Justine Devin / gammalib

Commit e962fb3f67003a695a4038a6252a9bce5dd4310e

Authored by Jürgen Knödlseder
1 parent f606364b

Use column names instead of hard-wired indices in CTA response table instrument response functions.

Inline Side-by-side
Showing 10 changed files with 517 additions and 298 deletions
ChangeLog
  View file @e962fb3
1   -2016-02-12 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-02-13 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 =========================
5 5  
  6 + Use column names for CTA response table access (#1656)
6 7 Add non-const GObservation::events() method (#1664)
7 8 Remove event fetching from GCTAObservation class (#1648)
8 9 Implement event fetching and disposal in GCTAEventList class (#1648)
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 1.1.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -12 February 2016
  5 +13 February 2016
6 6  
7 7  
8 8 1. Introduction
... ... @@ -15,6 +15,14 @@ made since the last release of GammaLib.
15 15 ---------------------------------
16 16 The following methods have been added to the public interface:
17 17 - GObservation::events() (non-const variant)
  18 +- GCTAResponseTable::axis(std::string)
  19 +- GCTAResponseTable::table(std::string)
  20 +- GCTAResponseTable::has_axis(std::string)
  21 +- GCTAResponseTable::has_table(std::string)
  22 +
  23 +The following methods have been renamed:
  24 +- GCTAResponseTable::tables(), was formerly GCTAResponseTable::size()
  25 +- GCTAResponseTable::axis_bins(), was formerly GCTAResponseTable::axis()
18 26  
19 27 The return value has been changed for the following methods:
20 28 - GFits::filename() returns now GFilename instead of std::string
... ... @@ -116,6 +124,10 @@ None
116 124  
117 125 20. CTA interface
118 126 -----------------
  127 +The response table instrument response functions now rely on the FITS table
  128 +column names to determine the information instead of hard-wired column indices
  129 +(#1656).
  130 +
119 131 Event fetching and disposal has been implemented in the GCTAEventList class by
120 132 additing the fetch() and dispose() methods. Now, events are only loaded upon
121 133 request and can be dispatched after user to save memory. The event fetching and
... ...
inst/cta/include/GCTABackground3D.hpp
  View file @e962fb3
1 1 /***************************************************************************
2 2 * GCTABackground3D.hpp - CTA 3D background class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -92,6 +92,10 @@ private:
92 92 GCTAResponseTable m_background; //!< Background response table
93 93 double m_mc_max_bin; //!< Maximum spatial binsize for MC
94 94 double m_mc_max_logE; //!< Maximum log energy binsize for MC
  95 + int m_inx_detx; //!< DETX index
  96 + int m_inx_dety; //!< DETY index
  97 + int m_inx_energy; //!< Energy index
  98 + int m_inx_bgd; //!< Background index
95 99  
96 100 // Monte Carlo cache
97 101 mutable std::vector<double> m_mc_cache; //!< Monte Carlo cache
... ...
inst/cta/include/GCTAEdisp2D.hpp
  View file @e962fb3
1 1 /***************************************************************************
2 2 * GCTAEdisp2D.hpp - CTA 2D energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015 by Florent Forest *
  4 + * copyright (C) 2015-2016 by Florent Forest *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -116,22 +116,26 @@ private:
116 116 // Protected classes
117 117 class edisp_kern : public GFunction {
118 118 public:
119   - edisp_kern(const GCTAEdisp2D* parent,
120   - const double& logEsrc,
121   - const double& theta) :
  119 + edisp_kern(const GCTAEdisp2D* parent,
  120 + const double& logEsrc,
  121 + const double& theta) :
122 122 m_parent(parent),
123 123 m_logEsrc(logEsrc),
124 124 m_theta(theta) { }
125 125 double eval(const double& x);
126 126 protected:
127   - const GCTAEdisp2D* m_parent; //!< Pointer to parent class
128   - double m_logEsrc; //!< True photon energy
129   - double m_theta; //!< Offset angle
  127 + const GCTAEdisp2D* m_parent; //!< Pointer to parent class
  128 + double m_logEsrc; //!< True photon energy
  129 + double m_theta; //!< Offset angle
130 130 };
131 131  
132 132 // Members
133   - std::string m_filename; //!< Name of Edisp response file
134   - GCTAResponseTable m_edisp; //!< Edisp response table
  133 + std::string m_filename; //!< Name of Edisp response file
  134 + GCTAResponseTable m_edisp; //!< Edisp response table
  135 + int m_inx_etrue; //!< True energy index
  136 + int m_inx_migra; //!< Migration index
  137 + int m_inx_theta; //!< Theta index
  138 + int m_inx_matrix; //!< Matrix
135 139  
136 140 // Computation cache
137 141 mutable bool m_ebounds_obs_computed;
... ...
inst/cta/include/GCTAPsf2D.hpp
  View file @e962fb3
... ... @@ -112,6 +112,13 @@ private:
112 112 // Members
113 113 std::string m_filename; //!< Name of Aeff response file
114 114 GCTAResponseTable m_psf; //!< PSF response table
  115 + int m_inx_energy; //!< Energy index
  116 + int m_inx_theta; //!< Theta index
  117 + int m_inx_sigma1; //!< 1st Gaussian sigma
  118 + int m_inx_ampl2; //!< 2nd Gaussian relative amplitude
  119 + int m_inx_sigma2; //!< 2nd Gaussian sigma
  120 + int m_inx_ampl3; //!< 3nd Gaussian relative amplitude
  121 + int m_inx_sigma3; //!< 3nd Gaussian sigma
115 122  
116 123 // Precomputation cache
117 124 mutable double m_par_logE; //!< Cache energy
... ...
inst/cta/include/GCTAPsfKing.hpp
  View file @e962fb3
1 1 /***************************************************************************
2 2 * GCTAPsfKing.hpp - King profile CTA point spread function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Michael Mayer *
  4 + * copyright (C) 2013-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -111,6 +111,10 @@ private:
111 111 // Members
112 112 std::string m_filename; //!< Name of Aeff response file
113 113 GCTAResponseTable m_psf; //!< PSF response table
  114 + int m_inx_energy; //!< Energy index
  115 + int m_inx_theta; //!< Theta index
  116 + int m_inx_gamma; //!< Gamma
  117 + int m_inx_sigma; //!< Sigma
114 118  
115 119 // Evaluation cache
116 120 mutable double m_par_logE; //!< Cache energy
... ...
inst/cta/src/GCTABackground3D.cpp
  View file @e962fb3
... ... @@ -61,6 +61,8 @@
61 61  
62 62 /***********************************************************************//**
63 63 * @brief Void constructor
  64 + *
  65 + * Constructs empty background.
64 66 ***************************************************************************/
65 67 GCTABackground3D::GCTABackground3D(void) : GCTABackground()
66 68 {
... ... @@ -97,6 +99,8 @@ GCTABackground3D::GCTABackground3D(const std::string&amp; filename) :
97 99 * @brief Copy constructor
98 100 *
99 101 * @param[in] bgd Background.
  102 + *
  103 + * Constructs background by copying from another background.
100 104 ***************************************************************************/
101 105 GCTABackground3D::GCTABackground3D(const GCTABackground3D& bgd) :
102 106 GCTABackground(bgd)
... ... @@ -114,6 +118,8 @@ GCTABackground3D::GCTABackground3D(const GCTABackground3D&amp; bgd) :
114 118  
115 119 /***********************************************************************//**
116 120 * @brief Destructor
  121 + *
  122 + * Destructs point spread function.
117 123 ***************************************************************************/
118 124 GCTABackground3D::~GCTABackground3D(void)
119 125 {
... ... @@ -136,6 +142,8 @@ GCTABackground3D::~GCTABackground3D(void)
136 142 *
137 143 * @param[in] bgd Background.
138 144 * @return Background.
  145 + *
  146 + * Assigns background.
139 147 ***************************************************************************/
140 148 GCTABackground3D& GCTABackground3D::operator=(const GCTABackground3D& bgd)
141 149 {
... ... @@ -190,9 +198,9 @@ double GCTABackground3D::operator()(const double&amp; logE,
190 198 // Get background rate
191 199 #if defined(G_LOG_INTERPOLATION)
192 200 // Retrieve node arrays
193   - const GNodeArray& detx_nodes = m_background.axis_nodes(0);
194   - const GNodeArray& dety_nodes = m_background.axis_nodes(1);
195   - const GNodeArray& energy_nodes = m_background.axis_nodes(2);
  201 + const GNodeArray& detx_nodes = m_background.axis_nodes(m_inx_detx);
  202 + const GNodeArray& dety_nodes = m_background.axis_nodes(m_inx_dety);
  203 + const GNodeArray& energy_nodes = m_background.axis_nodes(m_inx_energy);
196 204  
197 205 // Set values for node arrays
198 206 detx_nodes.set_value(detx);
... ... @@ -200,7 +208,7 @@ double GCTABackground3D::operator()(const double&amp; logE,
200 208 energy_nodes.set_value(logE);
201 209  
202 210 // Compute offsets of DETY in DETX-DETY plane
203   - int size1 = m_background.axis(0);
  211 + int size1 = m_background.axis_bins(m_inx_detx);
204 212 int offset_left = dety_nodes.inx_left() * size1;
205 213 int offset_right = dety_nodes.inx_right() * size1;
206 214  
... ... @@ -225,20 +233,20 @@ double GCTABackground3D::operator()(const double&amp; logE,
225 233 double wgt_emax = energy_nodes.wgt_right();
226 234  
227 235 // Compute offsets in energy dimension
228   - int npixels = m_background.axis_bins(0) *
229   - m_background.axis_bins(1);
  236 + int npixels = m_background.axis_bins(m_inx_detx) *
  237 + m_background.axis_bins(m_inx_dety);
230 238 int offset_emin = inx_emin * npixels;
231 239 int offset_emax = inx_emax * npixels;
232 240  
233 241 // Bi-linear interpolation the rates in both energy layers
234   - double rate_emin = wgt_ll * m_background(0, inx_ll + offset_emin) +
235   - wgt_lr * m_background(0, inx_lr + offset_emin) +
236   - wgt_rl * m_background(0, inx_rl + offset_emin) +
237   - wgt_rr * m_background(0, inx_rr + offset_emin);
238   - double rate_emax = wgt_ll * m_background(0, inx_ll + offset_emax) +
239   - wgt_lr * m_background(0, inx_lr + offset_emax) +
240   - wgt_rl * m_background(0, inx_rl + offset_emax) +
241   - wgt_rr * m_background(0, inx_rr + offset_emax);
  242 + double rate_emin = wgt_ll * m_background(m_inx_bgd, inx_ll + offset_emin) +
  243 + wgt_lr * m_background(m_inx_bgd, inx_lr + offset_emin) +
  244 + wgt_rl * m_background(m_inx_bgd, inx_rl + offset_emin) +
  245 + wgt_rr * m_background(m_inx_bgd, inx_rr + offset_emin);
  246 + double rate_emax = wgt_ll * m_background(m_inx_bgd, inx_ll + offset_emax) +
  247 + wgt_lr * m_background(m_inx_bgd, inx_lr + offset_emax) +
  248 + wgt_rl * m_background(m_inx_bgd, inx_rl + offset_emax) +
  249 + wgt_rr * m_background(m_inx_bgd, inx_rr + offset_emax);
242 250  
243 251 // If both rates are positive then perform a logarithmic
244 252 // interpolation in energy
... ... @@ -252,7 +260,11 @@ double GCTABackground3D::operator()(const double&amp; logE,
252 260 rate = wgt_emin * rate_emin + wgt_emax * rate_emax;
253 261 }
254 262 #else
255   - rate = m_background(0, detx, dety, logE);
  263 + double arg[3];
  264 + arg[m_inx_detx] = detx;
  265 + arg[m_inx_dety] = dety;
  266 + arg[m_inx_energy] = logE;
  267 + rate = m_background(m_inx_bgd, arg[0], arg[1], arg[2]);
256 268 #endif
257 269  
258 270 // Make sure that background rate is not negative
... ... @@ -276,11 +288,11 @@ double GCTABackground3D::operator()(const double&amp; logE,
276 288 /***********************************************************************//**
277 289 * @brief Clear background
278 290 *
279   - * This method properly resets the background to an initial state.
  291 + * Clears background.
280 292 ***************************************************************************/
281 293 void GCTABackground3D::clear(void)
282 294 {
283   - // Free class members (base and derived classes, derived class first)
  295 + // Free class members
284 296 free_members();
285 297 this->GCTABackground::free_members();
286 298  
... ... @@ -296,7 +308,9 @@ void GCTABackground3D::clear(void)
296 308 /***********************************************************************//**
297 309 * @brief Clone background
298 310 *
299   - * @return Pointer to deep copy of background.
  311 + * @return Deep copy of background.
  312 + *
  313 + * Returns a pointer to a deep copy of the background.
300 314 ***************************************************************************/
301 315 GCTABackground3D* GCTABackground3D::clone(void) const
302 316 {
... ... @@ -310,13 +324,21 @@ GCTABackground3D* GCTABackground3D::clone(void) const
310 324 * @param[in] table FITS table.
311 325 *
312 326 * @exception GException::invalid_value
313   - * FITS file format differs from expectation.
  327 + * Response table is not three-dimensional.
314 328 *
315   - * Reads the background from the FITS @p table.
  329 + * Reads the background form the FITS @p table. The following column names
  330 + * are mandatory:
316 331 *
317   - * The data are stored in m_background which is of type GCTAResponseTable.
318   - * The DETX and DETY axes will be set to radians, the energy axis will be set
319   - * to log10.
  332 + * DETX_LO - DETX lower bin boundaries
  333 + * DETX_HI - DETX upper bin boundaries
  334 + * DETY_LO - DETY lower bin boundaries
  335 + * DETY_HI - DETY upper bin boundaries
  336 + * ENERG_LO - Energy lower bin boundaries
  337 + * ENERG_HI - Energy upper bin boundaries
  338 + * BGD - Background template
  339 + *
  340 + * The data are stored in the m_background member. The DETX and DETY axes
  341 + * will be set to radians, the energy axis will be set to log10.
320 342 ***************************************************************************/
321 343 void GCTABackground3D::read(const GFitsTable& table)
322 344 {
... ... @@ -326,35 +348,28 @@ void GCTABackground3D::read(const GFitsTable&amp; table)
326 348 // Read background table
327 349 m_background.read(table);
328 350  
329   - // Check that axis names comply to format
330   - if (m_background.axis_lo_name(0) != "DETX_LO" ||
331   - m_background.axis_hi_name(0) != "DETX_HI") {
332   - std::string msg = "Background response table does not contain"
333   - " \"DETX_LO\" and \"DETX_HI\" columns as the"
334   - " first axis.";
335   - throw GException::invalid_value(G_READ, msg);
336   - }
337   - if (m_background.axis_lo_name(1) != "DETY_LO" ||
338   - m_background.axis_hi_name(1) != "DETY_HI") {
339   - std::string msg = "Background response table does not contain"
340   - " \"DETY_LO\" and \"DETY_HI\" columns as the"
341   - " second axis.";
342   - throw GException::invalid_value(G_READ, msg);
343   - }
344   - if (m_background.axis_lo_name(2) != "ENERG_LO" ||
345   - m_background.axis_hi_name(2) != "ENERG_HI") {
346   - std::string msg = "Background response table does not contain"
347   - " \"ENERG_LO\" and \"ENERG_HI\" columns as the"
348   - " third axis.";
  351 + // Get mandatory indices (throw exception if not found)
  352 + m_inx_detx = m_background.axis("DETX");
  353 + m_inx_dety = m_background.axis("DETY");
  354 + m_inx_energy = m_background.axis("ENERG");
  355 + m_inx_bgd = m_background.table("BGD");
  356 +
  357 + // Throw an exception if the table is not three-dimensional
  358 + if (m_background.axes() != 3) {
  359 + std::string msg = "Expected three-dimensional background "
  360 + "response table but found "+
  361 + gammalib::str(m_background.axes())+
  362 + " dimensions. Please specify a three-dimensional "
  363 + "background.";
349 364 throw GException::invalid_value(G_READ, msg);
350 365 }
351 366  
352 367 // Set DETX and DETY axis to radians
353   - m_background.axis_radians(0);
354   - m_background.axis_radians(1);
  368 + m_background.axis_radians(m_inx_detx);
  369 + m_background.axis_radians(m_inx_dety);
355 370  
356 371 // Set energy axis to logarithmic scale
357   - m_background.axis_log10(2);
  372 + m_background.axis_log10(m_inx_energy);
358 373  
359 374 // Return
360 375 return;
... ... @@ -490,15 +505,17 @@ GCTAInstDir GCTABackground3D::mc(const GEnergy&amp; energy,
490 505 index = int((logE - m_mc_logE_min) / m_mc_logE_bin);
491 506 }
492 507 if (index < 0 || index >= m_mc_nmaps) {
493   - std::string msg = "The specified energy "+energy.print()+" does"
494   - " not fall in any of the energy boundaries of"
495   - " the background response cube.\n"
496   - "Please select an energy that is comprised"
497   - " in the range ["+
498   - gammalib::str(m_background.axis_lo(2,0))+
499   - " "+m_background.axis_lo_unit(2)+" - "+
500   - gammalib::str(m_background.axis_hi(2,m_background.axis_bins(2)-1))+
501   - " "+m_background.axis_hi_unit(2)+"].";
  508 + double emin = m_background.axis_lo(m_inx_energy,0);
  509 + double emax = m_background.axis_hi(m_inx_energy,
  510 + m_background.axis_bins(m_inx_energy)-1);
  511 + std::string msg = "The specified energy "+energy.print()+" does "
  512 + "not fall in any of the energy boundaries of "
  513 + "the background response cube. Please select "
  514 + "an energy that is comprised in the range ["+
  515 + gammalib::str(emin)+" "+
  516 + m_background.axis_lo_unit(m_inx_energy)+" - "+
  517 + gammalib::str(emax)+" "+
  518 + m_background.axis_hi_unit(m_inx_energy)+"].";
502 519 throw GException::invalid_value(G_MC, msg);
503 520 }
504 521  
... ... @@ -564,24 +581,27 @@ std::string GCTABackground3D::print(const GChatter&amp; chatter) const
564 581 if (is_valid()) {
565 582  
566 583 // Compute DETX boundaries in deg
567   - double detx_min = m_background.axis_lo(0,0);
568   - double detx_max = m_background.axis_hi(0,m_background.axis_bins(0)-1);
  584 + double detx_min = m_background.axis_lo(m_inx_detx,0);
  585 + double detx_max = m_background.axis_hi(m_inx_detx,
  586 + m_background.axis_bins(m_inx_detx)-1);
569 587  
570 588 // Compute DETY boundaries in deg
571   - double dety_min = m_background.axis_lo(1,0);
572   - double dety_max = m_background.axis_hi(1,m_background.axis_bins(1)-1);
  589 + double dety_min = m_background.axis_lo(m_inx_dety,0);
  590 + double dety_max = m_background.axis_hi(m_inx_dety,
  591 + m_background.axis_bins(m_inx_dety)-1);
573 592  
574 593 // Compute energy boundaries in TeV
575   - double emin = m_background.axis_lo(2,0);
576   - double emax = m_background.axis_hi(2,m_background.axis_bins(2)-1);
  594 + double emin = m_background.axis_lo(m_inx_energy,0);
  595 + double emax = m_background.axis_hi(m_inx_energy,
  596 + m_background.axis_bins(m_inx_energy)-1);
577 597  
578 598 // Append information
579 599 result.append("\n"+gammalib::parformat("Number of DETX bins") +
580   - gammalib::str(m_background.axis_bins(0)));
  600 + gammalib::str(m_background.axis_bins(m_inx_detx)));
581 601 result.append("\n"+gammalib::parformat("Number of DETY bins") +
582   - gammalib::str(m_background.axis_bins(1)));
  602 + gammalib::str(m_background.axis_bins(m_inx_dety)));
583 603 result.append("\n"+gammalib::parformat("Number of energy bins") +
584   - gammalib::str(m_background.axis_bins(2)));
  604 + gammalib::str(m_background.axis_bins(m_inx_energy)));
585 605 result.append("\n"+gammalib::parformat("DETX range"));
586 606 result.append(gammalib::str(detx_min)+" - "+gammalib::str(detx_max)+" deg");
587 607 result.append("\n"+gammalib::parformat("DETX range"));
... ... @@ -628,6 +648,10 @@ void GCTABackground3D::init_members(void)
628 648 m_background.clear();
629 649 m_mc_max_bin = 0.05; //!< Spatial binning not worse than 0.05 deg
630 650 m_mc_max_logE = 0.02; //!< Spectral binning not worse than 0.02^10 TeV
  651 + m_inx_detx = 0;
  652 + m_inx_dety = 1;
  653 + m_inx_energy = 2;
  654 + m_inx_bgd = 0;
631 655  
632 656 // Initialise MC cache
633 657 m_mc_cache.clear();
... ... @@ -663,6 +687,10 @@ void GCTABackground3D::copy_members(const GCTABackground3D&amp; bgd)
663 687 m_background = bgd.m_background;
664 688 m_mc_max_bin = bgd.m_mc_max_bin;
665 689 m_mc_max_logE = bgd.m_mc_max_logE;
  690 + m_inx_detx = bgd.m_inx_detx;
  691 + m_inx_dety = bgd.m_inx_dety;
  692 + m_inx_energy = bgd.m_inx_energy;
  693 + m_inx_bgd = bgd.m_inx_bgd;
666 694  
667 695 // Copy MC cache
668 696 m_mc_cache = bgd.m_mc_cache;
... ... @@ -724,14 +752,14 @@ void GCTABackground3D::init_mc_cache(void) const
724 752 }
725 753  
726 754 // Determine number of response cube pixels and maps
727   - m_mc_nx = m_background.axis_bins(0);
728   - m_mc_ny = m_background.axis_bins(1);
729   - m_mc_nmaps = m_background.axis_bins(2);
  755 + m_mc_nx = m_background.axis_bins(m_inx_detx);
  756 + m_mc_ny = m_background.axis_bins(m_inx_dety);
  757 + m_mc_nmaps = m_background.axis_bins(m_inx_energy);
730 758  
731 759 // Compute DETX boundaries and binsize in deg. If DETX binsize is too
732 760 // coarse then request a finer sampling
733   - m_mc_detx_min = m_background.axis_lo(0,0);
734   - m_mc_detx_max = m_background.axis_hi(0,m_mc_nx-1);
  761 + m_mc_detx_min = m_background.axis_lo(m_inx_detx,0);
  762 + m_mc_detx_max = m_background.axis_hi(m_inx_detx,m_mc_nx-1);
735 763 m_mc_detx_bin = (m_mc_detx_max - m_mc_detx_min) / m_mc_nx;
736 764 if (m_mc_detx_bin > m_mc_max_bin) {
737 765 m_mc_nx = int((m_mc_detx_max - m_mc_detx_min) / m_mc_max_bin + 1.0);
... ... @@ -740,8 +768,8 @@ void GCTABackground3D::init_mc_cache(void) const
740 768  
741 769 // Compute DETY boundaries and binsize in deg. If DETY binsize is too
742 770 // coarse then request a finer sampling
743   - m_mc_dety_min = m_background.axis_lo(1,0);
744   - m_mc_dety_max = m_background.axis_hi(1,m_mc_ny-1);
  771 + m_mc_dety_min = m_background.axis_lo(m_inx_dety,0);
  772 + m_mc_dety_max = m_background.axis_hi(m_inx_dety,m_mc_ny-1);
745 773 m_mc_dety_bin = (m_mc_dety_max - m_mc_dety_min) / m_mc_ny;
746 774 if (m_mc_dety_bin > m_mc_max_bin) {
747 775 m_mc_ny = int((m_mc_dety_max - m_mc_dety_min) / m_mc_max_bin + 1.0);
... ... @@ -750,8 +778,10 @@ void GCTABackground3D::init_mc_cache(void) const
750 778  
751 779 // Compute energy sampling. If energy sampling is too coarse then
752 780 // request a finer sampling
753   - GEnergy emin(m_background.axis_lo(2,0), m_background.axis_lo_unit(2));
754   - GEnergy emax(m_background.axis_hi(2,m_mc_nmaps-1), m_background.axis_hi_unit(2));
  781 + GEnergy emin(m_background.axis_lo(m_inx_energy,0),
  782 + m_background.axis_lo_unit(m_inx_energy));
  783 + GEnergy emax(m_background.axis_hi(m_inx_energy,m_mc_nmaps-1),
  784 + m_background.axis_hi_unit(m_inx_energy));
755 785 m_mc_logE_min = emin.log10TeV();
756 786 m_mc_logE_max = emax.log10TeV();
757 787 m_mc_logE_bin = (m_mc_logE_max - m_mc_logE_min) / m_mc_nmaps;
... ... @@ -820,7 +850,7 @@ void GCTABackground3D::init_mc_cache(void) const
820 850  
821 851 // Set energy value (unit independent)
822 852 GEnergy energy;
823   - energy.log10(logE, m_background.axis_lo_unit(2));
  853 + energy.log10(logE, m_background.axis_lo_unit(m_inx_energy));
824 854  
825 855 // Only append node if rate > 0
826 856 if (total_rate > 0.0) {
... ...
inst/cta/src/GCTAEdisp2D.cpp
  View file @e962fb3
1 1 /***************************************************************************
2 2 * GCTAEdisp2D.cpp - CTA 2D energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015 by Florent Forest *
  4 + * copyright (C) 2015-2016 by Florent Forest *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -40,6 +40,7 @@
40 40 #include "GCTAEdisp2D.hpp"
41 41  
42 42 /* __ Method name definitions ____________________________________________ */
  43 +#define G_READ "GCTAEdisp2D::read(GFitsTable&)"
43 44  
44 45 /* __ Macros _____________________________________________________________ */
45 46  
... ... @@ -58,6 +59,8 @@
58 59  
59 60 /***********************************************************************//**
60 61 * @brief Void constructor
  62 + *
  63 + * Constructs empty energy dispersion.
61 64 ***************************************************************************/
62 65 GCTAEdisp2D::GCTAEdisp2D(void) : GCTAEdisp()
63 66 {
... ... @@ -94,6 +97,8 @@ GCTAEdisp2D::GCTAEdisp2D(const std::string&amp; filename) : GCTAEdisp()
94 97 * @brief Copy constructor
95 98 *
96 99 * @param[in] edisp Energy dispersion.
  100 + *
  101 + * Constructs energy dispersion by copying from another energy dispersion.
97 102 ***************************************************************************/
98 103 GCTAEdisp2D::GCTAEdisp2D(const GCTAEdisp2D& edisp) : GCTAEdisp(edisp)
99 104 {
... ... @@ -110,6 +115,8 @@ GCTAEdisp2D::GCTAEdisp2D(const GCTAEdisp2D&amp; edisp) : GCTAEdisp(edisp)
110 115  
111 116 /***********************************************************************//**
112 117 * @brief Destructor
  118 + *
  119 + * Destructs energy dispersion.
113 120 ***************************************************************************/
114 121 GCTAEdisp2D::~GCTAEdisp2D(void)
115 122 {
... ... @@ -132,6 +139,8 @@ GCTAEdisp2D::~GCTAEdisp2D(void)
132 139 *
133 140 * @param[in] edisp Energy dispersion.
134 141 * @return Energy dispersion.
  142 + *
  143 + * Assigns energy dispersion.
135 144 ***************************************************************************/
136 145 GCTAEdisp2D& GCTAEdisp2D::operator=(const GCTAEdisp2D& edisp)
137 146 {
... ... @@ -184,8 +193,14 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
184 193 // Compute Eobs/Esrc
185 194 double EobsOverEsrc = std::exp((logEobs-logEsrc) * gammalib::ln10);
186 195  
  196 + // Setup argument vector
  197 + double arg[3];
  198 + arg[m_inx_etrue] = logEsrc;
  199 + arg[m_inx_migra] = EobsOverEsrc;
  200 + arg[m_inx_theta] = theta;
  201 +
187 202 // Compute edisp
188   - edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
  203 + edisp = m_edisp(m_inx_matrix, arg[0], arg[1], arg[2]);
189 204  
190 205 // Return
191 206 return edisp;
... ... @@ -199,13 +214,13 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
199 214 ==========================================================================*/
200 215  
201 216 /***********************************************************************//**
202   - * @brief Clear instance
  217 + * @brief Clear energy dispersion
203 218 *
204   - * This method properly resets the object to an initial state.
  219 + * Clears energy dispersion.
205 220 ***************************************************************************/
206 221 void GCTAEdisp2D::clear(void)
207 222 {
208   - // Free class members (base and derived classes, derived class first)
  223 + // Free class members
209 224 free_members();
210 225 this->GCTAEdisp::free_members();
211 226  
... ... @@ -219,9 +234,11 @@ void GCTAEdisp2D::clear(void)
219 234  
220 235  
221 236 /***********************************************************************//**
222   - * @brief Clone instance
  237 + * @brief Clone energy dispersion
  238 + *
  239 + * @return Deep copy of energy dispersion.
223 240 *
224   - * @return Deep copy of energy dispersion instance.
  241 + * Returns a pointer to a deep copy of the point spread function.
225 242 ***************************************************************************/
226 243 GCTAEdisp2D* GCTAEdisp2D::clone(void) const
227 244 {
... ... @@ -230,52 +247,26 @@ GCTAEdisp2D* GCTAEdisp2D::clone(void) const
230 247  
231 248  
232 249 /***********************************************************************//**
233   - * @brief Load energy dispersion from FITS file
234   - *
235   - * @param[in] filename FITS file name.
236   - *
237   - * Loads the energy dispersion from a FITS file.
238   - *
239   - * If no extension name is provided, the energy dispersion will be loaded
240   - * from the "ENERGY DISPERSION" extension.
241   - ***************************************************************************/
242   -void GCTAEdisp2D::load(const std::string& filename)
243   -{
244   - // Create file name
245   - GFilename fname(filename);
246   -
247   - // Allocate FITS file
248   - GFits file;
249   -
250   - // Open FITS file
251   - file.open(fname.filename());
252   -
253   - // Get energy dispersion table
254   - const GFitsTable& table = *file.table(fname.extname("ENERGY DISPERSION"));
255   -
256   - // Read energy dispersion from table
257   - read(table);
258   -
259   - // Close FITS file
260   - file.close();
261   -
262   - // Store filename
263   - m_filename = filename;
264   -
265   - // Return
266   - return;
267   -}
268   -
269   -
270   -/***********************************************************************//**
271 250 * @brief Read energy dispersion from FITS table
272 251 *
273 252 * @param[in] table FITS table.
274 253 *
275   - * Reads the energy dispersion form the FITS @p table.
  254 + * @exception GException::invalid_value
  255 + * Response table is not three-dimensional.
  256 + *
  257 + * Reads the energy dispersion form the FITS @p table. The following column
  258 + * names are mandatory:
  259 + *
  260 + * ETRUE_LO - True energy lower bin boundaries
  261 + * ETRUE_HI - True energy upper bin boundaries
  262 + * MIGRA_LO - Migration lower bin boundaries
  263 + * MIGRA_HI - Migration upper bin boundaries
  264 + * THETA_LO - Offset angle lower bin boundaries
  265 + * THETA_HI - Offset angle upper bin boundaries
  266 + * MATRIX - Migration matrix
276 267 *
277   - * The data are stored in m_edisp which is of type GCTAResponseTable. The
278   - * energy axis will be set to log10, the offset angle axis to radians.
  268 + * The data are stored in the m_edisp member. The energy axis will be set to
  269 + * log10, the offset angle axis to radians.
279 270 *
280 271 * The method assures that the energy dispersion is properly normalised.
281 272 ***************************************************************************/
... ... @@ -287,16 +278,32 @@ void GCTAEdisp2D::read(const GFitsTable&amp; table)
287 278 // Read energy dispersion table
288 279 m_edisp.read(table);
289 280  
  281 + // Get mandatory indices (throw exception if not found)
  282 + m_inx_etrue = m_edisp.axis("ETRUE");
  283 + m_inx_migra = m_edisp.axis("MIGRA");
  284 + m_inx_theta = m_edisp.axis("THETA");
  285 + m_inx_matrix = m_edisp.table("MATRIX");
  286 +
  287 + // Throw an exception if the table is not two-dimensional
  288 + if (m_edisp.axes() != 3) {
  289 + std::string msg = "Expected three-dimensional energy dispersion "
  290 + "response table but found "+
  291 + gammalib::str(m_edisp.axes())+
  292 + " dimensions. Please specify a three-dimensional "
  293 + "energy dispersion.";
  294 + throw GException::invalid_value(G_READ, msg);
  295 + }
  296 +
290 297 // Set true energy axis to logarithmic scale
291   - m_edisp.axis_log10(0);
  298 + m_edisp.axis_log10(m_inx_etrue);
292 299  
293 300 // Set offset angle axis to radians
294   - m_edisp.axis_radians(2);
  301 + m_edisp.axis_radians(m_inx_theta);
295 302  
296 303 // Get axes dimensions
297   - int etrue_size = m_edisp.axis_bins(0);
298   - int migra_size = m_edisp.axis_bins(1);
299   - int theta_size = m_edisp.axis_bins(2);
  304 + int etrue_size = m_edisp.axis_bins(m_inx_etrue);
  305 + int migra_size = m_edisp.axis_bins(m_inx_migra);
  306 + int theta_size = m_edisp.axis_bins(m_inx_theta);
300 307  
301 308 // Now normalize the migration matrix vectors. We do this by integrating
302 309 // for each true energy and offset angle the measured energy. We perform
... ... @@ -313,16 +320,16 @@ void GCTAEdisp2D::read(const GFitsTable&amp; table)
313 320 for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
314 321  
315 322 // Get offset angle (in radians)
316   - double theta = 0.5 * (m_edisp.axis_lo(2,i_theta) +
317   - m_edisp.axis_hi(2,i_theta)) *
  323 + double theta = 0.5 * (m_edisp.axis_lo(m_inx_theta,i_theta) +
  324 + m_edisp.axis_hi(m_inx_theta,i_theta)) *
318 325 gammalib::deg2rad;
319 326  
320 327 // Loop over true photon energy
321 328 for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue) {
322 329  
323 330 // Get energy
324   - double emin = std::log10(m_edisp.axis_lo(0,i_etrue));
325   - double emax = std::log10(m_edisp.axis_hi(0,i_etrue));
  331 + double emin = std::log10(m_edisp.axis_lo(m_inx_etrue,i_etrue));
  332 + double emax = std::log10(m_edisp.axis_hi(m_inx_etrue,i_etrue));
326 333 double logEsrc = 0.5*(emin+emax);
327 334  
328 335 // Initialise integration
... ... @@ -364,7 +371,7 @@ void GCTAEdisp2D::read(const GFitsTable&amp; table)
364 371 if (sum > 0.0) {
365 372 int offset = i_etrue + (i_theta*migra_size) * etrue_size;
366 373 for (int k = 0, i = offset; k < migra_size; ++k, i += etrue_size) {
367   - m_edisp(0,i) /= sum;
  374 + m_edisp(m_inx_matrix,i) /= sum;
368 375 }
369 376 }
370 377 }
... ... @@ -400,6 +407,44 @@ void GCTAEdisp2D::write(GFitsBinTable&amp; table) const
400 407  
401 408  
402 409 /***********************************************************************//**
  410 + * @brief Load energy dispersion from FITS file
  411 + *
  412 + * @param[in] filename FITS file name.
  413 + *
  414 + * Loads the energy dispersion from a FITS file.
  415 + *
  416 + * If no extension name is provided, the energy dispersion will be loaded
  417 + * from the "ENERGY DISPERSION" extension.
  418 + ***************************************************************************/
  419 +void GCTAEdisp2D::load(const std::string& filename)
  420 +{
  421 + // Create file name
  422 + GFilename fname(filename);
  423 +
  424 + // Allocate FITS file
  425 + GFits file;
  426 +
  427 + // Open FITS file
  428 + file.open(fname.filename());
  429 +
  430 + // Get energy dispersion table
  431 + const GFitsTable& table = *file.table(fname.extname("ENERGY DISPERSION"));
  432 +
  433 + // Read energy dispersion from table
  434 + read(table);
  435 +
  436 + // Close FITS file
  437 + file.close();
  438 +
  439 + // Store filename
  440 + m_filename = filename;
  441 +
  442 + // Return
  443 + return;
  444 +}
  445 +
  446 +
  447 +/***********************************************************************//**
403 448 * @brief Save energy dispersion table into FITS file
404 449 *
405 450 * @param[in] filename FITS file name.
... ... @@ -522,7 +567,7 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
522 567 int high = m_ebounds_obs.size() - 1;
523 568 while ((high-low) > 1) {
524 569 int mid = (low+high) / 2;
525   - double e = m_edisp.axis_lo(0, mid);
  570 + double e = m_edisp.axis_lo(m_inx_etrue, mid);
526 571 if (logEsrc < std::log10(e)) {
527 572 high = mid;
528 573 }
... ... @@ -584,7 +629,7 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
584 629 int high = m_ebounds_src.size() - 1;
585 630 while ((high-low) > 1) {
586 631 int mid = (low+high) / 2;
587   - double e = m_edisp.axis_lo(0, mid);
  632 + double e = m_edisp.axis_lo(m_inx_etrue, mid);
588 633 if (logEobs < std::log10(e)) {
589 634 high = mid;
590 635 }
... ... @@ -618,12 +663,19 @@ std::string GCTAEdisp2D::print(const GChatter&amp; chatter) const
618 663 if (chatter != SILENT) {
619 664  
620 665 // Compute energy boundaries in TeV
621   - double emin = m_edisp.axis_lo(0,0);
622   - double emax = m_edisp.axis_hi(0,m_edisp.axis_bins(0)-1);
  666 + double emin = m_edisp.axis_lo(m_inx_etrue,0);
  667 + double emax = m_edisp.axis_hi(m_inx_etrue,
  668 + m_edisp.axis_bins(m_inx_etrue)-1);
  669 +
  670 + // Compute mingration
  671 + double mmin = m_edisp.axis_lo(m_inx_migra,0);
  672 + double mmax = m_edisp.axis_hi(m_inx_migra,
  673 + m_edisp.axis_bins(m_inx_migra)-1);
623 674  
624 675 // Compute offset angle boundaries in deg
625   - double omin = m_edisp.axis_lo(1,0);
626   - double omax = m_edisp.axis_hi(1,m_edisp.axis_bins(1)-1);
  676 + double omin = m_edisp.axis_lo(m_inx_theta,0);
  677 + double omax = m_edisp.axis_hi(m_inx_theta,
  678 + m_edisp.axis_bins(m_inx_theta)-1);
627 679  
628 680 // Append header
629 681 result.append("=== GCTAEdisp2D ===");
... ... @@ -631,11 +683,15 @@ std::string GCTAEdisp2D::print(const GChatter&amp; chatter) const
631 683 // Append information
632 684 result.append("\n"+gammalib::parformat("Filename")+m_filename);
633 685 result.append("\n"+gammalib::parformat("Number of energy bins") +
634   - gammalib::str(m_edisp.axis_bins(0)));
  686 + gammalib::str(m_edisp.axis_bins(m_inx_etrue)));
  687 + result.append("\n"+gammalib::parformat("Number of migration bins") +
  688 + gammalib::str(m_edisp.axis_bins(m_inx_migra)));
635 689 result.append("\n"+gammalib::parformat("Number of offset bins") +
636   - gammalib::str(m_edisp.axis_bins(1)));
  690 + gammalib::str(m_edisp.axis_bins(m_inx_theta)));
637 691 result.append("\n"+gammalib::parformat("Log10(Energy) range"));
638 692 result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
  693 + result.append("\n"+gammalib::parformat("Migration range"));
  694 + result.append(gammalib::str(mmin)+" - "+gammalib::str(mmax));
639 695 result.append("\n"+gammalib::parformat("Offset angle range"));
640 696 result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
641 697  
... ... @@ -660,6 +716,10 @@ void GCTAEdisp2D::init_members(void)
660 716 // Initialise members
661 717 m_filename.clear();
662 718 m_edisp.clear();
  719 + m_inx_etrue = 0;
  720 + m_inx_migra = 1;
  721 + m_inx_theta = 2;
  722 + m_inx_matrix = 0;
663 723  
664 724 // Initialise computation cache
665 725 m_ebounds_obs_computed = false;
... ... @@ -687,8 +747,12 @@ void GCTAEdisp2D::init_members(void)
687 747 void GCTAEdisp2D::copy_members(const GCTAEdisp2D& edisp)
688 748 {
689 749 // Copy members
690   - m_filename = edisp.m_filename;
691   - m_edisp = edisp.m_edisp;
  750 + m_filename = edisp.m_filename;
  751 + m_edisp = edisp.m_edisp;
  752 + m_inx_etrue = edisp.m_inx_etrue;
  753 + m_inx_migra = edisp.m_inx_migra;
  754 + m_inx_theta = edisp.m_inx_theta;
  755 + m_inx_matrix = edisp.m_inx_matrix;
692 756  
693 757 // Copy computation cache
694 758 m_ebounds_obs_computed = edisp.m_ebounds_obs_computed;
... ... @@ -743,12 +807,16 @@ void GCTAEdisp2D::compute_ebounds_obs(const double&amp; theta,
743 807 // Set epsilon
744 808 const double eps = 1.0e-12;
745 809  
  810 + // Determine number of true energy and migration bins
  811 + int n_bins = m_edisp.axis_bins(m_inx_etrue);
  812 + int n_migra = m_edisp.axis_bins(m_inx_migra);
  813 +
746 814 // Loop over Esrc
747   - for (int isrc = 0; isrc < m_edisp.axis_bins(0); ++isrc) {
  815 + for (int isrc = 0; isrc < n_bins; ++isrc) {
748 816  
749 817 // Set Esrc
750   - double Esrc = std::sqrt(m_edisp.axis_hi(0,isrc) *
751   - m_edisp.axis_lo(0,isrc));
  818 + double Esrc = std::sqrt(m_edisp.axis_hi(m_inx_etrue,isrc) *
  819 + m_edisp.axis_lo(m_inx_etrue,isrc));
752 820 double logEsrc = std::log10(Esrc);
753 821  
754 822 // Initialise results
... ... @@ -757,18 +825,19 @@ void GCTAEdisp2D::compute_ebounds_obs(const double&amp; theta,
757 825 bool minFound = false;
758 826 bool maxFound = false;
759 827  
760   - // Determine number of MIGRA bins
761   - int n_migra = m_edisp.axis_bins(1);
762   -
763 828 // Find minimum boundary
764 829 for (int i = 0; i < n_migra; ++i) {
765 830  
766 831 // Compute EobsOverEsrc value
767   - double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1,i) +
768   - m_edisp.axis_lo(1,i));
  832 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(m_inx_migra,i) +
  833 + m_edisp.axis_lo(m_inx_migra,i));
769 834  
770 835 // Get matrix term
771   - double edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
  836 + double arg[3];
  837 + arg[m_inx_etrue] = logEsrc;
  838 + arg[m_inx_migra] = EobsOverEsrc;
  839 + arg[m_inx_theta] = theta;
  840 + double edisp = m_edisp(m_inx_matrix, arg[0], arg[1], arg[2]);
772 841  
773 842 // Find first non-negligible matrix term
774 843 if (edisp >= eps) {
... ... @@ -783,11 +852,15 @@ void GCTAEdisp2D::compute_ebounds_obs(const double&amp; theta,
783 852 for (int i = n_migra-1; i >= 0; i--) {
784 853  
785 854 // Compute EobsOverEsrc value
786   - double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1,i) +
787   - m_edisp.axis_lo(1,i));
  855 + double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(m_inx_migra,i) +
  856 + m_edisp.axis_lo(m_inx_migra,i));
788 857  
789 858 // Get matrix term
790   - double edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
  859 + double arg[3];
  860 + arg[m_inx_etrue] = logEsrc;
  861 + arg[m_inx_migra] = EobsOverEsrc;
  862 + arg[m_inx_theta] = theta;
  863 + double edisp = m_edisp(m_inx_matrix, arg[0], arg[1], arg[2]);
791 864  
792 865 // Find first non-negligible matrix term
793 866 if (edisp >= eps) {
... ... @@ -846,12 +919,15 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
846 919 // Set epsilon
847 920 const double eps = 1.0e-12;
848 921  
  922 + // Determine number of true energy bins
  923 + int n_bins = m_edisp.axis_bins(m_inx_etrue);
  924 +
849 925 // Loop over Eobs
850   - for (int iobs = 0; iobs < m_edisp.axis(0); ++iobs) {
  926 + for (int iobs = 0; iobs < n_bins; ++iobs) {
851 927  
852 928 // Set Eobs
853   - double Eobs = std::sqrt(m_edisp.axis_hi(0,iobs) *
854   - m_edisp.axis_lo(0,iobs));
  929 + double Eobs = std::sqrt(m_edisp.axis_hi(m_inx_etrue,iobs) *
  930 + m_edisp.axis_lo(m_inx_etrue,iobs));
855 931 double logEobs = std::log10(Eobs);
856 932  
857 933 // Initialise results
... ... @@ -860,15 +936,12 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
860 936 bool minFound = false;
861 937 bool maxFound = false;
862 938  
863   - // Determine number of true energy bins
864   - int n_bins = m_edisp.axis(0);
865   -
866 939 // Find minimum boundary
867 940 for (int isrc = 0; isrc < n_bins; ++isrc) {
868 941  
869 942 // Set Esrc
870   - double Esrc = std::sqrt(m_edisp.axis_hi(0,isrc) *
871   - m_edisp.axis_lo(0,isrc));
  943 + double Esrc = std::sqrt(m_edisp.axis_hi(m_inx_etrue,isrc) *
  944 + m_edisp.axis_lo(m_inx_etrue,isrc));
872 945 double logEsrc = std::log10(Esrc);
873 946  
874 947 // Get matrix term
... ... @@ -887,8 +960,8 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
887 960 for (int isrc = n_bins-1; isrc >= 0; isrc--) {
888 961  
889 962 // Set Esrc
890   - double Esrc = std::sqrt(m_edisp.axis_hi(0,isrc) *
891   - m_edisp.axis_lo(0,isrc));
  963 + double Esrc = std::sqrt(m_edisp.axis_hi(m_inx_etrue,isrc) *
  964 + m_edisp.axis_lo(m_inx_etrue,isrc));
892 965 double logEsrc = std::log10(Esrc);
893 966  
894 967 // Get matrix term
... ... @@ -936,7 +1009,7 @@ void GCTAEdisp2D::set_max_edisp(void) const
936 1009  
937 1010 // Loop over all response table elements
938 1011 for (int i = 0; i < m_edisp.elements(); ++i) {
939   - double value = m_edisp(0,i);
  1012 + double value = m_edisp(m_inx_matrix,i);
940 1013 if (value > m_max_edisp) {
941 1014 m_max_edisp = value;
942 1015 }
... ...
inst/cta/src/GCTAPsf2D.cpp
  View file @e962fb3
... ... @@ -36,7 +36,6 @@
36 36 #include "GFits.hpp"
37 37 #include "GFitsBinTable.hpp"
38 38 #include "GCTAPsf2D.hpp"
39   -#include "GCTAException.hpp"
40 39  
41 40 /* __ Method name definitions ____________________________________________ */
42 41 #define G_READ "GCTAPsf2D::read(GFitsTable&)"
... ... @@ -61,6 +60,8 @@
61 60  
62 61 /***********************************************************************//**
63 62 * @brief Void constructor
  63 + *
  64 + * Constructs empty point spread function.
64 65 ***************************************************************************/
65 66 GCTAPsf2D::GCTAPsf2D(void) : GCTAPsf()
66 67 {
... ... @@ -96,6 +97,9 @@ GCTAPsf2D::GCTAPsf2D(const std::string&amp; filename) : GCTAPsf()
96 97 * @brief Copy constructor
97 98 *
98 99 * @param[in] psf Point spread function.
  100 + *
  101 + * Constructs point spread function by copying from another point spread
  102 + * function.
99 103 ***************************************************************************/
100 104 GCTAPsf2D::GCTAPsf2D(const GCTAPsf2D& psf) : GCTAPsf(psf)
101 105 {
... ... @@ -112,6 +116,8 @@ GCTAPsf2D::GCTAPsf2D(const GCTAPsf2D&amp; psf) : GCTAPsf(psf)
112 116  
113 117 /***********************************************************************//**
114 118 * @brief Destructor
  119 + *
  120 + * Destructs point spread function.
115 121 ***************************************************************************/
116 122 GCTAPsf2D::~GCTAPsf2D(void)
117 123 {
... ... @@ -134,6 +140,8 @@ GCTAPsf2D::~GCTAPsf2D(void)
134 140 *
135 141 * @param[in] psf Point spread function.
136 142 * @return Point spread function.
  143 + *
  144 + * Assigns point spread function.
137 145 ***************************************************************************/
138 146 GCTAPsf2D& GCTAPsf2D::operator=(const GCTAPsf2D& psf)
139 147 {
... ... @@ -242,13 +250,13 @@ double GCTAPsf2D::operator()(const double&amp; delta,
242 250 ==========================================================================*/
243 251  
244 252 /***********************************************************************//**
245   - * @brief Clear instance
  253 + * @brief Clear point spread function
246 254 *
247   - * This method properly resets the object to an initial state.
  255 + * Clears point spread function.
248 256 ***************************************************************************/
249 257 void GCTAPsf2D::clear(void)
250 258 {
251   - // Free class members (base and derived classes, derived class first)
  259 + // Free class members
252 260 free_members();
253 261 this->GCTAPsf::free_members();
254 262  
... ... @@ -262,9 +270,11 @@ void GCTAPsf2D::clear(void)
262 270  
263 271  
264 272 /***********************************************************************//**
265   - * @brief Clone instance
  273 + * @brief Clone point spread functions
  274 + *
  275 + * @return Deep copy of point spread function.
266 276 *
267   - * @return Deep copy of point spread function instance.
  277 + * Returns a pointer to a deep copy of the point spread function.
268 278 ***************************************************************************/
269 279 GCTAPsf2D* GCTAPsf2D::clone(void) const
270 280 {
... ... @@ -278,13 +288,23 @@ GCTAPsf2D* GCTAPsf2D::clone(void) const
278 288 * @param[in] table FITS table.
279 289 *
280 290 * @exception GException::invalid_value
281   - * FITS file format differs from expectation.
282   - *
283   - * Reads point spread function from the FITS @p table.
284   - *
285   - * The data are stored in m_psf which is of type GCTAResponseTable. The
286   - * energy axis will be set to log10. The offset angle axis and sigma
287   - * parameter columns will be set to radians.
  291 + * Response table is not two-dimensional.
  292 + *
  293 + * Reads the point spread function form the FITS @p table. The following
  294 + * column names are mandatory:
  295 + *
  296 + * ENERG_LO - Energy lower bin boundaries
  297 + * ENERG_HI - Energy upper bin boundaries
  298 + * THETA_LO - Offset angle lower bin boundaries
  299 + * THETA_HI - Offset angle upper bin boundaries
  300 + * SIGMA_1 - 1st Gaussian sigma
  301 + * AMPL_2 - 2nd Gaussian relative amplitude
  302 + * SIGMA_2 - 2nd Gaussian sigma
  303 + * AMPL_3 - 3rd Gaussian relative amplitude
  304 + * SIGMA_3 - 3rd Gaussian sigma
  305 + *
  306 + * The data are stored in the m_psf member. The energy axis will be set to
  307 + * log10, the offset angle axis to radians.
288 308 ***************************************************************************/
289 309 void GCTAPsf2D::read(const GFitsTable& table)
290 310 {
... ... @@ -294,32 +314,35 @@ void GCTAPsf2D::read(const GFitsTable&amp; table)
294 314 // Read PSF table
295 315 m_psf.read(table);
296 316  
297   - // Check that axis names comply to format
298   - if (m_psf.axis_lo_name(0) != "ENERG_LO" ||
299   - m_psf.axis_hi_name(0) != "ENERG_HI") {
300   - std::string msg = "Point spread function response table does not"
301   - " contain \"ENERG_LO\" and \"ENERG_HI\" columns"
302   - " as the first axis.";
303   - throw GException::invalid_value(G_READ, msg);
304   - }
305   - if (m_psf.axis_lo_name(1) != "THETA_LO" ||
306   - m_psf.axis_hi_name(1) != "THETA_HI") {
307   - std::string msg = "Point spread function response table does not"
308   - " contain \"THETA_LO\" and \"THETA_HI\" columns"
309   - " as the second axis.";
  317 + // Get mandatory indices (throw exception if not found)
  318 + m_inx_energy = m_psf.axis("ENERG");
  319 + m_inx_theta = m_psf.axis("THETA");
  320 + m_inx_sigma1 = m_psf.table("SIGMA_1");
  321 + m_inx_ampl2 = m_psf.table("AMPL_2");
  322 + m_inx_sigma2 = m_psf.table("SIGMA_2");
  323 + m_inx_ampl3 = m_psf.table("AMPL_3");
  324 + m_inx_sigma3 = m_psf.table("SIGMA_3");
  325 +
  326 + // Throw an exception if the table is not two-dimensional
  327 + if (m_psf.axes() != 2) {
  328 + std::string msg = "Expected two-dimensional point spread function "
  329 + "response table but found "+
  330 + gammalib::str(m_psf.axes())+
  331 + " dimensions. Please specify a two-dimensional "
  332 + "point spread function.";
310 333 throw GException::invalid_value(G_READ, msg);
311 334 }
312 335  
313 336 // Set energy axis to logarithmic scale
314   - m_psf.axis_log10(0);
  337 + m_psf.axis_log10(m_inx_energy);
315 338  
316 339 // Set offset angle axis to radians
317   - m_psf.axis_radians(1);
  340 + m_psf.axis_radians(m_inx_theta);
318 341  
319 342 // Convert sigma parameters to radians
320   - m_psf.scale(1, gammalib::deg2rad);
321   - m_psf.scale(3, gammalib::deg2rad);
322   - m_psf.scale(5, gammalib::deg2rad);
  343 + m_psf.scale(m_inx_sigma1, gammalib::deg2rad);
  344 + m_psf.scale(m_inx_sigma2, gammalib::deg2rad);
  345 + m_psf.scale(m_inx_sigma3, gammalib::deg2rad);
323 346  
324 347 // Return
325 348 return;
... ... @@ -341,11 +364,11 @@ void GCTAPsf2D::write(GFitsBinTable&amp; table) const
341 364 GCTAResponseTable psf(m_psf);
342 365  
343 366 // Convert sigma parameters back to degrees
344   - psf.scale(1, gammalib::rad2deg);
345   - psf.scale(3, gammalib::rad2deg);
346   - psf.scale(5, gammalib::rad2deg);
  367 + psf.scale(m_inx_sigma1, gammalib::rad2deg);
  368 + psf.scale(m_inx_sigma2, gammalib::rad2deg);
  369 + psf.scale(m_inx_sigma3, gammalib::rad2deg);
347 370  
348   - // Write background table
  371 + // Write response table
349 372 psf.write(table);
350 373  
351 374 // Return
... ... @@ -626,12 +649,14 @@ std::string GCTAPsf2D::print(const GChatter&amp; chatter) const
626 649 if (chatter != SILENT) {
627 650  
628 651 // Compute energy boundaries in TeV
629   - double emin = m_psf.axis_lo(0,0);
630   - double emax = m_psf.axis_hi(0,m_psf.axis_bins(0)-1);
  652 + double emin = m_psf.axis_lo(m_inx_energy,0);
  653 + double emax = m_psf.axis_hi(m_inx_energy,
  654 + m_psf.axis_bins(m_inx_energy)-1);
631 655  
632 656 // Compute offset angle boundaries in deg
633   - double omin = m_psf.axis_lo(1,0);
634   - double omax = m_psf.axis_hi(1,m_psf.axis_bins(1)-1);
  657 + double omin = m_psf.axis_lo(m_inx_theta,0);
  658 + double omax = m_psf.axis_hi(m_inx_theta,
  659 + m_psf.axis_bins(m_inx_theta)-1);
635 660  
636 661 // Append header
637 662 result.append("=== GCTAPsf2D ===");
... ... @@ -639,9 +664,9 @@ std::string GCTAPsf2D::print(const GChatter&amp; chatter) const
639 664 // Append information
640 665 result.append("\n"+gammalib::parformat("Filename")+m_filename);
641 666 result.append("\n"+gammalib::parformat("Number of energy bins") +
642   - gammalib::str(m_psf.axis_bins(0)));
  667 + gammalib::str(m_psf.axis_bins(m_inx_energy)));
643 668 result.append("\n"+gammalib::parformat("Number of offset bins") +
644   - gammalib::str(m_psf.axis_bins(1)));
  669 + gammalib::str(m_psf.axis_bins(m_inx_theta)));
645 670 result.append("\n"+gammalib::parformat("Log10(Energy) range"));
646 671 result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
647 672 result.append("\n"+gammalib::parformat("Offset angle range"));
... ... @@ -668,17 +693,24 @@ void GCTAPsf2D::init_members(void)
668 693 // Initialise members
669 694 m_filename.clear();
670 695 m_psf.clear();
671   - m_par_logE = -1.0e30;
672   - m_par_theta = -1.0;
673   - m_norm = 1.0;
674   - m_norm2 = 0.0;
675   - m_norm3 = 0.0;
676   - m_sigma1 = 0.0;
677   - m_sigma2 = 0.0;
678   - m_sigma3 = 0.0;
679   - m_width1 = 0.0;
680   - m_width2 = 0.0;
681   - m_width3 = 0.0;
  696 + m_inx_energy = 0;
  697 + m_inx_theta = 1;
  698 + m_inx_sigma1 = 0;
  699 + m_inx_ampl2 = 1;
  700 + m_inx_sigma2 = 2;
  701 + m_inx_ampl3 = 3;
  702 + m_inx_sigma3 = 4;
  703 + m_par_logE = -1.0e30;
  704 + m_par_theta = -1.0;
  705 + m_norm = 1.0;
  706 + m_norm2 = 0.0;
  707 + m_norm3 = 0.0;
  708 + m_sigma1 = 0.0;
  709 + m_sigma2 = 0.0;
  710 + m_sigma3 = 0.0;
  711 + m_width1 = 0.0;
  712 + m_width2 = 0.0;
  713 + m_width3 = 0.0;
682 714  
683 715 // Return
684 716 return;
... ... @@ -693,19 +725,26 @@ void GCTAPsf2D::init_members(void)
693 725 void GCTAPsf2D::copy_members(const GCTAPsf2D& psf)
694 726 {
695 727 // Copy members
696   - m_filename = psf.m_filename;
697   - m_psf = psf.m_psf;
698   - m_par_logE = psf.m_par_logE;
699   - m_par_theta = psf.m_par_theta;
700   - m_norm = psf.m_norm;
701   - m_norm2 = psf.m_norm2;
702   - m_norm3 = psf.m_norm3;
703   - m_sigma1 = psf.m_sigma1;
704   - m_sigma2 = psf.m_sigma2;
705   - m_sigma3 = psf.m_sigma3;
706   - m_width1 = psf.m_width1;
707   - m_width2 = psf.m_width2;
708   - m_width3 = psf.m_width3;
  728 + m_filename = psf.m_filename;
  729 + m_psf = psf.m_psf;
  730 + m_inx_energy = psf.m_inx_energy;
  731 + m_inx_theta = psf.m_inx_theta;
  732 + m_inx_sigma1 = psf.m_inx_sigma1;
  733 + m_inx_ampl2 = psf.m_inx_ampl2;
  734 + m_inx_sigma2 = psf.m_inx_sigma2;
  735 + m_inx_ampl3 = psf.m_inx_ampl3;
  736 + m_inx_sigma3 = psf.m_inx_sigma3;
  737 + m_par_logE = psf.m_par_logE;
  738 + m_par_theta = psf.m_par_theta;
  739 + m_norm = psf.m_norm;
  740 + m_norm2 = psf.m_norm2;
  741 + m_norm3 = psf.m_norm3;
  742 + m_sigma1 = psf.m_sigma1;
  743 + m_sigma2 = psf.m_sigma2;
  744 + m_sigma3 = psf.m_sigma3;
  745 + m_width1 = psf.m_width1;
  746 + m_width2 = psf.m_width2;
  747 + m_width3 = psf.m_width3;
709 748  
710 749 // Return
711 750 return;
... ... @@ -740,12 +779,13 @@ void GCTAPsf2D::update(const double&amp; logE, const double&amp; theta) const
740 779 m_par_theta = theta;
741 780  
742 781 // Interpolate response parameters
743   - std::vector<double> pars = m_psf(logE, theta);
  782 + std::vector<double> pars = (m_inx_energy == 0) ? m_psf(logE, theta) :
  783 + m_psf(theta, logE);
744 784  
745 785 // Set Gaussian sigmas
746   - m_sigma1 = pars[1];
747   - m_sigma2 = pars[3];
748   - m_sigma3 = pars[5];
  786 + m_sigma1 = pars[m_inx_sigma1];
  787 + m_sigma2 = pars[m_inx_sigma2];
  788 + m_sigma3 = pars[m_inx_sigma3];
749 789  
750 790 // Set width parameters
751 791 double sigma1 = m_sigma1 * m_sigma1;
... ... @@ -763,7 +803,7 @@ void GCTAPsf2D::update(const double&amp; logE, const double&amp; theta) const
763 803 // Compute Gaussian 2
764 804 if (sigma2 > 0.0) {
765 805 m_width2 = -0.5 / sigma2;
766   - m_norm2 = pars[2];
  806 + m_norm2 = pars[m_inx_ampl2];
767 807 }
768 808 else {
769 809 m_width2 = 0.0;
... ... @@ -773,7 +813,7 @@ void GCTAPsf2D::update(const double&amp; logE, const double&amp; theta) const
773 813 // Compute Gaussian 3
774 814 if (sigma3 > 0.0) {
775 815 m_width3 = -0.5 / sigma3;
776   - m_norm3 = pars[4];
  816 + m_norm3 = pars[m_inx_ampl3];
777 817 }
778 818 else {
779 819 m_width3 = 0.0;
... ... @@ -781,7 +821,8 @@ void GCTAPsf2D::update(const double&amp; logE, const double&amp; theta) const
781 821 }
782 822  
783 823 // Compute global normalization parameter
784   - double integral = gammalib::twopi * (sigma1 + sigma2*m_norm2 + sigma3*m_norm3);
  824 + double integral = gammalib::twopi *
  825 + (sigma1 + sigma2*m_norm2 + sigma3*m_norm3);
785 826 m_norm = (integral > 0.0) ? 1.0 / integral : 0.0;
786 827  
787 828 }
... ...
inst/cta/src/GCTAPsfKing.cpp
  View file @e962fb3
1 1 /***************************************************************************
2 2 * GCTAPsfKing.cpp - King profile CTA point spread function class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2013-2015 by Michael Mayer *
  4 + * copyright (C) 2013-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -65,6 +65,8 @@ const double r_max = 0.7 * gammalib::deg2rad; // Maximum delta for fixed
65 65  
66 66 /***********************************************************************//**
67 67 * @brief Void constructor
  68 + *
  69 + * Constructs empty point spread function.
68 70 ***************************************************************************/
69 71 GCTAPsfKing::GCTAPsfKing(void) : GCTAPsf()
70 72 {
... ... @@ -100,6 +102,9 @@ GCTAPsfKing::GCTAPsfKing(const std::string&amp; filename) : GCTAPsf()
100 102 * @brief Copy constructor
101 103 *
102 104 * @param[in] psf Point spread function.
  105 + *
  106 + * Constructs point spread function by copying from another point spread
  107 + * function.
103 108 ***************************************************************************/
104 109 GCTAPsfKing::GCTAPsfKing(const GCTAPsfKing& psf) : GCTAPsf(psf)
105 110 {
... ... @@ -116,6 +121,8 @@ GCTAPsfKing::GCTAPsfKing(const GCTAPsfKing&amp; psf) : GCTAPsf(psf)
116 121  
117 122 /***********************************************************************//**
118 123 * @brief Destructor
  124 + *
  125 + * Destructs point spread function.
119 126 ***************************************************************************/
120 127 GCTAPsfKing::~GCTAPsfKing(void)
121 128 {
... ... @@ -138,6 +145,8 @@ GCTAPsfKing::~GCTAPsfKing(void)
138 145 *
139 146 * @param[in] psf Point spread function.
140 147 * @return Point spread function.
  148 + *
  149 + * Assigns point spread function.
141 150 ***************************************************************************/
142 151 GCTAPsfKing& GCTAPsfKing::operator=(const GCTAPsfKing& psf)
143 152 {
... ... @@ -244,9 +253,9 @@ double GCTAPsfKing::operator()(const double&amp; delta,
244 253 ==========================================================================*/
245 254  
246 255 /***********************************************************************//**
247   - * @brief Clear instance
  256 + * @brief Clear point spread function
248 257 *
249   - * This method properly resets the object to an initial state.
  258 + * Clears point spread function.
250 259 ***************************************************************************/
251 260 void GCTAPsfKing::clear(void)
252 261 {
... ... @@ -264,9 +273,11 @@ void GCTAPsfKing::clear(void)
264 273  
265 274  
266 275 /***********************************************************************//**
267   - * @brief Clone instance
  276 + * @brief Clone point spread functions
268 277 *
269   - * @return Deep copy of point spread function instance.
  278 + * @return Deep copy of point spread function.
  279 + *
  280 + * Returns a pointer to a deep copy of the point spread function.
270 281 ***************************************************************************/
271 282 GCTAPsfKing* GCTAPsfKing::clone(void) const
272 283 {
... ... @@ -280,13 +291,20 @@ GCTAPsfKing* GCTAPsfKing::clone(void) const
280 291 * @param[in] table FITS table.
281 292 *
282 293 * @exception GException::invalid_value
283   - * FITS file format differs from expectation.
  294 + * Response table is not two-dimensional.
  295 + *
  296 + * Reads the point spread function form the FITS @p table. The following
  297 + * column names are mandatory:
284 298 *
285   - * Reads the point spread function from the FITS @p table.
  299 + * ENERG_LO - Energy lower bin boundaries
  300 + * ENERG_HI - Energy upper bin boundaries
  301 + * THETA_LO - Offset angle lower bin boundaries
  302 + * THETA_HI - Offset angle upper bin boundaries
  303 + * GAMMA - Gamma
  304 + * SIGMA - Sigma
286 305 *
287   - * The data are stored in m_psf which is of type GCTAResponseTable. The
288   - * energy axis will be set to log10. The offset angle axis and sigma
289   - * parameter columns will be set to radians.
  306 + * The data are stored in the m_psf member. The energy axis will be set to
  307 + * log10, the offset angle axis to radians.
290 308 ***************************************************************************/
291 309 void GCTAPsfKing::read(const GFitsTable& table)
292 310 {
... ... @@ -296,30 +314,30 @@ void GCTAPsfKing::read(const GFitsTable&amp; table)
296 314 // Read PSF table
297 315 m_psf.read(table);
298 316  
299   - // Check that axis names comply to format
300   - if (m_psf.axis_lo_name(0) != "ENERG_LO" ||
301   - m_psf.axis_hi_name(0) != "ENERG_HI") {
302   - std::string msg = "Point spread function response table does not"
303   - " contain \"ENERG_LO\" and \"ENERG_HI\" columns"
304   - " as the first axis.";
305   - throw GException::invalid_value(G_READ, msg);
306   - }
307   - if (m_psf.axis_lo_name(1) != "THETA_LO" ||
308   - m_psf.axis_hi_name(1) != "THETA_HI") {
309   - std::string msg = "Point spread function response table does not"
310   - " contain \"THETA_LO\" and \"THETA_HI\" columns"
311   - " as the second axis.";
  317 + // Get mandatory indices (throw exception if not found)
  318 + m_inx_energy = m_psf.axis("ENERG");
  319 + m_inx_theta = m_psf.axis("THETA");
  320 + m_inx_gamma = m_psf.table("GAMMA");
  321 + m_inx_sigma = m_psf.table("SIGMA");
  322 +
  323 + // Throw an exception if the table is not two-dimensional
  324 + if (m_psf.axes() != 2) {
  325 + std::string msg = "Expected two-dimensional point spread function "
  326 + "response table but found "+
  327 + gammalib::str(m_psf.axes())+
  328 + " dimensions. Please specify a two-dimensional "
  329 + "point spread function.";
312 330 throw GException::invalid_value(G_READ, msg);
313 331 }
314 332  
315 333 // Set energy axis to logarithmic scale
316   - m_psf.axis_log10(0);
  334 + m_psf.axis_log10(m_inx_energy);
317 335  
318 336 // Set offset angle axis to radians
319   - m_psf.axis_radians(1);
  337 + m_psf.axis_radians(m_inx_theta);
320 338  
321 339 // Convert sigma parameters to radians
322   - m_psf.scale(1, gammalib::deg2rad);
  340 + m_psf.scale(m_inx_sigma, gammalib::deg2rad);
323 341  
324 342 // Return
325 343 return;
... ... @@ -337,8 +355,14 @@ void GCTAPsfKing::read(const GFitsTable&amp; table)
337 355 ***************************************************************************/
338 356 void GCTAPsfKing::write(GFitsBinTable& table) const
339 357 {
340   - // Write background table
341   - m_psf.write(table);
  358 + // Create a copy of the response table
  359 + GCTAResponseTable psf(m_psf);
  360 +
  361 + // Convert sigma parameters back to degrees
  362 + psf.scale(m_inx_sigma, gammalib::rad2deg);
  363 +
  364 + // Write response table
  365 + psf.write(table);
342 366  
343 367 // Return
344 368 return;
... ... @@ -558,19 +582,41 @@ double GCTAPsfKing::containment_radius(const double&amp; fraction,
558 582 * @brief Print point spread function information
559 583 *
560 584 * @return Content of point spread function instance.
  585 + * @return String containing point spread function information.
561 586 ***************************************************************************/
562 587 std::string GCTAPsfKing::print(const GChatter& chatter) const
563 588 {
564 589 // Initialise result string
565 590 std::string result;
566 591  
  592 + // Continue only if chatter is not silent
567 593 if (chatter != SILENT) {
568 594  
569   - // Append information
570   - result.append("=== GCTAPsfKing ===");
571   - result.append("\n"+gammalib::parformat("Filename")+m_filename);
572   - result.append("\n"+m_psf.print(chatter));
573   - }
  595 + // Compute energy boundaries in TeV
  596 + double emin = m_psf.axis_lo(m_inx_energy,0);
  597 + double emax = m_psf.axis_hi(m_inx_energy,
  598 + m_psf.axis_bins(m_inx_energy)-1);
  599 +
  600 + // Compute offset angle boundaries in deg
  601 + double omin = m_psf.axis_lo(m_inx_theta,0);
  602 + double omax = m_psf.axis_hi(m_inx_theta,
  603 + m_psf.axis_bins(m_inx_theta)-1);
  604 +
  605 + // Append header
  606 + result.append("=== GCTAPsfKing ===");
  607 +
  608 + // Append information
  609 + result.append("\n"+gammalib::parformat("Filename")+m_filename);
  610 + result.append("\n"+gammalib::parformat("Number of energy bins") +
  611 + gammalib::str(m_psf.axis_bins(m_inx_energy)));
  612 + result.append("\n"+gammalib::parformat("Number of offset bins") +
  613 + gammalib::str(m_psf.axis_bins(m_inx_theta)));
  614 + result.append("\n"+gammalib::parformat("Log10(Energy) range"));
  615 + result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
  616 + result.append("\n"+gammalib::parformat("Offset angle range"));
  617 + result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
  618 +
  619 + } // endif: chatter was not silent
574 620  
575 621 // Return result
576 622 return result;
... ... @@ -591,6 +637,10 @@ void GCTAPsfKing::init_members(void)
591 637 // Initialise members
592 638 m_filename.clear();
593 639 m_psf.clear();
  640 + m_inx_energy = 0;
  641 + m_inx_theta = 1;
  642 + m_inx_gamma = 0;
  643 + m_inx_sigma = 1;
594 644 m_par_logE = -1.0e30;
595 645 m_par_theta = -1.0;
596 646 m_par_norm = 0.0;
... ... @@ -613,6 +663,10 @@ void GCTAPsfKing::copy_members(const GCTAPsfKing&amp; psf)
613 663 // Copy members
614 664 m_filename = psf.m_filename;
615 665 m_psf = psf.m_psf;
  666 + m_inx_energy = psf.m_inx_energy;
  667 + m_inx_theta = psf.m_inx_theta;
  668 + m_inx_gamma = psf.m_inx_gamma;
  669 + m_inx_sigma = psf.m_inx_sigma;
616 670 m_par_logE = psf.m_par_logE;
617 671 m_par_theta = psf.m_par_theta;
618 672 m_par_norm = psf.m_par_norm;
... ... @@ -658,23 +712,12 @@ void GCTAPsfKing::update(const double&amp; logE, const double&amp; theta) const
658 712 m_par_theta = theta;
659 713  
660 714 // Determine sigma and gamma by interpolating between nodes
661   - std::vector<double> pars = m_psf(logE,theta);
662   -
663   - // Throw an exception if there are not 2 parameters
664   - if (pars.size() != 2) {
665   - std::string msg = gammalib::str(pars.size()) + " parameters have"
666   - " been found in the response table of the"
667   - " King profile response function while 2"
668   - " parameters are expected.\n"
669   - "Possibly, the point spread function information"
670   - " has not yet been loaded. Please load the point"
671   - " spread function before using it.";
672   - throw GException::invalid_value(G_UPDATE, msg);
673   - }
  715 + std::vector<double> pars = (m_inx_energy == 0) ? m_psf(logE, theta) :
  716 + m_psf(theta, logE);
674 717  
675 718 // Set parameters
676   - m_par_gamma = pars[0];
677   - m_par_sigma = pars[1];
  719 + m_par_gamma = pars[m_inx_gamma];
  720 + m_par_sigma = pars[m_inx_sigma];
678 721 m_par_sigma2 = m_par_sigma * m_par_sigma;
679 722  
680 723 // Check for parameter sanity
... ...