Samuel Flis / gammalib

Commit 5e4a21b944bd127331ca99122cfc72ac746d1d84

Authored by Jürgen Knödlseder
1 parent 63e6870b

Add GCTABackground::rate_ebin() method (#2308)

Inline Side-by-side
Showing 13 changed files with 801 additions and 321 deletions
ChangeLog
  View file @5e4a21b
1   -2018-02-06
  1 +2018-02-09
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + Add GCTABackground::rate_ebin() method (#2308)
  7 +
6 8  
7 9 2018-01-23
8 10  
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -8 February 2018
  3 +9 February 2018
4 4  
5 5  
6 6 1. Introduction
... ... @@ -25,6 +25,9 @@ The following methods have been added:
25 25 - gammalib::plaw_integral(double&, double&, double&, double&)
26 26 - GCTAInstDir::theta()
27 27 - GCTAInstDir::phi()
  28 +- GCTABackground::rate_ebin()
  29 +- GCTABackgroundPerfTable::rate_ebin()
  30 +- GCTABackground3D::rate_ebin()
28 31  
29 32 The following methods have been removed:
30 33 - none
... ... @@ -133,6 +136,11 @@ None
133 136 A GCTAInstDir::theta() and GCTAInstDir::phi() methods were added that return the
134 137 distance to the camera centre in radians and the azimuth angle, respectively (#2308).
135 138  
  139 +A rate_ebin() method was added to the interface of the GCTABackground class and
  140 +implemented for the GCTABackgroundPerfTable and GCTABackground3D classes. The method
  141 +computes for a given instrument direction the background rate integrated over an
  142 +energy bin (#2308).
  143 +
136 144  
137 145 21. COMPTEL interface
138 146 ---------------------
... ...
doc/source/admin/release_history/1.6.rst
  View file @5e4a21b
... ... @@ -22,7 +22,8 @@ Bug fixes
22 22 Improvements
23 23 ------------
24 24  
25   -* ...
  25 +* [`2308 <https://cta-redmine.irap.omp.eu/issues/2308>`_] -
  26 + Add `GCTABackground::rate_ebin()` method
26 27  
27 28  
28 29 New features
... ...
inst/cta/include/GCTABackground.hpp
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackground.hpp - CTA background model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -30,11 +30,14 @@
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
32 32 #include "GBase.hpp"
33   -#include "GFilename.hpp"
34   -#include "GModelSpectralNodes.hpp"
35   -#include "GCTAInstDir.hpp"
36 33  
37 34 /* __ Forward declarations _______________________________________________ */
  35 +class GRan;
  36 +class GEnergy;
  37 +class GTime;
  38 +class GFilename;
  39 +class GCTAInstDir;
  40 +class GModelSpectralNodes;
38 41  
39 42  
40 43 /***********************************************************************//**
... ... @@ -68,6 +71,9 @@ public:
68 71 const GTime& time,
69 72 GRan& ran) const = 0;
70 73 virtual const GModelSpectralNodes& spectrum(void) const = 0;
  74 + virtual double rate_ebin(const GCTAInstDir& dir,
  75 + const GEnergy& emin,
  76 + const GEnergy& emax) const = 0;
71 77 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
72 78  
73 79 protected:
... ...
inst/cta/include/GCTABackground3D.hpp
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackground3D.hpp - CTA 3D background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -29,7 +29,9 @@
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
  32 +#include "GEnergies.hpp"
32 33 #include "GFilename.hpp"
  34 +#include "GModelSpectralNodes.hpp"
33 35 #include "GCTABackground.hpp"
34 36 #include "GCTAResponseTable.hpp"
35 37  
... ... @@ -75,6 +77,9 @@ public:
75 77 const GTime& time,
76 78 GRan& ran) const;
77 79 const GModelSpectralNodes& spectrum(void) const;
  80 + double rate_ebin(const GCTAInstDir& dir,
  81 + const GEnergy& emin,
  82 + const GEnergy& emax) const;
78 83 std::string print(const GChatter& chatter = NORMAL) const;
79 84  
80 85 // Methods
... ... @@ -91,17 +96,19 @@ private:
91 96 void init_members(void);
92 97 void copy_members(const GCTABackground3D& bgd);
93 98 void free_members(void);
94   - void set_limits(void);
  99 + void set_members(void);
95 100 int index(const int& idetx, const int& idety, const int& iebin) const;
96 101 void init_mc_cache(void) const;
97 102 void init_mc_max_rate(void) const;
98 103 double solid_angle(const double& detx1, const double& dety1,
99 104 const double& detx2, const double& dety2,
100 105 const double& detx3, const double& dety3) const;
  106 + double rate(const int& iebin, const double& detx, const double& dety) const;
101 107  
102 108 // Members
103 109 GFilename m_filename; //!< Name of background response file
104 110 GCTAResponseTable m_background; //!< Background response table
  111 + GEnergies m_energy; //!< Vector of energies
105 112 int m_inx_detx; //!< DETX index
106 113 int m_inx_dety; //!< DETY index
107 114 int m_inx_energy; //!< Energy index
... ... @@ -199,6 +206,7 @@ inline
199 206 void GCTABackground3D::table(const GCTAResponseTable& table)
200 207 {
201 208 m_background = table;
  209 + set_members();
202 210 }
203 211  
204 212 #endif /* GCTABACKGROUND3D_HPP */
... ...
inst/cta/include/GCTABackgroundPerfTable.hpp
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackgroundPerfTable.hpp - CTA performance table background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -30,8 +30,10 @@
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
32 32 #include "GFilename.hpp"
  33 +#include "GEnergies.hpp"
  34 +#include "GNodeArray.hpp"
  35 +#include "GModelSpectralNodes.hpp"
33 36 #include "GCTABackground.hpp"
34   -#include "GCTAResponseTable.hpp"
35 37  
36 38 /* __ Forward declarations _______________________________________________ */
37 39  
... ... @@ -68,6 +70,9 @@ public:
68 70 const GTime& time,
69 71 GRan& ran) const;
70 72 const GModelSpectralNodes& spectrum(void) const;
  73 + double rate_ebin(const GCTAInstDir& dir,
  74 + const GEnergy& emin,
  75 + const GEnergy& emax) const;
71 76 std::string print(const GChatter& chatter = NORMAL) const;
72 77  
73 78 // Methods
... ... @@ -82,6 +87,7 @@ private:
82 87 void free_members(void);
83 88 double solidangle(void) const;
84 89 void init_mc_cache(void) const;
  90 + double rate(const double& logE) const;
85 91  
86 92 // Radial integration class (used by solidangle() method). Note that
87 93 // the sigma parameter is given in rad^2
... ... @@ -99,10 +105,12 @@ private:
99 105 };
100 106  
101 107 // Members
102   - GFilename m_filename; //!< Name of background response file
103   - GNodeArray m_logE; //!< log(E) nodes for background interpolation
104   - std::vector<double> m_background; //!< Background rate
105   - double m_sigma; //!< Sigma for offset angle computation (0=none)
  108 + GFilename m_filename; //!< Name of background response file
  109 + GEnergies m_energy; //!< Vector of energies
  110 + GNodeArray m_log10_energy; //!< log10(E) nodes for background interpolation
  111 + std::vector<double> m_background; //!< Background rate
  112 + std::vector<double> m_log_background; //!< log(background rate)
  113 + double m_sigma; //!< Sigma for offset angle computation (0=none)
106 114  
107 115 // Monte Carlo cache
108 116 mutable GModelSpectralNodes m_mc_spectrum; //!< Background spectrum
... ... @@ -142,7 +150,7 @@ GFilename GCTABackgroundPerfTable::filename(void) const
142 150 inline
143 151 int GCTABackgroundPerfTable::size(void) const
144 152 {
145   - return (m_logE.size());
  153 + return (m_energy.size());
146 154 }
147 155  
148 156  
... ...
inst/cta/pyext/GCTABackground.i
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackground.i - CTA background model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -57,6 +57,9 @@ public:
57 57 const GTime& time,
58 58 GRan& ran) const = 0;
59 59 virtual const GModelSpectralNodes& spectrum(void) const = 0;
  60 + virtual double rate_ebin(const GCTAInstDir& dir,
  61 + const GEnergy& emin,
  62 + const GEnergy& emax) const = 0;
60 63 };
61 64  
62 65  
... ...
inst/cta/pyext/GCTABackground3D.i
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackground3D.i - CTA 3D background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -58,6 +58,9 @@ public:
58 58 const GTime& time,
59 59 GRan& ran) const;
60 60 const GModelSpectralNodes& spectrum(void) const;
  61 + double rate_ebin(const GCTAInstDir& dir,
  62 + const GEnergy& emin,
  63 + const GEnergy& emax) const;
