GammaLib  2.0.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMResponse.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMResponse.cpp - COMPTEL Response class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2020 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 GCOMResponse.cpp
23  * @brief COMPTEL response class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <string>
32 #include <typeinfo>
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GFits.hpp"
36 #include "GCaldb.hpp"
37 #include "GEvent.hpp"
38 #include "GPhoton.hpp"
39 #include "GSource.hpp"
40 #include "GEnergy.hpp"
41 #include "GTime.hpp"
42 #include "GObservation.hpp"
43 #include "GFitsImage.hpp"
44 #include "GFitsImageFloat.hpp"
45 #include "GModelSky.hpp"
47 #include "GModelSpectralConst.hpp"
48 #include "GCOMResponse.hpp"
49 #include "GCOMObservation.hpp"
50 #include "GCOMEventBin.hpp"
51 #include "GCOMInstDir.hpp"
52 
53 /* __ Method name definitions ____________________________________________ */
54 #define G_IRF "GCOMResponse::irf(GEvent&, GPhoton&, GObservation&)"
55 #define G_IRF_SPATIAL "GCOMResponse::irf_spatial(GEvent&, GSource&, "\
56  "GObservation&)"
57 #define G_NROI "GCOMResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
58  "GObservation&)"
59 #define G_EBOUNDS "GCOMResponse::ebounds(GEnergy&)"
60 #define G_IRF_DRM "GCOMResponse::irf_drm(GModelSky&, GObservation&)"
61 #define G_IRF_RADIAL "GCOMResponse::irf_radial(GModelSky&, GObservation&, "\
62  "GMatrix*)"
63 #define G_IRF_ELLIPTICAL "GCOMResponse::irf_elliptical(GModelSky&, "\
64  " GObservation&, GMatrix*)"
65 #define G_IRF_DIFFUSE "GCOMResponse::irf_diffuse(GModelSky&, GObservation&, "\
66  "GMatrix*)"
67 
68 /* __ Macros _____________________________________________________________ */
69 
70 /* __ Coding definitions _________________________________________________ */
71 #define G_USE_DRM_CUBE
72 
73 /* __ Debug definitions __________________________________________________ */
74 
75 /* __ Constants __________________________________________________________ */
76 
77 
78 /*==========================================================================
79  = =
80  = Constructors/destructors =
81  = =
82  ==========================================================================*/
83 
84 /***********************************************************************//**
85  * @brief Void constructor
86  *
87  * Creates an empty COMPTEL response.
88  ***************************************************************************/
90 {
91  // Initialise members
92  init_members();
93 
94  // Return
95  return;
96 }
97 
98 
99 /***********************************************************************//**
100  * @brief Copy constructor
101  *
102  * @param[in] rsp COM response.
103  **************************************************************************/
105 {
106  // Initialise members
107  init_members();
108 
109  // Copy members
110  copy_members(rsp);
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Response constructor
119  *
120  * @param[in] caldb Calibration database.
121  * @param[in] iaqname IAQ file name.
122  *
123  * Create COMPTEL response by loading an IAQ file from a calibration
124  * database.
125  ***************************************************************************/
127  const std::string& iaqname) : GResponse()
128 {
129  // Initialise members
130  init_members();
131 
132  // Set calibration database
133  this->caldb(caldb);
134 
135  // Load IRF
136  this->load(iaqname);
137 
138  // Return
139  return;
140 }
141 
142 
143 /***********************************************************************//**
144  * @brief Destructor
145  *
146  * Destroys instance of COMPTEL response object.
147  ***************************************************************************/
149 {
150  // Free members
151  free_members();
152 
153  // Return
154  return;
155 }
156 
157 
158 /*==========================================================================
159  = =
160  = Operators =
161  = =
162  ==========================================================================*/
163 
164 /***********************************************************************//**
165  * @brief Assignment operator
166  *
167  * @param[in] rsp COMPTEL response.
168  * @return COMPTEL response.
169  *
170  * Assigns COMPTEL response object to another COMPTEL response object. The
171  * assignment performs a deep copy of all information, hence the original
172  * object from which the assignment has been performed can be destroyed after
173  * this operation without any loss of information.
174  ***************************************************************************/
176 {
177  // Execute only if object is not identical
178  if (this != &rsp) {
179 
180  // Copy base class members
181  this->GResponse::operator=(rsp);
182 
183  // Free members
184  free_members();
185 
186  // Initialise members
187  init_members();
188 
189  // Copy members
190  copy_members(rsp);
191 
192  } // endif: object was not identical
193 
194  // Return this object
195  return *this;
196 }
197 
198 
199 /*==========================================================================
200  = =
201  = Public methods =
202  = =
203  ==========================================================================*/
204 
205 /***********************************************************************//**
206  * @brief Clear instance
207  *
208  * Clears COMPTEL response object by resetting all members to an initial
209  * state. Any information that was present in the object before will be lost.
210  ***************************************************************************/
212 {
213  // Free class members (base and derived classes, derived class first)
214  free_members();
215  this->GResponse::free_members();
216 
217  // Initialise members
218  this->GResponse::init_members();
219  init_members();
220 
221  // Return
222  return;
223 }
224 
225 
226 /***********************************************************************//**
227  * @brief Clone instance
228  *
229  * @return Pointer to deep copy of COMPTEL response.
230  ***************************************************************************/
232 {
233  return new GCOMResponse(*this);
234 }
235 
236 
237 /***********************************************************************//**
238  * @brief Return value of instrument response function
239  *
240  * @param[in] event Observed event.
241  * @param[in] photon Incident photon.
242  * @param[in] obs Observation.
243  * @return Instrument response function (cm^2 sr^-1)
244  *
245  * @exception GException::invalid_argument
246  * Observation is not a COMPTEL observation.
247  * Event is not a COMPTEL event bin.
248  *
249  * Returns the instrument response function for a given observed photon
250  * direction as function of the assumed true photon direction. The result
251  * is given by
252  * \f[IRF = \frac{IAQ \times DRG \times DRX}{ontime \times ewidth}\f]
253  * where
254  * \f$IRF\f$ is the instrument response function,
255  * \f$IAQ\f$ is the COMPTEL response matrix (sr-1),
256  * \f$DRG\f$ is the geometry factor (cm2),
257  * \f$DRX\f$ is the exposure (s),
258  * \f$ontime\f$ is the ontime (s), and
259  * \f$ewidth\f$ is the energy width (MeV).
260  *
261  * The observed photon direction is spanned by the 3 values (Chi,Psi,Phibar).
262  * (Chi,Psi) is the scatter direction of the event, given in sky coordinates.
263  * Phibar is the Compton scatter angle, computed from the energy deposits.
264  ***************************************************************************/
265 double GCOMResponse::irf(const GEvent& event,
266  const GPhoton& photon,
267  const GObservation& obs) const
268 {
269  // Extract COMPTEL observation
270  const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
271  if (observation == NULL) {
272  std::string cls = std::string(typeid(&obs).name());
273  std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
274  "observations. Please specify a COMPTEL observation "
275  "as argument.";
277  }
278 
279  // Extract COMPTEL event bin
280  const GCOMEventBin* bin = dynamic_cast<const GCOMEventBin*>(&event);
281  if (bin == NULL) {
282  std::string cls = std::string(typeid(&event).name());
283  std::string msg = "Event of type \""+cls+"\" is not a COMPTEL event. "
284  "Please specify a COMPTEL event as argument.";
286  }
287 
288  // Extract event parameters
289  const GCOMInstDir& obsDir = bin->dir();
290 
291  // Extract photon parameters
292  const GSkyDir& srcDir = photon.dir();
293  const GTime& srcTime = photon.time();
294 
295  // Compute angle between true photon arrival direction and scatter
296  // direction (Chi,Psi)
297  double phigeo = srcDir.dist_deg(obsDir.dir());
298 
299  // Compute scatter angle index
300  int iphibar = int(obsDir.phibar() / m_phibar_bin_size);
301 
302  // Initialise IRF
303  double iaq = 0.0;
304  double irf = 0.0;
305 
306  // Extract IAQ value by linear inter/extrapolation in Phigeo
307  if (iphibar < m_phibar_bins) {
308  double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
309  int iphigeo = int(phirat); // index into which Phigeo falls
310  double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
311  if (iphigeo < m_phigeo_bins) {
312  int i = iphibar * m_phigeo_bins + iphigeo;
313  if (eps < 0.0) { // interpolate towards left
314  if (iphigeo > 0) {
315  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
316  }
317  else {
318  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
319  }
320  }
321  else { // interpolate towards right
322  if (iphigeo < m_phigeo_bins-1) {
323  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
324  }
325  else {
326  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
327  }
328  }
329  }
330  }
331 
332  // Continue only if IAQ is positive
333  if (iaq > 0.0) {
334 
335  // Get DRG value (units: probability)
336  double drg = observation->drg()(obsDir.dir(), iphibar);
337 
338  // Get DRX value (units: cm^2 sec)
339  double drx = observation->drx()(srcDir);
340 
341  // Get ontime
342  double ontime = observation->ontime(); // sec
343 
344  // Compute IRF value
345  irf = iaq * drg * drx / ontime;
346 
347  // Apply deadtime correction
348  irf *= obs.deadc(srcTime);
349 
350  // Make sure that IRF is positive
351  if (irf < 0.0) {
352  irf = 0.0;
353  }
354 
355  // Compile option: Check for NaN/Inf
356  #if defined(G_NAN_CHECK)
358  std::cout << "*** ERROR: GCOMResponse::irf:";
359  std::cout << " NaN/Inf encountered";
360  std::cout << " (irf=" << irf;
361  std::cout << ", iaq=" << iaq;
362  std::cout << ", drg=" << drg;
363  std::cout << ", drx=" << drx;
364  std::cout << ")";
365  std::cout << std::endl;
366  }
367  #endif
368 
369  } // endif: IAQ was positive
370 
371  // Return IRF value
372  return irf;
373 }
374 
375 
376 /***********************************************************************//**
377  * @brief Return instrument response
378  *
379  * @param[in] event Event.
380  * @param[in] source Source.
381  * @param[in] obs Observation.
382  * @return Instrument response.
383  *
384  * @exception GException::invalid_argument
385  * Invalid observation or event type specified.
386  * @exception GException::feature_not_implemented
387  * Computation for specified spatial model not implemented.
388  *
389  * Returns the instrument response for a given event, source and observation.
390  * If all spatial model parameters are fixed a DRM cube stored in the COMPTEL
391  * observation is used. The computation of this cube is handled by the
392  * GCOMObservation::drm() method that also will update the cube in case that
393  * any of the parameters changes.
394  ***************************************************************************/
395 double GCOMResponse::irf_spatial(const GEvent& event,
396  const GSource& source,
397  const GObservation& obs) const
398 {
399  // Initialise IRF value
400  double irf = 0.0;
401 
402  // Determine if all spatial parameters are fixed
403  #if defined(G_USE_DRM_CUBE)
404  bool fixed = true;
405  for (int i = 0; i < source.model()->size(); ++i) {
406  if ((*(source.model()))[i].is_free()) {
407  fixed = false;
408  break;
409  }
410  }
411 
412  // If all parameters are fixed then use the DRM cache
413  if (fixed) {
414 
415  // Get pointer to COMPTEL observation
416  const GCOMObservation* comobs = dynamic_cast<const GCOMObservation*>(&obs);
417  if (comobs == NULL) {
418  std::string cls = std::string(typeid(&obs).name());
419  std::string msg = "Observation of type \""+cls+"\" is not a "
420  "COMPTEL observation. Please specify a COMPTEL "
421  "observation as argument.";
423  }
424 
425  // Get pointer to COMPTEL event bin
426  if (!event.is_bin()) {
427  std::string msg = "The current event is not a COMPTEL event bin. "
428  "This method only works on binned COMPTEL data. "
429  "Please make sure that a COMPTEL observation "
430  "containing binned data is provided.";
432  }
433  const GCOMEventBin* bin = static_cast<const GCOMEventBin*>(&event);
434 
435  // Setup sky model for DRM computation
436  GModelSky model(*(source.model()), GModelSpectralConst());
437 
438  // Get IRF value
439  irf = comobs->drm(model)[bin->index()];
440 
441  } // endif: all spatial parameters were fixed
442 
443  // ... otherwise compute IRF directly
444  else {
445  #endif
446 
447  // Select IRF depending on the spatial model type
448  switch (source.model()->code()) {
450  {
451  const GModelSpatialPointSource* src =
452  static_cast<const GModelSpatialPointSource*>(source.model());
453  GPhoton photon(src->dir(), source.energy(), source.time());
454  irf = this->irf(event, photon, obs);
455  }
456  break;
460  {
461  std::string msg = "Response computation not yet implemented "
462  "for spatial model type \""+
463  source.model()->type()+"\".";
465  }
466  break;
467  default:
468  break;
469  }
470 
471  #if defined(G_USE_DRM_CUBE)
472  } // endelse: used IRF computation directly
473  #endif
474 
475  // Return IRF value
476  return irf;
477 }
478 
479 
480 /***********************************************************************//**
481  * @brief Return integral of event probability for a given sky model over ROI
482  *
483  * @param[in] model Sky model.
484  * @param[in] obsEng Observed photon energy.
485  * @param[in] obsTime Observed photon arrival time.
486  * @param[in] obs Observation.
487  * @return 0.0
488  *
489  * @exception GException::feature_not_implemented
490  * Method is not implemented.
491  ***************************************************************************/
492 double GCOMResponse::nroi(const GModelSky& model,
493  const GEnergy& obsEng,
494  const GTime& obsTime,
495  const GObservation& obs) const
496 {
497  // Method is not implemented
498  std::string msg = "Spatial integration of sky model over the data space "
499  "is not implemented.";
501 
502  // Return Npred
503  return (0.0);
504 }
505 
506 
507 /***********************************************************************//**
508  * @brief Return true energy boundaries for a specific observed energy
509  *
510  * @param[in] obsEnergy Observed Energy.
511  * @return True energy boundaries for given observed energy.
512  *
513  * @exception GException::feature_not_implemented
514  * Method is not implemented.
515  ***************************************************************************/
516 GEbounds GCOMResponse::ebounds(const GEnergy& obsEnergy) const
517 {
518  // Initialise an empty boundary object
520 
521  // Throw an exception
522  std::string msg = "Energy dispersion not implemented.";
524 
525  // Return energy boundaries
526  return ebounds;
527 }
528 
529 
530 /***********************************************************************//**
531  * @brief Load COMPTEL response.
532  *
533  * @param[in] rspname COMPTEL response name.
534  *
535  * Loads the COMPTEL response with specified name @p rspname.
536  *
537  * The method first attempts to interpret @p rspname as a file name and to
538  * load the corresponding response.
539  *
540  * If @p rspname is not a FITS file the method searches for an appropriate
541  * response in the calibration database. If no appropriate response is found,
542  * the method takes the CALDB root path and response name to build the full
543  * path to the response file, and tries to load the response from these
544  * paths.
545  *
546  * If also this fails an exception is thrown.
547  ***************************************************************************/
548 void GCOMResponse::load(const std::string& rspname)
549 {
550  // Clear instance but conserve calibration database
551  GCaldb caldb = m_caldb;
552  clear();
553  m_caldb = caldb;
554 
555  // Save response name
556  m_rspname = rspname;
557 
558  // Interpret response name as a FITS file name
559  GFilename filename(rspname);
560 
561  // If the filename does not exist the try getting the response from the
562  // calibration database
563  if (!filename.is_fits()) {
564 
565  // Get GCaldb response
566  filename = m_caldb.filename("","","IAQ","","",rspname);
567 
568  // If filename is empty then build filename from CALDB root path
569  // and response name
570  if (filename.is_empty()) {
571  filename = gammalib::filepath(m_caldb.rootdir(), rspname);
572  if (!filename.exists()) {
573  GFilename testname = filename + ".fits";
574  if (testname.exists()) {
575  filename = testname;
576  }
577  }
578  }
579 
580  } // endif: response name is not a FITS file
581 
582  // Open FITS file
583  GFits fits(filename);
584 
585  // Get IAQ image
586  const GFitsImage& iaq = *fits.image(0);
587 
588  // Read IAQ
589  read(iaq);
590 
591  // Close IAQ FITS file
592  fits.close();
593 
594  // Return
595  return;
596 }
597 
598 
599 /***********************************************************************//**
600  * @brief Read COMPTEL response from FITS image.
601  *
602  * @param[in] image FITS image.
603  *
604  * Read the COMPTEL response from IAQ FITS file and convert the IAQ values
605  * into a probability per steradian.
606  ***************************************************************************/
607 void GCOMResponse::read(const GFitsImage& image)
608 {
609  // Continue only if there are two image axis
610  if (image.naxis() == 2) {
611 
612  // Store IAQ dimensions
613  m_phigeo_bins = image.naxes(0);
614  m_phibar_bins = image.naxes(1);
615 
616  // Store IAQ axes definitions
617  m_phigeo_ref_value = image.real("CRVAL1");
618  m_phigeo_ref_pixel = image.real("CRPIX1");
619  m_phigeo_bin_size = image.real("CDELT1");
620  m_phibar_ref_value = image.real("CRVAL2");
621  m_phibar_ref_pixel = image.real("CRPIX2");
622  m_phibar_bin_size = image.real("CDELT2");
623 
624  // Get axes minima (values of first bin)
627 
628  // Compute IAQ size. Continue only if size is positive
629  int size = m_phigeo_bins * m_phibar_bins;
630  if (size > 0) {
631 
632  // Allocate memory for IAQ
633  m_iaq.assign(size, 0.0);
634 
635  // Copy over IAQ values
636  for (int i = 0; i < size; ++i) {
637  m_iaq[i] = image.pixel(i);
638  }
639 
640  } // endif: size was positive
641 
642  // Convert IAQ matrix from probability per Phigeo bin into a
643  // probability per steradian
644  double omega0 = gammalib::fourpi *
646  for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
647  double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
648  double omega = omega0 * std::sin(phigeo * gammalib::deg2rad);
649  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
650  m_iaq[iphigeo+iphibar*m_phigeo_bins] /= omega;
651  }
652  }
653 
654  } // endif: image had two axes
655 
656  // Return
657  return;
658 }
659 
660 
661 /***********************************************************************//**
662  * @brief Write COMPTEL response into FITS image.
663  *
664  * @param[in] image FITS image.
665  *
666  * Writes the COMPTEL response into an IAQ FITS file.
667  ***************************************************************************/
669 {
670  // Continue only if response is not empty
671  if (m_phigeo_bins > 0 && m_phibar_bins > 0) {
672 
673  // Initialise image
675 
676  // Convert IAQ matrix from probability per steradian into
677  //probability per Phigeo bin
678  double omega0 = gammalib::fourpi *
680  for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
681  double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
682  double omega = omega0 * std::sin(phigeo * gammalib::deg2rad);
683  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
684  int i = iphigeo + iphibar*m_phigeo_bins;
685  image(iphigeo, iphibar) = m_iaq[i] * omega;
686  }
687  }
688 
689  // Set header keywords
690  image.card("CTYPE1", "Phigeo", "Geometrical scatter angle");
691  image.card("CRVAL1", m_phigeo_ref_value,
692  "[deg] Geometrical scatter angle for reference bin");
693  image.card("CDELT1", m_phigeo_bin_size,
694  "[deg] Geometrical scatter angle bin size");
695  image.card("CRPIX1", m_phigeo_ref_pixel,
696  "Reference bin of geometrical scatter angle");
697  image.card("CTYPE2", "Phibar", "Compton scatter angle");
698  image.card("CRVAL2", m_phibar_ref_value,
699  "[deg] Compton scatter angle for reference bin");
700  image.card("CDELT2", m_phibar_bin_size, "[deg] Compton scatter angle bin size");
701  image.card("CRPIX2", m_phibar_ref_pixel,
702  "Reference bin of Compton scatter angle");
703  image.card("BUNIT", "Probability per bin", "Unit of bins");
704  image.card("TELESCOP", "GRO", "Name of telescope");
705  image.card("INSTRUME", "COMPTEL", "Name of instrument");
706  image.card("ORIGIN", "GammaLib", "Origin of FITS file");
707  image.card("OBSERVER", "Unknown", "Observer that created FITS file");
708 
709  } // endif: response was not empty
710 
711  // Return
712  return;
713 }
714 
715 
716 /***********************************************************************//**
717  * @brief Print COMPTEL response information
718  *
719  * @param[in] chatter Chattiness.
720  * @return String containing COMPTEL response information.
721  ***************************************************************************/
722 std::string GCOMResponse::print(const GChatter& chatter) const
723 {
724  // Initialise result string
725  std::string result;
726 
727  // Continue only if chatter is not silent
728  if (chatter != SILENT) {
729 
730  // Append header
731  result.append("=== GCOMResponse ===");
732 
733  // Append response name
734  result.append("\n"+gammalib::parformat("Response name")+m_rspname);
735 
736  // EXPLICIT: Append detailed information
737  if (chatter >= EXPLICIT) {
738 
739  // Append information
740  result.append("\n"+gammalib::parformat("Number of Phigeo bins"));
741  result.append(gammalib::str(m_phigeo_bins));
742  result.append("\n"+gammalib::parformat("Number of Phibar bins"));
743  result.append(gammalib::str(m_phibar_bins));
744  result.append("\n"+gammalib::parformat("Phigeo reference value"));
745  result.append(gammalib::str(m_phigeo_ref_value)+" deg");
746  result.append("\n"+gammalib::parformat("Phigeo reference pixel"));
747  result.append(gammalib::str(m_phigeo_ref_pixel));
748  result.append("\n"+gammalib::parformat("Phigeo bin size"));
749  result.append(gammalib::str(m_phigeo_bin_size)+" deg");
750  result.append("\n"+gammalib::parformat("Phigeo first bin value"));
751  result.append(gammalib::str(m_phigeo_min)+" deg");
752  result.append("\n"+gammalib::parformat("Phibar reference value"));
753  result.append(gammalib::str(m_phibar_ref_value)+" deg");
754  result.append("\n"+gammalib::parformat("Phibar reference pixel"));
755  result.append(gammalib::str(m_phibar_ref_pixel));
756  result.append("\n"+gammalib::parformat("Phibar bin size"));
757  result.append(gammalib::str(m_phibar_bin_size)+" deg");
758  result.append("\n"+gammalib::parformat("Phibar first bin value"));
759  result.append(gammalib::str(m_phibar_min)+" deg");
760 
761  }
762 
763  // VERBOSE: Append calibration database
764  if (chatter == VERBOSE) {
765  result.append("\n"+m_caldb.print(chatter));
766  }
767 
768  } // endif: chatter was not silent
769 
770  // Return result
771  return result;
772 }
773 
774 
775 /*==========================================================================
776  = =
777  = Private methods =
778  = =
779  ==========================================================================*/
780 
781 /***********************************************************************//**
782  * @brief Initialise class members
783  ***************************************************************************/
785 {
786  // Initialise members
787  m_caldb.clear();
788  m_rspname.clear();
789  m_iaq.clear();
790  m_phigeo_bins = 0;
791  m_phibar_bins = 0;
792  m_phigeo_ref_value = 0.0;
793  m_phigeo_ref_pixel = 0.0;
794  m_phigeo_bin_size = 0.0;
795  m_phigeo_min = 0.0;
796  m_phibar_ref_value = 0.0;
797  m_phibar_ref_pixel = 0.0;
798  m_phibar_bin_size = 0.0;
799  m_phibar_min = 0.0;
800 
801  // Return
802  return;
803 }
804 
805 
806 /***********************************************************************//**
807  * @brief Copy class members
808  *
809  * @param[in] rsp COMPTEL response.
810  ***************************************************************************/
812 {
813  // Copy attributes
814  m_caldb = rsp.m_caldb;
815  m_rspname = rsp.m_rspname;
816  m_iaq = rsp.m_iaq;
827 
828  // Return
829  return;
830 }
831 
832 
833 /***********************************************************************//**
834  * @brief Delete class members
835  ***************************************************************************/
837 {
838  // Return
839  return;
840 }
841 
842 
843 /***********************************************************************//**
844  * @brief Return instrument response
845  *
846  * @param[in] event Event.
847  * @param[in] source Source.
848  * @param[in] obs Observation.
849  * @return Instrument response.
850  *
851  * @exception GException::invalid_argument
852  * Invalid observation or event type specified.
853  * @exception GException::feature_not_implemented
854  * Computation for specified spatial model not implemented.
855  *
856  * Returns the instrument response for a given event, source and observation.
857  * If all spatial model parameters are fixed a DRM cube stored in the COMPTEL
858  * observation is used. The computation of this cube is handled by the
859  * GCOMObservation::drm() method that also will update the cube in case that
860  * any of the parameters changes.
861  ***************************************************************************/
863  const GObservation& obs) const
864 {
865  // Get number of events
866  int nevents = obs.events()->size();
867 
868  // Initialise result
869  GVector irfs(nevents);
870 
871  // Get pointer to COMPTEL observation
872  const GCOMObservation* comobs = dynamic_cast<const GCOMObservation*>(&obs);
873  if (comobs == NULL) {
874  std::string cls = std::string(typeid(&obs).name());
875  std::string msg = "Observation of type \""+cls+"\" is not a "
876  "COMPTEL observation. Please specify a COMPTEL "
877  "observation as argument.";
879  }
880 
881  // Get reference to DRM
882  const GCOMDri& drm = comobs->drm(model);
883 
884  // Extract model values from DRM cache
885  for (int k = 0; k < nevents; ++k) {
886  irfs[k] = drm[k];
887  }
888 
889  // Return
890  return irfs;
891 }
892 
893 
894 /***********************************************************************//**
895  * @brief Return instrument response to point source sky model
896  *
897  * @param[in] model Sky model.
898  * @param[in] obs Observation.
899  * @param[in] gradients Gradients matrix.
900  * @return Instrument response to point source sky model.
901  *
902  * Returns the instrument response to a point source sky model for all
903  * events.
904  ***************************************************************************/
906  const GObservation& obs,
907  GMatrix* gradients) const
908 {
909  // Get number of events
910  int nevents = obs.events()->size();
911 
912  // Initialise result
913  GVector irfs(nevents);
914 
915  // Check if spatial model has free parameters
916  bool has_free_pars = model.spatial()->has_free_pars();
917 
918  // If model has free parameters then compute IRF value for all events
919  if (has_free_pars) {
920 
921  // Get point source direction
922  GSkyDir srcDir =
923  static_cast<const GModelSpatialPointSource*>(model.spatial())->dir();
924 
925  // Loop over events
926  for (int k = 0; k < nevents; ++k) {
927 
928  // Get reference to event
929  const GEvent& event = *((*obs.events())[k]);
930 
931  // Get source energy and time (no dispersion)
932  GEnergy srcEng = event.energy();
933  GTime srcTime = event.time();
934 
935  // Setup photon
936  GPhoton photon(srcDir, srcEng, srcTime);
937 
938  // Get IRF value
939  irfs[k] = this->irf(event, photon, obs);
940 
941  } // endfor: looped over events
942 
943  } // endif: model had free parameters
944 
945  // ... otherwise extract model from DRM cache
946  else {
947  irfs = irf_drm(model, obs);
948  }
949 
950  // Return IRF value
951  return irfs;
952 }
953 
954 
955 /***********************************************************************//**
956  * @brief Return instrument response to radial source sky model
957  *
958  * @param[in] model Sky model.
959  * @param[in] obs Observation.
960  * @param[in] gradients Gradients matrix.
961  * @return Instrument response to radial source sky model.
962  *
963  * @todo Implement method.
964  ***************************************************************************/
966  const GObservation& obs,
967  GMatrix* gradients) const
968 {
969  // Throw exception
970  std::string msg = "Response computation not yet implemented "
971  "for spatial model type \""+
972  model.spatial()->type()+"\".";
974 
975  // Get number of events
976  int nevents = obs.events()->size();
977 
978  // Initialise result
979  GVector irfs(nevents);
980 
981  // Return IRF value
982  return irfs;
983 }
984 
985 
986 /***********************************************************************//**
987  * @brief Return instrument response to elliptical source sky model
988  *
989  * @param[in] model Sky model.
990  * @param[in] obs Observation.
991  * @param[in] gradients Gradients matrix.
992  * @return Instrument response to elliptical source sky model.
993  *
994  * @todo Implement method.
995  ***************************************************************************/
997  const GObservation& obs,
998  GMatrix* gradients) const
999 {
1000  // Throw exception
1001  std::string msg = "Response computation not yet implemented "
1002  "for spatial model type \""+
1003  model.spatial()->type()+"\".";
1005 
1006  // Get number of events
1007  int nevents = obs.events()->size();
1008 
1009  // Initialise result
1010  GVector irfs(nevents);
1011 
1012  // Return IRF value
1013  return irfs;
1014 }
1015 
1016 
1017 /***********************************************************************//**
1018  * @brief Return instrument response to diffuse source sky model
1019  *
1020  * @param[in] model Sky model.
1021  * @param[in] obs Observation.
1022  * @param[in] gradients Gradients matrix.
1023  * @return Instrument response to diffuse source sky model.
1024  *
1025  * @todo Implement method.
1026  ***************************************************************************/
1028  const GObservation& obs,
1029  GMatrix* gradients) const
1030 {
1031  // Throw exception
1032  std::string msg = "Response computation not yet implemented "
1033  "for spatial model type \""+
1034  model.spatial()->type()+"\".";
1036 
1037  // Get number of events
1038  int nevents = obs.events()->size();
1039 
1040  // Initialise result
1041  GVector irfs(nevents);
1042 
1043  // Return IRF value
1044  return irfs;
1045 }
virtual GCOMResponse * clone(void) const
Clone instance.
void phibar(const double &phibar)
Set event Compton scatter angle.
virtual GVector irf_radial(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to radial source sky model.
double m_phibar_ref_value
Phigeo reference value (deg)
const std::string & rspname(void) const
Return response name.
GCaldb m_caldb
Calibration database.
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:284
#define G_EBOUNDS
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
const int & index(void) const
Return bin index.
const int & naxis(void) const
Return dimension of image.
Definition: GFitsImage.hpp:156
virtual std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL response information.
virtual void clear(void)
Clear instance.
Point source spatial model class interface definition.
double m_phigeo_bin_size
Phigeo binsize (deg)
Energy value class definition.
int m_phibar_bins
Number of Phibar bins.
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.
std::string m_rspname
Response name.
std::vector< double > m_iaq
IAQ array.
int size(void) const
Return number of parameters.
Single precision FITS image class definition.
bool has_free_pars(void) const
Checks if the spatial model has free parameters.
double m_phibar_min
Phigeo value of first bin (deg)
#define G_IRF_DIFFUSE
void init_members(void)
Initialise class members.
virtual double ontime(void) const
Return ontime.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
COMPTEL observation class interface definition.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
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
virtual GVector irf_ptsrc(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to point source sky model.
#define G_IRF_DRM
Interface for the COMPTEL instrument response function.
double m_phigeo_ref_value
Phigeo reference value (deg)
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
virtual double irf_spatial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response.
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:754
COMPTEL event bin class.
double m_phibar_ref_pixel
Phigeo reference pixel (starting from 1)
Source class definition.
virtual GClassCode code(void) const =0
COMPTEL event bin class interface definition.
Abstract FITS image base class definition.
double m_phigeo_min
Phigeo value of first bin (deg)
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
virtual GCOMResponse & operator=(const GCOMResponse &rsp)
Assignment operator.
Class that handles photons.
Definition: GPhoton.hpp:47
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
virtual bool is_bin(void) const =0
const double deg2rad
Definition: GMath.hpp:43
Calibration database class.
Definition: GCaldb.hpp:66
Energy boundaries container class.
Definition: GEbounds.hpp:60
COMPTEL instrument response function class interface definition.
void dir(const GSkyDir &dir)
Set event scatter direction.
Filename class.
Definition: GFilename.hpp:62
bool exists(void) const
Checks whether file exists.
Definition: GFilename.cpp:223
Constant spectral model class interface definition.
virtual double pixel(const int &ix) const =0
Abstract event base class definition.
void write(GFitsImageFloat &image) const
Write COMPTEL response into FITS image.
virtual double deadc(const GTime &time=GTime()) const =0
virtual GVector irf_diffuse(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to diffuse source sky model.
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
const GTime & time(void) const
Return photon arrival time.
Definition: GSource.hpp:156
GChatter
Definition: GTypemaps.hpp:33
const GCOMDri & drm(const GModelSky &model) const
Return model DRM cube.
bool is_fits(void) const
Checks whether file is a FITS file.
Definition: GFilename.cpp:313
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
std::string type(void) const
Return model type.
Single precision FITS image class.
Photon class definition.
Abstract observation base class interface definition.
#define G_IRF_SPATIAL
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const
Return value of instrument response function.
Calibration database class interface definition.
Point source spatial model.
void free_members(void)
Delete class members.
Definition: GResponse.cpp:802
int naxes(const int &axis) const
Return dimension of an image axis.
Definition: GFitsImage.cpp:344
virtual GResponse & operator=(const GResponse &rsp)
Assignment operator.
Definition: GResponse.cpp:138
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
const GSkyMap & drx(void) const
Return exposure.
const GCaldb & caldb(void) const
Return calibration database.
void free_members(void)
Delete class members.
double m_phigeo_ref_pixel
Phigeo reference pixel (starting from 1)
GSkyDir dir(void) const
Return position of point source.
Sky model class interface definition.
Sky model class.
Definition: GModelSky.hpp:121
#define G_IRF_ELLIPTICAL
virtual const GCOMInstDir & dir(void) const
Return instrument direction of event bin.
virtual GEbounds ebounds(const GEnergy &obsEnergy) const
Return true energy boundaries for a specific observed energy.
GCOMResponse(void)
Void constructor.
Interface class for COMPTEL observations.
#define G_NROI
int m_phigeo_bins
Number of Phigeo bins.
double m_phibar_bin_size
Phigeo binsize (deg)
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1233
virtual GEvents * events(void)
Return events.
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:640
COMPTEL Data Space class.
Definition: GCOMDri.hpp:60
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:250
const double fourpi
Definition: GMath.hpp:37
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
std::string filepath(const std::string &pathname, const std::string &filename)
Build file path from path name and file name.
Definition: GTools.cpp:388
virtual GVector irf_elliptical(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to elliptical source sky model.
void copy_members(const GCOMResponse &rsp)
Copy class members.
COMPTEL instrument direction class definition.
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:1033
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
virtual ~GCOMResponse(void)
Destructor.
#define G_IRF
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Sky direction class.
Definition: GSkyDir.hpp:62
GVector irf_drm(const GModelSky &model, const GObservation &obs) const
Return instrument response.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
#define G_IRF_RADIAL
Time class interface definition.
void clear(void)
Clear calibration database.
Definition: GCaldb.cpp:194
Constant spectral model class.
void read(const GFitsImage &hdu)
Read COMPTEL response from FITS image.
std::string rootdir(void) const
Return path to CALDB root directory.
Definition: GCaldb.cpp:259
Interface for the COMPTEL instrument direction class.
Definition: GCOMInstDir.hpp:45
Mathematical function definitions.
const GSkyMap & drg(void) const
Return geometry factors.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
virtual int size(void) const =0
void load(const std::string &rspname)
Load COMPTEL response.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415