C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        JET version 4.0
C
C            Version to calculate the two jet cross section
C                       d sigma/ d XA d XB d YSTAR
C
C                   S.D. Ellis, Z. Kunszt, D.E. Soper 
C
C                           14 June 1996
C
C                    Version 3.0,  23 February 1991
C                    Version 3.01, 14 March 1991
C                    Version 3.1,  27 February 1992
C                    Version 4.0,  27 September 1994
C                     6  October 1994: Set PSCALE = 0.05D0 * RTS
C                     10 October 1994: Slightly reduce max mu to avoid
C                                       out of range in PARTON.
C                     29 October 1994: Speed up BINIT
C                        14 June 1996: Incorporates generic parton.f 
C                                      version 1.1, made available
C                                      on World Wide Web
C
C This program needs jetsubs4_0.f. The program makes a table of
C K factors, which can then by read by readjet4_0.f.
C
C***********************************************************************
C************************BEGIN MAIN PROGRAM*****************************
C
      PROGRAM JET
C 
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
      DOUBLE PRECISION AUV,ACO
      COMMON /MUCOEFS/ AUV,ACO
C
      CHARACTER DATE*(28)
      REAL CPUTIME,DTIME,TIMEARRAY(2)
C
C  Switches for turning on parts of the integrand
C     .TRUE.  turns the corresponding part on
C     .FALSE. sets it to zero
C  In common /SW2TO3/,  the switch S2TO2 turns ALL of the 2 to 2 parts
C  off (independent of the 2 to 2 switches) if it is set to .FALSE.
C
      LOGICAL STYPEA,STYPEB,STYPEC,STYPED
      COMMON /SWTYPE/ STYPEA,STYPEB,STYPEC,STYPED
C
      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
C
      LOGICAL S2TO2
      LOGICAL SAFINITE,SBFINITE,S1FINITE,S2FINITE
      COMMON /SW2TO3/S2TO2,SAFINITE,SBFINITE,S1FINITE,S2FINITE
C
C  Switch for p-pbar collisions (if true) or p-p collisions (if false)
C
      LOGICAL PPBAR
      COMMON /BEAMS/ PPBAR
C
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
C   Parton distribution function
C
      DOUBLE PRECISION PARTON,PRTN
      CHARACTER*20 PARTONFILE
      COMMON /PARTONCHOICE/ PARTONFILE
C
C   RENO variables
C
      INTEGER SEED,NRENO
C
C   Name for output file and variable name for the output
C
      CHARACTER*10 FILENAME,VARNAME
C
C   Code to indentify version of this program in output file
C
      INTEGER VERSION
C
      DOUBLE PRECISION YSCALE,PSCALE
      COMMON /RENOSCALES/ YSCALE,PSCALE
C
      INTEGER NL,NYSTAR
      DOUBLE PRECISION LGRID(1:16),YSTARGRID(1:16)
      COMMON /GRIDS/  LGRID,YSTARGRID,NL,NYSTAR
C
      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
      DOUBLE PRECISION YSTARMAX,DELTAYSTAR
      DOUBLE PRECISION X1,XN,DELTAL
C
C   Results
C
      DOUBLE PRECISION LA,LB,YSTAR,XA,XB
      DOUBLE PRECISION KFACTORM1,KERR
C
C   Smearing and unsmearing parameters
C
      DOUBLE PRECISION  LSMEAR,YSMEAR
      COMMON /SMEAR/ LSMEAR,YSMEAR
C
C   Number of terms in jet defining function
C
      INTEGER NJETDEF
      COMMON /DEFJET/ NJETDEF
C
C  Bounds for jet variables
C
      DOUBLE PRECISION ETMIN,ETMAX,YBOUND
      COMMON /JETBOUNDS/ ETMIN,ETMAX,YBOUND
C
C   Integers for DO loops
C
      INTEGER I,IYSTAR,ILA,ILB
C
      DOUBLE PRECISION PI
      PARAMETER ( PI = 3.141592654 D0)
      DOUBLE PRECISION NBGEV2
      PARAMETER (NBGEV2 = 0.38938573 D6)
C
C------------INITIALIZATION-------------------
C
C---Timing function for Sun workstations; remove for others
C
      CPUTIME = DTIME(TIMEARRAY)
C---
C
C--------Adjustable parameters----------------
C
C
C  Switch for p-pbar collisions (if true) or p-p collisions (if false)
C
      PPBAR = .TRUE.
C
C  Lambda_MSbar^(NFL): set in initialization of PARTON
C
C  Input and output units.
C
      NIN  = 5
      NOUT = 6
      NWRT = 6
C
C  Number of terms in S-function defining the observable (normally 3).
C
      NJETDEF = 3
C
C  Jet physics is normally in the range of five flavors, so we
C  set the number of flavors to 5. (This could be changed for special
C  purposes, but the program works with a fixed NFL. Scheme switching
C  for as mu becomes greater than a heavy quark mass is not 
C  implemented.)
C
      NFLAVOR = 5.0D0
      NFL = 5
C
C  Switches (should all be .TRUE.)
C
      STYPEA    = .TRUE.
      STYPEB    = .TRUE.
      STYPEC    = .TRUE.
      STYPED    = .TRUE.
C
C  Physics parameters:
C
      WRITE(NOUT,*)' Give RTS,R'
      READ(NIN,*) RTS,R
C
C  Grid points in Log( x/(1-x) )
C
      WRITE(NOUT,*)' Give X grid: Xmin, Xmax, Number of points'
      READ(NIN,*) X1,XN,NL
      IF (NL.GT.16) THEN
        WRITE(NOUT,*) 'Need number of points .LE. 16 please.'
        STOP
      ENDIF
C
C  Grid points in Y
C
      WRITE(NOUT,*)' Give YSTAR grid: YSTARMAX, Number of points'
      READ(NIN,*) YSTARMAX,NYSTAR
      IF (NYSTAR.GT.16) THEN
        WRITE(NOUT,*) 'Need number of points .LE. 16 please.'
        STOP
      ENDIF
C
C  Array limits
C
      SIZE1 = NL
      SIZE2 = NL
      SIZE3 = NYSTAR
C
      WRITE(NOUT,*)' Give smearing distances in ln(x/(1-x)) and ystar.'
      READ(NIN,*) LSMEAR,YSMEAR
C
C  Parton set desired:
C
      WRITE(NOUT,*) ' Give file name of parton distribution data.'
      WRITE(NOUT,*) ' (Limit is 20 characters.)'
      WRITE(NOUT,*) ' (Typical example: cteq3m.ptn)'
      READ(*,*) PARTONFILE
