GCTAEdisp2D::operator() now returns energy dispersion per MeV (#2704)

The previous code returned the energy dispersion as function of the base 10 logarithm of the energy. It is not really clear why this was the case, as the GCTAResponseIrf::edisp() class returns energy dispersion per MeV, and all high-level code assumes that energy dispersion is per MeV. Therefore, to avoid confusion and renormalization of energy dispersion information, also the GCTAEdisp2D::operator() now also returns the information per MeV. The data is stored per migration bin, and hence a simple integration over migration should be unity for a give for a given offset angle and energy. The class was adapted so that this is now the case.

GCTAEdisp2D::operator() now returns energy dispersion per MeV (#2704)
The previous code returned the energy dispersion as function of the base 10 logarithm of the energy. It is not really clear why this was the case, as the GCTAResponseIrf::edisp() class returns energy dispersion per MeV, and all high-level code assumes that energy dispersion is per MeV. Therefore, to avoid confusion and renormalization of energy dispersion information, also the GCTAEdisp2D::operator() now also returns the information per MeV. The data is stored per migration bin, and hence a simple integration over migration should be unity for a give for a given offset angle and energy. The class was adapted so that this is now the case.
Jürgen Knödlseder
1 parent 2de42d5d
Showing 2 changed files with 123 additions and 147 deletions
inst/cta/include/GCTAEdisp2D.hpp
inst/cta/src/GCTAEdisp2D.cpp
@@ -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 {
@@ -135,7 +168,6 @@ private:
     void   set_table(void);
     void   set_boundaries(void);
     void   set_max_edisp(void) const;
-    void   convert_dm2dlogE(void);
     void   normalize_table(void);
     int    table_index(const int& ietrue,
                        const int& imigra,
@@ -188,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
@@ -178,7 +178,7 @@ 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).
@@ -186,19 +186,20 @@ GCTAEdisp2D&amp; GCTAEdisp2D::operator=(const GCTAEdisp2D&amp; edisp)
  * @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.
+ * @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, 
@@ -237,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
@@ -325,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.
  ***************************************************************************/
@@ -405,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
@@ -901,7 +915,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;
@@ -923,7 +937,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
@@ -1329,8 +1343,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
@@ -1367,9 +1391,6 @@ void GCTAEdisp2D::set_table(void)
     // Set table boundaries
     set_boundaries();
-    // Convert table from dP/dm to dP/dlogE
-    convert_dm2dlogE();
-
     // Smooth energy dispersion table (kludge only for CTA)
     #if defined(G_SMOOTH_EDISP_KLUDGE)
     if (gammalib::strip_whitespace(m_edisp.telescope()) == "CTA") {
@@ -1451,21 +1472,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
@@ -1502,7 +1519,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);
@@ -1516,9 +1533,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);
@@ -1561,106 +1578,6 @@ void GCTAEdisp2D::normalize_table(void)
 /***********************************************************************//**
- * @brief Convert table from dP/dm to dP/dlog10E
- *
- * The energy dispersion table when read-in is normalized so that
- *
- * \f[
- *    \int_{{\rm migra}_{\rm min}}^{{\rm migra}_{\rm max}}
- *    E_{\rm disp}({\rm migra}) \, d{\rm migra} = 1
- * \f]
- *
- * while the class needs values that are normalize according to
- *
- * \f[
- *    \int_{\log_{10} E_{\rm reco}^{\rm min}}^{\log_{10} E_{\rm reco}^{\rm max}}
- *    E_{\rm disp}(\log_{10} E_{\rm reco}) \, d(\log_{10} E_{\rm reco}) = 1
- * \f]
- *
- * The method thus multiplies each matrix element with
- *
- * \f[
- *    \frac{{\rm migra}_{\rm max} - {\rm migra}_{\rm min}}
- *         {\log_{10} E_{\rm reco}^{\rm max} - \log_{10} E_{\rm reco}^{\rm min}}
- * \f]
- *
- * where
- * \f${\rm migra}_{\rm min}\f$ and \f${\rm migra}_{\rm max}\f$ are the
- * minimum and maximum migration values of a given bin and
- * \f$E_{\rm reco}^{\rm min}\f$ and \f$E_{\rm reco}^{\rm max}\f$ are the
- * minimum and maximum reconstructed energy values for the same bin.
- ***************************************************************************/
-void GCTAEdisp2D::convert_dm2dlogE(void)
-{
-    // Optional dump header
-    #if defined(G_DEBUG_HESS_RENORMALIZATION)
-    std::cout << "GCTAEdisp2D::hess_renormalization:" << std::endl;
-    #endif
-
-    // Get axes dimensions
-    int etrue_size = m_edisp.axis_bins(m_inx_etrue);
-    int migra_size = m_edisp.axis_bins(m_inx_migra);
-    int theta_size = m_edisp.axis_bins(m_inx_theta);
-
-    // Perform renormalization
-    for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue) {
-
-        // Get true energy in TeV
-        double etrue = std::sqrt(m_edisp.axis_lo(m_inx_etrue, i_etrue) *
-                                 m_edisp.axis_hi(m_inx_etrue, i_etrue));
-
-        // Get offset angle
-        for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
-
-            // Intialize control sum
-            #if defined(G_DEBUG_HESS_RENORMALIZATION)
-            double sum = 0.0;
-            #endif
-
-            // Loop over all migration bins
-            int offset = table_index(i_etrue, 0, i_theta);
-            int stride = table_stride(m_inx_migra);
-            for (int i_migra = 0, i = offset; i_migra < migra_size; ++i_migra, i += stride) {
-            
-                // Get migration bin size
-                double migra_bin = m_edisp.axis_hi(m_inx_migra, i_migra) -
-                                   m_edisp.axis_lo(m_inx_migra, i_migra);
-
-                // Get dlog10(E_reco) bin size
-                double ereco_min = m_edisp.axis_lo(m_inx_migra, i_migra) * etrue;
-                double ereco_max = m_edisp.axis_hi(m_inx_migra, i_migra) * etrue;
-                double dlogE_bin = std::log10(ereco_max) - std::log10(ereco_min);
-
-                // Compute normalization factor
-                double norm = migra_bin / dlogE_bin;
-
-                // Multiply matrix element
-                m_edisp(m_inx_matrix, i) *= norm;
-
-                // Update test sum
-                #if defined(G_DEBUG_HESS_RENORMALIZATION)
-                sum += m_edisp(m_inx_matrix, i) * dlogE_bin;
-                #endif
-
-            } // endfor: looped over migration bins
-
-            // Optional test dump
-            #if defined(G_DEBUG_HESS_RENORMALIZATION)
-            std::cout << "etrue=" << etrue << " TeV ";
-            std::cout << "i_theta=" << i_theta << ": ";
-            std::cout << sum << std::endl;
-            #endif
-
-        } // endfor: looped over offset angle
-
-    } // endfor: looped over true energies
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
  * @brief Return index of response table element
  *
  * @param[in] ietrue True energy index.
@@ -1811,6 +1728,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)
 {
@@ -2178,23 +2098,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] x Function value.
+ * @param[in] logEobs Base 10 logarithm of observation energy (\f$\log_{10}\f$ MeV).
+ * @return Energy dispersion (\f$(\log_{10}\f$ MeV\f$)^{-1}\f$).
  *
- * This method implements the integration kernel needed for the edisp_kern()
- * class.
+ * 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]
+ *
+ * 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
+ *
+ * \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;