GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GLATLtCube.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GLATLtCube.cpp - Fermi/LAT livetime cube 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 GLATLtCube.cpp
23  * @brief Fermi/LAT livetime cube 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 "GSkyDir.hpp"
34 #include "GEnergy.hpp"
35 #include "GLATAeff.hpp"
36 #include "GLATPsf.hpp"
37 #include "GLATLtCube.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 
41 /* __ Macros _____________________________________________________________ */
42 
43 /* __ Coding definitions _________________________________________________ */
44 
45 /* __ Debug definitions __________________________________________________ */
46 
47 /* __ Constants __________________________________________________________ */
48 
49 
50 /*==========================================================================
51  = =
52  = Constructors/destructors =
53  = =
54  ==========================================================================*/
55 
56 /***********************************************************************//**
57  * @brief Void constructor
58  ***************************************************************************/
60 {
61  // Initialise class members
62  init_members();
63 
64  // Return
65  return;
66 }
67 
68 
69 /***********************************************************************//**
70  * @brief File constructor
71  *
72  * @param[in] filename Livetime cube filename.
73  ***************************************************************************/
75 {
76  // Initialise class members
77  init_members();
78 
79  // Load livetime cube from file
80  load(filename);
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief Copy constructor
89  *
90  * @param[in] cube Livetime cube.
91  ***************************************************************************/
93 {
94  // Initialise class members
95  init_members();
96 
97  // Copy members
98  copy_members(cube);
99 
100  // Return
101  return;
102 }
103 
104 
105 /***********************************************************************//**
106  * @brief Destructor
107  ***************************************************************************/
109 {
110  // Free members
111  free_members();
112 
113  // Return
114  return;
115 }
116 
117 
118 /*==========================================================================
119  = =
120  = Operators =
121  = =
122  ==========================================================================*/
123 
124 /***********************************************************************//**
125  * @brief Assignment operator
126  *
127  * @param[in] cube Livetime cube.
128  * @return Livetime cube.
129  ***************************************************************************/
131 {
132  // Execute only if object is not identical
133  if (this != &cube) {
134 
135  // Free members
136  free_members();
137 
138  // Initialise private members
139  init_members();
140 
141  // Copy members
142  copy_members(cube);
143 
144  } // endif: object was not identical
145 
146  // Return this object
147  return *this;
148 }
149 
150 
151 /***********************************************************************//**
152  * @brief Sum function multiplied by efficiency corrected livetime over
153  * zenith angle
154  *
155  * @param[in] dir Sky direction.
156  * @param[in] energy Energy.
157  * @param[in] fct Function to evaluate.
158  *
159  * Computes
160  * \f[\sum_{\cos \theta} T_{\rm corr.~live}(\cos \theta) f(\cos \theta)\f]
161  * where
162  * \f$T_{\rm corr.~live}(\cos \theta)\f$ is the efficiency corrected
163  * livetime as a function of the cosine of the zenith angle, and
164  * \f$f(\cos \theta)\f$ is a function that depends on the cosine of the
165  * zenith angle.
166  *
167  * Note that no efficieny correction is implemented for this method as the
168  * efficiency factors are stored together with the effective area. If we
169  * want efficiency correction here, we should think about passing this
170  * information to the method.
171  ***************************************************************************/
172 double GLATLtCube::operator()(const GSkyDir& dir, const GEnergy& energy,
173  _ltcube_ctheta fct) const
174 {
175  // Compute exposure
176  double exposure = m_exposure(dir, fct);
177 
178  // Optionally compute livetime factors for trigger rate- and
179  // energy-dependent efficiency corrections
180  /*
181  if (aeff.has_efficiency()) {
182 
183  // Compute correction factors
184  double f1 = aeff.efficiency_factor1(energy);
185  double f2 = aeff.efficiency_factor2(energy);
186 
187  // Compute correction
188  double correction = m_weighted_exposure(dir, fct);
189 
190  // Set exposure
191  exposure = f1 * exposure + f2 * correction;
192 
193  } // endif: corrections requested
194  */
195 
196  // Return exposure
197  return exposure;
198 }
199 
200 
201 /***********************************************************************//**
202  * @brief Sum function multiplied by efficiency corrected livetime over
203  * zenith and azimuth angles
204  *
205  * @param[in] dir Sky direction.
206  * @param[in] energy Energy.
207  * @param[in] fct Function to evaluate.
208  *
209  * Computes
210  * \f[\sum_{\cos \theta, \phi} T_{\rm corr.~live}(\cos \theta, \phi)
211  * f(\cos \theta, \phi)\f]
212  * where
213  * \f$T_{\rm corr.~live}(\cos \theta, \phi)\f$ is the efficiency corrected
214  * livetime as a function of the cosine of the zenith and of the azimuth
215  * angle, and
216  * \f$f(\cos \theta, \phi)\f$ is a function that depends on the cosine of
217  * the zenith angle and of the azimuth angle.
218  *
219  * Note that no efficieny correction is implemented for this method as the
220  * efficiency factors are stored together with the effective area. If we
221  * want efficiency correction here, we should think about passing this
222  * information to the method.
223  ***************************************************************************/
224 double GLATLtCube::operator()(const GSkyDir& dir, const GEnergy& energy,
225  _ltcube_ctheta_phi fct) const
226 {
227  // Compute exposure
228  double exposure = m_exposure(dir, fct);
229 
230  // Optionally compute livetime factors for trigger rate- and
231  // energy-dependent efficiency corrections
232  /*
233  if (aeff.has_efficiency()) {
234 
235  // Compute correction factors
236  double f1 = aeff.efficiency_factor1(energy);
237  double f2 = aeff.efficiency_factor2(energy);
238 
239  // Compute correction
240  double correction = m_weighted_exposure(dir, fct);
241 
242  // Set exposure
243  exposure = f1 * exposure + f2 * correction;
244 
245  } // endif: corrections requested
246  */
247 
248  // Return exposure
249  return exposure;
250 }
251 
252 
253 /***********************************************************************//**
254  * @brief Sum effective area multiplied by efficiency corrected livetime
255  * over zenith and (optionally) azimuth angles
256  *
257  * @param[in] dir Sky direction.
258  * @param[in] energy Energy.
259  * @param[in] aeff Effective area.
260  *
261  * Computes
262  * \f[\sum_{\cos \theta, \phi} T_{\rm corr.~live}(\cos \theta, \phi)
263  * A_{\rm eff}(\cos \theta, \phi)\f]
264  * where
265  * \f$T_{\rm corr.~live}(\cos \theta, \phi)\f$ is the efficiency corrected
266  * livetime as a function of the cosine of the zenith and of the azimuth
267  * angle, and
268  * \f$A_{\rm eff}(\cos \theta, \phi)\f$ is the effective area that depends
269  * on the cosine of the zenith angle and (optionally) of the azimuth angle.
270  ***************************************************************************/
271 double GLATLtCube::operator()(const GSkyDir& dir, const GEnergy& energy,
272  const GLATAeff& aeff) const
273 {
274  // Compute exposure
275  double exposure = m_exposure(dir, energy, aeff);
276 
277  // Optionally compute livetime factors for trigger rate- and
278  // energy-dependent efficiency corrections
279  if (aeff.has_efficiency()) {
280 
281  // Compute correction factors
282  double f1 = aeff.efficiency_factor1(energy);
283  double f2 = aeff.efficiency_factor2(energy);
284 
285  // Compute correction
286  double correction = m_weighted_exposure(dir, energy, aeff);
287 
288  // Set exposure
289  exposure = f1 * exposure + f2 * correction;
290 
291  } // endif: corrections requested
292 
293  // Return exposure
294  return exposure;
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Sum point spread function multiplied by efficiency corrected
300  * livetime over zenith angles
301  *
302  * @param[in] dir Sky direction.
303  * @param[in] energy Energy.
304  * @param[in] offset Offset from true direction (deg).
305  * @param[in] psf Point spread function.
306  * @param[in] aeff Effective area.
307  *
308  * Computes
309  * \f[\sum_{\cos \theta, \phi} T_{\rm corr.~live}(\cos \theta, \phi)
310  * PSF(\log E, \delta, \cos \theta) A_{\rm eff}(\cos \theta, \phi)\f]
311  * where
312  * \f$T_{\rm corr.~live}(\cos \theta, \phi)\f$ is the efficiency corrected
313  * livetime as a function of the cosine of the zenith and of the azimuth
314  * angle,
315  * \f$PSF(\log E, \delta, \cos \theta)\f$ is the point spread function that
316  * depends on
317  * the log10 of the energy (in MeV),
318  * the offset angle from the true direction (in degrees), and
319  * the cosine of the zenith angle, and
320  * \f$A_{\rm eff}(\cos \theta, \phi)\f$ is the effective area that depends
321  * on the cosine of the zenith angle and (optionally) of the azimuth angle.
322  ***************************************************************************/
323 double GLATLtCube::operator()(const GSkyDir& dir, const GEnergy& energy,
324  const double& offset, const GLATPsf& psf,
325  const GLATAeff& aeff) const
326 {
327  // Compute exposure
328  double exposure = m_exposure(dir, energy, offset, psf, aeff);
329 
330  // Optionally compute livetime factors for trigger rate- and
331  // energy-dependent efficiency corrections
332  if (aeff.has_efficiency()) {
333 
334  // Compute correction factors
335  double f1 = aeff.efficiency_factor1(energy);
336  double f2 = aeff.efficiency_factor2(energy);
337 
338  // Compute correction
339  double correction = m_weighted_exposure(dir, energy, offset, psf, aeff);
340 
341  // Set exposure
342  exposure = f1 * exposure + f2 * correction;
343 
344  } // endif: corrections requested
345 
346  // Return exposure
347  return exposure;
348 }
349 
350 
351 /*==========================================================================
352  = =
353  = Public methods =
354  = =
355  ==========================================================================*/
356 
357 /***********************************************************************//**
358  * @brief Clear livetime cube
359  ***************************************************************************/
361 {
362  // Free class members
363  free_members();
364 
365  // Initialise members
366  init_members();
367 
368  // Return
369  return;
370 }
371 
372 
373 /***********************************************************************//**
374  * @brief Clone livetime cube
375  *
376  * @return Pointer to deep copy of livetime cube.
377  ***************************************************************************/
379 {
380  return new GLATLtCube(*this);
381 }
382 
383 
384 /***********************************************************************//**
385  * @brief Load livetime cube from FITS file
386  *
387  * @param[in] filename FITS file name.
388  ***************************************************************************/
389 void GLATLtCube::load(const GFilename& filename)
390 {
391  // Clear object
392  clear();
393 
394  // Open livetime cube FITS file
395  GFits fits(filename);
396 
397  // Read livetime cube
398  GLATLtCube::read(fits);
399 
400  // Close FITS file
401  fits.close();
402 
403  // Return
404  return;
405 }
406 
407 
408 /***********************************************************************//**
409  * @brief Save livetime cube into FITS file
410  *
411  * @param[in] filename FITS file name.
412  * @param[in] clobber Overwrite existing file? (default: false)
413  *
414  * @todo Not yet implemented.
415  ***************************************************************************/
416 void GLATLtCube::save(const GFilename& filename, const bool& clobber) const
417 {
418  // Create FITS file
419  GFits fits;
420 
421  // Write effective area into file
422  write(fits);
423 
424  // Close FITS file
425  fits.saveto(filename, clobber);
426 
427  // Return
428  return;
429 }
430 
431 
432 /***********************************************************************//**
433  * @brief Read livetime cube from FITS file
434  *
435  * @param[in] fits FITS.
436  *
437  * Reads a livetime cube from a FITS file.
438  *
439  * @todo Reading of cos theta boundaries not yet implemented. This is not
440  * critical since they are not really needed. We just need them once
441  * we want to implement also saving.
442  ***************************************************************************/
443 void GLATLtCube::read(const GFits& fits)
444 {
445  // Clear object
446  clear();
447 
448  // Get HDUs
449  const GFitsTable& hdu_exposure = *fits.table(gammalib::extname_lat_exposure);
450  const GFitsTable& hdu_weighted_exposure = *fits.table(gammalib::extname_lat_wgtexposure);
451  //const GFitsTable& hdu_cthetabounds = *fits.table(gammalib::extname_lat_cthetabounds);
452  const GFitsTable& hdu_gti = *fits.table(gammalib::extname_gti);
453 
454  // Load exposure
455  m_exposure.read(hdu_exposure);
456 
457  // Load weighted exposure
458  m_weighted_exposure.read(hdu_weighted_exposure);
459 
460  // Load cos theta boundaries
461 
462  // Load GTIs
463  m_gti.read(hdu_gti);
464 
465  // Return
466  return;
467 }
468 
469 
470 /***********************************************************************//**
471  * @brief Write livetime cube into FITS file
472  *
473  * @param[in] fits FITS file.
474  *
475  * Writes the livetime cube into a FITS file.
476  ***************************************************************************/
477 void GLATLtCube::write(GFits& fits) const
478 {
479  // Write exposure
481 
482  // Write weighted exposure
484 
485  // Write cos theta boundaries
486 
487  // Write GTIs
489 
490  // Return
491  return;
492 }
493 
494 
495 /***********************************************************************//**
496  * @brief Print livetime cube information
497  *
498  * @param[in] chatter Chattiness.
499  * @return String containing livetime cube information.
500  ***************************************************************************/
501 std::string GLATLtCube::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("=== GLATLtCube ===");
511 
512  // Append detailed information
513  GChatter reduced_chatter = gammalib::reduce(chatter);
514  if (reduced_chatter > SILENT) {
515  result.append("\n"+m_exposure.print(reduced_chatter));
516  result.append("\n"+m_weighted_exposure.print(reduced_chatter));
517  result.append("\n"+m_gti.print(reduced_chatter));
518  }
519 
520  } // endif: chatter was not silent
521 
522  // Return result
523  return result;
524 }
525 
526 
527 /*==========================================================================
528  = =
529  = Private methods =
530  = =
531  ==========================================================================*/
532 
533 /***********************************************************************//**
534  * @brief Initialise class members
535  ***************************************************************************/
537 {
538  // Initialise members
539  m_exposure.clear();
541  m_gti.clear();
542 
543  // Return
544  return;
545 }
546 
547 
548 /***********************************************************************//**
549  * @brief Copy class members
550  *
551  * @param[in] cube Livetime cube.
552  ***************************************************************************/
554 {
555  // Copy members
556  m_exposure = cube.m_exposure;
558  m_gti = cube.m_gti;
559 
560  // Return
561  return;
562 }
563 
564 
565 /***********************************************************************//**
566  * @brief Delete class members
567  ***************************************************************************/
569 {
570  // Return
571  return;
572 }
double(* _ltcube_ctheta_phi)(const double &costheta, const double &phi)
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
void read(const GFits &file)
Read livetime cube from FITS file.
Definition: GLATLtCube.cpp:443
Energy value class definition.
double operator()(const GSkyDir &dir, const GEnergy &energy, _ltcube_ctheta fct) const
Sum function multiplied by efficiency corrected livetime over zenith angle.
Definition: GLATLtCube.cpp:172
void load(const GFilename &filename)
Load livetime cube from FITS file.
Definition: GLATLtCube.cpp:389
Sky direction class interface definition.
void write(GFits &file) const
Write livetime cube into FITS file.
Definition: GLATLtCube.cpp:477
GLATLtCubeMap m_weighted_exposure
Definition: GLATLtCube.hpp:99
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
const std::string extname_gti
Definition: GGti.hpp:44
GLATLtCube * clone(void) const
Clone livetime cube.
Definition: GLATLtCube.cpp:378
std::string print(const GChatter &chatter=NORMAL) const
Print livetime cube information.
Definition: GLATLtCube.cpp:501
std::string print(const GChatter &chatter=NORMAL) const
Print Good Time Intervals.
Definition: GGti.cpp:1106
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
const std::string extname_lat_exposure
Definition: GLATLtCube.hpp:45
void clear(void)
Clear Good Time Intervals.
Definition: GGti.cpp:237
Interface for the Fermi LAT point spread function.
Definition: GLATPsf.hpp:54
void write(GFits &fits, const std::string &extname=gammalib::extname_lat_ltcubemap) const
Write livetime cube map into FITS file.
Filename class.
Definition: GFilename.hpp:62
Fermi/LAT livetime cube class definition.
void save(const GFilename &filename, const bool &clobber=false) const
Save livetime cube into FITS file.
Definition: GLATLtCube.cpp:416
virtual ~GLATLtCube(void)
Destructor.
Definition: GLATLtCube.cpp:108
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
Fermi LAT point spread function class definition.
void init_members(void)
Initialise class members.
Definition: GLATLtCube.cpp:536
double efficiency_factor1(const GEnergy &srcEng) const
Returns efficiency factor 1.
Definition: GLATAeff.cpp:403
GLATLtCube(void)
Void constructor.
Definition: GLATLtCube.cpp:59
void clear(void)
Clear livetime cube map.
void free_members(void)
Delete class members.
Definition: GLATLtCube.cpp:568
GLATLtCubeMap m_exposure
Definition: GLATLtCube.hpp:98
void read(const GFitsTable &table)
Read Good Time Intervals and time reference from FITS table.
Definition: GGti.cpp:753
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
void copy_members(const GLATLtCube &cube)
Copy class members.
Definition: GLATLtCube.cpp:553
const std::string extname_lat_wgtexposure
Definition: GLATLtCube.hpp:46
double(* _ltcube_ctheta)(const double &costheta)
Fermi LAT effective area class definition.
GLATLtCube & operator=(const GLATLtCube &cube)
Assignment operator.
Definition: GLATLtCube.cpp:130
Sky direction class.
Definition: GSkyDir.hpp:62
void write(GFits &fits, const std::string &extname=gammalib::extname_gti) const
Write Good Time Intervals and time reference into FITS object.
Definition: GGti.cpp:796
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
void clear(void)
Clear livetime cube.
Definition: GLATLtCube.cpp:360
void read(const GFitsTable &table)
Load livetime cube from FITS file.
bool has_efficiency(void) const
Signals whether efficiency factors are present.
Definition: GLATAeff.cpp:384
Interface for the Fermi LAT livetime cube.
Definition: GLATLtCube.hpp:59
Filename class interface definition.
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
double efficiency_factor2(const GEnergy &srcEng) const
Returns efficiency factor 2.
Definition: GLATAeff.cpp:432