GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAPsf2D.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAPsf2D.cpp - CTA 2D point spread function 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 GCTAPsf2D.hpp
23  * @brief CTA 2D 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 <cmath>
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GException.hpp"
35 #include "GFilename.hpp"
36 #include "GRan.hpp"
37 #include "GFits.hpp"
38 #include "GFitsBinTable.hpp"
39 #include "GCTAPsf2D.hpp"
40 #include "GCTASupport.hpp"
41 
42 /* __ Method name definitions ____________________________________________ */
43 #define G_READ "GCTAPsf2D::read(GFitsTable&)"
44 #define G_CONTAINMENT_RADIUS "GCTAPsf2D::containment_radius(double&,"\
45  " double&, double&, double&, double&, double&, bool&)"
46 
47 /* __ Macros _____________________________________________________________ */
48 
49 /* __ Coding definitions _________________________________________________ */
50 #define G_SMOOTH_PSF
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 FITS file name.
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). Defaults to 0.0.
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.
186  ***************************************************************************/
187 double GCTAPsf2D::operator()(const double& delta,
188  const double& logE,
189  const double& theta,
190  const double& phi,
191  const double& zenith,
192  const double& azimuth,
193  const bool& etrue) const
194 {
195  #if defined(G_SMOOTH_PSF)
196  // Compute offset so that PSF goes to 0 at 5 times the sigma value. This
197  // is a kluge to get a PSF that smoothly goes to zero at the edge, which
198  // prevents steps or kinks in the log-likelihood function.
199  static const double offset = std::exp(-0.5*5.0*5.0);
200  #endif
201 
202  // Initialise PSF value
203  double psf = 0.0;
204 
205  // Update the parameter cache
206  update(logE, theta);
207 
208  // Continue only if normalization is positive
209  if (m_norm > 0.0) {
210 
211  // Compute distance squared
212  double delta2 = delta * delta;
213 
214  // Compute Psf value
215  #if defined(G_SMOOTH_PSF)
216  psf = std::exp(m_width1 * delta2) - offset;
217  if (m_norm2 > 0.0) {
218  psf += (std::exp(m_width2 * delta2)-offset) * m_norm2;
219  }
220  if (m_norm3 > 0.0) {
221  psf += (std::exp(m_width3 * delta2)-offset) * m_norm3;
222  }
223  #else
224  psf = std::exp(m_width1 * delta2);
225  if (m_norm2 > 0.0) {
226  psf += std::exp(m_width2 * delta2) * m_norm2;
227  }
228  if (m_norm3 > 0.0) {
229  psf += std::exp(m_width3 * delta2) * m_norm3;
230  }
231  #endif
232  psf *= m_norm;
233 
234  #if defined(G_SMOOTH_PSF)
235  // Make sure that PSF is non-negative
236  if (psf < 0.0) {
237  psf = 0.0;
238  }
239  #endif
240 
241  } // endif: normalization was positive
242 
243  // Return PSF
244  return psf;
245 }
246 
247 
248 /*==========================================================================
249  = =
250  = Public methods =
251  = =
252  ==========================================================================*/
253 
254 /***********************************************************************//**
255  * @brief Clear point spread function
256  *
257  * Clears point spread function.
258  ***************************************************************************/
260 {
261  // Free class members
262  free_members();
263  this->GCTAPsf::free_members();
264 
265  // Initialise members
266  this->GCTAPsf::init_members();
267  init_members();
268 
269  // Return
270  return;
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Clone point spread functions
276  *
277  * @return Deep copy of point spread function.
278  *
279  * Returns a pointer to a deep copy of the point spread function.
280  ***************************************************************************/
282 {
283  return new GCTAPsf2D(*this);
284 }
285 
286 
287 /***********************************************************************//**
288  * @brief Read point spread function from FITS table
289  *
290  * @param[in] table FITS table.
291  *
292  * @exception GException::invalid_value
293  * Response table is not two-dimensional.
294  *
295  * Reads the point spread function form the FITS @p table. The following
296  * column names are mandatory:
297  *
298  * ENERG_LO - Energy lower bin boundaries
299  * ENERG_HI - Energy upper bin boundaries
300  * THETA_LO - Offset angle lower bin boundaries
301  * THETA_HI - Offset angle upper bin boundaries
302  * SIGMA_1 - 1st Gaussian sigma
303  * AMPL_2 - 2nd Gaussian relative amplitude
304  * SIGMA_2 - 2nd Gaussian sigma
305  * AMPL_3 - 3rd Gaussian relative amplitude
306  * SIGMA_3 - 3rd Gaussian sigma
307  *
308  * The data are stored in the m_psf member. The energy axis will be set to
309  * log10, the offset angle axis to radians.
310  ***************************************************************************/
311 void GCTAPsf2D::read(const GFitsTable& table)
312 {
313  // Clear response table
314  m_psf.clear();
315 
316  // Read PSF table
317  m_psf.read(table);
318 
319  // Get mandatory indices (throw exception if not found)
320  m_inx_energy = m_psf.axis("ENERG");
321  m_inx_theta = m_psf.axis("THETA");
322  m_inx_sigma1 = m_psf.table("SIGMA_1");
323  m_inx_ampl2 = m_psf.table("AMPL_2");
324  m_inx_sigma2 = m_psf.table("SIGMA_2");
325  m_inx_ampl3 = m_psf.table("AMPL_3");
326  m_inx_sigma3 = m_psf.table("SIGMA_3");
327 
328  // Throw an exception if the table is not two-dimensional
329  if (m_psf.axes() != 2) {
330  std::string msg = "Expected two-dimensional point spread function "
331  "response table but found "+
333  " dimensions. Please specify a two-dimensional "
334  "point spread function.";
335  throw GException::invalid_value(G_READ, msg);
336  }
337 
338  // Set energy axis to logarithmic scale
340 
341  // Set offset angle axis to radians
343 
344  // Convert sigma parameters to radians
348 
349  // Return
350  return;
351 }
352 
353 
354 /***********************************************************************//**
355  * @brief Write point spread function into FITS binary table
356  *
357  * @param[in] table FITS binary table.
358  *
359  * Writes point spread function into FITS binary @p table.
360  *
361  * @todo Add necessary keywords.
362  ***************************************************************************/
363 void GCTAPsf2D::write(GFitsBinTable& table) const
364 {
365  // Create a copy of the response table
367 
368  // Convert sigma parameters back to degrees (only if there are tables)
369  if (psf.tables() > 0) {
373  }
374 
375  // Write response table
376  psf.write(table);
377 
378  // Return
379  return;
380 }
381 
382 
383 /***********************************************************************//**
384  * @brief Load point spread function from FITS file
385  *
386  * @param[in] filename FITS file name.
387  *
388  * Loads the point spread function from a FITS file.
389  *
390  * If no extension name is given the method scans the `HDUCLASS` keywords
391  * of all extensions and loads the background from the first extension
392  * for which `HDUCLAS4=PSF_3GAUSS`.
393  *
394  * Otherwise, the background will be loaded from the `POINT SPREAD FUNCTION`
395  * extension.
396  ***************************************************************************/
397 void GCTAPsf2D::load(const GFilename& filename)
398 {
399  // Open FITS file
400  GFits fits(filename);
401 
402  // Get the default extension name. If no GADF compliant name was found
403  // then set the default extension name to "POINT SPREAD FUNCTION".
404  std::string extname = gammalib::gadf_hduclas4(fits, "PSF_3GAUSS");
405  if (extname.empty()) {
406  extname = gammalib::extname_cta_psf2d;
407  }
408 
409  // Get PSF table
410  const GFitsTable& table = *fits.table(filename.extname(extname));
411 
412  // Read PSF from table
413  read(table);
414 
415  // Close FITS file
416  fits.close();
417 
418  // Store filename
420 
421  // Return
422  return;
423 }
424 
425 
426 /***********************************************************************//**
427  * @brief Save point spread function table into FITS file
428  *
429  * @param[in] filename FITS file name.
430  * @param[in] clobber Overwrite existing file?
431  *
432  * Save the point spread function into a FITS file.
433  *
434  * If no extension name is provided, the point spread function will be saved
435  * into the `POINT SPREAD FUNCTION` extension.
436  ***************************************************************************/
437 void GCTAPsf2D::save(const GFilename& filename, const bool& clobber) const
438 {
439  // Get extension name
440  std::string extname = filename.extname(gammalib::extname_cta_psf2d);
441 
442  // Open or create FITS file (without extension name since the requested
443  // extension may not yet exist in the file)
444  GFits fits(filename.url(), true);
445 
446  // Remove extension if it exists already
447  if (fits.contains(extname)) {
448  fits.remove(extname);
449  }
450 
451  // Create binary table
453 
454  // Write the background table
455  write(table);
456 
457  // Set binary table extension name
458  table.extname(extname);
459 
460  // Append table to FITS file
461  fits.append(table);
462 
463  // Save to file
464  fits.save(clobber);
465 
466  // Return
467  return;
468 }
469 
470 
471 /***********************************************************************//**
472  * @brief Simulate PSF offset (radians)
473  *
474  * @param[in] ran Random number generator.
475  * @param[in] logE Log10 of the true photon energy (TeV).
476  * @param[in] theta Offset angle in camera system (rad).
477  * @param[in] phi Azimuth angle in camera system (rad). Not used.
478  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
479  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
480  * @param[in] etrue Use true energy (true/false). Not used.
481  *
482  * Draws a random offset for a three component Gaussian PSF.
483  ***************************************************************************/
484 double GCTAPsf2D::mc(GRan& ran,
485  const double& logE,
486  const double& theta,
487  const double& phi,
488  const double& zenith,
489  const double& azimuth,
490  const bool& etrue) const
491 {
492  // Update the parameter cache
493  update(logE, theta);
494 
495  // Select in which Gaussian we are
496  double sigma = m_sigma1;
497  double sum1 = m_sigma1;
498  double sum2 = m_sigma2 * m_norm2;
499  double sum3 = m_sigma3 * m_norm3;
500  double sum = sum1 + sum2 + sum3;
501  double u = ran.uniform() * sum;
502  if (sum2 > 0.0 && u >= sum2) {
503  sigma = m_sigma3;
504  }
505  else if (sum1 > 0.0 && u >= sum1) {
506  sigma = m_sigma2;
507  }
508 
509  // Now draw from the selected Gaussian
510  double delta = sigma * ran.chisq2();
511 
512  // Return PSF offset
513  return delta;
514 }
515 
516 
517 /***********************************************************************//**
518  * @brief Return maximum size of PSF (radians)
519  *
520  * @param[in] logE Log10 of the true photon energy (TeV).
521  * @param[in] theta Offset angle in camera system (rad). Not used.
522  * @param[in] phi Azimuth angle in camera system (rad). Not used.
523  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
524  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
525  * @param[in] etrue Use true energy (true/false). Not used.
526  *
527  * Determine the radius beyond which the PSF becomes negligible. This radius
528  * is set by this method to \f$5 \times \sigma\f$, where \f$\sigma\f$ is the
529  * Gaussian width of the largest PSF component.
530  ***************************************************************************/
531 double GCTAPsf2D::delta_max(const double& logE,
532  const double& theta,
533  const double& phi,
534  const double& zenith,
535  const double& azimuth,
536  const bool& etrue) const
537 {
538  // Update the parameter cache
539  update(logE, theta);
540 
541  // Compute maximum sigma
542  double sigma = m_sigma1;
543  if (m_sigma2 > sigma) sigma = m_sigma2;
544  if (m_sigma3 > sigma) sigma = m_sigma3;
545 
546  // Compute maximum PSF radius
547  double radius = 5.0 * sigma;
548 
549  // Return maximum PSF radius
550  return radius;
551 }
552 
553 
554 /***********************************************************************//**
555  * @brief Return the radius that contains a fraction of the events (radians)
556  *
557  * @param[in] fraction of events (0.0 < fraction < 1.0)
558  * @param[in] logE Log10 of the true photon energy (TeV).
559  * @param[in] theta Offset angle in camera system (rad).
560  * @param[in] phi Azimuth angle in camera system (rad). Not used.
561  * @param[in] zenith Zenith angle in Earth system (rad). Not used.
562  * @param[in] azimuth Azimuth angle in Earth system (rad). Not used.
563  * @param[in] etrue Use true energy (true/false). Not used.
564  * @return Containment radius (radians).
565  *
566  * @exception GException::invalid_argument
567  * Invalid fraction specified.
568  *
569  * Uses the Newton-Raphson method to find which 'a' solves the equation:
570  *
571  * \f[
572  *
573  * fraction = m\_norm * \left( e^{m\_width1 * a^2} +
574  * m\_norm2 * e^{m\_width2 * a^2} +
575  * m\_norm3 * e^{m\_width3 * a^2} \right)
576  *
577  * \f]
578  *
579  * Calculate the radius from the center that contains @p fraction of the
580  * events (@p fraction * 100. = Containment % ). Fraction must be
581  *
582  * \f[
583  * 0.0 < fraction < 1.0
584  * \f]
585  ***************************************************************************/
586 double GCTAPsf2D::containment_radius(const double& fraction,
587  const double& logE,
588  const double& theta,
589  const double& phi,
590  const double& zenith,
591  const double& azimuth,
592  const bool& etrue) const
593 {
594  // Set maximum number of Newton-Raphson loops before giving up and
595  // required accuracy
596  const int itermax = 20;
597  const double convergence = 1.0e-6;
598 
599  // Check input argument
600  if (fraction <= 0.0 || fraction >= 1.0) {
601  std::string message = "Containment fraction "+
602  gammalib::str(fraction)+" must be between " +
603  "0.0 and 1.0, not inclusive.";
605  }
606 
607  // Update the parameter cache
608  update(logE, theta);
609 
610  // Set minimum and maximum radius (avoid starting at zero radius)
611  double amax = delta_max(logE, theta, phi, zenith, azimuth, etrue);
612  double amin = 0.000001 * amax;
613 
614  // Set initial radius to start Newton's method
615  double a = 0.20 * amax;
616 
617  // Set normalisation constants
618  double norm1 = (m_width1 != 0.0) ? -gammalib::pi/m_width1 : 0.0;
619  double norm2 = (m_width2 != 0.0) ? -gammalib::pi/m_width2 * m_norm2 : 0.0;
620  double norm3 = (m_width3 != 0.0) ? -gammalib::pi/m_width3 * m_norm3 : 0.0;
621 
622  // Do the Newton-Raphson loops
623  int iter = 0;
624  for (; iter < itermax; ++iter) {
625 
626  // Compute square of test radius
627  double a2 = a*a;
628 
629  // Precompute exponentials
630  double exp1 = std::exp(m_width1 * a2);
631  double exp2 = std::exp(m_width2 * a2);
632  double exp3 = std::exp(m_width3 * a2);
633 
634  // Calculate f(a)
635  double fa = 0.0;
636  fa += norm1 * (1.0 - exp1);
637  fa += norm2 * (1.0 - exp2);
638  fa += norm3 * (1.0 - exp3);
639  fa *= m_norm;
640  fa -= fraction;
641 
642  // If we've met the desired accuracy then exit the loop
643  if (std::abs(fa) < convergence) {
644  break;
645  }
646 
647  // Calculate f'(a)
648  double fp = 0.0;
649  fp += norm1 * exp1 * m_width1;
650  fp += norm2 * exp2 * m_width2;
651  fp += norm3 * exp3 * m_width3;
652  fp *= -2.0 * a * m_norm;
653 
654  // Calculate next point via x+1 = x - f(a) / f'(a)
655  if (fp != 0.0) {
656  a -= fa / fp;
657  if (a < amin) {
658  a = amin;
659  }
660  if (a > amax) {
661  a = amax;
662  }
663  }
664  else {
665  break;
666  }
667 
668  } // endfor: Newton-Raphson loops
669 
670  // Warn the user if we didn't converge
671  if (iter == itermax-1) {
672  std::string msg = "Unable to converge within " +
673  gammalib::str(convergence)+" of the root in less "
674  "than "+gammalib::str(itermax)+" iterations." ;
676  }
677 
678  // Return containment radius
679  return a;
680 }
681 
682 
683 /***********************************************************************//**
684  * @brief Print point spread function information
685  *
686  * @param[in] chatter Chattiness.
687  * @return String containing point spread function information.
688  ***************************************************************************/
689 std::string GCTAPsf2D::print(const GChatter& chatter) const
690 {
691  // Initialise result string
692  std::string result;
693 
694  // Continue only if chatter is not silent
695  if (chatter != SILENT) {
696 
697  // Append header
698  result.append("=== GCTAPsf2D ===");
699  result.append("\n"+gammalib::parformat("Filename")+m_filename);
700 
701  // Initialise information
702  int nebins = 0;
703  int nthetabins = 0;
704  double emin = 0.0;
705  double emax = 0.0;
706  double omin = 0.0;
707  double omax = 0.0;
708 
709  // Extract information if there are axes in the response table
710  if (m_psf.axes() > 0) {
711  nebins = m_psf.axis_bins(m_inx_energy);
712  nthetabins = m_psf.axis_bins(m_inx_theta);
713  emin = m_psf.axis_lo(m_inx_energy,0);
714  emax = m_psf.axis_hi(m_inx_energy,nebins-1);
715  omin = m_psf.axis_lo(m_inx_theta,0);
716  omax = m_psf.axis_hi(m_inx_theta,nthetabins-1);
717  }
718 
719  // Append information
720  result.append("\n"+gammalib::parformat("Number of energy bins") +
721  gammalib::str(nebins));
722  result.append("\n"+gammalib::parformat("Number of offset bins") +
723  gammalib::str(nthetabins));
724  result.append("\n"+gammalib::parformat("Energy range"));
725  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
726  result.append("\n"+gammalib::parformat("Offset angle range"));
727  result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
728 
729  } // endif: chatter was not silent
730 
731  // Return result
732  return result;
733 }
734 
735 
736 /*==========================================================================
737  = =
738  = Private methods =
739  = =
740  ==========================================================================*/
741 
742 /***********************************************************************//**
743  * @brief Initialise class members
744  ***************************************************************************/
746 {
747  // Initialise members
748  m_filename.clear();
749  m_psf.clear();
750  m_inx_energy = 0;
751  m_inx_theta = 1;
752  m_inx_sigma1 = 0;
753  m_inx_ampl2 = 1;
754  m_inx_sigma2 = 2;
755  m_inx_ampl3 = 3;
756  m_inx_sigma3 = 4;
757  m_par_logE = -1.0e30;
758  m_par_theta = -1.0;
759  m_norm = 1.0;
760  m_norm2 = 0.0;
761  m_norm3 = 0.0;
762  m_sigma1 = 0.0;
763  m_sigma2 = 0.0;
764  m_sigma3 = 0.0;
765  m_width1 = 0.0;
766  m_width2 = 0.0;
767  m_width3 = 0.0;
768 
769  // Return
770  return;
771 }
772 
773 
774 /***********************************************************************//**
775  * @brief Copy class members
776  *
777  * @param[in] psf Point spread function.
778  ***************************************************************************/
780 {
781  // Copy members
782  m_filename = psf.m_filename;
783  m_psf = psf.m_psf;
785  m_inx_theta = psf.m_inx_theta;
787  m_inx_ampl2 = psf.m_inx_ampl2;
789  m_inx_ampl3 = psf.m_inx_ampl3;
791  m_par_logE = psf.m_par_logE;
792  m_par_theta = psf.m_par_theta;
793  m_norm = psf.m_norm;
794  m_norm2 = psf.m_norm2;
795  m_norm3 = psf.m_norm3;
796  m_sigma1 = psf.m_sigma1;
797  m_sigma2 = psf.m_sigma2;
798  m_sigma3 = psf.m_sigma3;
799  m_width1 = psf.m_width1;
800  m_width2 = psf.m_width2;
801  m_width3 = psf.m_width3;
802 
803  // Return
804  return;
805 }
806 
807 
808 /***********************************************************************//**
809  * @brief Delete class members
810  ***************************************************************************/
812 {
813  // Return
814  return;
815 }
816 
817 
818 /***********************************************************************//**
819  * @brief Update PSF parameter cache
820  *
821  * @param[in] logE Log10 of the true photon energy (TeV).
822  * @param[in] theta Offset angle in camera system (rad).
823  *
824  * This method updates the PSF parameter cache.
825  ***************************************************************************/
826 void GCTAPsf2D::update(const double& logE, const double& theta) const
827 {
828  // Only compute PSF parameters if arguments have changed
829  if (logE != m_par_logE || theta != m_par_theta) {
830 
831  // Save parameters
832  m_par_logE = logE;
833  m_par_theta = theta;
834 
835  // Interpolate response parameters
836  std::vector<double> pars = (m_inx_energy == 0) ? m_psf(logE, theta) :
837  m_psf(theta, logE);
838 
839  // Set Gaussian sigmas
840  m_sigma1 = pars[m_inx_sigma1];
841  m_sigma2 = pars[m_inx_sigma2];
842  m_sigma3 = pars[m_inx_sigma3];
843 
844  // Set width parameters
845  double sigma1 = m_sigma1 * m_sigma1;
846  double sigma2 = m_sigma2 * m_sigma2;
847  double sigma3 = m_sigma3 * m_sigma3;
848 
849  // Compute Gaussian 1
850  if (sigma1 > 0.0) {
851  m_width1 = -0.5 / sigma1;
852  }
853  else {
854  m_width1 = 0.0;
855  }
856 
857  // Compute Gaussian 2
858  if (sigma2 > 0.0) {
859  m_width2 = -0.5 / sigma2;
860  m_norm2 = pars[m_inx_ampl2];
861  }
862  else {
863  m_width2 = 0.0;
864  m_norm2 = 0.0;
865  }
866 
867  // Compute Gaussian 3
868  if (sigma3 > 0.0) {
869  m_width3 = -0.5 / sigma3;
870  m_norm3 = pars[m_inx_ampl3];
871  }
872  else {
873  m_width3 = 0.0;
874  m_norm3 = 0.0;
875  }
876 
877  // Compute global normalization parameter
878  double integral = gammalib::twopi *
879  (sigma1 + sigma2*m_norm2 + sigma3*m_norm3);
880  m_norm = (integral > 0.0) ? 1.0 / integral : 0.0;
881 
882  }
883 
884  // Return
885  return;
886 }
void free_members(void)
Delete class members.
Definition: GCTAPsf.cpp:164
#define G_CONTAINMENT_RADIUS
Definition: GCTAPsf2D.cpp:44
double chisq2(void)
Returns Chi2 deviates for 2 degrees of freedom.
Definition: GRan.cpp:370
virtual ~GCTAPsf2D(void)
Destructor.
Definition: GCTAPsf2D.cpp:124
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
GCTAResponseTable m_psf
PSF response table.
Definition: GCTAPsf2D.hpp:119
GCTAPsf2D & operator=(const GCTAPsf2D &psf)
Assignment operator.
Definition: GCTAPsf2D.cpp:148
Random number generator class definition.
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1386
GFilename filename(void) const
Return filename.
Definition: GCTAPsf2D.hpp:161
const double pi
Definition: GMath.hpp:35
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
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)
Definition: GCTAPsf2D.cpp:187
void read(const GFitsTable &table)
Read point spread function from FITS table.
Definition: GCTAPsf2D.cpp:311
GCTAPsf2D(void)
Void constructor.
Definition: GCTAPsf2D.cpp:68
Abstract base class for the CTA point spread function.
Definition: GCTAPsf.hpp:47
int m_inx_energy
Energy index.
Definition: GCTAPsf2D.hpp:120
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:944
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
double m_width1
Gaussian 1 width.
Definition: GCTAPsf2D.hpp:137
Random number generator class.
Definition: GRan.hpp:44
double m_norm
Global normalization.
Definition: GCTAPsf2D.hpp:131
const int & tables(void) const
Return number of tables.
double m_sigma2
Gaussian 2 sigma.
Definition: GCTAPsf2D.hpp:135
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 m_par_theta
Cache offset angle.
Definition: GCTAPsf2D.hpp:130
int m_inx_ampl3
3nd Gaussian relative amplitude
Definition: GCTAPsf2D.hpp:125
void read(const GFitsTable &table)
Read response table from FITS table HDU.
CTA 2D point spread function class.
Definition: GCTAPsf2D.hpp:54
void init_members(void)
Initialise class members.
Definition: GCTAPsf.cpp:142
void copy_members(const GCTAPsf2D &psf)
Copy class members.
Definition: GCTAPsf2D.cpp:779
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
Definition of support function used by CTA classes.
const GCTAResponseTable & table(void) const
Return response table.
Definition: GCTAPsf2D.hpp:172
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 update(const double &logE, const double &theta) const
Update PSF parameter cache.
Definition: GCTAPsf2D.cpp:826
int m_inx_ampl2
2nd Gaussian relative amplitude
Definition: GCTAPsf2D.hpp:123
#define G_READ
Definition: GCTAPsf2D.cpp:43
int m_inx_sigma3
3nd Gaussian sigma
Definition: GCTAPsf2D.hpp:126
const double deg2rad
Definition: GMath.hpp:43
const std::string extname_cta_psf2d
Definition: GCTAPsf2D.hpp:42
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 m_width2
Gaussian 2 width.
Definition: GCTAPsf2D.hpp:138
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function information.
Definition: GCTAPsf2D.cpp:689
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)
Definition: GCTAPsf2D.cpp:484
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GChatter
Definition: GTypemaps.hpp:33
double m_norm3
Gaussian 3 normalization.
Definition: GCTAPsf2D.hpp:133
double m_sigma3
Gaussian 3 sigma.
Definition: GCTAPsf2D.hpp:136
CTA 2D point spread function class 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)
Definition: GCTAPsf2D.cpp:531
void save(const GFilename &filename, const bool &clobber=false) const
Save point spread function table into FITS file.
Definition: GCTAPsf2D.cpp:437
void load(const GFilename &filename)
Load point spread function from FITS file.
Definition: GCTAPsf2D.cpp:397
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.
void init_members(void)
Initialise class members.
Definition: GCTAPsf2D.cpp:745
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
int m_inx_sigma1
1st Gaussian sigma
Definition: GCTAPsf2D.hpp:122
int m_inx_sigma2
2nd Gaussian sigma
Definition: GCTAPsf2D.hpp:124
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
void clear(void)
Clear point spread function.
Definition: GCTAPsf2D.cpp:259
void write(GFitsBinTable &table) const
Write point spread function into FITS binary table.
Definition: GCTAPsf2D.cpp:363
FITS binary table class.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
int m_inx_theta
Theta index.
Definition: GCTAPsf2D.hpp:121
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)
Definition: GCTAPsf2D.cpp:586
double m_par_logE
Cache energy.
Definition: GCTAPsf2D.hpp:129
Exception handler interface definition.
void free_members(void)
Delete class members.
Definition: GCTAPsf2D.cpp:811
FITS binary table class definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
CTA response table class.
GCTAPsf2D * clone(void) const
Clone point spread functions.
Definition: GCTAPsf2D.cpp:281
const double twopi
Definition: GMath.hpp:36
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
double m_sigma1
Gaussian 1 sigma.
Definition: GCTAPsf2D.hpp:134
double m_width3
Gaussian 3 width.
Definition: GCTAPsf2D.hpp:139
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
const int & axes(void) const
Return number of axes of the tables.
GFilename m_filename
Name of Aeff response file.
Definition: GCTAPsf2D.hpp:118
Filename class interface definition.
Mathematical function definitions.
double m_norm2
Gaussian 2 normalization.
Definition: GCTAPsf2D.hpp:132
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489