mod_JHUGenMELA.F90

Go to the documentation of this file.

MODULE modjhugenmela
   use modparameters
   use modkinematics
   use modmisc
#if compiler==1
   use ifport
#endif
   implicit none
 
! JHUGenMELA-specific subroutines and functions
 
public :: setewparameters
public :: sethiggsmasswidth
public :: setdecaymodes
public :: settopdecays
public :: sethdk
public :: setmurenfac
public :: resetmubarhgabarh
public :: resetamplitudeincludes
public :: setdistinguishwwcouplingsflag
public :: setspinzeroggcouplings,setspinzeroqqcouplings,setspinzerovvcouplings
public :: setspinonecouplings
public :: setspintwocouplings
public :: setvprimecontactcouplings
public :: setazffcouplings
 
public :: getmvgv,getmvprimegvprime
public :: getalphasalphasmz
public :: getpdfconstants
public :: getdecaycouplings
 
 
contains
 
 
!=====================================================
! Subroutines visible to the user
!=====================================================
 
subroutine setewparameters(inMZ, inMW, inGf, inalpha_QED, inxw)
   implicit none
   real(8), intent(in) :: inmz, inmw, ingf, inalpha_qed, inxw
   m_z = inmz
   m_w = inmw
   gf = ingf
   alpha_qed = inalpha_qed
   xw = inxw
   call computeewvariables()
end subroutine setewparameters
 
 
subroutine sethiggsmasswidth(mass,width)
   implicit none
   real(8), intent(in) :: mass, width
   call setmass(mass,hig_)
   call setdecaywidth(width,hig_)
   return
end subroutine sethiggsmasswidth
 
subroutine setdecaymodes(idfirst,idsecond)
   implicit none
   integer, intent(in) :: idfirst(1:2)
   integer, intent(in) :: idsecond(1:2)
   integer :: idv(1:2)
 
   if( idfirst(1).eq.pho_ .or. idfirst(2).eq.pho_ ) then
      decaymode1 = 7
      idv(1)=pho_
   else
      idv(1)=coupledvertex(idfirst,-1)
      if( idv(1).eq.wp_ .or. idv(1).eq.wm_ ) then
         decaymode1=11
      elseif( idv(1).eq.z0_ ) then
         decaymode1=9
      endif
   endif
   if( idsecond(1).eq.pho_ .or. idsecond(2).eq.pho_ ) then
      decaymode2 = 7
      idv(2)=pho_
   else
      idv(2)=coupledvertex(idsecond,-1)
      if( idv(2).eq.wp_ .or. idv(2).eq.wm_ ) then
         decaymode2=11
      elseif( idv(2).eq.z0_ ) then
         decaymode2=9
      endif
   endif
   call setmvgv()
   return
end subroutine setdecaymodes
 
subroutine settopdecays(flag)
implicit none
integer, intent(in) :: flag
   topdecays=flag
end subroutine
 
subroutine sethdk(flag)
implicit none
logical, intent(in) :: flag
   h_dk=flag
end subroutine
 
subroutine setmurenfac(muren,mufac)
implicit none
real(8), intent(in) :: muren,mufac
    mu_ren = muren
   mu_fact = mufac
end subroutine
 
subroutine resetmubarhgabarh()
implicit none
   mubarh=-999d0
   gabarh=-999d0
end subroutine
 
subroutine resetamplitudeincludes()
implicit none
   includevprime=.false.
   includegammastar=.false.
end subroutine
 
