This document is a step-by-step to repeat the analysis in Huppenkothen, Watts & Levin (submitted). The steps below require the basic data from RXTE and RHESSI (pre-processed) as well as some intermediate data products to avoid long runtimes (where long may mean weeks).
The code runs on my system (MacOS 10.9.4, python 2.7.8, ipython 2.7.8, numpy 1.8.1, scipy 0.14.0), but I haven't tested it on any others. Feel free to drop me a note if you find it works/doesn't work for you, or file an issue on github (preferred), or issue a pull request (even more preferred!).
The code assumes the basic structure of the repository: code is in giantflare-paper/code, all data products are in giantflare-paper/data, and documents and plots are in giantflare-paper/documents. Make sure to copy all data files that are not part of the repository itself and that you do not want to run yourself into giantflare-paper/data before you start!
The code is best run from within the directory giantflare/documents, but running it from other directories is possible. Note that in this code, plots will also be saved in this directory, because I didn't want to hard-code the directory structure into the code.
To load either the RXTE or RHESSI data set, we can perform the following command:
import sys
sys.path.insert(0, '../code/') ## add file path to PYTHONPATH
import giantflare
import giantflare_paper_analysis
print giantflare_paper_analysis.__file__ ## check that it uses the right script
#tnew_rxte = giantflare_paper_analysis.load_rxte_data(datadir="../data/", tstart=198.0)
#tnew_rhessi = giantflare_paper_analysis.load_rhessi_data(datadir="../data/")
#print(tnew_rxte)
#print(tnew_rhessi)
../code/giantflare_paper_analysis.py
/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/__init__.py:1155: UserWarning: This call to matplotlib.use() has no effect because the backend has already been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot, or matplotlib.backends is imported for the first time. warnings.warn(_use_error_msg)
This automatically loads the data in the right energy channels and in the right time interval for the analysis of the paper. The start time can be changed using the keyword tstart
, the range of energy channels using the keyword climits
followed by a list of two elements, the lower and upper channel.
One can compute the p-values for the RXTE data across a number of averaged cycles simply using the function rxte_pvalues
in giantflare_paper_analysis.py
. The code will load the RXTE data file and perform a similar analysis to Strohmayer & Watts (2006), but taking into account both periodograms in single cycles as well as periodograms averaged at the same rotational phase of the neutron star over consecutive rotational cycles (see Huppenkothen et al (submitted) for details).
You can invoke the command from inside ipython with
allstack, pvals, pvals_err = giantflare_paper_analysis.rxte_pvalues()
Loading data file ../data/1806.dat. Loading data file ../data/1806.dat. I am here in search_singlepulse Loading simulation data files in ../data/sgr1806_rxte_simulated_maxpowers.txt
--------------------------------------------------------------------------- IOError Traceback (most recent call last) <ipython-input-2-d6e76047b73e> in <module>() ----> 1 allstack, pvals, pvals_err = giantflare_paper_analysis.rxte_pvalues() /Users/danielahuppenkothen/repositories/giantflare-paper/code/giantflare_paper_analysis.py in rxte_pvalues(datadir, maxpowfile, savgallroot) 110 if not maxpowfile is None: 111 print("Loading simulation data files in %s"%maxpowfile) --> 112 maxp_all = np.loadtxt(maxpowfile) 113 114 /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/lib/npyio.pyc in loadtxt(fname, dtype, comments, delimiter, converters, skiprows, usecols, unpack, ndmin) 732 fh = iter(bz2.BZ2File(fname)) 733 elif sys.version_info[0] == 2: --> 734 fh = iter(open(fname, 'U')) 735 else: 736 fh = iter(open(fname)) IOError: [Errno 2] No such file or directory: '../data/sgr1806_rxte_simulated_maxpowers.txt'
This assumes that a file sgr1806_rxte_simulated_maxpowers.txt
is present in ../data
/. This file contains the maximum powers at 625 Hz across all segments for $10^5$ single-cycle and cycle-averaged simulations. Files with data of this type can be specified with the maxpowfile
keyword. Alternatively, the output of the function make_rxte_sims
produces a file with these simulations and can be passed into the savgallfile
keyword. Here, I've provided the file sgr1806_rxte_simulated_maxpowers.txt
only, contact me if you want the raw output of make_rxte_sims
.
To run your own simulations, invoke make_rxte_sims
like this:
nsims=3 ## change the number of simulations, small here not to make the program run forever
savgall = giantflare_paper_analysis.make_rxte_sims(tnew=None, nsims=nsims,save=True, froot="../data/sgr1806_rxte")
maxpowers = giantflare.compute_maxpowers(savgall, 10, 10, froot="../data/sgr1806_rxte")
Loading data file ../data/1806.dat. I am on simulation 0 I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse Saving RXTE simulation results in ../data/sgr1806_rxte_savgall.dat Saving maximum powers across all segments in ../data/sgr1806_rxte_simulated_maxpowers.txt
If you are in an interactive ipython session (like I am most times I do this), you probably have the RXTE arrival times loaded in tnew_rxte
or a similar name. In order to save a little on execution time, you can pass it into the tnew
keyword.
The resulting file can be passed to the savgallfile
keyword of rxte_pvalues
. For running many simulations on a multi-core processor, it may be useful to break the run up into smaller segments and start them independently. Threading is not supported. Sorry.
In order to compare the observed QPO with simulations that include a signal, we can make simulated light curves by smoothing out all variability >100 Hz, and then injecting back a signal. In the paper, we inject a periodic signal into a short (0.5s) segment into the cycle where we observe the highest power in the data, and the cycle exactly before that. This can be done in the function rxte_qpo_sims_singlecycle
, which takes the number of simulations as input (parameter nsims
):
nsims=3 ## change the number of simulations, small here not to make the program run forever
mid_coarse_all, savg_coarse_all, mid_fine_all, savg_fine_all, pvals_all = \
giantflare_paper_analysis.rxte_qpo_sims_singlecycle(nsims)
Loading data file ../data/1806.dat. min(counts_qpo): 0.000000 min(counts_qpo): 0.000000 min(counts_qpo): 0.000000 Loading simulated maximum powers across all segments, ../data/sgr1806_rxte_simulated_maxpowers.txt. I am on simulation 0 I am here in search_singlepulse I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse q_sims: [ 0.15 1.5 2.85] Saving data output for cycle with the strongest QPO signal in a bunch of filesof type ../data/sgr1806_rxte_strongestcycle* Saving figure in f6.pdf
The result is Figure 6 in the paper. This basically concludes the RXTE analysis as performed in the paper. Of course, you can fiddle with the details and parameters to your heart's content.
The analysis of the RHESSI data proceeds in much the same vein as for the RXTE data, but is a little more involved because of the lower data quality of the RHESSI data. For the RXTE data, it is fairly obvious that the signal is confined to a single bin.
For the RHESSI data, we perform a somewhat more complex analysis, where we search segments of different size for signals.
This requires four separate simulation runs, one for each of the tested segment sizes (bins = [0.5, 1.0, 2.0, 2.5]
). Note that for the 2-second-segments, this run is very, very big ($10^6$ simulations; $10^4$ for the other segments), which took me about 3 weeks on a 20-core machine. If you don't want to run this yourself, I suggest using the files I made.
Below is nevertheless the code that makes these simulations:
tnew_rhessi = giantflare_paper_analysis.load_rhessi_data(datadir="../data/") ## load RHESSI data
nsims = 5 ## change number of simulations, small here not to make the program run forever
savgall_05 = giantflare_paper_analysis.make_rhessi_sims(tnew=tnew_rhessi, tseg_all=[0.5], nsims=nsims,save=True, froot="../data/sgr1806_rhessi_tseg=0.5s")
savgall_1 = giantflare_paper_analysis.make_rhessi_sims(tnew=tnew_rhessi, tseg_all=[1.0], nsims=nsims,save=True, froot="../data/sgr1806_rhessi_tseg=1.0s")
savgall_2 = giantflare_paper_analysis.make_rhessi_sims(tnew=tnew_rhessi, tseg_all=[2.0], nsims=nsims,save=True, froot="../data/sgr1806_rhessi_tseg=2.0s")
savgall_25 = giantflare_paper_analysis.make_rhessi_sims(tnew=tnew_rhessi, tseg_all=[2.5], nsims=nsims,save=True, froot="../data/sgr1806_rhessi_tseg=2.5s")
Loading RHESSI data from ../data/1806_rhessi.dat I am on simulation 0 I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse I am on simulation 3 I am here in search_singlepulse I am on simulation 4 I am here in search_singlepulse Saving output of simulations in ../data/sgr1806_rhessi_tseg=0.5s_tseg=0.50_df=2.00_savgall.txt I am on simulation 0 I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse I am on simulation 3 I am here in search_singlepulse I am on simulation 4 I am here in search_singlepulse Saving output of simulations in ../data/sgr1806_rhessi_tseg=1.0s_tseg=1.00_df=1.00_savgall.txt I am on simulation 0 I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse I am on simulation 3 I am here in search_singlepulse I am on simulation 4 I am here in search_singlepulse Saving output of simulations in ../data/sgr1806_rhessi_tseg=2.0s_tseg=2.00_df=0.50_savgall.txt I am on simulation 0 I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse I am on simulation 3 I am here in search_singlepulse I am on simulation 4 I am here in search_singlepulse Saving output of simulations in ../data/sgr1806_rhessi_tseg=2.5s_tseg=2.50_df=0.40_savgall.txt
This will create four numpy arrays with the power at 625 Hz in every segment. The array is also saved in a simple ascii file of filename "%s_savgall.dat"%froot
. Once you have run simulations for all your required segments (the ones in the example above are the default values I used, but you can fiddle with them as necessary), they can then be directly used to compute the p-values for each segment, and make Figure 7 of the paper:
pvals_all = giantflare_paper_analysis.rhessi_simulations_results(tnew=None, tseg_all=[0.5, 1.0, 2.0, 2.5], df_all=[2.0, 1.0, 1.0, 1.0],
freq = [627.0, 627.0, 626.5,626.0], classical=False,maxpowers=True,
froot_in="../data/sgr1806_rhessi", froot_out="../data/sgr1806_rhessi", plotdist=False)
Loading RHESSI data from ../data/1806_rhessi.dat I am here in search_singlepulse index: 0, f: 623.624, P: 15.062 index: 1, f: 625.626, P: 12.986 index: 2, f: 627.628, P: 18.808 index: 3, f: 629.630, P: 13.041 index: 4, f: 631.632, P: 22.583 ../data/sgr1806_rhessi*_tseg=0.5*maxpowers.txt ['../data/sgr1806_rhessi_tseg=0.5_simulated_maxpowers.txt'] sims: 10000 ind_sims: 9608 sims: 10000 ind_sims: 9597 sims: 10000 ind_sims: 8542 sims: 10000 ind_sims: 7194 sims: 10000 ind_sims: 7925 sims: 10000 ind_sims: 7729 sims: 10000 ind_sims: 8912 sims: 10000 ind_sims: 7980 sims: 10000 ind_sims: 8268 sims: 10000 ind_sims: 8329 sims: 10000 ind_sims: 8202 sims: 10000 ind_sims: 8319 sims: 10000 ind_sims: 8583 sims: 10000 ind_sims: 8709 sims: 10000 ind_sims: 9079 sims: 10000 ind_sims: 9581 sims: 10000 ind_sims: 9733 sims: 10000 ind_sims: 9763 sims: 10000 ind_sims: 9784 I am here in search_singlepulse index: 0, f: 625.813, P: 11.825 index: 1, f: 626.813, P: 12.182 index: 2, f: 627.814, P: 18.729 index: 3, f: 628.814, P: 16.987 index: 4, f: 629.815, P: 13.398 ../data/sgr1806_rhessi*_tseg=1.0*maxpowers.txt ['../data/sgr1806_rhessi_tseg=1.0_simulated_maxpowers.txt'] sims: 10000 ind_sims: 9616 sims: 10000 ind_sims: 9600 sims: 10000 ind_sims: 9374 sims: 10000 ind_sims: 8789 sims: 10000 ind_sims: 8733 sims: 10000 ind_sims: 9481 sims: 10000 ind_sims: 9267 sims: 10000 ind_sims: 8914 sims: 10000 ind_sims: 8456 sims: 10000 ind_sims: 8548 sims: 10000 ind_sims: 8116 sims: 10000 ind_sims: 6952 sims: 10000 ind_sims: 7228 sims: 10000 ind_sims: 7090 sims: 10000 ind_sims: 5841 sims: 10000 ind_sims: 4793 sims: 10000 ind_sims: 6084 sims: 10000 ind_sims: 5366 sims: 10000 ind_sims: 6860 I am here in search_singlepulse index: 0, f: 625.656, P: 12.789 index: 1, f: 626.157, P: 14.604 index: 2, f: 626.657, P: 21.511 index: 3, f: 627.157, P: 15.516 index: 4, f: 627.657, P: 20.697 ../data/sgr1806_rhessi*_tseg=2.0*maxpowers.txt ['../data/sgr1806_rhessi_tseg=2.0_simulated_maxpowers.txt'] sims: 1075000 ind_sims: 1074989 sims: 1075000 ind_sims: 1074992 sims: 1075000 ind_sims: 1075000 sims: 1075000 ind_sims: 1075000 sims: 1075000 ind_sims: 1074999 sims: 1075000 ind_sims: 1074998 sims: 1075000 ind_sims: 1075000 sims: 1075000 ind_sims: 1074999 sims: 1075000 ind_sims: 1074996 sims: 1075000 ind_sims: 1074982 sims: 1075000 ind_sims: 1074965 sims: 1075000 ind_sims: 1074922 sims: 1075000 ind_sims: 1074997 sims: 1075000 ind_sims: 1074999 sims: 1075000 ind_sims: 1075000 sims: 1075000 ind_sims: 1075000 sims: 1075000 ind_sims: 1075000 sims: 1075000 ind_sims: 1075000 sims: 1075000 ind_sims: 1075000 I am here in search_singlepulse index: 0, f: 624.725, P: 9.921 index: 1, f: 625.525, P: 8.835 index: 2, f: 626.325, P: 10.790 index: 3, f: 627.125, P: 9.189 index: 4, f: 627.926, P: 8.936 ../data/sgr1806_rhessi*_tseg=2.5*maxpowers.txt ['../data/sgr1806_rhessi_tseg=2.5_simulated_maxpowers.txt'] sims: 10600 ind_sims: 9515 sims: 10600 ind_sims: 7691 sims: 10600 ind_sims: 8236 sims: 10600 ind_sims: 8973 sims: 10600 ind_sims: 9603 sims: 10600 ind_sims: 9300 sims: 10600 ind_sims: 9825 sims: 10600 ind_sims: 8493 sims: 10600 ind_sims: 8897 sims: 10600 ind_sims: 7959 sims: 10600 ind_sims: 8332 sims: 10600 ind_sims: 8839 sims: 10600 ind_sims: 8872 sims: 10600 ind_sims: 10027 sims: 10600 ind_sims: 10169 sims: 10600 ind_sims: 10139 sims: 10600 ind_sims: 10265 sims: 10600 ind_sims: 10186 sims: 10600 ind_sims: 10294 Saving p-values in Figure f7.pdf
If you have the list of photon arrival times loaded, you can put it directly into the tnew
keyword and save a bit of execution time. tseg_all
needs to be a list of all time segments you are interested in (and have savgall-files for!). In df_all
, you can set a frequency resolution, but be warned that the new frequency resolution must be an integer multiple of the original frequency resolution! freq
defines the frequency of the QPO, i.e. the frequency we are interested in. Note that this changes for different segments here, because the frequency resolution is different, shifting the bulk of the signal into slightly different bins.
If classical
is True, then no simulations are used, and the probabilities for the maximum powers at the QPO frequency are computed using analytical equations from Groth (1975) appropriate for detection of periodic signals in white noise. This can serve as a first check of what we expect the significance to be. I mostly used it to roughly estimate the number of simulations required for each segment. DO NOT DO THIS FOR LOW FREQUENCIES! At low frequencies, the periodgrams are dominated by the underlying pulse profile, and classical p-values will vastly overstate the signficance of a signal!
If maxpowers=True
, the code assumes that you have a file with maximum powers across all segments for each segment size. If you need to re-run the function several times, using a file with the maximum powers at the QPO frequency for all segments across segment sizes and averaged cycles will save you a lot of time, because converting a savgall-file into a maxpowers file, while easily done, takes a long time to do (think a couple of hours for 10000 simulations), while loading the file from disc doesn't. So, run the first time with maxpowers=False
, if you have your own simulations (a maxpowers file is provided from my own runs), to have the code convert your simulation output files automatically in the right format for next time. For subsequent runs, set maxpowers=True
. The data products I supply on figshare are largely the maxpowers-files, because they have manageable size. Note that because the maxpowers-file for $t_{\mathrm{seg}}=2\,\mathrm{s}$ for the RHESSI data is too big for figshare to accept, I've split it in three (sgr1806_rhessi_tseg=2.0_simulated_maxpowers[1,2,3].txt). You'll have to stitch it together like this if you want to use it in its entirety:
files = ["sgr1806_rhessi_tseg=2.0_simulated_maxpowers1.txt", "sgr1806_rhessi_tseg=2.0_simulated_maxpowers2.txt",
"sgr1806_rhessi_tseg=2.0_simulated_maxpowers3.txt"]
maxp_all = []
for f in files:
data = np.atleast_2d(np.loadtxt(f))
maxp_all.extend(data)
savetxt("sgr1806_rhessi_tseg=2.0_simulated_maxpowers.txt", maxp_all)
In froot_in
, specify details of the input filename, either for your direct output of the simulations, or the processed maxpowers-data file. The filename needs to be of type "%s*_tseg=%.1f*savgall*"%(froot_in, tseg)
, where tseg
is an element of tseg_all
. The program will glob the file names in your directory. If you're like me, then you've probably run make_rhessi_sims
on a multi-core machine in several stages (I used increments of 5000, using 10000 and running several in parallel gets quite memory-intensive), so that I have files, for example, ["sgr1806_rhessi1_tseg=0.5_savgall.dat", "sgr1806_rhessi2_tseg=0.5_savgall.dat", ..., "sgr1806_rhessi25_tseg=0.5_savgall.dat"]
, which glob ought to find fairly easily. The code will then compute the averages over 2
to n
cycles (19 for the RHESSI data), find the maximum power across all segments, and save the resulting numpy array in a text-file of the type "%s_tseg=%.1f_simulated_maxpowers.txt"%(froot_out, tseg)
, which can be used in future runs.
If plotdist=True
, the code will produce a number of plots (one for each cycle average and each segment size), comparing the theoretical expected distribution of maximum powers across segments (derived from draws from a $\chi^2_{2n}$ distribution, where $n$ is the number of averaged cycles) with that derived from simulations of the light curve. This provides a test of how different the actual distributions are from the theoretical expectations.
In order to understand the behaviour of the QPO in the RHESSI data better, and to see whether the observations require a signal being present throughout all 19 cycles (as opposed to a short-lived, possibly re-excited signal), we performed simulations with a periodic signal injected into the smoothed light curve, as we did for the RXTE data. We made two different assumptions: (1) a QPO that is present at a fairly low level in all cycles, (2) a signal that is strongly present in 2 cycles, and re-excited a few cycles later, but at a lower level. The code also allows for a third possibility not included in the paper: a periodic signal that is re-excited in every cycle, but with an onset that varies randomly in rotational phase (of the neutron star).
In order to test the two different hypothesis, we can create simulations much like in the function make_rhessi_sims
, but with an additional periodic signal. Each of these simulations is then analysed in much the same way as the observed data, by comparing them to our initial run of simulations without periodic signal, for all segment sizes. Note that this function requires a maxpower-file for each segment size of interest, otherwise it'll never finish.
You can run it like this:
nsims=3 ## change number of simulations, small here not to make the program run forever
giantflare_paper_analysis.rhessi_qpo_sims_allcycles(nsims=nsims, froot="../data/sgr1806_rhessi")
giantflare_paper_analysis.rhessi_qpo_sims_singlecycle(nsims=nsims, froot="../data/sgr1806_rhessi")
## giantflare_paper_analysis.rhessi_qpo_sims_allcycles_randomised(nsims=nsims, froot="../data/sgr1806_rhessi")
Loading RHESSI data from ../data/1806_rhessi.dat min(counts_qpo): 0.000000 min(counts_qpo): 0.000000 min(counts_qpo): 0.000000 I am on simulation 0 I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse Saving data in ../data/sgr1806_rhessi_allcycles_tseg=0.5_savgall.txt Saving data in ../data/sgr1806_rhessi_allcycles_tseg=1.0_savgall.txt Saving data in ../data/sgr1806_rhessi_allcycles_tseg=2.0_savgall.txt Saving data in ../data/sgr1806_rhessi_allcycles_tseg=2.5_savgall.txt Loading RHESSI data from ../data/1806_rhessi.dat min(counts_qpo): 0.000000 min(counts_qpo): 0.000000 min(counts_qpo): 0.000000 I am on simulation 0 I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am on simulation 1 I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am on simulation 2 I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse I am here in search_singlepulse Saving data in ../data/sgr1806_rhessi_singlecycle_tseg=0.5_savgall.txt Saving data in ../data/sgr1806_rhessi_singlecycle_tseg=1.0_savgall.txt Saving data in ../data/sgr1806_rhessi_singlecycle_tseg=2.0_savgall.txt Saving data in ../data/sgr1806_rhessi_singlecycle_tseg=2.5_savgall.txt
Or you can use the files I provide. Note that the second (commented out) command is for case (3) mentioned above that is not part of the paper, but one of the possibilities I played around with.
I hard-coded the details of the periodic signal (where it appears in each cycle, in how many cycles it appears, at which frequency and amplitude) into these three functions, but of course, you can play around with the code as you wish. Note that the name rhessi_qpo_sims_singlecycle
is kind of misleading, because it actually contains a signal that's re-excited several times (but is not present in every cycle).
If you have run these simulations, or decided to use my runs, you can evaluate them like this:
pvals_all, pvals_hist_all = giantflare_paper_analysis.rhessi_qpo_sims_images(tseg_all=[0.5,1.0,2.0,2.5],
df_all=[2.0, 1.0, 1.0, 1.0],
nbins=30,
froot_in="../data/sgr1806_rhessi",
froot_sims="allcycles",
classical=False,
freq=[627.0, 627.0, 626.5,626.0], plotn=8)
Loading RHESSI data from ../data/1806_rhessi.dat froot_in: ../data/sgr1806_rhessi I am here in search_singlepulse qpofiles: ['../data/sgr1806_rhessi2_allcycles_tseg=0.5_savgall.txt', '../data/sgr1806_rhessi3_allcycles_tseg=0.5_savgall.txt', '../data/sgr1806_rhessi_allcycles_tseg=0.5_savgall.txt'] ../data/sgr1806_rhessi*_allcycles_tseg=0.5*savgall.txt len(pvals_data): 19 log10(pvals_data): [-1.40671393 -1.39469495 -0.83624248 -0.55191233 -0.6829819 -0.64378287 -0.9633711 -0.69464863 -0.76145211 -0.77702355 -0.74521031 -0.77443229 -0.84863015 -0.88907376 -1.03574037 -1.37778598 -1.57348874 -1.62525165 -1.66554625] pd: 0.0392 pd: 0.0403 pd: 0.1458 pd: 0.2806 pd: 0.2075 pd: 0.2271 pd: 0.1088 pd: 0.202 pd: 0.1732 pd: 0.1671 pd: 0.1798 pd: 0.1681 pd: 0.1417 pd: 0.1291 pd: 0.0921 pd: 0.0419 pd: 0.0267 pd: 0.0237 pd: 0.0216 I am here in search_singlepulse qpofiles: ['../data/sgr1806_rhessi2_allcycles_tseg=1.0_savgall.txt', '../data/sgr1806_rhessi3_allcycles_tseg=1.0_savgall.txt', '../data/sgr1806_rhessi_allcycles_tseg=1.0_savgall.txt'] ../data/sgr1806_rhessi*_allcycles_tseg=1.0*savgall.txt len(pvals_data): 19 log10(pvals_data): [-1.41566878 -1.39794001 -1.20342567 -0.91685586 -0.89722339 -1.28483264 -1.13489603 -0.96417017 -0.8113527 -0.83803338 -0.7249191 -0.51598504 -0.55720677 -0.53610701 -0.38101108 -0.28341242 -0.40715732 -0.33404397 -0.50307035] pd: 0.0384 pd: 0.04 pd: 0.0626 pd: 0.1211 pd: 0.1267 pd: 0.0519 pd: 0.0733 pd: 0.1086 pd: 0.1544 pd: 0.1452 pd: 0.1884 pd: 0.3048 pd: 0.2772 pd: 0.291 pd: 0.4159 pd: 0.5207 pd: 0.3916 pd: 0.4634 pd: 0.314 I am here in search_singlepulse qpofiles: ['../data/sgr1806_rhessi2_allcycles_tseg=2.0_savgall.txt', '../data/sgr1806_rhessi3_allcycles_tseg=2.0_savgall.txt', '../data/sgr1806_rhessi_allcycles_tseg=2.0_savgall.txt'] ../data/sgr1806_rhessi*_allcycles_tseg=2.0*savgall.txt len(pvals_data): 19 log10(pvals_data): [-4.99001578 -5.12831848 -inf -inf -6.03140846 -5.73037847 -inf -6.03140846 -5.42934847 -4.77613596 -4.48734042 -4.13931386 -5.55428721 -6.03140846 -inf -inf -inf -inf -inf] pd: 1.02325581395e-05 pd: 7.44186046512e-06 pd: 0.0 pd: 0.0 pd: 9.3023255814e-07 pd: 1.86046511628e-06 pd: 0.0 pd: 9.3023255814e-07 pd: 3.72093023256e-06 pd: 1.67441860465e-05 pd: 3.25581395349e-05 pd: 7.25581395349e-05 pd: 2.79069767442e-06 pd: 9.3023255814e-07 pd: 0.0 pd: 0.0 pd: 0.0 pd: 0.0 pd: 0.0 I am here in search_singlepulse qpofiles: ['../data/sgr1806_rhessi2_allcycles_tseg=2.5_savgall.txt', '../data/sgr1806_rhessi3_allcycles_tseg=2.5_savgall.txt', '../data/sgr1806_rhessi_allcycles_tseg=2.5_savgall.txt'] ../data/sgr1806_rhessi*_allcycles_tseg=2.5*savgall.txt len(pvals_data): 19 log10(pvals_data): [-0.98987613 -0.56156214 -0.65165839 -0.81391831 -1.02661071 -0.91136251 -1.13600416 -0.70164133 -0.79409122 -0.60353746 -0.66966282 -0.77954651 -0.78776213 -1.26715124 -1.3908286 -1.36160494 -1.50026106 -1.40830552 -1.53958444] pd: 0.102358490566 pd: 0.274433962264 pd: 0.223018867925 pd: 0.153490566038 pd: 0.0940566037736 pd: 0.122641509434 pd: 0.0731132075472 pd: 0.198773584906 pd: 0.160660377358 pd: 0.249150943396 pd: 0.213962264151 pd: 0.166132075472 pd: 0.163018867925 pd: 0.0540566037736 pd: 0.0406603773585 pd: 0.0434905660377 pd: 0.0316037735849 pd: 0.0390566037736 pd: 0.0288679245283 Saving plot for p-values of simulations with periodic signal in allcycles in file f8.pdf
Here, tseg_all
and df_all
should be the same as for the observed data, otherwise the results will not be comparable!
The code will produce a histograms for the simulations with injected signal (for each cycle average) and plot these as heat map, with the observed p-values overplotted. The parameter nbins
sets the number of bins of the histogram.
froot_in
should match froot
from the function that made the simulations with QPO, froot_sims
is a keyword that's either "allcycles"
, "allcyclesrandomised"
or "singlecycle"
, depending on which function you used above. As before classical
sets a Boolean variable that decides whether a classical (analytical) significance test is used rather than simulations of the light curve, and the list freq
sets the frequency at which the code looks for the QPO.
In the case above, I run the code for simulations with a periodic signal in all cycles of the type "%s*_%s_tseg=%.1f*savgall.txt"%(froot_in,froot_sims, tseg)
, where the code again uses glob
to find all files matching this description (this means you can run these simulations in increments across multiple processors, too).
Doing the same for the other two types of simulations is straightforward.
## The first 4 plots are just illustrations of how the analysis works
#giantflare_paper_analysis.plot_averaging_example() ## plot Figure 1
#giantflare_paper_analysis.maxpower_plot() ## plot Figure 2
#giantflare_paper_analysis.make_analysis_schematic() ## plot Figure 3
#giantflare_paper_analysis.plot_lightcurves(datadir="../data/") ## plot Figure 4
## RXTE Analysis
#giantflare_paper_analysis.plot_rxte_pvalues(nsims=100000, froot = "../data/sgr1806_rxte") ## plot Figure 5
#giantflare_paper_analysis.plot_rxte_sims_singlecycle(froot="../data/sgr1806_rxte_strongestcycle") ## plot Figure 6
## RHESSI Analysis
giantflare_paper_analysis.plot_rhessi_pvalues(filename="../data/sgr1806_rhessi_pvals_all.txt", tseg=[0.5,1.0,2.0,2.5]) ## plot Figure 7
#giantflare_paper_analysis.plot_rhessi_qpo_simulations(tnew=None, tseg_all=[0.5,1.0,2.0,2.5], df_all = [2.0, 1.0, 1.0, 1.0, 1.0],
# freq=[627.0, 627.0, 626.5,626.0])
shape pvals_all: (5, 19)
../code/utils.py:114: RuntimeWarning: invalid value encountered in sqrt perr_val = np.sqrt(np.array(pval)*(1.0-np.array(pval))/np.float(n))
If really all you want to do is re-make the plots from the paper, and not in re-running the entire analysis, you can run the plotting functions at the bottom of the script only:
And that's it! You're done! :)