GammaLib

I changed GCTAEdisp2D so that the matrix is now normalized like the H.E.S.S. matrix, i.e. the energy dispersion integrated over the migration values is unity. I also changed the GCTAEdisp2D::operator() so that it now returns a density per MeV, as is needed later for the response computation. The GCTAEdisp2D::prob_erecobin() was also adapted to that it returns the correct probabilities for the new definition. I tested the new code using the test_edisp.py script.

I then adapted GCTAResponseIrf::edisp() which no longer needs a transformation. The method returns a energy dispersion density per MeV, which is now directly provided by the GCTAEdisp2D::operator().

I added some debug code to GCTAResponseIrf, switched on using the G_DEBUG_NROI_EDISP definition, that compares the computation with and without energy dispersion. The code is used in GCTAResponseIrf::nroi which compares the nroi computations. I first switched off the weighting with the spectra, hence in principle the values without and with energy dispersion should be identical. Here is what I get for the first observation of RX J1713:

GCTAResponseIrf::nroi: obs_id=020326 energy=400 GeV nroi(noedisp)=1.16801e+09 nroi(edisp)=1.19095e+09 ratio(edisp/noedisp)=1.01964
GCTAResponseIrf::nroi: obs_id=020326 energy=50 TeV nroi(noedisp)=1.96924e+09 nroi(edisp)=2.13709e+09 ratio(edisp/noedisp)=1.08523
GCTAResponseIrf::nroi: obs_id=020326 energy=4.47213595499958 TeV nroi(noedisp)=2.48134e+09 nroi(edisp)=2.51083e+09 ratio(edisp/noedisp)=1.01188
GCTAResponseIrf::nroi: obs_id=020326 energy=1.33748060995284 TeV nroi(noedisp)=2.04281e+09 nroi(edisp)=2.05556e+09 ratio(edisp/noedisp)=1.00624
GCTAResponseIrf::nroi: obs_id=020326 energy=14.9534878122122 TeV nroi(noedisp)=2.41627e+09 nroi(edisp)=2.56464e+09 ratio(edisp/noedisp)=1.06141
GCTAResponseIrf::nroi: obs_id=020326 energy=731.43163999183 GeV nroi(noedisp)=1.61344e+09 nroi(edisp)=1.62867e+09 ratio(edisp/noedisp)=1.00944
GCTAResponseIrf::nroi: obs_id=020326 energy=2.4456890899877 TeV nroi(noedisp)=2.34658e+09 nroi(edisp)=2.37909e+09 ratio(edisp/noedisp)=1.01385
GCTAResponseIrf::nroi: obs_id=020326 energy=8.17765433957943 TeV nroi(noedisp)=2.49022e+09 nroi(edisp)=2.48173e+09 ratio(edisp/noedisp)=0.996591
GCTAResponseIrf::nroi: obs_id=020326 energy=27.3436352852105 TeV nroi(noedisp)=2.25699e+09 nroi(edisp)=2.27037e+09 ratio(edisp/noedisp)=1.00593
GCTAResponseIrf::nroi: obs_id=020326 energy=540.899857641627 GeV nroi(noedisp)=1.38996e+09 nroi(edisp)=1.37825e+09 ratio(edisp/noedisp)=0.991575
GCTAResponseIrf::nroi: obs_id=020326 energy=989.078174865405 GeV nroi(noedisp)=1.83682e+09 nroi(edisp)=1.81106e+09 ratio(edisp/noedisp)=0.985975
GCTAResponseIrf::nroi: obs_id=020326 energy=1.80860767880482 TeV nroi(noedisp)=2.21561e+09 nroi(edisp)=2.20908e+09 ratio(edisp/noedisp)=0.997054
GCTAResponseIrf::nroi: obs_id=020326 energy=3.30718220152506 TeV nroi(noedisp)=2.43407e+09 nroi(edisp)=2.4275e+09 ratio(edisp/noedisp)=0.997303
GCTAResponseIrf::nroi: obs_id=020326 energy=6.04744425353316 TeV nroi(noedisp)=2.49812e+09 nroi(edisp)=2.52044e+09 ratio(edisp/noedisp)=1.00893
GCTAResponseIrf::nroi: obs_id=020326 energy=11.0582301703023 TeV nroi(noedisp)=2.46243e+09 nroi(edisp)=2.45639e+09 ratio(edisp/noedisp)=0.99755
GCTAResponseIrf::nroi: obs_id=020326 energy=20.2208485721784 TeV nroi(noedisp)=2.34885e+09 nroi(edisp)=2.41567e+09 ratio(edisp/noedisp)=1.02845
GCTAResponseIrf::nroi: obs_id=020326 energy=36.9754210829371 TeV nroi(noedisp)=2.13018e+09 nroi(edisp)=2.26431e+09 ratio(edisp/noedisp)=1.06297
GCTAResponseIrf::nroi: obs_id=020326 energy=465.145077429237 GeV nroi(noedisp)=1.27497e+09 nroi(edisp)=1.27984e+09 ratio(edisp/noedisp)=1.00382
GCTAResponseIrf::nroi: obs_id=020326 energy=628.992265410444 GeV nroi(noedisp)=1.50057e+09 nroi(edisp)=1.5158e+09 ratio(edisp/noedisp)=1.01015
GCTAResponseIrf::nroi: obs_id=020326 energy=850.554567045483 GeV nroi(noedisp)=1.72893e+09 nroi(edisp)=1.69352e+09 ratio(edisp/noedisp)=0.979515
GCTAResponseIrf::nroi: obs_id=020326 energy=1.15016211057834 TeV nroi(noedisp)=1.94609e+09 nroi(edisp)=1.93318e+09 ratio(edisp/noedisp)=0.993365
GCTAResponseIrf::nroi: obs_id=020326 energy=1.55530630469155 TeV nroi(noedisp)=2.13237e+09 nroi(edisp)=2.16003e+09 ratio(edisp/noedisp)=1.01297
GCTAResponseIrf::nroi: obs_id=020326 energy=2.10316239699195 TeV nroi(noedisp)=2.28559e+09 nroi(edisp)=2.30368e+09 ratio(edisp/noedisp)=1.00792
GCTAResponseIrf::nroi: obs_id=020326 energy=2.84400060282542 TeV nroi(noedisp)=2.39412e+09 nroi(edisp)=2.42404e+09 ratio(edisp/noedisp)=1.0125
GCTAResponseIrf::nroi: obs_id=020326 energy=3.84579880300243 TeV nroi(noedisp)=2.46163e+09 nroi(edisp)=2.42719e+09 ratio(edisp/noedisp)=0.986013
GCTAResponseIrf::nroi: obs_id=020326 energy=5.20048006265587 TeV nroi(noedisp)=2.49336e+09 nroi(edisp)=2.54699e+09 ratio(edisp/noedisp)=1.02151
GCTAResponseIrf::nroi: obs_id=020326 energy=7.03234731389669 TeV nroi(noedisp)=2.49685e+09 nroi(edisp)=2.54483e+09 ratio(edisp/noedisp)=1.01922
GCTAResponseIrf::nroi: obs_id=020326 energy=9.50948915243301 TeV nroi(noedisp)=2.47846e+09 nroi(edisp)=2.42323e+09 ratio(edisp/noedisp)=0.977717
GCTAResponseIrf::nroi: obs_id=020326 energy=12.859203321989 TeV nroi(noedisp)=2.44232e+09 nroi(edisp)=2.52239e+09 ratio(edisp/noedisp)=1.03278
GCTAResponseIrf::nroi: obs_id=020326 energy=17.3888531156215 TeV nroi(noedisp)=2.3863e+09 nroi(edisp)=2.42444e+09 ratio(edisp/noedisp)=1.01598
GCTAResponseIrf::nroi: obs_id=020326 energy=23.5140704369771 TeV nroi(noedisp)=2.30605e+09 nroi(edisp)=2.41935e+09 ratio(edisp/noedisp)=1.04913
GCTAResponseIrf::nroi: obs_id=020326 energy=31.7968933798403 TeV nroi(noedisp)=2.19582e+09 nroi(edisp)=2.34322e+09 ratio(edisp/noedisp)=1.06712
GCTAResponseIrf::nroi: obs_id=020326 energy=42.9973377565041 TeV nroi(noedisp)=2.05282e+09 nroi(edisp)=2.31166e+09 ratio(edisp/noedisp)=1.12609
GCTAResponseIrf::nroi: obs_id=020326 energy=431.344445857015 GeV nroi(noedisp)=1.2213e+09 nroi(edisp)=1.23612e+09 ratio(edisp/noedisp)=1.01213
GCTAResponseIrf::nroi: obs_id=020326 energy=501.594364167081 GeV nroi(noedisp)=1.33102e+09 nroi(edisp)=1.33424e+09 ratio(edisp/noedisp)=1.00242
GCTAResponseIrf::nroi: obs_id=020326 energy=583.285373396415 GeV nroi(noedisp)=1.44509e+09 nroi(edisp)=1.46001e+09 ratio(edisp/noedisp)=1.01033

