Eirik Hatlen / gammalib

Commit 2110ec8a8fe69e078469f1a7e5f389e756d0452d

Authored by Jürgen Knödlseder
2 parents 29ed2bd1 c4911769

Merge branch 'bugfix-1.7.3' into integration 

# Conflicts:
#	ChangeLog
#	NEWS
#	README.md
#	configure.ac
#	doc/Doxyfile
#	doc/source/admin/release_history/index.rst
#	doc/source/conf.py
#	sonar-project.properties
Inline Side-by-side
Showing 12 changed files with 133 additions and 44 deletions
ChangeLog
  View file @2110ec8
... ... @@ -26,6 +26,14 @@
26 26 Implement pre-computation cache in GModelSpatialGauss (#1300)
27 27  
28 28  
  29 +2020-12-02
  30 +
  31 + * Bug fix version 1.7.3 released
  32 + ================================
  33 +
  34 + Fix inaccuracy in CTA background model evaluation (#3472)
  35 +
  36 +
29 37 2020-10-15
30 38  
31 39 * Bug fix version 1.7.2 released
... ...
dev/conda.recipe/conda_build_config.yaml
  View file @2110ec8
... ... @@ -3,3 +3,4 @@ python:
3 3 - 3.5
4 4 - 3.6
5 5 - 3.7
  6 + - 3.8
... ...
doc/source/admin/download.rst
  View file @2110ec8
... ... @@ -4,13 +4,14 @@ Download
4 4 ========
5 5  
6 6 The latest GammaLib release is
7   -`gammalib-1.7.1 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.1.tar.gz>`_
8   -(7 September 2020).
  7 +`gammalib-1.7.2 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.2.tar.gz>`_
  8 +(15 October 2020).
9 9  
10 10 Below a list of older GammaLib releases. Please read the :ref:`release` to
11 11 learn more about new features and corrected bugs in a given release.
12 12  
13 13 To download an older release, click on the corresponding release number:
  14 +`1.7.1 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.1.tar.gz>`_
14 15 `1.7.0 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.7.0.tar.gz>`_
15 16 `1.6.3 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.6.3.tar.gz>`_
16 17 `1.6.2 <http://cta.irap.omp.eu/ctools/releases/gammalib/gammalib-1.6.2.tar.gz>`_
... ...
doc/source/admin/release_history/1.7.2.rst 0 → 100644
  View file @2110ec8
  1 +.. _1.7.2:
  2 +
  3 +GammaLib 1.7.2 release (15 October 2020)
  4 +========================================
  5 +
  6 +Overview
  7 +--------
  8 +
  9 +GammaLib 1.7.2 is a bug fix release for GammaLib release 1.7.
  10 +
  11 +The following bugs were fixed:
  12 +
  13 +* [`3402 <https://cta-redmine.irap.omp.eu/issues/3402>`_] -
  14 + Fix copy of base class members in copy constructor of some spatial models.
... ...
doc/source/admin/release_history/1.7.3.rst 0 → 100644
  View file @2110ec8
  1 +.. _1.7.3:
  2 +
  3 +GammaLib 1.7.3 release (2 December 2020)
  4 +========================================
  5 +
  6 +Overview
  7 +--------
  8 +
  9 +GammaLib 1.7.3 is a bug fix release for GammaLib release 1.7.
  10 +
  11 +The following bugs were fixed:
  12 +
  13 +* [`3472 <https://cta-redmine.irap.omp.eu/issues/3472>`_] -
  14 + Fix inaccuracy in CTA background model evaluation
... ...
doc/source/admin/release_history/index.rst
  View file @2110ec8
... ... @@ -7,6 +7,8 @@ Release History
7 7 :maxdepth: 1
8 8  
9 9 GammaLib 2.0 <2.0>
  10 + GammaLib 1.7.3 <1.7.3>
  11 + GammaLib 1.7.2 <1.7.2>
10 12 GammaLib 1.7.1 <1.7.1>
11 13 GammaLib 1.7 <1.7>
12 14 GammaLib 1.6.3 <1.6.3>
... ...
inst/cta/include/GCTAEventBin.hpp
  View file @2110ec8
1 1 /***************************************************************************
2 2 * GCTAEventBin.hpp - CTA event bin class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -80,6 +80,8 @@ public:
80 80 const int& ipix(void) const;
81 81 const int& ieng(void) const;
82 82 const double& solidangle(void) const;
  83 + GEnergy emin(void) const;
  84 + GEnergy emax(void) const;
83 85 const GEnergy& ewidth(void) const;
84 86 const double& ontime(void) const;
85 87 const double& weight(void) const;
... ... @@ -185,4 +187,18 @@ void GCTAEventBin::ieng(const int&amp; ieng)
185 187 return;
186 188 }
187 189  
  190 +
  191 +/***********************************************************************//**
  192 + * @brief Return maximum energy of event bin
  193 + *
  194 + * @return Maximum energy of event bin.
  195 + *
  196 + * Returns the maximum energy of event bin.
  197 + ***************************************************************************/
  198 +inline
  199 +GEnergy GCTAEventBin::emax(void) const
  200 +{
  201 + return (emin() + ewidth());
  202 +}
  203 +
188 204 #endif /* GCTAEVENTBIN_HPP */
... ...
inst/cta/pyext/GCTAEventBin.i
  View file @2110ec8
1 1 /***************************************************************************
2 2 * GCTAEventBin.i - CTA event bin class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -61,6 +61,8 @@ public:
61 61 const int& ipix(void) const;
62 62 const int& ieng(void) const;
63 63 const double& solidangle(void) const;
  64 + GEnergy emin(void) const;
  65 + GEnergy emax(void) const;
64 66 const GEnergy& ewidth(void) const;
65 67 const double& ontime(void) const;
66 68 const double& weight(void) const;
... ...
inst/cta/src/GCTACubeBackground.cpp
  View file @2110ec8
... ... @@ -466,44 +466,23 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
466 466 // DETX/DETY coordinates
467 467 GCTAInstDir dir = cta->pointing().instdir(bin->dir().dir());
468 468  
469   - // Initialise values for integration over energy bin
470   - double f1 = 0.0;
471   -
472 469 // Loop over all energy bins of the cube
473 470 for (int iebin = 0, ibin = ipix; iebin < nebins; ++iebin, ibin += npix) {
474 471  
475   - // If this is the first bin then compute f1 at the lower bin edge
476   - if (iebin == 0) {
477   - GCTAEventBin* bin = eventcube[ibin];
478   - bin->dir(dir); // Set instrument direction
479   - bin->energy(m_ebounds.emin(iebin));
480   - f1 = obs.models().eval(*bin, *cta);
481   - }
482   -
483   - // Compute f2 at the upper bin edge
484   - GCTAEventBin* bin = eventcube[ibin];
485   - bin->dir(dir); // Set instrument direction
486   - bin->energy(m_ebounds.emax(iebin));
487   - double f2 = obs.models().eval(*bin, *cta);
488   -
489   - // Add background model value to cube if the energy bin is
490   - // fully container in the observation interval and if both f1
491   - // and f2 are positive. Add a small margin here to the energy
492   - // bin boundaries to take provision for rounding errors.
  472 + // Compute background model value if the energy bin is fully
  473 + // contained in the observation interval. Add a small margin
  474 + // here to the energy bin boundaries to take provision for
  475 + // rounding errors.
493 476 if (obs_ebounds.contains(m_ebounds.emin(iebin) + energy_margin,
494   - m_ebounds.emax(iebin) - energy_margin) &&
495   - (f1 > 0.0 && f2 > 0.0)) {
  477 + m_ebounds.emax(iebin) - energy_margin)) {
496 478  
497   - // Compute background model value by integrating over the
498   - // bin and deviding the result by the bin width
499   - double x1 = m_ebounds.emin(iebin).MeV();
500   - double x2 = m_ebounds.emax(iebin).MeV();
501   - double model = gammalib::plaw_integral(x1, f1, x2, f2) / (x2 - x1);
  479 + // Set energy bin
  480 + GCTAEventBin* bin = eventcube[ibin];
  481 + bin->dir(dir); // Set instrument direction
  482 + bin->energy(m_ebounds.elogmean(iebin));
502 483  
503   - // Multiply by ontime to get the correct weighting for each
504   - // observation. We divide by the total ontime later to get
505   - // the background model in units of counts/(MeV s sr).
506   - model *= ontime;
  484 + // Compute model
  485 + double model = obs.models().eval(*bin, *cta) * ontime;
507 486  
508 487 // Add existing number of counts
509 488 model += bin->counts();
... ... @@ -513,10 +492,7 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
513 492  
514 493 } // endif: added background model value to cube
515 494  
516   - // Set f1 to f2 for next energy layer
517   - f1 = f2;
518   -
519   - } // endfor: looped over energy layers
  495 + } // endfor: looped over energy bins
520 496  
521 497 } // endfor: looped over spatial pixels
522 498  
... ...
inst/cta/src/GCTAEventBin.cpp
  View file @2110ec8
1 1 /***************************************************************************
2 2 * GCTAEventBin.cpp - CTA event bin class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -342,6 +342,40 @@ const double&amp; GCTAEventBin::solidangle(void) const
342 342  
343 343  
344 344 /***********************************************************************//**
  345 + * @brief Return minimum energy of event bin
  346 + *
  347 + * @return Minimum energy of event bin
  348 + *
  349 + * Returns minimum energy of event bin, computed using
  350 + *
  351 + * \f[
  352 + * E_{\rm min} = \frac{-\Delta E + sqrt{\Delta E^2 + 4 E^2}}{2}
  353 + * \f]
  354 + *
  355 + * where
  356 + * \f$\Delta E\f$ is the energy bin width, returned by ewidth() and
  357 + * \f$E\f$ is the energy, returned by energy().
  358 + ***************************************************************************/
  359 +GEnergy GCTAEventBin::emin(void) const
  360 +{
  361 + // Get energy and energy width in MeV
  362 + double energy = this->energy().MeV();
  363 + double ewidth = this->ewidth().MeV();
  364 +
  365 + // Compute minimum energy in MeV
  366 + double arg = ewidth * ewidth + 4.0 * energy * energy;
  367 + double e_min = 0.5 * (-ewidth + std::sqrt(arg));
  368 +
  369 + // Set minimum energy
  370 + GEnergy emin;
  371 + emin.MeV(e_min);
  372 +
  373 + // Return minimum energy
  374 + return emin;
  375 +}
  376 +
  377 +
  378 +/***********************************************************************//**
345 379 * @brief Return energy width of event bin
346 380 *
347 381 * @return Energy width of event bin
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @2110ec8
... ... @@ -40,6 +40,7 @@
40 40 #include "GCTAObservation.hpp"
41 41 #include "GCTAResponseIrf.hpp"
42 42 #include "GCTABackground.hpp"
  43 +#include "GCTAEventBin.hpp"
43 44 #include "GCTASupport.hpp"
44 45  
45 46 /* __ Globals ____________________________________________________________ */
... ... @@ -58,9 +59,10 @@ const GModelRegistry g_cta_inst_background_registry(&amp;g_cta_inst_backgrou
58 59  
59 60 /* __ Coding definitions _________________________________________________ */
60 61 #define G_USE_NPRED_CACHE
  62 +#define G_USE_RATE_EBIN
61 63  
62 64 /* __ Debug definitions __________________________________________________ */
63   -//#define G_DUMP_MC
  65 +#define G_DUMP_MC
64 66 //#define G_DEBUG_NPRED
65 67  
66 68 /* __ Constants __________________________________________________________ */
... ... @@ -349,10 +351,27 @@ double GCTAModelIrfBackground::eval(const GEvent&amp; event,
349 351  
350 352 // Set DETX and DETY in instrument direction
351 353 GCTAInstDir inst_dir = pnt.instdir(dir.dir());
352   -
353   - // Evaluate function
  354 +
  355 + // Evaluate spatial component
  356 + #if defined(G_USE_RATE_EBIN)
  357 + double spat = 0.0;
  358 + if (event.is_atom()) {
  359 + double logE = event.energy().log10TeV();
  360 + spat = bgd(logE, inst_dir.detx(), inst_dir.dety());
  361 + }
  362 + else {
  363 + const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
  364 + GEnergy emin = bin->emin();
  365 + GEnergy emax = emin + bin->ewidth();
  366 + double norm = 1.0 / bin->ewidth().MeV();
  367 + spat = bgd.rate_ebin(inst_dir, emin, emax) * norm;
  368 + }
  369 + #else
354 370 double logE = event.energy().log10TeV();
355 371 double spat = bgd(logE, inst_dir.detx(), inst_dir.dety());
  372 + #endif
  373 +
  374 + // Evaluate function
356 375 double spec = (spectral() != NULL)
357 376 ? spectral()->eval(event.energy(), event.time(), gradients)
358 377 : 1.0;
... ...
inst/cta/test/test_CTA.cpp
  View file @2110ec8
... ... @@ -2294,6 +2294,8 @@ void TestGCTAObservation::test_event_bin(void)
2294 2294 test_value(bin.ieng(), 7, "Check ieng() for filled bin");
2295 2295 test_value(bin.solidangle(), 3.14, 1.0e-10,
2296 2296 "Check solidangle() for filled bin");
  2297 + test_value(bin.emin().TeV(), 1.996503, "Check emin() for filled bin");
  2298 + test_value(bin.emax().TeV(), 2.003503, "Check emax() for filled bin");
2297 2299 test_value(bin.ewidth().GeV(), 7.0, 1.0e-10,
2298 2300 "Check ewidth() for filled bin");
2299 2301 test_value(bin.ontime(), 101.0, 1.0e-10, "Check ontime() for filled bin");
... ...