C23456789012345678901234567890123456789012345678901234567890123456789012 C C Version 1.2 C 7 March 1997 C C Routines for calculating parton distributions from a table C This set of routines allows the use of more than one parton C distribution set at once. C C There are three parts C 1) PROGRAM PARTONLOOK C a small routine to demonstrate how the other two parts work. C 2) SUBROUTINE OPENPARTON(PARTONFILE,SET) C called to set up the needed tables in memory. C 3) FUNCTION PARTON(NPARTON,XIN,SCALE,SET) C called to compute the desired parton distribution function C at a particular value of x and scale. C C This is version 1.2 of this package. C Version 1.1 had a bug in OPENPARTON. C C23456789012345678901234567890123456789012345678901234567890123456789012 C 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 & PA, 30 January 1996 C Updated: 11 April 1996 C DOUBLE PRECISION PARTON, X, MU(10), RESULT INTEGER NPARTON(10) C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL INTEGER NIN,NOUT,NWRT COMMON /IOUNIT/ NIN,NOUT,NWRT C CHARACTER*100 PARTONFILE C INTEGER I,IMAX,J DOUBLE PRECISION XMIN,XMID,XMAX,B,C INTEGER SETNUMB 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(*,*)'Number of parton distributions?' READ(*,*) SETNUMB DO J = 1, SETNUMB WRITE(*,*) WRITE(*,*)'PARTON DISTRIBUTION ',J,' :' WRITE(*,*)'Check parton distributions' WRITE(*,*)'Give data file: cteq2m.ptn, cteq3m.ptn, etc' READ(*,*) PARTONFILE WRITE(*,*)'SCALE?' READ(*,*) MU(J) WRITE(*,*)'PARTON NUMBER?' READ(*,*) NPARTON(J) IF ((NPARTON(J).GT.6).OR.(NPARTON(J).LT.-6) ) THEN STOP ENDIF CALL OPENPARTON(PARTONFILE,J) ENDDO WRITE(*,*)'************************' DO J = 1,SETNUMB WRITE(*,*) WRITE(*,*)'PARTON DISTRIBUTION ',J,' :' WRITE(*,*)' X (X**2)*Parton Distribution' DO I = 0,IMAX X = XMIN * EXP(B*I) RESULT = X**2 * PARTON(NPARTON(J),X,MU(J),J) WRITE(*,2)X, RESULT ENDDO DO I = 1,IMAX X = XMID + C * I RESULT = X**2 * PARTON(NPARTON(J),X,MU(J),J) WRITE(*,2)X, RESULT ENDDO ENDDO 2 FORMAT(F10.5,F12.8) 3 FORMAT(F8.6,G20.4) END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE OPENPARTON(PARTONFILE,SET) CHARACTER*100 PARTONFILE INTEGER SET C C Setup subroutine to support function PARTON(NPARTON,XIN,SCALE,SET), C a function to interpolate parton distributions from a data C file. This function allows the use of more than one parton C distribution set at once. The set is specified by the integer SET. C This setup subroutine sets up the data tables for each set. C C This version corrects a bug in version 1.1. C Old, wrong, code: C LMUMIN(SET) = DLOG(MUMIN) C LMUMAX = DLOG(MUMAX) C New, right, code: C LMUMIN(SET) = DLOG(MUMIN(SET)) C LMUMAX = DLOG(MUMAX(SET)) C C Version 1.2 C P. Anandam and D. E. Soper C 7 March 1997 C C Previous version, 1.1 C P. Anandam and D. E. Soper C 11 April 1996 C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL INTEGER NIN,NOUT,NWRT COMMON /IOUNIT/ NIN,NOUT,NWRT C INTEGER NAMAX,NBMAX,NSET PARAMETER (NAMAX=64) PARAMETER (NBMAX=32) PARAMETER (NSET=6) C DOUBLE PRECISION XMIN(NSET),XMAX(NSET) DOUBLE PRECISION LXMIN(NSET),LXMAX,DELTALX(NSET) DOUBLE PRECISION MUMIN(NSET),MUMAX(NSET) DOUBLE PRECISION LMUMIN(NSET),LMUMAX,DELTALMU(NSET) INTEGER N,NA(NSET),NB(NSET),IA,IB COMMON /DATASET/ XMIN, XMAX, LXMIN, DELTALX, > MUMIN, MUMAX, LMUMIN, DELTALMU, NA, NB C DOUBLE PRECISION H(NSET,-2:6,NAMAX,NBMAX) DOUBLE PRECISION HA(NSET,-2:6,NAMAX,NBMAX) DOUBLE PRECISION HB(NSET,-2:6,NAMAX,NBMAX) DOUBLE PRECISION HAB(NSET,-2:6,NAMAX,NBMAX) COMMON /HVALUES/ H,HA,HB,HAB DOUBLE PRECISION HBL,HBR,HABL,HABR C CHARACTER*100 HEADER CHARACTER*100 VARNAME CHARACTER*100 DUMMY DOUBLE PRECISION VERSION INTEGER LAMBDAFLAVORS C C 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 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(SET) 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(SET) 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(SET) 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(SET).GT.NAMAX) THEN WRITE(NOUT,*)' N_XPOINTS too big for declared array size' ENDIF READ(28,*) VARNAME,DUMMY,MUMIN(SET) 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(SET) 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(SET) 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(SET).GT.NBMAX) THEN WRITE(NOUT,*)' N_MUPOINTS too big for declared array size' ENDIF DO 20 IA=1,NA(SET) DO 20 IB=1,NB(SET) READ(28,*)H(SET,-2,IA,IB),H(SET,-1,IA,IB),H(SET,0,IA,IB), > H(SET,1,IA,IB), H(SET,2,IA,IB), H(SET,3,IA,IB), > H(SET,4,IA,IB), H(SET,5,IA,IB), H(SET,6,IA,IB) 20 CONTINUE C Calculate lattice info. C LXMIN(SET) = DLOG(XMIN(SET)/(1.0D0 - XMIN(SET))) LXMAX = DLOG(XMAX(SET)/(1.0D0 - XMAX(SET))) DELTALX(SET) = (LXMAX - LXMIN(SET))/(NA(SET)-1) LMUMIN(SET) = DLOG(MUMIN(SET)) LMUMAX = DLOG(MUMAX(SET)) DELTALMU(SET) = (LMUMAX - LMUMIN(SET))/(NB(SET)-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(SET)-1 DO IB = 1,NB(SET) HA(SET,N,IA,IB) =0.5D0*(H(SET,N,IA+1,IB)-H(SET,N,IA-1,IB)) ENDDO ENDDO DO IB = 1,NB(SET) HA(SET,N,1,IB) = H(SET,N,2,IB) - H(SET,N,1,IB) HA(SET,N,NA(SET),IB) = H(SET,N,NA(SET),IB) > - H(SET,N,NA(SET)-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(SET) DO IB = 2,NB(SET)-1 HBL = H(SET,N,IA,IB) - H(SET,N,IA,IB-1) HBR = H(SET,N,IA,IB+1) - H(SET,N,IA,IB) IF (DABS(HBL).LT.DABS(HBR)) THEN HB(SET,N,IA,IB) = HBL ELSE HB(SET,N,IA,IB) = HBR ENDIF ENDDO ENDDO DO IA = 1,NA(SET) HB(SET,N,IA,1) = H(SET,N,IA,2) - H(SET,N,IA,1) HB(SET,N,IA,NB(SET)) = H(SET,N,IA,NB(SET)) > - H(SET,N,IA,NB(SET)-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(SET)-1 DO IB = 2,NB(SET)-1 HABL = 0.5D0*( H(SET,N,IA+1,IB) - H(SET,N,IA-1,IB) > - H(SET,N,IA+1,IB-1) + H(SET,N,IA-1,IB-1) ) HABR = 0.5D0*( H(SET,N,IA+1,IB+1) - H(SET,N,IA-1,IB+1) > - H(SET,N,IA+1,IB) + H(SET,N,IA-1,IB) ) IF (DABS(HABL).LT.DABS(HABR)) THEN HAB(SET,N,IA,IB) = HABL ELSE HAB(SET,N,IA,IB) = HABR ENDIF ENDDO ENDDO DO IA = 2,NA(SET)-1 HAB(SET,N,IA,1) = 0.5D0*( H(SET,N,IA+1,2) - H(SET,N,IA-1,2) > - H(SET,N,IA+1,1) + H(SET,N,IA-1,1) ) HAB(SET,N,IA,NB(SET)) = 0.5D0*( H(SET,N,IA+1,NB(SET)) > -H(SET,N,IA-1,NB(SET)) - H(SET,N,IA+1,NB(SET)-1) + > H(SET,N,IA-1,NB(SET)-1)) ENDDO DO IB = 2,NB(SET)-1 HABL = H(SET,N,2,IB) - H(SET,N,1,IB) > - H(SET,N,2,IB-1) + H(SET,N,1,IB-1) HABR = H(SET,N,2,IB+1) - H(SET,N,1,IB+1) > - H(SET,N,2,IB) + H(SET,N,1,IB) IF (DABS(HABL).LT.DABS(HABR)) THEN HAB(SET,N,1,IB) = HABL ELSE HAB(SET,N,1,IB) = HABR ENDIF HABL = H(SET,N,NA(SET),IB) - H(SET,N,NA(SET)-1,IB) > - H(SET,N,NA(SET),IB-1) + H(SET,N,NA(SET)-1,IB-1) HABR = H(SET,N,NA(SET),IB+1) - H(SET,N,NA(SET)-1,IB+1) > - H(SET,N,NA(SET),IB) + H(SET,N,NA(SET)-1,IB) IF (DABS(HABL).LT.DABS(HABR)) THEN HAB(SET,N,NA(SET),IB) = HABL ELSE HAB(SET,N,NA(SET),IB) = HABR ENDIF ENDDO HAB(SET,N,1,1) = H(SET,N,2,2) - H(SET,N,1,2) > - H(SET,N,2,1) + H(SET,N,1,1) HAB(SET,N,NA(SET),1) = H(SET,N,NA(SET),2) > - H(SET,N,NA(SET)-1,2) - H(SET,N,NA(SET),1) + > H(SET,N,NA(SET)-1,1) HAB(SET,N,1,NB(SET)) = H(SET,N,2,NB(SET)) > - H(SET,N,1,NB(SET)) - H(SET,N,2,NB(SET)-1) + > H(SET,N,1,NB(SET)-1) HAB(SET,N,NA(SET),NB(SET)) = H(SET,N,NA(SET),NB(SET)) - > H(SET,N,NA(SET)-1,NB(SET)) - H(SET,N,NA(SET),NB(SET)-1) + > H(SET,N,NA(SET)-1,NB(SET)-1) C C End DO N=1,9 C ENDDO C C 101 FORMAT(A80) C C WRITE(NOUT,601) PARTONFILE 601 FORMAT('Read file ',A20,' with header:') WRITE(NOUT,602) HEADER 602 FORMAT(A70) WRITE(NOUT,603) VERSION 603 FORMAT('The parton data file version is: ',F10) WRITE(NOUT,604)LAMBDA,NFL 604 FORMAT('These partons have LAMBDA =',F6.3, > ' for NFL =',I3) CLOSE(28) GOTO 755 899 WRITE(NOUT,*) 'Error opening parton data file!' STOP 199 WRITE(NOUT,*) 'Wrong format of parton data file!' STOP 755 CONTINUE END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C DOUBLE PRECISION FUNCTION PARTON(NPARTON,XIN,SCALE,SET) INTEGER NPARTON, SET 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.2 C P. Anandam and D. E. Soper C 7 March 1997 C C The function PARTON(...) here is identical with the previous version. C Only the setup subroutine OPENPARTON is changed between versions C 1.1 and 1.2. C C Version 1.1 C P. Anandam and D. E. Soper C 11 April 1996 C Based on function PARTON(NPARTON,XIN,SCALE) for a single C parton distribution set, 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 C INTEGER NAMAX,NBMAX,NSET PARAMETER (NAMAX=64) PARAMETER (NBMAX=32) PARAMETER (NSET=6) C DOUBLE PRECISION X,XMIN(NSET),XMAX(NSET) DOUBLE PRECISION LX,LXMIN(NSET),DELTALX(NSET) DOUBLE PRECISION MU,MUMIN(NSET),MUMAX(NSET) DOUBLE PRECISION LMUMIN(NSET),LMU,DELTALMU(NSET) INTEGER N,NA(NSET),NB(NSET),IA,IB DOUBLE PRECISION NORMALIZEDLX,NORMALIZEDLMU COMMON /DATASET/ XMIN, XMAX, LXMIN, DELTALX, > MUMIN, MUMAX, LMUMIN, DELTALMU, NA, NB C DOUBLE PRECISION H(NSET,-2:6,NAMAX,NBMAX) DOUBLE PRECISION HA(NSET,-2:6,NAMAX,NBMAX) DOUBLE PRECISION HB(NSET,-2:6,NAMAX,NBMAX) DOUBLE PRECISION HAB(NSET,-2:6,NAMAX,NBMAX) COMMON /HVALUES/ H,HA,HB,HAB 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 OLDSET C C Note that we save the variables! C SAVE C DATA OLDSET/0/ C C C ------------------- Initialization -------------------------- C IF(OLDSET.EQ.SET)GOTO 10 OLDSET=SET C C Reset the X and MU values C X=-1D0 MU=-1D0 C 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(SET) ).OR.(MU.GT.MUMAX(SET) )) THEN WRITE(NOUT,*)'PARTON called out of range, MU =',MU STOP ENDIF LMU = DLOG(MU) NORMALIZEDLMU = (LMU - LMUMIN(SET))/DELTALMU(SET) + 1.0D0 IB = INT(NORMALIZEDLMU) IF((IB.LT.1).OR.(IB.GT.(NB(SET)-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(SET))/DELTALX(SET) + 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(SET)) THEN EXTRAPOLATE = .TRUE. IA = NA(SET) 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(SET,N,IA, IB ) * B0 NET = NET + H(SET,N,IA, IB+1) * B1 NET = NET + HA(SET,N,IA ,IB ) * A * B0 NET = NET + HA(SET,N,IA ,IB+1) * A * B1 NET = NET + HB(SET,N,IA ,IB ) * B0B NET = NET + HB(SET,N,IA ,IB+1) * B1B NET = NET + HAB(SET,N,IA ,IB ) * A * B0B NET = NET + HAB(SET,N,IA ,IB+1) * A * B1B C ELSE C NET = H(SET,N,IA, IB ) * A0 * B0 NET = NET + H(SET,N,IA+1,IB ) * A1 * B0 NET = NET + H(SET,N,IA, IB+1) * A0 * B1 NET = NET + H(SET,N,IA+1,IB+1) * A1 * B1 C NET = NET + HA(SET,N,IA ,IB ) * A0A * B0 NET = NET + HA(SET,N,IA+1,IB ) * A1A * B0 NET = NET + HA(SET,N,IA ,IB+1) * A0A * B1 NET = NET + HA(SET,N,IA+1,IB+1) * A1A * B1 C NET = NET + HB(SET,N,IA ,IB ) * A0 * B0B NET = NET + HB(SET,N,IA+1,IB ) * A1 * B0B NET = NET + HB(SET,N,IA ,IB+1) * A0 * B1B NET = NET + HB(SET,N,IA+1,IB+1) * A1 * B1B C NET = NET + HAB(SET,N,IA ,IB ) * A0A * B0B NET = NET + HAB(SET,N,IA+1,IB ) * A1A * B0B NET = NET + HAB(SET,N,IA ,IB+1) * A0A * B1B NET = NET + HAB(SET,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 C END C C*********************************************************************** C***********************************************************************