mod_VHiggs.F90

Go to the documentation of this file.

!--YaofuZhou-----------------------------------------
module modvhiggs
  use modparameters
  use modmisc
  implicit none
  private
 
 
  !----- notation for subroutines
  public :: evalamp_vhiggs
 
contains
 
subroutine evalamp_vhiggs(id,helicity,MomExt,me2)
      integer, intent(in) :: id(9)
      real(8), intent(in) :: helicity(9)
      real(8), intent(in) :: momext(1:4,1:9)
      real(8), intent(out) :: me2
      real(8) :: mass(3:5,1:2)
      integer :: i
      complex(8) amplitude, a_vv(1:4), amptest
      integer :: idin(9)
      real(8) :: helin(9)
      real(8) :: pin(1:4,1:9)
 
      idin(:)=id(:)
      helin(:)=helicity(:)
      pin(:,:)=momext(:,:)
      if(id(2).eq.convertlhe(pho_)) then
         call swap(idin(1),idin(2))
         call swap(helin(1),helin(2))
         call swap(pin(:,1),pin(:,2))
      endif
      if(id(7).eq.convertlhe(pho_)) then
         call swap(idin(6),idin(7))
         call swap(helin(6),helin(7))
         call swap(pin(:,6),pin(:,7))
      endif
 
      do i=3,5
        mass(i,1) = getmass(convertlhereverse(idin(i)))
        mass(i,2) = getdecaywidth(convertlhereverse(idin(i)))
      enddo
 
      a_vv(:)=czero
      if(idin(1).ne.convertlhe(pho_) .and. idin(6).ne.convertlhe(pho_)) then
         !print *,"Case1:"
         a_vv(1)=matrixelement0(pin,mass,helin,idin,(/.false., .false./))
         if(includegammastar) then
            !print *,"Case2:"
            a_vv(2) = matrixelement0(pin,mass,helin,idin,(/.false., .true./))
            !print *,"Case3:"
            a_vv(3) = matrixelement0(pin,mass,helin,idin,(/.true., .false./))
            !print *,"Case4:"
            a_vv(4) = matrixelement0(pin,mass,helin,idin,(/.true., .true./))
         endif
      else if(idin(1).eq.convertlhe(pho_) .and. idin(6).eq.convertlhe(pho_)) then
         !print *,"Case5:"
         a_vv(1)=matrixelement0(pin,mass,helin,idin,(/.true., .true./))
      else if(idin(1).eq.convertlhe(pho_)) then
         !print *,"Case6:"
         a_vv(1)=matrixelement0(pin,mass,helin,idin,(/.true., .false./))
         if(includegammastar) then
            !print *,"Case7:"
            a_vv(2) = matrixelement0(pin,mass,helin,idin,(/.true., .true./))
         endif
      else !if(idin(6).eq.convertLHE(Pho_)) then
         !print *,"Case8:"
         a_vv(1)=matrixelement0(pin,mass,helin,idin,(/.false., .true./))
         if(includegammastar) then
            !print *,"Case9:"
            a_vv(2) = matrixelement0(pin,mass,helin,idin,(/.true., .true./))
         endif
      endif
      amplitude = a_vv(1)+a_vv(2)+a_vv(3)+a_vv(4)
 
      ! XCHECK FROM DECAY ME
      !print *,pin(:,1)
      !print *,pin(:,2)
      !print *,pin(:,6)
      !print *,pin(:,7)
      !print *,idin(1)
      !print *,idin(2)
      !print *,idin(3)
      !print *,idin(4)
      !print *,idin(6)
      !print *,idin(7)
      !print *,helin(1)
      !print *,helin(2)
      !print *,helin(6)
      !print *,helin(7)
      !print *,amplitude
      !amptest = MATRIXELEMENT02(pin,mass,helin,idin)
      !print *,amptest
      !pause
 
      me2=dble(amplitude*dconjg(amplitude))
      return
end subroutine evalamp_vhiggs
 
