modhiggsjj Module Reference

Functions/Subroutines
subroutine, public	evalamp_sbfh_unsymm_sa (p, res)
subroutine, public	evalamp_sbfh_unsymm_sa_select (p, iSel, jSel, flav_tag, iflip, res)
subroutine, public	evalamp_sbfh_unsymm_sa_select_exact (p, iSel, jSel, rSel, sSel, res)
subroutine, public	evalamp_wbfh_unsymm_sa (p, res)
subroutine, public	evalamp_wbfh_unsymm_sa_select (p, iSel, jSel, zz_fusion, iflip, res)
subroutine, public	evalamp_wbfh_unsymm_sa_select_exact (p, iSel, jSel, rSel, sSel, res)
integer function, dimension(2)	flip (i, a1, a2)
complex(dp) function, dimension(-1:1,-1:1)	a0_vv_4f (j1, j2, j3, j4, za, zb, sprod, mv, ga_v, useWWcoupl, Wpm_flip)
complex(dp) function, dimension(-1:1,-1:1)	a0_zz_4f (j1, j2, j3, j4, za, zb, sprod, iSel, jSel)
complex(dp) function, dimension(-1:1,-1:1)	a0_ww_4f (j1, j2, j3, j4, za, zb, sprod, iSel, jSel, useWWcoupl, Wpm_flip)
subroutine	me2_ggggh (j1, j2, j3, j4, za, zb, sprod, res)
subroutine	me2_qbqggh (j1, j2, j3, j4, za, zb, sprod, res)
subroutine	me2_qbqqbq (j1, j2, j3, j4, za, zb, sprod, res_diff, res_id)
complex(dp) function, dimension(1:2,-1:1,-1:1,-1:1,-1:1)	a0phigggg_xxxx (j1, j2, j3, j4, za, zb, sprod)
complex(dp) function, dimension(1:2,-1:1,-1:1,-1:1)	a0phiqbqgg_xxx (j1, j2, j3, j4, za, zb, sprod)
complex(dp) function, dimension(1:2,-1:1,-1:1)	a0phiqbqqbq_xx (j1, j2, j3, j4, za, zb, sprod)
complex(dp) function	a0phiggggpmmm (j1, j2, j3, j4, za, zb, sprod)
complex(dp) function	a0phiggggmpmp (j1, j2, j3, j4, za, zb, sprod)
complex(dp) function	a0phiggggmmpp (j1, j2, j3, j4, za, zb, sprod)
complex(dp) function	a0phiggggmmmm (j1, j2, j3, j4, za, zb, sprod)
subroutine	spinoru2 (n, p, za, zb, s)

Variables
real(dp), parameter, public	tag1 = 1.0_dp
real(dp), parameter, public	tag2 = 1.0_dp
real(dp), parameter, public	tagbot = 1.0_dp

Function/Subroutine Documentation

a0_vv_4f()

complex(dp) function, dimension(-1:1,-1:1) modhiggsjj::a0_vv_4f ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod, real(dp)  mv, real(dp)  ga_v, logical, optional  useWWcoupl, logical, optional  Wpm_flip  )

private

