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