pylal.imr

1 # Copyright (C) 2012 Chad Hanna 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 import sys 18 from glue.ligolw import lsctables 19 from glue.ligolw import dbtables 20 from glue import segments 21 from glue import segmentsUtils 22 from glue.ligolw import table 23 from pylal import db_thinca_rings 24 from lal import rate 25 import numpy 26 import math 27 import copy 28 from glue.ligolw.utils import search_summary as ligolw_search_summary 29 from glue.ligolw.utils import segments as ligolw_segments 30 from glue.ligolw.utils import process 31 from lalsimulation import SimInspiralTaylorF2ReducedSpinComputeChi, SimIMRPhenomBComputeChi 32 33 import sqlite3 34 35

36 -def allowed_analysis_table_names():

37 return (dbtables.lsctables.MultiBurstTable.tableName, dbtables.lsctables.CoincInspiralTable.tableName, dbtables.lsctables.CoincRingdownTable.tableName)

38 39

40 -def make_sim_inspiral_row_from_columns_in_db(connection):

41 """ 42 get the unique mapping of a sim inspiral row from columns in this 43 database 44 """ 45 return lsctables.SimInspiralTable.get_table(dbtables.get_xml(connection)).row_from_cols

46 47

48 -def time_within_segments(geocent_end_time, geocent_end_time_ns, zero_lag_segments = None):

49 """ 50 Return True if injection was made in the given segmentlist, if no 51 segments just return True 52 """ 53 if zero_lag_segments is None: 54 return True 55 else: 56 return lsctables.LIGOTimeGPS(geocent_end_time, geocent_end_time_ns) in zero_lag_segments

57 58

59 -def get_min_far_inspiral_injections(connection, segments = None, table_name = "coinc_inspiral"):

60 """ 61 This function returns the found injections from a database and the 62 minimum far associated with them as tuple of the form (far, sim). It also tells 63 you all of the injections that should have been injected. Subtracting the two 64 outputs should tell you the missed injections 65 """ 66 67 if table_name == dbtables.lsctables.CoincInspiralTable.tableName: 68 found_query = 'SELECT sim_inspiral.*, coinc_inspiral.combined_far FROM sim_inspiral JOIN coinc_event_map AS mapA ON mapA.event_id == sim_inspiral.simulation_id JOIN coinc_event_map AS mapB ON mapB.coinc_event_id == mapA.coinc_event_id JOIN coinc_inspiral ON coinc_inspiral.coinc_event_id == mapB.event_id JOIN coinc_event on coinc_event.coinc_event_id == coinc_inspiral.coinc_event_id WHERE mapA.table_name = "sim_inspiral" AND mapB.table_name = "coinc_event" AND injection_in_segments(sim_inspiral.geocent_end_time, sim_inspiral.geocent_end_time_ns)' 69 70 elif table_name == dbtables.lsctables.CoincRingdownTable.tableName: 71 found_query = 'SELECT sim_inspiral.*, coinc_ringdown.false_alarm_rate FROM sim_inspiral JOIN coinc_event_map AS mapA ON mapA.event_id == sim_inspiral.simulation_id JOIN coinc_event_map AS mapB ON mapB.coinc_event_id == mapA.coinc_event_id JOIN coinc_ringdown ON coinc_ringdown.coinc_event_id == mapB.event_id JOIN coinc_event on coinc_event.coinc_event_id == coinc_ringdown.coinc_event_id WHERE mapA.table_name = "sim_inspiral" AND mapB.table_name = "coinc_event" AND injection_in_segments(sim_inspiral.geocent_end_time, sim_inspiral.geocent_end_time_ns)' 72 73 elif table_name == dbtables.lsctables.MultiBurstTable.tableName: 74 found_query = 'SELECT sim_inspiral.*, multi_burst.false_alarm_rate FROM sim_inspiral JOIN coinc_event_map AS mapA ON mapA.event_id == sim_inspiral.simulation_id JOIN coinc_event_map AS mapB ON mapB.coinc_event_id == mapA.coinc_event_id JOIN multi_burst ON multi_burst.coinc_event_id == mapB.event_id JOIN coinc_event on coinc_event.coinc_event_id == multi_burst.coinc_event_id WHERE mapA.table_name = "sim_inspiral" AND mapB.table_name = "coinc_event" AND injection_in_segments(sim_inspiral.geocent_end_time, sim_inspiral.geocent_end_time_ns)' 75 76 else: 77 raise ValueError("table must be in " + " ".join(allowed_analysis_table_names())) 78 79 def injection_was_made(end_time, end_time_ns, segments = segments): 80 return time_within_segments(end_time, end_time_ns, segments)

