Correct diffuse map cube upper limit computation in ctulimit (#2298)

Jürgen Knödlseder
1 parent 11528349
Showing 6 changed files with 96 additions and 21 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.5.rst
doc/source/users/reference_manual/ctulimit.rst
src/ctulimit/ctulimit.cpp
src/ctulimit/ctulimit.hpp
-2018-01-24
+2018-01-25
         * Version 1.5.0 released
         ========================
+        Correct diffuse map cube upper limit computation in ctulimit (#2298)
         Add tmin and tmax parameters to csobsselect (#2297)
         Change DATE_OBS, DATE_END, TIME_OBS and TIME_END FITS keywords into DATA-OBS,
           DATE-END, TIME-OBS and TIME-END
-24 January 2018
+25 January 2018
 Introduction
@@ -142,6 +142,10 @@ Add &quot;statistic&quot; parameter to allow setting of optimizer statistic (#2237).
 ctulimit - Compute upper limit
 ------------------------------
+The tool now also computes correct upper flux limits for diffuse map cubes
+for which the spectral component is a normalization of the spectral information
+available in the map cube (#2298).
+
 Add "statistic" parameter to allow setting of optimizer statistic (#2237).
@@ -20,6 +20,8 @@ In particular, this release provides:
 Bug fixes
 ---------
+* [`2298 <https://cta-redmine.irap.omp.eu/issues/2298>`_] -
+  Correct diffuse map cube upper limit computation in :ref:`ctulimit`
 * Change ``DATE_OBS``, ``DATE_END``, ``TIME_OBS`` and ``TIME_END`` FITS keywords
   into ``DATA-OBS``, ``DATE-END``, ``TIME-OBS`` and ``TIME-END``
 * [`2287 <https://cta-redmine.irap.omp.eu/issues/2287>`_] -
@@ -19,7 +19,8 @@ using the hidden ``confidence`` parameter.
 :ref:`ctulimit` writes the differential upper flux limit at a given reference
 energy (specified by the hidden parameter ``eref``) and the integrated 
 upper flux limit (specified by the hidden parameters ``emin`` and ``emax``)
-into the log file.
+into the log file. This information can also be accessed using the ``diff_ulimit()``,
+``flux_ulimit()``, and ``eflux_ulimit()`` methods.
 General parameters
 /***************************************************************************
  *                   ctulimit - Upper limit calculation tool               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2015-2017 by Michael Mayer                               *
+ *  copyright (C) 2015-2018 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -302,17 +302,8 @@ void ctulimit::run(void)
     // Write final parameter into logger
     log_value(NORMAL, "Final parameter", m_model_par->value());
-    // Get reference energy for differential upper limit
-    GEnergy eref = GEnergy(m_eref, "TeV");
-
-    // Create energy range for flux limits
-    GEnergy emin = GEnergy(m_emin, "TeV");
-    GEnergy emax = GEnergy(m_emax, "TeV");
-
     // Compute upper limit intensity and fluxes
-    m_diff_ulimit  = m_skymodel->spectral()->eval(eref, GTime());
-    m_flux_ulimit  = m_skymodel->spectral()->flux(emin, emax);
-    m_eflux_ulimit = m_skymodel->spectral()->eflux(emin, emax);
+    compute_ulimit();
     // Write header into logger
     log_header1(TERSE, "Upper limit results");
@@ -640,3 +631,78 @@ void ctulimit::ulimit_bisection(const double&amp; min, const double&amp; max)
     return;
 }
+
+
+/***********************************************************************//**
+ * @brief Compute upper limit intensity and fluxes
+ *
+ * Compute upper limit intensity and fluxes, including a correct computation
+ * for the diffuse map cube models.
+ ***************************************************************************/
+void ctulimit::compute_ulimit(void)
+{
+    // Initialise upper limit intensity and fluxes
+    m_diff_ulimit  = 0.0;
+    m_flux_ulimit  = 0.0;
+    m_eflux_ulimit = 0.0;
+
+    // Get reference energy for differential upper limit
+    GEnergy eref = GEnergy(m_eref, "TeV");
+
+    // Create energy range for flux limits
+    GEnergy emin = GEnergy(m_emin, "TeV");
+    GEnergy emax = GEnergy(m_emax, "TeV");
+
+    // Set pointer to spectral model and initialise spectral nodes model
+    // pointer
+    GModelSpectral*      spectral = m_skymodel->spectral();
+    GModelSpectralNodes* nodes    = NULL;
+
+    // If the spectral model is a diffuse cube then create a node function
+    // spectral model that is the product of the diffuse cube node function
+    // and the spectral model evaluated at the energies of the node function
+    GModelSpatialDiffuseCube* cube =
+        dynamic_cast<GModelSpatialDiffuseCube*>(m_skymodel->spatial());
+    if (cube != NULL) {
+
+
+        // Set MC cone to the entire sky. This method call is needed to
+        // set-up the cube spectrum
+        cube->set_mc_cone(GSkyDir(), 180.0);
+
+        // Allocate node function to replace the spectral component
+        nodes = new GModelSpectralNodes(cube->spectrum());
+        for (int i = 0; i < nodes->nodes(); ++i) {
+            GEnergy energy    = nodes->energy(i);
+            double  intensity = nodes->intensity(i);
+            double  value     = spectral->eval(energy);
+            nodes->intensity(i, value*intensity);
+        }
+
+        // Set the spectral model pointer to the node function. If there are
+        // no nodes the spectral pointer is set to NULL since in this case the
+        // spectral->flux method will throw an exception
+        if (nodes->nodes() > 0) {
+            spectral = nodes;
+        }
+        else {
+            spectral = NULL;
+        }
+
+    } // endif: spatial model was a diffuse cube
+
+    // Compute upper limit intensity and fluxes
+    if (spectral != NULL) {
+        m_diff_ulimit  = spectral->eval(eref);
+        m_flux_ulimit  = spectral->flux(emin, emax);
+        m_eflux_ulimit = spectral->eflux(emin, emax);
+    }
+
+    // Free spectral nodes if they have been allocated
+    if (nodes != NULL) {
+        delete nodes;
+    }
+
+    // Return
+    return;
+}
 /***************************************************************************
  *                   ctulimit - Upper limit calculation tool               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2015-2016 by Michael Mayer                               *
+ *  copyright (C) 2015-2018 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -70,12 +70,13 @@ public:
 protected:
     // Protected methods
-    void   init_members(void);
-    void   copy_members(const ctulimit& app);
-    void   free_members(void);
-    void   get_parameters(void);
-    void   get_model_parameter(void);
-    void   ulimit_bisection(const double& parmin, const double& parmax);
+    void init_members(void);
+    void copy_members(const ctulimit& app);
+    void free_members(void);
+    void get_parameters(void);
+    void get_model_parameter(void);
+    void ulimit_bisection(const double& parmin, const double& parmax);
+    void compute_ulimit(void);
     // User parameters
     std::string   m_srcname;      //!< Name of source which is moved around