GammaLib  2.1.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-2022 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 "GVector.hpp"
36 #include "GMatrix.hpp"
37 #include "GIntegrals.hpp"
38 #include "GFits.hpp"
39 #include "GCaldb.hpp"
40 #include "GEvent.hpp"
41 #include "GPhoton.hpp"
42 #include "GSource.hpp"
43 #include "GEnergy.hpp"
44 #include "GTime.hpp"
45 #include "GObservation.hpp"
46 #include "GFitsImage.hpp"
47 #include "GFitsImageFloat.hpp"
48 #include "GFilename.hpp"
49 #include "GModelSky.hpp"
51 #include "GModelSpatialRadial.hpp"
53 #include "GModelSpectralConst.hpp"
54 #include "GCOMResponse.hpp"
55 #include "GCOMObservation.hpp"
56 #include "GCOMEventCube.hpp"
57 #include "GCOMEventBin.hpp"
58 #include "GCOMInstDir.hpp"
60 
61 
62 /* __ Method name definitions ____________________________________________ */
63 #define G_IRF "GCOMResponse::irf(GEvent&, GPhoton&, GObservation&)"
64 #define G_IRF_SPATIAL "GCOMResponse::irf_spatial(GEvent&, GSource&, "\
65  "GObservation&)"
66 #define G_NROI "GCOMResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
67  "GObservation&)"
68 #define G_EBOUNDS "GCOMResponse::ebounds(GEnergy&)"
69 #define G_IRF_PTSRC "GCOMResponse::irf_ptsrc(GModelSky&, GObservation&, "\
70  "GMatrix*)"
71 #define G_IRF_RADIAL "GCOMResponse::irf_radial(GModelSky&, GObservation&, "\
72  "GMatrix*)"
73 #define G_IRF_ELLIPTICAL "GCOMResponse::irf_elliptical(GModelSky&, "\
74  " GObservation&, GMatrix*)"
75 #define G_IRF_DIFFUSE "GCOMResponse::irf_diffuse(GModelSky&, GObservation&, "\
76  "GMatrix*)"
77 
78 /* __ Macros _____________________________________________________________ */
79 
80 /* __ Coding definitions _________________________________________________ */
81 
82 /* __ Debug definitions __________________________________________________ */
83 
84 /* __ Constants __________________________________________________________ */
85 
86 
87 /*==========================================================================
88  = =
89  = Constructors/destructors =
90  = =
91  ==========================================================================*/
92 
93 /***********************************************************************//**
94  * @brief Void constructor
95  *
96  * Creates an empty COMPTEL response.
97  ***************************************************************************/
99 {
100  // Initialise members
101  init_members();
102 
103  // Return
104  return;
105 }
106 
107 
108 /***********************************************************************//**
109  * @brief Copy constructor
110  *
111  * @param[in] rsp COM response.
112  **************************************************************************/
114 {
115  // Initialise members
116  init_members();
117 
118  // Copy members
119  copy_members(rsp);
120 
121  // Return
122  return;
123 }
124 
125 
126 /***********************************************************************//**
127  * @brief Response constructor
128  *
129  * @param[in] caldb Calibration database.
130  * @param[in] iaqname IAQ file name.
131  *
132  * Create COMPTEL response by loading an IAQ file from a calibration
133  * database.
134  ***************************************************************************/
136  const std::string& iaqname) : GResponse()
137 {
138  // Initialise members
139  init_members();
140 
141  // Set calibration database
142  this->caldb(caldb);
143 
144  // Load IRF
145  this->load(iaqname);
146 
147  // Return
148  return;
149 }
150 
151 
152 /***********************************************************************//**
153  * @brief Destructor
154  *
155  * Destroys instance of COMPTEL response object.
156  ***************************************************************************/
158 {
159  // Free members
160  free_members();
161 
162  // Return
163  return;
164 }
165 
166 
167 /*==========================================================================
168  = =
169  = Operators =
170  = =
171  ==========================================================================*/
172 
173 /***********************************************************************//**
174  * @brief Assignment operator
175  *
176  * @param[in] rsp COMPTEL response.
177  * @return COMPTEL response.
178  *
179  * Assigns COMPTEL response object to another COMPTEL response object. The
180  * assignment performs a deep copy of all information, hence the original
181  * object from which the assignment has been performed can be destroyed after
182  * this operation without any loss of information.
183  ***************************************************************************/
185 {
186  // Execute only if object is not identical
187  if (this != &rsp) {
188 
189  // Copy base class members
190  this->GResponse::operator=(rsp);
191 
192  // Free members
193  free_members();
194 
195  // Initialise members
196  init_members();
197 
198  // Copy members
199  copy_members(rsp);
200 
201  } // endif: object was not identical
202 
203  // Return this object
204  return *this;
205 }
206 
207 
208 /*==========================================================================
209  = =
210  = Public methods =
211  = =
212  ==========================================================================*/
213 
214 /***********************************************************************//**
215  * @brief Clear instance
216  *
217  * Clears COMPTEL response object by resetting all members to an initial
218  * state. Any information that was present in the object before will be lost.
219  ***************************************************************************/
221 {
222  // Free class members (base and derived classes, derived class first)
223  free_members();
224  this->GResponse::free_members();
225 
226  // Initialise members
227  this->GResponse::init_members();
228  init_members();
229 
230  // Return
231  return;
232 }
233 
234 
235 /***********************************************************************//**
236  * @brief Clone instance
237  *
238  * @return Pointer to deep copy of COMPTEL response.
239  ***************************************************************************/
241 {
242  return new GCOMResponse(*this);
243 }
244 
245 
246 /***********************************************************************//**
247  * @brief Return value of instrument response function
248  *
249  * @param[in] event Observed event.
250  * @param[in] photon Incident photon.
251  * @param[in] obs Observation.
252  * @return Instrument response function (\f$cm^2 sr^{-1}\f$)
253  *
254  * @exception GException::invalid_argument
255  * Observation is not a COMPTEL observation.
256  * Event is not a COMPTEL event bin.
257  * @exception GException::invalid_value
258  * Response not initialised with a valid IAQ
259  *
260  * Returns the instrument response function for a given observed photon
261  * direction as function of the assumed true photon direction. The result
262  * is given by
263  *
264  * \f[
265  * {\tt IRF} = \frac{{\tt IAQ} \times {\tt DRG} \times {\tt DRX}}
266  * {T \times {\tt TOFCOR}} \times {\tt PHASECOR}
267  * \f]
268  *
269  * where
270  * - \f${\tt IRF}\f$ is the instrument response function (\f$cm^2 sr^{-1}\f$),
271  * - \f${\tt IAQ}\f$ is the COMPTEL response matrix (\f$sr^{-1}\f$),
272  * - \f${\tt DRG}\f$ is the geometry factor (probability),
273  * - \f${\tt DRX}\f$ is the exposure (\f$cm^2 s\f$),
274  * - \f$T\f$ is the ontime (\f$s\f$), and
275  * - \f${\tt TOFCOR}\f$ is a correction that accounts for the Time of Flight
276  * selection window.
277  * - \f${\tt PHASECOR}\f$ is a correction that accounts for pulsar phase
278  * selection.
279  *
280  * The observed photon direction is spanned by the three values \f$\Chi\f$,
281  * \f$\Psi\f$, and \f$\bar{\varphi})\f$. \f$\Chi\f$ and \f$\Psi\f$ is the
282  * scatter direction of the event, given in sky coordinates.
283  * \f$\bar{\varphi}\f$ is the Compton scatter angle, computed from the
284  * energy deposits in the two detector planes.
285  ***************************************************************************/
286 double GCOMResponse::irf(const GEvent& event,
287  const GPhoton& photon,
288  const GObservation& obs) const
289 {
290  // Extract COMPTEL observation
291  const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
292  if (observation == NULL) {
293  std::string cls = std::string(typeid(&obs).name());
294  std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
295  "observations. Please specify a COMPTEL observation "
296  "as argument.";
298  }
299 
300  // Extract COMPTEL event cube
301  const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(observation->events());
302  if (cube == NULL) {
303  std::string cls = std::string(typeid(&cube).name());
304  std::string msg = "Event cube of type \""+cls+"\" is not a COMPTEL "
305  "event cube. Please specify a COMPTEL event cube "
306  "as argument.";
308  }
309 
310  // Extract COMPTEL event bin
311  const GCOMEventBin* bin = dynamic_cast<const GCOMEventBin*>(&event);
312  if (bin == NULL) {
313  std::string cls = std::string(typeid(&event).name());
314  std::string msg = "Event of type \""+cls+"\" is not a COMPTEL event. "
315  "Please specify a COMPTEL event as argument.";
317  }
318 
319  // Throw an exception if COMPTEL response is not set or if
320  if (m_iaq.empty()) {
321  std::string msg = "COMPTEL response is empty. Please initialise the "
322  "response with an \"IAQ\".";
323  throw GException::invalid_value(G_IRF, msg);
324  }
325  else if (m_phigeo_bin_size == 0.0) {
326  std::string msg = "COMPTEL response has a zero Phigeo bin size. "
327  "Please initialise the response with a valid "
328  "\"IAQ\".";
329  throw GException::invalid_value(G_IRF, msg);
330  }
331 
332  // Extract event parameters
333  const GCOMInstDir& obsDir = bin->dir();
334 
335  // Extract photon parameters
336  const GSkyDir& srcDir = photon.dir();
337  const GTime& srcTime = photon.time();
338 
339  // Compute angle between true photon arrival direction and scatter
340  // direction (Chi,Psi)
341  double phigeo = srcDir.dist_deg(obsDir.dir());
342 
343  // Compute scatter angle index
344  int iphibar = int(obsDir.phibar() / m_phibar_bin_size);
345 
346  // Initialise IRF
347  double iaq = 0.0;
348  double irf = 0.0;
349 
350  // Extract IAQ value by linear inter/extrapolation in Phigeo
351  if (iphibar < m_phibar_bins) {
352  double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
353  int iphigeo = int(phirat); // index into which Phigeo falls
354  double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
355  if (iphigeo < m_phigeo_bins) {
356  int i = iphibar * m_phigeo_bins + iphigeo;
357  if (eps < 0.0) { // interpolate towards left
358  if (iphigeo > 0) {
359  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
360  }
361  else {
362  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
363  }
364  }
365  else { // interpolate towards right
366  if (iphigeo < m_phigeo_bins-1) {
367  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
368  }
369  else {
370  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
371  }
372  }
373  }
374  }
375 
376  // Continue only if IAQ is positive
377  if (iaq > 0.0) {
378 
379  // Get DRG value (unit: probability)
380  double drg = observation->drg().map()(obsDir.dir(), iphibar);
381 
382  // Get DRX value (unit: cm^2 sec)
383  double drx = observation->drx().map()(srcDir);
384 
385  // Get ontime (unit: s)
386  double ontime = observation->ontime();
387 
388  // Get ToF correction
389  double tofcor = cube->dre().tof_correction();
390 
391  // Get pulsar phase correction
392  double phasecor = cube->dre().phase_correction();
393 
394  // Compute IRF value
395  irf = iaq * drg * drx / (ontime * tofcor) * phasecor;
396 
397  // Apply deadtime correction
398  irf *= obs.deadc(srcTime);
399 
400  // Make sure that IRF is positive
401  if (irf < 0.0) {
402  irf = 0.0;
403  }
404 
405  // Compile option: Check for NaN/Inf
406  #if defined(G_NAN_CHECK)
408  std::cout << "*** ERROR: GCOMResponse::irf:";
409  std::cout << " NaN/Inf encountered";
410  std::cout << " (irf=" << irf;
411  std::cout << ", iaq=" << iaq;
412  std::cout << ", drg=" << drg;
413  std::cout << ", drx=" << drx;
414  std::cout << ")";
415  std::cout << std::endl;
416  }
417  #endif
418 
419  } // endif: IAQ was positive
420 
421  // Return IRF value
422  return irf;
423 }
424 
425 
426 /***********************************************************************//**
427  * @brief Return integral of event probability for a given sky model over ROI
428  *
429  * @param[in] model Sky model.
430  * @param[in] obsEng Observed photon energy.
431  * @param[in] obsTime Observed photon arrival time.
432  * @param[in] obs Observation.
433  * @return 0.0
434  *
435  * @exception GException::feature_not_implemented
436  * Method is not implemented.
437  ***************************************************************************/
438 double GCOMResponse::nroi(const GModelSky& model,
439  const GEnergy& obsEng,
440  const GTime& obsTime,
441  const GObservation& obs) const
442 {
443  // Method is not implemented
444  std::string msg = "Spatial integration of sky model over the data space "
445  "is not implemented.";
447 
448  // Return Npred
449  return (0.0);
450 }
451 
452 
453 /***********************************************************************//**
454  * @brief Return true energy boundaries for a specific observed energy
455  *
456  * @param[in] obsEnergy Observed Energy.
457  * @return True energy boundaries for given observed energy.
458  *
459  * @exception GException::feature_not_implemented
460  * Method is not implemented.
461  ***************************************************************************/
462 GEbounds GCOMResponse::ebounds(const GEnergy& obsEnergy) const
463 {
464  // Initialise an empty boundary object
466 
467  // Throw an exception
468  std::string msg = "Energy dispersion not implemented.";
470 
471  // Return energy boundaries
472  return ebounds;
473 }
474 
475 
476 /***********************************************************************//**
477  * @brief Load COMPTEL response.
478  *
479  * @param[in] rspname COMPTEL response name.
480  *
481  * Loads the COMPTEL response with specified name @p rspname.
482  *
483  * The method first attempts to interpret @p rspname as a file name and to
484  * load the corresponding response.
485  *
486  * If @p rspname is not a FITS file the method searches for an appropriate
487  * response in the calibration database. If no appropriate response is found,
488  * the method takes the CALDB root path and response name to build the full
489  * path to the response file, and tries to load the response from these
490  * paths.
491  *
492  * If also this fails an exception is thrown.
493  ***************************************************************************/
494 void GCOMResponse::load(const std::string& rspname)
495 {
496  // Clear instance but conserve calibration database
497  GCaldb caldb = m_caldb;
498  clear();
499  m_caldb = caldb;
500 
501  // Save response name
502  m_rspname = rspname;
503 
504  // Interpret response name as a FITS file name
505  GFilename filename(rspname);
506 
507  // If the filename does not exist the try getting the response from the
508  // calibration database
509  if (!filename.is_fits()) {
510 
511  // Get GCaldb response
512  filename = m_caldb.filename("","","IAQ","","",rspname);
513 
514  // If filename is empty then build filename from CALDB root path
515  // and response name
516  if (filename.is_empty()) {
517  filename = gammalib::filepath(m_caldb.rootdir(), rspname);
518  if (!filename.exists()) {
519  GFilename testname = filename + ".fits";
520  if (testname.exists()) {
521  filename = testname;
522  }
523  }
524  }
525 
526  } // endif: response name is not a FITS file
527 
528  // Open FITS file
529  GFits fits(filename);
530 
531  // Get IAQ image
532  const GFitsImage& iaq = *fits.image(0);
533 
534  // Read IAQ
535  read(iaq);
536 
537  // Close IAQ FITS file
538  fits.close();
539 
540  // Return
541  return;
542 }
543 
544 
545 /***********************************************************************//**
546  * @brief Read COMPTEL response from FITS image.
547  *
548  * @param[in] image FITS image.
549  *
550  * Read the COMPTEL response from IAQ FITS file and convert the IAQ values
551  * into a probability per steradian.
552  *
553  * The IAQ values are divided by the solid angle of a Phigeo bin which is
554  * given by
555  *
556  * \f{eqnarray*}{
557  * \Omega & = & 2 \pi \left[
558  * \left(
559  * 1 - \cos \left( \varphi_{\rm geo} +
560  * \frac{1}{2} \Delta \varphi_{\rm geo} \right)
561  * \right) -
562  * \left(
563  * 1 - \cos \left( \varphi_{\rm geo} -
564  * \frac{1}{2} \Delta \varphi_{\rm geo} \right)
565  * \right) \right] \\
566  * &=& 2 \pi \left[
567  * \cos \left( \varphi_{\rm geo} -
568  * \frac{1}{2} \Delta \varphi_{\rm geo} \right) -
569  * \cos \left( \varphi_{\rm geo} +
570  * \frac{1}{2} \Delta \varphi_{\rm geo} \right)
571  * \right] \\
572  * &=& 4 \pi \sin \left( \varphi_{\rm geo} \right)
573  * \sin \left( \frac{1}{2} \Delta \varphi_{\rm geo} \right)
574  * \f}
575  ***************************************************************************/
576 void GCOMResponse::read(const GFitsImage& image)
577 {
578  // Continue only if there are two image axis
579  if (image.naxis() == 2) {
580 
581  // Store IAQ dimensions
582  m_phigeo_bins = image.naxes(0);
583  m_phibar_bins = image.naxes(1);
584 
585  // Store IAQ axes definitions
586  m_phigeo_ref_value = image.real("CRVAL1");
587  m_phigeo_ref_pixel = image.real("CRPIX1");
588  m_phigeo_bin_size = image.real("CDELT1");
589  m_phibar_ref_value = image.real("CRVAL2");
590  m_phibar_ref_pixel = image.real("CRPIX2");
591  m_phibar_bin_size = image.real("CDELT2");
592 
593  // Get axes minima (values of first bin)
598 
599  // Compute IAQ size. Continue only if size is positive
600  int size = m_phigeo_bins * m_phibar_bins;
601  if (size > 0) {
602 
603  // Allocate memory for IAQ
604  m_iaq.assign(size, 0.0);
605 
606  // Copy over IAQ values
607  for (int i = 0; i < size; ++i) {
608  m_iaq[i] = image.pixel(i);
609  }
610 
611  } // endif: size was positive
612 
613  // Precompute variable for conversion
614  double phigeo_min = m_phigeo_min * gammalib::deg2rad;
615  double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
616  double omega0 = gammalib::fourpi * std::sin(0.5 * phigeo_bin);
617 
618  // Convert IAQ matrix from probability per Phigeo bin into a
619  // probability per steradian
620  for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
621  double phigeo = iphigeo * phigeo_bin + phigeo_min;
622  double omega = omega0 * std::sin(phigeo);
623  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
624  m_iaq[iphigeo+iphibar*m_phigeo_bins] /= omega;
625  }
626  }
627 
628  } // endif: image had two axes
629 
630  // Return
631  return;
632 }
633 
634 
635 /***********************************************************************//**
636  * @brief Write COMPTEL response into FITS image.
637  *
638  * @param[in] image FITS image.
639  *
640  * Writes the COMPTEL response into an IAQ FITS file.
641  ***************************************************************************/
643 {
644  // Continue only if response is not empty
645  if (m_phigeo_bins > 0 && m_phibar_bins > 0) {
646 
647  // Initialise image
649 
650  // Convert IAQ matrix from probability per steradian into
651  //probability per Phigeo bin
652  double omega0 = gammalib::fourpi *
654  for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
655  double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
656  double omega = omega0 * std::sin(phigeo * gammalib::deg2rad);
657  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
658  int i = iphigeo + iphibar*m_phigeo_bins;
659  image(iphigeo, iphibar) = m_iaq[i] * omega;
660  }
661  }
662 
663  // Set header keywords
664  image.card("CTYPE1", "Phigeo", "Geometrical scatter angle");
665  image.card("CRVAL1", m_phigeo_ref_value,
666  "[deg] Geometrical scatter angle for reference bin");
667  image.card("CDELT1", m_phigeo_bin_size,
668  "[deg] Geometrical scatter angle bin size");
669  image.card("CRPIX1", m_phigeo_ref_pixel,
670  "Reference bin of geometrical scatter angle");
671  image.card("CTYPE2", "Phibar", "Compton scatter angle");
672  image.card("CRVAL2", m_phibar_ref_value,
673  "[deg] Compton scatter angle for reference bin");
674  image.card("CDELT2", m_phibar_bin_size, "[deg] Compton scatter angle bin size");
675  image.card("CRPIX2", m_phibar_ref_pixel,
676  "Reference bin of Compton scatter angle");
677  image.card("BUNIT", "Probability per bin", "Unit of bins");
678  image.card("TELESCOP", "GRO", "Name of telescope");
679  image.card("INSTRUME", "COMPTEL", "Name of instrument");
680  image.card("ORIGIN", "GammaLib", "Origin of FITS file");
681  image.card("OBSERVER", "Unknown", "Observer that created FITS file");
682 
683  } // endif: response was not empty
684 
685  // Return
686  return;
687 }
688 
689 
690 /***********************************************************************//**
691  * @brief Load response cache.
692  *
693  * @param[in] filename Response cache filename.
694  *
695  * Loads response cache from FITS file.
696  ***************************************************************************/
697 void GCOMResponse::load_cache(const GFilename& filename)
698 {
699  // Load response vector cache
700  m_irf_vector_cache.load(filename);
701 
702  // Return
703  return;
704 }
705 
706 
707 /***********************************************************************//**
708  * @brief Save response cache.
709  *
710  * @param[in] filename Response cache filename.
711  *
712  * Saves response cache from FITS file.
713  ***************************************************************************/
714 void GCOMResponse::save_cache(const GFilename& filename) const
715 {
716  // Save response vector cache
717  m_irf_vector_cache.save(filename);
718 
719  // Return
720  return;
721 }
722 
723 
724 /***********************************************************************//**
725  * @brief Print COMPTEL response information
726  *
727  * @param[in] chatter Chattiness.
728  * @return String containing COMPTEL response information.
729  ***************************************************************************/
730 std::string GCOMResponse::print(const GChatter& chatter) const
731 {
732  // Initialise result string
733  std::string result;
734 
735  // Continue only if chatter is not silent
736  if (chatter != SILENT) {
737 
738  // Append header
739  result.append("=== GCOMResponse ===");
740 
741  // Append response name
742  result.append("\n"+gammalib::parformat("Response name")+m_rspname);
743 
744  // EXPLICIT: Append detailed information
745  if (chatter >= EXPLICIT) {
746 
747  // Append information
748  result.append("\n"+gammalib::parformat("Number of Phigeo bins"));
749  result.append(gammalib::str(m_phigeo_bins));
750  result.append("\n"+gammalib::parformat("Number of Phibar bins"));
751  result.append(gammalib::str(m_phibar_bins));
752  result.append("\n"+gammalib::parformat("Phigeo reference value"));
753  result.append(gammalib::str(m_phigeo_ref_value)+" deg");
754  result.append("\n"+gammalib::parformat("Phigeo reference pixel"));
755  result.append(gammalib::str(m_phigeo_ref_pixel));
756  result.append("\n"+gammalib::parformat("Phigeo bin size"));
757  result.append(gammalib::str(m_phigeo_bin_size)+" deg");
758  result.append("\n"+gammalib::parformat("Phigeo first bin value"));
759  result.append(gammalib::str(m_phigeo_min)+" deg");
760  result.append("\n"+gammalib::parformat("Phibar reference value"));
761  result.append(gammalib::str(m_phibar_ref_value)+" deg");
762  result.append("\n"+gammalib::parformat("Phibar reference pixel"));
763  result.append(gammalib::str(m_phibar_ref_pixel));
764  result.append("\n"+gammalib::parformat("Phibar bin size"));
765  result.append(gammalib::str(m_phibar_bin_size)+" deg");
766  result.append("\n"+gammalib::parformat("Phibar first bin value"));
767  result.append(gammalib::str(m_phibar_min)+" deg");
768 
769  }
770 
771  // VERBOSE: Append calibration database
772  if (chatter == VERBOSE) {
773  result.append("\n"+m_caldb.print(chatter));
774  }
775 
776  } // endif: chatter was not silent
777 
778  // Return result
779  return result;
780 }
781 
782 
783 /*==========================================================================
784  = =
785  = Private methods =
786  = =
787  ==========================================================================*/
788 
789 /***********************************************************************//**
790  * @brief Initialise class members
791  ***************************************************************************/
793 {
794  // Initialise members
795  m_caldb.clear();
796  m_rspname.clear();
797  m_iaq.clear();
798  m_phigeo_bins = 0;
799  m_phibar_bins = 0;
800  m_phigeo_ref_value = 0.0;
801  m_phigeo_ref_pixel = 0.0;
802  m_phigeo_bin_size = 0.0;
803  m_phigeo_min = 0.0;
804  m_phibar_ref_value = 0.0;
805  m_phibar_ref_pixel = 0.0;
806  m_phibar_bin_size = 0.0;
807  m_phibar_min = 0.0;
808 
809  // Return
810  return;
811 }
812 
813 
814 /***********************************************************************//**
815  * @brief Copy class members
816  *
817  * @param[in] rsp COMPTEL response.
818  ***************************************************************************/
820 {
821  // Copy attributes
822  m_caldb = rsp.m_caldb;
823  m_rspname = rsp.m_rspname;
824  m_iaq = rsp.m_iaq;
835 
836  // Return
837  return;
838 }
839 
840 
841 /***********************************************************************//**
842  * @brief Delete class members
843  ***************************************************************************/
845 {
846  // Return
847  return;
848 }
849 
850 
851 /***********************************************************************//**
852  * @brief Return instrument response to point source
853  *
854  * @param[in] model Sky model.
855  * @param[in] obs Observation.
856  * @param[out] gradients Gradients matrix.
857  * @return Instrument response to point source for all events in
858  * observation (\f$cm^2\f$).
859  *
860  * @exception GException::invalid_argument
861  * Observation is not a COMPTEL observation.
862  * Event is not a COMPTEL event bin.
863  * @exception GException::invalid_value
864  * Response not initialised with a valid IAQ
865  * Incompatible IAQ encountered
866  *
867  * Returns the instrument response to a point source for all events in the
868  * observations.
869  *
870  * @p gradients is an optional matrix where the number of rows corresponds
871  * to the number of events in the observation and the number of columns
872  * corresponds to the number of spatial model parameters. Since for point
873  * sources no gradients are computed, the method does not alter the
874  * content of @p gradients.
875  ***************************************************************************/
877  const GObservation& obs,
878  GMatrix* gradients) const
879 {
880  // Extract COMPTEL observation
881  const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
882  if (obs_ptr == NULL) {
883  std::string cls = std::string(typeid(&obs).name());
884  std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
885  "observations. Please specify a COMPTEL observation "
886  "as argument.";
888  }
889 
890  // Extract COMPTEL event cube
891  const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(obs_ptr->events());
892  if (cube == NULL) {
893  std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
894  "not contain a COMPTEL event cube. Please specify "
895  "a COMPTEL observation containing and event cube.";
897  }
898 
899  // Throw an exception if COMPTEL response is not set or if
900  if (m_iaq.empty()) {
901  std::string msg = "COMPTEL response is empty. Please initialise the "
902  "response with an \"IAQ\".";
904  }
905  else if (m_phigeo_bin_size == 0.0) {
906  std::string msg = "COMPTEL response has a zero Phigeo bin size. "
907  "Please initialise the response with a valid "
908  "\"IAQ\".";
910  }
911 
912  // Get number of Chi/Psi pixels, Phibar layers and event bins
913  int npix = cube->naxis(0) * cube->naxis(1);
914  int nphibar = cube->naxis(2);
915  int nevents = cube->size();
916 
917  // Throw an exception if the number of Phibar bins does not match the
918  // response
919  if (nphibar != m_phibar_bins) {
920  std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
921  "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
922  "bins. Please specify a compatible IAQ.";
924  }
925 
926  // Initialise result
927  GVector irfs(nevents);
928 
929  // Get point source direction
930  GSkyDir srcDir =
931  static_cast<const GModelSpatialPointSource*>(model.spatial())->dir();
932 
933  // Get IAQ normalisation (cm2): DRX (cm2 s) * DEADC / ONTIME (s)
934  double iaq_norm = obs_ptr->drx().map()(srcDir) * obs_ptr->deadc() /
935  (obs_ptr->ontime() * cube->dre().tof_correction()) *
936  cube->dre().phase_correction();
937 
938  // Get pointer to DRG pixels
939  const double* drg = obs_ptr->drg().map().pixels();
940 
941  // Loop over Chi and Psi
942  for (int ipix = 0; ipix < npix; ++ipix) {
943 
944  // Get pointer to event bin
945  const GCOMEventBin* bin = (*cube)[ipix];
946 
947  // Get reference to instrument direction
948  const GCOMInstDir& obsDir = bin->dir();
949 
950  // Get reference to Phigeo sky direction
951  const GSkyDir& phigeoDir = obsDir.dir();
952 
953  // Compute angle between true photon arrival direction and scatter
954  // direction (Chi,Psi)
955  double phigeo = srcDir.dist_deg(phigeoDir);
956 
957  // Compute interpolation factors
958  double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
959  int iphigeo = int(phirat); // index into which Phigeo falls
960  double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
961 
962  // If Phigeo is in range then compute the IRF value for all
963  // Phibar layers
964  if (iphigeo < m_phigeo_bins) {
965 
966  // Loop over Phibar
967  for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
968 
969  // Get IAQ index
970  int i = iphibar * m_phigeo_bins + iphigeo;
971 
972  // Initialise IAQ
973  double iaq = 0.0;
974 
975  // Compute IAQ
976  if (eps < 0.0) { // interpolate towards left
977  if (iphigeo > 0) {
978  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
979  }
980  else {
981  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
982  }
983  }
984  else { // interpolate towards right
985  if (iphigeo < m_phigeo_bins-1) {
986  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
987  }
988  else {
989  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
990  }
991  }
992 
993  // Continue only if IAQ is positive
994  if (iaq > 0.0) {
995 
996  // Get DRI index
997  int idri = ipix + iphibar * npix;
998 
999  // Compute IRF value
1000  double irf = iaq * drg[idri] * iaq_norm;
1001 
1002  // Make sure that IRF is positive
1003  if (irf < 0.0) {
1004  irf = 0.0;
1005  }
1006 
1007  // Store IRF value
1008  irfs[idri] = irf;
1009 
1010  } // endif: IAQ was positive
1011 
1012  } // endfor: looped over Phibar
1013 
1014  } // endif: Phigeo was in valid range
1015 
1016  } // endfor: looped over Chi and Psi
1017 
1018  // Return IRF vector
1019  return irfs;
1020 }
1021 
1022 
1023 /***********************************************************************//**
1024  * @brief Return instrument response to radial source
1025  *
1026  * @param[in] model Sky model.
1027  * @param[in] obs Observation.
1028  * @param[out] gradients Gradients matrix.
1029  * @return Instrument response to radial source for all events in
1030  * observation (\f$cm^2\f$).
1031  *
1032  * @exception GException::invalid_argument
1033  * Observation is not a COMPTEL observation.
1034  * Event is not a COMPTEL event bin.
1035  * @exception GException::invalid_value
1036  * Response not initialised with a valid IAQ
1037  * Incompatible IAQ encountered
1038  *
1039  * Returns the instrument response to a radial source for all events in
1040  * the observations. See GCOMResponse::irf_extended for more information.
1041  ***************************************************************************/
1043  const GObservation& obs,
1044  GMatrix* gradients) const
1045 {
1046  // Set constants
1047  static const int iter_rho = 5;
1048  static const int iter_omega = 5;
1049 
1050  // Extract COMPTEL observation
1051  const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
1052  if (obs_ptr == NULL) {
1053  std::string cls = std::string(typeid(&obs).name());
1054  std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
1055  "observations. Please specify a COMPTEL observation "
1056  "as argument.";
1058  }
1059 
1060  // Extract COMPTEL event cube
1061  const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(obs_ptr->events());
1062  if (cube == NULL) {
1063  std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
1064  "not contain a COMPTEL event cube. Please specify "
1065  "a COMPTEL observation containing and event cube.";
1067  }
1068 
1069  // Throw an exception if COMPTEL response is not set or if
1070  if (m_iaq.empty()) {
1071  std::string msg = "COMPTEL response is empty. Please initialise the "
1072  "response with an \"IAQ\".";
1074  }
1075  else if (m_phigeo_bin_size == 0.0) {
1076  std::string msg = "COMPTEL response has a zero Phigeo bin size. "
1077  "Please initialise the response with a valid "
1078  "\"IAQ\".";
1080  }
1081 
1082  // Get pointer to radial model
1083  const GModelSpatialRadial* radial =
1084  static_cast<const GModelSpatialRadial*>(model.spatial());
1085 
1086  // Get references to model direction and maximum model radius. Reduce
1087  // the maximum model radius by a small amount to avoid rounding
1088  // problems at the boundary of a sharp edged model.
1089  const GSkyDir& model_dir = radial->dir();
1090  const double& theta_max = radial->theta_max() - 1.0e-12;
1091 
1092  // Get number of Chi/Psi pixels, Phibar layers and event bins
1093  int npix = cube->naxis(0) * cube->naxis(1);
1094  int nphibar = cube->naxis(2);
1095  int nevents = cube->size();
1096 
1097  // Throw an exception if the number of Phibar bins does not match the
1098  // response
1099  if (nphibar != m_phibar_bins) {
1100  std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
1101  "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
1102  "bins. Please specify a compatible IAQ.";
1104  }
1105 
1106  // Initialise result
1107  GVector irfs(nevents);
1108 
1109  // Initialise some variables
1110  double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
1111  double phigeo_min = m_phigeo_min * gammalib::deg2rad - 0.5 * phigeo_bin;
1112  double phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins;
1113  const GSkyMap& drx = obs_ptr->drx().map();
1114  const double* drg = obs_ptr->drg().map().pixels();
1115 
1116  // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
1117  double iaq_norm = obs_ptr->deadc() /
1118  (obs_ptr->ontime() * cube->dre().tof_correction()) *
1119  cube->dre().phase_correction();
1120 
1121  // Setup rotation matrix for conversion towards sky coordinates
1122  GMatrix ry;
1123  GMatrix rz;
1124  ry.eulery(model_dir.dec_deg() - 90.0);
1125  rz.eulerz(-model_dir.ra_deg());
1126  GMatrix rot = (ry * rz).transpose();
1127 
1128  // Loop over Chi and Psi
1129  for (int ipix = 0; ipix < npix; ++ipix) {
1130 
1131  // Get pointer to event bin
1132  const GCOMEventBin* bin = (*cube)[ipix];
1133 
1134  // Get reference to instrument direction
1135  const GCOMInstDir& obsDir = bin->dir();
1136 
1137  // Get reference to Phigeo sky direction
1138  const GSkyDir& phigeoDir = obsDir.dir();
1139 
1140  // Compute angle between model centre and Phigeo sky direction (radians)
1141  // and position angle of model with respect to Phigeo sky direction (radians)
1142  double zeta = phigeoDir.dist(model_dir);
1143 
1144  // Set radial model zenith angle range
1145  double rho_min = (zeta > phigeo_max) ? zeta - phigeo_max : 0.0;
1146  double rho_max = zeta + phigeo_max;
1147  if (rho_min > theta_max) {
1148  rho_min = theta_max;
1149  }
1150  if (rho_max > theta_max) {
1151  rho_max = theta_max;
1152  }
1153 
1154  // Perform model zenith angle integration if interval is valid
1155  if (rho_max > rho_min) {
1156 
1157  // Initialise IRF vector
1158  GVector irf(nphibar);
1159 
1160  // Setup integration kernel
1161  com_radial_kerns_rho integrands(m_iaq,
1162  *radial,
1163  irf,
1164  bin,
1165  rot,
1166  drx,
1167  phigeo_bin,
1168  m_phigeo_bins,
1169  nphibar,
1170  iter_omega);
1171 
1172  // Setup integrator
1173  GIntegrals integral(&integrands);
1174  integral.fixed_iter(iter_rho);
1175 
1176  // Integrate over Phigeo
1177  irf = integral.romberg(rho_min, rho_max, iter_rho);
1178 
1179  // Add IRF to result
1180  for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1181  int idri = ipix + iphibar * npix;
1182  irfs[idri] = irf[iphibar] * drg[idri] * iaq_norm;
1183  }
1184 
1185  } // endif: Phigeo angle interval was valid
1186 
1187  } // endfor: looped over Chi and Psi pixels
1188 
1189  // Return IRF vector
1190  return irfs;
1191 }
1192 
1193 
1194 /***********************************************************************//**
1195  * @brief Return instrument response to elliptical source
1196  *
1197  * @param[in] model Sky model.
1198  * @param[in] obs Observation.
1199  * @param[out] gradients Gradients matrix.
1200  * @return Instrument response to elliptical source for all events in
1201  * observation (\f$cm^2\f$).
1202  *
1203  * @exception GException::invalid_argument
1204  * Observation is not a COMPTEL observation.
1205  * Event is not a COMPTEL event bin.
1206  * @exception GException::invalid_value
1207  * Response not initialised with a valid IAQ
1208  * Incompatible IAQ encountered
1209  *
1210  * Returns the instrument response to an elliptical source for all events in
1211  * the observations. See GCOMResponse::irf_extended for more information.
1212  ***************************************************************************/
1214  const GObservation& obs,
1215  GMatrix* gradients) const
1216 {
1217  // Set constants
1218  static const int iter_rho = 5;
1219  static const int iter_omega = 5;
1220 
1221  // Extract COMPTEL observation
1222  const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
1223  if (obs_ptr == NULL) {
1224  std::string cls = std::string(typeid(&obs).name());
1225  std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
1226  "observations. Please specify a COMPTEL observation "
1227  "as argument.";
1229  }
1230 
1231  // Extract COMPTEL event cube
1232  const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(obs_ptr->events());
1233  if (cube == NULL) {
1234  std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
1235  "not contain a COMPTEL event cube. Please specify "
1236  "a COMPTEL observation containing and event cube.";
1238  }
1239 
1240  // Throw an exception if COMPTEL response is not set or if
1241  if (m_iaq.empty()) {
1242  std::string msg = "COMPTEL response is empty. Please initialise the "
1243  "response with an \"IAQ\".";
1245  }
1246  else if (m_phigeo_bin_size == 0.0) {
1247  std::string msg = "COMPTEL response has a zero Phigeo bin size. "
1248  "Please initialise the response with a valid "
1249  "\"IAQ\".";
1251  }
1252 
1253  // Get pointer to elliptical model
1254  const GModelSpatialElliptical* elliptical =
1255  static_cast<const GModelSpatialElliptical*>(model.spatial());
1256 
1257  // Get references to model direction and maximum model radius
1258  const GSkyDir& model_dir = elliptical->dir();
1259  const double& theta_max = elliptical->theta_max();
1260 
1261  // Get number of Chi/Psi pixels, Phibar layers and event bins
1262  int npix = cube->naxis(0) * cube->naxis(1);
1263  int nphibar = cube->naxis(2);
1264  int nevents = cube->size();
1265 
1266  // Throw an exception if the number of Phibar bins does not match the
1267  // response
1268  if (nphibar != m_phibar_bins) {
1269  std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
1270  "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
1271  "bins. Please specify a compatible IAQ.";
1273  }
1274 
1275  // Initialise result
1276  GVector irfs(nevents);
1277 
1278  // Initialise some variables
1279  double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
1280  double phigeo_min = m_phigeo_min * gammalib::deg2rad - 0.5 * phigeo_bin;
1281  double phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins;
1282  const GSkyMap& drx = obs_ptr->drx().map();
1283  const double* drg = obs_ptr->drg().map().pixels();
1284 
1285  // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
1286  double iaq_norm = obs_ptr->deadc() /
1287  (obs_ptr->ontime() * cube->dre().tof_correction()) *
1288  cube->dre().phase_correction();
1289 
1290  // Setup rotation matrix for conversion towards sky coordinates
1291  GMatrix ry;
1292  GMatrix rz;
1293  ry.eulery(model_dir.dec_deg() - 90.0);
1294  rz.eulerz(-model_dir.ra_deg());
1295  GMatrix rot = (ry * rz).transpose();
1296 
1297  // Loop over Chi and Psi
1298  for (int ipix = 0; ipix < npix; ++ipix) {
1299 
1300  // Get pointer to event bin
1301  const GCOMEventBin* bin = (*cube)[ipix];
1302 
1303  // Get reference to instrument direction
1304  const GCOMInstDir& obsDir = bin->dir();
1305 
1306  // Get reference to Phigeo sky direction
1307  const GSkyDir& phigeoDir = obsDir.dir();
1308 
1309  // Compute angle between model centre and Phigeo sky direction (radians)
1310  // and position angle of model with respect to Phigeo sky direction (radians)
1311  double zeta = phigeoDir.dist(model_dir);
1312 
1313  // Set radial model zenith angle range
1314  double rho_min = (zeta > phigeo_max) ? zeta - phigeo_max : 0.0;
1315  double rho_max = zeta + phigeo_max;
1316  if (rho_min > theta_max) {
1317  rho_min = theta_max;
1318  }
1319  if (rho_max > theta_max) {
1320  rho_max = theta_max;
1321  }
1322 
1323  // Perform model zenith angle integration if interval is valid
1324  if (rho_max > rho_min) {
1325 
1326  // Initialise IRF vector
1327  GVector irf(nphibar);
1328 
1329  // Setup integration kernel
1330  com_elliptical_kerns_rho integrands(m_iaq,
1331  model,
1332  irf,
1333  bin,
1334  rot,
1335  drx,
1336  phigeo_bin,
1337  m_phigeo_bins,
1338  nphibar,
1339  iter_omega);
1340 
1341  // Setup integrator
1342  GIntegrals integral(&integrands);
1343  integral.fixed_iter(iter_rho);
1344 
1345  // Integrate over Phigeo
1346  irf = integral.romberg(rho_min, rho_max, iter_rho);
1347 
1348  // Add IRF to result
1349  for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1350  int idri = ipix + iphibar * npix;
1351  irfs[idri] = irf[iphibar] * drg[idri] * iaq_norm;
1352  }
1353 
1354  } // endif: Phigeo angle interval was valid
1355 
1356  } // endfor: looped over Chi and Psi pixels
1357 
1358  // Return IRF vector
1359  return irfs;
1360 }
1361 
1362 
1363 /***********************************************************************//**
1364  * @brief Return instrument response to diffuse source
1365  *
1366  * @param[in] model Sky model.
1367  * @param[in] obs Observation.
1368  * @param[out] gradients Gradients matrix.
1369  * @return Instrument response to diffuse source for all events in
1370  * observation (\f$cm^2\f$).
1371  *
1372  * @exception GException::invalid_argument
1373  * Observation is not a COMPTEL observation.
1374  * @exception GException::invalid_value
1375  * Response not initialised with a valid IAQ
1376  * Incompatible IAQ encountered
1377  *
1378  * Returns the instrument response to a diffuse source for all events in
1379  * the observations. The diffuse source may be energy dependent.
1380  *
1381  * The computation is done by integrating the diffuse model for each pixel
1382  * in Chi and Psi over a circular region centred on the Chi/Psi pixel with
1383  * a radius equal to the maximum Phigeo value. The radial integration is
1384  * done by looping over all Phigeo bins of the response. For each Phigeo
1385  * value, the azimuthal integration is done by stepping with an angular step
1386  * size that corresponds to the Phigeo step size (which typically is 1 deg).
1387  *
1388  * @p gradients is an optional matrix where the number of rows corresponds
1389  * to the number of events in the observation and the number of columns
1390  * corresponds to the number of spatial model parameters. Since for point
1391  * sources no gradients are computed, the method does not alter the
1392  * content of @p gradients.
1393  ***************************************************************************/
1395  const GObservation& obs,
1396  GMatrix* gradients) const
1397 {
1398  // Extract COMPTEL observation
1399  const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
1400  if (obs_ptr == NULL) {
1401  std::string cls = std::string(typeid(&obs).name());
1402  std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
1403  "observations. Please specify a COMPTEL observation "
1404  "as argument.";
1406  }
1407 
1408  // Extract COMPTEL event cube
1409  const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(obs_ptr->events());
1410  if (cube == NULL) {
1411  std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
1412  "not contain a COMPTEL event cube. Please specify "
1413  "a COMPTEL observation containing and event cube.";
1415  }
1416 
1417  // Throw an exception if COMPTEL response is not set or if
1418  if (m_iaq.empty()) {
1419  std::string msg = "COMPTEL response is empty. Please initialise the "
1420  "response with an \"IAQ\".";
1422  }
1423  else if (m_phigeo_bin_size == 0.0) {
1424  std::string msg = "COMPTEL response has a zero Phigeo bin size. "
1425  "Please initialise the response with a valid "
1426  "\"IAQ\".";
1428  }
1429 
1430  // Get number of Chi/Psi pixels, Phibar layers and event bins
1431  int npix = cube->naxis(0) * cube->naxis(1);
1432  int nphibar = cube->naxis(2);
1433  int nevents = cube->size();
1434 
1435  // Throw an exception if the number of Phibar bins does not match the
1436  // response
1437  if (nphibar != m_phibar_bins) {
1438  std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
1439  "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
1440  "bins. Please specify a compatible IAQ.";
1442  }
1443 
1444  // Initialise result
1445  GVector irfs(nevents);
1446 
1447  // Initialise some variables
1448  double phigeo_min = m_phigeo_min * gammalib::deg2rad;
1449  double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
1450  double omega0 = 2.0 * std::sin(0.5 * phigeo_bin);
1451  const GSkyMap& drx = obs_ptr->drx().map();
1452  const double* drg = obs_ptr->drg().map().pixels();
1453 
1454  // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
1455  double iaq_norm = obs_ptr->deadc() /
1456  (obs_ptr->ontime() * cube->dre().tof_correction()) *
1457  cube->dre().phase_correction();
1458 
1459  // Loop over Chi and Psi
1460  for (int ipix = 0; ipix < npix; ++ipix) {
1461 
1462  // Get pointer to event bin
1463  const GCOMEventBin* bin = (*cube)[ipix];
1464 
1465  // Get reference to instrument direction
1466  const GCOMInstDir& obsDir = bin->dir();
1467 
1468  // Get reference to Phigeo sky direction
1469  const GSkyDir& phigeoDir = obsDir.dir();
1470 
1471  // Loop over Phigeo
1472  for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
1473 
1474  // Determine Phigeo value in radians. Note that phigeo_min is the value
1475  // for bin zero.
1476  double phigeo = phigeo_min + double(iphigeo) * phigeo_bin;
1477  double sin_phigeo = std::sin(phigeo);
1478 
1479  // Determine number of azimuthal integration steps and step size
1480  // in radians
1481  double length = gammalib::twopi * sin_phigeo;
1482  int naz = int(length / phigeo_bin + 0.5);
1483  if (naz < 2) {
1484  naz = 2;
1485  }
1486  double az_step = gammalib::twopi / double(naz);
1487 
1488  // Computes solid angle of integration bin multiplied by Jaccobian
1489  double omega = omega0 * az_step * sin_phigeo;
1490 
1491  // Loop over azimuth angle
1492  double az = 0.0;
1493  for (int iaz = 0; iaz < naz; ++iaz, az += az_step) {
1494 
1495  // Get sky direction
1496  GSkyDir skyDir = phigeoDir;
1497  skyDir.rotate(az, phigeo);
1498 
1499  // Fall through if sky direction is not contained in DRX
1500  if (!drx.contains(skyDir)) {
1501  continue;
1502  }
1503 
1504  // Set photon
1505  GPhoton photon(skyDir, bin->energy(), bin->time());
1506 
1507  // Get model sky intensity for photon (unit: sr^-1)
1508  double intensity = model.spatial()->eval(photon);
1509 
1510  // Fall through if intensity is zero
1511  if (intensity == 0.0) {
1512  continue;
1513  }
1514 
1515  // Multiply intensity by DRX value (unit: cm^2 s sr^-1)
1516  intensity *= drx(skyDir);
1517 
1518  // Multiply intensity by solid angle (unit: cm^2 s)
1519  intensity *= omega;
1520 
1521  // Loop over Phibar
1522  for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1523 
1524  // Get IAQ index
1525  int i = iphibar * m_phigeo_bins + iphigeo;
1526 
1527  // Get IAQ value
1528  double iaq = m_iaq[i];
1529 
1530  // Fall through if IAQ is not positive
1531  if (iaq <= 0.0) {
1532  continue;
1533  }
1534 
1535  // Get DRI index
1536  int idri = ipix + iphibar * npix;
1537 
1538  // Compute IRF value (unit: cm^2)
1539  double irf = iaq * drg[idri] * iaq_norm * intensity;
1540 
1541  // Add IRF value if it is positive
1542  if (irf > 0.0) {
1543  irfs[idri] += irf;
1544  }
1545 
1546  } // endfor: looped over Phibar
1547 
1548  } // endfor: looped over azimuth angle around locus of IAQ
1549 
1550  } // endfor: looped over Phigeo angles
1551 
1552  } // endfor: looped over Chi and Psi pixels
1553 
1554  // Return IRF vector
1555  return irfs;
1556 }
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.
double m_phibar_ref_value
Phigeo reference value (deg)
void save(const GFilename &filename, const bool &clobber=false) const
Save the response vector cache into FITS file.
const std::string & rspname(void) const
Return response name.
GCaldb m_caldb
Calibration database.
Sky map class.
Definition: GSkyMap.hpp:89
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:286
#define G_EBOUNDS
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
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.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:256
Point source spatial model class interface definition.
double m_phigeo_bin_size
Phigeo binsize (deg)
Energy value class definition.
Abstract elliptical spatial model base class.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const =0
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.
Single precision FITS image class definition.
Generic matrix class definition.
const GCOMDri & drx(void) const
Return exposure.
double m_phibar_min
Phigeo value of first bin (deg)
Abstract radial spatial model base class interface definition.
#define G_IRF_DIFFUSE
void init_members(void)
Initialise class members.
virtual int size(void) const
Return number of bins in event cube.
virtual double ontime(void) const
Return ontime.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
Kernel for rho angle integration of elliptical models.
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:524
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegrals.hpp:179
const double & tof_correction(void) const
Return ToF correction factor.
Definition: GCOMDri.hpp:398
const GSkyDir & dir(void) const
Return position of elliptical spatial model.
virtual GVector irf_ptsrc(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to point source.
Interface for the COMPTEL instrument response function.
double m_phigeo_ref_value
Phigeo reference value (deg)
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
virtual const GTime & time(void) const
Return time of event bin.
FITS file class interface definition.
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:787
COMPTEL event bin class.
double m_phibar_ref_pixel
Phigeo reference pixel (starting from 1)
virtual double deadc(const GTime &time=GTime()) const
Return deadtime correction factor.
Source class definition.
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1157
virtual double theta_max(void) const =0
GResponseVectorCache m_irf_vector_cache
Definition: GResponse.hpp:326
void id(const std::string &id)
Set observation identifier.
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:201
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:184
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
const double & phase_correction(void) const
Return pulsar phase correction factor.
Definition: GCOMDri.hpp:429
const double deg2rad
Definition: GMath.hpp:43
const GCOMDri & drg(void) const
Return geometry factors.
Calibration database class.
Definition: GCaldb.hpp:66
const GSkyDir & dir(void) const
Return position of radial spatial model.
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
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1194
Integration class for set of functions interface definition.
const GSkyMap & map(void) const
Return DRI sky map.
Definition: GCOMDri.hpp:278
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
void rotate(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:378
Abstract event base class definition.
void write(GFitsImageFloat &image) const
Write COMPTEL response into FITS image.
bool contains(const GSkyDir &dir) const
Checks if sky direction falls in map.
Definition: GSkyMap.cpp:2023
virtual double deadc(const GTime &time=GTime()) const =0
void load_cache(const GFilename &filename)
Load response cache.
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:229
virtual GVector irf_diffuse(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to diffuse source.
#define G_IRF_PTSRC
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
GChatter
Definition: GTypemaps.hpp:33
void save_cache(const GFilename &filename) const
Save response cache.
bool is_fits(void) const
Checks whether file is a FITS file.
Definition: GFilename.cpp:313
Abstract observation base class.
Vector class interface definition.
Single precision FITS image class.
COMPTEL event bin container 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:835
Defintion of COMPTEL helper classes for vector response.
Abstract elliptical spatial model base class interface definition.
int naxes(const int &axis) const
Return dimension of an image axis.
Definition: GFitsImage.cpp:344
const GCOMDri & dre(void) const
Return reference to DRE data.
virtual GResponse & operator=(const GResponse &rsp)
Assignment operator.
Definition: GResponse.cpp:139
void name(const std::string &name)
Set observation name.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:368
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)
Sky model class interface definition.
Sky model class.
Definition: GModelSky.hpp:122
#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.
Kernel for rho angle integration of radial models.
#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:1316
virtual GEvents * events(void)
Return events.
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:657
GVector romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegrals.cpp:210
virtual double theta_max(void) const =0
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:271
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:295
const double fourpi
Definition: GMath.hpp:37
Abstract radial spatial model base class.
const double twopi
Definition: GMath.hpp:36
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:393
virtual GVector irf_elliptical(const GModelSky &model, const GObservation &obs, GMatrix *gradients=NULL) const
Return instrument response to elliptical source.
void copy_members(const GCOMResponse &rsp)
Copy class members.
COMPTEL instrument direction class definition.
Integration class for set of functions.
Definition: GIntegrals.hpp:47
Abstract instrument response base class.
Definition: GResponse.hpp:77
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
COMPTEL event bin container class interface definition.
virtual ~GCOMResponse(void)
Destructor.
#define G_IRF
int iter_rho(const double &rho_max, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of radial Romberg iterations.
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
Sky direction class.
Definition: GSkyDir.hpp:62
void load(const GFilename &filename)
Load response vector cache from FITS file.
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
virtual int naxis(const int &axis) const
Return number of bins in axis.
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
const double * pixels(void) const
Returns pointer to pixel data.
Definition: GSkyMap.hpp:477
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:257
Filename class interface definition.
Interface for the COMPTEL instrument direction class.
Definition: GCOMInstDir.hpp:45
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
void load(const std::string &rspname)
Load COMPTEL response.
virtual const GEnergy & energy(void) const
Return energy of event bin.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489