[u/mrichter/AliRoot.git] / PYTHIA / pythia / rkbbv.F

 
C*********************************************************************
 
      SUBROUTINE RKBBV(AK1,AK2,AP1,AP2,ALEP1,ALEP2,IMC,RESULT)
 
C...The following routines have been written by Ronald Kleiss,
C...to evaluate the matrix element for g + g -> Z + q + qbar,
C...with massive quarks (e.g. q = b).
C...They have been modified, so that all routines and commonblocks
C...have names beginning with RK, and so that some unnecessary
C...initialization information is not printed. Further, COMPLEX*16
C...has been changed to COMPLEX and REAL*8 to DOUBLE PRECISION
C...(in a few cases to REAL), so as to make the program better
C...transportable.
 
* THE CROSS SECTION FOR
* G(K1) + G(K2) ---> Z(QV) + B(P1) + B_BAR(P2)
*                     |
*                     +---> L(LEP1) + LEP_BAR(LEP2)
* THE B QUARKS HAVE TO BE ON-SHELL, THE LEPTONS MASSLESS
* THE OPTION IMC=0 PERFORMS THE STANDARD SPIN SUM
* THE OPTION IMC=1 PERFORMS THE CALCULATION FOR 'NMC' RANDOMLY
* CHOSEN HELICITY STATES WHICH IMPROVES THE
* SPEED BY A FACTOR 32/NMC
      SAVE
 
      REAL AK1(0:3),AK2(0:3),AP1(0:3),AP2(0:3),ALEP1(0:3),ALEP2(0:3)
      DOUBLE PRECISION K1(0:4),K2(0:4),P1(0:4),P2(0:4),LEP1(0:4),
     &LEP2(0:4)
      REAL RMQ,RMV,RGV,GSTR,VB,AB,VL,AL
      INTEGER INIT
      INTEGER J1,J2,J3,J4,J5
      INTEGER K,IMC,KLOW,KUPP,NMC,OLDIMC
      DOUBLE PRECISION RKRAND,RKDOT,MULT,RMB
C      INTEGER CHKGL1,CHKGL2
      DOUBLE PRECISION QV(0:4),R1(0:4),R2(0:4),Q1(0:4),Q2(0:4)
      DOUBLE PRECISION PP2(0:4)
      DOUBLE PRECISION CROSS
      INTEGER LG1,LG2,LV,L1,L2,HELIX,HELI
      COMPLEX ZFACV,ZFAC1,ZFAC2
      DOUBLE PRECISION ZFACS,ZFACB,ZFACBB,ZFACL
      COMPLEX RKZSF
      COMPLEX ZFAC
      DOUBLE PRECISION VPA,VMA
      DOUBLE PRECISION RR1(0:4),RR2(0:4)
      DOUBLE PRECISION ZD12V,ZD21V,ZD1V2,ZD2V1,ZDV12,ZDV21
      COMPLEX RKZF,ZN12V,ZN21V,ZN1V2,ZN2V1,ZNV12,ZNV21
      COMPLEX ZDIA1,ZDIA2,ZDIA3,ZDIA4,ZDIA5,ZDIA6,ZDIA7,ZDIA8
      COMPLEX ZC12V,ZC21V,ZCV12,ZCV21
      DOUBLE PRECISION S,ZD11,ZD22
      COMPLEX ZABEL,ZNABEL,ZNABEM
      REAL RESULT
      DOUBLE PRECISION THIS1
      COMPLEX ANSS(-1:1,1:4,-1:1,1:4)
      INTEGER DONS(-1:1,1:4,-1:1,1:4)
      COMPLEX ANSF(-1:1,1:4,1:8,-1:1,1:4)
      INTEGER DONF(-1:1,1:4,1:8,-1:1,1:4)
 
C      PARAMETER(CHKGL1=0,CHKGL2=0)
      PARAMETER(NMC=1)
 
      COMMON/RKZSCO/ANSS,DONS
      COMMON/RKZFCO/ANSF,DONF
      COMMON/RKBBVC/RMQ,RMV,RGV,VB,AB,VL,AL
      DATA INIT/0/
 
