gtpsmap: check data/model agreement in Fermi-LAT analysis
This python script compares data and model 3D maps by computing a PS map.
It has been written by Philippe Bruel (LLR, CNRS/IN2P3, IPP Paris, France), based on a new method fully described in:
https://www.aanda.org/articles/aa/full_html/2021/12/aa41553-21/aa41553-21.html
For each spatial bin, the algorithm first computes the data and model count spectra integrated over an energy dependent region (following the PSF energy dependence) and then computes the p-value that the data count spectrum is compatible with the model count spectrum (using the likelihood statistics). The absolute value of PS is -log10(p-value) and its sign corresponds to the sign of the sum over the spectrum of the residuals in sigma. The algorithm also provides a PS map in sigma units.
WARNING: contrary to the previous version, the script is now able to handle multi-component analyses. This modification required that the user has to provide the energy range and binning.
The main inputs of the script are:
- the data 3D count map (output of gtbin) or the list of data 3D count maps for multi-component analyses
- the model 3D count map (output of gtmodel or gta.write_model_map() for fermipy users) or the list of model 3D count maps for multi-component analyses
- if you use weighted likelihood to perform the fit, then you also have to provide the 3D weight map
The user also has to provide the energy interval and binning (emin, emax, nbinloge) that are used to compute PS (which can be different from the 3D map binning). The recommendation is to set nbinloge in order to have a logE bin width in the range 0.2-0.3.
Usage:
gtpsmap.py --cmap --mmap --wmap --emin --emax --nbinloge --outfile