Correct transition point for shell model (was in degrees but the method expected…

… radians; this led to convergence problems with the shell model).

Correct transition point for shell model (was in degrees but the method expected…
… radians; this led to convergence problems with the shell model).
Jürgen Knödlseder
1 parent 7764d813
Showing 5 changed files with 145 additions and 1097 deletions
inst/cta/include/GCTAResponseCube.hpp
inst/cta/src/GCTAResponseCube.cpp
inst/cta/src/GCTAResponseIrf.cpp
inst/cta/src/GCTAResponse_helpers.cpp
inst/cta/src/GCTAResponse_helpers.hpp
@@ -111,7 +111,7 @@ private:
     void   free_members(void);
     int    cache_index(const std::string& name) const;
     double psf_radial(const GModelSpatialRadial* model,
-                      const double&              zeta,
+                      const double&              rho_obs,
                       const GSkyDir&             obsDir,
                       const GEnergy&             srcEng,
                       const GTime&               srcTime) const;
@@ -37,6 +37,7 @@
 #include "GCTASourceCubeDiffuse.hpp"
 #include "GModelSpatialPointSource.hpp"
 #include "GModelSpatialRadial.hpp"
+#include "GModelSpatialRadialShell.hpp"
 #include "GModelSpatialElliptical.hpp"
 #include "GModelSpatialDiffuse.hpp"
 #include "GPhoton.hpp"
@@ -63,8 +64,6 @@
 /* __ Macros _____________________________________________________________ */
 /* __ Coding definitions _________________________________________________ */
-#define G_NEW_PSF_RADIAL
-#define G_NEW_PSF_ELLIPTICAL
 /* __ Debug definitions __________________________________________________ */
