C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        JET version 4_0
C
C                Version to read the output file for
C              calculation of the two jet cross section
C                     d sigma/ d XA d XB d YSTAR
C
C                           14 June 1996
C
C                 S.D. Ellis, Z. Kunszt, D.E. Soper
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                      14 June 1996: Incorporates generic parton.f 
C                                    version 1.1, made available
C                                    on World Wide Web
C
C 
C***********************************************************************
C************************BEGIN MAIN PROGRAM*****************************
C
      PROGRAM JET
C
      INTEGER NFL
      DOUBLE PRECISION NFLAVOR,LAMBDA
      COMMON /QCDPARS/  NFLAVOR,LAMBDA,NFL
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
      LOGICAL STYPEA,STYPEB,STYPEC,STYPED
      COMMON /SWTYPE/ STYPEA,STYPEB,STYPEC,STYPED
C
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
C   Variables gathered by READDATA.
C
      DOUBLE PRECISION KM1ARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION ERRARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION SMEARARRAY(1:16,1:16,1:16)
      COMMON /ARRAYS/ KM1ARRAY,ERRARRAY,SMEARARRAY
C
      DOUBLE PRECISION LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR
      INTEGER NL,NYSTAR
      COMMON /GRIDINFO/ LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR,NL,NYSTAR
C
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
      LOGICAL PPBAR
      COMMON /BEAMS/ PPBAR
C
      DOUBLE PRECISION PARTON,PRTN
      CHARACTER*20 PARTONFILE
      COMMON /PARTONCHOICE/ PARTONFILE
C
      DOUBLE PRECISION AUV,ACO
      COMMON /MUCOEFS/ AUV,ACO
C
      DOUBLE PRECISION  LSMEAR,YSMEAR
      COMMON /SMEAR/ LSMEAR,YSMEAR
C
      INTEGER SEED,NRENO
      COMMON /RENOINFO/ SEED,NRENO
C
C   Ranges for output
C
      DOUBLE PRECISION OUTYSTARMAX,OUTXMIN,OUTXMAX
      INTEGER OUTNYSTAR,OUTNX
C
C   Results
C
      DOUBLE PRECISION LA,LB,YSTAR,XA,XB,YBOOST,Y1,Y2,P2
      DOUBLE PRECISION KFACTOR,KM1,KERR,SMEAR,BORN,BORNXSECT,XSECT
      DOUBLE PRECISION MUUV,MUA,MUB
C
C   Integers for DO loops
C
      INTEGER IYSTAR,ILA,ILB
C
      DOUBLE PRECISION FEXP
      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  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  Lambda_MSbar^(NFL): set in initialization of PARTON
C
C  Input and output units.
C
      NIN  = 5
      NOUT = 6
      NWRT = 6
C
C  Switches (should all be .TRUE.)
C
      STYPEA    = .TRUE.
      STYPEB    = .TRUE.
      STYPEC    = .TRUE.
      STYPED    = .TRUE.
C
C  Get the data; generate SMEARARRAY.
C
      CALL READDATA
      CALL MAKESMEARARRAY
C
C  Get desired ranges from user:
C
      WRITE(NOUT,*)' Give maximum YSTAR and number of ystar points'
      WRITE(NOUT,*)' Maximum Ystar .LE.',YSTARMAX
      READ(NIN,*) OUTYSTARMAX,OUTNYSTAR
      WRITE(NOUT,*)OUTYSTARMAX,OUTNYSTAR
      WRITE(NOUT,*)' Give min and max X and number of X points'
      WRITE(NOUT,*)' Xs must be in the range',
     >               1.0D0/(1.0D0 + FEXP(-LMIN)),
     >               1.0D0/(1.0D0 + FEXP(-LMAX))
      READ(NIN,*) OUTXMIN,OUTXMAX,OUTNX
      WRITE(NOUT,*) OUTXMIN,OUTXMAX,OUTNX
C
C   WRITE THE PARAMETERS
C
      CALL DAYTIME(DATE)
      WRITE(NOUT,101)DATE
101   FORMAT(//
     >    ,' Program READJET version 4.0 ',A,/,80('-'))
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(80('-'))
C
      IF (PPBAR) THEN
        WRITE(NOUT,*) 'Proton-antiproton collisions'
      ELSE
        WRITE(NOUT,*) 'Proton-proton collisions'
      ENDIF
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 used',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
      WRITE(NOUT,735)
735     FORMAT(/,' Data Grid: ')
      WRITE(NOUT,736)LMIN,LMAX,DELTAL
