cta-observatory / magic-cta-pipe

Pipeline for the analysis of MAGIC and LST1 data, and more.
https://magic-cta-pipe.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
5 stars 9 forks source link

Gammapy and energy bias #130

Open Elisa-Visentin opened 1 year ago

Elisa-Visentin commented 1 year ago

Somebody knows WITH CERTAINTY if (and where) the energy bias is corrected by Gammapy in check_flux_from_dl3.ipynb? We (me and some other LST member, they can reveal themselves if they want to) modified a real DL3 FITS file so that each event has E=0.12 TeV (but no changes to migration matrix, effective area...) and, by using this new FITS file in the notebook, there seems to be no energy fixing / correction through the cells up to SED. Maybe I made a mistake but I think we should know (or check) this. notebook: /fefs/aswg/workspace/elisa.visentin/myscripts/notebooks/bias_check_2.ipynb image

jsitarek commented 1 year ago

I do not know with certainty, but I expect that it works in a different way since it is using a joint likelihood fit . The main problem is not the bias but the energy resolution, so you cannot simply shift all the energies by a given factor. Instead when you make a fit you start with a particular shape of the spectrum, and fold it with the energy dispersion (that takes into account automatically the energy bias and energy resolution) to see how it is distributed in estimated energies and compare it with the actually measured distribution of events in estimated energy bins. Then you start modifying the fit parameters to minimize the difference and hence get the final fitted values.

Then if you want to get the points: https://docs.gammapy.org/0.18.2/tutorials/spectrum_analysis.html#Compute-Flux-Points the code is used to "compute flux points by fitting the norm of the global model in energy bands" so it seems that it is using the same kind of procedure (including the energy resolution and bias correction via migration matrix). I think this should be estimated energy, so then you try to reproduce the normalization of the previously derived spectral model that gives the total number of events in that estimated energy bin that is measured. If it works like this than of course it can happen that with the bias you use ~ 80 GeV events to estimate the flux at 100 GeV, but as long as the spectral model is correct the estimation is still correct (i.e. the points shift along the spectrum of the source). I do not know (and do not think) that there is any additional correction shifting the points again back along the spectral line of the source (like we do in flute) - this would probably show up in your test. If you still have some doubts it is best to ask it directly to gammapy developers.