pylal.geom_aligned_bank

1 from __future__ import division 2 import os,sys,math,copy 3 import numpy 4 import matplotlib 5 matplotlib.use('Agg') 6 import pylab 7 from lal import PI as LAL_PI 8 from lal import MTSUN_SI as LAL_MTSUN_SI 9 10 # FIXME: Use lal/pylals own variables here 11 mtsun = LAL_MTSUN_SI 12 LAL_GAMMA = 0.5772156649015328606065120900824024 13 14 # This function is taken from Stackoverflow: 15 # http://stackoverflow.com/questions/377017/test-if-executable-exists-in-python/377028#377028

16 -def which(program):

17 def is_exe(fpath): 18 return os.path.isfile(fpath) and os.access(fpath, os.X_OK)

19 20 fpath, fname = os.path.split(program) 21 if fpath: 22 if is_exe(program): 23 return program 24 else: 25 for path in os.environ["PATH"].split(os.pathsep): 26 exe_file = os.path.join(path, program) 27 if is_exe(exe_file): 28 return exe_file 29 30 return None 31

32 -def determine_eigen_directions(psd,order,f0,f_low,f_upper,delta_f,\ 33 return_gs=False,verbose=False,elapsed_time=None,\ 34 vary_fmax=False,vary_density=25):

35 36 evals = {} 37 evecs = {} 38 metric = {} 39 if verbose: 40 print >>sys.stdout,"Beginning to calculate moments at %d." %(elapsed_time()) 41 42 moments = get_moments(psd,f0,f_low,f_upper,delta_f,vary_fmax=vary_fmax,\ 43 vary_density=vary_density) 44 45 if verbose: 46 print >>sys.stdout,"Moments calculated, transforming to metric at %d." \ 47 %(elapsed_time()) 48 49 list = [] 50 if vary_fmax: 51 for t_fmax in numpy.arange(f_low+vary_density,f_upper,vary_density): 52 list.append(t_fmax) 53 else: 54 list.append('fixed') 55 56 for item in list: 57 Js = [] 58 for i in range(18): 59 Js.append(moments['J%d'%(i)][item]) 60 Js.append(moments['J%d'%(-1)][item]) 61 62 logJs = [] 63 for i in range(18): 64 logJs.append(moments['log%d'%(i)][item]) 65 logJs.append(moments['log%d'%(-1)][item]) 66 67 loglogJs = [] 68 for i in range(18): 69 loglogJs.append(moments['loglog%d'%(i)][item]) 70 loglogJs.append(moments['loglog%d'%(-1)][item]) 71 72 logloglogJs = [] 73 for i in range(18): 74 logloglogJs.append(moments['logloglog%d'%(i)][item]) 75 logloglogJs.append(moments['logloglog%d'%(-1)][item]) 76 77 loglogloglogJs = [] 78 for i in range(18): 79 loglogloglogJs.append(moments['loglogloglog%d'%(i)][item]) 80 loglogloglogJs.append(moments['loglogloglog%d'%(-1)][item]) 81 82 mapping = {} 83 84 if order == 'twoPN': 85 maxLen = 4 86 gs = numpy.matrix(numpy.zeros(shape=(maxLen,maxLen),dtype=float)) 87 mapping['Lambda0'] = 0 88 mapping['Lambda2'] = 1 89 mapping['Lambda3'] = 2 90 mapping['Lambda4'] = 3 91 elif order == 'threePointFivePN': 92 maxLen = 8 93 gs = numpy.matrix(numpy.zeros(shape=(maxLen,maxLen),dtype=float)) 94 mapping['Lambda0'] = 0 95 mapping['Lambda2'] = 1 96 mapping['Lambda3'] = 2 97 mapping['Lambda4'] = 3 98 mapping['LogLambda5'] = 4 99 mapping['Lambda6'] = 5 100 mapping['Lambda7'] = 6 101 mapping['LogLambda6'] = 7 102 elif order == 'taylorF4_45PN': 103 maxLen = 12 104 gs = numpy.matrix(numpy.zeros(shape=(maxLen,maxLen),dtype=float)) 105 mapping['Lambda0'] = 0 106 mapping['Lambda2'] = 1 107 mapping['Lambda3'] = 2 108 mapping['Lambda4'] = 3 109 mapping['LogLambda5'] = 4 110 mapping['Lambda6'] = 5 111 mapping['Lambda7'] = 6 112 mapping['LogLambda6'] = 7 113 mapping['LogLambda8'] = 8 114 mapping['LogLogLambda8'] = 9 115 mapping['Lambda9'] = 10 116 mapping['LogLambda9'] = 11 117 else: 118 raise BrokenError 119 120 for i in range(16): 121 for j in range(16): 122 # Normal terms 123 if mapping.has_key('Lambda%d'%i) and mapping.has_key('Lambda%d'%j): 124 gs[mapping['Lambda%d'%i],mapping['Lambda%d'%j]] = 0.5 * (Js[17-i-j] - Js[12-i]*Js[12-j] - (Js[9-i] - Js[4]*Js[12-i]) * (Js[9-j] - Js[4] * Js[12-j])/(Js[1] - Js[4]*Js[4])) 125 # Normal,log cross terms 126 if mapping.has_key('Lambda%d'%i) and mapping.has_key('LogLambda%d'%j): 127 gammaij = logJs[17-i-j] - logJs[12-j] * Js[12-i] 128 gamma0i = (Js[9-i] - Js[4] * Js[12-i]) 129 gamma0j = logJs[9-j] - logJs[12-j] * Js[4] 130 gs[mapping['Lambda%d'%i],mapping['LogLambda%d'%j]] = \ 131 gs[mapping['LogLambda%d'%j],mapping['Lambda%d'%i]] = \ 132 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4])) 133 # Log,log terms 134 if mapping.has_key('LogLambda%d'%i) and mapping.has_key('LogLambda%d'%j): 135 gammaij = loglogJs[17-i-j] - logJs[12-j] * logJs[12-i] 136 gamma0i = (logJs[9-i] - Js[4] * logJs[12-i]) 137 gamma0j = logJs[9-j] - logJs[12-j] * Js[4] 138 gs[mapping['LogLambda%d'%i],mapping['LogLambda%d'%j]] = \ 139 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4])) 140 # Normal,loglog cross terms 141 if mapping.has_key('Lambda%d'%i) and mapping.has_key('LogLogLambda%d'%j): 142 gammaij = loglogJs[17-i-j] - loglogJs[12-j] * Js[12-i] 143 gamma0i = (Js[9-i] - Js[4] * Js[12-i]) 144 gamma0j = loglogJs[9-j] - loglogJs[12-j] * Js[4] 145 gs[mapping['Lambda%d'%i],mapping['LogLogLambda%d'%j]] = \ 146 gs[mapping['LogLogLambda%d'%j],mapping['Lambda%d'%i]] = \ 147 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4])) 148 # log,loglog cross terms 149 if mapping.has_key('LogLambda%d'%i) and mapping.has_key('LogLogLambda%d'%j): 150 gammaij = logloglogJs[17-i-j] - loglogJs[12-j] * logJs[12-i] 151 gamma0i = (logJs[9-i] - Js[4] * logJs[12-i]) 152 gamma0j = loglogJs[9-j] - loglogJs[12-j] * Js[4] 153 gs[mapping['LogLambda%d'%i],mapping['LogLogLambda%d'%j]] = \ 154 gs[mapping['LogLogLambda%d'%j],mapping['LogLambda%d'%i]] = \ 155 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4])) 156 # Loglog,loglog terms 157 if mapping.has_key('LogLogLambda%d'%i) and mapping.has_key('LogLogLambda%d'%j): 158 gammaij = loglogloglogJs[17-i-j] - loglogJs[12-j] * loglogJs[12-i] 159 gamma0i = (loglogJs[9-i] - Js[4] * loglogJs[12-i]) 160 gamma0j = loglogJs[9-j] - loglogJs[12-j] * Js[4] 161 gs[mapping['LogLogLambda%d'%i],mapping['LogLogLambda%d'%j]] = \ 162 0.5 * (gammaij - gamma0i*gamma0j/(Js[1] - Js[4]*Js[4])) 163 164 evals[item],evecs[item] = numpy.linalg.eig(gs) 165 metric[item] = numpy.matrix(gs) 166 167 for i in range(len(evals[item])): 168 if evals[item][i] < 0: 169 print "WARNING: Negative eigenvalue %e. Setting as positive." %(evals[item][i]) 170 evals[item][i] = -evals[item][i] 171 if evecs[item][i,i] < 0: 172 # We demand a convention that all diagonal terms in the matrix 173 # of eigenvalues are positive. 174 evecs[item][:,i] = - evecs[item][:,i] 175 176 if verbose: 177 print >>sys.stdout,"Metric and eigenvalues calculated at %d." \ 178 %(elapsed_time()) 179 180 if return_gs: 181 return evals,evecs,gs 182 183 return evals,evecs

