Merge branch '1304-improve-ellipse-fitting' into integration

Conflicts: ChangeLog inst/cta/src/GCTAResponse_helpers.cpp inst/cta/src/GCTASupport.hpp

Merge branch '1304-improve-ellipse-fitting' into integration
Conflicts: ChangeLog inst/cta/src/GCTAResponse_helpers.cpp inst/cta/src/GCTASupport.hpp
Jürgen Knödlseder
2 parents 668e2cd5 0c5d33d6
Showing 22 changed files with 3259 additions and 1300 deletions
ChangeLog
NEWS
include/GModelSpatialElliptical.hpp
include/GModelSpatialEllipticalDisk.hpp
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
inst/cta/test/benchmark_ml_fitting.py
inst/cta/test/test_irf_elliptical.py
inst/cta/test/test_irf_radial.py
inst/cta/test/test_likelihood_profile.py
inst/cta/test/test_npred_elliptical.py
inst/cta/test/test_psf_radial.py
inst/lat/test/test_LAT.cpp
inst/mwl/test/test_MWL.cpp
pyext/GModelSpatialElliptical.i
pyext/GModelSpatialEllipticalDisk.i
src/model/GModelSpatialElliptical.cpp
-2014-10-26  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-10-30  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 0.9.0 released
         ========================
@@ -19,8 +19,8 @@
         Read/write altitude and azimuth attributes of CTA data
         Refactor CTA response class, add GCTAResponseCube class
         Improve CTA response computation for extended models
-        Implement CTA cube-style analysis
-        Optimize CTA binned and cube-style analysis
+        Implement CTA stacked analysis
+        Optimize CTA binned and stacked analysis
         Load events only when needed into GCTAObservation
         Implement full automatic typecasting in Python module
         Fill DETX and DETY values in GCTAResponseIrf::mc 
@@ -236,7 +236,7 @@ The type of the GCTAModelCubeBackground class has been set to GCTACubeBackground
 the type of the GCTAModelIrfBackground class to GCTAIrfBackground. XML files need
 to be adapted accordingly.
-For cube-style analysis, the GCTAExposure and GCTAMeanPsf classes have been added.
+For stacked analysis, the GCTAExposure and GCTAMeanPsf classes have been added.
 GCTAExposure represents the summed exposure for a list of CTA observations
 (computed as the sum of effective area times livetime), GCTAMeanPsf represents
 the average PSF, computed by weighing the PSF of each observation by the respective
@@ -257,10 +257,10 @@ GCTAResponseIrf derives. The class GCTAResponseCube, which also derives from
 GCTAResponse, has been added to handle cube-style analysis. GCTAResponseCube
 holds an exposure cube (GCTAExposure) and a PSF cube (GCTAMeanPsf) to describe
 the response. GCTAObservation has been modified to allocate either IRF-style or
-cube-style response, dependent on the information found in the XML file.
+stacked response, dependent on the information found in the XML file.
 The abstract GCTASourceCube class has been added to handle the pre-computation
-cache for cube-style analysis. A first implementation of this cache is provided
+cache for stacked analysis. A first implementation of this cache is provided
 by the GCTASourceCubePointSource class for point source analysis. The cache has
 been added to GCTAResponseCube who takes care of initializing and using the cache.
@@ -95,6 +95,10 @@ public:
     void    dec(const double& dec);
     double  posangle(void) const;
     void    posangle(const double& posangle);
+    double  semiminor(void) const;
+    double  semimajor(void) const;
+    void    semiminor(const double& semiminor);
+    void    semimajor(const double& semimajor);
     GSkyDir dir(void) const;
     void    dir(const GSkyDir& dir);
@@ -105,9 +109,11 @@ protected:
     void free_members(void);
     // Proteced members
-    GModelPar m_ra;       //!< Right Ascension (deg)
-    GModelPar m_dec;      //!< Declination (deg)
-    GModelPar m_posangle; //!< Position angle from North, counterclockwise (deg)
+    GModelPar m_ra;        //!< Right Ascension (deg)
+    GModelPar m_dec;       //!< Declination (deg)
+    GModelPar m_posangle;  //!< Position angle from North, counterclockwise (deg)
+    GModelPar m_semiminor; //!< Semi-minor axis of ellipse (deg)
+    GModelPar m_semimajor; //!< Semi-major axis of ellipse (deg)
 };
