1
subroutine jwwwnx(w1,w2,w3,gwwa,gwwz,wmass,wwidth , jwww)
3
c This subroutine computes an off-shell W+/W- current from the four-
4
c point gauge boson coupling. The vector propagators for the output
5
c W and the internal Z bosons are given in unitary gauge, and that of
6
c the internal photon is given in Feynman gauge.
9
c complex w1(6) : first vector w1
10
c complex w2(6) : second vector w2
11
c complex w3(6) : third vector w3
12
c real gwwa : coupling constant of W and A gwwa
13
c real gwwz : coupling constant of W and Z gwwz
14
c real zmass : mass of internal Z
15
c real zwidth : width of internal Z
16
c real wmass : mass of output W
17
c real wwidth : width of output W
19
c the possible sets of the inputs are as follows:
20
c -------------------------------------------------------------------
21
c | w1 | w2 | w3 |gwwa|gwwz|zmass|zwidth|wmass|wwidth || jwww |
22
c -------------------------------------------------------------------
23
c | W- | W+ | W- |gwwa|gwwz|zmass|zwidth|wmass|wwidth || W+ |
24
c | W+ | W- | W+ |gwwa|gwwz|zmass|zwidth|wmass|wwidth || W- |
25
c -------------------------------------------------------------------
26
c where all the bosons are defined by the flowing-OUT quantum number.
29
c complex jwww(6) : W current j^mu(w':w1,w2,w3)
32
double complex w1(6),w2(6),w3(6),jwww(6)
33
double complex dw1(0:3),dw2(0:3),dw3(0:3),jj(0:3)
34
double complex dw,w12,w32,w13,jq
35
double complex cm2 ! mass**2- I Gamma mass (Fabio)
36
double precision gwwa,gwwz,wmass,wwidth
37
double precision p1(0:3),p2(0:3),p3(0:3),q(0:3)
38
double precision dgwwa2,dgwwz2,dgw2,dmw,dww,mw2,q2
40
double precision rZero, rOne, rTwo
41
parameter( rZero = 0.0d0, rOne = 1.0d0, rTwo = 2.0d0 )
49
if ( abs(w1(1))+abs(w1(2))+abs(w1(3))+abs(w1(4)).eq.rZero ) then
50
write(stdo,*) ' helas-warn : w1 in jwwwxx is zero vector'
52
if ( abs(w1(5))+abs(w1(6)).eq.rZero ) then
54
& ' helas-error : w1 in jwwwxx has zero momentum'
56
if ( abs(w2(1))+abs(w2(2))+abs(w2(3))+abs(w2(4)).eq.rZero ) then
57
write(stdo,*) ' helas-warn : w2 in jwwwxx is zero vector'
59
if ( abs(w2(5))+abs(w2(6)).eq.rZero ) then
61
& ' helas-error : w2 in jwwwxx has zero momentum'
63
if ( abs(w3(1))+abs(w3(2))+abs(w3(3))+abs(w3(4)).eq.rZero ) then
64
write(stdo,*) ' helas-warn : w3 in jwwwxx is zero vector'
66
if ( abs(w3(5))+abs(w3(6)).eq.rZero ) then
68
& ' helas-error : w3 in jwwwxx has zero momentum'
70
c Neil edited the following to allow 3-site couplings.
71
c if ( gwwa.eq.rZero ) then
72
c write(stdo,*) ' helas-error : gwwa in jwwwxx is zero coupling'
74
c if ( gwwz.eq.rZero ) then
75
c write(stdo,*) ' helas-error : gwwz in jwwwxx is zero coupling'
77
c if ( gwwa.lt.rZero .or. gwwa.ge.gwwz ) then
79
c & ' helas-warn : gwwa/gwwz in jwwwxx are non-standard couplings'
81
c & ' : gwwa = ',gwwa,' gwwz = ',gwwz
84
if ( wmass.le.rZero ) then
85
write(stdo,*) ' helas-error : wmass in jwwwxx is not positive'
86
write(stdo,*) ' : wmass = ',wmass
88
if ( wwidth.lt.rZero ) then
89
write(stdo,*) ' helas-error : wwidth in jwwwxx is negative'
90
write(stdo,*) ' : wwidth = ',wwidth
97
jwww(5) = w1(5)+w2(5)+w3(5)
98
jwww(6) = w1(6)+w2(6)+w3(6)
100
dw1(0) = dcmplx(w1(1))
101
dw1(1) = dcmplx(w1(2))
102
dw1(2) = dcmplx(w1(3))
103
dw1(3) = dcmplx(w1(4))
104
dw2(0) = dcmplx(w2(1))
105
dw2(1) = dcmplx(w2(2))
106
dw2(2) = dcmplx(w2(3))
107
dw2(3) = dcmplx(w2(4))
108
dw3(0) = dcmplx(w3(1))
109
dw3(1) = dcmplx(w3(2))
110
dw3(2) = dcmplx(w3(3))
111
dw3(3) = dcmplx(w3(4))
114
p1(2) = dble(dimag(w1(6)))
115
p1(3) = dble(dimag(w1(5)))
118
p2(2) = dble(dimag(w2(6)))
119
p2(3) = dble(dimag(w2(5)))
122
p3(2) = dble(dimag(w3(6)))
123
p3(3) = dble(dimag(w3(5)))
124
q(0) = -(p1(0)+p2(0)+p3(0))
125
q(1) = -(p1(1)+p2(1)+p3(1))
126
q(2) = -(p1(2)+p2(2)+p3(2))
127
q(3) = -(p1(3)+p2(3)+p3(3))
128
q2 = q(0)**2 -(q(1)**2 +q(2)**2 +q(3)**2)
129
c Neil edited this file to allow implementation of 3-site model.
130
c Now, only gwwa matters and is the full coupling squared.
131
c dgwwa2 = dble(gwwa)**2
132
c dgwwz2 = dble(gwwz)**2
133
c dgw2 = dgwwa2+dgwwz2
137
c Benj's modif in order to have FR running
139
if(gwwa.eq.rZero) dgw2=gwwz
140
if(gwwz.eq.rZero) dgw2=gwwa
141
c End of Benj'S modif
150
if ( abs(jwww(5))+abs(jwww(6)).eq.rZero ) then
152
& ' helas-error : jwww in jwwwxx has zero momentum'
154
if ( wwidth.eq.rZero .and. q2.eq.mw2 ) then
156
& ' helas-error : jwww in jwwwxx is on wmass pole'
159
& real(q(0)),real(q(1)),real(q(2)),real(q(3))
161
& ' : abs(q)= ',sqrt(abs(real(q2)))
170
c Start of Claude'S modif (Removed minus Sign)
171
dw = rOne/dcmplx( q2-mw2, dmw*dww )
172
c End of Claude'S modif
174
c For the running width, use below instead of the above dw.
175
c dw = -rOne/dcmplx( q2-mw2 , max(dww*q2/dmw,rZero) )
177
w12=dw1(0)*dw2(0)-dw1(1)*dw2(1)-dw1(2)*dw2(2)-dw1(3)*dw2(3)
178
w32=dw3(0)*dw2(0)-dw3(1)*dw2(1)-dw3(2)*dw2(2)-dw3(3)*dw2(3)
180
w13=dw1(0)*dw3(0)-dw1(1)*dw3(1)-dw1(2)*dw3(2)-dw1(3)*dw3(3)
182
jj(0) = dgw2*( dw1(0)*w32 + dw3(0)*w12 - rTwo*dw2(0)*w13 )
183
jj(1) = dgw2*( dw1(1)*w32 + dw3(1)*w12 - rTwo*dw2(1)*w13 )
184
jj(2) = dgw2*( dw1(2)*w32 + dw3(2)*w12 - rTwo*dw2(2)*w13 )
185
jj(3) = dgw2*( dw1(3)*w32 + dw3(3)*w12 - rTwo*dw2(3)*w13 )
187
c Fabio's implementation of the fixed width
188
cm2=dcmplx( mw2, -dmw*dww )
189
c jq = (jj(0)*q(0)-jj(1)*q(1)-jj(2)*q(2)-jj(3)*q(3))/mw2
190
jq = (jj(0)*q(0)-jj(1)*q(1)-jj(2)*q(2)-jj(3)*q(3))/cm2
192
jwww(1) = dcmplx( (jj(0)-jq*q(0))*dw )
193
jwww(2) = dcmplx( (jj(1)-jq*q(1))*dw )
194
jwww(3) = dcmplx( (jj(2)-jq*q(2))*dw )
195
jwww(4) = dcmplx( (jj(3)-jq*q(3))*dw )