Merge branch '1039-allow-energy-thresholds-via-IRF-and-XML' into integration

Jürgen Knödlseder
2 parents 76e21e25 6f4762d7
Showing 9 changed files with 216 additions and 42 deletions
ChangeLog
NEWS
inst/cta/include/GCTAObservation.hpp
inst/cta/include/GCTAResponse.hpp
inst/cta/include/GCTAResponseIrf.hpp
inst/cta/pyext/GCTAObservation.i
inst/cta/src/GCTAObservation.cpp
inst/cta/src/GCTAResponseIrf.cpp
src/obs/GObservations.cpp
-2014-10-31  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-11-07  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * 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
@@ -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
 ---------------------
@@ -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&amp; 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 */
@@ -93,6 +93,7 @@ protected:
                                         const GObservation& obs) const;
     const GCTAInstDir&     retrieve_dir(const std::string& origin,
                                         const GEvent&      event) const;
+
 };
 #endif /* GCTARESPONSE_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 */
@@ -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;
 };
@@ -638,7 +638,10 @@ GEbounds GCTAObservation::ebounds(void) const
  *       ...
  *     </observation>
  *
- * for a binned observation.
+ * for a binned observation. Optional user energy thresholds can be
+ * specified by adding the emin and emax attributes:
+ *
+ *     <observation name="..." id="..." instrument="..." emin="..." emax="...">
  *
  * 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&amp; 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&amp; xml)
  *       ...
  *     </observation>
  *
- * 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:
+ *
+ *     <observation name="..." id="..." instrument="..." emin="..." emax="...">
+ *
  ***************************************************************************/
 void GCTAObservation::write(GXmlElement& xml) const
 {
@@ -813,6 +829,14 @@ void GCTAObservation::write(GXmlElement&amp; 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&amp; 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&amp; 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;
@@ -948,15 +948,47 @@ 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")) {
+
+            // 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&amp; 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;
@@ -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);