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