Implement GCOMDris::compute_drws() method

Jürgen Knödlseder
1 parent d0bfcda6
Showing 5 changed files with 316 additions and 4 deletions
inst/com/include/GCOMDri.hpp
inst/com/include/GCOMDris.hpp
inst/com/pyext/GCOMDris.i
inst/com/src/GCOMDri.cpp
inst/com/src/GCOMDris.cpp
 /***************************************************************************
  *                  GCOMDri.hpp - COMPTEL Data Space class                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017-2022 by Juergen Knoedlseder                         *
+ *  copyright (C) 2017-2023 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -60,6 +60,8 @@ namespace gammalib {
  ***************************************************************************/
 class GCOMDri : public GBase {
  
+    friend class GCOMDris;
+
 public:
     // Constructors and destructors
     GCOMDri(void);
@@ -34,8 +34,6 @@
 #include "GContainer.hpp"
  
 /* __ Forward declarations _______________________________________________ */
-//class GFilename;
-//class GFitsTable;
  
 /* __ Constants __________________________________________________________ */
  
@@ -75,6 +73,9 @@ public:
     void           remove(const int& index);
     void           reserve(const int& num);
     void           extend(const GCOMDris& dris);
+    void           compute_drws(const GCOMObservation& obs,
+                                const GCOMSelection&   select = GCOMSelection(),
+                                const double&          zetamin = 5.0);
     std::string    print(const GChatter& chatter = NORMAL) const;
  
 protected:
@@ -53,6 +53,9 @@ public:
     void        remove(const int& index);
     void        reserve(const int& num);
     void        extend(const GCOMDris& dris);
+    void        compute_drws(const GCOMObservation& obs,
+                             const GCOMSelection&   select = GCOMSelection(),
+                             const double&          zetamin = 5.0);
 };
  
  
@@ -45,7 +45,6 @@
 #include "GCOMOads.hpp"
 #include "GCOMTim.hpp"
 #include "GCOMStatus.hpp"
-#include "GCOMSupport.hpp"
 #include "GCOMObservation.hpp"
 #include "GCOMEventList.hpp"
 #include "GCOMSelection.hpp"
@@ -32,18 +32,28 @@
 #include "GMath.hpp"
 #include "GCOMTools.hpp"
 #include "GCOMSupport.hpp"
+#include "GCOMDri.hpp"
 #include "GCOMDris.hpp"
+#include "GCOMOad.hpp"
+#include "GCOMOads.hpp"
+#include "GCOMTim.hpp"
+#include "GCOMStatus.hpp"
+#include "GCOMObservation.hpp"
+#include "GCOMSelection.hpp"
  
 /* __ Method name definitions ____________________________________________ */
 #define G_AT                                             "GCOMDris::at(int&)"
 #define G_INSERT                           "GCOMDris::insert(int&, GCOMDri&)"
 #define G_REMOVE                                     "GCOMDris::remove(int&)"
+#define G_COMPUTE_DRWS            "GCOMDris::compute_drws(GCOMObservation&, "\
+                                                   "GCOMSelection&, double&)"
  
 /* __ Macros _____________________________________________________________ */
  
 /* __ Coding definitions _________________________________________________ */
  
 /* __ Debug definitions __________________________________________________ */
