Here I prototype some PSF models we are considering in Python to make a few plots for the internal note.

Some things are already implemented in gammalib: http://gammalib.sourceforge.net/doxygen/classGCTAPsf2D.html http://gammalib.sourceforge.net/doxygen/classGLATPsf.html

Everything we actually end up using should eventually be implemented properly in gammalib.

TODO:

Plot cumulative distributions, i.e. containment fraction
Plot containment radius, which is like containment fraction, but with x- and y-axis exchanged.

In [1]:

# Execute this if you didn't start this notebook with --pylab=inline
%pylab inline

Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

In [2]:

import numpy as np
from numpy import pi, sort, exp, log
import matplotlib.pyplot as plt

Definition of models

We define PSF models as probability density functions (PDFs) in 2D, i.e. as $\frac{dP}{dx dy}$. This is explained in more detail below.

TODO:

Implement as classes? Compute normalization integral only once!
Implement multi-Gauss and multi-King with correct normalisations.
Implement models in ROOT so that we can fit and draw random numbers easily!

In [3]:

def gauss(r, sigma, norm=1):
    """2D Gauss PDF dP / (dx dy) at r = sqrt(x ** 2 + y ** 2)
    Reference: TODO
    """
    N = norm / (2 * pi * sigma ** 2)
    return N * exp(-0.5 * r ** 2 / sigma ** 2)

In [4]:

def king(r, sigma, gamma, norm=1):
    """2D King PDF dP / (dx dy) at  r = sqrt(x ** 2 + y ** 2)
    Reference: This is the Fermi LAT PSF, see Equation (36) in [1]
    """
    N = norm / (2 * pi * sigma ** 2) * (1 - 1. / gamma)
    return N * (1 + 1. / (2 * gamma) * (r ** 2 / sigma ** 2)) ** (-gamma)

Transformation of probability densities

Definitions

There are three equivalent useful ways to describe a two-dimensional, radially symmetric PDF like a PSF.

$A(r) = \frac{dP}{dx dy}$ is the 2D probability density (unit: deg$^{-2}$) in $r = \sqrt{x^2 + y^2}$ (unit: deg). This is the distribution we used above to define our models, but it is rarely shown. A distribution with this shape will be found if histogramming events in $x$ for a small $y$-width stripe along the $x$-axis, or by evaluating the 2D PDF $P(x,y)$ at $P(x, y=0)$.
$B(t) = \frac{dP}{dt}$ is the 1D probability density (unit: deg$^2$) in $t = r^2$ (unit: deg$^2$). This is the theta square distribution that is often shown or fitted. A distribution with this shape will be found if histogramming events in $t$.
$C(r) = \frac{dP}{dr}$ is the 1D probability density (unit: deg) in $r= \sqrt{x^2 + y^2}$ (unit: deg). This distribution is not often shown or used, but is better suited to visualize differences between different distributions, because the peak is not at the edge. A distribution with this shape will be found if histogramming events in $r$.

Here we derive the relationship between $A$, $B$ and $C$ and illustrate them for the Gauss and King models. When implementing the fit function $B(t)$ for the theta square-distribution it might be better to re-code the formula for $B(t)$ directly.

Relations

Given $A(r)$, how to obtain / evaluate $B(t)$?

$t = r^2$, so $dt = 2 r dr$.
$dP = 2\pi r dr A(r)$ follows from the definition of $A$ as the probability $dP$ in a ring at radius $r$ of width $dr$.
$B(t) = \frac{dP}{dt} = \frac{2\pi r dr A(r)}{2 r dr} = \pi A(r)$, i.e. using the formula for $A(r)$ above we obtain the formula for $B(t)$ by substituting $r \to \sqrt{t}$ and multiplying by $\pi$.
$C(r) = \frac{dP}{dr} = \frac{2 \pi r dr A(r)}{dr} = 2 \pi r A(r)$, i.e. using the formula for $A(r)$ above we obtain the formula for $C(r)$ by multiplying with $2 \pi r$.
The relation between $B$ and $C$ is also simple, should we need it: $C(r) = 2 r B(t)$.