Definition at line 2267 of file mod_HiggsJJ.F90.

  use modmisc
  implicit none
    complex(dp) :: A0_VV_4f(-1:1,-1:1)
    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: mv, ga_v
    logical,optional :: useWWcoupl,Wpm_flip
    real(dp) :: sprod(4,4),q2Wplus,q2Wminus
    real(dp) :: mhsq, q1q2, kcoupl
    complex(dp) :: a1, a2, a3, struc1, struc2, struc3
    complex(dp) :: zab2
    complex(dp) :: iprop12, iprop34
    complex(dp) :: vvcoupl_prime(4)
    integer :: vv_it
    integer :: i,j,k,l
 
    zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
    !The previous line works, and assigns the whole zab2 correctly.
    !I have no idea how it works.
    !If you don't believe me, please uncomment the following lines:
    !do i=1,4
    ! do j=1,4
    !  do k=1,4
    !   do l=1,4
    !    print *,i,j,k,l,zab2(i,j,k,l),za(i,j)*zb(j,l) + za(i,k)*zb(k,l)
    !   enddo
    !  enddo
    ! enddo
    !enddo
    !print *,j1,j2,j3,j4
    !pause
 
    a0_vv_4f = czero
 
    q1q2 = (sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4))/two
    mhsq = two * q1q2 + sprod(j1,j2) + sprod(j3,j4)
 
    kcoupl = q1q2/lambda**2
 
    q2wplus  = sprod(j1,j2)
    q2wminus = sprod(j3,j4)
    if( present(wpm_flip) ) then
       if( wpm_flip ) call swapr(q2wplus,q2wminus)
    endif
    if( .not.present(usewwcoupl) ) then
       do vv_it=1,4
          vvcoupl_prime(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq)
       enddo
    else
       do vv_it=1,4
          vvcoupl_prime(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
       enddo
    endif
 
 
    a1 = vvcoupl_prime(1) * mv**2/mhsq + vvcoupl_prime(2) * two * q1q2/mhsq + vvcoupl_prime(3) * kcoupl * q1q2/mhsq
    a2 = -two * vvcoupl_prime(2) - kcoupl * vvcoupl_prime(3)
    a3 = -two * vvcoupl_prime(4)
 
    struc1 = two * (a1 * mhsq - ci * a3 * q1q2)
    struc2 = a2 + ci * a3
    struc3 = two * ci * a3
 
 
    a0_vv_4f(-1,-1) = za(j1,j3)*zb(j4,j2) * struc1 + &
         zab2(j1,j3,j4,j2)*zab2(j3,j1,j2,j4) * struc2 + &
         za(j1,j2)*za(j3,j4)*zb(j4,j2)**2 * struc3
 
    a0_vv_4f(-1,+1) = za(j1,j4)*zb(j3,j2) * struc1 + &
         zab2(j1,j3,j4,j2)*zab2(j4,j1,j2,j3) * struc2 + &
         za(j1,j2)*za(j4,j3)*zb(j3,j2)**2 * struc3
 
    a0_vv_4f(+1,-1) = za(j2,j3)*zb(j4,j1) * struc1 + &
         zab2(j2,j3,j4,j1)*zab2(j3,j1,j2,j4) * struc2 + &
         za(j2,j1)*za(j3,j4)*zb(j4,j1)**2 * struc3
 
    a0_vv_4f(+1,+1) = za(j2,j4)*zb(j3,j1) * struc1 + &
         zab2(j2,j3,j4,j1)*zab2(j4,j1,j2,j3) * struc2 + &
         za(j2,j1)*za(j4,j3)*zb(j3,j1)**2 * struc3
 
    iprop12 = sprod(j1,j2) - mv**2 + ci * mv * ga_v
    iprop34 = sprod(j3,j4) - mv**2 + ci * mv * ga_v
 
    a0_vv_4f = a0_vv_4f/vev /iprop12/iprop34
 
    return
 

a0_ww_4f()

complex(dp) function, dimension(-1:1,-1:1) modhiggsjj::a0_ww_4f	(	integer	j1,
		integer	j2,
		integer	j3,
		integer	j4,
		complex(dp), dimension(4,4)	za,
		complex(dp), dimension(4,4)	zb,
		real(dp), dimension(4,4)	sprod,
		integer	iSel,
		integer	jSel,
		logical, optional	useWWcoupl,
		logical, optional	Wpm_flip
	)

Definition at line 2640 of file mod_HiggsJJ.F90.

  use modmisc
  implicit none
    complex(dp) :: A0_WW_4f(-1:1,-1:1)
    integer :: j1,j2,j3,j4,iSel,jSel,line1,line2
    complex(dp) :: za(4,4), zb(4,4)
    logical,optional :: useWWcoupl,Wpm_flip
    real(dp) :: sprod(4,4),q2Wplus,q2Wminus
    real(dp) :: mhsq, q1q2, kcoupl
    complex(dp) :: vvcoupl_prime_ww(4)
    complex(dp) :: a1_ww, a2_ww, a3_ww, struc_ww(1:3)
    complex(dp) :: vvcoupl_prime_wwp(4)
    complex(dp) :: a1_wwp, a2_wwp, a3_wwp, struc_wwp(1:3)
    complex(dp) :: vvcoupl_prime_wpw(4)
    complex(dp) :: a1_wpw, a2_wpw, a3_wpw, struc_wpw(1:3)
    complex(dp) :: vvcoupl_prime_wpwp(4)
    complex(dp) :: a1_wpwp, a2_wpwp, a3_wpwp, struc_wpwp(1:3)
    complex(dp) :: zab2
    complex(dp) :: iprop12, iprop34, wpprop12, wpprop34
    complex(dp) :: helcoup(1:3,-1:1,-1:1)
    real(dp), dimension(5) :: LCT,RCT
    integer :: vv_it
    integer :: i,j,k,l
 
    zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
 
    a0_ww_4f = czero
    helcoup(:,:,:)=czero
 
    if( abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_ ) return
    if( abs(isel).eq.pdfglu_ .or. abs(jsel).eq.pdfglu_ ) return
 
    line1=abs(isel)
    line2=abs(jsel)
 
    q1q2 = (sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4))/two
    mhsq = two * q1q2 + sprod(j1,j2) + sprod(j3,j4)
 
    kcoupl = q1q2/lambda**2
 
    q2wplus  = sprod(j1,j2)
    q2wminus = sprod(j3,j4)
    if( present(wpm_flip) ) then
       if( wpm_flip ) then
          call swap(q2wplus,q2wminus)
          call swap(line1,line2)
       endif
    endif
    if( .not.present(usewwcoupl) ) then
       do vv_it=1,4
          vvcoupl_prime_ww(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq)
       enddo
    else
       do vv_it=1,4
          vvcoupl_prime_ww(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
       enddo
    endif
 
    ! SM contribution
    iprop12 = q2wplus - m_w**2 + ci * m_w * ga_w
    iprop34 = q2wminus - m_w**2 + ci * m_w * ga_w
 
    a1_ww = vvcoupl_prime_ww(1) * m_w**2/mhsq + vvcoupl_prime_ww(2) * two * q1q2/mhsq + vvcoupl_prime_ww(3) * kcoupl * q1q2/mhsq
    a2_ww = -two * vvcoupl_prime_ww(2) - kcoupl * vvcoupl_prime_ww(3)
    a3_ww = -two * vvcoupl_prime_ww(4)
 
    struc_ww(1) = two * (a1_ww * mhsq - ci * a3_ww * q1q2)
    struc_ww(2) = a2_ww + ci * a3_ww
    struc_ww(3) = two * ci * a3_ww
    struc_ww(:) = struc_ww(:) * couplwffsq
 
    helcoup(1:3,-1,-1) = struc_ww(1:3)/iprop12/iprop34
 
    if ( includevprime ) then
 
       lct = (/ewp_up_left,ewp_up_left,ewp_chm_left,ewp_chm_left,ewp_top_left/)
       rct = (/ewp_up_right,ewp_up_right,ewp_chm_right,ewp_chm_right,ewp_top_right/)
       ! Adjust couplings convention
       lct(:) = lct(:)/bl
       rct(:) = rct(:)/bl
 
       do vv_it=1,4
          vvcoupl_prime_wwp(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
          vvcoupl_prime_wpw(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q2wminus,q2wplus,mhsq,trywwcoupl=usewwcoupl)
          vvcoupl_prime_wpwp(vv_it) = hvvspinzerodynamiccoupling(15+vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
       enddo
 
 
       if (m_wprime.gt.0d0) then
          wpprop12 = q2wplus - m_wprime**2 + ci * m_wprime * ga_wprime
          wpprop34 = q2wminus - m_wprime**2 + ci * m_wprime * ga_wprime
       elseif (m_wprime.eq.0d0) then
          wpprop12 = q2wplus
          wpprop34 = q2wminus
       else
          wpprop12 = m_w**2
          wpprop34 = m_w**2
       endif
 
 
       a1_wwp = vvcoupl_prime_wwp(1) * m_w**2/mhsq + vvcoupl_prime_wwp(2) * two * q1q2/mhsq + vvcoupl_prime_wwp(3) * kcoupl * q1q2/mhsq
       a2_wwp = -two * vvcoupl_prime_wwp(2) - kcoupl * vvcoupl_prime_wwp(3)
       a3_wwp = -two * vvcoupl_prime_wwp(4)
       a1_wpw = vvcoupl_prime_wpw(1) * m_w**2/mhsq + vvcoupl_prime_wpw(2) * two * q1q2/mhsq + vvcoupl_prime_wpw(3) * kcoupl * q1q2/mhsq
       a2_wpw = -two * vvcoupl_prime_wpw(2) - kcoupl * vvcoupl_prime_wpw(3)
       a3_wpw = -two * vvcoupl_prime_wpw(4)
       a1_wpwp = vvcoupl_prime_wpwp(1) * m_w**2/mhsq + vvcoupl_prime_wpwp(2) * two * q1q2/mhsq + vvcoupl_prime_wpwp(3) * kcoupl * q1q2/mhsq
       a2_wpwp = -two * vvcoupl_prime_wpwp(2) - kcoupl * vvcoupl_prime_wpwp(3)
       a3_wpwp = -two * vvcoupl_prime_wpwp(4)
 
       struc_wwp(1) = two * (a1_wwp * mhsq - ci * a3_wwp * q1q2)
       struc_wwp(2) = (a2_wwp + ci * a3_wwp)
       struc_wwp(3) = two * ci * a3_wwp
       struc_wpw(1) = two * (a1_wpw * mhsq - ci * a3_wpw * q1q2)
       struc_wpw(2) = (a2_wpw + ci * a3_wpw)
       struc_wpw(3) = two * ci * a3_wpw
       struc_wpwp(1) = two * (a1_wpwp * mhsq - ci * a3_wpwp * q1q2)
       struc_wpwp(2) = (a2_wpwp + ci * a3_wpwp)
       struc_wpwp(3) = two * ci * a3_wpwp
 
       struc_wwp(:) = struc_wwp(:) * couplwffsq
       struc_wpw(:) = struc_wpw(:) * couplwffsq
       struc_wpwp(:) = struc_wpwp(:) * couplwffsq
 
       helcoup(1:3,-1,-1) = helcoup(1:3,-1,-1) + &
                            struc_wpwp(1:3) * lct(line1) * lct(line2)/wpprop12/wpprop34 + &
                            struc_wpw(1:3)  * lct(line1)             /wpprop12/iprop34  + &
                            struc_wwp(1:3)  *              lct(line2)/iprop12/wpprop34
 
       helcoup(1:3,-1,+1) = helcoup(1:3,-1,+1) + &
                            struc_wpwp(1:3) * lct(line1) * rct(line2)/wpprop12/wpprop34 + &
                            struc_wpw(1:3)  * lct(line1) * czero     /wpprop12/iprop34  + &
                            struc_wwp(1:3)               * rct(line2)/iprop12/wpprop34
 
       helcoup(1:3,+1,-1) = helcoup(1:3,+1,-1) + &
                            struc_wpwp(1:3) * rct(line1) * lct(line2)/wpprop12/wpprop34 + &
                            struc_wpw(1:3)  * rct(line1)             /wpprop12/iprop34  + &
                            struc_wwp(1:3)  * czero      * lct(line2)/iprop12/wpprop34
 
       helcoup(1:3,+1,+1) = helcoup(1:3,+1,+1) + &
                            struc_wpwp(1:3) * rct(line1) * rct(line2)/wpprop12/wpprop34 + &
                            struc_wpw(1:3)  * rct(line1) * czero     /wpprop12/iprop34  + &
                            struc_wwp(1:3)  * czero      * rct(line2)/iprop12/wpprop34
    endif
 
 
    a0_ww_4f(-1,-1) = za(j1,j3)*zb(j4,j2) * helcoup(1,-1,-1) + &
         zab2(j1,j3,j4,j2)*zab2(j3,j1,j2,j4) * helcoup(2,-1,-1) + &
         za(j1,j2)*za(j3,j4)*zb(j4,j2)**2 * helcoup(3,-1,-1)
 
    if ( includevprime ) then
       a0_ww_4f(-1,+1) = za(j1,j4)*zb(j3,j2) * helcoup(1,-1,+1) + &
            zab2(j1,j3,j4,j2)*zab2(j4,j1,j2,j3) * helcoup(2,-1,+1) + &
            za(j1,j2)*za(j4,j3)*zb(j3,j2)**2 * helcoup(3,-1,+1)
 
       a0_ww_4f(+1,-1) = za(j2,j3)*zb(j4,j1) * helcoup(1,+1,-1) + &
            zab2(j2,j3,j4,j1)*zab2(j3,j1,j2,j4) * helcoup(2,+1,-1) + &
            za(j2,j1)*za(j3,j4)*zb(j4,j1)**2 * helcoup(3,+1,-1)
 
       a0_ww_4f(+1,+1) = za(j2,j4)*zb(j3,j1) * helcoup(1,+1,+1) + &
            zab2(j2,j3,j4,j1)*zab2(j4,j1,j2,j3) * helcoup(2,+1,+1) + &
            za(j2,j1)*za(j4,j3)*zb(j3,j1)**2 * helcoup(3,+1,+1)
    endif
 
 
    a0_ww_4f = a0_ww_4f/vev
 
 
    return
 

a0_zz_4f()

complex(dp) function, dimension(-1:1,-1:1) modhiggsjj::a0_zz_4f	(	integer	j1,
		integer	j2,
		integer	j3,
		integer	j4,
		complex(dp), dimension(4,4)	za,
		complex(dp), dimension(4,4)	zb,
		real(dp), dimension(4,4)	sprod,
		integer	iSel,
		integer	jSel
	)

Definition at line 2358 of file mod_HiggsJJ.F90.

    use modmisc
    implicit none
    real(dp), dimension(5) :: Lz,Rz,LCT,RCT
    real(dp), parameter, dimension(5) :: La = (/qdl,qul,qdl,qul,qdl/)
    real(dp), parameter, dimension(5) :: Ra = (/qdr,qur,qdr,qur,qdr/)
    complex(dp) :: A0_ZZ_4f(-1:1,-1:1)
    integer :: j1,j2,j3,j4,iSel,jSel,line1,line2
    complex(dp) :: za(4,4),zb(4,4)
    real(dp) :: sprod(4,4)
    real(dp) :: mhsq,q1q2,kcoupl,q12sq,q34sq
    complex(dp) :: a1_zz,a2_zz,a3_zz,struc_zz(3),fac_zz(-1:1,-1:1)
    complex(dp) :: a1_az,a2_az,a3_az,struc_az(3),fac_az(-1:1,-1:1)
    complex(dp) :: a1_za,a2_za,a3_za,struc_za(3),fac_za(-1:1,-1:1)
    complex(dp) :: a1_aa,a2_aa,a3_aa,struc_aa(3),fac_aa(-1:1,-1:1)
    complex(dp) :: a1_zzp,a2_zzp,a3_zzp,struc_zzp(3),fac_zzp(-1:1,-1:1)
    complex(dp) :: a1_zpz,a2_zpz,a3_zpz,struc_zpz(3),fac_zpz(-1:1,-1:1)
    complex(dp) :: a1_azp,a2_azp,a3_azp,struc_azp(3),fac_azp(-1:1,-1:1)
    complex(dp) :: a1_zpa,a2_zpa,a3_zpa,struc_zpa(3),fac_zpa(-1:1,-1:1)
    complex(dp) :: a1_zpzp,a2_zpzp,a3_zpzp,struc_zpzp(3),fac_zpzp(-1:1,-1:1)
    complex(dp) :: helcoup(1:3,-1:1,-1:1)
    complex(dp) :: zab2
    complex(dp) :: iprop12,iprop34, zpprop12,zpprop34
    complex(dp) :: vvcoupl_prime_zz(4),vvcoupl_prime_az(4),vvcoupl_prime_za(4),vvcoupl_prime_aa(2:4)
    complex(dp) :: vvcoupl_prime_zzp(4),vvcoupl_prime_zpz(4),vvcoupl_prime_zpzp(4)
    complex(dp) :: vvcoupl_prime_azp(4),vvcoupl_prime_zpa(4)
    integer :: vv_it
    integer :: i,j,k,l
 
    zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
 
    a0_zz_4f = czero
    helcoup(:,:,:)=czero
    struc_zz(:)=czero
    struc_az(:)=czero
    struc_za(:)=czero
    struc_aa(:)=czero
    struc_zzp(:)=czero
    struc_zpz(:)=czero
    struc_azp(:)=czero
    struc_zpa(:)=czero
    struc_zpzp(:)=czero
    fac_zz(:,:)=czero
    fac_az(:,:)=czero
    fac_za(:,:)=czero
    fac_aa(:,:)=czero
    fac_zzp(:,:)=czero
    fac_zpz(:,:)=czero
    fac_azp(:,:)=czero
    fac_zpa(:,:)=czero
    fac_zpzp(:,:)=czero
 
    lz = (/al_qdn,al_qup,al_qdn,al_qup,al_qdn/)
    rz = (/ar_qdn,ar_qup,ar_qdn,ar_qup,ar_qdn/)
 
    if( abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_ ) return
    if( abs(isel).eq.pdfglu_ .or. abs(jsel).eq.pdfglu_ ) return
 
    line1=abs(isel)
    line2=abs(jsel)
 
    q1q2 = (sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4))/two
    mhsq = two * q1q2 + sprod(j1,j2) + sprod(j3,j4)
 
    kcoupl = q1q2/lambda**2
 
    q12sq = sprod(j1,j2)
    q34sq = sprod(j3,j4)
 
    iprop12 = q12sq - m_z**2 + ci * m_z * ga_z
    iprop34 = q34sq - m_z**2 + ci * m_z * ga_z
 
    !-- set up couplings
    do vv_it=1,4
       vvcoupl_prime_zz(vv_it) = hvvspinzerodynamiccoupling(vv_it,q12sq,q34sq,mhsq)
    enddo
 
    a1_zz = vvcoupl_prime_zz(1) * m_z**2/mhsq + vvcoupl_prime_zz(2) * two * q1q2/mhsq + vvcoupl_prime_zz(3) * kcoupl * q1q2/mhsq
    a2_zz = -two * vvcoupl_prime_zz(2) - kcoupl * vvcoupl_prime_zz(3)
    a3_zz = -two * vvcoupl_prime_zz(4)
 
    struc_zz(1) = two * (a1_zz * mhsq - ci * a3_zz * q1q2)
    struc_zz(2) = (a2_zz + ci * a3_zz)
    struc_zz(3) = two * ci * a3_zz
    struc_zz(:) = struc_zz(:) * couplzffsq
 
    fac_zz(-1,-1) = lz(line1) * lz(line2)/iprop12/iprop34
    fac_zz(-1,+1) = lz(line1) * rz(line2)/iprop12/iprop34
    fac_zz(+1,-1) = rz(line1) * lz(line2)/iprop12/iprop34
    fac_zz(+1,+1) = rz(line1) * rz(line2)/iprop12/iprop34
 
    if( includegammastar ) then
 
       do vv_it=1,4
          vvcoupl_prime_az(vv_it) = hvvspinzerodynamiccoupling(vv_it+4,0d0,q12sq,mhsq)
          vvcoupl_prime_za(vv_it) = hvvspinzerodynamiccoupling(vv_it+4,0d0,q34sq,mhsq)
       enddo
       do vv_it=1,3
          vvcoupl_prime_aa(vv_it+1) = hvvspinzerodynamiccoupling(vv_it+8,q12sq,q34sq,mhsq)
       enddo
 
       a1_aa = vvcoupl_prime_aa(2) * two * q1q2/mhsq + vvcoupl_prime_aa(3) * kcoupl * q1q2/mhsq
       a2_aa = -two * vvcoupl_prime_aa(2) - kcoupl * vvcoupl_prime_aa(3)
       a3_aa = -two * vvcoupl_prime_aa(4)
 
       a1_za = vvcoupl_prime_za(1) * m_z**2/mhsq + vvcoupl_prime_za(2) * two * q1q2/mhsq + vvcoupl_prime_za(3) * kcoupl * q1q2/mhsq
       a2_za = -two * vvcoupl_prime_za(2) - kcoupl * vvcoupl_prime_za(3)
       a3_za = -two * vvcoupl_prime_za(4)
 
       a1_az = vvcoupl_prime_az(1) * m_z**2/mhsq + vvcoupl_prime_az(2) * two * q1q2/mhsq + vvcoupl_prime_az(3) * kcoupl * q1q2/mhsq
       a2_az = -two * vvcoupl_prime_az(2) - kcoupl * vvcoupl_prime_az(3)
       a3_az = -two * vvcoupl_prime_az(4)
 
       struc_aa(1) = two * (a1_aa * mhsq - ci * a3_aa * q1q2)
       struc_aa(2) = (a2_aa + ci * a3_aa)
       struc_aa(3) = two * ci * a3_aa
       struc_aa(:) = struc_aa(:) * couplaffsq
 
       struc_az(1) = two * (a1_az * mhsq - ci * a3_az * q1q2)
       struc_az(2) = (a2_az + ci * a3_az)
       struc_az(3) = two * ci * a3_az
       struc_az(:) = struc_az(:) * couplazff
 
       struc_za(1) = two * (a1_za * mhsq - ci * a3_za * q1q2)
       struc_za(2) = (a2_za + ci * a3_za)
       struc_za(3) = two * ci * a3_za
       struc_za(:) = struc_za(:) * couplazff
 
       fac_aa(-1,-1) = la(line1) * la(line2)/q12sq/q34sq
       fac_aa(-1,+1) = la(line1) * ra(line2)/q12sq/q34sq
       fac_aa(+1,-1) = ra(line1) * la(line2)/q12sq/q34sq
       fac_aa(+1,+1) = ra(line1) * ra(line2)/q12sq/q34sq
 
       fac_az(-1,-1) = la(line1) * lz(line2)/q12sq/iprop34
       fac_az(-1,+1) = la(line1) * rz(line2)/q12sq/iprop34
       fac_az(+1,-1) = ra(line1) * lz(line2)/q12sq/iprop34
       fac_az(+1,+1) = ra(line1) * rz(line2)/q12sq/iprop34
 
       fac_za(-1,-1) = lz(line1) * la(line2)/iprop12/q34sq
       fac_za(-1,+1) = lz(line1) * ra(line2)/iprop12/q34sq
       fac_za(+1,-1) = rz(line1) * la(line2)/iprop12/q34sq
       fac_za(+1,+1) = rz(line1) * ra(line2)/iprop12/q34sq
 
    endif
 
    if ( includevprime ) then
 
       lct = (/ezp_dn_left,ezp_up_left,ezp_str_left,ezp_chm_left,ezp_bot_left/)
       rct = (/ezp_dn_right,ezp_up_right,ezp_str_right,ezp_chm_right,ezp_bot_right/)
 
       do vv_it=1,4
          vvcoupl_prime_zzp(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q12sq,q34sq,mhsq)
          vvcoupl_prime_zpz(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q34sq,q12sq,mhsq)
          vvcoupl_prime_zpzp(vv_it) = hvvspinzerodynamiccoupling(15+vv_it,q12sq,q34sq,mhsq)
       enddo
 
       if (m_zprime.gt.0d0) then
          zpprop12 = q12sq - m_zprime**2 + ci * m_zprime * ga_zprime
          zpprop34 = q34sq - m_zprime**2 + ci * m_zprime * ga_zprime
       elseif (m_zprime.eq.0d0) then
          zpprop12 = q12sq
          zpprop34 = q34sq
       else
          zpprop12 = m_z**2
          zpprop34 = m_z**2
       endif
 
 
       a1_zzp = vvcoupl_prime_zzp(1) * m_z**2/mhsq + vvcoupl_prime_zzp(2) * two * q1q2/mhsq + vvcoupl_prime_zzp(3) * kcoupl * q1q2/mhsq
       a2_zzp = -two * vvcoupl_prime_zzp(2) - kcoupl * vvcoupl_prime_zzp(3)
       a3_zzp = -two * vvcoupl_prime_zzp(4)
       a1_zpz = vvcoupl_prime_zpz(1) * m_z**2/mhsq + vvcoupl_prime_zpz(2) * two * q1q2/mhsq + vvcoupl_prime_zpz(3) * kcoupl * q1q2/mhsq
       a2_zpz = -two * vvcoupl_prime_zpz(2) - kcoupl * vvcoupl_prime_zpz(3)
       a3_zpz = -two * vvcoupl_prime_zpz(4)
       a1_zpzp = vvcoupl_prime_zpzp(1) * m_z**2/mhsq + vvcoupl_prime_zpzp(2) * two * q1q2/mhsq + vvcoupl_prime_zpzp(3) * kcoupl * q1q2/mhsq
       a2_zpzp = -two * vvcoupl_prime_zpzp(2) - kcoupl * vvcoupl_prime_zpzp(3)
       a3_zpzp = -two * vvcoupl_prime_zpzp(4)
 
       struc_zzp(1) = two * (a1_zzp * mhsq - ci * a3_zzp * q1q2)
       struc_zzp(2) = (a2_zzp + ci * a3_zzp)
       struc_zzp(3) = two * ci * a3_zzp
       struc_zpz(1) = two * (a1_zpz * mhsq - ci * a3_zpz * q1q2)
       struc_zpz(2) = (a2_zpz + ci * a3_zpz)
       struc_zpz(3) = two * ci * a3_zpz
       struc_zpzp(1) = two * (a1_zpzp * mhsq - ci * a3_zpzp * q1q2)
       struc_zpzp(2) = (a2_zpzp + ci * a3_zpzp)
       struc_zpzp(3) = two * ci * a3_zpzp
 
       struc_zzp(:) = struc_zzp(:) * couplzffsq
       struc_zpz(:) = struc_zpz(:) * couplzffsq
       struc_zpzp(:) = struc_zpzp(:) * couplzffsq
 
       fac_zpzp(-1,-1) = lct(line1) * lct(line2)/zpprop12/zpprop34
       fac_zpzp(-1,+1) = lct(line1) * rct(line2)/zpprop12/zpprop34
       fac_zpzp(+1,-1) = rct(line1) * lct(line2)/zpprop12/zpprop34
       fac_zpzp(+1,+1) = rct(line1) * rct(line2)/zpprop12/zpprop34
 
       fac_zpz(-1,-1) = lct(line1) * lz(line2)/zpprop12/iprop34
       fac_zpz(-1,+1) = lct(line1) * rz(line2)/zpprop12/iprop34
       fac_zpz(+1,-1) = rct(line1) * lz(line2)/zpprop12/iprop34
       fac_zpz(+1,+1) = rct(line1) * rz(line2)/zpprop12/iprop34
 
       fac_zzp(-1,-1) = lz(line1) * lct(line2)/iprop12/zpprop34
       fac_zzp(-1,+1) = lz(line1) * rct(line2)/iprop12/zpprop34
       fac_zzp(+1,-1) = rz(line1) * lct(line2)/iprop12/zpprop34
       fac_zzp(+1,+1) = rz(line1) * rct(line2)/iprop12/zpprop34
 
 
       if( includegammastar ) then
          do vv_it=1,4
             vvcoupl_prime_azp(vv_it) = hvvspinzerodynamiccoupling(19+vv_it,q12sq,q34sq,mhsq)
             vvcoupl_prime_zpa(vv_it) = hvvspinzerodynamiccoupling(19+vv_it,q34sq,q12sq,mhsq)
          enddo
 
          a1_azp = vvcoupl_prime_azp(1) * m_z**2/mhsq + vvcoupl_prime_azp(2) * two * q1q2/mhsq + vvcoupl_prime_azp(3) * kcoupl * q1q2/mhsq
          a2_azp = -two * vvcoupl_prime_azp(2) - kcoupl * vvcoupl_prime_azp(3)
          a3_azp = -two * vvcoupl_prime_azp(4)
 
          a1_zpa = vvcoupl_prime_zpa(1) * m_z**2/mhsq + vvcoupl_prime_zpa(2) * two * q1q2/mhsq + vvcoupl_prime_zpa(3) * kcoupl * q1q2/mhsq
          a2_zpa = -two * vvcoupl_prime_zpa(2) - kcoupl * vvcoupl_prime_zpa(3)
          a3_zpa = -two * vvcoupl_prime_zpa(4)
 
          struc_azp(1) = two * (a1_azp * mhsq - ci * a3_azp * q1q2)
          struc_azp(2) = (a2_azp + ci * a3_azp)
          struc_azp(3) = two * ci * a3_azp
          struc_azp(:) = struc_azp(:) * couplazff
 
          struc_zpa(1) = two * (a1_zpa * mhsq - ci * a3_zpa * q1q2)
          struc_zpa(2) = (a2_zpa + ci * a3_zpa)
          struc_zpa(3) = two * ci * a3_zpa
          struc_zpa(:) = struc_zpa(:) * couplazff
 
          fac_zpa(-1,-1) = lct(line1) * la(line2)/zpprop12/q34sq
          fac_zpa(-1,+1) = lct(line1) * ra(line2)/zpprop12/q34sq
          fac_zpa(+1,-1) = rct(line1) * la(line2)/zpprop12/q34sq
          fac_zpa(+1,+1) = rct(line1) * ra(line2)/zpprop12/q34sq
 
          fac_azp(-1,-1) = la(line1) * lct(line2)/q12sq/zpprop34
          fac_azp(-1,+1) = la(line1) * rct(line2)/q12sq/zpprop34
          fac_azp(+1,-1) = ra(line1) * lct(line2)/q12sq/zpprop34
          fac_azp(+1,+1) = ra(line1) * rct(line2)/q12sq/zpprop34
       endif
 
    endif
 
    do i=-1,1,2
    do j=-1,1,2
       helcoup(1:3,i,j) = struc_zz(1:3) * fac_zz(i,j) + &
                          struc_az(1:3) * fac_az(i,j) + &
                          struc_za(1:3) * fac_za(i,j) + &
                          struc_zpz(1:3) * fac_zpz(i,j) + &
                          struc_zzp(1:3) * fac_zzp(i,j) + &
                          struc_zpa(1:3) * fac_zpa(i,j) + &
                          struc_azp(1:3) * fac_azp(i,j) + &
                          struc_aa(1:3) * fac_aa(i,j) + &
                          struc_zpzp(1:3) * fac_zpzp(i,j)
    enddo
    enddo
 
 
    a0_zz_4f(-1,-1) = za(j1,j3)*zb(j4,j2) * helcoup(1,-1,-1) + &
         zab2(j1,j3,j4,j2)*zab2(j3,j1,j2,j4) * helcoup(2,-1,-1) + &
         za(j1,j2)*za(j3,j4)*zb(j4,j2)**2 * helcoup(3,-1,-1)
 
    a0_zz_4f(-1,+1) = za(j1,j4)*zb(j3,j2) * helcoup(1,-1,+1) + &
         zab2(j1,j3,j4,j2)*zab2(j4,j1,j2,j3) * helcoup(2,-1,+1) + &
         za(j1,j2)*za(j4,j3)*zb(j3,j2)**2 * helcoup(3,-1,+1)
 
    a0_zz_4f(+1,-1) = za(j2,j3)*zb(j4,j1) * helcoup(1,+1,-1) + &
         zab2(j2,j3,j4,j1)*zab2(j3,j1,j2,j4) * helcoup(2,+1,-1) + &
         za(j2,j1)*za(j3,j4)*zb(j4,j1)**2 * helcoup(3,+1,-1)
 
    a0_zz_4f(+1,+1) = za(j2,j4)*zb(j3,j1) * helcoup(1,+1,+1) + &
         zab2(j2,j3,j4,j1)*zab2(j4,j1,j2,j3) * helcoup(2,+1,+1) + &
         za(j2,j1)*za(j4,j3)*zb(j3,j1)**2 * helcoup(3,+1,+1)
 
    a0_zz_4f = a0_zz_4f/vev
 
    return
 

a0phigggg_xxxx()

complex(dp) function, dimension(1:2,-1:1,-1:1,-1:1,-1:1) modhiggsjj::a0phigggg_xxxx ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod  )

private

Definition at line 2946 of file mod_HiggsJJ.F90.

    complex(dp) :: A0phigggg_xxxx(1:2,-1:1,-1:1,-1:1,-1:1)
    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4),zb(4,4)
    real(dp) :: sprod(4,4)
 
    a0phigggg_xxxx = czero
 
    a0phigggg_xxxx(1,+1,+1,-1,-1) = a0phiggggmmpp(j3,j4,j1,j2,za,zb,sprod)
    a0phigggg_xxxx(1,+1,-1,+1,-1) = a0phiggggmpmp(j2,j3,j4,j1,za,zb,sprod)
    a0phigggg_xxxx(1,+1,-1,-1,+1) = a0phiggggmmpp(j2,j3,j4,j1,za,zb,sprod)
    a0phigggg_xxxx(1,+1,-1,-1,-1) = a0phiggggpmmm(j1,j2,j3,j4,za,zb,sprod)
    a0phigggg_xxxx(1,-1,+1,+1,-1) = a0phiggggmmpp(j4,j1,j2,j3,za,zb,sprod)
    a0phigggg_xxxx(1,-1,+1,-1,+1) = a0phiggggmpmp(j1,j2,j3,j4,za,zb,sprod)
    a0phigggg_xxxx(1,-1,+1,-1,-1) = a0phiggggpmmm(j2,j3,j4,j1,za,zb,sprod)
    a0phigggg_xxxx(1,-1,-1,+1,+1) = a0phiggggmmpp(j1,j2,j3,j4,za,zb,sprod)
    a0phigggg_xxxx(1,-1,-1,+1,-1) = a0phiggggpmmm(j3,j4,j1,j2,za,zb,sprod)
    a0phigggg_xxxx(1,-1,-1,-1,+1) = a0phiggggpmmm(j4,j1,j2,j3,za,zb,sprod)
    a0phigggg_xxxx(1,-1,-1,-1,-1) = a0phiggggmmmm(j1,j2,j3,j4,za,zb,sprod)
 
    a0phigggg_xxxx(2,-1,-1,+1,+1) = conjg(a0phigggg_xxxx(1,+1,+1,-1,-1))
    a0phigggg_xxxx(2,-1,+1,-1,+1) = conjg(a0phigggg_xxxx(1,+1,-1,+1,-1))
    a0phigggg_xxxx(2,-1,+1,+1,-1) = conjg(a0phigggg_xxxx(1,+1,-1,-1,+1))
    a0phigggg_xxxx(2,-1,+1,+1,+1) = conjg(a0phigggg_xxxx(1,+1,-1,-1,-1))
    a0phigggg_xxxx(2,+1,-1,-1,+1) = conjg(a0phigggg_xxxx(1,-1,+1,+1,-1))
    a0phigggg_xxxx(2,+1,-1,+1,-1) = conjg(a0phigggg_xxxx(1,-1,+1,-1,+1))
    a0phigggg_xxxx(2,+1,-1,+1,+1) = conjg(a0phigggg_xxxx(1,-1,+1,-1,-1))
    a0phigggg_xxxx(2,+1,+1,-1,-1) = conjg(a0phigggg_xxxx(1,-1,-1,+1,+1))
    a0phigggg_xxxx(2,+1,+1,-1,+1) = conjg(a0phigggg_xxxx(1,-1,-1,+1,-1))
    a0phigggg_xxxx(2,+1,+1,+1,-1) = conjg(a0phigggg_xxxx(1,-1,-1,-1,+1))
    a0phigggg_xxxx(2,+1,+1,+1,+1) = conjg(a0phigggg_xxxx(1,-1,-1,-1,-1))
 
    return
 

