Modified as per issue description

L. Tibaldo · Jürgen Knödlseder
1 parent 035428e5
Showing 3 changed files with 46 additions and 6 deletions
src/ctselect/ctselect.cpp
src/ctselect/ctselect.hpp
src/ctselect/ctselect.par
@@ -478,6 +478,7 @@ void ctselect::init_members(void)
     m_phases.clear();
     m_select_energy = false;
     m_select_phase  = false;
+    m_forcesel = false;
  
     // Return
     return;
@@ -506,6 +507,7 @@ void ctselect::copy_members(const ctselect&amp; app)
     m_phases        = app.m_phases;
     m_select_energy = app.m_select_energy;
     m_select_phase  = app.m_select_phase;
+    m_forcesel  = app.m_forcesel;
  
     // Return
     return;
@@ -548,6 +550,9 @@ void ctselect::get_parameters(void)
     // Query the RoI
     GCTAPointing pnt;
     GCTARoi      roi = get_roi(pnt);
+    // Signal whether the RoI selection should be enforced
+    // even when leading to invalid RoI
+    m_forcesel = (has_par("forcesel") && (*this)["forcesel"].boolean());
  
     // Query Good Time Intervals. The GammaLib time reference is specified
     // as a dummy argument since the relevant intervals will be queried later
@@ -730,6 +735,9 @@ void ctselect::select_events(GCTAObservation*   obs,
     // Signal that RoI selection should be performed
     bool select_roi = roi.is_valid();
  
+    //Create boolean to remember if RoI selection is enforced regardless of its validity
+    bool enforced_roi = false;
+
     // Set time selection interval. We make sure here that the time selection
     // interval cannot be wider than the GTIs covering the data. This is done
     // using GGti's reduce() method.
@@ -913,8 +921,38 @@ void ctselect::select_events(GCTAObservation*   obs,
             GSkyDir centre;
             centre.radec_deg(ra, dec);
             double distance = centre.dist_deg(roi_centre);
-            if (distance + rad > roi_radius) {
-                rad = roi_radius - distance;
+            // If the requested RoI and selected RoI are centred at the same position
+            // and the requested RoI is bigger then the original
+            // keep original radius
+            // RoIs same centre = distance < 1.e-4
+            // to cope with numerical precision limitations
+            if (distance < 1.e-4 && rad > roi_radius) {
+                rad = roi_radius;
+            }
+            //Otherwise, check if requested RoI is enclosed within original RoI
+            else {
+                if (distance + rad > roi_radius) {
+                    //If selection is enforced log warning and store
+                    if (m_forcesel){
+                    enforced_roi = true;
+                    log_string(TERSE, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
+                    log_string(TERSE,"Enforced RoI selection leading to inconsistent RoI. "\
+                                "RoI of data: Centre(RA,DEC)=("+gammalib::str(ra)+", "+
+                                gammalib::str(dec)+") deg, Radius="+gammalib::str(rad)+
+                                " deg. Requested RoI: Centre(RA,DEC)=("+
+                                gammalib::str(roi_centre.ra_deg())+", "+
+                                gammalib::str(roi_centre.dec_deg())+") deg, Radius="+
+                                gammalib::str(roi_radius)+" deg. The resulting event list "\
+                                "file is not usable for likelihood analysis");
+                    log_string(TERSE, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
+                    }
+                    // Otherwise, throw an exception
+                    else {
+                    std::string message;
+                    message = "Invalid RoI selection: the new RoI must be enclosed in the original RoI";
+                    throw GException::invalid_value(G_RUN, message);
+                    }
+                }
             }
         }
  
@@ -1011,8 +1049,8 @@ void ctselect::select_events(GCTAObservation*   obs,
     // closed later
     list->fetch();
  
-    // If RoI selection has been applied then set the event list RoI
-    if (select_roi) {
+    // If RoI selection has been applied without forcing then set the event list RoI
+    if (select_roi && !enforced_roi) {
         GCTAInstDir instdir;
         instdir.dir().radec_deg(ra, dec);
         list->roi(GCTARoi(instdir, rad));
@@ -87,6 +87,7 @@ protected:
     std::string m_expr;       //!< Selection expression
     std::string m_usethres;   //!< Energy threshold type
     GChatter    m_chatter;    //!< Chattiness
+    bool        m_forcesel;   //!< Enforce RoI selection
  
     // Protected members
     std::vector<std::string> m_infiles;       //!< Input event filenames
@@ -23,10 +23,11 @@ prefix, s, h, selected_,,, &quot;Prefix for event lists in observation definition XML
 #
 # Selection parameters
 #=====================
-usepnt,  b, h, no,,, "Use pointing instead of RA/DEC parameters?"
+usepnt,  b, h, yes,,, "Use pointing instead of RA/DEC parameters?"
 ra,      r, a, 83.63,0,360, "RA for ROI centre (degrees)"
 dec,     r, a, 22.51,-90,90, "Dec for ROI centre (degrees)"
-rad,     r, a, 3.0,0,180, "Radius of ROI (degrees)"
+rad,     r, a, 3.0,0,180, "Radius of ROI around pointing or specified RA/DEC (degrees)"
+forcesel,b, h, no,,, "Force RoI selection even if it produces invalid RoI"
 tmin,    t, a, NONE,,, "Start time (UTC string, JD, MJD or MET in seconds)"
 tmax,    t, a, NONE,,, "Stop time (UTC string, JD, MJD or MET in seconds)"
 emin,    r, a, 0.1,,, "Lower energy limit (TeV)"