C
C  Ratio of renormalization scale MUUV and factorization
C  scale MUCO to PT:
C
      WRITE(NOUT,*) ' Give scale choices: A_uv, A_co'
      READ(NIN,*) AUV,ACO
C
C  Integration parameters:
C
      WRITE(NOUT,*)' Give SEED ( 0 <= SEED < 259200)'
      READ(NIN,*) SEED
C
      CALL RANDOMINIT(SEED)
C
      WRITE(NOUT,*)' Give number of thousands of RENO points'
      READ(NIN,*) NRENO
C
      WRITE(NOUT,*)' Give name for output file'
      READ(NIN,*) FILENAME
C
C  END OF INPUT
C
      OPEN(100,FILE=FILENAME,STATUS='NEW',ERR=99)
C
C  Make ystar-grid.
C
      IF(NYSTAR.GT.1)THEN
        DELTAYSTAR = YSTARMAX/(NYSTAR-1)
        DO  I = 1,NYSTAR
         YSTARGRID(I) =  (I - 1) * DELTAYSTAR
        ENDDO
      ELSE
        YSTARGRID(1) = 0.0D0
      ENDIF      
C
C   Make L-grid.
C
      IF(NL.GT.1)THEN
        DELTAL = (  DLOG(XN/(1.0D0 - XN)) 
     >            - DLOG(X1/(1.0D0 - X1)) )/(NL-1)
      ELSE
        DELTAL = 1.0D0
      ENDIF
      DO  I = 1,NL
         LGRID(I) = DLOG(X1/(1.0D0 - X1)) + (I-1)*DELTAL
      ENDDO
C
C  Scales to tell RENO where to put rapidities and transverse
C  momenta for partons 1 and 2.
C
      YSCALE = 2.0D0
      PSCALE = 0.05D0 * RTS
C
C   WRITE THE PARAMETERS
C
      CALL DAYTIME(DATE)
      WRITE(NOUT,101)DATE
101   FORMAT(//
     >    ,' Program JET version 4.0 ',A,/,50('-'))
C
      WRITE(NOUT,*)'Two jet cross section d sigma / d XA d XB d YSTAR'
      WRITE(NOUT,*)'The jet variables are (sums over partons in jets)'
      WRITE(NOUT,*)'   XA = Sum_i (Pi/RTS) EXP( Yi),'
      WRITE(NOUT,*)'   XB = Sum_i (Pi/RTS) EXP(-Yi),'
      WRITE(NOUT,*)'   YSTAR = |YJ1-YJ2|/2.'
      WRITE(NOUT,87)
87    FORMAT(50('-'))
C
      CALL PRINTVERSION
C
      WRITE(100,102)DATE
102   FORMAT(' File written by JET version 4.0 ',A,/,50('-'))
      WRITE(100,*)' d sigma/ d XA d XB dYstar'
      IF (PPBAR) THEN
        WRITE(100,103)
103     FORMAT(' Proton-antiproton collisions')
      ELSE
        WRITE(100,104)
104     FORMAT(' Proton-proton collisions    ')
      ENDIF
      WRITE(100,*)' '
      VERSION = 940923
      WRITE(100,202)'   VERSION',VERSION
      VARNAME = '       RTS'
      WRITE(100,201)VARNAME,RTS
      VARNAME = '         R'
      WRITE(100,201)VARNAME,R
      VARNAME = '        X1'
      WRITE(100,201)VARNAME,X1
      VARNAME = '        XN'
      WRITE(100,201)VARNAME,XN
      VARNAME = '        NL'
      WRITE(100,202)VARNAME,NL
      VARNAME = '  YSTARMAX'
      WRITE(100,201)VARNAME,YSTARMAX
      VARNAME = '    NYSTAR'
      WRITE(100,202)VARNAME,NYSTAR
      VARNAME = '    LSMEAR'
      WRITE(100,201)VARNAME,LSMEAR
      VARNAME = '    YSMEAR'
      WRITE(100,201)VARNAME,YSMEAR
      VARNAME = '   PARTONS'
      WRITE(100,204)VARNAME,PARTONFILE
      VARNAME = '       AUV'
      WRITE(100,201)VARNAME,AUV
      VARNAME = '       ACO'
      WRITE(100,201)VARNAME,ACO
      VARNAME = '      SEED'
      WRITE(100,202)VARNAME,SEED
      VARNAME = '     NRENO'
      WRITE(100,202)VARNAME,NRENO
201   FORMAT(A10,' = ',G20.10)
202   FORMAT(A10,' = ',I10)
203   FORMAT(A10,' = ',A10)
204   FORMAT(A10,' = ',A20)
      WRITE(100,*)''
      WRITE(100,*)'Below are: YSTAR, XA, XB, KFACTORM1, KERR'
      WRITE(100,*)''
C
        S2TO2     = .TRUE.
         SBORN    = .FALSE.
         SVIRTUAL = .TRUE.
         SACOLL   = .TRUE.
         SASOFT   = .TRUE.
         SBCOLL   = .TRUE.
         SBSOFT   = .TRUE.
         S1COLL   = .TRUE.
         S1SOFT   = .TRUE.
         S2COLL   = .TRUE.
         S2SOFT   = .TRUE.
        SAFINITE  = .TRUE.
        SBFINITE  = .TRUE.
        S1FINITE  = .TRUE.
        S2FINITE  = .TRUE.
        WRITE(NOUT,*)' '
        WRITE(NOUT,*)'Calculate (K-factor - 1.0)'
        WRITE(NOUT,*)' '
C
      IF (PPBAR) THEN
        WRITE(NOUT,*) 'Proton-antiproton collisions'
      ELSE
        WRITE(NOUT,*) 'Proton-proton collisions'
      ENDIF
        WRITE(NOUT,*)' '
C
C  This call of PARTON initializes it and creates some output.
C
      PRTN =  PARTON(0,0.01D0,100.0D0)
      WRITE(NOUT,617) PRTN
617     FORMAT(' e.g. PARTON(0,0.01D0,100.0D0) =',G15.4,/)
C
      WRITE(NOUT,613)RTS,NFL,LAMBDA,R
613     FORMAT('           collision energy (GeV):', F10.1,/,
     >         '                          nflavor:', I10 ,/,
     >         '                     lambda (GeV):', F10.3 ,/,
     >         '                      cone radius:', F10.3,/)
C
      WRITE(NOUT,654)AUV,ACO    
654     FORMAT(' Choice of scales for mu:',/,
     > '   mu_A  = A_co SQRT(XA XB S)/[4 COSH((1-XA) YSTAR)]',/,
     > '   mu_B  = A_co SQRT(XA XB S)/[4 COSH((1-XB) YSTAR)]',/,
     > '   mu_uv = A_uv SQRT(XA XB S)/[4 COSH(0.7 YSTAR)], ',/,
     > '   where',/,
     > '                             A_uv:',F10.3,/,
     > '                             A_co:',F10.3,/)