subroutine setspinzerovvcouplings(vvcoupl, vvpcoupl, vpvpcoupl, cqsq, Lambda_qsq, useWWcoupl)
   implicit none
   complex(8), intent(in) :: vvcoupl(39)
   complex(8), intent(in) :: vvpcoupl(39)
   complex(8), intent(in) :: vpvpcoupl(39)
   integer, intent(in) :: cqsq(3)
   real(8), intent(in) :: lambda_qsq(1:3,1:4)
   logical, intent(in) :: usewwcoupl
 
   includevprime = includevprime .or. (                                   &
                    (any(vvpcoupl.ne.czero) .or. any(vpvpcoupl.ne.czero)) &
                   )
 
   if(.not.usewwcoupl) then
      cz_q1sq = cqsq(1)
      lambda_z11 = lambda_qsq(1,1)
      lambda_z21 = lambda_qsq(1,2)
      lambda_z31 = lambda_qsq(1,3)
      lambda_z41 = lambda_qsq(1,4)
      cz_q2sq = cqsq(2)
      lambda_z12 = lambda_qsq(2,1)
      lambda_z22 = lambda_qsq(2,2)
      lambda_z32 = lambda_qsq(2,3)
      lambda_z42 = lambda_qsq(2,4)
      cz_q12sq = cqsq(3)
      lambda_z10 = lambda_qsq(3,1)
      lambda_z20 = lambda_qsq(3,2)
      lambda_z30 = lambda_qsq(3,3)
      lambda_z40 = lambda_qsq(3,4)
 
      ghz1 =  vvcoupl(1)
      ghz2 =  vvcoupl(2)
      ghz3 =  vvcoupl(3)
      ghz4 =  vvcoupl(4)
 
      ghzgs2  = vvcoupl(5)
      ghzgs3  = vvcoupl(6)
      ghzgs4  = vvcoupl(7)
      ghgsgs2 = vvcoupl(8)
      ghgsgs3 = vvcoupl(9)
      ghgsgs4 = vvcoupl(10)
 
      ghz1_prime = vvcoupl(11)
      ghz1_prime2= vvcoupl(12)
      ghz1_prime3= vvcoupl(13)
      ghz1_prime4= vvcoupl(14)
      ghz1_prime5= vvcoupl(15)
 
      ghz2_prime = vvcoupl(16)
      ghz2_prime2= vvcoupl(17)
      ghz2_prime3= vvcoupl(18)
      ghz2_prime4= vvcoupl(19)
      ghz2_prime5= vvcoupl(20)
 
      ghz3_prime = vvcoupl(21)
      ghz3_prime2= vvcoupl(22)
      ghz3_prime3= vvcoupl(23)
      ghz3_prime4= vvcoupl(24)
      ghz3_prime5= vvcoupl(25)
 
      ghz4_prime = vvcoupl(26)
      ghz4_prime2= vvcoupl(27)
      ghz4_prime3= vvcoupl(28)
      ghz4_prime4= vvcoupl(29)
      ghz4_prime5= vvcoupl(30)
 
      ghzgs1_prime2= vvcoupl(31)
 
      ghz1_prime6  = vvcoupl(32)
      ghz1_prime7  = vvcoupl(33)
      ghz2_prime6  = vvcoupl(34)
      ghz2_prime7  = vvcoupl(35)
      ghz3_prime6  = vvcoupl(36)
      ghz3_prime7  = vvcoupl(37)
      ghz4_prime6  = vvcoupl(38)
      ghz4_prime7  = vvcoupl(39)
 
 
      ghzzp1 =  vvpcoupl(1)
      ghzzp2 =  vvpcoupl(2)
      ghzzp3 =  vvpcoupl(3)
      ghzzp4 =  vvpcoupl(4)
 
      ghzpgs2  = vvpcoupl(5)
      ghzpgs3  = vvpcoupl(6)
      ghzpgs4  = vvpcoupl(7)
 
      ghzzp1_prime = vvpcoupl(11)
      ghzzp1_prime2= vvpcoupl(12)
      ghzzp1_prime3= vvpcoupl(13)
      ghzzp1_prime4= vvpcoupl(14)
      ghzzp1_prime5= vvpcoupl(15)
 
      ghzzp2_prime = vvpcoupl(16)
      ghzzp2_prime2= vvpcoupl(17)
      ghzzp2_prime3= vvpcoupl(18)
      ghzzp2_prime4= vvpcoupl(19)
      ghzzp2_prime5= vvpcoupl(20)
 
      ghzzp3_prime = vvpcoupl(21)
      ghzzp3_prime2= vvpcoupl(22)
      ghzzp3_prime3= vvpcoupl(23)
      ghzzp3_prime4= vvpcoupl(24)
      ghzzp3_prime5= vvpcoupl(25)
 
      ghzzp4_prime = vvpcoupl(26)
      ghzzp4_prime2= vvpcoupl(27)
      ghzzp4_prime3= vvpcoupl(28)
      ghzzp4_prime4= vvpcoupl(29)
      ghzzp4_prime5= vvpcoupl(30)
 
      ghzpgs1_prime2= vvpcoupl(31)
 
      ghzzp1_prime6  = vvpcoupl(32)
      ghzzp1_prime7  = vvpcoupl(33)
      ghzzp2_prime6  = vvpcoupl(34)
      ghzzp2_prime7  = vvpcoupl(35)
      ghzzp3_prime6  = vvpcoupl(36)
      ghzzp3_prime7  = vvpcoupl(37)
      ghzzp4_prime6  = vvpcoupl(38)
      ghzzp4_prime7  = vvpcoupl(39)
 
      ghzpzp1 =  vpvpcoupl(1)
      ghzpzp2 =  vpvpcoupl(2)
      ghzpzp3 =  vpvpcoupl(3)
      ghzpzp4 =  vpvpcoupl(4)
 
      ghzpzp1_prime = vpvpcoupl(11)
      ghzpzp1_prime2= vpvpcoupl(12)
      ghzpzp1_prime3= vpvpcoupl(13)
      ghzpzp1_prime4= vpvpcoupl(14)
      ghzpzp1_prime5= vpvpcoupl(15)
 
      ghzpzp2_prime = vpvpcoupl(16)
      ghzpzp2_prime2= vpvpcoupl(17)
      ghzpzp2_prime3= vpvpcoupl(18)
      ghzpzp2_prime4= vpvpcoupl(19)
      ghzpzp2_prime5= vpvpcoupl(20)
 
      ghzpzp3_prime = vpvpcoupl(21)
      ghzpzp3_prime2= vpvpcoupl(22)
      ghzpzp3_prime3= vpvpcoupl(23)
      ghzpzp3_prime4= vpvpcoupl(24)
      ghzpzp3_prime5= vpvpcoupl(25)
 
      ghzpzp4_prime = vpvpcoupl(26)
      ghzpzp4_prime2= vpvpcoupl(27)
      ghzpzp4_prime3= vpvpcoupl(28)
      ghzpzp4_prime4= vpvpcoupl(29)
      ghzpzp4_prime5= vpvpcoupl(30)
 
      ghzpzp1_prime6  = vpvpcoupl(32)
      ghzpzp1_prime7  = vpvpcoupl(33)
      ghzpzp2_prime6  = vpvpcoupl(34)
      ghzpzp2_prime7  = vpvpcoupl(35)
      ghzpzp3_prime6  = vpvpcoupl(36)
      ghzpzp3_prime7  = vpvpcoupl(37)
      ghzpzp4_prime6  = vpvpcoupl(38)
      ghzpzp4_prime7  = vpvpcoupl(39)
 
 
      ! Set includeGammaStar based on the actual coupling values
      includegammastar = (          &
         ghzgs1_prime2.ne.czero .or. &
         ghzgs2.ne.czero .or.        &
         ghzgs3.ne.czero .or.        &
         ghzgs4.ne.czero .or.        &
         ghgsgs2.ne.czero .or.       &
         ghgsgs3.ne.czero .or.       &
         ghgsgs4.ne.czero .or.       &
         ghzpgs1_prime2.ne.czero .or. &
         ghzpgs2.ne.czero .or.        &
         ghzpgs3.ne.czero .or.        &
         ghzpgs4.ne.czero             &
         )
 
   else
      cw_q1sq = cqsq(1)
      lambda_w11 = lambda_qsq(1,1)
      lambda_w21 = lambda_qsq(1,2)
      lambda_w31 = lambda_qsq(1,3)
      lambda_w41 = lambda_qsq(1,4)
      cw_q2sq = cqsq(2)
      lambda_w12 = lambda_qsq(2,1)
      lambda_w22 = lambda_qsq(2,2)
      lambda_w32 = lambda_qsq(2,3)
      lambda_w42 = lambda_qsq(2,4)
      cw_q12sq = cqsq(3)
      lambda_w10 = lambda_qsq(3,1)
      lambda_w20 = lambda_qsq(3,2)
      lambda_w30 = lambda_qsq(3,3)
      lambda_w40 = lambda_qsq(3,4)
 
      ghw1 =  vvcoupl(1)
      ghw2 =  vvcoupl(2)
      ghw3 =  vvcoupl(3)
      ghw4 =  vvcoupl(4)
 
      ghw1_prime = vvcoupl(11)
      ghw1_prime2= vvcoupl(12)
      ghw1_prime3= vvcoupl(13)
      ghw1_prime4= vvcoupl(14)
      ghw1_prime5= vvcoupl(15)
 
      ghw2_prime = vvcoupl(16)
      ghw2_prime2= vvcoupl(17)
      ghw2_prime3= vvcoupl(18)
      ghw2_prime4= vvcoupl(19)
      ghw2_prime5= vvcoupl(20)
 
      ghw3_prime = vvcoupl(21)
      ghw3_prime2= vvcoupl(22)
      ghw3_prime3= vvcoupl(23)
      ghw3_prime4= vvcoupl(24)
      ghw3_prime5= vvcoupl(25)
 
      ghw4_prime = vvcoupl(26)
      ghw4_prime2= vvcoupl(27)
      ghw4_prime3= vvcoupl(28)
      ghw4_prime4= vvcoupl(29)
      ghw4_prime5= vvcoupl(30)
 
      ghw1_prime6  = vvcoupl(32)
      ghw1_prime7  = vvcoupl(33)
      ghw2_prime6  = vvcoupl(34)
      ghw2_prime7  = vvcoupl(35)
      ghw3_prime6  = vvcoupl(36)
      ghw3_prime7  = vvcoupl(37)
      ghw4_prime6  = vvcoupl(38)
      ghw4_prime7  = vvcoupl(39)
 
      ghwwp1 =  vvpcoupl(1)
      ghwwp2 =  vvpcoupl(2)
      ghwwp3 =  vvpcoupl(3)
      ghwwp4 =  vvpcoupl(4)
 
      ghwwp1_prime = vvpcoupl(11)
      ghwwp1_prime2= vvpcoupl(12)
      ghwwp1_prime3= vvpcoupl(13)
      ghwwp1_prime4= vvpcoupl(14)
      ghwwp1_prime5= vvpcoupl(15)
 
      ghwwp2_prime = vvpcoupl(16)
      ghwwp2_prime2= vvpcoupl(17)
      ghwwp2_prime3= vvpcoupl(18)
      ghwwp2_prime4= vvpcoupl(19)
      ghwwp2_prime5= vvpcoupl(20)
 
      ghwwp3_prime = vvpcoupl(21)
      ghwwp3_prime2= vvpcoupl(22)
      ghwwp3_prime3= vvpcoupl(23)
      ghwwp3_prime4= vvpcoupl(24)
      ghwwp3_prime5= vvpcoupl(25)
 
      ghwwp4_prime = vvpcoupl(26)
      ghwwp4_prime2= vvpcoupl(27)
      ghwwp4_prime3= vvpcoupl(28)
      ghwwp4_prime4= vvpcoupl(29)
      ghwwp4_prime5= vvpcoupl(30)
 
      ghwwp1_prime6  = vvpcoupl(32)
      ghwwp1_prime7  = vvpcoupl(33)
      ghwwp2_prime6  = vvpcoupl(34)
      ghwwp2_prime7  = vvpcoupl(35)
      ghwwp3_prime6  = vvpcoupl(36)
      ghwwp3_prime7  = vvpcoupl(37)
      ghwwp4_prime6  = vvpcoupl(38)
      ghwwp4_prime7  = vvpcoupl(39)
 
      ghwpwp1 =  vpvpcoupl(1)
      ghwpwp2 =  vpvpcoupl(2)
      ghwpwp3 =  vpvpcoupl(3)
      ghwpwp4 =  vpvpcoupl(4)
 
      ghwpwp1_prime = vpvpcoupl(11)
      ghwpwp1_prime2= vpvpcoupl(12)
      ghwpwp1_prime3= vpvpcoupl(13)
      ghwpwp1_prime4= vpvpcoupl(14)
      ghwpwp1_prime5= vpvpcoupl(15)
 
      ghwpwp2_prime = vpvpcoupl(16)
      ghwpwp2_prime2= vpvpcoupl(17)
      ghwpwp2_prime3= vpvpcoupl(18)
      ghwpwp2_prime4= vpvpcoupl(19)
      ghwpwp2_prime5= vpvpcoupl(20)
 
      ghwpwp3_prime = vpvpcoupl(21)
      ghwpwp3_prime2= vpvpcoupl(22)
      ghwpwp3_prime3= vpvpcoupl(23)
      ghwpwp3_prime4= vpvpcoupl(24)
      ghwpwp3_prime5= vpvpcoupl(25)
 
      ghwpwp4_prime = vpvpcoupl(26)
      ghwpwp4_prime2= vpvpcoupl(27)
      ghwpwp4_prime3= vpvpcoupl(28)
      ghwpwp4_prime4= vpvpcoupl(29)
      ghwpwp4_prime5= vpvpcoupl(30)
 
      ghwpwp1_prime6  = vpvpcoupl(32)
      ghwpwp1_prime7  = vpvpcoupl(33)
      ghwpwp2_prime6  = vpvpcoupl(34)
      ghwpwp2_prime7  = vpvpcoupl(35)
      ghwpwp3_prime6  = vpvpcoupl(36)
      ghwpwp3_prime7  = vpvpcoupl(37)
      ghwpwp4_prime6  = vpvpcoupl(38)
      ghwpwp4_prime7  = vpvpcoupl(39)
   endif
   return
