Merge branch 'devel' into 1948-implement-smoothbrokenpowerlaw

Josh Cardenzana
2 parents 76bc331f e0d833e6
Showing 8 changed files with 273 additions and 22 deletions
ChangeLog
NEWS
inst/cta/include/GCTAEdisp2D.hpp
inst/cta/src/GCTAEdisp2D.cpp
inst/cta/src/GCTAObservation.cpp
inst/cta/test/dev/test_background_integration.py
src/numerics/GIntegral.cpp
src/obs/GObservation.cpp
-2017-03-18
+2017-03-22
  
         * Version 1.3.0 released
         =========================
  
+        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)
         Add GCTAResponseIrf::rspname() method
-New Features and Important Changes in GammaLib 1.3.0
+New Features and Important Changes in GammaLib 1.3.0.dev1
  
-18 March 2017
+22 March 2017
  
  
 1. Introduction
@@ -24,6 +24,7 @@ The following classes have been renamed:
 The following methods have been added:
 - GCTAResponseIrf::rspname
 - GCTARoi::contains
+- GCTAEdisp2D::fetch
  
 The following methods have been removed:
 - ...
@@ -141,6 +142,10 @@ GCTAResponseIrf::read() method (#1969).
 Add GCTARoi::contains() method that checks whether an instrument direction
 is contained within the region of interest.
  
+Lazy loading of the energy dispersion in GCTAEdisp2D has been implemented by
+adding a fetch() method to the class. Now the energy dispersion is only loaded
+when needed, for example for function evaluation or Monte Carlo sampling (#1940).
+
  
 21. COMPTEL interface
 ---------------------
 /***************************************************************************
  *            GCTAEdisp2D.hpp - CTA 2D energy dispersion class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2015-2016 by Florent Forest                              *
+ *  copyright (C) 2015-2017 by Florent Forest                              *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -90,12 +90,14 @@ public:
     std::string  print(const GChatter& chatter = NORMAL) const;
  
     // Methods
+    void                     fetch(void) const;
     const GCTAResponseTable& table(void) const;
     void                     table(const GCTAResponseTable& table);
     void                     read(const GFitsTable& table);
     void                     write(GFitsBinTable& table) const;
     void                     save(const GFilename& filename,
                                   const bool&      clobber = false) const;
+
 private:
     // Methods
     void init_members(void);
@@ -133,6 +135,7 @@ private:
     // Members
     mutable GFilename m_filename;   //!< Name of Edisp response file
     GCTAResponseTable m_edisp;      //!< Edisp response table
+    mutable bool      m_fetched;    //!< Signals that Edisp has been fetched
     int               m_inx_etrue;  //!< True energy index
     int               m_inx_migra;  //!< Migration index
     int               m_inx_theta;  //!< Theta index
@@ -45,6 +45,7 @@
 #define G_READ                               "GCTAEdisp2D::read(GFitsTable&)"
 #define G_MC    "GCTAEdisp2D::mc(GRan&, double&, double&, double&, double&, "\
                                                                    "double&)"
+#define G_FETCH                                        "GCTAEdisp2D::fetch()"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -191,6 +192,9 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
                                const double& zenith,
                                const double& azimuth) const
 {
+    // Make sure that energy dispersion is online
+    fetch();
+
     // Initalize edisp
     double edisp = 0.0;
  
@@ -407,6 +411,9 @@ void GCTAEdisp2D::read(const GFitsTable&amp; table)
  ***************************************************************************/
 void GCTAEdisp2D::write(GFitsBinTable& table) const
 {
+    // Make sure that energy dispersion is online
+    fetch();
+
     // Write background table
     m_edisp.write(table);
  
@@ -427,17 +434,8 @@ void GCTAEdisp2D::write(GFitsBinTable&amp; table) const
  ***************************************************************************/
 void GCTAEdisp2D::load(const GFilename& filename)
 {
-    // Open FITS file
-    GFits fits(filename);
-
-    // Get energy dispersion table
-    const GFitsTable& table = *fits.table(filename.extname("ENERGY DISPERSION"));
-
-    // Read energy dispersion from table
-    read(table);
-
-    // Close FITS file
-    fits.close();
+    // Clear object
+    clear();
  
     // Store filename
     m_filename = filename;
@@ -521,6 +519,9 @@ GEnergy GCTAEdisp2D::mc(GRan&amp;         ran,
                         const double& zenith,
                         const double& azimuth) const
 {
+    // Make sure that energy dispersion is online
+    fetch();
+
     // Get boundaries for observed energy
     GEbounds ebounds = ebounds_obs(logEsrc, theta, phi, zenith, azimuth);
     double   emin    = ebounds.emin().log10TeV();
@@ -582,6 +583,9 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
                                   const double& zenith,
                                   const double& azimuth) const
 {
+    // Make sure that energy dispersion is online
+    fetch();
+
     // Compute only if parameters changed
     if (!m_ebounds_obs_computed || theta != m_last_theta_obs) {
  
@@ -644,6 +648,9 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
                                   const double& zenith,
                                   const double& azimuth) const
 {
+    // Make sure that energy dispersion is online
+    fetch();
+
     // Compute only if parameters changed
     if (!m_ebounds_src_computed || theta != m_last_theta_src) {
  
@@ -687,9 +694,101 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
  
  
 /***********************************************************************//**
+ * @brief Fetch energy dispersion
+ *
+ * @exception GException::file_error
+ *            File not found.
+ *            Unable to load energy dispersion.
+ * @exception GException::invalid_value
+ *            No file name has been specified.
+ *
+ * Fetches the energy dispersion by reading it from a FITS file. This method
+ * does nothing if the energy dispersion is already loaded, if there is
+ * nothing to fetch, or if the m_filename member is empty.
+ *
+ * The method is thread save. The method checks whether the file from which
+ * the energy dispersion should be loaded actually exists.
+ ***************************************************************************/
+void GCTAEdisp2D::fetch(void) const
+{
+    // Continue only if energy dispersion has not yet been fetched
+    if (!m_fetched) {
+
+        // Continue only if the file name is not empty
+        if (!m_filename.is_empty()) {
+
+            // Throw an exception if the file does not exist
+            if (!m_filename.exists()) {
+                std::string msg = "File \""+m_filename+"\" not found. Cannot "
+                                  "fetch energy dispersion. Maybe the file has "
+                                  "been deleted in the meantime.";
+                GException::file_error(G_FETCH, msg);
+            }
+
+            // Signal that energy dispersion will be fetched (has to come
+            // before reading since the ebounds_obs() and ebounds_src() methods
+            // will otherwise call the fetch() method recursively.
+            m_fetched = true;
+
+            // Initialise exception flag
+            bool has_exception = false;
+
+            // Load energy dispersion. Catch any exception. Put the code into
+            // a critical zone as it might be called from within a parallelized
+            // thread.
+            #pragma omp critical
+            {
+            try {
+
+
+                // Open FITS file
+                GFits fits(m_filename);
+
+                // Initialise energy dispersion extension name
+                std::string extname = m_filename.extname("ENERGY DISPERSION");
+
+                // Get energy dispersion table
+                const GFitsTable& table = *fits.table(extname);
+
+                // Read energy dispersion from table
+                const_cast<GCTAEdisp2D*>(this)->read(table);
+
+                // Close FITS file
+                fits.close();
+
+            }
+            catch (...) {
+                has_exception = true;
+            }
+            }
+
+            // Throw an exception if an exception has occured
+            if (has_exception) {
+                std::string msg = "Unable to load energy dispersion from "
+                                  "file \""+m_filename+"\".";
+                throw GException::file_error(G_FETCH, msg);
+            }
+
+        } // endif: filename was not empty
+
+        // Throw an exception if the FITS file name is not known
+        else {
+            std::string msg = "Unable to fetch energy dispersion since no "
+                              "filename is specified.";
+            throw GException::invalid_value(G_FETCH, msg);
+        }
+    
+    } // endif: energy dispersion had not yet been fetched
+    
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Print energy dispersion information
  *
- * @param[in] chatter Chattiness (defaults to NORMAL).
+ * @param[in] chatter Chattiness.
  * @return String containing energy dispersion information.
  ***************************************************************************/
 std::string GCTAEdisp2D::print(const GChatter& chatter) const
@@ -700,6 +799,9 @@ std::string GCTAEdisp2D::print(const GChatter&amp; chatter) const
     // Continue only if chatter is not silent
     if (chatter != SILENT) {
  
+        // Make sure that energy dispersion is online
+        fetch();
+
         // Compute energy boundaries in TeV
         double emin = m_edisp.axis_lo(m_inx_etrue,0);
         double emax = m_edisp.axis_hi(m_inx_etrue,
@@ -754,6 +856,7 @@ void GCTAEdisp2D::init_members(void)
     // Initialise members
     m_filename.clear();
     m_edisp.clear();
+    m_fetched    = false;
     m_inx_etrue  = 0;
     m_inx_migra  = 1;
     m_inx_theta  = 2;
@@ -787,6 +890,7 @@ void GCTAEdisp2D::copy_members(const GCTAEdisp2D&amp; edisp)
     // Copy members
     m_filename   = edisp.m_filename;
     m_edisp      = edisp.m_edisp;
+    m_fetched    = edisp.m_fetched;
     m_inx_etrue  = edisp.m_inx_etrue;
     m_inx_migra  = edisp.m_inx_migra;
     m_inx_theta  = edisp.m_inx_theta;
@@ -1441,7 +1441,7 @@ std::string GCTAObservation::print(const GChatter&amp; chatter) const
             if (m_events != NULL) {
                 result.append("\n"+m_events->print(reduced_chatter));
             }
-        
+
         } // endif: appended detailed information
  
     } // endif: chatter was not silent
+#! /usr/bin/env python
+# ==========================================================================
+# This energy integration of GCTABackground3D
+#
+# Copyright (C) 2017 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import os
+import math
+import gammalib
+import matplotlib.pyplot as plt
+
+
+# ========================================= #
+# Integrate using Trapezoid rule in log-log #
+# ========================================= #
+def integrate_trapezoid_loglog(bkg, logE1, logE2):
+    """
+    Integrate using Trapezoid rule in log-log
+    """
+    # Get function values
+    logf1 = math.log(bkg(logE1, 0.0, 0.0))
+    logf2 = math.log(bkg(logE2, 0.0, 0.0))
+    
+    # Compute ln of energy boundaries
+    lnE1 = logE1 * gammalib.ln10
+    lnE2 = logE2 * gammalib.ln10
+    
+    # Compute integral using Trapezoid rule in log-log
+    integral = (lnE2 - lnE1) * (math.exp(lnE1 + logf1) +
+                                math.exp(lnE2 + logf2)) * 0.5
+
+    # Return integral
+    return integral
+
+
+# ========================================= #
+# Integrate using Trapezoid rule in lin-lin #
+# ========================================= #
+def integrate_trapezoid_linlin(bkg, logE1, logE2):
+    """
+    Integrate using Trapezoid rule in lin-lin
+    """
+    # Set number of energy bins
+    nbins = 1000
+
+    # Compute log energy bin size
+    dlogE  = (logE2-logE1)/float(nbins-1)
+
+    # Perform integration
+    integral = 0.0
+    for i in range(nbins-1):
+        logE_low  = logE1    + i*dlogE
+        logE_high = logE_low + dlogE
+        dE        = math.pow(10.0, logE_high) - math.pow(10.0, logE_low)
+        integral += dE * (bkg(logE_low, 0.0, 0.0) + bkg(logE_high, 0.0, 0.0)) * 0.5
+
+    # Return integral
+    return integral
+
+
+# =============== #
+# Plot background #
+# =============== #
+def plot_background(bkg, emin=0.020, emax=120.0):
+    """
+    Plot background
+    """
+    # Generate energy vector
+    energies = gammalib.GEnergies(100, gammalib.GEnergy(emin,'TeV'),
+                                       gammalib.GEnergy(emax,'TeV'))
+
+    # Set energy values
+    x = [energy.TeV() for energy in energies]
+    y = [bkg(energy.log10TeV(), 0.0, 0.0) for energy in energies]
+
+    # Get node energies
+    i_eng   = bkg.table().axis('ENERG')
+    e_nodes = bkg.table().axis_nodes(i_eng)
+    xo      = [math.pow(10.0,logE) for logE in e_nodes]
+    yo      = [bkg(logE, 0.0, 0.0) for logE in e_nodes]
+
+    # Compute integrals
+    for i in range(len(xo)-1):
+        if yo[i] > 0 and yo[i+1] > 0:
+            integral_loglog = integrate_trapezoid_loglog(bkg, e_nodes[i], e_nodes[i+1])
+            integral_linlin = integrate_trapezoid_linlin(bkg, e_nodes[i], e_nodes[i+1])
+            print('%f-%f: %e %e (%f)' % (xo[i], xo[i+1], integral_loglog,
+                                        integral_linlin,
+                                        integral_loglog/integral_linlin))
+
+    # Create figure
+    plt.figure(1)
+    plt.title('CTA Background')
+
+    # Plot background
+    plt.loglog(x, y, 'r-')
+    plt.loglog(xo, yo, 'ro')
+
+    # Set axes
+    plt.xlabel("Energy (TeV)")
+    plt.ylabel("Background (events/s/MeV/sr)")
+
+    # Show plot
+    plt.show()
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+    # Load background
+    bkg = gammalib.GCTABackground3D('../../caldb/data/cta/prod2/bcf/South_50h/irf_file.fits')
+    print(bkg)
+
+    # Plot background
+    plot_background(bkg)
@@ -385,11 +385,18 @@ double GIntegral::romberg(std::vector&lt;double&gt; bounds, const int&amp; order)
     // Initialise integral
     double value = 0.0;
  
+    // Initialise integration status information
+    int calls = 0;
+
     // Add integral of all intervals
     for (int i = 0; i < bounds.size()-1; ++i) {
         value += romberg(bounds[i], bounds[i+1], order);
+        calls += m_calls;
     }
  
+    // Set integration status information
+    m_calls = calls;
+
     // Return value
     return value;
 }
@@ -352,17 +352,17 @@ double GObservation::model(const GModels&amp; models,
  *
  * @returns Total number of observed events.
  *
- * Returns the total number of observed events.
+ * Returns the total number of observed events. If the observation does not
+ * contain any events the method returns zero.
  ***************************************************************************/
 int GObservation::nobserved(void) const
 {
     // Initialise number of observed events
     int nobserved = 0;
  
-    // Extract number of observed events from event contained
-    const GEvents* events = this->events();
-    if (events != NULL) {
-        nobserved = events->number();
+    // Extract number of observed events
+    if (m_events != NULL) {
+        nobserved = m_events->number();
     }
  
     // Return number of observed events