MODULE GS09ST C GS module which stores the dPDF grids and sample T and X values. integer :: NX INTEGER :: NT INTEGER, PARAMETER :: NP = 8 INTEGER, PARAMETER :: NXMAX = 100 INTEGER, PARAMETER :: NTMAX = 100 DOUBLE PRECISION :: XX(NXMAX) DOUBLE PRECISION :: TT(NTMAX) DOUBLE PRECISION,DIMENSION(:,:,:,:,:),ALLOCATABLE, TARGET :: FF DOUBLE PRECISION,DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: JJBMAT DOUBLE PRECISION, POINTER, DIMENSION(:,:,:) :: PDFARY => NULL() LOGICAL :: FSTTIM = .TRUE. DOUBLE PRECISION, PARAMETER :: EPS=1d-8 DOUBLE PRECISION :: XMAX INTEGER NMCOR END MODULE GS09ST C---------------------------------------------------------------------- C Fortran interpolation code for GS dPDFs. Returns value of dPDF C with parton indices given by F1 and F2 at scale Q (in GeV) and momentum C arguments X1 and X2. F1 and F2 must be supplied as integers using C the PDG mapping between flavours and integers - i.e. C -6, -5, -4, -3, -2, -1,0,1,2,3,4,5,6 C = tbar,bbar,cbar,sbar,ubar,dbar,g,d,u,s,c,b,t. C Comments to Jo Gaunt . C---------------------------------------------------------------------- SUBROUTINE GS09(X1,X2,Q,F1,F2,DPDF) USE GS09ST IMPLICIT NONE INTEGER DMYINT DOUBLE PRECISION DMYDUB DOUBLE PRECISION, INTENT(IN) :: X1, X2, Q DOUBLE PRECISION, INTENT(OUT) :: DPDF INTEGER, INTENT(IN) :: F1,F2 INTEGER REALF1, REALF2 DOUBLE PRECISION T LOGICAL JJSWIT EXTERNAL OURNOT INTEGER OURNOT DOUBLE PRECISION QSQ INTEGER JJINT DOUBLE PRECISION X1I, X2I INTEGER F1I, F2I X1I = X1 X2I = X2 F1I = F1 F2I = F2 C Adjust parameters if outside grid. IF(X1I .LT. 1D-6) THEN X1I = 1E-6 END IF IF(X2I .LT. 1D-6) THEN X2I = 1E-6 END IF QSQ = Q*Q IF(QSQ .GT. 1D9) THEN QSQ = 1D9 END IF IF(QSQ .LT. 1D0) THEN QSQ = 1D0 END IF T = LOG(QSQ) C Return zero if user requests dPDF too near kinematic C boundary. IF(X1I+X2I.GT.1d0-1d-8) THEN DPDF = 0d0 RETURN END IF JJSWIT = .FALSE. C Initialise if first time call. IF(FSTTIM) THEN CALL GSINIT() FSTTIM = .FALSE. END IF C Special cases - ssbar, ccbar, bbbar distributions C are NOT equal to ss etc. For these, set PDFARY C to point to the relevant jjbar grid. IF(F1I+F2I.EQ.0.AND.ABS(F1I)>2.AND.ABS(F1I)<6) THEN IF(F2I .GT. 0) THEN DMYDUB = X1I X1I=X2I X2I=DMYDUB END IF JJINT = ABS(F1I)-2 PDFARY => JJBMAT(JJINT, 1:NX,1:NX,1:NT) JJSWIT = .TRUE. C Otherwise can use symmetry of sea for charm, bottom and C strange (nb we have set s = sbar here unlike MSTW 2008). C So if user requests ucbar for example, can print value of C uc as the two are equal. ELSE IF(F1I<-2) F1I = ABS(F1I) IF(F2I<-2) F2I = ABS(F2I) END IF C If user requests any dPDF with a top (or antitop) index, return C zero. IF(F1I.EQ.6 .OR. F2I.EQ.6) THEN DPDF = 0d0 RETURN END IF C We only have grid files for parton combinations with REALF1 integers) - this is to avoid redundancy. Here C we check whether the requested dPDF has this property - if not, swap C F1 <-> F2 and X1 <-> X2 simultaneously. IF(.NOT.JJSWIT) THEN F1I = OURNOT(F1I) F2I = OURNOT(F2I) IF(F2I .LT. F1I) THEN DMYINT = F1I F1I = F2I F2I = DMYINT DMYDUB = X1I X1I=X2I X2I=DMYDUB END IF C Set PDFARY to point to appropriate FF grid, unless it's already pointing C to a jjbar grid. PDFARY => FF(F1I,F2I,1:NX,1:NX,1:NT) END IF C Do the interpolation! CALL GSINTR(X1I,X2I,T,DPDF) NULLIFY(PDFARY) END SUBROUTINE GS09 FUNCTION OURNOT(PDGNOT) C Translates between our internal mapping between parton flavours and integers, C and that used by the PDG. IMPLICIT NONE INTEGER OURNOT, PDGNOT IF(PDGNOT.GT.2) OURNOT = PDGNOT IF(PDGNOT.EQ.1) OURNOT = 2 IF(PDGNOT.EQ.2) OURNOT = 1 IF(PDGNOT.EQ.-1) OURNOT = 7 IF(PDGNOT.EQ.-2) OURNOT = 6 IF(PDGNOT.EQ.0) OURNOT = 8 END FUNCTION OURNOT SUBROUTINE GSINTR(X1,X2,T,DPDF) C GS interpolation subroutine USE GS09ST IMPLICIT NONE DOUBLE PRECISION X1, X2, Q, DPDF INTEGER F1, F2 DOUBLE PRECISION FFLCL(2,2,2) EXTERNAL CORINT, LININT, GSLOCX, SIMPLN DOUBLE PRECISION LININT, CORINT, SIMPLN INTEGER GSLOCX INTEGER I, J, K INTEGER N1, N2, M DOUBLE PRECISION FRCX1, FRCX2, DELT, FRCT DOUBLE PRECISION FFX(2,2) DOUBLE PRECISION FFX2(2), FFQ(2) DOUBLE PRECISION T DOUBLE PRECISION ZERO, FRCX1LO, FRCX1HI DOUBLE PRECISION FRCX2LO, FRCX2HI ZERO = 0d0 M = GSLOCX(TT,NT,T) DELT = TT(M+1)-TT(M) FRCT = (T-TT(M))/DELT C If (x1,x2) point is right in the corner of the triangle (i.e C has an x1 or x2 value higher than the largest sampling value), C do a special corner interpolation (basically triangulation). IF(X1.GE.XMAX) THEN DO I = 1,2 FFQ(I) = CORINT(PDFARY(NX,1:NMCOR,M+I-1),X1,X2) END DO DPDF = LININT(FRCT, FFQ(1), FFQ(2)) RETURN ELSEIF(X2.GE.XMAX) THEN DO I = 1,2 FFQ(I) = CORINT(PDFARY(1:NMCOR,NX,M+I-1),X2,X1) END DO DPDF = LININT(FRCT, FFQ(1), FFQ(2)) RETURN END IF N1 = GSLOCX(XX,NX,X1) IF(N1.LT.1) N1 = 1 N2 = GSLOCX(XX,NX,X2) IF(N2.LT.1) N2 = 1 C Construct cube containing nearest 8 values. DO I = 1,2 DO J = 1,2 DO K = 1,2 FFLCL(I,J,K) = PDFARY(N1+I-1,N2+J-1,M+K-1) END DO END DO END DO C Interpolate in t. DO I = 1,2 DO J = 1,2 FFX(I,J) = LININT(FRCT, FFLCL(I,J,1),FFLCL(I,J,2)) END DO END DO C Check if point is right next to kinematic boundary - if so, C do special boundary interpolation. IF(XX(N1+1)+XX(N2+1) .GT. 1d0+EPS) THEN IF(XX(N1+1)+XX(N2) .LE. 1d0+eps) THEN FRCX1LO = (X1-XX(N1))/(XX(N1+1)-XX(N1)) FFX2(1) = LININT(FRCX1LO,FFX(1,1),FFX(2,1)) IF(X1 .GT. 1-XX(N2+1)) THEN FRCX2 = (X2-XX(N2))/(1-X1-XX(N2)) DPDF = LININT(FRCX2,FFX2(1),ZERO) RETURN ELSE FRCX1HI = (X1-XX(N1))/(1-XX(N2+1)-XX(N1)) FFX2(2) = LININT(FRCX1HI, FFX(1,2), ZERO) FRCX2 = (X2-XX(N2))/(XX(N2+1)-XX(N2)) DPDF = LININT(FRCX2, FFX2(1),FFX2(2)) RETURN END IF ELSEIF(XX(N2+1)+XX(N1) .LE. 1d0+eps) THEN FRCX2LO = (X2-XX(N2))/(XX(N2+1)-XX(N2)) FFX2(1) = LININT(FRCX2LO,FFX(1,1),FFX(1,2)) IF(X2 .GT. 1-XX(N1+1)) THEN FRCX1 = (X1-XX(N1))/(1-X2-XX(N1)) DPDF = LININT(FRCX1, FFX2(1),ZERO) RETURN ELSE FRCX2HI = (X2-XX(N2))/(1-XX(N1+1)-XX(N2)) FFX2(2) = LININT(FRCX2HI, FFX(2,1), ZERO) FRCX1 = (X1-XX(N1))/(XX(N1+1)-XX(N1)) DPDF = LININT(FRCX1, FFX2(1),FFX2(2)) RETURN END IF END IF FRCX1 = (X1-XX(N1))/(1-XX(N2)-XX(N1)) FRCX2 = (X2-XX(N2))/(1-XX(N1)-XX(N2)) DPDF = FFX(1,1)*SIMPLN(FRCX1,FRCX2) RETURN END IF C Do interpolation in x1 IF(XX(N1+1) .LT. 0.1d0) THEN FRCX1 = (LOG(X1)-LOG(XX(N1)))/(LOG(XX(N1+1))-LOG(XX(N1))) DO J = 1,2 IF(FFX(1,J).GT.0d0 .AND. FFX(2,J) .GT. 0d0) THEN FFX2(J) = LININT(FRCX1, LOG(FFX(1,J)), LOG(FFX(2,J))) FFX2(J) = EXP(FFX2(J)) ELSE FFX2(J) = LININT(FRCX1, FFX(1,J), FFX(2,J)) END IF END DO ELSE FRCX1 = (X1-XX(N1))/(XX(N1+1)-XX(N1)) DO J = 1,2 FFX2(J) = LININT(FRCX1, FFX(1,J), FFX(2,J)) END DO END IF C Do interpolation in x2 IF(XX(N2+1) .LT. 0.1d0) THEN FRCX2 = (LOG(X2)-LOG(XX(N2)))/(LOG(XX(N2+1))-LOG(XX(N2))) IF(FFX2(1) .GT. 0d0 .AND. FFX2(2) .GT. 0d0) THEN DPDF = LININT(FRCX2, LOG(FFX2(1)), LOG(FFX2(2))) DPDF = EXP(DPDF) ELSE DPDF = LININT(FRCX2, FFX2(1), FFX2(2)) END IF ELSE FRCX2 = (X2-XX(N2))/(XX(N2+1)-XX(N2)) DPDF = LININT(FRCX2, FFX2(1), FFX2(2)) END IF END SUBROUTINE GSINTR SUBROUTINE GSINIT() C GS initialisation subroutine. USE GS09ST IMPLICIT NONE INTEGER I, J, K, PART1, PART2 DOUBLE PRECISION X1, X2 INTEGER GSLOCX EXTERNAL GSLOCX INTEGER IO DOUBLE PRECISION DMYDUB OPEN(FILE = "dpdfgrids/XVALUES.txt", UNIT = 37) OPEN(FILE = "dpdfgrids/TVALUES.txt", UNIT = 38) C Put the sampling x values in the XX array NX = 1 DO READ(37,*,IOSTAT=IO) XX(NX) IF(IO .LT. 0) EXIT NX = NX + 1 END DO NX = NX - 1 C Put the sampling T values in the TT array. NT = 1 DO READ(38,*,IOSTAT=IO) TT(NT) IF(IO .LT. 0) EXIT NT = NT + 1 END DO NT = NT - 1 CLOSE(37) CLOSE(38) ALLOCATE(FF(8,8,NX,NX,NT)) ALLOCATE(JJBMAT(3,NX,NX,NT)) C Read in grid files to FF and JJBMAT arrays. OPEN(FILE = "dpdfgrids/gs2009lo.u.dat", UNIT = 21) OPEN(FILE = "dpdfgrids/gs2009lo.d.dat", UNIT = 22) OPEN(FILE = "dpdfgrids/gs2009lo.s.dat", UNIT = 23) OPEN(FILE = "dpdfgrids/gs2009lo.c.dat", UNIT = 24) OPEN(FILE = "dpdfgrids/gs2009lo.b.dat", UNIT = 25) OPEN(FILE = "dpdfgrids/gs2009lo.ub.dat", UNIT = 26) OPEN(FILE = "dpdfgrids/gs2009lo.db.dat", UNIT = 27) OPEN(FILE = "dpdfgrids/gs2009lo.g.dat", UNIT = 28) OPEN(FILE = "dpdfgrids/gs2009lo.jjb.dat", UNIT = 29) DO K = 1,NT DO I = 1, NX X1 = XX(I) DO J = 1,NX X2 = XX(J) DO PART1 = 1, NP IF(X1+X2-1 .GT. 0d0) THEN DO PART2 = PART1, NP FF(PART1, PART2, I, J, K) =0d0 END DO ELSE READ(20+PART1,*) & (FF(PART1, PART2, I, J, K), & PART2 = PART1,NP) END IF END DO IF(X1+X2-1 .GT. 0d0) THEN JJBMAT(1,I,J,K) = 0d0 JJBMAT(2,I,J,K) = 0d0 JJBMAT(3,I,J,K) = 0d0 ELSE READ(29,*) (JJBMAT(PART1,I,J,K), PART1 = 1,3) END IF END DO END DO END DO C Check to make sure number of entries in each grid file is as C expected. DO K = 1,NT DO PART1 = 1, NP READ(20+PART1,*,IOSTAT = IO ) DMYDUB IF(IO .EQ. 0) THEN WRITE(*,*) "GS Error: Wrong number & of entries in grid files." STOP END IF END DO READ(29,*, IOSTAT = IO) DMYDUB IF(IO .EQ. 0) THEN WRITE(*,*) "GS Error: Wrong number & of entries in grid file." STOP END IF END DO C Largest sampling X value - required for determining C whether one should do a corner interpolation or not. XMAX = XX(NX) C Number of samples there are in the x2 direction which C have x1 equal to the maximum sampling value. Also C used in corner interpolation. NMCOR = GSLOCX(XX, NX, 1-XMAX+EPS) DO I = 1,9 CLOSE(20+I) END DO END SUBROUTINE GSINIT FUNCTION SIMPLN(XS,YS) C For a plane which passes through (0,0,1), (1,0,0), (0,1,0), finds C the z value corresponding to x = xs, y = ys. Used when interpolating C near the kinematic boundary. IMPLICIT NONE DOUBLE PRECISION XS,YS, SIMPLN SIMPLN = 1d0 - XS - YS END FUNCTION SIMPLN SUBROUTINE FULPLN(X,Y,Z,XS,YS,ZOUT) C Finds the z value of a point which has x = xs, y = ys, and which C also lies on the plane passing through the vectors x(3), y(3), and C z(3). Used in the corner interpolation. IMPLICIT NONE DOUBLE PRECISION X(3), Y(3), Z(3) DOUBLE PRECISION XS, YS, ZOUT DOUBLE PRECISION A, B, C, D A = Y(1)*(Z(2)-Z(3))+Y(2)*(Z(3)-Z(1))+Y(3)*(Z(1)-Z(2)) B = Z(1)*(X(2)-X(3))+Z(2)*(X(3)-X(1))+Z(3)*(X(1)-X(2)) C = X(1)*(Y(2)-Y(3))+X(2)*(Y(3)-Y(1))+X(3)*(Y(1)-Y(2)) D = -X(1)*(Y(2)*Z(3)-Y(3)*Z(2))-X(2)*(Y(3)*Z(1)-Y(1)*Z(3)) & -X(3)*(Y(1)*Z(2)-Y(2)*Z(1)) ZOUT = (-D-A*XS-B*YS)/C END SUBROUTINE FULPLN FUNCTION LININT(XS, Z0, ZF) C One-dimensional linear interpolation. Interpolates linearly C between the points (0,Z0) and (1,ZF) to give the value of C Z at X = XS. IMPLICIT NONE DOUBLE PRECISION XS, Z0, ZF, LININT LININT = ZF*XS + (1d0-XS)*Z0 END FUNCTION LININT INTEGER FUNCTION FDPLPT(XS, YS) C Used during the corner interpolation. During this interpolation, the C corner of the kinematic region is divided into triangles by drawing lines radiating from C the corner to the sampled points with maximum x1/x2, and then joining up all C these points. This function essentially finds the triangle containing the point. USE GS09ST IMPLICIT NONE DOUBLE PRECISION YCOMP, XS, YS INTEGER I DO I = 1,NMCOR YCOMP = (1d0-XS)/(1-XMAX)*XX(I) IF (YS .LT. YCOMP) THEN FDPLPT = I-1 EXIT END IF END DO IF(YS .GT. YCOMP) THEN FDPLPT = -10 END IF END FUNCTION FDPLPT INTEGER FUNCTION GSLOCX(XX,NX,X) c-- returns an integer j such that x lies inbetween xx(j) and xx(j+1). c-- nx is the length of the array with xx(nx) the highest element. IMPLICIT NONE INTEGER NX,JL,JU,JM DOUBLE PRECISION X,XX(NX) IF(X.EQ.XX(1)) THEN GSLOCX=1 RETURN ENDIF IF(X.GE.XX(NX)) THEN GSLOCX=NX-1 RETURN ENDIF JU=NX+1 JL=0 1 IF((JU-JL).LE.1) GO TO 2 JM=(JU+JL)/2 IF(X.GE.XX(JM)) THEN JL=JM ELSE JU=JM ENDIF GO TO 1 2 GSLOCX=JL RETURN END FUNCTION CORINT(FCOR,X,Y) C Corner interpolation routine. In this interpolation, the C corner of the kinematic region is divided into triangles by drawing lines radiating from C the corner to the sampled points with maximum x1/x2, and then joining up all C these points. The points of the triangle that contains the requested (x1,x2) are C found with the help of FDPLPT. A plane is drawn through this triangle via FULPLN, and the Z value C of the plane at (x1,x2) is used as the interpolated value. USE GS09ST IMPLICIT NONE DOUBLE PRECISION CORINT DOUBLE PRECISION FCOR(NMCOR) INTEGER IL, FDPLPT DOUBLE PRECISION XVALS(3), YVALS(3), ZVALS(3) DOUBLE PRECISION X,Y IL = FDPLPT(X,Y) IF (IL == -10) THEN CORINT = 0d0 RETURN END IF XVALS = (/ XMAX, XMAX, 1d0 /) YVALS = (/ XX(IL), XX(IL+1), 0d0 /) ZVALS = (/ FCOR(IL), FCOR(IL+1), 0d0 /) CALL FULPLN(XVALS, YVALS, ZVALS, X,Y, CORINT) END FUNCTION CORINT