1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 """
28 Inspiral injection identification library. Contains code providing the
29 capacity to search a list of sngl_inspiral candidates for events
30 matching entries in a sim_inspiral list of software injections, recording the
31 matches as inspiral <--> injection coincidences using the standard coincidence
32 infrastructure. Also, any pre-recorded inspiral <--> inspiral coincidences are
33 checked for cases where all the inspiral events in a coincidence match an
34 injection, and these are recorded as coinc <--> injection coincidences,
35 again using the standard coincidence infrastructure.
36 """
37
38
39 import bisect
40 import sys
41
42
43 from glue.ligolw import lsctables
44 from glue.ligolw.utils import coincs as ligolw_coincs
45 from glue.ligolw.utils import process as ligolw_process
46 from pylal import git_version
47 from lalinspiral import thinca
48 from pylal import ligolw_rinca
49 from lalburst import timeslides as ligolw_tisi
50 from pylal import SimInspiralUtils
51 from pylal import SnglInspiralUtils
52
53
54 __author__ = "Kipp Cannon <kipp.cannon@ligo.org>"
55 __version__ = "git id %s" % git_version.id
56 __date__ = git_version.date
57
58
59
60
61
62
63
64
65
66
67
69
70 return cmp(self.end_time, other.seconds) or cmp(self.end_time_ns, other.nanoseconds)
71
72
73 lsctables.SnglInspiral.__cmp__ = sngl_inspiral___cmp__
74
76
77 return cmp(self.start_time, other.seconds) or cmp(self.start_time_ns, other.nanoseconds)
78
79
80 lsctables.SnglRingdown.__cmp__ = sngl_ringdown___cmp__
81
82
83
84
85
86
87
88
89
90
91 InspiralSICoincDef = lsctables.CoincDef(search = u"inspiral", search_coinc_type = 1, description = u"sim_inspiral<-->sngl_inspiral coincidences")
92 InspiralSCNearCoincDef = lsctables.CoincDef(search = u"inspiral", search_coinc_type = 2, description = u"sim_inspiral<-->coinc_event coincidences (nearby)")
93
94 RingdownSICoincDef = lsctables.CoincDef(search = u"ringdown", search_coinc_type = 1, description = u"sim_ringdown<-->sngl_ringdown coincidences")
95 RingdownSCNearCoincDef= lsctables.CoincDef(search = u"ringdown", search_coinc_type = 2, description = u"sim_ringdown<-->coinc_event coincidences (nearby)")
96
97 -class DocContents(object):
98 """
99 A wrapper interface to the XML document.
100 """
101 - def __init__(self, xmldoc, bbdef, sbdef, scndef, process, sngl_type, sim_type, get_sngl_time):
102
103
104
105
106 self.process = process
107
108
109
110
111
112 self.sngltable = sngl_type.get_table(xmldoc)
113 try:
114 self.simtable = sim_type.get_table(xmldoc)
115 except ValueError:
116 self.simtable = lsctables.SimInspiralTable.get_table(xmldoc)
117 print >>sys.stderr,"No SimRingdownTable, use SimInspiralTable instead!"
118
119
120
121
122
123
124
125
126
127
128 self.tisi_id = ligolw_tisi.get_time_slide_id(xmldoc, {}.fromkeys(self.sngltable.getColumnByName("ifo"), 0.0), create_new = process)
129
130
131
132
133
134
135
136 self.sb_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, sbdef.search, sbdef.search_coinc_type, create_new = True, description = sbdef.description)
137
138
139
140
141
142
143
144
145 try:
146 bb_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, bbdef.search, bbdef.search_coinc_type, create_new = False)
147 except KeyError:
148 bb_coinc_def_id = None
149 self.scn_coinc_def_id = None
150 else:
151 self.scn_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, scndef.search, scndef.search_coinc_type, create_new = True, description = scndef.description)
152
153
154
155
156
157 try:
158 self.coinctable = lsctables.CoincTable.get_table(xmldoc)
159 except ValueError:
160 self.coinctable = lsctables.New(lsctables.CoincTable)
161 xmldoc.childNodes[0].appendChild(self.coinctable)
162 self.coinctable.sync_next_id()
163
164
165
166
167
168 try:
169 self.coincmaptable = lsctables.CoincMapTable.get_table(xmldoc)
170 except ValueError:
171 self.coincmaptable = lsctables.New(lsctables.CoincMapTable)
172 xmldoc.childNodes[0].appendChild(self.coincmaptable)
173
174
175
176
177
178
179
180
181
182
183 index = {}
184 for row in self.sngltable:
185 index[row.event_id] = row
186
187 self.coincs = {}
188 for coinc in self.coinctable:
189 if coinc.coinc_def_id == bb_coinc_def_id:
190 self.coincs[coinc.coinc_event_id] = []
191
192 for row in self.coincmaptable:
193 if row.coinc_event_id in self.coincs:
194 self.coincs[row.coinc_event_id].append(index[row.event_id])
195 del index
196
197
198
199 for coinc_event_id, events in self.coincs.iteritems():
200 events.sort(key=get_sngl_time)
201 self.coincs[coinc_event_id]=tuple(events)
202
203
204 self.coincs = self.coincs.items()
205
206
207
208
209
210
211
212
213
214 self.sngltable.sort(key=get_sngl_time)
215 self.coincs.sort(key=lambda(id,a): get_sngl_time(a[0]))
216
217
218
219
220
221
222
223
224
225
226 self.search_time_window = 1.0
227 self.coinc_time_window = 1.0
228
229
230 - def type_near_time(self, t):
231 """
232 Return a list of the inspiral events whose peak times are
233 within self.search_time_window of t.
234 """
235 return self.sngltable[bisect.bisect_left(self.sngltable, t - self.search_time_window):bisect.bisect_right(self.sngltable, t + self.search_time_window)]
236
237 - def coincs_near_time(self, t, get_time):
238 """
239 Return a list of the (coinc_event_id, event list) tuples in
240 which at least one type event's end time is within
241 self.coinc_time_window of t.
242 """
243
244
245
246
247 return [(coinc_event_id, search_type) for coinc_event_id, search_type in self.coincs if (t - self.coinc_time_window <= get_time(search_type[-1])) and (get_time(search_type[0]) <= t + self.coinc_time_window)]
248
250 """
251 Sort the sngl_table's rows by ID (tidy-up document
252 for output).
253 """
254 self.sngltable.sort(lambda a, b: cmp(a.event_id, b.event_id))
255
256 - def new_coinc(self, coinc_def_id):
257 """
258 Construct a new coinc_event row attached to the given
259 process, and belonging to the set of coincidences defined
260 by the given coinc_def_id.
261 """
262 coinc = lsctables.Coinc()
263 coinc.process_id = self.process.process_id
264 coinc.coinc_def_id = coinc_def_id
265 coinc.coinc_event_id = self.coinctable.get_next_id()
266 coinc.time_slide_id = self.tisi_id
267 coinc.set_instruments(None)
268 coinc.nevents = 0
269 coinc.likelihood = None
270 self.coinctable.append(coinc)
271 return coinc
272
273
274
275
276
277
278
279
280
281
282
283 process_program_name = "lalapps_cbc_injfind"
284
285
287 """
288 Convenience wrapper for adding process metadata to the document.
289 """
290 process = ligolw_process.append_process(xmldoc, program = process_program_name, version = __version__, cvs_repository = u"lscsoft", cvs_entry_time = __date__, comment = comment)
291
292 params = [(u"--match-algorithm", u"lstring", match_algorithm)]
293 ligolw_process.append_process_params(xmldoc, process, params)
294
295 return process
296
297
298
299
300
301
302
303
304
305
311
312
318
320 tdiff = abs(get_sngl_time(sngl) - get_sim_time(sim))
321 if tdiff < twindow:
322 return 0
323 else:
324 return cmp(get_sngl_time(sngl), get_sim_time(sim))
325
326
327
328
329
330
331
332
333
334
335
337 """
338 Scan the type table for triggers matching sim.
339 """
340 return [type_table for type_table in contents.type_near_time(get_sim_time(sim)) if not comparefunc(sim, type_table, get_sim_time, get_sngl_time)]
341
342
344 """
345 Create a coinc_event in the coinc table, and add arcs in the
346 coinc_event_map table linking the sim_type row and the list of
347 sngl_type rows to the new coinc_event row.
348 """
349 coinc = contents.new_coinc(contents.sb_coinc_def_id)
350 coinc.set_instruments(set(event.ifo for event in types))
351 coinc.nevents = len(types)
352
353 coincmap = lsctables.CoincMap()
354 coincmap.coinc_event_id = coinc.coinc_event_id
355 coincmap.table_name = sim.simulation_id.table_name
356 coincmap.event_id = sim.simulation_id
357 contents.coincmaptable.append(coincmap)
358
359 for event in types:
360 coincmap = lsctables.CoincMap()
361 coincmap.coinc_event_id = coinc.coinc_event_id
362 coincmap.table_name = event.event_id.table_name
363 coincmap.event_id = event.event_id
364 contents.coincmaptable.append(coincmap)
365
366
367
368
369
370
371
372
373
374
375
377 """
378 Return a list of the coinc_event_ids of the inspiral<-->inspiral coincs
379 in which all inspiral events match sim.
380 """
381
382
383
384 return [coinc_event_id for coinc_event_id, inspirals in coincs if True not in [bool(comparefunc(sim, inspiral)) for inspiral in inspirals]]
385
387 """
388 Return a list of the coinc_event_ids of the type<-->type coincs
389 in which at least one type event matches sim.
390 """
391
392
393
394 return [coinc_event_id for coinc_event_id, coinc_types in coincs if False in [bool(comparefunc(sim, coinc_type, get_sim_time, get_sngl_time)) for coinc_type in coinc_types]]
395
396
398 """
399 Create a coinc_event in the coinc table, and add arcs in the
400 coinc_event_map table linking the sim_type row and the list of
401 coinc_event rows to the new coinc_event row.
402 """
403 coinc = contents.new_coinc(coinc_def_id)
404 coinc.nevents = len(coinc_event_ids)
405
406 coincmap = lsctables.CoincMap()
407 coincmap.coinc_event_id = coinc.coinc_event_id
408 coincmap.table_name = sim.simulation_id.table_name
409 coincmap.event_id = sim.simulation_id
410 contents.coincmaptable.append(coincmap)
411
412 for coinc_event_id in coinc_event_ids:
413 coincmap = lsctables.CoincMap()
414 coincmap.coinc_event_id = coinc.coinc_event_id
415 coincmap.table_name = coinc_event_id.table_name
416 coincmap.event_id = coinc_event_id
417 contents.coincmaptable.append(coincmap)
418
419
420
421
422
423
424
425
426
427
428
429 -def lalapps_cbc_injfind(xmldoc, process, search, snglcomparefunc, nearcoinccomparefunc, verbose = False):
430
431
432
433
434 if verbose:
435 print >>sys.stderr, "indexing ..."
436
437 bbdef = {"inspiral": thinca.InspiralCoincDef, "ringdown": ligolw_rinca.RingdownCoincDef}[search]
438 sbdef = {"inspiral": InspiralSICoincDef, "ringdown": RingdownSICoincDef}[search]
439 scndef = {"inspiral": InspiralSCNearCoincDef, "ringdown": RingdownSCNearCoincDef}[search]
440 sngl_type = {"inspiral": lsctables.SnglInspiralTable, "ringdown": lsctables.SnglRingdownTable}[search]
441 sim_type = {"inspiral": lsctables.SimInspiralTable, "ringdown": lsctables.SimRingdownTable}[search]
442 get_sngl_time = {"inspiral": lsctables.SnglInspiral.get_end, "ringdown": lsctables.SnglRingdown.get_start}[search]
443
444 contents = DocContents(xmldoc = xmldoc, bbdef = bbdef, sbdef = sbdef, scndef = scndef, process = process, sngl_type = sngl_type, sim_type = sim_type, get_sngl_time = get_sngl_time)
445 N = len(contents.simtable)
446
447 if isinstance(contents.simtable, lsctables.SimInspiralTable):
448 get_sim_time = lsctables.SimInspiral.get_end
449 elif isinstance(contents.simtable, lsctables.SimRingdownTable):
450 get_sim_time = lsctables.SimRingdown.get_start
451 else:
452 raise ValueError, "Unknown sim table"
453
454
455
456
457
458 if verbose:
459 print >>sys.stderr, "constructing %s:" % sbdef.description
460 for n, sim in enumerate(contents.simtable):
461 if verbose:
462 print >>sys.stderr, "\t%.1f%%\r" % (100.0 * n / N),
463 types = find_sngl_type_matches(contents, sim, snglcomparefunc, get_sim_time, get_sngl_time)
464 if types:
465 add_sim_type_coinc(contents, sim, types)
466 if verbose:
467 print >>sys.stderr, "\t100.0%"
468
469
470
471
472
473 if contents.scn_coinc_def_id:
474 if verbose:
475 print >>sys.stderr, "constructing %s:" % (scndef.description)
476 for n, sim in enumerate(contents.simtable):
477 if verbose:
478 print >>sys.stderr, "\t%.1f%%\r" % (100.0 * n / N),
479 coincs = contents.coincs_near_time(get_sim_time(sim), get_sngl_time)
480 coinc_event_ids = find_near_coinc_matches(coincs, sim, nearcoinccomparefunc, get_sim_time, get_sngl_time)
481 if coinc_event_ids:
482 add_sim_coinc_coinc(contents, sim, coinc_event_ids, contents.scn_coinc_def_id)
483 if verbose:
484 print >>sys.stderr, "\t100.0%"
485
486
487
488
489
490
491
492 if verbose:
493 print >>sys.stderr, "Checking for both SimInspiralTable and SimRingdownTable"
494 if isinstance(contents.simtable, lsctables.SimRingdownTable):
495 try:
496 insp_simtable = lsctables.SimInspiralTable.get_table(xmldoc)
497 except ValueError:
498 print >>sys.stderr,"No SimInspiralTable, only SimRingdownTable present"
499 insp_simtable = None
500 if insp_simtable is not None:
501 if verbose:
502 print >> sys.stderr, "found SimInspiralTable, creating maps to SimRingdownTable"
503
504 sim_ring_time_map = dict([ [(str(row.process_id), row.geocent_start_time), str(row.simulation_id)] for row in contents.simtable])
505
506 sim_ring_ceid_map = dict([ [str(row.event_id), row.coinc_event_id] for row in contents.coincmaptable if row.table_name == "sim_ringdown" ])
507
508 for this_sim_insp in insp_simtable:
509 if (str(this_sim_insp.process_id), this_sim_insp.geocent_end_time) in sim_ring_time_map:
510 this_sim_ring_id = sim_ring_time_map[( str(this_sim_insp.process_id), this_sim_insp.geocent_end_time )]
511 else:
512 continue
513 if str(this_sim_ring_id) in sim_ring_ceid_map:
514 this_ceid = sim_ring_ceid_map[ str(this_sim_ring_id) ]
515 else:
516 continue
517
518 new_mapping = lsctables.CoincMap()
519 new_mapping.table_name = "sim_inspiral"
520 new_mapping.event_id = this_sim_insp.simulation_id
521 new_mapping.coinc_event_id = this_ceid
522 contents.coincmaptable.append(new_mapping)
523
524
525
526
527
528 if verbose:
529 print >>sys.stderr, "finishing ..."
530 contents.sort_triggers_by_id()
531
532
533
534
535
536 return xmldoc
537