GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAAeffArf.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAAeffArf.cpp - CTA ARF effective area class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2017 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCTAAeffArf.hpp
23  * @brief CTA ARF effective area class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GMath.hpp"
33 #include "GIntegral.hpp"
34 #include "GFilename.hpp"
35 #include "GFitsTable.hpp"
36 #include "GFitsTableCol.hpp"
37 #include "GArf.hpp"
38 #include "GCTAAeffArf.hpp"
39 #include "GCTAResponseIrf.hpp"
40 #include "GCTAResponse_helpers.hpp"
41 #include "GCTAException.hpp"
42 
43 /* __ Method name definitions ____________________________________________ */
44 
45 /* __ Macros _____________________________________________________________ */
46 
47 /* __ Coding definitions _________________________________________________ */
48 
49 /* __ Debug definitions __________________________________________________ */
50 //#define G_DEBUG_APPLY_THETACUT //!< Debug thetacut application
51 
52 /* __ Constants __________________________________________________________ */
53 
54 
55 /*==========================================================================
56  = =
57  = Constructors/destructors =
58  = =
59  ==========================================================================*/
60 
61 /***********************************************************************//**
62  * @brief Void constructor
63  ***************************************************************************/
65 {
66  // Initialise class members
67  init_members();
68 
69  // Return
70  return;
71 }
72 
73 
74 /***********************************************************************//**
75  * @brief File constructor
76  *
77  * @param[in] filename ARF FITS file name.
78  *
79  * Constructs effective area from an ARF FITS file.
80  ***************************************************************************/
82 {
83  // Initialise class members
84  init_members();
85 
86  // Load ARF
87  load(filename);
88 
89  // Return
90  return;
91 }
92 
93 
94 /***********************************************************************//**
95  * @brief Copy constructor
96  *
97  * @param[in] aeff Effective area.
98  ***************************************************************************/
100 {
101  // Initialise class members
102  init_members();
103 
104  // Copy members
105  copy_members(aeff);
106 
107  // Return
108  return;
109 }
110 
111 
112 /***********************************************************************//**
113  * @brief Destructor
114  ***************************************************************************/
116 {
117  // Free members
118  free_members();
119 
120  // Return
121  return;
122 }
123 
124 
125 /*==========================================================================
126  = =
127  = Operators =
128  = =
129  ==========================================================================*/
130 
131 /***********************************************************************//**
132  * @brief Assignment operator
133  *
134  * @param[in] aeff Effective area.
135  * @return Effective area.
136  ***************************************************************************/
138 {
139  // Execute only if object is not identical
140  if (this != &aeff) {
141 
142  // Copy base class members
143  this->GCTAAeff::operator=(aeff);
144 
145  // Free members
146  free_members();
147 
148  // Initialise private members
149  init_members();
150 
151  // Copy members
152  copy_members(aeff);
153 
154  } // endif: object was not identical
155 
156  // Return this object
157  return *this;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief Return effective area in units of cm2
163  *
164  * @param[in] logE Log10 of the true photon energy (TeV).
165  * @param[in] theta Offset angle in camera system (rad).
166  * @param[in] phi Azimuth angle in camera system (rad). Not used.
167  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
168  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
169  * @param[in] etrue Use true energy (true/false). Not used.
170  * @return Effective area in cm2.
171  *
172  * Returns the effective area in units of cm2 for a given energy and
173  * offset angle. The effective area is linearily interpolated in
174  * log10(energy). The method assures that the effective area value never
175  * becomes negative.
176  *
177  * Outside the energy range that is covered by the ARF vector the effective
178  * area will be set to zero.
179  ***************************************************************************/
180 double GCTAAeffArf::operator()(const double& logE,
181  const double& theta,
182  const double& phi,
183  const double& zenith,
184  const double& azimuth,
185  const bool& etrue) const
186 {
187  // Initialise effective area
188  double aeff = 0.0;
189 
190  // Continue only if logE is in validity range
191  if ((logE >= m_logE_min) && (logE <= m_logE_max)) {
192 
193  // Get effective area value in cm2
194  aeff = m_logE.interpolate(logE, m_aeff);
195 
196  // Make sure that effective area is not negative
197  if (aeff < 0.0) {
198  aeff = 0.0;
199  }
200 
201  // Optionally add in Gaussian offset angle dependence
202  if (m_sigma != 0.0) {
203  double offset = theta * gammalib::rad2deg;
204  double arg = offset * offset / m_sigma;
205  double scale = exp(-0.5 * arg * arg);
206  aeff *= scale;
207  }
208 
209  } // endif: logE in validity range
210 
211  // Return effective area value
212  return aeff;
213 }
214 
215 
216 /*==========================================================================
217  = =
218  = Public methods =
219  = =
220  ==========================================================================*/
221 
222 /***********************************************************************//**
223  * @brief Clear instance
224  *
225  * This method properly resets the object to an initial state.
226  ***************************************************************************/
228 {
229  // Free class members (base and derived classes, derived class first)
230  free_members();
231  this->GCTAAeff::free_members();
232 
233  // Initialise members
234  this->GCTAAeff::init_members();
235  init_members();
236 
237  // Return
238  return;
239 }
240 
241 
242 /***********************************************************************//**
243  * @brief Clone instance
244  *
245  * @return Deep copy of effective area instance.
246  ***************************************************************************/
248 {
249  return new GCTAAeffArf(*this);
250 }
251 
252 
253 /***********************************************************************//**
254  * @brief Load effective area from ARF FITS file
255  *
256  * @param[in] filename ARF FITS file name.
257  *
258  * Loads the effective area from an ARF FITS file.
259  *
260  * If no extension name is provided, the effective area will be loaded from
261  * the `SPECRESP` extension.
262  ***************************************************************************/
263 void GCTAAeffArf::load(const GFilename& filename)
264 {
265  // Open FITS file
266  GFits fits(filename);
267 
268  // Get ARF table
269  const GFitsTable& table = *fits.table(filename.extname(gammalib::extname_arf));
270 
271  // Read ARF from table
272  read(table);
273 
274  // Close FITS file
275  fits.close();
276 
277  // Store filename
279 
280  // Return
281  return;
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Read CTA ARF vector
287  *
288  * @param[in] table FITS table.
289  *
290  * Reads a CTA ARF vector from the FITS @p table.
291  *
292  * The energies are converted to TeV and the effective area is converted to
293  * cm2. Conversion is done based on the units provided for the energy and
294  * effective area columns. Units that are recognized are 'keV', 'MeV', 'GeV',
295  * 'TeV', 'm^2', 'm2', 'cm^2' and 'cm^2' (case independent).
296  *
297  * @todo Assign appropriate theta angle for PSF. So far we use onaxis.
298  * For appropriate theta angle assignment, we would need this
299  * information in the response header.
300  ***************************************************************************/
301 void GCTAAeffArf::read(const GFitsTable& table)
302 {
303  // Clear arrays
304  m_logE.clear();
305  m_aeff.clear();
306  m_ebounds.clear();
307 
308  // Get pointers to table columns
309  const GFitsTableCol* energy_lo = table["ENERG_LO"];
310  const GFitsTableCol* energy_hi = table["ENERG_HI"];
311  const GFitsTableCol* specresp = table["SPECRESP"];
312 
313  // Determine unit conversion factors (default: TeV and cm^2)
314  std::string u_specresp = gammalib::tolower(
315  gammalib::strip_whitespace(specresp->unit()));
316  double c_specresp = 1.0;
317  if (u_specresp == "m^2" || u_specresp == "m2") {
318  c_specresp = 10000.0;
319  }
320 
321  // Extract number of energy bins
322  int num = energy_lo->nrows();
323 
324  // Set nodes
325  for (int i = 0; i < num; ++i) {
326 
327  // Compute energy boundaries
328  GEnergy emin(energy_lo->real(i), energy_lo->unit());
329  GEnergy emax(energy_hi->real(i), energy_hi->unit());
330 
331  // Compute log10 of mean energy in TeV
332  double logE = 0.5 * (emin.log10TeV() + emax.log10TeV());
333 
334  // Initialise scale factor
335  double scale = m_scale;
336 
337  // Compute effective area in cm2
338  double aeff = specresp->real(i) * c_specresp * scale;
339 
340  // Store log10 mean energy and effective area value
341  m_logE.append(logE);
342  m_aeff.push_back(aeff);
343 
344  } // endfor: looped over nodes
345 
346  // Set energy boundaries
347  GEnergy emin(energy_lo->real(0), energy_lo->unit());
348  GEnergy emax(energy_hi->real(num-1), energy_hi->unit());
349  m_logE_min = emin.log10TeV();
350  m_logE_max = emax.log10TeV();
351  m_ebounds.append(emin, emax);
352 
353  // Disable offset angle dependence
354  m_sigma = 0.0;
355 
356  // Return
357  return;
358 }
359 
360 
361 /***********************************************************************//**
362  * @brief Remove thetacut
363  *
364  * @param[in] rsp CTA response.
365  *
366  * Removes thetacut from Aeff values read from a FITS file. Note that this
367  * method should only be called once directly after loading all response
368  * components.
369  ***************************************************************************/
371 {
372  // Set number of iterations for Romberg integration.
373  static const int iter = 6;
374 
375  // Continue only if thetacut value has been set
376  if (m_thetacut > 0.0) {
377 
378  // Get maximum integration radius
379  double rmax = m_thetacut * gammalib::deg2rad;
380 
381  // Loop over Aeff array
382  for (int i = 0; i < m_aeff.size(); ++i) {
383 
384  // Setup integration kernel for on-axis PSF
385  cta_npsf_kern_rad_azsym integrand(&rsp,
386  rmax,
387  0.0,
388  m_logE[i],
389  0.0,
390  0.0,
391  0.0,
392  0.0);
393 
394  // Setup integration
395  GIntegral integral(&integrand);
396 
397  // Set fixed number of iterations
398  integral.fixed_iter(iter);
399 
400  // Perform integration
401  double fraction = integral.romberg(0.0, rmax);
402 
403  // Set scale factor
404  double scale = 1.0;
405  if (fraction > 0.0) {
406  scale /= fraction;
407  #if defined(G_DEBUG_APPLY_THETACUT)
408  std::cout << "GCTAAeffArf::apply_thetacut:";
409  std::cout << " logE=" << m_logE[i];
410  std::cout << " scale=" << scale;
411  std::cout << " fraction=" << fraction;
412  std::cout << std::endl;
413  #endif
414  }
415  else {
416  std::cout << "WARNING: GCTAAeffArf::apply_thetacut:";
417  std::cout << " Non-positive integral occured in";
418  std::cout << " PSF integration.";
419  std::cout << std::endl;
420  }
421 
422  // Apply scaling factor
423  if (scale != 1.0) {
424  m_aeff[i] *= scale;
425  }
426 
427  } // endfor: looped over Aeff array
428 
429  } // endif: thetacut value was set
430 
431  // Return
432  return;
433 }
434 
435 
436 /***********************************************************************//**
437  * @brief Return maximum effective area at a given energy
438  *
439  * @param[in] logE Log10 of the true photon energy (TeV).
440  * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
441  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
442  * @param[in] etrue Use true energy (true/false). Not used.
443  * @return Maximum effective area (cm2).
444  ***************************************************************************/
445 double GCTAAeffArf::max(const double& logE,
446  const double& zenith,
447  const double& azimuth,
448  const bool& etrue) const
449 {
450  // Get effective area value in cm2
451  double aeff_max = m_logE.interpolate(logE, m_aeff);
452 
453  // Make sure that effective area is not negative
454  if (aeff_max < 0.0) {
455  aeff_max = 0.0;
456  }
457 
458  // Return result
459  return aeff_max;
460 }
461 
462 
463 /***********************************************************************//**
464  * @brief Print effective area information
465  *
466  * @param[in] chatter Chattiness.
467  * @return String containing effective area information.
468  ***************************************************************************/
469 std::string GCTAAeffArf::print(const GChatter& chatter) const
470 {
471  // Initialise result string
472  std::string result;
473 
474  // Continue only if chatter is not silent
475  if (chatter != SILENT) {
476 
477  // Append header
478  result.append("=== GCTAAeffArf ===");
479 
480  // Append information
481  result.append("\n"+gammalib::parformat("Filename")+m_filename);
482  result.append("\n"+gammalib::parformat("Number of energy bins") +
483  gammalib::str(size()));
484  result.append("\n"+gammalib::parformat("Energy range"));
485  result.append(m_ebounds.emin().print() + " - " +
486  m_ebounds.emax().print());
487 
488  // Append offset angle dependence
489  if (m_sigma == 0) {
490  result.append("\n"+gammalib::parformat("Offset angle dependence") +
491  "none");
492  }
493  else {
494  std::string txt = "Fixed sigma=" + gammalib::str(m_sigma);
495  result.append("\n"+gammalib::parformat("Offset angle dependence") +
496  txt);
497  }
498 
499  } // endif: chatter was not silent
500 
501  // Return result
502  return result;
503 }
504 
505 
506 /*==========================================================================
507  = =
508  = Private methods =
509  = =
510  ==========================================================================*/
511 
512 /***********************************************************************//**
513  * @brief Initialise class members
514  ***************************************************************************/
516 {
517  // Initialise members
518  m_filename.clear();
519  m_logE.clear();
520  m_aeff.clear();
521  m_ebounds.clear();
522  m_sigma = 0.0;
523  m_thetacut = 0.0;
524  m_scale = 1.0;
525  m_logE_min = 0.0;
526  m_logE_max = 0.0;
527 
528  // Return
529  return;
530 }
531 
532 
533 /***********************************************************************//**
534  * @brief Copy class members
535  *
536  * @param[in] aeff Effective area.
537  ***************************************************************************/
539 {
540  // Copy members
541  m_filename = aeff.m_filename;
542  m_logE = aeff.m_logE;
543  m_aeff = aeff.m_aeff;
544  m_ebounds = aeff.m_ebounds;
545  m_sigma = aeff.m_sigma;
546  m_thetacut = aeff.m_thetacut;
547  m_scale = aeff.m_scale;
548  m_logE_min = aeff.m_logE_min;
549  m_logE_max = aeff.m_logE_max;
550 
551  // Return
552  return;
553 }
554 
555 
556 /***********************************************************************//**
557  * @brief Delete class members
558  ***************************************************************************/
560 {
561  // Return
562  return;
563 }
void unit(const std::string &unit)
Set column unit.
GFilename filename(void) const
Return filename.
CTA response helper classes definition.
void remove_thetacut(const GCTAResponseIrf &rsp)
Remove thetacut.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
GCTAAeffArf & operator=(const GCTAAeffArf &aeff)
Assignment operator.
const double & scale(void) const
Return effective area scaling factor.
GNodeArray m_logE
log(E) nodes for Aeff interpolation
GEbounds m_ebounds
Energy boundaries.
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:302
GCTAAeffArf * clone(void) const
Clone instance.
std::vector< double > m_aeff
Effective area in cm2.
GFilename m_filename
Name of Aeff response file.
CTA instrument response function class definition.
void nrows(const int &nrows)
Set number of rows in column.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
void clear(void)
Clear instance.
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
CTA ARF effective area class definition.
double max(const double &logE, const double &zenith, const double &azimuth, const bool &etrue=true) const
Return maximum effective area at a given energy.
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
FITS file class.
Definition: GFits.hpp:63
int size(void) const
Return number of node energies in response.
GCTAAeff & operator=(const GCTAAeff &aeff)
Assignment operator.
Definition: GCTAAeff.cpp:106
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
virtual ~GCTAAeffArf(void)
Destructor.
FITS table column abstract base class definition.
void free_members(void)
Delete class members.
Definition: GCTAAeff.cpp:164
double operator()(const double &logE, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0, const bool &etrue=true) const
Return effective area in units of cm2.
GCTAAeffArf(void)
Void constructor.
Definition: GCTAAeffArf.cpp:64
const double deg2rad
Definition: GMath.hpp:43
const std::string extname_arf
Definition: GArf.hpp:45
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
Filename class.
Definition: GFilename.hpp:62
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
Abstract interface for FITS table column.
XSPEC Auxiliary Response File class definition.
void init_members(void)
Initialise class members.
Definition: GCTAAeff.cpp:142
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
double m_thetacut
Theta cut for ARF.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void free_members(void)
Delete class members.
GChatter
Definition: GTypemaps.hpp:33
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:269
CTA exception handler interface definition.
void copy_members(const GCTAAeffArf &aeff)
Copy class members.
double m_sigma
Sigma for offset angle computation (0=none)
CTA instrument response function class.
Abstract base class for the CTA effective area.
Definition: GCTAAeff.hpp:54
void read(const GFitsTable &table)
Read CTA ARF vector.
void load(const GFilename &filename)
Load effective area from ARF FITS file.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
virtual double real(const int &row, const int &inx=0) const =0
double m_scale
Scale for ARF.
double m_logE_max
Maximum logE (log10(E/TeV))
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:834
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
void init_members(void)
Initialise class members.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
double m_logE_min
Minimum logE (log10(E/TeV))
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
Integration kernel for npsf() method.
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
Integration class interface definition.
std::string print(const GChatter &chatter=NORMAL) const
Print effective area information.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
Filename class interface definition.
Mathematical function definitions.
CTA ARF effective area class.
Definition: GCTAAeffArf.hpp:51
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.