CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C A program for calculating the thrust and C C heavy jet mass distributions in e+e- annihilation C C using dressed gluon exponentiation C C C C D G E S H A P E C C C C version 0.1, March 12 2002 C C C C Authors: C C Einan Gardi, CERN theory division, Einan.Gardi@cern.ch C C Johan Rathsman, Uppsala University, Johan.Rathsman@tsl.uu.se C C C C Availability: on request or via the WWW on C C http://www3.tsl.uu.se/~rathsman/dgeshape/ C C C C References: C C [1] E. Gardi and J. Rathsman, Nuclear Physics B 609 (2001) 123 C C [2] E. Gardi and J. Rathsman, The thrust and heavy-jet mass C C distributions in the two-jet region, hep-ph/0201019 C C C C Please report any problems or suggestions for improvements. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C----------------------------------------------------------------------- PROGRAM DGESHAPE IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,YVAL,DSIG,R,LOGR INTEGER I REAL*8 DSIGDY,RY REAL*8 LOGRTFOPTNLO,LOGRTDGENLO,LOGRTDGENLOML REAL*8 LOGRHFOPTNLO,LOGRHDGENLO,LOGRHDGENLOML C-- Short instructions for use: C-- C-- Compile with the cernlib library mathlib C-- C-- The program is initiated by reading the results of the C-- Maple calculation for the observable that is requested. C-- This is done by calling the subroutine READCALC('obsble') C-- where the argument 'obsble' is one of the following C-- 'thrust' - thrust C-- 'thrust-lm' - ditto with modified logR matching, L=log(1/t-1)) C-- 'heavy' - heavy jet mass C-- 'heavy-lm' - ditto with modified logR matching, L=log(1/rho_H-1)) C-- C-- The following DOUBLE PRECISION functions are implemented: C-- DSIGDY - the differential cross-section C-- RY - the integrated cross-section (from 0 to y) C-- both are functions of the DOUBLE PRECISION variables C-- ALPHAMZ - alpha_s(M_Z) in the MSbar scheme C-- QVAL - centre-of-mass energy in GeV C-- YVAL - value of the variable, C-- either t=1-thrust or rho_H=m_H^2/Q^2 C-- C-- Comments: C-- 1) The functions DISGDY and RY will give either the thrust C-- or the heavy jet mass cross-sections depending which of C-- them that has been read with READCALC C-- 2) The resummed cross-sections have been matched with the C-- exact LO and NLO coefficients using so called logR C-- (or modified logR) matching (cf. Eq (106) in [1]) C-- C-- The NLO results for logR using the exact coefficients or the C-- expanded resummed result (with or without the modified log) C-- are avialable in the following DOUBLE PRECISION functions with C-- the same arguments as the above. C-- for the thrust: C-- LOGRTFOPTNLO - exact LO and NLO coefficients C-- LOGRTDGENLO - DGE expanded to NLO (cf. Eq (107) in [1]) C-- LOGRTDGENLOLM - DGE expanded to NLO with modified log, L=log(1/t-1)) C-- and the heavy jet mass: C-- LOGRHFOPTNLO - exact LO and NLO coefficients C-- LOGRHDGENLO - DGE expanded to NLO C-- LOGRHDGENLOLM - DGE expanded to NLO with modified log, L=log(1/rho_H-1)) C-- These functions can for example be used to undo the logR-matching C-- Example 1: write the differential and integrated C-- cross-section for thrust with a modified logR matching C-- as a function of y(=1-thrust) CALL READCALC('thrust-lm') QVAL=91.2D0 ALPHAMZ=0.110D0 WRITE(*,*) 'Results for thrust with a modfied log:' WRITE(*,*) ' t=1-thrust dsigma/dt R_t ' DO 10 I=0,35 YVAL=0.01D0*I DSIG=DSIGDY(ALPHAMZ,QVAL,YVAL) R=RY(ALPHAMZ,QVAL,YVAL) WRITE(*,1000) YVAL,DSIG,R 10 CONTINUE 1000 FORMAT(' ',F5.2,' ',F8.4,' ',F7.4) C-- Example 2: compare the values of logR for the heavy jet mass C-- with the NLO results using the exact coefficients or the C-- expanded resummed result CALL READCALC('heavy') QVAL=91.2D0 ALPHAMZ=0.110D0 WRITE(*,*) 'Results for heavy jet mass:' WRITE(*,*) ' rho_H logR_H' WRITE(*,*) ' matched FOPTNLO DGENLO ' DO 20 I=1,35 YVAL=0.01D0*I LOGR=LOG(RY(ALPHAMZ,QVAL,YVAL)) WRITE(*,2000) YVAL,LOGR, & LOGRHFOPTNLO(ALPHAMZ,QVAL,YVAL),LOGRHDGENLO(ALPHAMZ,QVAL,YVAL) 20 CONTINUE 2000 FORMAT(' ',F5.2,' ',F8.4,' ',F8.4,' ',F8.4) END C----------------------------------------------------------------------- C-- Below follows the routines that can be called from the main C-- program as explained above and some internally used ones C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DSIGDY(ALPHAMZ,QVAL,YVAL) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,YVAL INTEGER NLOW REAL*8 HVAL,LVAL REAL*8 DPT(401),RPT(401) COMMON /RESULT/ DPT,RPT CALL CALCDSIG(ALPHAMZ,QVAL,HVAL) LVAL=LOG(1.D0/YVAL) NLOW=1+INT((LVAL-0.8D0)/HVAL+1.D-9) C-- interpolate in log(1/y) to get DSIGDY IF(NLOW.LE.400 .AND. NLOW.GE.1) THEN DSIGDY=DPT(NLOW)+ + ((LVAL-0.8D0)-HVAL*(NLOW-1))/HVAL*(DPT(NLOW+1)-DPT(NLOW)) ELSE DSIGDY=0.D0 ENDIF END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION RY(ALPHAMZ,QVAL,YVAL) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,YVAL INTEGER NLOW REAL*8 HVAL,LVAL REAL*8 DPT(401),RPT(401) COMMON /RESULT/ DPT,RPT CALL CALCDSIG(ALPHAMZ,QVAL,HVAL) LVAL=LOG(1.D0/YVAL) NLOW=1+INT((LVAL-0.8D0)/HVAL+1.D-9) C-- interpolate in log(1/y) to get RY IF(NLOW.LE.400 .AND. NLOW.GE.1) THEN RY=RPT(NLOW)+ + ((LVAL-0.8D0)-HVAL*(NLOW-1))/HVAL*(RPT(NLOW+1)-RPT(NLOW)) ELSE RY=0.D0 ENDIF END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION S2FCN(ALPHAMZ,QVAL) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL REAL*8 MZ,PI,BETA0,LAMBDA1LOOP,DELTA,S1MZ,LAMBDA2LOOP C-- calculate Lambda_2loop from alpha_s at MZ MZ=91.2D0 PI=3.14159265358979323D0 BETA0=11.D0/4.D0-5.D0/6.D0 LAMBDA1LOOP=MZ*EXP(-PI/2/ALPHAMZ/BETA0) DELTA=(51.D0-95.D0/3.D0)/8.D0/BETA0**2 S1MZ=LOG(MZ**2/LAMBDA1LOOP**2) LAMBDA2LOOP=LAMBDA1LOOP*EXP(DELTA/2.D0*LOG(1.D0+S1MZ/DELTA)) C-- calculate s2 for this Q S2FCN=LOG(QVAL**2/LAMBDA2LOOP**2) END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION ALPHAFCN(ALPHAMZ,QVAL) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL REAL*8 ALPHAVAL(500:2500) COMMON /ALPHAMSBAR/ ALPHAVAL REAL*8 S2FCN REAL*8 AINV,S2LOW,S2THIS INTEGER NLOW C-- get s2 for this Q S2THIS=S2FCN(ALPHAMZ,QVAL) C-- interpolate to get alpha_s NLOW=INT(S2THIS*100.D0+1.D-10) S2LOW=NLOW/100.D0 IF(NLOW.LE.2500 .AND. NLOW.GE.500) THEN AINV=1.D0/ALPHAVAL(NLOW)+ +(1.D0/ALPHAVAL(NLOW+1)-1.D0/ALPHAVAL(NLOW))*(S2THIS-S2LOW)/0.01D0 ELSE WRITE(*,*) 'alpha_s out of range' WRITE(*,*) 'If this error occurs please contact the authors' STOP ENDIF ALPHAFCN=1.D0/AINV END C----------------------------------------------------------------------- SUBROUTINE CALCDSIG(ALPHAMZ,QVAL,HLOW) IMPLICIT NONE C-- input REAL*8 ALPHAMZ,QVAL C-- output REAL*8 HLOW REAL*8 DPT(401),RPT(401) COMMON /RESULT/ DPT,RPT REAL*8 S2VAL(80:200),DELTAL(80:200) INTEGER S1MIN,S1MAX COMMON /SDATA/ S2VAL,DELTAL,S1MIN,S1MAX REAL*8 LRPT(80:200,401),LDPT(80:200,401) COMMON /RESUMPT/ LRPT,LDPT REAL*8 S2FCN REAL*8 LRPTLOW(401),LRPTHIGH(401),RPTHIS(401) REAL*8 LDPTLOW(401),LDPTHIGH(401),DPTHIS(401) REAL*8 S2THIS,S2LOW,S2HIGH,HHIGH,LVAL,LSHIFT INTEGER S1LOW,S1HIGH INTEGER I,NLOW C-- get s2 for this Q S2THIS=S2FCN(ALPHAMZ,QVAL) C-- find right columns for interpolation in s2 IF (S2THIS.LT.S2VAL(S1MIN)) THEN S1LOW=S1MIN S1HIGH=S1MIN+1 HLOW=DELTAL(S1MIN) HHIGH=DELTAL(S1MIN+1) S2LOW=S2VAL(S1MIN) S2HIGH=S2VAL(S1MIN+1) WRITE(*,*) WRITE(*,*) 'WARNING: extrapolating instead of interpolating' WRITE(*,*) 'If this error persists please contact the authors' ENDIF IF (S2THIS.GT.S2VAL(S1MAX-1)) THEN S1LOW=S1MAX-1 S1HIGH=S1MAX HLOW=DELTAL(S1MAX-1) HHIGH=DELTAL(S1MAX) S2LOW=S2VAL(S1MAX-1) S2HIGH=S2VAL(S1MAX) WRITE(*,*) WRITE(*,*) 'WARNING: extrapolating instead of interpolating' WRITE(*,*) 'If this error persists please contact the authors' ENDIF DO 10 I=S1MIN,S1MAX-1 IF (S2THIS.GE.S2VAL(I) .AND. S2THIS.LE.S2VAL(I+1)) THEN S1LOW=I S1HIGH=I+1 HLOW=DELTAL(I) HHIGH=DELTAL(I+1) S2LOW=S2VAL(I) S2HIGH=S2VAL(I+1) ENDIF 10 CONTINUE IF(S1LOW.LT.S1MIN .OR. S1LOW.GT.S1MAX .OR. / S1HIGH.LT.S1MIN .OR. S1HIGH.GT.S1MAX ) THEN WRITE(*,*) S1MIN,S1MAX WRITE(*,*) S1LOW,S1HIGH WRITE(*,*) S2THIS,S2VAL(S1MIN),S2VAL(S1MAX) ENDIF DO 20 I=1,401 LRPTLOW(I)=LRPT(S1LOW,I) LDPTLOW(I)=LDPT(S1LOW,I) 20 CONTINUE C-- interpolate in ln(1/t) for S1HIGH DO 30 I=1,401 LVAL=0.8D0+(I-1)*HLOW NLOW=1+INT((LVAL-0.8D0)/HHIGH+1.D-9) IF(NLOW.LE.400 .AND. NLOW.GE.1) THEN LRPTHIGH(I)=LRPT(S1HIGH,NLOW)+ + ((LVAL-0.8D0)-HHIGH*(NLOW-1))/HHIGH* * (LRPT(S1HIGH,NLOW+1)-LRPT(S1HIGH,NLOW)) LDPTHIGH(I)=LDPT(S1HIGH,NLOW)+ + ((LVAL-0.8D0)-HHIGH*(NLOW-1))/HHIGH* * (LDPT(S1HIGH,NLOW+1)-LDPT(S1HIGH,NLOW)) ELSE LRPTHIGH(I)=0.D0 LDPTHIGH(I)=0.D0 ENDIF 30 CONTINUE C-- interpolate in s2 DO 40 I=1,401 RPTHIS(I)=EXP(LRPTLOW(I)+ + (S2THIS-S2LOW)/(S2HIGH-S2LOW)*(LRPTHIGH(I)-LRPTLOW(I))) RPT(I)=RPTHIS(I) DPTHIS(I)=EXP(LDPTLOW(I)+ + (S2THIS-S2LOW)/(S2HIGH-S2LOW)*(LDPTHIGH(I)-LDPTLOW(I))) DPT(I)=DPTHIS(I) 40 CONTINUE END C----------------------------------------------------------------------- SUBROUTINE READCALC(OBSBLE) IMPLICIT NONE CHARACTER*(*) OBSBLE REAL*8 S2VAL(80:200),DELTAL(80:200) INTEGER S1MIN,S1MAX COMMON /SDATA/ S2VAL,DELTAL,S1MIN,S1MAX REAL*8 LRPT(80:200,401),LDPT(80:200,401) COMMON /RESUMPT/ LRPT,LDPT REAL*8 ALPHAVAL(500:2500) COMMON /ALPHAMSBAR/ ALPHAVAL REAL*8 S1,S2,LT(401),LR,LDS CHARACTER*20 FILE(80:200) INTEGER I,J S1MIN=80 S1MAX=200 FILE(80)='80' FILE(81)='81' FILE(82)='82' FILE(83)='83' FILE(84)='84' FILE(85)='85' FILE(86)='86' FILE(87)='87' FILE(88)='88' FILE(89)='89' FILE(90)='90' FILE(91)='91' FILE(92)='92' FILE(93)='93' FILE(94)='94' FILE(95)='95' FILE(96)='96' FILE(97)='97' FILE(98)='98' FILE(99)='99' FILE(100)='100' FILE(101)='101' FILE(102)='102' FILE(103)='103' FILE(104)='104' FILE(105)='105' FILE(106)='106' FILE(107)='107' FILE(108)='108' FILE(109)='109' FILE(110)='110' FILE(111)='111' FILE(112)='112' FILE(113)='113' FILE(114)='114' FILE(115)='115' FILE(116)='116' FILE(117)='117' FILE(118)='118' FILE(119)='119' FILE(120)='120' FILE(121)='121' FILE(122)='122' FILE(123)='123' FILE(124)='124' FILE(125)='125' FILE(126)='126' FILE(127)='127' FILE(128)='128' FILE(129)='129' FILE(130)='130' FILE(131)='131' FILE(132)='132' FILE(133)='133' FILE(134)='134' FILE(135)='135' FILE(136)='136' FILE(137)='137' FILE(138)='138' FILE(139)='139' FILE(140)='140' FILE(141)='141' FILE(142)='142' FILE(143)='143' FILE(144)='144' FILE(145)='145' FILE(146)='146' FILE(147)='147' FILE(148)='148' FILE(149)='149' FILE(150)='150' FILE(151)='151' FILE(152)='152' FILE(153)='153' FILE(154)='154' FILE(155)='155' FILE(156)='156' FILE(157)='157' FILE(158)='158' FILE(159)='159' FILE(160)='160' FILE(161)='161' FILE(162)='162' FILE(163)='163' FILE(164)='164' FILE(165)='165' FILE(166)='166' FILE(167)='167' FILE(168)='168' FILE(169)='169' FILE(170)='170' FILE(171)='171' FILE(172)='172' FILE(173)='173' FILE(174)='174' FILE(175)='175' FILE(176)='176' FILE(177)='177' FILE(178)='178' FILE(179)='179' FILE(180)='180' FILE(181)='181' FILE(182)='182' FILE(183)='183' FILE(184)='184' FILE(185)='185' FILE(186)='186' FILE(187)='187' FILE(188)='188' FILE(189)='189' FILE(190)='190' FILE(191)='191' FILE(192)='192' FILE(193)='193' FILE(194)='194' FILE(195)='195' FILE(196)='196' FILE(197)='197' FILE(198)='198' FILE(199)='199' FILE(200)='200' C-- check that the user has made a proper choice IF(OBSBLE.EQ.'thrust') THEN WRITE(*,*) 'Reading calculation for thrust' ELSEIF (OBSBLE.EQ.'thrust-lm') THEN WRITE(*,*) 'Reading calculation for thrust with modified log' ELSEIF (OBSBLE.EQ.'heavy') THEN WRITE(*,*) 'Reading calculation for heavy jet mass' ELSEIF (OBSBLE.EQ.'heavy-lm') THEN WRITE(*,*) & 'Reading calculation for heavy jet mass with modified log' ELSE WRITE(*,*) 'Observable ',OBSBLE,' not recognized:' WRITE(*,*) 'choose from thrust, thrust-lm, heavy or heavy-lm' STOP ENDIF C-- read appropriate calculation DO 200 J=S1MIN,S1MAX IF(OBSBLE.EQ.'thrust') THEN OPEN (UNIT=1,NAME='thrust/res_7_'//FILE(J),STATUS='OLD') ELSEIF (OBSBLE.EQ.'thrust-lm') THEN OPEN (UNIT=1, & NAME='thrust/res_7_logmod'//FILE(J),STATUS='OLD') ELSEIF (OBSBLE.EQ.'heavy') THEN OPEN (UNIT=1,NAME='heavy/hres_7_'//FILE(J),STATUS='OLD') ELSEIF (OBSBLE.EQ.'heavy-lm') THEN OPEN (UNIT=1, & NAME='heavy/hres_7_logmod'//FILE(J),STATUS='OLD') ENDIF DO 210 I=401,1,-1 READ(1,*) S1,S2,LT(I),LR,LDS LRPT(J,I)=LR LDPT(J,I)=LDS 210 CONTINUE S2VAL(J)=S2 DELTAL(J)=LT(2)-LT(1) CLOSE(1) 200 CONTINUE C-- read alpha_s in MSbar scheme OPEN (UNIT=1,NAME='alphamsbar',STATUS='OLD') DO 300 I=500,2500 READ(1,*) S2,ALPHAVAL(I) 300 CONTINUE END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION LOGRTDGENLO(ALPHAMZ,QVAL,TVAL) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,TVAL REAL*8 ALPHAFCN REAL*8 ALPHAS,LVAL,t0 ALPHAS=ALPHAFCN(ALPHAMZ,QVAL) LVAL=LOG(1.D0/TVAL) t0 = (0.6964364875294484D0*LVAL-0.6317084250556305D0*LVAL**2.D0 #-0.258931913752641D0*LVAL**3.D0+0.2187561762473097D0)*ALPHAS**2 #+(0.6366197723675813D0*LVAL- #0.4244131815783876D0*LVAL**2.D0)*ALPHAS LOGRTDGENLO=t0 END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION LOGRTDGENLOLM(ALPHAMZ,QVAL,TVAL) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,TVAL REAL*8 ALPHAFCN REAL*8 ALPHAS,LVAL,t0 ALPHAS=ALPHAFCN(ALPHAMZ,QVAL) LVAL=LOG(1.D0/TVAL-1.D0) t0 = (0.6964364875294484D0*LVAL-0.6317084250556305D0*LVAL**2.D0 #-0.258931913752641D0*LVAL**3.D0+0.2187561762473097D0)*ALPHAS**2 #+(0.6366197723675813D0*LVAL- #0.4244131815783876D0*LVAL**2.D0)*ALPHAS LOGRTDGENLOLM=t0 END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION LOGRTFOPTNLO(ALPHAMZ,QVAL,TVAL) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,TVAL REAL*8 ALPHAFCN,MDILOG,HEAVISIDE REAL*8 ALPHAS,t0,s1,s2,s3,s4,s5,s6,s7 ALPHAS=ALPHAFCN(ALPHAMZ,QVAL) s1 = 0.3183098861837907D0*(-0.3859912089130969D1-0.266666666666666 #7D1*mdilog(1.D0/TVAL-1.D0)-0.2666666666666667D1*LOG(1.D0/TVAL-2.D0 #)*LOG(1.D0/TVAL-1.D0)-2.D0*LOG(1.D0/TVAL)+2.D0*LOG(1.D0/TVAL-2. #D0)-4.D0*LOG(1.D0/TVAL-2.D0)*TVAL+3.D0*TVAL**2+4.D0*TVAL-2.D0*LOG #(TVAL))*Heaviside(0.3333333333333333D0-TVAL)*ALPHAS s3 = 0.1013211836423378D0 s6 = Heaviside(0.3333333333333333D0-TVAL)*(2.D0*LOG(1.D0/TVAL)-0. #2837267235D3*TVAL**2+0.1529609156666667D4*TVAL**3-0.3021640375D4*T #VAL**4+0.1813300976D4*TVAL**5-0.790045865D1*LOG(TVAL)-2.D0*LOG(1 #.D0/TVAL-2.D0)*(1.D0/TVAL-2.D0)*TVAL-0.5222222222222222D1*LOG(1.D #0/TVAL)**3.D0-0.6270150593197792D1*LOG(1.D0/TVAL)**2.D0+0.2666666 #666666667D1*LOG(1.D0/TVAL-2.D0)*LOG(1.D0/TVAL-1.D0)+0.8888888888 #888889D0*LOG(1.D0/TVAL)**4.D0+0.2666666666666667D1*mdilog(1.D0/TVA #L-1.D0)+0.170214422447D2-0.240492465D2*TVAL) s7 = Heaviside(TVAL-0.3333333333333333D0)*(-0.6062331071503D10*TVA #L+0.3066571862015677D11*TVAL**2-0.1204954418833919D12*TVAL**3+0.34 #19722762046396D12*TVAL**4-0.6892410427568841D12*TVAL**5+0.96337902 #35407135D12*TVAL**6-0.8891300934944143D12*TVAL**7+0.48795734240613 #41D12*TVAL**8-0.1207297744444444D12*TVAL**9+0.265995433912651D9*LO #G(TVAL)+0.9998040145502484D9)-0.5D0*(-0.3859912089130969D1-0.2666 #666666666667D1*mdilog(1.D0/TVAL-1.D0)-0.2666666666666667D1*LOG(1.D #0/TVAL-2.D0)*LOG(1.D0/TVAL-1.D0)-2.D0*LOG(1.D0/TVAL)+2.D0*LOG(1 #.D0/TVAL-2.D0)-4.D0*LOG(1.D0/TVAL-2.D0)*TVAL+3.D0*TVAL**2+4.D0*TV #AL-2.D0*LOG(TVAL))**2.D0*Heaviside(0.3333333333333333D0-TVAL)**2. #D0 s5 = s6+s7 s6 = ALPHAS**2 s4 = s5*s6 s2 = s3*s4 t0 = s1+s2 LOGRTFOPTNLO=t0 END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION HEAVISIDE(X) IMPLICIT NONE REAL*8 X IF (X.GT.0.D0) THEN HEAVISIDE=1.D0 ELSE HEAVISIDE=0.D0 ENDIF END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION MDILOG(X) IMPLICIT NONE REAL*8 X,DDILOG C-- convert from Maple definition to Cernlib MDILOG=DDILOG(1.D0-X) END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION LOGRHDGENLO(ALPHAMZ,QVAL,RHO) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,RHO REAL*8 ALPHAFCN REAL*8 ALPHAS,LVAL,t0 ALPHAS=ALPHAFCN(ALPHAMZ,QVAL) LVAL=LOG(1.D0/RHO) t0 = (0.6850367656932638D0*LVAL+0.1093780881236548D0 #-0.258931913752641D0*LVAL**3.D0-0.3354121287593342D0*LVAL**2.D0) #*ALPHAS**2+(0.6366197723675813D0*LVAL #-0.4244131815783876D0*LVAL**2.D0)*ALPHAS LOGRHDGENLO=t0 END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION LOGRHDGENLOLM(ALPHAMZ,QVAL,RHO) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,RHO REAL*8 ALPHAFCN REAL*8 ALPHAS,LVAL,t0 ALPHAS=ALPHAFCN(ALPHAMZ,QVAL) LVAL=LOG(1.D0/RHO-1.D0) t0 = (0.6850367656932638D0*LVAL+0.1093780881236548D0 #-0.258931913752641D0*LVAL**3.D0-0.3354121287593342D0*LVAL**2.D0) #*ALPHAS**2+(0.6366197723675813D0*LVAL #-0.4244131815783876D0*LVAL**2.D0)*ALPHAS LOGRHDGENLOLM=t0 END C----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION LOGRHFOPTNLO(ALPHAMZ,QVAL,RHO) IMPLICIT NONE REAL*8 ALPHAMZ,QVAL,RHO REAL*8 ALPHAFCN,MDILOG,HEAVISIDE REAL*8 ALPHAS,t0,s1,s2,s3,s4,s5,s6,s7 ALPHAS=ALPHAFCN(ALPHAMZ,QVAL) s4 = 0.1013211836423378D0 s6 = Heaviside(0.3333333333333333D0-RHO) s7 = -0.3116265944444444D6*RHO**9-0.4623191475D5*RHO**8+0.19218225 #71428571D5*RHO**7+0.138786606D5*RHO**6+0.243708168D4*RHO**5-0.1511 #9318175D4*RHO**4-0.7198404833333333D3*RHO**3+0.4302372855D3*RHO**2 #-0.89606722D1*LOG(RHO)+0.20268364D1*LOG(1.D0/RHO)**3.D0*RHO+0.26 #66666666666667D1*LOG(-2.D0+1.D0/RHO)*LOG(1.D0/RHO-1.D0)+2.D0*LOG #(1.D0/RHO)+0.144142630850583D2+0.60805092D1*RHO*LOG(1.D0/RHO)**2 #.D0+0.121610184D2*RHO*LOG(1.D0/RHO)+0.26906478D7*RHO**11-0.522222 #2222222222D1*LOG(1.D0/RHO)**3.D0-0.569233561D6*RHO**10+0.88888888 #88888889D0*LOG(1.D0/RHO)**4.D0-0.3345823363245389D1*LOG(1.D0/RHO #)**2.D0+0.2666666666666667D1*mdilog(1.D0/RHO-1.D0)-0.1364245696D3* #RHO-2.D0*LOG(-2.D0+1.D0/RHO)*(-2.D0+1.D0/RHO)*RHO s5 = s6*s7 s3 = s4*s5 s4 = 0.1013211836423378D0*Heaviside(RHO-0.3333333333333333D0)*(-0. #6554658366666667D10*RHO**9+0.2730131976425D11*RHO**8-0.51275211157 #13143D11*RHO**7+0.572704836803356D11*RHO**6-0.4224043541556267D11* #RHO**5+0.2160669063308841D11*RHO**4-0.7848960119805455D10*RHO**3+0 #.2059364236095161D10*RHO**2-0.4197062062677662D9*RHO+0.18984013106 #05831D8*LOG(RHO)+0.7077878235786739D8)-0.5066059182116889D-1*(-0. #3859912089130969D1-0.2666666666666667D1*mdilog(1.D0/RHO-1.D0)-0.26 #66666666666667D1*LOG(-2.D0+1.D0/RHO)*LOG(1.D0/RHO-1.D0)-2.D0*LOG #(1.D0/RHO)-4.D0*LOG(-2.D0+1.D0/RHO)*RHO+2.D0*LOG(-2.D0+1.D0/RHO) #+3.D0*RHO**2+4.D0*RHO-2.D0*LOG(RHO))**2.D0*Heaviside(0.3333333333 #333333D0-RHO)**2.D0 s2 = s3+s4 s3 = ALPHAS**2 s1 = s2*s3 s2 = 0.3183098861837907D0*(-0.3859912089130969D1-0.266666666666666 #7D1*mdilog(1.D0/RHO-1.D0)-0.2666666666666667D1*LOG(-2.D0+1.D0/RHO) #*LOG(1.D0/RHO-1.D0)-2.D0*LOG(1.D0/RHO)-4.D0*LOG(-2.D0+1.D0/RHO) #*RHO+2.D0*LOG(-2.D0+1.D0/RHO)+3.D0*RHO**2+4.D0*RHO-2.D0*LOG(RHO) #)*Heaviside(0.3333333333333333D0-RHO)*ALPHAS t0 = s1+s2 LOGRHFOPTNLO=t0 END C-----------------------------------------------------------------------