GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMEventList.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMEventList.cpp - COMPTEL event list class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-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 GCOMEventList.hpp
23  * @brief COMPTEL event list class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <typeinfo>
32 #include "GMath.hpp"
33 #include "GFits.hpp"
34 #include "GException.hpp"
35 #include "GCOMTools.hpp"
36 #include "GCOMSupport.hpp"
37 #include "GCOMEventList.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_OPERATOR "GCOMEventList::operator[](int&)"
41 #define G_ROI "GCOMEventList::roi(GRoi&)"
42 
43 /* __ Macros _____________________________________________________________ */
44 
45 /* __ Coding definitions _________________________________________________ */
46 
47 /* __ Debug definitions __________________________________________________ */
48 //#define G_DEBUG_TOFCOR //!< Debug TOF correction
49 
50 
51 /*==========================================================================
52  = =
53  = Constructors/destructors =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Void constructor
59  *
60  * Creates an empty COMPTEL event list.
61  ***************************************************************************/
63 {
64  // Initialise class members for clean destruction
65  init_members();
66 
67  // Return
68  return;
69 }
70 
71 
72 /***********************************************************************//**
73  * @brief File name constructor
74  *
75  * @param[in] filename COMPTEL event list filename.
76  *
77  * Construct COMPTEL event list object by loading the events from a
78  * FITS file.
79  ***************************************************************************/
81 {
82  // Initialise members
83  init_members();
84 
85  // Load event list
86  load(filename);
87 
88  // Return
89  return;
90 }
91 
92 
93 /***********************************************************************//**
94  * @brief Copy constructor
95  *
96  * @param[in] list COMPTEL event list.
97  ***************************************************************************/
99 {
100  // Initialise members
101  init_members();
102 
103  // Copy members
104  copy_members(list);
105 
106  // Return
107  return;
108 }
109 
110 
111 /***********************************************************************//**
112  * @brief Destructor
113  ***************************************************************************/
115 {
116  // Free members
117  free_members();
118 
119  // Return
120  return;
121 }
122 
123 
124 /*==========================================================================
125  = =
126  = Operators =
127  = =
128  ==========================================================================*/
129 
130 /***********************************************************************//**
131  * @brief Assignment operator
132  *
133  * @param[in] list COMPTEL event list.
134  * @return COMPTEL event list.
135  ***************************************************************************/
137 {
138  // Execute only if object is not identical
139  if (this != &list) {
140 
141  // Copy base class members
142  this->GEventList::operator=(list);
143 
144  // Free members
145  free_members();
146 
147  // Initialise members
148  init_members();
149 
150  // Copy members
151  copy_members(list);
152 
153  } // endif: object was not identical
154 
155  // Return this object
156  return *this;
157 }
158 
159 
160 /***********************************************************************//**
161  * @brief COMPTEL event atom access operator
162  *
163  * @param[in] index Event index [0,...,size()-1].
164  * @return Pointer to COMPTEL event atom.
165  *
166  * @exception GException::out_of_range
167  * Event index outside valid range.
168  *
169  * Returns pointer to a COMPTEL event atom.
170  ***************************************************************************/
172 {
173  // Optionally check if the index is valid
174  #if defined(G_RANGE_CHECK)
175  if (index < 0 || index >= size()) {
176  throw GException::out_of_range(G_OPERATOR, "Event index", index, size());
177  }
178  #endif
179 
180  // Return pointer
181  return (&(m_events[index]));
182 }
183 
184 
185 /***********************************************************************//**
186  * @brief COMPTEL event atom access operator
187  *
188  * @param[in] index Event index [0,...,size()-1].
189  * @return Pointer to COMPTEL event atom.
190  *
191  * @exception GException::out_of_range
192  * Event index outside valid range.
193  *
194  * Returns pointer to a COMPTEL event atom.
195  ***************************************************************************/
196 const GCOMEventAtom* GCOMEventList::operator[](const int& index) const
197 {
198  // Optionally check if the index is valid
199  #if defined(G_RANGE_CHECK)
200  if (index < 0 || index >= size()) {
201  throw GException::out_of_range(G_OPERATOR, "Event index", index, size());
202  }
203  #endif
204 
205  // Return pointer
206  return (&(m_events[index]));
207 }
208 
209 
210 /*==========================================================================
211  = =
212  = Public methods =
213  = =
214  ==========================================================================*/
215 
216 /***********************************************************************//**
217  * @brief Clear COMPTEL event list
218  *
219  * Clears COMPTEL event list by resetting all class members to an
220  * initial state. Any information that was present before will be lost.
221  ***************************************************************************/
223 {
224  // Free class members (base and derived classes, derived class first)
225  free_members();
226  this->GEventList::free_members();
227  this->GEvents::free_members();
228 
229  // Initialise members
230  this->GEvents::init_members();
231  this->GEventList::init_members();
232  init_members();
233 
234  // Return
235  return;
236 }
237 
238 
239 /***********************************************************************//**
240  * @brief Clone event list
241  *
242  * @return Pointer to deep copy of COMPTEL event list.
243  ***************************************************************************/
245 {
246  return new GCOMEventList(*this);
247 }
248 
249 
250 /***********************************************************************//**
251  * @brief Load COMPTEL events from FITS file
252  *
253  * @param[in] filename COMPTEL event list FITS file name.
254  *
255  * Loads COMPTEL events from a FITS file into the event list.
256  ***************************************************************************/
257 void GCOMEventList::load(const GFilename& filename)
258 {
259  // Open FITS file
260  GFits fits(filename);
261 
262  // Read event list from FITS file
263  read(fits);
264 
265  // Close FITS file
266  fits.close();
267 
268  // Return
269  return;
270 }
271 
272 
273 /***********************************************************************//**
274  * @brief Save COMPTEL events
275  *
276  * @param[in] filename COMPTEL event list FITS file name.
277  * @param[in] clobber Overwrite existing FITS file?
278  *
279  * Save COMPTEL events from a FITS file into the event list.
280  ***************************************************************************/
281 void GCOMEventList::save(const GFilename& filename,
282  const bool& clobber) const
283 {
284  // Open FITS file
285  GFits fits;
286 
287  // Write events
288  write(fits);
289 
290  // Save events
291  fits.saveto(filename, clobber);
292 
293  // Close FITS file
294  fits.close();
295 
296  // Return
297  return;
298 }
299 
300 
301 /***********************************************************************//**
302  * @brief Read COMPTEL events from FITS file.
303  *
304  * @param[in] file FITS file.
305  *
306  * Read the COMPTEL event list from a FITS file object.
307  ***************************************************************************/
308 void GCOMEventList::read(const GFits& file)
309 {
310  // Clear object
311  clear();
312 
313  // Get HDU (pointer is always valid)
314  const GFitsTable& hdu = *file.table(1);
315 
316  // Read event data
317  read_events(hdu);
318 
319  // If there are events then get start and stop time from the first
320  // and last event and append them as single time interval to the GTI
321  if (size() > 0) {
322  GTime start = m_events[0].time();
323  GTime stop = m_events[size()-1].time();
324  m_gti.append(start, stop);
325  }
326 
327  // Return
328  return;
329 }
330 
331 
332 /***********************************************************************//**
333  * @brief Write COMPTEL event list into FITS file.
334  *
335  * @param[in] file FITS file.
336  *
337  * Write the COMPTEL event list into FITS file.
338  *
339  * @todo Implement method.
340  ***************************************************************************/
341 void GCOMEventList::write(GFits& file) const
342 {
343  // TODO: You need to implement an interface to the FITS file that so
344  // that you can save your events.
345 
346  // Return
347  return;
348 }
349 
350 
351 /***********************************************************************//**
352  * @brief Set region of interest
353  *
354  * @param[in] roi Region of interest.
355  *
356  * @exception GException::invalid_argument
357  * Specified RoI is not a COMPTEL RoI.
358  *
359  * Sets the region of interest for the observation.
360  ***************************************************************************/
361 void GCOMEventList::roi(const GRoi& roi)
362 {
363  // Get pointer on COMPTEL region of interest
364  const GCOMRoi* comroi = dynamic_cast<const GCOMRoi*>(&roi);
365  if (comroi == NULL) {
366  std::string cls = std::string(typeid(&roi).name());
367  std::string msg = "Region of interest of type \""+cls+"\" is "
368  "not a COMPTEL RoI. Please specify a "
369  "COMPTEL RoI as argument.";
371  }
372 
373  // Set RoI
374  m_roi = *comroi;
375 
376  // Return
377  return;
378 }
379 
380 
381 /***********************************************************************//**
382  * @brief Append event to event list
383  *
384  * @param[in] event Event.
385  *
386  * Appends an event to the end of the event list.
387  ***************************************************************************/
389 {
390  // Append event
391  m_events.push_back(event);
392 
393  // Return
394  return;
395 }
396 
397 
398 /***********************************************************************//**
399  * @brief Remove events from event list
400  *
401  * @param[in] index Index from which on events should be removed.
402  * @param[in] number Number of event to remove (default: 1).
403  *
404  * Removes events from the event list. This method does nothing if @p index
405  * points beyond the event list. The method does also gently accept
406  * @p number arguments where @p index + @p number reach beyond the event
407  * list. In that case, all events from event @p index on will be removed.
408  ***************************************************************************/
409 void GCOMEventList::remove(const int& index, const int& number)
410 {
411  // Continue only if index is valid
412  if (index < size()) {
413 
414  // Determine number of elements to remove
415  int n_remove = (index + number > size()) ? size() - index : number;
416 
417  // Remove events
418  m_events.erase(m_events.begin() + index,
419  m_events.begin() + index + n_remove);
420 
421  } // endif: index was valid
422 
423  // Return
424  return;
425 }
426 
427 
428 /***********************************************************************//**
429  * @brief Print COMPTEL event list information
430  *
431  * @param[in] chatter Chattiness.
432  * @return String containing COMPTEL event list information.
433  ***************************************************************************/
434 std::string GCOMEventList::print(const GChatter& chatter) const
435 {
436  // Initialise result string
437  std::string result;
438 
439  // Continue only if chatter is not silent
440  if (chatter != SILENT) {
441 
442  // Append header
443  result.append("=== GCOMEventList ===");
444 
445  // Append information
446  result.append("\n"+gammalib::parformat("Number of events") +
447  gammalib::str(number()));
448 
449  // Append time information
450  if (gti().size() > 0) {
451  result.append("\n"+gammalib::parformat("TJD interval"));
452  result.append(gammalib::str(gammalib::com_tjd(tstart()))+":");
453  result.append(gammalib::str(gammalib::com_tics(tstart()))+" - ");
454  result.append(gammalib::str(gammalib::com_tjd(tstop()))+":");
455  result.append(gammalib::str(gammalib::com_tics(tstop())));
456  result.append("\n"+gammalib::parformat("MJD interval"));
457  result.append(gammalib::str(tstart().mjd())+" - ");
458  result.append(gammalib::str(tstop().mjd())+" days");
459  result.append("\n"+gammalib::parformat("UTC interval"));
460  result.append(tstart().utc()+" - ");
461  result.append(tstop().utc());
462  }
463  else {
464  result.append("\n"+gammalib::parformat("MJD interval"));
465  result.append("not defined");
466  }
467 
468  // Append energy information
469  result.append("\n"+gammalib::parformat("Energy interval"));
470  if (ebounds().size() > 0) {
471  result.append(gammalib::str(ebounds().emin().MeV()));
472  result.append(" - ");
473  result.append(gammalib::str(ebounds().emax().MeV())+" MeV");
474  }
475  else {
476  result.append("not defined");
477  }
478 
479  } // endif: chatter was not silent
480 
481  // Return result
482  return result;
483 }
484 
485 
486 /*==========================================================================
487  = =
488  = Private methods =
489  = =
490  ==========================================================================*/
491 
492 /***********************************************************************//**
493  * @brief Initialise class members
494  ***************************************************************************/
496 {
497  // Initialise members
498  m_roi.clear();
499  m_events.clear();
500 
501  // Return
502  return;
503 }
504 
505 
506 /***********************************************************************//**
507  * @brief Copy class members
508  *
509  * @param[in] list COMPTEL event list.
510  ***************************************************************************/
512 {
513  // Copy members
514  m_roi = list.m_roi;
515  m_events = list.m_events;
516 
517  // Return
518  return;
519 }
520 
521 
522 /***********************************************************************//**
523  * @brief Delete class members
524  ***************************************************************************/
526 {
527  // Return
528  return;
529 }
530 
531 
532 /***********************************************************************//**
533  * @brief Read COMPTEL events from FITS table
534  *
535  * @param[in] table FITS table.
536  *
537  * Reads COMPTEL events from a FITS table. The method expects the data in the
538  * format that is provided by the HEASARC archive.
539  *
540  * The method also sets the energy boundaries and Region of Interest of the
541  * data set as this information is not provided in the FITS table header.
542  * For that purpose it scans all data and determines the range that is
543  * covered by the data.
544  *
545  * The method assumes that no events exist so far.
546  ***************************************************************************/
548 {
549  // Extract number of events in FITS file
550  int num = table.nrows();
551 
552  // If there are events then load them
553  if (num > 0) {
554 
555  // Use pointing direction as roi centre
556  GSkyDir centre;
557  centre.lb_deg(table.real("GLON_SCZ"), table.real("GLAT_SCZ"));
558 
559  // Get data set representation version
560  int verno = table.integer("DSD_REP");
561  #if defined(G_DEBUG_TOFCOR)
562  std::cout << "GCOMEventList::read_events: DSD_REP = ";
563  std::cout << verno << std::endl;
564  #endif
565 
566  // Initialise boundaries
567  double radius_max = 0.0;
568  double phibar_min = 0.0;
569  double phibar_max = 0.0;
570 
571  // Reserve data
572  m_events.reserve(num);
573 
574  // Get column pointers
575  const GFitsTableCol* ptr_tjd = table["TJD"]; // days
576  const GFitsTableCol* ptr_tics = table["TICS"]; // ticks
577  const GFitsTableCol* ptr_glon_scat = table["GLON_SCAT"]; // rad
578  const GFitsTableCol* ptr_glat_scat = table["GLAT_SCAT"]; // rad
579  const GFitsTableCol* ptr_phi_scat = table["AZIMUTH_SCAT"]; // rad
580  const GFitsTableCol* ptr_theta_scat = table["ZENITH_SCAT"]; // rad
581  const GFitsTableCol* ptr_phibar = table["PHIBAR"]; // rad
582  const GFitsTableCol* ptr_eha = table["EARTH_HORIZON"]; // rad
583  const GFitsTableCol* ptr_e1 = table["E_D1"]; // keV
584  const GFitsTableCol* ptr_e2 = table["E_D2"]; // keV
585  const GFitsTableCol* ptr_x2 = table["X_D2"]; // mm
586  const GFitsTableCol* ptr_y2 = table["Y_D2"]; // mm
587  const GFitsTableCol* ptr_psd = table["PSD"]; // channel
588  const GFitsTableCol* ptr_tof = table["TOF"]; // channel
589  const GFitsTableCol* ptr_modcom = table["MODCOM"]; // id
590  const GFitsTableCol* ptr_reflag = table["RC_REFLAG"]; //
591  const GFitsTableCol* ptr_veto = table["RC_VETO"]; //
592 
593  // Initialise boundaries
594  GEnergy emin;
595  GEnergy emax;
596 
597  // Debug TOFCOR: initialise flag for first notice
598  #if defined(G_DEBUG_TOFCOR)
599  bool noticed = false;
600  #endif
601 
602  // Copy data from columns into GCOMEventAtom objects
603  for (int i = 0; i < num; ++i) {
604 
605  // Allocate event
606  GCOMEventAtom event;
607 
608  // Set instrument direction. Note that contrary to what is
609  // specified in the FITS file, angles are provided in radians.
610  // Also, GLON and GLAT are inverted.
611  GCOMInstDir inst_dir;
612  GSkyDir sky_dir;
613  sky_dir.lb(ptr_glat_scat->real(i), ptr_glon_scat->real(i));
614  inst_dir.dir(sky_dir);
615  inst_dir.phibar(ptr_phibar->real(i) * gammalib::rad2deg);
616 
617  // Set total energy
618  GEnergy etot;
619  etot.keV(ptr_e1->real(i) + ptr_e2->real(i));
620 
621  // Set phibar value (deg)
622  double phibar = ptr_phibar->real(i) * gammalib::rad2deg;
623 
624  // Set X and Y values of D2 interaction
625  double x_d2 = double(ptr_x2->integer(i))*0.03125;
626  double y_d2 = double(ptr_y2->integer(i))*0.03125;
627 
628  // Set PSD and TOF values. Note that the values are read unscaled,
629  // hence we need to divide explicitly by the 1/TSCAL factor. This
630  // emulates the behaviour of the evpdal13.pevpsr.f function.
631  double psd = double(ptr_psd->integer(i))/128.0;
632  double tof = double(ptr_tof->integer(i))/128.0;
633 
634  // Correct TOF if data representation flag is less than 3 and
635  // rejection flag is >=4 (see dalaaa19.pevpsr.f)
636  if (verno < 3) {
637  if (ptr_reflag->integer(i) >= 4) {
638  #if defined(G_DEBUG_TOFCOR)
639  double tof_raw = tof;
640  #endif
641  tof = tofcor(ptr_e1->real(i), ptr_e2->real(i), tof);
642  #if defined(G_DEBUG_TOFCOR)
643  if (!noticed) {
644  std::cout << "GCOMEventList::read_events: ";
645  std::cout << "apply TOF correction (";
646  std::cout << tof_raw;
647  std::cout << " => ";
648  std::cout << tof;
649  std::cout << ")" << std::endl;
650  noticed = true;
651  }
652  #endif
653  }
654  }
655 
656  // Set event information
657  event.time(ptr_tjd->integer(i), ptr_tics->integer(i));
658  event.energy(etot);
659  event.dir(inst_dir);
660  event.phi(ptr_phi_scat->real(i) * gammalib::rad2deg); // rad -> deg
661  event.theta(ptr_theta_scat->real(i) * gammalib::rad2deg); // rad -> deg
662  event.phibar(phibar); // rad -> deg
663  event.eha(ptr_eha->real(i) * gammalib::rad2deg); // rad -> deg
664  event.e1(ptr_e1->real(i) * 1.0e-3); // keV -> MeV
665  event.e2(ptr_e2->real(i) * 1.0e-3); // keV -> MeV
666  event.x_d2(x_d2);
667  event.y_d2(y_d2);
668  event.psd(psd);
669  event.tof(tof);
670  event.modcom(ptr_modcom->integer(i));
671  event.reflag(ptr_reflag->integer(i));
672  event.veto(ptr_veto->integer(i));
673 
674  // Append event
675  m_events.push_back(event);
676 
677  // Compute distance from pointing direction
678  double radius = centre.dist_deg(sky_dir);
679 
680  // Update boundaries
681  if (i == 0) {
682  emin = etot;
683  emax = etot;
684  phibar_min = phibar;
685  phibar_max = phibar;
686  radius_max = radius;
687  }
688  else {
689  if (etot < emin) {
690  emin = etot;
691  }
692  if (etot > emax) {
693  emax = etot;
694  }
695  if (phibar < phibar_min) {
696  phibar_min = phibar;
697  }
698  if (phibar > phibar_max) {
699  phibar_max = phibar;
700  }
701  if (radius > radius_max) {
702  radius_max = radius;
703  }
704  }
705 
706  } // endfor: looped over all events
707 
708  // Set region of interest
709  double phibar = 0.5 * (phibar_min + phibar_max);
710  m_roi = GCOMRoi(GCOMInstDir(centre, phibar), radius_max,
711  phibar_min, phibar_max);
712 
713  // Set energy boundaries
714  m_ebounds = GEbounds(emin, emax);
715 
716  } // endif: there were events
717 
718  // Return
719  return;
720 }
721 
722 
723 /***********************************************************************//**
724  * @brief Compute TOF correction
725  *
726  * @param[in] d1e D1 energy deposit (keV).
727  * @param[in] d2e D2 energy deposit (keV).
728  * @param[in] tof TOF value from EVP dataset (version < 3).
729  *
730  * Computes the TOF correction for EVP data sets with representation version
731  * < 3.
732  ***************************************************************************/
733 double GCOMEventList::tofcor(const double& d1e, const double& d2e, double tof) const
734 {
735  // Begin and end of D1E-range for which fit is valid
736  const double d1sta = 50.0;
737  const double d1end = 20000.0;
738 
739  // Begin and end of D2E-range for which fit is valid
740  const double d2sta = 600.0;
741  const double d2end = 30000.0;
742 
743  // Degree of polynomial fit to ToF(D1E) for which fit is valid
744  const int ndegd1 = 6;
745 
746  // D1-energy 'breaking' the D1E-fit range in 2 intervals
747  const double d1ebrk = 2250.0;
748 
749  // Degree of polynomial fit to ToF(D2E) for which fit is valid
750  const int ndegd2 = 4;
751 
752  // D2-energy 'breaking' the D2E-fit range in 3 intervals
753  const double d2ebrk1 = 1400.0;
754  const double d2ebrk2 = 5500.0;
755 
756  // Reference ToF: the ToF corrections are calculated with respect to this TOF
757  const double tofref = 118.30;
758 
759  // ToF the corrected peak will be centered at
760  const double tofcen = 120.0;
761 
762  // Polynomial coefficients of ToF(D1E) for interval 1: D1E < D1EBRK
763  const double d11cof[] = {111.74858,
764  28.280247,
765  -45.024305,
766  35.183210,
767  -14.639463,
768  3.1342536,
769  -0.2711735};
770 
771  // Polynomial coefficients of ToF(D1E) for interval 2: D1E > D1EBRK
772  const double d12cof[] = {116.25374,
773  0.500407092,
774  0.38182720,
775  -0.080145513,
776  0.0065569790,
777  -0.00024650067,
778  3.5077240e-6};
779 
780  // Polynomial coefficients of ToF(D2E) for interval 1: D2E < D2EBRK1
781  const double d21cof[] = { 181.77024,
782  -252.41070,
783  371.09898,
784  -232.83985,
785  52.918785};
786 
787  // Polynomial coefficients of ToF(D2E) for interval 2: D2EBRK1 < D2E < D2EBRK2
788  const double d22cof[] = { 120.91608,
789  -0.15048490,
790  -0.45526025,
791  0.11710009,
792  -0.0082172427};
793 
794  // Polynomial coefficients of ToF(D2E) for interval 3: D2EBRK2 < D2E
795  const double d23cof[] = { 119.24278,
796  -0.43134699,
797  0.060183080,
798  -0.0026847790,
799  3.7720986e-5};
800 
801  // Make sure that D1 and D2 energies are within fit range
802  double d1ener = d1e;
803  double d2ener = d2e;
804  if (d1ener < d1sta) {
805  d1ener = d1sta;
806  }
807  if (d1ener > d1end) {
808  d1ener = d1end;
809  }
810  if (d2ener < d2sta) {
811  d2ener = d2sta;
812  }
813  if (d2ener > d2end) {
814  d2ener = d2end;
815  }
816 
817  // Compute D1 correction
818  double d1ecor = -tofref;
819  double d1emev = d1ener * 1.0e-3;
820  double d1edum = 1.0;
821  if (d1ener <= d1ebrk) {
822  for (int i = 0; i <= ndegd1; ++i) {
823  d1ecor += d11cof[i] * d1edum;
824  d1edum *= d1emev;
825  }
826  }
827  else {
828  for (int i = 0; i <= ndegd1; ++i) {
829  d1ecor += d12cof[i] * d1edum;
830  d1edum *= d1emev;
831  }
832  }
833 
834  // Compute D2 correction
835  double d2ecor = -tofref;
836  double d2emev = d2ener * 1.0e-3;
837  double d2edum = 1.0;
838  if (d2ener <= d2ebrk1) {
839  for (int i = 0; i <= ndegd2; ++i) {
840  d2ecor += d21cof[i] * d2edum;
841  d2edum *= d2emev;
842  }
843  }
844  else if (d2ener > d2ebrk2) {
845  for (int i = 0; i <= ndegd2; ++i) {
846  d2ecor += d23cof[i] * d2edum;
847  d2edum *= d2emev;
848  }
849  }
850  else {
851  for (int i = 0; i <= ndegd2; ++i) {
852  d2ecor += d22cof[i] * d2edum;
853  d2edum *= d2emev;
854  }
855  }
856 
857  // Compute total correction
858  tof += tofcen - (tofref + d1ecor + d2ecor);
859 
860  // Return corrected TOF
861  return tof;
862 }
void phibar(const double &phibar)
Set event Compton scatter angle.
#define G_OPERATOR
void free_members(void)
Delete class members.
#define G_ROI
const GGti & gti(void) const
Return Good Time Intervals.
Definition: GEvents.hpp:134
double dist_deg(const GSkyDir &dir) const
Compute angular distance between sky directions in degrees.
Definition: GSkyDir.hpp:286
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:482
std::string number(const std::string &noun, const int &number)
Convert singular noun into number noun.
Definition: GTools.cpp:1167
Definition of COMPTEL tools.
virtual const GCOMRoi & roi(void) const
Return Region of Interest.
void lb(const double &l, const double &b)
Set galactic sky direction (radians)
Definition: GSkyDir.cpp:251
void init_members(void)
Initialise class members.
Definition: GEvents.cpp:176
Abstract event atom container class.
Definition: GEventList.hpp:49
COMPTEL event list class definition.
void init_members(void)
Initialise class members.
Definition: GEventList.cpp:143
const GTime & tstop(void) const
Return stop time.
Definition: GEvents.hpp:158
GEbounds m_ebounds
Energy boundaries covered by events.
Definition: GEvents.hpp:111
void read_events(const GFitsTable &table)
Read COMPTEL events from FITS table.
Time class.
Definition: GTime.hpp:55
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
void append(const GTime &tstart, const GTime &tstop)
Append Good Time Interval.
Definition: GGti.cpp:269
Implementation of support function used by COMPTEL classes.
std::string centre(const std::string &s, const int &n, const char &c= ' ')
Centre string to achieve a length of n characters.
Definition: GTools.cpp:1118
virtual GCOMEventList & operator=(const GCOMEventList &list)
Assignment operator.
GCOMRoi m_roi
Region of interest.
virtual GCOMEventList * clone(void) const
Clone event list.
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
virtual void load(const GFilename &filename)
Load COMPTEL events from FITS file.
virtual int size(void) const
Return number of events in list.
Energy boundaries container class.
Definition: GEbounds.hpp:60
void dir(const GSkyDir &dir)
Set event scatter direction.
Filename class.
Definition: GFilename.hpp:62
void free_members(void)
Delete class members.
Definition: GEvents.cpp:206
Abstract interface for FITS table column.
virtual void save(const GFilename &filename, const bool &clobber=false) const
Save COMPTEL events.
int com_tjd(const GTime &time)
Convert GTime in COMPTEL TJD.
Definition: GCOMTools.cpp:91
void remove(const int &index, const int &number=1)
Remove events from event list.
virtual GEventList & operator=(const GEventList &list)
Assignment operator.
Definition: GEventList.cpp:104
virtual void clear(void)
Clear region of interest.
Definition: GCOMRoi.cpp:166
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
GGti m_gti
Good time intervals covered by events.
Definition: GEvents.hpp:112
GChatter
Definition: GTypemaps.hpp:33
int integer(const std::string &keyname) const
Return card value as integer.
Definition: GFitsHDU.hpp:436
void append(const GCOMEventAtom &event)
Append event to event list.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL event list information.
void copy_members(const GCOMEventList &list)
Copy class members.
COMPTEL region of interest class.
Definition: GCOMRoi.hpp:44
const int & nrows(void) const
Return number of rows in table.
Definition: GFitsTable.hpp:119
virtual void write(GFits &file) const
Write COMPTEL event list into FITS file.
void free_members(void)
Delete class members.
Definition: GEventList.cpp:165
virtual GCOMEventAtom * operator[](const int &index)
COMPTEL event atom access operator.
int com_tics(const GTime &time)
Convert GTime in COMPTEL tics.
Definition: GCOMTools.cpp:125
virtual int integer(const int &row, const int &inx=0) const =0
virtual void clear(void)
Clear COMPTEL event list.
virtual double real(const int &row, const int &inx=0) const =0
const GEnergy & emin(void) const
Return minimum energy.
Definition: GEvents.hpp:170
std::vector< GCOMEventAtom > m_events
Events.
double keV(void) const
Return energy in keV.
Definition: GEnergy.cpp:306
GCOMEventList(void)
Void constructor.
void lb_deg(const double &l, const double &b)
Set galactic sky direction (degrees)
Definition: GSkyDir.cpp:278
virtual int number(void) const
Return number of events in list.
Exception handler interface definition.
COMPTEL event list class.
const GEbounds & ebounds(void) const
Return energy boundaries.
Definition: GEvents.hpp:122
virtual void read(const GFits &file)
Read COMPTEL events from FITS file.
Interface for the region of interest classes.
Definition: GRoi.hpp:48
const double rad2deg
Definition: GMath.hpp:44
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
const GTime & tstart(void) const
Return start time.
Definition: GEvents.hpp:146
Sky direction class.
Definition: GSkyDir.hpp:62
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
double tofcor(const double &d1e, const double &d2e, double tof) const
Compute TOF correction.
const GEnergy & emax(void) const
Return maximum energy.
Definition: GEvents.hpp:182
void init_members(void)
Initialise class members.
Interface for the COMPTEL instrument direction class.
Definition: GCOMInstDir.hpp:45
Mathematical function definitions.
virtual ~GCOMEventList(void)
Destructor.
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:489