JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
Functions/Subroutines
modhiggs Module Reference

Functions/Subroutines

subroutine, public evalamp_gg_h_vv (p, MY_IDUP, res)
 
subroutine, public calchelamp (ordering, VVMode, MY_IDUP, p, i1, i2, i3, i4, A)
 
subroutine gghzzampl (VVMode, p, sp, res)
 
subroutine, public evalamp_h_vv (p, MY_IDUP, res)
 
subroutine, public calchelamp2 (ordering, VVMode, MY_IDUP, p, i3, i4, A)
 
subroutine hzzampl (VVMode, p, sp, res)
 
subroutine, public evalamp_h_ff (pin, mass_F, res)
 
subroutine, public evalamp_h_tt_decay (pin, mass_F, ga_F, res)
 
subroutine getdecay_couplings_spinors_props (VVMode, idordered, pordered, h3
 
subroutine getdecay_vvmode_ordering (MY_IDUP, VVMode, ordering, VVMode_swap
 

Function/Subroutine Documentation

◆ calchelamp()

subroutine, public modhiggs::calchelamp ( integer, dimension(1:4)  ordering,
integer  VVMode,
integer, dimension(6: real(dp) :: p(1:4,1:6) complex(dp) :: propg real(dp) :: pin(4,4) complex(dp) :: a(1:1), sp(4,4) l1=ordering(1) l2=ordering(2) l3=ordering(3) l4=ordering(4) !print *,"l1=",l1 !print *,"l2=",l2 !print *,"l3=",l3 !print *,"l4=",l4 !print *,"id(6)=",convertlhe(my_idup(l1)) !print *,"id(7)=",convertlhe(my_idup(l2)) !print *,"id(8)=",convertlhe(my_idup(l3)) !print *,"id(9)=",convertlhe(my_idup(l4)) !print *,"p(6)=",(p(:,l1)) !print *,"p(7)=",(p(:,l2)) !print *,"p(8)=",(p(:,l3)) !print *,"p(9)=",(p(:,l4)) propg = one/dcmplx(2d0*scr(p(:,1),p(:,2)) - m_reso**2,m_reso*ga_re pin(1,:) = p(:,1) pin(2,:) = p(:,2) sp(1,:) = pol_mless2(dcmplx(p(:,1)),-3+2*i1,'in')  MY_IDUP,
  p,
integer  i1,
integer  i2,
integer  i3,
integer  i4,
  A 
)

Definition at line 168 of file mod_Higgs.F90.

168  implicit none
169  integer :: ordering(1:4),VVMode,i1,i2,i3,i4,l1,l2,l3,l4,MY_IDUP(6:9)
170  real(dp) :: p(1:4,1:6)
171  complex(dp) :: propG
172  real(dp) :: pin(4,4)
173  complex(dp) :: A(1:1), sp(4,4)
174 
175  l1=ordering(1)
176  l2=ordering(2)
177  l3=ordering(3)
178  l4=ordering(4)
179 
180  !print *,"l1=",l1
181  !print *,"l2=",l2

◆ calchelamp2()

subroutine, public modhiggs::calchelamp2 ( integer, dimension(1:4)  ordering,
integer  VVMode,
integer, dimension(6:9)  MY_IDUP,
real(dp), dimension(1:4,1:6)  p,
integer  i3,
integer  i4,
complex(dp), dimension(1:1)  A 
)

Definition at line 482 of file mod_Higgs.F90.

482 
483 
484 ! MADGRAPH CHECK
485 ! call coupsm(0)
486 ! if( (MY_IDUP(6).eq.MY_IDUP(8)) .and. (MY_IDUP(7).eq.MY_IDUP(9)) ) then
487 ! call SH_EMEPEMEP((/-P(1:4,1)-P(1:4,2),P(1:4,3),P(1:4,4),P(1:4,5),P(1:4,6)/)/GeV,res2)
488 ! else
489 ! call SH_EMEPEMEP_NOINT((/-P(1:4,1)-P(1:4,2),P(1:4,3),P(1:4,4),P(1:4,5),P(1:4,6)/)/GeV,res2)
490 ! endif
491 ! res2=res2* cdabs( (0d0,1d0)/dcmplx(2d0*scr(p(:,1),p(:,2))-M_Reso**2,M_Reso*Ga_Reso) * dconjg((0d0,1d0)/dcmplx(2d0*scr(p(:,1),p(:,2))-M_Reso**2,M_Reso*Ga_Reso)) )
492 ! res2=res2*GeV**4
493 ! pause
494 ! res=res2; RETURN
495 
496 ! res= res*prefactor
497 ! print *, "checker 1",res
498 ! print *, "checker 2",res2
499 ! print *, "checker 1/2",res/res2
500 ! pause
501 
502 

◆ evalamp_gg_h_vv()

subroutine, public modhiggs::evalamp_gg_h_vv ( real(dp), dimension(4,6), intent(in)  p,
integer, dimension(6:9), intent(in)  MY_IDUP,
real(dp), intent(out)  res 
)

Definition at line 21 of file mod_Higgs.F90.

21  implicit none
22  real(dp), intent(out) :: res
23  real(dp), intent(in) :: p(4,6)
24  integer, intent(in) :: MY_IDUP(6:9)
25  complex(dp) :: A_VV(1:18), A0_VV(1:2)
26  integer :: i1,i2,i3,i4,VVMode,VVmode_swap
27  real(dp) :: prefactor!,res2
28  real(dp) :: intcolfac
29  integer :: ordering(1:4),ordering_swap(1:4)
30  logical :: doInterference
31 
32  if(isaquark(my_idup(6)) .and. isaquark(my_idup(8))) then
33  intcolfac=1.0_dp/3.0_dp
34  else
35  intcolfac=1.0_dp
36  endif
37 
38  call getdecay_vvmode_ordering(my_idup(6:9),vvmode,ordering,vvmode_swap,ordering_swap)
39 
40 ! Global normalization
41  if( vvmode.eq.zzmode ) then! Z decay
42  prefactor = 8d0*overallcouplvffsq**2
43  elseif( vvmode.eq.wwmode ) then ! W decay
44  prefactor = 8d0*overallcouplvffsq**2
45  elseif( vvmode.eq.zgmode ) then ! Z+photon "decay"
46  prefactor = 8d0*overallcouplvffsq ! Only single powers
47  elseif( vvmode.eq.ggmode ) then ! photon "decay"
48  prefactor = 8d0
49  else
50  prefactor = 0d0
51  endif
52  prefactor = prefactor * (alphas/(3.0_dp*pi*vev))**2
53 
54 
55 ! ! MADGRAPH CHECK
56 ! res=0d0
57 ! if (MY_IDUP(6).ne.MY_IDUP(8) ) return
58 ! if (MY_IDUP(7).ne.MY_IDUP(9) ) return
59 ! if (MY_IDUP(6).eq.MY_IDUP(8) ) return
60 ! if (MY_IDUP(7).eq.MY_IDUP(9) ) return
61 ! print *, "MY COUPL",al1*dsqrt(overallCouplVffsq),ar1*dsqrt(overallCouplVffsq)
62 ! ar1=-0.158480099490745d0! this is MadGraphs gzl(R) , the MG couplings differ by a global minus sign, relative differences are because of different input parameters
63 ! al1=+0.210150647402957d0! this is MadGraphs gzl(L)
64 ! al2=al1
65 ! ar2=ar1
66 ! print *, "MG COUPL",al1,ar1
67 ! pause
68 
69  res = zero
70  a_vv(:) = 0d0
71  dointerference = includeinterference .and. ( &
72  ((vvmode.eq.zzmode) .and. (vvmode_swap.eq.zzmode)) &
73  )
74  if ( includevprime .and. .not.(vvmode.eq.zzmode .or. vvmode.eq.zgmode .or. vvmode.eq.wwmode) ) then
75  call error("Contact terms only for WW, ZZ or Zg!")
76  endif
77  do i1=1,2; do i2=1,2; do i3=1,2; do i4=1,2! sum over helicities
78  call calchelamp(ordering,vvmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(1))
79  if( vvmode.eq.zzmode ) then
80  if( includegammastar ) then
81  call calchelamp(ordering,zgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(3))
82  call calchelamp(ordering,gszmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(5))
83  call calchelamp(ordering,gsgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(7))
84  endif
85  if( includevprime ) then
86  call calchelamp(ordering,zzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(9))
87  call calchelamp(ordering,zpzmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(11))
88  call calchelamp(ordering,zpzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(13))
89  endif
90  if( includegammastar .and. includevprime ) then
91  call calchelamp(ordering,gszpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(15))
92  call calchelamp(ordering,zpgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(17))
93  endif
94  elseif( vvmode.eq.zgmode ) then
95  if(includegammastar) then
96  call calchelamp(ordering,gsgmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(3))
97  endif
98  if( includevprime ) then
99  call calchelamp(ordering,zpgmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(5))
100  endif
101  elseif( vvmode.eq.wwmode .and. includevprime ) then
102  call calchelamp(ordering,wwpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(9))
103  call calchelamp(ordering,wpwmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(11))
104  call calchelamp(ordering,wpwpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(13))
105  endif
106 
107  if( dointerference ) then
108  call calchelamp(ordering_swap,vvmode_swap,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(2))
109  if( includegammastar ) then
110  call calchelamp(ordering_swap,zgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(4))
111  call calchelamp(ordering_swap,gszmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(6))
112  call calchelamp(ordering_swap,gsgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(8))
113  endif
114  if( includevprime ) then
115  call calchelamp(ordering_swap,zzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(10))
116  call calchelamp(ordering_swap,zpzmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(12))
117  call calchelamp(ordering_swap,zpzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(14))
118  endif
119  if( includegammastar .and. includevprime ) then
120  call calchelamp(ordering_swap,gszpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(16))
121  call calchelamp(ordering_swap,zpgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(18))
122  endif
123  endif
124 
125  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
126  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
127  res = res + dreal(a0_vv(1)*dconjg(a0_vv(1)))
128  res = res + dreal(a0_vv(2)*dconjg(a0_vv(2)))
129  if( dointerference .and. (i3.eq.i4) ) then! interfere the 3456 with 5436 pieces
130  res = res - 2d0*intcolfac*dreal( a0_vv(1)*dconjg(a0_vv(2)) ) ! minus from Fermi statistics
131  endif
132  enddo; enddo; enddo; enddo
133 
134 
135 ! MADGRAPH CHECK
136 ! call coupsm(0)
137 ! if( (MY_IDUP(6).eq.MY_IDUP(8)) .and. (MY_IDUP(7).eq.MY_IDUP(9)) ) then
138 ! call SH_EMEPEMEP((/-P(1:4,1)-P(1:4,2),P(1:4,3),P(1:4,4),P(1:4,5),P(1:4,6)/)/GeV,res2)
139 ! else
140 ! call SH_EMEPEMEP_NOINT((/-P(1:4,1)-P(1:4,2),P(1:4,3),P(1:4,4),P(1:4,5),P(1:4,6)/)/GeV,res2)
141 ! call SH_TAMTAPTAMTAP_NOINT((/-P(1:4,1)-P(1:4,2),P(1:4,3),P(1:4,4),P(1:4,5),P(1:4,6)/)/GeV,res2)
142 ! endif
143 ! res2=res2* cdabs( (0d0,1d0)/dcmplx(2d0*scr(p(:,1),p(:,2))-M_Reso**2,M_Reso*Ga_Reso) * dconjg((0d0,1d0)/dcmplx(2d0*scr(p(:,1),p(:,2))-M_Reso**2,M_Reso*Ga_Reso)) )
144 ! res2=res2*GeV**4
145 ! pause
146 ! res=res2; RETURN
147 
148 ! res= res*prefactor
149 ! print *, "checker 1",res
150 ! print *, "checker 2",res2
151 ! print *, "checker 1/2",res/res2
152 ! pause
153 
154 
155  res = res*prefactor
156  if( (vvmode.eq.zzmode) .and. dointerference ) res = res * symmfac
157 
158  !print *,"VVmode:",VVmode
159  !print *,"VVmode_swap:",VVmode_swap
160  !print *,"ids:",MY_IDUP
161  !print *,"res:",res
162  !pause
163 
164  RETURN

◆ evalamp_h_ff()

subroutine, public modhiggs::evalamp_h_ff ( real(dp), dimension(1:4,1:2), intent(in)  pin,
real(dp), intent(in)  mass_F,
real(dp), intent(out)  res 
)

Definition at line 641 of file mod_Higgs.F90.

641  yyy3 = ahz3
642  elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.gszmode) .OR. (vvmode.eq.zgsmode) ) then! decay (Z-photon) OR (gamma*-Z) OR (Z-gamma*)
643  yyy1 = ahz1
644  yyy2 = -2*ahz1*q_q/(q_q-q3_q3)
645  yyy3 = ahz3
646  else
647  yyy1=czero
648  yyy2=czero
649  yyy3=czero
650  print *,"VVMode",vvmode,"not implemented in generate_as"
651  endif
652  endif

◆ evalamp_h_tt_decay()

subroutine, public modhiggs::evalamp_h_tt_decay ( real(dp), dimension(1:4,1:6), intent(in)  pin,
real(dp), intent(in)  mass_F,
real(dp), intent(in)  ga_F,
real(dp), intent(out)  res 
)

Definition at line 664 of file mod_Higgs.F90.

664 ! or H --> tbar(p1) + top(p2)
665 ! with tau/top controlled by value of mass_F
666 ! Since this is an isotropic scalar decay, just provide a matrix element SQUARED
667 ! R.Rontsch July 2015
668  SUBROUTINE evalamp_h_ff(pin,mass_F,res)
669  implicit none
670  real(dp), intent(out) :: res
671  complex(dp) :: amp2
672  real(dp), intent(in) :: pin(1:4,1:2),mass_F
673  real(dp) :: s12
674 
675  s12=2d0*(pin(1,1)*pin(1,2)-pin(2,1)*pin(2,2)-pin(3,1)*pin(3,2)-pin(4,1)*pin(4,2)) + 2d0*mass_f**2
676  amp2 = 2d0*s12*(kappa_tilde**2 + kappa**2) - 8d0*mass_f**2*kappa**2
677  amp2 = amp2*mass_f**2/vev**2
678  res = cdabs(amp2)
679 
680  RETURN
681  END SUBROUTINE
682 
683 
684 ! Decay amplitude H --> tau^-(-->l^-(p1)+nubar(p2)+nutau(p3)) + tau^+(-->nu(p4)+l^+(p5)+nutaubar(p6))
685 ! or H --> tbar^-(-->l^-(p1)+nubar(p2)+bbar(p3))+top(-->nu(p4)+l^+(p5)+b(p6))
686 ! Note: 1. Value of m allows to change between H->tau^+tau^- and H->ttbar
687 ! 2. Full kinematics of decaying particles -- no NWA.
688 ! 3. Final state b quark is massless.
689 ! 4. At present, no width in the tau/top propagator
690 ! R. Rontsch July 2015
691  SUBROUTINE evalamp_h_tt_decay(pin,mass_F,ga_F,res)
692  implicit none
693  real(dp), intent(out) :: res
694  real(dp), intent(in) :: pin(1:4,1:6),mass_F,ga_F

◆ evalamp_h_vv()

subroutine, public modhiggs::evalamp_h_vv ( real(dp), dimension(4,6), intent(in)  p,
integer, dimension(6:9), intent(in)  MY_IDUP,
real(dp), intent(out)  res 
)

Definition at line 341 of file mod_Higgs.F90.

341  yyy2 = -2*ahz1 !ahz2 ! gauge invariance fixes ahz2 in this case
342  yyy3 = ahz3
343  elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.gszmode) .OR. (vvmode.eq.zgsmode) ) then! decay (Z-photon) OR (gamma*-Z) OR (Z-gamma*)
344  yyy1 = ahz1
345  yyy2 = -2*ahz1*q_q/(q_q-q3_q3)
346  yyy3 = ahz3
347  else
348  yyy1=czero
349  yyy2=czero
350  yyy3=czero
351  print *,"VVMode",vvmode,"not implemented in generate_as"
352  endif
353  endif
354 
355  res= e1_e2*e3_e4*q_q**2*yyy1*xxx1 &
356  + e1_e2*e3_q4*e4_q3*q_q*yyy2*xxx1 &
357  + et1(e1,e2,q1,q2)*e3_e4*q_q*yyy1*xxx3 &
358  + et1(e1,e2,q1,q2)*e3_q4*e4_q3*yyy2*xxx3 &
359  + et1(e1,e2,q1,q2)*et1(e3,e4,q3,q4)*yyy3*xxx3 &
360  + et1(e3,e4,q3,q4)*e1_e2*q_q*yyy3*xxx1
361 
362  END SUBROUTINE gghzzampl
363 
364 
365 !----- a subroutinefor H -> (Z+gamma*)(Z+gamma*)/WW/gammagamma
366 !----- all outgoing convention and the following momentum assignment
367 !----- 0 -> Higgs(p1) + e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
368  subroutine evalamp_h_vv(p,MY_IDUP,res)
369  implicit none
370  real(dp), intent(out) :: res
371  real(dp), intent(in) :: p(4,6)
372  integer, intent(in) :: MY_IDUP(6:9)
373  complex(dp) :: A_VV(1:18),A0_VV(1:2)
374  integer :: i3,i4,VVMode,VVmode_swap
375  real(dp) :: prefactor
376  real(dp) :: intcolfac
377  integer :: ordering(1:4),ordering_swap(1:4)
378  logical :: doInterference
379 
380  if(isaquark(my_idup(6)) .and. isaquark(my_idup(8))) then
381  intcolfac=1.0_dp/3.0_dp
382  else
383  intcolfac=1.0_dp
384  endif
385 
386  call getdecay_vvmode_ordering(my_idup(6:9),vvmode,ordering,vvmode_swap,ordering_swap)
387 
388 ! Global normalization
389  if( vvmode.eq.zzmode ) then! Z decay
390  prefactor = overallcouplvffsq**2
391  elseif( vvmode.eq.wwmode ) then ! W decay
392  prefactor = overallcouplvffsq**2
393  elseif( vvmode.eq.zgmode ) then ! Z+photon "decay"
394  prefactor = overallcouplvffsq ! Only single powers
395  elseif( vvmode.eq.ggmode ) then ! photon "decay"
396  prefactor = 1d0
397  else
398  prefactor = 0d0
399  endif
400 
401 ! ! MADGRAPH CHECK
402 ! res=0d0
403 ! if (MY_IDUP(6).ne.MY_IDUP(8) ) return
404 ! if (MY_IDUP(7).ne.MY_IDUP(9) ) return
405 ! if (MY_IDUP(6).eq.MY_IDUP(8) ) return
406 ! if (MY_IDUP(7).eq.MY_IDUP(9) ) return
407 ! print *, "MY COUPL",al1*dsqrt(overallCouplVffsq),ar1*dsqrt(overallCouplVffsq)
408 ! ar1=-0.158480099490745d0! this is MadGraphs gzl(R) , the MG couplings differ by a global minus sign, relative differences are because of different input parameters
409 ! al1=+0.210150647402957d0! this is MadGraphs gzl(L)
410 ! al2=al1
411 ! ar2=ar1
412 ! print *, "MG COUPL",al1,ar1
413 ! pause
414 
415 
416  res = zero
417  a_vv(:) = 0d0
418  dointerference = includeinterference .and. ( &
419  ((vvmode.eq.zzmode) .and. (vvmode_swap.eq.zzmode)) &
420  )
421  if ( includevprime .and. .not.(vvmode.eq.zzmode .or. vvmode.eq.zgmode .or. vvmode.eq.wwmode) ) then
422  call error("Contact terms only for WW, ZZ or Zg!")
423  endif
424  do i3=1,2; do i4=1,2! sum over helicities
425  call calchelamp2(ordering,vvmode,my_idup,p(1:4,1:6),i3,i4,a_vv(1))
426  if( vvmode.eq.zzmode ) then
427  if( includegammastar ) then
428  call calchelamp2(ordering,zgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(3))
429  call calchelamp2(ordering,gszmode,my_idup,p(1:4,1:6),i3,i4,a_vv(5))
430  call calchelamp2(ordering,gsgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(7))
431  endif
432  if( includevprime ) then
433  call calchelamp2(ordering,zzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(9))
434  call calchelamp2(ordering,zpzmode,my_idup,p(1:4,1:6),i3,i4,a_vv(11))
435  call calchelamp2(ordering,zpzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(13))
436  endif
437  if( includegammastar .and. includevprime ) then
438  call calchelamp2(ordering,gszpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(15))
439  call calchelamp2(ordering,zpgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(17))
440  endif
441  elseif( vvmode.eq.zgmode ) then
442  if(includegammastar) then
443  call calchelamp2(ordering,gsgmode,my_idup,p(1:4,1:6),i3,i4,a_vv(3))
444  endif
445  if( includevprime ) then
446  call calchelamp2(ordering,zpgmode,my_idup,p(1:4,1:6),i3,i4,a_vv(5))
447  endif
448  elseif( vvmode.eq.wwmode .and. includevprime ) then
449  call calchelamp2(ordering,wwpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(9))
450  call calchelamp2(ordering,wpwmode,my_idup,p(1:4,1:6),i3,i4,a_vv(11))
451  call calchelamp2(ordering,wpwpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(13))
452  endif
453 
454  if( dointerference ) then
455  call calchelamp2(ordering_swap,vvmode,my_idup,p(1:4,1:6),i3,i4,a_vv(2))
456  if( includegammastar ) then
457  call calchelamp2(ordering_swap,zgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(4))
458  call calchelamp2(ordering_swap,gszmode,my_idup,p(1:4,1:6),i3,i4,a_vv(6))
459  call calchelamp2(ordering_swap,gsgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(8))
460  endif
461  if( includevprime ) then
462  call calchelamp2(ordering_swap,zzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(10))
463  call calchelamp2(ordering_swap,zpzmode,my_idup,p(1:4,1:6),i3,i4,a_vv(12))
464  call calchelamp2(ordering_swap,zpzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(14))
465  endif
466  if( includegammastar .and. includevprime ) then
467  call calchelamp2(ordering_swap,gszpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(16))
468  call calchelamp2(ordering_swap,zpgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(18))
469  endif
470  endif
471 
472  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
473  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
474  res = res + dreal(a0_vv(1)*dconjg(a0_vv(1)))
475  res = res + dreal(a0_vv(2)*dconjg(a0_vv(2)))
476  if( dointerference .and. (i3.eq.i4) ) then! interfere the 3456 with 5436 pieces
477  res = res - 2d0*intcolfac*dreal( a0_vv(1)*dconjg(a0_vv(2)) ) ! minus from Fermi statistics
478  endif

