← Back to branch summary

~chirality-flow/chiralityflow/ChiralityFlowMG

« back to all changes in this revision

Viewing changes to aloha/template_files/wavefunctions.py

  • Committer: andrew.lifson at lu
  • Date: 2021-09-02 13:57:34 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210902135734-4eybgli0iljkax9b
added fresh copy of MG5_aMC_v3.2.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
aloha/template_files/wavefunctions.py
 
1
from __future__ import division 
 
2
from __future__ import absolute_import
 
3
import math
 
4
from math import sqrt, pow
 
5
from itertools import product
 
6
from six.moves import range
 
7
 
 
8
class WaveFunction(list):
 
9
    """a objet for a WaveFunction"""
 
10
    
 
11
    spin_to_size={0:1,
 
12
                  1:3,
 
13
                  2:6,
 
14
                  3:6,
 
15
                  4:18,
 
16
                  5:18}
 
17
    
 
18
    def __init__(self, spin= None, size=None):
 
19
        """Init the list with zero value"""
 
20
        
 
21
        if spin:
 
22
            size = self.spin_to_size[spin]
 
23
        list.__init__(self, [0]*size)
 
24
        
 
25
 
 
26
def ixxxxx(p,fmass,nhel,nsf):
 
27
    """Defines an inflow fermion."""
 
28
    
 
29
    fi = WaveFunction(2)
 
30
    
 
31
    fi[0] = complex(-p[0]*nsf,-p[3]*nsf)
 
32
    fi[1] = complex(-p[1]*nsf,-p[2]*nsf) 
 
33
    nh = nhel*nsf 
 
34
    if (fmass != 0.):
 
35
        pp = min(p[0],sqrt(p[1]**2 + p[2]**2 + p[3]**2 ))
 
36
        if (pp == 0.): 
 
37
            sqm = [sqrt(abs(fmass))]
 
38
            sqm.append(sign(sqm[0],fmass)) 
 
39
            ip = (1+nh)//2 
 
40
            im = (1-nh)//2 
 
41
 
 
42
            fi[2] = ip*sqm[ip]
 
43
            fi[3] = im*nsf*sqm[ip]
 
44
            fi[4] = ip*nsf*sqm[im]
 
45
            fi[5] = im*sqm[im]
 
46
 
 
47
        else:
 
48
            sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
 
49
            omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))]
 
50
            ip = (1+nh)//2
 
51
            im = (1-nh)//2
 
52
            sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
 
53
            pp3 = max(pp+p[3],0.)
 
54
            if (pp3 == 0.):
 
55
                chi1 = complex(-nh,0.) 
 
56
            else:
 
57
                chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\
 
58
                p[2]/sqrt(2.*pp*pp3))
 
59
            chi = [complex(sqrt(pp3*0.5/pp)),chi1]
 
60
 
 
61
            fi[2] = sfomeg[0]*chi[im]
 
62
            fi[3] = sfomeg[0]*chi[ip]
 
63
            fi[4] = sfomeg[1]*chi[im]
 
64
            fi[5] = sfomeg[1]*chi[ip] 
 
65
    
 
66
    else: 
 
67
        sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf
 
68
        if (sqp0p3 == 0.):
 
69
            chi1 = complex(-nhel*sqrt(2.*p[0]),0.)
 
70
        else:
 
71
            chi1 = complex(nh*p[1]/sqp0p3,p[2]/sqp0p3)
 
72
        chi = [complex(sqp0p3,0.),chi1]
 
73
        if (nh == 1):
 
74
            fi[2] = complex(0.,0.)
 
75
            fi[3] = complex(0.,0.)
 
76
            fi[4] = chi[0]
 
77
            fi[5] = chi[1] 
 
78
        else:
 
79
            fi[2] = chi[1]
 
80
            fi[3] = chi[0]
 
81
            fi[4] = complex(0.,0.)
 
82
            fi[5] = complex(0.,0.) 
 
83
    
 
84
    return fi 
 
85
 
 
86
def oxxxxx(p,fmass,nhel,nsf):
 
87
    """ initialize an outgoing fermion"""
 
88
    
 
89
    fo = WaveFunction(2)
 
90
    fo[0] = complex(p[0]*nsf,p[3]*nsf)
 
91
    fo[1] = complex(p[1]*nsf,p[2]*nsf)
 
92
    nh = nhel*nsf
 
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
            sqm.append( sign(sqm[0], fmass)) 
 
