Merge branch '2694-hess-background-models' into integration

Jürgen Knödlseder
2 parents b86b2e29 5ed3c029
Showing 24 changed files with 734 additions and 374 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.6.rst
include/GResponse.hpp
inst/cta/include/GCTAEdisp2D.hpp
inst/cta/include/GCTAEdispPerfTable.hpp
inst/cta/include/GCTAEdispRmf.hpp
inst/cta/include/GCTAModelSpatialGradient.hpp
inst/cta/include/GCTAModelSpatialMultiplicative.hpp
inst/cta/pyext/GCTAModelSpatialGradient.i
inst/cta/pyext/GCTAModelSpatialMultiplicative.i
inst/cta/pyext/cta.i
inst/cta/src/GCTACubeBackground.cpp
inst/cta/src/GCTAEdisp2D.cpp
inst/cta/src/GCTAEdispPerfTable.cpp
inst/cta/src/GCTAEdispRmf.cpp
inst/cta/src/GCTAModelSpatialGradient.cpp
inst/cta/src/GCTAResponseIrf.cpp
inst/cta/src/GCTAResponse_helpers.cpp
inst/cta/test/test_CTA.cpp
-2018-10-17
+2018-10-25
  
         * Version 1.6.0 released
         ========================
  
+        GCTAEdisp returns now energy dispersion values per MeV (#2704) 
         Protect GCTAInstDir from access to invalid information (#2694)
         Implement general CTA background model (#2694)
         Fix bug that potentially wrote the wrong exposure time in a PHA file (#2683)
 New Features and Important Changes in GammaLib 1.6.0
  
-17 October 2018
+25 October 2018
  
  
 1. Introduction
@@ -300,6 +300,12 @@ None
  
 20. CTA interface
 -----------------
+The GCTAEdisp::operator() now returns energy dispersion values per MeV. The
+units are now consistent with the expectation in GResponse, and this changes
+assurance the coherence between the various classes. In particular, no re-
+normalization is needed anymore to digest the H.E.S.S. energy dispersion, and
+also the CTA energy dispersion is handled correctly (#2704).
+
 The GCTAInstDir class now assures that only valid information can be accessed.
 Since the class can hold a sky direction and detector coordinates, flags were
 added that signal what information actually exists in the class. If a non-
@@ -45,6 +45,8 @@ Bug fixes
 Improvements
 ------------
  
+* [`2704 <https://cta-redmine.irap.omp.eu/issues/2704>`_] -
+  ``GCTAEdisp`` returns now energy dispersion values per MeV
 * [`2694 <https://cta-redmine.irap.omp.eu/issues/2694>`_] -
   Protect ``GCTAInstDir`` from access to invalid information
 * [`2677 <https://cta-redmine.irap.omp.eu/issues/2677>`_] -
 /***************************************************************************
  *               GResponse.hpp - Abstract response base class              *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2008-2015 by Juergen Knoedlseder                         *
+ *  copyright (C) 2008-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -127,7 +127,7 @@ protected:
                    m_event(event),
                    m_srcTime(srcTime),
                    m_grad(grad) { }
-        double eval(const double& x);
+        double eval(const double& logEsrc);
     protected:
         const GResponse*    m_parent;  //!< Response
         const GObservation* m_obs;     //!< Observation
@@ -52,14 +52,47 @@ namespace gammalib {
  * @brief CTA 2D energy dispersion class
  *
  * This class implements the energy dispersion for the CTA 2D response. The
- * energy dispersion is defined as
+ * CTA 2D energy dispersion is in fact a 3-dimensional function, where the
+ * energy dispersion is given as function of true energy \f$E_{\rm true}\f$,
+ * migration value \f$m = E_{\rm reco}/E_{\rm true}\f$ and offset angle
+ * \f$\theta\f$
  *
  * \f[
- *    E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}, \theta) =
- *    \frac{dP}{d\log_{10} E_{\rm reco}}
+ *    E_{\rm disp}(E_{\rm true}, m, \theta)
  * \f]
  *
- * and given in units of \f$(\log_{10} MeV)^{-1}\f$.
+ * with
+ *
+ * \f[
+ *    \int_{m_{\rm min}}^{m_{\rm max}}
+ *    E_{\rm disp}(E_{\rm true}, m, \theta) \, dm = 1
+ * \f]
+ *
+ * The energy dispersion as function of reconstructed energy
+ * \f$E_{\rm reco}\f$ is given by
+ *
+ * \f[
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}, \theta) =
+ *    \frac{E_{\rm disp}(E_{\rm true}, m, \theta)}{E_{\rm true}}
+ * \f]
+ *
+ * the energy dispersion as function of the natural logarithm of
+ * reconstructed energy \f$\log E_{\rm reco}\f$ is given by
+ *
+ * \f[
+ *    E_{\rm disp}(E_{\rm true}, \log E_{\rm reco}, \theta) =
+ *    E_{\rm disp}(E_{\rm true}, m, \theta) \times
+ *    \frac{E_{\rm reco}}{E_{\rm true}}
+ * \f]
+ *
+ * and the energy dispersion as function of the base 10 logarithm of
+ * reconstructed energy \f$\log_{10} E_{\rm reco}\f$ is given by
+ *
+ * \f[
+ *    E_{\rm disp}(E_{\rm true}, \log_{10} E_{\rm reco}, \theta) =
+ *    E_{\rm disp}(E_{\rm true}, m, \theta) \times
+ *    \frac{\log 10 \times E_{\rm reco}}{E_{\rm true}}
+ * \f]
  ***************************************************************************/
 class GCTAEdisp2D : public GCTAEdisp {
  
@@ -187,7 +220,7 @@ private:
                    m_parent(parent),
                    m_logEsrc(logEsrc),
                    m_theta(theta) { }
-        double eval(const double& x);
+        double eval(const double& logEobs);
     protected:
         const GCTAEdisp2D* m_parent;  //!< Pointer to parent class
         double             m_logEsrc; //!< True photon energy
@@ -44,18 +44,25 @@ class GRan;
  * @brief CTA performance table energy dispersion class
  *
  * This class implements the CTA energy dispersion response as function
- * of energy as determined from a performance table. The performance table is
- * an ASCII file that specifies the CTA performance parameters in a simple
- * way.
+ * of true energy as determined from a performance table. A performance
+ * table is an ASCII file that represent a table, specifying the CTA
+ * performance parameters as function of true photon energy.
  *
  * The energy dispersion is defined as
  *
  * \f[
- *    E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}) =
- *    \frac{dP}{d\log_{10} E_{\rm reco}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}) =
+ *    \frac{1}{\sqrt{2\pi}\sigma(E_{\rm true})}
+ *    \exp \left(\frac{-(\log_{10} E_{\rm reco} - \log_{10} E_{\rm true})^2}
+ *                    {2 \sigma(E_{\rm true})^2} \right) \times
+ *    \frac{1}{\log_{10} E_{\rm reco}}
  * \f]
  *
- * and given in units of \f$(\log_{10} MeV)^{-1}\f$.
+ * and given in units of MeV\f$^{-1}\f$, where
+ * \f$E_{\rm reco}\f$ is the reconstructed energy in units of MeV,
+ * \f$E_{\rm true}\f$ is the true energy in units of MeV, and
+ * \f$\sigma(E_{\rm true})\f$ is the standard deviation of the energy
+ * dispersion that depends on the true photon energy.
  ***************************************************************************/
 class GCTAEdispPerfTable : public GCTAEdisp {
  
@@ -50,11 +50,14 @@ class GEbounds;
  * The energy dispersion is defined as
  *
  * \f[
- *    E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}) =
- *    \frac{dP}{d\log_{10} E_{\rm reco}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}) =
+ *    \frac{E_{\rm disp}(E_{\rm true}, \log_{10} E_{\rm reco}, \theta)}
+ *         {\log_{10} E_{\rm reco}}
  * \f]
  *
- * and given in units of \f$(\log_{10} MeV)^{-1}\f$.
+ * in units of MeV\f$^{-1}\f$ where
+ * \f$E_{\rm reco}\f$ is the reconstructed energy in units of MeV, and
+ * \f$E_{\rm true}\f$ is the true energy in units of MeV.
  ***************************************************************************/
 class GCTAEdispRmf : public GCTAEdisp {
  
@@ -113,7 +116,7 @@ private:
     void set_matrix(void);
     void set_cache(void) const;
     void set_max_edisp(void) const;
-    void update(const double& arg1, const double& arg2) const;
+    void update(const double& logEsrc, const double& logEobs) const;
     void compute_ebounds_obs(const double& theta = 0.0,
                              const double& phi = 0.0,
                              const double& zenith = 0.0,
@@ -50,6 +50,8 @@ class GCTAModelSpatialGradient  : public GCTAModelSpatial {
 public:
     // Constructors and destructors
     GCTAModelSpatialGradient(void);
+    GCTAModelSpatialGradient(const double& detx_gradient,
+                             const double& dety_gradient);
     explicit GCTAModelSpatialGradient(const GXmlElement& xml);
     GCTAModelSpatialGradient(const GCTAModelSpatialGradient& model);
     virtual ~GCTAModelSpatialGradient(void);
@@ -30,11 +30,6 @@
 /* __ Includes ___________________________________________________________ */
 #include <string>
 #include <vector>
-//#include "GBase.hpp"
-//#include "GFunction.hpp"
-//#include "GEnergy.hpp"
-//#include "GTime.hpp"
-//#include "GModelPar.hpp"
 #include "GCTAModelSpatial.hpp"
  
 /* __ Forward declarations _______________________________________________ */
+/***************************************************************************
+ *       GCTAModelSpatialGradient.i - Spatial gradient CTA model class     *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2018 by 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/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GCTAModelSpatialGradient.i
+ * @brief Spatial gradient CTA interface definition
+ * @author Juergen Knoedlseder
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "GCTAModelSpatialGradient.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class GCTAModelSpatialGradient
+ *
+ * @brief Spatial gradient CTA model class
+ ***************************************************************************/
+class GCTAModelSpatialGradient  : public GCTAModelSpatial {
+
+public:
+    // Constructors and destructors
+    GCTAModelSpatialGradient(void);
+    GCTAModelSpatialGradient(const double& detx_gradient,
+                             const double& dety_gradient);
+    explicit GCTAModelSpatialGradient(const GXmlElement& xml);
+    GCTAModelSpatialGradient(const GCTAModelSpatialGradient& model);
+    virtual ~GCTAModelSpatialGradient(void);
+
+    // Operators
+    virtual GCTAModelSpatialGradient& operator=(const GCTAModelSpatialGradient& model);
+
+    // Implemented pure virtual methods
+    virtual void                      clear(void);
+    virtual GCTAModelSpatialGradient* clone(void) const;
+    virtual std::string               classname(void) const;
+    virtual std::string               type(void) const;
+    virtual double                    eval(const GCTAInstDir& dir,
+                                           const GEnergy&     energy,
+                                           const GTime&       time,
+                                           const bool&        gradients = false) const;
+    virtual GCTAInstDir               mc(const GEnergy& energy,
+                                         const GTime&   time,
+                                         GRan& ran) const;
+    virtual void                      read(const GXmlElement& xml);
+    virtual void                      write(GXmlElement& xml) const;
+
+    // Other methods
+    double detx_gradient(void) const;
+    double dety_gradient(void) const;
+    void   detx_gradient(const double& detx_gradient);
+    void   dety_gradient(const double& dety_gradient);
+};
+
+
+/***********************************************************************//**
+ * @brief GCTAModelSpatialGradient class extension
+ ***************************************************************************/
+%extend GCTAModelSpatialGradient {
+    GCTAModelSpatialGradient copy() {
+        return (*self);
+    }
+%pythoncode {
+    def __getstate__(self):
+        xml = gammalib.GXmlElement()
+        self.write(xml)
+        state = (xml,)
+        return state
+    def __setstate__(self, state):
+        if state[0].elements('parameter') == 0:
+            self.__init__()
+        else:
+            self.__init__(state[0])
+}
+};
+/***************************************************************************
+ *  GCTAModelSpatialMultiplicative.i - Multiplicative spatial model class  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2018 by 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/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file GCTAModelSpatialMultiplicative.i
+ * @brief Multiplicative spatial model class interface definition
+ * @author Juergen Knoedlseder
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "GCTAModelSpatialMultiplicative.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class GCTAModelSpatialMultiplicative
+ *
+ * @brief Multiplicative spatial model class
+ ***************************************************************************/
+class GCTAModelSpatialMultiplicative : public GCTAModelSpatial {
+
+public:
+    // Constructors and destructors
+    GCTAModelSpatialMultiplicative(void);
+    explicit GCTAModelSpatialMultiplicative(const GXmlElement& xml);
+    GCTAModelSpatialMultiplicative(const GCTAModelSpatialMultiplicative& model);
+    virtual ~GCTAModelSpatialMultiplicative(void);
+
+    // Operators
+    virtual GCTAModelSpatialMultiplicative& operator=(const GCTAModelSpatialMultiplicative& model);
+
+    // Pure virtual methods
+    virtual void                            clear(void);
+    virtual GCTAModelSpatialMultiplicative* clone(void) const;
+    virtual std::string                     classname(void) const;
+    virtual std::string                     type(void) const;
+    virtual double                          eval(const GCTAInstDir& dir,
+                                                 const GEnergy&     energy,
+                                                 const GTime&       time,
+                                                 const bool&        gradients = false) const;
+    virtual GCTAInstDir                     mc(const GEnergy& energy,
+                                               const GTime&   time,
+                                               GRan& ran) const;
+    virtual void                            read(const GXmlElement& xml);
+    virtual void                            write(GXmlElement& xml) const;
+
+    // Other methods
+    void                    append(const GCTAModelSpatial& spatial,
+                                   const std::string&      name="");
+    int                     components(void) const;
+    const GCTAModelSpatial* component(const int& index) const;
+    const GCTAModelSpatial* component(const std::string& name) const;
+};
+
+
+/***********************************************************************//**
+ * @brief GCTAModelSpatialMultiplicative class extension
+ ***************************************************************************/
+%extend GCTAModelSpatialMultiplicative {
+    GCTAModelSpatialMultiplicative copy() {
+        return (*self);
+    }
+%pythoncode {
+    def __getstate__(self):
+        xml = gammalib.GXmlElement()
+        self.write(xml)
+        state = (xml,)
+        return state
+    def __setstate__(self, state):
+        self.__init__(state[0])
+}
+};
@@ -307,6 +307,8 @@
 %include "GCTAModelIrfBackground.i"
 %include "GCTAModelAeffBackground.i"
 %include "GCTAModelSpatial.i"
+%include "GCTAModelSpatialGradient.i"
+%include "GCTAModelSpatialMultiplicative.i"
 %include "GCTAModelSpatialRegistry.i"
 %include "GCTAModelRadial.i"
 %include "GCTAModelRadialRegistry.i"
@@ -153,15 +153,15 @@ GCTACubeBackground::GCTACubeBackground(const GCTACubeBackground&amp; bgd)
 /***********************************************************************//**
  * @brief Background cube constructor
  *
- * @param[in] wcs      World Coordinate System.
- * @param[in] coords   Coordinate System (CEL or GAL).
- * @param[in] x        X coordinate of sky map centre (deg).
- * @param[in] y        Y coordinate of sky map centre (deg).
- * @param[in] dx       Pixel size in x direction at centre (deg/pixel).
- * @param[in] dy       Pixel size in y direction at centre (deg/pixel).
- * @param[in] nx       Number of pixels in x direction.
- * @param[in] ny       Number of pixels in y direction.
- * @param[in] energies Energies.
+ * @param[in] wcs     World Coordinate System.
+ * @param[in] coords  Coordinate System (CEL or GAL).
+ * @param[in] x       X coordinate of sky map centre (deg).
+ * @param[in] y       Y coordinate of sky map centre (deg).
+ * @param[in] dx      Pixel size in x direction at centre (deg/pixel).
+ * @param[in] dy      Pixel size in y direction at centre (deg/pixel).
+ * @param[in] nx      Number of pixels in x direction.
+ * @param[in] ny      Number of pixels in y direction.
+ * @param[in] ebounds Energy boundaries.
  *
  * Constructs a background cube by specifying the sky map grid and the
  * energies.
@@ -57,6 +57,7 @@
 //#define G_SMOOTH_EDISP_SAVE_TEST //!< Save test matrix
  
 /* __ Debug definitions __________________________________________________ */
+//#define G_DEBUG_HESS_RENORMALIZATION
  
 /* __ Constants __________________________________________________________ */
  
@@ -177,27 +178,28 @@ GCTAEdisp2D&amp; GCTAEdisp2D::operator=(const GCTAEdisp2D&amp; edisp)
  
  
 /***********************************************************************//**
- * @brief Return energy dispersion in units of \f$(\log_{10} MeV)^{-1}\f$
+ * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
  *
- * @param[in] logEobs Log10 of the measured energy (TeV).
- * @param[in] logEsrc Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad).
- * @param[in] phi Azimuth angle in camera system (rad). Not used in this method.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
+ * @param[in] logEobs Log10 of the observed energy (\f$\log_{10}\f$ TeV).
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians).
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Energy dispersion (MeV\f$^{-1}\f$)
  *
  * Returns the energy dispersion
  *
  * \f[
- *    E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}, \theta) =
- *    \frac{dP}{d\log_{10} E_{\rm reco}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}, \theta) =
+ *    \frac{E_{\rm disp}(E_{\rm true}, m, \theta)}{E_{\rm true}}
  * \f]
  *
- * in units of \f$(\log_{10} MeV)^{-1}\f$ where
- * \f$\log_{10} E_{\rm reco}\f$ is the logarithm of the reconstructed energy
- * in TeV,
- * \f$\log_{10} E_{\rm true}\f$ is the logarithm of the true energy in TeV and
- * \f$\theta\f$ is the offset angle in radians.
+ * in units of MeV\f$^{-1}\f$ where
+ * \f$\log_{10} E_{\rm reco}\f$ is the base 10 logarithm of the reconstructed
+ * energy in TeV,
+ * \f$\log_{10} E_{\rm true}\f$ is the base 10 logarithm of the true energy
+ * in TeV and \f$\theta\f$ is the offset angle in radians.
  ***************************************************************************/
 double GCTAEdisp2D::operator()(const double& logEobs, 
                                const double& logEsrc, 
@@ -236,6 +238,19 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
                 edisp = 0.0;
             }
  
+            // If energy dispersion is positive then compute true energy
+            // in MeV and divide the energy dispersion value by this
+            // quantity to get the energy dispersion per MeV
+            else {
+
+                // Compute true energy
+                double etrue = std::pow(10.0, logEsrc+6.0);
+
+                // Compute energy dispersion per MeV
+                edisp /= etrue;
+
+            }
+
         } // endif: migration was in valid range
  
     } // endif: true energy and offset angle was in valid range
@@ -324,7 +339,7 @@ void GCTAEdisp2D::table(const GCTAResponseTable&amp; table)
  *     MATRIX   - Migration matrix
  *
  * The data are stored in the m_edisp member. The energy axis will be set to
- * log10, the offset angle axis to radians.
+ * \f$log_{10} E_{\rm true}\f$, the offset angle axis to radians.
  *
  * The method assures that the energy dispersion is properly normalised.
  ***************************************************************************/
@@ -404,7 +419,7 @@ void GCTAEdisp2D::load(const GFilename&amp; filename)
  * @brief Save energy dispersion table into FITS file
  *
  * @param[in] filename FITS file name.
- * @param[in] clobber Overwrite existing file? (default: false)
+ * @param[in] clobber Overwrite existing file?
  *
  * Saves energy dispersion into a FITS file. If a file with the given
  * @p filename does not yet exist it will be created, otherwise the method
@@ -454,11 +469,11 @@ void GCTAEdisp2D::save(const GFilename&amp; filename, const bool&amp; clobber) const
  * @brief Simulate energy dispersion
  *
  * @param[in] ran Random number generator.
- * @param[in] logEsrc Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad).
- * @param[in] phi Azimuth angle in camera system (rad).
- * @param[in] zenith Zenith angle in Earth system (rad).
- * @param[in] azimuth Azimuth angle in Earth system (rad).
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians).
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
  * @return Observed energy.
  *
  * @exception GException::invalid_return_value