61 64  
62 65 // Methods
63 66 bool is_valid(void) const;
... ...
inst/cta/pyext/GCTABackgroundPerfTable.i
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackgroundPerfTable.i - CTA performance table background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -58,6 +58,9 @@ public:
58 58 const GTime& time,
59 59 GRan& ran) const;
60 60 const GModelSpectralNodes& spectrum(void) const;
  61 + double rate_ebin(const GCTAInstDir& dir,
  62 + const GEnergy& emin,
  63 + const GEnergy& emax) const;
61 64  
62 65 // Methods
63 66 int size(void) const;
... ...
inst/cta/src/GCTABackground3D.cpp
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackground3D.cpp - CTA 3D background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,11 +35,13 @@
35 35 #include "GRan.hpp"
36 36 #include "GFits.hpp"
37 37 #include "GFitsBinTable.hpp"
  38 +#include "GCTAInstDir.hpp"
38 39 #include "GCTABackground3D.hpp"
39 40  
40 41 /* __ Method name definitions ____________________________________________ */
41   -#define G_READ "GCTABackground3D::read(GFitsTable&)"
42 42 #define G_MC "GCTABackground3D::mc(GEnergy&, GTime&, GRan&)"
  43 +#define G_SET_MEMBERS "GCTABackground3D::set_members()"
  44 +#define G_INIT_MC_CACHE "GCTABackground3D::init_mc_cache()"
43 45  
44 46 /* __ Macros _____________________________________________________________ */
45 47  
... ... @@ -195,68 +197,30 @@ double GCTABackground3D::operator()(const double&amp; logE,
195 197 if (is_valid() && (detx >= m_detx_min) && (detx <= m_detx_max) &&
196 198 (dety >= m_dety_min) && (dety <= m_dety_max)) {
197 199  
198   - // Retrieve references to node arrays
199   - const GNodeArray& detx_nodes = m_background.axis_nodes(m_inx_detx);
200   - const GNodeArray& dety_nodes = m_background.axis_nodes(m_inx_dety);
  200 + // Retrieve references to energy node array
201 201 const GNodeArray& energy_nodes = m_background.axis_nodes(m_inx_energy);
202 202  
203   - // Set values for node arrays
204   - detx_nodes.set_value(detx);
205   - dety_nodes.set_value(dety);
  203 + // Set energy for node array interpolation
206 204 energy_nodes.set_value(logE);
207 205  
208   - // Set weighting factors for bi-linear interpolation in DETX-DETY plane
209   - double wgt_ll = detx_nodes.wgt_left() * dety_nodes.wgt_left();
210   - double wgt_lr = detx_nodes.wgt_left() * dety_nodes.wgt_right();
211   - double wgt_rl = detx_nodes.wgt_right() * dety_nodes.wgt_left();
212   - double wgt_rr = detx_nodes.wgt_right() * dety_nodes.wgt_right();
213   -
214   - // Set weighting factors for energy interpolation
215   - double wgt_emin = energy_nodes.wgt_left();
216   - double wgt_emax = energy_nodes.wgt_right();
217   -
218   - // Get bin indices
219   - int inx_ll_emin = index(detx_nodes.inx_left(), dety_nodes.inx_left(),
220   - energy_nodes.inx_left());
221   - int inx_lr_emin = index(detx_nodes.inx_left(), dety_nodes.inx_right(),
222   - energy_nodes.inx_left());
223   - int inx_rl_emin = index(detx_nodes.inx_right(), dety_nodes.inx_left(),
224   - energy_nodes.inx_left());
225   - int inx_rr_emin = index(detx_nodes.inx_right(), dety_nodes.inx_right(),
226   - energy_nodes.inx_left());
227   - int inx_ll_emax = index(detx_nodes.inx_left(), dety_nodes.inx_left(),
228   - energy_nodes.inx_right());
229   - int inx_lr_emax = index(detx_nodes.inx_left(), dety_nodes.inx_right(),
230   - energy_nodes.inx_right());
231   - int inx_rl_emax = index(detx_nodes.inx_right(), dety_nodes.inx_left(),
232   - energy_nodes.inx_right());
233   - int inx_rr_emax = index(detx_nodes.inx_right(), dety_nodes.inx_right(),
234   - energy_nodes.inx_right());
235   -
236 206 // Bi-linearly interpolate the rates in both energy layers
237   - double rate_emin = wgt_ll * m_background(m_inx_bgd, inx_ll_emin) +
238   - wgt_lr * m_background(m_inx_bgd, inx_lr_emin) +
239   - wgt_rl * m_background(m_inx_bgd, inx_rl_emin) +
240   - wgt_rr * m_background(m_inx_bgd, inx_rr_emin);
241   - double rate_emax = wgt_ll * m_background(m_inx_bgd, inx_ll_emax) +
242   - wgt_lr * m_background(m_inx_bgd, inx_lr_emax) +
243   - wgt_rl * m_background(m_inx_bgd, inx_rl_emax) +
244   - wgt_rr * m_background(m_inx_bgd, inx_rr_emax);
  207 + double rate_emin = this->rate(energy_nodes.inx_left(), detx, dety);
  208 + double rate_emax = this->rate(energy_nodes.inx_right(), detx, dety);
245 209  
246 210 // If both rates are positive then perform a logarithmic interpolation
247 211 // in energy
248 212 if (rate_emin > 0.0 && rate_emax > 0.0) {
249   - rate = std::exp(wgt_emin * std::log(rate_emin) +
250   - wgt_emax * std::log(rate_emax));
251   - }
252   - else {
253   - rate = 0.0;
254   - }
255 213  
256   - // Make sure that background rate is not negative
257   - if (rate < 0.0) {
258   - rate = 0.0;
259   - }
  214 + // Perform energy interpolation
  215 + rate = std::exp(energy_nodes.wgt_left() * std::log(rate_emin) +
  216 + energy_nodes.wgt_right() * std::log(rate_emax));
  217 +
  218 + // Make sure that background rate is not negative
  219 + if (rate < 0.0) {
  220 + rate = 0.0;
  221 + }
  222 +
  223 + } // endif: both rates were positive
260 224  
261 225 } // endif: background model was valid
262 226  
... ... @@ -309,9 +273,6 @@ GCTABackground3D* GCTABackground3D::clone(void) const
309 273 *
310 274 * @param[in] table FITS table.
311 275 *
312   - * @exception GException::invalid_value
313   - * Response table is not three-dimensional.
314   - *
315 276 * Reads the background form the FITS @p table. The following column names
316 277 * are mandatory:
317 278 *
... ... @@ -334,41 +295,8 @@ void GCTABackground3D::read(const GFitsTable&amp; table)
334 295 // Read background table
335 296 m_background.read(table);
336 297  
337   - // Throw an exception if the table is not three-dimensional
338   - if (m_background.axes() != 3) {
339   - std::string msg = "Expected three-dimensional background "
340   - "response table but found "+
341   - gammalib::str(m_background.axes())+
342   - " dimensions. Please specify a three-dimensional "
343   - "background.";
344   - throw GException::invalid_value(G_READ, msg);
345   - }
346   -
347   - // Get mandatory indices (throw exception if not found)
348   - m_inx_detx = m_background.axis("DETX");
349   - m_inx_dety = m_background.axis("DETY");
350   - m_inx_energy = m_background.axis("ENERG");
351   - m_inx_bgd = m_background.table("BGD");
352   -
353   - // Get axes dimensions
354   - m_num_detx = m_background.axis_bins(m_inx_detx);
355   - m_num_dety = m_background.axis_bins(m_inx_dety);
356   - m_num_energy = m_background.axis_bins(m_inx_energy);
357   -
358   - // Set dimension array
359   - m_num[m_inx_detx] = m_num_detx;
360   - m_num[m_inx_dety] = m_num_dety;
361   - m_num[m_inx_energy] = m_num_energy;
362   -
363   - // Set DETX and DETY axis to radians
364   - m_background.axis_radians(m_inx_detx);
365   - m_background.axis_radians(m_inx_dety);
366   -
367   - // Set energy axis to logarithmic scale
368   - m_background.axis_log10(m_inx_energy);
369   -
370   - // Set DETX and DETY limits
371   - set_limits();
  298 + // Set class members
  299 + set_members();
