GammaLib  1.7.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-2018 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 "GCOMResponse.hpp"
48 #include "GCOMObservation.hpp"
49 #include "GCOMEventBin.hpp"
50 #include "GCOMInstDir.hpp"
51 
52 /* __ Method name definitions ____________________________________________ */
53 #define G_IRF "GCOMResponse::irf(GEvent&, GSource&, GObservation&)"
54 #define G_NROI "GCOMResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
55  "GObservation&)"
56 #define G_EBOUNDS "GCOMResponse::ebounds(GEnergy&)"
57 
58 /* __ Macros _____________________________________________________________ */
59 
60 /* __ Coding definitions _________________________________________________ */
61 #define G_USE_DRM_CUBE
62 
63 /* __ Debug definitions __________________________________________________ */
64 
65 /* __ Constants __________________________________________________________ */
66 
67 
68 /*==========================================================================
69  = =
70  = Constructors/destructors =
71  = =
72  ==========================================================================*/
73 
74 /***********************************************************************//**
75  * @brief Void constructor
76  *
77  * Creates an empty COMPTEL response.
78  ***************************************************************************/
80 {
81  // Initialise members
82  init_members();
83 
84  // Return
85  return;
86 }
87 
88 
89 /***********************************************************************//**
90  * @brief Copy constructor
91  *
92  * @param[in] rsp COM response.
93  **************************************************************************/
95 {
96  // Initialise members
97  init_members();
98 
99  // Copy members
100  copy_members(rsp);
101 
102  // Return
103  return;
104 }
105 
106 
107 /***********************************************************************//**
108  * @brief Response constructor
109  *
110  * @param[in] caldb Calibration database.
111  * @param[in] iaqname IAQ file name.
112  *
113  * Create COMPTEL response by loading an IAQ file from a calibration
114  * database.
115  ***************************************************************************/
117  const std::string& iaqname) : GResponse()
118 {
119  // Initialise members
120  init_members();
121 
122  // Set calibration database
123  this->caldb(caldb);
124 
125  // Load IRF
126  this->load(iaqname);
127 
128  // Return
129  return;
130 }
131 
132 
133 /***********************************************************************//**
134  * @brief Destructor
135  *
136  * Destroys instance of COMPTEL response object.
137  ***************************************************************************/
139 {
140  // Free members
141  free_members();
142 
143  // Return
144  return;
145 }
146 
147 
148 /*==========================================================================
149  = =
150  = Operators =
151  = =
152  ==========================================================================*/
153 
154 /***********************************************************************//**
155  * @brief Assignment operator
156  *
157  * @param[in] rsp COMPTEL response.
158  * @return COMPTEL response.
159  *
160  * Assigns COMPTEL response object to another COMPTEL response object. The
161  * assignment performs a deep copy of all information, hence the original
162  * object from which the assignment has been performed can be destroyed after
163  * this operation without any loss of information.
164  ***************************************************************************/
166 {
167  // Execute only if object is not identical
168  if (this != &rsp) {
169 
170  // Copy base class members
171  this->GResponse::operator=(rsp);
172 
173  // Free members
174  free_members();
175 
176  // Initialise members
177  init_members();
178 
179  // Copy members
180  copy_members(rsp);
181 
182  } // endif: object was not identical
183 
184  // Return this object
185  return *this;
186 }
187 
188 
189 /*==========================================================================
190  = =
191  = Public methods =
192  = =
193  ==========================================================================*/
194 
195 /***********************************************************************//**
196  * @brief Clear instance
197  *
198  * Clears COMPTEL response object by resetting all members to an initial
199  * state. Any information that was present in the object before will be lost.
200  ***************************************************************************/
202 {
203  // Free class members (base and derived classes, derived class first)
204  free_members();
205  this->GResponse::free_members();
206 
207  // Initialise members
208  this->GResponse::init_members();
209  init_members();
210 
211  // Return
212  return;
213 }
214 
215 
216 /***********************************************************************//**
217  * @brief Clone instance
218  *
219  * @return Pointer to deep copy of COMPTEL response.
220  ***************************************************************************/
222 {
223  return new GCOMResponse(*this);
224 }
225 
226 
227 /***********************************************************************//**
228  * @brief Return value of instrument response function
229  *
230  * @param[in] event Observed event.
231  * @param[in] photon Incident photon.
232  * @param[in] obs Observation.
233  * @return Instrument response function (cm^2 sr^-1)
234  *
235  * @exception GException::invalid_argument
236  * Observation is not a COMPTEL observation.
237  * Event is not a COMPTEL event bin.
238  *
239  * Returns the instrument response function for a given observed photon
240  * direction as function of the assumed true photon direction. The result
241  * is given by
242  * \f[IRF = \frac{IAQ \times DRG \times DRX}{ontime \times ewidth}\f]
243  * where
244  * \f$IRF\f$ is the instrument response function,
245  * \f$IAQ\f$ is the COMPTEL response matrix (sr-1),
246  * \f$DRG\f$ is the geometry factor (cm2),
247  * \f$DRX\f$ is the exposure (s),
248  * \f$ontime\f$ is the ontime (s), and
249  * \f$ewidth\f$ is the energy width (MeV).
250  *
251  * The observed photon direction is spanned by the 3 values (Chi,Psi,Phibar).
252  * (Chi,Psi) is the scatter direction of the event, given in sky coordinates.
253  * Phibar is the Compton scatter angle, computed from the energy deposits.
254  ***************************************************************************/
255 double GCOMResponse::irf(const GEvent& event,
256  const GPhoton& photon,
257  const GObservation& obs) const
258 {
259  // Extract COMPTEL observation
260  const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
261  if (observation == NULL) {
262  std::string cls = std::string(typeid(&obs).name());
263  std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
264  "observations. Please specify a COMPTEL observation "
265  "as argument.";
267  }
268 
269  // Extract COMPTEL event bin
270  const GCOMEventBin* bin = dynamic_cast<const GCOMEventBin*>(&event);
271  if (bin == NULL) {
272  std::string cls = std::string(typeid(&event).name());
273  std::string msg = "Event of type \""+cls+"\" is not a COMPTEL event. "
274  "Please specify a COMPTEL event as argument.";
276  }
277 
278  // Extract event parameters
279  const GCOMInstDir& obsDir = bin->dir();
280 
281  // Extract photon parameters
282  const GSkyDir& srcDir = photon.dir();
283  const GTime& srcTime = photon.time();
284 
285  // Compute angle between true photon arrival direction and scatter
286  // direction (Chi,Psi)
287  double phigeo = srcDir.dist_deg(obsDir.dir());
288 
289  // Compute scatter angle index
290  int iphibar = int(obsDir.phibar() / m_phibar_bin_size);
291 
292  // Initialise IRF
293  double iaq = 0.0;
294  double irf = 0.0;
295 
296  // Extract IAQ value by linear inter/extrapolation in Phigeo
297  if (iphibar < m_phibar_bins) {
298  double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
299  int iphigeo = int(phirat); // index into which Phigeo falls
300  double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
301  if (iphigeo < m_phigeo_bins) {
302  int i = iphibar * m_phigeo_bins + iphigeo;
303  if (eps < 0.0) { // interpolate towards left
304  if (iphigeo > 0) {
305  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
306  }
307  else {
308  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
309  }
310  }
311  else { // interpolate towards right
312  if (iphigeo < m_phigeo_bins-1) {
313  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
314  }
315  else {
316  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
317  }
318  }
319  }
320  }
321 
322  // Continue only if IAQ is positive
323  if (iaq > 0.0) {
324 
325  // Get DRG value (units: probability)
326  double drg = observation->drg()(obsDir.dir(), iphibar);
327 
328  // Get DRX value (units: cm^2 sec)
329  double drx = observation->drx()(srcDir);
330 
331  // Get ontime
332  double ontime = observation->ontime(); // sec
333 
334  // Compute IRF value
335  irf = iaq * drg * drx / ontime;
336 
337  // Apply deadtime correction
338  irf *= obs.deadc(srcTime);
339 
340  // Make sure that IRF is positive
341  if (irf < 0.0) {
342  irf = 0.0;
343  }
344 
345  // Compile option: Check for NaN/Inf
346  #if defined(G_NAN_CHECK)
348  std::cout << "*** ERROR: GCOMResponse::irf:";
349  std::cout << " NaN/Inf encountered";
350  std::cout << " (irf=" << irf;
351  std::cout << ", iaq=" << iaq;
352  std::cout << ", drg=" << drg;
353  std::cout << ", drx=" << drx;
354  std::cout << ")";
355  std::cout << std::endl;
356  }
357  #endif
358 
359  } // endif: IAQ was positive
360 
361  // Return IRF value
362  return irf;
363 }
364 
365 
366 /***********************************************************************//**
367  * @brief Return instrument response
368  *
369  * @param[in] event Event.
370  * @param[in] source Source.
371  * @param[in] obs Observation.
372  * @return Instrument response.
373  *
374  * @exception GException::invalid_argument
375  * Invalid observation or event type specified.
376  * @exception GException::feature_not_implemented
377  * Computation for specified spatial model not implemented.
378  *
379  * Returns the instrument response for a given event, source and observation.
380  * If all spatial model parameters are fixed a DRM cube stored in the COMPTEL
381  * observation is used. The computation of this cube is handled by the
382  * GCOMObservation::drm() method that also will update the cube in case that
383  * any of the parameters changes.
384  ***************************************************************************/
385 double GCOMResponse::irf(const GEvent& event,
386  const GSource& source,
387  const GObservation& obs) const
388 {
389  // Initialise IRF value
390  double irf = 0.0;
391 
392  // Determine if all spatial parameters are fixed
393  #if defined(G_USE_DRM_CUBE)
394  bool fixed = true;
395  for (int i = 0; i < source.model()->size(); ++i) {
396  if ((*(source.model()))[i].is_free()) {
397  fixed = false;
398  break;
399  }
400  }
401 
402  // If all parameters are fixed then use the DRM cache
403  if (fixed) {
404 
405  // Get pointer to COMPTEL observation
406  const GCOMObservation* comobs = dynamic_cast<const GCOMObservation*>(&obs);
407  if (comobs == NULL) {
408  std::string cls = std::string(typeid(&obs).name());
409  std::string msg = "Observation of type \""+cls+"\" is not a "
410  "COMPTEL observation. Please specify a COMPTEL "
411  "observation as argument.";
413  }
414 
415  // Get pointer to COMPTEL event bin
416  if (!event.is_bin()) {
417  std::string msg = "The current event is not a COMPTEL event bin. "
418  "This method only works on binned COMPTEL data. "
419  "Please make sure that a COMPTEL observation "
420  "containing binned data is provided.";
422  }
423  const GCOMEventBin* bin = static_cast<const GCOMEventBin*>(&event);
424 
425  // Get IRF value
426  irf = comobs->drm(source)[bin->index()];
427 
428  } // endif: all spatial parameters were fixed
429 
430  // ... otherwise compute IRF directly
431  else {
432  #endif
433 
434  // Select IRF depending on the spatial model type
435  switch (source.model()->code()) {
437  {
438  const GModelSpatialPointSource* src =
439  static_cast<const GModelSpatialPointSource*>(source.model());
440  GPhoton photon(src->dir(), source.energy(), source.time());
441  irf = this->irf(event, photon, obs);
442  }
443  break;
447  {
448  std::string msg = "Response computation not yet implemented "
449  "for spatial model type \""+
450  source.model()->type()+"\".";
452  }
453  break;
454  default:
455  break;
456  }
457 
458  #if defined(G_USE_DRM_CUBE)
459  } // endelse: used IRF computation directly
460  #endif
461 
462  // Return IRF value
463  return irf;
464 }
465 
466 
467 /***********************************************************************//**
468  * @brief Return integral of event probability for a given sky model over ROI
469  *
470  * @param[in] model Sky model.
471  * @param[in] obsEng Observed photon energy.
472  * @param[in] obsTime Observed photon arrival time.
473  * @param[in] obs Observation.
474  * @return 0.0
475  *
476  * @exception GException::feature_not_implemented
477  * Method is not implemented.
478  ***************************************************************************/
479 double GCOMResponse::nroi(const GModelSky& model,
480  const GEnergy& obsEng,
481  const GTime& obsTime,
482  const GObservation& obs) const
483 {
484  // Method is not implemented
485  std::string msg = "Spatial integration of sky model over the data space "
486  "is not implemented.";
488 
489  // Return Npred
490  return (0.0);
491 }
492 
493 
494 /***********************************************************************//**
495  * @brief Return true energy boundaries for a specific observed energy
496  *
497  * @param[in] obsEnergy Observed Energy.
498  * @return True energy boundaries for given observed energy.
499  *
500  * @exception GException::feature_not_implemented
501  * Method is not implemented.
502  ***************************************************************************/
503 GEbounds GCOMResponse::ebounds(const GEnergy& obsEnergy) const
504 {
505  // Initialise an empty boundary object
507 
508  // Throw an exception
509  std::string msg = "Energy dispersion not implemented.";
511 
512  // Return energy boundaries
513  return ebounds;
514 }
515 
516 
517 /***********************************************************************//**
518  * @brief Load COMPTEL response.
519  *
520  * @param[in] rspname COMPTEL response name.
521  *
522  * Loads the COMPTEL response with specified name @p rspname.
523  *
524  * The method first attempts to interpret @p rspname as a file name and to
525  * load the corresponding response.
526  *
527  * If @p rspname is not a FITS file the method searches for an appropriate
528  * response in the calibration database. If no appropriate response is found,
529  * the method takes the CALDB root path and response name to build the full
530  * path to the response file, and tries to load the response from these
531  * paths.
532  *
533  * If also this fails an exception is thrown.
534  ***************************************************************************/
535 void GCOMResponse::load(const std::string& rspname)
536 {
537  // Clear instance but conserve calibration database
538  GCaldb caldb = m_caldb;
539  clear();
540  m_caldb = caldb;
541 
542  // Save response name
543  m_rspname = rspname;
544 
545  // Interpret response name as a FITS file name
546  GFilename filename(rspname);
547 
548  // If the filename does not exist the try getting the response from the
549  // calibration database
550  if (!filename.is_fits()) {
551 
552  // Get GCaldb response
553  filename = m_caldb.filename("","","IAQ","","",rspname);
554 
555  // If filename is empty then build filename from CALDB root path
556  // and response name
557  if (filename.is_empty()) {
558  filename = gammalib::filepath(m_caldb.rootdir(), rspname);
559  if (!filename.exists()) {
560  GFilename testname = filename + ".fits";
561  if (testname.exists()) {
562  filename = testname;
563  }
564  }
565  }
566 
567  } // endif: response name is not a FITS file
568 
569  // Open FITS file
570  GFits fits(filename);
571 
572  // Get IAQ image
573  const GFitsImage& iaq = *fits.image(0);
574 
575  // Read IAQ
576  read(iaq);
577 
578  // Close IAQ FITS file
579  fits.close();
580 
581  // Return
582  return;
583 }
584 
585 
586 /***********************************************************************//**
587  * @brief Read COMPTEL response from FITS image.
588  *
589  * @param[in] image FITS image.
590  *
591  * Read the COMPTEL response from IAQ FITS file and convert the IAQ values
592  * into a probability per steradian.
593  ***************************************************************************/
594 void GCOMResponse::read(const GFitsImage& image)
595 {
596  // Continue only if there are two image axis
597  if (image.naxis() == 2) {
598 
599  // Store IAQ dimensions
600  m_phigeo_bins = image.naxes(0);
601  m_phibar_bins = image.naxes(1);
602 
603  // Store IAQ axes definitions
604  m_phigeo_ref_value = image.real("CRVAL1");
605  m_phigeo_ref_pixel = image.real("CRPIX1");
606  m_phigeo_bin_size = image.real("CDELT1");
607  m_phibar_ref_value = image.real("CRVAL2");
608  m_phibar_ref_pixel = image.real("CRPIX2");
609  m_phibar_bin_size = image.real("CDELT2");
610 
611  // Get axes minima (values of first bin)
614 
615  // Compute IAQ size. Continue only if size is positive
616  int size = m_phigeo_bins * m_phibar_bins;
617  if (size > 0) {
618 
619  // Allocate memory for IAQ
620  m_iaq.assign(size, 0.0);
621 
622  // Copy over IAQ values
623  for (int i = 0; i < size; ++i) {
624  m_iaq[i] = image.pixel(i);
625  }
626 
627  } // endif: size was positive
628 
629  // Convert IAQ matrix from probability per Phigeo bin into a
630  // probability per steradian
631  double omega0 = gammalib::fourpi *
633  for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
634  double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
635  double omega = omega0 * std::sin(phigeo * gammalib::deg2rad);
636  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
637  m_iaq[iphigeo+iphibar*m_phigeo_bins] /= omega;
638  }
639  }
640 
641  } // endif: image had two axes
642 
643  // Return
644  return;
645 }
646 
647 
648 /***********************************************************************//**
649  * @brief Write COMPTEL response into FITS image.
650  *
651  * @param[in] image FITS image.
652  *
653  * Writes the COMPTEL response into an IAQ FITS file.
654  ***************************************************************************/
656 {
657  // Continue only if response is not empty
658  if (m_phigeo_bins > 0 && m_phibar_bins > 0) {
659 
660  // Initialise image
662 
663  // Convert IAQ matrix from probability per steradian into
664  //probability per Phigeo bin
665  double omega0 = gammalib::fourpi *
667  for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
668  double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
669  double omega = omega0 * std::sin(phigeo * gammalib::deg2rad);
670  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
671  int i = iphigeo + iphibar*m_phigeo_bins;
672  image(iphigeo, iphibar) = m_iaq[i] * omega;
673  }
674  }
675 
676  // Set header keywords
677  image.card("CTYPE1", "Phigeo", "Geometrical scatter angle");
678  image.card("CRVAL1", m_phigeo_ref_value,
679  "[deg] Geometrical scatter angle for reference bin");
680  image.card("CDELT1", m_phigeo_bin_size,
681  "[deg] Geometrical scatter angle bin size");
682  image.card("CRPIX1", m_phigeo_ref_pixel,
683  "Reference bin of geometrical scatter angle");
684  image.card("CTYPE2", "Phibar", "Compton scatter angle");
685  image.card("CRVAL2", m_phibar_ref_value,
686  "[deg] Compton scatter angle for reference bin");
687  image.card("CDELT2", m_phibar_bin_size, "[deg] Compton scatter angle bin size");
688  image.card("CRPIX2", m_phibar_ref_pixel,
689  "Reference bin of Compton scatter angle");
690  image.card("BUNIT", "Probability per bin", "Unit of bins");
691  image.card("TELESCOP", "GRO", "Name of telescope");
692  image.card("INSTRUME", "COMPTEL", "Name of instrument");
693  image.card("ORIGIN", "GammaLib", "Origin of FITS file");
694  image.card("OBSERVER", "Unknown", "Observer that created FITS file");
695 
696  } // endif: response was not empty
697 
698  // Return
699  return;
700 }
701 
702 
703 /***********************************************************************//**
704  * @brief Print COMPTEL response information
705  *
706  * @param[in] chatter Chattiness.
707  * @return String containing COMPTEL response information.
708  ***************************************************************************/
709 std::string GCOMResponse::print(const GChatter& chatter) const
710 {
711  // Initialise result string
712  std::string result;
713 
714  // Continue only if chatter is not silent
715  if (chatter != SILENT) {
716 
717  // Append header
718  result.append("=== GCOMResponse ===");
719 
720  // Append response name
721  result.append("\n"+gammalib::parformat("Response name")+m_rspname);
722 
723  // EXPLICIT: Append detailed information
724  if (chatter >= EXPLICIT) {
725 
726  // Append information
727  result.append("\n"+gammalib::parformat("Number of Phigeo bins"));
728  result.append(gammalib::str(m_phigeo_bins));
729  result.append("\n"+gammalib::parformat("Number of Phibar bins"));
730  result.append(gammalib::str(m_phibar_bins));
731  result.append("\n"+gammalib::parformat("Phigeo reference value"));
732  result.append(gammalib::str(m_phigeo_ref_value)+" deg");
733  result.append("\n"+gammalib::parformat("Phigeo reference pixel"));
734  result.append(gammalib::str(m_phigeo_ref_pixel));
735  result.append("\n"+gammalib::parformat("Phigeo bin size"));
736  result.append(gammalib::str(m_phigeo_bin_size)+" deg");
737  result.append("\n"+gammalib::parformat("Phigeo first bin value"));
738  result.append(gammalib::str(m_phigeo_min)+" deg");
739  result.append("\n"+gammalib::parformat("Phibar reference value"));
740  result.append(gammalib::str(m_phibar_ref_value)+" deg");
741  result.append("\n"+gammalib::parformat("Phibar reference pixel"));
742  result.append(gammalib::str(m_phibar_ref_pixel));
743  result.append("\n"+gammalib::parformat("Phibar bin size"));
744  result.append(gammalib::str(m_phibar_bin_size)+" deg");
745  result.append("\n"+gammalib::parformat("Phibar first bin value"));
746  result.append(gammalib::str(m_phibar_min)+" deg");
747 
748  }
749 
750  // VERBOSE: Append calibration database
751  if (chatter == VERBOSE) {
752  result.append("\n"+m_caldb.print(chatter));
753  }
754 
755  } // endif: chatter was not silent
756 
757  // Return result
758  return result;
759 }
760 
761 
762 /*==========================================================================
763  = =
764  = Private methods =
765  = =
766  ==========================================================================*/
767 
768 /***********************************************************************//**
769  * @brief Initialise class members
770  ***************************************************************************/
772 {
773  // Initialise members
774  m_caldb.clear();
775  m_rspname.clear();
776  m_iaq.clear();
777  m_phigeo_bins = 0;
778  m_phibar_bins = 0;
779  m_phigeo_ref_value = 0.0;
780  m_phigeo_ref_pixel = 0.0;
781  m_phigeo_bin_size = 0.0;
782  m_phigeo_min = 0.0;
783  m_phibar_ref_value = 0.0;
784  m_phibar_ref_pixel = 0.0;
785  m_phibar_bin_size = 0.0;
786  m_phibar_min = 0.0;
787 
788  // Return
789  return;
790 }
791 
792 
793 /***********************************************************************//**
794  * @brief Copy class members
795  *
796  * @param[in] rsp COMPTEL response.
797  ***************************************************************************/
799 {
800  // Copy attributes
801  m_caldb = rsp.m_caldb;
802  m_rspname = rsp.m_rspname;
803  m_iaq = rsp.m_iaq;
814 
815  // Return
816  return;
817 }
818 
819 
820 /***********************************************************************//**
821  * @brief Delete class members
822  ***************************************************************************/
824 {
825  // Return
826  return;
827 }
virtual GCOMResponse * clone(void) const
Clone instance.
void phibar(const double &phibar)
Set event Compton scatter angle.
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:280
#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.
double m_phibar_min
Phigeo value of first bin (deg)
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
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.
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:260
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
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
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
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
Single precision FITS image class.
Photon class definition.
Abstract observation base class interface definition.
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:282
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:125
virtual std::string type(void) const =0
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:120
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:1226
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:640
const double fourpi
Definition: GMath.hpp:37
std::string filepath(const std::string &pathname, const std::string &filename)
Build file path from path name and file name.
Definition: GTools.cpp:386
void copy_members(const GCOMResponse &rsp)
Copy class members.
COMPTEL instrument direction 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
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
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
const GCOMDri & drm(const GSource &source) const
Return source model DRM cube.
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
Time class interface definition.
void clear(void)
Clear calibration database.
Definition: GCaldb.cpp:194
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:41
Mathematical function definitions.
const GSkyMap & drg(void) const
Return geometry factors.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
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:413