C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************
C***********************************************************************
C***********************************************************************
C
C                            JETSUBS
C                  Version 3.0,  23 February 1991
C                  Version 3.01, 14 March 1991:
C                     - misc small changes to enhance portability
C                  Version 3.1,  26 February 1992:
C                     - changes in soft integrals to match paper
C                     - a few other changes to match paper and simplify
C                     - common block /ARRAYS/ enlarged
C                     - switch for proton-proton collisions
C                  Version 4.0, 24 August 1994
C                                 
C***********************************************************************
C
C   Subroutines for JET.  Those that one is likely to want to modify,
C   which are consequently kept with the main program, are designated
C   with a (*).
C
C     SUBROUTINE RENO(NRENO)
C
C      SUBROUTINE INTEGRATE2TO2(Y1,Y2,PHI2,P2,JACOBIAN)
C        Virtual:
C         FUNCTION PSITILDE4(PROCESS,K)
C         FUNCTION PSITILDE6NS(PROCESS,K)
C        Collinear:
C         FUNCTION EFFPARTON(A,X,SCALE)
C            FUNCTION ALTARELLI(AOUT,AIN,Z)
C            FUNCTION ALTRLIPRIME(AOUT,AIN,Z)
C         FUNCTION GAMMA(NPARTON)
C         FUNCTION GAMMAPRIME(NPARTON)
C        Soft:
C         FUNCTION PSICTILDE(PROCESS,I1,I2,K)
C         FUNCTION LI2(X)
C
C      SUBROUTINE INTEGRATEA(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF)
C         FUNCTION FA(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB)
C      SUBROUTINE INTEGRATEB(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF)
C         FUNCTION FB(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB)
C      SUBROUTINE INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF)
C         FUNCTION F1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUA,MUB)
C      SUBROUTINE INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF)
C         FUNCTION F2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUA,MUB)
C
C Used in INTEGRATE2TO2,INTEGRATEA,INTEGRATEB,INTEGRATE1,INTEGRATE2.
C           SUBROUTINE JETDEFBORN(Y1,P2,Y2,PHI2,
C     >                       OK,JETVAR1,JETVAR2,JETVAR3,SVALUE) (*)
C           SUBROUTINE JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,
C     >                       OK,JETVAR1,JETVAR2,JETVAR3,SVALUE)  (*)
C           SUBROUTINE SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUA,MUB) (*)
C           SUBROUTINE BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND) (*)
C   The JETDEF subroutines calculate the 'measured' variables, like
C   PJ,YJ, given the parton momenta.  Then SETMU uses the measured
C   variables to set the mu values.  At this point the physics functions
C   can be calculated.  Finally, the results are smeared a bit over 
C   the physics variables and stored by BINIT.
C
C Used in FA,FB,F1,F2:
C      FUNCTION CROSSINGSIGN(PROCESS,PERMINV)
C      FUNCTION CONVERT(PHI)
C      FUNCTION RESIDUE(PROCESS,J,N,M,KIN)
C        SUBROUTINE RESIDUEINIT
C           LOGICAL FUNCTION MATCH(J,N,M,JJ,NN,MM)
C        SUBROUTINE PERMUTE
C      SUBROUTINE PARTONSIN(LMNSTY,XA,XB,MUA,MUB)
C      FUNCTION LUMINOSITY(LMNSTY,PROCESS,PERMINV)
C
C Auxiliary functions:
C      FUNCTION PARTON(NPARTON,X,SCALE)  Parton distributions (*)
C      FUNCTION ALPHAS(Q)  Alpha-s
C      FUNCTION WNUM(NPARTON)   Spin and color weight
C      FUNCTION COLOR(NPARTON)  Color charge
C      FUNCTION THETA(BOOL)  Theta function
C      FUNCTION FEXP(Y)  Exponential that doesn't crash for large y
C      SUBROUTINE PRINTVERSION
C
C Smearing:
C      FUNCTION UNSMEAR(M,J,DELTA)
C        FUNCTION DELTA(M,J)
C        FUNCTION FACTORIAL(N)
C Fitting:
C     SUBROUTINE FITJET(VALARRAY,ERRARRAY,
C    >                  V1GRID,V2GRID,V3GRID,NV1,NV2,NV3,
C    >                  V1SMEAR,V2SMEAR,V3SMEAR,
C    >                  J1MAX,J2MAX,J3MAX,CTILDE,C,NRENO) (*)
C        SUBROUTINE FITVAL(FUN,VAL,DEV,NPT,NPAR,PARAMS,CHISQ)
C          SUBROUTINE MMIN(FUNGRD,FUNHES,DELVEC,N)
C Random numbers:
C      FUNCTION RANDOM(DUMMY)
C        SUBROUTINE RANDOMINIT(IRAN)
C        SUBROUTINE NEWRAN
C
C***********************************************************************
C***********************************************************************
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE PRINTVERSION
C 
C
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
      WRITE(NOUT,1)
1     FORMAT(' JETSUBS version 4.0, 24 August 1994 '/,50('-'))
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE RENO(NRENO)
C 
      INTEGER NRENO
C
C   (Azimuthal angles are defined to lie between - pi and pi)
C
      DOUBLE PRECISION RANDOM
      DOUBLE PRECISION X1,X2,X3,X4,X5,X6,X7
      INTEGER ICALL1,ICALL2
C
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      DOUBLE PRECISION YSCALE,PSCALE
      COMMON /RENOSCALES/ YSCALE,PSCALE
C
      LOGICAL S2TO2
      LOGICAL SAFINITE,SBFINITE,S1FINITE,S2FINITE
      COMMON /SW2TO3/S2TO2,SAFINITE,SBFINITE,S1FINITE,S2FINITE
C
      DOUBLE PRECISION TINY
      DOUBLE PRECISION Y1,Y2,Y3,PHI1,PHI2,PHI3,P1,P2,P3,W,XI
      DOUBLE PRECISION COS0,PHI0,SIN0,X0,Y0,Z0
      DOUBLE PRECISION DY1DX1,DP2DX2,DY2DX3,DPHI2DX4,DP3DX5
      DOUBLE PRECISION DY2DX1,DP1DX2,DY1DX3,DPHI1DX4
      DOUBLE PRECISION DXIDX5,DWDX6,DPHI3DX7
      DOUBLE PRECISION DOMEGADX67,DY3PHI3DOMEGA
      DOUBLE PRECISION CONVERT
C
      DOUBLE PRECISION JACOBIAN,SINGULARITY,COEF
      DOUBLE PRECISION FEXP
C
      INTEGER PROCESS,NPERM,I,J
C
      DOUBLE PRECISION PI
      PARAMETER ( PI = 3.141592654 D0)
C
      INTEGER I1,I2,I3
      INTEGER SIZE1,SIZE2,SIZE3
      DOUBLE PRECISION TEMPARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION XSECTARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION ERRARRAY(1:16,1:16,1:16)
      COMMON /ARRAYS/ TEMPARRAY,XSECTARRAY,ERRARRAY,SIZE1,SIZE2,SIZE3
C
C  Arrays for RESIDUE, initialized by RESIDUEINIT
C
      INTEGER JA(5,5,5),JB(5,5,5),JC(5,5,5),JD(5,5,5)
      INTEGER SUBSCRIPTA(5,5,5),SUBSCRIPTB(5,5,5)
      INTEGER SUBSCRIPTC(5,5,5),SUBSCRIPTD(5,5,5)
      INTEGER PERMA(5,5,5,5),PERMB(5,5,5,5)
      INTEGER PERMC(5,5,5,5),PERMD(5,5,5,5)
      COMMON /RESIDUEHELP/ JA,JB,JC,JD,
     > SUBSCRIPTA,SUBSCRIPTB,SUBSCRIPTC,SUBSCRIPTD,
     > PERMA,PERMB,PERMC,PERMD
C
C   Permutations for each process.  The separate definitions for
C   each process, with an equivalence statement, are made in order
C   to keep under 20 continuation lines, for reasons of program
C   protability.
C
      INTEGER PERMS(30,5,4)
      INTEGER PERMSA(30,5),PERMSB(30,5),PERMSC(30,5),PERMSD(30,5)
      EQUIVALENCE (PERMS(1,1,1),PERMSA(1,1))
      EQUIVALENCE (PERMS(1,1,2),PERMSB(1,1))
      EQUIVALENCE (PERMS(1,1,3),PERMSC(1,1))
      EQUIVALENCE (PERMS(1,1,4),PERMSD(1,1))
C
C   The permutaions will be translated into these arrays for later use.
C
      INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4)
      COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR
C
C -- SET UP SUM OVER PERMUTATIONS FOR EACH PROCESS IN ORDER:
C
C -- FIRST PROCESS A:  USE (1,3) <-> (2,4) AND
C -- ALSO USE C INVARIANCE (1,2) <-> (3,4)
C -- TO MAKE PI(1)<PI(2),PI(1)<PI(3),PI(1)<PI(4)
C
      DATA PERMSA/
     >1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,
     >2,2,2,2,2,2,3,3,4,4,5,5,3,3,4,4,5,5,3,3,4,4,5,5,3,3,4,4,5,5,
     >3,3,4,4,5,5,2,2,2,2,2,2,4,5,3,5,3,4,4,5,3,5,3,4,4,5,3,5,3,4,
     >4,5,3,5,3,4,4,5,3,5,3,4,2,2,2,2,2,2,5,4,5,3,4,3,5,4,5,3,4,3,
     >5,4,5,3,4,3,5,4,5,3,4,3,5,4,5,3,4,3,2,2,2,2,2,2,1,1,1,1,1,1/
C
C -- FOR PROCESS B LET PI(3) < PI(4) ALSO USING EXTRA SYMMETRY
C
      DATA PERMSB/
     >1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >2,2,2,3,3,4,4,5,5,3,4,5,3,4,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >3,3,4,2,2,2,2,2,2,4,3,3,4,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >4,5,5,4,5,3,5,3,4,5,5,4,5,5,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >5,4,3,5,4,5,3,4,3,2,2,2,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
C
C -- FOR PROCESS C USE 345 SYMMETRY TO MAKE PI(3) < PI(4) < PI(5)
C    AND C INVARIANCE TO MAKE PI(1)<PI(2)
C
      DATA PERMSC/
     >1,1,1,1,2,2,2,3,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >2,3,4,5,3,4,5,4,5,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >3,2,2,2,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >4,4,3,3,4,3,3,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >5,5,5,4,5,5,4,5,4,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
C
C -- FOR PROCESS D SYMMETRIC IN EVERYTHING
C
      DATA PERMSD/
     >1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     >5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
C
C   Initializations.
C
      NPERMS(1) = 30
      NPERMS(2) = 15
      NPERMS(3) = 10
      NPERMS(4) = 1
C
C   Initialize logic arrays for RESIDUE
C
      CALL RESIDUEINIT
C
C   Initialize results arrays
C
      DO I1 = 1,SIZE1
      DO I2 = 1,SIZE2
      DO I3 = 1,SIZE3
          XSECTARRAY(I1,I2,I3) = 0.0D0
          ERRARRAY(I1,I2,I3) = 0.0D0
          TEMPARRAY(I1,I2,I3) = 0.0D0
      ENDDO
      ENDDO
      ENDDO
C
C   Initialize permutation arrays
C
      DO PROCESS = 1,4
       DO NPERM = 1,NPERMS(PROCESS)
        DO I = 1,5
          J = PERMS(NPERM,I,PROCESS)
          PERMINVARR(J,NPERM,PROCESS) = I
          PERMARR(I,NPERM,PROCESS) = J
        ENDDO
       ENDDO
      ENDDO
C
C   Set parameters for choosing kinematic variables
C
      TINY = 1.0 D-10
C
C   Calculation starts here.
C
C   Start loop over ICALL1, ICALL2.
C
      DO ICALL2 = 1,NRENO
      DO ICALL1 = 1,1000
C
C  ----------------------- Do 2-to-2 terms ----------------------------
C
      IF (S2TO2) THEN
C
      X1 = RANDOM(1.0)
      X2 = RANDOM(1.0)
      X3 = RANDOM(1.0)
      X4 = RANDOM(1.0)
C
      Y1 = YSCALE * DLOG(  (X1 + TINY) /(1.0D0 + TINY - X1) )
      DY1DX1 =  YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1))
C
      P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)
      DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS)
     >         /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2
C
      Y2 = YSCALE * DLOG(  (X3 + TINY) /(1.0D0 + TINY - X3) )
      DY2DX3 =  YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3))
C
      PHI2 = 2.0D0 * PI * X4 - PI
      DPHI2DX4 = 2.0D0 * PI
C
C
      JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4
C
      CALL INTEGRATE2TO2(Y1,Y2,PHI2,P2,JACOBIAN)
C
      ENDIF
C
C  -------------------------- Do Term A ----------------------------
C
      IF (SAFINITE) THEN
C
      X1 = RANDOM(1.0)
      X2 = RANDOM(1.0)
      X3 = RANDOM(1.0)
      X4 = RANDOM(1.0)
      X5 = RANDOM(1.0)
      X6 = RANDOM(1.0)
      X7 = RANDOM(1.0)
C
      Y1 = YSCALE * DLOG(  (X1 + TINY) /(1.0D0 + TINY - X1) )
      DY1DX1 =  YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1))
C
      P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)
      DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS)
     >         /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2
C
      Y2 = YSCALE * DLOG(  (X3 + TINY) /(1.0D0 + TINY - X3) )
      DY2DX3 =  YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3))
C
      PHI2 = 2.0D0 * PI * X4 - PI
      DPHI2DX4 = 2.0D0 * PI
C
      XI = X5**2
      DXIDX5 = 2.0D0 * X5
C
      W = RTS * X6**2 /(1.0D0 - X6)
      DWDX6 = RTS*X6*(2.0D0 - X6) /(1.0D0 - X6)**2
C
      PHI3 = 2.0D0 * PI * X7 - PI
      DPHI3DX7 = 2.0D0 * PI
C
      JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4*DXIDX5*DWDX6*DPHI3DX7
C
      SINGULARITY = 1.0D0/( XI * W )
C
C   Proceed as long as SINGULARITY is not too large:
C
      IF (SINGULARITY.LT.1.0D6) THEN
C
      COEF = SINGULARITY*JACOBIAN
C
        CALL INTEGRATEA(Y1,P2,Y2,PHI2,XI,   W,    PHI3, COEF)
      IF (W.LT.RTS) THEN
        CALL INTEGRATEA(Y1,P2,Y2,PHI2,XI,   0.0D0,PHI3,-COEF)
      ENDIF
      IF (XI.LT.(1.0D0 - ( FEXP(Y1) + FEXP(Y2) ) * P2 / RTS)) THEN
        CALL INTEGRATEA(Y1,P2,Y2,PHI2,0.0D0,W,    PHI3,-COEF)
       IF (W.LT.RTS) THEN
        CALL INTEGRATEA(Y1,P2,Y2,PHI2,0.0D0,0.0D0,PHI3, COEF)
       ENDIF
      ENDIF
C
C  End IF for SINGULARITY not too big.
C
      ENDIF
C
C  End IF for SAFINITE switch on
C
      ENDIF
C
C  -------------------------- Do Term B ----------------------------
C
      IF (SBFINITE) THEN
C
      X1 = RANDOM(1.0)
      X2 = RANDOM(1.0)
      X3 = RANDOM(1.0)
      X4 = RANDOM(1.0)
      X5 = RANDOM(1.0)
      X6 = RANDOM(1.0)
      X7 = RANDOM(1.0)
C
      Y1 = YSCALE * DLOG(  (X1 + TINY) /(1.0D0 + TINY - X1) )
      DY1DX1 =  YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1))
C
      P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)
      DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS)
     >         /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2
C
      Y2 = YSCALE * DLOG(  (X3 + TINY) /(1.0D0 + TINY - X3) )
      DY2DX3 =  YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3))
C
      PHI2 = 2.0D0 * PI * X4 - PI
      DPHI2DX4 = 2.0D0 * PI
C
      XI = X5**2
      DXIDX5 = 2.0D0 * X5
C
      W = RTS * X6**2 /(1.0D0 - X6)
      DWDX6 = RTS*X6*(2.0D0 - X6) /(1.0D0 - X6)**2
C
      PHI3 = 2.0D0 * PI * X7 - PI
      DPHI3DX7 = 2.0D0 * PI
C
      JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4*DXIDX5*DWDX6*DPHI3DX7
C
      SINGULARITY = 1.0D0/( XI * W )
C
C   Proceed as long as SINGULARITY is not too large:
C
      IF (SINGULARITY.LT.1.0D6) THEN
C
      COEF = SINGULARITY*JACOBIAN
C
        CALL INTEGRATEB(Y1,P2,Y2,PHI2,XI,   W,    PHI3, COEF)
      IF (W.LT.RTS) THEN
        CALL INTEGRATEB(Y1,P2,Y2,PHI2,XI,   0.0D0,PHI3,-COEF)
      ENDIF
      IF (XI.LT.(1.0D0 - ( FEXP(-Y1) + FEXP(-Y2) ) * P2 / RTS)) THEN
        CALL INTEGRATEB(Y1,P2,Y2,PHI2,0.0D0,W,    PHI3,-COEF)
       IF (W.LT.RTS) THEN
        CALL INTEGRATEB(Y1,P2,Y2,PHI2,0.0D0,0.0D0,PHI3, COEF)
       ENDIF
      ENDIF
C
C  End IF for SINGULARITY not too big
C
      ENDIF
C
C  End IF for SBFINITE switch on
C
      ENDIF
C
C  -------------------------- Do Term 1 ----------------------------
C
      IF (S1FINITE) THEN
