Eirik Hatlen / gammalib

Commit c4911769f899d44f0f2c60b68d67ffea2bbf06c7

Authored by Jürgen Knödlseder
1 parent cd207aa3

Use rate_ebin() method for CTA IRF background model evaluation (#3472) 

The GCTAModelIrfBackground::eval() method uses now the GCTABackground3D::rate_ebin() integrator method instead of the GCTABackground3D::operator() operator so that the background model evaluation is no longer sensitive to the number of energy bins.

The GCTACubeBackground::fill() method was simplified, removing the power law integration which was not accurate enough, and using directly the GCTAModelIrfBackground::eval() method which is now exact.
Inline Side-by-side
Showing 11 changed files with 124 additions and 51 deletions
ChangeLog
  View file @c491176
  1 +2020-12-02
  2 +
  3 + * Bug fix version 1.7.3 released
  4 + ================================
  5 +
  6 + Fix inaccuracy in CTA background model evaluation (#3472)
  7 +
  8 +
1 9 2020-10-15
2 10  
3 11 * Bug fix version 1.7.2 released
... ...
1   -New Features and Important Changes in GammaLib 1.7.2
  1 +New Features and Important Changes in GammaLib 1.7.3
2 2  
3   -15 October 2020
  3 +2 December 2020
4 4  
5 5  
6 6 1. Introduction
... ... @@ -63,11 +63,7 @@ None
63 63  
64 64 8. Model module
65 65 ---------------
66   -The copy constructors of the GModelSpatialDiffuse, GModelSpatialElliptical
67   -and GModelSpatialRadial classes did not properly copy the base class members
68   -in their copy constructors. This prevented specifically the fitting of
69   -diffuse maps in the stacked analysis. The code was modified to fix the
70   -issue.
  66 +None
71 67  
72 68  
73 69 9. Numerics module
... ... @@ -127,7 +123,12 @@ None
127 123  
128 124 20. CTA interface
129 125 -----------------
130   -None
  126 +The GCTAModelIrfBackground::eval() method uses now the GCTABackground3D::rate_ebin()
  127 +integrator method instead of the GCTABackground3D::operator() operator so that the
  128 +background model evaluation is no longer sensitive to the number of energy bins.
  129 +The GCTACubeBackground::fill() method was simplified, removing the power law
  130 +integration which was not accurate enough, and using directly the
  131 +GCTAModelIrfBackground::eval() method which is now exact (#3472).
131 132  
132 133  
133 134 21. COMPTEL interface
... ...
README.md
  View file @c491176
1 1 GammaLib information
2 2 ====================
3   -* Version: 1.7.2 (15 October 2020)
  3 +* Version: 1.7.3 (2 December 2020)
4 4  
5 5 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
6 6  
... ...
doc/source/admin/release_history/1.7.3.rst 0 → 100644
  View file @c491176
  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 @c491176
... ... @@ -6,6 +6,7 @@ Release History
6 6 .. toctree::
7 7 :maxdepth: 1
8 8  
  9 + GammaLib 1.7.3 <1.7.3>
9 10 GammaLib 1.7.2 <1.7.2>
10 11 GammaLib 1.7.1 <1.7.1>
11 12 GammaLib 1.7 <1.7>
... ...
inst/cta/include/GCTAEventBin.hpp
  View file @c491176
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 @c491176
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 @c491176
... ... @@ -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 @c491176
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 @c491176
... ... @@ -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 @c491176
... ... @@ -2229,6 +2229,8 @@ void TestGCTAObservation::test_event_bin(void)
2229 2229 test_value(bin.ieng(), 7, "Check ieng() for filled bin");
2230 2230 test_value(bin.solidangle(), 3.14, 1.0e-10,
2231 2231 "Check solidangle() for filled bin");
  2232 + test_value(bin.emin().TeV(), 1.996503, "Check emin() for filled bin");
  2233 + test_value(bin.emax().TeV(), 2.003503, "Check emax() for filled bin");
2232 2234 test_value(bin.ewidth().GeV(), 7.0, 1.0e-10,
2233 2235 "Check ewidth() for filled bin");
2234 2236 test_value(bin.ontime(), 101.0, 1.0e-10, "Check ontime() for filled bin");
... ...