modvhiggs Module Reference

Functions/Subroutines
subroutine, public	evalamp_vhiggs (id, helicity, MomExt, me2)
complex(8) function	matrixelement0 (MomExt, mass, helicity, id, useA)
complex(8) function	getvpffcoupling_vh (pdgid, hel, useWp)
subroutine	angles (sincos, vector)
subroutine	antisymmetric2 (p1, p2, epp)
real(8) function	antisymmetric (i, j, k, l)
subroutine	contra_field_tensor (POL, T_mu_nu)
subroutine	contra_outer (p1, p2, pp)
subroutine	covariant_field_tensor (POL, T_mu_nu)
subroutine	covariant_outer (p1, p2, pp)
complex(8) function	covariant_vector (p, mu)
complex(8) function	ffp (pdg_code1, p1, h1, pdg_code2, p2, h2)
subroutine	ffa (pdg_code1, p1, h1, pdg_code2, p2, h2, Acurrent)
complex(8) function	ffs (pdg_code1, p1, h1, pdg_code2, p2, h2)
subroutine	ffv (pdg_code1, p1, h1, pdg_code2, p2, h2, Vcurrent)
subroutine	inv_lorentz (vector, boost)
real(8) function	kronecker_delta (i, j)
real(8) function	metric (mu, nu)
subroutine	polarization (p, POL)
subroutine	polarization_single (p, lambda, POL)
subroutine	polarizationa (p, POL)
subroutine	polarizationx (p, POL)
complex(8) function	propagator (inv_mass, mass, width)
subroutine	vvp (p1, p2, epp)
subroutine	vvs1 (g_mu_nu)
subroutine	vvs2 (p1, p2, pp)
subroutine	spinoru2 (n, p, za, zb, s)
subroutine	matrixelement1 (p, FermFlav, UnPolSqAmp)
subroutine	getsme (p, FermFlav, SME)

Function/Subroutine Documentation

angles()

subroutine modvhiggs::angles ( real(8), dimension(4)  sincos, real(8), dimension(4)  vector  )

private

Definition at line 681 of file mod_VHiggs.F90.

      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

antisymmetric()

real(8) function modvhiggs::antisymmetric ( integer  i, integer  j, integer  k, integer  l  )

private

Definition at line 775 of file mod_VHiggs.F90.

 
      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

antisymmetric2()

subroutine modvhiggs::antisymmetric2 ( complex(8), dimension(4)  p1, complex(8), dimension(4)  p2, complex(8), dimension(4,4)  epp  )

private

Definition at line 733 of file mod_VHiggs.F90.

 
      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

contra_field_tensor()

subroutine modvhiggs::contra_field_tensor ( complex(8), dimension(3,4)  POL, complex(8), dimension(5,4,4)  T_mu_nu  )

private

Definition at line 792 of file mod_VHiggs.F90.

 
      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

contra_outer()

subroutine modvhiggs::contra_outer ( complex(8), dimension(4)  p1, complex(8), dimension(4)  p2, complex(8), dimension(4,4)  pp  )

private

Definition at line 835 of file mod_VHiggs.F90.

 
      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

covariant_field_tensor()

subroutine modvhiggs::covariant_field_tensor ( complex(8), dimension(3,4)  POL, complex(8), dimension(5,4,4)  T_mu_nu  )

private

Definition at line 857 of file mod_VHiggs.F90.

 
      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

covariant_outer()

subroutine modvhiggs::covariant_outer ( complex(8), dimension(4)  p1, complex(8), dimension(4)  p2, complex(8), dimension(4,4)  pp  )

private

Definition at line 899 of file mod_VHiggs.F90.

 
      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

covariant_vector()

complex(8) function modvhiggs::covariant_vector ( complex(8), dimension(4)  p, integer  mu  )

private

Definition at line 925 of file mod_VHiggs.F90.

 
      implicit none
 
      complex(8) p(4)
      integer mu
 
      if(mu.ne.1)then
        covariant_vector = -p(mu)
      else
        covariant_vector = p(mu)
      endif
 
      return

evalamp_vhiggs()

subroutine, public modvhiggs::evalamp_vhiggs	(	integer, dimension(9), intent(in)	id,
		real(8), dimension(9), intent(in)	helicity,
		real(8), dimension(1:4,1:9), intent(in)	MomExt,
		real(8), intent(out)	me2
	)

Definition at line 15 of file mod_VHiggs.F90.

      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

ffa()

subroutine modvhiggs::ffa ( integer  pdg_code1, real(8), dimension(4)  p1, real(8)  h1, integer  pdg_code2, real(8), dimension(4)  p2, real(8)  h2, complex(8), dimension(4)  Acurrent  )

private

Definition at line 983 of file mod_VHiggs.F90.

 
      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

ffp()

complex(8) function modvhiggs::ffp ( integer  pdg_code1, real(8), dimension(4)  p1, real(8)  h1, integer  pdg_code2, real(8), dimension(4)  p2, real(8)  h2  )