end subroutine setspinzerovvcouplings
 
subroutine setvprimecontactcouplings(Zpffcoupl, Wpffcoupl)
   implicit none
   complex(8), intent(in) :: zpffcoupl(20)
   complex(8), intent(in) :: wpffcoupl(20)
 
   ezp_el_left = zpffcoupl(1)
   ezp_el_right = zpffcoupl(2)
   ezp_mu_left = zpffcoupl(3)
   ezp_mu_right = zpffcoupl(4)
   ezp_ta_left = zpffcoupl(5)
   ezp_ta_right = zpffcoupl(6)
   ezp_nue_left = zpffcoupl(7)
   ezp_nue_right = zpffcoupl(8)
   ezp_dn_left = zpffcoupl(9)
   ezp_dn_right = zpffcoupl(10)
   ezp_up_left = zpffcoupl(11)
   ezp_up_right = zpffcoupl(12)
   ezp_str_left = zpffcoupl(13)
   ezp_str_right = zpffcoupl(14)
   ezp_chm_left = zpffcoupl(15)
   ezp_chm_right = zpffcoupl(16)
   ezp_bot_left = zpffcoupl(17)
   ezp_bot_right = zpffcoupl(18)
   ezp_top_left = zpffcoupl(19)
   ezp_top_right = zpffcoupl(20)
 
   ewp_el_left = wpffcoupl(1)
   ewp_el_right = wpffcoupl(2)
   ewp_mu_left = wpffcoupl(3)
   ewp_mu_right = wpffcoupl(4)
   ewp_ta_left = wpffcoupl(5)
   ewp_ta_right = wpffcoupl(6)
   ewp_up_left = wpffcoupl(11)
   ewp_up_right = wpffcoupl(12)
   ewp_chm_left = wpffcoupl(15)
   ewp_chm_right = wpffcoupl(16)
   ewp_top_left = wpffcoupl(19)
   ewp_top_right = wpffcoupl(20)
