pylal.skylocutils

1 #!/usr/bin/python 2 3 import git_version 4 5 __author__ = "Larry Price <larry.price@ligo.org> and Patrick Brady <patrick.brady@ligo.org>" 6 __version__ = "git id %s" % git_version.id 7 __date__ = git_version.date 8 9 import sys 10 import os 11 import operator 12 import gzip 13 from math import sqrt, sin, cos, modf 14 import numpy 15 from numpy import pi, linspace, interp, sum as npsum, exp, asarray 16 from bisect import bisect 17 18 from pylal import date 19 from pylal import CoincInspiralUtils, SnglInspiralUtils, SimInspiralUtils 20 from pylal import inject 21 import lal 22 from pylal.sphericalutils import angle_between_points 23 24 import glue.iterutils 25 from glue.ligolw import utils, table as tab, lsctables, ligolw 26 lsctables.use_in(ligolw.LIGOLWContentHandler) 27 28 29 30 31 ############################################################################## 32 # 33 # global variables 34 # 35 ############################################################################## 36 37 #set the detector locations 38 detector_locations = {} 39 detector_locations["L1"] = inject.cached_detector["LLO_4k"].location 40 detector_locations["H1"] = inject.cached_detector["LHO_4k"].location 41 detector_locations["V1"] = inject.cached_detector["VIRGO"].location 42 43 #set the detector responses 44 detector_responses = {} 45 detector_responses["L1"] = inject.cached_detector["LLO_4k"].response 46 detector_responses["H1"] = inject.cached_detector["LHO_4k"].response 47 detector_responses["V1"] = inject.cached_detector["VIRGO"].response 48 49 ############################################################################## 50 # 51 # function definitions 52 # 53 ############################################################################## 54

55 -def get_delta_t_rss(pt,coinc,reference_frequency=None):

56 """ 57 returns the rss timing error for a particular location in 58 the sky (longitude,latitude) 59 """ 60 latitude,longitude = pt 61 earth_center = (0.0,0.0,0.0) 62 tref = {} 63 tgeo={} 64 for ifo in coinc.ifo_list: 65 66 if reference_frequency: 67 tFromRefFreq = get_signal_duration(ifo,coinc,reference_frequency) 68 tref[ifo] = lal.LIGOTimeGPS(int(tFromRefFreq), 1.e9*(tFromRefFreq-int(tFromRefFreq))) 69 else: 70 tref[ifo] = 0.0 71 72 #compute the geocentric time from each trigger 73 tgeo[ifo] = coinc.gps[ifo] - tref[ifo] - \ 74 lal.LIGOTimeGPS(0,1.0e9*date.XLALArrivalTimeDiff(detector_locations[ifo],\ 75 earth_center,longitude,latitude,coinc.gps[ifo])) 76 77 #compute differences in these geocentric times 78 time={} 79 delta_t_rms = 0.0 80 for ifos in coinc.ifo_coincs: 81 time[ifos[0]+ifos[1]] = float(tgeo[ifos[0]] - tgeo[ifos[1]]) 82 delta_t_rms += time[ifos[0]+ifos[1]] * time[ifos[0]+ifos[1]] 83 84 return sqrt(delta_t_rms)

85

86 -def get_signal_duration(ifo,coinc,frequency):

87 """ 88 determine the amount by which the coalescence time must be translated to get the reference time 89 """ 90 M = coinc.mass1[ifo]+coinc.mass2[ifo] 91 mu = coinc.mass1[ifo]*coinc.mass2[ifo]/M 92 eta = mu/M 93 chirpM = pow(mu*mu*mu*M*M,1./5.) 94 M = M*4.92549095e-6 95 mu = mu*4.92549095e-6 96 chirpM = chirpM*4.92549095e-6 97 tau0 = 5./256. * pow(chirpM,-5./3.) * pow(pi*frequency,-8./3.) 98 tau1 = 5./(192.*mu*pi*pi*frequency*frequency) * (743./336. + 11./4.*eta) 99 tau1p5 = 1./(8.*mu) * pow(M/(pi*pi*pow(frequency,5.)),1./3.) 100 tau2 = 5./(128.*mu) * pow(M/(pi*pi*frequency*frequency),2./3.)\ 101 *(3058673./1016064. + 5429./1008.*eta + 617./144.*eta*eta) 102 duration = tau0 + tau1 - tau1p5 + tau2 103 104 return duration

105

106 -def get_delta_D_rss(pt,coinc):

107 """ 108 compute the rms difference in the ratio of the difference of the squares of Deff to 109 the sum of the squares of Deff between the measured values and a "marginalized" effective 110 distance this is just the squared Deff integrated over inclination and polarization which 111 is proportional to (F+^2 + Fx^2)^(-1) 112 """ 113 latitude,longitude = pt 114 gmst = {} 115 D_marg_sq = {} 116 F_plus = {} 117 F_cross = {} 118 for ifo in coinc.ifo_list: 119 gmst[ifo] = date.XLALGreenwichMeanSiderealTime(coinc.gps[ifo]) 120 F_plus[ifo], F_cross[ifo] = lal.ComputeDetAMResponse(detector_responses[ifo],\ 121 longitude,latitude,0,gmst[ifo]) 122 D_marg_sq[ifo] = 1/(F_plus[ifo]*F_plus[ifo]+F_cross[ifo]*F_cross[ifo]) 123 124 delta_D = {} 125 effD_diff = 0.0 126 effD_sum = 0.0 127 Dmarg_diff = 0.0 128 Dmarg_sum = 0.0 129 delta_D_rss = 0.0 130 for ifos in coinc.ifo_coincs: 131 effD_diff = coinc.eff_distances[ifos[0]] * coinc.eff_distances[ifos[0]]\ 132 - coinc.eff_distances[ifos[1]] * coinc.eff_distances[ifos[1]] 133 effD_sum = coinc.eff_distances[ifos[0]] * coinc.eff_distances[ifos[0]]\ 134 + coinc.eff_distances[ifos[1]] * coinc.eff_distances[ifos[1]] 135 Dmarg_diff = D_marg_sq[ifos[0]] - D_marg_sq[ifos[1]] 136 Dmarg_sum = D_marg_sq[ifos[0]] + D_marg_sq[ifos[1]] 137 delta_D[ifos[0]+ifos[1]] = (effD_diff/effD_sum) - (Dmarg_diff/Dmarg_sum) 138 delta_D_rss += delta_D[ifos[0]+ifos[1]]*delta_D[ifos[0]+ifos[1]] 139 140 return sqrt(delta_D_rss)

141

142 -def gridsky(resolution,shifted=False):

143 """ 144 grid the sky up into roughly square regions 145 resolution is the length of a side 146 the points get placed at the center of the squares and to 147 first order each square has an area of resolution^2 148 if shifted is true, the grids are reported with latitudes 149 in (0,pi). otherwise (default) they lie in (-pi/2,pi/2) 150 """ 151 points = [] 152 latitude = 0.0 153 longitude = 0.0 154 ds = pi*resolution/180.0 155 if shifted: 156 dlat = 0.0 157 else: 158 dlat = 0.5*pi 159 while latitude <= pi: 160 latitude += ds 161 longitude = 0.0 162 points.append((latitude-dlat, longitude)) 163 while longitude <= 2.0*pi: 164 longitude += ds / abs(sin(latitude)) 165 points.append((latitude-dlat, longitude)) 166 #add a point at the south pole 167 points.append((0.0-dlat,0.0)) 168 #there's some slop so get rid of it and only focus on points on the sphere 169 sphpts = [] 170 if shifted: 171 latmin = 0.0 172 latmax = pi 173 else: 174 latmin = -pi/2 175 latmax = pi/2 176 for pt in points: 177 if pt[0] > latmax or pt[0] < latmin or pt[1] > 2*pi or pt[1] < 0.0: 178 pass 179 else: 180 sphpts.append(pt) 181 return sphpts

182

183 -def map_grids(coarsegrid,finegrid,coarseres=4.0):

184 """ 185 takes the two grids (lists of lat/lon tuples) and returns a dictionary 186 where the points in the coarse grid are the keys and lists of tuples of 187 points in the fine grid are the values 188 """ 189 fgtemp = finegrid[:] 190 coarsedict = {} 191 192 ds = coarseres*pi/180.0 193 epsilon = ds/10.0 194 for cpt in coarsegrid: 195 flist = [] 196 for fpt in fgtemp: 197 if (cpt[0]-fpt[0])*(cpt[0]-fpt[0]) - ds*ds/4.0 <= epsilon and \ 198 (cpt[1]-fpt[1])*(cpt[1]-fpt[1])*sin(cpt[0])*sin(cpt[0]) - ds*ds/4.0 <= epsilon: 199 flist.append(fpt) 200 coarsedict[cpt] = flist 201 for rpt in flist: 202 fgtemp.remove(rpt) 203 first_column = [pt for pt in coarsegrid if pt[1] == 0.0] 204 for cpt in first_column: 205 flist = [] 206 for fpt in fgtemp: 207 if (cpt[0]-fpt[0])*(cpt[0]-fpt[0]) - ds*ds/4.0 <= epsilon and \ 208 (2*pi-fpt[1])*(2*pi-fpt[1])*sin(cpt[0])*sin(cpt[0]) - ds*ds/4.0 <= epsilon: 209 flist.append(fpt) 210 coarsedict[cpt] = flist 211 for rpt in flist: 212 fgtemp.remove(rpt) 213 214 return coarsedict, fgtemp

215

216 -def gaussian_kde(data,x,w):

217 """ 218 kernel density estimate of the pdf represented by data 219 at point x with bandwidth w 220 """ 221 N = float(len(data)) 222 223 return npsum([_gauss_kern(x,xn,w) for xn in data])/(N*w*sqrt(2.*pi))

224

225 -def _gauss_kern(x,xn,w):

226 """ 227 gaussian kernel for kernel density estimator 228 """ 229 a = x-xn 230 return exp(-a*a/(2.*w*w))

231

232 -def percentile(p,dat):

233 """ 234 compute p%-tile of data in dat 235 """ 236 N = len(dat) 237 n = float(p)/100*(N-1) + 1 238 values = dat[:] 239 values.sort() 240 241 if n == 1: 242 return values[0] 243 elif n == N: 244 return values[N-1] 245 else: 246 dpart, ipart = modf(n) 247 return values[int(ipart)-1] + dpart*(values[int(ipart)] - values[int(ipart)-1])

248 249

250 -def iqr(data):

251 """ 252 computes the interquartile range of data 253 useful for determing widths of bins in histograms or bandwidths in kdes 254 """ 255 return (percentile(75.,data) - percentile(25.,data))

256 257 #borrowed this from nick's ligolw_cbc_jitter_skyloc 258 #if/when it makes its way to sphericalutils use that instead

259 -def lonlat2polaz(lon, lat):

260 """ 261 Convert (longitude, latitude) in radians to (polar, azimuthal) angles 262 in radians. 263 """ 264 return lal.PI_2 - lat, lon

265

266 -def sbin(bins,pt,res=0.4):

267 """ 268 bin points on the sky 269 returns the index of the point in bin that is closest to pt 270 """ 271 #there's probably a better way to do this, but it works 272 #the assumption is that bins is generated from gridsky 273 ds = pi*res/180.0 274 #first pare down possibilities 275 xindex = bisect(bins,pt) 276 xminus = xindex - 1 277 newbins = [] 278 while 1: 279 #the 3 in the denominator should be a 4, but this allows for a little slop 280 #that helps you get closer in ra in the next step 281 if xindex < len(bins)-1 \ 282 and (bins[xindex][0]-pt[0])*(bins[xindex][0]-pt[0]) <= ds*ds/4.0: 283 newbins.append((bins[xindex][1],bins[xindex][0])) 284 xindex += 1 285 else: 286 break 287 while 1: 288 if (bins[xminus][0]-pt[0])*(bins[xminus][0]-pt[0]) <= ds*ds/4.0: 289 newbins.append((bins[xminus][1],bins[xminus][0])) 290 xminus -= 1 291 else: 292 break 293 if len(newbins) == 1: 294 return bins.index((newbins[0][1],newbins[0][0])) 295 #now break out the full spherical distance formula 296 newbins.sort() 297 rpt = (pt[1],pt[0]) 298 yindex = bisect(newbins,rpt) 299 finalbins = {} 300 #if it's the last index then work backwards from the bottom 301 if yindex > len(newbins)-1: 302 print yindex 303 print len(newbins)-1 304 mindist = angle_between_points(asarray(rpt),asarray(newbins[len(newbins)-1])) 305 finalbins[newbins[len(newbins)-1]] = mindist 306 i = 2 307 while 1: 308 angdist = angle_between_points(asarray(rpt),asarray(newbins[len(newbins)-i])) 309 if angdist <= mindist: 310 finalbins[newbins[len(newbins)-i]] = angdist 311 i += 1 312 else: 313 break 314 #make sure to cover the top, too 315 i = 0 316 while 1: 317 angdist = angle_between_points(asarray(rpt),asarray(newbins[i])) 318 if angdist <= mindist: 319 finalbins[newbins[i]] = angdist 320 i += 1 321 else: 322 break 323 else: 324 mindist = angle_between_points(asarray(rpt),asarray(newbins[yindex])) 325 finalbins[newbins[yindex]] = mindist 326 i = 1 327 while yindex + i < len(newbins) -1: 328 angdist = angle_between_points(asarray(rpt),asarray(newbins[yindex+i])) 329 if angdist <= mindist: 330 finalbins[newbins[yindex+i]] = angdist 331 i += 1 332 else: 333 break 334 i = 1 335 while yindex - i >= 0: 336 angdist = angle_between_points(asarray(rpt),asarray(newbins[yindex-i])) 337 if angdist <= mindist: 338 finalbins[newbins[yindex-i]] = angdist 339 i += 1 340 else: 341 break 342 mindist = min(finalbins.values()) 343 for key in finalbins.keys(): 344 if finalbins[key] == mindist: 345 sky_bin = key 346 347 return bins.index((sky_bin[1],sky_bin[0]))

348 349 350 351 ############################################################################## 352 # 353 # class definitions 354 # 355 ############################################################################## 356

357 -class Ranking(object):

358 """ 359 class for storing pdfs 360 """

361 - def __init__(self,xvals,yvals):

362 """ 363 storing pdfs 364 """ 365 self.x = xvals 366 self.y = yvals

367

368 - def get_rank(self,value):

369 """ 370 return the probability of value as obtained via linear interpolation 371 """ 372 return interp(value,self.x,self.y)

373 374

375 -class SkyPoints(list):

376 """ 377 useful class for doing sky localization. 378 assumes that it's a list of (latitude,longitude,L) tuples 379 and contains: 380 * a method for sorting those lists 381 * a method for writing itself to disk 382 """

383 - def nsort(self,n,rev=True):

384 """ 385 in place sort of (latitude,longitude,dt,dD...) tuples 386 according to the values in the nth column 387 """ 388 super(SkyPoints,self).sort(key=operator.itemgetter(n),reverse=rev)

389

390 - def normalize(self,n):

391 """ 392 in place normalization of the data in the n-th column 393 by the factor in normfac 394 """ 395 normfac = sum([pt[n] for pt in self]) 396 if normfac > 0.0: 397 for i in range(len(self)): 398 self[i][n] /= normfac 399 return normfac 400 else: 401 return -1

402

403 - def write(self,fname,post_dat,comment=None,debug=False,gz=False):

404 """ 405 write the grid to a text file 406 """ 407 self.nsort(1) 408 post_grid = '# normfac = ' + str(post_dat['normfac']) + '\n' 409 post_grid += 'snr = ' + str(post_dat['snr']) + '\n' 410 post_grid += 'FAR = ' + str(post_dat['FAR']) + '\n' 411 post_grid += 'gps = ' + str(post_dat['gps']) + '\n' 412 post_grid += '# ra' + '\t' + 'dec' + '\t' + 'probability (posterior)' + '\n' 413 post_grid_l = post_grid 414 for pt in self: 415 post_grid += str(pt[0][1]) + '\t' + str(pt[0][0]) + '\t' + str(pt[1]) + '\n' 416 for pt in self[:1000]: 417 post_grid_l += str(pt[0][1]) + '\t' + str(pt[0][0]) + '\t' + str(pt[1]) + '\n' 418 if comment: 419 post_grid += '# ' + comment + '\n' 420 post_grid_l += '# ' + comment + '\n' 421 if fname['galaxy']: 422 gal_grid = '# normfac = ' + str(post_dat['gnormfac']) + '\n' 423 gal_grid += 'snr = ' + str(post_dat['snr']) + '\n' 424 gal_grid += 'FAR = ' + str(post_dat['FAR']) + '\n' 425 gal_grid += 'gps = ' + str(post_dat['gps']) + '\n' 426 gal_grid += '# ra' + '\t' + 'dec' + '\t' + 'probability (posterior)' + '\n' 427 gal_grid_l = gal_grid 428 self.nsort(4) 429 for pt in self: 430 gal_grid += str(pt[0][1]) + '\t' + str(pt[0][0]) + '\t' + str(pt[4]) + '\n' 431 for pt in self[:1000]: 432 gal_grid_l += str(pt[0][1]) + '\t' + str(pt[0][0]) + '\t' + str(pt[4]) + '\n' 433 if gz: 434 fpost = gzip.open(fname['posterior']+'.gz', 'w') 435 fpost_l = gzip.open(fname['posterior'].replace('.txt','_lumin.txt')+'.gz', 'w') 436 if fname['galaxy']: 437 fgal = gzip.open(fname['galaxy']+'.gz', 'w') 438 fgal_l = gzip.open(fname['galaxy'].replace('.txt','_lumin.txt')+'.gz', 'w') 439 else: 440 fpost = open(fname['posterior'], 'w') 441 fpost_l = open(fname['posterior'].replace('.txt','_lumin.txt'), 'w') 442 if fname['galaxy']: 443 fgal = open(fname['galaxy'], 'w') 444 fgal_l = open(fname['galaxy'].replace('.txt','_lumin.txt'), 'w') 445 446 fpost.write(post_grid) 447 fpost_l.write(post_grid_l) 448 fpost.close() 449 fpost_l.close() 450 if fname['galaxy']: 451 fgal.write(gal_grid) 452 fgal_l.write(gal_grid_l) 453 fgal.close() 454 fgal_l.close()

455

456 -class CoincData(object):

457 """ 458 simple container for the information needed to run the sky localization code 459 """

460 - def __init__(self):

461 """ 462 here are all the things we need 463 """ 464 #start with data needed for every coinc 465 self.ifo_list = [] 466 self.ifo_coincs = [] 467 468 self.snr = {} 469 self.gps = {} 470 self.eff_distances = {} 471 self.mass1 = {} 472 self.mass2 = {} 473 474 self.time = None 475 self.FAR = 99 476 477 #this stuff is only needed for injections 478 self.is_injection = False 479 self.latitude_inj = None 480 self.longitude_inj = None 481 self.mass1_inj = None 482 self.mass2_inj = None 483 self.distance_inj = None 484 self.eff_distances_inj = {}

485 486

487 - def set_ifos(self,ifolist):

488 """ 489 set the ifo_list ans ifo_coincs from the list of ifos involved 490 in the coincidence 491 """ 492 self.ifo_list = ifolist 493 self.ifo_coincs.extend(list(glue.iterutils.choices(self.ifo_list,2)))

494

495 - def set_snr(self,snrdict):

496 self.snr = snrdict

497

498 - def set_gps(self,gpsdict):

499 self.gps = gpsdict 500 self.time = min(t for t in self.gps.values())

501

502 - def set_effDs(self,effDdict):

503 self.eff_distances = effDdict

504

505 - def set_masses(self,m1,m2):

506 self.mass1 = m1 507 self.mass2 = m2

508

509 - def set_inj_params(self,lat,lon,m1,m2,dist,effDs):

510 """ 511 set all of the injection parameters at once 512 """ 513 self.latitude_inj = lat 514 self.longitude_inj = lon 515 self.mass1_inj = m1 516 self.mass2_inj = m2 517 self.distance_inj = dist 518 self.eff_distances_inj = effDs

519

520 - def set_FAR(self,FAR_per_day):

521 self.FAR=FAR_per_day

522 523

524 -class Coincidences(list):

525 """ 526 takes input in either the form of coire files or coinc tables (xml format) 527 and produces a list of CoincData objects 528 """ 529

530 - def __init__(self,files,filetype='coinctable'):

531 """ 532 files is assumend to be a list of filenames 533 """ 534 if filetype == 'coinctable': 535 self.get_coincs_from_coinctable(files) 536 elif filetype == 'coire': 537 self.get_coincs_from_coire(files) 538 else: 539 print >>sys.stdout, 'Unknown input file type.' 540 sys.exit(0)

541

542 - def get_coincs_from_coinctable(self,files):

543 """ 544 read data from coinc tables (xml format) 545 546 FIXME: currently assumes one coinc per file!!! 547 """ 548 for file in files: 549 coinc = CoincData() 550 xmldoc = utils.load_filename(file, contenthandler=ligolw.LIGOLWContentHandler) 551 sngltab = lsctables.SnglInspiralTable.get_table(xmldoc) 552 coinc.set_snr(dict((row.ifo, row.snr) for row in sngltab)) 553 coinc.set_gps(dict((row.ifo, lal.LIGOTimeGPS(row.get_end())) for row in sngltab)) 554 #FIXME: this is put in place to deal with eff_distance = 0 555 # needs to be fixed upstream in the pipeline 556 effDs = list((row.ifo,row.eff_distance) for row in sngltab) 557 for eD in effDs: 558 if eD[1] == 0.: 559 effDs.append((eD[0],1.)) 560 effDs.remove(eD) 561 coinc.set_effDs(dict(effDs)) 562 # coinc.set_effDs(dict((row.ifo,row.eff_distance) for row in sngltab)) 563 coinc.set_masses(dict((row.ifo, row.mass1) for row in sngltab), \ 564 dict((row.ifo, row.mass2) for row in sngltab)) 565 ctab = lsctables.CoincInspiralTable.get_table(xmldoc) 566 #FIXME: ignoring H2 for now, but should be dealt in a better way 567 allifos = list(ctab[0].get_ifos()) 568 try: 569 allifos.remove('H2') 570 except ValueError: 571 pass 572 coinc.set_ifos(allifos) 573 if ctab[0].false_alarm_rate is not None: 574 coinc.set_FAR(ctab[0].false_alarm_rate) 575 576 try: 577 simtab = lsctables.SimInspiralTable.get_table(xmldoc) 578 row = siminsptab[0] 579 effDs_inj = {} 580 for ifo in coinc.ifo_list: 581 if ifo == 'H1': 582 effDs_inj[ifo] = row.eff_dist_h 583 elif ifo == 'L1': 584 effDs_inj[ifo] = row.eff_dist_l 585 elif ifo == 'V1': 586 effDs_inj[ifo] = row.eff_dist_v 587 dist_inj = row.distance 588 coinc.set_inj_params(row.latitude,row.longitude,row.mass1,row.mass2, \ 589 dist_inj,effDs_inj) 590 coinc.is_injection = True 591 #FIXME: name the exception! 592 except: 593 pass 594 595 self.append(coinc)

596

597 - def get_coincs_from_coire(self,files,stat='snr'):

598 """ 599 uses CoincInspiralUtils to get data from old-style (coire'd) coincs 600 """ 601 coincTrigs = CoincInspiralUtils.coincInspiralTable() 602 inspTrigs = SnglInspiralUtils.ReadSnglInspiralFromFiles(files, \ 603 mangle_event_id = True,verbose=None) 604 statistic = CoincInspiralUtils.coincStatistic(stat,None,None) 605 coincTrigs = CoincInspiralUtils.coincInspiralTable(inspTrigs,statistic) 606 607 try: 608 inspInj = SimInspiralUtils.ReadSimInspiralFromFiles(files) 609 coincTrigs.add_sim_inspirals(inspInj) 610 #FIXME: name the exception! 611 except: 612 pass 613 614 #now extract the relevant information into CoincData objects 615 for ctrig in coincTrigs: 616 coinc = CoincData() 617 coinc.set_ifos(ctrig.get_ifos()[1]) 618 coinc.set_gps(dict((trig.ifo,lal.LIGOTimeGPS(trig.get_end())) for trig in ctrig)) 619 coinc.set_snr(dict((trig.ifo,getattr(ctrig,trig.ifo).snr) for trig in ctrig)) 620 coinc.set_effDs(dict((trig.ifo,getattr(ctrig,trig.ifo).eff_distance) for trig in ctrig)) 621 coinc.set_masses(dict((trig.ifo,getattr(ctrig,trig.ifo).mass1) for trig in ctrig), \ 622 dict((trig.ifo,getattr(ctrig,trig.ifo).mass2) for trig in ctrig)) 623 624 try: 625 effDs_inj = {} 626 for ifo in coinc.ifo_list: 627 if ifo == 'H1': 628 effDs_inj[ifo] = getattr(ctrig,'sim').eff_dist_h 629 elif ifo == 'L1': 630 effDs_inj[ifo] = getattr(ctrig,'sim').eff_dist_l 631 elif ifo == 'V1': 632 effDs_inj[ifo] = getattr(ctrig,'sim').eff_dist_v 633 dist_inj = getattr(ctrig,'sim').distance 634 coinc.set_inj_params(getattr(ctrig,'sim').latitude,getattr(ctrig,'sim').longitude, \ 635 getattr(ctrig,'sim').mass1,getattr(ctrig,'sim').mass2,dist_inj,effDs_inj) 636 coinc.is_injection = True 637 #FIXME: name the exception! 638 except: 639 pass 640 641 self.append(coinc)

642 643 ############################################################################## 644 # 645 # table definitions and functions for populating them 646 # 647 ############################################################################## 648

649 -class SkyLocTable(tab.Table):

650 tableName = "SkyLoc" 651 validcolumns = { 652 "end_time": "int_4s", 653 "comb_snr": "real_4", 654 "ifos": "lstring", 655 "ra": "real_4", 656 "dec": "real_4", 657 "dt10": "real_4", 658 "dt20": "real_4", 659 "dt30": "real_4", 660 "dt40": "real_4", 661 "dt50": "real_4", 662 "dt60": "real_4", 663 "dt70": "real_4", 664 "dt80": "real_4", 665 "dt90": "real_4", 666 "min_eff_distance": "real_4", 667 "skymap": "lstring", 668 "galaxy_grid": "lstring", 669 "grid": "lstring" 670 }

671

672 -class SkyLocRow(object):

673 __slots__ = SkyLocTable.validcolumns.keys() 674

675 - def get_ifos(self):

676 """ 677 Return a set of the instruments for this row. 678 """ 679 return lsctables.instrument_set_from_ifos(self.ifos)

680

681 - def set_ifos(self, instruments):

682 """ 683 Serialize a sequence of instruments into the ifos 684 attribute. The instrument names must not contain the "," 685 character. 686 """ 687 self.ifos = lsctables.ifos_from_instrument_set(instruments)

688 689 SkyLocTable.RowType = SkyLocRow 690

691 -class SkyLocInjTable(tab.Table):

692 tableName = "SkyLocInj" 693 validcolumns = { 694 "end_time": "int_4s", 695 "ifos": "lstring", 696 "comb_snr": "real_4", 697 "h1_snr": "real_4", 698 "l1_snr": "real_4", 699 "v1_snr": "real_4", 700 "ra": "real_4", 701 "dec": "real_4", 702 "dt_area": "real_4", 703 "rank_area": "real_4", 704 "gal_area": "real_4", 705 "delta_t_rss": "real_8", 706 "delta_D_rss": "real_8", 707 "rank": "real_8", 708 "h1_eff_distance": "real_4", 709 "l1_eff_distance": "real_4", 710 "v1_eff_distance": "real_4", 711 "mass1": "real_4", 712 "mass2": "real_4", 713 "distance": "real_4", 714 "grid": "lstring", 715 "galaxy_grid": "lstring" 716 }

717

718 -class SkyLocInjRow(object):

719 __slots__ = SkyLocInjTable.validcolumns.keys() 720

721 - def get_ifos(self):

722 """ 723 Return a set of the instruments for this row. 724 """ 725 return lsctables.instrument_set_from_ifos(self.ifos)

726

727 - def set_ifos(self, instruments):

728 """ 729 Serialize a sequence of instruments into the ifos 730 attribute. The instrument names must not contain the "," 731 character. 732 """ 733 self.ifos = lsctables.ifos_from_instrument_set(instruments)

734 735 736 SkyLocInjTable.RowType = SkyLocInjRow 737

738 -class GalaxyTable(tab.Table):

739 tableName = "Galaxy" 740 validcolumns = { 741 "end_time": "int_4s", 742 "name": "lstring", 743 "ra": "real_8", 744 "dec": "real_8", 745 "distance_kpc": "real_8", 746 "distance_error": "real_8", 747 "luminosity_mwe": "real_8", 748 "metal_correction": "real_8", 749 "magnitude_error": "real_8" 750 }

751

752 -class GalaxyRow(object):

753 __slots__ = GalaxyTable.validcolumns.keys()

754 755 GalaxyTable.RowType = GalaxyRow 756

757 -def populate_SkyLocTable(skyloctable,coinc,grid,A,grid_fname,\ 758 skymap_fname=None,galmap_fname=None):

759 """ 760 populate a row in a skyloctable 761 """ 762 row = skyloctable.RowType() 763 764 row.end_time = coinc.time 765 row.set_ifos(coinc.ifo_list) 766 rhosquared = 0.0 767 for ifo in coinc.ifo_list: 768 rhosquared += coinc.snr[ifo]*coinc.snr[ifo] 769 row.comb_snr = sqrt(rhosquared) 770 row.dec,row.ra = grid[0][0] 771 #compute areas 772 def area(sp,pct,A,n): 773 return float(len([pt for pt in sp if pt[n] >= pct/100.]))*A

774 grid.nsort(2) 775 row.dt90 = area(grid,90.,A,2) 776 row.dt80 = area(grid,80.,A,2) 777 row.dt70 = area(grid,70.,A,2) 778 row.dt60 = area(grid,60.,A,2) 779 row.dt50 = area(grid,50.,A,2) 780 row.dt40 = area(grid,40.,A,2) 781 row.dt30 = area(grid,30.,A,2) 782 row.dt20 = area(grid,20.,A,2) 783 row.dt10 = area(grid,10.,A,2) 784 grid.nsort(1) 785 row.min_eff_distance = min(effD for effD in coinc.eff_distances.values()) 786 if skymap_fname: 787 row.skymap = os.path.basename(str(skymap_fname)) 788 else: 789 row.skymap = skymap_fname 790 791 row.grid = os.path.basename(str(grid_fname)) 792 if galmap_fname: 793 row.galaxy_grid = os.path.basename(str(galmap_fname)) 794 else: 795 row.galaxy_grid = galmap_fname 796 797 skyloctable.append(row) 798

799 -def populate_SkyLocInjTable(skylocinjtable,coinc,rank,area,\ 800 dtrss_inj,dDrss_inj,grid_fname,gal_grid):

801 """ 802 record injection data in a skylocinjtable 803 """ 804 row = skylocinjtable.RowType() 805 806 row.end_time = coinc.time 807 row.set_ifos(coinc.ifo_list) 808 row.rank = rank 809 row.distance = coinc.distance_inj 810 rhosquared = 0.0 811 for ifo in coinc.ifo_list: 812 rhosquared += coinc.snr[ifo]*coinc.snr[ifo] 813 row.comb_snr = sqrt(rhosquared) 814 try: 815 row.h1_snr = coinc.snr['H1'] 816 except: 817 row.h1_snr = None 818 try: 819 row.l1_snr = coinc.snr['L1'] 820 except: 821 row.l1_snr = None 822 try: 823 row.v1_snr = coinc.snr['V1'] 824 except: 825 row.v1_snr = None 826 row.ra = coinc.longitude_inj 827 row.dec = coinc.latitude_inj 828 row.dt_area = area['dt'] 829 row.rank_area = area['rank'] 830 if area['gal']: 831 row.gal_area = area['gal'] 832 else: 833 row.gal_area = None 834 row.delta_t_rss = dtrss_inj 835 row.delta_D_rss = dDrss_inj 836 try: 837 row.h1_eff_distance = coinc.eff_distances_inj['H1'] 838 except: 839 row.h1_eff_distance = None 840 try: 841 row.l1_eff_distance = coinc.eff_distances_inj['L1'] 842 except: 843 row.l1_eff_distance = None 844 try: 845 row.v1_eff_distance = coinc.eff_distances_inj['V1'] 846 except: 847 row.v1_eff_distance = None 848 row.mass1 = coinc.mass1_inj 849 row.mass2 = coinc.mass2_inj 850 row.grid = os.path.basename(str(grid_fname)) 851 if gal_grid: 852 row.galaxy_grid = os.path.basename(str(gal_grid)) 853 else: 854 row.galaxy_grid = gal_grid 855 skylocinjtable.append(row)

856

857 -def populate_GalaxyTable(galaxytable,coinc,galaxy):

858 """ 859 record galaxy data in a galaxytable 860 """ 861 row = galaxytable.RowType() 862 863 row.end_time = coinc.time 864 row.name = galaxy.name 865 row.ra = galaxy.ra 866 row.dec = galaxy.dec 867 row.distance_kpc = galaxy.distance_kpc 868 row.luminosity_mwe = galaxy.luminosity_mwe 869 row.magnitude_error = galaxy.magnitude_error 870 row.distance_error = galaxy.distance_error 871 row.metal_correction = galaxy.metal_correction 872 873 galaxytable.append(row)

874

Source Code for Module pylal.skylocutils