81 82 # restrict the found injections to only be within certain segments 83 connection.create_function("injection_in_segments", 2, injection_was_made) 84 85 # get the mapping of a record returned by the database to a sim 86 # inspiral row. Note that this is DB dependent potentially, so always 87 # do this! 88 make_sim_inspiral = make_sim_inspiral_row_from_columns_in_db(connection) 89 90 found_injections = {} 91 92 for values in connection.cursor().execute(found_query): 93 # all but the last column is used to build a sim inspiral object 94 sim = make_sim_inspiral(values[:-1]) 95 far = values[-1] 96 # update with the minimum far seen until now 97 this_inj = found_injections.setdefault(sim.simulation_id, (far, sim)) 98 if far < this_inj[0]: 99 found_injections[sim.simulation_id] = (far, sim) 100 101 total_query = 'SELECT * FROM sim_inspiral WHERE injection_in_segments(geocent_end_time, geocent_end_time_ns)' 102 103 total_injections = {} 104 # Missed injections start as a copy of the found injections 105 missed_injections = {} 106 for values in connection.cursor().execute(total_query): 107 sim = make_sim_inspiral(values) 108 total_injections[sim.simulation_id] = sim 109 missed_injections[sim.simulation_id] = sim 110 111 # now actually remove the missed injections 112 for k in found_injections: 113 del missed_injections[k] 114 115 116 return found_injections.values(), total_injections.values(), missed_injections.values() 117 118

119 -def get_max_snr_inspiral_injections(connection, segments = None, table_name = "coinc_inspiral"):

120 """ 121 Like get_min_far_inspiral_injections but uses SNR to rank injections. 122 """ 123 124 if table_name == dbtables.lsctables.CoincInspiralTable.tableName: 125 found_query = 'SELECT sim_inspiral.*, coinc_inspiral.snr FROM sim_inspiral JOIN coinc_event_map AS mapA ON mapA.event_id == sim_inspiral.simulation_id JOIN coinc_event_map AS mapB ON mapB.coinc_event_id == mapA.coinc_event_id JOIN coinc_inspiral ON coinc_inspiral.coinc_event_id == mapB.event_id JOIN coinc_event on coinc_event.coinc_event_id == coinc_inspiral.coinc_event_id WHERE mapA.table_name = "sim_inspiral" AND mapB.table_name = "coinc_event" AND injection_in_segments(sim_inspiral.geocent_end_time, sim_inspiral.geocent_end_time_ns)' 126 127 elif table_name in allowed_analysis_table_names(): 128 raise NotImplementedError("get_max_snr_inspiral_injections has not yet implemented querying against the table %s. Please consider submitting a patch. See get_min_far_inspiral_injections for how to construct your query." % table_name) 129 else: 130 raise ValueError("table must be in " + " ".join(allowed_analysis_table_names())) 131 132 133 def injection_was_made(end_time, end_time_ns, segments = segments): 134 return time_within_segments(end_time, end_time_ns, segments)

135 136 # restrict the found injections to only be within certain segments 137 connection.create_function("injection_in_segments", 2, injection_was_made) 138 139 # get the mapping of a record returned by the database to a sim 140 # inspiral row. Note that this is DB dependent potentially, so always 141 # do this! 142 make_sim_inspiral = make_sim_inspiral_row_from_columns_in_db(connection) 143 144 found_injections = {} 145 146 for values in connection.cursor().execute(found_query): 147 # all but the last column is used to build a sim inspiral object 148 sim = make_sim_inspiral(values[:-1]) 149 snr = values[-1] 150 # update with the minimum far seen until now 151 this_inj = found_injections.setdefault(sim.simulation_id, (snr, sim)) 152 if snr > this_inj[0]: 153 found_injections[sim.simulation_id] = (snr, sim) 154 155 total_query = 'SELECT * FROM sim_inspiral WHERE injection_in_segments(geocent_end_time, geocent_end_time_ns)' 156 157 total_injections = {} 158 # Missed injections start as a copy of the found injections 159 missed_injections = {} 160 for values in connection.cursor().execute(total_query): 161 sim = make_sim_inspiral(values) 162 total_injections[sim.simulation_id] = sim 163 missed_injections[sim.simulation_id] = sim 164 165 # now actually remove the missed injections 166 for k in found_injections: 167 del missed_injections[k] 168 169 return found_injections.values(), total_injections.values(), missed_injections.values() 170 171

