From 9878adf96920dd5f2e95fb3ccd7ba1c65ca32ee6 Mon Sep 17 00:00:00 2001
From: knodlseder <>
Date: Fri, 9 Mar 2012 16:34:08 +0000
Subject: [PATCH] Add option to read ANGRES40 column in GCTAResponse::read_psf.

---
 inst/cta/src/GCTAResponse.cpp | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/inst/cta/src/GCTAResponse.cpp b/inst/cta/src/GCTAResponse.cpp
index c8a84c8..0aa9fd1 100644
--- a/inst/cta/src/GCTAResponse.cpp
+++ b/inst/cta/src/GCTAResponse.cpp
@@ -833,6 +833,12 @@ void GCTAResponse::read_arf(const GFitsTable* hdu)
  * energies are converted to TeV. Conversion is done based on the units
  * provided for the energy columns. Units that are recognized are 'keV',
  * 'MeV', 'GeV', and 'TeV' (case independent).
+ *
+ * The Gaussian width parameter may be either given in 68% containment radius
+ * (R68) or in sigma (ANGRES40; corresponding to a 38% containment radius).
+ * The former format is used for CTA and H.E.S.S. data, the latter for
+ * MAGIC data. All these things should be more uniform once we have a well
+ * defined format.
  ***************************************************************************/
 void GCTAResponse::read_psf(const GFitsTable* hdu)
 {
@@ -843,7 +849,17 @@ void GCTAResponse::read_psf(const GFitsTable* hdu)
     // Get pointers to table columns
     const GFitsTableCol* energy_lo = &(*hdu)["ENERG_LO"];
     const GFitsTableCol* energy_hi = &(*hdu)["ENERG_HI"];
-    const GFitsTableCol* r68       = &(*hdu)["R68"];
+
+    // Handle various data formats (H.E.S.S. and MAGIC)
+    const GFitsTableCol* r68;
+    double               r68_scale = 1.0;
+    if (hdu->hascolumn("R68")) {
+        r68 = &(*hdu)["R68"];
+    }
+    else {
+        r68 = &(*hdu)["ANGRES40"];
+        r68_scale = 1.0 / 0.6624305; // MAGIC PSF is already 1 sigma
+    }
 
     // Determine unit conversion factors (default: TeV)
     std::string u_energy_lo = tolower(strip_whitespace(energy_lo->unit()));
@@ -880,8 +896,8 @@ void GCTAResponse::read_psf(const GFitsTable* hdu)
         double e_max = energy_hi->real(i) * c_energy_hi;
         double logE  = 0.5 * (log10(e_min) + log10(e_max));
         
-        // Extract r68 value
-        double r68_value = r68->real(i);
+        // Extract r68 value and scale as required
+        double r68_value = r68->real(i) * r68_scale;
         
         // Store log10 mean energy and r68 value
         m_psf_logE.append(logE);
--
libgit2 0.21.4