mod_Higgs.F90

Go to the documentation of this file.

      MODULE modhiggs
      use modparameters
      use modmisc
      implicit none
      private
 
 
!----- notation for subroutines
      public :: evalamp_gg_h_vv
      public :: evalamp_h_vv
      public :: evalamp_h_ff, evalamp_h_tt_decay
      public :: calchelamp,calchelamp2
 
      CONTAINS
 
 
!----- a subroutinefor gg -> H -> (Z+gamma*)(Z+gamma*)/WW/gammagamma
!----- all outgoing convention and the following momentum assignment
!-----  0 -> g(p1) + g(p2) + e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
      SUBROUTINE evalamp_gg_h_vv(p,MY_IDUP,res)
      implicit none
      real(dp), intent(out) ::  res
      real(dp), intent(in) :: p(4,6)
      integer, intent(in) :: my_idup(6:9)
      complex(dp) :: a_vv(1:18), a0_vv(1:2)
      integer :: i1,i2,i3,i4,vvmode,vvmode_swap
      real(dp) :: prefactor!,res2
      real(dp) :: intcolfac
      integer :: ordering(1:4),ordering_swap(1:4)
      logical :: dointerference
 
      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,vvmode_swap,ordering_swap)
 
! Global normalization
      if( vvmode.eq.zzmode ) then!  Z decay
         prefactor = 8d0*overallcouplvffsq**2
      elseif( vvmode.eq.wwmode ) then !  W decay
         prefactor = 8d0*overallcouplvffsq**2
      elseif( vvmode.eq.zgmode ) then !  Z+photon "decay"
         prefactor = 8d0*overallcouplvffsq ! Only single powers
      elseif( vvmode.eq.ggmode ) then !  photon "decay"
         prefactor = 8d0
      else
         prefactor = 0d0
      endif
      prefactor = prefactor * (alphas/(3.0_dp*pi*vev))**2
 
 
! ! MADGRAPH CHECK
! res=0d0
! if (MY_IDUP(6).ne.MY_IDUP(8) ) return
! if (MY_IDUP(7).ne.MY_IDUP(9) ) return
! if (MY_IDUP(6).eq.MY_IDUP(8) ) return
! if (MY_IDUP(7).eq.MY_IDUP(9) ) return
! print *, "MY COUPL",al1*dsqrt(overallCouplVffsq),ar1*dsqrt(overallCouplVffsq)
! 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
! al1=+0.210150647402957d0! this is MadGraphs gzl(L)
! al2=al1
! ar2=ar1
! print *, "MG COUPL",al1,ar1
! pause
 
          res = zero
          a_vv(:) = 0d0
          dointerference = includeinterference .and. (         &
            ((vvmode.eq.zzmode) .and. (vvmode_swap.eq.zzmode)) &
            )
          if ( includevprime .and. .not.(vvmode.eq.zzmode .or. vvmode.eq.zgmode .or. vvmode.eq.wwmode) ) then
              call error("Contact terms only for WW, ZZ or Zg!")
          endif
          do i1=1,2;  do i2=1,2;  do i3=1,2;  do i4=1,2!  sum over helicities
                  call calchelamp(ordering,vvmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(1))
                  if( vvmode.eq.zzmode ) then
                      if( includegammastar ) then
                         call calchelamp(ordering,zgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(3))
                         call calchelamp(ordering,gszmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(5))
                         call calchelamp(ordering,gsgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(7))
                      endif
                      if( includevprime ) then
                          call calchelamp(ordering,zzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(9))
                          call calchelamp(ordering,zpzmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(11))
                          call calchelamp(ordering,zpzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(13))
                      endif
                      if( includegammastar .and. includevprime ) then
                          call calchelamp(ordering,gszpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(15))
                          call calchelamp(ordering,zpgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(17))
                      endif
                  elseif( vvmode.eq.zgmode ) then
                      if(includegammastar) then
                          call calchelamp(ordering,gsgmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(3))
                      endif
                      if( includevprime ) then
                          call calchelamp(ordering,zpgmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(5))
                      endif
                  elseif( vvmode.eq.wwmode .and. includevprime ) then
                      call calchelamp(ordering,wwpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(9))
                      call calchelamp(ordering,wpwmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(11))
                      call calchelamp(ordering,wpwpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(13))
                  endif
 
                  if( dointerference ) then
                      call calchelamp(ordering_swap,vvmode_swap,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(2))
                      if( includegammastar ) then
                          call calchelamp(ordering_swap,zgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(4))
                          call calchelamp(ordering_swap,gszmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(6))
                          call calchelamp(ordering_swap,gsgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(8))
                      endif
                      if( includevprime ) then
                          call calchelamp(ordering_swap,zzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(10))
                          call calchelamp(ordering_swap,zpzmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(12))
                          call calchelamp(ordering_swap,zpzpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(14))
                      endif
                      if( includegammastar .and. includevprime ) then
                          call calchelamp(ordering_swap,gszpmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(16))
                          call calchelamp(ordering_swap,zpgsmode,my_idup,p(1:4,1:6),i1,i2,i3,i4,a_vv(18))
                      endif
                  endif
 
                  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
                  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
                  res = res + dreal(a0_vv(1)*dconjg(a0_vv(1)))
                  res = res + dreal(a0_vv(2)*dconjg(a0_vv(2)))
                  if( dointerference .and. (i3.eq.i4) ) then! interfere the 3456 with 5436 pieces
                      res = res - 2d0*intcolfac*dreal(  a0_vv(1)*dconjg(a0_vv(2))  ) ! minus from Fermi statistics
                  endif
          enddo;  enddo;  enddo;  enddo
 
 