@@ -215,6 +221,64 @@ void GModelSpatialElliptical::posangle(const double&amp; posangle)
 /***********************************************************************//**
+ * @brief Return semi-minor axis of ellipse
+ *
+ * @return Semi-minor axis of ellipse (degrees).
+ *
+ * Returns the semi-minor axis of the ellipse in degrees.
+ ***************************************************************************/
+inline
+double GModelSpatialElliptical::semiminor(void) const
+{
+    return (m_semiminor.value());
+}
+
+
+/***********************************************************************//**
+ * @brief Set semi-minor axis of ellipse
+ *
+ * @param[in] semiminor Semi-minor axis of ellipse (degrees)
+ *
+ * Sets the semi-minor axis of the ellipse in degrees.
+ ***************************************************************************/
+inline
+void GModelSpatialElliptical::semiminor(const double& semiminor)
+{
+    m_semiminor.value(semiminor);
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return semi-major axis of ellipse
+ *
+ * @return Semi-major axis of ellipse (degrees).
+ *
+ * Returns the semi-major axis of the ellipse in degrees.
+ ***************************************************************************/
+inline
+double GModelSpatialElliptical::semimajor(void) const
+{
+    return (m_semimajor.value());
+}
+
+
+/***********************************************************************//**
+ * @brief Set semi-major axis of ellipse
+ *
+ * @param[in] semimajor Semi-major axis of ellipse (degrees)
+ *
+ * Sets the semi-major axis of the ellipse in degrees.
+ ***************************************************************************/
+inline
+void GModelSpatialElliptical::semimajor(const double& semimajor)
+{
+    m_semimajor.value(semimajor);
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Return normalization of elliptical source
  *
  * @return Normalization.
@@ -81,11 +81,6 @@ public:
     virtual void                         write(GXmlElement& xml) const;
     virtual std::string                  print(const GChatter& chatter = NORMAL) const;
-    // Other methods
-    double semiminor(void) const;
-    double semimajor(void) const;
-    void   semiminor(const double& semiminor);
-    void   semimajor(const double& semimajor);
 protected:
     // Protected methods
@@ -94,10 +89,6 @@ protected:
     void free_members(void);
     void update(void) const;
-    // Protected members
-    GModelPar      m_semiminor;        //!< semi-minor axis of ellipse (degrees)
-    GModelPar      m_semimajor;        //!< semi-major axis of ellipse (degrees)
-
     // Cached members used for pre-computations
     mutable double m_last_semiminor;   //!< Last semi-minor axis
     mutable double m_last_semimajor;   //!< Last semi-major axis
@@ -132,62 +123,4 @@ std::string GModelSpatialEllipticalDisk::type(void) const
     return "EllipticalDisk";
 }
-
-/***********************************************************************//**
- * @brief Return semi-minor axis of ellipse
- *
- * @return Semi-minor axis of ellipse (degrees).
- *
- * Returns the semi-minor axis of the ellipse in degrees.
- ***************************************************************************/
-inline
-double GModelSpatialEllipticalDisk::semiminor(void) const
-{
-    return (m_semiminor.value());
-}
-
-
-/***********************************************************************//**
- * @brief Set semi-minor axis of ellipse
- *
- * @param[in] semiminor Semi-minor axis of ellipse (degrees)
- *
- * Sets the semi-minor axis of the ellipse in degrees.
- ***************************************************************************/
-inline
-void GModelSpatialEllipticalDisk::semiminor(const double& semiminor)
-{
-    m_semiminor.value(semiminor);
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Return semi-major axis of ellipse
- *
- * @return Semi-major axis of ellipse (degrees).
- *
- * Returns the semi-major axis of the ellipse in degrees.
- ***************************************************************************/
-inline
-double GModelSpatialEllipticalDisk::semimajor(void) const
-{
-    return (m_semimajor.value());
-}
-
-
-/***********************************************************************//**
- * @brief Set semi-major axis of ellipse
- *
- * @param[in] semimajor Semi-major axis of ellipse (degrees)
- *
- * Sets the semi-major axis of the ellipse in degrees.
- ***************************************************************************/
-inline
-void GModelSpatialEllipticalDisk::semimajor(const double& semimajor)
-{
-    m_semimajor.value(semimajor);
-    return;
-}
-
 #endif /* GMODELSPATIALELLIPTICALDISK_HPP */
@@ -112,13 +112,13 @@ 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;
     double psf_elliptical(const GModelSpatialElliptical* model,
-                          const double&                  zeta,
-                          const double&                  omega,
+                          const double&                  rho_obs,
+                          const double&                  posangle_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,6 +64,7 @@
 /* __ Macros _____________________________________________________________ */
 /* __ Coding definitions _________________________________________________ */
+#define G_RADIAL_PSF_BASED  //!< Use Psf-based integration for radial model
 /* __ Debug definitions __________________________________________________ */
@@ -236,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,
@@ -307,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,
@@ -347,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,
@@ -425,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,
@@ -462,19 +463,19 @@ 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);
         // Continue only if exposure is positive
         if (irf > 0.0) {
-            irf /= obs.livetime(); //!< Revover effective area from exposure
-            irf *= psf_radial(model, zeta, obsDir, obsEng, obsTime);
+            irf /= obs.livetime(); //!< Recover effective area from exposure
+            irf *= psf_radial(model, rho_obs, obsDir, obsEng, obsTime);
             irf *= obs.deadc(obsTime);
         }
@@ -496,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,
@@ -533,20 +533,20 @@ double GCTAResponseCube::irf_elliptical(const GEvent&amp;       event,
     // Compute angle between model centre and measured photon direction and
     // position angle (radians)
-    double zeta  = model->dir().dist(obsDir);
-    double omega = model->dir().posang(obsDir);
+    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
-    if (zeta <= model->theta_max()+psf().delta_max()) {
+    if (rho_obs <= model->theta_max()+psf().delta_max()) {
         // Get exposure
         irf = exposure()(obsDir, obsEng);
         // Continue only if exposure is positive
         if (irf > 0.0) {
-            irf /= obs.livetime(); //!< Revover effective area from exposure
-            irf *= psf_elliptical(model, zeta, omega, obsDir, obsEng, obsTime);
+            irf /= obs.livetime(); //!< Recover effective area from exposure
+            irf *= psf_elliptical(model, rho_obs, posangle_obs, obsDir, obsEng, obsTime);
             irf *= obs.deadc(obsTime);
         }
@@ -568,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,
@@ -651,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
@@ -668,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)
@@ -737,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
@@ -773,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
 {
@@ -818,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                             =
  =                                                                         =
  ==========================================================================*/
@@ -857,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)
 {
@@ -893,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
 {
@@ -921,142 +912,362 @@ 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] zeta Angle between model centre and measured photon direction
- *                 (radians).
+ * @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.
  *
- * 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_{\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${\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}(\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&              zeta,
+                                    const double&              delta_mod,
                                     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
+    // Set number of iterations for Romberg integration.
+    // 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;
+    static const int iter_delta = 5;
+    static const int iter_phi   = 6;
     // 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 maximum Psf radius (radians)
+    double psf_max = psf().delta_max();
-    // Get offset angle integration interval in radians
-    double delta_min = 0.0;
-    double delta_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,
-                                        sin_zeta, cos_zeta,
-                                        eps_phi, order_phi);
-
-    // Integrate over PSF delta angle
+    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.eps(eps_delta);
-    value = integral.romberg(delta_min, delta_max, order_delta);
+    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);
+    }
-    // Return PSF
+    // 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
+ *                    (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_{\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$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,
+                                    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/1299
+    static const int iter_rho = 5;
+    static const int iter_phi = 5;
+
+    // Initialise value
+    double irf = 0.0;
+
+    // Get maximum PSF radius (radians)
+    double delta_max = psf().delta_max();
+
+    // 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();
+    if (rho_max > src_max) {
+        rho_max = src_max;
+    }
+
+    // Perform zenith angle integration if interval is valid
+    if (rho_max > rho_min) {
+
+        // 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_rho integrand(this,
+                                          model,
+                                          obsDir,
+                                          srcEng,
+                                          srcTime,
+                                          rho_obs,
+                                          delta_max,
+                                          iter_phi);
+
+        // Integrate over model's zenith angle
+        GIntegral integral(&integrand);
+        integral.fixed_iter(iter_rho);
+
+        // Setup integration boundaries
+        std::vector<double> bounds;
+        bounds.push_back(rho_min);
+        bounds.push_back(rho_max);
+
+        // If the integration range includes a transition between full
+        // containment of model within Psf 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 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);
+
+        // Compile option: Check for NaN/Inf
+        #if defined(G_NAN_CHECK)
+        if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
+            std::cout << "*** ERROR: GCTAResponseCube::psf_radial:";
+            std::cout << " NaN/Inf encountered";
+            std::cout << " (irf=" << irf;
+            std::cout << ", rho_min=" << rho_min;
+            std::cout << ", rho_max=" << rho_max << ")";
+            std::cout << std::endl;
+        }
+        #endif
+
+    } // endif: integration interval is valid
+
+    // Return PSF
+    return irf;
+}
+#endif
 /***********************************************************************//**
- * @brief Integrate PSF over elliptical model
+ * @brief Integrate Psf over elliptical model
  *
  * @param[in] model Elliptical model.
- * @param[in] zeta Angle between model centre and measured photon direction
- *                 (radians).
- * @param[in] omega Position angle of measured photon direction with
- *                  respect to model centre (radians)
+ * @param[in] rho_obs Angle between model centre and measured photon direction
+ *                    (radians).
+ * @param[in] posangle_obs Position angle of measured photon direction with
+ *                         respect to model centre (radians).
  * @param[in] obsDir Observed event direction.
  * @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 
+ *    \int_{\omega} 
+ *    S_{\rm p}(\rho, \omega | E, t) \times
+ *    {\rm Psf}(\rho, \omega) d\omega d\rho
  * \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.
+ * where
+ * \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&                  zeta,
-                                        const double&                  omega,
+                                        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;
+    // Set number of iterations for Romberg integration.
+    // These values have been determined after careful testing, see
+    // https://cta-redmine.irap.omp.eu/issues/1299
+    static const int iter_rho = 5;
+    static const int iter_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);
+    double irf = 0.0;
-    // Get offset angle integration interval in radians
-    double delta_min = 0.0;
+    // Get maximum PSF radius (radians)
     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_zeta, cos_zeta, omega,
-                                            eps_phi, order_phi);
+    // Get the ellipse boundary (radians). Note that these are NOT the
+    // parameters of the ellipse but of a boundary ellipse that is used
+    // for computing the relevant omega angle intervals for a given angle
+    // rho. The boundary ellipse takes care of the possibility that the
+    // semiminor axis is larger than the semimajor axis
+    double semimajor;    // Will be the larger axis
+    double semiminor;    // Will be the smaller axis
+    double posangle;     // Will be the corrected position angle
+    double aspect_ratio; // Ratio between smaller/larger axis of model
+    if (model->semimajor() >= model->semiminor()) {
+        aspect_ratio = (model->semimajor() > 0.0) ?
+                        model->semiminor() / model->semimajor() : 0.0;
+        posangle     = model->posangle() * gammalib::deg2rad;
+    }
+    else {
+        aspect_ratio = (model->semiminor() > 0.0) ?
+                        model->semimajor() / model->semiminor() : 0.0;
+        posangle     = model->posangle() * gammalib::deg2rad + gammalib::pihalf;
+    }
+    semimajor = model->theta_max();
+    semiminor = semimajor * aspect_ratio;
+
+    // Set zenith angle integration range for elliptical model
+    double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
+    double rho_max = rho_obs + delta_max;
+    if (rho_max > semimajor) {
+        rho_max = semimajor;
+    }
-    // Integrate over PSF delta angle
-    GIntegral integral(&integrand);
-    integral.eps(eps_delta);
-    value = integral.romberg(delta_min, delta_max, order_delta);
+    // Perform zenith angle integration if interval is valid
+    if (rho_max > rho_min) {
+
+        // 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_rho integrand(this,
+                                              model,
+                                              semimajor,
+                                              semiminor,
+                                              posangle,
+                                              obsDir,
+                                              srcEng,
+                                              srcTime,
+                                              rho_obs,
+                                              posangle_obs,
+                                              delta_max,
+                                              iter_phi);
+
+        // Integrate over model's zenith angle
+        GIntegral integral(&integrand);
+        integral.fixed_iter(iter_rho);
+
+        // Setup integration boundaries
+        std::vector<double> bounds;
+        bounds.push_back(rho_min);
+        bounds.push_back(rho_max);
+
+        // If the integration range includes the semiminor boundary, then
+        // add an integration boundary at that location
+        if (semiminor > rho_min && semiminor < rho_max) {
+            bounds.push_back(semiminor);
+        }
+
+        // Integrate kernel
+        irf = integral.romberg(bounds, iter_rho);
+
+        // Compile option: Check for NaN/Inf
+        #if defined(G_NAN_CHECK)
+        if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
+            std::cout << "*** ERROR: GCTAResponseCube::psf_elliptical:";
+            std::cout << " NaN/Inf encountered";
+            std::cout << " (irf=" << irf;
+            std::cout << ", rho_min=" << rho_min;
+            std::cout << ", rho_max=" << rho_max << ")";
+            std::cout << std::endl;
+        }
+        #endif
+
+    } // endif: integration interval is valid
     // Return PSF
-    return value;
+    return irf;
 }
@@ -1267,7 +1267,7 @@ std::string GCTAResponseIrf::print(const GChatter&amp; chatter) const
 /*==========================================================================
  =                                                                         =
- =              Model type dependent CTA response methods                  =
+ =                Model type dependent CTA response methods                =
  =                                                                         =
  ==========================================================================*/
@@ -1281,192 +1281,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
- *
- * \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
+ * 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
@@ -1596,8 +1426,8 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp;       event,
         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
+        // containment of model within Psf and partial containment, then
+        // add a boundary at this location
         double transition_point = delta_max - zeta;
         if (transition_point > rho_min && transition_point < rho_max) {
             bounds.push_back(transition_point);
@@ -1608,7 +1438,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);
             }
@@ -1654,11 +1484,10 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp;       event,
     // Return IRF value
     return irf;
 }
-#endif
 /***********************************************************************//**
- * @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.
@@ -1675,17 +1504,16 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp;       event,
  * \f[
  *    \int_{\rho_{\rm min}}^{\rho_{\rm max}}
  *    \sin \rho \times
- *    \int_{\omega_{\rm min}}^{\omega_{\rm max}} 
+ *    \int_{\omega} 
  *    S_{\rm p}(\rho, \omega | E, t) \, IRF(\rho, \omega) d\omega d\rho
  * \f]
  *
  * where
- * - \f$S_{\rm p}(\rho, \omega | 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, \omega | E, t)\f$ is the elliptical model,
+ * \f$IRF(\rho, \omega)\f$ is the instrument response 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.
  *
  * The source model centre is located at \f$\vec{m}\f$, and a spherical
  * coordinate system is defined around this location with \f$(\rho,\omega)\f$
@@ -1731,26 +1559,16 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp;       event,
     // Determine angular distance between observed photon direction and model
     // centre and position angle of observed photon direction seen from the
     // model centre [radians]
-    double zeta     = centre.dist(obsDir);
-    double obsOmega = centre.posang(obsDir);
-
-    // Determine angular distance between measured photon direction and
-    // pointing direction [radians]
-    double eta = pnt.dir().dist(obsDir);
+    double rho_obs      = centre.dist(obsDir);
+    double posangle_obs = centre.posang(obsDir);
     // Determine angular distance between model centre and pointing direction
     // [radians]
-    double lambda = centre.dist(pnt.dir());
+    double rho_pnt      = centre.dist(pnt.dir());
+    double posangle_pnt = centre.posang(pnt.dir());
     // Compute azimuth angle of pointing in model coordinate system [radians]
-    // This azimuth angle is comprised in the interval [0,pi], and defines
-    // the zero point of the model coordinate system.
-    double omega0 = 0.0;
-    double denom  = std::sin(lambda) * std::sin(zeta);
-    if (denom != 0.0) {
-        double arg = (std::cos(eta) - std::cos(lambda) * std::cos(zeta))/denom;
-        omega0     = gammalib::acos(arg);
-    }
+    double omega_pnt = posangle_pnt - posangle_obs;
     // Get log10(E/TeV) of true photon energy
     double srcLogEng = srcEng.log10TeV();
@@ -1759,18 +1577,37 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp;       event,
     // angle (eta) as the true theta angle between the source and the pointing
     // directions. As we only use the angle to determine the maximum PSF size,
     // this should be sufficient.
-    double theta     = eta;
+    double theta     = pnt.dir().dist(obsDir);
     double phi       = 0.0; //TODO: Implement IRF Phi dependence
     double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
-    // Get maximum source model radius [radians]
-    double src_max = model->theta_max();
+    // Get the ellipse boundary (radians). Note that these are NOT the
+    // parameters of the ellipse but of a boundary ellipse that is used
+    // for computing the relevant omega angle intervals for a given angle
+    // rho. The boundary ellipse takes care of the possibility that the
+    // semiminor axis is larger than the semimajor axis
+    double semimajor;    // Will be the larger axis
+    double semiminor;    // Will be the smaller axis
+    double posangle;     // Will be the corrected position angle
+    double aspect_ratio; // Ratio between smaller/larger axis of model
+    if (model->semimajor() >= model->semiminor()) {
+        aspect_ratio = (model->semimajor() > 0.0) ?
+                        model->semiminor() / model->semimajor() : 0.0;
+        posangle     = model->posangle() * gammalib::deg2rad;
+    }
+    else {
+        aspect_ratio = (model->semiminor() > 0.0) ?
+                        model->semimajor() / model->semiminor() : 0.0;
+        posangle     = model->posangle() * gammalib::deg2rad + gammalib::pihalf;
+    }
+    semimajor = model->theta_max();
+    semiminor = semimajor * aspect_ratio;
     // Set zenith angle integration range for elliptical model
-    double rho_min = (zeta > delta_max) ? zeta - delta_max : 0.0;
-    double rho_max = zeta + delta_max;
-    if (rho_max > src_max) {
-        rho_max = src_max;
+    double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
+    double rho_max = rho_obs + delta_max;
+    if (rho_max > semimajor) {
+        rho_max = semimajor;
     }
     // Initialise IRF value
@@ -1782,16 +1619,19 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp;       event,
         // Setup integration kernel
         cta_irf_elliptical_kern_rho integrand(*this,
                                               *model,
+                                              semimajor,
+                                              semiminor,
+                                              posangle,
                                               zenith,
                                               azimuth,
                                               srcEng,
                                               srcTime,
                                               srcLogEng,
                                               obsEng,
-                                              zeta,
-                                              lambda,
-                                              obsOmega,
-                                              omega0,
+                                              rho_obs,
+                                              posangle_obs,
+                                              rho_pnt,
+                                              omega_pnt,
                                               delta_max,
                                               iter_phi);
@@ -1804,23 +1644,10 @@ 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 - zeta;
-        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
-        const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(model);
-        if (shell != NULL) {
-            double shell_radius = shell->radius();
-            if (shell_radius > rho_min && shell_radius < rho_max) {
-                bounds.push_back(shell_radius);
-            }
+        // If the integration range includes the semiminor boundary, then
+        // add an integration boundary at that location
+        if (semiminor > rho_min && semiminor < rho_max) {
+            bounds.push_back(semiminor);
         }
         // Integrate kernel
@@ -1833,8 +1660,7 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp;       event,
             std::cout << " NaN/Inf encountered";
             std::cout << " (irf=" << irf;
             std::cout << ", rho_min=" << rho_min;
-            std::cout << ", rho_max=" << rho_max;
-            std::cout << ", omega0=" << omega0 << ")";
+            std::cout << ", rho_max=" << rho_max << ")";
             std::cout << std::endl;
         }
         #endif
@@ -2188,19 +2014,27 @@ double GCTAResponseIrf::npred_radial(const GSource&amp; source,
         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
+        // containment of model within ROI and partial containment, then
+        // add a boundary at this location
         double transition_point = roi_psf_radius - roi_model_distance;
         if (transition_point > rho_min && transition_point < rho_max) {
             bounds.push_back(transition_point);
         }
+        // If the integration range includes a transition between full
+        // containment of ROI within model and partial containment, then
+        // add a boundary at this location
+        transition_point = roi_psf_radius + roi_model_distance;
+        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
         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);
             }
@@ -2327,23 +2161,44 @@ double GCTAResponseIrf::npred_elliptical(const GSource&amp; source,
     // to make sure we encompass the entire PSF.
     double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
-    // Extract ROI radius (radians)
-    double roi_radius = roi.radius() * gammalib::deg2rad;
-
-    // Compute distance between ROI and model centre (radians)
-    double roi_model_distance = roi.centre().dir().dist(centre);
-
-    // Compute the ROI radius plus maximum PSF radius (radians). Any photon
+    // Extract ROI radius plus maximum PSF radius (radians). Any photon
     // coming from beyond this radius will not make it in the dataspace and
     // thus can be neglected.
-    double roi_psf_radius = roi_radius + psf_max_radius;
+    double radius_roi = roi.radius() * gammalib::deg2rad + psf_max_radius;
+
+    // Compute distance between ROI and model centre (radians)
+    double rho_roi = roi.centre().dir().dist(centre);
+
+    // Get the ellipse boundary (radians). Note that these are NOT the
+    // parameters of the ellipse but of a boundary ellipse that is used
+    // for computing the relevant omega angle intervals for a given angle
+    // rho. The boundary ellipse takes care of the possibility that the
+    // semiminor axis is larger than the semimajor axis
+    double semimajor;    // Will be the larger axis
+    double semiminor;    // Will be the smaller axis
+    double posangle;     // Will be the corrected position angle
+    double aspect_ratio; // Ratio between smaller/larger axis of model
+    if (model->semimajor() >= model->semiminor()) {
+        aspect_ratio = (model->semimajor() > 0.0) ?
+                        model->semiminor() / model->semimajor() : 0.0;
+        posangle     = model->posangle() * gammalib::deg2rad;
+    }
+    else {
+        aspect_ratio = (model->semiminor() > 0.0) ?
+                        model->semimajor() / model->semiminor() : 0.0;
+        posangle     = model->posangle() * gammalib::deg2rad + gammalib::pihalf;
+    }
+    semimajor = model->theta_max();
+    semiminor = semimajor * aspect_ratio;
     // Set offset angle integration range. We take here the ROI+PSF into
     // account to make no integrations beyond the point where the
     // contribution drops to zero.
-    double rho_min = (roi_model_distance > roi_psf_radius)
-                     ? roi_model_distance - roi_psf_radius: 0.0;
-    double rho_max = model->theta_max();
+    double rho_min = (rho_roi > radius_roi) ? rho_roi - radius_roi: 0.0;
+    double rho_max = rho_roi + radius_roi;
+    if (rho_max > semimajor) {
+        rho_max = semimajor;
+    }
     // Perform offset angle integration only if interval is valid
     if (rho_max > rho_min) {
@@ -2358,18 +2213,21 @@ double GCTAResponseIrf::npred_elliptical(const GSource&amp; source,
         // Compute position angle of ROI centre with respect to model
         // centre (radians)
-        double omega0 = centre.posang(roi.centre().dir());
+        double posangle_roi = centre.posang(roi.centre().dir());
         // Setup integration kernel
         cta_npred_elliptical_kern_rho integrand(*this,
                                                 *model,
-                                                source.energy(),
-                                                source.time(),
+                                                semimajor,
+                                                semiminor,
+                                                posangle,
+                                                srcEng,
+                                                srcTime,
                                                 cta,
                                                 rot,
-                                                roi_model_distance,
-                                                roi_psf_radius,
-                                                omega0,
+                                                rho_roi,
+                                                posangle_roi,
+                                                radius_roi,
                                                 iter_phi);
         // Integrate over model's zenith angle
@@ -2381,23 +2239,10 @@ double GCTAResponseIrf::npred_elliptical(const GSource&amp; source,
         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 = roi_psf_radius - roi_model_distance;
-        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
-        const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(model);
-        if (shell != NULL) {
-            double shell_radius = shell->radius();
-            if (shell_radius > rho_min && shell_radius < rho_max) {
-                bounds.push_back(shell_radius);
-            }
+        // If the integration range includes the semiminor boundary, then
+        // add an integration boundary at that location
+        if (semiminor > rho_min && semiminor < rho_max) {
+            bounds.push_back(semiminor);
         }
         // Integrate kernel
@@ -46,12 +46,111 @@
 /* __ Debug definitions __________________________________________________ */
 //#define G_DEBUG_INTEGRAL                             //!< Debug integration
+#define G_DEBUG_MODEL_ZERO           //!< Debug check for zero model values
 /* __ Constants __________________________________________________________ */
 /*==========================================================================
  =                                                                         =
+ =                              Helper functions                           =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Limit omega interval
+ *
+ * @param[in] min Interval minimum (radians).
+ * @param[in] max Interval maximum (radians).
+ * @param[in] domega Half length of interval (radians).
+ * @return Vector of intervals.
+ *
+ * Limits an omega interval [@p min,@p max] to the interval specified by
+ * [-@p domega,@p domega]. This may lead to a split of [@p min,@p max] in
+ * several intervals. The method thus returns a vector of intervals that
+ * overlap with [-@p domega,@p domega]. If there is no overlap with the
+ * interval, the method returns an empty vector.
+ *
+ * The method takes care of wrap arounds. It is assumed that on input
+ * [@p min,@p max] is contained within [-2pi,+2pi].
+ ***************************************************************************/
+cta_omega_intervals gammalib::limit_omega(const double& min,
+                                          const double& max,
+                                          const double& domega)
+{
+    // Allocate intervals
+    cta_omega_intervals intervals;
+
+    // Continue only if domega is smaller than pi
+    if (domega < gammalib::pi) {
+
+        // Set limiting intervals. To take care of a possible wrap around in
+        // omega we consider also intervals that are shifted by +/-2pi
+        double omega_min       = -domega;
+        double omega_max       = +domega;
+        double omega_min_plus  = omega_min + gammalib::twopi;
+        double omega_max_plus  = omega_max + gammalib::twopi;
+        double omega_min_minus = omega_min - gammalib::twopi;
+        double omega_max_minus = omega_max - gammalib::twopi;
+
+        // If the [min,max] interval overlaps with the unshifted
+        // [-domega,domega] interval then constrain the interval
+        if (max > omega_min && min < omega_max) {
+            double interval_min = min;
+            double interval_max = max;
+            if (interval_min < omega_min) {
+                interval_min = omega_min;
+            }
+            if (interval_max > omega_max) {
+                interval_max = omega_max;
+            }
+            intervals.push_back(std::make_pair(interval_min,interval_max));
+        }
+
+        // If the [min,max] interval overlaps with the [-domega,domega]
+        // interval shifted by +2pi then constrain the interval using the
+        // shifted interval
+        if (max > omega_min_plus && min < omega_max_plus) {
+            double interval_min = min;
+            double interval_max = max;
+            if (interval_min < omega_min_plus) {
+                interval_min = omega_min_plus;
+            }
+            if (interval_max > omega_max_plus) {
+                interval_max = omega_max_plus;
+            }
+            intervals.push_back(std::make_pair(interval_min,interval_max));
+        }
+ 
+        // If the [min,max] interval overlaps with the [-domega,domega]
+        // interval shifted by -2pi then constrain the interval using the
+        // shifted interval
+        if (max > omega_min_minus && min < omega_max_minus) {
+            double interval_min = min;
+            double interval_max = max;
+            if (interval_min < omega_min_minus) {
+                interval_min = omega_min_minus;
+            }
+            if (interval_max > omega_max_minus) {
+                interval_max = omega_max_minus;
+            }
+            intervals.push_back(std::make_pair(interval_min,interval_max));
+        }
+
+    } // endif: interval was not the full circle
+
+    // ... otherwise append the provided interval
+    else {
+        intervals.push_back(std::make_pair(min,max));
+    }
+
+    // Return intervals
+    return intervals;
+}
+
+
+/*==========================================================================
+ =                                                                         =
  =              Helper class methods for response computation              =
  =                                                                         =
  ==========================================================================*/
@@ -184,8 +283,26 @@ 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)
+            if (model <= 0.0) {
+                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 << " radians" << std::endl;
+            }
+            #endif
             // Continue only if model is positive
             if (model > 0.0) {
@@ -336,7 +453,7 @@ double cta_irf_radial_kern_omega::eval(const double&amp; omega)
  * \f]
  * 
  * The azimuth angle integration range 
