From c7a508a50fb0ccd40bb19837f1f285edd0cbbeed Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Thu, 6 Nov 2014 16:46:51 +0100 Subject: [PATCH] Added functionality to read energy boundaries from IRF Fits file as well from observation definition XML files --- inst/cta/include/GCTAObservation.hpp | 32 ++++++++++++++++++++++++++++++++ inst/cta/include/GCTAResponse.hpp | 1 + inst/cta/include/GCTAResponseIrf.hpp | 31 +++++++++++++++++++++++++++++++ inst/cta/src/GCTAObservation.cpp | 23 +++++++++++++++++++++++ inst/cta/src/GCTAResponseIrf.cpp | 30 ++++++++++++++++++++++++++++++ src/obs/GObservations.cpp | 4 ++-- 6 files changed, 119 insertions(+), 2 deletions(-) diff --git a/inst/cta/include/GCTAObservation.hpp b/inst/cta/include/GCTAObservation.hpp index c80d9c0..0254ff0 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); + double lo_user_thres(void) const; + double hi_user_thres(void) const; protected: // Protected methods @@ -143,6 +145,10 @@ protected: 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 + // Special protected member for GCTAModelCubeBackground friend std::string m_bgdfile; //!< Background filename }; @@ -387,4 +393,30 @@ const std::string& GCTAObservation::eventfile(void) const return m_eventfile; } + +/***********************************************************************//** + * @brief Return low energy threshold from IRF + * + * @return return low energy threshold + ***************************************************************************/ +inline +double GCTAObservation::lo_user_thres(void) const +{ + return m_lo_user_thres; +} + +/***********************************************************************//** + * @brief Return high energy threshold from IRF + * + * @return return low energy threshold + ***************************************************************************/ +inline +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..8567952 100644 --- a/inst/cta/include/GCTAResponseIrf.hpp +++ b/inst/cta/include/GCTAResponseIrf.hpp @@ -128,6 +128,9 @@ public: const GCTABackground* background(void) const; void background(GCTABackground* background); + double lo_save_thres(void) const; + double hi_save_thres(void) const; + // Low-level response methods double aeff(const double& theta, const double& phi, @@ -178,6 +181,10 @@ private: 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 + // XML response filename std::string m_xml_caldb; //!< Calibration database string in XML file std::string m_xml_rspname; //!< Response name in XML file @@ -398,4 +405,28 @@ void GCTAResponseIrf::background(GCTABackground* background) return; } +/***********************************************************************//** + * @brief Return low energy threshold from IRF + * + * @return return low energy threshold + ***************************************************************************/ +inline +double GCTAResponseIrf::lo_save_thres(void) const +{ + return m_lo_save_thres; +} + +/***********************************************************************//** + * @brief Return high energy threshold from IRF + * + * @return return low energy threshold + ***************************************************************************/ +inline +double GCTAResponseIrf::hi_save_thres(void) const +{ + return m_hi_save_thres; +} + + + #endif /* GCTARESPONSEIRF_HPP */ diff --git a/inst/cta/src/GCTAObservation.cpp b/inst/cta/src/GCTAObservation.cpp index 3c8ce81..3311757 100644 --- a/inst/cta/src/GCTAObservation.cpp +++ b/inst/cta/src/GCTAObservation.cpp @@ -652,6 +652,17 @@ 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")); + } // Determine number of parameter nodes in XML element int npars = xml.elements("parameter"); @@ -813,6 +824,13 @@ void GCTAObservation::write(GXmlElement& xml) const throw GException::invalid_value(G_WRITE, msg); } + 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 +1051,7 @@ void GCTAObservation::load(const std::string& cntcube, ***************************************************************************/ void GCTAObservation::save(const std::string& filename, const bool& clobber) const { + // Create FITS file GFits fits; @@ -1227,6 +1246,8 @@ void GCTAObservation::init_members(void) 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; @@ -1252,6 +1273,8 @@ void GCTAObservation::copy_members(const GCTAObservation& obs) 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..ba2ecc0 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -948,6 +948,17 @@ 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")) { + + const GFitsHDU* hdu = file.at("EFFECTIVE AREA"); + + // Read save energy threshold 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"); + } + file.close(); m_aeff = new GCTAAeff2D(filename); } @@ -955,6 +966,17 @@ void GCTAResponseIrf::load_aeff(const std::string& filename) // ... else if file contains a "SPECRESP" extension then load it // as ARF else if (file.contains("SPECRESP")) { + + const GFitsHDU* hdu = file.at("SPECRESP"); + + // Read save energy threshold if available + if (hdu->has_card("LO_THRES")) { + m_lo_save_thres = file["SPECRESP"]->real("LO_THRES"); + } + if (hdu->has_card("LO_THRES")) { + m_hi_save_thres = file["SPECRESP"]->real("HI_THRES"); + } + file.close(); m_aeff = new GCTAAeffArf(filename); } @@ -2911,6 +2933,10 @@ void GCTAResponseIrf::init_members(void) m_background = NULL; m_apply_edisp = false; //!< Switched off by default + // Energy thresholds + m_lo_save_thres = 0.0; + m_hi_save_thres = 0.0; + // XML response filenames m_xml_caldb.clear(); m_xml_rspname.clear(); @@ -2942,6 +2968,10 @@ void GCTAResponseIrf::copy_members(const GCTAResponseIrf& rsp) m_rspname = rsp.m_rspname; m_apply_edisp = rsp.m_apply_edisp; + // Copy energy thresholds + 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; m_xml_rspname = rsp.m_xml_rspname; 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); -- libgit2 0.21.4