Use column names instead of hard-wired indices in GCTAAeff2D (#1656).

Jürgen Knödlseder
1 parent 67894a21
Showing 2 changed files with 118 additions and 67 deletions
inst/cta/include/GCTAAeff2D.hpp
inst/cta/src/GCTAAeff2D.cpp
@@ -87,8 +87,12 @@ private:
     void free_members(void);
     // Members
-    std::string       m_filename;  //!< Name of Aeff response file
-    GCTAResponseTable m_aeff;      //!< Aeff response table
+    std::string       m_filename;      //!< Name of Aeff response file
+    GCTAResponseTable m_aeff;          //!< Aeff response table
+    int               m_inx_energy;    //!< Energy index
+    int               m_inx_theta;     //!< Theta index
+    int               m_inx_aeff;      //!< Effective area (true energy)
+    int               m_inx_aeff_reco; //!< Effective area (reconstructed energy)
 };
 /***************************************************************************
- *                 GCTAAeff2D.cpp - CTA 2D effective area class            *
+ *                GCTAAeff2D.cpp - CTA 2D effective area class             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012-2015 by Juergen Knoedlseder                         *
+ *  copyright (C) 2012-2016 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -56,6 +56,8 @@
 /***********************************************************************//**
  * @brief Void constructor
+ *
+ * Constructs empty effective area.
  ***************************************************************************/
 GCTAAeff2D::GCTAAeff2D(void) : GCTAAeff()
 {
@@ -91,6 +93,8 @@ GCTAAeff2D::GCTAAeff2D(const std::string&amp; filename) : GCTAAeff()
  * @brief Copy constructor
  *
  * @param[in] aeff Effective area.
+ *
+ * Constructs effective area by copying from another effective area.
  ***************************************************************************/
 GCTAAeff2D::GCTAAeff2D(const GCTAAeff2D& aeff) : GCTAAeff(aeff)
 {
@@ -107,6 +111,8 @@ GCTAAeff2D::GCTAAeff2D(const GCTAAeff2D&amp; aeff) : GCTAAeff(aeff)
 /***********************************************************************//**
  * @brief Destructor
+ *
+ * Destructs effective area.
  ***************************************************************************/
 GCTAAeff2D::~GCTAAeff2D(void)
 {
@@ -129,8 +135,10 @@ GCTAAeff2D::~GCTAAeff2D(void)
  *
  * @param[in] aeff Effective area.
  * @return Effective area.
+ *
+ * Assigns effective area.
  ***************************************************************************/
-GCTAAeff2D& GCTAAeff2D::operator= (const GCTAAeff2D& aeff)
+GCTAAeff2D& GCTAAeff2D::operator=(const GCTAAeff2D& aeff)
 {
     // Execute only if object is not identical
     if (this != &aeff) {
@@ -158,14 +166,15 @@ GCTAAeff2D&amp; GCTAAeff2D::operator= (const GCTAAeff2D&amp; aeff)
  * @brief Return effective area in units of cm2
  *
  * @param[in] logE Log10 of the true photon energy (TeV).
- * @param[in] theta Offset angle in camera system (rad). Defaults to 0.0.
- * @param[in] phi Azimuth angle in camera system (rad). Not used in this method.
- * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
- * @param[in] etrue Use true energy (true/false). Defaults to true.
+ * @param[in] theta Offset angle in camera system (rad) (default: 0.0).
+ * @param[in] phi Azimuth angle in camera system (rad).
+ * @param[in] zenith Zenith angle in Earth system (rad).
+ * @param[in] azimuth Azimuth angle in Earth system (rad).
+ * @param[in] etrue Use true energy (default: true).
+ * @return Effective area in cm2.
  *
  * Returns the effective area in units of cm2 for a given energy and
- * offset angle. The effective area is bi-linearily interpolated in the
+ * offset angle. The effective area is bi-linearly interpolated in the
  * log10(energy) - offset angle plane. The method assures that the effective
  * area value never becomes negative.
  *
@@ -182,10 +191,11 @@ double GCTAAeff2D::operator()(const double&amp; logE,
                               const bool&   etrue) const
 {
     // Set parameter index
-    int index = (etrue) ? 0 : 1;
+    int index = (etrue) ? m_inx_aeff : m_inx_aeff_reco;
     // Get effective area value in cm2
-    double aeff = m_aeff(index, logE, theta);
+    double aeff = (m_inx_energy == 0) ? m_aeff(index, logE, theta) :
+                                        m_aeff(index, theta, logE);
     // Make sure that effective area is not negative
     if (aeff < 0.0) {
@@ -204,13 +214,13 @@ double GCTAAeff2D::operator()(const double&amp; logE,
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Clear instance
+ * @brief Clear effective area.
  *
- * This method properly resets the object to an initial state.
+ * Clears effective area.
  ***************************************************************************/
 void GCTAAeff2D::clear(void)
 {
-    // Free class members (base and derived classes, derived class first)
+    // Free class members
     free_members();
     this->GCTAAeff::free_members();
@@ -224,9 +234,11 @@ void GCTAAeff2D::clear(void)
 /***********************************************************************//**
- * @brief Clone instance
+ * @brief Clone effective area
+ *
+ * @return Deep copy of effective area.
  *
- * @return Deep copy of effective area instance.
+ * Returns a pointer to a deep copy of the effective area.
  ***************************************************************************/
 GCTAAeff2D* GCTAAeff2D::clone(void) const
 {
@@ -240,12 +252,26 @@ GCTAAeff2D* GCTAAeff2D::clone(void) const
  * @param[in] table FITS table.
  *
  * @exception GException::invalid_value
- *            FITS file format differs from expectation.
+ *            Response table is not two-dimensional.
+ *
+ * Reads the effective area form the FITS @p table. The following column
+ * names are mandatory:
  *
- * Reads the effective area form the FITS @p table.
+ *     ENERG_LO - Energy lower bin boundaries
+ *     ENERG_HI - Energy upper bin boundaries
+ *     THETA_LO - Offset angle lower bin boundaries
+ *     THETA_HI - Offset angle upper bin boundaries
+ *     EFFAREA  - Effective area
  *
- * The data are stored in m_aeff which is of type GCTAResponseTable. The
- * energy axis will be set to log10, the offset angle axis to radians.
+ * In addition, the following column names are optional:
+ *
+ *     EFFAREA_RECO - Effective area as function of reconstructed energy
+ *
+ * The data are stored in the m_aeff member. The energy axis will be set to
+ * log10, the offset angle axis to radians.
+ *
+ * @todo Analyse the unit of the parameter axis to determine the conversion
+ * factor for the effective areas. For the moment they are hard wired. 
  ***************************************************************************/
 void GCTAAeff2D::read(const GFitsTable& table)
 {
@@ -255,31 +281,34 @@ void GCTAAeff2D::read(const GFitsTable&amp; table)
     // Read effective area table
     m_aeff.read(table);
-    // Check that axis names comply to format
-    if (m_aeff.axis_lo_name(0) != "ENERG_LO" ||
-        m_aeff.axis_hi_name(0) != "ENERG_HI") {
-        std::string msg = "Effective area response table does not contain"
-                          " \"ENERG_LO\" and \"ENERG_HI\" columns as the"
-                          " first axis.";
-        throw GException::invalid_value(G_READ, msg);
-    }
-    if (m_aeff.axis_lo_name(1) != "THETA_LO" ||
-        m_aeff.axis_hi_name(1) != "THETA_HI") {
-        std::string msg = "Effective area response table does not contain"
-                          " \"THETA_LO\" and \"THETA_HI\" columns as the"
-                          " second axis.";
+    // Get mandatory indices (throw exception if not found)
+    m_inx_energy = m_aeff.axis("ENERG");
+    m_inx_theta  = m_aeff.axis("THETA");
+    m_inx_aeff   = m_aeff.table("EFFAREA");
+
+    // Get optional index (use "EFFAREA" if "EFFAREA_RECO" does not exist)
+    m_inx_aeff_reco = (m_aeff.has_table("EFFAREA_RECO")) ?
+                       m_aeff.table("EFFAREA_RECO") : m_inx_aeff;
+
+    // Throw an exception if the table is not two-dimensional
+    if (m_aeff.axes() != 2) {
+        std::string msg = "Expected two-dimensional effective area response "
+                          "table but found "+gammalib::str(m_aeff.axes())+
+                          " dimensions. Please specify a two-dimensional "
+                          "effective area.";
         throw GException::invalid_value(G_READ, msg);
     }
     // Set energy axis to logarithmic scale
-    m_aeff.axis_log10(0);
+    m_aeff.axis_log10(m_inx_energy);
     // Set offset angle axis to radians
-    m_aeff.axis_radians(1);
+    m_aeff.axis_radians(m_inx_theta);
     // Convert effective areas from m2 to cm2
-    m_aeff.scale(0, 1.0e4);
-    m_aeff.scale(1, 1.0e4);
+    for (int i = 0; i < m_aeff.tables(); ++i) {
+        m_aeff.scale(i, 1.0e4);
+    }
     // Return
     return;
@@ -291,7 +320,7 @@ void GCTAAeff2D::read(const GFitsTable&amp; table)
  *
  * @param[in] table FITS binary table.
  *
- * Writes effective area into the FITS binary @p table.
+ * Writes effective area into a FITS binary @p table.
  *
  * @todo Add keywords.
  ***************************************************************************/
@@ -301,8 +330,9 @@ void GCTAAeff2D::write(GFitsBinTable&amp; table) const
     GCTAResponseTable aeff(m_aeff);
     // Convert area from cm2 to m2
-    aeff.scale(0, 1.0e-4);
-    aeff.scale(1, 1.0e-4);
+    for (int i = 0; i < aeff.tables(); ++i) {
+        aeff.scale(i, 1.0e-4);
+    }
     // Write response table
     aeff.write(table);
@@ -384,44 +414,51 @@ void GCTAAeff2D::save(const std::string&amp; filename, const bool&amp; clobber) const
 /***********************************************************************//**
- * @brief Return maximum effective area at a given energy
+ * @brief Return maximum effective area at a given energy in cm2
  *
  * @param[in] logE Log10 of the true photon energy (TeV).
- * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
- * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
- * @param[in] etrue Use true energy (true/false). Not used.
+ * @param[in] zenith Zenith angle in Earth system (rad).
+ * @param[in] azimuth Azimuth angle in Earth system (rad).
+ * @param[in] etrue Use true energy (default: true).
  * @return Maximum effective area (cm2).
+ *
+ * Returns the maximum effective area for a given energy, zenith and azimuth
+ * angle in units of cm2.
  ***************************************************************************/
 double GCTAAeff2D::max(const double& logE,
                        const double& zenith,
                        const double& azimuth,
                        const bool&   etrue) const
 {
-    // Set parameter index
-    int index = (etrue) ? 0 : 1;
-
-    // Get number of theta bins
-    int n_theta = m_aeff.axis(1);
-
     // Initialise maximum effective area
     double max_aeff = 0.0;
-    // Compute stepsize
-    double theta_step = m_aeff.axis_hi(1, n_theta - 1) / (double(n_theta) - 1.0);
+    // Set parameter index
+    int index = (etrue) ? m_inx_aeff : m_inx_aeff_reco;
+
+    // Get number of theta bins
+    int n_theta = m_aeff.axis_bins(m_inx_theta);
     // Loop over theta values
     for (int i = 0; i < n_theta; ++i) {
-        // Compute theta
-        double theta = i * theta_step;
+        // Compute lower and upper theta bin values
+        double theta_lo = m_aeff.axis_lo(m_inx_theta, i);
+        double theta_hi = m_aeff.axis_hi(m_inx_theta, i);
         // Get effective area value in cm2
-        double aeff = m_aeff(index, logE, theta);
+        double aeff_lo = (m_inx_energy == 0) ? m_aeff(index, logE, theta_lo) :
+                                               m_aeff(index, theta_lo, logE);
+        double aeff_hi = (m_inx_energy == 0) ? m_aeff(index, logE, theta_hi) :
+                                               m_aeff(index, theta_hi, logE);
         // Update maximum effective area if larger than current maximum
         // effective area
-        if (aeff > max_aeff) {
-            max_aeff = aeff;
+        if (aeff_lo > max_aeff) {
+            max_aeff = aeff_lo;
+        }
+        if (aeff_hi > max_aeff) {
+            max_aeff = aeff_hi;
         }
     } // endfor: loop over theta values
@@ -446,12 +483,14 @@ std::string GCTAAeff2D::print(const GChatter&amp; chatter) const
     if (chatter != SILENT) {
         // Compute energy boundaries in TeV
-        double emin = m_aeff.axis_lo(0,0);
-        double emax = m_aeff.axis_hi(0,m_aeff.axis(0)-1);
+        double emin = m_aeff.axis_lo(m_inx_energy,0);
+        double emax = m_aeff.axis_hi(m_inx_energy,
+                                     m_aeff.axis_bins(m_inx_energy)-1);
         // Compute offset angle boundaries in deg
-        double omin = m_aeff.axis_lo(1,0);
-        double omax = m_aeff.axis_hi(1,m_aeff.axis(1)-1);
+        double omin = m_aeff.axis_lo(m_inx_theta,0);
+        double omax = m_aeff.axis_hi(m_inx_theta,
+                                     m_aeff.axis_bins(m_inx_theta)-1);
         // Append header
         result.append("=== GCTAAeff2D ===");
@@ -459,9 +498,9 @@ std::string GCTAAeff2D::print(const GChatter&amp; chatter) const
         // Append information
         result.append("\n"+gammalib::parformat("Filename")+m_filename);
         result.append("\n"+gammalib::parformat("Number of energy bins") +
-                      gammalib::str(m_aeff.axis(0)));
+                      gammalib::str(m_aeff.axis_bins(m_inx_energy)));
         result.append("\n"+gammalib::parformat("Number of offset bins") +
-                      gammalib::str(m_aeff.axis(1)));
+                      gammalib::str(m_aeff.axis_bins(m_inx_theta)));
         result.append("\n"+gammalib::parformat("Log10(Energy) range"));
         result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
         result.append("\n"+gammalib::parformat("Offset angle range"));
@@ -488,6 +527,10 @@ void GCTAAeff2D::init_members(void)
     // Initialise members
     m_filename.clear();
     m_aeff.clear();
+    m_inx_energy    = 0;
+    m_inx_theta     = 1;
+    m_inx_aeff      = 0;
+    m_inx_aeff_reco = 1;
     // Return
     return;
@@ -502,8 +545,12 @@ void GCTAAeff2D::init_members(void)
 void GCTAAeff2D::copy_members(const GCTAAeff2D& aeff)
 {
     // Copy members
-    m_filename = aeff.m_filename;
-    m_aeff     = aeff.m_aeff;
+    m_filename      = aeff.m_filename;
+    m_aeff          = aeff.m_aeff;
+    m_inx_energy    = aeff.m_inx_energy;
+    m_inx_theta     = aeff.m_inx_theta;
+    m_inx_aeff      = aeff.m_inx_aeff;
+    m_inx_aeff_reco = aeff.m_inx_aeff_reco;
     // Return
     return;