C23456789012345678901234567890123456789012345678901234567890123456789012 PROGRAM PARTONLOOK C C This small program is a "front end" for the fortran function PARTON, C enabling the user to inspect the parton distribution interpolated C from a given data file. C C DES, 30 January 1996 C DOUBLE PRECISION PARTON, X, MU, RESULT INTEGER NPARTON 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 INTEGER I,IMAX DOUBLE PRECISION XMIN,XMID,XMAX,B,C C NIN = 5 NOUT = 6 C NFLAVOR = 5.0D0 NFL = 5 IMAX = 10 XMIN = 0.0001D0 XMAX = 0.99 XMID = 0.1 B = DLOG(XMID/XMIN)/IMAX C = ( XMAX - XMID ) / DBLE(IMAX) C WRITE(*,*)'Check parton distributions' WRITE(*,*)'Give data file: cteq2m.ptn, cteq3m.ptn, etc' READ(*,*) PARTONFILE WRITE(*,*)'SCALE?' READ(*,*) MU WRITE(*,*)'PARTON NUMBER?' READ(*,*) NPARTON IF ((NPARTON.GT.6).OR.(NPARTON.LT.-6) ) THEN STOP ENDIF DO I = 0,IMAX X = XMIN * EXP(B*I) RESULT = X**2 * PARTON(NPARTON,X,MU) WRITE(*,2)X, RESULT ENDDO DO I = 1,IMAX X = XMID + C * I RESULT = X**2 * PARTON(NPARTON,X,MU) WRITE(*,2)X, RESULT ENDDO 2 FORMAT(F10.5,F12.8) 3 FORMAT(F8.6,G20.4) END C C23456789012345678901234567890123456789012345678901234567890123456789012 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 Small revisions on input and output for compatability with alpha C workstations. 18 May 1999, DES 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,801) VARNAME 801 FORMAT(A80) ENDDO IF ((INDEX(VARNAME,'PARTON').EQ.0).AND. > (INDEX(VARNAME,'Parton').EQ.0).AND. > (INDEX(VARNAME,'parton').EQ.0)) GOTO 199 READ(28,802) VARNAME,DUMMY,VERSION 802 FORMAT(A10,A3,I13) 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,803) VARNAME,DUMMY,LAMBDA 803 FORMAT(A10,A3,G20.10) IF ((INDEX(VARNAME,'LAMBDA').EQ.0).AND. > (INDEX(VARNAME,'Lambda').EQ.0).AND. > (INDEX(VARNAME,'lambda').EQ.0)) GOTO 199 READ(28,804) VARNAME,DUMMY,LAMBDAFLAVORS 804 FORMAT(A10,A3,I10) 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,805) VARNAME,DUMMY,XMIN 805 FORMAT(A10,A3,G20.10) IF ((INDEX(VARNAME,'XMIN').EQ.0).AND. > (INDEX(VARNAME,'Xmin').EQ.0).AND. > (INDEX(VARNAME,'xmin').EQ.0)) GOTO 199 READ(28,806) VARNAME,DUMMY,XMAX 806 FORMAT(A10,A3,G20.10) IF ((INDEX(VARNAME,'XMAX').EQ.0).AND. > (INDEX(VARNAME,'Xmax').EQ.0).AND. > (INDEX(VARNAME,'xmax').EQ.0)) GOTO 199 READ(28,807) VARNAME,DUMMY,NA 807 FORMAT(A10,A3,I10) 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,808) VARNAME,DUMMY,MUMIN 808 FORMAT(A10,A3,G20.10) IF ((INDEX(VARNAME,'MUMIN').EQ.0).AND. > (INDEX(VARNAME,'MuMin').EQ.0).AND. > (INDEX(VARNAME,'mumin').EQ.0)) GOTO 199 READ(28,809) VARNAME,DUMMY,MUMAX 809 FORMAT(A10,A3,G20.10) IF ((INDEX(VARNAME,'MUMAX').EQ.0).AND. > (INDEX(VARNAME,'MuMax').EQ.0).AND. > (INDEX(VARNAME,'mumax').EQ.0)) GOTO 199 READ(28,810) VARNAME,DUMMY,NB 810 FORMAT(A10,A3,I10) 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***********************************************************************