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