372 300  
373 301 // Return
374 302 return;
... ... @@ -580,9 +508,102 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy&amp; energy,
580 508  
581 509  
582 510 /***********************************************************************//**
  511 + * @brief Returns background count rate integrated over energy interval
  512 + *
  513 + * @param[in] dir Instrument direction.
  514 + * @param[in] emin Minimum energy of energy interval.
  515 + * @param[in] emax Maximum energy of energy interval.
  516 + * @return Integrated background count rate (counts/src/sr).
  517 + *
  518 + * Returns the background count rate for a given instrument direction that
  519 + * is integrated over a specified energy interval.
  520 + *
  521 + * If the energy interval is not positive, a zero background rate is
  522 + * returned.
  523 + ***************************************************************************/
  524 +double GCTABackground3D::rate_ebin(const GCTAInstDir& dir,
  525 + const GEnergy& emin,
  526 + const GEnergy& emax) const
  527 +{
  528 + // Initialise rate
  529 + double rate = 0.0;
  530 +
  531 + // Continue only if the background template is valid, the DETX and DETY
  532 + // values are comprised in template range and the energy interval is
  533 + // positive
  534 + if (is_valid() && (dir.detx() >= m_detx_min) &&
  535 + (dir.detx() <= m_detx_max) &&
  536 + (dir.dety() >= m_dety_min) &&
  537 + (dir.dety() <= m_dety_max) &&
  538 + (emax > emin)) {
  539 +
  540 + // Retrieve references to node arrays
  541 + //const GNodeArray& detxs = m_background.axis_nodes(m_inx_detx);
  542 + //const GNodeArray& detys = m_background.axis_nodes(m_inx_dety);
  543 +
  544 + // Initialise first and second node
  545 + double x1 = emin.MeV();
  546 + double x2 = 0.0;
  547 + double f1 = (*this)(emin.log10TeV(), dir.detx(), dir.dety());
  548 + double f2 = 0.0;
  549 +
  550 + // Loop over all nodes
  551 + for (int i = 0; i < m_energy.size(); ++i) {
  552 +
  553 + // If node energy is below x1 then skip node
  554 + if (m_energy[i].MeV() <= x1) {
  555 + continue;
  556 + }
  557 +
  558 + // If node energy is above emax then use emax as energy
  559 + if (m_energy[i] > emax) {
  560 + x2 = emax.MeV();
  561 + f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety());
  562 + }
  563 +
  564 + // ... otherwise use node energy
  565 + else {
  566 + x2 = m_energy[i].MeV();
  567 + f2 = this->rate(i, dir.detx(), dir.dety());
  568 + }
  569 +
  570 + // Compute integral
  571 + if (f1 > 0.0 && f2 > 0.0) {
  572 + rate += gammalib::plaw_integral(x1, f1, x2, f2);
  573 + }
  574 +
  575 + // Set second node as first node
  576 + x1 = x2;
  577 + f1 = f2;
  578 +
  579 + // If node energy is above emax then break now
  580 + if (m_energy[i] > emax) {
  581 + break;
  582 + }
  583 +
  584 + } // endfor: looped over all nodes
  585 +
  586 + // If last node energy is below emax then compute last part of
  587 + // integral up to emax
  588 + if (x1 < emax.MeV()) {
  589 + x2 = emax.MeV();
  590 + f2 = (*this)(emax.log10TeV(), dir.detx(), dir.dety());
  591 + if (f1 > 0.0 && f2 > 0.0) {
  592 + rate += gammalib::plaw_integral(x1, f1, x2, f2);
  593 + }
  594 + }
  595 +
  596 + } // endif: energy interval was positive
  597 +
  598 + // Return background rate
  599 + return rate;
  600 +}
  601 +
  602 +
  603 +/***********************************************************************//**
583 604 * @brief Print background information
584 605 *
585   - * @param[in] chatter Chattiness (defaults to NORMAL).
  606 + * @param[in] chatter Chattiness.
586 607 * @return String containing background information.
587 608 ***************************************************************************/
588 609 std::string GCTABackground3D::print(const GChatter& chatter) const
... ... @@ -661,6 +682,7 @@ void GCTABackground3D::init_members(void)
661 682 // Initialise members
662 683 m_filename.clear();
663 684 m_background.clear();
  685 + m_energy.clear();
664 686 m_inx_detx = 0;
665 687 m_inx_dety = 1;
666 688 m_inx_energy = 2;
... ... @@ -698,6 +720,7 @@ void GCTABackground3D::copy_members(const GCTABackground3D&amp; bgd)
698 720 // Copy members
699 721 m_filename = bgd.m_filename;
700 722 m_background = bgd.m_background;
  723 + m_energy = bgd.m_energy;
701 724 m_inx_detx = bgd.m_inx_detx;
702 725 m_inx_dety = bgd.m_inx_dety;
703 726 m_inx_energy = bgd.m_inx_energy;
... ... @@ -726,12 +749,59 @@ void GCTABackground3D::copy_members(const GCTABackground3D&amp; bgd)
726 749  
727 750  
728 751 /***********************************************************************//**
729   - * @brief Set DETX, DETY and log10(E/TeV) limits
  752 + * @brief Delete class members
  753 + ***************************************************************************/
  754 +void GCTABackground3D::free_members(void)
  755 +{
  756 + // Return
  757 + return;
  758 +}
  759 +
  760 +
  761 +/***********************************************************************//**
  762 + * @brief Set members from background table
730 763 *
731   - * Set the DETX, DETY and log10(E/TeV) limits.
  764 + * @exception GException::invalid_value
  765 + * Response table is not three-dimensional.
  766 + *
  767 + * Set class members based on the background table. The DETX and DETY axes
  768 + * will be set to radians, the energy axis will be set to log10.
732 769 ***************************************************************************/
733   -void GCTABackground3D::set_limits(void)
  770 +void GCTABackground3D::set_members(void)