736     FORMAT('  LA, LB from ',F8.3,' to ',F8.3,' in steps of ',F8.3)
      WRITE(NOUT,737)1.0D0/(1.0D0 + FEXP(-LMIN)),
     >               1.0D0/(1.0D0 + FEXP(-LMAX))
737     FORMAT('  XA, XB  from ',F12.6,' to ',F8.3)
      WRITE(NOUT,738)YSTARMAX,DELTAYSTAR
738     FORMAT('  YSTAR from 0 to ',F8.3,' in steps of ',F8.3)
C
C
C----------------------------------------------
C
      WRITE(NOUT,*)' '
      WRITE(NOUT,601)
601    FORMAT('----------------------------------------',
     >        '----------------------------------------')
      WRITE(NOUT,*)' '
      WRITE(NOUT,602)
602    FORMAT('                                    RESULTS')
      WRITE(NOUT,*)' '
C
      DO IYSTAR = 1,OUTNYSTAR
C
      WRITE(NOUT,601)
      WRITE(NOUT,603)
603    FORMAT('Ystar   XA     XB     ET     Y1    Y2   ',
     >        '   Xsect    Kerror Ksmear     Born      ')
      WRITE(NOUT,*)' '
C  
C---------  Loop over grid points, printing results for each.
C
      DO ILA = 1,OUTNX
      DO ILB = ILA,OUTNX
C
C  Kinematics.
C
      LA = (OUTNX - ILA)/(OUTNX-1.0D0) * DLOG(OUTXMIN/(1.0D0 - OUTXMIN))
     >     +(ILA - 1)/(OUTNX-1.0D0) * DLOG(OUTXMAX/(1.0D0 - OUTXMAX))
      LB = (OUTNX - ILB)/(OUTNX-1.0D0) * DLOG(OUTXMIN/(1.0D0 - OUTXMIN))
     >     +(ILB - 1)/(OUTNX-1.0D0) * DLOG(OUTXMAX/(1.0D0 - OUTXMAX))
      YSTAR = (IYSTAR - 1)/(OUTNYSTAR - 1.0D0) * OUTYSTARMAX
C
      XA = DEXP(LA) /( 1.0D0 + DEXP(LA) )
      XB = DEXP(LB) /( 1.0D0 + DEXP(LB) )
      YBOOST = 0.5D0 * DLOG(XA/XB)
      Y1 = YBOOST + YSTAR
      Y2 = YBOOST - YSTAR
      P2 =  0.5D0 * DSQRT(XA * XB) * RTS / DCOSH(YSTAR)
C
C  The K-factor, with its error.
C
      CALL INTERPOLATE(LA,LB,YSTAR,KM1,KERR,SMEAR)
      KFACTOR = KM1 + 1.0D0
C
C  We calculate the Born cross section, defining parton 1 to
C  be the parton with the greater value of y. We multiply by
C  the jacobian dY1 dY2 dP2 / dXA dXB d YSTAR = P2/(XA XB) so
C  that the Born cross section is d sigma / dXA dXB d YSTAR.
C
      CALL SETMU(LA,LB,YSTAR,MUUV,MUA,MUB)
      BORNXSECT = (P2/(XA*XB))* BORN(Y1,P2,Y2,MUUV,MUA,MUB)
C
C   Cross section is K-factor times Born cross section. 
C
      XSECT    = KFACTOR   * BORNXSECT 
C
C   Change units from GeV^-2 to nb
C
      XSECT     = NBGEV2 * XSECT
      BORNXSECT = NBGEV2 * BORNXSECT
C
C   Write results
C
      WRITE(NOUT,15)YSTAR,XA,XB,P2,Y1,Y2,XSECT,KERR,SMEAR,BORNXSECT 
15    FORMAT(F5.3,2F7.3,F7.1,2F6.2,G12.3,2F7.2,G12.3)
C
      ENDDO
      WRITE(*,*)' '
      ENDDO
      ENDDO
C
C----------------------------------------------
C 
      CALL DAYTIME(DATE)
      WRITE(NOUT,103)DATE
103   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
      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***********************************************************************
C
      SUBROUTINE READDATA
C
C Reads this information from file, puts information in common blocks.
C
      DOUBLE PRECISION KM1ARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION ERRARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION SMEARARRAY(1:16,1:16,1:16)
      COMMON /ARRAYS/ KM1ARRAY,ERRARRAY,SMEARARRAY