And here the values with the full nroi computation. Here a difference is expected due to the weighting with the energy spectrum. I was surprised that the difference is so large. The values below are ordered by increasing reconstructed energy, which shows that the difference is larger for lower energies, where the energy dispersion is larger, and smaller for higher energies, where the energy dispersion is smaller.

GCTAResponseIrf::nroi: obs_id=020326 energy=400 GeV              nroi(noedisp)=1.1139e-12  nroi(edisp)=9.26693e-11 ratio(edisp/noedisp)=83.1936
GCTAResponseIrf::nroi: obs_id=020326 energy=431.344445857015 GeV nroi(noedisp)=5.47747e-13 nroi(edisp)=5.22089e-11 ratio(edisp/noedisp)=95.3157
GCTAResponseIrf::nroi: obs_id=020326 energy=465.145077429237 GeV nroi(noedisp)=2.68915e-13 nroi(edisp)=5.55211e-11 ratio(edisp/noedisp)=206.463
GCTAResponseIrf::nroi: obs_id=020326 energy=501.594364167081 GeV nroi(noedisp)=1.32025e-13 nroi(edisp)=2.15619e-11 ratio(edisp/noedisp)=163.317
GCTAResponseIrf::nroi: obs_id=020326 energy=540.899857641627 GeV nroi(noedisp)=6.48387e-14 nroi(edisp)=1.81986e-11 ratio(edisp/noedisp)=280.675
GCTAResponseIrf::nroi: obs_id=020326 energy=628.992265410444 GeV nroi(noedisp)=1.54811e-14 nroi(edisp)=8.01247e-12 ratio(edisp/noedisp)=517.564
GCTAResponseIrf::nroi: obs_id=020326 energy=731.43163999183 GeV  nroi(noedisp)=3.68143e-15 nroi(edisp)=2.19206e-12 ratio(edisp/noedisp)=595.439
GCTAResponseIrf::nroi: obs_id=020326 energy=850.554567045483 GeV nroi(noedisp)=8.72482e-16 nroi(edisp)=3.46074e-13 ratio(edisp/noedisp)=396.654
GCTAResponseIrf::nroi: obs_id=020326 energy=989.078174865405 GeV nroi(noedisp)=2.05003e-16 nroi(edisp)=7.03383e-14 ratio(edisp/noedisp)=343.108
GCTAResponseIrf::nroi: obs_id=020326 energy=1.15016211057834 TeV nroi(noedisp)=4.80367e-17 nroi(edisp)=1.73202e-14 ratio(edisp/noedisp)=360.562
GCTAResponseIrf::nroi: obs_id=020326 energy=1.33748060995284 TeV nroi(noedisp)=1.1152e-17  nroi(edisp)=3.54512e-15 ratio(edisp/noedisp)=317.89
GCTAResponseIrf::nroi: obs_id=020326 energy=1.55530630469155 TeV nroi(noedisp)=2.57456e-18 nroi(edisp)=1.14239e-15 ratio(edisp/noedisp)=443.723
GCTAResponseIrf::nroi: obs_id=020326 energy=1.80860767880482 TeV nroi(noedisp)=5.9163e-19  nroi(edisp)=5.05359e-16 ratio(edisp/noedisp)=854.181
GCTAResponseIrf::nroi: obs_id=020326 energy=2.10316239699195 TeV nroi(noedisp)=1.3498e-19  nroi(edisp)=3.93526e-17 ratio(edisp/noedisp)=291.543
GCTAResponseIrf::nroi: obs_id=020326 energy=2.4456890899877 TeV  nroi(noedisp)=3.06495e-20 nroi(edisp)=8.25308e-18 ratio(edisp/noedisp)=269.273
GCTAResponseIrf::nroi: obs_id=020326 energy=2.84400060282542 TeV nroi(noedisp)=6.91589e-21 nroi(edisp)=1.39553e-18 ratio(edisp/noedisp)=201.787
GCTAResponseIrf::nroi: obs_id=020326 energy=3.84579880300243 TeV nroi(noedisp)=3.47822e-22 nroi(edisp)=5.57081e-20 ratio(edisp/noedisp)=160.163
GCTAResponseIrf::nroi: obs_id=020326 energy=4.47213595499958 TeV nroi(noedisp)=7.75419e-23 nroi(edisp)=1.28853e-20 ratio(edisp/noedisp)=166.172
GCTAResponseIrf::nroi: obs_id=020326 energy=3.30718220152506 TeV nroi(noedisp)=1.55508e-21 nroi(edisp)=3.95907e-19 ratio(edisp/noedisp)=254.59
GCTAResponseIrf::nroi: obs_id=020326 energy=5.20048006265587 TeV nroi(noedisp)=1.72326e-23 nroi(edisp)=2.6995e-21  ratio(edisp/noedisp)=156.651
GCTAResponseIrf::nroi: obs_id=020326 energy=6.04744425353316 TeV nroi(noedisp)=3.81852e-24 nroi(edisp)=4.37556e-22 ratio(edisp/noedisp)=114.588
GCTAResponseIrf::nroi: obs_id=020326 energy=7.03234731389669 TeV nroi(noedisp)=8.44092e-25 nroi(edisp)=7.62392e-23 ratio(edisp/noedisp)=90.3209
GCTAResponseIrf::nroi: obs_id=020326 energy=8.17765433957943 TeV nroi(noedisp)=1.86187e-25 nroi(edisp)=1.95125e-23 ratio(edisp/noedisp)=104.8
GCTAResponseIrf::nroi: obs_id=020326 energy=9.50948915243301 TeV nroi(noedisp)=4.09835e-26 nroi(edisp)=4.68426e-24 ratio(edisp/noedisp)=114.296
GCTAResponseIrf::nroi: obs_id=020326 energy=11.0582301703023 TeV nroi(noedisp)=9.00548e-27 nroi(edisp)=8.28975e-25 ratio(edisp/noedisp)=92.0523
GCTAResponseIrf::nroi: obs_id=020326 energy=12.859203321989 TeV  nroi(noedisp)=1.97543e-27 nroi(edisp)=2.23564e-25 ratio(edisp/noedisp)=113.172
GCTAResponseIrf::nroi: obs_id=020326 energy=14.9534878122122 TeV nroi(noedisp)=4.32235e-28 nroi(edisp)=3.68689e-26 ratio(edisp/noedisp)=85.2984
GCTAResponseIrf::nroi: obs_id=020326 energy=17.3888531156215 TeV nroi(noedisp)=9.44097e-29 nroi(edisp)=4.87566e-27 ratio(edisp/noedisp)=51.6437
GCTAResponseIrf::nroi: obs_id=020326 energy=20.2208485721784 TeV nroi(noedisp)=2.05524e-29 nroi(edisp)=3.43327e-28 ratio(edisp/noedisp)=16.705
GCTAResponseIrf::nroi: obs_id=020326 energy=23.5140704369771 TeV nroi(noedisp)=4.46263e-30 nroi(edisp)=1.14197e-28 ratio(edisp/noedisp)=25.5895
GCTAResponseIrf::nroi: obs_id=020326 energy=27.3436352852105 TeV nroi(noedisp)=9.65976e-31 nroi(edisp)=9.28774e-30 ratio(edisp/noedisp)=9.61488
GCTAResponseIrf::nroi: obs_id=020326 energy=31.7968933798403 TeV nroi(noedisp)=2.0785e-31  nroi(edisp)=1.37108e-30 ratio(edisp/noedisp)=6.59647
GCTAResponseIrf::nroi: obs_id=020326 energy=36.9754210829371 TeV nroi(noedisp)=4.45949e-32 nroi(edisp)=3.62126e-31 ratio(edisp/noedisp)=8.12034
GCTAResponseIrf::nroi: obs_id=020326 energy=42.9973377565041 TeV nroi(noedisp)=9.50461e-33 nroi(edisp)=8.54399e-32 ratio(edisp/noedisp)=8.98931
GCTAResponseIrf::nroi: obs_id=020326 energy=50 TeV               nroi(noedisp)=2.01651e-33 nroi(edisp)=9.74015e-33 ratio(edisp/noedisp)=4.83021