184

185 -def get_moments(psd_file,f0,f_low,f_high,deltaF,vary_fmax=False,vary_density=25):

186 psd = numpy.loadtxt(psd_file) 187 psd_f = psd[:,0] 188 psd_amp = psd[:,1] 189 psd_amp = psd_amp * psd_amp 190 new_f,new_amp = interpolate_psd(psd_f,psd_amp,deltaF) 191 192 # Need I7 193 funct = lambda x: 1 194 I7 = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,vary_fmax=vary_fmax,vary_density=vary_density) 195 196 # Do all the J moments 197 moments = {} 198 for i in range(-1,18): 199 funct = lambda x: x**((-i+7)/3.) 200 moments['J%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density) 201 202 # Do the logx multiplied by some power terms 203 for i in range(-1,18): 204 funct = lambda x: (numpy.log(x**(1./3.))) * x**((-i+7)/3.) 205 moments['log%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density) 206 207 # Do the loglog term 208 for i in range(-1,18): 209 funct = lambda x: (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * x**((-i+7)/3.) 210 moments['loglog%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density) 211 212 # Do the logloglog term 213 for i in range(-1,18): 214 funct = lambda x: (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * x**((-i+7)/3.) 215 moments['logloglog%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density) 216 217 # Do the logloglog term 218 for i in range(-1,18): 219 funct = lambda x: (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * (numpy.log(x**(1./3.))) * x**((-i+7)/3.) 220 moments['loglogloglog%d' %(i)] = calculate_moment(new_f,new_amp,f_low,f_high,f0,funct,norm=I7,vary_fmax=vary_fmax,vary_density=vary_density) 221 222 return moments

223

224 -def interpolate_psd(psd_f,psd_amp,deltaF):

225 new_psd_f = [] 226 new_psd_amp = [] 227 fcurr = psd_f[0] 228 229 for i in range(len(psd_f) - 1): 230 f_low = psd_f[i] 231 f_high = psd_f[i+1] 232 amp_low = psd_amp[i] 233 amp_high = psd_amp[i+1] 234 while(1): 235 if fcurr > f_high: 236 break 237 new_psd_f.append(fcurr) 238 gradient = (amp_high - amp_low) / (f_high - f_low) 239 fDiff = fcurr - f_low 240 new_psd_amp.append(amp_low + fDiff * gradient) 241 fcurr = fcurr + deltaF 242 return numpy.asarray(new_psd_f),numpy.asarray(new_psd_amp)

