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