C
      DOUBLE PRECISION LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR
      INTEGER NL,NYSTAR
      COMMON /GRIDINFO/ LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR,NL,NYSTAR
C
      DOUBLE PRECISION RTS,R
      COMMON /PHYSICS/ RTS,R
C
C  Switch for p-pbar collisions (if true) or p-p collisions (if false)
C
      LOGICAL PPBAR
      COMMON /BEAMS/ PPBAR
C
C   Parton distribution function
C
      CHARACTER*20 PARTONFILE
      COMMON /PARTONCHOICE/ PARTONFILE
C
      DOUBLE PRECISION AUV,ACO
      COMMON /MUCOEFS/ AUV,ACO
C
      DOUBLE PRECISION  LSMEAR,YSMEAR
      COMMON /SMEAR/ LSMEAR,YSMEAR
C
      INTEGER SEED,NRENO
      COMMON /RENOINFO/ SEED,NRENO
C
C Other variables.
C
C IO unit numbers
C
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
C Min and max values of X.
C
      DOUBLE PRECISION X1,XN
C
C The data as read, with XA and XB translated to LA and LB.
C
      DOUBLE PRECISION YSTAR,XA,XB,KFACTORM1,KERR
      DOUBLE PRECISION LA,LB
C
C Grids for YSTAR, LA, and LB
C
      DOUBLE PRECISION YSTARGRID(1:16),LGRID(1:16)
C
C Name for input file and variable name for the input
C
      CHARACTER*10 FILENAME,VARNAME
      CHARACTER*3 DUMMY
C
C  Code to indentify version of the program that wrote file
C
      INTEGER VERSION,SHOULDBE
C
C   Identification of pp or ppbar
C
      CHARACTER*29 BEAMTYPE
C
C Variables for DO loops
C
      INTEGER I,IYSTAR,ILA,ILB
C
C  Begin ---------
C
C  Output file to read:
C
      WRITE(NOUT,*)' Give name of file to read'
      READ(NIN,*) FILENAME
C
      OPEN(100,FILE=FILENAME,STATUS='OLD',ERR=99)
C
      READ(100,*)
      READ(100,*)
      READ(100,*)
C
C  Switch for p-pbar collisions (if true) or p-p collisions (if false)
C
      READ(100,204)BEAMTYPE
      IF (BEAMTYPE.EQ.' Proton-antiproton collisions') THEN
        PPBAR = .TRUE.
      ELSEIF (BEAMTYPE.EQ.' Proton-proton collisions    ') THEN
        PPBAR = .FALSE.
      ELSE
        VARNAME = '  BEAMTYPE'
        GOTO 199
      ENDIF
      READ(100,*)
      READ(100,202)VARNAME,DUMMY,VERSION
      IF(VARNAME.NE.'   VERSION') GOTO 199
      SHOULDBE = 940923
      IF(VERSION.NE.SHOULDBE) GOTO 199
      READ(100,201)VARNAME,DUMMY,RTS
      IF(VARNAME.NE.'       RTS') GOTO 199
      READ(100,201)VARNAME,DUMMY,R
      IF(VARNAME.NE.'         R') GOTO 199
      READ(100,201)VARNAME,DUMMY,X1
      IF(VARNAME.NE.'        X1') GOTO 199
      READ(100,201)VARNAME,DUMMY,XN
      IF(VARNAME.NE.'        XN') GOTO 199
      READ(100,202)VARNAME,DUMMY,NL
      IF(VARNAME.NE.'        NL') GOTO 199
      READ(100,201)VARNAME,DUMMY,YSTARMAX
      IF(VARNAME.NE.'  YSTARMAX') GOTO 199
      READ(100,202)VARNAME,DUMMY,NYSTAR
      IF(VARNAME.NE.'    NYSTAR') GOTO 199
      READ(100,201)VARNAME,DUMMY,LSMEAR
      IF(VARNAME.NE.'    LSMEAR') GOTO 199
      READ(100,201)VARNAME,DUMMY,YSMEAR
      IF(VARNAME.NE.'    YSMEAR') GOTO 199
      READ(100,203)VARNAME,DUMMY,PARTONFILE
      IF(VARNAME.NE.'   PARTONS') GOTO 199
      READ(100,201)VARNAME,DUMMY,AUV
      IF(VARNAME.NE.'       AUV') GOTO 199
      READ(100,201)VARNAME,DUMMY,ACO
      IF(VARNAME.NE.'       ACO') GOTO 199
      READ(100,202)VARNAME,DUMMY,SEED
      IF(VARNAME.NE.'      SEED') GOTO 199
      READ(100,202)VARNAME,DUMMY,NRENO
      IF(VARNAME.NE.'     NRENO') GOTO 199
