C23456789012345678901234567890123456789012345678901234567890123456789012 C C ---------------------- C beowulf.f Version 0.7 C ---------------------- C C Program for numerical integration of jet cross sections in C electron-positron annihilation. C C Comment lines with CM are for machine dependent code C Comment lines for special code for testing purposes: CT C C D. E. Soper C PROGRAM EECALC c implicit none INTEGER SIZE PARAMETER (SIZE = 3) C Seed for random number generator. INTEGER SEED C Number of groups of points for RENO. INTEGER GROUPSIZE,NRENO COMMON /MONTECARLO/GROUPSIZE,NRENO C Global size and scale parameters. INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX REAL*8 ENERGYSCALE COMMON /MSCALE/ ENERGYSCALE C RENO results. REAL*8 SUMR,ERRORR,SUMI,ERRORI REAL*8 SUMBISR,ERRORBISR,SUMBISI,ERRORBISI REAL*8 SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI INTEGER INCLUDED,EXTRACOUNT,OMITTED INTEGER NVALUE(-9:6) REAL*8 VALUEMAX REAL*8 KBAD(0:3*SIZE-1,0:3) INTEGER BADGRAPHNUMBER C Limits. REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT C Newpoint parameters. REAL*8 AS,BS,AM,BM,AH,BH,AUV,BUV,CEPS COMMON /POINTS/ AS,BS,AM,BM,AH,BH,AUV,BUV,CEPS C Deform parameters. REAL*8 DFRM1,DFRM2,DFRM3,DFRM4 COMMON /DEFORMSCALES/DFRM1,DFRM2,DFRM3,DFRM4 C Smear parameters. REAL*8 SMEARFCTR INTEGER LOWPWR,HIGHPWR COMMON /SMEARPARMS/ SMEARFCTR,LOWPWR,HIGHPWR C Renormaliation parameters. REAL*8 MUOVERRTS COMMON /RENORMALIZE/ MUOVERRTS C Graphs to include LOGICAL USEGRAPH(10) COMMON /WHICHGRAPHS/ USEGRAPH C Diagnostic data instructions. LOGICAL REPORT,DETAILS COMMON /CALCULOOK/ REPORT,DETAILS C Hrothgar dummy variables: REAL*8 KCUT0(SIZE+1,0:3) C CHARACTER DATE*(28) REAL CPUTIME,DTIME,TIMEARRAY(2) C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT C INTEGER N,MU C C----------------------------------------------------------------------- C C Input and ouput units. NIN = 5 NOUT = 6 C Parameters for subroutine CALCULATE for throwing away bad points. BADNESSLIMIT = 1.0D8 CANCELLIMIT = 1.0D4 C Cutoff on 1 - THRUST. THRUSTCUT = 0.0D0 C How many points in a RENO group GROUPSIZE = 100 C Parameters for subroutine NEWPOINT AH = 3.0D0 BH = 3.0D0 AM = 2.0D0 BM = 2.0D0 AS = 1.0D0 BS = 2.0D0 AUV = 3.0D0 BUV = 1.0D0 CEPS = 1.0D0 C Parameters for subroutine DEFORM. DFRM1 = 0.7D0 DFRM2 = 1.0D0 DFRM3 = 1.0D0 DFRM4 = 1.0D0 C Parmameters for function SMEAR. SMEARFCTR = 3.0D0 LOWPWR = 9 HIGHPWR = 18 C Renormalization parameters. MUOVERRTS = 1.0D0 C Graphs to include. USEGRAPH(1) = .TRUE. USEGRAPH(2) = .TRUE. USEGRAPH(3) = .TRUE. USEGRAPH(4) = .TRUE. USEGRAPH(5) = .TRUE. USEGRAPH(6) = .TRUE. USEGRAPH(7) = .TRUE. USEGRAPH(8) = .TRUE. USEGRAPH(9) = .TRUE. USEGRAPH(10) = .TRUE. C Parameters to control diagnostic reporting. REPORT = .FALSE. DETAILS = .FALSE. C Global size parameters: Number of loops, number of propagators C number of vertices, max number of cut propagators. NLOOPS = 3 NPROPS = 3 * NLOOPS - 1 NVERTS = 2 * NLOOPS CUTMAX = NLOOPS + 1 C Energy scale. ENERGYSCALE = 1.0D0 C Call Hrothgar to tell him to get ready. C Here KCUT0 are just dummy variables for Hrothgar. DO N = 1,SIZE+1 DO MU = 0,3 KCUT0(N,MU) = 1.0D0 ENDDO ENDDO CALL HROTHGAR(1,KCUT0,1.0D0,1.0D0,'STARTCALC ') C C---Timing function for Sun workstations; remove for others C CPUTIME = DTIME(TIMEARRAY) C--- WRITE(NOUT,*)'Please give number of groups of Monte Carlo points.' READ(NIN,*) NRENO WRITE(NOUT,*)'Please give seed (0 ,' Beowulf version 0.7 ',A,/,50('-')) C C**************************************************** WRITE(NOUT,*)'Latest revision 9 April 1998' C**************************************************** C CALL VERSION C WRITE(NOUT,*)' ' WRITE(NOUT,*)'I put as much faith in my martial' WRITE(NOUT,*)'might and power as Grendel puts in his.' WRITE(NOUT,*)'Therefore Ill not slay him with a sword...' WRITE(NOUT,*)'... no sword on earth' WRITE(NOUT,*)'even the best of iron war-blades,' WRITE(NOUT,*)'could make a dent in that miscreant; ' WRITE(NOUT,*)'for he had worked a spell on weapons' WRITE(NOUT,*)'to blunt their edge...' WRITE(NOUT,*)' - Beowulf, translated by Stanley B. Greenfield' WRITE(NOUT,*)' ' WRITE(NOUT,102) NRENO,SEED 102 FORMAT('Number of groups of points', I7,' Seed',I8) WRITE(NOUT,111) 111 FORMAT('Parameters for choosing points:') WRITE(NOUT,112)AH,BH,AM,BM,AS,BS,AUV,BUV,CEPS 112 FORMAT(' Hard map: A =',1P G9.2,' B = ',1P G9.2,/, > ' Medium map: A =',1P G9.2,' B = ',1P G9.2,/, > ' Soft map: A =',1P G9.2,' B = ',1P G9.2,/, > ' UV map: A =',1P G9.2,' B = ',1P G9.2,/, > ' Soft angles: C =',1P G9.2) WRITE(NOUT,113) 113 FORMAT('Parameters for contour deformation:') WRITE(NOUT,114)DFRM1,DFRM2,DFRM3,DFRM4 114 FORMAT(' Maximum of C:',1P G9.2,/, > ' Small P scale:',1P G9.2,/, > ' Medium P scale:',1P G9.2,/, > ' Large P scale:',1P G9.2) WRITE(NOUT,115) 115 FORMAT('Parameters for smearing in RTS:') WRITE(NOUT,116)SMEARFCTR,LOWPWR,HIGHPWR 116 FORMAT(' Scale factor for RTS/E0:',1P G9.2,/, > ' Small RTS power: ',I3,/, > ' Large RTS power: ',I3) WRITE(NOUT,117) 117 FORMAT('Renormalization parameter:') WRITE(NOUT,118)MUOVERRTS 118 FORMAT(' mu /Sqrt(s):',1P G9.2) WRITE(NOUT,119) 119 FORMAT('Cutoff parameters:') WRITE(NOUT,120)BADNESSLIMIT,CANCELLIMIT,THRUSTCUT 120 FORMAT(' Badness limit:',1P G9.2,/, > ' Cancellation limit:',1P G9.2,/, > ' Cut on 1 - thrust:',1P G9.2,/) WRITE(NOUT,121) 121 FORMAT('Graphs included: 1234567890') WRITE(NOUT,122)USEGRAPH 122 FORMAT(' ',10L1) C WRITE(NOUT,*)' ' C CALL RENO( > SUMR,ERRORR,SUMI,ERRORI, > SUMBISR,ERRORBISR,SUMBISI,ERRORBISI, > SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI, > INCLUDED,EXTRACOUNT,OMITTED, > NVALUE,VALUEMAX,KBAD,BADGRAPHNUMBER) C WRITE(NOUT,600) 600 FORMAT('Results for average of (1 - thrust)**2:') WRITE(NOUT,*)' ' WRITE(NOUT,800)INCLUDED WRITE(NOUT,801)EXTRACOUNT WRITE(NOUT,802)OMITTED 800 FORMAT(' Points included in result:',I8) 801 FORMAT(' Points in extra contribution:',I8) 802 FORMAT(' Points dropped entirely:',I8) WRITE(NOUT,*)' ' C WRITE(NOUT,601)SUMR,ERRORR 601 FORMAT(' Re(Result) =',1P G15.5,' +/-',1P G10.2) WRITE(NOUT,602)SUMI,ERRORI 602 FORMAT(' Im(Result) =',1P G15.5,' +/-',1P G10.2) WRITE(NOUT,603)SUMBISR,ERRORBISR 603 FORMAT(' Re(Extra Result) =',1P G15.5,' +/-',1P G10.2) WRITE(NOUT,604)SUMBISI,ERRORBISI 604 FORMAT(' Im(Extra Result) =',1P G15.5,' +/-',1P G10.2) WRITE(NOUT,605)SUMCHKR,ERRORCHKR 605 FORMAT(' Re(Check) =',1P G15.5,' +/-',1P G10.2) WRITE(NOUT,606)SUMCHKI,ERRORCHKI 606 FORMAT(' Im(Check) =',1P G15.5,' +/-',1P G10.2) C WRITE(NOUT,*)' ' C C Call Hrothgar to tell him that we are done with the calculation. C CALL HROTHGAR(1,KCUT0,1.0D0,'CALCDONE ') C WRITE(NOUT,*)' ' WRITE(NOUT,*)' ' WRITE(NOUT,700) 700 FORMAT('***********************') WRITE(NOUT,701) 701 FORMAT('Diagnostic information: ') WRITE(NOUT,*)' ' C DO N = -9,6 WRITE(NOUT,702)N,N+1,NVALUE(N) 702 FORMAT('Number of points with ',I2,' < log_10(|v|) <', I2, > ' is ',I8) ENDDO C WRITE(NOUT,*)' ' WRITE(NOUT,703)VALUEMAX 703 FORMAT('Biggest contribution was',1P G12.3) CALL DIAGNOSTIC(KBAD,BADGRAPHNUMBER) C C---------------------------------------------- C CALL DAYTIME(DATE) WRITE(NOUT,103)DATE 103 FORMAT(/ > ,' DONE ',A,/) C C---Timing function for Sun workstations; remove for others C CPUTIME = DTIME(TIMEARRAY) WRITE(NOUT,*)' CPUTIME = ',CPUTIME C--- C STOP END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE HROTHGAR(NCUT,KCUT,WEIGHT,WHATTODO) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER NCUT REAL*8 KCUT(SIZE+1,0:3) REAL*8 WEIGHT CHARACTER*10 WHATTODO C C Accumulates the results. Action depends on value of WHATTODO. C C 'STARTCALC ': Initialize variables at start of run. C 'STARTGROUP': Initialize variables at start of new group of points. C 'STARTPOINT': Initialize variables at start of a new point. C 'BADPOINT ': The point is too near a singularity. C Set switch BADPOINTQ to true. C 'XTRAPOINT ': The point involves a lot of cancellation, but not C so much as to trigger a BADPOINT signal. C Set switch XTRAPOINTQ to true. C 'NEWRESULT ': We have final state momentum variables NCUT, KCUT C and corresponding WEIGHT for a given final state cut. C We compute the corresponding INTEGRAND = WEIGHT*CALSVAL C and add it to VALUE for the current point. In Hrothgar, C WEIGHT and VALUE are real. C 'POINTDONE ': We now have the accumulated VALUE for this point. C *If we we earlier set BADPOINTQ to true, then we C simply throw away this point. If GOODPOINTQ is true C we proceed according to whether there was too much C cancellation or not: C *If XTRAPOINTQ is false we are ok and C we add VALUE to INTEGAL C *If XTRAPOINTQ is true C we add VALUE to INTEGALBIS C 'GROUPDONE ': We increment the sums for this group. C 'CALCDONE ': Compute the results and print them. C C 26 February 1998 C 9 March 1998 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT INTEGER GROUPSIZE,NRENO COMMON /MONTECARLO/GROUPSIZE,NRENO C SAVE SUM,SUMBIS,SQRSUM,SQRSUMBIS,INTEGRAL,INTEGRALBIS, > VALUE,MAXPART,BADPOINTQ,XTRAPOINTQ C REAL*8 SUM(30),SUMBIS(30),SQRSUM(30),SQRSUMBIS(30) REAL*8 INTEGRAL(30),INTEGRALBIS(30) REAL*8 VALUE(30),MAXPART(30) LOGICAL BADPOINTQ,XTRAPOINTQ REAL*8 CALSVALUE,CALS REAL*8 ERROR,ERRORBIS INTEGER NBINS PARAMETER (NBINS = 9) INTEGER N REAL*8 KN C IF (WHATTODO .EQ. 'STARTCALC ') THEN DO N = 1,NBINS SUM(N) = 0.0D0 SUMBIS(N) = 0.0D0 SQRSUM(N) = 0.0D0 SQRSUMBIS(N) = 0.0D0 ENDDO ELSEIF (WHATTODO .EQ. 'STARTGROUP') THEN DO N = 1,NBINS INTEGRAL(N) = 0.0D0 INTEGRALBIS(N) = 0.0D0 ENDDO ELSEIF (WHATTODO .EQ. 'STARTPOINT') THEN DO N = 1,NBINS VALUE(N) = 0.0D0 ENDDO BADPOINTQ = .FALSE. XTRAPOINTQ = .FALSE. ELSEIF (WHATTODO .EQ. 'BADPOINT ') THEN BADPOINTQ = .TRUE. ELSEIF (WHATTODO .EQ. 'XTRAPOINT ') THEN XTRAPOINTQ = .TRUE. ELSEIF (WHATTODO .EQ. 'NEWRESULT ') THEN DO N = 1,NBINS CALSVALUE = CALS(KCUT,NCUT,N) VALUE(N) = VALUE(N) + WEIGHT*CALSVALUE ENDDO ELSEIF (WHATTODO .EQ. 'POINTDONE ') THEN IF (.NOT.BADPOINTQ) THEN IF (XTRAPOINTQ) THEN DO N = 1,NBINS INTEGRALBIS(N) = INTEGRALBIS(N) + VALUE(N) ENDDO ELSE DO N = 1,NBINS INTEGRAL(N) = INTEGRAL(N) + VALUE(N) ENDDO ENDIF ENDIF ELSEIF (WHATTODO .EQ. 'GROUPDONE ') THEN DO N = 1,NBINS INTEGRAL(N) = INTEGRAL(N)/GROUPSIZE INTEGRALBIS(N) = INTEGRALBIS(N)/GROUPSIZE SUM(N) = SUM(N)+ INTEGRAL(N) SUMBIS(N) = SUMBIS(N) + INTEGRALBIS(N) SQRSUM(N) = SQRSUM(N) + INTEGRAL(N)**2 SQRSUMBIS(N) = SQRSUMBIS(N) + INTEGRALBIS(N)**2 ENDDO ELSEIF (WHATTODO .EQ. 'CALCDONE ') THEN WRITE(NOUT,001) DO N = 1,NBINS SUM(N) = SUM(N)/NRENO SUMBIS(N) = SUMBIS(N)/NRENO SQRSUM(N) = SQRSUM(N)/NRENO SQRSUMBIS(N) = SQRSUMBIS(N)/NRENO ERROR = SQRT((SQRSUM(N) - SUM(N)**2)/(NRENO - 1)) ERRORBIS = SQRT((SQRSUMBIS(N) - SUMBIS(N)**2)/(NRENO - 1)) WRITE(NOUT,002)(0.68D0 + N * 0.03D0),KN(0.68D0 + N * 0.03D0) WRITE(NOUT,003)SUM(N),ERROR WRITE(NOUT,004)SUMBIS(N),ERRORBIS ENDDO ENDIF RETURN 001 FORMAT(/,'Hrothgar reports',/) 002 FORMAT('T =',1P G10.3,' Kunszt-Nason function =',1P G10.2) 003 FORMAT(' Result =',1P G15.5,' +/-',1P G10.2) 004 FORMAT(' Added Result =',1P G15.5,' +/-',1P G10.2,/) END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION CALS(KCUT,NCUT,N) C INTEGER SIZE PARAMETER (SIZE = 3) C In: REAL*8 KCUT(SIZE+1,0:3) INTEGER NCUT,N C C Average value of 1 - thrust in bins. C C 9 March 1998 C REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT C REAL*8 RTS REAL*8 THRUST,T1,T2,T REAL*8 RHO,KN INTEGER I REAL*8 DELTATHRUST,MINTHRUST PARAMETER (DELTATHRUST = 0.03D0) PARAMETER (MINTHRUST = 0.68D0) C CALS = 0.0D0 C RTS = 0.0D0 DO I = 1,NCUT RTS = RTS + KCUT(I,0) ENDDO C T1 = MINTHRUST + (N-1) * DELTATHRUST T2 = MINTHRUST + (N+1) * DELTATHRUST T = THRUST(NCUT,KCUT,RTS) IF ( (T.GT.T1).AND.(T.LT.T2) ) THEN RHO = 15.0D0/(16.0D0*DELTATHRUST**5) * ((T-T1)*(T2-T))**2 CALS = (1.0D0 - T)*RHO/KN(T) ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION KN(T) C In: REAL*8 T C REAL*8 MUOVERRTS COMMON /RENORMALIZE/ MUOVERRTS REAL*8 NF PARAMETER (NF = 5.0D0) REAL*8 L,NONSING,SING,LO C L = LOG(1.0D0/(1.0D0 - T)) C NONSING = -16100.550195D0 * (0.46789210556 + T) > * (1.1184765748D0 - 2.1058049743D0 * T + T**2) > * (0.4925129711D0 - 1.3620871313D0 * T + T**2) C SING = -14.222222222D0 * L * (-5.098072232D0 + L ) > *( 0.691822232D0 + L) C LO = 2.0D0 * (3.0D0 * T**2 - 3.0D0 * T + 2.0D0) /T LO = LO * LOG( (2.0D0 * T - 1.0D0)/(1.0D0 - T) ) LO = LO - 3.0D0 * (3.0D0 * T - 2.0D0) * (2.0D0 - T) LO = LO * 4.0D0/3.0D0 C KN = SING + NONSING > + (11.0D0 - 2.0D0*NF/3.0D0)*LOG(MUOVERRTS)*LO C C We divide by 4 because KN report the coefficient of (alpha_s/2Pi)^2 C while I want the coefficient of (alpha_s/Pi)^2. C KN = KN/4.0D0 C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION CALS0(KCUT,NCUT) C INTEGER SIZE PARAMETER (SIZE = 3) C In: REAL*8 KCUT(SIZE+1,0:3) INTEGER NCUT C C Average value of (1 - thrust)**2. C C 10 March 1998 C REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT C REAL*8 RTS REAL*8 THRUST,ONEMINUSTHRUST INTEGER I C RTS = 0.0D0 DO I = 1,NCUT RTS = RTS + KCUT(I,0) ENDDO C ONEMINUSTHRUST = 1.0D0 - THRUST(NCUT,KCUT,RTS) C CALS0 = ONEMINUSTHRUST**2 C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION THRUST(NCUT,KCUT,RTS) C INTEGER SIZE PARAMETER (SIZE = 3) C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT c INTEGER NCUT REAL*8 KCUT(SIZE+1,0:3),RTS C C Calculates the value of the thrust. We calculate thrust as C 2.0D0*EJETMAX/RTS where EJETMAX is the maximum over subsets of the C final state particles of the absolute value of the sum of the momenta C of these particles. For three particles in the final state, EJETMAX is C the largest of the absolute values of the momenta of the three C particles. For four particles in the final state, EJETMAX is either C the largest of the absolute values of the momenta of the four C particles or else the largest of the absolute values of the C sum of the momenta of a pair of particles. Since the complementary C pair of paraticles has just the opposite momentum, we choose the pair C to contain the particle with the largest momentum. C C 6 March 1988 C INTEGER J,JMAX,MU REAL*8 EMAX,EMAXSQ,EJETSQ,EJETSQMAX,EJETMAX,DOT C EMAX = 0.0D0 DO J = 1,NCUT IF (KCUT(J,0).GT.EMAX) THEN EMAX = KCUT(J,0) JMAX = J ENDIF ENDDO C C What to do for NCUT = 3. C IF (NCUT.LE.3) THEN THRUST = 2.0D0*EMAX/RTS RETURN ELSEIF (NCUT.GT.4) THEN WRITE(NOUT,*)'Oops, there is a problem in THRUST' STOP ENDIF C C What to do for NCUT = 4. C EMAXSQ = EMAX**2 EJETSQMAX = EMAXSQ DO J = 1,NCUT IF (J.NE.JMAX) THEN DOT = 0.0D0 DO MU = 1,3 DOT = DOT + KCUT(J,MU)*KCUT(JMAX,MU) ENDDO EJETSQ = EMAXSQ + KCUT(J,0)**2 + 2.0D0*DOT IF (EJETSQ.GT.EJETSQMAX) THEN EJETSQMAX = EJETSQ ENDIF ENDIF ENDDO C EJETMAX = SQRT(EJETSQMAX) THRUST = 2.0D0*EJETMAX/RTS RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012