! MADGRAPH CHECK
! call coupsm(0)
! if( (MY_IDUP(6).eq.MY_IDUP(8)) .and. (MY_IDUP(7).eq.MY_IDUP(9)) ) then
!       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)
! else
!       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)
!       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)
! endif
! 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)) )
! res2=res2*GeV**4
! pause
! res=res2; RETURN
 
! res= res*prefactor
! print *, "checker 1",res
! print *, "checker 2",res2
! print *, "checker 1/2",res/res2
! pause
 
 
          res = res*prefactor
          if( (vvmode.eq.zzmode) .and. dointerference ) res = res * symmfac
 
         !print *,"VVmode:",VVmode
         !print *,"VVmode_swap:",VVmode_swap
         !print *,"ids:",MY_IDUP
         !print *,"res:",res
         !pause
 
      RETURN
      END SUBROUTINE
 
   subroutine calchelamp(ordering,VVMode,MY_IDUP,p,i1,i2,i3,i4,A)
      implicit none
      integer :: ordering(1:4),vvmode,i1,i2,i3,i4,l1,l2,l3,l4,my_idup(6:9)
      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_reso)
 
      pin(1,:) = p(:,1)
      pin(2,:) = p(:,2)
      sp(1,:) = pol_mless2(dcmplx(p(:,1)),-3+2*i1,'in')  ! gluon
      sp(2,:) = pol_mless2(dcmplx(p(:,2)),-3+2*i2,'in')  ! gluon
      !sp(1,1:4)=pin(1,1:4);print *, "this checks IS gauge invariance"
      !sp(2,1:4)=pin(2,1:4);print *, "this checks IS gauge invariance"
      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 gghzzampl(vvmode,pin,sp,a(1))
      a(1) = a(1) * propg
   end subroutine
 
      SUBROUTINE gghzzampl(VVMode,p,sp,res)
      implicit none
      integer, intent(in) :: VVMode
      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) :: 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) :: xxx1,xxx2,xxx3,yyy1,yyy2,yyy3,yyy4
      complex(dp) :: ghz1_dyn,ghz2_dyn,ghz3_dyn,ghz4_dyn
      real(dp) :: q34
      real(dp) :: q_q, q3_q3, q4_q4
 
 
      res = 0d0
      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)
      q3_q3 = sc(q3,q3)
      q4_q4 = sc(q4,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)
 
 
      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
 