C
      X1 = RANDOM(1.0)
      X2 = RANDOM(1.0)
      X3 = RANDOM(1.0)
      X4 = RANDOM(1.0)
      X5 = RANDOM(1.0)
      X6 = RANDOM(1.0)
      X7 = RANDOM(1.0)
C
      Y2 = YSCALE * DLOG(  (X1 + TINY) /(1.0D0 + TINY - X1) )
      DY2DX1 =  YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1))
C
      P1 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)
      DP1DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS)
     >         /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2
C
      Y1 = YSCALE * DLOG(  (X3 + TINY) /(1.0D0 + TINY - X3) )
      DY1DX3 =  YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3))
C
      PHI1 = 2.0D0 * PI * X4 - PI
      DPHI1DX4 = 2.0D0 * PI
C
      P3 = P1 * X5**2
      DP3DX5 = 2.0D0 * P1 * X5
C
      COS0 = 1.0D0 - 2.0D0 * X6**2
      PHI0 =  2.0D0 * PI * X7 - PI
      DOMEGADX67 = 8.0D0 * PI * X6
C
      SIN0 = DSQRT(1.0D0 - COS0**2)
      X0 = COS0
      Y0 = SIN0 * DCOS(PHI0)
      Z0 = SIN0 * DSIN(PHI0)
C
      Y3 = Y1 + 2.0D0 * DLOG( (1.0D0 + Z0)/(1.0D0 - Z0) )
      PHI3 = CONVERT( PHI1 + DATAN2(Y0,X0))
C
      DY3PHI3DOMEGA = 4.0D0/(1.0D0 - Z0**2)
C
      JACOBIAN = DY2DX1*DP1DX2*DY1DX3*DPHI1DX4*DP3DX5
     >           * DOMEGADX67 * DY3PHI3DOMEGA
C
      SINGULARITY = 1.0D0/(P3 * ( DCOSH(Y1-Y3) - DCOS(PHI1-PHI3) ))
C
C   Proceed as long as SINGULARITY is not too large:
C
      IF (SINGULARITY.LT.1.0D6) THEN
C
      COEF = SINGULARITY*JACOBIAN
C
        CALL INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF)
        CALL INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y1,PHI1,-COEF)
      IF (P3.LT.0.5D0*P1) THEN
        CALL INTEGRATE1(Y2,P1,Y1,PHI1,0.0D0,Y3,PHI3,-COEF)
        CALL INTEGRATE1(Y2,P1,Y1,PHI1,0.0D0,Y1,PHI1,COEF)
      ENDIF
C
C  End IF for SINGULARITY not too big
C
      ENDIF
C
C  End IF for S1FINITE switch on
C
      ENDIF
C
C  -------------------------- Do Term 2 ----------------------------
C
      IF (S2FINITE) THEN
C
      X1 = RANDOM(1.0)
      X2 = RANDOM(1.0)
      X3 = RANDOM(1.0)
      X4 = RANDOM(1.0)
      X5 = RANDOM(1.0)
      X6 = RANDOM(1.0)
      X7 = RANDOM(1.0)
C
      Y1 = YSCALE * DLOG(  (X1 + TINY) /(1.0D0 + TINY - X1) )
      DY1DX1 =  YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1))
C
      P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)
      DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS)
     >         /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2
C
      Y2 = YSCALE * DLOG(  (X3 + TINY) /(1.0D0 + TINY - X3) )
      DY2DX3 =  YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3))
C
      PHI2 = 2.0D0 * PI * X4 - PI
      DPHI2DX4 = 2.0D0 * PI
C
      P3 = P2 * X5**2
      DP3DX5 = 2.0D0 * P2 * X5
C
      COS0 = 1.0D0 - 2.0D0 * X6**2
      PHI0 =  2.0D0 * PI * X7 - PI
      DOMEGADX67 = 8.0D0 * PI * X6
C
      SIN0 = DSQRT(1.0D0 - COS0**2)
      X0 = COS0
      Y0 = SIN0 * DCOS(PHI0)
      Z0 = SIN0 * DSIN(PHI0)
C
      Y3 = Y2 + 2.0D0 * DLOG( (1.0D0 + Z0)/(1.0D0 - Z0) )
      PHI3 = CONVERT( PHI2 + DATAN2(Y0,X0))
C
      DY3PHI3DOMEGA = 4.0D0/(1.0D0 - Z0**2)
C
      JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4*DP3DX5
     >           * DOMEGADX67 * DY3PHI3DOMEGA
C
      SINGULARITY = 1.0D0/(P3 * ( DCOSH(Y2-Y3) - DCOS(PHI2-PHI3) ))
C
C   Proceed as long as SINGULARITY is not too large:
C
      IF (SINGULARITY.LT.1.0D6) THEN
C
      COEF = SINGULARITY*JACOBIAN
C
        CALL INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF)
        CALL INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y2,PHI2,-COEF)
      IF (P3.LT.0.5D0*P2) THEN
        CALL INTEGRATE2(Y1,P2,Y2,PHI2,0.0D0,Y3,PHI3,-COEF)
        CALL INTEGRATE2(Y1,P2,Y2,PHI2,0.0D0,Y2,PHI2,COEF)
      ENDIF
C
C  End IF for SINGULARITY not too big
C
      ENDIF
C
C  End IF for S2FINITE switch on
C
      ENDIF
C
C  End loop over ICALL1 = 1,1000
C
      ENDDO
C
C   Update arrays
C
      DO I1 = 1,SIZE1
      DO I2 = 1,SIZE2
      DO I3 = 1,SIZE3
       XSECTARRAY(I1,I2,I3) = XSECTARRAY(I1,I2,I3)
     >             + (TEMPARRAY(I1,I2,I3)/1000)/NRENO
       ERRARRAY(I1,I2,I3) = ERRARRAY(I1,I2,I3)
     >             + (TEMPARRAY(I1,I2,I3)/1000)**2/NRENO
       TEMPARRAY(I1,I2,I3)  = 0.0D0
      ENDDO
      ENDDO
      ENDDO
C
C  End loop over ICALL2 = 1,NRENO
C
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
C                   Routines for 2-to-2 integrals
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE INTEGRATE2TO2(Y1,Y2,PHI2,P2,JACOBIAN)
C 
      DOUBLE PRECISION Y1,Y2,PHI2,P2,JACOBIAN
C
C     This function calculates the integrand multiplying 
C     DY1 DP2 DY2 DPHI2
C     in the calculation of the 2 -> 2 cross section
C
      DOUBLE PRECISION INTEGRAND2TO2
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C  Jet parameters
      LOGICAL OK
      DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE
C  Switches
      LOGICAL SBORN,SVIRTUAL
      LOGICAL SACOLL,SASOFT
      LOGICAL SBCOLL,SBSOFT
      LOGICAL S1COLL,S1SOFT
      LOGICAL S2COLL,S2SOFT
      COMMON /SW2TO2/ SBORN,SVIRTUAL,
     >   SACOLL,SASOFT,SBCOLL,SBSOFT,
     >   S1COLL,S1SOFT,S2COLL,S2SOFT
      LOGICAL STYPE(4)
      COMMON /SWTYPE/ STYPE
C  Switch for p-pbar collisions (if true) or p-p collisions (if false)
      LOGICAL PPBAR
      COMMON /BEAMS/ PPBAR
C   Common block for mu values for 2-to-2 subroutines
      DOUBLE PRECISION MUUV,MUA,MUB
      COMMON /MUVALS/ MUUV,MUA,MUB
C   Loop controls,permutations,parton labels
      INTEGER PERMS(12,4,4),NPERMS(4),PERMINV(4),PERM(4)
      INTEGER PROCESS,NPERM,FLAVORS1,FLAVORS2,A1,A2
      INTEGER N1,N2,I,J
C   Parton distributions
      DOUBLE PRECISION PRTNA(-6:6),PRTNB(-6:6)
      DOUBLE PRECISION EFFPRTNA(-6:6),EFFPRTNB(-6:6)
C   Useful parameters
      DOUBLE PRECISION PI,V,N,LOG2
      PARAMETER ( PI = 3.141592654 D0)
      PARAMETER ( V = 8.0 D0)
      PARAMETER ( N = 3.0 D0)
      PARAMETER ( LOG2 = 0.693147181 D0)
C   Units for error reports:
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C   Variables for collinear integration
      DOUBLE PRECISION EFFPARTON
C   Parton flavors
      INTEGER AA,AB
C   Kinematical variables
      DOUBLE PRECISION SHAT,UHAT,THAT
      DOUBLE PRECISION XA,XB
      DOUBLE PRECISION PIN(4,4),K(4,4)
      DOUBLE PRECISION FEXP,EY1,EY2
C   Ellis and Sexton scale parameter
      DOUBLE PRECISION QES
      COMMON /KEITHSQ/ QES
C   QCD functions
      DOUBLE PRECISION D,D0,ALPHAS
      DOUBLE PRECISION L,LEFFA,LEFFB
      DOUBLE PRECISION PSI4,PSITILDE4,SIGN
      DOUBLE PRECISION PSI6NS,PSITILDE6NS
      DOUBLE PRECISION GAMMAPRIME
      DOUBLE PRECISION PARTON,WNUM,COLOR,GAMMA
      DOUBLE PRECISION JAB,JA1,JA2,JBA,JB1,JB2,J1A,J1B,J12,J2A,J2B,J21
      DOUBLE PRECISION PSICAB,PSICA1,PSICA2,PSICB1,PSICB2,PSIC12
      DOUBLE PRECISION PSICTILDE
      DOUBLE PRECISION LI2
C   Precalculated logs
      DOUBLE PRECISION LOGS,LOGT,LOGU,LOGMUASQ,LOGMUBSQ,LOG16P2SQ
      DOUBLE PRECISION LOGXA,LOGXB
C
C      Each column is one of the permutations we need.
C
      DATA PERMS/1,1,1,1,1,1,2,2,2,2,3,3,
     >           2,2,3,3,4,4,3,3,4,4,4,4,
     >           3,4,2,4,2,3,1,4,1,3,1,2,
     >           4,3,4,2,3,2,4,1,3,1,2,1,
C
     >           1,1,1,2,2,3,0,0,0,0,0,0,
     >           2,3,4,3,4,4,0,0,0,0,0,0,
     >           3,2,2,1,1,1,0,0,0,0,0,0,
     >           4,4,3,4,3,2,0,0,0,0,0,0,
C
     >           1,1,1,2,2,2,3,3,3,4,4,4,
     >           2,3,4,1,3,4,1,2,4,1,2,3,
     >           3,2,2,3,1,1,2,1,1,2,1,1,
     >           4,4,3,4,4,3,4,4,2,3,3,2,
C
     >           1,0,0,0,0,0,0,0,0,0,0,0,
     >           2,0,0,0,0,0,0,0,0,0,0,0,
     >           3,0,0,0,0,0,0,0,0,0,0,0,
     >           4,0,0,0,0,0,0,0,0,0,0,0/
C
      DATA NPERMS/12,6,12,1/
C
C -- Initialize computation ------------
C
C
      CALL  JETDEFBORN(Y1,P2,Y2,PHI2,OK,JETVAR1,JETVAR2,JETVAR3,SVALUE)
      IF(.NOT.OK) THEN
       RETURN
      ENDIF
      CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUA,MUB)
C
      EY1 = FEXP(Y1)
      EY2 = FEXP(Y2)
      XA = (EY1 + EY2)* P2/RTS
      XB = (1.0D0/EY1 + 1.0D0/EY2)* P2/RTS
      SHAT =  P2**2 * ( 2.0D0 + EY2/EY1 + EY1/EY2 )
      THAT = -  P2**2 * (1.0D0 + EY2/EY1)
      UHAT = -  P2**2 * (1.0D0 + EY1/EY2)
C
      IF( (XA.GT.1.0D0).OR.(XB.GT.1.0D0))THEN
        RETURN
      ENDIF
C
      QES = RTS/10.0D0
C
C  Momentum matrix (incoming)
C
      PIN(1,1) =  0.0D0
      PIN(1,2) =  SHAT /2.0D0
      PIN(1,3) =  THAT /2.0D0
      PIN(1,4) =  UHAT /2.0D0
      PIN(2,1) =  SHAT /2.0D0
      PIN(2,2) =  0.0D0
      PIN(2,3) =  UHAT /2.0D0
      PIN(2,4) =  THAT /2.0D0
      PIN(3,1) =  THAT /2.0D0
      PIN(3,2) =  UHAT /2.0D0
      PIN(3,3) =  0.0D0
      PIN(3,4) =  SHAT /2.0D0
      PIN(4,1) =  UHAT /2.0D0
      PIN(4,2) =  THAT /2.0D0
      PIN(4,3) =  SHAT /2.0D0
      PIN(4,4) =  0.0D0
C
      D0 = ALPHAS(MUUV)**2 * P2 /2.0D0 /RTS**4
      D = D0 * ALPHAS(MUUV)/2.0D0/PI
C
      LOGS = DLOG(SHAT /QES**2)
      LOGT = DLOG(-THAT /QES**2)
      LOGU = DLOG(-UHAT /QES**2)
      LOGMUASQ = DLOG(MUA**2 /QES**2)
      LOGMUBSQ = DLOG(MUB**2 /QES**2)
      LOGXA = DLOG( (XA/(1.0D0 - XA))**2 )
      LOGXB = DLOG( (XB/(1.0D0 - XB))**2 )
      LOG16P2SQ = DLOG(16.0D0*P2**2/QES**2)
C
C  The tilde J_MN soft integrals
C
      J21 =  0.25D0 * (2.0D0 * LOG2 - LOGS)**2 
     >     - 0.25D0 * (2.0D0 * LOGS - LOGT - LOGU)**2
     >     + (LOGS - LOGT)*(LOGS - LOGU)
C
      J2A =  0.25D0 * (2.0D0 * LOG2 - LOGU)**2
     >     - 0.25D0 * (LOGS - LOGT)**2
     >     - LOG2 * DLOG(2.0D0 + UHAT/SHAT)
C
      J2B =  0.25D0 * (2.0D0 * LOG2 - LOGT)**2
     >     - 0.25D0 * (LOGS - LOGU)**2
     >     - LOG2 * DLOG(2.0D0 + THAT/SHAT)
C
      J12 = J21
      J1A = J2B
      J1B = J2A
C
      JAB =  0.25D0 * (  LOGXA - LOGS )**2 + PI**2 / 12.0D0
C
      JA1 =  0.25D0 * (LOGXA - LOGT)**2
     >     - 0.25D0 * DLOG(0.5D0 + 0.25D0 * THAT / SHAT)**2
     >     + LOG2**2  + PI**2 / 12.0D0
     >     - 0.5D0 * LI2( (SHAT + THAT)/(2.0D0*SHAT + THAT) )
C
      JA2 =  0.25D0 * (LOGXA - LOGU)**2
     >     - 0.25D0 * DLOG(0.5D0 + 0.25D0 * UHAT / SHAT)**2
     >     + LOG2**2  + PI**2 / 12.0D0
     >     - 0.5D0 * LI2( (SHAT + UHAT)/(2.0D0*SHAT + UHAT) )
C
C
      JBA =  0.25D0 * (  LOGXB - LOGS )**2 + PI**2 / 12.0D0
C
      JB1 =  JA2 - 0.25D0 * (LOGXA - LOGU)**2
     >           + 0.25D0 * (LOGXB - LOGU)**2
C
      JB2 =  JA1 - 0.25D0 * (LOGXA - LOGT)**2
     >           + 0.25D0 * (LOGXB - LOGT)**2
C
C  ---- Calculate parton distributions. For p-pbar collisions, exchange
C       partons -> antipartons for hadron B.
C
      DO AA = -NFL,NFL
        PRTNA(AA) = PARTON(AA,XA,MUA)/(XA*WNUM(AA))
        EFFPRTNA(AA) =  EFFPARTON(AA,XA,MUA)/(XA*WNUM(AA))
      ENDDO
C
      IF (PPBAR) THEN
       DO AB = -NFL,NFL
        PRTNB(AB) = PARTON(-AB,XB,MUB)/(XB*WNUM(AB))
        EFFPRTNB(AB) =  EFFPARTON(-AB,XB,MUB)/(XB*WNUM(AB))
       ENDDO
      ELSE
       DO AB = -NFL,NFL
        PRTNB(AB) = PARTON(AB,XB,MUB)/(XB*WNUM(AB))
        EFFPRTNB(AB) =  EFFPARTON(AB,XB,MUB)/(XB*WNUM(AB))
       ENDDO
      ENDIF
C
C  ----- Initialize sum:
C
      INTEGRAND2TO2 = 0.0D0
C
C  ----- Sum over processes A,B,C,D (PROCESS = 1,2,3,4) -------
C
      DO PROCESS = 1,4
      IF (STYPE(PROCESS)) THEN
C
C  ----- Sum over independent permutations the partons
C
      DO NPERM = 1,NPERMS(PROCESS)
C
C    Calculate inverse permutation
C
      DO I = 1,4
        J = PERMS(NPERM,I,PROCESS)
        PERMINV(J) = I
        PERM(I) = J
       ENDDO
C
C      Set K matrix as permutation of P matrix:
C      K(I) = P(J) with J = PERM(I).
C
           DO I = 1,4
            DO J = 1,4
             K(I,J) = PIN(PERMS(NPERM,I,PROCESS),
     >                    PERMS(NPERM,J,PROCESS))
            ENDDO
           ENDDO
