GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GLATPsfV1.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GLATPsfV1.cpp - Fermi/LAT point spread function version 1 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 GLATPsfV1.cpp
23  * @brief Fermi/LAT point spread function version 1 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 "GMath.hpp"
33 #include "GFitsBinTable.hpp"
34 #include "GFitsTableFloatCol.hpp"
35 #include "GIntegral.hpp"
36 #include "GLATPsfV1.hpp"
37 
38 /* __ Method name definitions ____________________________________________ */
39 #define G_READ "GLATPsfV1::read(GFitsTable&)"
40 
41 /* __ Macros _____________________________________________________________ */
42 
43 /* __ Coding definitions _________________________________________________ */
44 
45 /* __ Debug definitions __________________________________________________ */
46 #define G_CHECK_PSF_NORM 0 //!< Check PSF normalization
47 
48 /* __ Constants __________________________________________________________ */
49 
50 
51 /*==========================================================================
52  = =
53  = Constructors/destructors =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Void constructor
59  ***************************************************************************/
61 {
62  // Initialise class members
63  init_members();
64 
65  // Return
66  return;
67 }
68 
69 
70 /***********************************************************************//**
71  * @brief Copy constructor
72  *
73  * @param[in] psf Point spread function.
74  ***************************************************************************/
76 {
77  // Initialise class members
78  init_members();
79 
80  // Copy members
81  copy_members(psf);
82 
83  // Return
84  return;
85 }
86 
87 
88 /***********************************************************************//**
89  * @brief Destructor
90  ***************************************************************************/
92 {
93  // Free members
94  free_members();
95 
96  // Return
97  return;
98 }
99 
100 
101 /*==========================================================================
102  = =
103  = Operators =
104  = =
105  ==========================================================================*/
106 
107 /***********************************************************************//**
108  * @brief Assignment operator
109  *
110  * @param[in] psf Point spread function.
111  * @return Point spread function.
112  ***************************************************************************/
114 {
115  // Execute only if object is not identical
116  if (this != &psf) {
117 
118  // Copy base class members
119  this->GLATPsfBase::operator=(psf);
120 
121  // Free members
122  free_members();
123 
124  // Initialise private members
125  init_members();
126 
127  // Copy members
128  copy_members(psf);
129 
130  } // endif: object was not identical
131 
132  // Return this object
133  return *this;
134 }
135 
136 
137 /*==========================================================================
138  = =
139  = Public methods =
140  = =
141  ==========================================================================*/
142 
143 /***********************************************************************//**
144  * @brief Clear point spread function
145  ***************************************************************************/
147 {
148  // Free class members (base and derived classes, derived class first)
149  free_members();
151 
152  // Initialise members
154  init_members();
155 
156  // Return
157  return;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief Clone point spread function
163  *
164  * @return Pointer to deep copy of point spread function.
165  ***************************************************************************/
167 {
168  return new GLATPsfV1(*this);
169 }
170 
171 
172 /***********************************************************************//**
173  * @brief Read point spread function from FITS table
174  *
175  * @param[in] table FITS table.
176  *
177  * @exception GException::invalid_argument
178  * Inconsistent response table encountered
179  *
180  * Reads point spread function information from FITS HDU. In addition to the
181  * energy and costheta binning information, 4 columns are expected:
182  * NCORE, SIGMA, GCORE, and GTAIL.
183  ***************************************************************************/
184 void GLATPsfV1::read(const GFitsTable& table)
185 {
186  // Clear arrays
187  m_ncore.clear();
188  m_sigma.clear();
189  m_gcore.clear();
190  m_gtail.clear();
191 
192  // Get energy and cos theta binning
193  m_rpsf_bins.read(table);
194 
195  // Set minimum cos(theta)
197 
198  // Continue only if there are bins
199  int size = m_rpsf_bins.size();
200  if (size > 0) {
201 
202  // Allocate arrays
203  m_ncore.reserve(size);
204  m_sigma.reserve(size);
205  m_gcore.reserve(size);
206  m_gtail.reserve(size);
207 
208  // Get pointer to columns
209  const GFitsTableCol* ncore = table["NCORE"];
210  const GFitsTableCol* sigma = table["SIGMA"];
211  const GFitsTableCol* gcore = table["GCORE"];
212  const GFitsTableCol* gtail = table["GTAIL"];
213 
214  // Check consistency of columns
215  if (ncore->number() != size) {
216  std::string msg = "Number of elements in \"NCORE\" column ("+
217  gammalib::str(ncore->number())+") is incompatible "
218  "with the expected size ("+
219  gammalib::str(size)+"). Please specify a valid "
220  "point spread function table.";
222  }
223  if (sigma->number() != size) {
224  std::string msg = "Number of elements in \"SIGMA\" column ("+
225  gammalib::str(sigma->number())+") is incompatible "
226  "with the expected size ("+
227  gammalib::str(size)+"). Please specify a valid "
228  "point spread function table.";
230  }
231  if (gcore->number() != size) {
232  std::string msg = "Number of elements in \"GCORE\" column ("+
233  gammalib::str(gcore->number())+") is incompatible "
234  "with the expected size ("+
235  gammalib::str(size)+"). Please specify a valid "
236  "point spread function table.";
238  }
239  if (gtail->number() != size) {
240  std::string msg = "Number of elements in \"GTAIL\" column ("+
241  gammalib::str(gtail->number())+") is incompatible "
242  "with the expected size ("+
243  gammalib::str(size)+"). Please specify a valid "
244  "point spread function table.";
246  }
247 
248  // Copy data
249  for (int i = 0; i < size; ++i) {
250  m_ncore.push_back(ncore->real(0,i));
251  m_sigma.push_back(sigma->real(0,i));
252  m_gcore.push_back(gcore->real(0,i));
253  m_gtail.push_back(gtail->real(0,i));
254  }
255 
256  } // endif: there were bins
257 
258  // Return
259  return;
260 }
261 
262 
263 /***********************************************************************//**
264  * @brief Write point spread function into FITS file
265  *
266  * @param[in] file FITS file.
267  *
268  * Writes the PSF into the extension "RPSF" of a FITS file. This method
269  * does not check if a "RPSF" extension exists so far, it simply adds one
270  * each time it is called.
271  *
272  * Nothing is done if the PSF size is 0.
273  *
274  * @todo Check if a RPSF extension exists already in FITS file
275  ***************************************************************************/
276 void GLATPsfV1::write(GFits& file) const
277 {
278  // Continue only if there are bins
279  int size = m_rpsf_bins.size();
280  if (size > 0) {
281 
282  // Create new binary table
283  GFitsBinTable* hdu_rpsf = new GFitsBinTable;
284 
285  // Set table attributes
286  hdu_rpsf->extname("RPSF");
287 
288  // Write boundaries into table
289  m_rpsf_bins.write(*hdu_rpsf);
290 
291  // Allocate floating point vector columns
292  GFitsTableFloatCol col_ncore = GFitsTableFloatCol("NCORE", 1, size);
293  GFitsTableFloatCol col_sigma = GFitsTableFloatCol("SIGMA", 1, size);
294  GFitsTableFloatCol col_gcore = GFitsTableFloatCol("GCORE", 1, size);
295  GFitsTableFloatCol col_gtail = GFitsTableFloatCol("GTAIL", 1, size);
296 
297  // Fill columns
298  for (int i = 0; i < size; ++i) {
299  col_ncore(0,i) = m_ncore[i];
300  col_sigma(0,i) = m_sigma[i];
301  col_gcore(0,i) = m_gcore[i];
302  col_gtail(0,i) = m_gtail[i];
303  }
304 
305  // Append columns to table
306  hdu_rpsf->append(col_ncore);
307  hdu_rpsf->append(col_sigma);
308  hdu_rpsf->append(col_gcore);
309  hdu_rpsf->append(col_gtail);
310 
311  // Set detector section
312  std::string detnam = (front()) ? "FRONT" : "BACK";
313 
314  // Set header keywords
315  hdu_rpsf->card("PSFVER", 1, "File format version");
316  hdu_rpsf->card("DETNAM", detnam, "Detector section");
317 
318  // Append HDU to FITS file
319  file.append(*hdu_rpsf);
320 
321  // Free binary table
322  delete hdu_rpsf;
323 
324  } // endif: there were data to write
325 
326  // Return
327  return;
328 }
329 
330 
331 /***********************************************************************//**
332  * @brief Return point spread function value
333  *
334  * @param[in] offset Offset angle (deg).
335  * @param[in] logE Log10 of the true photon energy (MeV).
336  * @param[in] ctheta Cosine of zenith angle.
337  *
338  * @todo Some optimisation could be done as in many cases gcore==gtail,
339  * and for this special case ntail=ncore, hence things become a little
340  * simpler.
341  ***************************************************************************/
342 double GLATPsfV1::psf(const double& offset, const double& logE,
343  const double& ctheta)
344 {
345  // Set constants
346  const double ub = 10.0;
347 
348  // Initialise response
349  double psf = 0.0;
350 
351  // Compute point spread function
352  if (ctheta >= m_min_ctheta) {
353 
354  // Get response parameters
355  double ncore = m_rpsf_bins.interpolate(logE, ctheta, m_ncore);
356  double sigma = m_rpsf_bins.interpolate(logE, ctheta, m_sigma);
357  double gcore = m_rpsf_bins.interpolate(logE, ctheta, m_gcore);
358  double gtail = m_rpsf_bins.interpolate(logE, ctheta, m_gtail);
359 
360  // Compute energy in MeV
361  double energy = pow(10.0, logE);
362 
363  // Rescale the sigma value after interpolation
364  sigma *= scale_factor(energy);
365 
366  // Compute base function argument
367  double r = gammalib::deg2rad * offset / sigma;
368  double u = 0.5 * r * r;
369 
370  // Compute normalization of tail
371  double ntail = ncore * (base_fct(ub, gcore) / base_fct(ub, gtail));
372 
373  // Ensure that PSF integrates to unity. For small energies perform
374  // a numerical integration over the solid angle, while for larger
375  // energies use a small angle approximation.
376  if (energy < 120.0) {
377  GLATPsfV1::base_integrand integrand(ncore, ntail, sigma, gcore, gtail);
378  GIntegral integral(&integrand);
379  ncore /= integral.romberg(0.0, gammalib::pihalf) * gammalib::twopi;
380  }
381  else {
382  double rmax = gammalib::pihalf / sigma;
383  double umax = 0.5 * rmax * rmax;
384  double norm = ncore*base_int(umax, gcore) + ntail*base_int(umax, gtail);
385  ncore /= norm * gammalib::twopi * sigma * sigma;
386  }
387 
388  // Re-compute normalization of tail
389  ntail = ncore * (base_fct(ub, gcore) / base_fct(ub, gtail));
390 
391  // Compute PSF value
392  psf = ncore * base_fct(u, gcore) + ntail * base_fct(u, gtail);
393 
394  // Compile option: check PSF normalization
395  #if G_CHECK_PSF_NORM
396  GLATPsfV1::base_integrand integrand(ncore, ntail, sigma, gcore, gtail);
397  GIntegral integral(&integrand);
398  double sum = integral.romberg(0.0, gammalib::pihalf) * gammalib::twopi;
399  std::cout << "Energy=" << energy;
400  std::cout << " Offset=" << offset;
401  std::cout << " cos(theta)=" << ctheta;
402  std::cout << " error=" << sum-1.0 << std::endl;
403  #endif
404  }
405 
406  // Return point spread function
407  return psf;
408 }
409 
410 
411 /***********************************************************************//**
412  * @brief Print point spread function
413  *
414  * @param[in] chatter Chattiness (defaults to NORMAL).
415  * @return String containing point spread function information.
416  ***************************************************************************/
417 std::string GLATPsfV1::print(const GChatter& chatter) const
418 {
419  // Initialise result string
420  std::string result;
421 
422  // Continue only if chatter is not silent
423  if (chatter != SILENT) {
424 
425  // Append header
426  result.append("=== GLATPsfV1 ===");
427 
428  } // endif: chatter was not silent
429 
430  // Return result
431  return result;
432 }
433 
434 
435 /*==========================================================================
436  = =
437  = Private methods =
438  = =
439  ==========================================================================*/
440 
441 /***********************************************************************//**
442  * @brief Initialise class members
443  ***************************************************************************/
445 {
446  // Initialise members
447  m_ncore.clear();
448  m_sigma.clear();
449  m_gcore.clear();
450  m_gtail.clear();
451 
452  // Return
453  return;
454 }
455 
456 
457 /***********************************************************************//**
458  * @brief Copy class members
459  *
460  * @param[in] psf Point spread function.
461  ***************************************************************************/
463 {
464  // Copy members
465  m_ncore = psf.m_ncore;
466  m_sigma = psf.m_sigma;
467  m_gcore = psf.m_gcore;
468  m_gtail = psf.m_gtail;
469 
470  // Return
471  return;
472 }
473 
474 
475 /***********************************************************************//**
476  * @brief Delete class members
477  ***************************************************************************/
479 {
480  // Return
481  return;
482 }
483 
484 
485 /***********************************************************************//**
486  * @brief Return point spread base function value
487  *
488  * @param[in] u Function argument.
489  * @param[in] gamma Index.
490  *
491  * The version 1 PSF base function is given by
492  * \f[\left(1 - \frac{1}{\Gamma} \right)
493  * \left(1 + \frac{u}{\Gamma} \right)^{-\Gamma}\f]
494  ***************************************************************************/
495 double GLATPsfV1::base_fct(const double& u, const double& gamma)
496 {
497  // Get base function value. The special case of gamma==1 is a ugly
498  // kluge because of sloppy programming in handoff response when
499  // setting boundaries of fit parameters for the PSF.
500  double base = (gamma == 1)
501  ? (1.0 - 1.0/1.001) * std::pow(1.0 + u/1.001, -1.001)
502  : (1.0 - 1.0/gamma) * std::pow(1.0 + u/gamma, -gamma);
503 
504  // Return base function
505  return base;
506 }
507 
508 
509 /***********************************************************************//**
510  * @brief Return approximation of point spread base function integral
511  *
512  * @param[in] u Function argument.
513  * @param[in] gamma Index.
514  *
515  * The version 1 PSF base function integral is approximated by
516  * \f[1 - \left(1 + \frac{u}{\Gamma} \right)^{1-\Gamma}\f]
517  * which is valid for small angles \f$u\f$. For larger angles a numerical
518  * integration of the base function has to be performed.
519  *
520  * @todo Verify that 1+u/gamma is not negative
521  ***************************************************************************/
522 double GLATPsfV1::base_int(const double& u, const double& gamma)
523 {
524  // Compute integral of base function
525  double integral = 1.0 - std::pow(1.0 + u/gamma, 1.0 - gamma);
526 
527  // Return integral
528  return integral;
529 }
void read(const GFitsTable &table)
Read point spread function from FITS table.
Definition: GLATPsfV1.cpp:184
void number(const int &number)
Set number of elements in column.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
std::vector< double > m_sigma
PSF sigma parameter.
Definition: GLATPsfV1.hpp:103
Fermi/LAT point spread function version 1 class definition.
void read(const GFitsTable &hdu)
Read response table from FITS table HDU.
Fermi/LAT point spread function version 1 class.
Definition: GLATPsfV1.hpp:47
void copy_members(const GLATPsfV1 &psf)
Copy class members.
Definition: GLATPsfV1.cpp:462
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:944
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
void write(GFitsTable &hdu) const
Write response table into FITS table.
void write(GFits &file) const
Write point spread function into FITS file.
Definition: GLATPsfV1.cpp:276
std::vector< double > m_ncore
PSF ncore parameter.
Definition: GLATPsfV1.hpp:102
FITS table float column class interface definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
FITS file class.
Definition: GFits.hpp:63
void init_members(void)
Initialise class members.
Definition: GLATPsfV1.cpp:444
GLATPsfV1 * clone(void) const
Clone point spread function.
Definition: GLATPsfV1.cpp:166
std::vector< double > m_gtail
PSF gtail parameter.
Definition: GLATPsfV1.hpp:105
void clear(void)
Clear point spread function.
Definition: GLATPsfV1.cpp:146
std::string print(const GChatter &chatter=NORMAL) const
Print point spread function.
Definition: GLATPsfV1.cpp:417
double scale_factor(const double &energy) const
Return scale factor for energy (in MeV)
double interpolate(const double &logE, const double &ctheta, const std::vector< double > &array)
Perform bi-linear interpolation of 2D array.
double psf(const double &offset, const double &logE, const double &ctheta)
Return point spread function value.
Definition: GLATPsfV1.cpp:342
const double pihalf
Definition: GMath.hpp:38
const double deg2rad
Definition: GMath.hpp:43
void free_members(void)
Delete class members.
Abstract Fermi/LAT point spread function base class.
Definition: GLATPsfBase.hpp:47
virtual ~GLATPsfV1(void)
Destructor.
Definition: GLATPsfV1.cpp:91
GLATPsfV1(void)
Void constructor.
Definition: GLATPsfV1.cpp:60
Abstract interface for FITS table column.
GLATPsfV1 & operator=(const GLATPsfV1 &psf)
Assignment operator.
Definition: GLATPsfV1.cpp:113
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GChatter
Definition: GTypemaps.hpp:33
void free_members(void)
Delete class members.
Definition: GLATPsfV1.cpp:478
void init_members(void)
Initialise class members.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
std::vector< double > m_gcore
PSF gcore parameter.
Definition: GLATPsfV1.hpp:104
double m_min_ctheta
Minimum valid cos(theta)
Definition: GLATPsfBase.hpp:97
GLATPsfBase & operator=(const GLATPsfBase &psf)
Assignment operator.
static double base_fct(const double &u, const double &gamma)
Return point spread base function value.
Definition: GLATPsfV1.cpp:495
static double base_int(const double &u, const double &gamma)
Return approximation of point spread base function integral.
Definition: GLATPsfV1.cpp:522
GLATResponseTable m_rpsf_bins
PSF energy and cos theta binning.
Definition: GLATPsfBase.hpp:93
virtual double real(const int &row, const int &inx=0) const =0
FITS binary table class.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
int size(void) const
Return number of bins in point spread function.
Exception handler interface definition.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:678
FITS binary table class definition.
double costheta_lo(const int &inx) const
Return lower bin cos theta [.
FITS table float column.
#define G_READ
Definition: GLATPsfV1.cpp:39
const double twopi
Definition: GMath.hpp:36
Integration class interface definition.
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
const bool & front(void) const
Signal that point spread function is for front section.
int size(void) const
Return number of bins in response table.
Mathematical function definitions.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489