GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GLATPsf.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GLATPsf.cpp - Fermi LAT point spread function *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2008-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 GLATPsf.cpp
23  * @brief Fermi LAT 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 "GException.hpp"
32 #include "GTools.hpp"
33 #include "GFilename.hpp"
34 #include "GFits.hpp"
35 #include "GFitsBinTable.hpp"
36 #include "GFitsTableFloatCol.hpp"
37 #include "GLATPsf.hpp"
38 #include "GLATPsfV1.hpp"
39 #include "GLATPsfV3.hpp"
40 
41 /* __ Method name definitions ____________________________________________ */
42 #define G_READ "GLATPsf::read(GFits&)"
43 #define G_WRITE "GLATPsf::write(GFits&)"
44 
45 /* __ Macros _____________________________________________________________ */
46 
47 /* __ Coding definitions _________________________________________________ */
48 
49 /* __ Debug definitions __________________________________________________ */
50 
51 /* __ Constants __________________________________________________________ */
52 
53 
54 /*==========================================================================
55  = =
56  = Constructors/destructors =
57  = =
58  ==========================================================================*/
59 
60 /***********************************************************************//**
61  * @brief Void constructor
62  *
63  * Constructs an empty Fermi LAT point spread function.
64  ***************************************************************************/
66 {
67  // Initialise class members
68  init_members();
69 
70  // Return
71  return;
72 }
73 
74 
75 /***********************************************************************//**
76  * @brief File constructor
77  *
78  * @param[in] filename FITS file name.
79  * @param[in] evtype Event type.
80  *
81  * Constructs a Fermi LAT point spread function by loading the point spread
82  * function for a given event type from a FITS response file.
83  ***************************************************************************/
84 GLATPsf::GLATPsf(const GFilename& filename, const std::string& evtype)
85 {
86  // Initialise class members
87  init_members();
88 
89  // Load PSF from FITS file
90  load(filename, evtype);
91 
92  // Return
93  return;
94 }
95 
96 
97 /***********************************************************************//**
98  * @brief Copy constructor
99  *
100  * @param[in] psf Point spread function.
101  *
102  * Constructs point spread function by copying point spread function from
103  * another object.
104  ***************************************************************************/
106 {
107  // Initialise class members
108  init_members();
109 
110  // Copy members
111  copy_members(psf);
112 
113  // Return
114  return;
115 }
116 
117 
118 /***********************************************************************//**
119  * @brief Destructor
120  ***************************************************************************/
122 {
123  // Free members
124  free_members();
125 
126  // Return
127  return;
128 }
129 
130 
131 /*==========================================================================
132  = =
133  = Operators =
134  = =
135  ==========================================================================*/
136 
137 /***********************************************************************//**
138  * @brief Assignment operator
139  *
140  * @param[in] psf Point spread function.
141  * @return Point spread function.
142  ***************************************************************************/
144 {
145  // Execute only if object is not identical
146  if (this != &psf) {
147 
148  // Free members
149  free_members();
150 
151  // Initialise private members
152  init_members();
153 
154  // Copy members
155  copy_members(psf);
156 
157  } // endif: object was not identical
158 
159  // Return this object
160  return *this;
161 }
162 
163 
164 /***********************************************************************//**
165  * @brief Return point spread function value
166  *
167  * @param[in] offset Offset angle (deg).
168  * @param[in] logE Log10 of the true photon energy (MeV).
169  * @param[in] ctheta Cosine of zenith angle.
170  *
171  * Returns the PSF value as function of the offset angle, the base 10
172  * logarithm of the energy, and the cosine of the zenith angle. This method
173  * calls the version dependent method.
174  *
175  * Returns 0 is no PSF has been allocated.
176  ***************************************************************************/
177 double GLATPsf::operator()(const double& offset, const double& logE,
178  const double& ctheta)
179 {
180  // Get PSF value
181  double psf = (m_psf != NULL) ? m_psf->psf(offset, logE, ctheta) : 0.0;
182 
183  // Return PSF value
184  return psf;
185 }
186 
187 
188 /*==========================================================================
189  = =
190  = Public methods =
191  = =
192  ==========================================================================*/
193 
194 /***********************************************************************//**
195  * @brief Clear instance
196  *
197  * This method properly resets the object to an initial state.
198  ***************************************************************************/
199 void GLATPsf::clear(void)
200 {
201  // Free class members
202  free_members();
203 
204  // Initialise members
205  init_members();
206 
207  // Return
208  return;
209 }
210 
211 
212 /***********************************************************************//**
213  * @brief Clone instance
214  *
215  * @return Pointer to deep copy of point spread function.
216  ***************************************************************************/
218 {
219  return new GLATPsf(*this);
220 }
221 
222 
223 /***********************************************************************//**
224  * @brief Load point spread function from FITS file
225  *
226  * @param[in] filename FITS file.
227  * @param[in] evtype Event type.
228  *
229  * Loads Fermi/LAT point spread function from FITS file.
230  ***************************************************************************/
231 void GLATPsf::load(const GFilename& filename, const std::string& evtype)
232 {
233  // Open FITS file
234  GFits fits(filename);
235 
236  // Read point spread function from file
237  read(fits, evtype);
238 
239  // Return
240  return;
241 }
242 
243 
244 /***********************************************************************//**
245  * @brief Save point spread function into FITS file
246  *
247  * @param[in] filename FITS file.
248  * @param[in] clobber Overwrite existing file?
249  *
250  * Saves Fermi/LAT point spread function into FITS file.
251  ***************************************************************************/
252 void GLATPsf::save(const GFilename& filename, const bool& clobber)
253 {
254  // Create FITS file
255  GFits fits;
256 
257  // Write point spread function into file
258  write(fits);
259 
260  // Close FITS file
261  fits.saveto(filename, clobber);
262 
263  // Return
264  return;
265 }
266 
267 
268 /***********************************************************************//**
269  * @brief Read point spread function from FITS file
270  *
271  * @param[in] fits FITS file.
272  * @param[in] evtype Event type.
273  *
274  * @exception GException::invalid_value
275  * Unsupported response version found.
276  *
277  * Reads the Fermi LAT point spread function from FITS file. The method
278  * determines the PSF version from the information found in the FITS file
279  * and allocates the proper PSF version class. It reads the PSF information
280  * from the FITS file.
281  *
282  * @todo Implement PSF version 2.
283  ***************************************************************************/
284 void GLATPsf::read(const GFits& fits, const std::string& evtype)
285 {
286  // Clear instance
287  clear();
288 
289  // Store event type
290  m_evtype = evtype;
291 
292  // Set extension names
293  std::string rpsf = "RPSF";
294  std::string scaling = "PSF_SCALING_PARAMS";
295  if (!fits.contains(rpsf)) {
296  rpsf += "_" + m_evtype;
297  }
298  if (!fits.contains(scaling)) {
299  scaling += "_" + m_evtype;
300  }
301 
302  // Get PSF parameters table
303  const GFitsTable& hdu_rpsf = *fits.table(rpsf);
304 
305  // Get PSF scaling parameters table
306  const GFitsTable& hdu_scale = *fits.table(scaling);
307 
308  // Determine PSF version (default version is version 1)
309  int version = (hdu_rpsf.has_card("PSFVER")) ? hdu_rpsf.integer("PSFVER") : 1;
310 
311  // Determine PSF type
312  bool front = true;
313  std::string detnam =
315  if (detnam == "FRONT") {
316  front = true;
317  }
318  else if (detnam == "BACK") {
319  front = false;
320  }
321  else {
322  front = false; // If we are here we should have a Pass 8 response.
323  // For Pass 8 the setting of this flag is irrelevant.
324  }
325 
326  // Allocate point spread function
327  switch (version) {
328  case 1:
329  m_psf = new GLATPsfV1;
330  break;
331  case 3:
332  m_psf = new GLATPsfV3;
333  break;
334  default:
335  std::string msg = "Unsupported response function version "+
336  gammalib::str(version)+". Please specify either a "
337  "version 1 or 3 point spread function.";
338  throw GException::invalid_value(G_READ, msg);
339  break;
340  }
341 
342  // Set PSF attributes (needed before reading scaling parameters for
343  // Pass 6 and Pass 7 response functions)
344  m_psf->front(front);
345 
346  // Read PSF scaling parameters
347  m_psf->read_scale(hdu_scale);
348 
349  // Read PSF
350  m_psf->read(hdu_rpsf);
351 
352  // Return
353  return;
354 }
355 
356 
357 /***********************************************************************//**
358  * @brief Write point spread function into FITS file
359  *
360  * @param[in] fits FITS file.
361  *
362  * Writes the version dependent point spread function into the FITS file. The
363  * method does nothing if no PSF has been allocated.
364  *
365  * @todo Implement PSF versions 2 and 3.
366  ***************************************************************************/
367 void GLATPsf::write(GFits& fits) const
368 {
369  // Continue only if PSF is valid
370  if (m_psf != NULL) {
371 
372  // Write point spread function
373  m_psf->write(fits);
374 
375  // Write scaling parameters
376  m_psf->write_scale(fits);
377 
378  }
379 
380  // Return
381  return;
382 }
383 
384 
385 /***********************************************************************//**
386  * @brief Returns size of PSF
387  *
388  * @return Size of point spread function.
389  *
390  * Returns the size of the PSF which is defined as the number of energy bins
391  * times the number of cos(theta) bins. If no PSF has been allocated, 0 is
392  * returned.
393  ***************************************************************************/
394 int GLATPsf::size(void) const
395 {
396  // Return size of PSF
397  return (nenergies()*ncostheta());
398 }
399 
400 
401 /***********************************************************************//**
402  * @brief Returns number of energy bins in PSF
403  *
404  * @return Number of energy bins.
405  *
406  * Returns the number of energy bins in PSF. If no PSF has been allocated,
407  * 0 is returned.
408  ***************************************************************************/
409 int GLATPsf::nenergies(void) const
410 {
411  // Retrieve number of energy bins
412  int nenergies = (m_psf != NULL) ? m_psf->nenergies() : 0;
413 
414  // Return number of energy bins
415  return nenergies;
416 }
417 
418 
419 /***********************************************************************//**
420  * @brief Returns number of cos(theta) bins in PSF
421  *
422  * @return Number of cos(theta) bins.
423  *
424  * Returns the number of cos(theta) bins in PSF. If no PSF has been
425  * allocated, 0 is returned.
426  ***************************************************************************/
427 int GLATPsf::ncostheta(void) const
428 {
429  // Retrieve number of cos(theta) bins
430  int ncostheta = (m_psf != NULL) ? m_psf->ncostheta() : 0;
431 
432  // Return number of cos(theta) bins
433  return ncostheta;
434 }
435 
436 
437 /***********************************************************************//**
438  * @brief Returns minimum cos(theta) angle
439  *
440  * @return Minimum cos(theta) angle.
441  *
442  * Returns the minimum cos(theta) angle for point spread function access. If
443  * no PSF has been allocated, 0 is returned.
444  ***************************************************************************/
445 double GLATPsf::costhetamin(void) const
446 {
447  // Retrieve number of energy bins
448  double costhetamin = (m_psf != NULL) ? m_psf->costhetamin() : 0.0;
449 
450  // Return minimum cos(theta)
451  return costhetamin;
452 }
453 
454 
455 /***********************************************************************//**
456  * @brief Set minimum cos(theta) angle for point spread function access
457  *
458  * @param[in] ctheta Cosine of maximum zenith angle.
459  *
460  * Sets the minimum cos(theta) angle. No verification of the specified
461  * maximum cos(theta) value is done.
462  ***************************************************************************/
463 void GLATPsf::costhetamin(const double& ctheta)
464 {
465  // Continue only if PSF is valid
466  if (m_psf != NULL) {
467 
468  // Set minimum cos(theta) value
469  m_psf->costhetamin(ctheta);
470 
471  }
472 
473  // Return
474  return;
475 }
476 
477 
478 /***********************************************************************//**
479  * @brief Returns PSF version
480  *
481  * Returns the PSF version number.
482  *
483  * Returns 0 if no PSF has been allocated.
484  ***************************************************************************/
485 int GLATPsf::version(void) const
486 {
487  // Retrieve version number
488  int version = (m_psf != NULL) ? m_psf->version() : 0;
489 
490  // Return version number
491  return version;
492 }
493 
494 
495 /***********************************************************************//**
496  * @brief Print point spread function information
497  *
498  * @param[in] chatter Chattiness.
499  * @return String containing point spread function information.
500  ***************************************************************************/
501 std::string GLATPsf::print(const GChatter& chatter) const
502 {
503  // Initialise result string
504  std::string result;
505 
506  // Continue only if chatter is not silent
507  if (chatter != SILENT) {
508 
509  // Append header
510  result.append("=== GLATPsf ===");
511 
512  // No PSF has been loaded ...
513  if (m_psf == NULL) {
514  result.append("\n"+gammalib::parformat("Version"));
515  result.append("No PSF loaded");
516  }
517 
518  // ... PSF has been loaded
519  else {
520  result.append("\n"+gammalib::parformat("Version"));
521  result.append(gammalib::str(version()));
522  result.append("\n"+gammalib::parformat("Event type"));
523  result.append(evtype());
524  result.append("\n"+gammalib::parformat("Energy scaling"));
525  result.append("sqrt(");
526  result.append("("+gammalib::str(m_psf->scale_par1())+"*(E/100)^");
527  result.append(gammalib::str(m_psf->scale_index())+")^2");
528  result.append(" + ("+gammalib::str(m_psf->scale_par2())+")^2)");
529  }
530 
531  } // endif: chatter was not silent
532 
533  // Return result
534  return result;
535 }
536 
537 
538 /*==========================================================================
539  = =
540  = Private methods =
541  = =
542  ==========================================================================*/
543 
544 /***********************************************************************//**
545  * @brief Initialise class members
546  ***************************************************************************/
548 {
549  // Initialise members
550  m_evtype.clear();
551  m_psf = NULL;
552 
553  // Return
554  return;
555 }
556 
557 
558 /***********************************************************************//**
559  * @brief Copy class members
560  *
561  * @param[in] psf Point spread function.
562  ***************************************************************************/
564 {
565  // Copy members
566  m_evtype = psf.m_evtype;
567 
568  // Clone members
569  if (psf.m_psf != NULL) {
570  m_psf = psf.m_psf->clone();
571  }
572  else {
573  m_psf = NULL;
574  }
575 
576  // Return
577  return;
578 }
579 
580 
581 /***********************************************************************//**
582  * @brief Delete class members
583  ***************************************************************************/
585 {
586  // Free PSF
587  if (m_psf != NULL) delete m_psf;
588 
589  // Signal that PSF is free
590  m_psf = NULL;
591 
592  // Return
593  return;
594 }
double operator()(const double &offset, const double &logE, const double &ctheta)
Return point spread function value.
Definition: GLATPsf.cpp:177
void read_scale(const GFitsTable &hdu)
Read PSF scale factors from FITS table.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
virtual GLATPsfBase * clone(void) const =0
Clones object.
virtual double psf(const double &offset, const double &logE, const double &ctheta)=0
void clear(void)
Clear instance.
Definition: GLATPsf.cpp:199
GLATPsf * clone(void) const
Clone instance.
Definition: GLATPsf.cpp:217
Fermi/LAT point spread function version 1 class definition.
const double & scale_index(void) const
Return scaling index.
Fermi/LAT point spread function version 1 class.
Definition: GLATPsfV1.hpp:47
void read(const GFits &file, const std::string &evtype)
Read point spread function from FITS file.
Definition: GLATPsf.cpp:284
const double & costhetamin(void) const
Return cosine theta minimum.
virtual void write(GFits &file) const =0
void save(const GFilename &filename, const bool &clobber=false)
Save point spread function into FITS file.
Definition: GLATPsf.cpp:252
void copy_members(const GLATPsf &psf)
Copy class members.
Definition: GLATPsf.cpp:563
const double & scale_par1(void) const
Return first scaling parameter.
int ncostheta(void) const
Return number of cosine theta bins in point spread function.
Gammalib tools definition.
FITS table float column class interface definition.
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:80
const std::string & evtype(void) const
Return event type.
Definition: GLATPsf.hpp:120
int nenergies(void) const
Returns number of energy bins in PSF.
Definition: GLATPsf.cpp:409
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
void write_scale(GFits &file) const
Write PSF scale factors.
Fermi/LAT point spread function version 3 class definition.
Interface for the Fermi LAT point spread function.
Definition: GLATPsf.hpp:54
Filename class.
Definition: GFilename.hpp:62
int version(void) const
Returns PSF version.
Definition: GLATPsf.cpp:485
double costhetamin(void) const
Returns minimum cos(theta) angle.
Definition: GLATPsf.cpp:445
#define G_READ
Definition: GLATPsf.cpp:42
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
virtual ~GLATPsf(void)
Destructor.
Definition: GLATPsf.cpp:121
int size(void) const
Returns size of PSF.
Definition: GLATPsf.cpp:394
GChatter
Definition: GTypemaps.hpp:33
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
Fermi LAT point spread function class definition.
void write(GFits &file) const
Write point spread function into FITS file.
Definition: GLATPsf.cpp:367
GLATPsfBase * m_psf
Pointer to versioned point spread function.
Definition: GLATPsf.hpp:96
std::string m_evtype
Event type.
Definition: GLATPsf.hpp:95
GLATPsf(void)
Void constructor.
Definition: GLATPsf.cpp:65
virtual int version(void) const =0
virtual void read(const GFitsTable &table)=0
int nenergies(void) const
Return number of energies in point spread function.
void load(const GFilename &filename, const std::string &evtype)
Load point spread function from FITS file.
Definition: GLATPsf.cpp:231
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
Exception handler interface definition.
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:941
FITS binary table class definition.
const double & scale_par2(void) const
Return second scaling parameter.
int ncostheta(void) const
Returns number of cos(theta) bins in PSF.
Definition: GLATPsf.cpp:427
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void init_members(void)
Initialise class members.
Definition: GLATPsf.cpp:547
const bool & front(void) const
Signal that point spread function is for front section.
Filename class interface definition.
GLATPsf & operator=(const GLATPsf &psf)
Assignment operator.
Definition: GLATPsf.cpp:143
Fermi/LAT point spread function version 3 class.
Definition: GLATPsfV3.hpp:51
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function information.
Definition: GLATPsf.cpp:501
void free_members(void)
Delete class members.
Definition: GLATPsf.cpp:584