C
      WRITE(NOUT,612) NRENO,SEED
612     FORMAT(' RENO will use',I6,' thousand points',
     >         ' with seed',I7,/)
      WRITE(NOUT,432)LSMEAR,YSMEAR
432     FORMAT(' Smearing information:',/,
     >         ' Smearing distance in ln(X/(1-X)):',F10.3,/,
     >         '       Smearing distance in YSTAR:',F10.3,/)
C
C
C---------------BEGIN CALCULATION----------------------
C
C
C   Calculate  cross section
C
      CALL RENO(NRENO)
C
C   Results are now in XSECTARRAY,ERRARRAY.  These results represent
C   the cross section (not including the Born cross section) over
C   BORN(Y1,P2,Y2), smeared over YSTAR and LA = Log ( XA/(1-XA) ) and
C   LB = Log ( XB/(1-XB) ).
C  
C---------  Loop over grid points, storing results for each.
C
      DO IYSTAR = 1,NYSTAR
      DO ILA = 1,NL
      DO ILB = 1,NL
C
        KFACTORM1 =  XSECTARRAY(ILA,ILB,IYSTAR)
C
C   Calculate the standard deviation.
C
      KERR = ERRARRAY(ILA,ILB,IYSTAR)
      KERR = DSQRT((KERR - XSECTARRAY(ILA,ILB,IYSTAR)**2)/(NRENO - 1))
C
      LA = LGRID(ILA)
      LB = LGRID(ILB)
      YSTAR = YSTARGRID(IYSTAR)
      XA = DEXP(LA) /( 1.0D0 + DEXP(LA) )
      XB = DEXP(LB) /( 1.0D0 + DEXP(LB) )
C
C   Write results
C
      WRITE(100,15)YSTAR,XA,XB,KFACTORM1,KERR 
15    FORMAT(5G20.10)
C
      ENDDO
      ENDDO
      ENDDO
C
      CLOSE(100)
C
C----------------------------------------------
C 
      CALL DAYTIME(DATE)
      WRITE(NOUT,567)DATE
567   FORMAT(/
     >    ,' DONE ',A,/)
C
C---Timing function for Sun workstations; remove for others
C
      CPUTIME = DTIME(TIMEARRAY)
      WRITE(NOUT,*)'  CPUTIME = ',CPUTIME
C---
C
      STOP
C
99    WRITE(NOUT,*)'Error opening data file for output'
      STOP
C
      END
C
C***********************************************************************
C***********************************************************************
C
C Sun version for DAYTIME
C
      SUBROUTINE DAYTIME(DATE)
      CHARACTER*28 DATE
C 
      CHARACTER*30 SUNDATE
      CALL FDATE(SUNDATE)
      DATE = SUNDATE(1:28)
      RETURN
      END
C
C***********************************************************************
C********************** Parton Distributions ***************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION PARTON(NPARTON,XIN,SCALE)
      INTEGER NPARTON
      DOUBLE PRECISION XIN,SCALE
C
C Fortran function to interpolate parton distributions from a data
C file. Uses cubic spline interpolation in the variables
C LX = ln(X/(1-X)) and LMU = ln(MU).
C
C Version 1.1
C 9 May 1996
C D. E. Soper and P. Anandam
C Version 1.0
C 25 August 1994
C D. E. Soper and S. D. Ellis
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
      CHARACTER*20 PARTONFILE
      COMMON /PARTONCHOICE/ PARTONFILE
C
      CHARACTER*100 HEADER
      CHARACTER*100 VARNAME
      CHARACTER*100 DUMMY
      INTEGER VERSION
      INTEGER LAMBDAFLAVORS
C
      DOUBLE PRECISION X,XMIN,XMAX
      DOUBLE PRECISION LX,LXMIN,LXMAX,DELTALX
      DOUBLE PRECISION MU,MUMIN,MUMAX
      DOUBLE PRECISION LMUMIN,LMU,LMUMAX,DELTALMU
      INTEGER N,NA,NB,IA,IB
      DOUBLE PRECISION NORMALIZEDLX,NORMALIZEDLMU
C
      INTEGER NAMAX,NBMAX
      PARAMETER(NAMAX=64)
      PARAMETER(NBMAX=32)
      DOUBLE PRECISION  H(-2:6,NAMAX,NBMAX), HA(-2:6,NAMAX,NBMAX)
      DOUBLE PRECISION HB(-2:6,NAMAX,NBMAX),HAB(-2:6,NAMAX,NBMAX)
      DOUBLE PRECISION HBL,HBR,HABL,HABR
C
      DOUBLE PRECISION A,AA,ASQ,AASQ,A0,A0A,A1,A1A
      DOUBLE PRECISION B,BB,BSQ,BBSQ,B0,B0B,B1,B1B
      LOGICAL EXTRAPOLATE
      DOUBLE PRECISION NET
C
      INTEGER INIT
C      
      DATA INIT/0/
C
C Note that we save the variables!
C
      SAVE
C
C ------------------- Initialization --------------------------
C
      IF(INIT.NE.0) GOTO 10
      INIT=1