243 244

245 -def calculate_moment(psd_f,psd_amp,fmin,fmax,f0,funct,norm=None,vary_fmax=False,vary_density=25):

246 # Must ensure deltaF in psd_f is constant 247 psd_x = psd_f / f0 248 deltax = psd_x[1] - psd_x[0] 249 250 comps = (psd_x)**(-7./3.) * funct(psd_x) * deltax/ psd_amp 251 moment = {} 252 logica = numpy.logical_and(psd_f > fmin, psd_f < fmax) 253 comps_red = comps[logica] 254 psdf_red = psd_f[logica] 255 moment['fixed'] = comps_red.sum() 256 if norm: 257 moment['fixed'] = moment['fixed']/norm['fixed'] 258 if vary_fmax: 259 for t_fmax in numpy.arange(fmin+vary_density,fmax,vary_density): 260 comps_red2 = comps_red[psdf_red < t_fmax] 261 moment[t_fmax] = comps_red2.sum() 262 if norm: 263 moment[t_fmax] = moment[t_fmax]/norm[t_fmax] 264 return moment

265

266 -def estimate_mass_range_slimline(numiters,order,evals,evecs,maxmass1,minmass1,maxmass2,minmass2,maxspin,f0,covary=True,maxBHspin=None,evecsCV=None,vary_fmax=False,maxmass=None):

267 out = [] 268 valsF = get_random_mass_slimline(numiters,minmass1,maxmass1,minmass2,maxmass2,maxspin,maxBHspin = maxBHspin,return_spins=True,maxmass=maxmass) 269 valsF = numpy.array(valsF) 270 mass = valsF[0] 271 eta = valsF[1] 272 beta = valsF[2] 273 sigma = valsF[3] 274 gamma = valsF[4] 275 chis = 0.5*(valsF[5] + valsF[6]) 276 if covary: 277 lambdas = get_cov_params(mass,eta,beta,sigma,gamma,chis,f0,evecs,evals,evecsCV,order) 278 else: 279 lambdas = get_conv_params(mass,eta,beta,sigma,gamma,chis,f0,evecs,evals,order,vary_fmax=vary_fmax) 280 281 return numpy.array(lambdas)

282

283 -def get_random_mass_slimline(N,minmass1,maxmass1,minmass2,maxmass2,maxspin,maxBHspin = None,return_spins=False,qm_scalar_fac=1,maxmass=None):

284 # WARNING: We expect mass1 > mass2 ALWAYS 285 minmass = minmass1 + minmass2 286 if not maxmass: 287 maxmass = maxmass1 + maxmass2 288 mincompmass = minmass2 289 maxcompmass = maxmass1 290 291 mass = numpy.random.random(N) * (minmass**(-5./3.)-maxmass**(-5./3.)) + maxmass**(-5./3.) 292 mass = mass**(-3./5.) 293 maxmass2 = numpy.minimum(mass/2.,maxmass2) 294 minmass1 = numpy.maximum(minmass1,mass/2.) 295 mineta = numpy.maximum(mincompmass * (mass-mincompmass)/(mass*mass), maxcompmass*(mass-maxcompmass)/(mass*mass)) 296 maxeta = numpy.minimum(0.25,maxmass2 * (mass - maxmass2) / (mass*mass)) 297 maxeta = numpy.minimum(maxeta,minmass1 * (mass - minmass1) / (mass*mass)) 298 if (maxeta < mineta).any(): 299 print "WARNING: Max eta is smaller than min eta!!" 300 eta = numpy.random.random(N) * (maxeta - mineta) + mineta 301 diff = (mass*mass * (1-4*eta))**0.5 302 mass1 = (mass + diff)/2. 303 mass2 = (mass - diff)/2. 304 if (mass1 > maxmass1).any() or (mass1 < minmass1).any(): 305 print "WARNING: Mass1 outside of mass range" 306 if (mass2 > maxmass2).any() or (mass2 < minmass2).any(): 307 print "WARNING: Mass1 outside of mass range" 308 if maxspin > 0: 309 mspin = numpy.zeros(len(mass1)) 310 mspin += maxspin 311 if maxBHspin: 312 mspin[mass1 > 3] = maxBHspin 313 spin1z = numpy.random.random(N) * mspin*2 - mspin 314 mspin = numpy.zeros(len(mass2)) 315 mspin += maxspin 316 if maxBHspin: 317 mspin[mass2 > 3] = maxBHspin 318 spin2z = numpy.random.random(N) * mspin*2 - mspin 319 320 spinspin = spin1z*spin2z 321 else: 322 spinspin = numpy.zeros(N,dtype=float) 323 spin1z = numpy.zeros(N,dtype=float) 324 spin2z = numpy.zeros(N,dtype=float) 325 326 chiS = 0.5 * (spin1z + spin2z) 327 chiA = 0.5 * (spin1z - spin2z) 328 delta = (mass1 - mass2) / (mass1 + mass2) 329 330 beta = (113. / 12. - 19./3. * eta) * chiS 331 beta += 113. / 12. * delta * chiA 332 sigma = eta / 48. * (474 * spinspin) 333 gamma = numpy.zeros(len(sigma)) 334 for spinA,massA in zip([spin1z,spin2z],[mass1,mass2]): 335 sigmaFac = 1. / 96. * (massA / mass)**2 336 sigmaFac2 = (720 * qm_scalar_fac -1) * spinA * spinA 337 sigmaFac3 = (240 * qm_scalar_fac - 7) * spinA * spinA 338 sigma += sigmaFac * (sigmaFac2 - sigmaFac3) 339 gamma = (732985./2268. - 24260./81.*eta - 340./9.*eta*eta)*chiS 340 gamma += (732985. / 2268. + 140./9. * eta) * delta * chiA 341 342 if return_spins: 343 return mass,eta,beta,sigma,gamma,spin1z,spin2z 344 else: 345 return mass,eta,beta,sigma,gamma

