GammaLib

I checked whether the residuals change with a changed number of energy bins. The figures below show the latitude residuals for 20, 30, 40 and 60 energy bins. Apparently the residuals do not reduce with an increase in the number of energy bins.

20 bins	30 bins
40 bins	60 bins

I did a single simulation of 50h of background only on the Galactic centre, generated response and background cubes, computed a model cube and displayed the residual between data and model. I selected an energy range of 70 GeV - 199 TeV and used the prod3b-v2 IRF South_z20_50h. I used 20 energy bins.

Below the residuals counts cube - model cube for the first energy bin (70-104.2 GeV). There is obviously an excess in the residuals that is not expected.

For comparison the residual that is obtained when the IRF response instead of the stacked response is used. This looks better, but is not perfect.

Here the same result for 100 energy bins, done using the IRF response computation. This is now a better, although there is still some structure visible, but maybe related to poor statistic.

I made a modification in GCTAModelIrfBackground::eval() using now the GCTABackground3D::rate_ebin() integrator method instead of the GCTABackground3D::operator() operator. This takes a bit more time (ctmodel takes 4.9 seconds instead of 3.6 seconds before, which is actually not so bad). The residuals are now clean.

I simplified the GCTACubeBackground::fill() method, removing the power law integration, and using directly the GCTAModelIrfBackground::eval() method which is now exact. Recomputing the background cube using ctbkgcube and then computing a model using the stacked response using ctmodel produced the following plot. It looks perfect, and actually identical to the plot generated using the IRF response.

I redid the GPS anticentre analysis. Below the result, top is original run, bottom the modified code. Things look now much better!

Also note that the background model fit parameters are now much closer to nominal. The normalisation of the IEM is also impacted, not being consistent with the nominal value of 1.

Before:

2020-11-30T16:23:18: === GCTAModelCubeBackground ===
2020-11-30T16:23:18:  Name ......................: BackgroundModel
2020-11-30T16:23:18:  Instruments ...............: CTA, HESS, MAGIC, VERITAS
2020-11-30T16:23:18:  Observation identifiers ...: all
2020-11-30T16:23:18:  Model type ................: "PowerLaw" * "Constant" 
2020-11-30T16:23:18:  Number of parameters ......: 4
2020-11-30T16:23:18:  Number of spectral par's ..: 3
2020-11-30T16:23:18:   Prefactor ................: 0.95425308484162 +/- 0.00401027199310809 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)
2020-11-30T16:23:18:   Index ....................: -0.0242787890275514 +/- 0.00147495295160353 [-5,5]  (free,scale=1,gradient)
2020-11-30T16:23:18:   PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
2020-11-30T16:23:18:  Number of temporal par's ..: 1
2020-11-30T16:23:18:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2020-11-30T16:23:18: === GCTAModelSkyCube ===
2020-11-30T16:23:18:  Name ......................: IEM
2020-11-30T16:23:18:  Instruments ...............: all
2020-11-30T16:23:18:  Observation identifiers ...: all
2020-11-30T16:23:18:  Model type ................: "ConstantValue" * "Constant" 
2020-11-30T16:23:18:  Number of parameters ......: 3
2020-11-30T16:23:18:  Number of spectral par's ..: 1
2020-11-30T16:23:18:   Value ....................: 16.6818659128961 +/- 1.41465490793952 [1e-05,100000]  (free,scale=1,gradient)
2020-11-30T16:23:18:  Number of temporal par's ..: 1
2020-11-30T16:23:18:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

2020-12-01T23:09:55: === GCTAModelCubeBackground ===
2020-12-01T23:09:55:  Name ......................: BackgroundModel
2020-12-01T23:09:55:  Instruments ...............: CTA, HESS, MAGIC, VERITAS
2020-12-01T23:09:55:  Observation identifiers ...: all
2020-12-01T23:09:55:  Model type ................: "PowerLaw" * "Constant" 
2020-12-01T23:09:55:  Number of parameters ......: 4
2020-12-01T23:09:55:  Number of spectral par's ..: 3
2020-12-01T23:09:55:   Prefactor ................: 0.998666663749527 +/- 0.00463828188617879 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)
2020-12-01T23:09:55:   Index ....................: -0.000905912390236462 +/- 0.00158187474876209 [-5,5]  (free,scale=1,gradient)
2020-12-01T23:09:55:   PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
2020-12-01T23:09:55:  Number of temporal par's ..: 1
2020-12-01T23:09:55:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2020-12-01T23:09:55: === GCTAModelSkyCube ===
2020-12-01T23:09:55:  Name ......................: IEM
2020-12-01T23:09:55:  Instruments ...............: all
2020-12-01T23:09:55:  Observation identifiers ...: all
2020-12-01T23:09:55:  Model type ................: "ConstantValue" * "Constant" 
2020-12-01T23:09:55:  Number of parameters ......: 3
2020-12-01T23:09:55:  Number of spectral par's ..: 1
2020-12-01T23:09:55:   Value ....................: 1.83811471079993 +/- 1.65034103271708 [1e-05,100000]  (free,scale=1,gradient)
2020-12-01T23:09:55:  Number of temporal par's ..: 1
2020-12-01T23:09:55:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

I moved the issue to the bugfix-1.7.3 release.

Note that I had also added GCTAEventBin::emin() and GCTAEventBin::emax() methods to compute the event boundaries for an event bin from the log mean energy and energy bin width. I kept these methods in release 2.0.0 but removed them from release 1.7.3 since a bugfix release should have no interface change.