In [5]:

def compute_B(A, r2, *args):
    """B(r2) for a given A(r)
    See derivation above."""
    r = sqrt(r2)
    return pi * A(r, *args)

In [6]:

def compute_C(A, r, *args):
    """C(r) for a given A(r)
    See derivation above."""
    return 2 *pi * r * A(r, *args)

For the Gauss and King models we implement the $B$ and $C$ PDFs here explicitly in addition to the $A$ PDF given above.

This is useful for fitting and for comparing with other code that might use e.g. $B$ directly.

In [7]:

def gauss_B(r2, sigma, norm=1):
    """2D Gauss PDF dP(r2) / dr2 at r2 = x ** 2 + y ** 2"""
    N = norm / (2 * sigma ** 2)
    return N * exp(-0.5 * r2 / sigma**2)

In [8]:

def gauss_C(r, sigma, norm=1):
    """2D Gauss PDF dP(r) / dr at r = sqrt(x ** 2 + y ** 2)"""
    N = norm * r / sigma ** 2
    return N * exp(-0.5 * r ** 2 / sigma ** 2)

In [9]:

def king_B(r, sigma, gamma, norm=1):
    """2D King PDF dP(r2) / dr2 at r2 = x ** 2 + y ** 2"""
    N = norm / (2 * sigma ** 2) * (1 - 1. / gamma)
    return N * (1 + 1. / (2 * gamma) * (r2 / sigma ** 2)) ** (-gamma)

In [10]:

def king_C(r, sigma, gamma, norm=1):
    """2D King PDF dP(r) / dr at r = sqrt(x ** 2 + y ** 2)"""
    N = norm * r / sigma ** 2 * (1 - 1. / gamma)
    return N * (1 + 1. / (2 * gamma) * (r ** 2 / sigma ** 2)) ** (-gamma)

In [11]:

# Check normalizations
r_step = 0.01 # deg

for r_max in [1, 100]: # deg
    r = np.arange(0, r_max, r_step)
    pdf_gauss = gauss(r, sigma=0.2)
    pdf_king_2 = king(r, sigma=0.2, gamma=1.5)
    pdf_king_3 = king(r, sigma=0.2, gamma=3)
    for _ in [pdf_gauss, pdf_king_2, pdf_king_3]:
        norm = (2 * pi * r * r_step * _).sum()
        print norm

0.999787429515
0.67113834294
0.961792040763
0.999791640617
0.996466292832
0.999861093167

In [12]:

# Plot PDF dP / (dx dy)
# Note that this only contains 67 % of the events for the King profile (see previous cell)
r_max, r_step = 1, 0.01 # deg
r = np.arange(0, r_max, r_step)

plt.plot(r, gauss(r, sigma=0.2), label='Gauss$(\sigma=0.2)$');
plt.plot(r, king(r, sigma=0.1, gamma=1.5), label='King$(\sigma=0.1, \gamma=1.5)$');
plt.plot(r, king(r, sigma=0.2, gamma=3), label='King$(\sigma=0.2, \gamma=3)$');
plt.xlabel('r (deg)')
plt.ylabel('dP / (dx dy) (deg^-2)')
plt.legend(loc='best');

In [13]:

# Same plot as previous cell, but with log y axis scale
# With a log scale for y it can be clearly seen that at large r
# the Gauss profile is a parabola and the King profile is a line,
# i.e. has a power-law decrease.
plt.plot(r, gauss(r, sigma=0.2), label='Gauss$(\sigma=0.2)$');
plt.plot(r, king(r, sigma=0.1, gamma=1.5), label='King$(\sigma=0.1, \gamma=1.5)$');
plt.plot(r, king(r, sigma=0.2, gamma=3), label='King$(\sigma=0.2, \gamma=3)$');
plt.xlabel('r (deg)')
plt.ylabel('dP(r) / (dx dy) (deg^-2)')
plt.legend(loc='best')
plt.ylim(1e-3, None)
plt.semilogy();

