Normalise phase curves by default (#2154)

Jürgen Knödlseder
1 parent 86c0af65
Showing 7 changed files with 159 additions and 32 deletions
ChangeLog
NEWS
doc/source/users/user_manual/modules/model.rst
include/GModelTemporalPhaseCurve.hpp
pyext/GModelTemporalPhaseCurve.i
src/model/GModelTemporalPhaseCurve.cpp
src/obs/GObservations.cpp
-2017-07-06
+2017-07-18
  
         * Version 1.4.0 released
         =========================
  
+        Normalise phase curves by default (#2154)
         Add GSkyRegions::overlaps() and GSkyMap::overlaps() classes (#1980)
  
  
 New Features and Important Changes in GammaLib 1.4.0
  
-6 July 2017
+18 July 2017
  
  
 1. Introduction
@@ -23,6 +23,7 @@ The following classes have been renamed:
 The following methods have been added:
 - GSkyRegion::overlaps(GSkyRegions&)
 - GSkyMap::overlaps(GSkyRegion&)
+- GModelTemporalPhaseCurve::normalise()
  
 The following methods have been removed:
 - n.a.
@@ -64,7 +65,10 @@ None
  
 8. Model module
 ---------------
-None
+The GModelTemporalPhaseCurve models are now normalised by default unless a
+normalize="0" attribute is explicitly added into the XML file. Thus the spectral
+component of a sky model with a GModelTemporalPhaseCurve temporal component now
+gives the phase-averaged flux (#2154).
  
  
 9. Numerics module
@@ -1084,6 +1084,18 @@ between 0 and 1.
  
    *Structure of phase curve FITS file*
  
+By default, the ``NORM`` values are recomputed internally so that the phase
+averaged normalisation is one. In that case, the spectral component corresponds
+to the phase-averaged spectrum. If the internal normalisation should be disabled
+the ``normalize="0"`` attribute needs to be added to the temporal tag, i.e.
+
+.. code-block:: xml
+
+   <temporal type="PhaseCurve" file="model_temporal_phasecurve.fits" normalize="0">
+
+In that case the ``NORM`` values are directly multiplied with the spectral
+component.
+
  
 Background models
 ~~~~~~~~~~~~~~~~~
@@ -85,7 +85,8 @@ public:
                              const double&    f0,
                              const double&    f1,
                              const double&    f2,
-                             const double&    norm = 1.0);
+                             const double&    norm = 1.0,
+                             const bool&      normalize = true);
     GModelTemporalPhaseCurve(const GModelTemporalPhaseCurve& model);
     virtual ~GModelTemporalPhaseCurve(void);
  
@@ -121,6 +122,7 @@ public:
     void             f1(const double& f1);
     double           f2(void) const;
     void             f2(const double& f2);
+    bool             normalize(void) const;
  
 protected:
     // Protected methods
@@ -128,14 +130,17 @@ protected:
     void copy_members(const GModelTemporalPhaseCurve& model);
     void free_members(void);
     void load_nodes(const GFilename& filename);
+    void normalise_nodes(void);
  
     // Protected parameter members
-    GModelPar           m_norm;     //!< Normalization factor
-    GModelPar           m_mjd;      //!< Reference MJD
-    GModelPar           m_phase;    //!< Phase at reference MJD
-    GModelPar           m_f0;       //!< Frequency at reference MJD
-    GModelPar           m_f1;       //!< First freq. derivative at reference MJD
-    GModelPar           m_f2;       //!< Second freq. derivative at reference MJD
+    GModelPar           m_norm;          //!< Normalization factor
+    GModelPar           m_mjd;           //!< Reference MJD
+    GModelPar           m_phase;         //!< Phase at reference MJD
+    GModelPar           m_f0;            //!< Frequency at reference MJD
+    GModelPar           m_f1;            //!< First freq. derivative at reference MJD
+    GModelPar           m_f2;            //!< Second freq. derivative at reference MJD
+    bool                m_normalize;     //!< Normalize phase curve (default: true)
+    bool                m_has_normalize; //!< XML has normalize attribute
  
     // Protected members
     mutable GNodeArray  m_nodes;    //!< Phase values of nodes
@@ -166,7 +171,7 @@ std::string GModelTemporalPhaseCurve::classname(void) const
 inline
 std::string GModelTemporalPhaseCurve::type(void) const
 {
-    return "PhaseCurve";
+    return ("PhaseCurve");
 }
  
  
@@ -374,4 +379,18 @@ void GModelTemporalPhaseCurve::f2(const double&amp; f2)
     return;
 }
  
+
+/***********************************************************************//**
+ * @brief Return normalization flag
+ *
+ * @return True if the phase curve has been normalized, false otherwise.
+ *
+ * Signals whether a phase curve has been normalized or not.
+ ***************************************************************************/
+inline
+bool GModelTemporalPhaseCurve::normalize(void) const
+{
+    return (m_normalize);
+}
+
 #endif /* GMODELTEMPORALPHASECURVE_HPP */
@@ -46,7 +46,8 @@ public:
                              const double&    f0,
                              const double&    f1,
                              const double&    f2,
-                             const double&    norm = 1.0);
+                             const double&    norm = 1.0,
+                             const bool&      normalize = true);
     GModelTemporalPhaseCurve(const GModelTemporalPhaseCurve& model);
     virtual ~GModelTemporalPhaseCurve(void);
  
