Fitting for a point source polarization with LeakageLib

Weights used:

• Spatial

• Spectral

This jupyter notebook fits for the polarization of a single simulated point source, with a photon and particle background. See other notebooks for how to use the more advanced weights.

import leakagelib

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

First, we need to load the data. See the next notebook for how to download and prepare real IXPE data. This notebook uses simulated data I made, which is located in <LEAKAGELIB>/docs/examples/data/ps. This notebook also comes downloaded with LeakageLib, located at <LEAKAGELIB>/docs/examples/fit-point-source.ipynb. If you run that notebook, the following load command should work.

datas = leakagelib.IXPEData.load_all_detectors_with_path("data", "ps")

>>> Reading (in memory) /opt/homebrew/anaconda3/lib/python3.12/site-packages/ixpeobssim/caldb/ixpe/xrt/bcf/vign/ixpe_d1_obssim20240101_vign_v013.fits...

Alternatives to IXPEData.load_all_detectors_with_path

• 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.

• 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.

• IXPEData. You must pass specific file names.

Next we center the data and cut to 2-8 keV.

for data in datas:
    data.iterative_centroid_center()
    data.retain(data.evt_energies > 2)
    data.retain(data.evt_energies < 8)

Alternatives to IXPEData.iterative_centroid_center

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

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

Now we have to tell the fitter how we want to do the fit. All this information is stored in the FitSettings object. First, we tell it to only consider events within 280 arcseconds of the center. More distant events will automatically be removed from the data sets.

settings = leakagelib.FitSettings(datas)
settings.apply_circular_roi(280)

6282 events were cut for being outside the region of interest.

Alternatives to FitSettings.apply_circular_roi

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

• IXPEData.apply_roi_cut: Adds a region file to the ROI (Or excludes it if you set exclude=True)

• 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.

You might be curious how many events are left. You can check this here:

sum([len(data) for data in datas])

15894

For this simplified scenario there are two components: a point source and a photon background. To create the point source,

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 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. These background sources are assumed to be spatially uniform and already have spatial weights applied.

settings.add_background() # Named "bkg" by default
settings.fix_qu("bkg", (0, 0)) # Assume an unpolarized background

To use energy weights, we must assign a spectrum. We create functions which take in the energy in keV and return the photon spectrum (dN/dE). Normalization doesn’t matter. The ARF and RMF are automatically included, but you can disable them with the use_arf and use_rmf arguments.

settings.set_spectrum("bkg", lambda e: e**-2.5) # Gamma=2.5, unabsorbed powerlaw
settings.set_spectrum("src", lambda 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...
>>> Using cached xEffectiveArea object at /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. Now we are ready to create the fitter…

fitter = leakagelib.Fitter(settings)

Data set ps DU 1 had no exposure map loaded. Please load an exposure map if you are fitting to events in the vignetted portion.

… and run the fit

result = fitter.fit()
result

FitResult:
        q (src) = 0.4697 +/- 0.0725
        u (src) = -0.0260 +/- 0.0721
        f (bkg) = 2.1813 +/- 0.0423

Polarization:
        PD (src): 0.4704 +/- 0.0725
        PA (src): -1.5836 deg +/- 4.3936
Likelihood 17886.60119278403, dof 15891
Optimization terminated successfully.

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

Alternatives to Fitter.fit

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

• IXPEData.fit_mcmc: Do a full MCMC fit.

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. The covariance matrix entries are stored in the same order as the result.parameter_names array.

print("===Results as human-readable dictionaries===")
print("Parameters", result.params)
print("Uncertainties", result.sigmas)
print()
print("===Results as arrays===")
print("Parameter names", result.parameter_names)
print("Best fit values", result.best_params)
print("Covariance", result.cov)

===Results as human-readable dictionaries===
Parameters {('q', 'src'): 0.46965223111815946, ('u', 'src'): -0.025987502310235725, ('f', 'bkg'): 2.1813398855067225}
Uncertainties {('q', 'src'): 0.07247836727152544, ('u', 'src'): 0.07213767036302343, ('f', 'bkg'): 0.0423476862090614}

===Results as arrays===
Parameter names [('q', 'src'), ('u', 'src'), ('f', 'bkg')]
Best fit values [ 0.46965223 -0.0259875   2.18133989]
Covariance [[ 5.25311372e-03 -1.73106799e-05  0.00000000e+00]
 [-1.73106799e-05  5.20384349e-03  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  1.79332653e-03]]