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/00046361/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 pvalue that the data count spectrum is compatible with the model count spectrum (using the likelihood statistics).
The absolute value of PS is log10(pvalue) 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.
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, specifying the data and model 3d map files, as well as the energy binning
(with the arguments emin
, emax
and nbinloge
). It is recommended to set nbinloge
so that the logE bin width is between 0.2 and 0.3.
# Write the source model map (after performing the fit)
gta.write_model_map(model_name="model01")
# Generate the PS map
psmap = gta.psmap(ccmap='ccube_00.fits',mmap='mcube_model01_00.fits', make_plots=True, emin=100, emax=100000, nbinloge=15)
In the case of an analysis with multiple components, the user has to provide the list of data and model files (separated by colons) corresponding to the components that the user wants to include in the PS computation:
# Write the source model map (after performing the fit)
gta.write_model_map(model_name="model01")
# Generate the PS map using the first three components
psmap = gta.psmap(ccmap='ccube_00.fits:ccube_01.fits:ccube_02.fits',mmap='mcube_model01_00.fits:mcube_model01_01.fits:mcube_model01_02.fits', make_plots=True, emin=100, emax=100000, nbinloge=15)
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 pointsource 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).
In the case of a multiple component analysis, it is possible to use the argument ``psfpar0lst
to provide the list of psfpar0 of each component (separated by colons).
Option 
Default 
Description 


1 
Output verbosity. 

Data 3d map file or list of files (separated by colons). 


1 
PS computation: subROI halfwidth in pixel (ignoring central pixel; dpix=0 corresponds to the central pixel only). 

1000000.0 
PS computation: maximum energy/MeV. 

100 
PS computation: minimum energy/MeV. 

1.0 
Spatial integration: fixed radius (deg). If negative (default), the PSFlike parameters are used. 

1 
PS computation: subROI central pixel iaxis position. 

1 
PS computation: subROI central pixel jaxis position. 

False 
Generate diagnostic plots. 

100 
LL computation: number of counts up to which Poisson statistics is considered. 

Model 3d map file or list of files (separated by colons). 


20 
PS computation: number of log10 energy bins. It is recommended to set nbinloge so that the log10(E) bin width is between 0.2 and 0.3. 

50 
LL computation: number of bins. 

Output fits file, which contains the PS map (hdu0: log10(deviation probability), hdu1: in sigma units). 


1e07 
LL computation: precision parameter. 

4.0 
Spatial integration: PSFlike parameter 0 (deg) of the parameterization sqrt(p0^2*pow(E/p1,2*p2) + p3^2). 

Spatial integration: list of the PSFlike parameters 0 of all components (separated by colons). By default, psfpar0 is used for all components. 


100 
Spatial integration: PSFlike parameter 1 (MeV) of the parameterization sqrt(p0^2*pow(E/p1,2*p2) + p3^2). 

0.9 
Spatial integration: PSFlike parameter 2 of the parameterization sqrt(p0^2*pow(E/p1,2*p2) + p3^2). 

0.1 
Spatial integration: PSFlike parameter 3 (deg) of the parameterization sqrt(p0^2*pow(E/p1,2*p2) + p3^2). 

20 
LL computation: scale axis. 

Weight 3d map file or list of files (separated by colons). 


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 
psmaxsigma 
float 
maximum of the ps map in sigma units 
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 sigma, at %s=%.2f, %s=%.2f' %(psmap['psmaxsigma'],
psmap['coordname1'],float(psmap['coordx']),
psmap['coordname2'],float(psmap['coordy'])))
Maximum PSvalue=3.80 sigma, at GLONAIT=86.75, GLATAIT=38.62
Diagnostic plots can be generated by setting make_plots=True
or by
passing the output dictionary to make_psmap_plots
:
gta.write_model_map(model_name="model01")
psmap = gta.psmap(cmap='ccube_00.fits',mmap='mcube_model01_00.fits', 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 2.57,4.20,6.24 (corresponding to 3,4,5 sigma) indicating values above this threshold.image_psmap_sigma
: Map of PS values converted in sigma. The color map is truncated at 5 sigma with isocontours at 3,4,5 indicating values above this threshold.image_ps_hist
: Left: Histogram of PS (in sigma) for all points in the map. Superimposed is the reference distribution for a gaussian with mean 0 and sigma=1. Right: Histogram of the absolute value of PS (in sigma) for all points in the map. Superimposed is the reference distribution for the absolute value of a gaussian with mean 0 and sigma=1. It must be noted that, because the PS computation parameters are set to provide a good precision for low probability deviations and because of the sign assignment method, the central part of the PS distribution likely has a notch around zero. This notch is much less visible in the distribution of the absolute value of PS. What is really interesting in terms of significant deviation is the distribution for PS above 2 sigma.
PS Map 
PS [Sigma] Map 
PS [Sigma] 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 requires a model map to be computed, therefore the user must run gta.write_model_map(model_name=”model01”) before running psmap = gta.psmap(cmap=’ccube_00.fits’,mmap=’mcube_model01_00.fits’,emin=100,emax=100000,nbinloge=15)