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 <data 3D map> --mmap <model 3D map> --wmap <weights 3D map (optional)> --emin <minimum energy/MeV> --emax <maximum energy/MeV> --nbinloge <number of log10 energy bins> --outfile <output file>

Example for one component between 100 MeV and 100 GeV:
gtpsmap.py --cmap ccube_00.fits --mmap mcube_00.fits --emin 100 --emax 100000 --nbinloge 15 --outfile psmap.fits
Example for four components between 100 MeV and 100 GeV:
gtpsmap.py --cmap ccube_00.fits:ccube_01.fits:ccube_02.fits:ccube_03.fits --mmap mcube_00.fits:mcube_01.fits:mcube_02.fits:mcube_03.fits --emin 100 --emax 100000 --nbinloge 15 --outfile psmap.fits

The output file contains:
- primary hdu: 2D map of PS in sigma (with sigma = sqrt(2) * ErfInv(1-10^(-|PS|)))
- hdu #1: 2D map of PS

The definition of PS at each pixel of the map is:
- |PS| = -log10(deviation probability), where the deviation probability is measured on the PSF-like integrated count data and model spectra
- sign(PS) = sign of sum(data_k-model_k)/sqrt(max(1,model_k)/weight_k), where the sum runs over the spectral bins

Here are some useful PS<->sigma conversions:
- PS -> sigma: 3/4/5/6/7/8 -> 3.29/3.89/4.42/4.89/5.33/5.73
- sigma -> PS: 3/4/5/6 -> 2.57/4.20/6.24/8.70

The other (optional) parameters of the python script are:
--psfpar0 <PSF-like par0; 4.0 deg> --psfpar1 <PSF-like par1; 100 MeV>  --psfpar2 <PSF-like par2; 0.9>  --psfpar3 <PSF-like par3; 0.1deg>
  these parameters define the PSF-like integration radius = sqrt(p0^2*pow(E/p1,-2*p2) + p3^2)
  The default values (4.0,100,0.9,0.1) are optimal for a point-like deviation. For potential extended deviations,
  it is advised to set psfpar3 to about the size of the extension (e.g. 0.5 or 1 deg)
--psfpar0lst: in the case of multiple components, the user can set psfpar0 independently for each component, e.g. --psfpar0lst 2.0:3.0:4.0:5.0
--chatter <output verbosity; 1>
--prob_epsilon <precision parameter; 1e-7>
  This parameter defines the k-interval over which the Poisson probability of each spectral bin is considered
  when computing the log-likelihood distribution (the k-interval is such that ΣPoisson(k,model_i) = 1-prob_epsilon).
  The default value 1e-7 provides a 1% precision on PS up to PS=20. If you want a better precision for larger PS values,
  you need a smaller prob_epsilon (1e-11 provides a 5% precision up to 250) but that makes the script slower.