In [14]:

# Check normalizations
r2_step = 0.001 # deg^2

for r2_max in [1, 1000]: # deg^2
    r2 = np.arange(0, r2_max, r2_step)
    pdf_gauss = gauss_B(r2, sigma=0.2)
    pdf_king_2 = king_B(r2, sigma=0.2, gamma=1.5)
    pdf_king_3 = king_B(r2, sigma=0.2, gamma=3)
    for _ in [pdf_gauss, pdf_king_2, pdf_king_3]:
        norm = (r2_step * _).sum()
        print norm

1.00625927081
0.674687757878
0.966684146376
1.0062630208
0.991133876887
1.0041752896

In [15]:

# Plot PDF dP(r^2) / dr^2
r2_max, r2_step = 1, 0.001 # deg
r2 = np.arange(0, r2_max, r2_step)

plt.plot(r2, gauss_B(r2, sigma=0.2), label='Gauss$(\sigma=0.2)$');
plt.plot(r2, king_B(r2, sigma=0.1, gamma=1.5), label='King$(\sigma=0.1, \gamma=1.5)$');
plt.plot(r2, king_B(r2, sigma=0.2, gamma=3), label='King$(\sigma=0.2, \gamma=3)$');
plt.xlabel('$r^2$ (deg$^2$)')
plt.ylabel('dP(r^2) / dr^2 (deg^-2)')
plt.legend(loc='best');

In [16]:

# Plot PDF dP(r^2) / dr^2
r2_max, r2_step = 1, 0.001 # deg
r2 = np.arange(0, r2_max, r2_step)

plt.plot(r2, gauss_B(r2, sigma=0.2), label='Gauss$(\sigma=0.2)$');
plt.plot(r2, king_B(r2, sigma=0.1, gamma=1.5), label='King$(\sigma=0.1, \gamma=1.5)$');
plt.plot(r2, king_B(r2, sigma=0.2, gamma=3), label='King$(\sigma=0.2, \gamma=3)$');
plt.xlabel('$r^2$ (deg$^2$)')
plt.ylabel('dP(r^2) / dr^2 (deg^-2)')
plt.legend(loc='best')
plt.ylim(1e-3, None)
plt.semilogy();

In [17]:

# Check normalizations
r_step = 0.01 # deg

for r_max in [1, 100]: # deg
    r = np.arange(0, r_max, r_step)
    pdf_gauss = gauss_C(r, sigma=0.2)
    pdf_king_2 = king_C(r, sigma=0.2, gamma=1.5)
    pdf_king_3 = king_C(r, sigma=0.2, gamma=3)
    for _ in [pdf_gauss, pdf_king_2, pdf_king_3]:
        norm = (r_step * _).sum()
        print norm

0.999787429515
0.67113834294
0.961792040763
0.999791640617
0.996466292832
0.999861093167

In [18]:

# Plot PDF dP(r) / dr
r_max, r_step = 1, 0.01 # deg
r = np.arange(0, r_max, r_step)

plt.plot(r, gauss_C(r, sigma=0.2), label='Gauss$(\sigma=0.2)$');
plt.plot(r, king_C(r, sigma=0.1, gamma=1.5), label='King$(\sigma=0.1, \gamma=1.5)$');
plt.plot(r, king_C(r, sigma=0.2, gamma=3), label='King$(\sigma=0.2, \gamma=3)$');
plt.xlabel('r (deg)')
plt.ylabel('dP / dr (deg^-1)')
plt.legend(loc='best');

In [19]:

# Plot PDF dP(r) / dr
r_max, r_step = 1, 0.01 # deg
r = np.arange(0, r_max, r_step)

