diff --git a/ChangeLog b/ChangeLog
index f3599a6..cef5278 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -3,6 +3,7 @@
         * Version 1.6.0 released
         ========================
 
+        Revise GCTAEdisp class interface (#2704) 
         GCTAEdisp returns now energy dispersion values per MeV (#2704) 
         Protect GCTAInstDir from access to invalid information (#2694)
         Implement general CTA background model (#2694)
diff --git a/NEWS b/NEWS
index e01fc9c..787f7dd 100644
--- a/NEWS
+++ b/NEWS
@@ -138,6 +138,8 @@ The following methods have been removed:
 
 The following methods have been renamed:
 - GPhoton::mc_id() from GPhoton::mcid()
+- GCTAEdisp::ereco_bounds() from GCTAEdisp::ebounds_obs()
+- GCTAEdisp::etrue_bounds() from GCTAEdisp::ebounds_src()
 
 The arguments for the following methods have been changed:
 - GCTAOnOffObservation(GCTAObservation&, GModelSpatial&, GEbounds&, GEbounds&, GSkyRegions&, GSkyRegions&, bool&)
@@ -310,6 +312,10 @@ 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 GCTAEdisp class now uses GEnergy objects instead of log10(E) to deal with
+true and reconstructed energies. This avoids computational overhead when going
+from one unit to another (#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-
diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst
index 19fa1cf..77e785d 100644
--- a/doc/source/admin/release_history/1.6.rst
+++ b/doc/source/admin/release_history/1.6.rst
@@ -46,7 +46,8 @@ Improvements
 ------------
 
 * [`2704 <https://cta-redmine.irap.omp.eu/issues/2704>`_] -
-  ``GCTAEdisp`` returns now energy dispersion values per MeV
+  ``GCTAEdisp`` returns now energy dispersion values per MeV.
+  The ``GCTAEdisp`` interface was revised.
 * [`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>`_] -
diff --git a/inst/cta/include/GCTAEdisp.hpp b/inst/cta/include/GCTAEdisp.hpp
index b19591c..ccb4f83 100644
--- a/inst/cta/include/GCTAEdisp.hpp
+++ b/inst/cta/include/GCTAEdisp.hpp
@@ -77,16 +77,16 @@ public:
                            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,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const = 0;
-    virtual GEbounds    ebounds_src(const double& logEobs,
-                                    const double& theta = 0.0,
-                                    const double& phi = 0.0,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const = 0;
+    virtual GEbounds    ereco_bounds(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    etrue_bounds(const GEnergy& ereco,
+                                     const double&  theta = 0.0,
+                                     const double&  phi = 0.0,
+                                     const double&  zenith = 0.0,
+                                     const double&  azimuth = 0.0) const = 0;
     virtual double      prob_erecobin(const GEnergy& ereco_min,
                                       const GEnergy& ereco_max,
                                       const GEnergy& etrue,
diff --git a/inst/cta/include/GCTAEdisp2D.hpp b/inst/cta/include/GCTAEdisp2D.hpp
index 2ce8133..f988ef2 100644
--- a/inst/cta/include/GCTAEdisp2D.hpp
+++ b/inst/cta/include/GCTAEdisp2D.hpp
@@ -125,16 +125,16 @@ public:
                     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,
-                             const double& zenith = 0.0,
-                             const double& azimuth = 0.0) const;
-    GEbounds     ebounds_src(const double& logEobs,
-                             const double& theta = 0.0,
-                             const double& phi = 0.0,
-                             const double& zenith = 0.0,
-                             const double& azimuth = 0.0) const;
+    GEbounds     ereco_bounds(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     etrue_bounds(const GEnergy& ereco,
+                              const double&  theta = 0.0,
+                              const double&  phi = 0.0,
+                              const double&  zenith = 0.0,
+                              const double&  azimuth = 0.0) const;
     double       prob_erecobin(const GEnergy& ereco_min,
                                const GEnergy& ereco_max,
                                const GEnergy& etrue,
@@ -248,8 +248,8 @@ private:
     mutable bool                  m_ebounds_src_computed;
     mutable double                m_last_theta_obs;
     mutable double                m_last_theta_src;
-    mutable double                m_last_logEsrc;
-    mutable double                m_last_logEobs;
+    mutable GEnergy               m_last_etrue;
+    mutable GEnergy               m_last_ereco;
     mutable int                   m_index_obs;
     mutable int                   m_index_src;
     mutable double                m_max_edisp;
diff --git a/inst/cta/include/GCTAEdispPerfTable.hpp b/inst/cta/include/GCTAEdispPerfTable.hpp
index 24825f8..75595ef 100644
--- a/inst/cta/include/GCTAEdispPerfTable.hpp
+++ b/inst/cta/include/GCTAEdispPerfTable.hpp
@@ -94,16 +94,16 @@ public:
                            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,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const;
-    GEbounds            ebounds_src(const double& logEobs,
-                                    const double& theta = 0.0,
-                                    const double& phi = 0.0,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const;
+    GEbounds            ereco_bounds(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            etrue_bounds(const GEnergy& ereco,
+                                     const double&  theta = 0.0,
+                                     const double&  phi = 0.0,
+                                     const double&  zenith = 0.0,
+                                     const double&  azimuth = 0.0) const;
     double              prob_erecobin(const GEnergy& ereco_min,
                                       const GEnergy& ereco_max,
                                       const GEnergy& etrue,
diff --git a/inst/cta/include/GCTAEdispRmf.hpp b/inst/cta/include/GCTAEdispRmf.hpp
index a5a93a9..6ea2741 100644
--- a/inst/cta/include/GCTAEdispRmf.hpp
+++ b/inst/cta/include/GCTAEdispRmf.hpp
@@ -89,16 +89,16 @@ public:
                      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,
-                              const double& zenith = 0.0,
-                              const double& azimuth = 0.0) const;
-    GEbounds      ebounds_src(const double& logEobs,
-                              const double& theta = 0.0,
-                              const double& phi = 0.0,
-                              const double& zenith = 0.0,
-                              const double& azimuth = 0.0) const;
+    GEbounds      ereco_bounds(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      etrue_bounds(const GEnergy& ereco,
+                               const double&  theta = 0.0,
+                               const double&  phi = 0.0,
+                               const double&  zenith = 0.0,
+                               const double&  azimuth = 0.0) const;
     double        prob_erecobin(const GEnergy& ereco_min,
                                 const GEnergy& ereco_max,
                                 const GEnergy& etrue,
@@ -149,8 +149,8 @@ private:
     mutable double                m_max_edisp;
     mutable double                m_last_theta_obs;
     mutable double                m_last_theta_src;
-    mutable double                m_last_logEsrc;
-    mutable double                m_last_logEobs;
+    mutable GEnergy               m_last_etrue_ereco_bounds;
+    mutable GEnergy               m_last_ereco_etrue_bounds;
     mutable int                   m_index_obs;
     mutable int                   m_index_src;
     mutable bool                  m_ebounds_obs_computed;
diff --git a/inst/cta/pyext/GCTAEdisp.i b/inst/cta/pyext/GCTAEdisp.i
index 69e3f86..b6628c7 100644
--- a/inst/cta/pyext/GCTAEdisp.i
+++ b/inst/cta/pyext/GCTAEdisp.i
@@ -62,16 +62,16 @@ public:
                            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,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const = 0;
-    virtual GEbounds    ebounds_src(const double& logEobs,
-                                    const double& theta = 0.0,
-                                    const double& phi = 0.0,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const = 0;
+    virtual GEbounds    ereco_bounds(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    etrue_bounds(const GEnergy& ereco,
+                                     const double&  theta = 0.0,
+                                     const double&  phi = 0.0,
+                                     const double&  zenith = 0.0,
+                                     const double&  azimuth = 0.0) const = 0;
     virtual double      prob_erecobin(const GEnergy& ereco_min,
                                       const GEnergy& ereco_max,
                                       const GEnergy& etrue,
diff --git a/inst/cta/pyext/GCTAEdisp2D.i b/inst/cta/pyext/GCTAEdisp2D.i
index 0787a88..fae8f48 100644
--- a/inst/cta/pyext/GCTAEdisp2D.i
+++ b/inst/cta/pyext/GCTAEdisp2D.i
@@ -63,16 +63,16 @@ public:
                     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,
-                             const double& zenith = 0.0,
-                             const double& azimuth = 0.0) const;
-    GEbounds     ebounds_src(const double& logEobs,
-                             const double& theta = 0.0,
-                             const double& phi = 0.0,
-                             const double& zenith = 0.0,
-                             const double& azimuth = 0.0) const;
+    GEbounds     ereco_bounds(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     etrue_bounds(const GEnergy& ereco,
+                              const double&  theta = 0.0,
+                              const double&  phi = 0.0,
+                              const double&  zenith = 0.0,
+                              const double&  azimuth = 0.0) const;
     double       prob_erecobin(const GEnergy& ereco_min,
                                const GEnergy& ereco_max,
                                const GEnergy& etrue,
diff --git a/inst/cta/pyext/GCTAEdispPerfTable.i b/inst/cta/pyext/GCTAEdispPerfTable.i
index 8f0eea8..8b904f1 100644
--- a/inst/cta/pyext/GCTAEdispPerfTable.i
+++ b/inst/cta/pyext/GCTAEdispPerfTable.i
@@ -63,16 +63,16 @@ public:
                            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,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const;
-    GEbounds            ebounds_src(const double& logEobs,
-                                    const double& theta = 0.0,
-                                    const double& phi = 0.0,
-                                    const double& zenith = 0.0,
-                                    const double& azimuth = 0.0) const;
+    GEbounds            ereco_bounds(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            etrue_bounds(const GEnergy& ereco,
+                                     const double&  theta = 0.0,
+                                     const double&  phi = 0.0,
+                                     const double&  zenith = 0.0,
+                                     const double&  azimuth = 0.0) const;
     double              prob_erecobin(const GEnergy& ereco_min,
                                       const GEnergy& ereco_max,
                                       const GEnergy& etrue,
diff --git a/inst/cta/pyext/GCTAEdispRmf.i b/inst/cta/pyext/GCTAEdispRmf.i
index 6c76d08..d1305e0 100644
--- a/inst/cta/pyext/GCTAEdispRmf.i
+++ b/inst/cta/pyext/GCTAEdispRmf.i
@@ -63,16 +63,16 @@ public:
                      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,
-                              const double& zenith = 0.0,
-                              const double& azimuth = 0.0) const;
-    GEbounds      ebounds_src(const double& logEobs,
-                              const double& theta = 0.0,
-                              const double& phi = 0.0,
-                              const double& zenith = 0.0,
-                              const double& azimuth = 0.0) const;
+    GEbounds      ereco_bounds(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      etrue_bounds(const GEnergy& ereco,
+                               const double&  theta = 0.0,
+                               const double&  phi = 0.0,
+                               const double&  zenith = 0.0,
+                               const double&  azimuth = 0.0) const;
     double        prob_erecobin(const GEnergy& ereco_min,
                                 const GEnergy& ereco_max,
                                 const GEnergy& etrue,
diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp
index e4525a9..afea25d 100644
--- a/inst/cta/src/GCTAEdisp2D.cpp
+++ b/inst/cta/src/GCTAEdisp2D.cpp
@@ -496,7 +496,7 @@ GEnergy GCTAEdisp2D::mc(GRan&          ran,
     GEnergy ereco;
 
     // Get boundaries for observed energy
-    GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth);
+    GEbounds ebounds = ereco_bounds(etrue, theta, phi, zenith, azimuth);
     double   emin    = ebounds.emin().log10TeV();
     double   emax    = ebounds.emax().log10TeV();
 
@@ -576,7 +576,7 @@ GEnergy GCTAEdisp2D::mc(GRan&          ran,
 /***********************************************************************//**
  * @brief Return observed energy interval that contains the energy dispersion.
  *
- * @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.
@@ -587,11 +587,11 @@ GEnergy GCTAEdisp2D::mc(GRan&          ran,
  * dispersion becomes negligible for a given true energy @p logEsrc and
  * offset angle @p theta.
  ***************************************************************************/
-GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc,
-                                  const double& theta,
-                                  const double& phi,
-                                  const double& zenith,
-                                  const double& azimuth) const
+GEbounds GCTAEdisp2D::ereco_bounds(const GEnergy& etrue,
+                                   const double&  theta,
+                                   const double&  phi,
+                                   const double&  zenith,
+                                   const double&  azimuth) const
 {
     // Make sure that energy dispersion is online
     fetch();
@@ -602,18 +602,21 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc,
         // Set computation flag
         m_ebounds_obs_computed = true;
         m_last_theta_obs       = theta;
-        m_last_logEsrc         = -30.0; // force update
+        m_last_etrue.log10TeV(-30.0); // force update
 
         // Compute ebounds_obs
         compute_ebounds_obs(theta, phi, zenith, azimuth);
 
     }
 
-    // Search index only if logEsrc has changed
-    if (logEsrc != m_last_logEsrc) {
+    // Search index only if true photon energy has changed
+    if (etrue != m_last_etrue) {
+
+        // Store last true photon energy
+        m_last_etrue = etrue;
 
-        // Store last log(Esrc)
-        m_last_logEsrc = logEsrc;
+        // Get true energy in TeV
+        double etrue_TeV = etrue.TeV();
 
         // Find right index with bisection
         int low  = 0;
@@ -621,7 +624,7 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc,
         while ((high-low) > 1) {
             int  mid = (low+high) / 2;
             double e = m_edisp.axis_lo(m_inx_etrue, mid);
-            if (logEsrc < std::log10(e)) {
+            if (etrue_TeV < e) {
                 high = mid;
             }
             else {
@@ -642,7 +645,7 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc,
 /***********************************************************************//**
  * @brief Return true energy interval that contains the energy dispersion.
  *
- * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
+ * @param[in] ereco Reconstructed event 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.
@@ -653,11 +656,11 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc,
  * dispersion becomes negligible for a given observed energy @p logEobs and
  * offset angle @p theta.
  ***************************************************************************/
-GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs,
-                                  const double& theta,
-                                  const double& phi,
-                                  const double& zenith,
-                                  const double& azimuth) const
+GEbounds GCTAEdisp2D::etrue_bounds(const GEnergy& ereco,
+                                   const double&  theta,
+                                   const double&  phi,
+                                   const double&  zenith,
+                                   const double&  azimuth) const
 {
     // Make sure that energy dispersion is online
     fetch();
@@ -668,18 +671,21 @@ GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs,
         // Set computation flag
         m_ebounds_src_computed = true;
         m_last_theta_src       = theta;
-        m_last_logEobs         = -30.0; // force update
+        m_last_ereco.log10TeV(-30.0); // force update
 
         // Compute ebounds_src
         compute_ebounds_src(theta, phi, zenith, azimuth);
 
     }
 
-    // Search index only if logEobs has changed
-    if (logEobs != m_last_logEobs) {
+    // Search index only if reconstructed event energy has changed
+    if (ereco != m_last_ereco) {
+
+        // Store reconstructed event energy
+        m_last_ereco = ereco;
 
-        // Store observed energy
-        m_last_logEobs = logEobs;
+        // Get reconstructed event energy in TeV
+        double ereco_TeV = ereco.TeV();
 
         // Find right index with bisection
         int low  = 0;
@@ -687,7 +693,7 @@ GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs,
         while ((high-low) > 1) {
             int  mid = (low+high) / 2;
             double e = m_edisp.axis_lo(m_inx_etrue, mid);
-            if (logEobs < std::log10(e)) {
+            if (ereco_TeV < e) {
                 high = mid;
             }
             else {
@@ -1046,8 +1052,8 @@ void GCTAEdisp2D::init_members(void)
     m_ebounds_src_computed = false;
     m_last_theta_obs       =  -1.0;
     m_last_theta_src       =  -1.0;
-    m_last_logEsrc         = -30.0;
-    m_last_logEobs         = -30.0;
+    m_last_etrue.log10TeV(-30.0);
+    m_last_ereco.log10TeV(-30.0);
     m_index_obs            =     0;
     m_index_src            =     0;
     m_max_edisp            =   0.0;
@@ -1086,8 +1092,8 @@ void GCTAEdisp2D::copy_members(const GCTAEdisp2D& edisp)
     m_ebounds_src_computed = edisp.m_ebounds_src_computed;
     m_last_theta_obs       = edisp.m_last_theta_obs;
     m_last_theta_src       = edisp.m_last_theta_src;
-    m_last_logEsrc         = edisp.m_last_logEsrc;
-    m_last_logEobs         = edisp.m_last_logEobs;
+    m_last_etrue           = edisp.m_last_etrue;
+    m_last_ereco           = edisp.m_last_ereco;
     m_index_obs            = edisp.m_index_obs;
     m_index_src            = edisp.m_index_src;
     m_max_edisp            = edisp.m_max_edisp;
@@ -1532,7 +1538,7 @@ void GCTAEdisp2D::normalize_table(void)
                 double sum = 0.0;
 
                 // Get integration boundaries
-                GEbounds ebounds = ebounds_obs(logEsrc, theta);
+                GEbounds ebounds = ereco_bounds(etrue, theta);
 
                 // Loop over all energy intervals
                 for (int i = 0; i < ebounds.size(); ++i) {
diff --git a/inst/cta/src/GCTAEdispPerfTable.cpp b/inst/cta/src/GCTAEdispPerfTable.cpp
index 1979e05..0396400 100644
--- a/inst/cta/src/GCTAEdispPerfTable.cpp
+++ b/inst/cta/src/GCTAEdispPerfTable.cpp
@@ -360,7 +360,7 @@ GEnergy GCTAEdispPerfTable::mc(GRan&          ran,
 /***********************************************************************//**
  * @brief Return observed energy interval that contains the energy dispersion.
  *
- * @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.
@@ -372,23 +372,26 @@ GEnergy GCTAEdispPerfTable::mc(GRan&          ran,
  * 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_obs(const double& logEsrc,
-                                         const double& theta,
-                                         const double& phi,
-                                         const double& zenith,
-                                         const double& azimuth) const
+GEbounds GCTAEdispPerfTable::ereco_bounds(const GEnergy& etrue,
+                                          const double&  theta,
+                                          const double&  phi,
+                                          const double&  zenith,
+                                          const double&  azimuth) const
 {
     // Set energy band constant
     const double number_of_sigmas = 5.0;
 
+    // Get log10 of true photon energy in TeV
+    double etrue_log10TeV = etrue.log10TeV();
+
     // Get energy dispersion sigma
-    double sigma = m_logE.interpolate(logEsrc, m_sigma);
+    double sigma = m_logE.interpolate(etrue_log10TeV, m_sigma);
 
     // Compute energy boundaries
     GEnergy emin;
     GEnergy emax;
-    emin.log10TeV(logEsrc - number_of_sigmas * sigma);
-    emax.log10TeV(logEsrc + number_of_sigmas * sigma);
+    emin.log10TeV(etrue_log10TeV - number_of_sigmas * sigma);
+    emax.log10TeV(etrue_log10TeV + number_of_sigmas * sigma);
 
     // Return energy boundaries
     return (GEbounds(emin, emax));
@@ -398,7 +401,7 @@ GEbounds GCTAEdispPerfTable::ebounds_obs(const double& logEsrc,
 /***********************************************************************//**
  * @brief Return true energy interval that contains the energy dispersion.
  *
- * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
+ * @param[in] ereco Reconstructed event 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.
@@ -410,23 +413,26 @@ GEbounds GCTAEdispPerfTable::ebounds_obs(const double& logEsrc,
  * 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,
-                                         const double& theta,
-                                         const double& phi,
-                                         const double& zenith,
-                                         const double& azimuth) const
+GEbounds GCTAEdispPerfTable::etrue_bounds(const GEnergy& ereco,
+                                          const double&  theta,
+                                          const double&  phi,
+                                          const double&  zenith,
+                                          const double&  azimuth) const
 {
     // Set energy band constant
     const double number_of_sigmas = 5.0;
 
+    // Get log10 of reconstructed event energy in TeV
+    double ereco_log10TeV = ereco.log10TeV();
+
     // Get energy dispersion sigma
-    double sigma = m_logE.interpolate(logEobs, m_sigma);
+    double sigma = m_logE.interpolate(ereco_log10TeV, m_sigma);
 
     // Compute energy boundaries
     GEnergy emin;
     GEnergy emax;
-    emin.log10TeV(logEobs - number_of_sigmas * sigma);
-    emax.log10TeV(logEobs + number_of_sigmas * sigma);
+    emin.log10TeV(ereco_log10TeV - number_of_sigmas * sigma);
+    emax.log10TeV(ereco_log10TeV + number_of_sigmas * sigma);
 
     // Return energy boundaries
     return (GEbounds(emin, emax));
diff --git a/inst/cta/src/GCTAEdispRmf.cpp b/inst/cta/src/GCTAEdispRmf.cpp
index 0e9f7d5..7df06b2 100644
--- a/inst/cta/src/GCTAEdispRmf.cpp
+++ b/inst/cta/src/GCTAEdispRmf.cpp
@@ -307,7 +307,7 @@ GEnergy GCTAEdispRmf::mc(GRan&          ran,
     GEnergy ereco;
 
     // Get boundaries for rejection method
-    GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth);
+    GEbounds ebounds = ereco_bounds(etrue, theta, phi, zenith, azimuth);
     double   emin    = ebounds.emin().log10TeV();
     double   emax    = ebounds.emax().log10TeV();
 
@@ -333,7 +333,7 @@ GEnergy GCTAEdispRmf::mc(GRan&          ran,
 /***********************************************************************//**
  * @brief Return observed energy interval that contains the energy dispersion.
  *
- * @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.
@@ -341,13 +341,13 @@ GEnergy GCTAEdispRmf::mc(GRan&          ran,
  * @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.
+ * dispersion becomes negligible for a given true energy @p etrue.
  ***************************************************************************/
-GEbounds GCTAEdispRmf::ebounds_obs(const double& logEsrc,
-                                   const double& theta,
-                                   const double& phi,
-                                   const double& zenith,
-                                   const double& azimuth) const
+GEbounds GCTAEdispRmf::ereco_bounds(const GEnergy& etrue,
+                                    const double&  theta,
+                                    const double&  phi,
+                                    const double&  zenith,
+                                    const double&  azimuth) const
 {
     // Compute only if parameters changed
     if (!m_ebounds_obs_computed || theta != m_last_theta_obs) {
@@ -355,6 +355,7 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double& logEsrc,
         // Store last theta value
         m_ebounds_obs_computed = true;
         m_last_theta_obs       = theta;
+        m_last_etrue_ereco_bounds.log10TeV(-30.0); // force update
 
         // Compute ebounds_obs
         compute_ebounds_obs(theta, phi, zenith, azimuth);
@@ -362,18 +363,21 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double& logEsrc,
     }
 
     // Search index only if logEsrc has changed
-    if (logEsrc != m_last_logEsrc) {
+    if (etrue != m_last_etrue_ereco_bounds) {
 
         // Store last photon energy
-        m_last_logEsrc = logEsrc;
+        m_last_etrue_ereco_bounds = etrue;
+
+        // Get true photon energy in TeV
+        double etrue_TeV = etrue.TeV();
 
         // Find right index with bisection
         int low  = 0;
         int high = m_ebounds_obs.size() - 1;
         while ((high-low) > 1) {
             int  mid = (low+high) / 2;
-            double e = m_rmf.etrue().emin(mid).log10TeV();
-            if (logEsrc < e) {
+            double e = m_rmf.etrue().emin(mid).TeV();
+            if (etrue_TeV < e) {
                 high = mid;
             }
             else {
@@ -394,7 +398,7 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double& logEsrc,
 /***********************************************************************//**
  * @brief Return true energy interval that contains the energy dispersion.
  *
- * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
+ * @param[in] ereco Reconstructed event 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.
@@ -404,11 +408,11 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double& logEsrc,
  * Returns the band of true photon energies outside of which the energy
  * dispersion becomes negligible for a given observed energy @p logEobs.
  ***************************************************************************/
-GEbounds GCTAEdispRmf::ebounds_src(const double& logEobs,
-                                   const double& theta,
-                                   const double& phi,
-                                   const double& zenith,
-                                   const double& azimuth) const
+GEbounds GCTAEdispRmf::etrue_bounds(const GEnergy& ereco,
+                                    const double&  theta,
+                                    const double&  phi,
+                                    const double&  zenith,
+                                    const double&  azimuth) const
 {
     // Compute only if parameters changed
     if (!m_ebounds_src_computed || theta != m_last_theta_src) {
@@ -416,25 +420,29 @@ GEbounds GCTAEdispRmf::ebounds_src(const double& logEobs,
         // Store last value
         m_ebounds_src_computed = true;
         m_last_theta_src       = theta;
+        m_last_ereco_etrue_bounds.log10TeV(-30.0); // force update
 
         // Compute ebounds_obs
         compute_ebounds_src(theta, phi, zenith, azimuth);
 
     }
 
-    // Search index only if logEobs has changed
-    if (logEobs != m_last_logEobs) {
+    // Search index only if ereco has changed
+    if (ereco != m_last_ereco_etrue_bounds) {
 
         // Store last value
-        m_last_logEobs = logEobs;
+        m_last_ereco_etrue_bounds = ereco;
+
+        // Get reconstructed energy in TeV
+        double ereco_TeV = ereco.TeV();
 
         // Find right index with bisection
         int low  = 0;
         int high = m_ebounds_src.size() - 1;
         while ((high-low) > 1) {
             int  mid = (low+high) / 2;
-            double e = m_rmf.emeasured().emin(mid).log10TeV();
-            if (logEobs < e) {
+            double e = m_rmf.emeasured().emin(mid).TeV();
+            if (ereco_TeV < e) {
                 high = mid;
             }
             else {
@@ -627,8 +635,8 @@ void GCTAEdispRmf::init_members(void)
     m_max_edisp            = 0.0;
     m_last_theta_obs       = -1.0;
     m_last_theta_src       = -1.0;
-    m_last_logEsrc         = -30.0;
-    m_last_logEobs         = -30.0;
+    m_last_etrue_ereco_bounds.log10TeV(-30.0);
+    m_last_ereco_etrue_bounds.log10TeV(-30.0);
     m_index_obs            = 0;
     m_index_src            = 0;
     m_ebounds_obs_computed = false;
@@ -668,17 +676,17 @@ void GCTAEdispRmf::copy_members(const GCTAEdispRmf& edisp)
     m_wgt4           = edisp.m_wgt4;
 
     // Copy Monte Carlo cache
-    m_max_edisp            = edisp.m_max_edisp;
-    m_last_theta_obs       = edisp.m_last_theta_obs;
-    m_last_theta_src       = edisp.m_last_theta_src;
-    m_last_logEsrc         = edisp.m_last_logEsrc;
-    m_last_logEobs         = edisp.m_last_logEobs;
-    m_index_obs            = edisp.m_index_obs;
-    m_index_src            = edisp.m_index_src;
-    m_ebounds_obs_computed = edisp.m_ebounds_obs_computed;
-    m_ebounds_src_computed = edisp.m_ebounds_src_computed;
-    m_ebounds_obs          = edisp.m_ebounds_obs;
-    m_ebounds_src          = edisp.m_ebounds_src;
+    m_max_edisp               = edisp.m_max_edisp;
+    m_last_theta_obs          = edisp.m_last_theta_obs;
+    m_last_theta_src          = edisp.m_last_theta_src;
+    m_last_etrue_ereco_bounds = edisp.m_last_etrue_ereco_bounds;
+    m_last_ereco_etrue_bounds = edisp.m_last_ereco_etrue_bounds;
+    m_index_obs               = edisp.m_index_obs;
+    m_index_src               = edisp.m_index_src;
+    m_ebounds_obs_computed    = edisp.m_ebounds_obs_computed;
+    m_ebounds_src_computed    = edisp.m_ebounds_src_computed;
+    m_ebounds_obs             = edisp.m_ebounds_obs;
+    m_ebounds_src             = edisp.m_ebounds_src;
 
     // Return
     return;
@@ -721,11 +729,8 @@ void GCTAEdispRmf::set_matrix(void)
     // Normalize matrix elements
     for (int itrue = 0; itrue < rows; ++itrue) {
 
-        // Get true photon energy
-        double logEsrc = m_rmf.etrue().elogmean(itrue).log10TeV();
-
         // Get integration boundaries
-        GEbounds ebounds = ebounds_obs(logEsrc, 0.0);
+        GEbounds ebounds = ereco_bounds(m_rmf.etrue().elogmean(itrue), 0.0);
 
         // Integrate contributions over energy boundaries
         double sum = 0.0;
diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp
index 620d5c8..5df5a3b 100644
--- a/inst/cta/src/GCTAResponseIrf.cpp
+++ b/inst/cta/src/GCTAResponseIrf.cpp
@@ -450,8 +450,7 @@ double GCTAResponseIrf::nroi(const GModelSky&    model,
     if (use_edisp()) {
 
         // Retrieve true energy boundaries
-        double   obsLogEng = obsEng.log10TeV();
-        GEbounds ebounds   = edisp()->ebounds_src(obsLogEng);
+        GEbounds ebounds = edisp()->etrue_bounds(obsEng);
 
         // Loop over all boundaries
         for (int i = 0; i < ebounds.size(); ++i) {
@@ -550,8 +549,7 @@ GEbounds GCTAResponseIrf::ebounds(const GEnergy& obsEnergy) const
 
     // If energy dispersion is available then set the energy boundaries
     if (edisp() != NULL) {
-        double obsLogEng = obsEnergy.log10TeV();
-        ebounds          = edisp()->ebounds_src(obsLogEng); // Requires TeV
+        ebounds = edisp()->etrue_bounds(obsEnergy);
     }
 
     // Return energy boundaries
diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp
index 25f0b9c..74ec6f9 100644
--- a/inst/cta/test/test_CTA.cpp
+++ b/inst/cta/test/test_CTA.cpp
@@ -890,12 +890,12 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
                "Test filename after loading perfromance table");
 
     // Test ebounds_obs method
-    test_value(edisp2.ebounds_obs(-1.0).emin().GeV(),  31.6795127349936);
-    test_value(edisp2.ebounds_obs(-1.0).emax().GeV(), 315.661420794325);
+    test_value(edisp2.ereco_bounds(GEnergy(0.1,"TeV")).emin().GeV(),  31.6795127349936);
+    test_value(edisp2.ereco_bounds(GEnergy(0.1,"TeV")).emax().GeV(), 315.661420794325);
 
     // Test ebounds_src method
-    test_value(edisp2.ebounds_src(-1.0).emin().GeV(),  31.6795127349936);
-    test_value(edisp2.ebounds_src(-1.0).emax().GeV(), 315.661420794325);
+    test_value(edisp2.etrue_bounds(GEnergy(0.1,"TeV")).emin().GeV(),  31.6795127349936);
+    test_value(edisp2.etrue_bounds(GEnergy(0.1,"TeV")).emax().GeV(), 315.661420794325);
 
     // Test if non-diagonal element (below diagonal) is zero
     test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(10.0,"TeV")), 0.0);
@@ -947,12 +947,12 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
                "Test filename after loading RMF file");
 
     // Test ebounds_obs method
-    test_value(edisp2.ebounds_obs(1.0).emin().TeV(),  5.75439929962158);
-    test_value(edisp2.ebounds_obs(1.0).emax().TeV(), 15.8489322662354);
+    test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emin().TeV(),  5.75439929962158);
+    test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 15.8489322662354);
 
     // Test ebounds_src method
-    test_value(edisp2.ebounds_src(1.0).emin().TeV(),  5.75439929962158);
-    test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8489322662354);
+    test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emin().TeV(),  5.75439929962158);
+    test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 15.8489322662354);
 
     // Test if non-diagonal element (below diagonal) is zero
     test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(1.0,"TeV")), 0.0);
@@ -1002,12 +1002,12 @@ void TestGCTAResponse::test_response_edisp_2D(void)
                "Test filename after loading response table file");
 
     // Test ebounds_obs method
-    test_value(edisp2.ebounds_obs(1.0).emin().TeV(),  6.82816152379073);
-    test_value(edisp2.ebounds_obs(1.0).emax().TeV(), 13.2011115155942);
+    test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emin().TeV(),  6.82816152379073);
+    test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 13.2011115155942);
 
     // Test ebounds_src method
-    test_value(edisp2.ebounds_src(1.0).emin().TeV(),  7.7813808018878);
-    test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8059295273666);
+    test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emin().TeV(),  7.7813808018878);
+    test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 15.8059295273666);
 
     // Test if non-diagonal element (below diagonal) is zero
     test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(1.0,"TeV")), 0.0);
@@ -1427,10 +1427,9 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf& rs
 
 	        // Compute log10 of true energy
             GEnergy etrue(e_src,"TeV");
-	        double  log10_e_src = etrue.log10TeV();
 
-	        // Retrieve boundaries in observed energy
-	        GEbounds ebounds = rsp.edisp()->ebounds_obs(log10_e_src);
+	        // Retrieve observed energy boundaries
+	        GEbounds ebounds = rsp.edisp()->ereco_bounds(etrue);
 	        double   emin    = ebounds.emin().MeV();
 	        double   emax    = ebounds.emax().MeV();
 
@@ -1480,10 +1479,9 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&   edisp,
 
         // Compute log10 of true energy
         GEnergy etrue(e_src,"TeV");
-        double  log10_e_src = etrue.log10TeV();
 
-        // Retrieve boundaries in observed energy
-        GEbounds ebounds  = edisp.ebounds_obs(log10_e_src);
+        // Retrieve observed energy boundaries
+        GEbounds ebounds = edisp.ereco_bounds(etrue);
 
         // Skip this energy if the boundaries are empty
         if (ebounds.is_empty()) {