Note
Click here to download the full example code
Targeted grid search with line-robust BSGL statistic
Search for a monochromatic (no spindown) signal using a parameter space grid (i.e. no MCMC) and the line-robust BSGL statistic to distinguish an astrophysical signal from an artifact in a single detector.
11 import pyfstat
12 import numpy as np
13 import os
14
15 label = "PyFstat_example_grid_search_BSGL"
16 outdir = os.path.join("PyFstat_example_data", label)
17
18 F0 = 30.0
19 F1 = 0
20 F2 = 0
21 Alpha = 1.0
22 Delta = 1.5
23
24 # Properties of the GW data - first we make data for two detectors,
25 # both including Gaussian noise and a coherent 'astrophysical' signal.
26 depth = 70
27 sqrtS = "1e-23"
28 h0 = float(sqrtS) / depth
29 cosi = 0
30 IFOs = "H1,L1"
31 sqrtSX = ",".join(np.repeat(sqrtS, len(IFOs.split(","))))
32 tstart = 1000000000
33 duration = 100 * 86400
34 tend = tstart + duration
35 tref = 0.5 * (tstart + tend)
36
37 data = pyfstat.Writer(
38 label=label,
39 outdir=outdir,
40 tref=tref,
41 tstart=tstart,
42 duration=duration,
43 F0=F0,
44 F1=F1,
45 F2=F2,
46 Alpha=Alpha,
47 Delta=Delta,
48 h0=h0,
49 cosi=cosi,
50 sqrtSX=sqrtSX,
51 detectors=IFOs,
52 SFTWindowType="tukey",
53 SFTWindowBeta=0.001,
54 Band=1,
55 )
56 data.make_data()
57
58 # Now we add an additional single-detector artifact to H1 only.
59 # For simplicity, this is modelled here as a fully modulated CW-like signal,
60 # just restricted to the single detector.
61 SFTs_H1 = data.sftfilepath.split(";")[0]
62 extra_writer = pyfstat.Writer(
63 label=label,
64 outdir=outdir,
65 tref=tref,
66 F0=F0 + 0.01,
67 F1=F1,
68 F2=F2,
69 Alpha=Alpha,
70 Delta=Delta,
71 h0=10 * h0,
72 cosi=cosi,
73 sqrtSX=0, # don't add yet another set of Gaussian noise
74 noiseSFTs=SFTs_H1,
75 SFTWindowType="tukey",
76 SFTWindowBeta=0.001,
77 )
78 extra_writer.make_data()
79
80 # set up search parameter ranges
81 dF0 = 0.0001
82 DeltaF0 = 1000 * dF0
83 F0s = [F0 - DeltaF0 / 2.0, F0 + DeltaF0 / 2.0, dF0]
84 F1s = [F1]
85 F2s = [F2]
86 Alphas = [Alpha]
87 Deltas = [Delta]
88
89 # first search: standard F-statistic
90 # This should show a weak peak from the coherent signal
91 # and a larger one from the "line artifact" at higher frequency.
92 searchF = pyfstat.GridSearch(
93 label + "_twoF",
94 outdir,
95 os.path.join(outdir, "*" + label + "*sft"),
96 F0s,
97 F1s,
98 F2s,
99 Alphas,
100 Deltas,
101 tref,
102 tstart,
103 tend,
104 )
105 searchF.run()
106
107 print("Plotting 2F(F0)...")
108 searchF.plot_1D(xkey="F0")
109
110 # second search: line-robust statistic BSGL activated
111 searchBSGL = pyfstat.GridSearch(
112 label + "_BSGL",
113 outdir,
114 os.path.join(outdir, "*" + label + "*sft"),
115 F0s,
116 F1s,
117 F2s,
118 Alphas,
119 Deltas,
120 tref,
121 tstart,
122 tend,
123 BSGL=True,
124 )
125 searchBSGL.run()
126
127 # The actual output statistic is log10BSGL.
128 # The peak at the higher frequency from the "line artifact" should now
129 # be massively suppressed.
130 print("Plotting log10BSGL(F0)...")
131 searchBSGL.plot_1D(xkey="F0")
Total running time of the script: ( 0 minutes 0.000 seconds)