Leonardo Di Venere / gammalib

Commit 44d3f497eb672a29adbe814218ff208dd1035bf0

Authored by Jürgen Knödlseder
1 parent f4cec345

Use CTA response cache in GCTAResponse and GCTAResponseIRF and remove former dif… 

…fuse model cache (#2715)
Inline Side-by-side
Showing 5 changed files with 297 additions and 325 deletions
... ... @@ -312,7 +312,12 @@ None
312 312 20. CTA interface
313 313 -----------------
314 314 The GCTAResponseCache class was added to implement a general response cache for
315   -CTA (#2715).
  315 +CTA. The response cache is now used in GCTAResponse::irf(GEvent&, GSource&,
  316 +GObservation&) and in GCTAResponseIrf::nroi(GModelSky&, GEnergy&, GTime&,
  317 +GEnergy&, GTime&, GObservation&) in case that the spatial model has no free
  318 +parameters. This means that now all other spatial models with frozen parameters
  319 +will be cached. The use of the cache can be toggled by internal data members
  320 +(#2715).
316 321  
317 322 GCTABackground3D::init_mc_cache now pre-computes the Monte Carlo cache at energy
318 323 nodes of the background cube, and in addition, inserts 8 logarithmically spaced
... ...
inst/cta/include/GCTAResponse.hpp
  View file @44d3f49
... ... @@ -30,6 +30,7 @@
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
32 32 #include "GResponse.hpp"
  33 +#include "GCTAResponseCache.hpp"
33 34  
34 35 /* __ Type definitions ___________________________________________________ */
35 36  
... ... @@ -113,6 +114,12 @@ protected:
113 114 virtual double irf_composite(const GEvent& event,
114 115 const GSource& source,
115 116 const GObservation& obs) const;
  117 +
  118 + // Protected members
  119 + bool m_use_irf_cache; //!< Control usage of IRF cache
  120 +
  121 + // Cache for irf(GEvent&, GSource&, GObservation&) computation
  122 + mutable GCTAResponseCache m_irf_cache;
116 123 };
117 124  
118 125 #endif /* GCTARESPONSE_HPP */
... ...
inst/cta/include/GCTAResponseIrf.hpp
  View file @44d3f49
... ... @@ -34,6 +34,7 @@
34 34 #include "GTime.hpp"
35 35 #include "GCaldb.hpp"
36 36 #include "GCTAResponse.hpp"
  37 +#include "GCTAResponseCache.hpp"
37 38 #include "GCTAAeff.hpp"
38 39 #include "GCTAPsf.hpp"
39 40 #include "GCTAEdisp.hpp"
... ... @@ -218,16 +219,15 @@ private:
218 219 mutable bool m_apply_edisp; //!< Apply energy dispersion
219 220 double m_lo_save_thres; //!< Save low energy threshold
220 221 double m_hi_save_thres; //!< Save high energy threshold
  222 + bool m_use_nroi_cache; //!< Control usage of nroi cache
221 223  
222 224 // XML response filename
223 225 std::string m_xml_caldb; //!< Calibration database string in XML file
224 226 std::string m_xml_rspname; //!< Response name in XML file
225 227  
226   - // Npred cache
227   - mutable std::vector<std::string> m_npred_names; //!< Model names
228   - mutable std::vector<GEnergy> m_npred_energies; //!< Model energy
229   - mutable std::vector<GTime> m_npred_times; //!< Model time
230   - mutable std::vector<double> m_npred_values; //!< Model values
  228 + // Cache for nroi(GModelSky&, GEnergy&, GTime&, GEnergy&, GTime&, GObservation&)
  229 + // computation
  230 + mutable GCTAResponseCache m_nroi_cache;
231 231 };
232 232  
233 233  
... ...
inst/cta/src/GCTAResponse.cpp
  View file @44d3f49
... ... @@ -34,6 +34,7 @@
34 34 #include "GModelSpatial.hpp"
35 35 #include "GModelSpatialComposite.hpp"
36 36 #include "GCTAResponse.hpp"
  37 +#include "GCTAResponse.hpp"
37 38  
38 39 /* __ Method name definitions ____________________________________________ */
39 40  
... ... @@ -166,25 +167,48 @@ double GCTAResponse::irf(const GEvent&amp; event,
166 167 // Initialise IRF value
167 168 double irf = 0.0;
168 169  
169   - // Select IRF depending on the spatial model type
170   - switch (source.model()->code()) {
171   - case GMODEL_SPATIAL_POINT_SOURCE:
172   - irf = irf_ptsrc(event, source, obs);
173   - break;
174   - case GMODEL_SPATIAL_RADIAL:
175   - irf = irf_radial(event, source, obs);
176   - break;
177   - case GMODEL_SPATIAL_ELLIPTICAL:
178   - irf = irf_elliptical(event, source, obs);
179   - break;
180   - case GMODEL_SPATIAL_DIFFUSE:
181   - irf = irf_diffuse(event, source, obs);
182   - break;
183   - case GMODEL_SPATIAL_COMPOSITE:
184   - irf = irf_composite(event, source, obs);
185   - break;
186   - default:
187   - break;
  170 + // Set IRF value attributes
  171 + std::string name = obs.id() + "::" + source.name();
  172 + const GEnergy& ereco = event.energy();
  173 + const GEnergy& etrue = source.energy();
  174 +
  175 + // Signal if spatial model has free parameters
  176 + bool has_free_pars = source.model()->has_free_pars();
  177 +
  178 + // If the spatial model component has free parameters, or the response
  179 + // cache should not be used, or the cache does not contain the requested
  180 + // IRF value then compute the IRF value for the spatial model.
  181 + if (has_free_pars ||
  182 + !m_use_irf_cache ||
  183 + !m_irf_cache.contains(name, ereco, etrue, &irf)) {
  184 +
  185 + // Compute IRF for spatial model
  186 + switch (source.model()->code()) {
  187 + case GMODEL_SPATIAL_POINT_SOURCE:
  188 + irf = irf_ptsrc(event, source, obs);
  189 + break;
  190 + case GMODEL_SPATIAL_RADIAL:
  191 + irf = irf_radial(event, source, obs);
  192 + break;
  193 + case GMODEL_SPATIAL_ELLIPTICAL:
  194 + irf = irf_elliptical(event, source, obs);
  195 + break;
  196 + case GMODEL_SPATIAL_DIFFUSE:
  197 + irf = irf_diffuse(event, source, obs);
  198 + break;
  199 + case GMODEL_SPATIAL_COMPOSITE:
  200 + irf = irf_composite(event, source, obs);
  201 + break;
  202 + default:
  203 + break;
  204 + }
  205 +
  206 + } // endif: computed spatial model
  207 +
  208 + // If the spatial model has no free parameters and the response cache
  209 + // should be used then put the IRF value in the response cache.
  210 + if (!has_free_pars && m_use_irf_cache) {
  211 + m_irf_cache.set(name, ereco, etrue, irf);
188 212 }
189 213  
190 214 // Return IRF value
... ... @@ -203,6 +227,10 @@ double GCTAResponse::irf(const GEvent&amp; event,
203 227 ***************************************************************************/
204 228 void GCTAResponse::init_members(void)
205 229 {
  230 + // Initialize members
  231 + m_use_irf_cache = true;
  232 + m_irf_cache.clear();
  233 +
206 234 // Return
207 235 return;
208 236 }
... ... @@ -215,6 +243,10 @@ void GCTAResponse::init_members(void)
215 243 ***************************************************************************/
216 244 void GCTAResponse::copy_members(const GCTAResponse& rsp)
217 245 {
  246 + // Copy members
  247 + m_use_irf_cache = rsp.m_use_irf_cache;
  248 + m_irf_cache = rsp.m_irf_cache;
  249 +
218 250 // Return
219 251 return;
220 252 }
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @44d3f49
... ... @@ -107,8 +107,6 @@
107 107 /* __ Macros _____________________________________________________________ */
108 108  
109 109 /* __ Coding definitions _________________________________________________ */
110   -#define G_USE_IRF_CACHE //!< Use IRF cache in irf_diffuse method
111   -#define G_USE_NPRED_CACHE //!< Use Npred cache in npred_diffuse method
112 110 //#define G_USE_PSF_SYSTEM //!< Do radial Irf integrations in Psf system
113 111  
114 112 /* __ Debug definitions __________________________________________________ */
... ... @@ -1604,8 +1602,8 @@ std::string GCTAResponseIrf::print(const GChatter&amp; chatter) const
1604 1602 }
1605 1603 }
1606 1604  
1607   - // Append save energy threshold information
1608   - result.append("\n"+gammalib::parformat("Save energy range"));
  1605 + // Append safe energy threshold information
  1606 + result.append("\n"+gammalib::parformat("Safe energy range"));
1609 1607 if (m_lo_save_thres > 0.0 && m_hi_save_thres) {
1610 1608 result.append(gammalib::str(m_lo_save_thres));
1611 1609 result.append(" - ");
... ... @@ -1626,9 +1624,11 @@ std::string GCTAResponseIrf::print(const GChatter&amp; chatter) const
1626 1624 result.append("undefined");
1627 1625 }
1628 1626  
1629   - // Append detailed information
  1627 + // Get reduced chatter level
1630 1628 GChatter reduced_chatter = gammalib::reduce(chatter);
1631   - if (reduced_chatter > SILENT) {
  1629 +
  1630 + // Append detailed information
  1631 + if (chatter >= NORMAL) {
1632 1632  
1633 1633 // Append calibration database
1634 1634 result.append("\n"+m_caldb.print(reduced_chatter));
... ... @@ -1655,19 +1655,11 @@ std::string GCTAResponseIrf::print(const GChatter&amp; chatter) const
1655 1655  
1656 1656 }
1657 1657  
1658   - // EXPLICIT: Append Npred cache information
  1658 + // Append cache information
1659 1659 if (chatter >= EXPLICIT) {
1660   - if (!m_npred_names.empty()) {
1661   - for (int i = 0; i < m_npred_names.size(); ++i) {
1662   - result.append("\n"+gammalib::parformat("Npred cache " +
1663   - gammalib::str(i)));
1664   - result.append(m_npred_names[i]+", ");
1665   - result.append(m_npred_energies[i].print()+", ");
1666   - result.append(m_npred_times[i].print()+" = ");
1667   - result.append(gammalib::str(m_npred_values[i]));
1668   - }
1669   - }
1670   - } // endif: chatter was explicit
  1660 + result.append("\n"+m_irf_cache.print(reduced_chatter));
  1661 + result.append("\n"+m_nroi_cache.print(reduced_chatter));
  1662 + }
1671 1663  
1672 1664 } // endif: chatter was not silent
1673 1665  
... ... @@ -1864,25 +1856,46 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
1864 1856 // Initialise response value
1865 1857 double nroi = 0.0;
1866 1858  
1867   - // Select method depending on the spatial model type
1868   - switch (model.spatial()->code()) {
1869   - case GMODEL_SPATIAL_POINT_SOURCE:
1870   - nroi = nroi_ptsrc(model, srcEng, srcTime, obsEng, obsTime, obs);
1871   - break;
1872   - case GMODEL_SPATIAL_RADIAL:
1873   - nroi = nroi_radial(model, srcEng, srcTime, obsEng, obsTime, obs);
1874   - break;
1875   - case GMODEL_SPATIAL_ELLIPTICAL:
1876   - nroi = nroi_elliptical(model, srcEng, srcTime, obsEng, obsTime, obs);
1877   - break;
1878   - case GMODEL_SPATIAL_DIFFUSE:
1879   - nroi = nroi_diffuse(model, srcEng, srcTime, obsEng, obsTime, obs);
1880   - break;
1881   - case GMODEL_SPATIAL_COMPOSITE:
1882   - nroi = nroi_composite(model, srcEng, srcTime, obsEng, obsTime, obs);
1883   - break;
1884   - default:
1885   - break;
  1859 + // Set response value attributes
  1860 + std::string name = obs.id() + "::" + model.name();
  1861 +
  1862 + // Signal if spatial model has free parameters
  1863 + bool has_free_pars = model.spatial()->has_free_pars();
  1864 +
  1865 + // If the spatial model component has free parameters, or the response
  1866 + // cache should not be used, or the cache does not contain the requested
  1867 + // IRF value then compute the IRF value for the spatial model.
  1868 + if (has_free_pars ||
  1869 + !m_use_nroi_cache ||
  1870 + !m_nroi_cache.contains(name, obsEng, srcEng, &nroi)) {
  1871 +
  1872 + // Select method depending on the spatial model type
  1873 + switch (model.spatial()->code()) {
  1874 + case GMODEL_SPATIAL_POINT_SOURCE:
  1875 + nroi = nroi_ptsrc(model, srcEng, srcTime, obsEng, obsTime, obs);
  1876 + break;
  1877 + case GMODEL_SPATIAL_RADIAL:
  1878 + nroi = nroi_radial(model, srcEng, srcTime, obsEng, obsTime, obs);
  1879 + break;
  1880 + case GMODEL_SPATIAL_ELLIPTICAL:
  1881 + nroi = nroi_elliptical(model, srcEng, srcTime, obsEng, obsTime, obs);
  1882 + break;
  1883 + case GMODEL_SPATIAL_DIFFUSE:
  1884 + nroi = nroi_diffuse(model, srcEng, srcTime, obsEng, obsTime, obs);
  1885 + break;
  1886 + case GMODEL_SPATIAL_COMPOSITE:
  1887 + nroi = nroi_composite(model, srcEng, srcTime, obsEng, obsTime, obs);
  1888 + break;
  1889 + default:
  1890 + break;
  1891 + }
  1892 +
  1893 + } // endif: computed spatial model
  1894 +
  1895 + // If the spatial model has no free parameters and the response cache
  1896 + // should be used then put the IRF value in the response cache.
  1897 + if (!has_free_pars && m_use_nroi_cache) {
  1898 + m_nroi_cache.set(name, obsEng, srcEng, nroi);
1886 1899 }
1887 1900  
1888 1901 // Return response value
... ... @@ -2105,23 +2118,21 @@ void GCTAResponseIrf::init_members(void)
2105 2118 // Initialise members
2106 2119 m_caldb.clear();
2107 2120 m_rspname.clear();
2108   - m_aeff = NULL;
2109   - m_psf = NULL;
2110   - m_edisp = NULL;
2111   - m_background = NULL;
2112   - m_apply_edisp = false; //!< Switched off by default
2113   - m_lo_save_thres = 0.0;
2114   - m_hi_save_thres = 0.0;
  2121 + m_aeff = NULL;
  2122 + m_psf = NULL;
  2123 + m_edisp = NULL;
  2124 + m_background = NULL;
  2125 + m_apply_edisp = false; //!< Switched off by default
  2126 + m_lo_save_thres = 0.0;
  2127 + m_hi_save_thres = 0.0;
  2128 + m_use_nroi_cache = true; //!< Switched on by default
2115 2129  
2116 2130 // XML response filenames
2117 2131 m_xml_caldb.clear();
2118 2132 m_xml_rspname.clear();
2119 2133  
2120   - // Initialise Npred cache
2121   - m_npred_names.clear();
2122   - m_npred_energies.clear();
2123   - m_npred_times.clear();
2124   - m_npred_values.clear();
  2134 + // Initialise nroi cache
  2135 + m_nroi_cache.clear();
2125 2136  
2126 2137 // Return
2127 2138 return;
... ... @@ -2136,21 +2147,19 @@ void GCTAResponseIrf::init_members(void)
2136 2147 void GCTAResponseIrf::copy_members(const GCTAResponseIrf& rsp)
2137 2148 {
2138 2149 // Copy members
2139   - m_caldb = rsp.m_caldb;
2140   - m_rspname = rsp.m_rspname;
2141   - m_apply_edisp = rsp.m_apply_edisp;
2142   - m_lo_save_thres = rsp.m_lo_save_thres;
2143   - m_hi_save_thres = rsp.m_hi_save_thres;
  2150 + m_caldb = rsp.m_caldb;
  2151 + m_rspname = rsp.m_rspname;
  2152 + m_apply_edisp = rsp.m_apply_edisp;
  2153 + m_lo_save_thres = rsp.m_lo_save_thres;
  2154 + m_hi_save_thres = rsp.m_hi_save_thres;
  2155 + m_use_nroi_cache = rsp.m_use_nroi_cache;
2144 2156  
2145 2157 // Copy response filenames
2146 2158 m_xml_caldb = rsp.m_xml_caldb;
2147 2159 m_xml_rspname = rsp.m_xml_rspname;
2148 2160  
2149   - // Copy cache
2150   - m_npred_names = rsp.m_npred_names;
2151   - m_npred_energies = rsp.m_npred_energies;
2152   - m_npred_times = rsp.m_npred_times;
2153   - m_npred_values = rsp.m_npred_values;
  2161 + // Copy nroi cache
  2162 + m_nroi_cache = rsp.m_nroi_cache;
2154 2163  
2155 2164 // Clone members
2156 2165 m_aeff = (rsp.m_aeff != NULL) ? rsp.m_aeff->clone() : NULL;
... ... @@ -2715,152 +2724,117 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp; event,
2715 2724 static const int max_iter_phi = 8;
2716 2725  
2717 2726 // Initialise IRF value
2718   - bool has_irf = false;
2719   - double irf = 0.0;
  2727 + double irf = 0.0;
2720 2728  
2721 2729 // Retrieve CTA pointing
2722 2730 const GCTAPointing& pnt = gammalib::cta_pnt(G_IRF_DIFFUSE, obs);
2723 2731  
2724   - // Try getting the IRF value from cache. Do not do this for energy
2725   - // dispersion since the the spectral model will change the IRF values
2726   - #if defined(G_USE_IRF_CACHE)
2727   - const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(obs.events());
2728   - const GCTAEventAtom* atom = dynamic_cast<const GCTAEventAtom*>(&event);
2729   - if (list != NULL && atom != NULL && !use_edisp()) {
2730   - irf = list->irf_cache(source.name(), atom->index());
2731   - if (irf >= 0.0) {
2732   - has_irf = true;
2733   - #if defined(G_DEBUG_IRF_DIFFUSE)
2734   - std::cout << "GCTAResponseIrf::irf_diffuse:";
2735   - std::cout << " cached irf=" << irf << std::endl;
2736   - #endif
2737   - }
2738   - else {
2739   - irf = 0.0;
2740   - }
2741   - }
2742   - #endif
2743   -
2744   - // Continue only if we have no IRF value
2745   - if (!has_irf) {
  2732 + // Get CTA instrument direction
  2733 + const GCTAInstDir& dir = gammalib::cta_dir(G_IRF_ELLIPTICAL, event);
2746 2734  
2747   - // Get CTA instrument direction
2748   - const GCTAInstDir& dir = gammalib::cta_dir(G_IRF_ELLIPTICAL, event);
2749   -
2750   - // Get pointer on spatial model
2751   - const GModelSpatial* model =
2752   - dynamic_cast<const GModelSpatial*>(source.model());
2753   - if (model == NULL) {
2754   - throw GCTAException::bad_model_type(G_IRF_DIFFUSE);
2755   - }
2756   -
2757   - // Get resolution of spatial model
2758   - double resolution = gammalib::resolution(model);
2759   -
2760   - // Get event attributes
2761   - const GEnergy& obsEng = event.energy();
  2735 + // Get pointer on spatial model
  2736 + const GModelSpatial* model =
  2737 + dynamic_cast<const GModelSpatial*>(source.model());
  2738 + if (model == NULL) {
  2739 + throw GCTAException::bad_model_type(G_IRF_DIFFUSE);
  2740 + }
2762 2741  
2763   - // Get source attributes
2764   - const GEnergy& srcEng = source.energy();
2765   - const GTime& srcTime = source.time();
  2742 + // Get resolution of spatial model
  2743 + double resolution = gammalib::resolution(model);
2766 2744  
2767   - // Get pointing direction zenith angle and azimuth [radians]
2768   - double zenith = pnt.zenith();
2769   - double azimuth = pnt.azimuth();
  2745 + // Get event attributes
  2746 + const GEnergy& obsEng = event.energy();
2770 2747  
2771   - // Determine angular distance between measured photon direction and
2772   - // pointing direction [radians]
2773   - double eta = pnt.dir().dist(dir.dir());
  2748 + // Get source attributes
  2749 + const GEnergy& srcEng = source.energy();
  2750 + const GTime& srcTime = source.time();
2774 2751  
2775   - // Get log10(E/TeV) of true photon energy
2776   - double srcLogEng = srcEng.log10TeV();
  2752 + // Get pointing direction zenith angle and azimuth [radians]
  2753 + double zenith = pnt.zenith();
  2754 + double azimuth = pnt.azimuth();
2777 2755  
2778   - // Assign the observed theta angle (eta) as the true theta angle
2779   - // between the source and the pointing directions. This is a (not
2780   - // too bad) approximation which helps to speed up computations.
2781   - // If we want to do this correctly, however, we would need to move
2782   - // the psf_dummy_sigma down to the integration kernel, and we would
2783   - // need to make sure that psf_delta_max really gives the absolute
2784   - // maximum (this is certainly less critical)
2785   - double theta = eta;
2786   - double phi = 0.0; //TODO: Implement Phi dependence
  2756 + // Determine angular distance between measured photon direction and
  2757 + // pointing direction [radians]
  2758 + double eta = pnt.dir().dist(dir.dir());
2787 2759  
2788   - // Get maximum PSF radius in radians
2789   - double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
  2760 + // Get log10(E/TeV) of true photon energy
  2761 + double srcLogEng = srcEng.log10TeV();
2790 2762  
2791   - // Perform zenith angle integration if interval is valid
2792   - if (delta_max > 0.0) {
  2763 + // Assign the observed theta angle (eta) as the true theta angle
  2764 + // between the source and the pointing directions. This is a (not
  2765 + // too bad) approximation which helps to speed up computations.
  2766 + // If we want to do this correctly, however, we would need to move
  2767 + // the psf_dummy_sigma down to the integration kernel, and we would
  2768 + // need to make sure that psf_delta_max really gives the absolute
  2769 + // maximum (this is certainly less critical)
  2770 + double theta = eta;
  2771 + double phi = 0.0; //TODO: Implement Phi dependence
2793 2772  
2794   - // Compute rotation matrix to convert from coordinates (theta,phi)
2795   - // in the reference frame of the observed arrival direction into
2796   - // celestial coordinates
2797   - GMatrix ry;
2798   - GMatrix rz;
2799   - ry.eulery(dir.dir().dec_deg() - 90.0);
2800   - rz.eulerz(-dir.dir().ra_deg());
2801   - GMatrix rot = (ry * rz).transpose();
  2773 + // Get maximum PSF radius in radians
  2774 + double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
2802 2775  
2803   - // Setup integration kernel
2804   - cta_irf_diffuse_kern_theta integrand(this,
2805   - model,
2806   - &rot,
2807   - theta,
2808   - phi,
2809   - zenith,
2810   - azimuth,
2811   - srcEng,
2812   - srcTime,
2813   - srcLogEng,
2814   - obsEng,
2815   - eta,
2816   - min_iter_phi,
2817   - max_iter_phi,
2818   - resolution);
2819   -
2820   - // Set number of radial integration iterations
2821   - int iter = gammalib::iter_rho(delta_max, resolution,
2822   - min_iter_rho, max_iter_rho);
2823   -
2824   - // Integrate over Psf delta angle
2825   - GIntegral integral(&integrand);
2826   - integral.fixed_iter(iter);
2827   - irf = integral.romberg(0.0, delta_max);
  2776 + // Perform zenith angle integration if interval is valid
  2777 + if (delta_max > 0.0) {
2828 2778  
2829   - // Compile option: Check for NaN/Inf
2830   - #if defined(G_NAN_CHECK)
2831   - if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
2832   - std::cout << "*** ERROR: GCTAResponseIrf::irf_diffuse:";
2833   - std::cout << " NaN/Inf encountered";
2834   - std::cout << " (irf=" << irf;
2835   - std::cout << ", delta_max=" << delta_max << ")";
2836   - std::cout << std::endl;
2837   - }
2838   - #endif
2839   - }
  2779 + // Compute rotation matrix to convert from coordinates (theta,phi)
  2780 + // in the reference frame of the observed arrival direction into
  2781 + // celestial coordinates
  2782 + GMatrix ry;
  2783 + GMatrix rz;
  2784 + ry.eulery(dir.dir().dec_deg() - 90.0);
  2785 + rz.eulerz(-dir.dir().ra_deg());
  2786 + GMatrix rot = (ry * rz).transpose();
2840 2787  
2841   - // Apply deadtime correction
2842   - irf *= obs.deadc(srcTime);
  2788 + // Setup integration kernel
  2789 + cta_irf_diffuse_kern_theta integrand(this,
  2790 + model,
  2791 + &rot,
  2792 + theta,
  2793 + phi,
  2794 + zenith,
  2795 + azimuth,
  2796 + srcEng,
  2797 + srcTime,
  2798 + srcLogEng,
  2799 + obsEng,
  2800 + eta,
  2801 + min_iter_phi,
  2802 + max_iter_phi,
  2803 + resolution);
  2804 +
  2805 + // Set number of radial integration iterations
  2806 + int iter = gammalib::iter_rho(delta_max, resolution,
  2807 + min_iter_rho, max_iter_rho);
  2808 +
  2809 + // Integrate over Psf delta angle
  2810 + GIntegral integral(&integrand);
  2811 + integral.fixed_iter(iter);
  2812 + irf = integral.romberg(0.0, delta_max);
2843 2813  
2844   - // Put IRF value in cache. Only do this in no energy dispersion is
2845   - // requested, since for energy dispersion the IRF values will
2846   - // depend on the spectrum
2847   - #if defined(G_USE_IRF_CACHE)
2848   - if (list != NULL && atom != NULL && !use_edisp()) {
2849   - list->irf_cache(source.name(), atom->index(), irf);
  2814 + // Compile option: Check for NaN/Inf
  2815 + #if defined(G_NAN_CHECK)
  2816 + if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
  2817 + std::cout << "*** ERROR: GCTAResponseIrf::irf_diffuse:";
  2818 + std::cout << " NaN/Inf encountered";
  2819 + std::cout << " (irf=" << irf;
  2820 + std::cout << ", delta_max=" << delta_max << ")";
  2821 + std::cout << std::endl;
2850 2822 }
2851 2823 #endif
  2824 + }
2852 2825  
2853   - // Compile option: Show integration results
2854   - #if defined(G_DEBUG_IRF_DIFFUSE)
2855   - std::cout << "GCTAResponseIrf::irf_diffuse:";
2856   - std::cout << " srcLogEng=" << srcLogEng;
2857   - std::cout << " obsLogEng=" << obsLogEng;
2858   - std::cout << " eta=" << eta;
2859   - std::cout << " delta_max=" << delta_max;
2860   - std::cout << " irf=" << irf << std::endl;
2861   - #endif
  2826 + // Apply deadtime correction
  2827 + irf *= obs.deadc(srcTime);
2862 2828  
2863   - } // endif: has no IRF
  2829 + // Compile option: Show integration results
  2830 + #if defined(G_DEBUG_IRF_DIFFUSE)
  2831 + std::cout << "GCTAResponseIrf::irf_diffuse:";
  2832 + std::cout << " srcLogEng=" << srcLogEng;
  2833 + std::cout << " obsLogEng=" << obsLogEng;
  2834 + std::cout << " eta=" << eta;
  2835 + std::cout << " delta_max=" << delta_max;
  2836 + std::cout << " irf=" << irf << std::endl;
  2837 + #endif
2864 2838  
2865 2839 // Return IRF value
2866 2840 return irf;
... ... @@ -3326,137 +3300,91 @@ double GCTAResponseIrf::nroi_diffuse(const GModelSky&amp; model,
3326 3300 static const int iter_phi = 9;
3327 3301  
3328 3302 // Initialise Nroi value
3329   - double nroi = 0.0;
3330   - bool has_nroi = false;
3331   -
3332   - // Build unique identifier
3333   - std::string id = model.name() + "::" + obs.id();
3334   -
3335   - // Check if Nroi value is already in cache. Do not use the Npred cache
3336   - // for energy dispersion since the actual value will depend on the
3337   - // spectral model
3338   - #if defined(G_USE_NPRED_CACHE)
3339   - if (!m_npred_names.empty() && !use_edisp()) {
3340   -
3341   - // Search for unique identifier, and if found, recover Npred value
3342   - // and break
3343   - for (int i = 0; i < m_npred_names.size(); ++i) {
3344   - if (m_npred_names[i] == id &&
3345   - m_npred_energies[i] == srcEng &&
3346   - m_npred_times[i] == srcTime) {
3347   - nroi = m_npred_values[i];
3348   - has_nroi = true;
3349   - #if defined(G_DEBUG_NROI_DIFFUSE)
3350   - std::cout << "GCTAResponseIrf::nroi_diffuse:";
3351   - std::cout << " cache=" << i;
3352   - std::cout << " nroi=" << nroi << std::endl;
3353   - #endif
3354   - break;
3355   - }
3356   - }
3357   -
3358   - } // endif: there were values in the Npred cache
3359   - #endif
3360   -
3361   - // Continue only if no Npred cache value was found
3362   - if (!has_nroi) {
  3303 + double nroi = 0.0;
3363 3304  
3364   - // Retrieve CTA observation, ROI and pointing
3365   - const GCTAObservation& cta = gammalib::cta_obs(G_NROI_DIFFUSE, obs);
3366   - const GCTARoi& roi = gammalib::cta_event_list(G_NROI_DIFFUSE, obs).roi();
3367   - const GCTAPointing& pnt = cta.pointing();
  3305 + // Retrieve CTA observation, ROI and pointing
  3306 + const GCTAObservation& cta = gammalib::cta_obs(G_NROI_DIFFUSE, obs);
  3307 + const GCTARoi& roi = gammalib::cta_event_list(G_NROI_DIFFUSE, obs).roi();
  3308 + const GCTAPointing& pnt = cta.pointing();
3368 3309  
3369   - // Get pointer on spatial model
3370   - const GModelSpatial* spatial =
3371   - dynamic_cast<const GModelSpatial*>(model.spatial());
3372   - if (spatial == NULL) {
3373   - throw GCTAException::bad_model_type(G_NROI_DIFFUSE);
3374   - }
  3310 + // Get pointer on spatial model
  3311 + const GModelSpatial* spatial =
  3312 + dynamic_cast<const GModelSpatial*>(model.spatial());
  3313 + if (spatial == NULL) {
  3314 + throw GCTAException::bad_model_type(G_NROI_DIFFUSE);
  3315 + }
3375 3316  
3376   - // Get pointing direction zenith angle and azimuth [radians]
3377   - double zenith = pnt.zenith();
3378   - double azimuth = pnt.azimuth();
  3317 + // Get pointing direction zenith angle and azimuth [radians]
  3318 + double zenith = pnt.zenith();
  3319 + double azimuth = pnt.azimuth();
3379 3320  
3380   - // Get log10(E/TeV) of true photon energy
3381   - double srcLogEng = srcEng.log10TeV();
  3321 + // Get log10(E/TeV) of true photon energy
  3322 + double srcLogEng = srcEng.log10TeV();
3382 3323  
3383   - // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
3384   - // as this allows us doing this computation in the outer loop. This
3385   - // should be sufficient here, unless the offaxis PSF becomes much worse
3386   - // than the onaxis PSF. In this case, we may add a safety factor here
3387   - // to make sure we encompass the entire PSF.
3388   - double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
  3324 + // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
  3325 + // as this allows us doing this computation in the outer loop. This
  3326 + // should be sufficient here, unless the offaxis PSF becomes much worse
  3327 + // than the onaxis PSF. In this case, we may add a safety factor here
  3328 + // to make sure we encompass the entire PSF.
  3329 + double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
3389 3330  
3390   - // Extract ROI radius (radians)
3391   - double roi_radius = roi.radius() * gammalib::deg2rad;
  3331 + // Extract ROI radius (radians)
  3332 + double roi_radius = roi.radius() * gammalib::deg2rad;
3392 3333  
3393   - // Compute the ROI radius plus maximum PSF radius (radians). Any photon
3394   - // coming from beyond this radius will not make it in the dataspace and
3395   - // thus can be neglected.
3396   - double roi_psf_radius = roi_radius + psf_max_radius;
  3334 + // Compute the ROI radius plus maximum PSF radius (radians). Any photon
  3335 + // coming from beyond this radius will not make it in the dataspace and
  3336 + // thus can be neglected.
  3337 + double roi_psf_radius = roi_radius + psf_max_radius;
3397 3338  
3398   - // Perform offset angle integration only if interval is valid
3399   - if (roi_psf_radius > 0.0) {
  3339 + // Perform offset angle integration only if interval is valid
  3340 + if (roi_psf_radius > 0.0) {
3400 3341  
3401   - // Compute rotation matrix to convert from native ROI coordinates,
3402   - // given by (theta,phi), into celestial coordinates.
3403   - GMatrix ry;
3404   - GMatrix rz;
3405   - ry.eulery(roi.centre().dir().dec_deg() - 90.0);
3406   - rz.eulerz(-roi.centre().dir().ra_deg());
3407   - GMatrix rot = (ry * rz).transpose();
  3342 + // Compute rotation matrix to convert from native ROI coordinates,
  3343 + // given by (theta,phi), into celestial coordinates.
  3344 + GMatrix ry;
  3345 + GMatrix rz;
  3346 + ry.eulery(roi.centre().dir().dec_deg() - 90.0);
  3347 + rz.eulerz(-roi.centre().dir().ra_deg());
  3348 + GMatrix rot = (ry * rz).transpose();
3408 3349  
3409   - // Setup integration kernel
3410   - cta_nroi_diffuse_kern_theta integrand(this,
3411   - &cta,
3412   - spatial,
3413   - &rot,
3414   - srcEng,
3415   - srcTime,
3416   - obsEng,
3417   - obsTime,
3418   - iter_phi);
3419   -
3420   - // Integrate over model's zenith angle
3421   - GIntegral integral(&integrand);
3422   - integral.fixed_iter(iter_rho);
3423   - nroi = integral.romberg(0.0, roi_psf_radius);
3424   -
3425   - // Compile option: Show integration results
3426   - #if defined(G_DEBUG_NROI_DIFFUSE)
3427   - std::cout << "GCTAResponseIrf::nroi_diffuse:";
3428   - std::cout << " roi_psf_radius=" << roi_psf_radius;
3429   - std::cout << " nroi=" << nroi;
3430   - std::cout << " Etrue=" << srcEng;
3431   - std::cout << " Ereco=" << obsEng;
3432   - std::cout << " id=" << id << std::endl;
3433   - #endif
  3350 + // Setup integration kernel
  3351 + cta_nroi_diffuse_kern_theta integrand(this,
  3352 + &cta,
  3353 + spatial,
  3354 + &rot,
  3355 + srcEng,
  3356 + srcTime,
  3357 + obsEng,
  3358 + obsTime,
  3359 + iter_phi);
3434 3360  
3435   - } // endif: offset angle range was valid
  3361 + // Integrate over model's zenith angle
  3362 + GIntegral integral(&integrand);
  3363 + integral.fixed_iter(iter_rho);
  3364 + nroi = integral.romberg(0.0, roi_psf_radius);
3436 3365  
3437   - // Store result in Npred cache in case that no energy dispersion
3438   - // is requested
3439   - #if defined(G_USE_NPRED_CACHE)
3440   - if (!use_edisp()) {
3441   - m_npred_names.push_back(id);
3442   - m_npred_energies.push_back(srcEng);
3443   - m_npred_times.push_back(srcTime);
3444   - m_npred_values.push_back(nroi);
3445   - }
  3366 + // Compile option: Show integration results
  3367 + #if defined(G_DEBUG_NROI_DIFFUSE)
  3368 + std::cout << "GCTAResponseIrf::nroi_diffuse:";
  3369 + std::cout << " roi_psf_radius=" << roi_psf_radius;
  3370 + std::cout << " nroi=" << nroi;
  3371 + std::cout << " Etrue=" << srcEng;
  3372 + std::cout << " Ereco=" << obsEng;
  3373 + std::cout << " id=" << id << std::endl;
3446 3374 #endif
3447 3375  
3448   - // Debug: Check for NaN
3449   - #if defined(G_NAN_CHECK)
3450   - if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
3451   - std::cout << "*** ERROR: GCTAResponseIrf::nroi_diffuse:";
3452   - std::cout << " NaN/Inf encountered";
3453   - std::cout << " (nroi=" << nroi;
3454   - std::cout << ", roi_psf_radius=" << roi_psf_radius;
3455   - std::cout << ")" << std::endl;
3456   - }
3457   - #endif
  3376 + } // endif: offset angle range was valid
3458 3377  
3459   - } // endif: Nroi computation required
  3378 + // Debug: Check for NaN
  3379 + #if defined(G_NAN_CHECK)
  3380 + if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
  3381 + std::cout << "*** ERROR: GCTAResponseIrf::nroi_diffuse:";
  3382 + std::cout << " NaN/Inf encountered";
  3383 + std::cout << " (nroi=" << nroi;
  3384 + std::cout << ", roi_psf_radius=" << roi_psf_radius;
  3385 + std::cout << ")" << std::endl;
  3386 + }
  3387 + #endif
3460 3388  
3461 3389 // Return Nroi
3462 3390 return nroi;
... ...