!---- data that defines couplings
      if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode)  ) then! decay ZZ's or WW's
         ghz1_dyn = hvvspinzerodynamiccoupling(1,q3_q3,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(2,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(3,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(4,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.gszmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q3_q3,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q3_q3,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q3_q3,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q3_q3,q_q)
      elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.zgsmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q4_q4,q_q)
      elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode)  .or. (vvmode.eq.gsgmode) ) then
         ghz1_dyn = czero
         ghz2_dyn = hvvspinzerodynamiccoupling(9,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(10,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(11,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.zzpmode) .or. (vvmode.eq.wwpmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(12,q3_q3,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(13,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(14,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(15,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.zpzmode) .or. (vvmode.eq.wpwmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(12,q4_q4,q3_q3,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(13,q4_q4,q3_q3,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(14,q4_q4,q3_q3,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(15,q4_q4,q3_q3,q_q)
      elseif( (vvmode.eq.zpzpmode) .or. (vvmode.eq.wpwpmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(16,q3_q3,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(17,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(18,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(19,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.gszpmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q3_q3,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q3_q3,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q3_q3,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q3_q3,q_q)
      elseif( (vvmode.eq.zpgmode) .OR. (vvmode.eq.zpgsmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q4_q4,q_q)
      else
         ghz1_dyn = czero
         ghz2_dyn = czero
         ghz3_dyn = czero
         ghz4_dyn = czero
         print *,"VVMode",vvmode,"not implemented"
      endif
 
      if( .not. generate_as ) then
         xxx1 = ghg2+ghg3/4d0/lambda**2*q_q
         xxx3 = -2d0*ghg4
         yyy1 = ghz1_dyn*m_v**2/q_q &  ! in this line M_V is indeed correct, not a misprint
              + ghz2_dyn*(q_q-q3_q3-q4_q4)/q_q &
              + ghz3_dyn/lambda**2*(q_q-q3_q3-q4_q4)*(q_q-q4_q4-q3_q3)/4d0/q_q
         yyy2 = -2d0*ghz2_dyn-ghz3_dyn/2d0/lambda**2*(q_q-q3_q3-q4_q4)
         yyy3 = -2d0*ghz4_dyn
      else
         xxx1 = ahg1
         xxx3 = ahg3
         if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode)  ) then! decay ZZ's or WW's
            yyy1 = ahz1
            yyy2 = ahz2
            yyy3 = ahz3
         elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode)  .or. (vvmode.eq.gsgmode) ) then! decay (gamma-gamma) OR (gamma*-gamma*) OR (gamma*-gamma)
            yyy1 = ahz1
            yyy2 = -2*ahz1 !ahz2  ! gauge invariance fixes ahz2 in this case
            yyy3 = ahz3
         elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.gszmode) .OR. (vvmode.eq.zgsmode) ) then! decay (Z-photon) OR (gamma*-Z) OR (Z-gamma*)
            yyy1 = ahz1
            yyy2 = -2*ahz1*q_q/(q_q-q3_q3)
            yyy3 = ahz3
         else
            yyy1=czero
            yyy2=czero
            yyy3=czero
            print *,"VVMode",vvmode,"not implemented in generate_as"
         endif
      endif
 
      res= e1_e2*e3_e4*q_q**2*yyy1*xxx1                  &
         + e1_e2*e3_q4*e4_q3*q_q*yyy2*xxx1               &
         + et1(e1,e2,q1,q2)*e3_e4*q_q*yyy1*xxx3          &
         + et1(e1,e2,q1,q2)*e3_q4*e4_q3*yyy2*xxx3        &
         + et1(e1,e2,q1,q2)*et1(e3,e4,q3,q4)*yyy3*xxx3   &
         + et1(e3,e4,q3,q4)*e1_e2*q_q*yyy3*xxx1
 
  END SUBROUTINE gghzzampl
 
 
!----- a subroutinefor H -> (Z+gamma*)(Z+gamma*)/WW/gammagamma
!----- all outgoing convention and the following momentum assignment
!-----  0 -> Higgs(p1) + e-(p3) + e+(p4) +mu-(p5) +mu+(p6)
      subroutine evalamp_h_vv(p,MY_IDUP,res)
      implicit none
      real(dp), intent(out) ::  res
      real(dp), intent(in) :: p(4,6)
      integer, intent(in) :: my_idup(6:9)
      complex(dp) :: a_vv(1:18),a0_vv(1:2)
      integer :: i3,i4,vvmode,vvmode_swap
      real(dp) :: prefactor
      real(dp) :: intcolfac
      integer :: ordering(1:4),ordering_swap(1:4)
      logical :: dointerference
 
      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,vvmode_swap,ordering_swap)
 
! Global normalization
         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
 
! ! MADGRAPH CHECK
! res=0d0
! if (MY_IDUP(6).ne.MY_IDUP(8) ) return
! if (MY_IDUP(7).ne.MY_IDUP(9) ) return
! if (MY_IDUP(6).eq.MY_IDUP(8) ) return
! if (MY_IDUP(7).eq.MY_IDUP(9) ) return
! print *, "MY COUPL",al1*dsqrt(overallCouplVffsq),ar1*dsqrt(overallCouplVffsq)
! 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
! al1=+0.210150647402957d0! this is MadGraphs gzl(L)
! al2=al1
! ar2=ar1
! print *, "MG COUPL",al1,ar1
! pause
 
 
          res = zero
          a_vv(:) = 0d0
          dointerference = includeinterference .and. (         &
            ((vvmode.eq.zzmode) .and. (vvmode_swap.eq.zzmode)) &
            )
          if ( includevprime .and. .not.(vvmode.eq.zzmode .or. vvmode.eq.zgmode .or. vvmode.eq.wwmode) ) then
              call error("Contact terms only for WW, ZZ or Zg!")
          endif
          do i3=1,2;  do i4=1,2!  sum over helicities
                  call calchelamp2(ordering,vvmode,my_idup,p(1:4,1:6),i3,i4,a_vv(1))
                  if( vvmode.eq.zzmode ) then
                      if( includegammastar ) then
                         call calchelamp2(ordering,zgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(3))
                         call calchelamp2(ordering,gszmode,my_idup,p(1:4,1:6),i3,i4,a_vv(5))
                         call calchelamp2(ordering,gsgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(7))
                      endif
                      if( includevprime ) then
                          call calchelamp2(ordering,zzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(9))
                          call calchelamp2(ordering,zpzmode,my_idup,p(1:4,1:6),i3,i4,a_vv(11))
                          call calchelamp2(ordering,zpzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(13))
                      endif
                      if( includegammastar .and. includevprime ) then
                          call calchelamp2(ordering,gszpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(15))
                          call calchelamp2(ordering,zpgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(17))
                      endif
                  elseif( vvmode.eq.zgmode ) then
                      if(includegammastar) then
                          call calchelamp2(ordering,gsgmode,my_idup,p(1:4,1:6),i3,i4,a_vv(3))
                      endif
                      if( includevprime ) then
                          call calchelamp2(ordering,zpgmode,my_idup,p(1:4,1:6),i3,i4,a_vv(5))
                      endif
                  elseif( vvmode.eq.wwmode .and. includevprime ) then
                      call calchelamp2(ordering,wwpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(9))
                      call calchelamp2(ordering,wpwmode,my_idup,p(1:4,1:6),i3,i4,a_vv(11))
                      call calchelamp2(ordering,wpwpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(13))
                  endif
 
                  if( dointerference ) then
                      call calchelamp2(ordering_swap,vvmode,my_idup,p(1:4,1:6),i3,i4,a_vv(2))
                      if( includegammastar ) then
                          call calchelamp2(ordering_swap,zgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(4))
                          call calchelamp2(ordering_swap,gszmode,my_idup,p(1:4,1:6),i3,i4,a_vv(6))
                          call calchelamp2(ordering_swap,gsgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(8))
                      endif
                      if( includevprime ) then
                          call calchelamp2(ordering_swap,zzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(10))
                          call calchelamp2(ordering_swap,zpzmode,my_idup,p(1:4,1:6),i3,i4,a_vv(12))
                          call calchelamp2(ordering_swap,zpzpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(14))
                      endif
                      if( includegammastar .and. includevprime ) then
                          call calchelamp2(ordering_swap,gszpmode,my_idup,p(1:4,1:6),i3,i4,a_vv(16))
                          call calchelamp2(ordering_swap,zpgsmode,my_idup,p(1:4,1:6),i3,i4,a_vv(18))
                      endif
                  endif
 
                  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
                  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
                  res = res + dreal(a0_vv(1)*dconjg(a0_vv(1)))
                  res = res + dreal(a0_vv(2)*dconjg(a0_vv(2)))
                  if( dointerference .and. (i3.eq.i4) ) then! interfere the 3456 with 5436 pieces
                      res = res - 2d0*intcolfac*dreal(  a0_vv(1)*dconjg(a0_vv(2))  ) ! minus from Fermi statistics
                  endif
          enddo;  enddo
 
 
 
 
! MADGRAPH CHECK
! call coupsm(0)
! if( (MY_IDUP(6).eq.MY_IDUP(8)) .and. (MY_IDUP(7).eq.MY_IDUP(9)) ) then
!       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)
! else
!       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)
! endif
! 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)) )
! res2=res2*GeV**4
! pause
! res=res2; RETURN
 
! res= res*prefactor
! print *, "checker 1",res
! print *, "checker 2",res2
! print *, "checker 1/2",res/res2
! pause
 
 
          res = res*prefactor
          if( (vvmode.eq.zzmode) .and. dointerference ) res = res * symmfac
 
      RETURN
      END SUBROUTINE
 
   subroutine calchelamp2(ordering,VVMode,MY_IDUP,p,i3,i4,A)
   implicit none
   integer :: ordering(1:4),VVMode,i3,i4,l1,l2,l3,l4,MY_IDUP(6:9)
   real(dp) :: p(1:4,1:6)
   real(dp) :: pin(4,4)
   complex(dp) :: A(1:1), sp(3:4,4)
 
      l1=ordering(1)
      l2=ordering(2)
      l3=ordering(3)
      l4=ordering(4)
 
      pin(1,:) = p(:,1)
 
      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 hzzampl(vvmode,pin,sp,a(1))
   end subroutine
 
      subroutine hzzampl(VVMode,p,sp,res)
      implicit none
      integer, intent(in) :: VVMode
      real(dp), intent(in) :: p(4,4)
      complex(dp), intent(in) :: sp(3:4,4)
      complex(dp), intent(out) :: res
      complex(dp) :: e3_e4
      complex(dp) :: e3_q4,e4_q3
      complex(dp) :: q1(4),q3(4),q4(4),q(4)
      complex(dp) :: e3(4),e4(4)
      complex(dp) :: yyy1,yyy2,yyy3,yyy4
      real(dp) :: q34
      real(dp) :: q_q, q3_q3, q4_q4
      complex(dp) :: ghz1_dyn,ghz2_dyn,ghz3_dyn,ghz4_dyn
 
 
 
      res = 0d0
      q1 = dcmplx(p(1,:),0d0)
      q3 = dcmplx(p(3,:),0d0)
      q4 = dcmplx(p(4,:),0d0)
 
 
      e3 = sp(3,:)
      e4 = sp(4,:)
 
      q = -q1
 
      q_q =sc(q,q)
      q3_q3 = sc(q3,q3)
      q4_q4 = sc(q4,q4)
 
      e3_e4 = sc(e3,e4)
      e3_q4 = sc(e3,q4)
      e4_q3 = sc(e4,q3)
 
 
      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
 
 
!---- data that defines couplings
      if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode)  ) then! decay ZZ's or WW's
         ghz1_dyn = hvvspinzerodynamiccoupling(1,q3_q3,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(2,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(3,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(4,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.gszmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q3_q3,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q3_q3,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q3_q3,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q3_q3,q_q)
      elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.zgsmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q4_q4,q_q)
      elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode)  .or. (vvmode.eq.gsgmode) ) then
         ghz1_dyn = czero
         ghz2_dyn = hvvspinzerodynamiccoupling(9,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(10,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(11,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.zzpmode) .or. (vvmode.eq.wwpmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(12,q3_q3,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(13,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(14,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(15,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.zpzmode) .or. (vvmode.eq.wpwmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(12,q4_q4,q3_q3,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(13,q4_q4,q3_q3,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(14,q4_q4,q3_q3,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(15,q4_q4,q3_q3,q_q)
      elseif( (vvmode.eq.zpzpmode) .or. (vvmode.eq.wpwpmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(16,q3_q3,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(17,q3_q3,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(18,q3_q3,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(19,q3_q3,q4_q4,q_q)
      elseif( (vvmode.eq.gszpmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q3_q3,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q3_q3,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q3_q3,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q3_q3,q_q)
      elseif( (vvmode.eq.zpgmode) .OR. (vvmode.eq.zpgsmode) ) then
         ghz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q4_q4,q_q)
         ghz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q4_q4,q_q)
         ghz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q4_q4,q_q)
         ghz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q4_q4,q_q)
      else
         ghz1_dyn = czero
         ghz2_dyn = czero
         ghz3_dyn = czero
         ghz4_dyn = czero
         print *,"VVMode",vvmode,"not implemented"
      endif
 
      if( .not. generate_as ) then
         yyy1 = ghz1_dyn*m_v**2/q_q &  ! in this line M_V is indeed correct, not a misprint
              + ghz2_dyn*(q_q-q3_q3-q4_q4)/q_q &
              + ghz3_dyn/lambda**2*(q_q-q3_q3-q4_q4)*(q_q-q4_q4-q3_q3)/4d0/q_q
         yyy2 = -2d0*ghz2_dyn-ghz3_dyn/2d0/lambda**2*(q_q-q3_q3-q4_q4)
         yyy3 = -2d0*ghz4_dyn
      else
         if( (vvmode.eq.zzmode) .or. (vvmode.eq.wwmode)  ) then! decay ZZ's or WW's
            yyy1 = ahz1
            yyy2 = ahz2
            yyy3 = ahz3
         elseif( (vvmode.eq.ggmode) .or. (vvmode.eq.gsgsmode)  .or. (vvmode.eq.gsgmode) ) then! decay (gamma-gamma) OR (gamma*-gamma*) OR (gamma*-gamma)
            yyy1 = ahz1
            yyy2 = -2*ahz1 !ahz2  ! gauge invariance fixes ahz2 in this case
            yyy3 = ahz3
         elseif( (vvmode.eq.zgmode) .OR. (vvmode.eq.gszmode) .OR. (vvmode.eq.zgsmode) ) then! decay (Z-photon) OR (gamma*-Z) OR (Z-gamma*)
            yyy1 = ahz1
            yyy2 = -2*ahz1*q_q/(q_q-q3_q3)
            yyy3 = ahz3
         else
            yyy1=czero
            yyy2=czero
            yyy3=czero
            print *,"VVMode",vvmode,"not implemented in generate_as"
         endif
      endif
 
      res = e3_e4*q_q*yyy1                  &
          + e3_q4*e4_q3*yyy2                      &
          + et1(e3,e4,q3,q4)*yyy3
 
 
   END SUBROUTINE hzzampl
 
 
! Higgs decay to tau^+ tau^-   or    top anti-top
! Decay amplitude H --> tau^-(p1) + tau^+(p2)
! or              H --> tbar(p1)  + top(p2)
! with tau/top controlled by value of mass_F
! Since this is an isotropic scalar decay, just provide a matrix element SQUARED
! R.Rontsch July 2015
   SUBROUTINE evalamp_h_ff(pin,mass_F,res)
   implicit none
   real(dp), intent(out) ::  res
   complex(dp) :: amp2
   real(dp), intent(in) :: pin(1:4,1:2),mass_f
   real(dp)             :: s12
 
      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
      amp2 =   2d0*s12*(kappa_tilde**2 + kappa**2) - 8d0*mass_f**2*kappa**2
      amp2 = amp2*mass_f**2/vev**2
      res = cdabs(amp2)
 
   RETURN
   END SUBROUTINE
 
 
! Decay amplitude H --> tau^-(-->l^-(p1)+nubar(p2)+nutau(p3))  +  tau^+(-->nu(p4)+l^+(p5)+nutaubar(p6))
! or H --> tbar^-(-->l^-(p1)+nubar(p2)+bbar(p3))+top(-->nu(p4)+l^+(p5)+b(p6))
! Note: 1. Value of m allows to change between H->tau^+tau^- and H->ttbar
!       2. Full kinematics of decaying particles -- no NWA.
!       3. Final state b quark is massless.
!       4. At present, no width in the tau/top propagator
! R. Rontsch July 2015
   SUBROUTINE evalamp_h_tt_decay(pin,mass_F,ga_F,res)
   implicit none
   real(dp), intent(out) ::  res
   real(dp), intent(in) :: pin(1:4,1:6),mass_f,ga_f
   integer              :: j
   real(dp)             :: p(1:6,1:4),s12,s45,s123,s456,s(6,6)
   complex(dp)          :: za(6,6),zb(6,6),amp,kl,kr
 
!       p=pin
      do j=1,6
          call convert_to_mcfm(pin(1:4,j),p(j,1:4))
      enddo
!       call spinoru(6,p,za,zb,s)
      call spinoru(p,za,zb,s)
 
      s12=s(1,2)
      s45=s(4,5)
      s123=s(1,2)+s(1,3)+s(2,3)
      s456=s(4,5)+s(4,6)+s(5,6)
 
      kl=-mass_f/vev*( kappa -(0d0,1d0)*kappa_tilde )
      kr=-mass_f/vev*( kappa +(0d0,1d0)*kappa_tilde )
 
      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)) &
             + 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))
 
! overall factors and propagators
      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)
      amp=amp*16d0*ci*mass_f*gwsq**2
      res = cdabs(amp)**2
 
   RETURN
   END SUBROUTINE
 
 
 
 
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
   complex(dp) :: propV(1:2), 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
      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
 
 
   elseif( vvmode.eq.zpzmode ) then
   !        Z'Z DECAYS
      al1 = vpffcoupling(idordered(6),-1,.false.)
      ar1 = vpffcoupling(idordered(6),+1,.false.)
      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))
      if( m_vprime .gt. 0d0 ) then
        propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(1) = 1d0
      else
        propv(1) = s/m_z**2
      endif
      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.zzpmode ) 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
      al2 = vpffcoupling(idordered(8),-1,.false.)
      ar2 = vpffcoupling(idordered(8),+1,.false.)
 
      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))
      if( m_vprime .gt. 0d0 ) then
        propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(2) = 1d0
      else
        propv(2) = s/m_z**2
      endif
 
 
   elseif( vvmode.eq.zpzpmode ) then
   !        Z'Z' DECAYS
      al1 = vpffcoupling(idordered(6),-1,.false.)
      ar1 = vpffcoupling(idordered(6),+1,.false.)
      al2 = vpffcoupling(idordered(8),-1,.false.)
      ar2 = vpffcoupling(idordered(8),+1,.false.)
 
      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))
      if( m_vprime .gt. 0d0 ) then
        propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(1) = 1d0
      else
        propv(1) = s/m_z**2
      endif
      s = scr(pordered(:,8)+pordered(:,9),pordered(:,8)+pordered(:,9))
      if( m_vprime .gt. 0d0 ) then
        propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(2) = 1d0
      else
        propv(2) = s/m_z**2
      endif
 
 
   elseif( vvmode.eq.zpgsmode ) then
   !        Z' gamma* DECAYS
      al1 = vpffcoupling(idordered(6),-1,.false.)
      ar1 = vpffcoupling(idordered(6),+1,.false.)
      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))
      if( m_vprime .gt. 0d0 ) then
        propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(1) = 1d0
      else
        propv(1) = s/m_z**2
      endif
      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.gszpmode ) 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
      al2 = vpffcoupling(idordered(8),-1,.false.)
      ar2 = vpffcoupling(idordered(8),+1,.false.)
 
      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))
      if( m_vprime .gt. 0d0 ) then
        propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(2) = 1d0
      else
        propv(2) = s/m_z**2
      endif
 
 
   elseif( vvmode.eq.zpgmode ) then
   !        Z' gamma DECAYS
      al1 = vpffcoupling(idordered(6),-1,.false.)
      ar1 = vpffcoupling(idordered(6),+1,.false.)
      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
      s = scr(pordered(:,6)+pordered(:,7),pordered(:,6)+pordered(:,7))
      if( m_vprime .gt. 0d0 ) then
        propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(1) = 1d0
      else
        propv(1) = s/m_z**2
      endif
      propv(2)=1d0
 
 
   elseif( vvmode.eq.wpwmode ) then
   !        W'W DECAYS
      al1 = vpffcoupling(idordered(6),-1,.true.)
      ar1 = vpffcoupling(idordered(6),+1,.true.)
      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(pv(3,:),pv(3,:))
      if( m_vprime .gt. 0d0 ) then
        propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(1) = 1d0
      else
        propv(1) = s/m_w**2
      endif
      s = scr(pv(4,:),pv(4,:))
      propv(2) = s/dcmplx(s - m_v**2,m_v*ga_v)
 
 
   elseif( vvmode.eq.wwpmode ) 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
      al2 = vpffcoupling(idordered(8),-1,.true.)
      ar2 = vpffcoupling(idordered(8),+1,.true.)
 
      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(pv(3,:),pv(3,:))
      propv(1) = s/dcmplx(s - m_v**2,m_v*ga_v)
      s = scr(pv(4,:),pv(4,:))
      if( m_vprime .gt. 0d0 ) then
        propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(2) = 1d0
      else
        propv(2) = s/m_w**2
      endif
 
 
   elseif( vvmode.eq.wpwpmode ) then
   !        W'W' DECAYS
      al1 = vpffcoupling(idordered(6),-1,.true.)
      ar1 = vpffcoupling(idordered(6),+1,.true.)
      al2 = vpffcoupling(idordered(8),-1,.true.)
      ar2 = vpffcoupling(idordered(8),+1,.true.)
 
      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(pv(3,:),pv(3,:))
      if( m_vprime .gt. 0d0 ) then
        propv(1) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(1) = 1d0
      else
        propv(1) = s/m_w**2
      endif
      s = scr(pv(4,:),pv(4,:))
      if( m_vprime .gt. 0d0 ) then
        propv(2) = s/dcmplx(s - m_vprime**2,m_vprime*ga_vprime)
      elseif( m_vprime .eq. 0d0 ) then
        propv(2) = 1d0
      else
        propv(2) = s/m_w**2
      endif
 
 
   else
      call error("Unsupported decay modes")
   endif
 
   !print *,"sp(3)=",sp(3,:)
   !print *,"propV(1)=",propV(1)
   !print *,"aL1,aR1=",aL1,aR1
   !print *,"sp(4)=",sp(4,:)
   !print *,"propV(2)=",propV(2)
   !print *,"aL2,aR2=",aL2,aR2
 
   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,VVMode_swap,ordering_swap)
   implicit none
   integer, intent(in) :: MY_IDUP(6:9)
   integer, intent(out) :: VVMode,ordering(1:4),VVMode_swap,ordering_swap(1:4)
   integer :: idV(1:2),idV_swap(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))
 
   idv_swap(1) = coupledvertex( (/ my_idup(3+ordering_swap(1)), my_idup(3+ordering_swap(2)) /), -1)
   idv_swap(2) = coupledvertex( (/ my_idup(3+ordering_swap(3)), my_idup(3+ordering_swap(4)) /), -1)
   if ( (idv_swap(1).eq.wm_) .and. (idv_swap(2).eq.wp_) ) then
      call swap(ordering_swap(1),ordering_swap(3))
      call swap(ordering_swap(2),ordering_swap(4))
      call swap(idv_swap(1),idv_swap(2))
   endif
 
   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 mode")
   endif
 
   vvmode_swap=invalidmode
   if(idv_swap(1).eq.z0_ .and. idv_swap(2).eq.z0_) then
      vvmode_swap=zzmode
   elseif(idv_swap(1).eq.wp_ .and. idv_swap(2).eq.wm_) then
      vvmode_swap=wwmode
   endif
   return
end subroutine
 
 
END MODULE