C
      OPEN(28, FILE=PARTONFILE, STATUS='OLD', ERR=899)
      READ(28,101) HEADER
      HEADER = HEADER(INDEX(HEADER,'!')+1:)
      VARNAME='!'
      DO WHILE (INDEX(VARNAME,'!').NE.0)
         READ (28,*) VARNAME
      ENDDO
      IF ((INDEX(VARNAME,'PARTON').EQ.0).AND.
     >    (INDEX(VARNAME,'Parton').EQ.0).AND.
     >    (INDEX(VARNAME,'parton').EQ.0)) GOTO 199
      READ(28,*) VARNAME,DUMMY,VERSION
      IF ((INDEX(VARNAME,'VERSION').EQ.0).AND.
     >    (INDEX(VARNAME,'Version').EQ.0).AND.
     >    (INDEX(VARNAME,'version').EQ.0)) GOTO 199
      IF (VERSION.NE.110496) THEN
        WRITE(NOUT,*) 'Wanted parton file with version number 110496'
        WRITE(NOUT,*) 'Found version number',VERSION
        WRITE(NOUT,*) 'Given possible data format problems, I stop.'
        STOP
      ENDIF
      READ(28,*) VARNAME,DUMMY,LAMBDA
      IF ((INDEX(VARNAME,'LAMBDA').EQ.0).AND.
     >    (INDEX(VARNAME,'Lambda').EQ.0).AND.
     >    (INDEX(VARNAME,'lambda').EQ.0)) GOTO 199
      READ(28,*) VARNAME,DUMMY,LAMBDAFLAVORS
      IF ((INDEX(VARNAME,'NFL_LAMBDA').EQ.0).AND.
     >    (INDEX(VARNAME,'NFL_Lambda').EQ.0).AND.
     >    (INDEX(VARNAME,'nfl_lambda').EQ.0)) GOTO 199
      IF(LAMBDAFLAVORS.NE.NFL) THEN
        WRITE(NOUT,*)'Oops, NFL_Lambda is not NFL'
        STOP
      ENDIF
      READ(28,*) VARNAME,DUMMY,XMIN
      IF ((INDEX(VARNAME,'XMIN').EQ.0).AND.
     >    (INDEX(VARNAME,'Xmin').EQ.0).AND.
     >    (INDEX(VARNAME,'xmin').EQ.0)) GOTO 199
      READ(28,*) VARNAME,DUMMY,XMAX
      IF ((INDEX(VARNAME,'XMAX').EQ.0).AND.
     >    (INDEX(VARNAME,'Xmax').EQ.0).AND.
     >    (INDEX(VARNAME,'xmax').EQ.0)) GOTO 199
      READ(28,*) VARNAME,DUMMY,NA
      IF ((INDEX(VARNAME,'N_XPOINTS').EQ.0).AND.
     >    (INDEX(VARNAME,'N_XPoints').EQ.0).AND.
     >    (INDEX(VARNAME,'n_xpoints').EQ.0)) GOTO 199
      IF(NA.GT.NAMAX) THEN
        WRITE(NOUT,*)' N_XPOINTS too big for declared array size'
        STOP
      ENDIF
      READ(28,*) VARNAME,DUMMY,MUMIN
      IF ((INDEX(VARNAME,'MUMIN').EQ.0).AND.
     >    (INDEX(VARNAME,'MuMin').EQ.0).AND.
     >    (INDEX(VARNAME,'mumin').EQ.0)) GOTO 199
      READ(28,*) VARNAME,DUMMY,MUMAX
      IF ((INDEX(VARNAME,'MUMAX').EQ.0).AND.
     >    (INDEX(VARNAME,'MuMax').EQ.0).AND.
     >    (INDEX(VARNAME,'mumax').EQ.0)) GOTO 199
      READ(28,*) VARNAME,DUMMY,NB
      IF ((INDEX(VARNAME,'N_MUPOINTS').EQ.0).AND.
     >    (INDEX(VARNAME,'N_MuPoints').EQ.0).AND.
     >    (INDEX(VARNAME,'n_mupoints').EQ.0)) GOTO 199
      IF(NB.GT.NBMAX) THEN
        WRITE(NOUT,*)' N_MUPOINTS too big for declared array size'
        STOP
      ENDIF
      DO 20 IA=1,NA
      DO 20 IB=1,NB
      READ(28,*)H(-2,IA,IB),H(-1,IA,IB),H(0,IA,IB),
     >            H(1,IA,IB), H(2,IA,IB), H(3,IA,IB),
     >            H(4,IA,IB), H(5,IA,IB), H(6,IA,IB)
  20  CONTINUE
      CLOSE(28)
C
101   FORMAT(A80)
C
      WRITE(NOUT,601) PARTONFILE
601   FORMAT(' Read file ',A20,' with header:')
      WRITE(NOUT,602) HEADER
602   FORMAT(A80)
      WRITE(NOUT,603)LAMBDA,NFL
603   FORMAT(' These partons have LAMBDA =',F10.3,
     >       ' for NFL =',I3)    
C
C Calculate lattice info.
C
      LXMIN  = DLOG(XMIN/(1.0D0 - XMIN))
      LXMAX  = DLOG(XMAX/(1.0D0 - XMAX))
      DELTALX = (LXMAX - LXMIN)/(NA-1)
      LMUMIN = DLOG(MUMIN)
      LMUMAX = DLOG(MUMAX)
      DELTALMU = (LMUMAX - LMUMIN)/(NB-1)
C
C Initialize the derivatives
C Define the derivatives at the edges of the lattice too.
C
      DO N=-2,6
C
C Derivatives with respect to A.
C
      DO IA = 2,NA-1
      DO IB = 1,NB
        HA(N,IA,IB) = 0.5D0*(H(N,IA+1,IB ) - H(N,IA-1,IB ))
      ENDDO
      ENDDO
      DO IB = 1,NB
        HA(N,1,IB)  = H(N,2,IB) - H(N,1,IB)
        HA(N,NA,IB) = H(N,NA,IB) - H(N,NA-1,IB)
      ENDDO
C
C Derivatives with respect to B. We choose the smaller in absolute
C value of the derivatives to the left and to the right because there
C is a very big derivative locally when a heavy quark "turns on" and
C we want to avoid wiggles in the interpolation.
C
      DO IA = 1,NA
      DO IB = 2,NB-1
        HBL = H(N,IA,IB) - H(N,IA,IB-1)
        HBR = H(N,IA,IB+1) - H(N,IA,IB)
        IF (DABS(HBL).LT.DABS(HBR)) THEN
          HB(N,IA,IB) = HBL
        ELSE
          HB(N,IA,IB) = HBR
        ENDIF
      ENDDO
      ENDDO
      DO IA = 1,NA
        HB(N,IA,1)  = H(N,IA,2) - H(N,IA,1)
        HB(N,IA,NB) = H(N,IA,NB) - H(N,IA,NB-1)
      ENDDO
