mod_TTBHiggs.F90

Go to the documentation of this file.

MODULE modttbhiggs
use modparameters
implicit none
 
 
public :: evalxsec_pp_ttbh, evalxsec_pp_bbbh, evalamp_gg_ttbh, evalamp_qqb_ttbh
public :: initprocess_ttbh
public :: exitprocess_ttbh
private
 
integer,parameter :: colorlesstag = 1
 
 
type :: particle
   integer  :: parttype
   integer  :: extref
   integer  :: helicity
   real(8)  :: mass
   real(8)  :: mass2
   complex(8) :: mom(1:4)
   complex(8) :: pol(1:4)
end type
 
type :: ptrtoparticle
   integer,pointer  :: parttype
   integer,pointer  :: extref
   real(8),pointer  :: mass
   real(8),pointer  :: mass2
   integer ,pointer :: helicity
   complex(8),pointer :: mom(:)
   complex(8),pointer :: pol(:)
end type
 
 
 
type :: treeprocess
   integer :: numpart
   integer :: numqua
   integer :: numsca
   integer :: numw
   integer :: numv
   integer :: bosonvertex
   integer,allocatable :: numglu(:)
   integer,allocatable :: partref(:)
   integer,allocatable :: parttype(:)
   type(ptrtoparticle),allocatable :: quarks(:)
   type(ptrtoparticle),allocatable :: gluons(:)
   type(ptrtoparticle) :: boson
   type(ptrtoparticle),allocatable :: scalars(:)
end type
 
 
INTERFACE OPERATOR (.ndot.)
   module procedure fourvecdot
END INTERFACE OPERATOR (.Ndot.)
 
 
real(8), parameter :: propcut = 1.0d-8
integer, parameter :: dv=4,ds=4
 
type(particle),save    :: extparticles(1:7)
type(treeprocess),save :: thetreeamps_gg_ttbh(1:2)
type(treeprocess),save :: thetreeamps_qqb_ttbh(1:1)
 
 
complex(8) :: couplhtt_right_dyn,couplhtt_left_dyn
 
 
 
 
   CONTAINS
 
 
SUBROUTINE evalxsec_pp_ttbh(Mom,SelectProcess,Res)
implicit none
real(8), intent(in) :: mom(1:4,1:13)
integer,intent(in) :: selectprocess! 0=gg, 1=qqb, 2=all
real(8), intent(out) :: res(-5:5,-5:5)
real(8) :: matelsq_gg,matelsq_qqb,matelsq_qbq
integer :: iq
 
   call initprocess_ttbh()
 
   matelsq_qqb = 0d0
   matelsq_qbq = 0d0
   matelsq_gg  = 0d0
   if( selectprocess.eq.0 ) then
      call evalamp_gg_ttbh(mom(1:4,1:13),matelsq_gg)
   elseif( selectprocess.eq.1 ) then
      call evalamp_qqb_ttbh(mom(1:4,1:13),matelsq_qqb)
      matelsq_qbq = matelsq_qqb
   else
      call evalamp_gg_ttbh(mom(1:4,1:13),matelsq_gg)
      call evalamp_qqb_ttbh(mom(1:4,1:13),matelsq_qqb)
      matelsq_qbq = matelsq_qqb
   endif
   do iq=0,5
      if(iq.eq.pdfglu_) then
         res(iq,iq) = matelsq_gg
      else
         res(iq,-iq) = matelsq_qqb
         res(-iq,iq) = matelsq_qbq
      endif
   enddo
 
   call exitprocess_ttbh()
   RETURN
END SUBROUTINE
 
SUBROUTINE evalxsec_pp_bbbh(Mom,SelectProcess,Res)
implicit none
real(8), intent(in) :: mom(1:4,1:13)
integer,intent(in) :: selectprocess! 0=gg, 1=qqb, 2=all
real(8), intent(out) :: res(-5:5,-5:5)
real(8) :: tmptopmass
integer :: tmptopdec
   tmptopmass = m_top
   tmptopdec = topdecays
 
   m_top = m_bot
   topdecays=0
   call evalxsec_pp_ttbh(mom,selectprocess,res)
 
   m_top = tmptopmass
   topdecays = tmptopdec
   RETURN
END SUBROUTINE
 
 
 
SUBROUTINE initprocess_ttbh()
implicit none
integer :: numquarks,numgluons,numboson
integer :: numtrees
integer :: itree,numparticles
 
 
! gg->ttbar+H
  numquarks=2; numgluons=2; numboson=1;
  numparticles=numquarks+numgluons+numboson
  numtrees=2
  do itree=1,numtrees
      thetreeamps_gg_ttbh(itree)%NumPart=numparticles
      thetreeamps_gg_ttbh(itree)%NumQua=numquarks
      allocate( thetreeamps_gg_ttbh(itree)%NumGlu(0:numquarks+1) )
      allocate( thetreeamps_gg_ttbh(itree)%PartRef(1:numparticles) )
      allocate( thetreeamps_gg_ttbh(itree)%PartType(1:numparticles) )
      allocate( thetreeamps_gg_ttbh(itree)%Quarks(1:numquarks) )
      allocate( thetreeamps_gg_ttbh(itree)%Gluons(1:numgluons) )
      thetreeamps_gg_ttbh(itree)%NumGlu(0) = numgluons
      thetreeamps_gg_ttbh(itree)%NumSca = 0
  enddo
 
! qqbar->ttbar+H
  numquarks=4; numgluons=0; numboson=1;
  numparticles=numquarks+numgluons+numboson
  numtrees=1
  do itree=1,numtrees
      thetreeamps_qqb_ttbh(itree)%NumPart=numparticles
      thetreeamps_qqb_ttbh(itree)%NumQua=numquarks
      allocate( thetreeamps_qqb_ttbh(itree)%NumGlu(0:numquarks+1) )
      allocate( thetreeamps_qqb_ttbh(itree)%PartRef(1:numparticles) )
      allocate( thetreeamps_qqb_ttbh(itree)%PartType(1:numparticles) )
      allocate( thetreeamps_qqb_ttbh(itree)%Quarks(1:numquarks) )
      allocate( thetreeamps_qqb_ttbh(itree)%Gluons(1:numgluons) )
      thetreeamps_qqb_ttbh(itree)%NumGlu(0) = numgluons
      thetreeamps_qqb_ttbh(itree)%NumSca = 0
  enddo
 
  extparticles(1)%PartType = atop_
  extparticles(1)%ExtRef   = 1
  extparticles(1)%Mass = m_top
  extparticles(1)%Mass2= extparticles(1)%Mass**2
  extparticles(1)%Helicity = 0
 
  extparticles(2)%PartType = top_
  extparticles(2)%ExtRef   = 2
  extparticles(2)%Mass = m_top
  extparticles(2)%Mass2= extparticles(2)%Mass**2
  extparticles(2)%Helicity = 0
 
  extparticles(3)%PartType = glu_
  extparticles(3)%ExtRef   = 3
  extparticles(3)%Mass = 0d0
  extparticles(3)%Mass2= 0d0
  extparticles(3)%Helicity = 0
 
  extparticles(4)%PartType = glu_
  extparticles(4)%ExtRef   = 4
  extparticles(4)%Mass = 0d0
  extparticles(4)%Mass2= 0d0
  extparticles(4)%Helicity = 0
 
  extparticles(5)%PartType = astr_
  extparticles(5)%ExtRef   = 5
  extparticles(5)%Mass = 0d0
  extparticles(5)%Mass2= 0d0
  extparticles(5)%Helicity = 0
 
  extparticles(6)%PartType = str_
  extparticles(6)%ExtRef   = 6
  extparticles(6)%Mass = 0d0
  extparticles(6)%Mass2= 0d0
  extparticles(6)%Helicity = 0
 
  extparticles(7)%PartType = hig_
  extparticles(7)%ExtRef   = 7
  extparticles(7)%Mass = m_reso
  extparticles(7)%Mass2= extparticles(5)%Mass**2
  extparticles(7)%Helicity = 0
 
 
  thetreeamps_gg_ttbh(1)%PartRef(1:5) = (/3,4,1,7,2/)! (/1,7,2,3,4/)
  thetreeamps_gg_ttbh(2)%PartRef(1:5) = (/3,1,7,2,4/)! (/1,7,2,4,3/)
  do itree=1,2
      call linktreeparticles(thetreeamps_gg_ttbh(itree),extparticles(1:7))
  enddo
 
  thetreeamps_qqb_ttbh(1)%PartRef(1:5) = (/5,6,1,7,2/)! (/1,7,2,5,6/)
  call linktreeparticles(thetreeamps_qqb_ttbh(1),extparticles(1:7))
 
 
  couplhtt_right_dyn = m_top/vev/2d0 * ( kappa + (0d0,1d0)*kappa_tilde )
  couplhtt_left_dyn  = m_top/vev/2d0 * ( kappa - (0d0,1d0)*kappa_tilde )
 
 
RETURN
END SUBROUTINE
 
 
SUBROUTINE exitprocess_ttbh()
implicit none
integer :: numtrees, itree
 
 
! gg->ttbar+H
  numtrees=2
  do itree=1,numtrees
      !call UnLinkTreeParticles(TheTreeAmps_GG_TTBH(iTree))
 
      if(allocated( thetreeamps_gg_ttbh(itree)%NumGlu )) then
         deallocate( thetreeamps_gg_ttbh(itree)%NumGlu )
      endif
      if(allocated( thetreeamps_gg_ttbh(itree)%PartRef )) then
         deallocate( thetreeamps_gg_ttbh(itree)%PartRef )
      endif
      if(allocated( thetreeamps_gg_ttbh(itree)%PartType )) then
         deallocate( thetreeamps_gg_ttbh(itree)%PartType )
      endif
      if(allocated( thetreeamps_gg_ttbh(itree)%Quarks )) then
         deallocate( thetreeamps_gg_ttbh(itree)%Quarks )
      endif
      if(allocated( thetreeamps_gg_ttbh(itree)%Gluons )) then
         deallocate( thetreeamps_gg_ttbh(itree)%Gluons )
      endif
      if(allocated( thetreeamps_gg_ttbh(itree)%Scalars )) then
         deallocate( thetreeamps_gg_ttbh(itree)%Scalars )
      endif
  enddo
 
  numtrees=1
  do itree=1,numtrees
      !call UnLinkTreeParticles(TheTreeAmps_QQB_TTBH(iTree))
 
      if(allocated( thetreeamps_qqb_ttbh(itree)%NumGlu )) then
         deallocate( thetreeamps_qqb_ttbh(itree)%NumGlu )
      endif
      if(allocated( thetreeamps_qqb_ttbh(itree)%PartRef )) then
         deallocate( thetreeamps_qqb_ttbh(itree)%PartRef )
      endif
      if(allocated( thetreeamps_qqb_ttbh(itree)%PartType )) then
         deallocate( thetreeamps_qqb_ttbh(itree)%PartType )
      endif
      if(allocated( thetreeamps_qqb_ttbh(itree)%Quarks )) then
         deallocate( thetreeamps_qqb_ttbh(itree)%Quarks )
      endif
      if(allocated( thetreeamps_qqb_ttbh(itree)%Gluons )) then
         deallocate( thetreeamps_qqb_ttbh(itree)%Gluons )
      endif
      if(allocated( thetreeamps_qqb_ttbh(itree)%Scalars )) then
         deallocate( thetreeamps_qqb_ttbh(itree)%Scalars )
      endif
  enddo
 
 
RETURN
END SUBROUTINE
 
 
 
SUBROUTINE evalamp_gg_ttbh(Mom,SqAmp)
use modtopdecay
implicit none
real(8) :: mom(1:4,1:13),sqamp
complex(8) :: resoffsh(1:4,1:2),res(1:2,1:2)
complex(8) :: glupol(1:4,1:2,1:2)
integer :: hel4,tophel1,tophel2,nhel
real(8),parameter :: c_aa=64.d0/3.d0, c_ab=-8.d0/3.d0
integer, parameter :: inleft=1,inright=2,hbos=3,tbar=4,t=5,  bbar=6,wm=7,lepm=8,nubar=9,  b=10,wp=11,lepp=12,nu=13
sqamp = 0d0
 
 
     extparticles(1)%Mom(1:4) = mom(1:4,tbar)
     extparticles(2)%Mom(1:4) = mom(1:4,t)
     extparticles(3)%Mom(1:4) =-mom(1:4,inleft)
     extparticles(4)%Mom(1:4) =-mom(1:4,inright)
     extparticles(7)%Mom(1:4) = mom(1:4,hbos)
 
     if( topdecays.ne.0 ) then
        call topdecay(atop_,(/mom(1:4,bbar),mom(1:4,lepm),mom(1:4,nubar)/),extparticles(1)%Pol(1:4))
        call topdecay(top_,(/mom(1:4,b),mom(1:4,lepp),mom(1:4,nu)/),extparticles(2)%Pol(1:4))
     endif
     extparticles(7)%Pol(1:4) = 1d0
!    call HDecay(ExtParticles(7),DK_LO,MomExt(1:4,12:13))
     glupol(1:4,1,1) = pol_mless(extparticles(3)%Mom(1:4),+1,outgoing=.true.)
     glupol(1:4,1,2) = pol_mless(extparticles(3)%Mom(1:4),-1,outgoing=.true.)
     glupol(1:4,2,1) = pol_mless(extparticles(4)%Mom(1:4),+1,outgoing=.true.)
     glupol(1:4,2,2) = pol_mless(extparticles(4)%Mom(1:4),-1,outgoing=.true.)
!      GluPol(1:4,1,1) = ExtParticles(3)%Mom(1:4);  GluPol(1:4,1,2) = ExtParticles(3)%Mom(1:4); print *, "checking gauge invariance"
 
 
     nhel=-1
     if( topdecays.EQ.0 ) nhel=+1
     do tophel1=-1,nhel,2
     do tophel2=-1,nhel,2
     if( topdecays.eq.0 ) then
             call ubarspi_dirac(extparticles(2)%Mom(1:4),m_top,tophel1,extparticles(2)%Pol(1:4))
             call    vspi_dirac(extparticles(1)%Mom(1:4),m_top,tophel2,extparticles(1)%Pol(1:4))
     endif
     do hel4=1,2
 
        extparticles(4)%Pol(1:4) = glupol(1:4,2,hel4)
!         ExtParticles(4)%Pol(1:4) = ExtParticles(4)%Mom(1:4); print *, "checking gauge invariance"
        call new_calc_ampl(0,0,thetreeamps_gg_ttbh(1),resoffsh(1:4,1))
        call new_calc_ampl(0,0,thetreeamps_gg_ttbh(2),resoffsh(1:4,2))
 
        res(1,1) = (resoffsh(1:4,1).ndot.glupol(1:4,1,1))! col1 hel+
        res(2,1) = (resoffsh(1:4,2).ndot.glupol(1:4,1,1))! col2 hel+
        res(1,2) = (resoffsh(1:4,1).ndot.glupol(1:4,1,2))! col1 hel-
        res(2,2) = (resoffsh(1:4,2).ndot.glupol(1:4,1,2))! col2 hel-
 
        sqamp = sqamp   &
              + c_aa * dreal( res(1,1)*dconjg(res(1,1)) + res(1,2)*dconjg(res(1,2)) )  &
              + c_ab * dreal( res(1,1)*dconjg(res(2,1)) + res(1,2)*dconjg(res(2,2)) )  &
              + c_ab * dreal( res(2,1)*dconjg(res(1,1)) + res(2,2)*dconjg(res(1,2)) )  &
              + c_aa * dreal( res(2,1)*dconjg(res(2,1)) + res(2,2)*dconjg(res(2,2)) )
    enddo
    enddo
    enddo
 
    sqamp = sqamp * spinavg * gluoncolavg**2 * (4d0*pi*alphas)**2  !* (4d0*pi*alpha_QED) * (m_top/(2d0*sitW*M_W))**2
 
 
RETURN
END SUBROUTINE
 
 
 
 
 
 
 
SUBROUTINE evalamp_qqb_ttbh(Mom,SqAmp)
use modtopdecay
implicit none
real(8) :: mom(1:4,1:13),sqamp
complex(8) :: resoffsh(1:4),res(1:2)
complex(8) :: quapol(1:4,1:2,1:2)
integer :: hel4,tophel1,tophel2,nhel
real(8),parameter :: c_aa=8.0d0
integer, parameter :: inleft=1,inright=2,hbos=3,tbar=4,t=5,  bbar=6,wm=7,lepm=8,nubar=9,  b=10,wp=11,lepp=12,nu=13
sqamp = 0d0
 
     extparticles(1)%Mom(1:4) = mom(1:4,tbar)
     extparticles(2)%Mom(1:4) = mom(1:4,t)
     extparticles(5)%Mom(1:4) =-mom(1:4,inleft)
     extparticles(6)%Mom(1:4) =-mom(1:4,inright)
     extparticles(7)%Mom(1:4) = mom(1:4,hbos)
 
     if( topdecays.ne.0 ) then
        call topdecay(atop_,(/mom(1:4,bbar),mom(1:4,lepm),mom(1:4,nubar)/),extparticles(1)%Pol(1:4))
        call topdecay(top_,(/mom(1:4,b),mom(1:4,lepp),mom(1:4,nu)/),extparticles(2)%Pol(1:4))
     endif
     extparticles(7)%Pol(1:4) = 1d0
!    call HDecay(ExtParticles(7),DK_LO,MomExt(1:4,12:13))
     call ubarspi_dirac(extparticles(6)%Mom(1:4),0d0,-1,quapol(1:4,1,1))
     call ubarspi_dirac(extparticles(6)%Mom(1:4),0d0,+1,quapol(1:4,1,2))
     call    vspi_dirac(extparticles(5)%Mom(1:4),0d0,-1,quapol(1:4,2,1))
     call    vspi_dirac(extparticles(5)%Mom(1:4),0d0,+1,quapol(1:4,2,2))
 
 
 
     nhel=-1
     if( topdecays.EQ.0 ) nhel=+1
     do tophel1=-1,nhel,2
     do tophel2=-1,nhel,2
     if( topdecays.eq.0 ) then
             call ubarspi_dirac(extparticles(2)%Mom(1:4),m_top,tophel1,extparticles(2)%Pol(1:4))
             call    vspi_dirac(extparticles(1)%Mom(1:4),m_top,tophel2,extparticles(1)%Pol(1:4))
     endif
     do hel4=1,2
 
        extparticles(6)%Pol(1:4) = quapol(1:4,2,hel4)
!         ExtParticles(4)%Pol(1:4) = ExtParticles(4)%Mom(1:4); print *, "checking gauge invariance"
        call new_calc_ampl(0,0,thetreeamps_qqb_ttbh(1),resoffsh(1:4))
 
        res(1) = psp1_(resoffsh(1:4),quapol(1:4,1,1))! hel+
        res(2) = psp1_(resoffsh(1:4),quapol(1:4,1,2))! hel-
 
        sqamp = sqamp   &
              + c_aa * dreal( res(1)*dconjg(res(1)) + res(2)*dconjg(res(2)) )
    enddo
    enddo
    enddo
    sqamp = sqamp * spinavg * quarkcolavg**2 * (4d0*pi*alphas)**2  !* (4d0*pi*alpha_QED) * (m_top/(2d0*sitW*M_W))**2
 
RETURN
END SUBROUTINE
 
 
 
 
 
 
 
 
SUBROUTINE new_calc_ampl(tag_f,tag_Z,TreeProc,Res)
implicit none
integer :: tag_f,tag_Z,n
complex(8) :: Res(1:Ds)
type(treeprocess) :: TreeProc
logical,parameter :: Boson=.true.
integer :: i,j,Order(1:6)
 
 
      if( treeproc%NumQua.eq.2 .and. treeproc%NumSca.eq.0 ) then!  2 quarks and no scalars
          if ( treeproc%PartType(1).eq.glu_ .and. boson ) then
             res(1:dv) = cur_g_2fv( treeproc%Gluons(2:treeproc%NumGlu(0)),treeproc%Quarks(1:treeproc%NumQua),treeproc%Boson,treeproc%NumGlu(0:3) )
          elseif( isaquark(treeproc%PartType(1)) .and. boson ) then
             if( treeproc%NumV.eq.1 ) then
                res(1:ds) = cur_f_2fv(treeproc%Gluons(1:treeproc%NumGlu(0)),treeproc%Quarks(2:2),treeproc%Quarks(1)%PartType,treeproc%Boson,treeproc%NumGlu(0:2))
             else
                print *, "requested current with a boson is not available"
                stop
             endif
          else
             print *, "requested current is not available 2q"
             stop
          endif
 
      elseif( treeproc%NumQua.eq.4  .and. treeproc%NumSca.eq.0) then!  4 quarks, no scalars
          if( treeproc%PartType(1).eq.glu_ ) then
                 res(1:dv) = cur_g_4fv( treeproc%Gluons(2:treeproc%NumGlu(0)),treeproc%Quarks(1:treeproc%NumQua),treeproc%Boson,treeproc%BosonVertex,treeproc%NumGlu(0:5) )
          elseif( isaquark(treeproc%PartType(1)) ) then
                 res(1:ds) = cur_f_4fv( treeproc%Gluons(1:treeproc%NumGlu(0)),treeproc%Quarks(2:4),treeproc%Quarks(1)%PartType,treeproc%Boson,treeproc%BosonVertex,treeproc%NumGlu(0:4),tag_f,tag_z )
          else
             print *, "requested current is not available 4q"
             stop
          endif
      else
           print *, "requested current is not available xx"
           stop
      endif
 
return
END SUBROUTINE
 
 
 
 
SUBROUTINE linktreeparticles(TheTreeAmp,TheParticles)
implicit none
type(treeprocess) :: TheTreeAmp
type(particle),target :: TheParticles(:)
integer :: iPart,PartRef,PartType,ig,iq,ib,NPart,counterQ,counterG,LastQuark,QuarkPos(1:6)
 
 
            thetreeamp%NumW = 0
            thetreeamp%NumV = 0
            counterq = 0
            counterg = 0
 
            do npart=1,thetreeamp%NumPart
                  thetreeamp%PartType(npart) = theparticles( thetreeamp%PartRef(npart) )%PartType
                  if( isaquark(thetreeamp%PartType(npart)) ) then
                     !TheTreeAmp%NumQua = TheTreeAmp%NumQua + 1   ! this is suppoed to be done outside this subroutine
                     counterq = counterq + 1
                     quarkpos(counterq) = counterq + counterg
                     lastquark = counterq! only required for BosonVertex below
!                   elseif( IsAScalar(TheTreeAmp%PartType(NPart)) ) then
! !                      TheTreeAmp%NumSca = TheTreeAmp%NumSca + 1   ! this is suppoed to be done outside this subroutine
!                      counterQ = counterQ + 1!     treat the scalar like a quark here because this is only to determine NumGlu
!                      QuarkPos(counterQ) = counterQ + counterG
                  elseif( thetreeamp%PartType(npart).eq.glu_ ) then
                     counterg = counterg + 1
                  elseif( isaboson(thetreeamp%PartType(npart)) ) then! careful: bosons should only be places *between* same flavor quark lines
                     if( npart.eq.1 ) print *, "Vector boson should not be the first particle."
                     if( abs(thetreeamp%PartType(npart)).eq.abs(wp_) ) thetreeamp%NumW = thetreeamp%NumW + 1
                     if( abs(thetreeamp%PartType(npart)).eq.abs(z0_) ) thetreeamp%NumV = thetreeamp%NumV + 1
                     if( abs(thetreeamp%PartType(npart)).eq.abs(pho_)) thetreeamp%NumV = thetreeamp%NumV + 1
