JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
mod_Graviton.F90
Go to the documentation of this file.
1  module modgraviton
2  use modparameters
3  use modmisc
4  implicit none
5  private
6 
7 
8 !----- notation for subroutines
9  public :: evalamp_gg_g_vv,evalamp_qqb_g_vv,evalamp_g_vv
10 
11  contains
12 
13 
14 !----- a subroutinefor gg -> G -> ZZ/WW
15 !----- all outgoing convention and the following momentum assignment
16 !----- 0 -> g(p1) + g(p2) + e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
17  subroutine evalamp_gg_g_vv(p,MY_IDUP,res)
18  implicit none
19  real(dp), intent(out) :: res
20  real(dp), intent(in) :: p(4,6)
21  integer, intent(in) :: my_idup(6:9)
22  complex(dp) :: a_vv(1:18), a0_vv(1:2)
23  integer :: i1,i2,i3,i4,vvmode,vvmode_swap
24  real(dp) :: prefactor
25  real(dp) :: intcolfac
26  integer :: ordering(1:4),ordering_swap(1:4)
27  logical :: dointerference
28 
29  if(isaquark(my_idup(6)) .and. isaquark(my_idup(8))) then
30  intcolfac=1.0_dp/3.0_dp
31  else
32  intcolfac=1.0_dp
33  endif
34 
35  call getdecay_vvmode_ordering(my_idup(6:9),vvmode,ordering,vvmode_swap,ordering_swap)
36 
37  if( vvmode.eq.zzmode ) then! Z decay
38  prefactor = 8d0*overallcouplvffsq**2
39  elseif( vvmode.eq.wwmode ) then ! W decay
40  prefactor = 8d0*overallcouplvffsq**2
41  elseif( vvmode.eq.zgmode ) then ! Z+photon "decay"
42  prefactor = 8d0*overallcouplvffsq ! Only single powers
43  elseif( vvmode.eq.ggmode ) then ! photon "decay"
44  prefactor = 8d0
45  else
46  prefactor = 0d0
47  endif
48 
49 
50  res = zero
51  a_vv(:) = 0d0
52  dointerference = includeinterference .and. ( &
53  ((vvmode.eq.zzmode) .and. (vvmode_swap.eq.zzmode)) &
54  )
55  if ( includevprime .and. .not.(vvmode.eq.zzmode .or. vvmode.eq.zgmode .or. vvmode.eq.wwmode) ) then
56  call error("Contact terms only for WW, ZZ or Zg!")
57  endif
58  do i1=1,2; do i2=1,2; do i3=1,2; do i4=1,2! sum over helicities
59  call calchelamp_gg(ordering,vvmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(1))
60  if( vvmode.eq.zzmode ) then
61  if( includegammastar ) then
62  call calchelamp_gg(ordering,zgsmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(3))
63  call calchelamp_gg(ordering,gszmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(5))
64  call calchelamp_gg(ordering,gsgsmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(7))
65  endif
66  if( includevprime ) then
67  call calchelamp_gg(ordering,zzpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(9))
68  call calchelamp_gg(ordering,zpzmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(11))
69  call calchelamp_gg(ordering,zpzpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(13))
70  endif
71  if( includegammastar .and. includevprime ) then
72  call calchelamp_gg(ordering,gszpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(15))
73  call calchelamp_gg(ordering,zpgsmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(17))
74  endif
75  elseif( vvmode.eq.zgmode ) then
76  if(includegammastar) then
77  call calchelamp_gg(ordering,gsgmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(3))
78  endif
79  if( includevprime ) then
80  call calchelamp_gg(ordering,zpgmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(5))
81  endif
82  elseif( vvmode.eq.wwmode .and. includevprime ) then
83  call calchelamp_gg(ordering,wwpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(9))
84  call calchelamp_gg(ordering,wpwmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(11))
85  call calchelamp_gg(ordering,wpwpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(13))
86  endif
87 
88  if( dointerference ) then
89  call calchelamp_gg(ordering_swap,vvmode_swap,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(2))
90  if( includegammastar ) then
91  call calchelamp_gg(ordering_swap,zgsmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(4))
92  call calchelamp_gg(ordering_swap,gszmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(6))
93  call calchelamp_gg(ordering_swap,gsgsmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(8))
94  endif
95  if( includevprime ) then
96  call calchelamp_gg(ordering_swap,zzpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(10))
97  call calchelamp_gg(ordering_swap,zpzmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(12))
98  call calchelamp_gg(ordering_swap,zpzpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(14))
99  endif
100  if( includegammastar .and. includevprime ) then
101  call calchelamp_gg(ordering_swap,gszpmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(16))
102  call calchelamp_gg(ordering_swap,zpgsmode,p(1:4,1:6),my_idup,i1,i2,i3,i4,a_vv(18))
103  endif
104  endif
105 
106  a0_vv(1) = a_vv(1)+a_vv(3)+a_vv(5)+a_vv(7)+a_vv(9)+a_vv(11)+a_vv(13)+a_vv(15)+a_vv(17) ! 3456 pieces
107  a0_vv(2) = a_vv(2)+a_vv(4)+a_vv(6)+a_vv(8)+a_vv(10)+a_vv(12)+a_vv(14)+a_vv(16)+a_vv(18) ! 5436 pieces
108  res = res + dreal(a0_vv(1)*dconjg(a0_vv(1)))
109  res = res + dreal(a0_vv(2)*dconjg(a0_vv(2)))
110  if( dointerference .and. (i3.eq.i4) ) then! interfere the 3456 with 5436 pieces
111  res = res - 2d0*intcolfac*dreal( a0_vv(1)*dconjg(a0_vv(2)) ) ! minus from Fermi statistics
112  endif
113  enddo; enddo; enddo; enddo
114 
115  res = res*prefactor
116  if( (vvmode.eq.zzmode) .and. dointerference ) res = res * symmfac
117 
118  end subroutine
119 
120  subroutine calchelamp_gg(ordering,VVMode,p,MY_IDUP,i1,i2,i3,i4,A)
121  implicit none
122  integer :: ordering(1:4),i1,i2,i3,i4,l1,l2,l3,l4,MY_IDUP(6:9),VVMode
123  real(dp) :: p(1:4,1:6)
124  complex(dp) :: propG
125  real(dp) :: pin(4,4)
126  complex(dp) :: A(1:1), sp(4,4)
127 
128 
129  l1=ordering(1)
130  l2=ordering(2)
131  l3=ordering(3)
132  l4=ordering(4)
133 
134  propg = one/dcmplx(2d0*scr(p(:,1),p(:,2)) - m_reso**2,m_reso*ga_reso)
135 
136  pin(1,:) = p(:,1)
137  pin(2,:) = p(:,2)
138  sp(1,:) = pol_mless2(dcmplx(p(:,1)),-3+2*i1,'in') ! gluon
139  sp(2,:) = pol_mless2(dcmplx(p(:,2)),-3+2*i2,'in') ! gluon
140  !sp(1,1:4)=pin(1,1:4);print *, "this checks IS gauge invariance"
141  !sp(2,1:4)=pin(2,1:4);print *, "this checks IS gauge invariance"
142  !careful: for gluon gauge invariance check the terms ~c3,c4 are needed because e1.q2 is not zero for e1-->q1
143  call getdecay_couplings_spinors_props( &
144  vvmode, &
145  (/my_idup(l1+3),my_idup(l2+3),my_idup(l3+3),my_idup(l4+3)/), &
146  (/p(:,l1),p(:,l2),p(:,l3),p(:,l4)/), &
147  -3+2*i3,-3+2*i4, &
148  sp(3:4,:),pin(3:4,:) &
149  )
150  call gggzzampl(vvmode,pin,sp,a(1))
151  a(1) = a(1)*propg
152 
153  end subroutine
154 
155  subroutine gggzzampl(VVMode,p,sp,res)
156  implicit none
157  integer, intent(in) :: VVMode
158  real(dp), intent(in) :: p(4,4)
159  complex(dp), intent(in) :: sp(4,4)
160  complex(dp), intent(out) :: res
161  complex(dp) :: e1_e2, e1_e3, e1_e4
162  complex(dp) :: e2_e3, e2_e4
163  complex(dp) :: e3_e4
164  complex(dp) :: q_q
165  complex(dp) :: q1_q2,q1_q3,q1_q4,q3_q3,q4_q4
166  complex(dp) :: q2_q3,q2_q4
167  complex(dp) :: q3_q4
168  complex(dp) :: q1_e3,q1_e4,q2_e3,q2_e4
169  complex(dp) :: e1_q3,e1_q4,e2_q3,e2_q4, q1_e1,q1_e2,q2_e1,q2_e2,q4_e1,q3_e1
170  complex(dp) :: e3_q4,e4_q3
171  complex(dp) :: q1(4),q2(4),q3(4),q4(4),q(4)
172  complex(dp) :: e1(4),e2(4),e3(4),e4(4)
173  complex(dp) :: xxx1,xxx2,xxx3,xxx4,yyy1,yyy2,yyy3,yyy4,yyy41,yyy42,yyy5,yyy6
174  complex(dp) :: yyy7,abr1
175  complex(dp) :: b_dyn(1:10)
176  real(dp) :: q34,MZ3,MZ4,MG
177  logical :: new
178  real(dp) :: rr_gam, rr
179 
180 
181  new = .true.
182 
183  b_dyn(:)=czero
184 
185  q1 = dcmplx(p(1,:),0d0)
186  q2 = dcmplx(p(2,:),0d0)
187  q3 = dcmplx(p(3,:),0d0)
188  q4 = dcmplx(p(4,:),0d0)
189 
190 
191  e1 = sp(1,:)
192  e2 = sp(2,:)
193  e3 = sp(3,:)
194  e4 = sp(4,:)
195 
196  q = -q1-q2
197  q_q = sc(q,q)
198  q3_q3 = sc(q3,q3)
199  q4_q4 = sc(q4,q4)
200 
201 
202  q1_q2 = sc(q1,q2)
203  q1_q3 = sc(q1,q3)
204  q1_q4 = sc(q1,q4)
205  q2_q3 = sc(q2,q3)
206  q2_q4 = sc(q2,q4)
207  q3_q4 = sc(q3,q4)
208 
209  e1_e2 = sc(e1,e2)
210  e1_e3 = sc(e1,e3)
211  e1_e4 = sc(e1,e4)
212  e2_e3 = sc(e2,e3)
213  e2_e4 = sc(e2,e4)
214  e3_e4 = sc(e3,e4)
215 
216 
217  q1_e3 = sc(q1,e3)
218  q1_e4 = sc(q1,e4)
219  q2_e3 = sc(q2,e3)
220  q2_e4 = sc(q2,e4)
221  e1_q3 = sc(e1,q3)
222  e1_q4 = sc(e1,q4)
223  e2_q3 = sc(e2,q3)
224  e2_q4 = sc(e2,q4)
225  e3_q4 = sc(e3,q4)
226  e4_q3 = sc(e4,q3)
227 
228  mz3=dsqrt(cdabs(q3_q3))
229  mz4=dsqrt(cdabs(q4_q4))
230  if( use_dynamic_mg ) then
231  mg = dsqrt(cdabs(q_q))
232  else
233  mg = m_reso
234  endif
235 
236 !---- define couplings
237  q34 = (mg**2-mz3**2-mz4**2)/2d0! = s = pV1.pV2 = (q_q-MZ3^2-MZ4^2)/2
238 
239  rr_gam = q_q/two/lambda**2! kappa for IS
240  xxx1 = (a1 + a2*rr_gam) ! those a's correspond to g's in eq.(7) for the IS
241  xxx2 = -a1/two + (a3+two*a4)*rr_gam
242  xxx3 = 4d0*a5*rr_gam * mg**2/q_q
243  ! for gluon gauge invariance check the terms ~c3,c4 are needed because e1.q2 is not zero for e1-->q1
244 
245  if (generate_bis) then
246  rr = q34/lambda**2! kappa for FS
247 
248  if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode) ) then! decay ZZ's or WW's
249  b_dyn(1)=b1
250  b_dyn(2)=b2
251  b_dyn(3)=b3
252  b_dyn(4)=b4
253  b_dyn(5)=b5
254  b_dyn(6)=b6
255  b_dyn(7)=b7
256  b_dyn(8)=b8
257  b_dyn(9)=b9
258  b_dyn(10)=b10
259  elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.gszmode) .OR. (vvmode.eq.zgsmode) ) then
260  b_dyn(1)=bzgs1
261  b_dyn(2)=bzgs2
262  b_dyn(3)=bzgs3
263  b_dyn(4)=bzgs4
264  b_dyn(8)=bzgs8
265  elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode) .or. (vvmode.eq.gsgmode) ) then
266  b_dyn(1)=bgsgs1
267  b_dyn(2)=bgsgs2
268  b_dyn(3)=bgsgs3
269  b_dyn(4)=bgsgs4
270  b_dyn(8)=bgsgs8
271  elseif( (vvmode.eq.zzpmode) .or. (vvmode.eq.wwpmode) .or. (vvmode.eq.zpzmode) .or. (vvmode.eq.wpwmode) ) then
272  b_dyn(1)=bzzp1
273  b_dyn(2)=bzzp2
274  b_dyn(3)=bzzp3
275  b_dyn(4)=bzzp4
276  b_dyn(5)=bzzp5
277  b_dyn(6)=bzzp6
278  b_dyn(7)=bzzp7
279  b_dyn(8)=bzzp8
280  b_dyn(9)=bzzp9
281  b_dyn(10)=bzzp10
282  elseif( (vvmode.eq.zpzpmode) .or. (vvmode.eq.wpwpmode) ) then
283  b_dyn(1)=bzpzp1
284  b_dyn(2)=bzpzp2
285  b_dyn(3)=bzpzp3
286  b_dyn(4)=bzpzp4
287  b_dyn(5)=bzpzp5
288  b_dyn(6)=bzpzp6
289  b_dyn(7)=bzpzp7
290  b_dyn(8)=bzpzp8
291  b_dyn(9)=bzpzp9
292  b_dyn(10)=bzpzp10
293  elseif( (vvmode.eq.zpgmode) .OR. (vvmode.eq.gszpmode) .OR. (vvmode.eq.zpgsmode) ) then
294  b_dyn(1)=bzpgs1
295  b_dyn(2)=bzpgs2
296  b_dyn(3)=bzpgs3
297  b_dyn(4)=bzpgs4
298  b_dyn(8)=bzpgs8
299  else
300  print *,"VVMode",vvmode,"not implemented"
301  endif
302 
303  yyy1 = q34*( b_dyn(1) + b_dyn(2)*rr*(one+mz3**2/q34)*(one+mz4**2/q34) ) + b_dyn(5)*m_v**2
304  yyy2 = -b_dyn(1)/two + b_dyn(3)*rr*(1d0-(mz3**2+mz4**2)/(2d0*q34)) + two*b_dyn(4)*rr + b_dyn(7)*rr*m_v**2/q34
305  yyy3 = (-b_dyn(2)/two - b_dyn(3)- two*b_dyn(4))*rr/q34
306  yyy41 = -b_dyn(1) - b_dyn(2)*(q34+mz3**2)/lambda**2 - b_dyn(3)*mz4**2/lambda**2 - 2d0*b_dyn(6)*m_v**2/lambda**2
307  yyy42 = -b_dyn(1) - b_dyn(2)*(q34+mz4**2)/lambda**2 - b_dyn(3)*mz3**2/lambda**2 - 2d0*b_dyn(6)*m_v**2/lambda**2
308  yyy5 = two*b_dyn(8)*rr*mg**2/q34
309  yyy6 = b_dyn(9) * m_v**2/lambda**2
310  yyy7 = b_dyn(10) * mg**2 * m_v**2/lambda**4
311 
312  else
313  yyy1 = q34*c1/2d0
314  yyy2 = c2
315  yyy3 = c3/mg**2
316  yyy41 = c41
317  yyy42 = c42
318  yyy5 = c5
319  yyy6 = czero
320  yyy7 = czero
321  if(vvmode.eq.zzmode .or. vvmode.eq.wwmode) then
322  yyy6 = c6
323  yyy7 = c7
324  endif
325  endif
326 
327  res = czero
328 
329 
330  if (new) then
331 
332 
333 ! res = &
334 ! + 8.D0*q1_e3*q1_e4*e1_e2*yyy1*xxx2 - 8.D0*q1_e3*q1_q3*e1_e2* &
335 ! e4_q3*yyy4*xxx2 + 4.D0*q1_e3*e1_e2*e4_q3*yyy1*xxx2 + 2.D0* &
336 ! q1_e3*e1_e2*e4_q3*MZ4**2*yyy4*xxx2 - 2.D0*q1_e3*e1_e2*e4_q3* &
337 ! MZ3**2*yyy4*xxx2 - 2.D0*q1_e3*e1_e2*e4_q3*MG**2*yyy4*xxx2 + 8.d0 &
338 ! *q1_e4*q1_q3*e1_e2*e3_q4*yyy4*xxx2 + 4.D0*q1_e4*e1_e2*e3_q4 &
339 ! *yyy1*xxx2 - 2.D0*q1_e4*e1_e2*e3_q4*MZ4**2*yyy4*xxx2 + 2.D0* &
340 ! q1_e4*e1_e2*e3_q4*MZ3**2*yyy4*xxx2 + 2.D0*q1_e4*e1_e2*e3_q4* &
341 ! MG**2*yyy4*xxx2 - 8.D0*q1_q3*e1_e2*e3_e4*MZ4**2*yyy2*xxx2 + 8.d0 &
342 ! *q1_q3*e1_e2*e3_e4*MZ3**2*yyy2*xxx2 + 8.D0*q1_q3*e1_e2* &
343 ! e3_e4*MG**2*yyy2*xxx2 - 8.D0*q1_q3*e1_e2*e3_q4*e4_q3*MZ4**2* &
344 ! yyy3*xxx2 + 8.D0*q1_q3*e1_e2*e3_q4*e4_q3*MZ3**2*yyy3*xxx2 + 8.d0 &
345 ! *q1_q3*e1_e2*e3_q4*e4_q3*MG**2*yyy3*xxx2 + 16.D0*q1_q3**2* &
346 ! e1_e2*e3_e4*yyy2*xxx2 + 16.D0*q1_q3**2*e1_e2*e3_q4*e4_q3*yyy3 &
347 ! *xxx2 + 2.D0/3.D0*e1_e2*e3_e4*MZ4**4*yyy2*xxx2 + 1.D0/3.D0* &
348 ! e1_e2*e3_e4*MZ4**4*yyy2*xxx1
349 ! res = res - 4.D0/3.D0*e1_e2*e3_e4*MZ3**2*MZ4**2*yyy2*xxx2 - 2.D0/ &
350 ! 3.D0*e1_e2*e3_e4*MZ3**2*MZ4**2*yyy2*xxx1 + 2.D0/3.D0*e1_e2* &
351 ! e3_e4*MZ3**4*yyy2*xxx2 + 1.D0/3.D0*e1_e2*e3_e4*MZ3**4*yyy2* &
352 ! xxx1 + 2.D0/3.D0*e1_e2*e3_e4*MG**2*yyy1*xxx2 - 2.D0/3.D0* &
353 ! e1_e2*e3_e4*MG**2*yyy1*xxx1 - 4.D0/3.D0*e1_e2*e3_e4*MG**2* &
354 ! MZ4**2*yyy2*xxx2 - 2.D0/3.D0*e1_e2*e3_e4*MG**2*MZ4**2*yyy2* &
355 ! xxx1 + 8.D0/3.D0*e1_e2*e3_e4*MG**2*MZ3**2*yyy2*xxx2 - 2.D0/3.D0 &
356 ! *e1_e2*e3_e4*MG**2*MZ3**2*yyy2*xxx1 + 2.D0/3.D0*e1_e2*e3_e4* &
357 ! MG**4*yyy2*xxx2 + 1.D0/3.D0*e1_e2*e3_e4*MG**4*yyy2*xxx1 + 4.D0 &
358 ! /3.D0*e1_e2*e3_q4*e4_q3*yyy1*xxx2 + 2.D0/3.D0*e1_e2*e3_q4* &
359 ! e4_q3*yyy1*xxx1 + 2.D0/3.D0*e1_e2*e3_q4*e4_q3*MZ4**4*yyy3* &
360 ! xxx2 + 1.D0/3.D0*e1_e2*e3_q4*e4_q3*MZ4**4*yyy3*xxx1 - 4.D0/3.D0 &
361 ! *e1_e2*e3_q4*e4_q3*MZ3**2*MZ4**2*yyy3*xxx2 - 2.D0/3.D0*e1_e2 &
362 ! *e3_q4*e4_q3*MZ3**2*MZ4**2*yyy3*xxx1 + 2.D0/3.D0*e1_e2*e3_q4* &
363 ! e4_q3*MZ3**4*yyy3*xxx2 + 1.D0/3.D0*e1_e2*e3_q4*e4_q3*MZ3**4* &
364 ! yyy3*xxx1
365 ! res = res + 2.D0/3.D0*e1_e2*e3_q4*e4_q3*MG**2*yyy4*xxx2 - 2.D0/3.D0 &
366 ! *e1_e2*e3_q4*e4_q3*MG**2*yyy4*xxx1 - 4.D0/3.D0*e1_e2*e3_q4* &
367 ! e4_q3*MG**2*MZ4**2*yyy3*xxx2 - 2.D0/3.D0*e1_e2*e3_q4*e4_q3* &
368 ! MG**2*MZ4**2*yyy3*xxx1 + 8.D0/3.D0*e1_e2*e3_q4*e4_q3*MG**2* &
369 ! MZ3**2*yyy3*xxx2 - 2.D0/3.D0*e1_e2*e3_q4*e4_q3*MG**2*MZ3**2* &
370 ! yyy3*xxx1 + 2.D0/3.D0*e1_e2*e3_q4*e4_q3*MG**4*yyy3*xxx2 + 1.D0 &
371 ! /3.D0*e1_e2*e3_q4*e4_q3*MG**4*yyy3*xxx1 + e1_e3*e2_e4*MG**2* &
372 ! yyy1*xxx1 - e1_e3*e2_q3*e4_q3*MG**2*yyy4*xxx1 + e1_e4*e2_e3* &
373 ! MG**2*yyy1*xxx1 + e1_e4*e2_q3*e3_q4*MG**2*yyy4*xxx1 - e1_q3* &
374 ! e2_e3*e4_q3*MG**2*yyy4*xxx1 + e1_q3*e2_e4*e3_q4*MG**2*yyy4* &
375 ! xxx1 + 4.D0*e1_q3*e2_q3*e3_e4*MG**2*yyy2*xxx1 + 4.D0*e1_q3* &
376 ! e2_q3*e3_q4*e4_q3*MG**2*yyy3*xxx1 + 8.D0*et1(q1,q2,e1,e2)* &
377 ! q1_e3*q1_e4*MG**(-2)*yyy1*xxx3 - 8.D0*et1(q1,q2,e1,e2)*q1_e3* &
378 ! q1_q3*e4_q3*MG**(-2)*yyy4*xxx3 + 4.D0*et1(q1,q2,e1,e2)*q1_e3* &
379 ! e4_q3*MG**(-2)*yyy1*xxx3 + 2.D0*et1(q1,q2,e1,e2)*q1_e3*e4_q3* &
380 ! MG**(-2)*MZ4**2*yyy4*xxx3
381 ! res = res - 2.D0*et1(q1,q2,e1,e2)*q1_e3*e4_q3*MG**(-2)*MZ3**2* &
382 ! yyy4*xxx3 - 2.D0*et1(q1,q2,e1,e2)*q1_e3*e4_q3*yyy4*xxx3 + 8.D0* &
383 ! et1(q1,q2,e1,e2)*q1_e4*q1_q3*e3_q4*MG**(-2)*yyy4*xxx3 + 4.D0* &
384 ! et1(q1,q2,e1,e2)*q1_e4*e3_q4*MG**(-2)*yyy1*xxx3 - 2.D0*et1(q1 &
385 ! ,q2,e1,e2)*q1_e4*e3_q4*MG**(-2)*MZ4**2*yyy4*xxx3 + 2.D0*et1( &
386 ! q1,q2,e1,e2)*q1_e4*e3_q4*MG**(-2)*MZ3**2*yyy4*xxx3 + 2.D0* &
387 ! et1(q1,q2,e1,e2)*q1_e4*e3_q4*yyy4*xxx3 - 8.D0*et1(q1,q2,e1,e2 &
388 ! )*q1_q3*e3_e4*MG**(-2)*MZ4**2*yyy2*xxx3 + 8.D0*et1(q1,q2,e1, &
389 ! e2)*q1_q3*e3_e4*MG**(-2)*MZ3**2*yyy2*xxx3 + 8.D0*et1(q1,q2,e1 &
390 ! ,e2)*q1_q3*e3_e4*yyy2*xxx3 - 8.D0*et1(q1,q2,e1,e2)*q1_q3* &
391 ! e3_q4*e4_q3*MG**(-2)*MZ4**2*yyy3*xxx3 + 8.D0*et1(q1,q2,e1,e2) &
392 ! *q1_q3*e3_q4*e4_q3*MG**(-2)*MZ3**2*yyy3*xxx3 + 8.D0*et1(q1,q2 &
393 ! ,e1,e2)*q1_q3*e3_q4*e4_q3*yyy3*xxx3 + 16.D0*et1(q1,q2,e1,e2)* &
394 ! q1_q3**2*e3_e4*MG**(-2)*yyy2*xxx3 + 16.D0*et1(q1,q2,e1,e2)* &
395 ! q1_q3**2*e3_q4*e4_q3*MG**(-2)*yyy3*xxx3 + 2.D0/3.D0*et1(q1,q2 &
396 ! ,e1,e2)*e3_e4*MG**(-2)*MZ4**4*yyy2*xxx3
397 ! res = res - 4.D0/3.D0*et1(q1,q2,e1,e2)*e3_e4*MG**(-2)*MZ3**2* &
398 ! MZ4**2*yyy2*xxx3 + 2.D0/3.D0*et1(q1,q2,e1,e2)*e3_e4*MG**(-2)* &
399 ! MZ3**4*yyy2*xxx3 + 2.D0/3.D0*et1(q1,q2,e1,e2)*e3_e4*yyy1*xxx3 &
400 ! - 4.D0/3.D0*et1(q1,q2,e1,e2)*e3_e4*MZ4**2*yyy2*xxx3 + 8.D0/3.D0 &
401 ! *et1(q1,q2,e1,e2)*e3_e4*MZ3**2*yyy2*xxx3 + 2.D0/3.D0*et1(q1 &
402 ! ,q2,e1,e2)*e3_e4*MG**2*yyy2*xxx3 + 4.D0/3.D0*et1(q1,q2,e1,e2) &
403 ! *e3_q4*e4_q3*MG**(-2)*yyy1*xxx3 + 2.D0/3.D0*et1(q1,q2,e1,e2)* &
404 ! e3_q4*e4_q3*MG**(-2)*MZ4**4*yyy3*xxx3 - 4.D0/3.D0*et1(q1,q2, &
405 ! e1,e2)*e3_q4*e4_q3*MG**(-2)*MZ3**2*MZ4**2*yyy3*xxx3 + 2.D0/3.D0 &
406 ! *et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**(-2)*MZ3**4*yyy3*xxx3 + 2.D0 &
407 ! /3.D0*et1(q1,q2,e1,e2)*e3_q4*e4_q3*yyy4*xxx3 - 4.D0/3.D0* &
408 ! et1(q1,q2,e1,e2)*e3_q4*e4_q3*MZ4**2*yyy3*xxx3 + 8.D0/3.D0* &
409 ! et1(q1,q2,e1,e2)*e3_q4*e4_q3*MZ3**2*yyy3*xxx3 + 2.D0/3.D0* &
410 ! et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**2*yyy3*xxx3 - et1(q1,q2,e1, &
411 ! e2)*et1(q1,q,e3,e4)*MG**(-2)*MZ4**2*yyy6*xxx3 + et1(q1,q2,e1, &
412 ! e2)*et1(q1,q,e3,e4)*MG**(-2)*MZ3**2*yyy6*xxx3
413 ! res = res + et1(q1,q2,e1,e2)*et1(q1,q,e3,e4)*yyy6*xxx3 + 4.D0* &
414 ! et1(q1,q2,e1,e2)*et1(q1,q,e3,e4)*q1_q3*MG**(-2)*yyy6*xxx3 - 4.0d0 &
415 ! *et1(q1,q2,e1,e2)*et1(q1,q,e3,q3)*q1_q3*e4_q3*MG**(-4)*yyy7 &
416 ! *xxx3 + et1(q1,q2,e1,e2)*et1(q1,q,e3,q3)*e4_q3*MG**(-4)* &
417 ! MZ4**2*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e3,q3)*e4_q3* &
418 ! MG**(-4)*MZ3**2*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e3,q3)* &
419 ! e4_q3*MG**(-2)*yyy7*xxx3 + 4.D0*et1(q1,q2,e1,e2)*et1(q1,q,e3, &
420 ! q4)*q1_q3*e4_q3*MG**(-4)*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1, &
421 ! q,e3,q4)*e4_q3*MG**(-4)*MZ4**2*yyy7*xxx3 + et1(q1,q2,e1,e2)* &
422 ! et1(q1,q,e3,q4)*e4_q3*MG**(-4)*MZ3**2*yyy7*xxx3 + et1(q1,q2, &
423 ! e1,e2)*et1(q1,q,e3,q4)*e4_q3*MG**(-2)*yyy7*xxx3 - 4.D0*et1(q1 &
424 ! ,q2,e1,e2)*et1(q1,q,e4,q3)*q1_q3*e3_q4*MG**(-4)*yyy7*xxx3 + &
425 ! et1(q1,q2,e1,e2)*et1(q1,q,e4,q3)*e3_q4*MG**(-4)*MZ4**2*yyy7* &
426 ! xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e4,q3)*e3_q4*MG**(-4)*MZ3**2 &
427 ! *yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e4,q3)*e3_q4*MG**(-2)* &
428 ! yyy7*xxx3
429 ! res = res + 4.D0*et1(q1,q2,e1,e2)*et1(q1,q,e4,q4)*q1_q3*e3_q4* &
430 ! MG**(-4)*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e4,q4)*e3_q4* &
431 ! MG**(-4)*MZ4**2*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q1,q,e4,q4)* &
432 ! e3_q4*MG**(-4)*MZ3**2*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q1,q, &
433 ! e4,q4)*e3_q4*MG**(-2)*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q2,q, &
434 ! e3,e4)*MG**(-2)*MZ4**2*yyy6*xxx3 - et1(q1,q2,e1,e2)*et1(q2,q, &
435 ! e3,e4)*MG**(-2)*MZ3**2*yyy6*xxx3 - et1(q1,q2,e1,e2)*et1(q2,q, &
436 ! e3,e4)*yyy6*xxx3 - 4.D0*et1(q1,q2,e1,e2)*et1(q2,q,e3,e4)* &
437 ! q1_q3*MG**(-2)*yyy6*xxx3 + 4.D0*et1(q1,q2,e1,e2)*et1(q2,q,e3, &
438 ! q3)*q1_q3*e4_q3*MG**(-4)*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q2, &
439 ! q,e3,q3)*e4_q3*MG**(-4)*MZ4**2*yyy7*xxx3 + et1(q1,q2,e1,e2)* &
440 ! et1(q2,q,e3,q3)*e4_q3*MG**(-4)*MZ3**2*yyy7*xxx3 + et1(q1,q2, &
441 ! e1,e2)*et1(q2,q,e3,q3)*e4_q3*MG**(-2)*yyy7*xxx3 - 4.D0*et1(q1 &
442 ! ,q2,e1,e2)*et1(q2,q,e3,q4)*q1_q3*e4_q3*MG**(-4)*yyy7*xxx3 + &
443 ! et1(q1,q2,e1,e2)*et1(q2,q,e3,q4)*e4_q3*MG**(-4)*MZ4**2*yyy7* &
444 ! xxx3
445 ! res = res - et1(q1,q2,e1,e2)*et1(q2,q,e3,q4)*e4_q3*MG**(-4)* &
446 ! MZ3**2*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q2,q,e3,q4)*e4_q3* &
447 ! MG**(-2)*yyy7*xxx3 + 4.D0*et1(q1,q2,e1,e2)*et1(q2,q,e4,q3)* &
448 ! q1_q3*e3_q4*MG**(-4)*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q2,q,e4 &
449 ! ,q3)*e3_q4*MG**(-4)*MZ4**2*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1( &
450 ! q2,q,e4,q3)*e3_q4*MG**(-4)*MZ3**2*yyy7*xxx3 + et1(q1,q2,e1,e2 &
451 ! )*et1(q2,q,e4,q3)*e3_q4*MG**(-2)*yyy7*xxx3 - 4.D0*et1(q1,q2, &
452 ! e1,e2)*et1(q2,q,e4,q4)*q1_q3*e3_q4*MG**(-4)*yyy7*xxx3 + et1( &
453 ! q1,q2,e1,e2)*et1(q2,q,e4,q4)*e3_q4*MG**(-4)*MZ4**2*yyy7*xxx3 &
454 ! - et1(q1,q2,e1,e2)*et1(q2,q,e4,q4)*e3_q4*MG**(-4)*MZ3**2* &
455 ! yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q2,q,e4,q4)*e3_q4*MG**(-2)* &
456 ! yyy7*xxx3 - 1.D0/3.D0*et1(q1,q2,e1,e2)*et1(q,e3,e4,q3)*yyy6* &
457 ! xxx3 + 1.D0/3.D0*et1(q1,q2,e1,e2)*et1(q,e3,e4,q4)*yyy6*xxx3 &
458 ! + 2.D0/3.D0*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*MG**(-4)* &
459 ! MZ4**4*yyy5*xxx3 - 4.D0/3.D0*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4 &
460 ! )*MG**(-4)*MZ3**2*MZ4**2*yyy5*xxx3
461 ! res = res + 2.D0/3.D0*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*MG**(-4)* &
462 ! MZ3**4*yyy5*xxx3 - 4.D0/3.D0*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)* &
463 ! MG**(-2)*MZ4**2*yyy5*xxx3 + 8.D0/3.D0*et1(q1,q2,e1,e2)*et1(e3 &
464 ! ,e4,q3,q4)*MG**(-2)*MZ3**2*yyy5*xxx3 + 2.D0/3.D0*et1(q1,q2,e1 &
465 ! ,e2)*et1(e3,e4,q3,q4)*yyy5*xxx3 - 8.D0*et1(q1,q2,e1,e2)*et1( &
466 ! e3,e4,q3,q4)*q1_q3*MG**(-4)*MZ4**2*yyy5*xxx3 + 8.D0*et1(q1,q2 &
467 ! ,e1,e2)*et1(e3,e4,q3,q4)*q1_q3*MG**(-4)*MZ3**2*yyy5*xxx3 + 8.D0 &
468 ! *et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*q1_q3*MG**(-2)*yyy5*xxx3 &
469 ! + 16.D0*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*q1_q3**2*MG**(-4)* &
470 ! yyy5*xxx3 + 4.D0*et1(q1,q,e3,e4)*q1_q3*e1_e2*yyy6*xxx2 - et1( &
471 ! q1,q,e3,e4)*e1_e2*MZ4**2*yyy6*xxx2 + et1(q1,q,e3,e4)*e1_e2* &
472 ! MZ3**2*yyy6*xxx2 + et1(q1,q,e3,e4)*e1_e2*MG**2*yyy6*xxx2 - 4.D0 &
473 ! *et1(q1,q,e3,q3)*q1_q3*e1_e2*e4_q3*MG**(-2)*yyy7*xxx2 + et1( &
474 ! q1,q,e3,q3)*e1_e2*e4_q3*MG**(-2)*MZ4**2*yyy7*xxx2 - et1(q1,q, &
475 ! e3,q3)*e1_e2*e4_q3*MG**(-2)*MZ3**2*yyy7*xxx2 - et1(q1,q,e3,q3 &
476 ! )*e1_e2*e4_q3*yyy7*xxx2
477 ! res = res + 4.D0*et1(q1,q,e3,q4)*q1_q3*e1_e2*e4_q3*MG**(-2)*yyy7* &
478 ! xxx2 - et1(q1,q,e3,q4)*e1_e2*e4_q3*MG**(-2)*MZ4**2*yyy7*xxx2 + &
479 ! et1(q1,q,e3,q4)*e1_e2*e4_q3*MG**(-2)*MZ3**2*yyy7*xxx2 + et1( &
480 ! q1,q,e3,q4)*e1_e2*e4_q3*yyy7*xxx2 - 4.D0*et1(q1,q,e4,q3)* &
481 ! q1_q3*e1_e2*e3_q4*MG**(-2)*yyy7*xxx2 + et1(q1,q,e4,q3)*e1_e2* &
482 ! e3_q4*MG**(-2)*MZ4**2*yyy7*xxx2 - et1(q1,q,e4,q3)*e1_e2*e3_q4 &
483 ! *MG**(-2)*MZ3**2*yyy7*xxx2 - et1(q1,q,e4,q3)*e1_e2*e3_q4*yyy7 &
484 ! *xxx2 + 4.D0*et1(q1,q,e4,q4)*q1_q3*e1_e2*e3_q4*MG**(-2)*yyy7* &
485 ! xxx2 - et1(q1,q,e4,q4)*e1_e2*e3_q4*MG**(-2)*MZ4**2*yyy7*xxx2 &
486 ! + et1(q1,q,e4,q4)*e1_e2*e3_q4*MG**(-2)*MZ3**2*yyy7*xxx2 + &
487 ! et1(q1,q,e4,q4)*e1_e2*e3_q4*yyy7*xxx2 - 4.D0*et1(q2,q,e3,e4)* &
488 ! q1_q3*e1_e2*yyy6*xxx2 + et1(q2,q,e3,e4)*e1_e2*MZ4**2*yyy6* &
489 ! xxx2 - et1(q2,q,e3,e4)*e1_e2*MZ3**2*yyy6*xxx2 - et1(q2,q,e3, &
490 ! e4)*e1_e2*MG**2*yyy6*xxx2 + 4.D0*et1(q2,q,e3,q3)*q1_q3*e1_e2* &
491 ! e4_q3*MG**(-2)*yyy7*xxx2 - et1(q2,q,e3,q3)*e1_e2*e4_q3* &
492 ! MG**(-2)*MZ4**2*yyy7*xxx2
493 ! res = res + et1(q2,q,e3,q3)*e1_e2*e4_q3*MG**(-2)*MZ3**2*yyy7*xxx2 &
494 ! + et1(q2,q,e3,q3)*e1_e2*e4_q3*yyy7*xxx2 - 4.D0*et1(q2,q,e3, &
495 ! q4)*q1_q3*e1_e2*e4_q3*MG**(-2)*yyy7*xxx2 + et1(q2,q,e3,q4)* &
496 ! e1_e2*e4_q3*MG**(-2)*MZ4**2*yyy7*xxx2 - et1(q2,q,e3,q4)*e1_e2 &
497 ! *e4_q3*MG**(-2)*MZ3**2*yyy7*xxx2 - et1(q2,q,e3,q4)*e1_e2* &
498 ! e4_q3*yyy7*xxx2 + 4.D0*et1(q2,q,e4,q3)*q1_q3*e1_e2*e3_q4* &
499 ! MG**(-2)*yyy7*xxx2 - et1(q2,q,e4,q3)*e1_e2*e3_q4*MG**(-2)* &
500 ! MZ4**2*yyy7*xxx2 + et1(q2,q,e4,q3)*e1_e2*e3_q4*MG**(-2)* &
501 ! MZ3**2*yyy7*xxx2 + et1(q2,q,e4,q3)*e1_e2*e3_q4*yyy7*xxx2 - 4.D0 &
502 ! *et1(q2,q,e4,q4)*q1_q3*e1_e2*e3_q4*MG**(-2)*yyy7*xxx2 + et1( &
503 ! q2,q,e4,q4)*e1_e2*e3_q4*MG**(-2)*MZ4**2*yyy7*xxx2 - et1(q2,q, &
504 ! e4,q4)*e1_e2*e3_q4*MG**(-2)*MZ3**2*yyy7*xxx2 - et1(q2,q,e4,q4 &
505 ! )*e1_e2*e3_q4*yyy7*xxx2 - et1(q,e1,e3,e4)*e2_q3*MG**2*yyy6* &
506 ! xxx1 + et1(q,e1,e3,q3)*e2_q3*e4_q3*yyy7*xxx1 - et1(q,e1,e3,q4 &
507 ! )*e2_q3*e4_q3*yyy7*xxx1 + et1(q,e1,e4,q3)*e2_q3*e3_q4*yyy7* &
508 ! xxx1
509 ! res = res - et1(q,e1,e4,q4)*e2_q3*e3_q4*yyy7*xxx1 - et1(q,e2,e3, &
510 ! e4)*e1_q3*MG**2*yyy6*xxx1 + et1(q,e2,e3,q3)*e1_q3*e4_q3*yyy7* &
511 ! xxx1 - et1(q,e2,e3,q4)*e1_q3*e4_q3*yyy7*xxx1 + et1(q,e2,e4,q3 &
512 ! )*e1_q3*e3_q4*yyy7*xxx1 - et1(q,e2,e4,q4)*e1_q3*e3_q4*yyy7* &
513 ! xxx1 - 1.D0/3.D0*et1(q,e3,e4,q3)*e1_e2*MG**2*yyy6*xxx2 + 1.D0/ &
514 ! 3.D0*et1(q,e3,e4,q3)*e1_e2*MG**2*yyy6*xxx1 + 1.D0/3.D0*et1(q, &
515 ! e3,e4,q4)*e1_e2*MG**2*yyy6*xxx2 - 1.D0/3.D0*et1(q,e3,e4,q4)* &
516 ! e1_e2*MG**2*yyy6*xxx1 - 8.D0*et1(e3,e4,q3,q4)*q1_q3*e1_e2* &
517 ! MG**(-2)*MZ4**2*yyy5*xxx2 + 8.D0*et1(e3,e4,q3,q4)*q1_q3*e1_e2 &
518 ! *MG**(-2)*MZ3**2*yyy5*xxx2 + 8.D0*et1(e3,e4,q3,q4)*q1_q3* &
519 ! e1_e2*yyy5*xxx2 + 16.D0*et1(e3,e4,q3,q4)*q1_q3**2*e1_e2* &
520 ! MG**(-2)*yyy5*xxx2 + 2.D0/3.D0*et1(e3,e4,q3,q4)*e1_e2* &
521 ! MG**(-2)*MZ4**4*yyy5*xxx2 + 1.D0/3.D0*et1(e3,e4,q3,q4)*e1_e2* &
522 ! MG**(-2)*MZ4**4*yyy5*xxx1 - 4.D0/3.D0*et1(e3,e4,q3,q4)*e1_e2* &
523 ! MG**(-2)*MZ3**2*MZ4**2*yyy5*xxx2 - 2.D0/3.D0*et1(e3,e4,q3,q4) &
524 ! *e1_e2*MG**(-2)*MZ3**2*MZ4**2*yyy5*xxx1
525 ! res = res + 2.D0/3.D0*et1(e3,e4,q3,q4)*e1_e2*MG**(-2)*MZ3**4*yyy5 &
526 ! *xxx2 + 1.D0/3.D0*et1(e3,e4,q3,q4)*e1_e2*MG**(-2)*MZ3**4*yyy5* &
527 ! xxx1 - 4.D0/3.D0*et1(e3,e4,q3,q4)*e1_e2*MZ4**2*yyy5*xxx2 - 2.D0 &
528 ! /3.D0*et1(e3,e4,q3,q4)*e1_e2*MZ4**2*yyy5*xxx1 + 8.D0/3.D0* &
529 ! et1(e3,e4,q3,q4)*e1_e2*MZ3**2*yyy5*xxx2 - 2.D0/3.D0*et1(e3,e4 &
530 ! ,q3,q4)*e1_e2*MZ3**2*yyy5*xxx1 + 2.D0/3.D0*et1(e3,e4,q3,q4)* &
531 ! e1_e2*MG**2*yyy5*xxx2 + 1.D0/3.D0*et1(e3,e4,q3,q4)*e1_e2* &
532 ! MG**2*yyy5*xxx1 + 4.D0*et1(e3,e4,q3,q4)*e1_q3*e2_q3*yyy5*xxx1
533 !
534 ! print *, "old res GG",res
535 
536 
537 ! this is the new code that includes couplings yyy41 and yyy42 instead of yyy4
538 ! res =&
539 ! & + 8.*q1_e3*q1_e4*e1_e2*yyy1*xxx2 - 8.*q1_e3*q1_q3*e1_e2*e4_q3*&
540 ! & yyy41*xxx2 + 4.*q1_e3*e1_e2*e4_q3*yyy1*xxx2 + 2.*q1_e3*e1_e2*&
541 ! & e4_q3*MZ4**2*yyy41*xxx2 - 2.*q1_e3*e1_e2*e4_q3*MZ3**2*yyy41*&
542 ! & xxx2 - 2.*q1_e3*e1_e2*e4_q3*MG**2*yyy41*xxx2 + 8.*q1_e4*q1_q3&
543 ! & *e1_e2*e3_q4*yyy42*xxx2 + 4.*q1_e4*e1_e2*e3_q4*yyy1*xxx2 - 2.&
544 ! & *q1_e4*e1_e2*e3_q4*MZ4**2*yyy42*xxx2 + 2.*q1_e4*e1_e2*e3_q4*&
545 ! & MZ3**2*yyy42*xxx2 + 2.*q1_e4*e1_e2*e3_q4*MG**2*yyy42*xxx2 - 8.&
546 ! & *q1_q3*e1_e2*e3_e4*MZ4**2*yyy2*xxx2 + 8.*q1_q3*e1_e2*e3_e4*&
547 ! & MZ3**2*yyy2*xxx2 + 8.*q1_q3*e1_e2*e3_e4*MG**2*yyy2*xxx2 + 4.*&
548 ! & q1_q3*e1_e2*e3_q4*e4_q3*yyy42*xxx2 - 4.*q1_q3*e1_e2*e3_q4*&
549 ! & e4_q3*yyy41*xxx2 - 8.*q1_q3*e1_e2*e3_q4*e4_q3*MZ4**2*yyy3*&
550 ! & xxx2 + 8.*q1_q3*e1_e2*e3_q4*e4_q3*MZ3**2*yyy3*xxx2 + 8.*q1_q3&
551 ! & *e1_e2*e3_q4*e4_q3*MG**2*yyy3*xxx2 + 16.*q1_q3**2*e1_e2*e3_e4&
552 ! & *yyy2*xxx2 + 16.*q1_q3**2*e1_e2*e3_q4*e4_q3*yyy3*xxx2 + 2./3.&
553 ! & *e1_e2*e3_e4*MZ4**4*yyy2*xxx2
554 ! res = res + 1./3.*e1_e2*e3_e4*MZ4**4*yyy2*xxx1 - 4./3.*e1_e2*&
555 ! & e3_e4*MZ3**2*MZ4**2*yyy2*xxx2 - 2./3.*e1_e2*e3_e4*MZ3**2*&
556 ! & MZ4**2*yyy2*xxx1 + 2./3.*e1_e2*e3_e4*MZ3**4*yyy2*xxx2 + 1./3.&
557 ! & *e1_e2*e3_e4*MZ3**4*yyy2*xxx1 + 2./3.*e1_e2*e3_e4*MG**2*yyy1*&
558 ! & xxx2 - 2./3.*e1_e2*e3_e4*MG**2*yyy1*xxx1 - 4./3.*e1_e2*e3_e4*&
559 ! & MG**2*MZ4**2*yyy2*xxx2 - 2./3.*e1_e2*e3_e4*MG**2*MZ4**2*yyy2*&
560 ! & xxx1 + 8./3.*e1_e2*e3_e4*MG**2*MZ3**2*yyy2*xxx2 - 2./3.*e1_e2&
561 ! & *e3_e4*MG**2*MZ3**2*yyy2*xxx1 + 2./3.*e1_e2*e3_e4*MG**4*yyy2*&
562 ! & xxx2 + 1./3.*e1_e2*e3_e4*MG**4*yyy2*xxx1 + 4./3.*e1_e2*e3_q4*&
563 ! & e4_q3*yyy1*xxx2 + 2./3.*e1_e2*e3_q4*e4_q3*yyy1*xxx1 - 2./3.*&
564 ! & e1_e2*e3_q4*e4_q3*MZ4**2*yyy42*xxx2 - 1./3.*e1_e2*e3_q4*e4_q3&
565 ! & *MZ4**2*yyy42*xxx1 + 2./3.*e1_e2*e3_q4*e4_q3*MZ4**2*yyy41*&
566 ! & xxx2 + 1./3.*e1_e2*e3_q4*e4_q3*MZ4**2*yyy41*xxx1 + 2./3.*&
567 ! & e1_e2*e3_q4*e4_q3*MZ4**4*yyy3*xxx2 + 1./3.*e1_e2*e3_q4*e4_q3*&
568 ! & MZ4**4*yyy3*xxx1 + 2./3.*e1_e2*e3_q4*e4_q3*MZ3**2*yyy42*xxx2&
569 ! & + 1./3.*e1_e2*e3_q4*e4_q3*MZ3**2*yyy42*xxx1
570 ! res = res - 2./3.*e1_e2*e3_q4*e4_q3*MZ3**2*yyy41*xxx2 - 1./3.*&
571 ! & e1_e2*e3_q4*e4_q3*MZ3**2*yyy41*xxx1 - 4./3.*e1_e2*e3_q4*e4_q3&
572 ! & *MZ3**2*MZ4**2*yyy3*xxx2 - 2./3.*e1_e2*e3_q4*e4_q3*MZ3**2*&
573 ! & MZ4**2*yyy3*xxx1 + 2./3.*e1_e2*e3_q4*e4_q3*MZ3**4*yyy3*xxx2&
574 ! & + 1./3.*e1_e2*e3_q4*e4_q3*MZ3**4*yyy3*xxx1 + 4./3.*e1_e2*&
575 ! & e3_q4*e4_q3*MG**2*yyy42*xxx2 - 1./3.*e1_e2*e3_q4*e4_q3*MG**2*&
576 ! & yyy42*xxx1 - 2./3.*e1_e2*e3_q4*e4_q3*MG**2*yyy41*xxx2 - 1./3.&
577 ! & *e1_e2*e3_q4*e4_q3*MG**2*yyy41*xxx1 - 4./3.*e1_e2*e3_q4*e4_q3&
578 ! & *MG**2*MZ4**2*yyy3*xxx2 - 2./3.*e1_e2*e3_q4*e4_q3*MG**2*&
579 ! & MZ4**2*yyy3*xxx1 + 8./3.*e1_e2*e3_q4*e4_q3*MG**2*MZ3**2*yyy3*&
580 ! & xxx2 - 2./3.*e1_e2*e3_q4*e4_q3*MG**2*MZ3**2*yyy3*xxx1 + 2./3.&
581 ! & *e1_e2*e3_q4*e4_q3*MG**4*yyy3*xxx2 + 1./3.*e1_e2*e3_q4*e4_q3*&
582 ! & MG**4*yyy3*xxx1 + e1_e3*e2_e4*MG**2*yyy1*xxx1 - e1_e3*e2_q3*&
583 ! & e4_q3*MG**2*yyy41*xxx1 + e1_e4*e2_e3*MG**2*yyy1*xxx1 + e1_e4*&
584 ! & e2_q3*e3_q4*MG**2*yyy42*xxx1 - e1_q3*e2_e3*e4_q3*MG**2*yyy41*&
585 ! & xxx1
586 ! res = res + e1_q3*e2_e4*e3_q4*MG**2*yyy42*xxx1 + 4.*e1_q3*e2_q3*&
587 ! & e3_e4*MG**2*yyy2*xxx1 + 4.*e1_q3*e2_q3*e3_q4*e4_q3*MG**2*yyy3&
588 ! & *xxx1 + 8.*et1(q1,q2,e1,e2)*q1_e3*q1_e4*MG**(-2)*yyy1*xxx3 - &
589 ! & 8.*et1(q1,q2,e1,e2)*q1_e3*q1_q3*e4_q3*MG**(-2)*yyy41*xxx3 + 4.&
590 ! & *et1(q1,q2,e1,e2)*q1_e3*e4_q3*MG**(-2)*yyy1*xxx3 + 2.*et1(q1,&
591 ! & q2,e1,e2)*q1_e3*e4_q3*MG**(-2)*MZ4**2*yyy41*xxx3 - 2.*et1(q1,&
592 ! & q2,e1,e2)*q1_e3*e4_q3*MG**(-2)*MZ3**2*yyy41*xxx3 - 2.*et1(q1,&
593 ! & q2,e1,e2)*q1_e3*e4_q3*yyy41*xxx3 + 8.*et1(q1,q2,e1,e2)*q1_e4*&
594 ! & q1_q3*e3_q4*MG**(-2)*yyy42*xxx3 + 4.*et1(q1,q2,e1,e2)*q1_e4*&
595 ! & e3_q4*MG**(-2)*yyy1*xxx3 - 2.*et1(q1,q2,e1,e2)*q1_e4*e3_q4*&
596 ! & MG**(-2)*MZ4**2*yyy42*xxx3 + 2.*et1(q1,q2,e1,e2)*q1_e4*e3_q4*&
597 ! & MG**(-2)*MZ3**2*yyy42*xxx3 + 2.*et1(q1,q2,e1,e2)*q1_e4*e3_q4*&
598 ! & yyy42*xxx3 - 8.*et1(q1,q2,e1,e2)*q1_q3*e3_e4*MG**(-2)*MZ4**2*&
599 ! & yyy2*xxx3 + 8.*et1(q1,q2,e1,e2)*q1_q3*e3_e4*MG**(-2)*MZ3**2*&
600 ! & yyy2*xxx3 + 8.*et1(q1,q2,e1,e2)*q1_q3*e3_e4*yyy2*xxx3 + 4.*&
601 ! & et1(q1,q2,e1,e2)*q1_q3*e3_q4*e4_q3*MG**(-2)*yyy42*xxx3
602 ! res = res - 4.*et1(q1,q2,e1,e2)*q1_q3*e3_q4*e4_q3*MG**(-2)*yyy41*&
603 ! & xxx3 - 8.*et1(q1,q2,e1,e2)*q1_q3*e3_q4*e4_q3*MG**(-2)*MZ4**2*&
604 ! & yyy3*xxx3 + 8.*et1(q1,q2,e1,e2)*q1_q3*e3_q4*e4_q3*MG**(-2)*&
605 ! & MZ3**2*yyy3*xxx3 + 8.*et1(q1,q2,e1,e2)*q1_q3*e3_q4*e4_q3*yyy3&
606 ! & *xxx3 + 16.*et1(q1,q2,e1,e2)*q1_q3**2*e3_e4*MG**(-2)*yyy2*&
607 ! & xxx3 + 16.*et1(q1,q2,e1,e2)*q1_q3**2*e3_q4*e4_q3*MG**(-2)*&
608 ! & yyy3*xxx3 + 2./3.*et1(q1,q2,e1,e2)*e3_e4*MG**(-2)*MZ4**4*yyy2&
609 ! & *xxx3 - 4./3.*et1(q1,q2,e1,e2)*e3_e4*MG**(-2)*MZ3**2*MZ4**2*&
610 ! & yyy2*xxx3 + 2./3.*et1(q1,q2,e1,e2)*e3_e4*MG**(-2)*MZ3**4*yyy2&
611 ! & *xxx3 + 2./3.*et1(q1,q2,e1,e2)*e3_e4*yyy1*xxx3 - 4./3.*et1(q1&
612 ! & ,q2,e1,e2)*e3_e4*MZ4**2*yyy2*xxx3 + 8./3.*et1(q1,q2,e1,e2)*&
613 ! & e3_e4*MZ3**2*yyy2*xxx3 + 2./3.*et1(q1,q2,e1,e2)*e3_e4*MG**2*&
614 ! & yyy2*xxx3 + 4./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**(-2)*yyy1*&
615 ! & xxx3 - 2./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**(-2)*MZ4**2*&
616 ! & yyy42*xxx3 + 2./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**(-2)*&
617 ! & MZ4**2*yyy41*xxx3
618 ! res = res + 2./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**(-2)*MZ4**4*&
619 ! & yyy3*xxx3 + 2./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**(-2)*MZ3**2*&
620 ! & yyy42*xxx3 - 2./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**(-2)*&
621 ! & MZ3**2*yyy41*xxx3 - 4./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*&
622 ! & MG**(-2)*MZ3**2*MZ4**2*yyy3*xxx3 + 2./3.*et1(q1,q2,e1,e2)*&
623 ! & e3_q4*e4_q3*MG**(-2)*MZ3**4*yyy3*xxx3 + 4./3.*et1(q1,q2,e1,e2&
624 ! & )*e3_q4*e4_q3*yyy42*xxx3 - 2./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3&
625 ! & *yyy41*xxx3 - 4./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MZ4**2*yyy3*&
626 ! & xxx3 + 8./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MZ3**2*yyy3*xxx3 + &
627 ! & 2./3.*et1(q1,q2,e1,e2)*e3_q4*e4_q3*MG**2*yyy3*xxx3 - et1(q1,&
628 ! & q2,e1,e2)*et1(q1,q,e3,e4)*MG**(-2)*MZ4**2*yyy6*xxx3 + et1(q1,&
629 ! & q2,e1,e2)*et1(q1,q,e3,e4)*MG**(-2)*MZ3**2*yyy6*xxx3 + et1(q1,&
630 ! & q2,e1,e2)*et1(q1,q,e3,e4)*yyy6*xxx3 + 4.*et1(q1,q2,e1,e2)*&
631 ! & et1(q1,q,e3,e4)*q1_q3*MG**(-2)*yyy6*xxx3 - 4.*et1(q1,q2,e1,e2&
632 ! & )*et1(q1,q,e3,q3)*q1_q3*e4_q3*MG**(-4)*yyy7*xxx3 + et1(q1,q2,&
633 ! & e1,e2)*et1(q1,q,e3,q3)*e4_q3*MG**(-4)*MZ4**2*yyy7*xxx3
634 ! res = res - et1(q1,q2,e1,e2)*et1(q1,q,e3,q3)*e4_q3*MG**(-4)*&
635 ! & MZ3**2*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e3,q3)*e4_q3*&
636 ! & MG**(-2)*yyy7*xxx3 + 4.*et1(q1,q2,e1,e2)*et1(q1,q,e3,q4)*&
637 ! & q1_q3*e4_q3*MG**(-4)*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e3&
638 ! & ,q4)*e4_q3*MG**(-4)*MZ4**2*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(&
639 ! & q1,q,e3,q4)*e4_q3*MG**(-4)*MZ3**2*yyy7*xxx3 + et1(q1,q2,e1,e2&
640 ! & )*et1(q1,q,e3,q4)*e4_q3*MG**(-2)*yyy7*xxx3 - 4.*et1(q1,q2,e1,&
641 ! & e2)*et1(q1,q,e4,q3)*q1_q3*e3_q4*MG**(-4)*yyy7*xxx3 + et1(q1,&
642 ! & q2,e1,e2)*et1(q1,q,e4,q3)*e3_q4*MG**(-4)*MZ4**2*yyy7*xxx3 - &
643 ! & et1(q1,q2,e1,e2)*et1(q1,q,e4,q3)*e3_q4*MG**(-4)*MZ3**2*yyy7*&
644 ! & xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e4,q3)*e3_q4*MG**(-2)*yyy7*&
645 ! & xxx3 + 4.*et1(q1,q2,e1,e2)*et1(q1,q,e4,q4)*q1_q3*e3_q4*&
646 ! & MG**(-4)*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q1,q,e4,q4)*e3_q4*&
647 ! & MG**(-4)*MZ4**2*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q1,q,e4,q4)*&
648 ! & e3_q4*MG**(-4)*MZ3**2*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q1,q,&
649 ! & e4,q4)*e3_q4*MG**(-2)*yyy7*xxx3
650 ! res = res + et1(q1,q2,e1,e2)*et1(q2,q,e3,e4)*MG**(-2)*MZ4**2*yyy6&
651 ! & *xxx3 - et1(q1,q2,e1,e2)*et1(q2,q,e3,e4)*MG**(-2)*MZ3**2*yyy6*&
652 ! & xxx3 - et1(q1,q2,e1,e2)*et1(q2,q,e3,e4)*yyy6*xxx3 - 4.*et1(q1&
653 ! & ,q2,e1,e2)*et1(q2,q,e3,e4)*q1_q3*MG**(-2)*yyy6*xxx3 + 4.*et1(&
654 ! & q1,q2,e1,e2)*et1(q2,q,e3,q3)*q1_q3*e4_q3*MG**(-4)*yyy7*xxx3&
655 ! & - et1(q1,q2,e1,e2)*et1(q2,q,e3,q3)*e4_q3*MG**(-4)*MZ4**2*&
656 ! & yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q2,q,e3,q3)*e4_q3*MG**(-4)*&
657 ! & MZ3**2*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q2,q,e3,q3)*e4_q3*&
658 ! & MG**(-2)*yyy7*xxx3 - 4.*et1(q1,q2,e1,e2)*et1(q2,q,e3,q4)*&
659 ! & q1_q3*e4_q3*MG**(-4)*yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q2,q,e3&
660 ! & ,q4)*e4_q3*MG**(-4)*MZ4**2*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(&
661 ! & q2,q,e3,q4)*e4_q3*MG**(-4)*MZ3**2*yyy7*xxx3 - et1(q1,q2,e1,e2&
662 ! & )*et1(q2,q,e3,q4)*e4_q3*MG**(-2)*yyy7*xxx3 + 4.*et1(q1,q2,e1,&
663 ! & e2)*et1(q2,q,e4,q3)*q1_q3*e3_q4*MG**(-4)*yyy7*xxx3 - et1(q1,&
664 ! & q2,e1,e2)*et1(q2,q,e4,q3)*e3_q4*MG**(-4)*MZ4**2*yyy7*xxx3 + &
665 ! & et1(q1,q2,e1,e2)*et1(q2,q,e4,q3)*e3_q4*MG**(-4)*MZ3**2*yyy7*&
666 ! & xxx3
667 ! res = res + et1(q1,q2,e1,e2)*et1(q2,q,e4,q3)*e3_q4*MG**(-2)*yyy7*&
668 ! & xxx3 - 4.*et1(q1,q2,e1,e2)*et1(q2,q,e4,q4)*q1_q3*e3_q4*MG**(-4)*&
669 ! & yyy7*xxx3 + et1(q1,q2,e1,e2)*et1(q2,q,e4,q4)*e3_q4*MG**(-4)*&
670 ! & MZ4**2*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q2,q,e4,q4)*e3_q4*&
671 ! & MG**(-4)*MZ3**2*yyy7*xxx3 - et1(q1,q2,e1,e2)*et1(q2,q,e4,q4)*&
672 ! & e3_q4*MG**(-2)*yyy7*xxx3 - 1./3.*et1(q1,q2,e1,e2)*et1(q,e3,e4&
673 ! & ,q3)*yyy6*xxx3 + 1./3.*et1(q1,q2,e1,e2)*et1(q,e3,e4,q4)*yyy6*&
674 ! & xxx3 + 2./3.*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*MG**(-4)*&
675 ! & MZ4**4*yyy5*xxx3 - 4./3.*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*&
676 ! & MG**(-4)*MZ3**2*MZ4**2*yyy5*xxx3 + 2./3.*et1(q1,q2,e1,e2)*&
677 ! & et1(e3,e4,q3,q4)*MG**(-4)*MZ3**4*yyy5*xxx3 - 4./3.*et1(q1,q2,&
678 ! & e1,e2)*et1(e3,e4,q3,q4)*MG**(-2)*MZ4**2*yyy5*xxx3 + 8./3.*&
679 ! & et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*MG**(-2)*MZ3**2*yyy5*xxx3&
680 ! & + 2./3.*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*yyy5*xxx3 - 8.*&
681 ! & et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*q1_q3*MG**(-4)*MZ4**2*yyy5*&
682 ! & xxx3
683 ! res = res + 8.*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*q1_q3*MG**(-4)*&
684 ! & MZ3**2*yyy5*xxx3 + 8.*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*q1_q3*&
685 ! & MG**(-2)*yyy5*xxx3 + 16.*et1(q1,q2,e1,e2)*et1(e3,e4,q3,q4)*&
686 ! & q1_q3**2*MG**(-4)*yyy5*xxx3 + 4.*et1(q1,q,e3,e4)*q1_q3*e1_e2*&
687 ! & yyy6*xxx2 - et1(q1,q,e3,e4)*e1_e2*MZ4**2*yyy6*xxx2 + et1(q1,q&
688 ! & ,e3,e4)*e1_e2*MZ3**2*yyy6*xxx2 + et1(q1,q,e3,e4)*e1_e2*MG**2*&
689 ! & yyy6*xxx2 - 4.*et1(q1,q,e3,q3)*q1_q3*e1_e2*e4_q3*MG**(-2)*&
690 ! & yyy7*xxx2 + et1(q1,q,e3,q3)*e1_e2*e4_q3*MG**(-2)*MZ4**2*yyy7*&
691 ! & xxx2 - et1(q1,q,e3,q3)*e1_e2*e4_q3*MG**(-2)*MZ3**2*yyy7*xxx2&
692 ! & - et1(q1,q,e3,q3)*e1_e2*e4_q3*yyy7*xxx2 + 4.*et1(q1,q,e3,q4)&
693 ! & *q1_q3*e1_e2*e4_q3*MG**(-2)*yyy7*xxx2 - et1(q1,q,e3,q4)*e1_e2&
694 ! & *e4_q3*MG**(-2)*MZ4**2*yyy7*xxx2 + et1(q1,q,e3,q4)*e1_e2*&
695 ! & e4_q3*MG**(-2)*MZ3**2*yyy7*xxx2 + et1(q1,q,e3,q4)*e1_e2*e4_q3&
696 ! & *yyy7*xxx2 - 4.*et1(q1,q,e4,q3)*q1_q3*e1_e2*e3_q4*MG**(-2)*&
697 ! & yyy7*xxx2 + et1(q1,q,e4,q3)*e1_e2*e3_q4*MG**(-2)*MZ4**2*yyy7*&
698 ! & xxx2
699 ! res = res - et1(q1,q,e4,q3)*e1_e2*e3_q4*MG**(-2)*MZ3**2*yyy7*xxx2&
700 ! & - et1(q1,q,e4,q3)*e1_e2*e3_q4*yyy7*xxx2 + 4.*et1(q1,q,e4,q4)&
701 ! & *q1_q3*e1_e2*e3_q4*MG**(-2)*yyy7*xxx2 - et1(q1,q,e4,q4)*e1_e2&
702 ! & *e3_q4*MG**(-2)*MZ4**2*yyy7*xxx2 + et1(q1,q,e4,q4)*e1_e2*&
703 ! & e3_q4*MG**(-2)*MZ3**2*yyy7*xxx2 + et1(q1,q,e4,q4)*e1_e2*e3_q4&
704 ! & *yyy7*xxx2 - 4.*et1(q2,q,e3,e4)*q1_q3*e1_e2*yyy6*xxx2 + et1(&
705 ! & q2,q,e3,e4)*e1_e2*MZ4**2*yyy6*xxx2 - et1(q2,q,e3,e4)*e1_e2*&
706 ! & MZ3**2*yyy6*xxx2 - et1(q2,q,e3,e4)*e1_e2*MG**2*yyy6*xxx2 + 4.&
707 ! & *et1(q2,q,e3,q3)*q1_q3*e1_e2*e4_q3*MG**(-2)*yyy7*xxx2 - et1(&
708 ! & q2,q,e3,q3)*e1_e2*e4_q3*MG**(-2)*MZ4**2*yyy7*xxx2 + et1(q2,q,&
709 ! & e3,q3)*e1_e2*e4_q3*MG**(-2)*MZ3**2*yyy7*xxx2 + et1(q2,q,e3,q3&
710 ! & )*e1_e2*e4_q3*yyy7*xxx2 - 4.*et1(q2,q,e3,q4)*q1_q3*e1_e2*&
711 ! & e4_q3*MG**(-2)*yyy7*xxx2 + et1(q2,q,e3,q4)*e1_e2*e4_q3*&
712 ! & MG**(-2)*MZ4**2*yyy7*xxx2 - et1(q2,q,e3,q4)*e1_e2*e4_q3*&
713 ! & MG**(-2)*MZ3**2*yyy7*xxx2 - et1(q2,q,e3,q4)*e1_e2*e4_q3*yyy7*&
714 ! & xxx2
715 ! res = res + 4.*et1(q2,q,e4,q3)*q1_q3*e1_e2*e3_q4*MG**(-2)*yyy7*&
716 ! & xxx2 - et1(q2,q,e4,q3)*e1_e2*e3_q4*MG**(-2)*MZ4**2*yyy7*xxx2 + &
717 ! & et1(q2,q,e4,q3)*e1_e2*e3_q4*MG**(-2)*MZ3**2*yyy7*xxx2 + et1(&
718 ! & q2,q,e4,q3)*e1_e2*e3_q4*yyy7*xxx2 - 4.*et1(q2,q,e4,q4)*q1_q3*&
719 ! & e1_e2*e3_q4*MG**(-2)*yyy7*xxx2 + et1(q2,q,e4,q4)*e1_e2*e3_q4*&
720 ! & MG**(-2)*MZ4**2*yyy7*xxx2 - et1(q2,q,e4,q4)*e1_e2*e3_q4*&
721 ! & MG**(-2)*MZ3**2*yyy7*xxx2 - et1(q2,q,e4,q4)*e1_e2*e3_q4*yyy7*&
722 ! & xxx2 - et1(q,e1,e3,e4)*e2_q3*MG**2*yyy6*xxx1 + et1(q,e1,e3,q3&
723 ! & )*e2_q3*e4_q3*yyy7*xxx1 - et1(q,e1,e3,q4)*e2_q3*e4_q3*yyy7*&
724 ! & xxx1 + et1(q,e1,e4,q3)*e2_q3*e3_q4*yyy7*xxx1 - et1(q,e1,e4,q4&
725 ! & )*e2_q3*e3_q4*yyy7*xxx1 - et1(q,e2,e3,e4)*e1_q3*MG**2*yyy6*&
726 ! & xxx1 + et1(q,e2,e3,q3)*e1_q3*e4_q3*yyy7*xxx1 - et1(q,e2,e3,q4&
727 ! & )*e1_q3*e4_q3*yyy7*xxx1 + et1(q,e2,e4,q3)*e1_q3*e3_q4*yyy7*&
728 ! & xxx1 - et1(q,e2,e4,q4)*e1_q3*e3_q4*yyy7*xxx1 - 1./3.*et1(q,e3&
729 ! & ,e4,q3)*e1_e2*MG**2*yyy6*xxx2 + 1./3.*et1(q,e3,e4,q3)*e1_e2*&
730 ! & MG**2*yyy6*xxx1
731 ! res = res + 1./3.*et1(q,e3,e4,q4)*e1_e2*MG**2*yyy6*xxx2 - 1./3.*&
732 ! & et1(q,e3,e4,q4)*e1_e2*MG**2*yyy6*xxx1 - 8.*et1(e3,e4,q3,q4)*&
733 ! & q1_q3*e1_e2*MG**(-2)*MZ4**2*yyy5*xxx2 + 8.*et1(e3,e4,q3,q4)*&
734 ! & q1_q3*e1_e2*MG**(-2)*MZ3**2*yyy5*xxx2 + 8.*et1(e3,e4,q3,q4)*&
735 ! & q1_q3*e1_e2*yyy5*xxx2 + 16.*et1(e3,e4,q3,q4)*q1_q3**2*e1_e2*&
736 ! & MG**(-2)*yyy5*xxx2 + 2./3.*et1(e3,e4,q3,q4)*e1_e2*MG**(-2)*&
737 ! & MZ4**4*yyy5*xxx2 + 1./3.*et1(e3,e4,q3,q4)*e1_e2*MG**(-2)*&
738 ! & MZ4**4*yyy5*xxx1 - 4./3.*et1(e3,e4,q3,q4)*e1_e2*MG**(-2)*&
739 ! & MZ3**2*MZ4**2*yyy5*xxx2 - 2./3.*et1(e3,e4,q3,q4)*e1_e2*&
740 ! & MG**(-2)*MZ3**2*MZ4**2*yyy5*xxx1 + 2./3.*et1(e3,e4,q3,q4)*&
741 ! & e1_e2*MG**(-2)*MZ3**4*yyy5*xxx2 + 1./3.*et1(e3,e4,q3,q4)*&
742 ! & e1_e2*MG**(-2)*MZ3**4*yyy5*xxx1 - 4./3.*et1(e3,e4,q3,q4)*&
743 ! & e1_e2*MZ4**2*yyy5*xxx2 - 2./3.*et1(e3,e4,q3,q4)*e1_e2*MZ4**2*&
744 ! & yyy5*xxx1 + 8./3.*et1(e3,e4,q3,q4)*e1_e2*MZ3**2*yyy5*xxx2 - 2.&
745 ! & /3.*et1(e3,e4,q3,q4)*e1_e2*MZ3**2*yyy5*xxx1 + 2./3.*et1(e3,e4&
746 ! & ,q3,q4)*e1_e2*MG**2*yyy5*xxx2
747 ! res = res + 1./3.*et1(e3,e4,q3,q4)*e1_e2*MG**2*yyy5*xxx1 + 4.*&
748 ! & et1(e3,e4,q3,q4)*e1_q3*e2_q3*yyy5*xxx1
749 
750 
751 ! print *, "new ",res
752 
753 
754 ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! SAME AS ABOVE BUT SHORTER ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
755 
756  abr1 = mg**2 + mz3**2 - mz4**2 + 4*q1_q3
757 
758  res = mg**2*xxx1*(e1_e4*(yyy1*e2_e3 + yyy42*e2_q3*e3_q4) + e1_e3*(yyy1*e2_e4 - yyy41*e2_q3*e4_q3) + &
759  e1_q3*(yyy42*e2_e4*e3_q4 - yyy41*e2_e3*e4_q3 + 4*e2_q3*(yyy2*e3_e4 + yyy3*e3_q4*e4_q3))) + &
760  e1_e2*((xxx1*((mg**4*yyy2 + (mz3**2 - mz4**2)**2*yyy2 - 2*mg**2*(yyy1 + (mz3**2 + mz4**2)*yyy2))*e3_e4 + &
761  (2*yyy1 + mg**4*yyy3 + (mz3 - mz4)*(mz3 + mz4)*((mz3 - mz4)*(mz3 + mz4)*yyy3 - yyy41 + yyy42) - mg**2*(2*(mz3**2 + mz4**2)*yyy3 + yyy41 + yyy42))* &
762  e3_q4*e4_q3))/3d0 + (2*xxx2*(3*q1_e3*((2*yyy1 - abr1*yyy41)*e4_q3 + 4*yyy1*q1_e4) + &
763  e3_e4*(abr1**2*yyy2 + abr1*(yyy1 + 2*mz3**2*yyy2) - (mz3 - mz4)*(mz3 + mz4)*(yyy1 + 2*mz3**2*yyy2) - &
764  4*q1_q3*(yyy1 - abr1*yyy2 + 2*mz3**2*yyy2 + 8*yyy2*q1_q3) + 24*yyy2*q1_q3**2) + &
765  e3_q4*(3*(2*yyy1 + abr1*yyy42)*q1_e4 + e4_q3* &
766  (2*yyy1 + abr1**2*yyy3 - (mz3 - mz4)*(mz3 + mz4)*(2*mz3**2*yyy3 + yyy42) + abr1*(2*mz3**2*yyy3 - yyy41 + 2*yyy42 ) - &
767  2*q1_q3*(-2*abr1*yyy3 + 4*mz3**2*yyy3 + yyy41 + yyy42 + 16*yyy3*q1_q3) + 24*yyy3*q1_q3**2))))/3d0) - &
768  mg**2*xxx1*yyy6*e2_q3*et1(q,e1,e3,e4) + xxx1*yyy7*e2_q3*e4_q3*et1(q,e1,e3,q3) - xxx1*yyy7*e2_q3*e4_q3*et1(q,e1,e3,q4) + &
769  xxx1*yyy7*e2_q3*e3_q4*et1(q,e1,e4,q3) - xxx1*yyy7*e2_q3*e3_q4*et1(q,e1,e4,q4) - mg**2*xxx1*yyy6*e1_q3*et1(q,e2,e3,e4) + &
770  xxx1*yyy7*e1_q3*e4_q3*et1(q,e2,e3,q3) - xxx1*yyy7*e1_q3*e4_q3*et1(q,e2,e3,q4) + xxx1*yyy7*e1_q3*e3_q4*et1(q,e2,e4,q3) - &
771  xxx1*yyy7*e1_q3*e3_q4*et1(q,e2,e4,q4) + et1(q,e3,e4,q3)* &
772  (((mg**2*xxx1*yyy6)/3d0 - (mg**2*xxx2*yyy6)/3d0)*e1_e2 - (xxx3*yyy6*et1(q1,q2,e1,e2))/3d0) + &
773  et1(q,e3,e4,q4)*((-(mg**2*xxx1*yyy6)/3d0 + (mg**2*xxx2*yyy6)/3d0)*e1_e2 + (xxx3*yyy6*et1(q1,q2,e1,e2))/3d0) + &
774  et1(q1,q,e3,e4)*(abr1*xxx2*yyy6*e1_e2 + (abr1*xxx3*yyy6*et1(q1,q2,e1,e2))/mg**2) + &
775  et1(q1,q,e4,q3)*(-((abr1*xxx2*yyy7*e1_e2*e3_q4)/mg**2) - (abr1*xxx3*yyy7*e3_q4*et1(q1,q2,e1,e2))/mg**4) + &
776  et1(q1,q,e4,q4)*((abr1*xxx2*yyy7*e1_e2*e3_q4)/mg**2 + (abr1*xxx3*yyy7*e3_q4*et1(q1,q2,e1,e2))/mg**4) + &
777  et1(q1,q,e3,q3)*(-((abr1*xxx2*yyy7*e1_e2*e4_q3)/mg**2) - (abr1*xxx3*yyy7*e4_q3*et1(q1,q2,e1,e2))/mg**4) + &
778  et1(q1,q,e3,q4)*((abr1*xxx2*yyy7*e1_e2*e4_q3)/mg**2 + (abr1*xxx3*yyy7*e4_q3*et1(q1,q2,e1,e2))/mg**4) + &
779  et1(e3,e4,q3,q4)*(4*xxx1*yyy5*e1_q3*e2_q3 + e1_e2* &
780  (((mg - mz3 - mz4)*(mg + mz3 - mz4)*(mg - mz3 + mz4)*(mg + mz3 + mz4)*xxx1*yyy5)/(3d0*mg**2) + &
781  (2*xxx2*yyy5*(abr1**2 + 2*abr1*mz3**2 - 2*mz3**4 + 2*mz3**2*mz4**2 + 4*(abr1 - 2*mz3**2 - 8*q1_q3)*q1_q3 + 24*q1_q3**2))/(3d0*mg**2)) + &
782  (2*xxx3*yyy5*(abr1**2 + 2*abr1*mz3**2 - 2*mz3**4 + 2*mz3**2*mz4**2 + 4*(abr1 - 2*mz3**2 - 8*q1_q3)*q1_q3 + 24*q1_q3**2)*et1(q1,q2,e1,e2))/ &
783  (3d0*mg**4)) - abr1*xxx2*yyy6*e1_e2*et1(q2,q,e3,e4) + (abr1*xxx2*yyy7*e1_e2*e4_q3*et1(q2,q,e3,q3))/mg**2 - &
784  (abr1*xxx2*yyy7*e1_e2*e4_q3*et1(q2,q,e3,q4))/mg**2 + (abr1*xxx2*yyy7*e1_e2*e3_q4*et1(q2,q,e4,q3))/mg**2 - &
785  (abr1*xxx2*yyy7*e1_e2*e3_q4*et1(q2,q,e4,q4))/mg**2 + &
786  et1(q1,q2,e1,e2)*((2*xxx3*(3*q1_e3*((2*yyy1 - abr1*yyy41)*e4_q3 + 4*yyy1*q1_e4) + &
787  e3_e4*(abr1**2*yyy2 + abr1*(yyy1 + 2*mz3**2*yyy2) - (mz3 - mz4)*(mz3 + mz4)*(yyy1 + 2*mz3**2*yyy2) - &
788  4*q1_q3*(yyy1 - abr1*yyy2 + 2*mz3**2*yyy2 + 8*yyy2*q1_q3) + 24*yyy2*q1_q3**2) + &
789  e3_q4*(3*(2*yyy1 + abr1*yyy42)*q1_e4 + e4_q3* &
790  (2*yyy1 + abr1**2*yyy3 - (mz3 - mz4)*(mz3 + mz4)*(2*mz3**2*yyy3 + yyy42) + abr1*(2*mz3**2*yyy3 - yyy41 + 2*yyy42) - &
791  2*q1_q3*(-2*abr1*yyy3 + 4*mz3**2*yyy3 + yyy41 + yyy42 + 16*yyy3*q1_q3) + 24*yyy3*q1_q3**2))))/(3d0*mg**2) - &
792  (abr1*xxx3*yyy6*et1(q2,q,e3,e4))/mg**2 + (abr1*xxx3*yyy7*e4_q3*et1(q2,q,e3,q3))/mg**4 - (abr1*xxx3*yyy7*e4_q3*et1(q2,q,e3,q4))/mg**4 + &
793  (abr1*xxx3*yyy7*e3_q4*et1(q2,q,e4,q3))/mg**4 - (abr1*xxx3*yyy7*e3_q4*et1(q2,q,e4,q4))/mg**4)
794 
795 
796 ! print *, "newer",res
797 ! pause
798 
799 ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
800 
801  else
802 
803 
804 
805  res = q1_e3*q2_e4*e1_e2 * ( m_reso**2 )
806  res = res + q1_e3*q2_q4*e1_e2*e4_q3 * ( - 2.0 )
807  res = res + q1_e4*q2_e3*e1_e2 * ( m_reso**2 )
808  res = res + q1_e4*q2_q3*e1_e2*e3_q4 * ( - 2.0 )
809  res = res + q1_q3*q2_e4*e1_e2*e3_q4 * ( - 2.0 )
810  res = res + q1_q3*q2_q4*e1_e2*e3_e4 * ( 2.0 )
811  res = res + q1_q3*e1_e2*e3_e4 * ( 1./2.*m_reso**2 )
812  res = res + q1_q3*e1_e2*e3_q4*e4_q3 * ( - 1. )
813  res = res + q1_q4*q2_e3*e1_e2*e4_q3 * ( - 2. )
814  res = res + q1_q4*q2_q3*e1_e2*e3_e4 * ( 2 )
815  res = res + q1_q4*e1_e2*e3_e4 * ( 1./2.*m_reso**2 )
816  res = res + q1_q4*e1_e2*e3_q4*e4_q3 * ( - 1. )
817  res = res + q2_q3*e1_e2*e3_e4 * ( 1./2.*m_reso**2 )
818  res = res + q2_q3*e1_e2*e3_q4*e4_q3 * ( - 1 )
819  res = res + q2_q4*e1_e2*e3_e4 * ( 1./2.*m_reso**2 )
820  res = res + q2_q4*e1_e2*e3_q4*e4_q3 * ( - 1. )
821  res = res + e1_e2*e3_e4 * ( m_reso**2*m_v**2 - 1./2.*m_reso**4 )
822  res = res + e1_e2*e3_q4*e4_q3 * ( m_reso**2 )
823  res = res + e1_e3*e2_e4 * ( 1./2.*m_reso**4 )
824  res = res + e1_e3*e2_q4*e4_q3 * ( - m_reso**2 )
825  res = res + e1_e4*e2_e3 * ( 1./2.*m_reso**4 )
826  res = res + e1_e4*e2_q3*e3_q4 * ( - m_reso**2 )
827  res = res + e1_q3*e2_e4*e3_q4 * ( - m_reso**2 )
828  res = res + e1_q3*e2_q4*e3_e4 * ( m_reso**2 )
829  res = res + e1_q4*e2_e3*e4_q3 * ( - m_reso**2 )
830  res = res + e1_q4*e2_q3*e3_e4 * ( m_reso**2 )
831 
832  print *,"this code should no longer be used"; stop 1
833 
834 !print *, "res GG old",res
835 !pause
836 
837 
838  endif
839 
840  end subroutine gggzzampl
841 
842 
843 !----- a subroutine for q qbar -> G -> Z -> lept + Z --> 2 lepts
844 !----- all outgoing convention and the following momentum assignment
845 !----- 0 -> bq(p1) + q(p2) + e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
846  subroutine evalamp_qqb_g_vv(p,MY_IDUP,res)
847  implicit none
848  real(dp), intent(out) :: res
849  real(dp), intent(in) :: p(4,6)
850  integer, intent(in) :: MY_IDUP(6:9)
851  real(dp) :: pin(4,4)
852  complex(dp) :: A_VV(1:18), A0_VV(1:2)
853  integer :: i1,i3,i4,VVMode,VVmode_swap
854  real(dp) :: prefactor
855  real(dp) :: intcolfac
856  integer :: ordering(1:4),ordering_swap(1:4)
857  logical :: doInterference
858 
859  if(isaquark(my_idup(6)) .and. isaquark(my_idup(8))) then
860  intcolfac=1.0_dp/3.0_dp
861  else
862  intcolfac=1.0_dp
863  endif
864 
865  call getdecay_vvmode_ordering(my_idup(6:9),vvmode,ordering,vvmode_swap,ordering_swap)
866 
867 !---- full prefactor; 3 is the color factor
868  if( vvmode.eq.zzmode ) then! Z decay
869  prefactor = 3d0*overallcouplvffsq**2
870  elseif( vvmode.eq.wwmode ) then ! W decay
871  prefactor = 3d0*overallcouplvffsq**2
872  elseif( vvmode.eq.zgmode ) then ! Z+photon "decay"
873  prefactor = 3d0*overallcouplvffsq ! Only single powers
874  elseif( vvmode.eq.ggmode ) then ! photon "decay"
875  prefactor = 3d0
876  else
877  prefactor = 0d0
878  endif
879 
880 
881  res = zero
882  a_vv(:) = 0d0
883  dointerference = includeinterference .and. ( &
884  ((vvmode.eq.zzmode) .and. (vvmode_swap.eq.zzmode)) &
885  )
886  if ( includevprime .and. .not.(vvmode.eq.zzmode .or. vvmode.eq.zgmode .or. vvmode.eq.wwmode) ) then
887  call error("Contact terms only for WW, ZZ or Zg!")
888  endif
889  do i1=1,2; do i3=1,2; do i4=1,2! sum over helicities
890  call calchelamp_qq(ordering,vvmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(1))
891  if( vvmode.eq.zzmode ) then
892  if( includegammastar ) then
893  call calchelamp_qq(ordering,zgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(3))
894  call calchelamp_qq(ordering,gszmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(5))
895  call calchelamp_qq(ordering,gsgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(7))
896  endif
897  if( includevprime ) then
898  call calchelamp_qq(ordering,zzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(9))
899  call calchelamp_qq(ordering,zpzmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(11))
900  call calchelamp_qq(ordering,zpzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(13))
901  endif
902  if( includegammastar .and. includevprime ) then
903  call calchelamp_qq(ordering,gszpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(15))
904  call calchelamp_qq(ordering,zpgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(17))
905  endif
906  elseif( vvmode.eq.zgmode ) then
907  if(includegammastar) then
908  call calchelamp_qq(ordering,gsgmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(3))
909  endif
910  if( includevprime ) then
911  call calchelamp_qq(ordering,zpgmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(5))
912  endif
913  elseif( vvmode.eq.wwmode .and. includevprime ) then
914  call calchelamp_qq(ordering,wwpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(9))
915  call calchelamp_qq(ordering,wpwmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(11))
916  call calchelamp_qq(ordering,wpwpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(13))
917  endif
918 
919  if( dointerference ) then
920  call calchelamp_qq(ordering_swap,vvmode_swap,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(2))
921  if( includegammastar ) then
922  call calchelamp_qq(ordering_swap,zgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(4))
923  call calchelamp_qq(ordering_swap,gszmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(6))
924  call calchelamp_qq(ordering_swap,gsgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(8))
925  endif
926  if( includevprime ) then
927  call calchelamp_qq(ordering_swap,zzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(10))
928  call calchelamp_qq(ordering_swap,zpzmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(12))
929  call calchelamp_qq(ordering_swap,zpzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(14))
930  endif
931  if( includegammastar .and. includevprime ) then
932  call calchelamp_qq(ordering_swap,gszpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(16))
933  call calchelamp_qq(ordering_swap,zpgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(18))
934  endif
935  endif
936 
937  a0_vv(1) = a_vv(1)+a_vv(3)+a_vv(5)+a_vv(7)+a_vv(9)+a_vv(11)+a_vv(13)+a_vv(15)+a_vv(17) ! 3456 pieces
938  a0_vv(2) = a_vv(2)+a_vv(4)+a_vv(6)+a_vv(8)+a_vv(10)+a_vv(12)+a_vv(14)+a_vv(16)+a_vv(18) ! 5436 pieces
939  res = res + dreal(a0_vv(1)*dconjg(a0_vv(1)))
940  res = res + dreal(a0_vv(2)*dconjg(a0_vv(2)))
941  if( dointerference .and. (i3.eq.i4) ) then! interfere the 3456 with 5436 pieces
942  res = res - 2d0*intcolfac*dreal( a0_vv(1)*dconjg(a0_vv(2)) ) ! minus from Fermi statistics
943  endif
944  enddo; enddo; enddo
945 
946  res = res*prefactor
947  if( (vvmode.eq.zzmode) .and. dointerference ) res = res * symmfac
948 
949  end subroutine
950 
951  subroutine calchelamp_qq(ordering,VVMode,p,MY_IDUP,i1,i3,i4,A)
952  implicit none
953  integer :: ordering(1:4),i1,i3,i4,l1,l2,l3,l4,MY_IDUP(6:9),VVMode
954  real(dp) :: p(1:4,1:6)
955  complex(dp) :: propG, propZ1, propZ2
956  real(dp) :: s, pin(4,4)
957  complex(dp) :: A(1:1), sp(4,4)
958 
959  l1=ordering(1)
960  l2=ordering(2)
961  l3=ordering(3)
962  l4=ordering(4)
963 
964  !print *,"p(6)=",(p(:,l1))
965  !print *,"p(7)=",(p(:,l2))
966  !print *,"p(8)=",(p(:,l3))
967  !print *,"p(9)=",(p(:,l4))
968  !pause
969 
970  s = 2d0*scr(p(:,1),p(:,2))
971  propg = s/dcmplx(s - m_reso**2,m_reso*ga_reso)
972  pin(1,:) = p(:,1)
973  pin(2,:) = p(:,2)
974  sp(1,:) = pol_dk2mom(dcmplx(p(:,2)),dcmplx(p(:,1)),-3+2*i1) !qbq
975  sp(2,:) = sp(1,:) !-- the same, isn't really needed but for uniform bookeeping
976  call getdecay_couplings_spinors_props( &
977  vvmode, &
978  (/my_idup(l1+3),my_idup(l2+3),my_idup(l3+3),my_idup(l4+3)/), &
979  (/p(:,l1),p(:,l2),p(:,l3),p(:,l4)/), &
980  -3+2*i3,-3+2*i4, &
981  sp(3:4,:),pin(3:4,:) &
982  )
983  call qqgzzampl(vvmode,pin,sp,a(1))
984 !---- chiral couplings of quarks to gravitons
985  if (i1.eq.1) then
986  a(1) = graviton_qq_left*a(1)
987  elseif(i1.eq.2) then
988  a(1) = graviton_qq_right*a(1)
989  endif
990  a(1) = a(1)*propg
991  end subroutine
992 
993  subroutine qqgzzampl(VVMode,p,sp,res)
994  implicit none
995  integer, intent(in) :: VVMode
996  real(dp), intent(in) :: p(4,4)
997  complex(dp), intent(in) :: sp(4,4)
998  complex(dp), intent(out) :: res
999  complex(dp) :: e1_e2, e1_e3, e1_e4
1000  complex(dp) :: e2_e3, e2_e4
1001  complex(dp) :: e3_e4
1002  complex(dp) :: q_q,q3_q3,q4_q4
1003  complex(dp) :: q1_q2,q1_q3,q1_q4
1004  complex(dp) :: q2_q3,q2_q4
1005  complex(dp) :: q3_q4
1006  complex(dp) :: q1_e3,q1_e4,q2_e3,q2_e4
1007  complex(dp) :: e1_q3,e1_q4,e2_q3,e2_q4
1008  complex(dp) :: e3_q4,e4_q3
1009  complex(dp) :: q1(4),q2(4),q3(4),q4(4),q(4)
1010  complex(dp) :: e1(4),e2(4),e3(4),e4(4),abr1
1011  complex(dp) :: yyy1,yyy2,yyy3,yyy41,yyy42,yyy5,yyy6,yyy7,yyy4
1012  complex(dp) :: b_dyn(1:10)
1013  real(dp) :: q34,MG,MZ3,MZ4
1014  real(dp) :: rr
1015 
1016  b_dyn(:)=czero
1017 
1018  q1 = dcmplx(p(1,:),0d0)
1019  q2 = dcmplx(p(2,:),0d0)
1020  q3 = dcmplx(p(3,:),0d0)
1021  q4 = dcmplx(p(4,:),0d0)
1022 
1023 
1024  e1 = sp(1,:)
1025  e2 = sp(2,:)
1026  e3 = sp(3,:)
1027  e4 = sp(4,:)
1028 
1029  q = -q1-q2
1030 
1031  q_q =sc(q,q)
1032  q3_q3 = sc(q3,q3)
1033  q4_q4 = sc(q4,q4)
1034 
1035 
1036  q1_q2 = sc(q1,q2)
1037  q1_q3 = sc(q1,q3)
1038  q1_q4 = sc(q1,q4)
1039  q2_q3 = sc(q2,q3)
1040  q2_q4 = sc(q2,q4)
1041  q3_q4 = sc(q3,q4)
1042 
1043  e1_e2 = sc(e1,e2)
1044  e1_e3 = sc(e1,e3)
1045  e1_e4 = sc(e1,e4)
1046 
1047  e2_e3 = sc(e2,e3)
1048  e2_e4 = sc(e2,e4)
1049 
1050  e3_e4 = sc(e3,e4)
1051 
1052 
1053  q1_e3 = sc(q1,e3)
1054  q1_e4 = sc(q1,e4)
1055  q2_e3 = sc(q2,e3)
1056  q2_e4 = sc(q2,e4)
1057  e1_q3 = sc(e1,q3)
1058  e1_q4 = sc(e1,q4)
1059  e2_q3 = sc(e2,q3)
1060  e2_q4 = sc(e2,q4)
1061  e3_q4 = sc(e3,q4)
1062  e4_q3 = sc(e4,q3)
1063 
1064  mz3=dsqrt(cdabs(q3_q3))
1065  mz4=dsqrt(cdabs(q4_q4))
1066  if( use_dynamic_mg ) then
1067  mg = dsqrt(cdabs(q_q))
1068  else
1069  mg = m_reso
1070  endif
1071 
1072  q34 = (mg**2-mz3**2-mz4**2)/2d0
1073 
1074  if (generate_bis) then
1075  rr = q34/lambda**2! kappa for FS
1076 
1077  if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode) ) then! decay ZZ's or WW's
1078  b_dyn(1)=b1
1079  b_dyn(2)=b2
1080  b_dyn(3)=b3
1081  b_dyn(4)=b4
1082  b_dyn(5)=b5
1083  b_dyn(6)=b6
1084  b_dyn(7)=b7
1085  b_dyn(8)=b8
1086  b_dyn(9)=b9
1087  b_dyn(10)=b10
1088  elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.gszmode) .OR. (vvmode.eq.zgsmode) ) then
1089  b_dyn(1)=bzgs1
1090  b_dyn(2)=bzgs2
1091  b_dyn(3)=bzgs3
1092  b_dyn(4)=bzgs4
1093  b_dyn(8)=bzgs8
1094  elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode) .or. (vvmode.eq.gsgmode) ) then
1095  b_dyn(1)=bgsgs1
1096  b_dyn(2)=bgsgs2
1097  b_dyn(3)=bgsgs3
1098  b_dyn(4)=bgsgs4
1099  b_dyn(8)=bgsgs8
1100  elseif( (vvmode.eq.zzpmode) .or. (vvmode.eq.wwpmode) .or. (vvmode.eq.zpzmode) .or. (vvmode.eq.wpwmode) ) then
1101  b_dyn(1)=bzzp1
1102  b_dyn(2)=bzzp2
1103  b_dyn(3)=bzzp3
1104  b_dyn(4)=bzzp4
1105  b_dyn(5)=bzzp5
1106  b_dyn(6)=bzzp6
1107  b_dyn(7)=bzzp7
1108  b_dyn(8)=bzzp8
1109  b_dyn(9)=bzzp9
1110  b_dyn(10)=bzzp10
1111  elseif( (vvmode.eq.zpzpmode) .or. (vvmode.eq.wpwpmode) ) then
1112  b_dyn(1)=bzpzp1
1113  b_dyn(2)=bzpzp2
1114  b_dyn(3)=bzpzp3
1115  b_dyn(4)=bzpzp4
1116  b_dyn(5)=bzpzp5
1117  b_dyn(6)=bzpzp6
1118  b_dyn(7)=bzpzp7
1119  b_dyn(8)=bzpzp8
1120  b_dyn(9)=bzpzp9
1121  b_dyn(10)=bzpzp10
1122  elseif( (vvmode.eq.zpgmode) .OR. (vvmode.eq.gszpmode) .OR. (vvmode.eq.zpgsmode) ) then
1123  b_dyn(1)=bzpgs1
1124  b_dyn(2)=bzpgs2
1125  b_dyn(3)=bzpgs3
1126  b_dyn(4)=bzpgs4
1127  b_dyn(8)=bzpgs8
1128  else
1129  print *,"VVMode",vvmode,"not implemented"
1130  endif
1131 
1132  yyy1 = q34*( b_dyn(1) + b_dyn(2)*rr*(one+mz3**2/q34)*(one+mz4**2/q34) ) + b_dyn(5)*m_v**2
1133  yyy2 = -b_dyn(1)/two + b_dyn(3)*rr*(1d0-(mz3**2+mz4**2)/(2d0*q34)) + two*b_dyn(4)*rr + b_dyn(7)*rr*m_v**2/q34
1134  yyy3 = (-b_dyn(2)/two - b_dyn(3)- two*b_dyn(4))*rr/q34
1135  yyy41 = -b_dyn(1) - b_dyn(2)*(q34+mz3**2)/lambda**2 - b_dyn(3)*mz4**2/lambda**2 - 2d0*b_dyn(6)*m_v**2/lambda**2
1136  yyy42 = -b_dyn(1) - b_dyn(2)*(q34+mz4**2)/lambda**2 - b_dyn(3)*mz3**2/lambda**2 - 2d0*b_dyn(6)*m_v**2/lambda**2
1137  yyy5 = two*b_dyn(8)*rr*mg**2/q34
1138  yyy6 = b_dyn(9) * m_v**2/lambda**2
1139  yyy7 = b_dyn(10) * mg**2 * m_v**2/lambda**4
1140  else
1141  yyy1 = q34*c1/2d0
1142  yyy2 = c2
1143  yyy3 = c3/mg**2
1144  yyy41 = c41
1145  yyy42 = c42
1146  yyy5 = c5
1147  yyy6 = czero
1148  yyy7 = czero
1149  if(vvmode.eq.zzmode .or. vvmode.eq.wwmode) then
1150  yyy6 = c6
1151  yyy7 = c7
1152  endif
1153  endif
1154 
1155  res = czero
1156 
1157 ! old code without c41 c42 coupligs
1158 ! res = &
1159 ! q1_e3*e1_e4*yyy1 - q1_e3*e1_q3*e4_q3*yyy4 + q1_e4*e1_e3*yyy1 + &
1160 ! q1_e4*e1_q3*e3_q4*yyy4 - q1_q3*e1_e3*e4_q3*yyy4 + q1_q3*e1_e4* &
1161 ! e3_q4*yyy4 + 4*q1_q3*e1_q3*e3_e4*yyy2 + 4*q1_q3*e1_q3*e3_q4* &
1162 ! e4_q3*yyy3 + 1./2.*e1_e3*e4_q3*yyy1 - 1./4.*e1_e3*e4_q3*MG**2* &
1163 ! yyy4 + 1./2.*e1_e4*e3_q4*yyy1 + 1./4.*e1_e4*e3_q4*MG**2*yyy4 + &
1164 ! e1_q3*e3_e4*MG**2*yyy2 + e1_q3*e3_q4*e4_q3*MG**2*yyy3
1165 !
1166 ! res = res + &
1167 ! 1./2.*et1(e3,e4,q1,q)*e1_q3*yyy6 - 1./2.*et1(e3,e4,q2,q)*e1_q3* &
1168 ! yyy6 + 1./4.*et1(e3,e4,e1,q)*MG**2*yyy6 + et1(e3,e4,e1,q)*q1_q3* &
1169 ! yyy6 + 4*et1(e3,e4,q3,q4)*q1_q3*e1_q3*MG**(-2)*yyy5 + et1(e3,e4, &
1170 ! q3,q4)*e1_q3*yyy5
1171 !
1172 ! res = res + &
1173 ! 1./2.*et1(q1,e3,q,q3)*e1_q3*e4_q3*MG**(-2)*yyy7 - 1./2.*et1(q1, &
1174 ! e3,q,q4)*e1_q3*e4_q3*MG**(-2)*yyy7 + 1./2.*et1(q1,e4,q,q3)*e1_q3 &
1175 ! *e3_q4*MG**(-2)*yyy7 - 1./2.*et1(q1,e4,q,q4)*e1_q3*e3_q4* &
1176 ! MG**(-2)*yyy7 - 1./2.*et1(q2,e3,q,q3)*e1_q3*e4_q3*MG**(-2)*yyy7 &
1177 ! + 1./2.*et1(q2,e3,q,q4)*e1_q3*e4_q3*MG**(-2)*yyy7 - 1./2.*et1( &
1178 ! q2,e4,q,q3)*e1_q3*e3_q4*MG**(-2)*yyy7 + 1./2.*et1(q2,e4,q,q4)* &
1179 ! e1_q3*e3_q4*MG**(-2)*yyy7 + et1(e1,e3,q,q3)*q1_q3*e4_q3*MG**(-2) &
1180 ! *yyy7 + 1./4.*et1(e1,e3,q,q3)*e4_q3*yyy7 - et1(e1,e3,q,q4)*q1_q3 &
1181 ! *e4_q3*MG**(-2)*yyy7 - 1./4.*et1(e1,e3,q,q4)*e4_q3*yyy7 + et1(e1 &
1182 ! ,e4,q,q3)*q1_q3*e3_q4*MG**(-2)*yyy7 + 1./4.*et1(e1,e4,q,q3)* &
1183 ! e3_q4*yyy7 - et1(e1,e4,q,q4)*q1_q3*e3_q4*MG**(-2)*yyy7 - 1./4.* &
1184 ! et1(e1,e4,q,q4)*e3_q4*yyy7
1185 ! print *, "res old QQB",res
1186 
1187 
1188 
1189 ! ! new code with c41 c42 coupligs
1190 ! res =&
1191 ! & q1_e3*e1_e4*yyy1 - q1_e3*e1_q3*e4_q3*yyy41 + q1_e4*e1_e3*yyy1 + &
1192 ! & q1_e4*e1_q3*e3_q4*yyy42 - q1_q3*e1_e3*e4_q3*yyy41 + q1_q3*e1_e4*&
1193 ! & e3_q4*yyy42 + 4.*q1_q3*e1_q3*e3_e4*yyy2 + 4.*q1_q3*e1_q3*e3_q4*&
1194 ! & e4_q3*yyy3 + 1./2.*e1_e3*e4_q3*yyy1 + 1./4.*e1_e3*e4_q3*MZ4**2*&
1195 ! & yyy41 - 1./4.*e1_e3*e4_q3*MZ3**2*yyy41 - 1./4.*e1_e3*e4_q3*MG**2&
1196 ! & *yyy41 + 1./2.*e1_e4*e3_q4*yyy1 - 1./4.*e1_e4*e3_q4*MZ4**2*yyy42&
1197 ! & + 1./4.*e1_e4*e3_q4*MZ3**2*yyy42 + 1./4.*e1_e4*e3_q4*MG**2*&
1198 ! & yyy42 - e1_q3*e3_e4*MZ4**2*yyy2 + e1_q3*e3_e4*MZ3**2*yyy2 + &
1199 ! & e1_q3*e3_e4*MG**2*yyy2 + 1./2.*e1_q3*e3_q4*e4_q3*yyy42 - 1./2.*&
1200 ! & e1_q3*e3_q4*e4_q3*yyy41 - e1_q3*e3_q4*e4_q3*MZ4**2*yyy3 + e1_q3*&
1201 ! & e3_q4*e4_q3*MZ3**2*yyy3 + e1_q3*e3_q4*e4_q3*MG**2*yyy3 + 1./2.*&
1202 ! & et1(q1,e3,q,q3)*e1_q3*e4_q3*MG**(-2)*yyy7 - 1./2.*et1(q1,e3,q,q4&
1203 ! & )*e1_q3*e4_q3*MG**(-2)*yyy7 + 1./2.*et1(q1,e4,q,q3)*e1_q3*e3_q4*&
1204 ! & MG**(-2)*yyy7 - 1./2.*et1(q1,e4,q,q4)*e1_q3*e3_q4*MG**(-2)*yyy7&
1205 ! & - 1./2.*et1(q2,e3,q,q3)*e1_q3*e4_q3*MG**(-2)*yyy7
1206 ! res = res + 1./2.*et1(q2,e3,q,q4)*e1_q3*e4_q3*MG**(-2)*yyy7 - 1./&
1207 ! & 2.*et1(q2,e4,q,q3)*e1_q3*e3_q4*MG**(-2)*yyy7 + 1./2.*et1(q2,e4,q&
1208 ! & ,q4)*e1_q3*e3_q4*MG**(-2)*yyy7 + et1(e1,e3,q,q3)*q1_q3*e4_q3*&
1209 ! & MG**(-2)*yyy7 - 1./4.*et1(e1,e3,q,q3)*e4_q3*MG**(-2)*MZ4**2*yyy7&
1210 ! & + 1./4.*et1(e1,e3,q,q3)*e4_q3*MG**(-2)*MZ3**2*yyy7 + 1./4.*et1(&
1211 ! & e1,e3,q,q3)*e4_q3*yyy7 - et1(e1,e3,q,q4)*q1_q3*e4_q3*MG**(-2)*&
1212 ! & yyy7 + 1./4.*et1(e1,e3,q,q4)*e4_q3*MG**(-2)*MZ4**2*yyy7 - 1./4.*&
1213 ! & et1(e1,e3,q,q4)*e4_q3*MG**(-2)*MZ3**2*yyy7 - 1./4.*et1(e1,e3,q,&
1214 ! & q4)*e4_q3*yyy7 + et1(e1,e4,q,q3)*q1_q3*e3_q4*MG**(-2)*yyy7 - 1./&
1215 ! & 4.*et1(e1,e4,q,q3)*e3_q4*MG**(-2)*MZ4**2*yyy7 + 1./4.*et1(e1,e4,&
1216 ! & q,q3)*e3_q4*MG**(-2)*MZ3**2*yyy7 + 1./4.*et1(e1,e4,q,q3)*e3_q4*&
1217 ! & yyy7 - et1(e1,e4,q,q4)*q1_q3*e3_q4*MG**(-2)*yyy7 + 1./4.*et1(e1,&
1218 ! & e4,q,q4)*e3_q4*MG**(-2)*MZ4**2*yyy7 - 1./4.*et1(e1,e4,q,q4)*&
1219 ! & e3_q4*MG**(-2)*MZ3**2*yyy7 - 1./4.*et1(e1,e4,q,q4)*e3_q4*yyy7 + &
1220 ! & 1./2.*et1(e3,e4,q1,q)*e1_q3*yyy6 - 1./2.*et1(e3,e4,q2,q)*e1_q3*&
1221 ! & yyy6
1222 ! res = res - 1./4.*et1(e3,e4,e1,q)*MZ4**2*yyy6 + 1./4.*et1(e3,e4,&
1223 ! & e1,q)*MZ3**2*yyy6 + 1./4.*et1(e3,e4,e1,q)*MG**2*yyy6 + et1(e3,e4&
1224 ! & ,e1,q)*q1_q3*yyy6 + 4.*et1(e3,e4,q3,q4)*q1_q3*e1_q3*MG**(-2)*&
1225 ! & yyy5 - et1(e3,e4,q3,q4)*e1_q3*MG**(-2)*MZ4**2*yyy5 + et1(e3,e4,&
1226 ! & q3,q4)*e1_q3*MG**(-2)*MZ3**2*yyy5 + et1(e3,e4,q3,q4)*e1_q3*yyy5
1227 
1228 ! print *, "res new QQB",res
1229 
1230 
1231 
1232 ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! SAME AS ABOVE BUT SHORTER ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
1233 
1234  abr1 = (mg**2 + mz3**2 - mz4**2 + 4*q1_q3)
1235 
1236  res = (e1_e3*(4*yyy1*q1_e4 + e4_q3*(2*yyy1 - (mg**2 + mz3**2 - mz4**2)*yyy41 - 4*yyy41*q1_q3)) + &
1237  e1_e4*(4*yyy1*q1_e3 + e3_q4*(2*yyy1 + (mg**2 + mz3**2 - mz4**2)*yyy42 + 4*yyy42*q1_q3)) + &
1238  2*e1_q3*(-2*yyy41*e4_q3*q1_e3 + 2*yyy2*e3_e4*abr1 + &
1239  e3_q4*(2*yyy42*q1_e4 + e4_q3*(2*(mg**2 + mz3**2 - mz4**2)*yyy3 - yyy41 + yyy42 + 8*yyy3*q1_q3))))/4d0 + &
1240  (yyy6*abr1*et1(e3,e4,e1,q))/4d0 + (yyy6*e1_q3*et1(e3,e4,q1,q))/2d0 &
1241  - (yyy6*e1_q3*et1(e3,e4,q2,q))/2d0 + &
1242  (yyy5*e1_q3*abr1*et1(e3,e4,q3,q4))/mg**2 + &
1243  yyy7*((e4_q3*abr1*et1(e1,e3,q,q3))/(4d0*mg**2) - &
1244  (e4_q3*abr1*et1(e1,e3,q,q4))/(4d0*mg**2) + &
1245  (e3_q4*abr1*et1(e1,e4,q,q3))/(4d0*mg**2) - &
1246  (e3_q4*abr1*et1(e1,e4,q,q4))/(4d0*mg**2) &
1247  + (e1_q3*e4_q3*et1(q1,e3,q,q3))/(2d0*mg**2) - &
1248  (e1_q3*e4_q3*et1(q1,e3,q,q4))/(2d0*mg**2) + (e1_q3*e3_q4*et1(q1,e4,q,q3))/(2d0*mg**2) - &
1249  (e1_q3*e3_q4*et1(q1,e4,q,q4))/(2d0*mg**2) - (e1_q3*e4_q3*et1(q2,e3,q,q3))/(2d0*mg**2) + &
1250  (e1_q3*e4_q3*et1(q2,e3,q,q4))/(2d0*mg**2) - (e1_q3*e3_q4*et1(q2,e4,q,q3))/(2d0*mg**2) + &
1251  (e1_q3*e3_q4*et1(q2,e4,q,q4))/(2d0*mg**2))
1252 
1253 ! print *, "res newer QQB ",res
1254 ! pause
1255 
1256  end subroutine qqgzzampl
1257 
1258 
1259 !----- a subroutine for G -> ZZ/WW/AA
1260 !----- all outgoing convention and the following momentum assignment
1261 !----- 0 -> G(p1) + e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
1262  subroutine evalamp_g_vv(p,MY_IDUP,res)
1263  implicit none
1264  real(dp), intent(out) :: res
1265  real(dp), intent(in) :: p(4,6)
1266  integer, intent(in) :: MY_IDUP(6:9)
1267  complex(dp) :: A_VV(1:18),A0_VV(1:2)
1268  integer :: i1,i3,i4,VVMode,VVmode_swap
1269  real(dp) :: prefactor
1270  real(dp) :: intcolfac
1271  integer :: ordering(1:4),ordering_swap(1:4)
1272  logical :: doInterference
1273 
1274  if(isaquark(my_idup(6)) .and. isaquark(my_idup(8))) then
1275  intcolfac=1.0_dp/3.0_dp
1276  else
1277  intcolfac=1.0_dp
1278  endif
1279 
1280  call getdecay_vvmode_ordering(my_idup(6:9),vvmode,ordering,vvmode_swap,ordering_swap)
1281 
1282 !---- full prefactor
1283  if( vvmode.eq.zzmode ) then! Z decay
1284  prefactor = overallcouplvffsq**2
1285  elseif( vvmode.eq.wwmode ) then ! W decay
1286  prefactor = overallcouplvffsq**2
1287  elseif( vvmode.eq.zgmode ) then ! Z+photon "decay"
1288  prefactor = overallcouplvffsq ! Only single powers
1289  elseif( vvmode.eq.ggmode ) then ! photon "decay"
1290  prefactor = 1d0
1291  else
1292  prefactor = 0d0
1293  endif
1294 
1295  res = zero
1296  a_vv(:) = 0d0
1297  dointerference = includeinterference .and. ( &
1298  ((vvmode.eq.zzmode) .and. (vvmode_swap.eq.zzmode)) &
1299  )
1300  if ( includevprime .and. .not.(vvmode.eq.zzmode .or. vvmode.eq.zgmode .or. vvmode.eq.wwmode) ) then
1301  call error("Contact terms only for WW, ZZ or Zg!")
1302  endif
1303  do i1 =-2,2; do i3=1,2; do i4=1,2! sum over helicities
1304  call calchelamp2(ordering,vvmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(1))
1305  if( vvmode.eq.zzmode ) then
1306  if( includegammastar ) then
1307  call calchelamp2(ordering,zgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(3))
1308  call calchelamp2(ordering,gszmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(5))
1309  call calchelamp2(ordering,gsgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(7))
1310  endif
1311  if( includevprime ) then
1312  call calchelamp2(ordering,zzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(9))
1313  call calchelamp2(ordering,zpzmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(11))
1314  call calchelamp2(ordering,zpzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(13))
1315  endif
1316  if( includegammastar .and. includevprime ) then
1317  call calchelamp2(ordering,gszpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(15))
1318  call calchelamp2(ordering,zpgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(17))
1319  endif
1320  elseif( vvmode.eq.zgmode ) then
1321  if(includegammastar) then
1322  call calchelamp2(ordering,gsgmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(3))
1323  endif
1324  if( includevprime ) then
1325  call calchelamp2(ordering,zpgmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(5))
1326  endif
1327  elseif( vvmode.eq.wwmode .and. includevprime ) then
1328  call calchelamp2(ordering,wwpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(9))
1329  call calchelamp2(ordering,wpwmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(11))
1330  call calchelamp2(ordering,wpwpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(13))
1331  endif
1332 
1333  if( dointerference ) then
1334  call calchelamp2(ordering_swap,vvmode_swap,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(2))
1335  if( includegammastar ) then
1336  call calchelamp2(ordering_swap,zgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(4))
1337  call calchelamp2(ordering_swap,gszmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(6))
1338  call calchelamp2(ordering_swap,gsgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(8))
1339  endif
1340  if( includevprime ) then
1341  call calchelamp2(ordering_swap,zzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(10))
1342  call calchelamp2(ordering_swap,zpzmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(12))
1343  call calchelamp2(ordering_swap,zpzpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(14))
1344  endif
1345  if( includegammastar .and. includevprime ) then
1346  call calchelamp2(ordering_swap,gszpmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(16))
1347  call calchelamp2(ordering_swap,zpgsmode,p(1:4,1:6),my_idup,i1,i3,i4,a_vv(18))
1348  endif
1349  endif
1350 
1351  a0_vv(1) = a_vv(1)+a_vv(3)+a_vv(5)+a_vv(7)+a_vv(9)+a_vv(11)+a_vv(13)+a_vv(15)+a_vv(17) ! 3456 pieces
1352  a0_vv(2) = a_vv(2)+a_vv(4)+a_vv(6)+a_vv(8)+a_vv(10)+a_vv(12)+a_vv(14)+a_vv(16)+a_vv(18) ! 5436 pieces
1353  res = res + dreal(a0_vv(1)*dconjg(a0_vv(1)))
1354  res = res + dreal(a0_vv(2)*dconjg(a0_vv(2)))
1355  if( dointerference .and. (i3.eq.i4) ) then! interfere the 3456 with 5436 pieces
1356  res = res - 2d0*intcolfac*dreal( a0_vv(1)*dconjg(a0_vv(2)) ) ! minus from Fermi statistics
1357  endif
1358  enddo; enddo; enddo
1359 
1360  res = res*prefactor
1361  if( (vvmode.eq.zzmode) .and. dointerference ) res = res * symmfac
1362 
1363  end subroutine
1364 
1365  subroutine calchelamp2(ordering,VVMode,p,MY_IDUP,i1,i3,i4,A)
1366  implicit none
1367  integer :: ordering(1:4),i1,i3,i4,l1,l2,l3,l4,mu,nu,MY_IDUP(6:9),VVMode
1368  real(dp) :: p(1:4,1:6)
1369  real(dp) :: pin(4,4)
1370  complex(dp) :: A(1:1), sp(0:4,1:4)
1371 
1372  l1=ordering(1)
1373  l2=ordering(2)
1374  l3=ordering(3)
1375  l4=ordering(4)
1376 
1377  pin(1,:) = p(:,1)
1378  pin(2,:) = 0d0 ! dummy
1379 
1380  sp(0,1:4) = pol_mass2(dcmplx(p(1:4,1)), 0,'in')
1381  sp(1,1:4) = pol_mass2(dcmplx(p(1:4,1)),-1,'in')
1382  sp(2,1:4) = pol_mass2(dcmplx(p(1:4,1)),+1,'in')
1383  call getdecay_couplings_spinors_props( &
1384  vvmode, &
1385  (/my_idup(l1+3),my_idup(l2+3),my_idup(l3+3),my_idup(l4+3)/), &
1386  (/p(:,l1),p(:,l2),p(:,l3),p(:,l4)/), &
1387  -3+2*i3,-3+2*i4, &
1388  sp(3:4,:),pin(3:4,:) &
1389  )
1390  call gzzampl(vvmode,pin,sp,i1,a(1))
1391  end subroutine
1392 
1393  subroutine gzzampl(VVMode,p,sp,i1,res)
1394  implicit none
1395  integer, intent(in) :: VVMode
1396  real(dp), intent(in) :: p(4,4)
1397  complex(dp), intent(in) :: sp(0:4,4)
1398  integer,intent(in) :: i1
1399  complex(dp), intent(out) :: res
1400  complex(dp) :: e1_e2, e1_e3, e1_e4
1401  complex(dp) :: e2_e3, e2_e4
1402  complex(dp) :: e3_e4
1403  complex(dp) :: q_q,q3_q3,q4_q4
1404  complex(dp) :: q1_q2,q1_q3,q1_q4
1405  complex(dp) :: q2_q3,q2_q4
1406  complex(dp) :: q3_q4
1407  complex(dp) :: q1_e3,q1_e4,q2_e3,q2_e4,e0_e3,e0_e4
1408  complex(dp) :: e1_q3,e1_q4,e2_q3,e2_q4,e0_q3,e0_q4,q_e3,q_e4
1409  complex(dp) :: e3_q4,e4_q3
1410  complex(dp) :: q1(4),q2(4),q3(4),q4(4),q(4)
1411  complex(dp) :: e1(4),e2(4),e3(4),e4(4),e0(4)
1412  complex(dp) :: yyy1,yyy2,yyy3,yyy41,yyy42,yyy5,yyy6,yyy7
1413  complex(dp) :: b_dyn(1:10)
1414  real(dp) :: q34,MG,MZ3,MZ4
1415  real(dp) :: rr
1416  real(dp),parameter :: sqrt6=dsqrt(6d0)
1417 
1418  b_dyn(:)=czero
1419 
1420  q1 = dcmplx(p(1,:),0d0)
1421  q2 = dcmplx(p(2,:),0d0)
1422  q3 = dcmplx(p(3,:),0d0)
1423  q4 = dcmplx(p(4,:),0d0)
1424 
1425 ! requirement: sp(0,:)= 0 pol
1426 ! requirement: sp(1,:)= - pol
1427 ! requirement: sp(2,:)= + pol
1428  if( i1.eq.+2 ) then
1429  e0 = 1d9! dummy
1430  e1 = sp(2,:)! +
1431  e2 = sp(2,:)! +
1432  elseif( i1.eq.+1 ) then
1433  e0 = 1d9! dummy
1434  e1 = sp(2,:)! +
1435  e2 = sp(0,:)! 0
1436  elseif( i1.eq.0 ) then
1437  e0 = sp(0,:)! 0
1438  e1 = sp(1,:)! -
1439  e2 = sp(2,:)! +
1440  elseif( i1.eq.-1 ) then
1441  e0 = 1d9! dummy
1442  e1 = sp(1,:)! -
1443  e2 = sp(0,:)! 0
1444  elseif( i1.eq.-2 ) then
1445  e0 = 1d9! dummy
1446  e1 = sp(1,:)! -
1447  e2 = sp(1,:)! -
1448  endif
1449 
1450  e3 = sp(3,:)
1451  e4 = sp(4,:)
1452 
1453  q = -q1-q2
1454 
1455  q_q =sc(q,q)
1456  q3_q3 = sc(q3,q3)
1457  q4_q4 = sc(q4,q4)
1458 
1459 
1460  q1_q2 = sc(q1,q2)
1461  q1_q3 = sc(q1,q3)
1462  q1_q4 = sc(q1,q4)
1463  q2_q3 = sc(q2,q3)
1464  q2_q4 = sc(q2,q4)
1465  q3_q4 = sc(q3,q4)
1466 
1467  e1_e2 = sc(e1,e2)
1468  e1_e3 = sc(e1,e3)
1469  e1_e4 = sc(e1,e4)
1470 
1471  e2_e3 = sc(e2,e3)
1472  e2_e4 = sc(e2,e4)
1473 
1474  e3_e4 = sc(e3,e4)
1475 
1476 ! new
1477  e0_e3 = sc(e0,e3)
1478  e0_e4 = sc(e0,e4)
1479  e0_q3 = sc(q3,e0)
1480  e0_q4 = sc(q4,e0)
1481  q_e3 = sc(q4,e3)
1482  q_e4 = sc(q3,e4)
1483 
1484  q1_e3 = sc(q1,e3)
1485  q1_e4 = sc(q1,e4)
1486  q2_e3 = sc(q2,e3)
1487  q2_e4 = sc(q2,e4)
1488  e1_q3 = sc(e1,q3)
1489  e1_q4 = sc(e1,q4)
1490  e2_q3 = sc(e2,q3)
1491  e2_q4 = sc(e2,q4)
1492  e3_q4 = sc(e3,q4)
1493  e4_q3 = sc(e4,q3)
1494 
1495  mz3=dsqrt(cdabs(q3_q3))
1496  mz4=dsqrt(cdabs(q4_q4))
1497  if( use_dynamic_mg ) then
1498  mg = dsqrt(cdabs(q_q))
1499  else
1500  mg = m_reso
1501  endif
1502 
1503  q34 = (mg**2-mz3**2-mz4**2)/2d0
1504 
1505  if (generate_bis) then
1506  rr = q34/lambda**2! kappa for FS
1507 
1508  if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode) ) then! decay ZZ's or WW's
1509  b_dyn(1)=b1
1510  b_dyn(2)=b2
1511  b_dyn(3)=b3
1512  b_dyn(4)=b4
1513  b_dyn(5)=b5
1514  b_dyn(6)=b6
1515  b_dyn(7)=b7
1516  b_dyn(8)=b8
1517  b_dyn(9)=b9
1518  b_dyn(10)=b10
1519  elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.gszmode) .OR. (vvmode.eq.zgsmode) ) then
1520  b_dyn(1)=bzgs1
1521  b_dyn(2)=bzgs2
1522  b_dyn(3)=bzgs3
1523  b_dyn(4)=bzgs4
1524  b_dyn(8)=bzgs8
1525  elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode) .or. (vvmode.eq.gsgmode) ) then
1526  b_dyn(1)=bgsgs1
1527  b_dyn(2)=bgsgs2
1528  b_dyn(3)=bgsgs3
1529  b_dyn(4)=bgsgs4
1530  b_dyn(8)=bgsgs8
1531  elseif( (vvmode.eq.zzpmode) .or. (vvmode.eq.wwpmode) .or. (vvmode.eq.zpzmode) .or. (vvmode.eq.wpwmode) ) then
1532  b_dyn(1)=bzzp1
1533  b_dyn(2)=bzzp2
1534  b_dyn(3)=bzzp3
1535  b_dyn(4)=bzzp4
1536  b_dyn(5)=bzzp5
1537  b_dyn(6)=bzzp6
1538  b_dyn(7)=bzzp7
1539  b_dyn(8)=bzzp8
1540  b_dyn(9)=bzzp9
1541  b_dyn(10)=bzzp10
1542  elseif( (vvmode.eq.zpzpmode) .or. (vvmode.eq.wpwpmode) ) then
1543  b_dyn(1)=bzpzp1
1544  b_dyn(2)=bzpzp2
1545  b_dyn(3)=bzpzp3
1546  b_dyn(4)=bzpzp4
1547  b_dyn(5)=bzpzp5
1548  b_dyn(6)=bzpzp6
1549  b_dyn(7)=bzpzp7
1550  b_dyn(8)=bzpzp8
1551  b_dyn(9)=bzpzp9
1552  b_dyn(10)=bzpzp10
1553  elseif( (vvmode.eq.zpgmode) .OR. (vvmode.eq.gszpmode) .OR. (vvmode.eq.zpgsmode) ) then
1554  b_dyn(1)=bzpgs1
1555  b_dyn(2)=bzpgs2
1556  b_dyn(3)=bzpgs3
1557  b_dyn(4)=bzpgs4
1558  b_dyn(8)=bzpgs8
1559  else
1560  print *,"VVMode",vvmode,"not implemented"
1561  endif
1562 
1563  yyy1 = q34*( b_dyn(1) + b_dyn(2)*rr*(one+mz3**2/q34)*(one+mz4**2/q34) ) + b_dyn(5)*m_v**2
1564  yyy2 = -b_dyn(1)/two + b_dyn(3)*rr*(1d0-(mz3**2+mz4**2)/(2d0*q34)) + two*b_dyn(4)*rr + b_dyn(7)*rr*m_v**2/q34
1565  yyy3 = (-b_dyn(2)/two - b_dyn(3)- two*b_dyn(4))*rr/q34
1566  yyy41 = -b_dyn(1) - b_dyn(2)*(q34+mz3**2)/lambda**2 - b_dyn(3)*mz4**2/lambda**2 - 2d0*b_dyn(6)*m_v**2/lambda**2
1567  yyy42 = -b_dyn(1) - b_dyn(2)*(q34+mz4**2)/lambda**2 - b_dyn(3)*mz3**2/lambda**2 - 2d0*b_dyn(6)*m_v**2/lambda**2
1568  yyy5 = two*b_dyn(8)*rr*mg**2/q34
1569  yyy6 = b_dyn(9) * m_v**2/lambda**2
1570  yyy7 = b_dyn(10) * mg**2 * m_v**2/lambda**4
1571 
1572  else
1573  yyy1 = q34*c1/2d0
1574  yyy2 = c2
1575  yyy3 = c3/mg**2
1576  yyy41 = c41
1577  yyy42 = c42
1578  yyy5 = c5
1579  yyy6 = czero
1580  yyy7 = czero
1581  if(vvmode.eq.zzmode .or. vvmode.eq.wwmode) then
1582  yyy6 = c6
1583  yyy7 = c7
1584  endif
1585  endif
1586 
1587 
1588  res = czero
1589 
1590 
1591  if( abs(i1).eq.2 ) then
1592 
1593  res =&
1594  & + yyy7 * ( et1(e1,e3,q,q3)*q_e4*e2_q3*mg**(-2) - et1(e1,e3,q,q3&
1595  & )*q_e4*e2_q4*mg**(-2) - et1(e1,e3,q,q4)*q_e4*e2_q3*mg**(-2)&
1596  & + et1(e1,e3,q,q4)*q_e4*e2_q4*mg**(-2) + et1(e1,e4,q,q3)*q_e3&
1597  & *e2_q3*mg**(-2) - et1(e1,e4,q,q3)*q_e3*e2_q4*mg**(-2) - et1(&
1598  & e1,e4,q,q4)*q_e3*e2_q3*mg**(-2) + et1(e1,e4,q,q4)*q_e3*e2_q4*&
1599  & mg**(-2) )
1600  res = res + yyy6 * ( et1(e3,e4,e1,q)*e2_q3 - et1(e3,e4,e1,q)*&
1601  & e2_q4 )
1602  res = res + yyy5 * ( et1(e3,e4,q3,q4)*e1_q3*e2_q3*mg**(-2) - et1(&
1603  & e3,e4,q3,q4)*e1_q3*e2_q4*mg**(-2) - et1(e3,e4,q3,q4)*e1_q4*&
1604  & e2_q3*mg**(-2) + et1(e3,e4,q3,q4)*e1_q4*e2_q4*mg**(-2) )
1605  res = res + yyy42 * ( e1_e4*e2_q3*e3_q4 + e1_q3*e2_e4*e3_q4 )
1606  res = res + yyy41 * ( e1_e3*e2_q4*e4_q3 + e1_q4*e2_e3*e4_q3 )
1607  res = res + yyy3 * ( e1_q3*e2_q3*e3_q4*e4_q3 - e1_q3*e2_q4*e3_q4*&
1608  & e4_q3 - e1_q4*e2_q3*e3_q4*e4_q3 + e1_q4*e2_q4*e3_q4*e4_q3 )
1609  res = res + yyy2 * ( e1_q3*e2_q3*e3_e4 - e1_q3*e2_q4*e3_e4 - &
1610  & e1_q4*e2_q3*e3_e4 + e1_q4*e2_q4*e3_e4 )
1611  res = res + yyy1 * ( e1_e3*e2_e4 + e1_e4*e2_e3 )
1612 
1613 
1614  elseif( abs(i1).eq.1 ) then
1615 
1616  res =&
1617  & + yyy7*sqrt2**(-1) * ( et1(e1,e3,q,q3)*q_e4*e2_q3*mg**(-2) - &
1618  & et1(e1,e3,q,q3)*q_e4*e2_q4*mg**(-2) - et1(e1,e3,q,q4)*q_e4*&
1619  & e2_q3*mg**(-2) + et1(e1,e3,q,q4)*q_e4*e2_q4*mg**(-2) + et1(e1&
1620  & ,e4,q,q3)*q_e3*e2_q3*mg**(-2) - et1(e1,e4,q,q3)*q_e3*e2_q4*&
1621  & mg**(-2) - et1(e1,e4,q,q4)*q_e3*e2_q3*mg**(-2) + et1(e1,e4,q,&
1622  & q4)*q_e3*e2_q4*mg**(-2) + et1(e2,e3,q,q3)*q_e4*e1_q3*mg**(-2)&
1623  & - et1(e2,e3,q,q3)*q_e4*e1_q4*mg**(-2) - et1(e2,e3,q,q4)*q_e4&
1624  & *e1_q3*mg**(-2) + et1(e2,e3,q,q4)*q_e4*e1_q4*mg**(-2) + et1(&
1625  & e2,e4,q,q3)*q_e3*e1_q3*mg**(-2) - et1(e2,e4,q,q3)*q_e3*e1_q4*&
1626  & mg**(-2) - et1(e2,e4,q,q4)*q_e3*e1_q3*mg**(-2) + et1(e2,e4,q,&
1627  & q4)*q_e3*e1_q4*mg**(-2) )
1628  res = res + yyy6*sqrt2**(-1) * ( et1(e3,e4,e1,q)*e2_q3 - et1(e3,&
1629  & e4,e1,q)*e2_q4 + et1(e3,e4,e2,q)*e1_q3 - et1(e3,e4,e2,q)*&
1630  & e1_q4 )
1631  res = res + yyy5*sqrt2**(-1) * ( 2.0d0*et1(e3,e4,q3,q4)*e1_q3*e2_q3*&
1632  & mg**(-2) - 2.0d0*et1(e3,e4,q3,q4)*e1_q3*e2_q4*mg**(-2) - 2.0d0*et1(&
1633  & e3,e4,q3,q4)*e1_q4*e2_q3*mg**(-2) + 2.0d0*et1(e3,e4,q3,q4)*e1_q4&
1634  & *e2_q4*mg**(-2) )
1635  res = res + yyy42*sqrt2**(-1) * ( 2.0d0*e1_e4*e2_q3*e3_q4 + 2.0d0*e1_q3&
1636  & *e2_e4*e3_q4 )
1637  res = res + yyy41*sqrt2**(-1) * ( 2.0d0*e1_e3*e2_q4*e4_q3 + 2.0d0*e1_q4&
1638  & *e2_e3*e4_q3 )
1639  res = res + yyy3*sqrt2**(-1) * ( 2.0d0*e1_q3*e2_q3*e3_q4*e4_q3 - 2.0d0*&
1640  & e1_q3*e2_q4*e3_q4*e4_q3 - 2.0d0*e1_q4*e2_q3*e3_q4*e4_q3 + 2.0d0*&
1641  & e1_q4*e2_q4*e3_q4*e4_q3 )
1642  res = res + yyy2*sqrt2**(-1) * ( 2.0d0*e1_q3*e2_q3*e3_e4 - 2.0d0*e1_q3*&
1643  & e2_q4*e3_e4 - 2.0d0*e1_q4*e2_q3*e3_e4 + 2.0d0*e1_q4*e2_q4*e3_e4 )
1644  res = res + yyy1*sqrt2**(-1) * ( 2.0d0*e1_e3*e2_e4 + 2.0d0*e1_e4*e2_e3&
1645  & )
1646 
1647 
1648  elseif( abs(i1).eq.0 ) then
1649  res =&
1650  & + yyy7*sqrt6**(-1) * ( - 2.0d0*et1(e0,e3,q,q3)*q_e4*e0_q3*&
1651  & mg**(-2) + 2.0d0*et1(e0,e3,q,q3)*q_e4*e0_q4*mg**(-2) + 2.0d0*et1(e0&
1652  & ,e3,q,q4)*q_e4*e0_q3*mg**(-2) - 2.0d0*et1(e0,e3,q,q4)*q_e4*e0_q4&
1653  & *mg**(-2) - 2.0d0*et1(e0,e4,q,q3)*q_e3*e0_q3*mg**(-2) + 2.0d0*et1(&
1654  & e0,e4,q,q3)*q_e3*e0_q4*mg**(-2) + 2.0d0*et1(e0,e4,q,q4)*q_e3*&
1655  & e0_q3*mg**(-2) - 2.0d0*et1(e0,e4,q,q4)*q_e3*e0_q4*mg**(-2) + &
1656  & et1(e1,e3,q,q3)*q_e4*e2_q3*mg**(-2) - et1(e1,e3,q,q3)*q_e4*&
1657  & e2_q4*mg**(-2) - et1(e1,e3,q,q4)*q_e4*e2_q3*mg**(-2) + et1(e1&
1658  & ,e3,q,q4)*q_e4*e2_q4*mg**(-2) + et1(e1,e4,q,q3)*q_e3*e2_q3*&
1659  & mg**(-2) - et1(e1,e4,q,q3)*q_e3*e2_q4*mg**(-2) - et1(e1,e4,q,&
1660  & q4)*q_e3*e2_q3*mg**(-2) + et1(e1,e4,q,q4)*q_e3*e2_q4*mg**(-2)&
1661  & + et1(e2,e3,q,q3)*q_e4*e1_q3*mg**(-2) - et1(e2,e3,q,q3)*q_e4&
1662  & *e1_q4*mg**(-2) - et1(e2,e3,q,q4)*q_e4*e1_q3*mg**(-2) + et1(&
1663  & e2,e3,q,q4)*q_e4*e1_q4*mg**(-2) + et1(e2,e4,q,q3)*q_e3*e1_q3*&
1664  & mg**(-2) - et1(e2,e4,q,q3)*q_e3*e1_q4*mg**(-2) - et1(e2,e4,q,&
1665  & q4)*q_e3*e1_q3*mg**(-2) )
1666  res = res + yyy7*sqrt6**(-1) * ( et1(e2,e4,q,q4)*q_e3*e1_q4*&
1667  & mg**(-2) )
1668  res = res + yyy6*sqrt6**(-1) * ( - 2.0d0*et1(e3,e4,e0,q)*e0_q3 + 2.0d0&
1669  & *et1(e3,e4,e0,q)*e0_q4 + et1(e3,e4,e1,q)*e2_q3 - et1(e3,e4,e1&
1670  & ,q)*e2_q4 + et1(e3,e4,e2,q)*e1_q3 - et1(e3,e4,e2,q)*e1_q4 )
1671  res = res + yyy5*sqrt6**(-1) * ( 4.0d0*et1(e3,e4,q3,q4)*e0_q3*e0_q4*&
1672  & mg**(-2) - 2.0d0*et1(e3,e4,q3,q4)*e0_q3**2*mg**(-2) - 2.0d0*et1(e3,&
1673  & e4,q3,q4)*e0_q4**2*mg**(-2) + 2.0d0*et1(e3,e4,q3,q4)*e1_q3*e2_q3&
1674  & *mg**(-2) - 2.0d0*et1(e3,e4,q3,q4)*e1_q3*e2_q4*mg**(-2) - 2.0d0*&
1675  & et1(e3,e4,q3,q4)*e1_q4*e2_q3*mg**(-2) + 2.0d0*et1(e3,e4,q3,q4)*&
1676  & e1_q4*e2_q4*mg**(-2) )
1677  res = res + yyy42*sqrt6**(-1) * ( - 4.0d0*e0_e4*e0_q3*e3_q4 + 2.0d0*&
1678  & e1_e4*e2_q3*e3_q4 + 2.0d0*e1_q3*e2_e4*e3_q4 )
1679  res = res + yyy41*sqrt6**(-1) * ( - 4.0d0*e0_e3*e0_q4*e4_q3 + 2.0d0*&
1680  & e1_e3*e2_q4*e4_q3 + 2.0d0*e1_q4*e2_e3*e4_q3 )
1681  res = res + yyy3*sqrt6**(-1) * ( 4.0d0*e0_q3*e0_q4*e3_q4*e4_q3 - 2.0d0*&
1682  & e0_q3**2*e3_q4*e4_q3 - 2.0d0*e0_q4**2*e3_q4*e4_q3 + 2.0d0*e1_q3*&
1683  & e2_q3*e3_q4*e4_q3 - 2.0d0*e1_q3*e2_q4*e3_q4*e4_q3 - 2.0d0*e1_q4*&
1684  & e2_q3*e3_q4*e4_q3 + 2.0d0*e1_q4*e2_q4*e3_q4*e4_q3 )
1685  res = res + yyy2*sqrt6**(-1) * ( 4.0d0*e0_q3*e0_q4*e3_e4 - 2.0d0*&
1686  & e0_q3**2*e3_e4 - 2.0d0*e0_q4**2*e3_e4 + 2.0d0*e1_q3*e2_q3*e3_e4 - 2.0d0&
1687  & *e1_q3*e2_q4*e3_e4 - 2.0d0*e1_q4*e2_q3*e3_e4 + 2.0d0*e1_q4*e2_q4*&
1688  & e3_e4 )
1689  res = res + yyy1*sqrt6**(-1) * ( - 4.0d0*e0_e3*e0_e4 + 2.0d0*e1_e3*&
1690  & e2_e4 + 2.0d0*e1_e4*e2_e3 )
1691 
1692 
1693  endif
1694 
1695  end subroutine gzzampl
1696 
1697 
1698 subroutine getdecay_couplings_spinors_props(VVMode,idordered,pordered,h3,h4, sp,pV)
1699  implicit none
1700  integer, intent(in) :: VVMode,idordered(6:9),h3,h4
1701  real(dp), intent(in) :: pordered(1:4,6:9)
1702  complex(dp), intent(out) :: sp(3:4,1:4)
1703  real(dp), intent(out) :: pV(3:4,1:4)
1704  real(dp) :: s, aL1,aR1,aL2,aR2
1705  complex(dp) :: propV(1:2)
1706 
1707  ! h3/h4 helicities: -1 == left, 1 == right
1708  if( vvmode.eq.zzmode ) then
1709  ! ZZ DECAYS
1710  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1711  al1=al_lep * dsqrt(scale_alpha_z_ll)
1712  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
1713  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1714  al1=al_lep * dsqrt(scale_alpha_z_tt)
1715  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
1716  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1717  al1=al_neu * dsqrt(scale_alpha_z_nn)
1718  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
1719  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1720  al1=al_qup * dsqrt(scale_alpha_z_uu)
1721  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
1722  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1723  al1=al_qdn * dsqrt(scale_alpha_z_dd)
1724  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
1725  else
1726  al1=0d0
1727  ar1=0d0
1728  endif
1729  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
1730  al2=al_lep * dsqrt(scale_alpha_z_ll)
1731  ar2=ar_lep * dsqrt(scale_alpha_z_ll)
1732  elseif( abs(idordered(8)).eq.abs(tam_) ) then
1733  al2=al_lep * dsqrt(scale_alpha_z_tt)
1734  ar2=ar_lep * dsqrt(scale_alpha_z_tt)
1735  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
1736  al2=al_neu * dsqrt(scale_alpha_z_nn)
1737  ar2=ar_neu * dsqrt(scale_alpha_z_nn)
1738  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
1739  al2=al_qup * dsqrt(scale_alpha_z_uu)
1740  ar2=ar_qup * dsqrt(scale_alpha_z_uu)
1741  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
1742  al2=al_qdn * dsqrt(scale_alpha_z_dd)
1743  ar2=ar_qdn * dsqrt(scale_alpha_z_dd)
1744  else
1745  al2=0d0
1746  ar2=0d0
1747  endif
1748  pv(3,:) = pordered(:,6)+pordered(:,7)
1749  pv(4,:) = pordered(:,8)+pordered(:,9)
1750  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1751  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1752  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1753  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1754  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1755  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
1756  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1757  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
1758 
1759 
1760 
1761  elseif( vvmode.eq.wwmode ) then
1762  ! WW DECAYS
1763  if( isaquark(idordered(6)) ) then
1764  al1 = bl * dsqrt(scale_alpha_w_ud)
1765  ar1 = br * dsqrt(scale_alpha_w_ud)! = 0
1766  elseif( &
1767  (abs(idordered(6)).eq.abs(elp_) .and. abs(idordered(7)).eq.abs(nue_)) .or. (abs(idordered(7)).eq.abs(elp_) .and. abs(idordered(6)).eq.abs(nue_)) .or. &
1768  (abs(idordered(6)).eq.abs(mup_) .and. abs(idordered(7)).eq.abs(num_)) .or. (abs(idordered(7)).eq.abs(mup_) .and. abs(idordered(6)).eq.abs(num_)) &
1769  ) then
1770  al1 = bl * dsqrt(scale_alpha_w_ln)
1771  ar1 = br * dsqrt(scale_alpha_w_ln)! = 0
1772  elseif( &
1773  (abs(idordered(6)).eq.abs(tap_) .and. abs(idordered(7)).eq.abs(nut_)) .or. (abs(idordered(7)).eq.abs(tap_) .and. abs(idordered(6)).eq.abs(nut_)) &
1774  ) then
1775  al1 = bl * dsqrt(scale_alpha_w_tn)
1776  ar1 = br * dsqrt(scale_alpha_w_tn)! = 0
1777  else
1778  al1=0d0
1779  ar1=0d0
1780  endif
1781  if( isaquark(idordered(8)) ) then
1782  al2 = bl * dsqrt(scale_alpha_w_ud)
1783  ar2 = br * dsqrt(scale_alpha_w_ud)! = 0
1784  elseif( &
1785  (abs(idordered(8)).eq.abs(elm_) .and. abs(idordered(9)).eq.abs(anue_)) .or. (abs(idordered(9)).eq.abs(elm_) .and. abs(idordered(8)).eq.abs(anue_)) .or. &
1786  (abs(idordered(8)).eq.abs(mum_) .and. abs(idordered(9)).eq.abs(anum_)) .or. (abs(idordered(9)).eq.abs(mum_) .and. abs(idordered(8)).eq.abs(anum_)) &
1787  ) then
1788  al2 = bl * dsqrt(scale_alpha_w_ln)
1789  ar2 = br * dsqrt(scale_alpha_w_ln)! = 0
1790  elseif( &
1791  (abs(idordered(8)).eq.abs(tam_) .and. abs(idordered(9)).eq.abs(anut_)) .or. (abs(idordered(9)).eq.abs(tam_) .and. abs(idordered(8)).eq.abs(anut_)) &
1792  ) then
1793  al2 = bl * dsqrt(scale_alpha_w_tn)
1794  ar2 = br * dsqrt(scale_alpha_w_tn)! = 0
1795  else
1796  al2=0d0
1797  ar2=0d0
1798  endif
1799  pv(3,:) = pordered(:,6)+pordered(:,7)
1800  pv(4,:) = pordered(:,8)+pordered(:,9)
1801  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1802  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1803  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1804  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1805  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1806  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
1807  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1808  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
1809 
1810 
1811  elseif( vvmode.eq.zgmode ) then
1812  ! Zgamma DECAYS
1813  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1814  al1=al_lep * dsqrt(scale_alpha_z_ll)
1815  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
1816  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1817  al1=al_lep * dsqrt(scale_alpha_z_tt)
1818  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
1819  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1820  al1=al_neu * dsqrt(scale_alpha_z_nn)
1821  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
1822  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1823  al1=al_qup * dsqrt(scale_alpha_z_uu)
1824  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
1825  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1826  al1=al_qdn * dsqrt(scale_alpha_z_dd)
1827  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
1828  else
1829  al1=0d0
1830  ar1=0d0
1831  endif
1832  al2=1d0
1833  ar2=1d0
1834  pv(3,:) = pordered(:,6)+pordered(:,7)
1835  pv(4,:) = pordered(:,8)
1836  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1837  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1838  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
1839  !sp(4,1:4)=pV(4,1:4); print *, "this checks gauge invariance"
1840  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1841  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
1842  propv(2)=1d0
1843 
1844 
1845  elseif( vvmode.eq.ggmode ) then
1846  ! gamma gamma DECAYS
1847  al1=1d0
1848  ar1=1d0
1849  al2=1d0
1850  ar2=1d0
1851  pv(3,:) = pordered(:,6)
1852  pv(4,:) = pordered(:,8)
1853  sp(3,:) = pol_mless2(dcmplx(pordered(:,6)),h3,'out') ! photon
1854  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
1855  !sp(3,1:4)=pV(3,1:4); print *, "this checks gauge invariance"
1856  !sp(4,1:4)=pV(4,1:4)
1857  propv(1)=1d0
1858  propv(2)=1d0
1859 
1860 
1861  elseif( vvmode.eq.gsgmode ) then
1862  ! gamma* gamma DECAYS
1863  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1864  al1=cl_lep * dsqrt(scale_alpha_z_ll)
1865  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
1866  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1867  al1=cl_lep * dsqrt(scale_alpha_z_tt)
1868  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
1869  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1870  al1=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
1871  ar1=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
1872  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1873  al1=cl_qup * dsqrt(scale_alpha_z_uu)
1874  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
1875  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1876  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
1877  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
1878  else
1879  al1=0d0
1880  ar1=0d0
1881  endif
1882  al2=1d0
1883  ar2=1d0
1884  pv(3,:) = pordered(:,6)+pordered(:,7)
1885  pv(4,:) = pordered(:,8)
1886  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1887  sp(3,:) = -sp(3,:) ! photon propagator
1888  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
1889  !sp(4,1:4)=pV(4,1:4); print *, "this checks gauge invariance"
1890  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1891  propv(1) = 1d0
1892  propv(2) = 1d0
1893  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
1894 
1895 
1896  elseif( vvmode.eq.gszmode ) then
1897  ! gamma* Z DECAYS
1898  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1899  al1=cl_lep * dsqrt(scale_alpha_z_ll)
1900  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
1901  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1902  al1=cl_lep * dsqrt(scale_alpha_z_tt)
1903  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
1904  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1905  al1=cl_neu * dsqrt(scale_alpha_z_nn)
1906  ar1=cr_neu * dsqrt(scale_alpha_z_nn)
1907  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1908  al1=cl_qup * dsqrt(scale_alpha_z_uu)
1909  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
1910  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1911  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
1912  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
1913  else
1914  al1=0d0
1915  ar1=0d0
1916  endif
1917  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
1918  al2=al_lep * dsqrt(scale_alpha_z_ll)
1919  ar2=ar_lep * dsqrt(scale_alpha_z_ll)
1920  elseif( abs(idordered(8)).eq.abs(tam_) ) then
1921  al2=al_lep * dsqrt(scale_alpha_z_tt)
1922  ar2=ar_lep * dsqrt(scale_alpha_z_tt)
1923  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
1924  al2=al_neu * dsqrt(scale_alpha_z_nn)
1925  ar2=ar_neu * dsqrt(scale_alpha_z_nn)
1926  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
1927  al2=al_qup * dsqrt(scale_alpha_z_uu)
1928  ar2=ar_qup * dsqrt(scale_alpha_z_uu)
1929  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
1930  al2=al_qdn * dsqrt(scale_alpha_z_dd)
1931  ar2=ar_qdn * dsqrt(scale_alpha_z_dd)
1932  else
1933  al2=0d0
1934  ar2=0d0
1935  endif
1936  pv(3,:) = pordered(:,6)+pordered(:,7)
1937  pv(4,:) = pordered(:,8)+pordered(:,9)
1938  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1939  sp(3,:) = -sp(3,:)
1940  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1941  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1942  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1943  propv(1) = 1d0! = s/dcmplx(s)
1944  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
1945  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1946  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
1947 
1948 
1949  elseif( vvmode.eq.zgsmode ) then
1950  ! Z gamma* DECAYS
1951  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1952  al1=al_lep * dsqrt(scale_alpha_z_ll)
1953  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
1954  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1955  al1=al_lep * dsqrt(scale_alpha_z_tt)
1956  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
1957  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1958  al1=al_neu * dsqrt(scale_alpha_z_nn)
1959  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
1960  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1961  al1=al_qup * dsqrt(scale_alpha_z_uu)
1962  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
1963  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1964  al1=al_qdn * dsqrt(scale_alpha_z_dd)
1965  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
1966  else
1967  al1=0d0
1968  ar1=0d0
1969  endif
1970  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
1971  al2=cl_lep * dsqrt(scale_alpha_z_ll)
1972  ar2=cr_lep * dsqrt(scale_alpha_z_ll)
1973  elseif( abs(idordered(8)).eq.abs(tam_) ) then
1974  al2=cl_lep * dsqrt(scale_alpha_z_tt)
1975  ar2=cr_lep * dsqrt(scale_alpha_z_tt)
1976  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
1977  al2=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
1978  ar2=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
1979  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
1980  al2=cl_qup * dsqrt(scale_alpha_z_uu)
1981  ar2=cr_qup * dsqrt(scale_alpha_z_uu)
1982  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
1983  al2=cl_qdn * dsqrt(scale_alpha_z_dd)
1984  ar2=cr_qdn * dsqrt(scale_alpha_z_dd)
1985  else
1986  al2=0d0
1987  ar2=0d0
1988  endif
1989  pv(3,:) = pordered(:,6)+pordered(:,7)
1990  pv(4,:) = pordered(:,8)+pordered(:,9)
1991  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1992  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1993  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1994  sp(4,:) = -sp(4,:)
1995  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1996  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
1997  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1998  propv(2) = 1d0 ! = s/dcmplx(s)
1999  if( s.lt.mphotoncutoff**2 ) propv(2)=czero
2000 
2001 
2002  elseif( vvmode.eq.gsgsmode ) then
2003  ! gamma* gamma* DECAYS
2004  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
2005  al1=cl_lep * dsqrt(scale_alpha_z_ll)
2006  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
2007  elseif( abs(idordered(6)).eq.abs(tam_) ) then
2008  al1=cl_lep * dsqrt(scale_alpha_z_tt)
2009  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
2010  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
2011  al1=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
2012  ar1=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
2013  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
2014  al1=cl_qup * dsqrt(scale_alpha_z_uu)
2015  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
2016  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
2017  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
2018  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
2019  else
2020  al1=0d0
2021  ar1=0d0
2022  endif
2023  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
2024  al2=cl_lep * dsqrt(scale_alpha_z_ll)
2025  ar2=cr_lep * dsqrt(scale_alpha_z_ll)
2026  elseif( abs(idordered(8)).eq.abs(tam_) ) then
2027  al2=cl_lep * dsqrt(scale_alpha_z_tt)
2028  ar2=cr_lep * dsqrt(scale_alpha_z_tt)
2029  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
2030  al2=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
2031  ar2=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
2032  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
2033  al2=cl_qup * dsqrt(scale_alpha_z_uu)
2034  ar2=cr_qup * dsqrt(scale_alpha_z_uu)
2035  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
2036  al2=cl_qdn * dsqrt(scale_alpha_z_dd)
2037  ar2=cr_qdn * dsqrt(scale_alpha_z_dd)
2038  else
2039  al2=0d0
2040  ar2=0d0
2041  endif
2042  pv(3,:) = pordered(:,6)+pordered(:,7)
2043  pv(4,:) = pordered(:,8)+pordered(:,9)
2044  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2045  sp(3,:) = -sp(3,:)
2046  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2047  sp(4,:) = -sp(4,:)
2048  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
2049  propv(1) = 1d0 ! = s/dcmplx(s)
2050  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
2051  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
2052  propv(2) = 1d0 ! = s/dcmplx(s)
2053  if( s.lt.mphotoncutoff**2 ) propv(2)=czero
2054 
2055 
2056  elseif( vvmode.eq.zpzmode ) then
2057  ! Z'Z DECAYS
2058  al1 = vpffcoupling(idordered(6),-1,.false.)
2059  ar1 = vpffcoupling(idordered(6),+1,.false.)
2060  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
2061  al2=al_lep * dsqrt(scale_alpha_z_ll)
2062  ar2=ar_lep * dsqrt(scale_alpha_z_ll)
2063  elseif( abs(idordered(8)).eq.abs(tam_) ) then
2064  al2=al_lep * dsqrt(scale_alpha_z_tt)
2065  ar2=ar_lep * dsqrt(scale_alpha_z_tt)
2066  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
2067  al2=al_neu * dsqrt(scale_alpha_z_nn)
2068  ar2=ar_neu * dsqrt(scale_alpha_z_nn)
2069  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
2070  al2=al_qup * dsqrt(scale_alpha_z_uu)
2071  ar2=ar_qup * dsqrt(scale_alpha_z_uu)
2072  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
2073  al2=al_qdn * dsqrt(scale_alpha_z_dd)
2074  ar2=ar_qdn * dsqrt(scale_alpha_z_dd)
2075  else
2076  al2=0d0
2077  ar2=0d0
2078  endif
2079 
2080  pv(3,:) = pordered(:,6)+pordered(:,7)
2081  pv(4,:) = pordered(:,8)+pordered(:,9)
2082  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2083  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2084  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2085  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
2086  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
2087  if( m_vprime .gt. 0d0 ) then
2088  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2089  elseif( m_vprime .eq. 0d0 ) then
2090  propv(1) = 1d0
2091  else
2092  propv(1) = s/m_z**2
2093  endif
2094  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
2095  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
2096 
2097 
2098  elseif( vvmode.eq.zzpmode ) then
2099  ! ZZ' DECAYS
2100  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
2101  al1=al_lep * dsqrt(scale_alpha_z_ll)
2102  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
2103  elseif( abs(idordered(6)).eq.abs(tam_) ) then
2104  al1=al_lep * dsqrt(scale_alpha_z_tt)
2105  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
2106  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
2107  al1=al_neu * dsqrt(scale_alpha_z_nn)
2108  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
2109  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
2110  al1=al_qup * dsqrt(scale_alpha_z_uu)
2111  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
2112  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
2113  al1=al_qdn * dsqrt(scale_alpha_z_dd)
2114  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
2115  else
2116  al1=0d0
2117  ar1=0d0
2118  endif
2119  al2 = vpffcoupling(idordered(8),-1,.false.)
2120  ar2 = vpffcoupling(idordered(8),+1,.false.)
2121 
2122  pv(3,:) = pordered(:,6)+pordered(:,7)
2123  pv(4,:) = pordered(:,8)+pordered(:,9)
2124  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2125  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2126  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2127  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
2128  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
2129  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
2130  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
2131  if( m_vprime .gt. 0d0 ) then
2132  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2133  elseif( m_vprime .eq. 0d0 ) then
2134  propv(2) = 1d0
2135  else
2136  propv(2) = s/m_z**2
2137  endif
2138 
2139 
2140  elseif( vvmode.eq.zpzpmode ) then
2141  ! Z'Z' DECAYS
2142  al1 = vpffcoupling(idordered(6),-1,.false.)
2143  ar1 = vpffcoupling(idordered(6),+1,.false.)
2144  al2 = vpffcoupling(idordered(8),-1,.false.)
2145  ar2 = vpffcoupling(idordered(8),+1,.false.)
2146 
2147  pv(3,:) = pordered(:,6)+pordered(:,7)
2148  pv(4,:) = pordered(:,8)+pordered(:,9)
2149  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2150  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2151  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2152  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
2153  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
2154  if( m_vprime .gt. 0d0 ) then
2155  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2156  elseif( m_vprime .eq. 0d0 ) then
2157  propv(1) = 1d0
2158  else
2159  propv(1) = s/m_z**2
2160  endif
2161  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
2162  if( m_vprime .gt. 0d0 ) then
2163  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2164  elseif( m_vprime .eq. 0d0 ) then
2165  propv(2) = 1d0
2166  else
2167  propv(2) = s/m_z**2
2168  endif
2169 
2170 
2171  elseif( vvmode.eq.zpgsmode ) then
2172  ! Z' gamma* DECAYS
2173  al1 = vpffcoupling(idordered(6),-1,.false.)
2174  ar1 = vpffcoupling(idordered(6),+1,.false.)
2175  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
2176  al2=cl_lep * dsqrt(scale_alpha_z_ll)
2177  ar2=cr_lep * dsqrt(scale_alpha_z_ll)
2178  elseif( abs(idordered(8)).eq.abs(tam_) ) then
2179  al2=cl_lep * dsqrt(scale_alpha_z_tt)
2180  ar2=cr_lep * dsqrt(scale_alpha_z_tt)
2181  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
2182  al2=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
2183  ar2=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
2184  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
2185  al2=cl_qup * dsqrt(scale_alpha_z_uu)
2186  ar2=cr_qup * dsqrt(scale_alpha_z_uu)
2187  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
2188  al2=cl_qdn * dsqrt(scale_alpha_z_dd)
2189  ar2=cr_qdn * dsqrt(scale_alpha_z_dd)
2190  else
2191  al2=0d0
2192  ar2=0d0
2193  endif
2194 
2195  pv(3,:) = pordered(:,6)+pordered(:,7)
2196  pv(4,:) = pordered(:,8)+pordered(:,9)
2197  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2198  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2199  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2200  sp(4,:) = -sp(4,:)
2201  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
2202  if( m_vprime .gt. 0d0 ) then
2203  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2204  elseif( m_vprime .eq. 0d0 ) then
2205  propv(1) = 1d0
2206  else
2207  propv(1) = s/m_z**2
2208  endif
2209  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
2210  propv(2) = 1d0 ! = s/dcmplx(s)
2211  if( s.lt.mphotoncutoff**2 ) propv(2)=czero
2212 
2213  elseif( vvmode.eq.gszpmode ) then
2214  ! gamma* Z' DECAYS
2215  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
2216  al1=cl_lep * dsqrt(scale_alpha_z_ll)
2217  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
2218  elseif( abs(idordered(6)).eq.abs(tam_) ) then
2219  al1=cl_lep * dsqrt(scale_alpha_z_tt)
2220  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
2221  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
2222  al1=cl_neu * dsqrt(scale_alpha_z_nn)
2223  ar1=cr_neu * dsqrt(scale_alpha_z_nn)
2224  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
2225  al1=cl_qup * dsqrt(scale_alpha_z_uu)
2226  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
2227  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
2228  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
2229  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
2230  else
2231  al1=0d0
2232  ar1=0d0
2233  endif
2234  al2 = vpffcoupling(idordered(8),-1,.false.)
2235  ar2 = vpffcoupling(idordered(8),+1,.false.)
2236 
2237  pv(3,:) = pordered(:,6)+pordered(:,7)
2238  pv(4,:) = pordered(:,8)+pordered(:,9)
2239  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2240  sp(3,:) = -sp(3,:)
2241  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2242  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
2243  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
2244  propv(1) = 1d0! = s/dcmplx(s)
2245  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
2246  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
2247  if( m_vprime .gt. 0d0 ) then
2248  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2249  elseif( m_vprime .eq. 0d0 ) then
2250  propv(2) = 1d0
2251  else
2252  propv(2) = s/m_z**2
2253  endif
2254 
2255 
2256  elseif( vvmode.eq.zpgmode ) then
2257  ! Z' gamma DECAYS
2258  al1 = vpffcoupling(idordered(6),-1,.false.)
2259  ar1 = vpffcoupling(idordered(6),+1,.false.)
2260  al2=1d0
2261  ar2=1d0
2262  pv(3,:) = pordered(:,6)+pordered(:,7)
2263  pv(4,:) = pordered(:,8)
2264  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2265  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2266  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
2267  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
2268  if( m_vprime .gt. 0d0 ) then
2269  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2270  elseif( m_vprime .eq. 0d0 ) then
2271  propv(1) = 1d0
2272  else
2273  propv(1) = s/m_z**2
2274  endif
2275  propv(2)=1d0
2276 
2277 
2278  elseif( vvmode.eq.wpwmode ) then
2279  ! W'W DECAYS
2280  al1 = vpffcoupling(idordered(6),-1,.true.)
2281  ar1 = vpffcoupling(idordered(6),+1,.true.)
2282  if( isaquark(idordered(8)) ) then
2283  al2 = bl * dsqrt(scale_alpha_w_ud)
2284  ar2 = br * dsqrt(scale_alpha_w_ud)! = 0
2285  elseif( &
2286  (abs(idordered(8)).eq.abs(elm_) .and. abs(idordered(9)).eq.abs(anue_)) .or. (abs(idordered(9)).eq.abs(elm_) .and. abs(idordered(8)).eq.abs(anue_)) .or. &
2287  (abs(idordered(8)).eq.abs(mum_) .and. abs(idordered(9)).eq.abs(anum_)) .or. (abs(idordered(9)).eq.abs(mum_) .and. abs(idordered(8)).eq.abs(anum_)) &
2288  ) then
2289  al2 = bl * dsqrt(scale_alpha_w_ln)
2290  ar2 = br * dsqrt(scale_alpha_w_ln)! = 0
2291  elseif( &
2292  (abs(idordered(8)).eq.abs(tam_) .and. abs(idordered(9)).eq.abs(anut_)) .or. (abs(idordered(9)).eq.abs(tam_) .and. abs(idordered(8)).eq.abs(anut_)) &
2293  ) then
2294  al2 = bl * dsqrt(scale_alpha_w_tn)
2295  ar2 = br * dsqrt(scale_alpha_w_tn)! = 0
2296  else
2297  al2=0d0
2298  ar2=0d0
2299  endif
2300  pv(3,:) = pordered(:,6)+pordered(:,7)
2301  pv(4,:) = pordered(:,8)+pordered(:,9)
2302  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2303  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2304  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2305  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
2306  s = scr(pv(3,:),pv(3,:))
2307  if( m_vprime .gt. 0d0 ) then
2308  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2309  elseif( m_vprime .eq. 0d0 ) then
2310  propv(1) = 1d0
2311  else
2312  propv(1) = s/m_w**2
2313  endif
2314  s = scr(pv(4,:),pv(4,:))
2315  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
2316 
2317 
2318  elseif( vvmode.eq.wwpmode ) then
2319  ! WW' DECAYS
2320  if( isaquark(idordered(6)) ) then
2321  al1 = bl * dsqrt(scale_alpha_w_ud)
2322  ar1 = br * dsqrt(scale_alpha_w_ud)! = 0
2323  elseif( &
2324  (abs(idordered(6)).eq.abs(elp_) .and. abs(idordered(7)).eq.abs(nue_)) .or. (abs(idordered(7)).eq.abs(elp_) .and. abs(idordered(6)).eq.abs(nue_)) .or. &
2325  (abs(idordered(6)).eq.abs(mup_) .and. abs(idordered(7)).eq.abs(num_)) .or. (abs(idordered(7)).eq.abs(mup_) .and. abs(idordered(6)).eq.abs(num_)) &
2326  ) then
2327  al1 = bl * dsqrt(scale_alpha_w_ln)
2328  ar1 = br * dsqrt(scale_alpha_w_ln)! = 0
2329  elseif( &
2330  (abs(idordered(6)).eq.abs(tap_) .and. abs(idordered(7)).eq.abs(nut_)) .or. (abs(idordered(7)).eq.abs(tap_) .and. abs(idordered(6)).eq.abs(nut_)) &
2331  ) then
2332  al1 = bl * dsqrt(scale_alpha_w_tn)
2333  ar1 = br * dsqrt(scale_alpha_w_tn)! = 0
2334  else
2335  al1=0d0
2336  ar1=0d0
2337  endif
2338  al2 = vpffcoupling(idordered(8),-1,.true.)
2339  ar2 = vpffcoupling(idordered(8),+1,.true.)
2340 
2341  pv(3,:) = pordered(:,6)+pordered(:,7)
2342  pv(4,:) = pordered(:,8)+pordered(:,9)
2343  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2344  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2345  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2346  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
2347  s = scr(pv(3,:),pv(3,:))
2348  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
2349  s = scr(pv(4,:),pv(4,:))
2350  if( m_vprime .gt. 0d0 ) then
2351  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2352  elseif( m_vprime .eq. 0d0 ) then
2353  propv(2) = 1d0
2354  else
2355  propv(2) = s/m_w**2
2356  endif
2357 
2358 
2359  elseif( vvmode.eq.wpwpmode ) then
2360  ! W'W' DECAYS
2361  al1 = vpffcoupling(idordered(6),-1,.true.)
2362  ar1 = vpffcoupling(idordered(6),+1,.true.)
2363  al2 = vpffcoupling(idordered(8),-1,.true.)
2364  ar2 = vpffcoupling(idordered(8),+1,.true.)
2365 
2366  pv(3,:) = pordered(:,6)+pordered(:,7)
2367  pv(4,:) = pordered(:,8)+pordered(:,9)
2368  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
2369  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
2370  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
2371  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
2372  s = scr(pv(3,:),pv(3,:))
2373  if( m_vprime .gt. 0d0 ) then
2374  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2375  elseif( m_vprime .eq. 0d0 ) then
2376  propv(1) = 1d0
2377  else
2378  propv(1) = s/m_w**2
2379  endif
2380  s = scr(pv(4,:),pv(4,:))
2381  if( m_vprime .gt. 0d0 ) then
2382  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
2383  elseif( m_vprime .eq. 0d0 ) then
2384  propv(2) = 1d0
2385  else
2386  propv(2) = s/m_w**2
2387  endif
2388 
2389 
2390  else
2391  call error("Unsupported decay modes")
2392  endif
2393 
2394  sp(3,:) = sp(3,:)*propv(1)
2395  sp(4,:) = sp(4,:)*propv(2)
2396  if (h3.eq.-1) then
2397  sp(3,:) = al1 * sp(3,:)
2398  elseif(h3.eq.1) then
2399  sp(3,:) = ar1 * sp(3,:)
2400  endif
2401  if (h4.eq.-1) then
2402  sp(4,:) = al2 * sp(4,:)
2403  elseif(h4.eq.1) then
2404  sp(4,:) = ar2 * sp(4,:)
2405  endif
2406 
2407  return
2408 end subroutine
2409 
2410 subroutine getdecay_vvmode_ordering(MY_IDUP, VVMode,ordering,VVMode_swap,ordering_swap)
2411  implicit none
2412  integer, intent(in) :: MY_IDUP(6:9)
2413  integer, intent(out) :: VVMode,ordering(1:4),VVMode_swap,ordering_swap(1:4)
2414  integer :: idV(1:2),idV_swap(1:2)
2415 
2416  ordering=(/3,4,5,6/)
2417  idv(1)=coupledvertex(my_idup(6:7),-1)
2418  idv(2)=coupledvertex(my_idup(8:9),-1)
2419  if(my_idup(6).eq.pho_ .or. my_idup(7).eq.pho_) idv(1)=pho_
2420  if(my_idup(8).eq.pho_ .or. my_idup(9).eq.pho_) idv(2)=pho_
2421  if(convertlhe(my_idup(6)).lt.0 .or. my_idup(6).eq.not_a_particle_) then
2422  call swap(ordering(1),ordering(2))
2423  endif
2424  if(convertlhe(my_idup(8)).lt.0 .or. my_idup(8).eq.not_a_particle_) then
2425  call swap(ordering(3),ordering(4))
2426  endif
2427  if( &
2428  (idv(1).eq.wm_ .and. idv(2).eq.wp_) .or. &
2429  (idv(2).eq.z0_ .and. idv(1).eq.pho_) &
2430  ) then
2431  call swap(ordering(1),ordering(3))
2432  call swap(ordering(2),ordering(4))
2433  call swap(idv(1),idv(2))
2434  endif
2435  ordering_swap(:)=ordering(:)
2436  call swap(ordering_swap(1),ordering_swap(3))
2437 
2438  idv_swap(1) = coupledvertex( (/ my_idup(3+ordering_swap(1)), my_idup(3+ordering_swap(2)) /), -1)
2439  idv_swap(2) = coupledvertex( (/ my_idup(3+ordering_swap(3)), my_idup(3+ordering_swap(4)) /), -1)
2440  if ( (idv_swap(1).eq.wm_) .and. (idv_swap(2).eq.wp_) ) then
2441  call swap(ordering_swap(1),ordering_swap(3))
2442  call swap(ordering_swap(2),ordering_swap(4))
2443  call swap(idv_swap(1),idv_swap(2))
2444  endif
2445 
2446  if(idv(1).eq.z0_ .and. idv(2).eq.z0_) then
2447  vvmode=zzmode
2448  elseif(idv(1).eq.z0_ .and. idv(2).eq.pho_) then
2449  vvmode=zgmode
2450  elseif(idv(1).eq.pho_ .and. idv(2).eq.pho_) then
2451  vvmode=ggmode
2452  elseif(idv(1).eq.wp_ .and. idv(2).eq.wm_) then
2453  vvmode=wwmode
2454  else
2455  print *,"idV=",idv
2456  call error("Unsupported decay mode")
2457  endif
2458 
2459  vvmode_swap=invalidmode
2460  if(idv_swap(1).eq.z0_ .and. idv_swap(2).eq.z0_) then
2461  vvmode_swap=zzmode
2462  elseif(idv_swap(1).eq.wp_ .and. idv_swap(2).eq.wm_) then
2463  vvmode_swap=wwmode
2464  endif
2465  return
2466 end subroutine
2467 
2468 
2469  end module
2470 
2471 
modparameters::b5
complex(8), public b5
Definition: mod_Parameters.F90:940
modparameters::b1
complex(8), public b1
Definition: mod_Parameters.F90:936
modparameters::bzpgs8
complex(8), public bzpgs8
Definition: mod_Parameters.F90:985
modparameters::m_vprime
real(8), public m_vprime
Definition: mod_Parameters.F90:78
modparameters::includegammastar
logical, public includegammastar
Definition: mod_Parameters.F90:213
modparameters::bzpgs3
complex(8), public bzpgs3
Definition: mod_Parameters.F90:983
modparameters::a5
complex(8), public a5
Definition: mod_Parameters.F90:927
modgraviton::getdecay_couplings_spinors_props
subroutine getdecay_couplings_spinors_props(VVMode, idordered, pordered, h3
Definition: mod_Graviton.F90:1699
modparameters::elm_
integer, target, public elm_
Definition: mod_Parameters.F90:1112
modparameters::zpgsmode
integer, parameter, public zpgsmode
Definition: mod_Parameters.F90:16
modparameters::mphotoncutoff
real(8), public mphotoncutoff
Definition: mod_Parameters.F90:215
modparameters::bzpzp1
complex(8), public bzpzp1
Definition: mod_Parameters.F90:958
modparameters::coupledvertex
integer function coupledvertex(id, hel, useAHcoupl)
Definition: mod_Parameters.F90:2467
modparameters::c7
complex(8), parameter, public c7
Definition: mod_Parameters.F90:995
modparameters::bzpzp6
complex(8), public bzpzp6
Definition: mod_Parameters.F90:963
modparameters::zgmode
integer, parameter, public zgmode
Definition: mod_Parameters.F90:13
modparameters::bzgs1
complex(8), public bzgs1
Definition: mod_Parameters.F90:969
modparameters::cl_neu
real(8), public cl_neu
Definition: mod_Parameters.F90:1071
modparameters::b10
complex(8), public b10
Definition: mod_Parameters.F90:945
modparameters::b8
complex(8), public b8
Definition: mod_Parameters.F90:943
modparameters::bzpgs2
complex(8), public bzpgs2
Definition: mod_Parameters.F90:982
modmisc::et1
double complex function et1(e1, e2, e3, e4)
Definition: mod_Misc.F90:109
modparameters::a2
complex(8), public a2
Definition: mod_Parameters.F90:924
modparameters::bzpzp8
complex(8), public bzpzp8
Definition: mod_Parameters.F90:965
modparameters::dn_
integer, target, public dn_
Definition: mod_Parameters.F90:1085
modparameters::bzgs4
complex(8), public bzgs4
Definition: mod_Parameters.F90:972
modparameters::bzzp2
complex(8), public bzzp2
Definition: mod_Parameters.F90:948
modparameters::pol_mless2
complex(dp) function, dimension(4) pol_mless2(p, i, out)
Definition: mod_Parameters.F90:3053
modgraviton::evalamp_gg_g_vv
subroutine, public evalamp_gg_g_vv(p, MY_IDUP, res)
Definition: mod_Graviton.F90:18
modparameters::scale_alpha_w_tn
real(8), public scale_alpha_w_tn
Definition: mod_Parameters.F90:335
modparameters::zpzpmode
integer, parameter, public zpzpmode
Definition: mod_Parameters.F90:15
modgraviton::calchelamp_gg
subroutine calchelamp_gg(ordering, VVMode, p, MY_IDUP, i1, i2, i3, i4, A)
Definition: mod_Graviton.F90:121
modparameters::one
real(8), parameter, public one
Definition: mod_Parameters.F90:83
modparameters::zgsmode
integer, parameter, public zgsmode
Definition: mod_Parameters.F90:13
modparameters::sqrt2
real(dp), parameter, public sqrt2
Definition: mod_Parameters.F90:80
modparameters::scale_alpha_w_ln
real(8), public scale_alpha_w_ln
Definition: mod_Parameters.F90:334
modparameters::bzpzp3
complex(8), public bzpzp3
Definition: mod_Parameters.F90:960
modparameters::two
real(8), parameter, public two
Definition: mod_Parameters.F90:84
modmisc::error
subroutine error(Message, ErrNum)
Definition: mod_Misc.F90:366
modparameters::includevprime
logical, public includevprime
Definition: mod_Parameters.F90:214
modparameters::c6
complex(8), parameter, public c6
Definition: mod_Parameters.F90:994
modparameters::lambda
real(8), parameter, public lambda
Definition: mod_Parameters.F90:354
modparameters::bgsgs1
complex(8), public bgsgs1
Definition: mod_Parameters.F90:975
modparameters::bgsgs8
complex(8), public bgsgs8
Definition: mod_Parameters.F90:979
modparameters::graviton_qq_left
complex(8), public graviton_qq_left
Definition: mod_Parameters.F90:928
modparameters::cr_qdn
real(8), public cr_qdn
Definition: mod_Parameters.F90:1074
modparameters::bzzp10
complex(8), public bzzp10
Definition: mod_Parameters.F90:956
modparameters::bgsgs4
complex(8), public bgsgs4
Definition: mod_Parameters.F90:978
modparameters::c42
complex(8), parameter, public c42
Definition: mod_Parameters.F90:992
modparameters::pol_mass2
complex(dp) function, dimension(4) pol_mass2(p, i, out)
Definition: mod_Parameters.F90:3176
modmisc::sc
double complex function sc(q1, q2)
Definition: mod_Misc.F90:127
modparameters::cr_neu
real(8), public cr_neu
Definition: mod_Parameters.F90:1070
modparameters::wwpmode
integer, parameter, public wwpmode
Definition: mod_Parameters.F90:14
modparameters::wm_
integer, target, public wm_
Definition: mod_Parameters.F90:1115
modparameters::zero
real(8), parameter, public zero
Definition: mod_Parameters.F90:85
modparameters::chm_
integer, target, public chm_
Definition: mod_Parameters.F90:1086
modparameters::up_
integer, target, public up_
Definition: mod_Parameters.F90:1084
modparameters::bot_
integer, target, public bot_
Definition: mod_Parameters.F90:1089
modparameters::cl_lep
real(8), public cl_lep
Definition: mod_Parameters.F90:1069
modparameters::cr_lep
real(8), public cr_lep
Definition: mod_Parameters.F90:1068
modparameters::a3
complex(8), public a3
Definition: mod_Parameters.F90:925
modparameters::ga_reso
real(8), public ga_reso
Definition: mod_Parameters.F90:231
modparameters::bzzp1
complex(8), public bzzp1
Definition: mod_Parameters.F90:947
modparameters::zzmode
integer, parameter, public zzmode
Definition: mod_Parameters.F90:13
modgraviton
Definition: mod_Graviton.F90:1
modparameters::b9
complex(8), public b9
Definition: mod_Parameters.F90:944
modparameters::ar_qdn
real(8), public ar_qdn
Definition: mod_Parameters.F90:1064
modparameters::bzpzp9
complex(8), public bzpzp9
Definition: mod_Parameters.F90:966
modparameters::c41
complex(8), parameter, public c41
Definition: mod_Parameters.F90:991
modparameters::bzgs2
complex(8), public bzgs2
Definition: mod_Parameters.F90:970
modparameters::bzgs8
complex(8), public bzgs8
Definition: mod_Parameters.F90:973
modparameters::gsgmode
integer, parameter, public gsgmode
Definition: mod_Parameters.F90:13
modparameters::scale_alpha_z_ll
real(8), public scale_alpha_z_ll
Definition: mod_Parameters.F90:329
modparameters::bzpzp5
complex(8), public bzpzp5
Definition: mod_Parameters.F90:962
modparameters::bl
real(8), public bl
Definition: mod_Parameters.F90:1066
modparameters::includeinterference
logical, public includeinterference
Definition: mod_Parameters.F90:77
modparameters::generate_bis
logical, parameter, public generate_bis
Definition: mod_Parameters.F90:932
modparameters::ar_neu
real(8), public ar_neu
Definition: mod_Parameters.F90:1060
modparameters::tam_
integer, target, public tam_
Definition: mod_Parameters.F90:1114
modparameters::scale_alpha_w_ud
real(8), public scale_alpha_w_ud
Definition: mod_Parameters.F90:332
modparameters::c5
complex(8), parameter, public c5
Definition: mod_Parameters.F90:993
modparameters::bzzp5
complex(8), public bzzp5
Definition: mod_Parameters.F90:951
modparameters::bzpgs4
complex(8), public bzpgs4
Definition: mod_Parameters.F90:984
modparameters::al_qdn
real(8), public al_qdn
Definition: mod_Parameters.F90:1065
modgraviton::calchelamp_qq
subroutine calchelamp_qq(ordering, VVMode, p, MY_IDUP, i1, i3, i4, A)
Definition: mod_Graviton.F90:952
modparameters::bzzp3
complex(8), public bzzp3
Definition: mod_Parameters.F90:949
modparameters::cl_qdn
real(8), public cl_qdn
Definition: mod_Parameters.F90:1075
modparameters::a4
complex(8), public a4
Definition: mod_Parameters.F90:926
modparameters::nue_
integer, target, public nue_
Definition: mod_Parameters.F90:1097
modparameters::bzgs3
complex(8), public bzgs3
Definition: mod_Parameters.F90:971
modparameters::ar_lep
real(8), public ar_lep
Definition: mod_Parameters.F90:1058
modparameters::b7
complex(8), public b7
Definition: mod_Parameters.F90:942
modgraviton::gzzampl
subroutine gzzampl(VVMode, p, sp, i1, res)
Definition: mod_Graviton.F90:1394
modparameters::bzzp8
complex(8), public bzzp8
Definition: mod_Parameters.F90:954
modparameters::zzpmode
integer, parameter, public zzpmode
Definition: mod_Parameters.F90:15
modparameters::czero
complex(8), parameter, public czero
Definition: mod_Parameters.F90:86
modparameters::gszmode
integer, parameter, public gszmode
Definition: mod_Parameters.F90:13
modparameters::bzpzp4
complex(8), public bzpzp4
Definition: mod_Parameters.F90:961
modparameters::c3
complex(8), parameter, public c3
Definition: mod_Parameters.F90:990
modparameters::elp_
integer, target, public elp_
Definition: mod_Parameters.F90:1090
modparameters::zpzmode
integer, parameter, public zpzmode
Definition: mod_Parameters.F90:15
modparameters::bzzp7
complex(8), public bzzp7
Definition: mod_Parameters.F90:953
modparameters::m_w
real(8), public m_w
Definition: mod_Parameters.F90:228
modparameters::bzpzp7
complex(8), public bzpzp7
Definition: mod_Parameters.F90:964
modparameters::scale_alpha_z_tt
real(8), public scale_alpha_z_tt
Definition: mod_Parameters.F90:330
modparameters::m_v
real(8), public m_v
Definition: mod_Parameters.F90:78
modparameters::anum_
integer, target, public anum_
Definition: mod_Parameters.F90:1117
modparameters::bzzp4
complex(8), public bzzp4
Definition: mod_Parameters.F90:950
modparameters::al_qup
real(8), public al_qup
Definition: mod_Parameters.F90:1063
modparameters::bzzp6
complex(8), public bzzp6
Definition: mod_Parameters.F90:952
modparameters::ga_v
real(8), public ga_v
Definition: mod_Parameters.F90:78
modparameters::anut_
integer, target, public anut_
Definition: mod_Parameters.F90:1118
modparameters::pol_dk2mom
complex(dp) function, dimension(4) pol_dk2mom(plepton, antilepton, i, outgoing)
Definition: mod_Parameters.F90:3068
modparameters::anue_
integer, target, public anue_
Definition: mod_Parameters.F90:1116
modgraviton::calchelamp2
subroutine calchelamp2(ordering, VVMode, p, MY_IDUP, i1, i3, i4, A)
Definition: mod_Graviton.F90:1366
modparameters
Definition: mod_Parameters.F90:1
modparameters::scale_alpha_z_nn
real(8), public scale_alpha_z_nn
Definition: mod_Parameters.F90:331
modparameters::wpwmode
integer, parameter, public wpwmode
Definition: mod_Parameters.F90:14
modmisc
Definition: mod_Misc.F90:1
modgraviton::qqgzzampl
subroutine qqgzzampl(VVMode, p, sp, res)
Definition: mod_Graviton.F90:994
modparameters::bzpzp10
complex(8), public bzpzp10
Definition: mod_Parameters.F90:967
modparameters::convertlhe
integer function convertlhe(Part)
Definition: mod_Parameters.F90:1798
modparameters::graviton_qq_right
complex(8), public graviton_qq_right
Definition: mod_Parameters.F90:929
modparameters::cr_qup
real(8), public cr_qup
Definition: mod_Parameters.F90:1072
modparameters::bzpzp2
complex(8), public bzpzp2
Definition: mod_Parameters.F90:959
modparameters::a1
complex(8), public a1
Definition: mod_Parameters.F90:923
modparameters::scale_alpha_z_dd
real(8), public scale_alpha_z_dd
Definition: mod_Parameters.F90:328
modparameters::symmfac
real(8), parameter, public symmfac
Definition: mod_Parameters.F90:91
modparameters::b2
complex(8), public b2
Definition: mod_Parameters.F90:937
modparameters::vpffcoupling
complex(8) function vpffcoupling(jhuid, hel, useWp)
Definition: mod_Parameters.F90:1370
modparameters::br
real(8), public br
Definition: mod_Parameters.F90:1067
modparameters::m_z
real(8), public m_z
Definition: mod_Parameters.F90:226
modparameters::nut_
integer, target, public nut_
Definition: mod_Parameters.F90:1099
modparameters::wwmode
integer, parameter, public wwmode
Definition: mod_Parameters.F90:13
modparameters::dp
integer, parameter, public dp
Definition: mod_Parameters.F90:11
modparameters::b6
complex(8), public b6
Definition: mod_Parameters.F90:941
modparameters::mum_
integer, target, public mum_
Definition: mod_Parameters.F90:1113
modparameters::b3
complex(8), public b3
Definition: mod_Parameters.F90:938
modparameters::c1
complex(8), parameter, public c1
Definition: mod_Parameters.F90:988
modparameters::wp_
integer, target, public wp_
Definition: mod_Parameters.F90:1096
modmisc::scr
double precision function scr(p1, p2)
Definition: mod_Misc.F90:135
modgraviton::gggzzampl
subroutine gggzzampl(VVMode, p, sp, res)
Definition: mod_Graviton.F90:156
modparameters::m_reso
real(8), public m_reso
Definition: mod_Parameters.F90:230
modparameters::str_
integer, target, public str_
Definition: mod_Parameters.F90:1087
modparameters::isaquark
logical function isaquark(PartType)
Definition: mod_Parameters.F90:2369
modparameters::bzpgs1
complex(8), public bzpgs1
Definition: mod_Parameters.F90:981
modparameters::scale_alpha_z_uu
real(8), public scale_alpha_z_uu
Definition: mod_Parameters.F90:327
modparameters::zpgmode
integer, parameter, public zpgmode
Definition: mod_Parameters.F90:16
modparameters::not_a_particle_
integer, parameter, public not_a_particle_
Definition: mod_Parameters.F90:1121
modparameters::bgsgs2
complex(8), public bgsgs2
Definition: mod_Parameters.F90:976
modparameters::b4
complex(8), public b4
Definition: mod_Parameters.F90:939
modparameters::ggmode
integer, parameter, public ggmode
Definition: mod_Parameters.F90:13
modparameters::al_lep
real(8), public al_lep
Definition: mod_Parameters.F90:1059
modparameters::num_
integer, target, public num_
Definition: mod_Parameters.F90:1098
modparameters::wpwpmode
integer, parameter, public wpwpmode
Definition: mod_Parameters.F90:14
modparameters::use_dynamic_mg
logical, parameter, public use_dynamic_mg
Definition: mod_Parameters.F90:933
modparameters::cl_qup
real(8), public cl_qup
Definition: mod_Parameters.F90:1073
modparameters::gszpmode
integer, parameter, public gszpmode
Definition: mod_Parameters.F90:16
modparameters::tap_
integer, target, public tap_
Definition: mod_Parameters.F90:1092
modparameters::bzzp9
complex(8), public bzzp9
Definition: mod_Parameters.F90:955
modparameters::invalidmode
integer, parameter, public invalidmode
Definition: mod_Parameters.F90:13
modparameters::gsgsmode
integer, parameter, public gsgsmode
Definition: mod_Parameters.F90:13
modparameters::overallcouplvffsq
real(8), public overallcouplvffsq
Definition: mod_Parameters.F90:1057
modparameters::ga_vprime
real(8), public ga_vprime
Definition: mod_Parameters.F90:78
modparameters::pho_
integer, target, public pho_
Definition: mod_Parameters.F90:1094
modgraviton::getdecay_vvmode_ordering
subroutine getdecay_vvmode_ordering(MY_IDUP, VVMode, ordering, VVMode_swap
Definition: mod_Graviton.F90:2411
modparameters::z0_
integer, target, public z0_
Definition: mod_Parameters.F90:1095
modparameters::bgsgs3
complex(8), public bgsgs3
Definition: mod_Parameters.F90:977
modparameters::c2
complex(8), parameter, public c2
Definition: mod_Parameters.F90:989
modmisc::swap
Definition: mod_Misc.F90:5
modparameters::al_neu
real(8), public al_neu
Definition: mod_Parameters.F90:1061
modparameters::mup_
integer, target, public mup_
Definition: mod_Parameters.F90:1091
modparameters::ar_qup
real(8), public ar_qup
Definition: mod_Parameters.F90:1062