346

347 -def rotate_vector(evecs,old_vector,rescale_factor,index,length):

348 temp = 0 349 for i in range(length): 350 temp += evecs[i,index] * old_vector[i] 351 temp *= rescale_factor 352 return temp

353

354 -def rotate_vector_inv(evecs,old_vector,rescale_factor,index,length):

355 temp = 0 356 for i in range(length): 357 temp += evecs[index,i] * old_vector[i] 358 temp *= rescale_factor 359 return temp

360

361 -def get_conv_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,order,vary_fmax=False):

362 363 lambdas = get_chirp_params(totmass,eta,beta,sigma,gamma,chis,f0,order) 364 365 lams = [] 366 if not vary_fmax: 367 length = len(evals) 368 for i in range(length): 369 lams.append(rotate_vector(evecs,lambdas,math.sqrt(evals[i]),i,length)) 370 return lams 371 else: 372 # Get the frequencies in the evecs/evals 373 fs = numpy.array(evals.keys(),dtype=float) 374 fs.sort() 375 # Get the frequencies of the input 376 fISCO = (1/6.)**(3./2.) / (LAL_PI * totmass * LAL_MTSUN_SI) 377 378 # INitialize output 379 length = len(evals[fs[0]]) 380 output=numpy.zeros([length,len(totmass)]) 381 lambdas = numpy.array(lambdas) 382 # We assume that the evecs are sampled at equal frequencies 383 for i in range(len(fs)): 384 if (i == 0): 385 logicArr = fISCO < ((fs[0] + fs[1])/2.) 386 if (i == (len(fs)-1)): 387 logicArr = fISCO > ((fs[-1] + fs[-1])/2.) 388 else: 389 logicArrA = fISCO > ((fs[i-1] + fs[i])/2.) 390 logicArrB = fISCO < ((fs[i] + fs[i+1])/2.) 391 logicArr = numpy.logical_and(logicArrA,logicArrB) 392 if logicArr.any(): 393 for j in range(length): 394 output[j,logicArr] = rotate_vector(evecs[fs[i]],lambdas[:,logicArr],math.sqrt(evals[fs[i]][j]),j,length) 395 # For now a list of arrays is returned so we convert 396 for i in range(length): 397 lams.append(output[i]) 398 return lams

399

400 -def get_chi_params(lambdas,f0,evecs,evals,order):

401 lams = [] 402 length = len(evals) 403 for i in range(length): 404 lams.append(rotate_vector(evecs,lambdas,math.sqrt(evals[i]),i,length)) 405 return lams

406

407 -def get_cov_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,evecsCV,order):

408 mus = get_conv_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,order) 409 xis = get_covaried_params(mus,evecsCV) 410 return xis

411

412 -def get_covaried_params(lambdas,evecsCV):

413 length = len(evecsCV) 414 lams = [] 415 for i in range(length): 416 lams.append(rotate_vector(evecsCV,lambdas,1.,i,length)) 417 return lams

418

419 -def get_chirp_params(totmass,eta,beta,sigma,gamma,chis,f0,order):