!                      TheTreeAmp%BosonVertex = TheTreeAmp%PartType(LastQuark)! this variable specifies the position of the boson wrt. to the quark lines
                     thetreeamp%BosonVertex = lastquark
                  endif
            enddo
 
 
!           set number of gluons between quark lines
!            if( IsAQuark( TheTreeAmp%PartType(1) ) .or. IsAScalar(TheTreeAmp%PartType(1)) ) then ! not a gluon and not a boson
           if( isaquark( thetreeamp%PartType(1) ) ) then ! not a gluon and not a boson
             if( thetreeamp%NumQua+thetreeamp%NumSca .eq. 2 ) then
                   thetreeamp%NumGlu(1) = quarkpos(2) - quarkpos(1) - 1
                   thetreeamp%NumGlu(2) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(2)
             endif
             if( thetreeamp%NumQua .eq. 4 ) then
                   thetreeamp%NumGlu(1) = quarkpos(2) - quarkpos(1) - 1
                   thetreeamp%NumGlu(2) = quarkpos(3) - quarkpos(2) - 1
                   thetreeamp%NumGlu(3) = quarkpos(4) - quarkpos(3) - 1
                   thetreeamp%NumGlu(4) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(4)
             endif
             if( thetreeamp%NumQua .eq. 6 ) then
                   thetreeamp%NumGlu(1) = quarkpos(2) - quarkpos(1) - 1
                   thetreeamp%NumGlu(2) = quarkpos(3) - quarkpos(2) - 1
                   thetreeamp%NumGlu(3) = quarkpos(4) - quarkpos(3) - 1
                   thetreeamp%NumGlu(4) = quarkpos(5) - quarkpos(4) - 1
                   thetreeamp%NumGlu(5) = quarkpos(6) - quarkpos(5) - 1
                   thetreeamp%NumGlu(6) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(6)
             endif
 
           elseif( thetreeamp%PartType(1).eq.10 ) then ! is a gluon
             if( thetreeamp%NumQua .eq. 2 ) then
                   thetreeamp%NumGlu(1) = quarkpos(1) - 2
                   thetreeamp%NumGlu(2) = quarkpos(2) - quarkpos(1) - 1
                   thetreeamp%NumGlu(3) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(2)
             endif
             if( thetreeamp%NumQua .eq. 4 ) then
                   thetreeamp%NumGlu(1) = quarkpos(1) - 2
                   thetreeamp%NumGlu(2) = quarkpos(2) - quarkpos(1) - 1
                   thetreeamp%NumGlu(3) = quarkpos(3) - quarkpos(2) - 1
                   thetreeamp%NumGlu(4) = quarkpos(4) - quarkpos(3) - 1
                   thetreeamp%NumGlu(5) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(4)
             endif
             if( thetreeamp%NumQua .eq. 6 ) then
                   thetreeamp%NumGlu(1) = quarkpos(1) - 2
                   thetreeamp%NumGlu(2) = quarkpos(2) - quarkpos(1) - 1
                   thetreeamp%NumGlu(3) = quarkpos(3) - quarkpos(2) - 1
                   thetreeamp%NumGlu(4) = quarkpos(4) - quarkpos(3) - 1
                   thetreeamp%NumGlu(5) = quarkpos(5) - quarkpos(4) - 1
                   thetreeamp%NumGlu(6) = quarkpos(6) - quarkpos(5) - 1
                   thetreeamp%NumGlu(7) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(6)
             endif
           endif
 
 
  ig=0; iq=0; ib=0;
  do ipart=1,thetreeamp%NumPart
     partref = thetreeamp%PartRef(ipart)
     parttype= theparticles(partref)%PartType
     if( parttype.eq.glu_ ) then
           ig=ig+1
           thetreeamp%Gluons(ig)%PartType => theparticles(partref)%PartType
           thetreeamp%Gluons(ig)%ExtRef   => theparticles(partref)%ExtRef
           thetreeamp%Gluons(ig)%Mass     => theparticles(partref)%Mass
           thetreeamp%Gluons(ig)%Mass2    => theparticles(partref)%Mass2
           thetreeamp%Gluons(ig)%Helicity => theparticles(partref)%Helicity
           thetreeamp%Gluons(ig)%Mom      => theparticles(partref)%Mom
           thetreeamp%Gluons(ig)%Pol      => theparticles(partref)%Pol
           if( parttype.ne.theparticles(partref)%PartType ) print *,"Error1 in LinkTreeParticles"
     elseif( isaquark(parttype) ) then ! PartType==Quark
           iq=iq+1
           thetreeamp%Quarks(iq)%PartType => theparticles(partref)%PartType
           thetreeamp%Quarks(iq)%ExtRef   => theparticles(partref)%ExtRef
           thetreeamp%Quarks(iq)%Mass     => theparticles(partref)%Mass
           thetreeamp%Quarks(iq)%Mass2    => theparticles(partref)%Mass2
           thetreeamp%Quarks(iq)%Helicity => theparticles(partref)%Helicity
           thetreeamp%Quarks(iq)%Mom      => theparticles(partref)%Mom
           thetreeamp%Quarks(iq)%Pol      => theparticles(partref)%Pol
           if( parttype.ne.theparticles(partref)%PartType ) print *,"Error2 in LinkTreeParticles"
     elseif( isaboson(parttype) ) then  ! PartType==Boson
           ib=ib+1
           if( ib.ge.2 ) print *, "Too many bosons in LinkTreeParticles"
           thetreeamp%Boson%PartType => theparticles(partref)%PartType
           thetreeamp%Boson%ExtRef   => theparticles(partref)%ExtRef
           thetreeamp%Boson%Mass     => theparticles(partref)%Mass
           thetreeamp%Boson%Mass2    => theparticles(partref)%Mass2
           thetreeamp%Boson%Helicity => theparticles(partref)%Helicity
           thetreeamp%Boson%Mom      => theparticles(partref)%Mom
           thetreeamp%Boson%Pol      => theparticles(partref)%Pol
           if( parttype.ne.theparticles(partref)%PartType ) print *,"Error2 in LinkTreeParticles"
     endif
  enddo
  if( ig.ne.thetreeamp%NumGlu(0) .OR. iq.ne.thetreeamp%NumQua+thetreeamp%NumSca .OR. ib.ne.thetreeamp%NumPart-thetreeamp%NumGlu(0)-thetreeamp%NumQua-thetreeamp%NumSca) print *,"Error3 in LinkTreeParticles"
 
 
RETURN
END SUBROUTINE
 
 
 
!SUBROUTINE UnLinkTreeParticles(TheTreeAmp)
!implicit none
!type(TreeProcess) :: TheTreeAmp
!integer :: iPart,PartType,PartRef,ig,iq,ib
 
!  ig=0; iq=0; ib=0;
!  do iPart=1,TheTreeAmp%NumPart
!     PartRef = TheTreeAmp%PartRef(iPart)
!     PartType = TheParticles(PartRef)%PartType
 
!     if( PartType.eq.Glu_ ) then
!           ig=ig+1
!           nullify(TheTreeAmp%Gluons(ig)%PartType)
!           nullify(TheTreeAmp%Gluons(ig)%ExtRef)
!           nullify(TheTreeAmp%Gluons(ig)%Mass)
!           nullify(TheTreeAmp%Gluons(ig)%Mass2)
!           nullify(TheTreeAmp%Gluons(ig)%Helicity)
!           nullify(TheTreeAmp%Gluons(ig)%Mom)
!           nullify(TheTreeAmp%Gluons(ig)%Pol)
!     elseif( IsAQuark(PartType) ) then ! PartType==Quark
!           iq=iq+1
!           nullify(TheTreeAmp%Quarks(iq)%PartType)
!           nullify(TheTreeAmp%Quarks(iq)%ExtRef)
!           nullify(TheTreeAmp%Quarks(iq)%Mass)
!           nullify(TheTreeAmp%Quarks(iq)%Mass2)
!           nullify(TheTreeAmp%Quarks(iq)%Helicity)
!           nullify(TheTreeAmp%Quarks(iq)%Mom)
!           nullify(TheTreeAmp%Quarks(iq)%Pol)
!     elseif( IsABoson(PartType) ) then  ! PartType==Boson
!           ib=ib+1
!           nullify(TheTreeAmp%Boson%PartType)
!           nullify(TheTreeAmp%Boson%ExtRef)
!           nullify(TheTreeAmp%Boson%Mass)
!           nullify(TheTreeAmp%Boson%Mass2)
!           nullify(TheTreeAmp%Boson%Helicity)
!           nullify(TheTreeAmp%Boson%Mom)
!           nullify(TheTreeAmp%Boson%Pol)
!     endif
!  enddo
 
!RETURN
!END SUBROUTINE
 
 
 
 
! ----------------------------------------------------
 
 
FUNCTION cur_f_2fv(Gluons,Quark,Quark1PartType,Boson,NumGlu) result(Res)           ! Quarks(:) DOES include the OFF-shell quark, in contrast to all other routines!
implicit none
complex(8) :: res(1:ds)
integer :: numglu(0:2),i,rin,rout,quark1parttype
type(ptrtoparticle) :: gluons(1:),boson,quark(2:2)
complex(8) :: glumom(1:dv,numglu(0)), quarkmom(1:dv)
complex(8) :: glupol(1:dv,numglu(0)), quarkpol(1:ds)
 
 
 
   do i=1,numglu(0)
    glumom(1:dv,i) = gluons(i)%Mom(1:dv)
    glupol(1:dv,i) = gluons(i)%Pol(1:dv)
   enddo
   quarkmom(1:dv) = quark(2)%Mom(1:dv)
   quarkpol(1:ds) = quark(2)%Pol(1:ds)
 
   if( quark1parttype.ne.-quark(2)%PartType ) print *, "Wrong quark flavors in cur_f_2fV"
   if( numglu(0)-numglu(1)-numglu(2).ne.0 ) print *, "Wrong NumGlu in cur_f_2fV",numglu(0)-numglu(1)-numglu(2)
 
   rin =1
   rout=numglu(0)
   if( quark(2)%PartType .gt.0 ) then      !    X----->----
      res(:) = fv(glupol(1:dv,rin:rout),glumom(1:dv,rin:rout),quarkpol(1:ds),quarkmom(1:dv),quark(2)%Mass,quark1parttype,boson%Pol(1:dv),boson%Mom(1:dv),numglu(1))
   else
      res(:) = bfv(glupol(1:dv,rin:rout),glumom(1:dv,rin:rout),quarkpol(1:ds),quarkmom(1:dv),quark(2)%Mass,quark1parttype,boson%Pol(1:dv),boson%Mom(1:dv),numglu(1))
   endif
 
 
return
END FUNCTION
 
 
 
 
 
 
 
      recursive function fv(e,k,sp,p,mass,QuarkFlavor,eV,kV,ms) result(res)
      implicit none
      complex(8), intent(in) :: e(:,:), k(:,:)
      complex(8), intent(in) :: sp(:), p(:)
      complex(8), intent(in) :: ev(:), kv(:)
      integer, intent(in) ::  ms,quarkflavor
      integer             :: ms1,m,ng1, ng2
      integer :: ngluon
      complex(8)             :: res(size(sp))
      complex(8)             :: tmp(size(sp))
      complex(8)             :: k1(size(p))
      complex(8)             :: k2(size(p))
      complex(8)             :: sp2(size(sp))
      complex(8)             :: sp3(size(sp))
      complex(8)             :: e1(size(e,dim=1))
      complex(8)             :: e2(size(e,dim=1))
      complex(8)  :: k1sq,k2sq,k3sq
      real(8) :: mass
      complex(8) :: couplvqq_left,couplvqq_right,couplvqq_left2,couplvqq_right2
      character,parameter :: ferfla*3="dum" ! dummy, only used for check of flavor consistency inside the functions f,bf
 
 
      ngluon = size(e,dim=2)
      ng1 = ms   !#gluons to the left of a f-line
      ng2 = ngluon - ms  !#gluons to the right of the f-line
      if (ng2 < 0) write(*,*) 'WRONG DEFINITION OF CURRENT fV'
 
      if( abs(quarkflavor).eq.top_ .or. abs(quarkflavor).eq.bot_ ) then!   note that Bot_ is treated as top quark in TOPAZ!
         couplvqq_left  = couplhtt_left_dyn
         couplvqq_right = couplhtt_right_dyn
      else
         couplvqq_left=0d0
         couplvqq_right=0d0
         print *, "this should not happen for the Higgs",quarkflavor,top_
      endif
 
 
 
if (ngluon == 0) then
         res = vbqv(sp,ev,couplvqq_left,couplvqq_right)
else
 
       res = (0d0,0d0)
       do m=0,ng2-1
           k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
           e1=g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
           k1sq=sc_(k1,k1)
 
           k2 = sum(k(:,1:ng1+m),dim=2)
           k2 = k2 + p + kv
           k2sq = sc_(k2,k2)-mass**2
           sp2 = fv(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,quarkflavor,ev,kv,ng1)
           sp2 = spb2_(sp2,k2)+mass*sp2
 
           tmp = vqg(sp2,e1)
           if (m < ng2-1)  then
              if(abs(k1sq) > propcut) then
                  tmp = -(0d0,1d0)/k1sq*tmp
              else
                  tmp = (0d0,0d0)
              endif
           endif
 
           if (abs(k2sq) > propcut) then
                  tmp =  (0d0,1d0)/k2sq*tmp
           else
                  tmp = (0d0,0d0)
           endif
           res = res + tmp
        enddo
 
 
 
 
        do m=1,ng1
           k1 = sum(k(:,1:m),dim=2)
           e1=g(e(:,1:m),k(:,1:m))
           k1sq = sc_(k1,k1)
 
           k2 = sum(k(:,m+1:ngluon),dim=2)
           k2 = k2 + p + kv
           k2sq = sc_(k2,k2) - mass**2
           ms1 = ng1 - m
           sp2=fv(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,quarkflavor,ev,kv,ms1)
 
           sp2 = spb2_(sp2,k2)+mass*sp2
 
           tmp = vgq(e1,sp2)
           if (m > 1) then
              if (abs(k1sq) > propcut) then
                  tmp=-(0d0,1d0)/k1sq*tmp
              else
                  tmp = (0d0,0d0)
              endif
           endif
 
           if (abs(k2sq) > propcut) then
              tmp=(0d0,1d0)/k2sq*tmp
           else
              tmp = (0d0,0d0)
           endif
 
           res = res + tmp
        enddo
 
 
        sp2 = f(e,k,sp,p,mass,ferfla,ferfla,ms)
        k2 = sum(k(:,1:ngluon),dim=2)
        k2 = k2 + p
        k2sq = sc_(k2,k2)  - mass**2
 
        sp2 = spb2_(sp2,k2)+ mass*sp2
 
        tmp = vbqv(sp2,ev,couplvqq_left,couplvqq_right)
        if (abs(k2sq) > propcut) then
               tmp = (0d0,1d0)/k2sq*tmp
        else
                tmp = (0d0,0d0)
        endif
        res = res + tmp
 
endif
 
end function fv
 
 
 
 
      recursive function bfv(e,k,sp,p,mass,QuarkFlavor,eV,kV,ms) result(res)
      implicit none
      complex(8), intent(in) :: e(:,:), k(:,:)
      complex(8), intent(in) :: sp(:), p(:)
      complex(8), intent(in) :: ev(:), kv(:)
      integer, intent(in) ::  ms,quarkflavor
      integer             :: ms1,m,ng1, ng2
      integer :: ngluon
      complex(8)             :: res(size(sp))
      complex(8)             :: tmp(size(sp))
      complex(8)             :: k1(size(p))
      complex(8)             :: k2(size(p))
      complex(8)             :: sp2(size(sp))
      complex(8)             :: sp3(size(sp))
      complex(8)             :: e1(size(e,dim=1))
      complex(8)             :: e2(size(e,dim=1))
      complex(8)  :: k1sq,k2sq,k3sq
      real(8) :: mass
      complex(8) :: couplvqq_left,couplvqq_right,couplvqq_left2,couplvqq_right2
      character,parameter :: ferfla*3="dum" ! dummy, only used for check of flavor consistency inside the functions f,bf
 
 
      ngluon = size(e,dim=2)
      ng1 = ms   !#gluons to the left of a f-line
      ng2 = ngluon - ms  !#gluons to the right of the f-line
 
      if (ng2 < 0) write(*,*) 'WRONG DEFINITION OF CURRENT fbV'
 
      if( abs(quarkflavor).eq.top_ .or. abs(quarkflavor).eq.bot_ ) then!   note that Bot_ is treated as top quark in TOPAZ!
          couplvqq_left  = couplhtt_left_dyn
          couplvqq_right = couplhtt_right_dyn
      else
          couplvqq_left=0d0
          couplvqq_right=0d0
          print *, "this should not happen for the Higgs",quarkflavor,top_
      endif
 
 
 
if (ngluon == 0) then
         res = vvq(ev,sp,couplvqq_left,couplvqq_right)
else
 
       res = (0d0,0d0)
       do m=0,ng2-1
           k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
           e1=g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
           k1sq=sc_(k1,k1)
 
           k2 = sum(k(:,1:ng1+m),dim=2)
           k2 = -k2 - p - kv
           k2sq = sc_(k2,k2)-mass**2
 
           sp2 = bfv(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,quarkflavor,ev,kv,ng1)
           sp2 = spi2_(k2,sp2)+mass*sp2
 
           tmp = vbqg(sp2,e1)
 
           if (m < ng2-1) then
            if (abs(k1sq) > propcut) then
                tmp = -(0d0,1d0)/k1sq*tmp
            else
                tmp = (0d0,0d0)
            endif
           endif
 
           if (abs(k2sq) > propcut) then
              tmp =  (0d0,1d0)/k2sq*tmp
           else
               tmp = (0d0,0d0)
           endif
 
           res = res + tmp
        enddo
 
 
        do m=1,ng1
 
           k1 = sum(k(:,1:m),dim=2)
           e1=g(e(:,1:m),k(:,1:m))
           k1sq = sc_(k1,k1)
 
           k2 = sum(k(:,m+1:ngluon),dim=2)
           k2 = -k2 - p - kv
           k2sq = sc_(k2,k2) - mass**2
           ms1 = ng1 - m
           sp2=bfv(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,quarkflavor,ev,kv,ms1)
 
           sp2 = spi2_(k2,sp2)+mass*sp2
 
           tmp = vgbq(e1,sp2)
 
           if (m > 1) then
              if (abs(k1sq) > propcut) then
                  tmp=-(0d0,1d0)/k1sq*tmp
              else
                  tmp = (0d0,0d0)
              endif
           endif
 
           if (abs(k2sq) > propcut) then
                 tmp=(0d0,1d0)/k2sq*tmp
              else
                 tmp = (0d0,0d0)
           endif
 
           res = res + tmp
 
        enddo
 
 
 
        sp2 = bf(e,k,sp,p,mass,ferfla,ferfla,ms)
        k2 = sum(k(:,1:ngluon),dim=2)
        k2 = -k2 - p
        k2sq = sc_(k2,k2)   -mass**2
 
!        sp2 = spb2_(sp2,k2) +mass*sp2
        sp2 = spi2_(k2,sp2) +mass*sp2
 
        tmp = vvq(ev,sp2,couplvqq_left,couplvqq_right)
        if (abs(k2sq) > propcut) then
               tmp = (0d0,1d0)/k2sq*tmp
        else
               tmp = (0d0,0d0)
        endif
        res = res + tmp
 
endif
 
 
 
 
      end function bfv
 
 
 
 
 
FUNCTION cur_f_4fv(Gluons,Quarks,Quark1PartType,Boson,BosonVertex,NumGlu,tag_f,tag_Z) result(res)           ! Quarks(:) does not include the OFF-shell quark
implicit none
integer :: numglu(0:4),quark1parttype
type(ptrtoparticle) :: gluons(1:),quarks(2:4),boson
integer :: tag_f,bosonvertex,tag_z
integer,target :: tmpextref
complex(8) :: res(1:ds),tmp(1:ds)
complex(8) :: ubar1(1:ds)
complex(8),target :: ubar0(1:ds)
complex(8) :: eps1(1:dv)
complex(8) :: eps2(1:dv)
type(ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(4)),tmpquark(1:1)
complex(8) :: propfac1,propfac2
complex(8),target :: pmom1(1:dv)
complex(8) :: pmom2(1:dv)
integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b
integer :: rin,rout,i,counter
 
 
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-numglu(1)-numglu(2)-numglu(3)-numglu(4).ne.0 ) print *, "wrong number of gluons in cur_f_4f"
    if(quarks(3)%PartType.eq.-quarks(4)%PartType .and. quark1parttype.ne.-quarks(2)%PartType ) print *,"wrong flavor in cur_f_4f (1)"
    if(quarks(2)%PartType.eq.-quarks(3)%PartType .and. quark1parttype.ne.-quarks(4)%PartType ) print *,"wrong flavor in cur_f_4f (2)"
