1 from __future__ import division
2
3 import sys
4 itertools = __import__("itertools")
5
6 import numpy
7
8 from glue import iterutils
9 from glue import segments, segmentsUtils
10 from glue.ligolw import ligolw
11 from glue.ligolw import table
12 from glue.ligolw import lsctables
13 from glue.ligolw import utils
14 from glue.ligolw.utils import ligolw_add
15 from glue.ligolw.utils import search_summary as ligolw_search_summary
16 from lal import rate
17 from pylal.SnglInspiralUtils import SnglInspiralID_old
18
19
20
21
22
24 """
25 Return a dictionary of sensitivity numbers for each detector, based on
26 a known sky location and an optional input dictionary of inspiral horizon
27 distances for a reference source of the user's choice.
28 If the horizons dictionary is specified, the returned values are interpreted
29 as inspiral horizons in that direction.
30 """
31
32 if type(gps_time)==int: gps_time=float(gps_time)
33
34 from pylal import antenna
35
36
37 if horizons is None:
38 horizons={}
39 for det in ifos:
40 horizons[det]=1.0
41 else:
42 assert len(ifos)==len(horizons)
43
44 resps={}
45
46 for det in ifos:
47 resps[det]=antenna.response(gps_time,RA,dec,0,0,'radians',det)[3]*horizons[det]
48
49 return resps
50
52 """
53 Return a set of detector thresholds adjusted for a particular
54 set of inspiral horizon distances (calculated with directional_horizon).
55 The min_threshold specified the minimum threshold which will be set
56 for all detectors less sensitive than the best one. The most sensitive
57 detector will have its threshold adjusted upward to a maximum of max_threshold.
58 """
59 assert min_threshold < max_threshold
60 threshs={}
61 worst_horizon=min(horizons.values())
62 best_horizon=max(horizons.values())
63
64 for det in horizons.keys():
65 if horizons[det]<best_horizon:
66 threshs[det]=min_threshold
67 else:
68 threshs[det]=min_threshold*(horizons[det]/worst_horizon)
69 if threshs[det]>max_threshold: threshs[det]=max_threshold
70 return threshs
71
72
73
74
75 sensitivity_dict = {"H1": 1, "L1": 2, "H2": 3, "V1": 4, "G1": 5}
77 """
78 Provide a comparison operator for IFOs such that they would get sorted
79 from most sensitive to least sensitive.
80 """
81 return cmp(sensitivity_dict[ifo1], sensitivity_dict[ifo2])
82
83
84
85
86
89 """
90 Compute and return the maximal off-source segment subject to the
91 following constraints:
92
93 1) The off-source segment is constrained to lie within a segment from the
94 analyzable segment list and to contain the on_source segment. If
95 no such segment exists, return None.
96 2) The off-source segment length is a multiple of the on-source segment
97 length. This multiple (minus one for the on-source segment) is called
98 the number of trials. By default, the number of trials is bounded
99 only by the availability of analyzable time.
100
101 Optionally:
102 3) padding_time is subtracted from the analyzable segments, but added
103 back to the off-source segment. This represents time that is thrown
104 away as part of the filtering process.
105 4) max_trials caps the number of trials that the off-source segment
106 can contain. The truncation is performed so that the resulting
107 off-source segment is as symmetric as possible.
108 5) symmetric being True will simply truncate the off-source segment to
109 be the symmetric about the on-source segment.
110 """
111 quantization_time = abs(on_source)
112
113 try:
114 super_seg = analyzable[analyzable.find(on_source)].contract(padding_time)
115 except ValueError:
116 return None
117
118
119 if on_source not in super_seg:
120 return None
121
122 nplus = (super_seg[1] - on_source[1]) // quantization_time
123 nminus = (on_source[0] - super_seg[0]) // quantization_time
124
125 if (max_trials is not None) and (nplus + nminus > max_trials):
126 half_max = max_trials // 2
127 if nplus < half_max:
128
129 remainder = max_trials - nplus
130 nminus = min(remainder, nminus)
131 elif nminus < half_max:
132
133 remainder = max_trials - nminus
134 nplus = min(remainder, nplus)
135 else:
136
137 nminus = half_max
138 nplus = max_trials - half_max
139
140 if symmetric:
141 nplus = nminus = min(nplus, nminus)
142
143 if (min_trials is not None) and (nplus + nminus < min_trials):
144 return None
145
146 return segments.segment((on_source[0] - nminus*quantization_time - padding_time,
147 on_source[1] + nplus*quantization_time + padding_time))
148
150 """
151 Return the off-source segment determined for multiple IFO times along with
152 the IFO combo that determined that segment. Calls
153 compute_offsource_segment as necessary, passing all kwargs as necessary.
154 """
155
156 new_analyzable_dict = segments.segmentlistdict()
157 for ifo, seglist in analyzable_dict.iteritems():
158 try:
159 ind = seglist.find(on_source)
160 except ValueError:
161 continue
162 new_analyzable_dict[ifo] = segments.segmentlist([seglist[ind]])
163 analyzable_ifos = new_analyzable_dict.keys()
164 analyzable_ifos.sort(sensitivity_cmp)
165
166
167
168 test_combos = itertools.chain( \
169 *itertools.imap(lambda n: iterutils.choices(analyzable_ifos, n),
170 xrange(len(analyzable_ifos), 1, -1)))
171
172 off_source_segment = None
173 the_ifo_combo = []
174 for ifo_combo in test_combos:
175 trial_seglist = new_analyzable_dict.intersection(ifo_combo)
176 temp_segment = compute_offsource_segment(trial_seglist, on_source,
177 **kwargs)
178 if temp_segment is not None:
179 off_source_segment = temp_segment
180 the_ifo_combo = list(ifo_combo)
181 the_ifo_combo.sort()
182 break
183
184 return off_source_segment, the_ifo_combo
185
187 """
188 Return the segments from a document
189 @param doc: document containing the desired segments
190 """
191
192
193 seg_dict = ligolw_search_summary.segmentlistdict_fromsearchsummary(doc)
194 segs = seg_dict.union(seg_dict.iterkeys()).coalesce()
195
196 segs = segments.segmentlist([segments.segment(int(seg[0]), int(seg[1]))\
197 for seg in segs])
198
199 return segs
200
201
203 """
204 Return a tuple of (off-source time bins, off-source veto mask,
205 index of trial that is on source).
206 The off-source veto mask is a one-dimensional boolean array where True
207 means vetoed.
208 @param onsource_doc: Document describing the on-source files
209 @param offsource_doc: Document describing the off-source files
210 @param veto_files: List of filenames containing vetoes
211 """
212
213
214 on_segs = get_segs_from_docs(onsource_doc)
215 off_segs = get_segs_from_docs(offsource_doc)
216
217 return get_exttrig_trials(on_segs, off_segs, veto_files)
218
220 """
221 Return a tuple of (off-source time bins, off-source veto mask,
222 index of trial that is on source).
223 The off-source veto mask is a one-dimensional boolean array where True
224 means vetoed.
225 @param on_segs: On-source segments
226 @param off_segs: Off-source segments
227 @param veto_files: List of filenames containing vetoes
228 """
229
230
231 trial_len = int(abs(on_segs))
232 if abs(off_segs) % trial_len != 0:
233 raise ValueError, "The provided file's analysis segment is not "\
234 "divisible by the fold time."
235 extent = (off_segs | on_segs).extent()
236
237
238 num_trials = int(abs(extent)) // trial_len
239 trial_bins = rate.LinearBins(extent[0], extent[1], num_trials)
240
241
242 trial_veto_mask = numpy.zeros(num_trials, dtype=numpy.bool8)
243 for veto_file in veto_files:
244 new_veto_segs = segmentsUtils.fromsegwizard(open(veto_file),
245 coltype=int)
246 if new_veto_segs.intersects(on_segs):
247 print >>sys.stderr, "warning: %s overlaps on-source segment" \
248 % veto_file
249 trial_veto_mask |= rate.bins_spanned(trial_bins, new_veto_segs).astype(bool)
250
251
252 onsource_mask = rate.bins_spanned(trial_bins, on_segs).astype(bool)
253 if sum(onsource_mask) != 1:
254 raise ValueError, "on-source segment spans more or less than one trial"
255 onsource_ind = numpy.arange(len(onsource_mask))[onsource_mask]
256
257 return trial_bins, trial_veto_mask, onsource_ind
258
260 """
261 Return the arithmetic average of the mchirps of all triggers in coinc.
262 """
263 return sum(t.mchirp for t in coinc) / coinc.numifos
264
265
266
267
268
269
271 doc = ligolw_add.ligolw_add(ligolw.Document(), [filename])
272 return lsctables.ExtTriggersTable.get_table(doc)
273
275 """
276 Create an empty LIGO_LW XML document, add a table of table_type,
277 insert the given rows, then write the document to a file.
278 """
279
280 xmldoc = ligolw.Document()
281 xmldoc.appendChild(ligolw.LIGO_LW())
282 tbl = lsctables.New(table_type)
283 xmldoc.childNodes[-1].appendChild(tbl)
284
285
286 tbl.extend(rows)
287
288
289 utils.write_filename(xmldoc, filename)
290
291 -def load_cache(xmldoc, cache, sieve_pattern, exact_match=False,
292 verbose=False):
293 """
294 Return a parsed and ligolw_added XML document from the files matched by
295 sieve_pattern in the given cache.
296 """
297 subcache = cache.sieve(description=sieve_pattern, exact_match=exact_match)
298 found, missed = subcache.checkfilesexist()
299 if len(found) == 0:
300 print >>sys.stderr, "warning: no files found for pattern %s" \
301 % sieve_pattern
302
303
304 old_id = lsctables.SnglInspiralTable.next_id
305 lsctables.SnglInspiralTable.next_id = SnglInspiralID_old(0)
306
307
308
309
310 urls = [c.url for c in found]
311 try:
312 xmldoc = ligolw_add.ligolw_add(ligolw.Document(), urls, verbose=verbose)
313 except ligolw.ElementError:
314
315 lsctables.SnglInspiralTable.validcolumns["event_id"] = "int_8s"
316 lsctables.SnglInspiralID = int
317 xmldoc = ligolw_add.ligolw_add(ligolw.Document(), urls, verbose=verbose)
318
319
320 lsctables.SnglInspiralTable.next_id = old_id
321
322 return xmldoc
323
325 """
326 Return the value of --num-slides found in the process_params table of
327 xmldoc. If no such entry is found, return 0.
328 """
329
330 for row in lsctables.ProcessParamsTable.get_table(xmldoc):
331 if row.param == "--num-slides":
332 return int(row.value)
333 return 0
334
335
336
337
339
340
341
342
343
344 rings = ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, program = "thinca").popitem()[1]
345
346
347
348
349
350
351 rings = segments.segmentlist(set(rings))
352 rings.sort()
353
354
355
356
357
358 for i in range(len(rings) - 1):
359 if rings[i].intersects(rings[i + 1]):
360 raise ValueError, "non-disjoint thinca rings detected in search_summary table"
361
362
363
364
365
366 for i, ring in enumerate(rings):
367 rings[i] = segments.segment(int(ring[0]), int(ring[1]))
368
369
370
371
372 return rings
373