C23456789012345678901234567890123456789012345678901234567890123456789012 C C JET version 3.4 C C Version to calculate the two jet cross section C d sigma/ d MJJ d y* C C S.D. Ellis, Z. Kunszt, D.E. Soper C C 13 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 C Function TRIAL revised to allow variation of RTS: 22 June 1992 C Size of DATE enlarged for new Sun operating system: 25 May 1994 C Small modification of output format 7 June 1994, 12 July 1994 C Smearing distances changed 8 June 1994, 12 July 1994 C C Version 3_3, 9 June 1996 C Uses generic parton program from jet 4.0 C Version 3.4, 13 June 1996 C Incorporates generic parton.f version 1.1 C Made available on World Wide Web C C C This is the master program, which is to be linked with the C subroutines in the file jetsubs3_4.f. This version of the C master program calculates the one jet cross section C C d sigma/ d MJJ d YSTAR C C where MJJ is the jet-jet mass and YSTAR = |Y_jet1 - Y_jet2|/2. C The cross section is integrated over a window of YCM = C |Y_jet1 + Y_jet2|/2. The cross section is given on a grid of C points MJJ, and YSTAR. At the end, it also integrates this over C a requested bin of MJJ. C Other master programs can calculate other C quantities of interest when linked with jetsubs3_4.f. C C Results from the version of this program that calculates C the one jet inclusive cross section were first presented in C C THE ONE JET INCLUSIVE CROSS-SECTION AT ORDER ALPHA-S**3 QUARKS C AND GLUONS. C By Stephen D. Ellis (Washington U., Seattle), Zoltan Kunszt C (Zurich,ETH), Davison E. Soper (Oregon U) C Phys.Rev.Lett.64:2121,1990. C C Results for the two jet cross section from this program were C published in C C TWO JET PRODUCTION IN HADRON COLLISIONS AT ORDER ALPHA-S**3 IN QCD. C By Stephen D. Ellis (Washington U., Seattle), Zoltan Kunszt C (Zurich, ETH), Davison E. Soper (Oregon U.) C Phys.Rev.Lett.69:1496-1499,1992. C C A technical description of the algorithm used in this program C is given in the paper C C CALCULATION OF JET CROSS-SECTIONS IN HADRON COLLISIONS C AT ORDER ALPHA-S**3. C By Zoltan Kunszt (Zurich, ETH), Davison E. Soper (Oregon U.). C Phys.Rev.D46:192-221,1992. C C An earlier version was described in C THE ONE JET INCLUSIVE CROSS-SECTION AT ORDER ALPHA-S**3. C 1. GLUONS ONLY. C By Stephen D. Ellis (Washington U., Seattle), Zoltan Kunszt C (Zurich, ETH), Davison E. Soper (Oregon U.). C Phys.Rev.D40:2188,1989. C C ----- C C There are three 'answers,' determined as follows. C C (1) For each grid point (MJJ,Y*), the program calculates the C 'SmearedXSECT,' S(MJJ,Y*), given in terms of the true cross section, C s(p*,y*), by C C S(MJJ,Y*) C C ~ f(MJJ,Y*) * \int dln(p) dy h[ln(MJJ)-ln(p*),Y* - y*] C * s(p*,y*)/f(p*,y*). C C where h is a gaussian smearing function and f (called TRIAL(MJJ,Y*) C in the program) is a fairly good fit to the cross section. (In C this way, we are smearing a function that is almost constant, so C smearing has little effect.) The calculation is not exact, hence C the "~", because of integration errors. C C (2) We fit a function, Poly(MJJ,Y*), which is a polynomial in C ln(MJJ) and Y*, to S(MJJ,Y*)/f(MJJ,Y*): C C S(MJJ,Y*) ~ f(MJJ,Y*) * Poly(MJJ,Y*). C C (3) We calculate the unsmeared polynomial poly(p*,y*) such that if C poly is smeared as above the result is (exactly) Poly: C C Poly(MJJ,Y*) = \int dln(p) dy h[ln(MJJ)-ln(p*),Y* - y*] poly(p*,y*). C C (4) The 'Unsmeared Fit XSECT' is C C f(MJJ,Y*) * poly(MJJ,Y*). C C (5) The 'Corrected XSECT' is C C S(MJJ,Y*) + f(MJJ,Y*)*( poly[MJJ,Y*] - Poly[MJJ,Y*] ). C C Each of the three answers has its advantages. To the extent that C the f(MJJ,Y*) times a polynomial of order, say 3 x 3, is able to C fit the true cross section very well, the 'Unsmeared Fit XSECT' C should be very good because it makes use of a wide range of C random points (p*,y*) chosen by the Monte Carlo program. On C the other hand, to the extent that the statistical errors C from the Monte Carlo integration are small, the 'Corrected C XSECT' should be good because it is not very dependent on the C quality of the fit (because the 'unsmearing correction' is C rather small). C C If f(MJJ,Y*) provides a very good fit to the cross section, then C the smearing correction will be smaller than the integration errors, C and the unsmearing proceedure is not needed. In this case, one can C just take the SmearedXSECT. C C To get a very good answer, one just has to use a large number of C Monte Carlo points and to either set the smearing distances C in the program to very small numbers or use a function f(MJJ,Y*) C that is a very good fit to S(MJJ,Y*). Also, it is possible to change C the number of parameters used in the fit. The default is JMAX = 2, C KMAX = 3 but with small statistical errors one could raise these. C C After this, the program integrates the unsmeared fit cross section C over the desired bin of MJJ for a set of Y* points. This is a simple C numerical integration using the fit function. The user should check C that the fit is accurate by comparing 'Unsmeared Fit XSECT' to C 'Corrected XSECT'. C C This program was mostly developed on a Unix (Sun) workstation, but C should be rather portable to any machine with a fortran compiler that C can deal with variable names that are longer than six characters. C To run on VAX computers: (1) use the VMS version of the function C DAYTIME, which is given below commented out, instead of the Sun C version; (2) erase the calls to the Sun CPU timing function DTIME in C this file; (3) compile with the G_FLOATING option to increase the C allowed range for floating point numbers. There is no need to change C anything in the file jetsubs3_4.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,BUV,BCO COMMON /MUCOEFS/ AUV,ACO,BUV,BCO 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 CHARACTER CALCTYPE*(1) 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 DOUBLE PRECISION YSCALE,PSCALE COMMON /RENOSCALES/ YSCALE,PSCALE C INTEGER NMJJ,NYSTAR DOUBLE PRECISION MJJGRID(1:30),YSTARGRID(1:30) COMMON /GRIDS/ MJJGRID,YSTARGRID,NMJJ,NYSTAR C DOUBLE PRECISION TEMPARRAY(1:30,1:30) DOUBLE PRECISION XSECTARRAY(1:30,1:30),ERRARRAY(1:30,1:30) COMMON /ARRAYS/ TEMPARRAY,XSECTARRAY,ERRARRAY C DOUBLE PRECISION Y1,YN,DELTAY DOUBLE PRECISION M1,MN,DELTAMJJ C DOUBLE PRECISION MJJSMALL,MJJLARGE C C Results C DOUBLE PRECISION XSECT,ERR,XSECTMID,XSECTOUT C C Smearing and unsmearing parameters C DOUBLE PRECISION LNMJJSMEAR,YSTARSMEAR COMMON /SMEAR/ LNMJJSMEAR,YSTARSMEAR DOUBLE PRECISION TRIAL INTEGER JMAX,KMAX DOUBLE PRECISION CMID(0:9,0:9),COUT(0:9,0:9) C C Number of terms in jet defining function C INTEGER NJETDEF COMMON /DEFJET/ NJETDEF C C Bounds for jet variables C DOUBLE PRECISION MJJMIN,MJJMAX,YSTARMAX,YCMMAX COMMON /JETBOUNDS/ MJJMIN,MJJMAX,YSTARMAX,YCMMAX C C Variables for integration over MJJ C DOUBLE PRECISION DELTAMJJ,DELTAYSTAR,MJJ,YSTAR,FACTOR INTEGER NY C C Integers for DO loops C INTEGER I,J,K,IMJJ,IYSTAR 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 Number of flavors: set in initialization of PARTON C C Lambda_MSbar^(NFL): set in initialization of PARTON C C Max powers of ln(MJJ)^2 and YSTAR^2, respectively, C in polynomial fit for unsmearing: C JMAX = 2 KMAX = 3 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 S2TO2 = .TRUE. SBORN = .TRUE. 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. C C USER INPUT C C Physics parameters: C WRITE(NOUT,*)' Give RTS,R' READ(NIN,*)RTS,R C C Option to compute Born cross section only C WRITE(NOUT,*)' Write f for full calculation, b for Born only' READ(NIN,*)CALCTYPE C C Grid points in MJJ: M1,M2,...,M(NMJJ)=MN C WRITE(NOUT,*)' Give MJJ grid: MJJmin, MJJmax, Number of points' READ(NIN,*) M1,MN,NMJJ IF (NMJJ.GT.30) THEN WRITE(NOUT,*) 'Need NMJJ .LE. 30 please.' STOP ENDIF WRITE(NOUT,*) ' Integrate over what MJJ bin? (MJJSMALL,MJJLARGE)' READ(NIN,*) MJJSMALL,MJJLARGE C C Grid points in YSTAR: Y1,Y2,...,Y(NYSTAR)=YN C WRITE(NOUT,*)' Give YSTAR grid: YSTARmin, YSTARmax,', > ' Number of points' READ(NIN,*) Y1,YN,NYSTAR IF (NYSTAR.GT.30) THEN WRITE(NOUT,*) 'Need NYSTAR .LE. 30 please.' STOP ENDIF C C Cut on YCM C WRITE(NOUT,*)' Give cut on YCM' READ(NIN,*) YCMMAX C 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: cetq3m.ptn)' READ(*,*) PARTONFILE C C C Ratio of renormalization scale MUUV and factorization C scale MUCO to PT: C WRITE(NOUT,*) ' Give scale choices: A_uv, A_co, B_uv, B_co ' READ(NIN,*) AUV,ACO,BUV,BCO 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 C END OF INPUT C C Make y-grid. We will use bins of |y|. C IF(NYSTAR.GT.1)THEN DELTAY = (YN - Y1)/(NYSTAR-1) ELSE DELTAY = 1.0D0 ENDIF DO I = 1,NYSTAR YSTARGRID(I) = Y1 + (I - 1) * DELTAY ENDDO C C Make MJJ-grid equally spaced in MJJ. C IF(NMJJ.GT.1)THEN DELTAMJJ = ( MN - M1 )/(NMJJ-1) ELSE DELTAMJJ = 1.0D0 ENDIF DO I = 1,NMJJ MJJGRID(I) = M1 + (I-1)*DELTAMJJ ENDDO C C Smearing distances in ln(MJJ) and YSTAR: C LNMJJSMEAR = 0.25 * DLOG(0.95 * RTS /MN ) IF(LNMJJSMEAR.GT.0.25) THEN LNMJJSMEAR = 0.25 ENDIF C YSTARSMEAR = 0.2D0 C C Throw away RENO points outside of these limits. C YSTARMAX = YN + 4.0D0 * YSTARSMEAR MJJMIN = M1 * DEXP(- 4.0D0 * LNMJJSMEAR) MJJMAX = MN * DEXP( 4.0D0 * LNMJJSMEAR) 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.5D0 * MJJSMALL C C WRITE THE PARAMETERS C CALL DAYTIME(DATE) WRITE(NOUT,101)DATE 101 FORMAT(// > ,' Program JET version 3.4 ',A,/,50('-')) C WRITE(NOUT,102) 102 FORMAT(' Two jet cross section d sigma / d MJJ d y*') WRITE(NOUT,600) 600 FORMAT('--------------------------------------------------') C CALL PRINTVERSION C IF ((CALCTYPE.EQ.'B').OR.(CALCTYPE.EQ.'b')) THEN SVIRTUAL = .FALSE. SACOLL = .FALSE. SASOFT = .FALSE. SBCOLL = .FALSE. SBSOFT = .FALSE. S1COLL = .FALSE. S1SOFT = .FALSE. S2COLL = .FALSE. S2SOFT = .FALSE. SAFINITE = .FALSE. SBFINITE = .FALSE. S1FINITE = .FALSE. S2FINITE = .FALSE. WRITE(NOUT,*)' ' WRITE(NOUT,*)'BORN Calculation only!' WRITE(NOUT,*)' ' ENDIF 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 WRITE(NOUT,*) ' ' WRITE(NOUT,*) 'Partons:' 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,YCMMAX 613 FORMAT( ' collision energy (GeV): ', F10.1,/, > ' nflavor: ', I10 ,/, > ' lambda (GeV): ', F10.3 ,/, > ' cone radius: ', F10.3,/, > ' cut on y_cm: ', F10.3,/) C WRITE(NOUT,654)AUV,ACO,BUV,BCO 654 FORMAT( ' Choice of scales:',/, > ' MU = A * MJJ / (2 COSH( B * YSTAR) )',/, > ' A_uv:',F10.3,/, > ' A_co:',F10.3,/, > ' B_uv:',F10.3,/, > ' B_co:',F10.3,/) C C WRITE(NOUT,612) NRENO,SEED 612 FORMAT(' RENO will use',I6,' thousand points', > ' with seed',I7,/) WRITE(NOUT,432)LNMJJSMEAR,YSTARSMEAR,JMAX,KMAX 432 FORMAT(' Smearing and unsmearing information:',/, > ' Smearing distance in ln(MJJ):',F10.4,/, > ' Smearing distance in YSTAR:',F10.4,/, > ' Max power of log(MJJ)^2:',I10,/, > ' Max power of YSTAR^2:',I10,/) 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 over TRIAL(MJJ,YSTAR), smeared over ln(MJJ) and C YSTAR. C C Write smeared Cross section/TRIAL(MJJ,YSTAR): C WRITE(NOUT,*)' ' WRITE(NOUT,*)' Smeared cross section / TRIAL(MJJ,YSTAR) ' WRITE(NOUT,*)' ' WRITE(NOUT,*)' MJJ YSTAR Value' C DO IMJJ = 1,NMJJ DO IYSTAR = 1,NYSTAR C WRITE(NOUT,14)MJJGRID(IMJJ),YSTARGRID(IYSTAR), > XSECTARRAY(IMJJ,IYSTAR)*NBGEV2 C 14 FORMAT(F8.1,F7.2,G20.10) C ENDDO WRITE(*,*)' ' ENDDO C C We need to unsmear by fitting the results to a polynomial and C applying an unsmearing transformation to the polynomial. C CALL FITJET(XSECTARRAY,ERRARRAY,MJJGRID,YSTARGRID,NMJJ,NYSTAR, > NRENO,LNMJJSMEAR,YSTARSMEAR,JMAX,KMAX,CMID,COUT) C C The coefficients of the polynomial fit to the smeared function are C now in CMID and the coefficients of the unsmeared polynomial are C now in COUT. C WRITE(NOUT,*)' Coefficients of polynomial fit:' WRITE(NOUT,*)' J K CMID(J,K) COUT(J,K) ' C DO J = 0,JMAX DO K = 0,KMAX WRITE(NOUT,88)J,K,CMID(J,K)*NBGEV2,COUT(J,K)*NBGEV2 88 FORMAT(2I3,2G15.5) ENDDO ENDDO C C---------------------------------------------- C WRITE(NOUT,*)' ' WRITE(NOUT,*)'----------------------------------------', > '---------------------------------------' WRITE(NOUT,*)' ' WRITE(NOUT,*)' RESULTS' WRITE(NOUT,*)' ' WRITE(NOUT,*)'----------------------------------------', > '---------------------------------------' WRITE(NOUT,*)' ' WRITE(NOUT,*)' MJJ YSTAR Smeared Integration ', > ' Unsmearing Corrected Unsmeared ' WRITE(NOUT,*)' XSECT error ', > ' Correction XSECT Fit XSECT ' WRITE(NOUT,*)' ' C C--------- Loop over grid points, printing results for each C DO IMJJ = 1,NMJJ DO IYSTAR = 1,NYSTAR C C XSECT = XSECTARRAY(IMJJ,IYSTAR) ERR = ERRARRAY(IMJJ,IYSTAR) C C Calculate the standard deviation C ERR = DSQRT((ERR - XSECT**2)/(NRENO - 1)) C C Calculate the polynomial functions at the grid point. (The K = 0 C case is done separately in order to avoid 1 = 0**0 if YSTAR = 0.) C XSECTMID = 0.0D0 XSECTOUT = 0.0D0 C DO J = 0,JMAX K = 0 XSECTMID = XSECTMID + DLOG(MJJGRID(IMJJ))**(2*J) * CMID(J,K) XSECTOUT = XSECTOUT + DLOG(MJJGRID(IMJJ))**(2*J) * COUT(J,K) DO K = 1,KMAX XSECTMID = XSECTMID > + DLOG(MJJGRID(IMJJ))**(2*J) > * YSTARGRID(IYSTAR)**(2*K) * CMID(J,K) XSECTOUT = XSECTOUT > + DLOG(MJJGRID(IMJJ))**(2*J) > * YSTARGRID(IYSTAR)**(2*K) * COUT(J,K) ENDDO ENDDO C C Restore factor TRIAL(MJJ,YSTAR) C XSECT = XSECT * TRIAL(MJJGRID(IMJJ),YSTARGRID(IYSTAR)) ERR = ERR * TRIAL(MJJGRID(IMJJ),YSTARGRID(IYSTAR)) XSECTMID = XSECTMID * TRIAL(MJJGRID(IMJJ),YSTARGRID(IYSTAR)) XSECTOUT = XSECTOUT * TRIAL(MJJGRID(IMJJ),YSTARGRID(IYSTAR)) C C Change units from GeV^-2 to nb C XSECT = NBGEV2 * XSECT ERR = NBGEV2 * ERR XSECTMID = NBGEV2 * XSECTMID XSECTOUT = NBGEV2 * XSECTOUT C C Write results C WRITE(NOUT,15)MJJGRID(IMJJ),YSTARGRID(IYSTAR),XSECT,ERR, > XSECTOUT-XSECTMID,XSECT+XSECTOUT-XSECTMID,XSECTOUT 15 FORMAT(F8.1,F7.2,G15.4,G10.2,G10.2,G15.4,G15.4) C ENDDO WRITE(*,*)' ' ENDDO C C Compute the integral over the requested bin in MJJ C WRITE(NOUT,*)' ' WRITE(NOUT,*)'--------------------------------------------' WRITE(NOUT,*)' ' WRITE(NOUT,50) MJJSMALL,MJJLARGE 50 FORMAT(' Cross section integrated from MJJ =',F8.2, > ' to MJJ =',F8.2,' (nb)',/) WRITE(NOUT,51) 51 FORMAT(' YSTAR XSECT XSECT/COSH(2 YSTAR)') C DELTAMJJ = (MJJLARGE - MJJSMALL)/1000.0D0 DELTAYSTAR = 0.1D0 NY = NINT(YN/DELTAYSTAR) C C Loop over desired YSTARs C DO IYSTAR = 0,NY YSTAR = IYSTAR * DELTAYSTAR C C Integrate over MJJ (trapezoidal rule) C XSECT = 0.0D0 DO IMJJ = 0,1000 MJJ = MJJSMALL + IMJJ * DELTAMJJ IF ((IMJJ.EQ.1).OR.(IMJJ.EQ.1000)) THEN FACTOR = 0.5D0 ELSE FACTOR = 1.0D0 ENDIF C C Fit cross section XSECTOUT(MJJ,YSTAR): C XSECTOUT = 0.0D0 DO J = 0,JMAX K = 0 XSECTOUT = XSECTOUT + DLOG(MJJ)**(2*J) * COUT(J,K) DO K = 1,KMAX XSECTOUT = XSECTOUT > + DLOG(MJJ)**(2*J) * YSTAR**(2*K) * COUT(J,K) ENDDO ENDDO XSECTOUT = XSECTOUT * TRIAL(MJJ,YSTAR) C C Done calculating XSECTOUT(MJJ,YSTAR), add to integral: C XSECT = XSECT + XSECTOUT * DELTAMJJ * FACTOR C C End of integration loop. C ENDDO C XSECT = XSECT * NBGEV2 WRITE(NOUT, 96) YSTAR, XSECT, XSECT/DCOSH(2.0D0*YSTAR) 96 FORMAT(F8.2,2G18.6) 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 DOUBLE PRECISION FUNCTION TRIAL(MJJ,YSTAR) DOUBLE PRECISION MJJ,YSTAR C C An approximate version of the cross section, used in the smearing C procedure. C DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C C TRIAL = DEXP( 1.81750D0 > - 4.38459D0 * DLOG(18.0D0*MJJ/RTS) > + 8.89467D0 * DLOG(1.0D0 - MJJ/RTS) ) > * DCOSH(2.0D0*YSTAR) C RETURN END C C*********************************************************************** C*********************************************************************** C C VMS version for DAYTIME C C SUBROUTINE DAYTIME(DT) C CHARACTER DT*(*) C IF (LEN(DT).LE.19) THEN C DT = 'TOO SHORT' C RETURN C ENDIF C CALL TIME(DT(1:8)) C CALL DATE(DT(11:19)) C DT(9:10) = '//' C RETURN C END C C*********************************************************************** C C Sun version for DAYTIME C SUBROUTINE DAYTIME(DATE) CHARACTER*28 DATE C CHARACTER*30 DT CALL FDATE(DT) DATE = DT(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(MJJ,YSTAR,YCM) C DOUBLE PRECISION MJJ,YSTAR,YCM C C DOUBLE PRECISION MJJMIN,MJJMAX,YSTARMAX,YCMMAX COMMON /JETBOUNDS/ MJJMIN,MJJMAX,YSTARMAX,YCMMAX C INBOUNDS = .FALSE. C IF ( (MJJMIN.LT.MJJ).AND.(MJJ.LT.MJJMAX) ) THEN IF ( DABS(YCM).LT.YCMMAX) THEN IF (YSTAR.LT.YSTARMAX) THEN INBOUNDS = .TRUE. ENDIF ENDIF ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE SETMU(MJJ,YSTAR,YCM,MUUV,MUCO) C C Input variables: DOUBLE PRECISION MJJ,YSTAR,YCM C Output variables: DOUBLE PRECISION MUUV,MUCO C C The input arguments are JETVAR1,JETVAR2,JETVAR3, which are the jet C parameters produced by JETDEF and JETDEFBORN. C DOUBLE PRECISION AUV,ACO,BUV,BCO COMMON /MUCOEFS/ AUV,ACO,BUV,BCO C MUUV = AUV * MJJ / (2.0D0 * DCOSH( BUV * YSTAR) ) MUCO = ACO * MJJ / (2.0D0 * DCOSH( BCO * YSTAR) ) C RETURN END C C*********************************************************************** C*********************************************************************** C SUBROUTINE JETDEFBORN(Y1,P2,Y2,PHI2,OK,MJJ,YSTAR,YCM,SVALUE) C C Input variables: DOUBLE PRECISION Y1,P2,Y2,PHI2 C Output variables: LOGICAL OK DOUBLE PRECISION MJJ,YSTAR,YCM,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) C C See JETDEF for defintions for this two jet inclusive cross section. C LOGICAL INBOUNDS C OK = .FALSE. YCM = 0.5D0*(Y1 + Y2) YSTAR = 0.5D0*DABS(Y1 - Y2) MJJ = DSQRT( 2.0D0 * P2**2 * (DCOSH(Y1-Y2) + 1.0D0)) SVALUE = 1.0D0 C IF ( INBOUNDS(MJJ,YSTAR,YCM) ) THEN OK = .TRUE. ENDIF C RETURN END C C*********************************************************************** C SUBROUTINE > JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,OK,MJJ,YSTAR,YCM,SVALUE) C C Input variables: INTEGER CASE DOUBLE PRECISION Y1,P2,Y2,PHI2,P3,Y3,PHI3 C Output variables: LOGICAL OK DOUBLE PRECISION MJJ,YSTAR,YCM,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.) C C The parameter NJETDEF is set in the main program. (Normally it is 3.) C C For two jet inclusive cross section, the variables are C MJJ = the exact invariant mass of the jet-jet system C YSTAR = | y(1) - y(2) |/2 C YCM = ( y(1) + y(2) )/2 C When there are three jets in the event, the two jets with the largest C transverse momenta are used. 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,YJET,CONVERT LOGICAL INBOUNDS DOUBLE PRECISION PI PARAMETER ( PI = 3.141592654 D0) C MJJ = 0.0D0 YSTAR = 0.0D0 YCM = 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 YCM = 0.5D0*(Y1 + Y2) YSTAR = 0.5D0*DABS(Y1 - Y2) MJJ = DSQRT( 2.0D0*P1*P2*(DCOSH(Y1-Y2) - DCOS(PHI12))) IF ( INBOUNDS(MJJ,YSTAR,YCM) ) 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 YJET = (P1*Y1 + P3*Y3)/(P1+P3) YCM = 0.5D0 * (YJET + Y2) YSTAR = 0.5D0 * DABS(YJET - Y2) MJJ = DSQRT( 2.0D0*P1*P2*(DCOSH(Y1-Y2) - DCOS(PHI12)) > + 2.0D0*P2*P3*(DCOSH(Y2-Y3) - DCOS(PHI23)) > + 2.0D0*P3*P1*(DCOSH(Y3-Y1) - DCOS(PHI31)) ) IF ( INBOUNDS(MJJ,YSTAR,YCM) ) 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 YJET = (P2*Y2 + P3*Y3)/(P2+P3) YCM = 0.5D0 * (YJET + Y1) YSTAR = 0.5D0 * DABS(YJET - Y1) MJJ = DSQRT( 2.0D0*P1*P2*(DCOSH(Y1-Y2) - DCOS(PHI12)) > + 2.0D0*P2*P3*(DCOSH(Y2-Y3) - DCOS(PHI23)) > + 2.0D0*P3*P1*(DCOSH(Y3-Y1) - DCOS(PHI31)) ) IF ( INBOUNDS(MJJ,YSTAR,YCM) ) THEN OK = .TRUE. ENDIF ENDIF RETURN C ELSE WRITE(NOUT,*) 'Fatal error in JETDEF' STOP ENDIF END C C C*********************************************************************** C SUBROUTINE BINIT(MJJ,YSTAR,UNUSED,INTEGRAND) C DOUBLE PRECISION MJJ,YSTAR,UNUSED,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 TRIAL(MJJ,YSTAR) and smearing with Gaussians in ln(MJJ) C and in YSTAR. In the case of YSTAR, we define YSTAR to be positive, C and we have grid points at positive YSTAR's only. However, we C wish to smear a smooth cross section that has been continued to C negative YSTAR by defining sigma(-YSTAR) = sigma(YSTAR). For this C reason, we add "artificial) contributions from the corresponding C negative values of YSTAR. IN the case of MJJ, note the Jacobian C factor 1/MJJ, which changes the integration from dMJJ to d ln(MJJ). C INTEGER NMJJ,NYSTAR DOUBLE PRECISION MJJGRID(1:30),YSTARGRID(1:30) COMMON /GRIDS/ MJJGRID,YSTARGRID,NMJJ,NYSTAR C DOUBLE PRECISION LNMJJSMEAR,YSTARSMEAR COMMON /SMEAR/ LNMJJSMEAR,YSTARSMEAR DOUBLE PRECISION TRIAL C DOUBLE PRECISION TEMPARRAY(1:30,1:30) DOUBLE PRECISION XSECTARRAY(1:30,1:30),ERRARRAY(1:30,1:30) COMMON /ARRAYS/ TEMPARRAY,XSECTARRAY,ERRARRAY C DOUBLE PRECISION INTEGRANDMOD INTEGER IMJJ,IYSTAR 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 INTEGRANDMOD = INTEGRAND/TRIAL(MJJ,YSTAR) C DO IMJJ = 1,NMJJ DO IYSTAR = 1,NYSTAR C TEMPARRAY(IMJJ,IYSTAR) = TEMPARRAY(IMJJ,IYSTAR) + INTEGRANDMOD > * DEXP(- DLOG(MJJGRID(IMJJ)/MJJ)**2 /(2.0D0*LNMJJSMEAR**2)) > /(DSQRT(2.0D0*PI)*LNMJJSMEAR * MJJ ) > *( DEXP(- (YSTAR-YSTARGRID(IYSTAR))**2/(2.0D0*YSTARSMEAR**2)) > + DEXP(- (YSTAR+YSTARGRID(IYSTAR))**2/(2.0D0*YSTARSMEAR**2))) > /(DSQRT(2.0D0*PI)*YSTARSMEAR) C ENDDO ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C Fitting Routine for Unsmearing Calculation C*********************************************************************** C*********************************************************************** C SUBROUTINE FITJET(VALARRAY,ERRARRAY,MJJGRID,YSTARGRID,NMJJ,NYSTAR, > NRENO,LNMJJSMEAR,YSTARSMEAR,JMAX,KMAX,CTILDE,C) DOUBLE PRECISION VALARRAY(1:30,1:30),ERRARRAY(1:30,1:30) DOUBLE PRECISION MJJGRID(1:30),YSTARGRID(1:30) DOUBLE PRECISION LNMJJSMEAR,YSTARSMEAR INTEGER NMJJ,NYSTAR,NRENO,JMAX,KMAX DOUBLE PRECISION C(0:9,0:9),CTILDE(0:9,0:9) C C Fits VALARRAY(IMJJ,IYSTAR) to a polynomial in log(MJJ) and YSTAR C using a modified version of P. Nason's subroutine FITVAL. C C The grid points in MJJ and YSTAR are in MJJGRID and YSTARGRID, C with sizes NMJJ and NYSTAR. NRENO is the number of RENO iterations, C used in computing the error. C C Then applies an 'unsmoothing' transformation to the polynomial. C The resulting polynomial has the form C C Sum C(J,K) * log(MJJ)^(2J) * YSTAR^(2K) C J,K C C The indices J and K run from 0 to JMAX and KMAX, respectively. C The parameters LNMJJSMEAR,YSTARSMEAR are the smearing lengths in C log(MJJ) and YSTAR, respecitively: the gaussian smearing function is C exp{- (x - xbar)^2 /(2 xSMEAR^2) } up to normalization. C C The parameters JMAX and KMAX are input. (They should not be large.) C The coefficients C(J,K) are the output. C C INTEGER MAXPAR,MAXPT PARAMETER (MAXPAR=100) PARAMETER (MAXPT=900) C DOUBLE PRECISION CTILDEJK(MAXPAR) DOUBLE PRECISION VAL(MAXPT),ERR(MAXPT) DOUBLE PRECISION FUN(MAXPAR,MAXPT) C INTEGER J,K,JK,M,N INTEGER IMJJ,IYSTAR,LPY,NPOINTS,NPAR DOUBLE PRECISION UNSMEAR,DELTA DOUBLE PRECISION CHISQ C C C We first map the given VALARRAY with two indices into a one C dimensional array VAL for the purpose of using FITVAL. The mapping C is L = IYSTAR + NYSTAR*(IMJJ - 1). We also compute a one dimensional C array ERR from ERRARRAY. Recall that ERRARRAY contains the average C of the squares of the Monte Carlo data; we compute in ERR the C corresponding standard deviations. C C We also compute the monomials to be fit at the grid points. The C monomials, log(MJJ)^(2J) * YSTAR^(2K), are labelled by a one C dimensionalindex JK, defined by JK = K + (KMAX + 1)*J + 1. C (The K = 0 case is calculated separately in order to C avoid 0**0 in case YSTAR = 0.) C DO IMJJ = 1,NMJJ DO IYSTAR = 1,NYSTAR LPY = IYSTAR + NYSTAR * (IMJJ - 1) VAL(LPY) = VALARRAY(IMJJ,IYSTAR) ERR(LPY) = (ERRARRAY(IMJJ,IYSTAR)-VALARRAY(IMJJ,IYSTAR)**2) > /(NRENO - 1) ERR(LPY) = DSQRT(ERR(LPY)) C-- DO J = 0,JMAX JK = (KMAX + 1)*J + 1 FUN(JK,LPY) = DLOG(MJJGRID(IMJJ))**(2*J) DO K = 1,KMAX JK = K + (KMAX + 1)*J + 1 FUN(JK,LPY) = DLOG(MJJGRID(IMJJ))**(2*J) > * YSTARGRID(IYSTAR)**(2*K) ENDDO ENDDO C-- ENDDO ENDDO C NPOINTS = NMJJ*NYSTAR NPAR = (JMAX + 1)*(KMAX + 1) C C call FITVAL(FUN,VAL,DEV,NPT,NPAR,PARAMS,CHISQ) C CALL FITVAL(FUN,VAL,ERR,NPOINTS,NPAR,CTILDEJK,CHISQ) C C Calculate C (and also CTILDE) from CTILDEJK C DO M = 0,JMAX DO N = 0,KMAX C C(M,N) = 0.0D0 CTILDE(M,N) = 0.0D0 C DO J = M,JMAX DO K = N,KMAX JK = K + (KMAX + 1)*J + 1 C C(M,N) = C(M,N) > + UNSMEAR(M,J,LNMJJSMEAR)*UNSMEAR(N,K,YSTARSMEAR)*CTILDEJK(JK) C CTILDE(M,N) = CTILDE(M,N) > + DELTA(M,J)*DELTA(N,K)*CTILDEJK(JK) ENDDO ENDDO ENDDO ENDDO C C Finished C RETURN END C C***********************************************************************