GammaLib

I did some additional tests:
- One pointing of 10,50,100h
- One point source in the field
- Unbinned analysis over 1-10TeV, also tried 100GeV-10TeV and 1-100TeV
- IRF .../data/cta/aar/bcf/DESY20140105_50h/irf_file.fits
- 50 fits in each case
- code version devel/bd092b9

Results are (model 1 is the point source):

True parameters (first background, then point source): [1.0, 0.0, 8.2e-19, -2.8]

100GeV-10TeV / 50h
Mean/deviation of background normalisation: 9.999464e-01 (9.115832e-04)
Mean/deviation of background index: 0.000000 (0.000000)
Mean/deviation of model 1 normalisation: 8.178346e-19 (1.263012e-20)
Mean/deviation of model 1 index: -2.801564 (0.011369)

1-100TeV / 50h
Mean/deviation of background normalisation: 1.000739e+00 (4.738377e-03)
Mean/deviation of background index: 0.000000 (0.000000)
Mean/deviation of model 1 normalisation: 8.105815e-19 (2.989788e-20)
Mean/deviation of model 1 index: -2.789380 (0.038785)

1-10TeV / 50h
Mean/deviation of background normalisation: 1.001022e+00 (4.984124e-03)
Mean/deviation of background index: 0.000000 (0.000000)
Mean/deviation of model 1 normalisation: 8.163553e-19 (3.259923e-20)
Mean/deviation of model 1 index: -2.800814 (0.046734)

I don’t see the problem anymore. There is no impact pivot energy of the background (1MeV or 1TeV, this is not used anyway). Not sure about what changed in the code since when I reported the problem. The mean parameters of the point source converges towards the true values if I increase the observing time:

1-10TeV / 100h
Mean/deviation of background normalisation: 1.000153e+00 (2.728596e-03)
Mean/deviation of background index: 0.000000 (0.000000)
Mean/deviation of model 1 normalisation: 8.201161e-19 (2.412416e-20)
Mean/deviation of model 1 index: -2.806473 (0.036199)

So overall, it looks like it is OK. I have not done more test than that but the case can probably be closed.

FORGET ABOUT THE PREVIOUS POST !

There must have been something wrong in the previous tests and I probably tested GCTAModelRadialAcceptance with ASCII IRF definition instead of GCTAModelInstBackground with FITS definition. With appropriate tests, I confirm there is an offset.

The setup now is:
- One pointing of 50h
- No source in the field
- Unbinned analysis over 1-10TeV
- ... also tried binned with 5 bins
- IRF .../data/cta/aar/bcf/DESY20140105_50h/irf_file.fits
- 50 fits in each case
- code version devel/bd092b9

The output is :
Mean/deviation of background normalisation: 1.004012e+00 (2.674310e-03)
... and the distribution of the 50 runs is shown in Unbinned_fits_50h.png

The situation is even worse for a binned analysis:
Mean/deviation of background normalisation: 1.020184e+00 (3.195677e-03)
... and the distribution also is worse, see Binned_fits_50h.png

Comparing that to other ways to set up the background:
- A gaussian for radial dependance and a power-law for spectrum: works perfectly for both binned and unbinned (tested only 50h)
- A gaussian for radial dependance and an ASCII file for spectrum: works for unbinned but not for binned (tested only 5 bins)
For the latter case, see plot Binned_file_50h.png

More testes needed...

Following up on that...

The test of background normalization when using a gaussian for radial dependance and an ASCII file for spectrum turns out to work well provided there are enough energy bins (something like 10-12 at least). See plot Binned_file_50h_12bins.png.

The problem of background normalisation when using GCTAModelInstBackground with the current set of instrument performances remains. It is most likely due to the coarse spatial binning of these responses, which leads to problems because the events simulation and the model-fitting handle this background differently: the former uses a cumulative distribution function, which has steps, while the latter uses an interpolation. When bins are too large, it is not surprising that the same number of counts is not obtained because there is no reason for it. A possible solution would be to finely rebin the background spatial distribution using the same function that is used for interpolation in model-fitting. Or wait for more accurate performance files.

Put this back into “In Progress” as Stefan Ohm was also stumbling over this in his simulations of Arp220. We should really find a fix for this.

Attached Stefan’s script that nicely produced the problem: simArp220.py

I recognized that the GCTAModelBackground constructor has optional parameters that allow for a rebinning. By default these parameters are 0, but you may add nx, ny, nenergies to the constructor to get a finer internal pixelisation of the background model, e.g.

    bck = gl.GCTAModelBackground(ob,filename,spec,100,100,100)

I did a couple of tests and come close to the true values. Higher resolution in nx and ny improves the source parameters, higher resolution in nenergies the background model parameters (see attached excel file).

I now implemented an internal rebinning to a finer grid into GCTAModelInstBackground. This gave reasonable results for the Arp220 test case. Attached an excel file and the test script used. Code is in the devel branch.

Please provide feedback and tell me if this also works for you.