GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAPsfTable.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAPsfTable.cpp - CTA point spread function table class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2016-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 GCTAPsfTable.cpp
23  * @brief CTA point spread function table 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 "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GException.hpp"
35 #include "GRan.hpp"
36 #include "GFilename.hpp"
37 #include "GFits.hpp"
38 #include "GFitsTable.hpp"
39 #include "GFitsBinTable.hpp"
40 #include "GCTAPsfTable.hpp"
41 #include "GCTASupport.hpp"
42 
43 /* __ Method name definitions ____________________________________________ */
44 #define G_READ "GCTAPsfTable::read(GFitsTable&)"
45 #define G_CONTAINMENT_RADIUS "GCTAPsfTable::containment_radius(double&,"\
46  " double&, double&, double&, double&, double&, bool&)"
47 
48 /* __ Macros _____________________________________________________________ */
49 
50 /* __ Coding definitions _________________________________________________ */
51 
52 /* __ Debug definitions __________________________________________________ */
53 
54 /* __ Constants __________________________________________________________ */
55 
56 
57 /*==========================================================================
58  = =
59  = Constructors/destructors =
60  = =
61  ==========================================================================*/
62 
63 /***********************************************************************//**
64  * @brief Void constructor
65  *
66  * Constructs empty point spread function.
67  ***************************************************************************/
69 {
70  // Initialise class members
71  init_members();
72 
73  // Return
74  return;
75 }
76 
77 
78 /***********************************************************************//**
79  * @brief File constructor
80  *
81  * @param[in] filename PSF FITS file.
82  *
83  * Constructs point spread function from a FITS file.
84  ***************************************************************************/
86 {
87  // Initialise class members
88  init_members();
89 
90  // Load point spread function from file
91  load(filename);
92 
93  // Return
94  return;
95 }
96 
97 
98 /***********************************************************************//**
99  * @brief Copy constructor
100  *
101  * @param[in] psf Point spread function.
102  *
103  * Constructs point spread function by copying from another point spread
104  * function.
105  ***************************************************************************/
107 {
108  // Initialise class members
109  init_members();
110 
111  // Copy members
112  copy_members(psf);
113 
114  // Return
115  return;
116 }
117 
118 
119 /***********************************************************************//**
120  * @brief Destructor
121  *
122  * Destructs point spread function.
123  ***************************************************************************/
125 {
126  // Free members
127  free_members();
128 
129  // Return
130  return;
131 }
132 
133 
134 /*==========================================================================
135  = =
136  = Operators =
137  = =
138  ==========================================================================*/
139 
140 /***********************************************************************//**
141  * @brief Assignment operator
142  *
143  * @param[in] psf Point spread function.
144  * @return Point spread function.
145  *
146  * Assigns point spread function.
147  ***************************************************************************/
149 {
150  // Execute only if object is not identical
151  if (this != &psf) {
152 
153  // Copy base class members
154  this->GCTAPsf::operator=(psf);
155 
156  // Free members
157  free_members();
158 
159  // Initialise private members
160  init_members();
161 
162  // Copy members
163  copy_members(psf);
164 
165  } // endif: object was not identical
166 
167  // Return this object
168  return *this;
169 }
170 
171 
172 /***********************************************************************//**
173  * @brief Return point spread function (in units of sr^-1)
174  *
175  * @param[in] delta Angular separation between true and measured photon
176  * directions (rad).
177  * @param[in] logE Log10 of the true photon energy (TeV).
178  * @param[in] theta Offset angle in camera system (rad).
179  * @param[in] phi Azimuth angle in camera system (rad). Not used.
180  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
181  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
182  * @param[in] etrue Use true energy (true/false). Not used.
183  *
184  * Returns the point spread function for a given angular separation in units
185  * of sr^-1 for a given energy and offset angle. The PSF value will be
186  * determined by trilinear interpolation (and extrapolation) in the PSF
187  * histograms. If the interpolation or extrapolation would lead to negative
188  * PSF values, the value of zero is returned. A zero value is also returned
189  * if the @p delta angle is larger than the largest angle in the response
190  * table.
191  ***************************************************************************/
192 double GCTAPsfTable::operator()(const double& delta,
193  const double& logE,
194  const double& theta,
195  const double& phi,
196  const double& zenith,
197  const double& azimuth,
198  const bool& etrue) const
199 {
200  // Initialise PSF value
201  double psf = 0.0;
202 
203  // Continue only if delta is not larger than delta_max
204  if (delta <= m_delta_max) {
205 
206  // Setup argument vector
207  double arg[3];
208  arg[m_inx_energy] = logE;
209  arg[m_inx_theta] = theta;
210  arg[m_inx_delta] = delta;
211 
212  // Compute point spread function by trilinear interpolation
213  psf = m_psf(m_inx_rpsf, arg[0], arg[1], arg[2]);
214 
215  // Constrain PSF value to non-negative values
216  if (psf < 0.0) {
217  psf = 0.0;
218  }
219 
220  } // endif: delta was within validity range
221 
222  // Return PSF
223  return psf;
224 }
225 
226 
227 /*==========================================================================
228  = =
229  = Public methods =
230  = =
231  ==========================================================================*/
232 
233 /***********************************************************************//**
234  * @brief Clear point spread function
235  *
236  * Clears point spread function.
237  ***************************************************************************/
239 {
240  // Free class members (base and derived classes, derived class first)
241  free_members();
242  this->GCTAPsf::free_members();
243 
244  // Initialise members
245  this->GCTAPsf::init_members();
246  init_members();
247 
248  // Return
249  return;
250 }
251 
252 
253 /***********************************************************************//**
254  * @brief Clone point spread functions
255  *
256  * @return Deep copy of point spread function.
257  *
258  * Returns a pointer to a deep copy of the point spread function.
259  ***************************************************************************/
261 {
262  return new GCTAPsfTable(*this);
263 }
264 
265 
266 /***********************************************************************//**
267  * @brief Read point spread function from FITS table
268  *
269  * @param[in] table FITS table.
270  *
271  * @exception GException::invalid_value
272  * Response table is not three-dimensional.
273  *
274  * Reads the point spread function form the FITS @p table. The following
275  * column names are mandatory:
276  *
277  * ENERG_LO - Energy lower bin boundaries
278  * ENERG_HI - Energy upper bin boundaries
279  * THETA_LO - Offset angle lower bin boundaries
280  * THETA_HI - Offset angle upper bin boundaries
281  * RAD_LO - Delta angle lower bin boundaries
282  * RAD_HI - Delta angle upper bin boundaries
283  * RPSF - PSF histogram
284  *
285  * The data are stored in the m_psf member. The energy axis will be set to
286  * log10, the offset and delta angle axes to radians.
287  ***************************************************************************/
288 void GCTAPsfTable::read(const GFitsTable& table)
289 {
290  // Clear response table
291  m_psf.clear();
292 
293  // Read PSF table
294  m_psf.read(table);
295 
296  // Get mandatory indices (throw exception if not found)
297  m_inx_energy = m_psf.axis("ENERG");
298  m_inx_theta = m_psf.axis("THETA");
299  m_inx_delta = m_psf.axis("RAD");
300  m_inx_rpsf = m_psf.table("RPSF");
301 
302  // Throw an exception if the table is not three-dimensional
303  if (m_psf.axes() != 3) {
304  std::string msg = "Expected three-dimensional point spread function "
305  "response table but found "+
307  " dimensions. Please specify a three-dimensional "
308  "point spread function.";
309  throw GException::invalid_value(G_READ, msg);
310  }
311 
312  // Set energy axis to logarithmic scale
314 
315  // Set offset angle axis to radians
317 
318  // Set delta angle axis to radians
320 
321  // Precompute PSF information
322  precompute();
323 
324  // Return
325  return;
326 }
327 
328 
329 /***********************************************************************//**
330  * @brief Write point spread function into FITS table
331  *
332  * @param[in] table FITS binary table.
333  *
334  * Writes point spread function into a FITS binary @p table.
335  *
336  * @todo Add keywords.
337  ***************************************************************************/
339 {
340  // Create a copy of the response table
342 
343  // Write response table
344  psf.write(table);
345 
346  // Return
347  return;
348 }
349 
350 
351 /***********************************************************************//**
352  * @brief Load point spread function from FITS file
353  *
354  * @param[in] filename FITS file name.
355  *
356  * Loads the point spread function from a FITS file.
357  *
358  * If no extension name is given the method scans the `HDUCLASS` keywords
359  * of all extensions and loads the background from the first extension
360  * for which `HDUCLAS4=PSF_TABLE`.
361  *
362  * Otherwise, the background will be loaded from the `POINT SPREAD FUNCTION`
363  * extension.
364  ***************************************************************************/
365 void GCTAPsfTable::load(const GFilename& filename)
366 {
367  // Open FITS file
368  GFits fits(filename);
369 
370  // Get the default extension name. If no GADF compliant name was found
371  // then set the default extension name to "POINT SPREAD FUNCTION".
372  std::string extname = gammalib::gadf_hduclas4(fits, "PSF_TABLE");
373  if (extname.empty()) {
375  }
376 
377  // Get PSF table
378  const GFitsTable& table = *fits.table(filename.extname(extname));
379 
380  // Read PSF from table
381  read(table);
382 
383  // Close FITS file
384  fits.close();
385 
386  // Store filename
388 
389  // Return
390  return;
391 }
392 
393 /***********************************************************************//**
394  * @brief Save point spread function table into FITS file
395  *
396  * @param[in] filename FITS file name.
397  * @param[in] clobber Overwrite existing file?
398  *
399  * Saves point spread function into a FITS file. If a file with the given
400  * @p filename does not yet exist it will be created, otherwise the method
401  * opens the existing file. The method will create a (or replace an existing)
402  * point spread function extension. The extension name can be specified as
403  * part of the @p filename, or if no extension name is given, is assumed to
404  * be `POINT SPREAD FUNCTION`.
405  *
406  * An existing file will only be modified if the @p clobber flag is set to
407  * true.
408  ***************************************************************************/
409 void GCTAPsfTable::save(const GFilename& filename, const bool& clobber) const
410 {
411  // Get extension name
412  std::string extname = filename.extname(gammalib::extname_cta_psftable);
413 
414  // Open or create FITS file (without extension name since the requested
415  // extension may not yet exist in the file)
416  GFits fits(filename.url(), true);
417 
418  // Remove extension if it exists already
419  if (fits.contains(extname)) {
420  fits.remove(extname);
421  }
422 
423  // Create binary table
425 
426  // Write the background table
427  write(table);
428 
429  // Set binary table extension name
430  table.extname(extname);
431 
432  // Append table to FITS file
433  fits.append(table);
434 
435  // Save to file
436  fits.save(clobber);
437 
438  // Return
439  return;
440 }
441 
442 
443 /***********************************************************************//**
444  * @brief Simulate PSF offset (radians)
445  *
446  * @param[in] ran Random number generator.
447  * @param[in] logE Log10 of the true photon energy (TeV).
448  * @param[in] theta Offset angle in camera system (rad).
449  * @param[in] phi Azimuth angle in camera system (rad). Not used.
450  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
451  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
452  * @param[in] etrue Use true energy (true/false). Not used.
453  *
454  * Draws a random offset from the point spread function using a rejection
455  * method.
456  ***************************************************************************/
457 double GCTAPsfTable::mc(GRan& ran,
458  const double& logE,
459  const double& theta,
460  const double& phi,
461  const double& zenith,
462  const double& azimuth,
463  const bool& etrue) const
464 {
465  // Initialise random offset
466  double delta = 0.0;
467 
468  // Simulate PSF only if maximum PSF radius is positive
469  if (m_delta_max > 0.0 && m_psf_max > 0.0) {
470 
471  // Pre-compute 1-cos(delta_max)
472  double cosrad = 1.0 - std::cos(m_delta_max);
473 
474  // Draw random positions
475  while (true) {
476 
477  // Simulate offset angle
478  delta = std::acos(1.0 - ran.uniform() * cosrad);
479 
480  // Get corresponding PSF value
481  double value = this->operator()(delta, logE, theta, phi, zenith,
482  azimuth, etrue);
483 
484  // Get uniform random number up to the maximum
485  double uniform = ran.uniform() * m_psf_max;
486 
487  // Exit loop if the random number is not larger than the PSF value
488  if (uniform <= value) {
489  break;
490  }
491 
492  } // endwhile: rejection method
493 
494  } // endif: PSF radius and maximum PSF value was positive
495 
496  // Return PSF offset
497  return delta;
498 }
499 
500 
501 /***********************************************************************//**
502  * @brief Return maximum size of PSF (radians)
503  *
504  * @param[in] logE Log10 of the true photon energy (TeV).
505  * @param[in] theta Offset angle in camera system (rad).
506  * @param[in] phi Azimuth angle in camera system (rad). Not used.
507  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
508  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
509  * @param[in] etrue Use true energy (true/false). Not used.
510  *
511  * Returns the maximum PSF radius.
512  ***************************************************************************/
513 double GCTAPsfTable::delta_max(const double& logE,
514  const double& theta,
515  const double& phi,
516  const double& zenith,
517  const double& azimuth,
518  const bool& etrue) const
519 {
520  // Return maximum PSF radius
521  return m_delta_max;
522 }
523 
524 
525 /***********************************************************************//**
526  * @brief Return the radius that contains a fraction of the events (radians)
527  *
528  * @param[in] fraction of events (0.0-1.0)
529  * @param[in] logE Log10 of the true photon energy (TeV).
530  * @param[in] theta Offset angle in camera system (rad). Not used.
531  * @param[in] phi Azimuth angle in camera system (rad). Not used.
532  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
533  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
534  * @param[in] etrue Use true energy (true/false). Not used.
535  * @return Containment radius (radians).
536  *
537  * @exception GException::invalid_argument
538  * Invalid fraction specified.
539  *
540  * Calculate the radius from the center that contains 'fraction' percent
541  * of the events. fraction * 100. = Containment %.
542  ***************************************************************************/
543 double GCTAPsfTable::containment_radius(const double& fraction,
544  const double& logE,
545  const double& theta,
546  const double& phi,
547  const double& zenith,
548  const double& azimuth,
549  const bool& etrue) const
550 {
551  // Check input argument
552  if (fraction <= 0.0 || fraction >= 1.0) {
553  std::string message = "Containment fraction "+
554  gammalib::str(fraction)+" must be between " +
555  "0.0 and 1.0, not inclusive.";
557  }
558 
559  // Determine number of delta bins
560  int ndelta = m_psf.axis_bins(m_inx_delta);
561 
562  // Initialise containment radius and PSF integral
563  double delta = 0.0;
564  double sum = 0.0;
565 
566  // Integrate delta array by computing the sum over the pixels times the
567  // pixel size times the delta angle
568  for (int idelta = 0; idelta < ndelta; ++idelta) {
569 
570  // Compute delta value (radians)
571  delta = m_psf.axis_nodes(m_inx_delta)[idelta];
572 
573  // Compute delta bin width (radians)
574  double width = (m_psf.axis_hi(m_inx_delta, idelta) -
575  m_psf.axis_lo(m_inx_delta, idelta)) *
577 
578  // Integrate PSF
579  sum += this->operator()(delta, logE, theta, phi, zenith, azimuth, etrue) *
580  std::sin(delta) * width * gammalib::twopi;
581 
582  // Break of the containment fraction is reached or exceeded
583  if (sum >= fraction) {
584  break;
585  }
586 
587  } // endfor: loop over delta array
588 
589  // Return containing radius
590  return delta;
591 }
592 
593 
594 /***********************************************************************//**
595  * @brief Print point spread function information
596  *
597  * @param[in] chatter Chattiness.
598  * @return String containing point spread function information.
599  ***************************************************************************/
600 std::string GCTAPsfTable::print(const GChatter& chatter) const
601 {
602  // Initialise result string
603  std::string result;
604 
605  // Continue only if chatter is not silent
606  if (chatter != SILENT) {
607 
608  // Append header
609  result.append("=== GCTAPsfTable ===");
610  result.append("\n"+gammalib::parformat("Filename")+m_filename);
611 
612  // Initialise information
613  int nebins = 0;
614  int nthetabins = 0;
615  int ndeltabins = 0;
616  double emin = 0.0;
617  double emax = 0.0;
618  double omin = 0.0;
619  double omax = 0.0;
620  double dmin = 0.0;
621  double dmax = 0.0;
622 
623  // Extract information if there are axes in the response table
624  if (m_psf.axes() > 0) {
625  nebins = m_psf.axis_bins(m_inx_energy);
626  nthetabins = m_psf.axis_bins(m_inx_theta);
627  ndeltabins = m_psf.axis_bins(m_inx_delta);
628  emin = m_psf.axis_lo(m_inx_energy,0);
629  emax = m_psf.axis_hi(m_inx_energy,nebins-1);
630  omin = m_psf.axis_lo(m_inx_theta,0);
631  omax = m_psf.axis_hi(m_inx_theta,nthetabins-1);
632  dmin = m_psf.axis_lo(m_inx_delta,0);
633  dmax = m_psf.axis_hi(m_inx_delta,ndeltabins-1);
634  }
635 
636 
637  // Append information
638  result.append("\n"+gammalib::parformat("Number of energy bins") +
639  gammalib::str(nebins));
640  result.append("\n"+gammalib::parformat("Number of offset bins") +
641  gammalib::str(nthetabins));
642  result.append("\n"+gammalib::parformat("Number of delta bins") +
643  gammalib::str(ndeltabins));
644  result.append("\n"+gammalib::parformat("Energy range"));
645  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
646  result.append("\n"+gammalib::parformat("Offset angle range"));
647  result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
648  result.append("\n"+gammalib::parformat("Delta angle range"));
649  result.append(gammalib::str(dmin)+" - "+gammalib::str(dmax)+" deg");
650 
651  } // endif: chatter was not silent
652 
653  // Return result
654  return result;
655 }
656 
657 
658 /*==========================================================================
659  = =
660  = Private methods =
661  = =
662  ==========================================================================*/
663 
664 /***********************************************************************//**
665  * @brief Initialise class members
666  ***************************************************************************/
668 {
669  // Initialise members
670  m_filename.clear();
671  m_psf.clear();
672  m_inx_energy = 0;
673  m_inx_theta = 1;
674  m_inx_delta = 2;
675  m_inx_rpsf = 0;
676  m_delta_max = 0.0;
677  m_psf_max = 0.0;
678 
679  // Return
680  return;
681 }
682 
683 
684 /***********************************************************************//**
685  * @brief Copy class members
686  *
687  * @param[in] psf Point spread function.
688  ***************************************************************************/
690 {
691  // Copy members
692  m_filename = psf.m_filename;
693  m_psf = psf.m_psf;
695  m_inx_theta = psf.m_inx_theta;
696  m_inx_delta = psf.m_inx_delta;
697  m_inx_rpsf = psf.m_inx_rpsf;
698  m_delta_max = psf.m_delta_max;
699  m_psf_max = psf.m_psf_max;
700 
701  // Return
702  return;
703 }
704 
705 
706 /***********************************************************************//**
707  * @brief Delete class members
708  ***************************************************************************/
710 {
711  // Return
712  return;
713 }
714 
715 
716 /***********************************************************************//**
717  * @brief Performs precomputations for point spread function
718  *
719  * Replaces any invalid PSF histogram by zeroes, normalises the PSF
720  * histograms to unity, and computes maximum PSF value and maximum delta
721  * value.
722  ***************************************************************************/
724 {
725  // Determine PSF dimensions
726  int neng = m_psf.axis_bins(m_inx_energy);
727  int ntheta = m_psf.axis_bins(m_inx_theta);
728  int ndelta = m_psf.axis_bins(m_inx_delta);
729 
730  // Determine delta element increment
731  int inx_inc = element(0,0,1) - element(0,0,0);
732 
733  // Replace any negative or invalid histogram values by zero values
734  for (int i = 0; i < m_psf.elements(); ++i) {
735  double element = m_psf(m_inx_rpsf, i);
736  if (element < 0.0 ||
737  gammalib::is_notanumber(element) ||
738  gammalib::is_infinite(element)) {
739  m_psf(m_inx_rpsf, i) = 0.0;
740  }
741  }
742 
743  // Normalise PSF to unity
744  for (int ieng = 0; ieng < neng; ++ieng) {
745  for (int itheta = 0; itheta < ntheta; ++itheta) {
746 
747  // Initialise PSF integral
748  double sum = 0.0;
749 
750  // Set start element
751  int inx = element(ieng, itheta, 0);
752 
753  // Integrate delta array by computing the sum over the pixels
754  // times the pixel size time the delta angle
755  for (int idelta = 0; idelta < ndelta; ++idelta, inx += inx_inc) {
756 
757  // Compute delta value (radians)
758  double delta = m_psf.axis_nodes(m_inx_delta)[idelta];
759 
760  // Compute delta bin width (radians)
761  double width = (m_psf.axis_hi(m_inx_delta, idelta) -
762  m_psf.axis_lo(m_inx_delta, idelta)) *
764 
765  // Integrate PSF
766  sum += m_psf(m_inx_rpsf, inx) * std::sin(delta) * width *
768 
769  }
770 
771  // If integral is positive then divide PSF by integral so that it
772  // normalises to unity
773  if (sum > 0.0) {
774 
775  // Set start element
776  int inx = element(ieng, itheta, 0);
777 
778  // Normalise PSF
779  for (int idelta = 0; idelta < ndelta; ++idelta, inx += inx_inc) {
780  m_psf(m_inx_rpsf, inx) /= sum;
781  }
782 
783  } // endif: PSF integral was positive
784 
785  } // endfor: looped over all theta angles
786  } // endfor: looped over all PSF values
787 
788  // Compute maximum PSF value
789  m_psf_max = 0.0;
790  for (int i = 0; i < m_psf.elements(); ++i) {
791  if (m_psf(m_inx_rpsf, i) > m_psf_max) {
792  m_psf_max = m_psf(m_inx_rpsf, i);
793  }
794  }
795 
796  // Set maximum PSF radius (radians)
799 
800  // Return
801  return;
802 }
803 
804 
805 /***********************************************************************//**
806  * @brief Return element index
807  *
808  * @param[in] ieng Energy index [0,...,neng-1]
809  * @param[in] itheta Offset angle index [0,...,ntheta-1]
810  * @param[in] idelta Delta angle index [0,...,ndelta-1]
811  ***************************************************************************/
812 int GCTAPsfTable::element(const int& ieng, const int& itheta, const int& idelta)
813 {
814  // Setup index vector for requested element
815  int index[3];
816  index[m_inx_energy] = ieng;
817  index[m_inx_theta] = itheta;
818  index[m_inx_delta] = idelta;
819 
820  // Get element index
821  int element= index[0] + (index[1] + index[2]*m_psf.axis_bins(1)) *
822  m_psf.axis_bins(0);
823 
824  // Return element index
825  return element;
826 }
void free_members(void)
Delete class members.
Definition: GCTAPsf.cpp:164
const GNodeArray & axis_nodes(const int &axis) const
Return axis nodes.
int m_inx_energy
Energy index.
GCTAPsfTable & operator=(const GCTAPsfTable &psf)
Assignment operator.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
double m_psf_max
Maximum PSF value.
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function information.
Random number generator class definition.
GCTAResponseTable m_psf
PSF response table.
int m_inx_theta
Theta index.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
const GCTAResponseTable & table(void) const
Return response table.
#define G_READ
Abstract base class for the CTA point spread function.
Definition: GCTAPsf.hpp:47
GCTAPsfTable(void)
Void constructor.
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:944
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
Random number generator class.
Definition: GRan.hpp:44
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
void axis_radians(const int &axis)
Set nodes for a radians axis.
GCTAPsf & operator=(const GCTAPsf &psf)
Assignment operator.
Definition: GCTAPsf.cpp:106
FITS file class interface definition.
void copy_members(const GCTAPsfTable &psf)
Copy class members.
void read(const GFitsTable &table)
Read response table from FITS table HDU.
void save(const GFilename &filename, const bool &clobber=false) const
Save point spread function table into FITS file.
int m_inx_delta
Delta index.
void init_members(void)
Initialise class members.
Definition: GCTAPsf.cpp:142
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
Definition of support function used by CTA classes.
void write(GFitsTable &table) const
Write response table into FITS table HDU.
int axis(const std::string &name) const
Determine index of an axis.
void free_members(void)
Delete class members.
const double deg2rad
Definition: GMath.hpp:43
int element(const int &ieng, const int &itheta, const int &idelta)
Return element index.
void read(const GFitsTable &table)
Read point spread function from FITS table.
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:862
void clear(void)
Clear response table.
Filename class.
Definition: GFilename.hpp:62
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)
GFilename filename(void) const
Return filename.
int m_inx_rpsf
PSF histogram.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
const std::string extname_cta_psftable
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)
#define G_CONTAINMENT_RADIUS
const int & elements(void) const
Return number of elements per table.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void load(const GFilename &filename)
Load point spread function from FITS file.
GChatter
Definition: GTypemaps.hpp:33
virtual ~GCTAPsfTable(void)
Destructor.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void axis_log10(const int &axis)
Set nodes for a logarithmic (base 10) axis.
int table(const std::string &name) const
Determine index of table.
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
int axis_bins(const int &axis) const
Return number bins in an axis.
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)
GCTAPsfTable * clone(void) const
Clone point spread functions.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void write(GFitsBinTable &table) const
Write point spread function into FITS table.
FITS binary table class.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
Exception handler interface definition.
CTA point spread function table class definition.
FITS binary table class definition.
GFilename m_filename
Name of Aeff response file.
CTA response table class.
const double twopi
Definition: GMath.hpp:36
CTA point spread function table class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void precompute(void)
Performs precomputations for point spread function.
double m_delta_max
Maximum delta angle (radians)
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
const int & axes(void) const
Return number of axes of the tables.
Filename class interface definition.
void clear(void)
Clear point spread function.
Mathematical function definitions.
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.
void init_members(void)
Initialise class members.