Dmitry Malyshev / gammalib

Commit 2f590d36e2d448ef379ae9622d67e54a656ba13b

Authored by Josh Cardenzana
1 parent fddb03f0

Removed old code from GModelSpectralSmoothBrokenPlaw (#1948) 

Removed all the deprecated code fragments from previous MC generation method.
Inline Side-by-side
Showing 3 changed files with 10 additions and 148 deletions
include/GModelSpectralSmoothBrokenPlaw.hpp
  View file @2f590d3
... ... @@ -221,11 +221,6 @@ protected:
221 221 mutable double m_mc_plaw_prefactor; //!< Prefactor for comparison p-laws
222 222 mutable double m_mc_exponentS; //!< Exponent (index+1) for softer index
223 223 mutable double m_mc_exponentH; //!< Exponent (index+1) for harder index
224   - mutable double m_mc_pow_ewidth_low; //!< width of energy range below break
225   - mutable double m_mc_norm; //!< Normalization term for energy generation
226   - mutable double m_mc_exponent_ratio; //!< Ratio of (m_mc_exponentS/m_mc_exponentH)
227   - mutable double m_mc_ebreak; //!< Break energy (depending on the energy range,
228   - //!< this may not actually be the break energy)
229 224 mutable GModelSpectralBrokenPlaw m_mc_brokenplaw;
230 225 };
231 226  
... ...
src/model/GModelSpectralSmoothBrokenPlaw.cpp
  View file @2f590d3
