diff --git a/ChangeLog b/ChangeLog index 6491204..d9a724e 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,9 +1,9 @@ -2017-03-18 +2017-03-21 * Version 1.3.0 released ======================== - Take into account RoI in ctskymap (#1975) + Take into account RoI in ctskymap and use numerical integration for IRF (#1975) Set local language to english (#1370) Enhance logfile output of csroot2caldb script Set default deadc value to 0.98 (CTA has <2% deadtime) diff --git a/NEWS b/NEWS index 2ffb2e8..2148b77 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.3.0.dev1 -18 March 2017 +21 March 2017 Introduction @@ -53,7 +53,10 @@ None ctskymap - CTA sky mapping -------------------------- Any RoI is now taken into account in ctskymap so that if IRF background -subtraction is requested the counts and background maps are coherent +subtraction is requested the counts and background maps are coherent. In +addition, the energy integration over the background template is now done +numerically so that strong variations in the background rate, such as +those occuring close to the threshold, are correctly taken into account (#1975). diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index 072bbac..ff5df67 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -30,7 +30,6 @@ #endif #include #include "ctskymap.hpp" -#include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ #define G_INIT_MAP "ctskymap::init_map(GCTAObservation* obs)" @@ -584,7 +583,8 @@ void ctskymap::map_background(GCTAObservation* obs) * No background template available in instrument response function. * * Estimates the background in the sky map using the IRF template for a given - * observation and adds this estimate to the background sky map. + * observation and adds this estimate to the background sky map. The IRF + * template is integrated numerically in energy. ***************************************************************************/ void ctskymap::map_background_irf(GCTAObservation* obs) { @@ -612,24 +612,16 @@ void ctskymap::map_background_irf(GCTAObservation* obs) throw GException::invalid_value(G_BKG_SUBTRACT_IRF, msg); } - // Setup energy range covered by data - GEnergy emin; - GEnergy emax; - GEbounds ebds; - emin.TeV(m_emin); - emax.TeV(m_emax); - - // Determine number of energy nodes, 20 per decade, at least 10 bins - int nebin = int((emax.log10TeV() - emin.log10TeV()) * 20.0); - if (nebin < 10) { - nebin = 10; - } + // Compute natural logarithm of energy range in MeV + double lnEmin = std::log(m_emin * 1.0e6); + double lnEmax = std::log(m_emax * 1.0e6); // Extract region of interest from observation GCTARoi roi = obs->roi(); - // Setup energy nodes for computation - GEnergies energies(nebin, emin, emax); // More complicated nbins later + // Initialise statistics + int calls = 0; + double total = 0.0; // Loop over all background map pixels for (int i = 0; i < m_bkgmap.npix(); ++i) { @@ -646,28 +638,57 @@ void ctskymap::map_background_irf(GCTAObservation* obs) continue; } - // Compute background value for each energy node - std::vector values(energies.size(), 0.0); - for (int k = 0; k < energies.size(); ++k) { - values[k] = (*bkg)(energies[k].log10TeV(), instdir.detx(), - instdir.dety()); - } + // Setup integration function + ctskymap::irf_kern integrand(bkg, &instdir); + GIntegral integral(&integrand); - // Integrate background values according to Trapezoid rule - double value = 0; - for (int k = 1; k < energies.size(); ++k) { - value += 0.5 * (values[k] + values[k-1]) * - (energies[k].MeV() - energies[k-1].MeV()); - } + // Set precision (has been carefully adjusted using a test simulation + // over the energy range 20 GeV - 120 TeV) + integral.eps(1.0e-7); + // Do Romberg integration + double value = integral.romberg(lnEmin, lnEmax); + + // Update number of background function calls + calls += integral.calls(); + // Multiply background rate with livetime and solid angle value *= obs->livetime() * m_bkgmap.solidangle(i); // Add number of background events to background map m_bkgmap(i,0) += value; + // Update total number of background events + total += value; + } // endfor: looped over background map pixels + + // Log background subtraction results + log_value(NORMAL, "Events in background", int(total+0.5)); + log_value(NORMAL, "Background evaluations", calls); // Return return; } + + +/***********************************************************************//** + * @brief Estimates the background in sky map based on IRF template + * + * @param[in] lnE Natural logarithm of energy in MeV. + ***************************************************************************/ +double ctskymap::irf_kern::eval(const double& lnE) +{ + // Get log10 of energy in TeV + double logE = lnE * gammalib::inv_ln10 - 6; + + // Get function value + double value = (*m_bgd)(logE, m_dir->detx(), m_dir->dety()); + + // Correct for variable substitution + value *= std::exp(lnE); + + // Return value + return value; +} + diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index c4efd25..67f0985 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -73,6 +73,19 @@ protected: void map_background(GCTAObservation* obs); void map_background_irf(GCTAObservation* obs); + // Background integration kernel + class irf_kern : public GFunction { + public: + irf_kern(const GCTABackground* bgd, + const GCTAInstDir* dir) : + m_bgd(bgd), + m_dir(dir) { } + double eval(const double& lnE); + protected: + const GCTABackground* m_bgd; //!< Pointer to background + const GCTAInstDir* m_dir; //!< Pointer to instrument direction + }; + // User parameters GFilename m_outmap; //!< Output file name double m_emin; //!< Minimum energy (TeV)