Justine Devin / gammalib

Commit ee42aa1f2fb6f284cedd2c2f8d8eaa10a9bd58b9

Authored by Jürgen Knödlseder
1 parent ea998dc4

Started to re-organise response interface. Code still needs to be cleaned up.

Inline Side-by-side
Showing 23 changed files with 1241 additions and 461 deletions
include/GModelSky.hpp
  View file @ee42aa1
1 1 /***************************************************************************
2 2 * GModelSky.hpp - Sky model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
include/GResponse.hpp
  View file @ee42aa1
... ... @@ -43,6 +43,7 @@
43 43  
44 44 /* __ Forward declarations _______________________________________________ */
45 45 class GObservation;
  46 +class GModelSky;
46 47  
47 48  
48 49 /***********************************************************************//**
... ... @@ -81,7 +82,8 @@ public:
81 82 const GObservation& obs) const = 0;
82 83 virtual double convolve(const GModelSky& model,
83 84 const GEvent& event,
84   - const GObservation& obs) const;
  85 + const GObservation& obs,
  86 + const bool& grad = true) const;
85 87 virtual double nroi(const GModelSky& model,
86 88 const GEnergy& obsEng,
87 89 const GTime& obsTime,
... ... @@ -89,8 +91,6 @@ public:
89 91 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
90 92  
91 93 // Old methods that will become obsolete
92   - virtual double npred(const GPhoton& photon,
93   - const GObservation& obs) const = 0;
94 94 virtual double irf(const GEvent& event,
95 95 const GSource& source,
96 96 const GObservation& obs) const;
... ... @@ -106,6 +106,7 @@ public:
106 106 virtual double irf_diffuse(const GEvent& event,
107 107 const GSource& source,
108 108 const GObservation& obs) const;
  109 + /*
109 110 virtual double npred(const GSource& source,
110 111 const GObservation& obs) const;
111 112 virtual double npred_ptsrc(const GSource& source,
... ... @@ -116,6 +117,7 @@ public:
116 117 const GObservation& obs) const;
117 118 virtual double npred_diffuse(const GSource& source,
118 119 const GObservation& obs) const;
  120 + */
119 121 virtual GEbounds ebounds_src(const GEnergy& obsEng) const;
120 122  
121 123 protected:
... ... @@ -125,6 +127,7 @@ protected:
125 127 void free_members(void);
126 128  
127 129 // Npred theta integration kernel for radial model
  130 + /*
128 131 class npred_radial_kern_theta : public GFunction {
129 132 public:
130 133 npred_radial_kern_theta(const GResponse& rsp,
... ... @@ -148,8 +151,10 @@ protected:
148 151 const GObservation& m_obs; //!< Observation
149 152 const GMatrix& m_rot; //!< Rotation matrix
150 153 };
  154 + */
151 155  
152 156 // Npred phi integration kernel for radial model
  157 + /*
153 158 class npred_radial_kern_phi : public GFunction {
154 159 public:
155 160 npred_radial_kern_phi(const GResponse& rsp,
... ... @@ -178,8 +183,10 @@ protected:
178 183 double m_cos_theta; //!< cosine of offset angle
179 184 const double& m_sin_theta; //!< Sine of offset angle
180 185 };
  186 + */
181 187  
182 188 // Npred theta integration kernel for elliptical model
  189 + /*
183 190 class npred_elliptical_kern_theta : public GFunction {
184 191 public:
185 192 npred_elliptical_kern_theta(const GResponse& rsp,
... ... @@ -203,8 +210,10 @@ protected:
203 210 const GObservation& m_obs; //!< Observation
204 211 const GMatrix& m_rot; //!< Rotation matrix
205 212 };
  213 + */
206 214  
207 215 // Npred phi integration kernel for elliptical model
  216 + /*
208 217 class npred_elliptical_kern_phi : public GFunction {
209 218 public:
210 219 npred_elliptical_kern_phi(const GResponse& rsp,
... ... @@ -236,6 +245,7 @@ protected:
236 245 double m_cos_theta; //!< cosine of offset angle
237 246 const double& m_sin_theta; //!< Sine of offset angle
238 247 };
  248 + */
239 249  
240 250 };
241 251  
... ...
inst/com/include/GCOMResponse.hpp
  View file @ee42aa1
... ... @@ -38,6 +38,7 @@
38 38 /* __ Type definitions ___________________________________________________ */
39 39  
40 40 /* __ Forward declaration ________________________________________________ */
  41 +class GModelSky;
41 42  
42 43  
43 44 /***********************************************************************//**
... ... @@ -66,9 +67,6 @@ public:
66 67 virtual double irf(const GEvent& event,
67 68 const GPhoton& photon,
68 69 const GObservation& obs) const;
69   - virtual double convolve(const GModelSky& model,
70   - const GEvent& event,
71   - const GObservation& obs) const;
72 70 virtual double nroi(const GModelSky& model,
73 71 const GEnergy& obsEng,
74 72 const GTime& obsTime,
... ...
inst/com/pyext/GCOMResponse.i
  View file @ee42aa1
... ... @@ -51,8 +51,10 @@ public:
51 51 virtual double irf(const GEvent& event,
52 52 const GPhoton& photon,
53 53 const GObservation& obs) const;
54   - virtual double npred(const GPhoton& photon,
55   - const GObservation& obs) const;
  54 + virtual double nroi(const GModelSky& model,
  55 + const GEnergy& obsEng,
  56 + const GTime& obsTime,
  57 + const GObservation& obs) const;
56 58  
57 59 // Other Methods
58 60 void caldb(const GCaldb& caldb);
... ...
inst/cta/include/GCTAResponse.hpp
  View file @ee42aa1
1 1 /***************************************************************************
2 2 * GCTAResponse.hpp - CTA response abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -74,8 +74,10 @@ public:
74 74 virtual double irf(const GEvent& event,
75 75 const GPhoton& photon,
76 76 const GObservation& obs) const = 0;
77   - virtual double npred(const GPhoton& photon,
78   - const GObservation& obs) const = 0;
  77 + virtual double nroi(const GModelSky& model,
  78 + const GEnergy& obsEng,
  79 + const GTime& obsTime,
  80 + const GObservation& obs) const = 0;
79 81 virtual void read(const GXmlElement& xml) = 0;
80 82 virtual void write(GXmlElement& xml) const = 0;
81 83 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
... ...
inst/cta/include/GCTAResponseCube.hpp
  View file @ee42aa1
... ... @@ -80,21 +80,14 @@ public:
80 80 virtual double irf(const GEvent& event,
81 81 const GPhoton& photon,
82 82 const GObservation& obs) const;
83   - virtual double npred(const GPhoton& photon,
84   - const GObservation& obs) const;
  83 + virtual double nroi(const GModelSky& model,
  84 + const GEnergy& obsEng,
  85 + const GTime& obsTime,
  86 + const GObservation& obs) const;
85 87 virtual void read(const GXmlElement& xml);
86 88 virtual void write(GXmlElement& xml) const;
87 89 virtual std::string print(const GChatter& chatter = NORMAL) const;
88 90  
89   - // New pure virtual methods
90   - virtual double convolve(const GModelSky& model,
91   - const GEvent& event,
92   - const GObservation& obs) const;
93   - virtual double nroi(const GModelSky& model,
94   - const GEnergy& obsEng,
95   - const GTime& obsTime,
96   - const GObservation& obs) const;
97   -
98 91 // Overload base class methods
99 92 virtual double irf(const GEvent& event,
100 93 const GSource& source,
... ...
inst/cta/include/GCTAResponseIrf.hpp
  View file @ee42aa1
1 1 /***************************************************************************
2 2 * GCTAResponseIrf.hpp - CTA instrument response function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -84,20 +84,53 @@ public:
84 84 virtual double irf(const GEvent& event,
85 85 const GPhoton& photon,
86 86 const GObservation& obs) const;
  87 + virtual double nroi(const GModelSky& model,
  88 + const GEnergy& obsEng,
  89 + const GTime& obsTime,
  90 + const GObservation& obs) const;
  91 + /*
87 92 virtual double npred(const GPhoton& photon,
88 93 const GObservation& obs) const;
  94 + */
89 95 virtual void read(const GXmlElement& xml);
90 96 virtual void write(GXmlElement& xml) const;
91 97 virtual std::string print(const GChatter& chatter = NORMAL) const;
92 98  
93   - // New pure virtual methods
94   - virtual double convolve(const GModelSky& model,
95   - const GEvent& event,
96   - const GObservation& obs) const;
97   - virtual double nroi(const GModelSky& model,
98   - const GEnergy& obsEng,
99   - const GTime& obsTime,
100   - const GObservation& obs) const;
  99 + // New support methods
  100 + double nroi(const GModelSky& model,
  101 + const GEnergy& srcEng,
  102 + const GTime& srcTime,
  103 + const GEnergy& obsEng,
  104 + const GTime& obsTime,
  105 + const GObservation& obs) const;
  106 + double nroi_ptsrc(const GModelSky& model,
  107 + const GEnergy& srcEng,
  108 + const GTime& srcTime,
  109 + const GEnergy& obsEng,
  110 + const GTime& obsTime,
  111 + const GObservation& obs) const;
  112 + double nroi_radial(const GModelSky& model,
  113 + const GEnergy& srcEng,
  114 + const GTime& srcTime,
  115 + const GEnergy& obsEng,
  116 + const GTime& obsTime,
  117 + const GObservation& obs) const;
  118 + double nroi_elliptical(const GModelSky& model,
  119 + const GEnergy& srcEng,
  120 + const GTime& srcTime,
  121 + const GEnergy& obsEng,
  122 + const GTime& obsTime,
  123 + const GObservation& obs) const;
  124 + double nroi_diffuse(const GModelSky& model,
  125 + const GEnergy& srcEng,
  126 + const GTime& srcTime,
  127 + const GEnergy& obsEng,
  128 + const GTime& obsTime,
  129 + const GObservation& obs) const;
  130 + double nroi(const GPhoton& photon,
  131 + const GEnergy& obsEng,
  132 + const GTime& obsTime,
  133 + const GObservation& obs) const;
101 134  
102 135 // Overload virtual base class methods
103 136 virtual double irf_radial(const GEvent& event,
... ... @@ -109,12 +142,14 @@ public:
109 142 virtual double irf_diffuse(const GEvent& event,
110 143 const GSource& source,
111 144 const GObservation& obs) const;
  145 + /*
112 146 virtual double npred_radial(const GSource& source,
113 147 const GObservation& obs) const;
114 148 virtual double npred_elliptical(const GSource& source,
115 149 const GObservation& obs) const;
116 150 virtual double npred_diffuse(const GSource& source,
117 151 const GObservation& obs) const;
  152 + */
118 153 virtual GEbounds ebounds_src(const GEnergy& obsEnergy) const;
119 154  
120 155 // Other Methods
... ...
inst/cta/pyext/GCTAResponse.i
  View file @ee42aa1
... ... @@ -50,8 +50,10 @@ public:
50 50 virtual double irf(const GEvent& event,
51 51 const GPhoton& photon,
52 52 const GObservation& obs) const = 0;
53   - virtual double npred(const GPhoton& photon,
54   - const GObservation& obs) const = 0;
  53 + virtual double nroi(const GModelSky& model,
  54 + const GEnergy& obsEng,
  55 + const GTime& obsTime,
  56 + const GObservation& obs) const = 0;
55 57 virtual void read(const GXmlElement& xml) = 0;
56 58 virtual void write(GXmlElement& xml) const = 0;
57 59 };
... ...
inst/cta/pyext/GCTAResponseCube.i
  View file @ee42aa1
... ... @@ -54,11 +54,10 @@ public:
54 54 virtual double irf(const GEvent& event,
55 55 const GPhoton& photon,
56 56 const GObservation& obs) const;
57   - virtual double irf(const GEvent& event,
58   - const GSource& source,
59   - const GObservation& obs) const;
60   - virtual double npred(const GPhoton& photon,
61   - const GObservation& obs) const;
  57 + virtual double nroi(const GModelSky& model,
  58 + const GEnergy& obsEng,
  59 + const GTime& obsTime,
  60 + const GObservation& obs) const;
62 61 virtual void read(const GXmlElement& xml);
63 62 virtual void write(GXmlElement& xml) const;
64 63  
... ...
inst/cta/pyext/GCTAResponseIrf.i
  View file @ee42aa1
... ... @@ -108,8 +108,10 @@ public:
108 108 virtual double irf(const GEvent& event,
109 109 const GPhoton& photon,
110 110 const GObservation& obs) const;
111   - virtual double npred(const GPhoton& photon,
112   - const GObservation& obs) const;
  111 + virtual double nroi(const GModelSky& model,
  112 + const GEnergy& obsEng,
  113 + const GTime& obsTime,
  114 + const GObservation& obs) const;
113 115 virtual void read(const GXmlElement& xml);
114 116 virtual void write(GXmlElement& xml) const;
115 117  
... ... @@ -123,12 +125,6 @@ public:
123 125 virtual double irf_diffuse(const GEvent& event,
124 126 const GSource& source,
125 127 const GObservation& obs) const;
126   - virtual double npred_radial(const GSource& source,
127   - const GObservation& obs) const;
128   - virtual double npred_elliptical(const GSource& source,
129   - const GObservation& obs) const;
130   - virtual double npred_diffuse(const GSource& source,
131   - const GObservation& obs) const;
132 128 virtual GEbounds ebounds_src(const GEnergy& obsEnergy) const;
133 129  
134 130 // Other Methods
... ...
inst/cta/src/GCTAResponse.cpp
  View file @ee42aa1
