1
C###############################################################################
3
C Copyright (c) 2010 The ALOHA Development team and Contributors
5
C This file is a part of the MadGraph5_aMC@NLO project, an application which
6
C automatically generates Feynman diagrams and matrix elements for arbitrary
7
C high-energy processes in the Standard Model and beyond.
9
C It is subject to the ALOHA license which should accompany this
12
C###############################################################################
13
subroutine ixxxxx(p, fmass, nhel, nsf ,fi)
15
c This subroutine computes a fermion wavefunction with the flowing-IN
19
c real p(0:3) : four-momentum of fermion
20
c real fmass : mass of fermion
21
c integer nhel = -1 or 1 : helicity of fermion
22
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
25
c complex fi(6) : fermion wavefunction |fi>
28
double complex fi(8),chi(2)
29
double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
30
& pp,pp3,sqp0p3,sqm(0:1)
31
integer nhel,nsf,ip,im,nh
33
double precision rZero, rHalf, rTwo
34
parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
38
c double precision epsi
39
c parameter( epsi = 2.0d-5 )
41
c parameter( stdo = 6 )
45
c pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
46
c if ( abs(p(0))+pp.eq.rZero ) then
48
c & ' helas-error : p(0:3) in ixxxxx is zero momentum'
50
c if ( p(0).le.rZero ) then
52
c & ' helas-error : p(0:3) in ixxxxx has non-positive energy'
56
c p2 = (p(0)-pp)*(p(0)+pp)
57
c if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
59
c & ' helas-error : p(0:3) in ixxxxx has inappropriate mass'
61
c & ' : p**2 = ',p2,' : fmass**2 = ',fmass**2
63
c if (abs(nhel).ne.1) then
64
c write(stdo,*) ' helas-error : nhel in ixxxxx is not -1,1'
65
c write(stdo,*) ' : nhel = ',nhel
67
c if (abs(nsf).ne.1) then
68
c write(stdo,*) ' helas-error : nsf in ixxxxx is not -1,1'
69
c write(stdo,*) ' : nsf = ',nsf
73
c Convention for trees
74
c fi(5) = dcmplx(p(0),p(3))*nsf
75
c fi(6) = dcmplx(p(1),p(2))*nsf
77
c Convention for loop computations
78
fi(1) = dcmplx(p(0),0.D0)*(-nsf)
79
fi(2) = dcmplx(p(1),0.D0)*(-nsf)
80
fi(3) = dcmplx(p(2),0.D0)*(-nsf)
81
fi(4) = dcmplx(p(3),0.D0)*(-nsf)
85
if ( fmass.ne.rZero ) then
87
pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
89
if ( pp.eq.rZero ) then
91
sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
92
sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
97
fi(6) = im*nsf * sqm(ip)
98
fi(7) = ip*nsf * sqm(im)
103
sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
104
sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
105
omega(1) = dsqrt(p(0)+pp)
106
omega(2) = fmass/omega(1)
109
sfomeg(1) = sf(1)*omega(ip)
110
sfomeg(2) = sf(2)*omega(im)
111
pp3 = max(pp+p(3),rZero)
112
chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
113
if ( pp3.eq.rZero ) then
114
chi(2) = dcmplx(-nh )
116
chi(2) = dcmplx( nh*p(1) , p(2) )/dsqrt(rTwo*pp*pp3)
119
fi(5) = sfomeg(1)*chi(im)
120
fi(6) = sfomeg(1)*chi(ip)
121
fi(7) = sfomeg(2)*chi(im)
122
fi(8) = sfomeg(2)*chi(ip)
128
if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
131
sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
133
chi(1) = dcmplx( sqp0p3 )
134
if ( sqp0p3.eq.rZero ) then
135
chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
137
chi(2) = dcmplx( nh*p(1), p(2) )/sqp0p3
140
fi(5) = dcmplx( rZero )
141
fi(6) = dcmplx( rZero )
147
fi(7) = dcmplx( rZero )
148
fi(8) = dcmplx( rZero )
156
subroutine ixxxso(p, fmass, nhel, nsf ,fi)
158
c This subroutine computes a fermion wavefunction with the flowing-IN
162
c real p(0:3) : four-momentum of fermion
163
c real fmass : mass of fermion
164
c integer nhel = -1 or 1 : helicity of fermion
165
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
168
c complex fi(6) : fermion wavefunction |fi>
171
double complex fi(4),chi(2)
172
double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
173
& pp,pp3,sqp0p3,sqm(0:1)
174
integer nhel,nsf,ip,im,nh
176
double precision rZero, rHalf, rTwo
177
parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
180
c double precision p2
181
c double precision epsi
182
c parameter( epsi = 2.0d-5 )
184
c parameter( stdo = 6 )
188
c pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
189
c if ( abs(p(0))+pp.eq.rZero ) then
191
c & ' helas-error : p(0:3) in ixxxxx is zero momentum'
193
c if ( p(0).le.rZero ) then
195
c & ' helas-error : p(0:3) in ixxxxx has non-positive energy'
197
c & ' : p(0) = ',p(0)
199
c p2 = (p(0)-pp)*(p(0)+pp)
200
c if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
202
c & ' helas-error : p(0:3) in ixxxxx has inappropriate mass'
204
c & ' : p**2 = ',p2,' : fmass**2 = ',fmass**2
206
c if (abs(nhel).ne.1) then
207
c write(stdo,*) ' helas-error : nhel in ixxxxx is not -1,1'
208
c write(stdo,*) ' : nhel = ',nhel
210
c if (abs(nsf).ne.1) then
211
c write(stdo,*) ' helas-error : nsf in ixxxxx is not -1,1'
212
c write(stdo,*) ' : nsf = ',nsf
216
c Convention for trees
217
c fi(5) = dcmplx(p(0),p(3))*nsf
218
c fi(6) = dcmplx(p(1),p(2))*nsf
220
c$$$c Convention for loop computations
221
c$$$ fi(1) = dcmplx(p(0),0.D0)*(-nsf)
222
c$$$ fi(2) = dcmplx(p(1),0.D0)*(-nsf)
223
c$$$ fi(3) = dcmplx(p(2),0.D0)*(-nsf)
224
c$$$ fi(4) = dcmplx(p(3),0.D0)*(-nsf)
228
if ( fmass.ne.rZero ) then
230
pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
232
if ( pp.eq.rZero ) then
234
sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
235
sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
240
fi(2) = im*nsf * sqm(ip)
241
fi(3) = ip*nsf * sqm(im)
246
sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
247
sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
248
omega(1) = dsqrt(p(0)+pp)
249
omega(2) = fmass/omega(1)
252
sfomeg(1) = sf(1)*omega(ip)
253
sfomeg(2) = sf(2)*omega(im)
254
pp3 = max(pp+p(3),rZero)
255
chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
256
if ( pp3.eq.rZero ) then
257
chi(2) = dcmplx(-nh )
259
chi(2) = dcmplx( nh*p(1) , p(2) )/dsqrt(rTwo*pp*pp3)
262
fi(1) = sfomeg(1)*chi(im)
263
fi(2) = sfomeg(1)*chi(ip)
264
fi(3) = sfomeg(2)*chi(im)
265
fi(4) = sfomeg(2)*chi(ip)
271
if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
274
sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
276
chi(1) = dcmplx( sqp0p3 )
277
if ( sqp0p3.eq.rZero ) then
278
chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
280
chi(2) = dcmplx( nh*p(1), p(2) )/sqp0p3
283
fi(1) = dcmplx( rZero )
284
fi(2) = dcmplx( rZero )
290
fi(3) = dcmplx( rZero )
291
fi(4) = dcmplx( rZero )
299
subroutine mp_ixxxxx(p, fmass, nhel, nsf ,fi)
301
c This subroutine computes a fermion wavefunction with the flowing-IN
302
c fermion number, in QUADRUPLE PRECISIOn
305
c real p(0:3) : four-momentum of fermion
306
c real fmass : mass of fermion
307
c integer nhel = -1 or 1 : helicity of fermion
308
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
311
c complex fi(6) : fermion wavefunction |fi>
314
complex*32 fi(8),chi(2)
315
real*16 p(0:3),sf(2),sfomeg(2),omega(2),fmass,
316
& pp,pp3,sqp0p3,sqm(0:1)
317
integer nhel,nsf,ip,im,nh
319
real*16 rZero, rHalf, rTwo
320
parameter( rZero = 0.0e0_16, rHalf = 0.5e0_16, rTwo = 2.0e0_16 )
321
c Convention for loop computations
322
fi(1) = cmplx(p(0),rZero,KIND=16)*(-nsf)
323
fi(2) = cmplx(p(1),rZero,KIND=16)*(-nsf)
324
fi(3) = cmplx(p(2),rZero,KIND=16)*(-nsf)
325
fi(4) = cmplx(p(3),rZero,KIND=16)*(-nsf)
329
if ( fmass.ne.rZero ) then
331
pp = min(p(0),sqrt(p(1)**2+p(2)**2+p(3)**2))
333
if ( pp.eq.rZero ) then
335
sqm(0) = sqrt(abs(fmass)) ! possibility of negative fermion masses
336
sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
341
fi(6) = im*nsf * sqm(ip)
342
fi(7) = ip*nsf * sqm(im)
347
sf(1) = REAL(1+nsf+(1-nsf)*nh,KIND=16)*rHalf
348
sf(2) = REAL(1+nsf-(1-nsf)*nh,KIND=16)*rHalf
349
omega(1) = sqrt(p(0)+pp)
350
omega(2) = fmass/omega(1)
353
sfomeg(1) = sf(1)*omega(ip)
354
sfomeg(2) = sf(2)*omega(im)
355
pp3 = max(pp+p(3),rZero)
356
chi(1) = cmplx( sqrt(pp3*rHalf/pp), KIND=16 )
357
if ( pp3.eq.rZero ) then
358
chi(2) = cmplx(-nh ,KIND=16)
360
chi(2) = cmplx( nh*p(1) , p(2),KIND=16)/sqrt(rTwo*pp*pp3)
363
fi(5) = sfomeg(1)*chi(im)
364
fi(6) = sfomeg(1)*chi(ip)
365
fi(7) = sfomeg(2)*chi(im)
366
fi(8) = sfomeg(2)*chi(ip)
372
if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
375
sqp0p3 = sqrt(max(p(0)+p(3),rZero))*nsf
377
chi(1) = cmplx( sqp0p3 ,KIND=16)
378
if ( sqp0p3.eq.rZero ) then
379
chi(2) = cmplx(-nhel ,KIND=16)*sqrt(rTwo*p(0))
381
chi(2) = cmplx( nh*p(1), p(2) ,KIND=16)/sqp0p3
384
fi(5) = cmplx( rZero ,KIND=16)
385
fi(6) = cmplx( rZero ,KIND=16)
391
fi(7) = cmplx( rZero ,KIND=16)
392
fi(8) = cmplx( rZero ,KIND=16)
399
subroutine oxxxxx(p,fmass,nhel,nsf , fo)
401
c This subroutine computes a fermion wavefunction with the flowing-OUT
405
c real p(0:3) : four-momentum of fermion
406
c real fmass : mass of fermion
407
c integer nhel = -1 or 1 : helicity of fermion
408
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
411
c complex fo(6) : fermion wavefunction <fo|
414
double complex fo(8),chi(2)
415
double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
416
& pp,pp3,sqp0p3,sqm(0:1)
417
integer nhel,nsf,nh,ip,im
419
double precision rZero, rHalf, rTwo
420
parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
423
c double precision p2
424
c double precision epsi
425
c parameter( epsi = 2.0d-5 )
427
c parameter( stdo = 6 )
431
c pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
432
c if ( abs(p(0))+pp.eq.rZero ) then
434
c & ' helas-error : p(0:3) in oxxxxx is zero momentum'
436
c if ( p(0).le.rZero ) then
438
c & ' helas-error : p(0:3) in oxxxxx has non-positive energy'
440
c & ' : p(0) = ',p(0)
442
c p2 = (p(0)-pp)*(p(0)+pp)
443
c if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
445
c & ' helas-error : p(0:3) in oxxxxx has inappropriate mass'
447
c & ' : p**2 = ',p2,' : fmass**2 = ',fmass**2
449
c if ( abs(nhel).ne.1 ) then
450
c write(stdo,*) ' helas-error : nhel in oxxxxx is not -1,1'
451
c write(stdo,*) ' : nhel = ',nhel
453
c if ( abs(nsf).ne.1 ) then
454
c write(stdo,*) ' helas-error : nsf in oxxxxx is not -1,1'
455
c write(stdo,*) ' : nsf = ',nsf
459
c Convention for trees
460
c fo(5) = dcmplx(p(0),p(3))*nsf
461
c fo(6) = dcmplx(p(1),p(2))*nsf
463
c Convention for loop computations
464
fo(1) = dcmplx(p(0),0.D0)*(nsf)
465
fo(2) = dcmplx(p(1),0.D0)*(nsf)
466
fo(3) = dcmplx(p(2),0.D0)*(nsf)
467
fo(4) = dcmplx(p(3),0.D0)*(nsf)
472
if ( fmass.ne.rZero ) then
474
pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
476
if ( pp.eq.rZero ) then
478
sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
479
sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
484
fo(6) = ip*nsf * sqm(im)
485
fo(7) = im*nsf * sqm(-ip)
486
fo(8) = ip * sqm(-ip)
490
pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
491
sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
492
sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
493
omega(1) = dsqrt(p(0)+pp)
494
omega(2) = fmass/omega(1)
497
sfomeg(1) = sf(1)*omega(ip)
498
sfomeg(2) = sf(2)*omega(im)
499
pp3 = max(pp+p(3),rZero)
500
chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
501
if ( pp3.eq.rZero ) then
502
chi(2) = dcmplx(-nh )
504
chi(2) = dcmplx( nh*p(1) , -p(2) )/dsqrt(rTwo*pp*pp3)
507
fo(5) = sfomeg(2)*chi(im)
508
fo(6) = sfomeg(2)*chi(ip)
509
fo(7) = sfomeg(1)*chi(im)
510
fo(8) = sfomeg(1)*chi(ip)
516
if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
519
sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
521
chi(1) = dcmplx( sqp0p3 )
522
if ( sqp0p3.eq.rZero ) then
523
chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
525
chi(2) = dcmplx( nh*p(1), -p(2) )/sqp0p3
530
fo(7) = dcmplx( rZero )
531
fo(8) = dcmplx( rZero )
533
fo(5) = dcmplx( rZero )
534
fo(6) = dcmplx( rZero )
544
subroutine oxxxso(p,fmass,nhel,nsf , fo)
546
c This subroutine computes a fermion wavefunction with the flowing-OUT
550
c real p(0:3) : four-momentum of fermion
551
c real fmass : mass of fermion
552
c integer nhel = -1 or 1 : helicity of fermion
553
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
556
c complex fo(6) : fermion wavefunction <fo|
559
double complex fo(4),chi(2)
560
double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
561
& pp,pp3,sqp0p3,sqm(0:1)
562
integer nhel,nsf,nh,ip,im
564
double precision rZero, rHalf, rTwo
565
parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
568
c double precision p2
569
c double precision epsi
570
c parameter( epsi = 2.0d-5 )
572
c parameter( stdo = 6 )
576
c pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
577
c if ( abs(p(0))+pp.eq.rZero ) then
579
c & ' helas-error : p(0:3) in oxxxxx is zero momentum'
581
c if ( p(0).le.rZero ) then
583
c & ' helas-error : p(0:3) in oxxxxx has non-positive energy'
585
c & ' : p(0) = ',p(0)
587
c p2 = (p(0)-pp)*(p(0)+pp)
588
c if ( abs(p2-fmass**2).gt.p(0)**2*epsi ) then
590
c & ' helas-error : p(0:3) in oxxxxx has inappropriate mass'
592
c & ' : p**2 = ',p2,' : fmass**2 = ',fmass**2
594
c if ( abs(nhel).ne.1 ) then
595
c write(stdo,*) ' helas-error : nhel in oxxxxx is not -1,1'
596
c write(stdo,*) ' : nhel = ',nhel
598
c if ( abs(nsf).ne.1 ) then
599
c write(stdo,*) ' helas-error : nsf in oxxxxx is not -1,1'
600
c write(stdo,*) ' : nsf = ',nsf
604
c Convention for trees
605
c fo(5) = dcmplx(p(0),p(3))*nsf
606
c fo(6) = dcmplx(p(1),p(2))*nsf
608
c$$$c Convention for loop computations
609
c$$$ fo(1) = dcmplx(p(0),0.D0)*(nsf)
610
c$$$ fo(2) = dcmplx(p(1),0.D0)*(nsf)
611
c$$$ fo(3) = dcmplx(p(2),0.D0)*(nsf)
612
c$$$ fo(4) = dcmplx(p(3),0.D0)*(nsf)
617
if ( fmass.ne.rZero ) then
619
pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
621
if ( pp.eq.rZero ) then
623
sqm(0) = dsqrt(abs(fmass)) ! possibility of negative fermion masses
624
sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
629
fo(2) = ip*nsf * sqm(im)
630
fo(3) = im*nsf * sqm(-ip)
631
fo(4) = ip * sqm(-ip)
635
pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
636
sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
637
sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
638
omega(1) = dsqrt(p(0)+pp)
639
omega(2) = fmass/omega(1)
642
sfomeg(1) = sf(1)*omega(ip)
643
sfomeg(2) = sf(2)*omega(im)
644
pp3 = max(pp+p(3),rZero)
645
chi(1) = dcmplx( dsqrt(pp3*rHalf/pp) )
646
if ( pp3.eq.rZero ) then
647
chi(2) = dcmplx(-nh )
649
chi(2) = dcmplx( nh*p(1) , -p(2) )/dsqrt(rTwo*pp*pp3)
652
fo(1) = sfomeg(2)*chi(im)
653
fo(2) = sfomeg(2)*chi(ip)
654
fo(3) = sfomeg(1)*chi(im)
655
fo(4) = sfomeg(1)*chi(ip)
661
if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
664
sqp0p3 = dsqrt(max(p(0)+p(3),rZero))*nsf
666
chi(1) = dcmplx( sqp0p3 )
667
if ( sqp0p3.eq.rZero ) then
668
chi(2) = dcmplx(-nhel )*dsqrt(rTwo*p(0))
670
chi(2) = dcmplx( nh*p(1), -p(2) )/sqp0p3
675
fo(3) = dcmplx( rZero )
676
fo(4) = dcmplx( rZero )
678
fo(1) = dcmplx( rZero )
679
fo(2) = dcmplx( rZero )
689
subroutine mp_oxxxxx(p,fmass,nhel,nsf , fo)
691
c This subroutine computes a fermion wavefunction with the flowing-OUT
692
c fermion number in quadruple precision.
695
c real p(0:3) : four-momentum of fermion
696
c real fmass : mass of fermion
697
c integer nhel = -1 or 1 : helicity of fermion
698
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
701
c complex fo(6) : fermion wavefunction <fo|
704
complex*32 fo(8),chi(2)
705
real*16 p(0:3),sf(2),sfomeg(2),omega(2),fmass,
706
& pp,pp3,sqp0p3,sqm(0:1)
707
integer nhel,nsf,nh,ip,im
709
real*16 rZero, rHalf, rTwo
710
parameter( rZero = 0.0e0_16, rHalf = 0.5e0_16, rTwo = 2.0e0_16 )
712
c Convention for loop computations
713
fo(1) = cmplx(p(0),rZero,KIND=16)*nsf
714
fo(2) = cmplx(p(1),rZero,KIND=16)*nsf
715
fo(3) = cmplx(p(2),rZero,KIND=16)*nsf
716
fo(4) = cmplx(p(3),rZero,KIND=16)*nsf
721
if ( fmass.ne.rZero ) then
723
pp = min(p(0),sqrt(p(1)**2+p(2)**2+p(3)**2))
725
if ( pp.eq.rZero ) then
727
sqm(0) = sqrt(abs(fmass)) ! possibility of negative fermion masses
728
sqm(1) = sign(sqm(0),fmass) ! possibility of negative fermion masses
733
fo(6) = ip*nsf * sqm(im)
734
fo(7) = im*nsf * sqm(-ip)
735
fo(8) = ip * sqm(-ip)
739
pp = min(p(0),sqrt(p(1)**2+p(2)**2+p(3)**2))
740
sf(1) = real(1+nsf+(1-nsf)*nh,KIND=16)*rHalf
741
sf(2) = real(1+nsf-(1-nsf)*nh,KIND=16)*rHalf
742
omega(1) = sqrt(p(0)+pp)
743
omega(2) = fmass/omega(1)
746
sfomeg(1) = sf(1)*omega(ip)
747
sfomeg(2) = sf(2)*omega(im)
748
pp3 = max(pp+p(3),rZero)
749
chi(1) = cmplx( sqrt(pp3*rHalf/pp) ,KIND=16)
750
if ( pp3.eq.rZero ) then
751
chi(2) = cmplx(-nh ,KIND=16)
753
chi(2) = cmplx( nh*p(1) , -p(2) , KIND=16)/
757
fo(5) = sfomeg(2)*chi(im)
758
fo(6) = sfomeg(2)*chi(ip)
759
fo(7) = sfomeg(1)*chi(im)
760
fo(8) = sfomeg(1)*chi(ip)
766
if(p(1).eq.0d0.and.p(2).eq.0d0.and.p(3).lt.0d0) then
769
sqp0p3 = sqrt(max(p(0)+p(3),rZero))*nsf
771
chi(1) = cmplx( sqp0p3 ,KIND=16)
772
if ( sqp0p3.eq.rZero ) then
773
chi(2) = cmplx(-nhel , KIND=16)*sqrt(rTwo*p(0))
775
chi(2) = cmplx( nh*p(1), -p(2) , KIND=16)/sqp0p3
780
fo(7) = cmplx( rZero , KIND=16)
781
fo(8) = cmplx( rZero , KIND=16)
783
fo(5) = cmplx( rZero , KIND=16)
784
fo(6) = cmplx( rZero , KIND=16)
794
subroutine pxxxxx(p,tmass,nhel,nst , tc)
798
c This subroutine computes a PSEUDOR wavefunction.
801
c real p(0:3) : four-momentum of tensor boson
802
c real tmass : mass of tensor boson
803
c integer nhel : helicity of tensor boson
804
c = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
805
c integer nst = -1 or 1 : +1 for final, -1 for initial
808
c complex tc(18) : PSEUDOR wavefunction epsilon^mu^nu(t)
811
double precision p(0:3), tmass
813
double complex tc(20)
815
double complex ft(6,4), ep(4), em(4), e0(4)
816
double precision pt, pt2, pp, pzpt, emp, sqh, sqs
819
double precision rZero, rHalf, rOne, rTwo
820
parameter( rZero = 0.0d0, rHalf = 0.5d0 )
821
parameter( rOne = 1.0d0, rTwo = 2.0d0 )
826
c Convention for trees
827
c tc(17) = dcmplx(p(0),p(3))*nst
828
c tc(18) = dcmplx(p(1),p(2))*nst
830
c Convention for loop computations
831
tc(1) = dcmplx(p(0),0.D0)*nst
832
tc(2) = dcmplx(p(1),0.D0)*nst
833
tc(3) = dcmplx(p(2),0.D0)*nst
834
tc(4) = dcmplx(p(3),0.D0)*nst
839
subroutine mp_pxxxxx(p,tmass,nhel,nst , tc)
843
c This subroutine computes a PSEUDOR wavefunction.
846
c real p(0:3) : four-momentum of tensor boson
847
c real tmass : mass of tensor boson
848
c integer nhel : helicity of tensor boson
849
c = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
850
c integer nst = -1 or 1 : +1 for final, -1 for initial
853
c complex tc(18) : PSEUDOR wavefunction epsilon^mu^nu(t)
856
real*16 p(0:3), tmass
860
complex*32 ft(6,4), ep(4), em(4), e0(4)
861
real*16 pt, pt2, pp, pzpt, emp, sqh, sqs
864
real*16 rZero, rHalf, rOne, rTwo
865
parameter( rZero = 0.0e0_16, rHalf = 0.5e0_16 )
866
parameter( rOne = 1.0e0_16, rTwo = 2.0e0_16 )
871
c Convention for trees
872
c tc(17) = dcmplx(p(0),p(3))*nst
873
c tc(18) = dcmplx(p(1),p(2))*nst
875
c Convention for loop computations
876
tc(1) = cmplx(p(0),0.0e0_16,KIND=16)*nst
877
tc(2) = cmplx(p(1),0.0e0_16,KIND=16)*nst
878
tc(3) = cmplx(p(2),0.0e0_16,KIND=16)*nst
879
tc(4) = cmplx(p(3),0.0e0_16,KIND=16)*nst
884
subroutine sxxxxx(p,nss , sc)
886
c This subroutine computes a complex SCALAR wavefunction.
889
c real p(0:3) : four-momentum of scalar boson
890
c integer nss = -1 or 1 : +1 for final, -1 for initial
893
c complex sc(3) : scalar wavefunction s
897
double precision p(0:3)
900
double precision rOne
901
parameter( rOne = 1.0d0 )
904
c double precision p2
905
c double precision epsi
906
c parameter( epsi = 2.0d-5 )
907
c double precision rZero
908
c parameter( rZero = 0.0d0 )
910
c parameter( stdo = 6 )
914
c if ( abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero ) then
916
c & ' helas-error : p(0:3) in sxxxxx is zero momentum'
918
c if ( p(0).le.rZero ) then
920
c & ' helas-error : p(0:3) in sxxxxx has non-positive energy'
922
c & ' : p(0) = ',p(0)
924
c p2 = p(0)**2-(p(1)**2+p(2)**2+p(3)**2)
925
c if ( p2.lt.-p(0)**2*epsi ) then
926
c write(stdo,*) ' helas-error : p(0:3) in sxxxxx is spacelike'
927
c write(stdo,*) ' : p**2 = ',p2
929
c if ( abs(nss).ne.1 ) then
930
c write(stdo,*) ' helas-error : nss in sxxxxx is not -1,1'
931
c write(stdo,*) ' : nss = ',nss
935
sc(5) = dcmplx( rOne )
937
c Convention for trees
938
c sc(2) = dcmplx(p(0),p(3))*nss
939
c sc(3) = dcmplx(p(1),p(2))*nss
941
c Convention for loop computations
942
sc(1) = dcmplx(p(0),0.D0)*nss
943
sc(2) = dcmplx(p(1),0.D0)*nss
944
sc(3) = dcmplx(p(2),0.D0)*nss
945
sc(4) = dcmplx(p(3),0.D0)*nss
950
subroutine mp_sxxxxx(p,nss , sc)
952
c This subroutine computes a complex SCALAR wavefunction.
953
c in quadrupole precision.
956
c real p(0:3) : four-momentum of scalar boson
957
c integer nss = -1 or 1 : +1 for final, -1 for initial
960
c complex sc(3) : scalar wavefunction s
968
parameter( rOne = 1.0e0_16 )
970
sc(5) = cmplx( rOne , KIND=16)
972
c Convention for trees
973
c sc(2) = dcmplx(p(0),p(3))*nss
974
c sc(3) = dcmplx(p(1),p(2))*nss
976
c Convention for loop computations
977
sc(1) = cmplx(p(0),0.0e0_16, KIND=16)*nss
978
sc(2) = cmplx(p(1),0.0e0_16, KIND=16)*nss
979
sc(3) = cmplx(p(2),0.0e0_16, KIND=16)*nss
980
sc(4) = cmplx(p(3),0.0e0_16, KIND=16)*nss
985
subroutine txxxxx(p,tmass,nhel,nst , tc)
987
c This subroutine computes a TENSOR wavefunction.
990
c real p(0:3) : four-momentum of tensor boson
991
c real tmass : mass of tensor boson
992
c integer nhel : helicity of tensor boson
993
c = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
994
c integer nst = -1 or 1 : +1 for final, -1 for initial
997
c complex tc(18) : tensor wavefunction epsilon^mu^nu(t)
1000
double precision p(0:3), tmass
1002
double complex tc(20)
1004
double complex ft(8,4), ep(4), em(4), e0(4)
1005
double precision pt, pt2, pp, pzpt, emp, sqh, sqs
1008
double precision rZero, rHalf, rOne, rTwo
1009
parameter( rZero = 0.0d0, rHalf = 0.5d0 )
1010
parameter( rOne = 1.0d0, rTwo = 2.0d0 )
1013
parameter( stdo = 6 )
1016
sqs = sqrt(rHalf/3.d0)
1018
pt2 = p(1)**2 + p(2)**2
1019
pp = min(p(0),sqrt(pt2+p(3)**2))
1020
pt = min(pp,sqrt(pt2))
1022
c Convention for trees
1023
c ft(5,1) = dcmplx(p(0),p(3))*nst
1024
c ft(6,1) = dcmplx(p(1),p(2))*nst
1026
c Convention for loop computations
1027
ft(5,1) = dcmplx(p(0),0.D0)*nst
1028
ft(6,1) = dcmplx(p(1),0.D0)*nst
1029
ft(7,1) = dcmplx(p(2),0.D0)*nst
1030
ft(8,1) = dcmplx(p(3),0.D0)*nst
1032
if ( nhel.ge.0 ) then
1034
if ( pp.eq.rZero ) then
1035
ep(1) = dcmplx( rZero )
1036
ep(2) = dcmplx( -sqh )
1037
ep(3) = dcmplx( rZero , nst*sqh )
1038
ep(4) = dcmplx( rZero )
1040
ep(1) = dcmplx( rZero )
1041
ep(4) = dcmplx( pt/pp*sqh )
1042
if ( pt.ne.rZero ) then
1043
pzpt = p(3)/(pp*pt)*sqh
1044
ep(2) = dcmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh )
1045
ep(3) = dcmplx( -p(2)*pzpt , nst*p(1)/pt*sqh )
1047
ep(2) = dcmplx( -sqh )
1048
ep(3) = dcmplx( rZero , nst*sign(sqh,p(3)) )
1053
if ( nhel.le.0 ) then
1055
if ( pp.eq.rZero ) then
1056
em(1) = dcmplx( rZero )
1057
em(2) = dcmplx( sqh )
1058
em(3) = dcmplx( rZero , nst*sqh )
1059
em(4) = dcmplx( rZero )
1061
em(1) = dcmplx( rZero )
1062
em(4) = dcmplx( -pt/pp*sqh )
1063
if ( pt.ne.rZero ) then
1064
pzpt = -p(3)/(pp*pt)*sqh
1065
em(2) = dcmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh )
1066
em(3) = dcmplx( -p(2)*pzpt , nst*p(1)/pt*sqh )
1068
em(2) = dcmplx( sqh )
1069
em(3) = dcmplx( rZero , nst*sign(sqh,p(3)) )
1074
if ( abs(nhel).le.1 ) then
1076
if ( pp.eq.rZero ) then
1077
e0(1) = dcmplx( rZero )
1078
e0(2) = dcmplx( rZero )
1079
e0(3) = dcmplx( rZero )
1080
e0(4) = dcmplx( rOne )
1082
emp = p(0)/(tmass*pp)
1083
e0(1) = dcmplx( pp/tmass )
1084
e0(4) = dcmplx( p(3)*emp )
1085
if ( pt.ne.rZero ) then
1086
e0(2) = dcmplx( p(1)*emp )
1087
e0(3) = dcmplx( p(2)*emp )
1089
e0(2) = dcmplx( rZero )
1090
e0(3) = dcmplx( rZero )
1095
if ( nhel.eq.2 ) then
1098
ft(i,j) = ep(i)*ep(j)
1101
else if ( nhel.eq.-2 ) then
1104
ft(i,j) = em(i)*em(j)
1107
else if (tmass.eq.0) then
1113
else if (tmass.ne.0) then
1114
if ( nhel.eq.1 ) then
1117
ft(i,j) = sqh*( ep(i)*e0(j) + e0(i)*ep(j) )
1120
else if ( nhel.eq.0 ) then
1123
ft(i,j) = sqs*( ep(i)*em(j) + em(i)*ep(j)
1124
& + rTwo*e0(i)*e0(j) )
1127
else if ( nhel.eq.-1 ) then
1130
ft(i,j) = sqh*( em(i)*e0(j) + e0(i)*em(j) )
1134
write(stdo,*) 'invalid helicity in TXXXXX'
1165
subroutine mp_txxxxx(p,tmass,nhel,nst , tc)
1167
c This subroutine computes a TENSOR wavefunction.
1170
c real p(0:3) : four-momentum of tensor boson
1171
c real tmass : mass of tensor boson
1172
c integer nhel : helicity of tensor boson
1173
c = -2,-1,0,1,2 : (0 is forbidden if tmass=0.0)
1174
c integer nst = -1 or 1 : +1 for final, -1 for initial
1177
c complex tc(18) : tensor wavefunction epsilon^mu^nu(t)
1180
real*16 p(0:3), tmass
1184
complex*32 ft(8,4), ep(4), em(4), e0(4)
1185
real*16 pt, pt2, pp, pzpt, emp, sqh, sqs
1188
real*16 rZero, rHalf, rOne, rTwo
1189
parameter( rZero = 0.0e0_16, rHalf = 0.5e0_16 )
1190
parameter( rOne = 1.0e0_16, rTwo = 2.0e0_16 )
1193
parameter( stdo = 6 )
1196
sqs = sqrt(rHalf/3.0e0_16)
1198
pt2 = p(1)**2 + p(2)**2
1199
pp = min(p(0),sqrt(pt2+p(3)**2))
1200
pt = min(pp,sqrt(pt2))
1202
c Convention for trees
1203
c ft(5,1) = dcmplx(p(0),p(3))*nst
1204
c ft(6,1) = dcmplx(p(1),p(2))*nst
1206
c Convention for loop computations
1207
ft(5,1) = cmplx(p(0),0.0e0_16)*nst
1208
ft(6,1) = cmplx(p(1),0.0e0_16)*nst
1209
ft(7,1) = cmplx(p(2),0.0e0_16)*nst
1210
ft(8,1) = cmplx(p(3),0.0e0_16)*nst
1212
if ( nhel.ge.0 ) then
1214
if ( pp.eq.rZero ) then
1215
ep(1) = cmplx( rZero ,KIND=16)
1216
ep(2) = cmplx( -sqh ,KIND=16)
1217
ep(3) = cmplx( rZero , nst*sqh ,KIND=16)
1218
ep(4) = cmplx( rZero ,KIND=16)
1220
ep(1) = cmplx( rZero ,KIND=16)
1221
ep(4) = cmplx( pt/pp*sqh ,KIND=16)
1222
if ( pt.ne.rZero ) then
1223
pzpt = p(3)/(pp*pt)*sqh
1224
ep(2) = cmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh ,KIND=16)
1225
ep(3) = cmplx( -p(2)*pzpt , nst*p(1)/pt*sqh ,KIND=16)
1227
ep(2) = cmplx( -sqh ,KIND=16)
1228
ep(3) = cmplx( rZero , nst*sign(sqh,p(3)) ,KIND=16)
1233
if ( nhel.le.0 ) then
1235
if ( pp.eq.rZero ) then
1236
em(1) = cmplx( rZero ,KIND=16)
1237
em(2) = cmplx( sqh ,KIND=16)
1238
em(3) = cmplx( rZero , nst*sqh ,KIND=16)
1239
em(4) = cmplx( rZero ,KIND=16)
1241
em(1) = cmplx( rZero ,KIND=16)
1242
em(4) = cmplx( -pt/pp*sqh ,KIND=16)
1243
if ( pt.ne.rZero ) then
1244
pzpt = -p(3)/(pp*pt)*sqh
1245
em(2) = cmplx( -p(1)*pzpt , -nst*p(2)/pt*sqh ,KIND=16)
1246
em(3) = cmplx( -p(2)*pzpt , nst*p(1)/pt*sqh ,KIND=16)
1248
em(2) = cmplx( sqh ,KIND=16)
1249
em(3) = cmplx( rZero , nst*sign(sqh,p(3)) ,KIND=16)
1254
if ( abs(nhel).le.1 ) then
1256
if ( pp.eq.rZero ) then
1257
e0(1) = cmplx( rZero ,KIND=16)
1258
e0(2) = cmplx( rZero ,KIND=16)
1259
e0(3) = cmplx( rZero ,KIND=16)
1260
e0(4) = cmplx( rOne ,KIND=16)
1262
emp = p(0)/(tmass*pp)
1263
e0(1) = cmplx( pp/tmass ,KIND=16)
1264
e0(4) = cmplx( p(3)*emp ,KIND=16)
1265
if ( pt.ne.rZero ) then
1266
e0(2) = cmplx( p(1)*emp ,KIND=16)
1267
e0(3) = cmplx( p(2)*emp ,KIND=16)
1269
e0(2) = cmplx( rZero ,KIND=16)
1270
e0(3) = cmplx( rZero ,KIND=16)
1275
if ( nhel.eq.2 ) then
1278
ft(i,j) = ep(i)*ep(j)
1281
else if ( nhel.eq.-2 ) then
1284
ft(i,j) = em(i)*em(j)
1287
else if (tmass.eq.0) then
1293
else if (tmass.ne.0) then
1294
if ( nhel.eq.1 ) then
1297
ft(i,j) = sqh*( ep(i)*e0(j) + e0(i)*ep(j) )
1300
else if ( nhel.eq.0 ) then
1303
ft(i,j) = sqs*( ep(i)*em(j) + em(i)*ep(j)
1304
& + rTwo*e0(i)*e0(j) )
1307
else if ( nhel.eq.-1 ) then
1310
ft(i,j) = sqh*( em(i)*e0(j) + e0(i)*em(j) )
1314
write(stdo,*) 'invalid helicity in TXXXXX'
1344
subroutine vxxxxx(p,vmass,nhel,nsv , vc)
1346
c This subroutine computes a VECTOR wavefunction.
1349
c real p(0:3) : four-momentum of vector boson
1350
c real vmass : mass of vector boson
1351
c integer nhel = -1, 0, 1: helicity of vector boson
1352
c (0 is forbidden if vmass=0.0)
1353
c integer nsv = -1 or 1 : +1 for final, -1 for initial
1356
c complex vc(6) : vector wavefunction epsilon^mu(v)
1359
double complex vc(8)
1360
double precision p(0:3),vmass,hel,hel0,pt,pt2,pp,pzpt,emp,sqh
1361
integer nhel,nsv,nsvahl
1363
double precision rZero, rHalf, rOne, rTwo
1364
parameter( rZero = 0.0d0, rHalf = 0.5d0 )
1365
parameter( rOne = 1.0d0, rTwo = 2.0d0 )
1368
c double precision p2
1369
c double precision epsi
1370
c parameter( epsi = 2.0d-5 )
1372
c parameter( stdo = 6 )
1376
c pp = sqrt(p(1)**2+p(2)**2+p(3)**2)
1377
c if ( abs(p(0))+pp.eq.rZero ) then
1379
c & ' helas-error : p(0:3) in vxxxxx is zero momentum'
1381
c if ( p(0).le.rZero ) then
1383
c & ' helas-error : p(0:3) in vxxxxx has non-positive energy'
1385
c & ' : p(0) = ',p(0)
1387
c p2 = (p(0)+pp)*(p(0)-pp)
1388
c if ( abs(p2-vmass**2).gt.p(0)**2*2.e-5 ) then
1390
c & ' helas-error : p(0:3) in vxxxxx has inappropriate mass'
1392
c & ' : p**2 = ',p2,' : vmass**2 = ',vmass**2
1394
c if ( vmass.ne.rZero ) then
1395
c if ( abs(nhel).gt.1 ) then
1396
c write(stdo,*) ' helas-error : nhel in vxxxxx is not -1,0,1'
1397
c write(stdo,*) ' : nhel = ',nhel
1400
c if ( abs(nhel).ne.1 ) then
1401
c write(stdo,*) ' helas-error : nhel in vxxxxx is not -1,1'
1402
c write(stdo,*) ' : nhel = ',nhel
1405
c if ( abs(nsv).ne.1 ) then
1406
c write(stdo,*) ' helas-error : nsv in vmxxxx is not -1,1'
1407
c write(stdo,*) ' : nsv = ',nsv
1413
nsvahl = nsv*dabs(hel)
1414
pt2 = p(1)**2+p(2)**2
1415
pp = min(p(0),dsqrt(pt2+p(3)**2))
1416
pt = min(pp,dsqrt(pt2))
1418
c Convention for trees
1419
c vc(5) = dcmplx(p(0),p(3))*nsv
1420
c vc(6) = dcmplx(p(1),p(2))*nsv
1422
c Convention for loop computations
1423
vc(1) = dcmplx(p(0),0.D0)*nsv
1424
vc(2) = dcmplx(p(1),0.D0)*nsv
1425
vc(3) = dcmplx(p(2),0.D0)*nsv
1426
vc(4) = dcmplx(p(3),0.D0)*nsv
1429
c nhel=4 option for scalar polarization
1430
c if( nhel.eq.4 ) then
1431
c if( vmass.eq.rZero ) then
1437
c vc(1) = p(0)/vmass
1438
c vc(2) = p(1)/vmass
1439
c vc(3) = p(2)/vmass
1440
c vc(4) = p(3)/vmass
1446
if ( vmass.ne.rZero ) then
1448
hel0 = rOne-dabs(hel)
1450
if ( pp.eq.rZero ) then
1452
vc(5) = dcmplx( rZero )
1453
vc(6) = dcmplx(-hel*sqh )
1454
vc(7) = dcmplx( rZero , nsvahl*sqh )
1455
vc(8) = dcmplx( hel0 )
1459
emp = p(0)/(vmass*pp)
1460
vc(5) = dcmplx( hel0*pp/vmass )
1461
vc(8) = dcmplx( hel0*p(3)*emp+hel*pt/pp*sqh )
1462
if ( pt.ne.rZero ) then
1463
pzpt = p(3)/(pp*pt)*sqh*hel
1464
vc(6) = dcmplx( hel0*p(1)*emp-p(1)*pzpt ,
1465
& -nsvahl*p(2)/pt*sqh )
1466
vc(7) = dcmplx( hel0*p(2)*emp-p(2)*pzpt ,
1467
& nsvahl*p(1)/pt*sqh )
1469
vc(6) = dcmplx( -hel*sqh )
1470
vc(7) = dcmplx( rZero , nsvahl*sign(sqh,p(3)) )
1478
pt = sqrt(p(1)**2+p(2)**2)
1479
vc(5) = dcmplx( rZero )
1480
vc(8) = dcmplx( hel*pt/pp*sqh )
1481
if ( pt.ne.rZero ) then
1482
pzpt = p(3)/(pp*pt)*sqh*hel
1483
vc(6) = dcmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh )
1484
vc(7) = dcmplx( -p(2)*pzpt , nsv*p(1)/pt*sqh )
1486
vc(6) = dcmplx( -hel*sqh )
1487
vc(7) = dcmplx( rZero , nsv*sign(sqh,p(3)) )
1495
subroutine mp_vxxxxx(p,vmass,nhel,nsv , vc)
1497
c This subroutine computes a VECTOR wavefunction in quadruple precision.
1500
c real p(0:3) : four-momentum of vector boson
1501
c real vmass : mass of vector boson
1502
c integer nhel = -1, 0, 1: helicity of vector boson
1503
c (0 is forbidden if vmass=0.0)
1504
c integer nsv = -1 or 1 : +1 for final, -1 for initial
1507
c complex vc(6) : vector wavefunction epsilon^mu(v)
1511
real*16 p(0:3),vmass,hel,hel0,pt,pt2,pp,pzpt,emp,sqh
1512
integer nhel,nsv,nsvahl
1514
real*16 rZero, rHalf, rOne, rTwo
1515
parameter( rZero = 0.0e0_16, rHalf = 0.5e0_16 )
1516
parameter( rOne = 1.0e0_16, rTwo = 2.0e0_16 )
1519
hel = real(nhel,KIND=16)
1520
nsvahl = nsv*abs(hel)
1521
pt2 = p(1)**2+p(2)**2
1522
pp = min(p(0),sqrt(pt2+p(3)**2))
1523
pt = min(pp,sqrt(pt2))
1525
c Convention for loop computations
1526
vc(1) = cmplx(p(0),0.e0_16, KIND=16)*nsv
1527
vc(2) = cmplx(p(1),0.e0_16, KIND=16)*nsv
1528
vc(3) = cmplx(p(2),0.e0_16, KIND=16)*nsv
1529
vc(4) = cmplx(p(3),0.e0_16, KIND=16)*nsv
1531
if ( vmass.ne.rZero ) then
1533
hel0 = rOne-abs(hel)
1535
if ( pp.eq.rZero ) then
1537
vc(5) = cmplx( rZero , KIND=16)
1538
vc(6) = cmplx(-hel*sqh , KIND=16)
1539
vc(7) = cmplx( rZero , nsvahl*sqh , KIND=16)
1540
vc(8) = cmplx( hel0 , KIND=16)
1544
emp = p(0)/(vmass*pp)
1545
vc(5) = cmplx( hel0*pp/vmass , KIND=16)
1546
vc(8) = cmplx( hel0*p(3)*emp+hel*pt/pp*sqh , KIND=16 )
1547
if ( pt.ne.rZero ) then
1548
pzpt = p(3)/(pp*pt)*sqh*hel
1549
vc(6) = cmplx( hel0*p(1)*emp-p(1)*pzpt ,
1550
& -nsvahl*p(2)/pt*sqh , KIND=16)
1551
vc(7) = cmplx( hel0*p(2)*emp-p(2)*pzpt ,
1552
& nsvahl*p(1)/pt*sqh , KIND=16)
1554
vc(6) = cmplx( -hel*sqh , KIND=16)
1555
vc(7) = cmplx( rZero , nsvahl*sign(sqh,p(3)) , KIND=16)
1563
pt = sqrt(p(1)**2+p(2)**2)
1564
vc(5) = cmplx( rZero , KIND=16)
1565
vc(8) = cmplx( hel*pt/pp*sqh , KIND=16)
1566
if ( pt.ne.rZero ) then
1567
pzpt = p(3)/(pp*pt)*sqh*hel
1568
vc(6) = cmplx( -p(1)*pzpt , -nsv*p(2)/pt*sqh , KIND=16)
1569
vc(7) = cmplx( -p(2)*pzpt , nsv*p(1)/pt*sqh , KIND=16)
1571
vc(6) = cmplx( -hel*sqh , KIND=16)
1572
vc(7) = cmplx( rZero , nsv*sign(sqh,p(3)), KIND=16 )
1580
subroutine boostx(p,q , pboost)
1582
c This subroutine performs the Lorentz boost of a four-momentum. The
1583
c momentum p is assumed to be given in the rest frame of q. pboost is
1584
c the momentum p boosted to the frame in which q is given. q must be a
1585
c timelike momentum.
1588
c real p(0:3) : four-momentum p in the q rest frame
1589
c real q(0:3) : four-momentum q in the boosted frame
1592
c real pboost(0:3) : four-momentum p in the boosted frame
1595
double precision p(0:3),q(0:3),pboost(0:3),pq,qq,m,lf
1597
double precision rZero
1598
parameter( rZero = 0.0d0 )
1600
qq = q(1)**2+q(2)**2+q(3)**2
1603
c if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
1605
c & ' helas-error : p(0:3) in boostx is zero momentum'
1607
c if (abs(q(0))+qq.eq.rZero) then
1609
c & ' helas-error : q(0:3) in boostx is zero momentum'
1611
c if (p(0).le.rZero) then
1613
c & ' helas-warn : p(0:3) in boostx has not positive energy'
1615
c & ' : p(0) = ',p(0)
1617
c if (q(0).le.rZero) then
1619
c & ' helas-error : q(0:3) in boostx has not positive energy'
1621
c & ' : q(0) = ',q(0)
1623
c pp=p(0)**2-p(1)**2-p(2)**2-p(3)**2
1624
c if (pp.lt.rZero) then
1626
c & ' helas-warn : p(0:3) in boostx is spacelike'
1630
c if (q(0)**2-qq.le.rZero) then
1632
c & ' helas-error : q(0:3) in boostx is not timelike'
1634
c & ' : q**2 = ',q(0)**2-qq
1636
c if (qq.eq.rZero) then
1638
c & ' helas-warn : q(0:3) in boostx has zero spacial components'
1642
if ( qq.ne.rZero ) then
1643
pq = p(1)*q(1)+p(2)*q(2)+p(3)*q(3)
1644
m = sqrt(q(0)**2-qq)
1645
lf = ((q(0)-m)*pq/qq+p(0))/m
1646
pboost(0) = (p(0)*q(0)+pq)/m
1647
pboost(1) = p(1)+q(1)*lf
1648
pboost(2) = p(2)+q(2)*lf
1649
pboost(3) = p(3)+q(3)*lf
1660
subroutine boostm(p,q,m, pboost)
1662
c This subroutine performs the Lorentz boost of a four-momentum. The
1663
c momentum p is assumed to be given in the rest frame of q. pboost is
1664
c the momentum p boosted to the frame in which q is given. q must be a
1665
c timelike momentum.
1668
c real p(0:3) : four-momentum p in the q rest frame
1669
c real q(0:3) : four-momentum q in the boosted frame
1670
c real m : mass of q (for numerical stability)
1673
c real pboost(0:3) : four-momentum p in the boosted frame
1676
double precision p(0:3),q(0:3),pboost(0:3),pq,qq,m,lf
1678
double precision rZero
1679
parameter( rZero = 0.0d0 )
1683
c parameter( stdo = 6 )
1684
c double precision pp
1687
qq = q(1)**2+q(2)**2+q(3)**2
1690
c if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
1692
c & ' helas-error : p(0:3) in boostx is zero momentum'
1694
c if (abs(q(0))+qq.eq.rZero) then
1696
c & ' helas-error : q(0:3) in boostx is zero momentum'
1698
c if (p(0).le.rZero) then
1700
c & ' helas-warn : p(0:3) in boostx has not positive energy'
1702
c & ' : p(0) = ',p(0)
1704
c if (q(0).le.rZero) then
1706
c & ' helas-error : q(0:3) in boostx has not positive energy'
1708
c & ' : q(0) = ',q(0)
1710
c pp=p(0)**2-p(1)**2-p(2)**2-p(3)**2
1711
c if (pp.lt.rZero) then
1713
c & ' helas-warn : p(0:3) in boostx is spacelike'
1717
c if (q(0)**2-qq.le.rZero) then
1719
c & ' helas-error : q(0:3) in boostx is not timelike'
1721
c & ' : q**2 = ',q(0)**2-qq
1723
c if (qq.eq.rZero) then
1725
c & ' helas-warn : q(0:3) in boostx has zero spacial components'
1729
if ( qq.ne.rZero ) then
1730
pq = p(1)*q(1)+p(2)*q(2)+p(3)*q(3)
1731
lf = ((q(0)-m)*pq/qq+p(0))/m
1732
pboost(0) = (p(0)*q(0)+pq)/m
1733
pboost(1) = p(1)+q(1)*lf
1734
pboost(2) = p(2)+q(2)*lf
1735
pboost(3) = p(3)+q(3)*lf
1746
subroutine momntx(energy,mass,costh,phi , p)
1748
c This subroutine sets up a four-momentum from the four inputs.
1751
c real energy : energy
1753
c real costh : cos(theta)
1754
c real phi : azimuthal angle
1757
c real p(0:3) : four-momentum
1760
double precision p(0:3),energy,mass,costh,phi,pp,sinth
1762
double precision rZero, rOne
1763
parameter( rZero = 0.0d0, rOne = 1.0d0 )
1766
c double precision rPi, rTwo
1767
c parameter( rPi = 3.14159265358979323846d0, rTwo = 2.d0 )
1769
c parameter( stdo = 6 )
1773
c if (energy.lt.mass) then
1775
c & ' helas-error : energy in momntx is less than mass'
1777
c & ' : energy = ',energy,' : mass = ',mass
1779
c if (mass.lt.rZero) then
1780
c write(stdo,*) ' helas-error : mass in momntx is negative'
1781
c write(stdo,*) ' : mass = ',mass
1783
c if (abs(costh).gt.rOne) then
1784
c write(stdo,*) ' helas-error : costh in momntx is improper'
1785
c write(stdo,*) ' : costh = ',costh
1787
c if (phi.lt.rZero .or. phi.gt.rTwo*rPi) then
1789
c & ' helas-warn : phi in momntx does not lie on 0.0 thru 2.0*pi'
1797
if ( energy.eq.mass ) then
1805
pp = sqrt((energy-mass)*(energy+mass))
1806
sinth = sqrt((rOne-costh)*(rOne+costh))
1808
if ( phi.eq.rZero ) then
1812
p(1) = pp*sinth*cos(phi)
1813
p(2) = pp*sinth*sin(phi)
1820
subroutine rotxxx(p,q , prot)
1822
c This subroutine performs the spacial rotation of a four-momentum.
1823
c the momentum p is assumed to be given in the frame where the spacial
1824
c component of q points the positive z-axis. prot is the momentum p
1825
c rotated to the frame where q is given.
1828
c real p(0:3) : four-momentum p in q(1)=q(2)=0 frame
1829
c real q(0:3) : four-momentum q in the rotated frame
1832
c real prot(0:3) : four-momentum p in the rotated frame
1835
double precision p(0:3),q(0:3),prot(0:3),qt2,qt,psgn,qq,p1
1837
double precision rZero, rOne
1838
parameter( rZero = 0.0d0, rOne = 1.0d0 )
1842
c parameter( stdo = 6 )
1847
qt2 = q(1)**2 + q(2)**2
1850
c if (abs(p(0))+abs(p(1))+abs(p(2))+abs(p(3)).eq.rZero) then
1852
c & ' helas-error : p(0:3) in rotxxx is zero momentum'
1854
c if (abs(q(0))+abs(q(3))+qt2.eq.rZero) then
1856
c & ' helas-error : q(0:3) in rotxxx is zero momentum'
1858
c if (qt2+abs(q(3)).eq.rZero) then
1860
c & ' helas-warn : q(0:3) in rotxxx has zero spacial momentum'
1864
if ( qt2.eq.rZero ) then
1865
if ( q(3).eq.rZero ) then
1870
psgn = dsign(rOne,q(3))
1876
qq = sqrt(qt2+q(3)**2)
1879
prot(1) = q(1)*q(3)/qq/qt*p1 -q(2)/qt*p(2) +q(1)/qq*p(3)
1880
prot(2) = q(2)*q(3)/qq/qt*p1 +q(1)/qt*p(2) +q(2)/qq*p(3)
1881
prot(3) = -qt/qq*p1 +q(3)/qq*p(3)
1887
subroutine mom2cx(esum,mass1,mass2,costh1,phi1 , p1,p2)
1889
c This subroutine sets up two four-momenta in the two particle rest
1893
c real esum : energy sum of particle 1 and 2
1894
c real mass1 : mass of particle 1
1895
c real mass2 : mass of particle 2
1896
c real costh1 : cos(theta) of particle 1
1897
c real phi1 : azimuthal angle of particle 1
1900
c real p1(0:3) : four-momentum of particle 1
1901
c real p2(0:3) : four-momentum of particle 2
1904
double precision p1(0:3),p2(0:3),
1905
& esum,mass1,mass2,costh1,phi1,md2,ed,pp,sinth1
1907
double precision rZero, rHalf, rOne, rTwo
1908
parameter( rZero = 0.0d0, rHalf = 0.5d0 )
1909
parameter( rOne = 1.0d0, rTwo = 2.0d0 )
1912
c double precision rPi
1913
c parameter( rPi = 3.14159265358979323846d0 )
1915
c parameter( stdo = 6 )
1919
c if (esum.lt.mass1+mass2) then
1921
c & ' helas-error : esum in mom2cx is less than mass1+mass2'
1923
c & ' : esum = ',esum,
1924
c & ' : mass1+mass2 = ',mass1,mass2
1926
c if (mass1.lt.rZero) then
1927
c write(stdo,*) ' helas-error : mass1 in mom2cx is negative'
1928
c write(stdo,*) ' : mass1 = ',mass1
1930
c if (mass2.lt.rZero) then
1931
c write(stdo,*) ' helas-error : mass2 in mom2cx is negative'
1932
c write(stdo,*) ' : mass2 = ',mass2
1934
c if (abs(costh1).gt.1.) then
1935
c write(stdo,*) ' helas-error : costh1 in mom2cx is improper'
1936
c write(stdo,*) ' : costh1 = ',costh1
1938
c if (phi1.lt.rZero .or. phi1.gt.rTwo*rPi) then
1940
c & ' helas-warn : phi1 in mom2cx does not lie on 0.0 thru 2.0*pi'
1942
c & ' : phi1 = ',phi1
1946
md2 = (mass1-mass2)*(mass1+mass2)
1948
if ( mass1*mass2.eq.rZero ) then
1949
pp = (esum-abs(ed))*rHalf
1951
pp = sqrt((md2/esum)**2-rTwo*(mass1**2+mass2**2)+esum**2)*rHalf
1953
sinth1 = sqrt((rOne-costh1)*(rOne+costh1))
1955
p1(0) = max((esum+ed)*rHalf,rZero)
1956
p1(1) = pp*sinth1*cos(phi1)
1957
p1(2) = pp*sinth1*sin(phi1)
1960
p2(0) = max((esum-ed)*rHalf,rZero)
1968
C###############################################################################
1969
C LOOP related universal subroutines
1970
C###############################################################################
1972
C===============================================================================
1973
C Subroutines to create the external wavefunctions of the L-cut particles
1974
C===============================================================================
1976
SUBROUTINE LCUT_F(Q,CFIG,W)
1982
CALL LCUT_V(Q,CFIG,W)
1985
SUBROUTINE LCUT_AF(Q,CFIG,W)
1991
CALL LCUT_V(Q,CFIG,W)
1994
SUBROUTINE LCUT_V(Q,CFIG,W)
2004
W(CFIG+4)=(1.d0,0.d0)
2013
SUBROUTINE LCUT_S(Q,CFIG,W)
2028
SUBROUTINE LCUT_AS(Q,CFIG,W)
2043
SUBROUTINE MP_LCUT_F(Q,CFIG,W)
2049
CALL MP_LCUT_V(Q,CFIG,W)
2052
SUBROUTINE MP_LCUT_AF(Q,CFIG,W)
2058
CALL MP_LCUT_V(Q,CFIG,W)
2061
SUBROUTINE MP_LCUT_V(Q,CFIG,W)
2066
COMPLEX*32 IONE, IZERO
2067
PARAMETER (IONE=(1.0e0_16,0.0e0_16))
2068
PARAMETER (IZERO=(0.0e0_16,0.0e0_16))
2083
SUBROUTINE MP_LCUT_AS(Q,CFIG,W)
2089
PARAMETER (IONE=(1.0e0_16,0.0e0_16))
2100
SUBROUTINE MP_LCUT_S(Q,CFIG,W)
2106
PARAMETER (IONE=(1.0e0_16,0.0e0_16))
2117
C===============================================================================
2118
C Subroutines to close the lorentz traces of loops,
2119
C===============================================================================
2121
SUBROUTINE CLOSE_4(AMPS,RES)
2126
RES=AMPS(1)+AMPS(2)+AMPS(3)+AMPS(4)
2130
SUBROUTINE CLOSE_1(AMPS,RES)
2139
SUBROUTINE MP_CLOSE_4(AMPS,RES)
2144
RES=AMPS(1)+AMPS(2)+AMPS(3)+AMPS(4)
2148
SUBROUTINE MP_CLOSE_1(AMPS,RES)
2157
C===============================================================================
2158
C OLD Subroutines to close the lorentz traces of loops,
2160
C===============================================================================
2162
SUBROUTINE CLOSE_V(Q,M,AMPS,RES)
2169
IF (M.NE.(0.D0,0.0D0)) THEN
2170
STOP 'Massive vector L-cut particle not supported'
2172
RES=AMPS(1)-AMPS(2)-AMPS(3)-AMPS(4)
2176
c This subroutine is to recreate the fermion propagator with 4 helicities
2177
c only. This has problems with certain configuration of the imaginary
2178
c momentum q, so it is not implemented yet.
2180
SUBROUTINE CLOSE_F4HEL(Q,M,AMPS,RES)
2183
COMPLEX*16 RES, QNORM
2190
write(*,*) 'PPM=',PPM
2191
write(*,*) 'PMM=',PMM
2193
QNORM=SQRT(Q(0)**2-Q(1)**2-Q(2)**2-Q(3)**2)
2195
write(*,*) 'QNORM=',QNORM
2197
RES=(0.D0,0.5D0)*((PPM+PMM)+(PPM-PMM)*(M/QNORM))
2198
write(*,*) 'RES1=',RES
2200
RES=(0.D0,0.5D0)*(PPM+PMM)
2201
write(*,*) 'RES2=',RES
2206
SUBROUTINE CLOSE_F(Q,M,AMPS,RES)
2214
RES=(Q(0)-Q(3))*AMPS(1)+
2215
& (-Q(1)+(0.0d0,1.0d0)*Q(2))*AMPS(2)+
2216
& (-Q(1)-(0.0d0,1.0d0)*Q(2))*AMPS(3)+
2217
& (Q(0)+Q(3))*AMPS(4)+
2218
& (Q(0)+Q(3))*AMPS(5)+
2219
& (Q(1)-(0.0d0,1.0d0)*Q(2))*AMPS(6)+
2220
& (Q(1)+(0.0d0,1.0d0)*Q(2))*AMPS(7)+
2221
& (Q(0)-Q(3))*AMPS(8)
2222
IF (M.NE.(0.D0,0.D0)) THEN
2223
RES=RES-M*(AMPS(9)+AMPS(10)+AMPS(11)+AMPS(12))
2228
SUBROUTINE CLOSE_S(Q,AMP,RES)
2238
c // QUAD PREC VERSIONS OF THE ABOVE //
2240
SUBROUTINE MP_CLOSE_V(Q,M,AMPS,RES)
2247
IF (M.NE.(0.0e0_16,0.0e0_16)) THEN
2248
STOP 'Massive vector L-cut particle not supported'
2250
RES=AMPS(1)-AMPS(2)-AMPS(3)-AMPS(4)
2254
SUBROUTINE MP_CLOSE_F(Q,M,AMPS,RES)
2261
RES=(0.0e0_16,0.0e0_16)
2262
RES=(Q(0)-Q(3))*AMPS(1)+
2263
& (-Q(1)+(0.0e0_16,1.0e0_16)*Q(2))*AMPS(2)+
2264
& (-Q(1)-(0.0e0_16,1.0e0_16)*Q(2))*AMPS(3)+
2265
& (Q(0)+Q(3))*AMPS(4)+
2266
& (Q(0)+Q(3))*AMPS(5)+
2267
& (Q(1)-(0.0e0_16,1.0e0_16)*Q(2))*AMPS(6)+
2268
& (Q(1)+(0.0e0_16,1.0e0_16)*Q(2))*AMPS(7)+
2269
& (Q(0)-Q(3))*AMPS(8)
2270
IF (M.NE.(0.0e0_16,0.0e0_16)) THEN
2271
RES=RES-M*(AMPS(9)+AMPS(10)+AMPS(11)+AMPS(12))
2276
SUBROUTINE MP_CLOSE_S(Q,AMP,RES)
2282
RES=(1.0e0_16,0.0e0_16)*AMP
2286
C===============================================================================
2287
C OLD Subroutines to create the external wavefunctions of the L-cut particles
2289
C===============================================================================
2291
c This subroutine is to recreate the fermion propagator with 4 helicities
2292
c only. This has problems with certain configuration of the imaginary
2293
c momentum q, so it is not implemented yet.
2295
SUBROUTINE LCUT_F4HEL(Q,M,CFIG,SCD,W)
2306
CALL ILXXXX(Q(0),M,-1,1,W(1))
2309
CALL ILXXXX(Q(0),M,-1,-1,W(1))
2311
write(*,*) 'Wcf(',j,',1)=',W(1)
2314
ELSEIF (CFIG.EQ.2) THEN
2317
CALL ILXXXX(Q(0),M,1,1,W(1))
2320
CALL ILXXXX(Q(0),M,1,-1,W(1))
2322
ELSEIF (CFIG.EQ.3) THEN
2325
CALL OLXXXX(Q(0),M,-1,-1,W(1))
2328
CALL OLXXXX(Q(0),M,-1,1,W(1))
2330
ELSEIF (CFIG.EQ.4) THEN
2333
CALL OLXXXX(Q(0),M,1,-1,W(1))
2336
CALL OLXXXX(Q(0),M,1,1,W(1))
2339
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT SPINORS
2349
SUBROUTINE OLD_LCUT_F(Q,M,CFIG,SCD,W)
2368
ELSEIF (CFIG.eq.2) then
2374
ELSEIF (CFIG.eq.3) then
2380
ELSEIF (CFIG.eq.4) then
2386
ELSEIF (CFIG.eq.5) then
2392
ELSEIF (CFIG.eq.6) then
2398
ELSEIF (CFIG.eq.7) then
2404
ELSEIF (CFIG.eq.8) then
2410
ELSEIF (CFIG.eq.9) then
2416
ELSEIF (CFIG.eq.10) then
2422
ELSEIF (CFIG.eq.11) then
2428
ELSEIF (CFIG.eq.12) then
2435
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT SPINORS
2450
SUBROUTINE OLD_MP_LCUT_F(Q,M,CFIG,SCD,W)
2458
PARAMETER (IZERO=(0.0e0_16,0.0e0_16))
2460
PARAMETER (IONE=(1.0e0_16,0.0e0_16))
2473
ELSEIF (CFIG.eq.2) then
2479
ELSEIF (CFIG.eq.3) then
2485
ELSEIF (CFIG.eq.4) then
2491
ELSEIF (CFIG.eq.5) then
2497
ELSEIF (CFIG.eq.6) then
2503
ELSEIF (CFIG.eq.7) then
2509
ELSEIF (CFIG.eq.8) then
2515
ELSEIF (CFIG.eq.9) then
2521
ELSEIF (CFIG.eq.10) then
2527
ELSEIF (CFIG.eq.11) then
2533
ELSEIF (CFIG.eq.12) then
2540
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT SPINORS
2555
SUBROUTINE OLD_LCUT_CF(Q,M,CFIG,SCD,W)
2566
CALL ICLXXX(Q(0),M,-1,1,W(1))
2569
CALL ICLXXX(Q(0),M,-1,-1,W(1))
2571
ELSEIF (CFIG.EQ.2) THEN
2574
CALL ICLXXX(Q(0),M,1,1,W(1))
2577
CALL ICLXXX(Q(0),M,1,-1,W(1))
2579
ELSEIF (CFIG.EQ.3) THEN
2582
CALL OCLXXX(Q(0),M,-1,-1,W(1))
2585
CALL OCLXXX(Q(0),M,-1,1,W(1))
2587
ELSEIF (CFIG.EQ.4) THEN
2590
CALL OCLXXX(Q(0),M,1,-1,W(1))
2593
CALL OCLXXX(Q(0),M,1,1,W(1))
2596
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT SPINORS
2606
SUBROUTINE OLD_LCUT_V(Q,M,CFIG,SCD,W)
2615
STOP 'Massive vector L-cut particle not supported'
2623
ELSEIF (CFIG.EQ.2) THEN
2625
ELSEIF (CFIG.EQ.3) THEN
2627
ELSEIF (CFIG.EQ.4) THEN
2630
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT VECTORS
2645
SUBROUTINE OLD_MP_LCUT_V(Q,M,CFIG,SCD,W)
2653
PARAMETER (IZERO=(0.0e0_16,0.0e0_16))
2655
PARAMETER (IONE=(1.0e0_16,0.0e0_16))
2658
IF (M.NE.IZERO) THEN
2659
STOP 'Massive vector L-cut particle not supported'
2667
ELSEIF (CFIG.EQ.2) THEN
2669
ELSEIF (CFIG.EQ.3) THEN
2671
ELSEIF (CFIG.EQ.4) THEN
2674
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND L-CUT VECTORS
2689
SUBROUTINE OLD_LCUT_S(Q,M,CFIG,SCD,W)
2699
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND SCALAR
2714
SUBROUTINE OLD_MP_LCUT_S(Q,M,CFIG,SCD,W)
2722
PARAMETER (IONE=(1.0e0_16,0.0e0_16))
2726
C REVERSE THE MOMENTUM IN THE WF FOR THE SECOND SCALAR
2741
C===============================================================================
2742
C Complex-momentum version of the subroutine to create on-shell
2743
C wavefunctions of particles with different spins. OBSOLETE
2744
C===============================================================================
2746
C The subroutine with charge conjugation are not yet implemented
2747
c Obsolete by now too
2749
subroutine oclxxx(p,ffmass,nhel,nsf,fo)
2752
double complex fo(8),p(0:3)
2753
double precision ffmass
2756
CALL olxxxx(p,ffmass,nhel,nsf,fo)
2760
subroutine iclxxx(p,ffmass,nhel,nsf,fi)
2763
double complex fi(8),p(0:3)
2764
double precision ffmass
2767
CALL ilxxxx(p,ffmass,nhel,nsf,fi)
2771
subroutine ilxxxx(p,ffmass,nhel,nsf,fi)
2773
c This subroutine computes a fermion wavefunction with the flowing-IN
2774
c fermion number and defined with complex ONSHELL momentium.
2777
c complex p(0:3) : four-momentum of fermion
2778
c real fmass : mass of fermion
2779
c integer nhel = -1 or 1 : helicity of fermion
2780
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
2783
c complex fi(8) : fermion wavefunction |fi>
2784
c Note: There are 4 components for the spinor and four for the
2787
double complex fi(8),chi(2), fmass
2788
c double precision p(0:3),sf(2),sfomeg(2),omega(2),fmass,
2789
c & pp,pp3,sqp0p3,sqm(0:1)
2790
double complex sqm(0:1)
2791
double precision sf(2),ffmass
2792
double complex p(0:3), sfomeg(2),omega(2),
2794
integer nhel,nsf,ip,im,nh
2796
double precision rZero, rHalf, rTwo
2797
parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
2801
c fi(5) = dcmplx(p(0),p(3))*nsf
2802
c fi(6) = dcmplx(p(1),p(2))*nsf
2810
fmass = sqrt(p(0)**2-p(1)**2-p(2)**2-p(3)**2)
2812
if ( ffmass.ne.rZero ) then
2813
c special treatment for massless particles.
2814
c pp = min(p(0),sqrt(p(1)**2+p(2)**2+p(3)**2))
2815
pp=sqrt(p(1)**2+p(2)**2+p(3)**2)
2816
c for time-like four-momenta we can always think of it as the p_vec^2
2817
if ( abs(pp).eq.rZero ) then
2819
sqm(0) = sqrt(fmass) ! possibility of negative fermion masses
2820
sqm(1) = sqm(0) ! possibility of negative fermion masses
2824
fi(1) = ip * sqm(ip)
2825
fi(2) = im*nsf * sqm(ip)
2826
fi(3) = ip*nsf * sqm(im)
2827
fi(4) = im * sqm(im)
2832
pp=sqrt(p(1)**2+p(2)**2+p(3)**2)
2833
write(*,*) 'ppre=',pp
2834
c if( (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) .or.
2835
c & (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) ) then
2838
sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
2839
c fermion spin using HELAS conventions.
2840
sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
2841
omega(1) = sqrt(p(0)+pp)
2842
c the omega of the definition.
2843
c omega(2) = fmass/omega(1)
2844
omega(2) = sqrt(p(0)-pp)
2848
sfomeg(1) = sf(1)*omega(ip)
2849
sfomeg(2) = sf(2)*omega(im)
2850
c pp3 = max(pp+p(3),rZero)
2852
chi(1) = sqrt(pp3*rHalf/pp)
2853
if ( abs(pp3).eq.rZero ) then
2854
chi(2) = dcmplx(-nh )
2856
chi(2) = ( (nh*p(1)) + ((0d0,1d0)*p(2)) )/
2861
c Write(*,*) 'Chi=',Chi(1),' and ',Chi(2)
2863
fi(1) = sfomeg(1)*chi(im)
2864
fi(2) = sfomeg(1)*chi(ip)
2865
c Write(*,*) 'fi(2)=',fi(2)
2866
fi(3) = sfomeg(2)*chi(im)
2867
c Write(*,*) 'fi(3)=',fi(3)
2868
fi(4) = sfomeg(2)*chi(ip)
2874
c if(zabs(p(1)).eq.0d0.and.zabs(p(2)).eq.0d0.and.
2875
c .zabs(p(3)).lt.0d0) then
2878
sqp0p3 = sqrt(p(0)+p(3))*nsf
2881
if ( abs(sqp0p3).eq.rZero ) then
2882
chi(2) = dcmplx(-nhel )*sqrt(rTwo*p(0))
2884
chi(2) = ( nh*p(1) + ((0d0,1d0)*p(2) ) )/sqp0p3
2887
fi(1) = dcmplx( rZero )
2888
fi(2) = dcmplx( rZero )
2894
fi(3) = dcmplx( rZero )
2895
fi(4) = dcmplx( rZero )
2902
subroutine olxxxx(p,ffmass,nhel,nsf,fo)
2904
c This subroutine computes a fermion wavefunction with the flowing-OUT
2905
c fermion number and defined with complex ONSHELL momentum.
2908
c complex p(0:3) : four-momentum of fermion
2909
c real fmass : mass of fermion
2910
c integer nhel = -1 or 1 : helicity of fermion
2911
c integer nsf = -1 or 1 : +1 for particle, -1 for anti-particle
2914
c complex fo(8) : fermion wavefunction <fo|
2917
double complex fo(8),chi(2), fmass
2919
double complex sqm(0:1)
2920
double precision sf(2), ffmass
2922
double complex p(0:3),sfomeg(2),omega(2),
2924
integer nhel,nsf,nh,ip,im
2926
double precision rZero, rHalf, rTwo
2927
parameter( rZero = 0.0d0, rHalf = 0.5d0, rTwo = 2.0d0 )
2930
c fo(5) = dcmplx(p(0),p(3))*nsf
2931
c fo(6) = dcmplx(p(1),p(2))*nsf
2939
fmass=sqrt(p(0)**2-p(1)**2-p(2)**2-p(3)**2)
2941
if ( ffmass.ne.rZero ) then
2943
c pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
2944
pp=sqrt(p(1)**2+p(2)**2+p(3)**2)
2946
if ( abs(pp).eq.rZero ) then
2948
sqm(0) = sqrt(fmass) ! possibility of negative fermion masses
2949
sqm(1) = sqm(0) ! possibility of negative fermion masses
2953
fo(1) = im * sqm(im)
2954
fo(2) = ip*nsf * sqm(im)
2955
fo(3) = im*nsf * sqm(-ip)
2956
fo(4) = ip * sqm(-ip)
2960
c pp = min(p(0),dsqrt(p(1)**2+p(2)**2+p(3)**2))
2961
pp=sqrt(p(1)**2+p(2)**2+p(3)**2) !repetition
2962
c if( (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) .or.
2963
c & (dble(p(0)) .lt. 0 .and. dble(pp) .gt. 0) ) then
2966
sf(1) = dble(1+nsf+(1-nsf)*nh)*rHalf
2967
sf(2) = dble(1+nsf-(1-nsf)*nh)*rHalf
2968
omega(1) = sqrt(p(0)+pp)
2969
c omega(2) = fmass/omega(1)
2970
omega(2) = sqrt(p(0)-pp)
2973
sfomeg(1) = sf(1)*omega(ip)
2974
sfomeg(2) = sf(2)*omega(im)
2976
chi(1) = sqrt(pp3*rHalf/pp)
2977
if ( abs(pp3).eq.rZero ) then
2978
chi(2) = dcmplx(-nh )
2980
chi(2) = ( nh*p(1) - ((0d0,1d0)*p(2)) )/
2984
fo(1) = sfomeg(2)*chi(im)
2985
fo(2) = sfomeg(2)*chi(ip)
2986
fo(3) = sfomeg(1)*chi(im)
2987
fo(4) = sfomeg(1)*chi(ip)
2993
c if(zabs(p(1)).eq.0d0.and.zabs(p(2)).eq.0d0.and.zabs(p(3)).lt.0d0) then
2996
sqp0p3 = sqrt(p(0)+p(3))*nsf
2999
if ( abs(sqp0p3).eq.rZero ) then
3000
chi(2) = dcmplx(-nhel )*sqrt(rTwo*p(0))
3002
chi(2) = ( nh*p(1)-((0d0,1d0)*p(2)) )/sqp0p3
3007
fo(3) = dcmplx( rZero )
3008
fo(4) = dcmplx( rZero )
3010
fo(1) = dcmplx( rZero )
3011
fo(2) = dcmplx( rZero )
added
removed
aloha/template_files/aloha_functions_loop.f