a0phiggggmmmm()

complex(dp) function modhiggsjj::a0phiggggmmmm ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod  )

private

Definition at line 3102 of file mod_HiggsJJ.F90.

    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: sprod(4,4),qsq
    complex(dp) :: A0phiggggmmmm
 
    qsq = sprod(j1,j2)+sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4)+sprod(j3,j4)
    a0phiggggmmmm=qsq**2/(zb(j1,j2)*zb(j2,j3)*zb(j3,j4)*zb(j4,j1))
 
    return
 

a0phiggggmmpp()

complex(dp) function modhiggsjj::a0phiggggmmpp ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod  )

private

Definition at line 3088 of file mod_HiggsJJ.F90.

    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: sprod(4,4)
    complex(dp) :: A0phiggggmmpp
 
    a0phiggggmmpp = za(j1,j2)**4 &
         /(za(j1,j2)*za(j2,j3)*za(j3,j4)*za(j4,j1))
 
    return
 

a0phiggggmpmp()

complex(dp) function modhiggsjj::a0phiggggmpmp ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod  )

private

Definition at line 3074 of file mod_HiggsJJ.F90.

    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: sprod(4,4)
    complex(dp) :: A0phiggggmpmp
 
    a0phiggggmpmp = za(j1,j3)**4 &
         /(za(j1,j2)*za(j2,j3)*za(j3,j4)*za(j4,j1))
 
    return
 

a0phiggggpmmm()

complex(dp) function modhiggsjj::a0phiggggpmmm ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod  )

private