1 1 /***************************************************************************
2 2 * GCTAResponse.cpp - CTA Response class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
inst/cta/src/GCTAResponseCube.cpp
  View file @ee42aa1
... ... @@ -57,7 +57,8 @@
57 57 " GObservation&)"
58 58 #define G_IRF_DIFFUSE "GCTAResponseCube::irf_diffuse(GEvent&, GSource&,"\
59 59 " GObservation&)"
60   -#define G_NPRED "GCTAResponseCube::npred(GPhoton&, GObservation&)"
  60 +#define G_NROI "GCTAResponseCube::nroi(GModelSky&, GEnergy&, GTime&, "\
  61 + "GObservation&)"
61 62 #define G_READ "GCTAResponseCube::read(GXmlElement&)"
62 63 #define G_WRITE "GCTAResponseCube::write(GXmlElement&)"
63 64  
... ... @@ -696,27 +697,45 @@ double GCTAResponseCube::irf_diffuse(const GEvent&amp; event,
696 697  
697 698  
698 699 /***********************************************************************//**
699   - * @brief Return spatial integral of point spread function
  700 + * @brief Return integral of event probability for a given sky model over ROI
700 701 *
701   - * @param[in] photon Incident photon.
  702 + * @param[in] model Incident photon.
702 703 * @param[in] obs Observation.
703   - * @return Spatial integral of point spread function.
  704 + * @return 0.0
704 705 *
705 706 * @exception GException::feature_not_implemented
706   - * Method not implemented.
  707 + * Method is not implemented.
  708 + *
  709 + * Computes the integral
  710 + *
  711 + * \f[
  712 + * N_{\rm ROI}(E',t') = \int_{\rm ROI} P(p',E',t') dp'
  713 + * \f]
  714 + *
  715 + * of the event probability
  716 + *
  717 + * \f[
  718 + * P(p',E',t') = \int \int \int
  719 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  720 + * \f]
  721 + *
  722 + * for a given sky model \f$S(p,E,t)\f$ and response function
  723 + * \f$R(p',E',t'|p,E,t)\f$ over the Region of Interest (ROI).
707 724 *
708   - * This method is a dummy method that is required for a class derived from
709   - * GResponse but that needs not to be implemented for a binned response.
  725 + * @todo Implement method (is maybe not really needed)
710 726 ***************************************************************************/
711   -double GCTAResponseCube::npred(const GPhoton& photon,
712   - const GObservation& obs) const
  727 +double GCTAResponseCube::nroi(const GModelSky& model,
  728 + const GEnergy& obsEng,
  729 + const GTime& obsTime,
  730 + const GObservation& obs) const
713 731 {
714   - // Feature not yet implemented
715   - throw GException::feature_not_implemented(G_NPRED,
716   - "Npred computation not implemented for cube analysis.");
  732 + // Method is not implemented
  733 + std::string msg = "Spatial integration of sky model over the data space "
  734 + "is not implemented.";
  735 + throw GException::feature_not_implemented(G_NROI, msg);
717 736  
718 737 // Return Npred
719   - return 0.0;
  738 + return (0.0);
720 739 }
721 740  
722 741  
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @ee42aa1
... ... @@ -36,6 +36,8 @@
36 36 #include "GMath.hpp"
37 37 #include "GIntegral.hpp"
38 38 #include "GCaldb.hpp"
  39 +#include "GModelSky.hpp"
  40 +#include "GModelSpatialPointSource.hpp"