C
C  Here we calculate functions that don't depend on the choice
C  of quark flavors. Note that these are the 'tilde' functions,
C  so we need a crossing sign when we use them.
C
      PSI4 = PSITILDE4(PROCESS,K)
      PSI6NS = PSITILDE6NS(PROCESS,K)
      PSICAB = PSICTILDE(PROCESS,PERMINV(1),PERMINV(2),K)
      PSICA1 = PSICTILDE(PROCESS,PERMINV(1),PERMINV(3),K)
      PSICA2 = PSICTILDE(PROCESS,PERMINV(1),PERMINV(4),K)
      PSICB1 = PSICTILDE(PROCESS,PERMINV(2),PERMINV(3),K)
      PSICB2 = PSICTILDE(PROCESS,PERMINV(2),PERMINV(4),K)
      PSIC12 = PSICTILDE(PROCESS,PERMINV(3),PERMINV(4),K)
C
C  Find how many flavors for incoming partons to sum over.
C
      IF (PROCESS.EQ.1) THEN
        FLAVORS1 = NFL
        FLAVORS2 = NFL
      ELSE IF (PROCESS.EQ.2) THEN
        FLAVORS1 = NFL
        FLAVORS2 = 1
      ELSE IF (PROCESS.EQ.3) THEN
        FLAVORS1 = NFL
        FLAVORS2 = 1
      ELSE IF (PROCESS.EQ.4) THEN
        FLAVORS1 = 1
        FLAVORS2 = 1
      ENDIF
C
C  ----- Sum over particular parton flavors to calculate luminosities
C
      L = 0.0D0
      LEFFA = 0.0D0
      LEFFB = 0.0D0
C
      DO N1 = 1,FLAVORS1
      DO N2 = 1,FLAVORS2
      IF (.NOT.((PROCESS.EQ.1).AND.(N1.EQ.N2))) THEN
C
C  Find flavors for incoming partons: a(J) = a_0(I) with J = PERM(I)
C  and set outgoing flavors to 0 (gluon) or 101 (generic
C  quark or antiquark).
C
       IF (PROCESS.EQ.1) THEN
             IF       (PERM(1).EQ.1) THEN
                 AA =   N1
              ELSE IF (PERM(2).EQ.1) THEN
                 AA =   N2
              ELSE IF (PERM(3).EQ.1) THEN
                 AA = - N1
              ELSE IF (PERM(4).EQ.1) THEN
                 AA = - N2
             END IF
             IF       (PERM(1).EQ.2) THEN
                 AB =   N1
              ELSE IF (PERM(2).EQ.2) THEN
                 AB =   N2
              ELSE IF (PERM(3).EQ.2) THEN
                 AB = - N1
              ELSE IF (PERM(4).EQ.2) THEN
                 AB = - N2
             END IF
             A1 = 101
             A2 = 101
       ELSE IF (PROCESS.EQ.2) THEN
             IF       (PERM(1).EQ.1) THEN
                 AA =   N1
              ELSE IF (PERM(2).EQ.1) THEN
                 AA =   N1
              ELSE IF (PERM(3).EQ.1) THEN
                 AA = - N1
              ELSE IF (PERM(4).EQ.1) THEN
                 AA = - N1
             END IF
             IF       (PERM(1).EQ.2) THEN
                 AB =   N1
              ELSE IF (PERM(2).EQ.2) THEN
                 AB =   N1
              ELSE IF (PERM(3).EQ.2) THEN
                 AB = - N1
              ELSE IF (PERM(4).EQ.2) THEN
                 AB = - N1
             END IF
             A1 = 101
             A2 = 101
       ELSE IF (PROCESS.EQ.3) THEN
             IF       (PERM(1).EQ.1) THEN
                 AA =   N1
              ELSE IF (PERM(2).EQ.1) THEN
                 AA = - N1
              ELSE IF (PERM(3).EQ.1) THEN
                 AA =   0
              ELSE IF (PERM(4).EQ.1) THEN
                 AA =   0
             END IF
             IF       (PERM(1).EQ.2) THEN
                 AB =   N1
              ELSE IF (PERM(2).EQ.2) THEN
                 AB = - N1
              ELSE IF (PERM(3).EQ.2) THEN
                 AB =   0
              ELSE IF (PERM(4).EQ.2) THEN
                 AB =   0
             END IF
             IF       (PERM(1).EQ.3) THEN
                 A1 =  101
              ELSE IF (PERM(2).EQ.3) THEN
                 A1 =  101
              ELSE IF (PERM(3).EQ.3) THEN
                 A1 =   0
              ELSE IF (PERM(4).EQ.3) THEN
                 A1 =   0
             END IF
             IF       (PERM(1).EQ.4) THEN
                 A2 =  101
              ELSE IF (PERM(2).EQ.4) THEN
                 A2 =  101
              ELSE IF (PERM(3).EQ.4) THEN
                 A2 =   0
              ELSE IF (PERM(4).EQ.4) THEN
                 A2 =   0
             END IF
       ELSE IF (PROCESS.EQ.4) THEN
             AA = 0
             AB = 0
             A1 = 0
             A2 = 0
      ENDIF
C
C  Calculate luminosity factors
C
       L = L + PRTNA(AA)*PRTNB(AB)
       LEFFA = LEFFA + EFFPRTNA(AA)*PRTNB(AB)
       LEFFB = LEFFB + PRTNA(AA)*EFFPRTNB(AB)
C
C  Here we close the loops for flavor sums, having calculated the
C  luminosity factors.  In the following parts of the calculation,
C  the indices AA,AB,A1,A2 are used, but the calculation needs only
C  to remember if Ax.EQ.0 for gluon or Ax.NE.0 for quark or antiquark.
C  As we exit these loops, the Ax are correctly set for this purpose.
C
C  ----- Close IF that omits the case N1=N2 in process A.
      ENDIF
C  ----- Close DO for sum over flavors N2.
      ENDDO
C  ----- Close DO for sum over flavors N1.
      ENDDO
C
C  Get crossing sign
C
      IF (((AA.EQ.0).AND.(AB.NE.0)).OR.
     >    ((AB.EQ.0).AND.(AA.NE.0))     ) THEN
        SIGN = - 1.0D0
       ELSE
        SIGN =   1.0D0
      ENDIF
C
C  Calculate contribution from Born graph
C
      IF (SBORN) THEN
        INTEGRAND2TO2 = INTEGRAND2TO2 + D0 * L * SIGN * PSI4
      ENDIF
C
C  Calculate contribution from virtual graphs
C
      IF (SVIRTUAL) THEN
        INTEGRAND2TO2 = INTEGRAND2TO2 + D * L  * SIGN * PSI6NS
      ENDIF
C
C  Calculate contribution from term A
C
      IF(SACOLL) THEN
C
        INTEGRAND2TO2 = INTEGRAND2TO2
     >  + D * SIGN * PSI4
     >    *( - L * LOGMUASQ
     >        *( GAMMA(AA) -  COLOR(AA) *LOGXA  )
     >  + LEFFA )
C
      ENDIF
C
      IF(SASOFT) THEN
        INTEGRAND2TO2 = INTEGRAND2TO2
     >   + D * L * SIGN
     >    *( PSICAB * JAB + PSICA1 * JA1 + PSICA2 * JA2)
C
      ENDIF
C
C  Calculate contribution from term B
C
      IF(SBCOLL) THEN
C
        INTEGRAND2TO2 = INTEGRAND2TO2
     >  + D * SIGN * PSI4
     >    *( - L * LOGMUBSQ
     >        *( GAMMA(AB) -  COLOR(AB) * LOGXB  )
     >  + LEFFB )
C
      ENDIF
C
      IF(SBSOFT) THEN
        INTEGRAND2TO2 = INTEGRAND2TO2
     >   + D * L * SIGN
     >    *( PSICAB * JBA + PSICB1 * JB1 + PSICB2 * JB2)
C
      ENDIF
C
C
C  Calculate contribution from term 1
C
C
      IF(S1COLL) THEN
C
        INTEGRAND2TO2 = INTEGRAND2TO2
     >  + D * L * SIGN * PSI4
     >    *( - LOG16P2SQ
     >           *( GAMMA(A1) - 2.0D0 * COLOR(A1) * LOG2  )
     >       + GAMMAPRIME(A1) - 2.0D0 * COLOR(A1) * LOG2**2 )
C
      ENDIF
C
      IF(S1SOFT) THEN
        INTEGRAND2TO2 = INTEGRAND2TO2
     >   + D * L * SIGN
     >    *( PSICA1 * J1A + PSICB1 * J1B + PSIC12 * J12)
C
      ENDIF
C
C
C  Calculate contribution from term 2
C
      IF(S2COLL) THEN
C
        INTEGRAND2TO2 = INTEGRAND2TO2
     >  + D * L * SIGN * PSI4
     >    *( - LOG16P2SQ
     >           *( GAMMA(A2) - 2.0D0 * COLOR(A2) * LOG2  )
     >       + GAMMAPRIME(A2) - 2.0D0 * COLOR(A2) * LOG2**2 )
C
      ENDIF
C
      IF(S2SOFT) THEN
        INTEGRAND2TO2 = INTEGRAND2TO2
     >   + D * L * SIGN
     >    *( PSICA2 * J2A + PSICB2 * J2B + PSIC12 * J21)
C
      ENDIF
C
C  ---- Done!
C
C  ----- Close DO for sum over permutations.
      ENDDO
C  ----- Close IF that omits process if switch is off.
      ENDIF
C  ----- Close DO for sum over processes A,B,C,D.
      ENDDO
C
      INTEGRAND2TO2 = INTEGRAND2TO2*JACOBIAN*SVALUE
      CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND2TO2)
      RETURN
      END
C
C***********************************************************************
C  --------------------------------------------
C  Functions for Born and Virtual Contributions
C  --------------------------------------------
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION PSITILDE4(PROCESS,K)
C 
      INTEGER PROCESS
      DOUBLE PRECISION K(4,4)
C
C      This function gives the psitilde^(4) functions from E.S.
C
      DOUBLE PRECISION S,T,U
      DOUBLE PRECISION N,V
      PARAMETER ( V = 8.0 D0)
      PARAMETER ( N = 3.0 D0)
C
      S = 2.0D0 * K(1,2)
      T = 2.0D0 * K(1,3)
      U = 2.0D0 * K(1,4)
C
      IF (PROCESS.EQ.1) THEN
       PSITILDE4 = 2.0D0 * V * (S**2 + U**2) / T**2
      ELSEIF (PROCESS.EQ.2) THEN
       PSITILDE4 = 2.0D0 * V * (S**2 + U**2) / T**2
     >           + 2.0D0 * V * (S**2 + T**2) / U**2
     >           - 4.0D0 * V / N * S**2 / U / T
      ELSEIF (PROCESS.EQ.3) THEN
       PSITILDE4 = 2.0D0 * V / N
     >             *( V / U / T - 2.0D0 * N**2 / S**2 )
     >             *( T**2 + U**2 )
      ELSEIF (PROCESS.EQ.4) THEN
       PSITILDE4 = 4.0D0 * V * N**2
     >          *( S**2 + T**2 + U**2 )
     >          /S**2 /T**2 /U**2
     >          *( S**4 + T**4 + U**4 )
      ELSE
       PSITILDE4 = 0.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION PSITILDE6NS(PROCESS,K)
C 
      INTEGER PROCESS
      DOUBLE PRECISION K(4,4)
C
C      This function gives the psitilde^(6)_NS functions calculated
C      by subtracting the singular pieces from from E.S.
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
      DOUBLE PRECISION MUUV,MUA,MUB
      COMMON /MUVALS/ MUUV,MUA,MUB
      DOUBLE PRECISION QES
      COMMON /KEITHSQ/ QES
      DOUBLE PRECISION PI,N,V
      PARAMETER ( PI = 3.141592654 D0)
      PARAMETER ( V = 8.0 D0)
      PARAMETER ( N = 3.0 D0)
      DOUBLE PRECISION S,T,U
      DOUBLE PRECISION THETA
      DOUBLE PRECISION LS,LT,LU,LMU,L2S,L2T,L2U
C
      S = 2.0D0 * K(1,2)
      T = 2.0D0 * K(1,3)
      U = 2.0D0 * K(1,4)
C
C  Here are the Ellis and Sexton log functions
C
      LS = DLOG(DABS(S/QES**2))
      LT = DLOG(DABS(T/QES**2))
      LU = DLOG(DABS(U/QES**2))
      LMU = DLOG(DABS(MUUV**2/QES**2))
      L2S = LS**2 - PI**2 * THETA(S.GT.0.0D0)
      L2T = LT**2 - PI**2 * THETA(T.GT.0.0D0)
      L2U = LU**2 - PI**2 * THETA(U.GT.0.0D0)
C
C  Find what type, then calculate
C
      IF (PROCESS.EQ.1) THEN
       PSITILDE6NS =
     >    V/(9.0D0*N*T**2)
     >    *(-  9*PI**2*(N**2-4 )*(S**2 - U**2)
     >      +  2*(72 + 13*N**2 - 10*N*NFLAVOR + 9*N**2*PI**2)
     >          *(S**2 + U**2)
     >      +  6*N*(11*N - 2*NFLAVOR)*(S**2 + U**2)*LMU
     >      - 36*T*U*LS
     >      +  3*(- 6*S**2 - 30*U**2 + 4*N*NFLAVOR*(S**2 + U**2)
     >                     - N**2*(7*S**2 + 3*T**2 + U**2))*LT
     >      - 18*S*T*(N**2 -2)*LU
     >      - 36*(3*S**2 + U**2)*LS*LT
     >      - 18*(N**2 -2)*(S**2 + 3*U**2)*LT*LU
     >      + 18*(S**2 - U**2)*L2S
     >      +  9*(N**2*(S**2 + 3*U**2) + 2*(3*S**2 - U**2))*L2T
     >      -  9*(N**2 - 2)*(S**2 - U**2)*L2U)
      ELSEIF (PROCESS.EQ.2) THEN
C
       PSITILDE6NS =
     >     40*N*NFLAVOR*S**2*T*U - 4*N**2*S**2*T*U*(13 + 9*PI**2)
     >   - 9*T*U*(32*S**2 + PI**2*(S**2 + 2*T*U))
     >   + 36*N*(S**2*(4 + PI**2) * (T**2 + U**2)
     >   + (4 - PI**2)  *  (T**4 + U**4))
     >   + N**3*(S**2*(26 + 9*PI**2)*(T**2 + U**2)
     >   + (26 + 27*PI**2)*(T**4 + U**4))
     >   - 20*N**2*NFLAVOR*(T**4 + U**4 + S**2*(T**2 + U**2))
C
       PSITILDE6NS =  PSITILDE6NS
     >  + 6*N*(11*N - 2*NFLAVOR)
     >     *(-2*S**2*T*U + N*(T**4 + U**4 + S**2*(T**2 + U**2))) *LMU
     >  - 36*N*T*U*(T**2 + U**2)*LS
     >  - 6*U*(- 2*N**2*S**2*T + 2*N*NFLAVOR*S**2*T
     >         - 2*N**2*NFLAVOR*U*(S**2 + U**2)
     >         + 6*S*T*(T + 2*U) + N**3*(-3*T**3 + 2*T**2*U
     >         + 7*T*U**2 + 4*U**3)
     >         + 3*N*(2*T**3 + 3*T**2*U + 2*T*U**2 + 6*U**3)) *LT
     >  - 6*T*(- 2*N**2*S**2*U + 2*N*NFLAVOR*S**2*U
     >         - 2*N**2*NFLAVOR*T*(S**2 + T**2)
     >         + 6*S*U*(2*T + U)
     >         + N**3*(4*T**3 + 7*T**2*U + 2*T*U**2 - 3*U**3)
     >         + 3*N*(6*T**3 + 2*T**2*U + 3*T*U**2 + 2*U**3)) * LU
C
       PSITILDE6NS =  PSITILDE6NS
     >  - 36*U*(- S**2*T - N**2*S**2*T + N*U*(3*S**2 + U**2)) *LS*LT
     >  - 36*T*(- S**2*U - N**2*S**2*U + N*T*(3*S**2 + T**2)) *LS*LU
     >  - 18*(  2*N*(-2 + N**2)*(T**2 - T*U + U**2)
     >          *(2*T**2 + 3*T*U + 2*U**2)
     >        + T*U*(3*T**2 + 4*T*U + 3*U**2)) * LT*LU
     >  + 18*N*T*U*(S**2 + T**2 + U**2) * L2S
     >  + 9*U*(- 2*N**2*S**2*T - 4*N*(S**3 - S*T**2 - T**3)
     >         - T*(3*T**2 + 8*T*U + 3*U**2)
     >         - 2*N**3*(T**3 - T*U**2 - 2*U**3)) * L2T
     >  + 9*T*(- 2*N**2*S**2*U - 4*N*(S**3 - S*U**2 - U**3)
     >         - U*(3*T**2 + 8*T*U + 3*U**2)
     >         + 2*N**3*(2*T**3 + T**2*U - U**3)) * L2U
C
       PSITILDE6NS = V/(9*N**2*T**2*U**2) * PSITILDE6NS
C
      ELSEIF (PROCESS.EQ.3) THEN
C
       PSITILDE6NS =
     >   3*S**2*( -7*(T**2 + U**2) + PI**2*(3*S**2 -2*T*U) )
     >   + 3*N**2*( 14*(T**4 + U**4) + T*U*(13*S**2 + 4*T*U )
     >            + PI**2*(T**3*U + 2*T**2*U**2 + T*U**3)    )
     >   - 3*N**4*( 7*(T**4 + U**4)+ T*U*(S**2 + 10*T*U)
     >            - PI**2*(T - U)**2*(T**2 + U**2)  )
     >   + 2*N*(11*N - 2*NFLAVOR)*(-S**2 + N**2*(T**2 + U**2))
     >         *(T**2 + U**2) * LMU
     >   - 3*T*U*(  4*T**2 + 8*T*U + 4*U**2
     >         + N**2*(-1 + N**2)*(T**2 - 10*T*U + U**2) ) * LS
     >   + 3*S*U*(S + N**2*U)*(2*T + 3*U + N**2*(2*T - 3*U)) * LT
     >   + 3*S*T*(S + N**2*T)*(3*T + 2*U + N**2*(-3*T + 2*U)) * LU
