GammaLib  1.7.0.dev
 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-2016 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 "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GRan.hpp"
36 #include "GCTAPsfPerfTable.hpp"
37 #include "GCTAException.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 GCTAExceptionHandler::file_open_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  throw GCTAException::file_open_error(G_LOAD, filename);
280  }
281 
282  // Read lines
283  while (std::fgets(line, n, fptr) != NULL) {
284 
285  // Split line in elements. Strip empty elements from vector.
286  std::vector<std::string> elements = gammalib::split(line, " ");
287  for (int i = elements.size()-1; i >= 0; i--) {
288  if (gammalib::strip_whitespace(elements[i]).length() == 0) {
289  elements.erase(elements.begin()+i);
290  }
291  }
292 
293  // Skip header
294  if (elements[0].find("log(E)") != std::string::npos) {
295  continue;
296  }
297 
298  // Break loop if end of data table has been reached
299  if (elements[0].find("----------") != std::string::npos) {
300  break;
301  }
302 
303  // Push elements in node array and vector
304  m_logE.append(gammalib::todouble(elements[0]));
305  m_r68.push_back(gammalib::todouble(elements[2]));
306  m_r80.push_back(gammalib::todouble(elements[3]));
307  m_sigma.push_back(gammalib::todouble(elements[2])*conv);
308 
309  } // endwhile: looped over lines
310 
311  // Close file
312  std::fclose(fptr);
313 
314  // Store filename
316 
317  // Return
318  return;
319 }
320 
321 
322 /***********************************************************************//**
323  * @brief Simulate PSF offset (radians)
324  *
325  * @param[in] ran Random number generator.
326  * @param[in] logE Log10 of the true photon energy (TeV).
327  * @param[in] theta Offset angle in camera system (rad). Not used.
328  * @param[in] phi Azimuth angle in camera system (rad). Not used.
329  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
330  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
331  * @param[in] etrue Use true energy (true/false). Not used.
332  ***************************************************************************/
334  const double& logE,
335  const double& theta,
336  const double& phi,
337  const double& zenith,
338  const double& azimuth,
339  const bool& etrue) const
340 {
341  // Update the parameter cache
342  update(logE);
343 
344  // Draw offset
345  double delta = m_par_sigma * ran.chisq2();
346 
347  // Return PSF offset
348  return delta;
349 }
350 
351 
352 /***********************************************************************//**
353  * @brief Return maximum size of PSF (radians)
354  *
355  * @param[in] logE Log10 of the true photon energy (TeV).
356  * @param[in] theta Offset angle in camera system (rad). Not used.
357  * @param[in] phi Azimuth angle in camera system (rad). Not used.
358  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
359  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
360  * @param[in] etrue Use true energy (true/false). Not used.
361  *
362  * Determine the radius beyond which the PSF becomes negligible. This radius
363  * is set by this method to \f$5 \times \sigma\f$, where \f$\sigma\f$ is the
364  * Gaussian width of the largest PSF component.
365  ***************************************************************************/
366 double GCTAPsfPerfTable::delta_max(const double& logE,
367  const double& theta,
368  const double& phi,
369  const double& zenith,
370  const double& azimuth,
371  const bool& etrue) const
372 {
373  // Update the parameter cache
374  update(logE);
375 
376  // Compute maximum PSF radius
377  double radius = 5.0 * m_par_sigma;
378 
379  // Return maximum PSF radius
380  return radius;
381 }
382 
383 
384 /***********************************************************************//**
385  * @brief Return the radius that contains a fraction of the events (radians)
386  *
387  * @param[in] fraction of events (0.0-1.0)
388  * @param[in] logE Log10 of the true photon energy (TeV).
389  * @param[in] theta Offset angle in camera system (rad). Not used.
390  * @param[in] phi Azimuth angle in camera system (rad). Not used.
391  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
392  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
393  * @param[in] etrue Use true energy (true/false). Not used.
394  * @return Containment radius (radians).
395  *
396  * @exception GException::invalid_argument
397  * Invalid fraction specified.
398  *
399  * Evaluates:
400  *
401  * \f[
402  *
403  * radius = \sqrt{ \frac{\ln{\left( \frac{ fraction * m\_par\_width}
404  * {\pi*m\_par\_scale} +1\right)}}{m\_par\_width} }
405  *
406  * \f]
407  *
408  * which is derived by integrating
409  *
410  * \f[
411  *
412  * fraction = \int_{0}^{2\pi}\int_{0}^{radius}r *
413  * e^{ m\_par\_width * r^{2}}dr d\phi
414  *
415  * \f]
416  *
417  * and solving for radius.
418  *
419  * Calculate the radius from the center that contains 'fraction' percent
420  * of the events. fraction * 100. = Containment % .
421  ***************************************************************************/
422 double GCTAPsfPerfTable::containment_radius(const double& fraction,
423  const double& logE,
424  const double& theta,
425  const double& phi,
426  const double& zenith,
427  const double& azimuth,
428  const bool& etrue) const
429 {
430  // Check input argument
431  if (fraction <= 0.0 || fraction >= 1.0) {
432  std::string message = "Containment fraction "+
433  gammalib::str(fraction)+" must be between " +
434  "0.0 and 1.0, not inclusive.";
436  }
437 
438  // Update the parameter cache
439  update(logE);
440 
441  // Compute radius
442  double arg = fraction * m_par_width / ( gammalib::pi * m_par_scale);
443  double radius = std::sqrt(std::log(arg + 1.0) / m_par_width);
444 
445  // Return containement radius
446  return radius;
447 }
448 
449 
450 /***********************************************************************//**
451  * @brief Print point spread function information
452  *
453  * @param[in] chatter Chattiness (defaults to NORMAL).
454  * @return String containing point spread function information.
455  ***************************************************************************/
456 std::string GCTAPsfPerfTable::print(const GChatter& chatter) const
457 {
458  // Initialise result string
459  std::string result;
460 
461  // Continue only if chatter is not silent
462  if (chatter != SILENT) {
463 
464  // Compute energy boundaries in TeV
465  int num = m_logE.size();
466  double emin = std::pow(10.0, m_logE[0]);
467  double emax = std::pow(10.0, m_logE[num-1]);
468 
469  // Append header
470  result.append("=== GCTAPsfPerfTable ===");
471 
472  // Append information
473  result.append("\n"+gammalib::parformat("Filename")+m_filename);
474  result.append("\n"+gammalib::parformat("Number of energy bins") +
475  gammalib::str(num));
476  result.append("\n"+gammalib::parformat("Log10(Energy) range"));
477  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
478 
479  } // endif: chatter was not silent
480 
481  // Return result
482  return result;
483 }
484 
485 
486 /*==========================================================================
487  = =
488  = Private methods =
489  = =
490  ==========================================================================*/
491 
492 /***********************************************************************//**
493  * @brief Initialise class members
494  ***************************************************************************/
496 {
497  // Initialise members
498  m_filename.clear();
499  m_logE.clear();
500  m_r68.clear();
501  m_r80.clear();
502  m_sigma.clear();
503  m_par_logE = -1.0e30;
504  m_par_scale = 1.0;
505  m_par_sigma = 0.0;
506  m_par_width = 0.0;
507 
508  // Return
509  return;
510 }
511 
512 
513 /***********************************************************************//**
514  * @brief Copy class members
515  *
516  * @param[in] psf Point spread function.
517  ***************************************************************************/
519 {
520  // Copy members
521  m_filename = psf.m_filename;
522  m_logE = psf.m_logE;
523  m_r68 = psf.m_r68;
524  m_r80 = psf.m_r80;
525  m_sigma = psf.m_sigma;
526  m_par_logE = psf.m_par_logE;
527  m_par_scale = psf.m_par_scale;
528  m_par_sigma = psf.m_par_sigma;
529  m_par_width = psf.m_par_width;
530 
531  // Return
532  return;
533 }
534 
535 
536 /***********************************************************************//**
537  * @brief Delete class members
538  ***************************************************************************/
540 {
541  // Return
542  return;
543 }
544 
545 
546 /***********************************************************************//**
547  * @brief Update PSF parameter cache
548  *
549  * @param[in] logE Log10 of the true photon energy (TeV).
550  *
551  * This method updates the PSF parameter cache. As the performance table PSF
552  * only depends on energy, the only parameter on which the cache values
553  * depend is the energy.
554  ***************************************************************************/
555 void GCTAPsfPerfTable::update(const double& logE) const
556 {
557  // Only compute PSF parameters if arguments have changed
558  if (logE != m_par_logE) {
559 
560  // Save energy
561  m_par_logE = logE;
562 
563  // Determine Gaussian sigma in radians
565 
566  // Derive width=-0.5/(sigma*sigma) and scale=1/(twopi*sigma*sigma)
567  double sigma2 = m_par_sigma * m_par_sigma;
568  m_par_scale = 1.0 / (gammalib::twopi * sigma2);
569  m_par_width = -0.5 / sigma2;
570 
571  }
572 
573  // Return
574  return;
575 }
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:188
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:862
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:73
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:1268
#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:1184
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 exception handler interface definition.
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:1332
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
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:1022
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:805
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
double m_par_logE
Energy for which precomputation is done.