39 41 #include "GModelSpatialRadial.hpp"
40 42 #include "GModelSpatialRadialShell.hpp"
41 43 #include "GModelSpatialElliptical.hpp"
... ... @@ -396,6 +398,843 @@ double GCTAResponseIrf::irf(const GEvent&amp; event,
396 398  
397 399  
398 400 /***********************************************************************//**
  401 + * @brief Return integral of event probability for a given sky model over ROI
  402 + *
  403 + * @param[in] model Incident photon.
  404 + * @param[in] obs Observation.
  405 + * @return 0.0
  406 + *
  407 + * @exception GException::feature_not_implemented
  408 + * Method is not implemented.
  409 + *
  410 + * Computes the integral
  411 + *
  412 + * \f[
  413 + * N_{\rm ROI}(E',t') = \int_{\rm ROI} P(p',E',t') dp'
  414 + * \f]
  415 + *
  416 + * of the event probability
  417 + *
  418 + * \f[
  419 + * P(p',E',t') = \int \int \int
  420 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  421 + * \f]
  422 + *
  423 + * for a given sky model \f$S(p,E,t)\f$ and response function
  424 + * \f$R(p',E',t'|p,E,t)\f$ over the Region of Interest (ROI).
  425 + ***************************************************************************/
  426 +double GCTAResponseIrf::nroi(const GModelSky& model,
  427 + const GEnergy& obsEng,
  428 + const GTime& obsTime,
  429 + const GObservation& obs) const
  430 +{
  431 + // Initialise Nroi value
  432 + double nroi = 0.0;
  433 +
  434 + // No time dispersion supported
  435 + const GTime& srcTime = obsTime;
  436 +
  437 + // If energy dispersion is requested then integrate over the relevant
  438 + // true photon energies ...
  439 + if (use_edisp()) {
  440 +
  441 + // Retrieve true energy boundaries
  442 + double obsLogEng = obsEng.log10TeV();
  443 + GEbounds ebounds = edisp()->ebounds_src(obsLogEng);
  444 +
  445 + // Loop over all boundaries
  446 + for (int i = 0; i < ebounds.size(); ++i) {
  447 +
  448 + // Get boundaries in MeV
  449 + double emin = ebounds.emin(i).MeV();
  450 + double emax = ebounds.emax(i).MeV();
  451 +
  452 + // Continue only if valid
  453 + if (emax > emin) {
  454 +
  455 + // Setup integration function
  456 + cta_nroi_kern integrand(model, *this, srcTime, obsEng, obsTime, obs);
  457 + GIntegral integral(&integrand);
  458 +
  459 + // Set integration precision
  460 + integral.eps(1.0e-3);
  461 +
  462 + // Do Romberg integration
  463 + emin = std::log(emin);
  464 + emax = std::log(emax);
  465 + nroi += integral.romberg(emin, emax);
  466 +
  467 + } // endif: interval was valid
  468 +
  469 + } // endfor: looped over energy intervals
  470 +
  471 + } // endif: energy dispersion requested
  472 +
  473 + // ... otherwise evaluate
  474 + else {
  475 +
  476 + // No energy dispersion
  477 + const GEnergy& srcEng = obsEng;
  478 +
  479 + // Compute response components
  480 + double npred_spatial = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
  481 + double npred_spectral = model.spectral()->eval(srcEng, srcTime);
  482 + double npred_temporal = model.temporal()->eval(srcTime);
  483 +
  484 + // Compute response
  485 + nroi = npred_spatial * npred_spectral * npred_temporal;
  486 +
  487 + } // endelse: no energy dispersion requested
  488 +
  489 + // If required, apply instrument specific model scaling
  490 + if (model.has_scales()) {
  491 + nroi *= model.scale(obs.instrument()).value();
  492 + }
  493 +
  494 + // Compile option: Check for NaN/Inf
  495 + #if defined(G_NAN_CHECK)
  496 + if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
  497 + std::cout << "*** ERROR: GCTAResponseIrf::nroi:";
  498 + std::cout << " NaN/Inf encountered";
  499 + std::cout << " (nroi=" << nroi;
  500 + std::cout << ", obsEng=" << obsEng;
  501 + std::cout << ", obsTime=" << obsTime;
  502 + std::cout << ")" << std::endl;
  503 + }
  504 + #endif
  505 +
  506 + // Return response value
  507 + return nroi;
  508 +}
  509 +double GCTAResponseIrf::nroi(const GModelSky& model,
  510 + const GEnergy& srcEng,
  511 + const GTime& srcTime,
  512 + const GEnergy& obsEng,
  513 + const GTime& obsTime,
  514 + const GObservation& obs) const
  515 +{
  516 + // Initialise response value
  517 + double nroi = 0.0;
  518 +
  519 + // Select method depending on the spatial model type
  520 + switch (model.spatial()->code()) {
  521 + case GMODEL_SPATIAL_POINT_SOURCE:
  522 + nroi = nroi_ptsrc(model, srcEng, srcTime, obsEng, obsTime, obs);
  523 + break;
  524 + case GMODEL_SPATIAL_RADIAL:
  525 + nroi = nroi_radial(model, srcEng, srcTime, obsEng, obsTime, obs);
  526 + break;
  527 + case GMODEL_SPATIAL_ELLIPTICAL:
  528 + nroi = nroi_elliptical(model, srcEng, srcTime, obsEng, obsTime, obs);
  529 + break;
  530 + case GMODEL_SPATIAL_DIFFUSE:
  531 + nroi = nroi_diffuse(model, srcEng, srcTime, obsEng, obsTime, obs);
  532 + break;
  533 + default:
  534 + break;
  535 + }
  536 +
  537 + // Return response value
  538 + return nroi;
  539 +}
  540 +double GCTAResponseIrf::nroi_ptsrc(const GModelSky& model,
  541 + const GEnergy& srcEng,
  542 + const GTime& srcTime,
  543 + const GEnergy& obsEng,
  544 + const GTime& obsTime,
  545 + const GObservation& obs) const
  546 +{
  547 + // Get point source spatial model
  548 + const GModelSpatialPointSource* src =
  549 + static_cast<const GModelSpatialPointSource*>(model.spatial());
  550 +
  551 + // Set Photon
  552 + GPhoton photon(src->dir(), srcEng, srcTime);
  553 +
  554 + // Compute Nroi
  555 + double nroi = this->nroi(photon, obsEng, obsTime, obs);
  556 +
  557 + // Return Nroi
  558 + return nroi;
  559 +}
  560 +
  561 +/***********************************************************************//**
  562 + * @brief Return spatial integral of radial source model over ROI
  563 + *
  564 + * @param[in] source Source.
  565 + * @param[in] obs Observation.
  566 + *
  567 + * @exception GCTAException::bad_model_type
  568 + * Model is not a radial model.
  569 + *
  570 + * Integrates the product of the radial model and Npred over the Region Of
  571 + * Interest using
  572 + *
  573 + * \f[
  574 + * \int_{\rho_{\rm min}}^{\rho_{\rm max}}
  575 + * \sin \rho \times S_{\rm p}(\rho | E, t) \times
  576 + * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
  577 + * N_{\rm pred}(\rho, \omega) d\omega
  578 + * d\rho
  579 + * \f]
  580 + *
  581 + * where
  582 + * - \f$S_{\rm p}(\rho | E, t)\f$ is the radial model,
  583 + * - \f$N_{\rm pred}(\rho,\omega)\f$ is the data space integral of the
  584 + * Instrument Response Function for a point spread function over the
  585 + * Region Of Interest,
  586 + * - \f$\rho\f$ is the distance from the model centre, and
  587 + * - \f$\omega\f$ is the azimuth angle is the position angle with respect to
  588 + * the connecting line between the model centre and the observed photon
  589 + * arrival direction.
  590 + *
  591 + * The integration is performed in a spherical coordinate system that is
  592 + * centred on the source model centre \f$\vec{m}\f$, with \f$(\rho,\omega)\f$
  593 + * being the zenith and azimuth angles, respectively.
  594 + *
  595 + * The zenith angle integration range \f$[\rho_{\rm min}, \rho_{\rm max}\f$
  596 + * and azimuth angle integration range
  597 + * \f$[\omega_{\rm min}, \omega_{\rm max}\f$
  598 + * are adjusted so that only coordinates within the circular region of
  599 + * interest will be considered.
  600 + *
  601 + * Note that we estimate the integration radius based on the size of the
  602 + * onaxis PSF. This should be fine as long as the offaxis PSF is not
  603 + * considerably larger than the onaxis PSF. We should verify this, however.
  604 + *
  605 + * @todo Verify that offaxis PSF is not considerably larger than onaxis
  606 + * PSF.
  607 + ***************************************************************************/
  608 +double GCTAResponseIrf::nroi_radial(const GModelSky& model,
  609 + const GEnergy& srcEng,
  610 + const GTime& srcTime,
  611 + const GEnergy& obsEng,
  612 + const GTime& obsTime,
  613 + const GObservation& obs) const
  614 +{
  615 + // Set number of iterations for Romberg integration.
  616 + // These values have been determined after careful testing, see
  617 + // https://cta-redmine.irap.omp.eu/issues/1299
  618 + static const int iter_rho = 6;
  619 + static const int iter_phi = 6;
  620 +
  621 + // Initialise Npred value
  622 + double npred = 0.0;
  623 +
  624 + // Retrieve CTA observation, ROI and pointing
  625 + const GCTAObservation& cta = retrieve_obs(G_NPRED_RADIAL, obs);
  626 + const GCTARoi& roi = retrieve_roi(G_NPRED_RADIAL, obs);
  627 + const GCTAPointing& pnt = cta.pointing();
  628 +
  629 + // Get pointer on radial model
  630 + const GModelSpatialRadial* radial =
  631 + dynamic_cast<const GModelSpatialRadial*>(&model);
  632 + if (radial == NULL) {
  633 + throw GCTAException::bad_model_type(G_NPRED_RADIAL);
  634 + }
  635 +
  636 + // Get source attributes
  637 + const GSkyDir& centre = radial->dir();
  638 + //const GEnergy& srcEng = source.energy();
  639 + //const GTime& srcTime = source.time();
  640 +
  641 + // Get pointing direction zenith angle and azimuth [radians]
  642 + double zenith = pnt.zenith();
  643 + double azimuth = pnt.azimuth();
  644 +
  645 + // Get log10(E/TeV) of true photon energy
  646 + double srcLogEng = srcEng.log10TeV();
  647 +
  648 + // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
  649 + // as this allows us doing this computation in the outer loop. This
  650 + // should be sufficient here, unless the offaxis PSF becomes much worse
  651 + // than the onaxis PSF. In this case, we may add a safety factor here
  652 + // to make sure we encompass the entire PSF.
  653 + double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
  654 +
  655 + // Extract ROI radius (radians)
  656 + double roi_radius = roi.radius() * gammalib::deg2rad;
  657 +
  658 + // Compute distance between ROI and model centre (radians)
  659 + double roi_model_distance = roi.centre().dir().dist(centre);
  660 +
  661 + // Compute the ROI radius plus maximum PSF radius (radians). Any photon
  662 + // coming from beyond this radius will not make it in the dataspace and
  663 + // thus can be neglected.
  664 + double roi_psf_radius = roi_radius + psf_max_radius;
  665 +
  666 + // Set offset angle integration range. We take here the ROI+PSF into
  667 + // account to make no integrations beyond the point where the
  668 + // contribution drops to zero.
  669 + double rho_min = (roi_model_distance > roi_psf_radius)
  670 + ? roi_model_distance - roi_psf_radius: 0.0;
  671 + double rho_max = radial->theta_max();
  672 +
  673 + // Perform offset angle integration only if interval is valid
  674 + if (rho_max > rho_min) {
  675 +
  676 + // Compute rotation matrix to convert from native model coordinates,
  677 + // given by (rho,omega), into celestial coordinates.
  678 + GMatrix ry;
  679 + GMatrix rz;
  680 + ry.eulery(radial->dec() - 90.0);
  681 + rz.eulerz(-radial->ra());
  682 + GMatrix rot = (ry * rz).transpose();
  683 +
  684 + // Compute position angle of ROI centre with respect to model
  685 + // centre (radians)
  686 + double omega0 = centre.posang(roi.centre().dir());
  687 +
  688 + // Setup integration kernel
  689 + cta_npred_radial_kern_rho integrand(*this,
  690 + *radial,
  691 + srcEng,
  692 + srcTime,
  693 + obsEng,
  694 + obsTime,
  695 + cta,
  696 + rot,
  697 + roi_model_distance,
  698 + roi_psf_radius,
  699 + omega0,
  700 + iter_phi);
  701 +
  702 + // Integrate over model's zenith angle
  703 + GIntegral integral(&integrand);
  704 + integral.fixed_iter(iter_rho);
  705 +
  706 + // Setup integration boundaries
  707 + std::vector<double> bounds;
  708 + bounds.push_back(rho_min);
  709 + bounds.push_back(rho_max);
  710 +
  711 + // If the integration range includes a transition between full
  712 + // containment of model within ROI and partial containment, then
  713 + // add a boundary at this location
  714 + double transition_point = roi_psf_radius - roi_model_distance;
  715 + if (transition_point > rho_min && transition_point < rho_max) {
  716 + bounds.push_back(transition_point);
  717 + }
  718 +
  719 + // If the integration range includes a transition between full
  720 + // containment of ROI within model and partial containment, then
  721 + // add a boundary at this location
  722 + transition_point = roi_psf_radius + roi_model_distance;
  723 + if (transition_point > rho_min && transition_point < rho_max) {
  724 + bounds.push_back(transition_point);
  725 + }
  726 +
  727 + // If we have a shell model then add an integration boundary for the
  728 + // shell radius as a function discontinuity will occur at this
  729 + // location
  730 + const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(radial);
  731 + if (shell != NULL) {
  732 + double shell_radius = shell->radius() * gammalib::deg2rad;
  733 + if (shell_radius > rho_min && shell_radius < rho_max) {
  734 + bounds.push_back(shell_radius);
  735 + }
  736 + }
  737 +
  738 + // Integrate kernel
  739 + npred = integral.romberg(bounds, iter_rho);
  740 +
  741 + // Compile option: Show integration results
  742 + #if defined(G_DEBUG_NPRED_RADIAL)
  743 + std::cout << "GCTAResponseIrf::npred_radial:";
  744 + std::cout << " rho_min=" << rho_min;
  745 + std::cout << " rho_max=" << rho_max;
  746 + std::cout << " npred=" << npred << std::endl;
  747 + #endif
  748 +
  749 + } // endif: offset angle range was valid
  750 +
  751 + // Debug: Check for NaN
  752 + #if defined(G_NAN_CHECK)
  753 + if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
  754 + std::cout << "*** ERROR: GCTAResponseIrf::npred_radial:";
  755 + std::cout << " NaN/Inf encountered";
  756 + std::cout << " (npred=" << npred;
  757 + std::cout << ", rho_min=" << rho_min;
  758 + std::cout << ", rho_max=" << rho_max;
  759 + std::cout << ")" << std::endl;
  760 + }
  761 + #endif
  762 +
  763 + // Return Npred
  764 + return npred;
  765 +}
  766 +
  767 +
  768 +/***********************************************************************//**
  769 + * @brief Return spatial integral of elliptical source model over ROI
  770 + *
  771 + * @param[in] source Source.
  772 + * @param[in] obs Observation.
  773 + *
  774 + * @exception GCTAException::bad_model_type
  775 + * Model is not an elliptical model.
  776 + *
  777 + * Integrates the product of the elliptical model and Npred over the Region
  778 + * Of Interest using
  779 + *
  780 + * \f[
  781 + * \int_{\rho_{\rm min}}^{\rho_{\rm max}}
  782 + * \sin \rho \times
  783 + * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
  784 + * S_{\rm p}(\rho,\omega | E, t) \,
  785 + * N_{\rm pred}(\rho,\omega) d\omega
  786 + * d\rho
  787 + * \f]
  788 + *
  789 + * where
  790 + * - \f$S_{\rm p}(\rho,\omega | E, t)\f$ is the elliptical model,
  791 + * - \f$N_{\rm pred}(\rho,\omega)\f$ is the data space integral of the
  792 + * Instrument Response Function for a point spread function over the
  793 + * Region Of Interest,
  794 + * - \f$\rho\f$ is the distance from the model centre, and
  795 + * - \f$\omega\f$ is the azimuth angle is the position angle with respect to
  796 + * the connecting line between the model centre and the observed photon
  797 + * arrival direction.
  798 + *
  799 + * The integration is performed in a spherical coordinate system that is
  800 + * centred on the source model centre \f$\vec{m}\f$, with \f$(\rho,\omega)\f$
  801 + * being the zenith and azimuth angles, respectively.
  802 + *
  803 + * The zenith angle integration range \f$[\rho_{\rm min}, \rho_{\rm max}\f$
  804 + * and azimuth angle integration range
  805 + * \f$[\omega_{\rm min}, \omega_{\rm max}\f$
  806 + * are adjusted so that only coordinates within the circular region of
  807 + * interest will be considered.
  808 + *
  809 + * Note that we estimate the integration radius based on the size of the
  810 + * onaxis PSF. This should be fine as long as the offaxis PSF is not
  811 + * considerably larger than the onaxis PSF. We should verify this, however.
  812 + *
  813 + * @todo Verify that offaxis PSF is not considerably larger than onaxis
  814 + * PSF.
  815 + ***************************************************************************/
  816 +double GCTAResponseIrf::nroi_elliptical(const GModelSky& model,
  817 + const GEnergy& srcEng,
  818 + const GTime& srcTime,
  819 + const GEnergy& obsEng,
  820 + const GTime& obsTime,
  821 + const GObservation& obs) const
  822 +{
  823 + // Set number of iterations for Romberg integration.
  824 + // These values have been determined after careful testing, see
  825 + // https://cta-redmine.irap.omp.eu/issues/1299
  826 + static const int iter_rho = 6;
  827 + static const int iter_phi = 6;
  828 +
  829 + // Initialise Npred value
  830 + double npred = 0.0;
  831 +
  832 + // Retrieve CTA observation, ROI and pointing
  833 + const GCTAObservation& cta = retrieve_obs(G_NPRED_ELLIPTICAL, obs);
  834 + const GCTARoi& roi = retrieve_roi(G_NPRED_ELLIPTICAL, obs);
  835 + const GCTAPointing& pnt = cta.pointing();
  836 +
  837 + // Get pointer on elliptical model
  838 + const GModelSpatialElliptical* elliptical =
  839 + dynamic_cast<const GModelSpatialElliptical*>(&model);
  840 + if (elliptical == NULL) {
  841 + throw GCTAException::bad_model_type(G_NPRED_ELLIPTICAL);
  842 + }
  843 +
  844 + // Get source attributes
  845 + const GSkyDir& centre = elliptical->dir();
  846 + //const GEnergy& srcEng = source.energy();
  847 + //const GTime& srcTime = source.time();
  848 +
  849 + // Get pointing direction zenith angle and azimuth [radians]
  850 + double zenith = pnt.zenith();
  851 + double azimuth = pnt.azimuth();
  852 +
  853 + // Get log10(E/TeV) of true photon energy
  854 + double srcLogEng = srcEng.log10TeV();
  855 +
  856 + // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
  857 + // as this allows us doing this computation in the outer loop. This
  858 + // should be sufficient here, unless the offaxis PSF becomes much worse
  859 + // than the onaxis PSF. In this case, we may add a safety factor here
  860 + // to make sure we encompass the entire PSF.
  861 + double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
  862 +
  863 + // Extract ROI radius plus maximum PSF radius (radians). Any photon
  864 + // coming from beyond this radius will not make it in the dataspace and
  865 + // thus can be neglected.
  866 + double radius_roi = roi.radius() * gammalib::deg2rad + psf_max_radius;
  867 +
  868 + // Compute distance between ROI and model centre (radians)
  869 + double rho_roi = roi.centre().dir().dist(centre);
  870 +
  871 + // Get the ellipse boundary (radians). Note that these are NOT the
  872 + // parameters of the ellipse but of a boundary ellipse that is used
  873 + // for computing the relevant omega angle intervals for a given angle
  874 + // rho. The boundary ellipse takes care of the possibility that the
  875 + // semiminor axis is larger than the semimajor axis
  876 + double semimajor; // Will be the larger axis
  877 + double semiminor; // Will be the smaller axis
  878 + double posangle; // Will be the corrected position angle
  879 + double aspect_ratio; // Ratio between smaller/larger axis of model
  880 + if (elliptical->semimajor() >= elliptical->semiminor()) {
  881 + aspect_ratio = (elliptical->semimajor() > 0.0) ?
  882 + elliptical->semiminor() / elliptical->semimajor() : 0.0;
  883 + posangle = elliptical->posangle() * gammalib::deg2rad;
  884 + }
  885 + else {
  886 + aspect_ratio = (elliptical->semiminor() > 0.0) ?
  887 + elliptical->semimajor() / elliptical->semiminor() : 0.0;
  888 + posangle = elliptical->posangle() * gammalib::deg2rad + gammalib::pihalf;
  889 + }
  890 + semimajor = elliptical->theta_max();
  891 + semiminor = semimajor * aspect_ratio;
  892 +
  893 + // Set offset angle integration range. We take here the ROI+PSF into
  894 + // account to make no integrations beyond the point where the
  895 + // contribution drops to zero.
  896 + double rho_min = (rho_roi > radius_roi) ? rho_roi - radius_roi: 0.0;
  897 + double rho_max = rho_roi + radius_roi;
  898 + if (rho_max > semimajor) {
  899 + rho_max = semimajor;
  900 + }
  901 +
  902 + // Perform offset angle integration only if interval is valid
  903 + if (rho_max > rho_min) {
  904 +
  905 + // Compute rotation matrix to convert from native model coordinates,
  906 + // given by (rho,omega), into celestial coordinates.
  907 + GMatrix ry;
  908 + GMatrix rz;
  909 + ry.eulery(elliptical->dec() - 90.0);
  910 + rz.eulerz(-elliptical->ra());
  911 + GMatrix rot = (ry * rz).transpose();
  912 +
  913 + // Compute position angle of ROI centre with respect to model
  914 + // centre (radians)
  915 + double posangle_roi = centre.posang(roi.centre().dir());
  916 +
  917 + // Setup integration kernel
  918 + cta_npred_elliptical_kern_rho integrand(*this,
  919 + *elliptical,
  920 + semimajor,
  921 + semiminor,
  922 + posangle,
  923 + srcEng,
  924 + srcTime,
  925 + obsEng,
  926 + obsTime,
  927 + cta,
  928 + rot,
  929 + rho_roi,
  930 + posangle_roi,
  931 + radius_roi,
  932 + iter_phi);
  933 +
  934 + // Integrate over model's zenith angle
  935 + GIntegral integral(&integrand);
  936 + integral.fixed_iter(iter_rho);
  937 +
  938 + // Setup integration boundaries
  939 + std::vector<double> bounds;
  940 + bounds.push_back(rho_min);
  941 + bounds.push_back(rho_max);
  942 +
  943 + // If the integration range includes the semiminor boundary, then
  944 + // add an integration boundary at that location
  945 + if (semiminor > rho_min && semiminor < rho_max) {
  946 + bounds.push_back(semiminor);
  947 + }
  948 +
  949 + // Integrate kernel
  950 + npred = integral.romberg(bounds, iter_rho);
  951 +
  952 + // Compile option: Show integration results
  953 + #if defined(G_DEBUG_NPRED_ELLIPTICAL)
  954 + std::cout << "GCTAResponseIrf::npred_elliptical:";
  955 + std::cout << " rho_min=" << rho_min;
  956 + std::cout << " rho_max=" << rho_max;
  957 + std::cout << " npred=" << npred << std::endl;
  958 + #endif
  959 +
  960 + } // endif: offset angle range was valid
  961 +
  962 + // Debug: Check for NaN
  963 + #if defined(G_NAN_CHECK)
  964 + if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
  965 + std::cout << "*** ERROR: GCTAResponseIrf::npred_elliptical:";
  966 + std::cout << " NaN/Inf encountered";
  967 + std::cout << " (npred=" << npred;
  968 + std::cout << ", rho_min=" << rho_min;
  969 + std::cout << ", rho_max=" << rho_max;
  970 + std::cout << ")" << std::endl;
  971 + }
  972 + #endif
  973 +
  974 + // Return Npred
  975 + return npred;
  976 +}
  977 +
  978 +
  979 +/***********************************************************************//**
  980 + * @brief Return spatial integral of diffuse source model over ROI
  981 + *
  982 + * @param[in] source Source.
  983 + * @param[in] obs Observation.
  984 + *
  985 + * @exception GCTAException::bad_model_type
  986 + * Model is not a radial model.
  987 + *
  988 + * Integrates the product of the diffuse model and Npred over the Region
  989 + * Of Interest using
  990 + *
  991 + * \f[
  992 + * \int_{0}^{\theta_{\rm max}}
  993 + * \sin \theta \times
  994 + * \int_{0}^{2\pi}
  995 + * S_{\rm p}(\theta, \phi | E, t) \,
  996 + * N_{\rm pred}(\theta, \phi) d\phi
  997 + * d\theta
  998 + * \f]
  999 + *
  1000 + * where
  1001 + * - \f$S_{\rm p}(\theta, \phi | E, t)\f$ is the diffuse model,
  1002 + * - \f$N_{\rm pred}(\theta, \phi)\f$ is the data space integral of the
  1003 + * Instrument Response Function for a point spread function over the
  1004 + * Region Of Interest in the reference frame of the diffuse source
  1005 + * model
  1006 + * - \f$\theta\f$ is the distance from the ROI centre, and
  1007 + * - \f$\phi\f$ is the azimuth angle.
  1008 + *
  1009 + * Note that the integration precision was adjusted trading-off between
  1010 + * computation time and computation precision. A value of 1e-4 was judged
  1011 + * appropriate.
  1012 + *
  1013 + * Note that we estimate the integration radius based on the size of the
  1014 + * onaxis PSF in this method. This should be fine as long as the offaxis
  1015 + * PSF is not considerably larger than the onaxis PSF. We should verify
  1016 + * this, however.
  1017 + *
  1018 + * @todo Verify that offaxis PSF is not considerably larger than onaxis PSF.
  1019 + ***************************************************************************/
  1020 +double GCTAResponseIrf::nroi_diffuse(const GModelSky& model,
  1021 + const GEnergy& srcEng,
  1022 + const GTime& srcTime,
  1023 + const GEnergy& obsEng,
  1024 + const GTime& obsTime,
  1025 + const GObservation& obs) const
  1026 +{
  1027 + // Set number of iterations for Romberg integration.
  1028 + // These values have been determined after careful testing, see
  1029 + // https://cta-redmine.irap.omp.eu/issues/1248
  1030 + static const int iter_rho = 9;
  1031 + static const int iter_phi = 9;
  1032 +
  1033 + // Initialise Npred value
  1034 + double npred = 0.0;
  1035 + bool has_npred = false;
  1036 +
  1037 + // Build unique identifier
  1038 + std::string id = model.name() + "::" + obs.id();
  1039 +
  1040 + // Check if Npred value is already in cache
  1041 + #if defined(G_USE_NPRED_CACHE)
  1042 + if (!m_npred_names.empty()) {
  1043 +
  1044 + // Search for unique identifier, and if found, recover Npred value
  1045 + // and break
  1046 + for (int i = 0; i < m_npred_names.size(); ++i) {
  1047 + if (m_npred_names[i] == id &&
  1048 + m_npred_energies[i] == srcEng &&
  1049 + m_npred_times[i] == srcTime) {
  1050 + npred = m_npred_values[i];
  1051 + has_npred = true;
  1052 + #if defined(G_DEBUG_NPRED_DIFFUSE)
  1053 + std::cout << "GCTAResponseIrf::npred_diffuse:";
  1054 + std::cout << " cache=" << i;
  1055 + std::cout << " npred=" << npred << std::endl;
  1056 + #endif
  1057 + break;
  1058 + }
  1059 + }
  1060 +
  1061 + } // endif: there were values in the Npred cache
  1062 + #endif
  1063 +
  1064 + // Continue only if no Npred cache value was found
  1065 + if (!has_npred) {
  1066 +
  1067 + // Retrieve CTA observation, ROI and pointing
  1068 + const GCTAObservation& cta = retrieve_obs(G_NPRED_DIFFUSE, obs);
  1069 + const GCTARoi& roi = retrieve_roi(G_NPRED_DIFFUSE, obs);
  1070 + const GCTAPointing& pnt = cta.pointing();
  1071 +
  1072 + // Get pointer on spatial model
  1073 + const GModelSpatial* spatial =
  1074 + dynamic_cast<const GModelSpatial*>(model.spatial());
  1075 + if (spatial == NULL) {
  1076 + throw GCTAException::bad_model_type(G_NPRED_DIFFUSE);
  1077 + }
  1078 +
  1079 + // Get source attributes
  1080 + //const GEnergy& srcEng = source.energy();
  1081 + //const GTime& srcTime = source.time();
  1082 +
  1083 + // Get pointing direction zenith angle and azimuth [radians]
  1084 + double zenith = pnt.zenith();
  1085 + double azimuth = pnt.azimuth();
  1086 +
  1087 + // Get log10(E/TeV) of true photon energy
  1088 + double srcLogEng = srcEng.log10TeV();
  1089 +
  1090 + // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
  1091 + // as this allows us doing this computation in the outer loop. This
  1092 + // should be sufficient here, unless the offaxis PSF becomes much worse
  1093 + // than the onaxis PSF. In this case, we may add a safety factor here
  1094 + // to make sure we encompass the entire PSF.
  1095 + double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
  1096 +
  1097 + // Extract ROI radius (radians)
  1098 + double roi_radius = roi.radius() * gammalib::deg2rad;
  1099 +
  1100 + // Compute the ROI radius plus maximum PSF radius (radians). Any photon
  1101 + // coming from beyond this radius will not make it in the dataspace and
  1102 + // thus can be neglected.
  1103 + double roi_psf_radius = roi_radius + psf_max_radius;
  1104 +
  1105 + // Perform offset angle integration only if interval is valid
  1106 + if (roi_psf_radius > 0.0) {
  1107 +
  1108 + // Compute rotation matrix to convert from native ROI coordinates,
  1109 + // given by (theta,phi), into celestial coordinates.
  1110 + GMatrix ry;
  1111 + GMatrix rz;
  1112 + ry.eulery(roi.centre().dir().dec_deg() - 90.0);
  1113 + rz.eulerz(-roi.centre().dir().ra_deg());
  1114 + GMatrix rot = (ry * rz).transpose();
  1115 +
  1116 + // Setup integration kernel
  1117 + cta_npred_diffuse_kern_theta integrand(*this,
  1118 + *spatial,
  1119 + srcEng,
  1120 + srcTime,
  1121 + obsEng,
  1122 + obsTime,
  1123 + cta,
  1124 + rot,
  1125 + iter_phi);
  1126 +
  1127 +
  1128 + // Integrate over model's zenith angle
  1129 + GIntegral integral(&integrand);
  1130 + integral.fixed_iter(iter_rho);
  1131 + npred = integral.romberg(0.0, roi_psf_radius);
  1132 +
  1133 + // Compile option: Show integration results
  1134 + #if defined(G_DEBUG_NPRED_DIFFUSE)
  1135 + std::cout << "GCTAResponseIrf::npred_diffuse:";
  1136 + std::cout << " roi_psf_radius=" << roi_psf_radius;
  1137 + std::cout << " npred=" << npred;
  1138 + std::cout << " id=" << id << std::endl;
  1139 + #endif
  1140 +
  1141 + } // endif: offset angle range was valid
  1142 +
  1143 + // Store result in Npred cache
  1144 + #if defined(G_USE_NPRED_CACHE)
  1145 + m_npred_names.push_back(id);
  1146 + m_npred_energies.push_back(srcEng);
  1147 + m_npred_times.push_back(srcTime);
  1148 + m_npred_values.push_back(npred);
  1149 + #endif
  1150 +
  1151 + // Debug: Check for NaN
  1152 + #if defined(G_NAN_CHECK)
  1153 + if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
  1154 + std::cout << "*** ERROR: GCTAResponseIrf::npred_diffuse:";
  1155 + std::cout << " NaN/Inf encountered";
  1156 + std::cout << " (npred=" << npred;
  1157 + std::cout << ", roi_psf_radius=" << roi_psf_radius;
  1158 + std::cout << ")" << std::endl;
  1159 + }
  1160 + #endif
  1161 +
  1162 + } // endif: Npred computation required
  1163 +
  1164 + // Return Npred
  1165 + return npred;
  1166 +}
  1167 +
  1168 +
  1169 +double GCTAResponseIrf::nroi(const GPhoton& photon,
  1170 + const GEnergy& obsEng,
  1171 + const GTime& obsTime,
  1172 + const GObservation& obs) const
  1173 +{
  1174 + // Retrieve CTA observation, ROI and pointing
  1175 + const GCTAObservation& cta = retrieve_obs(G_NPRED, obs);
  1176 + const GCTARoi& roi = retrieve_roi(G_NPRED, obs);
  1177 + const GCTAPointing& pnt = cta.pointing();
  1178 +
  1179 + // Get photon attributes
  1180 + const GSkyDir& srcDir = photon.dir();
  1181 + const GEnergy& srcEng = photon.energy();
  1182 + const GTime& srcTime = photon.time();
  1183 +
  1184 + // Get pointing direction zenith angle and azimuth [radians]
  1185 + double zenith = pnt.zenith();
  1186 + double azimuth = pnt.azimuth();
  1187 +
  1188 + // Get radial offset and polar angles of true photon in camera [radians]
  1189 + double theta = pnt.dir().dist(srcDir);
  1190 + double phi = 0.0; //TODO: Implement Phi dependence
  1191 +
  1192 + // Get log10(E/TeV) of true photon energy.
  1193 + double srcLogEng = srcEng.log10TeV();
  1194 +
  1195 + // Get effectve area components
  1196 + double nroi = aeff(theta, phi, zenith, azimuth, srcLogEng);
  1197 +
  1198 + // Multiply-in PSF
  1199 + if (nroi > 0.0) {
  1200 +
  1201 + // Get PSF
  1202 + nroi *= npsf(srcDir, srcLogEng, srcTime, pnt, roi);
  1203 +
  1204 + // Multiply-in energy dispersion
  1205 + if (use_edisp() && nroi > 0.0) {
  1206 +
  1207 + // Get log10(E/TeV) of observed photon energy.
  1208 + double obsLogEng = obsEng.log10TeV();
  1209 +
  1210 + // Get energy dispersion
  1211 + nroi *= edisp()->operator()(obsLogEng, srcLogEng, theta, phi, zenith, azimuth);
  1212 +
  1213 + } // endif: had energy dispersion
  1214 +
  1215 + // Apply deadtime correction
  1216 + nroi *= obs.deadc(srcTime);
  1217 +
  1218 + } // endif: had non-zero effective area
  1219 +
  1220 + // Compile option: Check for NaN/Inf
  1221 + #if defined(G_NAN_CHECK)
  1222 + if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
  1223 + std::cout << "*** ERROR: GCTAResponseIrf::nroi:";
  1224 + std::cout << " NaN/Inf encountered";
  1225 + std::cout << " (nroi=" << nroi;
  1226 + std::cout << ", theta=" << theta;
  1227 + std::cout << ", phi=" << phi << ")";
  1228 + std::cout << std::endl;
  1229 + }
  1230 + #endif
  1231 +
  1232 + // Return Nroi
  1233 + return nroi;
  1234 +}
  1235 +
  1236 +
  1237 +/***********************************************************************//**
399 1238 * @brief Return spatial integral of point spread function
400 1239 *
401 1240 * @param[in] photon Incident photon.
... ... @@ -404,6 +1243,7 @@ double GCTAResponseIrf::irf(const GEvent&amp; event,
404 1243 * @todo Set polar angle phi of photon in camera system
405 1244 * @todo Write method documentation
406 1245 ***************************************************************************/
  1246 +/*
407 1247 double GCTAResponseIrf::npred(const GPhoton& photon,
408 1248 const GObservation& obs) const
409 1249 {
... ... @@ -465,6 +1305,7 @@ double GCTAResponseIrf::npred(const GPhoton&amp; photon,
465 1305 // Return Npred
466 1306 return npred;
467 1307 }
  1308 +*/