!MATRIXELEMENT0.F
!VERSION 20160511
      complex(8) function matrixelement0(MomExt,mass,helicity,id,useA)
      implicit none
      real(8), intent(in) :: momext(1:4,1:9)
      real(8), intent(in) :: mass(3:5,1:2)
      real(8), intent(in) :: helicity(9)
      integer, intent(in) :: id(9)
      logical, intent(in) :: usea(2)
 
      integer mu3,mu4
      real(8) qq,q3_q3,q4_q4,q5_q5
      complex(8) pvvx0p
      complex(8) vcurrent1(4), acurrent1(4), current1(4), currentvp1(4)
      complex(8) vcurrent2(4), acurrent2(4), current2(4), currentvp2(4)
      complex(8) vpffcoupl(2,2)
      complex(8) pol1(3,4), pol2(3,4)
      complex(8) g_mu_nu(4,4), pp(4,4), epp(4,4)
      complex(8) prop1, prop2, prop3
      complex(8) prop_vp1, prop_vp2
      complex(8) gffz, gffa, gffw
      complex(8) gvvp, gvvs1, gvvs2
      complex(8) ghz1_dyn,ghz2_dyn,ghz3_dyn,ghz4_dyn
      complex(8) gvvpp, gvvps1, gvvps2
      complex(8) ghzzp1_dyn,ghzzp2_dyn,ghzzp3_dyn,ghzzp4_dyn
      complex(8) gvpvp, gvpvs1, gvpvs2
      complex(8) ghzpz1_dyn,ghzpz2_dyn,ghzpz3_dyn,ghzpz4_dyn
      complex(8) gvpvpp, gvpvps1, gvpvps2
      complex(8) ghzpzp1_dyn,ghzpzp2_dyn,ghzpzp3_dyn,ghzpzp4_dyn
 
      matrixelement0=czero
      if( &
         (id(1).ne.convertlhe(pho_) .and. usea(1) .and. .not.includegammastar) .or. &
         (id(6).ne.convertlhe(pho_) .and. usea(2) .and. .not.includegammastar) .or. &
         ((id(1)+id(2)).ne.0 .and. id(1).ne.convertlhe(pho_) .and. usea(1)) .or. &
         ((id(6)+id(7)).ne.0 .and. id(6).ne.convertlhe(pho_) .and. usea(2))      &
        ) then
         return
      endif
 
      vcurrent1 = czero
      acurrent1 = czero
      vcurrent2 = czero
      acurrent2 = czero
      vpffcoupl=czero
      currentvp1=czero
      currentvp2=czero
      prop_vp1=czero
      prop_vp2=czero
      ghzzp1_dyn = czero
      ghzzp2_dyn = czero
      ghzzp3_dyn = czero
      ghzzp4_dyn = czero
      ghzpz1_dyn = czero
      ghzpz2_dyn = czero
      ghzpz3_dyn = czero
      ghzpz4_dyn = czero
      ghzpzp1_dyn = czero
      ghzpzp2_dyn = czero
      ghzpzp3_dyn = czero
      ghzpzp4_dyn = czero
      gvvpp = czero
      gvvps1 = czero
      gvvps2 = czero
      gvpvp = czero
      gvpvs1 = czero
      gvpvs2 = czero
      gvpvpp = czero
      gvpvps1 = czero
      gvpvps2 = czero
 
      gffz = ci*2d0*dsqrt(couplzffsq) ! = sqrt(gwsq/(1.0_dp-xw))
      gffa = -ci*dsqrt(couplaffsq) ! = sqrt(gwsq*xw)
      gffw = ci*dsqrt(couplwffsq) ! = sqrt(gwsq/2.0_dp)
 
      qq = -scr(momext(:,3),momext(:,4))
      if (id(1).eq.convertlhe(pho_) .and. usea(1)) then
         q3_q3 = 0d0
      else
         q3_q3 = scr(momext(:,3),momext(:,3))
      endif
      if (id(6).eq.convertlhe(pho_) .and. usea(2)) then
         q4_q4 = 0d0
      else
         q4_q4 = scr(momext(:,4),momext(:,4))
      endif
      q5_q5 = scr(momext(:,5),momext(:,5))
      prop3 = propagator(dsqrt(q5_q5),mass(5,1),mass(5,2))
 
      if (includevprime) then
         if(.not.usea(1)) then
            vpffcoupl(1,1)=getvpffcoupling_vh(id(1), -1, ((id(1)+id(2)).ne.0))
            vpffcoupl(1,2)=getvpffcoupling_vh(id(1), +1, ((id(1)+id(2)).ne.0))
         endif
         if(.not.usea(2)) then
            vpffcoupl(2,1)=getvpffcoupling_vh(id(6), -1, ((id(6)+id(7)).ne.0))
            vpffcoupl(2,2)=getvpffcoupling_vh(id(6), +1, ((id(6)+id(7)).ne.0))
         endif
      endif
 
 
      if(.not.usea(1)) then
         prop1 = propagator(dsqrt(q3_q3),mass(3,1),mass(3,2))
         if(id(1).gt.0)then
            call ffv(id(2), momext(:,2), helicity(2), id(1), momext(:,1), helicity(1), vcurrent1)
            call ffa(id(2), momext(:,2), helicity(2), id(1), momext(:,1), helicity(1), acurrent1)
         else
            call ffv(id(1), momext(:,1), helicity(1), id(2), momext(:,2), helicity(2), vcurrent1)
            call ffa(id(1), momext(:,1), helicity(1), id(2), momext(:,2), helicity(2), acurrent1)
         endif
 
         ! Vpff current without the prefactor
         if (includevprime) then
            if((id(1)*helicity(1)).le.0d0)then
               currentvp1=( &
                  vcurrent1*(vpffcoupl(1,1)+vpffcoupl(1,2))*0.5 - &
                  acurrent1*(vpffcoupl(1,1)-vpffcoupl(1,2))*0.5   &
                  )
            else
               currentvp1=( &
                  vcurrent1*vpffcoupl(1,2)                        &
                  )
            endif
         endif
 
         !WH
         if((id(1)+id(2)).ne.0)then
            if (includevprime) then
               if (getmass(wppr_).ge.0d0) then
                  !print *,"Compute prop for Wppr"
                  prop_vp1 = propagator(dsqrt(q3_q3),getmass(wppr_),getdecaywidth(wppr_))
               else
                  prop_vp1 = propagator(m_w,0d0,0d0)
               endif
               currentvp1 = currentvp1*gffw*ckmbare(id(1),id(2))
            endif
            if((id(1)*helicity(1)).le.0d0)then
               current1=(vcurrent1-acurrent1)/2d0*gffw*ckmbare(id(1),id(2))
            else
               current1=0d0
            endif
         !ZH
         else
            if (includevprime) then
               if (getmass(zpr_).ge.0d0) then
                  !print *,"Compute prop for Zpr"
                  prop_vp1 = propagator(dsqrt(q3_q3),getmass(zpr_),getdecaywidth(zpr_))
               else
                  prop_vp1 = propagator(m_z,0d0,0d0)
               endif
               currentvp1 = currentvp1*gffz
            endif
            !e+ e- Z vertex for incoming states
            if((abs(id(1)).eq.11).or.(abs(id(1)).eq.13).or.(abs(id(1)).eq.15))then
              if((id(1)*helicity(1)).gt.0d0)then
                current1=(0.5d0*t3lr - qlr*sitw**2) *vcurrent1 -(0.5d0*t3lr)*acurrent1
              else
                current1=(0.5d0*t3ll - qll*sitw**2) *vcurrent1 -(0.5d0*t3ll)*acurrent1
              endif
              current1=current1*gffz
            !u u~ Z vertex for incoming states
            else if((abs(id(1)).eq.2).or.(abs(id(1)).eq.4))then
              if((id(1)*helicity(1)).gt.0d0)then
                current1=(0.5d0*t3ur - qur*sitw**2) *vcurrent1 -(0.5d0*t3ur)*acurrent1
              else
                current1=(0.5d0*t3ul - qul*sitw**2) *vcurrent1 -(0.5d0*t3ul)*acurrent1
              endif
              current1=current1*gffz
            !d d~ Z vertex for incoming states
            else if((abs(id(1)).eq.1).or.(abs(id(1)).eq.3).or.(abs(id(1)).eq.5))then
              if((id(1)*helicity(1)).gt.0d0)then
                current1=(0.5d0*t3dr - qdr*sitw**2) *vcurrent1 -(0.5d0*t3dr)*acurrent1
              else
                current1=(0.5d0*t3dl - qdl*sitw**2) *vcurrent1 -(0.5d0*t3dl)*acurrent1
              endif
              current1=current1*gffz
            else
              current1=0d0
              currentvp1=0d0
              print *, "invalid incoming state"
            endif
         endif
      else
         if(abs(id(1)).eq.convertlhe(pho_)) then
           prop1=cone
           if((id(1)*helicity(1)).gt.0d0) then
             call polarization_single(momext(:,3),+1,vcurrent1)
           else
             call polarization_single(momext(:,3),-1,vcurrent1)
           endif
         else
           prop1 = propagator(dsqrt(q3_q3),0d0,0d0)
           if(id(1).gt.0)then
             call ffv(id(2), momext(:,2), helicity(2), id(1), momext(:,1), helicity(1), vcurrent1)
           else
             call ffv(id(1), momext(:,1), helicity(1), id(2), momext(:,2), helicity(2), vcurrent1)
           endif
         endif
 
         !ZH
         if(abs(id(1)).eq.convertlhe(pho_))then
           current1=vcurrent1
         !e+ e- Z vertex for incoming states
         else if((abs(id(1)).eq.11).or.(abs(id(1)).eq.13).or.(abs(id(1)).eq.15))then
           if((id(1)*helicity(1)).gt.0d0)then
             current1 = qlr*vcurrent1
           else
             current1 = qll*vcurrent1
           endif
           current1=current1*gffa
         !u u~ Z vertex for incoming states
         else if((abs(id(1)).eq.2).or.(abs(id(1)).eq.4))then
           if((id(1)*helicity(1)).gt.0d0)then
             current1 = qur*vcurrent1
           else
             current1 = qul*vcurrent1
           endif
           current1=current1*gffa
         !d d~ Z vertex for incoming states
         else if((abs(id(1)).eq.1).or.(abs(id(1)).eq.3).or.(abs(id(1)).eq.5))then
           if((id(1)*helicity(1)).gt.0d0)then
             current1 = qdr*vcurrent1
           else
             current1 = qdl*vcurrent1
           endif
           current1=current1*gffa
         else
           current1=0d0
           print *, "invalid incoming state"
         endif
      endif
 
      if(.not.usea(2)) then
         prop2 = propagator(dsqrt(q4_q4),mass(4,1),mass(4,2))
 
         if(id(6).gt.0)then
           call ffv(id(6), momext(:,6), helicity(6), id(7), momext(:,7), helicity(7), vcurrent2)
           call ffa(id(6), momext(:,6), helicity(6), id(7), momext(:,7), helicity(7), acurrent2)
         else
           call ffv(id(7), momext(:,7), helicity(7), id(6), momext(:,6), helicity(6), vcurrent2)
           call ffa(id(7), momext(:,7), helicity(7), id(6), momext(:,6), helicity(6), acurrent2)
         endif
 
         ! Vpff current without the prefactor
         if (includevprime) then
            if((id(6)*helicity(6)).le.0d0)then
               currentvp2=( &
                  vcurrent2*(vpffcoupl(2,1)+vpffcoupl(2,2))*0.5 - &
                  acurrent2*(vpffcoupl(2,1)-vpffcoupl(2,2))*0.5   &
                  )
            else
               currentvp2=( &
                  vcurrent2*vpffcoupl(2,2)                        &
                  )
            endif
         endif
 
         !WH
         if((id(6)+id(7)).ne.0)then
           if (includevprime) then
             if (getmass(wppr_).ge.0d0) then
               prop_vp2 = propagator(dsqrt(q4_q4),getmass(wppr_),getdecaywidth(wppr_))
             else
               prop_vp2 = propagator(m_w,0d0,0d0)
             endif
             currentvp2 = currentvp2*gffw*ckmbare(id(6),id(7))
           endif
           if((id(6)*helicity(6)).le.0d0)then
             current2=(vcurrent2-acurrent2)/2d0*gffw*ckm(id(6),id(7))
           else
             current2=0d0
           endif
         !ZH
         else
           if (includevprime) then
             if (getmass(zpr_).ge.0d0) then
               prop_vp2 = propagator(dsqrt(q4_q4),getmass(zpr_),getdecaywidth(zpr_))
             else
               prop_vp2 = propagator(m_z,0d0,0d0)
             endif
             currentvp2 = currentvp2*gffz
           endif
           !l+ l- Z vertex for final state
           if((abs(id(6)).eq.11).or.(abs(id(6)).eq.13))then
             if((id(6)*helicity(6)).gt.0d0)then
               current2=(0.5d0*t3lr - qlr*sitw**2) *vcurrent2 -(0.5d0*t3lr)*acurrent2
             else
               current2=(0.5d0*t3ll - qll*sitw**2) *vcurrent2 -(0.5d0*t3ll)*acurrent2
             endif
             current2=current2*gffz*dsqrt(scale_alpha_z_ll)
           !tau+ tau- Z vertex for final state
           else if((abs(id(6)).eq.15))then
             if((id(6)*helicity(6)).gt.0d0)then
               current2=(0.5d0*t3lr - qlr*sitw**2) *vcurrent2 -(0.5d0*t3lr)*acurrent2
             else
               current2=(0.5d0*t3ll - qll*sitw**2) *vcurrent2 -(0.5d0*t3ll)*acurrent2
             endif
             current2=current2*gffz*dsqrt(scale_alpha_z_tt)
           !u u~ Z vertex for final state
           else if((abs(id(6)).eq.2).or.(abs(id(6)).eq.4))then
             if((id(6)*helicity(6)).gt.0d0)then
               current2=(0.5d0*t3ur - qur*sitw**2) *vcurrent2 -(0.5d0*t3ur)*acurrent2
             else
               current2=(0.5d0*t3ul - qul*sitw**2) *vcurrent2 -(0.5d0*t3ul)*acurrent2
             endif
             current2=current2*gffz*dsqrt(scale_alpha_z_uu)
           !d d~ Z vertex for final state
           else if((abs(id(6)).eq.1).or.(abs(id(6)).eq.3).or.(abs(id(6)).eq.5))then
             if((id(6)*helicity(6)).gt.0d0)then
               current2=(0.5d0*t3dr - qdr*sitw**2) *vcurrent2 -(0.5d0*t3dr)*acurrent2
             else
               current2=(0.5d0*t3dl - qdl*sitw**2) *vcurrent2 -(0.5d0*t3dl)*acurrent2
             endif
             current2=current2*gffz*dsqrt(scale_alpha_z_dd)
           !nu nu~ Z vertex for final state
           else if((abs(id(6)).eq.12).or.(abs(id(6)).eq.14).or.(abs(id(6)).eq.16))then
             current2=(0.5d0*t3nl - qnl*sitw**2) *vcurrent2 -(0.5d0*t3nl)*acurrent2
             current2=current2*gffz*dsqrt(scale_alpha_z_nn)
           else
             current2=0d0
             currentvp2 = 0d0
             print *, "invalid final state", id(6:7), helicity(6:7)
             stop
           endif
         endif
      else
         if(abs(id(6)).eq.convertlhe(pho_)) then
           prop2=cone
           if((id(6)*helicity(6)).gt.0d0) then
             call polarization_single(momext(:,4),+1,vcurrent2)
           else
             call polarization_single(momext(:,4),-1,vcurrent2)
           endif
           vcurrent2 = dconjg(vcurrent2)
         else
           prop2 = propagator(dsqrt(q4_q4),0d0,0d0)
           if(id(6).gt.0)then
             call ffv(id(6), momext(:,6), helicity(6), id(7), momext(:,7), helicity(7), vcurrent2)
           else
             call ffv(id(7), momext(:,7), helicity(7), id(6), momext(:,6), helicity(6), vcurrent2)
           endif
         endif
 
         !ZH
         if(abs(id(6)).eq.convertlhe(pho_)) then
           current2=vcurrent2
         !l+ l- Z vertex for final state
         else if((abs(id(6)).eq.11).or.(abs(id(6)).eq.13))then
           if((id(6)*helicity(6)).gt.0d0)then
             current2=qlr*vcurrent2
           else
             current2=qll*vcurrent2
           endif
           current2=current2*gffa*dsqrt(scale_alpha_z_ll)
         !tau+ tau- Z vertex for final state
         else if((abs(id(6)).eq.15))then
           if((id(6)*helicity(6)).gt.0d0)then
             current2=qlr*vcurrent2
           else
             current2=qll*vcurrent2
           endif
           current2=current2*gffa*dsqrt(scale_alpha_z_tt)
         !u u~ Z vertex for final state
         else if((abs(id(6)).eq.2).or.(abs(id(6)).eq.4))then
           if((id(6)*helicity(6)).gt.0d0)then
             current2=qur*vcurrent2
           else
             current2=qul*vcurrent2
           endif
           current2=current2*gffa*dsqrt(scale_alpha_z_uu)
         !d d~ Z vertex for final state
         else if((abs(id(6)).eq.1).or.(abs(id(6)).eq.3).or.(abs(id(6)).eq.5))then
           if((id(6)*helicity(6)).gt.0d0)then
             current2=qdr*vcurrent2
           else
             current2=qdl*vcurrent2
           endif
           current2=current2*gffa*dsqrt(scale_alpha_z_dd)
         !nu nu~ Z vertex for final state
         else if((abs(id(6)).eq.12).or.(abs(id(6)).eq.14).or.(abs(id(6)).eq.16))then
           current2=qnl*vcurrent2
           current2=current2*gffa*dsqrt(scale_alpha_z_nn)
         else
           current2=0d0
           print *, "invalid final state", id(6:7), helicity(6:7)
           stop
         endif
      endif
 
      if(.not.(usea(1) .and. abs(id(1)).eq.convertlhe(pho_))) then
         current1 = -current1 + scrc(momext(:,3),current1)/q3_q3
         currentvp1 = -currentvp1 + scrc(momext(:,3),currentvp1)/q3_q3
      endif
      if(.not.(usea(2) .and. abs(id(6)).eq.convertlhe(pho_))) then
         current2 = -current2 + scrc(momext(:,4),current2)/q4_q4
         currentvp2 = -currentvp2 + scrc(momext(:,4),currentvp2)/q4_q4
      endif
 
      !print *,"current1=",current1
      !print *,"currentVp1=",currentVp1
      !print *,"PROP1=",PROP1
      !print *,"PROP_Vp1=",PROP_Vp1
      !print *,"current2=",current2
      !print *,"currentVp2=",currentVp2
      !print *,"PROP2=",PROP2
      !print *,"PROP_Vp2=",PROP_Vp2
 
      current1 = current1*prop1
      current2 = current2*prop2
      currentvp1 = currentvp1*prop_vp1
      currentvp2 = currentvp2*prop_vp2
 
