diff --git a/NEWS b/NEWS
index 029e6be..e01fc9c 100644
--- a/NEWS
+++ b/NEWS
@@ -142,7 +142,8 @@ The following methods have been renamed:
 The arguments for the following methods have been changed:
 - GCTAOnOffObservation(GCTAObservation&, GModelSpatial&, GEbounds&, GEbounds&, GSkyRegions&, GSkyRegions&, bool&)
 - GCTAModelRadial::mc(GRan&)
-- GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&)  
+- GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&)
+- GCTAEdisp::mc(GRan&, GEnergy&, double&, double&, double&, double&)
 - GCTAResponseIrf::edisp(GEnergy&, GEnergy&, double&, double&, double&, double&)
 - GLATAeff(GFits&, std::string&)
 - GLATEdisp(GFits&, std::string&)
diff --git a/inst/cta/include/GCTAEdisp.hpp b/inst/cta/include/GCTAEdisp.hpp
index 142bab1..b19591c 100644
--- a/inst/cta/include/GCTAEdisp.hpp
+++ b/inst/cta/include/GCTAEdisp.hpp
@@ -71,12 +71,12 @@ public:
     virtual std::string classname(void) const = 0;
     virtual void        load(const GFilename& filename) = 0;
     virtual GFilename   filename(void) const = 0;
