GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAPsfVector.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAPsfVector.cpp - CTA point spread function vector 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 GCTAPsfVector.cpp
23  * @brief CTA point spread function vector class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GRan.hpp"
36 #include "GFilename.hpp"
37 #include "GFits.hpp"
38 #include "GFitsTable.hpp"
39 #include "GFitsTableCol.hpp"
40 #include "GCTAPsfVector.hpp"
41 
42 /* __ Method name definitions ____________________________________________ */
43 #define G_CONTAINMENT_RADIUS "GCTAPsfVector::containment_radius(double&,"\
44  " double&, double&, double&, double&, double&, bool&)"
45 
46 /* __ Macros _____________________________________________________________ */
47 
48 /* __ Coding definitions _________________________________________________ */
49 #define G_SMOOTH_PSF
50 
51 /* __ Debug definitions __________________________________________________ */
52 
53 /* __ Constants __________________________________________________________ */
54 
55 
56 /*==========================================================================
57  = =
58  = Constructors/destructors =
59  = =
60  ==========================================================================*/
61 
62 /***********************************************************************//**
63  * @brief Void constructor
64  ***************************************************************************/
66 {
67  // Initialise class members
68  init_members();
69 
70  // Return
71  return;
72 }
73 
74 
75 /***********************************************************************//**
76  * @brief File constructor
77  *
78  * @param[in] filename PSF FITS file.
79  *
80  * Constructs point spread function vector from a FITS file.
81  ***************************************************************************/
83 {
84  // Initialise class members
85  init_members();
86 
87  // Load point spread function from file
88  load(filename);
89 
90  // Return
91  return;
92 }
93 
94 
95 /***********************************************************************//**
96  * @brief Copy constructor
97  *
98  * @param[in] psf Point spread function.
99  ***************************************************************************/
101 {
102  // Initialise class members
103  init_members();
104 
105  // Copy members
106  copy_members(psf);
107 
108  // Return
109  return;
110 }
111 
112 
113 /***********************************************************************//**
114  * @brief Destructor
115  ***************************************************************************/
117 {
118  // Free members
119  free_members();
120 
121  // Return
122  return;
123 }
124 
125 
126 /*==========================================================================
127  = =
128  = Operators =
129  = =
130  ==========================================================================*/
131 
132 /***********************************************************************//**
133  * @brief Assignment operator
134  *
135  * @param[in] psf Point spread function.
136  * @return Point spread function.
137  ***************************************************************************/
139 {
140  // Execute only if object is not identical
141  if (this != &psf) {
142 
143  // Copy base class members
144  this->GCTAPsf::operator=(psf);
145 
146  // Free members
147  free_members();
148 
149  // Initialise private members
150  init_members();
151 
152  // Copy members
153  copy_members(psf);
154 
155  } // endif: object was not identical
156 
157  // Return this object
158  return *this;
159 }
160 
161 
162 /***********************************************************************//**
163  * @brief Return point spread function (in units of sr^-1)
164  *
165  * @param[in] delta Angular separation between true and measured photon
166  * directions (rad).
167  * @param[in] logE Log10 of the true photon energy (TeV).
168  * @param[in] theta Offset angle in camera system (rad). Not used.
169  * @param[in] phi Azimuth angle in camera system (rad). Not used.
170  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
171  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
172  * @param[in] etrue Use true energy (true/false). Not used.
173  *
174  * Returns the point spread function for a given angular separation in units
175  * of sr^-1 for a given energy.
176  ***************************************************************************/
177 double GCTAPsfVector::operator()(const double& delta,
178  const double& logE,
179  const double& theta,
180  const double& phi,
181  const double& zenith,
182  const double& azimuth,
183  const bool& etrue) const
184 {
185  #if defined(G_SMOOTH_PSF)
186  // Compute offset so that PSF goes to 0 at 5 times the sigma value. This
187  // is a kluge to get a PSF that smoothly goes to zero at the edge, which
188  // prevents steps or kinks in the log-likelihood function.
189  static const double offset = std::exp(-0.5*5.0*5.0);
190  #endif
191 
192  // Update the parameter cache
193  update(logE);
194 
195  // Compute PSF value
196  #if defined(G_SMOOTH_PSF)
197  double psf = m_par_scale * (std::exp(m_par_width * delta * delta) - offset);
198  #else
199  double psf = m_par_scale * (std::exp(m_par_width * delta * delta));
200  #endif
201 
202  #if defined(G_SMOOTH_PSF)
203  // Make sure that PSF is non-negative
204  if (psf < 0.0) {
205  psf = 0.0;
206  }
207  #endif
208 
209  // Return PSF
210  return psf;
211 }
212 
213 
214 /*==========================================================================
215  = =
216  = Public methods =
217  = =
218  ==========================================================================*/
219 
220 /***********************************************************************//**
221  * @brief Clear instance
222  *
223  * This method properly resets the object to an initial state.
224  ***************************************************************************/
226 {
227  // Free class members (base and derived classes, derived class first)
228  free_members();
229  this->GCTAPsf::free_members();
230 
231  // Initialise members
232  this->GCTAPsf::init_members();
233  init_members();
234 
235  // Return
236  return;
237 }
238 
239 
240 /***********************************************************************//**
241  * @brief Clone instance
242  *
243  * @return Deep copy of point spread function instance.
244  ***************************************************************************/
246 {
247  return new GCTAPsfVector(*this);
248 }
249 
250 
251 /***********************************************************************//**
252  * @brief Load point spread function from FITS file
253  *
254  * @param[in] filename FITS file name.
255  *
256  * Loads the point spread function from a FITS file.
257  *
258  * If no extension name is provided, the point spread function will be loaded
259  * from the "PSF" extension.
260  ***************************************************************************/
261 void GCTAPsfVector::load(const GFilename& filename)
262 {
263  // Open FITS file
264  GFits fits(filename);
265 
266  // Get PSFa table
267  const GFitsTable& table = *fits.table(filename.extname("PSF"));
268 
269  // Read PSF from table
270  read(table);
271 
272  // Close FITS file
273  fits.close();
274 
275  // Store filename
277 
278  // Return
279  return;
280 }
281 
282 
283 /***********************************************************************//**
284  * @brief Read point spread function vector from FITS table
285  *
286  * @param[in] table FITS table.
287  *
288  * Reads a point spread function vector from the FITS @p table.
289  *
290  * The energies are converted to TeV. Conversion is done based on the units
291  * provided for the energy columns. Units that are recognized are 'keV',
292  * 'MeV', 'GeV', and 'TeV' (case independent).
293  *
294  * The Gaussian width parameter may be either given in 68% containment radius
295  * (R68) or in sigma (ANGRES40; corresponding to a 38% containment radius).
296  * The former format is used for CTA and H.E.S.S. data, the latter for
297  * MAGIC data. All these things should be more uniform once we have a well
298  * defined format.
299  ***************************************************************************/
300 void GCTAPsfVector::read(const GFitsTable& table)
301 {
302  // Clear arrays
303  m_logE.clear();
304  m_r68.clear();
305  m_sigma.clear();
306 
307  // Set conversion factor from 68% containment radius to 1 sigma
308  const double conv = 0.6624305 * gammalib::deg2rad;
309 
310  // Get pointers to table columns
311  const GFitsTableCol* energy_lo = table["ENERG_LO"];
312  const GFitsTableCol* energy_hi = table["ENERG_HI"];
313 
314  // Handle various data formats (H.E.S.S. and MAGIC)
315  const GFitsTableCol* r68;
316  double r68_scale = 1.0;
317  if (table.contains("R68")) {
318  r68 = table["R68"];
319  }
320  else {
321  r68 = table["ANGRES40"];
322  r68_scale = 1.0 / 0.6624305; // MAGIC PSF is already 1 sigma
323  }
324 
325  // Determine unit conversion factors (default: TeV)
326  std::string u_energy_lo = gammalib::tolower(gammalib::strip_whitespace(energy_lo->unit()));
327  std::string u_energy_hi = gammalib::tolower(gammalib::strip_whitespace(energy_hi->unit()));
328  double c_energy_lo = 1.0;
329  double c_energy_hi = 1.0;
330  if (u_energy_lo == "kev") {
331  c_energy_lo = 1.0e-9;
332  }
333  else if (u_energy_lo == "mev") {
334  c_energy_lo = 1.0e-6;
335  }
336  else if (u_energy_lo == "gev") {
337  c_energy_lo = 1.0e-3;
338  }
339  if (u_energy_hi == "kev") {
340  c_energy_hi = 1.0e-9;
341  }
342  else if (u_energy_hi == "mev") {
343  c_energy_hi = 1.0e-6;
344  }
345  else if (u_energy_hi == "gev") {
346  c_energy_hi = 1.0e-3;
347  }
348 
349  // Extract number of energy bins
350  int num = energy_lo->nrows();
351 
352  // Set nodes
353  for (int i = 0; i < num; ++i) {
354 
355  // Compute log10 mean energy in TeV
356  double e_min = energy_lo->real(i) * c_energy_lo;
357  double e_max = energy_hi->real(i) * c_energy_hi;
358  double logE = 0.5 * (log10(e_min) + log10(e_max));
359 
360  // Extract r68 value and scale as required
361  double r68_value = r68->real(i) * r68_scale;
362 
363  // Store log10 mean energy and r68 value
364  m_logE.append(logE);
365  m_r68.push_back(r68_value);
366  m_sigma.push_back(r68_value*conv);
367 
368  } // endfor: looped over nodes
369 
370  // Return
371  return;
372 }
373 
374 
375 /***********************************************************************//**
376  * @brief Simulate PSF offset (radians)
377  *
378  * @param[in] ran Random number generator.
379  * @param[in] logE Log10 of the true photon energy (TeV).
380  * @param[in] theta Offset angle in camera system (rad). Not used.
381  * @param[in] phi Azimuth angle in camera system (rad). Not used.
382  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
383  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
384  * @param[in] etrue Use true energy (true/false). Not used.
385  ***************************************************************************/
387  const double& logE,
388  const double& theta,
389  const double& phi,
390  const double& zenith,
391  const double& azimuth,
392  const bool& etrue) const
393 {
394  // Update the parameter cache
395  update(logE);
396 
397  // Draw offset
398  double delta = m_par_sigma * ran.chisq2();
399 
400  // Return PSF offset
401  return delta;
402 }
403 
404 
405 /***********************************************************************//**
406  * @brief Return maximum size of PSF (radians)
407  *
408  * @param[in] logE Log10 of the true photon energy (TeV).
409  * @param[in] theta Offset angle in camera system (rad). Not used.
410  * @param[in] phi Azimuth angle in camera system (rad). Not used.
411  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
412  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
413  * @param[in] etrue Use true energy (true/false). Not used.
414  *
415  * Determine the radius beyond which the PSF becomes negligible. This radius
416  * is set by this method to \f$5 \times \sigma\f$, where \f$\sigma\f$ is the
417  * Gaussian width of the largest PSF component.
418  ***************************************************************************/
419 double GCTAPsfVector::delta_max(const double& logE,
420  const double& theta,
421  const double& phi,
422  const double& zenith,
423  const double& azimuth,
424  const bool& etrue) const
425 {
426  // Update the parameter cache
427  update(logE);
428 
429  // Compute maximum PSF radius
430  double radius = 5.0 * m_par_sigma;
431 
432  // Return maximum PSF radius
433  return radius;
434 }
435 
436 
437 /***********************************************************************//**
438  * @brief Return the radius that contains a fraction of the events (radians)
439  *
440  * @param[in] fraction of events (0.0-1.0)
441  * @param[in] logE Log10 of the true photon energy (TeV).
442  * @param[in] theta Offset angle in camera system (rad). Not used.
443  * @param[in] phi Azimuth angle in camera system (rad). Not used.
444  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
445  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
446  * @param[in] etrue Use true energy (true/false). Not used.
447  * @return Containment radius (radians).
448  *
449  * @exception GException::invalid_argument
450  * Invalid fraction specified.
451  *
452  * Evaluates:
453  *
454  * \f[
455  *
456  * radius = \sqrt{ \frac{\ln{\left( \frac{ fraction * m\_par\_width}
457  * {\pi*m\_par\_scale} +1\right)}}{m\_par\_width} }
458  *
459  * \f]
460  *
461  * which is derived by integrating
462  *
463  * \f[
464  *
465  * fraction = \int_{0}^{2\pi}\int_{0}^{radius}r *
466  * e^{ m\_par\_width * r^{2}}dr d\phi
467  *
468  * \f]
469  *
470  * and solving for radius.
471  *
472  * Calculate the radius from the center that contains 'fraction' percent
473  * of the events. fraction * 100. = Containment % .
474  ***************************************************************************/
475 double GCTAPsfVector::containment_radius(const double& fraction,
476  const double& logE,
477  const double& theta,
478  const double& phi,
479  const double& zenith,
480  const double& azimuth,
481  const bool& etrue) const
482 {
483  // Check input argument
484  if (fraction <= 0.0 || fraction >= 1.0) {
485  std::string msg = "Containment fraction "+gammalib::str(fraction)+
486  " must lie in interval ]0,1[. Please specify a "
487  "valid containment fraction";
489  }
490 
491  // Update the parameter cache
492  update(logE);
493 
494  // Compute radius for gaussian
495  double arg = fraction * m_par_width / (gammalib::pi * m_par_scale);
496  double radius = std::sqrt(std::log(arg + 1.0) / m_par_width);
497 
498  // Return containement radius
499  return radius;
500 }
501 
502 
503 /***********************************************************************//**
504  * @brief Print point spread function information
505  *
506  * @param[in] chatter Chattiness.
507  * @return String containing point spread function information.
508  ***************************************************************************/
509 std::string GCTAPsfVector::print(const GChatter& chatter) const
510 {
511  // Initialise result string
512  std::string result;
513 
514  // Continue only if chatter is not silent
515  if (chatter != SILENT) {
516 
517  // Compute energy boundaries in TeV
518  int num = m_logE.size();
519  double emin = std::pow(10.0, m_logE[0]);
520  double emax = std::pow(10.0, m_logE[num-1]);
521 
522  // Append header
523  result.append("=== GCTAPsfVector ===");
524 
525  // Append information
526  result.append("\n"+gammalib::parformat("Filename")+m_filename);
527  result.append("\n"+gammalib::parformat("Number of energy bins") +
528  gammalib::str(num));
529  result.append("\n"+gammalib::parformat("Energy range"));
530  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
531 
532  } // endif: chatter was not silent
533 
534  // Return result
535  return result;
536 }
537 
538 
539 /*==========================================================================
540  = =
541  = Private methods =
542  = =
543  ==========================================================================*/
544 
545 /***********************************************************************//**
546  * @brief Initialise class members
547  ***************************************************************************/
549 {
550  // Initialise members
551  m_filename.clear();
552  m_logE.clear();
553  m_r68.clear();
554  m_sigma.clear();
555  m_par_logE = -1.0e30;
556  m_par_scale = 1.0;
557  m_par_sigma = 0.0;
558  m_par_width = 0.0;
559 
560  // Return
561  return;
562 }
563 
564 
565 /***********************************************************************//**
566  * @brief Copy class members
567  *
568  * @param[in] psf Point spread function.
569  ***************************************************************************/
571 {
572  // Copy members
573  m_filename = psf.m_filename;
574  m_logE = psf.m_logE;
575  m_r68 = psf.m_r68;
576  m_sigma = psf.m_sigma;
577  m_par_logE = psf.m_par_logE;
578  m_par_scale = psf.m_par_scale;
579  m_par_sigma = psf.m_par_sigma;
580  m_par_width = psf.m_par_width;
581 
582  // Return
583  return;
584 }
585 
586 
587 /***********************************************************************//**
588  * @brief Delete class members
589  ***************************************************************************/
591 {
592  // Return
593  return;
594 }
595 
596 
597 /***********************************************************************//**
598  * @brief Update PSF parameter cache
599  *
600  * @param[in] logE Log10 of the true photon energy (TeV).
601  *
602  * This method updates the PSF parameter cache. As the performance table PSF
603  * only depends on energy, the only parameter on which the cache values
604  * depend is the energy.
605  ***************************************************************************/
606 void GCTAPsfVector::update(const double& logE) const
607 {
608  // Only compute PSF parameters if arguments have changed
609  if (logE != m_par_logE) {
610 
611  // Save energy
612  m_par_logE = logE;
613 
614  // Determine Gaussian sigma in radians
616 
617  // Derive width=-0.5/(sigma*sigma) and scale=1/(twopi*sigma*sigma)
618  double sigma2 = m_par_sigma * m_par_sigma;
619  m_par_scale = 1.0 / (gammalib::twopi * sigma2);
620  m_par_width = -0.5 / sigma2;
621 
622  }
623 
624  // Return
625  return;
626 }
double m_par_width
Gaussian width parameter.
void free_members(void)
Delete class members.
Definition: GCTAPsf.cpp:164
void unit(const std::string &unit)
Set column unit.
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
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
Random number generator class definition.
const double pi
Definition: GMath.hpp:35
CTA vector point spread function class.
void load(const GFilename &filename)
Load point spread function from FITS file.
bool contains(const std::string &colname) const
Checks the presence of a column in table.
Definition: GFitsTable.cpp:759
Abstract base class for the CTA point spread function.
Definition: GCTAPsf.hpp:47
std::vector< double > m_sigma
Sigma value of PSF in radians.
void nrows(const int &nrows)
Set number of rows in column.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
Random number generator class.
Definition: GRan.hpp:44
void free_members(void)
Delete class members.
Gammalib tools definition.
void clear(void)
Clear instance.
FITS file class.
Definition: GFits.hpp:63
GCTAPsf & operator=(const GCTAPsf &psf)
Assignment operator.
Definition: GCTAPsf.cpp:106
FITS file class interface definition.
double m_par_logE
Energy for which precomputation is done.
void copy_members(const GCTAPsfVector &psf)
Copy class members.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:80
GNodeArray m_logE
log(E) nodes for Aeff interpolation
FITS table column abstract base class definition.
void init_members(void)
Initialise class members.
Definition: GCTAPsf.cpp:142
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function information.
virtual ~GCTAPsfVector(void)
Destructor.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
void update(const double &logE) const
Update PSF parameter cache.
Filename class.
Definition: GFilename.hpp:62
Abstract interface for FITS table column.
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
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GChatter
Definition: GTypemaps.hpp:33
GFilename filename(void) const
Return filename.
#define G_CONTAINMENT_RADIUS
GCTAPsfVector * clone(void) const
Clone instance.
GFilename m_filename
Name of Aeff response file.
CTA point spread function vector class definition.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
GCTAPsfVector & operator=(const GCTAPsfVector &psf)
Assignment operator.
virtual double real(const int &row, const int &inx=0) const =0
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
void init_members(void)
Initialise class members.
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)
Exception handler interface definition.
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)
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)
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
double m_par_scale
Gaussian normalization.
const double twopi
Definition: GMath.hpp:36
void read(const GFitsTable &table)
Read point spread function vector from FITS table.
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
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
double m_par_sigma
Gaussian sigma (radians)
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)
std::vector< double > m_r68
68% containment radius of PSF in degrees
Filename class interface definition.
Mathematical function definitions.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1295
GCTAPsfVector(void)
Void constructor.
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.