Package pylal :: Module ligolw_cwb_to_coinc
[hide private]
[frames] | no frames]

Source Code for Module pylal.ligolw_cwb_to_coinc

  1  #!/usr/bin/python 
  2  # 
  3  # Copyright (C) 2011 Chris Pankow 
  4  # 
  5  # This program is free software; you can redistribute it and/or modify it 
  6  # under the terms of the GNU General Public License as published by the 
  7  # Free Software Foundation; either version 3 of the License, or (at your 
  8  # option) any later version. 
  9  # 
 10  # This program is distributed in the hope that it will be useful, but 
 11  # WITHOUT ANY WARRANTY; without even the implied warranty of 
 12  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
 13  # Public License for more details. 
 14  # 
 15  # You should have received a copy of the GNU General Public License along 
 16  # with this program; if not, write to the Free Software Foundation, Inc., 
 17  # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA. 
 18   
 19  # TODO: 
 20  # 2. Find a way to include other tables which require the IDs to be linked 
 21  # (i. e. coinc_events and such). 
 22   
 23  """ 
 24  Convert a ROOT file produced by cWB into a ligolw document. 
 25  """ 
 26   
 27  import sys 
 28  import os 
 29  import platform 
 30  import re 
 31   
 32  # ligo_lw xml handling modules 
 33  import glue 
 34   
 35  from glue import segments 
 36  from glue.ligolw import ligolw 
 37  from glue.ligolw import lsctables 
 38  from glue.ligolw import table 
 39  from glue.ligolw import types 
 40  from glue.ligolw import utils 
 41  from glue.ligolw import ilwd 
 42  from glue.ligolw.utils import time_slide 
 43   
 44  from glue.lal import LIGOTimeGPS 
 45  from glue.lal import CacheEntry 
 46  from lalburst import timeslides as ligolw_tisi 
 47  from pylal import llwapp 
 48  from pylal import git_version 
 49   
 50  # TODO: Reimplement this when the cwb_table module is included 
 51  #try: 
 52    # Import cwb module 
 53    #import cwb_table 
 54  #except: 
 55    #sys.exit("ERROR: Unable to import modules from cwb_table.") 
 56   
 57  try: 
 58    from ROOT import TFile, TChain, TTree 
 59    from ROOT import gDirectory 
 60  except ImportError: 
 61    print >>sys.stderr, "WARNING: Unable to import modules from ROOT. ROOT file processing will not be possible." 
 62   
 63  #CONVERT_ROOT_VERSION = "0.0.8" 
 64  __author__ = "Chris Pankow <chris.pankow@ligo.org>" 
 65  __version__ = "git id %s" % git_version.id 
 66  __date__ = git_version.date 
 67   
 68  # 
 69  # ============================================================================= 
 70  # 
 71  #                                  Helper Functions 
 72  # 
 73  # ============================================================================= 
 74  # 
 75   