468 1309  
469 1310  
470 1311 /***********************************************************************//**
... ... @@ -2031,6 +2872,7 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp; event,
2031 2872 * @todo Verify that offaxis PSF is not considerably larger than onaxis
2032 2873 * PSF.
2033 2874 ***************************************************************************/
  2875 +/*
2034 2876 double GCTAResponseIrf::npred_radial(const GSource& source,
2035 2877 const GObservation& obs) const
2036 2878 {
... ... @@ -2183,7 +3025,7 @@ double GCTAResponseIrf::npred_radial(const GSource&amp; source,
2183 3025 // Return Npred
2184 3026 return npred;
2185 3027 }
2186   -
  3028 +*/
2187 3029  
2188 3030 /***********************************************************************//**
2189 3031 * @brief Return spatial integral of elliptical source model over ROI
... ... @@ -2233,6 +3075,7 @@ double GCTAResponseIrf::npred_radial(const GSource&amp; source,
2233 3075 * @todo Verify that offaxis PSF is not considerably larger than onaxis
2234 3076 * PSF.
2235 3077 ***************************************************************************/
  3078 +/*
2236 3079 double GCTAResponseIrf::npred_elliptical(const GSource& source,
2237 3080 const GObservation& obs) const
2238 3081 {
... ... @@ -2388,7 +3231,7 @@ double GCTAResponseIrf::npred_elliptical(const GSource&amp; source,
2388 3231 // Return Npred
2389 3232 return npred;
2390 3233 }
2391   -
  3234 +*/
