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.
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
49import igwn_segments as segments
51from . import SimBurstUtils
54__author__ = "Kipp Cannon <kipp.cannon@ligo.org>"
55from .git_version import date as __date__
56from .git_version import version as __version__
60# =============================================================================
84 return self.
end < other
87 return self.
end <= other
90 return self.
end == other
93 return self.
end != other
96 return self.
end >= other
99 return self.
end > other
102lsctables.SnglBurstTable.RowType = lsctables.SnglBurst = SnglBurst
117lsctables.SimInspiralTable.RowType = lsctables.SimInspiral = SimInspiral
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)")
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)")
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)")
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)")
164 A wrapper interface to the XML document.
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):
187 timeslidetable = lsctables.TimeSlideTable.get_table(xmldoc)
189 timeslidetable =
None
190 if timeslidetable
is not None:
209 self.
seglists = ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, livetime_program).coalesce()
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)
224 sim.time_slide_id = time_slide_id
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)
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)
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)
251 b_b_coinc_def_id =
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)
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)
277 self.
coinctable = lsctables.New(lsctables.CoincTable)
278 xmldoc.childNodes[0].appendChild(self.
coinctable)
298 self.
coincs =
dict((row.coinc_event_id, [])
for row
in self.
coinctable if row.coinc_def_id == b_b_coinc_def_id)
302 self.
coincs[row.coinc_event_id].
append(index[row.event_id])
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)
327 self.
coincs.sort(key =
lambda coinc_id_events: coinc_id_events[1][0].peak)
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
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()), [])
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.
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])]
364 def new_coinc(self, coinc_def_id, time_slide_id):
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.
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)
377 coinc.likelihood =
None
391process_program_name =
"lalburst_injfind"
396 Convenience wrapper for adding process metadata to the document.
398 return ligolw_process.register_to_xmldoc(
400 program = process_program_name,
402 "match_algorithm": match_algorithm
404 version = __version__,
405 cvs_repository =
"lscsoft",
406 cvs_entry_time = __date__,
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.
426 tinj = sim.time_at_instrument(burst.ifo, offsetvector)
427 window = SimBurstUtils.stringcusp_autocorrelation_width / 2
429 return segments.segment(tinj - window, tinj + window).disjoint(burst.period)
434 Return False (injection matches event) if the peak time and centre
435 frequency of sim lie within the time-frequency tile of burst.
437 return (sim.time_at_instrument(burst.ifo, offsetvector)
not in burst.period)
or (sim.frequency
not in burst.band)
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
446 return abs(float(sim.time_at_instrument(burst.ifo, offsetvector) - burst.peak)) > delta_t
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
454 return abs(float(sim.time_at_instrument(burst.ifo, offsetvector) - burst.peak)) > delta_t
459 Return False (injection matches coinc) if the peak time of the sim
460 is "near" the burst event.
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)
469 Return False (injection matches coinc) if the peak time of the sim
470 is "near" the burst event.
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)
479 Return False (injection matches coinc) if the peak time of the sim
480 is "near" the burst event.
482 return OmegaSnglCompare(sim, burst, offsetvector, delta_t = 20.0 + burst.duration / 2)
486 Return False (injection matches coinc) if the peak time of the sim
487 is "near" the burst event.
489 return OmegaSnglCompare(sim, burst, offsetvector, delta_t = 20.0 + burst.duration / 2)
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.
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)]
514def add_sim_burst_coinc(contents, sim, events, coinc_def_id):
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.
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)
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)
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)
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.
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))
570 Return a set of the coinc_event_ids of the burst<-->burst coincs in
571 which at least one burst 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))
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.
591 coinc = contents.new_coinc(coinc_def_id, sim.time_slide_id)
592 coinc.nevents = len(coinc_event_ids)
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)
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)
619def binjfind(xmldoc, process, search, snglcomparefunc, nearcoinccomparefunc, verbose = False):
625 print(
"indexing ...", file=sys.stderr)
628 "StringCusp": burca.StringCuspBBCoincDef,
629 "excesspower": burca.ExcessPowerBBCoincDef,
630 "waveburst": CWBBBCoincDef,
631 "omega": OmegaBBCoincDef
634 "StringCusp": StringCuspSBBCoincDef,
635 "excesspower": ExcessPowerSBBCoincDef,
636 "waveburst": CWBSBBCoincDef,
637 "omega": OmegaSBBCoincDef
640 "StringCusp": StringCuspSIBCoincDef,
641 "excesspower": ExcessPowerSIBCoincDef,
642 "waveburst": CWBSIBCoincDef,
643 "omega": OmegaSIBCoincDef
646 "StringCusp": StringCuspSBCCoincDef,
647 "excesspower": ExcessPowerSBCCoincDef,
648 "waveburst": CWBSBCCoincDef,
649 "omega": OmegaSBCCoincDef
652 "StringCusp": StringCuspSBCNearCoincDef,
653 "excesspower": ExcessPowerSBCNearCoincDef,
654 "waveburst": CWBSBCNearCoincDef,
655 "omega": OmegaSBCNearCoincDef
658 "StringCusp": StringCuspSICCoincDef,
659 "excesspower": ExcessPowerSICCoincDef,
660 "waveburst": CWBSICCoincDef,
661 "omega": OmegaSICCoincDef
664 "StringCusp": StringCuspSICNearCoincDef,
665 "excesspower": ExcessPowerSICNearCoincDef,
666 "waveburst": CWBSICNearCoincDef,
667 "omega": OmegaSICNearCoincDef
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,
681 "StringCusp":
"StringSearch",
682 "excesspower":
"lalapps_power",
700 burst_peak_time_window = lal.REARTH_SI / lal.C_SI * 1.25
705 burst_peak_time_window += contents.longestduration
716 burst_peak_time_window += {
717 "StringCusp": SimBurstUtils.stringcusp_autocorrelation_width / 2,
728 coinc_peak_time_window = burst_peak_time_window + SimBurstUtils.burst_is_near_injection_window
734 if contents.simbursttable
is not None:
735 N = len(contents.simbursttable)
742 print(
"constructing %s:" % sb_b_def.description, file=sys.stderr)
743 for n, sim
in enumerate(contents.simbursttable):
745 print(
"\t%.1f%%\r" % (100.0 * n / N), end=
' ', file=sys.stderr)
750 print(
"\t100.0%", file=sys.stderr)
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):
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)
765 assert exact_coinc_event_ids.issubset(near_coinc_event_ids)
766 if exact_coinc_event_ids:
768 if near_coinc_event_ids:
771 print(
"\t100.0%", file=sys.stderr)
773 print(
"no %s table in document, skipping" % lsctables.SimBurstTable.tableName, file=sys.stderr)
779 if contents.siminspiraltable
is not None:
780 N = len(contents.siminspiraltable)
787 print(
"constructing %s:" % si_b_def.description, file=sys.stderr)
788 for n, sim
in enumerate(contents.siminspiraltable):
790 print(
"\t%.1f%%\r" % (100.0 * n / N), end=
' ', file=sys.stderr)
795 print(
"\t100.0%", file=sys.stderr)
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):
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)
810 assert exact_coinc_event_ids.issubset(near_coinc_event_ids)
811 if exact_coinc_event_ids:
813 if near_coinc_event_ids:
816 print(
"\t100.0%", file=sys.stderr)
818 print(
"no %s table in document, skipping" % lsctables.SimInspiralTable.tableName, file=sys.stderr)
static double max(double a, double b)
A wrapper interface to the XML document.
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...
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...
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)
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...
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...
def append_process(xmldoc, match_algorithm, comment)
Convenience wrapper for adding process metadata to the document.
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 ...
def StringCuspNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
def binjfind(xmldoc, process, search, snglcomparefunc, nearcoinccomparefunc, verbose=False)
def find_sngl_burst_matches(contents, sim, comparefunc, sieve_window)
Scan the burst table for triggers matching sim.
def StringCuspSnglCompare(sim, burst, offsetvector)
Return False (injection matches event) if an autocorrelation-width window centred on the injection is...
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...
def ExcessPowerSnglCompare(sim, burst, offsetvector)
Return False (injection matches event) if the peak time and centre frequency of sim lie within the ti...
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...
def ExcessPowerNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
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...
def OmegaNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
def CWBNearCoincCompare(sim, burst, offsetvector)
Return False (injection matches coinc) if the peak time of the sim is "near" the burst event.
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...