From f2695a7538b9d5983cc293f9c4bef44cd38d1167 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 24 Oct 2018 17:20:08 +0200 Subject: [PATCH] GCTAEdisp2D::operator() now returns energy dispersion per MeV (#2704) --- inst/cta/include/GCTAEdisp2D.hpp | 44 ++++++++++++++++++++++++++++++++++++++------ inst/cta/src/GCTAEdisp2D.cpp | 226 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------- 2 files changed, 123 insertions(+), 147 deletions(-) diff --git a/inst/cta/include/GCTAEdisp2D.hpp b/inst/cta/include/GCTAEdisp2D.hpp index d224bcf..cf1e027 100644 --- a/inst/cta/include/GCTAEdisp2D.hpp +++ b/inst/cta/include/GCTAEdisp2D.hpp @@ -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 diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp index 6812e9b..e7c813a 100644 --- a/inst/cta/src/GCTAEdisp2D.cpp +++ b/inst/cta/src/GCTAEdisp2D.cpp @@ -178,7 +178,7 @@ GCTAEdisp2D& GCTAEdisp2D::operator=(const GCTAEdisp2D& 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& GCTAEdisp2D::operator=(const GCTAEdisp2D& 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& 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& 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& 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& 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& 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& 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& 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& 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; -- libgit2 0.21.4