This code will search for a period in photon arrival times (as opposed to "countrate as a function of time" = "lightcurve") data in the context of X-ray and gamma-ray astronomy. This is a C implementation of the Hm test (de Jager, Raubenheimer & Swanepoel, 1989). The chance occurrence probability for a periodogram peak is calculated as Prob(> H) = exp (−0.4H) following equation (5) of de Jager & Büsching (2010) and corrected for multiple trials using an estimate of the number of independent frequencies as suggested by Schwarzenberg-Czerny (2003). While writing the code I was following the explanation of the Hm test by Kerr (2011). The code makes use of multiple processing cores thanks to OpenMP.
A detailed description of patpc is given in the Appendix by Sokolovsky et al. (2022).
In order to compile and run the code you'll need:
- A C compiler (tested with
gcc) GNU make(you'll need to install it if you are on Mac or something)- CFITSIO library (optional, needed to read OGIP FITS event files)
- GNU Scientific Library (optional)
- gnuplot and ImageMagick
convertfor plotting (optional)
Download and unpack the archive or clone this repository, then run make.
Specify the OGIP FITS event file or a single-column ASCII file listing the photon arrival times in seconds as the command line argument:
./patpc event_file.evt
or
./patpc photon_arrival_times.txt
This will print the parameters of the highest peak of the Hm periodogram and the power spectrum to the terminal. The code will also try to create plots of the power spectrum, Hm periodogram, binned lightcurve as a function of time, binned lightcurve as a function of phase (folded with the period corresponding to the highest Hm periodogram peak).
You may modify the plots by editing plot_everything.gnuplot and re-running gnuplot plot_everything.gnuplot.
You may also manually specify the search parameters:
./patpc event_file.evt 5000 5 0.1
where 5000 is the maximum trial period in seconds, 5 is the minimum trial period in seconds, 0.1 is the maximum phase shift between the first and the last point of the lightcurve (determines the width of the step between the trial frequencies).
Comparison of the power spectra computed with patpc, XRONOS/powspec
and z2n-periodogram codes from the same event file.
The slight difference with XRONOS/powspec
is because it computes a different thing: rather than follow the Discrete
Fourier Transform definition directly (Deeming, 1975)
as patpc and z2n-periodogram do,
powspec produces a binned ``lightcurve'' with a regular sampling where most of
the measurements are 0 (no photons arrived),
some measurements are equal to 1 (one photon arrived) and maybe some bins have more than one photon.
Then powspec runs an FFT on this lightcurve. Overall, this procedure is
similar, but not the same as computing the Discrete Fourier Transform over
the list of photon arrival times.
Bug reports, pull requests and feature suggestions are warmly welcome!