Sergio Hernández Cadena / gammalib

Commit bf1c6bc023c8a1e5a7c2a60f916e9638ef93d2ba

Authored by Jürgen Knödlseder
1 parent 546fff3d

Implement reading and writing of selection keywords and add ToF correction factor

Inline Side-by-side
Showing 6 changed files with 184 additions and 13 deletions
inst/com/include/GCOMDri.hpp
  View file @bf1c6bc
1 1 /***************************************************************************
2 2 * GCOMDri.hpp - COMPTEL Data Space class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2017-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2017-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -38,6 +38,7 @@
38 38 /* __ Forward declarations _______________________________________________ */
39 39 class GFilename;
40 40 class GFits;
  41 +class GFitsHDU;
41 42 class GFitsImage;
42 43 class GModel;
43 44 class GModelSky;
... ... @@ -94,6 +95,7 @@ public:
94 95 void gti(const GGti& gti);
95 96 const double& phimin(void) const;
96 97 const double& phibin(void) const;
  98 + const double& tof_correction(void) const;
97 99 void compute_dre(const GCOMObservation& obs,
98 100 const GCOMSelection& select = GCOMSelection(),
99 101 const double& zetamin = 5.0);
... ... @@ -112,6 +114,8 @@ public:
112 114 void read(const GFitsImage& image);
113 115 void write(GFits& fits,
114 116 const std::string& extname = "") const;
  117 + void read_attributes(const GFitsHDU& hdu);
  118 + void write_attributes(GFitsHDU& hdu) const;
115 119  
116 120 protected:
117 121 // Protected methods
... ... @@ -136,6 +140,7 @@ protected:
136 140 GGti m_gti; //!< Good Time Intervals of data cube
137 141 double m_phimin; //!< Phibar minimum (deg)
138 142 double m_phibin; //!< Phibar binsize (deg)
  143 + double m_tofcor; //!< ToF correction
139 144  
140 145 // Computation statistics
141 146 GTime m_tstart; //!< Selection start time
... ... @@ -145,6 +150,7 @@ protected:
145 150 int m_num_skipped_superpackets; //!< Number of skipped superpackets
146 151  
147 152 // Selection parameters
  153 + bool m_has_selection; //!< Signal that selection was applied
148 154 GCOMSelection m_selection; //!< Selection parameters
149 155 double m_zetamin; //!< Minimum zeta angle
150 156 };
... ... @@ -352,4 +358,19 @@ const double&amp; GCOMDri::phibin(void) const
352 358 return (m_phibin);
353 359 }
354 360  
  361 +
  362 +/***********************************************************************//**
  363 + * @brief Return ToF correction factor
  364 + *
  365 + * @return ToF correction factor.
  366 + *
  367 + * Returns the ToF correction factor that corrects for the event selection
  368 + * in a ToF window.
  369 + ***************************************************************************/
  370 +inline
  371 +const double& GCOMDri::tof_correction(void) const
  372 +{
  373 + return (m_tofcor);
  374 +}
  375 +
355 376 #endif /* GCOMDRI_HPP */
... ...
inst/com/include/GCOMSelection.hpp
  View file @bf1c6bc
1 1 /***************************************************************************
2 2 * GCOMSelection.hpp - COMPTEL selection set class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2017-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,6 +32,7 @@
32 32 #include "GBase.hpp"
33 33  
34 34 /* __ Forward declarations _______________________________________________ */
  35 +class GFitsHDU;
35 36 class GCOMEventAtom;
36 37  
37 38 /* __ Constants __________________________________________________________ */
... ... @@ -64,13 +65,15 @@ public:
64 65 // Other methods
65 66 void init_statistics(void) const;
66 67 bool use_event(const GCOMEventAtom& event) const;
  68 + void read(const GFitsHDU& hdu);
  69 + void write(GFitsHDU& hdu) const;
67 70  
68 71 protected:
69 72 // Protected methods
70 73 void init_members(void);
71 74 void copy_members(const GCOMSelection& select);
72 75 void free_members(void);
73   -
  76 +
74 77 // Protected members
75 78 double m_e1_min; //!< Minimum D1 energy deposit (MeV)
76 79 double m_e1_max; //!< Maximum D1 energy deposit (MeV)
... ... @@ -103,8 +106,8 @@ protected:
103 106 mutable int m_num_zeta_max; //!< Number of events above Zeta threshold
104 107 mutable int m_num_reflag_min; //!< Number of events below rejection flag threshold
105 108 mutable int m_num_reflag_max; //!< Number of events above rejection flag threshold
106   - mutable int m_num_vetoflag_min; //!< Number of events below rejection flag threshold
107   - mutable int m_num_vetoflag_max; //!< Number of events above rejection flag threshold
  109 + mutable int m_num_vetoflag_min; //!< Number of events below veto flag threshold
  110 + mutable int m_num_vetoflag_max; //!< Number of events above veto flag threshold
