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