◆ getdecay_couplings_spinors_props()

subroutine modhiggs::getdecay_couplings_spinors_props ( integer, intent(in)  VVMode,
integer, dimension(6:9), intent(in)  idordered,
real(dp), dimension(1:4,6:9), intent(in)  pordered,
integer, intent(in)  h3 
)
private

Definition at line 701 of file mod_Higgs.F90.

701  call convert_to_mcfm(pin(1:4,j),p(j,1:4))
702  enddo
703 ! call spinoru(6,p,za,zb,s)
704  call spinoru(p,za,zb,s)
705 
706  s12=s(1,2)
707  s45=s(4,5)
708  s123=s(1,2)+s(1,3)+s(2,3)
709  s456=s(4,5)+s(4,6)+s(5,6)
710 
711  kl=-mass_f/vev*( kappa -(0d0,1d0)*kappa_tilde )
712  kr=-mass_f/vev*( kappa +(0d0,1d0)*kappa_tilde )
713 
714  amp = + kr * ( za(1,3)*za(1,4)*zb(1,2)*zb(5,6)- za(1,3)*za(3,4)*zb(2,3)*zb(5,6)) &
715  + kl * (- za(1,3)*za(4,5)*zb(2,5)*zb(5,6)- za(1,3)*za(4,6)*zb(2,6)*zb(5,6))
716 
717 ! overall factors and propagators
718  amp=amp/(s123-mass_f**2+ci*mass_f*ga_f)/(s456-mass_f**2+ci*mass_f*ga_f)/(s12-m_w**2+ci*m_w*ga_w)/(s45-m_w**2+ci*m_w*ga_w)
719  amp=amp*16d0*ci*mass_f*gwsq**2
720  res = cdabs(amp)**2
721 
722  RETURN
723  END SUBROUTINE
724 
725 
726 
727 
728 subroutine getdecay_couplings_spinors_props(VVMode,idordered,pordered,h3,h4, sp,pV)
729  implicit none
730  integer, intent(in) :: VVMode,idordered(6:9),h3,h4
731  real(dp), intent(in) :: pordered(1:4,6:9)
732  complex(dp), intent(out) :: sp(3:4,1:4)
733  real(dp), intent(out) :: pV(3:4,1:4)
734  real(dp) :: s
735  complex(dp) :: propV(1:2), aL1,aR1,aL2,aR2
736 
737  ! h3/h4 helicities: -1 == left, 1 == right
738  if( vvmode.eq.zzmode ) then
739  ! ZZ DECAYS
740  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
741  al1=al_lep * dsqrt(scale_alpha_z_ll)
742  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
743  elseif( abs(idordered(6)).eq.abs(tam_) ) then
744  al1=al_lep * dsqrt(scale_alpha_z_tt)
745  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
746  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
747  al1=al_neu * dsqrt(scale_alpha_z_nn)
748  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
749  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
750  al1=al_qup * dsqrt(scale_alpha_z_uu)
751  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
752  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
753  al1=al_qdn * dsqrt(scale_alpha_z_dd)
754  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
755  else
756  al1=0d0
757  ar1=0d0
758  endif
759  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
760  al2=al_lep * dsqrt(scale_alpha_z_ll)
761  ar2=ar_lep * dsqrt(scale_alpha_z_ll)
762  elseif( abs(idordered(8)).eq.abs(tam_) ) then
763  al2=al_lep * dsqrt(scale_alpha_z_tt)
764  ar2=ar_lep * dsqrt(scale_alpha_z_tt)
765  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
766  al2=al_neu * dsqrt(scale_alpha_z_nn)
767  ar2=ar_neu * dsqrt(scale_alpha_z_nn)
768  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
769  al2=al_qup * dsqrt(scale_alpha_z_uu)
770  ar2=ar_qup * dsqrt(scale_alpha_z_uu)
771  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
772  al2=al_qdn * dsqrt(scale_alpha_z_dd)
773  ar2=ar_qdn * dsqrt(scale_alpha_z_dd)
774  else
775  al2=0d0
776  ar2=0d0
777  endif
778  pv(3,:) = pordered(:,6)+pordered(:,7)
779  pv(4,:) = pordered(:,8)+pordered(:,9)
780  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
781  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
782  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
783  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
784  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
785  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
786  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
787  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
788 
789 
790 
791  elseif( vvmode.eq.wwmode ) then
792  ! WW DECAYS
793  if( isaquark(idordered(6)) ) then
794  al1 = bl * dsqrt(scale_alpha_w_ud)
795  ar1 = br * dsqrt(scale_alpha_w_ud)! = 0
796  elseif( &
797  (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. &
798  (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_)) &
799  ) then
800  al1 = bl * dsqrt(scale_alpha_w_ln)
801  ar1 = br * dsqrt(scale_alpha_w_ln)! = 0
802  elseif( &
803  (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_)) &
804  ) then
805  al1 = bl * dsqrt(scale_alpha_w_tn)
806  ar1 = br * dsqrt(scale_alpha_w_tn)! = 0
807  else
808  al1=0d0
809  ar1=0d0
810  endif
811  if( isaquark(idordered(8)) ) then
812  al2 = bl * dsqrt(scale_alpha_w_ud)
813  ar2 = br * dsqrt(scale_alpha_w_ud)! = 0
814  elseif( &
815  (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. &
816  (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_)) &
817  ) then
818  al2 = bl * dsqrt(scale_alpha_w_ln)
819  ar2 = br * dsqrt(scale_alpha_w_ln)! = 0
820  elseif( &
821  (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_)) &
822  ) then
823  al2 = bl * dsqrt(scale_alpha_w_tn)
824  ar2 = br * dsqrt(scale_alpha_w_tn)! = 0
825  else
826  al2=0d0
827  ar2=0d0
828  endif
829  pv(3,:) = pordered(:,6)+pordered(:,7)
830  pv(4,:) = pordered(:,8)+pordered(:,9)
831  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
832  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
833  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
834  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
835  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
836  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
837  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
838  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
839 
840 
841  elseif( vvmode.eq.zgmode ) then
842  ! Zgamma DECAYS
843  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
844  al1=al_lep * dsqrt(scale_alpha_z_ll)
845  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
846  elseif( abs(idordered(6)).eq.abs(tam_) ) then
847  al1=al_lep * dsqrt(scale_alpha_z_tt)
848  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
849  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
850  al1=al_neu * dsqrt(scale_alpha_z_nn)
851  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
852  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
853  al1=al_qup * dsqrt(scale_alpha_z_uu)
854  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
855  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
856  al1=al_qdn * dsqrt(scale_alpha_z_dd)
857  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
858  else
859  al1=0d0
860  ar1=0d0
861  endif
862  al2=1d0
863  ar2=1d0
864  pv(3,:) = pordered(:,6)+pordered(:,7)
865  pv(4,:) = pordered(:,8)
866  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
867  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
868  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
869  !sp(4,1:4)=pV(4,1:4); print *, "this checks gauge invariance"
870  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
871  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
872  propv(2)=1d0
873 
874 
875  elseif( vvmode.eq.ggmode ) then
876  ! gamma gamma DECAYS
877  al1=1d0
878  ar1=1d0
879  al2=1d0
880  ar2=1d0
881  pv(3,:) = pordered(:,6)
882  pv(4,:) = pordered(:,8)
883  sp(3,:) = pol_mless2(dcmplx(pordered(:,6)),h3,'out') ! photon
884  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
885  !sp(3,1:4)=pV(3,1:4); print *, "this checks gauge invariance"
886  !sp(4,1:4)=pV(4,1:4)
887  propv(1)=1d0
888  propv(2)=1d0
889 
890 
891  elseif( vvmode.eq.gsgmode ) then
892  ! gamma* gamma DECAYS
893  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
894  al1=cl_lep * dsqrt(scale_alpha_z_ll)
895  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
896  elseif( abs(idordered(6)).eq.abs(tam_) ) then
897  al1=cl_lep * dsqrt(scale_alpha_z_tt)
898  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
899  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
900  al1=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
901  ar1=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
902  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
903  al1=cl_qup * dsqrt(scale_alpha_z_uu)
904  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
905  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
906  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
907  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
908  else
909  al1=0d0
910  ar1=0d0
911  endif
912  al2=1d0
913  ar2=1d0
914  pv(3,:) = pordered(:,6)+pordered(:,7)
915  pv(4,:) = pordered(:,8)
916  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
917  sp(3,:) = -sp(3,:) ! photon propagator
918  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
919  !sp(4,1:4)=pV(4,1:4); print *, "this checks gauge invariance"
920  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
921  propv(1) = 1d0
922  propv(2) = 1d0
923  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
924 
925 
926  elseif( vvmode.eq.gszmode ) then
927  ! gamma* Z DECAYS
928  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
929  al1=cl_lep * dsqrt(scale_alpha_z_ll)
930  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
931  elseif( abs(idordered(6)).eq.abs(tam_) ) then
932  al1=cl_lep * dsqrt(scale_alpha_z_tt)
933  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
934  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
935  al1=cl_neu * dsqrt(scale_alpha_z_nn)
936  ar1=cr_neu * dsqrt(scale_alpha_z_nn)
937  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
938  al1=cl_qup * dsqrt(scale_alpha_z_uu)
939  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
940  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
941  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
942  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
943  else
944  al1=0d0
945  ar1=0d0
946  endif
947  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
948  al2=al_lep * dsqrt(scale_alpha_z_ll)
949  ar2=ar_lep * dsqrt(scale_alpha_z_ll)
950  elseif( abs(idordered(8)).eq.abs(tam_) ) then
951  al2=al_lep * dsqrt(scale_alpha_z_tt)
952  ar2=ar_lep * dsqrt(scale_alpha_z_tt)
953  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
954  al2=al_neu * dsqrt(scale_alpha_z_nn)
955  ar2=ar_neu * dsqrt(scale_alpha_z_nn)
956  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
957  al2=al_qup * dsqrt(scale_alpha_z_uu)
958  ar2=ar_qup * dsqrt(scale_alpha_z_uu)
959  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
960  al2=al_qdn * dsqrt(scale_alpha_z_dd)
961  ar2=ar_qdn * dsqrt(scale_alpha_z_dd)
962  else
963  al2=0d0
964  ar2=0d0
965  endif
966  pv(3,:) = pordered(:,6)+pordered(:,7)
967  pv(4,:) = pordered(:,8)+pordered(:,9)
968  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
969  sp(3,:) = -sp(3,:)
970  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
971  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
972  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
973  propv(1) = 1d0! = s/dcmplx(s)
974  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
975  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
976  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
977 
978 
979  elseif( vvmode.eq.zgsmode ) then
980  ! Z gamma* DECAYS
981  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
982  al1=al_lep * dsqrt(scale_alpha_z_ll)
983  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
984  elseif( abs(idordered(6)).eq.abs(tam_) ) then
985  al1=al_lep * dsqrt(scale_alpha_z_tt)
986  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
987  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
988  al1=al_neu * dsqrt(scale_alpha_z_nn)
989  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
990  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
991  al1=al_qup * dsqrt(scale_alpha_z_uu)
992  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
993  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
994  al1=al_qdn * dsqrt(scale_alpha_z_dd)
995  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
996  else
997  al1=0d0
998  ar1=0d0
999  endif
1000  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
1001  al2=cl_lep * dsqrt(scale_alpha_z_ll)
1002  ar2=cr_lep * dsqrt(scale_alpha_z_ll)
1003  elseif( abs(idordered(8)).eq.abs(tam_) ) then
1004  al2=cl_lep * dsqrt(scale_alpha_z_tt)
1005  ar2=cr_lep * dsqrt(scale_alpha_z_tt)
1006  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
1007  al2=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
1008  ar2=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
1009  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
1010  al2=cl_qup * dsqrt(scale_alpha_z_uu)
1011  ar2=cr_qup * dsqrt(scale_alpha_z_uu)
1012  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
1013  al2=cl_qdn * dsqrt(scale_alpha_z_dd)
1014  ar2=cr_qdn * dsqrt(scale_alpha_z_dd)
1015  else
1016  al2=0d0
1017  ar2=0d0
1018  endif
1019  pv(3,:) = pordered(:,6)+pordered(:,7)
1020  pv(4,:) = pordered(:,8)+pordered(:,9)
1021  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1022  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1023  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1024  sp(4,:) = -sp(4,:)
1025  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1026  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
1027  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1028  propv(2) = 1d0 ! = s/dcmplx(s)
1029  if( s.lt.mphotoncutoff**2 ) propv(2)=czero
1030 
1031 
1032  elseif( vvmode.eq.gsgsmode ) then
1033  ! gamma* gamma* DECAYS
1034  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1035  al1=cl_lep * dsqrt(scale_alpha_z_ll)
1036  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
1037  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1038  al1=cl_lep * dsqrt(scale_alpha_z_tt)
1039  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
1040  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1041  al1=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
1042  ar1=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
1043  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1044  al1=cl_qup * dsqrt(scale_alpha_z_uu)
1045  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
1046  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1047  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
1048  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
1049  else
1050  al1=0d0
1051  ar1=0d0
1052  endif
1053  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
1054  al2=cl_lep * dsqrt(scale_alpha_z_ll)
1055  ar2=cr_lep * dsqrt(scale_alpha_z_ll)
1056  elseif( abs(idordered(8)).eq.abs(tam_) ) then
1057  al2=cl_lep * dsqrt(scale_alpha_z_tt)
1058  ar2=cr_lep * dsqrt(scale_alpha_z_tt)
1059  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
1060  al2=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
1061  ar2=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
1062  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
1063  al2=cl_qup * dsqrt(scale_alpha_z_uu)
1064  ar2=cr_qup * dsqrt(scale_alpha_z_uu)
1065  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
1066  al2=cl_qdn * dsqrt(scale_alpha_z_dd)
1067  ar2=cr_qdn * dsqrt(scale_alpha_z_dd)
1068  else
1069  al2=0d0
1070  ar2=0d0
1071  endif
1072  pv(3,:) = pordered(:,6)+pordered(:,7)
1073  pv(4,:) = pordered(:,8)+pordered(:,9)
1074  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1075  sp(3,:) = -sp(3,:)
1076  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1077  sp(4,:) = -sp(4,:)
1078  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1079  propv(1) = 1d0 ! = s/dcmplx(s)
1080  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
1081  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1082  propv(2) = 1d0 ! = s/dcmplx(s)
1083  if( s.lt.mphotoncutoff**2 ) propv(2)=czero
1084 
1085 
1086  elseif( vvmode.eq.zpzmode ) then
1087  ! Z'Z DECAYS
1088  al1 = vpffcoupling(idordered(6),-1,.false.)
1089  ar1 = vpffcoupling(idordered(6),+1,.false.)
1090  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
1091  al2=al_lep * dsqrt(scale_alpha_z_ll)
1092  ar2=ar_lep * dsqrt(scale_alpha_z_ll)
1093  elseif( abs(idordered(8)).eq.abs(tam_) ) then
1094  al2=al_lep * dsqrt(scale_alpha_z_tt)
1095  ar2=ar_lep * dsqrt(scale_alpha_z_tt)
1096  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
1097  al2=al_neu * dsqrt(scale_alpha_z_nn)
1098  ar2=ar_neu * dsqrt(scale_alpha_z_nn)
1099  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
1100  al2=al_qup * dsqrt(scale_alpha_z_uu)
1101  ar2=ar_qup * dsqrt(scale_alpha_z_uu)
1102  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
1103  al2=al_qdn * dsqrt(scale_alpha_z_dd)
1104  ar2=ar_qdn * dsqrt(scale_alpha_z_dd)
1105  else
1106  al2=0d0
1107  ar2=0d0
1108  endif
1109 
1110  pv(3,:) = pordered(:,6)+pordered(:,7)
1111  pv(4,:) = pordered(:,8)+pordered(:,9)
1112  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1113  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1114  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1115  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1116  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1117  if( m_vprime .gt. 0d0 ) then
1118  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1119  elseif( m_vprime .eq. 0d0 ) then
1120  propv(1) = 1d0
1121  else
1122  propv(1) = s/m_z**2
1123  endif
1124  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1125  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
1126 
1127 
1128  elseif( vvmode.eq.zzpmode ) then
1129  ! ZZ' DECAYS
1130  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1131  al1=al_lep * dsqrt(scale_alpha_z_ll)
1132  ar1=ar_lep * dsqrt(scale_alpha_z_ll)
1133  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1134  al1=al_lep * dsqrt(scale_alpha_z_tt)
1135  ar1=ar_lep * dsqrt(scale_alpha_z_tt)
1136  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1137  al1=al_neu * dsqrt(scale_alpha_z_nn)
1138  ar1=ar_neu * dsqrt(scale_alpha_z_nn)
1139  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1140  al1=al_qup * dsqrt(scale_alpha_z_uu)
1141  ar1=ar_qup * dsqrt(scale_alpha_z_uu)
1142  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1143  al1=al_qdn * dsqrt(scale_alpha_z_dd)
1144  ar1=ar_qdn * dsqrt(scale_alpha_z_dd)
1145  else
1146  al1=0d0
1147  ar1=0d0
1148  endif
1149  al2 = vpffcoupling(idordered(8),-1,.false.)
1150  ar2 = vpffcoupling(idordered(8),+1,.false.)
1151 
1152  pv(3,:) = pordered(:,6)+pordered(:,7)
1153  pv(4,:) = pordered(:,8)+pordered(:,9)
1154  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1155  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1156  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1157  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1158  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1159  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
1160  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1161  if( m_vprime .gt. 0d0 ) then
1162  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1163  elseif( m_vprime .eq. 0d0 ) then
1164  propv(2) = 1d0
1165  else
1166  propv(2) = s/m_z**2
1167  endif
1168 
1169 
1170  elseif( vvmode.eq.zpzpmode ) then
1171  ! Z'Z' DECAYS
1172  al1 = vpffcoupling(idordered(6),-1,.false.)
1173  ar1 = vpffcoupling(idordered(6),+1,.false.)
1174  al2 = vpffcoupling(idordered(8),-1,.false.)
1175  ar2 = vpffcoupling(idordered(8),+1,.false.)
1176 
1177  pv(3,:) = pordered(:,6)+pordered(:,7)
1178  pv(4,:) = pordered(:,8)+pordered(:,9)
1179  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1180  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1181  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1182  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1183  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1184  if( m_vprime .gt. 0d0 ) then
1185  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1186  elseif( m_vprime .eq. 0d0 ) then
1187  propv(1) = 1d0
1188  else
1189  propv(1) = s/m_z**2
1190  endif
1191  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1192  if( m_vprime .gt. 0d0 ) then
1193  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1194  elseif( m_vprime .eq. 0d0 ) then
1195  propv(2) = 1d0
1196  else
1197  propv(2) = s/m_z**2
1198  endif
1199 
1200 
1201  elseif( vvmode.eq.zpgsmode ) then
1202  ! Z' gamma* DECAYS
1203  al1 = vpffcoupling(idordered(6),-1,.false.)
1204  ar1 = vpffcoupling(idordered(6),+1,.false.)
1205  if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
1206  al2=cl_lep * dsqrt(scale_alpha_z_ll)
1207  ar2=cr_lep * dsqrt(scale_alpha_z_ll)
1208  elseif( abs(idordered(8)).eq.abs(tam_) ) then
1209  al2=cl_lep * dsqrt(scale_alpha_z_tt)
1210  ar2=cr_lep * dsqrt(scale_alpha_z_tt)
1211  elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
1212  al2=cl_neu * dsqrt(scale_alpha_z_nn)! = 0
1213  ar2=cr_neu * dsqrt(scale_alpha_z_nn)! = 0
1214  elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
1215  al2=cl_qup * dsqrt(scale_alpha_z_uu)
1216  ar2=cr_qup * dsqrt(scale_alpha_z_uu)
1217  elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
1218  al2=cl_qdn * dsqrt(scale_alpha_z_dd)
1219  ar2=cr_qdn * dsqrt(scale_alpha_z_dd)
1220  else
1221  al2=0d0
1222  ar2=0d0
1223  endif
1224 
1225  pv(3,:) = pordered(:,6)+pordered(:,7)
1226  pv(4,:) = pordered(:,8)+pordered(:,9)
1227  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1228  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1229  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1230  sp(4,:) = -sp(4,:)
1231  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1232  if( m_vprime .gt. 0d0 ) then
1233  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1234  elseif( m_vprime .eq. 0d0 ) then
1235  propv(1) = 1d0
1236  else
1237  propv(1) = s/m_z**2
1238  endif
1239  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1240  propv(2) = 1d0 ! = s/dcmplx(s)
1241  if( s.lt.mphotoncutoff**2 ) propv(2)=czero
1242 
1243  elseif( vvmode.eq.gszpmode ) then
1244  ! gamma* Z' DECAYS
1245  if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
1246  al1=cl_lep * dsqrt(scale_alpha_z_ll)
1247  ar1=cr_lep * dsqrt(scale_alpha_z_ll)
1248  elseif( abs(idordered(6)).eq.abs(tam_) ) then
1249  al1=cl_lep * dsqrt(scale_alpha_z_tt)
1250  ar1=cr_lep * dsqrt(scale_alpha_z_tt)
1251  elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
1252  al1=cl_neu * dsqrt(scale_alpha_z_nn)
1253  ar1=cr_neu * dsqrt(scale_alpha_z_nn)
1254  elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
1255  al1=cl_qup * dsqrt(scale_alpha_z_uu)
1256  ar1=cr_qup * dsqrt(scale_alpha_z_uu)
1257  elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
1258  al1=cl_qdn * dsqrt(scale_alpha_z_dd)
1259  ar1=cr_qdn * dsqrt(scale_alpha_z_dd)
1260  else
1261  al1=0d0
1262  ar1=0d0
1263  endif
1264  al2 = vpffcoupling(idordered(8),-1,.false.)
1265  ar2 = vpffcoupling(idordered(8),+1,.false.)
1266 
1267  pv(3,:) = pordered(:,6)+pordered(:,7)
1268  pv(4,:) = pordered(:,8)+pordered(:,9)
1269  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1270  sp(3,:) = -sp(3,:)
1271  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1272  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1273  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1274  propv(1) = 1d0! = s/dcmplx(s)
1275  if( s.lt.mphotoncutoff**2 ) propv(1)=czero
1276  s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
1277  if( m_vprime .gt. 0d0 ) then
1278  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1279  elseif( m_vprime .eq. 0d0 ) then
1280  propv(2) = 1d0
1281  else
1282  propv(2) = s/m_z**2
1283  endif
1284 
1285 
1286  elseif( vvmode.eq.zpgmode ) then
1287  ! Z' gamma DECAYS
1288  al1 = vpffcoupling(idordered(6),-1,.false.)
1289  ar1 = vpffcoupling(idordered(6),+1,.false.)
1290  al2=1d0
1291  ar2=1d0
1292  pv(3,:) = pordered(:,6)+pordered(:,7)
1293  pv(4,:) = pordered(:,8)
1294  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1295  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1296  sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out') ! photon
1297  s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
1298  if( m_vprime .gt. 0d0 ) then
1299  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1300  elseif( m_vprime .eq. 0d0 ) then
1301  propv(1) = 1d0
1302  else
1303  propv(1) = s/m_z**2
1304  endif
1305  propv(2)=1d0
1306 
1307 
1308  elseif( vvmode.eq.wpwmode ) then
1309  ! W'W DECAYS
1310  al1 = vpffcoupling(idordered(6),-1,.true.)
1311  ar1 = vpffcoupling(idordered(6),+1,.true.)
1312  if( isaquark(idordered(8)) ) then
1313  al2 = bl * dsqrt(scale_alpha_w_ud)
1314  ar2 = br * dsqrt(scale_alpha_w_ud)! = 0
1315  elseif( &
1316  (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. &
1317  (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_)) &
1318  ) then
1319  al2 = bl * dsqrt(scale_alpha_w_ln)
1320  ar2 = br * dsqrt(scale_alpha_w_ln)! = 0
1321  elseif( &
1322  (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_)) &
1323  ) then
1324  al2 = bl * dsqrt(scale_alpha_w_tn)
1325  ar2 = br * dsqrt(scale_alpha_w_tn)! = 0
1326  else
1327  al2=0d0
1328  ar2=0d0
1329  endif
1330  pv(3,:) = pordered(:,6)+pordered(:,7)
1331  pv(4,:) = pordered(:,8)+pordered(:,9)
1332  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1333  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1334  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1335  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1336  s = scr(pv(3,:),pv(3,:))
1337  if( m_vprime .gt. 0d0 ) then
1338  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1339  elseif( m_vprime .eq. 0d0 ) then
1340  propv(1) = 1d0
1341  else
1342  propv(1) = s/m_w**2
1343  endif
1344  s = scr(pv(4,:),pv(4,:))
1345  propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
1346 
1347 
1348  elseif( vvmode.eq.wwpmode ) then
1349  ! WW' DECAYS
1350  if( isaquark(idordered(6)) ) then
1351  al1 = bl * dsqrt(scale_alpha_w_ud)
1352  ar1 = br * dsqrt(scale_alpha_w_ud)! = 0
1353  elseif( &
1354  (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. &
1355  (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_)) &
1356  ) then
1357  al1 = bl * dsqrt(scale_alpha_w_ln)
1358  ar1 = br * dsqrt(scale_alpha_w_ln)! = 0
1359  elseif( &
1360  (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_)) &
1361  ) then
1362  al1 = bl * dsqrt(scale_alpha_w_tn)
1363  ar1 = br * dsqrt(scale_alpha_w_tn)! = 0
1364  else
1365  al1=0d0
1366  ar1=0d0
1367  endif
1368  al2 = vpffcoupling(idordered(8),-1,.true.)
1369  ar2 = vpffcoupling(idordered(8),+1,.true.)
1370 
1371  pv(3,:) = pordered(:,6)+pordered(:,7)
1372  pv(4,:) = pordered(:,8)+pordered(:,9)
1373  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1374  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1375  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1376  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1377  s = scr(pv(3,:),pv(3,:))
1378  propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
1379  s = scr(pv(4,:),pv(4,:))
1380  if( m_vprime .gt. 0d0 ) then
1381  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1382  elseif( m_vprime .eq. 0d0 ) then
1383  propv(2) = 1d0
1384  else
1385  propv(2) = s/m_w**2
1386  endif
1387 
1388 
1389  elseif( vvmode.eq.wpwpmode ) then
1390  ! W'W' DECAYS
1391  al1 = vpffcoupling(idordered(6),-1,.true.)
1392  ar1 = vpffcoupling(idordered(6),+1,.true.)
1393  al2 = vpffcoupling(idordered(8),-1,.true.)
1394  ar2 = vpffcoupling(idordered(8),+1,.true.)
1395 
1396  pv(3,:) = pordered(:,6)+pordered(:,7)
1397  pv(4,:) = pordered(:,8)+pordered(:,9)
1398  sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3) ! ubar(l1), v(l2)
1399  sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
1400  sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4) ! ubar(l3), v(l4)
1401  sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
1402  s = scr(pv(3,:),pv(3,:))
1403  if( m_vprime .gt. 0d0 ) then
1404  propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1405  elseif( m_vprime .eq. 0d0 ) then
1406  propv(1) = 1d0
1407  else
1408  propv(1) = s/m_w**2
1409  endif
1410  s = scr(pv(4,:),pv(4,:))
1411  if( m_vprime .gt. 0d0 ) then
1412  propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
1413  elseif( m_vprime .eq. 0d0 ) then
1414  propv(2) = 1d0
1415  else
1416  propv(2) = s/m_w**2