C
       PSITILDE6NS = PSITILDE6NS 
     >   - 6*(  N**4*U*(U - T)*(2*T**2 + T*U + U**2)
     >     + S*(S - N**2*T)*(2*T**2 + 2*T*U + U**2)  ) * LS*LT
     >   - 6*(  N**4*T*(T - U)*(T**2 + T*U + 2*U**2)
     >   + S*(S - N**2*U)*(T**2 + 2*T*U + 2*U**2)  ) * LS*LU
     >   + 12*N**2*S**2*(T**2 + U**2) * LT*LU
     >   + 3*( 2*S**4 - 2*N**4*T*U*(T**2 + U**2)
     >     + N**2*(T**2 + T*U + 2*U**2)*(2*T**2 + T*U + U**2) ) *L2S
     >   - 3*S*(- N**4*U*(2*T**2 - T*U + U**2)
     >       + (N**2*T - S)*(2*T**2 + 2*T*U + U**2) ) *L2T
     >   - 3*S*(- N**4*T*(2*U**2 - T*U + T**2 )
     >       + (N**2*U - S)*( 2*U**2 + 2*T*U + T**2)) *L2U
C
       PSITILDE6NS = V/(3.0D0*N**2*S**2*T*U) * PSITILDE6NS
C
      ELSEIF (PROCESS.EQ.4) THEN
C
       PSITILDE6NS =
     >  2*N**2*NFLAVOR*( - (66 + 27*PI**2)*S**2*T**2*U**2
     >           + 40*(S**6 + T**6 + U**6))
     >   + 2*N**3*(  6*(125 - 27*PI**2)*S**2*T**2*U**2
     >         - 4*(67 - 9*PI**2)*(S**6 + T**6 + U**6))
     >   + 6*N**2*(11*N - 2*NFLAVOR)*(S**2 + T**2 + U**2)**3 *LMU
     >   + 6*N**2*T*U*(S**2 + T**2 + U**2)
     >       *(  NFLAVOR *(5*T**2 + 2*T*U + 5*U**2)
     >         - 2*N*(7*T**2 - 8*T*U + 7*U**2) ) *LS
C
       PSITILDE6NS = PSITILDE6NS 
     >   + 6*N**2*S*U*(S**2 + T**2 + U**2)
     >       *(  NFLAVOR *(5*S**2 + 2*S*U + 5*U**2)
     >         - 2*N*(7*S**2 - 8*S*U + 7*U**2) ) *LT
     >   + 6*N**2*S*T*(S**2 + T**2 + U**2)
     >       *(  NFLAVOR* (5*S**2 + 2*S*T + 5*T**2)
     >         - 2*N*(7*S**2 - 8*S*T + 7*T**2) ) *LU
     >   - 36*N**2*U**2*( NFLAVOR*S*T*(S**2 + T**2)
     >             + 2*N*(2*S**4 + 2*S**3*T + 3*S**2*T**2 + 2*S*T**3
     >             + 2*T**4))*LS*LT
C
       PSITILDE6NS = PSITILDE6NS 
     >   - 36*N**2*S**2*(  NFLAVOR*T*U*(T**2 + U**2)
     >             + 2*N*(2*T**4 + 2*T**3*U + 3*T**2*U**2 + 2*T*U**3
     >             + 2*U**4))*LT*LU
     >   - 36*N**2*T**2*(  NFLAVOR*S*U*(S**2 + U**2)
     >             + 2*N*(2*S**4 + 2*S**3*U + 3*S**2*U**2 + 2*S*U**3
     >             + 2*U**4))*LS*LU
     >   + 18*N**2*S**2*T*U*(4*N*(T**2 + U**2)
     >             - NFLAVOR*(T**2 + 3*T*U + U**2))*L2S
     >   + 18*N**2*S*T**2*U*(4*N*(S**2 + U**2)
     >             - NFLAVOR*(S**2 + 3*S*U + U**2))*L2T
     >   + 18*N**2*S*T*U**2*(4*N*(S**2 + T**2)
     >             - NFLAVOR*(S**2 + 3*S*T + T**2))*L2U
C
       PSITILDE6NS = V/(9.0D0*S**2*T**2*U**2) * PSITILDE6NS
C
      ELSE
       PSITILDE6NS = 0.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C  -------------------------------------
C  Functions for Collinear Contributions
C  -------------------------------------
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION EFFPARTON(A,X,SCALE)
C
      INTEGER A
      DOUBLE PRECISION X,SCALE
C
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
C
      DOUBLE PRECISION RANDOM,Z,DZDX1
      INTEGER APRIME
      DOUBLE PRECISION FCOLL
      DOUBLE PRECISION PARTON,COLOR
      DOUBLE PRECISION ALTARELLI,ALTRLIPRIME
      DOUBLE PRECISION FACTOR,SUBTRACTION
C
C  Integrates FCOLL(z) from X to 1, Monte Carlo style,
C  if called many times.
C
      Z = X + (1.0D0 - X) * RANDOM(1.0)
      DZDX1 = 1.0D0 - X
C------
C  Calculate FCOLL(Z), the integrand for the 2 -> 2 coll pieces.
C
      FACTOR = 2.0D0 * DLOG(RTS*X*(1.0D0 - Z)/SCALE/Z)
      SUBTRACTION = FACTOR * 2.0D0 * COLOR(A)
     >              * PARTON(A,X,SCALE) /Z/(1.0D0 - Z)
C
      FCOLL = - SUBTRACTION
C
      IF (A.NE.0) THEN
       FCOLL =
     >   FCOLL
     >  + (FACTOR * ALTARELLI(A,A,Z) - ALTRLIPRIME(A,A,Z))
     >      * PARTON(A,X/Z,SCALE) /Z
     >  + (FACTOR * ALTARELLI(A,0,Z)- ALTRLIPRIME(A,0,Z))
     >      * PARTON(0,X/Z,SCALE) /Z
      ELSE
       DO APRIME = - NFL,+NFL
        FCOLL =
     >   FCOLL
     >  + (FACTOR * ALTARELLI(A,APRIME,Z)- ALTRLIPRIME(A,APRIME,Z))
     >      * PARTON(APRIME,X/Z,SCALE) /Z
C
       ENDDO
      ENDIF
C------
      EFFPARTON = FCOLL * DZDX1
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION ALTARELLI(AOUT,AIN,Z)
C 
      INTEGER AOUT,AIN
      DOUBLE PRECISION Z
      DOUBLE PRECISION V,N,TWON,CF
      PARAMETER( V = 8.0D0 )
      PARAMETER( N = 3.0D0 )
      PARAMETER( TWON = 6.0D0 )
      PARAMETER( CF = 4.0D0/3.0D0 )
C
C   The Altarelli-Parisi kernel P_{AOUT/AIN}(z)
C
      IF ( (Z.GT.1.0D0).OR.(Z.LT.0.0D0) ) THEN
        ALTARELLI = 0.0D0
        RETURN
      ENDIF
C
      IF ( (AOUT.EQ.0).AND.(AIN.EQ.0) ) THEN
        ALTARELLI = TWON
     >   *( (1.0D0 - Z) /Z + Z /(1.0D0 - Z) + Z*(1.0D0 - Z) )
      ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.AOUT) ) THEN
        ALTARELLI = CF *(1.0D0 + Z**2) /(1.0D0 - Z)
      ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.0) ) THEN
        ALTARELLI = 0.5D0 * ( Z**2 + (1.0D0 - Z)**2 )
      ELSE IF ( (AOUT.EQ.0).AND.(AIN.NE.0) ) THEN
        ALTARELLI = CF *( 1.0D0 + (1.0D0 - Z)**2 ) /Z
      ELSE
        ALTARELLI = 0.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION ALTRLIPRIME(AOUT,AIN,Z)
C 
      INTEGER AOUT,AIN
      DOUBLE PRECISION Z
      DOUBLE PRECISION V,N,TWON,CF
      PARAMETER( V = 8.0D0 )
      PARAMETER( N = 3.0D0 )
      PARAMETER( TWON = 6.0D0 )
      PARAMETER( CF = 4.0D0/3.0D0)
C
C   The derivative of the Altarelli-Parisi kernel
C   P_{AOUT/AIN}(z,\epsilon) w.r.t. \epsilon (\epsilon = 0)
C
      IF ( (Z.GT.1.0D0).OR.(Z.LT.0.0D0) ) THEN
        ALTRLIPRIME = 0.0D0
        RETURN
      ENDIF
C
      IF ( (AOUT.EQ.0).AND.(AIN.EQ.0) ) THEN
        ALTRLIPRIME = 0.0D0
      ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.AOUT) ) THEN
        ALTRLIPRIME = - CF * (1.0D0 - Z)
      ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.0) ) THEN
        ALTRLIPRIME = 0.5D0 * ( Z**2 + (1.0D0 - Z)**2 - 1.0D0 )
      ELSE IF ( (AOUT.EQ.0).AND.(AIN.NE.0) ) THEN
        ALTRLIPRIME = - CF * Z
      ELSE
        ALTRLIPRIME = 0.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION GAMMA(NPARTON)
C 
      INTEGER NPARTON
C
C   This function gives the A-P delta function coefficient for
C   parton NPARTON: 2 for quarks and antiquarks,
C   11/2 - Nfl/3 for gluons
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
C
      IF (NPARTON.EQ.0) THEN
        GAMMA = 11.0D0/2.0D0 - NFLAVOR/3.0D0
      ELSE
        GAMMA = 2.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION GAMMAPRIME(NPARTON)
C 
      INTEGER NPARTON
C
C   This function gives the remainder from the collinear integration
C   for the coll-1 and coll-2 contributions for parton NPARTON
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
C
      DOUBLE PRECISION COLOR
      DOUBLE PRECISION PI,LOG2,V,N
      PARAMETER ( PI = 3.141592654 D0)
      PARAMETER ( LOG2 = 0.693147181 D0)
      PARAMETER ( V = 8.0 D0)
      PARAMETER ( N = 3.0 D0)
C
      IF (NPARTON.EQ.0) THEN
        GAMMAPRIME = 67.0D0*N/9.0D0 - 23.0D0/18.0D0 * NFLAVOR
     >             -2.0D0 * PI**2 * COLOR(NPARTON) /3.0D0
      ELSE
        GAMMAPRIME = 13.0D0*V/(4.0D0*N)
     >             -2.0D0 * PI**2 * COLOR(NPARTON) /3.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C  --------------------------------
C  Functions for Soft Contributions
C  --------------------------------
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION PSICTILDE(PROCESS,I1,I2,K)
C 
      INTEGER PROCESS
      INTEGER I1,I2
      DOUBLE PRECISION K(4,4)
C
C      This function gives the lambdatilde(a1,a2,s,t,u) functions
C
      INTEGER II1,II2
      DOUBLE PRECISION S,T,U
      DOUBLE PRECISION HAT,HAU,HB,HCT,HCU,HCS,HDS,HDT,HDU
      DOUBLE PRECISION N,V,CF
      PARAMETER( N = 3.0D0 )
      PARAMETER( V = 8.0D0 )
      PARAMETER( CF = 4.0D0 / 3.0D0 )
C
      S = 2.0D0 * K(1,2)
      T = 2.0D0 * K(1,3)
      U = 2.0D0 * K(1,4)
C
C  Make it a symmetric function of I1,I2
C
      IF (I1.LT.I2) THEN
        II1 = I1
        II2 = I2
      ELSE
        II1 = I2
        II2 = I1
      ENDIF
C
      IF (PROCESS.EQ.1) THEN
        HAT = 2.0D0 * (S**2 + U**2) /T**2
        IF     ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN
          PSICTILDE = 4.0D0 * CF * HAT
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN
          PSICTILDE = - 2.0D0 * CF * HAT
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 2.0D0 * CF * (N**2 - 2.0D0) * HAT
        ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 4.0D0 * CF * HAT
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN
          PSICTILDE = - 2.0D0 * CF* HAT
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN
          PSICTILDE = 2.0D0 * CF * (N**2 - 2.0D0) * HAT
        ENDIF
C
      ELSEIF (PROCESS.EQ.2) THEN
        HAT = 2.0D0 * (S**2 + U**2) /T**2
        HAU = 2.0D0 * (S**2 + T**2) /U**2
        HB  = 2.0D0 * S**2 /T /U
        IF     ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN
          PSICTILDE = 4.0D0 * CF /N
     >       *( - HB + N*(HAT + HAU) - N**2 * HB )
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN
          PSICTILDE = 4.0D0 * CF /N
     >       *( HB - N*(HAU + 0.5D0*HAT) + 0.5D0 * N**3 * HAU )
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 4.0D0 * CF /N
     >       *( HB - N*(HAT + 0.5D0*HAU) + 0.5D0 * N**3 * HAT )
        ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 4.0D0 * CF /N
     >       *( - HB + N*(HAT + HAU) - N**2 * HB )
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 4.0D0 * CF /N
     >       *( HB - N*(HAU + 0.5D0*HAT) + 0.5D0 * N**3 * HAU )
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN
          PSICTILDE = 4.0D0 * CF /N
     >       *( HB - N*(HAT + 0.5D0*HAU) + 0.5D0 * N**3 * HAT )
        ENDIF
C
      ELSEIF (PROCESS.EQ.3) THEN
        HCT = 2.0D0 * (T**2 + U**2) /S**2 * T/U
        HCU = 2.0D0 * (T**2 + U**2) /S**2 * U/T
        HCS = 2.0D0 * (T**2 + U**2) /T /U
        IF     ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN
          PSICTILDE = V * ( - HCT - HCU + HCS + HCS /N**2 )
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN
          PSICTILDE = V * ( N**2 * HCU - HCS )
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN
          PSICTILDE = V * ( N**2 * HCT - HCS )
        ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN
          PSICTILDE =  V * N**2 * ( HCT + HCU )
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN
          PSICTILDE = V * ( N**2 * HCU - HCS )
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN
          PSICTILDE = V * ( N**2 * HCT - HCS )
        ENDIF
C
      ELSEIF (PROCESS.EQ.4) THEN
       HDS = 2.0D0 * (T**2 + U**2) /S**2 /T**2 /U**2
     >          * (S**4 + T**4 + U**4)
       HDT = 2.0D0 * (U**2 + S**2) /S**2 /T**2 /U**2
     >          * (S**4 + T**4 + U**4)
       HDU = 2.0D0 * (S**2 + T**2) /S**2 /T**2 /U**2
     >          * (S**4 + T**4 + U**4)
        IF     ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN
          PSICTILDE = 2.0D0 * V * N**3 * HDS
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN
          PSICTILDE = 2.0D0 * V * N**3 * HDT
        ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 2.0D0 * V * N**3 * HDU
        ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 2.0D0 * V * N**3 * HDS
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN
          PSICTILDE = 2.0D0 * V * N**3 * HDT
        ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN
          PSICTILDE = 2.0D0 * V * N**3 * HDU
        ENDIF
C
      ELSE
       PSICTILDE = 0.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION LI2(X)
C 
      DOUBLE PRECISION X
C
C   Let ERR be the target fractional error
C
      DOUBLE PRECISION XN
      INTEGER N
      DOUBLE PRECISION ERR
      PARAMETER(ERR = 1.0D-6)
C
      DOUBLE PRECISION PI
      PARAMETER ( PI = 3.141592654 D0)
C
      XN = 1.0D0
      LI2 = 0.0D0
C
      DO N=1,100
        XN = XN * X
        IF (XN.LT.ERR) THEN
         RETURN
        ENDIF
        LI2 = LI2 + XN/N**2
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
C  End of functions for 1- and 2-dimensional integrals
C
C***********************************************************************
C***********************************************************************
C
C
      SUBROUTINE INTEGRATEA(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF)
C 
      DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,COEF
C
C   Integrates Term A. First compute Y3 and P3. (For Y3, set 
C   Y3 = ln(\infty) = 30 if W is zero or very near to zero.)
C   We sum over terms in the jet-defining function, labelled by CASE.
C   For each term, JETDEF returns OK = .TRUE. if the theta functions
C   for that term are satisfied and then returns the values of PJ and
C   YJ determined by the delta functions for that term.  We then set
C   MUUV and MUA,MUB accordingly and call the BINIT subroutine, which
C   multiplies by gaussian smearing for each PJ,YJ grid point and
C   adds the resulting integrand to the integral for that grid point.
C
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      DOUBLE PRECISION Y3,P3
      DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE
      DOUBLE PRECISION MUUV,MUA,MUB
      DOUBLE PRECISION FA,INTEGRAND
      LOGICAL OK
      INTEGER CASE
C
      INTEGER NJETDEF
      COMMON /DEFJET/ NJETDEF
C
      IF(W.GT.1.0D-12)THEN
        Y3 = DLOG(RTS/W)
      ELSE
        Y3 = 30
      ENDIF
C
      P3 = XI*W
C
      DO CASE = 1,NJETDEF
      CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,
     >            OK,JETVAR1,JETVAR2,JETVAR3,SVALUE)
      IF(OK) THEN
         CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUA,MUB)
         INTEGRAND = COEF * SVALUE 
     >              * FA(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB)
         CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND)
      ENDIF
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE INTEGRATEB(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF)
C 
      DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,COEF
