Fitting for a point source polarization with LeakageLib

This jupyter notebook fits for the polarization of a single point source, with a photon and particle background. Spatial, particle, and spectral weights are used.

import numpy as np
import leakagelib
from ixpeobssim.irf import load_arf

>>> PyXSPEC is not installed, you will no be able to use it.

First, we need to load the data. LeakageLib needs both the level 2 event file and the housekeeping (hk) attitude file in order to run. (The attitude file contains the spacecraft orientation which is necessary to get the PSF model right). I made mock IXPE data for this example, which we’ll use.

source = leakagelib.source.Source.no_image(False) # Tells LeakageLib that the data was NOT NN reconstructed.
datas = [leakagelib.IXPEData(source, (
    "data/ps/event_l2/ixpeps_det1_evt2_v00.fits",
    "data/ps/hk/ixpeps_det1_att_v00.fits",
), energy_cut=(2,8))]

Alternatives to leakagelib.IXPEData

• leakagelib.IXPEData.load_all_detectors. Searches all directories listed in the DATA_DIRECTORIES variable defined in settings.py for the given obs id, and load all detectors.

• leakagelib.IXPEData.load_all_detectors_with_path. Provide a path to the observation and the function will load all detectors. This will not use the DATA_DIRECTORIES variable.

• leakagelib.IXPEData (this method). You must pass specific file names.

Next we will center the data and cut to data within 280 arcsec of the center.

for data in datas:
    data.iterative_centroid_center()
    data.retain(np.sqrt(data.evt_xs**2 + data.evt_ys**2) < 280)

Alternatives to leakagelib.IXPEData.iterative_centroid_center

• leakagelib.IXPEData.explicit_center: Center on a specific point, in pixels.

• leakagelib.IXPEData.iterative_centroid_center: Center by iteratively zooming in on the average event. This finds central PSs rather well.

Alternatives to leakagelib.IXPEData.retain

• leakagelib.IXPEData.retain_region: Pass in a region file to cut only events in the region. Note: the region must be ciao formatted in physical coordinates

• leakagelib.IXPEData.retain: Simply a list of bools, one per event, where True indicates the event should be cut.

Now we have to tell the fitter what the image looks like. There are two components: a point source and a photon background. To create the point source,

settings = leakagelib.FitSettings(datas)

settings.add_point_source() # Named "src" by default
settings.fix_flux("src", 1)
settings.set_initial_qu("src", (0.5,0))

Without loss of generality, we have defined the flux to be in units of the source flux.

You can use leakagelib.FitSettings.set_initial_qu to set the start point of the fitter. It’s important to set the start point as close to the true value as possible to avoid falling into the wrong maximum of the likelihood.

Now we create the background models:

settings.add_background() # Named "bkg" by default
settings.fix_qu("bkg", (0, 0)) # Assume an unpolarized background
settings.add_particle_background() # Named "pbkg" by default
settings.fix_qu("pbkg", (0, 0)) # Assume the particles have no direction dependence

These background sources are assumed to be spatially uniform and already have spatial weights applied. To use energy weights, we must assign a spectrum. To do this, we create functions which take in the energy in keV and return the photon spectrum (dN/dE). Normalization doesn’t matter, but remember to include the ARF.

arf = load_arf()
settings.set_spectrum("bkg", lambda e: arf(e) * e**-2.5) # Gamma=2.5, unabsorbed powerlaw
settings.set_spectrum("src", lambda e: arf(e) * e**-1.5) # Gamma=1.5, unabsorbed powerlaw

>>> Reading (in memory) /opt/homebrew/anaconda3/lib/python3.12/site-packages/ixpeobssim/caldb/ixpe/gpd/cpf/arf/ixpe_d1_obssim20240101_v013.arf...

We have assigned all the weights. We now have to tell the fitter which events we cut, so that it knows how to normalize the spatial weights. Since the ROI is circular, this is easy

settings.apply_circular_roi(280)

Alternatives to leakagelib.FitSettings.apply_circular_roi

• leakagelib.IXPEData.apply_circular_roi: Assumes a circular ROI centered on the origin

• leakagelib.IXPEData.apply_roi: Apply an arbitrary ROI, which you pass in as an image. The image must have the same pixel size and dimensions as the fit sources. To get these properties, check out settings.sources[0] (example below).

print(type(settings.sources[0]))
print(settings.sources[0].pixel_size)
print(settings.sources[0].source.shape)

<class 'leakagelib.source.Source'>
2.9729
(195, 195)

Now we are ready to run the fit.

fitter = leakagelib.Fitter(datas, settings)
print(fitter)
result = fitter.fit()

result

FITTED PARAMETERS:
Source  Param
src:    q
src:    u
bkg:    f
pbkg:   f

FIXED PARAMETERS:
Source  Param   Value
src:    f       1
bkg:    q       0
bkg:    u       0
pbkg:   q       0
pbkg:   u       0

FitResult:
        q (src) = 0.5092 +/- 0.0721
        u (src) = -0.0001 +/- 0.0723
        f (bkg) = 1.4602 +/- 0.0315
        f (pbkg) = 0.6852 +/- 0.0177

Polarization:
        PD (src): 0.5092 +/- 0.0721 (7.1 sig)
        PA (src): -0.0082 deg +/- 4.0700
Likelihood 23123.702880749537, dof 15691
Optimization terminated successfully.

Indeed, the true polarization of this mock source was q = 0.5, u = 0.0. The fit result agrees.

Alternatives to leakagelib.Fitter.fit

• leakagelib.Fitter.fit: Fit by numerically maximizing the likelihood and get Gaussian uncertainties by taking the Hessian and using Laplace’s approximation.

• leakagelib.IXPEData.fit_mcmc: Do a full MCMC fit. This code is less tested and you may need to edit it yourself.

You can access result.params and result.sigmas to get the best-fit values and uncertainties. Also result.cov to get the full covariance matrix.

print("Parameters", result.params)
print("Uncertainties", result.sigmas)
print("Covariance", result.cov)

Parameters {('q', 'src'): 0.5091997559169119, ('u', 'src'): -0.0001465637380886765, ('f', 'bkg'): 1.4602052476296876, ('f', 'pbkg'): 0.6851993225880586}
Uncertainties {('q', 'src'): 0.0720889216211086, ('u', 'src'): 0.072342012166437, ('f', 'bkg'): 0.03145924839584101, ('f', 'pbkg'): 0.017668913696187224}
Covariance [[ 5.19681262e-03 -8.08105348e-05  0.00000000e+00  0.00000000e+00]
 [-8.08105348e-05  5.23336672e-03  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  9.89684310e-04  1.96057259e-04]
 [ 0.00000000e+00  0.00000000e+00  1.96057259e-04  3.12190511e-04]]