GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GLATResponse.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GLATResponse.cpp - Fermi LAT response class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2008-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 GLATResponse.cpp
23  * @brief Fermi LAT response class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <unistd.h> // access() function
32 #include <typeinfo>
33 #include <cstdlib> // std::getenv() function
34 #include <string>
35 #include "GException.hpp"
36 #include "GFits.hpp"
37 #include "GTools.hpp"
38 #include "GCaldb.hpp"
39 #include "GSource.hpp"
41 #include "GLATInstDir.hpp"
42 #include "GLATResponse.hpp"
43 #include "GLATObservation.hpp"
44 #include "GLATEventAtom.hpp"
45 #include "GLATEventBin.hpp"
46 #include "GLATEventCube.hpp"
47 
48 /* __ Method name definitions ____________________________________________ */
49 #define G_CALDB "GLATResponse::caldb(std::string&)"
50 #define G_LOAD "GLATResponse::load(std::string&)"
51 #define G_IRF "GLATResponse::irf(GInstDir&, GEnergy&, GTime&, GSkyDir&,"\
52  " GEnergy&, GTime&, GObservation&)"
53 #define G_NROI "GLATResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
54  "GObservation&)"
55 #define G_EBOUNDS "GLATResponse::ebounds(GEnergy&)"
56 #define G_SAVE "GLATResponse::save(std::string&)"
57 #define G_AEFF "GLATResponse::aeff(int&)"
58 #define G_PSF "GLATResponse::psf(int&)"
59 #define G_EDISP "GLATResponse::edisp(int&)"
60 #define G_IRF_ATOM "GLATResponse::irf(GLATEventAtom&, GModel&, GEnergy&,"\
61  "GTime&, GObservation&)"
62 #define G_IRF_BIN "GLATResponse::irf(GLATEventBin&, GModel&, GEnergy&,"\
63  "GTime&, GObservation&)"
64 
65 /* __ Macros _____________________________________________________________ */
66 
67 /* __ Coding definitions _________________________________________________ */
68 
69 /* __ Debug definitions __________________________________________________ */
70 //#define G_DUMP_MEAN_PSF //!< Dump mean PSF allocation
71 //#define G_DEBUG_MEAN_PSF //!< Debug mean PSF computation
72 
73 /* __ Constants __________________________________________________________ */
74 
75 
76 /*==========================================================================
77  = =
78  = Constructors/destructors =
79  = =
80  ==========================================================================*/
81 
82 /***********************************************************************//**
83  * @brief Void constructor
84  ***************************************************************************/
86 {
87  // Initialise members
88  init_members();
89 
90  // Return
91  return;
92 }
93 
94 
95 /***********************************************************************//**
96  * @brief Copy constructor
97  *
98  * @param[in] rsp Response.
99  ***************************************************************************/
101 {
102  // Initialise members
103  init_members();
104 
105  // Copy members
106  copy_members(rsp);
107 
108  // Return
109  return;
110 }
111 
112 
113 /***********************************************************************//**
114  * @brief Destructor
115  ***************************************************************************/
117 {
118  // Free members
119  free_members();
120 
121  // Return
122  return;
123 }
124 
125 
126 /*==========================================================================
127  = =
128  = Operators =
129  = =
130  ==========================================================================*/
131 
132 /***********************************************************************//**
133  * @brief Assignment operator
134  *
135  * @param[in] rsp Response.
136  * @return Response.
137  ***************************************************************************/
139 {
140  // Execute only if object is not identical
141  if (this != &rsp) {
142 
143  // Copy base class members
144  this->GResponse::operator=(rsp);
145 
146  // Free members
147  free_members();
148 
149  // Initialise members
150  init_members();
151 
152  // Copy members
153  copy_members(rsp);
154 
155  } // endif: object was not identical
156 
157  // Return this object
158  return *this;
159 }
160 
161 
162 /*==========================================================================
163  = =
164  = Public methods =
165  = =
166  ==========================================================================*/
167 
168 /***********************************************************************//**
169  * @brief Clear response
170 ***************************************************************************/
172 {
173  // Free class members
174  free_members();
175  this->GResponse::free_members();
176 
177  // Initialise members
178  this->GResponse::init_members();
179  init_members();
180 
181  // Return
182  return;
183 }
184 
185 
186 /***********************************************************************//**
187  * @brief Clone response
188  *
189  * @return Pointer to deep copy of response.
190  ***************************************************************************/
192 {
193  return new GLATResponse(*this);
194 }
195 
196 
197 /***********************************************************************//**
198  * @brief Return value of point source IRF
199  *
200  * @param[in] event Observed event.
201  * @param[in] photon Incident photon.
202  * @param[in] obs Observation.
203  *
204  * @exception GException::invalid_argument
205  * Instrument direction is not a valid LAT instrument direction.
206  *
207  * @todo The IRF value is not divided by ontime of the event, but it is
208  * already time integrated.
209  ***************************************************************************/
210 double GLATResponse::irf(const GEvent& event,
211  const GPhoton& photon,
212  const GObservation& obs) const
213 {
214  // Get pointer to LAT instrument direction
215  const GLATInstDir* dir = dynamic_cast<const GLATInstDir*>(&(event.dir()));
216 
217  // If pointer is not valid then throw an exception
218  if (dir == NULL) {
219  std::string cls = std::string(typeid(&(event.dir())).name());
220  std::string msg = "Invalid instrument direction type \""+cls+
221  "\" specified. Please specify a \"GLATInstDir\" "
222  "instance as argument.";
224  }
225 
226  // Get photon attributes
227  const GSkyDir& srcDir = photon.dir();
228  const GEnergy& srcEng = photon.energy();
229 
230  // Search for mean PSF
231  int ipsf = -1;
232  for (int i = 0; i < m_ptsrc.size(); ++i) {
233  if (m_ptsrc[i]->dir() == srcDir) {
234  ipsf = i;
235  break;
236  }
237  }
238 
239  // If mean PSF has not been found then create it now
240  if (ipsf == -1) {
241 
242  // Allocate new mean PSF
243  GLATMeanPsf* psf = new GLATMeanPsf(srcDir, static_cast<const GLATObservation&>(obs));
244 
245  // Set source name
246  std::string name = "SRC("+gammalib::str(srcDir.ra_deg()) +
247  "," +
248  gammalib::str(srcDir.dec_deg())+")";
249  psf->name(name);
250 
251  // Push mean PSF on stack
252  const_cast<GLATResponse*>(this)->m_ptsrc.push_back(psf);
253 
254  // Set index of mean PSF
255  ipsf = m_ptsrc.size()-1;
256 
257  // Debug option: dump mean PSF
258  #if defined(G_DUMP_MEAN_PSF)
259  std::cout << "Added new mean PSF \""+name+"\"" << std::endl;
260  std::cout << *psf << std::endl;
261  #endif
262 
263  } // endif: created new mean PSF
264 
265  // Get IRF value
266  double offset = dir->dir().dist_deg(srcDir);
267  double irf = (*m_ptsrc[ipsf])(offset, srcEng.log10MeV());
268 
269  // Return IRF value
270  return irf;
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Return value of model instrument response function
276  *
277  * @param[in] event Event.
278  * @param[in] source Source.
279  * @param[in] obs Observation.
280  *
281  * This method returns the response of the instrument to a specific source
282  * model. The method handles both event atoms and event bins.
283  ***************************************************************************/
284 double GLATResponse::irf_spatial(const GEvent& event,
285  const GSource& source,
286  const GObservation& obs) const
287 {
288  // Get IRF value
289  double rsp;
290  if (event.is_atom()) {
291  rsp = irf_spatial_atom(static_cast<const GLATEventAtom&>(event), source, obs);
292  }
293  else {
294  rsp = irf_spatial_bin(static_cast<const GLATEventBin&>(event), source, obs);
295  }
296 
297  // Return IRF value
298  return rsp;
299 }
300 
301 
302 /***********************************************************************//**
303  * @brief Return value of model IRF for event atom
304  *
305  * @param[in] event Event atom.
306  * @param[in] source Source.
307  * @param[in] obs Observations.
308  *
309  * @exception GException::feature_not_implemented
310  * Method not yet implemented.
311  *
312  * @todo Not yet implemented.
313  ***************************************************************************/
315  const GSource& source,
316  const GObservation& obs) const
317 {
318  // Initialise IRF with "no response"
319  double irf = 0.0;
320 
321  // Dump warning that integration is not yet implemented
323 
324  // Return IRF value
325  return irf;
326 }
327 
328 
329 /***********************************************************************//**
330  * @brief Return value of model IRF for event bin
331  *
332  * @param[in] event Event bin.
333  * @param[in] source Source.
334  * @param[in] obs Observation.
335  *
336  * @exception GException::invalid_value
337  * Diffuse model not found.
338  *
339  * This method first searches for a corresponding source map, and if found,
340  * computes the response from the source map. If no source map is present
341  * it checks if the source is a point source. If this is the case a mean
342  * PSF is allocated for the source and the response is computed from the
343  * mean PSF. Otherwise an GException::invalid_value exception is thrown.
344  *
345  * @todo Extract event cube from observation. We do not need the cube
346  * pointer in the event anymore.
347  * @todo Implement more efficient method for response search (name search is
348  * not very rapid).
349  * @todo Instead of calling "offset = event.dir().dist_deg(srcDir)" we can
350  * precompute and store for each PSF the offsets. This should save
351  * quite some time since the distance computation is time
352  * consuming.
353  ***************************************************************************/
355  const GSource& source,
356  const GObservation& obs) const
357 {
358  // Initialise response value
359  double rsp = 0.0;
360 
361  // Get pointer to event cube
362  const GLATEventCube* cube = event.cube();
363 
364  // Get pointer on point source spatial model
365  const GModelSpatialPointSource* ptsrc =
366  dynamic_cast<const GModelSpatialPointSource*>(source.model());
367 
368  // Get source energy
369  GEnergy srcEng = source.energy();
370 
371  // Search for diffuse response in event cube
372  int idiff = -1;
373  for (int i = 0; i < cube->ndiffrsp(); ++i) {
374  if (cube->diffname(i) == source.name()) {
375  idiff = i;
376  break;
377  }
378  }
379 
380  // If diffuse response has been found then get response from source map
381  if (idiff != -1) {
382 
383  // Get srcmap indices and weighting factors
384  GNodeArray nodes = cube->enodes();
385  nodes.set_value(srcEng.log10MeV());
386 
387  // Compute diffuse response
388  GSkyMap* map = cube->diffrsp(idiff);
389  const double* pixels = map->pixels() + event.ipix();
390  rsp = nodes.wgt_left() * pixels[nodes.inx_left() * map->npix()] +
391  nodes.wgt_right() * pixels[nodes.inx_right() * map->npix()];
392 
393  // Divide by solid angle and ontime since source maps are given in units of
394  // counts/pixel/MeV.
395  rsp /= (event.solidangle() * event.ontime());
396 
397  } // endelse: diffuse response was present
398 
399  // ... otherwise check if model is a point source. If this is true
400  // then return response from mean PSF
401  if ((idiff == -1 || m_force_mean) && ptsrc != NULL) {
402 
403  // Search for mean PSF. If a mean PSF is found then check if the PSF
404  // needs to be updated, which may occur if the source direction
405  // changed by more then 0.01 deg.
406  int ipsf = -1;
407  bool update = false;
408  for (int i = 0; i < m_ptsrc.size(); ++i) {
409  if (m_ptsrc[i]->name() == source.name()) {
410  ipsf = i;
411  if (m_ptsrc[i]->dir().dist_deg(ptsrc->dir()) > 0.01) {
412  update = true;
413  }
414  break;
415  }
416  }
417 
418  // If mean PSF has not been found then create it now
419  if (ipsf == -1) {
420 
421  // Allocate new mean PSF
422  GLATMeanPsf* psf = new GLATMeanPsf(ptsrc->dir(),
423  static_cast<const GLATObservation&>(obs));
424 
425  // Set source name
426  psf->name(source.name());
427 
428  // Push mean PSF on stack
429  const_cast<GLATResponse*>(this)->m_ptsrc.push_back(psf);
430 
431  // Set index of mean PSF
432  ipsf = m_ptsrc.size()-1;
433 
434  // Debug option: dump mean PSF
435  #if defined(G_DUMP_MEAN_PSF)
436  std::cout << "Added new mean PSF \""+source.name();
437  std::cout << "\"" << std::endl << *psf << std::endl;
438  #endif
439 
440  } // endif: created new mean PSF
441 
442  // ... otherwise if an update is needed then update mean PSF now
443  else if (update) {
444 
445  // Allocate new mean PSF
446  GLATMeanPsf* psf = new GLATMeanPsf(ptsrc->dir(),
447  static_cast<const GLATObservation&>(obs));
448 
449  // Set source name
450  psf->name(source.name());
451 
452  // Delete existing mean PSF
453  if (m_ptsrc[ipsf] != NULL) delete m_ptsrc[ipsf];
454 
455  // Set new mean PSF
456  const_cast<GLATResponse*>(this)->m_ptsrc[ipsf] = psf;
457 
458  // Debug option: dump mean PSF
459  #if defined(G_DUMP_MEAN_PSF)
460  std::cout << "Replaced mean PSF \""+source.name();
461  std::cout << "\"" << std::endl << *psf << std::endl;
462  #endif
463 
464  } // endelse: Mean PSF update was needed
465 
466  // Get PSF value
467  GSkyDir srcDir = ptsrc->dir();
468  double offset = event.dir().dir().dist_deg(srcDir);
469  double mean_psf = (*m_ptsrc[ipsf])(offset, srcEng.log10MeV()) / (event.ontime());
470 
471  // Debug option: compare mean PSF to diffuse response
472  #if defined(G_DEBUG_MEAN_PSF)
473  std::cout << "Energy=" << srcEng.MeV();
474  std::cout << " MeanPsf=" << mean_psf;
475  std::cout << " DiffusePsf=" << rsp;
476  std::cout << " ratio(Mean/Diffuse)=" << mean_psf/rsp;
477  if (mean_psf/rsp < 0.99) {
478  std::cout << " <<< (1%)";
479  }
480  else if (mean_psf/rsp > 1.01) {
481  std::cout << " >>> (1%)";
482  }
483  std::cout << std::endl;
484  #endif
485 
486  // Set response
487  rsp = mean_psf;
488 
489  } // endif: model was point source
490 
491  // ... otherwise throw an exception
492  if ((idiff == -1) && ptsrc == NULL) {
493  std::string msg = "No diffuse source with name \""+source.name()+
494  "\" found and source is not a point source. So far "
495  "the LAT module can only handle diffuse sources and "
496  "point sources. Please adapt your model accordingly.";
498  }
499 
500  // Return IRF value
501  return rsp;
502 }
503 
504 
505 /***********************************************************************//**
506  * @brief Return integral of event probability for a given sky model over ROI
507  *
508  * @param[in] model Sky model.
509  * @param[in] obsEng Observed photon energy.
510  * @param[in] obsTime Observed photon arrival time.
511  * @param[in] obs Observation.
512  * @return 0
513  *
514  * @exception GException::feature_not_implemented
515  * Method is not implemented.
516  ***************************************************************************/
517 double GLATResponse::nroi(const GModelSky& model,
518  const GEnergy& obsEng,
519  const GTime& obsTime,
520  const GObservation& obs) const
521 {
522  // Method is not implemented
523  std::string msg = "Spatial integration of sky model over the data space "
524  "is not implemented.";
526 
527  // Return Npred
528  return (0.0);
529 }
530 
531 
532 /***********************************************************************//**
533  * @brief Return true energy boundaries for a specific observed energy
534  *
535  * @param[in] obsEnergy Observed Energy.
536  * @return True energy boundaries for given observed energy.
537  *
538  * @exception GException::feature_not_implemented
539  * Method is not implemented.
540  ***************************************************************************/
541 GEbounds GLATResponse::ebounds(const GEnergy& obsEnergy) const
542 {
543  // Initialise an empty boundary object
545 
546  // Throw an exception
547  std::string msg = "Energy dispersion not implemented.";
549 
550  // Return energy boundaries
551  return (ebounds);
552 }
553 
554 
555 /***********************************************************************//**
556  * @brief Load Fermi LAT response from calibration database
557  *
558  * @param[in] rspname Response name.
559  *
560  * @exception GException::invalid_argument
561  * Invalid response type encountered.
562  *
563  * Loads the specified Fermi LAT response from the calibration database.
564  * The following response names are supported (case insensitive):
565  *
566  * name (is equivalent to front+back)
567  * name::front
568  * name::back
569  * name::psf(0-3)
570  * name::edisp(0-3)
571  *
572  * where name is the response name (for example "P8R2_SOURCE_V6"). Note that
573  * the name is case sensitive, but the event typ is not case sensitive.
574  ***************************************************************************/
575 void GLATResponse::load(const std::string& rspname)
576 {
577  // Set possible response types
578  static const std::string type_names[] = {"FRONT",
579  "BACK",
580  "PSF0",
581  "PSF1",
582  "PSF2",
583  "PSF3",
584  "EDISP0",
585  "EDISP1",
586  "EDISP2",
587  "EDISP3"};
588 
589  // Clear instance
590  clear();
591 
592  // Allocate Fermi LAT calibration database
593  GCaldb caldb("glast", "lat");
594 
595  // Determine response types to be loaded. If no event type has been
596  // specified then load both front and back response.
597  int event_type = 0;
598  std::vector<std::string> array = gammalib::split(rspname, ":");
599  if (array.size() == 1) {
600  m_rspname = array[0];
601  event_type = 3; // 0x0000000011
602  }
603  else if (array.size() == 3) {
604  m_rspname = array[0];
605  std::string evtype = gammalib::strip_whitespace(gammalib::tolower(array[2]));
606  if (evtype == "front") {
607  event_type = 1; // 0x0000000001
608  }
609  else if (evtype == "back") {
610  event_type = 2; // 0x0000000010
611  }
612  else if (evtype == "psf0") {
613  event_type = 4; // 0x0000000100
614  }
615  else if (evtype == "psf1") {
616  event_type = 8; // 0x0000001000
617  }
618  else if (evtype == "psf2") {
619  event_type = 16; // 0x0000010000
620  }
621  else if (evtype == "psf3") {
622  event_type = 32; // 0x0000100000
623  }
624  else if (evtype == "psf") {
625  event_type = 60; // 0x0000111100
626  }
627  else if (evtype == "edisp0") {
628  event_type = 64; // 0x0001000000
629  }
630  else if (evtype == "edisp1") {
631  event_type = 128; // 0x0010000000
632  }
633  else if (evtype == "edisp2") {
634  event_type = 256; // 0x0100000000
635  }
636  else if (evtype == "edisp3") {
637  event_type = 512; // 0x1000000000
638  }
639  else if (evtype == "edisp") {
640  event_type = 960; // 0x1111000000
641  }
642  else {
643  std::string msg = "Invalid response type \""+array[2]+"\". "
644  "Expect one of \"FRONT\", \"BACK\", "
645  "\"PSF(0-3)\", or \"EDISP(0-3)\" "
646  "(case insensitive).";
648  }
649  }
650  else {
651  std::string msg = "Invalid response \""+m_rspname+"\".";
653  }
654 
655  // Loop over all possible response types and append the relevant
656  // response function components to the response
657  for (int i = 0; i < 10; ++i) {
658 
659  // Set bitmask for this response type
660  int bitmask = 1 << i;
661 
662  // Fall through if this response type is not requested
663  if ((event_type & bitmask) == 0) {
664  continue;
665  }
666 
667  // Get response using the GCaldb interface
668  std::string expr = "VERSION("+gammalib::toupper(m_rspname)+")";
669  GFilename aeffname = caldb.filename(type_names[i],"","EFF_AREA","","",expr);
670  GFilename psfname = caldb.filename(type_names[i],"","RPSF","","",expr);
671  GFilename edispname = caldb.filename(type_names[i],"","EDISP","","",expr);
672 
673  // Load IRF components
674  GLATAeff* aeff = new GLATAeff(aeffname, type_names[i]);
675  GLATPsf* psf = new GLATPsf(psfname, type_names[i]);
676  GLATEdisp* edisp = new GLATEdisp(edispname, type_names[i]);
677 
678  // Push IRF components on the response stack
679  m_aeff.push_back(aeff);
680  m_psf.push_back(psf);
681  m_edisp.push_back(edisp);
682 
683  } // endfor: looped over response types
684 
685  // Return
686  return;
687 }
688 
689 
690 /***********************************************************************//**
691  * @brief Save Fermi LAT response in calibration database
692  *
693  * @param[in] rspname Response name.
694  *
695  * @todo Not yet implemented.
696  ***************************************************************************/
697 void GLATResponse::save(const std::string& rspname) const
698 {
699  // Throw an exception
700  std::string msg = "Saving of LAT response not implemented.";
702 
703  // Return
704  return;
705 }
706 
707 
708 /***********************************************************************//**
709  * @brief Return pointer on effective area
710  *
711  * @param[in] index Response index [0,...,m_aeff.size()-1].
712  ***************************************************************************/
713 GLATAeff* GLATResponse::aeff(const int& index) const
714 {
715  // Optionally check if the index is valid
716  #if defined(G_RANGE_CHECK)
717  if (index < 0 || index >= m_aeff.size()) {
718  throw GException::out_of_range(G_AEFF, "Response index", index,
719  m_aeff.size());
720  }
721  #endif
722 
723  // Return pointer on effective area
724  return m_aeff[index];
725 }
726 
727 
728 /***********************************************************************//**
729  * @brief Return pointer on point spread function
730  *
731  * @param[in] index Response index [0,...,m_psf.size()-1].
732  ***************************************************************************/
733 GLATPsf* GLATResponse::psf(const int& index) const
734 {
735  // Optionally check if the index is valid
736  #if defined(G_RANGE_CHECK)
737  if (index < 0 || index >= m_psf.size()) {
738  throw GException::out_of_range(G_PSF, "Response index", index,
739  m_psf.size());
740  }
741  #endif
742 
743  // Return pointer on point spread function
744  return m_psf[index];
745 }
746 
747 
748 /***********************************************************************//**
749  * @brief Return pointer on energy dispersion
750  *
751  * @param[in] index Response index [0,...,m_edisp.size()-1].
752  ***************************************************************************/
753 GLATEdisp* GLATResponse::edisp(const int& index) const
754 {
755  // Optionally check if the index is valid
756  #if defined(G_RANGE_CHECK)
757  if (index < 0 || index >= m_edisp.size()) {
758  throw GException::out_of_range(G_EDISP, "Response index", index,
759  m_edisp.size());
760  }
761  #endif
762 
763  // Return pointer on energy dispersion
764  return m_edisp[index];
765 }
766 
767 
768 /***********************************************************************//**
769  * @brief Print Fermi-LAT response information
770  *
771  * @param[in] chatter Chattiness.
772  * @return String containing response information.
773  ***************************************************************************/
774 std::string GLATResponse::print(const GChatter& chatter) const
775 {
776  // Initialise result string
777  std::string result;
778 
779  // Continue only if chatter is not silent
780  if (chatter != SILENT) {
781 
782  // Append header
783  result.append("=== GLATResponse ===");
784 
785  // Append information
786  result.append("\n"+gammalib::parformat("Response name")+m_rspname);
787  result.append("\n"+gammalib::parformat("Event types"));
788  for (int i = 0; i < size(); ++i) {
789  result.append(m_aeff[i]->evtype()+" ");
790  }
791 
792  if (chatter > TERSE) {
793  for (int i = 0; i < size(); ++i) {
794  result.append("\n"+m_aeff[i]->print(gammalib::reduce(chatter)));
795  result.append("\n"+m_psf[i]->print(gammalib::reduce(chatter)));
796  result.append("\n"+m_edisp[i]->print(gammalib::reduce(chatter)));
797  }
798  for (int i = 0; i < m_ptsrc.size(); ++i) {
799  result.append("\n"+m_ptsrc[i]->print(gammalib::reduce(chatter)));
800  }
801  }
802 
803  } // endif: chatter was not silent
804 
805  // Return result
806  return result;
807 }
808 
809 
810 /*==========================================================================
811  = =
812  = Private methods =
813  = =
814  ==========================================================================*/
815 
816 /***********************************************************************//**
817  * @brief Initialise class members
818  ***************************************************************************/
820 {
821  // Initialise members
822  m_rspname.clear();
823  m_force_mean = false;
824  m_aeff.clear();
825  m_psf.clear();
826  m_edisp.clear();
827  m_ptsrc.clear();
828 
829  // Return
830  return;
831 }
832 
833 
834 /***********************************************************************//**
835  * @brief Copy class members
836  *
837  * @param[in] rsp Response.
838  ***************************************************************************/
840 {
841  // Copy members
842  m_rspname = rsp.m_rspname;
844 
845  // Clone Aeff response components
846  m_aeff.clear();
847  for (int i = 0; i < rsp.m_aeff.size(); ++i) {
848  m_aeff.push_back(rsp.m_aeff[i]->clone());
849  }
850 
851  // Clone Psf response components
852  m_psf.clear();
853  for (int i = 0; i < rsp.m_psf.size(); ++i) {
854  m_psf.push_back(rsp.m_psf[i]->clone());
855  }
856 
857  // Clone Edisp response components
858  m_edisp.clear();
859  for (int i = 0; i < rsp.m_edisp.size(); ++i) {
860  m_edisp.push_back(rsp.m_edisp[i]->clone());
861  }
862 
863  // Clone point sources
864  m_ptsrc.clear();
865  for (int i = 0; i < rsp.m_ptsrc.size(); ++i) {
866  m_ptsrc.push_back(rsp.m_ptsrc[i]->clone());
867  }
868 
869  // Return
870  return;
871 }
872 
873 
874 /***********************************************************************//**
875  * @brief Delete class members
876  ***************************************************************************/
878 {
879  // Free Aeff memory
880  for (int i = 0; i < m_aeff.size(); ++i) {
881  if (m_aeff[i] != NULL) delete m_aeff[i];
882  m_aeff[i] = NULL;
883  }
884  m_aeff.clear();
885 
886  // Free Psf memory
887  for (int i = 0; i < m_psf.size(); ++i) {
888  if (m_psf[i] != NULL) delete m_psf[i];
889  m_psf[i] = NULL;
890  }
891  m_psf.clear();
892 
893  // Free Edisp memory
894  for (int i = 0; i < m_edisp.size(); ++i) {
895  if (m_edisp[i] != NULL) delete m_edisp[i];
896  m_edisp[i] = NULL;
897  }
898  m_edisp.clear();
899 
900  // Free point source memory
901  for (int i = 0; i < m_ptsrc.size(); ++i) {
902  if (m_ptsrc[i] != NULL) delete m_ptsrc[i];
903  m_ptsrc[i] = NULL;
904  }
905  m_ptsrc.clear();
906 
907  // Return
908  return;
909 }
910 
911 
912 /***********************************************************************//**
913  * @brief Return instrument response to point source
914  *
915  * @param[in] event Observed event.
916  * @param[in] source Source.
917  * @param[in] obs Observation.
918  * @return Instrument response to point source.
919  *
920  * Returns the instrument response to a point source.
921  ***************************************************************************/
922 double GLATResponse::irf_ptsrc(const GEvent& event,
923  const GSource& source,
924  const GObservation& obs) const
925 {
926  // Initialise IRF value
927  double irf = 0.0;
928 
929  // Get IRF value
930  if (event.is_atom()) {
931  irf = irf_spatial_atom(static_cast<const GLATEventAtom&>(event), source, obs);
932  }
933  else {
934  irf = irf_spatial_bin(static_cast<const GLATEventBin&>(event), source, obs);
935  }
936 
937  // Return IRF value
938  return irf;
939 }
940 
941 
942 /***********************************************************************//**
943  * @brief Return instrument response to diffuse source
944  *
945  * @param[in] event Observed event.
946  * @param[in] source Source.
947  * @param[in] obs Observation.
948  * @return Instrument response to point source.
949  *
950  * Returns the instrument response to a diffuse source.
951  ***************************************************************************/
952 double GLATResponse::irf_diffuse(const GEvent& event,
953  const GSource& source,
954  const GObservation& obs) const
955 {
956  // Initialise IRF value
957  double irf = 0.0;
958 
959  // Get IRF value
960  if (event.is_atom()) {
961  irf = irf_spatial_atom(static_cast<const GLATEventAtom&>(event), source, obs);
962  }
963  else {
964  irf = irf_spatial_bin(static_cast<const GLATEventBin&>(event), source, obs);
965  }
966 
967  // Return IRF value
968  return irf;
969 }
virtual ~GLATResponse(void)
Destructor.
#define G_EBOUNDS
Fermi/LAT observation class.
#define G_IRF
Sky map class.
Definition: GSkyMap.hpp:89
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:286
GSkyMap * diffrsp(const int &index) const
Return diffuse response map.
std::string m_rspname
Name of the instrument response.
void dir(const GSkyDir &dir)
Set sky direction.
Fermi/LAT observation class definition.
virtual GEbounds ebounds(const GEnergy &obsEnergy) const
Return true energy boundaries for a specific observed energy.
Node array class.
Definition: GNodeArray.hpp:60
int size(void) const
Return number of event types.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:256
Point source spatial model class interface definition.
const std::string & name(void) const
Return source name for mean PSF.
const GSkyDir & dir(void) const
Return position of point source.
virtual double irf_spatial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return value of model instrument response function.
int ndiffrsp(void) const
Return number of diffuse model components.
Fermi/LAT instrument direction class.
Definition: GLATInstDir.hpp:48
Interface for the Fermi LAT energy dispersion.
Definition: GLATEdisp.hpp:55
Abstract interface for the event classes.
Definition: GEvent.hpp:71
GLATPsf * psf(const int &index) const
Return pointer on point spread function.
GFilename filename(const std::string &detector, const std::string &filter, const std::string &codename, const std::string &date, const std::string &time, const std::string &expr)
Return calibration file name based on selection parameters.
Definition: GCaldb.cpp:524
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:264
Fermi/LAT Response class.
void load(const std::string &rspname)
Load Fermi LAT response from calibration database.
std::vector< GLATMeanPsf * > m_ptsrc
Mean PSFs for point sources.
double irf_spatial_bin(const GLATEventBin &event, const GSource &source, const GObservation &obs) const
Return value of model IRF for event bin.
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:983
GLATResponse & operator=(const GLATResponse &rsp)
Assignment operator.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
FITS file class interface definition.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:587
virtual double irf_diffuse(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to diffuse source.
const std::string & name(void) const
Return model name.
Definition: GSource.hpp:114
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:787
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:80
Source class definition.
bool m_force_mean
Use mean PSF in any case.
void init_members(void)
Initialise class members.
virtual void clear(void)
Clear response.
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
GLATEdisp * edisp(const int &index) const
Return pointer on energy dispersion.
Class that handles photons.
Definition: GPhoton.hpp:47
Fermi/LAT event bin class.
Calibration database class.
Definition: GCaldb.hpp:66
Interface for the Fermi LAT point spread function.
Definition: GLATPsf.hpp:54
Fermi/LAT event bin class interface definition.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:279
Energy boundaries container class.
Definition: GEbounds.hpp:60
#define G_IRF_ATOM
Filename class.
Definition: GFilename.hpp:62
#define G_PSF
#define G_AEFF
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:229
virtual double irf_ptsrc(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to point source.
GLATResponse(void)
Void constructor.
GChatter
Definition: GTypemaps.hpp:33
Fermi/LAT event cube class definition.
void save(const std::string &rspname) const
Save Fermi LAT response in calibration database.
Interface for the Fermi/LAT effective area.
Definition: GLATAeff.hpp:59
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const
Return value of point source IRF.
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:235
void free_members(void)
Delete class members.
Fermi/LAT event cube class.
const GModelSpatial * model(void) const
Return spatial model component.
Definition: GSource.hpp:128
Abstract observation base class.
const GEnergy & energy(void) const
Return photon energy.
Definition: GSource.hpp:142
GLATAeff * aeff(const int &index) const
Return pointer on effective area.
#define G_EDISP
virtual std::string print(const GChatter &chatter=NORMAL) const
Print Fermi-LAT response information.
virtual GLATResponse * clone(void) const
Clone response.
const double & ontime(void) const
Return ontime of event bin.
Calibration database class interface definition.
std::vector< GLATAeff * > m_aeff
Effective areas.
Point source spatial model.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:249
Fermi/LAT event atom class definition.
void free_members(void)
Delete class members.
Definition: GResponse.cpp:835
#define G_NROI
virtual GResponse & operator=(const GResponse &rsp)
Assignment operator.
Definition: GResponse.cpp:139
double irf_spatial_atom(const GLATEventAtom &event, const GSource &source, const GObservation &obs) const
Return value of model IRF for event atom.
#define G_IRF_BIN
virtual double nroi(const GModelSky &model, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return integral of event probability for a given sky model over ROI.
virtual bool is_atom(void) const =0
Sky model class.
Definition: GModelSky.hpp:122
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
Fermi/LAT instrument direction class definition.
Exception handler interface definition.
#define G_SAVE
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:941
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:955
Fermi/LAT mean PSF class.
Definition: GLATMeanPsf.hpp:51
void enodes(const GNodeArray &enodes)
Set event cube energy nodes.
Fermi LAT Response class definition.
virtual const GInstDir & dir(void) const =0
Abstract instrument response base class.
Definition: GResponse.hpp:77
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:345
std::vector< GLATPsf * > m_psf
Point spread functions.
Sky direction class.
Definition: GSkyDir.hpp:62
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
const double * pixels(void) const
Returns pointer to pixel data.
Definition: GSkyMap.hpp:475
Fermi/LAT event atom class.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::vector< GLATEdisp * > m_edisp
Energy dispersions.
#define G_LOAD
std::string diffname(const int &index) const
Return name of diffuse model.
void copy_members(const GLATResponse &rsp)
Copy class members.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489