* CHECK ON EITHER FIRST CALL OR CHANGE IN IMC
      IF(INIT.EQ.0.OR.IMC.NE.OLDIMC) THEN
        OLDIMC=IMC
        INIT=1
* REPRODUCE INPUT DATA
C       WRITE(6,*) ' ----------------------------------------'
C       WRITE(6,*) ' BBV: G G ---> B B_BAR Z, Z ---> L L_BAR'
C       WRITE(6,*) ' B QUARK MASS      = ',RMB,' GEV'
C       WRITE(6,*) ' BOSON MASS        = ',RMV,' GEV'
C       WRITE(6,*) ' BOSON WIDTH       = ',RGV,' GEV'
C       WRITE(6,*) ' B VECTOR C.       = ',VB
C       WRITE(6,*) ' B AXIAL C.        = ',AB
C       WRITE(6,*) ' LEPTON VECTOR C.  = ',VL
C       WRITE(6,*) ' LEPTON AXIAL C.   = ',AL
        RMB=RMQ
* ADJUST STRONG COUPLING SO AS TO GIVE EFFECTIVELY ALPHA_S=1
        GSTR=4D0*DSQRT(DATAN(1D0))
C       WRITE(6,*) ' QCD COUPLING      = ',GSTR
* SEE WETHER GAUGE CHECKS ARE REQUIRED
C        IF(CHKGL1.EQ.1) THEN
C          WRITE(6,*) ' GAUGE CHECK ON GLUON 1'
C        ENDIF
C        IF(CHKGL2.EQ.1) THEN
C          WRITE(6,*) ' GAUGE CHECK ON GLUON 2'
C        ENDIF
* SEE WETHER HELICITY MONTE CARLO IS REQUIRED
        IF(IMC.EQ.0) THEN
          KLOW=1
          KUPP=32
          MULT=1D0
          WRITE(6,*) ' SUM OVER HELICITIES SELECTED'
        ELSEIF(IMC.EQ.1) THEN
          KLOW=1
          KUPP=NMC
          MULT=32D0/(1D0*NMC)
C         WRITE(6,*) ' MONTE CARLO OVER HELICITES SELECTED'
C         WRITE(6,*) ' WITH ',NMC,' HELICITY TRIALS'
C         WRITE(6,*) ' RESULT THEN MULTIPLIED BY ',MULT
        ELSE
          WRITE(6,*) ' ERROR: WRONG OPTION IMC=',IMC
        ENDIF
C       WRITE(6,*) ' THE RESULT IS BASED ON ALPHA_S=1,',
C    .  ' MUST BE MULTIPLIED BY ALPHA_S**2'
C       WRITE(6,*) ' ----------------------------------------'
C       WRITE(6,800)'NO.','LG1','LG2','LV','L1','L2','AMP**2'
C 800   FORMAT(' ',6A4,A10)
      ENDIF
 
* INITIALIZE THE ARRAYS ANSS,DONS
      DO 130 J1=-1,1,2
        DO 120 J2=1,4
          DO 110 J3=-1,1,2
            DO 100 J4=1,4
              ANSS(J1,J2,J3,J4)=(0.,0.)
              DONS(J1,J2,J3,J4)=0
  100       CONTINUE
  110     CONTINUE
  120   CONTINUE
  130 CONTINUE
 
* INITIALIZE THE ARRAYS ANSF,DONF
      DO 180 J1=-1,1,2
        DO 170 J2=1,4
          DO 160 J3=1,8
            DO 150 J4=-1,1,2
              DO 140 J5=1,4
                 ANSF(J1,J2,J3,J4,J5)=(0.,0.)
                 DONF(J1,J2,J3,J4,J5)=0
  140         CONTINUE
  150       CONTINUE
  160     CONTINUE
  170   CONTINUE
  180 CONTINUE
 
* EQUATE THE (0:4) INTERNAL MOMENTA TO THE (0:3) ARGUMENTS MOMENTA
      DO 190 K=0,3
        K1(K)=AK1(K)
        K2(K)=AK2(K)
        P1(K)=AP1(K)
        P2(K)=AP2(K)
        LEP1(K)=ALEP1(K)
        LEP2(K)=ALEP2(K)
  190 CONTINUE
 