201   FORMAT(A10,A3,G20.10)
202   FORMAT(A10,A3,I10)
203   FORMAT(A10,A3,A10)
204   FORMAT(A29)
      READ(100,*)
      READ(100,*)
      READ(100,*)
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
        LMIN = DLOG(X1/(1.0D0 - X1))
        LMAX = DLOG(XN/(1.0D0 - XN))
        DELTAL = ( LMAX - LMIN ) /(NL-1)
      ELSE
        DELTAL = 1.0D0
      ENDIF
      DO  I = 1,NL
         LGRID(I) = LMIN + (I-1) * DELTAL
      ENDDO
C
C  Read data.
C
      DO IYSTAR = 1,NYSTAR
      DO ILA = 1,NL
      DO ILB = 1,NL
C
      READ(100,16)YSTAR,XA,XB,KFACTORM1,KERR 
16    FORMAT(5G20.10)
      KM1ARRAY(ILA,ILB,IYSTAR) = KFACTORM1
      ERRARRAY(ILA,ILB,IYSTAR) = KERR
      LA = LGRID(ILA)
      LB = LGRID(ILB)
      IF( DABS(YSTAR - YSTARGRID(IYSTAR)) .GT. 1D-6 ) GOTO 299
      IF( DABS(DEXP(LA) /( 1.0D0 + DEXP(LA) )/XA - 1.0D0) .GT. 1D-6 )
     >   GOTO 299
      IF( DABS(DEXP(LB) /( 1.0D0 + DEXP(LB) )/XB - 1.0D0) .GT. 1D-6 ) 
     >   GOTO 299
C
      ENDDO
      ENDDO
      ENDDO
C
      CLOSE(100)
C
      RETURN
99    WRITE(NOUT,*)'Error opening data file for input'
      STOP
C
199   WRITE(NOUT,*)'Missmatch reading variable',VARNAME
      STOP
C
299   WRITE(NOUT,*)'Missmatch reading YSTAR,XA,XB'
      WRITE(NOUT,*) YSTAR, YSTARGRID(IYSTAR)
      WRITE(NOUT,*) XA,DEXP(LA) /( 1.0D0 + DEXP(LA) ),LA
      WRITE(NOUT,*) XB,DEXP(LB) /( 1.0D0 + DEXP(LB) ),LB
      STOP
C
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE MAKESMEARARRAY 
C
      DOUBLE PRECISION KM1ARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION ERRARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION SMEARARRAY(1:16,1:16,1:16)
      COMMON /ARRAYS/ KM1ARRAY,ERRARRAY,SMEARARRAY
C
      DOUBLE PRECISION LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR
      INTEGER NL,NYSTAR
      COMMON /GRIDINFO/ LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR,NL,NYSTAR
C
      DOUBLE PRECISION  LSMEAR,YSMEAR
      COMMON /SMEAR/ LSMEAR,YSMEAR
C
      DOUBLE PRECISION D2KDYSTAR2,D2KDLA2,D2KDLB2 
      INTEGER IYSTAR,ILA,ILB
C
C
      DO IYSTAR = 2,NYSTAR-1
      DO ILA = 2,NL-1
      DO ILB = 2,NL-1
C
      D2KDYSTAR2 = (  KM1ARRAY(ILA,ILB,IYSTAR+1)
     >              + KM1ARRAY(ILA,ILB,IYSTAR-1)
     >              - 2.0D0*KM1ARRAY(ILA,ILB,IYSTAR) )/DELTAYSTAR**2
      D2KDLA2    = (  KM1ARRAY(ILA+1,ILB,IYSTAR)
     >              + KM1ARRAY(ILA-1,ILB,IYSTAR)
     >              - 2.0D0*KM1ARRAY(ILA,ILB,IYSTAR) )/DELTAL**2
      D2KDLB2    = (  KM1ARRAY(ILA,ILB+1,IYSTAR) 
     >              + KM1ARRAY(ILA,ILB-1,IYSTAR)
     >              - 2.0D0*KM1ARRAY(ILA,ILB,IYSTAR) )/DELTAL**2
C
      SMEARARRAY(ILA,ILB,IYSTAR) =
     >                  - 0.5D0 * D2KDYSTAR2 * YSMEAR**2
     >                  - 0.5D0 * D2KDLA2    * LSMEAR**2
     >                  - 0.5D0 * D2KDLB2    * LSMEAR**2