Definition at line 3048 of file mod_HiggsJJ.F90.

    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: sprod(4,4)
    complex(dp) :: A0phiggggpmmm
 
    real(dp) :: s3
    complex(dp) :: zab2
 
    s3(j1,j2,j3)=sprod(j1,j2)+sprod(j2,j3)+sprod(j3,j1)
    zab2(j1,j2,j3,j4)=za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)
 
    a0phiggggpmmm = &
         +(zab2(j3,j2,j4,j1)*za(j2,j4))**2/(s3(j1,j2,j4)*sprod(j1,j2)*sprod(j1,j4)) &
         +(zab2(j4,j2,j3,j1)*za(j2,j3))**2/(s3(j1,j2,j3)*sprod(j1,j2)*sprod(j2,j3)) &
         +(zab2(j2,j3,j4,j1)*za(j3,j4))**2/(s3(j1,j3,j4)*sprod(j1,j4)*sprod(j3,j4)) &
         -za(j2,j4)/(za(j1,j2)*zb(j2,j3)*zb(j3,j4)*za(j4,j1)) &
         *(-sprod(j2,j3)*zab2(j2,j3,j4,j1)/zb(j4,j1) &
         -sprod(j3,j4)*zab2(j4,j2,j3,j1)/zb(j1,j2) &
         -s3(j2,j3,j4)*za(j2,j4))
 
    return
 

a0phiqbqgg_xxx()

complex(dp) function, dimension(1:2,-1:1,-1:1,-1:1) modhiggsjj::a0phiqbqgg_xxx ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod  )

private

Definition at line 2982 of file mod_HiggsJJ.F90.

    complex(dp) :: A0phiqbqgg_xxx(1:2,-1:1,-1:1,-1:1)
    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: sprod(4,4)
    real(dp) :: s3
    complex(dp) :: zab2
 
    s3(j1,j2,j3) = sprod(j1,j2) + sprod(j1,j3) + sprod(j2,j3)
    zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
 
    a0phiqbqgg_xxx = czero
 
    a0phiqbqgg_xxx(1,-1,+1,-1) = -za(j1,j4)**2*za(j2,j4)/(za(j1,j2)*za(j2,j3)*za(j3,j4))
    a0phiqbqgg_xxx(1,-1,-1,+1) = za(j1,j3)**3/(za(j1,j2)*za(j3,j4)*za(j4,j1))
    a0phiqbqgg_xxx(1,-1,-1,-1) =  -((za(j1,j3)*zab2(j4,j1,j3,j2)**2)/ &
         (s3(j1,j2,j3)*sprod(j1,j2)*zb(j2,j3))) - &
         ((1/sprod(j1,j2) + 1/sprod(j4,j1))*za(j4,j1)*zab2(j3,j1,j4,j2)**2)/ &
         (s3(j4,j1,j2)*zb(j2,j4)) + &
         zab2(j1,j3,j4,j2)**2/(za(j1,j2)*zb(j2,j3)*zb(j2,j4)*zb(j3,j4))
 
    a0phiqbqgg_xxx(2,-1,+1,-1) = -zb(j1,j3)*zb(j2,j3)**2/(zb(j1,j2)*zb(j3,j4)*zb(j4,j1))
    a0phiqbqgg_xxx(2,-1,-1,+1) = zb(j2,j4)**3/(zb(j1,j2)*zb(j2,j3)*zb(j3,j4))
    a0phiqbqgg_xxx(2,-1,+1,+1) = -(zab2(j1,j3,j4,j2)**2/(za(j1,j3)*za(j1,j4)*za(j3,j4)*zb(j1,j2))) - &
         ((1/sprod(j1,j2) + 1/sprod(j2,j3))*zab2(j1,j2,j3,j4)**2*zb(j2,j3))/ &
         (s3(j1,j2,j3)*za(j1,j3)) + &
         (zab2(j1,j2,j4,j3)**2*zb(j2,j4))/ &
         (s3(j4,j1,j2)*sprod(j1,j2)*za(j1,j4))
 
    a0phiqbqgg_xxx(1,+1,+1,-1) = conjg(a0phiqbqgg_xxx(2,-1,-1,+1))
    a0phiqbqgg_xxx(1,+1,-1,+1) = conjg(a0phiqbqgg_xxx(2,-1,+1,-1))
    a0phiqbqgg_xxx(1,+1,-1,-1) = conjg(a0phiqbqgg_xxx(2,-1,+1,+1))
 
    a0phiqbqgg_xxx(2,+1,+1,-1) = conjg(a0phiqbqgg_xxx(1,-1,-1,+1))
    a0phiqbqgg_xxx(2,+1,-1,+1) = conjg(a0phiqbqgg_xxx(1,-1,+1,-1))
    a0phiqbqgg_xxx(2,+1,+1,+1) = conjg(a0phiqbqgg_xxx(1,-1,-1,-1))
 
    return
 

a0phiqbqqbq_xx()

complex(dp) function, dimension(1:2,-1:1,-1:1) modhiggsjj::a0phiqbqqbq_xx ( integer  j1, integer  j2, integer  j3, integer  j4, complex(dp), dimension(4,4)  za, complex(dp), dimension(4,4)  zb, real(dp), dimension(4,4)  sprod  )

private

Definition at line 3023 of file mod_HiggsJJ.F90.

    complex(dp) :: A0phiqbqQBQ_xx(1:2,-1:1,-1:1)
    integer :: j1,j2,j3,j4
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: sprod(4,4)
 
    a0phiqbqqbq_xx = czero
 
    a0phiqbqqbq_xx(1,-1,+1) = za(j1,j4)**2/(za(j1,j2)*za(j3,j4))
    a0phiqbqqbq_xx(1,-1,-1) = za(j1,j3)**2/(za(j1,j2)*za(j4,j3))
 
    a0phiqbqqbq_xx(2,-1,+1) = zb(j2,j3)**2/(zb(j1,j2)*zb(j3,j4))
    a0phiqbqqbq_xx(2,-1,-1) = zb(j2,j4)**2/(zb(j1,j2)*zb(j4,j3))
 
    a0phiqbqqbq_xx(1,+1,-1) = conjg(a0phiqbqqbq_xx(2,-1,+1))
    a0phiqbqqbq_xx(1,+1,+1) = conjg(a0phiqbqqbq_xx(2,-1,-1))
 
    a0phiqbqqbq_xx(2,+1,-1) = conjg(a0phiqbqqbq_xx(1,-1,+1))
    a0phiqbqqbq_xx(2,+1,+1) = conjg(a0phiqbqqbq_xx(1,-1,-1))
 
    return
 

evalamp_sbfh_unsymm_sa()

subroutine, public modhiggsjj::evalamp_sbfh_unsymm_sa	(	real(dp), dimension(4,5), intent(in)	p,
		real(dp), dimension(-5:5,-5:5), intent(out)	res
	)

Definition at line 24 of file mod_HiggsJJ.F90.

    real(dp), intent(in) :: p(4,5)
    real(dp), intent(out) :: res(-5:5,-5:5)
    real(dp) :: sprod(4,4)
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: restmp, restmpid
    integer :: i, j
 
    res = zero
 
    call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
 
    !-- gg -> gg
    call me2_ggggh(1,2,3,4,za,zb,sprod,restmp)
    restmp = restmp * avegg * symmfac
    res(0,0) = res(0,0) + restmp
 
    !-- gg -> qqb
    call me2_qbqggh(4,3,1,2,za,zb,sprod,restmp)
    restmp = restmp * avegg
    res(0,0) = res(0,0) + restmp * nf
 
    !-- gq -> gq
    call me2_qbqggh(2,4,1,3,za,zb,sprod,restmp)
    restmp = restmp * aveqg
    do i = 1,5
       res(0,i) = res(0,i) + restmp
       res(0,-i) = res(0,-i) + restmp
    enddo
    call me2_qbqggh(1,4,2,3,za,zb,sprod,restmp)
    restmp = restmp * aveqg
    do i = 1,5
       res(i,0) = res(i,0) + restmp
       res(-i,0) = res(-i,0) + restmp
    enddo
 
    !-- qqb -> gg
    call me2_qbqggh(1,2,3,4,za,zb,sprod,restmp)
    restmp = restmp * aveqq * symmfac
    do i = 1,5
       res(i,-i) = res(i,-i) + restmp
       res(-i,i) = res(-i,i) + restmp
    enddo
 
    !-- qqb -> qqb
    call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
    restmp = restmpid * aveqq
    do i = 1,5
       res(i,-i) =res(i,-i) + restmp
    enddo
    call me2_qbqqbq(2,1,4,3,za,zb,sprod,restmp,restmpid)
    restmp = restmpid * aveqq
    do i = 1,5
       res(-i,i) =res(-i,i) + restmp
    enddo
 
    !-- qqb -> rrb
    call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
    restmp = restmp * aveqq
    do i = 1,5
       res(i,-i) =res(i,-i) + restmp * (nf-1.0_dp)
       res(-i,i) = res(-i,i) + restmp * (nf-1.0_dp)
    enddo
 
    !-- qrb -> qrb
    call me2_qbqqbq(1,3,4,2,za,zb,sprod,restmp,restmpid)
    restmp = restmp * aveqq
    do i = 1,5
       do j = 1,5
          if (i.ne.j) res(i,-j) = res(i,-j) + restmp
       enddo
    enddo
    call me2_qbqqbq(2,3,4,1,za,zb,sprod,restmp,restmpid)
    restmp = restmp * aveqq
    do i = 1,5
       do j = 1,5
          if (i.ne.j) res(-j,i) = res(-j,i) + restmp
       enddo
    enddo
 
    !-- qq -> qq
    call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
    restmp = restmpid * aveqq * symmfac
    do i = 1,5
       res(i,i) = res(i,i) + restmp
       res(-i,-i) = res(-i,-i) + restmp
    enddo
 
    !-- qr -> qr
    call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
    restmp = restmp * aveqq
    do i = 1,5
       do j = 1,5
          if (i.ne.j) then
             res(i,j) = res(i,j) + restmp
             res(-i,-j) = res(-i,-j) + restmp
          endif
       enddo
    enddo
 
    res(:,:) = res(:,:) * (2d0/3d0*alphas**2)**2
    return
 

evalamp_sbfh_unsymm_sa_select()

subroutine, public modhiggsjj::evalamp_sbfh_unsymm_sa_select	(	real(dp), dimension(4,5), intent(in)	p,
		integer, intent(in)	iSel,
		integer, intent(in)	jSel,
		integer, intent(in)	flav_tag,
		integer	iflip,
		real(dp), dimension(-5:5,-5:5), intent(out)	res
	)

Definition at line 135 of file mod_HiggsJJ.F90.

    real(dp), intent(in) :: p(4,5)
    integer, intent(in) :: iSel,jSel,flav_tag
    real(dp), intent(out) :: res(-5:5,-5:5)
    real(dp) :: sprod(4,4)
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: restmp, restmpid
    integer :: i, j,iflip
 
    res = zero
 
    call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
    iflip = 1
 
    !-- gg -> gg
    if( isel.eq.pdfglu_ .and. jsel.eq.pdfglu_ .and. flav_tag.eq.1 ) then
        call me2_ggggh(1,2,3,4,za,zb,sprod,restmp)
        restmp = restmp * avegg * symmfac
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
 
    !-- gg -> qqb
    if( isel.eq.pdfglu_ .and. jsel.eq.pdfglu_ .and. flav_tag.eq.2 ) then
        call me2_qbqggh(4,3,1,2,za,zb,sprod,restmp)
        restmp = restmp * avegg
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp * nf
        return
    endif
 
 
 
    !-- gq -> gq
    if( isel.eq.pdfglu_ .and. jsel.ne.0 ) then
        call me2_qbqggh(2,4,1,3,za,zb,sprod,restmp)
        restmp = restmp * aveqg
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
    if( jsel.eq.pdfglu_ .and. isel.ne.0 ) then
        iflip = 2
        call me2_qbqggh(1,4,2,3,za,zb,sprod,restmp)
        restmp = restmp * aveqg
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
 
 
    !-- qqb -> gg
    if( isel.ne.0 .and. jsel.eq.-isel .and. flav_tag.eq.1 ) then
        call me2_qbqggh(1,2,3,4,za,zb,sprod,restmp)
        restmp = restmp * aveqq * symmfac
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
 
    !-- qqb -> qqb
    if( isel.gt.0 .and. jsel.eq.-isel .and. flav_tag.eq.2 ) then
        call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
        restmp = restmpid * aveqq
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
    if( isel.lt.0 .and. jsel.eq.-isel .and. flav_tag.eq.2 ) then
        iflip = 2
        call me2_qbqqbq(2,1,4,3,za,zb,sprod,restmp,restmpid)
        restmp = restmpid * aveqq
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
 
    !-- qqb -> rrb
    if( isel.ne.0 .and. jsel.eq.-isel .and. flav_tag.eq.3 ) then
        call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
        restmp = restmp * aveqq
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp * (nf-1.0_dp)
        return
    endif
 
 
 
    !-- qrb -> qrb
    if( isel.gt.0 .and. jsel.lt.0 .and. isel.ne.jsel ) then
        call me2_qbqqbq(1,3,4,2,za,zb,sprod,restmp,restmpid)
        restmp = restmp * aveqq
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
    if( isel.lt.0 .and. jsel.gt.0 .and. isel.ne.jsel ) then
        iflip = 2
        call me2_qbqqbq(2,3,4,1,za,zb,sprod,restmp,restmpid)
        restmp = restmp * aveqq
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
 
 
 
    !-- qq -> qq
    if( isel.ne.0 .and. isel.eq.jsel ) then
        call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
        restmp = restmpid * aveqq * symmfac
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
 
 
    !-- qr -> qr
    if( isel.gt.0 .and. jsel.gt.0 .and. isel.ne.jsel ) then
        call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
        restmp = restmp * aveqq
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
    if( isel.lt.0 .and. jsel.lt.0 .and. isel.ne.jsel ) then
        call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
        restmp = restmp * aveqq
        restmp = restmp * (2d0/3d0*alphas**2)**2
        res(isel,jsel) = restmp
        return
    endif
 
 
 
  return
 

evalamp_sbfh_unsymm_sa_select_exact()

subroutine, public modhiggsjj::evalamp_sbfh_unsymm_sa_select_exact	(	real(dp), dimension(4,5), intent(in)	p,
		integer, intent(in)	iSel,
		integer, intent(in)	jSel,
		integer, intent(in)	rSel,
		integer, intent(in)	sSel,
		real(dp), intent(out)	res
	)

