Example 5: NPTF Correction for the Point Spread Function (PSF)

In this example we show how to account for the PSF correction using psf_correction.py

Fundamentally the presence of a non-zero PSF implies that the photons from any point source will be smeared out into some region around its true location. This effect must be accounted for in the NPTF. This is achieved via a function \(\rho(f)\). In the code we discretize \(\rho(f)\) as an approximation to the full function.

The two outputs of an instance of psf_correction are: 1. f_ary, an array of f values; and 2. df_rho_div_f_ary, an associated array of \(\Delta f \rho(f)/f\) values, where \(\Delta f\) is the width of the f_ary bins.

If the angular reconstruction of the data is perfect, then \(\rho(f) = \delta(f-1)\). In many situations, such as for the Fermi data at higher energies, a Gaussian approximation of the PSF will suffice. Even then there are a number of variables that go into evaluating the correction, as shown below. Finally we will show how the code can be used for the case of non-Gaussian PSFs.

As the calculation of \(\rho(f)\) can be time consuming, we always save the output to avoid recomputing the same correction twice. Consequently it can be convenient to have a common psf_dir where all PSF corrections for the runs are stored.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams

from NPTFit import psf_correction as pc # Module for determining the PSF correction

from __future__ import print_function

Example 1: Default Gaussian PSF

We start by showing the PSF correction for a Gaussian PSF - that is the PSF as a function of \(r\) is \(\exp \left[ -r^2 / (2\sigma^2) \right]\) - with \(\sigma\) set to the value of the 68% containment radius for the PSF of the Fermi dataset we will use in later examples.

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)
f_ary_1 = pc_inst.f_ary
df_rho_div_f_ary_1 = pc_inst.df_rho_div_f_ary

print('f_ary:', f_ary_1)
print('df_rho_div_f_ary:', df_rho_div_f_ary_1)

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='black', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.title('Gaussian PSF, $\sigma_\mathrm{PSF} = 0.1812$', y=1.04)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy
f_ary: [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]
df_rho_div_f_ary: [65.19815984  6.88897747  2.52908225  1.28920055  0.80522461  0.54317562
  0.09145394  0.          0.          0.        ]

Text(0.5,1.04,'Gaussian PSF, $\sigma_\mathrm{PSF} = 0.1812$')

Example 2: Impact of changing \(\sigma\)

Here we show the impact on the PSF of changing \(\sigma\). From the plot we can see that for a small PSF, \(\rho(f)\) approaches the no PSF case of \(\delta(f-1)\) implying that the flux fractions are concentrated at a single large value. As \(\sigma\) increases we move away from this idealized scenario and the flux becomes more spread out, leading to a \(\rho(f)\) peaked at lower flux values.

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.05)
f_ary_2 = pc_inst.f_ary
df_rho_div_f_ary_2 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.4)
f_ary_3 = pc_inst.f_ary
df_rho_div_f_ary_3 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='cornflowerblue',label='0.18', lw = 1.5)
plt.plot(f_ary_2,f_ary_2**2*df_rho_div_f_ary_2/(f_ary_2[1]-f_ary_2[0]),color='forestgreen',label='0.05', lw = 1.5)
plt.plot(f_ary_3,f_ary_3**2*df_rho_div_f_ary_3/(f_ary_3[1]-f_ary_3[0]),color='maroon',label='0.4', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='upper right', fancybox=True)
plt.title('Varying $\sigma_\mathrm{PSF}$', y=1.04)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.05_10_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.4_10_50000_1000_0.01.npy

Text(0.5,1.04,'Varying $\sigma_\mathrm{PSF}$')

Example 3: Changing the default options for determining \(\rho(f)\)

In this example we show how for a given PSF, the other parameters associated with how accurately we calculate \(\rho(f)\) can impact what we get back. The parameters that can be changed are are:

Argument	Defaults	Purpose
num_f_bins	10	number of f_bins used
n_psf	50000	number of PSFs placed down when calculating
n_pts_per_psf	1000	number of points to place per psf in the calculation
f_trunc	0.01	minimum flux fraction to keep track of
nside	128	nside of the map the PSF is used on

The default parameters have been chosen to be accurate enough for the Fermi analyses we will be performed later. But if the user changes the PSF (even just \(\sigma\)), it is important to be sure that the above parameters are chosen so that \(\rho(f)\) is evaluated accurately enough.

In general increasing num_f_bins, n_psf, and n_pts_per_psf, whilst decreasing f_trunc leads to a more accurate \(\rho(f)\). But each will also slow down the evaluation of \(\rho(f)\), and in the case of num_f_bin, slow down the subsequent non-Poissonian likelihood evaluation.

