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 0254ff0..2007dc6 100644 --- a/inst/cta/include/GCTAObservation.hpp +++ b/inst/cta/include/GCTAObservation.hpp @@ -120,8 +120,8 @@ public: void eventfile(const std::string& filename); const std::string& eventfile(void) const; void dispose_events(void); - double lo_user_thres(void) const; - double hi_user_thres(void) const; + const double& lo_user_thres(void) const; + const double& hi_user_thres(void) const; protected: // Protected methods @@ -133,21 +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) - - // User defined energy thresholds - double m_lo_user_thres; //!< user defined lower energy threshold - double m_hi_user_thres; //!< user defined upper energy boundary + 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 @@ -395,23 +393,23 @@ const std::string& GCTAObservation::eventfile(void) const /***********************************************************************//** - * @brief Return low energy threshold from IRF + * @brief Return user low energy threshold * - * @return return low energy threshold + * @return User low energy threshold. ***************************************************************************/ inline -double GCTAObservation::lo_user_thres(void) const +const double& GCTAObservation::lo_user_thres(void) const { return m_lo_user_thres; } /***********************************************************************//** - * @brief Return high energy threshold from IRF + * @brief Return user high energy threshold * - * @return return low energy threshold + * @return User high energy threshold. ***************************************************************************/ inline -double GCTAObservation::hi_user_thres(void) const +const double& GCTAObservation::hi_user_thres(void) const { return m_hi_user_thres; } diff --git a/inst/cta/include/GCTAResponseIrf.hpp b/inst/cta/include/GCTAResponseIrf.hpp index 8567952..a4d2bb5 100644 --- a/inst/cta/include/GCTAResponseIrf.hpp +++ b/inst/cta/include/GCTAResponseIrf.hpp @@ -127,9 +127,8 @@ public: void edisp(GCTAEdisp* edisp); const GCTABackground* background(void) const; void background(GCTABackground* background); - - double lo_save_thres(void) const; - double hi_save_thres(void) const; + const double& lo_save_thres(void) const; + const double& hi_save_thres(void) const; // Low-level response methods double aeff(const double& theta, @@ -180,10 +179,8 @@ private: GCTAEdisp* m_edisp; //!< Energy dispersion GCTABackground* m_background; //!< Energy dispersion mutable bool m_apply_edisp; //!< Apply energy dispersion - - // Energy thresholds - double m_lo_save_thres;//!< save energy threshold - double m_hi_save_thres;//!< save maximum energy + 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 @@ -408,10 +405,10 @@ void GCTAResponseIrf::background(GCTABackground* background) /***********************************************************************//** * @brief Return low energy threshold from IRF * - * @return return low energy threshold + * @return Low energy threshold from IRF. ***************************************************************************/ inline -double GCTAResponseIrf::lo_save_thres(void) const +const double& GCTAResponseIrf::lo_save_thres(void) const { return m_lo_save_thres; } @@ -419,14 +416,12 @@ double GCTAResponseIrf::lo_save_thres(void) const /***********************************************************************//** * @brief Return high energy threshold from IRF * - * @return return low energy threshold + * @return High energy threshold from IRF. ***************************************************************************/ inline -double GCTAResponseIrf::hi_save_thres(void) const +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 3311757..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 @@ -652,14 +655,11 @@ void GCTAObservation::read(const GXmlElement& xml) // Extract instrument name m_instrument = xml.attribute("instrument"); -// m_name = xml.attribute("name"); -// m_id = xml.attribute("id"); // 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")); } @@ -804,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 { @@ -824,11 +829,12 @@ 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)); + 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)); + xml.attribute("emax", gammalib::str(m_hi_user_thres)); } // If XML element has 0 nodes then add the required parameter nodes @@ -1192,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) { @@ -1240,12 +1268,12 @@ 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; @@ -1262,17 +1290,17 @@ 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; diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index ba2ecc0..1e49e79 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -949,9 +949,10 @@ void GCTAResponseIrf::load_aeff(const std::string& filename) // as CTA response table if (file.contains("EFFECTIVE AREA")) { + // Get HDU const GFitsHDU* hdu = file.at("EFFECTIVE AREA"); - // Read save energy threshold if available + // Read save energy thresholds if available if (hdu->has_card("LO_THRES")) { m_lo_save_thres = hdu->real("LO_THRES"); } @@ -959,26 +960,35 @@ void GCTAResponseIrf::load_aeff(const std::string& filename) 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 threshold if available + // Read save energy thresholds if available if (hdu->has_card("LO_THRES")) { - m_lo_save_thres = file["SPECRESP"]->real("LO_THRES"); + m_lo_save_thres = hdu->real("LO_THRES"); } if (hdu->has_card("LO_THRES")) { - m_hi_save_thres = file["SPECRESP"]->real("HI_THRES"); + m_hi_save_thres = hdu->real("HI_THRES"); } + // Close file file.close(); + + // Allocate Aeff from file m_aeff = new GCTAAeffArf(filename); + } } @@ -1243,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)); @@ -2927,13 +2959,11 @@ 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 - - // Energy thresholds + 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; @@ -2964,11 +2994,9 @@ 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; - - // Copy energy thresholds + 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;