!DEC$ ENDIF
    res(:)=(0d0,0d0)
   if ( quark1parttype.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType ) then
      do n2a=0,numglu(2)
         do n4a=0,numglu(4)
            n2b = numglu(2)-n2a
            n4b = numglu(4)-n4a
            rin =numglu(1)+n2a+1
            rout=numglu(1)+numglu(2)+numglu(3)+n4a
            if (bosonvertex .eq. 1 .or. bosonvertex .eq. 2 .or. bosonvertex .eq. 4) then
               eps2 = cur_g_2f(gluons(rin:rout),quarks(3:4),(/1+n2b+numglu(3)+n4a,n2b,numglu(3),n4a/))
               pmom1(:) = summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom
            elseif (bosonvertex .eq. 3) then
               eps2 = cur_g_2fv(gluons(rin:rout),quarks(3:4),boson,(/1+n2b+numglu(3)+n4a,n2b,numglu(3),n4a/))
               pmom1(:) = summom(gluons,rin,rout) + quarks(3)%Mom(:) + quarks(4)%Mom(:) + boson%Mom(:)
            endif
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
            do n1a=0,numglu(1)
               n1b = numglu(1)-n1a
               ! Fer2
               if (bosonvertex.eq.1 .or. bosonvertex.eq.2) then
                  rin =n1a+1
                  rout=numglu(1)+n2a
                  ubar1(:) = cur_f_2fv(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,boson,(/n2a+n1b,n1b,n2a/) )
                  pmom2(:) = quarks(2)%Mom(:) + summom(gluons,rin,rout) + boson%Mom(:)
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
 
                  if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) then
                     propfac2=(0d0,0d0)
                  endif
 
                  if( quarks(2)%PartType.lt.0 ) then
                     ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
                  else
                     ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
                  endif
 
                  if( quarks(2)%PartType.lt.0 ) then
                     ubar0(:) = vbqg(ubar1,eps2)
                  else
                     ubar0(:) = vqg(ubar1,eps2)
                  endif
 
                  pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a) + boson%Mom(:)
                  if(n1a.ge.1 .or. n4b.ge.1) then
                     propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
                     if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) then
                        propfac1=(0d0,0d0)
                     endif
                     if( quarks(2)%PartType.lt.0 ) then
                        ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
                     else
                        ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
                     endif
                  endif
 
                  tmpquark(1)%Mom  => pmom1(:)
                  tmpquark(1)%Pol  => ubar0(:)
                  tmpquark(1)%Mass => quarks(2)%Mass
                  tmpquark(1)%Mass2=> quarks(2)%Mass2
                  tmpextref = -1
                  tmpquark(1)%ExtRef => tmpextref
                  tmpquark(1)%PartType => quarks(2)%PartType
                  counter=1
                  rin =1
                  rout=n1a
                  do i=rin,rout
                     call copyparticleptr(gluons(i),tmpgluons(counter))
                     counter=counter+1
                  enddo
                  rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout=numglu(0)
                  do i=rin,rout
                     call copyparticleptr(gluons(i),tmpgluons(counter))
                     counter=counter+1
                  enddo
                  tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
                  res(:) = res(:) + tmp(:)
                endif
 
                ! Fer 3
                if (bosonvertex .eq. 1 .or. bosonvertex .eq. 4) then
                   rin =n1a+1
                   rout=numglu(1)+n2a
                   ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2a+n1b,n1b,n2a/) )
                   if(n1b.ge.1 .or. n2a.ge.1) then
                      pmom2(:) = quarks(2)%Mom + summom(gluons,rin,rout)  ! can be moved outside the n1a-loop
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
                      if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) cycle
                      if( quarks(2)%PartType.lt.0 ) then
                         ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
                      else
                         ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
                      endif
                   endif
                   if( quarks(2)%PartType.lt.0 ) then
                      ubar0(:) = vbqg(ubar1,eps2)
                   else
                      ubar0(:) = vqg(ubar1,eps2)
                   endif
 
                   pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
                   propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
                   if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) then
                      propfac1=(0d0,0d0)
                   endif
                   if( quarks(2)%PartType.lt.0 ) then
                      ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
                   else
                      ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
                   endif
 
                   tmpquark(1)%Mom  => pmom1(:)
                   tmpquark(1)%Pol  => ubar0(:)
                   tmpquark(1)%Mass => quarks(2)%Mass
                   tmpquark(1)%Mass2=> quarks(2)%Mass2
                   tmpextref = -1
                   tmpquark(1)%ExtRef => tmpextref
                   tmpquark(1)%PartType => quarks(2)%PartType
                   counter=1
                   rin =1
                   rout=n1a
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout=numglu(0)
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   tmp(:) = cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,boson,(/counter-1,n1a,n4b/) )
                   res(:) = res(:) + tmp(:)
                endif
 
                if (bosonvertex .eq. 3) then
                   rin =n1a+1
                   rout=numglu(1)+n2a
                   ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2a+n1b,n1b,n2a/) )
                   if(n1b.ge.1 .or. n2a.ge.1) then
                      pmom2(:) = quarks(2)%Mom + summom(gluons,rin,rout)  ! can be moved outside the n1a-loop
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
                      if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) then
                         propfac2=(0d0,0d0)
                      endif
                      if( quarks(2)%PartType.lt.0 ) then
                         ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
                      else
                         ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
                      endif
                   endif
                   if( quarks(2)%PartType.lt.0 ) then
                      ubar0(:) = vbqg(ubar1,eps2)
                   else
                      ubar0(:) = vqg(ubar1,eps2)
                   endif
 
                   pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)+boson%Mom
                   if(n1a.ge.1 .or. n4b.ge.1) then
                      propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
                      if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) cycle
                      if( quarks(2)%PartType.lt.0 ) then
                         ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
                      else
                         ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
                      endif
                   endif
 
                   tmpquark(1)%Mom  => pmom1(:)
                   tmpquark(1)%Pol  => ubar0(:)
                   tmpquark(1)%Mass => quarks(2)%Mass
                   tmpquark(1)%Mass2=> quarks(2)%Mass2
                   tmpextref = -1
                   tmpquark(1)%ExtRef => tmpextref
                   tmpquark(1)%PartType => quarks(2)%PartType
                   counter=1
                   rin =1
                   rout=n1a
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout=numglu(0)
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
                   res(:) = res(:) + tmp(:)
                endif
             enddo
          enddo
       enddo
    endif
 
 
    if( (quark1parttype.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType)      .AND.  &
         (quarks(4)%ExtRef.ne.-1 .or. tag_f.ne.1 .or. abs(quark1parttype).ne.abs(quarks(2)%PartType)) ) then
 
       do n1a=0,numglu(1)
          do n3a=0,numglu(3)
             n1b = numglu(1)-n1a
             n3b = numglu(3)-n3a
 
             rin =n1a+1
             rout=numglu(1)+numglu(2)+n3a
 
! This means that all the ext gluons are on this lines,
! and we must remove this to prevent color issues with a Z on the quark loop, see RR notes
 
             if (bosonvertex.eq.1 .or. bosonvertex.eq.3 .or. bosonvertex.eq.4)  then
                if ( tag_z .eq. 1 .and. n1b+numglu(2)+n3a .eq. numglu(0) ) then
!                     print  * , "cycle for tag_Z=1 in cur_f_4fV"
                    cycle
                endif
                eps2 = cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
                pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
             elseif (bosonvertex .eq. 2) then
                eps2 = cur_g_2fv(gluons(rin:rout),quarks(2:3),boson,(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
                pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
             endif
             propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
             if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
             eps2 = eps2*propfac1
 
             do n4a=0,numglu(4)
                n4b = numglu(4)-n4a
                ! radiate V off Fer4
                if (bosonvertex .eq. 3 .or. bosonvertex .eq. 4) then
                   rin =numglu(1)+numglu(2)+n3a+1
                   rout=numglu(1)+numglu(2)+numglu(3)+n4a
 
                   ubar1(:) = cur_f_2fv(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,boson,(/n3b+n4a,n3b,n4a/) )
                   pmom2(:) = quarks(4)%Mom(:) + summom(gluons,rin,rout) + boson%Mom(:)
                   propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(4)%Mass2)
                   if( abs(sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.propcut ) then
                      propfac2=(0d0,0d0)
                   endif
 
                   if( quarks(4)%PartType.lt.0 ) then
                      ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
                   else
                      ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
                   endif
                   if( quarks(4)%PartType.lt.0 ) then
                      ubar0(:) = vgbq(eps2,ubar1)
                   else
                      ubar0(:) = vgq(eps2,ubar1)
                   endif
 
                   pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a) + boson%Mom(:)
                   if(n1a.ge.1 .or. n4b.ge.1) then
                      propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(4)%Mass2)
                      if( abs(sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.propcut ) then
                         propfac1=(0d0,0d0)
                      endif
                      if( quarks(4)%PartType.lt.0 ) then
                         ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(4)%Mass*ubar0(:))*propfac1
                      else
                         ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(4)%Mass*ubar0(:))*propfac1
                      endif
                   endif
 
                   tmpquark(1)%Mom  => pmom1(:)
                   tmpquark(1)%Pol  => ubar0(:)
                   tmpquark(1)%Mass => quarks(4)%Mass
                   tmpquark(1)%Mass2=> quarks(4)%Mass2
                   tmpextref = -1
                   tmpquark(1)%ExtRef => tmpextref
                   tmpquark(1)%PartType => quarks(4)%PartType
                   counter=1
                   rin =1
                   rout=n1a
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout=numglu(0)
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
                   res(:) = res(:) + tmp(:)
                endif
 
                if (bosonvertex .eq. 1 .or. bosonvertex .eq. 4) then
                   ! radiate V off Fer1
                   rin =numglu(1)+numglu(2)+n3a+1
                   rout=numglu(1)+numglu(2)+numglu(3)+n4a
                   ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/n3b+n4a,n3b,n4a/) )
                   if(n3b.ge.1 .or. n4a.ge.1) then
                      pmom2(:) = quarks(4)%Mom + summom(gluons,rin,rout)  ! can be moved outside the n1a-loop
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.propcut ) then
                         propfac2=(0d0,0d0)
                      endif
 
                      if( quarks(4)%PartType.lt.0 ) then
                         ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
                      else
                         ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
                      endif
                   endif
                   if( quarks(4)%PartType.lt.0 ) then
                      ubar0(:) = vgbq(eps2,ubar1)
                   else
                      ubar0(:) = vgq(eps2,ubar1)
                   endif
 
                   pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
                   propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(4)%Mass2)
                   if( abs(sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.propcut ) then
                      propfac1=(0d0,0d0)
                   endif
                   if( quarks(4)%PartType.lt.0 ) then
                      ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(4)%Mass*ubar0(:))*propfac1
                   else
                      ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(4)%Mass*ubar0(:))*propfac1
                   endif
 
                   tmpquark(1)%Mom  => pmom1(:)
                   tmpquark(1)%Pol  => ubar0(:)
                   tmpquark(1)%Mass => quarks(4)%Mass
                   tmpquark(1)%Mass2=> quarks(4)%Mass2
                   tmpextref = -1
                   tmpquark(1)%ExtRef => tmpextref
                   tmpquark(1)%PartType => quarks(4)%PartType
                   counter=1
                   rin =1
                   rout=n1a
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout=numglu(0)
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   tmp(:) = cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,boson,(/counter-1,n1a,n4b/) )
                   res(:) = res(:) + tmp(:)
                endif
 
 
                if (bosonvertex .eq. 2) then
                   rin =numglu(1)+numglu(2)+n3a+1
                   rout=numglu(1)+numglu(2)+numglu(3)+n4a
                   ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/n3b+n4a,n3b,n4a/) )
                   if ( n3b .ge. 1 .or. n4a .ge. 1) then
                      pmom2(:) = quarks(4)%Mom(:) + summom(gluons,rin,rout)
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(4)%Mass2)
 
                      if( abs(sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.propcut ) then
                         propfac2=(0d0,0d0)
                      endif
 
                      if( quarks(4)%PartType.lt.0 ) then
                         ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
                      else
                         ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
                      endif
                   endif
                   if( quarks(4)%PartType.lt.0 ) then
                      ubar0(:) = vgbq(eps2,ubar1)
                   else
                      ubar0(:) = vgq(eps2,ubar1)
                   endif
                   pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)+boson%Mom
                   if(n1a.ge.1 .or. n4b.ge.1) then
                      propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(4)%Mass2)
                      if( abs(sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.propcut ) then
                         cycle
                      endif
                      if( quarks(4)%PartType.lt.0 ) then
                         ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(4)%Mass*ubar0(:))*propfac1
                      else
                         ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(4)%Mass*ubar0(:))*propfac1
                      endif
                   endif
 
                   tmpquark(1)%Mom  => pmom1(:)
                   tmpquark(1)%Pol  => ubar0(:)
                   tmpquark(1)%Mass => quarks(4)%Mass
                   tmpquark(1)%Mass2=> quarks(4)%Mass2
                   tmpextref = -1
                   tmpquark(1)%ExtRef => tmpextref
                   tmpquark(1)%PartType => quarks(4)%PartType
                   counter=1
                   rin =1
                   rout=n1a
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout=numglu(0)
                   do i=rin,rout
                      call copyparticleptr(gluons(i),tmpgluons(counter))
                      counter=counter+1
                   enddo
                   tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
                   res(:) = res(:) + tmp(:)
                endif
             enddo
          enddo
       enddo
 
    endif
 
    return
  END FUNCTION cur_f_4fv
 
 
 
 
 
 
FUNCTION cur_f_4f(Gluons,Quarks,Quark1PartType,NumGlu,tag_f,tag_Z_arg) result(res)           ! Quarks(:) does not include the OFF-shell quark
implicit none
integer :: numglu(0:4),quark1parttype
type(ptrtoparticle) :: gluons(1:),quarks(2:4)
integer :: tag_f,tag_z
integer, optional :: tag_z_arg
integer,target :: tmpextref
complex(8) :: res(1:ds),tmp(1:ds)
complex(8) :: ubar1(1:ds)
complex(8),target :: ubar0(1:ds)
complex(8) :: eps1(1:dv)
complex(8) :: eps2(1:dv)
type(ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(4)),tmpquark(1:1)
complex(8) :: propfac1,propfac2
complex(8),target :: pmom1(1:dv)
complex(8) :: pmom2(1:dv)
integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b
integer :: rin,rout,i,counter
 
 
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-numglu(1)-numglu(2)-numglu(3)-numglu(4).ne.0 ) print *, "wrong number of gluons in cur_f_4f"
    if(quarks(3)%PartType.eq.-quarks(4)%PartType .and. quark1parttype.ne.-quarks(2)%PartType ) print *,"wrong flavor in cur_f_4f (1)"
    if(quarks(2)%PartType.eq.-quarks(3)%PartType .and. quark1parttype.ne.-quarks(4)%PartType ) print *,"wrong flavor in cur_f_4f (2)"
!DEC$ ENDIF
 
   res(:)=(0d0,0d0)
 
   tag_z=0
   if( present(tag_z_arg) ) then
      tag_z=tag_z_arg
   endif
 
   if( quark1parttype.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType ) then
!     (I)
      do n2a=0,numglu(2)
      do n4a=0,numglu(4)
         n2b = numglu(2)-n2a
         n4b = numglu(4)-n4a
 
         rin =numglu(1)+n2a+1
         rout=numglu(1)+numglu(2)+numglu(3)+n4a
         eps2 = cur_g_2f(gluons(rin:rout),quarks(3:4),(/1+n2b+numglu(3)+n4a,n2b,numglu(3),n4a/))
         pmom1(:) = summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom
         propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
         if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
         eps2 = eps2*propfac1
         do n1a=0,numglu(1)
            n1b = numglu(1)-n1a
            ! Fer2
            rin =n1a+1
            rout=numglu(1)+n2a
            ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2a+n1b,n1b,n2a/) )
            if(n1b.ge.1 .or. n2a.ge.1) then
               pmom2(:) = quarks(2)%Mom + summom(gluons,rin,rout)  ! can be moved outside the n1a-loop
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
               if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
               else
                  ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
               endif
            endif
            if( quarks(2)%PartType.lt.0 ) then
               ubar0(:) = vbqg(ubar1,eps2)       ! re-checked
            else
               ubar0(:) = vqg(ubar1,eps2)        ! re-checked
            endif
 
            pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
            if(n1a.ge.1 .or. n4b.ge.1) then
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
               if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
               else
                  ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
               endif
            endif
 
            tmpquark(1)%Mom  => pmom1(:)
            tmpquark(1)%Pol  => ubar0(:)
            tmpquark(1)%Mass => quarks(2)%Mass
            tmpquark(1)%Mass2=> quarks(2)%Mass2
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpquark(1)%PartType => quarks(2)%PartType
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
            res(:) = res(:) + tmp(:)
         enddo
      enddo
      enddo
   endif
 
 
!    if( (Quark1PartType.eq.-Quarks(4)%PartType .and. Quarks(2)%PartType.eq.-Quarks(3)%PartType) .AND.  &
!        ((abs(Quark1PartType).ne.abs(Quarks(2)%PartType).and.(tag_f.ne.3)) .or. &
!         (abs(Quark1PartType).eq.abs(Quarks(2)%PartType).and.(tag_f.ne.1.and.tag_f.ne.3))) ) then
!    if( (Quark1PartType.eq.-Quarks(4)%PartType .and. Quarks(2)%PartType.eq.-Quarks(3)%PartType) .AND.  &
!        .not.(Quarks(4)%ExtRef.eq.-1 .and. tag_f.eq.1 .and. abs(Quark1PartType).eq.(Quarks(2)%PartType))  ) then
   if( (quark1parttype.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType) .AND.  &
       (quarks(4)%ExtRef.ne.-1 .or. tag_f.ne.1 .or. abs(quark1parttype).ne.abs(quarks(2)%PartType)) &
     ) then
!     (II)
      do n1a=0,numglu(1)
      do n3b=0,numglu(3)
         n1b = numglu(1)-n1a
         n3a = numglu(3)-n3b
         rin =n1a+1
         rout=numglu(1)+numglu(2)+n3a
 
 
! This prevents color issues with a Z on the quark loop, see RR notes
         if ( tag_z .eq. 1 .and. (n1b+numglu(2)+n3a .eq. numglu(0)) .and. quarks(4)%ExtRef .eq. -1) then
!             print  * , "cycle for tag_Z=1 in cur_f_4f",n1a,n1b
            cycle
         endif
 
 
 
 
         eps2 = cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
         pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
         propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
         if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
         eps2 = eps2*propfac1
 
         do n4a=0,numglu(4)
            n4b = numglu(4)-n4a
            ! Fer4
            rin =numglu(1)+numglu(2)+n3a+1
            rout=numglu(1)+numglu(2)+numglu(3)+n4a
            ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/numglu(3)+n4a-n3a,n3b,n4a/) )
            if(n3b.ge.1 .or. n4a.ge.1) then
               pmom2(:) = quarks(4)%Mom + summom(gluons,rin,rout)
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(4)%Mass2)
               if( abs(sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.propcut ) cycle
               if( quarks(4)%PartType.lt.0 ) then
                  ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
               else
                  ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
               endif
            endif
            if( quarks(4)%PartType.lt.0 ) then
               ubar0(:) = vgbq(eps2,ubar1)   !! changed from vqg(ubar1,eps2)       ! re-checked
            else
               ubar0(:) = vgq(eps2,ubar1)   !! changed from vbqg(ubar1,eps2)       ! re-checked
            endif
 
            pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
            if(n1a.ge.1 .or. n4b.ge.1) then
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(4)%Mass2)
               if( abs(sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.propcut ) cycle
               if( quarks(4)%PartType.lt.0 ) then
                  ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(4)%Mass*ubar0(:))*propfac1
               else
                  ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(4)%Mass*ubar0(:))*propfac1
               endif
            endif
            tmpquark(1)%Mom  => pmom1(:)
            tmpquark(1)%Pol  => ubar0(:)
            tmpquark(1)%Mass => quarks(4)%Mass
            tmpquark(1)%Mass2=> quarks(4)%Mass2
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpquark(1)%PartType => quarks(4)%PartType
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
            res(:) = res(:) + tmp(:)
         enddo
      enddo
      enddo
   endif
return
END FUNCTION
 
 
 
 
      recursive function g(e,k) result(res)
      implicit none
      complex(8), intent(in) :: e(:,:),k(:,:)
      complex(8)             :: e1(size(e,dim=1))
      complex(8)             :: res(size(e,dim=1))
      complex(8) :: k1(size(e,dim=1))
      complex(8) :: k2(size(e,dim=1))
      complex(8) :: k3(size(e,dim=1))
      complex(8) :: e2(size(e,dim=1))
      complex(8) :: e3(size(e,dim=1))
      complex(8) :: tmp(size(e,dim=1))
      complex(8) :: k1sq, k2sq, k3sq
      integer :: npart, m, m1
 
      npart = size(e,dim=2)
 
      if (npart == 1) then
         res = e(:,1)
 
      elseif (npart == 2) then
         res = vggg(e(:,1),k(:,1),e(:,2),k(:,2))
        else
 
           res = (0d0,0d0)
 
           do m=1,npart-1
              k1=sum(k(:,1:m),dim=2)
              k2=sum(k(:,m+1:npart),dim=2)
 
              e1 = g(e(:,1:m),k(:,1:m))
              e2 = g(e(:,m+1:npart),k(:,m+1:npart))
 
              tmp = vggg(e1,k1,e2,k2)
 
              if (m > 1) then
              k1sq = sc_(k1,k1)
               if (abs(k1sq) > propcut) then
              tmp = -(0d0,1d0)*tmp/k1sq
               else
                  tmp = (0d0,0d0)
               endif
              endif
 
              if (m + 1 < npart) then
              k2sq = sc_(k2,k2)
                if (abs(k2sq) > propcut) then
              tmp = -(0d0,1d0)*tmp/k2sq
                else
              tmp = (0d0,0d0)
                endif
              endif
 
              res = res + tmp
 
      if (m <= npart-2) then
 
        do m1=m+1,npart-1
!        e1 = g(e(:,1:m),k(:,1:m))  ! e1 is already computed
        e2=g(e(:,m+1:m1),k(:,m+1:m1))
        e3=g(e(:,m1+1:npart),k(:,m1+1:npart))
        k2 = sum(k(:,m+1:m1),dim=2)
        k3 = sum(k(:,m1+1:npart),dim=2)
        tmp = vgggg(e1,e2,e3)
        if (m > 1) then
        k1sq = sc_(k1,k1)
          if(abs(k1sq) > propcut) then
        tmp = -(0d0,1d0)*tmp/k1sq
         else
            tmp = (0d0,0d0)
          endif
        endif
        if (m+1 < m1) then
         k2sq = sc_(k2,k2)
           if(abs(k2sq) > propcut) then
           tmp = -(0d0,1d0)*tmp/k2sq
           else
            tmp = (0d0,0d0)
           endif
        endif
        if (m1+1 < npart) then
           k3sq = sc_(k3,k3)
           if(abs(k3sq) > propcut) then
           tmp = -(0d0,1d0)*tmp/k3sq
           else
           tmp = (0d0,0d0)
           endif
        endif
        res = res + tmp
        enddo
 
        endif
        enddo
        endif
 
      end function
 
 
 
 
 
FUNCTION cur_g(Gluons,NumGlu)!  note the off-shell gluon has be counted in NumGlu
implicit none
type(ptrtoparticle) :: gluons(1:)
complex(8) :: cur_g(1:dv)
integer :: ngluons,numglu
complex(8) :: glu_subcur(1:dv,1:36)  ! max. 8 gluons allowed
complex(8) :: mom_sum(1:dv,1:36), propfactor,propdenom
integer :: ind0,ind1,ind2,ind3,j,l,mu
integer :: a,b,i1,i2
 
!DEC$ IF (_DebugWriteCurrents==1)
character :: parts(20)*4
integer :: parti(20)
 
    do i1=1,numglu-1
       if( gluons(i1)%ExtRef.eq.-1 ) then
          exit
       else
          parti(i1)=gluons(i1)%ExtRef
       endif
       if(i1.eq.numglu-1) print *, parti(1:numglu-1)
!       write (parts(1:20), '(I20)') parti(1:NumGlu-1)
    enddo
!DEC$ ENDIF
 
   ngluons = numglu-1
 
   do a=0,ngluons-1
      do b=1,ngluons-a
 
         i1 = b
         i2 = a+b
         ind0 = linear_map(i1,i2,ngluons)
 
         if (i1.eq.i2) then
            glu_subcur(1:dv,ind0) = gluons(i1)%Pol(1:dv)
            mom_sum(1:dv,ind0)    = gluons(i1)%Mom(1:dv)
         else
            mom_sum(1:dv,ind0) = mom_sum(1:dv,ind0-ngluons+i2-i1-1) + gluons(i2)%Mom(1:dv)
            if ( i1 .ne. 1 .or. i2 .ne. ngluons ) then
               propdenom = mom_sum(1:dv,ind0).ndot.mom_sum(1:dv,ind0)
               if( abs(propdenom).lt.propcut ) cycle
               propfactor = (0d0,-1d0)/propdenom
            else
               propfactor = 1d0
            endif
            do mu=1,dv
               glu_subcur(mu,ind0) = 0d0
            enddo
            do j=i1,i2-1
               ind1 = linear_map(i1,j,ngluons)
               ind2 = linear_map(j+1,i2,ngluons)
               glu_subcur(1:dv,ind0) = glu_subcur(1:dv,ind0) +            &
                   eval_tripvert( mom_sum(1:dv,ind1),mom_sum(1:dv,ind2),  &
                                   glu_subcur(1:dv,ind1),glu_subcur(1:dv,ind2)) * propfactor
            enddo
            do j=i1,i2-2
               do l=j+1,i2-1
                  ind1 = linear_map(i1,j,ngluons)
                  ind2 = linear_map(j+1,l,ngluons)
                  ind3 = linear_map(l+1,i2,ngluons)
                  glu_subcur(1:dv,ind0) = glu_subcur(1:dv,ind0) +  &
                      eval_quadvert(                                       &
                        glu_subcur(1:dv,ind1),glu_subcur(1:dv,ind2),glu_subcur(1:dv,ind3)) * propfactor
               enddo
            enddo
         endif
      enddo
   enddo
   cur_g(1:dv) = glu_subcur(1:dv,ind0)
 