@@ -464,11 +463,11 @@ double GCTAResponseCube::irf_radial(const GEvent&amp;       event,
     // Compute angle between model centre and measured photon direction
     // (radians)
-    double zeta = model->dir().dist(obsDir);
+    double rho_obs = model->dir().dist(obsDir);
     // Continue only if we're sufficiently close to the model centre to get a
     // non-zero response
-    if (zeta <= model->theta_max()+psf().delta_max()) {
+    if (rho_obs <= model->theta_max()+psf().delta_max()) {
         // Get exposure
         irf = exposure()(obsDir, obsEng);
@@ -476,7 +475,7 @@ double GCTAResponseCube::irf_radial(const GEvent&amp;       event,
         // Continue only if exposure is positive
         if (irf > 0.0) {
             irf /= obs.livetime(); //!< Recover effective area from exposure
-            irf *= psf_radial(model, zeta, obsDir, obsEng, obsTime);
+            irf *= psf_radial(model, rho_obs, obsDir, obsEng, obsTime);
             irf *= obs.deadc(obsTime);
         }
@@ -535,8 +534,8 @@ double GCTAResponseCube::irf_elliptical(const GEvent&amp;       event,
     // Compute angle between model centre and measured photon direction and
     // position angle (radians)
-    double rho_obs      = model->dir().dist(obsDir);      // was zeta
-    double posangle_obs = model->dir().posang(obsDir);    // was omega
+    double rho_obs      = model->dir().dist(obsDir);
+    double posangle_obs = model->dir().posang(obsDir);
     // Continue only if we're sufficiently close to the model centre to get a
     // non-zero response
@@ -927,8 +926,8 @@ int GCTAResponseCube::cache_index(const std::string&amp; name) const
  * @brief Integrate PSF over radial model
  *
  * @param[in] model Radial model.
- * @param[in] zeta Angle between model centre and measured photon direction
- *                 (radians).
+ * @param[in] rho_obs 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.
@@ -949,9 +948,8 @@ int GCTAResponseCube::cache_index(const std::string&amp; name) const
  * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread
  * function.
  ***************************************************************************/
-#if defined(G_NEW_PSF_ELLIPTICAL)
 double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
-                                    const double&              rho_obs, // was zeta
+                                    const double&              rho_obs,
                                     const GSkyDir&             obsDir,
                                     const GEnergy&             srcEng,
                                     const GTime&               srcTime) const
@@ -968,7 +966,7 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
     // Get maximum PSF radius (radians)
     double delta_max = psf().delta_max();
-    // Set zenith angle integration range for elliptical model (radians)
+    // Set zenith angle integration range for radial model (radians)
     double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
     double rho_max = rho_obs + delta_max;
     double src_max = model->theta_max();
@@ -1016,6 +1014,17 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
             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
+        const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(model);
+        if (shell != NULL) {
+            double shell_radius = shell->radius() * gammalib::deg2rad;
+            if (shell_radius > rho_min && shell_radius < rho_max) {
+                bounds.push_back(shell_radius);
+            }
+        }
+
         // Integrate kernel
         irf = integral.romberg(bounds, iter_rho);
@@ -1036,49 +1045,6 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
     // Return PSF
     return irf;
 }
-#else
-double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
-                                    const double&              zeta,
-                                    const GSkyDir&             obsDir,
-                                    const GEnergy&             srcEng,
-                                    const GTime&               srcTime) const
-{
-    // Set integration precision and Romberg order. These values have been
-    // determined after careful testing, see
-    // https://cta-redmine.irap.omp.eu/issues/1291
-    const double eps_delta   = 1.0e-2;
-    const int    order_delta =      5;
-    const double eps_phi     = 1.0e-2;
-    const int    order_phi   =      5;
-
-    // Initialise value
-    double value = 0.0;
-
-    // Compute sin and cosine of angle between model centre and measured
-    // photon direction
-    double sin_zeta = std::sin(zeta);
-    double cos_zeta = std::cos(zeta);
-
-    // Get offset angle integration interval in radians
-    double delta_min = 0.0;
-    double delta_max = psf().delta_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,
-                                        sin_zeta, cos_zeta,
-                                        eps_phi, order_phi);
-
-    // Integrate over PSF delta angle
-    GIntegral integral(&integrand);
-    integral.eps(eps_delta);
-    value = integral.romberg(delta_min, delta_max, order_delta);
-
-    // Return PSF
-    return value;
-}
-#endif
 /***********************************************************************//**
@@ -1124,7 +1090,6 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model,
  * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread
  * function.
  ***************************************************************************/
-#if defined(G_NEW_PSF_ELLIPTICAL)
 double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model,
                                         const double&                  rho_obs,
                                         const double&                  posangle_obs,
@@ -1201,16 +1166,6 @@ double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model,
         bounds.push_back(rho_min);
         bounds.push_back(rho_max);
-        // If the integration range includes a transition between full
-        // containment and partial containment, then add a boundary at
-        // this location
-        /*
-        double transition_point = delta_max - rho_obs;
-        if (transition_point > rho_min && transition_point < rho_max) {
-            bounds.push_back(transition_point);
-        }
-        */
-
         // If the integration range includes the semiminor boundary, then
         // add an integration boundary at that location
         if (semiminor > rho_min && semiminor < rho_max) {
@@ -1237,47 +1192,3 @@ double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model,
     // Return PSF
     return irf;
 }
-#else
-double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model,
-                                        const double&                  rho_obs,
-                                        const double&                  posangle_obs,
-                                        const GSkyDir&                 obsDir,
-                                        const GEnergy&                 srcEng,
-                                        const GTime&                   srcTime) const
-{
-    // Set integration precision and Romberg order. These values have been
-    // determined after careful testing, see
-    // https://cta-redmine.irap.omp.eu/issues/1291
-    const double eps_delta   = 1.0e-2;
-    const int    order_delta =      5;
-    const double eps_phi     = 1.0e-2;
-    const int    order_phi   =      5;
-
-    // Initialise value
-    double value = 0.0;
-
-    // Compute sin and cosine of angle between model centre and measured
-    // photon direction
-    double sin_rho_obs = std::sin(rho_obs);
-    double cos_rho_obs = std::cos(rho_obs);
-
-    // Get offset angle integration interval in radians
-    double delta_min = 0.0;
-    double delta_max = psf().delta_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_elliptical_kern_delta integrand(this, model, obsDir, srcEng, srcTime,
-                                            sin_rho_obs, cos_rho_obs, posangle_obs,
-                                            eps_phi, order_phi);
-
-    // Integrate over PSF delta angle
-    GIntegral integral(&integrand);
-    integral.eps(eps_delta);
-    value = integral.romberg(delta_min, delta_max, order_delta);
-
-    // Return PSF
-    return value;
-}
-#endif
@@ -1282,174 +1282,6 @@ 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
- *
- * \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
- * \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.
- *
- * 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
- * from the symmetry of the source model.
- *
- * The source centre is located at \f$\vec{m}\f$, and a spherical system
- * is defined around this location with \f$(\omega,\rho)\f$ being the
- * azimuth and zenith angles, respectively. \f$\omega=0\f$ is defined
- * by the direction that connects the source centre \f$\vec{m}\f$ to the
- * measured photon direction \f$\vec{p'}\f$, and \f$\omega\f$ increases
- * counterclockwise.
- *
- * Note that this method approximates the true theta angle (angle between
- * incident photon and pointing direction) by the measured theta angle
- * (angle between the measured photon arrival direction and the pointing
- * direction). Given the slow variation of the PSF shape over the field of
- * view, this approximation should be fine. It helps in fact a lot in
- * speeding up the computations.
- ***************************************************************************/
-#if defined(G_USE_PSF_SYSTEM)
-double GCTAResponseIrf::irf_radial(const GEvent&       event,
-                                   const GSource&      source,
-                                   const GObservation& obs) const
-{
-    // Set integration precision and Romberg order. These values have been
-    // determined after careful testing, see
-    // https://cta-redmine.irap.omp.eu/issues/1299
-    static const int iter_delta = 5;
-    static const int iter_phi   = 5;
-
-    // Set a tiny angle (in radians) to be removed from the model boundary
-    // to avoid requesting values beyond the boundary due to round-off
-    // errors (this could produce a discontinuity that the integrator does
-    // not like)
-    static const double tiny = 1.0e-9;
-
-    // Initialise Irf
-    double irf = 0.0;
-
-    // Retrieve CTA pointing
-    const GCTAPointing& pnt = retrieve_pnt(G_IRF_RADIAL, obs);
-    const GCTAInstDir&  dir = retrieve_dir(G_IRF_RADIAL, event);
-
-    // Get pointer to radial model
-    const GModelSpatialRadial* model = static_cast<const GModelSpatialRadial*>(source.model());
-
-    // Get event attributes
-    const GSkyDir& obsDir = dir.dir();
-    const GEnergy& obsEng = event.energy();
-
-    // Get source attributes
-    const GEnergy& srcEng    = source.energy();
-    const GTime&   srcTime   = source.time();
-    double         srcLogEng = srcEng.log10TeV();
-
-    // Angle between measured photon and model centre (radians)
-    double zeta = obsDir.dist(model->dir());
-
-    // Angle between measured photon and pointing (radians)
-    double obsOffset = obsDir.dist(pnt.dir());
-
-    // Compute position angle of pointing in Psf system (radians).
-    // The angle Phi=0 of the Psf system is given by the connecting line
-    // between the measured photon direction and the model centre
-    double phi0 = obsDir.posang(pnt.dir()) - obsDir.posang(model->dir());
-    
-    // Get Psf delta integration interval [delta_min, delta_max] (radians)
-    double theta_max = model->theta_max() - tiny;
-    double delta_min = (zeta > theta_max) ? zeta - theta_max : 0.0;
-    double delta_max = psf_delta_max(obsOffset, 0.0,
-                                     pnt.zenith(), pnt.azimuth(),
-                                     srcLogEng);
-    double distance  = zeta + theta_max;
-    if (distance < delta_max) {
-        delta_max = distance;
-    }
-
-    // Continue only if interval has positive length
-    if (delta_max > delta_min) {
-
-        // Setup integration kernel.
-        cta_irf_radial_kern_delta integrand(this, model, pnt,
-                                            srcEng, srcLogEng, srcTime,
-                                            obsOffset, obsEng,
-                                            zeta, phi0, theta_max,
-                                            iter_phi);
-
-        // Setup integration over Psf delta 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 and partial containment, then add a boundary at
-        // this location
-        double transition_point = theta_max - zeta;
-        if (transition_point > delta_min && transition_point < delta_max) {
-            bounds.push_back(transition_point);
-        }
-
-        // Integrate kernel
-        irf = integral.romberg(bounds, iter_delta);
-
-        // Compile option: Check for NaN/Inf
-        #if defined(G_NAN_CHECK)
-        if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
-            std::cout << "*** ERROR: GCTAResponseIrf::irf_radial:";
-            std::cout << " NaN/Inf encountered";
-            std::cout << " (irf=" << irf;
-            std::cout << ", delta_min=" << delta_min;
-            std::cout << ", delta_max=" << delta_max << ")";
-            std::cout << std::endl;
-        }
-        #endif
-
-        // Apply deadtime correction
-        irf *= obs.deadc(srcTime);
-
-    } // endif: delta interval has positive length
-
-    // Compile option: Show integration results
-    #if defined(G_DEBUG_IRF_RADIAL)
-    std::cout << "GCTAResponseIrf::irf_radial:";
-    std::cout << " delta=[" << delta_min*gammalib::rad2deg;
-    std::cout << ", " << delta_max*gammalib::rad2deg << " deg;";
-    std::cout << " zeta=" << zeta*gammalib::rad2deg << " deg;";
-    std::cout << " phi0=" << phi0*gammalib::rad2deg << " deg;";
-    std::cout << " theta_max=" << theta_max*gammalib::rad2deg << " deg;";
-    std::cout << " obsDir=" << obsDir << ";";
-    std::cout << " modelDir=" << model->dir() << ";";
-    std::cout << " irf=" << irf << std::endl;
-    #endif
-
-    // Return Irf value
-    return irf;
-}
-#else
-/***********************************************************************//**
- * @brief Return IRF value for radial source model
- *
- * @param[in] event Observed event.
- * @param[in] source Source.
- * @param[in] obs Observation.
- *
- * @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
  *
@@ -1617,7 +1449,7 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp;       event,
         // location
         const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(model);
         if (shell != NULL) {
-            double shell_radius = shell->radius();
+            double shell_radius = shell->radius() * gammalib::deg2rad;
             if (shell_radius > rho_min && shell_radius < rho_max) {
                 bounds.push_back(shell_radius);
             }
@@ -1663,7 +1495,6 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp;       event,
     // Return IRF value
     return irf;
 }
-#endif
 /***********************************************************************//**
@@ -1824,16 +1655,6 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp;       event,
         bounds.push_back(rho_min);
         bounds.push_back(rho_max);
-        // If the integration range includes a transition between full
-        // containment and partial containment, then add a boundary at
-        // this location
-        /*
-        double transition_point = delta_max - rho_obs;
-        if (transition_point > rho_min && transition_point < rho_max) {
-            bounds.push_back(transition_point);
-        }
-        */
-
         // If the integration range includes the semiminor boundary, then
         // add an integration boundary at that location
         if (semiminor > rho_min && semiminor < rho_max) {
@@ -2224,7 +2045,7 @@ double GCTAResponseIrf::npred_radial(const GSource&amp; source,
         // location
         const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(model);
         if (shell != NULL) {
-            double shell_radius = shell->radius();
+            double shell_radius = shell->radius() * gammalib::deg2rad;
             if (shell_radius > rho_min && shell_radius < rho_max) {
                 bounds.push_back(shell_radius);
             }
@@ -283,8 +283,15 @@ double cta_irf_radial_kern_rho::eval(const double&amp; rho)
             double omega_min = -domega;
             double omega_max = +domega;
+            // Reduce rho by an infinite amount to avoid rounding errors
+            // at the boundary of a sharp edged model
+            double rho_kluge = rho - 1.0e-12;
+            if (rho_kluge < 0.0) {
+                rho_kluge = 0.0;
+            }
+
             // Evaluate sky model
-            double model = m_model.eval(rho, m_srcEng, m_srcTime);
+            double model = m_model.eval(rho_kluge, m_srcEng, m_srcTime);
             // Debug: test if model is non positive
             #if defined(G_DEBUG_MODEL_ZERO)
@@ -292,7 +299,7 @@ double cta_irf_radial_kern_rho::eval(const double&amp; rho)
                 std::cout << "*** WARNING: cta_irf_radial_kern_rho::eval";
                 std::cout << " zero model for (rho)=(";
                 std::cout << rho*gammalib::rad2deg << ")";
-                std::cout << " rho-r_model" << (rho-m_model.theta_max());
+                std::cout << " rho-r_model=" << (rho-m_model.theta_max());
                 std::cout << " radians" << std::endl;
             }
             #endif
@@ -474,8 +481,15 @@ double cta_npred_radial_kern_rho::eval(const double&amp; rho)
             double omega_min = m_omega0 - domega;
             double omega_max = m_omega0 + domega;
+            // Reduce rho by an infinite amount to avoid rounding errors
+            // at the boundary of a sharp edged model
+            double rho_kluge = rho - 1.0e-12;
+            if (rho_kluge < 0.0) {
+                rho_kluge = 0.0;
+            }
+
             // Get radial model value
-            double model = m_model.eval(rho, m_srcEng, m_srcTime);
+            double model = m_model.eval(rho_kluge, m_srcEng, m_srcTime);
             // Debug: test if model is non positive
             #if defined(G_DEBUG_MODEL_ZERO)
@@ -483,7 +497,7 @@ double cta_npred_radial_kern_rho::eval(const double&amp; rho)
                 std::cout << "*** WARNING: cta_npred_radial_kern_rho::eval";
                 std::cout << " zero model for (rho)=(";
                 std::cout << rho*gammalib::rad2deg << ")";
-                std::cout << " rho-r_model" << (rho-m_model.theta_max());
+                std::cout << " rho-r_model=" << (rho-m_model.theta_max());
                 std::cout << " radians" << std::endl;
             }
             #endif
@@ -1500,237 +1514,11 @@ double cta_npred_diffuse_kern_phi::eval(const double&amp; phi)
 }