C
      ENDDO
      ENDDO
      ENDDO
C
      DO ILA = 2,NL-1
      DO ILB = 2,NL-1
C
      SMEARARRAY(ILA,ILB,1)      = SMEARARRAY(ILA,ILB,2)
      SMEARARRAY(ILA,ILB,NYSTAR) = SMEARARRAY(ILA,ILB,NYSTAR-1)
C
      ENDDO
      ENDDO
C
      DO IYSTAR = 1,NYSTAR
      DO ILB = 2,NL-1
C
      SMEARARRAY( 1,ILB,IYSTAR) = SMEARARRAY(   2,ILB,IYSTAR)
      SMEARARRAY(NL,ILB,IYSTAR) = SMEARARRAY(NL-1,ILB,IYSTAR)
C
      ENDDO
      ENDDO
C
      DO IYSTAR = 1,NYSTAR
      DO ILA = 1,NL
C
      SMEARARRAY(ILA, 1,IYSTAR) = SMEARARRAY(ILA,   2,IYSTAR)
      SMEARARRAY(ILA,NL,IYSTAR) = SMEARARRAY(ILA,NL-1,IYSTAR)
C
      ENDDO
      ENDDO
C
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE INTERPOLATE(LA,LB,YSTAR,KM1,KERR,SMEAR)
C  In:
      DOUBLE PRECISION LA,LB,YSTAR
C  Out:
      DOUBLE PRECISION KM1,KERR,SMEAR
C
C  Interpolates from the tables KM1ARRAY and ERRARRAY to give the 
C  interpolated values of KM1 and ERR at the point specified by the 
C  variables LA,LB,YSTAR. The tables are taken from the common
C  block /ARRAYS/. The grids for LA,LB,YSTAR are constructed from
C  the data in the common block /GRIDINFO/.
C
      INTEGER NIN,NOUT,NWRT
      COMMON /IOUNIT/ NIN,NOUT,NWRT
C
      DOUBLE PRECISION KM1ARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION ERRARRAY(1:16,1:16,1:16)
      DOUBLE PRECISION SMEARARRAY(1:16,1:16,1:16)
      COMMON /ARRAYS/ KM1ARRAY,ERRARRAY,SMEARARRAY
C
      DOUBLE PRECISION LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR
      INTEGER NL,NYSTAR
      COMMON /GRIDINFO/ LMIN,LMAX,DELTAL,YSTARMAX,DELTAYSTAR,NL,NYSTAR
C
      INTEGER ILA,ILB,IYSTAR
      DOUBLE PRECISION ZLA,ZLB,ZYSTAR
      DOUBLE PRECISION XXLA,XXLB,XXYSTAR
      DOUBLE PRECISION YYLA,YYLB,YYYSTAR
C
C  Range checking. If we are just slightly out of range, we
C  put the variables into the range.
C
      IF ( LA.LE.LMIN ) THEN
        IF (LA.LT.LMIN - 1.0D-6) THEN
          WRITE(NOUT,*)'FXSECT called out of range, LA =',LA
          STOP
        ENDIF
        LA = LMIN + 1.0D-12
      ENDIF
      IF ( LA.GE.LMAX ) THEN
        IF (LA.GT.LMAX + 1.0D-6) THEN
          WRITE(NOUT,*)'FXSECT called out of range, LA =',LA
          STOP
        ENDIF
        LA = LMAX - 1.0D-12
      ENDIF
C
      IF ( LB.LE.LMIN ) THEN
        IF (LB.LT.LMIN - 1.0D-6) THEN
          WRITE(NOUT,*)'FXSECT called out of range, LB =',LB
          STOP
        ENDIF
        LB = LMIN + 1.0D-12
      ENDIF
      IF ( LB.GE.LMAX ) THEN
        IF (LB.GT.LMAX + 1.0D-6) THEN
          WRITE(NOUT,*)'FXSECT called out of range, LB =',LB
          STOP
        ENDIF
        LB = LMAX - 1.0D-12
      ENDIF
C
      IF ( YSTAR.LE.0.0D0 ) THEN
        IF (YSTAR.LT. -1.0D-6) THEN
          WRITE(NOUT,*)'FXSECT called out of range, YSTAR =',YSTAR
          STOP
        ENDIF
        YSTAR = 1.0D-12
      ENDIF
      IF ( YSTAR.GE.YSTARMAX ) THEN
        IF (YSTAR.GT.YSTARMAX + 1.0D-6) THEN
          WRITE(NOUT,*)'FXSECT called out of range, YSTAR =',YSTAR
          STOP
        ENDIF
        YSTAR = YSTARMAX - 1.0D-12
      ENDIF
