Added functionality to read energy boundaries from IRF Fits file as well from ob…

…servation definition XML files

Added functionality to read energy boundaries from IRF Fits file as well from ob…
…servation definition XML files
Michael Mayer
1 parent 4dc6a2e5
Showing 6 changed files with 119 additions and 2 deletions
inst/cta/include/GCTAObservation.hpp
inst/cta/include/GCTAResponse.hpp
inst/cta/include/GCTAResponseIrf.hpp
inst/cta/src/GCTAObservation.cpp
inst/cta/src/GCTAResponseIrf.cpp
src/obs/GObservations.cpp
@@ -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&amp; 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 */
@@ -93,6 +93,7 @@ protected:
                                         const GObservation& obs) const;
     const GCTAInstDir&     retrieve_dir(const std::string& origin,
                                         const GEvent&      event) const;
+
 };
 #endif /* GCTARESPONSE_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 */
@@ -652,6 +652,17 @@ void GCTAObservation::read(const GXmlElement&amp; 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&amp; 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&amp; 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&amp; 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;
@@ -948,6 +948,17 @@ void GCTAResponseIrf::load_aeff(const std::string&amp; 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&amp; 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&amp; 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;
@@ -646,7 +646,7 @@ void GObservations::write(GXml&amp; 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&amp; xml) const
             obs->attribute("id", m_obs[i]->id());
             obs->attribute("instrument", m_obs[i]->instrument());
         }
-    
+
         // Write now observation
         m_obs[i]->write(*obs);