From 85fdb567c623d5f5018a90f2f478a90a08be548c Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 24 Mar 2017 16:00:14 +0100 Subject: [PATCH] Integrate smoothly broken power law spectrum (#1948) --- AUTHORS | 1 + ChangeLog | 3 ++- NEWS | 6 +++++- include/GModelSpectralSmoothBrokenPlaw.hpp | 58 ++++++++++++++++++++++++++++------------------------------ src/model/GModelSpectralSmoothBrokenPlaw.cpp | 159 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------- test/data/legacy_point_smoothbplaw.xml | 17 +++++++++++++++++ test/data/model_point_smoothbplaw.xml | 42 ++++++++++++++---------------------------- test/test_GModel.cpp | 117 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------- test/test_GModel.hpp | 2 ++ 9 files changed, 231 insertions(+), 174 deletions(-) create mode 100644 test/data/legacy_point_smoothbplaw.xml diff --git a/AUTHORS b/AUTHORS index 80376bb..dc30cac 100644 --- a/AUTHORS +++ b/AUTHORS @@ -1,5 +1,6 @@ Rolf Buehler (rolf.buehler@desy.de) François Brun (brun@cenbg.in2p3.fr) +Joshua Cardenzana (jcardenzana@irap.omp.eu) Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr) Christoph Deil (Christoph.Deil@mpi-hd.mpg.de) Alex Donath (Axel.Donath@mpi-hd.mpg.de) diff --git a/ChangeLog b/ChangeLog index 84eac9b..9dcf5c0 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2017-03-22 +2017-03-24 * Version 1.3.0 released ========================= + Add GModelSpectralSmoothBrokenPlaw class (#1948) Implement lazy loading of energy dispersion in GCTAEdisp2D (#1940) Add GCTARoi::contains method Do no longer interpret a CTA calibration database as file path (#1969) diff --git a/NEWS b/NEWS index 953220f..6708d66 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.3.0.dev1 -22 March 2017 +24 March 2017 1. Introduction @@ -14,6 +14,7 @@ made since the last release of GammaLib. The following classes have been added: - GWcsGLS - GWcsSFL +- GModelSpectralSmoothBrokenPlaw The following classes have been removed: - ... @@ -69,6 +70,9 @@ None The GModelSpectralNodes class now throws an exception if the energy or intensity of a node is set to a non-positive value (#990). +The GModelSpectralSmoothBrokenPlaw class that implements a smoothly broken +power law spectrum has been added (#1948). + 9. Numerics module ------------------ diff --git a/include/GModelSpectralSmoothBrokenPlaw.hpp b/include/GModelSpectralSmoothBrokenPlaw.hpp index 2ecf031..17335e6 100644 --- a/include/GModelSpectralSmoothBrokenPlaw.hpp +++ b/include/GModelSpectralSmoothBrokenPlaw.hpp @@ -1,6 +1,6 @@ /*************************************************************************** - * GModelSpectralSmoothBrokenPlaw.cpp - * - * Smoothly broken power law spectrum class * + * GModelSpectralSmoothBrokenPlaw.hpp * + * Smoothly broken power law spectrum class * * ----------------------------------------------------------------------- * * copyright (C) 2017 by Josh Cardenzana * * ----------------------------------------------------------------------- * @@ -40,6 +40,7 @@ class GRan; class GTime; class GXmlElement; + /***********************************************************************//** * @class GModelSpectralSmoothBrokenPlaw * @@ -65,13 +66,13 @@ class GXmlElement; * Notes: * - The pivot energy should be set far away from the expected break energy * value. - * - When the two indices are close together, the beta parameter becomes poorly - * constrained. Since the beta parameter also scales the indices, this can - * cause very large errors in the estimates of the various spectral parameters. - * In this case, consider fixing beta. + * - When the two indices are close together, the beta parameter becomes + * poorly constrained. Since the beta parameter also scales the indices, this + * can cause very large errors in the estimates of the various spectral + * parameters. In this case, consider fixing beta. ***************************************************************************/ class GModelSpectralSmoothBrokenPlaw : public GModelSpectral { - + public: // Constructors and destructors GModelSpectralSmoothBrokenPlaw(void); @@ -116,14 +117,12 @@ public: virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - // Methods for getting the current values double prefactor(void) const; double index1(void) const; double index2(void) const; GEnergy pivot(void) const; GEnergy breakenergy(void) const; double beta(void) const; - // Methods for setting the parameter values void prefactor(const double& prefactor); void index1(const double& index1); void index2(const double& index2); @@ -149,15 +148,14 @@ protected: const double& index2, const GEnergy& breakenergy, const double& beta) : - m_prefactor(prefactor), - m_index1(index1), - m_index2(index2), - m_pivot(pivot), - m_breakenergy(breakenergy), - m_beta(beta) - {}; + m_prefactor(prefactor), + m_index1(index1), + m_index2(index2), + m_pivot(pivot), + m_breakenergy(breakenergy), + m_beta(beta) {}; - // Method + // Evaluation method double eval(const double& energy) { double epivot = energy / m_pivot.MeV(); double ebreak = energy / m_breakenergy.MeV(); @@ -185,10 +183,10 @@ protected: const double& index2, const GEnergy& breakenergy, const double& beta): - flux_kern(prefactor, index1, pivot, index2, breakenergy, beta) - {}; + flux_kern(prefactor, index1, pivot, index2, breakenergy, beta) + {}; - // Method + // Evaluation method double eval(const double& energy) { return energy * flux_kern::eval(energy); } @@ -204,18 +202,18 @@ protected: GModelPar m_beta; //!< Break smoothness // Cached members used for pre-computations - mutable GEnergy m_last_energy; //!< Last energy value - mutable double m_last_index1; //!< Last index1 parameter - mutable double m_last_index2; //!< Last index2 parameter - mutable double m_last_pivot; //!< Last pivot parameter - mutable double m_last_breakenergy; //!< Last breakenergy parameter - mutable double m_last_beta; //!< Last beta parameter - mutable double m_last_epivot_norm; //!< Last E/Epivot value - mutable double m_last_ebreak_norm; //!< Last E/Ebreakenergy value + mutable GEnergy m_last_energy; //!< Last energy value + mutable double m_last_index1; //!< Last index1 parameter + mutable double m_last_index2; //!< Last index2 parameter + mutable double m_last_pivot; //!< Last pivot parameter + mutable double m_last_breakenergy; //!< Last breakenergy parameter + mutable double m_last_beta; //!< Last beta parameter + mutable double m_last_epivot_norm; //!< Last E/Epivot value + mutable double m_last_ebreak_norm; //!< Last E/Ebreakenergy value mutable double m_last_log_epivot_norm; //!< Last ln(E/Epivot) value mutable double m_last_log_ebreak_norm; //!< Last ln(E/Ebreakenergy) value - mutable double m_last_epivot_pow; //!< Last pow(E/Epivot,index1) value - mutable double m_last_ebreak_pow; //!< Last pow(E/Ebreakenergy,(index1-index2)/beta) + mutable double m_last_epivot_pow; //!< Last pow(E/Epivot,index1) value + mutable double m_last_ebreak_pow; //!< Last pow(E/Ebreakenergy,(index1-index2)/beta) mutable double m_mc_emin; //!< Minimum energy mutable double m_mc_emax; //!< Maximum energy diff --git a/src/model/GModelSpectralSmoothBrokenPlaw.cpp b/src/model/GModelSpectralSmoothBrokenPlaw.cpp index 16ee40d..7c64139 100644 --- a/src/model/GModelSpectralSmoothBrokenPlaw.cpp +++ b/src/model/GModelSpectralSmoothBrokenPlaw.cpp @@ -1,6 +1,6 @@ /*************************************************************************** - * GModelSpectralSmoothBrokenPlaw.cpp - * - * Smoothly broken power law spectrum class * + * GModelSpectralSmoothBrokenPlaw.cpp * + * Smoothly broken power law spectrum class * * ----------------------------------------------------------------------- * * copyright (C) 2017 by Joshua Cardenzana * * ----------------------------------------------------------------------- * @@ -21,7 +21,7 @@ ***************************************************************************/ /** * @file GModelSpectralSmoothBrokenPlaw.cpp - * @brief Broken power law spectrum class implementation + * @brief Smoothly broken power law spectrum class implementation * @author Joshua Cardenzana */ @@ -47,7 +47,7 @@ const GModelSpectralSmoothBrokenPlaw g_spectral_sblaw_seed1("SmoothBrokenPowerLa "Index2", "BreakEnergy", "BreakSmoothness"); -const GModelSpectralRegistry g_spectral_sblaw_registry1(&g_spectral_sblaw_seed1); +const GModelSpectralRegistry g_spectral_sblaw_registry1(&g_spectral_sblaw_seed1); const GModelSpectralSmoothBrokenPlaw g_spectral_sblaw_seed2("SmoothBrokenPowerLaw", "Prefactor", "Index1", @@ -55,13 +55,13 @@ const GModelSpectralSmoothBrokenPlaw g_spectral_sblaw_seed2("SmoothBrokenPowerLa "Index2", "BreakValue", "Beta"); -const GModelSpectralRegistry g_spectral_sblaw_registry2(&g_spectral_sblaw_seed2); +const GModelSpectralRegistry g_spectral_sblaw_registry2(&g_spectral_sblaw_seed2); /* __ Method name definitions ____________________________________________ */ -#define G_MC "GModelSpectralSmoothBrokenPlaw::mc(GEnergy&, GEnergy&, GTime&,"\ -" GRan&)" -#define G_READ "GModelSpectralSmoothBrokenPlaw::read(GXmlElement&)" -#define G_WRITE "GModelSpectralSmoothBrokenPlaw::write(GXmlElement&)" +#define G_MC "GModelSpectralSmoothBrokenPlaw::mc(GEnergy&, GEnergy&, GTime&,"\ + " GRan&)" +#define G_READ "GModelSpectralSmoothBrokenPlaw::read(GXmlElement&)" +#define G_WRITE "GModelSpectralSmoothBrokenPlaw::write(GXmlElement&)" /* __ Macros _____________________________________________________________ */ @@ -80,7 +80,7 @@ const GModelSpectralRegistry g_spectral_sblaw_registry2(&g_spectral_sblaw_seed2 * @brief Void constructor ***************************************************************************/ GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw(void) : - GModelSpectral() + GModelSpectral() { // Initialise private members for clean destruction init_members(); @@ -93,13 +93,13 @@ GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw(void) : /***********************************************************************//** * @brief Model type and parameter name constructor * - * @param[in] type Model type. - * @param[in] prefactor Name of prefactor parameter. - * @param[in] index1 Name of index1 parameter. - * @param[in] pivot Name of pivot parameter. - * @param[in] index2 Name of index2 parameter. - * @param[in] breakenergy Name of breakenergy parameter. - * @param[in] beta Name of beta parameter. + * @param[in] type Model type. + * @param[in] prefactor Name of prefactor parameter. + * @param[in] index1 Name of index1 parameter. + * @param[in] pivot Name of pivot parameter. + * @param[in] index2 Name of index2 parameter. + * @param[in] breakenergy Name of breakenergy parameter. + * @param[in] beta Name of beta parameter. ***************************************************************************/ GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw( const std::string& type, @@ -109,7 +109,7 @@ GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw( const std::string& index2, const std::string& breakenergy, const std::string& beta) : -GModelSpectral() + GModelSpectral() { // Initialise members init_members(); @@ -133,14 +133,14 @@ GModelSpectral() /***********************************************************************//** * @brief Parameter constructor * - * @param[in] prefactor Smoothly broken Plaw pre factor (ph/cm2/s/MeV). - * @param[in] index1 Smoothly broken Plaw index1. - * @param[in] pivot Smoothly broken Plaw pivot energy - * @param[in] index2 Smoothly broken Plaw index1. - * @param[in] breakenergy Break energy. - * @param[in] beta Break smoothness parameter + * @param[in] prefactor Smoothly broken power law pre factor (ph/cm2/s/MeV). + * @param[in] index1 Smoothly broken power law index1. + * @param[in] pivot Smoothly broken power law pivot energy + * @param[in] index2 Smoothly broken power law index1. + * @param[in] breakenergy Break energy. + * @param[in] beta Break smoothness parameter * - * Constructs a spectral power law using the model parameters + * Constructs a smoothly broken power law using the model parameters * * power law @p prefactor (ph/cm2/s/MeV) * spectral @p index1 @@ -156,7 +156,7 @@ GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw( const double& index2, const GEnergy& breakenergy, const double& beta) : -GModelSpectral() + GModelSpectral() { // Initialise members init_members(); @@ -182,13 +182,12 @@ GModelSpectral() * * @param[in] xml XML element containing position information. * - * Constructs a power law spectral model by extracting information from an - * XML element. See the read() method for more information about the expected - * structure of the XML element. + * Constructs a smoothly broken power law spectral model by extracting + * information from an XML element. See the read() method for more + * information about the expected structure of the XML element. ***************************************************************************/ -GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw( - const GXmlElement& xml) : - GModelSpectral() +GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw(const GXmlElement& xml) : + GModelSpectral() { // Initialise members init_members(); @@ -208,7 +207,7 @@ GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw( ***************************************************************************/ GModelSpectralSmoothBrokenPlaw::GModelSpectralSmoothBrokenPlaw( const GModelSpectralSmoothBrokenPlaw& model) : - GModelSpectral(model) + GModelSpectral(model) { // Initialise members init_members(); @@ -409,6 +408,7 @@ double GModelSpectralSmoothBrokenPlaw::eval(const GEnergy& srcEng, ? value * m_index2.scale() * m_last_log_ebreak_norm * m_last_ebreak_pow / (1.0 + m_last_ebreak_pow) : 0.0; + // Compute pivot and break energy value gradients double g_pivot = (m_pivot.is_free()) ? -value * m_last_index1 / m_pivot.factor_value() @@ -417,6 +417,7 @@ double GModelSpectralSmoothBrokenPlaw::eval(const GEnergy& srcEng, ? value * (m_last_index1-m_last_index2) * m_last_ebreak_pow / ((1.0+m_last_ebreak_pow) * m_breakenergy.factor_value()) : 0.0; + // Compute beta gradient double g_beta = (m_beta.is_free()) ? value * m_beta.scale() * @@ -435,7 +436,7 @@ double GModelSpectralSmoothBrokenPlaw::eval(const GEnergy& srcEng, } // endif: gradient computation was requested // Compile option: Check for NaN/Inf -#if defined(G_NAN_CHECK) + #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { std::cout << "*** ERROR: GModelSpectralSmoothBrokenPlaw::eval"; std::cout << "(srcEng=" << srcEng; @@ -456,7 +457,7 @@ double GModelSpectralSmoothBrokenPlaw::eval(const GEnergy& srcEng, std::cout << ", m_last_ebreak_pow=" << m_last_ebreak_pow; std::cout << ")" << std::endl; } -#endif + #endif // Return return value; @@ -572,7 +573,8 @@ double GModelSpectralSmoothBrokenPlaw::eflux(const GEnergy& emin, * @exception GException::erange_invalid * Energy range is invalid (emin < emax required). * - * Returns Monte Carlo energy by randomly drawing from a smoothly broken power law. + * Returns Monte Carlo energy by randomly drawing from a smoothly broken + * power law. ***************************************************************************/ GEnergy GModelSpectralSmoothBrokenPlaw::mc(const GEnergy& emin, const GEnergy& emax, @@ -611,16 +613,19 @@ GEnergy GModelSpectralSmoothBrokenPlaw::mc(const GEnergy& emin, // Case 1: indices are the same (things are simple) if (m_mc_exponentH == m_mc_exponentS) { + // Case 1A: Both indices are -1 if (m_mc_exponentH == 0.0) { eng = std::exp(u * m_mc_norm + std::log(m_mc_emin)); } + // Case 1B: Indices are not equal to -1 else { eng = std::exp(std::log(u * m_mc_norm + std::pow(m_mc_emin,m_mc_exponentH)) / m_mc_exponentH); } + // There is no differentiation between the powerlaws when the // exponents are the same. e_norm = eng / breakenergy().MeV(); @@ -654,12 +659,14 @@ GEnergy GModelSpectralSmoothBrokenPlaw::mc(const GEnergy& emin, } // Case 2B: Random number suggests energy above breakenergy else { + // Case 2B-1: Soft index is not -1 if (m_mc_exponentS != 0.0) { eng = std::exp(std::log(u * m_mc_norm + std::pow(m_mc_ebreak, m_mc_exponentS) - m_mc_pow_ewidth_low) / m_mc_exponentS) ; } + // Case 2B-2: Soft index is -1 else { eng = std::exp(u * m_mc_norm + @@ -741,7 +748,7 @@ void GModelSpectralSmoothBrokenPlaw::write(GXmlElement& xml) const // Verify model type if (xml.attribute("type") != type()) { throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"), - "Spectral model is not of type \""+type()+"\"."); + "Spectral model is not of type \""+type()+"\"."); } // Get XML parameters @@ -768,7 +775,7 @@ void GModelSpectralSmoothBrokenPlaw::write(GXmlElement& xml) const /***********************************************************************//** * @brief Print smooth broken powerlaw information * - * @param[in] chatter Chattiness (defaults to NORMAL). + * @param[in] chatter Chattiness. * @return String containing model information. ***************************************************************************/ std::string GModelSpectralSmoothBrokenPlaw::print(const GChatter& chatter) const @@ -882,17 +889,17 @@ void GModelSpectralSmoothBrokenPlaw::init_members(void) // Initialise eval cache m_last_energy.clear(); - m_last_index1 = 1.0e30; - m_last_index2 = 1.0e30; - m_last_pivot = 1.0e30; - m_last_breakenergy = 1.0e30; - m_last_beta = 1.0e30; - m_last_epivot_norm = 1.0e30; - m_last_ebreak_norm = 1.0e30; + m_last_index1 = 1.0e30; + m_last_index2 = 1.0e30; + m_last_pivot = 1.0e30; + m_last_breakenergy = 1.0e30; + m_last_beta = 1.0e30; + m_last_epivot_norm = 1.0e30; + m_last_ebreak_norm = 1.0e30; m_last_log_epivot_norm = 1.0e30; m_last_log_ebreak_norm = 1.0e30; - m_last_epivot_pow = 1.0e30; - m_last_ebreak_pow = 1.0e30; + m_last_epivot_pow = 1.0e30; + m_last_ebreak_pow = 1.0e30; // Initialise MC cache m_mc_emin = 0.0; @@ -912,9 +919,10 @@ void GModelSpectralSmoothBrokenPlaw::init_members(void) /***********************************************************************//** * @brief Copy class members * - * @param[in] model GModelSpectralSmoothBrokenPlaw members which should be copied. + * @param[in] model Smooth broken power law. ***************************************************************************/ -void GModelSpectralSmoothBrokenPlaw::copy_members(const GModelSpectralSmoothBrokenPlaw& model) +void GModelSpectralSmoothBrokenPlaw::copy_members( + const GModelSpectralSmoothBrokenPlaw& model) { // Copy members m_type = model.m_type; @@ -935,18 +943,18 @@ void GModelSpectralSmoothBrokenPlaw::copy_members(const GModelSpectralSmoothBrok m_pars.push_back(&m_beta); // Copy eval cache - m_last_energy = model.m_last_energy; - m_last_index1 = model.m_last_index1; - m_last_index2 = model.m_last_index2; - m_last_pivot = model.m_last_pivot; - m_last_breakenergy = model.m_last_breakenergy; - m_last_beta = model.m_last_beta; - m_last_epivot_norm = model.m_last_epivot_norm; - m_last_ebreak_norm = model.m_last_ebreak_norm; + m_last_energy = model.m_last_energy; + m_last_index1 = model.m_last_index1; + m_last_index2 = model.m_last_index2; + m_last_pivot = model.m_last_pivot; + m_last_breakenergy = model.m_last_breakenergy; + m_last_beta = model.m_last_beta; + m_last_epivot_norm = model.m_last_epivot_norm; + m_last_ebreak_norm = model.m_last_ebreak_norm; m_last_log_epivot_norm = model.m_last_log_epivot_norm; m_last_log_ebreak_norm = model.m_last_log_ebreak_norm; - m_last_epivot_pow = model.m_last_epivot_pow; - m_last_ebreak_pow = model.m_last_ebreak_pow; + m_last_epivot_pow = model.m_last_epivot_pow; + m_last_ebreak_pow = model.m_last_ebreak_pow; // Copy MC cache m_mc_emin = model.m_mc_emin; @@ -991,24 +999,24 @@ void GModelSpectralSmoothBrokenPlaw::update_eval_cache(const GEnergy& energy) co // If the energy or one of the parameters index1, index2, breakenergy // energy, or beta has changed then recompute the cache - if ((m_last_energy != energy) || - (m_last_index1 != index1) || - (m_last_index2 != index2) || - (m_last_pivot != pivot_eng) || - (m_last_breakenergy != break_eng) || - (m_last_beta != beta)) { + if ((m_last_energy != energy) || + (m_last_index1 != index1) || + (m_last_index2 != index2) || + (m_last_pivot != pivot_eng) || + (m_last_breakenergy != break_eng) || + (m_last_beta != beta)) { // Store actual energy and parameter values - m_last_energy = energy; - m_last_index1 = index1; - m_last_index2 = index2; - m_last_pivot = pivot_eng; - m_last_breakenergy = break_eng; - m_last_beta = beta; + m_last_energy = energy; + m_last_index1 = index1; + m_last_index2 = index2; + m_last_pivot = pivot_eng; + m_last_breakenergy = break_eng; + m_last_beta = beta; // Compute and store value - double eng = energy.MeV(); - m_last_epivot_norm = eng / m_last_pivot; + double eng = energy.MeV(); + m_last_epivot_norm = eng / m_last_pivot; m_last_ebreak_norm = eng / m_last_breakenergy; m_last_log_epivot_norm = std::log(m_last_epivot_norm); m_last_log_ebreak_norm = std::log(m_last_ebreak_norm); @@ -1057,6 +1065,7 @@ void GModelSpectralSmoothBrokenPlaw::update_mc_cache(const GEnergy& emin, m_mc_exponentH = index1() + 1.0; // Hard index + 1 m_mc_exponentS = index2() + 1.0; // Soft index + 1 } + // CASE: index2 is harder, or index1==index2 else { m_mc_exponentH = index2() + 1.0; // Hard index + 1 @@ -1070,7 +1079,8 @@ void GModelSpectralSmoothBrokenPlaw::update_mc_cache(const GEnergy& emin, m_mc_ebreak = breakenergy().MeV(); if (m_mc_emax <= m_mc_ebreak) { m_mc_ebreak = m_mc_emax; - } else if (m_mc_emin >= m_mc_ebreak) { + } + else if (m_mc_emin >= m_mc_ebreak) { m_mc_ebreak = m_mc_emin; } @@ -1082,6 +1092,7 @@ void GModelSpectralSmoothBrokenPlaw::update_mc_cache(const GEnergy& emin, m_mc_pow_ewidth_low = 0.0 ; m_mc_norm = std::log(m_mc_emax) - std::log(m_mc_emin); } + // CASE: Only harder exponent is -1 else { m_mc_pow_ewidth_low = std::log(m_mc_ebreak/m_mc_emin); diff --git a/test/data/legacy_point_smoothbplaw.xml b/test/data/legacy_point_smoothbplaw.xml new file mode 100644 index 0000000..54328a2 --- /dev/null +++ b/test/data/legacy_point_smoothbplaw.xml @@ -0,0 +1,17 @@ + + + + + + + + + + + + + + + + + diff --git a/test/data/model_point_smoothbplaw.xml b/test/data/model_point_smoothbplaw.xml index c6e1c0a..a08a7d0 100644 --- a/test/data/model_point_smoothbplaw.xml +++ b/test/data/model_point_smoothbplaw.xml @@ -1,31 +1,17 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index 3541a9b..c03a38d 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -97,6 +97,7 @@ void TestGModel::set(void) m_xml_legacy_point_plaw = datadir + "/legacy_point_plaw.xml"; m_xml_legacy_point_plaw2 = datadir + "/legacy_point_plaw2.xml"; m_xml_legacy_point_eplaw = datadir + "/legacy_point_eplaw.xml"; + m_xml_legacy_point_smoothbplaw = datadir + "/legacy_point_smoothbplaw.xml"; m_xml_legacy_point_supeplaw = datadir + "/legacy_point_supeplaw.xml"; m_xml_legacy_point_logparabola = datadir + "/legacy_point_logparabola.xml"; @@ -204,6 +205,8 @@ void TestGModel::set(void) "Test GModelSpectralPlaw legacy model"); append(static_cast(&TestGModel::test_legacy_model_point_plaw2), "Test GModelSpectralPlaw2 legacy model"); + append(static_cast(&TestGModel::test_legacy_model_point_smoothbplaw), + "Test GModelSpectralSmoothBrokenPlaw legacy model"); append(static_cast(&TestGModel::test_legacy_model_point_eplaw), "Test GModelSpectralExpPlaw legacy model"); append(static_cast(&TestGModel::test_legacy_model_point_supeplaw), @@ -1861,63 +1864,51 @@ void TestGModel::test_smoothbplaw(void) test_value(model2.beta(), 2.0); - // Test XML constructor - 1 - GModels models(m_xml_model_point_smoothbplaw); - GModelSpectralSmoothBrokenPlaw* model3 = static_cast( - static_cast( - models[0])->spectral()); - test_value(model3->size(), 6); - test_value(model3->type(), "SmoothBrokenPowerLaw", "Check model type"); - test_value(model3->prefactor(), 5.7e-16); - test_value(model3->index1(), -2.48); - test_value(model3->pivot().TeV(), 1.0); - test_value(model3->index2(), -2.70); - test_value(model3->breakenergy().TeV(), 0.3); + // Test XML constructor + GXml xml(m_xml_model_point_smoothbplaw); + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0); + GModelSpectralSmoothBrokenPlaw model3(*element); + test_value(model3.size(), 6); + test_value(model3.type(), "SmoothBrokenPowerLaw", "Check model type"); + test_value(model3.prefactor(), 5.7e-16); + test_value(model3.index1(), -2.48); + test_value(model3.pivot().TeV(), 1.0); + test_value(model3.index2(), -2.70); + test_value(model3.breakenergy().TeV(), 0.3); - // Test XML constructor - 2 - GModelSpectralSmoothBrokenPlaw* model4 = static_cast( - static_cast( - models[1])->spectral()); - test_value(model4->size(), 6); - test_value(model4->type(), "SmoothBrokenPowerLaw", "Check model type"); - test_value(model4->prefactor(), 5.7e-16); - test_value(model4->index1(), -2.48); - test_value(model4->pivot().TeV(), 1.0); - test_value(model4->index2(), -2.70); - test_value(model4->breakenergy().TeV(), 0.3); // Test prefactor method - model3->prefactor(2.3e-16); - test_value(model3->prefactor(), 2.3e-16); + model3.prefactor(2.3e-16); + test_value(model3.prefactor(), 2.3e-16); // Test index1 method - model3->index1(-2.6); - test_value(model3->index1(), -2.6); + model3.index1(-2.6); + test_value(model3.index1(), -2.6); // Test pivot method - model3->pivot(GEnergy(0.5, "TeV")); - test_value(model3->pivot().TeV(), 0.5); + model3.pivot(GEnergy(0.5, "TeV")); + test_value(model3.pivot().TeV(), 0.5); // Test index2 method - model3->index2(-3.6); - test_value(model3->index2(), -3.6); + model3.index2(-3.6); + test_value(model3.index2(), -3.6); // Test breakenergy method - model3->breakenergy(GEnergy(0.5, "TeV")); - test_value(model3->breakenergy().TeV(), 0.5); + model3.breakenergy(GEnergy(0.5, "TeV")); + test_value(model3.breakenergy().TeV(), 0.5); // Test operator access const char* strarray[] = {"Prefactor", "Index1", "PivotEnergy", "Index2", "BreakEnergy", "BreakSmoothness"}; for (int i = 0; i < 4; ++i) { std::string keyname(strarray[i]); - (*model3)[keyname].remove_range(); // To allow setting of any value - (*model3)[keyname].value(2.1); - (*model3)[keyname].error(1.9); - (*model3)[keyname].gradient(0.8); - test_value((*model3)[keyname].value(), 2.1); - test_value((*model3)[keyname].error(), 1.9); - test_value((*model3)[keyname].gradient(), 0.8); + model3[keyname].remove_range(); // To allow setting of any value + model3[keyname].value(2.1); + model3[keyname].error(1.9); + model3[keyname].gradient(0.8); + test_value(model3[keyname].value(), 2.1); + test_value(model3[keyname].error(), 1.9); + test_value(model3[keyname].gradient(), 0.8); } // Exit test @@ -3192,6 +3183,52 @@ void TestGModel::test_legacy_model_point_plaw2(void) /***********************************************************************//** + * @brief Test smooth broken power law legacy model + ***************************************************************************/ +void TestGModel::test_legacy_model_point_smoothbplaw(void) +{ + // Load model from XML file + GModels models(m_xml_legacy_point_smoothbplaw); + + // Extract spectral component + GModelSpectralSmoothBrokenPlaw* spectral = + static_cast + (static_cast(models[0])->spectral()); + + + // Test model values + test_value(spectral->size(), 6); + test_value(spectral->type(), "SmoothBrokenPowerLaw"); + test_value(spectral->prefactor(), 5.7e-16); + test_value(spectral->index1(), -2.48); + test_value(spectral->pivot().TeV(), 1.0); + test_value(spectral->index2(), -2.70); + test_value(spectral->breakenergy().TeV(), 0.3); + + // Save file to disk and reload file (tests the proper saving and loading) + models.save("test_xml_legacy_point_smoothbplaw.xml"); + models.clear(); + models.load("test_xml_legacy_point_smoothbplaw.xml"); + + // Extract spectral component + spectral = static_cast + (static_cast(models[0])->spectral()); + + // Test model values + test_value(spectral->size(), 6); + test_value(spectral->type(), "SmoothBrokenPowerLaw"); + test_value(spectral->prefactor(), 5.7e-16); + test_value(spectral->index1(), -2.48); + test_value(spectral->pivot().TeV(), 1.0); + test_value(spectral->index2(), -2.70); + test_value(spectral->breakenergy().TeV(), 0.3); + + // Exit test + return; +} + + +/***********************************************************************//** * @brief Test exponentially cut off power law legacy model ***************************************************************************/ void TestGModel::test_legacy_model_point_eplaw(void) diff --git a/test/test_GModel.hpp b/test/test_GModel.hpp index cee3158..c1a2acd 100644 --- a/test/test_GModel.hpp +++ b/test/test_GModel.hpp @@ -92,6 +92,7 @@ public: void test_legacy_model_point_const(void); void test_legacy_model_point_plaw(void); void test_legacy_model_point_plaw2(void); + void test_legacy_model_point_smoothbplaw(void); void test_legacy_model_point_eplaw(void); void test_legacy_model_point_supeplaw(void); void test_legacy_model_point_logparabola(void); @@ -147,6 +148,7 @@ private: std::string m_xml_legacy_point_const; std::string m_xml_legacy_point_plaw; std::string m_xml_legacy_point_plaw2; + std::string m_xml_legacy_point_smoothbplaw; std::string m_xml_legacy_point_eplaw; std::string m_xml_legacy_point_supeplaw; std::string m_xml_legacy_point_logparabola; -- libgit2 0.21.4