76 -def branch_array_to_list(branch, len):
77 """ 78 Turn a ROOT TBranch into a python list 79 """ 80 81 list = [] 82 for i in range(0, len): 83 list.append( branch[i] ) 84 85 return list
86
87 -def get_ifos_from_index(indx):
88 """ 89 Index map copied from detector definitions in WaveBurst. Constructs a string with all relevant detector strings from a set of indexes given. 90 """ 91 92 if not isinstance(indx, list): 93 indx = [indx] 94 95 ifolist = [] 96 for i in indx : 97 if i == 1 : ifolist.append("L1") 98 if i == 2 : ifolist.append("H1") 99 if i == 3 : ifolist.append("H2") 100 if i == 4 : ifolist.append("G1") 101 if i == 5 : ifolist.append("T1") 102 if i == 6 : ifolist.append("V1") 103 if i == 7 : ifolist.append("A1") 104 105 if len(ifolist) == 0 : return None 106 if len(ifolist) == 1 : return ifolist[0] 107 return ifolist
108 109 # 110 # ============================================================================= 111 # 112 # XML Table Functions 113 # 114 # ============================================================================= 115 # 116
117 -class CWB2Coinc(object):
118 """ 119 Class to convert a set of rootfiles to a ligolw_document. 120 """ 121
122 - def __init__(self, joblist=None, start=None, end=None, instruments=None, 123 #cwbtables=False, 124 waveoffset=0, verbose=False):
125 self.job_list = joblist 126 self.start = start 127 self.end = end 128 #self.cwbtable = cwbtables 129 self.verbose = verbose 130 self.instruments = lsctables.ifos_from_instrument_set( instruments.split(",") ) 131 self.waveoffset = waveoffset
132
133 - def create_tables(self, xmldoc, rootfiles):
134 """ 135 Sets up table structures and calls populating methods. 136 """ 137 138 if os.path.splitext(rootfiles[0])[1] == ".root": 139 sim_tree = TChain("waveburst") 140 else: # If the file is (for example) text, use a proxy class 141 sim_tree = CWBTextConverter() 142 143 for rootfile in rootfiles: 144 sim_tree.Add(rootfile) 145 146 # Define tables 147 sngl_burst_table = lsctables.New(lsctables.SnglBurstTable, 148 ["peak_time_ns", "start_time_ns", "stop_time_ns", 149 "process_id", "ifo", "peak_time", "start_time", "stop_time", 150 "duration", "time_lag", "peak_frequency", "search", 151 "flow", "fhigh", "bandwidth", "tfvolume", "hrss", "event_id", "snr"]) 152 xmldoc.childNodes[0].appendChild(sngl_burst_table) 153 sngl_burst_table.sync_next_id() 154 155 coinc_event_table = lsctables.New(lsctables.CoincTable, 156 [ "process_id", "coinc_event_id", "nevents", "instruments", "time_slide_id", 157 "coinc_def_id", "likelihood"]) 158 xmldoc.childNodes[0].appendChild(coinc_event_table) 159 coinc_event_table.sync_next_id() 160 161 # TODO: Reimplement this when the cwb_table module is included 162 #if self.cwbtable: 163 #cohwb_table = lsctables.New(cwb_table.CoherentWaveburstTable, 164 #["ellipticity", "correlated_energy", "eff_correlated_energy", 165 #"coherent_network_amp", "penalty", "network_correlation", 166 #"energy_disbalance", "ellip_energy_disbalance", "process_id", 167 #"coinc_event_id", "cwb_id"]) 168 #xmldoc.childNodes[0].appendChild(cohwb_table) 169 #cohwb_table.sync_next_id() 170 171 multi_burst_table = lsctables.New(lsctables.MultiBurstTable, 172 ["process_id", "peak_time", "peak_time_ns", "coinc_event_id", "snr", "ifos", 173 # NOTE: Added to the table definition 174 "false_alarm_rate", "ligo_axis_ra", "ligo_axis_dec"]) 175 xmldoc.childNodes[0].appendChild(multi_burst_table) 176 177 coinc_event_map_table = lsctables.New(lsctables.CoincMapTable) 178 xmldoc.childNodes[0].appendChild(coinc_event_map_table) 179 180 jobsegment = None 181 if self.start and self.end : 182 jobsegment = segments.segment(LIGOTimeGPS(self.start), LIGOTimeGPS(self.end)) 183 184 if self.verbose: 185 print "Creating Process Table...", 186 187 if self.job_list: 188 self.do_process_table(xmldoc, sim_tree) 189 else : 190 self.do_process_table_from_segment(xmldoc, sim_tree, jobsegment) 191 192 process_index = dict((int(row.process_id), row) for row in lsctables.ProcessTable.get_table(xmldoc)) 193 194 if self.verbose: 195 print " done." 196 197 if self.verbose: 198 print "Creating Summary Table...", 199 # If we are provided a joblist, use it to generate the list 200 if self.job_list : 201 self.do_summary_table_from_joblist(xmldoc, sim_tree) 202 elif self.job_list == None and self.start and self.end : 203 self.do_summary_table_from_segment(xmldoc, jobsegment, sim_tree) 204 else : 205 self.do_summary_table(xmldoc, sim_tree) 206 if self.verbose: 207 print " done." 208 209 # create coinc_definer row 210 row = self.get_coinc_def_row(sim_tree) 211 coinc_def_id = llwapp.get_coinc_def_id(xmldoc, row.search, row.search_coinc_type, description = row.description) 212 213 entries = sim_tree.GetEntries() 214 for i in range(0, entries): 215 sim_tree.GetEntry(i) 216 if self.start != None : 217 if float(self.start) > sim_tree.start[0] : continue 218 if self.end != None : 219 if float(self.end) < sim_tree.stop[0] : continue 220 221 222 offset_vector = dict((get_ifos_from_index(instrument_index), offset) for instrument_index, offset in zip(sim_tree.ifo, sim_tree.lag)) 223 224 coinc_event = coinc_event_table.RowType() 225 coinc_event.process_id = process_index[sim_tree.run].process_id 226 coinc_event.coinc_event_id = coinc_event_table.get_next_id() 227 coinc_event.coinc_def_id = coinc_def_id 228 coinc_event.nevents = sim_tree.ndim 229 coinc_event.instruments = lsctables.ifos_from_instrument_set( get_ifos_from_index( branch_array_to_list ( sim_tree.ifo, sim_tree.ndim ) ) ) 230 coinc_event.time_slide_id = time_slide.get_time_slide_id(xmldoc, offset_vector, process_index[sim_tree.run]) 231 coinc_event.likelihood = sim_tree.likelihood 232 coinc_event_table.append(coinc_event) 233 234 for d in range(0, sim_tree.ndim): 235 sngl_burst = self.get_sngl_burst_row(sngl_burst_table, sim_tree, d) 236 sngl_burst.process_id = coinc_event.process_id 237 sngl_burst.event_id = sngl_burst_table.get_next_id() 238 sngl_burst_table.append(sngl_burst) 239 240 coinc_event_map = coinc_event_map_table.RowType() 241 coinc_event_map.event_id = sngl_burst.event_id 242 coinc_event_map.table_name = sngl_burst.event_id.table_name 243 coinc_event_map.coinc_event_id = coinc_event.coinc_event_id 244 coinc_event_map_table.append(coinc_event_map) 245 246 # TODO: Reimplement when cwb_table module is included 247 #if self.cwbtable: 248 #cwb = get_cwb_row(cohwb_table, sim_tree) 249 #cwb.process_id = coinc_event.process_id 250 #cwb.coinc_event_id = coinc_event.coinc_event_id 251 #cwb.cwb_id = cohwb_table.get_next_id() 252 #cohwb_table.append(cwb) 253 254 multi_burst = self.get_multi_burst_row(multi_burst_table, sim_tree) 255 multi_burst.process_id = coinc_event.process_id 256 multi_burst.coinc_event_id = coinc_event.coinc_event_id 257 # NOTE: Until we have an embedded cwb table definition, this will be 258 # copied here so official tools have a ranking statistic to use 259 try: 260 multi_burst.snr = sim_tree.rho[1] 261 except TypeError: # difference in definition between ROOT and text 262 multi_burst.snr = sim_tree.rho 263 # NOTE: To be filled in later by farburst 264 multi_burst.false_alarm_rate = -1.0 265 # Reconstructed right ascension and declination 266 multi_burst.ligo_axis_ra = sim_tree.phi[2] 267 multi_burst.ligo_axis_dec = sim_tree.theta[2] 268 multi_burst.ifos = lsctables.ifos_from_instrument_set( get_ifos_from_index( branch_array_to_list ( sim_tree.ifo, sim_tree.ndim ) ) ) 269 270 multi_burst_table.append(multi_burst)
271
272 - def get_coinc_def_row(self, sim_tree):
273 """ 274 Fill in a coinc definer entry for waveburst 275 """ 276 277 return lsctables.CoincDef(search = u"waveburst", search_coinc_type = 0, description = u"fully coherent search (burst events)")
278
279 - def get_cwb_row(self, cohwb_table, sim_tree):
280 """ 281 Fill in a custom cwb event table 282 """ 283 284 row = cohwb_table.RowType() 285 # ellipticity of reconstructed event 286 row.ellipticity = sim_tree.norm 287 # sum of correlated energy 288 row.correlated_energy = sim_tree.ecor 289 # sum of correlated energy weighted by likelihood matrix 290 row.eff_correlated_energy = sim_tree.ECOR 291 # Coherent Network Amplitude 292 # The currently used definition uses effective correlated energy which 293 # is the second entry in the array (the first is normal corrrelated energy) 294 row.coherent_network_amp = sim_tree.rho[1] 295 # Penalty factor 296 # is an expression of the mismatch between signal energy in the h(t) 297 # and the energy in the reconstructed signal, if the recon. signal energy is 298 # greater than the data signal energy -- otherwise this is identically unity. 299 row.penalty = sim_tree.penalty 300 # Network Correlation Coefficient 301 # Correlated energy divided by sum of null and correlated energy 302 # Note again, the first entry is used, corresponding to ecor rather than 303 # eff_ecor in this case 304 row.network_correlation = sim_tree.netcc[0] 305 # Energy Disbalance 306 # Summed difference of the reconstructed and data signal energies 307 # See Penaly factor for description of above 308 row.energy_disbalance = sim_tree.neted[0]/sim_tree.ecor 309 # ELliptical Energy Disbalance 310 # Same as above, but between the signal and phase shifted signal 311 row.ellip_energy_disbalance = sim_tree.neted[1]/sim_tree.ecor 312 313 return row
314
315 - def get_sngl_burst_row(self, sngl_burst_table, sim_tree, d):
316 """ 317 Fill in a sngl_burst row for a cwb event 318 """ 319 row = sngl_burst_table.RowType() 320 row.search = u"waveburst" 321 # Interferometer name -> ifo 322 row.ifo = get_ifos_from_index(sim_tree.ifo[d] ) 323 # Timing 324 peak = LIGOTimeGPS(sim_tree.time[d]) 325 seg = segments.segment(LIGOTimeGPS(sim_tree.start[d]), LIGOTimeGPS(sim_tree.stop[d])) 326 # Central time in the detector -> cent_time 327 row.set_peak(peak) 328 # Start time in the detector -> end_time 329 row.set_start(seg[0]) 330 # Stop time in the detector -> star_time 331 row.set_stop(seg[1]) 332 # Event duration 333 row.duration = abs(seg) 334 # TODO: Make sure this is right = Time lag used to shift detector -> lag 335 row.time_lag = sim_tree.lag[d] 336 # Frequency 337 # Central frequency in the detector -> frequency 338 row.peak_frequency = sim_tree.frequency[d] 339 try: 340 # Low frequency of the event in the detector -> flow 341 row.flow = sim_tree.low[d] 342 # High frequency of the event in the detector -> fhigh 343 row.fhigh = sim_tree.high[d] 344 except TypeError: 345 row.flow = sim_tree.low 346 row.fhigh = sim_tree.high 347 # Bandwidth 348 row.bandwidth = sim_tree.bandwidth[d] 349 # Shape 350 # number of pizels on the TF plane -> tfsize 351 row.tfvolume = sim_tree.size[d] 352 # Energy 353 # energy / noise variance -> snr 354 row.snr = sim_tree.snr[d]**(1./2.) 355 # TODO: What to do with this? GW strain 356 #row.strain = sim_tree.strain[d]) 357 # h _ root square sum 358 row.hrss = sim_tree.hrss[d] 359 360 return row
361
362 - def get_multi_burst_row(self, multi_burst_table, sim_tree):
363 """ 364 Fill in a multi_burst row for a cwb event -- these should exactly correlate to the sngl_burst events within a single coherent event 365 """ 366 row = multi_burst_table.RowType() 367 # TODO : Script chokes on the line below: we need a better way to handle this anyway 368 #row.set_peak(sum(LIGOTimeGPS(t) for t in sim_tree.time[0:3]) / int(sim_tree.ndim)) 369 row.set_peak(LIGOTimeGPS(sim_tree.time[0])) 370 return row
371
372 - def do_process_table_from_segment(self, xmldoc, sim_tree, segment, jobid=-1):
373 """ 374 Create the process_table for the cWB job(s) from a job segment. 375 """ 376 377 try: 378 process_table = lsctables.ProcessTable.get_table(xmldoc) 379 except ValueError: 380 process_table = lsctables.New(lsctables.ProcessTable, 381 ["process_id", "ifos", "comment", "program", "start_time", "jobid", 382 "end_time", 383 # Unused, but required for column synchronicity 384 "version", "cvs_repository", "cvs_entry_time", "is_online", "node", 385 "username", "unix_procid", "domain"]) 386 xmldoc.childNodes[0].appendChild(process_table) 387 388 389 # get run IDs 390 if(jobid < 0): 391 N = sim_tree.GetEntries() 392 if N == 0: 393 exit("There is no information available to parse from the ROOT file, and no external information was provided.") 394 runs = set() 395 for i in xrange(N): 396 sim_tree.GetEntry(i) 397 run = sim_tree.run 398 runs.add(int(run)) 399 else: 400 run = jobid 401 seg = segment 402 403 # Imstruments involved in the search 404 ifos = lsctables.ifos_from_instrument_set( get_ifos_from_index( branch_array_to_list ( sim_tree.ifo, sim_tree.ndim ) ) ) 405 406 if( ifos == None or len(ifos) == 0 ): 407 if( self.instruments ): 408 ifos = self.instruments 409 else: # Not enough information to completely fill out the table 410 sys.exit("Found a job with no IFOs on, or not enough to determine IFOs. Try specifying instruments directly.") 411 412 for run in runs: 413 row = process_table.RowType() 414 row.process_id = type(process_table.next_id)(run) 415 row.ifos = ifos 416 row.comment = u"waveburst" 417 row.program = u"waveburst" 418 row.start_time = None 419 row.end_time = None 420 row.jobid = run 421 row.version = __version__ 422 row.cvs_repository = u"lscsoft" 423 # TODO: Need a handle on where this comes from 424 row.cvs_entry_time = 0 425 row.is_online = 0 426 row.username = os.environ['USER'] 427 row.node = platform.node() 428 row.unix_procid = os.getpid() 429 row.domain = u"pylal" 430 process_table.append(row)
431
432 - def do_process_table(self, xmldoc, sim_tree):
433 """ 434 Create the process_table for the cWB job(s) 435 """ 436 437 try: 438 process_table = lsctables.ProcessTable.get_table(xmldoc) 439 except ValueError: 440 process_table = lsctables.New(lsctables.ProcessTable, 441 ["process_id", "ifos", "comment", "program", "start_time", "jobid", 442 "end_time", 443 # Unused, but required for column synchronicity 444 "version", "cvs_repository", "cvs_entry_time", "is_online", "node", 445 "username", "unix_procid", "domain"]) 446 xmldoc.childNodes[0].appendChild(process_table) 447 448 runids = set() 449 runid = dict() 450 451 itr = 0 452 for line in file(self.job_list) : 453 if line[0] == '#' : 454 continue # skip comments 455 456 # Determine the file type 457 # WARNING: Mixing the two types could cause overwrite behavior 458 line = line.split() 459 if len(line) == 2 : # start and stop 460 seg = segments.segment( map( LIGOTimeGPS, map(float, line[0:2]) ) ) 461 runid[itr] = seg 462 itr += 1 463 elif len(line) == 3 : # index start and stop 464 seg = segments.segment( map( LIGOTimeGPS, map(float, line[1:3]) ) ) 465 runid[int(line[0])] = seg 466 elif len(line) == 4 : # index start and stop and length 467 seg = segments.segment( map( LIGOTimeGPS, map(float, line[1:3]) ) ) 468 runid[int(line[0])] = seg 469 else : # dunno! 470 sys.exit("Unable to understand job list segment format. Why do you confused me so?") 471 472 473 sim_tree.GetEntry(0) 474 for run, seg in runid.iteritems() : 475 if self.start != None : 476 if float(self.start) > seg[0] : continue 477 if self.end != None : 478 if float(self.end) < seg[1] : continue 479 480 row = process_table.RowType() 481 row.process_id = type(process_table.next_id)(run) 482 #print "Run id %d was given process_id %d" % (sim_tree.run, row.process_id) 483 runids.add(sim_tree.run) 484 485 # Imstruments involved in the search 486 ifos = lsctables.ifos_from_instrument_set( get_ifos_from_index( branch_array_to_list ( sim_tree.ifo, sim_tree.ndim ) ) ) 487 488 if( ifos == None or len(ifos) == 0 ): 489 if( self.instruments ): 490 ifos = self.instruments 491 else: # Not enough information to completely fill out the table 492 sys.exit("Found a job with no IFOs on, or not enough to determine IFOs. Try specifying instruments directly.") 493 494 row.ifos = ifos 495 row.comment = u"waveburst" 496 row.program = u"waveburst" 497 498 row.start_time = None 499 row.end_time = None 500 row.jobid = run 501 502 row.version = __version__ 503 row.cvs_repository = u"lscsoft" 504 # TODO: Need a handle on where this comes from 505 row.cvs_entry_time = 0 506 row.is_online = 0 507 row.username = os.environ['USER'] 508 row.node = platform.node() 509 row.unix_procid = os.getpid() 510 row.domain = "pylal" 511 512 process_table.append(row)
513
514 - def do_summary_table_from_segment(self, xmldoc, segment, sim_tree, jobid=-1):
515 """ 516 Create the search_summary table for a cWB from a segment specified from the command line. The function will try to determine the proper job intervals from the waveoffset, if specified. 517 """ 518 519 try: 520 search_summary = lsctables.SearchSummaryTable.get_table(xmldoc) 521 except ValueError: 522 search_summary = lsctables.New(lsctables.SearchSummaryTable, 523 ["process_id", "nevents", "ifos", "comment", "in_start_time", 524 "in_start_time_ns", "out_start_time", "out_start_time_ns", 525 "in_end_time", "in_end_time_ns", "out_end_time", "out_end_time_ns"]) 526 xmldoc.childNodes[0].appendChild(search_summary) 527 528 process_id_type = type(lsctables.ProcessTable.get_table(xmldoc).next_id) 529 530 sim_tree.GetEntry(0) 531 532 if(jobid < 0): 533 run = sim_tree.run 534 else: run = jobid 535 seg = segment 536 537 # Search Summary Table 538 # events found in the run -> nevents 539 row = search_summary.RowType() 540 row.process_id = process_id_type(run) 541 row.nevents = sim_tree.GetEntries() 542 543 ifos = lsctables.ifos_from_instrument_set( get_ifos_from_index( branch_array_to_list ( sim_tree.ifo, sim_tree.ndim ) ) ) 544 # Imstruments involved in the search 545 if( ifos == None or len(ifos) == 0 ): 546 if( self.instruments ): 547 ifos = self.instruments 548 else: # Not enough information to completely fill out the table 549 sys.exit("Found a job with no IFOs on, or not enough to determine IFOs. Try specifying instruments directly.") 550 551 row.ifos = ifos 552 row.comment = "waveburst" 553 554 # Begin and end time of the segment 555 waveoffset = self.waveoffset 556 if waveoffset == None: waveoffset = 0 557 558 # in -- with waveoffset 559 row.set_in(seg) 560 # out -- without waveoffset 561 waveoffset = LIGOTimeGPS(waveoffset) 562 row.set_out(segments.segment(seg[0]+waveoffset, seg[1]-waveoffset)) 563 search_summary.append(row)
564
565 - def do_summary_table_from_joblist(self, xmldoc, sim_tree):
566 """ 567 Create the search_summary table for the cWB job(s) and a provided cWB joblist. The function will try to determine the proper job intervals from the waveoffset, if specified. 568 """ 569 570 try: 571 search_summary = lsctables.SearchSummaryTable.get_table(xmldoc) 572 except ValueError: 573 search_summary = lsctables.New(lsctables.SearchSummaryTable, 574 ["process_id", "nevents", "ifos", "comment", "in_start_time", 575 "in_start_time_ns", "out_start_time", "out_start_time_ns", 576 "in_end_time", "in_end_time_ns", "out_end_time", "out_end_time_ns"]) 577 xmldoc.childNodes[0].appendChild(search_summary) 578 579 process_id_type = type(lsctables.ProcessTable.get_table(xmldoc).next_id) 580 581 runid = dict() 582 itr = 0 583 for line in file(self.job_list) : 584 if line[0] == '#' : 585 continue # skip comments 586 587 # Determine the file type 588 # WARNING: Mixing the two types could cause overwrite behavior 589 line = line.split() 590 591 if len(line) == 2 : # start and stop 592 seg = segments.segment( map( LIGOTimeGPS, map(float, line[0:2]) ) ) 593 runid[itr] = seg 594 itr += 1 595 elif len(line) == 3 : # index start and stop 596 seg = segments.segment( map( LIGOTimeGPS, map(float, line[1:3]) ) ) 597 runid[int(line[0])] = seg 598 elif len(line) == 4 : # index start and stop and length 599 seg = segments.segment( map( LIGOTimeGPS, map(float, line[1:3]) ) ) 600 runid[int(line[0])] = seg 601 else : # dunno! 602 sys.exit("Unable to understand job list segment format.") 603 604 for run, seg in runid.iteritems() : 605 if self.start != None : 606 if float(self.start) > seg[0] : continue 607 if self.end != None : 608 if float(self.end) < seg[1] : continue 609 610 row = search_summary.RowType() 611 row.process_id = process_id_type(run) 612 613 # Search Summary Table 614 # events found in the run -> nevents 615 row.nevents = 0 616 #entries = sim_tree.GetEntries() 617 618 # Imstruments involved in the search 619 sim_tree.GetEntry(0) 620 ifos = lsctables.ifos_from_instrument_set( get_ifos_from_index( branch_array_to_list ( sim_tree.ifo, sim_tree.ndim ) ) ) 621 # Imstruments involved in the search 622 if( ifos == None or len(ifos) == 0 ): 623 if( self.instruments ): 624 ifos = self.instruments 625 else: # Not enough information to completely fill out the table 626 sys.exit("Found a job with no IFOs on, or not enough to determine IFOs. Try specifying instruments directly.") 627 628 row.ifos = ifos 629 row.comment = "waveburst" 630 631 # Begin and end time of the segment 632 # TODO: This is a typical offset on either side of the job for artifacts 633 # It can, and probably will change in the future, and shouldn't be hardcoded 634 #waveoffset, livetime = 8, -1 635 waveoffset, livetime = self.waveoffset, -1 636 if waveoffset == None: waveoffset = 0 637 livetime = abs(seg) 638 639 if livetime < 0: 640 print >>sys.stderr, "WARNING: Run %d will have zero livetime because no events are recorded in this run, and therefore livetime cannot be calculated." % run 641 # in -- with waveoffset 642 row.set_in(segments.segment(LIGOTimeGPS(0), LIGOTimeGPS(1))) 643 # out -- without waveoffset 644 row.set_out(segments.segment(LIGOTimeGPS(0), LIGOTimeGPS(1))) 645 else: 646 # in -- with waveoffset 647 row.set_in(seg) 648 # out -- without waveoffset 649 row.set_out(segments.segment(seg[0]+waveoffset, seg[1]-waveoffset)) 650 651 search_summary.append(row)
652
653 - def do_summary_table(self, xmldoc, sim_tree):
654 """ 655 Create the search_summary table for the cWB job(s). This function exists as a backup in case no job list exists. It will try and reconstruct the job segments from the event data, but this list will be incomplete in the case where no events were found during a job. 656 """ 657 658 try: 659 search_summary = lsctables.SearchSummaryTable.get_table(xmldoc) 660 except ValueError: 661 search_summary = lsctables.New(lsctables.SearchSummaryTable, 662 ["process_id", "nevents", "ifos", "comment", "in_start_time", 663 "in_start_time_ns", "out_start_time", "out_start_time_ns", 664 "in_end_time", "in_end_time_ns", "out_end_time", "out_end_time_ns"]) 665 xmldoc.childNodes[0].appendChild(search_summary) 666 667 process_id_type = type(lsctables.ProcessTable.get_table(xmldoc).next_id) 668 669 runids = set() 670 entries = sim_tree.GetEntries() 671 for i in range(0, entries) : 672 sim_tree.GetEntry(i) 673 674 if self.start != None : 675 if float(self.start) > sim_tree.start[0] : continue 676 if self.end != None : 677 if float(self.end) < sim_tree.stop[0] : continue 678 679 # Id for the run processed by WaveBurst -> process ID 680 run = sim_tree.run 681 if run in runids : 682 continue 683 684 row = search_summary.RowType() 685 row.process_id = process_id_type(run) 686 runids.add(run) 687 688 # Search Summary Table 689 # events found in the run -> nevents 690 # TODO: Destroy ROOT, because it hates me 691 #row.nevents = sim_tree.nevent 692 row.nevents = 0 693 694 # Imstruments involved in the search 695 ifos = lsctables.ifos_from_instrument_set( get_ifos_from_index( branch_array_to_list ( sim_tree.ifo, sim_tree.ndim ) ) ) 696 # Imstruments involved in the search 697 if( ifos == None or len(ifos) == 0 ): 698 if( self.instruments ): 699 ifos = self.instruments 700 else: # Not enough information to completely fill out the table 701 sys.exit("Found a job with no IFOs on, or not enough to determine IFOs. Try specifying instruments directly.") 702 703 row.ifos = ifos 704 row.comment = "waveburst" 705 706 # Begin and end time of the segment 707 if self.waveoffset != None : 708 waveoffset = self.waveoffset 709 else: 710 waveoffset = 0 711 712 livetime = -1 713 # WARNING: This will fail miserably if there are no events in analyzed 714 # see do_search_summary_from_joblist for a better solution 715 livetime = sim_tree.left[0] + sim_tree.duration[0] + sim_tree.right[0] 716 717 if livetime < 0: 718 print >>sys.stderr, "WARNING: Run %d will have zero livetime because no events are recorded in this run, and therefore livetime cannot be calculated." % run 719 # in -- with waveoffset 720 row.set_in(segments.segment(LIGOTimeGPS(0), LIGOTimeGPS(1))) 721 # out -- without waveoffset 722 row.set_out(segments.segment(LIGOTimeGPS(0), LIGOTimeGPS(1))) 723 else: 724 725 seg_start_with_offset = sim_tree.gps - waveoffset 726 seg_start_without_offset = sim_tree.gps 727 seg_end_with_offset = sim_tree.gps + waveoffset + livetime 728 seg_end_without_offset = sim_tree.gps + livetime 729 # in -- with waveoffset 730 row.set_in(segments.segment(LIGOTimeGPS(seg_start_with_offset), LIGOTimeGPS(seg_end_with_offset))) 731 # out -- without waveoffset 732 row.set_out(segments.segment(LIGOTimeGPS(seg_start_without_offset), LIGOTimeGPS(seg_end_without_offset))) 733 734 search_summary.append(row)
735
736 -class CWBTextConverter(object):
737 """ 738 Class to transparently emit data as if the input were a ROOT file. The input is an unstructured text file made by cWB at runtime. 739 """ 740 741 CWB_SEARCH_TYPES = { 'r': "unmodelled", 'i': "elliptical", 's': "linear", 'g': "circular", 'R': "unmodelled", 'I': "elliptical", 'S': "linear", 'G': "circular" } 742
743 - def __init__(self):
744 self.data = [] 745 self.this_data = None 746 self.version = None 747 self.online_version = None 748 self.search_type = None 749 self._skymap_reg = re.compile( "^map_lenght" ) # NOT A TYPO, THIS IS ACTUALLY IN THE FILE. 750 self._far_info_reg = re.compile( "^##" )
751
752 - def __getattr__(self, name):
753 """ 754 This overrides the attribute getter for this class so that when the proxy class is required to produce information based on an attribute (say 'rho') it responds just like a TBranch would. If no event is 'selected' via GetEntry, then None is returned, regardless of the attribute, because there's no way to know whether we're asking for a TBranch leaf or something else. 755 """ 756 if self.this_data is None: 757 return None 758 try: 759 return self.this_data[name] 760 except KeyError: 761 return getattr(self, name)
762
763 - def GetEntry(self, i):
764 """ 765 Emulate a TTree GetEntry() call by switching internal pointer to a row in the table. 766 """ 767 self.this_data = self.data[i]
768
769 - def GetEntries(self):
770 """ 771 Emulate a TTree GetEntries() call by returning the length of the internal data table. 772 """ 773 return len(self.data)
774
775 - def Add(self, filen):
776 """ 777 Add a file which will be parsed by the class and added to the internal table. 778 """ 779 self.parse_transcript( open(filen) )
780
781 - def get_far_segments( self ):
782 """ 783 Gets the segments for which various FARs have been defined. 784 """ 785 if self.this_data["far_data"] is None: 786 raise ValueError( "No FAR info for this event." ) 787 788 return self.this_data["far_data"].keys()
789 790
791 - def get_far_information( self, segment=None ):
792 """ 793 Returns the false alarm rate for a given segment, usually one of the segments retrieved from 'get_far_segments'. If segment is None (default), then the FAR from the biggest segment is returned. 794 """ 795 796 if self.this_data["far_data"] is None: 797 raise ValueError( "No FAR info for this event." ) 798 799 # return FAR from largest segment of data 800 if segment is None: 801 segs = sorted( self.get_far_segments(), key=lambda s: abs(s) ) 802 return self.this_data["far_data"][segs[-1]] 803 804 return self.this_data["far_data"][segment]
805 806
807 - def get_skymap_coords( self ):
808 """ 809 Get a list of the skymap coordinates (R.A. and dec) which cWB has probability values for. 810 """ 811 if self.this_data["skymap"] is None: 812 raise ValueError( "No skymap set for this event." ) 813 814 return self.this_data["skymap"].keys()
815
816 - def get_skymap_prob( self, ra, dec ):
817 """ 818 Return the probability for a given R.A. and dec from the skymap data. If there is no entry for this, then 0 is returned. This is the implied value for no R.A. and dec entry in the skymap data. 819 """ 820 if self.this_data["skymap"] is None: 821 raise ValueError( "No skymap set for this event." ) 822 823 try: 824 return self.this_data["skymap"][(ra, dec)] 825 except KeyError: 826 return 0
827
828 - def parse_with_type(self, val):
829 """ 830 Attempt to coerce data to the proper type. Note that this depends on python's finickiness with int / float parsing. It won't do string -> float -> int truncation with one call to float(). So, this check whether or not int() works, then tries float(). Then it gives up and gives the value (presumed a string) back. 831 """ 832 try: 833 return int(val) 834 except ValueError: 835 pass 836 try: 837 return float(val) 838 except ValueError: 839 pass 840 return val
841
842 - def parse_transcript(self, text):
843 """ 844 Parses a textual event transcript and adds the data to the internal table. 845 """ 846 row, mode = {}, "var" 847 row["skymap"] = {} 848 row["far_data"] = {} 849 850 for line in text.readlines(): 851 852 # ALl the unstructured stuff is at the end... 853 if "wat version" in line: 854 self.version = line.split("=")[1].strip() 855 continue 856 elif "online version" in line: 857 self.online_version = line.split("=")[1].strip() 858 continue 859 elif "search" in line: 860 self.search_type = self.CWB_SEARCH_TYPES[line.split("=")[1].strip()] 861 continue 862 elif re.match( self._far_info_reg, line ) is not None: 863 mode = "far" 864 continue 865 elif re.match( self._skymap_reg, line ) is not None: 866 mode = "skymap" 867 continue 868 elif len(line.strip()) == 0: 869 # a blank line is assumed to be a new event entry 870 row["gps"] = row["segment"][0] 871 self.data.append( row ) 872 row = {} 873 row["skymap"] = {} 874 row["far_data"] = {} 875 mode = "var" 876 continue 877 878 if line[0] == "#": continue 879 880 # Modal parsing logic 881 if mode == "var": 882 883 varn, val = line.split(":") 884 val = [ self.parse_with_type(v) for v in val.split() ] 885 if len(val) == 1: val = val[0] 886 887 row[varn] = val 888 elif mode == "skymap": 889 sid, theta, dec, step, phi, ra, step2, prob, cum = line.split() 890 row["skymap"][(float(ra), float(dec))] = float(prob) 891 elif mode == "far": 892 try: 893 rank, rate, segstart, segend, segdur = line.split() 894 except ValueError: 895 # There is now two lines at the end of the file which 896 # contain extraneous paths and URLs 897 continue 898 far_seg = segments.segment( float(segstart), float(segend) ) 899 row["far_data"][far_seg] = float(rate) 900 901 # Not currently in the file, but an alias is needed for compatibility 902 row["gps"] = row["segment"][0] 903 self.data.append( row ) 904 return row
905