nside should be set to the value of the map being analysed, but we also highlight the impact of changing it below. For an analysis on a non-HEALPix grid, the PSF can often be approximated by an appropriate HEALPix binning. If this is not the case, however, a different approach must be pursued in calculating \(\rho(f)\).

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,num_f_bins=20)
f_ary_4 = pc_inst.f_ary
df_rho_div_f_ary_4 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,n_psf=5000,n_pts_per_psf=100)
f_ary_5 = pc_inst.f_ary
df_rho_div_f_ary_5 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,f_trunc=0.1)
f_ary_6 = pc_inst.f_ary
df_rho_div_f_ary_6 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,nside=64)
f_ary_7 = pc_inst.f_ary
df_rho_div_f_ary_7 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='black',label=r'Default', lw=2.2)
plt.plot(f_ary_4,f_ary_4**2*df_rho_div_f_ary_4/(f_ary_4[1]-f_ary_4[0]),color='forestgreen',label=r'more f\_bins', lw = 1.5)
plt.plot(f_ary_5,f_ary_5**2*df_rho_div_f_ary_5/(f_ary_5[1]-f_ary_5[0]),color='cornflowerblue',label=r'fewer points', lw = 1.5)
plt.plot(f_ary_6,f_ary_6**2*df_rho_div_f_ary_6/(f_ary_6[1]-f_ary_6[0]),color='salmon',label=r'larger f\_trunc', lw = 1.5)
plt.plot(f_ary_7,f_ary_7**2*df_rho_div_f_ary_7/(f_ary_7[1]-f_ary_7[0]),color='orchid',label=r'lower nside', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), fancybox=True)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_20_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_5000_100_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.1.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_64_0.181_10_50000_1000_0.01.npy

<matplotlib.legend.Legend at 0x7fe4c00337b8>

Example 4: PSF on a Cartesian Grid

For some applications, particularly when analyzing smaller regions of the sky, it may be desirable to work with data on a Cartesian grid rather than a healpix map. Note generally for larger regions, in order to account for curvature on the sky a healpix pixelization is recommended. Code to convert from Cartesian grids to healpix can be found here: https://github.com/nickrodd/grid2healpix

In order to calculate the appropriate PSF correction for Cartesian maps the general syntax is the same, except now the healpix_map keyword should be set to False and the pixarea keyword set to the area in sr of each pixel of the Cartesian map. In addition the gridsize keyword determines how large the map is, and flux that falls outside the map is lost in the Cartesian case.

As an example of this syntax we calculate the PSF correction on a Cartesian map that has pixels the same size as an nside=128 healpix map, and compare the two PSF corrections. Note they are essentially identical.

