C--------------------------------------------------------------------------- C F2 Structure Function Calculation C with C Variable Flavor Number Scheme or Fixed Flavor Number Scheme C C Ref[1]: "IMPROVED PARTON DISTRIBUTIONS FROM GLOBAL ANALYSIS OF RECENT DEEP C INELASTIC SCATTERING AND INCLUSIVE JET DATA" C By: H.L. Lai, J. Huston, S. Kuhlmann, F. Olness, J. Owens, D. Soper C W.K. Tung, H. Weerts C Phys. Rev. D55, 1280 (1997) C This paper explains the properties of the traditional zero-mass parton C scheme CTEQ parton distributions listed as C Iset = 1 - 10 C in the table below. C C Ref[2]: "CHARM PRODUCTION AND PARTON DISTRIBUTIONS" C By: H.L. Lai and W.K. Tung C MSU-HEP-61222, CTEQ-622, e-Print Archive: hep-ph/9701256 C to appear in Z. Phys. C The e-Print hep-ph/9701256 describes the C CTEQ4HQ and CTEQ4F3,4 distributions. C The published version (Z.Phys.) contains a Note Added in Proof which C introduces a refined set CTEQ4HQ1 which is preferred over CTEQ4HQ. The C latter is kept only for historical reason. These distributions C correspond to Iset = 11 - 14 in the table below. C -------------------------------------------------------------------- C For charm production calculation including charm mass effects, the C scheme used to calculate the hard scattering cross-section must C match that used to generate the parton distributions. This program C is provided mainly to allow a user to C calculate the F2 structure function according to the generalized C factorization scheme formulated by C Aivazis, Collins, Olness, and Tung (ACOT). C C The ACOT scheme is formulated in a series of papers: C Nucl.Phys.B278:934,1986 C Nucl.Phys.B308:813,1988 C Phys.Rev.Lett.65:2339-2342,1990 C Phys.Rev.D50:3102-3118,1994. C C Since the CTEQ4HQ(1) distributions are the only ones generated in this C scheme, the current package is set up to use only the CTEQ parton C distribution sets listed below. The programs and data files for C calculating CTEQ4 parton distributions can be downloaded from C http://www.phys.psu.edu/~cteq/ C C For comparison purposes, we have also included options to calculate C F2 in the conventional zero-mass parton formalism (used by all C popular parton distributions so far, from EHLQ to MRS and CTEQ until C the latest sets) as well as in the C fixed-flavor-number (FFN) scheme (used by GRV and most "NLO" calculations C by the StonyBrook-Leiden collaboration -- J. Smith, W. van Neerven etal.). C Cf. Ref.[2] above for a brief review of these comparisons. C ------------------------------------------------------------------------- C Before calculating F2, it is necessary to initialize by the call: C Call SetF2(Isch, Iset) C where C Isch: The choice of Scheme C 0 --> Variable Flavor Number Scheme (zero-mass parton scheme) C 1 --> Variable Flavor Number Scheme (ACOT scheme) C 2 --> Fixed Flavor Number Scheme C Iset: The choice of PDF's C This package contains 13 sets of CTEQ4 PDF's. Details are: C --------------------------------------------------------------------------- C Iset PDF Description Alpha_s(Mz) Lam4 Lam5 Table_File C --------------------------------------------------------------------------- C Ref[1] C 1 CTEQ4M Standard MSbar scheme 0.116 298 202 cteq4m.tbl C 2 CTEQ4D Standard DIS scheme 0.116 298 202 cteq4d.tbl C 3 CTEQ4L Leading Order 0.132 236 181 cteq4l.tbl C 4 CTEQ4A1 Alpha_s series 0.110 215 140 cteq4a1.tbl C 5 CTEQ4A2 Alpha_s series 0.113 254 169 cteq4a2.tbl C 6 CTEQ4A3 ( same as CTEQ4M ) C 7 CTEQ4A4 Alpha_s series 0.119 346 239 cteq4a4.tbl C 8 CTEQ4A5 Alpha_s series 0.122 401 282 cteq4a5.tbl C 9 CTEQ4HJ High Jet 0.116 303 206 cteq4hj.tbl C 10 CTEQ4LQ Low Q0 0.114 261 174 cteq4lq.tbl C --------------------------------------------------------------------------- C Ref[2] C 11 CTEQ4HQ Heavy Quark 0.116 298 202 cteq4hq.tbl C 12 CTEQ4HQ1 Heavy Quark:Q0=1,Mc=1.3 0.116 298 202 cteq4hq1.tbl C (Improved version of CTEQ4HQ, recommended) C 13 CTEQ4F3 Nf=3 FixedFlavorNumber 0.106 (Lam3=385) cteq4f3.tbl C 14 CTEQ4F4 Nf=4 FixedFlavorNumber 0.111 292 XXX cteq4f4.tbl C --------------------------------------------------------------------------- C The valid combinations of (Isch, Iset) are: C (0, Iset) Iset = 1 - 10 (VFN: zero-mass parton scheme) C (1, Iset) Iset = 11, 12(preferred) (VFN: ACOT scheme) C (2, Iset) Iset = 13, 14 (fixed-flavor number scheme) C C The function F2(Mode, X, Q) C returns the value of F2 struction function of the proton in the C desired mode, where C Mode: The available mode for F2 calculation C 1 --> F2 Total C 2 --> F2charm C 3 --> F2bottom C X : Bjorken X C Q : Scale (GeV) C C If you have detailed questions, or if you find problems/bugs using C this package, direct inquires to C Hung-Liang Lai(Lai_H@pa.msu.edu) or Wu-Ki Tung(Tung@pa.msu.edu). C C Below is a simple test program as a sample. C--------------------------------------------------------------------------- Program F2Test IMPLICIT DOUBLE PRECISION (A-H, O-Z) Write(*,*)'Input Isch, Iset, Mode, X, Q' read(*,*) Isch, Iset, Mode, X, Q Call SetF2(Isch, Iset) ff2 = F2(Mode, X, Q) If (Mode.eq.1) Then Print*,'Total F2=', ff2 Elseif (Mode.eq.2) Then Print*,'F2charm=', ff2 Elseif (Mode.eq.3) Then Print*,'F2bottom=', ff2 Endif End C ========================================================================= Function F2(Mode, X, Q) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Parameter(Iset4F4=14) Common /Ischeme/ Isch, Iset Data Pmass /.938D0/ Q2 = Q * Q W = SQRT (-Q2 + Pmass + Q2 / X) If (W .lt. 2D0*Bmass(4)) then NEFF=3 Elseif (W .lt. 2D0*Bmass(5)) then NEFF=4 Elseif (W .lt. 2D0*Bmass(6)) then NEFF=5 Else NEFF=6 Endif If (Isch .eq. 0) then If (Mode.eq.1) then F2total=STlight(Neff, X, Q) Elseif(Mode.eq.2) then If (Neff.ge.4) then F2c=STlight(4, X, Q)-STlight(3, X, Q) Else F2c=0D0 Endif Elseif(Mode.eq.3) then If (Neff.ge.5) then F2b=STlight(5, X, Q)-STlight(4, X, Q) Else F2b=0D0 Endif Endif Else If ((Mode.eq.1 .or. Mode.eq.2) .and. Neff.ge.4) Then If (Iset .eq. Iset4F4) then F2c=STlight(4, X, Q)-STlight(3, X, Q) Else F2c=STheavy(4, X, Q) Endif Else F2c = 0D0 Endif If ((Mode.eq.1 .or. Mode.eq.3) .and. Neff.ge.5) Then F2b=STheavy(5, X, Q) Else F2b = 0D0 Endif If (Mode.eq.1 .and. Neff.ge.6) Then F2t=STheavy(6, X, Q) Else F2t = 0D0 Endif If (Mode.eq.1) F2total=STlight(3, X, Q)+F2c+F2b+F2t Endif If (Mode.eq.1) Then F2=F2total Elseif (Mode.eq.2) Then F2=F2c Elseif (Mode.eq.3) Then F2=F2b Else Print*, 'invalid mode for F2 calculation: Mode=',Mode Stop Endif Return End Subroutine SetF2(Isch0, Iset0) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (I0M=10, IVFN = 12, Mxsch=2, MxPDF=14) PARAMETER (MXQRK = 6) Common /Ischeme/ Isch, Iset Common /coupling/ Qcpln(MXQRK), SFNCF(MXQRK) Common /Nrd/ Ird Qcpln(1) = 0.66666667D0 Qcpln(2) =-0.33333333D0 Qcpln(3) = Qcpln(2) Qcpln(4) = Qcpln(1) Qcpln(5) = Qcpln(2) Qcpln(6) = Qcpln(1) DO 10 I=1,6 SFNCF(I)= Qcpln(I)**2 10 Continue Isch= Isch0 Iset= Iset0 If (Iset.eq.2 .or. Iset.eq.3) then Ird=1 Else Ird=2 Endif If (Iset .le. 0 .or. Iset.gt.MxPDF) Goto 100 If (Isch .lt. 0 .or. Isch.gt.Mxsch) Goto 101 If (Isch .eq. 0 .and. Iset.gt. I0M) Goto 99 If (Isch .eq. 1 .and. (Iset.gt. IVFN .or. Iset.le.I0M)) Goto 99 If (Isch .eq. 2 .and. Iset.le. IVFN) Goto 99 Call SetCtq4(Iset) Return 99 Print*,'The chosen PDF is not matched with the chosen scheme!' Stop 100 Print*,'The chosen PDF is not available!' Stop 101 Print*,'The chosen Scheme is not available!' Stop End Function Parton(Iparton, X, Q) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Common > / QCDtable / Alambda, Nfl, Iorder If (Iparton .lt. -Nfl .or. Iparton .gt. Nfl) Then Parton = 0D0 Return Endif Parton = Ctq4Pdf (Iparton, X, Q) Return End Function Alpi(Q) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Common > / QCDtable / Alambda, Nfl, Iorder Blambda= Alambda Nf= Nfl If (Q .lt. Bmass(Nf)) Blambda= 0.298D0 B0=12D0/(33D0-2D0*Nf) Aln=2D0*log(Q/Blambda) Alpi0=B0/Aln B1=6D0*(153D0-19D0*Nf)/(33D0-2D0*Nf)**2 Alpi=Alpi0*(1D0-B1*log(Aln)/Aln) Return End FUNCTION STheavy (IPRTN, X, Q) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXQRK = 6) Common /coupling/ Qcpln(MXQRK), SFNCF(MXQRK) F1m=Bmass(IPRTN) STheavy= X * Fnc(IPRTN,X,Q,F1m > ,QCPLN(IPRTN)/2,QCPLN(IPRTN)/2) RETURN END FUNCTION STlight(NEFF, X, Q) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXQRK = 6) Common /coupling/ Qcpln(MXQRK), SFNCF(MXQRK) Common /Nrd/ Ird NEFF0 = NEFF FLO = ST1 (NEFF0, X, Q) If(Ird.eq.1) then STlight = X * FLO Return Endif CPLGLU = 0.0 DO 10 IPRTN = 1, NEFF0 CPLGLU = CPLGLU + 2D0 * SFNCF(IPRTN) 10 CONTINUE FNLOQ = ST2 ( 1, NEFF0, X, Q) FNLOG = ST2 ( 0, NEFF0, X, Q) 98 FNLO = FNLOQ + FNLOG * CPLGLU 99 STlight = X* (FLO + FNLO) RETURN END FUNCTION ST1 (NEFF, X, Q) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MXQRK = 6) Common /coupling/ Qcpln(MXQRK), SFNCF(MXQRK) FLO = 0.0 DO 10 IPRTN = 1, NEFF FLO = FLO + SFNCF(IPRTN)* > (parton(IPRTN, X, Q) + parton(-IPRTN, X, Q)) 10 CONTINUE ST1 = FLO RETURN END FUNCTION ST2(IPRT, NEFF, X, Q) IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL C2Q0, FWINTG PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) PARAMETER (MXQRK = 6, NPOL = 3, NBSN = 4) COMMON 1 / STRFNC / XX, QQ, JPRT, NF DATA 1 CF, TR / 1.33333, 0.5 / 1 AER, RER / 0., 0.05 / NF = NEFF JPRT = IPRT XX = X QQ = Q AL2PI= ALPI(Q)/2. IF (IPRT.EQ.1) THEN F0 = CF * ADZINT(C2Q0, D0, X, AER, RER, ER, IR, 1, 1) 1 * ST1 (NF, X, Q) ELSE F0 = 0.0 ENDIF IF (IPRT .EQ. 0) THEN F1 = TR * ADZINT(FWINTG, X, D1, AER, RER, ER, IR, 2, 1) ELSEIF (IPRT .EQ. 1) THEN F1 = CF * ADZINT(FWINTG, X, D1, AER, RER, ER, IR, 2, 1) ENDIF STR2 = (F1 - F0) * AL2PI ST2 = STR2 RETURN END FUNCTION FWINTG (Y) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON 1 / STRFNC / XX, QQ, JPRT, NF C2Q(z) = ((1.+ Z**2) * ( LOG((1.- Z)/Z) - 3./4.) / (1. - z)) > + (9.+ 5.*Z) / 4. C2G(z) = (Z**2 + (1.- Z)**2) * LOG ((1.- Z)/Z) > - 1. + 8.* Z * (1.- Z) X = XX Q = QQ XY = X / Y IF (JPRT .EQ. 1) THEN TEM = C2Q(Y) * 1 ( ST1 (NF, XY, Q) / Y 1 - ST1 (NF, X, Q) ) ELSEIF (JPRT .EQ. 0) THEN TEM = C2G(Y) * parton(JPRT, XY, Q) / Y ENDIF FWINTG = TEM RETURN ENTRY C2Q0 (Y) C2Q0 = C2Q(Y) RETURN Entry Fintg(Y) X = XX Q = QQ XY = X / Y Fintg = C2G(Y) * parton(0, XY, Q) / Y Return END Double Precision Function Fnc > (IPRTN,X,Q,Fm,GLQ,GRQ) Double Precision X,Q,Fm,GLQ,GRQ Integer IPRTN Integer I Double Precision XHadNLo(-1:1) Double Precision XHadLo(-1:1),AHadLo(-1:1) Double Precision XHadSub(-1:1),AHadSub(-1:1) Double Precision Small,F1m,F2m,Hmass,Xmu,WHel2Ten Double Precision ScaleMu, parton, tempdf Integer Isch, Iset Common /Ischeme/ Isch, Iset Data Hmass,Small /.938D0, 1.D-3/ Fnc=0D0 F1m= MAX(Fm,Small) F2m= F1m Xmu= ScaleMu(F1m, F2m, Q) Call HADNLO (X,Q,F1m,F2m,GLQ,GRQ,Xmu,Hmass,XHadNLo) If (Isch.eq.1) then tempdf= parton(IPRTN,X,Xmu) If (tempdf .gt. 0D0) then Call HADLO (IPRTN,X,Q,F1m,F2m,GLQ,GRQ,Xmu,Hmass,XHadLo) Call HADLO (-IPRTN,X,Q,F1m,F2m,GRQ,GLQ,Xmu,Hmass,AHadLo) Call HADSUB (IPRTN,X,Q,F1m,F2m,GLQ,GRQ,Xmu,Hmass,XHadSub) Call HADSUB (-IPRTN,X,Q,F1m,F2m,GRQ,GLQ,Xmu,Hmass,AHadSub) Do 11 I=-1,1 XHadNLo(I)=XHadNLo(I)-XHadSub(I)-AHadSub(I) > + XHadLo(I)+AHadLo(I) 11 Continue Endif Endif Fnc=WHel2Ten(X,Q,Hmass,XHadNLo,2) Return End Function Bmass(Iparton) IMPLICIT DOUBLE PRECISION (A-H, O-Z) Common / Masstbl / Amass(6) Bmass=Amass(Iparton) Return End Double Precision Function ScaleMu(F1m, F2m, Q) Double Precision F1m, F2m, Q Double Precision Fm Double Precision scalc, scaln Common /ScalPar/ scalc, scaln Data scalc, scaln /.5D0, 2D0/ Fm= MAX(F1m, F2m) If (Q .gt. Fm) then ScaleMu= Sqrt(Fm**2 + scalc* Q**2 *(1D0- (Fm/Q)**2)**scaln) Else ScaleMu= Fm Endif Return End Double Precision Function WHel2Ten(X,Q,Hmass,XHel,ISF) Double Precision X,Q,Hmass,XHel(-1:1) Integer ISF Double Precision XH123(3) If (ISF.eq.0) then WHel2Ten=XHel(0) Return Endif Call HEL2TEN(X,Q,Hmass,XHel,XH123) If (ISF.eq.1) then WHel2Ten=XH123(1)*2 Elseif(ISF.eq.2) then WHel2Ten=XH123(2)/X Elseif(ISF.eq.3) then WHel2Ten=XH123(3) Else Print*,'invalid ISF=',ISF stop EndIf Return End SUBROUTINE HADLO (IPARTON,XBJ,Q,F1M,F2M,GLQ,GRQ,XMU,HMASS, > XHADLO) Implicit Double Precision (A-H, O-Z) Dimension OmGLO(-1:1), XHADLO(-1:1) DELTA(A,B,C) = SQRT(A**2 + B**2 + C**2 - 2.0*(A*B+B*C+C*A)) Q2=Q*Q Discr = 1.0 + (4.0*(HMASS*XBJ)**2)/Q2 ETA = 2.0* XBJ/(1.0 + Sqrt(Discr)) If (F1M .eq. F2M) then SHATTH= 4*F2M**2 Else SHATTH= F2M**2 Endif XITH=ETA*(Q2-F1M**2+SHATTH + > DELTA(-Q2,F1M**2,SHATTH))/(2.0*Q2) If (XITH.gt.1.) XITH=1.0 CALL PARLO (Q,F1M,F2M,GLQ,GRQ,OmGLO) PdfTMP = Parton(IPARTON, XITH, XMU) XHADLO( 1) = OmGLO( 1) * PDFTMP XHADLO( 0) = OmGLO( 0) * PDFTMP XHADLO(-1) = OmGLO(-1) * PDFTMP RETURN END SUBROUTINE PARLO (Q,F1M,F2M,GLQ,GRQ,OmGLO) Implicit Double Precision (A-H, O-Z) Dimension OmGLO(-1:1), OHelLO(-1:1,-1:1) Data iLeft, iLong, iRight / -1, 0, 1 / > igL2, igRgL, igR2 / -1, 0, 1 / DELTA(A,B,C) = SQRT(A**2 + B**2 + C**2 - 2.0*(A*B+B*C+C*A)) F1m2 = F1m**2 F2m2 = F2m**2 Q2 = Q**2 DEL = DELTA(-Q2,F1M2,F2M2) OHelLO(iRight, igR2) = (Q2+F1m2+F2m2 + DEL)/DEL OHelLO(iRight, igRgL) = -2*F1m*F2m/DEL OHelLO(iRight, igL2) = (Q2+F1m2+F2m2 - DEL)/DEL OHelLO(iLeft, igR2) = OHelLO(iRight, igL2) OHelLO(iLeft, igRgL) = OHelLO(iRight, igRgL) OHelLO(iLeft, igL2) = OHelLO(iRight, igR2) OHelLO(iLong, igR2) = ((F1m2+F2m2) + ((F1m2-F2m2)**2/Q2))/DEL OHelLO(iLong, igRgL) = +2*F1m*F2m/DEL OHelLO(iLong, igL2) = OHelLO(iLong, igR2) WR= GRQ**2 *OHelLO(iRight, igR2) + > 2.0*GRQ*GLQ*OHelLO(iRight, igRgL) + > GLQ**2 *OHelLO(iRight, igL2) WZ= GRQ**2 *OHelLO(iLong, igR2) + > 2.0*GRQ*GLQ*OHelLO(iLong, igRgL) + > GLQ**2 *OHelLO(iLong, igL2) WL= GRQ**2 *OHelLO(iLeft, igR2) + > 2.0*GRQ*GLQ*OHelLO(iLeft, igRgL) + > GLQ**2 *OHelLO(iLeft, igL2) OmGLO( 1)=WR OmGLO( 0)=WZ OmGLO(-1)=WL Return End SUBROUTINE HADSUB (IPARTON,XBJ,Q,F1M,F2M,GLQ,GRQ,XMU,HMASS, > XHADSUB) Implicit Double Precision (A-H, O-Z) Dimension OmGLO(-1:1), XHADSUB(-1:1) EXTERNAL XINTSUB Common / CINTSUB / XMUIN, XiTHIN DATA AERR, RERR, iActL, iActU / 0.D0, 1.0D-3, 1, 1 / DELTA(A,B,C) = SQRT(A**2 + B**2 + C**2 - 2.0*(A*B+B*C+C*A)) XHADSUB( 1) = 0.0 XHADSUB( 0) = 0.0 XHADSUB(-1) = 0.0 IF(XMU.LE.F1M) RETURN XMUIN =XMU Q2=Q*Q Discr = 1.0 + (4.0*(HMASS*XBJ)**2)/Q2 ETA = 2.0* XBJ/(1.0 + Sqrt(Discr)) If (F1M .eq. F2M) then SHATTH= 4*F2M**2 Else SHATTH= F2M**2 Endif XITH=ETA*(Q2-F1M**2+SHATTH + > DELTA(-Q2,F1M**2,SHATTH))/(2.0*Q2) XiTHIN =XITH ALPIX = Alpi(XMU) IF(ALPIX.GT.0.15) ALPIX=0.30 DivFct = ALPIX * Log(XMU/F1m) D1=1.0 If (XITH .lt. D1) then Convol = AdzInt (XINTSUB, XITH, D1, > AERR, RERR, ErrEst, iErr, iActL, iActU) Else Convol = 0. Endif PDFSUB= DIVFCT * CONVOL CALL PARLO (Q,F1M,F2M,GLQ,GRQ,OmGLO) XHADSUB( 1) = OmGLO( 1) * PDFSUB XHADSUB( 0) = OmGLO( 0) * PDFSUB XHADSUB(-1) = OmGLO(-1) * PDFSUB RETURN END Function XINTSUB(X) Implicit Double Precision (A-H, O-Z) Common / CINTSUB / XMU, XiTH iGluon=0 iHadr=1 XX=XiTH/X If (XX.le.1.) then GluPdf = Parton(iGluon,XX, XMU) SpltFn = 0.5 - x + x*x XINTSUB = SpltFn * GluPdf / X Else XINTSUB = 0. Endif Return End SUBROUTINE HADNLO(XBJ,Q,F1M,F2M,GLQ,GRQ,XMU,HMASS,XHLR) Implicit Double Precision (A-H, O-Z) Dimension XHLR(-1:1) Do 11 I=-1,1 XHLR(I)=FHADHEL ( I,XBJ,Q,F1M,F2M,GLQ,GRQ,XMU,HMASS) 11 Continue RETURN END Function FHADHEL (IHEL,XBJ,Q,F1M,F2M,GLQ,GRQ,XMU,HMASS) Implicit Double Precision (A-H, O-Z) COMMON /CXINTNLO/ ETAout,Qout,F1Mout,F2Mout,GLQout, > GRQout,XMUout,IRETout,IHELout EXTERNAL XINTNLO DATA AERR, RERR, iActL, iActU / 0.D0, 1.0D-3, 1, 1 / Q2 = Q*Q Discr = 1.0 + (4.0*(HMASS*XBJ)**2)/Q2 ETA = 2.0* XBJ/(1.0 + Sqrt(Discr)) IHELout=IHEL ETAout =ETA Qout =Q F1Mout =F1M F2Mout =F2M GLQout =GLQ GRQout =GRQ XMUout =XMU XIMIN= ETA*((F1M+F2M)**2+Q2)/Q2 XIMAX= 1.0 If (XIMIN.lt.XIMAX) then Conv = AdzInt (XINTNLO, XIMIN, XIMAX, > AERR, RERR, ErrEst, iErr, iActL, iActU) Else CONV=0. Endif FHADHEL = CONV RETURN END Function XINTNLO(XI) Implicit Double Precision (A-H, O-Z) Dimension OmGNLO(-1:1) COMMON /CXINTNLO/ ETA,Q,F1M,F2M,GLQ,GRQ,XMU,IRET,IHEL IGLUON=0 XHAT = ETA/XI IF(XI.le.1.) then CALL PARNLO (XHAT,Q,F1M,F2M,GLQ,GRQ,XMU,OmGNLO) PdfTMP = Parton(IGLUON, XI, XMU) XINTNLO= OmGNLO(IHEL) * PDFTMP / XI Else XINTNLO = 0. Endif RETURN END SUBROUTINE PARNLO (X,Q,F1M,F2M,GLQ,GRQ,XMU,OmGNLO) Implicit Double Precision (A-H, O-Z) Dimension OmGNLO(-1:1) DELTA(A,B,C) = SQRT(A**2 + B**2 + C**2 - 2.0*(A*B+B*C+C*A)) OmGNLO( 1)=0.0 OmGNLO( 0)=0.0 OmGNLO(-1)=0.0 Q2= Q*Q F1M2 = F1M ** 2 F2M2 = F2M ** 2 XMAX=Q2/((F1M+F2M)**2+Q2) IF((X.GE.1.0).OR.(X.GE.XMAX)) RETURN HMASS=0 !*** TEMP PATCH: THIS IS PARTON LEVEL SMIN=((F1M+F2M+HMASS)**2-HMASS**2) QMIN=SQRT( SMIN*X/(1.0-X) ) IF(Q.LE.QMIN) RETURN s = Q2 * (1./X - 1.) RS= SQRT(s) DEL = DELTA(s, F1M2, F2M2) AK = ( s + Q2) /(2.0*RS) Eq = ( s - Q2) /(2.0*RS) P = DEL /(2.0*RS) E1 = ( s - F2M2 + F1M2) /(2.0*RS) E2 = ( s + F2M2 - F1M2) /(2.0*RS) TLog = Log((s+F1M2-F2M2+Del)**2/(4*F1M2*s)) ULog = Log((s-F1M2+F2M2+Del)**2/(4*F2M2*s)) WsPlus = TLog * ( 0.5 + E1/AK * (E1/AK - 1.)) > + ULog * ( 0.5 + E2/AK * (E2/AK - 1.)) > - 2.*P/RS * ( Eq/AK )**2 WxPlus = (TLog + ULog)* (2.*F1M*F2M)/(Q2+s)**2 * (s-(F2M2+F1M2)) > - 8.*F1M*F2M*P*RS/(Q2+s)**2 WaPlus = TLog*(.5+E1/AK*(E1/AK-1.)+ F1M2*(F2M2-F1M2)/(2.*s*AK**2)) > - ULog*(.5+E2/AK*(E2/AK-1.)+ F2M2*(F1M2-F2M2)/(2.*s*AK**2)) > + P*(F2M2-F1M2)/(AK**2 * RS) WsZero = -(TLog + ULog) > *((F2M2+F1M2)*(Q2*Q2-(F2M2-F1M2)**2-2.*s*AK**2)/(4.*s*Q2*AK**2) > -Eq/(2.*AK**2 *RS)* (F2M2+F1M2-(F2M2-F1M2)**2/Q2) > +F2M2*F1M2/(AK**2 *s) ) > -(TLog-ULog)* Eq*(F2M2-F1M2)/(AK**2 *RS) > +P/(AK**2*Q2*RS)* ( (F2M2-F1M2)**2+ Q2*(2.*Q2-(F2M2+F1M2)) ) WxZero = -(TLog+ULog)*F1M*F2M > *( 1./Q2 + (s-(F2M2+F1M2))/(2.*AK**2*s) ) > + 2.*F1M*F2M*P/(AK**2 * RS) WaZero = 0. Xfac = Alpi(XMu)/2. Cs = GRQ*GRQ + GLQ*GLQ Ca = GRQ*GRQ - GLQ*GLQ Cx = 2.*GRQ*GLQ OmGNLO( 1)=Xfac*(Cs*WsPlus + Cx*WxPlus + Ca*WaPlus) OmGNLO( 0)=Xfac*(Cs*WsZero + Cx*WxZero + Ca*WaZero) OmGNLO(-1)=Xfac*(Cs*WsPlus + Cx*WxPlus - Ca*WaPlus) Return End SUBROUTINE HEL2TEN(XBJ,Q,HMASS,XHLR,XH123) Implicit Double Precision (A-H, O-Z) Dimension XH123(3), XHLR(-1:1) RHO=1.0 IF(HMASS.GT.0.0) THEN XNU=Q**2/(2.0*HMASS*XBJ) RHO=Sqrt(1.0+(Q/XNU)**2) ENDIF XH123(1)=(1./2.)*(XHLR(1)+XHLR(-1)) XH123(2)=(XBJ/RHO**2)*(XHLR(1)+XHLR(-1)+2.0*XHLR(0)) XH123(3)=(1./RHO)*(-XHLR(1)+XHLR(-1)) RETURN END FUNCTION ADZINT (F, A, B, AERR, RERR, ERREST, IER, IACTA, IACTB) IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL F PARAMETER (MAXINT = 1000) COMMON / ADZWRK / U(MAXINT), V(MAXINT), FU(MAXINT), ERS, RES, > FW(MAXINT), ERR(MAXINT), RESULT(MAXINT), FV(MAXINT), FA, FB, > ICTA, ICTB, NUMINT, IB SAVE / ADZWRK / DATA SMLL / 1E-20 / IER = 0 IF (AERR.LE.SMLL .AND. RERR.LE.SMLL) 1 STOP 'Both Aerr and Rerr are zero in ADZINT!' IF (IACTA.LT.0 .OR. IACTA.GT.2) THEN PRINT '(A, I4/ A)', ' Illegal value of IACT in ADZINT call', > 'IACTA =', IACTA, ' IACTA set for regular open-end option.' IACTA = 1 IER = 2 ENDIF IF (IACTB.LT.0 .OR. IACTB.GT.2) THEN PRINT '(A, I4/ A)', ' Illegal value of IACT in ADZINT call', > 'IACTB =', IACTB, ' IACTB set for regular open-end option.' IACTB = 1 IER = 3 ENDIF ICTA = IACTA ICTB = IACTB NUMINT = 3 DX = (B-A)/ NUMINT DO 10 I = 1, NUMINT IF (I .EQ. 1) THEN U(1) = A IF (IACTA .EQ. 0) THEN FU(1) = F(U(1)) ELSE FA = F(A+DX/2.) ENDIF ELSE U(I) = V(I-1) FU(I) = FV(I-1) ENDIF IF (I .EQ. NUMINT) THEN V(I) = B IF (IACTB .EQ. 0) THEN FV(I) = F(V(I)) ELSE IB = I FB = F(B-DX/2.) ENDIF ELSE V(I) = A + DX * I FV(I) = F(V(I)) ENDIF CALL ADZCAL(F,I) 10 CONTINUE CALL TOTALZ 30 TARGET = ABS(AERR) + ABS(RERR * RES) IF (ERS .GT. TARGET) THEN NUMOLD = NUMINT DO 40, I = 1, NUMINT IF (ERR(I)*NUMOLD .GT. TARGET) CALL ADZSPL(F,I,IER) 40 CONTINUE IF (IER.EQ.0 .AND. NUMINT.NE.NUMOLD) GOTO 30 ENDIF ADZINT = RES ERREST = ERS RETURN END SUBROUTINE ADZCAL (F,I) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D1 = 1.0, D2 = 2.0, HUGE = 1.E15) EXTERNAL F PARAMETER (MAXINT = 1000) COMMON / ADZWRK / U(MAXINT), V(MAXINT), FU(MAXINT), ERS, RES, > FW(MAXINT), ERR(MAXINT), RESULT(MAXINT), FV(MAXINT), FA, FB, > ICTA, ICTB, NUMINT, IB SAVE / ADZWRK / DX = V(I) - U(I) W = (U(I) + V(I)) / 2. IF (I .EQ. 1 .AND. ICTA .GT. 0) THEN FW(I) = FA FA = F (U(I) + DX / 4.) CALL SGLINT (ICTA, FA, FW(I), FV(I), DX, TEM, ER) ELSEIF (I .EQ. IB .AND. ICTB .GT. 0) THEN FW(I) = FB FB = F (V(I) - DX / 4.) CALL SGLINT (ICTB, FB, FW(I), FU(I), DX, TEM, ER) ELSE FW(I) = F(W) TEM = DX * (FU(I) + 4. * FW(I) + FV(I)) / 6. ER = DX * (FU(I) - 2. * FW(I) + FV(I)) / 12. ENDIF RESULT(I) = TEM ERR (I) = ABS (ER) RETURN END SUBROUTINE SGLINT (IACT, F1, F2, F3, DX, FINT, ESTER) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) DATA HUGE / 1.E20 / TEM = DX * (4.*F1 + 3.*F2 + 2.*F3) / 9. ER = DX * (4.*F1 - 6.*F2 + 2.*F3) / 9. IF (IACT .EQ. 2) THEN T1 = F2 - F1 T2 = F3 - F2 IF (T1*T2 .LE. 0.) GOTO 7 T3 = T2 - T1 IF (ABS(T3)*HUGE .LT. T1**2) GOTO 7 CC = LOG (T2/T1) / LOG(D2) IF (CC .LE. -D1) GOTO 7 BB = T1**2 / T3 AA = (F1*F3 - F2**2) / T3 TMP = DX * (AA + BB* 4.**CC / (CC + 1.)) ER = TEM - TMP TEM= TMP ENDIF 7 FINT = TEM ESTER= ER RETURN END SUBROUTINE ADZSPL (F, I, IER) IMPLICIT DOUBLE PRECISION (A-H, O-Z) EXTERNAL F PARAMETER (MAXINT = 1000) COMMON / ADZWRK / U(MAXINT), V(MAXINT), FU(MAXINT), ERS, RES, > FW(MAXINT), ERR(MAXINT), RESULT(MAXINT), FV(MAXINT), FA, FB, > ICTA, ICTB, NUMINT, IB SAVE / ADZWRK / DATA TINY / 1.D-20 / IF (NUMINT .GE. MAXINT) THEN IER = 1 RETURN ENDIF NUMINT = NUMINT + 1 IF (I .EQ. IB) IB = NUMINT U(NUMINT) = (U(I) + V(I)) / 2. V(NUMINT) = V(I) FU(NUMINT) = FW(I) FV(NUMINT) = FV(I) V(I) = U(NUMINT) FV(I) = FU(NUMINT) OLDRES = RESULT(I) OLDERR = ERR(I) CALL ADZCAL (F, I) CALL ADZCAL (F, NUMINT) DELRES = RESULT(I) + RESULT(NUMINT) - OLDRES RES = RES + DELRES GODERR = ABS(DELRES) ERS = ERS + GODERR - OLDERR SUMERR = ERR(I) + ERR(NUMINT) IF (SUMERR .GT. TINY) THEN FAC = GODERR / SUMERR ELSE FAC = 1. ENDIF ERR(I) = ERR(I) * FAC ERR(NUMINT) = ERR(NUMINT) * FAC RETURN END SUBROUTINE TOTALZ IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (MAXINT = 1000) COMMON / ADZWRK / U(MAXINT), V(MAXINT), FU(MAXINT), ERS, RES, > FW(MAXINT), ERR(MAXINT), RESULT(MAXINT), FV(MAXINT), FA, FB, > ICTA, ICTB, NUMINT, IB SAVE / ADZWRK / RES = 0. ERS = 0. DO 10 I = 1, NUMINT RES = RES + RESULT(I) ERS = ERS + ERR(I) 10 CONTINUE END