end subroutine
 
subroutine setdistinguishwwcouplingsflag(doAllow)
   implicit none
   logical, intent(in) :: doallow
   distinguish_hwwcouplings = doallow
   return
end subroutine setdistinguishwwcouplingsflag
 
subroutine setspinzeroggcouplings(ggcoupl)
   implicit none
   complex(8), intent(in) :: ggcoupl(1:3)
   ghg2 =  ggcoupl(1)
   ghg3 =  ggcoupl(2)
   ghg4 =  ggcoupl(3)
   return
end subroutine setspinzeroggcouplings
 
subroutine setspinzeroqqcouplings(qqcoupl)
   implicit none
   complex(8), intent(in) :: qqcoupl(1:2)
   kappa =  qqcoupl(1)
   kappa_tilde =  qqcoupl(2)
   return
end subroutine setspinzeroqqcouplings
 
subroutine setspinonecouplings(qqcoupl,vvcoupl)
   implicit none
   complex(8), intent(in) :: qqcoupl(1:2)
   complex(8), intent(in) :: vvcoupl(1:2)
 
   zprime_qq_left  = qqcoupl(1)
   zprime_qq_right = qqcoupl(2)
   zprime_zz_1 = vvcoupl(1)
   zprime_zz_2 = vvcoupl(2)
   return