172 -def get_instruments_from_coinc_event_table(connection):

173 """ 174 This function returns a list of the instruments analyzed according to the coinc_event_table 175 """ 176 instruments = [] 177 for ifos in connection.cursor().execute('SELECT DISTINCT(instruments) FROM coinc_event WHERE instruments!=""'): 178 # ignore null columns 179 if ifos[0]: 180 instruments.append(frozenset(lsctables.instrument_set_from_ifos(ifos[0]))) 181 return instruments

182 183

184 -def get_segments(connection, xmldoc, table_name, live_time_program, veto_segments_name = None, data_segments_name = "datasegments"):

185 segs = segments.segmentlistdict() 186 187 if table_name == dbtables.lsctables.CoincInspiralTable.tableName: 188 if live_time_program == "gstlal_inspiral": 189 segs = ligolw_segments.segmenttable_get_by_name(xmldoc, data_segments_name).coalesce() 190 segs &= ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, live_time_program).coalesce() 191 elif live_time_program == "thinca": 192 segs = db_thinca_rings.get_thinca_zero_lag_segments(connection, program_name = live_time_program).coalesce() 193 else: 194 raise ValueError("for burst tables livetime program must be one of gstlal_inspiral, thinca") 195 if veto_segments_name is not None: 196 veto_segs = db_thinca_rings.get_veto_segments(connection, veto_segments_name) 197 segs -= veto_segs 198 return segs 199 elif table_name == dbtables.lsctables.CoincRingdownTable.tableName: 200 segs = ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, live_time_program).coalesce() 201 if veto_segments_name is not None: 202 veto_segs = ligolw_segments.segmenttable_get_by_name(xmldoc, veto_segments_name).coalesce() 203 segs -= veto_segs 204 return segs 205 elif table_name == dbtables.lsctables.MultiBurstTable.tableName: 206 if live_time_program == "omega_to_coinc": 207 segs = ligolw_search_summary.segmentlistdict_fromsearchsummary(xmldoc, live_time_program).coalesce() 208 if veto_segments_name is not None: 209 veto_segs = ligolw_segments.segmenttable_get_by_name(xmldoc, veto_segments_name).coalesce() 210 segs -= veto_segs 211 elif live_time_program == "waveburst": 212 segs = db_thinca_rings.get_thinca_zero_lag_segments(connection, program_name = live_time_program).coalesce() 213 if veto_segments_name is not None: 214 veto_segs = db_thinca_rings.get_veto_segments(connection, veto_segments_name) 215 segs -= veto_segs 216 else: 217 raise ValueError("for burst tables livetime program must be one of omega_to_coinc, waveburst") 218 return segs 219 else: 220 raise ValueError("table must be in " + " ".join(allowed_analysis_table_names()))

221 222

223 -def get_event_fars(connection, table_name, segments = None):

224 """ 225 return the false alarm rate of the most rare zero-lag coinc by instruments 226 """ 227 228 def event_in_requested_segments(end_time, end_time_ns, segments = segments): 229 return time_within_segments(end_time, end_time_ns, segments)

