Johan Bregeon / gammalib

Commit 9878adf96920dd5f2e95fb3ccd7ba1c65ca32ee6

Authored by knodlseder
1 parent bff928b2

Add option to read ANGRES40 column in GCTAResponse::read_psf. 

This is needed to support the preliminary MAGIC data response format.
Inline Side-by-side
Showing 1 changed file with 19 additions and 3 deletions
inst/cta/src/GCTAResponse.cpp
  View file @9878adf
... ... @@ -833,6 +833,12 @@ void GCTAResponse::read_arf(const GFitsTable* hdu)
833 833 * energies are converted to TeV. Conversion is done based on the units
834 834 * provided for the energy columns. Units that are recognized are 'keV',
835 835 * 'MeV', 'GeV', and 'TeV' (case independent).
  836 + *
  837 + * The Gaussian width parameter may be either given in 68% containment radius
  838 + * (R68) or in sigma (ANGRES40; corresponding to a 38% containment radius).
  839 + * The former format is used for CTA and H.E.S.S. data, the latter for
  840 + * MAGIC data. All these things should be more uniform once we have a well
  841 + * defined format.
836 842 ***************************************************************************/
837 843 void GCTAResponse::read_psf(const GFitsTable* hdu)
838 844 {
... ... @@ -843,7 +849,17 @@ void GCTAResponse::read_psf(const GFitsTable* hdu)
843 849 // Get pointers to table columns
844 850 const GFitsTableCol* energy_lo = &(*hdu)["ENERG_LO"];
845 851 const GFitsTableCol* energy_hi = &(*hdu)["ENERG_HI"];
846   - const GFitsTableCol* r68 = &(*hdu)["R68"];
  852 +
  853 + // Handle various data formats (H.E.S.S. and MAGIC)
  854 + const GFitsTableCol* r68;
  855 + double r68_scale = 1.0;
  856 + if (hdu->hascolumn("R68")) {
  857 + r68 = &(*hdu)["R68"];
  858 + }
  859 + else {
  860 + r68 = &(*hdu)["ANGRES40"];
  861 + r68_scale = 1.0 / 0.6624305; // MAGIC PSF is already 1 sigma
  862 + }
847 863  
848 864 // Determine unit conversion factors (default: TeV)
849 865 std::string u_energy_lo = tolower(strip_whitespace(energy_lo->unit()));
... ... @@ -880,8 +896,8 @@ void GCTAResponse::read_psf(const GFitsTable* hdu)
880 896 double e_max = energy_hi->real(i) * c_energy_hi;
881 897 double logE = 0.5 * (log10(e_min) + log10(e_max));
882 898  
883   - // Extract r68 value
884   - double r68_value = r68->real(i);
  899 + // Extract r68 value and scale as required
  900 + double r68_value = r68->real(i) * r68_scale;
885 901  
886 902 // Store log10 mean energy and r68 value
887 903 m_psf_logE.append(logE);
... ...