@@ -565,17 +580,16 @@ GEnergy GCTAEdisp2D::mc(GRan&amp;         ran,
 /***********************************************************************//**
  * @brief Return observed energy interval that contains the energy dispersion.
  *
- * @param[in] logEsrc Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad).
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians).
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
  * @return Observed energy boundaries.
  *
  * Returns the band of observed energies outside of which the energy
  * dispersion becomes negligible for a given true energy @p logEsrc and
- * offset angle @p theta. An energy is considered negligible if inferior
- * to 1e-6.
+ * offset angle @p theta.
  ***************************************************************************/
 GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc,
                                   const double& theta,
@@ -632,17 +646,16 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
 /***********************************************************************//**
  * @brief Return true energy interval that contains the energy dispersion.
  *
- * @param[in] logEobs Log10 of the observed event energy (TeV).
- * @param[in] theta Offset angle in camera system (rad).
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians).
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
  * @return True energy boundaries.
  *
  * Returns the band of true photon energies outside of which the energy
  * dispersion becomes negligible for a given observed energy @p logEobs and
- * offset angle @p theta. An energy in considered negligible if inferior to
- * 1e-6.
+ * offset angle @p theta.
  ***************************************************************************/
 GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs,
                                   const double& theta,
@@ -799,12 +812,20 @@ void GCTAEdisp2D::fetch(void) const
  * @param[in] ereco_min Minimum of reconstructed energy interval.
  * @param[in] ereco_max Maximum of reconstructed energy interval.
  * @param[in] etrue True energy.
- * @param[in] theta Offset angle (rad).
+ * @param[in] theta Offset angle (radians).
  * @return Integrated energy dispersion probability.
  *
- * Integrates the energy dispersion probability over an interval in
- * reconstructed energy, defined by @p ereco_min and @p ereco_max, for a
- * given true energy @p etrue and offset angle @p theta.
+ * Computes
+ *
+ * \f[
+ *    \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}, \theta) \, dE_{\rm reco}
+ * \f]
+ *
+ * where
+ * \f$E_{\rm reco}\f$ is the reconstructed energy,
+ * \f$E_{\rm true}\f$ is the true energy and
+ * \f$\theta\f$ is the offset angle.
  *
  * The method takes into account that the energy dispersion data are stored
  * in a matrix and uses the trapezoidal rule for integration of the tabulated