C
C Second derivatives with respect to A and B. Same logic as before
C for dH/ dB.
C
      DO IA = 2,NA-1
      DO IB = 2,NB-1
        HABL = 0.5D0*(  H(N,IA+1,IB) - H(N,IA-1,IB)
     >                - H(N,IA+1,IB-1) + H(N,IA-1,IB-1) )
        HABR = 0.5D0*(  H(N,IA+1,IB+1) - H(N,IA-1,IB+1)
     >                - H(N,IA+1,IB) + H(N,IA-1,IB) )
        IF (DABS(HABL).LT.DABS(HABR)) THEN
          HAB(N,IA,IB) = HABL
        ELSE
          HAB(N,IA,IB) = HABR
        ENDIF
      ENDDO
      ENDDO
      DO IA = 2,NA-1
        HAB(N,IA,1)  = 0.5D0*(  H(N,IA+1,2) - H(N,IA-1,2)
     >                        - H(N,IA+1,1) + H(N,IA-1,1) )
        HAB(N,IA,NB) = 0.5D0*(  H(N,IA+1,NB)  -  H(N,IA-1,NB)
     >                        - H(N,IA+1,NB-1) + H(N,IA-1,NB-1) )
      ENDDO
      DO IB = 2,NB-1
        HABL =  H(N,2,IB)   - H(N,1,IB)
     >        - H(N,2,IB-1) + H(N,1,IB-1)
        HABR =  H(N,2,IB+1) - H(N,1,IB+1)
     >        - H(N,2,IB)   + H(N,1,IB)
        IF (DABS(HABL).LT.DABS(HABR)) THEN
          HAB(N,1,IB) = HABL
        ELSE
          HAB(N,1,IB) = HABR
        ENDIF
        HABL =  H(N,NA,IB)   - H(N,NA-1,IB)
     >        - H(N,NA,IB-1) + H(N,NA-1,IB-1)
        HABR =  H(N,NA,IB+1) - H(N,NA-1,IB+1)
     >        - H(N,NA,IB)   + H(N,NA-1,IB)
        IF (DABS(HABL).LT.DABS(HABR)) THEN
          HAB(N,NA,IB) = HABL
        ELSE
          HAB(N,NA,IB) = HABR
        ENDIF
      ENDDO
      HAB(N,1,1)   =  H(N,2,2)     - H(N,1,2)
     >              - H(N,2,1)     + H(N,1,1)
      HAB(N,NA,1)  =  H(N,NA,2)    - H(N,NA-1,2)
     >              - H(N,NA,1)    + H(N,NA-1,1)
      HAB(N,1,NB)  =  H(N,2,NB)    - H(N,1,NB) 
     >              - H(N,2,NB-1)  + H(N,1,NB-1)
      HAB(N,NA,NB) =  H(N,NA,NB)   - H(N,NA-1,NB)
     >              - H(N,NA,NB-1) + H(N,NA-1,NB-1)
C
C End DO N=1,9
C
      ENDDO
C
C --------------- Main function routine -----------------
C
  10  CONTINUE
C
C -------------------------------------------------------
C   Translate parton numbers.  In data table we have
C   dbar ubar glue  u   d   s   c   b   t
C    -2   -1    0   1   2   3   4   5   6
C   with sbar = s, cbar = c, bbar = b, tbar = b.
C -------------------------------------------------------
C
      IF ((NPARTON.GE.-2).AND.(NPARTON.LE.6)) THEN
         N = NPARTON
      ELSE IF (NPARTON.GE.-6) THEN
         N = -NPARTON
      ELSE
        WRITE(NOUT,*)' NPARTON OUT OF RANGE',NPARTON
        STOP
      ENDIF
C
C We interpolate in LMU = LOG(MU). We define a normalized
C version of LMU such that the lattice points are at 1,2,...,NB.
C We define B and IB such that the normalized LMU is IB + B with
C 0 < B < 1.
C
C Check whether the previous scale 'MU' equals the new one. If
C it does, we can just use a lot of the old variables.
C
      IF (SCALE.NE.MU) THEN
C
      MU = SCALE
      IF ( ( MU.LE.MUMIN ).OR.(MU.GT.MUMAX )) THEN
          WRITE(NOUT,*)'PARTON called out of range, MU =',MU
          STOP
      ENDIF
      LMU = DLOG(MU)
      NORMALIZEDLMU = (LMU - LMUMIN)/DELTALMU + 1.0D0
      IB = INT(NORMALIZEDLMU)  
      IF((IB.LT.1).OR.(IB.GT.(NB-1))) THEN
         WRITE(NOUT,*)' IB OUT OF RANGE, IB =',IB
         STOP
      ENDIF
      B = NORMALIZEDLMU - DBLE(IB)
      BB = 1.0D0 - B
      BSQ = B**2
      BBSQ = BB**2
C
C Functions of B that have
C f(0) = 1, f'(0) = 0, f(1) = 0, f'(1) = 0;
C f(0) = 0, f'(0) = 1, f(1) = 0, f'(1) = 0;
C f(0) = 0, f'(0) = 0, f(1) = 1, f'(1) = 0;
C f(0) = 0, f'(0) = 0, f(1) = 0, f'(1) = 1.
C
      B0  =   (1.0D0 + 2.0D0*B) * BBSQ
      B0B =   B   * BBSQ
      B1  =   BSQ * (1.0D0 + 2.0D0*BB)
      B1B = - BSQ * BB
C
C End IF (SCALE.NE.MU) THEN
C
      ENDIF
C
C We interpolate in LX = LOG(X/(1-X)). We define a normalized
C version of LX such that the lattice points are at 1,2,...,NA.
C We define A and IA such that the normalized LX is IA + AB with
C 0 < A < 1.
C
C Check whether the previous momentum fraction 'X' equals the new one.
C If it does, we can just use a lot of the old variables.
C
      IF (XIN.NE.X) THEN
C
      X = XIN
      IF (X.LE.0.0D0) THEN
        WRITE(NOUT,*) 'PARTON called for negative x!'
        STOP
      ENDIF
      IF (X.GE.1.0D0) THEN
        PARTON = 0.0D0
        RETURN
      ENDIF
C
      LX = LOG(X/(1.0D0-X))
      NORMALIZEDLX = (LX - LXMIN)/DELTALX + 1.0D0
      IA = INT(NORMALIZEDLX)
C
C Special case: if X is outside the range of the table, then we
C extrapolate LOG(PARTON) as a linear function of LOG(X/(1-X)).
C
      IF(IA.LT.1) THEN
        EXTRAPOLATE = .TRUE.
        IA = 1
      ELSEIF(IA.GE.NA) THEN
        EXTRAPOLATE = .TRUE.
        IA = NA
      ELSE
        EXTRAPOLATE = .FALSE.
      ENDIF
C
      A = NORMALIZEDLX - DBLE(IA)
      AA = 1.0D0 - A
      ASQ = A**2
      AASQ = AA**2 
      A0  =   (1.0D0 + 2.0D0*A) * AASQ
      A0A =   A   * AASQ
      A1  =   ASQ * (1.0D0 + 2.0D0*AA)
      A1A = - ASQ * AA
C
C End IF (XIN.NE.X) THEN
C
      ENDIF
C
C Use the magic functions to construct the interpolation; but
C first check if we were supposed to extrapolate instead of
C interpolate.
C
      IF (EXTRAPOLATE) THEN
C
        NET =       H(N,IA,  IB  ) * B0
        NET = NET + H(N,IA,  IB+1) * B1
        NET = NET + HA(N,IA  ,IB  )  * A * B0
        NET = NET + HA(N,IA  ,IB+1)  * A * B1
        NET = NET + HB(N,IA  ,IB  )  * B0B
        NET = NET + HB(N,IA  ,IB+1)  * B1B
        NET = NET + HAB(N,IA  ,IB  )  * A * B0B
        NET = NET + HAB(N,IA  ,IB+1)  * A * B1B
