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 described in a publication currently in preparation.
Meanwhile one can look at the corresponding Fermi symposium 2021 presentation:
https://indico.cern.ch/event/1010947/contributions/4278096/

The main inputs of the script are:
- the data 3D count map of the data (output of gtbin)
- the model 3D count map (output of gtmodel or gta.write_model_map() for fermipy users)
- if you use weighted likelihood to perform the fit, then you also have to provide the 3D weight map

Usage:
gtpsmap.py --cmap <data 3D map> --mmap <model 3D map> --wmap <weights 3D map (optional)> --outfile <output file>

The output file contains:
- primary hdu: 2D map of PS
- hdu #1 ('PS in sigma'): deviation in sigma units (with sigma = sqrt(2) * ErfInv(1-10^(-|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:
--emin <minimum energy/MeV; 1> --emax <maximum energy/MeV; 1e9>
  energy range to compute PS
--rebin <nb of bins merged when rebinning the energy axis; 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.
--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>
  PSF-like selection radius = sqrt(p0^2*pow(E/p1,-2*p2) + p3^2)
--chatter <output verbosity; 1>