diff --git a/inst/com/include/GCOMDri.hpp b/inst/com/include/GCOMDri.hpp index 687f97c..a1c705f 100644 --- a/inst/com/include/GCOMDri.hpp +++ b/inst/com/include/GCOMDri.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMDri.hpp - COMPTEL Data Space class * * ----------------------------------------------------------------------- * - * copyright (C) 2017-2019 by Juergen Knoedlseder * + * copyright (C) 2017-2021 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -38,6 +38,7 @@ /* __ Forward declarations _______________________________________________ */ class GFilename; class GFits; +class GFitsHDU; class GFitsImage; class GModel; class GModelSky; @@ -94,6 +95,7 @@ public: void gti(const GGti& gti); const double& phimin(void) const; const double& phibin(void) const; + const double& tof_correction(void) const; void compute_dre(const GCOMObservation& obs, const GCOMSelection& select = GCOMSelection(), const double& zetamin = 5.0); @@ -112,6 +114,8 @@ public: void read(const GFitsImage& image); void write(GFits& fits, const std::string& extname = "") const; + void read_attributes(const GFitsHDU& hdu); + void write_attributes(GFitsHDU& hdu) const; protected: // Protected methods @@ -136,6 +140,7 @@ protected: GGti m_gti; //!< Good Time Intervals of data cube double m_phimin; //!< Phibar minimum (deg) double m_phibin; //!< Phibar binsize (deg) + double m_tofcor; //!< ToF correction // Computation statistics GTime m_tstart; //!< Selection start time @@ -145,6 +150,7 @@ protected: int m_num_skipped_superpackets; //!< Number of skipped superpackets // Selection parameters + bool m_has_selection; //!< Signal that selection was applied GCOMSelection m_selection; //!< Selection parameters double m_zetamin; //!< Minimum zeta angle }; @@ -352,4 +358,19 @@ const double& GCOMDri::phibin(void) const return (m_phibin); } + +/***********************************************************************//** + * @brief Return ToF correction factor + * + * @return ToF correction factor. + * + * Returns the ToF correction factor that corrects for the event selection + * in a ToF window. + ***************************************************************************/ +inline +const double& GCOMDri::tof_correction(void) const +{ + return (m_tofcor); +} + #endif /* GCOMDRI_HPP */ diff --git a/inst/com/include/GCOMSelection.hpp b/inst/com/include/GCOMSelection.hpp index 621ff1d..123fdf7 100644 --- a/inst/com/include/GCOMSelection.hpp +++ b/inst/com/include/GCOMSelection.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMSelection.hpp - COMPTEL selection set class * * ----------------------------------------------------------------------- * - * copyright (C) 2017 by Juergen Knoedlseder * + * copyright (C) 2017-2021 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -32,6 +32,7 @@ #include "GBase.hpp" /* __ Forward declarations _______________________________________________ */ +class GFitsHDU; class GCOMEventAtom; /* __ Constants __________________________________________________________ */ @@ -64,13 +65,15 @@ public: // Other methods void init_statistics(void) const; bool use_event(const GCOMEventAtom& event) const; + void read(const GFitsHDU& hdu); + void write(GFitsHDU& hdu) const; protected: // Protected methods void init_members(void); void copy_members(const GCOMSelection& select); void free_members(void); - + // Protected members double m_e1_min; //!< Minimum D1 energy deposit (MeV) double m_e1_max; //!< Maximum D1 energy deposit (MeV) @@ -103,8 +106,8 @@ protected: mutable int m_num_zeta_max; //!< Number of events above Zeta threshold mutable int m_num_reflag_min; //!< Number of events below rejection flag threshold mutable int m_num_reflag_max; //!< Number of events above rejection flag threshold - mutable int m_num_vetoflag_min; //!< Number of events below rejection flag threshold - mutable int m_num_vetoflag_max; //!< Number of events above rejection flag threshold + mutable int m_num_vetoflag_min; //!< Number of events below veto flag threshold + mutable int m_num_vetoflag_max; //!< Number of events above veto flag threshold mutable int m_num_no_scatter; //!< Number of events without scatter angle mutable int m_num_invalid_modcom; //!< Number of events with invalid minitelescopes }; diff --git a/inst/com/pyext/GCOMDri.i b/inst/com/pyext/GCOMDri.i index 3328fe4..fafc4df 100644 --- a/inst/com/pyext/GCOMDri.i +++ b/inst/com/pyext/GCOMDri.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMDri.i - COMPTEL Data Space class * * ----------------------------------------------------------------------- * - * copyright (C) 2017-2019 by Juergen Knoedlseder * + * copyright (C) 2017-2021 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -65,6 +65,7 @@ public: void gti(const GGti& gti); const double& phimin(void) const; const double& phibin(void) const; + const double& tof_correction(void) const; void compute_dre(const GCOMObservation& obs, const GCOMSelection& select = GCOMSelection(), const double& zetamin = 5.0); @@ -83,6 +84,8 @@ public: void read(const GFitsImage& image); void write(GFits& fits, const std::string& extname = "") const; + void read_attributes(const GFitsHDU& hdu); + void write_attributes(GFitsHDU& hdu) const; }; diff --git a/inst/com/pyext/GCOMSelection.i b/inst/com/pyext/GCOMSelection.i index 6b93708..6b0b531 100644 --- a/inst/com/pyext/GCOMSelection.i +++ b/inst/com/pyext/GCOMSelection.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMSelection.i - COMPTEL selection set class * * ----------------------------------------------------------------------- * - * copyright (C) 2017 by Juergen Knoedlseder * + * copyright (C) 2017-2021 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -50,6 +50,8 @@ public: // Other methods void init_statistics(void) const; bool use_event(const GCOMEventAtom& event) const; + void read(const GFitsHDU& hdu); + void write(GFitsHDU& hdu) const; }; diff --git a/inst/com/src/GCOMDri.cpp b/inst/com/src/GCOMDri.cpp index 913521f..fb8beb2 100644 --- a/inst/com/src/GCOMDri.cpp +++ b/inst/com/src/GCOMDri.cpp @@ -293,6 +293,11 @@ void GCOMDri::compute_dre(const GCOMObservation& obs, // Initialise selection statistics select.init_statistics(); + // Signal that event selection was used (this will enable writing into + // the FITS HDU) and store Earth horizon cut + m_has_selection = true; + m_zetamin = zetamin; + // Initialise statistics int num_used_events = 0; int num_event_outside_sp = 0; @@ -379,7 +384,7 @@ void GCOMDri::compute_dre(const GCOMObservation& obs, continue; } - // Apply event selectio + // Apply event selection if (!select.use_event(*event)) { continue; } @@ -978,6 +983,10 @@ std::string GCOMDri::print(const GChatter& chatter) const result.append(gammalib::str(phimax)+" deg"); result.append("\n"+gammalib::parformat("Phibar bin size")); result.append(gammalib::str(m_phibin)+" deg"); + if (m_tofcor > 1.0) { + result.append("\n"+gammalib::parformat("ToF correction")); + result.append(gammalib::str(m_tofcor)); + } // Append energy boundaries result.append("\n"+m_ebounds.print(gammalib::reduce(chatter))); @@ -998,6 +1007,11 @@ std::string GCOMDri::print(const GChatter& chatter) const result.append("\n"+gammalib::parformat("Skipped superpackets")); result.append(gammalib::str(m_num_skipped_superpackets)); + // Append selection + if (m_has_selection) { + result.append("\n"+m_selection.print(gammalib::reduce(chatter))); + } + } // endif: chatter was not silent // Return result @@ -1023,11 +1037,13 @@ void GCOMDri::init_members(void) m_gti.clear(); m_phimin = 0.0; m_phibin = 0.0; + m_tofcor = 1.0; // Initialise statistics init_statistics(); // Initialise selection parameters + m_has_selection = false; m_selection.clear(); m_zetamin = -90.0; // Signals no selection @@ -1050,6 +1066,7 @@ void GCOMDri::copy_members(const GCOMDri& dri) m_gti = dri.m_gti; m_phimin = dri.m_phimin; m_phibin = dri.m_phibin; + m_tofcor = dri.m_tofcor; // Copy statistics m_tstart = dri.m_tstart; @@ -1059,8 +1076,9 @@ void GCOMDri::copy_members(const GCOMDri& dri) m_num_skipped_superpackets = dri.m_num_skipped_superpackets; // Copy selection parameters - m_selection = dri.m_selection; - m_zetamin = dri.m_zetamin; + m_has_selection = dri.m_has_selection; + m_selection = dri.m_selection; + m_zetamin = dri.m_zetamin; // Return return; @@ -1219,6 +1237,19 @@ void GCOMDri::read_attributes(const GFitsHDU* hdu) m_ebounds.clear(); } + // Read selection set + m_selection.read(*hdu); + + // Optionally read Earth horizon cut + if (hdu->has_card("EHAMIN")) { + m_zetamin = hdu->real("EHAMIN"); + } + + // Optionally read ToF correction + if (hdu->has_card("TOFCOR")) { + m_tofcor = hdu->real("TOFCOR"); + } + // Optionally read superpacket statistics if (hdu->has_card("NSPINP")) { m_num_superpackets = hdu->integer("NSPINP"); @@ -1282,6 +1313,21 @@ void GCOMDri::write_attributes(GFitsHDU* hdu) const hdu->card("E_MAX", m_ebounds.emax().MeV(), "[MeV] Upper bound of energy range"); } + // If there was an event seletion then write selection set + if (m_has_selection) { + m_selection.write(*hdu); + } + + // If there was a valid Earth horizon cut then write it + if (m_zetamin != -90.0) { + hdu->card("EHAMIN", m_zetamin, "[deg] Minimum Earth horizon - Phibar cut"); + } + + // If there is a ToF correction then write it + if (m_tofcor > 1.0) { + hdu->card("TOFCOR", m_tofcor, "ToF correction"); + } + // Write superpacket statistics hdu->card("NSPINP", m_num_superpackets, "Number of input superpackets"); hdu->card("NSPUSE", m_num_used_superpackets, "Number of used superpackets"); diff --git a/inst/com/src/GCOMSelection.cpp b/inst/com/src/GCOMSelection.cpp index e8bdbfb..751ae5c 100644 --- a/inst/com/src/GCOMSelection.cpp +++ b/inst/com/src/GCOMSelection.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCOMSelection.cpp - COMPTEL selection set class * * ----------------------------------------------------------------------- * - * copyright (C) 2017-2020 by Juergen Knoedlseder * + * copyright (C) 2017-2021 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -28,6 +28,7 @@ #ifdef HAVE_CONFIG_H #include #endif +#include "GFitsHDU.hpp" #include "GCOMSelection.hpp" #include "GCOMEventAtom.hpp" @@ -290,6 +291,101 @@ bool GCOMSelection::use_event(const GCOMEventAtom& event) const /***********************************************************************//** + * @brief Read selection set from FITS HDU keywords + * + * @param[in] hdu FITS HDU. + ***************************************************************************/ +void GCOMSelection::read(const GFitsHDU& hdu) +{ + // Read D1 energy selection + if (hdu.has_card("D1EMIN") && hdu.has_card("D1EMAX")) { + m_e1_min = hdu.real("D1EMIN"); + m_e1_max = hdu.real("D1EMAX"); + } + + // Read D2 energy selection + if (hdu.has_card("D2EMIN") && hdu.has_card("D2EMAX")) { + m_e2_min = hdu.real("D2EMIN"); + m_e2_max = hdu.real("D2EMAX"); + } + + // Read ToF selection + if (hdu.has_card("TOFMIN") && hdu.has_card("TOFMAX")) { + m_tof_min = hdu.real("TOFMIN"); + m_tof_max = hdu.real("TOFMAX"); + } + + // Read PSD selection + if (hdu.has_card("PSDMIN") && hdu.has_card("PSDMAX")) { + m_psd_min = hdu.real("PSDMIN"); + m_psd_max = hdu.real("PSDMAX"); + } + + + // Read Earth horizon angle - Phibar selection + if (hdu.has_card("ZETMIN") && hdu.has_card("ZETMAX")) { + m_zeta_min = hdu.real("ZETMIN"); + m_zeta_max = hdu.real("ZETMAX"); + } + + // Read rejection flag selection + if (hdu.has_card("RFLMIN") && hdu.has_card("RFLMAX")) { + m_reflag_min = hdu.integer("RFLMIN"); + m_reflag_max = hdu.integer("RFLMAX"); + } + + // Read veto flag selection + if (hdu.has_card("VFLMIN") && hdu.has_card("VFLMAX")) { + m_vetoflag_min = hdu.integer("VFLMIN"); + m_vetoflag_max = hdu.integer("VFLMAX"); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write selection set keywords into FITS HDU + * + * @param[in] hdu FITS HDU. + ***************************************************************************/ +void GCOMSelection::write(GFitsHDU& hdu) const +{ + // Write D1 energy selection + hdu.card("D1EMIN", m_e1_min, "[MeV] Minimum D1 energy deposit"); + hdu.card("D1EMAX", m_e1_max, "[MeV] Maximum D1 energy deposit"); + + // Write D2 energy selection + hdu.card("D2EMIN", m_e2_min, "[MeV] Minimum D2 energy deposit"); + hdu.card("D2EMAX", m_e2_max, "[MeV] Maximum D2 energy deposit"); + + // Write ToF selection + hdu.card("TOFMIN", m_tof_min, "ToF selection minimum"); + hdu.card("TOFMAX", m_tof_max, "ToF selection maximum"); + + // Write PSD selection + hdu.card("PSDMIN", m_psd_min, "PSD selection minimum"); + hdu.card("PSDMAX", m_psd_max, "PSD selection maximum"); + + // Write Earth horizon angle - Phibar selection + hdu.card("ZETMIN", m_zeta_min, "[deg] Earth horizon angle - Phibar minimum"); + hdu.card("ZETMAX", m_zeta_max, "[deg] Earth horizon angle - Phibar maximum"); + + // Write rejection flag selection + hdu.card("RFLMIN", m_reflag_min, "Rejection flag minimum"); + hdu.card("RFLMAX", m_reflag_max, "Rejection flag maximum"); + + // Write veto flag selection + hdu.card("VFLMIN", m_vetoflag_min, "Veto flag minimum"); + hdu.card("VFLMAX", m_vetoflag_max, "Veto flag maximum"); + + // Return + return; +} + + +/***********************************************************************//** * @brief Print COMPTEL selection set * * @param[in] chatter Chattiness. @@ -402,10 +498,10 @@ void GCOMSelection::init_members(void) m_reflag_max = 1000; //!< Maximum rejection flag m_vetoflag_min = 0; //!< Minimum veto flag m_vetoflag_max = 0; //!< Maximum veto flag - + // Initialise statistics init_statistics(); - + // Return return; }