return
END FUNCTION
 
 
 
 
 
      recursive function f(e,k,sp,p,mass,flout,flin,ms) result(res)
      implicit none
      complex(8), intent(in) :: e(:,:), k(:,:)
      complex(8), intent(in) :: sp(:), p(:)
      character,  intent(in) :: flin*3  ! flavor of off-shell f-line
      character,  intent(in)  :: flout*3 ! flavor of on-shell f-line
      integer, intent(in) ::  ms
      integer             :: ms1,m,ng1, ng2
      integer :: ngluon
      complex(8)             :: res(size(sp))
      complex(8)             :: tmp(size(sp))
      complex(8)             :: k1(size(p))
      complex(8)             :: k2(size(p))
      complex(8)             :: sp2(size(sp))
      complex(8)             :: sp3(size(sp))
      complex(8)             :: e1(size(e,dim=1))
      complex(8)             :: e2(size(e,dim=1))
      complex(8)  :: k1sq,k2sq,k3sq
      real(8) :: mass
 
      ngluon = size(e,dim=2)
      ng1 = ms   !#gluons to the left of a f-line
      ng2 = ngluon - ms  !#gluons to the right of the f-line
 
      if (flout.ne.flin) then
         res = (0d0,0d0)
      else
 
      if (ng2 < 0) write(*,*) 'WRONG DEFINITION OF CURRENT A'
 
      if (ngluon == 0) then
         res = sp
 
      else
 
       res = (0d0,0d0)
       do m=0,ng2-1
           k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
           e1=g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
           k1sq=sc_(k1,k1)
 
           k2 = sum(k(:,1:ng1+m),dim=2)
           k2 = k2 + p
           k2sq = sc_(k2,k2)-mass**2
           sp2 = f(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,flout,flin,ng1)
           if (ng1 >0.or.m>0) sp2 = spb2_(sp2,k2)+mass*sp2
           tmp = vqg(sp2,e1)
 
! print *, m,e1
! pause
           if (m < ng2-1)  then
              if(abs(k1sq) > propcut) then
                  tmp = -(0d0,1d0)/k1sq*tmp
              else
                  tmp = (0d0,0d0)
              endif
           endif
 
           if (ng1>0.or.m>0) then
              if (abs(k2sq) > propcut) then
                  tmp =  (0d0,1d0)/k2sq*tmp
              else
                  tmp = (0d0,0d0)
               endif
           endif
           res = res + tmp
 
        enddo
 
        do m=1,ng1
 
           k1 = sum(k(:,1:m),dim=2)
           e1=g(e(:,1:m),k(:,1:m))
           k1sq = sc_(k1,k1)
 
           k2 = sum(k(:,m+1:ngluon),dim=2)
           k2 = k2 + p
           k2sq = sc_(k2,k2) - mass**2
           ms1 = ng1 - m
           sp2=f(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,flout,flin,ms1)
 
           if (ng2 > 0.or.m < ng1) sp2 = spb2_(sp2,k2)+mass*sp2
           tmp = vgq(e1,sp2)
           if (m > 1) then
              if (abs(k1sq) > propcut) then
                  tmp=-(0d0,1d0)/k1sq*tmp
              else
                  tmp = (0d0,0d0)
              endif
           endif
 
           if (ng2 > 0.or. m < ng1) then
            if (abs(k2sq) > propcut) then
              tmp=(0d0,1d0)/k2sq*tmp
            else
              tmp = (0d0,0d0)
            endif
           endif
 
           res = res + tmp
           enddo
 
         endif
         endif ! endif for flavor consistency condition
 
      end function f
 
 
 
      recursive function bf(e,k,sp,p,mass,flout,flin,ms) result(res)
      implicit none
      complex(8), intent(in) :: e(:,:), k(:,:)
      complex(8), intent(in) :: sp(:), p(:)
      integer, intent(in) ::  ms
      character, intent(in) :: flout*3  ! flavor of on-shell f-line
      character, intent(in) :: flin*3   ! flavor of off-shell f-line
      integer             :: ms1,m,ng1, ng2
      integer :: ngluon
      complex(8)             :: res(size(sp))
      complex(8)             :: tmp(size(sp))
      complex(8)             :: k1(size(p))
      complex(8)             :: k2(size(p))
      complex(8)             :: sp2(size(sp))
      complex(8)             :: sp3(size(sp))
      complex(8)             :: e1(size(e,dim=1))
      complex(8)             :: e2(size(e,dim=1))
      complex(8)  :: k1sq,k2sq,k3sq
      real(8) :: mass
 
 
       if (flout.ne.flin) then
          res = (0d0,0d0)
       else
 
       ngluon = size(e,dim=2)
      ng1 = ms   !#gluons to the left of a f-line
      ng2 = ngluon - ms  !#gluons to the right of the f-line
 
      if (ng2 < 0) write(*,*) 'WRONG DEFINITION OF CURRENT B'
      if (ngluon == 0) then
         res = sp
      else
 
       res = (0d0,0d0)
       do m=0,ng2-1
           k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
           e1=g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
           k1sq=sc_(k1,k1)
 
           k2 = sum(k(:,1:ng1+m),dim=2)
           k2 = -k2 - p
           k2sq = sc_(k2,k2)-mass**2
           sp2 = bf(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,flout,flin,ng1)
           if (ng1>0.or.m>0) sp2 = spi2_(k2,sp2)+mass*sp2
 
            tmp = vbqg(sp2,e1)
 
           if (m < ng2-1) then
            if (abs(k1sq) > propcut) then
           tmp = -(0d0,1d0)/k1sq*tmp
            else
           tmp = (0d0,0d0)
            endif
            endif
           if (ng1>0.or.m>0) then
            if (abs(k2sq) > propcut) then
           tmp =  (0d0,1d0)/k2sq*tmp
            else
            tmp = (0d0,0d0)
             endif
             endif
 
           res = res + tmp
 
 
        enddo
 
 
        do m=1,ng1
 
           k1 = sum(k(:,1:m),dim=2)
           e1=g(e(:,1:m),k(:,1:m))
           k1sq = sc_(k1,k1)
 
           k2 = sum(k(:,m+1:ngluon),dim=2)
           k2 = -k2 - p
           k2sq = sc_(k2,k2) - mass**2
           ms1 = ng1 - m
           sp2=bf(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,flout,flin,ms1)
 
           if (ng2 > 0.or.m < ng1) sp2 = spi2_(k2,sp2)+mass*sp2
 
           tmp = vgbq(e1,sp2)
 
           if (m > 1) then
               if (abs(k1sq) > propcut) then
              tmp=-(0d0,1d0)/k1sq*tmp
              else
              tmp = (0d0,0d0)
              endif
           endif
 
           if (ng2 > 0.or. m < ng1) then
              if (abs(k2sq) > propcut) then
              tmp=(0d0,1d0)/k2sq*tmp
              else
              tmp = (0d0,0d0)
              endif
           endif
 
           res = res + tmp
 
           enddo
 
          endif
          endif ! endif for flavor consisency condition
      end function bf
 
 
 
 
 
 
FUNCTION cur_g_2f(Gluons,Quarks,NumGlu) result(Res)           ! Gluons(:) does not include the OFF-shell gluon,  however NumGlu(0) is the number of all gluons
implicit none
integer :: numglu(0:3),i,counter
type(ptrtoparticle) :: gluons(1:),quarks(1:2)
integer :: rin,rout,n1a,n1b,n2a,n2b,n3a,n3b
integer,target :: tmpextref
complex(8) :: res(1:dv)
complex(8) :: u1(1:ds),ubar2(1:ds)
complex(8),target :: eps1(1:dv)
complex(8) :: eps2(1:dv)
type(ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(3)+1)
complex(8) :: pmom1(1:dv),pmom2(1:dv),pmom4(1:dv)
complex(8),target :: pmom3(1:dv)
complex(8) :: propfac1,propfac2,propfac3,propfac4
integer :: partkey,helkey,currkey,hel_tmp
 
 
 
 
   res = (0d0,0d0)
   if( quarks(1)%PartType .ne. -quarks(2)%PartType ) return
   do n1a=0,numglu(1)
   do n3a=0,numglu(3)
   do n2a=0,numglu(2)
      n1b=numglu(1)-n1a
      n2b=numglu(2)-n2a
      n3b=numglu(3)-n3a
      ! Fer1
      rin=n1a+1
      rout=numglu(1)+n2a
      pmom1(:) = quarks(1)%Mom(:)+ summom(gluons,rin,rout)
      u1(:) = cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/n1b+n2a,n1b,n2a/))
      if(n1b.ge.1 .or. n2a.ge.1) then
         propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(1)%Mass2)
         if( abs(sc_(pmom1,pmom1)-quarks(1)%Mass2).lt.propcut ) cycle
         if( quarks(1)%PartType.lt.0 ) then
            u1(:) = (-spi2_(pmom1,u1)+quarks(1)%Mass*u1(:) )*propfac1
         else
            u1(:) = (+spb2_(u1,pmom1)+quarks(1)%Mass*u1(:) )*propfac1
         endif
      endif
 
      ! Fer2
      rin=numglu(1)+n2a+1
      rout=numglu(1)+numglu(2)+n3a
      pmom2(:) = quarks(2)%Mom(:)+ summom(gluons,rin,rout)
      ubar2(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2b+n3a,n2b,n3a/))
      if(n2b.ge.1 .or. n3a.ge.1) then
         propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
         if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) cycle
         if( quarks(2)%PartType.lt.0 ) then
            ubar2(:) = (-spi2_(pmom2,ubar2)+quarks(2)%Mass*ubar2(:) )*propfac2
         else
            ubar2(:) = (+spb2_(ubar2,pmom2)+quarks(2)%Mass*ubar2(:) )*propfac2
         endif
      endif
 
      if( quarks(1)%PartType.lt.0 ) then
         eps1(:)= -vbqq(dv,ubar2,u1)       ! re-checked
      else
         eps1(:)= +vbqq(dv,u1,ubar2)       ! re-checked
      endif
      pmom3(:) = quarks(1)%Mom(:)+quarks(2)%Mom(:) + summom(gluons,n1a+1,numglu(1)+numglu(2)+n3a)
      counter=1
      rin =1
      rout=n1a
      do i=rin,rout
         call copyparticleptr(gluons(i),tmpgluons(counter))
         counter=counter+1
      enddo
      tmpgluons(counter)%Mom => pmom3(:)
      tmpgluons(counter)%Pol => eps1(:)
      tmpextref = -1
      tmpgluons(counter)%ExtRef => tmpextref
      counter=counter+1
      rin =numglu(1)+numglu(2)+n3a+1
      rout=numglu(1)+numglu(2)+numglu(3)
      do i=rin,rout
         call copyparticleptr(gluons(i),tmpgluons(counter))
         counter=counter+1
      enddo
      eps2(:) = cur_g(tmpgluons(1:counter-1),1+n1a+n3b+1)
 
 
      if(n1a.ge.1 .or. n3b.ge.1) then
         propfac3 = (0d0,-1d0)/sc_(pmom3,pmom3)
         if( abs(sc_(pmom3,pmom3)).lt.propcut ) cycle
         eps2(:) = eps2(:)*propfac3
      endif
 
      res(:) = res(:) + eps2(:)
   enddo
   enddo
   enddo
 
 
 
return
END FUNCTION
 
 
 
 
 
FUNCTION cur_g_2fv(Gluons,Quarks,Boson,NumGlu) result(Res)           ! Gluons(:) does not include the OFF-shell gluon,  however NumGlu(0) is the number of all gluons
implicit none
integer :: numglu(0:3),i,counter
type(ptrtoparticle) :: gluons(1:),quarks(1:2),boson
integer :: rin,rout,n1a,n1b,n2a,n2b,n3a,n3b
integer,target :: tmpextref
complex(8) :: res(1:dv)
complex(8) :: u1(1:ds),ubar2(1:ds)
complex(8),target :: eps1(1:dv)
complex(8) :: eps2(1:dv)
type(ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(3)+1)
complex(8) :: pmom1(1:dv),pmom2(1:dv),pmom4(1:dv)
complex(8),target :: pmom3(1:dv)
complex(8) :: propfac1,propfac2,propfac3,propfac4
integer :: partkey,helkey,currkey,hel_tmp
 
 
 
 
!DEC$ IF (_DebugCheckMyImpl1==1)
   if(quarks(1)%PartType*quarks(2)%PartType.ge.0) print *,"Error in cur_g_2f: wrong PartTypes"
!DEC$ ENDIF
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-1-numglu(1)-numglu(2)-numglu(3).ne.0 ) print *, "wrong number of gluons in cur_g_2f"
!DEC$ ENDIF
 
   res = (0d0,0d0)
   if( quarks(1)%PartType .ne. -quarks(2)%PartType ) return
   do n1a=0,numglu(1)
   do n3a=0,numglu(3)
   do n2a=0,numglu(2)
      n1b=numglu(1)-n1a
      n2b=numglu(2)-n2a
      n3b=numglu(3)-n3a
 
      ! Fer1 and V coupling
      rin=n1a+1
      rout=numglu(1)+n2a
      pmom1(:) = quarks(1)%Mom(:) + summom(gluons,rin,rout) + boson%Mom(:)
      u1(:) = cur_f_2fv(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,boson,(/n1b+n2a,n1b,n2a/))
      !if(n1b.ge.1 .or. n2a.ge.1) then
         propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(1)%Mass2)
         if( abs(sc_(pmom1,pmom1)-quarks(1)%Mass2).lt.propcut ) then
            propfac1=(0d0,0d0)
         endif
         if( quarks(1)%PartType.lt.0 ) then
            u1(:) = (-spi2_(pmom1,u1)+quarks(1)%Mass*u1(:) )*propfac1
         else
            u1(:) = (+spb2_(u1,pmom1)+quarks(1)%Mass*u1(:) )*propfac1
         endif
      !endif
 
      ! Fer2
      rin=numglu(1)+n2a+1
      rout=numglu(1)+numglu(2)+n3a
      pmom2(:) = quarks(2)%Mom(:)+ summom(gluons,rin,rout)
      ubar2(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2b+n3a,n2b,n3a/))
      if(n2b.ge.1 .or. n3a.ge.1) then
         propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
         if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) then
            propfac2=(0d0,0d0)
         endif
         if( quarks(2)%PartType.lt.0 ) then
            ubar2(:) = (-spi2_(pmom2,ubar2)+quarks(2)%Mass*ubar2(:) )*propfac2
         else
            ubar2(:) = (+spb2_(ubar2,pmom2)+quarks(2)%Mass*ubar2(:) )*propfac2
         endif
      endif
 
 
      if( quarks(1)%PartType.lt.0 ) then
         eps1(:)= -vbqq(dv,ubar2,u1)
      else
         eps1(:)= +vbqq(dv,u1,ubar2)
      endif
 
 
 
 
 
      ! Fer1
      rin=n1a+1
      rout=numglu(1)+n2a
      pmom1(:) = quarks(1)%Mom(:) + summom(gluons,rin,rout)
      u1(:) = cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/n1b+n2a,n1b,n2a/))
      if(n1b.ge.1 .or. n2a.ge.1) then
         propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(1)%Mass2)
         if( abs(sc_(pmom1,pmom1)-quarks(1)%Mass2).lt.propcut ) cycle
         if( quarks(1)%PartType.lt.0 ) then
            u1(:) = (-spi2_(pmom1,u1)+quarks(1)%Mass*u1(:) )*propfac1
         else
            u1(:) = (+spb2_(u1,pmom1)+quarks(1)%Mass*u1(:) )*propfac1
         endif
      endif
 
      ! Fer2 and V coupling
      rin=numglu(1)+n2a+1
      rout=numglu(1)+numglu(2)+n3a
      pmom2(:) = quarks(2)%Mom(:)+ summom(gluons,rin,rout) + boson%Mom(:)
      ubar2(:) = cur_f_2fv(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,boson,(/n2b+n3a,n2b,n3a/))
      !if(n2b.ge.1 .or. n3a.ge.1) then
         propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
         if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) cycle
         if( quarks(2)%PartType.lt.0 ) then
            ubar2(:) = (-spi2_(pmom2,ubar2)+quarks(2)%Mass*ubar2(:) )*propfac2
         else
            ubar2(:) = (+spb2_(ubar2,pmom2)+quarks(2)%Mass*ubar2(:) )*propfac2
         endif
      !endif
 
      if( quarks(1)%PartType.lt.0 ) then
         eps1(:)= eps1(:) - vbqq(dv,ubar2,u1)
      else
         eps1(:)= eps1(:) + vbqq(dv,u1,ubar2)
      endif
!
 
      pmom3(:) = quarks(1)%Mom(:)+quarks(2)%Mom(:) + summom(gluons,n1a+1,numglu(1)+numglu(2)+n3a) + boson%Mom(:)
      counter=1
      rin =1
      rout=n1a
      do i=rin,rout
         call copyparticleptr(gluons(i),tmpgluons(counter))
         counter=counter+1
      enddo
      tmpgluons(counter)%Mom => pmom3(:)
      tmpgluons(counter)%Pol => eps1(:)
      tmpextref = -1
      tmpgluons(counter)%ExtRef => tmpextref
      counter=counter+1
      rin =numglu(1)+numglu(2)+n3a+1
      rout=numglu(1)+numglu(2)+numglu(3)
      do i=rin,rout
         call copyparticleptr(gluons(i),tmpgluons(counter))
         counter=counter+1
      enddo
      eps2(:) = cur_g(tmpgluons(1:counter-1),1+n1a+n3b+1)
 
 
      if(n1a.ge.1 .or. n3b.ge.1) then
         propfac3 = (0d0,-1d0)/sc_(pmom3,pmom3)
         if( abs(sc_(pmom3,pmom3)).lt.propcut ) cycle
         eps2(:) = eps2(:)*propfac3
      endif
 
      res(:) = res(:) + eps2(:)
 
   enddo
   enddo
   enddo
 
 
 
return
END FUNCTION
 
 
 
 
 
 
 
 
FUNCTION cur_f_2f(Gluons,Quarks,Quark1PartType,NumGlu) result(Res)           ! Quarks(:) does not include the OFF-shell quark
implicit none
complex(8) :: res(1:ds)
integer :: numglu(0:2),i,rin,rout,quark1parttype
type(ptrtoparticle) :: gluons(1:),quarks(2:2)      ! off-shell quark is not included in Quarks(:)
complex(8) :: glumom(1:dv,numglu(0)), quarkmom(1:dv)
complex(8) :: glupol(1:dv,numglu(0)), quarkpol(1:ds)
character :: ferfla1*3,ferfla2*3
integer :: partkey,helkey,currkey,hel_tmp
 
 
 
 
 
 
!DEC$ IF (_DebugGeneralChecks==1)
   if( quarks(2)%PartType .eq.0 .or. .not.isaquark(quarks(2)%PartType)) then
      print *, "Error in cur_f_2f"
      stop
   endif
!DEC$ ENDIF
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-numglu(1)-numglu(2).ne.0 ) print *, "wrong number of gluons in cur_f_2f"
    if( quarks(2)%PartType.ne.-quark1parttype ) print *, "unequal flavors in cur_f_2f"
!DEC$ ENDIF
 
 
   do i=1,numglu(0)
    glumom(1:dv,i) = gluons(i)%Mom(1:dv)
    glupol(1:dv,i) = gluons(i)%Pol(1:dv)
   enddo
   quarkmom(1:dv) = quarks(2)%Mom(1:dv)
   quarkpol(1:ds) = quarks(2)%Pol(1:ds)
 
   if( abs(quark1parttype).eq.5 ) then
    ferfla1="top"
   elseif( abs(quark1parttype).eq.6 ) then
    ferfla1="bot"
   elseif( abs(quark1parttype).eq.3 ) then
    ferfla1="chm"
   else
    ferfla1="str"
   endif
 
   if( abs(quarks(2)%PartType).eq.5 ) then
    ferfla2="top"
   elseif( abs(quarks(2)%PartType).eq.6 ) then
    ferfla2="bot"
   elseif( abs(quarks(2)%PartType).eq.3 ) then
    ferfla2="chm"
   else
    ferfla2="str"
   endif
 
   rin =1
   rout=numglu(0)
   if( quarks(2)%PartType .gt.0 ) then      !    X----->----
      res(:) = f(glupol(1:dv,rin:rout),glumom(1:dv,rin:rout),quarkpol(1:ds),quarkmom(1:dv),quarks(2)%Mass,ferfla2,ferfla1,numglu(1))
   else                                     !    X-----<----
      res(:) = bf(glupol(1:dv,rin:rout),glumom(1:dv,rin:rout),quarkpol(1:ds),quarkmom(1:dv),quarks(2)%Mass,ferfla2,ferfla1,numglu(1))
   endif
 
 
return
END FUNCTION
 
 
 
 
 
 
FUNCTION cur_g_4f(Gluons,Quarks,NumGlu) result(res)           ! Gluons(:) does not include the OFF-shell gluon,  however NumGlu is the number of all gluons
implicit none
integer,intent(in) :: numglu(0:5)
type(ptrtoparticle) :: gluons(1:),quarks(1:)
integer :: na,nb,nc,nd,ne,nf,ng,nh,ni,nj,nk
integer :: rin,rout
integer :: tag_f,counter,i
complex(8) :: res(dv)
type(ptrtoparticle) :: tmpgluons(1:2+numglu(1)+numglu(3)+numglu(5))
complex(8),target :: tmpmom1(1:dv),tmpmom2(1:dv)
integer,target :: tmpextref1,tmpextref2
complex(8),target :: eps1(1:dv)
complex(8),target :: eps2(1:dv)
complex(8) :: eps3(1:dv)
complex(8) :: u1(1:ds)
complex(8) :: ubar2(1:ds)
complex(8) :: propfac1,propfac2,propfac3,propfac4
complex(8) :: pmom1(1:dv)
complex(8) :: pmom2(1:dv)
complex(8) :: pmom3(1:dv)
complex(8) :: pmom4(1:dv)
 
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-1-numglu(1)-numglu(2)-numglu(3)-numglu(4)-numglu(5).ne.0 ) print *, "wrong number of gluons in cur_g_4f"
!DEC$ ENDIF
 
      res = (0d0,0d0)
      if( (quarks(1)%PartType.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType) &
     .OR. (quarks(1)%PartType.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType) ) then
