Finalize response optimization.

Jürgen Knödlseder
1 parent ab4461b6
Showing 7 changed files with 440 additions and 143 deletions
inst/cta/src/GCTAResponseCube.cpp
inst/cta/src/GCTAResponseIrf.cpp
inst/cta/src/GCTAResponse_helpers.cpp
inst/cta/src/GCTAResponse_helpers.hpp
inst/cta/test/test_irf_radial.py
inst/cta/test/test_likelihood_profile.py
inst/cta/test/test_npred_elliptical.py
@@ -64,6 +64,7 @@
 /* __ Macros _____________________________________________________________ */
 /* __ Coding definitions _________________________________________________ */
+#define G_RADIAL_PSF_BASED  //!< Use Psf-based integration for radial model
 /* __ Debug definitions __________________________________________________ */
@@ -237,11 +238,12 @@ GCTAResponseCube* GCTAResponseCube::clone(void) const
 /***********************************************************************//**
- * @brief Return value of instrument response function
+ * @brief Return instrument response
  *
  * @param[in] event Observed event.
  * @param[in] photon Incident photon.
  * @param[in] obs Observation.
+ * @return Instrument response.
  ***************************************************************************/
 double GCTAResponseCube::irf(const GEvent&       event,
                              const GPhoton&      photon,
@@ -308,14 +310,14 @@ double GCTAResponseCube::irf(const GEvent&amp;       event,
 /***********************************************************************//**
- * @brief Return value of instrument response function
+ * @brief Return instrument response
  *
  * @param[in] event Event.
  * @param[in] source Source.
  * @param[in] obs Observation.
+ * @return Instrument response.
  *
- * Returns the instrument response function for a given event, source and
- * observation.
+ * Returns the instrument response for a given event, source and observation.
  ***************************************************************************/
 double GCTAResponseCube::irf(const GEvent&       event,
                              const GSource&      source,
@@ -348,15 +350,14 @@ double GCTAResponseCube::irf(const GEvent&amp;       event,
 /***********************************************************************//**
- * @brief Return value of point source instrument response function
+ * @brief Return instrument response to point source
  *
  * @param[in] event Observed event.
  * @param[in] source Source.
  * @param[in] obs Observation.
- * @return Value of instrument response function for a point source.
+ * @return Instrument response to point source.
  *
- * This method returns the value of the instrument response function for a
- * point source.
+ * Returns the instrument response to a specified point source.
  ***************************************************************************/
 double GCTAResponseCube::irf_ptsrc(const GEvent&       event,
                                    const GSource&      source,
@@ -426,15 +427,14 @@ double GCTAResponseCube::irf_ptsrc(const GEvent&amp;       event,
 /***********************************************************************//**
- * @brief Return value of radial source instrument response function
+ * @brief Return instrument response to radial source
  *
  * @param[in] event Observed event.
  * @param[in] source Source.
  * @param[in] obs Observation.
- * @return Value of instrument response function for a radial source.
+ * @return Instrument response to radial source.
  *
- * This method returns the value of the instrument response function for a
- * radial source. Response values are computed on the fly.
+ * Returns the instrument response to a specified radial source.
  ***************************************************************************/
 double GCTAResponseCube::irf_radial(const GEvent&       event,
                                     const GSource&      source,
@@ -497,15 +497,14 @@ double GCTAResponseCube::irf_radial(const GEvent&amp;       event,
 /***********************************************************************//**
- * @brief Return value of elliptical source instrument response function
+ * @brief Return instrument response to elliptical source
  *
  * @param[in] event Observed event.
  * @param[in] source Source.
  * @param[in] obs Observation.
- * @return Value of instrument response function for an elliptical source.
+ * @return Instrument response to elliptical source.
  *
- * This method returns the value of the instrument response function for an
- * elliptical source. Response values are computed on the fly.
+ * Returns the instrument response to a specified elliptical source.
  ***************************************************************************/
 double GCTAResponseCube::irf_elliptical(const GEvent&       event,
                                         const GSource&      source,
@@ -569,21 +568,20 @@ double GCTAResponseCube::irf_elliptical(const GEvent&amp;       event,
 /***********************************************************************//**
- * @brief Return value of extended source instrument response function
+ * @brief Return instrument response to diffuse source
  *
  * @param[in] event Observed event.
  * @param[in] source Source.
  * @param[in] obs Observation.
- * @return Value of instrument response function for an extended source.
+ * @return Instrument response to diffuse source.
  *
- * This method returns the value of the instrument response function for an
- * extended source. It uses a pre-computation cache to store the IRF for
- * the spatial model component.
+ * Returns the instrument response to a specified diffuse source.
  *
- * The pre-computation cache is initialised if no cache has yet been
- * allocated, or if at the beginning of a scan over the events, the model
- * parameters have changed. The beginning of a scan is defined by an event
- * bin index of 0.
+ * The method uses a pre-computation cache to store the instrument response
+ * for the spatial model component. The pre-computation cache is initialised
+ * if no cache has yet been allocated, or if at the beginning of a scan over
+ * the events, the model parameters have changed. The beginning of a scan is
+ * defined by an event bin index of 0.
  ***************************************************************************/
 double GCTAResponseCube::irf_diffuse(const GEvent&       event,
                                      const GSource&      source,
@@ -652,9 +650,13 @@ double GCTAResponseCube::irf_diffuse(const GEvent&amp;       event,
  *
  * @param[in] photon Incident photon.
  * @param[in] obs Observation.
+ * @return Spatial integral of point spread function. 
  *
  * @exception GException::feature_not_implemented
  *            Method not implemented.
+ *
+ * This method is a dummy method that is required for a class derived from
+ * GResponse but that needs not to be implemented for a binned response.
  ***************************************************************************/
 double GCTAResponseCube::npred(const GPhoton&      photon,
                                const GObservation& obs) const
@@ -669,21 +671,21 @@ double GCTAResponseCube::npred(const GPhoton&amp;      photon,
 /***********************************************************************//**
- * @brief Read response
+ * @brief Read response information from Xml element
  *
  * @param[in] xml XML element.
  *
  * @exception GException::xml_invalid_parnames
  *            Invalid parameter names found in XML element.
  *
- * Reads information for a CTA observation from an XML element. The exposure
- * and PSF cubes are specified using
+ * Reads response information from an Xml element. The Exposure and Psf cubes
+ * are specified using
  *
- *     <observation name="..." id="..." instrument="...">
- *       ...
- *       <parameter name="ExposureCube" file="..."/>
- *       <parameter name="PsfCube"      file="..."/>
- *     </observation>
+ *      <observation name="..." id="..." instrument="...">
+ *        ...
+ *        <parameter name="ExposureCube" file="..."/>
+ *        <parameter name="PsfCube"      file="..."/>
+ *      </observation>
  *
  ***************************************************************************/
 void GCTAResponseCube::read(const GXmlElement& xml)
@@ -738,18 +740,18 @@ void GCTAResponseCube::read(const GXmlElement&amp; xml)
 /***********************************************************************//**
- * @brief Write response
+ * @brief Write response information into Xml element
  *
  * @param[in] xml XML element.
  *
- * Writes information for a CTA observation into an XML element. The exposure
- * and PSF cubes are specified using
+ * Writes response information into an Xml element. The Exposure and Psf
+ * cubes are specified using
  *
- *     <observation name="..." id="..." instrument="...">
- *       ...
- *       <parameter name="ExposureCube" file="..."/>
- *       <parameter name="PsfCube"      file="..."/>
- *     </observation>
+ *      <observation name="..." id="..." instrument="...">
+ *        ...
+ *        <parameter name="ExposureCube" file="..."/>
+ *        <parameter name="PsfCube"      file="..."/>
+ *      </observation>
  *
  ***************************************************************************/
 void GCTAResponseCube::write(GXmlElement& xml) const
@@ -774,10 +776,10 @@ void GCTAResponseCube::write(GXmlElement&amp; xml) const
 /***********************************************************************//**
- * @brief Print CTA response information
+ * @brief Print response information
  *
  * @param[in] chatter Chattiness (defaults to NORMAL).
- * @return String containing CTA response information.
+ * @return String containing response information.
  ***************************************************************************/
 std::string GCTAResponseCube::print(const GChatter& chatter) const
 {
@@ -819,20 +821,6 @@ std::string GCTAResponseCube::print(const GChatter&amp; chatter) const
 /*==========================================================================
  =                                                                         =
- =              Model type dependent CTA response methods                  =
- =                                                                         =
- ==========================================================================*/
-
-
-/*==========================================================================
- =                                                                         =
- =                    Low-level CTA response methods                       =
- =                                                                         =
- ==========================================================================*/
-
-
-/*==========================================================================
- =                                                                         =
  =                             Private methods                             =
  =                                                                         =
  ==========================================================================*/
@@ -858,7 +846,7 @@ void GCTAResponseCube::init_members(void)
 /***********************************************************************//**
  * @brief Copy class members
  *
- * @param[in] rsp Response to be copied
+ * @param[in] rsp Response.
  ***************************************************************************/
 void GCTAResponseCube::copy_members(const GCTAResponseCube& rsp)
 {
@@ -894,10 +882,12 @@ void GCTAResponseCube::free_members(void)
 /***********************************************************************//**
- * @brief Determines the cache index for a given model name
+ * @brief Determine cache index for model
  *
  * @param[in] name Model name.
  * @return Cache index (-1 if model has not been found).
+ *
+ * Determines the cache index for a given model @p name.
  ***************************************************************************/
 int GCTAResponseCube::cache_index(const std::string& name) const
 {
@@ -922,8 +912,114 @@ int GCTAResponseCube::cache_index(const std::string&amp; name) const
 }
+#if defined(G_RADIAL_PSF_BASED)
 /***********************************************************************//**
- * @brief Integrate PSF over radial model
+ * @brief Integrate Psf over radial model
+ *
+ * @param[in] model Radial model.
+ * @param[in] delta_mod Angle between model centre and measured photon
+ *                      direction (radians).
+ * @param[in] obsDir Observed event direction.
+ * @param[in] srcEng True photon energy.
+ * @param[in] srcTime True photon arrival time.
+ *
+ * Integrates the product of the spatial model and the point spread
+ * function over the true photon arrival direction using
+ * 
+ * \f[
+ *    \int_{\delta_{\rm min}}^{\delta_{\rm max}}
+ *    \sin \delta \times {\rm Psf}(\delta) \times
+ *    \int_{\phi_{\rm min}}^{\phi_{\rm max}} 
+ *    S_{\rm p}(\delta, \phi | E, t) d\phi d\delta
+ * \f]
+ *
+ * where
+ * \f$S_{\rm p}(\delta, \phi | E, t)\f$ is the radial spatial model,
+ * \f${\rm Psf}(\delta)\f$ is the point spread function,
+ * \f$\delta\f$ is angular distance between the true and the measured
+ * photon direction, and
+ * \f$\phi\f$ is the position angle around the observed photon direction
+ * measured counterclockwise from the connecting line between the model
+ * centre and the observed photon arrival direction.
+ ***************************************************************************/
+double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
+                                    const double&              delta_mod,
+                                    const GSkyDir&             obsDir,
+                                    const GEnergy&             srcEng,
+                                    const GTime&               srcTime) const
+{
+    // Set number of iterations for Romberg integration.
+    // These values have been determined after careful testing, see
+    // https://cta-redmine.irap.omp.eu/issues/1291
+    static const int iter_delta = 5;
+    static const int iter_phi   = 6;
+
+    // Initialise value
+    double value = 0.0;
+
+    // Get maximum Psf radius (radians)
+    double psf_max = psf().delta_max();
+
+    // Get maximum model radius (radians)
+    double theta_max = model->theta_max();
+
+    // Set offset angle integration range (radians)
+    double delta_min = (delta_mod > theta_max) ? delta_mod - theta_max : 0.0;
+    double delta_max = delta_mod + theta_max;
+    if (delta_max > psf_max) {
+        delta_max = psf_max;
+    }
+
+    // Setup integration kernel. We take here the observed photon arrival
+    // direction as the true photon arrival direction because the PSF does
+    // not vary significantly over a small region.
+    cta_psf_radial_kern_delta integrand(this,
+                                        model,
+                                        obsDir,
+                                        srcEng,
+                                        srcTime,
+                                        delta_mod,
+                                        theta_max,
+                                        iter_phi);
+
+    // Integrate over model's zenith angle
+    GIntegral integral(&integrand);
+    integral.fixed_iter(iter_delta);
+
+    // Setup integration boundaries
+    std::vector<double> bounds;
+    bounds.push_back(delta_min);
+    bounds.push_back(delta_max);
+
+    // If the integration range includes a transition between full
+    // containment of Psf within model and partial containment, then
+    // add a boundary at this location
+    double transition_point = theta_max - delta_mod;
+    if (transition_point > delta_min && transition_point < delta_max) {
+        bounds.push_back(transition_point);
+    }
+
+    // Integrate kernel
+    value = integral.romberg(bounds, iter_delta);
+
+    // Compile option: Check for NaN/Inf
+    #if defined(G_NAN_CHECK)
+    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+        std::cout << "*** ERROR: GCTAResponseCube::psf_radial:";
+        std::cout << " NaN/Inf encountered";
+        std::cout << " (value=" << value;
+        std::cout << ", delta_min=" << delta_min;
+        std::cout << ", delta_max=" << delta_max << ")";
+        std::cout << std::endl;
+    }
+    #endif
+
+    // Return integral
+    return value;
+}
+#else
+/***********************************************************************//**
+ * @brief Integrate Psf over radial model
  *
  * @param[in] model Radial model.
  * @param[in] rho_obs Angle between model centre and measured photon direction
@@ -932,21 +1028,23 @@ int GCTAResponseCube::cache_index(const std::string&amp; name) const
  * @param[in] srcEng True photon energy.
  * @param[in] srcTime True photon arrival time.
  *
- * Computes the integral
+ * Integrates the product of the spatial model and the point spread
+ * function over the true photon arrival direction using
  * 
  * \f[
- *    \int_0^{\delta_{\rm max}}
- *    {\rm PSF}(\delta) \times
- *    \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta
- *    {\rm d}\phi {\rm d}\delta
+ *    \int_{\rho_{\rm min}}^{\rho_{\rm max}}
+ *    \sin \rho \times S_{\rm p}(\rho | E, t) \times
+ *    \int_{\omega_{\rm min}}^{\omega_{\rm max}} 
+ *    {\rm Psf}(\rho, \omega) d\omega d\rho
  * \f]
  *
- * where \f${\rm M}(\delta, \phi)\f$ is a radial model in the coordinate
- * system of the point spread function, defined by the angle \f$\delta\f$
- * between the true and the measured photon direction and the azimuth angle
- * \f$\phi\f$ around the measured photon direction.
- * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread
- * function.
+ * where
+ * \f$S_{\rm p}(\rho | E, t)\f$ is the radial spatial model,
+ * \f${\rm Psf}(\rho, \omega)\f$ is the point spread function,
+ * \f$\rho\f$ is the radial distance from the model centre, and
+ * \f$\omega\f$ is the position angle around the model centre measured
+ * counterclockwise from the connecting line between the model centre and
+ * the observed photon arrival direction.
  ***************************************************************************/
 double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
                                     const double&              rho_obs,
@@ -1006,14 +1104,6 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
             bounds.push_back(transition_point);
         }
-        // If the integration range includes a transition between full
-        // containment of Psf within model and partial containment, then
-        // add a boundary at this location
-        transition_point = src_max - delta_max;
-        if (transition_point > rho_min && transition_point < rho_max) {
-            bounds.push_back(transition_point);
-        }
-
         // If we have a shell model then add an integration boundary for the
         // shell radius as a function discontinuity will occur at this
         // location
@@ -1045,10 +1135,11 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
     // Return PSF
     return irf;
 }
+#endif
 /***********************************************************************//**
- * @brief Integrate PSF over elliptical model
+ * @brief Integrate Psf over elliptical model
  *
  * @param[in] model Elliptical model.
  * @param[in] rho_obs Angle between model centre and measured photon direction
@@ -1059,36 +1150,24 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
  * @param[in] srcEng True photon energy.
  * @param[in] srcTime True photon arrival time.
  *
- * Computes the integral
+ * Integrates the product of the spatial model and the point spread
+ * function over the true photon arrival direction using
  * 
  * \f[
  *    \int_{\rho_{\rm min}}^{\rho_{\rm max}}
- *    \sin \rho \times
+ *    \sin \rho \times 
  *    \int_{\omega} 
- *    S_{\rm p}(\rho, \omega | E, t) \, PSF(\rho, \omega) d\omega d\rho
+ *    S_{\rm p}(\rho, \omega | E, t) \times
+ *    {\rm Psf}(\rho, \omega) d\omega d\rho
  * \f]
  *
  * where
- * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical model,
- * \f$PSF(\rho, \omega)\f$ is the point spread function,
- * \f$\rho\f$ is the distance from the model centre, and
- * \f$\omega\f$ is the position angle with respect to the connecting line
- * between the model centre and the observed photon arrival direction.
- ***************************************************************************/
-/*
- * \f[
- *    \int_0^{\delta_{\rm max}}
- *    {\rm PSF}(\delta) \times
- *    \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta
- *    {\rm d}\phi {\rm d}\delta
- * \f]
- *
- * where \f${\rm M}(\delta, \phi)\f$ is an elliptical model in the coordinate
- * system of the point spread function, defined by the angle \f$\delta\f$
- * between the true and the measured photon direction and the azimuth angle
- * \f$\phi\f$ around the measured photon direction.
- * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread
- * function.
+ * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical spatial model,
+ * \f${\rm Psf}(\rho, \omega)\f$ is the point spread function,
+ * \f$\rho\f$ is the radial distance from the model centre, and
+ * \f$\omega\f$ is the position angle around the model centre measured
+ * counterclockwise from the connecting line between the model centre and
+ * the observed photon arrival direction.
  ***************************************************************************/
 double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model,
                                         const double&                  rho_obs,
@@ -1268,7 +1268,7 @@ std::string GCTAResponseIrf::print(const GChatter&amp; chatter) const
 /*==========================================================================
  =                                                                         =
- =              Model type dependent CTA response methods                  =
+ =                Model type dependent CTA response methods                =
  =                                                                         =
  ==========================================================================*/
@@ -1282,24 +1282,22 @@ std::string GCTAResponseIrf::print(const GChatter&amp; chatter) const
  * @exception GCTAException::bad_model_type
  *            Model is not a radial model.
  *
- * Integrates the product of the model and the Irf over the true photon
- * arrival direction using
+ * Integrates the product of the spatial model and the instrument response
+ * function over the true photon arrival direction using
  *
  * \f[
  *    \int_{\rho_{\rm min}}^{\rho_{\rm max}}
  *    \sin \rho \times S_{\rm p}(\rho | E, t) \times
  *    \int_{\omega_{\rm min}}^{\omega_{\rm max}} 
- *    Irf(\rho, \omega) d\omega d\rho
+ *    {\rm Irf}(\rho, \omega) d\omega d\rho
  * \f]
  *
  * where
- *
- *     \f$S_{\rm p}(\rho | E, t)\f$ is the radial model,
- *     \f$Irf(\rho, \omega)\f$ is the IRF
- *     \f$\rho\f$ is the distance from the model centre, and
- *     \f$\omega\f$ is the azimuth angle is the position angle with respect
- *     to the connecting line between the model centre and the observed
- *     photon arrival direction.
+ * \f$S_{\rm p}(\rho | E, t)\f$ is the radial spatial model,
+ * \f${\rm Irf}(\rho, \omega)\f$ is the instrument response function,
+ * \f$\rho\f$ is the radial distance from the model centre, and
+ * \f$\omega\f$ is the position angle with respect to the connecting line
+ * between the model centre and the observed photon arrival direction.
  *
  * The integration is performed in the coordinate system of the source
  * model spanned by \f$\rho\f$ and \f$\omega\f$ which allows to benefit
@@ -1436,14 +1434,6 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp;       event,
             bounds.push_back(transition_point);
         }
-        // If the integration range includes a transition between full
-        // containment of Psf within model and partial containment, then
-        // add a boundary at this location
-        transition_point = src_max - delta_max;
-        if (transition_point > rho_min && transition_point < rho_max) {
-            bounds.push_back(transition_point);
-        }
-
         // If we have a shell model then add an integration boundary for the
         // shell radius as a function discontinuity will occur at this
         // location
@@ -1498,7 +1488,7 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp;       event,
 /***********************************************************************//**
- * @brief Return IRF value for elliptical source model
+ * @brief Return Irf value for elliptical source model
  *
  * @param[in] event Observed event.
  * @param[in] source Source.
@@ -692,7 +692,7 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
             // If the radius rho is not larger than the semiminor axis
             // boundary, the circle with that radius is fully contained in
             // the ellipse and we can just integrate over the relevant arc
-            if (rho < m_semiminor) {
+            if (rho <= m_semiminor) {
                 // Compute omega integration range
                 double omega_min = -domega;
@@ -984,7 +984,7 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
             // If the radius rho is not larger than the semiminor axis
             // boundary, the circle with that radius is fully contained in
             // the ellipse and we can just integrate over the relevant arc
-            if (rho < m_semiminor) {
+            if (rho <= m_semiminor) {
                 // Compute omega integration range
                 double omega_min = -domega;
@@ -1699,6 +1699,165 @@ double cta_psf_radial_kern_omega::eval(const double&amp; omega)
 /***********************************************************************//**
+ * @brief Kernel for PSF integration of radial model
+ *
+ * @param[in] delta PSF offset angle (radians).
+ * @return Azimuthally integrated product between PSF and radial model.
+ *
+ * Computes the azimuthally integrated product of point spread function and
+ * the radial model intensity. As the PSF is azimuthally symmetric, it is
+ * not included in the azimuthally integration, but just multiplied on the
+ * azimuthally integrated model. The method returns thus
+ *
+ * \f[
+ *    {\rm PSF}(\delta) \times
+ *    \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi
+ * \f]
+ *
+ * where \f${\rm M}(\delta, \phi)\f$ is the radial model in the coordinate
+ * system of the point spread function, defined by the angle \f$\delta\f$
+ * between the true and the measured photon direction and the azimuth angle
+ * \f$\phi\f$ around the measured photon direction.
+ ***************************************************************************/
+double cta_psf_radial_kern_delta::eval(const double& delta)
+{
+    // Initialise value
+    double value = 0.0;
+
+    // If we're at the Psf peak the model is zero (due to the sin(delta)
+    // term. We thus only integrate for positive deltas.
+    if (delta > 0.0) {
+
+        // Get Psf for this delta
+        double psf = m_rsp->psf()(m_srcDir, delta, m_srcEng);
+
+        // Continue only if Psf is positive
+        if (psf > 0.0) {
+
+            // Compute half length of the arc (in radians) from a circle with
+            // radius delta that intersects with the model, defined as a circle
+            // with maximum radius m_theta_max
+            double dphi = 0.5 * gammalib::cta_roi_arclength(delta,
+                                                            m_delta_mod,
+                                                            m_cos_delta_mod,
+                                                            m_sin_delta_mod,
+                                                            m_theta_max,
+                                                            m_cos_theta_max);
+
+            // Continue only if arc length is positive
+            if (dphi > 0.0) {
+
+                // Compute phi integration range
+                double phi_min = -dphi;
+                double phi_max = +dphi;
+
+                // Precompute cosine and sine terms for azimuthal integration
+                double sin_delta = std::sin(delta);
+                double cos_delta = std::cos(delta);
+                double sin_fact  = sin_delta * m_sin_delta_mod;
+                double cos_fact  = cos_delta * m_cos_delta_mod;
+
+                // Setup kernel for azimuthal integration of the spatial model
+                cta_psf_radial_kern_phi integrand(m_model,
+                                                  m_srcEng,
+                                                  m_srcTime,
+                                                  sin_fact,
+                                                  cos_fact);
+
+                // Setup integrator
+                GIntegral integral(&integrand);
+                integral.fixed_iter(m_iter);
+
+                // Integrate over azimuth
+                value = integral.romberg(phi_min, phi_max, m_iter) *
+                        psf * sin_delta;
+
+                // Debug: Check for NaN
+                #if defined(G_NAN_CHECK)
+                if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+                    std::cout << "*** ERROR: cta_psf_radial_kern_delta::eval";
+                    std::cout << "(delta=" << delta << "):";
+                    std::cout << " NaN/Inf encountered";
+                    std::cout << " (value=" << value;
+                    std::cout << ")" << std::endl;
+                }
+                #endif
+
+            } // endif: Psf value was positive
+
+        } // endif: arc length was positive
+
+    } // endif: delta was positive
+
+    // Return kernel value
+    return value;
+}
+
+
+/***********************************************************************//**
+ * @brief Kernel for azimuthal radial model integration
+ *
+ * @param[in] phi Azimuth angle (radians).
+ * @return Radial model value.
+ *
+ * Computes the value of the radial model at the position \f$(\delta,\phi)\f$
+ * given in point spread function coordinates. The \f$\theta\f$ angle of the
+ * radial model is computed using
+ *
+ * \f[
+ *    \theta = \arccos \left( \cos \delta \cos \zeta +
+ *                            \sin \delta \sin \zeta \cos \phi \right)
+ * \f]
+ *
+ * where \f$\delta\f$ is the angle between true and measured photon
+ * direction, \f$\zeta\f$ is the angle between model centre and measured
+ * photon direction, and \f$\phi\f$ is the azimuth angle with respect to the
+ * measured photon direction, where \f$\phi=0\f$ corresponds to the 
+ * connecting line between model centre and measured photon direction.
+ ***************************************************************************/
+double cta_psf_radial_kern_phi::eval(const double& phi)
+{
+    // Compute radial model theta angle
+    double theta = std::acos(m_cos_fact + m_sin_fact * std::cos(phi));
+
+    // Reduce theta by an infinite amount to avoid rounding errors at the
+    // boundary of a sharp edged model
+    double theta_kluge = theta - 1.0e-12;
+    if (theta_kluge < 0.0) {
+        theta_kluge = 0.0;
+    }
+
+    // Get radial model value
+    double value = m_model->eval(theta_kluge, m_srcEng, m_srcTime); 
+
+    // Debug: test if model is non positive
+    #if defined(G_DEBUG_MODEL_ZERO)
+    if (value <= 0.0) {
+        std::cout << "*** WARNING: cta_psf_radial_kern_phi::eval";
+        std::cout << " zero model for (phi)=(";
+        std::cout << phi*gammalib::rad2deg << ")";
+        std::cout << " theta-r_model=" << (theta-m_model->theta_max());
+        std::cout << " radians" << std::endl;
+    }
+    #endif
+
+    // Debug: Check for NaN
+    #if defined(G_NAN_CHECK)
+    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+        std::cout << "*** ERROR: cta_psf_radial_kern_phi::eval";
+        std::cout << "(phi=" << phi << "):";
+        std::cout << " NaN/Inf encountered";
+        std::cout << " (value=" << value;
+        std::cout << ")" << std::endl;
+    }
+    #endif
+
+    // Return kernel value
+    return value;
+}
+
+
+/***********************************************************************//**
  * @brief Kernel for elliptical model integration over zenith angle
  *
  * @param[in] rho Radial distance from model centre (radians).
@@ -1777,7 +1936,7 @@ double cta_psf_elliptical_kern_rho::eval(const double&amp; rho)
             // If the radius rho is not larger than the semiminor axis
             // boundary, the circle with that radius is fully contained in
             // the ellipse and we can just integrate over the relevant arc
-            if (rho < m_semiminor) {
+            if (rho <= m_semiminor) {
                 // Compute omega integration range
                 double omega_min = -domega;
@@ -1091,6 +1091,75 @@ public:
 /***********************************************************************//**
+ * @class cta_psf_radial_kern_delta
+ *
+ * @brief Kernel for Psf delta angle integration used for stacked analysis
+ ***************************************************************************/
+class cta_psf_radial_kern_delta : public GFunction {
+public:
+    cta_psf_radial_kern_delta(const GCTAResponseCube*    rsp,
+                              const GModelSpatialRadial* model,
+                              const GSkyDir&             srcDir,
+                              const GEnergy&             srcEng,
+                              const GTime&               srcTime,
+                              const double&              delta_mod,
+                              const double&              theta_max,
+                              const int&                 iter) :
+                              m_rsp(rsp),
+                              m_model(model),
+                              m_srcDir(srcDir),
+                              m_srcEng(srcEng),
+                              m_srcTime(srcTime),
+                              m_delta_mod(delta_mod),
+                              m_cos_delta_mod(std::cos(delta_mod)),
+                              m_sin_delta_mod(std::sin(delta_mod)),
+                              m_theta_max(theta_max),
+                              m_cos_theta_max(std::cos(theta_max)),
+                              m_iter(iter) { }
+    double eval(const double& delta);
+protected:
+    const GCTAResponseCube*    m_rsp;           //!< Response cube
+    const GModelSpatialRadial* m_model;         //!< Radial model
+    const GSkyDir&             m_srcDir;        //!< True photon arrival direction
+    const GEnergy&             m_srcEng;        //!< True photon energy
+    const GTime&               m_srcTime;       //!< True photon arrival time
+    const double&              m_delta_mod;     //!< Distance of model from Psf
+    double                     m_cos_delta_mod; //!< Cosine of m_delta_mod
+    double                     m_sin_delta_mod; //!< Sine of m_delta_mod
+    const double&              m_theta_max;     //!< Maximum model radius
+    double                     m_cos_theta_max; //!< Cosine of m_theta_max
+    const int&                 m_iter;          //!< Integration iterations
+};
+
+
+/***********************************************************************//**
+ * @class cta_psf_radial_kern_phi
+ *
+ * @brief Kernel for Psf phi angle integration used for stacked analysis
+ ***************************************************************************/
+class cta_psf_radial_kern_phi : public GFunction {
+public:
+    cta_psf_radial_kern_phi(const GModelSpatialRadial* model,
+                            const GEnergy&             srcEng,
+                            const GTime&               srcTime,
+                            const double&              sin_fact,
+                            const double&              cos_fact) :
+                            m_model(model),
+                            m_srcEng(srcEng),
+                            m_srcTime(srcTime),
+                            m_sin_fact(sin_fact),
+                            m_cos_fact(cos_fact) { }
+    double eval(const double& phi);
+protected:
+    const GModelSpatialRadial* m_model;     //!< Radial model
+    const GEnergy&             m_srcEng;    //!< True photon energy
+    const GTime&               m_srcTime;   //!< True photon arrival time
+    const double&              m_sin_fact;  //!< sin(delta)*sin(delta_mod)
+    const double&              m_cos_fact;  //!< cos(delta)*cos(delta_mod)
+};
+
+
+/***********************************************************************//**
  * @class cta_irf_elliptical_kern_rho
  *
  * @brief Kernel for elliptical model zenith angle integration
@@ -210,18 +210,19 @@ if __name__ == &#39;__main__&#39;:
         sys.exit()
     # Set Psf size
-    delta_max = 0.6
+    delta_max = 0.8
     # Set model parameters
     centre = gammalib.GSkyDir()
     centre.radec_deg(0.0, 0.0)
-    radius = 0.8
+    radius = 1.0
     # Set observed photon direction
     obsDir = gammalib.GSkyDir()
-    obsDir.radec_deg(0.4, -0.3)
+    #obsDir.radec_deg(0.4, -0.3)
     #obsDir.radec_deg(-0.4, -0.3)
     #obsDir.radec_deg(0.4, 0.4)
+    obsDir.radec_deg(0.0, 0.3)
     # Show elliptical response geometry
     show_response(centre, radius, obsDir, delta_max)
@@ -82,7 +82,7 @@ def show_likelihood_profile(result):
 # ========================== #
 # Compute likelihood profile #
 # ========================== #
-def compute_likelihood_profile(obs, parname="RA", scale=1.0e-4, steps=10):
+def compute_likelihood_profile(obs, parname="Radius", scale=1.0e-5, steps=10):
     """
     Compute the log-likelihood profile for a given parameter.
@@ -93,12 +93,12 @@ def compute_likelihood_profile(obs, parname=&quot;RA&quot;, scale=1.0e-4, steps=10):
     """
     # Free parameter for which profile is to be determined
     obs.models()[0][parname].free()
-    obs.models()[0]["DEC"].free()
     # Allocate Levenberg-Marquardt optimizer
     log = gammalib.GLog()
     log.cout(True)
     opt = gammalib.GOptimizerLM(log)
+    opt.eps(5.0e-3)
     # Optimize model parameters
     obs.optimize(opt)
@@ -256,8 +256,8 @@ if __name__ == &#39;__main__&#39;:
     # Set test
     #test = "point"
-    test = "gauss"
-    #test = "disk"
+    #test = "gauss"
+    test = "disk"
     #test = "ldisk"
     #test = "shell"
     #test = "ellipse"
@@ -309,10 +309,10 @@ if __name__ == &#39;__main__&#39;:
         cntref = 337.5 # +/- 1.8 (from ctobssim)
     # Perform unbinned computation
-    original_likelihood_profile(model, events, irf, gammalib.GCaldb(caldb), cntref)
+    #original_likelihood_profile(model, events, irf, gammalib.GCaldb(caldb), cntref)
     # Perform binned computation
     #original_likelihood_profile(model, cntmap, irf, gammalib.GCaldb(caldb), cntref)
     # Perform stacked computation
-    #stacked_likelihood_profile(model, cntmap, expcube, psfcube, cntref)
+    stacked_likelihood_profile(model, cntmap, expcube, psfcube, cntref)
@@ -279,7 +279,7 @@ def show_response(centre, semimajor, semiminor, posang, roiDir, pntDir, roi_radi
             if (domega > 0.0):
                 # Full containment within ellipse
-                if (rho < bound_minor):
+                if (rho <= bound_minor):
                     # Compute phi integration range
                     omega_min = -domega
@@ -400,17 +400,16 @@ if __name__ == &#39;__main__&#39;:
     # Set model parameters
     centre = gammalib.GSkyDir()
     centre.radec_deg(0.0, 0.0)
-    #semimajor = 1.0
-    #semiminor = 0.50
-    semiminor = 1.0
-    semimajor = 0.50
-    posang    = 310.0
+    semimajor = 1.0
+    semiminor = 0.50
+    posang    = 30.0
     # Set ROI direction
     roiDir = gammalib.GSkyDir()
-    #roiDir.radec_deg(-0.4, -0.3)
-    roiDir.radec_deg(+0.4, -0.3)
+    roiDir.radec_deg(-0.4, -0.3)
+    #roiDir.radec_deg(+0.4, -0.3)
     #roiDir.radec_deg(+0.4, +0.3)
+    #roiDir.radec_deg(0.0, 0.1)
     # Set pointing direction
     pntDir = gammalib.GSkyDir()