!XVV vertex
      if(id(3).eq.convertlhe(wp_))then
         call swap(q3_q3,q4_q4)
         call swap(current1,current2)
         call swap(currentvp1,currentvp2)
      endif
 
      if(.not.usea(1) .and. .not.usea(2)) then
         ghz1_dyn = hvvspinzerodynamiccoupling(1,q3_q3,q4_q4,q5_q5)
         ghz2_dyn = hvvspinzerodynamiccoupling(2,q3_q3,q4_q4,q5_q5)
         ghz3_dyn = hvvspinzerodynamiccoupling(3,q3_q3,q4_q4,q5_q5)
         ghz4_dyn = hvvspinzerodynamiccoupling(4,q3_q3,q4_q4,q5_q5)
 
         if (includevprime) then
            ghzzp1_dyn = hvvspinzerodynamiccoupling(12,q3_q3,q4_q4,q5_q5)
            ghzzp2_dyn = hvvspinzerodynamiccoupling(13,q3_q3,q4_q4,q5_q5)
            ghzzp3_dyn = hvvspinzerodynamiccoupling(14,q3_q3,q4_q4,q5_q5)
            ghzzp4_dyn = hvvspinzerodynamiccoupling(15,q3_q3,q4_q4,q5_q5)
 
            ghzpz1_dyn = hvvspinzerodynamiccoupling(12,q4_q4,q3_q3,q5_q5)
            ghzpz2_dyn = hvvspinzerodynamiccoupling(13,q4_q4,q3_q3,q5_q5)
            ghzpz3_dyn = hvvspinzerodynamiccoupling(14,q4_q4,q3_q3,q5_q5)
            ghzpz4_dyn = hvvspinzerodynamiccoupling(15,q4_q4,q3_q3,q5_q5)
 
            ghzpzp1_dyn = hvvspinzerodynamiccoupling(16,q3_q3,q4_q4,q5_q5)
            ghzpzp2_dyn = hvvspinzerodynamiccoupling(17,q3_q3,q4_q4,q5_q5)
            ghzpzp3_dyn = hvvspinzerodynamiccoupling(18,q3_q3,q4_q4,q5_q5)
            ghzpzp4_dyn = hvvspinzerodynamiccoupling(19,q3_q3,q4_q4,q5_q5)
         endif
      else if(usea(1) .and. usea(2)) then
         ghz1_dyn = czero
         ghz2_dyn = hvvspinzerodynamiccoupling(9,q3_q3,q4_q4,q5_q5)
         ghz3_dyn = hvvspinzerodynamiccoupling(10,q3_q3,q4_q4,q5_q5)
         ghz4_dyn = hvvspinzerodynamiccoupling(11,q3_q3,q4_q4,q5_q5)
      else if(usea(1)) then
         ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q3_q3,q5_q5)
         ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q3_q3,q5_q5)
         ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q3_q3,q5_q5)
         ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q3_q3,q5_q5)
 
         if (includevprime) then
            ghzzp1_dyn = hvvspinzerodynamiccoupling(20,0d0,q3_q3,q5_q5)
            ghzzp2_dyn = hvvspinzerodynamiccoupling(21,0d0,q3_q3,q5_q5)
            ghzzp3_dyn = hvvspinzerodynamiccoupling(22,0d0,q3_q3,q5_q5)
            ghzzp4_dyn = hvvspinzerodynamiccoupling(23,0d0,q3_q3,q5_q5)
         endif
      else !if(useA(2)) then
         ghz1_dyn = hvvspinzerodynamiccoupling(5,0d0,q4_q4,q5_q5)
         ghz2_dyn = hvvspinzerodynamiccoupling(6,0d0,q4_q4,q5_q5)
         ghz3_dyn = hvvspinzerodynamiccoupling(7,0d0,q4_q4,q5_q5)
         ghz4_dyn = hvvspinzerodynamiccoupling(8,0d0,q4_q4,q5_q5)
 
         if (includevprime) then
            ghzpz1_dyn = hvvspinzerodynamiccoupling(20,0d0,q4_q4,q5_q5)
            ghzpz2_dyn = hvvspinzerodynamiccoupling(21,0d0,q4_q4,q5_q5)
            ghzpz3_dyn = hvvspinzerodynamiccoupling(22,0d0,q4_q4,q5_q5)
            ghzpz4_dyn = hvvspinzerodynamiccoupling(23,0d0,q4_q4,q5_q5)
         endif
      endif
 
      gvvs1 = ghz1_dyn*(mass(3,1)**2) + qq * ( 2d0*ghz2_dyn + ghz3_dyn*qq/lambda**2 )
      gvvs2 = -( 2d0*ghz2_dyn + ghz3_dyn*qq/lambda**2 )
      gvvp = -2d0*ghz4_dyn
 
      if (includevprime) then
         if(.not.usea(1) .and. .not.usea(2)) then
            gvvps1 = ghzzp1_dyn*(mass(3,1)**2) + qq * ( 2d0*ghzzp2_dyn + ghzzp3_dyn*qq/lambda**2 )
            gvvps2 = -( 2d0*ghzzp2_dyn + ghzzp3_dyn*qq/lambda**2 )
            gvvpp = -2d0*ghzzp4_dyn
 
            gvpvs1 = ghzpz1_dyn*(mass(3,1)**2) + qq * ( 2d0*ghzpz2_dyn + ghzpz3_dyn*qq/lambda**2 )
            gvpvs2 = -( 2d0*ghzpz2_dyn + ghzpz3_dyn*qq/lambda**2 )
            gvpvp = -2d0*ghzpz4_dyn
 
            gvpvps1 = ghzpzp1_dyn*(mass(3,1)**2) + qq * ( 2d0*ghzpzp2_dyn + ghzpzp3_dyn*qq/lambda**2 )
            gvpvps2 = -( 2d0*ghzpzp2_dyn + ghzpzp3_dyn*qq/lambda**2 )
            gvpvpp = -2d0*ghzpzp4_dyn
         else if (usea(1)) then
            gvvps1 = ghzzp1_dyn*(mass(3,1)**2) + qq * ( 2d0*ghzzp2_dyn + ghzzp3_dyn*qq/lambda**2 )
            gvvps2 = -( 2d0*ghzzp2_dyn + ghzzp3_dyn*qq/lambda**2 )
            gvvpp = -2d0*ghzzp4_dyn
         else! if (useA(2)) then
            gvpvs1 = ghzpz1_dyn*(mass(3,1)**2) + qq * ( 2d0*ghzpz2_dyn + ghzpz3_dyn*qq/lambda**2 )
            gvpvs2 = -( 2d0*ghzpz2_dyn + ghzpz3_dyn*qq/lambda**2 )
            gvpvp = -2d0*ghzpz4_dyn
         endif
      endif
 
 
      call vvs1(g_mu_nu)
      call vvs2(momext(:,5),momext(:,5),pp)
      if(id(3).eq.convertlhe(wp_))then
         call vvp(momext(:,4),-momext(:,3),epp)
      else
         call vvp(-momext(:,3),momext(:,4),epp)
      endif
 