!        type (1)
         do na=0,numglu(1)
         do nc=0,numglu(2)
         do ne=0,numglu(5)
            nb=numglu(1)-na
            nd=numglu(2)-nc
            nf=numglu(5)-ne
            rin = numglu(1)+nc+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
            u1 = cur_f_4f(gluons(rin:rout),quarks(2:4),quarks(1)%PartType,(/nd+numglu(3)+numglu(4)+ne,nd,numglu(3),numglu(4),ne/),0,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
            propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(1)%Mass2)
            if( abs(sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.propcut ) cycle
            if( quarks(1)%PartType.lt.0 ) then
               u1 = ( spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
            else
               u1 = (-spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
            endif
 
            rin = na+1
            rout= numglu(1)+nc
            ubar2 = cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
            pmom3 = summom(gluons,rin,rout) + quarks(1)%Mom
            if( nb.ge.1 .or. nc.ge.1 ) then
               propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(1)%Mass2)
               if( abs(sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.propcut ) cycle
               if( quarks(1)%PartType.lt.0 ) then
                  ubar2 = (-spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
               else
                  ubar2 = (+spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
               endif
            endif
 
            if( quarks(1)%PartType.lt.0 ) then
                eps1 = -vbqq(dv,u1,ubar2)       ! re-checked
            else
                eps1 = +vbqq(dv,ubar2,u1)       ! re-checked
            endif
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpmom1(:) = pmom2(:)+pmom3(:)
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
            if( na.ge.1 .or. nf.ge.1 ) then
               propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
               if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
               eps2 = eps2*propfac1
            endif
 
            res = res + eps2
         enddo
         enddo
         enddo
 
 
!        type (2)
         do na=0,numglu(1)
         do nc=0,numglu(4)
         do ne=0,numglu(5)        ! can be replaced by above ne-loop
            nb=numglu(1)-na
            nd=numglu(4)-nc
            nf=numglu(5)-ne
 
            rin = na+1
            rout= numglu(1)+numglu(2)+numglu(3)+nc
            u1 = cur_f_4f(gluons(rin:rout),quarks(1:3),quarks(4)%PartType,(/nb+numglu(2)+numglu(3)+nc,nb,numglu(2),numglu(3),nc/),0,0)
            pmom2 =  summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom
            propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
            if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
            if( quarks(4)%PartType.lt.0 ) then
              u1 = (+spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
            else
              u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
            endif
            rin = numglu(1)+numglu(2)+numglu(3)+nc+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
            ubar2 = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/nd+ne,nd,ne/))
            pmom3  = summom(gluons,rin,rout) + quarks(4)%Mom
            if( nd.ge.1 .or. ne.ge.1 ) then
               propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(4)%Mass2)
               if( abs(sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.propcut ) cycle
               if( quarks(4)%PartType.lt.0 ) then
                  ubar2 = (-spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
               else
                  ubar2 = (+spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
               endif
            endif
 
            if( quarks(4)%PartType.lt.0 ) then
              eps1 = +vbqq(dv,u1,ubar2)       ! re-checked
            else
              eps1 = -vbqq(dv,ubar2,u1)       ! re-checked
            endif
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpmom1(:) = pmom2(:)+pmom3(:)
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
            if( na.ge.1 .or. nf.ge.1 ) then
               propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
               if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
               eps2 = eps2*propfac1
            endif
 
            res = res + eps2
         enddo
         enddo
         enddo
      endif
 
 
 
      if( quarks(1)%PartType.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType) then
!        type(3)
         do na=0,numglu(1)
         do nc=0,numglu(2)
         do ne=0,numglu(3)
         do nf=0,numglu(3)-ne
         do nh=0,numglu(4)   ! this loop can be placed after Eps1 has been calculated
         do nj=0,numglu(5)
            nb=numglu(1)-na
            nd=numglu(2)-nc
            ng=numglu(3)-ne-nf
            ni=numglu(4)-nh
            nk=numglu(5)-nj
 
            rin = na+1
            rout= numglu(1)+nc
            ubar2 = cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
            pmom1 = summom(gluons,rin,rout) + quarks(1)%Mom
            if( nb.ge.1 .or. nc.ge.1 ) then
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1) - quarks(1)%Mass2)
               if( abs(sc_(pmom1,pmom1) - quarks(1)%Mass2).lt.propcut ) cycle
               if( quarks(1)%PartType.lt.0 ) then
                  ubar2 = (-spi2_(pmom1,ubar2) + quarks(1)%Mass*ubar2 )*propfac1
               else
                  ubar2 = (+spb2_(ubar2,pmom1) + quarks(1)%Mass*ubar2 )*propfac1
               endif
            endif
            rin = numglu(1)+nc+1
            rout= numglu(1)+numglu(2)+ne
            u1 = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/nd+ne,nd,ne/))
            pmom2 = summom(gluons,rin,rout) + quarks(2)%Mom
            if( nd.ge.1 .or. ne.ge.1 ) then
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
               if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  u1 = (-spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
               else
                  u1 = (+spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
               endif
            endif
 
            if( quarks(2)%PartType.lt.0 ) then
              eps1 = +vbqq(dv,ubar2,u1)       ! re-checked
            else
              eps1 = -vbqq(dv,u1,ubar2)       ! re-checked
            endif
            tmpmom1 = pmom1 + pmom2
            propfac3 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
            if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
            eps1 = eps1*propfac3
 
 
            rin = numglu(1)+numglu(2)+ne+nf+1
            rout= numglu(1)+numglu(2)+numglu(3)+nh
            u1 = cur_f_2f(gluons(rin:rout),quarks(3:3),-quarks(3)%PartType,(/ng+nh,ng,nh/))
            pmom3 = summom(gluons,rin,rout) + quarks(3)%Mom
            if( ng.ge.1 .or. nh.ge.1 ) then
               propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(3)%Mass2)
               if( abs(sc_(pmom3,pmom3) - quarks(3)%Mass2).lt.propcut ) cycle
               if( quarks(3)%PartType.lt.0 ) then
                  u1 = (-spi2_(pmom3,u1) + quarks(3)%Mass*u1 )*propfac3
               else
                  u1 = (+spb2_(u1,pmom3) + quarks(3)%Mass*u1 )*propfac3
               endif
            endif
            rin = numglu(1)+numglu(2)+numglu(3)+nh+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj
            ubar2 = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/ni+nj,ni,nj/))
            pmom4 = summom(gluons,rin,rout) + quarks(4)%Mom
            if( ni.ge.1 .or. nj.ge.1 ) then
               propfac4 = (0d0,1d0)/(sc_(pmom4,pmom4) - quarks(4)%Mass2)
               if( abs(sc_(pmom4,pmom4) - quarks(4)%Mass2).lt.propcut ) cycle
               if( quarks(4)%PartType.lt.0 ) then
                  ubar2 = (-spi2_(pmom4,ubar2) + quarks(4)%Mass*ubar2 )*propfac4
               else
                  ubar2 = (+spb2_(ubar2,pmom4) + quarks(4)%Mass*ubar2 )*propfac4
               endif
            endif
 
            if( quarks(4)%PartType.lt.0 ) then
               eps2 = +vbqq(dv,u1,ubar2)       ! re-checked
            else
               eps2 = -vbqq(dv,ubar2,u1)       ! re-checked
            endif
            tmpmom2 = pmom3 + pmom4
            propfac1 = (0d0,-1d0)/sc_(tmpmom2,tmpmom2)
            if( abs(sc_(tmpmom2,tmpmom2)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+ne+1
            rout=numglu(1)+numglu(2)+ne+nf
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpgluons(counter)%Mom => tmpmom2(:)
            tmpgluons(counter)%Pol => eps2(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps3(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+nk+2)
            res = res + eps3
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
      endif
 
 
return
END FUNCTION
 
 
 
 
FUNCTION cur_f_6f(Gluons,Quarks,Quark1PartType,NumGlu,tag_f,tag_Z) result(res)           ! Quarks(:) does not include the OFF-shell quark
implicit none
integer :: numglu(0:6),quark1parttype,tag_f,tag_z
type(ptrtoparticle) :: gluons(1:),quarks(2:6)
integer,target :: tmpparttype,tmpextref
complex(8) :: res(1:ds),tmp(1:ds)
! complex(8) :: Res1(1:Ds),Res2(1:Ds),Res3(1:Ds),Res4(1:Ds)
complex(8) :: u1(1:ds),ubar1(1:ds)
complex(8),target :: ubar0(1:ds)
complex(8) :: eps1(1:dv)
complex(8) :: eps2(1:dv)
type(ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(6)),tmpquark(1:1)
complex(8) :: propfac1,propfac2
complex(8),target :: pmom1(1:dv)
complex(8),target :: pmom2(1:dv)
integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b,n5a,n5b,n6a,n6b
integer :: rin,rout,i,counter
 
 
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-numglu(1)-numglu(2)-numglu(3)-numglu(4)-numglu(5)-numglu(6).ne.0 ) print *, "wrong number of gluons in cur_f_6f"
!DEC$ ENDIF
 
    res = (0d0,0d0)
!     Res1=(0d0,0d0); Res2=(0d0,0d0); Res3=(0d0,0d0); Res4=(0d0,0d0);
 
 
!   (A)
    if( quark1parttype.eq.-quarks(2)%PartType .AND. (quarks(3)%PartType.eq.-quarks(4)%PartType .or. quarks(3)%PartType.eq.-quarks(6)%PartType) &
        .AND. (quarks(2)%ExtRef.ne.-1 .or. tag_f.ne.1) &
      ) then
 
      do n2a=0,numglu(2)
      do n6a=0,numglu(6)
         n2b = numglu(2)-n2a
         n6b = numglu(6)-n6a
 
         rin =numglu(1)+n2a+1
         rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
         eps2 = cur_g_4f(gluons(rin:rout),quarks(3:6),(/1+n2b+numglu(3)+numglu(4)+numglu(5)+n6a,n2b,numglu(3),numglu(4),numglu(5),n6a/))
         pmom1(:) = summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
         propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
         if( abs(sc_(pmom1,pmom1)).lt.propcut) cycle
         eps2 = eps2*propfac1
         do n1a=0,numglu(1)
            n1b = numglu(1)-n1a
            ! Fer2
            rin =n1a+1
            rout=numglu(1)+n2a
            ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n1b+n2a,n1b,n2a/) )
            if(n1b.ge.1 .or. n2a.ge.1) then
               pmom2(:) = quarks(2)%Mom + summom(gluons,rin,rout)
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
               if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
               else
                  ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
               endif
            endif
            if( quarks(2)%PartType.lt.0 ) then
               ubar0(:) = vbqg(ubar1,eps2)       ! re-checked
            else
               ubar0(:) = vqg(ubar1,eps2)       ! re-checked
            endif
 
            rin = n1a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            pmom1 = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom  ! can be simplified with PMom1(:)
            if(n1a.ge.1 .or. n6b.ge.1) then
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
               if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
               else
                  ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
               endif
            endif
 
            tmpquark(1)%Mom  => pmom1(:)
            tmpquark(1)%Pol  => ubar0(:)
            tmpquark(1)%Mass => quarks(2)%Mass
            tmpquark(1)%Mass2=> quarks(2)%Mass2
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpquark(1)%PartType => quarks(2)%PartType
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n6b/) )
 
            res(:) = res(:) + tmp(:)
!             Res1(:) = Res1(:) + tmp(:)
!             print *, "1",tmp(:)
         enddo
      enddo
      enddo
    endif
 
 
 
!   (B)
    if( quark1parttype.eq.-quarks(6)%PartType .AND. (quarks(2)%PartType.eq.-quarks(5)%PartType .or. quarks(2)%PartType.eq.-quarks(3)%PartType) &
        .AND. (quarks(6)%ExtRef.ne.-1 .or. tag_f.ne.1) &
      ) then
      do n1a=0,numglu(1)
      do n5b=0,numglu(5)
         n1b = numglu(1)-n1a
         n5a = numglu(5)-n5b
 
         rin =n1a+1
         rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+n5a
 
! preventing color issues with Z on the loop, see RR notes 05-22-2013
!note: only tested for zero gluons
         if ( tag_z .eq. 1 .and. (n1b+numglu(2)+numglu(3)+numglu(4)+n5a .eq. numglu(0)) ) then
!             print  * , "cycle for tag_Z=1 in cur_f_6f"
            cycle
         endif
 
         eps2 = cur_g_4f(gluons(rin:rout),quarks(2:5),(/1+n1b+numglu(2)+numglu(3)+numglu(4)+n5a,n1b,numglu(2),numglu(3),numglu(4),n5a/))
         pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom
         propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
         if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
         eps2 = eps2*propfac1
 
         do n6a=0,numglu(6)
            n6b = numglu(6)-n6a
            ! Fer6
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+n5a+1
            rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(6:6),-quarks(6)%PartType,(/n5b+n6a,n5b,n6a/) )
            if(n5b.ge.1 .or. n6a.ge.1) then
               pmom2(:) = summom(gluons,rin,rout) + quarks(6)%Mom
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(6)%Mass2)
               if( abs(sc_(pmom2,pmom2)-quarks(6)%Mass2).lt.propcut ) cycle
               if( quarks(6)%PartType.lt.0 ) then
                  ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(6)%Mass*ubar1(:))*propfac2
               else
                  ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(6)%Mass*ubar1(:))*propfac2
               endif
            endif
            if( quarks(6)%PartType.lt.0 ) then
               ubar0(:) = vgbq(eps2,ubar1)       ! re-checked
            else
               ubar0(:) = vgq(eps2,ubar1)       ! re-checked
            endif
 
            rin = n1a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            pmom1 = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
            if(n1a.ge.1 .or. n6b.ge.1) then
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(6)%Mass2)
               if( abs(sc_(pmom1,pmom1)-quarks(6)%Mass2).lt.propcut ) cycle
               if( quarks(6)%PartType.lt.0 ) then
                  ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(6)%Mass*ubar0(:))*propfac1
               else
                  ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(6)%Mass*ubar0(:))*propfac1
               endif
            endif
            tmpquark(1)%Mom => pmom1(:)
            tmpquark(1)%Pol => ubar0(:)
            tmpquark(1)%Mass => quarks(6)%Mass
            tmpquark(1)%Mass2=> quarks(6)%Mass2
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpquark(1)%PartType => quarks(6)%PartType
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
 
            res(:) = res(:) + tmp(:)
 
!             Res2(:) = Res2(:) + tmp(:)
 
!             print *, "2",tmp(:)
         enddo
      enddo
      enddo
    endif
 
 
 
 