@@ -78,6 +79,7 @@ public:
     void             f1(const double& f1);
     double           f2(void) const;
     void             f2(const double& f2);
+    bool             normalize(void) const;
 };
  
  
@@ -47,6 +47,7 @@ const GModelTemporalRegistry   g_temporal_phase_registry(&amp;g_temporal_phase_seed)
 #define G_READ                 "GModelTemporalPhaseCurve::read(GXmlElement&)"
 #define G_WRITE               "GModelTemporalPhaseCurve::write(GXmlElement&)"
 #define G_LOAD_NODES       "GModelTemporalPhaseCurve::load_nodes(GFilename&)"
+#define G_NORMALISE_NODES       "GModelTemporalPhaseCurve::normalise_nodes()"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -84,6 +85,7 @@ GModelTemporalPhaseCurve::GModelTemporalPhaseCurve(void) : GModelTemporal()
  * @param[in] f1 First frequency derivative at reference time.
  * @param[in] f2 Second frequency derivative at reference time.
  * @param[in] norm Normalization factor.
+ * @param[in] normalize Normalize phase curve?
  *
  * Constructs phase curve from a list of nodes that is found in the specified
  * FITS file, a reference time and phase information at the reference time.
@@ -96,7 +98,8 @@ GModelTemporalPhaseCurve::GModelTemporalPhaseCurve(const GFilename&amp; filename,
                                                    const double&    f0,
                                                    const double&    f1,
                                                    const double&    f2,
-                                                   const double&    norm) :
+                                                   const double&    norm,
+                                                   const bool&      normalize) :
                           GModelTemporal()
 {
     // Initialise members
@@ -110,6 +113,9 @@ GModelTemporalPhaseCurve::GModelTemporalPhaseCurve(const GFilename&amp; filename,
     this->f2(f2);
     this->norm(norm);
  
+    // Set normalization flag
+    m_normalize = normalize;
+
     // Load nodes
     load_nodes(filename);
  
@@ -385,6 +391,17 @@ void GModelTemporalPhaseCurve::read(const GXmlElement&amp; xml)
     m_f1.read(*f1);
     m_f2.read(*f2);
  
+    // Get optional normalization attribute
+    m_normalize     = true;
+    m_has_normalize = false;
+    if (xml.has_attribute("normalize")) {
+        m_has_normalize = true;
+        std::string arg = xml.attribute("normalize");
+        if (arg == "0" || gammalib::tolower(arg) == "false") {
+            m_normalize = false;
+        }
+    }
+
     // Load nodes from file
     load_nodes(gammalib::xml_file_expand(xml, xml.attribute("file")));
  
@@ -446,6 +463,16 @@ void GModelTemporalPhaseCurve::write(GXmlElement&amp; xml) const
     // Set file attribute
     xml.attribute("file", gammalib::xml_file_reduce(xml, m_filename));
  
+    // Set optional normalization attribute
+    if (m_has_normalize || !m_normalize) {
+        if (m_normalize) {
+            xml.attribute("normalize", "1");
+        }
+        else {
+            xml.attribute("normalize", "0");
+        }
+    }
+
     // Return
     return;
 }
@@ -514,19 +541,24 @@ std::string GModelTemporalPhaseCurve::print(const GChatter&amp; chatter) const
         // Append information
         result.append("\n"+gammalib::parformat("Function file"));
         result.append(m_filename.url());
+        if (normalize()) {
+            result.append(" [normalized]");
+        }
         result.append("\n"+gammalib::parformat("Number of parameters"));
         result.append(gammalib::str(size()));
         for (int i = 0; i < size(); ++i) {
             result.append("\n"+m_pars[i]->print(chatter));
         }
-
-        // Append node information
         result.append("\n"+gammalib::parformat("Number of phase nodes"));
         result.append(gammalib::str(m_nodes.size()));
-        for (int i = 0; i < m_nodes.size(); ++i) {
-            result.append("\n");
-            result.append(gammalib::parformat(" Phase "+gammalib::str(m_nodes[i])));
-            result.append(gammalib::str(m_values[i]));
+
+        // EXPLICIT: Append node values
+        if (chatter >= EXPLICIT) {
+            for (int i = 0; i < m_nodes.size(); ++i) {
+                result.append("\n");
+                result.append(gammalib::parformat(" Phase "+gammalib::str(m_nodes[i])));
+                result.append(gammalib::str(m_values[i]));
+            }
         }
  
     } // endif: chatter was not silent
@@ -626,7 +658,9 @@ void GModelTemporalPhaseCurve::init_members(void)
     m_nodes.clear();
     m_values.clear();
     m_filename.clear();
-    
+    m_normalize     = true;
+    m_has_normalize = false;
+
     // Return
     return;
 }
