From 272721f91ca423df0bcec0fd0bc8808fb5e9f76a Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sat, 8 Nov 2014 01:17:19 +0100 Subject: [PATCH] Integrate code. --- ChangeLog | 3 ++- NEWS | 26 +++++++++++++++++++++++++- src/ctselect/ctselect.cpp | 372 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------------------------- src/ctselect/ctselect.hpp | 5 +++-- src/ctselect/ctselect.par | 20 ++++++++++---------- 5 files changed, 254 insertions(+), 172 deletions(-) diff --git a/ChangeLog b/ChangeLog index 93248c7..1ba9210 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-10-30 Juergen Knoedlseder +2014-11-08 Juergen Knoedlseder * 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 diff --git a/NEWS b/NEWS index 71ace61..75df195 100644 --- a/NEWS +++ b/NEWS @@ -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 ---------------------------- diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 3daccb0..4dd0400 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -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(const_cast(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& filename) ra = pnt.dir().ra_deg(); dec = pnt.dir().dec_deg(); } - - // Get CTA event list pointer - GCTAEventList* list = - static_cast(const_cast(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& 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& 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(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(const_cast(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; - -} diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index ce558c4..8a85bf5 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -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;//!