end subroutine setspinonecouplings
 
subroutine setspintwocouplings(acoupl,vvcoupl,vvpcoupl,vpvpcoupl,qLR)
   implicit none
   integer, parameter :: indexgammabegin=11
   integer, parameter :: indexvvsize=20
   complex(8), intent(in) :: acoupl(1:5)
   complex(8), intent(in) :: vvcoupl(1:indexvvsize),vvpcoupl(1:indexvvsize),vpvpcoupl(1:indexvvsize)
   complex(8), intent(in) :: qlr(1:2)
 
   includevprime = (any(vvpcoupl.ne.czero) .or. any(vpvpcoupl.ne.czero))
   includegammastar = (                                                       &
                    any(vvcoupl(indexgammabegin:indexvvsize).ne.czero) .or.   &
                    any(vvpcoupl(indexgammabegin:indexvvsize).ne.czero)       &
                   )
 
   a1 = acoupl(1)
   a2 = acoupl(2)
   a3 = acoupl(3)
   a4 = acoupl(4)
   a5 = acoupl(5)
 
   graviton_qq_left  = qlr(1)
   graviton_qq_right = qlr(2)
 
   b1 = vvcoupl(1)
   b2 = vvcoupl(2)
   b3 = vvcoupl(3)
   b4 = vvcoupl(4)
   b5 = vvcoupl(5)
   b6 = vvcoupl(6)
   b7 = vvcoupl(7)
   b8 = vvcoupl(8)
   b9 = vvcoupl(9)
   b10 = vvcoupl(10)
 
   bzgs1 = vvcoupl(11)
   bzgs2 = vvcoupl(12)
   bzgs3 = vvcoupl(13)
   bzgs4 = vvcoupl(14)
   bzgs8 = vvcoupl(15)
 
   bgsgs1 = vvcoupl(16)
   bgsgs2 = vvcoupl(17)
   bgsgs3 = vvcoupl(18)
   bgsgs4 = vvcoupl(19)
   bgsgs8 = vvcoupl(20)
 
   bzzp1 = vvpcoupl(1)
   bzzp2 = vvpcoupl(2)
   bzzp3 = vvpcoupl(3)
   bzzp4 = vvpcoupl(4)
   bzzp5 = vvpcoupl(5)
   bzzp6 = vvpcoupl(6)
   bzzp7 = vvpcoupl(7)
   bzzp8 = vvpcoupl(8)
   bzzp9 = vvpcoupl(9)
   bzzp10 = vvpcoupl(10)
 
   bzpgs1 = vvpcoupl(11)
   bzpgs2 = vvpcoupl(12)
   bzpgs3 = vvpcoupl(13)
   bzpgs4 = vvpcoupl(14)
   bzpgs8 = vvpcoupl(15)
 
 
   bzpzp1 = vpvpcoupl(1)
   bzpzp2 = vpvpcoupl(2)
   bzpzp3 = vpvpcoupl(3)
   bzpzp4 = vpvpcoupl(4)
   bzpzp5 = vpvpcoupl(5)
   bzpzp6 = vpvpcoupl(6)
   bzpzp7 = vpvpcoupl(7)
   bzpzp8 = vpvpcoupl(8)
   bzpzp9 = vpvpcoupl(9)
   bzpzp10 = vpvpcoupl(10)
 
   return
