~maddevelopers/mg5amcnlo/3.0.2-alpha0

« back to all changes in this revision

Viewing changes to HELAS/jwwwnx.F

Added Template and HELAS into bzr

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
      subroutine jwwwnx(w1,w2,w3,gwwa,gwwz,wmass,wwidth , jwww)
 
2
c
 
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.
 
7
c
 
8
c input:
 
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
 
18
c
 
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.
 
27
c
 
28
c output:
 
29
c       complex jwww(6)        : W current             j^mu(w':w1,w2,w3)
 
30
c     
 
31
      implicit none
 
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
 
39
 
 
40
      double precision rZero, rOne, rTwo
 
41
      parameter( rZero = 0.0d0, rOne = 1.0d0, rTwo = 2.0d0 )
 
42
 
 
43
#ifdef HELAS_CHECK
 
44
      integer stdo
 
45
      parameter( stdo = 6 )
 
46
#endif
 
47
c
 
48
#ifdef HELAS_CHECK
 
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'
 
51
      endif
 
52
      if ( abs(w1(5))+abs(w1(6)).eq.rZero ) then
 
53
         write(stdo,*)
 
54
     &        ' helas-error : w1 in jwwwxx has zero momentum'
 
55
      endif
 
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'
 
58
      endif
 
59
      if ( abs(w2(5))+abs(w2(6)).eq.rZero ) then
 
60
         write(stdo,*)
 
61
     &        ' helas-error : w2 in jwwwxx has zero momentum'
 
62
      endif
 
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'
 
65
      endif
 
66
      if ( abs(w3(5))+abs(w3(6)).eq.rZero ) then
 
67
         write(stdo,*)
 
68
     &        ' helas-error : w3 in jwwwxx has zero momentum'
 
69
      endif 
 
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'
 
73
c      endif
 
74
c      if ( gwwz.eq.rZero ) then
 
75
c         write(stdo,*) ' helas-error : gwwz in jwwwxx is zero coupling'
 
76
c      endif
 
77
c      if ( gwwa.lt.rZero .or. gwwa.ge.gwwz ) then
 
78
c         write(stdo,*)
 
79
c     &  ' helas-warn  : gwwa/gwwz in jwwwxx are non-standard couplings'
 
80
c         write(stdo,*)
 
81
c     &  '             : gwwa = ',gwwa,'  gwwz = ',gwwz
 
82
c      endif
 
83
c     End Neil's edit
 
84
      if ( wmass.le.rZero ) then
 
85
         write(stdo,*) ' helas-error : wmass in jwwwxx is not positive'
 
86
         write(stdo,*) '             : wmass = ',wmass
 
87
      endif
 
88
      if ( wwidth.lt.rZero ) then
 
89
         write(stdo,*) ' helas-error : wwidth in jwwwxx is negative'
 
90
         write(stdo,*) '             : wwidth = ',wwidth
 
91
      endif
 
92
#endif
 
93
 
 
94
 
 
95
 
 
96
      
 
97
      jwww(5) = w1(5)+w2(5)+w3(5)
 
98
      jwww(6) = w1(6)+w2(6)+w3(6)
 
99
 
 
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))
 
112
      p1(0) = dble(      w1(5))
 
113
      p1(1) = dble(      w1(6))
 
114
      p1(2) = dble(dimag(w1(6)))
 
115
      p1(3) = dble(dimag(w1(5)))
 
116
      p2(0) = dble(      w2(5))
 
117
      p2(1) = dble(      w2(6))
 
118
      p2(2) = dble(dimag(w2(6)))
 
119
      p2(3) = dble(dimag(w2(5)))
 
120
      p3(0) = dble(      w3(5))
 
121
      p3(1) = dble(      w3(6))
 
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
 
134
c      dgw2 = gwwa
 
135
c     End Neil's edit
 
136
 
 
137
c Benj's modif in order to have FR running
 
138
      dgw2=rZero
 
139
      if(gwwa.eq.rZero) dgw2=gwwz
 
140
      if(gwwz.eq.rZero) dgw2=gwwa
 
141
c End of Benj'S modif
 
142
 
 
143
 
 
144
 
 
145
      dmw = dble(wmass)
 
146
      dww = dble(wwidth)
 
147
      mw2 = dmw**2
 
148
 
 
149
#ifdef HELAS_CHECK
 
150
      if ( abs(jwww(5))+abs(jwww(6)).eq.rZero ) then
 
151
         write(stdo,*)
 
152
     &        ' helas-error : jwww in jwwwxx has zero momentum'
 
153
      endif
 
154
      if ( wwidth.eq.rZero .and. q2.eq.mw2 ) then
 
155
         write(stdo,*)
 
156
     &        ' helas-error : jwww in jwwwxx is on wmass pole'
 
157
         write(stdo,*)
 
158
     &        '             : q     = ',
 
159
     &        real(q(0)),real(q(1)),real(q(2)),real(q(3))
 
160
         write(stdo,*)
 
161
     &        '             : abs(q)= ',sqrt(abs(real(q2)))
 
162
         jwww(1) = cZero
 
163
         jwww(2) = cZero
 
164
         jwww(3) = cZero
 
165
         jwww(4) = cZero
 
166
         return
 
167
      endif
 
168
#endif
 
169
 
 
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
 
173
 
 
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) )
 
176
 
 
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)
 
179
 
 
180
      w13=dw1(0)*dw3(0)-dw1(1)*dw3(1)-dw1(2)*dw3(2)-dw1(3)*dw3(3)
 
181
 
 
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 )
 
186
 
 
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
 
191
 
 
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 )
 
196
c
 
197
      return
 
198
      end