* ASSIGN LABELS TO THE MOMENTA FOR RECOGNITION
* THE MOMENTA K1,K2,LEP1,LEP2 (AND R1,R2) CAN OCCUR AS THE MASSLESS
* MOMENTA IN ARGUMENTS NO.2 AND 6 IN ZF, AND NO.2 AND 4 IN RKZSF
* R1,R2 AND Q1,Q2 ARE SOME OF THESE, AND CAN ALSO OCCUR
* AS ARGUMENTS NO.2 AND 6 IN ZF AND NO.2 AND 4 IN RKZSF
        K1(4)=1D0
        K2(4)=2D0
        LEP1(4)=3D0
        LEP2(4)=4D0
* THE OTHER MOMENTA P1,P2 AND THE VARIOUS RR1,RR2 CAN OCCUR ONLY
* AS ARGUMENT NO.3 IN ZF
        P1(4)=1D0
        P2(4)=2D0
 
* THE TOTAL BOSON MOMENTUM
* NO NEED TO ASSIGN 4TH COMPONENT LABEL SINCE IT IS NOT USED
      DO 200 K=0,3
        QV(K)=LEP1(K)+LEP2(K)
  200 CONTINUE
 
 
* DEFINE THE AUXILIARY VECTORS: THE RESULT SHOULD BE THE SAME
* FOR EVERY NON-SINGULAR CHOICE OF THE AUXILIARY VECTORS
* SINGULAR CHOICES ARE R1=K1 OR R2=K2
* THESE ARE OBTAINED BY PUTTING CHKGL1=1 OR CHKGL2=1
 
* AUXILIARY VECTOR FOR GLUON 1
* NEED TO ASSIGN ALSO 4TH COMPONENT LABELS HERE!
C      IF(CHKGL1.EQ.1) THEN
C        DO 210 K=0,4
C          R1(K)=K1(K)
C  210   CONTINUE
C      ELSE
        DO 210 K=0,4
          R1(K)=K2(K)
  210   CONTINUE
C      ENDIF
 
* AUXILIARY VECTOR FOR GLUON 2
C      IF(CHKGL2.EQ.1) THEN
C        DO 230 K=0,4
C          R2(K)=K2(K)
C  230   CONTINUE
C      ELSE
        DO 220 K=0,4
          R2(K)=K1(K)
  220   CONTINUE
C      ENDIF
 
* AUXILIARY VECTOR FOR THE B QUARK
      DO 230 K=0,4
        Q1(K)=LEP1(K)
  230 CONTINUE
 
* AUXILIARY VECTOR FOR THE B_BAR QUARK
      DO 240 K=0,4
        Q2(K)=LEP2(K)
  240 CONTINUE
 
* INITIALIZE THE CROSS SECTION TO ZERO
      CROSS=0D0
 
* SINCE P2 CORRESPONDS TO AN ANTIFERMION WE HAVE TO
* CHANGE ITS SIGN MOMENTARILY: PUT THE OLD RESULT IN PP2(0:3)
* BU MAKE SURE TO KEEP THE LABEL POSITIVE!
      DO 250 K=0,3
        PP2(K)=P2(K)
        P2(K)=-P2(K)
  250 CONTINUE
 
* COMPUTE OVERALL FACTORS: FOR EVERY SLASHED POLARIZATION THERE
* APPEARS A FACTOR OF 2 IN ADDITION TO THE NORMALIZATION
* FOLLOWING FROM THE CHISHOLM IDENTITY
* IN PRINCIPLE THE OVERALL FACTORS ARE DIFFERENT FOR EACH DIFFERENT
* HELICITY COMPBINATION BUT IN THIS CASE WE ARE ONLY INTERESTED IN
* THEIR ABSOLUTE VALUE (NO TRANSVERSE GLUON POLARIZATION ETC.)
* SO WE CAN TAKE THIS OUT OF THE LOOP, EXCEPT FOR THE NONTRIVIAL
* HELICITY DEPENDENCE IN 'ZFACV'
 