I tested the code on the H.E.S.S Crab data for the exponentially cut-off power law. Since energy dispersion leads at low energies to too high reconstructed energies (Ereco/Etrue > 1), neglecting energy dispersion makes the spectrum stepper, and taking it into account makes it flatter.

Below the spectra for the three cases:

The new code looks much better!

Thinks look okay and code has been merged in devel.

I also took the opportunity to rework the interface of the GCTAEdisp class, where all arguments that take a base 10 logarithm of an energy were replaced using GEnergy objects, which should reduce the number of computations. This affected the energy dispersion access operator, but also the mc() method as well as the energy boundary methods. These have also been renamed to ereco_bounds() and etrue_bounds(). Below the new interface. Please also note that all obs and src arguments were renamed into ereco and etrue, which is the new naming convention for reconstructed and true energies.

    // Constructors and destructors
    GCTAEdisp(void);
    GCTAEdisp(const GCTAEdisp& edisp);
    virtual ~GCTAEdisp(void);

    // Pure virtual operators
    virtual double operator()(const GEnergy& ereco,
                              const GEnergy& etrue,
                              const double&  theta = 0.0,
                              const double&  phi = 0.0,
                              const double&  zenith = 0.0,
                              const double&  azimuth = 0.0) const = 0;

    // Operators
    GCTAEdisp& operator=(const GCTAEdisp& edisp);

    // Pure virtual methods
    virtual void        clear(void) = 0;
    virtual GCTAEdisp*  clone(void) const = 0;
    virtual std::string classname(void) const = 0;
    virtual void        load(const GFilename& filename) = 0;
    virtual GFilename   filename(void) const = 0;
    virtual GEnergy     mc(GRan&          ran,
                           const GEnergy& etrue,
                           const double&  theta = 0.0,
                           const double&  phi = 0.0,
                           const double&  zenith = 0.0,
                           const double&  azimuth = 0.0) const = 0;
    virtual GEbounds    ereco_bounds(const GEnergy& etrue,
                                     const double&  theta = 0.0,
                                     const double&  phi = 0.0,
                                     const double&  zenith = 0.0,
                                     const double&  azimuth = 0.0) const = 0;
    virtual GEbounds    etrue_bounds(const GEnergy& ereco,
                                     const double&  theta = 0.0,
                                     const double&  phi = 0.0,
                                     const double&  zenith = 0.0,
                                     const double&  azimuth = 0.0) const = 0;
    virtual double      prob_erecobin(const GEnergy& ereco_min,
                                      const GEnergy& ereco_max,
                                      const GEnergy& etrue,
                                      const double&  theta) const = 0;
    virtual std::string print(const GChatter& chatter = NORMAL) const = 0;

The energy dispersion classes were revised, and I corrected the code so that the Monte Carlo simulations are compliant with the model. Below a number of figures that illustrate for the different energy dispersion types the consistency between Monte Carlo simulations and the model.

The comparisons between both a shown in the left panels, where red are Montre Carlo simulations and blue is the model. The vertical green line indicates the true photon energy.

The model for a given reconstructed energy as function of true photon energy is shown in the right panel. Here, the vertical green line indicates the reconstructed event energy.

Code merged into devel.