mod_Zprime.F90

Go to the documentation of this file.

      module modzprime
      use modparameters
      use modmisc
      implicit none
      private
 
 
!----- notation for subroutines
      public :: evalamp_qqb_zprime_vv,evalamp_zprime_vv
      public :: calchelamp,calchelamp2
 
      contains
 
!----- a subroutine for q qbar -> Zprime -> ZZ/WW
!----- all outgoing convention and the following momentum assignment
!-----  0 -> bq(p1) + q(p2) + e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
     subroutine evalamp_qqb_zprime_vv(p,MY_IDUP,sum)
      implicit none
      real(dp), intent(out) ::  sum
      real(dp), intent(in) :: p(4,6)
      integer, intent(in) :: MY_IDUP(6:9)
      real(dp) :: pin(4,4)
      complex(dp) :: A(2)
      integer :: i1,i2,i3,i4,ordering(1:4),ordering_swap(1:4),idV(1:2),VVmode
      real(dp) :: prefactor
      real(dp) :: intcolfac
 
      if(isaquark(my_idup(6)) .and. isaquark(my_idup(8))) then
         intcolfac=1.0_dp/3.0_dp
      else
         intcolfac=1.0_dp
      endif
 
      call getdecay_vvmode_ordering(my_idup(6:9),vvmode,ordering,ordering_swap)
 
!---- full prefactor; 3 is  the color factor
      if( vvmode.eq.zzmode ) then!  Z decay
         prefactor = 3d0*overallcouplvffsq**2
      elseif( vvmode.eq.wwmode ) then !  W decay
         prefactor = 3d0*overallcouplvffsq**2
      elseif( vvmode.eq.zgmode ) then !  Z+photon "decay"
         prefactor = 3d0*overallcouplvffsq ! Only single powers
      elseif( vvmode.eq.ggmode ) then !  photon "decay"
         prefactor = 3d0
      else
         prefactor = 0d0
      endif
 
 
      sum = zero
do i1=1,2
do i3 = 1,2
do i4 = 1,2
!do i3 = -1,1! on-shell check!
!do i4 = -1,1! on-shell check!
 
         call calchelamp(ordering,vvmode,p(1:4,1:6),my_idup,i1,i3,i4,a(1))
         if( includeinterference .and. (my_idup(6).eq.my_idup(8)) .and. (my_idup(7).eq.my_idup(9)) ) then
             call calchelamp(ordering_swap,vvmode,p(1:4,1:6),my_idup,i1,i3,i4,a(2))
             a(2) = -a(2) ! minus comes from fermi statistics
         endif
         if( includeinterference .and. (my_idup(6).eq.my_idup(8)) .and. (my_idup(7).eq.my_idup(9)) ) then
             sum = sum + symmfac * (cdabs( a(1)*dconjg(a(1)) ) + cdabs( a(2)*dconjg(a(2)) ))
             if( i3.eq.i4 ) sum = sum + symmfac * 2d0*intcolfac*dreal(a(1)*dconjg(a(2)))
         else
             sum = sum + cdabs( a(1)*dconjg(a(1)) )
         endif
 
enddo
enddo
enddo
 
!  sum ~ GeV~-8
!  prefactor ~ GeV^0
         sum = sum*prefactor
 
      end subroutine
 
   subroutine calchelamp(ordering,VVmode,p,MY_IDUP,i1,i3,i4,A)
      implicit none
      integer :: ordering(1:4),i1,i3,i4,l1,l2,l3,l4,vvmode,my_idup(6:9)
      real(dp) :: p(1:4,1:6)
      complex(dp) :: propg
      real(dp) :: s, pin(4,4)
      complex(dp) :: a(1:1), sp(4,4)
 
 
      l1=ordering(1)
      l2=ordering(2)
      l3=ordering(3)
      l4=ordering(4)
 
      !print *,"p(6)=",(p(:,l1))
      !print *,"p(7)=",(p(:,l2))
      !print *,"p(8)=",(p(:,l3))
      !print *,"p(9)=",(p(:,l4))
      !pause
 
      s  = two*scr(p(:,1),p(:,2))
      propg = s/dcmplx(s - m_reso**2,m_reso*ga_reso)