420 # Convert mass to seconds 421 totmass = totmass * mtsun 422 pi = numpy.pi 423 lambda0 = 3 / (128 * eta * (pi * totmass * f0)**(5/3)) 424 lambda2 = 5 / (96 * pi * eta * totmass * f0) * (743/336 + 11/4 * eta) 425 lambda3 = (-3 * pi**(1/3))/(8 * eta * (totmass*f0)**(2/3)) * (1 - beta/ (4 * pi)) 426 lambda4 = 15 / (64 * eta * (pi * totmass * f0)**(1/3)) * (3058673/1016064 + 5429/1008 * eta + 617/144 * eta**2 - sigma) 427 if order == 'twoPN': 428 return lambda0,lambda2,lambda3,lambda4 429 elif order[0:16] == 'threePointFivePN' or order[0:8] == 'taylorF4': 430 lambda5 = 3. * (38645.*pi/756. - 65.*pi*eta/9. - gamma) 431 lambda5 = lambda5 * (3./(128.*eta)) 432 lambda6 = 11583231236531./4694215680. - (640.*pi*pi)/3. - (6848.*LAL_GAMMA)/21. 433 lambda6 -= (6848./21.) * numpy.log(4 * (pi * totmass * f0)**(1./3.)) 434 lambda6 += (-15737765635/3048192. + 2255.*pi*pi/12.)*eta 435 lambda6 += (76055.*eta*eta)/1728. - (127825.*eta*eta*eta)/1296.; 436 lambda6 = lambda6 437 lambda6 = lambda6 * 3./(128.*eta) * (pi * totmass * f0)**(1/3.) 438 lambda7 = (77096675.*pi)/254016. + (378515.*pi*eta)/1512. - (74045.*pi*eta*eta)/756. 439 lambda7 = lambda7 440 lambda7 = lambda7 * 3./(128.*eta) * (pi * totmass * f0)**(2/3.) 441 lambda6log = -( 6848./21) 442 lambda6log = lambda6log * 3./(128.*eta) * (pi * totmass * f0)**(1/3.) 443 if order[0:16] == 'threePointFivePN': 444 return lambda0,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7,lambda6log 445 elif order[0:13] == 'taylorF4_45PN' or order[0:13] == 'taylorF4_35PN': 446 # Add 3PN spin corr term 447 lambda6spin = 502.6548245743669 * beta + 88.45238095238095 * sigma + \ 448 (110. * eta * sigma) - 20. * beta * beta 449 lambda6spin = lambda6spin * 3./(128.*eta) * (pi * totmass * f0)**(1/3.) 450 lambda6 += lambda6spin 451 # Add 3.5PN spin corr term 452 lambda7spin = -510.0603994773098*beta - 368.01525846326734*beta*eta + \ 453 1944.363555525455*chis*eta - 502.6548245743669*sigma + \ 454 40.*beta*sigma + 241.47615535889872*beta*eta*eta + \ 455 2961.654024441635*chis*eta*eta + 676.0619469026549*chis*eta*eta*eta 456 lambda7spin = lambda7spin * 3./(128.*eta) * (pi * totmass * f0)**(2/3.) 457 lambda7 += lambda7spin 458 # Add 4PN non spin term (this has log and loglog terms) 459 lambda8 = 342.6916926002232 + 2869.024558661873*eta - \ 460 3773.659169914512*eta*eta + 172.0609149438239*eta*eta*eta - \ 461 24.284336419753085*eta*eta*eta*eta 462 lambda8log = -1028.0750778006693 - 8607.073675985623*eta + \ 463 11320.977509743536*eta*eta - 516.1827448314717*eta*eta*eta + \ 464 72.85300925925927*eta*eta*eta*eta 465 lambda8loglog = 480.7316704459562 + 597.8412698412699*eta 466 # And the 4PN spin terms 467 lambda8spin = 936.7471880419974*beta - 311.03929625364435*beta*eta - \ 468 2455.4171568883194*chis*eta + 195.39588893571195*beta*chis*eta + \ 469 48.491201779065534*sigma + 101.92901234567901*eta*sigma - \ 470 58.81315259633844*beta*beta + \ 471 8.918387413962636*eta*beta*beta - 686.5167663065837*chis*eta*eta \ 472 + 54.631268436578175*beta*chis*eta*eta + \ 473 71.69753086419753*sigma*eta*eta - \ 474 4.444444444444445*sigma*sigma 475 lambda8logspin = -2810.241564125992*beta + 933.117888760933*beta*eta + \ 476 7366.251470664957*chis*eta - 586.1876668071359*beta*chis*eta - \ 477 145.4736053371966*sigma - 305.78703703703707*eta*sigma \ 478 + 176.4394577890153*beta*beta - \ 479 26.755162241887906*eta*beta*beta + \ 480 2059.5502989197507*chis*eta*eta - \ 481 163.89380530973452*beta*chis*eta*eta - \ 482 215.0925925925926*sigma*eta*eta + \ 483 13.333333333333334*sigma*sigma 484 # Construct the combined 4PN terms 485 lambda8 = (lambda8 + lambda8spin) * 3./(128.*eta) * \ 486 (pi * totmass * f0)**(3./3.) 487 lambda8log = (lambda8log + lambda8logspin) * 3./(128.*eta) * \ 488 (pi * totmass * f0)**(3./3.) 489 lambda8loglog = lambda8loglog * 3./(128.*eta) * \ 490 (pi * totmass * f0)**(3./3.) 491 # And now the 4.5PN non-spin term 492 lambda9 = 20021.24571514093 - 42141.161261993766*eta - \ 493 4047.211701119762*eta*eta - 2683.4848475303043*eta*eta*eta 494 lambda9log = -4097.833617482457 495 # And now the 4.5PN spin term 496 lambda9spin = 2105.9471080635244*beta + 3909.271818583914*beta*eta - \ 497 2398.354686411564*chis*eta + 1278.4225104920606*sigma - \ 498 198.6688790560472*beta*sigma + 589.0486225480862*eta*sigma - \ 499 62.43362831858406*beta*eta*sigma + \ 500 439.6407501053519*chis*eta*sigma - 376.99111843077515*beta*beta + \ 501 10.*beta*beta*beta - 202.1451909795383*beta*eta*eta - \ 502 5711.929102446965*chis*eta*eta + \ 503 122.9203539823009*chis*sigma*eta*eta - \ 504 493.00738145963066*beta*eta*eta*eta - \ 505 4955.659484448894*chis*eta*eta*eta - \ 506 991.4721607669617*chis*eta*eta*eta*eta 507 lambda9logspin = 326.0952380952381*beta 508 lambda9 = (lambda9 + lambda9spin) * 3./(128.*eta) * \ 509 (pi * totmass * f0)**(4./3.) 510 lambda9log = (lambda9log + lambda9logspin) * 3./(128.*eta) * \ 511 (pi * totmass * f0)**(4./3.) 512 if order[0:13] == 'taylorF4_45PN': 513 return lambda0,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7,lambda6log,lambda8log,lambda8loglog,lambda9,lambda9log 514 elif order[0:13] == 'taylorF4_35PN': 515 return lambda0,lambda2,lambda3,lambda4,lambda5,lambda6,lambda7,lambda6log 516 else: 517 raise BrokenError 518 else: 519 raise BrokenError

520

521 -def make_plots(a,b,c,d,aname,bname,cname,dname,paper_plots=False):

