diff --git a/ChangeLog b/ChangeLog index 293769d..06d5429 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-10-31 Juergen Knoedlseder +2014-11-07 Juergen Knoedlseder * Version 0.9.0 released ======================== @@ -34,6 +34,7 @@ Introduce class codes to avoid dynamic type casting Remove error computation from optimization Add GOptimize::errors() method for error computation + Add support for user and save CTA energy thresholds Several bug corrections diff --git a/NEWS b/NEWS index 6827ef3..5ee06b4 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-09-00 Author(s) : Juergen Knoedlseder -31 October 2014 +7 November 2014 1. Introduction @@ -307,6 +307,11 @@ Events are now only loaded when needed into GCTAObservation. This reduces memory the usage. The method dispose_events() has been added to explicitely allow the disposal of the events. +Support for user provided energy thresholds and save energy thresholds has been +implemented. Save energy thresholds will be extracted from the LO_THRES and +HI_THRES keywords if they are present in the effective area response FITS +file. The user may supply energy thresholds in the XML file. + 21. COMPTEL interface --------------------- diff --git a/inst/cta/include/GCTAObservation.hpp b/inst/cta/include/GCTAObservation.hpp index c80d9c0..2007dc6 100644 --- a/inst/cta/include/GCTAObservation.hpp +++ b/inst/cta/include/GCTAObservation.hpp @@ -120,6 +120,8 @@ public: void eventfile(const std::string& filename); const std::string& eventfile(void) const; void dispose_events(void); + const double& lo_user_thres(void) const; + const double& hi_user_thres(void) const; protected: // Protected methods @@ -131,17 +133,19 @@ protected: void set_event_type(void); // Protected members - std::string m_instrument; //!< Instrument name - std::string m_eventfile; //!< Event filename - std::string m_eventtype; //!< Event type (for XML file) - GCTAResponse* m_response; //!< Pointer to instrument response functions - GCTAPointing m_pointing; //!< Pointing direction - int m_obs_id; //!< Observation ID - double m_ontime; //!< Ontime (seconds) - double m_livetime; //!< Livetime (seconds) - double m_deadc; //!< Deadtime correction (livetime/ontime) - double m_ra_obj; //!< Right Ascension of object (degrees) - double m_dec_obj; //!< Declination of object (degrees) + std::string m_instrument; //!< Instrument name + std::string m_eventfile; //!< Event filename + std::string m_eventtype; //!< Event type (for XML file) + GCTAResponse* m_response; //!< Pointer to instrument response functions + GCTAPointing m_pointing; //!< Pointing direction + int m_obs_id; //!< Observation ID + double m_ontime; //!< Ontime (seconds) + double m_livetime; //!< Livetime (seconds) + double m_deadc; //!< Deadtime correction (livetime/ontime) + double m_ra_obj; //!< Right Ascension of object (degrees) + double m_dec_obj; //!< Declination of object (degrees) + double m_lo_user_thres; //!< User defined lower energy threshold + double m_hi_user_thres; //!< User defined upper energy boundary // Special protected member for GCTAModelCubeBackground friend std::string m_bgdfile; //!< Background filename @@ -387,4 +391,30 @@ const std::string& GCTAObservation::eventfile(void) const return m_eventfile; } + +/***********************************************************************//** + * @brief Return user low energy threshold + * + * @return User low energy threshold. + ***************************************************************************/ +inline +const double& GCTAObservation::lo_user_thres(void) const +{ + return m_lo_user_thres; +} + +/***********************************************************************//** + * @brief Return user high energy threshold + * + * @return User high energy threshold. + ***************************************************************************/ +inline +const double& GCTAObservation::hi_user_thres(void) const +{ + return m_hi_user_thres; +} + + + + #endif /* GCTAOBSERVATION_HPP */ diff --git a/inst/cta/include/GCTAResponse.hpp b/inst/cta/include/GCTAResponse.hpp index 3329420..0a65c20 100644 --- a/inst/cta/include/GCTAResponse.hpp +++ b/inst/cta/include/GCTAResponse.hpp @@ -93,6 +93,7 @@ protected: const GObservation& obs) const; const GCTAInstDir& retrieve_dir(const std::string& origin, const GEvent& event) const; + }; #endif /* GCTARESPONSE_HPP */ diff --git a/inst/cta/include/GCTAResponseIrf.hpp b/inst/cta/include/GCTAResponseIrf.hpp index 0850b96..a4d2bb5 100644 --- a/inst/cta/include/GCTAResponseIrf.hpp +++ b/inst/cta/include/GCTAResponseIrf.hpp @@ -127,6 +127,8 @@ public: void edisp(GCTAEdisp* edisp); const GCTABackground* background(void) const; void background(GCTABackground* background); + const double& lo_save_thres(void) const; + const double& hi_save_thres(void) const; // Low-level response methods double aeff(const double& theta, @@ -177,6 +179,8 @@ private: GCTAEdisp* m_edisp; //!< Energy dispersion GCTABackground* m_background; //!< Energy dispersion mutable bool m_apply_edisp; //!< Apply energy dispersion + double m_lo_save_thres; //!< Save low energy threshold + double m_hi_save_thres; //!< Save high energy threshold // XML response filename std::string m_xml_caldb; //!< Calibration database string in XML file @@ -398,4 +402,26 @@ void GCTAResponseIrf::background(GCTABackground* background) return; } +/***********************************************************************//** + * @brief Return low energy threshold from IRF + * + * @return Low energy threshold from IRF. + ***************************************************************************/ +inline +const double& GCTAResponseIrf::lo_save_thres(void) const +{ + return m_lo_save_thres; +} + +/***********************************************************************//** + * @brief Return high energy threshold from IRF + * + * @return High energy threshold from IRF. + ***************************************************************************/ +inline +const double& GCTAResponseIrf::hi_save_thres(void) const +{ + return m_hi_save_thres; +} + #endif /* GCTARESPONSEIRF_HPP */ diff --git a/inst/cta/pyext/GCTAObservation.i b/inst/cta/pyext/GCTAObservation.i index 08a3ce5..a42c6bb 100644 --- a/inst/cta/pyext/GCTAObservation.i +++ b/inst/cta/pyext/GCTAObservation.i @@ -93,6 +93,8 @@ public: void eventfile(const std::string& filename); const std::string& eventfile(void) const; void dispose_events(void); + const double& lo_user_thres(void) const; + const double& hi_user_thres(void) const; }; diff --git a/inst/cta/src/GCTAObservation.cpp b/inst/cta/src/GCTAObservation.cpp index 3c8ce81..71a983f 100644 --- a/inst/cta/src/GCTAObservation.cpp +++ b/inst/cta/src/GCTAObservation.cpp @@ -638,7 +638,10 @@ GEbounds GCTAObservation::ebounds(void) const * ... * * - * for a binned observation. + * for a binned observation. Optional user energy thresholds can be + * specified by adding the emin and emax attributes: + * + * * * The method does no load the events into memory but stores the file name * of the event file. The events are only loaded when required. This reduces @@ -653,6 +656,14 @@ void GCTAObservation::read(const GXmlElement& xml) // Extract instrument name m_instrument = xml.attribute("instrument"); + // Read in user defined energy boundaries of this observation + if (xml.attribute("emin") != "") { + m_lo_user_thres = gammalib::todouble(xml.attribute("emin")); + } + if (xml.attribute("emax") != "") { + m_hi_user_thres = gammalib::todouble(xml.attribute("emax")); + } + // Determine number of parameter nodes in XML element int npars = xml.elements("parameter"); @@ -793,7 +804,12 @@ void GCTAObservation::read(const GXmlElement& xml) * ... * * - * for a binned observation. + * for a binned observation. If user energy threshold are defined (i.e. + * threshold values are >0) the additional emin and emax attributes will + * be written to the file: + * + * + * ***************************************************************************/ void GCTAObservation::write(GXmlElement& xml) const { @@ -813,6 +829,14 @@ void GCTAObservation::write(GXmlElement& xml) const throw GException::invalid_value(G_WRITE, msg); } + // Add user energy threshold attributes (if set) + if (m_lo_user_thres > 0.0) { + xml.attribute("emin", gammalib::str(m_lo_user_thres)); + } + if (m_hi_user_thres > 0.0) { + xml.attribute("emax", gammalib::str(m_hi_user_thres)); + } + // If XML element has 0 nodes then add the required parameter nodes if (xml.elements() == 0) { std::string text = "parameter name=\""+m_eventtype+"\""; @@ -1033,6 +1057,7 @@ void GCTAObservation::load(const std::string& cntcube, ***************************************************************************/ void GCTAObservation::save(const std::string& filename, const bool& clobber) const { + // Create FITS file GFits fits; @@ -1173,6 +1198,28 @@ std::string GCTAObservation::print(const GChatter& chatter) const result.append("\n"+gammalib::parformat("Deadtime correction")); result.append(gammalib::str(m_deadc)); + // Append user energy threshold information + result.append("\n"+gammalib::parformat("User energy range")); + if (m_lo_user_thres > 0.0 && m_hi_user_thres) { + result.append(gammalib::str(m_lo_user_thres)); + result.append(" - "); + result.append(gammalib::str(m_hi_user_thres)); + result.append(" TeV"); + } + else if (m_lo_user_thres > 0.0) { + result.append("> "); + result.append(gammalib::str(m_lo_user_thres)); + result.append(" TeV"); + } + else if (m_hi_user_thres > 0.0) { + result.append("< "); + result.append(gammalib::str(m_hi_user_thres)); + result.append(" TeV"); + } + else { + result.append("undefined"); + } + // Append detailed information GChatter reduced_chatter = gammalib::reduce(chatter); if (reduced_chatter > SILENT) { @@ -1221,12 +1268,14 @@ void GCTAObservation::init_members(void) m_bgdfile.clear(); m_response = NULL; m_pointing.clear(); - m_obs_id = 0; - m_ontime = 0.0; - m_livetime = 0.0; - m_deadc = 0.0; - m_ra_obj = 0.0; - m_dec_obj = 0.0; + m_obs_id = 0; + m_ontime = 0.0; + m_livetime = 0.0; + m_deadc = 0.0; + m_ra_obj = 0.0; + m_dec_obj = 0.0; + m_lo_user_thres = 0.0; + m_hi_user_thres = 0.0; // Return return; @@ -1241,17 +1290,19 @@ void GCTAObservation::init_members(void) void GCTAObservation::copy_members(const GCTAObservation& obs) { // Copy members - m_instrument = obs.m_instrument; - m_eventfile = obs.m_eventfile; - m_eventtype = obs.m_eventtype; - m_bgdfile = obs.m_bgdfile; - m_pointing = obs.m_pointing; - m_obs_id = obs.m_obs_id; - m_ontime = obs.m_ontime; - m_livetime = obs.m_livetime; - m_deadc = obs.m_deadc; - m_ra_obj = obs.m_ra_obj; - m_dec_obj = obs.m_dec_obj; + m_instrument = obs.m_instrument; + m_eventfile = obs.m_eventfile; + m_eventtype = obs.m_eventtype; + m_bgdfile = obs.m_bgdfile; + m_pointing = obs.m_pointing; + m_obs_id = obs.m_obs_id; + m_ontime = obs.m_ontime; + m_livetime = obs.m_livetime; + m_deadc = obs.m_deadc; + m_ra_obj = obs.m_ra_obj; + m_dec_obj = obs.m_dec_obj; + m_lo_user_thres = obs.m_lo_user_thres; + m_hi_user_thres = obs.m_hi_user_thres; // Clone members m_response = (obs.m_response != NULL) ? obs.m_response->clone() : NULL; diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index 048a8a1..1e49e79 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -948,15 +948,47 @@ void GCTAResponseIrf::load_aeff(const std::string& filename) // If file contains an "EFFECTIVE AREA" extension then load it // as CTA response table if (file.contains("EFFECTIVE AREA")) { + + // Get HDU + const GFitsHDU* hdu = file.at("EFFECTIVE AREA"); + + // Read save energy thresholds if available + if (hdu->has_card("LO_THRES")) { + m_lo_save_thres = hdu->real("LO_THRES"); + } + if (hdu->has_card("HI_THRES")) { + m_hi_save_thres = hdu->real("HI_THRES"); + } + + // Close file file.close(); + + // Allocate Aeff from file m_aeff = new GCTAAeff2D(filename); + } // ... else if file contains a "SPECRESP" extension then load it // as ARF else if (file.contains("SPECRESP")) { + + // Get HDU + const GFitsHDU* hdu = file.at("SPECRESP"); + + // Read save energy thresholds if available + if (hdu->has_card("LO_THRES")) { + m_lo_save_thres = hdu->real("LO_THRES"); + } + if (hdu->has_card("LO_THRES")) { + m_hi_save_thres = hdu->real("HI_THRES"); + } + + // Close file file.close(); + + // Allocate Aeff from file m_aeff = new GCTAAeffArf(filename); + } } @@ -1221,6 +1253,28 @@ std::string GCTAResponseIrf::print(const GChatter& chatter) const } } + // Append save energy threshold information + result.append("\n"+gammalib::parformat("Save energy range")); + if (m_lo_save_thres > 0.0 && m_hi_save_thres) { + result.append(gammalib::str(m_lo_save_thres)); + result.append(" - "); + result.append(gammalib::str(m_hi_save_thres)); + result.append(" TeV"); + } + else if (m_lo_save_thres > 0.0) { + result.append("> "); + result.append(gammalib::str(m_lo_save_thres)); + result.append(" TeV"); + } + else if (m_hi_save_thres > 0.0) { + result.append("< "); + result.append(gammalib::str(m_hi_save_thres)); + result.append(" TeV"); + } + else { + result.append("undefined"); + } + // Append calibration database result.append("\n"+m_caldb.print(chatter)); @@ -2905,11 +2959,13 @@ void GCTAResponseIrf::init_members(void) // Initialise members m_caldb.clear(); m_rspname.clear(); - m_aeff = NULL; - m_psf = NULL; - m_edisp = NULL; - m_background = NULL; - m_apply_edisp = false; //!< Switched off by default + m_aeff = NULL; + m_psf = NULL; + m_edisp = NULL; + m_background = NULL; + m_apply_edisp = false; //!< Switched off by default + m_lo_save_thres = 0.0; + m_hi_save_thres = 0.0; // XML response filenames m_xml_caldb.clear(); @@ -2938,9 +2994,11 @@ void GCTAResponseIrf::init_members(void) void GCTAResponseIrf::copy_members(const GCTAResponseIrf& rsp) { // Copy members - m_caldb = rsp.m_caldb; - m_rspname = rsp.m_rspname; - m_apply_edisp = rsp.m_apply_edisp; + m_caldb = rsp.m_caldb; + m_rspname = rsp.m_rspname; + m_apply_edisp = rsp.m_apply_edisp; + m_lo_save_thres = rsp.m_lo_save_thres; + m_hi_save_thres = rsp.m_hi_save_thres; // Copy response filenames m_xml_caldb = rsp.m_xml_caldb; diff --git a/src/obs/GObservations.cpp b/src/obs/GObservations.cpp index a8b15fe..553ae62 100644 --- a/src/obs/GObservations.cpp +++ b/src/obs/GObservations.cpp @@ -646,7 +646,7 @@ void GObservations::write(GXml& xml) const // Write all observations into library for (int i = 0; i < size(); ++i) { - + // Initialise pointer on observation GXmlElement* obs = NULL; @@ -670,7 +670,7 @@ void GObservations::write(GXml& xml) const obs->attribute("id", m_obs[i]->id()); obs->attribute("instrument", m_obs[i]->instrument()); } - + // Write now observation m_obs[i]->write(*obs);