GammaLib

I had actually not recognised how bad the situation really is. Below a plot for the current kludge implementation (clipping and slight smooth in etrue) where I changed the number of integration steps using the iter parameter (recall: iter is the fixed number of iterations in a Romberg integration). The plot shows the fractional difference in the model computation as function of energy with respect to a very high precision computation done with iter=10.

The current value that is used is iter=6 (red line), and it turns out that for that value there are spikes that reach about 15% of the reference value. Only for iter>=8 the integration converges, and uncertainties are at the level of a few percent.

Just to illustrate how bad it was before: at least the kludge goes in the right direction.

A 5 sigma smoothing in true energy instead of a 3 sigma smoothing does not really change the picture.

Going to a smoothing sigma of 10 pixels in etrue improves the situation as expected, but the accuracy with iter=6 is still not better than about 10% at high energies.

Smoothing bin 3 bins in etrue and 5 bins in migra improves the situation. The migra smoothing obviously has a strong impact on the required integration precision. This is certainly the way to follow. The biggest residuals are at the highest energies, where the statistic is lowest in the migration matrix.

The following figure illustrates the issue with the energy dispersion matrix. The noise in a given true energy bin is considerable, even for the domain where the statistic is high.

Here an attempt to get rid of the noise (blue curves).

The algorithm first determines for each true energy bin the mean and rms migration value since the zero order shape of the energy migration is a Gaussian (that’s in fact the blue curve shown in the plot above). The Gaussian is then subtracted from the migration values and the residual is then smoothed by a Gaussian with a certain width, and the smoothed residual is then added back to the Gaussian. In the plot below, the smoothed residual is shown as green line.

The trick now is to change the Gaussian smoothing sigma as function of the noise. For that purpose, the mean number of events per filled migration pixel mean is determined. This is hampered by the fact that the migration values have been rescaled, but by assuming that the minimum non-zero pixel value corresponds to one event one can get an estimate of the mean number of events per filled pixel.

Then the number of pixels num_pixels that are covered by the migration values is computed, computing the difference between the last and first non-zero migration pixels. To handle the case of a single filled pixel, we assume that at least 10% of the migration pixels are covered by the migration matrix.

Finally, the sigma (in number of pixels) is computed using

sigma = 0.5 * num_pixels / mean

It turned out that the current level of integration is needed to give accurate results.