GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GLATMeanPsf.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GLATMeanPsf.cpp - Fermi/LAT mean PSF class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2016 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 GLATMeanPsf.cpp
23  * @brief Fermi/LAT mean PSF class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GLATMeanPsf.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GLATAeff.hpp"
36 #include "GLATPsf.hpp"
37 #include "GLATObservation.hpp"
38 #include "GLATEventCube.hpp"
39 #include "GLATException.hpp"
40 
41 /* __ Method name definitions ____________________________________________ */
42 #define G_SET "GLATMeanPsf::set(GSkyDir&, GLATObservation&)"
43 #define G_EXPOSURE "GLATMeanPsf::exposure(int&)"
44 
45 /* __ Macros _____________________________________________________________ */
46 
47 /* __ Coding definitions _________________________________________________ */
48 
49 /* __ Debug definitions __________________________________________________ */
50 #define G_SIGNAL_NEGATIVE_MEAN_PSF //!< Signal negative mean PSF value
51 //#define G_DEBUG_INTEGRAL //!< Debug mean PSF integration
52 //#define G_DEBUG_MAP_CORRECTIONS //!< Debug map corrections
53 
54 /* __ Constants __________________________________________________________ */
55 
56 /* __ Typedefs ___________________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  ***************************************************************************/
69 {
70  // Initialise class members
71  init_members();
72 
73  // Return
74  return;
75 }
76 
77 
78 /***********************************************************************//**
79  * @brief PSF constructor
80  *
81  * @param[in] dir Sky direction.
82  * @param[in] obs LAT observation.
83  ***************************************************************************/
85 {
86  // Initialise class members
87  init_members();
88 
89  // Setup PSF
90  set(dir, obs);
91 
92  // Return
93  return;
94 }
95 
96 
97 /***********************************************************************//**
98  * @brief Copy constructor
99  *
100  * @param[in] psf Mean PSF.
101  ***************************************************************************/
103 {
104  // Initialise class members
105  init_members();
106 
107  // Copy members
108  copy_members(psf);
109 
110  // Return
111  return;
112 }
113 
114 
115 /***********************************************************************//**
116  * @brief Destructor
117  ***************************************************************************/
119 {
120  // Free members
121  free_members();
122 
123  // Return
124  return;
125 }
126 
127 
128 /*==========================================================================
129  = =
130  = Operators =
131  = =
132  ==========================================================================*/
133 
134 /***********************************************************************//**
135  * @brief Assignment operator
136  *
137  * @param[in] psf Mean PSF.
138  * @return Mean PSF.
139  ***************************************************************************/
141 {
142  // Execute only if object is not identical
143  if (this != &psf) {
144 
145  // Free members
146  free_members();
147 
148  // Initialise private members
149  init_members();
150 
151  // Copy members
152  copy_members(psf);
153 
154  } // endif: object was not identical
155 
156  // Return this object
157  return *this;
158 }
159 
160 
161 /***********************************************************************//**
162  * @brief Return mean PSF*exposure value
163  *
164  * @param[in] offset Angular distance from true source direction (degrees).
165  * @param[in] logE log10 of energy in MeV.
166  *
167  * Computes
168  * \f[{\rm PSF}(\delta, \log E) \times {\rm exposure}(\log E)\f]
169  * by bilinear interpolation, where
170  * \f$\delta\f$ is the offset angle from the source direction in degrees, and
171  * \f$\log E\f$ is the logarithm of base 10 of the energy in MeV.
172  * A zero value is returned if the offset angle is equal or larger than
173  * 70 degrees or if \f$\log E\f$ is not positive.
174  ***************************************************************************/
175 double GLATMeanPsf::operator() (const double& offset, const double& logE)
176 {
177  // Initialise response
178  double value = 0.0;
179 
180  // Continue only if arguments are within valid range
181  if (offset < 70.0 && logE > 0.0) {
182 
183  // Flag no change of values
184  //bool change = false;
185 
186  // Set offset interpolation
187  //if (offset != m_last_offset) {
188  m_offset.set_value(offset);
189  // m_last_offset = offset;
190  // change = true;
191  //}
192 
193  // Set energy interpolation
194  //if (logE != m_last_energy) {
195  m_energy.set_value(logE);
196  // m_last_energy = logE;
197  // change = true;
198  //}
199 
200  // If change occured then update interpolation indices and weighting
201  // factors
202  //if (change) {
203 
204  // Set energy indices for exposure computation
207 
208  // Set energy indices for PSF computation
209  int inx_energy_left = m_inx1_exp * noffsets();
210  int inx_energy_right = m_inx2_exp * noffsets();
211 
212  // Set array indices for bi-linear interpolation
213  m_inx1 = m_offset.inx_left() + inx_energy_left;
214  m_inx2 = m_offset.inx_left() + inx_energy_right;
215  m_inx3 = m_offset.inx_right() + inx_energy_left;
216  m_inx4 = m_offset.inx_right() + inx_energy_right;
217 
218  // Set weighting factors for bi-linear interpolation
223 
224  //} // endif: logE or ctheta changed
225 
226  // Compute energy dependent exposure and map corrections
227  double fac_left = m_exposure[m_inx1_exp] * m_mapcorr[m_inx1_exp];
228  double fac_right = m_exposure[m_inx2_exp] * m_mapcorr[m_inx2_exp];
229 
230  // Perform bi-linear interpolation
231  value = m_wgt1 * m_psf[m_inx1] * fac_left +
232  m_wgt2 * m_psf[m_inx2] * fac_right +
233  m_wgt3 * m_psf[m_inx3] * fac_left +
234  m_wgt4 * m_psf[m_inx4] * fac_right;
235 
236  // Optionally check for negative values
237  #if defined(G_SIGNAL_NEGATIVE_MEAN_PSF)
238  if (value < 0.0) {
239  std::cout << "GLATMeanPsf::operator()(offset=" << offset;
240  std::cout << ", logE=" << logE << ") = " << value << std::endl;
241  }
242  #endif
243 
244  } // endif: arguments were valid
245 
246  // Return value
247  return value;
248 }
249 
250 
251 /*==========================================================================
252  = =
253  = Public methods =
254  = =
255  ==========================================================================*/
256 
257 /***********************************************************************//**
258  * @brief Clear mean PSF
259  ***************************************************************************/
261 {
262  // Free class members
263  free_members();
264 
265  // Initialise members
266  init_members();
267 
268  // Return
269  return;
270 }
271 
272 
273 /***********************************************************************//**
274  * @brief Clone mean PSF
275  *
276  * @return Pointer to deep copy of mean PSF.
277  ***************************************************************************/
279 {
280  return new GLATMeanPsf(*this);
281 }
282 
283 
284 /***********************************************************************//**
285  * @brief Compute mean PSF and exposure
286  *
287  * @param[in] dir Source location.
288  * @param[in] obs LAT observation.
289  *
290  * @exception GLATException::no_ltcube
291  * Livetime cube has not been defined.
292  *
293  * Computes the mean PSF and the energy dependent exposure for a source at
294  * a given sky location. The PSF is computed for all bin boundaries of the
295  * observation, hence if there are N energy bins there will be N+1 energies
296  * at which the mean PSF is computed.
297  ***************************************************************************/
298 void GLATMeanPsf::set(const GSkyDir& dir, const GLATObservation& obs)
299 {
300  // Clear PSF, exposure and energy arrays
301  m_psf.clear();
302  m_exposure.clear();
303  m_energy.clear();
304 
305  // Get pointers on response, livetime cube and energy boundaries
306  const GLATResponse* rsp = obs.response();
307 
308  // Get pointer on livetime cube
309  const GLATLtCube* ltcube = obs.ltcube();
310  if (ltcube == NULL) {
312  }
313 
314  // Get energy boundaries
315  GEbounds ebds = obs.events()->ebounds();
316 
317  // Store source direction
318  m_dir = dir;
319 
320  // Limit computation to zenith angles < m_theta_max (typically 70
321  // degrees - this is the hardwired value in the ST). For this purpose
322  // set the costhetamin parameter of Aeff to m_theta_max. Store also the
323  // original values in a vector for later restoration.
324  /*
325  std::vector<double> save_costhetamin;
326  for (int i = 0; i < rsp.size(); ++i) {
327  save_costhetamin.push_back(rsp.aeff(i)->costhetamin());
328  rsp.aeff(i)->costhetamin(cos(m_theta_max*gammalib::deg2rad));
329  }
330  */
331 
332  // Allocate room for arrays
333  m_psf.reserve(size());
334  m_exposure.reserve(m_energy.size());
335 
336  // Set energy nodes from the bin boundaries of the observations energy
337  // boundaries. Store the energy nodes locally as GEnergy objects and
338  // save them also in the class as log10 of energy in MeV.
339  std::vector<GEnergy> energy;
340  energy.reserve(ebds.size()+1);
341  energy.push_back(ebds.emin(0));
342  m_energy.append(ebds.emin(0).log10MeV());
343  for (int i = 0; i < ebds.size(); ++i) {
344  m_energy.append(ebds.emax(i).log10MeV());
345  energy.push_back(ebds.emax(i));
346  }
347 
348  // Loop over energies
349  for (int ieng = 0; ieng < energy.size(); ++ieng) {
350 
351  // Compute exposure by looping over the responses
352  double exposure = 0.0;
353  for (int i = 0; i < rsp->size(); ++i) {
354  exposure += (*ltcube)(dir, energy[ieng], *rsp->aeff(i));
355  }
356 
357  // Set exposure
358  m_exposure.push_back(exposure);
359 
360  // Loop over all offset angles
361  for (int ioffset = 0; ioffset < m_offset.size(); ++ioffset) {
362 
363  // Compute point spread function by looping over the responses
364  double psf = 0.0;
365  for (int i = 0; i < rsp->size(); ++i) {
366  psf += (*ltcube)(dir, energy[ieng], m_offset[ioffset],
367  *rsp->psf(i), *rsp->aeff(i));
368  }
369 
370  // Normalize PSF by exposure and clip when exposure drops to 0
371  psf = (exposure > 0.0) ? psf/exposure : 0.0;
372 
373  // Set PSF value
374  m_psf.push_back(psf);
375 
376  } // endfor: looped over offsets
377  } // endfor: looped over energies
378 
379  // Restore initial Aeff zenith angle restriction
380  /*
381  for (int i = 0; i < rsp.size(); ++i) {
382  rsp.aeff(i)->costhetamin(save_costhetamin[i]);
383  }
384  */
385 
386  // Compute map corrections
387  set_map_corrections(obs);
388 
389  // Return
390  return;
391 }
392 
393 
394 /***********************************************************************//**
395  * @brief Return mean PSF value
396  *
397  * @param[in] offset Angular distance from true source direction (degrees).
398  * @param[in] logE log10 of energy in MeV.
399  *
400  * Computes
401  * \f[{\rm PSF}(\delta, \log E)\f]
402  * by bilinear interpolation, where
403  * \f$\delta\f$ is the offset angle from the source direction in degrees, and
404  * \f$\log E\f$ is the logarithm of base 10 of the energy in MeV.
405  * A zero value is returned if the offset angle is equal or larger than
406  * 70 degrees or if \f$\log E\f$ is not positive.
407  ***************************************************************************/
408 double GLATMeanPsf::psf(const double& offset, const double& logE)
409 {
410  // Initialise response
411  double value = 0.0;
412 
413  // Continue only if arguments are within valid range
414  if (offset < 70.0 && logE > 0.0) {
415 
416  // Flag no change of values
417  //bool change = false;
418 
419  // Set offset interpolation
420  //if (offset != m_last_offset) {
421  m_offset.set_value(offset);
422  // m_last_offset = offset;
423  // change = true;
424  //}
425 
426  // Set energy interpolation
427  //if (logE != m_last_energy) {
428  m_energy.set_value(logE);
429  // m_last_energy = logE;
430  // change = true;
431  //}
432 
433  // If change occured then update interpolation indices and weighting
434  // factors
435  //if (change) {
436 
437  // Set energy indices for exposure computation
440 
441  // Set energy indices
442  int inx_energy_left = m_energy.inx_left() * noffsets();
443  int inx_energy_right = m_energy.inx_right() * noffsets();
444 
445  // Set array indices for bi-linear interpolation
446  m_inx1 = m_offset.inx_left() + inx_energy_left;
447  m_inx2 = m_offset.inx_left() + inx_energy_right;
448  m_inx3 = m_offset.inx_right() + inx_energy_left;
449  m_inx4 = m_offset.inx_right() + inx_energy_right;
450 
451  // Set weighting factors for bi-linear interpolation
456 
457  //} // endif: logE or ctheta changed
458 
459  // Compute energy dependentmap corrections
460  double fac_left = m_mapcorr[m_inx1_exp];
461  double fac_right = m_mapcorr[m_inx2_exp];
462 
463  // Perform bi-linear interpolation
464  value = m_wgt1 * m_psf[m_inx1] * fac_left +
465  m_wgt2 * m_psf[m_inx2] * fac_right +
466  m_wgt3 * m_psf[m_inx3] * fac_left +
467  m_wgt4 * m_psf[m_inx4] * fac_right;
468 
469  // Optionally check for negative values
470  #if defined(G_SIGNAL_NEGATIVE_MEAN_PSF)
471  if (value < 0.0) {
472  std::cout << "GLATMeanPsf::psf(offset=" << offset;
473  std::cout << ", logE=" << logE << ") = " << value << std::endl;
474  }
475  #endif
476 
477  } // endif: arguments were valid
478 
479  // Return value
480  return value;
481 }
482 
483 
484 /***********************************************************************//**
485  * @brief Return exposure value
486  *
487  * @param[in] logE log10 of energy in MeV.
488  *
489  * Computes
490  * \f[{\rm exposure}(\log E)\f]
491  * by linear interpolation, where
492  * \f$\log E\f$ is the logarithm of base 10 of the energy in MeV.
493  * A zero value is returned if \f$\log E\f$ is not positive.
494  ***************************************************************************/
495 double GLATMeanPsf::exposure(const double& logE)
496 {
497  // Initialise response
498  double value = 0.0;
499 
500  // Continue only if arguments are within valid range
501  if (logE > 0.0) {
502 
503  // Set energy interpolation
504  //if (logE != m_last_energy) {
505  m_energy.set_value(logE);
506  // m_last_energy = logE;
509  //}
510 
511  // Perform linear interpolation
512  value = m_energy.wgt_left() * m_exposure[m_inx1_exp] +
514 
515  } // endif: arguments were in valid range
516 
517  // Return value
518  return value;
519 }
520 
521 
522 /***********************************************************************//**
523  * @brief Print livetime cube information
524  *
525  * @param[in] chatter Chattiness (defaults to NORMAL).
526  * @return String containing ROI information.
527  ***************************************************************************/
528 std::string GLATMeanPsf::print(const GChatter& chatter) const
529 {
530  // Initialise result string
531  std::string result;
532 
533  // Continue only if chatter is not silent
534  if (chatter != SILENT) {
535 
536  // Compute exposure range
537  double min_exposure = 0.0;
538  double max_exposure = 0.0;
539  for (int i = 0; i < nenergies(); ++i) {
540  if (i == 0) {
541  min_exposure = m_exposure[i];
542  max_exposure = m_exposure[i];
543  }
544  else {
545  if (m_exposure[i] < min_exposure) {
546  min_exposure = m_exposure[i];
547  }
548  if (m_exposure[i] > max_exposure) {
549  max_exposure = m_exposure[i];
550  }
551  }
552  }
553 
554  // Append header
555  result.append("=== GLATMeanPsf ===");
556 
557  // Append information
558  result.append("\n"+gammalib::parformat("Source name")+name());
559  result.append("\n"+gammalib::parformat("Source direction"));
560  result.append(gammalib::str(m_dir.ra_deg()));
561  result.append(", ");
562  result.append(gammalib::str(m_dir.dec_deg()));
563  result.append("\n"+gammalib::parformat("Offset angles")+gammalib::str(noffsets()));
564  result.append("\n"+gammalib::parformat("Energy values")+gammalib::str(nenergies()));
565  result.append("\n"+gammalib::parformat("Exposure range"));
566  result.append(gammalib::str(min_exposure)+" - "+gammalib::str(max_exposure)+" s cm2");
567  for (int i = 0; i < nenergies(); ++i) {
568  GEnergy energy;
569  energy.log10MeV(m_energy[i]);
570  result.append("\n"+gammalib::parformat(energy.print()));
571  result.append(gammalib::str(m_exposure[i])+" s cm2");
572  if (m_mapcorr.size() == nenergies()) {
573  result.append(" (map correction="+gammalib::str(m_mapcorr[i])+")");
574  }
575  }
576 
577  } // endif: chatter was not silent
578 
579  // Return result
580  return result;
581 }
582 
583 
584 /*==========================================================================
585  = =
586  = Private methods =
587  = =
588  ==========================================================================*/
589 
590 /***********************************************************************//**
591  * @brief Initialise class members
592  ***************************************************************************/
594 {
595  // Initialise members
596  m_name.clear();
597  m_dir.clear();
598  m_psf.clear();
599  m_exposure.clear();
600  m_mapcorr.clear();
601  m_energy.clear();
602  m_offset.clear();
603  m_theta_max = 70.0; //!< Maximum zenith angle
604  m_last_energy = -1.0;
605  m_last_offset = -1.0;
606  m_inx1_exp = 0;
607  m_inx2_exp = 0;
608  m_inx1 = 0;
609  m_inx2 = 0;
610  m_inx3 = 0;
611  m_inx4 = 0;
612  m_wgt1 = 0.0;
613  m_wgt2 = 0.0;
614  m_wgt3 = 0.0;
615  m_wgt4 = 0.0;
616 
617  // Set offset array
618  set_offsets();
619 
620  // Return
621  return;
622 }
623 
624 
625 /***********************************************************************//**
626  * @brief Copy class members
627  *
628  * @param[in] psf Mean PSF.
629  ***************************************************************************/
631 {
632  // Copy members
633  m_name = psf.m_name;
634  m_dir = psf.m_dir;
635  m_psf = psf.m_psf;
636  m_exposure = psf.m_exposure;
637  m_mapcorr = psf.m_mapcorr;
638  m_energy = psf.m_energy;
639  m_offset = psf.m_offset;
640  m_theta_max = psf.m_theta_max;
643  m_inx1_exp = psf.m_inx1_exp;
644  m_inx2_exp = psf.m_inx2_exp;
645  m_inx1 = psf.m_inx1;
646  m_inx2 = psf.m_inx2;
647  m_inx3 = psf.m_inx3;
648  m_inx4 = psf.m_inx4;
649  m_wgt1 = psf.m_wgt1;
650  m_wgt2 = psf.m_wgt2;
651  m_wgt3 = psf.m_wgt3;
652  m_wgt4 = psf.m_wgt4;
653 
654  // Return
655  return;
656 }
657 
658 
659 /***********************************************************************//**
660  * @brief Delete class members
661  ***************************************************************************/
663 {
664  // Return
665  return;
666 }
667 
668 
669 /***********************************************************************//**
670  * @brief Set array of offset values in degrees
671  *
672  * The array of offset values defines the points at with the mean PSF will
673  * be evaluated and stored. The array is logarithmically spaced.
674  ***************************************************************************/
676 {
677  // Set array parameters
678  const double offset_min = 1.0e-4;
679  const double offset_max = 70.0;
680  const int offset_num = 200;
681 
682  // Clear offset array
683  m_offset.clear();
684 
685  // Set array
686  double step = log(offset_max/offset_min)/(offset_num - 1.0);
687  for (int i = 0; i < offset_num; ++i)
688  m_offset.append(offset_min*exp(i*step));
689 
690  // Return
691  return;
692 }
693 
694 
695 /***********************************************************************//**
696  * @brief Compute map corrections
697  *
698  * @param[in] obs LAT observation.
699  *
700  * Sets up a vector of energy dependent corrections that assures that the
701  * integral over the PSF in a specific event cube is correctly normalised.
702  * For this purpose the PSF is integrated over a radius that fully lies
703  * within the event cube and devided by the PSF pixel sum within this radius.
704  * The map corrections are mainly important at high energies where the
705  * PSF size can become comparable to the pixel size.
706  *
707  * See SourceMap::getMapCorrections for the original code in the Fermi/LAT
708  * ScienceTools.
709  *
710  * @todo We can also implement a method for event atoms, yet it is not clear
711  * whether we really need this.
712  ***************************************************************************/
714 {
715  // Initialise map corrections to unity vector
716  m_mapcorr.assign(m_energy.size(), 1.0);
717 
718  // Get pointer on event cube
719  const GLATEventCube* cube = dynamic_cast<const GLATEventCube*>(obs.events());
720 
721  // Continue only if we have an event cube
722  if (cube != NULL) {
723 
724  // Get maximum PSF radius
725  double radius = cube->maxrad(m_dir);
726  #if defined(G_DEBUG_MAP_CORRECTIONS)
727  std::cout << "GLATMeanPsf::set_map_corrections:";
728  std::cout << " radius=" << radius << std::endl;
729  #endif
730  if (radius > 0.0) {
731 
732  // Initialise sum over pixels
733  std::vector<double> sum(m_energy.size(), 0.0);
734 
735  // Loop over all event cube spatial pixels
736  for (int iy = 0; iy < cube->ny(); ++iy) {
737  for (int ix = 0; ix < cube->nx(); ++ix) {
738 
739  // Compute offset angle in degrees
740  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
741  double offset = cube->map().pix2dir(pixel).dist_deg(m_dir);
742  double solidangle = cube->map().solidangle(pixel);
743 
744  // Use only pixels within maximum PSF radius
745  if (offset <= radius) {
746 
747  // Accumulate energy dependent pixel sum
748  for (int ieng = 0; ieng < m_energy.size(); ++ieng) {
749  sum[ieng] += psf(offset, m_energy[ieng]) * solidangle;
750  }
751 
752  } // endif: pixel was within maximum PSF radius
753 
754  } // endfor: looped over latitude pixels
755  } // endfor: looped over longitude pixels
756 
757  // Compute map correction
758  for (int ieng = 0; ieng < m_energy.size(); ++ieng) {
759  if (sum[ieng] > 0.0) {
760  m_mapcorr[ieng] = integral(radius, m_energy[ieng]) / sum[ieng];
761  }
762  }
763 
764  } // endif: radius was positive
765  } // endif: observation had an event cube
766 
767  // Return
768  return;
769 }
770 
771 
772 /***********************************************************************//**
773  * @brief Compute integral over PSF
774  *
775  * @param[in] offsetmax Maximum offset angle.
776  * @param[in] logE log10 of energy in MeV.
777  ***************************************************************************/
778 double GLATMeanPsf::integral(const double& offsetmax, const double& logE)
779 {
780  // Set energy and offset interpolation weigths
781  m_energy.set_value(logE);
782  m_offset.set_value(offsetmax);
783 
784  // Get PSF array offsets
785  int inx_energy_left = m_energy.inx_left() * noffsets();
786  int inx_energy_right = m_energy.inx_right() * noffsets();
787 
788  // Initialise integrals
789  double int_left = 0.0;
790  double int_right = 0.0;
791 
792  // Loop over all offset angles
793  for (int i = 0; i < m_offset.size()-1; ++i) {
794 
795  // If interval is fully contained within offset then simply use
796  // tabulated values for integration
797  if (m_offset[i+1] <= offsetmax) {
798  double theta_min = m_offset[i] * gammalib::deg2rad;
799  double theta_max = m_offset[i+1] * gammalib::deg2rad;
800  int_left += 0.5 * (m_psf[inx_energy_left+i] * sin(theta_min) +
801  m_psf[inx_energy_left+i+1] * sin(theta_max)) *
802  (theta_max - theta_min);
803  int_right += 0.5 * (m_psf[inx_energy_right+i] * sin(theta_min) +
804  m_psf[inx_energy_right+i+1] * sin(theta_max)) *
805  (theta_max - theta_min);
806  }
807 
808  // ... otherwise interpolate the PSF value for theta_max. We can
809  // then exit the loop since we're done with the integration
810  else {
811  double theta_min = m_offset[i] * gammalib::deg2rad;
812  double theta_max = offsetmax * gammalib::deg2rad;
813  double psf_left = m_psf[inx_energy_left+i] * m_offset.wgt_left() +
814  m_psf[inx_energy_left+i+1] * m_offset.wgt_right();
815  double psf_right = m_psf[inx_energy_right+i] * m_offset.wgt_left() +
816  m_psf[inx_energy_right+i+1] * m_offset.wgt_right();
817  int_left += 0.5 * (m_psf[inx_energy_left+i] * sin(theta_min) +
818  psf_left * sin(theta_max)) *
819  (theta_max - theta_min);
820  int_right += 0.5 * (m_psf[inx_energy_right+i] * sin(theta_min) +
821  psf_right * sin(theta_max)) *
822  (theta_max - theta_min);
823 
824  // Debug option: Dump integral computation results
825  #if defined(G_DEBUG_INTEGRAL)
826  std::cout << "offsetmax=" << offsetmax;
827  std::cout << " offset.left=" << m_offset.inx_left();
828  std::cout << " offset.right=" << m_offset.inx_right();
829  std::cout << " i=" << i;
830  std::cout << " psf_left=" << psf_left;
831  std::cout << " psf_right=" << psf_right;
832  std::cout << " int_left=" << int_left*gammalib::twopi;
833  std::cout << " int_right=" << int_right*gammalib::twopi;
834  #endif
835  break;
836  }
837 
838  } // endfor: looped over offset angles
839 
840  // Interpolate now in energy
841  double integral = gammalib::twopi * (m_energy.wgt_left() * int_left +
842  m_energy.wgt_right() * int_right);
843 
844  // Debug option: Dump integral computation results
845  #if defined(G_DEBUG_INTEGRAL)
846  std::cout << " integral=" << integral << std::endl;
847  #endif
848 
849  // Return integral
850  return integral;
851 }
int m_inx3
Index 3.
Fermi/LAT observation class.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
std::string m_name
Source name for mean PSF.
Definition: GLATMeanPsf.hpp:93
const GSkyDir & dir(void) const
Return sky direction for mean PSF.
virtual void response(const GResponse &rsp)
Set response function.
Fermi/LAT observation class definition.
double m_theta_max
Maximum inclination angle (default 70 deg)
int size(void) const
Return number of event types.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:250
double operator()(const double &offset, const double &logE)
Return mean PSF*exposure value.
const std::string & name(void) const
Return source name for mean PSF.
double m_wgt4
Weighting factor 4.
void init_members(void)
Initialise class members.
const double & offset(const int &inx) const
Return offset angle for given bin.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
int m_inx4
Index 4.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
int m_inx2_exp
Exposure index 2.
double maxrad(const GSkyDir &dir) const
Computes the maximum radius (in degrees) around a given source direction that fits spatially into the...
GLATMeanPsf * clone(void) const
Clone mean PSF.
GLATPsf * psf(const int &index) const
Return pointer on point spread function.
int nx(void) const
Return number of bins in X direction.
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
Fermi/LAT Response class.
void free_members(void)
Delete class members.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
void set_offsets(void)
Set array of offset values in degrees.
Gammalib tools definition.
int size(void) const
Return number of bins in mean PSF.
#define G_SET
Definition: GLATMeanPsf.cpp:42
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
int m_inx2
Index 2.
GLATMeanPsf & operator=(const GLATMeanPsf &cube)
Assignment operator.
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
double exposure(const double &logE)
Return exposure value.
GLATMeanPsf(void)
Void constructor.
Definition: GLATMeanPsf.cpp:68
Sky map pixel class.
Definition: GSkyPixel.hpp:74
void clear(void)
Clear mean PSF.
const double deg2rad
Definition: GMath.hpp:43
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
double m_wgt1
Weighting factor 1.
Energy boundaries container class.
Definition: GEbounds.hpp:60
const GLATLtCube * ltcube(void) const
Return Fermi/LAT livetime cube.
GSkyDir m_dir
Source direction for mean PSF.
Definition: GLATMeanPsf.hpp:94
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
LAT exception handler interface definition.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
double m_last_energy
Last requested logE value.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:223
GChatter
Definition: GTypemaps.hpp:33
Fermi/LAT event cube class definition.
std::vector< double > m_mapcorr
Map corrections.
Definition: GLATMeanPsf.hpp:97
double m_last_offset
Last requested offset value.
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
Fermi LAT point spread function class definition.
Fermi/LAT event cube class.
GNodeArray m_offset
Offsets of mean PSF.
Definition: GLATMeanPsf.hpp:98
GLATAeff * aeff(const int &index) const
Return pointer on effective area.
void set(const GSkyDir &dir, const GLATObservation &obs)
Compute mean PSF and exposure.
double m_wgt3
Weighting factor 3.
std::vector< double > m_psf
Mean PSF values.
Definition: GLATMeanPsf.hpp:95
Fermi/LAT mean PSF class definition.
double psf(const double &offset, const double &logE)
Return mean PSF value.
double m_wgt2
Weighting factor 2.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
void map(const GSkyMap &map)
Set event cube from sky map.
void copy_members(const GLATMeanPsf &psf)
Copy class members.
double integral(const double &radmax, const double &logE)
Compute integral over PSF.
std::vector< double > m_exposure
Mean exposure.
Definition: GLATMeanPsf.hpp:96
void clear(void)
Clear sky direction.
Definition: GSkyDir.cpp:142
int nenergies(void) const
Return number of energy bins.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
virtual GEvents * events(void)
Return events.
Fermi/LAT mean PSF class.
Definition: GLATMeanPsf.hpp:51
std::string print(const GChatter &chatter=NORMAL) const
Print livetime cube information.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
int noffsets(void) const
Return number of offset bins.
const double twopi
Definition: GMath.hpp:36
GNodeArray m_energy
log10(energy) of mean PSF
Definition: GLATMeanPsf.hpp:99
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
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
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
int m_inx1_exp
Exposure index 1.
virtual ~GLATMeanPsf(void)
Destructor.
Sky direction class.
Definition: GSkyDir.hpp:62
const double & energy(const int &inx) const
Return energy for given bin.
void set_map_corrections(const GLATObservation &obs)
Compute map corrections.
int m_inx1
Index 1.
int ny(void) const
Return number of bins in Y direction.
Interface for the Fermi LAT livetime cube.
Definition: GLATLtCube.hpp:59
Mathematical function definitions.
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