end subroutine setspintwocouplings
 
subroutine setazffcouplings(azffcoupl)
   implicit none
   complex(8), intent(in) :: azffcoupl(1:12)
   ar_lep = azffcoupl(1) ! rez
   al_lep = azffcoupl(2) ! lez
   ar_qup = azffcoupl(3) ! rquz
   al_qup = azffcoupl(4) ! lquz
   ar_qdn = azffcoupl(5) ! rqdz
   al_qdn = azffcoupl(6) ! lqdz
   ar_neu = azffcoupl(7) ! rnz
   al_neu = azffcoupl(8) ! lnz
   return
end subroutine setazffcouplings
 
 
subroutine setmvgv()
implicit none
  if( isazdecay(decaymode1) ) then
    m_v = m_z
    ga_v= ga_z
    m_vprime = m_zprime
    ga_vprime= ga_zprime
  elseif( isawdecay(decaymode1) ) then
    m_v = m_w
    ga_v= ga_w
    m_vprime = m_wprime
    ga_vprime= ga_wprime
  else
    m_v = 0d0
    ga_v= 0d0
    m_vprime = -1d0
    ga_vprime= 0d0
  endif
end subroutine setmvgv
 
subroutine setmvgvfromvertex(idV)
implicit none
integer idV
  if( idv.eq.z0_ ) then
    m_v = m_z
    ga_v= ga_z
    m_vprime = m_zprime
    ga_vprime= ga_zprime
  else if( idv.eq.wp_ .or. idv.eq.wm_ ) then
    m_v = m_w
    ga_v= ga_w
    m_vprime = m_wprime
    ga_vprime= ga_wprime
  else
    m_v = 0d0
    ga_v= 0d0
    m_vprime = -1d0
    ga_vprime= 0d0
  endif
end subroutine
 
subroutine getmvgv(mv,gv)
   implicit none
   real(8), intent(out) :: mv,gv
   mv=m_v
   gv=ga_v
end subroutine getmvgv
 
subroutine getmvprimegvprime(mv,gv)
   implicit none
   real(8), intent(out) :: mv,gv
   mv=m_vprime
   gv=ga_vprime
end subroutine getmvprimegvprime
 
subroutine getalphasalphasmz(val_as, val_asmz)
implicit none
real(8), intent(out) :: val_as, val_asmz
   val_as=alphas
   val_asmz=alphas_mz
end subroutine
 
subroutine getpdfconstants(pdfzmass, pdfnloops, pdfnf)
implicit none
real(8), intent(out) :: pdfzmass
integer, intent(out) :: pdfnloops, pdfnf
   pdfzmass=zmass_pdf
   pdfnloops=nloops_pdf
   pdfnf=nqflavors_pdf