!-------- For fermions, -1 == left, 1 == right
      pin(1,:) = p(:,1)
      pin(2,:) = p(:,2)
      sp(1,:) = pol_dk2mom(dcmplx(p(:,2)),dcmplx(p(:,1)),-3+2*i1) !qbq
      sp(2,:) = sp(1,:)  !-- the same, isn't really needed but for uniform bookeeping
!-- on-shell check
!       sp(3,1:4) = pol_mass(dcmplx(p(1:4,3)+p(1:4,4)),dsqrt(2d0*scr(p(:,3),p(:,4))),i3)
!       sp(4,1:4) = pol_mass(dcmplx(p(1:4,5)+p(1:4,6)),dsqrt(2d0*scr(p(:,5),p(:,6))),i4)
 
!        if( i3.eq.0 .and. i4.eq.0 ) then
!             a(:) = 0d0
!             return
!        endif
!-- end: on-shell check
      call getdecay_couplings_spinors_props(                                                             &
                                       vvmode,                                                      &
                                       (/my_idup(l1+3),my_idup(l2+3),my_idup(l3+3),my_idup(l4+3)/), &
                                       (/p(:,l1),p(:,l2),p(:,l3),p(:,l4)/),                         &
                                       -3+2*i3,-3+2*i4,                                             &
                                       sp(3:4,:),pin(3:4,:)                                         &
                                      )
      call qqzprimezzampl(pin,sp,a(1))
!---- chiral couplings of quarks to Zprimes
      if (i1.eq.1) then
         a(1) = zprime_qq_left*a(1)
      else if(i1.eq.2) then
         a(1) = zprime_qq_right*a(1)
      endif
      a(1) = a(1)*propg
! A(1) ~ GeV^-2
! propG*propZ1*propZ2 ~  GeV^-2
   end subroutine
 
      subroutine qqzprimezzampl(p,sp,res)
      implicit none
      real(dp), intent(in) :: p(4,4)
      complex(dp), intent(in) :: sp(4,4)
      complex(dp), intent(out) :: res
      complex(dp) :: e1_e2, e1_e3, e1_e4
      complex(dp) :: e2_e3, e2_e4
      complex(dp) :: e3_e4
      complex(dp) :: q_q
      complex(dp) :: q1_q2,q1_q3,q1_q4
      complex(dp) :: q2_q3,q2_q4
      complex(dp) :: q3_q4
      complex(dp) :: q1_e3,q1_e4,q2_e3,q2_e4
      complex(dp) :: e1_q3,e1_q4,e2_q3,e2_q4
      complex(dp) :: e3_q4,e4_q3
      complex(dp) :: q1(4),q2(4),q3(4),q4(4),q(4)
      complex(dp) :: e1(4),e2(4),e3(4),e4(4)
      complex(dp) :: yyy1,yyy2,yyy3,yyy4,xxx1,epsZpr(1:4,-1:+1)
 
      q1 = dcmplx(p(1,:),0d0)
      q2 = dcmplx(p(2,:),0d0)
      q3 = dcmplx(p(3,:),0d0)
      q4 = dcmplx(p(4,:),0d0)
 
 
      e1 = sp(1,:)
      e2 = sp(2,:)
      e3 = sp(3,:)
      e4 = sp(4,:)
 
      q = -q1-q2
 
      q_q =sc(q,q)
 
 
      q1_q2 = sc(q1,q2)
      q1_q3 = sc(q1,q3)
      q1_q4 = sc(q1,q4)
      q2_q3 = sc(q2,q3)
      q2_q4 = sc(q2,q4)
      q3_q4 = sc(q3,q4)
      e1_e2 = sc(e1,e2)
 
      e1_e3 = sc(e1,e3)
      e1_e4 = sc(e1,e4)
 
 
 
 
 
 