2392 3235  
2393 3236 /***********************************************************************//**
2394 3237 * @brief Return spatial integral of diffuse source model over ROI
... ... @@ -2431,6 +3274,7 @@ double GCTAResponseIrf::npred_elliptical(const GSource&amp; source,
2431 3274 *
2432 3275 * @todo Verify that offaxis PSF is not considerably larger than onaxis PSF.
2433 3276 ***************************************************************************/
  3277 +/*
2434 3278 double GCTAResponseIrf::npred_diffuse(const GSource& source,
2435 3279 const GObservation& obs) const
2436 3280 {
... ... @@ -2572,7 +3416,7 @@ double GCTAResponseIrf::npred_diffuse(const GSource&amp; source,
2572 3416 // Return Npred
2573 3417 return npred;
2574 3418 }
2575   -
  3419 +*/
2576 3420  
2577 3421 /***********************************************************************//**
2578 3422 * @brief Return true energy boundaries for a specific observed energy
... ...
inst/cta/src/GCTAResponse_helpers.cpp
  View file @ee42aa1
... ... @@ -245,6 +245,32 @@ double cta_nedisp_kern::eval(const double&amp; logEobs)
245 245 }
246 246  
247 247  
  248 +/***********************************************************************//**
  249 + * @brief Integration kernel for nedisp() method
  250 + *
  251 + * @param[in] logEobs Logarithm of true photon energy.
  252 + * @return Nroi.
  253 + ***************************************************************************/
  254 +double cta_nroi_kern::eval(const double& logEsrc)
  255 +{
  256 + // Set true energy
  257 + GEnergy srcEng;
  258 + double expx = std::exp(logEsrc);
  259 + srcEng.MeV(expx);
  260 +
  261 + // Compute response components
  262 + double nroi_spatial = m_rsp.nroi(m_model, srcEng, m_srcTime, m_obsEng, m_obsTime, m_obs);
  263 + double nroi_spectral = m_model.spectral()->eval(srcEng, m_srcTime);
  264 + double nroi_temporal = m_model.temporal()->eval(m_srcTime);
  265 +
  266 + // Compute response
  267 + double nroi = nroi_spatial * nroi_spectral * nroi_temporal * expx;
  268 +
  269 + // Return response
  270 + return nroi;
  271 +}
  272 +
  273 +