522 if paper_plots: 523 paper_plot() 524 if not os.path.isdir('plots'): 525 os.makedirs('plots') 526 vals = [a,b,c,d] 527 names = [aname,bname,cname,dname] 528 for i in range(4): 529 for j in range(i+1,4): 530 pylab.plot(vals[i],vals[j],'b.') 531 pylab.xlabel(names[i]) 532 pylab.ylabel(names[j]) 533 # if i == 0: 534 # pylab.xlim([-700,0]) 535 # if j == 1: 536 # pylab.ylim([-3.5,0.5]) 537 # if j == 2: 538 # pylab.ylim([-2,1]) 539 # if j == 3: 540 # pylab.ylim([0.,0.1]) 541 pylab.savefig('plots/%s_vs_%s.png' % (names[i],names[j])) 542 pylab.clf()

543 544 # 3D plots 545 # for i in range(4): 546 # for j in range(i+1,4): 547 # for k in range(j+1,4): 548 # fig = pylab.figure() 549 # ax = Axes3D(fig) 550 # ax.plot(vals[i],vals[j],'b.',zs=vals[k]) 551 # Axes3D.set_xlabel(ax,names[i]) 552 # Axes3D.set_ylabel(ax,names[j]) 553 # Axes3D.set_zlabel(ax,names[k]) 554 # pylab.savefig('plots/%s_vs_%s_vs_%s.png' % (names[i],names[j],names[k])) 555 # pylab.clf() 556

557 -def generate_hexagonal_lattice(maxv1,minv1,maxv2,minv2,mindist):

558 # Place first point 559 v1s = [minv1] 560 v2s = [minv2] 561 initPoint = [minv1,minv2] 562 # Place first line 563 initLine = [initPoint] 564 tmpv1 = minv1 565 while (tmpv1 < maxv1): 566 tmpv1 = tmpv1 + (3 * mindist)**(0.5) 567 initLine.append([tmpv1,minv2]) 568 v1s.append(tmpv1) 569 v2s.append(minv2) 570 initLine = numpy.array(initLine) 571 initLine2 = copy.deepcopy(initLine) 572 initLine2[:,0] += 0.5 * (3*mindist)**0.5 573 initLine2[:,1] += 1.5 * (mindist)**0.5 574 for i in xrange(len(initLine2)): 575 v1s.append(initLine2[i,0]) 576 v2s.append(initLine2[i,1]) 577 tmpv2_1 = initLine[0,1] 578 tmpv2_2 = initLine2[0,1] 579 while tmpv2_1 < maxv2 and tmpv2_2 < maxv2: 580 tmpv2_1 = tmpv2_1 + 3.0 * (mindist)**0.5 581 tmpv2_2 = tmpv2_2 + 3.0 * (mindist)**0.5 582 initLine[:,1] = tmpv2_1 583 initLine2[:,1] = tmpv2_2 584 for i in xrange(len(initLine)): 585 v1s.append(initLine[i,0]) 586 v2s.append(initLine[i,1]) 587 for i in xrange(len(initLine2)): 588 v1s.append(initLine2[i,0]) 589 v2s.append(initLine2[i,1]) 590 v1s = numpy.array(v1s) 591 v2s = numpy.array(v2s) 592 return v1s,v2s

593

594 -def generate_anstar_3d_lattice(maxv1,minv1,maxv2,minv2,maxv3,minv3,mindist):

595 import lal 596 tiling = lal.CreateFlatLatticeTiling(3) 597 lal.SetFlatLatticeConstantBound(tiling,0,minv1,maxv1) 598 lal.SetFlatLatticeConstantBound(tiling,1,minv2,maxv2) 599 lal.SetFlatLatticeConstantBound(tiling,2,minv3,maxv3) 600 lal.SetFlatLatticeGenerator(tiling,lal.AnstarLatticeGeneratorPtr) 601 # Make a 3x3 Euclidean lattice 602 a = lal.gsl_matrix(3,3) 603 a.data[0,0] = 1 604 a.data[1,1] = 1 605 a.data[2,2] = 1 606 lal.SetFlatLatticeMetric(tiling,a,mindist) 607 608 vs1 = [] 609 vs2 = [] 610 vs3 = [] 611 count = 0 612 while (lal.NextFlatLatticePoint(tiling) >= 0): 613 count += 1 614 if not (count % 100000): 615 print "Now %d points" %(count) 616 p = lal.GetFlatLatticePoint(tiling) 617 vs1.append(p.data[0]) 618 vs2.append(p.data[1]) 619 vs3.append(p.data[2]) 620 return vs1,vs2,vs3

621

622 -def get_physical_covaried_masses(xis,bestMasses,bestXis,f0,temp_number,req_match,order,evecs,evals,evecsCV,maxmass1,minmass1,maxmass2,minmass2,maxNSspin,maxBHspin,return_smaller=False,nsbh_flag=False):