@@ -900,7 +921,7 @@ double GCTAEdisp2D::prob_erecobin(const GEnergy&amp; ereco_min,
             }
  
             // Compute integral
-            prob += 0.5 * (f1 + f2) * std::log10(x2/x1);
+            prob += 0.5 * (f1 + f2) * (x2 - x1);
  
             // Set second node as first node
             x1 = x2;
@@ -922,7 +943,7 @@ double GCTAEdisp2D::prob_erecobin(const GEnergy&amp; ereco_min,
                                 etrue_nodes.wgt_right(),
                                 theta_nodes.wgt_left(),
                                 theta_nodes.wgt_right(), x2);
-            prob += 0.5 * (f1 + f2) * std::log10(x2/x1);
+            prob += 0.5 * (f1 + f2) * (x2 - x1);
         }
  
     } // endif: logEsrc, theta and migration range were valid
@@ -1095,10 +1116,10 @@ void GCTAEdisp2D::free_members(void)
 /***********************************************************************//**
  * @brief Compute ebounds_obs vector
  *
- * @param[in] theta Offset angle (rad).
- * @param[in] phi Azimuth angle (rad).
- * @param[in] zenith Zenith angle (rad).
- * @param[in] azimuth Azimuth angle (rad).
+ * @param[in] theta Offset angle (radians).
+ * @param[in] phi Azimuth angle (radians).
+ * @param[in] zenith Zenith angle (radians).
+ * @param[in] azimuth Azimuth angle (radians).
  *
  * Computes for all true energies the energy boundaries of the observed
  * energies covered by valid migration matrix elements. Only matrix elements
@@ -1207,10 +1228,10 @@ void GCTAEdisp2D::compute_ebounds_obs(const double&amp; theta,
 /***********************************************************************//**
  * @brief Compute ebounds_src vector
  *
- * @param[in] theta Offset angle (rad).
- * @param[in] phi Azimuth angle (rad).
- * @param[in] zenith Zenith angle (rad).
- * @param[in] azimuth Azimuth angle (rad).
+ * @param[in] theta Offset angle (radians).
+ * @param[in] phi Azimuth angle (radians).
+ * @param[in] zenith Zenith angle (radians).
+ * @param[in] azimuth Azimuth angle (radians).
  *
  * Computes for all observed energies the energy boundaries of the true
  * energies covered by valid migration matrix elements. Only matrix elements
@@ -1328,8 +1349,18 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
  * @brief Set table
  *
  * After assigning or loading a response table, performs all necessary steps
- * to set the table. The method also smoothes the table get rid of numerical
- * fluctuations and normalises the table properly.
+ * to set the table. Specifically, the method normalizes the energy
+ * dispersion for each offset angle and true energy so that
+ *
+ * \f[
+ *    \int_{m_{\rm min}}^{m_{\rm max}}
+ *    E_{\rm disp}(E_{\rm true}, m, \theta) \, dm = 1
+ * \f]
+ *
+ * where \f$m=E_{\rm reco}/E_{\rm true}\f$.
+ *
+ * For CTA, the method applies a kludge that smoothes the energy dispersion
+ * table to get rid of numerical fluctuations.
  *
  * Sets the data members
  *      m_inx_etrue
@@ -1447,21 +1478,17 @@ void GCTAEdisp2D::set_max_edisp(void) const
  * Normalize the energy dispersion table using
  *
  * \f[
- *    \int_{\log_{10} E_{\rm reco}^{\rm min}}^{\log_{10} E_{\rm reco}^{\rm max}}
- *    E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}, \theta)
- *    \, d(\log_{10} E_{\rm reco}) = 1
+ *    \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}, \theta) \, dE_{\rm reco} = 1
  * \f]
  *
  * where
- * \f$E_{\rm reco}^{\rm min}\f$ is the minimum and
- * \f$E_{\rm reco}^{\rm max}\f$ the maximum reconstructed energy for a given
- * true energy, returned by the method ebounds_obs(),
- * \f$E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}, \theta)\f$
- * is the energy dispersion returned by the GCTAEdisp2D::operator(), given in
- * units of \f$(\log_{10} MeV)^{-1}\f$,
- * \f$\log_{10} E_{\rm reco}\f$ is the logarithm of the reconstructued energy
- * in TeV, and
- * \f$\log_{10} E_{\rm true}\f$ is the logarithm of the true energy in TeV.
+ * \f$E_{\rm reco}_{\rm min}\f$ and \f$E_{\rm reco}_{\rm max}\f$ is the
+ * minimum and maximum migration value for a given true energy as returned
+ * by the ebounds_obs() method, and
+ * \f$E_{\rm disp}(E_{\rm true}, E_{\rm reco}, \theta)\f$ is the energy
+ * dispersion returned by the GCTAEdisp2D::operator(), given in units of
+ * \f$MeV^{-1}\f$.
  *
  * The normalisation is performed for each \f$E_{\rm true}\f$ and
  * \f$\theta\f$ bin using a Romberg integration method. Since the
@@ -1498,7 +1525,7 @@ void GCTAEdisp2D::normalize_table(void)
             // Loop over true photon energy
             for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue) {
  
-                // Get energy
+                // Get the base 10 logarithm of the true energy
                 double emin    = std::log10(m_edisp.axis_lo(m_inx_etrue,i_etrue));
                 double emax    = std::log10(m_edisp.axis_hi(m_inx_etrue,i_etrue));
                 double logEsrc = 0.5*(emin+emax);
@@ -1512,9 +1539,9 @@ void GCTAEdisp2D::normalize_table(void)
                 // Loop over all energy intervals
                 for (int i = 0; i < ebounds.size(); ++i) {
  
-                    // Get energy boundaries
-                    emin = ebounds.emin(i).log10TeV();
-                    emax = ebounds.emax(i).log10TeV();
+                    // Get energy boundaries in log10 of energy in MeV
+                    emin = ebounds.emin(i).log10MeV();
+                    emax = ebounds.emax(i).log10MeV();
  
                     // Setup integration function
                     edisp_kern integrand(this, logEsrc, theta);
@@ -1707,6 +1734,9 @@ double GCTAEdisp2D::table_value(const int&amp;    base_ll,
  
 /***********************************************************************//**
  * @brief Smoothed energy dispersion table
+ *
+ * This method implements the kludge for CTA that reduces the statistical
+ * noise in the energy dispersion matrix.
  ***************************************************************************/
 void GCTAEdisp2D::smooth_table(void)
 {
@@ -2074,23 +2104,47 @@ GNdarray GCTAEdisp2D::gaussian_array(const double&amp; mean,
  
  
 /***********************************************************************//**
- * @brief Integration kernel for edisp_kern() class
+ * @brief Integration kernel for GCTAEdisp2D::edisp_kern class
+ *
+ * @param[in] logEobs Log10 of observation energy (\f$\log_{10}\f$ MeV).
+ * @return Energy dispersion (\f$(\log_{10}\f$ MeV\f$)^{-1}\f$).
+ *
+ * This method implements the function
+ *
+ * \f[
+ *    E_{\rm disp}(E_{\rm true}, \log_{10} E_{\rm reco}, \theta) =
+ *    E_{\rm disp}(E_{\rm true}, m, \theta) \times
+ *    \frac{\log 10 \times E_{\rm reco}}{E_{\rm true}}
+ * \f]
  *
- * @param[in] x Function value.
+ * which is the integration kernel needed for the GCTAEdisp2D::edisp_kern
+ * class. The class is used by GCTAEdisp2D::normalize_table to integrate the
+ * energy dispersion information using
  *
- * This method implements the integration kernel needed for the edisp_kern()
- * class.
+ * \f[
+ *    \int_{\log_{10} E_{\rm reco}^{\rm min}}^{\log_{10} E_{\rm reco}^{\rm max}}
+ *    E_{\rm disp}(E_{\rm true}, \log_{10} E_{\rm reco}, \theta) \,
+ *    d(\log_{10} E_{\rm reco})
+ * \f]
  ***************************************************************************/
-double GCTAEdisp2D::edisp_kern::eval(const double& x)
+double GCTAEdisp2D::edisp_kern::eval(const double& logEobs)
 {
-    // Get function value
-    double value = m_parent->operator()(x, m_logEsrc, m_theta);
+    // Get function value. Note that logEobs is in log10 MeV, while we need
+    // the value in log10 TeV.
+    double value = m_parent->operator()(logEobs-6.0, m_logEsrc, m_theta);
+
+    // Compute reconstructed energy in MeV. Note that logEobs is in
+    // log10 MeV!
+    double ereco = std::pow(10.0, logEobs);
+
+    // Normalize energy dispersion so that the units are per log10 MeV
+    value *= gammalib::ln10 * ereco;
  
     // Compile option: Check for NaN
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
         std::cout << "*** ERROR: GCTAEdisp2D::edisp_kern::eval";
-        std::cout << "(x=" << x << "): ";
+        std::cout << "(logEobs=" << logEobs << "): ";
         std::cout << " NaN/Inf encountered";
         std::cout << " (value=" << value;
         std::cout << ")" << std::endl;
@@ -158,34 +158,31 @@ GCTAEdispPerfTable&amp; GCTAEdispPerfTable::operator=(const GCTAEdispPerfTable&amp; edis
  
  
 /***********************************************************************//**
- * @brief Return energy dispersion in units of \f$(\log_{10} MeV)^{-1}\f$
+ * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
  *
- * @param[in] logEobs log10 of the observed photon energy (TeV).
- * @param[in] logEsrc log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEobs Log10 of the observed photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Energy dispersion (MeV\f$^{-1}\f$)
  *
  * Returns the energy dispersion
  *
  * \f[
- *    E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}) =
- *    \frac{dP}{d\log_{10} E_{\rm reco}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}) =
+ *    \frac{1}{\sqrt{2\pi}\sigma(E_{\rm true})}
+ *    \exp \left(\frac{-(\log_{10} E_{\rm reco} - \log_{10} E_{\rm true})^2}
+ *                    {2 \sigma(E_{\rm true})^2} \right) \times
+ *    \frac{1}{\log_{10} E_{\rm reco}}
  * \f]
  *
- * in units of \f$(\log_{10} MeV)^{-1}\f$ where
- * \f$\log_{10} E_{\rm reco}\f$ is the logarithm of the reconstructed energy
- * in TeV, and
- * \f$\log_{10} E_{\rm true}\f$ is the logarithm of the true energy in TeV.
- *
- * The operator evaluates
- *
- * \f[
- *    E_{\rm disp} = \frac{1}{\sqrt{2\pi}m\_sigma}
- *                   \exp(\frac{-(\log_{10} E_{\rm reco} -
- *                                \log_{10} E_{\rm true})^2}{2 m\_sigma^2})
- * \f]
+ * in units of MeV\f$^{-1}\f$ where
+ * \f$E_{\rm reco}\f$ is the reconstructed energy in units of MeV,
+ * \f$E_{\rm true}\f$ is the true energy in units of MeV, and
+ * \f$\sigma(E_{\rm true})\f$ is the standard deviation of the energy
+ * dispersion that depends on the true photon energy.
  ***************************************************************************/
 double GCTAEdispPerfTable::operator()(const double& logEobs,
                                       const double& logEsrc,
@@ -201,6 +198,12 @@ double GCTAEdispPerfTable::operator()(const double&amp; logEobs,
     double delta = logEobs - logEsrc;
     double edisp = m_par_scale * std::exp(m_par_width * delta * delta);
  
+    // Compute reconstructed energy in MeV
+    double ereco = std::pow(10.0, logEobs+6.0);
+
+    // Compute energy dispersion per MeV
+    edisp /= (gammalib::ln10 * ereco);
+
     // Return energy dispersion
     return edisp;
 }
@@ -318,11 +321,12 @@ void GCTAEdispPerfTable::load(const GFilename&amp; filename)
  * @brief Simulate energy dispersion
  *
  * @param[in] ran Random number generator.
- * @param[in] logEsrc Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Energy.
  *
  * Draws observed energy value from a normal distribution of width
  * m_par_sigma around @p logE.
@@ -352,11 +356,12 @@ GEnergy GCTAEdispPerfTable::mc(GRan&amp;         ran,
 /***********************************************************************//**
  * @brief Return observed energy interval that contains the energy dispersion.
  *
- * @param[in] logEsrc Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Boundaries of reconstruced energies.
  *
  * Returns the band of observed energies outside of which the energy
  * dispersion becomes negligible for a given true energy @p logEsrc. This
@@ -389,15 +394,16 @@ GEbounds GCTAEdispPerfTable::ebounds_obs(const double&amp; logEsrc,
 /***********************************************************************//**
  * @brief Return true energy interval that contains the energy dispersion.
  *
- * @param[in] logEobs Log10 of the observed event energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Boundaries of true energies.
  *
  * Returns the band of true photon energies outside of which the energy
- * dispersion becomes negligible for a given observed energy @p logEobs. This
- * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
+ * dispersion becomes negligible for a given observed energy @p logEobs.
+ * This band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
  * Gaussian width of the energy dispersion.
  ***************************************************************************/
 GEbounds GCTAEdispPerfTable::ebounds_src(const double& logEobs,
@@ -430,8 +436,19 @@ GEbounds GCTAEdispPerfTable::ebounds_src(const double&amp; logEobs,
  * @param[in] ereco_min Minimum of reconstructed energy interval.
  * @param[in] ereco_max Maximum of reconstructed energy interval.
  * @param[in] etrue True energy.
- * @param[in] theta Offset angle (not used).
+ * @param[in] theta Offset angle (radians). Not used.
  * @return Integrated energy dispersion probability.
+ *
+ * Computes
+ *
+ * \f[
+ *    \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}) \, dE_{\rm reco}
+ * \f]
+ *
+ * where
+ * \f$E_{\rm reco}\f$ is the reconstructed energy and
+ * \f$E_{\rm true}\f$ is the true energy.
  ***************************************************************************/
 double GCTAEdispPerfTable::prob_erecobin(const GEnergy& ereco_min,
                                          const GEnergy& ereco_max,
@@ -556,11 +573,12 @@ void GCTAEdispPerfTable::free_members(void)
 /***********************************************************************//**
  * @brief Update energy dispersion parameter cache
  *
- * @param[in] logE Log10 of the true photon energy (TeV).
+ * @param[in] logE Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
  *
- * This method updates the energy dispersion parameter cache.
- * As the performance table energy dispersion only depends on energy,
- * the only parameter on which the cache values depend is the energy.
+ * This method updates the energy dispersion parameter cache. As the
+ * performance table energy dispersion only depends on true photon energy,
+ * the only parameter on which the cache values depend is the true photon
+ * energy.
  ***************************************************************************/
 void GCTAEdispPerfTable::update(const double& logE) const
 {
@@ -160,27 +160,27 @@ GCTAEdispRmf&amp; GCTAEdispRmf::operator=(const GCTAEdispRmf&amp; edisp)
  
  
 /***********************************************************************//**
- * @brief Return energy dispersion in units of \f$(\log_{10} MeV)^{-1}\f$
+ * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
  *
- * @param[in] logEobs log10 of the observed photon energy (TeV).
- * @param[in] logEsrc log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
- * @return Energy dispersion
+ * @param[in] logEobs Log10 of the observed photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Energy dispersion (MeV\f$^{-1}\f$)
  *
  * Returns the energy dispersion
  *
  * \f[
- *    E_{\rm disp}(\log_{10} E_{\rm reco} | \log_{10} E_{\rm true}) =
- *    \frac{dP}{d\log_{10} E_{\rm reco}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}) =
+ *    \frac{E_{\rm disp}(E_{\rm true}, \log_{10} E_{\rm reco}, \theta)}
+ *         {\log_{10} E_{\rm reco}}
  * \f]
  *
- * in units of \f$(\log_{10} MeV)^{-1}\f$ where
- * \f$\log_{10} E_{\rm reco}\f$ is the logarithm of the reconstructed energy
- * in TeV, and
- * \f$\log_{10} E_{\rm true}\f$ is the logarithm of the true energy in TeV.
+ * in units of MeV\f$^{-1}\f$ where
+ * \f$E_{\rm reco}\f$ is the reconstructed energy in units of MeV, and
+ * \f$E_{\rm true}\f$ is the true energy in units of MeV.
  ***************************************************************************/
 double GCTAEdispRmf::operator()(const double& logEobs,
                                 const double& logEsrc,
@@ -203,6 +203,11 @@ double GCTAEdispRmf::operator()(const double&amp; logEobs,
         edisp = 0.0;
     }
  
+    // If the energy dispersion is positive, return it per MeV
+    else {
+        edisp /= (gammalib::ln10 * std::pow(10.0, logEobs+6.0));
+    }
+
     // Return energy dispersion
     return edisp;
 }
@@ -278,11 +283,12 @@ void GCTAEdispRmf::load(const GFilename&amp; filename)
  * @brief Simulate energy dispersion
  *
  * @param[in] ran Random number generator.
- * @param[in] logEsrc Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Energy.
  *
  * Draws observed energy value from RMF matrix.
  ***************************************************************************/
@@ -323,11 +329,12 @@ GEnergy GCTAEdispRmf::mc(GRan&amp;         ran,
 /***********************************************************************//**
  * @brief Return observed energy interval that contains the energy dispersion.
  *
- * @param[in] logEsrc Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return Reconstructed energy boundaries.
  *
  * Returns the band of observed energies outside of which the energy
  * dispersion becomes negligible for a given true energy @p logEsrc.
@@ -383,11 +390,12 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double&amp; logEsrc,
 /***********************************************************************//**
  * @brief Return true energy interval that contains the energy dispersion.
  *
- * @param[in] logEobs Log10 of the observed event energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
+ * @return True energy boundaries.
  *
  * Returns the band of true photon energies outside of which the energy
  * dispersion becomes negligible for a given observed energy @p logEobs.
@@ -447,8 +455,19 @@ GEbounds GCTAEdispRmf::ebounds_src(const double&amp; logEobs,
  * @param[in] ereco_min Minimum of reconstructed energy interval.
  * @param[in] ereco_max Maximum of reconstructed energy interval.
  * @param[in] etrue True energy.
- * @param[in] theta Offset angle (not used).
+ * @param[in] theta Offset angle. Not used.
  * @return Integrated energy dispersion probability.
+ *
+ * Computes
+ *
+ * \f[
+ *    \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
+ *    E_{\rm disp}(E_{\rm true}, E_{\rm reco}) \, dE_{\rm reco}
+ * \f]
+ *
+ * where
+ * \f$E_{\rm reco}\f$ is the reconstructed energy and
+ * \f$E_{\rm true}\f$ is the true energy.
  ***************************************************************************/
 double GCTAEdispRmf::prob_erecobin(const GEnergy& ereco_min,
                                    const GEnergy& ereco_max,
@@ -783,25 +802,25 @@ void GCTAEdispRmf::set_max_edisp(void) const
 /***********************************************************************//**
  * @brief Update cache
  *
- * @param[in] etrue Log10 of true energy in TeV.
- * @param[in] emeasured Log10 of measured energy in TeV.
+ * @param[in] logEsrc Log10 of true energy (\f$\log_{10}\f$ TeV).
+ * @param[in] logEobs Log10 of reconstructed energy (\f$\log_{10}\f$ TeV).
  *
  * Updates the interpolation cache. The interpolation cache is composed
  * of four indices and weights that define 4 data values of the RMF matrix
  * that are used for bilinear interpolation.
  ***************************************************************************/
-void GCTAEdispRmf::update(const double& etrue, const double& emeasured) const
+void GCTAEdispRmf::update(const double& logEsrc, const double& logEobs) const
 {
     // Update cache only of arguments have changed
-    if (etrue != m_last_etrue || emeasured != m_last_emeasured) {
+    if (logEsrc != m_last_etrue || logEobs != m_last_emeasured) {
  
         // Store actual values
-        m_last_etrue     = etrue;
-        m_last_emeasured = emeasured;
+        m_last_etrue     = logEsrc;
+        m_last_emeasured = logEobs;
  
         // Set values for node arrays
-        m_etrue.set_value(etrue);
-        m_emeasured.set_value(emeasured);
+        m_etrue.set_value(logEsrc);
+        m_emeasured.set_value(logEobs);
  
         // Set indices for bi-linear interpolation
         m_itrue1 = m_etrue.inx_left();
@@ -825,10 +844,10 @@ void GCTAEdispRmf::update(const double&amp; etrue, const double&amp; emeasured) const
 /***********************************************************************//**
  * @brief Compute ebounds_obs vector
  *
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
  ***************************************************************************/
 void GCTAEdispRmf::compute_ebounds_obs(const double& theta,
                                        const double& phi,
@@ -861,10 +880,10 @@ void GCTAEdispRmf::compute_ebounds_obs(const double&amp; theta,
 /***********************************************************************//**
  * @brief Compute ebounds_src vector
  *
- * @param[in] theta Offset angle in camera system (rad). Not used.
- * @param[in] phi Azimuth angle in camera system (rad). Not used.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
+ * @param[in] theta Offset angle in camera system (radians). Not used.
+ * @param[in] phi Azimuth angle in camera system (radians). Not used.
+ * @param[in] zenith Zenith angle in Earth system (radians). Not used.
+ * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
  ***************************************************************************/
 void GCTAEdispRmf::compute_ebounds_src(const double& theta,
                                        const double& phi,
@@ -75,6 +75,28 @@ GCTAModelSpatialGradient::GCTAModelSpatialGradient(void) : GCTAModelSpatial()
  
  
 /***********************************************************************//**
+ * @brief Gradient constructor
+ *
+ * @param[in] detx_gradient DETX gradient (degrees\f$^{-1}\f$).
+ * @param[in] dety_gradient DETY gradient (degrees\f$^{-1}\f$).
+ ***************************************************************************/
+GCTAModelSpatialGradient::GCTAModelSpatialGradient(const double& detx_gradient,
+                                                   const double& dety_gradient) :
+                          GCTAModelSpatial()
+{
+    // Initialise members
+    init_members();
+
+    // Assign gradients
+    this->detx_gradient(detx_gradient);
+    this->dety_gradient(dety_gradient);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief XML constructor
  *
  * @param[in] xml XML element.
@@ -318,7 +340,7 @@ void GCTAModelSpatialGradient::read(const GXmlElement&amp; xml)
  *
  * @param[in] xml XML element.
  *
- * @exceptino GException::invalid_value
+ * @exception GException::invalid_value
  *            Spatial model is not of valid type.
  *
  * Write the gradient spatial model information into an XML element.
@@ -115,6 +115,7 @@
 //#define G_DEBUG_IRF_RADIAL                     //!< Debug irf_radial method
 //#define G_DEBUG_IRF_DIFFUSE                   //!< Debug irf_diffuse method
 //#define G_DEBUG_IRF_ELLIPTICAL             //!< Debug irf_elliptical method
+//#define G_DEBUG_NROI_EDISP    //!< Debug nroi energy dispersion computation
 //#define G_DEBUG_NROI_RADIAL                  //!< Debug npred_radial method
 //#define G_DEBUG_NROI_ELLIPTICAL          //!< Debug npred_elliptical method
 //#define G_DEBUG_NROI_DIFFUSE                //!< Debug npred_diffuse method
@@ -455,12 +456,12 @@ double GCTAResponseIrf::nroi(const GModelSky&amp;    model,
         // Loop over all boundaries
         for (int i = 0; i < ebounds.size(); ++i) {
  
-            // Get boundaries in MeV
-            double emin = ebounds.emin(i).MeV();
-            double emax = ebounds.emax(i).MeV();
+            // Get true energy boundaries in log10 MeV
+            double log10_etrue_min = ebounds.emin(i).log10MeV();
+            double log10_etrue_max = ebounds.emax(i).log10MeV();
  
             // Continue only if valid
-            if (emax > emin) {
+            if (log10_etrue_max > log10_etrue_min) {
  
                 // Setup integration function
                 cta_nroi_kern integrand(this, &obs, &model, srcTime, obsEng, obsTime);
@@ -470,14 +471,32 @@ double GCTAResponseIrf::nroi(const GModelSky&amp;    model,
                 integral.fixed_iter(iter);
  
                 // Do Romberg integration
-                emin  = std::log(emin);
-                emax  = std::log(emax);
-                nroi += integral.romberg(emin, emax);
+                nroi += integral.romberg(log10_etrue_min, log10_etrue_max);
  
             } // endif: interval was valid
  
         } // endfor: looped over energy intervals
  
+        // Debug energy dispersion computation
+        #if defined(G_DEBUG_NROI_EDISP)
+        apply_edisp(false);
+        const GEnergy& srcEng = obsEng;
+        double nroi_spatial   = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
+        double nroi_spectral  = model.spectral()->eval(srcEng, srcTime);
+        double nroi_temporal  = model.temporal()->eval(srcTime);
+        double nroi_noedisp   = nroi_spatial * nroi_spectral * nroi_temporal;
+        apply_edisp(true);
+        std::cout << "GCTAResponseIrf::nroi:";
+        std::cout << " obs_id=" << obs.id();
+        std::cout << " ereco=" << obsEng;
+        std::cout << " nroi(noedisp)=" << nroi_noedisp;
+        std::cout << " nroi(edisp)=" << nroi;
+        if (nroi_noedisp > 0.0) {
+            std::cout << " ratio(edisp/noedisp)=" << nroi/nroi_noedisp;
+        }
+        std::cout << std::endl;
+        #endif
+
     } // endif: energy dispersion requested
  
     // ... otherwise evaluate
@@ -487,12 +506,12 @@ double GCTAResponseIrf::nroi(const GModelSky&amp;    model,
         const GEnergy& srcEng = obsEng;
  
         // Compute response components
-        double npred_spatial  = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
-        double npred_spectral = model.spectral()->eval(srcEng, srcTime);
-        double npred_temporal = model.temporal()->eval(srcTime);
+        double nroi_spatial  = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
+        double nroi_spectral = model.spectral()->eval(srcEng, srcTime);
+        double nroi_temporal = model.temporal()->eval(srcTime);
  
         // Compute response
-        nroi = npred_spatial * npred_spectral * npred_temporal;
+        nroi = nroi_spatial * nroi_spectral * nroi_temporal;
  
     } // endelse: no energy dispersion requested
  
@@ -1783,7 +1802,7 @@ double GCTAResponseIrf::psf_delta_max(const double&amp; theta,
  
  
 /***********************************************************************//**
- * @brief Return energy dispersion (in units of MeV^-1)
+ * @brief Return energy dispersion (in units of MeV\f$^{-1}\f$)
  *
  * @param[in] obsEng Measured event energy.
  * @param[in] theta Radial offset angle in camera (radians).
@@ -1791,6 +1810,7 @@ double GCTAResponseIrf::psf_delta_max(const double&amp; theta,
  * @param[in] zenith Zenith angle of telescope pointing (radians).
  * @param[in] azimuth Azimuth angle of telescope pointing (radians).
  * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
+ * @return Energy dispersion (MeV\f$^{-1}\f$).
  ***************************************************************************/
 double GCTAResponseIrf::edisp(const GEnergy& obsEng,
                               const double&  theta,
@@ -1799,13 +1819,11 @@ double GCTAResponseIrf::edisp(const GEnergy&amp; obsEng,
                               const double&  azimuth,
                               const double&  srcLogEng) const
 {
-    // Compute log10 energy in TeV and linear energy in MeV
+    // Compute log10 energy in TeV
     double obsLogEng = obsEng.log10TeV();
-    double energy    = obsEng.MeV();
  
     // Compute energy dispersion
-    double edisp = (*m_edisp)(obsLogEng, srcLogEng, theta, phi, zenith, azimuth) /
-                   (gammalib::ln10 * energy);
+    double edisp = (*m_edisp)(obsLogEng, srcLogEng, theta, phi, zenith, azimuth);
  
     // Return energy dispersion
     return edisp;
@@ -2707,11 +2725,12 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp;       event,
     // Retrieve CTA pointing
     const GCTAPointing& pnt = gammalib::cta_pnt(G_IRF_DIFFUSE, obs);
  
-    // Try getting the IRF value from cache
+    // Try getting the IRF value from cache. Do not do this for energy
+    // dispersion since the the spectral model will change the IRF values
     #if defined(G_USE_IRF_CACHE)
     const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(obs.events());
     const GCTAEventAtom* atom = dynamic_cast<const GCTAEventAtom*>(&event);
-    if (list != NULL && atom != NULL) {
+    if (list != NULL && atom != NULL && !use_edisp()) {
         irf = list->irf_cache(source.name(), atom->index());
         if (irf >= 0.0) {
             has_irf = true;
@@ -2826,9 +2845,11 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp;       event,
         // Apply deadtime correction
         irf *= obs.deadc(srcTime);
  
-        // Put IRF value in cache
+        // Put IRF value in cache. Only do this in no energy dispersion is
+        // requested, since for energy dispersion the IRF values will
+        // depend on the spectrum
         #if defined(G_USE_IRF_CACHE)
-        if (list != NULL && atom != NULL) {
+        if (list != NULL && atom != NULL && !use_edisp()) {
             list->irf_cache(source.name(), atom->index(), irf);
         }
         #endif
@@ -3315,9 +3336,11 @@ double GCTAResponseIrf::nroi_diffuse(const GModelSky&amp;    model,
     // Build unique identifier
     std::string id = model.name() + "::" + obs.id();
  
-    // Check if Nroi value is already in cache
+    // Check if Nroi value is already in cache. Do not use the Npred cache
+    // for energy dispersion since the actual value will depend on the
+    // spectral model
     #if defined(G_USE_NPRED_CACHE)
-    if (!m_npred_names.empty()) {
+    if (!m_npred_names.empty() && !use_edisp()) {
  
          // Search for unique identifier, and if found, recover Npred value
          // and break
@@ -3407,18 +3430,23 @@ double GCTAResponseIrf::nroi_diffuse(const GModelSky&amp;    model,
             #if defined(G_DEBUG_NROI_DIFFUSE)
             std::cout << "GCTAResponseIrf::nroi_diffuse:";
             std::cout << " roi_psf_radius=" << roi_psf_radius;
-            std::cout << " nroi=" << npred;
+            std::cout << " nroi=" << nroi;
+            std::cout << " Etrue=" << srcEng;
+            std::cout << " Ereco=" << obsEng;
             std::cout << " id=" << id << std::endl;
             #endif
  
         } // endif: offset angle range was valid
  
-        // Store result in Npred cache
+        // Store result in Npred cache in case that no energy dispersion
+        // is requested
         #if defined(G_USE_NPRED_CACHE)
-        m_npred_names.push_back(id);
-        m_npred_energies.push_back(srcEng);
-        m_npred_times.push_back(srcTime);
-        m_npred_values.push_back(nroi);
+        if (!use_edisp()) {
+            m_npred_names.push_back(id);
+            m_npred_energies.push_back(srcEng);
+            m_npred_times.push_back(srcTime);
+            m_npred_values.push_back(nroi);
+        }
         #endif
  
         // Debug: Check for NaN
@@ -133,7 +133,7 @@ cta_omega_intervals gammalib::limit_omega(const double&amp; min,
             }
             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
@@ -339,7 +339,7 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
 {
     // Initialise PSF value
     double value = 0.0;
-    
+
     // Get arclength for given radius in radians
     double phi = gammalib::cta_roi_arclength(delta,
                                              m_psf,
@@ -350,7 +350,7 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
  
     // If arclength is positive then compute the PSF value
     if (phi > 0) {
-    
+
         // Compute PSF value
         value = m_rsp->psf(delta, m_theta, m_phi, m_zenith, m_azimuth, m_logE) * 
                            phi * std::sin(delta);
@@ -366,7 +366,7 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
             std::cout << std::endl;
         }
         #endif
-        
+
     } // endif: arclength was positive
  
     // Return
@@ -375,25 +375,25 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
  
  
 /***********************************************************************//**
- * @brief Integration kernel for nedisp() method
+ * @brief Integration kernel for GCTAResponseIrf::nroi method
  *
- * @param[in] logEsrc Logarithm of true photon energy.
+ * @param[in] logEsrc Base 10 logarithm of true photon energy in MeV.
  * @return Nroi.
  ***************************************************************************/
 double cta_nroi_kern::eval(const double& logEsrc)
 {
     // Set true energy
     GEnergy srcEng;
-    double expx = std::exp(logEsrc);
-    srcEng.MeV(expx);
+    srcEng.log10MeV(logEsrc);
  
     // Compute response components
     double nroi_spatial  = m_rsp->nroi(*m_model, srcEng, m_srcTime, m_obsEng, m_obsTime, *m_obs);
     double nroi_spectral = m_model->spectral()->eval(srcEng, m_srcTime);
     double nroi_temporal = m_model->temporal()->eval(m_srcTime);
-    
+
     // Compute response
-    double nroi = nroi_spatial * nroi_spectral * nroi_temporal * expx;
+    double nroi = nroi_spatial * nroi_spectral * nroi_temporal *
+                  gammalib::ln10 * m_obsEng.MeV();
  
     // Return response
     return nroi;
@@ -567,10 +567,10 @@ double cta_irf_radial_kern_omega::eval(const double&amp; omega)
 {
     // Compute PSF offset angle [radians]
     double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
-    
+
     // Compute true photon offset angle in camera system [radians]
     double offset = std::acos(m_cos_ph + m_sin_ph * std::cos(m_omega0 - omega));
-    
+
     //TODO: Compute true photon azimuth angle in camera system [radians]
     double azimuth = 0.0;
  
@@ -582,7 +582,7 @@ double cta_irf_radial_kern_omega::eval(const double&amp; omega)
     if (m_rsp->use_edisp() && irf > 0.0) {
         irf *= m_rsp->edisp(m_obsEng, offset, azimuth, m_zenith, m_azimuth, m_srcLogEng);
     }
-    
+
     // Compile option: Check for NaN/Inf
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
@@ -944,7 +944,7 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
             #endif
  
         } // endif: arc length was positive
-    
+
     } // endif: rho was positive
  
     // Return result
@@ -1035,7 +1035,7 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
  
         // 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_omega_pnt - omega));
@@ -1242,7 +1242,7 @@ double cta_nroi_elliptical_kern_rho::eval(const double&amp; rho)
             #endif
  
         } // endif: arc length was positive
-    
+
     } // endif: rho was positive
  
     // Return Nroi
@@ -1302,10 +1302,10 @@ double cta_nroi_elliptical_kern_omega::eval(const double&amp; omega)
         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);
@@ -1339,7 +1339,7 @@ double cta_nroi_elliptical_kern_omega::eval(const double&amp; omega)
         #endif
  
     } // endif: sky intensity was positive
-    
+
     // Return Nroi
     return nroi;
 }
@@ -2210,7 +2210,7 @@ double cta_psf_elliptical_kern_rho::eval(const double&amp; rho)
             #endif
  
         } // endif: arc length was positive
-    
+
     } // endif: rho was positive
  
     // Return result
@@ -901,7 +901,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
     test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
  
     // Test that diagonal element is non-zero
-    test_value(edisp2(std::log10(30.0),std::log10(30.0)), 14.7025305813298);
+    test_value(edisp2(std::log10(30.0),std::log10(30.0)), 2.12840930049511e-07);
  
     // Test if non-diagonal element (above diagonal) is zero
     test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
@@ -911,7 +911,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
     GEnergy etrue(1.0, "TeV");
     GEnergy emin(0.1, "TeV");
     GEnergy emax(10.0, "TeV");
-    test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.0);
+    test_value(edisp2.prob_erecobin(emin, emax, etrue,  0.0), 1.0);
     test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.5);
     test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.5);
  
@@ -958,7 +958,7 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
     test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
  
     // Test that diagonal element is non-zero
-    test_value(edisp2(std::log10(30.0),std::log10(30.0)), 7.4480009344395);
+    test_value(edisp2(std::log10(30.0),std::log10(30.0)), 1.0785985662216e-07);
  
     // Test if non-diagonal element (above diagonal) is zero
     test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
@@ -967,13 +967,13 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
     GEnergy etrue(1.0, "TeV");
     GEnergy emin(0.1, "TeV");
     GEnergy emax(10.0, "TeV");
-    test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0),  0.999664174482241);
-    test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.494679101633454);
-    test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.504985072848788);
+    test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0),  1.00002576227557);
+    test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.42065218271377);
+    test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.43096191010834);
  
     // Test mc method
     GRan ran;
-    test_value(edisp2.mc(ran, 0.0).TeV(), 0.942937281664539);
+    test_value(edisp2.mc(ran, 0.0).TeV(), 0.930126476703161);
  
     // Test normalisation
     test_edisp_integration(edisp2, 0.5, 10.0);
@@ -1006,14 +1006,14 @@ void TestGCTAResponse::test_response_edisp_2D(void)
     test_value(edisp2.ebounds_obs(1.0).emax().TeV(), 13.2011115155942);
  
     // Test ebounds_src method
-    test_value(edisp2.ebounds_src(1.0).emin().TeV(),  7.72198100249182);
-    test_value(edisp2.ebounds_src(1.0).emax().TeV(), 19.8348918834691);
+    test_value(edisp2.ebounds_src(1.0).emin().TeV(),  7.7813808018878);
+    test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8059295273666);
  
     // Test if non-diagonal element (below diagonal) is zero
     test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
  
     // Test that diagonal element is non-zero
-    test_value(edisp2(std::log10(30.0),std::log10(30.0)), 11.7884028286084);
+    test_value(edisp2(std::log10(30.0),std::log10(30.0)), 1.74200329848547e-07);
  
     // Test if non-diagonal element (above diagonal) is zero
     test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
@@ -1022,13 +1022,13 @@ void TestGCTAResponse::test_response_edisp_2D(void)
     GEnergy etrue(1.0, "TeV");
     GEnergy emin(0.1, "TeV");
     GEnergy emax(10.0, "TeV");
-    test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0),  1.00008878464952);
-    test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.448839360815326);
-    test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.551249618415791);
+    test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0),  1.00002101900895);
+    test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.480451962185927);
+    test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.519569056823021);
  
     // Test mc method
     GRan ran;
-    test_value(edisp2.mc(ran, 0.0).TeV(), 0.97532927330811);
+    test_value(edisp2.mc(ran, 0.0).TeV(), 0.872387491302767);
  
     // Test normalisation
     test_edisp_integration(edisp2);
@@ -1411,9 +1411,9 @@ void TestGCTAResponse::test_response_edispcube(void)
 /***********************************************************************//**
  * @brief Utility function for energy dispersion tests
  *
- * @param[in] rsp Response
- * @param[in] e_src_min Minimum true energy (TeV) (default: 0.1)
- * @param[in] e_src_max Maximum true energy (TeV) (default: 10.0)
+ * @param[in] rsp Response.
+ * @param[in] e_src_min Minimum true energy (TeV).
+ * @param[in] e_src_max Maximum true energy (TeV).
  ***************************************************************************/
 void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf& rsp,
                                                        const double&          e_src_min,
@@ -1429,40 +1429,29 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf&amp; rs
 	        double log10_e_src = std::log10(e_src);
  
 	        // Retrieve boundaries in observed energy
-	        GEbounds ebounds  = rsp.edisp()->ebounds_obs(log10_e_src);
-	        GEnergy  emin     = ebounds.emin();
-	        GEnergy  emax     = ebounds.emax();
-	        double   logE_min = std::log10(emin.TeV());
-	        double   logE_max = std::log10(emax.TeV());
-
-	        // Compute step size for numerical integration
-	        const int steps = 1000;
-	        double    dlogE = (logE_max-logE_min)/steps;
+	        GEbounds ebounds = rsp.edisp()->ebounds_obs(log10_e_src);
+	        double   emin    = ebounds.emin().MeV();
+	        double   emax    = ebounds.emax().MeV();
  
 	        // Perform numerical integration by summing
-	        double sum      = 0.0;
-	        double logE_obs = logE_min;
-	        for (int i = 0; i < steps; ++i) {
-	            double dp_dlogE = (*rsp.edisp())(logE_obs, log10_e_src);
-	            sum            += dp_dlogE * dlogE;
-	            logE_obs       += dlogE;
-	        }
-	        GEnergy eng(e_src, "TeV");
-	        test_value(sum, 1.0, 0.001, "Energy Dispersion integration for "+eng.print());
-
-	        // And now in linear energies, test GCTAResponse::edisp() method
 	        const int nE = 10000;
-	        double E_min = emin.MeV();
-	        double E_max = emax.MeV();
-	        double dE    = (E_max-E_min)/nE;
+	        double dE    = (emax-emin)/nE;
+	        double sum1  = 0.0;
 	        double sum2  = 0.0;
-	        double E_obs = E_min;
-	        for (int i = 0; i < nE; ++i) {
-	            double dp_dE = rsp.edisp(GEnergy(E_obs,"MeV"), 0.0, 0.0, 0.0, 0.0, log10_e_src);
-	            sum2        += dp_dE * dE;
-	            E_obs       += dE;
+	        double E_obs = emin;
+	        for (int i = 0; i < nE; ++i, E_obs += dE) {
+                GEnergy ereco(E_obs,"MeV");
+	            double dp1  = (*rsp.edisp())(ereco.log10TeV(), log10_e_src);
+	            double dp2  = rsp.edisp(ereco, 0.0, 0.0, 0.0, 0.0, log10_e_src);
+	            sum1       += dp1 * dE;
+	            sum2       += dp2 * dE;
 	        }
-	        test_value(sum2, 1.0, 0.001, "Energy Dispersion integration for "+eng.print());
+	        GEnergy eng(e_src, "TeV");
+	        test_value(sum1, 1.0, 0.001, "Energy dispersion integration using "+
+                                         rsp.edisp()->classname()+"::operator() "
+                                         "for "+eng.print());
+	        test_value(sum2, 1.0, 0.001, "Energy dispersion integration using "
+                                         "GCTAResponseIrf::edisp() for "+eng.print());
 	    }
  
 	} // endif: energy dispersion was available
@@ -1475,15 +1464,17 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf&amp; rs
 /***********************************************************************//**
  * @brief Utility function for energy dispersion integration tests
  *
- * @param[in] edisp Energy dispersion
- * @param[in] e_src_min Minimum true energy (TeV) (default: 0.1)
- * @param[in] e_src_max Maximum true energy (TeV) (default: 10.0)
+ * @param[in] edisp Energy dispersion instance.
+ * @param[in] e_src_min Minimum true energy (TeV).
+ * @param[in] e_src_max Maximum true energy (TeV).
+ *
+ * Tests whether the energy dispersion integrated over observed energies
+ * gives unity.
  ***************************************************************************/
-void TestGCTAResponse::test_edisp_integration(const GCTAEdisp& edisp,
-                                              const double&    e_src_min,
-                                              const double&    e_src_max)
+void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&   edisp,
+                                              const double&      e_src_min,
+                                              const double&      e_src_max)
 {
-
     // Loop over source energies
     for (double e_src = e_src_min; e_src <= e_src_max; e_src *= 2.0) {
  
@@ -1498,30 +1489,22 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
             continue;
         }
  
-        // Set measured energy range
-        GEnergy  emin     = ebounds.emin();
-        GEnergy  emax     = ebounds.emax();
-        double   logE_min = std::log10(emin.TeV());
-        double   logE_max = std::log10(emax.TeV());
-
-        // Compute step size for numerical integration
-        const int steps = 1000;
-        double    dlogE = (logE_max-logE_min)/steps;
-
         // Perform numerical integration by summing
-        double sum      = 0.0;
-        double logE_obs = logE_min;
-        for (int i = 0; i <= steps; ++i) {
-            double dp_dlogE = edisp(logE_obs, log10_e_src);
-            sum            += dp_dlogE * dlogE;
-            logE_obs       += dlogE;
+        const int nE = 10000;
+        double emin  = ebounds.emin().MeV();
+        double emax  = ebounds.emax().MeV();
+        double dE    = (emax-emin)/nE;
+        double sum   = 0.0;
+        double E_obs = emin;
+        for (int i = 0; i < nE; ++i, E_obs += dE) {
+            GEnergy ereco(E_obs,"MeV");
+            double dp  = edisp(ereco.log10TeV(), log10_e_src);
+            sum       += dp * dE;
         }
-
-        // Set message string
-        GEnergy     eng(e_src, "TeV");
-        std::string msg = edisp.classname()+" integration for "+eng.print();
-        test_value(sum, 1.0, 0.01, msg);
-
+        GEnergy eng(e_src, "TeV");
+        test_value(sum, 1.0, 0.001, "Energy dispersion integration using "+
+                                    edisp.classname()+"::operator() for "+
+                                    eng.print());
     }
  
     // Return
@@ -95,9 +95,9 @@ public:
     void test_response_edisp_integration(const GCTAResponseIrf& rsp,
                                          const double&          e_src_min = 0.1,
                                          const double&          e_src_max = 10.0);
-    void test_edisp_integration(const GCTAEdisp& edisp,
-                                const double&    e_src_min = 0.1,
-                                const double&    e_src_max = 10.0);
+    void test_edisp_integration(const GCTAEdisp&   edisp,
+                                const double&      e_src_min = 0.1,
+                                const double&      e_src_max = 10.0);
 };
  
  
@@ -250,12 +250,12 @@ class Test(gammalib.GPythonTestSuite):
         edisp    = gammalib.GCTAEdispPerfTable(filename)
  
         # Test energy dispersion values
-        self.test_value(edisp(0.0, 0.0), 9.99019627861, 1.0e-6)
-        self.test_value(edisp(0.001, 0.0), 9.9870644077, 1.0e-6)
-        self.test_value(edisp(0.01, 0.0), 9.68182, 1.0e-6)
-        self.test_value(edisp(0.1, 0.0), 0.434382, 1.0e-6)
-        self.test_value(edisp(1.0, 1.0), 18.064868197, 1.0e-6)
-        self.test_value(edisp(1.001, 1.0, 0.0), 18.0463571212, 1.0e-6)
+        self.test_value(edisp(0.0, 0.0),        4.3386871e-06, 1.0e-6)
+        self.test_value(edisp(0.001, 0.0),      4.3273514e-06, 1.0e-6)
+        self.test_value(edisp(0.01, 0.0),       4.1090489e-06, 1.0e-6)
+        self.test_value(edisp(0.1, 0.0),        1.4984969e-07, 1.0e-6)
+        self.test_value(edisp(1.0, 1.0),        7.8454726e-07, 1.0e-6)
+        self.test_value(edisp(1.001, 1.0, 0.0), 7.8194077e-07, 1.0e-6)
  
         # Test GCTAResponseIrf file constructor
         db  = gammalib.GCaldb(self._caldb)
@@ -397,6 +397,7 @@ class Test(gammalib.GPythonTestSuite):
         test_support.pickeling(self, gammalib.GCTAEventCube())
         test_support.pickeling(self, gammalib.GCTAEventList())
         test_support.pickeling(self, gammalib.GCTAInstDir())
+        test_support.pickeling(self, gammalib.GCTAModelBackground())
         test_support.pickeling(self, gammalib.GCTAModelAeffBackground())
         test_support.pickeling(self, gammalib.GCTAModelCubeBackground())
         test_support.pickeling(self, gammalib.GCTAModelIrfBackground())
@@ -404,6 +405,8 @@ class Test(gammalib.GPythonTestSuite):
         test_support.pickeling(self, gammalib.GCTAModelRadialGauss())
         test_support.pickeling(self, gammalib.GCTAModelRadialPolynom())
         test_support.pickeling(self, gammalib.GCTAModelRadialProfile())
+        test_support.pickeling(self, gammalib.GCTAModelSpatialGradient())
+        test_support.pickeling(self, gammalib.GCTAModelSpatialMultiplicative())
         test_support.pickeling(self, gammalib.GCTAObservation())
         test_support.pickeling(self, gammalib.GCTAOnOffObservation())
         test_support.pickeling(self, gammalib.GCTAPointing())
@@ -467,6 +470,9 @@ class Test(gammalib.GPythonTestSuite):
         obs5      = gammalib.GCTAOnOffObservation(obs1, ptsrc, ebds, ebds, regs, regs)
         radgauss  = gammalib.GCTAModelRadialGauss(1.0)
         radacc    = gammalib.GCTAModelRadialAcceptance(radgauss, plaw)
+        multi     = gammalib.GCTAModelSpatialMultiplicative()
+        multi.append(radgauss, 'Src1')
+        multi.append(radgauss, 'Src2')
         kingname  = self._caldb + '/prod1_king.fits'
         hessaeff  = self._data  + '/irf_hess_aeff.fits.gz'
         hesspsf   = self._data  + '/irf_hess_psf.fits.gz'
@@ -498,6 +504,7 @@ class Test(gammalib.GPythonTestSuite):
         test_support.pickeling(self, gammalib.GCTAEventCube(cube))
         test_support.pickeling(self, gammalib.GCTAEventList(list))
         test_support.pickeling(self, gammalib.GCTAInstDir(instdir))
+        test_support.pickeling(self, gammalib.GCTAModelBackground(radgauss,plaw))
         test_support.pickeling(self, gammalib.GCTAModelAeffBackground(plaw))
         test_support.pickeling(self, gammalib.GCTAModelCubeBackground(plaw))
         test_support.pickeling(self, gammalib.GCTAModelIrfBackground(plaw))
@@ -505,6 +512,8 @@ class Test(gammalib.GPythonTestSuite):
         test_support.pickeling(self, gammalib.GCTAModelRadialGauss(1.0))
         test_support.pickeling(self, gammalib.GCTAModelRadialPolynom([1.0,2.0]))
         test_support.pickeling(self, gammalib.GCTAModelRadialProfile(1.0,2.0,3.0))
+        test_support.pickeling(self, gammalib.GCTAModelSpatialGradient(1.0,2.0))
+        test_support.pickeling(self, gammalib.GCTAModelSpatialMultiplicative(multi))
         test_support.pickeling(self, gammalib.GCTAObservation(obs1))
         test_support.pickeling(self, gammalib.GCTAObservation(obs2))
         test_support.pickeling(self, gammalib.GCTAObservation(obs3))
@@ -423,11 +423,11 @@ GEnergy GModelSpectralExponential::mc(const GEnergy&amp; emin,
  * Reads the spectral information from an XML element. The XML element shall
  * have the following format
  *
- *    <spectrum type="Exponential">
- *      <spectrum type="Constant">
- *        <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
- *	    </spectrum>
- *    </spectrum>
+ *     <spectrum type="Exponential">
+ *       <spectrum type="Constant">
+ *         <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
+ *	     </spectrum>
+ *     </spectrum>
  *
  ***************************************************************************/
 void GModelSpectralExponential::read(const GXmlElement& xml)
@@ -469,17 +469,17 @@ void GModelSpectralExponential::read(const GXmlElement&amp; xml)
  * Writes the spectral information into an XML element. The XML element
  * will have the following format:
  *
- *    <spectrum type="Exponential">
- *      <spectrum type="Constant">
- *        <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
- *	    </spectrum>
- *    </spectrum>
+ *     <spectrum type="Exponential">
+ *       <spectrum type="Constant">
+ *         <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
+ *	     </spectrum>
+ *     </spectrum>
  *
  * If no exponential model component is defined the method writes the
  * following XML structure
  *
- *    <spectrum type="Exponential">
- *    </spectrum>
+ *     <spectrum type="Exponential">
+ *     </spectrum>
  *
  ***************************************************************************/
 void GModelSpectralExponential::write(GXmlElement& xml) const
 /***************************************************************************
  *                GResponse.cpp - Abstract response base class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2008-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2008-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -57,6 +57,7 @@
 /* __ Coding definitions _________________________________________________ */
  
 /* __ Debug definitions __________________________________________________ */
+#define G_DEBUG_CONVOLVE_EDISP    //!< Debug convolve for energy dispersion
  
  
 /*==========================================================================
@@ -187,19 +188,19 @@ double GResponse::convolve(const GModelSky&amp;    model,
  
         // Case A: Integration
         if (use_edisp()) {
-    
+
             // Retrieve true energy boundaries
             GEbounds ebounds = this->ebounds(event.energy());
-    
+
             // Loop over all boundaries
             for (int i = 0; i < ebounds.size(); ++i) {
  
-                // Get boundaries in MeV
-                double emin = ebounds.emin(i).MeV();
-                double emax = ebounds.emax(i).MeV();
+                // Get true energy boundaries in log10 MeV
+                double log10_etrue_min = ebounds.emin(i).log10MeV();
+                double log10_etrue_max = ebounds.emax(i).log10MeV();
  
                 // Continue only if valid
-                if (emax > emin) {
+                if (log10_etrue_max > log10_etrue_min) {
  
                     // Setup integration function
                     edisp_kern integrand(this, &obs, &model, &event, srcTime, grad);
@@ -209,10 +210,8 @@ double GResponse::convolve(const GModelSky&amp;    model,
                     integral.fixed_iter(iter);
  
                     // Do Romberg integration
-                    emin  = std::log(emin);
-                    emax  = std::log(emax);
-                    prob += integral.romberg(emin, emax);
-    
+                    prob += integral.romberg(log10_etrue_min, log10_etrue_max);
+
                 } // endif: interval was valid
  
             } // endfor: looped over intervals
@@ -221,7 +220,7 @@ double GResponse::convolve(const GModelSky&amp;    model,
  
         // Case B: No integration (assume no energy dispersion)
         else {
-    
+
             // Get source energy (no dispersion)
             GEnergy srcEng  = event.energy();
  
@@ -410,39 +409,34 @@ double GResponse::eval_prob(const GModelSky&amp;    model,
  
  
 /***********************************************************************//**
- * @brief Integration kernel for edisp_kern() method
+ * @brief Integration kernel for GResponse::edisp_kern() class
  *
- * @param[in] x Function value.
+ * @param[in] logEsrc Base 10 logarithm of true photon energy in MeV.
  *
- * This method implements the integration kernel needed for the edisp_kern()
- * method.
+ * This method implements the integration kernel needed for the
+ * GResponse::edisp_kern() class.
  ***************************************************************************/
-double GResponse::edisp_kern::eval(const double& x)
+double GResponse::edisp_kern::eval(const double& logEsrc)
 {
-    // Set energy
-    GEnergy eng;
-    double expx = std::exp(x);
-    eng.MeV(expx);
+    // Set true energy
+    GEnergy srcEng;
+    srcEng.log10MeV(logEsrc);
  
     // Get function value
-    double value = m_parent->eval_prob(*m_model, *m_event, eng, m_srcTime, *m_obs, m_grad);
-
-    // Save value if needed
-    #if defined(G_NAN_CHECK)
-    double value_out = value;
-    #endif
+    double value = m_parent->eval_prob(*m_model, *m_event,
+                                       srcEng, m_srcTime, *m_obs, m_grad);
  
-    // Correct for variable substitution
-    value *= expx;
+    // Account for the fact that integration is done in log10 MeV and
+    // not in MeV
+    value *= gammalib::ln10 * m_event->energy().MeV();
  
     // Compile option: Check for NaN
     #if defined(G_NAN_CHECK)
     if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
         std::cout << "*** ERROR: GResponse::edisp_kern::eval";
-        std::cout << "(x=" << x << "): ";
+        std::cout << "(logEsrc=" << logEsrc << "): ";
         std::cout << " NaN/Inf encountered";
-        std::cout << " (value=" << value_out;
-        std::cout << " exp(x)=" << expx;
+        std::cout << " (value=" << value;
         std::cout << ")" << std::endl;
     }
     #endif