◆ getdecay_vvmode_ordering()

subroutine modhiggs::getdecay_vvmode_ordering ( integer, dimension(6:9), intent(in)  MY_IDUP,
integer, intent(out)  VVMode,
integer, dimension(1:4), intent(out)  ordering,
integer, intent(out)  VVMode_swap 
)
private

Definition at line 1420 of file mod_Higgs.F90.

1420  else
1421  call error("Unsupported decay modes")
1422  endif
1423 
1424  !print *,"sp(3)=",sp(3,:)
1425  !print *,"propV(1)=",propV(1)
1426  !print *,"aL1,aR1=",aL1,aR1
1427  !print *,"sp(4)=",sp(4,:)
1428  !print *,"propV(2)=",propV(2)
1429  !print *,"aL2,aR2=",aL2,aR2
1430 
1431  sp(3,:) = sp(3,:)*propv(1)
1432  sp(4,:) = sp(4,:)*propv(2)
1433  if (h3.eq.-1) then
1434  sp(3,:) = al1 * sp(3,:)
1435  elseif(h3.eq.1) then
1436  sp(3,:) = ar1 * sp(3,:)
1437  endif
1438  if (h4.eq.-1) then
1439  sp(4,:) = al2 * sp(4,:)
1440  elseif(h4.eq.1) then
1441  sp(4,:) = ar2 * sp(4,:)
1442  endif
1443 
1444  return
1445 end subroutine
1446 
1447 subroutine getdecay_vvmode_ordering(MY_IDUP, VVMode,ordering,VVMode_swap,ordering_swap)
1448  implicit none
1449  integer, intent(in) :: MY_IDUP(6:9)
1450  integer, intent(out) :: VVMode,ordering(1:4),VVMode_swap,ordering_swap(1:4)
1451  integer :: idV(1:2),idV_swap(1:2)
1452 
1453  ordering=(/3,4,5,6/)
1454  idv(1)=coupledvertex(my_idup(6:7),-1)
1455  idv(2)=coupledvertex(my_idup(8:9),-1)
1456  if(my_idup(6).eq.pho_ .or. my_idup(7).eq.pho_) idv(1)=pho_
1457  if(my_idup(8).eq.pho_ .or. my_idup(9).eq.pho_) idv(2)=pho_
1458  if(convertlhe(my_idup(6)).lt.0 .or. my_idup(6).eq.not_a_particle_) then
1459  call swap(ordering(1),ordering(2))
1460  endif
1461  if(convertlhe(my_idup(8)).lt.0 .or. my_idup(8).eq.not_a_particle_) then
1462  call swap(ordering(3),ordering(4))
1463  endif
1464  if( &
1465  (idv(1).eq.wm_ .and. idv(2).eq.wp_) .or. &
1466  (idv(2).eq.z0_ .and. idv(1).eq.pho_) &
1467  ) then
1468  call swap(ordering(1),ordering(3))
1469  call swap(ordering(2),ordering(4))
1470  call swap(idv(1),idv(2))
1471  endif
1472  ordering_swap(:)=ordering(:)
1473  call swap(ordering_swap(1),ordering_swap(3))
1474 