! assemble everything and get iM
      matrixelement0=(0d0,0d0)
      do mu3=1,4
      do mu4=1,4
         matrixelement0 = matrixelement0 +      &
            current1(mu3)*current2(mu4)*(       &
               gvvs1*g_mu_nu(mu3,mu4)         + &
               gvvs2*pp(mu3,mu4)         + &
               gvvp *epp(mu3,mu4)           &
            )                                   &
            +                                   &
            current1(mu3)*currentvp2(mu4)*(    &
               gvvps1*g_mu_nu(mu3,mu4)        + &
               gvvps2*pp(mu3,mu4)        + &
               gvvpp *epp(mu3,mu4)          &
            )                                   &
            +                                   &
            currentvp1(mu3)*current2(mu4)*(    &
               gvpvs1*g_mu_nu(mu3,mu4)        + &
               gvpvs2*pp(mu3,mu4)        + &
               gvpvp *epp(mu3,mu4)          &
            )                                   &
            +                                   &
            currentvp1(mu3)*currentvp2(mu4)*( &
               gvpvps1*g_mu_nu(mu3,mu4)       + &
               gvpvps2*pp(mu3,mu4)       + &
               gvpvpp *epp(mu3,mu4)         &
            )
      enddo !mu4
      enddo !mu3
      matrixelement0 = matrixelement0*ci/vev
 
      if(h_dk.eqv..false.)then
        matrixelement0=matrixelement0 *prop3
      else if(id(8).ne.not_a_particle_) then
        matrixelement0=matrixelement0 *prop3 &
        *(kappa*ffs(id(8), momext(:,8), helicity(8), id(9), momext(:,9), helicity(9)) &
         +kappa_tilde*ffp(id(8), momext(:,8), helicity(8), id(9), momext(:,9), helicity(9)))&
        *(-ci/vev*getmass(convertlhereverse(id(8))))
      else
        matrixelement0=czero
      endif
 
      !print *,"MATRIXELEMENT0=",MATRIXELEMENT0
 
      return
      END function
 
 
 
function getvpffcoupling_vh(pdgid, hel, useWp)
integer, intent(in) :: pdgid
integer, intent(in) :: hel
logical, intent(in) :: usewp
complex(8) :: getvpffcoupling_vh
   getvpffcoupling_vh=vpffcoupling_pdg(pdgid,hel,usewp)
   if(usewp) then
     getvpffcoupling_vh = getvpffcoupling_vh / bl ! Bc of the couplings formalism from decay
   else
     getvpffcoupling_vh = getvpffcoupling_vh / 2.0_dp
   endif
end function
 
 
!ANGLES.F
!VERSION 20130531
!
!A subroutine that calculates the polar and azimuthal angles of a
!given vector(4) in terms of their sin and cos, which will be
!returned by the array sincos(4).
      subroutine angles(sincos, vector)
      implicit none
!     real(8) Pi
      real(8) Twopi, Fourpi, epsilon
!     parameter( Pi = 3.14159265358979323846d0 )
      parameter( twopi = 2d0 * pi )
      parameter( fourpi = 4d0 * pi )
      parameter( epsilon = 1d-13 ) !a small quantity slightly above machine precision
      real(8) sincos(4), vector(4), abs3p, phi
!sincos(1)=cos(theta)
!sincos(2)=sin(theta)
!sincos(3)=cos(phi)
!sincos(4)=sin(phi)
 
!|3-momentum|
      abs3p = dsqrt(vector(2)**2+vector(3)**2+vector(4)**2)
 
!if |3-momentum|=0
      if(abs3p.lt.epsilon)then
        sincos(1)=1d0
        sincos(2)=0d0
      else
        sincos(1)=vector(4)/abs3p
        sincos(2)=dsqrt((1d0+sincos(1))*(1d0-sincos(1)))
      endif
 
!if colinear
      if(dabs(vector(3)).lt.epsilon)then
        phi=0d0
      else
        if(dabs(vector(2)).lt.epsilon)then
           phi=(twopi/2d0)/2d0 * dsign(1d0,vector(3))
        else
           phi=datan(vector(3)/vector(2))
        endif
      endif
!shift phi so that 0 < phi < 2Pi
      if(vector(2).lt.0d0)then
         phi=phi+pi
      endif
      if(phi.lt.0d0)then
         phi=phi+twopi
      endif
!     print *,phi
      sincos(3)=dcos(phi)
      sincos(4)=dsin(phi)
 
      return
      END subroutine angles
 
!ANTISYMMETRIC2.F
!VERSION 20130702
      subroutine antisymmetric2(p1,p2,epp)
 
      implicit none
!     real(8) Pi
      real(8) Twopi, Fourpi, epsilon
!     parameter( Pi = 3.14159265358979323846d0 )
      parameter( twopi = 2d0 * pi )
      parameter( fourpi = 4d0 * pi )
      parameter( epsilon = 1d-13 ) !a small quantity slightly above machine precision
 
      complex(8) p1(4), p2(4)
      complex(8) epp(4,4)
!      real(8) ANTISYMMETRIC
      integer i,j,k,l
 
!      external ANTISYMMETRIC
 
!     do i=1,4
!     do j=1,4
      epp(i,j)=0d0
!     enddo
!     enddo
 
      do i=1,4
      do j=1,4
      do k=1,4
      do l=1,4
          epp(i,j)=epp(i,j)+antisymmetric(i,j,k,l)*p1(k)*p2(l)
      enddo
      enddo
      enddo
      enddo
 
      return
      END subroutine antisymmetric2
 
!ANTISYMMETRIC.F
!VERSION 20130618
 
!returns the element of the rank-4 COVARIANT total antysymmetric
!tensor.
!ANTISYMMETRIC(0,1,2,3)=1.
      real(8) function ANTISYMMETRIC(i,j,k,l)
 
      implicit none
!     include '../COMMON.INI'
 
      integer i,j,k,l
 
      antisymmetric=dble((i-j)*(i-k)*(i-l)*(j-k)*(j-l)*(k-l))/12d0
 
      return
      END function antisymmetric
 
!CONTRA_FIELD_TENSOR.F
!VERSION 20130630
 
!in T_mu_nu(lambda, ^mu, ^nu) returns the contrvariant field tensors
!for given polarization vectors.
      subroutine contra_field_tensor(POL, T_mu_nu)
 
      implicit none
!     real(8) Pi
      real(8) Twopi, Fourpi, epsilon
!     parameter( Pi = 3.14159265358979323846d0 )
      parameter( twopi = 2d0 * pi )
      parameter( fourpi = 4d0 * pi )
      parameter( epsilon = 1d-13 ) !a small quantity slightly above machine precision
      complex(8) epep(4,4),emem(4,4),epe0(4,4),eme0(4,4),e0e0(4,4)
      complex(8) epem(4,4),e0ep(4,4),e0em(4,4),emep(4,4)
      complex(8) POL(3,4), T_mu_nu(5,4,4)
 
      call contra_outer(pol(1,:), pol(1,:), epep)
      call contra_outer(pol(2,:), pol(2,:), emem)
      call contra_outer(pol(1,:), pol(3,:), epe0)
      call contra_outer(pol(3,:), pol(1,:), e0ep)
      call contra_outer(pol(2,:), pol(3,:), eme0)
      call contra_outer(pol(3,:), pol(2,:), e0em)
      call contra_outer(pol(3,:), pol(3,:), e0e0)
      call contra_outer(pol(1,:), pol(2,:), epem)
      call contra_outer(pol(2,:), pol(1,:), emep)
 
!lambda = +2
      t_mu_nu(1,:,:)=epep
