Package aloha :: Package template_files :: Module wavefunctions
[hide private]
[frames] | no frames]

Source Code for Module aloha.template_files.wavefunctions

  1  from __future__ import division  
  2  import math 
  3  from math import sqrt, pow 
  4   
5 -class WaveFunction(list):
6 """a objet for a WaveFunction""" 7 8 spin_to_size={0:1, 9 1:3, 10 2:6, 11 3:6, 12 5:18} 13
14 - def __init__(self, spin= None, size=None):
15 """Init the list with zero value""" 16 17 if spin: 18 size = self.spin_to_size[spin] 19 list.__init__(self, [0]*size)
20 21
22 -def ixxxxx(p,fmass,nhel,nsf):
23 """Defines an inflow fermion.""" 24 25 fi = WaveFunction(2) 26 27 fi[4] = complex(p[0]*nsf,p[3]*nsf) 28 fi[5] = complex(p[1]*nsf,p[2]*nsf) 29 30 nh = nhel*nsf 31 32 if (fmass != 0.): 33 pp = min(p[0],sqrt(p[1]**2 + p[2]**2 + p[3]**2 )) 34 if (pp == 0.): 35 sqm = sqrt(abs(fmass)) 36 ip = (1+nh)/2 37 im = (1-nh)/2 38 39 fi[0] = ip*sqm 40 fi[1] = im*nsf*sqm 41 fi[2] = ip*nsf*sqm 42 fi[3] = im*sqm 43 44 else: 45 sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5] 46 omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))] 47 ip = (1+nh)//2 48 im = (1-nh)//2 49 sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]] 50 pp3 = max(pp+p[3],0.) 51 if (pp3 == 0.): 52 chi1 = complex(-nh,0.) 53 else: 54 chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\ 55 p[2]/sqrt(2.*pp*pp3)) 56 chi = [complex(sqrt(pp3*0.5/pp)),chi1] 57 58 fi[0] = sfomeg[0]*chi[im] 59 fi[1] = sfomeg[0]*chi[ip] 60 fi[2] = sfomeg[1]*chi[im] 61 fi[3] = sfomeg[1]*chi[ip] 62 63 else: 64 sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf 65 if (sqp0p3 == 0.): 66 chi1 = complex(-nhel*sqrt(2.*p[0]),0.) 67 else: 68 chi1 = complex(nh*p[1]/sqp0p3,p[2]/sqp0p3) 69 chi = [complex(sqp0p3,0.),chi1] 70 if (nh == 1): 71 fi[0] = complex(0.,0.) 72 fi[1] = complex(0.,0.) 73 fi[2] = chi[0] 74 fi[3] = chi[1] 75 else: 76 fi[0] = chi[1] 77 fi[1] = chi[0] 78 fi[2] = complex(0.,0.) 79 fi[3] = complex(0.,0.) 80 81 return fi
82
83 -def oxxxxx(p,fmass,nhel,nsf):
84 """ initialize an outgoing fermion""" 85 86 fo = WaveFunction(2) 87 88 fo[4] = complex(p[0]*nsf,p[3]*nsf) 89 fo[5] = complex(p[1]*nsf,p[2]*nsf) 90 91 nh = nhel*nsf 92 93 if (fmass != 0.): 94 pp = min(p[0],sqrt(p[1]**2 + p[2]**2 + p[3]**2 )) 95 if (pp == 0.): 96 sqm = sqrt(abs(fmass)) 97 ip = -((1-nh)/2) * nhel 98 im = (1+nh)/2 * nhel 99 100 fo[0] = im*sqm 101 fo[1] = ip*nsf*sqm 102 fo[2] = im*nsf*sqm 103 fo[3] = ip*sqm 104 105 else: 106 sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5] 107 omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))] 108 ip = (1+nh)//2 109 im = (1-nh)//2 110 sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]] 111 pp3 = max(pp+p[3],0.) 112 if (pp3 == 0.): 113 chi1 = complex(-nh,0.) 114 else: 115 chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\ 116 -p[2]/sqrt(2.*pp*pp3)) 117 chi = [complex(sqrt(pp3*0.5/pp)),chi1] 118 119 fo[0] = sfomeg[1]*chi[im] 120 fo[1] = sfomeg[1]*chi[ip] 121 fo[2] = sfomeg[0]*chi[im] 122 fo[3] = sfomeg[0]*chi[ip] 123 124 else: 125 sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf 126 if (sqp0p3 == 0.): 127 chi1 = complex(-nhel*sqrt(2.*p[0]),0.) 128 else: 129 chi1 = complex(nh*p[1]/sqp0p3,-p[2]/sqp0p3) 130 chi = [complex(sqp0p3,0.),chi1] 131 if (nh == 1): 132 fo[0] = chi[0] 133 fo[1] = chi[1] 134 fo[2] = complex(0.,0.) 135 fo[3] = complex(0.,0.) 136 else: 137 fo[0] = complex(0.,0.) 138 fo[1] = complex(0.,0.) 139 fo[2] = chi[1] 140 fo[3] = chi[0] 141 142 return fo
143
144 -def vxxxxx(p,vmass,nhel,nsv):
145 """ initialize a vector wavefunction. nhel=4 is for checking BRST""" 146 147 vc = WaveFunction(3) 148 149 sqh = sqrt(0.5) 150 nsvahl = nsv*abs(nhel) 151 pt2 = p[1]**2 + p[2]**2 152 pp = min(p[0],sqrt(pt2 + p[3]**2)) 153 pt = min(pp,sqrt(pt2)) 154 155 vc[4] = complex(p[0]*nsv,p[3]*nsv) 156 vc[5] = complex(p[1]*nsv,p[2]*nsv) 157 158 if (nhel == 4): 159 if (vmass == 0.): 160 vc[0] = 1. 161 vc[1]=p[1]/p[0] 162 vc[2]=p[2]/p[0] 163 vc[3]=p[3]/p[0] 164 else: 165 vc[0] = p[0]/vmass 166 vc[1] = p[1]/vmass 167 vc[2] = p[2]/vmass 168 vc[3] = p[3]/vmass 169 170 return vc 171 172 if (vmass != 0.): 173 hel0 = 1.-abs(nhel) 174 175 if (pp == 0.): 176 vc[0] = complex(0.,0.) 177 vc[1] = complex(-nhel*sqh,0.) 178 vc[2] = complex(0.,nsvahl*sqh) 179 vc[3] = complex(hel0,0.) 180 181 else: 182 emp = p[0]/(vmass*pp) 183 vc[0] = complex(hel0*pp/vmass,0.) 184 vc[3] = complex(hel0*p[3]*emp+nhel*pt/pp*sqh) 185 if (pt != 0.): 186 pzpt = p[3]/(pp*pt)*sqh*nhel 187 vc[1] = complex(hel0*p[1]*emp-p[1]*pzpt, \ 188 -nsvahl*p[2]/pt*sqh) 189 vc[2] = complex(hel0*p[2]*emp-p[2]*pzpt, \ 190 nsvahl*p[1]/pt*sqh) 191 else: 192 vc[1] = complex(-nhel*sqh,0.) 193 vc[2] = complex(0.,nsvahl*sign(sqh,p[3])) 194 else: 195 pp = p[0] 196 pt = sqrt(p[1]**2 + p[2]**2) 197 vc[0] = complex(0.,0.) 198 vc[3] = complex(nhel*pt/pp*sqh) 199 if (pt != 0.): 200 pzpt = p[3]/(pp*pt)*sqh*nhel 201 vc[1] = complex(-p[1]*pzpt,-nsv*p[2]/pt*sqh) 202 vc[2] = complex(-p[2]*pzpt,nsv*p[1]/pt*sqh) 203 else: 204 vc[1] = complex(-nhel*sqh,0.) 205 vc[2] = complex(0.,nsv*sign(sqh,p[3])) 206 207 return vc
208
209 -def sign(x,y):
210 """Fortran's sign transfer function""" 211 if (y < 0.): 212 return -abs(x) 213 else: 214 return abs(x)
215 216
217 -def sxxxxx(p,nss):
218 """initialize a scalar wavefunction""" 219 220 sc = WaveFunction(1) 221 222 sc[0] = complex(1.,0.) 223 sc[1] = complex(p[0]*nss,p[3]*nss) 224 sc[2] = complex(p[1]*nss,p[2]*nss) 225 return sc
226 227
228 -def txxxxx(p, tmass, nhel, nst):
229 """ initialize a tensor wavefunction""" 230 231 tc = WaveFunction(5) 232 233 sqh = sqrt(0.5) 234 sqs = sqrt(1/6) 235 236 pt2 = p[1]**2 + p[2]**2 237 pp = min(p[0],sqrt(pt2+p[3]**2)) 238 pt = min(pp,sqrt(pt2)) 239 240 ft = {} 241 ft[(4,0)] = complex(p[0], p[3]) * nst 242 ft[(5,0)] = complex(p[1], p[2]) * nst 243 244 if ( nhel >= 0 ): 245 #construct eps+ 246 ep = [0] * 4 247 248 if ( pp == 0 ): 249 #ep[0] = 0 250 ep[1] = -sqh 251 ep[2] = complex(0, nst*sqh) 252 #ep[3] = 0 253 else: 254 #ep[0] = 0 255 ep[3] = pt/pp*sqh 256 if (pt != 0): 257 pzpt = p[3]/(pp*pt)*sqh 258 ep[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh ) 259 ep[2] = complex( -p[2]*pzpt , nst*p[1]/pt*sqh ) 260 else: 261 ep[1] = -sqh 262 ep[2] = complex( 0 , nst*sign(sqh,p[3]) ) 263 264 265 266 if ( nhel <= 0 ): 267 #construct eps- 268 em = [0] * 4 269 if ( pp == 0 ): 270 #em[0] = 0 271 em[1] = sqh 272 em[2] = cpmplex( 0 , nst*sqh ) 273 #em[3] = 0 274 else: 275 #em[0] = 0 276 em[3] = -pt/pp*sqh 277 if pt: 278 pzpt = -p[3]/(pp*pt)*sqh 279 em[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh ) 280 em[2] = complex( -p[2]*pzpt , nst*p[1]/pt*sqh ) 281 else: 282 em[1] = sqh 283 em[2] = complex( 0 , nst*sign(sqh,p[3]) ) 284 285 286 if ( abs(nhel) <= 1 ): 287 #construct eps0 288 e0 = [0] * 4 289 if ( pp == 0 ): 290 #e0[0] = dcmplx( rZero ) 291 #e0[1] = dcmplx( rZero ) 292 #e0[2] = dcmplx( rZero ) 293 e0[3] = 1 294 else: 295 emp = p[0]/(tmass*pp) 296 e0[0] = pp/tmass 297 e0[3] = p[3]*emp 298 if pt: 299 e0[1] = p[1]*emp 300 e0[2] = p[2]*emp 301 #else: 302 # e0[1] = dcmplx( rZero ) 303 # e0[2] = dcmplx( rZero ) 304 305 if nhel == 2: 306 for j in range(4): 307 for i in range(4): 308 ft[(i,j)] = ep[i]*ep[j] 309 elif nhel == -2: 310 for j in range(4): 311 for i in range(4): 312 ft[(i,j)] = em[i]*em[j] 313 elif tmass == 0: 314 for j in range(4): 315 for i in range(4): 316 ft[(i,j)] = 0 317 elif nhel == 1: 318 for j in range(4): 319 for i in range(4): 320 ft[(i,j)] = sqh*( ep[i]*e0[j] + e0[i]*ep[j] ) 321 elif nhel == 0: 322 for j in range(4): 323 for i in range(4): 324 ft[(i,j)] = sqs*( ep[i]*em[j] + em[i]*ep[j] + 2 *e0[i]*e0[j] ) 325 elif nhel == -1: 326 for j in range(4): 327 for i in range(4): 328 ft[(i,j)] = sqh*( em[i]*e0[j] + e0[i]*em[j] ) 329 330 else: 331 raise Exception, 'invalid helicity TXXXXXX' 332 333 tc[0] = ft[(0,0)] 334 tc[1] = ft[(0,1)] 335 tc[2] = ft[(0,2)] 336 tc[3] = ft[(0,3)] 337 tc[4] = ft[(1,0)] 338 tc[5] = ft[(1,1)] 339 tc[6] = ft[(1,2)] 340 tc[7] = ft[(1,3)] 341 tc[8] = ft[(2,0)] 342 tc[9] = ft[(2,1)] 343 tc[10] = ft[(2,2)] 344 tc[11] = ft[(2,3)] 345 tc[12] = ft[(3,0)] 346 tc[13] = ft[(3,1)] 347 tc[14] = ft[(3,2)] 348 tc[15] = ft[(3,3)] 349 tc[16] = ft[(4,0)] 350 tc[17] = ft[(5,0)] 351 352 return tc
353