C
      ELSE
C     
        NET =       H(N,IA,  IB  ) * A0 * B0
        NET = NET + H(N,IA+1,IB  ) * A1 * B0
        NET = NET + H(N,IA,  IB+1) * A0 * B1
        NET = NET + H(N,IA+1,IB+1) * A1 * B1
C
        NET = NET + HA(N,IA  ,IB  )  * A0A * B0
        NET = NET + HA(N,IA+1,IB  )  * A1A * B0
        NET = NET + HA(N,IA  ,IB+1)  * A0A * B1
        NET = NET + HA(N,IA+1,IB+1)  * A1A * B1
C
        NET = NET + HB(N,IA  ,IB  )  * A0 * B0B
        NET = NET + HB(N,IA+1,IB  )  * A1 * B0B
        NET = NET + HB(N,IA  ,IB+1)  * A0 * B1B
        NET = NET + HB(N,IA+1,IB+1)  * A1 * B1B
C
        NET = NET + HAB(N,IA  ,IB  )  * A0A * B0B
        NET = NET + HAB(N,IA+1,IB  )  * A1A * B0B
        NET = NET + HAB(N,IA  ,IB+1)  * A0A * B1B
        NET = NET + HAB(N,IA+1,IB+1)  * A1A * B1B
C
      ENDIF
C
      PARTON = DEXP(NET) - 1.0D-16
      IF(PARTON.LT.(1.0D-16)) PARTON = 0.0D0 
C
      RETURN
C
C-------
C
899    WRITE(NOUT,*) 'Error opening parton data file!'
       STOP
199    WRITE(NOUT,*) 'Wrong format of parton data file!',VARNAME
       STOP
C
      END
C
C***********************************************************************
C***********************************************************************
C
      LOGICAL FUNCTION INBOUNDS(LA,LB,YSTAR)
C 
      DOUBLE PRECISION LA,LB,YSTAR
C
      INBOUNDS = .TRUE.
C
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE SETMU(LA,LB,YSTAR,MUUV,MUA,MUB)
C 
C  Input variables:
      DOUBLE PRECISION LA,LB,YSTAR
C  Output variables:
      DOUBLE PRECISION MUUV,MUA,MUB
C
C  The input arguments are JETVAR1,JETVAR2,JETVAR3, which are the jet
C  parameters produced by JETDEF and JETDEFBORN.
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION AUV,ACO
      COMMON /MUCOEFS/ AUV,ACO
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      DOUBLE PRECISION FEXP
      DOUBLE PRECISION ELA,ELB,XA,XB,RTSHAT,UVSCALE,ASCALE,BSCALE
C
      ELA = FEXP(LA)
      XA = ELA/(1.0D0 + ELA)
      ELB = FEXP(LB)
      XB = ELB/(1.0D0 + ELB)
      RTSHAT = DSQRT(XA*XB) * RTS
      UVSCALE = 0.25D0 * RTSHAT / DCOSH(0.7D0*YSTAR)
      ASCALE =  0.25D0 * RTSHAT / DCOSH((1.0D0 - XA)*YSTAR)
      BSCALE =  0.25D0 * RTSHAT / DCOSH((1.0D0 - XB)*YSTAR)
C
      MUUV = MAX(LAMBDA+0.01D0,AUV * UVSCALE)
      MUA  = MAX(2.24D0,ACO * ASCALE )
      MUB  = MAX(2.24D0,ACO * BSCALE )
      MUUV = MIN(MUUV,999.999D0)
      MUA  = MIN(MUA, 999.999D0)
      MUB  = MIN(MUB, 999.999D0)
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE JETDEFBORN(Y1,P2,Y2,PHI2,OK,LA,LB,YSTAR,SVALUE)
C
C  Input variables:
      DOUBLE PRECISION Y1,P2,Y2,PHI2
C  Output variables:
      LOGICAL OK
      DOUBLE PRECISION LA,LB,YSTAR,SVALUE
C
C  The output variables are JETVAR1,JETVAR2,JETVAR3,SVALUE, where
C  JETVAR1,JETVAR2,JETVAR3 are fed to BINIT and SETMU.  SVALUE is
C  the value of the coeffient of the delta and theta functions in
C  the jet-defining function.  (Normally, SVALUE = 1, obtained when
C  we treat the partons as indistinguishable.)
C
C  See JETDEF for defintions for this two jet inclusive cross section.
C
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      LOGICAL INBOUNDS
      DOUBLE PRECISION FEXP
      DOUBLE PRECISION XA,XB
C
      LA  = 0.0D0
      LB = 0.0D0
      YSTAR = 0.0D0
      OK = .FALSE.

C
      XA = (P2/RTS)*( FEXP( Y1) + FEXP( Y2) )
      XB = (P2/RTS)*( FEXP(-Y1) + FEXP(-Y2) )
C
      IF (XA.GE.1.0D0.OR.XB.GE.1.0D0) THEN
        RETURN
      ENDIF
C
      LA = DLOG( XA/(1.0D0 - XA) )
      LB = DLOG( XB/(1.0D0 - XB) )
      YSTAR = 0.5D0 * DABS(Y1-Y2)
C
      SVALUE = 1.0D0
C
      IF (  INBOUNDS(LA,LB,YSTAR) ) THEN
        OK = .TRUE.
      ENDIF
C
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE 
     >    JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,OK,LA,LB,YSTAR,SVALUE)
C
C  Input variables:
      INTEGER CASE
      DOUBLE PRECISION Y1,P2,Y2,PHI2,P3,Y3,PHI3
C  Output variables:
      LOGICAL OK
      DOUBLE PRECISION  LA,LB,YSTAR,SVALUE