* OVERALL FACTOR FOR THE BOSON CURRENT, WITH BREIT-WIGNER
      ZFACV=2./CMPLX(SNGL(RKDOT(QV,QV))-RMV**2,RMV*RGV)
 
* OVERALL FACTOR FOR GLUON 1
C      IF(CHKGL1.EQ.1) THEN
C        ZFAC1=(1.,0.)
C      ELSE
* ORIGINAL FORM: ZFAC1=2D0*LG1/(DSQRT(2D0)*RKZPR(-LG1,K1,R1))
        ZFAC1=DSQRT(2D0)/RKZSF(1,K1,-1,R1)
C      ENDIF
 
* OVERALL FACTOR FOR GLUON 2
C      IF(CHKGL2.EQ.1) THEN
C        ZFAC2=1D0
C      ELSE
* ORIGINAL FORM: ZFAC2=2D0*LG2/(DSQRT(2D0)*RKZPR(-LG2,K2,R2))
        ZFAC2=DSQRT(2D0)/RKZSF(1,K2,-1,R2)
C      ENDIF
 
* OVERALL FACTOR FOR QCD COUPLINGS
      ZFACS=GSTR**2
 
* OVERALL FACTOR FOR THE B QUARK
      ZFACB=1/DSQRT(2D0*RKDOT(P1,Q1))
 
* OVERALL FACTOR FOR THE B_BAR QUARK
      ZFACBB=1D0/DSQRT(2D0*RKDOT(PP2,Q2))
 
* FINAL OVERALL FACTOR
      ZFAC=ZFACV*ZFAC1*ZFAC2*ZFACS*ZFACB*ZFACBB
 
* DO A BIG LOOP OVER ALL HELICITIES OR A RANDOM CHOICE OF HELICITIES
* NB: FUNNY INDENTATION HERE!
* ALSO INITIALIZE COUNTERS FOR RKZSF AND ZF
 
      DO 340 HELIX=KLOW,KUPP
      IF(IMC.EQ.0) THEN
        CALL RKHLPK(HELIX,LG1,LG2,LV,L1,L2)
      ELSE
        HELI=IDINT(32D0*RKRAND(HELIX))+1
        CALL RKHLPK(HELI,LG1,LG2,LV,L1,L2)
      ENDIF
 
* DETERMINE THE 'LEFT-' AND 'RIGHT-'HANDED COUPLINGS OF THE B TO THE Z
      VPA=VB+LV*AB
      VMA=VB-LV*AB
* AND THE LEPTON HELICITY FACTOR
      ZFACL=(VL-LV*AL)
 