private

Definition at line 945 of file mod_VHiggs.F90.

 
      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

ffs()

complex(8) function modvhiggs::ffs ( integer  pdg_code1, real(8), dimension(4)  p1, real(8)  h1, integer  pdg_code2, real(8), dimension(4)  p2, real(8)  h2  )

private

Definition at line 1056 of file mod_VHiggs.F90.

 
      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

ffv()

subroutine modvhiggs::ffv ( integer  pdg_code1, real(8), dimension(4)  p1, real(8)  h1, integer  pdg_code2, real(8), dimension(4)  p2, real(8)  h2, complex(8), dimension(4)  Vcurrent  )

private

Definition at line 1094 of file mod_VHiggs.F90.

 
      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

getsme()

subroutine modvhiggs::getsme	(	real(8), dimension(1:4,1:9)	p,
		integer, dimension(1:6)	FermFlav,
		complex(8), dimension(1:3,-1:+1,-1:+1)	SME
	)

Definition at line 1682 of file mod_VHiggs.F90.

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

getvpffcoupling_vh()

complex(8) function modvhiggs::getvpffcoupling_vh ( integer, intent(in)  pdgid, integer, intent(in)  hel, logical, intent(in)  useWp  )

private

Definition at line 661 of file mod_VHiggs.F90.

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

inv_lorentz()

subroutine modvhiggs::inv_lorentz ( real(8), dimension(4)  vector, real(8), dimension(4)  boost  )

private

Definition at line 1166 of file mod_VHiggs.F90.

 
      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

kronecker_delta()

real(8) function modvhiggs::kronecker_delta ( integer  i, integer  j  )

private

Definition at line 1212 of file mod_VHiggs.F90.

      integer i,j
      if(i.eq.j)then
        kronecker_delta = 1d0
      else
        kronecker_delta = 0d0
      endif
 
      return

matrixelement0()

complex(8) function modvhiggs::matrixelement0 ( real(8), dimension(1:4,1:9), intent(in)  MomExt, real(8), dimension(3:5,1:2), intent(in)  mass, real(8), dimension(9), intent(in)  helicity, integer, dimension(9), intent(in)  id, logical, dimension(2), intent(in)  useA  )

private

Definition at line 104 of file mod_VHiggs.F90.

      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

matrixelement1()

subroutine modvhiggs::matrixelement1 ( real(8), dimension(1:4,1:9)  p, integer, dimension(1:6)  FermFlav, real(8)  UnPolSqAmp  )

private

Definition at line 1621 of file mod_VHiggs.F90.

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

metric()

real(8) function modvhiggs::metric ( integer  mu, integer  nu  )

private

Definition at line 1228 of file mod_VHiggs.F90.

 
      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

polarization()

subroutine modvhiggs::polarization ( real(8), dimension(4)  p, complex(8), dimension(3,4)  POL  )

private

Definition at line 1249 of file mod_VHiggs.F90.

 
      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

polarization_single()

subroutine modvhiggs::polarization_single ( real(8), dimension(4)  p, integer  lambda, complex(8), dimension(4)  POL  )

private

Definition at line 1262 of file mod_VHiggs.F90.

 
      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

polarizationa()

subroutine modvhiggs::polarizationa ( real(8), dimension(4)  p, complex(8), dimension(2,4)  POL  )

private

Definition at line 1321 of file mod_VHiggs.F90.

 
      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

polarizationx()

subroutine modvhiggs::polarizationx ( real(8), dimension(4)  p, complex(8), dimension(3,4)  POL  )

private

Definition at line 1339 of file mod_VHiggs.F90.

 
      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

propagator()

complex(8) function modvhiggs::propagator ( real(8)  inv_mass, real(8)  mass, real(8)  width  )

private

Definition at line 1357 of file mod_VHiggs.F90.

      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

spinoru2()

subroutine modvhiggs::spinoru2 ( integer, intent(in)  n, real(8), dimension(4,n), intent(in)  p, complex(8), dimension(n,n), intent(out)  za, complex(8), dimension(n,n), intent(out)  zb, real(8), dimension(n,n), intent(out)  s  )

private

Definition at line 1457 of file mod_VHiggs.F90.

       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
 

vvp()

subroutine modvhiggs::vvp ( real(8), dimension(4)  p1, real(8), dimension(4)  p2, complex(8), dimension(4,4)  epp  )

private

Definition at line 1383 of file mod_VHiggs.F90.

 
      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

vvs1()

subroutine modvhiggs::vvs1 ( complex(8), dimension(4,4)  g_mu_nu )

private

Definition at line 1416 of file mod_VHiggs.F90.

 
      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

vvs2()

subroutine modvhiggs::vvs2 ( real(8), dimension(4)  p1, real(8), dimension(4)  p2, complex(8), dimension(4,4)  pp  )

private

Definition at line 1435 of file mod_VHiggs.F90.

 
      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