Fitting for an extended source polarization with LeakageLib

This jupyter notebook fits for the polarization of an extended source, with a photon and particle background. Here, the source is a 60 arcsecond-radius uniform disk, though because of the IXPE PSF, it appears wider in the IXPE data.

The source is assumed to have constant polarization, though we provide guidance on variable polarization sources at the end. Spatial, particle, and spectral weights are used.

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

The first step is to load the mock Moments data. We’ll center it on pixel (300, 300), and cut to 280 arcseconds around that point. See the point source fit example for a description of these functions.

source = leakagelib.Source.no_image(False)

# Load the extended source data
datas = [leakagelib.IXPEData(source, (
    "data/extended/event_l2/ixpeextended_det1_evt2_v00.fits",
    "data/extended/hk/ixpeextended_det1_att_v00.fits",
), energy_cut=(2,8))]

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

Now we have to create a custom Source object for our extended source: a 60”-radius disk. The source object should not be blurred by the PSF; i.e. it should depict the source as seen by Chandra, not IXPE

pixel_size = 2.9729 # This is the size of the sky PSF pixels; use it for best results
source_pixel_edges = np.arange(0, 280+pixel_size, pixel_size)
source_pixel_edges = np.concatenate([-np.flip(source_pixel_edges)[:-1], source_pixel_edges])
# source_pixel_edges now ranges from -280 to 280, centered on zero, and spaced by `pixel_size`
source_pixel_centers = (source_pixel_edges[1:] + source_pixel_edges[:-1]) / 2

xs, ys = np.meshgrid(source_pixel_centers, source_pixel_centers)
image = (np.sqrt(xs**2 + ys**2) < 60).astype(float)

Now we give the source object to the fitter

source = leakagelib.Source(image, False, len(image), pixel_size)
settings = leakagelib.FitSettings(datas)
settings.add_source(source, "src-ext")
settings.fix_flux("src-ext", 1)

The background and spectral weights are added in the same way as the previous fit, but this time we won’t assume zero background polarization.

We are assuming no particles for convenience. This is what you should do if you did not generate particle weights with leakagelib_bkg.

settings.add_background("bkg")
settings.set_initial_flux("bkg", 0.1) # Use this function to set the initial value of the background flux

arf = load_arf()
settings.set_spectrum("bkg", lambda e: arf(e) * e**-2.5)
settings.set_spectrum("src-ext", lambda e: arf(e) * e**-1.5)

settings.apply_circular_roi(280)

fitter = leakagelib.Fitter(datas, settings)

>>> Using cached xEffectiveArea object at /opt/homebrew/anaconda3/lib/python3.12/site-packages/ixpeobssim/caldb/ixpe/gpd/cpf/arf/ixpe_d1_obssim20240101_v013.arf...

Let’s view what the sources look like:

fitter.display_sources(data_pixel_size=8);

It appears that the blurred source does indeed match the data well. We can also view what the fitter thinks the data should look like, based on the starting parameters

fitter.plot();

Finally, we perform the fit

result = fitter.fit()
result

FitResult:
        q (src-ext) = 0.4289 +/- 0.0751
        u (src-ext) = 0.0077 +/- 0.0751
        q (bkg) = -0.0202 +/- 0.0462
        u (bkg) = -0.0532 +/- 0.0467
        f (bkg) = 2.1462 +/- 0.0448

Polarization:
        PD (src-ext): 0.4289 +/- 0.0751 (5.7 sig)
        PA (src-ext): 0.5144 deg +/- 5.0130
        PD (bkg): 0.0569 +/- 0.0466 (1.2 sig)
        PA (bkg): -55.3973 deg +/- 23.2842
Likelihood 11694.843175584967, dof 15655
Optimization terminated successfully.

The true source polarization was q=0.5, u=0, which the fit agrees with.