diff --git a/ChangeLog b/ChangeLog index 8955ac3..84eac9b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -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 diff --git a/NEWS b/NEWS index 70def40..953220f 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ -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 --------------------- diff --git a/inst/cta/include/GCTAEdisp2D.hpp b/inst/cta/include/GCTAEdisp2D.hpp index ae98e52..6baf505 100644 --- a/inst/cta/include/GCTAEdisp2D.hpp +++ b/inst/cta/include/GCTAEdisp2D.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * 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 diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp index 67e3167..39a1bb4 100644 --- a/inst/cta/src/GCTAEdisp2D.cpp +++ b/inst/cta/src/GCTAEdisp2D.cpp @@ -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& 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& 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& 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& 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& 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& 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& 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(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& 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& 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; diff --git a/inst/cta/src/GCTAObservation.cpp b/inst/cta/src/GCTAObservation.cpp index a426c4b..b768730 100644 --- a/inst/cta/src/GCTAObservation.cpp +++ b/inst/cta/src/GCTAObservation.cpp @@ -1441,7 +1441,7 @@ std::string GCTAObservation::print(const GChatter& chatter) const if (m_events != NULL) { result.append("\n"+m_events->print(reduced_chatter)); } - + } // endif: appended detailed information } // endif: chatter was not silent diff --git a/inst/cta/test/dev/test_background_integration.py b/inst/cta/test/dev/test_background_integration.py new file mode 100755 index 0000000..a6330c8 --- /dev/null +++ b/inst/cta/test/dev/test_background_integration.py @@ -0,0 +1,131 @@ +#! /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 . +# +# ========================================================================== +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) diff --git a/src/numerics/GIntegral.cpp b/src/numerics/GIntegral.cpp index e360663..4bb87c5 100644 --- a/src/numerics/GIntegral.cpp +++ b/src/numerics/GIntegral.cpp @@ -385,11 +385,18 @@ double GIntegral::romberg(std::vector bounds, const int& 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; } diff --git a/src/obs/GObservation.cpp b/src/obs/GObservation.cpp index 8ba4e48..d358581 100644 --- a/src/obs/GObservation.cpp +++ b/src/obs/GObservation.cpp @@ -352,17 +352,17 @@ double GObservation::model(const GModels& 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