GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GLATLtCubeMap.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GLATLtCubeMap.cpp - Fermi LAT livetime cube map class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-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 GLATLtCubeMap.cpp
23  * @brief Fermi/LAT livetime cube map 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 "GMath.hpp"
33 #include "GEnergy.hpp"
34 #include "GFitsBinTable.hpp"
35 #include "GLATAeff.hpp"
36 #include "GLATPsf.hpp"
37 #include "GLATLtCubeMap.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_COSTHETA "GLATLtCubeMap::costheta(int&)"
41 #define G_PHI "GLATLtCubeMap::phi(int&)"
42 
43 /* __ Macros _____________________________________________________________ */
44 
45 /* __ Coding definitions _________________________________________________ */
46 
47 /* __ Debug definitions __________________________________________________ */
48 
49 /* __ Constants __________________________________________________________ */
50 
51 
52 /*==========================================================================
53  = =
54  = Constructors/destructors =
55  = =
56  ==========================================================================*/
57 
58 /***********************************************************************//**
59  * @brief Void constructor
60  ***************************************************************************/
62 {
63  // Initialise class members
64  init_members();
65 
66  // Return
67  return;
68 }
69 
70 
71 /***********************************************************************//**
72  * @brief Copy constructor
73  *
74  * @param map Livetime cube map.
75  ***************************************************************************/
77 {
78  // Initialise class members
79  init_members();
80 
81  // Copy members
82  copy_members(map);
83 
84  // Return
85  return;
86 }
87 
88 
89 /***********************************************************************//**
90  * @brief Destructor
91  ***************************************************************************/
93 {
94  // Free members
95  free_members();
96 
97  // Return
98  return;
99 }
100 
101 
102 /*==========================================================================
103  = =
104  = Operators =
105  = =
106  ==========================================================================*/
107 
108 /***********************************************************************//**
109  * @brief Assignment operator
110  *
111  * @param[in] map Livetime cube map.
112  * @return Livetime cube map.
113  ***************************************************************************/
115 {
116  // Execute only if object is not identical
117  if (this != &map) {
118 
119  // Free members
120  free_members();
121 
122  // Initialise private members
123  init_members();
124 
125  // Copy members
126  copy_members(map);
127 
128  } // endif: object was not identical
129 
130  // Return this object
131  return *this;
132 }
133 
134 
135 /***********************************************************************//**
136  * @brief Sum function multiplied by livetime over zenith angle
137  *
138  * @param[in] dir Sky direction.
139  * @param[in] fct Function to evaluate.
140  *
141  * Computes
142  * \f[\sum_{\cos \theta} T_{\rm live}(\cos \theta) f(\cos \theta)\f]
143  * where
144  * \f$T_{\rm live}(\cos \theta)\f$ is the livetime as a function of the
145  * cosine of the zenith angle, and
146  * \f$f(\cos \theta)\f$ is a function that depends on the cosine of the
147  * zenith angle.
148  * This method assumes that \f$T_{\rm live}(\cos \theta)\f$ is stored as a
149  * set of maps.
150  ***************************************************************************/
151 double GLATLtCubeMap::operator()(const GSkyDir& dir, _ltcube_ctheta fct) const
152 {
153  // Get map index
154  int pixel = m_map.dir2pix(dir);
155 
156  // Initialise sum
157  double sum = 0.0;
158 
159  // Loop over zenith angles
160  for (int i = 0; i < m_num_ctheta; ++i) {
161  sum += m_map(pixel, i) * (*fct)(costheta(i));
162  }
163 
164  // Return sum
165  return sum;
166 }
167 
168 
169 /***********************************************************************//**
170  * @brief Sum function multiplied by livetime over zenith and azimuth angles
171  *
172  * @param[in] dir Sky direction.
173  * @param[in] fct Function to evaluate.
174  *
175  * Computes
176  * \f[\sum_{\cos \theta, \phi} T_{\rm live}(\cos \theta, \phi)
177  * f(\cos \theta, \phi)\f]
178  * where
179  * \f$T_{\rm live}(\cos \theta, \phi)\f$ is the livetime as a function of
180  * the cosine of the zenith and of the azimuth angle, and
181  * \f$f(\cos \theta, \phi)\f$ is a function that depends on the cosine of
182  * the zenith angle and of the azimuth angle.
183  * This method assumes that \f$T_{\rm live}(\cos \theta, \phi)\f$ is
184  * stored as a set of maps in a 2D array with \f$\cos \theta\f$ being the
185  * most rapidely varying parameter and with the first map starting at
186  * index m_num_ctheta (the first m_num_ctheta maps are the livetime cube
187  * maps without any \f$\phi\f$ dependence).
188  ***************************************************************************/
190 {
191  // Get map index
192  int pixel = m_map.dir2pix(dir);
193 
194  // Initialise sum
195  double sum = 0.0;
196 
197  // Loop over azimuth and zenith angles. Note that the map index starts
198  // with m_num_ctheta as the first m_num_ctheta maps correspond to an
199  // evaluation without any phi-dependence.
200  for (int iphi = 0, i = m_num_ctheta; iphi < m_num_phi; ++iphi) {
201  double p = phi(iphi);
202  for (int itheta = 0; itheta < m_num_ctheta; ++itheta, ++i) {
203  sum += m_map(pixel, i) * (*fct)(costheta(itheta), p);
204  }
205  }
206 
207  // Return sum
208  return sum;
209 }
210 
211 
212 /***********************************************************************//**
213  * @brief Sum effective area multiplied by livetime over zenith and
214  * (optionally) azimuth angles
215  *
216  * @param[in] dir True sky direction.
217  * @param[in] energy True photon energy.
218  * @param[in] aeff Effective area.
219  *
220  * Computes
221  * \f[\sum_{\cos \theta, \phi} T_{\rm live}(\cos \theta, \phi)
222  * A_{\rm eff}(\log E, \cos \theta, \phi)\f]
223  * where
224  * \f$T_{\rm live}(\cos \theta, \phi)\f$ is the livetime as a function of
225  * the cosine of the zenith and the azimuth angle, and
226  * \f$A_{\rm eff}(\log E, \cos \theta, \phi)\f$ is the effective area that
227  * depends on
228  * the log10 of the energy (in MeV),
229  * the cosine of the zenith angle, and
230  * the azimuth angle.
231  * This method assumes that \f$T_{\rm live}(\cos \theta, \phi)\f$ is
232  * stored as a set of maps in a 2D array with \f$\cos \theta\f$ being the
233  * most rapidely varying parameter and with the first map starting at
234  * index m_num_ctheta (the first m_num_ctheta maps are the livetime cube
235  * maps without any \f$\phi\f$ dependence).
236  ***************************************************************************/
237 double GLATLtCubeMap::operator()(const GSkyDir& dir, const GEnergy& energy,
238  const GLATAeff& aeff) const
239 {
240  // Get map index
241  int pixel = m_map.dir2pix(dir);
242 
243  // Initialise sum
244  double sum = 0.0;
245 
246  // Circumvent const correctness
247  GLATAeff* fct = ((GLATAeff*)&aeff);
248 
249  // If livetime cube and response have phi dependence then sum over
250  // zenith and azimuth. Note that the map index starts with m_num_ctheta
251  // as the first m_num_ctheta maps correspond to an evaluation without
252  // any phi-dependence.
253  if (has_phi() && aeff.has_phi()) {
254  for (int iphi = 0, i = m_num_ctheta; iphi < m_num_phi; ++iphi) {
255  double p = phi(iphi);
256  for (int itheta = 0; itheta < m_num_ctheta; ++itheta, ++i) {
257  sum += m_map(pixel, i) * (*fct)(energy.log10MeV(), costheta(i), p);
258  }
259  }
260  }
261 
262  // ... otherwise sum only over zenith angle
263  else {
264  for (int i = 0; i < m_num_ctheta; ++i) {
265  sum += m_map(pixel, i) * (*fct)(energy.log10MeV(), costheta(i));
266  }
267  }
268 
269  // Return sum
270  return sum;
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Sum effective area multiplied by livetime over zenith angles
276  *
277  * @param[in] dir True sky direction.
278  * @param[in] energy True photon energy.
279  * @param[in] offset Offset from true direction (deg).
280  * @param[in] psf Point spread function.
281  * @param[in] aeff Effective area.
282  *
283  * Computes
284  * \f[\sum_{\cos \theta} T_{\rm live}(\cos \theta)
285  * PSF(\log E, \delta, \cos \theta) A_{\rm eff}(\cos \theta, \phi)\f]
286  * where
287  * \f$T_{\rm live}(\cos \theta)\f$ is the livetime as a function of the
288  * cosine of the zenith angle, and
289  * \f$PSF(\log E, \delta, \cos \theta)\f$ is the point spread function that
290  * depends on
291  * the log10 of the energy (in MeV),
292  * the offset angle from the true direction (in degrees), and
293  * the cosine of the zenith angle, and
294  * \f$A_{\rm eff}(\cos \theta, \phi)\f$ is the effective area that depends
295  * on the cosine of the zenith angle and (optionally) of the azimuth angle.
296  ***************************************************************************/
297 double GLATLtCubeMap::operator()(const GSkyDir& dir, const GEnergy& energy,
298  const double& offset, const GLATPsf& psf,
299  const GLATAeff& aeff) const
300 {
301  // Get map index
302  int pixel = m_map.dir2pix(dir);
303 
304  // Initialise sum
305  double sum = 0.0;
306 
307  // Circumvent const correctness
308  GLATPsf* fpsf = const_cast<GLATPsf*>(&psf);
309  GLATAeff* faeff = const_cast<GLATAeff*>(&aeff);
310 
311  // Get log10 of energy
312  double logE = energy.log10MeV();
313 
314  // If livetime cube and response have phi dependence then sum over
315  // zenith and azimuth. Note that the map index starts with m_num_ctheta
316  // as the first m_num_ctheta maps correspond to an evaluation without
317  // any phi-dependence.
318  if (has_phi() && aeff.has_phi()) {
319  for (int iphi = 0, i = m_num_ctheta; iphi < m_num_phi; ++iphi) {
320  double p = phi(iphi);
321  for (int itheta = 0; itheta < m_num_ctheta; ++itheta, ++i) {
322  sum += m_map(pixel, i) * (*faeff)(logE, costheta(i), p) *
323  (*fpsf)(offset, logE, costheta(i));
324  }
325  }
326  }
327 
328  // ... otherwise sum only over zenith angle
329  else {
330  for (int i = 0; i < m_num_ctheta; ++i) {
331  sum += m_map(pixel, i) * (*faeff)(logE, costheta(i)) *
332  (*fpsf)(offset, logE, costheta(i));
333  }
334  }
335 
336  // Return sum
337  return sum;
338 }
339 
340 
341 /*==========================================================================
342  = =
343  = Public methods =
344  = =
345  ==========================================================================*/
346 
347 /***********************************************************************//**
348  * @brief Clear livetime cube map
349  ***************************************************************************/
351 {
352  // Free class members
353  free_members();
354 
355  // Initialise members
356  init_members();
357 
358  // Return
359  return;
360 }
361 
362 
363 /***********************************************************************//**
364  * @brief Clone livetime cube map
365  *
366  * @return Pointer to deep copy of livetime cube map
367  ***************************************************************************/
369 {
370  return new GLATLtCubeMap(*this);
371 }
372 
373 
374 /***********************************************************************//**
375  * @brief Load livetime cube from FITS file
376  *
377  * @param[in] table FITS table.
378  ***************************************************************************/
379 void GLATLtCubeMap::read(const GFitsTable& table)
380 {
381  // Clear object
382  clear();
383 
384  // Read skymap
385  m_map.read(table);
386 
387  // Set costheta binning scheme
388  std::string scheme =
390  m_sqrt_bin = (scheme == "SQRT(1-COSTHETA)");
391 
392  // Read attributes
393  m_num_ctheta = table.integer("NBRBINS");
394  m_num_phi = table.integer("PHIBINS");
395  m_min_ctheta = table.real("COSMIN");
396 
397  // Return
398  return;
399 }
400 
401 
402 /***********************************************************************//**
403  * @brief Write livetime cube map into FITS file
404  *
405  * @param[in] fits FITS file.
406  * @param[in] extname Livetime cube map extension name.
407  *
408  * Writes livetime cube map into a FITS file.
409  ***************************************************************************/
410 void GLATLtCubeMap::write(GFits& fits, const std::string& extname) const
411 {
412  // If the FITS file contains already an extension with the same name
413  // then remove this extension now
414  if (fits.contains(extname)) {
415  fits.remove(extname);
416  }
417 
418  // Write sky map info FITS file
419  GFitsHDU* hdu = m_map.write(fits, extname);
420 
421  // If HDU is valid then set sky map keywords
422  if (hdu != NULL) {
423  hdu->card("THETABIN", costhetabin(), "Costheta binning scheme");
424  hdu->card("NBRBINS", m_num_ctheta, "Number of costheta bins");
425  hdu->card("PHIBINS", m_num_phi, "Number of phi bins");
426  hdu->card("COSMIN", m_min_ctheta, "Minimum costheta value");
427  }
428 
429  // Return
430  return;
431 }
432 
433 
434 /***********************************************************************//**
435  * @brief Return cos theta value for an index
436  *
437  * @param[in] index Bin index (starting from 0).
438  *
439  * @exception GException::out_of_range
440  * Bin index outside value range.
441  *
442  * The cos theta value is computed using
443  * \f[\cos \theta = 1 - \left( \frac{{\rm index} + 0.5}{\rm ncostheta}
444  * \right)^N (1 - {\rm costhetamin})\f]
445  * where
446  * \f${\rm index}\f$ is the bin index,
447  * \f${\rm ncostheta}\f$ is the number of cos theta bins,
448  * \f$N\f$ is either 1 or 2, and
449  * \f${\rm costhetamin}\f$ is the minimum cos theta value.
450  * Default values for LAT are \f$N=2\f$ and \f${\rm costhetamin}=0\f$.
451  ***************************************************************************/
452 double GLATLtCubeMap::costheta(const int& index) const
453 {
454  // Optionally check if the index is valid
455  #if defined(G_RANGE_CHECK)
456  if (index < 0 || index >= m_num_ctheta) {
458  }
459  #endif
460 
461  // Set cos theta scale
462  double f = (index+0.5)/m_num_ctheta;
463  if (m_sqrt_bin) {
464  f = f*f;
465  }
466 
467  // Set cos theta value
468  double costheta = 1.0 - f * (1.0 - m_min_ctheta);
469 
470  // Return costheta
471  return costheta;
472 }
473 
474 
475 /***********************************************************************//**
476  * @brief Return phi value (in radians) for an index
477  *
478  * @param[in] index Bin index (starting from 0).
479  *
480  * @exception GException::out_of_range
481  * Bin index outside value range.
482  *
483  * The phi value is computed using
484  * \f[\phi = \frac{{\rm index} + 0.5}{\rm nphi} \frac{\pi}{4}\f]
485  * where
486  * \f${\rm index}\f$ is the bin index, and
487  * \f${\rm nphi}\f$ is the number of phi bins.
488  ***************************************************************************/
489 double GLATLtCubeMap::phi(const int& index) const
490 {
491  // Optionally check if the index is valid
492  #if defined(G_RANGE_CHECK)
493  if (index < 0 || index >= m_num_phi) {
494  throw GException::out_of_range(G_PHI, index, 0, m_num_phi-1);
495  }
496  #endif
497 
498  // Set phi value
499  double phi = (index+0.5) / m_num_phi * gammalib::pi / 4.0;
500 
501  // Return phi
502  return phi;
503 }
504 
505 
506 /***********************************************************************//**
507  * @brief Return cos theta binning scheme
508  *
509  * Returns either "SQRT(1-COSTHETA)" or "COSTHETA".
510  ***************************************************************************/
511 std::string GLATLtCubeMap::costhetabin(void) const
512 {
513  // Set string
514  std::string result = (m_sqrt_bin) ? "SQRT(1-COSTHETA)" : "COSTHETA";
515 
516  // Return result
517  return result;
518 }
519 
520 
521 /***********************************************************************//**
522  * @brief Print livetime cube map information
523  *
524  * @param[in] chatter Chattiness (defaults to NORMAL).
525  * @return String containing livetime cube map information.
526  ***************************************************************************/
527 std::string GLATLtCubeMap::print(const GChatter& chatter) const
528 {
529  // Initialise result string
530  std::string result;
531 
532  // Continue only if chatter is not silent
533  if (chatter != SILENT) {
534 
535  // Append header
536  result.append("=== GLATLtCubeMap ===");
537 
538  // Append information
539  result.append("\n"+gammalib::parformat("Number of cos theta bins") +
541  result.append("\n"+gammalib::parformat("Number of phi bins") +
542  gammalib::str(nphi()));
543  result.append("\n"+gammalib::parformat("Cos theta binning"));
544  if (m_sqrt_bin) {
545  result.append("sqrt");
546  }
547  else {
548  result.append("linear");
549  }
550  result.append("\n"+gammalib::parformat("Minimum cos theta") +
552  result.append("\n"+m_map.print(gammalib::reduce(chatter)));
553 
554  } // endif: chatter was not silent
555 
556  // Return result
557  return result;
558 }
559 
560 
561 /*==========================================================================
562  = =
563  = Private methods =
564  = =
565  ==========================================================================*/
566 
567 /***********************************************************************//**
568  * @brief Initialise class members
569  ***************************************************************************/
571 {
572  // Initialise members
573  m_map.clear();
574  m_num_ctheta = 0;
575  m_num_phi = 0;
576  m_min_ctheta = 0.0;
577  m_sqrt_bin = true;
578 
579  // Return
580  return;
581 }
582 
583 
584 /***********************************************************************//**
585  * @brief Copy class members
586  *
587  * @param map Livetime cube map.
588  ***************************************************************************/
590 {
591  // Copy members
592  m_map = map.m_map;
594  m_num_phi = map.m_num_phi;
596  m_sqrt_bin = map.m_sqrt_bin;
597 
598  // Return
599  return;
600 }
601 
602 
603 /***********************************************************************//**
604  * @brief Delete class members
605  ***************************************************************************/
607 {
608  // Return
609  return;
610 }
void free_members(void)
Delete class members.
GLATLtCubeMap * clone(void) const
Clone livetime cube map.
bool m_sqrt_bin
Square root binning?
double(* _ltcube_ctheta_phi)(const double &costheta, const double &phi)
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
Energy value class definition.
GLATLtCubeMap & operator=(const GLATLtCubeMap &cube)
Assignment operator.
const double pi
Definition: GMath.hpp:35
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
double m_min_ctheta
Minimum cos theta value.
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
bool has_phi(void) const
Signal if livetime cube map has phi dependence.
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
int m_num_ctheta
Number of bins in cos theta.
const double & costhetamin(void) const
Return minimum of cosine theta.
#define G_PHI
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2557
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
Interface for the Fermi LAT point spread function.
Definition: GLATPsf.hpp:54
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
virtual ~GLATLtCubeMap(void)
Destructor.
void write(GFits &fits, const std::string &extname=gammalib::extname_lat_ltcubemap) const
Write livetime cube map into FITS file.
std::string costhetabin(void) const
Return cos theta binning scheme.
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2467
int m_num_phi
Number of bins in phi.
const int & nphi(void) const
Return number of phi bins.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GChatter
Definition: GTypemaps.hpp:33
Interface for the Fermi/LAT effective area.
Definition: GLATAeff.hpp:59
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
Fermi LAT point spread function class definition.
bool has_phi(void) const
Signal that effective area has Phi dependence.
Definition: GLATAeff.hpp:207
void copy_members(const GLATLtCubeMap &cube)
Copy class members.
Fermi LAT livetime cube map class definition.
void clear(void)
Clear livetime cube map.
double phi(const int &index) const
Return phi value (in radians) for an index.
GSkyPixel dir2pix(const GSkyDir &dir) const
Returns sky map pixel for a given sky direction.
Definition: GSkyMap.cpp:1462
GSkyMap m_map
Lifetime cube map.
double operator()(const GSkyDir &dir, _ltcube_ctheta fct) const
Sum function multiplied by livetime over zenith angle.
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
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 int & ncostheta(void) const
Return number of cosine theta bins.
GLATLtCubeMap(void)
Void constructor.
double(* _ltcube_ctheta)(const double &costheta)
Fermi LAT effective area class definition.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Sky direction class.
Definition: GSkyDir.hpp:62
void init_members(void)
Initialise class members.
#define G_COSTHETA
void read(const GFitsTable &table)
Load livetime cube from FITS file.
Fermi LAT livetime cube map class.
double costheta(const int &index) const
Return cos theta value for an index.
Mathematical function definitions.
std::string print(const GChatter &chatter=NORMAL) const
Print livetime cube map information.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413