98
            ip = int(-((1-nh)//2) * nhel)
 
99
            im = int((1+nh)//2 * nhel)
 
100
            
 
101
            fo[2] = im*sqm[abs(im)]
 
102
            fo[3] = ip*nsf*sqm[abs(im)]
 
103
            fo[4] = im*nsf*sqm[abs(ip)]
 
104
            fo[5] = ip*sqm[abs(ip)]
 
105
 
 
106
        else:
 
107
            sf = [(1+nsf+(1-nsf)*nh)*0.5,(1+nsf-(1-nsf)*nh)*0.5]
 
108
            omega = [sqrt(p[0]+pp),fmass/(sqrt(p[0]+pp))]
 
109
            ip = (1+nh)//2
 
110
            im = (1-nh)//2
 
111
            sfomeg = [sf[0]*omega[ip],sf[1]*omega[im]]
 
112
            pp3 = max(pp+p[3],0.)
 
113
            if (pp3 == 0.):
 
114
                chi1 = complex(-nh,0.) 
 
115
            else:
 
116
                chi1 = complex(nh*p[1]/sqrt(2.*pp*pp3),\
 
117
                -p[2]/sqrt(2.*pp*pp3))
 
118
            chi = [complex(sqrt(pp3*0.5/pp)),chi1]
 
119
 
 
120
            fo[2] = sfomeg[1]*chi[im]
 
121
            fo[3] = sfomeg[1]*chi[ip]
 
122
            fo[4] = sfomeg[0]*chi[im]
 
123
            fo[5] = sfomeg[0]*chi[ip] 
 
124
            
 
125
    else: 
 
126
        sqp0p3 = sqrt(max(p[0]+p[3],0.))*nsf
 
127
        if (sqp0p3 == 0.):
 
128
            chi1 = complex(-nhel*sqrt(2.*p[0]),0.)
 
129
        else:
 
130
            chi1 = complex(nh*p[1]/sqp0p3,-p[2]/sqp0p3)
 
131
        chi = [complex(sqp0p3,0.),chi1]
 
132
        if (nh == 1):
 
133
            fo[2] = chi[0]
 
134
            fo[3] = chi[1]
 
135
            fo[4] = complex(0.,0.)
 
136
            fo[5] = complex(0.,0.) 
 
137
        else:
 
138
            fo[2] = complex(0.,0.)
 
139
            fo[3] = complex(0.,0.)
 
140
            fo[4] = chi[1]
 
141
            fo[5] = chi[0] 
 
142
    
 
143
    return fo
 
144
 
 
145
def vxxxxx(p,vmass,nhel,nsv):
 
146
    """ initialize a vector wavefunction. nhel=4 is for checking BRST"""
 
147
    
 
148
    vc = WaveFunction(3)
 
149
    
 
150
    sqh = sqrt(0.5)
 
151
    nsvahl = nsv*abs(nhel)
 
152
    pt2 = p[1]**2 + p[2]**2 
 
153
    pp = min(p[0],sqrt(pt2 + p[3]**2))
 
154
    pt = min(pp,sqrt(pt2))
 
155
 
 
156
    vc[0] = complex(p[0]*nsv,p[3]*nsv)
 
157
    vc[1] = complex(p[1]*nsv,p[2]*nsv)
 
158
 
 
159
    if (nhel == 4):
 
160
        if (vmass == 0.):
 
161
            vc[2] = 1.
 
162
            vc[3]=p[1]/p[0]
 
163
            vc[4]=p[2]/p[0]
 
164
            vc[5]=p[3]/p[0]
 
165
        else:
 
166
            vc[2] = p[0]/vmass
 
167
            vc[3] = p[1]/vmass
 
168
            vc[4] = p[2]/vmass
 
169
            vc[5] = p[3]/vmass
 
170
        
 
171
        return vc 
 
172
 
 
173
    if (vmass != 0.):
 
174
        hel0 = 1.-abs(nhel) 
 
175
 
 
176
        if (pp == 0.):
 
177
            vc[2] = complex(0.,0.)
 
178
            vc[3] = complex(-nhel*sqh,0.)
 
179
            vc[4] = complex(0.,nsvahl*sqh) 
 
180
            vc[5] = complex(hel0,0.)
 
181
 
 
182
        else:
 
183
            emp = p[0]/(vmass*pp)
 
184
            vc[2] = complex(hel0*pp/vmass,0.)
 
185
            vc[5] = complex(hel0*p[3]*emp+nhel*pt/pp*sqh)
 
186
            if (pt != 0.):
 
187
                pzpt = p[3]/(pp*pt)*sqh*nhel
 
188
                vc[3] = complex(hel0*p[1]*emp-p[1]*pzpt, \
 
189
                    -nsvahl*p[2]/pt*sqh)
 
190
                vc[4] = complex(hel0*p[2]*emp-p[2]*pzpt, \
 
191
                    nsvahl*p[1]/pt*sqh) 
 
192
            else:
 
193
                vc[3] = complex(-nhel*sqh,0.)
 
194
                vc[4] = complex(0.,nsvahl*sign(sqh,p[3]))
 
195
    else: 
 
196
        pp = p[0]
 
197
        pt = sqrt(p[1]**2 + p[2]**2)
 
198
        vc[2] = complex(0.,0.)
 
199
        vc[5] = complex(nhel*pt/pp*sqh)
 
200
        if (pt != 0.):
 
201
            pzpt = p[3]/(pp*pt)*sqh*nhel
 
202
            vc[3] = complex(-p[1]*pzpt,-nsv*p[2]/pt*sqh)
 
203
            vc[4] = complex(-p[2]*pzpt,nsv*p[1]/pt*sqh)
 
204
        else:
 
205
            vc[3] = complex(-nhel*sqh,0.)
 
206
            vc[4] = complex(0.,nsv*sign(sqh,p[3]))
 
207
    
 
208
    return vc
 
209
 
 
210
def sign(x,y):
 
211
    """Fortran's sign transfer function"""
 
212
    try:
 
213
        cmp = (y < 0.)
 
214
    except TypeError:
 
215
        # y should be complex
 
216
        if abs(y.imag) < 1e-6 * abs(y.real):
 
217
            y = y.real
 
218
        else:
 
219
            raise
 
220
    finally:
 
221
        if (y < 0.):
 
222
            return -abs(x) 
 
223
        else:
 
224
            return abs(x) 
 
225
 
 
226
def sxxxxx(p,nss):
 
227
    """initialize a scalar wavefunction"""
 
228
    
 
229
    sc = WaveFunction(1)
 
230
    
 
231
    sc[2] = complex(1.,0.)
 
232
    sc[0] = complex(p[0]*nss,p[3]*nss)
 
233
    sc[1] = complex(p[1]*nss,p[2]*nss)
 
234
    return sc
 
235
 
 
236
 
 
237
def txxxxx(p, tmass, nhel, nst):
 
238
    """ initialize a tensor wavefunction"""
 
239
    
 
240
    tc = WaveFunction(5)
 
241
    
 
242
    sqh = sqrt(0.5)
 
243
    sqs = sqrt(1/6)
 
244
 
 
245
    pt2 = p[1]**2 + p[2]**2
 
246
    pp = min(p[0],sqrt(pt2+p[3]**2))
 
247
    pt = min(pp,sqrt(pt2))
 
248
 
 
249
    ft = {}
 
250
    ft[(4,0)] = complex(p[0], p[3]) * nst
 
251
    ft[(5,0)] = complex(p[1], p[2]) * nst
 
252
 
 
253
    if ( nhel >= 0 ): 
 
254
        #construct eps+
 
255
        ep = [0] * 4
 
256
        
 
257
        if ( pp == 0 ):
 
258
            #ep[0] = 0
 
259
            ep[1] = -sqh
 
260
            ep[2] = complex(0, nst*sqh)
 
261
            #ep[3] = 0
 
262
        else:            
 
263
            #ep[0] = 0
 
264
            ep[3] = pt/pp*sqh
 
265
            if (pt != 0):
 
266
               pzpt = p[3]/(pp*pt)*sqh
 
267
               ep[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh )
 
268
               ep[2] = complex( -p[2]*pzpt ,  nst*p[1]/pt*sqh )
 
269
            else:
 
270
               ep[1] = -sqh 
 
271
               ep[2] = complex( 0 , nst*sign(sqh,p[3]) )
 
272
            
 
273
         
 
274
     
 
275
    if ( nhel <= 0 ): 
 
276
        #construct eps-
 
277
        em = [0] * 4
 
278
        if ( pp == 0 ):
 
279
            #em[0] = 0
 
280
            em[1] = sqh 
 
281
            em[2] = complex( 0 , nst*sqh )
 
282
            #em[3] = 0
 
283
        else:
 
284
            #em[0] = 0
 
285
            em[3] = -pt/pp*sqh
 
286
            if pt:
 
287
               pzpt = -p[3]/(pp*pt)*sqh
 
288
               em[1] = complex( -p[1]*pzpt , -nst*p[2]/pt*sqh )
 
289
               em[2] = complex( -p[2]*pzpt ,  nst*p[1]/pt*sqh )
 
290
            else:
 
291
               em[1] = sqh
 
292
               em[2] = complex( 0 , nst*sign(sqh,p[3]) )
 
293
            
 
294
    
 
295
    if ( abs(nhel) <= 1 ):  
 
296
        #construct eps0
 
297
        e0 = [0] * 4
 
298
        if ( pp == 0 ):
 
299
            #e0[0] = complex( rZero )
 
300
            #e0[1] = complex( rZero )
 
301
            #e0[2] = complex( rZero )
 
302
            e0[3] = 1
 
303
        else:
 
304
            emp = p[0]/(tmass*pp)
 
305
            e0[0] = pp/tmass 
 
306
            e0[3] = p[3]*emp
 
307
            if pt:
 
308
               e0[1] = p[1]*emp 
 
309
               e0[2] = p[2]*emp 
 
310
            #else:
 
311
            #   e0[1] = complex( rZero )
 
312
            #   e0[2] = complex( rZero )
 
313
 
 
314
    if nhel == 2:
 
315
        for j in range(4):
 
316
            for i in range(4):         
 
317
                ft[(i,j)] = ep[i]*ep[j]
 
318
    elif nhel == -2:
 
319
        for j in range(4):
 
320
            for i in range(4):         
 
321
                ft[(i,j)] = em[i]*em[j]
 
322
    elif tmass == 0:
 
323
        for j in range(4):
 
324
            for i in range(4):         
 
325
                ft[(i,j)] = 0
 
326
    elif nhel == 1:
 
327
        for j in range(4):
 
328
            for i in range(4): 
 
329
                ft[(i,j)] = sqh*( ep[i]*e0[j] + e0[i]*ep[j] )
 
330
    elif nhel == 0:
 
331
        for j in range(4):
 
332
            for i in range(4):       
 
333
                ft[(i,j)] = sqs*( ep[i]*em[j] + em[i]*ep[j] + 2 *e0[i]*e0[j] )
 
334
    elif nhel == -1:
 
335
        for j in range(4):
 
336
            for i in range(4): 
 
337
                ft[(i,j)] = sqh*( em[i]*e0[j] + e0[i]*em[j] )
 
338
 
 
339
    else:
 
340
        raise Exception('invalid helicity TXXXXXX') 
 
341
 
 
342
 
 
343
 
 
344
    tc[2] = ft[(0,0)]
 
345
    tc[3] = ft[(0,1)]
 
346
    tc[4] = ft[(0,2)]
 
347
    tc[5] = ft[(0,3)]
 
348
    tc[6] = ft[(1,0)]
 
349
    tc[7] = ft[(1,1)]
 
350
    tc[8] = ft[(1,2)]
 
351
    tc[9] = ft[(1,3)]
 
352
    tc[10] = ft[(2,0)]
 
353
    tc[11] = ft[(2,1)]
 
354
    tc[12] = ft[(2,2)]
 
355
    tc[13] = ft[(2,3)]
 
356
    tc[14] = ft[(3,0)]
 
357
    tc[15] = ft[(3,1)]
 
358
    tc[16] = ft[(3,2)]
 
359
    tc[17] = ft[(3,3)]
 
360
    tc[0] = ft[(4,0)]
 
361
    tc[1] = ft[(5,0)]
 
362
 
 
363
    return tc
 
364
 
 
365
def irxxxx(p, mass, nhel, nsr):
 
366
    """ initialize a incoming spin 3/2 wavefunction."""
 
367
    
 
368
    # This routines is a python translation of a routine written by
 
369
    # K.Mawatari in fortran dated from the 2008/02/26
 
370
    
 
371
    ri = WaveFunction(4)
 
372
    
 
373
    sqh = sqrt(0.5)
 
374
    sq2 = sqrt(2)
 
375
    sq3 = sqrt(3)
 
376
    
 
377
    pt2 = p[1]**2 + p[2]**2
 
378
    pp = min([p[0], sqrt(pt2+p[3]**2)])
 
379
    pt = min([pp,sqrt(pt2)])
 
380
    
 
381
    rc = {}
 
382
    rc[(4,0)] = -1*complex(p[0],p[3])*nsr
 
383
    rc[(5,0)] = -1*complex(p[1],p[2])*nsr
 
384
    
 
385
 
 
386
    nsv = -nsr # nsv=+1 for final, -1 for initial   
 
387
        
 
388
    # Construct eps+
 
389
    if nhel > 0:
 
390
        ep = [0] * 4
 
391
        if pp:
 
392
            #ep[0] = 0
 
393
            ep[3] = pt/pp*sqh
 
394
            if pt:
 
395
                pzpt = p[3]/(pp*pt)*sqh
 
396
                ep[1] = complex( -p[1]*pzpt , -nsv*p[2]/pt*sqh )
 
397
                ep[2] = complex( -p[2]*pzpt ,  nsv*p[1]/pt*sqh )
 
398
            else:
 
399
                ep[1] = -sqh 
 
400
                ep[2] = complex( 0 , nsv*sign(sqh,p[3]) )
 
401
        else:
 
402
            #ep[0] = 0d0
 
403
            ep[1] = -sqh
 
404
            ep[2] = complex(0, nsv * sqh)
 
405
            #ep[3] = 0d0
 
406
 
 
407
         
 
408
    if ( nhel < 0 ): 
 
409
        #construct eps-
 
410
        em = [0] * 4
 
411
        if ( pp == 0 ):
 
412
            #em[0] = 0
 
413
            em[1] = sqh 
 
414
            em[2] = complex( 0 , nsv*sqh )
 
415
            #em[3] = 0
 
416
        else:
 
417
            #em[0] = 0
 
418
            em[3] = -pt/pp*sqh
 
419
            if pt:
 
420
                pzpt = -p[3]/(pp*pt)*sqh
 
421
                em[1] = complex( -p[1]*pzpt , -nsv*p[2]/pt*sqh )
 
422
                em[2] = complex( -p[2]*pzpt ,  nsv*p[1]/pt*sqh )
 
423
            else:
 
424
                em[1] = sqh
 
425
                em[2] = complex( 0 , nsv*sign(sqh,p[3]) )            
 
426
                
 
427
    if ( abs(nhel) <= 1 ):  
 
428
        #construct eps0
 
429
        e0 = [0] * 4
 
430
        if ( pp == 0 ):
 
431
            #e0[0] = complex( rZero )
 
432
            #e0[1] = complex( rZero )
 
433
            #e0[2] = complex( rZero )
 
434
            e0[3] = 1
 
435
        else:
 
436
            emp = p[0]/(mass*pp)
 
437
            e0[0] = pp/mass 
 
438
            e0[3] = p[3]*emp
 
439
            if pt:
 
440
               e0[1] = p[1]*emp 
 
441
               e0[2] = p[2]*emp 
 
442
            #else:
 
443
            #   e0[1] = complex( rZero )
 
444
            #   e0[2] = complex( rZero )
 
445
 
 
446
    
 
447
 
 
448
    if ( nhel >= -1 ):
 
449
        # constract spinor+ 
 
450
        fip = [0] * 4
 
451
        sf, omega, sfomeg, chi = [0, 0], [0,0], [0,0], [0,0]
 
452
        nh = nsr
 
453
        if  mass:
 
454
            pp = min([p[0],sqrt(p[1]**2+p[2]**2+p[3]**2)])
 
455
            if pp == 0:
 
456
                sqm = sqrt(mass)
 
457
                ip = (1+nh)//2
 
458
                im = (1-nh)//2
 
459
                fip[0] = ip * sqm
 
460
                fip[1] = im * nsr * sqm
 
461
                fip[2] = ip * nsr * sqm
 
462
                fip[3] = im * sqm
 
463
            else:
 
464
                sf[0] = float(1+nsr+(1-nsr)*nh)*0.5
 
465
                sf[1] = float(1+nsr-(1-nsr)*nh)*0.5
 
466
                omega[0] = sqrt(p[0]+pp)
 
467
                omega[1] = mass/omega[0]
 
468
                ip = ((3+nh)//2) -1 # -1 since they are index 
 
469
                im = ((3-nh)//2) -1 # -1 since they are index
 
470
                sfomeg[0] = sf[0]*omega[ip]
 
471
                sfomeg[1] = sf[1]*omega[im]
 
472
                pp3 = max([pp+p[3],0])
 
473
                chi[0] = sqrt(pp3*0.5/pp)
 
474
                if  pp3 ==0:
 
475
                    chi[1] = -nh
 
476
                else:
 
477
                    chi[1] = complex( nh*p[1] , p[2] )/sqrt(2*pp*pp3)
 
478
            
 
479
                fip[0] = sfomeg[0]*chi[im]
 
480
                fip[1] = sfomeg[0]*chi[ip]
 
481
                fip[2] = sfomeg[1]*chi[im]
 
482
                fip[3] = sfomeg[1]*chi[ip]
 
483
        else:
 
484
            sqp0p3 = sqrt(max([p[0]+p[3],0])) * nsr
 
485
            chi[0] = sqp0p3
 
486
            if  sqp0p3 == 0:
 
487
                chi[1] = -nhel *  sqrt(2*p[0])
 
488
            else:
 
489
                chi[1] = complex( nh*p[1], p[2] )/sqp0p3
 
490
            if  nh == 1:
 
491
                #fip[0] = complex( rZero )
 
492
                #fip[1] = complex( rZero )
 
493
                fip[2] = chi[0]
 
494
                fip[3] = chi[1]
 
495
            else:
 
496
                fip[0] = chi[1]
 
497
                fip[1] = chi[0]
 
498
                #fip(3) = complex( rZero )
 
499
                #fip(4) = complex( rZero )
 
500
            
 
501
    if ( nhel <= 1 ):
 
502
        # constract spinor- 
 
503
        fim = [0] * 4
 
504
        sf, omega, sfomeg, chi = [0, 0], [0,0], [0,0], [0,0]
 
505
        nh = -nsr
 
506
        if  mass:
 
507
            pp = min([p[0],sqrt(p[1]**2+p[2]**2+p[3]**2)])
 
508
            if pp == 0:
 
509
                sqm = sqrt(mass)
 
510
                ip = (1+nh)/2
 
511
                im = (1-nh)/2
 
512
                fim[0] = ip * sqm
 
513
                fim[1] = im * nsr * sqm
 
514
                fim[2] = ip * nsr * sqm
 
515
                fim[3] = im * sqm
 
516
            else:
 
517
                sf[0] = float(1+nsr+(1-nsr)*nh)*0.5
 
518
                sf[1] = float(1+nsr-(1-nsr)*nh)*0.5
 
519
                omega[0] = sqrt(p[0]+pp)
 
520
                omega[1] = mass/omega[0]
 
521
                ip = (3+nh)//2 -1
 
522
                im = (3-nh)//2 -1
 
523
                sfomeg[0] = sf[0]*omega[ip]
 
524
                sfomeg[1] = sf[1]*omega[im]
 
525
                pp3 = max([pp+p[3],0])
 
526
                chi[0] = sqrt(pp3*0.5/pp)
 
527
                if  pp3 ==0:
 
528
                    chi[1] = -nh
 
529
                else:
 
530
                    chi[1] = complex( nh*p[1] , p[2] )/sqrt(2*pp*pp3)
 
531
            
 
532
                fim[0] = sfomeg[0]*chi[im]
 
533
                fim[1] = sfomeg[0]*chi[ip]
 
534
                fim[2] = sfomeg[1]*chi[im]
 
535
                fim[3] = sfomeg[1]*chi[ip]
 
536
        else:
 
537
            sqp0p3 = sqrt(max([p[0]+p[3],0])) * nsr
 
538
            chi[0] = sqp0p3
 
539
            if  sqp0p3 == 0:
 
540
                chi[1] = -nhel *  sqrt(2*p[0])
 
541
            else:
 
542
                chi[1] = complex( nh*p[1], p[2] )/sqp0p3
 
543
            if  nh == 1:
 
544
                #fip[0] = complex( rZero )
 
545
                #fip[1] = complex( rZero )
 
546
                fim[2] = chi[0]
 
547
                fim[3] = chi[1]
 
548
            else:
 
549
                fim[0] = chi[1]
 
550
                fim[1] = chi[0]
 
551
                #fip(3) = complex( rZero )
 
552
                #fip(4) = complex( rZero )        
 
553
      
 
554
    
 
555
 
 
556
    # recurent relation put her for optimization
 
557
    cond1  = (pt == 0 and p[3] >= 0)
 
558
    cond2  = (pt == 0 and p[3] < 0)
 
559
    
 
560
    # spin-3/2 fermion wavefunction
 
561
    if nhel == 3:
 
562
        for i,j in product(list(range(4)), list(range(4))):
 
563
            rc[(i, j)] = ep[i] *fip[j]
 
564
    
 
565
    elif nhel == 1:
 
566
        for i,j in product(list(range(4)), list(range(4))):
 
567
            if cond1:
 
568
                rc[(i,j)] = sq2/sq3*e0[i]*fip[j] +1/sq3*ep[i]*fim[j]
 
569
            elif cond2:
 
570
                rc[(i,j)] = sq2/sq3*e0[i]*fip[j] -1/sq3*ep[i]*fim[j]
 
571
            else:
 
572
                rc[(i,j)] = sq2/sq3*e0[i]*fip[j] + \
 
573
                                   1/sq3*ep[i]*fim[j] *complex(p[1],nsr*p[2])/pt  
 
574
    elif nhel == -1:
 
575
        for i,j in product(list(range(4)), list(range(4))):
 
576
            if cond1:
 
577
                rc[(i,j)] = 1/sq3*em[i]*fip[j] +sq2/sq3*e0[i]*fim[j]
 
578
            elif cond2:
 
579
                rc[(i,j)] = 1/sq3*em[i]*fip[j] -sq2/sq3*e0[i]*fim[j]
 
580
            else:
 
581
                rc[(i,j)] = 1/sq3*em[i]*fip[j] + \
 
582
                                sq2/sq3*e0[i]*fim[j] *complex(p[1],nsr*p[2])/pt  
 
583
    else:
 
584
        for i,j in product(list(range(4)), list(range(4))):
 
585
            if cond1:
 
586
                rc[(i, j)] = em[i] *fim[j]
 
587
            elif cond2:
 
588
                rc[(i, j)] = -em[i] *fim[j]
 
589
            else:
 
590
                rc[(i, j)] = em[i]*fim[j] *complex(p[1],nsr*p[2])/pt 
 
591
                
 
592
    ri[2] = rc[(0,0)]
 
593
    ri[3] = rc[(0,1)]
 
594
    ri[4] = rc[(0,2)]
 
595
    ri[5] = rc[(0,3)]
 
596
    ri[6] = rc[(1,0)]
 
597
    ri[7] = rc[(1,1)]
 
598
    ri[8] = rc[(1,2)]
 
599
    ri[9] = rc[(1,3)]
 
600
    ri[10] = rc[(2,0)]
 
601
    ri[11] = rc[(2,1)]
 
602
    ri[12] = rc[(2,2)]
 
603
    ri[13] = rc[(2,3)]
 
604
    ri[14] = rc[(3,0)]
 
605
    ri[15] = rc[(3,1)]
 
606
    ri[16] = rc[(3,2)]
 
607
    ri[17] = rc[(3,3)]
 
608
    ri[0] = rc[(4,0)]
 
609
    ri[1] = rc[(5,0)]              
 
610
 
 
611
    return ri
 
612
 
 
613
def orxxxx(p, mass, nhel, nsr):
 
614
    """ initialize a incoming spin 3/2 wavefunction."""
 
615
    
 
616
    # This routines is a python translation of a routine written by
 
617
    # K.Mawatari in fortran dated from the 2008/02/26
 
618
 
 
619
   
 
620
    ro = WaveFunction(spin=4)
 
621
    
 
622
    sqh = sqrt(0.5)
 
623
    sq2 = sqrt(2)
 
624
    sq3 = sqrt(3)
 
625
    
 
626
    pt2 = p[1]**2 + p[2]**2
 
627
    pp = min([p[0], sqrt(pt2+p[3]**2)])
 
628
    pt = min([pp,sqrt(pt2)])
 
629
    rc = {}
 
630
    rc[(4,0)] = complex(p[0],p[3])*nsr
 
631
    rc[(5,0)] = complex(p[1],p[2])*nsr
 
632
    
 
633
 
 
634
    nsv = nsr # nsv=+1 for final, -1 for initial   
 
635
        
 
636
    # Construct eps+
 
637
    if nhel > 0:
 
638
        ep = [0] * 4
 
639
        if pp:
 
640
            #ep[0] = 0
 
641
            ep[3] = pt/pp*sqh
 
642
            if pt:
 
643
                pzpt = p[3]/(pp*pt)*sqh
 
644
                ep[1] = complex( -p[1]*pzpt , -nsv*p[2]/pt*sqh )
 
645
                ep[2] = complex( -p[2]*pzpt ,  nsv*p[1]/pt*sqh )
 
646
            else:
 
647
                ep[1] = -sqh 
 
648
                ep[2] = complex( 0 , nsv*sign(sqh,p[3]) )
 
649
        else:
 
650
            #ep[0] = 0d0
 
651
            ep[1] = -sqh
 
652
            ep[2] = complex(0, nsv * sqh)
 
653
            #ep[3] = 0d0
 
654
         
 
655
    if ( nhel < 0 ): 
 
656
        #construct eps-
 
657
        em = [0] * 4
 
658
        if ( pp == 0 ):
 
659
            #em[0] = 0
 
660
            em[1] = sqh 
 
661
            em[2] = complex( 0 , nsv*sqh )
 
662
            #em[3] = 0
 
663
        else:
 
664
            #em[0] = 0
 
665
            em[3] = -pt/pp*sqh
 
666
            if pt:
 
667
                pzpt = -p[3]/(pp*pt)*sqh
 
668
                em[1] = complex( -p[1]*pzpt , -nsv*p[2]/pt*sqh )
 
669
                em[2] = complex( -p[2]*pzpt ,  nsv*p[1]/pt*sqh )
 
670
            else:
 
671
                em[1] = sqh
 
672
                em[2] = complex( 0 , nsv*sign(sqh,p[3]) )            
 
673
                
 
674
    if ( abs(nhel) <= 1 ):  
 
675
        #construct eps0
 
676
        e0 = [0] * 4
 
677
        if ( pp == 0 ):
 
678
            #e0[0] = complex( rZero )
 
679
            #e0[1] = complex( rZero )
 
680
            #e0[2] = complex( rZero )
 
681
            e0[3] = 1
 
682
        else:
 
683
            emp = p[0]/(mass*pp)
 
684
            e0[0] = pp/mass 
 
685
            e0[3] = p[3]*emp
 
686
            if pt:
 
687
               e0[1] = p[1]*emp 
 
688
               e0[2] = p[2]*emp 
 
689
            #else:
 
690
            #   e0[1] = complex( rZero )
 
691
            #   e0[2] = complex( rZero )
 
692
 
 
693
    if nhel >= -1:
 
694
        #constract spinor+ 
 
695
        nh = nsr
 
696
        sqm, fop, omega, sf, sfomeg = [0]*2,[0]*4,[0]*2,[0]*2,[0]*2
 
697
        chi = [0]*2
 
698
        if mass:
 
699
            pp = min([p[0],sqrt(p[1]**2+p[2]**2+p[3]**2)])
 
700
            if ( pp == 0):
 
701
                sqm[0] = sqrt(abs(mass)) # possibility of negative fermion masses
 
702
                sqm[1] = sign(sqm[0],mass) # possibility of negative fermion masses
 
703
                ip = -((1+nh)/2)
 
704
                im =  (1-nh)/2
 
705
                fop[0] = im     * sqm[im]
 
706
                fop[1] = ip*nsr * sqm[im]
 
707
                fop[2] = im*nsr * sqm[-ip]
 
708
                fop[3] = ip     * sqm[-ip]
 
709
            else:
 
710
                pp = min(p[0],sqrt(p[1]**2+p[2]**2+p[3]**2))
 
711
                sf[0] = (1+nsr+(1-nsr)*nh)*0.5
 
712
                sf[1] = (1+nsr-(1-nsr)*nh)*0.5
 
713
                omega[0] = sqrt(p[0]+pp)
 
714
                omega[1] = mass/omega[0]
 
715
                ip = (3+nh)//2  -1 # -1 since this is index
 
716
                im = (3-nh)//2  -1 # -1 since this is index 
 
717
                sfomeg[0] = sf[0]*omega[ip]
 
718
                sfomeg[1] = sf[1]*omega[im]
 
719
                pp3 = max([pp+p[3],0])
 
720
                chi[0] = sqrt(pp3*0.5/pp)
 
721
                if pp3 == 0:
 
722
                    chi[1] = -nh 
 
723
                else:
 
724
                    chi[1] = complex( nh*p[1] , -p[2] )/sqrt(2*pp*pp3)
 
725
 
 
726
            
 
727
                fop[0] = sfomeg[1]*chi[im]
 
728
                fop[1] = sfomeg[1]*chi[ip]
 
729
                fop[2] = sfomeg[0]*chi[im]
 
730
                fop[3] = sfomeg[0]*chi[ip]
 
731
 
 
732
        else:
 
733
            if(p[1] == 0 and p[2] == 0 and p[3] < 0):
 
734
                sqp0p3 = 0
 
735
            else:
 
736
                sqp0p3 = sqrt(max(p[0]+p[3], 0))*nsr
 
737
                
 
738
            chi[0] =  sqp0p3
 
739
            if ( sqp0p3 == 0 ):
 
740
                chi[1] = complex(-nhel )*sqrt(2*p[0])
 
741
            else:
 
742
                chi[1] = complex( nh*p[1], -p[2] )/sqp0p3
 
743
         
 
744
            if ( nh == 1 ):
 
745
                fop[0] = chi[0]
 
746
                fop[1] = chi[1]
 
747
                #fop[2] = 0
 
748
                #fop[3] = 0
 
749
            else:
 
750
                #fop[0] = 0
 
751
                #fop[1] = 0
 
752
                fop[2] = chi[1]
 
753
                fop[3] = chi[0]
 
754
         
 
755
    
 
756
    if ( nhel < 2 ):
 
757
        # constract spinor+ 
 
758
        sqm, fom, omega, sf, sfomeg = [0]*2,[0]*4,[0]*2,[0]*2,[0]*2
 
759
        chi = [0]*2
 
760
 
 
761
        
 
762
        nh = -nsr
 
763
        if mass:
 
764
            pp = min([p[0],sqrt(p[1]**2+p[2]**2+p[3]**2)])
 
765
            if ( pp == 0):
 
766
                sqm[0] = sqrt(abs(mass)) # possibility of negative fermion masses
 
767
                sqm[1] = sign(sqm[0],mass) # possibility of negative fermion masses
 
768
                ip = -((1+nh)/2)
 
769
                im =  (1-nh)/2
 
770
            
 
771
                fom[0] = im     * sqm[im]
 
772
                fom[1] = ip*nsr * sqm[im]
 
773
                fom[2] = im*nsr * sqm[-ip]
 
774
                fom[3] = ip     * sqm[-ip]
 
775
            
 
776
            else:
 
777
                pp = min([p[0],sqrt(p[1]**2+p[2]**2+p[3]**2)])
 
778
                sf[0] = (1+nsr+(1-nsr)*nh)*0.5
 
779
                sf[1] = (1+nsr-(1-nsr)*nh)*0.5
 
780
                omega[0] = sqrt(p[0]+pp)
 
781
                omega[1] = mass/omega[0]
 
782
                ip = (3+nh)//2 -1 #-1 since ip is an index
 
783
                im = (3-nh)//2 -1 
 
784
                sfomeg[0] = sf[0]*omega[ip]
 
785
                sfomeg[1] = sf[1]*omega[im]
 
786
                pp3 = max([pp+p[3], 0])
 
787
                chi[0] = sqrt(pp3*0.5/pp)
 
788
                if ( pp3 == 0):
 
789
                    chi[1] = -nh
 
790
                else:
 
791
                    chi[1] = complex( nh*p[1] , -p[2] )/sqrt(2*pp*pp3)
 
792
            
 
793
            
 
794
                fom[0] = sfomeg[1]*chi[im]
 
795
                fom[1] = sfomeg[1]*chi[ip]
 
796
                fom[2] = sfomeg[0]*chi[im]
 
797
                fom[3] = sfomeg[0]*chi[ip]
 
798
        else:
 
799
            if(p[1] == 0 == p[2] and p[3] < 0):
 
800
                sqp0p3 = 0
 
801
            else:
 
802
                sqp0p3 = sqrt(max([p[0]+p[3],0]))*nsr
 
803
            chi[0] = sqp0p3
 
804
            if ( sqp0p3 == 0):
 
805
                chi[1] = complex(-nhel )*sqrt(2*p[0])
 
806
            else:
 
807
                chi[1] = complex( nh*p[1], -p[2] )/sqp0p3
 
808
            if ( nh == 1 ):
 
809
                fom[0] = chi[0]
 
810
                fom[1] = chi[1]
 
811
                #fom[2] = 0
 
812
                #fom[3] = 0
 
813
            else:
 
814
                #fom[1] = 0
 
815
                #fom[2] = 0
 
816
                fom[2] = chi[1]
 
817
                fom[3] = chi[0]
 
818
 
 
819
    cond1 = ( pt==0 and p[3]>=0)
 
820
    cond2= (pt==0 and p[3]<0)
 
821
 
 
822
   
 
823
    # spin-3/2 fermion wavefunction
 
824
    if nhel == 3:
 
825
        for i,j in product(list(range(4)), list(range(4))):
 
826
            rc[(i, j)] = ep[i] *fop[j]  
 
827
    
 
828
 
 
829
    elif nhel == 1:
 
830
        for i,j in product(list(range(4)), list(range(4))):
 
831
            if cond1:
 
832
                rc[(i,j)] = sq2/sq3*e0[i]*fop[j] + 1/sq3*ep[i]*fom[j]
 
833
            elif cond2:
 
834
                rc[(i,j)] = sq2/sq3*e0[i]*fop[j] - 1/sq3*ep[i]*fom[j]
 
835
            else:
 
836
                rc[(i,j)] = sq2/sq3*e0[i]*fop[j] + 1/sq3*ep[i]*fom[j] * \
 
837
                                                      complex(p[1],-nsr*p[2])/pt  
 
838
                
 
839
    elif nhel == -1:
 
840
        for i,j in product(list(range(4)), list(range(4))):
 
841
            if cond1:
 
842
                rc[(i,j)] = 1/sq3*em[i]*fop[j]+sq2/sq3*e0[i]*fom[j]
 
843
            elif cond2:
 
844
                rc[(i,j)] =1/sq3*em[i]*fop[j]-sq2/sq3*e0[i]*fom[j]
 
845
            else:
 
846
                rc[(i,j)] =  1/sq3*em[i]*fop[j] + sq2/sq3*e0[i]*fom[j] *\
 
847
                                                      complex(p[1],-nsr*p[2])/pt              
 
848
    else:
 
849
        for i,j in product(list(range(4)), list(range(4))):
 
850
            if cond1:
 
851
                rc[(i,j)] = em[i] * fom[j]
 
852
            elif cond2:
 
853
                rc[(i,j)] = - em[i] * fom[j]
 
854
            else:
 
855
                rc[(i,j)] = em[i] * fom[j] * complex(p[1],-nsr*p[2])/pt 
 
856
 
 
857
 
 
858
 
 
859
    ro[2] = rc[(0,0)]
 
860
    ro[3] = rc[(0,1)]
 
861
    ro[4] = rc[(0,2)]
 
862
    ro[5] = rc[(0,3)]
 
863
    ro[6] = rc[(1,0)]
 
864
    ro[7] = rc[(1,1)]
 
865
    ro[8] = rc[(1,2)]
 
866
    ro[9] = rc[(1,3)]
 
867
    ro[10] = rc[(2,0)]
 
868
    ro[11] = rc[(2,1)]
 
869
    ro[12] = rc[(2,2)]
 
870
    ro[13] = rc[(2,3)]
 
871
    ro[14] = rc[(3,0)]
 
872
    ro[15] = rc[(3,1)]
 
873
    ro[16] = rc[(3,2)]
 
874
    ro[17] = rc[(3,3)]
 
875
    ro[0] = rc[(4,0)]
 
876
    ro[1] = rc[(5,0)]
 
877
    
 
878
    return ro
 
879
    
 
880
 
 
881