Integrate code.

Jürgen Knödlseder
1 parent c7a508a5
Showing 7 changed files with 135 additions and 78 deletions
ChangeLog
NEWS
inst/cta/include/GCTAObservation.hpp
inst/cta/include/GCTAResponseIrf.hpp
inst/cta/pyext/GCTAObservation.i
inst/cta/src/GCTAObservation.cpp
inst/cta/src/GCTAResponseIrf.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,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&amp; 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;
 }
@@ -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 */
@@ -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
@@ -652,14 +655,11 @@ 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"));
     }
@@ -804,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
 {
@@ -824,11 +829,12 @@ 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));
+        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&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) {
@@ -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;
  
@@ -949,9 +949,10 @@ void GCTAResponseIrf::load_aeff(const std::string&amp; 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&amp; 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&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));
  
@@ -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;