plt.plot(r, gauss_C(r, sigma=0.2), label='Gauss$(\sigma=0.2)$');
plt.plot(r, king_C(r, sigma=0.1, gamma=1.5), label='King$(\sigma=0.1, \gamma=1.5)$');
plt.plot(r, king_C(r, sigma=0.2, gamma=3), label='King$(\sigma=0.2, \gamma=3)$');
plt.xlabel('r (deg)')
plt.ylabel('dP / dr (deg^-1)')
plt.legend(loc='best')
plt.ylim(1e-3, None)


plt.semilogy();

We want to use gammalib for most of our simulations / fitting.

At the moment the Fermi LAT PSF is fully implemented in gammalib, for HESS there is the GCTAPsf2D which implements a triple-Gauss model: http://gammalib.sourceforge.net/doxygen/classGCTAPsf2D.html

In [120]:

import gammalib

In [ ]:

# As far as I can see the PSF has to be initialized from a FITS file, which is not nice for simulations / checks ...
# This is an exmple PSF, I think for HD / Hillas and some Crab run, but I'm not sure ...
psf = gammalib.GCTAPsf2D('/Users/deil/code/gammalib/inst/cta/test/caldb/dc1/psf.fits')

In [ ]:

[np.degrees(psf.delta_max(logE)) for logE in np.arange(-1, 1.1, 0.5)]

In [ ]:

# This is how you can evaluate the psf, i.e. compute the PDF (in radians^(-2))
# at a given offset "delta" (in radians) and a given logE (in TeV?)
delta, logE = 0, 0
psf(delta, logE)

In [ ]:

# Let's plot the psf at 1 TeV, i.e. at 
x_max, x_step = 1, 0.01 # deg
x = np.arange(0, x_max, x_step)
y = np.zeros_like(x)
for ii in range(len(x)):
    # We use deg, gammalib uses radians, so convert before and after call
    y[ii] = (pi / 180) ** 2 * psf(np.radians(x[ii]), 0)

In [ ]:

# Check normalization
(2 * pi * x * x_step * y).sum()

In [ ]:

# Plot
plt.plot(x * x, y, label='gammalib');

plt.semilogy()
#plt.loglog()
plt.title('PDF of 2D distributions on x-axis')
plt.xlabel('x (deg)')
plt.ylabel('dP / dx (deg^-2)')
plt.legend(loc='best');

In [20]:

class Fermi_LAT_PSF(object):
    """Fermi LAT PSF model

    References:
    [1] http://adsabs.harvard.edu/abs/2012ApJS..203....4A
    [2] http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_PSF.html
    [3] http://gammalib.sourceforge.net/doxygen/classGLATPsfV3.html
    """

    def __init__(self, c0=3.32, c1=0.022, beta=0.80, sigma=0.1, gamma=0.2):
        self.scale_pars = dict(c0=c0, c1=c1, beta=beta)
        self.radial_pars = dict(sigma=sigma, gamma=gamma)

    @staticmethod
    def default():
        """Set default parameters"""
        # Front converting values from Table 13 in [1]
        c0, c1, beta = 3.32, 0.022, 0.80
        # Random values
        sigma, gamma = 0.1, 2
        return Fermi_LAT_PSF(c0, c1, beta, gamma, sigma)
        
    def scale_factor(self, energy):
        """Scale factor containing most of the energy dependence
        energy in MeV

        Reference: Equation (34) in [1].
        """
        p = self.scale_pars
        c0, c1, beta = p['c0'], p['c1'], p['beta']
        term1 = c0 * (energy / 100) ** (-beta)
        term2 = term1 ** 2 + c1 ** 2
        return sort(term2)

    def x(self, r, energy):
        """Scaled angular deviation x.
        Reference: Equation (35) in [1]
        """
        return r / scale_factor(energy)

    def K(self, r):
        """Radial PDF
        Reference: Equation (36) in [1]
        """
        p = self.radial_pars
        sigma, gamma = p['sigma'], p['gamma']
        return king(r, sigma, gamma)

PSF Models

Definition of models

Transformation of probability densities

Definitions

Relations

Checks / Plots

For dP(r) / (dx dy)

For dP(r^2) / dr^2

For dP(r) / dr

Gammalib

Backup Cells