GammaLib  2.0.0
 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-2021 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file 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 "GException.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 "GLATMeanPsf.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 GException::invalid_argument
291  * No livetime cube found in observation.
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) {
311  std::string msg = "Observation does not contain a livetime cube. "
312  "Please specify an observation containing a livetime "
313  "cube.";
315  }
316 
317  // Get energy boundaries
318  GEbounds ebds = obs.events()->ebounds();
319 
320  // Store source direction
321  m_dir = dir;
322 
323  // Limit computation to zenith angles < m_theta_max (typically 70
324  // degrees - this is the hardwired value in the ST). For this purpose
325  // set the costhetamin parameter of Aeff to m_theta_max. Store also the
326  // original values in a vector for later restoration.
327  /*
328  std::vector<double> save_costhetamin;
329  for (int i = 0; i < rsp.size(); ++i) {
330  save_costhetamin.push_back(rsp.aeff(i)->costhetamin());
331  rsp.aeff(i)->costhetamin(cos(m_theta_max*gammalib::deg2rad));
332  }
333  */
334 
335  // Allocate room for arrays
336  m_psf.reserve(size());
337  m_exposure.reserve(m_energy.size());
338 
339  // Set energy nodes from the bin boundaries of the observations energy
340  // boundaries. Store the energy nodes locally as GEnergy objects and
341  // save them also in the class as log10 of energy in MeV.
342  std::vector<GEnergy> energy;
343  energy.reserve(ebds.size()+1);
344  energy.push_back(ebds.emin(0));
345  m_energy.append(ebds.emin(0).log10MeV());
346  for (int i = 0; i < ebds.size(); ++i) {
347  m_energy.append(ebds.emax(i).log10MeV());
348  energy.push_back(ebds.emax(i));
349  }
350 
351  // Loop over energies
352  for (int ieng = 0; ieng < energy.size(); ++ieng) {
353 
354  // Compute exposure by looping over the responses
355  double exposure = 0.0;
356  for (int i = 0; i < rsp->size(); ++i) {
357  exposure += (*ltcube)(dir, energy[ieng], *rsp->aeff(i));
358  }
359 
360  // Set exposure
361  m_exposure.push_back(exposure);
362 
363  // Loop over all offset angles
364  for (int ioffset = 0; ioffset < m_offset.size(); ++ioffset) {
365 
366  // Compute point spread function by looping over the responses
367  double psf = 0.0;
368  for (int i = 0; i < rsp->size(); ++i) {
369  psf += (*ltcube)(dir, energy[ieng], m_offset[ioffset],
370  *rsp->psf(i), *rsp->aeff(i));
371  }
372 
373  // Normalize PSF by exposure and clip when exposure drops to 0
374  psf = (exposure > 0.0) ? psf/exposure : 0.0;
375 
376  // Set PSF value
377  m_psf.push_back(psf);
378 
379  } // endfor: looped over offsets
380  } // endfor: looped over energies
381 
382  // Restore initial Aeff zenith angle restriction
383  /*
384  for (int i = 0; i < rsp.size(); ++i) {
385  rsp.aeff(i)->costhetamin(save_costhetamin[i]);
386  }
387  */
388 
389  // Compute map corrections
390  set_map_corrections(obs);
391 
392  // Return
393  return;
394 }
395 
396 
397 /***********************************************************************//**
398  * @brief Return mean PSF value
399  *
400  * @param[in] offset Angular distance from true source direction (degrees).
401  * @param[in] logE log10 of energy in MeV.
402  *
403  * Computes
404  * \f[{\rm PSF}(\delta, \log E)\f]
405  * by bilinear interpolation, where
406  * \f$\delta\f$ is the offset angle from the source direction in degrees, and
407  * \f$\log E\f$ is the logarithm of base 10 of the energy in MeV.
408  * A zero value is returned if the offset angle is equal or larger than
409  * 70 degrees or if \f$\log E\f$ is not positive.
410  ***************************************************************************/
411 double GLATMeanPsf::psf(const double& offset, const double& logE)
412 {
413  // Initialise response
414  double value = 0.0;
415 
416  // Continue only if arguments are within valid range
417  if (offset < 70.0 && logE > 0.0) {
418 
419  // Flag no change of values
420  //bool change = false;
421 
422  // Set offset interpolation
423  //if (offset != m_last_offset) {
424  m_offset.set_value(offset);
425  // m_last_offset = offset;
426  // change = true;
427  //}
428 
429  // Set energy interpolation
430  //if (logE != m_last_energy) {
431  m_energy.set_value(logE);
432  // m_last_energy = logE;
433  // change = true;
434  //}
435 
436  // If change occured then update interpolation indices and weighting
437  // factors
438  //if (change) {
439 
440  // Set energy indices for exposure computation
443 
444  // Set energy indices
445  int inx_energy_left = m_energy.inx_left() * noffsets();
446  int inx_energy_right = m_energy.inx_right() * noffsets();
447 
448  // Set array indices for bi-linear interpolation
449  m_inx1 = m_offset.inx_left() + inx_energy_left;
450  m_inx2 = m_offset.inx_left() + inx_energy_right;
451  m_inx3 = m_offset.inx_right() + inx_energy_left;
452  m_inx4 = m_offset.inx_right() + inx_energy_right;
453 
454  // Set weighting factors for bi-linear interpolation
459 
460  //} // endif: logE or ctheta changed
461 
462  // Compute energy dependentmap corrections
463  double fac_left = m_mapcorr[m_inx1_exp];
464  double fac_right = m_mapcorr[m_inx2_exp];
465 
466  // Perform bi-linear interpolation
467  value = m_wgt1 * m_psf[m_inx1] * fac_left +
468  m_wgt2 * m_psf[m_inx2] * fac_right +
469  m_wgt3 * m_psf[m_inx3] * fac_left +
470  m_wgt4 * m_psf[m_inx4] * fac_right;
471 
472  // Optionally check for negative values
473  #if defined(G_SIGNAL_NEGATIVE_MEAN_PSF)
474  if (value < 0.0) {
475  std::cout << "GLATMeanPsf::psf(offset=" << offset;
476  std::cout << ", logE=" << logE << ") = " << value << std::endl;
477  }
478  #endif
479 
480  } // endif: arguments were valid
481 
482  // Return value
483  return value;
484 }
485 
486 
487 /***********************************************************************//**
488  * @brief Return exposure value
489  *
490  * @param[in] logE log10 of energy in MeV.
491  *
492  * Computes
493  * \f[{\rm exposure}(\log E)\f]
494  * by linear interpolation, where
495  * \f$\log E\f$ is the logarithm of base 10 of the energy in MeV.
496  * A zero value is returned if \f$\log E\f$ is not positive.
497  ***************************************************************************/
498 double GLATMeanPsf::exposure(const double& logE)
499 {
500  // Initialise response
501  double value = 0.0;
502 
503  // Continue only if arguments are within valid range
504  if (logE > 0.0) {
505 
506  // Set energy interpolation
507  //if (logE != m_last_energy) {
508  m_energy.set_value(logE);
509  // m_last_energy = logE;
512  //}
513 
514  // Perform linear interpolation
515  value = m_energy.wgt_left() * m_exposure[m_inx1_exp] +
517 
518  } // endif: arguments were in valid range
519 
520  // Return value
521  return value;
522 }
523 
524 
525 /***********************************************************************//**
526  * @brief Print livetime cube information
527  *
528  * @param[in] chatter Chattiness (defaults to NORMAL).
529  * @return String containing ROI information.
530  ***************************************************************************/
531 std::string GLATMeanPsf::print(const GChatter& chatter) const
532 {
533  // Initialise result string
534  std::string result;
535 
536  // Continue only if chatter is not silent
537  if (chatter != SILENT) {
538 
539  // Compute exposure range
540  double min_exposure = 0.0;
541  double max_exposure = 0.0;
542  for (int i = 0; i < nenergies(); ++i) {
543  if (i == 0) {
544  min_exposure = m_exposure[i];
545  max_exposure = m_exposure[i];
546  }
547  else {
548  if (m_exposure[i] < min_exposure) {
549  min_exposure = m_exposure[i];
550  }
551  if (m_exposure[i] > max_exposure) {
552  max_exposure = m_exposure[i];
553  }
554  }
555  }
556 
557  // Append header
558  result.append("=== GLATMeanPsf ===");
559 
560  // Append information
561  result.append("\n"+gammalib::parformat("Source name")+name());
562  result.append("\n"+gammalib::parformat("Source direction"));
563  result.append(gammalib::str(m_dir.ra_deg()));
564  result.append(", ");
565  result.append(gammalib::str(m_dir.dec_deg()));
566  result.append("\n"+gammalib::parformat("Offset angles")+gammalib::str(noffsets()));
567  result.append("\n"+gammalib::parformat("Energy values")+gammalib::str(nenergies()));
568  result.append("\n"+gammalib::parformat("Exposure range"));
569  result.append(gammalib::str(min_exposure)+" - "+gammalib::str(max_exposure)+" s cm2");
570  for (int i = 0; i < nenergies(); ++i) {
571  GEnergy energy;
572  energy.log10MeV(m_energy[i]);
573  result.append("\n"+gammalib::parformat(energy.print()));
574  result.append(gammalib::str(m_exposure[i])+" s cm2");
575  if (m_mapcorr.size() == nenergies()) {
576  result.append(" (map correction="+gammalib::str(m_mapcorr[i])+")");
577  }
578  }
579 
580  } // endif: chatter was not silent
581 
582  // Return result
583  return result;
584 }
585 
586 
587 /*==========================================================================
588  = =
589  = Private methods =
590  = =
591  ==========================================================================*/
592 
593 /***********************************************************************//**
594  * @brief Initialise class members
595  ***************************************************************************/
597 {
598  // Initialise members
599  m_name.clear();
600  m_dir.clear();
601  m_psf.clear();
602  m_exposure.clear();
603  m_mapcorr.clear();
604  m_energy.clear();
605  m_offset.clear();
606  m_theta_max = 70.0; //!< Maximum zenith angle
607  m_last_energy = -1.0;
608  m_last_offset = -1.0;
609  m_inx1_exp = 0;
610  m_inx2_exp = 0;
611  m_inx1 = 0;
612  m_inx2 = 0;
613  m_inx3 = 0;
614  m_inx4 = 0;
615  m_wgt1 = 0.0;
616  m_wgt2 = 0.0;
617  m_wgt3 = 0.0;
618  m_wgt4 = 0.0;
619 
620  // Set offset array
621  set_offsets();
622 
623  // Return
624  return;
625 }
626 
627 
628 /***********************************************************************//**
629  * @brief Copy class members
630  *
631  * @param[in] psf Mean PSF.
632  ***************************************************************************/
634 {
635  // Copy members
636  m_name = psf.m_name;
637  m_dir = psf.m_dir;
638  m_psf = psf.m_psf;
639  m_exposure = psf.m_exposure;
640  m_mapcorr = psf.m_mapcorr;
641  m_energy = psf.m_energy;
642  m_offset = psf.m_offset;
643  m_theta_max = psf.m_theta_max;
646  m_inx1_exp = psf.m_inx1_exp;
647  m_inx2_exp = psf.m_inx2_exp;
648  m_inx1 = psf.m_inx1;
649  m_inx2 = psf.m_inx2;
650  m_inx3 = psf.m_inx3;
651  m_inx4 = psf.m_inx4;
652  m_wgt1 = psf.m_wgt1;
653  m_wgt2 = psf.m_wgt2;
654  m_wgt3 = psf.m_wgt3;
655  m_wgt4 = psf.m_wgt4;
656 
657  // Return
658  return;
659 }
660 
661 
662 /***********************************************************************//**
663  * @brief Delete class members
664  ***************************************************************************/
666 {
667  // Return
668  return;
669 }
670 
671 
672 /***********************************************************************//**
673  * @brief Set array of offset values in degrees
674  *
675  * The array of offset values defines the points at with the mean PSF will
676  * be evaluated and stored. The array is logarithmically spaced.
677  ***************************************************************************/
679 {
680  // Set array parameters
681  const double offset_min = 1.0e-4;
682  const double offset_max = 70.0;
683  const int offset_num = 200;
684 
685  // Clear offset array
686  m_offset.clear();
687 
688  // Set array
689  double step = log(offset_max/offset_min)/(offset_num - 1.0);
690  for (int i = 0; i < offset_num; ++i)
691  m_offset.append(offset_min*exp(i*step));
692 
693  // Return
694  return;
695 }
696 
697 
698 /***********************************************************************//**
699  * @brief Compute map corrections
700  *
701  * @param[in] obs LAT observation.
702  *
703  * Sets up a vector of energy dependent corrections that assures that the
704  * integral over the PSF in a specific event cube is correctly normalised.
705  * For this purpose the PSF is integrated over a radius that fully lies
706  * within the event cube and devided by the PSF pixel sum within this radius.
707  * The map corrections are mainly important at high energies where the
708  * PSF size can become comparable to the pixel size.
709  *
710  * See SourceMap::getMapCorrections for the original code in the Fermi/LAT
711  * ScienceTools.
712  *
713  * @todo We can also implement a method for event atoms, yet it is not clear
714  * whether we really need this.
715  ***************************************************************************/
717 {
718  // Initialise map corrections to unity vector
719  m_mapcorr.assign(m_energy.size(), 1.0);
720 
721  // Get pointer on event cube
722  const GLATEventCube* cube = dynamic_cast<const GLATEventCube*>(obs.events());
723 
724  // Continue only if we have an event cube
725  if (cube != NULL) {
726 
727  // Get maximum PSF radius
728  double radius = cube->maxrad(m_dir);
729  #if defined(G_DEBUG_MAP_CORRECTIONS)
730  std::cout << "GLATMeanPsf::set_map_corrections:";
731  std::cout << " radius=" << radius << std::endl;
732  #endif
733  if (radius > 0.0) {
734 
735  // Initialise sum over pixels
736  std::vector<double> sum(m_energy.size(), 0.0);
737 
738  // Loop over all event cube spatial pixels
739  for (int iy = 0; iy < cube->ny(); ++iy) {
740  for (int ix = 0; ix < cube->nx(); ++ix) {
741 
742  // Compute offset angle in degrees
743  GSkyPixel pixel = GSkyPixel(double(ix), double(iy));
744  double offset = cube->map().pix2dir(pixel).dist_deg(m_dir);
745  double solidangle = cube->map().solidangle(pixel);
746 
747  // Use only pixels within maximum PSF radius
748  if (offset <= radius) {
749 
750  // Accumulate energy dependent pixel sum
751  for (int ieng = 0; ieng < m_energy.size(); ++ieng) {
752  sum[ieng] += psf(offset, m_energy[ieng]) * solidangle;
753  }
754 
755  } // endif: pixel was within maximum PSF radius
756 
757  } // endfor: looped over latitude pixels
758  } // endfor: looped over longitude pixels
759 
760  // Compute map correction
761  for (int ieng = 0; ieng < m_energy.size(); ++ieng) {
762  if (sum[ieng] > 0.0) {
763  m_mapcorr[ieng] = integral(radius, m_energy[ieng]) / sum[ieng];
764  }
765  }
766 
767  } // endif: radius was positive
768  } // endif: observation had an event cube
769 
770  // Return
771  return;
772 }
773 
774 
775 /***********************************************************************//**
776  * @brief Compute integral over PSF
777  *
778  * @param[in] offsetmax Maximum offset angle.
779  * @param[in] logE log10 of energy in MeV.
780  ***************************************************************************/
781 double GLATMeanPsf::integral(const double& offsetmax, const double& logE)
782 {
783  // Set energy and offset interpolation weigths
784  m_energy.set_value(logE);
785  m_offset.set_value(offsetmax);
786 
787  // Get PSF array offsets
788  int inx_energy_left = m_energy.inx_left() * noffsets();
789  int inx_energy_right = m_energy.inx_right() * noffsets();
790 
791  // Initialise integrals
792  double int_left = 0.0;
793  double int_right = 0.0;
794 
795  // Loop over all offset angles
796  for (int i = 0; i < m_offset.size()-1; ++i) {
797 
798  // If interval is fully contained within offset then simply use
799  // tabulated values for integration
800  if (m_offset[i+1] <= offsetmax) {
801  double theta_min = m_offset[i] * gammalib::deg2rad;
802  double theta_max = m_offset[i+1] * gammalib::deg2rad;
803  int_left += 0.5 * (m_psf[inx_energy_left+i] * sin(theta_min) +
804  m_psf[inx_energy_left+i+1] * sin(theta_max)) *
805  (theta_max - theta_min);
806  int_right += 0.5 * (m_psf[inx_energy_right+i] * sin(theta_min) +
807  m_psf[inx_energy_right+i+1] * sin(theta_max)) *
808  (theta_max - theta_min);
809  }
810 
811  // ... otherwise interpolate the PSF value for theta_max. We can
812  // then exit the loop since we're done with the integration
813  else {
814  double theta_min = m_offset[i] * gammalib::deg2rad;
815  double theta_max = offsetmax * gammalib::deg2rad;
816  double psf_left = m_psf[inx_energy_left+i] * m_offset.wgt_left() +
817  m_psf[inx_energy_left+i+1] * m_offset.wgt_right();
818  double psf_right = m_psf[inx_energy_right+i] * m_offset.wgt_left() +
819  m_psf[inx_energy_right+i+1] * m_offset.wgt_right();
820  int_left += 0.5 * (m_psf[inx_energy_left+i] * sin(theta_min) +
821  psf_left * sin(theta_max)) *
822  (theta_max - theta_min);
823  int_right += 0.5 * (m_psf[inx_energy_right+i] * sin(theta_min) +
824  psf_right * sin(theta_max)) *
825  (theta_max - theta_min);
826 
827  // Debug option: Dump integral computation results
828  #if defined(G_DEBUG_INTEGRAL)
829  std::cout << "offsetmax=" << offsetmax;
830  std::cout << " offset.left=" << m_offset.inx_left();
831  std::cout << " offset.right=" << m_offset.inx_right();
832  std::cout << " i=" << i;
833  std::cout << " psf_left=" << psf_left;
834  std::cout << " psf_right=" << psf_right;
835  std::cout << " int_left=" << int_left*gammalib::twopi;
836  std::cout << " int_right=" << int_right*gammalib::twopi;
837  #endif
838  break;
839  }
840 
841  } // endfor: looped over offset angles
842 
843  // Interpolate now in energy
844  double integral = gammalib::twopi * (m_energy.wgt_left() * int_left +
845  m_energy.wgt_right() * int_right);
846 
847  // Debug option: Dump integral computation results
848  #if defined(G_DEBUG_INTEGRAL)
849  std::cout << " integral=" << integral << std::endl;
850  #endif
851 
852  // Return integral
853  return integral;
854 }
int m_inx3
Index 3.
Fermi/LAT observation class.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:192
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:256
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:163
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:264
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:944
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:587
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:279
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
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
double m_last_energy
Last requested logE value.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:229
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:235
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:249
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:164
int nenergies(void) const
Return number of energy bins.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
virtual GEvents * events(void)
Return events.
Exception handler interface definition.
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:1232
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:199
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:1143
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:489