From de01eb121bb85f19300bad41b0d41bfba80bab14 Mon Sep 17 00:00:00 2001 From: Josh Cardenzana Date: Wed, 22 Mar 2017 10:20:30 +0100 Subject: [PATCH] Resolved issues in MC generation for GModelSpectralSmoothBrokenPlaw (#1948) --- include/GModelSpectralSmoothBrokenPlaw.hpp | 29 ++++++++++++++++------------- src/model/GModelSpectralSmoothBrokenPlaw.cpp | 61 +++++++++++++++++++++++++++++++++++-------------------------- 2 files changed, 51 insertions(+), 39 deletions(-) diff --git a/include/GModelSpectralSmoothBrokenPlaw.hpp b/include/GModelSpectralSmoothBrokenPlaw.hpp index f135da8..20e77a4 100644 --- a/include/GModelSpectralSmoothBrokenPlaw.hpp +++ b/include/GModelSpectralSmoothBrokenPlaw.hpp @@ -40,7 +40,6 @@ class GRan; class GTime; class GXmlElement; -// TODO: Updated class documentation /***********************************************************************//** * @class GModelSpectralSmoothBrokenPlaw * @@ -50,17 +49,18 @@ class GXmlElement; * defined by * * \f[ - * S_{\rm E}(E | t) = k_0 \times \left \{ - * \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise} - * \right . + * S_{\rm E}(E | t) = k_0 \left( \frac{E}{E_0} \right)^{\gamma_1} + * \left[ 1 + \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}} + * \right]^{-\beta} * \f] * * where * \f$k_0\f$ is the normalization or prefactor, * \f$\gamma_1\f$ is the spectral index before the break, - * \f$\gamma_2\f$ is the spectral index after the break, and - * \f$E_b\f$ is the break energy. - * \f$\beta\f$ defines the break smoothness + * \f$\gamma_2\f$ is the spectral index after the break, + * \f$E_0\f$ is the pivot energy, + * \f$E_b\f$ is the break energy, and + * \f$\beta\f$ defines the break smoothness (larger = smoother transition) ***************************************************************************/ class GModelSpectralSmoothBrokenPlaw : public GModelSpectral { @@ -159,12 +159,12 @@ protected: // Members protected: - double m_prefactor; //!< Normalization - double m_index1; //!< Spectral index1 - double m_index2; //!< Spectral index2 - GEnergy m_pivot; //!< Pivot energy - GEnergy m_breakenergy; //!< Break energy - double m_beta; //!< Break smoothness parameter + double m_prefactor; //!< Normalization + double m_index1; //!< Spectral index1 + double m_index2; //!< Spectral index2 + GEnergy m_pivot; //!< Pivot energy + GEnergy m_breakenergy; //!< Break energy + double m_beta; //!< Break smoothness parameter }; // Class to determine the integral energy flux, derivation of flux_kern @@ -216,6 +216,9 @@ protected: mutable double m_mc_exponentH; //!< Exponent (index+1) for harder index mutable double m_mc_pow_ewidth_low; //!< width of energy range below break mutable double m_mc_norm; //!< Normalization term for energy generation + mutable double m_mc_exponent_ratio; //!< Ratio of (m_mc_exponentS/m_mc_exponentH) + mutable double m_mc_ebreak; //!< Break energy (depending on the energy range, + //!< this may not actually be the break energy) }; diff --git a/src/model/GModelSpectralSmoothBrokenPlaw.cpp b/src/model/GModelSpectralSmoothBrokenPlaw.cpp index ab41189..fa82ff6 100644 --- a/src/model/GModelSpectralSmoothBrokenPlaw.cpp +++ b/src/model/GModelSpectralSmoothBrokenPlaw.cpp @@ -318,12 +318,9 @@ GModelSpectralSmoothBrokenPlaw* GModelSpectralSmoothBrokenPlaw::clone(void) cons * Evaluates * * \f[ - * S_{\rm E}(E | t) = k_0 \times \left \{ - * \begin{eqnarray} - * \left( \frac{E}{E_b} \right)^{\gamma_1} & {\rm if\,\,} E < E_b \\ - * \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise} - * \end{eqnarray} - * \right . + * S_{\rm E}(E | t) = k_0 \left( \frac{E}{E_0} \right)^{\gamma_1} + * \left[ 1 + \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}} + * \right]^{-\beta} * \f] * * where @@ -394,11 +391,6 @@ double GModelSpectralSmoothBrokenPlaw::eval(const GEnergy& srcEng, ? value * m_index1.scale() * (m_last_log_epivot_norm - ((m_last_ebreak_pow * m_last_log_ebreak_norm) / (m_last_ebreak_pow + 1.0))) - /* - value * m_index1.scale() / (1.0+m_last_ebreak_pow) * - (m_last_log_epivot_norm+(m_last_log_epivot_norm-m_last_log_ebreak_norm)* - m_last_ebreak_pow) - */ : 0.0; double g_index2 = (m_index2.is_free()) ? value * m_index2.scale() * m_last_log_ebreak_norm * @@ -622,15 +614,18 @@ GEnergy GModelSpectralSmoothBrokenPlaw::mc(const GEnergy& emin, plaw = m_mc_plaw_prefactor * std::pow(e_norm, m_mc_exponentH-1.0); } + // Case 2: indices are different (things are more complicated) else { + double break_u = m_mc_pow_ewidth_low/m_mc_norm; + // Case 2A: Random number suggests energy less than breakenergy if (u <= (m_mc_pow_ewidth_low/m_mc_norm)) { // Case 2A-1: Hard index is not -1 if (m_mc_exponentH != 0.0) { - eng = std::exp(std::log(u * m_mc_norm + + eng = std::exp(std::log(u * m_mc_norm / m_mc_exponent_ratio + std::pow(m_mc_emin, m_mc_exponentH)) / m_mc_exponentH); } @@ -649,13 +644,13 @@ GEnergy GModelSpectralSmoothBrokenPlaw::mc(const GEnergy& emin, // 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(breakenergy().MeV(), m_mc_exponentS) - + 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 + - std::log(breakenergy().MeV()) - + std::log(m_mc_ebreak) - m_mc_pow_ewidth_low); } @@ -894,6 +889,7 @@ void GModelSpectralSmoothBrokenPlaw::init_members(void) m_mc_exponentH = 0.0; m_mc_pow_ewidth_low = 0.0; m_mc_norm = 0.0; + m_mc_exponent_ratio = 1.0; // Return return; @@ -947,6 +943,7 @@ void GModelSpectralSmoothBrokenPlaw::copy_members(const GModelSpectralSmoothBrok m_mc_exponentH = model.m_mc_exponentH; m_mc_pow_ewidth_low = model.m_mc_pow_ewidth_low; m_mc_norm = model.m_mc_norm; + m_mc_exponent_ratio = model.m_mc_exponent_ratio; // Return return; @@ -1044,14 +1041,24 @@ void GModelSpectralSmoothBrokenPlaw::update_mc_cache(const GEnergy& emin, // CASE: index1 is harder if (index1() > index2()) { - m_mc_exponentH = index1() + 1.0 ; // Hard index + 1 - m_mc_exponentS = index2() + 1.0 ; // Soft index + 1 - + 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 - m_mc_exponentS = index1() + 1.0 ; // Soft index + 1 + m_mc_exponentH = index2() + 1.0; // Hard index + 1 + m_mc_exponentS = index1() + 1.0; // Soft index + 1 + } + + // Store the ratio of the two exponents + m_mc_exponent_ratio = m_mc_exponentS / m_mc_exponentH; + + // Create a break energy variable that select the appropriate range + 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) { + m_mc_ebreak = m_mc_emin; } // Now handle the cases where the hard exponent is -1 @@ -1064,29 +1071,31 @@ void GModelSpectralSmoothBrokenPlaw::update_mc_cache(const GEnergy& emin, } // CASE: Only harder exponent is -1 else { - m_mc_pow_ewidth_low = std::log(breakenergy().MeV()/m_mc_emin); + m_mc_pow_ewidth_low = std::log(m_mc_ebreak/m_mc_emin); m_mc_norm = m_mc_pow_ewidth_low + std::pow(m_mc_emax,m_mc_exponentS) - - std::pow(breakenergy().MeV(),m_mc_exponentS); + std::pow(m_mc_ebreak,m_mc_exponentS); } } // CASE: Only softer exponent is -1 else if (m_mc_exponentS == 0.0) { - m_mc_pow_ewidth_low = std::pow(breakenergy().MeV(),m_mc_exponentH) - + m_mc_pow_ewidth_low = std::pow(m_mc_ebreak,m_mc_exponentH) - std::pow(m_mc_emin, m_mc_exponentH) ; m_mc_norm = m_mc_pow_ewidth_low + std::log(m_mc_emax) - - std::log(breakenergy().MeV()) ; + std::log(m_mc_ebreak) ; } // CASE: Neither exponent is -1. else { + // Get the width of the energy range below the break energy - m_mc_pow_ewidth_low = std::pow(breakenergy().MeV(), m_mc_exponentH) - - std::pow(m_mc_emin, m_mc_exponentH); + m_mc_pow_ewidth_low = m_mc_exponent_ratio * + (std::pow(m_mc_ebreak, m_mc_exponentH) - + std::pow(m_mc_emin, m_mc_exponentH)); // Get the normalization term m_mc_norm = m_mc_pow_ewidth_low + (std::pow(m_mc_emax,m_mc_exponentS) - - std::pow(breakenergy().MeV(),m_mc_exponentS)); + std::pow(m_mc_ebreak,m_mc_exponentS)); } } // endif: Update was required -- libgit2 0.21.4