◆ gghzzampl()

subroutine modhiggs::gghzzampl ( integer, intent(in)  VVMode,
real(dp), dimension(4,4), intent(in)  p,
complex(dp), dimension(4,4), intent(in)  sp,
complex(dp), intent(out)  res 
)
private

Definition at line 185 of file mod_Higgs.F90.

185  !print *,"id(7)=",convertLHE(MY_IDUP(l2))
186  !print *,"id(8)=",convertLHE(MY_IDUP(l3))
187  !print *,"id(9)=",convertLHE(MY_IDUP(l4))
188  !print *,"p(6)=",(p(:,l1))
189  !print *,"p(7)=",(p(:,l2))
190  !print *,"p(8)=",(p(:,l3))
191  !print *,"p(9)=",(p(:,l4))
192 
193  propg = one/dcmplx(2d0*scr(p(:,1),p(:,2)) - m_reso**2,m_reso*ga_reso)
194 
195  pin(1,:) = p(:,1)
196  pin(2,:) = p(:,2)
197  sp(1,:) = pol_mless2(dcmplx(p(:,1)),-3+2*i1,'in') ! gluon
198  sp(2,:) = pol_mless2(dcmplx(p(:,2)),-3+2*i2,'in') ! gluon
199  !sp(1,1:4)=pin(1,1:4);print *, "this checks IS gauge invariance"
200  !sp(2,1:4)=pin(2,1:4);print *, "this checks IS gauge invariance"
201  call getdecay_couplings_spinors_props( &
202  vvmode, &
203  (/my_idup(l1+3),my_idup(l2+3),my_idup(l3+3),my_idup(l4+3)/), &
204  (/p(:,l1),p(:,l2),p(:,l3),p(:,l4)/), &
205  -3+2*i3,-3+2*i4, &
206  sp(3:4,:),pin(3:4,:) &
207  )
208  call gghzzampl(vvmode,pin,sp,a(1))
209  a(1) = a(1) * propg
210  end subroutine
211 
212  SUBROUTINE gghzzampl(VVMode,p,sp,res)
213  implicit none
214  integer, intent(in) :: VVMode
215  real(dp), intent(in) :: p(4,4)
216  complex(dp), intent(in) :: sp(4,4)
217  complex(dp), intent(out) :: res
218  complex(dp) :: e1_e2, e1_e3, e1_e4
219  complex(dp) :: e2_e3, e2_e4
220  complex(dp) :: e3_e4
221  complex(dp) :: q1_e3,q1_e4,q2_e3,q2_e4
222  complex(dp) :: e1_q3,e1_q4,e2_q3,e2_q4
223  complex(dp) :: e3_q4,e4_q3
224  complex(dp) :: q1(4),q2(4),q3(4),q4(4),q(4)
225  complex(dp) :: e1(4),e2(4),e3(4),e4(4)
226  complex(dp) :: xxx1,xxx2,xxx3,yyy1,yyy2,yyy3,yyy4
227  complex(dp) :: ghz1_dyn,ghz2_dyn,ghz3_dyn,ghz4_dyn
228  real(dp) :: q34
229  real(dp) :: q_q, q3_q3, q4_q4
230 
231 
232  res = 0d0
233  q1 = dcmplx(p(1,:),0d0)
234  q2 = dcmplx(p(2,:),0d0)
235  q3 = dcmplx(p(3,:),0d0)
236  q4 = dcmplx(p(4,:),0d0)
237 
238  e1 = sp(1,:)
239  e2 = sp(2,:)
240  e3 = sp(3,:)
241  e4 = sp(4,:)
242 
243  q = -q1-q2
244 
245  q_q =sc(q,q)
246  q3_q3 = sc(q3,q3)
247  q4_q4 = sc(q4,q4)
248 
249  e1_e2 = sc(e1,e2)
250  e1_e3 = sc(e1,e3)
251  e1_e4 = sc(e1,e4)
252  e2_e3 = sc(e2,e3)
253  e2_e4 = sc(e2,e4)
254  e3_e4 = sc(e3,e4)
255 
256  q1_e3 = sc(q1,e3)
257  q1_e4 = sc(q1,e4)
258  q2_e3 = sc(q2,e3)
259  q2_e4 = sc(q2,e4)
260  e1_q3 = sc(e1,q3)
261  e1_q4 = sc(e1,q4)
262  e2_q3 = sc(e2,q3)
263  e2_q4 = sc(e2,q4)
264  e3_q4 = sc(e3,q4)
265  e4_q3 = sc(e4,q3)
266 
267 
268  if (q_q.lt.-0.1d0 .or. q3_q3.lt.-0.1d0 .or. q4_q4.lt.-0.1d0) return ! if negative invariant masses return zero
269 
270 !---- data that defines couplings
271  if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode) ) then! decay ZZ's or WW's
272  ghz1_dyn = hvvspinzerodynamiccoupling(1,q3_q3,q4_q4,q_q)
273  ghz2_dyn = hvvspinzerodynamiccoupling(2,q3_q3,q4_q4,q_q)
274  ghz3_dyn = hvvspinzerodynamiccoupling(3,q3_q3,q4_q4,q_q)
275  ghz4_dyn = hvvspinzerodynamiccoupling(4,q3_q3,q4_q4,q_q)
276  elseif( (vvmode.eq.gszmode) ) then
277  ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q3_q3,q_q)
278  ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q3_q3,q_q)
279  ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q3_q3,q_q)
280  ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q3_q3,q_q)
281  elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.zgsmode) ) then
282  ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q4_q4,q_q)
283  ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q4_q4,q_q)
284  ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q4_q4,q_q)
285  ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q4_q4,q_q)
286  elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode) .or. (vvmode.eq.gsgmode) ) then
287  ghz1_dyn = czero
288  ghz2_dyn = hvvspinzerodynamiccoupling(9,q3_q3,q4_q4,q_q)
289  ghz3_dyn = hvvspinzerodynamiccoupling(10,q3_q3,q4_q4,q_q)
290  ghz4_dyn = hvvspinzerodynamiccoupling(11,q3_q3,q4_q4,q_q)
291  elseif( (vvmode.eq.zzpmode) .or. (vvmode.eq.wwpmode) ) then
292  ghz1_dyn = hvvspinzerodynamiccoupling(12,q3_q3,q4_q4,q_q)
293  ghz2_dyn = hvvspinzerodynamiccoupling(13,q3_q3,q4_q4,q_q)
294  ghz3_dyn = hvvspinzerodynamiccoupling(14,q3_q3,q4_q4,q_q)
295  ghz4_dyn = hvvspinzerodynamiccoupling(15,q3_q3,q4_q4,q_q)
296  elseif( (vvmode.eq.zpzmode) .or. (vvmode.eq.wpwmode) ) then
297  ghz1_dyn = hvvspinzerodynamiccoupling(12,q4_q4,q3_q3,q_q)
298  ghz2_dyn = hvvspinzerodynamiccoupling(13,q4_q4,q3_q3,q_q)
299  ghz3_dyn = hvvspinzerodynamiccoupling(14,q4_q4,q3_q3,q_q)
300  ghz4_dyn = hvvspinzerodynamiccoupling(15,q4_q4,q3_q3,q_q)
301  elseif( (vvmode.eq.zpzpmode) .or. (vvmode.eq.wpwpmode) ) then
302  ghz1_dyn = hvvspinzerodynamiccoupling(16,q3_q3,q4_q4,q_q)
303  ghz2_dyn = hvvspinzerodynamiccoupling(17,q3_q3,q4_q4,q_q)
304  ghz3_dyn = hvvspinzerodynamiccoupling(18,q3_q3,q4_q4,q_q)
305  ghz4_dyn = hvvspinzerodynamiccoupling(19,q3_q3,q4_q4,q_q)
306  elseif( (vvmode.eq.gszpmode) ) then
307  ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q3_q3,q_q)
308  ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q3_q3,q_q)
309  ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q3_q3,q_q)
310  ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q3_q3,q_q)
311  elseif( (vvmode.eq.zpgmode) .OR. (vvmode.eq.zpgsmode) ) then
312  ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q4_q4,q_q)
313  ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q4_q4,q_q)
314  ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q4_q4,q_q)
315  ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q4_q4,q_q)
316  else
317  ghz1_dyn = czero
318  ghz2_dyn = czero
319  ghz3_dyn = czero
320  ghz4_dyn = czero
321  print *,"VVMode",vvmode,"not implemented"
322  endif
323 
324  if( .not. generate_as ) then
325  xxx1 = ghg2+ghg3/4d0/lambda**2*q_q
326  xxx3 = -2d0*ghg4
327  yyy1 = ghz1_dyn*m_v**2/q_q & ! in this line M_V is indeed correct, not a misprint
328  + ghz2_dyn*(q_q-q3_q3-q4_q4)/q_q &
329  + ghz3_dyn/lambda**2*(q_q-q3_q3-q4_q4)*(q_q-q4_q4-q3_q3)/4d0/q_q
330  yyy2 = -2d0*ghz2_dyn-ghz3_dyn/2d0/lambda**2*(q_q-q3_q3-q4_q4)
331  yyy3 = -2d0*ghz4_dyn
332  else
333  xxx1 = ahg1

