1
2
3 from glue.ligolw import ligolw, lsctables, utils
4 from pylal.tools import XLALCalculateEThincaParameter
5 from glue import segments, lal, git_version
6 from glue.ligolw.utils import process
7 import copy
8 import lal
9 import numpy
10 import math
11 import glob
12 import optparse
13 import sys
14
16
17
18
19 a11 = row.Gamma0 / eMatch;
20 a12 = row.Gamma1 / eMatch;
21 a13 = row.Gamma2 / eMatch;
22 a22 = row.Gamma3 / eMatch;
23 a23 = row.Gamma4 / eMatch;
24 a33 = row.Gamma5 / eMatch;
25
26 x = (a23 * a23 - a22 * a33) * a22;
27 denom = (a12*a23 - a22*a13) * (a12*a23 - a22*a13) - (a23*a23 - a22*a33) * (a12*a12 - a22*a11);
28
29 return math.sqrt( x / denom ) + 2. * lal.REARTH_SI / lal.C_SI;
30
31 -def find_slide_coincs(htrigs, ltrigs, min_mchirp, max_mchirp,
32 ethinca, slide, threshold):
33 """ Calculate coincs from single detector triggers. We window based on
34 mchirp, apply a new snr threshold, and only keep triggers that lie on a
35 slide boundary.
36 """
37
38 coinc_sngls = lsctables.New(lsctables.SnglInspiralTable)
39 num_trig = 0
40 t_snrsq = float(threshold) ** 2
41
42
43
44 l_ends = []
45 lt_mc = []
46 lt_sn = []
47 lt_en = []
48 for l in ltrigs:
49 lt_mc.append(l.mchirp)
50 lt_sn.append(l.get_new_snr()**2)
51 lt_en.append(float(l.get_end()))
52 lt_mc = numpy.array(lt_mc)
53 lt_sn = numpy.array(lt_sn)
54 lt_en = numpy.array(lt_en)
55
56
57
58
59 loc = numpy.arange(0, len(ltrigs))
60
61
62 stats = []
63 htrigs.sort(lambda a, b: cmp(a.get_new_snr(), b.get_new_snr()), reverse=True)
64
65
66
67 for h in htrigs:
68 h_end = float(h.get_end())
69 h_mchirp = h.mchirp
70 h_max_td = max_dt(h, ethinca)
71 h_snrsq = h.get_new_snr() ** 2
72 r_snrsq = t_snrsq - h_snrsq
73 lind = 0
74
75
76
77 lt_snl = lt_sn[loc]
78 lt_mcl = lt_mc[loc]
79 lt_enl = lt_en[loc]
80
81
82
83 num_slides = numpy.round((h_end - lt_enl)/slide)
84 lt_en_tmp = lt_enl + num_slides * slide
85 td = abs(h_end - lt_en_tmp)
86
87
88
89 goods = (lt_snl >= r_snrsq)
90
91
92
93 good = goods & (td <=h_max_td) & (lt_mcl >= 2*min_mchirp - h_mchirp) & (lt_mcl <= 2*max_mchirp - h_mchirp)
94
95
96
97 for lind, ns in zip(loc[good], num_slides[good]):
98
99
100 if abs(slide * ns) < 0.1:
101 continue
102
103
104
105 l = ltrigs[lind]
106 l_end_old = l.get_end()
107 l_end_tmp = l_end_old + float( slide * ns )
108
109 l.set_end(l_end_tmp)
110 epar = XLALCalculateEThincaParameter(h, l)
111
112
113
114 if epar < ethinca:
115
116 if abs(slide*ns) < 0.1:
117 print 'I USED ZEROLAG!!'
118 exit(1)
119 hcopy = copy.deepcopy(h)
120 l.set_end(h.get_end())
121 lcopy = copy.deepcopy(l)
122 hcopy.event_id = lsctables.SnglInspiralID(num_trig)
123 lcopy.event_id = lsctables.SnglInspiralID(num_trig)
124 num_trig += 1
125 coinc_sngls.append(hcopy)
126 coinc_sngls.append(lcopy)
127 l.set_end(l_end_old)
128
129
130
131 loc = loc[goods]
132
133
134 if len(loc) == 0:
135 break
136 return( coinc_sngls )
137
139 slots = param_range.split(";")
140 bins = []
141 for slot in slots:
142 minb, maxb = slot[1:-1].split(',')
143 bins.append((float(minb), float(maxb)))
144 return bins
145
147 for mchirp_low, mchirp_high in mchirp_bins:
148 if mchirp >= mchirp_low and mchirp <= mchirp_high:
149 return mchirp_low, mchirp_high
150 return None, None
151
153
154 veto_seg = segments.segmentlist()
155 veto_seg.append(segments.segment(veto_start, veto_end))
156
157
158 veto_coincs = coincs.veto(veto_seg)
159 return veto_coincs
160
162 """
163 Returns the expected total time-slid coincident livetime for two detectors
164 having single-ifo livetimes lt1, lt2, when doing all possible slides with
165 time offsets being multiples of slide_step.
166 Specifying veto_window > 0 removes a given duration from each of the
167 single-ifo livetimes.
168
169 Derivation: Consider an analysis of length T with two ifos having duty
170 factors d1, d2. Then the expected coinc livetime per slide is T*d1*d2.
171 The maximum possible number of slides is T/s where s is the step. The
172 total time is the product which can be expressed as
173 (T*d1)(T*d2)/s = lt1*lt2/s.
174 """
175 lt1 = float(lt1)
176 lt2 = float(lt2)
177 slide_step = float(slide_step)
178 veto_window = float(veto_window)
179 return (lt1 - veto_window*2) * (lt2 - veto_window*2) / slide_step
180