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***********************************************************************