-/***********************************************************************//**
- * @brief Kernel for PSF integration of spatial model
- *
- * @param[in] delta PSF offset angle (radians).
- * @return Azimuthally integrated product between PSF and model.
- *
- * Computes the azimuthally integrated product of point spread function and
- * the spatial 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 spatial 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_diffuse_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
-        value = m_rsp->psf()(m_srcDir, delta, m_srcEng);
-
-        // Continue only if PSF is positive
-        if (value > 0.0) {
-
-            // Compute sine and cosine of delta
-            double sin_delta = std::sin(delta);
-            double cos_delta = std::cos(delta);
-
-            // Setup kernel for azimuthal integration of the diffuse model
-            cta_psf_diffuse_kern_phi integrand(m_model, m_srcEng, m_srcTime,
-                                               m_rot,
-                                               sin_delta, cos_delta);
-
-            // Azimuthally integrate model
-            GIntegral integral(&integrand);
-            integral.eps(m_eps);
-            value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
-
-        } // endif: PSF value was positive
-
-    } // endif: delta was positive
-
-    // Debug: Check for NaN
-    #if defined(G_NAN_CHECK)
-    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: cta_psf_diffuse_kern_delta::eval";
-        std::cout << "(delta=" << delta << "):";
-        std::cout << " NaN/Inf encountered";
-        std::cout << " (value=" << value;
-        std::cout << ")" << std::endl;
-    }
-    #endif
-
-    // Return kernel value
-    return value;
-}
-
-
-/***********************************************************************//**
- * @brief Kernel for map integration of spatial model
- *
- * @param[in] phi Azimuth angle (radians).
- * @return Spatial model value.
- *
- * Computes the value of the spatial model at the position (delta,phi) given
- * in point spread function coordinates. The transformation from point
- * spread function coordinates into sky coordinates is done using a rotation
- * matrix that is pre-computed on entry.
- ***************************************************************************/
-double cta_psf_diffuse_kern_phi::eval(const double& phi)
-{
-    // Compute sky direction vector in native coordinates
-    double  cos_phi = std::cos(phi);
-    double  sin_phi = std::sin(phi);
-    GVector native(-cos_phi*m_sin_delta, sin_phi*m_sin_delta, m_cos_delta);
-
-    // Rotate from native into celestial system
-    GVector cel = m_rot * native;
-
-    // Set sky direction
-    GSkyDir srcDir;
-    srcDir.celvector(cel);
-
-    // Set photon
-    GPhoton photon(srcDir, m_srcEng, m_srcTime);
-
-    // Compute map value this sky direction
-    double value = m_model->eval(photon); 
-
-    // Debug: Check for NaN
-    #if defined(G_NAN_CHECK)
-    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: cta_psf_diffuse_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 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
-        value = m_rsp->psf()(m_srcDir, delta, m_srcEng);
-
-        // Continue only if PSF is positive
-        if (value > 0.0) {
-
-            // Compute sin(delta), sin(delta)*sin(zeta), cos(delta)*cos(zeta)
-            double sin_delta = std::sin(delta);
-            double sin_fact  = sin_delta * m_sin_zeta;
-            double cos_fact  = std::cos(delta)* m_cos_zeta;
-
-            // 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);
-
-            // Azimuthally integrate model
-            GIntegral integral(&integrand);
-            integral.eps(m_eps);
-            value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
-
-        } // endif: PSF value was positive
-
-    } // endif: delta was positive
-
-    // 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
-
-    // 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
- * \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));
-
-    // Get radial model value
-    double value = m_model->eval(theta, m_srcEng, m_srcTime); 
-
-    // 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;
-}
-
+/*==========================================================================
+ =                                                                         =
+ =                Helper class methods for stacked analysis                =
+ =                                                                         =
+ ==========================================================================*/
 /***********************************************************************//**
  * @brief Kernel for radial model integration over zenith angle
@@ -1777,8 +1565,7 @@ double cta_psf_radial_kern_rho::eval(const double&amp; rho)
         if (domega > 0.0) {
             // Reduce rho by an infinite amount to avoid rounding errors
-            // at the boundary of a sharp edged model (e.g. an elliptical
-            // disk model)
+            // at the boundary of a sharp edged model
             double rho_kluge = rho - 1.0e-12;
             if (rho_kluge < 0.0) {
                 rho_kluge = 0.0;
@@ -1793,7 +1580,7 @@ double cta_psf_radial_kern_rho::eval(const double&amp; rho)
                 std::cout << "*** WARNING: cta_psf_radial_kern_rho::eval";
                 std::cout << " zero model for (rho)=(";
                 std::cout << rho*gammalib::rad2deg << ")";
-                std::cout << " rho-r_model" << (rho-m_model->theta_max());
+                std::cout << " rho-r_model=" << (rho-m_model->theta_max());
                 std::cout << " radians" << std::endl;
             }
             #endif
@@ -1912,133 +1699,6 @@ double cta_psf_radial_kern_omega::eval(const double&amp; omega)
 /***********************************************************************//**
- * @brief Kernel for PSF integration of elliptical model
- *
- * @param[in] delta PSF offset angle (radians).
- * @return Azimuthally integrated product between PSF and elliptical model.
- *
- * Computes the azimuthally integrated product of point spread function and
- * the elliptical 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 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.
- ***************************************************************************/
-double cta_psf_elliptical_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
-        value = m_rsp->psf()(m_srcDir, delta, m_srcEng);
-
-        // Continue only if PSF is positive
-        if (value > 0.0) {
-
-            // Compute sin(delta), sin(delta)*sin(zeta), cos(delta)*cos(zeta)
-            double sin_delta = std::sin(delta);
-            double sin_fact  = sin_delta * m_sin_zeta;
-            double cos_fact  = std::cos(delta)* m_cos_zeta;
-
-            // Setup kernel for azimuthal integration of the elliptical model
-            cta_psf_elliptical_kern_phi integrand(m_model, m_srcEng, m_srcTime,
-                                                  m_omega,
-                                                  sin_delta, sin_fact, cos_fact);
-
-            // Azimuthally integrate model
-            GIntegral integral(&integrand);
-            integral.eps(m_eps);
-            value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
-
-        } // endif: PSF value was positive
-
-    } // endif: delta was positive
-
-    // Debug: Check for NaN
-    #if defined(G_NAN_CHECK)
-    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: cta_psf_elliptical_kern_delta::eval";
-        std::cout << "(delta=" << delta << "):";
-        std::cout << " NaN/Inf encountered";
-        std::cout << " (value=" << value;
-        std::cout << ")" << std::endl;
-    }
-    #endif
-
-    // Return kernel value
-    return value;
-}
-
-
-/***********************************************************************//**
- * @brief Kernel for azimuthal elliptical model integration
- *
- * @param[in] phi Azimuth angle (radians).
- * @return Radial model value.
- *
- * Computes the value of the elliptical model at the position \f$(\delta,\phi)\f$
- * given in point spread function coordinates. The \f$\theta\f$ angle of the
- * elliptical 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.
- *
- * The position angle of the elliptical model is computed using
- *
- * \f[
- *    {\rm posangle} = \arcsin \left( \frac{\sin \phi \sin \delta}
- *                                         {\sin \theta} \right)
- * \f]
- ***************************************************************************/
-double cta_psf_elliptical_kern_phi::eval(const double& phi)
-{
-    // Compute radial model theta angle
-    double theta = std::acos(m_cos_fact + m_sin_fact * std::cos(phi));
-
-    // Compute radial model position angle
-    double posangle = std::asin((std::sin(phi) * m_sin_delta)/std::sin(theta)) +
-                      m_omega;
-
-    // Get radial model value
-    double value = m_model->eval(theta, posangle, m_srcEng, m_srcTime);
-
-    // Debug: Check for NaN
-    #if defined(G_NAN_CHECK)
-    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: cta_psf_elliptical_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).
@@ -2305,112 +1965,64 @@ double cta_psf_elliptical_kern_omega::eval(const double&amp; omega)
 }