* FIRST PART OF THE RESULT: THE ABELIAN TERMS
* COMPUTE THE NUMERATORS (ZN...) USING THE ZF FUNCTION
* AND THE DENOMINATORS (ZD...) THE STANDARD WAY
* THE INTERNAL FERMION MOMENTA ARE DIFFERENT IN EACH DIAGRAM
* AND ARE DENOTED BY RR1 AND RR2
* THE 4TH COMPONENT LABELS ARE NONTRIVIAL HERE: HAVING ALREADY
* P1(4)=1 AND P2(4)=2 WE ALSO DEFINE
* (P1-K1)(4)=3,
* (P1-K1-K2)(4)=(P1-K2-K1)(4)=4
* (P1-K2)(4)=5
* (P1-K1+QV)(4)=6
* (P1-K2+QV)(5)=7
* (P1+QV)(4)=8
* SO THAT IN THE VARIOUS DIAGRAMS WE HAVE
* IN ZN12V: RR1(4)=3, RR2(4)=4
* IN ZN21V: RR1(4)=5, RR2(4)=4
* IN ZN1V2: RR1(4)=3, RR2(4)=6
* IN ZN2V1: RR1(4)=5, RR2(4)=7
* IN ZNV12: RR1(4)=8, RR2(4)=6
* IN ZNV21: RR1(4)=8, RR2(4)=7
 
      DO 260 K=0,3
        RR1(K)=P1(K)-K1(K)
        RR2(K)=RR1(K)-K2(K)
  260 CONTINUE
      RR1(4)=3D0
      RR2(4)=4D0
      ZD12V=(RKDOT(RR1,RR1)-RMB**2)*(RKDOT(RR2,RR2)-RMB**2)
      ZN12V =
     . + RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
 
      DO 270 K=0,3
        RR1(K)=P1(K)-K2(K)
        RR2(K)=RR1(K)-K1(K)
  270 CONTINUE
      RR1(4)=5D0
      RR2(4)=4D0
      ZD21V=(RKDOT(RR1,RR1)-RMB**2)*(RKDOT(RR2,RR2)-RMB**2)
      ZN21V =
     .   RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
 
      DO 280 K=0,3
        RR1(K)=P1(K)-K1(K)
        RR2(K)=RR1(K)+QV(K)
  280 CONTINUE
      RR1(4)=3D0
      RR2(4)=6D0
      ZD1V2=(RKDOT(RR1,RR1)-RMB**2)*(RKDOT(RR2,RR2)-RMB**2)
      ZN1V2 =
     .   RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,LG2,R2)  *RKZF(LG2,K2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,-LG2,K2) *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,LG2,R2) *RKZF(LG2,K2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1)     *RKZF(LG1,K1,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,-LG2,K2)*RKZF(-LG2,R2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,LG2,R2)  *RKZF(LG2,K2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,-LG2,K2) *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,LG2,R2) *RKZF(LG2,K2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)    *RKZF(-LG1,R1,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,-LG2,K2)*RKZF(-LG2,R2,P2,RMB,L2,Q2)*VMA
 
      DO 290 K=0,3
        RR1(K)=P1(K)-K2(K)
        RR2(K)=RR1(K)+QV(K)
  290 CONTINUE
      RR1(4)=5D0
      RR2(4)=7D0
      ZD2V1=(RKDOT(RR1,RR1)-RMB**2)*(RKDOT(RR2,RR2)-RMB**2)
      ZN2V1 =
     .   RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,LG1,R1)  *RKZF(LG1,K1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,-LG1,K1) *RKZF(-LG1,R1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,LG1,R1) *RKZF(LG1,K1,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2)     *RKZF(LG2,K2,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,-LG1,K1)*RKZF(-LG1,R1,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,LG1,R1)  *RKZF(LG1,K1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,LV,LEP2)
     .  *RKZF(LV,LEP1,RR2,RMB,-LG1,K1) *RKZF(-LG1,R1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,LG1,R1) *RKZF(LG1,K1,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)    *RKZF(-LG2,R2,RR1,RMB,-LV,LEP1)
     .  *RKZF(-LV,LEP2,RR2,RMB,-LG1,K1)*RKZF(-LG1,R1,P2,RMB,L2,Q2)*VMA
 
      DO 300 K=0,3
        RR1(K)=P1(K)+QV(K)
        RR2(K)=RR1(K)-K1(K)
  300 CONTINUE
      RR1(4)=8D0
      RR2(4)=6D0
      ZDV12=(RKDOT(RR1,RR1)-RMB**2)*(RKDOT(RR2,RR2)-RMB**2)
      ZNV12 =
     .   RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,LG2,R2)  *RKZF(LG2,K2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,-LG2,K2) *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,LG2,R2) *RKZF(LG2,K2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,-LG2,K2)*RKZF(-LG2,R2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,LG2,R2)  *RKZF(LG2,K2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,LG1,R1)
     .  *RKZF(LG1,K1,RR2,RMB,-LG2,K2) *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,LG2,R2) *RKZF(LG2,K2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,-LG1,K1)
     .  *RKZF(-LG1,R1,RR2,RMB,-LG2,K2)*RKZF(-LG2,R2,P2,RMB,L2,Q2)*VMA
 
      DO 310 K=0,3
        RR1(K)=P1(K)+QV(K)
        RR2(K)=RR1(K)-K2(K)
  310 CONTINUE
      RR1(4)=8D0
      RR2(4)=7D0
      ZDV21=(RKDOT(RR1,RR1)-RMB**2)*(RKDOT(RR2,RR2)-RMB**2)
      ZNV21 =
     .   RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,LG1,R1)  *RKZF(LG1,K1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,-LG1,K1) *RKZF(-LG1,R1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,LG1,R1) *RKZF(LG1,K1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2)   *RKZF(LV,LEP1,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,-LG1,K1)*RKZF(-LG1,R1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,LG1,R1)  *RKZF(LG1,K1,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,LG2,R2)
     .  *RKZF(LG2,K2,RR2,RMB,-LG1,K1) *RKZF(-LG1,R1,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,LG1,R1) *RKZF(LG1,K1,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)  *RKZF(-LV,LEP2,RR1,RMB,-LG2,K2)
     .  *RKZF(-LG2,R2,RR2,RMB,-LG1,K1)*RKZF(-LG1,R1,P2,RMB,L2,Q2)*VMA
 
