C23456789012345678901234567890123456789012345678901234567890123456789012 C C ------------------------ C beowulf.f Version 1.1.1 C ------------------------ C C Program for numerical integration of jet cross sections in C electron-positron annihilation. -- D. E. Soper C C First code: 29 November 1992. C Latest revision: see revision date below. C C The main program and subroutines that a user might want to modify C are contained in this package, beowulf.f. In particular, a C user may very well want to modify parameter settings in the main C program and to change the observables calculated in the subroutine C HROTHGAR and in the functions CALStype(NCUT,KCUT,index). Subroutines C that can be modified only at extreme peril to the reliability of C the results are in the companion package beowulfsubs.f. C PROGRAM BEOWULF C CHARACTER*36 REVISIONDATE PARAMETER(REVISIONDATE = C ------------------------------------ >'Latest revision 16 August 2001 ') C ------------------------------------ C Array sizes: INTEGER SIZE,MAXGRAPHS,MAXCUTS,MAXMAPS PARAMETER (SIZE = 3) PARAMETER (MAXGRAPHS = 10) PARAMETER (MAXCUTS = 9) PARAMETER (MAXMAPS = 64) C Input and output units: INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT C Number of groups of points and which graphs to use for RENO: INTEGER NSETS INTEGER COUNTFACTOR(MAXGRAPHS,MAXMAPS) INTEGER GROUPSIZE(MAXGRAPHS,MAXMAPS) INTEGER GROUPSIZEGRAPH(MAXGRAPHS) INTEGER GROUPSIZETOTAL COMMON /MONTECARLO/GROUPSIZE,GROUPSIZEGRAPH,GROUPSIZETOTAL LOGICAL USEGRAPH(MAXGRAPHS) COMMON /WHICHGRAPHS/ USEGRAPH C Global size and scale parameters: INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX REAL*8 ENERGYSCALE COMMON /MSCALE/ ENERGYSCALE C Color factors: REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C Information on cut structure: INTEGER NCUTINFO(MAXGRAPHS,MAXCUTS) INTEGER ISIGNINFO(MAXGRAPHS,MAXCUTS,3*SIZE + 1) INTEGER CUTINDEXINFO(MAXGRAPHS,MAXCUTS,SIZE + 1) INTEGER CUTSIGNINFO(MAXGRAPHS,MAXCUTS,SIZE + 1) LOGICAL LEFTLOOPINFO(MAXGRAPHS,MAXCUTS) LOGICAL RIGHTLOOPINFO(MAXGRAPHS,MAXCUTS) INTEGER NINLOOPINFO(MAXGRAPHS,MAXCUTS) INTEGER LOOPINDEXINFO(MAXGRAPHS,MAXCUTS,SIZE+1) INTEGER LOOPSIGNINFO(MAXGRAPHS,MAXCUTS,SIZE+1) COMMON /CUTINFORMATION/ NCUTINFO,ISIGNINFO, > CUTINDEXINFO,CUTSIGNINFO,LEFTLOOPINFO,RIGHTLOOPINFO, > NINLOOPINFO,LOOPINDEXINFO,LOOPSIGNINFO C How many graphs and how many cuts and maps for each: INTEGER NUMBEROFGRAPHS INTEGER NUMBEROFCUTS(MAXGRAPHS) INTEGER NUMBEROFMAPS(MAXGRAPHS) COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS C Limits: REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT C Timing limit: REAL*8 HOURLIMIT REAL*8 TIMELIMIT COMMON /MAXTIME/ TIMELIMIT C Deform parameters: REAL*8 DEFORMALPHA,DEFORMBETA,DEFORMGAMMA COMMON /DEFORMSCALES/DEFORMALPHA,DEFORMBETA,DEFORMGAMMA C Smear parameters: REAL*8 SMEARFCTR INTEGER LOWPWR,HIGHPWR COMMON /SMEARPARMS/ SMEARFCTR,LOWPWR,HIGHPWR C Renormaliation parameters: REAL*8 MUOVERRTS COMMON /RENORMALIZE/ MUOVERRTS C Diagnostic data instructions: LOGICAL REPORT,DETAILS COMMON /CALCULOOK/ REPORT,DETAILS C RENO results: REAL*8 SUMR,ERRORR,SUMI,ERRORI REAL*8 SUMBIS,ERRORBIS REAL*8 SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI REAL*8 FLUCT(MAXGRAPHS,MAXMAPS) INTEGER*8 INCLUDED,EXTRACOUNT,OMITTED INTEGER NVALPT(-9:6) REAL*8 VALPTMAX REAL*8 KBAD(0:3*SIZE-1,0:3) INTEGER BADGRAPHNUMBER,BADMAPNUMBER INTEGER NRENO REAL*8 CPUTIME C Hrothgar dummy variables: REAL*8 KCUT0(SIZE+1,0:3) C Seed for random number generator. INTEGER SEED C Today's date: CHARACTER DATE*(28) C Loop control variables: INTEGER N,MU,GRAPHNUMBER,MAPNUMBER C Logical variables that control aspects of the output: LOGICAL ALLGRAPHS,WRITEPROGRAMINFO,WRITEFLUCT,WRITEDIAGNOSTIC C A temporary real number: REAL*8 TEMP C A temporary integer: INTEGER ITEMP C Average of FLUCT variables: REAL*8 AVFLUCT C C----------------------------------------------------------------------- C C Input and ouput units. NIN = 5 NOUT = 6 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 Color factors. NC = 3.0D0 NF = 5.0D0 C Energy scale. ENERGYSCALE = 1.0D0 C Initialize common /CUTINFORMATION/ and /GRAPHCOUNTS/. CALL MAKECUTINFO C Parameters for subroutine CALCULATE for throwing away bad points. BADNESSLIMIT = 1.0D4 CANCELLIMIT = 1.0D4 C Cutoff on 1 - THRUST. In the this version of beowulf.f, THRUSTCUT C has no effect. However, it is still included in the common block C /LIMITS/ since it is sometimes useful to use a cutoff that prevents C points that produce a final state with thrust very near 1 from being C accepted. THRUSTCUT = 0.0D0 C Parameters for subroutine DEFORM. DEFORMALPHA = 0.7D0 DEFORMBETA = 1.0D0 DEFORMGAMMA = 1.0D0 C Parmameters for function SMEAR. SMEARFCTR = 3.0D0 LOWPWR = 9 HIGHPWR = 18 C Renormalization parameters. MUOVERRTS = 0.3D0 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 C Here we set the count factors COUNTFACTOR(GRAPHNUMBER,MAPNUMBER). C These, multiplied by NSETS, give GROUPSIZE(GRAPHNUMBER,MAPNUMBER). C CALL GETCOUNTS(COUNTFACTOR) C C Call Hrothgar to tell him to get ready. C Here KCUT0 are just dummy variables for Hrothgar. C DO N = 1,SIZE+1 DO MU = 0,3 KCUT0(N,MU) = 1.0D0 ENDDO ENDDO CALL HROTHGAR(1,KCUT0,1.0D0,1,'STARTCALC ') C--- WRITE(NOUT,*)'Please give the approximate CPU time limit (hours).' READ(NIN,*) HOURLIMIT TIMELIMIT = HOURLIMIT * 3600 IF(HOURLIMIT.LT.1.0D0) THEN NSETS = 1 + NINT(10.0D0*HOURLIMIT) ELSE NSETS = 11 + NINT(HOURLIMIT) ENDIF C--- C To change the default value calculated above, use this. C C WRITE(NOUT,*)'Please give number of sets of points in each group' C READ(NIN,*) NSETS C--- GROUPSIZETOTAL = 0 DO GRAPHNUMBER = 1,MAXGRAPHS GROUPSIZEGRAPH(GRAPHNUMBER) = 0 IF (USEGRAPH(GRAPHNUMBER)) THEN DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER) ITEMP = NSETS * COUNTFACTOR(GRAPHNUMBER,MAPNUMBER) GROUPSIZE(GRAPHNUMBER,MAPNUMBER) = ITEMP GROUPSIZEGRAPH(GRAPHNUMBER) = GROUPSIZEGRAPH(GRAPHNUMBER)+ITEMP GROUPSIZETOTAL = GROUPSIZETOTAL + ITEMP ENDDO ENDIF ENDDO C-- C If you want to extract parts of the results with different color C and number-of-flavor factors, you may want to be able to C adjust the number of colors and flavors. But note that the C function KN(T) that contains the "standard" results C for the thrust distribution is for three colors and five flavors. C C WRITE(NOUT,*)'Please give the numbers of colors and flavors.' C READ(NIN,*)NC,NF C-- C This overrides the default value MUOVERRTS = 0.3 that was set above. C WRITE(NOUT,*) > 'Please give ratio of the renormalization scale to sqrt(s).' READ(NIN,*) MUOVERRTS C WRITE(NOUT,*)'Please give a seed (0 ,' beowulf version 1.1.1 ',A,/,50('-')) C C**************************************************** WRITE(NOUT,*) REVISIONDATE 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,101) HOURLIMIT,NSETS 101 FORMAT('Beowulf will work for',F8.2,' hours',/, > 'using groups of',I4,' sets of points.') WRITE(NOUT,102) SEED 102 FORMAT('The seed is',I8,'.') WRITE(NOUT,103)NC,NF 103 FORMAT('Use',1P G9.2,' colors and',1P G9.2,' flavors.') WRITE(NOUT,104) 104 FORMAT('Renormalization parameter:') WRITE(NOUT,105)MUOVERRTS 105 FORMAT(' mu /Sqrt(s):',1P G9.2) WRITE(NOUT,106) 106 FORMAT('Cutoff parameters:') WRITE(NOUT,107)BADNESSLIMIT,CANCELLIMIT 107 FORMAT(' Badness limit:',1P G9.2,/, > ' Cancellation limit:',1P G9.2,/) C ALLGRAPHS = .TRUE. DO N = 1,MAXGRAPHS IF (.NOT.USEGRAPH(N) ) THEN ALLGRAPHS = .FALSE. ENDIF ENDDO IF (ALLGRAPHS) THEN WRITE(NOUT,108) 108 FORMAT('All graphs are used.') ELSE DO N = 1,MAXGRAPHS WRITE(NOUT,109) N,USEGRAPH(N) 109 FORMAT('Use of graph ',I2,': ',L1) ENDDO ENDIF C WRITEPROGRAMINFO = .FALSE. IF (WRITEPROGRAMINFO) THEN WRITE(NOUT,123) 123 FORMAT(/,'Points per group:',/) DO GRAPHNUMBER = 1,NUMBEROFGRAPHS IF (USEGRAPH(GRAPHNUMBER)) THEN DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER) WRITE(NOUT,124) GRAPHNUMBER,MAPNUMBER, > GROUPSIZE(GRAPHNUMBER,MAPNUMBER) 124 FORMAT('Graph',I3,' Map',I3,' Points used',I6) ENDDO WRITE(NOUT,125)GROUPSIZEGRAPH(GRAPHNUMBER) 125 FORMAT(/,'Total points for this graph',I8,/) ENDIF ENDDO WRITE(NOUT,126)GROUPSIZETOTAL 126 FORMAT('Total points in each group:',I8,/) WRITE(NOUT,127) 127 FORMAT('Parameters for contour deformation:') WRITE(NOUT,128)DEFORMALPHA,DEFORMBETA,DEFORMGAMMA 128 FORMAT(' alpha:',1P G9.2,/, > ' beta :',1P G9.2,/, > ' gamma:',1P G9.2) WRITE(NOUT,129) 129 FORMAT('Parameters for smearing in RTS:') WRITE(NOUT,130)SMEARFCTR,LOWPWR,HIGHPWR 130 FORMAT(' Scale factor for RTS/E0:',1P G9.2,/, > ' Small RTS power: ',I3,/, > ' Large RTS power: ',I3) ENDIF C WRITE(NOUT,*)' ' C CALL RENO( > SUMR,ERRORR,SUMI,ERRORI, > SUMBIS,ERRORBIS, > SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI,FLUCT, > INCLUDED,EXTRACOUNT,OMITTED, > NVALPT,VALPTMAX,KBAD,BADGRAPHNUMBER,BADMAPNUMBER, > NRENO,CPUTIME) C WRITE(NOUT,598) CPUTIME 598 FORMAT('After ',F10.1,' CPU seconds, Beowulf is done.') WRITE(NOUT,599) NRENO 599 FORMAT('Beowulf used',I7,' groups of points.',/) 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:',I12) 801 FORMAT(' Points in cutoff correction:',I12) 802 FORMAT(' Points dropped entirely:',I12) WRITE(NOUT,*)' ' C WRITE(NOUT,601)SUMR,ERRORR 601 FORMAT(' Re(Result) =',1P G13.5,' +/-',1P G10.2) WRITE(NOUT,602)SUMI,ERRORI 602 FORMAT(' Im(Result) =',1P G13.5,' +/-',1P G10.2) WRITE(NOUT,603)SUMBIS,ERRORBIS 603 FORMAT(' Cutoff correction =',1P G13.5,' +/-',1P G10.2) C C Check integral. C WRITE(NOUT,604) 604 FORMAT(/,'Check integral, equal to 1.0 plus a cutoff error.',/, > 'For Badness limit = Cancellation limit = 1.00E+04,',/, > 'the check integral is approximately 0.95.') WRITE(NOUT,605)SUMCHKR,ERRORCHKR 605 FORMAT(' Re(Check) =',1P G13.5,' +/-',1P G10.2) WRITE(NOUT,606)SUMCHKI,ERRORCHKI 606 FORMAT(' Im(Check) =',1P G13.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,NRENO,'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,NVALPT(N) 702 FORMAT('Number of points with ',I2,' < log_10(|v|) <', I2, > ' is ',I10) ENDDO C WRITE(NOUT,*)' ' WRITE(NOUT,703)VALPTMAX 703 FORMAT('Biggest contribution was',1P G12.3) WRITE(NOUT,704)BADGRAPHNUMBER,BADMAPNUMBER 704 FORMAT('From graph',I3,', map',I3) C C Setting WRITEFLUCT to .TRUE. will result in the generation of C information where the fluctions in SUMR come from, including C a suggested set of COUNTFACTOR values that should make the C fluctions in SUMR smaller. C WRITEFLUCT = .FALSE. IF (WRITEFLUCT) THEN AVFLUCT = 0.0D0 DO GRAPHNUMBER = 1,NUMBEROFGRAPHS DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER) AVFLUCT = AVFLUCT > + FLUCT(GRAPHNUMBER,MAPNUMBER) > * GROUPSIZE(GRAPHNUMBER,MAPNUMBER)/GROUPSIZETOTAL ENDDO ENDDO WRITE(NOUT,*)' ' WRITE(NOUT,*)'Information on fluctuations' DO GRAPHNUMBER = 1,NUMBEROFGRAPHS DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER) WRITE(NOUT,705) GRAPHNUMBER,MAPNUMBER, > GROUPSIZE(GRAPHNUMBER,MAPNUMBER), > FLUCT(GRAPHNUMBER,MAPNUMBER)/AVFLUCT 705 FORMAT('Graph',I3,' Map',I3,' Points used',I6, > ' Fluctuation',G10.2) ENDDO ENDDO WRITE(NOUT,*)' ' WRITE(NOUT,*)'Recommended new point distribution' DO GRAPHNUMBER = 1,NUMBEROFGRAPHS DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER) TEMP = FLUCT(GRAPHNUMBER,MAPNUMBER)/AVFLUCT TEMP = SQRT(TEMP) * 0.95 TEMP = COUNTFACTOR(GRAPHNUMBER,MAPNUMBER)*TEMP ITEMP = NINT(TEMP) IF (ITEMP.EQ.0) THEN ITEMP = 1 ENDIF WRITE(NOUT,706) GRAPHNUMBER,MAPNUMBER,ITEMP 706 FORMAT(' COUNTFACTOR(',I2,',',I2,') = ',I4) ENDDO ENDDO WRITE(NOUT,*)' ' ENDIF C WRITEDIAGNOSTIC = .TRUE. IF (WRITEDIAGNOSTIC) THEN CALL DIAGNOSTIC(KBAD,BADGRAPHNUMBER) ENDIF C C---------------------------------------------- C CALL DAYTIME(DATE) WRITE(NOUT,707)DATE 707 FORMAT(/ > ,' DONE ',A,/) C C--- C STOP END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE HROTHGAR(NCUT,KCUT,WEIGHT,NRENO,WHATTODO) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER NCUT REAL*8 KCUT(SIZE+1,0:3) REAL*8 WEIGHT INTEGER NRENO CHARACTER*10 WHATTODO C C Accumulates the results. C C The results are stored as they accumulate in arrays of size NBINSTOT. C Each array element corresponds to a different observable -- for C example a certain moment of the thrust distribution. C The arrays are divided into subsets of size NBINS1, NBINS2, ... , C with NBINS1 + NBINS2 + ... = NBINSTOT. Each subset corresponds C to a different type of observable, such as "moments of the C thrust distribution." The calculations of the C observation function for each type of observable is done C by a function named CALStype(NCUT,KCUT,index), where the index C for type j takes values 1,...,NBINSj. Users can remove C observable types and can add new observable types by writing new C functions CALStype(NCUT,KCUT,index) and by modfifying the reporting C of the results, which is controlled by this subroutine. C C Note that, in addition to sending sets of final state momenta along C with weights to HROTHGAR, the other subroutines of beowulf compute C a single integral, which is defined by the observable in the C function CALS0(NCUT,KCUT) and is the average value of (1-thrust)^2 C in the default version of the beowulf.f subroutines. If there is C too much cancellation in the calculation of this observable among C the contributions from different final states corresponding to C a given point in the space of loop momenta, then HROTHGAR is sent C a BADPOINT or XTRAPOINT signal. C C The action of HROTHGAR 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 in the C space of loop momenta. C 'BADPOINT ': The point is too near a singularity C OR there is too much cancellation. C Set switch BADPOINTQ to true. C 'XTRAPOINT ': The point is too near a singularity C OR there is too much cancellation. C But not so bad 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 BADPOINTQ is false 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. This includes C sums for the observables being computed and also C sums to use for computing the statistical errors. C 'CALCDONE ': Compute the results and print them. C C The limits: C Cancel VALUE,MAXPART,BADPOINTQ,XTRAPOINTQ C INTEGER NBINS1,NBINS2,NBINS3,NBINSTOT PARAMETER (NBINS1=9,NBINS2=10,NBINS3=10,NBINSTOT = 29) REAL*8 SUM(NBINSTOT),SUMBIS(NBINSTOT) REAL*8 SUM2(NBINSTOT),SUM2BIS(NBINSTOT) REAL*8 SUM3(NBINSTOT),SUM4(NBINSTOT) REAL*8 INTEGRAL(NBINSTOT),INTEGRALBIS(NBINSTOT) REAL*8 VALUE(NBINSTOT),MAXPART(NBINSTOT) REAL*8 ERROR(NBINSTOT),ERROR4(NBINSTOT),ERRORBIS(NBINSTOT) C LOGICAL BADPOINTQ,XTRAPOINTQ REAL*8 CALSVALUE,CALSTHRUST,CALSTMOMENTS,CALSYMOMENTS INTEGER N,NTEMP REAL*8 KN C IF (WHATTODO .EQ. 'STARTCALC ') THEN DO N = 1,NBINSTOT SUM(N) = 0.0D0 SUMBIS(N) = 0.0D0 SUM2(N) = 0.0D0 SUM2BIS(N) = 0.0D0 SUM3(N) = 0.0D0 SUM4(N) = 0.0D0 ENDDO ELSEIF (WHATTODO .EQ. 'STARTGROUP') THEN DO N = 1,NBINSTOT INTEGRAL(N) = 0.0D0 INTEGRALBIS(N) = 0.0D0 ENDDO ELSEIF (WHATTODO .EQ. 'STARTPOINT') THEN DO N = 1,NBINSTOT VALUE(N) = 0.0D0 ENDDO BADPOINTQ = .FALSE. XTRAPOINTQ = .FALSE. ELSEIF (WHATTODO .EQ. 'BADPOINT ') THEN BADPOINTQ = .TRUE. ELSEIF (WHATTODO .EQ. 'XTRAPOINT ') THEN XTRAPOINTQ = .TRUE. C--------------------------------------------------------------- C Here is where we calculate the functions CALS for each type of C observable using functions CALStype(NCUT,KCUT,index). C ELSEIF (WHATTODO .EQ. 'NEWRESULT ') THEN N = 0 DO NTEMP = 1,NBINS1 N = N+1 CALSVALUE = CALSTHRUST(NCUT,KCUT,NTEMP) VALUE(N) = VALUE(N) + WEIGHT*CALSVALUE ENDDO DO NTEMP = 1,NBINS2 N = N+1 CALSVALUE = CALSTMOMENTS(NCUT,KCUT,NTEMP) VALUE(N) = VALUE(N) + WEIGHT*CALSVALUE ENDDO DO NTEMP = 1,NBINS3 N = N+1 CALSVALUE = CALSYMOMENTS(NCUT,KCUT,NTEMP) VALUE(N) = VALUE(N) + WEIGHT*CALSVALUE ENDDO C--------------------------------------------------------------- ELSEIF (WHATTODO .EQ. 'POINTDONE ') THEN IF (.NOT.BADPOINTQ) THEN IF (XTRAPOINTQ) THEN DO N = 1,NBINSTOT INTEGRALBIS(N) = INTEGRALBIS(N) + VALUE(N) ENDDO ELSE DO N = 1,NBINSTOT INTEGRAL(N) = INTEGRAL(N) + VALUE(N) ENDDO ENDIF ENDIF ELSEIF (WHATTODO .EQ. 'GROUPDONE ') THEN DO N = 1,NBINSTOT SUM(N) = SUM(N) + INTEGRAL(N) SUMBIS(N) = SUMBIS(N) + INTEGRALBIS(N) SUM2(N) = SUM2(N) + INTEGRAL(N)**2 SUM2BIS(N) = SUM2BIS(N) + INTEGRALBIS(N)**2 SUM3(N) = SUM3(N) + INTEGRAL(N)**3 SUM4(N) = SUM4(N) + INTEGRAL(N)**4 ENDDO ELSEIF (WHATTODO .EQ. 'CALCDONE ') THEN WRITE(NOUT,001) DO N = 1,NBINSTOT SUM(N) = SUM(N)/NRENO SUMBIS(N) = SUMBIS(N)/NRENO SUM2(N) = SUM2(N)/NRENO SUM2BIS(N) = SUM2BIS(N)/NRENO ERROR(N) = SQRT((SUM2(N) - SUM(N)**2)/(NRENO - 1)) ERRORBIS(N) = SQRT((SUM2BIS(N) - SUMBIS(N)**2)/(NRENO - 1)) SUM3(N) = SUM3(N)/NRENO SUM4(N) = SUM4(N)/NRENO ERROR4(N) = (SUM4(N) - 4*SUM3(N)*SUM(N) + 6*SUM2(N)*SUM(N)**2 > - 3*SUM(N)**4)/(3.0D0*(NRENO-1)**2) ERROR4(N) = ABS(ERROR4(N))**0.25 ENDDO C------------------------------------------------------- C Write the labels for each subset of our results arrays C and write the results themselves. C N = 0 DO NTEMP = 1,NBINS1 N = N+1 WRITE(NOUT,002)(0.68D0 + NTEMP*0.03D0), > KN(0.68D0 + NTEMP*0.03D0) WRITE(NOUT,003)SUM(N),ERROR(N) WRITE(NOUT,004)SUMBIS(N),ERRORBIS(N) WRITE(NOUT,005)ERROR4(N) ENDDO WRITE(NOUT,006) DO NTEMP = 1,NBINS2 N = N+1 WRITE(NOUT,007)1.0D0+NTEMP*0.5D0 WRITE(NOUT,003)SUM(N),ERROR(N) WRITE(NOUT,004)SUMBIS(N),ERRORBIS(N) WRITE(NOUT,005)ERROR4(N) ENDDO WRITE(NOUT,008) DO NTEMP = 1,NBINS3 N = N+1 WRITE(NOUT,007)1.0D0+NTEMP*0.5D0 WRITE(NOUT,003)SUM(N),ERROR(N) WRITE(NOUT,004)SUMBIS(N),ERRORBIS(N) WRITE(NOUT,005)ERROR4(N) ENDDO C------------------------------------------------------- ENDIF RETURN C C FORMAT statements for the output. C 001 FORMAT(/, ' --- Hrothgar reports ---',/,/, > 'First, the order alpha_s^2 contribution to the',/, > 'thrust distribution d sigma /d T divided by the',/, > 'same quantity as reported in Kunszt and Nason.'/) 002 FORMAT('T =',F6.3,' Kunszt-Nason function =',1P G10.2) 003 FORMAT(' Result =',1P G13.5,' +/-',1P G10.2) 004 FORMAT(' Cutoff correction =',1P G13.5,' +/-',1P G10.2) 005 FORMAT(' Alternative error =',1P G10.2,/) 006 FORMAT(/,'Next, moments <(1-t)^n>.'/) 007 FORMAT('n =',F6.3) 008 FORMAT(/,'Next, moments <(y_cut)^n>.'/) END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION CALSTHRUST(NCUT,KCUT,N) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER NCUT,N REAL*8 KCUT(SIZE+1,0:3) C C Average value of 1 - thrust in bins. C C 9 March 1998 C 5 August 1999 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 CALSTHRUST = 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 CALSTHRUST = (1.0D0 - T)*RHO/KN(T) ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION CALSTMOMENTS(NCUT,KCUT,N) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER NCUT,N REAL*8 KCUT(SIZE+1,0:3) C C Moments of 1 - thrust. C 20 September 1999 C 29 September 1999 C REAL*8 RTS REAL*8 THRUST,T,POWER INTEGER I REAL*8 TINY PARAMETER (TINY = 1.0D-12) C CALSTMOMENTS = 0.0D0 C RTS = 0.0D0 DO I = 1,NCUT RTS = RTS + KCUT(I,0) ENDDO C T = THRUST(NCUT,KCUT,RTS) POWER = 1.0D0 + 0.5D0*N IF ( (1.0D0 - T) .GT. TINY ) THEN CALSTMOMENTS = (1.0D0 - T)**POWER ELSE CALSTMOMENTS = 0.0D0 ENDIF 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 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 C The coefficient of (alpha_s/Pi)^2 in the thrust distribution as C given by Kunszt and Nason. C Note that this function is only for NC = 3, NF = 5. C C 9 April 1998 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 CALSYMOMENTS(NCUT,KCUT,N) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER NCUT,N REAL*8 KCUT(SIZE+1,0:3) C C Moments of Yjet. C 20 September 1999 C 30 September 1999 C REAL*8 Y2JET,Y4JET,POWER REAL*8 TINY PARAMETER (TINY = 1.0D-12) C CALSYMOMENTS = 0.0D0 C CALL COMBINEJETS(NCUT,KCUT,Y4JET,Y2JET) C IF (Y2JET.GT.Y4JET) THEN POWER = 1.0D0 + 0.5D0*N CALSYMOMENTS = 0.0D0 IF (Y2JET.GT.TINY) THEN CALSYMOMENTS = CALSYMOMENTS + Y2JET**POWER ENDIF IF (Y4JET.GT.TINY) THEN CALSYMOMENTS = CALSYMOMENTS - Y4JET**POWER ENDIF ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE COMBINEJETS(NCUT,KCUT,Y4JET,Y2JET) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER NCUT REAL*8 KCUT(SIZE+1,0:3) C Out: REAL*8 Y4JET,Y2JET C C Find Y4jet and Y2jet for combining three or four partons to jets C according to the Durham algorithm. C C 4 August 1999 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT C REAL*8 RTS REAL*8 DOTIJ,COSIJ,YIJ REAL*8 KJET(SIZE,0:3),ABSKJET(SIZE) REAL*8 KJETSQ INTEGER I,J,ISTAR,JSTAR,MU LOGICAL OK C RTS = 0.0D0 DO I = 1,NCUT RTS = RTS + KCUT(I,0) ENDDO C IF (NCUT.EQ.4) THEN C C What to do for NCUT = 4. C C 1) Find the smallest Yij. C Y4JET = 1.0D0 OK = .FALSE. DO I = 2,4 DO J = 1,I-1 DOTIJ = 0.0D0 DO MU = 1,3 DOTIJ = DOTIJ + KCUT(I,MU)*KCUT(J,MU) ENDDO COSIJ = DOTIJ/KCUT(I,0)/KCUT(J,0) YIJ = 1.0D0 - COSIJ YIJ = YIJ * 2.0D0 * MIN(KCUT(I,0),KCUT(J,0))**2 YIJ = YIJ / RTS**2 IF (YIJ.LE.Y4JET) THEN Y4JET = YIJ ISTAR = I JSTAR = J OK = .TRUE. ENDIF ENDDO ENDDO IF (.NOT.OK) THEN WRITE(NOUT,*) 'Problem in jet combination.' STOP ENDIF C C 2) Combine the momenta. C KJET(1,0) = KCUT(ISTAR,0) + KCUT(JSTAR,0) KJETSQ = 0.0D0 DO MU = 1,3 KJET(1,MU) = KCUT(ISTAR,MU) + KCUT(JSTAR,MU) KJETSQ = KJETSQ + KJET(1,MU)**2 ENDDO ABSKJET(1) = SQRT(KJETSQ) C C 3) Make the momenta for the other protojets. C I = 1 DO J = 1,4 IF ((J.NE.ISTAR).AND.(J.NE.JSTAR)) THEN I = I + 1 DO MU = 0,3 KJET(I,MU) = KCUT(J,MU) ENDDO ABSKJET(I) = KJET(I,0) ENDIF ENDDO IF (I.NE.3) THEN WRITE(NOUT,*) 'Messed up jet combination.' STOP ENDIF C ELSEIF (NCUT.EQ.3) THEN C Y4JET = 0.0D0 DO I = 1,3 DO MU = 0,3 KJET(I,MU) = KCUT(I,MU) ENDDO ABSKJET(I) = KJET(I,0) ENDDO C C End NCUT cases. C ELSE WRITE(NOUT,*) 'NCUT not 3 or 4 in jet combination.' STOP ENDIF C C At this point we have three protojets and we find the YIJ to C combine them into two. C Y2JET = 1.0D0 OK = .FALSE. DO I = 2,3 DO J = 1,I-1 DOTIJ = 0.0D0 DO MU = 1,3 DOTIJ = DOTIJ + KJET(I,MU)*KJET(J,MU) ENDDO COSIJ = DOTIJ/ABSKJET(I)/ABSKJET(J) YIJ = 1.0D0 - COSIJ YIJ = YIJ * 2.0D0 * MIN(KJET(I,0),KJET(J,0))**2 YIJ = YIJ / RTS**2 IF (YIJ.LE.Y2JET) THEN Y2JET = YIJ OK = .TRUE. ENDIF ENDDO ENDDO IF (.NOT.OK) THEN WRITE(NOUT,*) 'Problem in second jet combination.' STOP ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION CALS0(NCUT,KCUT) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER NCUT REAL*8 KCUT(SIZE+1,0:3) C C Average value of (1 - thrust)**2. C C 10 March 1998 C 20 August 1999 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 C Use the code that is commented out below to add ten times the average C value of y_cut^2 for the differential three jet cross section with C the Durham algorithm to the computed result. C C REAL*8 Y4JET,Y2JET C C CALL COMBINEJETS(NCUT,KCUT,Y4JET,Y2JET) C IF (Y2JET.GT.Y4JET) THEN C CALS0 = CALS0 + 10.0D0 * (Y2JET**2 - Y4JET**2) C ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE DAYTIME(DATE) CHARACTER*28 DATE C C Returns DATE = today's date and time as a character string. C This version for Sun computers. Replace if necessary, for C instance with C date = 'today' C CHARACTER*30 DT CALL FDATE(DT) DATE = DT(1:28) RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE TIMING(DELTATIME) REAL*8 DELTATIME C C You call TIMING multiple times. Each time that you call it, C DELTATIME gets set to the elapsed time since the previous time you C called it. This version for Sun computers. Replace if necessary. C REAL*4 DTIME,TIMEARRAY(2) DELTATIME = DTIME(TIMEARRAY) RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE GETCOUNTS(COUNTFACTOR) C C Out: INTEGER MAXGRAPHS,MAXMAPS PARAMETER (MAXGRAPHS = 10) PARAMETER (MAXMAPS = 64) INTEGER COUNTFACTOR(MAXGRAPHS,MAXMAPS) C C This subroutine contains information on how many points to devote C to each of the ten graphs and, within each graph, how many to C devote to each of the maps from the random numbers x to the loop C momenta. C C Latest revision 29 September 1999. C INTEGER GRAPHNUMBER,MAPNUMBER C DO GRAPHNUMBER = 1,MAXGRAPHS DO MAPNUMBER = 1,MAXMAPS COUNTFACTOR(GRAPHNUMBER,MAPNUMBER) = 10 ENDDO ENDDO C C-------- COUNTFACTOR( 1, 1) = 57 COUNTFACTOR( 1, 2) = 58 COUNTFACTOR( 2, 1) = 29 COUNTFACTOR( 2, 2) = 26 COUNTFACTOR( 3, 1) = 86 COUNTFACTOR( 3, 2) = 89 COUNTFACTOR( 4, 1) = 13 COUNTFACTOR( 4, 2) = 13 COUNTFACTOR( 4, 3) = 13 COUNTFACTOR( 4, 4) = 14 COUNTFACTOR( 5, 1) = 51 COUNTFACTOR( 5, 2) = 34 COUNTFACTOR( 5, 3) = 49 COUNTFACTOR( 5, 4) = 17 COUNTFACTOR( 5, 5) = 4 COUNTFACTOR( 6, 1) = 52 COUNTFACTOR( 6, 2) = 18 COUNTFACTOR( 6, 3) = 3 COUNTFACTOR( 6, 4) = 37 COUNTFACTOR( 6, 5) = 44 COUNTFACTOR( 7, 1) = 247 COUNTFACTOR( 7, 2) = 195 COUNTFACTOR( 7, 3) = 33 COUNTFACTOR( 7, 4) = 23 COUNTFACTOR( 7, 5) = 255 COUNTFACTOR( 7, 6) = 192 COUNTFACTOR( 7, 7) = 36 COUNTFACTOR( 7, 8) = 28 COUNTFACTOR( 8, 1) = 55 COUNTFACTOR( 8, 2) = 2 COUNTFACTOR( 8, 3) = 68 COUNTFACTOR( 8, 4) = 19 COUNTFACTOR( 8, 5) = 43 COUNTFACTOR( 8, 6) = 35 COUNTFACTOR( 8, 7) = 52 COUNTFACTOR( 8, 8) = 1 COUNTFACTOR( 9, 1) = 3 COUNTFACTOR( 9, 2) = 3 COUNTFACTOR( 9, 3) = 5 COUNTFACTOR( 9, 4) = 4 COUNTFACTOR(10, 1) = 5 COUNTFACTOR(10, 2) = 1 COUNTFACTOR(10, 3) = 14 COUNTFACTOR(10, 4) = 1 COUNTFACTOR(10, 5) = 4 COUNTFACTOR(10, 6) = 13 COUNTFACTOR(10, 7) = 10 COUNTFACTOR(10, 8) = 1 COUNTFACTOR(10, 9) = 15 COUNTFACTOR(10,10) = 1 COUNTFACTOR(10,11) = 8 COUNTFACTOR(10,12) = 11 C-------- C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012