Definition at line 286 of file mod_HiggsJJ.F90.

    real(dp), intent(in) :: p(4,5)
    integer, intent(in) :: iSel,jSel,rSel,sSel!,flavor_tag ! flavor_tag for TEST
    real(dp), intent(out) :: res
    real(dp) :: sprod(4,4)
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: restmp, restmpid
    integer :: i, j
    integer :: j1,j2,k1,k2
    logical :: isGGGG,isQQGG,isQQQQ,isQQQQ_idQ
 
    restmp=0.0_dp
    restmpid=0.0_dp
    res=0.0_dp
    isgggg=.false.
    isqqgg=.false.
    isqqqq=.false.
    isqqqq_idq=.false.
    j1 = 0
    j2 = 0
    k1 = 0
    k2 = 0
    if( (abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_) .or. (abs(rsel).eq.pdftop_ .or. abs(ssel).eq.pdftop_) ) return
 
    !print *, "Begin EvalAmp_SBFH_UnSymm_SA_Select_exact"
    !print *, "iSel: ",iSel,", jSel: ",jSel," rSel: ",rSel," sSel: ",sSel
 
    if(isel.eq.pdfglu_) then ! 1==g
       if(jsel.eq.pdfglu_) then ! gg
          if(rsel.eq.pdfglu_ .and. ssel .eq. pdfglu_) then ! gg->gg
             isgggg=.true.
            j1=1
          j2=2
             k1=3
             k2=4
 
          elseif(ssel.eq.(-rsel)) then ! gg->qqb/qbq
             isqqgg=.true.
            if(rsel.gt.0) then ! gg->qqb
               j1=4
                j2=3
            else ! gg->qbq
               j1=3
                j2=4
             endif
               k1=1
              k2=2
          endif
 
       elseif(abs(jsel).eq.pdfup_ .or. abs(jsel).eq.pdfchm_ .or. abs(jsel).eq.pdfdn_ .or. abs(jsel).eq.pdfstr_ .or. abs(jsel).eq.pdfbot_) then ! gq/gqb
           isqqgg=.true.
          k1=1
          if(jsel.gt.0) then ! gq
             j1=2
             if(rsel.eq.isel .and. ssel.eq.jsel) then ! gq->gq
                k2=3
                  j2=4
             elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! gq->qg
                 j2=3
                  k2=4
             endif
          else ! gqb
             j2=2
            if(rsel.eq.isel .and. ssel.eq.jsel) then ! gqb->gqb
                k2=3
                  j1=4
             elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! gqb->qbg
                 j1=3
                  k2=4
             endif
          endif
       endif
    elseif(abs(isel).eq.pdfup_ .or. abs(isel).eq.pdfchm_ .or. abs(isel).eq.pdfdn_ .or. abs(isel).eq.pdfstr_ .or. abs(isel).eq.pdfbot_) then ! q/qb?->?
       if(jsel.eq.pdfglu_) then ! qg->?
          isqqgg=.true.
          k1=2
          if(isel.gt.0) then
             j1=1
             if(rsel.eq.isel .and. ssel.eq.jsel) then ! qg->qg
                k2=4
                j2=3
             elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! qg->gq
                j2=4
                k2=3
             endif
          else ! qbg
             j2=1
             if(rsel.eq.isel .and. ssel.eq.jsel) then ! qbg->qbg
                k2=4
                  j1=3
             elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! qbg->gqb
                j1=4
                 k2=3
             endif
          endif
 
       elseif(abs(jsel).eq.pdfup_ .or. abs(jsel).eq.pdfchm_ .or. abs(jsel).eq.pdfdn_ .or. abs(jsel).eq.pdfstr_ .or. abs(jsel).eq.pdfbot_) then ! qq/qqb/qbq/qbqb->?
          if(isel.gt.0) then ! qq'/qqb'->?
             j1=1
          else ! qbq'/qbqb'->?
             j2=1
          endif
 
          if(jsel.eq.(-isel)) then ! qqb/qbq->?
             if(jsel.gt.0) then ! qbq->?
                j1=2
             else ! qqb->?
                j2=2
             endif
 
             if(rsel.eq.pdfglu_ .and. ssel .eq. pdfglu_) then ! qqb/qbq->gg
                isqqgg=.true.
                k1=3
                k2=4
             elseif(ssel.eq.(-rsel)) then ! qqb->q'qb'/qb'q'
                isqqqq=.true.
                if(rsel.gt.0) then ! ->q'qb'
                   if(rsel .eq. isel .or. rsel.eq.jsel) isqqqq_idq=.true.
                   k1=4
                   k2=3
                else ! ->qb'q'
                   if(rsel .eq. isel .or. rsel.eq.jsel) isqqqq_idq=.true.
                   k1=3
                   k2=4
                endif
             endif
 
          elseif(jsel.eq.isel) then ! qq/qbqb->?
             isqqqq=.true.
             if(jsel.gt.0) then ! qq
                k1=2
             else ! qbqb
                k2=2
             endif
             if(rsel.gt.0 .and. rsel.eq.ssel .and. rsel.eq.isel) then ! qq->qq
                isqqqq_idq=.true.
                j2=3
                k2=4
             elseif(rsel.lt.0 .and. rsel.eq.ssel .and. rsel.eq.isel) then ! qbqb->qbqb
                isqqqq_idq=.true.
                j1=3
                k1=4
             endif
 
          elseif(jsel.eq.sign(jsel,isel)) then ! qq'/qbqb'->?
             isqqqq=.true.
             if(jsel.gt.0) then ! qq'
                k1=2 ! j1=1
             else ! qbqb'
                k2=2 ! j2=1
             endif
             if(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.isel) then ! ->qq'
                j2=3
                k2=4
             elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.isel) then ! ->qbqb'
                j1=3
                k1=4
             elseif(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.jsel) then ! ->q'q
                j2=4
                k2=3
             elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.jsel) then ! ->qb'qb
                j1=4
                k1=3
             endif
 
          elseif(jsel.eq.(-sign(jsel,isel))) then ! qqb'/qbq'->?
             isqqqq=.true.
             if(jsel.gt.0) then ! qbq'
                k1=2 ! j2=1
             else ! qqb'
                k2=2 ! j1=1
             endif
             if(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.isel) then ! qqb'->qqb'
                j2=3
                k1=4
             elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.isel) then ! qbq'->qbq'
                j1=3
                k2=4
             elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.jsel) then ! qqb'->qb'q
                j2=4
                k1=3
             elseif(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.jsel) then ! qbq'->q'qb
                j1=4
                k2=3
             endif
          endif
 
       endif
    endif
 
    if(j1.eq.0 .or. j2.eq.0 .or. k1.eq.0 .or. k2.eq.0) then
       print *,"Could not recognize the incoming flavors!"
       print *,"iSel: ",isel,", jSel: ",jsel," rSel: ",rsel," sSel: ",ssel
       print *,"j1: ",j1,", j2: ",j2," k1: ",k1," k2: ",k2
        if(isgggg) then
            print *,"is gggg"
        endif
        if(isqqgg) then
            print *,"is qbqgg"
        endif
        if(isqqqq) then
            print *,"is qbqQBQ"
        endif
        if(isqqqq_idq) then
            print *,"is qbqQBQ id"
        endif
    endif
 
    call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
 
    if(isgggg) then
       call me2_ggggh(j1,j2,k1,k2,za,zb,sprod,restmp)
       restmp = restmp * avegg * symmfac
    elseif(isqqqq) then
        call me2_qbqqbq(j1,j2,k1,k2,za,zb,sprod,restmp,restmpid)
        if(isqqqq_idq) then
           restmp = restmpid * aveqq
           if(isel.eq.jsel) restmp = restmp * symmfac
        else
           restmp = restmp * aveqq
           if(isel .eq. (-jsel) .and. rsel .eq. (-ssel) .and. abs(isel) .ne. abs(rsel)) restmp = restmp * (nf-1.0_dp)
        endif
    elseif(isqqgg) then
        call me2_qbqggh(j1,j2,k1,k2,za,zb,sprod,restmp)
        if( isel.eq.pdfglu_ .and. jsel.eq.pdfglu_) then
           restmp = restmp * avegg * nf
        elseif( isel.ne.pdfglu_ .and. jsel.ne.pdfglu_) then
           restmp = restmp * aveqq * symmfac
        else
           restmp = restmp * aveqg
        endif
    endif
 
   restmp = restmp * (2d0/3d0*alphas**2)**2
    res = restmp
    return

evalamp_wbfh_unsymm_sa()

subroutine, public modhiggsjj::evalamp_wbfh_unsymm_sa	(	real(dp), dimension(4,5), intent(in)	p,
		real(dp), dimension(-5:5,-5:5), intent(out)	res
	)

Definition at line 526 of file mod_HiggsJJ.F90.

    real(dp), intent(in) :: p(4,5)
    real(dp), intent(out) :: res(-5:5,-5:5)
    complex(dp) :: amp_z(-1:1,-1:1), amp_z_b(-1:1,-1:1)
    complex(dp) :: amp_w(-1:1,-1:1)
    real(dp) :: sprod(4,4)
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: Lu, Ru
    real(dp) :: Ld, Rd
    real(dp) :: couplz
    real(dp) :: couplw
    real(dp) :: restmp=0d0
    integer :: i, j, j1, j2, iflip, pdfindex(2)
 
    lu = al_qup**2
    ru = ar_qup**2
    ld = al_qdn**2
    rd = ar_qdn**2
    couplz = gwsq * xw/twosc**2
    couplw = gwsq/two
 
    res = zero
 
    call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
 
    !-- qq->qq, up
    amp_z = a0_vv_4f(4,1,3,2,za,zb,sprod,m_z,ga_z)
    amp_z_b = -a0_vv_4f(3,1,4,2,za,zb,sprod,m_z,ga_z)
 