C
C  Calculate the factors for interpolation. The grid points are at
C  integer values of ZLA, ZLB, and ZYSTAR from 1 to NL or NYSTAR.
C
      ZLA = (LA - LMIN)/DELTAL + 1.0D0
      ILA = INT(ZLA)
      XXLA = ZLA - DBLE(ILA)
      YYLA = 1.0D0 - XXLA
C
      ZLB = (LB - LMIN)/DELTAL + 1.0D0
      ILB = INT(ZLB)
      XXLB = ZLB - DBLE(ILB)
      YYLB = 1.0D0 - XXLB
C
      ZYSTAR = YSTAR/DELTAYSTAR + 1.0D0
      IYSTAR = INT(ZYSTAR)
      XXYSTAR = ZYSTAR - DBLE(IYSTAR)
      YYYSTAR = 1.0D0 - XXYSTAR
C
C  The interpolation.
C
      KM1 = YYLA*YYLB*YYYSTAR*KM1ARRAY(ILA,  ILB,  IYSTAR)
     >   +  XXLA*YYLB*YYYSTAR*KM1ARRAY(ILA+1,ILB,  IYSTAR)
     >   +  YYLA*XXLB*YYYSTAR*KM1ARRAY(ILA,  ILB+1,IYSTAR)
     >   +  XXLA*XXLB*YYYSTAR*KM1ARRAY(ILA+1,ILB+1,IYSTAR)
     >   +  YYLA*YYLB*XXYSTAR*KM1ARRAY(ILA,  ILB,  IYSTAR+1)
     >   +  XXLA*YYLB*XXYSTAR*KM1ARRAY(ILA+1,ILB,  IYSTAR+1)
     >   +  YYLA*XXLB*XXYSTAR*KM1ARRAY(ILA,  ILB+1,IYSTAR+1)
     >   +  XXLA*XXLB*XXYSTAR*KM1ARRAY(ILA+1,ILB+1,IYSTAR+1)
C
      KERR = YYLA*YYLB*YYYSTAR*ERRARRAY(ILA,  ILB,  IYSTAR)
     >   +   XXLA*YYLB*YYYSTAR*ERRARRAY(ILA+1,ILB,  IYSTAR)
     >   +   YYLA*XXLB*YYYSTAR*ERRARRAY(ILA,  ILB+1,IYSTAR)
     >   +   XXLA*XXLB*YYYSTAR*ERRARRAY(ILA+1,ILB+1,IYSTAR)
     >   +   YYLA*YYLB*XXYSTAR*ERRARRAY(ILA,  ILB,  IYSTAR+1)
     >   +   XXLA*YYLB*XXYSTAR*ERRARRAY(ILA+1,ILB,  IYSTAR+1)
     >   +   YYLA*XXLB*XXYSTAR*ERRARRAY(ILA,  ILB+1,IYSTAR+1)
     >   +   XXLA*XXLB*XXYSTAR*ERRARRAY(ILA+1,ILB+1,IYSTAR+1)
C
      SMEAR = YYLA*YYLB*YYYSTAR*SMEARARRAY(ILA,  ILB,  IYSTAR)
     >   +    XXLA*YYLB*YYYSTAR*SMEARARRAY(ILA+1,ILB,  IYSTAR)
     >   +    YYLA*XXLB*YYYSTAR*SMEARARRAY(ILA,  ILB+1,IYSTAR)
     >   +    XXLA*XXLB*YYYSTAR*SMEARARRAY(ILA+1,ILB+1,IYSTAR)
     >   +    YYLA*YYLB*XXYSTAR*SMEARARRAY(ILA,  ILB,  IYSTAR+1)
     >   +    XXLA*YYLB*XXYSTAR*SMEARARRAY(ILA+1,ILB,  IYSTAR+1)
     >   +    YYLA*XXLB*XXYSTAR*SMEARARRAY(ILA,  ILB+1,IYSTAR+1)
     >   +    XXLA*XXLB*XXYSTAR*SMEARARRAY(ILA+1,ILB+1,IYSTAR+1)
C
      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
      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,1000D0)
      MUA  = MIN(MUA,1000D0)
      MUB  = MIN(MUB,1000D0)
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***********************************************************************
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***********************************************************************
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 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***********************************************************************
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***********************************************************************
C***********************************************************************