Add ctool::get_roi() method (#1828)

The ctool::get_roi() replaces some duplicated code in ctselect and ctcubemask. The method extracts the RoI from the user parameters "ra", "dec" and "rad". In case the "usepnt" parameter is present and set to "yes", the "ra" and "dec" parameters will not be read, and it's up to the client to set the RoI centre properly.

Add ctool::get_roi() method (#1828)
The ctool::get_roi() replaces some duplicated code in ctselect and ctcubemask. The method extracts the RoI from the user parameters "ra", "dec" and "rad". In case the "usepnt" parameter is present and set to "yes", the "ra" and "dec" parameters will not be read, and it's up to the client to set the RoI centre properly.
Jürgen Knödlseder
1 parent f9a30a18
Showing 7 changed files with 85 additions and 68 deletions
pyext/ctool.i
src/ctcubemask/ctcubemask.cpp
src/ctcubemask/ctcubemask.hpp
src/ctselect/ctselect.cpp
src/ctselect/ctselect.hpp
src/support/ctool.cpp
src/support/ctool.hpp
@@ -97,6 +97,7 @@ public:
     %rename(_require_inobs_nocube)     require_inobs_nocube;
     %rename(_log_observations)         log_observations;
     %rename(_log_models)               log_models;
+    %rename(_get_roi)                  get_roi;
     %rename(_set_response)             set_response;
     %rename(_set_edisp)                set_edisp;
     %rename(_restore_edisp)            restore_edisp;
@@ -140,6 +141,9 @@ public:
                                const GModels&     models,
                                const std::string& what = "Model");
+    // Protected methods that extract user parameters
+    GCTARoi           get_roi(void);
+
     // Protected support methods
     void              set_response(GObservations& obs);
     std::vector<bool> set_edisp(GObservations& obs, const bool& edisp) const;
@@ -377,9 +377,7 @@ void ctcubemask::init_members(void)
     m_outcube.clear();
 	m_prefix.clear();
     m_usepnt  = false;
-    m_ra      = 0.0;
-    m_dec     = 0.0;
-    m_rad     = 0.0;
+    m_roi.clear();
     m_emin    = 0.0;
     m_emax    = 0.0;
     m_publish = false;
@@ -407,9 +405,7 @@ void ctcubemask::copy_members(const ctcubemask&amp; app)
     m_outcube = app.m_outcube;
     m_prefix  = app.m_prefix;
     m_usepnt  = app.m_usepnt;
-    m_ra      = app.m_ra;
-    m_dec     = app.m_dec;
-    m_rad     = app.m_rad;
+    m_roi     = app.m_roi;
     m_emin    = app.m_emin;
     m_emax    = app.m_emax;
     m_publish = app.m_publish;
@@ -458,26 +454,10 @@ void ctcubemask::get_parameters(void)
     // Get parameters
     m_regfile = (*this)["regfile"].filename();
     m_usepnt  = (*this)["usepnt"].boolean();
-    if (!m_usepnt) {
-        // Check RA/DEC parameters for validity to read
-        if ((*this)["ra"].is_valid() && (*this)["dec"].is_valid()) {
-           m_ra         = (*this)["ra"].real();
-           m_dec        = (*this)["dec"].real();
-           m_select_roi = true;
-        }
-        else {
-           m_select_roi = false;
-        }
-    }
-
-    // Check if radius is valid for a RoI selection
-    if (m_select_roi && (*this)["rad"].is_valid()) {
-       m_rad = (*this)["rad"].real();
-    }
-    else {
-       m_select_roi = false;
-    }
+    // Get the RoI and enable RoI selection if the RoI is valid
+    m_roi        = get_roi();
+    m_select_roi = m_roi.is_valid();
     // Check for sanity of energy selection parameters
     if ((*this)["emin"].is_valid() && (*this)["emax"].is_valid()) {
@@ -587,7 +567,7 @@ void ctcubemask::apply_mask(GCTAObservation* obs)
         // Set all pixels inside selected energy bands but outside RoI
         // to -1.0 if requested
         if (m_select_roi) {
-            GSkyRegionCircle roi(m_ra, m_dec, m_rad);
+            GSkyRegionCircle roi(m_roi.centre().dir(), m_roi.radius());
             for (int i = e_idx1; i <= e_idx2; ++i) {
                 for (int pixel = 0; pixel < npix; ++pixel) {
                     GSkyDir dir = map.inx2dir(pixel);
@@ -599,9 +579,9 @@ void ctcubemask::apply_mask(GCTAObservation* obs)
             // Log selected RoI
             log_value(NORMAL, "Selected RoI",
-                              "RA="+gammalib::str(m_ra)+" deg, "+
-                              "DEC="+gammalib::str(m_dec)+" deg, "+
-                              "Radius="+gammalib::str(m_rad)+" deg");
+                      "RA="+gammalib::str(roi.centre().ra_deg())+" deg, "+
+                      "DEC="+gammalib::str(roi.centre().dec_deg())+" deg, "+
+                      "Radius="+gammalib::str(roi.radius())+" deg");
         } // endif: applied RoI selection
@@ -81,9 +81,7 @@ protected:
     GFilename   m_outcube;    //!< Output event list or XML file
 	std::string m_prefix;     //!< Prefix for multiple counts maps
     bool        m_usepnt;     //!< Use pointing instead of RA/DEC parameters
-    double      m_ra;         //!< RA of ROI centre
-    double      m_dec;        //!< DEC of ROI centre
-    double      m_rad;        //!< ROI radius
+    GCTARoi     m_roi;        //!< RoI selection
     double      m_emin;       //!< Lower energy
     double      m_emax;       //!< Upper energy
     bool        m_publish;    //!< Publish counts cube?
@@ -408,9 +408,7 @@ void ctselect::init_members(void)
     m_outobs.clear();
     m_prefix.clear();
     m_usepnt = false;
-    m_ra     = -1.0;
-    m_dec    = -1.0;
-    m_rad    = -1.0;
+    m_roi.clear();
     m_tmin   = 0.0;
     m_tmax   = 0.0;
     m_emin   = 0.0;
@@ -446,9 +444,7 @@ void ctselect::copy_members(const ctselect&amp; app)
     m_outobs   = app.m_outobs;
     m_prefix   = app.m_prefix;
     m_usepnt   = app.m_usepnt;
-    m_ra       = app.m_ra;
-    m_dec      = app.m_dec;
-    m_rad      = app.m_rad;
+    m_roi      = app.m_roi;
     m_tmin     = app.m_tmin;
     m_tmax     = app.m_tmax;
     m_emin     = app.m_emin;
@@ -504,28 +500,12 @@ void ctselect::get_parameters(void)
     // information and do not accept counts cubes.
     setup_observations(m_obs, false, true, false);
-    // Get parameters
+    // Get some parameters
     m_usepnt = (*this)["usepnt"].boolean();
-    if (!m_usepnt) {
-        // Check RA/DEC parameters for validity to read
-        if ((*this)["ra"].is_valid() && (*this)["dec"].is_valid()) {
-            m_ra         = (*this)["ra"].real();
-            m_dec        = (*this)["dec"].real();
-            m_select_roi = true;
-        }
-        else {
-            m_select_roi = false;
-        }
-    }
-
-    // Check if radius is vaild for a RoI selection
-    if (m_select_roi && (*this)["rad"].is_valid()) {
-        m_rad = (*this)["rad"].real();
-    }
-    else {
-        m_select_roi = false;
-    }
+    // Get the RoI and enable RoI selection if the RoI is valid
+    m_roi        = get_roi();
+    m_select_roi = m_roi.is_valid();
     // Check for sanity of time selection parameters
     if ((*this)["tmin"].is_valid() && (*this)["tmax"].is_valid()) {
@@ -640,11 +620,11 @@ void ctselect::select_events(GCTAObservation*   obs,
     // Analyse expression to see if a selection is required
     bool select_expr = (gammalib::strip_whitespace(m_expr).length() > 0);
-    // Set RA/DEC selection. If the "usepnt" parameter is set to true then
-    // use the pointing direction as the ROI centre.
-    double ra  = m_ra;
-    double dec = m_dec;
-    double rad = m_rad;
+    // Set RoI selection. If the "usepnt" parameter is set to true then
+    // use the pointing direction as the RoI centre.
+    double ra  = m_roi.centre().dir().ra_deg();
+    double dec = m_roi.centre().dir().dec_deg();
+    double rad = m_roi.radius();
     if (m_usepnt) {
         const GCTAPointing &pnt = obs->pointing();
         ra  = pnt.dir().ra_deg();
@@ -751,7 +731,7 @@ void ctselect::select_events(GCTAObservation*   obs,
         // Log the requested RoI
         log_value(NORMAL, "Requested RoI",
                   "Centre(RA,DEC)=("+gammalib::str(ra)+", "+
-                  gammalib::str(dec)+") deg, Radius="+gammalib::str(m_rad)+
+                  gammalib::str(dec)+") deg, Radius="+gammalib::str(rad)+
                   " deg");
         // If we have already an RoI then make sure that the selected
@@ -92,9 +92,7 @@ protected:
     std::string m_outobs;     //!< Output event list or XML file
     std::string m_prefix;     //!< Prefix for multiple event lists
     bool        m_usepnt;     //!< Use pointing instead of RA/DEC parameters
-    double      m_ra;         //!< RA of ROI centre
-    double      m_dec;        //!< DEC of ROI centre
-    double      m_rad;        //!< ROI radius
+    GCTARoi     m_roi;        //!< RoI selection
     double      m_tmin;       //!< Start time
     double      m_tmax;       //!< Stop time
     double      m_emin;       //!< Lower energy
@@ -111,7 +109,7 @@ protected:
     GTime                    m_timemin;       //!< Earliest time
     GTime                    m_timemax;       //!< Latest time
     bool                     m_select_energy; //!< Perform energy selection
-    bool                     m_select_roi;    //!< Perform ROI selection
+    bool                     m_select_roi;    //!< Perform RoI selection
     bool                     m_select_time;   //!< Perform time selection
 };
@@ -997,6 +997,60 @@ void ctool::log_models(const GChatter&amp;    chatter,
 /***********************************************************************//**
+ * @brief Return RoI from user parameters
+ *
+ * @param[in] obs Observation container.
+ * @return Region of interest.
+ *
+ * Returns region of interest from the user parameters "ra", "dec" and "rad".
+ * If the "usepnt" parameter exists and is "yes" then postpone then assume
+ * that the Right Ascension and Declination of the RoI centre will be set
+ * from the pointing direction somewhere else.
+ ***************************************************************************/
+GCTARoi ctool::get_roi(void)
+{
+    // Initialise an empty RoI
+    GCTARoi roi;
+
+    // Signal whether the pointing instead of the user parameters should be
+    // use for the RoI centre
+    bool usepnt = (has_par("usepnt") && (*this)["usepnt"].boolean());
+
+    // If pointing should not bet used the extract the RoI centre from the
+    // "ra" and "dec" parameters
+    if (!usepnt) {
+
+        // Read "ra" and "dec" parameters only if both are valid
+        if ((*this)["ra"].is_valid() && (*this)["dec"].is_valid()) {
+
+            // Read parameters
+            double ra  = (*this)["ra"].real();
+            double dec = (*this)["dec"].real();
+
+            // Set RoI
+            GCTAInstDir instdir;
+            instdir.dir().radec_deg(ra, dec);
+            roi.centre(instdir);
+
+            // Signal that we have a pointing direction
+            usepnt = true;
+
+        } // endif: "ra" and "dec" were valid
+
+    } // endif: pointing was not used for RoI centre
+
+    // If we have a RoI centre and if the "rad" parameter is valid then
+    // extract the RoI radius from the parameter
+    if (usepnt && (*this)["rad"].is_valid()) {
+        roi.radius((*this)["rad"].real());
+    }
+
+    // Return RoI
+    return roi;
+}
+
+
+/***********************************************************************//**
  * @brief Set response for all CTA observations in container
  *
  * @param[in,out] obs Observation container
@@ -103,6 +103,9 @@ protected:
                                const GModels&     models,
                                const std::string& what = "Model");
+    // Protected methods that extract user parameters
+    GCTARoi           get_roi(void);
+
     // Protected support methods
     void              set_response(GObservations& obs);
     std::vector<bool> set_edisp(GObservations& obs, const bool& edisp) const;