C
C   Integrates Term A. First compute Y3 and P3. (For Y3, set 
C   Y3 = -ln(\infty) = -30 if W is zero or very near to zero.)
C   We sum over terms in the jet-defining function, labelled by CASE.
C   For each term, JETDEF returns OK = .TRUE. if the theta functions
C   for that term are satisfied and then returns the values of PJ and
C   YJ determined by the delta functions for that term.  We then set
C   MUUV and MUA,MUB accordingly and call the BINIT subroutine, which
C   multiplies by gaussian smearing for each PJ,YJ grid point and
C   adds the resulting integrand to the integral for that grid point.
C
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      DOUBLE PRECISION Y3,P3
      DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE
      DOUBLE PRECISION MUUV,MUA,MUB
      DOUBLE PRECISION FB,INTEGRAND
      LOGICAL OK
      INTEGER CASE
C
      INTEGER NJETDEF
      COMMON /DEFJET/ NJETDEF
C
      IF(W.GT.1.0D-12)THEN
        Y3 = - DLOG(RTS/W)
      ELSE
        Y3 = - 30
      ENDIF
C
      P3 = XI*W
C
      DO CASE = 1,NJETDEF
      CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,
     >            OK,JETVAR1,JETVAR2,JETVAR3,SVALUE)
      IF(OK) THEN
         CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUA,MUB)
         INTEGRAND = COEF * SVALUE
     >              * FB(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB)
         CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND)
      ENDIF
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF)
C 
      DOUBLE PRECISION Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF
C
C   Integrates Term 1. First we compute new variables P2 and PHI2.
C   We sum over terms in the jet-defining function, labelled by CASE.
C   For each term, JETDEF returns OK = .TRUE. if the theta functions
C   for that term are satisfied and then returns the values of PJ and
C   YJ determined by the delta functions for that term.  We then set
C   MUUV and MUA,MUB accordingly and call the BINIT subroutine, which
C   multiplies by gaussian smearing for each PJ,YJ grid point and
C   adds the resulting integrand to the integral for that grid point.
C
      DOUBLE PRECISION P2X,P2Y,P2,PHI2
      DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE
      DOUBLE PRECISION MUUV,MUA,MUB
      DOUBLE PRECISION F1,INTEGRAND
      LOGICAL OK
      INTEGER CASE
C
      INTEGER NJETDEF
      COMMON /DEFJET/ NJETDEF
C
      P2X = - P1 * DCOS(PHI1) - P3 * DCOS(PHI3)
      P2Y = - P1 * DSIN(PHI1) - P3 * DSIN(PHI3)
      P2 = DSQRT(P2X**2 + P2Y**2)
      PHI2 = DATAN2(P2Y,P2X)
C
      DO CASE = 1,NJETDEF
      CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,
     >            OK,JETVAR1,JETVAR2,JETVAR3,SVALUE)
      IF(OK) THEN
         CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUA,MUB)
         INTEGRAND = COEF * SVALUE
     >               * F1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUA,MUB)
         CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND)
      ENDIF
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF)
C 
      DOUBLE PRECISION Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF
C
C   Integrates Term 2.
C   We sum over terms in the jet-defining function, labelled by CASE.
C   For each term, JETDEF returns OK = .TRUE. if the theta functions
C   for that term are satisfied and then returns the values of PJ and
C   YJ determined by the delta functions for that term.  We then set
C   MUUV and MUA,MUB accordingly and call the BINIT subroutine, which
C   multiplies by gaussian smearing for each PJ,YJ grid point and
C   adds the resulting integrand to the integral for that grid point. 
C
      DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE
      DOUBLE PRECISION MUUV,MUA,MUB
      DOUBLE PRECISION F2,INTEGRAND
      LOGICAL OK
      INTEGER CASE
C
      INTEGER NJETDEF
      COMMON /DEFJET/ NJETDEF
C
      DO CASE = 1,NJETDEF
      CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,
     >            OK,JETVAR1,JETVAR2,JETVAR3,SVALUE)
      IF(OK) THEN
         CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUA,MUB)
         INTEGRAND = COEF * SVALUE
     >               * F2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUA,MUB)
         CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND)
      ENDIF
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION
     >              FA(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB)
C 
      DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      LOGICAL STYPE(4)
      COMMON /SWTYPE/ STYPE
C
      INTEGER PERM(5),PERMINV(5)
      INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4)
      COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR
C
      DOUBLE PRECISION PIN(5,5),K(5,5)
      DOUBLE PRECISION PAB,PA1,PA2,PA3MOD,
     >                 PB1,PB2,PB3MOD,P12,P13MOD,P23MOD
      DOUBLE PRECISION P1,XA,XB
      DOUBLE PRECISION COS23,EY1,EY2,P3,FEXP
      DOUBLE PRECISION ALPHAS,FACOEF,FASUM,LUMINOSITY
      DOUBLE PRECISION LMNSTY(7)
      INTEGER I,J,PROCESS,NPERM
C
      DOUBLE PRECISION SIGN,CROSSINGSIGN
      DOUBLE PRECISION PSIAB,PSIA1,PSIA2,RESIDUE
C
      DOUBLE PRECISION PISQUARE
      PARAMETER (PISQUARE = 9.869604404D0)
C
C
      COS23 = DCOS(PHI2-PHI3)
      EY1 = FEXP(Y1)
      EY2 = FEXP(Y2)
      P3 = XI*W
      P1 = DSQRT( P2**2 + P3**2 + 2.0D0*P2*P3*COS23 )
      XA = (1/RTS)*( P1*EY1 + P2*EY2 + XI*RTS )
      XB = (1/RTS)*( P1/EY1 + P2/EY2 + XI * W**2 /RTS )
C
C  Test for particle identification conditions and XA,XB:
C
      IF (P3.GT.P1) THEN
        FA = 0.0D0
        RETURN
      ELSE IF (P3.GT.P2) THEN
        FA = 0.0D0
        RETURN
      ELSE IF (XA.GT.1.0D0) THEN
        FA = 0.0D0
        RETURN
      ELSE IF (XB.GT.1.0D0) THEN
        FA = 0.0D0
        RETURN
      ENDIF
C
C   Construct the dot products of the P_J vectors, J = A,B,1,2,3.
C   The 'MOD' variables are the dot products divided by XI.
C   This division is done algebraicly because XI may be zero.
C
      PAB =  XA*XB*RTS**2 /2.0D0
      PA1 =  XA*RTS*P1 /EY1 /2.0D0
      PA2 =  XA*RTS*P2 /EY2 /2.0D0
      PB1 =  XB*RTS*P1 *EY1 /2.0D0
      PB2 =  XB*RTS*P2 *EY2 /2.0D0
      P12 =  P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P2**2 + P2*P3*COS23
      PA3MOD   = XA * W**2 /2.0D0
      PB3MOD   = XB * RTS**2 /2.0D0
      P13MOD   = 0.5D0 * P1 * RTS /EY1 + 0.5D0 * P1 * W**2 * EY1 /RTS
     >                   + XI * W**2 + P2 * W * COS23
      P23MOD   = 0.5D0 * P2 * RTS /EY2 + 0.5D0 * P2 * W**2 * EY2 /RTS
     >                               - P2 * W * COS23
C
C  Construct the PIN matrix (incoming momenta)
C
      PIN(1,2) =  PAB
      PIN(1,3) = -PA1
      PIN(1,4) = -PA2
      PIN(1,5) = -XI*PA3MOD
      PIN(2,3) = -PB1
      PIN(2,4) = -PB2
      PIN(2,5) = -XI*PB3MOD
      PIN(3,4) =  P12
      PIN(3,5) =  XI*P13MOD
      PIN(4,5) =  XI*P23MOD
C
      DO I = 1,5
       DO J = 1,I-1
        PIN(I,J) = PIN(J,I)
       ENDDO
      ENDDO
C
C  Calculate coefficient
C
      FACOEF = ALPHAS(MUUV)**3 * P2 /(4.0D0 * PISQUARE * XA * RTS**4)
C
C  Generate LMNSTY array for the given XA,XB
C
      CALL PARTONSIN(LMNSTY,XA,XB,MUA,MUB)
C
C  Initialize FA
C
      FA = 0.0D0
C
C  Main loops over processes and permutations
C
      DO PROCESS = 1,4
C
      IF ( .NOT.STYPE(PROCESS) ) GOTO 300
C
      DO NPERM = 1,NPERMS(PROCESS)
C
C  Initialize for each process and permutation
C
      DO I = 1,5
         PERM(I) = PERMARR(I,NPERM,PROCESS)
         PERMINV(I) = PERMINVARR(I,NPERM,PROCESS)
      ENDDO
C
      DO I = 1,5
       DO J = 1,5
        K(I,J) = PIN(PERM(I),PERM(J))
       ENDDO
      ENDDO
C
C   Complete calculation
C
      SIGN = CROSSINGSIGN(PROCESS,PERMINV)
C
      PSIAB = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(2),K)
      PSIA1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(3),K)
      PSIA2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(4),K)
C
      FASUM =    PAB/(PA3MOD + PB3MOD) * PSIAB
     >         + PA1/(PA3MOD + P13MOD) * PSIA1
     >         + PA2/(PA3MOD + P23MOD) * PSIA2 
C
      FA = FA + FACOEF * FASUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV)
C
C  Close loop over NPERM
      ENDDO
C  Close loop over PROCESS
300   ENDDO
C
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION
     >              FB(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB)
C 
      DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUA,MUB
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      LOGICAL STYPE(4)
      COMMON /SWTYPE/ STYPE
C
 
      INTEGER PERM(5),PERMINV(5)
      INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4)
      COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR
C
      DOUBLE PRECISION PIN(5,5),K(5,5)
      DOUBLE PRECISION PAB,PA1,PA2,PA3MOD,
     >                 PB1,PB2,PB3MOD,P12,P13MOD,P23MOD
      DOUBLE PRECISION P1,XA,XB
      DOUBLE PRECISION COS23,EY1,EY2,P3,FEXP
      DOUBLE PRECISION ALPHAS,FBCOEF,FBSUM,LUMINOSITY
      DOUBLE PRECISION LMNSTY(7)
      INTEGER I,J,PROCESS,NPERM
C
      DOUBLE PRECISION SIGN,CROSSINGSIGN
      DOUBLE PRECISION PSIAB,PSIB1,PSIB2,RESIDUE
C
      DOUBLE PRECISION PISQUARE
      PARAMETER (PISQUARE = 9.869604404D0)
C
C
      COS23 = DCOS(PHI2-PHI3)
      EY1 = FEXP(Y1)
      EY2 = FEXP(Y2)
      P3 = XI*W
      P1 = DSQRT( P2**2 + P3**2 + 2.0D0*P2*P3*COS23 )
      XA = (1/RTS)*( P1*EY1 + P2*EY2 + XI * W**2 / RTS)
      XB = (1/RTS)*( P1/EY1 + P2/EY2 + XI*RTS )
C
C  Test for particle identification conditions and XA,XB:
C
      IF (P3.GT.P1) THEN
        FB = 0.0D0
        RETURN
      ELSE IF (P3.GT.P2) THEN
        FB = 0.0D0
        RETURN
      ELSE IF (XA.GT.1.0D0) THEN
        FB = 0.0D0
        RETURN
      ELSE IF (XB.GT.1.0D0) THEN
        FB = 0.0D0
        RETURN
      ENDIF
C
C   Construct the dot products of the P_J vectors, J = A,B,1,2,3.
C   The 'MOD' variables are the dot products divided by XI.
C   This division is done algebraicly because XI may be zero.
C
      PAB = XA*XB*RTS**2 /2.0D0
      PA1 = XA*RTS*P1 /EY1 /2.0D0
      PA2 = XA*RTS*P2 /EY2 /2.0D0
      PB1 = XB*RTS*P1* EY1 /2.0D0
      PB2 = XB*RTS*P2* EY2 /2.0D0
      P12 = P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P2**2 + P2*P3*COS23
      PA3MOD = XA * RTS**2 /2.0D0
      PB3MOD = XB * W**2 /2.0D0
      P13MOD = 0.5D0 * P1 * RTS * EY1 + 0.5D0 * P1* W**2 / EY1 / RTS
     >                   + XI * W**2 + P2 * W * COS23
      P23MOD = 0.5D0 * P2 * RTS * EY2 + 0.5D0 * P2 * W**2 / EY2 / RTS
     >                             - P2 * W * COS23
C
C  Construct the PIN matrix (incoming momenta):
C
      PIN(1,2) =  PAB
      PIN(1,3) = -PA1
      PIN(1,4) = -PA2
      PIN(1,5) = -XI*PA3MOD
      PIN(2,3) = -PB1
      PIN(2,4) = -PB2
      PIN(2,5) = -XI*PB3MOD
      PIN(3,4) =  P12
      PIN(3,5) =  XI*P13MOD
      PIN(4,5) =  XI*P23MOD
C
      DO I = 1,5
       DO J = 1,I-1
        PIN(I,J) = PIN(J,I)
       ENDDO
      ENDDO
C
C  Calculate coefficient for factor HB
C
      FBCOEF = ALPHAS(MUUV)**3 * P2 /(4.0D0 * PISQUARE * XB * RTS**4)
C
C  Generate LMNSTY array for the given XA,XB
C
      CALL PARTONSIN(LMNSTY,XA,XB,MUA,MUB)
C
C  Initialize FB
C
      FB = 0.0D0
C
C  Main loops over processes and permutations
C
      DO PROCESS = 1,4
C
      IF ( .NOT.STYPE(PROCESS) ) GOTO 300
C
      DO NPERM = 1,NPERMS(PROCESS)
C
C  Initialize for each process and permutation
C
      DO I = 1,5
         PERM(I) = PERMARR(I,NPERM,PROCESS)
         PERMINV(I) = PERMINVARR(I,NPERM,PROCESS)
      ENDDO
C
      DO I = 1,5
       DO J = 1,5
         K(I,J) = PIN(PERM(I),PERM(J))
       ENDDO
      ENDDO
C
C   Complete calculation:
C
      SIGN = CROSSINGSIGN(PROCESS,PERMINV)
C
      PSIAB = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(2),K)
      PSIB1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(3),K)
      PSIB2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(4),K)
C
      FBSUM =    PAB/(PA3MOD + PB3MOD) * PSIAB
     >         + PB1/(PB3MOD + P13MOD) * PSIB1
     >         + PB2/(PB3MOD + P23MOD) * PSIB2
C
      FB = FB + FBCOEF * FBSUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV)
C
C  Close loop over NPERM
      ENDDO
C  Close loop over PROCESS
300   ENDDO
C
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION
     >              F1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUA,MUB)
C 
      DOUBLE PRECISION Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUA,MUB
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      LOGICAL STYPE(4)
      COMMON /SWTYPE/ STYPE
C
      INTEGER PERM(5),PERMINV(5)
      INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4)
      COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR
C
      DOUBLE PRECISION PIN(5,5),K(5,5)
      DOUBLE PRECISION PAB,PA1,PA2,PA3MOD,
     >                 PB1,PB2,PB3MOD,P12,P13MOD,P23MOD
      DOUBLE PRECISION P2,XA,XB
      DOUBLE PRECISION COS13,EY1,EY2,EY3,FEXP
      DOUBLE PRECISION ALPHAS,F1COEF,F1SUM,LUMINOSITY
      DOUBLE PRECISION LMNSTY(7)
      INTEGER I,J,PROCESS,NPERM
C
      DOUBLE PRECISION SIGN,CROSSINGSIGN
      DOUBLE PRECISION PSIA1,PSIB1,PSI12,RESIDUE
C
      DOUBLE PRECISION PISQUARE
      PARAMETER (PISQUARE = 9.869604404D0)
C
      COS13 = DCOS(PHI1-PHI3)
      EY1 = FEXP(Y1)
      EY2 = FEXP(Y2)
      EY3 = FEXP(Y3)
      P2 = DSQRT( P1**2 + P3**2 + 2.0D0*P1*P3*COS13 )
      XA = (1/RTS)*( P1*EY1 + P2*EY2 + P3*EY3 )
      XB = (1/RTS)*( P1/EY1 + P2/EY2 + P3/EY3 )
C
C  Test for particle identification conditions and XA,XB:
C
      IF (P3.GT.P1) THEN
        F1 = 0.0D0
        RETURN
      ELSE IF (P3.GT.P2) THEN
        F1 = 0.0D0
        RETURN
      ELSE IF (XA.GT.1.0D0) THEN
        F1 = 0.0D0
        RETURN
      ELSE IF (XB.GT.1.0D0) THEN
        F1 = 0.0D0
        RETURN
      ENDIF
C
C   Construct the dot products of the P_J vectors, J = A,B,1,2,3.
C   The 'MOD' variables are the dot products divided by P3.
C   This division is done algebraicly because P3 may be zero.
C
      PAB    = XA*XB*(RTS)**2 /2.0D0
      PA1    = XA*RTS*P1 /EY1 /2.0D0
      PA2    = XA*RTS*P2 /EY2 /2.0D0
      PA3MOD = XA*RTS    /EY3 /2.0D0
      PB1    = XB*RTS*P1* EY1 /2.0D0
      PB2    = XB*RTS*P2* EY2 /2.0D0
      PB3MOD = XB*RTS   * EY3 /2.0D0
      P12    = P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P1**2 + P1*P3*COS13
      P13MOD = P1   *0.5D0*(EY1/EY3 + EY3/EY1) - P1   *COS13
      P23MOD = P2   *0.5D0*(EY2/EY3 + EY3/EY2) + P3    + P1   *COS13
