GammaLib  1.7.0.dev
 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-2018 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 "GTools.hpp"
32 #include "GFilename.hpp"
33 #include "GFits.hpp"
34 #include "GFitsBinTable.hpp"
35 #include "GFitsTableFloatCol.hpp"
36 #include "GLATPsf.hpp"
37 #include "GLATPsfV1.hpp"
38 #include "GLATPsfV3.hpp"
39 #include "GLATException.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_response
275  * Invalid response type or 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:
336  "Unsupported response function version "+gammalib::str(version)+".");
337  break;
338  }
339 
340  // Set PSF attributes (needed before reading scaling parameters for
341  // Pass 6 and Pass 7 response functions)
342  m_psf->front(front);
343 
344  // Read PSF scaling parameters
345  m_psf->read_scale(hdu_scale);
346 
347  // Read PSF
348  m_psf->read(hdu_rpsf);
349 
350  // Return
351  return;
352 }
353 
354 
355 /***********************************************************************//**
356  * @brief Write point spread function into FITS file
357  *
358  * @param[in] fits FITS file.
359  *
360  * Writes the version dependent point spread function into the FITS file. The
361  * method does nothing if no PSF has been allocated.
362  *
363  * @todo Implement PSF versions 2 and 3.
364  ***************************************************************************/
365 void GLATPsf::write(GFits& fits) const
366 {
367  // Continue only if PSF is valid
368  if (m_psf != NULL) {
369 
370  // Write point spread function
371  m_psf->write(fits);
372 
373  // Write scaling parameters
374  m_psf->write_scale(fits);
375 
376  }
377 
378  // Return
379  return;
380 }
381 
382 
383 /***********************************************************************//**
384  * @brief Returns size of PSF
385  *
386  * @return Size of point spread function.
387  *
388  * Returns the size of the PSF which is defined as the number of energy bins
389  * times the number of cos(theta) bins. If no PSF has been allocated, 0 is
390  * returned.
391  ***************************************************************************/
392 int GLATPsf::size(void) const
393 {
394  // Return size of PSF
395  return (nenergies()*ncostheta());
396 }
397 
398 
399 /***********************************************************************//**
400  * @brief Returns number of energy bins in PSF
401  *
402  * @return Number of energy bins.
403  *
404  * Returns the number of energy bins in PSF. If no PSF has been allocated,
405  * 0 is returned.
406  ***************************************************************************/
407 int GLATPsf::nenergies(void) const
408 {
409  // Retrieve number of energy bins
410  int nenergies = (m_psf != NULL) ? m_psf->nenergies() : 0;
411 
412  // Return number of energy bins
413  return nenergies;
414 }
415 
416 
417 /***********************************************************************//**
418  * @brief Returns number of cos(theta) bins in PSF
419  *
420  * @return Number of cos(theta) bins.
421  *
422  * Returns the number of cos(theta) bins in PSF. If no PSF has been
423  * allocated, 0 is returned.
424  ***************************************************************************/
425 int GLATPsf::ncostheta(void) const
426 {
427  // Retrieve number of cos(theta) bins
428  int ncostheta = (m_psf != NULL) ? m_psf->ncostheta() : 0;
429 
430  // Return number of cos(theta) bins
431  return ncostheta;
432 }
433 
434 
435 /***********************************************************************//**
436  * @brief Returns minimum cos(theta) angle
437  *
438  * @return Minimum cos(theta) angle.
439  *
440  * Returns the minimum cos(theta) angle for point spread function access. If
441  * no PSF has been allocated, 0 is returned.
442  ***************************************************************************/
443 double GLATPsf::costhetamin(void) const
444 {
445  // Retrieve number of energy bins
446  double costhetamin = (m_psf != NULL) ? m_psf->costhetamin() : 0.0;
447 
448  // Return minimum cos(theta)
449  return costhetamin;
450 }
451 
452 
453 /***********************************************************************//**
454  * @brief Set minimum cos(theta) angle for point spread function access
455  *
456  * @param[in] ctheta Cosine of maximum zenith angle.
457  *
458  * Sets the minimum cos(theta) angle. No verification of the specified
459  * maximum cos(theta) value is done.
460  ***************************************************************************/
461 void GLATPsf::costhetamin(const double& ctheta)
462 {
463  // Continue only if PSF is valid
464  if (m_psf != NULL) {
465 
466  // Set minimum cos(theta) value
467  m_psf->costhetamin(ctheta);
468 
469  }
470 
471  // Return
472  return;
473 }
474 
475 
476 /***********************************************************************//**
477  * @brief Returns PSF version
478  *
479  * Returns the PSF version number.
480  *
481  * Returns 0 if no PSF has been allocated.
482  ***************************************************************************/
483 int GLATPsf::version(void) const
484 {
485  // Retrieve version number
486  int version = (m_psf != NULL) ? m_psf->version() : 0;
487 
488  // Return version number
489  return version;
490 }
491 
492 
493 /***********************************************************************//**
494  * @brief Print point spread function information
495  *
496  * @param[in] chatter Chattiness.
497  * @return String containing point spread function information.
498  ***************************************************************************/
499 std::string GLATPsf::print(const GChatter& chatter) const
500 {
501  // Initialise result string
502  std::string result;
503 
504  // Continue only if chatter is not silent
505  if (chatter != SILENT) {
506 
507  // Append header
508  result.append("=== GLATPsf ===");
509 
510  // No PSF has been loaded ...
511  if (m_psf == NULL) {
512  result.append("\n"+gammalib::parformat("Version"));
513  result.append("No PSF loaded");
514  }
515 
516  // ... PSF has been loaded
517  else {
518  result.append("\n"+gammalib::parformat("Version"));
519  result.append(gammalib::str(version()));
520  result.append("\n"+gammalib::parformat("Event type"));
521  result.append(evtype());
522  result.append("\n"+gammalib::parformat("Energy scaling"));
523  result.append("sqrt(");
524  result.append("("+gammalib::str(m_psf->scale_par1())+"*(E/100)^");
525  result.append(gammalib::str(m_psf->scale_index())+")^2");
526  result.append(" + ("+gammalib::str(m_psf->scale_par2())+")^2)");
527  }
528 
529  } // endif: chatter was not silent
530 
531  // Return result
532  return result;
533 }
534 
535 
536 /*==========================================================================
537  = =
538  = Private methods =
539  = =
540  ==========================================================================*/
541 
542 /***********************************************************************//**
543  * @brief Initialise class members
544  ***************************************************************************/
546 {
547  // Initialise members
548  m_evtype.clear();
549  m_psf = NULL;
550 
551  // Return
552  return;
553 }
554 
555 
556 /***********************************************************************//**
557  * @brief Copy class members
558  *
559  * @param[in] psf Point spread function.
560  ***************************************************************************/
562 {
563  // Copy members
564  m_evtype = psf.m_evtype;
565 
566  // Clone members
567  if (psf.m_psf != NULL) {
568  m_psf = psf.m_psf->clone();
569  }
570  else {
571  m_psf = NULL;
572  }
573 
574  // Return
575  return;
576 }
577 
578 
579 /***********************************************************************//**
580  * @brief Delete class members
581  ***************************************************************************/
583 {
584  // Free PSF
585  if (m_psf != NULL) delete m_psf;
586 
587  // Signal that PSF is free
588  m_psf = NULL;
589 
590  // Return
591  return;
592 }
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:472
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:561
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:73
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:407
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1267
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
LAT exception handler interface definition.
Filename class.
Definition: GFilename.hpp:62
int version(void) const
Returns PSF version.
Definition: GLATPsf.cpp:483
double costhetamin(void) const
Returns minimum cos(theta) angle.
Definition: GLATPsf.cpp:443
#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:392
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:365
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
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:820
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:425
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void init_members(void)
Initialise class members.
Definition: GLATPsf.cpp:545
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:413
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function information.
Definition: GLATPsf.cpp:499
void free_members(void)
Delete class members.
Definition: GLATPsf.cpp:582