734 771 {
  772 + // Throw an exception if the table is not three-dimensional
  773 + if (m_background.axes() != 3) {
  774 + std::string msg = "Expected three-dimensional background "
  775 + "response table but found "+
  776 + gammalib::str(m_background.axes())+
  777 + " dimensions. Please specify a three-dimensional "
  778 + "background.";
  779 + throw GException::invalid_value(G_SET_MEMBERS, msg);
  780 + }
  781 +
  782 + // Get mandatory indices (throw exception if not found)
  783 + m_inx_detx = m_background.axis("DETX");
  784 + m_inx_dety = m_background.axis("DETY");
  785 + m_inx_energy = m_background.axis("ENERG");
  786 + m_inx_bgd = m_background.table("BGD");
  787 +
  788 + // Get axes dimensions
  789 + m_num_detx = m_background.axis_bins(m_inx_detx);
  790 + m_num_dety = m_background.axis_bins(m_inx_dety);
  791 + m_num_energy = m_background.axis_bins(m_inx_energy);
  792 +
  793 + // Set dimension array
  794 + m_num[m_inx_detx] = m_num_detx;
  795 + m_num[m_inx_dety] = m_num_dety;
  796 + m_num[m_inx_energy] = m_num_energy;
  797 +
  798 + // Set DETX and DETY axis to radians
  799 + m_background.axis_radians(m_inx_detx);
  800 + m_background.axis_radians(m_inx_dety);
  801 +
  802 + // Set energy axis to logarithmic scale
  803 + m_background.axis_log10(m_inx_energy);
  804 +
735 805 // Compute DETX boundaries in radians
736 806 m_detx_min = m_background.axis_lo(m_inx_detx, 0) *
737 807 gammalib::deg2rad;
... ... @@ -752,6 +822,17 @@ void GCTABackground3D::set_limits(void)
752 822 m_logE_min = emin.log10TeV();
753 823 m_logE_max = emax.log10TeV();
754 824  
  825 + // Setup energy vector
  826 + m_energy.clear();
  827 + for (int i = 0; i < m_num_energy; ++i) {
  828 + GEnergy elo(m_background.axis_lo(m_inx_energy, i),
  829 + m_background.axis_lo_unit(m_inx_energy));
  830 + GEnergy ehi(m_background.axis_hi(m_inx_energy, i),
  831 + m_background.axis_hi_unit(m_inx_energy));
  832 + GEnergy energy(std::sqrt(elo.TeV()*ehi.TeV()), "TeV");
  833 + m_energy.append(energy);
  834 + }
  835 +
755 836 // Return
756 837 return;
757 838 }
... ... @@ -786,16 +867,6 @@ int GCTABackground3D::index(const int&amp; idetx,
786 867  
787 868  
788 869 /***********************************************************************//**
789   - * @brief Delete class members
790   - ***************************************************************************/
791   -void GCTABackground3D::free_members(void)
792   -{
793   - // Return
794   - return;
795   -}
796   -
797   -
798   -/***********************************************************************//**
799 870 * @brief Initialise Monte Carlo cache
800 871 *
801 872 * @exception GException::invalid_value
... ... @@ -923,7 +994,7 @@ void GCTABackground3D::init_mc_cache(void) const
923 994 if (m_mc_spectrum.nodes() == 0) {
924 995 std::string msg = "Background response table is empty. Please provide "
925 996 "a valid three-dimensional background template.";
926   - throw GException::invalid_value(G_READ, msg);
  997 + throw GException::invalid_value(G_INIT_MC_CACHE, msg);
927 998 }
928 999  
929 1000 // Return
... ... @@ -1045,3 +1116,47 @@ double GCTABackground3D::solid_angle(const double&amp; detx1, const double&amp; dety1,
1045 1116 // Return solid angle
1046 1117 return solidangle;
1047 1118 }
  1119 +
  1120 +
  1121 +/***********************************************************************//**
  1122 + * @brief Return background rate for a given energy bin and DETX-DETY value
  1123 + * (events/s/MeV/sr)
  1124 + *
  1125 + * @param[in] iebin Energy index.
  1126 + * @param[in] detx Tangential X coordinate in nominal system (radians).
  1127 + * @param[in] dety Tangential Y coordinate in nominal system (radians).
  1128 + * @return Background rate (events/s/MeV/sr)
  1129 + ***************************************************************************/
  1130 +double GCTABackground3D::rate(const int& iebin, const double& detx,
  1131 + const double& dety) const
  1132 +{
  1133 + // Retrieve references to node arrays
  1134 + const GNodeArray& detxs = m_background.axis_nodes(m_inx_detx);
  1135 + const GNodeArray& detys = m_background.axis_nodes(m_inx_dety);
  1136 +
  1137 + // Set DETX and DETY values for node array interpolation
  1138 + detxs.set_value(detx);
  1139 + detys.set_value(dety);
  1140 +
  1141 + // Get bin indices
  1142 + int inx_ll = index(detxs.inx_left(), detys.inx_left(), iebin);
  1143 + int inx_lr = index(detxs.inx_left(), detys.inx_right(), iebin);
  1144 + int inx_rl = index(detxs.inx_right(), detys.inx_left(), iebin);
  1145 + int inx_rr = index(detxs.inx_right(), detys.inx_right(), iebin);
  1146 +
  1147 + // Bi-linearly interpolate the rates
  1148 + double rate = detxs.wgt_left() *
  1149 + (m_background(m_inx_bgd, inx_ll) * detys.wgt_left() +
  1150 + m_background(m_inx_bgd, inx_lr) * detys.wgt_right()) +
  1151 + detxs.wgt_right() *
  1152 + (m_background(m_inx_bgd, inx_rl) * detys.wgt_left() +
  1153 + m_background(m_inx_bgd, inx_rr) * detys.wgt_right());
  1154 +
  1155 + // Make sure that background rate is not negative
  1156 + if (rate < 0.0) {
  1157 + rate = 0.0;
  1158 + }
  1159 +
  1160 + // Return
  1161 + return rate;
  1162 +}
... ...
inst/cta/src/GCTABackgroundPerfTable.cpp
  View file @5e4a21b
1 1 /***************************************************************************
2 2 * GCTABackgroundPerfTable.cpp - CTA performance table background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -33,6 +33,7 @@
33 33 #include "GMath.hpp"
34 34 #include "GIntegral.hpp"
35 35 #include "GRan.hpp"
  36 +#include "GCTAInstDir.hpp"
36 37 #include "GCTABackgroundPerfTable.hpp"
37 38  
38 39 /* __ Method name definitions ____________________________________________ */
... ... @@ -45,7 +46,7 @@
45 46  
46 47 /* __ Debug definitions __________________________________________________ */
47 48 //#define G_DEBUG_MC //!< Debug Monte Carlo method
48   -#define G_LOG_INTERPOLATION //!< Energy interpolate log10(background rate)
  49 +#define G_LOG_INTERPOLATION //!< Energy interpolate log(background rate)
49 50  
50 51 /* __ Constants __________________________________________________________ */
51 52  
... ... @@ -163,27 +164,23 @@ GCTABackgroundPerfTable&amp; GCTABackgroundPerfTable::operator=(const GCTABackground
163 164 * @brief Return background rate in units of events/s/MeV/sr
164 165 *
165 166 * @param[in] logE Log10 of the true photon energy (TeV).
166   - * @param[in] detx Tangential coord in nominal sys (rad).
167   - * @param[in] dety Tangential coord in nominal sys (rad).
  167 + * @param[in] detx Tangential coordinate in nominal system (rad).
  168 + * @param[in] dety Tangential coordinate in nominal system (rad).
  169 + * @return Background rate (events/s/MeV/sr)
168 170 *
169 171 * Returns the background rate in units of events/s/MeV/sr for a given energy
170 172 * and detector coordinates. The method assures that the background rate
171 173 * never becomes negative.
  174 + *
  175 + * If the performance table contains less than 2 nodes an empty value is
  176 + * returned.
172 177 ***************************************************************************/
173 178 double GCTABackgroundPerfTable::operator()(const double& logE,
174 179 const double& detx,
175 180 const double& dety) const
176 181 {
177   - // Get background rate
178   - double rate = m_logE.interpolate(logE, m_background);
179   - #if defined(G_LOG_INTERPOLATION)
180   - rate = std::pow(10.0, rate);
181   - #endif
182   -
183   - // Make sure that background rate is not negative
184   - if (rate < 0.0) {
185   - rate = 0.0;
186   - }
  182 + // Initialise rate
  183 + double rate = this->rate(logE);
187 184  
188 185 // Optionally add in Gaussian offset angle dependence
189 186 if (m_sigma != 0.0) {
... ... @@ -192,7 +189,7 @@ double GCTABackgroundPerfTable::operator()(const double&amp; logE,
192 189 double scale = std::exp(-0.5 * arg * arg);
193 190 rate *= scale;
194 191 }
195   -
  192 +
196 193 // Return background rate
197 194 return rate;
198 195 }
... ... @@ -248,8 +245,10 @@ GCTABackgroundPerfTable* GCTABackgroundPerfTable::clone(void) const
248 245 void GCTABackgroundPerfTable::load(const GFilename& filename)
249 246 {
250 247 // Clear arrays
251   - m_logE.clear();
  248 + m_energy.clear();
  249 + m_log10_energy.clear();
252 250 m_background.clear();
  251 + m_log_background.clear();
253 252  
254 253 // Allocate line buffer
255 254 const int n = 1000;
... ... @@ -297,13 +296,15 @@ void GCTABackgroundPerfTable::load(const GFilename&amp; filename)
297 296 bgrate = 0.0;
298 297 }
299 298  
300   - // Push elements in node array and vector
301   - m_logE.append(logE);
302   - #if defined(G_LOG_INTERPOLATION)
303   - m_background.push_back(std::log10(bgrate));
304   - #else
  299 + // Set energy
  300 + GEnergy energy;
  301 + energy.log10TeV(logE);
  302 +
  303 + // Push elements in node array and vectors
  304 + m_energy.append(energy);
  305 + m_log10_energy.append(logE);
305 306 m_background.push_back(bgrate);
306   - #endif
  307 + m_log_background.push_back(std::log(bgrate));
307 308  
308 309 } // endwhile: looped over lines
309 310  
... ... @@ -376,9 +377,96 @@ GCTAInstDir GCTABackgroundPerfTable::mc(const GEnergy&amp; energy,
376 377  
377 378  
378 379 /***********************************************************************//**
  380 + * @brief Returns background count rate integrated over energy interval
  381 + *
  382 + * @param[in] dir Instrument direction.
  383 + * @param[in] emin Minimum energy of energy interval.
  384 + * @param[in] emax Maximum energy of energy interval.
  385 + * @return Integrated background count rate (counts/src/sr).
  386 + *
  387 + * Returns the background count rate for a given instrument direction that
  388 + * is integrated over a specified energy interval.
  389 + *
  390 + * If the energy interval is not positive, a zero background rate is
  391 + * returned.
  392 + ***************************************************************************/
  393 +double GCTABackgroundPerfTable::rate_ebin(const GCTAInstDir& dir,
  394 + const GEnergy& emin,
  395 + const GEnergy& emax) const
  396 +{
  397 + // Initialise rate
  398 + double rate = 0.0;
  399 +
  400 + // Continue only if energy interval is positive
  401 + if (emax > emin) {
  402 +
  403 + // Initialise first and second node
  404 + double x1 = emin.MeV();
  405 + double x2 = 0.0;
  406 + double f1 = this->rate(emin.log10TeV());
  407 + double f2 = 0.0;
  408 +
  409 + // Loop over all nodes
  410 + for (int i = 0; i < size(); ++i) {
  411 +
  412 + // If node energy is below x1 then skip node
  413 + if (m_energy[i].MeV() <= x1) {
  414 + continue;
  415 + }
  416 +
  417 + // If node energy is above emax then use emax as energy
  418 + if (m_energy[i] > emax) {
  419 + x2 = emax.MeV();
  420 + f2 = this->rate(emax.log10TeV());
  421 + }
  422 +
  423 + // ... otherwise use node energy
  424 + else {
  425 + x2 = m_energy[i].MeV();
  426 + f2 = m_background[i];
  427 + }
  428 +
  429 + // Compute integral
  430 + rate += gammalib::plaw_integral(x1, f1, x2, f2);
  431 +
  432 + // Set second node as first node
  433 + x1 = x2;
  434 + f1 = f2;
  435 +
  436 + // If node energy is above emax then break now
  437 + if (m_energy[i] > emax) {
  438 + break;
  439 + }
  440 +
  441 + } // endfor: looped over all nodes
  442 +
  443 + // If last node energy is below emax then compute last part of
  444 + // integral up to emax
  445 + if (x1 < emax.MeV()) {
  446 + x2 = emax.MeV();
  447 + f2 = this->rate(emax.log10TeV());
  448 + rate += gammalib::plaw_integral(x1, f1, x2, f2);
  449 + }
  450 +
  451 + // Optionally add in Gaussian offset angle dependence
  452 + if (m_sigma != 0.0) {
  453 + double theta = dir.theta() * gammalib::rad2deg;
  454 + double arg = theta * theta / m_sigma;
  455 + double scale = std::exp(-0.5 * arg * arg);
  456 + rate *= scale;
  457 + }
  458 +
  459 + } // endif: energy interval was positive
  460 +
  461 + // Return background rate
  462 + return rate;
  463 +}
  464 +
  465 +
  466 +/***********************************************************************//**
379 467 * @brief Print background information
380 468 *
381   - * @param[in] chatter Chattiness (defaults to NORMAL).
  469 + * @param[in] chatter Chattiness.
382 470 * @return String containing background information.
383 471 ***************************************************************************/
384 472 std::string GCTABackgroundPerfTable::print(const GChatter& chatter) const
... ... @@ -390,8 +478,8 @@ std::string GCTABackgroundPerfTable::print(const GChatter&amp; chatter) const
390 478 if (chatter != SILENT) {
391 479  
392 480 // Compute energy boundaries in TeV
393   - double emin = std::pow(10.0, m_logE[0]);
394   - double emax = std::pow(10.0, m_logE[size()-1]);
  481 + double emin = (size() > 0) ? m_energy[0].TeV() : 0.0;
  482 + double emax = (size() > 1) ? m_energy[size()-1].TeV() : 0.0;
395 483  
396 484 // Append header
397 485 result.append("=== GCTABackgroundPerfTable ===");
... ... @@ -400,7 +488,7 @@ std::string GCTABackgroundPerfTable::print(const GChatter&amp; chatter) const
400 488 result.append("\n"+gammalib::parformat("Filename")+m_filename);
401 489 result.append("\n"+gammalib::parformat("Number of energy bins") +
402 490 gammalib::str(size()));
403   - result.append("\n"+gammalib::parformat("Log10(Energy) range"));
  491 + result.append("\n"+gammalib::parformat("Energy range"));
404 492 result.append(gammalib::str(emin) +
405 493 " - " +
406 494 gammalib::str(emax) +
... ... @@ -437,8 +525,10 @@ void GCTABackgroundPerfTable::init_members(void)
437 525 {
438 526 // Initialise members
439 527 m_filename.clear();
440   - m_logE.clear();
  528 + m_energy.clear();
  529 + m_log10_energy.clear();
441 530 m_background.clear();
  531 + m_log_background.clear();
442 532 m_sigma = 3.0;
443 533  
444 534 // Initialise MC cache
... ... @@ -457,10 +547,12 @@ void GCTABackgroundPerfTable::init_members(void)
457 547 void GCTABackgroundPerfTable::copy_members(const GCTABackgroundPerfTable& bgd)
458 548 {
459 549 // Copy members
460   - m_filename = bgd.m_filename;
461   - m_logE = bgd.m_logE;
462   - m_background = bgd.m_background;
463   - m_sigma = bgd.m_sigma;
  550 + m_filename = bgd.m_filename;
  551 + m_energy = bgd.m_energy;
  552 + m_log10_energy = bgd.m_log10_energy;
  553 + m_background = bgd.m_background;
  554 + m_log_background = bgd.m_log_background;
  555 + m_sigma = bgd.m_sigma;
464 556  
465 557 // Copy MC cache
466 558 m_mc_spectrum = bgd.m_mc_spectrum;
... ... @@ -539,20 +631,12 @@ void GCTABackgroundPerfTable::init_mc_cache(void) const
539 631 // Loop over nodes
540 632 for (int i = 0; i < size(); ++i) {
541 633  
542   - // Set energy
543   - GEnergy energy;
544   - energy.log10TeV(m_logE[i]);
545   -
546 634 // Compute total rate
547   - #if defined(G_LOG_INTERPOLATION)
548   - double total_rate = std::pow(10.0, m_background[i]) * solidangle;
549   - #else
550 635 double total_rate = m_background[i] * solidangle;
551   - #endif
552 636  
553 637 // Set node
554 638 if (total_rate > 0.0) {
555   - m_mc_spectrum.append(energy, total_rate);
  639 + m_mc_spectrum.append(m_energy[i], total_rate);
556 640 }
557 641  
558 642 }
... ... @@ -560,3 +644,36 @@ void GCTABackgroundPerfTable::init_mc_cache(void) const
560 644 // Return
561 645 return;
562 646 }
  647 +
  648 +
  649 +/***********************************************************************//**
  650 + * @brief Return background rate for a given energy (events/s/MeV/sr)
  651 + *
  652 + * @param[in] logE Log10 of the true photon energy (TeV).
  653 + * @return Background rate (events/s/MeV/sr)
  654 + ***************************************************************************/
  655 +double GCTABackgroundPerfTable::rate(const double& logE) const
  656 +{
  657 + // Initialise rate
  658 + double rate = 0.0;
  659 +
  660 + // Continue only if there are at least two nodes
  661 + if (size() > 1) {
  662 +
  663 + // Get background rate
  664 + #if defined(G_LOG_INTERPOLATION)
  665 + rate = std::exp(m_log10_energy.interpolate(logE, m_log_background));
  666 + #else
  667 + rate = m_log10_energy.interpolate(logE, m_background);
  668 + #endif
  669 +
  670 + // Make sure that background rate is not negative
  671 + if (rate < 0.0) {
  672 + rate = 0.0;
  673 + }
  674 +
  675 + } // endif: there were enough rates
  676 +
  677 + // Return
  678 + return rate;
  679 +}
... ...
inst/cta/test/test_CTA.cpp
  View file @5e4a21b
... ... @@ -64,9 +64,10 @@ const std::string cta_aeff_bgd_xml = datadir+&quot;/cta_model_aeff_bgd.xml&quot;;
64 64 const std::string cta_caldb_king = cta_caldb+"/data/cta/e/bcf/IFAE20120510_50h_King";
65 65 const std::string cta_irf_king = "irf_file.fits";
66 66 const std::string cta_psf_table = caldbdir+"/psf_table.fits[PSF_2D_TABLE]";
67   -const std::string cta_edisp_perf = caldbdir+"/cta_dummy_irf.dat";
  67 +const std::string cta_perf_table = caldbdir+"/cta_dummy_irf.dat";
68 68 const std::string cta_edisp_rmf = caldbdir+"/dc1/rmf.fits";
69 69 const std::string cta_edisp_2D = caldbdir+"/edisp_matrix.fits";
  70 +const std::string cta_bgd_3D = caldbdir+"/edisp_matrix.fits";
70 71 const std::string cta_modbck_fit = datadir+"/bg_test.fits";
71 72 const std::string cta_point_table = datadir+"/crab_pointing.fits";
72 73  
... ... @@ -84,6 +85,9 @@ const std::string cta_onoff_obs = datadir+&quot;/onoff_obs.xml&quot;;
84 85 const std::string cta_onoff_model = datadir+"/onoff_model.xml";
85 86 const std::string cta_onoff_onreg = datadir+"/onoff_region_on.reg";
86 87  
  88 +/* __ Debug definitions __________________________________________________ */
  89 +//#define G_COMPUTE_REF_RATE_EBIN //!< Compute reference rate
  90 +
87 91  
88 92 /***********************************************************************//**
89 93 * @brief Set miscellaneous CTA classes test methods
... ... @@ -129,26 +133,30 @@ void TestGCTAResponse::set(void)
129 133 "Test Table point spread function");
130 134 append(static_cast<pfunction>(&TestGCTAResponse::test_response_npsf),
131 135 "Test integrated point spread function");
132   - append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp),
133   - "Test energy dispersion");
134   - append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp_PerfTable),
135   - "Test Performance Table energy dispersion");
136   - append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp_RMF),
137   - "Test RMF energy dispersion");
138   - append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp_2D),
139   - "Test 2D energy dispersion");
140 136 append(static_cast<pfunction>(&TestGCTAResponse::test_response_irf_diffuse),
141 137 "Test diffuse IRF");
142 138 append(static_cast<pfunction>(&TestGCTAResponse::test_response_npred_diffuse),
143 139 "Test diffuse IRF integration");
  140 + append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp),
  141 + "Test GCTAEdisp class");
  142 + append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp_PerfTable),
  143 + "Test GCTAEdispPerfTable class");
  144 + append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp_RMF),
  145 + "Test GCTAEdispRmf class");
  146 + append(static_cast<pfunction>(&TestGCTAResponse::test_response_edisp_2D),
  147 + "Test GCTAEdisp2D class");
  148 + append(static_cast<pfunction>(&TestGCTAResponse::test_response_bgd_PerfTable),
  149 + "Test GCTABackgroundPerfTable class");
  150 + append(static_cast<pfunction>(&TestGCTAResponse::test_response_bgd_3D),
  151 + "Test GCTABackground3D class");
144 152 append(static_cast<pfunction>(&TestGCTAResponse::test_response_expcube),
145   - "Test exposure cube");
  153 + "Test GCTACubeExposure class");
146 154 append(static_cast<pfunction>(&TestGCTAResponse::test_response_psfcube),
147   - "Test point spread function cube");
  155 + "Test GCTACubePsf class");
148 156 append(static_cast<pfunction>(&TestGCTAResponse::test_response_edispcube),
149   - "Test energy dispersion cube");
  157 + "Test GCTACubeEdisp class");
150 158 append(static_cast<pfunction>(&TestGCTAResponse::test_response_bkgcube),
151   - "Test background cube");
  159 + "Test GCTACubeBackground class");
152 160  
153 161 // Return
154 162 return;
... ... @@ -608,6 +616,165 @@ void TestGCTAResponse::test_response_npsf(void)
608 616 return;
609 617 }
610 618  
  619 +
  620 +/***********************************************************************//**
  621 + * @brief Test CTA IRF computation for diffuse source model
  622 + *
  623 + * Tests the IRF computation for the diffuse source model. This is done
  624 + * by calling the GCTAObservation::model method which in turn calls the
  625 + * GCTAResponse::irf_diffuse method. The test is done for a small counts
  626 + * map to keep the test executing reasonably fast.
  627 + ***************************************************************************/
  628 +void TestGCTAResponse::test_response_irf_diffuse(void)
  629 +{
  630 + // Set reference value
  631 + //double ref = 13803.800313356;
  632 + //const double ref = 13803.5186374;
  633 + //const double ref = 13803.6932774; // after GWcs::solidangle improvement
  634 + //const double ref = 13803.3453994; // after Irf computation improvements
  635 + //const double ref = 13803.3461768; // using skymap flux() method
  636 + const double ref = 13803.5251338; // after correcting contains() method
  637 +
  638 + // Set parameters
  639 + double src_ra = 201.3651;
  640 + double src_dec = -43.0191;
  641 + int nebins = 5;
  642 +
  643 + // Setup pointing on Cen A
  644 + GSkyDir skyDir;
  645 + skyDir.radec_deg(src_ra, src_dec);
  646 + GCTAPointing pnt;
  647 + pnt.dir(skyDir);
  648 +
  649 + // Setup skymap (10 energy layers)
  650 + GSkyMap map("CAR", "CEL", src_ra, src_dec, 0.5, 0.5, 10, 10, nebins);
  651 +
  652 + // Setup time interval
  653 + GGti gti;
  654 + GTime tstart(0.0);
  655 + GTime tstop(1800.0);
  656 + gti.append(tstart, tstop);
  657 +
  658 + // Setup energy boundaries
  659 + GEnergy emin(0.1, "TeV");
  660 + GEnergy emax(100.0, "TeV");
  661 + GEbounds ebounds(nebins, emin, emax);
  662 +
  663 + // Setup event cube centreed on Cen A
  664 + GCTAEventCube cube(map, ebounds, gti);
  665 +
  666 + // Setup dummy CTA observation
  667 + GCTAObservation obs;
  668 + obs.ontime(1800.0);
  669 + obs.livetime(1600.0);
  670 + obs.deadc(1600.0/1800.0);
  671 + obs.response(cta_irf, GCaldb(cta_caldb));
  672 + obs.events(cube);
  673 + obs.pointing(pnt);
  674 +
  675 + // Load model for IRF computation
  676 + GModels models(cta_rsp_xml);
  677 +
  678 + // Reset sum
  679 + double sum = 0.0;
  680 +
  681 + // Iterate over all bins in event cube
  682 + for (int i = 0; i < obs.events()->size(); ++i) {
  683 +
  684 + // Get event pointer
  685 + const GEventBin* bin = (*(static_cast<const GEventCube*>(obs.events())))[i];
  686 +
  687 + // Get model and add to sum
  688 + double model = obs.model(models, *bin, NULL) * bin->size();
  689 + sum += model;
  690 +
  691 + }
  692 +
  693 + // Test sum
  694 + test_value(sum, ref, 0.2, "Diffuse IRF computation");
  695 +
  696 + // Return
  697 + return;
  698 +}
  699 +
  700 +
  701 +/***********************************************************************//**
  702 + * @brief Test CTA Npred computation
  703 + *
  704 + * Tests the Npred computation for the diffuse source model. This is done
  705 + * by loading the model from the XML file and by calling the
  706 + * GCTAObservation::npred method which in turn calls the
  707 + * GCTAResponse::npred_diffuse method. The test takes a few seconds.
  708 + ***************************************************************************/
  709 +void TestGCTAResponse::test_response_npred_diffuse(void)
  710 +{
  711 + // Set reference value
  712 + //const double ref = 11212.26274; // npred_spec precision of 1e-6
  713 + //const double ref = 11212.437464; // npred_spec precision of 1e-5
  714 + //const double ref = 11212.4370702; // After GWcs::solidangle improvement
  715 + //const double ref = 12644.3391902; // After correcting for deadtime bug
  716 + //const double ref = 12643.6159142; // npred_spec precision of 1e-6
  717 + //const double ref = 10997.0851335; // After response computation change
  718 + //const double ref = 11212.7161901; // npred_diffuse precision to iter=8
  719 + //const double ref = 11238.6928076; // npred_diffuse precision to iter=9
  720 + //const double ref = 11238.6933328; // using skymap flux() method
  721 + const double ref = 11238.8390288; // after correcting contains() method
  722 +
  723 + // Set parameters
  724 + double src_ra = 201.3651;
  725 + double src_dec = -43.0191;
  726 + double roi_rad = 4.0;
  727 +
  728 + // Setup ROI centred on Cen A with a radius of 4 deg
  729 + GCTARoi roi;
  730 + GCTAInstDir instDir;
  731 + instDir.dir().radec_deg(src_ra, src_dec);
  732 + roi.centre(instDir);
  733 + roi.radius(roi_rad);
  734 +
  735 + // Setup pointing on Cen A
  736 + GSkyDir skyDir;
  737 + skyDir.radec_deg(src_ra, src_dec);
  738 + GCTAPointing pnt;
  739 + pnt.dir(skyDir);
  740 +
  741 + // Setup dummy event list
  742 + GGti gti;
  743 + GEbounds ebounds;
  744 + GTime tstart(0.0);
  745 + GTime tstop(1800.0);
  746 + GEnergy emin(0.1, "TeV");
  747 + GEnergy emax(100.0, "TeV");
  748 + gti.append(tstart, tstop);
  749 + ebounds.append(emin, emax);
  750 + GCTAEventList events;
  751 + events.roi(roi);
  752 + events.gti(gti);
  753 + events.ebounds(ebounds);
  754 +
  755 + // Setup dummy CTA observation
  756 + GCTAObservation obs;
  757 + obs.ontime(1800.0);
  758 + obs.livetime(1600.0);
  759 + obs.deadc(1600.0/1800.0);
  760 + obs.response(cta_irf, GCaldb(cta_caldb));
  761 + obs.events(events);
  762 + obs.pointing(pnt);
  763 +
  764 + // Load models for Npred computation
  765 + GModels models(cta_rsp_xml);
  766 +
  767 + // Perform Npred computation
  768 + double npred = obs.npred(models, NULL);
  769 +
  770 + // Test Npred
  771 + test_value(npred, ref, 0.2, "Diffuse Npred computation");
  772 +
  773 + // Return
  774 + return;
  775 +}
  776 +
  777 +
611 778 /***********************************************************************//**
612 779 * @brief Test CTA Energy Dispersion computation
613 780 *
... ... @@ -641,7 +808,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
641 808 // Test GCTAEdispPerfTable constructor
642 809 test_try("GCTAEdispPerfTable constructor");
643 810 try {
644   - GCTAEdispPerfTable test(cta_edisp_perf);
  811 + GCTAEdispPerfTable test(cta_perf_table);
645 812 test_try_success();
646 813 }
647 814 catch (std::exception &e) {
... ... @@ -649,7 +816,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
649 816 }
650 817  
651 818 // Load energy dispersion
652   - GCTAEdispPerfTable edisp(cta_edisp_perf);
  819 + GCTAEdispPerfTable edisp(cta_perf_table);
653 820  
654 821 // Test ebounds_obs and ebounds_src methods
655 822 test_try("GCTAEdispPerfTable ebounds_obs and ebounds_src methods");
... ... @@ -687,8 +854,6 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
687 854 }
688 855  
689 856  
690   -
691   -
692 857 /***********************************************************************//**
693 858 * @brief Test CTA RMF energy dispersion class
694 859 ***************************************************************************/
... ... @@ -783,80 +948,99 @@ void TestGCTAResponse::test_response_edisp_2D(void)
783 948  
784 949  
785 950 /***********************************************************************//**
786   - * @brief Test CTA IRF computation for diffuse source model
787   - *
788   - * Tests the IRF computation for the diffuse source model. This is done
789   - * by calling the GCTAObservation::model method which in turn calls the
790   - * GCTAResponse::irf_diffuse method. The test is done for a small counts
791   - * map to keep the test executing reasonably fast.
  951 + * @brief Test GCTABackgroundPerfTable class
792 952 ***************************************************************************/
793   -void TestGCTAResponse::test_response_irf_diffuse(void)
  953 +void TestGCTAResponse::test_response_bgd_PerfTable(void)
794 954 {
795   - // Set reference value
796   - //double ref = 13803.800313356;
797   - //const double ref = 13803.5186374;
798   - //const double ref = 13803.6932774; // after GWcs::solidangle improvement
799   - //const double ref = 13803.3453994; // after Irf computation improvements
800   - //const double ref = 13803.3461768; // using skymap flux() method
801   - const double ref = 13803.5251338; // after correcting contains() method
802   -
803   - // Set parameters
804   - double src_ra = 201.3651;
805   - double src_dec = -43.0191;
806   - int nebins = 5;
807   -
808   - // Setup pointing on Cen A
809   - GSkyDir skyDir;
810   - skyDir.radec_deg(src_ra, src_dec);
811   - GCTAPointing pnt;
812   - pnt.dir(skyDir);
813   -
814   - // Setup skymap (10 energy layers)
815   - GSkyMap map("CAR", "CEL", src_ra, src_dec, 0.5, 0.5, 10, 10, nebins);
816   -
817   - // Setup time interval
818   - GGti gti;
819   - GTime tstart(0.0);
820   - GTime tstop(1800.0);
821   - gti.append(tstart, tstop);
822   -
823   - // Setup energy boundaries
824   - GEnergy emin(0.1, "TeV");
825   - GEnergy emax(100.0, "TeV");
826   - GEbounds ebounds(nebins, emin, emax);
827   -
828   - // Setup event cube centreed on Cen A
829   - GCTAEventCube cube(map, ebounds, gti);
830   -
831   - // Setup dummy CTA observation
832   - GCTAObservation obs;
833   - obs.ontime(1800.0);
834   - obs.livetime(1600.0);
835   - obs.deadc(1600.0/1800.0);
836   - obs.response(cta_irf, GCaldb(cta_caldb));
837   - obs.events(cube);
838   - obs.pointing(pnt);
839   -
840   - // Load model for IRF computation
841   - GModels models(cta_rsp_xml);
842   -
843   - // Reset sum
844   - double sum = 0.0;
845   -
846   - // Iterate over all bins in event cube
847   - for (int i = 0; i < obs.events()->size(); ++i) {
848   -
849   - // Get event pointer
850   - const GEventBin* bin = (*(static_cast<const GEventCube*>(obs.events())))[i];
851   -
852   - // Get model and add to sum
853   - double model = obs.model(models, *bin, NULL) * bin->size();
854   - sum += model;
855   -
  955 + // Set reference values
  956 + const double ref_rate_below = 0.000612295557409636;
  957 + const double ref_rate_above = 5.74072516267095e-11;
  958 + const double ref_rate_onaxis = 3.62782669235078e-05;
  959 + const double ref_rate_offaxis = 3.12350491143832e-05;
  960 + const double ref_rate_ebin_within = 44.9181299290927;
  961 + const double ref_rate_ebin_below = 3676.68195289697;
  962 + const double ref_rate_ebin_above = 76.2328310391366;
  963 +
  964 + // Set energies for energy range for rate integration
  965 + const GEnergy ebelow(10.0, "GeV");
  966 + const GEnergy emin(1.0, "TeV");
  967 + const GEnergy emax(10.0, "TeV");
  968 + const GEnergy eabove(200.0, "TeV");
  969 +
  970 + // Test GCTABackgroundPerfTable constructor
  971 + test_try("GCTABackgroundPerfTable constructor");
  972 + try {
  973 + GCTABackgroundPerfTable bgd1(cta_perf_table);
  974 + test_try_success();
  975 + }
  976 + catch (std::exception &e) {
  977 + test_try_failure(e);
856 978 }
857 979  
858   - // Test sum
859   - test_value(sum, ref, 0.2, "Diffuse IRF computation");
  980 + // Check content of empty instance
  981 + GCTABackgroundPerfTable bgd1;
  982 + test_value(bgd1.size(), 0, "Nodes of empty instance");
  983 + test_value(bgd1.sigma(), 3.0, "Sigma value of empty instance");
  984 + test_value(bgd1.filename(), "", "Filename of empty instance");
  985 + test_value(bgd1(0.0, 0.0, 0.0), 0.0,
  986 + "Onaxis operator value of empty instance");
  987 + test_value(bgd1(0.0, 0.01, 0.02), 0.0,
  988 + "Offaxis operator value of empty instance");
  989 + test_value(bgd1.rate_ebin(GCTAInstDir(), emin, emax), 0.0,
  990 + "Integrated rate value of empty instance");
  991 +
  992 + // Check content of filled instance
  993 + GCTABackgroundPerfTable bgd2(cta_perf_table);
  994 + test_value(bgd2.size(), 20, "Nodes of filled instance");
  995 + test_value(bgd2.sigma(), 3.0, "Sigma value of filled instance");
  996 + test_value(bgd2.filename(), cta_perf_table, "Filename of filled instance");
  997 + test_value(bgd2(-2.0, 0.0, 0.0), ref_rate_below,
  998 + "Operator value of filled instance extrapolated below lowest energy");
  999 + test_value(bgd2(3.0, 0.0, 0.0), ref_rate_above,
  1000 + "Operator value of filled instance extrapolated above highest energy");
  1001 + test_value(bgd2(0.0, 0.0, 0.0), ref_rate_onaxis,
  1002 + "Onaxis operator value of filled instance");
  1003 + test_value(bgd2(0.0, 0.01, 0.02), ref_rate_offaxis,
  1004 + "Offaxis operator value of filled instance");
  1005 + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within,
  1006 + "Integrated rate value of filled instance (within covered energies)");
  1007 + test_value(bgd2.rate_ebin(GCTAInstDir(), ebelow, emax), ref_rate_ebin_below,
  1008 + "Integrated rate value of filled instance (extrapolate below lowest energy)");
  1009 + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, eabove), ref_rate_ebin_above,
  1010 + "Integrated rate value of filled instance (extrapolate above highest energy)");
  1011 + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emin), 0.0,
  1012 + "Integrated rate value of filled instance for zero interval");
  1013 +
  1014 + // Check content of copied instance
  1015 + GCTABackgroundPerfTable bgd3(bgd2);
  1016 + test_value(bgd3.size(), 20, "Nodes of copied instance");
  1017 + test_value(bgd3.sigma(), 3.0, "Sigma value of copied instance");
  1018 + test_value(bgd3.filename(), cta_perf_table, "Filename of copied instance");
  1019 + test_value(bgd3(-2.0, 0.0, 0.0), ref_rate_below,
  1020 + "Operator value of copied instance extrapolated below lowest energy");
  1021 + test_value(bgd3(3.0, 0.0, 0.0), ref_rate_above,
  1022 + "Operator value of copied instance extrapolated above highest energy");
  1023 + test_value(bgd3(0.0, 0.0, 0.0), ref_rate_onaxis,
  1024 + "Onaxis operator value of copied instance");
  1025 + test_value(bgd3(0.0, 0.01, 0.02), ref_rate_offaxis,
  1026 + "Offaxis operator value of copied instance");
  1027 + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within,
  1028 + "Integrated rate value of copied instance");
  1029 + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emin), 0.0,
  1030 + "Integrated rate value of copied instance for zero interval");
  1031 +
  1032 + // Compute reference rate
  1033 + #if defined(G_COMPUTE_REF_RATE_EBIN)
  1034 + double rate = 0.0;
  1035 + int nbins = 10000;
  1036 + GEnergy E = emin;
  1037 + GEnergy dE = (eabove-emin)/nbins;
  1038 + for (int i = 0; i < nbins; ++i) {
  1039 + rate += bgd2(E.log10TeV(), 0.0, 0.0) * dE.MeV();
  1040 + E += dE;
  1041 + }
  1042 + std::cout << rate << std::endl;
  1043 + #endif