230 231 connection.create_function("event_in_requested_segments", 2, event_in_requested_segments) 232 233 if table_name == dbtables.lsctables.CoincInspiralTable.tableName: 234 query = 'SELECT coinc_event.instruments, coinc_inspiral.combined_far AS combined_far, EXISTS(SELECT * FROM time_slide WHERE time_slide.time_slide_id == coinc_event.time_slide_id AND time_slide.offset != 0) FROM coinc_inspiral JOIN coinc_event ON (coinc_inspiral.coinc_event_id == coinc_event.coinc_event_id) WHERE event_in_requested_segments(coinc_inspiral.end_time, coinc_inspiral.end_time_ns);' 235 236 elif table_name == dbtables.lsctables.MultiBurstTable.tableName: 237 query = 'SELECT coinc_event.instruments, multi_burst.false_alarm_rate AS combined_far, EXISTS(SELECT * FROM time_slide WHERE time_slide.time_slide_id == coinc_event.time_slide_id AND time_slide.offset != 0) FROM multi_burst JOIN coinc_event ON (multi_burst.coinc_event_id == coinc_event.coinc_event_id) WHERE event_in_requested_segments(multi_burst.peak_time, multi_burst.peak_time_ns);' 238 239 elif table_name == dbtables.lsctables.CoincRingdownTable.tableName: 240 query = 'SELECT coinc_event.instruments, coinc_ringdown.false_alarm_rate AS combined_far, EXISTS(SELECT * FROM time_slide WHERE time_slide.time_slide_id == coinc_event.time_slide_id AND time_slide.offset != 0) FROM coinc_ringdown JOIN coinc_event ON (coinc_ringdown.coinc_event_id == coinc_event.coinc_event_id) WHERE event_in_requested_segments(coinc_ringdown.start_time, coinc_ringdown.start_time_ns);' 241 242 else: 243 raise ValueError("table must be in " + " ".join(allowed_analysis_table_names())) 244 245 for inst, far, ts in connection.cursor().execute(query): 246 inst = frozenset(lsctables.instrument_set_from_ifos(inst)) 247 yield (inst, far, ts) 248 249

250 -def compute_search_efficiency_in_bins(found, total, ndbins, sim_to_bins_function = lambda sim: (sim.distance,)):

251 """ 252 This program creates the search efficiency in the provided ndbins. The 253 first dimension of ndbins must be the distance. You also must provide a 254 function that maps a sim inspiral row to the correct tuple to index the ndbins. 255 """ 256 257 num = rate.BinnedArray(ndbins) 258 den = rate.BinnedArray(ndbins) 259 260 # increment the numerator with the found injections 261 for sim in found: 262 num[sim_to_bins_function(sim)] += 1 263 264 # increment the denominator with the total injections 265 for sim in total: 266 den[sim_to_bins_function(sim)] += 1 267 268 # sanity check 269 assert (num.array <= den.array).all(), "some bins have more found injections than were made" 270 271 # regularize by setting empty bins to zero efficiency 272 den.array[numpy.logical_and(num.array == 0, den.array == 0)] = 1 273 274 # pull out the efficiency array, it is the ratio 275 eff = rate.BinnedArray(rate.NDBins(ndbins), array = num.array / den.array) 276 277 # compute binomial uncertainties in each bin 278 k = num.array 279 N = den.array 280 eff_lo_arr = ( N*(2*k + 1) - numpy.sqrt(4*N*k*(N - k) + N**2) ) / (2*N*(N + 1)) 281 eff_hi_arr = ( N*(2*k + 1) + numpy.sqrt(4*N*k*(N - k) + N**2) ) / (2*N*(N + 1)) 282 283 eff_lo = rate.BinnedArray(rate.NDBins(ndbins), array = eff_lo_arr) 284 eff_hi = rate.BinnedArray(rate.NDBins(ndbins), array = eff_hi_arr) 285 286 return eff_lo, eff, eff_hi

287 288

289 -def compute_search_volume(eff):

290 """ 291 Integrate efficiency to get search volume. 292 """ 293 # get distance bins 294 ndbins = eff.bins 295 dx = ndbins[0].upper() - ndbins[0].lower() 296 r = ndbins[0].centres() 297 298 # we have one less dimension on the output 299 vol = rate.BinnedArray(rate.NDBins(ndbins[1:])) 300 301 # integrate efficiency to obtain volume 302 vol.array = numpy.trapz(eff.array.T * 4. * numpy.pi * r**2, r, dx) 303 304 return vol

305 306