-/*==========================================================================
- =                                                                         =
- =          Helper class methods for response computation testing          =
- =                                                                         =
- = WARNING: These helper classes are for testing purposes only and not for =
- =          production. These classes may be removed at any moment without =
- =          any warning.                                                   =
- ==========================================================================*/
-
-
 /***********************************************************************//**
- * @brief Kernel for integration of radial model in Psf system
+ * @brief Kernel for PSF integration of spatial model
  *
- * @param[in] delta Psf offset angle (radians).
- * @return Azimuthally integrated product between Psf and radial model.
+ * @param[in] delta PSF offset angle (radians).
+ * @return Azimuthally integrated product between PSF and model.
  *
  * Computes the azimuthally integrated product of point spread function and
- * the radial model intensity. As the Psf is azimuthally symmetric, it is
+ * the spatial 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
+ *    {\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
+ * where \f${\rm M}(\delta, \phi)\f$ is the spatial 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_irf_radial_kern_delta::eval(const double& delta)
+double cta_psf_diffuse_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)
+    // 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. We use here the event offset instead of
-        // the true photon offset because the Psf varies little over the
-        // dimension of the Psf, hence it should basically be invariant
-        // over such small scales
-        value = m_rsp->psf(delta, m_obsOffset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng);
-
-        // Compile option: Get Aeff now using the event offset instead of
-        // computing Aeff later in the azimuthal integration kernel where
-        // the true photon offset is available. This is an approximation and
-        // assumes that Aeff varies little over the dimension of the Psf.
-        // This option is only for testing purposes and should not be used
-        // for production.
-        #if defined(G_USE_OBSDIR_FOR_AEFF)
-        value *= m_rsp->aeff(m_obsOffset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng);
-        #endif
+        // Get PSF for this delta
+        value = m_rsp->psf()(m_srcDir, delta, m_srcEng);
-        // Continue only if Psf is positive
+        // Continue only if PSF is positive
         if (value > 0.0) {
-            // Compute half length of arc that lies within model (radians)
-            double dphi = 0.5 * gammalib::cta_roi_arclength(delta,
-                                                            m_zeta,
-                                                            m_cos_zeta,
-                                                            m_sin_zeta,
-                                                            m_theta_max,
-                                                            m_cos_theta_max);
-
-            // Continue only if arc length is positive
-            if (dphi > 0.0) {
-
-                // Compute azimuth integration range [phi_min,phi_max].
-                double phi_min = -dphi;
-                double phi_max = +dphi;
-
-                // Pre-compute terms for angular distance computation in
-                // azimuth integration kernel.
-                double sin_delta      = std::sin(delta);
-                double cos_delta      = std::cos(delta);
-                double sin_fact_model = sin_delta * m_sin_zeta;
-                double cos_fact_model = cos_delta * m_cos_zeta;
-                double sin_fact_inst  = sin_delta * m_sin_obs_offset;
-                double cos_fact_inst  = cos_delta * m_cos_obs_offset;
-
-                // Setup kernel for azimuthal integration of the spatial model
-                cta_irf_radial_kern_phi integrand(m_rsp, m_model, m_pnt,
-                                                  m_srcEng, m_srcLogEng, m_srcTime,
-                                                  m_obsEng,
-                                                  m_phi0,
-                                                  sin_fact_model, cos_fact_model,
-                                                  sin_fact_inst,  cos_fact_inst);
-
-                // Integrate kernel
-                GIntegral integral(&integrand);
-                integral.fixed_iter(m_iter);
-                value *= integral.romberg(phi_min, phi_max, m_iter) * sin_delta;
+            // Compute sine and cosine of delta
+            double sin_delta = std::sin(delta);
+            double cos_delta = std::cos(delta);
+
+            // Setup kernel for azimuthal integration of the diffuse model
+            cta_psf_diffuse_kern_phi integrand(m_model, m_srcEng, m_srcTime,
+                                               m_rot,
+                                               sin_delta, cos_delta);
-            } // endif: arc length was positive
+            // Azimuthally integrate model
+            GIntegral integral(&integrand);
+            integral.eps(m_eps);
+            value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
-        } // endif: Psf value was positive
+        } // endif: PSF value was positive
     } // endif: delta was positive
     // Debug: Check for NaN
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: cta_irf_radial_kern_delta::eval";
+        std::cout << "*** ERROR: cta_psf_diffuse_kern_delta::eval";
         std::cout << "(delta=" << delta << "):";
         std::cout << " NaN/Inf encountered";
         std::cout << " (value=" << value;
@@ -2424,79 +2036,40 @@ double cta_irf_radial_kern_delta::eval(const double&amp; delta)
 /***********************************************************************//**
- * @brief Kernel for azimuthal radial model integration
+ * @brief Kernel for map integration of spatial model
  *
  * @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
- * \f]
+ * @return Spatial model value.
  *
- * 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.
+ * Computes the value of the spatial model at the position (delta,phi) given
+ * in point spread function coordinates. The transformation from point
+ * spread function coordinates into sky coordinates is done using a rotation
+ * matrix that is pre-computed on entry.
  ***************************************************************************/