!lambda = -2
      t_mu_nu(2,:,:)=emem
!lambda = +3
      t_mu_nu(3,:,:)=(epe0+e0ep)/dsqrt(2d0)
!lambda = -3
      t_mu_nu(4,:,:)=(eme0+e0em)/dsqrt(2d0)
!lambda = 0
      t_mu_nu(5,:,:)=(epem+emep)/dsqrt(6d0) + e0e0/dsqrt(1.5d0)
 
 
      return
      END subroutine contra_field_tensor
 
!CONTRA_OUTER.F
!VERSION 20130620
 
!in pp(_mu, _nu) returns the CONTRAVARIANT tensor of
!p1 p2 outer product.
      subroutine contra_outer(p1,p2,pp)
 
      implicit none
!     include '../COMMON.INI'
      complex(8) p1(4), p2(4)
      complex(8) pp(4,4)
      integer mu, nu
 
      do mu=1,4
        do nu=1,4
          pp(mu,nu)=p1(mu)*p2(nu)
        enddo
      enddo
 
      return
      END subroutine contra_outer
 
!COVARIANT_FIELD_TENSOR.F
!VERSION 20130630
 
!in T_mu_nu(lambda, _mu, _nu) returns the covariant field tensors
!for given polarization vectors.
      subroutine covariant_field_tensor(POL, T_mu_nu)
 
      implicit none
!     real(8) Pi
      real(8) Twopi, Fourpi, epsilon
!     parameter( Pi = 3.14159265358979323846d0 )
      parameter( twopi = 2d0 * pi )
      parameter( fourpi = 4d0 * pi )
      parameter( epsilon = 1d-13 ) !a small quantity slightly above machine precision
      complex(8) epep(4,4),emem(4,4),epe0(4,4),eme0(4,4),e0e0(4,4)
      complex(8) epem(4,4),e0ep(4,4),e0em(4,4),emep(4,4)
      complex(8) POL(3,4), T_mu_nu(5,4,4)
 
      call covariant_outer(pol(1,:), pol(1,:), epep)
      call covariant_outer(pol(2,:), pol(2,:), emem)
      call covariant_outer(pol(1,:), pol(3,:), epe0)
      call covariant_outer(pol(3,:), pol(1,:), e0ep)
      call covariant_outer(pol(2,:), pol(3,:), eme0)
      call covariant_outer(pol(3,:), pol(2,:), e0em)
      call covariant_outer(pol(3,:), pol(3,:), e0e0)
      call covariant_outer(pol(1,:), pol(2,:), epem)
      call covariant_outer(pol(2,:), pol(1,:), emep)
 
!lambda = +2
      t_mu_nu(1,:,:)=epep
!lambda = -2
      t_mu_nu(2,:,:)=emem
!lambda = +3
      t_mu_nu(3,:,:)=(epe0+e0ep)/dsqrt(2d0)
!lambda = -3
      t_mu_nu(4,:,:)=(eme0+e0em)/dsqrt(2d0)
!lambda = 0
      t_mu_nu(5,:,:)=(epem+emep)/dsqrt(6d0) + e0e0/dsqrt(1.5d0)
 
      return
      END subroutine covariant_field_tensor
 
!COVARIANT_OUTER.F
!VERSION 20130620
 
!in pp(_mu, _nu) returns the COVARIANT tensor of p1 p2
!outer product.
      subroutine covariant_outer(p1,p2,pp)
 
      implicit none
!     include '../COMMON.INI'
      complex(8) p1(4), p2(4)
      complex(8) pp(4,4)
      integer mu, nu
 
      do mu=1,4
        do nu=1,4
          pp(mu,nu)=p1(mu)*p2(nu)
          if( ( (mu.ne.1).and.(nu.eq.1) ).or. &
              ( (mu.eq.1).and.(nu.ne.1) ) )then
            pp(mu,nu)=-pp(mu,nu)
          endif
        enddo
      enddo
 
      return
      END subroutine covariant_outer
 
!COVARIANT_VECTOR.F
!VERSION 20130703
 
!returns the component of the COVARIANT vector for given 4-vector
!and Lorentz index.
      complex(8) function covariant_vector(p,mu)
 
      implicit none
 
      complex(8) p(4)
      integer mu
 
      if(mu.ne.1)then
        covariant_vector = -p(mu)
      else
        covariant_vector = p(mu)
      endif
 
      return
      END function covariant_vector
 
!FFP.A
!VERSION 20130522
 
!returns i.Psi~(p1,s1).gamma5.Psi(p2,s2) for massless states
      complex(8) function ffp(pdg_code1, p1, h1, pdg_code2, p2, h2)
 
      implicit none
      real(8), parameter :: epsilon = 1d-13 !a small quantity slightly above machine precision
 
      real(8) p1(4), p2(4), h1, h2
      integer pdg_code1, pdg_code2
      real(8) sqrt_pp1dpp2
 
      if( ( dble(pdg_code1) *h1* dble(pdg_code2) *h2 ).gt.0d0)then
        ffp=0d0
 
      else if( ( dabs( p1(1)+p1(4) ).lt.epsilon ).and. (   dabs( p2(1)+p2(4) ).lt.epsilon ) )then
        ffp=0d0
      else if(   dabs( p1(1)+p1(4) ).lt.epsilon )then
        ffp=-dsqrt(2d0*p1(1)*(p2(1)+p2(4)))
      else if(   dabs( p2(1)+p2(4) ).lt.epsilon )then
        ffp= dsqrt(2d0*p2(1)*(p1(1)+p1(4)))
      else
        sqrt_pp1dpp2 = dsqrt((p1(1)+p1(4))/(p2(1)+p2(4)))
        ffp=(p2(2)-(0d0,1d0)*p2(3))*sqrt_pp1dpp2- (p1(2)-(0d0,1d0)*p1(3))/sqrt_pp1dpp2
      endif
 
      ffp=ffp*(0d0,-1d0)
 
      if( (dble(pdg_code1)*h1) .lt. 0d0)then
        ffp=-dconjg(ffp)
 
      endif
 
      return
      END function ffp
 
!FFA.F
!VERSION 20130523
 
!in Acurrent(4) returns Psi~(p1,s1).gamma^mu.gamma5.Psi(p2,s2) for massless
!states.
      subroutine ffa(pdg_code1, p1, h1, pdg_code2, p2, h2, Acurrent)
 
      implicit none
      real(8), parameter :: epsilon = 1d-13 !a small quantity slightly above machine precision
 
      real(8) p1(4), p2(4), h1, h2
      integer pdg_code1, pdg_code2
      real(8) sqrt_pp1Dpp2, sqrt_pp1Xpp2
      complex(8) Acurrent(4)
      integer mu
 
      acurrent = (0d0,0d0)
 
      if( ( dble(pdg_code1) *h1* dble(pdg_code2) *h2 ).lt.0d0)then
        do mu=1,4
          acurrent(mu)=0d0
        enddo
 
      else if( ( dabs( p1(1)+p1(4) ).lt.epsilon ).and. &
               ( dabs( p2(1)+p2(4) ).lt.epsilon ) )then
 
        acurrent(1)= 2d0*dsqrt(p1(1)*p2(1))
        acurrent(2)= 0d0
        acurrent(3)= 0d0
        acurrent(4)=-acurrent(1)
 
      else if(   dabs( p1(1)+p1(4) ).lt.epsilon )then
 
        acurrent(1)= dsqrt( 2d0*p1(1) / ( p2(1)+p2(4) ) ) &
                    *( p2(2) + (0d0,1d0)*p2(3) )
        acurrent(2)= dsqrt( 2d0*p1(1) * ( p2(1)+p2(4) ) )
        acurrent(3)= (0d0,1d0)*acurrent(2)
        acurrent(4)=-acurrent(1)
 
      else if(   dabs( p2(1)+p2(4) ).lt.epsilon )then
 
        acurrent(1)= dsqrt( 2d0*p2(1) / ( p1(1)+p1(4) ) ) &
                    *( p1(2) - (0d0,1d0)*p1(3) )
        acurrent(2)= dsqrt( 2d0*p2(1) * ( p1(1)+p1(4) ) )
        acurrent(3)=-(0d0,1d0)*acurrent(2)
        acurrent(4)=-acurrent(1)
 
      else
 
        sqrt_pp1dpp2= dsqrt( (p1(1)+p1(4)) / (p2(1)+p2(4)) )
        sqrt_pp1xpp2= dsqrt( (p1(1)+p1(4)) * (p2(1)+p2(4)) )
        acurrent(1)= sqrt_pp1xpp2         &
                    +( p1(2) - (0d0,1d0)*p1(3) )     &
                    *( p2(2) + (0d0,1d0)*p2(3) )/sqrt_pp1xpp2
        acurrent(2)= ( p1(2) - (0d0,1d0)*p1(3) )/sqrt_pp1dpp2   &
                    +( p2(2) + (0d0,1d0)*p2(3) )*sqrt_pp1dpp2
        acurrent(3)= ( (0d0,1d0)*p1(2) + p1(3) )/sqrt_pp1dpp2   &
                    -( (0d0,1d0)*p2(2) - p2(3) )*sqrt_pp1dpp2
        acurrent(4)=sqrt_pp1xpp2          &
                    -( p1(2) - (0d0,1d0)*p1(3) )                &
                    *( p2(2) + (0d0,1d0)*p2(3) )/sqrt_pp1xpp2
      endif
 
