GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAPsfKing.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAPsfKing.cpp - King profile CTA point spread function class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2018 by Michael Mayer *
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 GCTAPsfKing.hpp
23  * @brief King profile CTA point spread function class definition
24  * @author Michael Mayer
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 "GCTAPsfKing.hpp"
41 #include "GCTAException.hpp"
42 #include "GCTASupport.hpp"
43 
44 /* __ Method name definitions ____________________________________________ */
45 #define G_READ "GCTAPsfKing::read(GFitsTable&)"
46 #define G_CONTAINMENT_RADIUS "GCTAPsfKing::containment_radius(double&,"\
47  " double&, double&, double&, double&, double&, bool&)"
48 #define G_UPDATE "GCTAPsfKing::update(double&, double&)"
49 
50 /* __ Macros _____________________________________________________________ */
51 
52 /* __ Coding definitions _________________________________________________ */
53 #define G_SMOOTH_PSF
54 
55 /* __ Debug definitions __________________________________________________ */
56 
57 /* __ Constants __________________________________________________________ */
58 
59 
60 /*==========================================================================
61  = =
62  = Constructors/destructors =
63  = =
64  ==========================================================================*/
65 
66 /***********************************************************************//**
67  * @brief Void constructor
68  *
69  * Constructs empty point spread function.
70  ***************************************************************************/
72 {
73  // Initialise class members
74  init_members();
75 
76  // Return
77  return;
78 }
79 
80 
81 /***********************************************************************//**
82  * @brief File constructor
83  *
84  * @param[in] filename PSF FITS file.
85  *
86  * Constructs point spread function from a FITS file.
87  ***************************************************************************/
89 {
90  // Initialise class members
91  init_members();
92 
93  // Load point spread function from file
94  load(filename);
95 
96  // Return
97  return;
98 }
99 
100 
101 /***********************************************************************//**
102  * @brief Copy constructor
103  *
104  * @param[in] psf Point spread function.
105  *
106  * Constructs point spread function by copying from another point spread
107  * function.
108  ***************************************************************************/
110 {
111  // Initialise class members
112  init_members();
113 
114  // Copy members
115  copy_members(psf);
116 
117  // Return
118  return;
119 }
120 
121 
122 /***********************************************************************//**
123  * @brief Destructor
124  *
125  * Destructs point spread function.
126  ***************************************************************************/
128 {
129  // Free members
130  free_members();
131 
132  // Return
133  return;
134 }
135 
136 
137 /*==========================================================================
138  = =
139  = Operators =
140  = =
141  ==========================================================================*/
142 
143 /***********************************************************************//**
144  * @brief Assignment operator
145  *
146  * @param[in] psf Point spread function.
147  * @return Point spread function.
148  *
149  * Assigns point spread function.
150  ***************************************************************************/
152 {
153  // Execute only if object is not identical
154  if (this != &psf) {
155 
156  // Copy base class members
157  this->GCTAPsf::operator=(psf);
158 
159  // Free members
160  free_members();
161 
162  // Initialise private members
163  init_members();
164 
165  // Copy members
166  copy_members(psf);
167 
168  } // endif: object was not identical
169 
170  // Return this object
171  return *this;
172 }
173 
174 
175 /***********************************************************************//**
176  * @brief Return point spread function (in units of sr^-1)
177  *
178  * @param[in] delta Angular separation between true and measured photon
179  * directions (rad).
180  * @param[in] logE Log10 of the true photon energy (TeV).
181  * @param[in] theta Offset angle in camera system (rad).
182  * @param[in] phi Azimuth angle in camera system (rad). Not used.
183  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
184  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
185  * @param[in] etrue Use true energy (true/false). Not used.
186  *
187  * Returns the point spread function for a given angular separation in units
188  * of sr^-1 for a given energy. If the King profile parameters are invalid
189  * the method returns zero.
190  ***************************************************************************/
191 double GCTAPsfKing::operator()(const double& delta,
192  const double& logE,
193  const double& theta,
194  const double& phi,
195  const double& zenith,
196  const double& azimuth,
197  const bool& etrue) const
198 {
199  // Initialise PSF value
200  double psf = 0.0;
201 
202  // Update the parameter cache
203  update(logE, theta);
204 
205  // Continue only if delta is smaller than PSF radius and normalization is
206  // positive
207  if ((delta <= m_par_rmax) && (m_par_norm > 0.0)) {
208 
209  // Compute PSF value
210  double arg = delta / m_par_sigma;
211  double arg2 = arg * arg;
212  psf = m_par_norm *
213  std::pow((1.0 + 1.0 / (2.0 * m_par_gamma) * arg2), -m_par_gamma);
214 
215  // If we are at large offset angles, add a smooth ramp down to
216  // avoid steps in the log-likelihood computation
217  #if defined(G_SMOOTH_PSF)
218  double ramp_down = 0.95 * m_par_rmax;
219  double norm_down = 1.0 / (m_par_rmax - ramp_down);
220  if (delta > ramp_down) {
221  double x = norm_down * (delta - ramp_down);
222  psf *= 1.0 - x * x;
223  }
224  #endif
225 
226  } // endif: PSF radius was smaller than PSF radius and normalization was > 0
227 
228  // Return PSF
229  return psf;
230 }
231 
232 
233 /*==========================================================================
234  = =
235  = Public methods =
236  = =
237  ==========================================================================*/
238 
239 /***********************************************************************//**
240  * @brief Clear point spread function
241  *
242  * Clears point spread function.
243  ***************************************************************************/
245 {
246  // Free class members (base and derived classes, derived class first)
247  free_members();
248  this->GCTAPsf::free_members();
249 
250  // Initialise members
251  this->GCTAPsf::init_members();
252  init_members();
253 
254  // Return
255  return;
256 }
257 
258 
259 /***********************************************************************//**
260  * @brief Clone point spread functions
261  *
262  * @return Deep copy of point spread function.
263  *
264  * Returns a pointer to a deep copy of the point spread function.
265  ***************************************************************************/
267 {
268  return new GCTAPsfKing(*this);
269 }
270 
271 
272 /***********************************************************************//**
273  * @brief Read point spread function from FITS table
274  *
275  * @param[in] table FITS table.
276  *
277  * @exception GException::invalid_value
278  * Response table is not two-dimensional.
279  *
280  * Reads the point spread function form the FITS @p table. The following
281  * column names are mandatory:
282  *
283  * ENERG_LO - Energy lower bin boundaries
284  * ENERG_HI - Energy upper bin boundaries
285  * THETA_LO - Offset angle lower bin boundaries
286  * THETA_HI - Offset angle upper bin boundaries
287  * GAMMA - Gamma
288  * SIGMA - Sigma
289  *
290  * The data are stored in the m_psf member. The energy axis will be set to
291  * log10, the offset angle axis to radians.
292  ***************************************************************************/
293 void GCTAPsfKing::read(const GFitsTable& table)
294 {
295  // Clear response table
296  m_psf.clear();
297 
298  // Read PSF table
299  m_psf.read(table);
300 
301  // Get mandatory indices (throw exception if not found)
302  m_inx_energy = m_psf.axis("ENERG");
303  m_inx_theta = m_psf.axis("THETA");
304  m_inx_gamma = m_psf.table("GAMMA");
305  m_inx_sigma = m_psf.table("SIGMA");
306 
307  // Throw an exception if the table is not two-dimensional
308  if (m_psf.axes() != 2) {
309  std::string msg = "Expected two-dimensional point spread function "
310  "response table but found "+
312  " dimensions. Please specify a two-dimensional "
313  "point spread function.";
314  throw GException::invalid_value(G_READ, msg);
315  }
316 
317  // Set energy axis to logarithmic scale
319 
320  // Set offset angle axis to radians
322 
323  // Convert sigma parameters to radians
325 
326  // Return
327  return;
328 }
329 
330 
331 /***********************************************************************//**
332  * @brief Write point spread function into FITS table
333  *
334  * @param[in] table FITS binary table.
335  *
336  * Writes point spread function into a FITS binary @p table.
337  *
338  * @todo Add keywords.
339  ***************************************************************************/
341 {
342  // Create a copy of the response table
344 
345  // Convert sigma parameters back to degrees (only if there are tables)
346  if (psf.tables() > 0) {
348  }
349 
350  // Write response table
351  psf.write(table);
352 
353  // Return
354  return;
355 }
356 
357 
358 /***********************************************************************//**
359  * @brief Load point spread function from FITS file
360  *
361  * @param[in] filename FITS file name.
362  *
363  * Loads the point spread function from a FITS file.
364  *
365  * If no extension name is given the method scans the `HDUCLASS` keywords
366  * of all extensions and loads the background from the first extension
367  * for which `HDUCLAS4=PSF_KING`.
368  *
369  * Otherwise, the background will be loaded from the `POINT SPREAD FUNCTION`
370  * extension.
371  ***************************************************************************/
372 void GCTAPsfKing::load(const GFilename& filename)
373 {
374  // Open FITS file
375  GFits fits(filename);
376 
377  // Get the default extension name. If no GADF compliant name was found
378  // then set the default extension name to "POINT SPREAD FUNCTION".
379  std::string extname = gammalib::gadf_hduclas4(fits, "PSF_KING");
380  if (extname.empty()) {
382  }
383 
384  // Get PSF table
385  const GFitsTable& table = *fits.table(filename.extname(extname));
386 
387  // Read PSF from table
388  read(table);
389 
390  // Close FITS file
391  fits.close();
392 
393  // Store filename
395 
396  // Return
397  return;
398 }
399 
400 /***********************************************************************//**
401  * @brief Save point spread function table into FITS file
402  *
403  * @param[in] filename FITS file name.
404  * @param[in] clobber Overwrite existing file?
405  *
406  * Saves point spread function into a FITS file. If a file with the given
407  * @p filename does not yet exist it will be created, otherwise the method
408  * opens the existing file. The method will create a (or replace an existing)
409  * point spread function extension. The extension name can be specified as
410  * part of the @p filename, or if no extension name is given, is assumed to
411  * be `POINT SPREAD FUNCTION`.
412  *
413  * An existing file will only be modified if the @p clobber flag is set to
414  * true.
415  ***************************************************************************/
416 void GCTAPsfKing::save(const GFilename& filename, const bool& clobber) const
417 {
418  // Get extension name
419  std::string extname = filename.extname(gammalib::extname_cta_psfking);
420 
421  // Open or create FITS file (without extension name since the requested
422  // extension may not yet exist in the file)
423  GFits fits(filename.url(), true);
424 
425  // Remove extension if it exists already
426  if (fits.contains(extname)) {
427  fits.remove(extname);
428  }
429 
430  // Create binary table
432 
433  // Write the background table
434  write(table);
435 
436  // Set binary table extension name
437  table.extname(extname);
438 
439  // Append table to FITS file
440  fits.append(table);
441 
442  // Save to file
443  fits.save(clobber);
444 
445  // Return
446  return;
447 }
448 
449 
450 /***********************************************************************//**
451  * @brief Simulate PSF offset (radians)
452  *
453  * @param[in] ran Random number generator.
454  * @param[in] logE Log10 of the true photon energy (TeV).
455  * @param[in] theta Offset angle in camera system (rad).
456  * @param[in] phi Azimuth angle in camera system (rad). Not used.
457  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
458  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
459  * @param[in] etrue Use true energy (true/false). Not used.
460  *
461  * Draws a random offset angle for the King Profile. If the King profile
462  * parameters are invalid the method returns a zero offset angle.
463  ***************************************************************************/
464 double GCTAPsfKing::mc(GRan& ran,
465  const double& logE,
466  const double& theta,
467  const double& phi,
468  const double& zenith,
469  const double& azimuth,
470  const bool& etrue) const
471 {
472  // Initialise random offset
473  double delta = 0.0;
474 
475  // Update the parameter cache
476  update(logE, theta);
477 
478  // Continue only if normalization is positive
479  if (m_par_norm > 0.0) {
480 
481  // Compute exponent
482  double exponent = 1.0 / (1.0-m_par_gamma);
483 
484  // Sample until delta <= m_par_rmax
485  do {
486 
487  // Get uniform random number
488  double u = ran.uniform();
489 
490  // Draw random offset using inversion sampling
491  double u_max = (std::pow((1.0 - u), exponent) - 1.0) * m_par_gamma;
492  delta = m_par_sigma * std::sqrt(2.0 * u_max);
493 
494  } while (delta > m_par_rmax);
495 
496  } // endif: normalization was positive
497 
498  // Return PSF offset
499  return delta;
500 }
501 
502 
503 /***********************************************************************//**
504  * @brief Return maximum size of PSF (radians)
505  *
506  * @param[in] logE Log10 of the true photon energy (TeV).
507  * @param[in] theta Offset angle in camera system (rad).
508  * @param[in] phi Azimuth angle in camera system (rad). Not used.
509  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
510  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
511  * @param[in] etrue Use true energy (true/false). Not used.
512  *
513  * Determine the radius beyond which the PSF becomes negligible. This radius
514  * is set by this method to where the containment fraction become 99.9%. If
515  * the King profile parameters are invalid the method returns zero.
516  ***************************************************************************/
517 double GCTAPsfKing::delta_max(const double& logE,
518  const double& theta,
519  const double& phi,
520  const double& zenith,
521  const double& azimuth,
522  const bool& etrue) const
523 {
524  // Initialise PSF radius
525  double radius = 0.0;
526 
527  // Update the parameter cache
528  update(logE, theta);
529 
530  // Compute radius if normalization is positive
531  if (m_par_norm > 0.0) {
532  radius = r_max(logE, theta);
533  }
534 
535  // Return maximum PSF radius
536  return radius;
537 }
538 
539 
540 /***********************************************************************//**
541  * @brief Return the radius that contains a fraction of the events (radians)
542  *
543  * @param[in] fraction of events (0.0-1.0)
544  * @param[in] logE Log10 of the true photon energy (TeV).
545  * @param[in] theta Offset angle in camera system (rad). Not used.
546  * @param[in] phi Azimuth angle in camera system (rad). Not used.
547  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
548  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
549  * @param[in] etrue Use true energy (true/false). Not used.
550  * @return Containment radius (radians).
551  *
552  * @exception GException::invalid_argument
553  * Invalid fraction specified.
554  *
555  * Calculate the radius from the center that contains 'fraction' percent
556  * of the events. fraction * 100. = Containment %. If the King profile
557  * parameters are invalid the method returns zero.
558  ***************************************************************************/
559 double GCTAPsfKing::containment_radius(const double& fraction,
560  const double& logE,
561  const double& theta,
562  const double& phi,
563  const double& zenith,
564  const double& azimuth,
565  const bool& etrue) const
566 {
567  // Initialise containment radius
568  double radius = 0.0;
569 
570  // Check input argument
571  if (fraction <= 0.0 || fraction >= 1.0) {
572  std::string message = "Containment fraction "+
573  gammalib::str(fraction)+" must be between " +
574  "0.0 and 1.0, not inclusive.";
576  }
577 
578  // Update the parameter cache
579  update(logE, theta);
580 
581  // Continue only if normalization is positive
582  if (m_par_norm > 0.0) {
583 
584  // Use analytic calculation
585  double arg1 = std::pow(1.0 - fraction, 1.0/(1.0-m_par_gamma));
586  double arg2 = 2.0 * m_par_gamma * (arg1 - 1.0);
587  radius = (arg2 > 0) ? m_par_sigma * std::sqrt(arg2) : 0.0;
588 
589  } // endif: King profile parameters are valid
590 
591  // Return radius containing fraction of events
592  return radius;
593 }
594 
595 
596 /***********************************************************************//**
597  * @brief Print point spread function information
598  *
599  * @param[in] chatter Chattiness.
600  * @return String containing point spread function information.
601  ***************************************************************************/
602 std::string GCTAPsfKing::print(const GChatter& chatter) const
603 {
604  // Initialise result string
605  std::string result;
606 
607  // Continue only if chatter is not silent
608  if (chatter != SILENT) {
609 
610  // Append header
611  result.append("=== GCTAPsfKing ===");
612  result.append("\n"+gammalib::parformat("Filename")+m_filename);
613 
614  // Initialise information
615  int nebins = 0;
616  int nthetabins = 0;
617  double emin = 0.0;
618  double emax = 0.0;
619  double omin = 0.0;
620  double omax = 0.0;
621 
622  // Extract information if there are axes in the response table
623  if (m_psf.axes() > 0) {
624  nebins = m_psf.axis_bins(m_inx_energy);
625  nthetabins = m_psf.axis_bins(m_inx_theta);
626  emin = m_psf.axis_lo(m_inx_energy,0);
627  emax = m_psf.axis_hi(m_inx_energy,nebins-1);
628  omin = m_psf.axis_lo(m_inx_theta,0);
629  omax = m_psf.axis_hi(m_inx_theta,nthetabins-1);
630  }
631 
632  // Append information
633  result.append("\n"+gammalib::parformat("Number of energy bins") +
634  gammalib::str(nebins));
635  result.append("\n"+gammalib::parformat("Number of offset bins") +
636  gammalib::str(nthetabins));
637  result.append("\n"+gammalib::parformat("Log10(Energy) range"));
638  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
639  result.append("\n"+gammalib::parformat("Offset angle range"));
640  result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
641 
642  } // endif: chatter was not silent
643 
644  // Return result
645  return result;
646 }
647 
648 
649 /*==========================================================================
650  = =
651  = Private methods =
652  = =
653  ==========================================================================*/
654 
655 /***********************************************************************//**
656  * @brief Initialise class members
657  ***************************************************************************/
659 {
660  // Initialise members
661  m_filename.clear();
662  m_psf.clear();
663  m_inx_energy = 0;
664  m_inx_theta = 1;
665  m_inx_gamma = 0;
666  m_inx_sigma = 1;
667  m_par_logE = -1.0e30;
668  m_par_theta = -1.0;
669  m_par_norm = 0.0;
670  m_par_gamma = 0.0;
671  m_par_sigma = 0.0;
672  m_par_sigma2 = 0.0;
673  m_par_rmax = 0.0;
674 
675  // Return
676  return;
677 }
678 
679 
680 /***********************************************************************//**
681  * @brief Copy class members
682  *
683  * @param[in] psf Point spread function.
684  ***************************************************************************/
686 {
687  // Copy members
688  m_filename = psf.m_filename;
689  m_psf = psf.m_psf;
691  m_inx_theta = psf.m_inx_theta;
692  m_inx_gamma = psf.m_inx_gamma;
693  m_inx_sigma = psf.m_inx_sigma;
694  m_par_logE = psf.m_par_logE;
695  m_par_theta = psf.m_par_theta;
696  m_par_norm = psf.m_par_norm;
697  m_par_gamma = psf.m_par_gamma;
698  m_par_sigma = psf.m_par_sigma;
700  m_par_rmax = psf.m_par_rmax;
701 
702  // Return
703  return;
704 }
705 
706 
707 /***********************************************************************//**
708  * @brief Delete class members
709  ***************************************************************************/
711 {
712  // Return
713  return;
714 }
715 
716 
717 /***********************************************************************//**
718  * @brief Update PSF parameter cache
719  *
720  * @param[in] logE Log10 of the true photon energy (TeV).
721  * @param[in] theta Offset angle.
722  *
723  * This method updates the PSF parameter cache. As the performance table PSF
724  * only depends on energy, the only parameter on which the cache values
725  * depend is the energy. If the PSF parameters are invalid the m_par_norm
726  * member will be set to zero. Valid PSF parameters are \f$\gamma > 1\f$ and
727  * \f$\sigma > 0\f$.
728  ***************************************************************************/
729 void GCTAPsfKing::update(const double& logE, const double& theta) const
730 {
731  // Maximum PSF radius in radians. We set the maximum radius to 2 degrees
732  // since it's hard to imaging that a PSF will ever be larger.
733  const double r_max = 2.0 * gammalib::deg2rad;
734 
735  // Only compute PSF parameters if arguments have changed
736  if (logE != m_par_logE || theta != m_par_theta) {
737 
738  // Save parameters
739  m_par_logE = logE;
740  m_par_theta = theta;
741 
742  // Determine sigma and gamma by interpolating between nodes
743  std::vector<double> pars = (m_inx_energy == 0) ? m_psf(logE, theta) :
744  m_psf(theta, logE);
745 
746  // Set parameters
747  m_par_gamma = pars[m_inx_gamma];
748  m_par_sigma = pars[m_inx_sigma];
750 
751  // Check for parameter sanity. The King function is not defined for
752  // gamma values equal or smaller than one and non-positive sigma
753  // values. If this is the case we set the normalization and the
754  // maximum radius to zero ...
755  if (m_par_gamma <= 1.0 || m_par_sigma <= 0.0) {
756  m_par_norm = 0.0;
757  m_par_rmax = 0.0;
758  }
759 
760  // ... otherwise we compute the normalization
761  else {
762 
763  // Determine normalisation for given parameters
764  m_par_norm = 1.0 / gammalib::twopi * (1.0 - 1.0 / m_par_gamma) /
765  m_par_sigma2;
766 
767  // Determine maximum PSF radius
768  m_par_rmax = this->r_max(logE, theta);
769 
770  // Make sure that radius is smaller than 2 degrees
771  if (m_par_rmax > r_max) {
772 
773  // Restrict PSF radius to 2 degrees
774  m_par_rmax = r_max;
775 
776  // Correct PSF normalization for restriction
777  double u_max = (m_par_rmax*m_par_rmax) / (2.0 * m_par_sigma2);
778  double norm = 1.0 - std::pow((1.0 + u_max/m_par_gamma),
779  (1.0-m_par_gamma));
780  m_par_norm /= norm;
781 
782  } // endif: restricted PSF radius
783 
784  // If maximum PSF radius is zero then also set the normalization
785  // to zero
786  else if (m_par_rmax == 0.0) {
787  m_par_norm = 0.0;
788  }
789 
790  } // endif: King profile parameters were valid
791 
792  } // endif: PSF parameters have changed
793 
794  // Return
795  return;
796 }
797 
798 
799 /***********************************************************************//**
800  * @brief Return maximum size of PSF (radians)
801  *
802  * @param[in] logE Log10 of the true photon energy (TeV).
803  * @param[in] theta Offset angle in camera system (radians).
804  * @return Maximum size of point spread function (radians).
805  *
806  * Determine the radius beyond which the PSF becomes negligible. This radius
807  * is set by this method to where the containment fraction become 99.9%.
808  *
809  * This method requires the m_par_gamma and m_par_sigma to be set to valid
810  * values. In case that the parameters lead to an infinite result the method
811  * will return a maximum size of zero.
812  ***************************************************************************/
813 double GCTAPsfKing::r_max(const double& logE,
814  const double& theta) const
815 {
816  // Set 99.9% containment fraction
817  const double F = 0.999;
818 
819  // Compute maximum PSF radius for containment fraction
820  double u_max = (std::pow((1.0 - F), (1.0/(1.0-m_par_gamma))) - 1.0) *
821  m_par_gamma;
822  double radius = (gammalib::is_infinite(u_max)) ? 0.0 :
823  m_par_sigma * std::sqrt(2.0 * u_max);
824 
825  // Return maximum PSF radius
826  return radius;
827 }
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)
void init_members(void)
Initialise class members.
GCTAPsfKing(void)
Void constructor.
Definition: GCTAPsfKing.cpp:71
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
Random number generator class definition.
GCTAPsfKing * clone(void) const
Clone point spread functions.
const GCTAResponseTable & table(void) const
Return response table.
void load(const GFilename &filename)
Load point spread function from FITS file.
double m_par_rmax
Maximum PSF radius.
Abstract base class for the CTA point spread function.
Definition: GCTAPsf.hpp:47
CTA point spread function class with a King profile.
Definition: GCTAPsfKing.hpp:54
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
Random number generator class.
Definition: GRan.hpp:44
const int & tables(void) const
Return number of tables.
int m_inx_gamma
Gamma.
void scale(const int &table, const double &scale)
Scale table.
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.
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)
void read(const GFitsTable &table)
Read response table from FITS table HDU.
void read(const GFitsTable &table)
Read point spread function from FITS table.
void init_members(void)
Initialise class members.
Definition: GCTAPsf.cpp:142
GFilename filename(void) const
Return filename.
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)
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
GCTAResponseTable m_psf
PSF response table.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
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.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
void clear(void)
Clear point spread function.
const std::string extname_cta_psfking
Definition: GCTAPsfKing.hpp:42
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
void clear(void)
Clear response table.
Filename class.
Definition: GFilename.hpp:62
#define G_READ
Definition: GCTAPsfKing.cpp:45
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
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)
double r_max(const double &logE, const double &theta) const
Return maximum size of PSF (radians)
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void copy_members(const GCTAPsfKing &psf)
Copy class members.
GChatter
Definition: GTypemaps.hpp:33
void write(GFitsBinTable &table) const
Write point spread function into FITS table.
GCTAPsfKing & operator=(const GCTAPsfKing &psf)
Assignment operator.
GFilename m_filename
Name of Aeff response file.
void save(const GFilename &filename, const bool &clobber=false) const
Save point spread function table into FITS file.
King profile CTA point spread function class definition.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
int m_inx_energy
Energy index.
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.
CTA exception handler interface definition.
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.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
double m_par_sigma
King profile sigma (radians)
double m_par_sigma2
King profile sigma squared.
double m_par_norm
King profile normalization.
#define G_CONTAINMENT_RADIUS
Definition: GCTAPsfKing.cpp:46
void free_members(void)
Delete class members.
FITS binary table class.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function information.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
virtual ~GCTAPsfKing(void)
Destructor.
double m_par_theta
Cache offset angle.
Exception handler interface definition.
void update(const double &logE, const double &theta) const
Update PSF parameter cache.
FITS binary table class definition.
CTA response table class.
double m_par_gamma
King profile gamma parameter.
const double twopi
Definition: GMath.hpp:36
const double rad2deg
Definition: GMath.hpp:44
double m_par_logE
Cache energy.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
int m_inx_theta
Theta index.
const int & axes(void) const
Return number of axes of the tables.
Filename class interface definition.
int m_inx_sigma
Sigma.
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.