end subroutine getpdfconstants
 
 
! This subroutine is slightly different form the one in the decay MEs in the sense that onshell photon returns 0,0 instead of 1,1
subroutine getdecaycouplings(VVMode,idordered,aL1,aR1,aL2,aR2)
   implicit none
   integer, intent(in) :: vvmode,idordered(6:9)
   real(dp), intent(out) :: al1,ar1,al2,ar2
 
   !        h3/h4 helicities: -1 == left, 1 == right
   if( vvmode.eq.zzmode ) then
   !        ZZ DECAYS
      if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_)  ) then
         al1=al_lep    * dsqrt(scale_alpha_z_ll)
         ar1=ar_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(6)).eq.abs(tam_) ) then
         al1=al_lep    * dsqrt(scale_alpha_z_tt)
         ar1=ar_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
         al1=al_neu    * dsqrt(scale_alpha_z_nn)
         ar1=ar_neu    * dsqrt(scale_alpha_z_nn)
      elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
         al1=al_qup    * dsqrt(scale_alpha_z_uu)
         ar1=ar_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
         al1=al_qdn    * dsqrt(scale_alpha_z_dd)
         ar1=ar_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al1=0d0
         ar1=0d0
      endif
      if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_)  ) then
         al2=al_lep    * dsqrt(scale_alpha_z_ll)
         ar2=ar_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(8)).eq.abs(tam_) ) then
         al2=al_lep    * dsqrt(scale_alpha_z_tt)
         ar2=ar_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
         al2=al_neu    * dsqrt(scale_alpha_z_nn)
         ar2=ar_neu    * dsqrt(scale_alpha_z_nn)
      elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
         al2=al_qup    * dsqrt(scale_alpha_z_uu)
         ar2=ar_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
         al2=al_qdn    * dsqrt(scale_alpha_z_dd)
         ar2=ar_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al2=0d0
         ar2=0d0
      endif
 
   elseif( vvmode.eq.wwmode ) then
   !        WW DECAYS
      if( isaquark(idordered(6)) ) then
         al1 = bl * dsqrt(scale_alpha_w_ud)
         ar1 = br * dsqrt(scale_alpha_w_ud)! = 0
      elseif( &
               (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. &
               (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_))      &
            ) then
         al1 = bl * dsqrt(scale_alpha_w_ln)
         ar1 = br * dsqrt(scale_alpha_w_ln)! = 0
      elseif( &
               (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_))      &
            ) then
         al1 = bl * dsqrt(scale_alpha_w_tn)
         ar1 = br * dsqrt(scale_alpha_w_tn)! = 0
      else
         al1=0d0
         ar1=0d0
      endif
      if( isaquark(idordered(8)) ) then
         al2 = bl * dsqrt(scale_alpha_w_ud)
         ar2 = br * dsqrt(scale_alpha_w_ud)! = 0
      elseif( &
               (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. &
               (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_))      &
            ) then
         al2 = bl * dsqrt(scale_alpha_w_ln)
         ar2 = br * dsqrt(scale_alpha_w_ln)! = 0
      elseif( &
               (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_))      &
            ) then
         al2 = bl * dsqrt(scale_alpha_w_tn)
         ar2 = br * dsqrt(scale_alpha_w_tn)! = 0
      else
         al2=0d0
         ar2=0d0
      endif
 
   elseif( vvmode.eq.zgmode ) then
   !        Zgamma DECAYS
      if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
         al1=al_lep    * dsqrt(scale_alpha_z_ll)
         ar1=ar_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(6)).eq.abs(tam_) ) then
         al1=al_lep    * dsqrt(scale_alpha_z_tt)
         ar1=ar_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
         al1=al_neu    * dsqrt(scale_alpha_z_nn)
         ar1=ar_neu    * dsqrt(scale_alpha_z_nn)
      elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
         al1=al_qup    * dsqrt(scale_alpha_z_uu)
         ar1=ar_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
         al1=al_qdn    * dsqrt(scale_alpha_z_dd)
         ar1=ar_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al1=0d0
         ar1=0d0
      endif
      al2=0d0
      ar2=0d0
 
   elseif( vvmode.eq.ggmode ) then
   !        gamma gamma DECAYS
      al1=0d0
      ar1=0d0
      al2=0d0
      ar2=0d0
 
   elseif( vvmode.eq.gsgmode ) then
   !        gamma* gamma DECAYS
      if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
         al1=cl_lep    * dsqrt(scale_alpha_z_ll)
         ar1=cr_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(6)).eq.abs(tam_) ) then
         al1=cl_lep    * dsqrt(scale_alpha_z_tt)
         ar1=cr_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
         al1=cl_neu    * dsqrt(scale_alpha_z_nn)! = 0
         ar1=cr_neu    * dsqrt(scale_alpha_z_nn)! = 0
      elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
         al1=cl_qup    * dsqrt(scale_alpha_z_uu)
         ar1=cr_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
         al1=cl_qdn    * dsqrt(scale_alpha_z_dd)
         ar1=cr_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al1=0d0
         ar1=0d0
      endif
      al2=0d0
      ar2=0d0
 
   elseif( vvmode.eq.gszmode ) then
   !        gamma* Z DECAYS
      if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_)  ) then
         al1=cl_lep    * dsqrt(scale_alpha_z_ll)
         ar1=cr_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(6)).eq.abs(tam_) ) then
         al1=cl_lep    * dsqrt(scale_alpha_z_tt)
         ar1=cr_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
         al1=cl_neu    * dsqrt(scale_alpha_z_nn)
         ar1=cr_neu    * dsqrt(scale_alpha_z_nn)
      elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
         al1=cl_qup    * dsqrt(scale_alpha_z_uu)
         ar1=cr_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
         al1=cl_qdn    * dsqrt(scale_alpha_z_dd)
         ar1=cr_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al1=0d0
         ar1=0d0
      endif
      if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
         al2=al_lep    * dsqrt(scale_alpha_z_ll)
         ar2=ar_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(8)).eq.abs(tam_) ) then
         al2=al_lep    * dsqrt(scale_alpha_z_tt)
         ar2=ar_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
         al2=al_neu    * dsqrt(scale_alpha_z_nn)
         ar2=ar_neu    * dsqrt(scale_alpha_z_nn)
      elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
         al2=al_qup    * dsqrt(scale_alpha_z_uu)
         ar2=ar_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
         al2=al_qdn    * dsqrt(scale_alpha_z_dd)
         ar2=ar_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al2=0d0
         ar2=0d0
      endif
 
   elseif( vvmode.eq.zgsmode ) then
   !        Z gamma* DECAYS
      if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
         al1=al_lep    * dsqrt(scale_alpha_z_ll)
         ar1=ar_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(6)).eq.abs(tam_) ) then
         al1=al_lep    * dsqrt(scale_alpha_z_tt)
         ar1=ar_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
         al1=al_neu    * dsqrt(scale_alpha_z_nn)
         ar1=ar_neu    * dsqrt(scale_alpha_z_nn)
      elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
         al1=al_qup    * dsqrt(scale_alpha_z_uu)
         ar1=ar_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
         al1=al_qdn    * dsqrt(scale_alpha_z_dd)
         ar1=ar_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al1=0d0
         ar1=0d0
      endif
      if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
         al2=cl_lep    * dsqrt(scale_alpha_z_ll)
         ar2=cr_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(8)).eq.abs(tam_) ) then
         al2=cl_lep    * dsqrt(scale_alpha_z_tt)
         ar2=cr_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
         al2=cl_neu    * dsqrt(scale_alpha_z_nn)! = 0
         ar2=cr_neu    * dsqrt(scale_alpha_z_nn)! = 0
      elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
         al2=cl_qup    * dsqrt(scale_alpha_z_uu)
         ar2=cr_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
         al2=cl_qdn    * dsqrt(scale_alpha_z_dd)
         ar2=cr_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al2=0d0
         ar2=0d0
      endif
 
   elseif( vvmode.eq.gsgsmode ) then
   !        gamma* gamma* DECAYS
      if( abs(idordered(6)).eq.abs(elm_) .or. abs(idordered(6)).eq.abs(mum_) ) then
         al1=cl_lep    * dsqrt(scale_alpha_z_ll)
         ar1=cr_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(6)).eq.abs(tam_) ) then
         al1=cl_lep    * dsqrt(scale_alpha_z_tt)
         ar1=cr_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(6)).eq.abs(nue_) .or. abs(idordered(6)).eq.abs(num_) .or. abs(idordered(6)).eq.abs(nut_) ) then
         al1=cl_neu    * dsqrt(scale_alpha_z_nn)! = 0
         ar1=cr_neu    * dsqrt(scale_alpha_z_nn)! = 0
      elseif( abs(idordered(6)).eq.abs(up_) .or. abs(idordered(6)).eq.abs(chm_) ) then
         al1=cl_qup    * dsqrt(scale_alpha_z_uu)
         ar1=cr_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(6)).eq.abs(dn_) .or. abs(idordered(6)).eq.abs(str_) .or. abs(idordered(6)).eq.abs(bot_) ) then
         al1=cl_qdn    * dsqrt(scale_alpha_z_dd)
         ar1=cr_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al1=0d0
         ar1=0d0
      endif
      if( abs(idordered(8)).eq.abs(elm_) .or. abs(idordered(8)).eq.abs(mum_) ) then
         al2=cl_lep    * dsqrt(scale_alpha_z_ll)
         ar2=cr_lep    * dsqrt(scale_alpha_z_ll)
      elseif( abs(idordered(8)).eq.abs(tam_) ) then
         al2=cl_lep    * dsqrt(scale_alpha_z_tt)
         ar2=cr_lep    * dsqrt(scale_alpha_z_tt)
      elseif( abs(idordered(8)).eq.abs(nue_) .or. abs(idordered(8)).eq.abs(num_) .or. abs(idordered(8)).eq.abs(nut_) ) then
         al2=cl_neu    * dsqrt(scale_alpha_z_nn)! = 0
         ar2=cr_neu    * dsqrt(scale_alpha_z_nn)! = 0
      elseif( abs(idordered(8)).eq.abs(up_) .or. abs(idordered(8)).eq.abs(chm_) ) then
         al2=cl_qup    * dsqrt(scale_alpha_z_uu)
         ar2=cr_qup    * dsqrt(scale_alpha_z_uu)
      elseif( abs(idordered(8)).eq.abs(dn_) .or. abs(idordered(8)).eq.abs(str_) .or. abs(idordered(8)).eq.abs(bot_) ) then
         al2=cl_qdn    * dsqrt(scale_alpha_z_dd)
         ar2=cr_qdn    * dsqrt(scale_alpha_z_dd)
      else
         al2=0d0
         ar2=0d0
      endif
 
   else
      al1=0d0
      ar1=0d0
      al2=0d0
      ar2=0d0
   endif
 
   return
end subroutine
 
 
END MODULE modjhugenmela