GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAPsfPerfTable.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAPsfPerfTable.cpp - CTA performance table PSF 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 GCTAPsfPerfTable.hpp
23  * @brief CTA performance table point spread function class implementation
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 <cmath>
33 #include "GException.hpp"
34 #include "GTools.hpp"
35 #include "GMath.hpp"
36 #include "GRan.hpp"
37 #include "GCTAPsfPerfTable.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_LOAD "GCTAPsfPerfTable::load(GFilename&)"
41 #define G_CONTAINMENT_RADIUS "GCTAPsfPerfTable::containment_radius(double&,"\
42  " double&, double&, double&, double&, double&, bool&)"
43 
44 /* __ Macros _____________________________________________________________ */
45 
46 /* __ Coding definitions _________________________________________________ */
47 #define G_SMOOTH_PSF
48 
49 /* __ Debug definitions __________________________________________________ */
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 Performance table file name.
77  *
78  * Construct instance by loading the point spread function information from
79  * an ASCII performance table.
80  ***************************************************************************/
82 {
83  // Initialise class members
84  init_members();
85 
86  // Load point spread function from file
87  load(filename);
88 
89  // Return
90  return;
91 }
92 
93 
94 /***********************************************************************//**
95  * @brief Copy constructor
96  *
97  * @param[in] psf Point spread function.
98  ***************************************************************************/
100 {
101  // Initialise class members
102  init_members();
103 
104  // Copy members
105  copy_members(psf);
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] psf Point spread function.
135  * @return Point spread function.
136  ***************************************************************************/
138 {
139  // Execute only if object is not identical
140  if (this != &psf) {
141 
142  // Copy base class members
143  this->GCTAPsf::operator=(psf);
144 
145  // Free members
146  free_members();
147 
148  // Initialise private members
149  init_members();
150 
151  // Copy members
152  copy_members(psf);
153 
154  } // endif: object was not identical
155 
156  // Return this object
157  return *this;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief Return point spread function (in units of sr^-1)
163  *
164  * @param[in] delta Angular separation between true and measured photon
165  * directions (rad).
166  * @param[in] logE Log10 of the true photon energy (TeV).
167  * @param[in] theta Offset angle in camera system (rad). Not used.
168  * @param[in] phi Azimuth angle in camera system (rad). Not used.
169  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
170  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
171  * @param[in] etrue Use true energy (true/false). Not used.
172  *
173  * Returns the point spread function for a given angular separation in units
174  * of sr^-1 for a given energy.
175  ***************************************************************************/
176 double GCTAPsfPerfTable::operator()(const double& delta,
177  const double& logE,
178  const double& theta,
179  const double& phi,
180  const double& zenith,
181  const double& azimuth,
182  const bool& etrue) const
183 {
184  #if defined(G_SMOOTH_PSF)
185  // Compute offset so that PSF goes to 0 at 5 times the sigma value. This
186  // is a kluge to get a PSF that smoothly goes to zero at the edge, which
187  // prevents steps or kinks in the log-likelihood function.
188  static const double offset = std::exp(-0.5*5.0*5.0);
189  #endif
190 
191  // Update the parameter cache
192  update(logE);
193 
194  // Compute PSF value
195  #if defined(G_SMOOTH_PSF)
196  double psf = m_par_scale * (std::exp(m_par_width * delta * delta) - offset);
197  #else
198  double psf = m_par_scale * (std::exp(m_par_width * delta * delta));
199  #endif
200 
201  #if defined(G_SMOOTH_PSF)
202  // Make sure that PSF is non-negative
203  if (psf < 0.0) {
204  psf = 0.0;
205  }
206  #endif
207 
208  // Return PSF
209  return psf;
210 }
211 
212 
213 /*==========================================================================
214  = =
215  = Public methods =
216  = =
217  ==========================================================================*/
218 
219 /***********************************************************************//**
220  * @brief Clear instance
221  *
222  * This method properly resets the object to an initial state.
223  ***************************************************************************/
225 {
226  // Free class members (base and derived classes, derived class first)
227  free_members();
228  this->GCTAPsf::free_members();
229 
230  // Initialise members
231  this->GCTAPsf::init_members();
232  init_members();
233 
234  // Return
235  return;
236 }
237 
238 
239 /***********************************************************************//**
240  * @brief Clone instance
241  *
242  * @return Deep copy of point spread function instance.
243  ***************************************************************************/
245 {
246  return new GCTAPsfPerfTable(*this);
247 }
248 
249 
250 /***********************************************************************//**
251  * @brief Load point spread function from performance table
252  *
253  * @param[in] filename Performance table file name.
254  *
255  * @exception GException::file_error
256  * File could not be opened for read access.
257  *
258  * This method loads the point spread function information from an ASCII
259  * performance table.
260  ***************************************************************************/
261 void GCTAPsfPerfTable::load(const GFilename& filename)
262 {
263  // Set conversion factor from 68% containment radius to 1 sigma
264  const double conv = 0.6624305 * gammalib::deg2rad;
265 
266  // Clear arrays
267  m_logE.clear();
268  m_r68.clear();
269  m_r80.clear();
270  m_sigma.clear();
271 
272  // Allocate line buffer
273  const int n = 1000;
274  char line[n];
275 
276  // Open performance table readonly
277  FILE* fptr = std::fopen(filename.url().c_str(), "r");
278  if (fptr == NULL) {
279  std::string msg = "Point spread function file \""+filename.url()+
280  "\" not found or readable. Please specify a valid "
281  "and readable point spread function file.";
282  throw GException::file_error(G_LOAD, msg);
283  }
284 
285  // Read lines
286  while (std::fgets(line, n, fptr) != NULL) {
287 
288  // Split line in elements. Strip empty elements from vector.
289  std::vector<std::string> elements = gammalib::split(line, " ");
290  for (int i = elements.size()-1; i >= 0; i--) {
291  if (gammalib::strip_whitespace(elements[i]).length() == 0) {
292  elements.erase(elements.begin()+i);
293  }
294  }
295 
296  // Skip header
297  if (elements[0].find("log(E)") != std::string::npos) {
298  continue;
299  }
300 
301  // Break loop if end of data table has been reached
302  if (elements[0].find("----------") != std::string::npos) {
303  break;
304  }
305 
306  // Push elements in node array and vector
307  m_logE.append(gammalib::todouble(elements[0]));
308  m_r68.push_back(gammalib::todouble(elements[2]));
309  m_r80.push_back(gammalib::todouble(elements[3]));
310  m_sigma.push_back(gammalib::todouble(elements[2])*conv);
311 
312  } // endwhile: looped over lines
313 
314  // Close file
315  std::fclose(fptr);
316 
317  // Store filename
319 
320  // Return
321  return;
322 }
323 
324 
325 /***********************************************************************//**
326  * @brief Simulate PSF offset (radians)
327  *
328  * @param[in] ran Random number generator.
329  * @param[in] logE Log10 of the true photon energy (TeV).
330  * @param[in] theta Offset angle in camera system (rad). Not used.
331  * @param[in] phi Azimuth angle in camera system (rad). Not used.
332  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
333  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
334  * @param[in] etrue Use true energy (true/false). Not used.
335  ***************************************************************************/
337  const double& logE,
338  const double& theta,
339  const double& phi,
340  const double& zenith,
341  const double& azimuth,
342  const bool& etrue) const
343 {
344  // Update the parameter cache
345  update(logE);
346 
347  // Draw offset
348  double delta = m_par_sigma * ran.chisq2();
349 
350  // Return PSF offset
351  return delta;
352 }
353 
354 
355 /***********************************************************************//**
356  * @brief Return maximum size of PSF (radians)
357  *
358  * @param[in] logE Log10 of the true photon energy (TeV).
359  * @param[in] theta Offset angle in camera system (rad). Not used.
360  * @param[in] phi Azimuth angle in camera system (rad). Not used.
361  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
362  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
363  * @param[in] etrue Use true energy (true/false). Not used.
364  *
365  * Determine the radius beyond which the PSF becomes negligible. This radius
366  * is set by this method to \f$5 \times \sigma\f$, where \f$\sigma\f$ is the
367  * Gaussian width of the largest PSF component.
368  ***************************************************************************/
369 double GCTAPsfPerfTable::delta_max(const double& logE,
370  const double& theta,
371  const double& phi,
372  const double& zenith,
373  const double& azimuth,
374  const bool& etrue) const
375 {
376  // Update the parameter cache
377  update(logE);
378 
379  // Compute maximum PSF radius
380  double radius = 5.0 * m_par_sigma;
381 
382  // Return maximum PSF radius
383  return radius;
384 }
385 
386 
387 /***********************************************************************//**
388  * @brief Return the radius that contains a fraction of the events (radians)
389  *
390  * @param[in] fraction of events (0.0-1.0)
391  * @param[in] logE Log10 of the true photon energy (TeV).
392  * @param[in] theta Offset angle in camera system (rad). Not used.
393  * @param[in] phi Azimuth angle in camera system (rad). Not used.
394  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
395  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
396  * @param[in] etrue Use true energy (true/false). Not used.
397  * @return Containment radius (radians).
398  *
399  * @exception GException::invalid_argument
400  * Invalid fraction specified.
401  *
402  * Evaluates:
403  *
404  * \f[
405  *
406  * radius = \sqrt{ \frac{\ln{\left( \frac{ fraction * m\_par\_width}
407  * {\pi*m\_par\_scale} +1\right)}}{m\_par\_width} }
408  *
409  * \f]
410  *
411  * which is derived by integrating
412  *
413  * \f[
414  *
415  * fraction = \int_{0}^{2\pi}\int_{0}^{radius}r *
416  * e^{ m\_par\_width * r^{2}}dr d\phi
417  *
418  * \f]
419  *
420  * and solving for radius.
421  *
422  * Calculate the radius from the center that contains 'fraction' percent
423  * of the events. fraction * 100. = Containment % .
424  ***************************************************************************/
425 double GCTAPsfPerfTable::containment_radius(const double& fraction,
426  const double& logE,
427  const double& theta,
428  const double& phi,
429  const double& zenith,
430  const double& azimuth,
431  const bool& etrue) const
432 {
433  // Check input argument
434  if (fraction <= 0.0 || fraction >= 1.0) {
435  std::string message = "Containment fraction "+
436  gammalib::str(fraction)+" must be between " +
437  "0.0 and 1.0, not inclusive.";
439  }
440 
441  // Update the parameter cache
442  update(logE);
443 
444  // Compute radius
445  double arg = fraction * m_par_width / ( gammalib::pi * m_par_scale);
446  double radius = std::sqrt(std::log(arg + 1.0) / m_par_width);
447 
448  // Return containement radius
449  return radius;
450 }
451 
452 
453 /***********************************************************************//**
454  * @brief Print point spread function information
455  *
456  * @param[in] chatter Chattiness.
457  * @return String containing point spread function information.
458  ***************************************************************************/
459 std::string GCTAPsfPerfTable::print(const GChatter& chatter) const
460 {
461  // Initialise result string
462  std::string result;
463 
464  // Continue only if chatter is not silent
465  if (chatter != SILENT) {
466 
467  // Compute energy boundaries in TeV
468  int num = m_logE.size();
469  double emin = std::pow(10.0, m_logE[0]);
470  double emax = std::pow(10.0, m_logE[num-1]);
471 
472  // Append header
473  result.append("=== GCTAPsfPerfTable ===");
474 
475  // Append information
476  result.append("\n"+gammalib::parformat("Filename")+m_filename);
477  result.append("\n"+gammalib::parformat("Number of energy bins") +
478  gammalib::str(num));
479  result.append("\n"+gammalib::parformat("Energy range"));
480  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
481 
482  } // endif: chatter was not silent
483 
484  // Return result
485  return result;
486 }
487 
488 
489 /*==========================================================================
490  = =
491  = Private methods =
492  = =
493  ==========================================================================*/
494 
495 /***********************************************************************//**
496  * @brief Initialise class members
497  ***************************************************************************/
499 {
500  // Initialise members
501  m_filename.clear();
502  m_logE.clear();
503  m_r68.clear();
504  m_r80.clear();
505  m_sigma.clear();
506  m_par_logE = -1.0e30;
507  m_par_scale = 1.0;
508  m_par_sigma = 0.0;
509  m_par_width = 0.0;
510 
511  // Return
512  return;
513 }
514 
515 
516 /***********************************************************************//**
517  * @brief Copy class members
518  *
519  * @param[in] psf Point spread function.
520  ***************************************************************************/
522 {
523  // Copy members
524  m_filename = psf.m_filename;
525  m_logE = psf.m_logE;
526  m_r68 = psf.m_r68;
527  m_r80 = psf.m_r80;
528  m_sigma = psf.m_sigma;
529  m_par_logE = psf.m_par_logE;
530  m_par_scale = psf.m_par_scale;
531  m_par_sigma = psf.m_par_sigma;
532  m_par_width = psf.m_par_width;
533 
534  // Return
535  return;
536 }
537 
538 
539 /***********************************************************************//**
540  * @brief Delete class members
541  ***************************************************************************/
543 {
544  // Return
545  return;
546 }
547 
548 
549 /***********************************************************************//**
550  * @brief Update PSF parameter cache
551  *
552  * @param[in] logE Log10 of the true photon energy (TeV).
553  *
554  * This method updates the PSF parameter cache. As the performance table PSF
555  * only depends on energy, the only parameter on which the cache values
556  * depend is the energy.
557  ***************************************************************************/
558 void GCTAPsfPerfTable::update(const double& logE) const
559 {
560  // Only compute PSF parameters if arguments have changed
561  if (logE != m_par_logE) {
562 
563  // Save energy
564  m_par_logE = logE;
565 
566  // Determine Gaussian sigma in radians
568 
569  // Derive width=-0.5/(sigma*sigma) and scale=1/(twopi*sigma*sigma)
570  double sigma2 = m_par_sigma * m_par_sigma;
571  m_par_scale = 1.0 / (gammalib::twopi * sigma2);
572  m_par_width = -0.5 / sigma2;
573 
574  }
575 
576  // Return
577  return;
578 }
void free_members(void)
Delete class members.
Definition: GCTAPsf.cpp:164
double containment_radius(const double &fraction, 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 the radius that contains a fraction of the events (radians)
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:192
double chisq2(void)
Returns Chi2 deviates for 2 degrees of freedom.
Definition: GRan.cpp:370
void update(const double &logE) const
Update PSF parameter cache.
Random number generator class definition.
virtual ~GCTAPsfPerfTable(void)
Destructor.
const double pi
Definition: GMath.hpp:35
double mc(GRan &ran, 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
Simulate PSF offset (radians)
double m_par_scale
Gaussian normalization.
GCTAPsfPerfTable * clone(void) const
Clone instance.
Abstract base class for the CTA point spread function.
Definition: GCTAPsf.hpp:47
void load(const GFilename &filename)
Load point spread function from performance table.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
Random number generator class.
Definition: GRan.hpp:44
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:983
Gammalib tools definition.
GFilename m_filename
Name of Aeff response file.
GCTAPsf & operator=(const GCTAPsf &psf)
Assignment operator.
Definition: GCTAPsf.cpp:106
double m_par_sigma
Gaussian sigma (radians)
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:80
void init_members(void)
Initialise class members.
Definition: GCTAPsf.cpp:142
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function information.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
#define G_CONTAINMENT_RADIUS
void free_members(void)
Delete class members.
Filename class.
Definition: GFilename.hpp:62
void clear(void)
Clear instance.
GNodeArray m_logE
log(E) nodes for Aeff interpolation
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
GChatter
Definition: GTypemaps.hpp:33
void copy_members(const GCTAPsfPerfTable &psf)
Copy class members.
std::vector< double > m_r80
80% containment radius of PSF in degrees
CTA performance table point spread function class definition.
std::vector< double > m_r68
68% containment radius of PSF in degrees
#define G_LOAD
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void init_members(void)
Initialise class members.
GFilename filename(void) const
Return filename.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
double operator()(const double &delta, 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 point spread function (in units of sr^-1)
GCTAPsfPerfTable & operator=(const GCTAPsfPerfTable &psf)
Assignment operator.
const double twopi
Definition: GMath.hpp:36
CTA performance table point spread function class.
double delta_max(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 maximum size of PSF (radians)
std::vector< double > m_sigma
Sigma value of PSF in radians.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
GCTAPsfPerfTable(void)
Void constructor.
double m_par_width
Gaussian width parameter.
Mathematical function definitions.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:926
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
double m_par_logE
Energy for which precomputation is done.