!       epsZpr(1:4,-1) = pol_mass(q,m_reso,-1)
!       epsZpr(1:4, 0) = pol_mass(q,m_reso, 0)
!       epsZpr(1:4,+1) = pol_mass(q,m_reso,+1)
!       e1_e3 = sc(e1,epsZpr(1:4,-1)) * sc(e3,dconjg(epsZpr(1:4,-1)))  &
!             + sc(e1,epsZpr(1:4,+1)) * sc(e3,dconjg(epsZpr(1:4,+1)))  &
!             + sc(e1,epsZpr(1:4, 0)) * sc(e3,dconjg(epsZpr(1:4, 0)))
!       e1_e4 = sc(e1,epsZpr(1:4,-1)) * sc(e4,dconjg(epsZpr(1:4,-1)))  &
!             + sc(e1,epsZpr(1:4,+1)) * sc(e4,dconjg(epsZpr(1:4,+1)))  &
!             + sc(e1,epsZpr(1:4, 0)) * sc(e4,dconjg(epsZpr(1:4, 0)))
 
 
 
 
      e2_e3 = sc(e2,e3)
      e2_e4 = sc(e2,e4)
 
      e3_e4 = sc(e3,e4)
 
 
      q1_e3 = sc(q1,e3)
      q1_e4 = sc(q1,e4)
      q2_e3 = sc(q2,e3)
      q2_e4 = sc(q2,e4)
      e1_q3 = sc(e1,q3)
      e1_q4 = sc(e1,q4)
      e2_q3 = sc(e2,q3)
      e2_q4 = sc(e2,q4)
      e3_q4 = sc(e3,q4)
      e4_q3 = sc(e4,q3)
 
 
      res = czero
 
      xxx1 = (1d0,0d0)  !  different possibilities for fermion couplings
                        !  to zprime are accounted in the amplitude call
 
 
      yyy1 = zprime_zz_1
      yyy2 = zprime_zz_2
 
 
       res= - e1_e3*e4_q3*xxx1*yyy1    &
            - e1_e4*e3_q4*xxx1*yyy1    &
            - et1(e1,e3,e4,q3)*xxx1*yyy2 + et1(e1,e3,e4,q4)*xxx1*yyy2
 
 
!       e3 = pol_mass(q3,dreal(cdsqrt(sc(q3,q3))),-1)
!       e4 = pol_mass(q4,dreal(cdsqrt(sc(q4,q4))),-1)
!       e1_e3 = sc(e1,e3)
!       e1_e4 = sc(e1,e4)
!       e3_q4 = sc(e3,q4)
!       e4_q3 = sc(e4,q3)
!print *, dreal(cdsqrt(sc(q3,q3))),dreal(cdsqrt(sc(q4,q4)))
!print * , "e3x",q3,dreal(cdsqrt(sc(q3,q3))),-1
!print *, - e1_e3*e4_q3 - e1_e4*e3_q4
!pause
 
 
 
      end subroutine qqzprimezzampl
 
 
 
!----- a subroutine for Zprime -> ZZ/WW
!----- all outgoing convention and the following momentum assignment
!-----  0 --> Zprime(p1) e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
     subroutine evalamp_zprime_vv(p,MY_IDUP,sum)
      implicit none
      real(dp), intent(out) ::  sum
      real(dp), intent(in) :: p(4,6)
      integer, intent(in) :: MY_IDUP(6:9)
      real(dp) :: pin(4,4)
      complex(dp) :: A(2)
      integer :: i1,i2,i3,i4,ordering(1:4),ordering_swap(1:4),idV(1:2),VVmode
      real(dp) :: aL1,aR1,aL2,aR2
      real(dp) :: prefactor
      real(dp) :: intcolfac
 
      if(isaquark(my_idup(6)) .and. isaquark(my_idup(8))) then
         intcolfac=1.0_dp/3.0_dp
      else
         intcolfac=1.0_dp
      endif
 
      call getdecay_vvmode_ordering(my_idup(6:9),vvmode,ordering,ordering_swap)
 
!---- full prefactor
      if( vvmode.eq.zzmode ) then!  Z decay
         prefactor = overallcouplvffsq**2
      elseif( vvmode.eq.wwmode ) then !  W decay
         prefactor = overallcouplvffsq**2
      elseif( vvmode.eq.zgmode ) then !  Z+photon "decay"
         prefactor = overallcouplvffsq ! Only single powers
      elseif( vvmode.eq.ggmode ) then !  photon "decay"
         prefactor = 1d0
      else
         prefactor = 0d0
      endif
 
 
      sum = zero
 