... ... @@ -597,86 +597,12 @@ GEnergy GModelSpectralSmoothBrokenPlaw::mc(const GEnergy&amp; emin,
597 597 double acceptance_fraction(0.0);
598 598  
599 599 // Use rejection method to draw a random energy. We first draw
600   - // analytically from a power law, and then compare the power law
601   - // at the drawn energy to the curved function. This
602   - // gives an acceptance fraction, and we accept the energy only if
603   - // a uniform random number is <= the acceptance fraction.
  600 + // analytically from a broken power law, and then compare the power law
  601 + // at the drawn energy to the curved function. This gives an acceptance
  602 + // fraction, and we accept the energy only if a uniform random number
  603 + // is <= the acceptance fraction.
604 604 do {
605   - /*
606   - double e_norm(0.0), plaw(0.0);
607   -
608   - // Get uniform random number
609   - double u = ran.uniform();
610   -
611   - // Case 1: indices are the same (things are simple)
612   - if (m_mc_exponentH == m_mc_exponentS) {
613   -
614   - // Case 1A: Both indices are -1
615   - if (m_mc_exponentH == 0.0) {
616   - eng = std::exp(u * m_mc_norm + std::log(m_mc_emin));
617   -
618   - }
619   -
620   - // Case 1B: Indices are not equal to -1
621   - else {
622   - eng = std::exp(std::log(u * m_mc_norm +
623   - std::pow(m_mc_emin,m_mc_exponentH)) / m_mc_exponentH);
624   - }
625   -
626   - // There is no differentiation between the powerlaws when the
627   - // exponents are the same.
628   - e_norm = eng / breakenergy().MeV();
629   - plaw = m_mc_plaw_prefactor * std::pow(e_norm, m_mc_exponentH-1.0);
630   -
631   - }
632   -
633   - // Case 2: indices are different (things are more complicated)
634   - else {
635   -
636   - // Case 2A: Random number suggests energy less than breakenergy
637   - if (u <= (m_mc_pow_ewidth_low/m_mc_norm)) {
638   -
639   - // Case 2A-1: Hard index is not -1
640   - if (m_mc_exponentH != 0.0) {
641   - eng = std::exp(std::log(u * m_mc_norm / m_mc_exponent_ratio +
642   - std::pow(m_mc_emin, m_mc_exponentH)) / m_mc_exponentH);
643   - }
644   -
645   - // Case 2A-2: Hard index is -1
646   - else {
647   - eng = std::exp(u * m_mc_norm + std::log(m_mc_emin));
648   - }
649   -
650   - // Compute powerlaw at given energy
651   - e_norm = eng / breakenergy().MeV();
652   - plaw = m_mc_plaw_prefactor * std::pow(e_norm, m_mc_exponentH-1.0);
653   -
654   - }
655   - // Case 2B: Random number suggests energy above breakenergy
656   - else {
657   -
658   - // Case 2B-1: Soft index is not -1
659   - if (m_mc_exponentS != 0.0) {
660   - eng = std::exp(std::log(u * m_mc_norm +
661   - std::pow(m_mc_ebreak, m_mc_exponentS) -
662   - m_mc_pow_ewidth_low) / m_mc_exponentS) ;
663   - }
664   -
665   - // Case 2B-2: Soft index is -1
666   - else {
667   - eng = std::exp(u * m_mc_norm +
668   - std::log(m_mc_ebreak) -
669   - m_mc_pow_ewidth_low);
670   - }
671   -
672   - // Compute powerlaw at given energy
673   - e_norm = eng / breakenergy().MeV();
674   - plaw = m_mc_plaw_prefactor * std::pow(e_norm, m_mc_exponentS-1.0);
675   -
676   - }
677   - }
678   - */
679   -
  605 + // Generate an energy from the broken power law distribution
680 606 energy = m_mc_brokenplaw.mc(emin, emax, time, ran);
681 607  
682 608 // Compute acceptance fraction
... ... @@ -898,9 +824,6 @@ void GModelSpectralSmoothBrokenPlaw::init_members(void)
898 824 m_mc_plaw_prefactor = 0.0;
899 825 m_mc_exponentS = 0.0;
900 826 m_mc_exponentH = 0.0;
901   - m_mc_pow_ewidth_low = 0.0;
902   - m_mc_norm = 0.0;
903   - m_mc_exponent_ratio = 1.0;
904 827 m_mc_brokenplaw = GModelSpectralBrokenPlaw();
905 828  
906 829 // Return
... ... @@ -954,10 +877,7 @@ void GModelSpectralSmoothBrokenPlaw::copy_members(
954 877 m_mc_plaw_prefactor = model.m_mc_plaw_prefactor;
955 878 m_mc_exponentS = model.m_mc_exponentS;
956 879 m_mc_exponentH = model.m_mc_exponentH;
957   - m_mc_pow_ewidth_low = model.m_mc_pow_ewidth_low;
958   - m_mc_norm = model.m_mc_norm;
959   - m_mc_exponent_ratio = model.m_mc_exponent_ratio;
960   - m_mc_brokenplaw = model.m_mc_brokenplaw;
  880 + m_mc_brokenplaw = model.m_mc_brokenplaw;
961 881  
962 882 // Return
963 883 return;
... ... @@ -1045,7 +965,7 @@ void GModelSpectralSmoothBrokenPlaw::update_mc_cache(const GEnergy&amp; emin,
1045 965  
1046 966 /* Predefine some variables */
1047 967  
1048   - // Prefactor for comparison power law functions
  968 + // Prefactor for comparison broken power law function
1049 969 m_mc_plaw_prefactor = prefactor() * std::pow(breakenergy().MeV()/pivot().MeV(),
1050 970 index1());
1051 971  
... ... @@ -1069,60 +989,7 @@ void GModelSpectralSmoothBrokenPlaw::update_mc_cache(const GEnergy&amp; emin,
1069 989 m_mc_exponentH-1.0,
1070 990 breakenergy(),
1071 991 m_mc_exponentS-1.0);
1072   - /*
1073   - // Create a break energy variable that select the appropriate range
1074   - m_mc_ebreak = breakenergy().MeV();
1075 992  
1076   - // Store the ratio of the two exponents
1077   - m_mc_exponent_ratio = (m_mc_exponentS / m_mc_exponentH) *
1078   - std::pow(m_mc_ebreak, m_mc_exponentS - m_mc_exponentH);
1079   -
1080   -
1081   - if (m_mc_emax <= m_mc_ebreak) {
1082   - m_mc_ebreak = m_mc_emax;
1083   - }
1084   - else if (m_mc_emin >= m_mc_ebreak) {
1085   - m_mc_ebreak = m_mc_emin;
1086   - }
1087   -
1088   - // Now handle the cases where the hard exponent is -1
1089   - if (m_mc_exponentH == 0.0) {
1090   -
1091   - // CASE: both exponents are -1
1092   - if (m_mc_exponentS == 0.0) {
1093   - m_mc_pow_ewidth_low = 0.0 ;
1094   - m_mc_norm = std::log(m_mc_emax) - std::log(m_mc_emin);
1095   - }
1096   -
1097   - // CASE: Only harder exponent is -1
1098   - else {
1099   - m_mc_pow_ewidth_low = std::log(m_mc_ebreak/m_mc_emin);
1100   - m_mc_norm = m_mc_pow_ewidth_low + std::pow(m_mc_emax,m_mc_exponentS) -
1101   - std::pow(m_mc_ebreak,m_mc_exponentS);
1102   - }
1103   - }
1104   -
1105   - // CASE: Only softer exponent is -1
1106   - else if (m_mc_exponentS == 0.0) {
1107   - m_mc_pow_ewidth_low = std::pow(m_mc_ebreak,m_mc_exponentH) -
1108   - std::pow(m_mc_emin, m_mc_exponentH) ;
1109   - m_mc_norm = m_mc_pow_ewidth_low + std::log(m_mc_emax) -
1110   - std::log(m_mc_ebreak) ;
1111   - }
1112   -
1113   - // CASE: Neither exponent is -1.
1114   - else {
1115   -
1116   - // Get the width of the energy range below the break energy
1117   - m_mc_pow_ewidth_low = m_mc_exponent_ratio *
1118   - (std::pow(m_mc_ebreak, m_mc_exponentH) -
1119   - std::pow(m_mc_emin, m_mc_exponentH));
1120   -
1121   - // Get the normalization term
1122   - m_mc_norm = m_mc_pow_ewidth_low + (std::pow(m_mc_emax,m_mc_exponentS) -
1123   - std::pow(m_mc_ebreak,m_mc_exponentS));
1124   - }
1125   - */
1126 993 } // endif: Update was required
1127 994  
1128 995 // Return
... ...
test/dev/test_model_spectral.py
  View file @2f590d3
... ... @@ -88,11 +88,11 @@ if __name__ == &#39;__main__&#39;:
88 88  
89 89 # Allocate spectral model
90 90 model = gammalib.GModelSpectralSmoothBrokenPlaw(5.7e-16,
91   - -2.48,
  91 + -1.5,
92 92 gammalib.GEnergy(0.3, 'TeV'),
93   - -3.5,
  93 + -4.5,
94 94 gammalib.GEnergy(1.0, 'TeV'),
95   - 0.2)
  95 + 2.0)
96 96  
97 97 # Evaluate model
98 98 eval_model(model, ntrials=1000000)
... ...