From 8fc4e141a422254ee4bdd05ae677d35d762d4f35 Mon Sep 17 00:00:00 2001 From: L. Tibaldo Date: Fri, 15 Jun 2018 12:04:47 +0200 Subject: [PATCH] Modified as per issue description --- src/ctselect/ctselect.cpp | 46 ++++++++++++++++++++++++++++++++++++++++++---- src/ctselect/ctselect.hpp | 1 + src/ctselect/ctselect.par | 5 +++-- 3 files changed, 46 insertions(+), 6 deletions(-) diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 242cbfa..49a3e68 100755 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -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& 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)); diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index d76ce9b..7256c22 100755 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -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 m_infiles; //!< Input event filenames diff --git a/src/ctselect/ctselect.par b/src/ctselect/ctselect.par index 2e4d1b8..0ce9387 100755 --- a/src/ctselect/ctselect.par +++ b/src/ctselect/ctselect.par @@ -23,10 +23,11 @@ prefix, s, h, selected_,,, "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)" -- libgit2 0.21.4