307 -def guess_nd_bins(sims, bin_dict = {"distance": (200, rate.LinearBins)}):

308 """ 309 Given a dictionary of bin counts and bin objects keyed by sim 310 attribute, come up with a sensible NDBins scheme 311 """ 312 return rate.NDBins([bintup[1](min([getattr(sim, attr) for sim in sims]), max([getattr(sim, attr) for sim in sims]) + sys.float_info.min, bintup[0]) for attr, bintup in bin_dict.items()])

313 314

315 -def guess_distance_mass1_mass2_bins_from_sims(sims, mass1bins = 11, mass2bins = 11, distbins = 200):

316 """ 317 Given a list of the injections, guess at the mass1, mass2 and distance 318 bins. 319 """ 320 return guess_nd_bins(sims, bin_dict = {"distance": (distbins, rate.LinearBins), "mass1": (mass1bins, rate.LinearBins), "mass2": (mass2bins, rate.LinearBins)})

321 322

323 -def guess_distance_spin1z_spin2z_bins_from_sims(sims, spin1bins = 11, spin2bins = 11, distbins = 200):

324 """ 325 Given a list of the injections, guess at the spin1, spin2 and distance 326 bins. 327 """ 328 return guess_nd_bins(sims, bin_dict = {"distance": (distbins, rate.LinearBins), "spin1z": (spin1bins, rate.LinearBins), "spin2z": (spin2bins, rate.LinearBins)})

329 330

331 -def guess_distance_effective_spin_parameter_bins_from_sims(sims, chibins = 11, distbins = 200):

332 """ 333 Given a list of the injections, guess at the chi = (m1*s1z + 334 m2*s2z)/(m1+m2) and distance bins. 335 """ 336 dist_chi_vals = map(sim_to_distance_effective_spin_parameter_bins_function, sims) 337 338 distances = [tup[0] for tup in dist_chi_vals] 339 chis = [tup[1] for tup in dist_chi_vals] 340 341 return rate.NDBins([rate.LinearBins(min(distances), max(distances), distbins), rate.LinearBins(min(chis), max(chis), chibins)])

342 343

344 -def guess_distance_mass_ratio_bins_from_sims(sims, qbins = 11, distbins = 200):

345 """ 346 Given a list of the injections, guess at the chi and distance 347 bins. 348 """ 349 dist_mratio_vals = map(sim_to_distance_mass_ratio_bins_function, sims) 350 351 distances = [tup[0] for tup in dist_mratio_vals] 352 mratios = [tup[1] for tup in dist_mratio_vals] 353 354 return rate.NDBins([rate.LinearBins(min(distances), max(distances), distbins), rate.LinearBins(min(mratios), max(mratios), qbins)])

355 356

357 -def guess_distance_chirp_mass_bins_from_sims(sims, mbins = 11, distbins = 200):

358 """ 359 Given a list of the injections, guess at the chirp mass and distance 360 bins. 361 """ 362 dist_mchirp_vals = map(sim_to_distance_chirp_mass_bins_function, sims) 363 364 distances = [tup[0] for tup in dist_mchirp_vals] 365 mchirps = [tup[1] for tup in dist_mchirp_vals] 366 367 return rate.NDBins([rate.LinearBins(min(distances), max(distances), distbins), rate.LinearBins(min(mchirps), max(mchirps), mbins)])

368 369

370 -def guess_distance_total_mass_bins_from_sims(sims, nbins = 11, distbins = 200):

371 """ 372 Given a list of the injections, guess at the mass1, mass2 and distance 373 bins. Floor and ceil will be used to round down to the nearest integers. 374 """ 375 376 total_lo = numpy.floor(min([sim.mass1 + sim.mass2 for sim in sims])) 377 total_hi = numpy.ceil(max([sim.mass1 + sim.mass2 for sim in sims])) 378 mindist = numpy.floor(min([sim.distance for sim in sims])) 379 maxdist = numpy.ceil(max([sim.distance for sim in sims])) 380 381 return rate.NDBins((rate.LinearBins(mindist, maxdist, distbins), rate.LinearBins(total_lo, total_hi, nbins)))

382 383

384 -def sim_to_distance_mass1_mass2_bins_function(sim):