- * \f$[\omega_{\rm min}, \omega_{\rm max}\f$
+ * \f$[\omega_{\rm min}, \omega_{\rm max}]\f$
  * is limited to an arc around the vector connecting the model centre to
  * the ROI centre. This limitation assures that the integration converges
  * properly.
@@ -364,8 +481,26 @@ 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)
+            if (model <= 0.0) {
+                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 << " radians" << std::endl;
+            }
+            #endif
             // Continue only if we have a positive model value
             if (model > 0.0) {
@@ -470,7 +605,8 @@ double cta_npred_radial_kern_omega::eval(const double&amp; omega)
 /***********************************************************************//**
  * @brief Kernel for elliptical model integration over model's zenith angle
  *
- * @param[in] rho Zenith angle with respect to model centre [radians].
+ * @param[in] rho Radial distance from model centre [radians].
+ * @return Radial IRF integration kernel.
  *
  * Computes
  *
@@ -480,6 +616,20 @@ double cta_npred_radial_kern_omega::eval(const double&amp; omega)
  *                     S_{\rm p}(\rho, \omega | E, t) \, IRF(\rho, \omega)
  *                     d\omega
  * \f]
+ *
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the observed photon direction, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical source model
+ * for a given true photon energy and photon arrival time,
+ * \f$IRF(\rho, \omega)\f$ is the instrument response function.
+ *
+ * The method performs the required coordinate transformations from the
+ * model system, spanned by \f$(\rho, \omega)\f$, to the system needed for
+ * IRF computations. Furthermore, the method limits the integration range
+ * to area where the ellipse intersects the IRF.
  ***************************************************************************/
 double cta_irf_elliptical_kern_rho::eval(const double& rho)
 {
@@ -489,29 +639,34 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
     // Continue only if rho is positive
     if (rho > 0.0) {
-        // Compute half length of arc that lies within PSF validity circle
-        // (in radians)
+        // Compute half length of the arc (in radians) from a circle with
+        // radius rho that intersects with the point spread function, defined
+        // as a circle with maximum radius m_delta_max
         double domega = 0.5 * gammalib::cta_roi_arclength(rho,
-                                                          m_zeta,
-                                                          m_cos_zeta,
-                                                          m_sin_zeta,
+                                                          m_rho_obs,
+                                                          m_cos_rho_obs,
+                                                          m_sin_rho_obs,
                                                           m_delta_max,
                                                           m_cos_delta_max);
         // Continue only if arc length is positive
         if (domega > 0.0) {
-            // Compute omega integration range
-            double omega_min = -domega;
-            double omega_max = +domega;
-
             // Precompute cosine and sine terms for azimuthal integration
             double cos_rho = std::cos(rho);
             double sin_rho = std::sin(rho);
-            double cos_psf = cos_rho*m_cos_zeta;
-            double sin_psf = sin_rho*m_sin_zeta;
-            double cos_ph  = cos_rho*m_cos_lambda;
-            double sin_ph  = sin_rho*m_sin_lambda;
+            double cos_psf = cos_rho * m_cos_rho_obs;
+            double sin_psf = sin_rho * m_sin_rho_obs;
+            double cos_ph  = cos_rho * m_cos_rho_pnt;
+            double sin_ph  = sin_rho * m_sin_rho_pnt;
+
+            // Reduce rho by an infinite amount to avoid rounding errors
+            // at the boundary of a sharp edged model (e.g. an elliptical
+            // disk model)
+            double rho_kluge = rho - 1.0e-12;
+            if (rho_kluge < 0.0) {
+                rho_kluge = 0.0;
+            }
             // Setup integration kernel
             cta_irf_elliptical_kern_omega integrand(m_rsp,
@@ -522,18 +677,93 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
                                                     m_srcTime,
                                                     m_srcLogEng,
                                                     m_obsEng,
-                                                    m_obsOmega,
-                                                    m_omega0,
-                                                    rho,
+                                                    m_posangle_obs,
+                                                    m_omega_pnt,
+                                                    rho_kluge,
                                                     cos_psf,
                                                     sin_psf,
                                                     cos_ph,
                                                     sin_ph);
-            // Integrate over phi
+            // Setup integrator
             GIntegral integral(&integrand);
             integral.fixed_iter(m_iter);
-            irf = integral.romberg(omega_min, omega_max, m_iter) * sin_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) {
+
+                // Compute omega integration range
+                double omega_min = -domega;
+                double omega_max = +domega;
+
+                // Integrate over omega
+                irf = integral.romberg(omega_min, omega_max, m_iter) *
+                      sin_rho;
+
+            } // endif: circle comprised in ellipse
+
+            // ... otherwise there are arcs that intersect with the Psf circle
+            else {
+
+                // Compute half the arc length (in radians) of a circle of
+                // radius rho, centred on the model, that intersects with
+                // the ellipse boundary
+                double arg1 = 1.0 - (m_semiminor*m_semiminor) / (rho*rho);
+                double arg2 = 1.0 - (m_semiminor*m_semiminor) /
+                                    (m_semimajor*m_semimajor);
+                double omega_width = std::acos(std::sqrt(arg1/arg2));
+
+                // Continue only if the arclength is positive
+                if (omega_width > 0.0) {
+
+                    // Compute azimuth angle difference between ellipse
+                    // position angle and position angle of observed
+                    // photon in the model system. This angle will define
+                    // the reference point around which the circle arc. Make
+                    // sure that omega_0 is within [-pi,pi] thus that the omega
+                    // intervals are between [-2pi,2pi]
+                    double omega_0 = m_posangle - m_posangle_obs;
+                    if (omega_0 > gammalib::pi) {
+                        omega_0 -= gammalib::pi;
+                    }
+                    else if (omega_0 < -gammalib::pi) {
+                        omega_0 += gammalib::pi;
+                    }
+
+                    // Compute azimuth angle intervals
+                    double omega1_min = omega_0    - omega_width;
+                    double omega1_max = omega_0    + omega_width;
+                    double omega2_min = omega1_min + gammalib::pi;
+                    double omega2_max = omega1_max + gammalib::pi;
+
+                    // Limit intervals to the intersection of the ellipse with
+                    // the Psf circle. This may lead to a split of intervals,
+                    // and we gather all these intervals in a special interval
+                    // pair containers
+                    cta_omega_intervals intervals1 = 
+                        gammalib::limit_omega(omega1_min, omega1_max, domega);
+                    cta_omega_intervals intervals2 = 
+                        gammalib::limit_omega(omega2_min, omega2_max, domega);
+
+                    // Integrate over all intervals for omega1
+                    for (int i = 0; i < intervals1.size(); ++i) {
+                        double min = intervals1[i].first;
+                        double max = intervals1[i].second;
+                        irf       += integral.romberg(min, max, m_iter) * sin_rho;
+                    }
+
+                    // Integrate over all intervals for omega2
+                    for (int i = 0; i < intervals2.size(); ++i) {
+                        double min = intervals2[i].first;
+                        double max = intervals2[i].second;
+                        irf       += integral.romberg(min, max, m_iter) * sin_rho;
+                    }
+
+                } // endif: arc length was positive
+
+            } // endelse: circle was not comprised in ellipse
             // Compile option: Check for NaN/Inf
             #if defined(G_NAN_CHECK)
@@ -565,12 +795,18 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
  * Computes
  *
  * \f[
- *    S_{\rm p}(\rho, \omega | E, t) \, IRF(\rho, \omega)
+ *    S_{\rm p}(\omega | \rho, E, t) \, IRF(\omega | \rho)
  * \f]
  *
- * From the model coordinates \f$(\rho,\omega)\f$, the method computes the
- * angle between the true (\f$\vec{p}\f$) and observed (\f$\vec{p'}\f$) 
- * photon arrival direction using
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the observed photon direction, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ *
+ * From the coordinates \f$(\rho,\omega)\f$ in the model system, the method
+ * computes the angle between the true (\f$\vec{p}\f$) and observed
+ * (\f$\vec{p'}\f$) photon arrival direction using
  * 
  * \f[\delta = \arccos(\cos \rho \cos \zeta + 
  *                     \sin \rho \sin \zeta \cos \omega)\f]
@@ -598,18 +834,42 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
     // Initialise IRF value
     double irf = 0.0;
+    // Compute azimuth angle in model coordinate system (radians)
+    double omega_model = omega + m_posangle_obs;
+
     // Evaluate sky model
-    double model = m_model.eval(m_rho, omega + m_obsOmega, m_srcEng, m_srcTime);
+    double model = m_model.eval(m_rho, omega_model, m_srcEng, m_srcTime);
+
+    // Debug: test if model is non positive
+    #if defined(G_DEBUG_MODEL_ZERO)
+    if (model <= 0.0) {
+        double m_semiminor_rad = m_model.semiminor() * gammalib::deg2rad;
+        double m_semimajor_rad = m_model.semimajor() * gammalib::deg2rad;
+        double diff_angle      = omega_model - m_model.posangle() * gammalib::deg2rad;
+        double cosinus         = std::cos(diff_angle);
+        double sinus           = std::sin(diff_angle);
+        double arg1            = m_semiminor_rad * cosinus;
+        double arg2            = m_semimajor_rad * sinus;
+        double r_ellipse       = m_semiminor_rad * m_semimajor_rad /
+                                 std::sqrt(arg1*arg1 + arg2*arg2);
+        std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval";
+        std::cout << " zero model for (rho,omega)=(";
+        std::cout << m_rho*gammalib::rad2deg << ",";
+        std::cout << omega*gammalib::rad2deg << ")";
+        std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
+        std::cout << std::endl;
+    }
+    #endif
     // Continue only if model is positive
     if (model > 0.0) {
-        // Compute PSF offset angle [radians]
+        // Compute Psf offset angle [radians]
         double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
         // Compute true photon offset and azimuth angle in camera system
         // [radians]
-        double theta = std::acos(m_cos_ph + m_sin_ph * std::cos(m_omega0 - omega));
+        double theta = std::acos(m_cos_ph + m_sin_ph * std::cos(m_omega_pnt - omega));
         double phi   = 0.0; //TODO: Implement IRF Phi dependence
         // Evaluate IRF * model
@@ -659,10 +919,20 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
  *                     N_{\rm pred}(\rho,\omega) d\omega
  * \f]
  *
- * The azimuth angle integration range 
- * \f$[\omega_{\rm min}, \omega_{\rm max}\f$
- * is limited to an arc around the vector connecting the model centre to
- * the ROI centre. This limitation assures proper converges properly.
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the ROI centre, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical source model
+ * for a given true photon energy and photon arrival time,
+ * \f$N_{\rm pred}(\rho,\omega)\f$ is the data space integral over the
+ * response function.
+ *
+ * The method performs the required coordinate transformations from the
+ * model system, spanned by \f$(\rho, \omega)\f$, to the system needed for
+ * Npred computations. Furthermore, the method limits the integration range
+ * to area where the ellipse intersects the ROI.
  ***************************************************************************/
 double cta_npred_elliptical_kern_rho::eval(const double& rho)
 {
@@ -674,23 +944,27 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
         // Compute half length of arc that lies within ROI+PSF radius (radians)
         double domega = 0.5 * gammalib::cta_roi_arclength(rho,
-                                                          m_dist,
-                                                          m_cos_dist,
-                                                          m_sin_dist,
-                                                          m_radius,
-                                                          m_cos_radius);
+                                                          m_rho_roi,
+                                                          m_cos_rho_roi,
+                                                          m_sin_rho_roi,
+                                                          m_radius_roi,
+                                                          m_cos_radius_roi);
         // Continue only if arc length is positive
         if (domega > 0.0) {
-            // Compute phi integration range
-            double omega_min = m_omega0 - domega;
-            double omega_max = m_omega0 + domega;
-
-            // Compute sine and cosine of offset angle
+            // Compute sine and cosine terms for azimuthal integration
             double sin_rho = std::sin(rho);
             double cos_rho = std::cos(rho);
+            // Reduce rho by an infinite amount to avoid rounding errors
+            // at the boundary of a sharp edged model (e.g. an elliptical
+            // disk model)
+            double rho_kluge = rho - 1.0e-12;
+            if (rho_kluge < 0.0) {
+                rho_kluge = 0.0;
+            }
+
             // Setup phi integration kernel
             cta_npred_elliptical_kern_omega integrand(m_rsp,
                                                       m_model,
@@ -698,13 +972,90 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
                                                       m_srcTime,
                                                       m_obs,
                                                       m_rot,
+                                                      rho_kluge,
                                                       sin_rho,
-                                                      cos_rho);
+                                                      cos_rho,
+                                                      m_posangle_roi);
-            // Integrate over phi
+            // Setup integrator
             GIntegral integral(&integrand);
             integral.fixed_iter(m_iter);
-            npred = integral.romberg(omega_min, omega_max, m_iter) * sin_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) {
+
+                // Compute omega integration range
+                double omega_min = -domega;
+                double omega_max = +domega;
+
+                // Integrate over omega
+                npred = integral.romberg(omega_min, omega_max, m_iter) *
+                        sin_rho;
+
+            } // endif: circle comprised in ellipse
+
+            // ... otherwise there are arcs that intersect with the ROI
+            // circle
+            else {
+
+                // Compute half the arc length (in radians) of a circle of
+                // radius rho, centred on the model, that intersects with
+                // the ellipse boundary
+                double arg1 = 1.0 - (m_semiminor*m_semiminor) / (rho*rho);
+                double arg2 = 1.0 - (m_semiminor*m_semiminor) /
+                                    (m_semimajor*m_semimajor);
+                double omega_width = std::acos(std::sqrt(arg1/arg2));
+
+                // Continue only if the arclength is positive
+                if (omega_width > 0.0) {
+
+                    // Compute azimuth angle difference between ellipse
+                    // position angle and position angle of ROI in the model
+                    // system. This angle will define the reference point for
+                    // the circle arcs. Make sure that omega_0 is within [-pi,pi]
+                    // thus that the omega intervals are between [-2pi,2pi]
+                    double omega_0 = m_posangle - m_posangle_roi;
+                    if (omega_0 > gammalib::pi) {
+                        omega_0 -= gammalib::pi;
+                    }
+                    else if (omega_0 < -gammalib::pi) {
+                        omega_0 += gammalib::pi;
+                    }
+
+                    // Compute azimuth angle intervals
+                    double omega1_min = omega_0    - omega_width;
+                    double omega1_max = omega_0    + omega_width;
+                    double omega2_min = omega1_min + gammalib::pi;
+                    double omega2_max = omega1_max + gammalib::pi;
+
+                    // Limit intervals to the intersection of the ellipse with
+                    // the Psf circle. This may lead to a split of intervals,
+                    // and we gather all these intervals in a special interval
+                    // pair containers
+                    cta_omega_intervals intervals1 = 
+                        gammalib::limit_omega(omega1_min, omega1_max, domega);
+                    cta_omega_intervals intervals2 = 
+                        gammalib::limit_omega(omega2_min, omega2_max, domega);
+
+                    // Integrate over all intervals for omega1
+                    for (int i = 0; i < intervals1.size(); ++i) {
+                        double min = intervals1[i].first;
+                        double max = intervals1[i].second;
+                        npred     += integral.romberg(min, max, m_iter) * sin_rho;
+                    }
+
+                    // Integrate over all intervals for omega2
+                    for (int i = 0; i < intervals2.size(); ++i) {
+                        double min = intervals2[i].first;
+                        double max = intervals2[i].second;
+                        npred     += integral.romberg(min, max, m_iter) * sin_rho;
+                    }
+
+                } // endif: arc length was positive
+
+            } // endelse: circle was not comprised in ellipse
             // Debug: Check for NaN
             #if defined(G_NAN_CHECK)
@@ -713,7 +1064,6 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
                 std::cout << "(rho=" << rho << "):";
                 std::cout << " NaN/Inf encountered";
                 std::cout << " (npred=" << npred;
-                std::cout << ", omega=[" << omega_min << "," << omega_max << "]";
                 std::cout << ", sin_rho=" << sin_rho;
                 std::cout << ", cos_rho=" << cos_rho;
                 std::cout << ")" << std::endl;
@@ -737,34 +1087,68 @@ double cta_npred_elliptical_kern_rho::eval(const double&amp; rho)
  * Computes
  *
  * \f[
- *    S_{\rm p}(\rho,\omega | E, t) \, N_{\rm pred}(\rho,\omega)
+ *    S_{\rm p}(\omega | \rho, E, t) \, N_{\rm pred}(\omega | \rho)
  * \f]
+ *
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the centre of the Region of Interest (ROI), and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ *
+ * @todo Npred computation goes over sky coordinates. This can maybe be
+ *       optimized to reduce the number of coordinate transformations.
+ * @todo Check whether the Npred omega argument is the right one.
  ***************************************************************************/
 double cta_npred_elliptical_kern_omega::eval(const double& omega)
 {
     // Initialise Npred value
     double npred = 0.0;
-    // Compute sky direction vector in native coordinates
-    double  cos_omega = std::cos(omega);
-    double  sin_omega = std::sin(omega);
-    GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
-
-    // Rotate from native into celestial system
-    GVector cel = m_rot * native;
+    // Compute azimuth angle in model coordinate system (radians)
+    double omega_model = omega + m_posangle_roi;
-    // Set sky direction
-    GSkyDir srcDir;
-    srcDir.celvector(cel);
+    // Evaluate sky model
+    double model = m_model.eval(m_rho, omega_model, m_srcEng, m_srcTime);
+
+    // Debug: test if model is non positive
+    #if defined(G_DEBUG_MODEL_ZERO)
+    if (model <= 0.0) {
+        double m_semiminor_rad = m_model.semiminor() * gammalib::deg2rad;
+        double m_semimajor_rad = m_model.semimajor() * gammalib::deg2rad;
+        double diff_angle      = omega_model - m_model.posangle() * gammalib::deg2rad;
+        double cosinus         = std::cos(diff_angle);
+        double sinus           = std::sin(diff_angle);
+        double arg1            = m_semiminor_rad * cosinus;
+        double arg2            = m_semimajor_rad * sinus;
+        double r_ellipse       = m_semiminor_rad * m_semimajor_rad /
+                                 std::sqrt(arg1*arg1 + arg2*arg2);
+        std::cout << "*** WARNING: cta_npred_elliptical_kern_omega::eval";
+        std::cout << " zero model for (rho,omega)=(";
+        std::cout << m_rho*gammalib::rad2deg << ",";
+        std::cout << omega*gammalib::rad2deg << ")";
+        std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
+        std::cout << std::endl;
+    }
+    #endif
+    
+    // Continue only if model is positive
+    if (model > 0.0) {
+    
+        // Compute sky direction vector in native coordinates
+        double  cos_omega = std::cos(omega_model);
+        double  sin_omega = std::sin(omega_model);
+        GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
-    // Set Photon
-    GPhoton photon(srcDir, m_srcEng, m_srcTime);
+        // Rotate from native into celestial system
+        GVector cel = m_rot * native;
-    // Get model value for this photon
-    double model = m_model.eval(photon);
+        // Set sky direction
+        GSkyDir srcDir;
+        srcDir.celvector(cel);
-    // Continue only if sky intensity is positive
-    if (model > 0.0) {
+        // Set Photon
+        GPhoton photon(srcDir, m_srcEng, m_srcTime);
         // Compute Npred for this sky direction
         npred = m_rsp.npred(photon, m_obs) * model;
@@ -1130,136 +1514,203 @@ double cta_npred_diffuse_kern_phi::eval(const double&amp; phi)
 }
+/*==========================================================================
+ =                                                                         =
+ =                Helper class methods for stacked analysis                =
+ =                                                                         =
+ ==========================================================================*/
+
 /***********************************************************************//**
- * @brief Kernel for PSF integration of spatial model
+ * @brief Kernel for radial model integration over zenith angle
  *
- * @param[in] delta PSF offset angle (radians).
- * @return Azimuthally integrated product between PSF and model.
+ * @param[in] rho Radial distance from model centre (radians).
+ * @return Integration kernel.
  *
- * 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
+ * Computes the integration kernel
  *
  * \f[
- *    {\rm PSF}(\delta) \times
- *    \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi
+ *    K(\rho | E, t) = \sin \rho \times
+ *                     S_{\rm p}(\rho | E, t) \times
+ *                     \int_{\omega} PSF(\rho, \omega) d\omega
  * \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.
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the observed photon direction, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ * \f$S_{\rm p}(\rho | E, t)\f$ is the radial source model
+ * for a given true photon energy and photon arrival time,
+ * \f$PSF(\rho, \omega)\f$ is the point spread function.
  ***************************************************************************/
-double cta_psf_diffuse_kern_delta::eval(const double& delta)
+double cta_psf_radial_kern_rho::eval(const double& rho)
 {
-    // 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);
+    // Initialise result
+    double irf = 0.0;
-        // Continue only if PSF is positive
-        if (value > 0.0) {
+    // Continue only if rho is positive
+    if (rho > 0.0) {
-            // Compute sine and cosine of delta
-            double sin_delta = std::sin(delta);
-            double cos_delta = std::cos(delta);
+        // Compute half length of the arc (in radians) from a circle with
+        // radius rho that intersects with the point spread function, defined
+        // as a circle with maximum radius m_delta_max
+        double domega = 0.5 * gammalib::cta_roi_arclength(rho,
+                                                          m_rho_obs,
+                                                          m_cos_rho_obs,
+                                                          m_sin_rho_obs,
+                                                          m_delta_max,
+                                                          m_cos_delta_max);
-            // 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);
+        // Continue only if arc length is positive
+        if (domega > 0.0) {
-            // Azimuthally integrate model
-            GIntegral integral(&integrand);
-            integral.eps(m_eps);
-            value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
+            // 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;
+            }
-        } // endif: PSF value was positive
+            // Evaluate sky model
+            double model = m_model->eval(rho_kluge, m_srcEng, m_srcTime);
+
+            // Debug: test if model is non positive
+            #if defined(G_DEBUG_MODEL_ZERO)
+            if (model <= 0.0) {
+                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 << " radians" << std::endl;
+            }
+            #endif
-    } // endif: delta was positive
+            // Continue only if model is positive
+            if (model > 0.0) {
-    // 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
+                // Compute omega integration range
+                double omega_min = -domega;
+                double omega_max = +domega;
-    // Return kernel value
-    return value;
-}
+                // Precompute cosine and sine terms for azimuthal integration
+                double cos_rho = std::cos(rho);
+                double sin_rho = std::sin(rho);
+                double cos_psf = cos_rho * m_cos_rho_obs;
+                double sin_psf = sin_rho * m_sin_rho_obs;
+                // Setup integration kernel
+                cta_psf_radial_kern_omega integrand(m_rsp,
+                                                    m_model,
+                                                    m_srcDir,
+                                                    m_srcEng,
+                                                    m_srcTime,
+                                                    cos_psf,
+                                                    sin_psf);
-/***********************************************************************//**
- * @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);
+                // Setup integrator
+                GIntegral integral(&integrand);
+                integral.fixed_iter(m_iter);
-    // Rotate from native into celestial system
-    GVector cel = m_rot * native;
+                // Integrate over omega
+                irf = integral.romberg(omega_min, omega_max, m_iter) *
+                      model * sin_rho;
-    // Set sky direction
-    GSkyDir srcDir;
-    srcDir.celvector(cel);
+                // Compile option: Check for NaN/Inf
+                #if defined(G_NAN_CHECK)
+                if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
+                    std::cout << "*** ERROR: cta_psf_radial_kern_rho";
+                    std::cout << "(rho=" << rho << "):";
+                    std::cout << " NaN/Inf encountered";
+                    std::cout << " (irf=" << irf;
+                    std::cout << ", domega=" << domega;
+                    std::cout << ", model=" << model;
+                    std::cout << ", sin_rho=" << sin_rho << ")";
+                    std::cout << std::endl;
+                }
+                #endif
-    // Set photon
-    GPhoton photon(srcDir, m_srcEng, m_srcTime);
+            } // endif: model was positive
-    // Compute map value this sky direction
-    double value = m_model->eval(photon); 
+        } // endif: arclength 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_phi::eval";
-        std::cout << "(phi=" << phi << "):";
-        std::cout << " NaN/Inf encountered";
-        std::cout << " (value=" << value;
-        std::cout << ")" << std::endl;
-    }
-    #endif
+    } // endif: rho was positive
-    // Return kernel value
-    return value;
+    // Return result
+    return irf;
 }
 /***********************************************************************//**
- * @brief Kernel for PSF integration of radial model
+ * @brief Kernel for radial model integration over azimuth angle
  *
- * @param[in] delta PSF offset angle (radians).
- * @return Azimuthally integrated product between PSF and radial model.
+ * @param[in] omega Azimuth angle (radians).
  *
- * 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
+ * Computes
  *
  * \f[
- *    {\rm PSF}(\delta) \times
+ *    K(\omega | \rho, E, t) = PSF(\omega | \rho)
+ * \f]
+ *
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the observed photon direction, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ *
+ * From the coordinates \f$(\rho,\omega)\f$ in the model system, the method
+ * computes the angle between the true (\f$\vec{p}\f$) and observed
+ * (\f$\vec{p'}\f$) photon arrival direction using
+ * 
+ * \f[\delta = \arccos(\cos \rho \cos \rho_{\rm obs} + 
+ *                     \sin \rho \sin \rho_{\rm obs} \cos \omega)\f]
+ *
+ * where
+ * \f$\rho_{\rm obs}\f$ is the angular distance between the observed photon
+ * arrival direction \f$\vec{p'}\f$ and the model centre \f$\vec{m}\f$.
+ * \f$\delta\f$ is used to compute the value of the point spread function.
+ ***************************************************************************/
+double cta_psf_radial_kern_omega::eval(const double& omega)
+{
+    // Initialise Irf value
+    double irf = 0.0;
+
+    // Compute Psf offset angle (radians)
+    double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
+
+    // Evaluate Psf * model for this delta
+    irf = m_rsp->psf()(m_srcDir, delta, m_srcEng);
+
+    // Compile option: Check for NaN/Inf
+    #if defined(G_NAN_CHECK)
+    if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
+        std::cout << "*** ERROR: cta_psf_radial_kern_omega::eval";
+        std::cout << "(omega=" << omega << "):";
+        std::cout << " NaN/Inf encountered";
+        std::cout << " (irf=" << irf;
+        std::cout << ", delta=" << delta;
+        std::cout << ", omega=" << omega << ")";
+        std::cout << std::endl;
+    }
+    #endif
+
+    // Return
+    return irf;
+}
+
+
+/***********************************************************************//**
+ * @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]
  *
@@ -1273,44 +1724,70 @@ double cta_psf_radial_kern_delta::eval(const double&amp; 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
-        value = m_rsp->psf()(m_srcDir, delta, m_srcEng);
+        // Get Psf for this delta
+        double psf = m_rsp->psf()(m_srcDir, delta, m_srcEng);
-        // Continue only if PSF is positive
-        if (value > 0.0) {
+        // Continue only if Psf is positive
+        if (psf > 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;
+            // 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);
-            // 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);
+            // Continue only if arc length is positive
+            if (dphi > 0.0) {
-            // Azimuthally integrate model
-            GIntegral integral(&integrand);
-            integral.eps(m_eps);
-            value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
+                // Compute phi integration range
+                double phi_min = -dphi;
+                double phi_max = +dphi;
-        } // endif: PSF value was positive
+                // 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;
-    } // endif: delta was positive
+                // 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);
-    // 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
+                // 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;
@@ -1343,8 +1820,26 @@ double cta_psf_radial_kern_phi::eval(const double&amp; 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, m_srcEng, m_srcTime); 
+    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)
@@ -1363,238 +1858,330 @@ double cta_psf_radial_kern_phi::eval(const double&amp; phi)
 /***********************************************************************//**
- * @brief Kernel for PSF integration of elliptical model
+ * @brief Kernel for elliptical model integration over zenith angle
  *
- * @param[in] delta PSF offset angle (radians).
- * @return Azimuthally integrated product between PSF and elliptical model.
+ * @param[in] rho Radial distance from model centre (radians).
+ * @return Integration kernel.
  *
- * 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
+ * Computes the integration kernel
  *
  * \f[
- *    {\rm PSF}(\delta) \times
- *    \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi
+ *    K(\rho | E, t) = \sin \rho \times
+ *                     \int_{\omega} 
+ *                     S_{\rm p}(\rho, \omega | E, t) \, PSF(\rho, \omega)
+ *                     d\omega
  * \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.
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the observed photon direction, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical source model
+ * for a given true photon energy and photon arrival time,
+ * \f$PSF(\rho, \omega)\f$ is the point spread function.
+ * The integration is over the \f$\omega\f$ values that are comprised within
+ * the ellipse boundaries.
  ***************************************************************************/
-double cta_psf_elliptical_kern_delta::eval(const double& delta)
+double cta_psf_elliptical_kern_rho::eval(const double& rho)
 {
-    // Initialise value
-    double value = 0.0;
+    // Initialise result
+    double irf = 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) {
+    // Continue only if rho is positive
+    if (rho > 0.0) {
-        // Get PSF for this delta
-        value = m_rsp->psf()(m_srcDir, delta, m_srcEng);
+        // Compute half length of the arc (in radians) from a circle with
+        // radius rho that intersects with the point spread function, defined
+        // as a circle with maximum radius m_delta_max
+        double domega = 0.5 * gammalib::cta_roi_arclength(rho,
+                                                          m_rho_obs,
+                                                          m_cos_rho_obs,
+                                                          m_sin_rho_obs,
+                                                          m_delta_max,
+                                                          m_cos_delta_max);
-        // Continue only if PSF is positive
-        if (value > 0.0) {
+        // Continue only if arc length is positive
+        if (domega > 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;
+            // Precompute cosine and sine terms for azimuthal integration
+            double cos_rho = std::cos(rho);
+            double sin_rho = std::sin(rho);
+            double cos_psf = cos_rho * m_cos_rho_obs;
+            double sin_psf = sin_rho * m_sin_rho_obs;
+
+            // Reduce rho by an infinite amount to avoid rounding errors
+            // at the boundary of a sharp edged model (e.g. an elliptical
+            // disk model)
+            double rho_kluge = rho - 1.0e-12;
+            if (rho_kluge < 0.0) {
+                rho_kluge = 0.0;
+            }
-            // 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);
+            // Setup integration kernel
+            cta_psf_elliptical_kern_omega integrand(m_rsp,
+                                                    m_model,
+                                                    m_srcDir,
+                                                    m_srcEng,
+                                                    m_srcTime,
+                                                    m_posangle_obs,
+                                                    rho_kluge,
+                                                    cos_psf,
+                                                    sin_psf);
-            // Azimuthally integrate model
+            // Setup integrator
             GIntegral integral(&integrand);
-            integral.eps(m_eps);
-            value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
+            integral.fixed_iter(m_iter);
-        } // endif: PSF value was positive
+            // 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) {
-    } // endif: delta was positive
+                // Compute omega integration range
+                double omega_min = -domega;
+                double omega_max = +domega;
-    // 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
+                // Integrate over omega
+                irf = integral.romberg(omega_min, omega_max, m_iter) *
+                      sin_rho;
+
+            } // endif: circle comprised in ellipse
+
+            // ... otherwise there are arcs that intersect with the Psf circle
+            else {
+
+                // Compute half the arc length (in radians) of a circle of
+                // radius rho, centred on the model, that intersects with
+                // the ellipse boundary
+                double arg1 = 1.0 - (m_semiminor*m_semiminor) / (rho*rho);
+                double arg2 = 1.0 - (m_semiminor*m_semiminor) /
+                                    (m_semimajor*m_semimajor);
+                double omega_width = std::acos(std::sqrt(arg1/arg2));
+
+                // Continue only if the arclength is positive
+                if (omega_width > 0.0) {
+
+                    // Compute azimuth angle difference between ellipse
+                    // position angle and position angle of observed
+                    // photon in the model system. This angle will define
+                    // the reference point around which the circle arc. Make
+                    // sure that omega_0 is within [-pi,pi] thus that the omega
+                    // intervals are between [-2pi,2pi]
+                    double omega_0 = m_posangle - m_posangle_obs;
+                    if (omega_0 > gammalib::pi) {
+                        omega_0 -= gammalib::pi;
+                    }
+                    else if (omega_0 < -gammalib::pi) {
+                        omega_0 += gammalib::pi;
+                    }
+
+                    // Compute azimuth angle intervals
+                    double omega1_min = omega_0    - omega_width;
+                    double omega1_max = omega_0    + omega_width;
+                    double omega2_min = omega1_min + gammalib::pi;
+                    double omega2_max = omega1_max + gammalib::pi;
+
+                    // Limit intervals to the intersection of the ellipse with
+                    // the Psf circle. This may lead to a split of intervals,
+                    // and we gather all these intervals in a special interval
+                    // pair containers
+                    cta_omega_intervals intervals1 = 
+                        gammalib::limit_omega(omega1_min, omega1_max, domega);
+                    cta_omega_intervals intervals2 = 
+                        gammalib::limit_omega(omega2_min, omega2_max, domega);
+
+                    // Integrate over all intervals for omega1
+                    for (int i = 0; i < intervals1.size(); ++i) {
+                        double min = intervals1[i].first;
+                        double max = intervals1[i].second;
+                        irf       += integral.romberg(min, max, m_iter) * sin_rho;
+                    }
+
+                    // Integrate over all intervals for omega2
+                    for (int i = 0; i < intervals2.size(); ++i) {
+                        double min = intervals2[i].first;
+                        double max = intervals2[i].second;
+                        irf       += integral.romberg(min, max, m_iter) * sin_rho;
+                    }
+
+                } // endif: arc length was positive
+
+            } // endelse: circle was not comprised in ellipse
-    // Return kernel value
-    return value;
+            // Compile option: Check for NaN/Inf
+            #if defined(G_NAN_CHECK)
+            if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
+                std::cout << "*** ERROR: cta_psf_elliptical_kern_rho";
+                std::cout << "(rho=" << rho << "):";
+                std::cout << " NaN/Inf encountered";
+                std::cout << " (irf=" << irf;
+                std::cout << ", domega=" << domega;
+                std::cout << ", sin_rho=" << sin_rho << ")";
+                std::cout << std::endl;
+            }
+            #endif
+
+        } // endif: arc length was positive
+    
+    } // endif: rho was positive
+
+    // Return result
+    return irf;
 }
 /***********************************************************************//**
- * @brief Kernel for azimuthal elliptical model integration
+ * @brief Kernel for elliptical model integration over azimuth angle
  *
- * @param[in] phi Azimuth angle (radians).
- * @return Radial model value.
+ * @param[in] omega Azimuth angle (radians).
  *
- * 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
+ * Computes
  *
  * \f[
- *    \theta = \arccos \left( \cos \delta \cos \zeta +
- *                            \sin \delta \sin \zeta \cos \phi \right)
+ *    K(\omega | \rho, E, t) = S_{\rm p}(\omega | \rho, E, t) \,
+ *                             PSF(\omega | \rho)
  * \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
+ * where
+ * \f$\omega\f$ is the azimuth angle with respect to the model centre,
+ * counted counterclockwise from the vector connecting the model centre
+ * to the observed photon direction, and
+ * \f$\rho\f$ is the radial distance from the model centre.
+ *
+ * From the coordinates \f$(\rho,\omega)\f$ in the model system, the method
+ * computes the angle between the true (\f$\vec{p}\f$) and observed
+ * (\f$\vec{p'}\f$) photon arrival direction using
+ * 
+ * \f[\delta = \arccos(\cos \rho \cos \rho_{\rm obs} + 
+ *                     \sin \rho \sin \rho_{\rm obs} \cos \omega)\f]
  *
- * \f[
- *    {\rm posangle} = \arcsin \left( \frac{\sin \phi \sin \delta}
- *                                         {\sin \theta} \right)
- * \f]
+ * where
+ * \f$\rho_{\rm obs}\f$ is the angular distance between the observed photon
+ * arrival direction \f$\vec{p'}\f$ and the model centre \f$\vec{m}\f$.
+ * \f$\delta\f$ is used to compute the value of the point spread function.
  ***************************************************************************/
-double cta_psf_elliptical_kern_phi::eval(const double& phi)
+double cta_psf_elliptical_kern_omega::eval(const double& omega)
 {
-    // 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;
+    // Initialise Irf value
+    double irf = 0.0;
-    // Get radial model value
-    double value = m_model->eval(theta, posangle, m_srcEng, m_srcTime);
+    // Compute azimuth angle in model coordinate system (radians)
+    double omega_model = omega + m_posangle_obs;
-    // 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;
+    // Evaluate sky model
+    double model = m_model->eval(m_rho, omega_model, m_srcEng, m_srcTime);
+
+    // Debug: test if model is non positive
+    #if defined(G_DEBUG_MODEL_ZERO)
+    if (model <= 0.0) {
+        double semiminor_rad = m_model->semiminor() * gammalib::deg2rad;
+        double semimajor_rad = m_model->semimajor() * gammalib::deg2rad;
+        double diff_angle    = omega_model - m_model->posangle() * gammalib::deg2rad;
+        double cosinus       = std::cos(diff_angle);
+        double sinus         = std::sin(diff_angle);
+        double arg1          = semiminor_rad * cosinus;
+        double arg2          = semimajor_rad * sinus;
+        double r_ellipse     = semiminor_rad * semimajor_rad /
+                               std::sqrt(arg1*arg1 + arg2*arg2);
+        std::cout << "*** WARNING: cta_psf_elliptical_kern_omega::eval";
+        std::cout << " zero model for (rho,omega)=";
+        std::cout << m_rho*gammalib::rad2deg << ",";
+        std::cout << omega*gammalib::rad2deg << ")";
+        std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
+        std::cout << std::endl;
     }
     #endif
-    // Return kernel value
-    return value;
-}
+    // Continue only if model is positive
+    if (model > 0.0) {
+        // Compute Psf offset angle (radians)
+        double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(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.                                                   =
- ==========================================================================*/
+        // Evaluate Psf * model for this delta
+        irf = m_rsp->psf()(m_srcDir, delta, m_srcEng) * model;
+
+        // Compile option: Check for NaN/Inf
+        #if defined(G_NAN_CHECK)
+        if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
+            std::cout << "*** ERROR: cta_psf_elliptical_kern_omega::eval";
+            std::cout << "(omega=" << omega << "):";
+            std::cout << " NaN/Inf encountered";
+            std::cout << " (irf=" << irf;
+            std::cout << ", model=" << model;
+            std::cout << ", delta=" << delta;
+            std::cout << ", rho=" << m_rho;
+            std::cout << ", omega=" << omega << ")";
+            std::cout << std::endl;
+        }
+        #endif
+
+    } // endif: model is positive
+
+    // Return
+    return irf;
+}
 /***********************************************************************//**
- * @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;
+            // 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 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;
+            // 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;
@@ -1608,79 +2195,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 \right)
- * \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;
@@ -29,6 +29,8 @@
 /* __ Includes ___________________________________________________________ */
 #include <cmath>
+#include <vector>
+#include <utility>
 #include "GCTAResponseIrf.hpp"
 #include "GCTAObservation.hpp"
 #include "GMatrix.hpp"
@@ -44,6 +46,16 @@
 /* __ Forward declaration ________________________________________________ */
+/* __ Typedefs ___________________________________________________________ */
+typedef std::vector<std::pair<double,double> > cta_omega_intervals;
+
+/* __ Prototypes _________________________________________________________ */
+namespace gammalib {
+    cta_omega_intervals limit_omega(const double& min,
+                                    const double& max,
+                                    const double& domega);
+}
+
 /***********************************************************************//**
  * @class cta_npsf_kern_rad_azsym
@@ -433,55 +445,64 @@ class cta_irf_elliptical_kern_rho : public GFunction {
 public:
     cta_irf_elliptical_kern_rho(const GCTAResponseIrf&         rsp,
                                 const GModelSpatialElliptical& model,
+                                const double&                  semimajor,
+                                const double&                  semiminor,
+                                const double&                  posangle,
                                 const double&                  zenith,
                                 const double&                  azimuth,
                                 const GEnergy&                 srcEng,
                                 const GTime&                   srcTime,
                                 const double&                  srcLogEng,
                                 const GEnergy&                 obsEng,
-                                const double&                  zeta,
-                                const double&                  lambda,
-                                const double&                  obsOmega,
-                                const double&                  omega0,
+                                const double&                  rho_obs,
+                                const double&                  posangle_obs,
+                                const double&                  rho_pnt,
+                                const double&                  omega_pnt,
                                 const double&                  delta_max,
                                 const int&                     iter) :
                                 m_rsp(rsp),
                                 m_model(model),
+                                m_semimajor(semimajor),
+                                m_semiminor(semiminor),
+                                m_posangle(posangle),
                                 m_zenith(zenith),
                                 m_azimuth(azimuth),
                                 m_srcEng(srcEng),
                                 m_srcTime(srcTime),
                                 m_srcLogEng(srcLogEng),
                                 m_obsEng(obsEng),
-                                m_zeta(zeta),
-                                m_cos_zeta(std::cos(zeta)),
-                                m_sin_zeta(std::sin(zeta)),
-                                m_lambda(lambda),
-                                m_cos_lambda(std::cos(lambda)),
-                                m_sin_lambda(std::sin(lambda)),
-                                m_obsOmega(obsOmega),
-                                m_omega0(omega0),
+                                m_rho_obs(rho_obs),
+                                m_cos_rho_obs(std::cos(rho_obs)),
+                                m_sin_rho_obs(std::sin(rho_obs)),
+                                m_posangle_obs(posangle_obs),
+                                m_rho_pnt(rho_pnt),
+                                m_cos_rho_pnt(std::cos(rho_pnt)),
+                                m_sin_rho_pnt(std::sin(rho_pnt)),
+                                m_omega_pnt(omega_pnt),
                                 m_delta_max(delta_max),
                                 m_cos_delta_max(std::cos(delta_max)),
                                 m_iter(iter) { }
     double eval(const double& rho);
 public:
     const GCTAResponseIrf&         m_rsp;           //!< CTA response
-    const GModelSpatialElliptical& m_model;         //!< Spatial model
+    const GModelSpatialElliptical& m_model;         //!< Elliptical model
+    const double&                  m_semimajor;     //!< Ellipse boundary semimajor axis
+    const double&                  m_semiminor;     //!< Ellipse boundary semiminor axis
+    const double&                  m_posangle;      //!< Ellipse boundary position angle
     const double&                  m_zenith;        //!< Zenith angle
     const double&                  m_azimuth;       //!< Azimuth angle
     const GEnergy&                 m_srcEng;        //!< True photon energy
     const GTime&                   m_srcTime;       //!< True photon time
     const double&                  m_srcLogEng;     //!< True photon log energy
     const GEnergy&                 m_obsEng;        //!< Measured event energy
-    const double&                  m_zeta;          //!< Distance model centre - measured photon
-    double                         m_cos_zeta;      //!< Cosine of zeta
-    double                         m_sin_zeta;      //!< Sine of zeta
-    const double&                  m_lambda;        //!< Distance model centre - pointing
-    double                         m_cos_lambda;    //!< Cosine of lambda
-    double                         m_sin_lambda;    //!< Sine of lambda
-    const double&                  m_obsOmega;      //!< Measured photon position angle from model centre
-    const double&                  m_omega0;        //!< Azimuth of pointing in model system
+    const double&                  m_rho_obs;       //!< Distance of model centre from measured photon
+    double                         m_cos_rho_obs;   //!< Cosine of m_rho_obs
+    double                         m_sin_rho_obs;   //!< Sine of m_rho_obs
+    const double&                  m_posangle_obs;  //!< Photon position angle measured from model centre
+    const double&                  m_rho_pnt;       //!< Distance of model centre from pointing
+    double                         m_cos_rho_pnt;   //!< Cosine of m_rho_pnt
+    double                         m_sin_rho_pnt;   //!< Sine of m_rho_pnt
+    const double&                  m_omega_pnt;     //!< Azimuth of pointing in model system
     const double&                  m_delta_max;     //!< Maximum PSF radius
     double                         m_cos_delta_max; //!< Cosine of maximum PSF radius
     const int&                     m_iter;          //!< Integration iterations
@@ -517,8 +538,8 @@ public:
                                   const GTime&                   srcTime,
                                   const double&                  srcLogEng,
                                   const GEnergy&                 obsEng,
-                                  const double&                  obsOmega,
-                                  const double&                  omega0,
+                                  const double&                  posangle_obs,
+                                  const double&                  omega_pnt,
                                   const double&                  rho,
                                   const double&                  cos_psf,
                                   const double&                  sin_psf,
@@ -532,8 +553,8 @@ public:
                                   m_srcTime(srcTime),
                                   m_srcLogEng(srcLogEng),
                                   m_obsEng(obsEng),
-                                  m_obsOmega(obsOmega),
-                                  m_omega0(omega0),
+                                  m_posangle_obs(posangle_obs),
+                                  m_omega_pnt(omega_pnt),
                                   m_rho(rho),
                                   m_cos_psf(cos_psf),
                                   m_sin_psf(sin_psf),
@@ -541,21 +562,21 @@ public:
                                   m_sin_ph(sin_ph) { }
     double eval(const double& omega);
 public:
-    const GCTAResponseIrf&         m_rsp;        //!< CTA response
-    const GModelSpatialElliptical& m_model;      //!< Spatial model
-    const double&                  m_zenith;     //!< Zenith angle
-    const double&                  m_azimuth;    //!< Azimuth angle
-    const GEnergy&                 m_srcEng;     //!< True photon energy
-    const GTime&                   m_srcTime;    //!< True photon time
-    const double&                  m_srcLogEng;  //!< True photon log energy
-    const GEnergy&                 m_obsEng;     //!< Measured event energy
-    const double&                  m_obsOmega;   //!< Measured photon position angle from model centre
-    const double&                  m_omega0;     //!< Azimuth of pointing in model system
-    const double&                  m_rho;        //!< Model zenith angle
-    const double&                  m_cos_psf;    //!< Cosine term for PSF offset angle computation
-    const double&                  m_sin_psf;    //!< Sine term for PSF offset angle computation
-    const double&                  m_cos_ph;     //!< Cosine term for photon offset angle computation
-    const double&                  m_sin_ph;     //!< Sine term for photon offset angle computation
+    const GCTAResponseIrf&         m_rsp;          //!< CTA response
+    const GModelSpatialElliptical& m_model;        //!< Spatial model
+    const double&                  m_zenith;       //!< Zenith angle
+    const double&                  m_azimuth;      //!< Azimuth angle
+    const GEnergy&                 m_srcEng;       //!< True photon energy
+    const GTime&                   m_srcTime;      //!< True photon time
+    const double&                  m_srcLogEng;    //!< True photon log energy
+    const GEnergy&                 m_obsEng;       //!< Measured event energy
+    const double&                  m_posangle_obs; //!< Measured photon position angle from model centre
+    const double&                  m_omega_pnt;    //!< Azimuth of pointing in model system
+    const double&                  m_rho;          //!< Model zenith angle
+    const double&                  m_cos_psf;      //!< Cosine term for PSF offset angle computation
+    const double&                  m_sin_psf;      //!< Sine term for PSF offset angle computation
+    const double&                  m_cos_ph;       //!< Cosine term for photon offset angle computation
+    const double&                  m_sin_ph;       //!< Sine term for photon offset angle computation
 };
@@ -594,42 +615,51 @@ class cta_npred_elliptical_kern_rho : public GFunction {
 public:
     cta_npred_elliptical_kern_rho(const GCTAResponseIrf&         rsp,
                                   const GModelSpatialElliptical& model,
+                                  const double&                  semimajor,
+                                  const double&                  semiminor,
+                                  const double&                  posangle,
                                   const GEnergy&                 srcEng,
                                   const GTime&                   srcTime,
                                   const GCTAObservation&         obs,
                                   const GMatrix&                 rot,
-                                  const double&                  dist,
-                                  const double&                  radius,
-                                  const double&                  omega0,
+                                  const double&                  rho_roi,
+                                  const double&                  posangle_roi,
+                                  const double&                  radius_roi,
                                   const int&                     iter) :
                                   m_rsp(rsp),
                                   m_model(model),
+                                  m_semimajor(semimajor),
+                                  m_semiminor(semiminor),
+                                  m_posangle(posangle),
                                   m_srcEng(srcEng),
                                   m_srcTime(srcTime),
                                   m_obs(obs),
                                   m_rot(rot),
-                                  m_dist(dist),
-                                  m_cos_dist(std::cos(dist)),
-                                  m_sin_dist(std::sin(dist)),
-                                  m_radius(radius),
-                                  m_cos_radius(std::cos(radius)),
-                                  m_omega0(omega0),
+                                  m_rho_roi(rho_roi),
+                                  m_cos_rho_roi(std::cos(rho_roi)),
+                                  m_sin_rho_roi(std::sin(rho_roi)),
+                                  m_posangle_roi(posangle_roi),
+                                  m_radius_roi(radius_roi),
+                                  m_cos_radius_roi(std::cos(radius_roi)),
                                   m_iter(iter) { }
     double eval(const double& rho);
 protected:
-    const GCTAResponseIrf&         m_rsp;        //!< CTA response
-    const GModelSpatialElliptical& m_model;      //!< Spatial model
-    const GEnergy&                 m_srcEng;     //!< True photon energy
-    const GTime&                   m_srcTime;    //!< True photon arrival time
-    const GCTAObservation&         m_obs;        //!< CTA observation
-    const GMatrix&                 m_rot;        //!< Rotation matrix
-    const double&                  m_dist;       //!< Distance model-ROI centre
-    double                         m_cos_dist;   //!< Cosine of distance model-ROI centre
-    double                         m_sin_dist;   //!< Sine of distance model-ROI centre
-    const double&                  m_radius;     //!< ROI+PSF radius
-    double                         m_cos_radius; //!< Cosine of ROI+PSF radius
-    const double&                  m_omega0;     //!< Position angle of ROI
-    const int&                     m_iter;       //!< Integration iterations
+    const GCTAResponseIrf&         m_rsp;            //!< CTA response
+    const GModelSpatialElliptical& m_model;          //!< Elliptical model
+    const double&                  m_semimajor;      //!< Ellipse boundary semimajor axis
+    const double&                  m_semiminor;      //!< Ellipse boundary semiminor axis
+    const double&                  m_posangle;       //!< Ellipse boundary position angle
+    const GEnergy&                 m_srcEng;         //!< True photon energy
+    const GTime&                   m_srcTime;        //!< True photon arrival time
+    const GCTAObservation&         m_obs;            //!< CTA observation
+    const GMatrix&                 m_rot;            //!< Rotation matrix
+    const double&                  m_rho_roi;        //!< Distance between model and ROI centre
+    double                         m_cos_rho_roi;    //!< Cosine of m_rho_roi
+    double                         m_sin_rho_roi;    //!< Sine of m_rho_roi
+    const double&                  m_posangle_roi;   //!< Position angle of ROI
+    const double&                  m_radius_roi;     //!< ROI+PSF radius
+    double                         m_cos_radius_roi; //!< Cosine of m_radius_roi
+    const int&                     m_iter;           //!< Integration iterations
 };
@@ -663,26 +693,32 @@ public:
                                     const GTime&                   srcTime,
                                     const GCTAObservation&         obs,
                                     const GMatrix&                 rot,
+                                    const double&                  rho,
                                     const double&                  sin_rho,
-                                    const double&                  cos_rho) :
+                                    const double&                  cos_rho,
+                                    const double&                  posangle_roi) :
                                     m_rsp(rsp),
                                     m_model(model),
                                     m_srcEng(srcEng),
                                     m_srcTime(srcTime),
                                     m_obs(obs),
                                     m_rot(rot),
+                                    m_rho(rho),
                                     m_sin_rho(sin_rho),
-                                    m_cos_rho(cos_rho) { }
+                                    m_cos_rho(cos_rho),
+                                    m_posangle_roi(posangle_roi) { }
     double eval(const double& omega);
 protected:
-    const GCTAResponseIrf&         m_rsp;      //!< CTA response
-    const GModelSpatialElliptical& m_model;    //!< Model
-    const GEnergy&                 m_srcEng;   //!< True photon energy
-    const GTime&                   m_srcTime;  //!< True photon arrival time
-    const GCTAObservation&         m_obs;      //!< Pointer to observation
-    const GMatrix&                 m_rot;      //!< Rotation matrix
-    const double&                  m_sin_rho;  //!< Sine of offset angle
-    const double&                  m_cos_rho;  //!< Cosine of offset angle
+    const GCTAResponseIrf&         m_rsp;          //!< CTA response
+    const GModelSpatialElliptical& m_model;        //!< Model
+    const GEnergy&                 m_srcEng;       //!< True photon energy
+    const GTime&                   m_srcTime;      //!< True photon arrival time
+    const GCTAObservation&         m_obs;          //!< Pointer to observation
+    const GMatrix&                 m_rot;          //!< Rotation matrix
+    const double&                  m_rho;
+    const double&                  m_sin_rho;      //!< Sine of offset angle
+    const double&                  m_cos_rho;      //!< Cosine of offset angle
+    const double&                  m_posangle_roi; //!< Position angle of ROI
 };
@@ -946,70 +982,111 @@ protected:
 /***********************************************************************//**
- * @class cta_psf_diffuse_kern_delta
+ * @class cta_irf_radial_kern_rho
  *
- * @brief Kernel for Psf delta angle integration used for stacked analysis
+ * @brief Kernel for radial model zenith angle integration
+ *
+ * This class implements the integration kernel \f$K(\rho)\f$ for the
+ * integration
+ *
+ * \f[
+ *    \int_{\rho_{\rm min}}^{\rho_{\rm max}} K(\rho | E, t) d\rho
+ * \f]
+ *
+ * of radial spatial models. The eval() method computes
+ *
+ * \f[
+ *    K(\rho | E, t) = \sin \rho \times
+ *                     S_{\rm p}(\rho | E, t) \times
+ *                     \int_{\omega} PSF(\rho, \omega) d\omega
+ * \f]
+ *
+ * where
+ * \f$S_{\rm p}(\rho | E, t)\f$ is the radial 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.
  ***************************************************************************/
-class cta_psf_diffuse_kern_delta : public GFunction {
+class cta_psf_radial_kern_rho : 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
+    cta_psf_radial_kern_rho(const GCTAResponseCube*    rsp,
+                            const GModelSpatialRadial* model,
+                            const GSkyDir&             srcDir,
+                            const GEnergy&             srcEng,
+                            const GTime&               srcTime,
+                            const double&              rho_obs,
+                            const double&              delta_max,
+                            const int&                 iter) :
+                            m_rsp(rsp),
+                            m_model(model),
+                            m_srcDir(srcDir),
+                            m_srcEng(srcEng),
+                            m_srcTime(srcTime),
+                            m_rho_obs(rho_obs),
+                            m_cos_rho_obs(std::cos(rho_obs)),
+                            m_sin_rho_obs(std::sin(rho_obs)),
+                            m_delta_max(delta_max),
+                            m_cos_delta_max(std::cos(delta_max)),
+                            m_iter(iter) { }
+    double eval(const double& rho);
+public:
+    const GCTAResponseCube*    m_rsp;           //!< CTA response
+    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 time
+    const double&              m_rho_obs;       //!< Distance of model centre from measured photon
+    double                     m_cos_rho_obs;   //!< Cosine of m_rho_obs
+    double                     m_sin_rho_obs;   //!< Sine of m_rho_obs
+    const double&              m_delta_max;     //!< Maximum PSF radius
+    double                     m_cos_delta_max; //!< Cosine of maximum PSF radius
+    const int&                 m_iter;          //!< Integration iterations
 };
 /***********************************************************************//**
- * @class cta_psf_diffuse_kern_phi
+ * @class cta_irf_radial_kern_omega
  *
- * @brief Kernel for Psf phi angle integration used for stacked analysis
+ * @brief Kernel for radial model azimuth angle integration
+ *
+ * This class implements the computation of
+ *
+ * \f[
+ *    K(\omega | \rho, E, t) = PSF(\omega | \rho)
+ * \f]
+ *
+ * where
+ * \f$PSF(\omega | \rho)\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.
  ***************************************************************************/
-class cta_psf_diffuse_kern_phi : public GFunction {
+class cta_psf_radial_kern_omega : 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)
+    cta_psf_radial_kern_omega(const GCTAResponseCube*    rsp,
+                              const GModelSpatialRadial* model,
+                              const GSkyDir&             srcDir,
+                              const GEnergy&             srcEng,
+                              const GTime&               srcTime,
+                              const double&              cos_psf,
+                              const double&              sin_psf) :
+                              m_rsp(rsp),
+                              m_model(model),
+                              m_srcDir(srcDir),
+                              m_srcEng(srcEng),
+                              m_srcTime(srcTime),
+                              m_cos_psf(cos_psf),
+                              m_sin_psf(sin_psf) { }
+    double eval(const double& omega);
+public:
+    const GCTAResponseCube*    m_rsp;     //!< CTA response
+    const GModelSpatialRadial* m_model;   //!< Radial model
+    const GSkyDir&             m_srcDir;  //!< True photon sky direction
+    const GEnergy&             m_srcEng;  //!< True photon energy
+    const GTime&               m_srcTime; //!< True photon time
+    const double&              m_cos_psf; //!< Cosine term for PSF offset angle computation
+    const double&              m_sin_psf; //!< Sine term for PSF offset angle computation
 };
@@ -1025,32 +1102,33 @@ public:
                               const GSkyDir&             srcDir,
                               const GEnergy&             srcEng,
                               const GTime&               srcTime,
-                              const double&              sin_zeta,
-                              const double&              cos_zeta,
-                              const double&              eps,
-                              const int&                 order) :
+                              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_sin_zeta(sin_zeta),
-                              m_cos_zeta(cos_zeta),
-                              m_eps(eps),
-                              m_order (order),
-                              m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { }
+                              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_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
+    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
 };
@@ -1076,197 +1154,205 @@ 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)
+    const double&              m_sin_fact;  //!< sin(delta)*sin(delta_mod)
+    const double&              m_cos_fact;  //!< cos(delta)*cos(delta_mod)
 };
 /***********************************************************************//**
- * @class cta_psf_elliptical_kern_delta
+ * @class cta_irf_elliptical_kern_rho
  *
- * @brief Kernel for Psf delta angle integration used for stacked analysis
+ * @brief Kernel for elliptical model zenith angle integration
+ *
+ * This class implements the integration kernel \f$K(\rho)\f$ for the
+ * integration
+ *
+ * \f[
+ *    \int_{\rho_{\rm min}}^{\rho_{\rm max}} K(\rho | E, t) d\rho
+ * \f]
+ *
+ * of elliptical spatial models. The eval() method computes
+ *
+ * \f[
+ *    K(\rho | E, t) = \sin \rho \times
+ *                     \int_{\omega} 
+ *                     S_{\rm p}(\rho, \omega | E, t) \, PSF(\rho, \omega)
+ *                     d\omega
+ * \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.
  ***************************************************************************/
