Merge branch '1996-ctselect-add-phase-selection' into integration

# Conflicts: # src/ctselect/ctselect.par

Merge branch '1996-ctselect-add-phase-selection' into integration
# Conflicts: # src/ctselect/ctselect.par
Jürgen Knödlseder
2 parents b1bd4caa 2e1aefc0
Showing 8 changed files with 322 additions and 27 deletions
ChangeLog
NEWS
cscripts/Makefile.am
doc/source/users/reference_manual/ctselect.rst
src/ctselect/ctselect.cpp
src/ctselect/ctselect.hpp
src/ctselect/ctselect.par
test/test_ctselect.py
-2017-04-11
+2017-04-13
  
         * Version 1.3.0 released
         ========================
  
+        Add "phase" parameter to ctselect (#1996)
         Add ctphase tool (#1982)
         Add ctprob tool (#1991)
         ctlike stores covariance matrix in either CSV or FITS format (#1863)
 New Features and Important Changes in ctools 1.3.0.dev1
  
-11 April 2017
+13 April 2017
  
  
 Introduction
@@ -47,7 +47,12 @@ Now, for map cubes outside the RoI the flux is set to zero.
  
 ctselect - CTA event selection
 ------------------------------
-None
+A "phase" parameter has been added to ctselect to allow event selection
+according to their phase information. The phase parameter accepts a
+string that provides a comma separated list of phase intervals in the
+format "0.1:0.4,0.6:0.8". Wrap around is supported (e.g. "0.8:0.1"
+selects events with phases in the intervals [0.8,1.0] and [0.0,0.1]
+(#1996).
  
  
 ctbin - CTA event binning
@@ -127,6 +127,7 @@ testscripts = $(srcdir)/tests/__init__.py \
               $(top_srcdir)/test/test_csobsdef.py \
               $(top_srcdir)/test/test_csobsinfo.py \
               $(top_srcdir)/test/test_csobsselect.py \
+              $(top_srcdir)/test/test_csphasecrv.py \
               $(top_srcdir)/test/test_cspull.py \
               $(top_srcdir)/test/test_csresmap.py \
               $(top_srcdir)/test/test_csroot2caldb.py \
@@ -52,39 +52,48 @@ General parameters
  
 ``ra [real]``
     Right Ascension of acceptance cone (or ROI) centre (J2000, in degrees).
-    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, not ROI
+    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, no ROI
     selection will be performed.
  
 ``dec [real]``
     Declination of acceptance cone (or ROI) centre (J2000, in degrees).
-    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, not ROI
+    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, no ROI
     selection will be performed.
  
 ``rad [real]``
     Radius of acceptance cone (or ROI) centre (in degrees).
-    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, not ROI
+    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, no ROI
     selection will be performed.
  
 ``tmin [real]``
     CTA mission elapsed start time defining the earliest event (in seconds).
-    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, not time
+    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, no time
     selection will be performed.
  
 ``tmax [real]``
     CTA mission elapsed stop time defining the latest event (in seconds).
-    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, not time
+    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, no time
     selection will be performed.
  
 ``emin [real]``
     Lower energy limit of events (in TeV).
-    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, not energy
+    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, no energy
     selection will be performed.
  
 ``emax [real]``
     Upper energy limit of events (in TeV).
-    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, not energy
+    If INDEF, NONE, UNDEF or UNDEFINED is passed as value, no energy
     selection will be performed.
  
+``(phase = NONE) [string]``
+    String to apply a phase selection. The string must contain the boundaries 
+    of the phase interval to be selected separated be a colon. More than one
+    interval can be specified at the same time. In this case intervals must be 
+    separated by a comma. Examples of valid strings are: phase = "0.3:0.6",
+    phase = "0.3:0.6,0.8:0.9" or phase = "0.8:0.2". In the last case, events
+    with phases in the intervals [0.8,1.0] and [0.0,0.2] are selected. If NONE
+    is passed as value, no phase selection will be performed.
+ 	 	 
 ``(expr = "") [string]``
     Additional event selection expression (cfitsio syntax).
  
@@ -34,6 +34,7 @@
  
 /* __ Method name definitions ____________________________________________ */
 #define G_RUN                                               "ctselect::run()"
+#define G_GET_PARAMETERS                         "ctselect::get_parameters()"
 #define G_SELECT_EVENTS          "ctselect::select_events(GCTAObservation*, "\
                                   "std::string&, std::string&, std::string&)"
 #define G_SET_EBOUNDS    "ctselect::set_ebounds(GCTAObservation*, GEbounds&)"
@@ -481,9 +482,12 @@ void ctselect::init_members(void)
     m_gtiname.clear();
     m_timemin.clear();
     m_timemax.clear();
-    m_select_energy = true;
-    m_select_roi    = true;
-    m_select_time   = true;
+    m_phasemin.clear();
+    m_phasemax.clear();
+    m_select_energy = false;
+    m_select_roi    = false;
+    m_select_time   = false;
+    m_select_phase  = false;
  
     // Return
     return;
@@ -516,9 +520,12 @@ void ctselect::copy_members(const ctselect&amp; app)
     m_gtiname       = app.m_gtiname;
     m_timemin       = app.m_timemin;
     m_timemax       = app.m_timemax;
+    m_phasemin      = app.m_phasemin;
+    m_phasemax      = app.m_phasemax;
     m_select_energy = app.m_select_energy;
     m_select_roi    = app.m_select_roi;
     m_select_time   = app.m_select_time;
+    m_select_phase  = app.m_select_phase;
  
     // Return
     return;
@@ -548,9 +555,10 @@ void ctselect::free_members(void)
 void ctselect::get_parameters(void)
 {
     // Initialise selection flags
-    m_select_energy = true;
-    m_select_roi    = true;
-    m_select_time   = true;
+    m_select_energy = false;
+    m_select_roi    = false;
+    m_select_time   = false;
+    m_select_phase  = false;
  
     // Setup observations from "inobs" parameter. Do not request response
     // information and do not accept counts cubes.
@@ -563,7 +571,7 @@ void ctselect::get_parameters(void)
     m_roi        = get_roi();
     m_select_roi = m_roi.is_valid();
  
-    // Check for sanity of time selection parameters
+    // Get time selection parameters
     if ((*this)["tmin"].is_valid() && (*this)["tmax"].is_valid()) {
  
         // Get User parameters
@@ -578,19 +586,109 @@ void ctselect::get_parameters(void)
             m_select_time = true;
         }
     }
-    else {
-        m_select_time = false;
-    }
  
-    // Check for sanity of energy selection parameters
+    // Get energy selection parameters
     if ((*this)["emin"].is_valid() && (*this)["emax"].is_valid()) {
         m_emin          = (*this)["emin"].real();
         m_emax          = (*this)["emax"].real();
         m_select_energy = true;
     }
-    else {
-        m_select_energy = false;
-    }
+
+    // Get phase selection string
+    std::string phase_expr = gammalib::strip_whitespace(gammalib::tolower(
+                                       (*this)["phase"].string()));
+
+    // Get phase selection parameters
+    if ((phase_expr.length() > 0) && (phase_expr != "none")) {
+
+        // Get phase selection string
+        std::string phase_expr = (*this)["phase"].string();
+
+        // Extract phase bins from selection string
+        std::vector<std::string> phase_splits = gammalib::split(phase_expr, ",");
+        for (int i = 0; i < phase_splits.size(); ++i) {
+
+            // Get minimum and maximum of phase bin
+            std::vector<std::string> phase_toks =
+                                     gammalib::split(phase_splits[i], ":");
+
+            // Throw an exception if there are not two elements
+            if (phase_toks.size() != 2) {
+                std::string msg = "Invalid phase selection string \""+
+                                  phase_splits[i]+"\" encountered. The phase "
+                                  "selection string requires a minimum and "
+                                  "maximum value separated by a colon, e.g. "
+                                  "\"0.3:0.5\".";
+                throw GException::invalid_value(G_GET_PARAMETERS, msg);
+            }
+
+            // Throw an exception if one of the elements is empty
+            if (gammalib::strip_whitespace(phase_toks[0]).length() == 0) {
+                std::string msg = "Invalid phase selection string \""+
+                                  phase_splits[i]+"\" encountered. The phase "
+                                  "selection string requires a minimum and "
+                                  "maximum value separated by a colon, e.g. "
+                                  "\"0.3:0.5\".";
+                throw GException::invalid_value(G_GET_PARAMETERS, msg);
+            }
+            if (gammalib::strip_whitespace(phase_toks[1]).length() == 0) {
+                std::string msg = "Invalid phase selection string \""+
+                                  phase_splits[i]+"\" encountered. The phase "
+                                  "selection string requires a minimum and "
+                                  "maximum value separated by a colon, e.g. "
+                                  "\"0.3:0.5\".";
+                throw GException::invalid_value(G_GET_PARAMETERS, msg);
+            }
+
+            // Get phase bin boundaries
+            double phase_min = gammalib::todouble(phase_toks[0]);
+            double phase_max = gammalib::todouble(phase_toks[1]);
+
+            // Throw an exception if phase minimum or maximum is outside the
+            // range [0,1]
+            if (phase_min < 0.0 || phase_min > 1.0) {
+                std::string msg = "Phase minimum "+gammalib::str(phase_min)+
+                                  " outside the valid range [0,1]. Please "
+                                  "specify phase interval boundaries comprised "
+                                  "within 0 and 1.";
+                throw GException::invalid_value(G_GET_PARAMETERS, msg);
+            }
+            if (phase_max < 0.0 || phase_max > 1.0) {
+                std::string msg = "Phase maximum "+gammalib::str(phase_max)+
+                                  " outside the valid range [0,1]. Please "
+                                  "specify phase interval boundaries comprised "
+                                  "within 0 and 1.";
+                throw GException::invalid_value(G_GET_PARAMETERS, msg);
+            }
+
+            // If phase minimum is smaller than phase maximum then append a
+            // single interval
+            if (phase_min < phase_max) {
+                m_phasemin.push_back(phase_min);
+                m_phasemax.push_back(phase_max);
+            }
+
+            // Otherwise if the phase minimum is larger than the maximum then
+            // consider this a wrap around interval and append one interval
+            // from phase minimum to 1 and another interval from 0 to phase
+            // maximum
+            else if (phase_min > phase_max) {
+                if (phase_min < 1.0) {
+                    m_phasemin.push_back(phase_min);
+                    m_phasemax.push_back(1.0);
+                }
+                if (phase_max > 0.0) {
+                    m_phasemin.push_back(0.0);
+                    m_phasemax.push_back(phase_max);
+                }
+            }
+
+        } // endfor: looped over phase selection string
+
+        // Signal that a phase selection information is available
+        m_select_phase = (m_phasemin.size() > 0);
+
+    } // endif: phase selection parameters were valid
  
     // Get other User parameters
     m_expr     = (*this)["expr"].string();
@@ -722,6 +820,54 @@ void ctselect::select_events(GCTAObservation*   obs,
  
     } // endif: made time selection
  
+    // Make phase selection
+    if (m_select_phase) {
+
+        // Check if event_list has phase
+        if (list->has_phase()) { 
+
+            // Loop over phase selections
+            for (int i = 0; i < m_phasemin.size(); ++i) {
+
+                // Check if phasemax is larger than phasemin
+                if (i == 0) {
+                    selection += add + "( ";
+                }
+                else {
+                    selection += add;
+                }
+
+                // Format phase with sufficient accuracy and add to selection string
+                char cmin[80];
+                char cmax[80];
+                sprintf(cmin, "%.8f", m_phasemin[i]);
+                sprintf(cmax, "%.8f", m_phasemax[i]);
+                selection += "(PHASE >= "+std::string(cmin)+" && PHASE <= "+
+                             std::string(cmax)+")";
+                add        = " || ";
+                log_value(NORMAL, "Phase range "+gammalib::str(i+1),
+                                  gammalib::str(m_phasemin[i])+" - "+
+                                  gammalib::str(m_phasemax[i]));
+
+            } // end for loop
+
+            // If phase intervals have been appended then close the selection
+            // string and take provision for appending another selection
+            if (m_phasemin.size() > 0) {
+                selection += " )";
+                add        = " && ";
+            }
+
+        } // endif: there was PHASE information in the event list
+
+        // Otherwise signal that no PHASE information is available
+        else {
+            log_value(NORMAL, "Phase range", "Event list has no \"PHASE\" "
+                              "column. Phase selection skipped.");
+        }
+
+    } // endif: made phase selection
+
     // Make energy selection
     if (m_select_energy) {
  
@@ -99,10 +99,12 @@ protected:
     std::vector<std::string> m_gtiname;       //!< GTI extension names
     GTime                    m_timemin;       //!< Earliest time
     GTime                    m_timemax;       //!< Latest time
+    std::vector<double>      m_phasemin;      //!< Phase bin minima
+    std::vector<double>      m_phasemax;      //!< Phase bin maxima
     bool                     m_select_energy; //!< Perform energy selection
     bool                     m_select_roi;    //!< Perform RoI selection
     bool                     m_select_time;   //!< Perform time selection
+    bool                     m_select_phase;  //!< Perform phase selection
 };
  
-
 #endif /* CTSELECT_HPP */
@@ -31,6 +31,7 @@ tmin,    r, a, 0.0,,, &quot;Start time (CTA MET in seconds)&quot;
 tmax,    r, a, 0.0,,, "End time (CTA MET in seconds)"
 emin,    r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,    r, a, 100.0,,, "Upper energy limit (TeV)"
+phase,   s, h, NONE,,, "Phase expression in the format phasemin0:phasemax0,phasemin1:phasemax1,..."
 expr,    s, h, ,,, "Additional selection expression"
 usethres,s, h, NONE,NONE|DEFAULT|USER,, "Energy threshold type"
  
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the ctselect tool.
 #
-# Copyright (C) 2014-2016 Juergen Knoedlseder
+# Copyright (C) 2014-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -40,6 +40,7 @@ class Test(test):
         test.__init__(self)
  
         # Set test data
+        self._phased_events  = self._datadir + '/phased_events.fits'
         self._invalid_events = self._datadir + '/invalid_event_list.fits'
  
         # Return
@@ -355,7 +356,136 @@ class Test(test):
         select.execute()
  
         # Check result file
-        self._check_result_file('ctselect_py10.fits', nevents=0, dec=-22.01, rad=5.0)
+        self._check_result_file('ctselect_py10.fits', nevents=0, dec=-22.01,
+                                rad=5.0)
+
+        # Now test phase cut
+        select = ctools.ctselect()
+        select['inobs']   = self._phased_events
+        select['ra']      = 83.63
+        select['dec']     = 22.01
+        select['rad']     = 5.0
+        select['tmin']    = 'INDEF'
+        select['tmax']    = 'INDEF'
+        select['emin']    = 0.1
+        select['emax']    = 100.0
+        select['phase']   = '0.1:0.4,0.6:0.8'
+        select['outobs']  = 'ctselect_py11.fits'
+        select['logfile'] = 'ctselect_py11.log'
+        select['chatter'] = 2
+
+        # Run ctselect tool
+        select.logFileOpen()   # Make sure we get a log file
+        select.execute()
+
+        # Check result file
+        self._check_result_file('ctselect_py11.fits', nevents=1909, rad=5.0)
+
+        # Now test anoother phase cut
+        select = ctools.ctselect()
+        select['inobs']   = self._phased_events
+        select['ra']      = 83.63
+        select['dec']     = 22.01
+        select['rad']     = 5.0
+        select['tmin']    = 'INDEF'
+        select['tmax']    = 'INDEF'
+        select['emin']    = 0.1
+        select['emax']    = 100.0
+        select['phase']   = '0.8:0.3'
+        select['outobs']  = 'ctselect_py12.fits'
+        select['logfile'] = 'ctselect_py12.log'
+        select['chatter'] = 2
+
+        # Run ctselect tool
+        select.logFileOpen()   # Make sure we get a log file
+        select.execute()
+
+        # Check result file
+        self._check_result_file('ctselect_py12.fits', nevents=1995, rad=5.0)
+
+        # Test invalid minimum phase value
+        self.test_try('Test invalid minimum phase value')
+        try:
+            select['phase']   = '-0.1:0.3'
+            select['outobs']  = 'ctselect_py13.fits'
+            select['logfile'] = 'ctselect_py13.log'
+            select.logFileOpen()
+            select.execute()
+            self.test_try_failure('Exception not thrown for "-0.1:0.3"')
+        except ValueError:
+            self.test_try_success()
+
+        # Test invalid minimum phase value
+        self.test_try('Test invalid minimum phase value')
+        try:
+            select['phase']   = '1.1:0.3'
+            select['outobs']  = 'ctselect_py14.fits'
+            select['logfile'] = 'ctselect_py14.log'
+            select.logFileOpen()
+            select.execute()
+            self.test_try_failure('Exception not thrown for "1.1:0.3"')
+        except ValueError:
+            self.test_try_success()
+
+        # Test invalid maximum phase value
+        self.test_try('Test invalid maximum phase value')
+        try:
+            select['phase']   = '0.1:-0.3'
+            select['outobs']  = 'ctselect_py15.fits'
+            select['logfile'] = 'ctselect_py15.log'
+            select.logFileOpen()
+            select.execute()
+            self.test_try_failure('Exception not thrown for "0.1:-0.3"')
+        except ValueError:
+            self.test_try_success()
+
+        # Test invalid maximum phase value
+        self.test_try('Test invalid maximum phase value')
+        try:
+            select['phase']   = '0.1:1.1'
+            select['outobs']  = 'ctselect_py16.fits'
+            select['logfile'] = 'ctselect_py16.log'
+            select.logFileOpen()
+            select.execute()
+            self.test_try_failure('Exception not thrown for "0.1:1.1"')
+        except ValueError:
+            self.test_try_success()
+
+        # Test invalid phase string
+        self.test_try('Test invalid phase string')
+        try:
+            select['phase']   = '0.1:'
+            select['outobs']  = 'ctselect_py17.fits'
+            select['logfile'] = 'ctselect_py17.log'
+            select.logFileOpen()
+            select.execute()
+            self.test_try_failure('Exception not thrown for "0.1:"')
+        except ValueError:
+            self.test_try_success()
+
+        # Test invalid phase string
+        self.test_try('Test invalid phase string')
+        try:
+            select['phase']   = ':1.0'
+            select['outobs']  = 'ctselect_py18.fits'
+            select['logfile'] = 'ctselect_py18.log'
+            select.logFileOpen()
+            select.execute()
+            self.test_try_failure('Exception not thrown for ":1.0"')
+        except ValueError:
+            self.test_try_success()
+
+        # Test invalid phase string
+        self.test_try('Test invalid phase string')
+        try:
+            select['phase']   = '0.1-1.0'
+            select['outobs']  = 'ctselect_py19.fits'
+            select['logfile'] = 'ctselect_py19.log'
+            select.logFileOpen()
+            select.execute()
+            self.test_try_failure('Exception not thrown for "0.1-1.0"')
+        except ValueError:
+            self.test_try_success()
  
         # Return
         return