C
C  Construct the PIN matrix (incoming momenta):
C
      PIN(1,2) =  PAB
      PIN(1,3) = -PA1
      PIN(1,4) = -PA2
      PIN(1,5) = -P3*PA3MOD
      PIN(2,3) = -PB1
      PIN(2,4) = -PB2
      PIN(2,5) = -P3*PB3MOD
      PIN(3,4) =  P12
      PIN(3,5) =  P3*P13MOD
      PIN(4,5) =  P3*P23MOD
C
      DO I = 1,5
       DO J = 1,I-1
        PIN(I,J) = PIN(J,I)
       ENDDO
      ENDDO
C
C  Calculate coefficient for factor H1
C
      F1COEF =  ALPHAS(MUUV)**3 /(8.0D0 * PISQUARE * RTS**4)
C
C  Generate LMNSTY array for the given XA,XB
C
      CALL PARTONSIN(LMNSTY,XA,XB,MUA,MUB)
C
C  Initialize psi variables
C
      F1 = 0.0D0
C
C  Main loops over processes and permutations
C
      DO PROCESS = 1,4
C
      IF ( .NOT.STYPE(PROCESS) ) GOTO 300
C
      DO NPERM = 1,NPERMS(PROCESS)
C
C  Initialize for each process and permutation
C
      DO I = 1,5
         PERM(I) = PERMARR(I,NPERM,PROCESS)
         PERMINV(I) = PERMINVARR(I,NPERM,PROCESS)
      ENDDO
C
      DO I = 1,5
       DO J = 1,5
        K(I,J) = PIN(PERM(I),PERM(J))
       ENDDO
      ENDDO
C
C   Complete calculation:
C
      SIGN = CROSSINGSIGN(PROCESS,PERMINV)
C
      PSIA1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(3),K)
      PSIB1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(3),K)
      PSI12 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(3),PERMINV(4),K)
C
      F1SUM =    PA1/(PA3MOD + P13MOD) * PSIA1
     >         + PB1/(PB3MOD + P13MOD) * PSIB1
     >         + P12/(P13MOD + P23MOD) * PSI12
C
      F1 = F1 + F1COEF * F1SUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV)
C
C  Close loop over NPERM
      ENDDO
C  Close loop over PROCESS
300   ENDDO
C
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION
     >              F2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUA,MUB)
C 
      DOUBLE PRECISION Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUA,MUB
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      LOGICAL STYPE(4)
      COMMON /SWTYPE/ STYPE
C
      INTEGER PERM(5),PERMINV(5)
      INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4)
      COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR
C
      DOUBLE PRECISION PIN(5,5),K(5,5)
      DOUBLE PRECISION PAB,PA1,PA2,PA3MOD,
     >                 PB1,PB2,PB3MOD,P12,P13MOD,P23MOD
      DOUBLE PRECISION P1,XA,XB
      DOUBLE PRECISION COS23,EY1,EY2,EY3,FEXP
      DOUBLE PRECISION ALPHAS,F2COEF,F2SUM,LUMINOSITY
      DOUBLE PRECISION LMNSTY(7)
      INTEGER I,J,PROCESS,NPERM
C
      DOUBLE PRECISION SIGN,CROSSINGSIGN
      DOUBLE PRECISION PSIA2,PSIB2,PSI12,RESIDUE
C
      DOUBLE PRECISION PISQUARE
      PARAMETER (PISQUARE = 9.869604404D0)
C
      COS23 = DCOS(PHI2-PHI3)
      EY1 = FEXP(Y1)
      EY2 = FEXP(Y2)
      EY3 = FEXP(Y3)
      P1 = DSQRT( P2**2 + P3**2 + 2.0D0*P2*P3*COS23 )
      XA = (1/RTS)*( P1*EY1 + P2*EY2 + P3*EY3 )
      XB = (1/RTS)*( P1/EY1 + P2/EY2 + P3/EY3 )
C
C  Test for particle identification conditions and XA,XB:
C
      IF (P3.GT.P1) THEN
        F2 = 0.0D0
        RETURN
      ELSE IF (P3.GT.P2) THEN
        F2 = 0.0D0
        RETURN
      ELSE IF (XA.GT.1.0D0) THEN
        F2 = 0.0D0
        RETURN
      ELSE IF (XB.GT.1.0D0) THEN
        F2 = 0.0D0
        RETURN
      ENDIF
C
C   Construct the dot products of the P_J vectors, J = A,B,1,2,3.
C   The 'MOD' variables are the dot products divided by P3.
C   This division is done algebraicly because P3 may be zero.
C
      PAB    = XA*XB*(RTS)**2 /2.0D0
      PA1    = XA*RTS*P1 /EY1 /2.0D0
      PA2    = XA*RTS*P2 /EY2 /2.0D0
      PA3MOD = XA*RTS    /EY3 /2.0D0
      PB1    = XB*RTS*P1* EY1 /2.0D0
      PB2    = XB*RTS*P2* EY2 /2.0D0
      PB3MOD = XB*RTS   * EY3 /2.0D0
      P12    = P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P2**2 + P2*P3*COS23
      P13MOD = P1   *0.5D0*(EY1/EY3 + EY3/EY1) + P3    + P2   *COS23
      P23MOD = P2   *0.5D0*(EY2/EY3 + EY3/EY2) - P2  *COS23
C
C  Construct the PIN matrix (incoming momenta):
C
      PIN(1,2) =  PAB
      PIN(1,3) = -PA1
      PIN(1,4) = -PA2
      PIN(1,5) = -P3*PA3MOD
      PIN(2,3) = -PB1
      PIN(2,4) = -PB2
      PIN(2,5) = -P3*PB3MOD
      PIN(3,4) =  P12
      PIN(3,5) =  P3*P13MOD
      PIN(4,5) =  P3*P23MOD
C
      DO I = 1,5
       DO J = 1,I-1
        PIN(I,J) = PIN(J,I)
       ENDDO
      ENDDO
C
C  Calculate coefficient for factor H2
C
      F2COEF = ALPHAS(MUUV)**3 /(8.0D0 * PISQUARE * RTS**4)
C
C  Generate LMNSTY array for the given XA,XB
C
      CALL PARTONSIN(LMNSTY,XA,XB,MUA,MUB)
C
C  Initialize F2
C
      F2 = 0.0D0
C
C  Main loops over processes and permutations
C
      DO PROCESS = 1,4
C
      IF ( .NOT.STYPE(PROCESS) ) GOTO 300
C
      DO NPERM = 1,NPERMS(PROCESS)
C
C  Initialize for each process and permutation
C
      DO I = 1,5
         PERM(I) = PERMARR(I,NPERM,PROCESS)
         PERMINV(I) = PERMINVARR(I,NPERM,PROCESS)
      ENDDO
C
      DO I = 1,5
       DO J = 1,5
        K(I,J) = PIN(PERM(I),PERM(J))
       ENDDO
      ENDDO
C
C   Complete calculation:
C
      SIGN = CROSSINGSIGN(PROCESS,PERMINV)
C
      PSIA2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(4),K)
      PSIB2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(4),K)
      PSI12 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(3),PERMINV(4),K)
C
      F2SUM =    PA2/(PA3MOD + P23MOD) * PSIA2
     >         + PB2/(PB3MOD + P23MOD) * PSIB2
     >         + P12/(P13MOD + P23MOD) * PSI12
C
      F2 = F2 + F2COEF * F2SUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV)
C
C  Close loop over NPERM
      ENDDO
C  Close loop over PROCESS
300   ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION CROSSINGSIGN(PROCESS,PERMINV)
C 
      INTEGER PROCESS,PERMINV(5)
C
C    This function calculated the crossing sign needed in the
C    the four dimensional integrals.  SIGN is -1 if one of
C    particles A and B is a gluon while the other is not.  Otherwise
C    SIGN = +1.
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
C
      INTEGER A,B
C
C   Start with a plus sign, then change it if necessary:
C
      CROSSINGSIGN = 1.0D0
      A = PERMINV(1)
      B = PERMINV(2)
C
C-----------------------------
C
      IF ((PROCESS.EQ.1).OR.(PROCESS.EQ.2)) THEN
C
        IF ((A.EQ.5).OR.(B.EQ.5)) THEN
          CROSSINGSIGN = - CROSSINGSIGN
        ENDIF
C
      ELSE IF (PROCESS.EQ.3) THEN
C
        IF ((A.EQ.1).OR.(A.EQ.2)) THEN
          CROSSINGSIGN = - CROSSINGSIGN
        ENDIF
        IF ((B.EQ.1).OR.(B.EQ.2)) THEN
          CROSSINGSIGN = - CROSSINGSIGN
        ENDIF
C
      ENDIF
C
C------------------------
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION CONVERT(PHI)
      DOUBLE PRECISION PHI
C
C  Gets PHI + 2*pi*N such that the result lies between - pi and pi.
C  Here IDNINT( x ) is the nearest integer to x.
C
      DOUBLE PRECISION TWOPI
      PARAMETER ( TWOPI = 6.283185308 D0)
C
      CONVERT = PHI - TWOPI * IDNINT( PHI/TWOPI )
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION RESIDUE(PROCESS,J,N,M,KIN)
C 
      INTEGER PROCESS,J,N,M
      DOUBLE PRECISION KIN(5,5)
C
C   This function calculates the function tilde-Psi^j_{nm},
C   which is the coefficient of (k_n*k_m) / (k_n*k_j)(k_m*k_j))
C   in the expansion of tilde-\Psi(K) for the PROCESS.
C
C   In the computation,only certain standard values of j,n,m are used.
C   Then tilde-\Psi for other values of j,n,m are obtained by
C   using symmetries.  The relation is
C
C   Psi^j_{n m}(KIN) = Psi^jj_{nn mm}(K)
C
C   with j = PERM(jj), n = PERM(nn), m = PERM(nm)
C   KIN(I,J) = K(PERM(I),PERM(J)).
C   The proper value of JJ for process A is given in JA(J,N,M).
C   The proper value of NN,MM for process A is in SUBSCRIPTA(J,N,M).
C   as a two digit integer SUBSCRIPTA = 10*min(NN,MM) + max(NN,MM).
C   (Recall that  Psi^j_{n m} = Psi^j_{m n} ).
C   The proper value of PERM(I) for process A is in PERMA(I,J,N,M).
C   Simarly for processes B,C,D.
C   These values are put into COMMON /RESIDUEHELP/  by the
C   subroutine RESIDUEINIT at the start of the program.
C
      INTEGER I
      INTEGER JJ,SUBSCRIPT,PERM(5)
C
      INTEGER JA(5,5,5),JB(5,5,5),JC(5,5,5),JD(5,5,5)
      INTEGER SUBSCRIPTA(5,5,5),SUBSCRIPTB(5,5,5)
      INTEGER SUBSCRIPTC(5,5,5),SUBSCRIPTD(5,5,5)
      INTEGER PERMA(5,5,5,5),PERMB(5,5,5,5)
      INTEGER PERMC(5,5,5,5),PERMD(5,5,5,5)
      COMMON /RESIDUEHELP/ JA,JB,JC,JD,
     > SUBSCRIPTA,SUBSCRIPTB,SUBSCRIPTC,SUBSCRIPTD,
     > PERMA,PERMB,PERMC,PERMD
C
C  Variables for the calculation of all processes
C
      DOUBLE PRECISION K12,K13,K14,K15,K23,K24,K25,K34,K35,K45
      DOUBLE PRECISION NCOLOR,V,CF,TWOCF,FOURCF
      PARAMETER( NCOLOR = 3.0D0)
      PARAMETER( V = 8.0D0)
      PARAMETER( CF = V/2.0D0/NCOLOR)
      PARAMETER(TWOCF = 2.0D0 * CF)
      PARAMETER(FOURCF = 4.0D0 * CF)
      DOUBLE PRECISION NCF,TWONCF,FOURNCF2,FOURCF2
      PARAMETER(NCF = NCOLOR * CF)
      PARAMETER(TWONCF = 2.0D0 * NCOLOR * CF)
      PARAMETER(FOURNCF2 = 4.0D0 * NCOLOR * CF**2)
      PARAMETER(FOURCF2 = 4.0D0 * CF**2)
C
C  Variables for process A:
C
      DOUBLE PRECISION COEF512,COEF513,COEF514,FCTOR
      PARAMETER(COEF512 =   4.0D0 * CF)
      PARAMETER(COEF513 = - 2.0D0 * CF)
      PARAMETER(COEF514 = 2.0D0 * CF * ( 2.0D0 * NCOLOR * CF - 1.0D0 ))
      DOUBLE PRECISION Q1,Q2,Q3
C
C  Variables for process B:
C
      DOUBLE PRECISION FCTOR1,FCTOR2,FCTOR3,COMBO,COMBO2
      DOUBLE PRECISION H1,H2,H3,H4,H5,H6,H7
      DOUBLE PRECISION SFACTOR,TFACTOR,UFACTOR,STUFACTOR
      DOUBLE PRECISION ONELESS2NCF,TWOPLUS2NCF
      PARAMETER(ONELESS2NCF = 1.0D0 - 2.0D0*NCOLOR*CF)
      PARAMETER(TWOPLUS2NCF = 2.0D0 + 2.0D0*NCOLOR*CF)
C
C  Variables for process C:
C
      DOUBLE PRECISION VOVERN2,ONEPLUSN2,N2,N2OVER3,N2COMBO
      PARAMETER(VOVERN2 = V/NCOLOR**2)
      PARAMETER(ONEPLUSN2 = 1.0D0 + NCOLOR**2)
      PARAMETER(N2 = NCOLOR**2)
      PARAMETER(N2OVER3 = NCOLOR**2/3.0D0)
      PARAMETER(N2COMBO = (1.0D0 + NCOLOR**2)/NCOLOR**2)
      DOUBLE PRECISION H123,H124
      DOUBLE PRECISION U123,U124,V124
      DOUBLE PRECISION W123,W124,W125
      DOUBLE PRECISION S1314,S1525,S1535,S1545
      DOUBLE PRECISION S2324,S2334,S2434,S2535,S2545,S3545
      DOUBLE PRECISION N2OVERK12
      DOUBLE PRECISION X123,X124,X125,Y123,Y124,Y125
      DOUBLE PRECISION L2345,L3435,L3445,L3545
C
C  Variables for process D:
C
      DOUBLE PRECISION SUM4,DENOMS,TWOVN3
      PARAMETER(TWOVN3 = 2.0D0 * V * NCOLOR**3)
C
C--------------
C
C  Program must change this!
C
      RESIDUE = 1.0D15
C
      IF (PROCESS.EQ.1) THEN
C
        JJ = JA(J,N,M)
        SUBSCRIPT = SUBSCRIPTA(J,N,M)
C
        DO I = 1,5
         PERM(I) = PERMA(I,J,N,M)
        ENDDO
C
      IF ( JJ.EQ.5 ) THEN
C
      IF ( SUBSCRIPT.EQ.12 ) THEN
C
C  The function tilde-Psi^5_{12} for process A:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          FCTOR = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24)
          RESIDUE = COEF512 * FCTOR
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.13 ) THEN
C
C  The function tilde-Psi^5_{13} for process A:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          FCTOR = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24)
          RESIDUE = COEF513 * FCTOR
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.14 ) THEN
C
C  The function tilde-Psi^5_{14} for process A:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          FCTOR = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24)
          RESIDUE = COEF514 * FCTOR
C--
          RETURN
C
C  End SUBSCRIPT test, continue JJ test
C
      ENDIF
      ELSEIF ( JJ.EQ.1 ) THEN
C
      IF ( SUBSCRIPT.EQ.35 ) THEN
C
C  The function tilde-Psi^1_{35} for process A:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          Q1 = K14*K23 - K12*K34
          Q2 = K13*K24
          Q3 = K12**2 + K14**2 + K23**2 + K34**2
          RESIDUE = FOURCF * Q3 /(K24 * K25 * K35**2 * K45)
     >             *( -  K14*K23*(K13 + K24)
     >                -  K13*K24*(K14 + K23)
     >                + (K12 + K34)*(Q1 + Q2)/2.0D0
     >                + NCF * (  K14*(K12*K13 + K24*K34)
     >                         + K23*(K12*K24 + K13*K34)
     >                         - (K14 + K23)*(Q1 - Q2)   ))
C--
          RETURN
C
      ELSEIF((SUBSCRIPT.EQ.23).OR.(SUBSCRIPT.EQ.24).OR.
     >       (SUBSCRIPT.EQ.25).OR.(SUBSCRIPT.EQ.34).OR.
     >       (SUBSCRIPT.EQ.45) ) THEN
C
C  The function tilde-Psi^1_{23} for process A:
C  The function tilde-Psi^1_{24} for process A:
C  The function tilde-Psi^1_{25} for process A:
C  The function tilde-Psi^1_{34} for process A:
C  The function tilde-Psi^1_{45} for process A:
C--
          RESIDUE = 0.0D0
C--
          RETURN
C
C  End  SUBSCRIPT test and JJ test
C
       ENDIF
       ENDIF
C
C--------------
C
      ELSE IF (PROCESS.EQ.2) THEN
C
        JJ = JB(J,N,M)
        SUBSCRIPT = SUBSCRIPTB(J,N,M)
C
        DO I = 1,5
         PERM(I) = PERMB(I,J,N,M)
        ENDDO
C
      IF ( JJ.EQ.5 ) THEN
C
      IF ( SUBSCRIPT.EQ.12 ) THEN