do i1 =-1,1! Z' boson
do i3 = 1,2! lepton string1
do i4 = 1,2! lepton string2
 
         call calchelamp2(ordering,vvmode,p(1:4,1:6),my_idup,i1,i3,i4,a(1))
         if( includeinterference .and. (my_idup(6).eq.my_idup(8)) .and. (my_idup(7).eq.my_idup(9)) ) then
             call calchelamp2(ordering_swap,vvmode,p(1:4,1:6),my_idup,i1,i3,i4,a(2))
             a(2) = -a(2) ! minus comes from fermi statistics
         endif
         if( includeinterference .and. (my_idup(6).eq.my_idup(8)) .and. (my_idup(7).eq.my_idup(9)) ) then
             sum = sum + symmfac * (cdabs( a(1)*dconjg(a(1)) ) + cdabs( a(2)*dconjg(a(2)) ))
             if( i3.eq.i4 ) sum = sum + symmfac * 2d0*intcolfac*dreal(a(1)*dconjg(a(2)))
         else
             sum = sum + cdabs( a(1)*dconjg(a(1)) )
         endif
 
enddo
enddo
enddo
 
         sum = sum*prefactor
 
      end subroutine
 
   subroutine calchelamp2(ordering,VVMode,p,MY_IDUP,i1,i3,i4,A)
      implicit none
      integer :: ordering(1:4),i1,i3,i4,l1,l2,l3,l4,vvmode,my_idup(6:9)
      real(dp) :: p(1:4,1:6)
      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)
 
      pin(1,:) = p(:,1)
      sp(1,:) = pol_mass2(dcmplx(p(1:4,1)),i1,'in')
      pin(2,:) = 0d0  ! dummy
      sp(2,:)  = czero
      call getdecay_couplings_spinors_props(                                                             &
                                       vvmode,                                                      &
                                       (/my_idup(l1+3),my_idup(l2+3),my_idup(l3+3),my_idup(l4+3)/), &
                                       (/p(:,l1),p(:,l2),p(:,l3),p(:,l4)/),                         &
                                       -3+2*i3,-3+2*i4,                                             &
                                       sp(3:4,:),pin(3:4,:)                                         &
                                      )
      call zprimezzampl(pin,sp,a(1))
   end subroutine
 
      subroutine zprimezzampl(p,sp,res)
      implicit none
      real(dp), intent(in) :: p(4,4)
      complex(dp), intent(in) :: sp(4,4)
      complex(dp), intent(out) :: res
      complex(dp) :: e1_e2, e1_e3, e1_e4
      complex(dp) :: e2_e3, e2_e4
      complex(dp) :: e3_e4
      complex(dp) :: q_q,MZ3,MZ4
      complex(dp) :: q1_q2,q1_q3,q1_q4
      complex(dp) :: q2_q3,q2_q4
      complex(dp) :: q3_q4,q3_q3,q4_q4
      complex(dp) :: q1_e3,q1_e4,q2_e3,q2_e4
      complex(dp) :: e1_q3,e1_q4,e2_q3,e2_q4
      complex(dp) :: e3_q4,e4_q3
      complex(dp) :: q1(4),q2(4),q3(4),q4(4),q(4)
      complex(dp) :: e1(4),e2(4),e3(4),e4(4)
      complex(dp) :: yyy1,yyy2,yyy3,yyy4,epsZpr(1:4,-1:+1)
 
 
      q1 = dcmplx(p(1,:),0d0)
      q2 = dcmplx(p(2,:),0d0)
      q3 = dcmplx(p(3,:),0d0)
      q4 = dcmplx(p(4,:),0d0)
 
 
      e1 = sp(1,:)
      e2 = sp(2,:)
      e3 = sp(3,:)
      e4 = sp(4,:)
 
      q = -q1-q2
 
      q_q =sc(q,q)
 
 
      q1_q2 = sc(q1,q2)
      q3_q3 = sc(q3,q3)
      q4_q4 = sc(q4,q4)
      q1_q3 = sc(q1,q3)
      q1_q4 = sc(q1,q4)
      q2_q3 = sc(q2,q3)
      q2_q4 = sc(q2,q4)
      q3_q4 = sc(q3,q4)
      e1_e2 = sc(e1,e2)
 
      e1_e3 = sc(e1,e3)
      e1_e4 = sc(e1,e4)
 
 
 
 
      e2_e3 = sc(e2,e3)
      e2_e4 = sc(e2,e4)
 
      e3_e4 = sc(e3,e4)
 
 
      q1_e3 = sc(q1,e3)
      q1_e4 = sc(q1,e4)
      q2_e3 = sc(q2,e3)
      q2_e4 = sc(q2,e4)
      e1_q3 = sc(e1,q3)
      e1_q4 = sc(e1,q4)
      e2_q3 = sc(e2,q3)
      e2_q4 = sc(e2,q4)
      e3_q4 = sc(e3,q4)
      e4_q3 = sc(e4,q3)
 
 
      mz3=dsqrt(cdabs(q3_q3))
      mz4=dsqrt(cdabs(q4_q4))
 
 
      res = czero
 
 
      yyy1 = zprime_zz_1
      yyy2 = zprime_zz_2
 
 
       res= yyy1*( q1_e3*e1_e4 + q1_e4*e1_e3 )  +   yyy2*( et1(e1,e3,e4,q3) - et1(e1,e3,e4,q4) )
 
 
      end subroutine zprimezzampl
 
 