623 # TUNABLE PARAMETERS GO HERE! 624 origScaleFactor = 1 625 626 # Set up 627 xi_size = len(xis) 628 scaleFactor = origScaleFactor 629 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.) 630 count = 0 631 unFixedCount = 0 632 currDist = 100000000000000000 633 while(1): 634 # if not (count % 100): 635 # print '\rTemplate %d Distance %e Count %d Scale Factor %d ' %(temp_number,currDist,count,scaleFactor), 636 # If we are a long way away we use larger jumps 637 if count: 638 if currDist > 1 and scaleFactor == origScaleFactor: 639 scaleFactor = origScaleFactor*10 640 chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis = get_mass_distribution(bestChirpmass,bestMasses[1],bestMasses[2],bestMasses[3],scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = nsbh_flag) 641 cDist = (new_xis[0] - xis[0])**2 642 for j in range(1,xi_size): 643 cDist += (new_xis[j] - xis[j])**2 644 if (cDist.min() < req_match): 645 idx = cDist.argmin() 646 scaleFactor = origScaleFactor 647 return mass1[idx],mass2[idx],spin1z[idx],spin2z[idx],count,cDist.min(),new_xis[0][idx],new_xis[1][idx],new_xis[2][idx],new_xis[3][idx] 648 if (cDist.min() < currDist): 649 idx = cDist.argmin() 650 bestMasses[0] = totmass[idx] 651 bestMasses[1] = eta[idx] 652 bestMasses[2] = spin1z[idx] 653 bestMasses[3] = spin2z[idx] 654 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.) 655 currDist = cDist.min() 656 unFixedCount = 0 657 scaleFactor = origScaleFactor 658 count += 1 659 unFixedCount += 1 660 if unFixedCount > 5000: 661 if return_smaller: 662 diff = (bestMasses[0]*bestMasses[0] * (1-4*bestMasses[1]))**0.5 663 mass1 = (bestMasses[0] + diff)/2. 664 mass2 = (bestMasses[0] - diff)/2. 665 return mass1,mass2,bestMasses[2],bestMasses[3],count,currDist,new_xis[0][0],new_xis[1][0],new_xis[2][0],new_xis[3][0] 666 # Give up 667 else: 668 raise BrokenError 669 if not unFixedCount % 100: 670 scaleFactor *= 2 671 if scaleFactor > 64: 672 scaleFactor = 1 673 # Shouldn't be here! 674 raise BrokenError

675

676 -def get_mass_distribution(bestChirpmass,bestEta,bestSpin1z,bestSpin2z,scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = False):

677 chirpmass = bestChirpmass * ( 1 - (numpy.random.random(100) - 0.5) * 0.0001 * scaleFactor) 678 chirpmass[0] = bestChirpmass 679 eta = bestEta * ( 1 - (numpy.random.random(100) - 0.5) * 0.01 * scaleFactor) 680 eta[0] = bestEta 681 # eta = bestMasses[2] + ( (numpy.random.random(100) - 0.5) * 0.01 * scaleFactor) 682 eta[eta > 0.25] = 0.25 683 eta[eta < 0.0001] = 0.0001 684 totmass = chirpmass / (eta**(3./5.)) 685 spin1z = bestSpin1z + ( (numpy.random.random(100) - 0.5) * 0.01 * scaleFactor) 686 spin1z[0] = bestSpin1z 687 spin2z = bestSpin2z + ( (numpy.random.random(100) - 0.5) * 0.01 * scaleFactor) 688 spin2z[0] = bestSpin2z 689 beta,sigma,gamma,chis = get_beta_sigma_from_aligned_spins(totmass,eta,spin1z,spin2z) 690 691 diff = (totmass*totmass * (1-4*eta))**0.5 692 mass1 = (totmass + diff)/2. 693 mass2 = (totmass - diff)/2. 694 695 numploga1 = numpy.logical_and(abs(spin1z) < maxBHspin,mass1 > 2.99) 696 if nsbh_flag: 697 numploga = numploga1 698 else: 699 numploga2 = numpy.logical_and(abs(spin1z) < maxNSspin,mass1 < 3.01) 700 numploga = numpy.logical_or(numploga1,numploga2) 701 numplogb2 = numpy.logical_and(abs(spin2z) < maxNSspin,mass2 < 3.01) 702 if nsbh_flag: 703 numplogb = numplogb2 704 else: 705 numplogb1 = numpy.logical_and(abs(spin2z) < maxBHspin,mass2 > 2.99) 706 numplogb = numpy.logical_or(numplogb1,numplogb2) 707 numplog1 = numpy.logical_and(numploga,numplogb) 708 numplog = numpy.logical_not(numplog1) 709 beta[numplog] = 0 710 sigma[numplog] = 0 711 gamma[numplog] = 0 712 chis[numplog] = 0 713 spin1z[numplog] = 0 714 spin2z[numplog] = 0 715 716 totmass[mass1 < minmass1] = 0.0001 717 totmass[mass1 > maxmass1] = 0.0001 718 totmass[mass2 < minmass2] = 0.0001 719 totmass[mass2 > maxmass2] = 0.0001 720 721 new_xis = get_cov_params(totmass,eta,beta,sigma,gamma,chis,f0,evecs,evals,evecsCV,order) 722 return chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis

723

724 -def stack_xi_direction_brute(xis,bestMasses,bestXis,f0,temp_number,direction_num,req_match,order,evecs,evals,evecsCV,maxmass1,minmass1,maxmass2,minmass2,maxNSspin,maxBHspin,nsbh_flag=False):