108 111 mutable int m_num_no_scatter; //!< Number of events without scatter angle
109 112 mutable int m_num_invalid_modcom; //!< Number of events with invalid minitelescopes
110 113 };
... ...
inst/com/pyext/GCOMDri.i
  View file @bf1c6bc
1 1 /***************************************************************************
2 2 * GCOMDri.i - COMPTEL Data Space class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2017-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2017-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -65,6 +65,7 @@ public:
65 65 void gti(const GGti& gti);
66 66 const double& phimin(void) const;
67 67 const double& phibin(void) const;
  68 + const double& tof_correction(void) const;
68 69 void compute_dre(const GCOMObservation& obs,
69 70 const GCOMSelection& select = GCOMSelection(),
70 71 const double& zetamin = 5.0);
... ... @@ -83,6 +84,8 @@ public:
83 84 void read(const GFitsImage& image);
84 85 void write(GFits& fits,
85 86 const std::string& extname = "") const;
  87 + void read_attributes(const GFitsHDU& hdu);
  88 + void write_attributes(GFitsHDU& hdu) const;
86 89 };
87 90  
88 91  
... ...
inst/com/pyext/GCOMSelection.i
  View file @bf1c6bc
1 1 /***************************************************************************
2 2 * GCOMSelection.i - COMPTEL selection set class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2017-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -50,6 +50,8 @@ public:
50 50 // Other methods
51 51 void init_statistics(void) const;
52 52 bool use_event(const GCOMEventAtom& event) const;
  53 + void read(const GFitsHDU& hdu);
  54 + void write(GFitsHDU& hdu) const;
53 55 };
54 56  
55 57  
... ...
inst/com/src/GCOMDri.cpp
  View file @bf1c6bc
... ... @@ -293,6 +293,11 @@ void GCOMDri::compute_dre(const GCOMObservation&amp; obs,
293 293 // Initialise selection statistics
294 294 select.init_statistics();
295 295  
  296 + // Signal that event selection was used (this will enable writing into
  297 + // the FITS HDU) and store Earth horizon cut
  298 + m_has_selection = true;
  299 + m_zetamin = zetamin;
  300 +
296 301 // Initialise statistics
297 302 int num_used_events = 0;
298 303 int num_event_outside_sp = 0;
... ... @@ -379,7 +384,7 @@ void GCOMDri::compute_dre(const GCOMObservation&amp; obs,
379 384 continue;
380 385 }
381 386  
382   - // Apply event selectio
  387 + // Apply event selection
383 388 if (!select.use_event(*event)) {
384 389 continue;
385 390 }
... ... @@ -978,6 +983,10 @@ std::string GCOMDri::print(const GChatter&amp; chatter) const
978 983 result.append(gammalib::str(phimax)+" deg");
979 984 result.append("\n"+gammalib::parformat("Phibar bin size"));
980 985 result.append(gammalib::str(m_phibin)+" deg");
  986 + if (m_tofcor > 1.0) {
  987 + result.append("\n"+gammalib::parformat("ToF correction"));
  988 + result.append(gammalib::str(m_tofcor));
  989 + }
981 990  
982 991 // Append energy boundaries
983 992 result.append("\n"+m_ebounds.print(gammalib::reduce(chatter)));
... ... @@ -998,6 +1007,11 @@ std::string GCOMDri::print(const GChatter&amp; chatter) const
998 1007 result.append("\n"+gammalib::parformat("Skipped superpackets"));
999 1008 result.append(gammalib::str(m_num_skipped_superpackets));
1000 1009  
  1010 + // Append selection
  1011 + if (m_has_selection) {
  1012 + result.append("\n"+m_selection.print(gammalib::reduce(chatter)));
  1013 + }
  1014 +
1001 1015 } // endif: chatter was not silent
1002 1016  
1003 1017 // Return result
... ... @@ -1023,11 +1037,13 @@ void GCOMDri::init_members(void)
1023 1037 m_gti.clear();
1024 1038 m_phimin = 0.0;
1025 1039 m_phibin = 0.0;
  1040 + m_tofcor = 1.0;
1026 1041  
1027 1042 // Initialise statistics
1028 1043 init_statistics();
1029 1044  
1030 1045 // Initialise selection parameters
  1046 + m_has_selection = false;
1031 1047 m_selection.clear();
1032 1048 m_zetamin = -90.0; // Signals no selection
1033 1049  
... ... @@ -1050,6 +1066,7 @@ void GCOMDri::copy_members(const GCOMDri&amp; dri)
1050 1066 m_gti = dri.m_gti;
1051 1067 m_phimin = dri.m_phimin;
1052 1068 m_phibin = dri.m_phibin;
  1069 + m_tofcor = dri.m_tofcor;
1053 1070  
1054 1071 // Copy statistics
1055 1072 m_tstart = dri.m_tstart;
... ... @@ -1059,8 +1076,9 @@ void GCOMDri::copy_members(const GCOMDri&amp; dri)
1059 1076 m_num_skipped_superpackets = dri.m_num_skipped_superpackets;
1060 1077  
1061 1078 // Copy selection parameters
1062   - m_selection = dri.m_selection;
1063   - m_zetamin = dri.m_zetamin;
  1079 + m_has_selection = dri.m_has_selection;
  1080 + m_selection = dri.m_selection;
  1081 + m_zetamin = dri.m_zetamin;
1064 1082  
1065 1083 // Return
1066 1084 return;
... ... @@ -1219,6 +1237,19 @@ void GCOMDri::read_attributes(const GFitsHDU* hdu)
1219 1237 m_ebounds.clear();
1220 1238 }
1221 1239  
  1240 + // Read selection set
  1241 + m_selection.read(*hdu);
  1242 +
  1243 + // Optionally read Earth horizon cut
  1244 + if (hdu->has_card("EHAMIN")) {
  1245 + m_zetamin = hdu->real("EHAMIN");
  1246 + }
  1247 +
  1248 + // Optionally read ToF correction
  1249 + if (hdu->has_card("TOFCOR")) {
  1250 + m_tofcor = hdu->real("TOFCOR");
  1251 + }
  1252 +
1222 1253 // Optionally read superpacket statistics
1223 1254 if (hdu->has_card("NSPINP")) {
1224 1255 m_num_superpackets = hdu->integer("NSPINP");
... ... @@ -1282,6 +1313,21 @@ void GCOMDri::write_attributes(GFitsHDU* hdu) const
1282 1313 hdu->card("E_MAX", m_ebounds.emax().MeV(), "[MeV] Upper bound of energy range");
1283 1314 }
1284 1315  
  1316 + // If there was an event seletion then write selection set
  1317 + if (m_has_selection) {
  1318 + m_selection.write(*hdu);
  1319 + }
  1320 +
  1321 + // If there was a valid Earth horizon cut then write it
  1322 + if (m_zetamin != -90.0) {
  1323 + hdu->card("EHAMIN", m_zetamin, "[deg] Minimum Earth horizon - Phibar cut");
  1324 + }
  1325 +
  1326 + // If there is a ToF correction then write it
  1327 + if (m_tofcor > 1.0) {
  1328 + hdu->card("TOFCOR", m_tofcor, "ToF correction");
  1329 + }
  1330 +
1285 1331 // Write superpacket statistics
1286 1332 hdu->card("NSPINP", m_num_superpackets, "Number of input superpackets");
1287 1333 hdu->card("NSPUSE", m_num_used_superpackets, "Number of used superpackets");
... ...
inst/com/src/GCOMSelection.cpp
  View file @bf1c6bc
1 1 /***************************************************************************
2 2 * GCOMSelection.cpp - COMPTEL selection set class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2017-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2017-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -28,6 +28,7 @@
28 28 #ifdef HAVE_CONFIG_H
29 29 #include <config.h>
30 30 #endif
  31 +#include "GFitsHDU.hpp"
31 32 #include "GCOMSelection.hpp"
32 33 #include "GCOMEventAtom.hpp"
33 34  
... ... @@ -290,6 +291,101 @@ bool GCOMSelection::use_event(const GCOMEventAtom&amp; event) const
290 291  
291 292  
292 293 /***********************************************************************//**
  294 + * @brief Read selection set from FITS HDU keywords
  295 + *
  296 + * @param[in] hdu FITS HDU.
  297 + ***************************************************************************/
  298 +void GCOMSelection::read(const GFitsHDU& hdu)
  299 +{
  300 + // Read D1 energy selection
  301 + if (hdu.has_card("D1EMIN") && hdu.has_card("D1EMAX")) {
  302 + m_e1_min = hdu.real("D1EMIN");
  303 + m_e1_max = hdu.real("D1EMAX");
  304 + }
  305 +
  306 + // Read D2 energy selection
  307 + if (hdu.has_card("D2EMIN") && hdu.has_card("D2EMAX")) {
  308 + m_e2_min = hdu.real("D2EMIN");
  309 + m_e2_max = hdu.real("D2EMAX");
  310 + }
  311 +
  312 + // Read ToF selection
  313 + if (hdu.has_card("TOFMIN") && hdu.has_card("TOFMAX")) {
  314 + m_tof_min = hdu.real("TOFMIN");
  315 + m_tof_max = hdu.real("TOFMAX");
  316 + }
  317 +
  318 + // Read PSD selection
  319 + if (hdu.has_card("PSDMIN") && hdu.has_card("PSDMAX")) {
  320 + m_psd_min = hdu.real("PSDMIN");
  321 + m_psd_max = hdu.real("PSDMAX");
  322 + }
  323 +
  324 +
  325 + // Read Earth horizon angle - Phibar selection
  326 + if (hdu.has_card("ZETMIN") && hdu.has_card("ZETMAX")) {
  327 + m_zeta_min = hdu.real("ZETMIN");
  328 + m_zeta_max = hdu.real("ZETMAX");
  329 + }
  330 +
  331 + // Read rejection flag selection
  332 + if (hdu.has_card("RFLMIN") && hdu.has_card("RFLMAX")) {
  333 + m_reflag_min = hdu.integer("RFLMIN");
  334 + m_reflag_max = hdu.integer("RFLMAX");
  335 + }
  336 +
  337 + // Read veto flag selection
  338 + if (hdu.has_card("VFLMIN") && hdu.has_card("VFLMAX")) {
  339 + m_vetoflag_min = hdu.integer("VFLMIN");
  340 + m_vetoflag_max = hdu.integer("VFLMAX");
  341 + }
  342 +
  343 + // Return
  344 + return;
  345 +}
  346 +
  347 +
  348 +/***********************************************************************//**
  349 + * @brief Write selection set keywords into FITS HDU
  350 + *
  351 + * @param[in] hdu FITS HDU.
  352 + ***************************************************************************/
  353 +void GCOMSelection::write(GFitsHDU& hdu) const
  354 +{
  355 + // Write D1 energy selection
  356 + hdu.card("D1EMIN", m_e1_min, "[MeV] Minimum D1 energy deposit");
  357 + hdu.card("D1EMAX", m_e1_max, "[MeV] Maximum D1 energy deposit");
  358 +
  359 + // Write D2 energy selection
  360 + hdu.card("D2EMIN", m_e2_min, "[MeV] Minimum D2 energy deposit");
  361 + hdu.card("D2EMAX", m_e2_max, "[MeV] Maximum D2 energy deposit");
  362 +
  363 + // Write ToF selection
  364 + hdu.card("TOFMIN", m_tof_min, "ToF selection minimum");
  365 + hdu.card("TOFMAX", m_tof_max, "ToF selection maximum");
  366 +
  367 + // Write PSD selection
  368 + hdu.card("PSDMIN", m_psd_min, "PSD selection minimum");
  369 + hdu.card("PSDMAX", m_psd_max, "PSD selection maximum");
  370 +
  371 + // Write Earth horizon angle - Phibar selection
  372 + hdu.card("ZETMIN", m_zeta_min, "[deg] Earth horizon angle - Phibar minimum");
  373 + hdu.card("ZETMAX", m_zeta_max, "[deg] Earth horizon angle - Phibar maximum");
  374 +
  375 + // Write rejection flag selection
  376 + hdu.card("RFLMIN", m_reflag_min, "Rejection flag minimum");
  377 + hdu.card("RFLMAX", m_reflag_max, "Rejection flag maximum");
  378 +
  379 + // Write veto flag selection
  380 + hdu.card("VFLMIN", m_vetoflag_min, "Veto flag minimum");
  381 + hdu.card("VFLMAX", m_vetoflag_max, "Veto flag maximum");
  382 +
  383 + // Return
  384 + return;
  385 +}
  386 +
  387 +
  388 +/***********************************************************************//**
293 389 * @brief Print COMPTEL selection set
294 390 *
295 391 * @param[in] chatter Chattiness.
... ... @@ -402,10 +498,10 @@ void GCOMSelection::init_members(void)
402 498 m_reflag_max = 1000; //!< Maximum rejection flag
403 499 m_vetoflag_min = 0; //!< Minimum veto flag
404 500 m_vetoflag_max = 0; //!< Maximum veto flag
405   -
  501 +
406 502 // Initialise statistics
407 503 init_statistics();
408   -
  504 +
409 505 // Return
410 506 return;
411 507 }
... ...