C23456789012345678901234567890123456789012345678901234567890123456789012
C
C JET version 3.4
C
C Version to calculate the one jet inclusive cross section
C d sigma/ d ET
C
C S.D. Ellis, Z. Kunszt, D.E. Soper
C
C 18 March 1997
C
C Version 3.0, 23 February 1991
C Version 3.01, 14 March 1991
C Version 3.1, 28 February 1992
C ET distribution integrated over yJ, 10 Mar 1992
C Version 3_3, 15 August 1995
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 Additions to include Rsep, 12 December 1996
C Modification to TRIAL, 18 March 1997
C Allows better performance with RTS other than
C 1800 GeV. Results for RTS = 1800 GeV should
C be unchanged.
C 18 August 1997 version
C Version to calculate d sigma/ d ET1 dy1 dy2
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 ET1 dy1 dy2 >
C
C Here there are two jets plus anything and jets 1 and 2 are the
C two jets in the event with the highest ET. (ET1,y1) are the transverse
C energy and rapidity of jet 1; y2 is the rapidity of jet 2. We
C average the cross section over intervals in y1 and y2 specified
C by the user. The cross section is given on a grid of points in ET1.
C
C Other master programs can calculate other quantities of interest
C when linked with jetsubs3_4.f.
C
C Results from this program were first presented in
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 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
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 PJ, the program calculates the
C 'SmearedXSECT,' S(PJ), given in terms of the true cross section,
C s(pj), by
C
C S(PJ) ~ f(PJ) * \int dln(pj) h[ln(PJ)-ln(pj)] s(pj)/f(pj).
C
C where h is a gaussian smearing function and f (called TRIAL(PJ)
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(PJ), which is a polynomial in
C ln(PJ), to S(PJ)/f(PJ):
C
C S(PJ) ~ f(PJ) * Poly(PJ).
C
C (3) We calculate the unsmeared polynomial poly(pj) such that if
C poly is smeared as above the result is (exactly) Poly:
C
C Poly(PJ) = \int dln(p) h[ln(PJ)-ln(pj)] poly(pj).
C
C (4) The 'Unsmeared Fit XSECT' is
C
C f(PJ) * poly(PJ).
C
C (5) The 'Corrected XSECT' is
C
C S(PJ) + f(PJ) * [ Poly(PJ) - poly(PJ) ].
C
C Each of the three answers has its advantages. To the extent that
C the f(PJ) times a polynomial of order, say, 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 pj 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(PJ) 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(PJ)
C that is a very good fit to S(PJ). Also, it is possible to change
C the number of parameters used in the fit. The default is JMAX = 2
C but with small statistical errors one could raise JMAX.
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 MUUVOVERPJ,MUCOOVERPJ
COMMON /MUCOEFS/ MUUVOVERPJ,MUCOOVERPJ
DOUBLE PRECISION RSEP,RSEPOVERR
COMMON /CONE/ RSEP
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 NPJ,NYJ
DOUBLE PRECISION PJGRID(1:30),YJGRID(1:30)
COMMON /GRIDS/ PJGRID,YJGRID,NPJ,NYJ
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 P1,PN,DELTAPJ
C
DOUBLE PRECISION PJMIN,PJMAX,YJ1MIN,YJ1MAX,YJ2MIN,YJ2MAX
COMMON /JETBOUNDS/ PJMIN,PJMAX,YJ1MIN,YJ1MAX,YJ2MIN,YJ2MAX
C
C Results
C
DOUBLE PRECISION XSECT,ERR,XSECTMID,XSECTOUT
C
C Smearing and unsmearing parameters
C
DOUBLE PRECISION LNPJSMEAR,YJSMEAR
COMMON /SMEAR/ LNPJSMEAR,YJSMEAR
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 Integers for DO loops
C
INTEGER I,J,K,IPJ,IYJ
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(PJ)^2 and YJ^2, respectively, in polynomial fit for
C unsmearing:
C
JMAX = 4
KMAX = 0
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,Rsep/R'
READ(NIN,*)RTS,R,RSEPOVERR
RSEP = RSEPOVERR * 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 PJ: P1,P2,...,P(NPJ)=PN
C
WRITE(NOUT,*)' Give PJ grid: PJmin, PJmax, Number of points'
READ(NIN,*) P1,PN,NPJ
IF (NPJ.GT.30) THEN
WRITE(NOUT,*) 'Need NPJ .LE. 30 please.'
STOP
ENDIF
C
C YJ cuts
C
WRITE(NOUT,*)' Give ABS(YJ) cuts for jet 1: YJ1min, YJ1max'
READ(NIN,*) YJ1MIN,YJ1MAX
WRITE(NOUT,*)' Give ABS(YJ) cuts for jet 2: YJ2min, YJ2max'
READ(NIN,*) YJ2MIN,YJ2MAX
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 PJ:
C
WRITE(NOUT,*) ' Give scale choices: mu_uv / ET1 and mu_co / ET1'
READ(NIN,*) MUUVOVERPJ,MUCOOVERPJ
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
C YJ-grid (not really used, dummy arguments for fit routine)
C
NYJ = 1
YJGRID(1) = 0.0D0
C
C Make PJ-grid equally spaced in PJ.
C
IF(NPJ.GT.1)THEN
DELTAPJ = ( PN - P1 )/(NPJ-1)
ELSE
DELTAPJ = 1.0D0
ENDIF
DO I = 1,NPJ
PJGRID(I) = P1 + (I-1)*DELTAPJ
ENDDO
C
C Scales to tell RENO where to put rapidities and transverse
C momenta for partons 1 and 2.
C
YSCALE = 0.6D0 * YJ2MAX
PSCALE = 0.5D0 * DSQRT(P1*PN)
C
C Smearing distances in ln(PJ) and YJ (YJSMEAR is just a dummy value):
C
IF (PN.LT.(DEXP(-YJ2MAX)*RTS/3.0D0)) THEN
LNPJSMEAR = 0.20D0
ELSE
LNPJSMEAR = 0.10D0
ENDIF
C
YJSMEAR = 1.0D0
C
C Throw away RENO points outside of these limits.
C
PJMIN = P1 * DEXP(- 4.0D0 * LNPJSMEAR)
PJMAX = PN * DEXP( 4.0D0 * LNPJSMEAR)
C
C--------------------------
C
C WRITE THE PARAMETERS
C
CALL DAYTIME(DATE)
WRITE(NOUT,101)DATE
101 FORMAT(//
> ,' Program JET version 3.4 (18 August 1997) ',/,
> ' Run on ',A,/,50('-'),/,
> ,' This version includes Rsep.',/)
C
WRITE(NOUT,102)
102 FORMAT(' Jet ET Distribution',/,
> ' d sigma / dET1 dy1 dy2, averaged over y1 and y2')
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,RSEPOVERR,MUUVOVERPJ,MUCOOVERPJ
613 FORMAT( ' collision energy (GeV): ', F10.1,/,
> ' nflavor: ', I10 ,/,
> ' lambda (GeV): ', F10.3 ,/,
> ' cone radius R: ', F10.3,/,
> ' Rsep/R (Snowmass = 2): ', F10.3,/,
> ' mu_ultraviolet/PJ: ', F10.3,/,
> ' mu_collinear/PJ: ', F10.3,/)
C
WRITE(NOUT,433)YJ1MIN,YJ1MAX,YJ2MIN,YJ2MAX
433 FORMAT( ' Cuts on absolute value of YJ:',/,
> ' YJ1min:',F10.4,/,
> ' YJ1max:',F10.4,/,
> ' YJ2min:',F10.4,/,
> ' YJ2max:',F10.4,/)
C
WRITE(NOUT,612) NRENO,SEED
612 FORMAT(' RENO will use',I6,' thousand points',
> ' with seed',I7,/)
C
WRITE(NOUT,432)LNPJSMEAR,JMAX
432 FORMAT(' Smearing and unsmearing information:',/,
> ' Smearing distance in ln(PJ):',F10.4,/,
> ' Max power of log(PJ)^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(PJ,YJ), smeared over ln(PJ) and
C YJ.
C
C Write smeared Cross section/TRIAL(PJ,YJ):
C
WRITE(NOUT,*)' '
WRITE(NOUT,*)
> ' Smeared cross section / TRIAL(PJ)'
WRITE(NOUT,*)' '
WRITE(NOUT,*)' PJ Value'
C
DO IPJ = 1,NPJ
C
WRITE(NOUT,14)PJGRID(IPJ), XSECTARRAY(IPJ,1)
> *NBGEV2 /( 4.0D0*(YJ1MAX-YJ1MIN)*(YJ2MAX-YJ2MIN) )
C
14 FORMAT(F8.1,G18.10)
C
ENDDO
WRITE(NOUT,*)' '
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,PJGRID,YJGRID,NPJ,NYJ,NRENO,
> LNPJSMEAR,YJSMEAR,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,601)
601 FORMAT('----------------------------------------',
> '----------------------------------------')
WRITE(NOUT,*)' '
WRITE(NOUT,602)
602 FORMAT(' Cross section <d sigma / d ET1 d Y1 d Y2>'
> ' (nb/GeV) ',/, ' (averaged over Y1 and Y2) ')
WRITE(NOUT,*)' '
WRITE(NOUT,601)
WRITE(NOUT,*)' '
WRITE(NOUT,603)
603 FORMAT(' ET1 Smeared Integration',
> ' Unsmearing Corrected Unsmeared ')
WRITE(NOUT,604)
604 FORMAT(' XSECT error ',
> ' Correction XSECT Fit XSECT ')
WRITE(NOUT,*)' '
C
C--------- Loop over grid points, printing results for each
C
IYJ = 1
DO IPJ = 1,NPJ
C
C
XSECT = XSECTARRAY(IPJ,IYJ)
ERR = ERRARRAY(IPJ,IYJ)
C
C Calculate the standard deviation
C
ERR = DSQRT((ERR - XSECT**2)/(NRENO - 1))
C
C Calculate the polynomial functions at the grid point.
C
XSECTMID = 0.0D0
XSECTOUT = 0.0D0
C
K = 0
DO J = 0,JMAX
XSECTMID = XSECTMID
> + DLOG(PJGRID(IPJ))**(2*J) * CMID(J,K)
XSECTOUT = XSECTOUT
> + DLOG(PJGRID(IPJ))**(2*J) * COUT(J,K)
ENDDO
C
C Restore factor TRIAL(PJ)
C
XSECT = XSECT * TRIAL(PJGRID(IPJ))
ERR = ERR * TRIAL(PJGRID(IPJ))
XSECTMID = XSECTMID * TRIAL(PJGRID(IPJ))
XSECTOUT = XSECTOUT * TRIAL(PJGRID(IPJ))
C
C Factor for average over YJs instead of integral over YJs:
C
XSECT = XSECT /( 4.0D0*(YJ1MAX-YJ1MIN)*(YJ2MAX-YJ2MIN) )
ERR = ERR /( 4.0D0*(YJ1MAX-YJ1MIN)*(YJ2MAX-YJ2MIN) )
XSECTMID = XSECTMID /( 4.0D0*(YJ1MAX-YJ1MIN)*(YJ2MAX-YJ2MIN) )
XSECTOUT = XSECTOUT /( 4.0D0*(YJ1MAX-YJ1MIN)*(YJ2MAX-YJ2MIN) )
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)PJGRID(IPJ),XSECT,ERR,
> XSECTOUT-XSECTMID,XSECT+XSECTOUT-XSECTMID,XSECTOUT
15 FORMAT(F8.2,G15.4,G12.2,G12.2,G15.4,G15.4)
C
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(PJ)
DOUBLE PRECISION PJ
C
C An approximate version of the cross section, used in the smearing
C procedure. The first part of this is a careful fit for the
C ordinary one jet inclusive cross section. This is modified
C to work reasonably well for d sigma/ dET1 dy1 dy2. Evidently
C the modification is rather ad hoc.
C
DOUBLE PRECISION PJMIN,PJMAX,YJ1MIN,YJ1MAX,YJ2MIN,YJ2MAX
COMMON /JETBOUNDS/ PJMIN,PJMAX,YJ1MIN,YJ1MAX,YJ2MIN,YJ2MAX
DOUBLE PRECISION RTS,R
COMMON /PHYSICS/ RTS,R
C
DOUBLE PRECISION LL,LH,LHY
C
LL = DLOG(PJ*18.0D0/RTS)
LH = DLOG(1.0D0 - 2.0D0 * PJ/RTS)
LHY = - 3.0 * DLOG((1.0D0 - PJ/RTS * DEXP(YJ2MIN))**2 + 0.07D0)
C
TRIAL = DEXP( - 3.31589D0
> - 7.57548D0 * LL
> - 22.7385D0 * LH
> - 0.875663D0 * LL**2
> - 1.77124D0 * LH**2
> + 11.3259D0 * LL * LH )
C
TRIAL = TRIAL * DEXP(- 0.2D0 * YJ2MIN * (LL + 0.7D0)**2 - LHY)
C
RETURN
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***********************************************************************
C23456789012345678901234567890123456789012345678901234567890123456789012
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(PJ,YJ1,YJ2)
C
DOUBLE PRECISION PJ,YJ1,YJ2
C
DOUBLE PRECISION PJMIN,PJMAX,YJ1MIN,YJ1MAX,YJ2MIN,YJ2MAX
COMMON /JETBOUNDS/ PJMIN,PJMAX,YJ1MIN,YJ1MAX,YJ2MIN,YJ2MAX
C
DOUBLE PRECISION YJ1ABS,YJ2ABS
C
YJ1ABS = DABS(YJ1)
YJ2ABS = DABS(YJ2)
C
INBOUNDS = .FALSE.
C
IF ( (PJMIN.LT.PJ).AND.(PJ.LT.PJMAX) ) THEN
IF ( (YJ1MIN.LT.YJ1ABS).AND.(YJ1ABS.LT.YJ1MAX) ) THEN
IF ( (YJ2MIN.LT.YJ2ABS).AND.(YJ2ABS.LT.YJ2MAX) ) THEN
INBOUNDS = .TRUE.
ENDIF
ENDIF
ENDIF
C
RETURN
END
C
C***********************************************************************
C
SUBROUTINE SETMU(PJ,YJ1,YJ2,MUUV,MUCO)
C
C Input variables:
DOUBLE PRECISION PJ,YJ1,YJ2
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 MUUVOVERPJ,MUCOOVERPJ
COMMON /MUCOEFS/ MUUVOVERPJ,MUCOOVERPJ
C
MUUV = PJ * MUUVOVERPJ
MUCO = PJ * MUCOOVERPJ
C
RETURN
END
C
C***********************************************************************
C
SUBROUTINE JETDEFBORN(Y1,P2,Y2,PHI2,OK,PJ,YJ1,YJ2,SVALUE)
C
C Input variables:
DOUBLE PRECISION Y1,P2,Y2,PHI2
C Output variables:
LOGICAL OK
DOUBLE PRECISION PJ,YJ1,YJ2,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 = 2 to account
C for the breaking of 1<-->2 symmetry.)
C
C We identify jet 1 with parton 2 and jet 2 with parton 1.
C
LOGICAL INBOUNDS
C
OK = .FALSE.
PJ = P2
YJ1 = Y2
YJ2 = Y1
SVALUE = 2.0D0
C
IF ( INBOUNDS(PJ,YJ1,YJ2) ) THEN
OK = .TRUE.
ENDIF
C
RETURN
END
C
C***********************************************************************
C
SUBROUTINE
> JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3,OK,PJ,YJ1,YJ2,SVALUE)
C
C Input variables:
INTEGER CASE
DOUBLE PRECISION Y1,P2,Y2,PHI2,P3,Y3,PHI3
C Output variables:
LOGICAL OK
DOUBLE PRECISION PJ,YJ1,YJ2,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 = 2 to account
C for the breaking of 1<-->2 symmetry.)
C
C The parameter NJETDEF is set in the main program. (Normally it is 3.)
C
C
DOUBLE PRECISION RTS,R
COMMON /PHYSICS/ RTS,R
DOUBLE PRECISION RSEP
COMMON /CONE/ RSEP
INTEGER NIN,NOUT,NWRT
COMMON /IOUNIT/ NIN,NOUT,NWRT
DOUBLE PRECISION PHI32
LOGICAL INBOUNDS
DOUBLE PRECISION CONVERT
DOUBLE PRECISION RLIMITSQ,DELTASQ
C
PJ = 0.0D0
YJ1 = 0.0D0
YJ2 = 0.0D0
OK = .FALSE.
C
SVALUE = 2.0D0
C
PHI32 = CONVERT(PHI3 - PHI2)
RLIMITSQ = MIN( RSEP**2, ( (P2+P3)/P2 * R )**2 )
DELTASQ = (Y3-Y2)**2 + PHI32**2
C
C CASE 1: 2 and 3 are the jet
C
IF (CASE.EQ.1) THEN
IF( DELTASQ.LT.RLIMITSQ ) THEN
PJ = P2 + P3
YJ1 = ( P2*Y2 + P3*Y3 )/PJ
YJ2 = Y1
IF ( INBOUNDS(PJ,YJ1,YJ2) ) THEN
OK = .TRUE.
ENDIF
ENDIF
RETURN
C
C CASE 2: 2 is the jet
C
ELSEIF (CASE.EQ.2) THEN
IF( DELTASQ.GE.RLIMITSQ ) THEN
PJ = P2
YJ1 = Y2
YJ2 = Y1
IF ( INBOUNDS(PJ,YJ1,YJ2) ) THEN
OK = .TRUE.
ENDIF
ENDIF
RETURN
C
C CASE 3: 3 is the jet
C WE DO NOT ALLOW THIS, SINCE IT HAS THE SMALLEST ET!
C
ELSEIF (CASE.EQ.3) THEN
RETURN
C
ELSE
WRITE(NOUT,*) 'Fatal error in JETDEF'
STOP
ENDIF
END
C
C***********************************************************************
C
SUBROUTINE BINIT(PJ,UNUSED1,UNUSED2,INTEGRAND)
C
DOUBLE PRECISION PJ,UNUSED1,UNUSED2,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(PJ) and smearing with Gaussians in ln(PJ).
C Note the Jacobian factor 1/PJ, which changes the integration
C from dPJ to d ln(PJ).
C
C
INTEGER NPJ,NYJ
DOUBLE PRECISION PJGRID(1:30),YJGRID(1:30)
COMMON /GRIDS/ PJGRID,YJGRID,NPJ,NYJ
C
DOUBLE PRECISION LNPJSMEAR,YJSMEAR
COMMON /SMEAR/ LNPJSMEAR,YJSMEAR
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 IPJ,IYJ
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(PJ)
C
IYJ = 1
DO IPJ = 1,NPJ
C
TEMPARRAY(IPJ,IYJ) = TEMPARRAY(IPJ,IYJ) + INTEGRANDMOD
> * DEXP(- DLOG(PJGRID(IPJ)/PJ)**2 /(2.0D0*LNPJSMEAR**2))
> /(DSQRT(2.0D0*PI)*LNPJSMEAR * PJ )
C
ENDDO
C
RETURN
END
C
C***********************************************************************
C***********************************************************************
C Fitting Routine for Unsmearing Calculation
C***********************************************************************
C***********************************************************************
C
SUBROUTINE FITJET(VALARRAY,ERRARRAY,PJGRID,YJGRID,NPJ,NYJ,NRENO,
> LNPJSMEAR,YJSMEAR,JMAX,KMAX,CTILDE,C)
DOUBLE PRECISION VALARRAY(1:30,1:30),ERRARRAY(1:30,1:30)
DOUBLE PRECISION PJGRID(1:30),YJGRID(1:30)
DOUBLE PRECISION LNPJSMEAR,YJSMEAR
INTEGER NPJ,NYJ,NRENO,JMAX,KMAX
DOUBLE PRECISION C(0:9,0:9),CTILDE(0:9,0:9)
C
C Fits VALARRAY(IPJ,IYJ) to a polynomial in log(PJ) and YJ using
C a modified version of P. Nason's subroutine FITVAL.
C
C The grid points in PJ and YJ are in PJGRID and YJGRID,
C with sizes NPJ and NYJ. 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(PJ)^(2J) * YJ^(2K)
C J,K
C
C The indices J and K run from 0 to JMAX and KMAX, respectively.
C The parameters LNPJSMEAR,YJSMEAR are the smearing lengths in
C log(PJ) and YJ, 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 IPJ,IYJ,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 = IYJ + NYJ*(IPJ - 1). We also compute a one dimensional array
C ERR from ERRARRAY. Recall that ERRARRAY contains the average of the
C 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(PJ)^(2J) * YJ^(2K), are labelled by a one dimensional
C index JK, defined by JK = K + (KMAX + 1)*J + 1. (The K = 0 case
C is calculated separately in order to avoid 0**0 in case YJ = 0.)
C
DO IPJ = 1,NPJ
DO IYJ = 1,NYJ
LPY = IYJ + NYJ * (IPJ - 1)
VAL(LPY) = VALARRAY(IPJ,IYJ)
ERR(LPY) = (ERRARRAY(IPJ,IYJ)-VALARRAY(IPJ,IYJ)**2)
> /(NRENO - 1)
ERR(LPY) = DSQRT(ERR(LPY))
C--
DO J = 0,JMAX
JK = (KMAX + 1)*J + 1
FUN(JK,LPY) = DLOG(PJGRID(IPJ))**(2*J)
DO K = 1,KMAX
JK = K + (KMAX + 1)*J + 1
FUN(JK,LPY) = DLOG(PJGRID(IPJ))**(2*J) * YJGRID(IYJ)**(2*K)
ENDDO
ENDDO
C--
ENDDO
ENDDO
C
NPOINTS = NPJ*NYJ
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,LNPJSMEAR)*UNSMEAR(N,K,YJSMEAR)*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***********************************************************************