!   (C)
!     if( Quarks(5)%PartType.eq.-Quarks(6)%PartType .AND. (Quark1PartType.eq.-Quarks(2)%PartType .or. Quark1PartType.eq.-Quarks(4)%PartType) &
!         .AND. .not.(Quarks(6)%ExtRef.eq.-1 .and. tag_f.eq.1) &
!       ) then
!     if( Quark1PartType.eq.-Quarks(2)%PartType .and. Quarks(5)%PartType.eq.-Quarks(6)%PartType .and..not.(Quarks(2)%ExtRef.eq.-1 .and. tag_f.eq.1) &
!    .OR. Quark1PartType.eq.-Quarks(4)%PartType .and. Quarks(5)%PartType.eq.-Quarks(6)%PartType .and..not.(Quarks(4)%ExtRef.eq.-1 .and. tag_f.eq.1) &
!       ) then
 
    if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
        ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
    .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
      ) then
 
      do n1a=0,numglu(1)
      do n4a=0,numglu(4)
      do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n4b = numglu(4)-n4a
         n6b = numglu(6)-n6a
 
            rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            eps2 = cur_g_2f(gluons(rin:rout),quarks(5:6),(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
            pmom1(:) = summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            rin = n1a+1
            rout= numglu(1)+numglu(2)+numglu(3)+n4a
            u1 = cur_f_4f(gluons(rin:rout),quarks(2:4),quark1parttype,(/n1b+numglu(2)+numglu(3)+n4a,n1b,numglu(2),numglu(3),n4a/),0,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
 
            if( quark1parttype.eq.-quarks(2)%PartType) then
                if( quarks(2)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
                  ubar0 = vbqg(u1,eps2)       ! re-checked
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom   ! this PMom2 will be re-used below
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(2)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(2)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
                  ubar0 = vqg(u1,eps2)       ! re-checked
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(2)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(2)%Mass
                tmpquark(1)%Mass2=> quarks(2)%Mass2
            elseif( quark1parttype.eq.-quarks(4)%PartType) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vbqg(u1,eps2)       ! re-checked
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vqg(u1,eps2)       ! re-checked
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(4)%Mass
                tmpquark(1)%Mass2=> quarks(4)%Mass2
            endif
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
 
            res(:) = res(:) + tmp(:)
!             Res3(:) = Res3(:) + tmp(:)
!             print *, "3",tmp(:)
      enddo
      enddo
      enddo
    endif
 
 
!   (D)
    if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
        (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
   .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1))  )) then
!     if( Quark1PartType.eq.-Quarks(4)%PartType .and. Quarks(2)%PartType.eq.-Quarks(3)%PartType .and..not.(Quarks(4)%ExtRef.eq.-1 .and. tag_f.eq.1) &
!     .OR.Quark1PartType.eq.-Quarks(6)%PartType .and. Quarks(2)%PartType.eq.-Quarks(3)%PartType .and..not.(Quarks(6)%ExtRef.eq.-1 .and. tag_f.eq.1) &
!       ) then
 
      do n1a=0,numglu(1)
      do n3a=0,numglu(3)
      do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n3b = numglu(3)-n3a
         n6b = numglu(6)-n6a
 
            counter=1
            rin = n1a+1
            rout= numglu(1)+numglu(2)+n3a
 
 
            eps2 = cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
            pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            rin = numglu(1)+numglu(2)+n3a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            u1 = cur_f_4f(gluons(rin:rout),quarks(4:6),quark1parttype,(/n3b+numglu(4)+numglu(5)+n6a,n3b,numglu(4),numglu(5),n6a/),tag_f,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
 
            if( quark1parttype.eq.-quarks(4)%PartType ) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom   !this PMom2 will be re-used below  ! CAN be written as PMom2=PMom2+PMom1
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(4)%Mass
                tmpquark(1)%Mass2=> quarks(4)%Mass2
            elseif(quark1parttype.eq.-quarks(6)%PartType) then
                if( quarks(6)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(6)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(6)%Mass
                tmpquark(1)%Mass2=> quarks(6)%Mass2
            endif
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/))
 
            res(:)  = res(:) + tmp(:)
!             Res4(:) = Res4(:) + Tmp(:)
!             print *, "4",tmp(:)
      enddo
      enddo
      enddo
    endif
!    print *, 'res', res
!    pause
!     Res(:) = Res1(:)+Res2(:)+Res3(:)+Res4(:)
!     print *, "Res1:",Res1(1:4)
!     print *, "Res2:",Res2(1:4)
!     print *, "Res3:",Res3(1:4)
!     print *, "Res4:",Res4(1:4)
return
END FUNCTION
 
 
 
 
 
 
 
 
FUNCTION cur_g_4fv(Gluons,Quarks,Boson,BosonVertex,NumGlu) result(res)           ! Gluons(:) does not include the OFF-shell gluon,  however NumGlu is the number of all gluons
implicit none
integer,intent(in) :: numglu(0:5),bosonvertex
type(ptrtoparticle) :: gluons(1:),quarks(1:),boson
integer :: na,nb,nc,nd,ne,nf,ng,nh,ni,nj,nk,bosonvertex_mod
integer :: rin,rout
integer :: tag_f,counter,i
complex(8) :: res(dv)
type(ptrtoparticle) :: tmpgluons(1:2+numglu(1)+numglu(3)+numglu(5))
complex(8),target :: tmpmom1(1:dv),tmpmom2(1:dv)
integer,target :: tmpextref1,tmpextref2
complex(8),target :: eps1(1:dv)
complex(8),target :: eps2(1:dv)
complex(8) :: eps3(1:dv)
complex(8) :: u1(1:ds)
complex(8) :: ubar2(1:ds)
complex(8) :: propfac1,propfac2,propfac3,propfac4
complex(8) :: pmom1(1:dv)
complex(8) :: pmom2(1:dv)
complex(8) :: pmom3(1:dv)
complex(8) :: pmom4(1:dv)
 
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-1-numglu(1)-numglu(2)-numglu(3)-numglu(4)-numglu(5).ne.0 ) print *, "wrong number of gluons in cur_g_4f"
!DEC$ ENDIF
 
    res = (0d0,0d0)
    if (quarks(1)%PartType.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType) then
       ! for this flavor structure, the BosonVertex implies the following:
       ! BV = 1 : Boson on quark 1 only
       ! BV = 2 : Boson on both quarks 2 and 3
       ! BV = 3 : Boson on quark 4 only
       ! See RR notes
 
       if( bosonvertex.eq.1)  then
          !        type (1)
         do na=0,numglu(1)
         do nc=0,numglu(2)
         do ne=0,numglu(5)
            nb=numglu(1)-na
            nd=numglu(2)-nc
            nf=numglu(5)-ne
            rin = numglu(1)+nc+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
            u1 = cur_f_4f(gluons(rin:rout),quarks(2:4),quarks(1)%PartType,(/nd+numglu(3)+numglu(4)+ne,nd,numglu(3),numglu(4),ne/),0,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
            propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(1)%Mass2)
            if( abs(sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.propcut ) cycle
            if( quarks(1)%PartType.lt.0 ) then
               u1 = ( spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
            else
               u1 = (-spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
            endif
 
            rin = na+1
            rout= numglu(1)+nc
            ubar2 = cur_f_2fv(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,boson,(/nb+nc,nb,nc/))
            pmom3 = summom(gluons,rin,rout) + quarks(1)%Mom + boson%Mom
            propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(1)%Mass2)
            if( abs(sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.propcut ) cycle
            if( quarks(1)%PartType.lt.0 ) then
               ubar2 = (-spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
            else
               ubar2 = (+spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
            endif
 
            if( quarks(1)%PartType.lt.0 ) then
                eps1 = -vbqq(dv,u1,ubar2)       ! re-checked
            else
                eps1 = +vbqq(dv,ubar2,u1)       ! re-checked
            endif
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpmom1(:) = pmom2(:)+pmom3(:)
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
            if( na.ge.1 .or. nf.ge.1 ) then
               propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
               if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
               eps2 = eps2*propfac1
            endif
 
            res = res + eps2
         enddo
         enddo
         enddo
      endif! Bosonvertex
 
      if( bosonvertex.eq.1 .or. bosonvertex.eq.2 .or. bosonvertex .eq. 3)  then
!        type (2)
         do na=0,numglu(1)
         do nc=0,numglu(4)
         do ne=0,numglu(5)
            nb=numglu(1)-na
            nd=numglu(4)-nc
            nf=numglu(5)-ne
 
            rin = na+1
            rout= numglu(1)+numglu(2)+numglu(3)+nc
            bosonvertex_mod = bosonvertex  + 1
!            if (BosonVertex .eq. 1) BosonVertex_mod=BosonVertex
            u1 = cur_f_4fv(gluons(rin:rout),quarks(1:3),quarks(4)%PartType,boson,bosonvertex_mod,(/nb+numglu(2)+numglu(3)+nc,nb,numglu(2),numglu(3),nc/),0,0)
            pmom2 =  summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
            propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
            if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
            if( quarks(4)%PartType.lt.0 ) then
              u1 = (+spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
            else
              u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
            endif
 
            rin = numglu(1)+numglu(2)+numglu(3)+nc+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
            ubar2 = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/nd+ne,nd,ne/))
            pmom3  = summom(gluons,rin,rout) + quarks(4)%Mom
            if( nd.ge.1 .or. ne.ge.1 ) then
               propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(4)%Mass2)
               if( abs(sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.propcut ) cycle
               if( quarks(4)%PartType.lt.0 ) then
                  ubar2 = (-spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
               else
                  ubar2 = (+spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
               endif
            endif
 
            if( quarks(4)%PartType.lt.0 ) then
              eps1 = +vbqq(dv,u1,ubar2)
            else
              eps1 = -vbqq(dv,ubar2,u1)
            endif
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpmom1(:) = pmom2(:)+pmom3(:)
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
            if( na.ge.1 .or. nf.ge.1 ) then
               propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
               if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
               eps2 = eps2*propfac1
            endif
 
            res = res + eps2
 
         enddo
         enddo
         enddo
 
      endif ! BosonVertex
 
 
      if (bosonvertex .eq. 1 .or. bosonvertex .eq. 2 .or. bosonvertex .eq. 3 ) then
 
         do na=0,numglu(1)
         do nc=0,numglu(2)
         do ne=0,numglu(5)
            nb=numglu(1)-na
            nd=numglu(2)-nc
            nf=numglu(5)-ne
            rin = numglu(1)+nc+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
            bosonvertex_mod=bosonvertex
            u1 = cur_f_4fv(gluons(rin:rout),quarks(2:4),quarks(1)%PartType,boson,bosonvertex_mod,(/nd+numglu(3)+numglu(4)+ne,nd,numglu(3),numglu(4),ne/),0,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + boson%Mom
            propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(1)%Mass2)
            if( abs(sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.propcut ) cycle
            if( quarks(1)%PartType.lt.0 ) then
               u1 = ( spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
            else
               u1 = (-spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
            endif
 
            rin = na+1
            rout= numglu(1)+nc
            ubar2 = cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
            pmom3 = summom(gluons,rin,rout) + quarks(1)%Mom
            propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(1)%Mass2)
            if ( nb .ge. 1 .or. nc .ge. 1) then
               if( abs(sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.propcut ) cycle
               if( quarks(1)%PartType.lt.0 ) then
                  ubar2 = (-spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
               else
                  ubar2 = (+spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
               endif
            endif
 
            if( quarks(1)%PartType.lt.0 ) then
                eps1 = -vbqq(dv,u1,ubar2)       ! re-checked
            else
                eps1 = +vbqq(dv,ubar2,u1)       ! re-checked
            endif
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpmom1(:) = pmom2(:)+pmom3(:)
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
            if( na.ge.1 .or. nf.ge.1 ) then
               propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
               if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
               eps2 = eps2*propfac1
            endif
 
            res = res + eps2
         enddo
         enddo
         enddo
      endif! Bosonvertex
 
 
      if( bosonvertex.eq.3)  then
!        type (2)
         do na=0,numglu(1)
         do nc=0,numglu(4)
         do ne=0,numglu(5)
            nb=numglu(1)-na
            nd=numglu(4)-nc
            nf=numglu(5)-ne
 
            rin = na+1
            rout= numglu(1)+numglu(2)+numglu(3)+nc
            u1 = cur_f_4f(gluons(rin:rout),quarks(1:3),quarks(4)%PartType,(/nb+numglu(2)+numglu(3)+nc,nb,numglu(2),numglu(3),nc/),0,0)
            pmom2 =  summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom
            propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
            if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
            if( quarks(4)%PartType.lt.0 ) then
              u1 = (+spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
            else
              u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
            endif
 
            rin = numglu(1)+numglu(2)+numglu(3)+nc+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
            bosonvertex_mod = bosonvertex - 3
            ubar2 = cur_f_2fv(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,boson,(/nd+ne,nd,ne/))
            pmom3  = summom(gluons,rin,rout) + quarks(4)%Mom + boson%Mom
            propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(4)%Mass2)
            if( abs(sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.propcut ) cycle
            if( quarks(4)%PartType.lt.0 ) then
               ubar2 = (-spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
            else
               ubar2 = (+spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
            endif
 
            if( quarks(4)%PartType.lt.0 ) then
              eps1 = +vbqq(dv,u1,ubar2)       ! re-checked
            else
              eps1 = -vbqq(dv,ubar2,u1)       ! re-checked
            endif
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpmom1(:) = pmom2(:)+pmom3(:)
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
            if( na.ge.1 .or. nf.ge.1 ) then
               propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
               if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
               eps2 = eps2*propfac1
            endif
 
            res = res + eps2
         enddo
         enddo
         enddo
      endif! Bosonvertex
 
 
      if (bosonvertex .ne. 1 .and. bosonvertex .ne. 2 .and. bosonvertex .ne. 3) then
         print *, 'cur_g_4fV not implemented for this flavor choice and BosonVertex=', bosonvertex
         stop
      endif
 
      endif! Flavor check
 
      if( quarks(1)%PartType.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType) then
 
         if (bosonvertex .eq. 1 .or. bosonvertex .eq. 3) then
             do na=0,numglu(1)
             do nc=0,numglu(2)
             do ne=0,numglu(5)
 
                nb=numglu(1)-na
                nd=numglu(2)-nc
                nf=numglu(5)-ne
                rin = numglu(1)+nc+1
                rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
                u1 = cur_f_4fv(gluons(rin:rout),quarks(2:4),quarks(1)%PartType,boson, bosonvertex,(/nd+numglu(3)+numglu(4)+ne,nd,numglu(3),numglu(4),ne/),0,0)
                pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + boson%Mom
                propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(1)%Mass2)
                if( abs(sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.propcut ) cycle
                if( quarks(1)%PartType.lt.0 ) then
                   u1 = ( spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
                else
                   u1 = (-spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
                endif
 
                rin = na+1
                rout= numglu(1)+nc
                ubar2 = cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
                pmom3 = summom(gluons,rin,rout) + quarks(1)%Mom
                if( nb.ge.1 .or. nc.ge.1 ) then
                   propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(1)%Mass2)
                   if( abs(sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.propcut ) cycle
                   if( quarks(1)%PartType.lt.0 ) then
                      ubar2 = (-spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
                   else
                      ubar2 = (+spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
                   endif
                endif
 
                if( quarks(1)%PartType.lt.0 ) then
                   eps1 = -vbqq(dv,u1,ubar2)       ! re-checked
                else
                   eps1 = +vbqq(dv,ubar2,u1)       ! re-checked
                endif
 
                counter=1
                rin =1
                rout=na
                do i=rin,rout
                   call copyparticleptr(gluons(i),tmpgluons(counter))
                   counter=counter+1
                enddo
                tmpmom1(:) = pmom2(:)+pmom3(:)
                tmpextref1 = -1
                tmpgluons(counter)%Mom => tmpmom1(:)
                tmpgluons(counter)%Pol => eps1(:)
                tmpgluons(counter)%ExtRef => tmpextref1
                counter=counter+1
                rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
                rout=numglu(0)-1
                do i=rin,rout
                   call copyparticleptr(gluons(i),tmpgluons(counter))
                   counter=counter+1
                enddo
                eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
                if( na.ge.1 .or. nf.ge.1 ) then
                   propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
                   if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
                   eps2 = eps2*propfac1
                endif
 
                res = res + eps2
             enddo
             enddo
          enddo
       endif
 
       if (bosonvertex .eq. 1) then
          do na=0,numglu(1)
          do nc=0,numglu(2)
          do ne=0,numglu(5)
             nb=numglu(1)-na
             nd=numglu(2)-nc
             nf=numglu(5)-ne
             rin = numglu(1)+nc+1
             rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
             u1 = cur_f_4f(gluons(rin:rout),quarks(2:4),quarks(1)%PartType,(/nd+numglu(3)+numglu(4)+ne,nd,numglu(3),numglu(4),ne/),0,0)
             pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
             propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(1)%Mass2)
             if( abs(sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.propcut ) cycle
             if( quarks(1)%PartType.lt.0 ) then
                u1 = ( spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
             else
                u1 = (-spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
             endif
 
             rin = na+1
             rout= numglu(1)+nc
             ubar2 = cur_f_2fv(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,boson,(/nb+nc,nb,nc/))
             pmom3 = summom(gluons,rin,rout) + quarks(1)%Mom +  boson%Mom
             propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(1)%Mass2)
             if( abs(sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.propcut ) cycle
             if( quarks(1)%PartType.lt.0 ) then
                ubar2 = (-spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
             else
                ubar2 = (+spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
             endif
 
             if( quarks(1)%PartType.lt.0 ) then
                eps1 = -vbqq(dv,u1,ubar2)       ! re-checked
             else
                eps1 = +vbqq(dv,ubar2,u1)       ! re-checked
             endif
 
             counter=1
             rin =1
             rout=na
             do i=rin,rout
                call copyparticleptr(gluons(i),tmpgluons(counter))
                counter=counter+1
             enddo
             tmpmom1(:) = pmom2(:)+pmom3(:)
             tmpextref1 = -1
             tmpgluons(counter)%Mom => tmpmom1(:)
             tmpgluons(counter)%Pol => eps1(:)
             tmpgluons(counter)%ExtRef => tmpextref1
             counter=counter+1
             rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
             rout=numglu(0)-1
             do i=rin,rout
                call copyparticleptr(gluons(i),tmpgluons(counter))
                counter=counter+1
             enddo
             eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
             if( na.ge.1 .or. nf.ge.1 ) then
                propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
                if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
                eps2 = eps2*propfac1
             endif
 
             res = res + eps2
          enddo
       enddo
    enddo
 endif
 
 
 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 3) then
!        type (3)
    do na=0,numglu(1)
    do nc=0,numglu(4)
    do ne=0,numglu(5)        ! can be replaced by above ne-loop
       nb=numglu(1)-na
       nd=numglu(4)-nc
       nf=numglu(5)-ne
 
       rin = na+1
       rout= numglu(1)+numglu(2)+numglu(3)+nc
       bosonvertex_mod=bosonvertex+1
       u1 = cur_f_4fv(gluons(rin:rout),quarks(1:3),quarks(4)%PartType,boson,bosonvertex_mod,(/nb+numglu(2)+numglu(3)+nc,nb,numglu(2),numglu(3),nc/),0,0)
       pmom2 =  summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
       propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
       if ( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) then
          cycle
       endif
       if( quarks(4)%PartType.lt.0 ) then
          u1 = (+spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
       else
          u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
       endif
       rin = numglu(1)+numglu(2)+numglu(3)+nc+1
       rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
       ubar2 = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/nd+ne,nd,ne/))
       pmom3  = summom(gluons,rin,rout) + quarks(4)%Mom
       if( nd.ge.1 .or. ne.ge.1 ) then
          propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(4)%Mass2)
          if( abs(sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.propcut ) cycle
          if( quarks(4)%PartType.lt.0 ) then
             ubar2 = (-spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
          else
             ubar2 = (+spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
          endif
       endif
 
       if( quarks(4)%PartType.lt.0 ) then
          eps1 = +vbqq(dv,u1,ubar2)       ! re-checked
       else
          eps1 = -vbqq(dv,ubar2,u1)       ! re-checked
       endif
 
       counter=1
       rin =1
       rout=na
       do i=rin,rout
          call copyparticleptr(gluons(i),tmpgluons(counter))
          counter=counter+1
       enddo
       tmpmom1(:) = pmom2(:)+pmom3(:)
       tmpextref1 = -1
       tmpgluons(counter)%Mom => tmpmom1(:)
       tmpgluons(counter)%Pol => eps1(:)
       tmpgluons(counter)%ExtRef => tmpextref1
       counter=counter+1
       rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
       rout=numglu(0)-1
       do i=rin,rout
          call copyparticleptr(gluons(i),tmpgluons(counter))
          counter=counter+1
       enddo
       eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
       if( na.ge.1 .or. nf.ge.1 ) then
          propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
          if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
          eps2 = eps2*propfac1
       endif
 
       res = res + eps2
    enddo
    enddo
    enddo
   endif
 
   if (bosonvertex .eq. 1) then
! type(4)
      do na=0,numglu(1)
      do nc=0,numglu(2)
      do ne=0,numglu(3)
      do nf=0,numglu(3)-ne
      do nh=0,numglu(4)   ! this loop can be placed after Eps1 has been calculated
      do nj=0,numglu(5)
         nb=numglu(1)-na
         nd=numglu(2)-nc
         ng=numglu(3)-ne-nf
         ni=numglu(4)-nh
         nk=numglu(5)-nj
 
         rin = na+1
         rout= numglu(1)+numglu(2)+ne
         eps1 = cur_g_2fv(gluons(rin:rout),quarks(1:2),boson,(/1+nb+nc+nd+ne,nc,nc+nd,ne/) )
         tmpmom1 =  summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + boson%Mom
         propfac3 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
         if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
         eps1 = eps1*propfac3
 
         rin = numglu(1)+numglu(2)+ne+nf+1
         rout= numglu(1)+numglu(2)+numglu(3)+nh
         u1 = cur_f_2f(gluons(rin:rout),quarks(3:3),-quarks(3)%PartType,(/ng+nh,ng,nh/))
         pmom3 = summom(gluons,rin,rout) + quarks(3)%Mom
         if( ng.ge.1 .or. nh.ge.1 ) then
            propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(3)%Mass2)
            if( abs(sc_(pmom3,pmom3) - quarks(3)%Mass2).lt.propcut ) cycle
            if( quarks(3)%PartType.lt.0 ) then
               u1 = (-spi2_(pmom3,u1) + quarks(3)%Mass*u1 )*propfac3
            else
               u1 = (+spb2_(u1,pmom3) + quarks(3)%Mass*u1 )*propfac3
            endif
         endif
         rin = numglu(1)+numglu(2)+numglu(3)+nh+1
         rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj
         ubar2 = cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/ni+nj,ni,nj/))
         pmom4 = summom(gluons,rin,rout) + quarks(4)%Mom
         if( ni.ge.1 .or. nj.ge.1 ) then
            propfac4 = (0d0,1d0)/(sc_(pmom4,pmom4) - quarks(4)%Mass2)
            if( abs(sc_(pmom4,pmom4) - quarks(4)%Mass2).lt.propcut ) cycle
            if( quarks(4)%PartType.lt.0 ) then
               ubar2 = (-spi2_(pmom4,ubar2) + quarks(4)%Mass*ubar2 )*propfac4
            else
               ubar2 = (+spb2_(ubar2,pmom4) + quarks(4)%Mass*ubar2 )*propfac4
            endif
         endif
 
         if( quarks(4)%PartType.lt.0 ) then
            eps2 = +vbqq(dv,u1,ubar2)       ! re-checked
         else
            eps2 = -vbqq(dv,ubar2,u1)       ! re-checked
         endif
         tmpmom2 = pmom3 + pmom4
         propfac1 = (0d0,-1d0)/sc_(tmpmom2,tmpmom2)
         if( abs(sc_(tmpmom2,tmpmom2)).lt.propcut ) cycle
         eps2 = eps2*propfac1
 
 
         counter=1
         rin =1
         rout=na
         do i=rin,rout
            call copyparticleptr(gluons(i),tmpgluons(counter))
            counter=counter+1
         enddo
         tmpextref1 = -1
         tmpgluons(counter)%Mom => tmpmom1(:)
         tmpgluons(counter)%Pol => eps1(:)
         tmpgluons(counter)%ExtRef => tmpextref1
         counter=counter+1
         rin =numglu(1)+numglu(2)+ne+1
         rout=numglu(1)+numglu(2)+ne+nf
         do i=rin,rout
            call copyparticleptr(gluons(i),tmpgluons(counter))
            counter=counter+1
         enddo
         tmpgluons(counter)%Mom => tmpmom2(:)
         tmpgluons(counter)%Pol => eps2(:)
         tmpgluons(counter)%ExtRef => tmpextref1
         counter=counter+1
         rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj+1
         rout=numglu(0)-1
         do i=rin,rout
            call copyparticleptr(gluons(i),tmpgluons(counter))
            counter=counter+1
         enddo
         eps3(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+nk+2)
         res = res + eps3
         enddo
         enddo
         enddo
         enddo
         enddo
         enddo
      endif
 
      if (bosonvertex .eq. 3) then
 
         do na=0,numglu(1)
         do nc=0,numglu(4)
         do ne=0,numglu(5)        ! can be replaced by above ne-loop
            nb=numglu(1)-na
            nd=numglu(4)-nc
            nf=numglu(5)-ne
 
            rin = na+1
            rout= numglu(1)+numglu(2)+numglu(3)+nc
            u1 = cur_f_4f(gluons(rin:rout),quarks(1:3),quarks(4)%PartType,(/nb+numglu(2)+numglu(3)+nc,nb,numglu(2),numglu(3),nc/),0,0)
            pmom2 =  summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom
            propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
            if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
            if( quarks(4)%PartType.lt.0 ) then
              u1 = (+spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
            else
              u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
            endif
            rin = numglu(1)+numglu(2)+numglu(3)+nc+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne
            ubar2 = cur_f_2fv(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,boson,(/nd+ne,nd,ne/))
            pmom3  = summom(gluons,rin,rout) + quarks(4)%Mom + boson%Mom
            propfac3 = (0d0,1d0)/(sc_(pmom3,pmom3) - quarks(4)%Mass2)
            if( abs(sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.propcut ) cycle
            if( quarks(4)%PartType.lt.0 ) then
               ubar2 = (-spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
            else
               ubar2 = (+spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
            endif
 
            if( quarks(4)%PartType.lt.0 ) then
              eps1 = +vbqq(dv,u1,ubar2)       ! re-checked
            else
              eps1 = -vbqq(dv,ubar2,u1)       ! re-checked
            endif
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpmom1(:) = pmom2(:)+pmom3(:)
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+ne+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps2(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+1)
 
            if( na.ge.1 .or. nf.ge.1 ) then
               propfac1 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
               if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
               eps2 = eps2*propfac1
            endif
 
            res = res + eps2
         enddo
         enddo
         enddo
      endif
 
      if (bosonvertex .eq. 3) then
 
         do na=0,numglu(1)
         do nc=0,numglu(2)
         do ne=0,numglu(3)
         do nf=0,numglu(3)-ne
         do nh=0,numglu(4)   ! this loop can be placed after Eps1 has been calculated
         do nj=0,numglu(5)
            nb=numglu(1)-na
            nd=numglu(2)-nc
            ng=numglu(3)-ne-nf
            ni=numglu(4)-nh
            nk=numglu(5)-nj
 
            rin = na+1
            rout= numglu(1)+nc
            ubar2 = cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
            pmom1 = summom(gluons,rin,rout) + quarks(1)%Mom
            if( nb.ge.1 .or. nc.ge.1 ) then
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1) - quarks(1)%Mass2)
               if( abs(sc_(pmom1,pmom1) - quarks(1)%Mass2).lt.propcut ) cycle
               if( quarks(1)%PartType.lt.0 ) then
                  ubar2 = (-spi2_(pmom1,ubar2) + quarks(1)%Mass*ubar2 )*propfac1
               else
                  ubar2 = (+spb2_(ubar2,pmom1) + quarks(1)%Mass*ubar2 )*propfac1
               endif
            endif
            rin = numglu(1)+nc+1
            rout= numglu(1)+numglu(2)+ne
            u1 = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/nd+ne,nd,ne/))
            pmom2 = summom(gluons,rin,rout) + quarks(2)%Mom
            if( nd.ge.1 .or. ne.ge.1 ) then
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
               if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  u1 = (-spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
               else
                  u1 = (+spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
               endif
            endif
 
            if( quarks(2)%PartType.lt.0 ) then
              eps1 = +vbqq(dv,ubar2,u1)       ! re-checked
            else
              eps1 = -vbqq(dv,u1,ubar2)       ! re-checked
            endif
            tmpmom1 = pmom1 + pmom2
            propfac3 = (0d0,-1d0)/sc_(tmpmom1,tmpmom1)
            if( abs(sc_(tmpmom1,tmpmom1)).lt.propcut ) cycle
            eps1 = eps1*propfac3
 
            rin = numglu(1)+numglu(2)+ne+nf+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj
            eps2=cur_g_2fv(gluons(rin:rout),quarks(3:4),boson,(/1+ng+nh+ni+nj,ng,nh,ni,nj /) )
            tmpmom2 = summom(gluons,rin,rout) +quarks(3)%Mom + quarks(4)%Mom + boson%Mom
            propfac1 = (0d0,-1d0)/sc_(tmpmom2,tmpmom2)
            if( abs(sc_(tmpmom2,tmpmom2)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            counter=1
            rin =1
            rout=na
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpextref1 = -1
            tmpgluons(counter)%Mom => tmpmom1(:)
            tmpgluons(counter)%Pol => eps1(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+ne+1
            rout=numglu(1)+numglu(2)+ne+nf
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmpgluons(counter)%Mom => tmpmom2(:)
            tmpgluons(counter)%Pol => eps2(:)
            tmpgluons(counter)%ExtRef => tmpextref1
            counter=counter+1
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj+1
            rout=numglu(0)-1
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            eps3(:) = cur_g(tmpgluons(1:counter-1),1+na+nf+nk+2)
            res = res + eps3
 
         enddo
      enddo
   enddo
enddo
enddo
enddo
 
endif   ! BosonVertex
 
 
if (bosonvertex .ne. 1 .and. bosonvertex .ne. 3) then
   print *, 'Not implemented for for cur_g_4fV for this flavor structure and BosonVertex=', bosonvertex
   stop
endif
endif
 
return
END FUNCTION
 
 
 
 
 
 
 
 
 
 
 
FUNCTION cur_f_6fv(Gluons,Quarks,Quark1PartType,Boson,BosonVertex,NumGlu,tag_f) result(res)           ! Quarks(:) does not include the OFF-shell quark
implicit none
integer :: numglu(0:6),quark1parttype,tag_f,bosonvertex,bosonvertex_mod
type(ptrtoparticle) :: gluons(1:),quarks(2:6),boson
integer,target :: tmpparttype,tmpextref
complex(8) :: res(1:ds),tmp(1:ds)
! complex(8) :: Res1(1:Ds),Res2(1:Ds),Res3(1:Ds),Res4(1:Ds)
complex(8) :: u1(1:ds),ubar1(1:ds)
complex(8),target :: ubar0(1:ds)
complex(8) :: eps1(1:dv)
complex(8) :: eps2(1:dv)
type(ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(6)),tmpquark(1:1)
complex(8) :: propfac1,propfac2
complex(8),target :: pmom1(1:dv)
complex(8),target :: pmom2(1:dv)
integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b,n5a,n5b,n6a,n6b
integer :: rin,rout,i,counter
 
 
!DEC$ IF (_DebugCheckMyImpl1==1)
    if( numglu(0)-numglu(1)-numglu(2)-numglu(3)-numglu(4)-numglu(5)-numglu(6).ne.0 ) print *, "wrong number of gluons in cur_f_6fV"
!DEC$ ENDIF
 
    res = (0d0,0d0)
 
 
! -- not all possible choices of flavors and BosonVertex values are used in ttb+Z calculation. This warns you if you try to use something new and potentially buggy.
    if (quark1parttype.eq.-quarks(2)%PartType) then
       if ( ((quarks(3)%PartType.eq.-quarks(4)%PartType) .AND. (bosonvertex .eq. 2 .OR.bosonvertex .eq. 4 .OR.bosonvertex .eq. 6)) .OR. &
         ((quarks(3)%PartType.eq.-quarks(6)%PartType) .AND. (bosonvertex .eq. 2 .OR.bosonvertex .eq. 6)) ) then
          print *, 'WARNING :  cur_f_6fV with this flavor structure and this choice of BosonVertex is implemented, but not checked!'
          print *, 'Quark flavors: ',quark1parttype, quarks(2)%PartType,quarks(3)%PartType,quarks(4)%PartType,quarks(5)%PartType,quarks(6)%PartType
          print *, 'BosonVertex =', bosonvertex
       endif
    elseif (quark1parttype.eq.-quarks(6)%PartType ) then
       if ( ((quark1parttype.eq.-quarks(2)%PartType) .AND. (bosonvertex .eq. 2 .or. bosonvertex .eq. 4 .or. bosonvertex .eq. 6)) .OR. &
            & ((quark1parttype.eq.-quarks(4)%PartType) .AND. (bosonvertex .eq. 4 .or. bosonvertex .eq. 6)) ) then
          print *, 'WARNING :  cur_f_6fV with this flavor structure and this choice of BosonVertex is implemented, but not checked!'
          print *, 'Quark flavors: ',quark1parttype, quarks(2)%PartType,quarks(3)%PartType,quarks(4)%PartType,quarks(5)%PartType,quarks(6)%PartType
          print *, 'BosonVertex =', bosonvertex
       endif
    elseif (quarks(2)%PartType.eq.-quarks(3)%PartType) then
       if ( ((quark1parttype.eq.-quarks(4)%PartType) .AND.  (bosonvertex .eq. 4 .or. bosonvertex .eq. 6) ) .OR. &
            ((quark1parttype.eq.-quarks(6)%PartType) .AND.  (bosonvertex .eq. 4 .or. bosonvertex .eq. 6) ) ) then
          print *, 'WARNING :  cur_f_6fV with this flavor structure and this choice of BosonVertex is implemented, but not checked!'
          print *, 'Quark flavors: ',quark1parttype, quarks(2)%PartType,quarks(3)%PartType,quarks(4)%PartType,quarks(5)%PartType,quarks(6)%PartType
          print *, 'BosonVertex =', bosonvertex
       endif
    endif
 
 
! probably some tag_Z checks are needed here: not yet implemented
 
 
!   (A)
    if( quark1parttype.eq.-quarks(2)%PartType .AND. (quarks(3)%PartType.eq.-quarks(4)%PartType .or. quarks(3)%PartType.eq.-quarks(6)%PartType) &
         .AND. (quarks(2)%ExtRef.ne.-1 .or. tag_f.ne.1) &
         .AND. ( bosonvertex.eq.1 .OR. bosonvertex.eq.2) &
         ) then
!       print *, 'A-1'
      do n2a=0,numglu(2)
      do n6a=0,numglu(6)
         n2b = numglu(2)-n2a
         n6b = numglu(6)-n6a
 
         rin =numglu(1)+n2a+1
         rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
         eps2 = cur_g_4f(gluons(rin:rout),quarks(3:6),(/1+n2b+numglu(3)+numglu(4)+numglu(5)+n6a,n2b,numglu(3),numglu(4),numglu(5),n6a/))
         pmom1(:) = summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
         propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
         if( abs(sc_(pmom1,pmom1)).lt.propcut) cycle
         eps2 = eps2*propfac1
         do n1a=0,numglu(1)
            n1b = numglu(1)-n1a
            ! Fer2 couple V on the top
            rin =n1a+1
            rout=numglu(1)+n2a
            if (bosonvertex .eq. 1 .or. bosonvertex .eq. 2) then
               ubar1(:) = cur_f_2fv(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,boson,(/n1b+n2a,n1b,n2a/) )
               pmom2(:) = quarks(2)%Mom + summom(gluons,rin,rout) + boson%Mom
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
               if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) then
                  propfac2=(0d0,0d0)
               endif
 
               if( quarks(2)%PartType.lt.0 ) then
                  ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
               else
                  ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
               endif
               if( quarks(2)%PartType.lt.0 ) then
                  ubar0(:) = vbqg(ubar1,eps2)
               else
                  ubar0(:) = vqg(ubar1,eps2)
               endif
 
               rin = n1a+1
               rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
               pmom1 = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom + boson%Mom  ! can be simplified with PMom1(:)
               if(n1a.ge.1 .or. n6b.ge.1) then
                  propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
                  if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) then
                     propfac1=(0d0,0d0)
                  endif
 
                  if( quarks(2)%PartType.lt.0 ) then
                     ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
                  else
                     ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
                  endif
               endif
 
               tmpquark(1)%Mom  => pmom1(:)
               tmpquark(1)%Pol  => ubar0(:)
               tmpquark(1)%Mass => quarks(2)%Mass
               tmpquark(1)%Mass2=> quarks(2)%Mass2
               tmpextref = -1
               tmpquark(1)%ExtRef => tmpextref
               tmpquark(1)%PartType => quarks(2)%PartType
               counter=1
               rin =1
               rout=n1a
               do i=rin,rout
                  call copyparticleptr(gluons(i),tmpgluons(counter))
                  counter=counter+1
               enddo
               rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
               rout=numglu(0)
               do i=rin,rout
                  call copyparticleptr(gluons(i),tmpgluons(counter))
                  counter=counter+1
               enddo
               tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n6b/) )
               res(:) = res(:) + tmp(:)
            endif
            if (bosonvertex .eq. 1 .or. bosonvertex .eq. 6) then
               ! Fer2  couple V on the bottom
               rin =n1a+1
               rout=numglu(1)+n2a
               ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n1b+n2a,n1b,n2a/) )
               if(n1b.ge.1 .or. n2a.ge.1) then
                  pmom2(:) = quarks(2)%Mom + summom(gluons,rin,rout)
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
                  if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) cycle
                  if( quarks(2)%PartType.lt.0 ) then
                     ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
                  else
                     ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
                  endif
               endif
               if( quarks(2)%PartType.lt.0 ) then
                  ubar0(:) = vbqg(ubar1,eps2)       ! re-checked
               else
                  ubar0(:) = vqg(ubar1,eps2)       ! re-checked
               endif
 
               rin = n1a+1
               rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
               pmom1 = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
 
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
               if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
               else
                  ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
                  endif
 
                  tmpquark(1)%Mom  => pmom1(:)
                  tmpquark(1)%Pol  => ubar0(:)
                  tmpquark(1)%Mass => quarks(2)%Mass
                  tmpquark(1)%Mass2=> quarks(2)%Mass2
                  tmpextref = -1
                  tmpquark(1)%ExtRef => tmpextref
                  tmpquark(1)%PartType => quarks(2)%PartType
                  counter=1
                  rin =1
                  rout=n1a
                  do i=rin,rout
                     call copyparticleptr(gluons(i),tmpgluons(counter))
                     counter=counter+1
                  enddo
                  rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
                  rout=numglu(0)
                  do i=rin,rout
                     call copyparticleptr(gluons(i),tmpgluons(counter))
                     counter=counter+1
                  enddo
                  tmp(:) = cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,boson,(/counter-1,n1a,n6b/) )
                  res(:) = res(:) + tmp(:)
               endif
            enddo
         enddo
      enddo
   endif
 
 
 
!   (A)
    if( quark1parttype.eq.-quarks(2)%PartType .AND. (quarks(3)%PartType.eq.-quarks(4)%PartType .or. quarks(3)%PartType.eq.-quarks(6)%PartType) &
        .AND. (quarks(2)%ExtRef.ne.-1 .or. tag_f.ne.1) &
        .AND. ( bosonvertex.eq.3 .OR. bosonvertex .eq. 4 .OR. bosonvertex.eq.5 ) &
      ) then
!       print *, 'A-2'
      do n2a=0,numglu(2)
      do n6a=0,numglu(6)
         n2b = numglu(2)-n2a
         n6b = numglu(6)-n6a
 
         rin =numglu(1)+n2a+1
         rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
         bosonvertex_mod = bosonvertex - 2
         eps2 = cur_g_4fv(gluons(rin:rout),quarks(3:6),boson,bosonvertex_mod,(/1+n2b+numglu(3)+numglu(4)+numglu(5)+n6a,n2b,numglu(3),numglu(4),numglu(5),n6a/))
         pmom1(:) = summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
         propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
         if( abs(sc_(pmom1,pmom1)).lt.propcut) cycle
         eps2 = eps2*propfac1
         do n1a=0,numglu(1)
            n1b = numglu(1)-n1a
            ! Fer2
            rin =n1a+1
            rout=numglu(1)+n2a
            ubar1(:) = cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n1b+n2a,n1b,n2a/) )
            if(n1b.ge.1 .or. n2a.ge.1) then
               pmom2(:) = quarks(2)%Mom + summom(gluons,rin,rout)
               propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2)-quarks(2)%Mass2)
               if( abs(sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  ubar1(:) = (-spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
               else
                  ubar1(:) = (+spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
               endif
            endif
            if( quarks(2)%PartType.lt.0 ) then
               ubar0(:) = vbqg(ubar1,eps2)
            else
               ubar0(:) = vqg(ubar1,eps2)
            endif
 
            rin = n1a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            pmom1 = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
            if(n1a.ge.1 .or. n6b.ge.1) then
               propfac1 = (0d0,1d0)/(sc_(pmom1,pmom1)-quarks(2)%Mass2)
               if( abs(sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.propcut ) cycle
               if( quarks(2)%PartType.lt.0 ) then
                  ubar0(:) = (-spi2_(pmom1,ubar0)+quarks(2)%Mass*ubar0(:))*propfac1
               else
                  ubar0(:) = (+spb2_(ubar0,pmom1)+quarks(2)%Mass*ubar0(:))*propfac1
               endif
            endif
 
            tmpquark(1)%Mom  => pmom1(:)
            tmpquark(1)%Pol  => ubar0(:)
            tmpquark(1)%Mass => quarks(2)%Mass
            tmpquark(1)%Mass2=> quarks(2)%Mass2
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpquark(1)%PartType => quarks(2)%PartType
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n6b/) )
            res(:) = res(:) + tmp(:)
 
         enddo
      enddo
      enddo
   endif
 
 
 
!   (B)
!! RR -- I think that only times this will be called in ttbZ will be color-zero in any case (tadpoles), so I'm going to just comment it out.
!! I imagine that if implemented correctly, then the it would return zero for ttbZ for kinematic reasons, but there doesn't seem to be any poin in coding it up for now...
 
!    if( Quark1PartType.eq.-Quarks(6)%PartType .AND. (Quarks(2)%PartType.eq.-Quarks(5)%PartType .or. Quarks(2)%PartType.eq.-Quarks(3)%PartType) &
!        .AND. (Quarks(6)%ExtRef.ne.-1 .or. tag_f.ne.1) &
!      ) then
!      call Error("This 6fV(B) current is not yet implemented")
!
!    endif
 
!   (C)
    if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
        ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
    .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
    .AND. ( bosonvertex.eq.1 .OR. bosonvertex.eq.6 )  &
      ) then
!       print *, 'C-1'
      do n1a=0,numglu(1)
      do n4a=0,numglu(4)
      do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n4b = numglu(4)-n4a
         n6b = numglu(6)-n6a
 
            rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            eps2 = cur_g_2f(gluons(rin:rout),quarks(5:6),(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
            pmom1(:) = summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            rin = n1a+1
            rout= numglu(1)+numglu(2)+numglu(3)+n4a
            u1 = cur_f_4f(gluons(rin:rout),quarks(2:4),quark1parttype,(/n1b+numglu(2)+numglu(3)+n4a,n1b,numglu(2),numglu(3),n4a/),0,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
 
            if( quark1parttype.eq.-quarks(2)%PartType) then
                 if( quarks(2)%PartType.lt.0 ) then
                   propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                   if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                   u1 = (-spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
                   ubar0 = vbqg(u1,eps2)       ! re-checked
                   rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                   pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom   ! this PMom2 will be re-used below
 
                   if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(2)%Mass*ubar0 )*propfac2
                   endif
                elseif( quarks(2)%PartType.gt.0 ) then
                   propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                   if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                   u1 = ( spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
                   ubar0 = vqg(u1,eps2)       ! re-checked
                   rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                   pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
 
                   if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                       if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                       ubar0 = ( spb2_(ubar0,pmom2) + quarks(2)%Mass*ubar0 )*propfac2
                    endif
                 endif
                 tmpquark(1)%Mom  => pmom2(:)
                 tmpquark(1)%Pol  => ubar0(:)
                 tmpquark(1)%Mass => quarks(2)%Mass
                 tmpquark(1)%Mass2=> quarks(2)%Mass2
            elseif( quark1parttype.eq.-quarks(4)%PartType) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vbqg(u1,eps2)
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vqg(u1,eps2)
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(4)%Mass
                tmpquark(1)%Mass2=> quarks(4)%Mass2
            endif
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,boson,(/counter-1,n1a,n6b/) )
 
            res(:) = res(:) + tmp(:)
 
      enddo
      enddo
      enddo
    endif
 
 
 
    if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
        ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
    .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
    .AND. ( bosonvertex.eq.1 .OR. bosonvertex.eq.2 .OR. bosonvertex.eq.3 .OR. bosonvertex.eq.4  )  &
      ) then
!       print *, 'C-2'
      do n1a=0,numglu(1)
      do n4a=0,numglu(4)
      do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n4b = numglu(4)-n4a
         n6b = numglu(6)-n6a
 
            rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            eps2 = cur_g_2f(gluons(rin:rout),quarks(5:6),(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
            pmom1(:) = summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            rin = n1a+1
            rout= numglu(1)+numglu(2)+numglu(3)+n4a
            bosonvertex_mod = bosonvertex
            u1 = cur_f_4fv(gluons(rin:rout),quarks(2:4),quark1parttype,boson,bosonvertex_mod,(/n1b+numglu(2)+numglu(3)+n4a,n1b,numglu(2),numglu(3),n4a/),0,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + boson%Mom
 
            if( quark1parttype.eq.-quarks(2)%PartType) then
                 if( quarks(2)%PartType.lt.0 ) then
                   propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                   if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                   u1 = (-spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
                   ubar0 = vbqg(u1,eps2)       ! re-checked
                   rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                   pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom   ! this PMom2 will be re-used below
                   if( n1a.ge.1 .or. n6b.ge.1 ) then
                       propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                       if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                       ubar0 = (-spi2_(pmom2,ubar0) + quarks(2)%Mass*ubar0 )*propfac2
                   endif
                 elseif( quarks(2)%PartType.gt.0 ) then
                   propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                   if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                   u1 = ( spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
                   ubar0 = vqg(u1,eps2)       ! re-checked
                   rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                   pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                   if( n1a.ge.1 .or. n6b.ge.1 ) then
                       propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                       if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                       ubar0 = ( spb2_(ubar0,pmom2) + quarks(2)%Mass*ubar0 )*propfac2
                   endif
                 endif
                 tmpquark(1)%Mom  => pmom2(:)
                 tmpquark(1)%Pol  => ubar0(:)
                 tmpquark(1)%Mass => quarks(2)%Mass
                 tmpquark(1)%Mass2=> quarks(2)%Mass2
            elseif( quark1parttype.eq.-quarks(4)%PartType) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vbqg(u1,eps2)
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vqg(u1,eps2)
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(4)%Mass
                tmpquark(1)%Mass2=> quarks(4)%Mass2
            endif
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
 
            res(:) = res(:) + tmp(:)
 
      enddo
      enddo
      enddo
    endif
 
 
!   (C)
    if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
        ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
    .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
    .AND. ( bosonvertex.eq.5 )  &
      ) then
!       print *, 'C-3'
       do n1a=0,numglu(1)
       do n4a=0,numglu(4)
       do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n4b = numglu(4)-n4a
         n6b = numglu(6)-n6a
 
            rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            eps2 = cur_g_2fv(gluons(rin:rout),quarks(5:6),boson,(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
            pmom1(:) = summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            rin = n1a+1
            rout= numglu(1)+numglu(2)+numglu(3)+n4a
            u1 = cur_f_4f(gluons(rin:rout),quarks(2:4),quark1parttype,(/n1b+numglu(2)+numglu(3)+n4a,n1b,numglu(2),numglu(3),n4a/),0,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
 
            if( quark1parttype.eq.-quarks(2)%PartType) then
 
                 if( quarks(2)%PartType.lt.0 ) then
                   propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                   if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                   u1 = (-spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
                   ubar0 = vbqg(u1,eps2)       ! re-checked
                   rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                   pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom   + boson%Mom! this PMom2 will be re-used below
                   if( n1a.ge.1 .or. n6b.ge.1 ) then
                       propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                       if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                       ubar0 = (-spi2_(pmom2,ubar0) + quarks(2)%Mass*ubar0 )*propfac2
                   endif
                 elseif( quarks(2)%PartType.gt.0 ) then
                    propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                    if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                    u1 = ( spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
                    ubar0 = vqg(u1,eps2)       ! re-checked
                   rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                   rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                   pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom   + boson%Mom! this PMom2 will be re-used below
                   if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(2)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.propcut ) cycle
                      u1 = ( spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
                   endif
                endif
                 tmpquark(1)%Mom  => pmom2(:)
                 tmpquark(1)%Pol  => ubar0(:)
                 tmpquark(1)%Mass => quarks(2)%Mass
                 tmpquark(1)%Mass2=> quarks(2)%Mass2
            elseif( quark1parttype.eq.-quarks(4)%PartType) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vbqg(u1,eps2)
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vqg(u1,eps2)
                  rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
                  rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(4)%Mass
                tmpquark(1)%Mass2=> quarks(4)%Mass2
            endif
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
 
            res(:) = res(:) + tmp(:)
         enddo
      enddo
   enddo
 
 
endif
 
 
 
!   (D)
    if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
        (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
   .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1))  )  &
   .AND. ( bosonvertex.eq.1  )  &
     ) then
!       print *, 'D-1'
      do n1a=0,numglu(1)
      do n3a=0,numglu(3)
      do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n3b = numglu(3)-n3a
         n6b = numglu(6)-n6a
 
            counter=1
            rin = n1a+1
            rout= numglu(1)+numglu(2)+n3a
            eps2 = cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
            pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            rin = numglu(1)+numglu(2)+n3a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            u1 = cur_f_4f(gluons(rin:rout),quarks(4:6),quark1parttype,(/n3b+numglu(4)+numglu(5)+n6a,n3b,numglu(4),numglu(5),n6a/),tag_f,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
 
            if( quark1parttype.eq.-quarks(4)%PartType ) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom   !this PMom2 will be re-used below  ! CAN be written as PMom2=PMom2+PMom1
                  if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1  .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(4)%Mass
                tmpquark(1)%Mass2=> quarks(4)%Mass2
            elseif(quark1parttype.eq.-quarks(6)%PartType) then
               if( quarks(6)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(6)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1  .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                endif
            endif
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,boson,(/counter-1,n1a,n6b/))
 
            res(:)  = res(:) + tmp(:)
 
      enddo
      enddo
      enddo
   endif
 
 
 
!   (D)
    if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
        (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
   .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1))  ) &
   .AND. ( bosonvertex.eq.2)  &
     ) then
!       print *, 'D-2'
      do n1a=0,numglu(1)
      do n3a=0,numglu(3)
      do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n3b = numglu(3)-n3a
         n6b = numglu(6)-n6a
 
            counter=1
            rin = n1a+1
            rout= numglu(1)+numglu(2)+n3a
            eps2 = cur_g_2fv(gluons(rin:rout),quarks(2:3),boson,(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
            pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom +  boson%Mom
            propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
            if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
            eps2 = eps2*propfac1
 
            rin = numglu(1)+numglu(2)+n3a+1
            rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
            bosonvertex_mod = bosonvertex
            u1 = cur_f_4f(gluons(rin:rout),quarks(4:6),quark1parttype,(/n3b+numglu(4)+numglu(5)+n6a,n3b,numglu(4),numglu(5),n6a/),tag_f,0)
            pmom2  = summom(gluons,rin,rout)  + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
 
            if( quark1parttype.eq.-quarks(4)%PartType ) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom +  boson%Mom   !this PMom2 will be re-used below  ! CAN be written as PMom2=PMom2+PMom1
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
            elseif(quark1parttype.eq.-quarks(6)%PartType) then
               if( quarks(6)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom+  boson%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(6)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                endif
            endif
            tmpquark(1)%Mom  => pmom2(:)
            tmpquark(1)%Pol  => ubar0(:)
            tmpquark(1)%Mass => quarks(4)%Mass
            tmpquark(1)%Mass2=> quarks(4)%Mass2
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/))
 
            res(:)  = res(:) + tmp(:)
      enddo
      enddo
      enddo
    endif
 
 
!   (D)
    if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
        (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
   .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1))  ) &
   .AND. ( bosonvertex .ge. 3 .and. bosonvertex .le. 6  )  ) then
!       print *, 'D-3'
       do n1a=0,numglu(1)
       do n3a=0,numglu(3)
       do n6a=0,numglu(6)
         n1b = numglu(1)-n1a
         n3b = numglu(3)-n3a
         n6b = numglu(6)-n6a
 
         counter=1
         rin = n1a+1
         rout= numglu(1)+numglu(2)+n3a
         eps2 = cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
         pmom1(:) = summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
         propfac1 = (0d0,-1d0)/sc_(pmom1,pmom1)
         if( abs(sc_(pmom1,pmom1)).lt.propcut ) cycle
         eps2 = eps2*propfac1
 
         rin = numglu(1)+numglu(2)+n3a+1
         rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
         bosonvertex_mod=bosonvertex-2
         u1 = cur_f_4fv(gluons(rin:rout),quarks(4:6),quark1parttype,boson,bosonvertex_mod,(/n3b+numglu(4)+numglu(5)+n6a,n3b,numglu(4),numglu(5),n6a/),tag_f,0)
         pmom2  = summom(gluons,rin,rout)  + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom+boson%Mom
 
            if( quark1parttype.eq.-quarks(4)%PartType ) then
                if( quarks(4)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom   !this PMom2 will be re-used below  ! CAN be written as PMom2=PMom2+PMom1
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(4)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(4)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(4)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(4)%Mass
                tmpquark(1)%Mass2=> quarks(4)%Mass2
            elseif(quark1parttype.eq.-quarks(6)%PartType) then
                if( quarks(6)%PartType.lt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = (-spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgbq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = (-spi2_(pmom2,ubar0) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                elseif( quarks(6)%PartType.gt.0 ) then
                  propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                  if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                  u1 = ( spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
                  ubar0 = vgq(eps2,u1)       ! re-checked
                  rin = n1a+1
                  rout= numglu(1)+numglu(2)+n3a
                  pmom2 = pmom2 + summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
                  if( n1a.ge.1 .or. n6b.ge.1 ) then
                      propfac2 = (0d0,1d0)/(sc_(pmom2,pmom2) - quarks(6)%Mass2)
                      if( abs(sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.propcut ) cycle
                      ubar0 = ( spb2_(ubar0,pmom2) + quarks(6)%Mass*ubar0 )*propfac2
                  endif
                endif
                tmpquark(1)%Mom  => pmom2(:)
                tmpquark(1)%Pol  => ubar0(:)
                tmpquark(1)%Mass => quarks(6)%Mass
                tmpquark(1)%Mass2=> quarks(6)%Mass2
            endif
            tmpextref = -1
            tmpquark(1)%ExtRef => tmpextref
            tmpparttype = -quark1parttype
            tmpquark(1)%PartType => tmpparttype
            counter=1
            rin =1
            rout=n1a
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
            rout=numglu(0)
            do i=rin,rout
              call copyparticleptr(gluons(i),tmpgluons(counter))
              counter=counter+1
            enddo
            tmp(:) = cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/))
 
            res(:)  = res(:) + tmp(:)
      enddo
      enddo
      enddo
    endif
 
 
return
END FUNCTION
 
 
 
 
 
 
!---------------------------------------
 
 
FUNCTION summom(Particles,i1,i2)
implicit none
complex(8) :: summom(1:dv)
type(ptrtoparticle) :: particles(:)
integer :: i1,i2,j
 
   summom(1:dv)= (0d0,0d0)
   if (i2.ge.i1) then
      do j=i1,i2
         summom(1:dv) = summom(1:dv) + particles(j)%Mom(1:dv)
      enddo
   endif
END FUNCTION
 
 
 
FUNCTION fourvecdot(p1,p2)
implicit none
complex(8), intent(in) :: p1(1:dv),p2(1:dv)
complex(8)  :: fourvecdot
integer :: mu
 
   fourvecdot = p1(1)*p2(1)
!DEC$ UNROLL
   do mu=2,dv
      fourvecdot = fourvecdot - p1(mu)*p2(mu)
   enddo
return
END FUNCTION fourvecdot
 
FUNCTION eval_tripvert(k1,k2,v1,v2)
implicit none
complex(8) :: eval_tripvert(1:dv)
complex(8) :: k1(1:dv),k2(1:dv),v1(1:dv),v2(1:dv)
complex(8), parameter :: ioversqrt2=(0d0,1d0)/dsqrt(2d0)
 
   eval_tripvert(1:dv) = ioversqrt2 * ( (k1(1:dv)-k2(1:dv))*(v1.ndot.v2)  &
                           - 2d0*v1(1:dv)*(k1.ndot.v2)  &
                           + 2d0*v2(1:dv)*(k2.ndot.v1) )
return
END FUNCTION eval_tripvert
 
FUNCTION eval_quadvert(k1,k2,k3)
implicit none
complex(8) :: eval_quadvert(1:dv)
complex(8) :: k1(1:dv),k2(1:dv),k3(1:dv)
complex(8), parameter :: i=(0d0,1d0)
 
   eval_quadvert(1:dv) = i * (-k1(1:dv)*(k2.ndot.k3)*0.5d0  &
                           + k2(1:dv)*(k1.ndot.k3)  &
                           - k3(1:dv)*(k1.ndot.k2)*0.5d0 )
return
END FUNCTION eval_quadvert
 
 
 
 
!DEC$ ATTRIBUTES INLINE :: linear_map
FUNCTION linear_map(i1,i2,Ngluons)
implicit none
integer :: linear_map,i1,i2,ngluons
 
   linear_map = i2+ngluons*(i2-i1)-((i2-i1)*(i2-i1+1))/2
return
END FUNCTION linear_map
 
 
 
 
 
 
SUBROUTINE copyparticleptr(InPointer,OutPointer)
implicit none
type(ptrtoparticle), intent(in) :: InPointer
type(ptrtoparticle), intent(out):: OutPointer
 
   outpointer%PartType => inpointer%PartType
   outpointer%ExtRef => inpointer%ExtRef
   outpointer%Mass => inpointer%Mass
   outpointer%Mass2 => inpointer%Mass2
   outpointer%Helicity => inpointer%Helicity
   outpointer%Mom => inpointer%Mom
   outpointer%Pol => inpointer%Pol
return
END SUBROUTINE
 
 
 
 
 
 
 
!---------------------------------------
 
 
 
 
!-----------modified procedure
!---------- Dv is the dimensionality of the vector space
!---------- Ds is the dimensionality of the spinorial representation
         subroutine spi2(Dv,Ds,v,sp,f)
         implicit none
         integer i,i1,i2,i3,imax,Dv,Ds
         double complex sp(Ds),v(Dv),f(Ds)
         double complex x0(4,4),xx(4,4),xy(4,4)
         double complex xz(4,4),x5(4,4)
         double complex y1,y2,y3,y4,bp,bm,cp,cm
         double complex test(Ds)
 
 
         imax = ds/4
 
           do i=1,imax
           i1= 1+4*(i-1)
           i2=i1+3
 
           y1=sp(i1)
           y2=sp(i1+1)
           y3=sp(i1+2)
           y4=sp(i1+3)
 
           x0(1,i)=y1
           x0(2,i)=y2
           x0(3,i)=-y3
           x0(4,i)=-y4
 
 
           xx(1,i) = y4
           xx(2,i) = y3
           xx(3,i) = -y2
           xx(4,i) = -y1
 
 
           xy(1,i)=dcmplx(0d0,-1d0)*y4
           xy(2,i)=dcmplx(0d0,1d0)*y3
           xy(3,i)=dcmplx(0d0,1d0)*y2
           xy(4,i)=dcmplx(0d0,-1d0)*y1
 
           xz(1,i)=y3
           xz(2,i)=-y4
           xz(3,i)=-y1
           xz(4,i)=y2
 
           x5(1,i)=y3
           x5(2,i)=y4
           x5(3,i)=y1
           x5(4,i)=y2
 
           enddo
 
           if(dv.eq.4) then
 
           do i=1,4
 
           f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)
           enddo
 
           endif
 
 
          if (dv.eq.6) then
           bp = (v(5)+dcmplx(0d0,1d0)*v(6))
           bm=(v(5)-dcmplx(0d0,1d0)*v(6))
 
 
           do i=1,4
 
           f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)-bp*x5(i,2)
 
            i1=i+4
 
            f(i1)=v(1)*x0(i,2)-v(2)*xx(i,2)-v(3)*xy(i,2)-v(4)*xz(i,2)+bm*x5(i,1)
 
            enddo
 
          endif
 
          if (dv.eq.8) then
 
           bp = (v(5)+dcmplx(0d0,1d0)*v(6))
           bm=(v(5)-dcmplx(0d0,1d0)*v(6))
           cp=(v(7)+dcmplx(0d0,1d0)*v(8))
           cm=(v(7)-dcmplx(0d0,1d0)*v(8))
 
 
 
           do i=1,4
 
           f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)-bp*x5(i,2)+ cp*x5(i,3)
 
            i1=i+4
 
            f(i1)=v(1)*x0(i,2)-v(2)*xx(i,2)-v(3)*xy(i,2)-v(4)*xz(i,2)+bm*x5(i,1)-cp*x5(i,4)
 
             i2=i1+4
 
            f(i2)=v(1)*x0(i,3)-v(2)*xx(i,3)-v(3)*xy(i,3)-v(4)*xz(i,3)-bp*x5(i,4)-cm*x5(i,1)
 
             i3=i2+4
 
             f(i3)=v(1)*x0(i,4)-v(2)*xx(i,4)-v(3)*xy(i,4)-v(4)*xz(i,4)+bm*x5(i,3)+cm*x5(i,2)
 
 
            enddo
 
            endif
 
! if(Ds.eq.16) then
!   test(:) = VSpiL(v,sp)
!   print *, ""
!   print *, "spi2",f(1:Ds)
!   print *, "test",test(1:Ds)
!   print *, "diff",test(1:Ds)-f(1:Ds)
!   if( any( abs(test(1:Ds)-f(1:Ds) ).gt.1d-10  ) ) pause
! endif
 
           return
           end subroutine
 
 
 
         function spi2_(v,sp)
         implicit none
         double complex, intent(in) :: sp(:),v(:)
         double complex :: spi2_(size(sp)) ,tmp(size(sp))
         integer :: dv,ds
 
          ds = size(sp)
          if (ds == 4) dv = 4
          if (ds == 8) dv = 6
          if (ds == 16) dv = 8
          call spi2(dv,ds,v,sp,spi2_)
 
!  tmp = VSpiL(v,sp)
!  print *, "diff1",tmp-spi2_
! pause
 
        return
        end function
 
 
 
      function spivl(sp,v)   ! SpiVL=sp.(v_mu*gamma^mu) =spb2(4,4,...)
      implicit none
      double complex :: sp(:),v(:)
      double complex :: spivl(size(sp))
 
            spivl(1) = sp(1)*v(1) + sp(4)*(v(2) + (0d0,1d0)*v(3)) + sp(3)*v(4)
            spivl(2) = sp(2)*v(1) + sp(3)*(v(2) - (0d0,1d0)*v(3)) - sp(4)*v(4)
            spivl(3) =-sp(3)*v(1) - sp(2)*(v(2) + (0d0,1d0)*v(3)) - sp(1)*v(4)
            spivl(4) =-sp(4)*v(1) - sp(1)*(v(2) - (0d0,1d0)*v(3)) + sp(2)*v(4)
 
      return
      end function
 
 
 
subroutine spb2(Dv,Ds,sp,v,f)
         implicit none
         integer i,i1,i2,i3,Dv,Ds,imax
         double complex sp(Ds),v(Dv),f(Ds)
         double complex x0(4,4),xx(4,4),xy(4,4)
         double complex xz(4,4),x5(4,4)
         double complex y1,y2,y3,y4,bp,bm,cp,cm
         double complex test(Ds)
 
           imax = ds/4
 
           do i=1,imax
           i1= 1+4*(i-1)
           i2=i1+3
 
           y1=sp(i1)
           y2=sp(i1+1)
           y3=sp(i1+2)
           y4=sp(i1+3)
 
           x0(1,i)=y1
           x0(2,i)=y2
           x0(3,i)=-y3
           x0(4,i)=-y4
 
           xx(1,i) = -y4
           xx(2,i) = -y3
           xx(3,i) = y2
           xx(4,i) = y1
 
           xy(1,i)=dcmplx(0d0,-1d0)*y4
           xy(2,i)=dcmplx(0d0,1d0)*y3
           xy(3,i)=dcmplx(0d0,1d0)*y2
           xy(4,i)=dcmplx(0d0,-1d0)*y1
 
           xz(1,i)=-y3
           xz(2,i)=y4
           xz(3,i)=y1
           xz(4,i)=-y2
 
           x5(1,i)=y3
           x5(2,i)=y4
           x5(3,i)=y1
           x5(4,i)=y2
 
           enddo
 
           if (dv.eq.4) then
 
           do i=1,4
 
           f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)
 
           enddo
 
           endif
 
           if (dv.eq.6) then
           bp = (v(5)+dcmplx(0d0,1d0)*v(6))
           bm=(v(5)-dcmplx(0d0,1d0)*v(6))
 
           do i=1,4
 
           f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)+bm*x5(i,2)
 
            i1 = i+4
 
            f(i1)= v(1)*x0(i,2)-v(2)*xx(i,2)-v(3)*xy(i,2)-v(4)*xz(i,2)-bp*x5(i,1)
 
 
            enddo
 
           endif
 
           if (dv.eq.8) then
           bp=(v(5)+dcmplx(0d0,1d0)*v(6))
           bm=(v(5)-dcmplx(0d0,1d0)*v(6))
           cp=(v(7)+dcmplx(0d0,1d0)*v(8))
           cm=(v(7)-dcmplx(0d0,1d0)*v(8))
 
           do i=1,4
 
           f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)+bm*x5(i,2)-cm*x5(i,3)
 
            i1 = i+4
 
            f(i1)= v(1)*x0(i,2)-v(2)*xx(i,2)-v(3)*xy(i,2)-v(4)*xz(i,2)-bp*x5(i,1)+cm*x5(i,4)
 
             i2 = i1+4
 
             f(i2)=v(1)*x0(i,3)-v(2)*xx(i,3)-v(3)*xy(i,3)-v(4)*xz(i,3)+bm*x5(i,4)+cp*x5(i,1)
 
              i3=i2+4
 
              f(i3)=v(1)*x0(i,4)-v(2)*xx(i,4)-v(3)*xy(i,4)-v(4)*xz(i,4)-bp*x5(i,3)-cp*x5(i,2)
 
              enddo
 
              endif
 
 
 
! if(Ds.eq.16) then
!   test(:) = SpiVL(sp,v)
!   print *, ""
!   print *, "spb2",f(1:Ds)
!   print *, "test",test(1:Ds)
!   print *, "diff",test(1:Ds)-f(1:Ds)
!   if( any( abs(test(1:Ds)-f(1:Ds) ).gt.1d-10  ) ) pause
! endif
 
 
 
 
               return
end subroutine
 
 
        function spb2_(sp,v)
        implicit none
        double complex, intent(in) :: sp(:),v(:)
        double complex :: spb2_(size(sp)) ,tmp(size(sp))
        integer :: dv,ds
 
          ds = size(sp)
          if (ds == 4) dv = 4
          if (ds == 8) dv = 6
          if (ds == 16) dv = 8
          call spb2(dv,ds,sp,v,spb2_)
 
!  tmp = SpiVL(sp(1:Ds),v(1:Dv))
!  print *, "diff2",tmp-spb2_
!  pause
        return
        end function
 
 
 
 
          function psp1_(sp1,sp2) result(res)
          implicit none
          complex(8), intent(in) :: sp1(:)
          complex(8), intent(in) :: sp2(:)
          complex(8) :: res
 
            res = sum(sp1(1:)*sp2(1:))
 
           end function
 
 
 
 
! RR -- new for D-dim chirality
             recursive function chir(sign,sp) result(res)
               implicit none
               logical :: sign
               double complex :: sp(:)
               double complex :: res(size(sp))
               integer        :: d
 
               d = size(sp)
               if ( d .eq. 4) then
                  if(sign) then !omega_+
                     res(1) = 0.5d0*(sp(1)+sp(3))
                     res(2) = 0.5d0*(sp(2)+sp(4))
                     res(3) = res(1)
                     res(4) = res(2)
                  else !omega_-
                     res(1) = 0.5d0*(sp(1)-sp(3))
                     res(2) = 0.5d0*(sp(2)-sp(4))
                     res(3) =-res(1)
                     res(4) =-res(2)
                  endif
               else
                  res(1:d/2)     =  chir(sign,sp(1:d/2))
                  res((d/2+1):d) =  chir(sign,sp( (d/2+1):d ))
               endif
 
             end function chir
 
 
 
 
             recursive function ichir(sign,sp) result(res)
! RR -- this function is needed for the electric and dipole moment like couplings,
!!      which have a 1 +\- i*gamma5 (see 0811.3842)
               implicit none
               logical :: sign
               double complex :: sp(:)
               double complex :: res(size(sp))
               double complex :: ci
               integer        :: d
 
               ci=(0d0,1d0)
               d = size(sp)
               if ( d .eq. 4) then
                  if(sign) then !omega_+
                     res(1) = 0.5d0*(sp(1)+ci*sp(3))
                     res(2) = 0.5d0*(sp(2)+ci*sp(4))
                     res(3) = 0.5d0*(ci*sp(1)+sp(3))
                     res(4) = 0.5d0*(ci*sp(2)+sp(4))
                  else !omega_-
                     res(1) = 0.5d0*(sp(1)-ci*sp(3))
                     res(2) = 0.5d0*(sp(2)-ci*sp(4))
                     res(3) = 0.5d0*(-ci*sp(1)+sp(3))
                     res(4) = 0.5d0*(-ci*sp(2)+sp(4))
                  endif
               else
                  res(1:d/2)     =  ichir(sign,sp(1:d/2))
                  res((d/2+1):d) =  ichir(sign,sp( (d/2+1):d ))
               endif
 
             end function ichir
 
 
 
 
      function vbqv(sp,e1,coupl_left,coupl_right)
      implicit none
      complex(8), intent(in) :: e1(:)
      complex(8), intent(in) :: sp(:)
      complex(8), intent(in) :: coupl_left,coupl_right
      complex(8) :: vbqv(size(sp))
 
!            vbqV = -(0d0,1d0)*( coupl_left*Chir(.false.,spb2_(sp,e1)) + coupl_right*Chir(.true.,spb2_(sp,e1)) )
           vbqv = -(0d0,1d0)*( coupl_left*chir(.false.,sp) + coupl_right*chir(.true.,sp) )
 
      return
      end function
 
 
 
      function vvq(e1,sp,coupl_left,coupl_right)
      implicit none
      complex(8), intent(in) :: e1(:)
      complex(8), intent(in) :: sp(:)
      complex(8), intent(in) :: coupl_left,coupl_right
      complex(8) :: vvq(size(sp))
 
!            vVq = -(0d0,1d0)*( coupl_left*Chir(.true., spi2_(e1,sp)) + coupl_right*Chir(.false., spi2_(e1,sp)) )
           vvq = -(0d0,1d0)*( coupl_left*chir(.false.,sp) + coupl_right*chir(.true.,sp) )
 
 
      return
      end function
 
 
 
 
          function sc_(p1,p2)
          implicit none
          complex(8) :: p1(:),p2(:)
          complex(8) :: sc_
          integer :: sizemin
 
              sizemin=min(size(p1),size(p2))
              call rsc_(sizemin,p1,p2,sc_)
 
          return
          end function
 
 
 
         subroutine rsc_(n,x,y,r)
         implicit none
         integer i,n
         complex(8) x(*),y(*)
         complex(8) r
 
            r = x(1)*y(1)
            do i=2, n
              r = r - x(i)*y(i)
            enddo
 
         return
         end subroutine
 
 
 
 
 
 
      function vggg(e1,k1,e2,k2)
      implicit none
      complex(8), intent(in) :: e1(:), e2(:)
      complex(8), intent(in) :: k1(:), k2(:)
      complex(8)             :: vggg(size(e1))
      complex(8):: sk1e2,se1e2,sk2e1,xx
      real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
 
 
          sk1e2=sc_(k1,e2)
          sk2e1=sc_(k2,e1)
          se1e2=sc_(e1,e2)
          xx=(0.0d0,1.0d0)*sqrt2
          vggg = xx*(-sk1e2*e1+sk2e1*e2+se1e2/2d0*(k1-k2))
 
       end function vggg
 
 
       function  vgggg(e1,e2,e3)
       implicit none
       complex(8), intent(in) :: e1(:),e2(:),e3(:)
       complex(8)             :: vgggg(size(e1))
       complex(8):: se1e3,se2e3,se1e2
 
          se1e3=sc_(e1,e3)
          se2e3=sc_(e2,e3)
          se1e2=sc_(e1,e2)
          vgggg = (0.0d0,1.0d0)*(e2*se1e3-0.5d0*(e1*se2e3+ e3*se1e2))
 
       end function vgggg
 
 
      function vqg(sp,e1)
      implicit none
      complex(8), intent(in) :: e1(:)
      complex(8), intent(in) :: sp(:)
      complex(8) :: vqg(size(sp))
      real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
 
          vqg = (0d0,1d0)/sqrt2*spb2_(sp,e1)
 
      end function vqg
 
 
 
 
      function vgq(e1,sp)
      implicit none
      complex(8), intent(in) :: e1(:)
      complex(8), intent(in) :: sp(:)
      complex(8) :: vgq(size(sp))
      real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
 
          vgq = (0d0,-1d0)/sqrt2*spb2_(sp,e1)
 
      end function vgq
 
 
 
 
 
      function vbqg(sp,e1)
      implicit none
      complex(8), intent(in) :: e1(:)
      complex(8), intent(in) :: sp(:)
      complex(8) :: vbqg(size(sp))
      real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
 
          vbqg = (0d0,-1d0)/sqrt2*spi2_(e1,sp)
 
      end function vbqg
 
 
 
      function vgbq(e1,sp)
      implicit none
      complex(8), intent(in) :: e1(:)
      complex(8), intent(in) :: sp(:)
      complex(8) :: vgbq(size(sp))
      real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
 
          vgbq = (0d0,1d0)/sqrt2*spi2_(e1,sp)
 
      end function vgbq
 
 
 
!       function vbqq(Dv,sp1,sp2)
!       implicit none
!       complex(8), intent(in) :: sp1(:), sp2(:)
!       integer, intent(in) ::  Dv
!       integer :: i
!       complex(8) :: vbqq(Dv)
!       complex(8) :: rr, va(Dv),sp1a(size(sp1))
!       real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
!
!           va=(0d0,0d0)
!           vbqq=(0d0,0d0)
!
!           do i=1,Dv
!              if (i.eq.1) then
!                va(1)=(1d0,0d0)
!              else
!                va(i)=(-1d0,0d0)
!              endif
!              sp1a=spb2_(sp1,va)
!
!              rr=(0d0,-1d0)/sqrt2*psp1_(sp1a,sp2)
!              if (i.eq.1) then
!                   vbqq = vbqq + rr*va
!               else
!                   vbqq = vbqq - rr*va
!              endif
!              va(i)=(0d0,0d0)
!           enddo
!
!       end function vbqq
 
 
 
 
 
      function vbqq(Dummy,sp1,sp2)! this is my own simpler 4-dim version
      implicit none
      complex(8), intent(in) :: sp1(1:4), sp2(1:4)
      integer :: i,dummy
      complex(8) :: vbqq(4)
      complex(8) :: sp1a(4)
      real(8) :: va(1:4,1:4)
      real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
 
         va(1,1:4)=(/+1d0,0d0,0d0,0d0/)
         va(2,1:4)=(/0d0,-1d0,0d0,0d0/)
         va(3,1:4)=(/0d0,0d0,-1d0,0d0/)
         va(4,1:4)=(/0d0,0d0,0d0,-1d0/)
 
          do i=1,4
             sp1a=spivl(sp1,dcmplx(va(i,1:4)))
             vbqq(i) = (sp1a(1)*sp2(1)+sp1a(2)*sp2(2)+sp1a(3)*sp2(3)+sp1a(4)*sp2(4)) * (0d0,-1d0)/sqrt2
          enddo
 
 
      end function vbqq
 
 
 
 
      function vvbqq(sp1,sp2)
      implicit none
      complex(8), intent(in) :: sp1(:), sp2(:)
      integer :: i,j
      complex(8) :: vvbqq(4,4)
      complex(8) :: sp1a(4)
      real(8) :: va(1:4,1:4)
 
         va(1,1:4)=(/+1d0,0d0,0d0,0d0/)
         va(2,1:4)=(/0d0,-1d0,0d0,0d0/)
         va(3,1:4)=(/0d0,0d0,-1d0,0d0/)
         va(4,1:4)=(/0d0,0d0,0d0,-1d0/)
 
         do i=1,4
          do j=1,4
             sp1a=spb2_(sp1, dcmplx(va(i,1:4)))
             sp1a=spb2_(sp1a,dcmplx(va(j,1:4)))
             vvbqq(i,j) = psp1_(sp1a,sp2)
          enddo
         enddo
 
      end function vvbqq
 
 
 
 
 
 
 
END MODULE