* COMPUTE THE DIAGRAMS SO FAR
      ZDIA1=ZN12V/ZD12V
      ZDIA2=ZN21V/ZD21V
      ZDIA3=ZN1V2/ZD1V2
      ZDIA4=ZN2V1/ZD2V1
      ZDIA5=ZNV12/ZDV12
      ZDIA6=ZNV21/ZDV21
 
* SECOND PART OF THE RESULT: THE NONABELIAN PART.
* THIS IS MADE UP PARTLY FROM THE ABELIAN PART AND PARTLY FROM
* NEW PIECES
* THE ASSIGNMENT OF THE 4TH COMPONENT LABELS IS NOW UNNECESSARY
* FOR RR1 SINCE IT DOES NOT OCCUR IN ANY ZF HERE
 
      S=2D0*RKDOT(K1,K2)
 
      DO 320 K=0,3
        RR1(K)=PP2(K)+QV(K)
  320 CONTINUE
      ZD11=S*(RKDOT(RR1,RR1)-RMB**2)
 
      ZC12V =
     . + RKZF(L1,Q1,P1,RMB,LG1,R1) *RKZSF(LG1,K1,LG2,R2)
     .  *RKZSF(LG2,K2,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1) *RKZSF(LG1,K1,LG2,R2)
     .  *RKZSF(LG2,K2,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1) *RKZSF(LG1,K1,-LG2,K2)
     .  *RKZSF(-LG2,R2,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG1,R1) *RKZSF(LG1,K1,-LG2,K2)
     .  *RKZSF(-LG2,R2,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)*RKZSF(-LG1,R1,LG2,R2)
     .  *RKZSF(LG2,K2,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)*RKZSF(-LG1,R1,LG2,R2)
     .  *RKZSF(LG2,K2,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)*RKZSF(-LG1,R1,-LG2,K2)
     .  *RKZSF(-LG2,R2,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG1,K1)*RKZSF(-LG1,R1,-LG2,K2)
     .  *RKZSF(-LG2,R2,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
 
      ZC21V =
     . + RKZF(L1,Q1,P1,RMB,LG2,R2) *RKZSF(LG2,K2,LG1,R1)
     .  *RKZSF(LG1,K1,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2) *RKZSF(LG2,K2,LG1,R1)
     .  *RKZSF(LG1,K1,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2) *RKZSF(LG2,K2,-LG1,K1)
     .  *RKZSF(-LG1,R1,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LG2,R2) *RKZSF(LG2,K2,-LG1,K1)
     .  *RKZSF(-LG1,R1,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)*RKZSF(-LG2,R2,LG1,R1)
     .  *RKZSF(LG1,K1,LV,LEP2)  *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)*RKZSF(-LG2,R2,LG1,R1)
     .  *RKZSF(LG1,K1,-LV,LEP1) *RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)*RKZSF(-LG2,R2,-LG1,K1)
     .  *RKZSF(-LG1,R1,LV,LEP2) *RKZF(LV,LEP1,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LG2,K2)*RKZSF(-LG2,R2,-LG1,K1)
     .  *RKZSF(-LG1,R1,-LV,LEP1)*RKZF(-LV,LEP2,P2,RMB,L2,Q2)*VMA
      ZDIA7=(-ZN12V+ZN21V)/ZD11-(ZC12V-ZC21V)/(2D0*S)
 
      DO 330 K=0,3
        RR1(K)=P1(K)+QV(K)
  330 CONTINUE
      ZD22=S*(RKDOT(RR1,RR1)-RMB**2)
 
      ZCV12 =
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2) *RKZSF(LV,LEP1,LG1,R1)
     .  *RKZSF(LG1,K1,LG2,R2)    *RKZF(LG2,K2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2) *RKZSF(LV,LEP1,LG1,R1)
     .  *RKZSF(LG1,K1,-LG2,K2)   *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2) *RKZSF(LV,LEP1,-LG1,K1)
     .  *RKZSF(-LG1,R1,LG2,R2)   *RKZF(LG2,K2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,LV,LEP2) *RKZSF(LV,LEP1,-LG1,K1)
     .  *RKZSF(-LG1,R1,-LG2,K2)  *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VPA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)*RKZSF(-LV,LEP2,LG1,R1)
     .  *RKZSF(LG1,K1,LG2,R2)    *RKZF(LG2,K2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)*RKZSF(-LV,LEP2,LG1,R1)
     .  *RKZSF(LG1,K1,-LG2,K2)   *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)*RKZSF(-LV,LEP2,-LG1,K1)
     .  *RKZSF(-LG1,R1,LG2,R2)   *RKZF(LG2,K2,P2,RMB,L2,Q2)*VMA
     . + RKZF(L1,Q1,P1,RMB,-LV,LEP1)*RKZSF(-LV,LEP2,-LG1,K1)
     .  *RKZSF(-LG1,R1,-LG2,K2)  *RKZF(-LG2,R2,P2,RMB,L2,Q2)*VMA
 