!     print *, pdg_code1,h1,dble(pdg_code1)*h1,'!'
      if( (dble(pdg_code1)*h1) .lt. 0d0)then
!       print *, Acurrent
        do mu=1,4
          acurrent(mu)=-dconjg(acurrent(mu))
        enddo
      endif
 
      return
      END subroutine ffa
 
!FFS.F
!VERSION 20130522
 
!returns Psi~(p1,s1).Psi(p2,s2) for massless states
      complex(8) function ffs(pdg_code1, p1, h1, pdg_code2, p2, h2)
 
      implicit none
      real(8), parameter :: epsilon = 1d-13 !a small quantity slightly above machine precision
 
      real(8) p1(4), p2(4), h1, h2
      integer pdg_code1, pdg_code2
      real(8) sqrt_pp1dpp2
 
      ffs = (0d0,0d0)
 
      if( ( dble(pdg_code1) *h1* dble(pdg_code2) *h2 ).gt.0d0)then
        ffs=0d0
 
      else if( ( dabs( p1(1)+p1(4) ).lt.epsilon ).and. (   dabs( p2(1)+p2(4) ).lt.epsilon ) )then
        ffs=0d0
      else if(   dabs( p1(1)+p1(4) ).lt.epsilon )then
        ffs=-dsqrt(2d0*p1(1)*(p2(1)+p2(4)))
      else if(   dabs( p2(1)+p2(4) ).lt.epsilon )then
        ffs= dsqrt(2d0*p2(1)*(p1(1)+p1(4)))
      else
        sqrt_pp1dpp2 = dsqrt((p1(1)+p1(4))/(p2(1)+p2(4)))
        ffs=(p2(2)-(0d0,1d0)*p2(3))*sqrt_pp1dpp2- (p1(2)-(0d0,1d0)*p1(3))/sqrt_pp1dpp2
      endif
 
      if( (dble(pdg_code1)*h1) .lt. 0d0)then
        ffs=-dconjg(ffs)
 
      endif
 
      return
      END function ffs
 
!FFV.F
!VERSION 20130523
 
!in Vcurrent(4) returns Psi~(p1,s1).gamma^mu.Psi(p2,s2) for massless
!states.
      subroutine ffv(pdg_code1, p1, h1, pdg_code2, p2, h2, Vcurrent)
 
      implicit none
      real(8), parameter :: epsilon = 1d-13 !a small quantity slightly above machine precision
 
      real(8) p1(4), p2(4), h1, h2
      integer pdg_code1, pdg_code2
      real(8) sqrt_pp1Dpp2, sqrt_pp1Xpp2
      complex(8) Vcurrent(4)
      integer mu
 
      vcurrent = (0d0,0d0)
 
      if( ( dble(pdg_code1) *h1* dble(pdg_code2) *h2 ).lt.0d0)then
        do mu=1,4
          vcurrent(mu)=0d0
        enddo
 
      else if( ( dabs( p1(1)+p1(4) ).lt.epsilon ).and. ( dabs( p2(1)+p2(4) ).lt.epsilon ) )then
 
        vcurrent(1)= 2d0*dsqrt(p1(1)*p2(1))
        vcurrent(2)= 0d0
        vcurrent(3)= 0d0
        vcurrent(4)=-vcurrent(1)
 
      else if(   dabs( p1(1)+p1(4) ).lt.epsilon )then
 
        vcurrent(1)= dsqrt( 2d0*p1(1) / ( p2(1)+p2(4) ) ) *( p2(2) + (0d0,1d0)*p2(3) )
        vcurrent(2)= dsqrt( 2d0*p1(1) * ( p2(1)+p2(4) ) )
        vcurrent(3)= (0d0,1d0)*vcurrent(2)
        vcurrent(4)=-vcurrent(1)
 
      else if(   dabs( p2(1)+p2(4) ).lt.epsilon )then
 
        vcurrent(1)= dsqrt( 2d0*p2(1) / ( p1(1)+p1(4) ) ) *( p1(2) - (0d0,1d0)*p1(3) )
        vcurrent(2)= dsqrt( 2d0*p2(1) * ( p1(1)+p1(4) ) )
        vcurrent(3)=-(0d0,1d0)*vcurrent(2)
        vcurrent(4)=-vcurrent(1)
 
      else
 
        sqrt_pp1dpp2= dsqrt( (p1(1)+p1(4)) / (p2(1)+p2(4)) )
        sqrt_pp1xpp2= dsqrt( (p1(1)+p1(4)) * (p2(1)+p2(4)) )
        vcurrent(1)= sqrt_pp1xpp2  &
                    +( p1(2) - (0d0,1d0)*p1(3) )   &
                    *( p2(2) + (0d0,1d0)*p2(3) )/sqrt_pp1xpp2
        vcurrent(2)= ( p1(2) - (0d0,1d0)*p1(3) )/sqrt_pp1dpp2  &
                    +( p2(2) + (0d0,1d0)*p2(3) )*sqrt_pp1dpp2
        vcurrent(3)= ( (0d0,1d0)*p1(2) + p1(3) )/sqrt_pp1dpp2  &
                    -( (0d0,1d0)*p2(2) - p2(3) )*sqrt_pp1dpp2
        vcurrent(4)=sqrt_pp1xpp2  &
                    -( p1(2) - (0d0,1d0)*p1(3) )   &
                    *( p2(2) + (0d0,1d0)*p2(3) )/sqrt_pp1xpp2
      endif
 
      if( (dble(pdg_code1)*h1) .lt. 0d0)then
        do mu=1,4
          vcurrent(mu)=dconjg(vcurrent(mu))
        enddo
      endif
 
      return
      END subroutine ffv
 
 
!INV_LORENTZ.F
!VERSION 20130602
!
!A subroutine that performs a general inverse boost to a four vector
!(vector) based on another four vector (boost). The primed and
!unprimed frames have their axes in parallel to one another.
!Rotation is not performed by this subroutine.
      subroutine inv_lorentz(vector, boost)
 
      implicit none
 
      real(8) vector(4), boost(4)
      real(8) lambda(4,4), vector_copy(4)
      real(8) beta(2:4), beta_sq, gamma
      integer i,j
 
      do i=2,4
        beta(i) = -boost(i)/boost(1)
      enddo
 
      beta_sq = beta(2)**2+beta(3)**2+beta(4)**2
 
      gamma = 1d0/dsqrt(1d0-beta_sq)
 
      lambda(1,1) = gamma
 
      do i=2,4
        lambda(1,i) = gamma*beta(i)
        lambda(i,1) = lambda(1,i)
      enddo
 
      do i=2,4
      do j=2,4
        lambda(i,j) = (gamma-1d0)*beta(i)*beta(j)/beta_sq + kronecker_delta(i,j)
      enddo
      enddo
 
!apply boost to vector1
      vector_copy = vector
      vector = 0d0
      do i=1,4
      do j=1,4
        vector(i) = vector(i) + lambda(i,j)*vector_copy(j)
      enddo
      enddo
 
      return
      END subroutine inv_lorentz
 
!KRONECKER_DELTA.F
 
!KRONECKER_DELTA(i,j)
!A function that returns 1 if i=j, and 0 otherwise.
      real(8) function KRONECKER_DELTA(i,j)
      integer i,j
      if(i.eq.j)then
        kronecker_delta = 1d0
      else
        kronecker_delta = 0d0
      endif
 
      return
      end function kronecker_delta
 
!METRIC.F
!VERSION 20130524
 
!in METRIC returns the element of the Minkovski metric, with
!signature (1,-1,-1,-1), for given (_mu, _nu) or given (^mu, ^nu).
      real(8) function metric(mu,nu)
 
      implicit none
      integer mu, nu
 
      if(mu.ne.nu)then
        metric=0d0
      else if(mu.eq.1)then
        metric=1d0
      else
        metric=-1d0
      endif
 
      return
      END function metric
 
