GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMDri.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMDri.cpp - COMPTEL Data Space class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-2019 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 GCOMDri.cpp
23  * @brief COMPTEL Data Space class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GMath.hpp"
33 #include "GWcs.hpp"
34 #include "GFits.hpp"
35 #include "GFitsImage.hpp"
36 #include "GModelSky.hpp"
37 #include "GModelSpatial.hpp"
39 #include "GCOMDri.hpp"
40 #include "GCOMOad.hpp"
41 #include "GCOMOads.hpp"
42 #include "GCOMTim.hpp"
43 #include "GCOMStatus.hpp"
44 #include "GCOMObservation.hpp"
45 #include "GCOMEventList.hpp"
46 #include "GCOMSelection.hpp"
47 #include "GCOMSupport.hpp"
48 
49 /* __ Method name definitions ____________________________________________ */
50 #define G_COMPUTE_DRE "GCOMDri::compute_dre(GCOMObservation&, "\
51  "GCOMSelection&, double&)"
52 #define G_COMPUTE_DRM "GCOMDri::compute_drm(GCOMObservation&, GModel&)"
53 #define G_COMPUTE_DRE_PTSRC "GCOMDri::compute_drm_ptsrc(GCOMObservation&, "\
54  "GModelSky&)"
55 
56 /* __ Macros _____________________________________________________________ */
57 
58 /* __ Coding definitions _________________________________________________ */
59 //#define G_CHECK_EHA_COMPUTATION
60 
61 /* __ Debug definitions __________________________________________________ */
62 //#define G_DEBUG_COMPUTE_DRE
63 //#define G_DEBUG_COMPUTE_DRG
64 //#define G_DEBUG_COMPUTE_DRX
65 
66 /* __ Constants __________________________________________________________ */
67 const double superpacket_duration = 16.384; // Duration of superpacket (s)
68 const double r1 = 13.8; // D1 module radius (cm)
69 
70 /* __ Derived constants __________________________________________________ */
71 const double r1sq = r1 * r1;
72 const double d1_area = 7.0 * gammalib::pi * r1sq;
73 
74 
75 /*==========================================================================
76  = =
77  = Constructors/destructors =
78  = =
79  ==========================================================================*/
80 
81 /***********************************************************************//**
82  * @brief Void constructor
83  ***************************************************************************/
85 {
86  // Initialise class members
87  init_members();
88 
89  // Return
90  return;
91 }
92 
93 
94 /***********************************************************************//**
95  * @brief File name constructor
96  *
97  * @param[in] filename COMPTEL Data Space FITS file name.
98  ***************************************************************************/
99 GCOMDri::GCOMDri(const GFilename& filename)
100 {
101  // Initialise class members
102  init_members();
103 
104  // Load DRI FITS file
105  load(filename);
106 
107  // Return
108  return;
109 }
110 
111 
112 /***********************************************************************//**
113  * @brief Sky map constructor
114  *
115  * @param[in] map Sky map defining the DRI cube.
116  * @param[in] phimin Minimum Phibar angle (deg).
117  * @param[in] phibin Bin size of Phibar angle (deg).
118  * @param[in] nphibin Number of Phibar bins.
119  *
120  * Constructs a DRI cube from a sky map and a Phibar binning definition.
121  ***************************************************************************/
123  const double& phimin,
124  const double& phibin,
125  const int& nphibin)
126 {
127  // Initialise class members
128  init_members();
129 
130  // Set sky map
131  m_dri = map;
132 
133  // Make sure that sky map can hold the required number of maps
134  if (nphibin > 0) {
135  m_dri.nmaps(nphibin);
136  }
137 
138  // Set all sky map pixels to zero
139  m_dri = 0;
140 
141  // Store minimum Phibar and bin size
142  m_phimin = phimin;
143  m_phibin = phibin;
144 
145  // Return
146  return;
147 }
148 
149 
150 /***********************************************************************//**
151  * @brief Copy constructor
152  *
153  * @param[in] dri COMPTEL Data Space.
154  ***************************************************************************/
156 {
157  // Initialise class members
158  init_members();
159 
160  // Copy members
161  copy_members(dri);
162 
163  // Return
164  return;
165 }
166 
167 
168 /***********************************************************************//**
169  * @brief Destructor
170  ***************************************************************************/
172 {
173  // Free members
174  free_members();
175 
176  // Return
177  return;
178 }
179 
180 
181 /*==========================================================================
182  = =
183  = Operators =
184  = =
185  ==========================================================================*/
186 
187 /***********************************************************************//**
188  * @brief Assignment operator
189  *
190  * @param[in] dri COMPTEL Data Space.
191  * @return COMPTEL Data Space.
192  ***************************************************************************/
194 {
195  // Execute only if object is not identical
196  if (this != &dri) {
197 
198  // Free members
199  free_members();
200 
201  // Initialise private members
202  init_members();
203 
204  // Copy members
205  copy_members(dri);
206 
207  } // endif: object was not identical
208 
209  // Return this object
210  return *this;
211 }
212 
213 
214 /*==========================================================================
215  = =
216  = Public methods =
217  = =
218  ==========================================================================*/
219 
220 /***********************************************************************//**
221  * @brief Clear COMPTEL Data Space
222  ***************************************************************************/
223 void GCOMDri::clear(void)
224 {
225  // Free members
226  free_members();
227 
228  // Initialise private members
229  init_members();
230 
231  // Return
232  return;
233 }
234 
235 
236 /***********************************************************************//**
237  * @brief Clone COMPTEL Data Space
238  *
239  * @return Pointer to deep copy of COMPTEL Data Space.
240  ***************************************************************************/
242 {
243  return new GCOMDri(*this);
244 }
245 
246 
247 /***********************************************************************//**
248  * @brief Compute event cube
249  *
250  * @param[in] obs COMPTEL observation.
251  * @param[in] select Selection set.
252  * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
253  *
254  * @exception GException::invalid_argument
255  * DRE cube has a non-positive Phibar bin size.
256  *
257  * Compute DRE event cube for a COMPTEL observation.
258  ***************************************************************************/
260  const GCOMSelection& select,
261  const double& zetamin)
262 {
263  // Debug
264  #if defined(G_DEBUG_COMPUTE_DRE)
265  std::cout << "GCOMDri::compute_dre" << std::endl;
266  std::cout << "====================" << std::endl;
267  #endif
268 
269  // Check if observation contains a COMPTEL event list
270  const GCOMEventList* events = dynamic_cast<const GCOMEventList*>(obs.events());
271  if (events == NULL) {
272  std::string msg = "No event list found in COMPTEL observation. Please "
273  "specify an observation that contains an event list.";
275  }
276 
277  // Check for positive Phibar bin size
278  if (m_phibin <= 0.0) {
279  std::string msg = "DRE cube has a non-positive Phibar bin size. Please "
280  "specify a DRE cube with a positive Phibar bin size.";
282  }
283 
284  // Initialise event counter
285  int i_evt = 0;
286 
287  // Initialise D1 & D2 module status
288  GCOMStatus status;
289 
290  // Initialise superpacket statistics
291  init_statistics();
292 
293  // Initialise selection statistics
294  select.init_statistics();
295 
296  // Initialise statistics
297  int num_used_events = 0;
298  int num_event_outside_sp = 0;
299  int num_energy_too_low = 0;
300  int num_energy_too_high = 0;
301  int num_eha_too_small = 0;
302  int num_phibar_too_low = 0;
303  int num_phibar_too_high = 0;
304  int num_outside_dre = 0;
305  int num_event_before_dre = 0;
306  int num_event_after_dre = 0;
307  int num_d1module_off = 0;
308  int num_d2module_off = 0;
309  int num_processed = 0;
310 
311  // Set all DRE bins to zero
312  init_cube();
313 
314  // Signal that loop should be terminated
315  bool terminate = false;
316 
317  // Loop over Orbit Aspect Data
318  for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
319 
320  // Get reference to Orbit Aspect Data of superpacket
321  const GCOMOad &oad = obs.oads()[i_oad];
322 
323  // Check superpacket usage
324  if (!use_superpacket(oad, obs.tim())) {
325  continue;
326  }
327 
328  // Prepare Earth horizon angle comparison
329  GSkyDir sky_geocentre;
330  double theta_geocentre = double(oad.gcel());
331  double phi_geocentre = double(oad.gcaz());
332  sky_geocentre.radec_deg(phi_geocentre, 90.0-theta_geocentre);
333 
334  // Collect all events within superpacket. Break if the end
335  // of the event list was reached.
336  for (; i_evt < events->size(); ++i_evt) {
337 
338  // Get pointer to event
339  const GCOMEventAtom* event = (*events)[i_evt];
340 
341  // Break loop if the end of the superpacket was reached
342  if (event->time() > oad.tstop()) {
343  break;
344  }
345 
346  // Increase number of processed events
347  num_processed++;
348 
349  // Check GTIs if the DRE has GTIs
350  if (m_gti.size() > 0) {
351 
352  // Skip event if it lies before the DRE start.
353  if (event->time() < m_gti.tstart()) {
354  num_event_before_dre++;
355  continue;
356  }
357 
358  // Break if event lies after the DRE stop
359  else if (event->time() > m_gti.tstop()) {
360  num_event_after_dre++;
361  terminate = true;
362  break;
363  }
364  }
365 
366  // Skip event if it lies before the superpacket start
367  if (event->time() < oad.tstart()) {
368  num_event_outside_sp++;
369  continue;
370  }
371 
372  // Skip event if it lies outside energy range
373  if (event->energy() < m_ebounds.emin()) {
374  num_energy_too_low++;
375  continue;
376  }
377  else if (event->energy() > m_ebounds.emax()) {
378  num_energy_too_high++;
379  continue;
380  }
381 
382  // Apply event selectio
383  if (!select.use_event(*event)) {
384  continue;
385  }
386 
387  // Extract module IDs from MODCOM
388  int id2 = (event->modcom()-1)/7 + 1; // [1-14]
389  int id1 = event->modcom() - (id2-1) * 7; // [1-7]
390 
391  // Check module status. Skip all events where module is
392  // signalled as not on.
393  if (status.d1status(oad.tjd(), id1) != 1) {
394  num_d1module_off++;
395  continue;
396  }
397  if (status.d2status(oad.tjd(), id2) != 1) {
398  num_d2module_off++;
399  continue;
400  }
401 
402  // Compute Compton scatter angle index. Skip if it's invalid.
403  int iphibar = (event->phibar() - m_phimin) / m_phibin;
404  if (iphibar < 0) {
405  num_phibar_too_low++;
406  continue;
407  }
408  else if (iphibar >= nphibar()) {
409  num_phibar_too_high++;
410  continue;
411  }
412 
413  // Option: Earth horizon angle comparison
414  #if defined(G_CHECK_EHA_COMPUTATION)
415  GSkyDir sky_event;
416  double theta_event = double(event->theta());
417  double phi_event = double(event->phi());
418  sky_event.radec_deg(-phi_event, 90.0-theta_event);
419  double eha = sky_geocentre.dist_deg(sky_event) - oad.georad();
420  if (std::abs(eha - event->eha()) > 1.5) {
421  std::string msg = "Earth horizon angle from EVP dataset ("+
422  gammalib::str(event->eha())+" deg) "
423  "differs from Earth horizon angle "
424  "computed from Orbit Aspect Data ("+
425  gammalib::str(eha)+" deg). Use the EVP "
426  "value.";
428  }
429  #endif
430 
431  // Check for Earth horizon angle. There is a constant EHA limit
432  // over a Phibar layer to be compliant with the DRG.
433  double ehamin = double(iphibar) * m_phibin + zetamin;
434  if (event->eha() < ehamin) {
435  num_eha_too_small++;
436  continue;
437  }
438 
439  // Now fill the event into the DRE. Put this in a try-catch
440  // block so that any invalid transformations are catched.
441  try {
442  GSkyPixel pixel = m_dri.dir2pix(event->dir().dir());
443  if (m_dri.contains(pixel)) {
444  int inx = m_dri.pix2inx(pixel) + iphibar * m_dri.npix();
445  (*this)[inx] += 1.0;
446  num_used_events++;
447  }
448  else {
449  num_outside_dre++;
450  }
451  }
453  num_outside_dre++;
454  }
455 
456  } // endfor: collected events
457 
458  // Break if termination was signalled or if there are no more events
459  if (terminate || i_evt >= events->size()) {
460  break;
461  }
462 
463  } // endfor: looped over Orbit Aspect Data
464 
465  // Set Good Time interval for DRE
466  if (m_num_used_superpackets > 0) {
468  }
469 
470  // Debug
471  #if defined(G_DEBUG_COMPUTE_DRE)
472  std::cout << "Total number of superpackets .: " << m_num_superpackets << std::endl;
473  std::cout << "Used superpackets ............: " << m_num_used_superpackets << std::endl;
474  std::cout << "Skipped superpackets .........: " << m_num_skipped_superpackets << std::endl;
475  std::cout << "Total number of events .......: " << events.size() << std::endl;
476  std::cout << "Processed events .............: " << num_processed << std::endl;
477  std::cout << "Used events ..................: " << num_used_events << std::endl;
478  std::cout << "Events outside superpacket ...: " << num_event_outside_sp << std::endl;
479  std::cout << "Events before DRE GTI ........: " << num_event_before_dre << std::endl;
480  std::cout << "Events after DRE GTI .........: " << num_event_after_dre << std::endl;
481  std::cout << "Events with D1 module off ....: " << num_d1module_off << std::endl;
482  std::cout << "Events with D2 module off ....: " << num_d2module_off << std::endl;
483  std::cout << "Energy too low ...............: " << num_energy_too_low << std::endl;
484  std::cout << "Energy too high ..............: " << num_energy_too_high << std::endl;
485  std::cout << "Phibar too low ...............: " << num_phibar_too_low << std::endl;
486  std::cout << "Phibar too high ..............: " << num_phibar_too_high << std::endl;
487  std::cout << "Outside DRE cube .............: " << num_outside_dre << std::endl;
488  std::cout << "Earth horizon angle too small : " << num_eha_too_small << std::endl;
489  std::cout << select << std::endl;
490  #endif
491 
492  // Return
493  return;
494 }
495 
496 
497 /***********************************************************************//**
498  * @brief Compute geometry cube
499  *
500  * @param[in] obs COMPTEL observation.
501  * @param[in] select Selection set (not used so far).
502  * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
503  *
504  * Compute DRG cube for a COMPTEL observation.
505  ***************************************************************************/
507  const GCOMSelection& select,
508  const double& zetamin)
509 {
510  // Debug
511  #if defined(G_DEBUG_COMPUTE_DRG)
512  std::cout << "GCOMDri::compute_drg" << std::endl;
513  std::cout << "====================" << std::endl;
514  #endif
515 
516  // Initialise D1 & D2 module status
517  GCOMStatus status;
518 
519  // Initialise superpacket statistics
520  init_statistics();
521 
522  // Initialise selection statistics
523  select.init_statistics();
524 
525  // Set all DRG bins to zero
526  init_cube();
527 
528  // Loop over Orbit Aspect Data
529  for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
530 
531  // Get reference to Orbit Aspect Data of superpacket
532  const GCOMOad &oad = obs.oads()[i_oad];
533 
534  // Check superpacket usage
535  if (!use_superpacket(oad, obs.tim())) {
536  continue;
537  }
538 
539  // Prepare Earth horizon angle computation. The celestial system
540  // is reinterpreted as the COMPTEL coordinate system, where the
541  // 90 degrees - zenith angle becomes the declination and the azimuth
542  // angle becomes Right Ascension. This allows us later to use the
543  // GSkyDir::dist method to compute the distance between the geocentre
544  // and the telescope Z-axis.
545  double theta_geocentre = double(oad.gcel());
546  double phi_geocentre = double(oad.gcaz());
547  GSkyDir geocentre_comptel;
548  geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
549 
550  // Compute geometrical factors for all (Chi, Psi)
551  for (int index = 0; index < m_dri.npix(); ++index) {
552 
553  // Get sky direction for (Chi, Psi)
554  GSkyDir sky = m_dri.inx2dir(index);
555 
556  // Convert sky direction to COMPTEL coordinates
557  double theta = oad.theta(sky);
558  double phi = oad.phi(sky);
559 
560  // Compute geometric factor summed over all D1, D2
561  double geometry = compute_geometry(oad.tjd(), theta, phi, status);
562 
563  // Compute Earth horizon angle as the distance between the Earth
564  // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
565  // COMPTEL coordinates minus the radius of the Earth.
566  GSkyDir chipsi_comptel;
567  chipsi_comptel.radec_deg(phi, 90.0-theta);
568  double eha = geocentre_comptel.dist_deg(chipsi_comptel) - oad.georad();
569 
570  // Loop over Phibar
571  for (int iphibar = 0; iphibar < nphibar(); ++iphibar) {
572 
573  // Compute minimum Earth horizon angle for Phibar layer
574  double ehamin = double(iphibar) * m_phibin + zetamin;
575 
576  // Add up geometry if the Earth horizon angle is equal or
577  // larger than the minimum
578  if (eha >= ehamin) {
579  int inx = index + iphibar * m_dri.npix();
580  (*this)[inx] += geometry;
581  }
582 
583  } // endfor: looped over Phibar
584 
585  } // endfor: looped over (Chi, Psi)
586 
587  } // endfor: looped over Orbit Aspect Data
588 
589  // Divide DRG by number of used superpackets
590  if (m_num_used_superpackets > 0) {
591  double norm = 1.0 / double(m_num_used_superpackets);
592  for (int i = 0; i < size(); ++i) {
593  (*this)[i] *= norm;
594  }
595  }
596 
597  // Clear energy boundaries for DRG since DRG does not depend on energy)
598  m_ebounds.clear();
599 
600  // Set the Good Time interval for DRG
601  if (m_num_used_superpackets > 0) {
603  }
604 
605  // Debug
606  #if defined(G_DEBUG_COMPUTE_DRG)
607  std::cout << "Total number of superpackets .: " << m_num_superpackets << std::endl;
608  std::cout << "Used superpackets ............: " << m_num_used_superpackets << std::endl;
609  std::cout << "Skipped superpackets .........: " << m_num_skipped_superpackets << std::endl;
610  #endif
611 
612  // Return
613  return;
614 }
615 
616 
617 /***********************************************************************//**
618  * @brief Compute DRX exposure map
619  *
620  * @param[in] obs COMPTEL observation.
621  *
622  * Compute DRG cube for a COMPTEL observation.
623  *
624  * Compute DRX exposure map for a COMPTEL observation.
625  *
626  * For a given superpacket, the exposure is computed using
627  *
628  * \f[
629  * X_i(\theta_c) = 7 \pi r_1^2 \cos \theta_c
630  * \frac{1 - \exp \left( -\tau \ \cos \theta_c \right)}
631  * {1 - \exp \left( -\tau \right)}
632  * \f]
633  *
634  * where
635  * \f$\tau=0.2\f$ is the typical thickness of a D1 module in radiation
636  * lengths,
637  * \f$r_1=13.8\f$ cm is the radius of a D1 module, and
638  * \f$\theta_c\f$ is the zenith angle in COMPTEL coordinates.
639  ***************************************************************************/
641 {
642  // Debug
643  #if defined(G_DEBUG_COMPUTE_DRX)
644  std::cout << "GCOMDri::compute_drx" << std::endl;
645  std::cout << "====================" << std::endl;
646  #endif
647 
648  // Initialise constants
649  const double tau = 0.2; // Thickness in radiation lenghts
650 
651  // Initialise superpacket statistics
652  init_statistics();
653 
654  // Set all DRX bins to zero
655  init_cube();
656 
657  // Loop over Orbit Aspect Data
658  for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
659 
660  // Get reference to Orbit Aspect Data of superpacket
661  const GCOMOad &oad = obs.oads()[i_oad];
662 
663  // Check superpacket usage
664  if (!use_superpacket(oad, obs.tim())) {
665  continue;
666  }
667 
668  // Loop over all DRX pixels
669  for (int index = 0; index < m_dri.npix(); ++index) {
670 
671  // Get sky direction for DRX pixel
672  GSkyDir sky = m_dri.inx2dir(index);
673 
674  // Compute zenith angle of pixel in COMPTEL coordinates
675  double theta = oad.theta(sky);
676 
677  // Initialise exposure
678  double exposure = 0.0;
679 
680  // Skip pixel if zenith angle is beyond 90 degrees
681  if (theta >= 90.0) {
682  continue;
683  }
684 
685  // ... otherwise compute the exposure
686  else {
687  double costheta = std::cos(theta * gammalib::deg2rad);
688  if (theta < 89.0) {
689  exposure = d1_area * costheta *
690  (1.0 - std::exp(-tau / costheta)) /
691  (1.0 - std::exp(-tau));
692  }
693  else {
694  exposure = d1_area * costheta / (1.0 - std::exp(-tau));
695  }
696  }
697 
698  // Accumulate exposure
699  (*this)[index] += exposure;
700 
701  } // endfor: looped over all DRX pixels
702 
703  } // endfor: looped over Orbit Aspect Data
704 
705  // Multiply by time per superpacket to give the result in cm^2 s
706  for (int i = 0; i < size(); ++i) {
707  (*this)[i] *= superpacket_duration;
708  }
709 
710  // Clear energy boundaries for DRX since DRX does not depend on energy)
711  m_ebounds.clear();
712 
713  // Set the Good Time interval for DRX
714  if (m_num_used_superpackets > 0) {
716  }
717 
718  // Debug
719  #if defined(G_DEBUG_COMPUTE_DRX)
720  std::cout << "Total number of superpackets .: " << m_num_superpackets << std::endl;
721  std::cout << "Used superpackets ............: " << m_num_used_superpackets << std::endl;
722  std::cout << "Skipped superpackets .........: " << m_num_skipped_superpackets << std::endl;
723  #endif
724 
725  // Return
726  return;
727 }
728 
729 
730 /***********************************************************************//**
731  * @brief Compute DRM model
732  *
733  * @param[in] obs COMPTEL observation.
734  * @param[in] model Model.
735  *
736  * @exception GException::feature_not_implemented
737  * Method is only implemented for a point source sky model.
738  *
739  * Compute DRM model cube for a COMPTEL observation.
740  ***************************************************************************/
742  const GModel& model)
743 {
744  // Check if model is a sky model
745  const GModelSky* skymodel = dynamic_cast<const GModelSky*>(&model);
746  if (skymodel == NULL) {
747  std::string msg = "Method is only implement for sky models.";
749  }
750 
751  // Extract spatial model component
752  const GModelSpatial* spatial = skymodel->spatial();
753 
754  // Select computation depending on the spatial model type
755  switch (spatial->code()) {
757  {
758  compute_drm_ptsrc(obs, *skymodel);
759  }
760  break;
764  {
765  std::string msg = "Method is not yet implemented for spatial model "
766  "type \""+spatial->type()+"\".";
768  }
769  break;
770  default:
771  break;
772  }
773 
774  // Return
775  return;
776 }
777 
778 
779 /***********************************************************************//**
780  * @brief Compute content in cone
781  *
782  * @param[in] dir Sky direction of cone apex.
783  * @param[in] armmin Minimum Angular Resolution Measure (deg).
784  * @param[in] armmax Maximum Angular Resolution Measure (deg).
785  * @return Content in cone.
786  *
787  * Compute the sum of the DRI bins within an event cone with apex at a given
788  * sky direction. All bins with an Angular Resolution Measure comprised
789  * between @p armmin (inclusive) and @p armmax (exclusive) will be
790  * considered. The bin centres will be used for the computation of the
791  * Angular Resolution Measure. The Angular Resolution Measure is defined as
792  * phibar - phigeo.
793  ***************************************************************************/
794 double GCOMDri::cone_content(const GSkyDir& dir,
795  const double& armmin,
796  const double& armmax) const
797 {
798  // Initialise content
799  double content = 0.0;
800 
801  // Create Phigeo map in degrees
802  GSkyMap phigeo = m_dri.extract(0);
803  for (int i = 0; i < phigeo.npix(); ++i) {
804  phigeo(i) = dir.dist_deg(phigeo.inx2dir(i));
805  }
806 
807  // Loop over phibar layers
808  for (int iphibar = 0; iphibar < this->nphibar(); ++iphibar) {
809 
810  // Compute phibar
811  double phibar = this->phimin() + (iphibar+0.5) * this->phibin();
812 
813  // Compute index offset
814  int offset = iphibar * phigeo.npix();
815 
816  // Loop over bins in phibar layer
817  for (int i = 0; i < phigeo.npix(); ++i) {
818 
819  // Compute ARM
820  double arm = phibar - phigeo(i);
821 
822  // If ARM is within specified interval then add bin content
823  if (arm >= armmin && arm < armmax) {
824  content += (*this)[i + offset];
825  }
826 
827  } // endfor: looped over bins in phibar layers
828 
829  } // endfor: looped over phibar layers
830 
831  // Return content
832  return content;
833 }
834 
835 
836 /***********************************************************************//**
837  * @brief Load COMPTEL Data Space from DRI FITS file
838  *
839  * @param[in] filename DRI FITS file name.
840  ***************************************************************************/
841 void GCOMDri::load(const GFilename& filename)
842 {
843  // Open FITS file
844  GFits fits(filename);
845 
846  // Get HDU (pointer is always valid)
847  const GFitsImage& hdu = *fits.image(0);
848 
849  // Read DRI file
850  read(hdu);
851 
852  // Close FITS file
853  fits.close();
854 
855  // Return
856  return;
857 }
858 
859 
860 /***********************************************************************//**
861  * @brief Save COMPTEL Data Space into DRI FITS file
862  *
863  * @param[in] filename DRI FITS file name.
864  * @param[in] clobber Overwrite existing file?
865  ***************************************************************************/
866 void GCOMDri::save(const GFilename& filename, const bool& clobber) const
867 {
868  // Create FITS file
869  GFits fits;
870 
871  // Write data space into FITS file
872  write(fits, filename.extname(gammalib::extname_dri));
873 
874  // Save FITS file
875  fits.saveto(filename, clobber);
876 
877  // Close FITS file
878  fits.close();
879 
880  // Return
881  return;
882 }
883 
884 
885 /***********************************************************************//**
886  * @brief Read COMPTEL Data Space from DRI FITS image
887  *
888  * @param[in] image DRI FITS image.
889  ***************************************************************************/
890 void GCOMDri::read(const GFitsImage& image)
891 {
892  // Clear
893  clear();
894 
895  // Read sky map
896  m_dri.read(image);
897 
898  // Correct WCS projection (HEASARC data format kluge)
900 
901  // Read attributes
902  read_attributes(&image);
903 
904  // Return
905  return;
906 }
907 
908 
909 /***********************************************************************//**
910  * @brief Write COMPTEL Data Space into FITS image
911  *
912  * @param[in] fits FITS file.
913  * @param[in] extname Extension name.
914  ***************************************************************************/
915 void GCOMDri::write(GFits& fits, const std::string& extname) const
916 {
917  // Write sky map into FITS file
918  GFitsHDU *image = m_dri.write(fits, extname);
919 
920  // Write DRI attributes
921  if (image != NULL) {
922  write_attributes(image);
923  }
924 
925  // Return
926  return;
927 }
928 
929 
930 /***********************************************************************//**
931  * @brief Print COMPTEL Data Space
932  *
933  * @param[in] chatter Chattiness.
934  * @return String containing COMPTEL Data Space information.
935  ***************************************************************************/
936 std::string GCOMDri::print(const GChatter& chatter) const
937 {
938  // Initialise result string
939  std::string result;
940 
941  // Continue only if chatter is not silent
942  if (chatter != SILENT) {
943 
944  // Compute scatter angle dimensions
945  double chimin = 0.0;
946  double chimax = 0.0;
947  double chibin = 0.0;
948  double psimin = 0.0;
949  double psimax = 0.0;
950  double psibin = 0.0;
951  const GWcs* wcs = dynamic_cast<const GWcs*>(m_dri.projection());
952  if (wcs != NULL) {
953  chibin = wcs->cdelt(0);
954  chimin = wcs->crval(0) - (wcs->crpix(0)-0.5) * chibin;
955  chimax = chimin + m_dri.nx() * chibin;
956  psibin = wcs->cdelt(1);
957  psimin = wcs->crval(1) - (wcs->crpix(1)-0.5) * psibin;
958  psimax = psimin + m_dri.ny() * psibin;
959  }
960 
961  // Compute Phibar maximum
962  double phimax = m_phimin + m_dri.nmaps() * m_phibin;
963 
964  // Append header
965  result.append("=== GCOMDri ===");
966 
967  // Append Phibar information
968  result.append("\n"+gammalib::parformat("Chi range"));
969  result.append(gammalib::str(chimin)+" - ");
970  result.append(gammalib::str(chimax)+" deg");
971  result.append("\n"+gammalib::parformat("Chi bin size"));
972  result.append(gammalib::str(chibin)+" deg");
973  result.append("\n"+gammalib::parformat("Psi range"));
974  result.append(gammalib::str(psimin)+" - ");
975  result.append(gammalib::str(psimax)+" deg");
976  result.append("\n"+gammalib::parformat("Psi bin size"));
977  result.append(gammalib::str(psibin)+" deg");
978  result.append("\n"+gammalib::parformat("Phibar range"));
979  result.append(gammalib::str(m_phimin)+" - ");
980  result.append(gammalib::str(phimax)+" deg");
981  result.append("\n"+gammalib::parformat("Phibar bin size"));
982  result.append(gammalib::str(m_phibin)+" deg");
983 
984  // Append energy boundaries
985  result.append("\n"+m_ebounds.print(gammalib::reduce(chatter)));
986 
987  // Append GTI
988  result.append("\n"+m_gti.print(gammalib::reduce(chatter)));
989 
990  // EXPLICIT: Append sky map
991  if (chatter >= EXPLICIT) {
992  result.append("\n"+m_dri.print(gammalib::reduce(chatter)));
993  }
994 
995  // Append computation statistics
996  result.append("\n"+gammalib::parformat("Input superpackets"));
997  result.append(gammalib::str(m_num_superpackets));
998  result.append("\n"+gammalib::parformat("Used superpackets"));
999  result.append(gammalib::str(m_num_used_superpackets));
1000  result.append("\n"+gammalib::parformat("Skipped superpackets"));
1001  result.append(gammalib::str(m_num_skipped_superpackets));
1002 
1003  } // endif: chatter was not silent
1004 
1005  // Return result
1006  return result;
1007 }
1008 
1009 
1010 /*==========================================================================
1011  = =
1012  = Private methods =
1013  = =
1014  ==========================================================================*/
1015 
1016 /***********************************************************************//**
1017  * @brief Initialise class members
1018  ***************************************************************************/
1020 {
1021  // Initialise members
1022  m_name.clear();
1023  m_dri.clear();
1024  m_ebounds.clear();
1025  m_gti.clear();
1026  m_phimin = 0.0;
1027  m_phibin = 0.0;
1028 
1029  // Initialise statistics
1030  init_statistics();
1031 
1032  // Initialise selection parameters
1033  m_selection.clear();
1034  m_zetamin = -90.0; // Signals no selection
1035 
1036  // Return
1037  return;
1038 }
1039 
1040 
1041 /***********************************************************************//**
1042  * @brief Copy class members
1043  *
1044  * @param[in] dri COMPTEL Data Space.
1045  ***************************************************************************/
1047 {
1048  // Copy members
1049  m_name = dri.m_name;
1050  m_dri = dri.m_dri;
1051  m_ebounds = dri.m_ebounds;
1052  m_gti = dri.m_gti;
1053  m_phimin = dri.m_phimin;
1054  m_phibin = dri.m_phibin;
1055 
1056  // Copy statistics
1057  m_tstart = dri.m_tstart;
1058  m_tstop = dri.m_tstop;
1062 
1063  // Copy selection parameters
1064  m_selection = dri.m_selection;
1065  m_zetamin = dri.m_zetamin;
1066 
1067  // Return
1068  return;
1069 }
1070 
1071 
1072 /***********************************************************************//**
1073  * @brief Delete class members
1074  ***************************************************************************/
1076 {
1077  // Return
1078  return;
1079 }
1080 
1081 
1082 /***********************************************************************//**
1083  * @brief Initialise DRI cube
1084  *
1085  * Sets all DRI cube bins to zero.
1086  ***************************************************************************/
1088 {
1089  // Set all cube bins to zero
1090  for (int i = 0; i < size(); ++i) {
1091  (*this)[i] = 0.0;
1092  }
1093 
1094  // Return
1095  return;
1096 }
1097 
1098 
1099 /***********************************************************************//**
1100  * @brief Initialise computation statistics
1101  ***************************************************************************/
1103 {
1104  // Initialise statistics
1105  m_tstart.clear();
1106  m_tstop.clear();
1107  m_num_superpackets = 0;
1110 
1111  // Return
1112  return;
1113 }
1114 
1115 
1116 /***********************************************************************//**
1117  * @brief Check if superpacket should be used
1118  *
1119  * @param[in] oad Orbit Aspect Data record (i.e. superpacket).
1120  * @param[in] tim Good Time Intervals.
1121  * @return True if superpacket should be used, false otherwise.
1122  *
1123  * Checks if a superpacket should be used and updated the superpacket
1124  * statistics and selected time interval. A superpacket will be used if
1125  * it is fully enclosed within the COMPTEL Good Time Intervals and the
1126  * Good Time Intervals of the DRI dataset.
1127  ***************************************************************************/
1128 bool GCOMDri::use_superpacket(const GCOMOad &oad, const GCOMTim& tim)
1129 {
1130  // Initialise usage flag
1131  bool use = true;
1132 
1133  // Increment superpacket counter
1135 
1136  // Skip superpacket if it is not fully enclosed within the COMPTEL
1137  // Good Time Intervals
1138  if (!(tim.contains(oad.tstart()) && tim.contains(oad.tstop()))) {
1140  use = false;
1141  }
1142 
1143  // Skip superpacket if it is not fully enclosed within the DRI Good
1144  // Time Intervals. Only check if there are Good Time Intervals in the
1145  // DRI.
1146  else if ((m_gti.size() > 0) &&
1147  !(m_gti.contains(oad.tstart()) && m_gti.contains(oad.tstop()))) {
1149  use = false;
1150  }
1151 
1152  // ... otherwise use superpacket
1153  else {
1155  }
1156 
1157  // Update selection validity interval
1158  if (m_num_used_superpackets == 1) {
1159  m_tstart = oad.tstart();
1160  m_tstop = oad.tstop();
1161  }
1162  else {
1163  m_tstop = oad.tstop();
1164  }
1165 
1166  // Return usage
1167  return use;
1168 }
1169 
1170 
1171 /***********************************************************************//**
1172  * @brief Read DRI attributes from FITS HDU
1173  *
1174  * @param[in] hdu FITS HDU pointer.
1175  *
1176  * Reads the time interval from the FITS header and sets the Phibar definiton
1177  * and energy boundaries from the header keywords if they are provided.
1178  ***************************************************************************/
1180 {
1181  // Get time attributes
1182  GTime tstart = gammalib::com_time(hdu->integer("VISDAY"), hdu->integer("VISTIM"));
1183  GTime tstop = gammalib::com_time(hdu->integer("VIEDAY"), hdu->integer("VIETIM"));
1184 
1185  // Set Good Time Intervals
1186  if (tstop > tstart) {
1187  m_gti = GGti(tstart, tstop);
1188  }
1189 
1190  // Optionally read Phibar attributes
1191  if (hdu->has_card("CDELT3") &&
1192  hdu->has_card("CRVAL3") &&
1193  hdu->has_card("CRPIX3")) {
1194 
1195  // Get phibar attributes
1196  m_phibin = hdu->real("CDELT3");
1197  m_phimin = hdu->real("CRVAL3") - (hdu->real("CRPIX3")-0.5) * m_phibin;
1198 
1199  }
1200 
1201  // ... otherwise set Phibar attributes to zero
1202  else {
1203  m_phimin = 0.0;
1204  m_phibin = 0.0;
1205  }
1206 
1207  // Optionally read energy attributes
1208  if (hdu->has_card("E_MIN") && hdu->has_card("E_MAX")) {
1209 
1210  // Get energy attributes
1211  GEnergy emin(hdu->real("E_MIN"), "MeV");
1212  GEnergy emax(hdu->real("E_MAX"), "MeV");
1213 
1214  // Set energy boundaries
1215  m_ebounds = GEbounds(emin, emax);
1216 
1217  }
1218 
1219  // ... otherwise clear energy boundaries
1220  else {
1221  m_ebounds.clear();
1222  }
1223 
1224  // Optionally read superpacket statistics
1225  if (hdu->has_card("NSPINP")) {
1226  m_num_superpackets = hdu->integer("NSPINP");
1227  }
1228  if (hdu->has_card("NSPUSE")) {
1229  m_num_used_superpackets = hdu->integer("NSPUSE");
1230  }
1231  if (hdu->has_card("NSPSKP")) {
1232  m_num_skipped_superpackets = hdu->integer("NSPSKP");
1233  }
1234 
1235  // Return
1236  return;
1237 }
1238 
1239 
1240 /***********************************************************************//**
1241  * @brief Write DRI attributes into FITS HDU
1242  *
1243  * @param[in] hdu FITS HDU pointer.
1244  ***************************************************************************/
1246 {
1247  // Set Phibar keywords
1248  double crval3 = m_phimin + 0.5 * m_phibin;
1249 
1250  // Write Phibar keywords
1251  hdu->card("CTYPE3", "Phibar", "Compton scatter angle");
1252  hdu->card("CRPIX3", 1.0, "Pixel coordinate of reference point (starting from 1)");
1253  hdu->card("CRVAL3", crval3, "[deg] Coordinate value at reference point");
1254  hdu->card("CDELT3", m_phibin, "[deg] Coordinate increment at reference point");
1255 
1256  // Write OGIP time keywords
1257  m_gti.reference().write(*hdu);
1258  hdu->card("TSTART", m_gti.tstart().secs(), "[s] Start time");
1259  hdu->card("TSTOP", m_gti.tstop().secs(), "[s] Stop time");
1260  hdu->card("DATE-OBS", m_gti.tstart().utc(), "Start of observation in UTC");
1261  hdu->card("DATE-END", m_gti.tstop().utc(), "Stop of observation in UTC");
1262 
1263  // Set time keywords
1264  int visday = gammalib::com_tjd(m_gti.tstart());
1265  int vistim = gammalib::com_tics(m_gti.tstart());
1266  int vieday = gammalib::com_tjd(m_gti.tstop());
1267  int vietim = gammalib::com_tics(m_gti.tstop());
1268 
1269  // Write COMPTEL time keywords
1270  hdu->card("VISDAY", visday, "[TJD] Data validity interval start day");
1271  hdu->card("VISTIM", vistim, "[tics] Data validity interval start time");
1272  hdu->card("VIEDAY", vieday, "[TJD] Data validity interval end day");
1273  hdu->card("VIETIM", vietim, "[tics] Data validity interval start time");
1274 
1275  // If there are energy boundaries then write them
1276  if (m_ebounds.size() > 0) {
1277  hdu->card("E_MIN", m_ebounds.emin().MeV(), "[MeV] Lower bound of energy range");
1278  hdu->card("E_MAX", m_ebounds.emax().MeV(), "[MeV] Upper bound of energy range");
1279  }
1280 
1281  // Write superpacket statistics
1282  hdu->card("NSPINP", m_num_superpackets, "Number of input superpackets");
1283  hdu->card("NSPUSE", m_num_used_superpackets, "Number of used superpackets");
1284  hdu->card("NSPSKP", m_num_skipped_superpackets, "Number of skipped superpackets");
1285 
1286  // Return
1287  return;
1288 }
1289 
1290 
1291 /***********************************************************************//**
1292  * @brief Compute DRG geometry factor
1293  *
1294  * @param[in] tjd TJD for module status
1295  * @param[in] theta Zenith angle in COMPTEL coordinates (deg).
1296  * @param[in] phi Azimuth angle in COMPTEL coordinates (deg).
1297  * @param[in] status D1 and D2 module status
1298  * @return Geometry factor.
1299  *
1300  * Computes the DRG geometry factor as function of zenith and azimuth angles
1301  * given in the COMPTEL coordinate system.
1302  ***************************************************************************/
1303 double GCOMDri::compute_geometry(const int& tjd,
1304  const double& theta,
1305  const double& phi,
1306  const GCOMStatus& status) const
1307 {
1308  // Set D1 module positions (from COM-RP-MPE-M10-123, Issue 1, Page 3-3)
1309  const double xd1[] = {0.0,-42.3,-26.0,26.0, 42.3,26.0,-26.0};
1310  const double yd1[] = {0.0, 0.0, 39.1,39.1, 0.0,-39.1,-39.1};
1311 
1312  // Set D2 module positions (from COM-RP-MPE-M10-123, Issue 1, Page 4-4)
1313  const double xd2[] = { 30.2, 0.0, -30.2, 45.3, 15.1, -15.1, -45.3,
1314  45.3, 15.1, -15.1, -45.3, 30.2, 0.0, -30.2};
1315  const double yd2[] = {-41.254,-41.254,-41.254,-15.1, -15.1, -15.1, -15.1,
1316  15.1, 15.1, 15.1 , 15.1, 41.254,41.254,41.254};
1317 
1318  // Set distance between D1 and D2 levels in cm (from COM-SP-MPE-M10-123,
1319  // Issue 1, page 8-5
1320  const double delz = 158.0;
1321 
1322  // Set D1 module radius in cm (from COM-TN-UNH-F70-051)
1323  //const double r1 = 13.8; // Defined globally
1324 
1325  // Set D2 module radius in cm (from SIM-AL-005)
1326  const double r2 = 14.085;
1327 
1328  // Derive some constants
1329  //const double r1sq = r1 * r1; // Defined globally
1330  const double r2sq = r2 * r2;
1331  const double drsq = r1sq - r2sq;
1332  const double norm_geo = 1.0 / (gammalib::pi * r1sq);
1333 
1334  // Initialise geometry factor
1335  double geometry = 0.0;
1336 
1337  // Precompute results
1338  double theta_rad = theta * gammalib::deg2rad;
1339  double phi_rad = phi * gammalib::deg2rad;
1340  double cosphi = std::cos(phi_rad);
1341  double sinphi = std::sin(phi_rad);
1342  double tantheta = std::tan(theta_rad);
1343 
1344  // Loop over all D1 modules
1345  for (int id1 = 0; id1 < 7; ++id1) {
1346 
1347  // Skip D1 module if it's off
1348  if (status.d1status(tjd, id1+1) != 1) {
1349  continue;
1350  }
1351 
1352  // Compute coordinates of D1 projected on D2 plane
1353  double xd1_prj = xd1[id1] - delz * tantheta * cosphi;
1354  double yd1_prj = yd1[id1] - delz * tantheta * sinphi;
1355 
1356  // Loop over all D2 modules
1357  for (int id2 = 0; id2 < 14; ++id2) {
1358 
1359  // Skip D2 module if it's off
1360  if (status.d2status(tjd, id2+1) != 1) {
1361  continue;
1362  }
1363 
1364  // Compute distance between center of D2 and projected centre
1365  // of D1
1366  double dx = xd2[id2] - xd1_prj;
1367  double dy = yd2[id2] - yd1_prj;
1368  double dsq = dx * dx + dy * dy;
1369  double d = std::sqrt(dsq);
1370 
1371  // If there is no overlap then skip the module
1372  if (d >= r1 + r2) {
1373  continue;
1374  }
1375 
1376  // If there is total overlap (within 0.1 cm) then add 1 to the
1377  // geometry factor
1378  else if (d <= r2 - r1 + 0.1) {
1379  geometry += 1.0;
1380  }
1381 
1382  // ... otherwise if there is a partial overlap then compute the
1383  // semiangle subtended by overlap sector of D2 and projected
1384  // D1 module
1385  else {
1386 
1387  // Cosine beta
1388  double d2 = 2.0 * d;
1389  double cbeta1 = (dsq + drsq) / (d2 * r1);
1390  double cbeta2 = (dsq - drsq) / (d2 * r2);
1391 
1392  // Sin beta
1393  double sbeta1 = std::sqrt(1.0 - cbeta1 * cbeta1);
1394  double sbeta2 = std::sqrt(1.0 - cbeta2 * cbeta2);
1395 
1396  // Beta
1397  double beta1 = std::acos(cbeta1);
1398  double beta2 = std::acos(cbeta2);
1399 
1400  // Projection
1401  geometry += (r1sq * (beta1 - sbeta1 * cbeta1) +
1402  r2sq * (beta2 - sbeta2 * cbeta2)) * norm_geo;
1403 
1404  } // endelse: there was partial overlap
1405 
1406  } // endfor: looped over D2 modules
1407 
1408  } // endfor: looped over D1 modules
1409 
1410  // Now divide by 7 since we want the probability of hitting a given
1411  // D1 module
1412  geometry /= 7.0;
1413 
1414  // Return geometry factor
1415  return geometry;
1416 }
1417 
1418 
1419 /***********************************************************************//**
1420  * @brief Compute DRM model for a point source
1421  *
1422  * @param[in] obs COMPTEL observation.
1423  * @param[in] model Sky model.
1424  *
1425  * @exception GException::invalid_argument
1426  * Model is not a point source model.
1427  *
1428  * Compute point source DRM model cube for a COMPTEL observation.
1429  ***************************************************************************/
1431  const GModelSky& model)
1432 {
1433  // Extract source direction from spatial model component
1434  const GModelSpatialPointSource* source =
1435  dynamic_cast<const GModelSpatialPointSource*>(model.spatial());
1436  if (source == NULL) {
1437  std::string msg = "Spatial component of model \""+model.name()+"\" "
1438  "is not a point source.";
1440  }
1441  const GSkyDir& srcDir = source->dir();
1442 
1443  // Extract response
1444  const GCOMResponse* rsp = obs.response();
1445 
1446  // Set all DRM cube bins to zero
1447  init_cube();
1448 
1449  // Get DRX value (units: cm^2 sec)
1450  double drx = obs.drx()(srcDir);
1451 
1452  // Get ontime
1453  double ontime = obs.ontime(); // sec
1454 
1455  // Get deadtime correction
1456  double deadc = obs.deadc(GTime());
1457 
1458  // Compute multiplicate IAQ normalisation
1459  double norm = drx * deadc / ontime;
1460 
1461  // Loop over all scatter angle pixels
1462  for (int index = 0; index < m_dri.npix(); ++index) {
1463 
1464  // Get sky direction of pixel
1465  GSkyDir obsDir = m_dri.inx2dir(index);
1466 
1467  // Compute angle between true photon arrival direction and scatter
1468  // direction (Chi,Psi)
1469  double phigeo = srcDir.dist_deg(obsDir);
1470 
1471  // Get Phigeo interpolation factor
1472  double phirat = phigeo / rsp->m_phigeo_bin_size; // 0.5 at bin centre
1473  int iphigeo = int(phirat); // index into which Phigeo falls
1474  double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre
1475 
1476  // Continue only if Phigeo is inside IAQ
1477  if (iphigeo < rsp->m_phigeo_bins) {
1478 
1479  // Initialise IAQ pixel index
1480  int i = iphigeo;
1481 
1482  // Interpolate towards left
1483  if (eps < 0.0) {
1484 
1485  // Not the first bin
1486  if (iphigeo > 0) {
1487  for (int iphibar = 0; iphibar < nphibar(); ++iphibar,
1488  i += rsp->m_phigeo_bins) {
1489  double iaq = obs.drg()(index, iphibar);
1490  iaq *= (1.0 + eps) * rsp->m_iaq[i] - eps * rsp->m_iaq[i-1];
1491  m_dri(index, iphibar) = iaq * norm;
1492  }
1493  }
1494  else {
1495  for (int iphibar = 0; iphibar < nphibar(); ++iphibar,
1496  i += rsp->m_phigeo_bins) {
1497  double iaq = obs.drg()(index, iphibar);
1498  iaq *= (1.0 - eps) * rsp->m_iaq[i] + eps * rsp->m_iaq[i+1];
1499  m_dri(index, iphibar) = iaq * norm;
1500  }
1501  }
1502 
1503  }
1504 
1505  // Interpolate towards right
1506  else {
1507 
1508  // Not the last IAQ bin
1509  if (iphigeo < rsp->m_phigeo_bins-1) {
1510  for (int iphibar = 0; iphibar < nphibar(); ++iphibar,
1511  i += rsp->m_phigeo_bins) {
1512  double iaq = obs.drg()(index, iphibar);
1513  iaq *= (1.0 - eps) * rsp->m_iaq[i] + eps * rsp->m_iaq[i+1];
1514  m_dri(index, iphibar) = iaq * norm;
1515  }
1516  }
1517  else {
1518  for (int iphibar = 0; iphibar < nphibar(); ++iphibar,
1519  i += rsp->m_phigeo_bins) {
1520  double iaq = obs.drg()(index, iphibar);
1521  iaq *= (1.0 + eps) * rsp->m_iaq[i] - eps * rsp->m_iaq[i-1];
1522  m_dri(index, iphibar) = iaq * norm;
1523  }
1524  }
1525 
1526  } // endfor: looped over all Compton scatter angles
1527 
1528  } // endif: Phigeo was inside IAQ
1529 
1530  } // endfor: looped over all scatter angle pixels
1531 
1532  // Return
1533  return;
1534 }
bool contains(const GTime &time) const
Checks whether Good Time Intervals contains time.
Definition: GGti.cpp:952
COMPTEL instrument status class.
Definition: GCOMStatus.hpp:49
void compute_drm_ptsrc(const GCOMObservation &obs, const GModelSky &model)
Compute DRM model for a point source.
Definition: GCOMDri.cpp:1430
COMPTEL instrument status class definition.
Abstract model class.
Definition: GModel.hpp:97
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:284
Abstract FITS image base class.
Definition: GFitsImage.hpp:43
void modcom(const int &modcom)
Set mini telescope.
const GCOMOads & oads(void) const
Return Orbit Aspect Data.
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
const int & ny(void) const
Returns number of pixels in y coordinate.
Definition: GSkyMap.hpp:366
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
int pix2inx(const GSkyPixel &pixel) const
Converts sky map pixel into pixel index.
Definition: GSkyMap.cpp:1397
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:414
Point source spatial model class interface definition.
double m_phigeo_bin_size
Phigeo binsize (deg)
#define G_COMPUTE_DRM
Definition: GCOMDri.cpp:52
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1276
GTime m_tstart
Selection start time.
Definition: GCOMDri.hpp:141
GCOMSelection m_selection
Selection parameters.
Definition: GCOMDri.hpp:148
const double & phimin(void) const
Return minimum Compton scatter angle of DRI cube.
Definition: GCOMDri.hpp:338
GSkyMap extract(const int &map, const int &nmaps=1) const
Extract maps into a new sky map object.
Definition: GSkyMap.cpp:2036
const double pi
Definition: GMath.hpp:35
std::vector< double > m_iaq
IAQ array.
const double superpacket_duration
Definition: GCOMDri.cpp:67
int size(void) const
Return number of bins.
Definition: GCOMDri.hpp:197
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
Abstract FITS extension base class.
Definition: GFitsHDU.hpp:51
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
GSkyMap m_dri
Data cube.
Definition: GCOMDri.hpp:134
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
COMPTEL Data Space class definition.
COMPTEL event list class definition.
double m_phibin
Phibar binsize (deg)
Definition: GCOMDri.hpp:138
void compute_drg(const GCOMObservation &obs, const GCOMSelection &select=GCOMSelection(), const double &zetamin=5.0)
Compute geometry cube.
Definition: GCOMDri.cpp:506
virtual double ontime(void) const
Return ontime.
double compute_geometry(const int &tjd, const double &theta, const double &phi, const GCOMStatus &status) const
Compute DRG geometry factor.
Definition: GCOMDri.cpp:1303
void write_attributes(GFitsHDU *hdu) const
Write DRI attributes into FITS HDU.
Definition: GCOMDri.cpp:1245
COMPTEL observation class interface definition.
const float & gcel(void) const
Return Geocentre zenith angle.
Definition: GCOMOad.hpp:276
void clear(void)
Clear time.
Definition: GTime.cpp:251
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
Abstract spatial model base class interface definition.
int m_num_superpackets
Number of superpackets.
Definition: GCOMDri.hpp:143
double crval(const int &inx) const
Return value of reference pixel.
Definition: GWcs.cpp:877
Interface for the COMPTEL instrument response function.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
double phi(const GSkyDir &sky) const
Return azimuth angle of sky direction in COMPTEL coordinates.
Definition: GCOMOad.hpp:410
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
const GTime & tstop(void) const
Return stop time of superpacket.
Definition: GCOMOad.hpp:160
int size(void) const
Return number of Orbit Aspect Data in container.
Definition: GCOMOads.hpp:141
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
GCOMDri(void)
Void constructor.
Definition: GCOMDri.cpp:84
GTime com_time(const int &tjd, const int &tics)
Convert TJD and COMPTEL ticks in GTime object.
COMPTEL selection set class.
COMPTEL Good Time Intervals class definition.
std::string print(const GChatter &chatter=NORMAL) const
Print Good Time Intervals.
Definition: GGti.cpp:988
const float & georad(void) const
Return apparent radius of Earth.
Definition: GCOMOad.hpp:305
virtual double deadc(const GTime &time=GTime()) const
Return deadtime correction factor.
void read(const GFitsHDU &hdu)
Read skymap from FITS HDU.
Definition: GSkyMap.cpp:2434
virtual GClassCode code(void) const =0
COMPTEL Orbit Aspect Data class.
Definition: GCOMOad.hpp:48
Implementation of support function used by COMPTEL classes.
int size(void) const
Return number of Good Time Intervals.
Definition: GGti.hpp:153
Abstract FITS image base class definition.
const GTime & tstart(void) const
Return start time of superpacket.
Definition: GCOMOad.hpp:129
COMPTEL selection set class definition.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1267
std::string print(const GChatter &chatter=NORMAL) const
Print sky map.
Definition: GSkyMap.cpp:2557
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
void copy_members(const GCOMDri &dri)
Copy class members.
Definition: GCOMDri.cpp:1046
Sky map pixel class.
Definition: GSkyPixel.hpp:74
void reference(const GTimeReference &ref)
Set time reference for Good Time Intervals.
Definition: GGti.hpp:256
const GCOMTim & tim(void) const
Return COMPTEL Good Time Intervals.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
double crpix(const int &inx) const
Return reference pixel.
Definition: GWcs.cpp:900
void clear(void)
Clear Good Time Intervals.
Definition: GGti.cpp:237
virtual GCOMDri * clone(void) const
Clone COMPTEL Data Space.
Definition: GCOMDri.cpp:241
const double d1_area
Definition: GCOMDri.cpp:72
double cone_content(const GSkyDir &dir, const double &armmin, const double &armmax) const
Compute content in cone.
Definition: GCOMDri.cpp:794
virtual int size(void) const
Return number of events in list.
bool contains(const GTime &time) const
Check if time is comprised in the Good Time Intervals.
Definition: GCOMTim.hpp:113
Energy boundaries container class.
Definition: GEbounds.hpp:60
const int & nmaps(void) const
Returns number of maps.
Definition: GSkyMap.hpp:378
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1321
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1289
void init_cube(void)
Initialise DRI cube.
Definition: GCOMDri.cpp:1087
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:253
GFitsHDU * write(GFits &file, const std::string &extname="") const
Write sky map into FITS file.
Definition: GSkyMap.cpp:2467
Filename class.
Definition: GFilename.hpp:62
GTime m_tstop
Selection stop time.
Definition: GCOMDri.hpp:142
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
void write(GFits &fits, const std::string &extname="") const
Write COMPTEL Data Space into FITS image.
Definition: GCOMDri.cpp:915
void compute_drm(const GCOMObservation &obs, const GModel &model)
Compute DRM model.
Definition: GCOMDri.cpp:741
void com_wcs_mer2car(GSkyMap &map)
Changes Mercator&#39;s projection to cartesian projection.
Definition: GCOMSupport.cpp:58
const GSkyMap & map(void) const
Return DRI sky map.
Definition: GCOMDri.hpp:245
int com_tjd(const GTime &time)
Convert GTime in COMPTEL TJD.
double cdelt(const int &inx) const
Return pixel size.
Definition: GWcs.cpp:923
void init_statistics(void) const
Initialise selection statistics.
bool contains(const GSkyDir &dir) const
Checks if sky direction falls in map.
Definition: GSkyMap.cpp:1879
int m_num_skipped_superpackets
Number of skipped superpackets.
Definition: GCOMDri.hpp:145
std::string print(const GChatter &chatter=NORMAL) const
Print energy boundaries.
Definition: GEbounds.cpp:1244
GChatter
Definition: GTypemaps.hpp:33
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
void free_members(void)
Delete class members.
Definition: GCOMDri.cpp:1075
void save(const GFilename &filename, const bool &clobber=false) const
Save COMPTEL Data Space into DRI FITS file.
Definition: GCOMDri.cpp:866
void clear(void)
Clear instance.
Definition: GSkyMap.cpp:1090
Good Time Interval class.
Definition: GGti.hpp:62
const double r1
Definition: GCOMDri.cpp:68
int d1status(const int &tjd, const int &module) const
Return D1 module status.
Definition: GCOMStatus.cpp:188
#define G_COMPUTE_DRE
Definition: GCOMDri.cpp:50
void radec_deg(const double &ra, const double &dec)
Set equatorial sky direction (degrees)
Definition: GSkyDir.cpp:224
const GTime & tstop(void) const
Returns latest stop time in Good Time Intervals.
Definition: GGti.hpp:206
virtual ~GCOMDri(void)
Destructor.
Definition: GCOMDri.cpp:171
const GTime & tstart(void) const
Returns earliest start time in Good Time Intervals.
Definition: GGti.hpp:193
void clear(void)
Clear energy boundaries.
Definition: GEbounds.cpp:269
Point source spatial model.
int nphibar(void) const
Return number of Phibar bins.
Definition: GCOMDri.hpp:233
GCOMDri & operator=(const GCOMDri &dri)
Assignment operator.
Definition: GCOMDri.cpp:193
bool use_event(const GCOMEventAtom &event) const
Check if event should be used.
const std::string extname_dri
Definition: GCOMDri.hpp:51
virtual std::string type(void) const =0
int d2status(const int &tjd, const int &module) const
Return D2 module status.
Definition: GCOMStatus.cpp:217
virtual void response(const GResponse &rsp)
Set response function.
GFitsImage * image(const int &extno)
Get pointer to image HDU.
Definition: GFits.cpp:360
const GSkyMap & drx(void) const
Return exposure.
const double & secs(void) const
Return time in seconds in native reference (TT)
Definition: GTime.hpp:153
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
int com_tics(const GTime &time)
Convert GTime in COMPTEL tics.
GSkyPixel dir2pix(const GSkyDir &dir) const
Returns sky map pixel for a given sky direction.
Definition: GSkyMap.cpp:1462
GSkyDir dir(void) const
Return position of point source.
void init_statistics(void)
Initialise computation statistics.
Definition: GCOMDri.cpp:1102
Sky model class interface definition.
Sky model class.
Definition: GModelSky.hpp:120
virtual void clear(void)
Clear COMPTEL Data Space.
Definition: GCOMDri.cpp:223
const int & tjd(void) const
Return Truncated Julian Days of Orbit Aspect Record.
Definition: GCOMOad.hpp:189
void read(const GFitsImage &image)
Read COMPTEL Data Space from DRI FITS image.
Definition: GCOMDri.cpp:890
void init_members(void)
Initialise class members.
Definition: GCOMDri.cpp:1019
GGti m_gti
Good Time Intervals of data cube.
Definition: GCOMDri.hpp:136
Interface class for COMPTEL observations.
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
int m_phigeo_bins
Number of Phigeo bins.
virtual void clear(void)
Clear COMPTEL selection set.
void load(const GFilename &filename)
Load COMPTEL Data Space from DRI FITS file.
Definition: GCOMDri.cpp:841
GEbounds m_ebounds
Energy boundaries of data cube.
Definition: GCOMDri.hpp:135
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
virtual GEvents * events(void)
Return events.
void compute_dre(const GCOMObservation &obs, const GCOMSelection &select=GCOMSelection(), const double &zetamin=5.0)
Compute event cube.
Definition: GCOMDri.cpp:259
COMPTEL event list class.
const int & nx(void) const
Returns number of pixels in x coordinate.
Definition: GSkyMap.hpp:350
COMPTEL Data Space class.
Definition: GCOMDri.hpp:60
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:246
int m_num_used_superpackets
Number of used superpackets.
Definition: GCOMDri.hpp:144
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
double m_zetamin
Minimum zeta angle.
Definition: GCOMDri.hpp:149
Abstract spatial model base class.
void compute_drx(const GCOMObservation &obs)
Compute DRX exposure map.
Definition: GCOMDri.cpp:640
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
COMPTEL Good Time Intervals class.
Definition: GCOMTim.hpp:50
bool use_superpacket(const GCOMOad &oad, const GCOMTim &tim)
Check if superpacket should be used.
Definition: GCOMDri.cpp:1128
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1033
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:334
const double & phibin(void) const
Return Compton scatter angle bin of DRI cube.
Definition: GCOMDri.hpp:350
void read_attributes(const GFitsHDU *hdu)
Read DRI attributes from FITS HDU.
Definition: GCOMDri.cpp:1179
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
virtual std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL Data Space.
Definition: GCOMDri.cpp:936
std::string utc(const int &precision=0) const
Return time as string in UTC time system.
Definition: GTime.cpp:464
double m_phimin
Phibar minimum (deg)
Definition: GCOMDri.hpp:137
const double r1sq
Definition: GCOMDri.cpp:71
Abstract world coordinate system base class definition.
double theta(const GSkyDir &sky) const
Return zenith angle of sky direction in COMPTEL coordinates.
Definition: GCOMOad.hpp:395
#define G_COMPUTE_DRE_PTSRC
Definition: GCOMDri.cpp:53
COMPTEL Orbit Aspect Data container class definition.
COMPTEL Orbit Aspect Data class definition.
Mathematical function definitions.
const float & gcaz(void) const
Return Geocentre azimuth angle.
Definition: GCOMOad.hpp:247
const GSkyMap & drg(void) const
Return geometry factors.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
COMPTEL event atom class.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415
std::string m_name
Data cube name.
Definition: GCOMDri.hpp:133