diff --git a/ChangeLog b/ChangeLog index d47060d..5da4e18 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,10 @@ -2018-02-06 +2018-02-09 * Version 1.6.0 released ======================== + Add GCTABackground::rate_ebin() method (#2308) + 2018-01-23 diff --git a/NEWS b/NEWS index cbdf5a7..dbcc7c8 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.6.0 -8 February 2018 +9 February 2018 1. Introduction @@ -25,6 +25,9 @@ The following methods have been added: - gammalib::plaw_integral(double&, double&, double&, double&) - GCTAInstDir::theta() - GCTAInstDir::phi() +- GCTABackground::rate_ebin() +- GCTABackgroundPerfTable::rate_ebin() +- GCTABackground3D::rate_ebin() The following methods have been removed: - none @@ -133,6 +136,11 @@ None A GCTAInstDir::theta() and GCTAInstDir::phi() methods were added that return the distance to the camera centre in radians and the azimuth angle, respectively (#2308). +A rate_ebin() method was added to the interface of the GCTABackground class and +implemented for the GCTABackgroundPerfTable and GCTABackground3D classes. The method +computes for a given instrument direction the background rate integrated over an +energy bin (#2308). + 21. COMPTEL interface --------------------- diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 8004414..2269bf2 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -22,7 +22,8 @@ Bug fixes Improvements ------------ -* ... +* [`2308 `_] - + Add `GCTABackground::rate_ebin()` method New features diff --git a/inst/cta/include/GCTABackground.hpp b/inst/cta/include/GCTABackground.hpp index 59f7cb7..4bb9023 100644 --- a/inst/cta/include/GCTABackground.hpp +++ b/inst/cta/include/GCTABackground.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackground.hpp - CTA background model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -30,11 +30,14 @@ /* __ Includes ___________________________________________________________ */ #include #include "GBase.hpp" -#include "GFilename.hpp" -#include "GModelSpectralNodes.hpp" -#include "GCTAInstDir.hpp" /* __ Forward declarations _______________________________________________ */ +class GRan; +class GEnergy; +class GTime; +class GFilename; +class GCTAInstDir; +class GModelSpectralNodes; /***********************************************************************//** @@ -68,6 +71,9 @@ public: const GTime& time, GRan& ran) const = 0; virtual const GModelSpectralNodes& spectrum(void) const = 0; + virtual double rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const = 0; virtual std::string print(const GChatter& chatter = NORMAL) const = 0; protected: diff --git a/inst/cta/include/GCTABackground3D.hpp b/inst/cta/include/GCTABackground3D.hpp index 0feb04a..dcdba7e 100644 --- a/inst/cta/include/GCTABackground3D.hpp +++ b/inst/cta/include/GCTABackground3D.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackground3D.hpp - CTA 3D background class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2017 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -29,7 +29,9 @@ /* __ Includes ___________________________________________________________ */ #include +#include "GEnergies.hpp" #include "GFilename.hpp" +#include "GModelSpectralNodes.hpp" #include "GCTABackground.hpp" #include "GCTAResponseTable.hpp" @@ -75,6 +77,9 @@ public: const GTime& time, GRan& ran) const; const GModelSpectralNodes& spectrum(void) const; + double rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const; std::string print(const GChatter& chatter = NORMAL) const; // Methods @@ -91,17 +96,19 @@ private: void init_members(void); void copy_members(const GCTABackground3D& bgd); void free_members(void); - void set_limits(void); + void set_members(void); int index(const int& idetx, const int& idety, const int& iebin) const; void init_mc_cache(void) const; void init_mc_max_rate(void) const; double solid_angle(const double& detx1, const double& dety1, const double& detx2, const double& dety2, const double& detx3, const double& dety3) const; + double rate(const int& iebin, const double& detx, const double& dety) const; // Members GFilename m_filename; //!< Name of background response file GCTAResponseTable m_background; //!< Background response table + GEnergies m_energy; //!< Vector of energies int m_inx_detx; //!< DETX index int m_inx_dety; //!< DETY index int m_inx_energy; //!< Energy index @@ -199,6 +206,7 @@ inline void GCTABackground3D::table(const GCTAResponseTable& table) { m_background = table; + set_members(); } #endif /* GCTABACKGROUND3D_HPP */ diff --git a/inst/cta/include/GCTABackgroundPerfTable.hpp b/inst/cta/include/GCTABackgroundPerfTable.hpp index 1cef8e6..349df26 100644 --- a/inst/cta/include/GCTABackgroundPerfTable.hpp +++ b/inst/cta/include/GCTABackgroundPerfTable.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackgroundPerfTable.hpp - CTA performance table background class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -30,8 +30,10 @@ /* __ Includes ___________________________________________________________ */ #include #include "GFilename.hpp" +#include "GEnergies.hpp" +#include "GNodeArray.hpp" +#include "GModelSpectralNodes.hpp" #include "GCTABackground.hpp" -#include "GCTAResponseTable.hpp" /* __ Forward declarations _______________________________________________ */ @@ -68,6 +70,9 @@ public: const GTime& time, GRan& ran) const; const GModelSpectralNodes& spectrum(void) const; + double rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const; std::string print(const GChatter& chatter = NORMAL) const; // Methods @@ -82,6 +87,7 @@ private: void free_members(void); double solidangle(void) const; void init_mc_cache(void) const; + double rate(const double& logE) const; // Radial integration class (used by solidangle() method). Note that // the sigma parameter is given in rad^2 @@ -99,10 +105,12 @@ private: }; // Members - GFilename m_filename; //!< Name of background response file - GNodeArray m_logE; //!< log(E) nodes for background interpolation - std::vector m_background; //!< Background rate - double m_sigma; //!< Sigma for offset angle computation (0=none) + GFilename m_filename; //!< Name of background response file + GEnergies m_energy; //!< Vector of energies + GNodeArray m_log10_energy; //!< log10(E) nodes for background interpolation + std::vector m_background; //!< Background rate + std::vector m_log_background; //!< log(background rate) + double m_sigma; //!< Sigma for offset angle computation (0=none) // Monte Carlo cache mutable GModelSpectralNodes m_mc_spectrum; //!< Background spectrum @@ -142,7 +150,7 @@ GFilename GCTABackgroundPerfTable::filename(void) const inline int GCTABackgroundPerfTable::size(void) const { - return (m_logE.size()); + return (m_energy.size()); } diff --git a/inst/cta/pyext/GCTABackground.i b/inst/cta/pyext/GCTABackground.i index a0ff302..a4bf3f9 100644 --- a/inst/cta/pyext/GCTABackground.i +++ b/inst/cta/pyext/GCTABackground.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackground.i - CTA background model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -57,6 +57,9 @@ public: const GTime& time, GRan& ran) const = 0; virtual const GModelSpectralNodes& spectrum(void) const = 0; + virtual double rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const = 0; }; diff --git a/inst/cta/pyext/GCTABackground3D.i b/inst/cta/pyext/GCTABackground3D.i index 7b131cc..e2c6d47 100644 --- a/inst/cta/pyext/GCTABackground3D.i +++ b/inst/cta/pyext/GCTABackground3D.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackground3D.i - CTA 3D background class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -58,6 +58,9 @@ public: const GTime& time, GRan& ran) const; const GModelSpectralNodes& spectrum(void) const; + double rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const; // Methods bool is_valid(void) const; diff --git a/inst/cta/pyext/GCTABackgroundPerfTable.i b/inst/cta/pyext/GCTABackgroundPerfTable.i index dd708f3..39fdec1 100644 --- a/inst/cta/pyext/GCTABackgroundPerfTable.i +++ b/inst/cta/pyext/GCTABackgroundPerfTable.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackgroundPerfTable.i - CTA performance table background class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -58,6 +58,9 @@ public: const GTime& time, GRan& ran) const; const GModelSpectralNodes& spectrum(void) const; + double rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const; // Methods int size(void) const; diff --git a/inst/cta/src/GCTABackground3D.cpp b/inst/cta/src/GCTABackground3D.cpp index d38887f..f2e294e 100644 --- a/inst/cta/src/GCTABackground3D.cpp +++ b/inst/cta/src/GCTABackground3D.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackground3D.cpp - CTA 3D background class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2017 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -35,11 +35,13 @@ #include "GRan.hpp" #include "GFits.hpp" #include "GFitsBinTable.hpp" +#include "GCTAInstDir.hpp" #include "GCTABackground3D.hpp" /* __ Method name definitions ____________________________________________ */ -#define G_READ "GCTABackground3D::read(GFitsTable&)" #define G_MC "GCTABackground3D::mc(GEnergy&, GTime&, GRan&)" +#define G_SET_MEMBERS "GCTABackground3D::set_members()" +#define G_INIT_MC_CACHE "GCTABackground3D::init_mc_cache()" /* __ Macros _____________________________________________________________ */ @@ -195,68 +197,30 @@ double GCTABackground3D::operator()(const double& logE, if (is_valid() && (detx >= m_detx_min) && (detx <= m_detx_max) && (dety >= m_dety_min) && (dety <= m_dety_max)) { - // Retrieve references to node arrays - const GNodeArray& detx_nodes = m_background.axis_nodes(m_inx_detx); - const GNodeArray& dety_nodes = m_background.axis_nodes(m_inx_dety); + // Retrieve references to energy node array const GNodeArray& energy_nodes = m_background.axis_nodes(m_inx_energy); - // Set values for node arrays - detx_nodes.set_value(detx); - dety_nodes.set_value(dety); + // Set energy for node array interpolation energy_nodes.set_value(logE); - // Set weighting factors for bi-linear interpolation in DETX-DETY plane - double wgt_ll = detx_nodes.wgt_left() * dety_nodes.wgt_left(); - double wgt_lr = detx_nodes.wgt_left() * dety_nodes.wgt_right(); - double wgt_rl = detx_nodes.wgt_right() * dety_nodes.wgt_left(); - double wgt_rr = detx_nodes.wgt_right() * dety_nodes.wgt_right(); - - // Set weighting factors for energy interpolation - double wgt_emin = energy_nodes.wgt_left(); - double wgt_emax = energy_nodes.wgt_right(); - - // Get bin indices - int inx_ll_emin = index(detx_nodes.inx_left(), dety_nodes.inx_left(), - energy_nodes.inx_left()); - int inx_lr_emin = index(detx_nodes.inx_left(), dety_nodes.inx_right(), - energy_nodes.inx_left()); - int inx_rl_emin = index(detx_nodes.inx_right(), dety_nodes.inx_left(), - energy_nodes.inx_left()); - int inx_rr_emin = index(detx_nodes.inx_right(), dety_nodes.inx_right(), - energy_nodes.inx_left()); - int inx_ll_emax = index(detx_nodes.inx_left(), dety_nodes.inx_left(), - energy_nodes.inx_right()); - int inx_lr_emax = index(detx_nodes.inx_left(), dety_nodes.inx_right(), - energy_nodes.inx_right()); - int inx_rl_emax = index(detx_nodes.inx_right(), dety_nodes.inx_left(), - energy_nodes.inx_right()); - int inx_rr_emax = index(detx_nodes.inx_right(), dety_nodes.inx_right(), - energy_nodes.inx_right()); - // Bi-linearly interpolate the rates in both energy layers - double rate_emin = wgt_ll * m_background(m_inx_bgd, inx_ll_emin) + - wgt_lr * m_background(m_inx_bgd, inx_lr_emin) + - wgt_rl * m_background(m_inx_bgd, inx_rl_emin) + - wgt_rr * m_background(m_inx_bgd, inx_rr_emin); - double rate_emax = wgt_ll * m_background(m_inx_bgd, inx_ll_emax) + - wgt_lr * m_background(m_inx_bgd, inx_lr_emax) + - wgt_rl * m_background(m_inx_bgd, inx_rl_emax) + - wgt_rr * m_background(m_inx_bgd, inx_rr_emax); + double rate_emin = this->rate(energy_nodes.inx_left(), detx, dety); + double rate_emax = this->rate(energy_nodes.inx_right(), detx, dety); // If both rates are positive then perform a logarithmic interpolation // in energy if (rate_emin > 0.0 && rate_emax > 0.0) { - rate = std::exp(wgt_emin * std::log(rate_emin) + - wgt_emax * std::log(rate_emax)); - } - else { - rate = 0.0; - } - // Make sure that background rate is not negative - if (rate < 0.0) { - rate = 0.0; - } + // Perform energy interpolation + rate = std::exp(energy_nodes.wgt_left() * std::log(rate_emin) + + energy_nodes.wgt_right() * std::log(rate_emax)); + + // Make sure that background rate is not negative + if (rate < 0.0) { + rate = 0.0; + } + + } // endif: both rates were positive } // endif: background model was valid @@ -309,9 +273,6 @@ GCTABackground3D* GCTABackground3D::clone(void) const * * @param[in] table FITS table. * - * @exception GException::invalid_value - * Response table is not three-dimensional. - * * Reads the background form the FITS @p table. The following column names * are mandatory: * @@ -334,41 +295,8 @@ void GCTABackground3D::read(const GFitsTable& table) // Read background table m_background.read(table); - // Throw an exception if the table is not three-dimensional - if (m_background.axes() != 3) { - std::string msg = "Expected three-dimensional background " - "response table but found "+ - gammalib::str(m_background.axes())+ - " dimensions. Please specify a three-dimensional " - "background."; - throw GException::invalid_value(G_READ, msg); - } - - // Get mandatory indices (throw exception if not found) - m_inx_detx = m_background.axis("DETX"); - m_inx_dety = m_background.axis("DETY"); - m_inx_energy = m_background.axis("ENERG"); - m_inx_bgd = m_background.table("BGD"); - - // Get axes dimensions - m_num_detx = m_background.axis_bins(m_inx_detx); - m_num_dety = m_background.axis_bins(m_inx_dety); - m_num_energy = m_background.axis_bins(m_inx_energy); - - // Set dimension array - m_num[m_inx_detx] = m_num_detx; - m_num[m_inx_dety] = m_num_dety; - m_num[m_inx_energy] = m_num_energy; - - // Set DETX and DETY axis to radians - m_background.axis_radians(m_inx_detx); - m_background.axis_radians(m_inx_dety); - - // Set energy axis to logarithmic scale - m_background.axis_log10(m_inx_energy); - - // Set DETX and DETY limits - set_limits(); + // Set class members + set_members(); // Return return; @@ -580,9 +508,102 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy& energy, /***********************************************************************//** + * @brief Returns background count rate integrated over energy interval + * + * @param[in] dir Instrument direction. + * @param[in] emin Minimum energy of energy interval. + * @param[in] emax Maximum energy of energy interval. + * @return Integrated background count rate (counts/src/sr). + * + * Returns the background count rate for a given instrument direction that + * is integrated over a specified energy interval. + * + * If the energy interval is not positive, a zero background rate is + * returned. + ***************************************************************************/ +double GCTABackground3D::rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise rate + double rate = 0.0; + + // Continue only if the background template is valid, the DETX and DETY + // values are comprised in template range and the energy interval is + // positive + if (is_valid() && (dir.detx() >= m_detx_min) && + (dir.detx() <= m_detx_max) && + (dir.dety() >= m_dety_min) && + (dir.dety() <= m_dety_max) && + (emax > emin)) { + + // Retrieve references to node arrays + //const GNodeArray& detxs = m_background.axis_nodes(m_inx_detx); + //const GNodeArray& detys = m_background.axis_nodes(m_inx_dety); + + // Initialise first and second node + double x1 = emin.MeV(); + double x2 = 0.0; + double f1 = (*this)(emin.log10TeV(), dir.detx(), dir.dety()); + double f2 = 0.0; + + // Loop over all nodes + for (int i = 0; i < m_energy.size(); ++i) { + + // If node energy is below x1 then skip node + if (m_energy[i].MeV() <= x1) { + continue; + } + + // If node energy is above emax then use emax as energy + if (m_energy[i] > emax) { + x2 = emax.MeV(); + f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety()); + } + + // ... otherwise use node energy + else { + x2 = m_energy[i].MeV(); + f2 = this->rate(i, dir.detx(), dir.dety()); + } + + // Compute integral + if (f1 > 0.0 && f2 > 0.0) { + rate += gammalib::plaw_integral(x1, f1, x2, f2); + } + + // Set second node as first node + x1 = x2; + f1 = f2; + + // If node energy is above emax then break now + if (m_energy[i] > emax) { + break; + } + + } // endfor: looped over all nodes + + // If last node energy is below emax then compute last part of + // integral up to emax + if (x1 < emax.MeV()) { + x2 = emax.MeV(); + f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety()); + if (f1 > 0.0 && f2 > 0.0) { + rate += gammalib::plaw_integral(x1, f1, x2, f2); + } + } + + } // endif: energy interval was positive + + // Return background rate + return rate; +} + + +/***********************************************************************//** * @brief Print background information * - * @param[in] chatter Chattiness (defaults to NORMAL). + * @param[in] chatter Chattiness. * @return String containing background information. ***************************************************************************/ std::string GCTABackground3D::print(const GChatter& chatter) const @@ -661,6 +682,7 @@ void GCTABackground3D::init_members(void) // Initialise members m_filename.clear(); m_background.clear(); + m_energy.clear(); m_inx_detx = 0; m_inx_dety = 1; m_inx_energy = 2; @@ -698,6 +720,7 @@ void GCTABackground3D::copy_members(const GCTABackground3D& bgd) // Copy members m_filename = bgd.m_filename; m_background = bgd.m_background; + m_energy = bgd.m_energy; m_inx_detx = bgd.m_inx_detx; m_inx_dety = bgd.m_inx_dety; m_inx_energy = bgd.m_inx_energy; @@ -726,12 +749,59 @@ void GCTABackground3D::copy_members(const GCTABackground3D& bgd) /***********************************************************************//** - * @brief Set DETX, DETY and log10(E/TeV) limits + * @brief Delete class members + ***************************************************************************/ +void GCTABackground3D::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Set members from background table * - * Set the DETX, DETY and log10(E/TeV) limits. + * @exception GException::invalid_value + * Response table is not three-dimensional. + * + * Set class members based on the background table. The DETX and DETY axes + * will be set to radians, the energy axis will be set to log10. ***************************************************************************/ -void GCTABackground3D::set_limits(void) +void GCTABackground3D::set_members(void) { + // Throw an exception if the table is not three-dimensional + if (m_background.axes() != 3) { + std::string msg = "Expected three-dimensional background " + "response table but found "+ + gammalib::str(m_background.axes())+ + " dimensions. Please specify a three-dimensional " + "background."; + throw GException::invalid_value(G_SET_MEMBERS, msg); + } + + // Get mandatory indices (throw exception if not found) + m_inx_detx = m_background.axis("DETX"); + m_inx_dety = m_background.axis("DETY"); + m_inx_energy = m_background.axis("ENERG"); + m_inx_bgd = m_background.table("BGD"); + + // Get axes dimensions + m_num_detx = m_background.axis_bins(m_inx_detx); + m_num_dety = m_background.axis_bins(m_inx_dety); + m_num_energy = m_background.axis_bins(m_inx_energy); + + // Set dimension array + m_num[m_inx_detx] = m_num_detx; + m_num[m_inx_dety] = m_num_dety; + m_num[m_inx_energy] = m_num_energy; + + // Set DETX and DETY axis to radians + m_background.axis_radians(m_inx_detx); + m_background.axis_radians(m_inx_dety); + + // Set energy axis to logarithmic scale + m_background.axis_log10(m_inx_energy); + // Compute DETX boundaries in radians m_detx_min = m_background.axis_lo(m_inx_detx, 0) * gammalib::deg2rad; @@ -752,6 +822,17 @@ void GCTABackground3D::set_limits(void) m_logE_min = emin.log10TeV(); m_logE_max = emax.log10TeV(); + // Setup energy vector + m_energy.clear(); + for (int i = 0; i < m_num_energy; ++i) { + GEnergy elo(m_background.axis_lo(m_inx_energy, i), + m_background.axis_lo_unit(m_inx_energy)); + GEnergy ehi(m_background.axis_hi(m_inx_energy, i), + m_background.axis_hi_unit(m_inx_energy)); + GEnergy energy(std::sqrt(elo.TeV()*ehi.TeV()), "TeV"); + m_energy.append(energy); + } + // Return return; } @@ -786,16 +867,6 @@ int GCTABackground3D::index(const int& idetx, /***********************************************************************//** - * @brief Delete class members - ***************************************************************************/ -void GCTABackground3D::free_members(void) -{ - // Return - return; -} - - -/***********************************************************************//** * @brief Initialise Monte Carlo cache * * @exception GException::invalid_value @@ -923,7 +994,7 @@ void GCTABackground3D::init_mc_cache(void) const if (m_mc_spectrum.nodes() == 0) { std::string msg = "Background response table is empty. Please provide " "a valid three-dimensional background template."; - throw GException::invalid_value(G_READ, msg); + throw GException::invalid_value(G_INIT_MC_CACHE, msg); } // Return @@ -1045,3 +1116,47 @@ double GCTABackground3D::solid_angle(const double& detx1, const double& dety1, // Return solid angle return solidangle; } + + +/***********************************************************************//** + * @brief Return background rate for a given energy bin and DETX-DETY value + * (events/s/MeV/sr) + * + * @param[in] iebin Energy index. + * @param[in] detx Tangential X coordinate in nominal system (radians). + * @param[in] dety Tangential Y coordinate in nominal system (radians). + * @return Background rate (events/s/MeV/sr) + ***************************************************************************/ +double GCTABackground3D::rate(const int& iebin, const double& detx, + const double& dety) const +{ + // Retrieve references to node arrays + const GNodeArray& detxs = m_background.axis_nodes(m_inx_detx); + const GNodeArray& detys = m_background.axis_nodes(m_inx_dety); + + // Set DETX and DETY values for node array interpolation + detxs.set_value(detx); + detys.set_value(dety); + + // Get bin indices + int inx_ll = index(detxs.inx_left(), detys.inx_left(), iebin); + int inx_lr = index(detxs.inx_left(), detys.inx_right(), iebin); + int inx_rl = index(detxs.inx_right(), detys.inx_left(), iebin); + int inx_rr = index(detxs.inx_right(), detys.inx_right(), iebin); + + // Bi-linearly interpolate the rates + double rate = detxs.wgt_left() * + (m_background(m_inx_bgd, inx_ll) * detys.wgt_left() + + m_background(m_inx_bgd, inx_lr) * detys.wgt_right()) + + detxs.wgt_right() * + (m_background(m_inx_bgd, inx_rl) * detys.wgt_left() + + m_background(m_inx_bgd, inx_rr) * detys.wgt_right()); + + // Make sure that background rate is not negative + if (rate < 0.0) { + rate = 0.0; + } + + // Return + return rate; +} diff --git a/inst/cta/src/GCTABackgroundPerfTable.cpp b/inst/cta/src/GCTABackgroundPerfTable.cpp index 6bfe2cb..ec8b4a4 100644 --- a/inst/cta/src/GCTABackgroundPerfTable.cpp +++ b/inst/cta/src/GCTABackgroundPerfTable.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTABackgroundPerfTable.cpp - CTA performance table background class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -33,6 +33,7 @@ #include "GMath.hpp" #include "GIntegral.hpp" #include "GRan.hpp" +#include "GCTAInstDir.hpp" #include "GCTABackgroundPerfTable.hpp" /* __ Method name definitions ____________________________________________ */ @@ -45,7 +46,7 @@ /* __ Debug definitions __________________________________________________ */ //#define G_DEBUG_MC //!< Debug Monte Carlo method -#define G_LOG_INTERPOLATION //!< Energy interpolate log10(background rate) +#define G_LOG_INTERPOLATION //!< Energy interpolate log(background rate) /* __ Constants __________________________________________________________ */ @@ -163,27 +164,23 @@ GCTABackgroundPerfTable& GCTABackgroundPerfTable::operator=(const GCTABackground * @brief Return background rate in units of events/s/MeV/sr * * @param[in] logE Log10 of the true photon energy (TeV). - * @param[in] detx Tangential coord in nominal sys (rad). - * @param[in] dety Tangential coord in nominal sys (rad). + * @param[in] detx Tangential coordinate in nominal system (rad). + * @param[in] dety Tangential coordinate in nominal system (rad). + * @return Background rate (events/s/MeV/sr) * * Returns the background rate in units of events/s/MeV/sr for a given energy * and detector coordinates. The method assures that the background rate * never becomes negative. + * + * If the performance table contains less than 2 nodes an empty value is + * returned. ***************************************************************************/ double GCTABackgroundPerfTable::operator()(const double& logE, const double& detx, const double& dety) const { - // Get background rate - double rate = m_logE.interpolate(logE, m_background); - #if defined(G_LOG_INTERPOLATION) - rate = std::pow(10.0, rate); - #endif - - // Make sure that background rate is not negative - if (rate < 0.0) { - rate = 0.0; - } + // Initialise rate + double rate = this->rate(logE); // Optionally add in Gaussian offset angle dependence if (m_sigma != 0.0) { @@ -192,7 +189,7 @@ double GCTABackgroundPerfTable::operator()(const double& logE, double scale = std::exp(-0.5 * arg * arg); rate *= scale; } - + // Return background rate return rate; } @@ -248,8 +245,10 @@ GCTABackgroundPerfTable* GCTABackgroundPerfTable::clone(void) const void GCTABackgroundPerfTable::load(const GFilename& filename) { // Clear arrays - m_logE.clear(); + m_energy.clear(); + m_log10_energy.clear(); m_background.clear(); + m_log_background.clear(); // Allocate line buffer const int n = 1000; @@ -297,13 +296,15 @@ void GCTABackgroundPerfTable::load(const GFilename& filename) bgrate = 0.0; } - // Push elements in node array and vector - m_logE.append(logE); - #if defined(G_LOG_INTERPOLATION) - m_background.push_back(std::log10(bgrate)); - #else + // Set energy + GEnergy energy; + energy.log10TeV(logE); + + // Push elements in node array and vectors + m_energy.append(energy); + m_log10_energy.append(logE); m_background.push_back(bgrate); - #endif + m_log_background.push_back(std::log(bgrate)); } // endwhile: looped over lines @@ -376,9 +377,96 @@ GCTAInstDir GCTABackgroundPerfTable::mc(const GEnergy& energy, /***********************************************************************//** + * @brief Returns background count rate integrated over energy interval + * + * @param[in] dir Instrument direction. + * @param[in] emin Minimum energy of energy interval. + * @param[in] emax Maximum energy of energy interval. + * @return Integrated background count rate (counts/src/sr). + * + * Returns the background count rate for a given instrument direction that + * is integrated over a specified energy interval. + * + * If the energy interval is not positive, a zero background rate is + * returned. + ***************************************************************************/ +double GCTABackgroundPerfTable::rate_ebin(const GCTAInstDir& dir, + const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise rate + double rate = 0.0; + + // Continue only if energy interval is positive + if (emax > emin) { + + // Initialise first and second node + double x1 = emin.MeV(); + double x2 = 0.0; + double f1 = this->rate(emin.log10TeV()); + double f2 = 0.0; + + // Loop over all nodes + for (int i = 0; i < size(); ++i) { + + // If node energy is below x1 then skip node + if (m_energy[i].MeV() <= x1) { + continue; + } + + // If node energy is above emax then use emax as energy + if (m_energy[i] > emax) { + x2 = emax.MeV(); + f2 = this->rate(emax.log10TeV()); + } + + // ... otherwise use node energy + else { + x2 = m_energy[i].MeV(); + f2 = m_background[i]; + } + + // Compute integral + rate += gammalib::plaw_integral(x1, f1, x2, f2); + + // Set second node as first node + x1 = x2; + f1 = f2; + + // If node energy is above emax then break now + if (m_energy[i] > emax) { + break; + } + + } // endfor: looped over all nodes + + // If last node energy is below emax then compute last part of + // integral up to emax + if (x1 < emax.MeV()) { + x2 = emax.MeV(); + f2 = this->rate(emax.log10TeV()); + rate += gammalib::plaw_integral(x1, f1, x2, f2); + } + + // Optionally add in Gaussian offset angle dependence + if (m_sigma != 0.0) { + double theta = dir.theta() * gammalib::rad2deg; + double arg = theta * theta / m_sigma; + double scale = std::exp(-0.5 * arg * arg); + rate *= scale; + } + + } // endif: energy interval was positive + + // Return background rate + return rate; +} + + +/***********************************************************************//** * @brief Print background information * - * @param[in] chatter Chattiness (defaults to NORMAL). + * @param[in] chatter Chattiness. * @return String containing background information. ***************************************************************************/ std::string GCTABackgroundPerfTable::print(const GChatter& chatter) const @@ -390,8 +478,8 @@ std::string GCTABackgroundPerfTable::print(const GChatter& chatter) const if (chatter != SILENT) { // Compute energy boundaries in TeV - double emin = std::pow(10.0, m_logE[0]); - double emax = std::pow(10.0, m_logE[size()-1]); + double emin = (size() > 0) ? m_energy[0].TeV() : 0.0; + double emax = (size() > 1) ? m_energy[size()-1].TeV() : 0.0; // Append header result.append("=== GCTABackgroundPerfTable ==="); @@ -400,7 +488,7 @@ std::string GCTABackgroundPerfTable::print(const GChatter& chatter) const result.append("\n"+gammalib::parformat("Filename")+m_filename); result.append("\n"+gammalib::parformat("Number of energy bins") + gammalib::str(size())); - result.append("\n"+gammalib::parformat("Log10(Energy) range")); + result.append("\n"+gammalib::parformat("Energy range")); result.append(gammalib::str(emin) + " - " + gammalib::str(emax) + @@ -437,8 +525,10 @@ void GCTABackgroundPerfTable::init_members(void) { // Initialise members m_filename.clear(); - m_logE.clear(); + m_energy.clear(); + m_log10_energy.clear(); m_background.clear(); + m_log_background.clear(); m_sigma = 3.0; // Initialise MC cache @@ -457,10 +547,12 @@ void GCTABackgroundPerfTable::init_members(void) void GCTABackgroundPerfTable::copy_members(const GCTABackgroundPerfTable& bgd) { // Copy members - m_filename = bgd.m_filename; - m_logE = bgd.m_logE; - m_background = bgd.m_background; - m_sigma = bgd.m_sigma; + m_filename = bgd.m_filename; + m_energy = bgd.m_energy; + m_log10_energy = bgd.m_log10_energy; + m_background = bgd.m_background; + m_log_background = bgd.m_log_background; + m_sigma = bgd.m_sigma; // Copy MC cache m_mc_spectrum = bgd.m_mc_spectrum; @@ -539,20 +631,12 @@ void GCTABackgroundPerfTable::init_mc_cache(void) const // Loop over nodes for (int i = 0; i < size(); ++i) { - // Set energy - GEnergy energy; - energy.log10TeV(m_logE[i]); - // Compute total rate - #if defined(G_LOG_INTERPOLATION) - double total_rate = std::pow(10.0, m_background[i]) * solidangle; - #else double total_rate = m_background[i] * solidangle; - #endif // Set node if (total_rate > 0.0) { - m_mc_spectrum.append(energy, total_rate); + m_mc_spectrum.append(m_energy[i], total_rate); } } @@ -560,3 +644,36 @@ void GCTABackgroundPerfTable::init_mc_cache(void) const // Return return; } + + +/***********************************************************************//** + * @brief Return background rate for a given energy (events/s/MeV/sr) + * + * @param[in] logE Log10 of the true photon energy (TeV). + * @return Background rate (events/s/MeV/sr) + ***************************************************************************/ +double GCTABackgroundPerfTable::rate(const double& logE) const +{ + // Initialise rate + double rate = 0.0; + + // Continue only if there are at least two nodes + if (size() > 1) { + + // Get background rate + #if defined(G_LOG_INTERPOLATION) + rate = std::exp(m_log10_energy.interpolate(logE, m_log_background)); + #else + rate = m_log10_energy.interpolate(logE, m_background); + #endif + + // Make sure that background rate is not negative + if (rate < 0.0) { + rate = 0.0; + } + + } // endif: there were enough rates + + // Return + return rate; +} diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index a7fca15..176a87d 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -64,9 +64,10 @@ const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml"; const std::string cta_caldb_king = cta_caldb+"/data/cta/e/bcf/IFAE20120510_50h_King"; const std::string cta_irf_king = "irf_file.fits"; const std::string cta_psf_table = caldbdir+"/psf_table.fits[PSF_2D_TABLE]"; -const std::string cta_edisp_perf = caldbdir+"/cta_dummy_irf.dat"; +const std::string cta_perf_table = caldbdir+"/cta_dummy_irf.dat"; const std::string cta_edisp_rmf = caldbdir+"/dc1/rmf.fits"; const std::string cta_edisp_2D = caldbdir+"/edisp_matrix.fits"; +const std::string cta_bgd_3D = caldbdir+"/edisp_matrix.fits"; const std::string cta_modbck_fit = datadir+"/bg_test.fits"; const std::string cta_point_table = datadir+"/crab_pointing.fits"; @@ -84,6 +85,9 @@ const std::string cta_onoff_obs = datadir+"/onoff_obs.xml"; const std::string cta_onoff_model = datadir+"/onoff_model.xml"; const std::string cta_onoff_onreg = datadir+"/onoff_region_on.reg"; +/* __ Debug definitions __________________________________________________ */ +//#define G_COMPUTE_REF_RATE_EBIN //!< Compute reference rate + /***********************************************************************//** * @brief Set miscellaneous CTA classes test methods @@ -129,26 +133,30 @@ void TestGCTAResponse::set(void) "Test Table point spread function"); append(static_cast(&TestGCTAResponse::test_response_npsf), "Test integrated point spread function"); - append(static_cast(&TestGCTAResponse::test_response_edisp), - "Test energy dispersion"); - append(static_cast(&TestGCTAResponse::test_response_edisp_PerfTable), - "Test Performance Table energy dispersion"); - append(static_cast(&TestGCTAResponse::test_response_edisp_RMF), - "Test RMF energy dispersion"); - append(static_cast(&TestGCTAResponse::test_response_edisp_2D), - "Test 2D energy dispersion"); append(static_cast(&TestGCTAResponse::test_response_irf_diffuse), "Test diffuse IRF"); append(static_cast(&TestGCTAResponse::test_response_npred_diffuse), "Test diffuse IRF integration"); + append(static_cast(&TestGCTAResponse::test_response_edisp), + "Test GCTAEdisp class"); + append(static_cast(&TestGCTAResponse::test_response_edisp_PerfTable), + "Test GCTAEdispPerfTable class"); + append(static_cast(&TestGCTAResponse::test_response_edisp_RMF), + "Test GCTAEdispRmf class"); + append(static_cast(&TestGCTAResponse::test_response_edisp_2D), + "Test GCTAEdisp2D class"); + append(static_cast(&TestGCTAResponse::test_response_bgd_PerfTable), + "Test GCTABackgroundPerfTable class"); + append(static_cast(&TestGCTAResponse::test_response_bgd_3D), + "Test GCTABackground3D class"); append(static_cast(&TestGCTAResponse::test_response_expcube), - "Test exposure cube"); + "Test GCTACubeExposure class"); append(static_cast(&TestGCTAResponse::test_response_psfcube), - "Test point spread function cube"); + "Test GCTACubePsf class"); append(static_cast(&TestGCTAResponse::test_response_edispcube), - "Test energy dispersion cube"); + "Test GCTACubeEdisp class"); append(static_cast(&TestGCTAResponse::test_response_bkgcube), - "Test background cube"); + "Test GCTACubeBackground class"); // Return return; @@ -608,6 +616,165 @@ void TestGCTAResponse::test_response_npsf(void) return; } + +/***********************************************************************//** + * @brief Test CTA IRF computation for diffuse source model + * + * Tests the IRF computation for the diffuse source model. This is done + * by calling the GCTAObservation::model method which in turn calls the + * GCTAResponse::irf_diffuse method. The test is done for a small counts + * map to keep the test executing reasonably fast. + ***************************************************************************/ +void TestGCTAResponse::test_response_irf_diffuse(void) +{ + // Set reference value + //double ref = 13803.800313356; + //const double ref = 13803.5186374; + //const double ref = 13803.6932774; // after GWcs::solidangle improvement + //const double ref = 13803.3453994; // after Irf computation improvements + //const double ref = 13803.3461768; // using skymap flux() method + const double ref = 13803.5251338; // after correcting contains() method + + // Set parameters + double src_ra = 201.3651; + double src_dec = -43.0191; + int nebins = 5; + + // Setup pointing on Cen A + GSkyDir skyDir; + skyDir.radec_deg(src_ra, src_dec); + GCTAPointing pnt; + pnt.dir(skyDir); + + // Setup skymap (10 energy layers) + GSkyMap map("CAR", "CEL", src_ra, src_dec, 0.5, 0.5, 10, 10, nebins); + + // Setup time interval + GGti gti; + GTime tstart(0.0); + GTime tstop(1800.0); + gti.append(tstart, tstop); + + // Setup energy boundaries + GEnergy emin(0.1, "TeV"); + GEnergy emax(100.0, "TeV"); + GEbounds ebounds(nebins, emin, emax); + + // Setup event cube centreed on Cen A + GCTAEventCube cube(map, ebounds, gti); + + // Setup dummy CTA observation + GCTAObservation obs; + obs.ontime(1800.0); + obs.livetime(1600.0); + obs.deadc(1600.0/1800.0); + obs.response(cta_irf, GCaldb(cta_caldb)); + obs.events(cube); + obs.pointing(pnt); + + // Load model for IRF computation + GModels models(cta_rsp_xml); + + // Reset sum + double sum = 0.0; + + // Iterate over all bins in event cube + for (int i = 0; i < obs.events()->size(); ++i) { + + // Get event pointer + const GEventBin* bin = (*(static_cast(obs.events())))[i]; + + // Get model and add to sum + double model = obs.model(models, *bin, NULL) * bin->size(); + sum += model; + + } + + // Test sum + test_value(sum, ref, 0.2, "Diffuse IRF computation"); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Test CTA Npred computation + * + * Tests the Npred computation for the diffuse source model. This is done + * by loading the model from the XML file and by calling the + * GCTAObservation::npred method which in turn calls the + * GCTAResponse::npred_diffuse method. The test takes a few seconds. + ***************************************************************************/ +void TestGCTAResponse::test_response_npred_diffuse(void) +{ + // Set reference value + //const double ref = 11212.26274; // npred_spec precision of 1e-6 + //const double ref = 11212.437464; // npred_spec precision of 1e-5 + //const double ref = 11212.4370702; // After GWcs::solidangle improvement + //const double ref = 12644.3391902; // After correcting for deadtime bug + //const double ref = 12643.6159142; // npred_spec precision of 1e-6 + //const double ref = 10997.0851335; // After response computation change + //const double ref = 11212.7161901; // npred_diffuse precision to iter=8 + //const double ref = 11238.6928076; // npred_diffuse precision to iter=9 + //const double ref = 11238.6933328; // using skymap flux() method + const double ref = 11238.8390288; // after correcting contains() method + + // Set parameters + double src_ra = 201.3651; + double src_dec = -43.0191; + double roi_rad = 4.0; + + // Setup ROI centred on Cen A with a radius of 4 deg + GCTARoi roi; + GCTAInstDir instDir; + instDir.dir().radec_deg(src_ra, src_dec); + roi.centre(instDir); + roi.radius(roi_rad); + + // Setup pointing on Cen A + GSkyDir skyDir; + skyDir.radec_deg(src_ra, src_dec); + GCTAPointing pnt; + pnt.dir(skyDir); + + // Setup dummy event list + GGti gti; + GEbounds ebounds; + GTime tstart(0.0); + GTime tstop(1800.0); + GEnergy emin(0.1, "TeV"); + GEnergy emax(100.0, "TeV"); + gti.append(tstart, tstop); + ebounds.append(emin, emax); + GCTAEventList events; + events.roi(roi); + events.gti(gti); + events.ebounds(ebounds); + + // Setup dummy CTA observation + GCTAObservation obs; + obs.ontime(1800.0); + obs.livetime(1600.0); + obs.deadc(1600.0/1800.0); + obs.response(cta_irf, GCaldb(cta_caldb)); + obs.events(events); + obs.pointing(pnt); + + // Load models for Npred computation + GModels models(cta_rsp_xml); + + // Perform Npred computation + double npred = obs.npred(models, NULL); + + // Test Npred + test_value(npred, ref, 0.2, "Diffuse Npred computation"); + + // Return + return; +} + + /***********************************************************************//** * @brief Test CTA Energy Dispersion computation * @@ -641,7 +808,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void) // Test GCTAEdispPerfTable constructor test_try("GCTAEdispPerfTable constructor"); try { - GCTAEdispPerfTable test(cta_edisp_perf); + GCTAEdispPerfTable test(cta_perf_table); test_try_success(); } catch (std::exception &e) { @@ -649,7 +816,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void) } // Load energy dispersion - GCTAEdispPerfTable edisp(cta_edisp_perf); + GCTAEdispPerfTable edisp(cta_perf_table); // Test ebounds_obs and ebounds_src methods test_try("GCTAEdispPerfTable ebounds_obs and ebounds_src methods"); @@ -687,8 +854,6 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void) } - - /***********************************************************************//** * @brief Test CTA RMF energy dispersion class ***************************************************************************/ @@ -783,80 +948,99 @@ void TestGCTAResponse::test_response_edisp_2D(void) /***********************************************************************//** - * @brief Test CTA IRF computation for diffuse source model - * - * Tests the IRF computation for the diffuse source model. This is done - * by calling the GCTAObservation::model method which in turn calls the - * GCTAResponse::irf_diffuse method. The test is done for a small counts - * map to keep the test executing reasonably fast. + * @brief Test GCTABackgroundPerfTable class ***************************************************************************/ -void TestGCTAResponse::test_response_irf_diffuse(void) +void TestGCTAResponse::test_response_bgd_PerfTable(void) { - // Set reference value - //double ref = 13803.800313356; - //const double ref = 13803.5186374; - //const double ref = 13803.6932774; // after GWcs::solidangle improvement - //const double ref = 13803.3453994; // after Irf computation improvements - //const double ref = 13803.3461768; // using skymap flux() method - const double ref = 13803.5251338; // after correcting contains() method - - // Set parameters - double src_ra = 201.3651; - double src_dec = -43.0191; - int nebins = 5; - - // Setup pointing on Cen A - GSkyDir skyDir; - skyDir.radec_deg(src_ra, src_dec); - GCTAPointing pnt; - pnt.dir(skyDir); - - // Setup skymap (10 energy layers) - GSkyMap map("CAR", "CEL", src_ra, src_dec, 0.5, 0.5, 10, 10, nebins); - - // Setup time interval - GGti gti; - GTime tstart(0.0); - GTime tstop(1800.0); - gti.append(tstart, tstop); - - // Setup energy boundaries - GEnergy emin(0.1, "TeV"); - GEnergy emax(100.0, "TeV"); - GEbounds ebounds(nebins, emin, emax); - - // Setup event cube centreed on Cen A - GCTAEventCube cube(map, ebounds, gti); - - // Setup dummy CTA observation - GCTAObservation obs; - obs.ontime(1800.0); - obs.livetime(1600.0); - obs.deadc(1600.0/1800.0); - obs.response(cta_irf, GCaldb(cta_caldb)); - obs.events(cube); - obs.pointing(pnt); - - // Load model for IRF computation - GModels models(cta_rsp_xml); - - // Reset sum - double sum = 0.0; - - // Iterate over all bins in event cube - for (int i = 0; i < obs.events()->size(); ++i) { - - // Get event pointer - const GEventBin* bin = (*(static_cast(obs.events())))[i]; - - // Get model and add to sum - double model = obs.model(models, *bin, NULL) * bin->size(); - sum += model; - + // Set reference values + const double ref_rate_below = 0.000612295557409636; + const double ref_rate_above = 5.74072516267095e-11; + const double ref_rate_onaxis = 3.62782669235078e-05; + const double ref_rate_offaxis = 3.12350491143832e-05; + const double ref_rate_ebin_within = 44.9181299290927; + const double ref_rate_ebin_below = 3676.68195289697; + const double ref_rate_ebin_above = 76.2328310391366; + + // Set energies for energy range for rate integration + const GEnergy ebelow(10.0, "GeV"); + const GEnergy emin(1.0, "TeV"); + const GEnergy emax(10.0, "TeV"); + const GEnergy eabove(200.0, "TeV"); + + // Test GCTABackgroundPerfTable constructor + test_try("GCTABackgroundPerfTable constructor"); + try { + GCTABackgroundPerfTable bgd1(cta_perf_table); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); } - // Test sum - test_value(sum, ref, 0.2, "Diffuse IRF computation"); + // Check content of empty instance + GCTABackgroundPerfTable bgd1; + test_value(bgd1.size(), 0, "Nodes of empty instance"); + test_value(bgd1.sigma(), 3.0, "Sigma value of empty instance"); + test_value(bgd1.filename(), "", "Filename of empty instance"); + test_value(bgd1(0.0, 0.0, 0.0), 0.0, + "Onaxis operator value of empty instance"); + test_value(bgd1(0.0, 0.01, 0.02), 0.0, + "Offaxis operator value of empty instance"); + test_value(bgd1.rate_ebin(GCTAInstDir(), emin, emax), 0.0, + "Integrated rate value of empty instance"); + + // Check content of filled instance + GCTABackgroundPerfTable bgd2(cta_perf_table); + test_value(bgd2.size(), 20, "Nodes of filled instance"); + test_value(bgd2.sigma(), 3.0, "Sigma value of filled instance"); + test_value(bgd2.filename(), cta_perf_table, "Filename of filled instance"); + test_value(bgd2(-2.0, 0.0, 0.0), ref_rate_below, + "Operator value of filled instance extrapolated below lowest energy"); + test_value(bgd2(3.0, 0.0, 0.0), ref_rate_above, + "Operator value of filled instance extrapolated above highest energy"); + test_value(bgd2(0.0, 0.0, 0.0), ref_rate_onaxis, + "Onaxis operator value of filled instance"); + test_value(bgd2(0.0, 0.01, 0.02), ref_rate_offaxis, + "Offaxis operator value of filled instance"); + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within, + "Integrated rate value of filled instance (within covered energies)"); + test_value(bgd2.rate_ebin(GCTAInstDir(), ebelow, emax), ref_rate_ebin_below, + "Integrated rate value of filled instance (extrapolate below lowest energy)"); + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, eabove), ref_rate_ebin_above, + "Integrated rate value of filled instance (extrapolate above highest energy)"); + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emin), 0.0, + "Integrated rate value of filled instance for zero interval"); + + // Check content of copied instance + GCTABackgroundPerfTable bgd3(bgd2); + test_value(bgd3.size(), 20, "Nodes of copied instance"); + test_value(bgd3.sigma(), 3.0, "Sigma value of copied instance"); + test_value(bgd3.filename(), cta_perf_table, "Filename of copied instance"); + test_value(bgd3(-2.0, 0.0, 0.0), ref_rate_below, + "Operator value of copied instance extrapolated below lowest energy"); + test_value(bgd3(3.0, 0.0, 0.0), ref_rate_above, + "Operator value of copied instance extrapolated above highest energy"); + test_value(bgd3(0.0, 0.0, 0.0), ref_rate_onaxis, + "Onaxis operator value of copied instance"); + test_value(bgd3(0.0, 0.01, 0.02), ref_rate_offaxis, + "Offaxis operator value of copied instance"); + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within, + "Integrated rate value of copied instance"); + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emin), 0.0, + "Integrated rate value of copied instance for zero interval"); + + // Compute reference rate + #if defined(G_COMPUTE_REF_RATE_EBIN) + double rate = 0.0; + int nbins = 10000; + GEnergy E = emin; + GEnergy dE = (eabove-emin)/nbins; + for (int i = 0; i < nbins; ++i) { + rate += bgd2(E.log10TeV(), 0.0, 0.0) * dE.MeV(); + E += dE; + } + std::cout << rate << std::endl; + #endif // Return return; @@ -864,76 +1048,96 @@ void TestGCTAResponse::test_response_irf_diffuse(void) /***********************************************************************//** - * @brief Test CTA Npred computation - * - * Tests the Npred computation for the diffuse source model. This is done - * by loading the model from the XML file and by calling the - * GCTAObservation::npred method which in turn calls the - * GCTAResponse::npred_diffuse method. The test takes a few seconds. + * @brief Test GCTABackground3D class ***************************************************************************/ -void TestGCTAResponse::test_response_npred_diffuse(void) +void TestGCTAResponse::test_response_bgd_3D(void) { - // Set reference value - //const double ref = 11212.26274; // npred_spec precision of 1e-6 - //const double ref = 11212.437464; // npred_spec precision of 1e-5 - //const double ref = 11212.4370702; // After GWcs::solidangle improvement - //const double ref = 12644.3391902; // After correcting for deadtime bug - //const double ref = 12643.6159142; // npred_spec precision of 1e-6 - //const double ref = 10997.0851335; // After response computation change - //const double ref = 11212.7161901; // npred_diffuse precision to iter=8 - //const double ref = 11238.6928076; // npred_diffuse precision to iter=9 - //const double ref = 11238.6933328; // using skymap flux() method - const double ref = 11238.8390288; // after correcting contains() method - - // Set parameters - double src_ra = 201.3651; - double src_dec = -43.0191; - double roi_rad = 4.0; - - // Setup ROI centred on Cen A with a radius of 4 deg - GCTARoi roi; - GCTAInstDir instDir; - instDir.dir().radec_deg(src_ra, src_dec); - roi.centre(instDir); - roi.radius(roi_rad); - - // Setup pointing on Cen A - GSkyDir skyDir; - skyDir.radec_deg(src_ra, src_dec); - GCTAPointing pnt; - pnt.dir(skyDir); - - // Setup dummy event list - GGti gti; - GEbounds ebounds; - GTime tstart(0.0); - GTime tstop(1800.0); - GEnergy emin(0.1, "TeV"); - GEnergy emax(100.0, "TeV"); - gti.append(tstart, tstop); - ebounds.append(emin, emax); - GCTAEventList events; - events.roi(roi); - events.gti(gti); - events.ebounds(ebounds); - - // Setup dummy CTA observation - GCTAObservation obs; - obs.ontime(1800.0); - obs.livetime(1600.0); - obs.deadc(1600.0/1800.0); - obs.response(cta_irf, GCaldb(cta_caldb)); - obs.events(events); - obs.pointing(pnt); - - // Load models for Npred computation - GModels models(cta_rsp_xml); + // Set reference values + const double ref_rate_below = 0.0; + const double ref_rate_above = 0.0; + const double ref_rate_onaxis = 9.7632195934345e-05; + const double ref_rate_offaxis = 13.5162501643753e-05; + const double ref_rate_ebin_within = 65.6069433285784; + const double ref_rate_ebin_below = 4647.37614735018; + const double ref_rate_ebin_above = 69.5331125123672; + + // Set energies for energy range for rate integration + const GEnergy ebelow(10.0, "GeV"); + const GEnergy emin(1.0, "TeV"); + const GEnergy emax(10.0, "TeV"); + const GEnergy eabove(200.0, "TeV"); + + // Test GCTABackground3D constructor + test_try("GCTABackground3D constructor"); + try { + GCTABackground3D bgd1(cta_bgd_3D); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } - // Perform Npred computation - double npred = obs.npred(models, NULL); + // Check content of empty instance + GCTABackground3D bgd1; + test_value(bgd1.table().axes(), 0, "Response table axes of empty instance"); + test_value(bgd1.filename(), "", "Filename of empty instance"); + test_value(bgd1(0.0, 0.0, 0.0), 0.0, + "Onaxis operator value of empty instance"); + test_value(bgd1(0.0, 0.01, 0.02), 0.0, + "Offaxis operator value of empty instance"); + test_value(bgd1.rate_ebin(GCTAInstDir(), emin, emax), 0.0, + "Integrated rate value of empty instance"); - // Test Npred - test_value(npred, ref, 0.2, "Diffuse Npred computation"); + // Check content of filled instance + GCTABackground3D bgd2(cta_bgd_3D); + test_value(bgd2.table().axes(), 3, "Response table axes of filled instance"); + test_value(bgd2.filename(), cta_bgd_3D, "Filename of filled instance"); + test_value(bgd2(-2.0, 0.0, 0.0), ref_rate_below, + "Operator value of filled instance extrapolated below lowest energy"); + test_value(bgd2(3.0, 0.0, 0.0), ref_rate_above, + "Operator value of filled instance extrapolated above highest energy"); + test_value(bgd2(0.0, 0.0, 0.0), ref_rate_onaxis, + "Onaxis operator value of filled instance"); + test_value(bgd2(0.0, 0.01, 0.02), ref_rate_offaxis, + "Offaxis operator value of filled instance"); + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within, + "Integrated rate value of filled instance"); + test_value(bgd2.rate_ebin(GCTAInstDir(), ebelow, emax), ref_rate_ebin_below, + "Integrated rate value of filled instance (extrapolate below lowest energy)"); + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, eabove), ref_rate_ebin_above, + "Integrated rate value of filled instance (extrapolate above highest energy)"); + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emin), 0.0, + "Integrated rate value of filled instance for zero interval"); + + // Check content of copied instance + GCTABackground3D bgd3(bgd2); + test_value(bgd3.table().axes(), 3, "Response table axes of copied instance"); + test_value(bgd3.filename(), cta_bgd_3D, "Filename of copied instance"); + test_value(bgd3(-2.0, 0.0, 0.0), ref_rate_below, + "Operator value of copied instance extrapolated below lowest energy"); + test_value(bgd3(3.0, 0.0, 0.0), ref_rate_above, + "Operator value of copied instance extrapolated above highest energy"); + test_value(bgd3(0.0, 0.0, 0.0), ref_rate_onaxis, + "Onaxis operator value of copied instance"); + test_value(bgd3(0.0, 0.01, 0.02), ref_rate_offaxis, + "Offaxis operator value of copied instance"); + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within, + "Integrated rate value of copied instance"); + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emin), 0.0, + "Integrated rate value of copied instance for zero interval"); + + // Compute reference rate + #if defined(G_COMPUTE_REF_RATE_EBIN) + double rate = 0.0; + int nbins = 10000; + GEnergy E = emin; + GEnergy dE = (eabove-E)/nbins; + for (int i = 0; i < nbins; ++i) { + rate += bgd2(E.log10TeV(), 0.0, 0.0) * dE.MeV(); + E += dE; + } + std::cout << rate << std::endl; + #endif // Return return; diff --git a/inst/cta/test/test_CTA.hpp b/inst/cta/test/test_CTA.hpp index f21aa0f..034234a 100644 --- a/inst/cta/test/test_CTA.hpp +++ b/inst/cta/test/test_CTA.hpp @@ -78,12 +78,14 @@ public: void test_response_psf_king(void); void test_response_psf_table(void); void test_response_npsf(void); + void test_response_irf_diffuse(void); + void test_response_npred_diffuse(void); void test_response_edisp(void); void test_response_edisp_PerfTable(void); void test_response_edisp_RMF(void); void test_response_edisp_2D(void); - void test_response_irf_diffuse(void); - void test_response_npred_diffuse(void); + void test_response_bgd_PerfTable(void); + void test_response_bgd_3D(void); void test_response_expcube(void); void test_response_psfcube(void); void test_response_bkgcube(void);