C
C  The function tilde-Psi^5_{12} for process B:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          FCTOR1 = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24)
          FCTOR2 = (K12**2 + K13**2 + K24**2 + K34**2)/(K14*K23)
          FCTOR3 = (K12**2 + K34**2)
     >             *(K12*K34 - K14*K23 - K13*K24)/(K13*K14*K23*K24)
          COMBO2 = FOURCF * (FCTOR1 +  FCTOR2 -  FCTOR3 /NCOLOR)
          RESIDUE = COMBO2 - FOURCF2 * FCTOR3
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.13 ) THEN
C
C  The function tilde-Psi^5_{13} for process B:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          FCTOR1 = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24)
          FCTOR2 = (K12**2 + K13**2 + K24**2 + K34**2)/(K14*K23)
          FCTOR3 = (K12**2 + K34**2)
     >             *(K12*K34 - K14*K23 - K13*K24)/(K13*K14*K23*K24)
          COMBO = TWOCF * (FCTOR1 +  FCTOR2 -  FCTOR3 /NCOLOR)
          RESIDUE = - COMBO + FOURNCF2 * FCTOR2
C--
        RETURN
C
C  End SUBSCRIPT test, continue JJ test
C
      ENDIF
      ELSEIF ( JJ.EQ.1 ) THEN
      IF ( SUBSCRIPT.EQ.23 ) THEN
C
C  The function tilde-Psi^1_{23} for process B:
C--
          RESIDUE = 0.0D0
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.34 ) THEN
C
C  The function tilde-Psi^1_{34} for process B:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          H2 = K24*K35
          H3 = K23*K45
          H4 = K35*K45
          H5 = K25*K34
          H1 = 2.0D0*H5 - H2 - H3
          SFACTOR = K12**2 + K34**2
          TFACTOR = K13**2 + K24**2
          UFACTOR = K14**2 + K23**2
          STUFACTOR = K12*K34 - K14*K23 - K13*K24
C
          RESIDUE = TWOCF/(H5*H4)
     >      *( ( 2.0D0*H4 - H1  + TWONCF*(H4-H5) )/(K23*K24*NCOLOR)
     >           *SFACTOR*STUFACTOR
     >        + (K13/K23)*(SFACTOR + TFACTOR)*(H1 + TWONCF*H2)
     >        + (K14/K24)*(SFACTOR + UFACTOR)*(H1 + TWONCF*H3) )
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.25 ) THEN
C
C  The function tilde-Psi^1_{25} for process B:
C--
          RESIDUE = 0.0D0
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.35 ) THEN
C
C  The function tilde-Psi^1_{35} for process B:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          H6 = K24*K45
          H7 = K24*K25
          SFACTOR = K12**2 + K34**2
          UFACTOR = K14**2 + K23**2
          STUFACTOR = K12*K34 - K14*K23 - K13*K24
C
          RESIDUE = ( TWOCF/K35)
     >    *((SFACTOR + UFACTOR)
     >        *( (2.0D0*K12 - K13*K25/K35)/H7 - K14*ONELESS2NCF/H6 )
     >    + SFACTOR*STUFACTOR/(K23*NCOLOR)*(1.0D0/H6 + TWOPLUS2NCF/H7))
C--
          RETURN
C
C  End  SUBSCRIPT test and JJ test
C
       ENDIF
       ENDIF
C
C--------------
C
      ELSE IF (PROCESS.EQ.3) THEN
C
        JJ = JC(J,N,M)
        SUBSCRIPT = SUBSCRIPTC(J,N,M)
C
        DO I = 1,5
         PERM(I) = PERMC(I,J,N,M)
        ENDDO
C
      IF ( JJ.EQ.5 ) THEN
C
      IF ( SUBSCRIPT.EQ.12 ) THEN
C
C  The function tilde-Psi^5_{12} for process C:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          H123 = K13**2 + K23**2
          H124 = K14**2 + K24**2
C
          RESIDUE = VOVERN2
     >     *( ONEPLUSN2*(H123/(K14*K24) + H124/(K13*K23))
     >       - (N2/K34/K12)
     >        *( H123*(K13/K14 + K23/K24)+ H124*(K24/K23 + K14/K13)))
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.13 ) THEN
C
C  The function tilde-Psi^5_{13} for process C:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          H123 = K13**2 + K23**2
          H124 = K14**2 + K24**2
          U123 = K13**2 + K23**2
          U124 = K14**2 + K24**2
          V124 =  K14*K24
C
          RESIDUE = V/K13
     >     *( - U124/K23 - K13*U123/V124
     >        + N2/(K12*K24*K34)
     >           *( V124*H124 + K13*K23*H123) )
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.34 ) THEN
C
C  The function tilde-Psi^5_{34} for process C:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          W123 = K13**2 + K23**2
          W124 = K14**2 + K24**2
          W125 = K15**2 + K25**2
          S1314 = K13*K14
          S1525 = K15*K25
          S1535 = K15*K35
          S1545 = K15*K45
          S2324 = K23*K24
          S2334 = K23*K34
          S2434 = K24*K34
          S2535 = K25*K35
          S2545 = K25*K45
          S3545 = K35*K45
          N2OVERK12 = N2/K12
C
          RESIDUE = V/K34 *( - W125
     >    *( (S1535/K23 + S1545/K24 )/S1314
     >     + (S2535/S2324 + S3545/S2434)/K13
     >     + (S2545/S2324 + S3545/S2334)/K14
     >     -  S3545*ONEPLUSN2 /(S1314*S2324*N2OVERK12)
     >     - N2OVERK12
     >        *(  (S1525/K23 + S1545/K34)/K14
     >          + (S1525/K24 + S1535/K34)/K13
     >          + S2535/S2334 + S2545/S2434 ) )
     >      + N2OVERK12
     >        *(  W123*(K13/K14 + K23/K24)
     >          + W124*(K14/K13 + K24/K23) )  )
C--
          RETURN
C
C  End SUBSCRIPT test, continue JJ test
C
      ENDIF
      ELSEIF ( JJ.EQ.1 ) THEN
C
      IF ( SUBSCRIPT.EQ.23 ) THEN
C
C  The function tilde-Psi^1_{23} for process C:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          X123 = K13**2 + K23**2
          X124 = K14**2 + K24**2
          X125 = K15**2 + K25**2
          Y123 = K23 * X123
          Y124 = K24 * X124
          Y125 = K25 * X125
          L2345 = K23 * K45
          L3435 = K34 * K35
          L3445 = K34 * K45
          L3545 = K35 * K45
C
          RESIDUE =  V/K23
     >       *(  N2*(K25*K34 + K24*K35)/(L3435*K45*K24*K25)
     >             *( K15*Y125 + K14*Y124)
     >          - K12/K24/K25/L2345*(  Y125*(L2345/K34 + K25)
     >                               + Y124*(L2345/K35 + K24)  )
     >          + K13*N2OVER3
     >               *(  (Y124/L3545 + Y125/L3445)/K23
     >                 + (Y125/L3435 + Y123/L3545)/K24
     >                 + (Y123/L3445 + Y124/L3435)/K25 ))
 
C--
          RETURN
C
      ELSEIF ( SUBSCRIPT.EQ.34 ) THEN
C
C  The function tilde-Psi^1_{34} for process C:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          RESIDUE = V*(K15**2 + K25**2)/K34/K23/K24
     >               *( K12*N2COMBO - K15*(K23/K35 + K24/K45))
C--
          RETURN
C
C  End  SUBSCRIPT test and JJ test
C
       ENDIF
       ENDIF
C--------------
C
      ELSE IF (PROCESS.EQ.4) THEN
C
        JJ = JD(J,N,M)
        SUBSCRIPT = SUBSCRIPTD(J,N,M)
C
        DO I = 1,5
         PERM(I) = PERMD(I,J,N,M)
        ENDDO
C
      IF ( JJ.EQ.5 ) THEN
C
      IF ( SUBSCRIPT.EQ.34 ) THEN
C
C  The function tilde-Psi^5_{34} for process D:
C
          CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C--
          SUM4 = K12**4 + K13**4 + K14**4 + K15**4
     >                  + K23**4 + K24**4 + K25**4
     >                           + K34**4 + K35**4
     >                                    + K45**4
C
          DENOMS = (1.0D0 /K13/K24 + 1.0D0 /K23/K14)/K12/K34
C
          RESIDUE = TWOVN3 * SUM4 * DENOMS
C--
          RETURN
C
C  End  SUBSCRIPT test and JJ test
C
       ENDIF
       ENDIF
C
C--------------
C
C  End PROCESS test
C
      ENDIF
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE PERMUTE
     >   (KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM)
C 
      DOUBLE PRECISION KIN(5,5)
      DOUBLE PRECISION K12,K13,K14,K15,K23,K24,K25,K34,K35,K45
      INTEGER PERM(5)
C
C  Produces K(I,J) = KIN(PERM(I),PERM(J))
C
      K12 = KIN( PERM(1),PERM(2) )
      K13 = KIN( PERM(1),PERM(3) )
      K14 = KIN( PERM(1),PERM(4) )
      K15 = KIN( PERM(1),PERM(5) )
      K23 = KIN( PERM(2),PERM(3) )
      K24 = KIN( PERM(2),PERM(4) )
      K25 = KIN( PERM(2),PERM(5) )
      K34 = KIN( PERM(3),PERM(4) )
      K35 = KIN( PERM(3),PERM(5) )
      K45 = KIN( PERM(4),PERM(5) )
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE RESIDUEINIT
C 
C
C   This is to initialize arrays for the function RESIDUE,
C   which calculates the function tilde-Psi^j_{nm},
C   which is the coefficient of (k_n*k_m) / (k_n*k_j)(k_m*k_j))
C   in the expansion of tilde-\Psi(K) for the PROCESS.
C
C   In the computation,only certain standard values of j,n,m are used.
C   Then tilde-\Psi for other values of j,n,m are obtained by
C   using symmetries.  The relation is
C
C   Psi^j_{n m}(KIN) = Psi^jj_{nn mm}(K)
C
C   with jj = PERM(j), nn = PERM(n), mm = PERM(n)
C   K(I,J) = KIN(PERM(I),PERM(J)).
C   The proper value of JJ for process A is given in JA(J,N,M).
C   The proper value of NN,MM for process A is in SUBSCRIPTA(J,N,M)
C   as a two digit integer SUBSCRIPTA = 10*min(NN,MM) + max(NN,MM).
C   (Recall that  Psi^j_{n m} = Psi^j_{m n} ).
C   The proper value of PERM(I) for process A is in PERMA(I,J,N,M).
C   Simarly for processes B,C,D.
C   These values are put into COMMON /RESIDUEHELP/  by the
C   subroutine RESIDUEINIT at the start of the program.
C
C
      INTEGER JA(5,5,5),JB(5,5,5),JC(5,5,5),JD(5,5,5)
      INTEGER SUBSCRIPTA(5,5,5),SUBSCRIPTB(5,5,5)
      INTEGER SUBSCRIPTC(5,5,5),SUBSCRIPTD(5,5,5)
      INTEGER PERMA(5,5,5,5),PERMB(5,5,5,5)
      INTEGER PERMC(5,5,5,5),PERMD(5,5,5,5)
      COMMON /RESIDUEHELP/ JA,JB,JC,JD,
     > SUBSCRIPTA,SUBSCRIPTB,SUBSCRIPTC,SUBSCRIPTD,
     > PERMA,PERMB,PERMC,PERMD
C
      INTEGER PROCESS,J,N,M
      LOGICAL MATCH
      INTEGER I,L,LL
      INTEGER A(4,5),B(8,5),C(12,5)
      DATA    A/1,2,3,4,
     >          2,1,4,3,
     >          3,4,1,2,
     >          4,3,2,1,
     >          5,5,5,5/
      DATA    B/1,1,2,2,3,3,4,4,
     >          2,2,1,1,4,4,3,3,
     >          3,4,3,4,1,2,1,2,
     >          4,3,4,3,2,1,2,1,
     >          5,5,5,5,5,5,5,5/
      DATA    C/1,1,1,1,1,1,2,2,2,2,2,2,
     >          2,2,2,2,2,2,1,1,1,1,1,1,
     >          3,3,4,4,5,5,3,3,4,4,5,5,
     >          4,5,3,5,3,4,4,5,3,5,3,4,
     >          5,4,5,3,4,3,5,4,5,3,4,3/
C
C--------------
C
      DO J = 1,5
      DO N = 1,5
      DO M = 1,5
C
      IF((N.NE.J).AND.(M.NE.J).AND.(M.NE.N) ) THEN
C
      PROCESS=1
C
      DO L = 1,4
C
      IF      (MATCH(J,N,M,A(L,5),A(L,1),A(L,2))) THEN
C
      JA(J,N,M) = 5
      SUBSCRIPTA(J,N,M) = 12
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,5),A(L,1),A(L,3))) THEN
C
      JA(J,N,M) = 5
      SUBSCRIPTA(J,N,M) = 13
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,5),A(L,1),A(L,4))) THEN
C
      JA(J,N,M) = 5
      SUBSCRIPTA(J,N,M) = 14
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,1),A(L,2),A(L,3))) THEN
C
      JA(J,N,M) = 1
      SUBSCRIPTA(J,N,M) = 23
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,1),A(L,2),A(L,4))) THEN
C
      JA(J,N,M) = 1
      SUBSCRIPTA(J,N,M) = 24
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,1),A(L,3),A(L,4))) THEN
C
      JA(J,N,M) = 1
      SUBSCRIPTA(J,N,M) = 34
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,1),A(L,2),A(L,5))) THEN
C
      JA(J,N,M) = 1
      SUBSCRIPTA(J,N,M) = 25
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,1),A(L,3),A(L,5))) THEN
C
      JA(J,N,M) = 1
      SUBSCRIPTA(J,N,M) = 35
      GOTO 101
C
      ELSE IF (MATCH(J,N,M,A(L,1),A(L,4),A(L,5))) THEN
C
      JA(J,N,M) = 1
      SUBSCRIPTA(J,N,M) = 45
      GOTO 101
C
C  If none of the tests was successful, try a different L
C
      ENDIF
      ENDDO
C
C  This gets executed when test is successful for this L
C
101   DO I = 1,5
         PERMA(I,J,N,M) = A(L,I)
      ENDDO
C
C
C--------------
C
      PROCESS=2
C
      DO L = 1,8
C
      IF      (MATCH(J,N,M,B(L,5),B(L,1),B(L,2))) THEN
C
      JB(J,N,M) = 5
      SUBSCRIPTB(J,N,M) = 12
      GOTO 102
C
      ELSE IF (MATCH(J,N,M,B(L,5),B(L,1),B(L,3))) THEN
C
      JB(J,N,M) = 5
      SUBSCRIPTB(J,N,M) = 13
      GOTO 102
C
      ELSE IF (MATCH(J,N,M,B(L,1),B(L,2),B(L,3))) THEN
C
      JB(J,N,M) = 1
      SUBSCRIPTB(J,N,M) = 23
      GOTO 102
C
      ELSE IF (MATCH(J,N,M,B(L,1),B(L,3),B(L,4))) THEN
C
      JB(J,N,M) = 1
      SUBSCRIPTB(J,N,M) = 34
      GOTO 102
C
      ELSE IF (MATCH(J,N,M,B(L,1),B(L,2),B(L,5))) THEN
C
      JB(J,N,M) = 1
      SUBSCRIPTB(J,N,M) = 25
      GOTO 102
C
      ELSE IF (MATCH(J,N,M,B(L,1),B(L,3),B(L,5))) THEN
C
      JB(J,N,M) = 1
      SUBSCRIPTB(J,N,M) = 35
      GOTO 102
C
C  If none of the tests was successful, try a different L
C
      ENDIF
      ENDDO
C
C  This gets executed when test is successful for this L
C
102   DO I = 1,5
         PERMB(I,J,N,M) = B(L,I)
      ENDDO
C
C--------------
C
      PROCESS=3
C
      DO L = 1,12
C
      IF      (MATCH(J,N,M,C(L,5),C(L,1),C(L,2))) THEN
C
      JC(J,N,M) = 5
      SUBSCRIPTC(J,N,M) = 12
      GOTO 103
C
      ELSE IF (MATCH(J,N,M,C(L,5),C(L,1),C(L,3))) THEN
C
      JC(J,N,M) = 5
      SUBSCRIPTC(J,N,M) = 13
      GOTO 103
C
      ELSE IF (MATCH(J,N,M,C(L,5),C(L,3),C(L,4))) THEN
C
      JC(J,N,M) = 5
      SUBSCRIPTC(J,N,M) = 34
      GOTO 103
C
      ELSE IF (MATCH(J,N,M,C(L,1),C(L,2),C(L,3))) THEN
C
      JC(J,N,M) = 1
      SUBSCRIPTC(J,N,M) = 23
      GOTO 103
C
      ELSE IF (MATCH(J,N,M,C(L,1),C(L,3),C(L,4))) THEN
C
      JC(J,N,M) = 1
      SUBSCRIPTC(J,N,M) = 34
      GOTO 103
C
C  If none of the tests was successful, try a different L
C
      ENDIF
      ENDDO
C
C  This gets executed when test is successful for this L
C
103   DO I = 1,5
         PERMC(I,J,N,M) = C(L,I)
      ENDDO
C
C--------------
C
      PROCESS=4
C
      JD(J,N,M) = 5
      SUBSCRIPTD(J,N,M) = 34
C
        PERMD(5,J,N,M) = J
        PERMD(4,J,N,M) = M
        PERMD(3,J,N,M) = N
C
C  Now assign PERMD(2)=LL, PERMD(1)=L where L and LL are chosen
C  from the remaining numbers (i.e. not J,N,M) with L < LL.
C  When done, get out of the Do loops with a GOTO.
C
        DO L = 1,4
          IF ((L.NE.J).AND.(L.NE.N).AND.(L.NE.M)) THEN
            PERMD(1,J,N,M) = L
            DO LL = L+1,5
               IF ((LL.NE.J).AND.(LL.NE.N).AND.(LL.NE.M)) THEN
                  PERMD(2,J,N,M) = LL
                  GOTO 702
               ENDIF
            ENDDO
          ENDIF
        ENDDO
