Add D1 and D2 module statistics

Jürgen Knödlseder
1 parent 22933023
Showing 8 changed files with 284 additions and 12 deletions
ChangeLog
NEWS
README.md
doc/source/admin/release_history/2.0.rst
inst/com/include/GCOMSelection.hpp
inst/com/pyext/GCOMSelection.i
inst/com/src/GCOMDri.cpp
inst/com/src/GCOMSelection.cpp
-2021-12-01
+2021-12-02
         * Version 2.0.0 released
         ========================
-        Add GCOMSelection::fpmtflag() method (#3917)
+	Add GCOMSelection::fpmtflag(), GCOMSelection::use_d1() and 
+          GCOMSelection::use_d2() methods (#3917)
         Add GCOMModelDRBPhibarBins class (#3913)
         Compute spectral model gradients for negative factor values
         Add likelihood fit event selection hook (#3912)
@@ -114,6 +114,8 @@ The following methods have been added:
 - GCOMSelection::phase_curve()
 - GCOMSelection::phases()
 - GCOMSelection::fpmtflag()
+- GCOMSelection::use_d1()
+- GCOMSelection::use_d2()
 - gammalib::tmpnam()
 - gammalib::compare()
 - gammalib::xml_check_parnum()
@@ -389,6 +391,11 @@ spatial model gradients (#3203).
 22. COMPTEL interface
 ---------------------
+Added GCOMSelection::use_d1() and GCOMSelection::use_d2() methods to control the usage
+of D1 and D2 modules. The control logic was implemented in GCOMDri and header keywords
+with per module information are now written into the DRE FITS files. A D1 and D2
+module usage string is also written into DRE and DRG files (#3917).
+
 Added GCOMSelection::fpmtflag() methods to select the method how D2 modules with
 failed PMTs will be handled. If the flag is set to 0, D2 modules with failed PMTs
 will be discarded, if set to 1 modules with failed PMT will be included, and if set
 GammaLib information
 ====================
-* Version: 2.0.0.dev (1 December 2021)
+* Version: 2.0.0.dev (2 December 2021)
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
 .. _2.0:
-GammaLib 2.0 release (30 November 2021)
-=======================================
+GammaLib 2.0 release (2 December 2021)
+======================================
 Overview
 --------
@@ -37,6 +37,8 @@ Bug fixes
 Improvements
 ------------
+* [`3917 <https://cta-redmine.irap.omp.eu/issues/3913>`_] -
+  Add ``GCOMSelection::fpmtflag()``, ``GCOMSelection::use_d1()`` and ``GCOMSelection::use_d2()`` methods
 * [`3913 <https://cta-redmine.irap.omp.eu/issues/3913>`_] -
   Add ``GCOMModelDRBPhibarBins`` class
 * [`3912 <https://cta-redmine.irap.omp.eu/issues/3912>`_] -
@@ -97,6 +97,10 @@ public:
     void                            vetoflag_max(const int& vetoflag_max);
     const int&                      fpmtflag(void) const;
     void                            fpmtflag(const int& fpmtflag);
+    const bool&                     use_d1(const int& id1) const;
+    void                            use_d1(const int& id1, const bool& use);
+    const bool&                     use_d2(const int& id2) const;
+    void                            use_d2(const int& id2, const bool& use);
     const GModelTemporalPhaseCurve& phase_curve(void) const;
     void                            phase_curve(const GModelTemporalPhaseCurve& phase_curve);
     const GPhases&                  phases(void) const;
@@ -128,6 +132,8 @@ protected:
     int                      m_fpmtflag;     //!< D2 PMT failures
                                              //!< (0: exclude modules, 1: include modules,
                                              //!<  2: use FPM information)
+    bool                     m_use_d1[7];    //!< D1 module usage
+    bool                     m_use_d2[14];   //!< D2 module usage
     GModelTemporalPhaseCurve m_phase_curve;  //!< Phase curve for phase selection
     GPhases                  m_phases;       //!< Phases for phase selection
@@ -154,6 +160,8 @@ protected:
     mutable int m_num_d1module_off;    //!< Number of events excluded since D1 module off
     mutable int m_num_d2module_off;    //!< Number of events excluded since D2 module off
     mutable int m_num_fpmt;            //!< Number of events excluded due to failed PMT
+    mutable int m_num_d1[7];           //!< Number of events per D1 module
+    mutable int m_num_d2[14];          //!< Number of events per D2 module
 };
@@ -80,6 +80,10 @@ public:
     void                            vetoflag_max(const int& vetoflag_max);
     const int&                      fpmtflag(void) const;
     void                            fpmtflag(const int& fpmtflag);
+    const bool&                     use_d1(const int& id1) const;
+    void                            use_d1(const int& id1, const bool& use);
+    const bool&                     use_d2(const int& id2) const;
+    void                            use_d2(const int& id2, const bool& use);
     const GModelTemporalPhaseCurve& phase_curve(void) const;
     void                            phase_curve(const GModelTemporalPhaseCurve& phase_curve);
     const GPhases&                  phases(void) const;
@@ -1309,12 +1309,33 @@ void GCOMDri::write_attributes(GFitsHDU* hdu) const
     }
     // If there was an event seletion then write selection set. Otherwise write
-    // handling of D2 PMT failures.
+    // handling D1 and D2 module usage and handling mode of D2 PMT failures.
     if (m_has_selection) {
         m_selection.write(*hdu);
     }
     else {
+
+        // Write D2 PMT failure handling flag
         hdu->card("D2FPMT", m_selection.fpmtflag(), "D2 PMT failure handling");
+
+        // Write D1 module usage
+        std::string d1use = "0000000";
+        for (int i = 0; i < 7; ++i) {
+            if (m_selection.use_d1(i)) {
+                d1use[i] = '1';
+            }
+        }
+        hdu->card("D1USE", d1use, "D1 module usage");
+
+        // Write D2 module usage
+        std::string d2use = "00000000000000";
+        for (int i = 0; i < 14; ++i) {
+            if (m_selection.use_d2(i)) {
+                d2use[i] = '1';
+            }
+        }
+        hdu->card("D2USE", d2use, "D2 module usage");
+
     }
     // If there was a valid Earth horizon cut then write it
@@ -1407,6 +1428,11 @@ double GCOMDri::compute_geometry(const int&amp;           tjd,
         // Loop over all D2 modules
         for (int id2 = 0; id2 < 14; ++id2) {
+            // Skip if coresponding modules should not be used
+            if (!select.use_d1(id1) || !select.use_d2(id2)) {
+                continue;
+            }
+
             // Get D2 module status
             int d2status = status.d2status(tjd, id2+1);
@@ -37,6 +37,10 @@
 /* __ Method name definitions ____________________________________________ */
 #define G_FPMTFLAG                            "GCOMSelection::fpmtflag(int&)"
+#define G_USE_D1_GET                            "GCOMSelection::use_d1(int&)"
+#define G_USE_D1_SET                     "GCOMSelection::use_d1(int&, bool&)"
+#define G_USE_D2_GET                            "GCOMSelection::use_d2(int&)"
+#define G_USE_D2_SET                     "GCOMSelection::use_d2(int&, bool&)"
 /* __ Macros _____________________________________________________________ */
@@ -190,6 +194,12 @@ void GCOMSelection::init_statistics(void) const
     m_num_d1module_off    = 0;
     m_num_d2module_off    = 0;
     m_num_fpmt            = 0;
+    for (int i = 0; i < 7; ++i) {
+        m_num_d1[i] = 0;
+    }
+    for (int i = 0; i < 14; ++i) {
+        m_num_d2[i] = 0;
+    }
     // Return
     return;
@@ -304,9 +314,19 @@ bool GCOMSelection::use_event(const GCOMEventAtom&amp; event) const
         int d1status = status.d1status(tjd, id1);
         int d2status = status.d2status(tjd, id2);
+        // Exclude event if the coresponding modules should not be used
+        if (!m_use_d1[id1-1]) {
+            m_num_d1module_off++;
+            use = false;
+        }
+        else if (!m_use_d2[id2-1]) {
+            m_num_d2module_off++;
+            use = false;
+        }
+
         // Exclude event if the corresponding modules are signalled
         // as off
-        if (d1status != 1) {
+        else if (d1status != 1) {
             m_num_d1module_off++;
             use = false;
         }
@@ -319,16 +339,18 @@ bool GCOMSelection::use_event(const GCOMEventAtom&amp; event) const
         // failure flag
         else if (d2status > 1) {
-            // If failure flag = 0 then reject event
+            // If failure flag = 0 then reject event and signal that the
+            // module was off
             if (m_fpmtflag == 0) {
-                m_num_fpmt++;
+                m_num_d2module_off++;
                 use = false;
             }
             // ... otherwise if failure flag = 2 and the module has a
             // valid exclusion circle then reject event if it is
             // within the exclusion circle; if the module has no
-            // exclusion region then reject the entire module.
+            // exclusion region then reject the entire module and
+            // signal that it was off.
             else if (m_fpmtflag == 2) {
                 double r = gammalib::com_exd2r(id2);
                 if (r >= 0.1) {
@@ -342,7 +364,7 @@ bool GCOMSelection::use_event(const GCOMEventAtom&amp; event) const
                     }
                 }
                 else {
-                    m_num_fpmt++;
+                    m_num_d2module_off++;
                     use = false;
                 }
             }
@@ -353,7 +375,16 @@ bool GCOMSelection::use_event(const GCOMEventAtom&amp; event) const
     // Update acceptance and rejection statistics
     if (use) {
+
+        // Extract module IDs from MODCOM
+        int id2 = (event.modcom()-1)/7 + 1;      // [1-14]
+        int id1 =  event.modcom() - (id2-1) * 7; // [1-7]
+
+        // Update statistics
         m_num_events_used++;
+        m_num_d1[id1-1]++;
+        m_num_d2[id2-1]++;
+
     }
     else {
         m_num_events_rejected++;
@@ -395,6 +426,96 @@ void GCOMSelection::fpmtflag(const int&amp; fpmtflag)
 /***********************************************************************//**
+ * @brief Return D1 module usage flag
+ *
+ * @param[in] id1 D1 module identifier [0,...,6].
+ * @return True if D1 module @p id1 should be used.
+ ***************************************************************************/
+const bool& GCOMSelection::use_d1(const int& id1) const
+{
+    // Check module validity range
+    if (id1 < 0 or id1 > 6) {
+        std::string msg = "Invalid D1 module identifier "+gammalib::str(id1)+
+                          " specified. The D1 module identifier needs to be "
+                          "comprised within 0 and 6.";
+        throw GException::invalid_argument(G_USE_D1_GET, msg);
+    }
+
+    // Return usage flag
+    return (m_use_d1[id1]);
+}
+
+
+/***********************************************************************//**
+ * @brief Set D1 module usage flag
+ *
+ * @param[in] id1 D1 module identifier [0,...,6].
+ * @param[in] use D1 module usage.
+ ***************************************************************************/
+void GCOMSelection::use_d1(const int& id1, const bool& use)
+{
+    // Check module validity range
+    if (id1 < 0 or id1 > 6) {
+        std::string msg = "Invalid D1 module identifier "+gammalib::str(id1)+
+                          " specified. The D1 module identifier needs to be "
+                          "comprised within 0 and 6.";
+        throw GException::invalid_argument(G_USE_D1_SET, msg);
+    }
+
+    // Set usage flag
+    m_use_d1[id1] = use;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return D2 module usage flag
+ *
+ * @param[in] id2 D2 module identifier [0,...,13].
+ * @return True if D2 module @p id2 should be used.
+ ***************************************************************************/
+const bool& GCOMSelection::use_d2(const int& id2) const
+{
+    // Check module validity range
+    if (id2 < 0 or id2 > 13) {
+        std::string msg = "Invalid D2 module identifier "+gammalib::str(id2)+
+                          " specified. The D2 module identifier needs to be "
+                          "comprised within 0 and 13.";
+        throw GException::invalid_argument(G_USE_D2_GET, msg);
+    }
+
+    // Return usage flag
+    return (m_use_d2[id2]);
+}
+
+
+/***********************************************************************//**
+ * @brief Set D2 module usage flag
+ *
+ * @param[in] id2 D2 module identifier [0,...,13].
+ * @param[in] use D2 module usage.
+ ***************************************************************************/
+void GCOMSelection::use_d2(const int& id2, const bool& use)
+{
+    // Check module validity range
+    if (id2 < 0 or id2 > 13) {
+        std::string msg = "Invalid D2 module identifier "+gammalib::str(id2)+
+                          " specified. The D2 module identifier needs to be "
+                          "comprised within 0 and 13.";
+        throw GException::invalid_argument(G_USE_D2_SET, msg);
+    }
+
+    // Set usage flag
+    m_use_d2[id2] = use;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Read selection set from FITS HDU keywords
  *
  * @param[in] hdu FITS HDU.
@@ -425,7 +546,6 @@ void GCOMSelection::read(const GFitsHDU&amp; hdu)
         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");
@@ -449,6 +569,26 @@ void GCOMSelection::read(const GFitsHDU&amp; hdu)
         m_fpmtflag = hdu.integer("D2FPMT");
     }
+    // Read D1 module usage flag
+    if (hdu.has_card("D1USE")) {
+        std::string use = gammalib::strip_whitespace(hdu.string("D1USE"));
+        if (use.length() == 7) {
+            for (int i = 0; i < 7; ++i) {
+                m_use_d1[i] = (use[i] == '1');
+            }
+        }
+    }
+
+    // Read D2 module usage flag
+    if (hdu.has_card("D2USE")) {
+        std::string use = gammalib::strip_whitespace(hdu.string("D2USE"));
+        if (use.length() == 14) {
+            for (int i = 0; i < 14; ++i) {
+                m_use_d2[i] = (use[i] == '1');
+            }
+        }
+    }
+
     // Read selection statistics
     if (hdu.has_card("NEVCHK")) {
         m_num_events_checked = hdu.integer("NEVCHK");
@@ -490,6 +630,18 @@ void GCOMSelection::read(const GFitsHDU&amp; hdu)
     if (hdu.has_card("NNOSCT")) {
         m_num_no_scatter = hdu.integer("NNOSCT");
     }
+    for (int i = 0; i < 7; ++i) {
+        std::string key = "ND1-"+gammalib::str(i+1,"%2.2d");
+        if (hdu.has_card(key)) {
+             m_num_d1[i] = hdu.integer(key);
+        }
+    }
+    for (int i = 0; i < 14; ++i) {
+        std::string key = "ND2-"+gammalib::str(i+1,"%2.2d");
+        if (hdu.has_card(key)) {
+             m_num_d2[i] = hdu.integer(key);
+        }
+    }
     if (hdu.has_card("ND1OFF")) {
         m_num_d1module_off = hdu.integer("ND1OFF");
     }
@@ -543,6 +695,24 @@ void GCOMSelection::write(GFitsHDU&amp; hdu) const
     // Write D2 PMT failure handling flag
     hdu.card("D2FPMT",  m_fpmtflag, "D2 PMT failure handling");
+    // Write D1 module usage
+    std::string d1use = "0000000";
+    for (int i = 0; i < 7; ++i) {
+        if (m_use_d1[i]) {
+            d1use[i] = '1';
+        }
+    }
+    hdu.card("D1USE", d1use, "D1 module usage");
+
+    // Write D2 module usage
+    std::string d2use = "00000000000000";
+    for (int i = 0; i < 14; ++i) {
+        if (m_use_d2[i]) {
+            d2use[i] = '1';
+        }
+    }
+    hdu.card("D2USE", d2use, "D2 module usage");
+
     // Write selection statistics
     hdu.card("NEVCHK",  m_num_events_checked,  "Number of checked events");
     hdu.card("NEVUSE",  m_num_events_used,     "Number of used events");
@@ -563,6 +733,16 @@ void GCOMSelection::write(GFitsHDU&amp; hdu) const
     hdu.card("NVFLHI",  m_num_vetoflag_max,    "Number of events > VFLMAX");
     hdu.card("NNOSCT",  m_num_no_scatter,      "Number of events w/o scatter angle");
     hdu.card("NINVMT",  m_num_invalid_modcom,  "Number of events with invalid minitelescope");
+    for (int i = 0; i < 7; ++i) {
+        std::string key = "ND1-"+gammalib::str(i+1,"%2.2d");
+        std::string com = "Number of events in D1 module "+gammalib::str(i+1);
+        hdu.card(key, m_num_d1[i], com);
+    }
+    for (int i = 0; i < 14; ++i) {
+        std::string key = "ND2-"+gammalib::str(i+1,"%2.2d");
+        std::string com = "Number of events in D2 module "+gammalib::str(i+1);
+        hdu.card(key, m_num_d2[i], com);
+    }
     hdu.card("ND1OFF",  m_num_d1module_off,    "Number of events with D1 module off");
     hdu.card("ND2OFF",  m_num_d2module_off,    "Number of events with D2 module off");
     hdu.card("ND2FPM",  m_num_fpmt,            "Number of events excluded due to failed PMTs");
@@ -646,6 +826,30 @@ std::string GCOMSelection::print(const GChatter&amp; chatter) const
         result.append(gammalib::str(m_vetoflag_max)+"] < ");
         result.append(gammalib::str(m_num_vetoflag_max));
+        // Append D1 module usage
+        for (int i = 0; i < 7; ++i) {
+            result.append("\n"+gammalib::parformat("D1 module "+gammalib::str(i+1)));
+            result.append(gammalib::str(m_num_d1[i]));
+            if (m_use_d1[i]) {
+                result.append(" [use]");
+            }
+            else {
+                result.append(" [don't use]");
+            }
+        }
+
+        // Append D2 module usage
+        for (int i = 0; i < 14; ++i) {
+            result.append("\n"+gammalib::parformat("D2 module "+gammalib::str(i+1)));
+            result.append(gammalib::str(m_num_d2[i]));
+            if (m_use_d2[i]) {
+                result.append(" [use]");
+            }
+            else {
+                result.append(" [don't use]");
+            }
+        }
+
         // Append D1 module off
         result.append("\n"+gammalib::parformat("D1 modules off"));
         result.append(gammalib::str(m_num_d1module_off));
@@ -707,6 +911,12 @@ void GCOMSelection::init_members(void)
     m_fpmtflag     =     0;  //!< D2 PMT failure flag
     m_phase_curve.clear();
     m_phases.clear();
+    for (int i = 0; i < 7; ++i) {
+        m_use_d1[i] = true;
+    }
+    for (int i = 0; i < 14; ++i) {
+        m_use_d2[i] = true;
+    }
     // Initialise statistics
     init_statistics();
@@ -741,6 +951,12 @@ void GCOMSelection::copy_members(const GCOMSelection&amp; select)
     m_fpmtflag     = select.m_fpmtflag;
     m_phase_curve  = select.m_phase_curve;
     m_phases       = select.m_phases;
+    for (int i = 0; i < 7; ++i) {
+        m_use_d1[i] = select.m_use_d1[i];
+    }
+    for (int i = 0; i < 14; ++i) {
+        m_use_d2[i] = select.m_use_d2[i];
+    }
     // Copy statistics
     m_num_events_checked  = select.m_num_events_checked;
@@ -766,6 +982,14 @@ void GCOMSelection::copy_members(const GCOMSelection&amp; select)
     m_num_d2module_off    = select.m_num_d2module_off;
     m_num_fpmt            = select.m_num_fpmt;
+    // Copy module usage
+    for (int i = 0; i < 7; ++i) {
+        m_num_d1[i] = select.m_num_d1[i];
+    }
+    for (int i = 0; i < 14; ++i) {
+        m_num_d2[i] = select.m_num_d2[i];
+    }
+
     // Return
     return;
 }