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