!`.F
!VERSION 20130524
 
!in POL(lambda,^mu) returns the polarization vectors for given
!4-momentum p
      subroutine polarization(p, POL)
 
      implicit none
      real(8) p(4)
      complex(8) POL(3,4)
 
      call polarization_single(p,+1,pol(1,:))
      call polarization_single(p,-1,pol(2,:))
      call polarization_single(p, 0,pol(3,:))
 
      return
      END subroutine polarization
 
      subroutine polarization_single(p, lambda, POL)
 
      implicit none
      real(8) p(4), sincos(4), inv_mass, abs3p
      complex(8) POL(4)
      integer lambda
 
      pol(:)=czero
 
      call angles(sincos, p)
      !sincos(1)=cos(theta)
      !sincos(2)=sin(theta)
      !sincos(3)=cos(phi)
      !sincos(4)=sin(phi)
 
!lambda = +1
      if(lambda.eq.1) then
         pol(1)= 0d0
         pol(2)= (sincos(3)*sincos(1)-(0d0,1d0)*sincos(4))/dsqrt(2d0)
         pol(3)= (sincos(4)*sincos(1)+(0d0,1d0)*sincos(3))/dsqrt(2d0)
         pol(4)= -sincos(2)/dsqrt(2d0)
!lambda = -1
      else if(lambda.eq.-1) then
         pol(1)= 0d0
         pol(2)= (sincos(3)*sincos(1)+(0d0,1d0)*sincos(4))/dsqrt(2d0)
         pol(3)= (sincos(4)*sincos(1)-(0d0,1d0)*sincos(3))/dsqrt(2d0)
         pol(4)= -sincos(2)/dsqrt(2d0)
!lambda = 0 (z)
      else if(lambda.eq.0) then
!|3-momentum|
         abs3p = dsqrt(p(2)**2+p(3)**2+p(4)**2)
!invariant mass
         inv_mass= dsqrt(p(1)**2-abs3p**2)
         pol(1)= abs3p/inv_mass
         pol(2)= sincos(3)*sincos(2)*p(1)/inv_mass
         pol(3)= sincos(4)*sincos(2)*p(1)/inv_mass
         pol(4)= sincos(1)*p(1)/inv_mass
!lambda = 2 (vec{p-hat})
      else if(lambda.eq.2) then
         pol(1)= 0d0
         pol(2)= sincos(3)*sincos(2)
         pol(3)= sincos(4)*sincos(2)
         pol(4)= sincos(1)
!lambda = -2 (-vec{p-hat})
      else if(lambda.eq.-2) then
         pol(1)= 0d0
         pol(2)= -sincos(3)*sincos(2)
         pol(3)= -sincos(4)*sincos(2)
         pol(4)= -sincos(1)
      endif
 
      return
      END subroutine polarization_single
 
!POLARIZATIONA.F
!VERSION 20130529
 
!in POL(lambda,^mu) returns the photon polarization vectors for given
!4-momentum p
      subroutine polarizationa(p, POL)
 
      implicit none
      real(8) p(4)
      complex(8) POL(2,4)
 
      call polarization_single(p,+1,pol(1,:))
      call polarization_single(p,-1,pol(2,:))
 
      return
      END subroutine polarizationa
 
!POLARIZATIONX.F
!VERSION 20130529
 
!in POL(lambda,^mu) returns the polarization vectors for given
!4-momentum p, with POL(L, ^mu) = (0, 0, 0, 1) when p is along the
!Z direction.
      subroutine polarizationx(p, POL)
 
      implicit none
      real(8) p(4)
      complex(8) POL(3,4)
 
      call polarization_single(p,+1,pol(1,:))
      call polarization_single(p,-1,pol(2,:))
      call polarization_single(p,+2,pol(3,:))
 
      return
      END subroutine polarizationx
!PROPAGATOR.F
!VERSION 20130522
!
!PROPAGATOR() returns the generi!complex-valued propagator
!without tensor structure (numerator), given mass, invariant mass
!and width.
      complex(8) function propagator(inv_mass, mass, width)
      implicit none
 
      real(8) inv_mass, mass, width
 
!not assuming auto-conversion
!     PROPAGATOR = (0d0,1d0)/(dcmplx(inv_mass**2,0d0)
!    &            -dcmplx(mass**2,0d0)+
!    &             (0d0,1d0)*dcmplx(mass,0d0)*dcmplx(width,0d0))
 
!assuming auto-conversion. works with gfortran
      !print *,"Called propagator with",inv_mass,mass,width
      if (mass.ge.0d0) then
         propagator = ci / ( inv_mass**2 - mass**2 + ci*mass*width )
      else
         propagator = ci / ( inv_mass**2 )
      endif
!     print *, PROPAGATOR
 
      return
      END function propagator
 
!VVP.F
!VERSION 20130524
 
!in epp(_mu, _nu) returns the
      subroutine vvp(p1,p2,epp)
 
      implicit none
      real(8) p1(4), p2(4)
      complex(8) epp(4,4)
!      real(8) ANTISYMMETRIC
      integer i,j,k,l
 
!      external ANTISYMMETRIC
 
      do i=1,4
      do j=1,4
        epp(i,j)=0d0
      enddo
      enddo
 
      do i=1,4
      do j=1,4
      do k=1,4
      do l=1,4
          epp(i,j)=epp(i,j)+antisymmetric(i,j,k,l)*p1(k)*p2(l)
      enddo
      enddo
      enddo
      enddo
 
      return
      END subroutine vvp
 
!VVS1.F
!VERSION 20130524
 
!in g_mu_nu(_mu, _nu) returns the VVS1 4*4 tensor.
      subroutine vvs1(g_mu_nu)
 
      implicit none
      complex(8) g_mu_nu(4,4)
 
      g_mu_nu = 0d0
 
      g_mu_nu(1,1) =  1d0
      g_mu_nu(2,2) = -1d0
      g_mu_nu(3,3) = -1d0
      g_mu_nu(4,4) = -1d0
 
      return
      END subroutine vvs1
 
!VVS2.F
!VERSION 20130524
 
!in pp(_mu, _nu) returns the tensor of p1_mu * p2_nu.
      subroutine vvs2(p1,p2,pp)
 
      implicit none
      real(8) p1(4), p2(4)
      complex(8) pp(4,4)
      integer mu, nu
 
      do mu=1,4
        do nu=1,4
          pp(mu,nu)=p1(mu)*p2(nu)
          if( ( (mu.ne.1).and.(nu.eq.1) ).or. &
              ( (mu.eq.1).and.(nu.ne.1) ) )then
            pp(mu,nu)=-pp(mu,nu)
          endif
        enddo
      enddo
 
      return
      END subroutine vvs2
 
  !-- generic functions below
  !- MCFM spinors in non-MCFM momentum convention
      subroutine spinoru2(n,p,za,zb,s)
       implicit none
       integer, intent(in) :: n
       real(8), intent(in) :: p(4,n)
       complex(8), intent(out) :: za(n,n), zb(n,n)
       real(8), intent(out) :: s(n,n)
       integer :: i,j
       complex(8) :: c23(n), f(n)
       real(8) :: rt(n)
 
       !---if one of the vectors happens to be zero this routine fails.
       do j=1,n
          za(j,j)=czero
          zb(j,j)=za(j,j)
 
          !-----positive energy case
          if (p(1,j) .gt. zero) then
             rt(j)=sqrt(abs(p(2,j)+p(1,j)))
             c23(j)=dcmplx(p(4,j),-p(3,j))
             f(j)=(one,zero)
          else
          !-----negative energy case
             rt(j)=sqrt(abs(-p(1,j)-p(2,j)))
             c23(j)=dcmplx(-p(4,j),p(3,j))
             f(j)=ci
          endif
       enddo
 
       do i=2,n
 
        do j=1,i-1
             s(i,j)=two*scr(p(:,i),p(:,j))
             za(i,j)=f(i)*f(j)  * ( c23(i)*dcmplx(rt(j)/(rt(i)+1d-16))-c23(j)*dcmplx(rt(i)/(rt(j)+1d-16)) )
 
             if (abs(s(i,j)).lt.1d-5) then
                zb(i,j)=-(f(i)*f(j))**2*conjg(za(i,j))
             else
                zb(i,j)=-dcmplx(s(i,j))/(za(i,j)+1d-16)
             endif
 
             za(j,i)=-za(i,j)
             zb(j,i)=-zb(i,j)
             s(j,i)=s(i,j)
 
          enddo
 
       enddo
 
       return
 
      end subroutine spinoru2
 
 
!!!!!!!!!!!!!!
! XCHECK MEs !
!!!!!!!!!!!!!!
 
!MATRIXELEMENT02.F
!VERSION 20160924
!complex(8) function MATRIXELEMENT02(p,mass,helicity,id)
!use ModHiggs
!implicit none
!real(8), intent(in) :: p(1:4,1:9)
!real(8), intent(in) :: mass(3:5,1:2)
!real(8), intent(in) :: helicity(9)
!integer, intent(in) :: id(9)
 
!complex(dp) :: A_VV(1:4), propH
!integer :: MY_IDUP(6:9),i3,i4,j,VVMode
!real(dp) :: pUsed(4,6)
!integer :: ordering(1:4),ordering_swap(1:4)
 
!! Set the ids and momenta
!MY_IDUP(:)=Not_a_particle_
!pUsed(:,:)=zero
!do j=1,4
!   if(id(j) .ne. Not_a_particle_) then
!      if(j.le.2) then
!         MY_IDUP(j+5) = -convertLHEreverse(id(j))
!         pUsed(:,j+2) = -p(:,j)
!      else
!         MY_IDUP(j+5) = convertLHEreverse(id(j+3))
!         pUsed(:,j+2) = p(:,j+3)
!      endif
!   endif
!   pUsed(:,1) = pUsed(:,1) + pUsed(:,j+2)
!enddo
!call getDecay_VVMode_Ordering(MY_IDUP, VVMode,ordering,ordering_swap)
 
!! Set the helicities
!if(id(1)*helicity(1).gt.0) then
!   i3=2
!else
!   i3=1
!endif
!if(id(6)*helicity(6).gt.0) then
!   i4=2
!else
!   i4=1
!endif
!if(ordering(1).eq.5 .or. ordering(2).eq.5) then
!   call swap(i3,i4)
!endif
 
!A_VV(:) = 0d0
!PROPH = PROPAGATOR(dsqrt(scr(p(:,5),p(:,5))),mass(5,1),mass(5,2))
!call calcHelAmp2(ordering,VVMode,MY_IDUP,pUsed,i3,i4,A_VV(1))
!if( (VVMode.eq.ZZMode) .and. includeGammaStar ) then
!    call calcHelAmp2(ordering,ZgsMode,MY_IDUP,pUsed,i3,i4,A_VV(2))
!    call calcHelAmp2(ordering,gsZMode,MY_IDUP,pUsed,i3,i4,A_VV(3))
!    call calcHelAmp2(ordering,gsgsMode,MY_IDUP,pUsed,i3,i4,A_VV(4))
!elseif( VVMode.eq.ZgMode .and. includeGammaStar ) then
!    call calcHelAmp2(ordering,gsgMode,MY_IDUP,pUsed,i3,i4,A_VV(2))
!endif
!A_VV(:) = -A_VV(:) * propH
 
!MATRIXELEMENT02 = A_VV(1)+A_VV(2)+A_VV(3)+A_VV(4)
!return
!END function
 
!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
 
 
SUBROUTINE matrixelement1(p,FermFlav,UnPolSqAmp)
use modparameters
use modmisc
implicit none
complex(8) :: SME(1:3,-1:+1,-1:+1),HelAmp
real(8) :: p(1:4,1:9),UnpolSqAmp,PreFac,IZis(-1:+1)
real(8) :: qsq_V1,qsq_V2,qsq_V1V2,qsq_H
complex(8) ghz1_dyn,ghz2_dyn,ghz3_dyn,ghz4_dyn
complex(8) :: a1HVV,a2HVV,a3HVV,Prop
integer :: ishel,fshel,FermFlav(1:6)! 12:IS, 34:ZDK, 56:HDK
 
    ! q1 qbar2 --> 3 --> 45 --> Z4-->f6 fbar7 + H5-->89
    call getsme(p,fermflav,sme)
    if( h_dk ) call error("Higgs decay not implemented")
 
    prop = (0d0,1d0)/(((p(1:4,4)+p(1:4,5)).dot.(p(1:4,4)+p(1:4,5))) -m_v**2 + (0d0,1d0)*m_v*ga_v )
    prefac = 4d0*pi*alpha_qed/4d0/sitw**2/(1d0-sitw**2)      ! gets squared below
 
    ! initial state couplings
    if( isawdecay(decaymode1) ) then
        izis(+1) = br   *ckm(fermflav(1),fermflav(2))
        izis(-1) = bl   *ckm(fermflav(1),fermflav(2))
    elseif( isazdecay(decaymode1) .and. (abs(fermflav(1)).eq.2 .or. abs(fermflav(1)).eq.4) ) then
        izis(+1) = ar_qup
        izis(-1) = al_qup
    elseif( isazdecay(decaymode1) .and. (abs(fermflav(1)).eq.1 .or. abs(fermflav(1)).eq.3 .or. abs(fermflav(1)).eq.5) ) then
        izis(+1) = ar_qdn
        izis(-1) = al_qdn
    endif
 
 
    ! anomalous HVV couplings
    qsq_v1  =  p(1:4,3).dot.p(1:4,3)
    qsq_v2  =  p(1:4,4).dot.p(1:4,4)
    qsq_v1v2=-(p(1:4,3).dot.p(1:4,4))
    qsq_h   =  p(1:4,5).dot.p(1:4,5)
 
    ghz1_dyn = hvvspinzerodynamiccoupling(1,qsq_v1,qsq_v2,qsq_h)
    ghz2_dyn = hvvspinzerodynamiccoupling(2,qsq_v1,qsq_v2,qsq_h)
    ghz3_dyn = hvvspinzerodynamiccoupling(3,qsq_v1,qsq_v2,qsq_h)
    ghz4_dyn = hvvspinzerodynamiccoupling(4,qsq_v1,qsq_v2,qsq_h)
 
    a1hvv = ghz1_dyn*m_v**2 + qsq_v1v2*( 2d0*ghz2_dyn + ghz3_dyn*qsq_v1v2/lambda**2 )
    a2hvv =-2d0*ghz2_dyn - ghz3_dyn*qsq_v1v2/lambda**2
    a3hvv =-2d0*ghz4_dyn
 
    unpolsqamp = 0d0
    do ishel=-1,+1,2
    do fshel=-1,+1,2
          helamp =  a3hvv * (  - izis(ishel)*sme(3,ishel,fshel)  )          &
                  + a2hvv * (  - izis(ishel)*sme(2,ishel,fshel)  )          &
                  + a1hvv * (  + izis(ishel)*sme(1,ishel,fshel)  )
          helamp = helamp * prefac/vev * prop
          unpolsqamp = unpolsqamp + dreal( helamp*dconjg(helamp) )
    enddo
    enddo
    unpolsqamp = unpolsqamp * cf
 
RETURN
END SUBROUTINE
 
SUBROUTINE getsme(p,FermFlav,SME)
use modparameters
use modmisc
implicit none
complex(8) :: SME(1:3,-1:+1,-1:+1)
real(8) :: sprod(9,9),p(1:4,1:9),IZfs(-1:+1)
complex(8) :: za(9,9), zb(9,9),Prop
integer :: FermFlav(1:6)
 
    call spinoru2(9,(/-p(1:4,1),-p(1:4,2),-p(1:4,1)-p(1:4,2),p(1:4,6)+p(1:4,7),p(1:4,8)+p(1:4,9),p(1:4,6),p(1:4,7),p(1:4,8),p(1:4,9)/),za,zb,sprod)
 
 
    ! Z-final state couplings
    if( isawdecay(decaymode1) ) then
        izfs(+1) = br   *ckm(fermflav(3),fermflav(4))
        izfs(-1) = bl   *ckm(fermflav(3),fermflav(4))
    elseif( abs(fermflav(3)).eq.11 .or. abs(fermflav(3)).eq.13 .or. abs(fermflav(3)).eq.15) then
         izfs(-1)=al_lep    * dsqrt(scale_alpha_z_ll)
         izfs(+1)=ar_lep    * dsqrt(scale_alpha_z_ll)
    elseif( abs(fermflav(3)).eq.12 .or. abs(fermflav(3)).eq.14 .or. abs(fermflav(3)).eq.16 ) then
         izfs(-1)=al_neu    * dsqrt(scale_alpha_z_nn)
         izfs(+1)=ar_neu    * dsqrt(scale_alpha_z_nn)
    elseif( abs(fermflav(3)).eq.2 .or. abs(fermflav(3)).eq.4 ) then
         izfs(-1)=al_qup    * dsqrt(scale_alpha_z_uu)
         izfs(+1)=ar_qup    * dsqrt(scale_alpha_z_uu)
    elseif( abs(fermflav(3)).eq.1 .or. abs(fermflav(3)).eq.3 .or. abs(fermflav(3)).eq.5 ) then
         izfs(-1)=al_qdn    * dsqrt(scale_alpha_z_dd)
         izfs(+1)=ar_qdn    * dsqrt(scale_alpha_z_dd)
    else
         call error("Wrong flavor in getSME",fermflav(3))
    endif
 
    sme(1,+1,+1) =  -2*izfs(1)*za(1,7)*zb(2,6)
    sme(2,+1,+1) =  izfs(1)*(za(1,6)*zb(2,6) + za(1,7)*zb(2,7))*(za(7,8)*zb(6,8) + za(7,9)*zb(6,9))
    sme(3,+1,+1) = ci*izfs(1)*(za(1,7)*(za(7,8)*zb(2,8) + za(7,9)*zb(2,9))*zb(6,7) + za(6,7)*zb(2,6)*(za(1,8)*zb(6,8) + za(1,9)*zb(6,9)))
    sme(1,+1,-1) =  -2*izfs(-1)*za(1,6)*zb(2,7)
    sme(2,+1,-1) =  izfs(-1)*(za(1,6)*zb(2,6) + za(1,7)*zb(2,7))*(za(6,8)*zb(7,8) + za(6,9)*zb(7,9))
    sme(3,+1,-1) = ci*izfs(-1)*(-(za(1,7)*zb(2,7)*(za(6,8)*zb(7,8) + za(6,9)*zb(7,9))) + za(1,6)*(za(6,8)*(-(zb(2,7)*zb(6,8)) + zb(2,6)*zb(7,8)) + zb(2,7)*(za(7,8)*zb(7,8) + za(7,9)*zb(7,9)) + za(6,9)*(-(zb(2,7)*zb(6,9)) + zb(2,6)*zb(7,9))))
    sme(1,-1,+1) = -2*izfs(1)*za(2,7)*zb(1,6)
    sme(2,-1,+1) = izfs(1)*(za(2,6)*zb(1,6) + za(2,7)*zb(1,7))*(za(7,8)*zb(6,8) + za(7,9)*zb(6,9))
    sme(3,-1,+1) = ci*izfs(1)*(za(2,7)*(za(7,8)*zb(1,8) + za(7,9)*zb(1,9))*zb(6,7) + za(6,7)*zb(1,6)*(za(2,8)*zb(6,8) + za(2,9)*zb(6,9)))
    sme(1,-1,-1) = -2*izfs(-1)*za(2,6)*zb(1,7)
    sme(2,-1,-1) =  izfs(-1)*(za(2,6)*zb(1,6) + za(2,7)*zb(1,7))*(za(6,8)*zb(7,8) + za(6,9)*zb(7,9))
    sme(3,-1,-1) =   -(ci*izfs(-1)*(za(2,6)*(za(6,8)*zb(1,8) + za(6,9)*zb(1,9))*zb(6,7) + za(6,7)*zb(1,7)*(za(2,8)*zb(7,8) + za(2,9)*zb(7,9))))
 
    prop = (0d0,1d0)/(2*(p(1:4,6).dot.p(1:4,7)) - m_v**2 + (0d0,1d0)*m_v*ga_v )
    sme(:,:,:) = sme(:,:,:) * prop
 
RETURN
END SUBROUTINE
 
 
 
 
end module modvhiggs
!!--YaofuZhou-----------------------------------------