C
C  The output variables are JETVAR1,JETVAR2,JETVAR3,SVALUE, where 
C  JETVAR1,JETVAR2,JETVAR3 are fed to BINIT and SETMU and SVALUE
C  is a factor contributed to the integrand.
C
C  The jet definition is contained in a function {\cal S}, which
C  is assumed to be a sum over CASE = 1,NJETDEF of terms, each of
C  which consists of
C
C  (1) delta functions setting JETVAR1,JETVAR2,JETVAR3 to some
C      function of the kinematic variables.
C  (2) theta functions, implemented by setting OK = .TRUE. if
C      the theta functions are satisfied (and the values of the
C      JETVARs are not too outrageous.
C  (3) a continuous function, SVALUE.  (Normally, SVALUE = 1, obtained
C      when we consider the partons as indistinguishable)
C
C  The parameter NJETDEF is set in the main program. (Normally it is 3.)
C
C  Let the two jets have rapidities YJ1 and YJ2 and transverse energies
C  PJ1 and PJ2 (computed according to the Snowmass convention).
C  When there are three jets in the event, the two jets with the largest
C  transverse energies are used. Note that we do not distinguish
C  between jets 1 and 2.
C
C  The jet variables are
C     LA = ln(XA/(1-XA)
C     LB = ln(XB/(1-XB)
C     YSTAR = |YJ1-YJ2|/2
C  where
C     XA = Sum_i (PJi/RTS) EXP( Yi)
C     XB = Sum_i (PJi/RTS) EXP(-Yi)
C  summed over i in jet.
C   
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
      DOUBLE PRECISION P1X,P1Y,P1,PHI1,PHI12,PHI23,PHI31,CONVERT
      DOUBLE PRECISION YJ1,YJ2
      DOUBLE PRECISION XA,XB
      DOUBLE PRECISION FEXP 
      LOGICAL INBOUNDS
      DOUBLE PRECISION PI
      PARAMETER ( PI = 3.141592654 D0)
C
      LA  = 0.0D0
      LB = 0.0D0
      YSTAR = 0.0D0
      OK = .FALSE.
C
C  Set S to 1
C
      SVALUE = 1.0D0
C
      P1X = - P2 * DCOS(PHI2) - P3 * DCOS(PHI3)
      P1Y = - P2 * DSIN(PHI2) - P3 * DSIN(PHI3)
      P1 = DSQRT(P1X**2 + P1Y**2)
      PHI1 = DATAN2(P1Y,P1X)
C
      PHI12 = CONVERT(PHI1 - PHI2)
      PHI23 = CONVERT(PHI2 - PHI3)
      PHI31 = CONVERT(PHI3 - PHI1)
C
C   CASE 1: 1 and 2 are the jets (so neither (1,3) nor (2,3) are jets)
C
      IF (CASE.EQ.1) THEN
       IF( ((Y3-Y1)**2 + PHI31**2).GE.( ((P1+P3)/P1*R)**2 )) THEN
       IF( ((Y3-Y2)**2 + PHI23**2).GE.( ((P2+P3)/P2*R)**2 )) THEN
         YJ1 = Y1
         YJ2 = Y2
C
         XA = ( P1 * FEXP( Y1) + P2 * FEXP( Y2) )/RTS
         XB = ( P1 * FEXP(-Y1) + P2 * FEXP(-Y2) )/RTS
         IF (XA.GE.1.0D0.OR.XB.GE.1.0D0) THEN
            RETURN
         ENDIF
         LA = DLOG( XA/(1.0D0 - XA) )
         LB = DLOG( XB/(1.0D0 - XB) )
         YSTAR = 0.5D0 * DABS(YJ1-YJ2)
C
         IF (  INBOUNDS(LA,LB,YSTAR) ) THEN
           OK = .TRUE.
         ENDIF
       ENDIF
       ENDIF
       RETURN
C
C   CASE 2: (1,3) and 2 are the jets
C
      ELSEIF (CASE.EQ.2) THEN
       IF( ((Y3-Y1)**2 + PHI31**2).LT.( ((P1+P3)/P1*R)**2 )) THEN
         YJ1 = (P1*Y1 + P3*Y3)/(P1+P3)
         YJ2 = Y2 
C
         XA = ( P1 * FEXP( Y1) + P2 * FEXP( Y2) + P3 * FEXP( Y3) )/RTS
         XB = ( P1 * FEXP(-Y1) + P2 * FEXP(-Y2) + P3 * FEXP(-Y3) )/RTS
         IF (XA.GE.1.0D0.OR.XB.GE.1.0D0) THEN
            RETURN
         ENDIF
         LA = DLOG( XA/(1.0D0 - XA) )
         LB = DLOG( XB/(1.0D0 - XB) )
         YSTAR = 0.5D0 * DABS(YJ1-YJ2)
C
         IF (  INBOUNDS(LA,LB,YSTAR) ) THEN
           OK = .TRUE.
         ENDIF
        ENDIF
       RETURN
C
C   CASE 3: (2,3) and 1 are the jets
C
      ELSEIF (CASE.EQ.3) THEN
       IF( ((Y3-Y2)**2 + PHI23**2).LT.( ((P2+P3)/P2*R)**2 )) THEN
         YJ2  = (P2*Y2 + P3*Y3)/(P2+P3)
         YJ1 = Y1
C
         XA = ( P1 * FEXP( Y1) + P2 * FEXP( Y2) + P3 * FEXP( Y3) )/RTS
         XB = ( P1 * FEXP(-Y1) + P2 * FEXP(-Y2) + P3 * FEXP(-Y3) )/RTS
         IF (XA.GE.1.0D0.OR.XB.GE.1.0D0) THEN
            RETURN
         ENDIF
         LA = DLOG( XA/(1.0D0 - XA) )
         LB = DLOG( XB/(1.0D0 - XB) )
         YSTAR = 0.5D0 * DABS(YJ1-YJ2)
C
         IF (  INBOUNDS(LA,LB,YSTAR) ) THEN
           OK = .TRUE.
         ENDIF
        ENDIF
        RETURN
C
      ELSE
        WRITE(NOUT,*) 'Fatal error in JETDEF'
        STOP
      ENDIF
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE BINIT(LA,LB,YSTAR,INTEGRAND)
C 
      DOUBLE PRECISION LA,LB,YSTAR,INTEGRAND
C
C  The arguments of binit are JETVAR1,JETVAR2,JETVAR3,INTEGRAND
C  where JETVAR1,JETVAR2,JETVAR3 are the jet parameters produced
C  by JETDEF and JETDEFBORN.
C
C  Adds result INTEGRAND to proper bin of TEMPARRAY after dividing
C  by a factor BORN and smearing with Gaussians in YSTAR
C  and in LA AND LB.
C
      INTEGER NL,NYSTAR
      DOUBLE PRECISION LGRID(1:16),YSTARGRID(1:16)
      COMMON /GRIDS/  LGRID,YSTARGRID,NL,NYSTAR
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
      DOUBLE PRECISION  LSMEAR,YSMEAR
      COMMON /SMEAR/ LSMEAR,YSMEAR
C
      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
      DOUBLE PRECISION BORN,BORNXSECT
      DOUBLE PRECISION INTEGRANDMOD
      DOUBLE PRECISION PREFACTOR,CY,CL
      DOUBLE PRECISION EXPONENT1,EXPONENT1BIS,EXPONENT2,EXPONENT3
      DOUBLE PRECISION EXPONENTIAL1,EXPONENTIAL1BIS
      DOUBLE PRECISION FACTOR1(1:16),FACTOR2(1:16),FACTOR3(1:16)
      DOUBLE PRECISION ADDITION
      DOUBLE PRECISION FEXP
      DOUBLE PRECISION ELA,ELB,XA,XB,ET,YBOOST,YJ1,YJ2
      DOUBLE PRECISION MUUV,MUA,MUB
      INTEGER IYSTAR,ILA,ILB
C
      DOUBLE PRECISION PI
      PARAMETER ( PI = 3.141592654 D0)
C
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
      IF (INTEGRAND.EQ.0.0D0) THEN
         RETURN
      ENDIF
C
C We divide by the Born cross section d sigma(B) / d LA dLB dYSTAR,
C including a jacobian factor
C           dYJ1 dYJ2 dET / dLA dLB dYSTAR = (1-XA)*(1-XB)*ET.
C
      ELA = FEXP(LA)
      XA = ELA/(1.0D0 + ELA)
      ELB = FEXP(LB)
      XB = ELB/(1.0D0 + ELB)
      ET = 0.5D0 * DSQRT(XA*XB) * RTS / DCOSH(YSTAR)
      YBOOST = 0.5D0 * DLOG(XA/XB)
      YJ1 = YBOOST + YSTAR
      YJ2 = YBOOST - YSTAR
C
      CALL SETMU(LA,LB,YSTAR,MUUV,MUA,MUB)
      BORNXSECT = (1-XA)*(1-XB)*ET*BORN(YJ1,ET,YJ2,MUUV,MUA,MUB)
      IF (BORNXSECT.EQ.0.0D0) THEN
         RETURN
      ENDIF
C
      INTEGRANDMOD = INTEGRAND /BORNXSECT
C
C The values of YSTAR are by definition positive.  However, we
C also add phantom events with YSTAR --> - YSTAR so that we are
C smearing a distribution that is smooth near YSTAR = 0.  We also 
C AVERAGE over contributions with XA interchanged with XB to get
C better statistics.  This is allowed because of P or CP symmetry.
C
      PREFACTOR = 0.5D0/(DSQRT(2.0D0*PI)**3 * YSMEAR * LSMEAR**2 )
      CY =  0.5D0 /YSMEAR**2
      CL =  0.5D0 /LSMEAR**2
C
C
      DO IYSTAR = 1,NYSTAR
        EXPONENT1 = CY * (YSTAR - YSTARGRID(IYSTAR))**2
        IF( EXPONENT1.LT.20 ) THEN
          EXPONENTIAL1 = DEXP(- EXPONENT1 )
        ELSE
          EXPONENTIAL1 = 0.0D0
        ENDIF
        EXPONENT1BIS = CY * (YSTAR + YSTARGRID(IYSTAR))**2
        IF( EXPONENT1BIS.LT.20 ) THEN
          EXPONENTIAL1BIS = DEXP(- EXPONENT1BIS )
        ELSE
          EXPONENTIAL1BIS = 0.0D0
        ENDIF
        FACTOR1(IYSTAR) = EXPONENTIAL1 + EXPONENTIAL1BIS
      ENDDO
C
      DO ILA = 1,NL
        EXPONENT2 = CL * (LA - LGRID(ILA))**2
        IF( EXPONENT2.LT.20 ) THEN
          FACTOR2(ILA) = DEXP(- EXPONENT2 )
        ELSE
          FACTOR2(ILA) = 0.0D0
        ENDIF
      ENDDO
C
      DO ILB = 1,NL
        EXPONENT3 = CL * (LB - LGRID(ILB))**2
        IF( EXPONENT3.LT.20 ) THEN
          FACTOR3(ILB) = DEXP(- EXPONENT3 )
        ELSE
          FACTOR3(ILB) = 0.0D0
        ENDIF
      ENDDO
C
      DO IYSTAR = 1,NYSTAR
      DO ILA = 1,NL
      DO ILB = 1,NL
C
      ADDITION = INTEGRANDMOD * PREFACTOR 
     >           * FACTOR1(IYSTAR) * FACTOR2(ILA) * FACTOR3(ILB)
C
      TEMPARRAY(ILA,ILB,IYSTAR) = TEMPARRAY(ILA,ILB,IYSTAR)
     >                            + ADDITION
      TEMPARRAY(ILB,ILA,IYSTAR) = TEMPARRAY(ILB,ILA,IYSTAR)
     >                            + ADDITION
C
      ENDDO
      ENDDO
      ENDDO
C
      RETURN
      END
C C***********************************************************************
C***********************************************************************
C
      DOUBLE PRECISION FUNCTION BORN(Y1,P2,Y2,MUUV,MUA,MUB)
C 
C
C     Calculates the Born cross section d sigma / d Y1 d P2 d Y2.
C     Note: we remove a factor of 1/2 from the normalization compared
C     to that in INTEGRATE2TO2: we want d sigma / dP2 d y1 dy2, not
C     1/(2!) times this.  We also multiply by 2 PI so that we get
C     d sigma / dP2 d y1 dy2 instead of d sigma / dP2 d y1 dy2 d phi2.
C     DES, 29 August 1993; 29 December 1993,28 January 1994.
C
      DOUBLE PRECISION Y1,P2,Y2,MUUV,MUA,MUB
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C  Switches
      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   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)
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   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 D0,ALPHAS
      DOUBLE PRECISION L
      DOUBLE PRECISION PSI4,PSITILDE4,SIGN
      DOUBLE PRECISION PARTON,WNUM
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
      BORN = 0.0D0
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
C Note: we remove a factor of 1/2 from D0 compared to the form
C in INTEGRATE2TO2 because we want d sigma / dP2 d y1 dy2, not
C 1/(2!) times this.  We also multiply by 2 PI so that we get
C d sigma / dP2 d y1 dy2 instead of d sigma / dP2 d y1 dy2 d phi2.
C
      D0 = 2.0D0*PI * ALPHAS(MUUV)**2 * P2 /RTS**4
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))
      ENDDO
C
      IF (PPBAR) THEN
       DO AB = -NFL,NFL
        PRTNB(AB) = PARTON(-AB,XB,MUB)/(XB*WNUM(AB))
       ENDDO
      ELSE
       DO AB = -NFL,NFL
        PRTNB(AB) = PARTON(AB,XB,MUB)/(XB*WNUM(AB))
       ENDDO
      ENDIF
C
C  ----- Initialize sum:
C
      BORN = 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)
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
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)
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
      BORN = BORN + D0 * L * SIGN * PSI4
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
      RETURN
      END
C
C***********************************************************************
C**************            End of File               *******************
C***********************************************************************