-class cta_psf_elliptical_kern_delta : public GFunction {
+class cta_psf_elliptical_kern_rho : public GFunction {
 public:
-    cta_psf_elliptical_kern_delta(const GCTAResponseCube*        rsp,
+    cta_psf_elliptical_kern_rho(const GCTAResponseCube*        rsp,
+                                const GModelSpatialElliptical* model,
+                                const double&                  semimajor,
+                                const double&                  semiminor,
+                                const double&                  posangle,
+                                const GSkyDir&                 srcDir,
+                                const GEnergy&                 srcEng,
+                                const GTime&                   srcTime,
+                                const double&                  rho_obs,
+                                const double&                  posangle_obs,
+                                const double&                  delta_max,
+                                const int&                     iter) :
+                                m_rsp(rsp),
+                                m_model(model),
+                                m_semimajor(semimajor),
+                                m_semiminor(semiminor),
+                                m_posangle(posangle),
+                                m_srcDir(srcDir),
+                                m_srcEng(srcEng),
+                                m_srcTime(srcTime),
+                                m_rho_obs(rho_obs),
+                                m_cos_rho_obs(std::cos(rho_obs)),
+                                m_sin_rho_obs(std::sin(rho_obs)),
+                                m_posangle_obs(posangle_obs),
+                                m_delta_max(delta_max),
+                                m_cos_delta_max(std::cos(delta_max)),
+                                m_iter(iter) { }
+    double eval(const double& rho);
+public:
+    const GCTAResponseCube*        m_rsp;           //!< CTA response
+    const GModelSpatialElliptical* m_model;         //!< Elliptical model
+    const double&                  m_semimajor;     //!< Ellipse boundary semimajor axis
+    const double&                  m_semiminor;     //!< Ellipse boundary semiminor axis
+    const double&                  m_posangle;      //!< Ellipse boundary position angle
+    const GSkyDir&                 m_srcDir;        //!< True photon arrival direction
+    const GEnergy&                 m_srcEng;        //!< True photon energy
+    const GTime&                   m_srcTime;       //!< True photon time
+    const double&                  m_rho_obs;       //!< Distance of model centre from measured photon
+    double                         m_cos_rho_obs;   //!< Cosine of m_rho_obs
+    double                         m_sin_rho_obs;   //!< Sine of m_rho_obs
+    const double&                  m_posangle_obs;  //!< Photon position angle measured from model centre
+    const double&                  m_delta_max;     //!< Maximum PSF radius
+    double                         m_cos_delta_max; //!< Cosine of maximum PSF radius
+    const int&                     m_iter;          //!< Integration iterations
+};
+
+
+/***********************************************************************//**
+ * @class cta_irf_elliptical_kern_omega
+ *
+ * @brief Kernel for elliptical model azimuth angle integration
+ *
+ * This class implements the computation of
+ *
+ * \f[
+ *    S_{\rm p}(\omega | \rho, E, t) \, PSF(\omega | \rho)
+ * \f]
+ *
+ * where
+ * \f$S_{\rm p}(\omega | \rho, E, t)\f$ is the elliptical model,
+ * \f$ PSF(\omega | \rho)\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.
+ ***************************************************************************/
+class cta_psf_elliptical_kern_omega : public GFunction {
+public:
+    cta_psf_elliptical_kern_omega(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) :
+                                  const double&                  posangle_obs,
+                                  const double&                  rho,
+                                  const double&                  cos_psf,
+                                  const double&                  sin_psf) :
                                   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 {
+                                  m_posangle_obs(posangle_obs),
+                                  m_rho(rho),
+                                  m_cos_psf(cos_psf),
+                                  m_sin_psf(sin_psf) { }
+    double eval(const double& omega);
 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)
+    const GCTAResponseCube*        m_rsp;          //!< CTA response
+    const GModelSpatialElliptical* m_model;        //!< Spatial model
+    const GSkyDir&                 m_srcDir;       //!< True photon sky direction
+    const GEnergy&                 m_srcEng;       //!< True photon energy
+    const GTime&                   m_srcTime;      //!< True photon time
+    const double&                  m_posangle_obs; //!< Measured photon position angle from model centre
+    const double&                  m_rho;          //!< Model zenith angle
+    const double&                  m_cos_psf;      //!< Cosine term for PSF offset angle computation
+    const double&                  m_sin_psf;      //!< Sine term for PSF offset angle computation
 };
 /***********************************************************************//**
- * @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 */
@@ -59,7 +59,10 @@ def unbinned_analysis(model, evtfile, irf, caldb):
     obs.models(model)
     # Allocate Levenberg-Marquardt optimizer
-    opt = gammalib.GOptimizerLM()
+    log = gammalib.GLog()
+    log.cout(True)
+    opt = gammalib.GOptimizerLM(log)
+    opt.eps(5.0e-3)
     # Optimize model parameters
     obs.optimize(opt)
@@ -113,7 +116,10 @@ def binned_analysis(model, cntmap, irf, caldb):
     obs.models(model)
     # Allocate Levenberg-Marquardt optimizer
-    opt = gammalib.GOptimizerLM()
+    log = gammalib.GLog()
+    log.cout(True)
+    opt = gammalib.GOptimizerLM(log)
+    opt.eps(5.0e-3)
     # Optimize model parameters
     obs.optimize(opt)
@@ -169,7 +175,10 @@ def stacked_analysis(model, cntmap, expcube, psfcube):
     obs.models(model)
     # Allocate Levenberg-Marquardt optimizer
-    opt = gammalib.GOptimizerLM()
+    log = gammalib.GLog()
+    log.cout(True)
+    opt = gammalib.GOptimizerLM(log)
+    opt.eps(5.0e-3)
     # Optimize model parameters
     obs.optimize(opt)
+#! /usr/bin/env python
+# ==========================================================================
+# This script tests the elliptical response computation. The script requires
+# matplotlib to be installed.
+#
+# --------------------------------------------------------------------------
+#
+# Copyright (C) 2014 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import math
+
+
+# ================= #
+# Compute arclength #
+# ================= #
+def cta_roi_arclength(rad, dist, roi):
+    """
+    Returns length of circular arc within circular ROI
+
+    Parameters:
+     rad  - Circle radius in radians (<pi)
+     dist - Circle centre distance to ROI centre (<pi)
+     roi  - Radius of ROI in radians
+    """
+    # Initialise
+    arclength = 0.0
+
+    # Handle special case that circle centre matches ROI centre
+    if dist == 0.0:
+        if rad > roi:
+            arclength = 0.0
+        else:
+            arclength = 2.0 * math.pi
+
+    # ... otherwise circle and ROI centres are not identical
+    else:
+
+        # Handle special case that we evaluate exactly at the circle
+        # centre. In this case we have in fact a point, and if this point
+        # falls within the ROI it has a formal arclength of 2pi.
+        if rad == 0.0:
+            if dist > roi:
+                arclength = 0.0
+            else:
+                arclength = 2.0 * math.pi
+
+        # ... otherwise we have to handle the general case
+        else:
+            d = roi - dist
+            if -rad >= d:
+                arclength = 0.0
+            elif rad <= d:
+                arclength = 2.0 * math.pi
+            else:
+                cosrad  = math.cos(rad)
+                sinrad  = math.sin(rad)
+                cosroi  = math.cos(roi)
+                cosdist = math.cos(dist)
+                sindist = math.sin(dist)
+                cosang  = (cosroi - cosdist * cosrad) / (sindist * sinrad)
+                arclength = 2.0 * gammalib.acos(cosang)
+
+    # Return
+    return arclength
+
+# ============== #
+# Limit interval #
+# ============== #
+def limit_omega(min, max, domega, log=False):
+    """
+    """
+    # Compute only if domega < pi
+    if (domega < gammalib.pi):
+
+        # Save initial interval
+        min_initial = min
+        max_initial = max
+
+        # Initialise intervals
+        intervals = []
+        
+        # Limit each of the intervals to [omega_min, omega_max]
+        omega_min       = -domega
+        omega_max       = +domega
+        omega_min_plus  = omega_min + gammalib.twopi
+        omega_max_plus  = omega_max + gammalib.twopi
+        omega_min_minus = omega_min - gammalib.twopi
+        omega_max_minus = omega_max - gammalib.twopi
+
+        # Check for overlaps
+        in_unshifted = False
+        in_plus      = False
+        in_minus     = False
+        if (max_initial > omega_min and min_initial < omega_max):
+            min = min_initial
+            max = max_initial
+            in_unshifted = True
+            if (min < omega_min):
+                min = omega_min
+            if (max > omega_max):
+                max = omega_max
+            intervals.append({'min': min, 'max': max})
+        if (max_initial > omega_min_plus and min_initial < omega_max_plus):
+            min = min_initial
+            max = max_initial
+            in_plus = True
+            if (min < omega_min_plus):
+                min = omega_min_plus
+            if (max > omega_max_plus):
+                max = omega_max_plus
+            intervals.append({'min': min, 'max': max})
+        if (max_initial > omega_min_minus and min_initial < omega_max_minus):
+            min = min_initial
+            max = max_initial
+            in_minus = True
+            if (min < omega_min_minus):
+                min = omega_min_minus
+            if (max > omega_max_minus):
+                max = omega_max_minus
+            intervals.append({'min': min, 'max': max})
+
+        # Debug
+        if log:
+            print("[%7.2f,%7.2f] => [%7.2f,%7.2f] (%i[%7.2f,%7.2f]  %i[%7.2f,%7.2f]  %i[%7.2f,%7.2f])" %
+              (min_initial*gammalib.rad2deg, max_initial*gammalib.rad2deg, \
+               min*gammalib.rad2deg, max*gammalib.rad2deg, \
+               in_unshifted, \
+               omega_min*gammalib.rad2deg, omega_max*gammalib.rad2deg, \
+               in_plus, \
+               omega_min_plus*gammalib.rad2deg, omega_max_plus*gammalib.rad2deg, \
+               in_minus, \
+               omega_min_minus*gammalib.rad2deg, omega_max_minus*gammalib.rad2deg))
+
+    # ... just copy over the input interval
+    else:
+        intervals = [{'min': min, 'max': max}]
+
+    # Return
+    return intervals
+    
+
+# ======================= #
+# Show ellptical response #
+# ======================= #
+def show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, delta_max, log=False):
+    """
+    Show ellptical response
+    """
+    # Create figure
+    plt.figure(1)
+    plt.title("Elliptical response")
+
+    # Get ellipse boundary (used for intersection computation)
+    if (semimajor >= semiminor):
+        bound_major  = semimajor
+        bound_minor  = semiminor
+        bound_pa     = posang
+    else:
+        bound_major  = semiminor
+        bound_minor  = semimajor
+        bound_pa     = posang + 90.0
+
+    # Compute some stuff
+    cos_pa       = math.cos(posang * gammalib.deg2rad)
+    sin_pa       = math.sin(posang * gammalib.deg2rad)
+    rho_obs      = centre.dist_deg(obsDir)
+    posangle_obs = centre.posang(obsDir)
+    rho_pnt      = centre.dist(pntDir)
+    posangle_pnt = centre.posang(pntDir)
+    omega_pnt    = posangle_pnt - posangle_obs
+    src_max      = semimajor
+
+    # Rho limits
+    rho_min = 0.0
+    rho_max = rho_obs + delta_max
+    if (rho_obs > delta_max):
+        rho_min = rho_obs - delta_max
+    if rho_max > src_max:
+        rho_max = src_max
+
+    # Convert to radians
+    delta_max *= gammalib.deg2rad
+    rho_obs   *= gammalib.deg2rad
+    src_max   *= gammalib.deg2rad
+
+    # Setup rhos
+    rhos = []
+    for i in range(51):
+        rho = i*0.02
+        if (rho >= rho_min and rho <= rho_max):
+            rhos.append(rho)
+
+    # Setup ellipse as vector
+    x = []
+    y = []
+    for angle in range(0, 360):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(semimajor * cos_radians * sin_pa - semiminor * sin_radians * cos_pa)
+        y.append(semimajor * cos_radians * cos_pa + semiminor * sin_radians * sin_pa)
+
+    # Plot ellipse
+    plt.plot(x, y, 'r-')
+    plt.xlim([ 1.5*bound_major,-1.5*bound_major])
+    plt.ylim([-1.5*bound_major, 1.5*bound_major])
+    plt.plot([0.0], [0.0], 'ro')
+
+    # Plot observed direction
+    x_obs = obsDir.ra_deg()  - centre.ra_deg()
+    y_obs = obsDir.dec_deg() - centre.dec_deg()
+    plt.plot([x_obs], [y_obs], 'bo')
+
+    # Draw connecting line model-photon'
+    n    = 10
+    step = rho_obs/n * gammalib.rad2deg
+    x    = []
+    y    = []
+    for i in range(n+1):
+        x.append(i*step * math.sin(posangle_obs))
+        y.append(i*step * math.cos(posangle_obs))
+    plt.plot(x, y, 'b-')
+
+    # Plot Psf validity circle
+    x = []
+    y = []
+    for angle in range(0, 360):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(delta_max * cos_radians * gammalib.rad2deg + x_obs)
+        y.append(delta_max * sin_radians * gammalib.rad2deg + y_obs)
+    plt.plot(x, y, 'b-')
+
+    # Plot pointing direction
+    x = pntDir.ra_deg()  - centre.ra_deg()
+    y = pntDir.dec_deg() - centre.dec_deg()
+    plt.plot([x], [y], 'go')
+
+    # Print some stuff
+    print("Distance model-event (rho_obs) ..........: %.2f deg" % (rho_obs*gammalib.rad2deg))
+    print("PA of event from model (posangle_obs) ...: %.2f deg" % (posangle_obs*gammalib.rad2deg))
+    print("Distance model-pointing (rho_pnt) .......: %.2f deg" % (rho_pnt*gammalib.rad2deg))
+    print("PA of pointing from model (posangle_pnt) : %.2f deg" % (posangle_pnt*gammalib.rad2deg))
+    print("Azimuth of pointing (omega_pnt) .........: %.2f deg" % (omega_pnt*gammalib.rad2deg))
+    print("Rho interval ............................: %.2f-%.2f deg" % (rho_min*gammalib.rad2deg,rho_max*gammalib.rad2deg))
+
+    # Plot rhos
+    for rho in rhos:
+
+        # Print header
+        if log:
+            print("Rho ..................................: %.2f deg" % (rho))
+
+        # Compute some stuff
+        rho_rad = rho   * gammalib.deg2rad
+
+        # Compute half length of arc that lies within Psf validity circle
+        domega = 0.5 * cta_roi_arclength(rho_rad, rho_obs, delta_max)
+
+        # Print domega
+        if log:
+            print("  domega .............................: %.2f deg" % (domega*gammalib.rad2deg))
+
+        # Continue only if arc length is positive
+        if (domega > 0.0):
+
+            # Precompute cosine and sine terms for azimuthal integration
+            cos_rho = math.cos(rho_rad);
+            sin_rho = math.sin(rho_rad);
+            cos_psf = cos_rho * math.cos(rho_obs)
+            sin_psf = sin_rho * math.sin(rho_obs)
+            cos_ph  = cos_rho * math.cos(rho_pnt)
+            sin_ph  = sin_rho * math.sin(rho_pnt)
+
+            # Full containment within ellipse
+            if (rho < bound_minor):
+
+                # Set Omega range
+                omega_min = -domega
+                omega_max = +domega
+
+                # Print
+                if log:
+                    print("  Omega range ........................: %.2f - %.2f deg" % (omega_min*gammalib.rad2deg,omega_max*gammalib.rad2deg))
+
+                # Plot segment
+                n    = 100
+                x    = []
+                y    = []
+                step = (omega_max - omega_min)/n
+                for i in range(n+1):
+                    radians = i*step + omega_min
+                    cos_radians = math.sin(radians+posangle_obs)
+                    sin_radians = math.cos(radians+posangle_obs)
+                    x.append(rho * cos_radians)
+                    y.append(rho * sin_radians)
+                plt.plot(x, y, 'k-')
+
+            # Partial containment within ellipse
+            else:
+
+                # Azimuth angle of intersection points between ellipse and
+                # circle with radius of rho
+                arg1        = 1.0 - (bound_minor*bound_minor) / (rho*rho)
+                arg2        = 1.0 - (bound_minor*bound_minor) / (bound_major*bound_major)
+                omega_width = math.acos(math.sqrt(arg1/arg2))
+
+                # Print Omega width
+                if log:
+                    print("  omega_width ........................: %.2f deg" % (omega_width*gammalib.rad2deg))
+
+                # 
+                if (omega_width > 0.0):
+                
+                    # Compute Omega intervals. Recall that Omega is with respect
+                    # to connecting line between model and observed photon. Make
+                    # sure that omega_0 is within [-pi,pi] thus that the omega
+                    # intervals are between [-2pi,2pi]
+                    omega_0 = bound_pa * gammalib.deg2rad - posangle_obs
+                    if (omega_0 > gammalib.pi):
+                        omega_0 -= gammalib.pi
+                    elif (omega_0 < -gammalib.pi):
+                        omega_0 += gammalib.pi
+
+                    # Compute intervals
+                    omega1_min = omega_0    - omega_width
+                    omega1_max = omega_0    + omega_width
+                    omega2_min = omega1_min + gammalib.pi
+                    omega2_max = omega1_max + gammalib.pi
+
+                    # Logging
+                    if log:
+                        print("  omega_0 ............................: %.2f deg" % (omega_0*gammalib.rad2deg))
+                        print("  Omega range 1 ......................: %.2f - %.2f deg" % (omega1_min*gammalib.rad2deg,omega1_max*gammalib.rad2deg))
+                        print("  Omega range 2 ......................: %.2f - %.2f deg" % (omega2_min*gammalib.rad2deg,omega2_max*gammalib.rad2deg))
+
+                    # Limit intervals
+                    intervals  = limit_omega(omega1_min, omega1_max, domega)
+                    intervals2 = limit_omega(omega2_min, omega2_max, domega)
+                    intervals.extend(intervals2)
+
+                    # Print results
+                    if log:
+                        for interval in intervals:
+                            print("  Omega interval .....................: %.2f - %.2f deg" % (interval['min']*gammalib.rad2deg,interval['max']*gammalib.rad2deg))
+
+                    # Draw Omega range lines and circle segments
+                    for interval in intervals:
+                        omega_min = interval['min']
+                        omega_max = interval['max']
+                        if (omega_min < omega_max):
+                            n    = 100
+                            x    = []
+                            y    = []
+                            step = (omega_max - omega_min)/n
+                            for i in range(n+1):
+                                radians = i*step + omega_min
+                                cos_radians = math.sin(radians+posangle_obs)
+                                sin_radians = math.cos(radians+posangle_obs)
+                                x.append(rho * cos_radians)
+                                y.append(rho * sin_radians)
+                            plt.plot(x, y, 'k-')
+
+
+    # Show plot
+    plt.show()
+
+    # Return
+    return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    Test elliptical response computation
+    """
+    # Dump header
+    print("")
+    print("****************************************")
+    print("* Test elliptical response computation *")
+    print("****************************************")
+
+    # Check if matplotlib is available
+    try:
+        import matplotlib.pyplot as plt
+    except ImportError:
+        print("Matplotlib is not installed on your system. Abort.")
+        sys.exit()
+
+    # Set Psf size
+    delta_max = 0.5
+
+    # Set model parameters
+    centre = gammalib.GSkyDir()
+    centre.radec_deg(0.0, 0.0)
+    semimajor = 1.0
+    semiminor = 0.50
+    #semimajor = 0.5
+    #semiminor = 1.0
+    posang    = 30.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.4)
+
+    # Set pointing direction
+    pntDir = gammalib.GSkyDir()
+    pntDir.radec_deg(0.1, -1.0)
+
+    # Show elliptical response geometry
+    show_response(centre, semimajor, semiminor, posang, obsDir, pntDir, delta_max)
+    
+#! /usr/bin/env python
+# ==========================================================================
+# This script tests the radial response computation. The script requires
+# matplotlib to be installed.
+#
+# --------------------------------------------------------------------------
+#
+# Copyright (C) 2014 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import math
+
+
+# ================= #
+# Compute arclength #
+# ================= #
+def cta_roi_arclength(rad, dist, roi):
+    """
+    Returns length of circular arc within circular ROI
+
+    Parameters:
+     rad  - Circle radius in radians (<pi)
+     dist - Circle centre distance to ROI centre (<pi)
+     roi  - Radius of ROI in radians
+    """
+    # Initialise
+    arclength = 0.0
+
+    # Handle special case that circle centre matches ROI centre
+    if dist == 0.0:
+        if rad > roi:
+            arclength = 0.0
+        else:
+            arclength = 2.0 * math.pi
+
+    # ... otherwise circle and ROI centres are not identical
+    else:
+
+        # Handle special case that we evaluate exactly at the circle
+        # centre. In this case we have in fact a point, and if this point
+        # falls within the ROI it has a formal arclength of 2pi.
+        if rad == 0.0:
+            if dist > roi:
+                arclength = 0.0
+            else:
+                arclength = 2.0 * math.pi
+
+        # ... otherwise we have to handle the general case
+        else:
+            d = roi - dist
+            if -rad >= d:
+                arclength = 0.0
+            elif rad <= d:
+                arclength = 2.0 * math.pi
+            else:
+                cosrad  = math.cos(rad)
+                sinrad  = math.sin(rad)
+                cosroi  = math.cos(roi)
+                cosdist = math.cos(dist)
+                sindist = math.sin(dist)
+                cosang  = (cosroi - cosdist * cosrad) / (sindist * sinrad)
+                arclength = 2.0 * gammalib.acos(cosang)
+
+    # Return
+    return arclength
+    
+
+# ==================== #
+# Show radial response #
+# ==================== #
+def show_response(centre, radius, obsDir, delta_max):
+    """
+    Show radial response
+    """
+    # Create figure
+    plt.figure(1)
+    plt.title("Radial response")
+
+    # Compute some stuff
+    rho_obs      = centre.dist_deg(obsDir)
+    posangle_obs = centre.posang(obsDir)
+
+    # Rho limits
+    rho_min = 0.0
+    rho_max = rho_obs + delta_max
+    if (rho_obs > delta_max):
+        rho_min = rho_obs - delta_max
+    if rho_max > radius:
+        rho_max = radius
+
+    # Convert to radians
+    delta_max *= gammalib.deg2rad
+    rho_obs   *= gammalib.deg2rad
+
+    # Setup rhos
+    rhos = []
+    for i in range(51):
+        rho = i*0.02
+        if (rho >= rho_min and rho <= rho_max):
+            rhos.append(rho)
+
+    # Plot radial model
+    x = []
+    y = []
+    for angle in range(0, 361):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(radius * cos_radians)
+        y.append(radius * sin_radians)
+    plt.plot(x, y, 'r-')
+    plt.xlim([ 2.0*radius,-2.0*radius])
+    plt.ylim([-2.0*radius, 2.0*radius])
+    plt.plot([0.0], [0.0], 'ro')
+
+    # Plot observed direction
+    x_obs = obsDir.ra_deg()  - centre.ra_deg()
+    y_obs = obsDir.dec_deg() - centre.dec_deg()
+    plt.plot([x_obs], [y_obs], 'bo')
+
+    # Draw connecting line model-photon'
+    n    = 10
+    step = rho_obs/n * gammalib.rad2deg
+    x    = []
+    y    = []
+    for i in range(n+1):
+        x.append(i*step * math.sin(posangle_obs))
+        y.append(i*step * math.cos(posangle_obs))
+    plt.plot(x, y, 'b-')
+
+    # Plot Psf validity circle
+    x = []
+    y = []
+    for angle in range(0, 361):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(delta_max * cos_radians * gammalib.rad2deg + x_obs)
+        y.append(delta_max * sin_radians * gammalib.rad2deg + y_obs)
+    plt.plot(x, y, 'b-')
+
+    # Plot rhos
+    for rho in rhos:
+
+        # Compute some stuff
+        rho_rad = rho * gammalib.deg2rad
+
+        # Compute half length of arc that lies within Psf validity circle
+        domega = 0.5 * cta_roi_arclength(rho_rad, rho_obs, delta_max)
+
+        # Continue only if arc length is positive
+        if (domega > 0.0):
+
+            # Set Omega range
+            omega_min = -domega
+            omega_max = +domega
+
+            # Plot segment
+            n    = 100
+            x    = []
+            y    = []
+            step = (omega_max - omega_min)/n
+            for i in range(n+1):
+                radians = i*step + omega_min
+                cos_radians = math.sin(radians+posangle_obs)
+                sin_radians = math.cos(radians+posangle_obs)
+                x.append(rho * cos_radians)
+                y.append(rho * sin_radians)
+            plt.plot(x, y, 'k-')
+
+    # Show plot
+    plt.show()
+
+    # Return
+    return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    Test radial response computation
+    """
+    # Dump header
+    print("")
+    print("************************************")
+    print("* Test radial response computation *")
+    print("************************************")
+
+    # Check if matplotlib is available
+    try:
+        import matplotlib.pyplot as plt
+    except ImportError:
+        print("Matplotlib is not installed on your system. Abort.")
+        sys.exit()
+
+    # Set Psf size
+    delta_max = 0.8
+
+    # Set model parameters
+    centre = gammalib.GSkyDir()
+    centre.radec_deg(0.0, 0.0)
+    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.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)
+#! /usr/bin/env python
+# ==========================================================================
+# This script tests the elliptical Npred computation. The script requires
+# matplotlib to be installed.
+#
+# --------------------------------------------------------------------------
+#
+# Copyright (C) 2014 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import math
+
+
+# ================= #
+# Compute arclength #
+# ================= #
+def cta_roi_arclength(rad, dist, roi):
+    """
+    Returns length of circular arc within circular ROI
+
+    Parameters:
+     rad  - Circle radius in radians (<pi)
+     dist - Circle centre distance to ROI centre (<pi)
+     roi  - Radius of ROI in radians
+    """
+    # Initialise
+    arclength = 0.0
+
+    # Handle special case that circle centre matches ROI centre
+    if dist == 0.0:
+        if rad > roi:
+            arclength = 0.0
+        else:
+            arclength = 2.0 * math.pi
+
+    # ... otherwise circle and ROI centres are not identical
+    else:
+
+        # Handle special case that we evaluate exactly at the circle
+        # centre. In this case we have in fact a point, and if this point
+        # falls within the ROI it has a formal arclength of 2pi.
+        if rad == 0.0:
+            if dist > roi:
+                arclength = 0.0
+            else:
+                arclength = 2.0 * math.pi
+
+        # ... otherwise we have to handle the general case
+        else:
+            d = roi - dist
+            if -rad >= d:
+                arclength = 0.0
+            elif rad <= d:
+                arclength = 2.0 * math.pi
+            else:
+                cosrad  = math.cos(rad)
+                sinrad  = math.sin(rad)
+                cosroi  = math.cos(roi)
+                cosdist = math.cos(dist)
+                sindist = math.sin(dist)
+                cosang  = (cosroi - cosdist * cosrad) / (sindist * sinrad)
+                arclength = 2.0 * gammalib.acos(cosang)
+
+    # Return
+    return arclength
+
+
+# ============== #
+# Limit interval #
+# ============== #
+def limit_omega(min, max, domega, log=True):
+    """
+    """
+    # Compute only if domega < pi
+    if (domega < gammalib.pi):
+
+        # Save initial interval
+        min_initial = min
+        max_initial = max
+
+        # Initialise intervals
+        intervals = []
+        
+        # Limit each of the intervals to [omega_min, omega_max]
+        omega_min       = -domega
+        omega_max       = +domega
+        omega_min_plus  = omega_min + gammalib.twopi
+        omega_max_plus  = omega_max + gammalib.twopi
+        omega_min_minus = omega_min - gammalib.twopi
+        omega_max_minus = omega_max - gammalib.twopi
+
+        # Check for overlaps
+        in_unshifted = False
+        in_plus      = False
+        in_minus     = False
+        if (max_initial > omega_min and min_initial < omega_max):
+            min = min_initial
+            max = max_initial
+            in_unshifted = True
+            if (min < omega_min):
+                min = omega_min
+            if (max > omega_max):
+                max = omega_max
+            intervals.append({'min': min, 'max': max})
+        if (max_initial > omega_min_plus and min_initial < omega_max_plus):
+            min = min_initial
+            max = max_initial
+            in_plus = True
+            if (min < omega_min_plus):
+                min = omega_min_plus
+            if (max > omega_max_plus):
+                max = omega_max_plus
+            intervals.append({'min': min, 'max': max})
+        if (max_initial > omega_min_minus and min_initial < omega_max_minus):
+            min = min_initial
+            max = max_initial
+            in_minus = True
+            if (min < omega_min_minus):
+                min = omega_min_minus
+            if (max > omega_max_minus):
+                max = omega_max_minus
+            intervals.append({'min': min, 'max': max})
+
+        # Debug
+        if log:
+            print("[%7.2f,%7.2f] => [%7.2f,%7.2f] (%i[%7.2f,%7.2f]  %i[%7.2f,%7.2f]  %i[%7.2f,%7.2f])" %
+              (min_initial*gammalib.rad2deg, max_initial*gammalib.rad2deg, \
+               min*gammalib.rad2deg, max*gammalib.rad2deg, \
+               in_unshifted, \
+               omega_min*gammalib.rad2deg, omega_max*gammalib.rad2deg, \
+               in_plus, \
+               omega_min_plus*gammalib.rad2deg, omega_max_plus*gammalib.rad2deg, \
+               in_minus, \
+               omega_min_minus*gammalib.rad2deg, omega_max_minus*gammalib.rad2deg))
+
+    # ... just copy over the input interval
+    else:
+        intervals = [{'min': min, 'max': max}]
+
+    # Return
+    return intervals
+
+
+# ======================= #
+# Show ellptical response #
+# ======================= #
+def show_response(centre, semimajor, semiminor, posang, roiDir, pntDir, roi_radius, psf_max_radius, log=False):
+    """
+    Show ellptical response
+    """
+    # Create figure
+    plt.figure(1)
+    plt.title("Elliptical Npred computation")
+
+    # Get ellipse boundary (used for intersection computation)
+    if (semimajor >= semiminor):
+        bound_major  = semimajor
+        bound_minor  = semiminor
+        bound_pa     = posang
+    else:
+        bound_major  = semiminor
+        bound_minor  = semimajor
+        bound_pa     = posang + 90.0
+
+    # Compute effective ROI radius
+    radius_roi = roi_radius + psf_max_radius
+
+    # Compute some stuff
+    cos_pa       = math.cos(posang * gammalib.deg2rad)
+    sin_pa       = math.sin(posang * gammalib.deg2rad)
+    rho_roi      = roiDir.dist_deg(centre)
+    theta_max    = bound_major
+    posangle_roi = centre.posang(roiDir) # omega0
+
+    # Rho limits
+    rho_min = 0.0
+    rho_max = rho_roi + radius_roi
+    if (rho_roi > radius_roi):
+        rho_min = rho_roi - radius_roi
+    if (rho_max > theta_max):
+        rho_max = theta_max
+
+    # Convert to radians
+    roi_radius     *= gammalib.deg2rad
+    psf_max_radius *= gammalib.deg2rad
+    radius_roi     *= gammalib.deg2rad
+    rho_roi        *= gammalib.deg2rad
+    theta_max      *= gammalib.deg2rad
+
+    # Setup rhos
+    rhos = []
+    for i in range(51):
+        rho = i*0.02
+        if (rho >= rho_min and rho <= rho_max):
+            rhos.append(rho)
+
+    # Setup ellipse as vector
+    x = []
+    y = []
+    for angle in range(0, 360):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(semimajor * cos_radians * sin_pa - semiminor * sin_radians * cos_pa)
+        y.append(semimajor * cos_radians * cos_pa + semiminor * sin_radians * sin_pa)
+
+    # Plot ellipse
+    plt.plot(x, y, 'r-')
+    plt.xlim([ 2.0*bound_major,-2.0*bound_major])
+    plt.ylim([-2.0*bound_major, 2.0*bound_major])
+    plt.plot([0.0], [0.0], 'ro')
+
+    # Plot ROI centre
+    x_roi = roiDir.ra_deg()  - centre.ra_deg()
+    y_roi = roiDir.dec_deg() - centre.dec_deg()
+    plt.plot([x_roi], [y_roi], 'bo')
+
+    # Plot ROI circle
+    x = []
+    y = []
+    for angle in range(0, 360):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(roi_radius * cos_radians * gammalib.rad2deg + x_roi)
+        y.append(roi_radius * sin_radians * gammalib.rad2deg + y_roi)
+    plt.plot(x, y, 'b-.')
+
+    # Plot ROI+Psf circle
+    x = []
+    y = []
+    for angle in range(0, 360):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(radius_roi * cos_radians * gammalib.rad2deg + x_roi)
+        y.append(radius_roi * sin_radians * gammalib.rad2deg + y_roi)
+    plt.plot(x, y, 'b-')
+
+    # Print some stuff
+    print("Distance model-ROI (rho_roi) ............: %.2f deg" % (rho_roi*gammalib.rad2deg))
+    print("Rho interval ............................: %.2f-%.2f deg" % (rho_min*gammalib.rad2deg,rho_max*gammalib.rad2deg))
+
+    # Continue only if rho range is valid
+    if (rho_max > rho_min):
+
+        # Plot rhos
+        for rho in rhos:
+
+            # Print header
+            if log:
+                print("Rho ..................................: %.2f deg" % (rho))
+
+            # Compute some stuff
+            rho_rad = rho * gammalib.deg2rad
+
+            # Compute half length of arc that lies within Psf validity circle
+            domega = 0.5 * cta_roi_arclength(rho_rad, rho_roi, radius_roi)
+
+            # Print domega
+            if log:
+                print("  domega .............................: %.2f deg" % (domega*gammalib.rad2deg))
+
+            #
+            if (domega > 0.0):
+
+                # Full containment within ellipse
+                if (rho <= bound_minor):
+            
+                    # Compute phi integration range
+                    omega_min = -domega
+                    omega_max = +domega
+
+                    # Print
+                    if log:
+                        print("  Omega range ........................: %.2f - %.2f deg" % (omega_min*gammalib.rad2deg,omega_max*gammalib.rad2deg))
+
+                    # Plot segment
+                    n    = 100
+                    x    = []
+                    y    = []
+                    step = (omega_max - omega_min)/n
+                    for i in range(n+1):
+                        radians = i*step + omega_min
+                        cos_radians = math.sin(radians+posangle_roi)
+                        sin_radians = math.cos(radians+posangle_roi)
+                        x.append(rho * cos_radians)
+                        y.append(rho * sin_radians)
+                    plt.plot(x, y, 'k-')
+
+                # ... otherwise compute the arcs
+                else:
+
+                    # Azimuth angle of intersection points between ellipse and
+                    # circle with radius of rho
+                    arg1        = 1.0 - (bound_minor*bound_minor) / (rho*rho)
+                    arg2        = 1.0 - (bound_minor*bound_minor) / (bound_major*bound_major)
+                    omega_width = math.acos(math.sqrt(arg1/arg2))
+
+                    # Print Omega width
+                    if log:
+                        print("  omega_width ........................: %.2f deg" % (omega_width*gammalib.rad2deg))
+
+                    # 
+                    if (omega_width > 0.0):
+                
+                        # Compute Omega intervals. Recall that Omega is with respect
+                        # to connecting line between model and observed photon. Make
+                        # sure that omega_0 is within [-pi,pi] thus that the omega
+                        # intervals are between [-2pi,2pi]
+                        omega_0 = bound_pa * gammalib.deg2rad - posangle_roi
+                        if (omega_0 > gammalib.pi):
+                            omega_0 -= gammalib.pi
+                        elif (omega_0 < -gammalib.pi):
+                            omega_0 += gammalib.pi
+
+                        # Compute intervals
+                        omega1_min = omega_0    - omega_width
+                        omega1_max = omega_0    + omega_width
+                        omega2_min = omega1_min + gammalib.pi
+                        omega2_max = omega1_max + gammalib.pi
+
+                        # Logging
+                        if log:
+                            print("  Omega range 1 (initial) ............: %.2f - %.2f deg" % (omega1_min*gammalib.rad2deg,omega1_max*gammalib.rad2deg))
+                            print("  Omega range 2 (initial) ............: %.2f - %.2f deg" % (omega2_min*gammalib.rad2deg,omega2_max*gammalib.rad2deg))
+
+                        # Limit intervals
+                        intervals  = limit_omega(omega1_min, omega1_max, domega)
+                        intervals2 = limit_omega(omega2_min, omega2_max, domega)
+                        intervals.extend(intervals2)
+
+                        # Print results
+                        if log:
+                            for interval in intervals:
+                                print("  Omega interval .....................: %.2f - %.2f deg" % (interval['min']*gammalib.rad2deg,interval['max']*gammalib.rad2deg))
+
+                        # Draw Omega range lines and circle segments
+                        for interval in intervals:
+                            omega_min = interval['min']
+                            omega_max = interval['max']
+                            if (omega_min < omega_max):
+                                n    = 100
+                                x    = []
+                                y    = []
+                                step = (omega_max - omega_min)/n
+                                for i in range(n+1):
+                                    radians = i*step + omega_min
+                                    cos_radians = math.sin(radians+posangle_roi)
+                                    sin_radians = math.cos(radians+posangle_roi)
+                                    x.append(rho * cos_radians)
+                                    y.append(rho * sin_radians)
+                                plt.plot(x, y, 'k-')
+
+    # Show plot
+    plt.show()
+
+    # Return
+    return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    Test elliptical response computation
+    """
+    # Dump header
+    print("")
+    print("*************************************")
+    print("* Test elliptical Npred computation *")
+    print("*************************************")
+
+    # Check if matplotlib is available
+    try:
+        import matplotlib.pyplot as plt
+    except ImportError:
+        print("Matplotlib is not installed on your system. Abort.")
+        sys.exit()
+
+    # Set Parameters
+    roi_radius = 0.5
+    delta_max  = 0.4
+
+    # Set model parameters
+    centre = gammalib.GSkyDir()
+    centre.radec_deg(0.0, 0.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.0, 0.1)
+
+    # Set pointing direction
+    pntDir = gammalib.GSkyDir()
+    pntDir.radec_deg(0.1, -1.0)
+
+    # Show elliptical response geometry
+    show_response(centre, semimajor, semiminor, posang, roiDir, pntDir, roi_radius, delta_max)
+    
+#! /usr/bin/env python
+# ==========================================================================
+# This script tests the radial response computation. The script requires
+# matplotlib to be installed.
+#
+# --------------------------------------------------------------------------
+#
+# Copyright (C) 2014 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import math
+
+
+# ================= #
+# Compute arclength #
+# ================= #
+def cta_roi_arclength(rad, dist, roi):
+    """
+    Returns length of circular arc within circular ROI
+
+    Parameters:
+     rad  - Circle radius in radians (<pi)
+     dist - Circle centre distance to ROI centre (<pi)
+     roi  - Radius of ROI in radians
+    """
+    # Initialise
+    arclength = 0.0
+
+    # Handle special case that circle centre matches ROI centre
+    if dist == 0.0:
+        if rad > roi:
+            arclength = 0.0
+        else:
+            arclength = 2.0 * math.pi
+
+    # ... otherwise circle and ROI centres are not identical
+    else:
+
+        # Handle special case that we evaluate exactly at the circle
+        # centre. In this case we have in fact a point, and if this point
+        # falls within the ROI it has a formal arclength of 2pi.
+        if rad == 0.0:
+            if dist > roi:
+                arclength = 0.0
+            else:
+                arclength = 2.0 * math.pi
+
+        # ... otherwise we have to handle the general case
+        else:
+            d = roi - dist
+            if -rad >= d:
+                arclength = 0.0
+            elif rad <= d:
+                arclength = 2.0 * math.pi
+            else:
+                cosrad  = math.cos(rad)
+                sinrad  = math.sin(rad)
+                cosroi  = math.cos(roi)
+                cosdist = math.cos(dist)
+                sindist = math.sin(dist)
+                cosang  = (cosroi - cosdist * cosrad) / (sindist * sinrad)
+                arclength = 2.0 * gammalib.acos(cosang)
+
+    # Return
+    return arclength
+    
+
+# ==================== #
+# Show radial response #
+# ==================== #
+def show_response(centre, theta_max, obsDir, psf_max):
+    """
+    Show radial response
+    """
+    # Create figure
+    plt.figure(1)
+    plt.title("Radial response")
+
+    # Compute some stuff
+    delta_mod    = obsDir.dist_deg(centre)
+    posangle_mod = obsDir.posang(centre)
+
+    # Delta limits
+    delta_min = 0.0
+    delta_max = delta_mod + theta_max
+    if (delta_mod > theta_max):
+        delta_min = delta_mod - theta_max
+    if delta_max > psf_max:
+        delta_max = psf_max
+
+    # Convert to radians
+    delta_mod *= gammalib.deg2rad
+    theta_max *= gammalib.deg2rad
+
+    # Setup deltas
+    deltas = []
+    for i in range(51):
+        delta = i*0.02
+        if (delta >= delta_min and delta <= delta_max):
+            deltas.append(delta)
+
+    # Plot Psf validity circle
+    x = []
+    y = []
+    for angle in range(0, 361):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(psf_max * cos_radians)
+        y.append(psf_max * sin_radians)
+    plt.plot(x, y, 'b-')
+    plt.xlim([ 2.0*psf_max,-2.0*psf_max])
+    plt.ylim([-2.0*psf_max, 2.0*psf_max])
+    plt.plot([0.0], [0.0], 'bo')
+
+    # Plot model direction
+    x_mod = centre.ra_deg()  - obsDir.ra_deg()
+    y_mod = centre.dec_deg() - obsDir.dec_deg()
+    plt.plot([x_mod], [y_mod], 'ro')
+
+    # Plot radial model
+    x = []
+    y = []
+    for angle in range(0, 361):
+        radians     = angle * gammalib.deg2rad
+        cos_radians = math.cos(radians)
+        sin_radians = math.sin(radians)
+        x.append(theta_max * cos_radians * gammalib.rad2deg + x_mod)
+        y.append(theta_max * sin_radians * gammalib.rad2deg + y_mod)
+    plt.plot(x, y, 'r-')
+
+    # Plot deltas
+    for delta in deltas:
+
+        # Compute some stuff
+        delta_rad = delta * gammalib.deg2rad
+
+        # Compute half length of arc that lies within model
+        dphi = 0.5 * cta_roi_arclength(delta_rad, delta_mod, theta_max)
+
+        # Continue only if arc length is positive
+        if (dphi > 0.0):
+
+            # Set Phi range
+            phi_min = -dphi
+            phi_max = +dphi
+
+            # Plot segment
+            n    = 100
+            x    = []
+            y    = []
+            step = (phi_max - phi_min)/n
+            for i in range(n+1):
+                radians = i*step + phi_min
+                cos_radians = math.sin(radians+posangle_mod)
+                sin_radians = math.cos(radians+posangle_mod)
+                x.append(delta * cos_radians)
+                y.append(delta * sin_radians)
+            plt.plot(x, y, 'k-')
+
+    # Show plot
+    plt.show()
+
+    # Return
+    return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    Test radial response computation
+    """
+    # Dump header
+    print("")
+    print("*******************************")
+    print("* Test radial Psf computation *")
+    print("*******************************")
+
+    # Check if matplotlib is available
+    try:
+        import matplotlib.pyplot as plt
+    except ImportError:
+        print("Matplotlib is not installed on your system. Abort.")
+        sys.exit()
+
+    # Set Psf size
+    psf_max = 1.0
+
+    # Set observed photon direction
+    obsDir = gammalib.GSkyDir()
+    obsDir.radec_deg(0.0, 0.0)
+
+    # Set model parameters
+    centre = gammalib.GSkyDir()
+    centre.radec_deg(0.0, 0.6)
+    theta_max = 0.5
+
+    # Show elliptical response geometry
+    show_response(centre, theta_max, obsDir, psf_max)
+    
@@ -617,23 +617,8 @@ void TestGLATObservation::test_one_binned_obs(const std::string&amp; datadir, const 
 void TestGLATOptimize::test_binned_optimizer_p6(void)
 {
     // Set expected fit results
-    /*
     double fit_results[] = {1, 0,
-                            3.214411831, 0.744781,
-                            1, 0,
-                            1, 0,
-                            0.8068908818, 0.08764042542,
-                            1, 0,
-                            83.6331, 0,
-                            22.0145, 0,
-                            2.103749294e-06, 1.560611406e-07,
-                            -2.126826057, 0.053618,
-                            100, 0,
-                            500000, 0,
-                            1, 0};
-    */
-    double fit_results[] = {1, 0,
-                            2.342467831, 0.7339917,
+                            2.34236, 0.7339917,
                             1, 0,
                             1, 0,
                             0.8167256818, 0.08569407542,
@@ -663,23 +648,8 @@ void TestGLATOptimize::test_binned_optimizer_p6(void)
 void TestGLATOptimize::test_binned_optimizer_p7(void)
 {
     // Set expected fit results
-    /*
-    double fit_results[] = {1, 0,
-                            2.896021569, 0.441842,
-                            1, 0,
-                            1, 0,
-                            0.8629520722, 0.06376048926,
-                            1, 0,
-                            83.6331, 0,
-                            22.0145, 0,
-                            1.922525774e-06, 1.209507237e-07,
-                            -2.038629122, 0.0509173084,
-                            100, 0,
-                            500000, 0,
-                            1, 0};
-    */
     double fit_results[] = {1, 0,
-                            2.374909569, 0.4548979,
+                            2.37468, 0.4548979,
                             1, 0,
                             1, 0,
                             0.8419824722, 0.06340531926,
@@ -687,7 +657,7 @@ void TestGLATOptimize::test_binned_optimizer_p7(void)
                             83.6331, 0,
                             22.0145, 0,
                             1.922525774e-06, 1.209507237e-07,
-                            -2.124145722, 0.0493105484,
+                            -2.12421, 0.0493105484,
                             100, 0,
                             500000, 0,
                             1, 0};
 /***************************************************************************
  *               test_MWL.cpp - Test multi-wavelength classes              *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2013 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -219,13 +219,13 @@ void TestGMWL::test_optimizer(void)
     double fit_results[] = {83.6331,             0,
                             22.0145,             0,
                             0.00231438040883405, 0.000191759607698647,
-                            -3.1111154491587,    0.125598,
+                            -3.12006,            0.128382,
                             1,                   0,
                             1,                   0,
                             83.6331,             0,
                             22.0145,             0,
                             4.5050091460345e-11, 1.95885803229878e-12,
-                            -2.13740092737976,   0.00578641,
+                            -2.13765,            0.00578641,
                             1000,                0,
                             1,                   0};
     test_try("Perform LM optimization");
@@ -73,6 +73,10 @@ public:
     void    dec(const double& dec);
     double  posangle(void) const;
     void    posangle(const double& posangle);
+    double  semiminor(void) const;
+    double  semimajor(void) const;
+    void    semiminor(const double& semiminor);
+    void    semimajor(const double& semimajor);
     GSkyDir dir(void) const;
     void    dir(const GSkyDir& dir);
 };
 /***************************************************************************
  *    GModelSpatialEllipticalDisk.i - Elliptical disk source model class   *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013 by Michael Mayer                                    *
+ *  copyright (C) 2013-2014 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -65,12 +65,6 @@ public:
     virtual double                       theta_max(void) const;
     virtual void                         read(const GXmlElement& xml);
     virtual void                         write(GXmlElement& xml) const;
-
-    // Other methods
-    double semiminor(void) const;
-    double semimajor(void) const;
-    void   semiminor(const double& semiminor);
-    void   semimajor(const double& semimajor);
 };
 /***************************************************************************
  *  GModelSpatialElliptical.cpp - Abstract elliptical spatial model class  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013 by Juergen Knoedlseder                              *
+ *  copyright (C) 2013-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -467,11 +467,35 @@ void GModelSpatialElliptical::init_members(void)
     m_posangle.gradient(0.0);
     m_posangle.has_grad(false);
+    // Initialise semi-major axis
+    m_semimajor.clear();
+    m_semimajor.name("MajorRadius");
+    m_semimajor.unit("deg");
+    m_semimajor.value(2.778e-4); // 1 arcsec
+    m_semimajor.min(2.778e-4);   // 1 arcsec
+    m_semimajor.free();
+    m_semimajor.scale(1.0);
+    m_semimajor.gradient(0.0);
+    m_semimajor.has_grad(false); // Elliptical components never have gradients
+
+    // Initialise semi-minor axis
+    m_semiminor.clear();
+    m_semiminor.name("MinorRadius");
+    m_semiminor.unit("deg");
+    m_semiminor.value(2.778e-4); // 1 arcsec
+    m_semiminor.min(2.778e-4);   // 1 arcsec
+    m_semiminor.free();
+    m_semiminor.scale(1.0);
+    m_semiminor.gradient(0.0);
+    m_semiminor.has_grad(false); // Elliptical components never have gradients
+
     // Set parameter pointer(s)
     m_pars.clear();
     m_pars.push_back(&m_ra);
     m_pars.push_back(&m_dec);
     m_pars.push_back(&m_posangle);
+    m_pars.push_back(&m_semimajor);
+    m_pars.push_back(&m_semiminor);
     // Return
     return;
@@ -486,15 +510,19 @@ void GModelSpatialElliptical::init_members(void)
 void GModelSpatialElliptical::copy_members(const GModelSpatialElliptical& model)
 {
     // Copy members
-    m_ra       = model.m_ra;
-    m_dec      = model.m_dec;
-    m_posangle = model.m_posangle;
+    m_ra        = model.m_ra;
+    m_dec       = model.m_dec;
+    m_posangle  = model.m_posangle;
+	m_semiminor = model.m_semiminor;
+	m_semimajor = model.m_semimajor;
     // Set parameter pointer(s)
     m_pars.clear();
     m_pars.push_back(&m_ra);
     m_pars.push_back(&m_dec);
     m_pars.push_back(&m_posangle);
+    m_pars.push_back(&m_semimajor);
+    m_pars.push_back(&m_semiminor);
     // Return
     return;
 /***************************************************************************
  *   GModelSpatialEllipticalDisk.cpp - Elliptical disk source model class  *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2013 by Michael Mayer                                    *
+ *  copyright (C) 2013-2014 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -405,7 +405,7 @@ GSkyDir GModelSpatialEllipticalDisk::mc(const GEnergy&amp; energy,
 /***********************************************************************//**
  * @brief Return maximum model radius (in radians)
  *
- * @return Returns maximum model radians.
+ * @return Returns maximum model radius.
  ***************************************************************************/
 double GModelSpatialEllipticalDisk::theta_max(void) const
 {
@@ -634,32 +634,6 @@ std::string GModelSpatialEllipticalDisk::print(const GChatter&amp; chatter) const
  ***************************************************************************/
 void GModelSpatialEllipticalDisk::init_members(void)
 {
-    // Initialise semi-minor axis
-    m_semiminor.clear();
-    m_semiminor.name("MinorRadius");
-    m_semiminor.unit("deg");
-    m_semiminor.value(2.778e-4); // 1 arcsec
-    m_semiminor.min(2.778e-4);   // 1 arcsec
-    m_semiminor.free();
-    m_semiminor.scale(1.0);
-    m_semiminor.gradient(0.0);
-    m_semiminor.has_grad(false);  // Elliptical components never have gradients
-
-    // Initialise semi-major axis
-    m_semimajor.clear();
-    m_semimajor.name("MajorRadius");
-    m_semimajor.unit("deg");
-    m_semimajor.value(2.778e-4); // 1 arcsec
-    m_semimajor.min(2.778e-4);   // 1 arcsec
-    m_semimajor.free();
-    m_semimajor.scale(1.0);
-    m_semimajor.gradient(0.0);
-    m_semimajor.has_grad(false);  // Elliptical components never have gradients
-
-    // Set parameter pointer(s)
-    m_pars.push_back(&m_semiminor);
-    m_pars.push_back(&m_semimajor);
-
     // Initialise precomputation cache. Note that zero values flag
     // uninitialised as a zero radius is not meaningful
     m_last_semiminor = 0.0;
@@ -684,10 +658,6 @@ void GModelSpatialEllipticalDisk::init_members(void)
  ***************************************************************************/
 void GModelSpatialEllipticalDisk::copy_members(const GModelSpatialEllipticalDisk& model)
 {
-    // Copy members
-	m_semiminor = model.m_semiminor;
-	m_semimajor = model.m_semimajor;
-
     // Copy precomputation cache
     m_last_semiminor = model.m_last_semiminor;
     m_last_semimajor = model.m_last_semimajor;
@@ -722,7 +692,7 @@ void GModelSpatialEllipticalDisk::free_members(void)
  ***************************************************************************/
 void GModelSpatialEllipticalDisk::update() const
 {
-    // Update if one axis has changed has changed
+    // Update if one axis has changed
     if (m_last_semiminor != semiminor() || m_last_semimajor != semimajor()) {
         // Store last values
@@ -635,10 +635,10 @@ void GOptimizerLM::init_members(void)
     m_lambda_start = 1.0e-3;
     m_lambda_inc   = 10.0;
     m_lambda_dec   = 0.1;
-    m_eps          = 1.0e-6;
-    m_max_iter     = 1000;
+    m_eps          = 5.0e-3;  //!< Changed on 30/10/2014 from 1.0e-6
+    m_max_iter     = 100;     //!< Changed on 30/10/2014 from 1000
     m_max_stall    = 10;
-    m_max_hit      = 3; //!< Maximum successive boundary hits before freeze
+    m_max_hit      = 3;       //!< Maximum successive boundary hits before freeze
     m_step_adjust  = true;
     // Initialise bookkeeping arrays