GammaLib  1.7.0.dev
 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-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 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 "GLATPsfV1.hpp"
32 #include "GLATException.hpp"
33 #include "GMath.hpp"
34 #include "GFitsBinTable.hpp"
35 #include "GFitsTableFloatCol.hpp"
36 #include "GIntegral.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 GLATException::inconsistent_response
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) {
217  ncore->number(), size);
218  }
219  if (sigma->number() != size) {
221  sigma->number(), size);
222  }
223  if (gcore->number() != size) {
225  gcore->number(), size);
226  }
227  if (gtail->number() != size) {
229  gtail->number(), size);
230  }
231 
232  // Copy data
233  for (int i = 0; i < size; ++i) {
234  m_ncore.push_back(ncore->real(0,i));
235  m_sigma.push_back(sigma->real(0,i));
236  m_gcore.push_back(gcore->real(0,i));
237  m_gtail.push_back(gtail->real(0,i));
238  }
239 
240  } // endif: there were bins
241 
242  // Return
243  return;
244 }
245 
246 
247 /***********************************************************************//**
248  * @brief Write point spread function into FITS file
249  *
250  * @param[in] file FITS file.
251  *
252  * Writes the PSF into the extension "RPSF" of a FITS file. This method
253  * does not check if a "RPSF" extension exists so far, it simply adds one
254  * each time it is called.
255  *
256  * Nothing is done if the PSF size is 0.
257  *
258  * @todo Check if a RPSF extension exists already in FITS file
259  ***************************************************************************/
260 void GLATPsfV1::write(GFits& file) const
261 {
262  // Continue only if there are bins
263  int size = m_rpsf_bins.size();
264  if (size > 0) {
265 
266  // Create new binary table
267  GFitsBinTable* hdu_rpsf = new GFitsBinTable;
268 
269  // Set table attributes
270  hdu_rpsf->extname("RPSF");
271 
272  // Write boundaries into table
273  m_rpsf_bins.write(*hdu_rpsf);
274 
275  // Allocate floating point vector columns
276  GFitsTableFloatCol col_ncore = GFitsTableFloatCol("NCORE", 1, size);
277  GFitsTableFloatCol col_sigma = GFitsTableFloatCol("SIGMA", 1, size);
278  GFitsTableFloatCol col_gcore = GFitsTableFloatCol("GCORE", 1, size);
279  GFitsTableFloatCol col_gtail = GFitsTableFloatCol("GTAIL", 1, size);
280 
281  // Fill columns
282  for (int i = 0; i < size; ++i) {
283  col_ncore(0,i) = m_ncore[i];
284  col_sigma(0,i) = m_sigma[i];
285  col_gcore(0,i) = m_gcore[i];
286  col_gtail(0,i) = m_gtail[i];
287  }
288 
289  // Append columns to table
290  hdu_rpsf->append(col_ncore);
291  hdu_rpsf->append(col_sigma);
292  hdu_rpsf->append(col_gcore);
293  hdu_rpsf->append(col_gtail);
294 
295  // Set detector section
296  std::string detnam = (front()) ? "FRONT" : "BACK";
297 
298  // Set header keywords
299  hdu_rpsf->card("PSFVER", 1, "File format version");
300  hdu_rpsf->card("DETNAM", detnam, "Detector section");
301 
302  // Append HDU to FITS file
303  file.append(*hdu_rpsf);
304 
305  // Free binary table
306  delete hdu_rpsf;
307 
308  } // endif: there were data to write
309 
310  // Return
311  return;
312 }
313 
314 
315 /***********************************************************************//**
316  * @brief Return point spread function value
317  *
318  * @param[in] offset Offset angle (deg).
319  * @param[in] logE Log10 of the true photon energy (MeV).
320  * @param[in] ctheta Cosine of zenith angle.
321  *
322  * @todo Some optimisation could be done as in many cases gcore==gtail,
323  * and for this special case ntail=ncore, hence things become a little
324  * simpler.
325  ***************************************************************************/
326 double GLATPsfV1::psf(const double& offset, const double& logE,
327  const double& ctheta)
328 {
329  // Set constants
330  const double ub = 10.0;
331 
332  // Initialise response
333  double psf = 0.0;
334 
335  // Compute point spread function
336  if (ctheta >= m_min_ctheta) {
337 
338  // Get response parameters
339  double ncore = m_rpsf_bins.interpolate(logE, ctheta, m_ncore);
340  double sigma = m_rpsf_bins.interpolate(logE, ctheta, m_sigma);
341  double gcore = m_rpsf_bins.interpolate(logE, ctheta, m_gcore);
342  double gtail = m_rpsf_bins.interpolate(logE, ctheta, m_gtail);
343 
344  // Compute energy in MeV
345  double energy = pow(10.0, logE);
346 
347  // Rescale the sigma value after interpolation
348  sigma *= scale_factor(energy);
349 
350  // Compute base function argument
351  double r = gammalib::deg2rad * offset / sigma;
352  double u = 0.5 * r * r;
353 
354  // Compute normalization of tail
355  double ntail = ncore * (base_fct(ub, gcore) / base_fct(ub, gtail));
356 
357  // Ensure that PSF integrates to unity. For small energies perform
358  // a numerical integration over the solid angle, while for larger
359  // energies use a small angle approximation.
360  if (energy < 120.0) {
361  GLATPsfV1::base_integrand integrand(ncore, ntail, sigma, gcore, gtail);
362  GIntegral integral(&integrand);
363  ncore /= integral.romberg(0.0, gammalib::pihalf) * gammalib::twopi;
364  }
365  else {
366  double rmax = gammalib::pihalf / sigma;
367  double umax = 0.5 * rmax * rmax;
368  double norm = ncore*base_int(umax, gcore) + ntail*base_int(umax, gtail);
369  ncore /= norm * gammalib::twopi * sigma * sigma;
370  }
371 
372  // Re-compute normalization of tail
373  ntail = ncore * (base_fct(ub, gcore) / base_fct(ub, gtail));
374 
375  // Compute PSF value
376  psf = ncore * base_fct(u, gcore) + ntail * base_fct(u, gtail);
377 
378  // Compile option: check PSF normalization
379  #if G_CHECK_PSF_NORM
380  GLATPsfV1::base_integrand integrand(ncore, ntail, sigma, gcore, gtail);
381  GIntegral integral(&integrand);
382  double sum = integral.romberg(0.0, gammalib::pihalf) * gammalib::twopi;
383  std::cout << "Energy=" << energy;
384  std::cout << " Offset=" << offset;
385  std::cout << " cos(theta)=" << ctheta;
386  std::cout << " error=" << sum-1.0 << std::endl;
387  #endif
388  }
389 
390  // Return point spread function
391  return psf;
392 }
393 
394 
395 /***********************************************************************//**
396  * @brief Print point spread function
397  *
398  * @param[in] chatter Chattiness (defaults to NORMAL).
399  * @return String containing point spread function information.
400  ***************************************************************************/
401 std::string GLATPsfV1::print(const GChatter& chatter) const
402 {
403  // Initialise result string
404  std::string result;
405 
406  // Continue only if chatter is not silent
407  if (chatter != SILENT) {
408 
409  // Append header
410  result.append("=== GLATPsfV1 ===");
411 
412  } // endif: chatter was not silent
413 
414  // Return result
415  return result;
416 }
417 
418 
419 /*==========================================================================
420  = =
421  = Private methods =
422  = =
423  ==========================================================================*/
424 
425 /***********************************************************************//**
426  * @brief Initialise class members
427  ***************************************************************************/
429 {
430  // Initialise members
431  m_ncore.clear();
432  m_sigma.clear();
433  m_gcore.clear();
434  m_gtail.clear();
435 
436  // Return
437  return;
438 }
439 
440 
441 /***********************************************************************//**
442  * @brief Copy class members
443  *
444  * @param[in] psf Point spread function.
445  ***************************************************************************/
447 {
448  // Copy members
449  m_ncore = psf.m_ncore;
450  m_sigma = psf.m_sigma;
451  m_gcore = psf.m_gcore;
452  m_gtail = psf.m_gtail;
453 
454  // Return
455  return;
456 }
457 
458 
459 /***********************************************************************//**
460  * @brief Delete class members
461  ***************************************************************************/
463 {
464  // Return
465  return;
466 }
467 
468 
469 /***********************************************************************//**
470  * @brief Return point spread base function value
471  *
472  * @param[in] u Function argument.
473  * @param[in] gamma Index.
474  *
475  * The version 1 PSF base function is given by
476  * \f[\left(1 - \frac{1}{\Gamma} \right)
477  * \left(1 + \frac{u}{\Gamma} \right)^{-\Gamma}\f]
478  ***************************************************************************/
479 double GLATPsfV1::base_fct(const double& u, const double& gamma)
480 {
481  // Get base function value. The special case of gamma==1 is a ugly
482  // kluge because of sloppy programming in handoff response when
483  // setting boundaries of fit parameters for the PSF.
484  double base = (gamma == 1)
485  ? (1.0 - 1.0/1.001) * std::pow(1.0 + u/1.001, -1.001)
486  : (1.0 - 1.0/gamma) * std::pow(1.0 + u/gamma, -gamma);
487 
488  // Return base function
489  return base;
490 }
491 
492 
493 /***********************************************************************//**
494  * @brief Return approximation of point spread base function integral
495  *
496  * @param[in] u Function argument.
497  * @param[in] gamma Index.
498  *
499  * The version 1 PSF base function integral is approximated by
500  * \f[1 - \left(1 + \frac{u}{\Gamma} \right)^{1-\Gamma}\f]
501  * which is valid for small angles \f$u\f$. For larger angles a numerical
502  * integration of the base function has to be performed.
503  *
504  * @todo Verify that 1+u/gamma is not negative
505  ***************************************************************************/
506 double GLATPsfV1::base_int(const double& u, const double& gamma)
507 {
508  // Compute integral of base function
509  double integral = 1.0 - std::pow(1.0 + u/gamma, 1.0 - gamma);
510 
511  // Return integral
512  return integral;
513 }
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:821
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
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:446
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
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:260
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:45
FITS file class.
Definition: GFits.hpp:63
void init_members(void)
Initialise class members.
Definition: GLATPsfV1.cpp:428
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:401
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:326
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
LAT exception handler interface definition.
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:462
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:479
static double base_int(const double &u, const double &gamma)
Return approximation of point spread base function integral.
Definition: GLATPsfV1.cpp:506
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:1332
int size(void) const
Return number of bins in point spread function.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:665
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.