subroutine getdecay_couplings_spinors_props(VVMode,idordered,pordered,h3,h4, sp,pV)
   implicit none
   integer, intent(in) :: VVMode,idordered(6:9),h3,h4
   real(dp), intent(in) :: pordered(1:4,6:9)
   complex(dp), intent(out) :: sp(3:4,1:4)
   real(dp), intent(out) :: pV(3:4,1:4)
   real(dp) :: s, aL1,aR1,aL2,aR2
   complex(dp) :: propV(1:2)
 
   !        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
      pv(3,:) = pordered(:,6)+pordered(:,7)
      pv(4,:) = pordered(:,8)+pordered(:,9)
      sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3)  ! ubar(l1), v(l2)
      sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
      sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4)  ! ubar(l3), v(l4)
      sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
      s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
      propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
 
   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
      pv(3,:) = pordered(:,6)+pordered(:,7)
      pv(4,:) = pordered(:,8)+pordered(:,9)
      sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3)  ! ubar(l1), v(l2)
      sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
      sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4)  ! ubar(l3), v(l4)
      sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
      s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
      propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
 
   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=1d0
      ar2=1d0
      pv(3,:) = pordered(:,6)+pordered(:,7)
      pv(4,:) = pordered(:,8)
      sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3)  ! ubar(l1), v(l2)
      sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
      sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out')  ! photon
      !sp(4,1:4)=pV(4,1:4); print *, "this checks gauge invariance"
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
      propv(2)=1d0
 
   elseif( vvmode.eq.ggmode ) then
   !        gamma gamma DECAYS
      al1=1d0
      ar1=1d0
      al2=1d0
      ar2=1d0
      pv(3,:) = pordered(:,6)
      pv(4,:) = pordered(:,8)
      sp(3,:) = pol_mless2(dcmplx(pordered(:,6)),h3,'out')  ! photon
      sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out')  ! photon
      !sp(3,1:4)=pV(3,1:4); print *, "this checks gauge invariance"
      !sp(4,1:4)=pV(4,1:4)
      propv(1)=1d0
      propv(2)=1d0
 
   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=1d0
      ar2=1d0
      pv(3,:) = pordered(:,6)+pordered(:,7)
      pv(4,:) = pordered(:,8)
      sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3)  ! ubar(l1), v(l2)
      sp(3,:) = -sp(3,:) ! photon propagator
      sp(4,:) = pol_mless2(dcmplx(pordered(:,8)),h4,'out')  ! photon
      !sp(4,1:4)=pV(4,1:4); print *, "this checks gauge invariance"
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      propv(1) = 1d0
      propv(2) = 1d0
      if( s.lt.mphotoncutoff**2 ) propv(1)=czero
 
   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
      pv(3,:) = pordered(:,6)+pordered(:,7)
      pv(4,:) = pordered(:,8)+pordered(:,9)
      sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3)  ! ubar(l1), v(l2)
      sp(3,:) = -sp(3,:)
      sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4)  ! ubar(l3), v(l4)
      sp(4,:) = -sp(4,:) + pv(4,:)*( sc(sp(4,:),dcmplx(pv(4,:))) )/scr(pv(4,:),pv(4,:))! full propagator numerator
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      propv(1) = 1d0! = s/dcmplx(s)
      if( s.lt.mphotoncutoff**2 ) propv(1)=czero
      s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
      propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
 
   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
      pv(3,:) = pordered(:,6)+pordered(:,7)
      pv(4,:) = pordered(:,8)+pordered(:,9)
      sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3)  ! ubar(l1), v(l2)
      sp(3,:) = -sp(3,:) + pv(3,:)*( sc(sp(3,:),dcmplx(pv(3,:))) )/scr(pv(3,:),pv(3,:))! full propagator numerator
      sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4)  ! ubar(l3), v(l4)
      sp(4,:) = -sp(4,:)
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
      s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
      propv(2) = 1d0 ! = s/dcmplx(s)
      if( s.lt.mphotoncutoff**2 ) propv(2)=czero
 
   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
      pv(3,:) = pordered(:,6)+pordered(:,7)
      pv(4,:) = pordered(:,8)+pordered(:,9)
      sp(3,:) = pol_dk2mom(dcmplx(pordered(:,6)),dcmplx(pordered(:,7)),h3)  ! ubar(l1), v(l2)
      sp(3,:) = -sp(3,:)
      sp(4,:) = pol_dk2mom(dcmplx(pordered(:,8)),dcmplx(pordered(:,9)),h4)  ! ubar(l3), v(l4)
      sp(4,:) = -sp(4,:)
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      propv(1) = 1d0 ! = s/dcmplx(s)
      if( s.lt.mphotoncutoff**2 ) propv(1)=czero
      s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
      propv(2) = 1d0 ! = s/dcmplx(s)
      if( s.lt.mphotoncutoff**2 ) propv(2)=czero
 
   else
      call error("Unsupported decay modes")
   endif
 
   sp(3,:) = sp(3,:)*propv(1)
   sp(4,:) = sp(4,:)*propv(2)
   if (h3.eq.-1) then
      sp(3,:) = al1 * sp(3,:)
   elseif(h3.eq.1) then
      sp(3,:) = ar1 * sp(3,:)
   endif
   if (h4.eq.-1) then
      sp(4,:) = al2 * sp(4,:)
   elseif(h4.eq.1) then
      sp(4,:) = ar2 * sp(4,:)
   endif
 
   return