+#define G_DEBUG_COMPUTE_DRWS
  
  
 /*==========================================================================
@@ -328,6 +338,303 @@ void GCOMDris::extend(const GCOMDris&amp; dris)
  
  
 /***********************************************************************//**
+ * @brief Compute backgroud weighting cubes
+ *
+ * @param[in] obs COMPTEL observation.
+ * @param[in] select Selection set.
+ * @param[in] zetamin Minimum Earth horizon - Phibar cut (deg).
+ *
+ * @exception GException::invalid_value
+ *            DRW with mismatching definition encountered in container.
+ * @exception GException::invalid_argument
+ *            No event list found in COMPTEL observation.
+ *
+ * Compute DRW cubes for a COMPTEL observation.
+ ***************************************************************************/
+void GCOMDris::compute_drws(const GCOMObservation& obs,
+                            const GCOMSelection&   select,
+                            const double&          zetamin)
+{
+    // Debug
+    #if defined(G_DEBUG_COMPUTE_DRWS)
+    std::cout << "GCOMDris::compute_drws" << std::endl;
+    std::cout << "======================" << std::endl;
+    #endif
+
+    // Set constants
+    const double dtime = 5.0 * 60.0; // Event rate integration time in seconds
+
+    // Continue only if there are DRWs in container
+    if (!is_empty()) {
+
+        // Initialise D1 & D2 module status
+        GCOMStatus status;
+
+        // Get pointer to first DRW in container (for convenient access to
+        // members of GCOMDri; we won't use it for modifications, hence we
+        // can declare it as constant)
+        const GCOMDri* dri = &(m_dris[0]);
+
+        // Initialise all DRWs
+        for (int i = 0; i < size(); ++i) {
+
+            // Throw an exception if the DRW definition differs
+            if (dri->m_dri != m_dris[i].m_dri) {
+                std::string msg = "Mismatch between definition of DRW "+
+                                  gammalib::str(i)+" and the first DRW. "
+                                  "Method only works on DRWs with identical "
+                                  "definitions.";
+                throw GException::invalid_value(G_COMPUTE_DRWS, msg);
+            }
+
+            // Initialise superpacket statistics
+            m_dris[i].init_statistics();
+
+            // Store selection set so that handling of failed PMT flag is
+            // correctly ritten into the FITS header
+            m_dris[i].m_selection = select;
+
+            // Set all DRG bins to zero
+            m_dris[i].init_cube();
+
+        } // endfor: looped over all DRWs
+
+        // Get pointer to event list. Throw an exception if the observation
+        // did not contain an event list
+        const GCOMEventList* events = dynamic_cast<const GCOMEventList*>(obs.events());
+        if (events == NULL) {
+            std::string msg = "No event list found in COMPTEL observation. Please "
+                              "specify an observation that contains an event list.";
+            throw GException::invalid_argument(G_COMPUTE_DRWS, msg);
+        }
+
+        // Initialise current time and vectors for computation of node function
+        // and vectors for energy-dependent rate interpolation
+        double              time;
+        double              rate0 = 0.0;
+        double              rate1 = 0.0;
+        double              rate2 = 0.0;
+        GNodeArray          times;
+        std::vector<double> rates0;
+        std::vector<double> rates1;
+        std::vector<double> rates2;
+
+        // Set time of first event
+        if (events->size() > 0) {
+            time = (*events)[0]->time().secs();
+        }
+
+        // Determine energy-dependent event rates. The standard event selection
+        // criteria and earth horizon cut is applied so that the event rate
+        // corresponds to the current selection, yet we do not consider any
+        // Good Time Intervals or OAD time information as we later only access
+        // relevant times (YET THERE COULD BE AN ISSUE WITH A GLITCH). Rates
+        // are computed for three hard-coded energy intervals.
+        for (int i_evt = 0; i_evt < events->size(); ++i_evt) {
+
+            // Get pointer to event
+            const GCOMEventAtom* event = (*events)[i_evt];
+
+            // Apply event selection
+            if (!select.use_event(*event)) {
+                continue;
+            }
+
+            // Apply Earth horizon cut
+            if ((event->eha() - event->phibar()) < zetamin) {
+                continue;
+            }
+
+            // If time exceeds time bin then add rates
+            while (event->time().secs() > (time + dtime)) {
+
+                // If rates are not all zero then append rates to vectors
+                if ((rate0 != 0.0) || (rate1 != 0.0) || (rate2 != 0.0)) {
+
+                    // Append time at mid-point of interval
+                    times.append(time + 0.5 * dtime);
+
+                    // Append rates
+                    rates0.push_back(rate0);
+                    rates1.push_back(rate1);
+                    rates2.push_back(rate2);
+
+                    // Debug
+                    #if defined(G_DEBUG_COMPUTE_DRWS)
+                    std::cout << "t=" << time + 0.5 * dtime;
+                    std::cout << " r0=" << rate0;
+                    std::cout << " r1=" << rate1;
+                    std::cout << " r2=" << rate2;
+                    std::cout << std::endl;
+                    #endif
+
+                    // Initialise rates for new accumulation
+                    rate0 = 0.0;
+                    rate1 = 0.0;
+                    rate2 = 0.0;
+
+                } // endif: there were rates to append
+
+                // Increment time
+                time += dtime;
+
+            } // endwhile: event time was after the current integration time stop
+
+            // Accumulate rates
+            double energy = event->energy().MeV();
+            if (energy < 1.8) {
+                rate0 += 1.0;
+            }
+            else if (energy < 4.3) {
+                rate1 += 1.0;
+            }
+            else {
+                rate2 += 1.0;
+            }
+
+        } // endfor: looped over events
+
+        // Get Good Time Intervals. If a pulsar selection is specified then
+        // reduce the Good Time Intervals to the validity intervals of the
+        // pulsar emphemerides.
+        GCOMTim tim = obs.tim();
+        if (select.has_pulsar()) {
+            tim.reduce(select.pulsar().validity());
+        }
+
+        // Loop over Orbit Aspect Data
+        for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
+
+            // Get reference to Orbit Aspect Data of superpacket
+            const GCOMOad &oad = obs.oads()[i_oad];
+
+            // Check superpacket usage for all DRWs so that their statistics
+            // get correctly updated
+            bool skip = false;
+            for (int i = 0; i < size(); ++i) {
+                if (!m_dris[i].use_superpacket(oad, tim, select)) {
+                    skip = true;
+                }
+            }
+            if (skip) {
+                continue;
+            }
+
+            // Get Orbit Aspect Data mid-time in seconds
+            double time = oad.tstart().secs() + 0.5 * (oad.tstop() - oad.tstart());
+
+            // Prepare Earth horizon angle computation. The celestial system
+            // is reinterpreted as the COMPTEL coordinate system, where the
+            // 90 degrees - zenith angle becomes the declination and the azimuth
+            // angle becomes Right Ascension. This allows us later to use the
+            // GSkyDir::dist method to compute the distance between the geocentre
+            // and the telescope Z-axis.
+            double  theta_geocentre = double(oad.gcel());
+            double  phi_geocentre   = double(oad.gcaz());
+            GSkyDir geocentre_comptel;
+            geocentre_comptel.radec_deg(phi_geocentre, 90.0-theta_geocentre);
+
+            // Compute geometrical factors for all (Chi, Psi). We assume here
+            // that their size and definition is identical, which checked on
+            // entry
+            for (int index = 0; index < dri->m_dri.npix(); ++index) {
+
+                // Get sky direction for (Chi, Psi)
+                GSkyDir sky = dri->m_dri.inx2dir(index);
+
+                // Convert sky direction to COMPTEL coordinates
+                double theta = oad.theta(sky);
+                double phi   = oad.phi(sky);
+
+                // Compute geometric factor summed over all D1, D2
+                double geometry = dri->compute_geometry(oad.tjd(), theta, phi, select, status);
+
+                // Compute Earth horizon angle as the distance between the Earth
+                // centre in COMPTEL coordinates and the (Chi, Psi) pixel in
+                // COMPTEL coordinates minus the radius of the Earth.
+                GSkyDir chipsi_comptel;
+                chipsi_comptel.radec_deg(phi, 90.0-theta);
+                double eha = geocentre_comptel.dist_deg(chipsi_comptel) - oad.georad();
+
+                // Loop over Phibar
+                for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
+
+                    // Compute minimum Earth horizon angle for Phibar layer
+                    double ehamin = double(iphibar) * dri->m_phibin + zetamin;
+
+                    // Add up geometry if the Earth horizon angle is equal or
+                    // larger than the minimum
+                    if (eha >= ehamin) {
+
+                        // Get data space index
+                        int inx = index + iphibar * dri->m_dri.npix();
+
+                        // Loop over DRWs and multiply with appopriate rates
+                        for (int i = 0; i < size(); ++i) {
+
+                            // Set rate by interpolating the precomputed vectors
+                            double rate = 0.0;
+                            if (m_dris[i].m_ebounds.emin().MeV() > 4.3) {
+                                rate = times.interpolate(time, rates2);
+                            }
+                            else if (m_dris[i].m_ebounds.emin().MeV() > 1.8) {
+                                rate = times.interpolate(time, rates1);
+                            }
+                            else {
+                                rate = times.interpolate(time, rates0);
+                            }
+
+                            // Set DRW
+                            m_dris[i][inx] += geometry * rate;
+
+                        } // endfor: looped over DRWs
+
+                    } // endif: Earth horizon cut passed
+
+                } // endfor: looped over Phibar
+
+            } // endfor: looped over (Chi, Psi)
+
+        } // endfor: looped over Orbit Aspect Data
+
+        // Normalisation Phibar layers of all DRWs to unity
+        for (int i = 0; i < size(); ++i) {
+
+            // Loop over Phibar layers
+            for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
+
+                // Compute sum in Phibar layer
+                double sum    = 0.0;
+                int    istart = iphibar * dri->m_dri.npix();
+                int    istop  = istart  + dri->m_dri.npix();
+                for (int inx = istart; inx < istop; ++inx) {
+                    sum += m_dris[i][inx];
+                }
+
+                // Normalise by sum
+                if (sum != 0.0) {
+                    for (int inx = istart; inx < istop; ++inx) {
+                        m_dris[i][inx] /= sum;
+                    }
+                }
+
+            } // endfor: looped over Phibar layers
+
+            // Set the Good Time interval for DRWs
+            if (m_dris[i].m_num_used_superpackets > 0) {
+                m_dris[i].m_gti = GGti(dri->m_tstart, dri->m_tstop);
+            }
+
+        } // endfor: looped over DRWs
+
+    } // endif: container was not empty
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Print Data Space container
  *
  * @param[in] chatter Chattiness.