248 274 /*==========================================================================
249 275 = =
250 276 = Helper class methods for radial models =
... ... @@ -521,6 +547,8 @@ double cta_npred_radial_kern_rho::eval(const double&amp; rho)
521 547 cta_npred_radial_kern_omega integrand(m_rsp,
522 548 m_srcEng,
523 549 m_srcTime,
  550 + m_obsEng,
  551 + m_obsTime,
524 552 m_obs,
525 553 m_rot,
526 554 sin_rho,
... ... @@ -584,7 +612,8 @@ double cta_npred_radial_kern_omega::eval(const double&amp; omega)
584 612 GPhoton photon(srcDir, m_srcEng, m_srcTime);
585 613  
586 614 // Compute point source Npred for this sky direction
587   - double npred = m_rsp.npred(photon, m_obs);
  615 + //double npred = m_rsp.npred(photon, m_obs);
  616 + double npred = m_rsp.nroi(photon, m_obsEng, m_obsTime, m_obs);
588 617  
589 618 // Debug: Check for NaN
590 619 #if defined(G_NAN_CHECK)
... ... @@ -978,6 +1007,8 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
978 1007 m_model,
979 1008 m_srcEng,
980 1009 m_srcTime,
  1010 + m_obsEng,
  1011 + m_obsTime,
981 1012 m_obs,
982 1013 m_rot,
983 1014 rho_kluge,
... ... @@ -1159,7 +1190,8 @@ double cta_npred_elliptical_kern_omega::eval(const double&amp; omega)
1159 1190 GPhoton photon(srcDir, m_srcEng, m_srcTime);
1160 1191  
1161 1192 // Compute Npred for this sky direction
1162   - npred = m_rsp.npred(photon, m_obs) * model;
  1193 + //npred = m_rsp.npred(photon, m_obs) * model;
  1194 + npred = m_rsp.nroi(photon, m_obsEng, m_obsTime, m_obs) * model;
1163 1195  
1164 1196 // Debug: Check for NaN
1165 1197 #if defined(G_NAN_CHECK)
... ... @@ -1421,6 +1453,8 @@ double cta_npred_diffuse_kern_theta::eval(const double&amp; theta)
1421 1453 m_model,
1422 1454 m_srcEng,
1423 1455 m_srcTime,
  1456 + m_obsEng,
  1457 + m_obsTime,
1424 1458 m_obs,
1425 1459 m_rot,
1426 1460 theta,
... ... @@ -1499,7 +1533,8 @@ double cta_npred_diffuse_kern_phi::eval(const double&amp; phi)
1499 1533 if (intensity > 0.0) {
1500 1534  
1501 1535 // Compute Npred for this sky direction
1502   - npred = m_rsp.npred(photon, m_obs) * intensity;
  1536 + //npred = m_rsp.npred(photon, m_obs) * intensity;
  1537 + npred = m_rsp.nroi(photon, m_obsEng, m_obsTime, m_obs) * intensity;
1503 1538  
1504 1539 // Debug: Check for NaN
1505 1540 #if defined(G_NAN_CHECK)
... ...
inst/cta/src/GCTAResponse_helpers.hpp
  View file @ee42aa1
1 1 /***************************************************************************
2 2 * GCTAResponse_helpers.hpp - CTA response helper classes *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -36,6 +36,7 @@
36 36 #include "GMatrix.hpp"
37 37 #include "GEnergy.hpp"
38 38 #include "GTime.hpp"
  39 +#include "GModelSky.hpp"
39 40 #include "GModelSpatial.hpp"
40 41 #include "GModelSpatialRadial.hpp"
41 42 #include "GCTAResponseCube.hpp"
... ... @@ -144,6 +145,41 @@ protected:
144 145  
145 146  
146 147 /***********************************************************************//**
  148 + * @class cta_nedisp_kern
  149 + *
  150 + * @brief Integration kernel for nedisp() method
  151 + *
  152 + * This class implements the integration kernel for the nedisp() method.
  153 + * The cta_nedisp_kern::eval method evaluates the energy dispersion for a
  154 + * given log10 of true photon energy as function of the log10 of observed
  155 + * event energy.
  156 + ***************************************************************************/
  157 +class cta_nroi_kern : public GFunction {
  158 +public:
  159 + cta_nroi_kern(const GModelSky& model,
  160 + const GCTAResponseIrf& rsp,
  161 + const GTime& srcTime,
  162 + const GEnergy& obsEng,
  163 + const GTime& obsTime,
  164 + const GObservation& obs) :
  165 + m_model(model),
  166 + m_rsp(rsp),
  167 + m_srcTime(srcTime),
  168 + m_obsEng(obsEng),
  169 + m_obsTime(obsTime),
  170 + m_obs(obs) {}
  171 + double eval(const double& logEsrc);
  172 +protected:
  173 + const GModelSky& m_model;
  174 + const GCTAResponseIrf& m_rsp; //!< CTA response function
  175 + const GTime& m_srcTime;
  176 + const GEnergy& m_obsEng;
  177 + const GTime& m_obsTime;
  178 + const GObservation& m_obs;
  179 +};
  180 +
  181 +
  182 +/***********************************************************************//**
147 183 * @class cta_irf_radial_kern_rho
148 184 *
149 185 * @brief Kernel for radial model zenith angle integration of IRF
... ... @@ -327,6 +363,8 @@ public:
327 363 const GModelSpatialRadial& model,
328 364 const GEnergy& srcEng,
329 365 const GTime& srcTime,
  366 + const GEnergy& obsEng,
  367 + const GTime& obsTime,
330 368 const GCTAObservation& obs,
331 369 const GMatrix& rot,
332 370 const double& dist,
... ... @@ -337,6 +375,8 @@ public:
337 375 m_model(model),
338 376 m_srcEng(srcEng),
339 377 m_srcTime(srcTime),
  378 + m_obsEng(obsEng),
  379 + m_obsTime(obsTime),
340 380 m_obs(obs),
341 381 m_rot(rot),
342 382 m_dist(dist),
... ... @@ -352,6 +392,8 @@ protected:
352 392 const GModelSpatialRadial& m_model; //!< Radial spatial model
353 393 const GEnergy& m_srcEng; //!< True photon energy
354 394 const GTime& m_srcTime; //!< True photon arrival time
  395 + const GEnergy& m_obsEng; //!< Observed photon energy
  396 + const GTime& m_obsTime; //!< Observed photon arrival time
355 397 const GCTAObservation& m_obs; //!< CTA observation
356 398 const GMatrix& m_rot; //!< Rotation matrix
357 399 const double& m_dist; //!< Distance model-ROI centre
... ... @@ -389,6 +431,8 @@ public:
389 431 cta_npred_radial_kern_omega(const GCTAResponseIrf& rsp,
390 432 const GEnergy& srcEng,
391 433 const GTime& srcTime,
  434 + const GEnergy& obsEng,
  435 + const GTime& obsTime,
392 436 const GCTAObservation& obs,
393 437 const GMatrix& rot,
394 438 double sin_rho,
... ... @@ -396,6 +440,8 @@ public:
396 440 m_rsp(rsp),
397 441 m_srcEng(srcEng),
398 442 m_srcTime(srcTime),
  443 + m_obsEng(obsEng),
  444 + m_obsTime(obsTime),
399 445 m_obs(obs),
400 446 m_rot(rot),
401 447 m_cos_rho(cos_rho),
... ... @@ -405,6 +451,8 @@ protected:
405 451 const GCTAResponseIrf& m_rsp; //!< CTA response
406 452 const GEnergy& m_srcEng; //!< True photon energy
407 453 const GTime& m_srcTime; //!< True photon arrival time
  454 + const GEnergy& m_obsEng; //!< Observed photon energy
  455 + const GTime& m_obsTime; //!< Observed photon arrival time
408 456 const GCTAObservation& m_obs; //!< CTA observation
409 457 const GMatrix& m_rot; //!< Rotation matrix
410 458 double m_cos_rho; //!< Cosine of offset angle
... ... @@ -620,6 +668,8 @@ public:
620 668 const double& posangle,
621 669 const GEnergy& srcEng,
622 670 const GTime& srcTime,
  671 + const GEnergy& obsEng,
  672 + const GTime& obsTime,
623 673 const GCTAObservation& obs,
624 674 const GMatrix& rot,
625 675 const double& rho_roi,
... ... @@ -633,6 +683,8 @@ public:
633 683 m_posangle(posangle),
634 684 m_srcEng(srcEng),
635 685 m_srcTime(srcTime),
  686 + m_obsEng(obsEng),
  687 + m_obsTime(obsTime),
636 688 m_obs(obs),
637 689 m_rot(rot),
638 690 m_rho_roi(rho_roi),
... ... @@ -651,6 +703,8 @@ protected:
651 703 const double& m_posangle; //!< Ellipse boundary position angle
652 704 const GEnergy& m_srcEng; //!< True photon energy
653 705 const GTime& m_srcTime; //!< True photon arrival time
  706 + const GEnergy& m_obsEng; //!< Observed photon energy
  707 + const GTime& m_obsTime; //!< Observed photon arrival time
654 708 const GCTAObservation& m_obs; //!< CTA observation
655 709 const GMatrix& m_rot; //!< Rotation matrix
656 710 const double& m_rho_roi; //!< Distance between model and ROI centre
... ... @@ -691,6 +745,8 @@ public:
691 745 const GModelSpatialElliptical& model,
692 746 const GEnergy& srcEng,
693 747 const GTime& srcTime,
  748 + const GEnergy& obsEng,
  749 + const GTime& obsTime,
694 750 const GCTAObservation& obs,
695 751 const GMatrix& rot,
696 752 const double& rho,
... ... @@ -701,6 +757,8 @@ public:
701 757 m_model(model),
702 758 m_srcEng(srcEng),
703 759 m_srcTime(srcTime),
  760 + m_obsEng(obsEng),
  761 + m_obsTime(obsTime),
704 762 m_obs(obs),
705 763 m_rot(rot),
706 764 m_rho(rho),
... ... @@ -713,6 +771,8 @@ protected:
713 771 const GModelSpatialElliptical& m_model; //!< Model
714 772 const GEnergy& m_srcEng; //!< True photon energy
715 773 const GTime& m_srcTime; //!< True photon arrival time
  774 + const GEnergy& m_obsEng; //!< Observed photon energy
  775 + const GTime& m_obsTime; //!< Observed photon arrival time
716 776 const GCTAObservation& m_obs; //!< Pointer to observation
717 777 const GMatrix& m_rot; //!< Rotation matrix
718 778 const double& m_rho;
... ... @@ -906,6 +966,8 @@ public:
906 966 const GModelSpatial& model,
907 967 const GEnergy& srcEng,
908 968 const GTime& srcTime,
  969 + const GEnergy& obsEng,
  970 + const GTime& obsTime,
909 971 const GCTAObservation& obs,
910 972 const GMatrix& rot,
911 973 const int& iter) :
... ... @@ -913,6 +975,8 @@ public:
913 975 m_model(model),
914 976 m_srcEng(srcEng),
915 977 m_srcTime(srcTime),
  978 + m_obsEng(obsEng),
  979 + m_obsTime(obsTime),
916 980 m_obs(obs),
917 981 m_rot(rot),
918 982 m_iter(iter) { }
... ... @@ -922,6 +986,8 @@ protected:
922 986 const GModelSpatial& m_model; //!< Spatial model
923 987 const GEnergy& m_srcEng; //!< True photon energy
924 988 const GTime& m_srcTime; //!< True photon arrival time
  989 + const GEnergy& m_obsEng; //!< Observed photon energy
  990 + const GTime& m_obsTime; //!< Observed photon arrival time
925 991 const GCTAObservation& m_obs; //!< CTA observation
926 992 const GMatrix& m_rot; //!< Rotation matrix
927 993 const int& m_iter; //!< Integration iterations
... ... @@ -954,6 +1020,8 @@ public:
954 1020 const GModelSpatial& model,
955 1021 const GEnergy& srcEng,
956 1022 const GTime& srcTime,
  1023 + const GEnergy& obsEng,
  1024 + const GTime& obsTime,
957 1025 const GCTAObservation& obs,
958 1026 const GMatrix& rot,
959 1027 const double& theta,
... ... @@ -962,6 +1030,8 @@ public:
962 1030 m_model(model),
963 1031 m_srcEng(srcEng),
964 1032 m_srcTime(srcTime),
  1033 + m_obsEng(obsEng),
  1034 + m_obsTime(obsTime),
965 1035 m_obs(obs),
966 1036 m_rot(rot),
967 1037 m_theta(theta),
... ... @@ -973,6 +1043,8 @@ protected:
973 1043 const GModelSpatial& m_model; //!< Spatial model
974 1044 const GEnergy& m_srcEng; //!< True photon energy
975 1045 const GTime& m_srcTime; //!< True photon arrival time
  1046 + const GEnergy& m_obsEng; //!< Observed photon energy
  1047 + const GTime& m_obsTime; //!< Observed photon arrival time
976 1048 const GCTAObservation& m_obs; //!< CTA observation
977 1049 const GMatrix& m_rot; //!< Rotation matrix
978 1050 const double& m_theta; //!< Offset angle (radians)
... ...
inst/lat/include/GLATResponse.hpp
  View file @ee42aa1
... ... @@ -67,19 +67,12 @@ public:
67 67 virtual double irf(const GEvent& event,
68 68 const GPhoton& photon,
69 69 const GObservation& obs) const;
70   - virtual double npred(const GPhoton& photon,
71   - const GObservation& obs) const;
  70 + virtual double nroi(const GModelSky& model,
  71 + const GEnergy& obsEng,
  72 + const GTime& obsTime,
  73 + const GObservation& obs) const;
72 74 virtual std::string print(const GChatter& chatter = NORMAL) const;
73 75  
74   - // New pure virtual methods
75   - virtual double convolve(const GModelSky& model,
76   - const GEvent& event,
77   - const GObservation& obs) const;
78   - virtual double nroi(const GModelSky& model,
79   - const GEnergy& obsEng,
80   - const GTime& obsTime,
81   - const GObservation& obs) const;
82   -
83 76 // Implemented virtual methods
84 77 virtual double irf(const GEvent& event,
85 78 const GSource& source,
... ...
inst/lat/pyext/GLATResponse.i
  View file @ee42aa1
... ... @@ -50,8 +50,10 @@ public:
50 50 virtual double irf(const GEvent& event,
51 51 const GPhoton& photon,
52 52 const GObservation& obs) const;
53   - virtual double npred(const GPhoton& photon,
54   - const GObservation& obs) const;
  53 + virtual double nroi(const GModelSky& model,
  54 + const GEnergy& obsEng,
  55 + const GTime& obsTime,
  56 + const GObservation& obs) const;
55 57  
56 58 // Implemented virtual methods
57 59 virtual double irf(const GEvent& event,
... ...
inst/lat/src/GLATResponse.cpp
  View file @ee42aa1
1 1 /***************************************************************************
2 2 * GLATResponse.cpp - Fermi/LAT response class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -49,6 +49,8 @@
49 49 #define G_LOAD "GLATResponse::load(std::string&)"
50 50 #define G_IRF "GLATResponse::irf(GInstDir&, GEnergy&, GTime&, GSkyDir&,"\
51 51 " GEnergy&, GTime&, GObservation&)"
  52 +#define G_NROI "GLATResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
  53 + "GObservation&)"
52 54 #define G_AEFF "GLATResponse::aeff(int&)"
53 55 #define G_PSF "GLATResponse::psf(int&)"
54 56 #define G_EDISP "GLATResponse::edisp(int&)"
... ... @@ -56,8 +58,6 @@
56 58 "GTime&, GObservation&)"
57 59 #define G_IRF_BIN "GLATResponse::irf(GLATEventBin&, GModel&, GEnergy&,"\
58 60 "GTime&, GObservation&)"
59   -#define G_NPRED "GLATResponse::npred(GSkyDir&, GEnergy&, GTime&,"\
60   - " GObservation&)"
61 61  
62 62 /* __ Macros _____________________________________________________________ */
63 63  
... ... @@ -487,24 +487,45 @@ double GLATResponse::irf(const GLATEventBin&amp; event,
487 487  
488 488  
489 489 /***********************************************************************//**
490   - * @brief Return integral of instrument response function.
  490 + * @brief Return integral of event probability for a given sky model over ROI
491 491 *
492   - * @param[in] photon Incident photon.
  492 + * @param[in] model Incident photon.
493 493 * @param[in] obs Observation.
  494 + * @return 0.0
494 495 *
495   - * @todo Not yet implemented.
  496 + * @exception GException::feature_not_implemented
  497 + * Method is not implemented.
  498 + *
  499 + * Computes the integral
  500 + *
  501 + * \f[
  502 + * N_{\rm ROI}(E',t') = \int_{\rm ROI} P(p',E',t') dp'
  503 + * \f]
  504 + *
  505 + * of the event probability
  506 + *
  507 + * \f[
  508 + * P(p',E',t') = \int \int \int
  509 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  510 + * \f]
  511 + *
  512 + * for a given sky model \f$S(p,E,t)\f$ and response function
  513 + * \f$R(p',E',t'|p,E,t)\f$ over the Region of Interest (ROI).
  514 + *
  515 + * @todo Implement method (is maybe not really needed)
496 516 ***************************************************************************/
497   -double GLATResponse::npred(const GPhoton& photon,
498   - const GObservation& obs) const
  517 +double GLATResponse::nroi(const GModelSky& model,
  518 + const GEnergy& obsEng,
  519 + const GTime& obsTime,
  520 + const GObservation& obs) const
499 521 {
500   - // Initialise
501   - double npred = 0.0;
502   -
503   - // Notify that method is not yet implemented
504   - throw GException::feature_not_implemented(G_NPRED);
  522 + // Method is not implemented
  523 + std::string msg = "Spatial integration of sky model over the data space "
  524 + "is not implemented.";
  525 + throw GException::feature_not_implemented(G_NROI, msg);
505 526  
506   - // Return integrated IRF value
507   - return npred;
  527 + // Return Npred
  528 + return (0.0);
508 529 }
509 530  
510 531  
... ...
inst/mwl/include/GMWLResponse.hpp
  View file @ee42aa1
... ... @@ -52,7 +52,7 @@ public:
52 52 virtual ~GMWLResponse(void);
53 53  
54 54 // Operators
55   - virtual GMWLResponse& operator= (const GMWLResponse & rsp);
  55 + virtual GMWLResponse& operator=(const GMWLResponse & rsp);
56 56  
57 57 // Implemented pure virtual methods
58 58 virtual void clear(void);
... ... @@ -63,19 +63,12 @@ public:
63 63 virtual double irf(const GEvent& event,
64 64 const GPhoton& photon,
65 65 const GObservation& obs) const;
66   - virtual double npred(const GPhoton& photon,
67   - const GObservation& obs) const;
  66 + virtual double nroi(const GModelSky& model,
  67 + const GEnergy& obsEng,
  68 + const GTime& obsTime,
  69 + const GObservation& obs) const;
68 70 virtual std::string print(const GChatter& chatter = NORMAL) const;
69 71  
70   - // New pure virtual methods
71   - virtual double convolve(const GModelSky& model,
72   - const GEvent& event,
73   - const GObservation& obs) const;
74   - virtual double nroi(const GModelSky& model,
75   - const GEnergy& obsEng,
76   - const GTime& obsTime,
77   - const GObservation& obs) const;
78   -
79 72 protected:
80 73 // Protected methods
81 74 void init_members(void);
... ... @@ -136,19 +129,4 @@ double GMWLResponse::irf(const GEvent&amp; event, const GPhoton&amp; photon,
136 129 return 1.0;
137 130 }
138 131  
139   -
140   -/***********************************************************************//**
141   - * @brief Return predicted number of events
142   - *
143   - * @param[in] photon Photon.
144   - * @param[in] obs Observation.
145   - *
146   - * @return Instrument response function (always 1).
147   - ***************************************************************************/
148   -inline
149   -double GMWLResponse::npred(const GPhoton& photon, const GObservation& obs) const
150   -{
151   - return 1.0;
152   -}
153   -
154 132 #endif /* GMWLRESPONSE_HPP */
... ...
inst/mwl/src/GMWLResponse.cpp
  View file @ee42aa1
1 1 /***************************************************************************
2 2 * GMWLResponse.cpp - Multi-wavelength response class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2015 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,6 +32,8 @@
32 32 #include "GMWLResponse.hpp"
33 33  
34 34 /* __ Method name definitions ____________________________________________ */
  35 +#define G_NROI "GMWLResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
  36 + "GObservation&)"
35 37  
36 38 /* __ Macros _____________________________________________________________ */
37 39  
... ... @@ -151,7 +153,7 @@ void GMWLResponse::clear(void)
151 153  
152 154 /***********************************************************************//**
153 155 * @brief Clone instance
154   -***************************************************************************/
  156 + ***************************************************************************/
155 157 GMWLResponse* GMWLResponse::clone(void) const
156 158 {
157 159 return new GMWLResponse(*this);
... ... @@ -159,6 +161,49 @@ GMWLResponse* GMWLResponse::clone(void) const
159 161  
160 162  
161 163 /***********************************************************************//**
  164 + * @brief Return integral of event probability for a given sky model over ROI
  165 + *
  166 + * @param[in] model Incident photon.
  167 + * @param[in] obs Observation.
  168 + * @return 0.0
  169 + *
  170 + * @exception GException::feature_not_implemented
  171 + * Method is not implemented.
  172 + *
  173 + * Computes the integral
  174 + *
  175 + * \f[
  176 + * N_{\rm ROI}(E',t') = \int_{\rm ROI} P(p',E',t') dp'
  177 + * \f]
  178 + *
  179 + * of the event probability
  180 + *
  181 + * \f[
  182 + * P(p',E',t') = \int \int \int
  183 + * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  184 + * \f]
  185 + *
  186 + * for a given sky model \f$S(p,E,t)\f$ and response function
  187 + * \f$R(p',E',t'|p,E,t)\f$ over the Region of Interest (ROI).
  188 + *
  189 + * @todo Implement method (is maybe not really needed)
  190 + ***************************************************************************/
  191 +double GMWLResponse::nroi(const GModelSky& model,
  192 + const GEnergy& obsEng,
  193 + const GTime& obsTime,
  194 + const GObservation& obs) const
  195 +{
  196 + // Method is not implemented
  197 + std::string msg = "Spatial integration of sky model over the data space "
  198 + "is not implemented.";
  199 + throw GException::feature_not_implemented(G_NROI, msg);
  200 +
  201 + // Return Npred
  202 + return (0.0);
  203 +}
  204 +
  205 +
  206 +/***********************************************************************//**
162 207 * @brief Print response information
163 208 *
164 209 * @param[in] chatter Chattiness (defaults to NORMAL).
... ...
pyext/GResponse.i
  View file @ee42aa1
... ... @@ -51,8 +51,14 @@ public:
51 51 virtual double irf(const GEvent& event,
52 52 const GPhoton& photon,
53 53 const GObservation& obs) const = 0;
54   - virtual double npred(const GPhoton& photon,
55   - const GObservation& obs) const = 0;
  54 + virtual double convolve(const GModelSky& model,
  55 + const GEvent& event,
  56 + const GObservation& obs,
  57 + const bool& grad = true) const;
  58 + virtual double nroi(const GModelSky& model,
  59 + const GEnergy& obsEng,
  60 + const GTime& obsTime,
  61 + const GObservation& obs) const = 0;
56 62  
57 63 // Virtual methods
58 64 virtual double irf(const GEvent& event,
... ... @@ -70,16 +76,6 @@ public:
70 76 virtual double irf_diffuse(const GEvent& event,
71 77 const GSource& source,
72 78 const GObservation& obs) const;
73   - virtual double npred(const GSource& source,
74   - const GObservation& obs) const;
75   - virtual double npred_ptsrc(const GSource& source,
76   - const GObservation& obs) const;
77   - virtual double npred_radial(const GSource& source,
78   - const GObservation& obs) const;
79   - virtual double npred_elliptical(const GSource& source,
80   - const GObservation& obs) const;
81   - virtual double npred_diffuse(const GSource& source,
82   - const GObservation& obs) const;
83 79 virtual GEbounds ebounds_src(const GEnergy& obsEnergy) const;
84 80 };
85 81  
... ...
src/model/GModelSky.cpp
  View file @ee42aa1
... ... @@ -562,9 +562,10 @@ double GModelSky::npred(const GEnergy&amp; obsEng, const GTime&amp; obsTime,
562 562 // Continue only if model is valid)
563 563 if (valid_model()) {
564 564  
565   - // Get response function
566   - const GResponse* rsp = obs.response();
  565 + // Compute Nroi
  566 + npred = obs.response()->nroi(*this, obsEng, obsTime, obs);
567 567  
  568 + /*
568 569 // Here we make the simplifying approximations
569 570 // srcEng=obsEng and srcTime=obsTime. To be fully correct we should
570 571 // integrate over true energy and true time here ... at least true
... ... @@ -587,6 +588,7 @@ double GModelSky::npred(const GEnergy&amp; obsEng, const GTime&amp; obsTime,
587 588 if (!m_scales.empty()) {
588 589 npred *= scale(obs.instrument()).value();
589 590 }
  591 + */
590 592  
591 593 // Compile option: Check for NaN/Inf
592 594 #if defined(G_NAN_CHECK)
... ... @@ -594,11 +596,8 @@ double GModelSky::npred(const GEnergy&amp; obsEng, const GTime&amp; obsTime,
594 596 std::cout << "*** ERROR: GModelSky::npred:";
595 597 std::cout << " NaN/Inf encountered";
596 598 std::cout << " (npred=" << npred;
597   - std::cout << ", npred_spatial=" << npred_spatial;
598   - std::cout << ", npred_spectral=" << npred_spectral;
599   - std::cout << ", npred_temporal=" << npred_temporal;
600   - std::cout << ", srcEng=" << srcEng;
601   - std::cout << ", srcTime=" << srcTime;
  599 + std::cout << ", obsEng=" << obsEng;
  600 + std::cout << ", obsTime=" << obsTime;
602 601 std::cout << ")" << std::endl;
603 602 }
604 603 #endif
... ...
src/obs/GResponse.cpp
  View file @ee42aa1
... ... @@ -38,6 +38,7 @@
38 38 #include "GException.hpp"
39 39 #include "GTools.hpp"
40 40 #include "GMath.hpp"
  41 +#include "GModelSky.hpp"
41 42 #include "GModelSpatialPointSource.hpp"
42 43 #include "GModelSpatialRadial.hpp"
43 44 #include "GModelSpatialElliptical.hpp"
... ... @@ -153,19 +154,13 @@ GResponse&amp; GResponse::operator=(const GResponse&amp; rsp)
153 154 = =
154 155 ==========================================================================*/
155 156  
156   -
157   -/*==========================================================================
158   - = =
159   - = Obsolete public methods =
160   - = =
161   - ==========================================================================*/
162   -
163 157 /***********************************************************************//**
164 158 * @brief Convolve sky model with the instrument response
165 159 *
166 160 * @param[in] model Sky model.
167 161 * @param[in] event Event.
168 162 * @param[in] obs Observation.
  163 + * @param[in] grad Should model gradients be computed? (default: true)
169 164 * @return Event probability.
170 165 *
171 166 * Computes the event probability
... ... @@ -181,48 +176,47 @@ GResponse&amp; GResponse::operator=(const GResponse&amp; rsp)
181 176 ***************************************************************************/
182 177 double GResponse::convolve(const GModelSky& model,
183 178 const GEvent& event,
184   - const GObservation& obs) const
  179 + const GObservation& obs,
  180 + const bool& grad) const
185 181 {
186 182 // Initialise result
187 183 double prob = 0.0;
188 184  
189   - // Get response function
190   - const GResponse* rsp = obs.response();
191   -
192 185 // Continue only if the model has a spatial component
193   - if (spatial() != NULL) {
  186 + if (model.spatial() != NULL) {
194 187  
195   - // Get response function
196   - const GResponse* rsp = obs.response();
  188 + // Get source energy and time (no dispersion)
  189 + GEnergy srcEng = event.energy();
  190 + GTime srcTime = event.time();
197 191  
198 192 // Set source
199   - GSource source(this->name(), m_spatial, srcEng, srcTime);
  193 + GSource source(model.name(), model.spatial(), srcEng, srcTime);
200 194  
201 195 // Get IRF value. This method returns the spatial component of the
202 196 // source model.
203   - double irf = rsp->irf(event, source, obs);
  197 + double irf = this->irf(event, source, obs);
204 198  
205 199 // If required, apply instrument specific model scaling
206   - if (!m_scales.empty()) {
207   - irf *= scale(obs.instrument()).value();
  200 + if (model.has_scales()) {
  201 + irf *= model.scale(obs.instrument()).value();
208 202 }
209 203  
210 204 // Case A: evaluate gradients
211 205 if (grad) {
212 206  
213 207 // Evaluate source model
214   - double spec = (spectral() != NULL) ? spectral()->eval_gradients(srcEng, srcTime) : 1.0;
215   - double temp = (temporal() != NULL) ? temporal()->eval_gradients(srcTime) : 1.0;
  208 + double spec = (model.spectral() != NULL) ? model.spectral()->eval_gradients(srcEng, srcTime) : 1.0;
  209 + double temp = (model.temporal() != NULL) ? model.temporal()->eval_gradients(srcTime) : 1.0;
216 210  
217   - // Set value
218   - value = spec * temp * irf;
  211 + // Set probability
  212 + prob = spec * temp * irf;
219 213  
220 214 // Compile option: Check for NaN/Inf
221 215 #if defined(G_NAN_CHECK)
222   - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
223   - std::cout << "*** ERROR: GModelSky::integrate_dir:";
  216 + if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) {
  217 + std::cout << "*** ERROR: GResponse::convolve:";
224 218 std::cout << " NaN/Inf encountered";
225   - std::cout << " (value=" << value;
  219 + std::cout << " (prob=" << prob;
226 220 std::cout << ", spec=" << spec;
227 221 std::cout << ", temp=" << temp;
228 222 std::cout << ", irf=" << irf;
... ... @@ -231,21 +225,21 @@ double GResponse::convolve(const GModelSky&amp; model,
231 225 #endif
232 226  
233 227 // Multiply factors to spectral gradients
234   - if (spectral() != NULL) {
  228 + if (model.spectral() != NULL) {
235 229 double fact = temp * irf;
236 230 if (fact != 1.0) {
237   - for (int i = 0; i < spectral()->size(); ++i) {
238   - (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
  231 + for (int i = 0; i < model.spectral()->size(); ++i) {
  232 + (*model.spectral())[i].factor_gradient((*model.spectral())[i].factor_gradient() * fact);
239 233 }
240 234 }
241 235 }
242 236  
243 237 // Multiply factors to temporal gradients
244   - if (temporal() != NULL) {
  238 + if (model.temporal() != NULL) {
245 239 double fact = spec * irf;
246 240 if (fact != 1.0) {
247   - for (int i = 0; i < temporal()->size(); ++i) {
248   - (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
  241 + for (int i = 0; i < model.temporal()->size(); ++i) {
  242 + (*model.temporal())[i].factor_gradient((*model.temporal())[i].factor_gradient() * fact);
249 243 }
250 244 }
251 245 }
... ... @@ -256,18 +250,18 @@ double GResponse::convolve(const GModelSky&amp; model,
256 250 else {
257 251  
258 252 // Evaluate source model
259   - double spec = (m_spectral != NULL) ? m_spectral->eval(srcEng, srcTime) : 1.0;
260   - double temp = (m_temporal != NULL) ? m_temporal->eval(srcTime) : 1.0;
  253 + double spec = (model.spectral() != NULL) ? model.spectral()->eval(srcEng, srcTime) : 1.0;
  254 + double temp = (model.temporal() != NULL) ? model.temporal()->eval(srcTime) : 1.0;
261 255  
262   - // Set value
263   - value = spec * temp * irf;
  256 + // Set probability
  257 + prob = spec * temp * irf;
264 258  
265 259 // Compile option: Check for NaN/Inf
266 260 #if defined(G_NAN_CHECK)
267   - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
268   - std::cout << "*** ERROR: GModelSky::integrate_dir:";
  261 + if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) {
  262 + std::cout << "*** ERROR: GResponse::convolve:";
269 263 std::cout << " NaN/Inf encountered";
270   - std::cout << " (value=" << value;
  264 + std::cout << " (prob=" << prob;
271 265 std::cout << ", spec=" << spec;
272 266 std::cout << ", temp=" << temp;
273 267 std::cout << ", irf=" << irf;
... ... @@ -279,22 +273,17 @@ double GResponse::convolve(const GModelSky&amp; model,
279 273  
280 274 } // endif: Gamma-ray source model had a spatial component
281 275  
282   - // Compile option: Check for NaN/Inf
283   - #if defined(G_NAN_CHECK)
284   - if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) {
285   - std::cout << "*** ERROR: GResponse::convolve:";
286   - std::cout << " NaN/Inf encountered";
287   - std::cout << " (value=" << prob;
288   - std::cout << ", event=" << event;
289   - std::cout << ")" << std::endl;
290   - }
291   - #endif
292   -
293 276 // Return event probability
294 277 return prob;
295 278 }
296 279  
297 280  
  281 +/*==========================================================================
  282 + = =
  283 + = Obsolete public methods =
  284 + = =
  285 + ==========================================================================*/
  286 +
298 287 /***********************************************************************//**
299 288 * @brief Return value of instrument response function
300 289 *
... ... @@ -448,6 +437,7 @@ double GResponse::irf_diffuse(const GEvent&amp; event,
448 437 * The method applies the deadtime correction, so that the result can be
449 438 * directly multiplied by the exposure time (also known as ontime).
450 439 ***************************************************************************/
  440 +/*
451 441 double GResponse::npred(const GSource& source, const GObservation& obs) const
452 442 {
453 443 // Initialise Npred value
... ... @@ -474,7 +464,7 @@ double GResponse::npred(const GSource&amp; source, const GObservation&amp; obs) const
474 464 // Return response value
475 465 return npred;
476 466 }
477   -
  467 +*/
478 468  
479 469 /***********************************************************************//**
480 470 * @brief Return ROI integral of point source model
... ... @@ -487,6 +477,7 @@ double GResponse::npred(const GSource&amp; source, const GObservation&amp; obs) const
487 477 * region of interest. The method assumes that source.model() is of type
488 478 * GModelSpatialPointSource.
489 479 ***************************************************************************/
  480 +/*
490 481 double GResponse::npred_ptsrc(const GSource& source,
491 482 const GObservation& obs) const
492 483 {
... ... @@ -503,7 +494,7 @@ double GResponse::npred_ptsrc(const GSource&amp; source,
503 494 // Return Npred
504 495 return npred;
505 496 }
506   -
  497 +*/
507 498  
508 499 /***********************************************************************//**
509 500 * @brief Return spatial integral of radial source model
... ... @@ -516,6 +507,7 @@ double GResponse::npred_ptsrc(const GSource&amp; source,
516 507 * the region of interest. The method assumes that source.model() is of type
517 508 * GModelSpatialRadial.
518 509 ***************************************************************************/
  510 +/*
519 511 double GResponse::npred_radial(const GSource& source,
520 512 const GObservation& obs) const
521 513 {
... ... @@ -578,7 +570,7 @@ double GResponse::npred_radial(const GSource&amp; source,
578 570 // Return Npred
579 571 return npred;
580 572 }
581   -
  573 +*/
582 574  
583 575 /***********************************************************************//**
584 576 * @brief Return spatial integral of elliptical source model
... ... @@ -590,6 +582,7 @@ double GResponse::npred_radial(const GSource&amp; source,
590 582 * the region of interest. The method assumes that source.model() is of type
591 583 * GModelSpatialElliptical.
592 584 ***************************************************************************/
  585 +/*
593 586 double GResponse::npred_elliptical(const GSource& source,
594 587 const GObservation& obs) const
595 588 {
... ... @@ -652,7 +645,7 @@ double GResponse::npred_elliptical(const GSource&amp; source,
652 645 // Return Npred
653 646 return npred;
654 647 }
655   -
  648 +*/
656 649  
657 650 /***********************************************************************//**
658 651 * @brief Return spatial integral of diffuse source model
... ... @@ -663,6 +656,7 @@ double GResponse::npred_elliptical(const GSource&amp; source,
663 656 * @exception GException::feature_not_implemented
664 657 * Method not yet implemented.
665 658 ***************************************************************************/
  659 +/*
666 660 double GResponse::npred_diffuse(const GSource& source,
667 661 const GObservation& obs) const
668 662 {
... ... @@ -673,270 +667,7 @@ double GResponse::npred_diffuse(const GSource&amp; source,
673 667 // Return Npred
674 668 return 0.0;
675 669 }
676   -
677   -
678   -/***********************************************************************//**
679   - * @brief Return data space integral of instrument response function
680   - *
681   - * @param[in] source Source.
682   - * @param[in] obsEng Observed energy.
683   - * @param[in] obsTime Observed time.
684   - * @param[in] obs Observation.
685   - * @return Integral of model over ROI.
686   - *
687   - * Returns the data space integral of the instrument response function for
688   - * a given source and a particular observation. This method is needed for
689   - * an unbinned maximum likelihood analysis.
690   - *
691   - * The method applies the deadtime correction, so that the result can be
692   - * directly multiplied by the exposure time (also known as ontime).
693   - ***************************************************************************/
694   -double GResponse::npred(const GSource& source,
695   - const GEnergy& obsEng,
696   - const GTime& obsTime,
697   - const GObservation& obs) const
698   -{
699   - // Initialise Npred value
700   - double npred = 0.0;
701   -
702   - // Select NPRED depending on the spatial model type
703   - switch (source.model()->code()) {
704   - case GMODEL_SPATIAL_POINT_SOURCE:
705   - npred = npred_ptsrc(source, obsEng, obsTime, obs);
706   - break;
707   - case GMODEL_SPATIAL_RADIAL:
708   - npred = npred_radial(source, obsEng, obsTime, obs);
709   - break;
710   - case GMODEL_SPATIAL_ELLIPTICAL:
711   - npred = npred_elliptical(source, obsEng, obsTime, obs);
712   - break;
713   - case GMODEL_SPATIAL_DIFFUSE:
714   - npred = npred_diffuse(source, obsEng, obsTime, obs);
715   - break;
716   - default:
717   - break;
718   - }
719   -
720   - // Return response value
721   - return npred;
722   -}
723   -
724   -
725   -/***********************************************************************//**
726   - * @brief Return ROI integral of point source model
727   - *
728   - * @param[in] source Source.
729   - * @param[in] obsEng Observed energy.
730   - * @param[in] obsTime Observed time.
731   - * @param[in] obs Observation.
732   - * @return Integral of point source model over ROI.
733   - *
734   - * This method returns the spatial integral of a point source model over the
735   - * region of interest. The method assumes that source.model() is of type
736   - * GModelSpatialPointSource.
737   - ***************************************************************************/
738   -double GResponse::npred_ptsrc(const GSource& source,
739   - const GEnergy& obsEng,
740   - const GTime& obsTime,
741   - const GObservation& obs) const
742   -{
743   - // Get point source spatial model
744   - const GModelSpatialPointSource* src =
745   - static_cast<const GModelSpatialPointSource*>(source.model());
746   -
747   - // Set Photon
748   - GPhoton photon(src->dir(), source.energy(), source.time());
749   -
750   - // Compute Npred
751   - double npred = this->npred(photon, obsEng, obsTime, obs);
752   -
753   - // Return Npred
754   - return npred;
755   -}
756   -
757   -
758   -/***********************************************************************//**
759   - * @brief Return spatial integral of radial source model
760   - *
761   - * @param[in] source Source.
762   - * @param[in] obsEng Observed energy.
763   - * @param[in] obsTime Observed time.
764   - * @param[in] obs Observation.
765   - * @return Integral of radial source model over ROI.
766   - *
767   - * This method returns the spatial integral of a radial source model over
768   - * the region of interest. The method assumes that source.model() is of type
769   - * GModelSpatialRadial.
770   - ***************************************************************************/
771   -double GResponse::npred_radial(const GSource& source,
772   - const GEnergy& obsEng,
773   - const GTime& obsTime,
774   - const GObservation& obs) const
775   -{
776   - // Initialise Npred value
777   - double npred = 0.0;
778   -
779   - // Get radial spatial model
780   - const GModelSpatialRadial* spatial =
781   - static_cast<const GModelSpatialRadial*>(source.model());
782   -
783   - // Compute rotation matrix to convert from native coordinates given
784   - // by (theta,phi) into celestial coordinates.
785   - GMatrix ry;
786   - GMatrix rz;
787   - ry.eulery(spatial->dec() - 90.0);
788   - rz.eulerz(-spatial->ra());
789   - GMatrix rot = (ry * rz).transpose();
790   -
791   - // Set offset angle integration range
792   - double theta_min = 0.0;
793   - double theta_max = spatial->theta_max();
794   -
795   - // Perform offset angle integration if interval is valid
796   - if (theta_max > theta_min) {
797   -
798   - // Setup integration kernel
799   - GResponse::npred_radial_kern_theta integrand(*this,
800   - *spatial,
801   - source.energy(),
802   - source.time(),
803   - obs,
804   - rot);
805   -
806   - // Integrate over theta
807   - GIntegral integral(&integrand);
808   - npred = integral.romberg(theta_min, theta_max);
809   -
810   - // Compile option: Show integration results
811   - #if defined(G_DEBUG_NPRED_RADIAL)
812   - std::cout << "GResponse::npred_radial:";
813   - std::cout << " theta_min=" << theta_min;
814   - std::cout << " theta_max=" << theta_max;
815   - std::cout << " npred=" << npred << std::endl;
816   - #endif
817   -
818   - } // endif: offset angle range was valid
819   -
820   - // Debug: Check for NaN
821   - #if defined(G_NAN_CHECK)
822   - if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
823   - std::cout << "*** ERROR: GResponse::npred_radial:";
824   - std::cout << " NaN/Inf encountered";
825   - std::cout << " (npred=" << npred;
826   - std::cout << ", theta_min=" << theta_min;
827   - std::cout << ", theta_max=" << theta_max;
828   - std::cout << ")" << std::endl;
829   - }
830   - #endif
831   -
832   - // Return Npred
833   - return npred;
834   -}
835   -
836   -
837   -/***********************************************************************//**
838   - * @brief Return spatial integral of elliptical source model
839   - *
840   - * @param[in] source Source.
841   - * @param[in] obsEng Observed energy.
842   - * @param[in] obsTime Observed time.
843   - * @param[in] obs Observation.
844   - * @return Integral of elliptical source model over ROI.
845   - *
846   - * This method returns the spatial integral of a radial source model over
847   - * the region of interest. The method assumes that source.model() is of type
848   - * GModelSpatialElliptical.
849   - ***************************************************************************/
850   -double GResponse::npred_elliptical(const GSource& source,
851   - const GEnergy& obsEng,
852   - const GTime& obsTime,
853   - const GObservation& obs) const
854   -{
855   - // Initialise Npred value
856   - double npred = 0.0;
857   -
858   - // Get elliptical spatial model
859   - const GModelSpatialElliptical* spatial =
860   - static_cast<const GModelSpatialElliptical*>(source.model());
861   -
862   - // Compute rotation matrix to convert from native coordinates given
863   - // by (theta,phi) into celestial coordinates.
864   - GMatrix ry;
865   - GMatrix rz;
866   - ry.eulery(spatial->dec() - 90.0);
867   - rz.eulerz(-spatial->ra());
868   - GMatrix rot = (ry * rz).transpose();
869   -
870   - // Set offset angle integration range
871   - double theta_min = 0.0;
872   - double theta_max = spatial->theta_max();
873   -
874   - // Perform offset angle integration if interval is valid
875   - if (theta_max > theta_min) {
876   -
877   - // Setup integration kernel
878   - GResponse::npred_elliptical_kern_theta integrand(*this,
879   - *spatial,
880   - source.energy(),
881   - source.time(),
882   - obs,
883   - rot);
884   -
885   - // Integrate over theta
886   - GIntegral integral(&integrand);
887   - npred = integral.romberg(theta_min, theta_max);
888   -
889   - // Compile option: Show integration results
890   - #if defined(G_DEBUG_NPRED_ELLIPTICAL)
891   - std::cout << "GResponse::npred_elliptical:";
892   - std::cout << " theta_min=" << theta_min;
893   - std::cout << " theta_max=" << theta_max;
894   - std::cout << " npred=" << npred << std::endl;
895   - #endif
896   -
897   - } // endif: performed offset angle integration
898   -
899   - // Debug: Check for NaN
900   - #if defined(G_NAN_CHECK)
901   - if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
902   - std::cout << "*** ERROR: GResponse::npred_elliptical:";
903   - std::cout << " NaN/Inf encountered";
904   - std::cout << " (npred=" << npred;
905   - std::cout << ", theta_min=" << theta_min;
906   - std::cout << ", theta_max=" << theta_max;
907   - std::cout << ")" << std::endl;
908   - }
909   - #endif
910   -
911   - // Return Npred
912   - return npred;
913   -}
914   -
915   -
916   -/***********************************************************************//**
917   - * @brief Return spatial integral of diffuse source model
918   - *
919   - * @param[in] source Source.
920   - * @param[in] obsEng Observed energy.
921   - * @param[in] obsTime Observed time.
922   - * @param[in] obs Observation.
923   - * @return Integral of diffuse source model over ROI.
924   - *
925   - * @exception GException::feature_not_implemented
926   - * Method not yet implemented.
927   - ***************************************************************************/
928   -double GResponse::npred_diffuse(const GSource& source,
929   - const GEnergy& obsEng,
930   - const GTime& obsTime,
931   - const GObservation& obs) const
932   -{
933   - // Feature not yet implemented
934   - throw GException::feature_not_implemented(G_NPRED_DIFFUSE,
935   - "Npred computation not implemented for diffuse models.");
936   -
937   - // Return Npred
938   - return 0.0;
939   -}
  670 +*/
940 671  
941 672  
942 673 /***********************************************************************//**
... ... @@ -1006,6 +737,7 @@ void GResponse::free_members(void)
1006 737 *
1007 738 * @param[in] theta Radial model offset angle (radians).
1008 739 ***************************************************************************/
  740 +/*
1009 741 double GResponse::npred_radial_kern_theta::eval(const double& theta)
1010 742 {
1011 743 // Initialise Npred value
... ... @@ -1046,6 +778,7 @@ double GResponse::npred_radial_kern_theta::eval(const double&amp; theta)
1046 778 // Return Npred
1047 779 return npred;
1048 780 }
  781 +*/
1049 782  
1050 783  
1051 784 /***********************************************************************//**
... ... @@ -1053,6 +786,7 @@ double GResponse::npred_radial_kern_theta::eval(const double&amp; theta)
1053 786 *
1054 787 * @param[in] phi Azimuth angle (radians).
1055 788 ***************************************************************************/
  789 +/*
1056 790 double GResponse::npred_radial_kern_phi::eval(const double& phi)
1057 791 {
1058 792 // Compute sky direction vector in native coordinates
... ... @@ -1089,6 +823,7 @@ double GResponse::npred_radial_kern_phi::eval(const double&amp; phi)
1089 823 // Return Npred
1090 824 return npred;
1091 825 }
  826 +*/
1092 827  
1093 828  
1094 829 /***********************************************************************//**
... ... @@ -1096,6 +831,7 @@ double GResponse::npred_radial_kern_phi::eval(const double&amp; phi)
1096 831 *
1097 832 * @param[in] theta Radial model offset angle (radians).
1098 833 ***************************************************************************/
  834 +/*
1099 835 double GResponse::npred_elliptical_kern_theta::eval(const double& theta)
1100 836 {
1101 837 // Initialise Npred value
... ... @@ -1133,6 +869,7 @@ double GResponse::npred_elliptical_kern_theta::eval(const double&amp; theta)
1133 869 // Return Npred
1134 870 return npred;
1135 871 }
  872 +*/
1136 873  
1137 874  
1138 875 /***********************************************************************//**
... ... @@ -1140,6 +877,7 @@ double GResponse::npred_elliptical_kern_theta::eval(const double&amp; theta)
1140 877 *
1141 878 * @param[in] phi Azimuth angle (radians).
1142 879 ***************************************************************************/
  880 +/*
1143 881 double GResponse::npred_elliptical_kern_phi::eval(const double& phi)
1144 882 {
1145 883 // Compute sky direction vector in native coordinates
... ... @@ -1180,3 +918,4 @@ double GResponse::npred_elliptical_kern_phi::eval(const double&amp; phi)
1180 918 // Return Npred
1181 919 return npred;
1182 920 }
  921 +*/
... ...