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