PS Map¶
psmap()
quantifies the 3d data/model agreement by computing the PS
at each spatial bin in the ROI according to the algorithm described in Bruel P. (2021), A&A, 656, A81.
(doi:10.1051/0004-6361/202141553). 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.
Caveat: the algorithm is currently not able to run in the context of a joint analysis with several event types.
Examples¶
In order to run psmap, the user must first compute the source model map using the write_model_map'
function specifying the name of the model with the parameter ``model_name`()
.
The psmap()
will then be called with the same model_name
.
# Write the source model map (after performing the fit)
gta.write_model_map(model_name="model01")
# Generate the PS map
psmap = gta.psmap(model_name='model01', make_plots=True)
Configuration¶
The default configuration of the method is controlled with the psmap section of the configuration file. The default configuration can be overriden by passing the option as a kwargs argument to the method.
If the analysis uses likelihood weights, the user can specify the likelihood weight file with the argument wmap
.
The user can set the parameters defining the energy dependent region used in the count spectrum integration step. The integration radius is defined by the function sqrt(psfpar0^2*pow(E/psfpar1,-2*psfpar2) + psfpar3^2), with E in MeV. The default parameters (``psfpar0``=4.0,``psfpar1``=100,``psfpar2``=0.9,``psfpar3``=0.1) are optimized to look for point-source like deviations. When looking for extended deviations (~1deg scale), it is recommended to use (``psfpar0``=5.0,``psfpar1``=100,``psfpar2``=0.8,``psfpar3``=1.0).
The user can set the energy range with the arguments emin
and emax
. Depending on the energy binning, it can be optimal
to rebin the count spectra thanks to the argument rebin
. For an analysis with 10 bins per decade, it is recommended to set rebin
to 3.
Option |
Default |
Description |
---|---|---|
|
1 |
output verbosity |
|
1000000000.0 |
maximum energy/MeV |
|
1.0 |
minimum energy/MeV |
|
-1.0 |
Fixed radius |
|
-1 |
number of pixel i axis |
|
-1 |
number of pixel j axis |
|
False |
Generate diagnostic plots. |
|
100 |
Maximum number of poisson counts |
|
None |
Model name |
|
50 |
Number of bin of the PSF |
|
psmap.fits |
outfile name |
|
1e-07 |
precision parameter |
|
4.0 |
PSF parameter 0 |
|
100 |
PSF parameter 1 |
|
0.9 |
PSF parameter 2 |
|
0.1 |
PSF parameter 3 |
|
1 |
Rebin |
|
20 |
SCale axis |
|
weight 3d map |
|
|
True |
Write the output to a FITS file. |
Output¶
psmap()
returns a dictionary containing the following variables:
Key |
Type |
Description |
---|---|---|
psmax |
float |
maximum of the ps map |
coordname1 |
str |
Name of the coordinate of the ps maximum |
coordname2 |
str |
Name of the coordinaste of the ps maximum |
coordx |
float |
Value of the X coordinate of the ps maximum |
coordy |
float |
Value of the Y coordinate of the ps maximum |
ipix |
int |
Value of the i pixel of the of the ps maximum |
jpix |
int |
Value of the j pixel of the of the ps maximum |
wcs2d |
WCS Keywords |
WCS of the ps map |
psmap |
np.array |
PSMAP |
psmapsigma |
np.array |
PSMAP in sigma units |
name |
str |
NAmke of the model |
ps_map |
WcsNDMap PSMAP |
|
pssigma_map |
WcsNDMap PSMAP in sigma units |
|
config |
dict |
Dictionary of the input configuration |
file |
str |
Name of the output file |
file_name |
str |
Full path of the output file |
The write_fits
option can used to write the output to a FITS or numpy file. The value of the maximum of the PS map
can be retrieved from the output dictionary:
print('PS maximum value=%.2f, at %s=%.2f, %s=%.2f' %(psmap['psmax'],
psmap['coordname1'],float(psmap['coordx']),
psmap['coordname2'],float(psmap['coordy'])))
PS maximum value=3.85, at GLON-AIT=86.75, GLAT-AIT=38.62
Diagnostic plots can be generated by setting make_plots=True
or by
passing the output dictionary to make_psmap_plots
:
psmap = gta.psmap(model_name='model01', make_plots=True)
//equivalent to:
gta.plotter.make_tsmap_plots(psmap, roi=gta.roi)
This will generate the following plots:
image_psmap
: Map of PS values. The color map is truncated at 5 sigma with isocontours at 3,4,5 PS intervals indicating values above this threshold.image_pssigma
: Map of PS values converted in sigma. The color map is truncated at 5 sigma with isocontours at 3,4,5 PS intervals indicating values above this threshold.image_ps_hist
: Histogram of PS values for all points in the map. Overplotted is the reference distribution for a gaussian with mean 0 and sigma=1.
PS Map |
Sigma (PS) Map |
PS Histogram |
---|---|---|
Reference/API¶
- GTAnalysis.psmap(prefix='', **kwargs)
Generate a spatial PS map for evaluate the data/model comparison The PS map will have the same geometry as the ROI. The output of this method is a dictionary containing
Map
objects with the PS amplitude and sigma equivalent. By default this method will also save maps to FITS files and render them as image files. psmap requirtes a model map to be compute, therefore the user must run gta.write_model_map(model_name=”model01”) before running psmap = gta.psmap(model_name=’model01’)