C
C
C  End of loops over J,N,M and ENDIF for (N.NE.J),(M.NE.J),(M.NE.N)
C
C
702   ENDIF
      ENDDO
      ENDDO
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C
      LOGICAL FUNCTION MATCH(J,N,M,JJ,NN,MM)
C 
      INTEGER J,N,M,JJ,NN,MM
C
      IF (J.EQ.JJ) THEN
        IF      ( (N.EQ.NN).AND.(M.EQ.MM) ) THEN
          MATCH = .TRUE.
        ELSE IF ( (N.EQ.MM).AND.(M.EQ.NN) ) THEN
          MATCH = .TRUE.
        ELSE
          MATCH = .FALSE.
        ENDIF
      ELSE
        MATCH = .FALSE.
      ENDIF
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION LUMINOSITY(LMNSTY,PROCESS,PERMINV)
C 
      DOUBLE PRECISION LMNSTY(7)
      INTEGER PROCESS,PERMINV(5)
C
C    This function calculated the luminosity function needed in the
C    the four dimensional integrals.  LUMINOSITY is the sum over
C    quark flavors of L(\vec a_A, \vec a_B, x_A, x_B), where the
C    sum runs over flavors a_A and a_B of the initial state particles
C    and over flavors a_1, a_2, a_3 of the final state particles,
C    consistent with PROCESS and PERMINV. For each case, the charge
C    conjugate luninosity is included also. When the sum over final
C    state particles is nontrivial, it produces a factor NFLAVOR -1
C    (for process A) or 2*NFLAVOR (for process C).
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
C
      INTEGER A,B
C
C   Initialize:
C
      A = PERMINV(1)
      B = PERMINV(2)
C
C-----------------------------
C
      IF (PROCESS.EQ.1) THEN
C
C  Possible cases are (A,B) = (1,2),(1,3),(1,4),(1,5),(5,1).
C
       IF (A.EQ.1) THEN
C
        IF      (B.EQ.2) THEN
C
C  Luminosity for (PROCESS;A,B) = (1;1,2) i.e. (A,B) = (q,Q):
C  We use WNUM(101) to get the counting factor for a quark or antiquark.
C
          LUMINOSITY = LMNSTY(1)
C
        ELSE IF (B.EQ.3) THEN
C
C  Luminosity for (PROCESS;A,B) = (1;1,3) i.e. (A,B) = (q,qbar):
C
          LUMINOSITY = LMNSTY(4)
C
          LUMINOSITY = (NFLAVOR - 1.0D0) * LUMINOSITY
C
        ELSE IF (B.EQ.4) THEN
C
C  Luminosity for (PROCESS;A,B) = (1;1,4) i.e. (A,B) = (q,Qbar):
C
          LUMINOSITY = LMNSTY(2)
C
        ELSE IF (B.EQ.5) THEN
C
C  Luminosity for (PROCESS;A,B) = (1;1,5) i.e. (A,B) = (q,g):
C
          LUMINOSITY = LMNSTY(5)
C
          LUMINOSITY = (NFLAVOR - 1.0D0) * LUMINOSITY
C
        ENDIF
C
       ELSE IF ((A.EQ.5).AND.(B.EQ.1)) THEN
C
C  Luminosity for (PROCESS;A,B) = (1;5,1) i.e. (A,B) = (g,q):
C
          LUMINOSITY = LMNSTY(6)
C
          LUMINOSITY = (NFLAVOR - 1.0D0) * LUMINOSITY
C
       ENDIF
C
C------------------
C
      ELSE IF (PROCESS.EQ.2) THEN
C
C  Possible cases are (A,B) = (1,2),(1,3),(1,5),(5,1).
C
       IF (A.EQ.1) THEN
C
        IF      (B.EQ.2) THEN
C
C  Luminosity for (PROCESS;A,B) = (2;1,2) i.e. (A,B) = (q,q):
C
          LUMINOSITY = LMNSTY(3)
C
        ELSE IF (B.EQ.3) THEN
C
C  Luminosity for (PROCESS;A,B) = (2;1,3) i.e. (A,B) = (q,qbar):
C
          LUMINOSITY = LMNSTY(4)
C
        ELSE IF (B.EQ.5) THEN
C
C  Luminosity for (PROCESS;A,B) = (2;1,5) i.e. (A,B) = (q,g):
C
          LUMINOSITY = LMNSTY(5)
C
        ENDIF
C
       ELSE IF ((A.EQ.5).AND.(B.EQ.1)) THEN
C
C  Luminosity for (PROCESS;A,B) = (2;5,1) i.e. (A,B) = (g,q):
C
          LUMINOSITY = LMNSTY(6)
C
       ENDIF
C
C------------------
C
      ELSE IF (PROCESS.EQ.3) THEN
C
C  Possible cases are (A,B) = (1,2),(1,3),(3,1),(3,4).
C
       IF (A.EQ.1) THEN
C
        IF      (B.EQ.2) THEN
C
C  Luminosity for (PROCESS;A,B) = (3;1,2) i.e. (A,B) = (q,qbar):
C
          LUMINOSITY = LMNSTY(4)
C
        ELSE IF (B.EQ.3) THEN
C
C  Luminosity for (PROCESS;A,B) = (3;1,3) i.e. (A,B) = (q,g):
C
          LUMINOSITY = LMNSTY(5)
C
        ENDIF
C
       ELSE IF(A.EQ.3) THEN
C
        IF      (B.EQ.1) THEN
C
C  Luminosity for (PROCESS;A,B) = (3;3,1) i.e. (A,B) = (g,q):
C
          LUMINOSITY = LMNSTY(6)
C
        ELSE IF (B.EQ.4) THEN
C
C  Luminosity for (PROCESS;A,B) = (3;3,4) i.e. (A,B) = (g,g):
C
          LUMINOSITY = LMNSTY(7)
C
          LUMINOSITY = 2.0D0 * NFLAVOR * LUMINOSITY
C
        ENDIF
C
       ENDIF
C
C------------------------
C
      ELSE IF (PROCESS.EQ.4) THEN
C
C  Only case is (A,B) = (g,g).
C
          LUMINOSITY = LMNSTY(7)
C
C------------------------
C
      ENDIF
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE PARTONSIN(LMNSTY,XA,XB,MUA,MUB)
C 
      DOUBLE PRECISION LMNSTY(7)
      DOUBLE PRECISION XA,XB,MUA,MUB
C
C  This function computes the product, LMNSTY(KIND), of a parton
C  distribution function for A and a parton distribution function
C  for B, summed over flavors according to KIND:
C       A            B
C  1 = quark     + quark      (different flavors)
C      antiquark + antiquark  (different flavors)
C
C  2 = quark     + antiquark  (different flavors)
C      antiquark + quark      (different flavors)
C
C  3 = quark     + quark      (same flavor)
C      antiquark + antiquark  (same flavor)
C
C  4 = quark     + antiquark  (same flavor)
C      antiquark + quark      (same flavor)
C
C  5 = quark     + gluon
C      antiquark + gluon
C
C  6 = gluon     + quark
C      gluon     + antiquark
C
C  7 = gluon     + gluon
C
C
      INTEGER I,J,KIND
      DOUBLE PRECISION PARTON,PRTNS,WNUM
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
C  Switch for p-pbar collisions (if true) or p-p collisions (if false)
      LOGICAL PPBAR
      COMMON /BEAMS/ PPBAR
C
      DOUBLE PRECISION PRTNA(-6:+6),PRTNB(-6:+6)
C
C  For p-pbar collisions, replace partons -> antipartons for hadron B.
C
      IF (PPBAR) THEN
       DO I = -NFL,+NFL
        PRTNA(I) = PARTON( I,XA,MUA)/XA/WNUM(I)
        PRTNB(I) = PARTON(-I,XB,MUB)/XB/WNUM(I)
       ENDDO
      ELSE
       DO I = -NFL,+NFL
        PRTNA(I) = PARTON( I,XA,MUA)/XA/WNUM(I)
        PRTNB(I) = PARTON( I,XB,MUB)/XB/WNUM(I)
       ENDDO
      ENDIF
C
      DO KIND = 1,7
C
       PRTNS = 0.0D0
C
      IF (KIND.EQ.1) THEN
C                          i.e. (A,B) = (q,Q):
C
          DO I = 1,NFL
           DO J = 1,NFL
            IF (J.NE.I) THEN
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB( J)
            ENDIF
           ENDDO
          ENDDO
C
          DO I = -NFL,-1
           DO J = -NFL,-1
            IF (J.NE.I) THEN
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB( J)
            ENDIF
           ENDDO
          ENDDO
C
      ELSE IF (KIND.EQ.2) THEN
C                             i.e. (A,B) = (q,Qbar):
C
          DO I = 1,NFL
           DO J = 1,NFL
            IF (J.NE.I) THEN
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB(-J)
            ENDIF
           ENDDO
          ENDDO
C
          DO I = -NFL,-1
           DO J = -NFL,-1
            IF (J.NE.I) THEN
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB(-J)
            ENDIF
           ENDDO
          ENDDO
C
      ELSE IF (KIND.EQ.3) THEN
C                                i.e. (A,B) = (q,q):
C
          DO I = 1,NFL
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB( I)
          ENDDO
C
          DO I = -NFL,-1
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB( I)
          ENDDO
C
      ELSE IF (KIND.EQ.4) THEN
C                              i.e. (A,B) = (q,qbar):
C
          DO I = 1,NFL
              PRTNS =   PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB(-I)
          ENDDO
C
          DO I = -NFL,-1
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB(-I)
          ENDDO
C
      ELSE IF (KIND.EQ.5) THEN
C                             i.e. (A,B) = (q,g):
C
          DO I = 1,NFL
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB( 0)
          ENDDO
C
          DO I = -NFL,-1
              PRTNS =  PRTNS
     >                   + PRTNA( I)
     >                    *PRTNB( 0)
          ENDDO
C
      ELSE IF (KIND.EQ.6) THEN
C                                i.e. (A,B) = (g,q):
C
          DO I = 1,NFL
              PRTNS =  PRTNS
     >                   + PRTNA( 0)
     >                    *PRTNB( I)
          ENDDO
C
          DO I = -NFL,-1
              PRTNS =  PRTNS
     >                   + PRTNA( 0)
     >                    *PRTNB( I)
          ENDDO
C
      ELSE IF (KIND.EQ.7) THEN
C                             i.e. (A,B) = (g,g):
C
           PRTNS =      PRTNA( 0)
     >                    *PRTNB( 0)
C
      ENDIF
C
      LMNSTY(KIND) =  PRTNS
C
C  End loop over KIND
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C             END OF SUBROUTINES FOR 4 DIMENSIONAL INTEGRATIONS
C***********************************************************************
C
C***********************************************************************
C       MISCELLANEOUS LIBRARY ROUTINES FOR JET CALCULATION
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION ALPHAS(Q)
C 
      DOUBLE PRECISION Q
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
C
      DOUBLE PRECISION PI
      PARAMETER ( PI = 3.141592654 D0)
      DOUBLE PRECISION FACTOR,LOG
C
      FACTOR = 33 - 2.0D0*NFLAVOR
      LOG = DLOG( Q**2 / LAMBDA**2)
      ALPHAS = 12.0D0 * PI /FACTOR /LOG
     >    - 72.0D0 * PI * (153.0D0 - 19.0D0*NFLAVOR) * DLOG(LOG)
     >          /FACTOR**3 /LOG**2
C
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION WNUM(NPARTON)
C 
      INTEGER NPARTON
C
C   This function gives the spin and color weighting for
C   parton NPARTON: 6 for quarks and antiquarks, 16 for gluons
C
      IF (NPARTON.EQ.0) THEN
        WNUM = 16.0D0
      ELSE
        WNUM = 6.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION COLOR(NPARTON)
C 
      INTEGER NPARTON
C
C   This function gives the color charge for
C   parton NPARTON: 4/3 for quarks and antiquarks, 3 for gluons
C
      DOUBLE PRECISION CA,CF
      PARAMETER( CA = 3.0D0 )
      PARAMETER( CF = 4.0D0/3.0D0)
C
      IF (NPARTON.EQ.0) THEN
        COLOR = CA
      ELSE
        COLOR = CF
      ENDIF
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION THETA(BOOL)
C 
      LOGICAL BOOL
      IF (BOOL) THEN
          THETA = 1.0D0
      ELSE
          THETA = 0.0D0
      ENDIF
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION FEXP(Y)
C 
      DOUBLE PRECISION Y
C  THE PURPOSE OF THIS FUNCTION IS JUST TO PREVENT A PROGRAM CRASH
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
      IF (Y.GT.35.0D0) THEN
          FEXP=DEXP(35.0D0)
          WRITE(NWRT,101)Y
      ELSE IF (Y.LT.-35.0D0) THEN
          FEXP=DEXP(-35.0D0)
          WRITE(NWRT,101)Y
      ELSE
          FEXP=DEXP(Y)
      ENDIF
      RETURN
101   FORMAT(' Y = ',G10.2,'IN FEXP(Y); NO HARM IF THIS IS RARE.')
      END
C
C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION UNSMEAR(M,J,DELTA)
      INTEGER M,J
      DOUBLE PRECISION DELTA
C
C  Calculates the 'unsmearing coefficients'
C
C  UNSMEAR(M,J,DELTA) = (-1 DELTA^2 /2)^(J-M)  (2J)! / (2M)! /(J-M)!
C
C  where DELTA is the smearing length.
C
      DOUBLE PRECISION FACTORIAL
C
      UNSMEAR = ( - DELTA**2 / 2.0D0 )**(J-M) 
     >         * FACTORIAL(2*J) / FACTORIAL(2*M) / FACTORIAL(J-M)
C
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION DELTA(M,J)
      INTEGER M,J
C
C  Kroneker delta
C
      IF(M.EQ.J) THEN
       DELTA = 1.0D0
      ELSE
       DELTA = 0.0D0
      ENDIF
C
      RETURN
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION FACTORIAL(N)
      INTEGER N
C
      INTEGER J
C
      FACTORIAL = 1.0D0
      DO J = 1,N
        FACTORIAL = FACTORIAL * J
      ENDDO
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE RANDOMINIT(IRAN)
      INTEGER IRAN
C
C  Initialize the shift-register random number generator from a random
C  seed  IRAN, 0 <= IRAN < 259200, with the help of a portable "quick 
C  and dirty" generator
C
      INTEGER LBIT
      DOUBLE PRECISION FAC
      PARAMETER(LBIT=31,FAC= 2.0D0**(-LBIT))
C  Here LBIT is the number of bits used in the shift register generator
      INTEGER J,IR(250)
      DOUBLE PRECISION RR(250)
      COMMON/RANDO/ RR,J,IR
      SAVE /RANDO/
      INTEGER I,LB,JRAN,IFAC1,ISK,IDI,I1S
C  Configuration of the shift register generator
      INTEGER IM,IA,IC
      DATA IM,IA,IC /259200,7141,54773/
C
C
      IF(IRAN.LT.0.OR.IRAN.GE.IM)STOP
     > 'RINI: IRAN OUT OF RANGE'
      JRAN=IRAN
C  Warm up the auxiliary generator a little
      DO I=1,10
       JRAN=MOD(JRAN*IA+IC,IM)
      ENDDO
      IFAC1=((2**(LBIT-1)-1)*2+1)/IM
      DO I=2,250
       JRAN=MOD(JRAN*IA+IC,IM)
       IR(I)=IFAC1*JRAN
      ENDDO
C  Guarantee LBIT linearly independent (over the field of 2 el.)
C  elements in IR(I):
      IR(1)=1
      I=1
      ISK=250/LBIT
      IDI=1
      I1S=1
      DO LB=1,LBIT-1
       I=I+ISK
       IDI=IDI*2
       I1S=I1S+IDI
       IR(I)=IR(I).OR.IDI
       IR(I)=IR(I).AND.I1S
      ENDDO
      CALL NEWRAN
      END
C
C***********************************************************************
C
      SUBROUTINE NEWRAN
C
C  Fills IR(I),RR(I) with 250 new random numbers, resets J=0.
C  Increment J before use!
C
      INTEGER LBIT
      DOUBLE PRECISION FAC
      PARAMETER(LBIT=31,FAC= 2.0D0**(-LBIT))
      INTEGER J,IR(250)
      DOUBLE PRECISION RR(250)
      COMMON/RANDO/ RR,J,IR
      SAVE /RANDO/
C
      INTEGER N
C
      DO N=1,103
       IR(N)=IR(N+147).XOR.IR(N)
       RR(N)=FAC*(DBLE(IR(N)) + 0.5D0)
      ENDDO
C
      DO N=104,250
       IR(N)=IR(N-103).XOR.IR(N)
       RR(N)=FAC*(DBLE(IR(N)) + 0.5D0)
      ENDDO
C
      J=0
      END
C
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION RANDOM(DUMMY)
      REAL DUMMY
C
C    Random number between 2**(-32) and 1 - 2**(-32):
C
      INTEGER J,IR(250)
      DOUBLE PRECISION RR(250)
      COMMON/RANDO/ RR,J,IR
      SAVE /RANDO/
C
      IF(J.GE.250)CALL NEWRAN
      J=J+1
      RANDOM = RR(J)
      END
C
C***********************************************************************
C             END OF LIBRARY ROUTINES FOR JET CALCULATION
C***********************************************************************