! !     adding contract terms
!     iprop12 = sprod(4,1) - mv**2 + ci * mv * ga_v
!     iprop34 = sprod(3,2) - mv**2 + ci * mv * ga_v
!
!     amp_z(-1,-1) = amp_z(-1,-1) + amp_z(-1,-1)*iprop12/(aL_QUp*couplz)*ehz_L_U  &
!                                 + amp_z(-1,-1)*iprop34/(aL_QUp*couplz)*ehz_L_U  &
!                                 + amp_z(-1,-1)*iprop12/(aL_QUp*couplz)*ehz_L_U  &
!                                               *iprop34/(aL_QUp*couplz)*ehz_L_U
!
!     amp_z(+1,-1) = amp_z(+1,-1) + amp_z(+1,-1)*iprop12/(aR_QUp*couplz)*ehz_R_U  &
!                                 + amp_z(+1,-1)*iprop34/(aL_QUp*couplz)*ehz_L_U  &
!                                 + amp_z(+1,-1)*iprop12/(aR_QUp*couplz)*ehz_R_U  &
!                                               *iprop34/(aL_QUp*couplz)*ehz_L_U
!
!     amp_z(-1,+1) = amp_z(-1,+1) + amp_z(-1,+1)*iprop12/(aL_QUp*couplz)*ehz_L_U  &
!                                 + amp_z(-1,+1)*iprop34/(aR_QUp*couplz)*ehz_R_U  &
!                                 + amp_z(-1,+1)*iprop12/(aL_QUp*couplz)*ehz_L_U  &
!                                               *iprop34/(aR_QUp*couplz)*ehz_R_U
!
!     amp_z(+1,+1) = amp_z(+1,+1) + amp_z(+1,+1)*iprop12/(aR_QUp*couplz)*ehz_R_U  &
!                                 + amp_z(+1,+1)*iprop34/(aR_QUp*couplz)*ehz_R_U  &
!                                 + amp_z(+1,+1)*iprop12/(aR_QUp*couplz)*ehz_R_U  &
!                                               *iprop34/(aR_QUp*couplz)*ehz_R_U
 
 
!     Lu = aL_QUp + iprop12/couplz*ehz_L_U
!     Ru = aR_QUp + iprop12/couplz*ehz_R_U
!     amp_z(-1,-1) = amp_z(-1,-1)*Lu*Lu
!     amp_z(-1,+1) = amp_z(-1,+1)*Lu*Ru
!     amp_z(+1,-1) = amp_z(+1,-1)*Ru*Lu
!     amp_z(+1,+1) = amp_z(+1,+1)*Ru*Ru
!
!     amp_z_b(-1,-1) = amp_z_b(-1,-1)*Lu*Lu
!     amp_z_b(-1,+1) = amp_z_b(-1,+1)*Lu*Ru
!     amp_z_b(+1,-1) = amp_z_b(+1,-1)*Ru*Lu
!     amp_z_b(+1,+1) = amp_z_b(+1,+1)*Ru*Ru
 
    restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * lu**2 + &
         (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * lu * ru + &
         (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * lu * ru + &
         (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * ru**2) * xn**2
 
    restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * lu**2 + &
                       two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * ru**2) * xn
 
    restmp = restmp * symmfac * aveqq * couplz**2
 
    res(pdfup_,pdfup_) = restmp
    res(pdfchm_,pdfchm_) = restmp
 
 
 
 
 
    !-- qq->qq, down
    restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * ld**2 + &
         (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * ld * rd + &
         (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * ld * rd + &
         (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * rd**2) * xn**2
 
    restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * ld**2 + &
            two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * rd**2) * xn
 
    restmp = restmp * symmfac * aveqq * couplz**2
 
    res(pdfdn_,pdfdn_) = restmp
    res(pdfstr_,pdfstr_) = restmp
    res(pdfbot_,pdfbot_) = restmp * tagbot
 
    !-- qbqb->qbqb, aup
    amp_z = a0_vv_4f(1,4,2,3,za,zb,sprod,m_z,ga_z)
    amp_z_b = -a0_vv_4f(1,3,2,4,za,zb,sprod,m_z,ga_z)
 
    restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * lu**2 + &
         (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * lu * ru + &
         (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * lu * ru + &
         (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * ru**2) * xn**2
 
    restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * lu**2 + &
            two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * ru**2) * xn
 
    restmp = restmp * symmfac * aveqq * couplz**2
 
    res(pdfaup_,pdfaup_) = restmp
    res(pdfachm_,pdfachm_) = restmp
 
    !-- qbqb->qbqb, adn
    restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * ld**2 + &
         (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * ld * rd + &
         (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * ld * rd + &
         (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * rd**2) * xn**2
 
    restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * ld**2 + &
            two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * rd**2) * xn
 
    restmp = restmp * symmfac * aveqq * couplz**2
 
    res(pdfadn_,pdfadn_) = restmp
    res(pdfastr_,pdfastr_) = restmp
    res(pdfabot_,pdfabot_) = restmp * tagbot
 
    !-- W/Z interference
    j1 = 1
    j2 = 2
    do iflip = 1, 2
       !-- ud -> ud
       amp_z =  a0_vv_4f(4,j2,3,j1,za,zb,sprod,m_z,ga_z)
       amp_w = -a0_vv_4f(4,j1,3,j2,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
 
       restmp = ((abs(amp_z(-1,-1))**2) * ld * lu + &
                 (abs(amp_z(-1,+1))**2) * ld * ru + &
                 (abs(amp_z(+1,-1))**2) * rd * lu + &
                 (abs(amp_z(+1,+1))**2) * rd * ru) * couplz**2 * xn**2
 
       restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * xn**2
 
       restmp = restmp + two * real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) * &
            al_qup * al_qdn * couplz * couplw * xn
 
       restmp = restmp * aveqq
 
       pdfindex = flip(iflip,pdfup_,pdfdn_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfchm_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       !-- ubdb -> ubdb
       amp_z = a0_vv_4f(j2,4,j1,3,za,zb,sprod,m_z,ga_z)
       amp_w = -a0_vv_4f(j1,4,j2,3,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
 
       restmp = ((abs(amp_z(-1,-1))**2) * ld * lu + &
            (abs(amp_z(-1,+1))**2) * ld * ru + &
            (abs(amp_z(+1,-1))**2) * rd * lu + &
            (abs(amp_z(+1,+1))**2) * rd * ru) * couplz**2 * xn**2
 
       restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * xn**2
 
       restmp = restmp + two * real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) * &
            al_qup * al_qdn * couplz * couplw * xn
 
       restmp = restmp * aveqq
 
       pdfindex = flip(iflip,pdfaup_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfachm_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       j1 = 2
       j2 = 1
    enddo
 
    !-- non-interfering diagrams below
    j1 = 1
    j2 = 2
    do iflip = 1,2
 
       !-- qqb processes
 
       !--uub -> uub // ddb
       amp_z = a0_vv_4f(3,j1,j2,4,za,zb,sprod,m_z,ga_z)
       amp_w = a0_vv_4f(3,j1,j2,4,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
 
       restmp = ((abs(amp_z(-1,-1))**2) * lu * lu + &
            (abs(amp_z(-1,+1))**2) * lu * ru + &
            (abs(amp_z(+1,-1))**2) * ru * lu + &
            (abs(amp_z(+1,+1))**2) * ru * ru) * couplz**2 * spinavg * tag1
       restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
 
       pdfindex = flip(iflip,pdfup_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfchm_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfup_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfchm_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       !--udb -> udb
       restmp = ((abs(amp_z(-1,-1))**2) * lu * ld + &
            (abs(amp_z(-1,+1))**2) * lu * rd + &
            (abs(amp_z(+1,-1))**2) * ru * ld + &
            (abs(amp_z(+1,+1))**2) * ru * rd) * couplz**2 * spinavg
 
       pdfindex = flip(iflip,pdfup_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfchm_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfup_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfchm_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfup_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfchm_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       !--dub -> dub
       restmp = ((abs(amp_z(-1,-1))**2) * ld * lu + &
            (abs(amp_z(-1,+1))**2) * ld * ru + &
            (abs(amp_z(+1,-1))**2) * rd * lu + &
            (abs(amp_z(+1,+1))**2) * rd * ru) * couplz**2 * spinavg
 
       pdfindex = flip(iflip,pdfdn_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfstr_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfdn_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfstr_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfbot_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfbot_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       !--ddb -> uub/ddb
       restmp = ((abs(amp_z(-1,-1))**2) * ld * ld + &
            (abs(amp_z(-1,+1))**2) * ld * rd + &
            (abs(amp_z(+1,-1))**2) * rd * ld + &
            (abs(amp_z(+1,+1))**2) * rd * rd) * couplz**2 * spinavg * tag1
 
       pdfindex = flip(iflip,pdfbot_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfdn_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfstr_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfbot_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfbot_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
 
       pdfindex = flip(iflip,pdfdn_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfstr_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfdn_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfstr_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       !-- non-symmetric qq processes
       amp_z = a0_vv_4f(3,j1,4,j2,za,zb,sprod,m_z,ga_z)
       amp_w = a0_vv_4f(3,j2,4,j1,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
 
       !--uc -> uc
       restmp = ((abs(amp_z(-1,-1))**2) * lu * lu + &
            (abs(amp_z(-1,+1))**2) * lu * ru + &
            (abs(amp_z(+1,-1))**2) * ru * lu + &
            (abs(amp_z(+1,+1))**2) * ru * ru) * couplz**2 * spinavg
 
       pdfindex = flip(iflip,pdfup_,pdfchm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       !--us -> us/cd
       restmp = ((abs(amp_z(-1,-1))**2) * lu * ld + &
            (abs(amp_z(-1,+1))**2) * lu * rd + &
            (abs(amp_z(+1,-1))**2) * ru * ld + &
            (abs(amp_z(+1,+1))**2) * ru * rd) * couplz**2 * spinavg * tag1
 
       pdfindex = flip(iflip,pdfup_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfchm_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
 
       pdfindex = flip(iflip,pdfup_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfchm_,pdfdn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       !--ds -> ds
       restmp = ((abs(amp_z(-1,-1))**2) * ld * ld + &
            (abs(amp_z(-1,+1))**2) * ld * rd + &
            (abs(amp_z(+1,-1))**2) * rd * ld + &
            (abs(amp_z(+1,+1))**2) * rd * rd) * couplz**2 * spinavg
 
       pdfindex = flip(iflip,pdfdn_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfdn_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfstr_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       !-- qbqb processes
       amp_z = a0_vv_4f(j1,3,j2,4,za,zb,sprod,m_z,ga_z)
       amp_w = a0_vv_4f(j1,4,j2,3,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
 
       !--ubcb -> ubcb
       restmp = ((abs(amp_z(-1,-1))**2) * lu * lu + &
            (abs(amp_z(-1,+1))**2) * lu * ru + &
            (abs(amp_z(+1,-1))**2) * ru * lu + &
            (abs(amp_z(+1,+1))**2) * ru * ru) * couplz**2 * spinavg
 
       pdfindex = flip(iflip,pdfaup_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       !--ubsb -> ubsb//cbdb
       restmp = ((abs(amp_z(-1,-1))**2) * lu * ld + &
            (abs(amp_z(-1,+1))**2) * lu * rd + &
            (abs(amp_z(+1,-1))**2) * ru * ld + &
            (abs(amp_z(+1,+1))**2) * ru * rd) * couplz**2 * spinavg * tag1
 
       pdfindex = flip(iflip,pdfaup_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfachm_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
 
       pdfindex = flip(iflip,pdfaup_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfachm_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       !--dbsb -> dbsb
       restmp = ((abs(amp_z(-1,-1))**2) * ld * ld + &
            (abs(amp_z(-1,+1))**2) * ld * rd + &
            (abs(amp_z(+1,-1))**2) * rd * ld + &
            (abs(amp_z(+1,+1))**2) * rd * rd) * couplz**2 * spinavg
 
       pdfindex = flip(iflip,pdfadn_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfadn_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       pdfindex = flip(iflip,pdfastr_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       j1 = 2
       j2 = 1
 
    enddo
 
    return
 

evalamp_wbfh_unsymm_sa_select()

subroutine, public modhiggsjj::evalamp_wbfh_unsymm_sa_select	(	real(dp), dimension(4,5), intent(in)	p,
		integer, intent(in)	iSel,
		integer, intent(in)	jSel,
		logical, intent(in)	zz_fusion,
		integer	iflip,
		real(dp), dimension(-5:5,-5:5), intent(out)	res
	)

Definition at line 898 of file mod_HiggsJJ.F90.

    implicit none
    real(dp), intent(in) :: p(4,5)
    real(dp), intent(out) :: res(-5:5,-5:5)
    logical, intent(in) :: zz_fusion
    integer, intent(in) :: iSel,jSel
    complex(dp) :: amp_z(-1:1,-1:1), amp_z_b(-1:1,-1:1)
    complex(dp) :: amp_w(-1:1,-1:1)
    real(dp) :: sprod(4,4)
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: restmp
    integer :: i, j, j1, j2, iflip, pdfindex(2)
 
    call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
    iflip = 1
 
    !-- qq->qq, up
    if( (isel.eq.pdfup_ .and. jsel.eq.pdfup_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfchm_) ) then
 
       amp_z = a0_zz_4f(4,1,3,2,za,zb,sprod,isel,jsel)
       amp_z_b = -a0_zz_4f(3,1,4,2,za,zb,sprod,isel,jsel)
 
       restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
 
       restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
            two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
 
       restmp = restmp * symmfac * aveqq
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       res(pdfup_,pdfup_) = restmp
       res(pdfchm_,pdfchm_) = restmp
 
       return
    endif
 
 
   !-- qq->qq, down
    if( (isel.eq.pdfdn_ .and. jsel.eq.pdfdn_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfstr_) .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfbot_) ) then
       amp_z = a0_zz_4f(4,1,3,2,za,zb,sprod,isel,jsel)
       amp_z_b = -a0_zz_4f(3,1,4,2,za,zb,sprod,isel,jsel)
 
       restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
 
       restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
            two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
 
       restmp = restmp * symmfac * aveqq
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       res(pdfdn_,pdfdn_) = restmp
       res(pdfstr_,pdfstr_) = restmp
       res(pdfbot_,pdfbot_) = restmp * tagbot
       return
    endif
 
 
    !-- qbqb->qbqb, aup
    if( (isel.eq.pdfaup_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfachm_) ) then
       amp_z = a0_zz_4f(1,4,2,3,za,zb,sprod,isel,jsel)
       amp_z_b = -a0_zz_4f(1,3,2,4,za,zb,sprod,isel,jsel)
 
       restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
 
       restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
            two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
 
       restmp = restmp * symmfac * aveqq
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       res(pdfaup_,pdfaup_) = restmp
       res(pdfachm_,pdfachm_) = restmp
       return
    endif
 
 
    !-- qbqb->qbqb, adn
    if( (isel.eq.pdfadn_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfastr_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfabot_ .and. jsel.eq.pdfabot_) ) then
       amp_z = a0_zz_4f(1,4,2,3,za,zb,sprod,isel,jsel)
       amp_z_b = -a0_zz_4f(1,3,2,4,za,zb,sprod,isel,jsel)
       restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
 
       restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
            two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
 
       restmp = restmp * symmfac * aveqq
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       res(pdfadn_,pdfadn_) = restmp
       res(pdfastr_,pdfastr_) = restmp
       res(pdfabot_,pdfabot_) = restmp * tagbot
       return
    endif
 
 
 
    ! NEED TO ADJSUT J1,J2 in W AMPS HERE
    !-- W/Z interference
    j1 = 1
    j2 = 2
    iflip = 1
 
    !-- ud -> ud
    if( (isel.eq.pdfup_ .and. jsel.eq.pdfdn_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfstr_) ) then
       amp_z = a0_zz_4f(4,j2,3,j1,za,zb,sprod,jsel,isel)
       amp_w = -a0_ww_4f(4,j1,3,j2,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.true.)
       if( zz_fusion ) then
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * xn**2
       else
          restmp= ( &
                   abs(amp_w(-1,-1))**2 + &
                   abs(amp_w(-1,+1))**2 + &
                   abs(amp_w(+1,-1))**2 + &
                   abs(amp_w(+1,+1))**2   &
                  ) * xn**2
          restmp = restmp + two * ( &
             real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
             real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
             real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
             real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp)   &
             ) * xn
 
       endif
 
       restmp = restmp * aveqq
 
       pdfindex = flip(iflip,pdfup_,pdfdn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !-- ubdb -> ubdb
    if( (isel.eq.pdfaup_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfastr_) ) then
       amp_z = a0_zz_4f(j2,4,j1,3,za,zb,sprod,jsel,isel)
       amp_w = -a0_ww_4f(j1,4,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
 
       if( zz_fusion ) then
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * xn**2
       else
          restmp= ( &
                   abs(amp_w(-1,-1))**2 + &
                   abs(amp_w(-1,+1))**2 + &
                   abs(amp_w(+1,-1))**2 + &
                   abs(amp_w(+1,+1))**2   &
                  ) * xn**2
          restmp = restmp + two * ( &
             real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
             real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
             real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
             real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp)   &
             ) * xn
 
       endif
 
       restmp = restmp * aveqq
 
       pdfindex = flip(iflip,pdfaup_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
       pdfindex = flip(iflip,pdfachm_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    ! NEED TO ADJSUT J1,J2 in W AMPS HERE
    !-- W/Z interference
    j1 = 2
    j2 = 1
    iflip = 2
 
    !-- ud -> ud
    if( (isel.eq.pdfdn_ .and. jsel.eq.pdfup_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfchm_) ) then
       amp_z = a0_zz_4f(4,j2,3,j1,za,zb,sprod,isel,jsel)
       amp_w = -a0_ww_4f(4,j1,3,j2,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.true.)
       if( zz_fusion ) then
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * xn**2
       else
          restmp= ( &
                   abs(amp_w(-1,-1))**2 + &
                   abs(amp_w(-1,+1))**2 + &
                   abs(amp_w(+1,-1))**2 + &
                   abs(amp_w(+1,+1))**2   &
                  ) * xn**2
          restmp = restmp + two * ( &
             real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
             real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
             real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
             real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp)   &
             ) * xn
 
       endif
 
       restmp = restmp * aveqq
 
       pdfindex = flip(iflip,pdfup_,pdfdn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !-- ubdb -> ubdb
    if( (isel.eq.pdfadn_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfastr_ .and. jsel.eq.pdfachm_) ) then
       amp_z = a0_zz_4f(j2,4,j1,3,za,zb,sprod,isel,jsel)
       amp_w = -a0_ww_4f(j1,4,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
 
       if( zz_fusion ) then
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * xn**2
       else
          restmp= ( &
                   abs(amp_w(-1,-1))**2 + &
                   abs(amp_w(-1,+1))**2 + &
                   abs(amp_w(+1,-1))**2 + &
                   abs(amp_w(+1,+1))**2   &
                  ) * xn**2
          restmp = restmp + two * ( &
             real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
             real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
             real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
             real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp)   &
             ) * xn
       endif
 
       restmp = restmp * aveqq
 
       pdfindex = flip(iflip,pdfaup_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfachm_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !-- non-interfering diagrams below
    j1 = 1
    j2 = 2
    iflip = 1
 
    !-- qqb processes
 
    !--uub -> uub // ddb
    if( (isel.eq.pdfup_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfachm_) .or. (isel.eq.pdfup_ .and. jsel.eq.pdfachm_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfaup_) ) then
 
       if( zz_fusion ) then
          amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.true.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfup_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfup_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !--udb -> udb
    if( (isel.eq.pdfup_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfup_ .and. jsel.eq.pdfastr_) .or. &
         (isel.eq.pdfchm_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfup_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfabot_) ) then
       amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfup_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfup_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfup_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfchm_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    !--dub -> dub
    if( (isel.eq.pdfdn_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfachm_)  .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfachm_)  .or. &
         (isel.eq.pdfstr_ .and. jsel.eq.pdfaup_)  .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfaup_)  .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfachm_) ) then
       amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfdn_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfdn_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfbot_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfbot_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    !--ddb -> uub/ddb
    if( (isel.eq.pdfbot_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfabot_) .or.&
         (isel.eq.pdfbot_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfastr_) ) then
       amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
 
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg * tag1
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfbot_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfdn_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfstr_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfbot_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfbot_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    if( (isel.eq.pdfdn_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfastr_)  .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfastr_)  .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfadn_) ) then
       if( zz_fusion ) then
          amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfdn_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfdn_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    if( (isel.eq.pdfup_ .and. jsel.eq.pdfchm_)  ) then
 
       !-- non-symmetric qq processes
       amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
       !--uc -> uc
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfup_,pdfchm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !--us -> us/cd
    if( (isel.eq.pdfup_ .and. jsel.eq.pdfbot_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfbot_) ) then
 
       amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg * tag1
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfup_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfchm_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    if( (isel.eq.pdfup_ .and. jsel.eq.pdfstr_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfdn_) ) then
 
       if( zz_fusion ) then
          amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(3,j2,4,j1,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfup_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfdn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !--ds -> ds
    if( (isel.eq.pdfdn_ .and. jsel.eq.pdfstr_) .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfbot_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfbot_) ) then
 
       amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfdn_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfdn_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfstr_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    !-- qbqb processes
    if( (isel.eq.pdfaup_ .and. jsel.eq.pdfachm_) ) then
       amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
       !--ubcb -> ubcb
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
 
       pdfindex = flip(iflip,pdfaup_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
       return
    endif
 
 
    !--ubsb -> ubsb//cbdb
    if( (isel.eq.pdfaup_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfabot_) ) then
       amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg * tag1
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfaup_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfachm_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    if( (isel.eq.pdfaup_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfadn_) ) then
 
       if( zz_fusion ) then
          amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(j1,4,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfaup_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfachm_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
 
       return
    endif
 
 
    !--dbsb -> dbsb
    if( (isel.eq.pdfadn_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfadn_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfastr_ .and. jsel.eq.pdfabot_) ) then
 
       amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
 
       pdfindex = flip(iflip,pdfadn_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfadn_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfastr_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
    !-------------
 
    j1 = 2
    j2 = 1
    iflip = 2
    !-- qqb processes
 
    !--uub -> uub // ddb
    if( (jsel.eq.pdfup_ .and. isel.eq.pdfaup_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfachm_) .or. (jsel.eq.pdfup_ .and. isel.eq.pdfachm_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfaup_) ) then
 
       if( zz_fusion ) then
          amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.true.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfup_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfup_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
 
       return
    endif
 
 
    !--udb -> udb
    if( (jsel.eq.pdfup_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfup_ .and. isel.eq.pdfastr_) .or. &
         (jsel.eq.pdfchm_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfup_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfabot_) ) then
       amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
 
       pdfindex = flip(iflip,pdfup_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfup_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfup_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfchm_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    !--dub -> dub
    if( (jsel.eq.pdfdn_ .and. isel.eq.pdfaup_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfachm_)  .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfachm_)  .or. &
         (jsel.eq.pdfstr_ .and. isel.eq.pdfaup_)  .or. (jsel.eq.pdfbot_ .and. isel.eq.pdfaup_)  .or. (jsel.eq.pdfbot_ .and. isel.eq.pdfachm_) ) then
       amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
 
       pdfindex = flip(iflip,pdfdn_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfdn_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfbot_,pdfaup_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfbot_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    !--ddb -> uub/ddb
    if( (jsel.eq.pdfbot_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfabot_) .or.&
         (jsel.eq.pdfbot_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfbot_ .and. isel.eq.pdfastr_) ) then
       amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
 
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg * tag1
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
 
       pdfindex = flip(iflip,pdfbot_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfdn_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfstr_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfbot_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfbot_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
 
 
    if( (jsel.eq.pdfdn_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfastr_)  .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfastr_)  .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfadn_) ) then
       if( zz_fusion ) then
          amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfdn_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfdn_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfstr_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    if( (jsel.eq.pdfup_ .and. isel.eq.pdfchm_) ) then
 
       !-- non-symmetric qq processes
       amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
       !--uc -> uc
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfup_,pdfchm_)
       res(pdfindex(1),pdfindex(2)) = restmp
       return
    endif
 
 
    !--us -> us/cd
    if( (jsel.eq.pdfup_ .and. isel.eq.pdfbot_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfbot_) ) then
 
       amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg * tag1
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfup_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfchm_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    if( (jsel.eq.pdfup_ .and. isel.eq.pdfstr_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfdn_) ) then
 
       if( zz_fusion ) then
          amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(3,j2,4,j1,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfup_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfchm_,pdfdn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !--ds -> ds
    if( (jsel.eq.pdfdn_ .and. isel.eq.pdfstr_) .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfbot_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfbot_) ) then
 
       amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfdn_,pdfstr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfdn_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfstr_,pdfbot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    !-- qbqb processes
    if( (jsel.eq.pdfaup_ .and. isel.eq.pdfachm_) ) then
       amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
       !--ubcb -> ubcb
       restmp = ((abs(amp_z(-1,-1))**2)+ &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
 
       pdfindex = flip(iflip,pdfaup_,pdfachm_)
       res(pdfindex(1),pdfindex(2)) = restmp
       return
    endif
 
 
    !--ubsb -> ubsb//cbdb
    if( (jsel.eq.pdfaup_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfachm_ .and. isel.eq.pdfabot_) ) then
       amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg * tag1
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
 
       pdfindex = flip(iflip,pdfaup_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfachm_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
 
    if( (jsel.eq.pdfaup_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfachm_ .and. isel.eq.pdfadn_) ) then
       if( zz_fusion ) then
          amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
          restmp = ((abs(amp_z(-1,-1))**2) + &
               (abs(amp_z(-1,+1))**2) + &
               (abs(amp_z(+1,-1))**2) + &
               (abs(amp_z(+1,+1))**2)) * spinavg * tag1
       else
          amp_w = a0_ww_4f(j1,4,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
          restmp= ( &
                  abs(amp_w(-1,-1))**2 + &
                  abs(amp_w(-1,+1))**2 + &
                  abs(amp_w(+1,-1))**2 + &
                  abs(amp_w(+1,+1))**2   &
                  ) * spinavg * tag2
       endif
 
       pdfindex = flip(iflip,pdfaup_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfachm_,pdfadn_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       return
    endif
 
 
    !--dbsb -> dbsb
    if( (jsel.eq.pdfadn_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfadn_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfastr_ .and. isel.eq.pdfabot_) ) then
 
       amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
       restmp = ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * spinavg
 
       if( .not. zz_fusion ) restmp = 0.0d0
 
       pdfindex = flip(iflip,pdfadn_,pdfastr_)
       res(pdfindex(1),pdfindex(2)) = restmp
 
       pdfindex = flip(iflip,pdfadn_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
 
       pdfindex = flip(iflip,pdfastr_,pdfabot_)
       res(pdfindex(1),pdfindex(2)) = restmp * tagbot
       return
    endif
 
    RETURN

evalamp_wbfh_unsymm_sa_select_exact()

subroutine, public modhiggsjj::evalamp_wbfh_unsymm_sa_select_exact	(	real(dp), dimension(4,5), intent(in)	p,
		integer, intent(in)	iSel,
		integer, intent(in)	jSel,
		integer, intent(in)	rSel,
		integer, intent(in)	sSel,
		real(dp), intent(out)	res
	)

Definition at line 1872 of file mod_HiggsJJ.F90.

    implicit none
    real(dp), intent(in) :: p(4,5)
    real(dp), intent(out) :: res
    integer, intent(in) :: iSel,jSel,rSel,sSel
    complex(dp) :: amp_z(-1:1,-1:1), amp_z_b(-1:1,-1:1)
    complex(dp) :: amp_w(-1:1,-1:1)
    real(dp) :: sprod(4,4)
    complex(dp) :: za(4,4), zb(4,4)
    real(dp) :: restmp
    real(dp) :: ckm_wfactor
    real(dp) :: kfactor_z,kfactor_w
    integer :: i, j, jz1, jz2, jw1, jw2, pdfindex(2)
    integer :: line_i,line_j
    integer :: kz1, kz2, kw1, kw2
    logical :: ZZ_fusion,WW_fusion
 
    !print *, "Begin EvalAmp_WBFH_UnSymm_SA_Select_exact"
    !print *, "iSel: ",iSel,", jSel: ",jSel," rSel: ",rSel," sSel: ",sSel
 
    res=0.0_dp
 
    zz_fusion=.false.
    ww_fusion=.false.
 
    ckm_wfactor=1.0_dp
    kfactor_z=1.0_dp
    kfactor_w=1.0_dp
    restmp=0.0_dp
 
    line_i=0
    line_j=0
    jz1 = 1
    jz2 = 2
    jw1 = jz1
    jw2 = jz2
 
    if( (abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_) .or. (abs(rsel).eq.pdftop_ .or. abs(ssel).eq.pdftop_) ) return
    if( (abs(isel).eq.pdfglu_ .or. abs(jsel).eq.pdfglu_) .or. (abs(rsel).eq.pdfglu_ .or. abs(ssel).eq.pdfglu_) ) return
 
    if( &
         (isel.eq.rsel .and. jsel.eq.ssel) &
         .or. &
         (isel.eq.ssel .and. jsel.eq.rsel) &
         ) then
       zz_fusion=.true.
       if( isel.eq.ssel .and. jsel.eq.rsel ) then
          kz1 = 4
          kz2 = 3
          kfactor_z = 1.0_dp ! dsqrt(ScaleFactor(iSel,sSel)*ScaleFactor(jSel,rSel))
          !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isZZ and is-jr"
       else
          kz1 = 3
          kz2 = 4
          kfactor_z = 1.0_dp ! dsqrt(ScaleFactor(iSel,rSel)*ScaleFactor(jSel,sSel))
          !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isZZ and ir-js"
       endif
    endif
 
    if( &
         ( (sign(isel,rsel).eq.isel .and. sign(jsel,ssel).eq.jsel) .and. (abs(isel-rsel).eq.1 .or. abs(isel-rsel).eq.3 .or. abs(isel-rsel).eq.5) .and. (abs(jsel-ssel).eq.1 .or. abs(jsel-ssel).eq.3 .or. abs(jsel-ssel).eq.5) ) &
         .or. &
         ( (sign(isel,ssel).eq.isel .and. sign(jsel,rsel).eq.jsel) .and. (abs(isel-ssel).eq.1 .or. abs(isel-ssel).eq.3 .or. abs(isel-ssel).eq.5) .and. (abs(jsel-rsel).eq.1 .or. abs(jsel-rsel).eq.3 .or. abs(jsel-rsel).eq.5) ) &
         ) then
       ww_fusion=.true.
       ! W_is W_jr fusion
       if( (sign(isel,ssel).eq.isel .and. sign(jsel,rsel).eq.jsel) .and. (abs(isel-ssel).eq.1 .or. abs(isel-ssel).eq.3 .or. abs(isel-ssel).eq.5) .and. (abs(jsel-rsel).eq.1 .or. abs(jsel-rsel).eq.3 .or. abs(jsel-rsel).eq.5) ) then
          kw1 = 4
          kw2 = 3
          kfactor_w = 1.0_dp ! dsqrt(ScaleFactor(iSel,sSel)*ScaleFactor(jSel,rSel))
 
          !if(ZZ_fusion) then ! Special treatment for WW+ZZ interference, not included through phasespace
          ckm_wfactor = ckmbare(isel,ssel)*ckmbare(jsel,rsel)
          !print *, "iSel, sSel: ",iSel," ",sSel, "; jSel, rSel: ",jSel," ",rSel, ", ckm: ",ckm_wfactor
          !endif
          !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isWW and is-jr"
       else ! W_ir W_js fusion
          kw1 = 3
          kw2 = 4
          kfactor_w = 1.0_dp ! dsqrt(ScaleFactor(iSel,rSel)*ScaleFactor(jSel,sSel))
 
          ckm_wfactor = ckmbare(isel,rsel)*ckmbare(jsel,ssel)
          !print *, "iSel, rSel: ",iSel," ",rSel, "; jSel, sSel: ",jSel," ",sSel, ", ckm: ",ckm_wfactor
          !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isWW and ir-js"
       endif
    endif
 
 
    if( .not.(zz_fusion .or. ww_fusion) ) return
    if(isel.lt.0) then
       if(zz_fusion) call swapi(jz1,kz1)
       if(ww_fusion) call swapi(jw1,kw1)
    endif
    if(jsel.lt.0) then
       if(zz_fusion) call swapi(jz2,kz2)
       if(ww_fusion) call swapi(jw2,kw2)
    endif
 
    line_i = isel
    line_j = jsel
 
    if(ww_fusion) then
       if (isel.eq.pdfup_ .or. isel.eq.pdfchm_ .or. isel.eq.pdfadn_ .or. isel.eq.pdfastr_ .or. isel.eq.pdfabot_) then ! W+ should be passed as the second set of partons
          call swapi(jw1,jw2)
          call swapi(kw1,kw2)
          call swapi(line_i,line_j)
          if(zz_fusion) then ! If also ZZ fusion, swap it as well
             call swapi(jz1,jz2)
             call swapi(kz1,kz2)
          endif
       endif
    endif
 
    call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
 
    !if(ZZ_fusion) print *, "jz1: ",jz1,", jz2: ",jz2," kz1: ",kz1," kz2: ",kz2
    !if(WW_fusion) print *, "jw1: ",jw1,", jw2: ",jw2," kw1: ",kw1," kw2: ",kw2
 
    !ckm_wfactor = 1.0_dp ! TEST!!!
    if(zz_fusion) then
 
       amp_z = a0_zz_4f(kz1,jz1,kz2,jz2,za,zb,sprod,line_i,line_j)
 
       restmp = restmp + &
            ((abs(amp_z(-1,-1))**2) + &
            (abs(amp_z(-1,+1))**2) + &
            (abs(amp_z(+1,-1))**2) + &
            (abs(amp_z(+1,+1))**2)) * xn**2
 
       if(isel.eq.jsel) then
          amp_z_b = -a0_zz_4f(kz2,jz1,kz1,jz2,za,zb,sprod,line_i,line_j)
          restmp = restmp + &
               ((abs(amp_z_b(-1,-1))**2) + &
               (abs(amp_z_b(-1,+1))**2) + &
               (abs(amp_z_b(+1,-1))**2) + &
               (abs(amp_z_b(+1,+1))**2)) * xn**2
          restmp = restmp + &
               (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
               two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
 
          restmp = restmp * symmfac
       endif
       restmp = restmp * kfactor_z**2
    endif
 
    if(ww_fusion) then
       amp_w = -a0_ww_4f(kw1,jw1,kw2,jw2,za,zb,sprod,line_i,line_j,usewwcoupl=.true.)
       restmp = restmp + &
                (        &
                abs(amp_w(-1,-1))**2 + &
                abs(amp_w(-1,+1))**2 + &
                abs(amp_w(+1,-1))**2 + &
                abs(amp_w(+1,+1))**2   &
                ) * xn**2 * kfactor_w**2 * abs(ckm_wfactor)**2
       if(zz_fusion) then
          restmp = restmp + &
                   two * (  &
                   real(amp_w(-1,-1)*ckm_wfactor*conjg(amp_z(-1,-1)),kind=dp) + &
                   real(amp_w(-1,+1)*ckm_wfactor*conjg(amp_z(-1,+1)),kind=dp) + &
                   real(amp_w(+1,-1)*ckm_wfactor*conjg(amp_z(+1,-1)),kind=dp) + &
                   real(amp_w(+1,+1)*ckm_wfactor*conjg(amp_z(+1,+1)),kind=dp)   &
                   ) * xn * kfactor_z * kfactor_w
       endif
    endif
 
    restmp = restmp * aveqq
    if(abs(isel).eq.pdfbot_ .or. abs(jsel).eq.pdfbot_) restmp = restmp * tagbot
    res = restmp
 
    RETURN

flip()

integer function, dimension(2) modhiggsjj::flip ( integer  i, integer  a1, integer  a2  )

private

Definition at line 2252 of file mod_HiggsJJ.F90.

    integer :: flip(2)
    integer :: i, a1, a2
 
    if (i .eq. 1) flip = (/a1,a2/)
    if (i .eq. 2) flip = (/a2,a1/)
 
    return
 

me2_ggggh()

subroutine modhiggsjj::me2_ggggh	(	integer, intent(in)	j1,
		integer, intent(in)	j2,
		integer, intent(in)	j3,
		integer, intent(in)	j4,
		complex(dp), dimension(4,4), intent(in)	za,
		complex(dp), dimension(4,4), intent(in)	zb,
		real(dp), dimension(4,4), intent(in)	sprod,
		real(dp), intent(out)	res
	)

Definition at line 2814 of file mod_HiggsJJ.F90.

    integer, intent(in) :: j1,j2,j3,j4
    complex(dp), intent(in) :: za(4,4), zb(4,4)
    real(dp), intent(in) :: sprod(4,4)
    real(dp), intent(out) :: res
    complex(dp) :: a1234(-1:1,-1:1,-1:1,-1:1), a1324(-1:1,-1:1,-1:1,-1:1)
    complex(dp) :: aphi1234(1:2,-1:1,-1:1,-1:1,-1:1), aphi1324(1:2,-1:1,-1:1,-1:1,-1:1)
    real(dp), parameter :: col = 8.0_dp * ca**3 * cf
    integer :: i1,i2,i3,i4
    complex(dp) :: scalar, pseudo
 
    res = zero
 
    scalar = ghg2
    pseudo = -ghg4
 
    aphi1234 = a0phigggg_xxxx(j1,j2,j3,j4,za,zb,sprod)
    aphi1324 = a0phigggg_xxxx(j1,j3,j2,j4,za,zb,sprod)
 
    a1234(:,:,:,:) = scalar * (aphi1234(1,:,:,:,:) + aphi1234(2,:,:,:,:)) + &
         pseudo * (-ci) * (aphi1234(1,:,:,:,:) - aphi1234(2,:,:,:,:))
 
    a1324(:,:,:,:) = scalar * (aphi1324(1,:,:,:,:) + aphi1324(2,:,:,:,:)) + &
         pseudo * (-ci) * (aphi1324(1,:,:,:,:) - aphi1324(2,:,:,:,:))
 
    do i1 = -1,1,2
    do i2 = -1,1,2
    do i3 = -1,1,2
    do i4 = -1,1,2
 
       res = res + real(a1234(i1,i2,i3,i4)*conjg(a1234(i1,i2,i3,i4)),kind=dp)
       res = res + real(a1324(i1,i2,i3,i4)*conjg(a1324(i1,i2,i3,i4)),kind=dp)
       res = res + real(a1234(i1,i2,i3,i4)*conjg(a1324(i1,i3,i2,i4)),kind=dp)
 
    enddo
    enddo
    enddo
    enddo
 
    res = res * col / vev**2
 
    return
 

me2_qbqggh()

subroutine modhiggsjj::me2_qbqggh ( integer, intent(in)  j1, integer, intent(in)  j2, integer, intent(in)  j3, integer, intent(in)  j4, complex(dp), dimension(4,4), intent(in)  za, complex(dp), dimension(4,4), intent(in)  zb, real(dp), dimension(4,4), intent(in)  sprod, real(dp), intent(out)  res  )

private

Definition at line 2859 of file mod_HiggsJJ.F90.

    integer, intent(in) :: j1,j2,j3,j4
    complex(dp), intent(in) :: za(4,4), zb(4,4)
    real(dp), intent(in) :: sprod(4,4)
    real(dp), intent(out) :: res
    complex(dp) :: a1234(-1:1,-1:1,-1:1), a1243(-1:1,-1:1,-1:1)
    complex(dp) :: aphi1234(1:2,-1:1,-1:1,-1:1), aphi1243(1:2,-1:1,-1:1,-1:1)
    real(dp), parameter :: col1 = 4.0_dp * xn * cf**2
    real(dp), parameter :: col2 = 2.0_dp * xn * cf * (2.0_dp * cf - ca)
    integer :: i12,i3,i4
    complex(dp) :: scalar, pseudo
 
    res = zero
 
    scalar = ghg2
    pseudo = -ghg4
 
    aphi1234 = a0phiqbqgg_xxx(j1,j2,j3,j4,za,zb,sprod)
    aphi1243 = a0phiqbqgg_xxx(j1,j2,j4,j3,za,zb,sprod)
 
    a1234(:,:,:) = scalar * (aphi1234(1,:,:,:) + aphi1234(2,:,:,:)) + &
         pseudo * (-ci) * (aphi1234(1,:,:,:) - aphi1234(2,:,:,:))
 
    a1243(:,:,:) = scalar * (aphi1243(1,:,:,:) + aphi1243(2,:,:,:)) + &
         pseudo * (-ci) * (aphi1243(1,:,:,:) - aphi1243(2,:,:,:))
 
    do i12 = -1,1,2
    do i3 = -1,1,2
    do i4 = -1,1,2
 
       res = res + real(a1234(i12,i3,i4)*conjg(a1234(i12,i3,i4)),kind=dp) * col1
       res = res + real(a1243(i12,i4,i3)*conjg(a1243(i12,i4,i3)),kind=dp) * col1
       res = res + two * real(a1234(i12,i3,i4)*conjg(a1243(i12,i4,i3)),kind=dp) * col2
 
    enddo
    enddo
    enddo
 
    res = res / vev**2
 
    return
 

me2_qbqqbq()

subroutine modhiggsjj::me2_qbqqbq ( integer, intent(in)  j1, integer, intent(in)  j2, integer, intent(in)  j3, integer, intent(in)  j4, complex(dp), dimension(4,4), intent(in)  za, complex(dp), dimension(4,4), intent(in)  zb, real(dp), dimension(4,4), intent(in)  sprod, real(dp), intent(out)  res_diff, real(dp), intent(out)  res_id  )

private

Definition at line 2903 of file mod_HiggsJJ.F90.

    integer, intent(in) :: j1,j2,j3,j4
    complex(dp), intent(in) :: za(4,4), zb(4,4)
    real(dp), intent(in) :: sprod(4,4)
    real(dp), intent(out) :: res_diff, res_id
    complex(dp) :: amp_a(-1:1,-1:1), amp_b(-1:1,-1:1)
    complex(dp) :: aphi_a(1:2,-1:1,-1:1), aphi_b(1:2,-1:1,-1:1)
    real(dp), parameter :: col = xn**2 - one
    integer :: h12,h34
    complex(dp) :: scalar, pseudo
 
    res_diff = zero
    res_id = zero
 
    scalar = ghg2
    pseudo = -ghg4
 
    aphi_a = a0phiqbqqbq_xx(j1,j2,j3,j4,za,zb,sprod)
    aphi_b = -a0phiqbqqbq_xx(j1,j4,j3,j2,za,zb,sprod)
 
    amp_a(:,:) = scalar * (aphi_a(1,:,:) + aphi_a(2,:,:)) + &
         pseudo * (-ci) * (aphi_a(1,:,:) - aphi_a(2,:,:))
 
    amp_b(:,:) = scalar * (aphi_b(1,:,:) + aphi_b(2,:,:)) + &
         pseudo * (-ci) * (aphi_b(1,:,:) - aphi_b(2,:,:))
 
    do h12 = -1,1,2
    do h34 = -1,1,2
       res_diff = res_diff + real(amp_a(h12,h34)*conjg(amp_a(h12,h34)),kind=dp)
       res_id = res_id + real(amp_a(h12,h34)*conjg(amp_a(h12,h34)),kind=dp)
       res_id = res_id + real(amp_b(h12,h34)*conjg(amp_b(h12,h34)),kind=dp)
    enddo
    res_id = res_id - two/xn * real(amp_a(h12,h12)*conjg(amp_b(h12,h12)),kind=dp)
    enddo
 
    res_id = res_id * col / vev**2
    res_diff = res_diff * col / vev**2
 
    return
 

spinoru2()

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

private

Definition at line 3118 of file mod_HiggsJJ.F90.

    implicit none
    integer, intent(in) :: n
    real(dp), intent(in) :: p(4,n)
    complex(dp), intent(out) :: za(n,n), zb(n,n)
    real(dp), intent(out) :: s(n,n)
    integer :: i,j
    complex(dp) :: c23(n), f(n)
    real(dp) :: 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)=cmplx(p(4,j),-p(3,j),kind=dp)
          f(j)=(one,zero)
       else
       !-----negative energy case
          rt(j)=sqrt(abs(-p(1,j)-p(2,j)))
          c23(j)=cmplx(-p(4,j),p(3,j),kind=dp)
          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)*cmplx(rt(j)/(rt(i)+1d-16),kind=dp)-c23(j)*cmplx(rt(i)/(rt(j)+1d-16),kind=dp) )
 
          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)=-cmplx(s(i,j),kind=dp)/(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
 

Variable Documentation

tag1

real(dp), parameter, public modhiggsjj::tag1 = 1.0_dp

Definition at line 12 of file mod_HiggsJJ.F90.

   real(dp), public, parameter :: tag1 = 1.0_dp

tag2

real(dp), parameter, public modhiggsjj::tag2 = 1.0_dp

Definition at line 13 of file mod_HiggsJJ.F90.

   real(dp), public, parameter :: tag2 = 1.0_dp

tagbot

real(dp), parameter, public modhiggsjj::tagbot = 1.0_dp

Definition at line 14 of file mod_HiggsJJ.F90.

   real(dp), public, parameter :: tagbot = 1.0_dp