385 """ 386 create a function to map a sim to a distance, mass1, mass2 NDBins based object 387 """ 388 389 return (sim.distance, sim.mass1, sim.mass2)

390 391

392 -def sim_to_distance_total_mass_bins_function(sim):

393 """ 394 create a function to map a sim to a distance, total mass NDBins based object 395 """ 396 return (sim.distance, sim.mass1 + sim.mass2)

397 398

399 -def sim_to_distance_spin1z_spin2z_bins_function(sim):

400 """ 401 create a function to map a sim to a distance, spin1z, spin2z NDBins based object 402 """ 403 404 return (sim.distance, sim.spin1z, sim.spin2z)

405 406

407 -def sim_to_distance_effective_spin_parameter_bins_function(sim):

408 """ 409 Map a sim_inspiral row to a distance, "chi" spin parameter 410 bin. For IMR waveforms, "chi" refers to the effective spin, 411 412 chi = (m1*s1z + m2*s2z)/(m1 + m2) 413 414 where s1z, s2z are the components of the spins along the 415 direction of the total angular momentum. For inspiral 416 waveforms, "chi" refers to the reduced spin, 417 418 chi_red = chi_s + delta*chi_a - 76.*eta/113*chi_s, 419 420 where chi_s and chi_a are the symmetric and anti-symmetric 421 combinations of the spins, and delta=(m1-m2)/(m1+m2). Some 422 waveforms, e.g., SpinTaylorT4, use different coordinate 423 conventions and require a coordinate transformation before 424 applying these definitions. 425 """ 426 427 if sim.waveform.startswith("SpinTaylorT4"): 428 chi1 = sim.spin1x * math.sin(sim.inclination) + sim.spin1z * math.cos(sim.inclination) 429 chi2 = sim.spin2x * math.sin(sim.inclination) + sim.spin2z * math.cos(sim.inclination) 430 chi = SimInspiralTaylorF2ReducedSpinComputeChi(sim.mass1, sim.mass2, chi1, chi2) 431 432 elif sim.waveform.startswith("SpinTaylorT5"): 433 chi1 = sim.spin1z 434 chi2 = sim.spin2z 435 chi = SimInspiralTaylorF2ReducedSpinComputeChi(sim.mass1, sim.mass2, chi1, chi2) 436 437 elif sim.waveform.startswith("IMRPhenomB") or sim.waveform.startswith("IMRPhenomC") or sim.waveform.startswith("SEOBNR"): 438 chi = SimIMRPhenomBComputeChi(sim.mass1, sim.mass2, sim.spin1z, sim.spin2z) 439 440 else: 441 raise ValueError(sim.waveform) 442 443 return (sim.distance, chi)

444 445

446 -def sim_to_distance_mass_ratio_bins_function(sim):

447 """ 448 create a function to map a sim to a distance, mass ratio NDBins based object 449 """ 450 # note that if you use symmetrize_sims() below, m2/m1 > 1 451 # which just strikes me as more intuitive 452 return (sim.distance, sim.mass2/sim.mass1)

453 454

455 -def sim_to_distance_chirp_mass_bins_function(sim):

456 """ 457 create a function to map a sim to a distance, chirp mass NDBins based object 458 """ 459 return (sim.distance, sim.mchirp)

460

461 -def symmetrize_sims(sims, col1, col2):

462 """ 463 symmetrize by two columns that should be symmetric. For example mass1 and mass2 464 """ 465 for sim in sims: 466 c1 = getattr(sim, col1) 467 c2 = getattr(sim, col2) 468 if c1 > c2: 469 setattr(sim, col1, c2) 470 setattr(sim, col2, c1) 471 return sims

472

473 -class DataBaseSummary(object):

474 """ 475 This class stores summary information gathered across the databases 476 """ 477

478 - def __init__(self, filelist, live_time_program = None, veto_segments_name = None, data_segments_name = "datasegments", tmp_path = None, verbose = False):

