Fourier Amplitude Difference correction in Stingray

[1]:
import numpy as np
import matplotlib.pyplot as plt
from stingray.lightcurve import Lightcurve
from stingray.deadtime.fad import calculate_FAD_correction
from stingray.deadtime.filters import filter_for_deadtime
from stingray.crossspectrum import AveragedCrossspectrum
from stingray.powerspectrum import AveragedPowerspectrum

import matplotlib.pyplot as plt

Dead time affects most counting experiments. While the instrument is busy processing one event, it is “dead” to other photons/particles hitting the detector. This is usually not an issue if the count rate is low enough, or the processing time (dead time) is small enough. However, at high count rate dead time affects greatly the statistical properties of the data, to a point where a standard periodicity search based on the periodogram/power density spectrum (PDS) cannot be carried out.

The Fourier Amplitude Difference (FAD) correction is described in Bachetti & Huppenkothen, 2018, ApJ, 853L, 21, and is able to correct precisely deadtime affected PDSs if we have at least two identical and independent detectors. This is common in new generation X-ray timing instruments, often based on multiple-detector configurations (e.g. NuSTAR, NICER, AstroSAT, etc.).

In the code below, we calculate the PDS of light curves without dead time, after applying a dead time filter, and after applying the FAD to the dead-time affected dataset.

[2]:
def generate_events(length, ncounts):
    ev = np.random.uniform(0, length, ncounts)
    ev.sort()
    return ev


def generate_deadtime_lc(ev, dt, tstart=0, tseg=None, deadtime=2.5e-3):
    ev = filter_for_deadtime(ev, deadtime)
    return Lightcurve.make_lightcurve(ev, dt=dt, tstart=tstart, tseg=tseg, gti=np.array([[tstart, tseg]]))

[3]:
ctrate = 500
dt = 0.001
deadtime = 2.5e-3
tstart = 0
length = 25600
segment_size = 256.
ncounts = np.int(ctrate * length)
ev1 = generate_events(length, ncounts)
ev2 = generate_events(length, ncounts)
lc1 = Lightcurve.make_lightcurve(ev1, dt=dt, tstart=0, tseg=length,
                                 gti=np.array([[tstart, length]]))
lc2 = Lightcurve.make_lightcurve(ev2, dt=dt, tstart=0, tseg=length,
                                 gti=np.array([[tstart, length]]))

pds1 = AveragedPowerspectrum(lc1, segment_size=segment_size, norm='leahy')
pds2 = AveragedPowerspectrum(lc2, segment_size=segment_size, norm='leahy')
ptot = AveragedPowerspectrum(lc1 + lc2, segment_size=segment_size, norm='leahy')
cs = AveragedCrossspectrum(lc1, lc2, segment_size=segment_size, norm='leahy')
/home/mbachett/devel/stingray/build/lib/stingray/utils.py:102: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)
/home/mbachett/devel/stingray/build/lib/stingray/utils.py:102: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)
/home/mbachett/devel/stingray/build/lib/stingray/utils.py:102: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)
[4]:
lc1_dt = generate_deadtime_lc(ev1, dt, tstart=0, tseg=length, deadtime=deadtime)
lc2_dt = generate_deadtime_lc(ev2, dt, tstart=0, tseg=length, deadtime=deadtime)

pds1_dt = AveragedPowerspectrum(lc1_dt, segment_size=segment_size, norm='leahy')
pds2_dt = AveragedPowerspectrum(lc2_dt, segment_size=segment_size, norm='leahy')
ptot_dt = AveragedPowerspectrum(lc1_dt + lc2_dt, segment_size=segment_size, norm='leahy')
cs_dt = AveragedCrossspectrum(lc1_dt, lc2_dt, segment_size=segment_size, norm='leahy')
INFO:astropy:filter_for_deadtime: 7111455/12800000 events rejected
INFO: filter_for_deadtime: 7111455/12800000 events rejected [stingray.deadtime.filters]
INFO:astropy:filter_for_deadtime: 7109916/12800000 events rejected
INFO: filter_for_deadtime: 7109916/12800000 events rejected [stingray.deadtime.filters]
/home/mbachett/devel/stingray/build/lib/stingray/utils.py:102: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)
/home/mbachett/devel/stingray/build/lib/stingray/utils.py:102: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)
/home/mbachett/devel/stingray/build/lib/stingray/utils.py:102: UserWarning: SIMON says: Errorbars on cross spectra are not thoroughly tested. Please report any inconsistencies.
  warnings.warn("SIMON says: {0}".format(message), **kwargs)
[5]:
results = \
    calculate_FAD_correction(lc1_dt, lc2_dt, segment_size, plot=False,
                      smoothing_alg='gauss',
                      smoothing_length=segment_size*2,
                      strict=True, verbose=False,
                      tolerance=0.05, all_leahy=True)

freq_f = results['freq']
pds1_f = results['pds1']
pds2_f = results['pds2']
cs_f = results['cs']
ptot_f = results['ptot']

[6]:
for spec, spec_dt, spec_f, label in zip(
        [pds1, pds1, ptot, cs],
        [pds1_dt, pds2_dt, ptot_dt, cs_dt],
        [pds1_f, pds2_f, ptot_f, cs_f],
        ['PDS from light curve 1', 'PDS from light curve 2', 'PDS from lcs 1+2', 'cospectrum']
        ):
    plt.figure(figsize=(10, 8))
    plt.title(label)
    plt.plot(spec.freq, spec.power, label='No dead time', alpha=0.5)
    plt.plot(spec_dt.freq, spec_dt.power, label='Dead time-affected', alpha=0.5)
    plt.plot(freq_f, spec_f, label='FAD-corrected', alpha=0.5)
    plt.legend()
/home/mbachett/virtualenvs/py37/lib/python3.7/site-packages/IPython/core/events.py:88: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  func(*args, **kwargs)
/home/mbachett/virtualenvs/py37/lib/python3.7/site-packages/IPython/core/pylabtools.py:128: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
../../_images/notebooks_Deadtime_Check_FAD_correction_in_Stingray_7_1.png
../../_images/notebooks_Deadtime_Check_FAD_correction_in_Stingray_7_2.png
../../_images/notebooks_Deadtime_Check_FAD_correction_in_Stingray_7_3.png
../../_images/notebooks_Deadtime_Check_FAD_correction_in_Stingray_7_4.png

As can be seen above, all power density and co- spectra have been corrected accurately in their basic property (the white noise level). See Bachetti & Huppenkothen 2019 for more information.

[ ]: