1
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
34
35
36
37
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
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
52
53
54
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
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
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
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
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
167 points.append((0.0-dlat,0.0))
168
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
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
226 """
227 gaussian kernel for kernel density estimator
228 """
229 a = x-xn
230 return exp(-a*a/(2.*w*w))
231
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
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
258
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
272
273 ds = pi*res/180.0
274
275 xindex = bisect(bins,pt)
276 xminus = xindex - 1
277 newbins = []
278 while 1:
279
280
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
296 newbins.sort()
297 rpt = (pt[1],pt[0])
298 yindex = bisect(newbins,rpt)
299 finalbins = {}
300
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
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
354
355
356
358 """
359 class for storing pdfs
360 """
362 """
363 storing pdfs
364 """
365 self.x = xvals
366 self.y = yvals
367
369 """
370 return the probability of value as obtained via linear interpolation
371 """
372 return interp(value,self.x,self.y)
373
374
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
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
457 """
458 simple container for the information needed to run the sky localization code
459 """
461 """
462 here are all the things we need
463 """
464
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
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
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
497
499 self.gps = gpsdict
500 self.time = min(t for t in self.gps.values())
501
503 self.eff_distances = effDdict
504
506 self.mass1 = m1
507 self.mass2 = m2
508
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
522
523
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'):
541
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
555
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
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
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
592 except:
593 pass
594
595 self.append(coinc)
596
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
611 except:
612 pass
613
614
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
638 except:
639 pass
640
641 self.append(coinc)
642
643
644
645
646
647
648
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
673 __slots__ = SkyLocTable.validcolumns.keys()
674
676 """
677 Return a set of the instruments for this row.
678 """
679 return lsctables.instrument_set_from_ifos(self.ifos)
680
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
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
719 __slots__ = SkyLocInjTable.validcolumns.keys()
720
722 """
723 Return a set of the instruments for this row.
724 """
725 return lsctables.instrument_set_from_ifos(self.ifos)
726
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
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
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
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
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
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