479 480 self.segments = segments.segmentlistdict() 481 self.instruments = set() 482 self.table_name = None 483 self.found_injections_by_instrument_set = {} 484 self.missed_injections_by_instrument_set = {} 485 self.total_injections_by_instrument_set = {} 486 self.zerolag_fars_by_instrument_set = {} 487 self.ts_fars_by_instrument_set = {} 488 self.numslides = set() 489 490 for f in filelist: 491 if verbose: 492 print >> sys.stderr, "Gathering stats from: %s...." % (f,) 493 working_filename = dbtables.get_connection_filename(f, tmp_path = tmp_path, verbose = verbose) 494 connection = sqlite3.connect(working_filename) 495 xmldoc = dbtables.get_xml(connection) 496 497 sim = False 498 499 # look for a sim inspiral table. This is IMR work we have to have one of these :) 500 try: 501 sim_inspiral_table = table.get_table(xmldoc, dbtables.lsctables.SimInspiralTable.tableName) 502 sim = True 503 except ValueError: 504 pass 505 506 # look for the relevant table for analyses 507 for table_name in allowed_analysis_table_names(): 508 try: 509 setattr(self, table_name, table.get_table(xmldoc, table_name)) 510 if self.table_name is None or self.table_name == table_name: 511 self.table_name = table_name 512 else: 513 raise ValueError("detected more than one table type out of " + " ".join(allowed_analysis_table_names())) 514 except ValueError: 515 setattr(self, table_name, None) 516 517 # the non simulation databases are where we get information about segments 518 if not sim: 519 self.numslides.add(connection.cursor().execute('SELECT count(DISTINCT(time_slide_id)) FROM time_slide').fetchone()[0]) 520 [self.instruments.add(ifos) for ifos in get_instruments_from_coinc_event_table(connection)] 521 # save a reference to the segments for this file, needed to figure out the missed and found injections 522 self.this_segments = get_segments(connection, xmldoc, self.table_name, live_time_program, veto_segments_name, data_segments_name = data_segments_name) 523 # FIXME we don't really have any reason to use playground segments, but I put this here as a reminder 524 # self.this_playground_segments = segmentsUtils.S2playground(self.this_segments.extent_all()) 525 self.segments += self.this_segments 526 527 # get the far thresholds for the loudest events in these databases 528 for (instruments_set, far, ts) in get_event_fars(connection, self.table_name): 529 if not ts: 530 self.zerolag_fars_by_instrument_set.setdefault(instruments_set, []).append(far) 531 else: 532 self.ts_fars_by_instrument_set.setdefault(instruments_set, []).append(far) 533 # get the injections 534 else: 535 # We need to know the segments in this file to determine which injections are found 536 self.this_injection_segments = get_segments(connection, xmldoc, self.table_name, live_time_program, veto_segments_name, data_segments_name = data_segments_name) 537 self.this_injection_instruments = [] 538 distinct_instruments = connection.cursor().execute('SELECT DISTINCT(instruments) FROM coinc_event WHERE instruments!=""').fetchall() 539 for instruments, in distinct_instruments: 540 instruments_set = frozenset(lsctables.instrument_set_from_ifos(instruments)) 541 self.this_injection_instruments.append(instruments_set) 542 segments_to_consider_for_these_injections = self.this_injection_segments.intersection(instruments_set) - self.this_injection_segments.union(set(self.this_injection_segments.keys()) - instruments_set) 543 found, total, missed = get_min_far_inspiral_injections(connection, segments = segments_to_consider_for_these_injections, table_name = self.table_name) 544 if verbose: 545 print >> sys.stderr, "%s total injections: %d; Found injections %d: Missed injections %d" % (instruments, len(total), len(found), len(missed)) 546 self.found_injections_by_instrument_set.setdefault(instruments_set, []).extend(found) 547 self.total_injections_by_instrument_set.setdefault(instruments_set, []).extend(total) 548 self.missed_injections_by_instrument_set.setdefault(instruments_set, []).extend(missed) 549 550 # All done 551 dbtables.discard_connection_filename(f, working_filename, verbose = verbose) 552 if len(self.numslides) > 1: 553 raise ValueError('number of slides differs between input files') 554 elif self.numslides: 555 self.numslides = min(self.numslides) 556 else: 557 self.numslides = 0

558 559 # FIXME 560 # Things left to do 561 # 1) summarize the far threshold over the entire dataset 562

Source Code for Module pylal.imr_utils