end subroutine
 
subroutine getdecay_vvmode_ordering(MY_IDUP, VVMode,ordering,ordering_swap)
   implicit none
   integer, intent(in) :: MY_IDUP(6:9)
   integer, intent(out) :: VVMode,ordering(1:4),ordering_swap(1:4)
   integer :: idV(1:2)
 
   ordering=(/3,4,5,6/)
   idv(1)=coupledvertex(my_idup(6:7),-1)
   idv(2)=coupledvertex(my_idup(8:9),-1)
   if(my_idup(6).eq.pho_ .or. my_idup(7).eq.pho_) idv(1)=pho_
   if(my_idup(8).eq.pho_ .or. my_idup(9).eq.pho_) idv(2)=pho_
   if(convertlhe(my_idup(6)).lt.0 .or. my_idup(6).eq.not_a_particle_) then
      call swap(ordering(1),ordering(2))
   endif
   if(convertlhe(my_idup(8)).lt.0 .or. my_idup(8).eq.not_a_particle_) then
      call swap(ordering(3),ordering(4))
   endif
   if( &
         (idv(1).eq.wm_ .and. idv(2).eq.wp_) .or. &
         (idv(2).eq.z0_ .and. idv(1).eq.pho_) &
     ) then
      call swap(ordering(1),ordering(3))
      call swap(ordering(2),ordering(4))
      call swap(idv(1),idv(2))
   endif
   ordering_swap(:)=ordering(:)
   call swap(ordering_swap(1),ordering_swap(3))
 
   if(idv(1).eq.z0_ .and. idv(2).eq.z0_) then
      vvmode=zzmode
   elseif(idv(1).eq.z0_ .and. idv(2).eq.pho_) then
      vvmode=zgmode
   elseif(idv(1).eq.pho_ .and. idv(2).eq.pho_) then
      vvmode=ggmode
   elseif(idv(1).eq.wp_ .and. idv(2).eq.wm_) then
      vvmode=wwmode
   else
      print *,"idV=",idv
      call error("Unsupported decay Modes")
   endif
   return
end subroutine
 
 
       end module