725 # TUNING PARAMETERS GO HERE 726 origScaleFactor = 0.8 727 numTestPoints = 3000 728 729 # Setup 730 xi_size = len(xis) 731 origMasses = copy.deepcopy(bestMasses) 732 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.) 733 count = 0 734 unFixedCount = 0 735 xi3min = 10000000000 736 xi3max = -100000000000 737 738 for i in range(numTestPoints): 739 # Evaluate upper extent of xi3 740 scaleFactor = origScaleFactor 741 chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis = get_mass_distribution(bestChirpmass,bestMasses[1],bestMasses[2],bestMasses[3],scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = nsbh_flag) 742 cDist = (new_xis[0] - xis[0])**2 743 for j in range(1,xi_size): 744 cDist += (new_xis[j] - xis[j])**2 745 redCDist = cDist[cDist < req_match] 746 redXis = (new_xis[direction_num])[cDist < req_match] 747 748 if len(redCDist): 749 new_xis[direction_num][cDist > req_match] = -10000000 750 maxXi3 = (new_xis[direction_num]).max() 751 idx = (new_xis[direction_num]).argmax() 752 if maxXi3 > xi3max: 753 xi3max = maxXi3 754 bestMasses[0] = totmass[idx] 755 bestMasses[1] = eta[idx] 756 bestMasses[2] = spin1z[idx] 757 bestMasses[3] = spin2z[idx] 758 m1 = mass1[idx] 759 m2 = mass2[idx] 760 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.) 761 762 bestMasses = copy.deepcopy(origMasses) 763 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.) 764 765 for i in range(numTestPoints): 766 # Evaluate lower extent of xi3 767 scaleFactor = origScaleFactor 768 chirpmass,totmass,eta,spin1z,spin2z,diff,mass1,mass2,beta,sigma,gamma,chis,new_xis = get_mass_distribution(bestChirpmass,bestMasses[1],bestMasses[2],bestMasses[3],scaleFactor,order,evecs,evals,evecsCV,maxmass1,minmass1,minmass2,maxmass2,maxNSspin,maxBHspin,f0,nsbh_flag = nsbh_flag) 769 770 cDist = (new_xis[0] - xis[0])**2 771 for j in range(1,xi_size): 772 cDist += (new_xis[j] - xis[j])**2 773 redCDist = cDist[cDist < req_match] 774 redXis = (new_xis[direction_num])[cDist < req_match] 775 776 if len(redCDist): 777 new_xis[direction_num][cDist > req_match] = 10000000 778 maxXi3 = (new_xis[direction_num]).min() 779 idx = (new_xis[direction_num]).argmin() 780 if maxXi3 < xi3min: 781 xi3min = maxXi3 782 bestMasses[0] = totmass[idx] 783 bestMasses[1] = eta[idx] 784 bestMasses[2] = spin1z[idx] 785 bestMasses[3] = spin2z[idx] 786 m1 = mass1[idx] 787 m2 = mass2[idx] 788 bestChirpmass = bestMasses[0] * (bestMasses[1])**(3./5.) 789 790 return xi3min,xi3max,xi3max-xi3min

791

792 -def get_beta_sigma_from_aligned_spins(mass,eta,spin1z,spin2z):

793 diff = (mass*mass * (1-4*eta))**0.5 794 mass1 = (mass + diff)/2. 795 mass2 = (mass - diff)/2. 796 chiS = 0.5 * (spin1z + spin2z) 797 chiA = 0.5 * (spin1z - spin2z) 798 delta = (mass1 - mass2) / (mass1 + mass2) 799 spinspin = spin1z*spin2z 800 801 beta = (113. / 12. - 19./3. * eta) * chiS 802 beta += 113. / 12. * delta * chiA 803 sigma = eta / 48. * (474 * spinspin) 804 sigma += (1 - 2*eta) * (81./16. * (chiS*chiS + chiA*chiA)) 805 sigma += delta * (81. / 8. * (chiS*chiA)) 806 gamma = (732985./2268. - 24260./81.*eta - 340./9.*eta*eta)*chiS 807 gamma += (732985. / 2268. + 140./9. * eta) * delta * chiA 808 return beta,sigma,gamma,chiS

809

810 -def test_point_distance(point1,point2,evals,evecs,evecsCV,order,f0,return_xis=False):

811 # Note: I think this will work if one of these inputs is an array, but not if both are 812 aMass1 = point1[0] 813 aMass2 = point1[1] 814 aSpin1 = point1[2] 815 aSpin2 = point1[3] 816 try: 817 leng = len(aMass1) 818 aArray = True 819 except: 820 aArray = False 821 822 bMass1 = point2[0] 823 bMass2 = point2[1] 824 bSpin1 = point2[2] 825 bSpin2 = point2[3] 826 try: 827 leng = len(bMass1) 828 bArray = True 829 except: 830 bArray = False 831 832 if aArray and bArray: 833 print "I cannot take point1 and point2 as arrays" 834 835 aTotMass = aMass1 + aMass2 836 aEta = (aMass1 * aMass2) / (aTotMass * aTotMass) 837 aCM = aTotMass * aEta**(3./5.) 838 839 bTotMass = bMass1 + bMass2 840 bEta = (bMass1 * bMass2) / (bTotMass * bTotMass) 841 bCM = bTotMass * bEta**(3./5.) 842 843 abeta,asigma,agamma,achis = get_beta_sigma_from_aligned_spins(aTotMass,aEta,aSpin1,aSpin2) 844 bbeta,bsigma,bgamma,bchis = get_beta_sigma_from_aligned_spins(bTotMass,bEta,bSpin1,bSpin2) 845 846 aXis = get_cov_params(aTotMass,aEta,abeta,asigma,agamma,achis,f0,evecs,evals,evecsCV,order) 847 if return_xis and not aArray: 848 xis1 = aXis 849 850 bXis = get_cov_params(bTotMass,bEta,bbeta,bsigma,bgamma,bchis,f0,evecs,evals,evecsCV,order) 851 if return_xis and not bArray: 852 xis2 = bXis 853 854 dist = (aXis[0] - bXis[0])**2 855 for i in range(1,len(aXis)): 856 dist += (aXis[i] - bXis[i])**2 857 858 if aArray and return_xis: 859 aXis = numpy.array(aXis) 860 xis1 = aXis[:,dist.argmin()] 861 if bArray and return_xis: 862 bXis = numpy.array(bXis) 863 xis2 = bXis[:,dist.argmin()] 864 865 if return_xis: 866 return xis1,xis2 867 868 return dist

869

Source Code for Module pylal.geom_aligned_bank_utils