@@ -640,15 +674,17 @@ void GModelTemporalPhaseCurve::init_members(void)
 void GModelTemporalPhaseCurve::copy_members(const GModelTemporalPhaseCurve& model)
 {
     // Copy members
-    m_norm     = model.m_norm;
-    m_mjd      = model.m_mjd;
-    m_phase    = model.m_phase;
-    m_f0       = model.m_f0;
-    m_f1       = model.m_f1;
-    m_f2       = model.m_f2;
-    m_nodes    = model.m_nodes;
-    m_values   = model.m_values;
-    m_filename = model.m_filename;
+    m_norm          = model.m_norm;
+    m_mjd           = model.m_mjd;
+    m_phase         = model.m_phase;
+    m_f0            = model.m_f0;
+    m_f1            = model.m_f1;
+    m_f2            = model.m_f2;
+    m_nodes         = model.m_nodes;
+    m_values        = model.m_values;
+    m_filename      = model.m_filename;
+    m_normalize     = model.m_normalize;
+    m_has_normalize = model.m_has_normalize;
  
     // Set parameter pointer(s)
     m_pars.clear();
@@ -797,6 +833,59 @@ void GModelTemporalPhaseCurve::load_nodes(const GFilename&amp; filename)
     // Close FITS file
     fits.close();
  
+    // Optionally normalise nodes
+    if (m_normalize) {
+        normalise_nodes();
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Normalise nodes
+ *
+ * @exception GException::invalid_value
+ *            Phase curve FITS file is invalid
+ *
+ * Normalise the node values so that the integral over the phase interval
+ * [0,1] gives an average normalisation of 1.
+ ***************************************************************************/
+void GModelTemporalPhaseCurve::normalise_nodes(void)
+{
+    // Initialise node integral
+    double sum = 0.0;
+
+    // Loop over all nodes-1
+    for (int i = 0; i < m_nodes.size()-1; ++i) {
+
+        // Compute phase values
+        double a = (m_nodes[i]   < 0.0) ? 0.0 : m_nodes[i];
+        double b = (m_nodes[i+1] > 1.0) ? 1.0 : m_nodes[i+1];
+
+        // Compute normalisation at phase values
+        double fa = m_nodes.interpolate(a, m_values);
+        double fb = m_nodes.interpolate(b, m_values);
+
+        // Compute intergal for this interval using the Trapezoid rule
+        sum += 0.5 * (b - a) * (fa + fb);
+
+    }
+
+    // Throw an exception if the integral is not positive
+    if (sum <= 0.0) {
+        std::string msg = "Integral over phase curve is not positive ("+
+                          gammalib::str(sum)+"). Please provide a valid phase "
+                          "curve file.";
+        throw GException::invalid_value(G_NORMALISE_NODES, msg);
+    }
+
+    // Normalise the node values
+    for (int i = 0; i < m_values.size(); ++i) {
+        m_values[i] /= sum;
+    }
+
     // Return
     return;
 }
@@ -880,7 +880,7 @@ int GObservations::nobserved(void) const
 /***********************************************************************//**
  * @brief Print observation list information
  *
- * @param[in] chatter Chattiness (defaults to NORMAL).
+ * @param[in] chatter Chattiness.
  * @return String containing observation list information
  ***************************************************************************/
 std::string GObservations::print(const GChatter& chatter) const
@@ -910,12 +910,12 @@ std::string GObservations::print(const GChatter&amp; chatter) const
                 result.append("\n");
                 result.append((*this)[i]->print());
             }
-        } // endif: chatter was explicit
+        }
  
         // VERBOSE: Append models
         if (chatter == VERBOSE) {
             result.append("\n"+m_models.print());
-        } // endif: chatter was verbose
+        }
  
     } // endif: chatter was not silent