GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAAeffPerfTable.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAAeffPerfTable.hpp - CTA performance table 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 GCTAAeffPerfTable.hpp
23  * @brief CTA performance table effective area class definition
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdio> // std::fopen, std::fgets, and std::fclose
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GFitsTable.hpp"
35 #include "GCTAAeffPerfTable.hpp"
36 #include "GCTAException.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_LOAD "GCTAAeffPerfTable::load(std::string&)"
40 
41 /* __ Macros _____________________________________________________________ */
42 
43 /* __ Coding definitions _________________________________________________ */
44 
45 /* __ Debug definitions __________________________________________________ */
46 
47 /* __ Constants __________________________________________________________ */
48 
49 
50 /*==========================================================================
51  = =
52  = Constructors/destructors =
53  = =
54  ==========================================================================*/
55 
56 /***********************************************************************//**
57  * @brief Void constructor
58  ***************************************************************************/
60 {
61  // Initialise class members
62  init_members();
63 
64  // Return
65  return;
66 }
67 
68 
69 /***********************************************************************//**
70  * @brief File constructor
71  *
72  * @param[in] filename Performance table file name.
73  *
74  * Construct instance by loading the effective area information from an
75  * ASCII performance table.
76  ***************************************************************************/
78 {
79  // Initialise class members
80  init_members();
81 
82  // Load effective area from file
83  load(filename);
84 
85  // Return
86  return;
87 }
88 
89 
90 /***********************************************************************//**
91  * @brief Copy constructor
92  *
93  * @param[in] aeff Effective area.
94  ***************************************************************************/
96 {
97  // Initialise class members
98  init_members();
99 
100  // Copy members
101  copy_members(aeff);
102 
103  // Return
104  return;
105 }
106 
107 
108 /***********************************************************************//**
109  * @brief Destructor
110  ***************************************************************************/
112 {
113  // Free members
114  free_members();
115 
116  // Return
117  return;
118 }
119 
120 
121 /*==========================================================================
122  = =
123  = Operators =
124  = =
125  ==========================================================================*/
126 
127 /***********************************************************************//**
128  * @brief Assignment operator
129  *
130  * @param[in] aeff Effective area.
131  * @return Effective area.
132  ***************************************************************************/
134 {
135  // Execute only if object is not identical
136  if (this != &aeff) {
137 
138  // Copy base class members
139  this->GCTAAeff::operator=(aeff);
140 
141  // Free members
142  free_members();
143 
144  // Initialise private members
145  init_members();
146 
147  // Copy members
148  copy_members(aeff);
149 
150  } // endif: object was not identical
151 
152  // Return this object
153  return *this;
154 }
155 
156 
157 /***********************************************************************//**
158  * @brief Return effective area in units of cm2
159  *
160  * @param[in] logE Log10 of the true photon energy (TeV).
161  * @param[in] theta Offset angle in camera system (rad).
162  * @param[in] phi Azimuth angle in camera system (rad). Not used.
163  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
164  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
165  * @param[in] etrue Use true energy (true/false). Not used.
166  *
167  * Returns the effective area in units of cm2 for a given energy and
168  * offset angle. The effective area is linearily interpolated in
169  * log10(energy). The method assures that the effective area value never
170  * becomes negative.
171  *
172  * Outside the energy range that is covered by the performance table the
173  * effective area will be set to zero.
174  ***************************************************************************/
175 double GCTAAeffPerfTable::operator()(const double& logE,
176  const double& theta,
177  const double& phi,
178  const double& zenith,
179  const double& azimuth,
180  const bool& etrue) const
181 {
182  // Initialise effective area
183  double aeff = 0.0;
184 
185  // Continue only if logE is in validity range
186  if ((logE >= m_logE_min) && (logE <= m_logE_max)) {
187 
188  // Get effective area value in cm2
189  aeff = m_logE.interpolate(logE, m_aeff);
190 
191  // Make sure that effective area is not negative
192  if (aeff < 0.0) {
193  aeff = 0.0;
194  }
195 
196  // Optionally add in Gaussian offset angle dependence
197  if (m_sigma != 0.0) {
198  double offset = theta * gammalib::rad2deg;
199  double arg = offset * offset / m_sigma;
200  double scale = std::exp(-0.5 * arg * arg);
201  aeff *= scale;
202  }
203 
204  } // endif: logE in validity range
205 
206  // Return effective area value
207  return aeff;
208 }
209 
210 
211 /*==========================================================================
212  = =
213  = Public methods =
214  = =
215  ==========================================================================*/
216 
217 /***********************************************************************//**
218  * @brief Clear instance
219  *
220  * This method properly resets the object to an initial state.
221  ***************************************************************************/
223 {
224  // Free class members (base and derived classes, derived class first)
225  free_members();
226  this->GCTAAeff::free_members();
227 
228  // Initialise members
229  this->GCTAAeff::init_members();
230  init_members();
231 
232  // Return
233  return;
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Clone instance
239  *
240  * @return Deep copy of effective area instance.
241  ***************************************************************************/
243 {
244  return new GCTAAeffPerfTable(*this);
245 }
246 
247 
248 /***********************************************************************//**
249  * @brief Load effective area from performance table
250  *
251  * @param[in] filename Performance table file name.
252  *
253  * @exception GCTAExceptionHandler::file_open_error
254  * File could not be opened for read access.
255  *
256  * This method loads the effective area information from an ASCII
257  * performance table.
258  ***************************************************************************/
259 void GCTAAeffPerfTable::load(const GFilename& filename)
260 {
261  // Clear arrays
262  m_logE.clear();
263  m_aeff.clear();
264 
265  // Allocate line buffer
266  const int n = 1000;
267  char line[n];
268 
269  // Open performance table readonly
270  FILE* fptr = std::fopen(filename.url().c_str(), "r");
271  if (fptr == NULL) {
272  throw GCTAException::file_open_error(G_LOAD, filename.url());
273  }
274 
275  // Read lines
276  while (std::fgets(line, n, fptr) != NULL) {
277 
278  // Split line in elements. Strip empty elements from vector.
279  std::vector<std::string> elements = gammalib::split(line, " ");
280  for (int i = elements.size()-1; i >= 0; i--) {
281  if (gammalib::strip_whitespace(elements[i]).length() == 0) {
282  elements.erase(elements.begin()+i);
283  }
284  }
285 
286  // Skip header
287  if (elements[0].find("log(E)") != std::string::npos) {
288  continue;
289  }
290 
291  // Break loop if end of data table has been reached
292  if (elements[0].find("----------") != std::string::npos) {
293  break;
294  }
295 
296  // Push elements in node array and vector
297  m_logE.append(gammalib::todouble(elements[0]));
298  m_aeff.push_back(gammalib::todouble(elements[1])*10000.0);
299 
300  } // endwhile: looped over lines
301 
302  // Close file
303  std::fclose(fptr);
304 
305  // Store filename
307 
308  // Set the energy boundaries
309  set_boundaries();
310 
311  // Return
312  return;
313 }
314 
315 
316 /***********************************************************************//**
317  * @brief Return maximum effective area at a given energy
318  *
319  * @param[in] logE Log10 of the true photon energy (TeV).
320  * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
321  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
322  * @param[in] etrue Use true energy (true/false). Not used.
323  * @return Maximum effective area (cm2).
324  ***************************************************************************/
325 double GCTAAeffPerfTable::max(const double& logE,
326  const double& zenith,
327  const double& azimuth,
328  const bool& etrue) const
329 {
330  // Get effective area value in cm2
331  double aeff_max = m_logE.interpolate(logE, m_aeff);
332 
333  // Make sure that effective area is not negative
334  if (aeff_max < 0.0) {
335  aeff_max = 0.0;
336  }
337 
338  // Return result
339  return aeff_max;
340 }
341 
342 
343 /***********************************************************************//**
344  * @brief Print effective area information
345  *
346  * @param[in] chatter Chattiness.
347  * @return String containing effective area information.
348  ***************************************************************************/
349 std::string GCTAAeffPerfTable::print(const GChatter& chatter) const
350 {
351  // Initialise result string
352  std::string result;
353 
354  // Continue only if chatter is not silent
355  if (chatter != SILENT) {
356 
357  // Append header
358  result.append("=== GCTAAeffPerfTable ===");
359 
360  // Append information
361  result.append("\n"+gammalib::parformat("Filename")+m_filename);
362  result.append("\n"+gammalib::parformat("Number of energy bins") +
363  gammalib::str(size()));
364  result.append("\n"+gammalib::parformat("Energy range"));
365  result.append(m_ebounds.emin().print() + " - " +
366  m_ebounds.emax().print());
367 
368  // Append offset angle dependence
369  if (m_sigma == 0) {
370  result.append("\n"+gammalib::parformat("Offset angle dependence") +
371  "none");
372  }
373  else {
374  std::string txt = "Fixed sigma=" + gammalib::str(m_sigma);
375  result.append("\n"+gammalib::parformat("Offset angle dependence") +
376  txt);
377  }
378 
379  } // endif: chatter was not silent
380 
381  // Return result
382  return result;
383 }
384 
385 
386 /*==========================================================================
387  = =
388  = Private methods =
389  = =
390  ==========================================================================*/
391 
392 /***********************************************************************//**
393  * @brief Initialise class members
394  ***************************************************************************/
396 {
397  // Initialise members
398  m_filename.clear();
399  m_logE.clear();
400  m_aeff.clear();
401  m_ebounds.clear();
402  m_sigma = 3.0;
403  m_logE_min = 0.0;
404  m_logE_max = 0.0;
405 
406  // Return
407  return;
408 }
409 
410 
411 /***********************************************************************//**
412  * @brief Copy class members
413  *
414  * @param[in] aeff Effective area.
415  ***************************************************************************/
417 {
418  // Copy members
419  m_filename = aeff.m_filename;
420  m_logE = aeff.m_logE;
421  m_aeff = aeff.m_aeff;
422  m_ebounds = aeff.m_ebounds;
423  m_sigma = aeff.m_sigma;
424  m_logE_min = aeff.m_logE_min;
425  m_logE_max = aeff.m_logE_max;
426 
427  // Return
428  return;
429 }
430 
431 
432 /***********************************************************************//**
433  * @brief Delete class members
434  ***************************************************************************/
436 {
437  // Return
438  return;
439 }
440 
441 
442 /***********************************************************************//**
443  * @brief Set effective area boundaries
444  *
445  * Sets the data members m_ebounds, m_logE_min and m_logE_max that define
446  * the validity range of the effective area.
447  ***************************************************************************/
449 {
450  // Clear energy boundaries
451  m_ebounds.clear();
452 
453  // Set log10 of minimum and maximum energies. Since the energy values are
454  // given at the bin centre we subtract half of the distance to the second
455  // bin from the minimum energy and we add half of the distance to the before
456  // last bin to the maximum energy
457  m_logE_min = m_logE[0] - 0.5*(m_logE[1] - m_logE[0]);
458  m_logE_max = m_logE[m_logE.size()-1] + 0.5*(m_logE[m_logE.size()-1] -
459  m_logE[m_logE.size()-2]);
460 
461  // Set energy boundaries
462  GEnergy emin;
463  GEnergy emax;
464  emin.log10TeV(m_logE_min);
465  emax.log10TeV(m_logE_max);
466  m_ebounds.append(emin, emax);
467 
468  // Return
469  return;
470 }
void load(const GFilename &filename)
Load effective area from performance table.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
GCTAAeffPerfTable & operator=(const GCTAAeffPerfTable &aeff)
Assignment operator.
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.
CTA performance table effective area class.
GEbounds m_ebounds
Energy boundaries.
double m_logE_min
Minimum logE (log10(E/TeV))
void set_boundaries(void)
Set effective area boundaries.
void free_members(void)
Delete class members.
std::string print(const GChatter &chatter=NORMAL) const
Print effective area information.
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:302
GNodeArray m_logE
log(E) nodes for Aeff interpolation
double m_logE_max
Maximum logE (log10(E/TeV))
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
virtual ~GCTAAeffPerfTable(void)
Destructor.
GCTAAeffPerfTable * clone(void) const
Clone instance.
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:862
Gammalib tools definition.
GCTAAeff & operator=(const GCTAAeff &aeff)
Assignment operator.
Definition: GCTAAeff.cpp:106
void copy_members(const GCTAAeffPerfTable &aeff)
Copy class members.
double max(const double &logE, const double &zenith, const double &azimuth, const bool &etrue=true) const
Return maximum effective area at a given energy.
GFilename filename(void) const
Return filename.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
void free_members(void)
Delete class members.
Definition: GCTAAeff.cpp:164
#define G_LOAD
CTA performance table effective area class definition.
double m_sigma
Sigma for offset angle computation (0=none)
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
GFilename m_filename
Name of Aeff response file.
Filename class.
Definition: GFilename.hpp:62
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
void clear(void)
Clear instance.
void init_members(void)
Initialise class members.
Definition: GCTAAeff.cpp:142
void init_members(void)
Initialise class members.
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
GChatter
Definition: GTypemaps.hpp:33
GCTAAeffPerfTable(void)
Void constructor.
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:269
CTA exception handler interface definition.
Abstract base class for the CTA effective area.
Definition: GCTAAeff.hpp:54
std::vector< double > m_aeff
Effective area in cm2.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
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
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
int size(void) const
Return number of node energies in response.
Mathematical function definitions.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:805
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.