Fitting for variable sources

Weights used:

• Spatial

• Temporal

• Spectral

This simulated example source is pulsed and has a sweeping PA. This notebook determines the PD and EVPA of the sweep.

import numpy as np
import leakagelib

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

We’ll load the data as usual, and cut to within 280 arcsec of the center

datas = leakagelib.IXPEData.load_all_detectors_with_path("data", "pulse")
for data in datas:
    data.iterative_centroid_center()
    data.retain(data.evt_energies > 2)
    data.retain(data.evt_energies < 8)

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

We plan to implement phase weights, but LeakageLib records times, not phases. You should calculate them yourself (e.g. with PINT, or load phases from a separate file), then replace the evt_times field with these phases. For this mock data, the frequency is exactly 10 Hz, so the phases are just \(10 t\).

# WARNING: only run this box once.

for data in datas:
    # It's convenient to overwrite the "times" list with a list of phases. This source sweeps with frequency of 10 Hz, so multiplying by 10 gives the phase
    data.evt_times *= 10
    data.evt_times = np.fmod(data.evt_times, 1)

Now we create the point source and background sources as usual, setting their spectra. For this fit, we’ll fix the background to be unpolarized.

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

settings.add_point_source("src")
settings.fix_flux("src", 1)

settings.add_background("bkg")
settings.fix_qu("bkg", (0, 0))
settings.set_initial_flux("bkg", 1)

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

6462 events were cut for being outside the region of interest.
>>> 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...

To implement phase weights, we need to tell LeakageLib the source’s light curve. This source was simulated with a sine squared light curve. Note that normalization doesn’t matter.

Warning

If you cut based on phase e.g. for an on-off fit, you will need to supply the duty_cycle argument in set_lightcurve.

settings.set_lightcurve("src", lambda ph: np.sin(ph * np.pi)**2)

Now we can perform the fit and see what we get.

fitter = leakagelib.Fitter(settings)
print(fitter)
const_result = fitter.fit()
const_result

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

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

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

FitResult:
        q (src) = -0.0033 +/- 0.0721
        u (src) = -0.0027 +/- 0.0721
        f (bkg) = 2.1789 +/- 0.0417

Polarization:
        PD (src): 0.0043 +/- 0.0721
        PA (src): -70.5491 deg +/- 481.0701
Likelihood 19526.389495555064, dof 15927
Optimization terminated successfully.

The best-fit polarization degree was quite low! This is because the source’s true PA sweeps, and we modeled it as constant. We need to make a sweeping PA model.

Fitting a custom polarization model

We’ll test a model where the EVPA rotates by 2 pi every pulse, and PD is constant. This model will have two parameters: a constant PD and a EVPA at phase zero. Let’s create those parameters. The add_param function allows us to do this, and give initial fit values and bounds.

settings.add_param("sweep-PD", 0.1, [0, 1])
settings.add_param("sweep-PA", 0, [-np.pi, np.pi])

Now we need to create a function which takes in the phase and gives the expected Q and U polarization, and tell the fitter to use this function.

def model_fn(ph, fit_data, param_array):
    pd = fit_data.param_to_value(param_array, "sweep-PD")
    pa = fit_data.param_to_value(param_array, "sweep-PA")
    q = pd * np.cos(2 * (ph * 2*np.pi + pa))
    u = pd * np.sin(2 * (ph * 2*np.pi + pa))
    return q, u
settings.set_model_fn("src", model_fn)

Note that set_model_fn automatically tells the fitter to no longer use the src q and src u params. You can see that by printing the fitter object, which displays the free parameters.

fitter = leakagelib.Fitter(settings)
print(fitter)

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

FITTED PARAMETERS:
Source  Param
bkg:    f
None:   sweep-PD
None:   sweep-PA

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

Running the fit will only fit for f, sweep-PD, and sweep-PA.

sweep_result = fitter.fit()
sweep_result

FitResult:
        f (bkg) = 2.1787 +/- 0.0417
        sweep-PD = 0.4722 +/- 0.0711
        sweep-PA = 0.0244 +/- 0.0763

Polarization:
Likelihood 19548.00863653022, dof 15927
Optimization terminated successfully.

The PD of the sweeping fit was higher. The likelihood was also higher, indicating the sweep model provided a better fit. The true source polarization had PD = 0.5 and PA = 0, which the fit agrees with.