* THE FOURTH COMBINATION CAN BE GOTTEN FROM
* THE FIRST THREE USING DIRAC ALGEBRA:
* EPS1*EPS2*EPVS+EPS2*EPS1*EPSV = 2(EPS1.EPS2)*EPSV ETC.
      ZCV21=ZC12V+ZC21V-ZCV12
 
      ZDIA8=(-ZNV12+ZNV21)/ZD22-(ZCV12-ZCV21)/(2D0*S)
 
* CONSTRUCT THE ABELIAN AND NONABELIAN PART
 
      ZABEL= ZDIA1+ZDIA2+ZDIA3+ZDIA4+ZDIA5+ZDIA6
      ZNABEL=ZDIA1-ZDIA2+ZDIA3-ZDIA4+ZDIA5-ZDIA6
      ZNABEM=2D0*ZDIA7+2D0*ZDIA8
      ZNABEL=ZNABEL-ZNABEM
      ZABEL=ZABEL*ZFAC*ZFACL
      ZNABEL=ZNABEL*ZFAC*ZFACL
 
* INCLUDE COLOUR FACTORS:
* (N**2-1)*(N**2-2)/(8*N) = 7/3 FOR THE ABELIAN PART
* N*(N**2-1)/8 = 3 FOR THE NONABELIAN PART
* AND ADD THE RESULT TO THE CROSS SECTION
      THIS1=7D0/3D0*ABS(ZABEL)**2+3D0*ABS(ZNABEL)**2
CC    WRITE(6,801)HELIX,LG1,LG2,LV,L1,L2,THIS1
CC801 FORMAT(' ',6I4,D30.20)
      CROSS=CROSS+THIS1
 
* END OF THE BIG LOOP OVER HELICITIES
  340 CONTINUE
 
* DO NOT FORGET TO PUT P2 BACK TO ITS ORIGINAL VALUE IN PP2!
      DO 350 K=0,3
        P2(K)=PP2(K)
  350 CONTINUE
 
* ADD AVERAGING FACTORS:
* 1/2 FOR EACH GLUON SPIN, 1/8 FOR EACH GLUON COLOUR
      CROSS=CROSS/256D0
 
* TAKE INTO ACCOUNT A POSSIBLE FACTOR FOR THE HELICITY SUM OPTION
* AND RETURN THE FINAL RESULT
      IF(IMC.EQ.1) CROSS=CROSS*MULT
      RESULT=CROSS
      END