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.

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(cmap='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(cmap='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 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). 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).

psmap Options

Option, Default, Description

chatter

1

Output verbosity.

cmap

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

dpix

-1

PS computation: sub-ROI half-width in pixel (ignoring central pixel; dpix=0 corresponds to the central pixel only).

emax

1000000.0

PS computation: maximum energy/MeV.

emin

100

PS computation: minimum energy/MeV.

fixedradius

-1.0

Spatial integration: fixed radius (deg). If negative (default), the PSF-like parameters are used.

ipix

-1

PS computation: sub-ROI central pixel i-axis position.

jpix

-1

PS computation: sub-ROI central pixel j-axis position.

make_plots

False

Generate diagnostic plots.

maxpoissoncount

100

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

mmap

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

nbinloge

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.

nbinpdf

50

LL computation: number of bins.

outfile

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

prob_epsilon

1e-07

LL computation: precision parameter.

psfpar0

4.0

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

psfpar0lst

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

psfpar1

100

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

psfpar2

0.9

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

psfpar3

0.1

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

scaleaxis

20

LL computation: scale axis.

wmap

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

write_fits

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

Map

WcsNDMap PSMAP

pssigma_map

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 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:

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

image_psmap

image_psmap_sigma

image_ps_hist_sigma

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)

Parameters:

  • prefix (str) – Optional string that will be prepended to all output files.

  • kwargs (Any) – these will override the psmap configuration file

  • {options}

Returns:

psmap – A dictionary containing the Map objects for PS amplitude and sigma equivalent.

Return type:

dict