Integrate code.

Jürgen Knödlseder
1 parent 8e535189
Showing 5 changed files with 254 additions and 172 deletions
ChangeLog
NEWS
src/ctselect/ctselect.cpp
src/ctselect/ctselect.hpp
src/ctselect/ctselect.par
-2014-10-30  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-11-08  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
  
         * Version 0.8.0 released
         ========================
@@ -17,6 +17,7 @@
         ctmodel now always creates a single model cube
         Support energy dispersion in ctobssim, ctlike and ctmodel	
         Correct ctobssim for thread safety
+        Support save and user energy thresholds in ctselect
         Correctly apply time references in ctselect
         Reduce memory usage
         Move obsutils to ctools namespace
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
  
 Author(s) : Juergen Knoedlseder
  
-30 October 2014
+8 November 2014
  
  
 1. Introduction
@@ -65,6 +65,30 @@ tmin and tmax parameters are given in CTA reference time, they are
 internally converted to the reference time of the GTI before event
 selection (#1329).
  
+ctselect now supports save and user energy thresholds (#1039). There
+are three different possibilities for specifying a threshold. They
+are toggled using the hidden usethres parameter:
+
+NONE: Use the emin and emax parameters to select an energy range.
+      This is the default mode which is identical to the previous
+      workings. The only difference is now that if the emin and/or
+      emax parameter is set to 0.0, the threshold will be ignored.
+      This allows to select all events above or below a given
+      threshold.
+
+DEFAULT: Use the save threshold that is found in the effective
+      area FITS header keywords LO_THRES and HI_THRES. If the
+      emin and emax parameters are non-zero, they will provide
+      an additional energy constraint.
+
+USER: Use the user specified threshold that are found in the
+      XML observation definition file. If the emin and emax
+      parameters are non-zero, they will provide an additional
+      energy constraint.
+
+Note that any existing energy boundaries will be used as envelope,
+i.e. a wider energy threshold specification will have no effect.
+
  
 7. ctbin - CTA event binning
 ----------------------------
@@ -36,9 +36,8 @@
 #define G_RUN                                               "ctselect::run()"
 #define G_SELECT_EVENTS          "ctselect::select_events(GCTAObservation*, "\
                                                               "std::string&)"
-#define G_FIND_EBOUNDS          "ctselect::find_ebounds(GCTAObservation*, "\
-                                                              "GEbounds&)"
-#define G_GET_PARAMETERS                      "ctselect::get_parameters()"
+#define G_SET_EBOUNDS    "ctselect::set_ebounds(GCTAObservation*, GEbounds&)"
+#define G_GET_PARAMETERS                         "ctselect::get_parameters()"
  
 /* __ Debug definitions __________________________________________________ */
  
@@ -299,7 +298,6 @@ void ctselect::run(void)
             // Load observation from temporary file, including event selection
             select_events(obs, filename);
  
-
             // Close temporary file
             #if G_USE_MKSTEMP
             close(fileid);
@@ -408,7 +406,7 @@ void ctselect::init_members(void)
     m_infiles.clear();
     m_timemin.clear();
     m_timemax.clear();
-    m_use_xml    = false;
+    m_use_xml = false;
  
     // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD,
     // which is defined in GCTALib.hpp. This is somehow a kluge. We need
@@ -428,18 +426,18 @@ void ctselect::init_members(void)
 void ctselect::copy_members(const ctselect& app)
 {
     // Copy parameters
-    m_infile  = app.m_infile;
-    m_outfile = app.m_outfile;
-    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_tmin    = app.m_tmin;
-    m_tmax    = app.m_tmax;
-    m_emin    = app.m_emin;
-    m_emax    = app.m_emax;
-    m_expr    = app.m_expr;
+    m_infile   = app.m_infile;
+    m_outfile  = app.m_outfile;
+    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_tmin     = app.m_tmin;
+    m_tmax     = app.m_tmax;
+    m_emin     = app.m_emin;
+    m_emax     = app.m_emax;
+    m_expr     = app.m_expr;
     m_usethres = app.m_usethres;
  
     // Copy protected members
@@ -520,15 +518,19 @@ void ctselect::get_parameters(void)
         m_ra  = (*this)["ra"].real();
         m_dec = (*this)["dec"].real();
     }
-    m_rad  = (*this)["rad"].real();
-    m_tmin = (*this)["tmin"].real();
-    m_tmax = (*this)["tmax"].real();
-    m_emin = (*this)["emin"].real();
-    m_emax = (*this)["emax"].real();
-    m_expr = (*this)["expr"].string();
+    m_rad      = (*this)["rad"].real();
+    m_tmin     = (*this)["tmin"].real();
+    m_tmax     = (*this)["tmax"].real();
+    m_emin     = (*this)["emin"].real();
+    m_emax     = (*this)["emax"].real();
+    m_expr     = (*this)["expr"].string();
     m_usethres = (*this)["usethres"].string();
+
+    // Check parameters (Kluge as the parfile interface does not yet check
+    // parameters; this can be removed once this is done properly)
     if (m_usethres != "NONE" && m_usethres != "DEFAULT" && m_usethres != "USER") {
-        std::string msg = "Parameter \"usethres\" must be one of \"NONE\", \"DEFAULT\" or \"USER\"";
+        std::string msg = "Parameter \"usethres\" must be one of \"NONE\", "
+                          "\"DEFAULT\" or \"USER\".";
         throw GException::invalid_value(G_GET_PARAMETERS, msg);
     }
  
@@ -574,23 +576,31 @@ void ctselect::get_parameters(void)
  ***************************************************************************/
 void ctselect::select_events(GCTAObservation* obs, const std::string& filename)
 {
-
     // Initialise selection string
     std::string selection;
  
     // Initialise selection flags
     bool remove_all = false;
  
-    // Analyse task parameters to see what selections are required
+    // Get CTA event list pointer
+    GCTAEventList* list =
+        static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
+
+    // Get energy boundaries
+    double   emin    = 0.0;
+    double   emax    = 0.0;
+    GEbounds ebounds = set_ebounds(obs, list->ebounds());
+    if (ebounds.size() >= 1) {
+        emin = ebounds.emin(0).TeV();
+        emax = ebounds.emax(0).TeV();
+    }
+
+    // Analyse parameters to see what selections are required
     bool select_time   = (m_tmin != 0.0 || m_tmax != 0.0);
-    bool select_energy = true;
+    bool select_energy = (emin != 0.0 || emax != 0.0);
     bool select_roi    = ((m_dec != -1.0 && m_ra != -1.0 && m_rad != -1.0) || m_usepnt);
     bool select_expr   = (gammalib::strip_whitespace(m_expr).length() > 0);
  
-    // Set energy selection
-    double emin = m_emin;
-    double emax = m_emax;
-    
     // 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;
@@ -601,10 +611,6 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
         ra  = pnt.dir().ra_deg();
         dec = pnt.dir().dec_deg();
     }
-
-    // Get CTA event list pointer
-    GCTAEventList* list =
-        static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
  
     // 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
@@ -650,55 +656,58 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
             selection += " && ";
         }
  
-        // Log the requested energy range
+        // Log the requested energy selection
         if (logTerse()) {
-            log << gammalib::parformat("Requested energy range");
-            log << m_emin << " - " << m_emax << " TeV" << std::endl;
-        }
-
-        // Set the energy boundaries. Make sure that they overlap with any
-        // existing energy boundaries
-
-        // Find the energy boundaries, taking into account any previous selections
-        // and also user or default energy thresholds
-        GEbounds ebounds = find_ebounds(obs,list->ebounds());
-
-        // Check if energy bounds are valid
-        if (ebounds.size() >= 1) {
-
-            // Write ebounds to log
-            if (logTerse()) {
-                log << emin << " - " << emax << " TeV" << std::endl;
-            }
-
-            // Set emin and emax values
-            emin = ebounds.emin(0).TeV();
-            emax = ebounds.emax(0).TeV();
-
-        } //endif: ebounds were valid
-
-        else {
-
-            // Write log that energy range is not existent
-            if (logTerse()) {
-                log << gammalib::parformat("Selected energy range");
+            log << gammalib::parformat("Selected energy range");
+            if (emin > 0.0 && emax > 0.0) {
                 if (emin >= emax) {
                     log << "None. There is no overlap between existing ";
                     log << "and requested energy range." << std::endl;
                 }
+                else {
+                    log << emin << " - " << emax << " TeV" << std::endl;
+                }
+            }
+            else if (emin > 0.0) {
+                log << "> " << emin << " TeV" << std::endl;
+            }
+            else if (emax > 0.0) {
+                log << "< " << emax << " TeV" << std::endl;
             }
+            else {
+                log << "None" << std::endl;
+            }
+        }
  
-            // Signal to remove all events
+        // Apply the requested energy selection
+        if (emin > 0.0 && emax > 0.0) {
+            if (emin >= emax) {
+                remove_all = true;
+            }
+            else {
+                char cmin[80];
+                char cmax[80];
+                sprintf(cmin, "%.8f", emin);
+                sprintf(cmax, "%.8f", emax);
+                selection += "ENERGY >= "+std::string(cmin)+" && ENERGY <= "+
+                             std::string(cmax);
+            }
+        }
+        else if (emin > 0.0) {
+            char cmin[80];
+            sprintf(cmin, "%.8f", emin);
+            selection += "ENERGY >= "+std::string(cmin);
+        }
+        else if (emax > 0.0) {
+            char cmax[80];
+            sprintf(cmax, "%.8f", emax);
+            selection += "ENERGY <= "+std::string(cmax);
+        }
+        else {
+            log << "None" << std::endl;
             remove_all = true;
         }
  
-        // Format energy selection
-        char cmin[80];
-        char cmax[80];
-        sprintf(cmin, "%.8f", emin);
-        sprintf(cmax, "%.8f", emax);
-        selection += "ENERGY >= "+std::string(cmin)+" && ENERGY <= "+std::string(cmax);
-
     } // endif: made energy selection
  
     // Make ROI selection
@@ -860,6 +869,137 @@ void ctselect::select_events(GCTAObservation* obs, const std::string&amp; filename)
  
  
 /***********************************************************************//**
+ * @brief Return energy boundaries for a given observation
+ *
+ * @param[in] obs Pointer to CTA observation.
+ * @param[in] ebounds Current energy boundaries.
+ * @return Energy boundaries.
+ *
+ * Returns the energy boundaries for a given observation. Depending on the
+ * value of the usethres parameter, the following values will be
+ * returned:
+ *
+ *     NONE:    [max(emin,emin_exist),min(emax,emax_exist)]
+ *     DEFAULT: [max(emin,emin_exist,emin_save),min(emax,emax_exist,emax_save)]
+ *     USER:    [max(emin,emin_exist,emin_user),min(emax,emax_exist,emax_user)]
+ *
+ * where
+ *
+ *     emin is the value of the emin parameter
+ *     emax is the value of the emax parameter
+ *     emin_exist is the value of any existing minimum boundary
+ *     emax_exist is the value of any existing maximum boundary
+ *     emin_save is the lower save threshold
+ *     emax_save is the upper save threshold
+ *     emin_user is the lower user threshold
+ *     emax_user is the upper user threshold
+ *
+ * Any threshold value of 0 will be ignored.
+ ***************************************************************************/
+GEbounds ctselect::set_ebounds(GCTAObservation* obs, const GEbounds& ebounds) const
+{
+    // Set emin and emax values from user parameters
+    double emin = m_emin;
+    double emax = m_emax;
+
+    // Check if we have already energy boundaries
+    if (ebounds.size() > 0) {
+
+        // Raise minimum energy to lower boundary
+        if ((emin == 0.0) ||
+            ((ebounds.emin(0).TeV() > 0.0) && (emin < ebounds.emin(0).TeV()))) {
+            emin = ebounds.emin(0).TeV();
+        }
+
+        // Lower maximum energy to upper boundary
+        if ((emax == 0.0) ||
+            ((ebounds.emax(0).TeV() > 0.0) && (emax > ebounds.emax(0).TeV()))) {
+            emax = ebounds.emax(0).TeV();
+        }
+
+    } // endif: there were already energy boundaries
+
+    // Check if default threshold (the one from the IRF, also known as save
+    // threshold) should be applied
+    if (m_usethres == "DEFAULT") {
+
+        // Get CTA IRF respsonse pointer and throw exception if we don't
+        //  find IRF
+        const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
+        if (rsp == NULL) {
+            std::string msg = "No IRF response attached to given observation.";
+            throw GException::invalid_value(G_SET_EBOUNDS, msg);
+        }
+
+        // Retrieve energy range from response information
+        double lo_save_thres = rsp->lo_save_thres();
+        double hi_save_thres = rsp->hi_save_thres();
+
+        // Check threshold information for validity
+        if ((lo_save_thres > 0.0) &&
+            (hi_save_thres > 0.0) &&
+            (lo_save_thres >= hi_save_thres)) {
+            std::string msg = "IRF \""+obs->name()+"\" contains an invalid "
+                              "energy range ["+gammalib::str(lo_save_thres)+","+
+                              gammalib::str(hi_save_thres)+"] TeV.";
+            throw GException::invalid_value(G_SET_EBOUNDS, msg);
+        }
+
+        // Raise minimum energy to lower threshold (if lower threshold exists)
+        if ((emin == 0.0) ||
+            ((lo_save_thres > 0.0) && (emin < lo_save_thres))) {
+            emin = lo_save_thres;
+        }
+
+        // Lower maximum energy to upper threshold (if upper threshold exists)
+        if ((emax == 0.0) ||
+            ((hi_save_thres > 0.0) && (emax > hi_save_thres))) {
+            emax = hi_save_thres;
+        }
+
+    } //endif: usethres was "DEFAULT"
+
+    // ... otherwise check if a user specified threshold should be applied
+    else if (m_usethres == "USER") {
+
+        // Retrieve thresholds from observation
+        double lo_user_thres = obs->lo_user_thres();
+        double hi_user_thres = obs->hi_user_thres();
+
+        // Check threshold information for validity
+        if ((lo_user_thres > 0.0) &&
+            (hi_user_thres > 0.0) &&
+            (lo_user_thres >= hi_user_thres)) {
+            std::string msg = "User energy range ["+gammalib::str(lo_user_thres)+
+                              ","+gammalib::str(hi_user_thres)+"] TeV is invalid.";
+            throw GException::invalid_value(G_SET_EBOUNDS, msg);
+        }
+
+        // Raise minimum energy to lower threshold (if lower threshold exists)
+        if ((emin == 0.0) ||
+            ((lo_user_thres > 0.0) && (emin < lo_user_thres))) {
+            emin = lo_user_thres;
+        }
+
+        // Lower maximum energy to upper threshold (if upper threshold exists)
+        if ((emax == 0.0) ||
+            ((hi_user_thres > 0.0) && (emax > hi_user_thres))) {
+            emax = hi_user_thres;
+        }
+
+    } // endif: usethres was "USER
+
+    // Set selection energy boundaries
+    GEbounds result;
+    result.append(GEnergy(emin, "TeV"), GEnergy(emax, "TeV"));
+
+    // Return result
+    return result;
+
+}
+
+
+/***********************************************************************//**
  * @brief Check input filename
  *
  * @param[in] filename File name.
@@ -1079,87 +1219,3 @@ void ctselect::save_event_list(const GCTAObservation* obs,
     // Return
     return;
 }
-
-GEbounds ctselect::find_ebounds(GCTAObservation* obs, const GEbounds& list_ebounds) const
-{
-    // Initialise ebounds
-    GEbounds ebounds;
-
-    // Set emin and emax values from user parameters
-    double emin = m_emin;
-    double emax = m_emax;
-
-    // Check if list already has ebounds
-    if (list_ebounds.size() >= 1) {
-
-            // Make sure the user parameters don't condradict previous selections
-            if (emin < list_ebounds.emin(0).TeV()) {
-                emin = list_ebounds.emin(0).TeV();
-            }
-            if (emax > list_ebounds.emax(0).TeV()) {
-                emax = list_ebounds.emax(0).TeV();
-            }
-        } //endif: ebounds were valid
-
-    // Check if threshold should be applied on default
-    if (m_usethres == "DEFAULT") {
-
-        // Get CTA IRF respsonse pointer and throw exception
-        // if we don't find IRF
-        GCTAResponseIrf* rsp =  static_cast<GCTAResponseIrf*>(const_cast<GCTAResponse*>(obs->response()));
-        if (rsp == NULL) {
-            std::string msg = "No IRF response attached to given observation";
-            throw GException::invalid_value(G_FIND_EBOUNDS, msg);
-        }
-
-        // retrieve energy range from response information
-        double lo_save_thres = rsp->lo_save_thres();
-        double hi_save_thres = rsp->hi_save_thres();
-
-        // Check threshold information for validity
-        // throw exception if not valid
-        if (lo_save_thres >= hi_save_thres){
-            std::string msg = "Energy range provided by response from \""+obs->name()+"\" not valid";
-            throw GException::invalid_value(G_FIND_EBOUNDS, msg);
-        }
-
-        // Overwrite energy threshold definitions if necessary
-        if (emin < lo_save_thres) {
-            emin = lo_save_thres;
-        }
-        if (emax > hi_save_thres) {
-            emax = hi_save_thres;
-        }
-
-    } //endif: m_usethres was "DEFAULT"
-
-    else if (m_usethres == "USER") {
-
-        // Retrieve parameters from observation
-        // these should have been set via the xml file
-        double lo_user_thres = obs->lo_user_thres();
-        double hi_user_thres = obs->hi_user_thres();
-
-        // Check provided energy range for validity
-        if (lo_user_thres >= hi_user_thres) {
-                std::string msg = "User defined energy range which is provided in observation \""+obs->name()+"\" not valid";
-                throw GException::invalid_value(G_FIND_EBOUNDS, msg);
-        }
-
-        // Overwrite energy range definition if necessary
-        if (emin < lo_user_thres) {
-            emin = lo_user_thres;
-        }
-        if (emax > hi_user_thres) {
-            emax = hi_user_thres;
-        }
-
-    } // endif: usethres was "USER
-
-    // set ebounds for selection
-    ebounds.append(GEnergy(emin, "TeV"), GEnergy(emax, "TeV"));
-
-    // Return ebounds
-    return ebounds;
-
-}
@@ -71,6 +71,8 @@ protected:
     void        free_members(void);
     void        get_parameters(void);
     void        select_events(GCTAObservation* obs, const std::string& filename);
+    GEbounds    set_ebounds(GCTAObservation* obs,
+                            const GEbounds& ebounds) const;
     std::string check_infile(const std::string& filename) const;
     std::string set_outfile_name(const std::string& filename) const;
     void        save_fits(void);
@@ -78,7 +80,6 @@ protected:
     void        save_event_list(const GCTAObservation* obs,
                                 const std::string&     infile,
                                 const std::string&     outfile) const;
-    GEbounds find_ebounds(GCTAObservation* obs, const GEbounds& list_ebounds) const;
  
     // User parameters
     std::string m_infile;     //!< Input event list or XML file
@@ -93,7 +94,7 @@ protected:
     double      m_emin;       //!< Lower energy
     double      m_emax;       //!< Upper energy
     std::string m_expr;       //!< Selection expression
-    std::string m_usethres;//!<Choice of energy threshold usage
+    std::string m_usethres;   //!< Energy threshold type
  
     // Protected members
     GObservations            m_obs;        //!< Observations container
@@ -23,16 +23,16 @@ prefix, s, h, &quot;selected_&quot;,,,&quot;Prefix for event lists in observation definition fi
 #
 # Selection parameters
 #=====================
-usepnt, b, h,            no,,, "Use pointing instead of RA/DEC parameters?"
-ra,     r, a,  83.63,   0,360, "RA for ROI centre (degrees)"
-dec,    r, a,  22.01, -90, 90, "Dec for ROI centre (degrees)"
-rad,    r, a,   3.0,    0,180, "Radius of ROI (degrees)"
-tmin,   r, a,   0.0,    0,   , "Start time (CTA MET in seconds)"
-tmax,   r, a,   0.0,    0,   , "End time (CTA MET in seconds)"
-emin,   r, a,   0.1,    0,   , "Lower energy limit (TeV)"
-emax,   r, a, 100.0,    0,   , "Upper energy limit (TeV)"
-expr,   s, h,    "",     ,   , "Additional selection expression"
-usethres,s,h,"NONE",NONE|DEFAULT|USER, ,"Threshold applying algorithm"
+usepnt,  b, h,            no,,, "Use pointing instead of RA/DEC parameters?"
+ra,      r, a,  83.63,   0,360, "RA for ROI centre (degrees)"
+dec,     r, a,  22.01, -90, 90, "Dec for ROI centre (degrees)"
+rad,     r, a,   3.0,    0,180, "Radius of ROI (degrees)"
+tmin,    r, a,   0.0,    0,   , "Start time (CTA MET in seconds)"
+tmax,    r, a,   0.0,    0,   , "End time (CTA MET in seconds)"
+emin,    r, a,   0.1,    0,   , "Lower energy limit (TeV)"
+emax,    r, a, 100.0,    0,   , "Upper energy limit (TeV)"
+expr,    s, h,    "",     ,   , "Additional selection expression"
+usethres,s, h,"NONE",NONE|DEFAULT|USER,,"Energy threshold type (NONE|DEFAULT|USER)"
 #
 # Standard parameters
 #====================