Justine Devin / gammalib

Commit f332ce84d6427d0a6025f343c55c91517d567367

Authored by Michael Mayer
1 parent 5116e38a

Fixed a mc()-bug and added correct reading of fermi-style xml-files

Inline Side-by-side
Showing 3 changed files with 31 additions and 17 deletions
include/GModelSpectralLogParabola.hpp
  View file @f332ce8
... ... @@ -142,6 +142,7 @@ protected:
142 142 mutable double m_mc_exponent; //!< Exponent (index+1)
143 143 mutable double m_mc_pow_emin; //!< Power of minimum energy
144 144 mutable double m_mc_pow_ewidth; //!< Power of energy width
  145 + mutable double m_mc_norm; //!< Norm of powerlaw model at logparabola pivot energy
145 146 };
146 147  
147 148 #endif /* GMODELSPECTRALLOGPARABOLA_HPP */
... ...
pyext/GModelSpectralLogParabola.i
  View file @f332ce8
... ... @@ -65,6 +65,7 @@ public:
65 65 double index(void) const;
66 66 double curvature(void) const;
67 67 double pivot(void) const;
  68 +
68 69 };
69 70  
70 71  
... ...
src/model/GModelSpectralLogParabola.cpp
  View file @f332ce8
... ... @@ -438,7 +438,6 @@ double GModelSpectralLogParabola::eflux(const GEnergy&amp; emin,
438 438 *
439 439 * Returns Monte Carlo energy by randomly drawing from a power law.
440 440 *
441   - * @todo To be implemented.
442 441 ***************************************************************************/
443 442 GEnergy GModelSpectralLogParabola::mc(const GEnergy& emin,
444 443 const GEnergy& emax,
... ... @@ -489,10 +488,10 @@ GEnergy GModelSpectralLogParabola::mc(const GEnergy&amp; emin,
489 488  
490 489 // Compute powerlaw at given energy
491 490 double e_norm = eng / pivot();
492   - double plaw = std::pow(e_norm, m_mc_exponent-1.0);
  491 + double plaw = m_mc_norm*std::pow(eng/pivot(), m_mc_exponent-1.0);
493 492  
494 493 // Compute logparabola at given energy
495   - double logparabola = plaw * std::pow(e_norm,curvature()*std::log(e_norm));
  494 + double logparabola = norm()*std::pow(e_norm,index()+curvature()*std::log(e_norm));
496 495  
497 496 // Compute acceptance fraction
498 497 acceptance_fraction = logparabola / plaw;
... ... @@ -596,18 +595,28 @@ void GModelSpectralLogParabola::read(const GXmlElement&amp; xml)
596 595 }
597 596  
598 597 // Handle index
599   - else if (par->attribute("name") == "Index" ||
600   - par->attribute("name") == "alpha") {
  598 + else if (par->attribute("name") == "Index"){
601 599 m_index.read(*par);
602 600 npar[1]++;
603 601 }
  602 + // change sign if index is defined Fermi-like
  603 + else if(par->attribute("name") == "alpha") {
  604 + m_index.read(*par);
  605 + m_index.scale(-m_index.scale());
  606 + npar[1]++;
  607 + }
604 608  
605   - // Handle index
606   - else if (par->attribute("name") == "Curvature" ||
607   - par->attribute("name") == "beta") {
  609 + // Handle curvature
  610 + else if (par->attribute("name") == "Curvature") {
608 611 m_curvature.read(*par);
609 612 npar[2]++;
610 613 }
  614 + // change sign if curvature is defined Fermi-like
  615 + else if(par->attribute("name") == "beta") {
  616 + m_curvature.read(*par);
  617 + m_curvature.scale(-m_curvature.scale());
  618 + npar[2]++;
  619 + }
611 620  
612 621 // Handle pivot energy
613 622 else if (par->attribute("name") == "Scale" ||
... ... @@ -726,7 +735,7 @@ void GModelSpectralLogParabola::write(GXmlElement&amp; xml) const
726 735  
727 736  
728 737 /***********************************************************************//**
729   - * @brief Print powerlaw information
  738 + * @brief Print LogParabola information
730 739 ***************************************************************************/
731 740 std::string GModelSpectralLogParabola::print(void) const
732 741 {
... ... @@ -812,6 +821,7 @@ void GModelSpectralLogParabola::init_members(void)
812 821 m_mc_exponent = 0.0;
813 822 m_mc_pow_emin = 0.0;
814 823 m_mc_pow_ewidth = 0.0;
  824 + m_mc_norm = 0.0;
815 825  
816 826 // Return
817 827 return;
... ... @@ -821,7 +831,7 @@ void GModelSpectralLogParabola::init_members(void)
821 831 /***********************************************************************//**
822 832 * @brief Copy class members
823 833 *
824   - * @param[in] model GModelSpectralPlaw members which should be copied.
  834 + * @param[in] model GModelSpectralLogParabola members which should be copied.
825 835 ***************************************************************************/
826 836 void GModelSpectralLogParabola::copy_members(const GModelSpectralLogParabola& model)
827 837 {
... ... @@ -844,6 +854,7 @@ void GModelSpectralLogParabola::copy_members(const GModelSpectralLogParabola&amp; mo
844 854 m_mc_exponent = model.m_mc_exponent;
845 855 m_mc_pow_emin = model.m_mc_pow_emin;
846 856 m_mc_pow_ewidth = model.m_mc_pow_ewidth;
  857 + m_mc_norm = model.m_mc_norm;
847 858  
848 859 // Return
849 860 return;
... ... @@ -882,19 +893,20 @@ void GModelSpectralLogParabola::update_mc_cache(const GEnergy&amp; emin,
882 893 // checking the sign of spectrum curvature
883 894 if(curvature() < 0){
884 895  
885   - // Use the spectral index of the logarithmic center of the interval
886   - double e_center = std::pow(10,(std::log10(emin.MeV())+std::log10(emax.MeV()))/2);
887   - index_pl = index()+curvature()*std::log(e_center/pivot());
  896 + // Use the spectral index at the pivot energy of the LogParabola
  897 + index_pl = index();
  898 + m_mc_norm = norm();
888 899 }
889 900  
890 901 else{
891   -
892 902 // Use a power law which connects the ends of the convex, curved model
893   - double y1 = std::pow(emin.MeV()/pivot(),index()+curvature()*std::log(emin.MeV()/pivot()));
894   - double y2 = std::pow(emax.MeV()/pivot(),index()+curvature()*std::log(emax.MeV()/pivot()));
895 903  
896 904 // Plaw index defined by the slope of a straight line in the log-log-plane
897   - index_pl = (std::log(y2)-std::log(y1))/(std::log(emax.MeV())-std::log(emin.MeV()));
  905 + index_pl = std::log(eval(emin)/eval(emax))/std::log(emin.MeV()/emax.MeV());
  906 +
  907 + // Plaw norm defined such that Plaw = LogParabola at emin
  908 + m_mc_norm = eval(emin)/std::pow(emin.MeV()/pivot(),index_pl);
  909 +
898 910 }
899 911  
900 912 if(index_pl != -1.0){
... ...