Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALBurst 2.0.7.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Macros Modules Pages
binjfind.py
Go to the documentation of this file.
1# Copyright (C) 2006-2021 Kipp Cannon
2#
3# This program is free software; you can redistribute it and/or modify it
4# under the terms of the GNU General Public License as published by the
5# Free Software Foundation; either version 2 of the License, or (at your
6# option) any later version.
7#
8# This program is distributed in the hope that it will be useful, but
9# WITHOUT ANY WARRANTY; without even the implied warranty of
10# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11# Public License for more details.
12#
13# You should have received a copy of the GNU General Public License along
14# with this program; if not, write to the Free Software Foundation, Inc.,
15# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
16
17
18#
19# =============================================================================
20#
21# Preamble
22#
23# =============================================================================
24#
25
26
27"""
28Burst injection identification library. Contains code providing the
29capacity to search a list of sngl_burst burst candidates for events
30matching entries in a sim_burst list of software injections, recording the
31matches as burst <--> injection coincidences using the standard coincidence
32infrastructure. Also, any pre-recorded burst <--> burst coincidences are
33checked for cases where all the burst events in a coincidence match an
34injection, and these are recorded as coinc <--> injection coincidences,
35again using the standard coincidence infrastructure.
36"""
37
38
39import bisect
40import sys
41
42
43from igwn_ligolw import lsctables
44from igwn_ligolw.utils import coincs as ligolw_coincs
45from igwn_ligolw.utils import process as ligolw_process
46from igwn_ligolw.utils import search_summary as ligolw_search_summary
47from igwn_ligolw.utils import time_slide as ligolw_time_slide
48import lal
49import igwn_segments as segments
50from . import burca
51from . import SimBurstUtils
52
53
54__author__ = "Kipp Cannon <kipp.cannon@ligo.org>"
55from .git_version import date as __date__
56from .git_version import version as __version__
57
58
59#
60# =============================================================================
61#
62# Speed Hacks
63#
64# =============================================================================
65#
66
67
68#
69# Construct a subclass of the sngl_burst row class with the methods that
70# are needed
71#
72
73
74class SnglBurst(lsctables.SnglBurst):
75 __slots__ = ()
76
77 #
78 # compare self's end time to the LIGOTimeGPS instance other.
79 # allows bisection searches by GPS time to find ranges of triggers
80 # quickly
81 #
82
83 def __lt__(self, other):
84 return self.end < other
85
86 def __le__(self, other):
87 return self.end <= other
88
89 def __eq__(self, other):
90 return self.end == other
91
92 def __ne__(self, other):
93 return self.end != other
94
95 def __ge__(self, other):
96 return self.end >= other
97
98 def __gt__(self, other):
99 return self.end > other
100
101
102lsctables.SnglBurstTable.RowType = lsctables.SnglBurst = SnglBurst
103
104
105#
106# place-holder class for sim_inspiral table's row class so that
107# time_slide_id attributes can be assigned to it (this would normally be
108# prevented by the __slots__ mechanism). FIXME: delete when the
109# sim_inspiral table has a time_slide_id column
110#
111
112
113class SimInspiral(lsctables.SimInspiral):
114 pass
115
116
117lsctables.SimInspiralTable.RowType = lsctables.SimInspiral = SimInspiral
118
119
120#
121# =============================================================================
122#
123# Document Interface
124#
125# =============================================================================
126#
127
128
129ExcessPowerSBBCoincDef = lsctables.CoincDef(search = "excesspower", search_coinc_type = 1, description = "sim_burst<-->sngl_burst coincidences")
130ExcessPowerSIBCoincDef = lsctables.CoincDef(search = "excesspower", search_coinc_type = 4, description = "sim_inspiral<-->sngl_burst coincidences")
131ExcessPowerSBCCoincDef = lsctables.CoincDef(search = "excesspower", search_coinc_type = 2, description = "sim_burst<-->coinc_event coincidences (exact)")
132ExcessPowerSBCNearCoincDef = lsctables.CoincDef(search = "excesspower", search_coinc_type = 3, description = "sim_burst<-->coinc_event coincidences (nearby)")
133ExcessPowerSICCoincDef = lsctables.CoincDef(search = "excesspower", search_coinc_type = 5, description = "sim_inspiral<-->coinc_event coincidences (exact)")
134ExcessPowerSICNearCoincDef = lsctables.CoincDef(search = "excesspower", search_coinc_type = 6, description = "sim_inspiral<-->coinc_event coincidences (nearby)")
135
136
137StringCuspSBBCoincDef = lsctables.CoincDef(search = "StringCusp", search_coinc_type = 1, description = "sim_burst<-->sngl_burst coincidences")
138StringCuspSIBCoincDef = lsctables.CoincDef(search = "StringCusp", search_coinc_type = 4, description = "sim_inspiral<-->sngl_burst coincidences")
139StringCuspSBCCoincDef = lsctables.CoincDef(search = "StringCusp", search_coinc_type = 2, description = "sim_burst<-->coinc_event coincidences (exact)")
140StringCuspSBCNearCoincDef = lsctables.CoincDef(search = "StringCusp", search_coinc_type = 3, description = "sim_burst<-->coinc_event coincidences (nearby)")
141StringCuspSICCoincDef = lsctables.CoincDef(search = "StringCusp", search_coinc_type = 5, description = "sim_inspiral<-->coinc_event coincidences (exact)")
142StringCuspSICNearCoincDef = lsctables.CoincDef(search = "StringCusp", search_coinc_type = 6, description = "sim_inspiral<-->coinc_event coincidences (nearby)")
143
144
145OmegaBBCoincDef = lsctables.CoincDef(search = "omega", search_coinc_type = 0, description = "sngl_burst<-->sngl_burst coincidences")
146OmegaSBBCoincDef = lsctables.CoincDef(search = "omega", search_coinc_type = 1, description = "sim_burst<-->sngl_burst coincidences")
147OmegaSIBCoincDef = lsctables.CoincDef(search = "omega", search_coinc_type = 4, description = "sim_inspiral<-->sngl_burst coincidences")
148OmegaSBCCoincDef = lsctables.CoincDef(search = "omega", search_coinc_type = 2, description = "sim_burst<-->coinc_event coincidences (exact)")
149OmegaSBCNearCoincDef = lsctables.CoincDef(search = "omega", search_coinc_type = 3, description = "sim_burst<-->coinc_event coincidences (nearby)")
150OmegaSICCoincDef = lsctables.CoincDef(search = "omega", search_coinc_type = 5, description = "sim_inspiral<-->coinc_event coincidences (exact)")
151OmegaSICNearCoincDef = lsctables.CoincDef(search = "omega", search_coinc_type = 6, description = "sim_inspiral<-->coinc_event coincidences (nearby)")
152
153CWBBBCoincDef = lsctables.CoincDef(search = "waveburst", search_coinc_type = 0, description = "sngl_burst<-->sngl_burst coincidences")
154CWBSBBCoincDef = lsctables.CoincDef(search = "waveburst", search_coinc_type = 1, description = "sim_burst<-->sngl_burst coincidences")
155CWBSIBCoincDef = lsctables.CoincDef(search = "waveburst", search_coinc_type = 4, description = "sim_inspiral<-->sngl_burst coincidences")
156CWBSBCCoincDef = lsctables.CoincDef(search = "waveburst", search_coinc_type = 2, description = "sim_burst<-->coinc_event coincidences (exact)")
157CWBSBCNearCoincDef = lsctables.CoincDef(search = "waveburst", search_coinc_type = 3, description = "sim_burst<-->coinc_event coincidences (nearby)")
158CWBSICCoincDef = lsctables.CoincDef(search = "waveburst", search_coinc_type = 5, description = "sim_inspiral<-->coinc_event coincidences (exact)")
159CWBSICNearCoincDef = lsctables.CoincDef(search = "waveburst", search_coinc_type = 6, description = "sim_inspiral<-->coinc_event coincidences (nearby)")
160
161
162class DocContents(object):
163 """
164 A wrapper interface to the XML document.
165 """
166 def __init__(self, xmldoc, b_b_def, sb_b_def, si_b_def, sb_c_e_def, sb_c_n_def, si_c_e_def, si_c_n_def, process, livetime_program):
167 #
168 # store the process row
169 #
170
171 self.process = process
172
173 #
174 # locate the sngl_burst, time_slide and injection tables
175 #
176
177 self.snglbursttable = lsctables.SnglBurstTable.get_table(xmldoc)
178 try:
179 self.simbursttable = lsctables.SimBurstTable.get_table(xmldoc)
180 except ValueError:
181 self.simbursttable = None
182 try:
183 self.siminspiraltable = lsctables.SimInspiralTable.get_table(xmldoc)
184 except ValueError:
185 self.siminspiraltable = None
186 try:
187 timeslidetable = lsctables.TimeSlideTable.get_table(xmldoc)
188 except ValueError:
189 timeslidetable = None
190 if timeslidetable is not None:
191 self.offsetvectors = timeslidetable.as_dict()
192 else:
193 self.offsetvectors = {}
194
195 #
196 # store the longest duration of a burst event
197 #
198
199 if self.snglbursttable:
200 self.longestduration = max(self.snglbursttable.getColumnByName("duration"))
201 else:
202 self.longestduration = 0.0
203
204 #
205 # get a segmentlistdict indicating the times when the jobs
206 # that produced burst events could have produced output
207 #
208
209 self.seglists = ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, livetime_program).coalesce()
210
211 #
212 # FIXME: in the future, the sim_inspiral table should
213 # indicate time slide at which the injection was done, like
214 # the sim_burst table does. for now, fake it by giving
215 # each one a time_slide_id attribute. this is the only
216 # reason why a place-holder class is required for the
217 # SimInspiral row type; that class can be deleted when
218 # this is no longer needed
219 #
220
221 if self.siminspiraltable is not None:
222 time_slide_id = ligolw_time_slide.get_time_slide_id(xmldoc, {}.fromkeys(self.seglists, 0.0), create_new = process, superset_ok = True, nonunique_ok = False)
223 for sim in self.siminspiraltable:
224 sim.time_slide_id = time_slide_id
225
226 #
227 # get coinc_definer rows for sim_* <--> sngl_burst coincs
228 # for whichever sim_* tables are present; this creates a
229 # coinc_definer table if the document doesn't have one
230 #
231
232 if self.simbursttable is not None:
233 self.sb_b_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, sb_b_def.search, sb_b_def.search_coinc_type, create_new = True, description = sb_b_def.description)
234 else:
235 self.sb_b_coinc_def_id = None
236 if self.siminspiraltable is not None:
237 self.si_b_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, si_b_def.search, si_b_def.search_coinc_type, create_new = True, description = si_b_def.description)
238 else:
239 self.si_b_coinc_def_id = None
240
241 #
242 # get coinc_def_id's for sngl_burst <--> sngl_burst, and
243 # both kinds of sim_* <--> coinc_event coincs. set all to
244 # None if this document does not contain any sngl_burst
245 # <--> sngl_burst coincs.
246 #
247
248 try:
249 b_b_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, b_b_def.search, b_b_def.search_coinc_type, create_new = False)
250 except KeyError:
251 b_b_coinc_def_id = None
256 else:
257 if self.simbursttable is not None:
258 self.sb_c_e_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, sb_c_e_def.search, sb_c_e_def.search_coinc_type, create_new = True, description = sb_c_e_def.description)
259 self.sb_c_n_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, sb_c_n_def.search, sb_c_n_def.search_coinc_type, create_new = True, description = sb_c_n_def.description)
260 else:
261 self.sb_c_e_coinc_def_id = None
262 self.sb_c_n_coinc_def_id = None
263 if self.siminspiraltable is not None:
264 self.si_c_e_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, si_c_e_def.search, si_c_e_def.search_coinc_type, create_new = True, description = si_c_e_def.description)
265 self.si_c_n_coinc_def_id = ligolw_coincs.get_coinc_def_id(xmldoc, si_c_n_def.search, si_c_n_def.search_coinc_type, create_new = True, description = si_c_n_def.description)
266 else:
267 self.si_c_e_coinc_def_id = None
268 self.si_c_n_coinc_def_id = None
269
270 #
271 # get coinc table, create one if needed
272 #
273
274 try:
275 self.coinctable = lsctables.CoincTable.get_table(xmldoc)
276 except ValueError:
277 self.coinctable = lsctables.New(lsctables.CoincTable)
278 xmldoc.childNodes[0].appendChild(self.coinctable)
279 self.coinctable.sync_next_id()
280
281 #
282 # get coinc_map table, create one if needed
283 #
284
285 try:
286 self.coincmaptable = lsctables.CoincMapTable.get_table(xmldoc)
287 except ValueError:
288 self.coincmaptable = lsctables.New(lsctables.CoincMapTable)
289 xmldoc.childNodes[0].appendChild(self.coincmaptable)
290
291 #
292 # index the document
293 #
294
295 # index sngl_burst table
296 index = dict((row.event_id, row) for row in self.snglbursttable)
297 # find IDs of burst<-->burst coincs
298 self.coincs = dict((row.coinc_event_id, []) for row in self.coinctable if row.coinc_def_id == b_b_coinc_def_id)
299 # construct event list for each burst<-->burst coinc
300 for row in self.coincmaptable:
301 try:
302 self.coincs[row.coinc_event_id].append(index[row.event_id])
303 except KeyError:
304 pass
305 del index
306 # sort the event list for each coin by peak time and
307 # convert to tuples for speed
308 for coinc_event_id, events in self.coincs.items():
309 events.sort(key = lambda event: event.peak)
310 self.coincs[coinc_event_id] = tuple(events)
311 # convert dictionary to a list of (coinc_event_id, events)
312 # tuples and create a coinc_event_id to offset vector
313 # look-up table
314 self.coincs = self.coincs.items()
315 self.coincoffsets = dict((row.coinc_event_id, self.offsetvectors[row.time_slide_id]) for row in self.coinctable if row.coinc_def_id == b_b_coinc_def_id)
316
317 #
318 # represent the sngl_burst table as a dictionary indexed by
319 # instrument whose values are the event lists for those
320 # instruments sorted by peak time. sort the coincs list by
321 # the peak time of the first (earliest) event in each coinc
322 # (recall that the event tuple for each coinc has been
323 # time-ordered)
324 #
325
326 self.snglbursttable = dict((instrument, sorted((event for event in self.snglbursttable if event.ifo == instrument), key = lambda event: event.peak)) for instrument in set(self.snglbursttable.getColumnByName("ifo")))
327 self.coincs.sort(key = lambda coinc_id_events: coinc_id_events[1][0].peak)
328
329 def bursts_near_peaktime(self, t, window, offsetvector):
330 """
331 Return a list of the burst events whose peak times (with
332 offsetvector added) are within window seconds of t. This
333 is not used to define any coincidences, only to provide a
334 short list of burst events for use in more costly
335 comparison tests.
336 """
337 # instead of adding the offsets to each burst trigger, the
338 # offsets are subtracted from the times against which the
339 # bursts are compared
340 return sum((events[bisect.bisect_left(events, t - offsetvector[instrument] - window):bisect.bisect_right(events, t - offsetvector[instrument] + window)] for instrument, events in self.snglbursttable.items()), [])
341
342 def coincs_near_peaktime(self, t, window, offsetvector):
343 """
344 Return a list of the (coinc_event_id, event list) tuples in
345 which at least one burst event's peak time is within window
346 seconds of t. This is not used to define any coincidences,
347 only to provide a short list of coinc events for use in
348 more costly comparison tests.
349 """
350 #
351 # burst events within window of t
352 #
353
354 near_events = set(self.bursts_near_peaktime(t, window, offsetvector))
355
356 #
357 # coincs that involve at least one of those burst events
358 # and were found following the application of an offset
359 # vector consistent with the injection offset vector
360 #
361
362 return [(coinc_event_id, events) for coinc_event_id, events in self.coincs if set(events) & near_events and offsetvector.contains(self.coincoffsets[coinc_event_id])]
363
364 def new_coinc(self, coinc_def_id, time_slide_id):
365 """
366 Construct a new coinc_event row attached to the given
367 process, and belonging to the set of coincidences defined
368 by the given coinc_def_id.
369 """
370 coinc = self.coinctable.RowType()
371 coinc.process_id = self.process.process_id
372 coinc.coinc_def_id = coinc_def_id
373 coinc.coinc_event_id = self.coinctable.get_next_id()
374 coinc.time_slide_id = time_slide_id
375 coinc.set_instruments(None)
376 coinc.nevents = 0
377 coinc.likelihood = None
378 self.coinctable.append(coinc)
379 return coinc
380
381
382#
383# =============================================================================
384#
385# Add Process Information
386#
387# =============================================================================
388#
389
390
391process_program_name = "lalburst_injfind"
392
393
394def append_process(xmldoc, match_algorithm, comment):
395 """
396 Convenience wrapper for adding process metadata to the document.
397 """
398 return ligolw_process.register_to_xmldoc(
399 xmldoc,
400 program = process_program_name,
401 paramdict = {
402 "match_algorithm": match_algorithm
403 },
404 version = __version__,
405 cvs_repository = "lscsoft",
406 cvs_entry_time = __date__,
407 comment = comment
408 )
409
410
411#
412# =============================================================================
413#
414# Injection <--> Burst Event Comparison Tests
415#
416# =============================================================================
417#
418
419
420def StringCuspSnglCompare(sim, burst, offsetvector):
421 """
422 Return False (injection matches event) if an autocorrelation-width
423 window centred on the injection is continuous with the time
424 interval of the burst.
425 """
426 tinj = sim.time_at_instrument(burst.ifo, offsetvector)
427 window = SimBurstUtils.stringcusp_autocorrelation_width / 2
428 # uncomment last part of expression to impose an amplitude cut
429 return segments.segment(tinj - window, tinj + window).disjoint(burst.period) #or abs(sim.amplitude / SimBurstUtils.string_amplitude_in_instrument(sim, burst.ifo, offsetvector)) > 3
430
431
432def ExcessPowerSnglCompare(sim, burst, offsetvector):
433 """
434 Return False (injection matches event) if the peak time and centre
435 frequency of sim lie within the time-frequency tile of burst.
436 """
437 return (sim.time_at_instrument(burst.ifo, offsetvector) not in burst.period) or (sim.frequency not in burst.band)
438
439
440def OmegaSnglCompare(sim, burst, offsetvector, delta_t = 10.0):
441 """
442 Return False (injection matches event) if the time of the sim and
443 the peak time of the burst event differ by less than or equal to
444 delta_t seconds.
445 """
446 return abs(float(sim.time_at_instrument(burst.ifo, offsetvector) - burst.peak)) > delta_t
447
448def CWBSnglCompare(sim, burst, offsetvector, delta_t = 10.0):
449 """
450 Return False (injection matches event) if the time of the sim and
451 the peak time of the burst event differ by less than or equal to
452 delta_t seconds.
453 """
454 return abs(float(sim.time_at_instrument(burst.ifo, offsetvector) - burst.peak)) > delta_t
455
456
457def StringCuspNearCoincCompare(sim, burst, offsetvector):
458 """
459 Return False (injection matches coinc) if the peak time of the sim
460 is "near" the burst event.
461 """
462 tinj = sim.time_at_instrument(burst.ifo, offsetvector)
463 window = SimBurstUtils.stringcusp_autocorrelation_width / 2 + SimBurstUtils.burst_is_near_injection_window
464 return segments.segment(tinj - window, tinj + window).disjoint(burst.period)
465
466
467def ExcessPowerNearCoincCompare(sim, burst, offsetvector):
468 """
469 Return False (injection matches coinc) if the peak time of the sim
470 is "near" the burst event.
471 """
472 tinj = sim.time_at_instrument(burst.ifo, offsetvector)
473 window = SimBurstUtils.burst_is_near_injection_window
474 return segments.segment(tinj - window, tinj + window).disjoint(burst.period)
475
476
477def OmegaNearCoincCompare(sim, burst, offsetvector):
478 """
479 Return False (injection matches coinc) if the peak time of the sim
480 is "near" the burst event.
481 """
482 return OmegaSnglCompare(sim, burst, offsetvector, delta_t = 20.0 + burst.duration / 2)
483
484def CWBNearCoincCompare(sim, burst, offsetvector):
485 """
486 Return False (injection matches coinc) if the peak time of the sim
487 is "near" the burst event.
488 """
489 return OmegaSnglCompare(sim, burst, offsetvector, delta_t = 20.0 + burst.duration / 2)
490
491
492#
493# =============================================================================
494#
495# Build sim_* <--> sngl_burst Coincidences
496#
497# =============================================================================
498#
499
500
501def find_sngl_burst_matches(contents, sim, comparefunc, sieve_window):
502 """
503 Scan the burst table for triggers matching sim. sieve_window is
504 used in a bisection search to quickly identify burst events within
505 that many seconds of the injection's peak time at the geocentre;
506 it should be larger than the greatest time difference that can
507 separate a burst event's peak time from an injection's peak time at
508 the geocentre and the two still be considered a match.
509 """
510 offsetvector = contents.offsetvectors[sim.time_slide_id]
511 return [burst for burst in contents.bursts_near_peaktime(sim.time_geocent, sieve_window, offsetvector) if not comparefunc(sim, burst, offsetvector)]
512
513
514def add_sim_burst_coinc(contents, sim, events, coinc_def_id):
515 """
516 Create a coinc_event in the coinc table, and add arcs in the
517 coinc_event_map table linking the sim_burst row and the list of
518 sngl_burst rows to the new coinc_event row.
519 """
520 # the coinc always carries the same time_slide_id as the injection
521 coinc = contents.new_coinc(coinc_def_id, sim.time_slide_id)
522 coinc.set_instruments(set(event.ifo for event in events))
523 coinc.nevents = len(events)
524
525 coincmap = contents.coincmaptable.RowType()
526 coincmap.coinc_event_id = coinc.coinc_event_id
527 coincmap.table_name = sim.simulation_id.table_name
528 coincmap.event_id = sim.simulation_id
529 contents.coincmaptable.append(coincmap)
530
531 for event in events:
532 coincmap = contents.coincmaptable.RowType()
533 coincmap.coinc_event_id = coinc.coinc_event_id
534 coincmap.table_name = event.event_id.table_name
535 coincmap.event_id = event.event_id
536 contents.coincmaptable.append(coincmap)
537
538 return coinc
539
540
541#
542# =============================================================================
543#
544# Build sim_burst <--> coinc Coincidences
545#
546# =============================================================================
547#
548
549
550def find_exact_coinc_matches(coincs, sim, comparefunc, seglists, offsetvector):
551 """
552 Return a set of the coinc_event_ids of the burst<-->burst coincs in
553 which all burst events match sim and to which all instruments on at
554 the time of the sim contributed events.
555 """
556 # note: this doesn't check that the coinc and the sim share
557 # compatible offset vectors, it is assumed this condition was
558 # applied in the .coincs_near_peaktime() method that was used to
559 # assemble the list of candidate coincs
560
561 # comparefunc() returns False --> event matches sim
562 # any() --> at least one compare returns True == not all events match
563 # not any() --> all events match sim
564 on_instruments = SimBurstUtils.on_instruments(sim, seglists, offsetvector)
565 return set(coinc_event_id for coinc_event_id, events in coincs if on_instruments.issubset(set(event.ifo for event in events)) and not any(comparefunc(sim, event, offsetvector) for event in events))
566
567
568def find_near_coinc_matches(coincs, sim, comparefunc, offsetvector):
569 """
570 Return a set of the coinc_event_ids of the burst<-->burst coincs in
571 which at least one burst event matches sim.
572 """
573 # note: this doesn't check that the coinc and the sim share
574 # compatible offset vectors, it is assumed this condition was
575 # applied in the .coincs_near_peaktime() method that was used to
576 # assemble the list of candidate coincs
577
578 # comparefunc() returns False --> event matches sim
579 # all() --> all compares return True == no events match
580 # not all() --> at least one event matches sim
581 return set(coinc_event_id for coinc_event_id, events in coincs if not all(comparefunc(sim, event, offsetvector) for event in events))
582
583
584def add_sim_coinc_coinc(contents, sim, coinc_event_ids, coinc_def_id):
585 """
586 Create a coinc_event in the coinc table, and add arcs in the
587 coinc_event_map table linking the sim_burst row and the list of
588 coinc_event rows to the new coinc_event row.
589 """
590 # the coinc always carries the same time_slide_id as the injection
591 coinc = contents.new_coinc(coinc_def_id, sim.time_slide_id)
592 coinc.nevents = len(coinc_event_ids)
593
594 coincmap = contents.coincmaptable.RowType()
595 coincmap.coinc_event_id = coinc.coinc_event_id
596 coincmap.table_name = sim.simulation_id.table_name
597 coincmap.event_id = sim.simulation_id
598 contents.coincmaptable.append(coincmap)
599
600 for coinc_event_id in coinc_event_ids:
601 coincmap = contents.coincmaptable.RowType()
602 coincmap.coinc_event_id = coinc.coinc_event_id
603 coincmap.table_name = coinc_event_id.table_name
604 coincmap.event_id = coinc_event_id
605 contents.coincmaptable.append(coincmap)
606
607 return coinc
608
609
610#
611# =============================================================================
612#
613# Library API
614#
615# =============================================================================
616#
617
618
619def binjfind(xmldoc, process, search, snglcomparefunc, nearcoinccomparefunc, verbose = False):
620 #
621 # Analyze the document's contents.
622 #
623
624 if verbose:
625 print("indexing ...", file=sys.stderr)
626
627 b_b_def = {
628 "StringCusp": burca.StringCuspBBCoincDef,
629 "excesspower": burca.ExcessPowerBBCoincDef,
630 "waveburst": CWBBBCoincDef,
631 "omega": OmegaBBCoincDef
632 }[search]
633 sb_b_def = {
634 "StringCusp": StringCuspSBBCoincDef,
635 "excesspower": ExcessPowerSBBCoincDef,
636 "waveburst": CWBSBBCoincDef,
637 "omega": OmegaSBBCoincDef
638 }[search]
639 si_b_def = {
640 "StringCusp": StringCuspSIBCoincDef,
641 "excesspower": ExcessPowerSIBCoincDef,
642 "waveburst": CWBSIBCoincDef,
643 "omega": OmegaSIBCoincDef
644 }[search]
645 sb_c_e_def = {
646 "StringCusp": StringCuspSBCCoincDef,
647 "excesspower": ExcessPowerSBCCoincDef,
648 "waveburst": CWBSBCCoincDef,
649 "omega": OmegaSBCCoincDef
650 }[search]
651 sb_c_n_def = {
652 "StringCusp": StringCuspSBCNearCoincDef,
653 "excesspower": ExcessPowerSBCNearCoincDef,
654 "waveburst": CWBSBCNearCoincDef,
655 "omega": OmegaSBCNearCoincDef
656 }[search]
657 si_c_e_def = {
658 "StringCusp": StringCuspSICCoincDef,
659 "excesspower": ExcessPowerSICCoincDef,
660 "waveburst": CWBSICCoincDef,
661 "omega": OmegaSICCoincDef
662 }[search]
663 si_c_n_def = {
664 "StringCusp": StringCuspSICNearCoincDef,
665 "excesspower": ExcessPowerSICNearCoincDef,
666 "waveburst": CWBSICNearCoincDef,
667 "omega": OmegaSICNearCoincDef
668 }[search]
669
670 contents = DocContents(
671 xmldoc = xmldoc,
672 b_b_def = b_b_def,
673 sb_b_def = sb_b_def,
674 si_b_def = si_b_def,
675 sb_c_e_def = sb_c_e_def,
676 sb_c_n_def = sb_c_n_def,
677 si_c_e_def = si_c_e_def,
678 si_c_n_def = si_c_n_def,
679 process = process,
680 livetime_program = {
681 "StringCusp": "StringSearch",
682 "excesspower": "lalapps_power",
683 "omega": None, # FIXME: this causes all segments in the search_summary table to be retrieved
684 "waveburst": None # FIXME: this causes all segments in the search_summary table to be retrieved
685 }[search]
686 )
687
688 #
689 # set the window for the sieve in bursts_near_peaktime(). this
690 # window is the amount of time such that if an injection's peak
691 # time and a burst event's peak time differ by more than this it is
692 # *impossible* for them to match one another.
693 #
694
695 # the radius of Earth in light seconds == the most an injection's
696 # peak time column can differ from the time it peaks in an
697 # instrument. 1.25 = add 25% for good luck (we're not being
698 # careful with asphericity here, so a bit of padding is needed
699 # anyway, just make sure it's enough).
700 burst_peak_time_window = lal.REARTH_SI / lal.C_SI * 1.25
701
702 # add the duration of the longest burst event (the most a burst
703 # event's peak time could differ from either the start or stop time
704 # of the event)
705 burst_peak_time_window += contents.longestduration
706
707 # add a search-specific padding
708 #
709 # for the string search this is 1/2 the typical width of a whitened
710 # template's autocorrelation function (because this is added to a
711 # time window somewhere)
712 #
713 # for omega, the number should be larger than the largest delta_t
714 # passed to OmegaSnglCompare() not including the duration of the
715 # burst event (already taken into account above)
716 burst_peak_time_window += {
717 "StringCusp": SimBurstUtils.stringcusp_autocorrelation_width / 2,
718 "excesspower": 0.0,
719 "omega": 30.0,
720 "waveburst": 30.0
721 }[search]
722
723 #
724 # set the window for identifying coincs near a peak time
725 # FIXME: this is kinda specific to the excess power search.
726 #
727
728 coinc_peak_time_window = burst_peak_time_window + SimBurstUtils.burst_is_near_injection_window
729
730 #
731 # Search for sim_burst <--> * coincidences
732 #
733
734 if contents.simbursttable is not None:
735 N = len(contents.simbursttable)
736
737 #
738 # Find sim_burst <--> sngl_burst coincidences.
739 #
740
741 if verbose:
742 print("constructing %s:" % sb_b_def.description, file=sys.stderr)
743 for n, sim in enumerate(contents.simbursttable):
744 if verbose:
745 print("\t%.1f%%\r" % (100.0 * n / N), end=' ', file=sys.stderr)
746 events = find_sngl_burst_matches(contents, sim, snglcomparefunc, burst_peak_time_window)
747 if events:
748 add_sim_burst_coinc(contents, sim, events, contents.sb_b_coinc_def_id)
749 if verbose:
750 print("\t100.0%", file=sys.stderr)
751
752 #
753 # Find sim_burst <--> coinc_event coincidences.
754 #
755
756 if verbose:
757 print("constructing %s and %s:" % (sb_c_e_def.description, sb_c_n_def.description), file=sys.stderr)
758 for n, sim in enumerate(contents.simbursttable):
759 if verbose:
760 print("\t%.1f%%\r" % (100.0 * n / N), end=' ', file=sys.stderr)
761 offsetvector = contents.offsetvectors[sim.time_slide_id]
762 coincs = contents.coincs_near_peaktime(sim.time_geocent, coinc_peak_time_window, offsetvector)
763 exact_coinc_event_ids = find_exact_coinc_matches(coincs, sim, snglcomparefunc, contents.seglists, offsetvector)
764 near_coinc_event_ids = find_near_coinc_matches(coincs, sim, nearcoinccomparefunc, offsetvector)
765 assert exact_coinc_event_ids.issubset(near_coinc_event_ids)
766 if exact_coinc_event_ids:
767 add_sim_coinc_coinc(contents, sim, exact_coinc_event_ids, contents.sb_c_e_coinc_def_id)
768 if near_coinc_event_ids:
769 add_sim_coinc_coinc(contents, sim, near_coinc_event_ids, contents.sb_c_n_coinc_def_id)
770 if verbose:
771 print("\t100.0%", file=sys.stderr)
772 elif verbose:
773 print("no %s table in document, skipping" % lsctables.SimBurstTable.tableName, file=sys.stderr)
774
775 #
776 # Search for sim_inspiral <--> * coincidences
777 #
778
779 if contents.siminspiraltable is not None:
780 N = len(contents.siminspiraltable)
781
782 #
783 # Find sim_inspiral <--> sngl_burst coincidences.
784 #
785
786 if verbose:
787 print("constructing %s:" % si_b_def.description, file=sys.stderr)
788 for n, sim in enumerate(contents.siminspiraltable):
789 if verbose:
790 print("\t%.1f%%\r" % (100.0 * n / N), end=' ', file=sys.stderr)
791 events = find_sngl_burst_matches(contents, sim, snglcomparefunc, burst_peak_time_window)
792 if events:
793 add_sim_burst_coinc(contents, sim, events, contents.si_b_coinc_def_id)
794 if verbose:
795 print("\t100.0%", file=sys.stderr)
796
797 #
798 # Find sim_inspiral <--> coinc_event coincidences.
799 #
800
801 if verbose:
802 print("constructing %s and %s:" % (si_c_e_def.description, si_c_n_def.description), file=sys.stderr)
803 for n, sim in enumerate(contents.siminspiraltable):
804 if verbose:
805 print("\t%.1f%%\r" % (100.0 * n / N), end=' ', file=sys.stderr)
806 offsetvector = contents.offsetvectors[sim.time_slide_id]
807 coincs = contents.coincs_near_peaktime(sim.time_geocent, coinc_peak_time_window, offsetvector)
808 exact_coinc_event_ids = find_exact_coinc_matches(coincs, sim, snglcomparefunc, contents.seglists, offsetvector)
809 near_coinc_event_ids = find_near_coinc_matches(coincs, sim, nearcoinccomparefunc, offsetvector)
810 assert exact_coinc_event_ids.issubset(near_coinc_event_ids)
811 if exact_coinc_event_ids:
812 add_sim_coinc_coinc(contents, sim, exact_coinc_event_ids, contents.si_c_e_coinc_def_id)
813 if near_coinc_event_ids:
814 add_sim_coinc_coinc(contents, sim, near_coinc_event_ids, contents.si_c_n_coinc_def_id)
815 if verbose:
816 print("\t100.0%", file=sys.stderr)
817 elif verbose:
818 print("no %s table in document, skipping" % lsctables.SimInspiralTable.tableName, file=sys.stderr)
819
820 #
821 # Done.
822 #
823
824 return xmldoc
static double max(double a, double b)
Definition: EPFilters.c:43
A wrapper interface to the XML document.
Definition: binjfind.py:165
def new_coinc(self, coinc_def_id, time_slide_id)
Construct a new coinc_event row attached to the given process, and belonging to the set of coincidenc...
Definition: binjfind.py:369
def bursts_near_peaktime(self, t, window, offsetvector)
Return a list of the burst events whose peak times (with offsetvector added) are within window second...
Definition: binjfind.py:336
def __init__(self, xmldoc, b_b_def, sb_b_def, si_b_def, sb_c_e_def, sb_c_n_def, si_c_e_def, si_c_n_def, process, livetime_program)
Definition: binjfind.py:166
def coincs_near_peaktime(self, t, window, offsetvector)
Return a list of the (coinc_event_id, event list) tuples in which at least one burst event's peak tim...
Definition: binjfind.py:349
def __gt__(self, other)
Definition: binjfind.py:98
def __eq__(self, other)
Definition: binjfind.py:89
def __ge__(self, other)
Definition: binjfind.py:95
def __lt__(self, other)
Definition: binjfind.py:83
def __ne__(self, other)
Definition: binjfind.py:92
def __le__(self, other)
Definition: binjfind.py:86
def OmegaSnglCompare(sim, burst, offsetvector, delta_t=10.0)
Return False (injection matches event) if the time of the sim and the peak time of the burst event di...
Definition: binjfind.py:445
def append_process(xmldoc, match_algorithm, comment)
Convenience wrapper for adding process metadata to the document.
Definition: binjfind.py:397
def find_exact_coinc_matches(coincs, sim, comparefunc, seglists, offsetvector)
Return a set of the coinc_event_ids of the burst<-->burst coincs in which all burst events match sim ...
Definition: binjfind.py:555
def StringCuspNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
Definition: binjfind.py:461
def binjfind(xmldoc, process, search, snglcomparefunc, nearcoinccomparefunc, verbose=False)
Definition: binjfind.py:619
def find_sngl_burst_matches(contents, sim, comparefunc, sieve_window)
Scan the burst table for triggers matching sim.
Definition: binjfind.py:509
def StringCuspSnglCompare(sim, burst, offsetvector)
Return False (injection matches event) if an autocorrelation-width window centred on the injection is...
Definition: binjfind.py:425
def add_sim_coinc_coinc(contents, sim, coinc_event_ids, coinc_def_id)
Create a coinc_event in the coinc table, and add arcs in the coinc_event_map table linking the sim_bu...
Definition: binjfind.py:589
def ExcessPowerSnglCompare(sim, burst, offsetvector)
Return False (injection matches event) if the peak time and centre frequency of sim lie within the ti...
Definition: binjfind.py:436
def CWBSnglCompare(sim, burst, offsetvector, delta_t=10.0)
Return False (injection matches event) if the time of the sim and the peak time of the burst event di...
Definition: binjfind.py:453
def ExcessPowerNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
Definition: binjfind.py:471
def find_near_coinc_matches(coincs, sim, comparefunc, offsetvector)
Return a set of the coinc_event_ids of the burst<-->burst coincs in which at least one burst event ma...
Definition: binjfind.py:572
def OmegaNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
Definition: binjfind.py:481
def CWBNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
Definition: binjfind.py:488
def add_sim_burst_coinc(contents, sim, events, coinc_def_id)
Create a coinc_event in the coinc table, and add arcs in the coinc_event_map table linking the sim_bu...
Definition: binjfind.py:519