diff --git a/ChangeLog b/ChangeLog index 56f6353..9bdb0cf 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -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 diff --git a/NEWS b/NEWS index ac69d02..921fdff 100644 --- a/NEWS +++ b/NEWS @@ -1,4 +1,4 @@ -24 January 2018 +25 January 2018 Introduction @@ -142,6 +142,10 @@ Add "statistic" 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). diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst index 6687787..4d593de 100644 --- a/doc/source/admin/release_history/1.5.rst +++ b/doc/source/admin/release_history/1.5.rst @@ -20,6 +20,8 @@ In particular, this release provides: Bug fixes --------- +* [`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 `_] - diff --git a/doc/source/users/reference_manual/ctulimit.rst b/doc/source/users/reference_manual/ctulimit.rst index 1df15c5..749f3ba 100644 --- a/doc/source/users/reference_manual/ctulimit.rst +++ b/doc/source/users/reference_manual/ctulimit.rst @@ -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 diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index ddf4cfa..ef37d81 100644 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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& min, const double& 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(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; +} diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp index f942418..dd4ec60 100644 --- a/src/ctulimit/ctulimit.hpp +++ b/src/ctulimit/ctulimit.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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