pixarea = 4*np.pi/(12*128*128)
pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812, healpix_map=False, pixarea=pixarea, gridsize=100)
f_ary_8 = pc_inst.f_ary
df_rho_div_f_ary_8 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='black', label=r'healpix', lw = 1.5)
plt.plot(f_ary_8,f_ary_8**2*df_rho_div_f_ary_8/(f_ary_8[1]-f_ary_8[0]),color='forestgreen', label=r'cartesian', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_0.21_0.181_10_50000_1000_0.01.npy

Text(0,0.5,'$f \times \rho(f)$')

Example 5: Custom PSF

In addition to the default Gausian PSF, psf_correction.py also has the option of taking in a custom PSF. In order to use this ability, the user needs to initialise psf_correction with delay_compute=True, manually define the parameters that define the PSF and then call make_or_load_psf_corr.

The variables that need to be redefined in the instance of psf_correction are:

Argument	Purpose
psf_r_func	the psf as a function of r, distance in radians from the center of the point source
sample_psf_max	maximum distance to sample the psf from the center, should be larger for psfs with long tails
psf_samples	number of samples to make from the psf (linearly spaced) from 0 to sample_psf_max, should be large enough to adequately represent the full psf
psf_tag	label the psf is saved with

As an example of a more complicated PSF we consider the full Fermi-LAT PSF. The PSF of Fermi is approximately Gaussian near the core, but has larger tails. To model this a pair of King functions are used to describe the radial distribution. Below we show a comparison between the Gaussian approximation and the full PSF, for two different energies. As shown, for low energies where the Fermi PSF is larger, the difference between the two can be significant. For higher energies where the PSF becomes smaller, however, the difference is marginal.

For the full details of the Fermi-LAT PSF, see: http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_PSF.html

# Fermi-LAT PSF at 2 GeV

# Calculate the appropriate Gaussian approximation to the PSF for 2 GeV
pc_inst = pc.PSFCorrection(psf_sigma_deg=0.2354)
f_ary_9 = pc_inst.f_ary
df_rho_div_f_ary_9 = pc_inst.df_rho_div_f_ary

# Define parameters that specify the Fermi-LAT PSF at 2 GeV
fcore = 0.748988248179
score = 0.428653790656
gcore = 7.82363229341
stail = 0.715962650769
gtail = 3.61883748683
spe = 0.00456544262478

# Define the full PSF in terms of two King functions
def king_fn(x, sigma, gamma):
    return 1./(2.*np.pi*sigma**2.)*(1.-1./gamma)*(1.+(x**2./(2.*gamma*sigma**2.)))**(-gamma)

def Fermi_PSF(r):
    return fcore*king_fn(r/spe,score,gcore) + (1-fcore)*king_fn(r/spe,stail,gtail)

# Modify the relevant parameters in pc_inst and then make or load the PSF
pc_inst = pc.PSFCorrection(delay_compute=True)
pc_inst.psf_r_func = lambda r: Fermi_PSF(r)
pc_inst.sample_psf_max = 10.*spe*(score+stail)/2.
pc_inst.psf_samples = 10000
pc_inst.psf_tag = 'Fermi_PSF_2GeV'
pc_inst.make_or_load_psf_corr()

# Extract f_ary and df_rho_div_f_ary as usual
f_ary_10 = pc_inst.f_ary
df_rho_div_f_ary_10 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_9,f_ary_9**2*df_rho_div_f_ary_9/(f_ary_9[1]-f_ary_9[0]),color='maroon',label='Gauss PSF', lw = 1.5)
plt.plot(f_ary_10,f_ary_10**2*df_rho_div_f_ary_10/(f_ary_10[1]-f_ary_10[0]),color='forestgreen',label='Fermi PSF', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='upper right', fancybox=True)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.235_10_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/Fermi_PSF_2GeV.npy

<matplotlib.legend.Legend at 0x7fe4b5217e10>

# Fermi-LAT PSF at 20 GeV

# Calculate the appropriate Gaussian approximation to the PSF for 20 GeV
pc_inst = pc.PSFCorrection(psf_sigma_deg=0.05529)
f_ary_11 = pc_inst.f_ary
df_rho_div_f_ary_11 = pc_inst.df_rho_div_f_ary

# Define parameters that specify the Fermi-LAT PSF at 20 GeV
fcore = 0.834725201378
score = 0.498192326976
gcore = 6.32075520959
stail = 1.06648424558
gtail = 4.49677834267
spe = 0.000943339426754

# Define the full PSF in terms of two King functions
def king_fn(x, sigma, gamma):
    return 1./(2.*np.pi*sigma**2.)*(1.-1./gamma)*(1.+(x**2./(2.*gamma*sigma**2.)))**(-gamma)

def Fermi_PSF(r):
    return fcore*king_fn(r/spe,score,gcore) + (1-fcore)*king_fn(r/spe,stail,gtail)

# Modify the relevant parameters in pc_inst and then make or load the PSF
pc_inst = pc.PSFCorrection(delay_compute=True)
pc_inst.psf_r_func = lambda r: Fermi_PSF(r)
pc_inst.sample_psf_max = 10.*spe*(score+stail)/2.
pc_inst.psf_samples = 10000
pc_inst.psf_tag = 'Fermi_PSF_20GeV'
pc_inst.make_or_load_psf_corr()

# Extract f_ary and df_rho_div_f_ary as usual
f_ary_12 = pc_inst.f_ary
df_rho_div_f_ary_12 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_11,f_ary_11**2*df_rho_div_f_ary_11/(f_ary_11[1]-f_ary_11[0]),color='maroon',label='Gauss PSF', lw = 1.5)
plt.plot(f_ary_12,f_ary_12**2*df_rho_div_f_ary_12/(f_ary_12[1]-f_ary_12[0]),color='forestgreen',label='Fermi PSF', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='upper left', fancybox=True)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.055_10_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/Fermi_PSF_20GeV.npy

<matplotlib.legend.Legend at 0x7fe4b51cccc0>

The above example also serves as a tutorial on how to combine various PSFs into a single PSF. In the case of the Fermi PSF the full radial dependence is the sum of two King functions. More generally if the full PSF is a combination of multiple individual ones (for example from multiple energy bins), then this can be formed by just adding these functions with an appropriate weighting to get a single psf_r_func.