860 1044  
861 1045 // Return
862 1046 return;
... ... @@ -864,76 +1048,96 @@ void TestGCTAResponse::test_response_irf_diffuse(void)
864 1048  
865 1049  
866 1050 /***********************************************************************//**
867   - * @brief Test CTA Npred computation
868   - *
869   - * Tests the Npred computation for the diffuse source model. This is done
870   - * by loading the model from the XML file and by calling the
871   - * GCTAObservation::npred method which in turn calls the
872   - * GCTAResponse::npred_diffuse method. The test takes a few seconds.
  1051 + * @brief Test GCTABackground3D class
873 1052 ***************************************************************************/
874   -void TestGCTAResponse::test_response_npred_diffuse(void)
  1053 +void TestGCTAResponse::test_response_bgd_3D(void)
875 1054 {
876   - // Set reference value
877   - //const double ref = 11212.26274; // npred_spec precision of 1e-6
878   - //const double ref = 11212.437464; // npred_spec precision of 1e-5
879   - //const double ref = 11212.4370702; // After GWcs::solidangle improvement
880   - //const double ref = 12644.3391902; // After correcting for deadtime bug
881   - //const double ref = 12643.6159142; // npred_spec precision of 1e-6
882   - //const double ref = 10997.0851335; // After response computation change
883   - //const double ref = 11212.7161901; // npred_diffuse precision to iter=8
884   - //const double ref = 11238.6928076; // npred_diffuse precision to iter=9
885   - //const double ref = 11238.6933328; // using skymap flux() method
886   - const double ref = 11238.8390288; // after correcting contains() method
887   -
888   - // Set parameters
889   - double src_ra = 201.3651;
890   - double src_dec = -43.0191;
891   - double roi_rad = 4.0;
892   -
893   - // Setup ROI centred on Cen A with a radius of 4 deg
894   - GCTARoi roi;
895   - GCTAInstDir instDir;
896   - instDir.dir().radec_deg(src_ra, src_dec);
897   - roi.centre(instDir);
898   - roi.radius(roi_rad);
899   -
900   - // Setup pointing on Cen A
901   - GSkyDir skyDir;
902   - skyDir.radec_deg(src_ra, src_dec);
903   - GCTAPointing pnt;
904   - pnt.dir(skyDir);
905   -
906   - // Setup dummy event list
907   - GGti gti;
908   - GEbounds ebounds;
909   - GTime tstart(0.0);
910   - GTime tstop(1800.0);
911   - GEnergy emin(0.1, "TeV");
912   - GEnergy emax(100.0, "TeV");
913   - gti.append(tstart, tstop);
914   - ebounds.append(emin, emax);
915   - GCTAEventList events;
916   - events.roi(roi);
917   - events.gti(gti);
918   - events.ebounds(ebounds);
919   -
920   - // Setup dummy CTA observation
921   - GCTAObservation obs;
922   - obs.ontime(1800.0);
923   - obs.livetime(1600.0);
924   - obs.deadc(1600.0/1800.0);
925   - obs.response(cta_irf, GCaldb(cta_caldb));
926   - obs.events(events);
927   - obs.pointing(pnt);
928   -
929   - // Load models for Npred computation
930   - GModels models(cta_rsp_xml);
  1055 + // Set reference values
  1056 + const double ref_rate_below = 0.0;
  1057 + const double ref_rate_above = 0.0;
  1058 + const double ref_rate_onaxis = 9.7632195934345e-05;
  1059 + const double ref_rate_offaxis = 13.5162501643753e-05;
  1060 + const double ref_rate_ebin_within = 65.6069433285784;
  1061 + const double ref_rate_ebin_below = 4647.37614735018;
  1062 + const double ref_rate_ebin_above = 69.5331125123672;
  1063 +
  1064 + // Set energies for energy range for rate integration
  1065 + const GEnergy ebelow(10.0, "GeV");
  1066 + const GEnergy emin(1.0, "TeV");
  1067 + const GEnergy emax(10.0, "TeV");
  1068 + const GEnergy eabove(200.0, "TeV");
  1069 +
  1070 + // Test GCTABackground3D constructor
  1071 + test_try("GCTABackground3D constructor");
  1072 + try {
  1073 + GCTABackground3D bgd1(cta_bgd_3D);
  1074 + test_try_success();
  1075 + }
  1076 + catch (std::exception &e) {
  1077 + test_try_failure(e);
  1078 + }
931 1079  
932   - // Perform Npred computation
933   - double npred = obs.npred(models, NULL);
  1080 + // Check content of empty instance
  1081 + GCTABackground3D bgd1;
  1082 + test_value(bgd1.table().axes(), 0, "Response table axes of empty instance");
  1083 + test_value(bgd1.filename(), "", "Filename of empty instance");
  1084 + test_value(bgd1(0.0, 0.0, 0.0), 0.0,
  1085 + "Onaxis operator value of empty instance");
  1086 + test_value(bgd1(0.0, 0.01, 0.02), 0.0,
  1087 + "Offaxis operator value of empty instance");
  1088 + test_value(bgd1.rate_ebin(GCTAInstDir(), emin, emax), 0.0,
  1089 + "Integrated rate value of empty instance");
934 1090  
935   - // Test Npred
936   - test_value(npred, ref, 0.2, "Diffuse Npred computation");
  1091 + // Check content of filled instance
  1092 + GCTABackground3D bgd2(cta_bgd_3D);
  1093 + test_value(bgd2.table().axes(), 3, "Response table axes of filled instance");
  1094 + test_value(bgd2.filename(), cta_bgd_3D, "Filename of filled instance");
  1095 + test_value(bgd2(-2.0, 0.0, 0.0), ref_rate_below,
  1096 + "Operator value of filled instance extrapolated below lowest energy");
  1097 + test_value(bgd2(3.0, 0.0, 0.0), ref_rate_above,
  1098 + "Operator value of filled instance extrapolated above highest energy");
  1099 + test_value(bgd2(0.0, 0.0, 0.0), ref_rate_onaxis,
  1100 + "Onaxis operator value of filled instance");
  1101 + test_value(bgd2(0.0, 0.01, 0.02), ref_rate_offaxis,
  1102 + "Offaxis operator value of filled instance");
  1103 + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within,
  1104 + "Integrated rate value of filled instance");
  1105 + test_value(bgd2.rate_ebin(GCTAInstDir(), ebelow, emax), ref_rate_ebin_below,
  1106 + "Integrated rate value of filled instance (extrapolate below lowest energy)");
  1107 + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, eabove), ref_rate_ebin_above,
  1108 + "Integrated rate value of filled instance (extrapolate above highest energy)");
  1109 + test_value(bgd2.rate_ebin(GCTAInstDir(), emin, emin), 0.0,
  1110 + "Integrated rate value of filled instance for zero interval");
  1111 +
  1112 + // Check content of copied instance
  1113 + GCTABackground3D bgd3(bgd2);
  1114 + test_value(bgd3.table().axes(), 3, "Response table axes of copied instance");
  1115 + test_value(bgd3.filename(), cta_bgd_3D, "Filename of copied instance");
  1116 + test_value(bgd3(-2.0, 0.0, 0.0), ref_rate_below,
  1117 + "Operator value of copied instance extrapolated below lowest energy");
  1118 + test_value(bgd3(3.0, 0.0, 0.0), ref_rate_above,
  1119 + "Operator value of copied instance extrapolated above highest energy");
  1120 + test_value(bgd3(0.0, 0.0, 0.0), ref_rate_onaxis,
  1121 + "Onaxis operator value of copied instance");
  1122 + test_value(bgd3(0.0, 0.01, 0.02), ref_rate_offaxis,
  1123 + "Offaxis operator value of copied instance");
  1124 + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emax), ref_rate_ebin_within,
  1125 + "Integrated rate value of copied instance");
  1126 + test_value(bgd3.rate_ebin(GCTAInstDir(), emin, emin), 0.0,
  1127 + "Integrated rate value of copied instance for zero interval");
  1128 +
  1129 + // Compute reference rate
  1130 + #if defined(G_COMPUTE_REF_RATE_EBIN)
  1131 + double rate = 0.0;
  1132 + int nbins = 10000;
  1133 + GEnergy E = emin;
  1134 + GEnergy dE = (eabove-E)/nbins;
  1135 + for (int i = 0; i < nbins; ++i) {
  1136 + rate += bgd2(E.log10TeV(), 0.0, 0.0) * dE.MeV();
  1137 + E += dE;
  1138 + }
  1139 + std::cout << rate << std::endl;
  1140 + #endif
937 1141  
938 1142 // Return
939 1143 return;
... ...
inst/cta/test/test_CTA.hpp
  View file @5e4a21b
... ... @@ -78,12 +78,14 @@ public:
78 78 void test_response_psf_king(void);
79 79 void test_response_psf_table(void);
80 80 void test_response_npsf(void);
  81 + void test_response_irf_diffuse(void);
  82 + void test_response_npred_diffuse(void);
81 83 void test_response_edisp(void);
82 84 void test_response_edisp_PerfTable(void);
83 85 void test_response_edisp_RMF(void);
84 86 void test_response_edisp_2D(void);
85   - void test_response_irf_diffuse(void);
86   - void test_response_npred_diffuse(void);
  87 + void test_response_bgd_PerfTable(void);
  88 + void test_response_bgd_3D(void);
87 89 void test_response_expcube(void);
88 90 void test_response_psfcube(void);
89 91 void test_response_bkgcube(void);
... ...