-double cta_irf_radial_kern_phi::eval(const double& phi)
+double cta_psf_diffuse_kern_phi::eval(const double& phi)
 {
-    // Compute radial model theta angle
-    double theta = gammalib::acos(m_cos_fact_model +
-                                  m_sin_fact_model * std::cos(phi));
-
-    // Get radial model value
-    double value = m_model->eval(theta, m_srcEng, m_srcTime);
-
-    // For debugging: Check boundary violation
-    if (value == 0.0) {
-        std::cout << "Model 0 at theta=" << theta*gammalib::rad2deg;
-        std::cout << " deg; phi=";
-        std::cout << phi*gammalib::rad2deg << " deg";
-        std::cout << std::endl;
-        throw;
-    }
-
-    // If model is positive the multiply it with the effective area and
-    // optionally fold in energy dispersion
-    #if defined(G_USE_OBSDIR_FOR_AEFF)
-    if (m_rsp->use_edisp() && value > 0.0) {
-
-        // Compute offset angle of true photon in camera system (radians)
-        double offset = gammalib::acos(m_cos_fact_inst +
-                                       m_sin_fact_inst * std::cos(m_phi0 - phi));
-
-        // Take energy dispersion into account
-        value *= m_rsp->edisp(m_obsEng, offset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng);
+    // Compute sky direction vector in native coordinates
+    double  cos_phi = std::cos(phi);
+    double  sin_phi = std::sin(phi);
+    GVector native(-cos_phi*m_sin_delta, sin_phi*m_sin_delta, m_cos_delta);
-    } // endif: model value was positive
-    #else
-    if (value > 0.0) {
+    // Rotate from native into celestial system
+    GVector cel = m_rot * native;
-        // Compute offset angle of true photon in camera system (radians)
-        double offset = gammalib::acos(m_cos_fact_inst +
-                                       m_sin_fact_inst * std::cos(m_phi0 - phi));
+    // Set sky direction
+    GSkyDir srcDir;
+    srcDir.celvector(cel);
-        // Multiply with effective area
-        value *= m_rsp->aeff(offset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng);
+    // Set photon
+    GPhoton photon(srcDir, m_srcEng, m_srcTime);
-        // Optionally take energy dispersion into account
-        if (m_rsp->use_edisp() && value > 0.0) {
-            value *= m_rsp->edisp(m_obsEng, offset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng);
-        }
+    // Compute map value this sky direction
+    double value = m_model->eval(photon); 
-    } // endif: model value was positive
-    #endif
-    
     // Debug: Check for NaN
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
-        std::cout << "*** ERROR: cta_irf_radial_kern_phi::eval";
+        std::cout << "*** ERROR: cta_psf_diffuse_kern_phi::eval";
         std::cout << "(phi=" << phi << "):";
         std::cout << " NaN/Inf encountered";
         std::cout << " (value=" << value;
@@ -982,142 +982,6 @@ protected:
 /***********************************************************************//**
- * @class cta_psf_diffuse_kern_delta
- *
- * @brief Kernel for Psf delta angle integration used for stacked analysis
- ***************************************************************************/
-class cta_psf_diffuse_kern_delta : public GFunction {
-public:
-    cta_psf_diffuse_kern_delta(const GCTAResponseCube* rsp,
-                               const GModelSpatial*    model,
-                               const GSkyDir&          srcDir,
-                               const GEnergy&          srcEng,
-                               const GTime&            srcTime,
-                               const GMatrix&          rot,
-                               const double&           eps,
-                               const int&              order) :
-                               m_rsp(rsp),
-                               m_model(model),
-                               m_srcDir(srcDir),
-                               m_srcEng(srcEng),
-                               m_srcTime(srcTime),
-                               m_rot(rot),
-                               m_eps(eps),
-                               m_order (order),
-                               m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { }
-    double eval(const double& delta);
-protected:
-    const GCTAResponseCube* m_rsp;     //!< Response cube
-    const GModelSpatial*    m_model;   //!< Spatial 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 GMatrix&          m_rot;     //!< Rotation matrix
-    const double&           m_eps;     //!< Integration precision
-    const int&              m_order;   //!< Romberg order
-    double                  m_psf_max; //!< Maximum PSF value
-};
-
-
-/***********************************************************************//**
- * @class cta_psf_diffuse_kern_phi
- *
- * @brief Kernel for Psf phi angle integration used for stacked analysis
- ***************************************************************************/
-class cta_psf_diffuse_kern_phi : public GFunction {
-public:
-    cta_psf_diffuse_kern_phi(const GModelSpatial* model,
-                             const GEnergy&       srcEng,
-                             const GTime&         srcTime,
-                             const GMatrix&       rot,
-                             const double&        sin_delta,
-                             const double&        cos_delta) :
-                             m_model(model),
-                             m_srcEng(srcEng),
-                             m_srcTime(srcTime),
-                             m_rot(rot),
-                             m_sin_delta(sin_delta),
-                             m_cos_delta(cos_delta) { }
-    double eval(const double& phi);
-protected:
-    const GModelSpatial* m_model;     //!< Spatial model
-    const GEnergy&       m_srcEng;    //!< True photon energy
-    const GTime&         m_srcTime;   //!< True photon arrival time
-    const GMatrix&       m_rot;       //!< Rotation matrix
-    const double&        m_sin_delta; //!< sin(delta)
-    const double&        m_cos_delta; //!< cos(delta)
-};
-
-
-/***********************************************************************//**
- * @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&              sin_zeta,
-                              const double&              cos_zeta,
-                              const double&              eps,
-                              const int&                 order) :
-                              m_rsp(rsp),
-                              m_model(model),
-                              m_srcDir(srcDir),
-                              m_srcEng(srcEng),
-                              m_srcTime(srcTime),
-                              m_sin_zeta(sin_zeta),
-                              m_cos_zeta(cos_zeta),
-                              m_eps(eps),
-                              m_order (order),
-                              m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { }
-    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_sin_zeta; //!< sin(zeta)
-    const double&              m_cos_zeta; //!< cos(zeta)
-    const double&              m_eps;      //!< Integration precision
-    const int&                 m_order;    //!< Romberg order
-    double                     m_psf_max;  //!< Maximum PSF value
-};
-
-
-/***********************************************************************//**
- * @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(zeta)
-    const double&              m_cos_fact;  //!< cos(delta)*cos(zeta)
-};
-
-
-/***********************************************************************//**
  * @class cta_irf_radial_kern_rho
  *
  * @brief Kernel for radial model zenith angle integration
@@ -1227,83 +1091,6 @@ public:
 /***********************************************************************//**
- * @class cta_psf_elliptical_kern_delta
- *
- * @brief Kernel for Psf delta angle integration used for stacked analysis
- ***************************************************************************/
-class cta_psf_elliptical_kern_delta : public GFunction {
-public:
-    cta_psf_elliptical_kern_delta(const GCTAResponseCube*        rsp,
-                                  const GModelSpatialElliptical* model,
-                                  const GSkyDir&                 srcDir,
-                                  const GEnergy&                 srcEng,
-                                  const GTime&                   srcTime,
-                                  const double&                  sin_zeta,
-                                  const double&                  cos_zeta,
-                                  const double&                  omega,
-                                  const double&                  eps,
-                                  const int&                     order) :
-                                  m_rsp(rsp),
-                                  m_model(model),
-                                  m_srcDir(srcDir),
-                                  m_srcEng(srcEng),
-                                  m_srcTime(srcTime),
-                                  m_sin_zeta(sin_zeta),
-                                  m_cos_zeta(cos_zeta),
-                                  m_omega(omega),
-                                  m_eps(eps),
-                                  m_order (order),
-                                  m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { }
-    double eval(const double& delta);
-protected:
-    const GCTAResponseCube*        m_rsp;       //!< Response cube
-    const GModelSpatialElliptical* m_model;     //!< Elliptical 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_sin_zeta;  //!< sin(zeta)
-    const double&                  m_cos_zeta;  //!< cos(zeta)
-    const double&                  m_omega;     //!< Position angle of photon w/r model
-    const double&                  m_eps;       //!< Integration precision
-    const int&                     m_order;     //!< Romberg order
-    double                         m_psf_max;   //!< Maximum PSF value
-};
-
-
-/***********************************************************************//**
- * @class cta_psf_elliptical_kern_delta
- *
- * @brief Kernel for Psf phi angle integration used for stacked analysis
- ***************************************************************************/
-class cta_psf_elliptical_kern_phi : public GFunction {
-public:
-    cta_psf_elliptical_kern_phi(const GModelSpatialElliptical* model,
-                                const GEnergy&                 srcEng,
-                                const GTime&                   srcTime,
-                                const double&                  omega,
-                                const double&                  sin_delta,
-                                const double&                  sin_fact,
-                                const double&                  cos_fact) :
-                                m_model(model),
-                                m_srcEng(srcEng),
-                                m_srcTime(srcTime),
-                                m_omega(omega),
-                                m_sin_delta(sin_delta),
-                                m_sin_fact(sin_fact),
-                                m_cos_fact(cos_fact) { }
-    double eval(const double& phi);
-protected:
-    const GModelSpatialElliptical* m_model;     //!< Elliptical model
-    const GEnergy&                 m_srcEng;    //!< True photon energy
-    const GTime&                   m_srcTime;   //!< True photon arrival time
-    const double&                  m_omega;     //!< Position angle of photon w/r model
-    const double&                  m_sin_delta; //!< sin(delta)
-    const double&                  m_sin_fact;  //!< sin(delta)*sin(zeta)
-    const double&                  m_cos_fact;  //!< cos(delta)*cos(zeta)
-};
-
-
-/***********************************************************************//**
  * @class cta_irf_elliptical_kern_rho
  *
  * @brief Kernel for elliptical model zenith angle integration
@@ -1433,114 +1220,70 @@ public:
 /***********************************************************************//**
- * @class cta_irf_radial_kern_delta
- *
- * @brief Test kernel for Irf delta angle integration
+ * @class cta_psf_diffuse_kern_delta
  *
- * This class is used for Psf system based integration. It is only here for
- * testing and this class may be removed in the future. 
+ * @brief Kernel for Psf delta angle integration used for stacked analysis
  ***************************************************************************/
-class cta_irf_radial_kern_delta : public GFunction {
+class cta_psf_diffuse_kern_delta : public GFunction {
 public:
-    cta_irf_radial_kern_delta(const GCTAResponseIrf*     rsp,
-                              const GModelSpatialRadial* model,
-                              const GCTAPointing&        pnt,
-                              const GEnergy&             srcEng,
-                              const double&              srcLogEng,
-                              const GTime&               srcTime,
-                              const double&              obsOffset,
-                              const GEnergy&             obsEng,
-                              const double&              zeta,
-                              const double&              phi0,
-                              const double&              theta_max,
-                              const int&                 iter) :
-                              m_rsp(rsp),
-                              m_model(model),
-                              m_pnt(pnt),
-                              m_srcEng(srcEng),
-                              m_srcLogEng(srcLogEng),
-                              m_srcTime(srcTime),
-                              m_obsOffset(obsOffset),
-                              m_sin_obs_offset(std::sin(obsOffset)),
-                              m_cos_obs_offset(std::cos(obsOffset)),
-                              m_obsEng(obsEng),
-                              m_zeta(zeta),
-                              m_phi0(phi0),
-                              m_sin_zeta(std::sin(zeta)),
-                              m_cos_zeta(std::cos(zeta)),
-                              m_theta_max(theta_max),
-                              m_cos_theta_max(std::cos(theta_max)),
-                              m_iter(iter) { }
+    cta_psf_diffuse_kern_delta(const GCTAResponseCube* rsp,
+                               const GModelSpatial*    model,
+                               const GSkyDir&          srcDir,
+                               const GEnergy&          srcEng,
+                               const GTime&            srcTime,
+                               const GMatrix&          rot,
+                               const double&           eps,
+                               const int&              order) :
+                               m_rsp(rsp),
+                               m_model(model),
+                               m_srcDir(srcDir),
+                               m_srcEng(srcEng),
+                               m_srcTime(srcTime),
+                               m_rot(rot),
+                               m_eps(eps),
+                               m_order (order),
+                               m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { }
     double eval(const double& delta);
 protected:
-    const GCTAResponseIrf*     m_rsp;       //!< Response
-    const GModelSpatialRadial* m_model;     //!< Radial model
-    const GCTAPointing&        m_pnt;       //!< CTA pointing
-    const GEnergy&             m_srcEng;    //!< True photon energy
-    const double&              m_srcLogEng; //!< Log10(true photon energy)
-    const GTime&               m_srcTime;   //!< True photon arrival time
-    const double&              m_obsOffset; //!< Observed photon offset angle
-    const double               m_sin_obs_offset;
-    const double               m_cos_obs_offset;
-    const GEnergy&             m_obsEng;    //!< Measured event energy
-    const double&              m_zeta;      //!< zeta
-    const double&              m_phi0;      //!< phi0
-    const double               m_sin_zeta;  //!< sin(zeta)
-    const double               m_cos_zeta;  //!< cos(zeta)
-    const double&              m_theta_max; //!< Model radius
-    const double               m_cos_theta_max;
-    const int&                 m_iter;      //!< Integration iterations
+    const GCTAResponseCube* m_rsp;     //!< Response cube
+    const GModelSpatial*    m_model;   //!< Spatial 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 GMatrix&          m_rot;     //!< Rotation matrix
+    const double&           m_eps;     //!< Integration precision
+    const int&              m_order;   //!< Romberg order
+    double                  m_psf_max; //!< Maximum PSF value
 };
 /***********************************************************************//**
- * @class cta_irf_radial_kern_phi
- *
- * @brief Test kernel for Irf phi angle integration
+ * @class cta_psf_diffuse_kern_phi
  *
- * This class is used for Psf system based integration. It is only here for
- * testing and this class may be removed in the future. 
+ * @brief Kernel for Psf phi angle integration used for stacked analysis
  ***************************************************************************/
-class cta_irf_radial_kern_phi : public GFunction {
+class cta_psf_diffuse_kern_phi : public GFunction {
 public:
-    cta_irf_radial_kern_phi(const GCTAResponseIrf*     rsp,
-                            const GModelSpatialRadial* model,
-                            const GCTAPointing&        pnt,
-                            const GEnergy&             srcEng,
-                            const double&              srcLogEng,
-                            const GTime&               srcTime,
-                            const GEnergy&             obsEng,
-                            const double&              phi0,
-                            const double&              sin_fact_model,
-                            const double&              cos_fact_model,
-                            const double&              sin_fact_inst,
-                            const double&              cos_fact_inst) :
-                            m_rsp(rsp),
-                            m_model(model),
-                            m_pnt(pnt),
-                            m_srcEng(srcEng),
-                            m_srcLogEng(srcLogEng),
-                            m_srcTime(srcTime),
-                            m_obsEng(obsEng),
-                            m_phi0(phi0),
-                            m_sin_fact_model(sin_fact_model),
-                            m_cos_fact_model(cos_fact_model),
-                            m_sin_fact_inst(sin_fact_inst),
-                            m_cos_fact_inst(cos_fact_inst) { }
+    cta_psf_diffuse_kern_phi(const GModelSpatial* model,
+                             const GEnergy&       srcEng,
+                             const GTime&         srcTime,
+                             const GMatrix&       rot,
+                             const double&        sin_delta,
+                             const double&        cos_delta) :
+                             m_model(model),
+                             m_srcEng(srcEng),
+                             m_srcTime(srcTime),
+                             m_rot(rot),
+                             m_sin_delta(sin_delta),
+                             m_cos_delta(cos_delta) { }
     double eval(const double& phi);
 protected:
-    const GCTAResponseIrf*     m_rsp;            //!< Response
-    const GModelSpatialRadial* m_model;          //!< Radial model
-    const GCTAPointing&        m_pnt;            //!< CTA pointing
-    const GEnergy&             m_srcEng;         //!< True photon energy
-    const double&              m_srcLogEng;      //!< Log10(true photon energy)
-    const GTime&               m_srcTime;        //!< True photon arrival time
-    const GEnergy&             m_obsEng;         //!< Measured event energy
-    const double&              m_phi0;           //!< phi0
-    const double&              m_sin_fact_model; //!< sin(delta)*sin(zeta)
-    const double&              m_cos_fact_model; //!< cos(delta)*cos(zeta)
-    const double&              m_sin_fact_inst;  //!< sin(delta)*sin(obsOffset)
-    const double&              m_cos_fact_inst;  //!< cos(delta)*cos(obsOffset)
+    const GModelSpatial* m_model;     //!< Spatial model
+    const GEnergy&       m_srcEng;    //!< True photon energy
+    const GTime&         m_srcTime;   //!< True photon arrival time
+    const GMatrix&       m_rot;       //!< Rotation matrix
+    const double&        m_sin_delta; //!< sin(delta)
+    const double&        m_cos_delta; //!< cos(delta)
 };
 #endif /* GCTARESPONSE_HELPERS_HPP */