diff --git a/ChangeLog b/ChangeLog index 504da29..b66b010 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -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) diff --git a/NEWS b/NEWS index 4313ee1..792ae2e 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ 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 diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 3a82f60..2adc17d 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -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 \ diff --git a/doc/source/users/reference_manual/ctselect.rst b/doc/source/users/reference_manual/ctselect.rst index b9b4a9f..f9d94ae 100644 --- a/doc/source/users/reference_manual/ctselect.rst +++ b/doc/source/users/reference_manual/ctselect.rst @@ -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). diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 42e0b37..1497d07 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -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& 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 phase_splits = gammalib::split(phase_expr, ","); + for (int i = 0; i < phase_splits.size(); ++i) { + + // Get minimum and maximum of phase bin + std::vector 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) { diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index 250b3c9..1d5abd2 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -99,10 +99,12 @@ protected: std::vector m_gtiname; //!< GTI extension names GTime m_timemin; //!< Earliest time GTime m_timemax; //!< Latest time + std::vector m_phasemin; //!< Phase bin minima + std::vector 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 */ diff --git a/src/ctselect/ctselect.par b/src/ctselect/ctselect.par index 6240c80..1b13c02 100644 --- a/src/ctselect/ctselect.par +++ b/src/ctselect/ctselect.par @@ -31,6 +31,7 @@ tmin, r, a, 0.0,,, "Start time (CTA MET in seconds)" 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" diff --git a/test/test_ctselect.py b/test/test_ctselect.py index 7b135c3..e8f9bf6 100755 --- a/test/test_ctselect.py +++ b/test/test_ctselect.py @@ -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