Integrate smoothly broken power law spectrum (#1948)

Jürgen Knödlseder
1 parent 5fbb6c3e
Showing 9 changed files with 231 additions and 174 deletions
AUTHORS
ChangeLog
NEWS
include/GModelSpectralSmoothBrokenPlaw.hpp
src/model/GModelSpectralSmoothBrokenPlaw.cpp
test/data/legacy_point_smoothbplaw.xml
test/data/model_point_smoothbplaw.xml
test/test_GModel.cpp
test/test_GModel.hpp
 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)
-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)
 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
 ------------------
 /***************************************************************************
- *               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
 /***************************************************************************
- *               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(&quot;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(&quot;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(&amp;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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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);
+<?xml version="1.0" standalone="no"?>
+<source_library title="source library">
+  <source name="Crab" type="PointSource">
+    <spectrum type="SmoothBrokenPowerLaw">
+      <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index1"      scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="Scale"       scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
+      <parameter name="Index2"      scale="-1"    value="2.70" min="0.01"  max="+5.0"   free="1"/>
+      <parameter name="BreakValue"  scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
+      <parameter name="Beta"        scale="1.0"   value="0.2"  min="0.01"  max="10.0"     free="1"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
+      <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
+    </spatialModel>
+  </source>
+</source_library>
 <?xml version="1.0" standalone="no"?>
 <source_library title="source library">
-    <source name="Crab1" type="PointSource">
-        <spectrum type="SmoothBrokenPowerLaw">
-            <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
-            <parameter name="Index1"      scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
-            <parameter name="PivotEnergy" scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
-            <parameter name="Index2"      scale="-1"    value="2.70" min="0.01"  max="+5.0"   free="1"/>
-            <parameter name="BreakEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
-            <parameter name="BreakSmoothness" scale="1.0" value="0.2" min="0.01" max="10.0"     free="1"/>
-        </spectrum>
-        <spatialModel type="PointSource">
-            <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
-            <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
-        </spatialModel>
-    </source>
-    <source name="Crab2" type="PointSource">
-        <spectrum type="SmoothBrokenPowerLaw">
-            <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
-            <parameter name="Index1"      scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
-            <parameter name="Scale"       scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
-            <parameter name="Index2"      scale="-1"    value="2.70" min="0.01"  max="+5.0"   free="1"/>
-            <parameter name="BreakValue"  scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
-            <parameter name="Beta"        scale="1.0"   value="0.2"  min="0.01"  max="10.0"     free="1"/>
-        </spectrum>
-        <spatialModel type="PointSource">
-            <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
-            <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
-        </spatialModel>
-    </source>
+  <source name="Crab" type="PointSource">
+    <spectrum type="SmoothBrokenPowerLaw">
+      <parameter name="Prefactor"       scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index1"          scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="PivotEnergy"     scale="1e6"   value="1.0"  min="0.01"  max="1000.0" free="0"/>
+      <parameter name="Index2"          scale="-1"    value="2.70" min="0.01"  max="+5.0"   free="1"/>
+      <parameter name="BreakEnergy"     scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="1"/>
+      <parameter name="BreakSmoothness" scale="1.0" value="0.2" min="0.01" max="10.0"     free="1"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
+      <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
+    </spatialModel>
+  </source>
 </source_library>
@@ -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<pfunction>(&TestGModel::test_legacy_model_point_plaw2),
            "Test GModelSpectralPlaw2 legacy model");
+    append(static_cast<pfunction>(&TestGModel::test_legacy_model_point_smoothbplaw),
+           "Test GModelSpectralSmoothBrokenPlaw legacy model");
     append(static_cast<pfunction>(&TestGModel::test_legacy_model_point_eplaw),
            "Test GModelSpectralExpPlaw legacy model");
     append(static_cast<pfunction>(&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<GModelSpectralSmoothBrokenPlaw*>(
-                                            static_cast<GModelSky*>(
-                                                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<GModelSpectralSmoothBrokenPlaw*>(
-                                            static_cast<GModelSky*>(
-                                                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<GModelSpectralSmoothBrokenPlaw*>
+                        (static_cast<GModelSky*>(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<GModelSpectralSmoothBrokenPlaw*>
+               (static_cast<GModelSky*>(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)
@@ -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;