-    virtual GEnergy     mc(GRan&         ran,
-                           const double& logE,
-                           const double& theta = 0.0,
-                           const double& phi = 0.0,
-                           const double& zenith = 0.0,
-                           const double& azimuth = 0.0) const = 0;
+    virtual GEnergy     mc(GRan&          ran,
+                           const GEnergy& etrue,
+                           const double&  theta = 0.0,
+                           const double&  phi = 0.0,
+                           const double&  zenith = 0.0,
+                           const double&  azimuth = 0.0) const = 0;
     virtual GEbounds    ebounds_obs(const double& logEsrc,
                                     const double& theta = 0.0,
                                     const double& phi = 0.0,
diff --git a/inst/cta/include/GCTAEdisp2D.hpp b/inst/cta/include/GCTAEdisp2D.hpp
index 9401e64..2ce8133 100644
--- a/inst/cta/include/GCTAEdisp2D.hpp
+++ b/inst/cta/include/GCTAEdisp2D.hpp
@@ -119,12 +119,12 @@ public:
     std::string  classname(void) const;
     void         load(const GFilename& filename);
     GFilename    filename(void) const;
-    GEnergy      mc(GRan&         ran,
-                    const double& logE,
-                    const double& theta = 0.0,
-                    const double& phi = 0.0,
-                    const double& zenith = 0.0,
-                    const double& azimuth = 0.0) const;
+    GEnergy      mc(GRan&          ran,
+                    const GEnergy& etrue,
+                    const double&  theta = 0.0,
+                    const double&  phi = 0.0,
+                    const double&  zenith = 0.0,
+                    const double&  azimuth = 0.0) const;
     GEbounds     ebounds_obs(const double& logEsrc,
                              const double& theta = 0.0,
                              const double& phi = 0.0,
diff --git a/inst/cta/include/GCTAEdispPerfTable.hpp b/inst/cta/include/GCTAEdispPerfTable.hpp
index 6737738..24825f8 100644
--- a/inst/cta/include/GCTAEdispPerfTable.hpp
+++ b/inst/cta/include/GCTAEdispPerfTable.hpp
@@ -88,12 +88,12 @@ public:
     std::string         classname(void) const;
     void                load(const GFilename& filename);
     GFilename           filename(void) const;
-    GEnergy             mc(GRan&         ran,
-                           const double& logE,
-                           const double& theta = 0.0,
-                           const double& phi = 0.0,
-                           const double& zenith = 0.0,
-                           const double& azimuth = 0.0) const;
+    GEnergy             mc(GRan&          ran,
+                           const GEnergy& etrue,
+                           const double&  theta = 0.0,
+                           const double&  phi = 0.0,
+                           const double&  zenith = 0.0,
+                           const double&  azimuth = 0.0) const;
     GEbounds            ebounds_obs(const double& logEsrc,
                                     const double& theta = 0.0,
                                     const double& phi = 0.0,
diff --git a/inst/cta/include/GCTAEdispRmf.hpp b/inst/cta/include/GCTAEdispRmf.hpp
index d8d654d..a5a93a9 100644
--- a/inst/cta/include/GCTAEdispRmf.hpp
+++ b/inst/cta/include/GCTAEdispRmf.hpp
@@ -84,11 +84,11 @@ public:
     void          load(const GFilename& filename);
     GFilename     filename(void) const;
     GEnergy       mc(GRan& ran,
-                     const double& logEsrc,
-                     const double& theta = 0.0,
-                     const double& phi = 0.0,
-                     const double& zenith = 0.0,
-                     const double& azimuth = 0.0) const;
+                     const GEnergy& etrue,
+                     const double&  theta = 0.0,
+                     const double&  phi = 0.0,
+                     const double&  zenith = 0.0,
+                     const double&  azimuth = 0.0) const;
     GEbounds      ebounds_obs(const double& logEsrc,
                               const double& theta = 0.0,
                               const double& phi = 0.0,
diff --git a/inst/cta/pyext/GCTAEdisp.i b/inst/cta/pyext/GCTAEdisp.i
index ed0a240..69e3f86 100644
--- a/inst/cta/pyext/GCTAEdisp.i
+++ b/inst/cta/pyext/GCTAEdisp.i
@@ -56,12 +56,12 @@ public:
     virtual std::string classname(void) const = 0;
     virtual void        load(const GFilename& filename) = 0;
     virtual GFilename   filename(void) const = 0;
-    virtual GEnergy     mc(GRan&         ran,
-                           const double& logE,
-                           const double& theta = 0.0,
-                           const double& phi = 0.0,
-                           const double& zenith = 0.0,
-                           const double& azimuth = 0.0) const = 0;
+    virtual GEnergy     mc(GRan&          ran,
+                           const GEnergy& etrue,
+                           const double&  theta = 0.0,
+                           const double&  phi = 0.0,
+                           const double&  zenith = 0.0,
+                           const double&  azimuth = 0.0) const = 0;
     virtual GEbounds    ebounds_obs(const double& logEsrc,
                                     const double& theta = 0.0,
                                     const double& phi = 0.0,
diff --git a/inst/cta/pyext/GCTAEdisp2D.i b/inst/cta/pyext/GCTAEdisp2D.i
index 29a7419..0787a88 100644
--- a/inst/cta/pyext/GCTAEdisp2D.i
+++ b/inst/cta/pyext/GCTAEdisp2D.i
@@ -57,12 +57,12 @@ public:
     std::string  classname(void) const;
     void         load(const GFilename& filename);
     GFilename    filename(void) const;
-    GEnergy      mc(GRan&         ran,
-                    const double& logE,
-                    const double& theta = 0.0,
-                    const double& phi = 0.0,
-                    const double& zenith = 0.0,
-                    const double& azimuth = 0.0) const;
+    GEnergy      mc(GRan&          ran,
+                    const GEnergy& etrue,
+                    const double&  theta = 0.0,
+                    const double&  phi = 0.0,
+                    const double&  zenith = 0.0,
+                    const double&  azimuth = 0.0) const;
     GEbounds     ebounds_obs(const double& logEsrc,
                              const double& theta = 0.0,
                              const double& phi = 0.0,
diff --git a/inst/cta/pyext/GCTAEdispPerfTable.i b/inst/cta/pyext/GCTAEdispPerfTable.i
index 59563b9..8f0eea8 100644
--- a/inst/cta/pyext/GCTAEdispPerfTable.i
+++ b/inst/cta/pyext/GCTAEdispPerfTable.i
@@ -57,12 +57,12 @@ public:
     std::string         classname(void) const;
     void                load(const GFilename& filename);
     GFilename           filename(void) const;
-    GEnergy             mc(GRan&         ran,
-                           const double& logE,
-                           const double& theta = 0.0,
-                           const double& phi = 0.0,
-                           const double& zenith = 0.0,
-                           const double& azimuth = 0.0) const;
+    GEnergy             mc(GRan&          ran,
+                           const GEnergy& etrue,
+                           const double&  theta = 0.0,
+                           const double&  phi = 0.0,
+                           const double&  zenith = 0.0,
+                           const double&  azimuth = 0.0) const;
     GEbounds            ebounds_obs(const double& logEsrc,
                                     const double& theta = 0.0,
                                     const double& phi = 0.0,
diff --git a/inst/cta/pyext/GCTAEdispRmf.i b/inst/cta/pyext/GCTAEdispRmf.i
index 4a30602..6c76d08 100644
--- a/inst/cta/pyext/GCTAEdispRmf.i
+++ b/inst/cta/pyext/GCTAEdispRmf.i
@@ -57,12 +57,12 @@ public:
     std::string   classname(void) const;
     void          load(const GFilename& filename);
     GFilename     filename(void) const;
-    GEnergy       mc(GRan& ran,
-                     const double& logE,
-                     const double& theta = 0.0,
-                     const double& phi = 0.0,
-                     const double& zenith = 0.0,
-                     const double& azimuth = 0.0) const;
+    GEnergy       mc(GRan&          ran,
+                     const GEnergy& etrue,
+                     const double&  theta = 0.0,
+                     const double&  phi = 0.0,
+                     const double&  zenith = 0.0,
+                     const double&  azimuth = 0.0) const;
     GEbounds      ebounds_obs(const double& logEsrc,
                               const double& theta = 0.0,
                               const double& phi = 0.0,
diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp
index 8160663..e4525a9 100644
--- a/inst/cta/src/GCTAEdisp2D.cpp
+++ b/inst/cta/src/GCTAEdisp2D.cpp
@@ -467,7 +467,7 @@ void GCTAEdisp2D::save(const GFilename& filename, const bool& clobber) const
  * @brief Simulate energy dispersion
  *
  * @param[in] ran Random number generator.
- * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] etrue True photon energy.
  * @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.
@@ -482,23 +482,21 @@ void GCTAEdisp2D::save(const GFilename& filename, const bool& clobber) const
  * angle @p theta. If no energy dispersion information is available the
  * method will return the true photon energy.
  ***************************************************************************/
-GEnergy GCTAEdisp2D::mc(GRan&         ran,
-                        const double& logEsrc,
-                        const double& theta,
-                        const double& phi,
-                        const double& zenith,
-                        const double& azimuth) const
+GEnergy GCTAEdisp2D::mc(GRan&          ran,
+                        const GEnergy& etrue,
+                        const double&  theta,
+                        const double&  phi,
+                        const double&  zenith,
+                        const double&  azimuth) const
 {
     // Make sure that energy dispersion is online
     fetch();
 
-    // Initialise energies
-    GEnergy etrue;
+    // Initialise reconstructed photon energy
     GEnergy ereco;
-    etrue.log10TeV(logEsrc);
 
     // Get boundaries for observed energy
-    GEbounds ebounds = ebounds_obs(logEsrc, theta, phi, zenith, azimuth);
+    GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth);
     double   emin    = ebounds.emin().log10TeV();
     double   emax    = ebounds.emax().log10TeV();
 
@@ -506,9 +504,8 @@ GEnergy GCTAEdisp2D::mc(GRan&         ran,
     // energy (this means that no energy dispersion information is available
     // for that energy)
     if (emin >= emax) {
-        double      eng = std::pow(10.0, logEsrc);
         std::string msg = "No valid energy dispersion information available "
-                          "for true photon energy "+gammalib::str(eng)+" TeV,"
+                          "for true photon energy "+etrue.print()+","
                           "offset angle "+
                           gammalib::str(theta*gammalib::rad2deg)+" deg "
                           "and azimuth angle "+
@@ -527,7 +524,7 @@ GEnergy GCTAEdisp2D::mc(GRan&         ran,
 
     // Initialise rejection method
     double    ewidth               = emax - emin;
-    double    logEobs              = logEsrc;
+    double    logEobs              = etrue.log10TeV();
     double    f                    = 0.0;
     double    ftest                = 1.0;
     int       zeros                = 0;
@@ -549,10 +546,9 @@ GEnergy GCTAEdisp2D::mc(GRan&         ran,
             if (f == 0.0) {
                 zeros++;
                 if (zeros > max_subsequent_zeros) {
-                    double      eng = std::pow(10.0, logEsrc);
                     std::string msg = "No valid energy dispersion information "
                                       "available  for true photon energy "+
-                                      gammalib::str(eng)+" TeV, offset angle "+
+                                      etrue.print()+", offset angle "+
                                       gammalib::str(theta*gammalib::rad2deg)+
                                       " deg and azimuth angle "+
                                       gammalib::str(phi*gammalib::rad2deg)+
diff --git a/inst/cta/src/GCTAEdispPerfTable.cpp b/inst/cta/src/GCTAEdispPerfTable.cpp
index 7503f24..1979e05 100644
--- a/inst/cta/src/GCTAEdispPerfTable.cpp
+++ b/inst/cta/src/GCTAEdispPerfTable.cpp
@@ -322,7 +322,7 @@ void GCTAEdispPerfTable::load(const GFilename& filename)
  * @brief Simulate energy dispersion
  *
  * @param[in] ran Random number generator.
- * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] etrue True photon energy.
  * @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.
@@ -332,13 +332,16 @@ void GCTAEdispPerfTable::load(const GFilename& filename)
  * Draws observed energy value from a normal distribution of width
  * m_par_sigma around @p logE.
  ***************************************************************************/
-GEnergy GCTAEdispPerfTable::mc(GRan&         ran,
-                               const double& logEsrc,
-                               const double& theta,
-                               const double& phi,
-                               const double& zenith,
-                               const double& azimuth) const
+GEnergy GCTAEdispPerfTable::mc(GRan&          ran,
+                               const GEnergy& etrue,
+                               const double&  theta,
+                               const double&  phi,
+                               const double&  zenith,
+                               const double&  azimuth) const
 {
+    // Get log10 TeV of true photon energy
+    double logEsrc = etrue.log10TeV();
+
     // Update the parameter cache
     update(logEsrc);
 
diff --git a/inst/cta/src/GCTAEdispRmf.cpp b/inst/cta/src/GCTAEdispRmf.cpp
index a26f780..0e9f7d5 100644
--- a/inst/cta/src/GCTAEdispRmf.cpp
+++ b/inst/cta/src/GCTAEdispRmf.cpp
@@ -287,7 +287,7 @@ void GCTAEdispRmf::load(const GFilename& filename)
  * @brief Simulate energy dispersion
  *
  * @param[in] ran Random number generator.
- * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
+ * @param[in] etrue True photon energy.
  * @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.
@@ -296,20 +296,18 @@ void GCTAEdispRmf::load(const GFilename& filename)
  *
  * Draws observed energy value from RMF matrix.
  ***************************************************************************/
-GEnergy GCTAEdispRmf::mc(GRan&         ran,
-                         const double& logEsrc,
-                         const double& theta,
-                         const double& phi,
-                         const double& zenith,
-                         const double& azimuth) const
+GEnergy GCTAEdispRmf::mc(GRan&          ran,
+                         const GEnergy& etrue,
+                         const double&  theta,
+                         const double&  phi,
+                         const double&  zenith,
+                         const double&  azimuth) const
 {
     // Initialise energies
-    GEnergy etrue;
     GEnergy ereco;
-    etrue.log10TeV(logEsrc);
 
     // Get boundaries for rejection method
-    GEbounds ebounds = ebounds_obs(logEsrc, theta, phi, zenith, azimuth);
+    GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth);
     double   emin    = ebounds.emin().log10TeV();
     double   emax    = ebounds.emax().log10TeV();
 
diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp
index 6ca0902..620d5c8 100644
--- a/inst/cta/src/GCTAResponseIrf.cpp
+++ b/inst/cta/src/GCTAResponseIrf.cpp
@@ -578,8 +578,10 @@ GEbounds GCTAResponseIrf::ebounds(const GEnergy& obsEnergy) const
  * @todo Set polar angle phi of photon in camera system
  * @todo Implement energy dispersion
  ***************************************************************************/
-GCTAEventAtom* GCTAResponseIrf::mc(const double& area, const GPhoton& photon,
-                                   const GObservation& obs, GRan& ran) const
+GCTAEventAtom* GCTAResponseIrf::mc(const double&       area,
+                                   const GPhoton&      photon,
+                                   const GObservation& obs,
+                                   GRan&               ran) const
 {
     // Initialise event
     GCTAEventAtom* event = NULL;
@@ -639,7 +641,8 @@ GCTAEventAtom* GCTAResponseIrf::mc(const double& area, const GPhoton& photon,
                 // Set measured photon energy
                 GEnergy energy = photon.energy();
                 if (use_edisp()) {
-                    energy = edisp()->mc(ran, srcLogEng, theta, phi, zenith, azimuth);
+                    energy = edisp()->mc(ran, photon.energy(), theta, phi,
+                                         zenith, azimuth);
                 }
 
                 // Allocate event
diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp
index e10a8e9..25f0b9c 100644
--- a/inst/cta/test/test_CTA.cpp
+++ b/inst/cta/test/test_CTA.cpp
@@ -917,7 +917,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
 
     // Test mc method
     GRan ran;
-    test_value(edisp2.mc(ran, 0.0).TeV(), 0.88667418590027);
+    test_value(edisp2.mc(ran, etrue).TeV(), 0.88667418590027);
 
     // Test normalisation
     test_edisp_integration(edisp2);
@@ -973,7 +973,7 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
 
     // Test mc method
     GRan ran;
-    test_value(edisp2.mc(ran, 0.0).TeV(), 0.930126476703161);
+    test_value(edisp2.mc(ran, etrue).TeV(), 0.930126476703161);
 
     // Test normalisation
     test_edisp_integration(edisp2, 0.5, 10.0);
@@ -1028,7 +1028,7 @@ void TestGCTAResponse::test_response_edisp_2D(void)
 
     // Test mc method
     GRan ran;
-    test_value(edisp2.mc(ran, 0.0).TeV(), 0.872387491302767);
+    test_value(edisp2.mc(ran, etrue).TeV(), 0.872387491302767);
 
     // Test normalisation
     test_edisp_integration(edisp2);