Add GModelSpectralFunc spectral constructor (#3626)

Jürgen Knödlseder
1 parent 0c45c848
Showing 9 changed files with 78 additions and 21 deletions
ChangeLog
NEWS
README.md
doc/source/admin/release_history/2.0.rst
include/GModelSpectralFunc.hpp
pyext/GModelSpectralFunc.i
src/model/GModelSpectralFunc.cpp
src/model/GModelSpectralNodes.cpp
test/test_GModel.cpp
-2021-04-29
+2021-04-30
         * Version 2.0.0 released
         ========================
         Add ASCL shield to "About" page (#3630)
-        Support spectral line fitting (#3386)
+        Add GModelSpectralFunc spectral constructor (#3626)
         Add GModelSpectralBins class (#3606)
         Add carbon footprint estimation to GApplication (#3592)
         Add phase selection to GCOMSelection and GCOMDri (#3564)
@@ -34,6 +34,7 @@
         Add XML value constructors and access methods to GXmlElement class (#3431)
         Correct application elapsed CPU time for OpenMP (#3417) 
         Implement sparse vector response cache (#3403)
+        Support spectral line fitting (#3386)
         Implement GMatrixSparse::row() and GMatrixSparse::add_to_row() methods
         Integrate analytical gradients for energy dispersion computation (#3333)
         Add GFunctions abstract base class (#3328)
 New Features and Important Changes in GammaLib 2.0.0
-29 April 2021
+30 April 2021
 1. Introduction
@@ -56,6 +56,7 @@ The following methods have been added:
 - GModelSpatialPointSource::flux(GSkyRegion&, GEnergy&, GTime&)
 - GModelSpatialRadial::is_energy_dependent()
 - GModelSpatialRadial::is_time_dependent()
+- GModelSpectralFunc(GModelSpectral&, GEnergies&)
 - GModelSpectralFunc::nodes()
 - GModelSpectralFunc::is_empty()
 - GModelSpectralFunc::append()
@@ -166,6 +167,9 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods.
 9. Model module
 ---------------
+Added GModelSpectralFunc spectral model constructor that enables creating a file
+function from an arbitrary spectral model (#3626).
+
 Added GModelSpectralBins class that implements statistically independent spectral
 bins as the spectral model component. In addition, added GModelSpectralNodes::error()
 method so that intensity errors can be accessed directly (#3606).
 GammaLib information
 ====================
-* Version: 2.0.0.dev (29 April 2021)
+* Version: 2.0.0.dev (30 April 2021)
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
@@ -181,7 +181,8 @@ refer to the section &quot;Known problems&quot; for more information.
 Windows Installation
 ====================
-There have been no efforts so far to compile GammaLib under Windows.
+On Windows GammaLib needs to be installed into a virtual machine running
+a Linux distribution.
 Testing
@@ -202,8 +203,8 @@ or
     ============================================================================
     Testsuite summary for gammalib 2.0.0.dev
     ============================================================================
-    # TOTAL: 21
-    # PASS:  21
+    # TOTAL: 22
+    # PASS:  22
     # SKIP:  0
     # XFAIL: 0
     # FAIL:  0
 .. _2.0:
-GammaLib 2.0 release (23 April 2021)
+GammaLib 2.0 release (30 April 2021)
 ====================================
 Overview
@@ -35,8 +35,8 @@ Improvements
 * [`3630 <https://cta-redmine.irap.omp.eu/issues/3630>`_] -
   Add ASCL shield to **About** page
-* [`3386 <https://cta-redmine.irap.omp.eu/issues/3386>`_] -
-  Support spectral line fitting
+* [`3626 <https://cta-redmine.irap.omp.eu/issues/3626>`_] -
+  Add ``GModelSpectralFunc`` spectral constructor
 * [`3592 <https://cta-redmine.irap.omp.eu/issues/3592>`_] -
   Add carbon footprint estimation to ``GApplication``
 * [`3564 <https://cta-redmine.irap.omp.eu/issues/3564>`_] -
@@ -58,6 +58,8 @@ Improvements
   Add GModelSpectralFunc methods to build file function in memory
 * [`3501 <https://cta-redmine.irap.omp.eu/issues/3501>`_] -
   Implement vectorised response computation for COMPTEL
+* [`3386 <https://cta-redmine.irap.omp.eu/issues/3386>`_] -
+  Support spectral line fitting
 * [`3333 <https://cta-redmine.irap.omp.eu/issues/3333>`_] -
   Integrate analytical gradients for energy dispersion computation
 * [`3324 <https://cta-redmine.irap.omp.eu/issues/3324>`_] -
@@ -39,6 +39,7 @@
 /* __ Forward declarations _______________________________________________ */
 class GRan;
 class GTime;
+class GEnergies;
 class GXmlElement;
@@ -63,6 +64,7 @@ public:
     // Constructors and destructors
     GModelSpectralFunc(void);
     GModelSpectralFunc(const GFilename& filename, const double& norm);
+    GModelSpectralFunc(const GModelSpectral& model, const GEnergies& energies);
     explicit GModelSpectralFunc(const GXmlElement& xml);
     GModelSpectralFunc(const GModelSpectralFunc& model);
     virtual ~GModelSpectralFunc(void);
@@ -132,7 +134,7 @@ protected:
     mutable std::vector<double> m_epivot;    //!< Power-law pivot energies
     mutable std::vector<double> m_flux;      //!< Photon fluxes
     mutable std::vector<double> m_eflux;     //!< Energy fluxes
-    
+
     // Cached members for MC
     mutable GEnergy             m_mc_emin;   //!< Minimum energy
     mutable GEnergy             m_mc_emax;   //!< Maximum energy
@@ -39,6 +39,7 @@ public:
     // Constructors and destructors
     GModelSpectralFunc(void);
     GModelSpectralFunc(const GFilename& filename, const double& norm);
+    GModelSpectralFunc(const GModelSpectral& model, const GEnergies& energies);
     explicit GModelSpectralFunc(const GXmlElement& xml);
     GModelSpectralFunc(const GModelSpectralFunc& model);
     virtual ~GModelSpectralFunc(void);
@@ -33,6 +33,7 @@
 #include "GTools.hpp"
 #include "GCsv.hpp"
 #include "GRan.hpp"
+#include "GEnergies.hpp"
 #include "GModelSpectralFunc.hpp"
 #include "GModelSpectralRegistry.hpp"
@@ -135,6 +136,35 @@ GModelSpectralFunc::GModelSpectralFunc(const GXmlElement&amp; xml) :
 /***********************************************************************//**
+ * @brief Spectral model constructor
+ *
+ * @param[in] model Spectral model.
+ * @param[in] energies File function node energies.
+ *
+ * Constructs a spectral file function model from any spectral model.
+ * The file function normalisation will be set to unity.
+ ***************************************************************************/
+GModelSpectralFunc::GModelSpectralFunc(const GModelSpectral& model,
+                                       const GEnergies&      energies) :
+                    GModelSpectral()
+{
+    // Initialise members
+    init_members();
+
+    // Reserve space for file function nodes
+    reserve(energies.size());
+
+    // Append nodes for all energies
+    for (int i = 0; i < energies.size(); ++i) {
+        append(energies[i], model.eval(energies[i]));
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Copy constructor
  *
  * @param[in] model File function model.
@@ -178,7 +208,7 @@ GModelSpectralFunc::~GModelSpectralFunc(void)
  * @param[in] model File function model.
  * @return File function model.
  ***************************************************************************/
-GModelSpectralFunc& GModelSpectralFunc::operator= (const GModelSpectralFunc& model)
+GModelSpectralFunc& GModelSpectralFunc::operator=(const GModelSpectralFunc& model)
 {
     // Execute only if object is not identical
     if (this != &model) {
@@ -1016,7 +1046,7 @@ void GModelSpectralFunc::save(const GFilename&amp; filename, const bool&amp; clobber) co
         csv(i,0) = gammalib::str(m_lin_nodes[i]);
         csv(i,1) = gammalib::str(m_lin_values[i]);
     }
-    
+
     // Save CSV file
     csv.save(filename, " ", clobber);
@@ -100,6 +100,9 @@ GModelSpectralNodes::GModelSpectralNodes(const GModelSpectral&amp; model,
     // Initialise members
     init_members();
+    // Reserve space for nodes
+    reserve(energies.size());
+
     // Append nodes for all energies
     for (int i = 0; i < energies.size(); ++i) {
         append(energies[i], model.eval(energies[i]));
@@ -2721,6 +2721,19 @@ void TestGModel::test_filefct(void)
     test_value(model2.norm(), 2.0, "Check normalisation of loaded model");
     test_value(model2.filename().url(), m_filefct, "Check filename of loaded model");
+    // Test spectral model constructor
+    GModelSpectralPlaw plaw(5.7e-16, -2.48, GEnergy(0.3, "TeV"));
+    GEnergies          energies(3, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV"));
+    GModelSpectralFunc model8(plaw, energies);
+    test_value(model8.size(), 1);
+    test_value(model8.nodes(), 3);
+    test_value(model8.energy(0).TeV(), 0.1);
+    test_value(model8.intensity(0), 8.69231722336369e-15);
+    test_value(model8.energy(1).TeV(), 3.16227766016838);
+    test_value(model8.intensity(1), 1.65628690171949e-18);
+    test_value(model8.energy(2).TeV(), 100.0);
+    test_value(model8.intensity(2), 3.15598962890354e-22);
+
     // Test XML constructor
     GXml               xml(m_xml_model_point_filefct);
     GXmlElement*       element = xml.element(0)->element(0)->element("spectrum", 0);
@@ -2981,7 +2994,7 @@ void TestGModel::test_temp_const(void)
     // Test value constructor
     GModelTemporalConst model2(3.0);
     test_value(model2.norm(), 3.0);
-    
+
     // Test XML constructor
     GModelTemporalConst model3;
@@ -3021,7 +3034,7 @@ void TestGModel::test_temp_lightcurve(void)
     test_value(model2.filename().url(), m_temp_lightcurve,
                "Check light curve data file name");
     test_value(model2.norm(), 2.0);
-   
+
     // Test XML constructor
     GXml         xml(m_xml_model_point_temp_lightcurve);
     GXmlElement* element = xml.element(0)->element(0)->element("temporal", 0);
@@ -3385,7 +3398,7 @@ void TestGModel::test_models(void)
                     "Model \"2FGL J0005.7+3815\" found but not expected.");
         test_assert(!models.is_empty(), "Model container is empty.");
         test_value(models.size(), 1);
-        
+
         // Append model with same name
         test_try("Append model with same name");
         try {
@@ -3398,7 +3411,7 @@ void TestGModel::test_models(void)
         catch (std::exception &e) {
             test_try_failure(e);
         }
-        
+
         // Append model with different name
         model->name("Appended model");
         models.append(*model);
@@ -3460,7 +3473,7 @@ void TestGModel::test_models(void)
         models.load(m_xml_file);
         test_assert(!models.is_empty(), "Model container is empty.");
         test_value(models.size(), 1);
-        
+
         // Append identical container
         test_try("Append identical container");
         try {
@@ -3490,7 +3503,7 @@ void TestGModel::test_models(void)
     catch (std::exception &e) {
         test_try_failure(e);
     }
-    
+
     // Setup Crab model
     GModelSky crab;
@@ -3549,7 +3562,7 @@ void TestGModel::test_model_registry(void)
             test_assert(ptr != NULL, "Model pointer for \"" +
                                      registry.name(i)+"\" is NULL");
             if (ptr != NULL) {
-            
+
                 // Test model type
                 test_assert(ptr->type() == registry.name(i),
                             "Expected \""+registry.name(i)+"\" instead"
@@ -4214,7 +4227,7 @@ void TestGModel::test_flux(void)
     // Test point source outside ROI
     test_value(ps.flux(circle1), 0.0, 1.e-7, "Check point source not contained in region");
-    
+
     // Exit test
     return;
 }