◆ hzzampl()

subroutine modhiggs::hzzampl ( integer, intent(in)  VVMode,
real(dp), dimension(4,4), intent(in)  p,
complex(dp), dimension(3:4,4), intent(in)  sp,
complex(dp), intent(out)  res 
)
private

Definition at line 506 of file mod_Higgs.F90.

506  RETURN
507  END SUBROUTINE
508 
509  subroutine calchelamp2(ordering,VVMode,MY_IDUP,p,i3,i4,A)
510  implicit none
511  integer :: ordering(1:4),VVMode,i3,i4,l1,l2,l3,l4,MY_IDUP(6:9)
512  real(dp) :: p(1:4,1:6)
513  real(dp) :: pin(4,4)
514  complex(dp) :: A(1:1), sp(3:4,4)
515 
516  l1=ordering(1)
517  l2=ordering(2)
518  l3=ordering(3)
519  l4=ordering(4)
520 
521  pin(1,:) = p(:,1)
522 
523  call getdecay_couplings_spinors_props( &
524  vvmode, &
525  (/my_idup(l1+3),my_idup(l2+3),my_idup(l3+3),my_idup(l4+3)/), &
526  (/p(:,l1),p(:,l2),p(:,l3),p(:,l4)/), &
527  -3+2*i3,-3+2*i4, &
528  sp(3:4,:),pin(3:4,:) &
529  )
530  call hzzampl(vvmode,pin,sp,a(1))
531  end subroutine
532 
533  subroutine hzzampl(VVMode,p,sp,res)
534  implicit none
535  integer, intent(in) :: VVMode
536  real(dp), intent(in) :: p(4,4)
537  complex(dp), intent(in) :: sp(3:4,4)
538  complex(dp), intent(out) :: res
539  complex(dp) :: e3_e4
540  complex(dp) :: e3_q4,e4_q3
541  complex(dp) :: q1(4),q3(4),q4(4),q(4)
542  complex(dp) :: e3(4),e4(4)
543  complex(dp) :: yyy1,yyy2,yyy3,yyy4
544  real(dp) :: q34
545  real(dp) :: q_q, q3_q3, q4_q4
546  complex(dp) :: ghz1_dyn,ghz2_dyn,ghz3_dyn,ghz4_dyn
547 
548 
549 
550  res = 0d0
551  q1 = dcmplx(p(1,:),0d0)
552  q3 = dcmplx(p(3,:),0d0)
553  q4 = dcmplx(p(4,:),0d0)
554 
555 
556  e3 = sp(3,:)
557  e4 = sp(4,:)
558 
559  q = -q1
560 
561  q_q =sc(q,q)
562  q3_q3 = sc(q3,q3)
563  q4_q4 = sc(q4,q4)
564 
565  e3_e4 = sc(e3,e4)
566  e3_q4 = sc(e3,q4)
567  e4_q3 = sc(e4,q3)
568 
569 
570  if ((q_q).lt.-0.1d0 .or. (q3_q3).lt.-0.1d0 .or. (q4_q4).lt.-0.1d0) return ! if negative invariant masses return zero
571 
572 
573 !---- data that defines couplings
574  if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode) ) then! decay ZZ's or WW's
575  ghz1_dyn = hvvspinzerodynamiccoupling(1,q3_q3,q4_q4,q_q)
576  ghz2_dyn = hvvspinzerodynamiccoupling(2,q3_q3,q4_q4,q_q)
577  ghz3_dyn = hvvspinzerodynamiccoupling(3,q3_q3,q4_q4,q_q)
578  ghz4_dyn = hvvspinzerodynamiccoupling(4,q3_q3,q4_q4,q_q)
579  elseif( (vvmode.eq.gszmode) ) then
580  ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q3_q3,q_q)
581  ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q3_q3,q_q)
582  ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q3_q3,q_q)
583  ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q3_q3,q_q)
584  elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.zgsmode) ) then
585  ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q4_q4,q_q)
586  ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q4_q4,q_q)
587  ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q4_q4,q_q)
588  ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q4_q4,q_q)
589  elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode) .or. (vvmode.eq.gsgmode) ) then
590  ghz1_dyn = czero
591  ghz2_dyn = hvvspinzerodynamiccoupling(9,q3_q3,q4_q4,q_q)
592  ghz3_dyn = hvvspinzerodynamiccoupling(10,q3_q3,q4_q4,q_q)
593  ghz4_dyn = hvvspinzerodynamiccoupling(11,q3_q3,q4_q4,q_q)
594  elseif( (vvmode.eq.zzpmode) .or. (vvmode.eq.wwpmode) ) then
595  ghz1_dyn = hvvspinzerodynamiccoupling(12,q3_q3,q4_q4,q_q)
596  ghz2_dyn = hvvspinzerodynamiccoupling(13,q3_q3,q4_q4,q_q)
597  ghz3_dyn = hvvspinzerodynamiccoupling(14,q3_q3,q4_q4,q_q)
598  ghz4_dyn = hvvspinzerodynamiccoupling(15,q3_q3,q4_q4,q_q)
599  elseif( (vvmode.eq.zpzmode) .or. (vvmode.eq.wpwmode) ) then
600  ghz1_dyn = hvvspinzerodynamiccoupling(12,q4_q4,q3_q3,q_q)
601  ghz2_dyn = hvvspinzerodynamiccoupling(13,q4_q4,q3_q3,q_q)
602  ghz3_dyn = hvvspinzerodynamiccoupling(14,q4_q4,q3_q3,q_q)
603  ghz4_dyn = hvvspinzerodynamiccoupling(15,q4_q4,q3_q3,q_q)
604  elseif( (vvmode.eq.zpzpmode) .or. (vvmode.eq.wpwpmode) ) then
605  ghz1_dyn = hvvspinzerodynamiccoupling(16,q3_q3,q4_q4,q_q)
606  ghz2_dyn = hvvspinzerodynamiccoupling(17,q3_q3,q4_q4,q_q)
607  ghz3_dyn = hvvspinzerodynamiccoupling(18,q3_q3,q4_q4,q_q)
608  ghz4_dyn = hvvspinzerodynamiccoupling(19,q3_q3,q4_q4,q_q)
609  elseif( (vvmode.eq.gszpmode) ) then
610  ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q3_q3,q_q)
611  ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q3_q3,q_q)
612  ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q3_q3,q_q)
613  ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q3_q3,q_q)
614  elseif( (vvmode.eq.zpgmode) .OR. (vvmode.eq.zpgsmode) ) then
615  ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q4_q4,q_q)
616  ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q4_q4,q_q)
617  ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q4_q4,q_q)
618  ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q4_q4,q_q)
619  else
620  ghz1_dyn = czero
621  ghz2_dyn = czero
622  ghz3_dyn = czero
623  ghz4_dyn = czero
624  print *,"VVMode",vvmode,"not implemented"
625  endif
626 
627  if( .not. generate_as ) then
628  yyy1 = ghz1_dyn*m_v**2/q_q & ! in this line M_V is indeed correct, not a misprint
629  + ghz2_dyn*(q_q-q3_q3-q4_q4)/q_q &
630  + ghz3_dyn/lambda**2*(q_q-q3_q3-q4_q4)*(q_q-q4_q4-q3_q3)/4d0/q_q
hto_units::one
real *8, parameter one
Definition: CALLING_cpHTO.f:184
hto_units::ci
real *8, dimension(1:2) ci
Definition: CALLING_cpHTO.f:189
hto_units::zero
real *8, parameter zero
Definition: CALLING_cpHTO.f:185