20 complex(8) :: mom(1:4)
21 complex(8) :: pol(1:4)
25 integer,
pointer :: parttype
26 integer,
pointer :: extref
27 real(8),
pointer :: mass
28 real(8),
pointer :: mass2
29 integer ,
pointer :: helicity
30 complex(8),
pointer :: mom(:)
31 complex(8),
pointer :: pol(:)
42 integer :: bosonvertex
43 integer,
allocatable :: numglu(:)
44 integer,
allocatable :: partref(:)
45 integer,
allocatable :: parttype(:)
53 INTERFACE OPERATOR (.ndot.)
55 END INTERFACE OPERATOR (.Ndot.)
59 integer,
parameter ::
dv=4,
ds=4
76 real(8),
intent(in) :: mom(1:4,1:13)
77 integer,
intent(in) :: selectprocess
78 real(8),
intent(out) :: res(-5:5,-5:5)
79 real(8) :: matelsq_gg,matelsq_qqb,matelsq_qbq
87 if( selectprocess.eq.0 )
then
89 elseif( selectprocess.eq.1 )
then
91 matelsq_qbq = matelsq_qqb
95 matelsq_qbq = matelsq_qqb
99 res(iq,iq) = matelsq_gg
101 res(iq,-iq) = matelsq_qqb
102 res(-iq,iq) = matelsq_qbq
112 real(8),
intent(in) :: mom(1:4,1:13)
113 integer,
intent(in) :: selectprocess
114 real(8),
intent(out) :: res(-5:5,-5:5)
115 real(8) :: tmptopmass
133 integer :: numquarks,numgluons,numboson
135 integer :: itree,numparticles
139 numquarks=2; numgluons=2; numboson=1;
140 numparticles=numquarks+numgluons+numboson
155 numquarks=4; numgluons=0; numboson=1;
156 numparticles=numquarks+numgluons+numboson
233 integer :: numtrees, itree
294 real(8) :: mom(1:4,1:13),sqamp
295 complex(8) :: resoffsh(1:4,1:2),res(1:2,1:2)
296 complex(8) :: glupol(1:4,1:2,1:2)
297 integer :: hel4,tophel1,tophel2,nhel
298 real(8),
parameter :: c_aa=64.d0/3.d0, c_ab=-8.d0/3.d0
299 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
337 res(1,1) = (resoffsh(1:4,1).ndot.glupol(1:4,1,1))
338 res(2,1) = (resoffsh(1:4,2).ndot.glupol(1:4,1,1))
339 res(1,2) = (resoffsh(1:4,1).ndot.glupol(1:4,1,2))
340 res(2,2) = (resoffsh(1:4,2).ndot.glupol(1:4,1,2))
343 + c_aa * dreal( res(1,1)*dconjg(res(1,1)) + res(1,2)*dconjg(res(1,2)) ) &
344 + c_ab * dreal( res(1,1)*dconjg(res(2,1)) + res(1,2)*dconjg(res(2,2)) ) &
345 + c_ab * dreal( res(2,1)*dconjg(res(1,1)) + res(2,2)*dconjg(res(1,2)) ) &
346 + c_aa * dreal( res(2,1)*dconjg(res(2,1)) + res(2,2)*dconjg(res(2,2)) )
366 real(8) :: mom(1:4,1:13),sqamp
367 complex(8) :: resoffsh(1:4),res(1:2)
368 complex(8) :: quapol(1:4,1:2,1:2)
369 integer :: hel4,tophel1,tophel2,nhel
370 real(8),
parameter :: c_aa=8.0d0
371 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
407 res(1) =
psp1_(resoffsh(1:4),quapol(1:4,1,1))
408 res(2) =
psp1_(resoffsh(1:4),quapol(1:4,1,2))
411 + c_aa * dreal( res(1)*dconjg(res(1)) + res(2)*dconjg(res(2)) )
429 integer :: tag_f,tag_Z,n
430 complex(8) :: Res(1:Ds)
432 logical,
parameter :: Boson=.true.
433 integer :: i,j,Order(1:6)
436 if( treeproc%NumQua.eq.2 .and. treeproc%NumSca.eq.0 ) then
437 if ( treeproc%PartType(1).eq.
glu_ .and. boson )
then
438 res(1:
dv) =
cur_g_2fv( treeproc%Gluons(2:treeproc%NumGlu(0)),treeproc%Quarks(1:treeproc%NumQua),treeproc%Boson,treeproc%NumGlu(0:3) )
439 elseif(
isaquark(treeproc%PartType(1)) .and. boson )
then
440 if( treeproc%NumV.eq.1 )
then
441 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))
443 print *,
"requested current with a boson is not available"
447 print *,
"requested current is not available 2q"
451 elseif( treeproc%NumQua.eq.4 .and. treeproc%NumSca.eq.0) then
452 if( treeproc%PartType(1).eq.
glu_ )
then
453 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) )
454 elseif(
isaquark(treeproc%PartType(1)) )
then
455 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 )
457 print *,
"requested current is not available 4q"
461 print *,
"requested current is not available xx"
474 type(
particle),
target :: TheParticles(:)
475 integer :: iPart,PartRef,PartType,ig,iq,ib,NPart,counterQ,counterG,LastQuark,QuarkPos(1:6)
483 do npart=1,thetreeamp%NumPart
484 thetreeamp%PartType(npart) = theparticles( thetreeamp%PartRef(npart) )%PartType
485 if(
isaquark(thetreeamp%PartType(npart)) )
then
487 counterq = counterq + 1
488 quarkpos(counterq) = counterq + counterg
494 elseif( thetreeamp%PartType(npart).eq.
glu_ )
then
495 counterg = counterg + 1
496 elseif(
isaboson(thetreeamp%PartType(npart)) ) then
497 if( npart.eq.1 ) print *,
"Vector boson should not be the first particle."
498 if( abs(thetreeamp%PartType(npart)).eq.abs(
wp_) ) thetreeamp%NumW = thetreeamp%NumW + 1
499 if( abs(thetreeamp%PartType(npart)).eq.abs(
z0_) ) thetreeamp%NumV = thetreeamp%NumV + 1
500 if( abs(thetreeamp%PartType(npart)).eq.abs(
pho_)) thetreeamp%NumV = thetreeamp%NumV + 1
502 thetreeamp%BosonVertex = lastquark
509 if(
isaquark( thetreeamp%PartType(1) ) )
then
510 if( thetreeamp%NumQua+thetreeamp%NumSca .eq. 2 )
then
511 thetreeamp%NumGlu(1) = quarkpos(2) - quarkpos(1) - 1
512 thetreeamp%NumGlu(2) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(2)
514 if( thetreeamp%NumQua .eq. 4 )
then
515 thetreeamp%NumGlu(1) = quarkpos(2) - quarkpos(1) - 1
516 thetreeamp%NumGlu(2) = quarkpos(3) - quarkpos(2) - 1
517 thetreeamp%NumGlu(3) = quarkpos(4) - quarkpos(3) - 1
518 thetreeamp%NumGlu(4) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(4)
520 if( thetreeamp%NumQua .eq. 6 )
then
521 thetreeamp%NumGlu(1) = quarkpos(2) - quarkpos(1) - 1
522 thetreeamp%NumGlu(2) = quarkpos(3) - quarkpos(2) - 1
523 thetreeamp%NumGlu(3) = quarkpos(4) - quarkpos(3) - 1
524 thetreeamp%NumGlu(4) = quarkpos(5) - quarkpos(4) - 1
525 thetreeamp%NumGlu(5) = quarkpos(6) - quarkpos(5) - 1
526 thetreeamp%NumGlu(6) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(6)
529 elseif( thetreeamp%PartType(1).eq.10 )
then
530 if( thetreeamp%NumQua .eq. 2 )
then
531 thetreeamp%NumGlu(1) = quarkpos(1) - 2
532 thetreeamp%NumGlu(2) = quarkpos(2) - quarkpos(1) - 1
533 thetreeamp%NumGlu(3) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(2)
535 if( thetreeamp%NumQua .eq. 4 )
then
536 thetreeamp%NumGlu(1) = quarkpos(1) - 2
537 thetreeamp%NumGlu(2) = quarkpos(2) - quarkpos(1) - 1
538 thetreeamp%NumGlu(3) = quarkpos(3) - quarkpos(2) - 1
539 thetreeamp%NumGlu(4) = quarkpos(4) - quarkpos(3) - 1
540 thetreeamp%NumGlu(5) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(4)
542 if( thetreeamp%NumQua .eq. 6 )
then
543 thetreeamp%NumGlu(1) = quarkpos(1) - 2
544 thetreeamp%NumGlu(2) = quarkpos(2) - quarkpos(1) - 1
545 thetreeamp%NumGlu(3) = quarkpos(3) - quarkpos(2) - 1
546 thetreeamp%NumGlu(4) = quarkpos(4) - quarkpos(3) - 1
547 thetreeamp%NumGlu(5) = quarkpos(5) - quarkpos(4) - 1
548 thetreeamp%NumGlu(6) = quarkpos(6) - quarkpos(5) - 1
549 thetreeamp%NumGlu(7) = thetreeamp%NumSca+thetreeamp%NumQua+thetreeamp%NumGlu(0) - quarkpos(6)
555 do ipart=1,thetreeamp%NumPart
556 partref = thetreeamp%PartRef(ipart)
557 parttype= theparticles(partref)%PartType
558 if( parttype.eq.
glu_ )
then
560 thetreeamp%Gluons(ig)%PartType => theparticles(partref)%PartType
561 thetreeamp%Gluons(ig)%ExtRef => theparticles(partref)%ExtRef
562 thetreeamp%Gluons(ig)%Mass => theparticles(partref)%Mass
563 thetreeamp%Gluons(ig)%Mass2 => theparticles(partref)%Mass2
564 thetreeamp%Gluons(ig)%Helicity => theparticles(partref)%Helicity
565 thetreeamp%Gluons(ig)%Mom => theparticles(partref)%Mom
566 thetreeamp%Gluons(ig)%Pol => theparticles(partref)%Pol
567 if( parttype.ne.theparticles(partref)%PartType ) print *,
"Error1 in LinkTreeParticles"
570 thetreeamp%Quarks(iq)%PartType => theparticles(partref)%PartType
571 thetreeamp%Quarks(iq)%ExtRef => theparticles(partref)%ExtRef
572 thetreeamp%Quarks(iq)%Mass => theparticles(partref)%Mass
573 thetreeamp%Quarks(iq)%Mass2 => theparticles(partref)%Mass2
574 thetreeamp%Quarks(iq)%Helicity => theparticles(partref)%Helicity
575 thetreeamp%Quarks(iq)%Mom => theparticles(partref)%Mom
576 thetreeamp%Quarks(iq)%Pol => theparticles(partref)%Pol
577 if( parttype.ne.theparticles(partref)%PartType ) print *,
"Error2 in LinkTreeParticles"
580 if( ib.ge.2 ) print *,
"Too many bosons in LinkTreeParticles"
581 thetreeamp%Boson%PartType => theparticles(partref)%PartType
582 thetreeamp%Boson%ExtRef => theparticles(partref)%ExtRef
583 thetreeamp%Boson%Mass => theparticles(partref)%Mass
584 thetreeamp%Boson%Mass2 => theparticles(partref)%Mass2
585 thetreeamp%Boson%Helicity => theparticles(partref)%Helicity
586 thetreeamp%Boson%Mom => theparticles(partref)%Mom
587 thetreeamp%Boson%Pol => theparticles(partref)%Pol
588 if( parttype.ne.theparticles(partref)%PartType ) print *,
"Error2 in LinkTreeParticles"
591 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"
648 FUNCTION cur_f_2fv(Gluons,Quark,Quark1PartType,Boson,NumGlu)
result(Res)
650 complex(8) :: res(1:
ds)
651 integer :: numglu(0:2),i,rin,rout,quark1parttype
653 complex(8) :: glumom(1:
dv,numglu(0)), quarkmom(1:
dv)
654 complex(8) :: glupol(1:
dv,numglu(0)), quarkpol(1:
ds)
659 glumom(1:
dv,i) = gluons(i)%Mom(1:
dv)
660 glupol(1:
dv,i) = gluons(i)%Pol(1:
dv)
662 quarkmom(1:
dv) = quark(2)%Mom(1:
dv)
663 quarkpol(1:
ds) = quark(2)%Pol(1:
ds)
665 if( quark1parttype.ne.-quark(2)%PartType ) print *,
"Wrong quark flavors in cur_f_2fV"
666 if( numglu(0)-numglu(1)-numglu(2).ne.0 ) print *,
"Wrong NumGlu in cur_f_2fV",numglu(0)-numglu(1)-numglu(2)
670 if( quark(2)%PartType .gt.0 )
then
671 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))
673 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))
686 recursive function fv(e,k,sp,p,mass,QuarkFlavor,eV,kV,ms)
result(res)
688 complex(8),
intent(in) :: e(:,:), k(:,:)
689 complex(8),
intent(in) :: sp(:), p(:)
690 complex(8),
intent(in) ::
ev(:), kv(:)
691 integer,
intent(in) :: ms,quarkflavor
692 integer :: ms1,m,ng1, ng2
694 complex(8) :: res(size(sp))
695 complex(8) :: tmp(size(sp))
696 complex(8) :: k1(size(p))
697 complex(8) :: k2(size(p))
698 complex(8) :: sp2(size(sp))
699 complex(8) :: sp3(size(sp))
700 complex(8) :: e1(size(e,dim=1))
701 complex(8) :: e2(size(e,dim=1))
702 complex(8) :: k1sq,k2sq,k3sq
704 complex(8) :: couplvqq_left,couplvqq_right,couplvqq_left2,couplvqq_right2
705 character,
parameter :: ferfla*3=
"dum"
708 ngluon =
size(e,dim=2)
711 if (ng2 < 0)
write(*,*)
'WRONG DEFINITION OF CURRENT fV'
713 if( abs(quarkflavor).eq.
top_ .or. abs(quarkflavor).eq.
bot_ ) then
719 print *,
"this should not happen for the Higgs",quarkflavor,
top_
724 if (ngluon == 0)
then
725 res =
vbqv(sp,
ev,couplvqq_left,couplvqq_right)
730 k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
731 e1=
g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
734 k2 = sum(k(:,1:ng1+m),dim=2)
736 k2sq =
sc_(k2,k2)-mass**2
737 sp2 =
fv(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,quarkflavor,
ev,kv,ng1)
738 sp2 =
spb2_(sp2,k2)+mass*sp2
743 tmp = -(0d0,1d0)/k1sq*tmp
750 tmp = (0d0,1d0)/k2sq*tmp
761 k1 = sum(k(:,1:m),dim=2)
762 e1=
g(e(:,1:m),k(:,1:m))
765 k2 = sum(k(:,m+1:ngluon),dim=2)
767 k2sq =
sc_(k2,k2) - mass**2
769 sp2=
fv(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,quarkflavor,
ev,kv,ms1)
771 sp2 =
spb2_(sp2,k2)+mass*sp2
776 tmp=-(0d0,1d0)/k1sq*tmp
783 tmp=(0d0,1d0)/k2sq*tmp
792 sp2 =
f(e,k,sp,p,mass,ferfla,ferfla,ms)
793 k2 = sum(k(:,1:ngluon),dim=2)
795 k2sq =
sc_(k2,k2) - mass**2
797 sp2 =
spb2_(sp2,k2)+ mass*sp2
799 tmp =
vbqv(sp2,
ev,couplvqq_left,couplvqq_right)
801 tmp = (0d0,1d0)/k2sq*tmp
814 recursive function bfv(e,k,sp,p,mass,QuarkFlavor,eV,kV,ms)
result(res)
816 complex(8),
intent(in) :: e(:,:), k(:,:)
817 complex(8),
intent(in) :: sp(:), p(:)
818 complex(8),
intent(in) ::
ev(:), kv(:)
819 integer,
intent(in) :: ms,quarkflavor
820 integer :: ms1,m,ng1, ng2
822 complex(8) :: res(size(sp))
823 complex(8) :: tmp(size(sp))
824 complex(8) :: k1(size(p))
825 complex(8) :: k2(size(p))
826 complex(8) :: sp2(size(sp))
827 complex(8) :: sp3(size(sp))
828 complex(8) :: e1(size(e,dim=1))
829 complex(8) :: e2(size(e,dim=1))
830 complex(8) :: k1sq,k2sq,k3sq
832 complex(8) :: couplvqq_left,couplvqq_right,couplvqq_left2,couplvqq_right2
833 character,
parameter :: ferfla*3=
"dum"
836 ngluon =
size(e,dim=2)
840 if (ng2 < 0)
write(*,*)
'WRONG DEFINITION OF CURRENT fbV'
842 if( abs(quarkflavor).eq.
top_ .or. abs(quarkflavor).eq.
bot_ ) then
848 print *,
"this should not happen for the Higgs",quarkflavor,
top_
853 if (ngluon == 0)
then
854 res =
vvq(
ev,sp,couplvqq_left,couplvqq_right)
859 k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
860 e1=
g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
863 k2 = sum(k(:,1:ng1+m),dim=2)
865 k2sq =
sc_(k2,k2)-mass**2
867 sp2 =
bfv(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,quarkflavor,
ev,kv,ng1)
868 sp2 =
spi2_(k2,sp2)+mass*sp2
874 tmp = -(0d0,1d0)/k1sq*tmp
881 tmp = (0d0,1d0)/k2sq*tmp
892 k1 = sum(k(:,1:m),dim=2)
893 e1=
g(e(:,1:m),k(:,1:m))
896 k2 = sum(k(:,m+1:ngluon),dim=2)
898 k2sq =
sc_(k2,k2) - mass**2
900 sp2=
bfv(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,quarkflavor,
ev,kv,ms1)
902 sp2 =
spi2_(k2,sp2)+mass*sp2
908 tmp=-(0d0,1d0)/k1sq*tmp
915 tmp=(0d0,1d0)/k2sq*tmp
926 sp2 =
bf(e,k,sp,p,mass,ferfla,ferfla,ms)
927 k2 = sum(k(:,1:ngluon),dim=2)
929 k2sq =
sc_(k2,k2) -mass**2
932 sp2 =
spi2_(k2,sp2) +mass*sp2
934 tmp =
vvq(
ev,sp2,couplvqq_left,couplvqq_right)
936 tmp = (0d0,1d0)/k2sq*tmp
953 FUNCTION cur_f_4fv(Gluons,Quarks,Quark1PartType,Boson,BosonVertex,NumGlu,tag_f,tag_Z)
result(res)
955 integer :: numglu(0:4),quark1parttype
957 integer :: tag_f,bosonvertex,tag_z
958 integer,
target :: tmpextref
959 complex(8) :: res(1:
ds),tmp(1:
ds)
960 complex(8) :: ubar1(1:
ds)
961 complex(8),
target ::
ubar0(1:
ds)
962 complex(8) :: eps1(1:
dv)
963 complex(8) :: eps2(1:
dv)
964 type(
ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(4)),tmpquark(1:1)
965 complex(8) :: propfac1,propfac2
966 complex(8),
target :: pmom1(1:
dv)
967 complex(8) :: pmom2(1:
dv)
968 integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b
969 integer :: rin,rout,i,counter
973 if( numglu(0)-numglu(1)-numglu(2)-numglu(3)-numglu(4).ne.0 ) print *,
"wrong number of gluons in cur_f_4f"
974 if(quarks(3)%PartType.eq.-quarks(4)%PartType .and. quark1parttype.ne.-quarks(2)%PartType ) print *,
"wrong flavor in cur_f_4f (1)"
975 if(quarks(2)%PartType.eq.-quarks(3)%PartType .and. quark1parttype.ne.-quarks(4)%PartType ) print *,
"wrong flavor in cur_f_4f (2)"
978 if ( quark1parttype.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType )
then
984 rout=numglu(1)+numglu(2)+numglu(3)+n4a
985 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 2 .or. bosonvertex .eq. 4)
then
986 eps2 =
cur_g_2f(gluons(rin:rout),quarks(3:4),(/1+n2b+numglu(3)+n4a,n2b,numglu(3),n4a/))
987 pmom1(:) =
summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom
988 elseif (bosonvertex .eq. 3)
then
989 eps2 =
cur_g_2fv(gluons(rin:rout),quarks(3:4),boson,(/1+n2b+numglu(3)+n4a,n2b,numglu(3),n4a/))
990 pmom1(:) =
summom(gluons,rin,rout) + quarks(3)%Mom(:) + quarks(4)%Mom(:) + boson%Mom(:)
992 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
998 if (bosonvertex.eq.1 .or. bosonvertex.eq.2)
then
1001 ubar1(:) =
cur_f_2fv(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,boson,(/n2a+n1b,n1b,n2a/) )
1002 pmom2(:) = quarks(2)%Mom(:) +
summom(gluons,rin,rout) + boson%Mom(:)
1003 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
1005 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut )
then
1009 if( quarks(2)%PartType.lt.0 )
then
1010 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
1012 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
1015 if( quarks(2)%PartType.lt.0 )
then
1021 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a) + boson%Mom(:)
1022 if(n1a.ge.1 .or. n4b.ge.1)
then
1023 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
1024 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut )
then
1027 if( quarks(2)%PartType.lt.0 )
then
1034 tmpquark(1)%Mom => pmom1(:)
1035 tmpquark(1)%Pol =>
ubar0(:)
1036 tmpquark(1)%Mass => quarks(2)%Mass
1037 tmpquark(1)%Mass2=> quarks(2)%Mass2
1039 tmpquark(1)%ExtRef => tmpextref
1040 tmpquark(1)%PartType => quarks(2)%PartType
1048 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1054 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
1055 res(:) = res(:) + tmp(:)
1059 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 4)
then
1062 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2a+n1b,n1b,n2a/) )
1063 if(n1b.ge.1 .or. n2a.ge.1)
then
1064 pmom2(:) = quarks(2)%Mom +
summom(gluons,rin,rout)
1065 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
1066 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut ) cycle
1067 if( quarks(2)%PartType.lt.0 )
then
1068 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
1070 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
1073 if( quarks(2)%PartType.lt.0 )
then
1079 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
1080 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
1081 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut )
then
1084 if( quarks(2)%PartType.lt.0 )
then
1090 tmpquark(1)%Mom => pmom1(:)
1091 tmpquark(1)%Pol =>
ubar0(:)
1092 tmpquark(1)%Mass => quarks(2)%Mass
1093 tmpquark(1)%Mass2=> quarks(2)%Mass2
1095 tmpquark(1)%ExtRef => tmpextref
1096 tmpquark(1)%PartType => quarks(2)%PartType
1104 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1110 tmp(:) =
cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,boson,(/counter-1,n1a,n4b/) )
1111 res(:) = res(:) + tmp(:)
1114 if (bosonvertex .eq. 3)
then
1117 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2a+n1b,n1b,n2a/) )
1118 if(n1b.ge.1 .or. n2a.ge.1)
then
1119 pmom2(:) = quarks(2)%Mom +
summom(gluons,rin,rout)
1120 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
1121 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut )
then
1124 if( quarks(2)%PartType.lt.0 )
then
1125 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
1127 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
1130 if( quarks(2)%PartType.lt.0 )
then
1136 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)+boson%Mom
1137 if(n1a.ge.1 .or. n4b.ge.1)
then
1138 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
1139 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut ) cycle
1140 if( quarks(2)%PartType.lt.0 )
then
1147 tmpquark(1)%Mom => pmom1(:)
1148 tmpquark(1)%Pol =>
ubar0(:)
1149 tmpquark(1)%Mass => quarks(2)%Mass
1150 tmpquark(1)%Mass2=> quarks(2)%Mass2
1152 tmpquark(1)%ExtRef => tmpextref
1153 tmpquark(1)%PartType => quarks(2)%PartType
1161 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1167 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
1168 res(:) = res(:) + tmp(:)
1176 if( (quark1parttype.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType) .AND. &
1177 (quarks(4)%ExtRef.ne.-1 .or. tag_f.ne.1 .or. abs(quark1parttype).ne.abs(quarks(2)%PartType)) )
then
1185 rout=numglu(1)+numglu(2)+n3a
1190 if (bosonvertex.eq.1 .or. bosonvertex.eq.3 .or. bosonvertex.eq.4)
then
1191 if ( tag_z .eq. 1 .and. n1b+numglu(2)+n3a .eq. numglu(0) )
then
1195 eps2 =
cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
1196 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
1197 elseif (bosonvertex .eq. 2)
then
1198 eps2 =
cur_g_2fv(gluons(rin:rout),quarks(2:3),boson,(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
1199 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
1201 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
1203 eps2 = eps2*propfac1
1208 if (bosonvertex .eq. 3 .or. bosonvertex .eq. 4)
then
1209 rin =numglu(1)+numglu(2)+n3a+1
1210 rout=numglu(1)+numglu(2)+numglu(3)+n4a
1212 ubar1(:) =
cur_f_2fv(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,boson,(/n3b+n4a,n3b,n4a/) )
1213 pmom2(:) = quarks(4)%Mom(:) +
summom(gluons,rin,rout) + boson%Mom(:)
1214 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(4)%Mass2)
1215 if( abs(
sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.
propcut )
then
1219 if( quarks(4)%PartType.lt.0 )
then
1220 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
1222 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
1224 if( quarks(4)%PartType.lt.0 )
then
1230 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a) + boson%Mom(:)
1231 if(n1a.ge.1 .or. n4b.ge.1)
then
1232 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(4)%Mass2)
1233 if( abs(
sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.
propcut )
then
1236 if( quarks(4)%PartType.lt.0 )
then
1243 tmpquark(1)%Mom => pmom1(:)
1244 tmpquark(1)%Pol =>
ubar0(:)
1245 tmpquark(1)%Mass => quarks(4)%Mass
1246 tmpquark(1)%Mass2=> quarks(4)%Mass2
1248 tmpquark(1)%ExtRef => tmpextref
1249 tmpquark(1)%PartType => quarks(4)%PartType
1257 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1263 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
1264 res(:) = res(:) + tmp(:)
1267 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 4)
then
1269 rin =numglu(1)+numglu(2)+n3a+1
1270 rout=numglu(1)+numglu(2)+numglu(3)+n4a
1271 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/n3b+n4a,n3b,n4a/) )
1272 if(n3b.ge.1 .or. n4a.ge.1)
then
1273 pmom2(:) = quarks(4)%Mom +
summom(gluons,rin,rout)
1274 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(4)%Mass2)
1275 if( abs(
sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.
propcut )
then
1279 if( quarks(4)%PartType.lt.0 )
then
1280 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
1282 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
1285 if( quarks(4)%PartType.lt.0 )
then
1291 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
1292 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(4)%Mass2)
1293 if( abs(
sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.
propcut )
then
1296 if( quarks(4)%PartType.lt.0 )
then
1302 tmpquark(1)%Mom => pmom1(:)
1303 tmpquark(1)%Pol =>
ubar0(:)
1304 tmpquark(1)%Mass => quarks(4)%Mass
1305 tmpquark(1)%Mass2=> quarks(4)%Mass2
1307 tmpquark(1)%ExtRef => tmpextref
1308 tmpquark(1)%PartType => quarks(4)%PartType
1316 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1322 tmp(:) =
cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,boson,(/counter-1,n1a,n4b/) )
1323 res(:) = res(:) + tmp(:)
1327 if (bosonvertex .eq. 2)
then
1328 rin =numglu(1)+numglu(2)+n3a+1
1329 rout=numglu(1)+numglu(2)+numglu(3)+n4a
1330 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/n3b+n4a,n3b,n4a/) )
1331 if ( n3b .ge. 1 .or. n4a .ge. 1)
then
1332 pmom2(:) = quarks(4)%Mom(:) +
summom(gluons,rin,rout)
1333 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(4)%Mass2)
1335 if( abs(
sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.
propcut )
then
1339 if( quarks(4)%PartType.lt.0 )
then
1340 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
1342 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
1345 if( quarks(4)%PartType.lt.0 )
then
1350 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)+boson%Mom
1351 if(n1a.ge.1 .or. n4b.ge.1)
then
1352 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(4)%Mass2)
1353 if( abs(
sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.
propcut )
then
1356 if( quarks(4)%PartType.lt.0 )
then
1363 tmpquark(1)%Mom => pmom1(:)
1364 tmpquark(1)%Pol =>
ubar0(:)
1365 tmpquark(1)%Mass => quarks(4)%Mass
1366 tmpquark(1)%Mass2=> quarks(4)%Mass2
1368 tmpquark(1)%ExtRef => tmpextref
1369 tmpquark(1)%PartType => quarks(4)%PartType
1377 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1383 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
1384 res(:) = res(:) + tmp(:)
1400 FUNCTION cur_f_4f(Gluons,Quarks,Quark1PartType,NumGlu,tag_f,tag_Z_arg)
result(res)
1402 integer :: numglu(0:4),quark1parttype
1404 integer :: tag_f,tag_z
1405 integer,
optional :: tag_z_arg
1406 integer,
target :: tmpextref
1407 complex(8) :: res(1:
ds),tmp(1:
ds)
1408 complex(8) :: ubar1(1:
ds)
1409 complex(8),
target ::
ubar0(1:
ds)
1410 complex(8) :: eps1(1:
dv)
1411 complex(8) :: eps2(1:
dv)
1412 type(
ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(4)),tmpquark(1:1)
1413 complex(8) :: propfac1,propfac2
1414 complex(8),
target :: pmom1(1:
dv)
1415 complex(8) :: pmom2(1:
dv)
1416 integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b
1417 integer :: rin,rout,i,counter
1421 if( numglu(0)-numglu(1)-numglu(2)-numglu(3)-numglu(4).ne.0 ) print *,
"wrong number of gluons in cur_f_4f"
1422 if(quarks(3)%PartType.eq.-quarks(4)%PartType .and. quark1parttype.ne.-quarks(2)%PartType ) print *,
"wrong flavor in cur_f_4f (1)"
1423 if(quarks(2)%PartType.eq.-quarks(3)%PartType .and. quark1parttype.ne.-quarks(4)%PartType ) print *,
"wrong flavor in cur_f_4f (2)"
1429 if(
present(tag_z_arg) )
then
1433 if( quark1parttype.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType )
then
1440 rin =numglu(1)+n2a+1
1441 rout=numglu(1)+numglu(2)+numglu(3)+n4a
1442 eps2 =
cur_g_2f(gluons(rin:rout),quarks(3:4),(/1+n2b+numglu(3)+n4a,n2b,numglu(3),n4a/))
1443 pmom1(:) =
summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom
1444 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
1446 eps2 = eps2*propfac1
1452 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2a+n1b,n1b,n2a/) )
1453 if(n1b.ge.1 .or. n2a.ge.1)
then
1454 pmom2(:) = quarks(2)%Mom +
summom(gluons,rin,rout)
1455 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
1456 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut ) cycle
1457 if( quarks(2)%PartType.lt.0 )
then
1458 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
1460 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
1463 if( quarks(2)%PartType.lt.0 )
then
1469 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
1470 if(n1a.ge.1 .or. n4b.ge.1)
then
1471 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
1472 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut ) cycle
1473 if( quarks(2)%PartType.lt.0 )
then
1480 tmpquark(1)%Mom => pmom1(:)
1481 tmpquark(1)%Pol =>
ubar0(:)
1482 tmpquark(1)%Mass => quarks(2)%Mass
1483 tmpquark(1)%Mass2=> quarks(2)%Mass2
1485 tmpquark(1)%ExtRef => tmpextref
1486 tmpquark(1)%PartType => quarks(2)%PartType
1494 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1500 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
1501 res(:) = res(:) + tmp(:)
1513 if( (quark1parttype.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType) .AND. &
1514 (quarks(4)%ExtRef.ne.-1 .or. tag_f.ne.1 .or. abs(quark1parttype).ne.abs(quarks(2)%PartType)) &
1522 rout=numglu(1)+numglu(2)+n3a
1526 if ( tag_z .eq. 1 .and. (n1b+numglu(2)+n3a .eq. numglu(0)) .and. quarks(4)%ExtRef .eq. -1)
then
1534 eps2 =
cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
1535 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
1536 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
1538 eps2 = eps2*propfac1
1543 rin =numglu(1)+numglu(2)+n3a+1
1544 rout=numglu(1)+numglu(2)+numglu(3)+n4a
1545 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/numglu(3)+n4a-n3a,n3b,n4a/) )
1546 if(n3b.ge.1 .or. n4a.ge.1)
then
1547 pmom2(:) = quarks(4)%Mom +
summom(gluons,rin,rout)
1548 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(4)%Mass2)
1549 if( abs(
sc_(pmom2,pmom2)-quarks(4)%Mass2).lt.
propcut ) cycle
1550 if( quarks(4)%PartType.lt.0 )
then
1551 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(4)%Mass*ubar1(:))*propfac2
1553 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(4)%Mass*ubar1(:))*propfac2
1556 if( quarks(4)%PartType.lt.0 )
then
1562 pmom1 = quarks(2)%Mom+quarks(3)%Mom+quarks(4)%Mom+
summom(gluons,n1a+1,numglu(1)+numglu(2)+numglu(3)+n4a)
1563 if(n1a.ge.1 .or. n4b.ge.1)
then
1564 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(4)%Mass2)
1565 if( abs(
sc_(pmom1,pmom1)-quarks(4)%Mass2).lt.
propcut ) cycle
1566 if( quarks(4)%PartType.lt.0 )
then
1572 tmpquark(1)%Mom => pmom1(:)
1573 tmpquark(1)%Pol =>
ubar0(:)
1574 tmpquark(1)%Mass => quarks(4)%Mass
1575 tmpquark(1)%Mass2=> quarks(4)%Mass2
1577 tmpquark(1)%ExtRef => tmpextref
1578 tmpquark(1)%PartType => quarks(4)%PartType
1586 rin =numglu(1)+numglu(2)+numglu(3)+n4a+1
1592 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n4b/) )
1593 res(:) = res(:) + tmp(:)
1604 recursive function g(e,k)
result(res)
1606 complex(8),
intent(in) :: e(:,:),k(:,:)
1607 complex(8) :: e1(size(e,dim=1))
1608 complex(8) :: res(size(e,dim=1))
1609 complex(8) :: k1(size(e,dim=1))
1610 complex(8) :: k2(size(e,dim=1))
1611 complex(8) :: k3(size(e,dim=1))
1612 complex(8) :: e2(size(e,dim=1))
1613 complex(8) :: e3(size(e,dim=1))
1614 complex(8) :: tmp(size(e,dim=1))
1615 complex(8) :: k1sq, k2sq, k3sq
1616 integer :: npart, m, m1
1618 npart =
size(e,dim=2)
1620 if (npart == 1)
then
1623 elseif (npart == 2)
then
1624 res =
vggg(e(:,1),k(:,1),e(:,2),k(:,2))
1630 k1=sum(k(:,1:m),dim=2)
1631 k2=sum(k(:,m+1:npart),dim=2)
1633 e1 =
g(e(:,1:m),k(:,1:m))
1634 e2 =
g(e(:,m+1:npart),k(:,m+1:npart))
1636 tmp =
vggg(e1,k1,e2,k2)
1641 tmp = -(0d0,1d0)*tmp/k1sq
1647 if (m + 1 < npart)
then
1650 tmp = -(0d0,1d0)*tmp/k2sq
1658 if (m <= npart-2)
then
1662 e2=
g(e(:,m+1:m1),k(:,m+1:m1))
1663 e3=
g(e(:,m1+1:npart),k(:,m1+1:npart))
1664 k2 = sum(k(:,m+1:m1),dim=2)
1665 k3 = sum(k(:,m1+1:npart),dim=2)
1666 tmp =
vgggg(e1,e2,e3)
1670 tmp = -(0d0,1d0)*tmp/k1sq
1678 tmp = -(0d0,1d0)*tmp/k2sq
1683 if (m1+1 < npart)
then
1686 tmp = -(0d0,1d0)*tmp/k3sq
1704 FUNCTION cur_g(Gluons,NumGlu)
1708 integer :: ngluons,numglu
1709 complex(8) :: glu_subcur(1:
dv,1:36)
1710 complex(8) :: mom_sum(1:
dv,1:36), propfactor,propdenom
1711 integer :: ind0,ind1,ind2,ind3,j,l,mu
1712 integer :: a,b,i1,i2
1715 character :: parts(20)*4
1716 integer :: parti(20)
1719 if( gluons(i1)%ExtRef.eq.-1 )
then
1722 parti(i1)=gluons(i1)%ExtRef
1724 if(i1.eq.numglu-1) print *, parti(1:numglu-1)
1739 glu_subcur(1:
dv,ind0) = gluons(i1)%Pol(1:
dv)
1740 mom_sum(1:
dv,ind0) = gluons(i1)%Mom(1:
dv)
1742 mom_sum(1:
dv,ind0) = mom_sum(1:
dv,ind0-ngluons+i2-i1-1) + gluons(i2)%Mom(1:
dv)
1743 if ( i1 .ne. 1 .or. i2 .ne. ngluons )
then
1744 propdenom = mom_sum(1:
dv,ind0).ndot.mom_sum(1:
dv,ind0)
1745 if( abs(propdenom).lt.
propcut ) cycle
1746 propfactor = (0d0,-1d0)/propdenom
1751 glu_subcur(mu,ind0) = 0d0
1756 glu_subcur(1:
dv,ind0) = glu_subcur(1:
dv,ind0) + &
1758 glu_subcur(1:
dv,ind1),glu_subcur(1:
dv,ind2)) * propfactor
1765 glu_subcur(1:
dv,ind0) = glu_subcur(1:
dv,ind0) + &
1767 glu_subcur(1:
dv,ind1),glu_subcur(1:
dv,ind2),glu_subcur(1:
dv,ind3)) * propfactor
1782 recursive function f(e,k,sp,p,mass,flout,flin,ms)
result(res)
1784 complex(8),
intent(in) :: e(:,:), k(:,:)
1785 complex(8),
intent(in) :: sp(:), p(:)
1786 character,
intent(in) :: flin*3
1787 character,
intent(in) :: flout*3
1788 integer,
intent(in) :: ms
1789 integer :: ms1,m,ng1, ng2
1791 complex(8) :: res(size(sp))
1792 complex(8) :: tmp(size(sp))
1793 complex(8) :: k1(size(p))
1794 complex(8) :: k2(size(p))
1795 complex(8) :: sp2(size(sp))
1796 complex(8) :: sp3(size(sp))
1797 complex(8) :: e1(size(e,dim=1))
1798 complex(8) :: e2(size(e,dim=1))
1799 complex(8) :: k1sq,k2sq,k3sq
1802 ngluon =
size(e,dim=2)
1806 if (flout.ne.flin)
then
1810 if (ng2 < 0)
write(*,*)
'WRONG DEFINITION OF CURRENT A'
1812 if (ngluon == 0)
then
1819 k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
1820 e1=
g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
1823 k2 = sum(k(:,1:ng1+m),dim=2)
1825 k2sq =
sc_(k2,k2)-mass**2
1826 sp2 =
f(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,flout,flin,ng1)
1827 if (ng1 >0.or.m>0) sp2 =
spb2_(sp2,k2)+mass*sp2
1834 tmp = -(0d0,1d0)/k1sq*tmp
1840 if (ng1>0.or.m>0)
then
1842 tmp = (0d0,1d0)/k2sq*tmp
1853 k1 = sum(k(:,1:m),dim=2)
1854 e1=
g(e(:,1:m),k(:,1:m))
1857 k2 = sum(k(:,m+1:ngluon),dim=2)
1859 k2sq =
sc_(k2,k2) - mass**2
1861 sp2=
f(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,flout,flin,ms1)
1863 if (ng2 > 0.or.m < ng1) sp2 =
spb2_(sp2,k2)+mass*sp2
1867 tmp=-(0d0,1d0)/k1sq*tmp
1873 if (ng2 > 0.or. m < ng1)
then
1875 tmp=(0d0,1d0)/k2sq*tmp
1891 recursive function bf(e,k,sp,p,mass,flout,flin,ms)
result(res)
1893 complex(8),
intent(in) :: e(:,:), k(:,:)
1894 complex(8),
intent(in) :: sp(:), p(:)
1895 integer,
intent(in) :: ms
1896 character,
intent(in) :: flout*3
1897 character,
intent(in) :: flin*3
1898 integer :: ms1,m,ng1, ng2
1900 complex(8) :: res(size(sp))
1901 complex(8) :: tmp(size(sp))
1902 complex(8) :: k1(size(p))
1903 complex(8) :: k2(size(p))
1904 complex(8) :: sp2(size(sp))
1905 complex(8) :: sp3(size(sp))
1906 complex(8) :: e1(size(e,dim=1))
1907 complex(8) :: e2(size(e,dim=1))
1908 complex(8) :: k1sq,k2sq,k3sq
1912 if (flout.ne.flin)
then
1916 ngluon =
size(e,dim=2)
1920 if (ng2 < 0)
write(*,*)
'WRONG DEFINITION OF CURRENT B'
1921 if (ngluon == 0)
then
1927 k1 = sum(k(:,ng1+1+m:ngluon),dim=2)
1928 e1=
g(e(:,ng1+1+m:ngluon),k(:,ng1+1+m:ngluon))
1931 k2 = sum(k(:,1:ng1+m),dim=2)
1933 k2sq =
sc_(k2,k2)-mass**2
1934 sp2 =
bf(e(:,1:ng1+m),k(:,1:ng1+m),sp,p,mass,flout,flin,ng1)
1935 if (ng1>0.or.m>0) sp2 =
spi2_(k2,sp2)+mass*sp2
1941 tmp = -(0d0,1d0)/k1sq*tmp
1946 if (ng1>0.or.m>0)
then
1948 tmp = (0d0,1d0)/k2sq*tmp
1962 k1 = sum(k(:,1:m),dim=2)
1963 e1=
g(e(:,1:m),k(:,1:m))
1966 k2 = sum(k(:,m+1:ngluon),dim=2)
1968 k2sq =
sc_(k2,k2) - mass**2
1970 sp2=
bf(e(:,m+1:ngluon),k(:,m+1:ngluon),sp,p,mass,flout,flin,ms1)
1972 if (ng2 > 0.or.m < ng1) sp2 =
spi2_(k2,sp2)+mass*sp2
1978 tmp=-(0d0,1d0)/k1sq*tmp
1984 if (ng2 > 0.or. m < ng1)
then
1986 tmp=(0d0,1d0)/k2sq*tmp
2005 FUNCTION cur_g_2f(Gluons,Quarks,NumGlu)
result(Res)
2007 integer :: numglu(0:3),i,counter
2009 integer :: rin,rout,n1a,n1b,n2a,n2b,n3a,n3b
2010 integer,
target :: tmpextref
2011 complex(8) :: res(1:
dv)
2012 complex(8) :: u1(1:
ds),ubar2(1:
ds)
2013 complex(8),
target :: eps1(1:
dv)
2014 complex(8) :: eps2(1:
dv)
2016 complex(8) :: pmom1(1:
dv),pmom2(1:
dv),pmom4(1:
dv)
2017 complex(8),
target :: pmom3(1:
dv)
2018 complex(8) :: propfac1,propfac2,propfac3,propfac4
2019 integer :: partkey,helkey,currkey,hel_tmp
2025 if( quarks(1)%PartType .ne. -quarks(2)%PartType )
return
2035 pmom1(:) = quarks(1)%Mom(:)+
summom(gluons,rin,rout)
2036 u1(:) =
cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/n1b+n2a,n1b,n2a/))
2037 if(n1b.ge.1 .or. n2a.ge.1)
then
2038 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(1)%Mass2)
2039 if( abs(
sc_(pmom1,pmom1)-quarks(1)%Mass2).lt.
propcut ) cycle
2040 if( quarks(1)%PartType.lt.0 )
then
2041 u1(:) = (-
spi2_(pmom1,u1)+quarks(1)%Mass*u1(:) )*propfac1
2043 u1(:) = (+
spb2_(u1,pmom1)+quarks(1)%Mass*u1(:) )*propfac1
2049 rout=numglu(1)+numglu(2)+n3a
2050 pmom2(:) = quarks(2)%Mom(:)+
summom(gluons,rin,rout)
2051 ubar2(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2b+n3a,n2b,n3a/))
2052 if(n2b.ge.1 .or. n3a.ge.1)
then
2053 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
2054 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut ) cycle
2055 if( quarks(2)%PartType.lt.0 )
then
2056 ubar2(:) = (-
spi2_(pmom2,ubar2)+quarks(2)%Mass*ubar2(:) )*propfac2
2058 ubar2(:) = (+
spb2_(ubar2,pmom2)+quarks(2)%Mass*ubar2(:) )*propfac2
2062 if( quarks(1)%PartType.lt.0 )
then
2063 eps1(:)= -
vbqq(
dv,ubar2,u1)
2065 eps1(:)= +
vbqq(
dv,u1,ubar2)
2067 pmom3(:) = quarks(1)%Mom(:)+quarks(2)%Mom(:) +
summom(gluons,n1a+1,numglu(1)+numglu(2)+n3a)
2075 tmpgluons(counter)%Mom => pmom3(:)
2076 tmpgluons(counter)%Pol => eps1(:)
2078 tmpgluons(counter)%ExtRef => tmpextref
2080 rin =numglu(1)+numglu(2)+n3a+1
2081 rout=numglu(1)+numglu(2)+numglu(3)
2086 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+n1a+n3b+1)
2089 if(n1a.ge.1 .or. n3b.ge.1)
then
2090 propfac3 = (0d0,-1d0)/
sc_(pmom3,pmom3)
2092 eps2(:) = eps2(:)*propfac3
2095 res(:) = res(:) + eps2(:)
2109 FUNCTION cur_g_2fv(Gluons,Quarks,Boson,NumGlu)
result(Res)
2111 integer :: numglu(0:3),i,counter
2113 integer :: rin,rout,n1a,n1b,n2a,n2b,n3a,n3b
2114 integer,
target :: tmpextref
2115 complex(8) :: res(1:
dv)
2116 complex(8) :: u1(1:
ds),ubar2(1:
ds)
2117 complex(8),
target :: eps1(1:
dv)
2118 complex(8) :: eps2(1:
dv)
2120 complex(8) :: pmom1(1:
dv),pmom2(1:
dv),pmom4(1:
dv)
2121 complex(8),
target :: pmom3(1:
dv)
2122 complex(8) :: propfac1,propfac2,propfac3,propfac4
2123 integer :: partkey,helkey,currkey,hel_tmp
2129 if(quarks(1)%PartType*quarks(2)%PartType.ge.0) print *,
"Error in cur_g_2f: wrong PartTypes"
2132 if( numglu(0)-1-numglu(1)-numglu(2)-numglu(3).ne.0 ) print *,
"wrong number of gluons in cur_g_2f"
2136 if( quarks(1)%PartType .ne. -quarks(2)%PartType )
return
2147 pmom1(:) = quarks(1)%Mom(:) +
summom(gluons,rin,rout) + boson%Mom(:)
2148 u1(:) =
cur_f_2fv(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,boson,(/n1b+n2a,n1b,n2a/))
2150 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(1)%Mass2)
2151 if( abs(
sc_(pmom1,pmom1)-quarks(1)%Mass2).lt.
propcut )
then
2154 if( quarks(1)%PartType.lt.0 )
then
2155 u1(:) = (-
spi2_(pmom1,u1)+quarks(1)%Mass*u1(:) )*propfac1
2157 u1(:) = (+
spb2_(u1,pmom1)+quarks(1)%Mass*u1(:) )*propfac1
2163 rout=numglu(1)+numglu(2)+n3a
2164 pmom2(:) = quarks(2)%Mom(:)+
summom(gluons,rin,rout)
2165 ubar2(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n2b+n3a,n2b,n3a/))
2166 if(n2b.ge.1 .or. n3a.ge.1)
then
2167 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
2168 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut )
then
2171 if( quarks(2)%PartType.lt.0 )
then
2172 ubar2(:) = (-
spi2_(pmom2,ubar2)+quarks(2)%Mass*ubar2(:) )*propfac2
2174 ubar2(:) = (+
spb2_(ubar2,pmom2)+quarks(2)%Mass*ubar2(:) )*propfac2
2179 if( quarks(1)%PartType.lt.0 )
then
2180 eps1(:)= -
vbqq(
dv,ubar2,u1)
2182 eps1(:)= +
vbqq(
dv,u1,ubar2)
2192 pmom1(:) = quarks(1)%Mom(:) +
summom(gluons,rin,rout)
2193 u1(:) =
cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/n1b+n2a,n1b,n2a/))
2194 if(n1b.ge.1 .or. n2a.ge.1)
then
2195 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(1)%Mass2)
2196 if( abs(
sc_(pmom1,pmom1)-quarks(1)%Mass2).lt.
propcut ) cycle
2197 if( quarks(1)%PartType.lt.0 )
then
2198 u1(:) = (-
spi2_(pmom1,u1)+quarks(1)%Mass*u1(:) )*propfac1
2200 u1(:) = (+
spb2_(u1,pmom1)+quarks(1)%Mass*u1(:) )*propfac1
2206 rout=numglu(1)+numglu(2)+n3a
2207 pmom2(:) = quarks(2)%Mom(:)+
summom(gluons,rin,rout) + boson%Mom(:)
2208 ubar2(:) =
cur_f_2fv(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,boson,(/n2b+n3a,n2b,n3a/))
2210 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
2211 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut ) cycle
2212 if( quarks(2)%PartType.lt.0 )
then
2213 ubar2(:) = (-
spi2_(pmom2,ubar2)+quarks(2)%Mass*ubar2(:) )*propfac2
2215 ubar2(:) = (+
spb2_(ubar2,pmom2)+quarks(2)%Mass*ubar2(:) )*propfac2
2219 if( quarks(1)%PartType.lt.0 )
then
2220 eps1(:)= eps1(:) -
vbqq(
dv,ubar2,u1)
2222 eps1(:)= eps1(:) +
vbqq(
dv,u1,ubar2)
2226 pmom3(:) = quarks(1)%Mom(:)+quarks(2)%Mom(:) +
summom(gluons,n1a+1,numglu(1)+numglu(2)+n3a) + boson%Mom(:)
2234 tmpgluons(counter)%Mom => pmom3(:)
2235 tmpgluons(counter)%Pol => eps1(:)
2237 tmpgluons(counter)%ExtRef => tmpextref
2239 rin =numglu(1)+numglu(2)+n3a+1
2240 rout=numglu(1)+numglu(2)+numglu(3)
2245 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+n1a+n3b+1)
2248 if(n1a.ge.1 .or. n3b.ge.1)
then
2249 propfac3 = (0d0,-1d0)/
sc_(pmom3,pmom3)
2251 eps2(:) = eps2(:)*propfac3
2254 res(:) = res(:) + eps2(:)
2272 FUNCTION cur_f_2f(Gluons,Quarks,Quark1PartType,NumGlu)
result(Res)
2274 complex(8) :: res(1:
ds)
2275 integer :: numglu(0:2),i,rin,rout,quark1parttype
2277 complex(8) :: glumom(1:
dv,numglu(0)), quarkmom(1:
dv)
2278 complex(8) :: glupol(1:
dv,numglu(0)), quarkpol(1:
ds)
2279 character :: ferfla1*3,ferfla2*3
2280 integer :: partkey,helkey,currkey,hel_tmp
2288 if( quarks(2)%PartType .eq.0 .or. .not.
isaquark(quarks(2)%PartType))
then
2289 print *,
"Error in cur_f_2f"
2294 if( numglu(0)-numglu(1)-numglu(2).ne.0 ) print *,
"wrong number of gluons in cur_f_2f"
2295 if( quarks(2)%PartType.ne.-quark1parttype ) print *,
"unequal flavors in cur_f_2f"
2300 glumom(1:
dv,i) = gluons(i)%Mom(1:
dv)
2301 glupol(1:
dv,i) = gluons(i)%Pol(1:
dv)
2303 quarkmom(1:
dv) = quarks(2)%Mom(1:
dv)
2304 quarkpol(1:
ds) = quarks(2)%Pol(1:
ds)
2306 if( abs(quark1parttype).eq.5 )
then
2308 elseif( abs(quark1parttype).eq.6 )
then
2310 elseif( abs(quark1parttype).eq.3 )
then
2316 if( abs(quarks(2)%PartType).eq.5 )
then
2318 elseif( abs(quarks(2)%PartType).eq.6 )
then
2320 elseif( abs(quarks(2)%PartType).eq.3 )
then
2328 if( quarks(2)%PartType .gt.0 )
then
2329 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))
2331 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))
2343 FUNCTION cur_g_4f(Gluons,Quarks,NumGlu)
result(res)
2345 integer,
intent(in) :: numglu(0:5)
2347 integer :: na,nb,nc,nd,
ne,
nf,ng,nh,ni,nj,nk
2349 integer :: tag_f,counter,i
2350 complex(8) :: res(
dv)
2351 type(
ptrtoparticle) :: tmpgluons(1:2+numglu(1)+numglu(3)+numglu(5))
2352 complex(8),
target :: tmpmom1(1:
dv),tmpmom2(1:
dv)
2353 integer,
target :: tmpextref1,tmpextref2
2354 complex(8),
target :: eps1(1:
dv)
2355 complex(8),
target :: eps2(1:
dv)
2356 complex(8) :: eps3(1:
dv)
2357 complex(8) :: u1(1:
ds)
2358 complex(8) :: ubar2(1:
ds)
2359 complex(8) :: propfac1,propfac2,propfac3,propfac4
2360 complex(8) :: pmom1(1:
dv)
2361 complex(8) :: pmom2(1:
dv)
2362 complex(8) :: pmom3(1:
dv)
2363 complex(8) :: pmom4(1:
dv)
2366 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"
2370 if( (quarks(1)%PartType.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType) &
2371 .OR. (quarks(1)%PartType.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType) )
then
2379 rin = numglu(1)+nc+1
2380 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
2381 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)
2382 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
2383 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(1)%Mass2)
2384 if( abs(
sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.
propcut ) cycle
2385 if( quarks(1)%PartType.lt.0 )
then
2386 u1 = (
spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
2388 u1 = (-
spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
2393 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
2394 pmom3 =
summom(gluons,rin,rout) + quarks(1)%Mom
2395 if( nb.ge.1 .or. nc.ge.1 )
then
2396 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(1)%Mass2)
2397 if( abs(
sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.
propcut ) cycle
2398 if( quarks(1)%PartType.lt.0 )
then
2399 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
2401 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
2405 if( quarks(1)%PartType.lt.0 )
then
2406 eps1 = -
vbqq(
dv,u1,ubar2)
2408 eps1 = +
vbqq(
dv,ubar2,u1)
2418 tmpmom1(:) = pmom2(:)+pmom3(:)
2420 tmpgluons(counter)%Mom => tmpmom1(:)
2421 tmpgluons(counter)%Pol => eps1(:)
2422 tmpgluons(counter)%ExtRef => tmpextref1
2424 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
2430 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
2432 if( na.ge.1 .or.
nf.ge.1 )
then
2433 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
2434 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
2435 eps2 = eps2*propfac1
2453 rout= numglu(1)+numglu(2)+numglu(3)+nc
2454 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)
2455 pmom2 =
summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom
2456 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
2457 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
2458 if( quarks(4)%PartType.lt.0 )
then
2459 u1 = (+
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
2461 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
2463 rin = numglu(1)+numglu(2)+numglu(3)+nc+1
2464 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
2465 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/nd+
ne,nd,
ne/))
2466 pmom3 =
summom(gluons,rin,rout) + quarks(4)%Mom
2467 if( nd.ge.1 .or.
ne.ge.1 )
then
2468 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(4)%Mass2)
2469 if( abs(
sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.
propcut ) cycle
2470 if( quarks(4)%PartType.lt.0 )
then
2471 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
2473 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
2477 if( quarks(4)%PartType.lt.0 )
then
2478 eps1 = +
vbqq(
dv,u1,ubar2)
2480 eps1 = -
vbqq(
dv,ubar2,u1)
2490 tmpmom1(:) = pmom2(:)+pmom3(:)
2492 tmpgluons(counter)%Mom => tmpmom1(:)
2493 tmpgluons(counter)%Pol => eps1(:)
2494 tmpgluons(counter)%ExtRef => tmpextref1
2496 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
2502 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
2504 if( na.ge.1 .or.
nf.ge.1 )
then
2505 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
2506 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
2507 eps2 = eps2*propfac1
2518 if( quarks(1)%PartType.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType)
then
2523 do nf=0,numglu(3)-
ne
2534 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
2535 pmom1 =
summom(gluons,rin,rout) + quarks(1)%Mom
2536 if( nb.ge.1 .or. nc.ge.1 )
then
2537 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1) - quarks(1)%Mass2)
2538 if( abs(
sc_(pmom1,pmom1) - quarks(1)%Mass2).lt.
propcut ) cycle
2539 if( quarks(1)%PartType.lt.0 )
then
2540 ubar2 = (-
spi2_(pmom1,ubar2) + quarks(1)%Mass*ubar2 )*propfac1
2542 ubar2 = (+
spb2_(ubar2,pmom1) + quarks(1)%Mass*ubar2 )*propfac1
2545 rin = numglu(1)+nc+1
2546 rout= numglu(1)+numglu(2)+
ne
2547 u1 =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/nd+
ne,nd,
ne/))
2548 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom
2549 if( nd.ge.1 .or.
ne.ge.1 )
then
2550 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
2551 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
2552 if( quarks(2)%PartType.lt.0 )
then
2553 u1 = (-
spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
2555 u1 = (+
spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
2559 if( quarks(2)%PartType.lt.0 )
then
2560 eps1 = +
vbqq(
dv,ubar2,u1)
2562 eps1 = -
vbqq(
dv,u1,ubar2)
2564 tmpmom1 = pmom1 + pmom2
2565 propfac3 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
2566 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
2567 eps1 = eps1*propfac3
2570 rin = numglu(1)+numglu(2)+
ne+
nf+1
2571 rout= numglu(1)+numglu(2)+numglu(3)+nh
2572 u1 =
cur_f_2f(gluons(rin:rout),quarks(3:3),-quarks(3)%PartType,(/ng+nh,ng,nh/))
2573 pmom3 =
summom(gluons,rin,rout) + quarks(3)%Mom
2574 if( ng.ge.1 .or. nh.ge.1 )
then
2575 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(3)%Mass2)
2576 if( abs(
sc_(pmom3,pmom3) - quarks(3)%Mass2).lt.
propcut ) cycle
2577 if( quarks(3)%PartType.lt.0 )
then
2578 u1 = (-
spi2_(pmom3,u1) + quarks(3)%Mass*u1 )*propfac3
2580 u1 = (+
spb2_(u1,pmom3) + quarks(3)%Mass*u1 )*propfac3
2583 rin = numglu(1)+numglu(2)+numglu(3)+nh+1
2584 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj
2585 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/ni+nj,ni,nj/))
2586 pmom4 =
summom(gluons,rin,rout) + quarks(4)%Mom
2587 if( ni.ge.1 .or. nj.ge.1 )
then
2588 propfac4 = (0d0,1d0)/(
sc_(pmom4,pmom4) - quarks(4)%Mass2)
2589 if( abs(
sc_(pmom4,pmom4) - quarks(4)%Mass2).lt.
propcut ) cycle
2590 if( quarks(4)%PartType.lt.0 )
then
2591 ubar2 = (-
spi2_(pmom4,ubar2) + quarks(4)%Mass*ubar2 )*propfac4
2593 ubar2 = (+
spb2_(ubar2,pmom4) + quarks(4)%Mass*ubar2 )*propfac4
2597 if( quarks(4)%PartType.lt.0 )
then
2598 eps2 = +
vbqq(
dv,u1,ubar2)
2600 eps2 = -
vbqq(
dv,ubar2,u1)
2602 tmpmom2 = pmom3 + pmom4
2603 propfac1 = (0d0,-1d0)/
sc_(tmpmom2,tmpmom2)
2604 if( abs(
sc_(tmpmom2,tmpmom2)).lt.
propcut ) cycle
2605 eps2 = eps2*propfac1
2616 tmpgluons(counter)%Mom => tmpmom1(:)
2617 tmpgluons(counter)%Pol => eps1(:)
2618 tmpgluons(counter)%ExtRef => tmpextref1
2620 rin =numglu(1)+numglu(2)+
ne+1
2621 rout=numglu(1)+numglu(2)+
ne+
nf
2626 tmpgluons(counter)%Mom => tmpmom2(:)
2627 tmpgluons(counter)%Pol => eps2(:)
2628 tmpgluons(counter)%ExtRef => tmpextref1
2630 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj+1
2636 eps3(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+nk+2)
2653 FUNCTION cur_f_6f(Gluons,Quarks,Quark1PartType,NumGlu,tag_f,tag_Z)
result(res)
2655 integer :: numglu(0:6),quark1parttype,tag_f,tag_z
2657 integer,
target :: tmpparttype,tmpextref
2658 complex(8) :: res(1:
ds),tmp(1:
ds)
2660 complex(8) :: u1(1:
ds),ubar1(1:
ds)
2661 complex(8),
target ::
ubar0(1:
ds)
2662 complex(8) :: eps1(1:
dv)
2663 complex(8) :: eps2(1:
dv)
2664 type(
ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(6)),tmpquark(1:1)
2665 complex(8) :: propfac1,propfac2
2666 complex(8),
target :: pmom1(1:
dv)
2667 complex(8),
target :: pmom2(1:
dv)
2668 integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b,n5a,n5b,n6a,n6b
2669 integer :: rin,rout,i,counter
2673 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"
2681 if( quark1parttype.eq.-quarks(2)%PartType .AND. (quarks(3)%PartType.eq.-quarks(4)%PartType .or. quarks(3)%PartType.eq.-quarks(6)%PartType) &
2682 .AND. (quarks(2)%ExtRef.ne.-1 .or. tag_f.ne.1) &
2690 rin =numglu(1)+n2a+1
2691 rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2692 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/))
2693 pmom1(:) =
summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
2694 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
2696 eps2 = eps2*propfac1
2702 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n1b+n2a,n1b,n2a/) )
2703 if(n1b.ge.1 .or. n2a.ge.1)
then
2704 pmom2(:) = quarks(2)%Mom +
summom(gluons,rin,rout)
2705 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
2706 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut ) cycle
2707 if( quarks(2)%PartType.lt.0 )
then
2708 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
2710 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
2713 if( quarks(2)%PartType.lt.0 )
then
2720 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2721 pmom1 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
2722 if(n1a.ge.1 .or. n6b.ge.1)
then
2723 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
2724 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut ) cycle
2725 if( quarks(2)%PartType.lt.0 )
then
2732 tmpquark(1)%Mom => pmom1(:)
2733 tmpquark(1)%Pol =>
ubar0(:)
2734 tmpquark(1)%Mass => quarks(2)%Mass
2735 tmpquark(1)%Mass2=> quarks(2)%Mass2
2737 tmpquark(1)%ExtRef => tmpextref
2738 tmpquark(1)%PartType => quarks(2)%PartType
2746 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
2752 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n6b/) )
2754 res(:) = res(:) + tmp(:)
2765 if( quark1parttype.eq.-quarks(6)%PartType .AND. (quarks(2)%PartType.eq.-quarks(5)%PartType .or. quarks(2)%PartType.eq.-quarks(3)%PartType) &
2766 .AND. (quarks(6)%ExtRef.ne.-1 .or. tag_f.ne.1) &
2774 rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+n5a
2778 if ( tag_z .eq. 1 .and. (n1b+numglu(2)+numglu(3)+numglu(4)+n5a .eq. numglu(0)) )
then
2783 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/))
2784 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom
2785 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
2787 eps2 = eps2*propfac1
2792 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+n5a+1
2793 rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2794 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(6:6),-quarks(6)%PartType,(/n5b+n6a,n5b,n6a/) )
2795 if(n5b.ge.1 .or. n6a.ge.1)
then
2796 pmom2(:) =
summom(gluons,rin,rout) + quarks(6)%Mom
2797 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(6)%Mass2)
2798 if( abs(
sc_(pmom2,pmom2)-quarks(6)%Mass2).lt.
propcut ) cycle
2799 if( quarks(6)%PartType.lt.0 )
then
2800 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(6)%Mass*ubar1(:))*propfac2
2802 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(6)%Mass*ubar1(:))*propfac2
2805 if( quarks(6)%PartType.lt.0 )
then
2812 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2813 pmom1 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
2814 if(n1a.ge.1 .or. n6b.ge.1)
then
2815 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(6)%Mass2)
2816 if( abs(
sc_(pmom1,pmom1)-quarks(6)%Mass2).lt.
propcut ) cycle
2817 if( quarks(6)%PartType.lt.0 )
then
2823 tmpquark(1)%Mom => pmom1(:)
2824 tmpquark(1)%Pol =>
ubar0(:)
2825 tmpquark(1)%Mass => quarks(6)%Mass
2826 tmpquark(1)%Mass2=> quarks(6)%Mass2
2828 tmpquark(1)%ExtRef => tmpextref
2829 tmpquark(1)%PartType => quarks(6)%PartType
2837 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
2843 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
2845 res(:) = res(:) + tmp(:)
2866 if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
2867 ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
2868 .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
2878 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
2879 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2880 eps2 =
cur_g_2f(gluons(rin:rout),quarks(5:6),(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
2881 pmom1(:) =
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
2882 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
2884 eps2 = eps2*propfac1
2887 rout= numglu(1)+numglu(2)+numglu(3)+n4a
2888 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)
2889 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
2891 if( quark1parttype.eq.-quarks(2)%PartType)
then
2892 if( quarks(2)%PartType.lt.0 )
then
2893 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
2894 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
2895 u1 = (-
spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
2897 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
2898 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2899 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
2900 if( n1a.ge.1 .or. n6b.ge.1 )
then
2901 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
2902 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
2905 elseif( quarks(2)%PartType.gt.0 )
then
2906 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
2907 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
2908 u1 = (
spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
2910 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
2911 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2912 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
2913 if( n1a.ge.1 .or. n6b.ge.1 )
then
2914 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
2915 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
2919 tmpquark(1)%Mom => pmom2(:)
2920 tmpquark(1)%Pol =>
ubar0(:)
2921 tmpquark(1)%Mass => quarks(2)%Mass
2922 tmpquark(1)%Mass2=> quarks(2)%Mass2
2923 elseif( quark1parttype.eq.-quarks(4)%PartType)
then
2924 if( quarks(4)%PartType.lt.0 )
then
2925 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
2926 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
2927 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
2929 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
2930 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2931 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
2932 if( n1a.ge.1 .or. n6b.ge.1 )
then
2933 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
2934 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
2937 elseif( quarks(4)%PartType.gt.0 )
then
2938 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
2939 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
2940 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
2942 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
2943 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
2944 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
2945 if( n1a.ge.1 .or. n6b.ge.1 )
then
2946 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
2947 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
2951 tmpquark(1)%Mom => pmom2(:)
2952 tmpquark(1)%Pol =>
ubar0(:)
2953 tmpquark(1)%Mass => quarks(4)%Mass
2954 tmpquark(1)%Mass2=> quarks(4)%Mass2
2957 tmpquark(1)%ExtRef => tmpextref
2958 tmpparttype = -quark1parttype
2959 tmpquark(1)%PartType => tmpparttype
2967 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
2973 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
2975 res(:) = res(:) + tmp(:)
2985 if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
2986 (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
2987 .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1)) ))
then
3001 rout= numglu(1)+numglu(2)+n3a
3004 eps2 =
cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
3005 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
3006 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
3008 eps2 = eps2*propfac1
3010 rin = numglu(1)+numglu(2)+n3a+1
3011 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
3012 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)
3013 pmom2 =
summom(gluons,rin,rout) + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
3015 if( quark1parttype.eq.-quarks(4)%PartType )
then
3016 if( quarks(4)%PartType.lt.0 )
then
3017 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3018 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
3019 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
3022 rout= numglu(1)+numglu(2)+n3a
3023 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
3024 if( n1a.ge.1 .or. n6b.ge.1 )
then
3025 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3026 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
3029 elseif( quarks(4)%PartType.gt.0 )
then
3030 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3031 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
3032 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
3035 rout= numglu(1)+numglu(2)+n3a
3036 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
3037 if( n1a.ge.1 .or. n6b.ge.1 )
then
3038 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3039 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
3043 tmpquark(1)%Mom => pmom2(:)
3044 tmpquark(1)%Pol =>
ubar0(:)
3045 tmpquark(1)%Mass => quarks(4)%Mass
3046 tmpquark(1)%Mass2=> quarks(4)%Mass2
3047 elseif(quark1parttype.eq.-quarks(6)%PartType)
then
3048 if( quarks(6)%PartType.lt.0 )
then
3049 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
3050 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
3051 u1 = (-
spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
3054 rout= numglu(1)+numglu(2)+n3a
3055 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
3056 if( n1a.ge.1 .or. n6b.ge.1 )
then
3057 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
3058 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
3061 elseif( quarks(6)%PartType.gt.0 )
then
3062 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
3063 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
3064 u1 = (
spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
3067 rout= numglu(1)+numglu(2)+n3a
3068 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
3069 if( n1a.ge.1 .or. n6b.ge.1 )
then
3070 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
3071 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
3075 tmpquark(1)%Mom => pmom2(:)
3076 tmpquark(1)%Pol =>
ubar0(:)
3077 tmpquark(1)%Mass => quarks(6)%Mass
3078 tmpquark(1)%Mass2=> quarks(6)%Mass2
3081 tmpquark(1)%ExtRef => tmpextref
3082 tmpparttype = -quark1parttype
3083 tmpquark(1)%PartType => tmpparttype
3091 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
3097 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/))
3099 res(:) = res(:) + tmp(:)
3123 FUNCTION cur_g_4fv(Gluons,Quarks,Boson,BosonVertex,NumGlu)
result(res)
3125 integer,
intent(in) :: numglu(0:5),bosonvertex
3127 integer :: na,nb,nc,nd,
ne,
nf,ng,nh,ni,nj,nk,bosonvertex_mod
3129 integer :: tag_f,counter,i
3130 complex(8) :: res(
dv)
3131 type(
ptrtoparticle) :: tmpgluons(1:2+numglu(1)+numglu(3)+numglu(5))
3132 complex(8),
target :: tmpmom1(1:
dv),tmpmom2(1:
dv)
3133 integer,
target :: tmpextref1,tmpextref2
3134 complex(8),
target :: eps1(1:
dv)
3135 complex(8),
target :: eps2(1:
dv)
3136 complex(8) :: eps3(1:
dv)
3137 complex(8) :: u1(1:
ds)
3138 complex(8) :: ubar2(1:
ds)
3139 complex(8) :: propfac1,propfac2,propfac3,propfac4
3140 complex(8) :: pmom1(1:
dv)
3141 complex(8) :: pmom2(1:
dv)
3142 complex(8) :: pmom3(1:
dv)
3143 complex(8) :: pmom4(1:
dv)
3146 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"
3150 if (quarks(1)%PartType.eq.-quarks(4)%PartType .and. quarks(2)%PartType.eq.-quarks(3)%PartType)
then
3157 if( bosonvertex.eq.1)
then
3165 rin = numglu(1)+nc+1
3166 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3167 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)
3168 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
3169 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(1)%Mass2)
3170 if( abs(
sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.
propcut ) cycle
3171 if( quarks(1)%PartType.lt.0 )
then
3172 u1 = (
spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
3174 u1 = (-
spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
3179 ubar2 =
cur_f_2fv(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,boson,(/nb+nc,nb,nc/))
3180 pmom3 =
summom(gluons,rin,rout) + quarks(1)%Mom + boson%Mom
3181 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(1)%Mass2)
3182 if( abs(
sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.
propcut ) cycle
3183 if( quarks(1)%PartType.lt.0 )
then
3184 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
3186 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
3189 if( quarks(1)%PartType.lt.0 )
then
3190 eps1 = -
vbqq(
dv,u1,ubar2)
3192 eps1 = +
vbqq(
dv,ubar2,u1)
3202 tmpmom1(:) = pmom2(:)+pmom3(:)
3204 tmpgluons(counter)%Mom => tmpmom1(:)
3205 tmpgluons(counter)%Pol => eps1(:)
3206 tmpgluons(counter)%ExtRef => tmpextref1
3208 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3214 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3216 if( na.ge.1 .or.
nf.ge.1 )
then
3217 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3218 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3219 eps2 = eps2*propfac1
3228 if( bosonvertex.eq.1 .or. bosonvertex.eq.2 .or. bosonvertex .eq. 3)
then
3238 rout= numglu(1)+numglu(2)+numglu(3)+nc
3239 bosonvertex_mod = bosonvertex + 1
3241 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)
3242 pmom2 =
summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
3243 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3244 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
3245 if( quarks(4)%PartType.lt.0 )
then
3246 u1 = (+
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
3248 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
3251 rin = numglu(1)+numglu(2)+numglu(3)+nc+1
3252 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3253 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/nd+
ne,nd,
ne/))
3254 pmom3 =
summom(gluons,rin,rout) + quarks(4)%Mom
3255 if( nd.ge.1 .or.
ne.ge.1 )
then
3256 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(4)%Mass2)
3257 if( abs(
sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.
propcut ) cycle
3258 if( quarks(4)%PartType.lt.0 )
then
3259 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
3261 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
3265 if( quarks(4)%PartType.lt.0 )
then
3266 eps1 = +
vbqq(
dv,u1,ubar2)
3268 eps1 = -
vbqq(
dv,ubar2,u1)
3278 tmpmom1(:) = pmom2(:)+pmom3(:)
3280 tmpgluons(counter)%Mom => tmpmom1(:)
3281 tmpgluons(counter)%Pol => eps1(:)
3282 tmpgluons(counter)%ExtRef => tmpextref1
3284 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3290 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3292 if( na.ge.1 .or.
nf.ge.1 )
then
3293 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3294 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3295 eps2 = eps2*propfac1
3307 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 2 .or. bosonvertex .eq. 3 )
then
3315 rin = numglu(1)+nc+1
3316 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3317 bosonvertex_mod=bosonvertex
3318 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)
3319 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + boson%Mom
3320 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(1)%Mass2)
3321 if( abs(
sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.
propcut ) cycle
3322 if( quarks(1)%PartType.lt.0 )
then
3323 u1 = (
spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
3325 u1 = (-
spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
3330 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
3331 pmom3 =
summom(gluons,rin,rout) + quarks(1)%Mom
3332 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(1)%Mass2)
3333 if ( nb .ge. 1 .or. nc .ge. 1)
then
3334 if( abs(
sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.
propcut ) cycle
3335 if( quarks(1)%PartType.lt.0 )
then
3336 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
3338 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
3342 if( quarks(1)%PartType.lt.0 )
then
3343 eps1 = -
vbqq(
dv,u1,ubar2)
3345 eps1 = +
vbqq(
dv,ubar2,u1)
3355 tmpmom1(:) = pmom2(:)+pmom3(:)
3357 tmpgluons(counter)%Mom => tmpmom1(:)
3358 tmpgluons(counter)%Pol => eps1(:)
3359 tmpgluons(counter)%ExtRef => tmpextref1
3361 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3367 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3369 if( na.ge.1 .or.
nf.ge.1 )
then
3370 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3371 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3372 eps2 = eps2*propfac1
3382 if( bosonvertex.eq.3)
then
3392 rout= numglu(1)+numglu(2)+numglu(3)+nc
3393 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)
3394 pmom2 =
summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom
3395 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3396 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
3397 if( quarks(4)%PartType.lt.0 )
then
3398 u1 = (+
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
3400 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
3403 rin = numglu(1)+numglu(2)+numglu(3)+nc+1
3404 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3405 bosonvertex_mod = bosonvertex - 3
3406 ubar2 =
cur_f_2fv(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,boson,(/nd+
ne,nd,
ne/))
3407 pmom3 =
summom(gluons,rin,rout) + quarks(4)%Mom + boson%Mom
3408 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(4)%Mass2)
3409 if( abs(
sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.
propcut ) cycle
3410 if( quarks(4)%PartType.lt.0 )
then
3411 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
3413 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
3416 if( quarks(4)%PartType.lt.0 )
then
3417 eps1 = +
vbqq(
dv,u1,ubar2)
3419 eps1 = -
vbqq(
dv,ubar2,u1)
3429 tmpmom1(:) = pmom2(:)+pmom3(:)
3431 tmpgluons(counter)%Mom => tmpmom1(:)
3432 tmpgluons(counter)%Pol => eps1(:)
3433 tmpgluons(counter)%ExtRef => tmpextref1
3435 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3441 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3443 if( na.ge.1 .or.
nf.ge.1 )
then
3444 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3445 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3446 eps2 = eps2*propfac1
3456 if (bosonvertex .ne. 1 .and. bosonvertex .ne. 2 .and. bosonvertex .ne. 3)
then
3457 print *,
'cur_g_4fV not implemented for this flavor choice and BosonVertex=', bosonvertex
3463 if( quarks(1)%PartType.eq.-quarks(2)%PartType .and. quarks(3)%PartType.eq.-quarks(4)%PartType)
then
3465 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 3)
then
3473 rin = numglu(1)+nc+1
3474 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3475 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)
3476 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + boson%Mom
3477 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(1)%Mass2)
3478 if( abs(
sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.
propcut ) cycle
3479 if( quarks(1)%PartType.lt.0 )
then
3480 u1 = (
spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
3482 u1 = (-
spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
3487 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
3488 pmom3 =
summom(gluons,rin,rout) + quarks(1)%Mom
3489 if( nb.ge.1 .or. nc.ge.1 )
then
3490 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(1)%Mass2)
3491 if( abs(
sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.
propcut ) cycle
3492 if( quarks(1)%PartType.lt.0 )
then
3493 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
3495 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
3499 if( quarks(1)%PartType.lt.0 )
then
3500 eps1 = -
vbqq(
dv,u1,ubar2)
3502 eps1 = +
vbqq(
dv,ubar2,u1)
3512 tmpmom1(:) = pmom2(:)+pmom3(:)
3514 tmpgluons(counter)%Mom => tmpmom1(:)
3515 tmpgluons(counter)%Pol => eps1(:)
3516 tmpgluons(counter)%ExtRef => tmpextref1
3518 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3524 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3526 if( na.ge.1 .or.
nf.ge.1 )
then
3527 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3528 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3529 eps2 = eps2*propfac1
3538 if (bosonvertex .eq. 1)
then
3545 rin = numglu(1)+nc+1
3546 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3547 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)
3548 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
3549 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(1)%Mass2)
3550 if( abs(
sc_(pmom2,pmom2) - quarks(1)%Mass2).lt.
propcut ) cycle
3551 if( quarks(1)%PartType.lt.0 )
then
3552 u1 = (
spb2_(u1,pmom2) + quarks(1)%Mass*u1 )*propfac2
3554 u1 = (-
spi2_(pmom2,u1) + quarks(1)%Mass*u1 )*propfac2
3559 ubar2 =
cur_f_2fv(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,boson,(/nb+nc,nb,nc/))
3560 pmom3 =
summom(gluons,rin,rout) + quarks(1)%Mom + boson%Mom
3561 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(1)%Mass2)
3562 if( abs(
sc_(pmom3,pmom3) - quarks(1)%Mass2).lt.
propcut ) cycle
3563 if( quarks(1)%PartType.lt.0 )
then
3564 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(1)%Mass*ubar2 )*propfac3
3566 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(1)%Mass*ubar2 )*propfac3
3569 if( quarks(1)%PartType.lt.0 )
then
3570 eps1 = -
vbqq(
dv,u1,ubar2)
3572 eps1 = +
vbqq(
dv,ubar2,u1)
3582 tmpmom1(:) = pmom2(:)+pmom3(:)
3584 tmpgluons(counter)%Mom => tmpmom1(:)
3585 tmpgluons(counter)%Pol => eps1(:)
3586 tmpgluons(counter)%ExtRef => tmpextref1
3588 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3594 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3596 if( na.ge.1 .or.
nf.ge.1 )
then
3597 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3598 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3599 eps2 = eps2*propfac1
3609 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 3)
then
3619 rout= numglu(1)+numglu(2)+numglu(3)+nc
3620 bosonvertex_mod=bosonvertex+1
3621 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)
3622 pmom2 =
summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
3623 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3624 if ( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut )
then
3627 if( quarks(4)%PartType.lt.0 )
then
3628 u1 = (+
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
3630 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
3632 rin = numglu(1)+numglu(2)+numglu(3)+nc+1
3633 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3634 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/nd+
ne,nd,
ne/))
3635 pmom3 =
summom(gluons,rin,rout) + quarks(4)%Mom
3636 if( nd.ge.1 .or.
ne.ge.1 )
then
3637 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(4)%Mass2)
3638 if( abs(
sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.
propcut ) cycle
3639 if( quarks(4)%PartType.lt.0 )
then
3640 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
3642 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
3646 if( quarks(4)%PartType.lt.0 )
then
3647 eps1 = +
vbqq(
dv,u1,ubar2)
3649 eps1 = -
vbqq(
dv,ubar2,u1)
3659 tmpmom1(:) = pmom2(:)+pmom3(:)
3661 tmpgluons(counter)%Mom => tmpmom1(:)
3662 tmpgluons(counter)%Pol => eps1(:)
3663 tmpgluons(counter)%ExtRef => tmpextref1
3665 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3671 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3673 if( na.ge.1 .or.
nf.ge.1 )
then
3674 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3675 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3676 eps2 = eps2*propfac1
3685 if (bosonvertex .eq. 1)
then
3690 do nf=0,numglu(3)-
ne
3700 rout= numglu(1)+numglu(2)+
ne
3701 eps1 =
cur_g_2fv(gluons(rin:rout),quarks(1:2),boson,(/1+nb+nc+nd+
ne,nc,nc+nd,
ne/) )
3702 tmpmom1 =
summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + boson%Mom
3703 propfac3 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3704 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3705 eps1 = eps1*propfac3
3707 rin = numglu(1)+numglu(2)+
ne+
nf+1
3708 rout= numglu(1)+numglu(2)+numglu(3)+nh
3709 u1 =
cur_f_2f(gluons(rin:rout),quarks(3:3),-quarks(3)%PartType,(/ng+nh,ng,nh/))
3710 pmom3 =
summom(gluons,rin,rout) + quarks(3)%Mom
3711 if( ng.ge.1 .or. nh.ge.1 )
then
3712 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(3)%Mass2)
3713 if( abs(
sc_(pmom3,pmom3) - quarks(3)%Mass2).lt.
propcut ) cycle
3714 if( quarks(3)%PartType.lt.0 )
then
3715 u1 = (-
spi2_(pmom3,u1) + quarks(3)%Mass*u1 )*propfac3
3717 u1 = (+
spb2_(u1,pmom3) + quarks(3)%Mass*u1 )*propfac3
3720 rin = numglu(1)+numglu(2)+numglu(3)+nh+1
3721 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj
3722 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,(/ni+nj,ni,nj/))
3723 pmom4 =
summom(gluons,rin,rout) + quarks(4)%Mom
3724 if( ni.ge.1 .or. nj.ge.1 )
then
3725 propfac4 = (0d0,1d0)/(
sc_(pmom4,pmom4) - quarks(4)%Mass2)
3726 if( abs(
sc_(pmom4,pmom4) - quarks(4)%Mass2).lt.
propcut ) cycle
3727 if( quarks(4)%PartType.lt.0 )
then
3728 ubar2 = (-
spi2_(pmom4,ubar2) + quarks(4)%Mass*ubar2 )*propfac4
3730 ubar2 = (+
spb2_(ubar2,pmom4) + quarks(4)%Mass*ubar2 )*propfac4
3734 if( quarks(4)%PartType.lt.0 )
then
3735 eps2 = +
vbqq(
dv,u1,ubar2)
3737 eps2 = -
vbqq(
dv,ubar2,u1)
3739 tmpmom2 = pmom3 + pmom4
3740 propfac1 = (0d0,-1d0)/
sc_(tmpmom2,tmpmom2)
3741 if( abs(
sc_(tmpmom2,tmpmom2)).lt.
propcut ) cycle
3742 eps2 = eps2*propfac1
3753 tmpgluons(counter)%Mom => tmpmom1(:)
3754 tmpgluons(counter)%Pol => eps1(:)
3755 tmpgluons(counter)%ExtRef => tmpextref1
3757 rin =numglu(1)+numglu(2)+
ne+1
3758 rout=numglu(1)+numglu(2)+
ne+
nf
3763 tmpgluons(counter)%Mom => tmpmom2(:)
3764 tmpgluons(counter)%Pol => eps2(:)
3765 tmpgluons(counter)%ExtRef => tmpextref1
3767 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj+1
3773 eps3(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+nk+2)
3783 if (bosonvertex .eq. 3)
then
3793 rout= numglu(1)+numglu(2)+numglu(3)+nc
3794 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)
3795 pmom2 =
summom(gluons,rin,rout) + quarks(1)%Mom + quarks(2)%Mom + quarks(3)%Mom
3796 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
3797 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
3798 if( quarks(4)%PartType.lt.0 )
then
3799 u1 = (+
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
3801 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
3803 rin = numglu(1)+numglu(2)+numglu(3)+nc+1
3804 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne
3805 ubar2 =
cur_f_2fv(gluons(rin:rout),quarks(4:4),-quarks(4)%PartType,boson,(/nd+
ne,nd,
ne/))
3806 pmom3 =
summom(gluons,rin,rout) + quarks(4)%Mom + boson%Mom
3807 propfac3 = (0d0,1d0)/(
sc_(pmom3,pmom3) - quarks(4)%Mass2)
3808 if( abs(
sc_(pmom3,pmom3) - quarks(4)%Mass2).lt.
propcut ) cycle
3809 if( quarks(4)%PartType.lt.0 )
then
3810 ubar2 = (-
spi2_(pmom3,ubar2) + quarks(4)%Mass*ubar2 )*propfac3
3812 ubar2 = (+
spb2_(ubar2,pmom3) + quarks(4)%Mass*ubar2 )*propfac3
3815 if( quarks(4)%PartType.lt.0 )
then
3816 eps1 = +
vbqq(
dv,u1,ubar2)
3818 eps1 = -
vbqq(
dv,ubar2,u1)
3828 tmpmom1(:) = pmom2(:)+pmom3(:)
3830 tmpgluons(counter)%Mom => tmpmom1(:)
3831 tmpgluons(counter)%Pol => eps1(:)
3832 tmpgluons(counter)%ExtRef => tmpextref1
3834 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+
ne+1
3840 eps2(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+1)
3842 if( na.ge.1 .or.
nf.ge.1 )
then
3843 propfac1 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3844 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3845 eps2 = eps2*propfac1
3854 if (bosonvertex .eq. 3)
then
3859 do nf=0,numglu(3)-
ne
3870 ubar2 =
cur_f_2f(gluons(rin:rout),quarks(1:1),-quarks(1)%PartType,(/nb+nc,nb,nc/))
3871 pmom1 =
summom(gluons,rin,rout) + quarks(1)%Mom
3872 if( nb.ge.1 .or. nc.ge.1 )
then
3873 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1) - quarks(1)%Mass2)
3874 if( abs(
sc_(pmom1,pmom1) - quarks(1)%Mass2).lt.
propcut ) cycle
3875 if( quarks(1)%PartType.lt.0 )
then
3876 ubar2 = (-
spi2_(pmom1,ubar2) + quarks(1)%Mass*ubar2 )*propfac1
3878 ubar2 = (+
spb2_(ubar2,pmom1) + quarks(1)%Mass*ubar2 )*propfac1
3881 rin = numglu(1)+nc+1
3882 rout= numglu(1)+numglu(2)+
ne
3883 u1 =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/nd+
ne,nd,
ne/))
3884 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom
3885 if( nd.ge.1 .or.
ne.ge.1 )
then
3886 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
3887 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
3888 if( quarks(2)%PartType.lt.0 )
then
3889 u1 = (-
spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
3891 u1 = (+
spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
3895 if( quarks(2)%PartType.lt.0 )
then
3896 eps1 = +
vbqq(
dv,ubar2,u1)
3898 eps1 = -
vbqq(
dv,u1,ubar2)
3900 tmpmom1 = pmom1 + pmom2
3901 propfac3 = (0d0,-1d0)/
sc_(tmpmom1,tmpmom1)
3902 if( abs(
sc_(tmpmom1,tmpmom1)).lt.
propcut ) cycle
3903 eps1 = eps1*propfac3
3905 rin = numglu(1)+numglu(2)+
ne+
nf+1
3906 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj
3907 eps2=
cur_g_2fv(gluons(rin:rout),quarks(3:4),boson,(/1+ng+nh+ni+nj,ng,nh,ni,nj /) )
3908 tmpmom2 =
summom(gluons,rin,rout) +quarks(3)%Mom + quarks(4)%Mom + boson%Mom
3909 propfac1 = (0d0,-1d0)/
sc_(tmpmom2,tmpmom2)
3910 if( abs(
sc_(tmpmom2,tmpmom2)).lt.
propcut ) cycle
3911 eps2 = eps2*propfac1
3921 tmpgluons(counter)%Mom => tmpmom1(:)
3922 tmpgluons(counter)%Pol => eps1(:)
3923 tmpgluons(counter)%ExtRef => tmpextref1
3925 rin =numglu(1)+numglu(2)+
ne+1
3926 rout=numglu(1)+numglu(2)+
ne+
nf
3931 tmpgluons(counter)%Mom => tmpmom2(:)
3932 tmpgluons(counter)%Pol => eps2(:)
3933 tmpgluons(counter)%ExtRef => tmpextref1
3935 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+nj+1
3941 eps3(:) =
cur_g(tmpgluons(1:counter-1),1+na+
nf+nk+2)
3954 if (bosonvertex .ne. 1 .and. bosonvertex .ne. 3)
then
3955 print *,
'Not implemented for for cur_g_4fV for this flavor structure and BosonVertex=', bosonvertex
3973 FUNCTION cur_f_6fv(Gluons,Quarks,Quark1PartType,Boson,BosonVertex,NumGlu,tag_f)
result(res)
3975 integer :: numglu(0:6),quark1parttype,tag_f,bosonvertex,bosonvertex_mod
3977 integer,
target :: tmpparttype,tmpextref
3978 complex(8) :: res(1:
ds),tmp(1:
ds)
3980 complex(8) :: u1(1:
ds),ubar1(1:
ds)
3981 complex(8),
target ::
ubar0(1:
ds)
3982 complex(8) :: eps1(1:
dv)
3983 complex(8) :: eps2(1:
dv)
3984 type(
ptrtoparticle) :: tmpgluons(1:numglu(1)+numglu(6)),tmpquark(1:1)
3985 complex(8) :: propfac1,propfac2
3986 complex(8),
target :: pmom1(1:
dv)
3987 complex(8),
target :: pmom2(1:
dv)
3988 integer :: n1a,n1b,n2a,n2b,n3a,n3b,n4a,n4b,n5a,n5b,n6a,n6b
3989 integer :: rin,rout,i,counter
3993 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"
4000 if (quark1parttype.eq.-quarks(2)%PartType)
then
4001 if ( ((quarks(3)%PartType.eq.-quarks(4)%PartType) .AND. (bosonvertex .eq. 2 .OR.bosonvertex .eq. 4 .OR.bosonvertex .eq. 6)) .OR. &
4002 ((quarks(3)%PartType.eq.-quarks(6)%PartType) .AND. (bosonvertex .eq. 2 .OR.bosonvertex .eq. 6)) )
then
4003 print *,
'WARNING : cur_f_6fV with this flavor structure and this choice of BosonVertex is implemented, but not checked!'
4004 print *,
'Quark flavors: ',quark1parttype, quarks(2)%PartType,quarks(3)%PartType,quarks(4)%PartType,quarks(5)%PartType,quarks(6)%PartType
4005 print *,
'BosonVertex =', bosonvertex
4007 elseif (quark1parttype.eq.-quarks(6)%PartType )
then
4008 if ( ((quark1parttype.eq.-quarks(2)%PartType) .AND. (bosonvertex .eq. 2 .or. bosonvertex .eq. 4 .or. bosonvertex .eq. 6)) .OR. &
4009 & ((quark1parttype.eq.-quarks(4)%PartType) .AND. (bosonvertex .eq. 4 .or. bosonvertex .eq. 6)) )
then
4010 print *,
'WARNING : cur_f_6fV with this flavor structure and this choice of BosonVertex is implemented, but not checked!'
4011 print *,
'Quark flavors: ',quark1parttype, quarks(2)%PartType,quarks(3)%PartType,quarks(4)%PartType,quarks(5)%PartType,quarks(6)%PartType
4012 print *,
'BosonVertex =', bosonvertex
4014 elseif (quarks(2)%PartType.eq.-quarks(3)%PartType)
then
4015 if ( ((quark1parttype.eq.-quarks(4)%PartType) .AND. (bosonvertex .eq. 4 .or. bosonvertex .eq. 6) ) .OR. &
4016 ((quark1parttype.eq.-quarks(6)%PartType) .AND. (bosonvertex .eq. 4 .or. bosonvertex .eq. 6) ) )
then
4017 print *,
'WARNING : cur_f_6fV with this flavor structure and this choice of BosonVertex is implemented, but not checked!'
4018 print *,
'Quark flavors: ',quark1parttype, quarks(2)%PartType,quarks(3)%PartType,quarks(4)%PartType,quarks(5)%PartType,quarks(6)%PartType
4019 print *,
'BosonVertex =', bosonvertex
4028 if( quark1parttype.eq.-quarks(2)%PartType .AND. (quarks(3)%PartType.eq.-quarks(4)%PartType .or. quarks(3)%PartType.eq.-quarks(6)%PartType) &
4029 .AND. (quarks(2)%ExtRef.ne.-1 .or. tag_f.ne.1) &
4030 .AND. ( bosonvertex.eq.1 .OR. bosonvertex.eq.2) &
4038 rin =numglu(1)+n2a+1
4039 rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4040 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/))
4041 pmom1(:) =
summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
4042 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4044 eps2 = eps2*propfac1
4050 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 2)
then
4051 ubar1(:) =
cur_f_2fv(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,boson,(/n1b+n2a,n1b,n2a/) )
4052 pmom2(:) = quarks(2)%Mom +
summom(gluons,rin,rout) + boson%Mom
4053 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
4054 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut )
then
4058 if( quarks(2)%PartType.lt.0 )
then
4059 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
4061 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
4063 if( quarks(2)%PartType.lt.0 )
then
4070 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4071 pmom1 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
4072 if(n1a.ge.1 .or. n6b.ge.1)
then
4073 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
4074 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut )
then
4078 if( quarks(2)%PartType.lt.0 )
then
4085 tmpquark(1)%Mom => pmom1(:)
4086 tmpquark(1)%Pol =>
ubar0(:)
4087 tmpquark(1)%Mass => quarks(2)%Mass
4088 tmpquark(1)%Mass2=> quarks(2)%Mass2
4090 tmpquark(1)%ExtRef => tmpextref
4091 tmpquark(1)%PartType => quarks(2)%PartType
4099 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4105 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n6b/) )
4106 res(:) = res(:) + tmp(:)
4108 if (bosonvertex .eq. 1 .or. bosonvertex .eq. 6)
then
4112 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n1b+n2a,n1b,n2a/) )
4113 if(n1b.ge.1 .or. n2a.ge.1)
then
4114 pmom2(:) = quarks(2)%Mom +
summom(gluons,rin,rout)
4115 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
4116 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut ) cycle
4117 if( quarks(2)%PartType.lt.0 )
then
4118 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
4120 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
4123 if( quarks(2)%PartType.lt.0 )
then
4130 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4131 pmom1 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
4133 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
4134 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut ) cycle
4135 if( quarks(2)%PartType.lt.0 )
then
4141 tmpquark(1)%Mom => pmom1(:)
4142 tmpquark(1)%Pol =>
ubar0(:)
4143 tmpquark(1)%Mass => quarks(2)%Mass
4144 tmpquark(1)%Mass2=> quarks(2)%Mass2
4146 tmpquark(1)%ExtRef => tmpextref
4147 tmpquark(1)%PartType => quarks(2)%PartType
4155 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4161 tmp(:) =
cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,boson,(/counter-1,n1a,n6b/) )
4162 res(:) = res(:) + tmp(:)
4172 if( quark1parttype.eq.-quarks(2)%PartType .AND. (quarks(3)%PartType.eq.-quarks(4)%PartType .or. quarks(3)%PartType.eq.-quarks(6)%PartType) &
4173 .AND. (quarks(2)%ExtRef.ne.-1 .or. tag_f.ne.1) &
4174 .AND. ( bosonvertex.eq.3 .OR. bosonvertex .eq. 4 .OR. bosonvertex.eq.5 ) &
4182 rin =numglu(1)+n2a+1
4183 rout=numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4184 bosonvertex_mod = bosonvertex - 2
4185 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/))
4186 pmom1(:) =
summom(gluons,rin,rout) + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
4187 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4189 eps2 = eps2*propfac1
4195 ubar1(:) =
cur_f_2f(gluons(rin:rout),quarks(2:2),-quarks(2)%PartType,(/n1b+n2a,n1b,n2a/) )
4196 if(n1b.ge.1 .or. n2a.ge.1)
then
4197 pmom2(:) = quarks(2)%Mom +
summom(gluons,rin,rout)
4198 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2)-quarks(2)%Mass2)
4199 if( abs(
sc_(pmom2,pmom2)-quarks(2)%Mass2).lt.
propcut ) cycle
4200 if( quarks(2)%PartType.lt.0 )
then
4201 ubar1(:) = (-
spi2_(pmom2,ubar1)+quarks(2)%Mass*ubar1(:))*propfac2
4203 ubar1(:) = (+
spb2_(ubar1,pmom2)+quarks(2)%Mass*ubar1(:))*propfac2
4206 if( quarks(2)%PartType.lt.0 )
then
4213 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4214 pmom1 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
4215 if(n1a.ge.1 .or. n6b.ge.1)
then
4216 propfac1 = (0d0,1d0)/(
sc_(pmom1,pmom1)-quarks(2)%Mass2)
4217 if( abs(
sc_(pmom1,pmom1)-quarks(2)%Mass2).lt.
propcut ) cycle
4218 if( quarks(2)%PartType.lt.0 )
then
4225 tmpquark(1)%Mom => pmom1(:)
4226 tmpquark(1)%Pol =>
ubar0(:)
4227 tmpquark(1)%Mass => quarks(2)%Mass
4228 tmpquark(1)%Mass2=> quarks(2)%Mass2
4230 tmpquark(1)%ExtRef => tmpextref
4231 tmpquark(1)%PartType => quarks(2)%PartType
4239 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4245 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),-tmpquark(1)%PartType,(/counter-1,n1a,n6b/) )
4246 res(:) = res(:) + tmp(:)
4267 if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
4268 ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
4269 .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
4270 .AND. ( bosonvertex.eq.1 .OR. bosonvertex.eq.6 ) &
4280 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4281 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4282 eps2 =
cur_g_2f(gluons(rin:rout),quarks(5:6),(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
4283 pmom1(:) =
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4284 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4286 eps2 = eps2*propfac1
4289 rout= numglu(1)+numglu(2)+numglu(3)+n4a
4290 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)
4291 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
4293 if( quark1parttype.eq.-quarks(2)%PartType)
then
4294 if( quarks(2)%PartType.lt.0 )
then
4295 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4296 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4297 u1 = (-
spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
4299 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4300 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4301 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4303 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4304 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4305 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4308 elseif( quarks(2)%PartType.gt.0 )
then
4309 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4310 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4311 u1 = (
spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
4313 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4314 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4315 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4317 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4318 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4319 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4323 tmpquark(1)%Mom => pmom2(:)
4324 tmpquark(1)%Pol =>
ubar0(:)
4325 tmpquark(1)%Mass => quarks(2)%Mass
4326 tmpquark(1)%Mass2=> quarks(2)%Mass2
4327 elseif( quark1parttype.eq.-quarks(4)%PartType)
then
4328 if( quarks(4)%PartType.lt.0 )
then
4329 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4330 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4331 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
4333 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4334 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4335 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4336 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4337 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4338 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4341 elseif( quarks(4)%PartType.gt.0 )
then
4342 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4343 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4344 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
4346 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4347 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4348 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4349 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4350 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4351 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4355 tmpquark(1)%Mom => pmom2(:)
4356 tmpquark(1)%Pol =>
ubar0(:)
4357 tmpquark(1)%Mass => quarks(4)%Mass
4358 tmpquark(1)%Mass2=> quarks(4)%Mass2
4361 tmpquark(1)%ExtRef => tmpextref
4362 tmpparttype = -quark1parttype
4363 tmpquark(1)%PartType => tmpparttype
4371 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4377 tmp(:) =
cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,boson,(/counter-1,n1a,n6b/) )
4379 res(:) = res(:) + tmp(:)
4388 if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
4389 ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
4390 .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
4391 .AND. ( bosonvertex.eq.1 .OR. bosonvertex.eq.2 .OR. bosonvertex.eq.3 .OR. bosonvertex.eq.4 ) &
4401 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4402 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4403 eps2 =
cur_g_2f(gluons(rin:rout),quarks(5:6),(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
4404 pmom1(:) =
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4405 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4407 eps2 = eps2*propfac1
4410 rout= numglu(1)+numglu(2)+numglu(3)+n4a
4411 bosonvertex_mod = bosonvertex
4412 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)
4413 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom + boson%Mom
4415 if( quark1parttype.eq.-quarks(2)%PartType)
then
4416 if( quarks(2)%PartType.lt.0 )
then
4417 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4418 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4419 u1 = (-
spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
4421 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4422 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4423 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4424 if( n1a.ge.1 .or. n6b.ge.1 )
then
4425 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4426 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4429 elseif( quarks(2)%PartType.gt.0 )
then
4430 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4431 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4432 u1 = (
spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
4434 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4435 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4436 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4437 if( n1a.ge.1 .or. n6b.ge.1 )
then
4438 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4439 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4443 tmpquark(1)%Mom => pmom2(:)
4444 tmpquark(1)%Pol =>
ubar0(:)
4445 tmpquark(1)%Mass => quarks(2)%Mass
4446 tmpquark(1)%Mass2=> quarks(2)%Mass2
4447 elseif( quark1parttype.eq.-quarks(4)%PartType)
then
4448 if( quarks(4)%PartType.lt.0 )
then
4449 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4450 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4451 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
4453 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4454 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4455 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4456 if( n1a.ge.1 .or. n6b.ge.1 )
then
4457 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4458 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4461 elseif( quarks(4)%PartType.gt.0 )
then
4462 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4463 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4464 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
4466 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4467 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4468 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4469 if( n1a.ge.1 .or. n6b.ge.1 )
then
4470 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4471 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4475 tmpquark(1)%Mom => pmom2(:)
4476 tmpquark(1)%Pol =>
ubar0(:)
4477 tmpquark(1)%Mass => quarks(4)%Mass
4478 tmpquark(1)%Mass2=> quarks(4)%Mass2
4481 tmpquark(1)%ExtRef => tmpextref
4482 tmpparttype = -quark1parttype
4483 tmpquark(1)%PartType => tmpparttype
4491 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4497 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
4499 res(:) = res(:) + tmp(:)
4508 if( quarks(5)%PartType.eq.-quarks(6)%PartType .AND. &
4509 ((quark1parttype.eq.-quarks(2)%PartType .and. (quarks(2)%ExtRef.ne.-1.or.tag_f.ne.1) ) &
4510 .OR. (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1) ))&
4511 .AND. ( bosonvertex.eq.5 ) &
4521 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4522 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4523 eps2 =
cur_g_2fv(gluons(rin:rout),quarks(5:6),boson,(/1+n4b+numglu(5)+n6a,n4b,numglu(5),n6a/))
4524 pmom1(:) =
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
4525 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4527 eps2 = eps2*propfac1
4530 rout= numglu(1)+numglu(2)+numglu(3)+n4a
4531 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)
4532 pmom2 =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + quarks(4)%Mom
4534 if( quark1parttype.eq.-quarks(2)%PartType)
then
4536 if( quarks(2)%PartType.lt.0 )
then
4537 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4538 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4539 u1 = (-
spi2_(pmom2,u1) + quarks(2)%Mass*u1 )*propfac2
4541 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4542 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4543 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
4544 if( n1a.ge.1 .or. n6b.ge.1 )
then
4545 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4546 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4549 elseif( quarks(2)%PartType.gt.0 )
then
4550 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4551 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4552 u1 = (
spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
4554 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4555 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4556 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom + boson%Mom
4557 if( n1a.ge.1 .or. n6b.ge.1 )
then
4558 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(2)%Mass2)
4559 if( abs(
sc_(pmom2,pmom2) - quarks(2)%Mass2).lt.
propcut ) cycle
4560 u1 = (
spb2_(u1,pmom2) + quarks(2)%Mass*u1 )*propfac2
4563 tmpquark(1)%Mom => pmom2(:)
4564 tmpquark(1)%Pol =>
ubar0(:)
4565 tmpquark(1)%Mass => quarks(2)%Mass
4566 tmpquark(1)%Mass2=> quarks(2)%Mass2
4567 elseif( quark1parttype.eq.-quarks(4)%PartType)
then
4568 if( quarks(4)%PartType.lt.0 )
then
4569 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4570 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4571 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
4573 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4574 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4575 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4576 if( n1a.ge.1 .or. n6b.ge.1 )
then
4577 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4578 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4581 elseif( quarks(4)%PartType.gt.0 )
then
4582 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4583 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4584 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
4586 rin = numglu(1)+numglu(2)+numglu(3)+n4a+1
4587 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4588 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(5)%Mom + quarks(6)%Mom
4589 if( n1a.ge.1 .or. n6b.ge.1 )
then
4590 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4591 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4595 tmpquark(1)%Mom => pmom2(:)
4596 tmpquark(1)%Pol =>
ubar0(:)
4597 tmpquark(1)%Mass => quarks(4)%Mass
4598 tmpquark(1)%Mass2=> quarks(4)%Mass2
4601 tmpquark(1)%ExtRef => tmpextref
4602 tmpparttype = -quark1parttype
4603 tmpquark(1)%PartType => tmpparttype
4611 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4617 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/) )
4619 res(:) = res(:) + tmp(:)
4630 if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
4631 (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
4632 .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1)) ) &
4633 .AND. ( bosonvertex.eq.1 ) &
4645 rout= numglu(1)+numglu(2)+n3a
4646 eps2 =
cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
4647 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4648 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4650 eps2 = eps2*propfac1
4652 rin = numglu(1)+numglu(2)+n3a+1
4653 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4654 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)
4655 pmom2 =
summom(gluons,rin,rout) + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
4657 if( quark1parttype.eq.-quarks(4)%PartType )
then
4658 if( quarks(4)%PartType.lt.0 )
then
4659 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4660 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4661 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
4664 rout= numglu(1)+numglu(2)+n3a
4665 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4666 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4667 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4668 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4671 elseif( quarks(4)%PartType.gt.0 )
then
4672 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4673 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4674 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
4677 rout= numglu(1)+numglu(2)+n3a
4678 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4679 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4680 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4681 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4685 tmpquark(1)%Mom => pmom2(:)
4686 tmpquark(1)%Pol =>
ubar0(:)
4687 tmpquark(1)%Mass => quarks(4)%Mass
4688 tmpquark(1)%Mass2=> quarks(4)%Mass2
4689 elseif(quark1parttype.eq.-quarks(6)%PartType)
then
4690 if( quarks(6)%PartType.lt.0 )
then
4691 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4692 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4693 u1 = (-
spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
4696 rout= numglu(1)+numglu(2)+n3a
4697 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4698 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4699 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4700 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4703 elseif( quarks(6)%PartType.gt.0 )
then
4704 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4705 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4706 u1 = (
spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
4709 rout= numglu(1)+numglu(2)+n3a
4710 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4711 if( n1a.ge.1 .or. n6b.ge.1 .or. bosonvertex.eq.1 .or. bosonvertex.eq.6 )
then
4712 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4713 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4719 tmpquark(1)%ExtRef => tmpextref
4720 tmpparttype = -quark1parttype
4721 tmpquark(1)%PartType => tmpparttype
4729 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4735 tmp(:) =
cur_f_2fv(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,boson,(/counter-1,n1a,n6b/))
4737 res(:) = res(:) + tmp(:)
4747 if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
4748 (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
4749 .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1)) ) &
4750 .AND. ( bosonvertex.eq.2) &
4762 rout= numglu(1)+numglu(2)+n3a
4763 eps2 =
cur_g_2fv(gluons(rin:rout),quarks(2:3),boson,(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
4764 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
4765 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4767 eps2 = eps2*propfac1
4769 rin = numglu(1)+numglu(2)+n3a+1
4770 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4771 bosonvertex_mod = bosonvertex
4772 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)
4773 pmom2 =
summom(gluons,rin,rout) + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom
4775 if( quark1parttype.eq.-quarks(4)%PartType )
then
4776 if( quarks(4)%PartType.lt.0 )
then
4777 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4778 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4779 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
4782 rout= numglu(1)+numglu(2)+n3a
4783 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
4784 if( n1a.ge.1 .or. n6b.ge.1 )
then
4785 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4786 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4789 elseif( quarks(4)%PartType.gt.0 )
then
4790 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4791 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4792 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
4795 rout= numglu(1)+numglu(2)+n3a
4796 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
4797 if( n1a.ge.1 .or. n6b.ge.1 )
then
4798 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4799 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4803 elseif(quark1parttype.eq.-quarks(6)%PartType)
then
4804 if( quarks(6)%PartType.lt.0 )
then
4805 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4806 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4807 u1 = (-
spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
4810 rout= numglu(1)+numglu(2)+n3a
4811 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom+ boson%Mom
4812 if( n1a.ge.1 .or. n6b.ge.1 )
then
4813 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4814 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4817 elseif( quarks(6)%PartType.gt.0 )
then
4818 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4819 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4820 u1 = (
spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
4823 rout= numglu(1)+numglu(2)+n3a
4824 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom + boson%Mom
4825 if( n1a.ge.1 .or. n6b.ge.1 )
then
4826 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4827 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4832 tmpquark(1)%Mom => pmom2(:)
4833 tmpquark(1)%Pol =>
ubar0(:)
4834 tmpquark(1)%Mass => quarks(4)%Mass
4835 tmpquark(1)%Mass2=> quarks(4)%Mass2
4837 tmpquark(1)%ExtRef => tmpextref
4838 tmpparttype = -quark1parttype
4839 tmpquark(1)%PartType => tmpparttype
4847 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4853 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/))
4855 res(:) = res(:) + tmp(:)
4863 if( quarks(2)%PartType.eq.-quarks(3)%PartType .AND. ( &
4864 (quark1parttype.eq.-quarks(4)%PartType .and. (quarks(4)%ExtRef.ne.-1.or.tag_f.ne.1)) &
4865 .OR. (quark1parttype.eq.-quarks(6)%PartType .and. (quarks(6)%ExtRef.ne.-1.or.tag_f.ne.1)) ) &
4866 .AND. ( bosonvertex .ge. 3 .and. bosonvertex .le. 6 ) )
then
4877 rout= numglu(1)+numglu(2)+n3a
4878 eps2 =
cur_g_2f(gluons(rin:rout),quarks(2:3),(/1+n1b+numglu(2)+n3a,n1b,numglu(2),n3a/))
4879 pmom1(:) =
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4880 propfac1 = (0d0,-1d0)/
sc_(pmom1,pmom1)
4882 eps2 = eps2*propfac1
4884 rin = numglu(1)+numglu(2)+n3a+1
4885 rout= numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a
4886 bosonvertex_mod=bosonvertex-2
4887 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)
4888 pmom2 =
summom(gluons,rin,rout) + quarks(4)%Mom + quarks(5)%Mom + quarks(6)%Mom+boson%Mom
4890 if( quark1parttype.eq.-quarks(4)%PartType )
then
4891 if( quarks(4)%PartType.lt.0 )
then
4892 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4893 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4894 u1 = (-
spi2_(pmom2,u1) + quarks(4)%Mass*u1 )*propfac2
4897 rout= numglu(1)+numglu(2)+n3a
4898 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4899 if( n1a.ge.1 .or. n6b.ge.1 )
then
4900 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4901 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4904 elseif( quarks(4)%PartType.gt.0 )
then
4905 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4906 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4907 u1 = (
spb2_(u1,pmom2) + quarks(4)%Mass*u1 )*propfac2
4910 rout= numglu(1)+numglu(2)+n3a
4911 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4912 if( n1a.ge.1 .or. n6b.ge.1 )
then
4913 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(4)%Mass2)
4914 if( abs(
sc_(pmom2,pmom2) - quarks(4)%Mass2).lt.
propcut ) cycle
4918 tmpquark(1)%Mom => pmom2(:)
4919 tmpquark(1)%Pol =>
ubar0(:)
4920 tmpquark(1)%Mass => quarks(4)%Mass
4921 tmpquark(1)%Mass2=> quarks(4)%Mass2
4922 elseif(quark1parttype.eq.-quarks(6)%PartType)
then
4923 if( quarks(6)%PartType.lt.0 )
then
4924 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4925 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4926 u1 = (-
spi2_(pmom2,u1) + quarks(6)%Mass*u1 )*propfac2
4929 rout= numglu(1)+numglu(2)+n3a
4930 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4931 if( n1a.ge.1 .or. n6b.ge.1 )
then
4932 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4933 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4936 elseif( quarks(6)%PartType.gt.0 )
then
4937 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4938 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4939 u1 = (
spb2_(u1,pmom2) + quarks(6)%Mass*u1 )*propfac2
4942 rout= numglu(1)+numglu(2)+n3a
4943 pmom2 = pmom2 +
summom(gluons,rin,rout) + quarks(2)%Mom + quarks(3)%Mom
4944 if( n1a.ge.1 .or. n6b.ge.1 )
then
4945 propfac2 = (0d0,1d0)/(
sc_(pmom2,pmom2) - quarks(6)%Mass2)
4946 if( abs(
sc_(pmom2,pmom2) - quarks(6)%Mass2).lt.
propcut ) cycle
4950 tmpquark(1)%Mom => pmom2(:)
4951 tmpquark(1)%Pol =>
ubar0(:)
4952 tmpquark(1)%Mass => quarks(6)%Mass
4953 tmpquark(1)%Mass2=> quarks(6)%Mass2
4956 tmpquark(1)%ExtRef => tmpextref
4957 tmpparttype = -quark1parttype
4958 tmpquark(1)%PartType => tmpparttype
4966 rin =numglu(1)+numglu(2)+numglu(3)+numglu(4)+numglu(5)+n6a+1
4972 tmp(:) =
cur_f_2f(tmpgluons(1:counter-1),tmpquark(1:1),quark1parttype,(/counter-1,n1a,n6b/))
4974 res(:) = res(:) + tmp(:)
4992 FUNCTION summom(Particles,i1,i2)
5010 complex(8),
intent(in) :: p1(1:
dv),p2(1:
dv)
5025 complex(8) :: k1(1:
dv),k2(1:
dv),v1(1:
dv),v2(1:
dv)
5026 complex(8),
parameter :: ioversqrt2=(0d0,1d0)/dsqrt(2d0)
5029 - 2d0*v1(1:
dv)*(k1.ndot.v2) &
5030 + 2d0*v2(1:
dv)*(k2.ndot.v1) )
5037 complex(8) :: k1(1:
dv),k2(1:
dv),k3(1:
dv)
5038 complex(8),
parameter :: i=(0d0,1d0)
5041 + k2(1:
dv)*(k1.ndot.k3) &
5042 - k3(1:
dv)*(k1.ndot.k2)*0.5d0 )
5054 linear_map = i2+ngluons*(i2-i1)-((i2-i1)*(i2-i1+1))/2
5068 outpointer%PartType => inpointer%PartType
5069 outpointer%ExtRef => inpointer%ExtRef
5070 outpointer%Mass => inpointer%Mass
5071 outpointer%Mass2 => inpointer%Mass2
5072 outpointer%Helicity => inpointer%Helicity
5073 outpointer%Mom => inpointer%Mom
5074 outpointer%Pol => inpointer%Pol
5092 subroutine spi2(Dv,Ds,v,sp,f)
5094 integer i,i1,i2,i3,imax,Dv,Ds
5095 double complex sp(Ds),v(Dv),f(Ds)
5096 double complex x0(4,4),xx(4,4),xy(4,4)
5097 double complex xz(4,4),x5(4,4)
5098 double complex y1,y2,y3,y4,bp,bm,cp,cm
5099 double complex test(Ds)
5125 xy(1,i)=dcmplx(0d0,-1d0)*y4
5126 xy(2,i)=dcmplx(0d0,1d0)*y3
5127 xy(3,i)=dcmplx(0d0,1d0)*y2
5128 xy(4,i)=dcmplx(0d0,-1d0)*y1
5146 f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)
5153 bp = (v(5)+dcmplx(0d0,1d0)*v(6))
5154 bm=(v(5)-dcmplx(0d0,1d0)*v(6))
5159 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)
5163 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)
5171 bp = (v(5)+dcmplx(0d0,1d0)*v(6))
5172 bm=(v(5)-dcmplx(0d0,1d0)*v(6))
5173 cp=(v(7)+dcmplx(0d0,1d0)*v(8))
5174 cm=(v(7)-dcmplx(0d0,1d0)*v(8))
5180 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)
5184 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)
5188 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)
5192 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)
5213 function spi2_(v,sp)
5215 double complex,
intent(in) :: sp(:),v(:)
5216 double complex ::
spi2_(size(sp)) ,tmp(size(sp))
5222 if (
ds == 16)
dv = 8
5234 function spivl(sp,v)
5236 double complex :: sp(:),v(:)
5237 double complex ::
spivl(size(sp))
5239 spivl(1) = sp(1)*v(1) + sp(4)*(v(2) + (0d0,1d0)*v(3)) + sp(3)*v(4)
5240 spivl(2) = sp(2)*v(1) + sp(3)*(v(2) - (0d0,1d0)*v(3)) - sp(4)*v(4)
5241 spivl(3) =-sp(3)*v(1) - sp(2)*(v(2) + (0d0,1d0)*v(3)) - sp(1)*v(4)
5242 spivl(4) =-sp(4)*v(1) - sp(1)*(v(2) - (0d0,1d0)*v(3)) + sp(2)*v(4)
5249 subroutine spb2(Dv,Ds,sp,v,f)
5251 integer i,i1,i2,i3,Dv,Ds,imax
5252 double complex sp(Ds),v(Dv),f(Ds)
5253 double complex x0(4,4),xx(4,4),xy(4,4)
5254 double complex xz(4,4),x5(4,4)
5255 double complex y1,y2,y3,y4,bp,bm,cp,cm
5256 double complex test(Ds)
5279 xy(1,i)=dcmplx(0d0,-1d0)*y4
5280 xy(2,i)=dcmplx(0d0,1d0)*y3
5281 xy(3,i)=dcmplx(0d0,1d0)*y2
5282 xy(4,i)=dcmplx(0d0,-1d0)*y1
5300 f(i)=v(1)*x0(i,1)-v(2)*xx(i,1)-v(3)*xy(i,1)-v(4)*xz(i,1)
5307 bp = (v(5)+dcmplx(0d0,1d0)*v(6))
5308 bm=(v(5)-dcmplx(0d0,1d0)*v(6))
5312 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)
5316 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)
5324 bp=(v(5)+dcmplx(0d0,1d0)*v(6))
5325 bm=(v(5)-dcmplx(0d0,1d0)*v(6))
5326 cp=(v(7)+dcmplx(0d0,1d0)*v(8))
5327 cm=(v(7)-dcmplx(0d0,1d0)*v(8))
5331 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)
5335 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)
5339 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)
5343 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)
5367 function spb2_(sp,v)
5369 double complex,
intent(in) :: sp(:),v(:)
5370 double complex ::
spb2_(size(sp)) ,tmp(size(sp))
5376 if (
ds == 16)
dv = 8
5388 function psp1_(sp1,sp2)
result(res)
5390 complex(8),
intent(in) :: sp1(:)
5391 complex(8),
intent(in) :: sp2(:)
5394 res = sum(sp1(1:)*sp2(1:))
5402 recursive function chir(sign,sp)
result(res)
5405 double complex :: sp(:)
5406 double complex :: res(size(sp))
5412 res(1) = 0.5d0*(sp(1)+sp(3))
5413 res(2) = 0.5d0*(sp(2)+sp(4))
5417 res(1) = 0.5d0*(sp(1)-sp(3))
5418 res(2) = 0.5d0*(sp(2)-sp(4))
5423 res(1:d/2) =
chir(sign,sp(1:d/2))
5424 res((d/2+1):d) =
chir(sign,sp( (d/2+1):d ))
5432 recursive function ichir(sign,sp)
result(res)
5437 double complex :: sp(:)
5438 double complex :: res(size(sp))
5439 double complex ::
ci
5446 res(1) = 0.5d0*(sp(1)+
ci*sp(3))
5447 res(2) = 0.5d0*(sp(2)+
ci*sp(4))
5448 res(3) = 0.5d0*(
ci*sp(1)+sp(3))
5449 res(4) = 0.5d0*(
ci*sp(2)+sp(4))
5451 res(1) = 0.5d0*(sp(1)-
ci*sp(3))
5452 res(2) = 0.5d0*(sp(2)-
ci*sp(4))
5453 res(3) = 0.5d0*(-
ci*sp(1)+sp(3))
5454 res(4) = 0.5d0*(-
ci*sp(2)+sp(4))
5457 res(1:d/2) =
ichir(sign,sp(1:d/2))
5458 res((d/2+1):d) =
ichir(sign,sp( (d/2+1):d ))
5466 function vbqv(sp,e1,coupl_left,coupl_right)
5468 complex(8),
intent(in) :: e1(:)
5469 complex(8),
intent(in) :: sp(:)
5470 complex(8),
intent(in) :: coupl_left,coupl_right
5471 complex(8) ::
vbqv(size(sp))
5474 vbqv = -(0d0,1d0)*( coupl_left*
chir(.false.,sp) + coupl_right*
chir(.true.,sp) )
5481 function vvq(e1,sp,coupl_left,coupl_right)
5483 complex(8),
intent(in) :: e1(:)
5484 complex(8),
intent(in) :: sp(:)
5485 complex(8),
intent(in) :: coupl_left,coupl_right
5486 complex(8) ::
vvq(size(sp))
5489 vvq = -(0d0,1d0)*( coupl_left*
chir(.false.,sp) + coupl_right*
chir(.true.,sp) )
5500 complex(8) :: p1(:),p2(:)
5504 sizemin=min(
size(p1),
size(p2))
5512 subroutine rsc_(n,x,y,r)
5515 complex(8) x(*),y(*)
5531 function vggg(e1,k1,e2,k2)
5533 complex(8),
intent(in) :: e1(:), e2(:)
5534 complex(8),
intent(in) :: k1(:), k2(:)
5535 complex(8) ::
vggg(size(e1))
5536 complex(8):: sk1e2,se1e2,sk2e1,xx
5537 real(8),
parameter ::
sqrt2 = 1.4142135623730950488016887242096980786d0
5543 xx=(0.0d0,1.0d0)*
sqrt2
5544 vggg = xx*(-sk1e2*e1+sk2e1*e2+se1e2/2d0*(k1-k2))
5549 function vgggg(e1,e2,e3)
5551 complex(8),
intent(in) :: e1(:),e2(:),e3(:)
5552 complex(8) ::
vgggg(size(e1))
5553 complex(8):: se1e3,se2e3,se1e2
5558 vgggg = (0.0d0,1.0d0)*(e2*se1e3-0.5d0*(e1*se2e3+ e3*se1e2))
5565 complex(8),
intent(in) :: e1(:)
5566 complex(8),
intent(in) :: sp(:)
5567 complex(8) ::
vqg(size(sp))
5568 real(8),
parameter ::
sqrt2 = 1.4142135623730950488016887242096980786d0
5579 complex(8),
intent(in) :: e1(:)
5580 complex(8),
intent(in) :: sp(:)
5581 complex(8) ::
vgq(size(sp))
5582 real(8),
parameter ::
sqrt2 = 1.4142135623730950488016887242096980786d0
5592 function vbqg(sp,e1)
5594 complex(8),
intent(in) :: e1(:)
5595 complex(8),
intent(in) :: sp(:)
5596 complex(8) ::
vbqg(size(sp))
5597 real(8),
parameter ::
sqrt2 = 1.4142135623730950488016887242096980786d0
5605 function vgbq(e1,sp)
5607 complex(8),
intent(in) :: e1(:)
5608 complex(8),
intent(in) :: sp(:)
5609 complex(8) ::
vgbq(size(sp))
5610 real(8),
parameter ::
sqrt2 = 1.4142135623730950488016887242096980786d0
5653 function vbqq(Dummy,sp1,sp2)
5655 complex(8),
intent(in) :: sp1(1:4), sp2(1:4)
5657 complex(8) ::
vbqq(4)
5658 complex(8) :: sp1a(4)
5659 real(8) :: va(1:4,1:4)
5660 real(8),
parameter ::
sqrt2 = 1.4142135623730950488016887242096980786d0
5662 va(1,1:4)=(/+1d0,0d0,0d0,0d0/)
5663 va(2,1:4)=(/0d0,-1d0,0d0,0d0/)
5664 va(3,1:4)=(/0d0,0d0,-1d0,0d0/)
5665 va(4,1:4)=(/0d0,0d0,0d0,-1d0/)
5668 sp1a=
spivl(sp1,dcmplx(va(i,1:4)))
5669 vbqq(i) = (sp1a(1)*sp2(1)+sp1a(2)*sp2(2)+sp1a(3)*sp2(3)+sp1a(4)*sp2(4)) * (0d0,-1d0)/
sqrt2
5678 function vvbqq(sp1,sp2)
5680 complex(8),
intent(in) :: sp1(:), sp2(:)
5682 complex(8) ::
vvbqq(4,4)
5683 complex(8) :: sp1a(4)
5684 real(8) :: va(1:4,1:4)
5686 va(1,1:4)=(/+1d0,0d0,0d0,0d0/)
5687 va(2,1:4)=(/0d0,-1d0,0d0,0d0/)
5688 va(3,1:4)=(/0d0,0d0,-1d0,0d0/)
5689 va(4,1:4)=(/0d0,0d0,0d0,-1d0/)
5693 sp1a=
spb2_(sp1, dcmplx(va(i,1:4)))
5694 sp1a=
spb2_(sp1a,dcmplx(va(j,1:4)))