!2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! Module for program beowulf module beowulf_helpers implicit none integer, parameter :: dbl = kind(1.0d0) integer, parameter :: long = selected_int_kind(15) integer, parameter :: nin = 5 integer, parameter :: nout = 6 integer, parameter :: size = 3 integer, parameter :: maxgraphs = 12 integer, parameter :: maxcuts = 9 integer, parameter :: maxmaps = 64 integer, parameter :: maxparticles = 500 end module beowulf_helpers ! Note that for many compilers, the modules have to come before the ! programs in which they are used. !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! ! ------------------------ ! beowulf.f Version 3.1 ! ------------------------ ! ! Program for numerical integration of jet cross sections in ! electron-positron annihilation. -- D. E. Soper ! ! First code: 29 November 1992. ! Latest revision: see revision date below. ! ! The main program and subroutines that a user might want to modify ! are contained in this package, beowulf.f. In particular, a ! user may very well want to modify parameter settings in the main ! program and to change the observables calculated in the subroutine ! HROTHGAR and in the functions CALStype(nparticles,kf,index). Subroutines ! that can be modified only at extreme peril to the reliability of ! the results are in the companion package beowulfsubs.f. ! program beowulf ! use beowulf_helpers implicit none character(len=36), parameter :: & revisiondate = 'Latest revision 2 October 2003 ' ! ------------------------------------ ! Input and output units: 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 ! Global size and scale parameters: integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 ! Color factors: real(kind=dbl) :: nc,nf,cf common /colorfactors/ nc,nf,cf ! Choice of gauge: character(len=7) :: gauge common /gaugechoice/ gauge ! What the program should do character(len=6) :: mode common /programmode/ mode ! Physics data: real(kind=dbl) :: alphasofmz,mz,externalrts common /physicsdata/ alphasofmz,mz,externalrts ! How many graphs and how many maps for each: integer :: numberofgraphs integer :: numberofmaps(maxgraphs) common /graphcounts/ numberofgraphs,numberofmaps ! Limits: real(kind=dbl) :: badnesslimit,cancellimit,thrustcut common /limits/ badnesslimit,cancellimit,thrustcut real(kind=dbl) :: showercut common /showercutinfo/ showercut real(kind=dbl) :: showerendratio common /showerend/ showerendratio logical :: dosoft common /softswitch/ dosoft ! Timing limit: real(kind=dbl) :: hourlimit real(kind=dbl) :: timelimit common /maxtime/ timelimit ! Deform parameters: real(kind=dbl) :: deformalpha,deformbeta,deformgamma common /deformscales/deformalpha,deformbeta,deformgamma ! Smear parameters: real(kind=dbl) :: smearfctr integer :: lowpwr,highpwr common /smearparms/ smearfctr,lowpwr,highpwr ! Renormaliation parameters: real(kind=dbl) :: muoverrts common /renormalize/ muoverrts ! Diagnostics on showers: logical :: showercheck common /showercheckinfo/ showercheck ! Diagnostic data instructions: logical :: report,details common /calculook/ report,details ! RENO results: real(kind=dbl) :: sumr,errorr,sumi,errori real(kind=dbl) :: sumbis,errorbis real(kind=dbl) :: sumchkr,errorchkr,sumchki,errorchki real(kind=dbl) :: fluct(maxgraphs,maxmaps) integer(kind=long) :: included,extracount,omitted integer :: nvalpt(-9:6) real(kind=dbl) :: valptmax type badpointstate integer :: jr integer :: ir(250) real(kind=dbl) :: rr(250) integer :: graphnumber integer :: mapnumber real(kind=dbl) :: k(0:3*size-1,0:3) end type badpointstate type(badpointstate) :: badpointinfo integer :: nreno real(kind=dbl) :: cputime ! Hrothgar dummy variables: real(kind=dbl) :: kf0(maxparticles,0:3) ! Seed for random number generator. integer :: seed ! Today's date: character(len=28) :: date ! Loop control variables: integer :: n,mu,graphnumber,mapnumber ! Logical variables that control aspects of the output: logical :: allgraphs,writeprograminfo,writefluct,writediagnostic ! A temporary real number: real(kind=dbl) :: temp ! A temporary integer: integer :: itemp ! Average of FLUCT variables: real(kind=dbl) :: avfluct ! alpha_s/pi real(kind=dbl)alpi real (kind=dbl), parameter :: pi = 3.141592653589793239d0 ! !----------------------------------------------------------------------- ! ! Parameters to control diagnostic reporting. report = .false. details = .false. ! Generate diagnostics on showers. showercheck = .false. ! Global size parameters: Number of loops, number of propagators ! number of vertices, max number of cut propagators. nloops1 = 2 nprops1 = 3 * nloops1 - 1 nverts1 = 2 * nloops1 cutmax1 = nloops1 + 1 nloops2 = 3 nprops2 = 3 * nloops2 - 1 nverts2 = 2 * nloops2 cutmax2 = nloops2 + 1 ! Color factors. nc = 3.0d0 nf = 5.0d0 cf = (nc**2 - 1)/2.0d0/nc ! Physics data. mz = 91.1876d0 alphasofmz = 0.118d0 externalrts = mz ! Initialize common /GRAPHCOUNTS/. call graphcountinit ! Parameters for subroutine CALCULATE for throwing away bad points. badnesslimit = 1.0d4 cancellimit = 1.0d4 ! Gradual cutoff on 1 - THRUST before showering in case of mode 'shower'. ! IF (onemthrust .LT. thrustcut) THEN ! value = onemthrust/thrustcut*value ! weight = onemthrust/thrustcut*weight ! END IF thrustcut = 0.00d0 ! Cut, in shower mode, for virtuality treated as parton splitting from ! a Born graph: qbarsq < showercut * calqsq. showercut = 1.0d0 ! Minimum value of qbarsq/rts0^2 in a secondary parton shower (ie showerII). ! If qbarsq/rts0^2 comes out less than this, the splitting is vetoed and the ! mother parton is not split. showerendratio = 3.0d-4 ! Switch for including soft gluon effects. Setting this to .false. gives ! correct NLO results, but with showering from the order alpha_s^3 graphs ! turned off and soft gluon radiation from the Born graphs turned off. dosoft = .true. ! Parameters for subroutine DEFORM. deformalpha = 0.7d0 deformbeta = 1.0d0 deformgamma = 1.0d0 ! Parmameters for function SMEAR. smearfctr = 3.0d0 lowpwr = 9 highpwr = 18 ! Renormalization parameters. muoverrts = 0.16667d0 ! 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. usegraph(11) = .true. usegraph(12) = .true. ! write(nout,*) & 'Please give program mode (born, nlo, hocoef, shower, showr0).' read(nin,98)mode 98 format(a6) ! ! Depending on the mode, we have to set USEGRAPH to .false. for those ! graphs not used in that mode. ! IF (mode.EQ.'born ') THEN usegraph(1) = .false. usegraph(2) = .false. usegraph(3) = .false. usegraph(4) = .false. usegraph(5) = .false. usegraph(6) = .false. usegraph(7) = .false. usegraph(8) = .false. usegraph(9) = .false. usegraph(10) = .false. ELSE IF (mode.EQ.'hocoef') THEN usegraph(11) = .false. usegraph(12) = .false. ELSE IF (mode.EQ.'nlo ') THEN continue ELSE IF (mode.EQ.'shower') THEN continue ELSE IF (mode.EQ.'showr0') THEN usegraph(1) = .false. usegraph(2) = .false. usegraph(3) = .false. usegraph(4) = .false. usegraph(5) = .false. usegraph(6) = .false. usegraph(7) = .false. usegraph(8) = .false. usegraph(9) = .false. usegraph(10) = .false. mode = 'shower' ELSE write(nout,*)'Not prepared for mode ',mode,'.' stop END IF ! write(nout,*)'Please give alpha_s/0.118' read(nin,*)temp alphasofmz = 0.118d0*temp ! IF (mode.EQ.'shower') THEN gauge = 'coulomb' ELSE write(nout,*) & 'Please give gauge choice (coulomb or feynman).' read(nin,99)gauge 99 format(a7) END IF ! ! Here we set the count factors COUNTFACTOR(GRAPHNUMBER,MAPNUMBER). ! These, multiplied by NSETS, give GROUPSIZE(GRAPHNUMBER,MAPNUMBER). ! IF (mode.EQ.'shower') THEN call getcountsshwr(countfactor) ELSE IF (gauge.EQ.'feynman') THEN call getcountsfeyn(countfactor) ELSE IF (gauge.EQ.'coulomb') THEN call getcountscoul(countfactor) ELSE write(*,*)'Not programmed for gauge choice ', gauge STOP END IF ! ! Call Hrothgar to tell him to get ready. ! Here kf0 are just dummy variables for Hrothgar. ! DO n = 1,maxparticles DO mu = 0,3 kf0(n,mu) = 0.0d0 END DO END DO call hrothgar(1,kf0,1.0d0,1,'startcalc ') !--- 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) END IF !--- ! To change the default value calculated above, use this. ! ! WRITE(NOUT,*)'Please give number of sets of points in each group' ! READ(NIN,*) NSETS !--- 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 END DO END IF END DO !-- ! If you want to extract parts of the results with different color ! and number-of-flavor factors, you may want to be able to ! adjust the number of colors and flavors. But note that the ! function KN(T) that contains the "standard" results ! for the thrust distribution is for three colors and five flavors. ! ! WRITE(NOUT,*)'Please give the numbers of colors and flavors.' ! READ(NIN,*)NC,NF !-- ! This overrides the default value MUOVERRTS = 0.3 that was set above. ! write(nout,*) & 'Please give ratio of the renormalization scale to sqrt(s).' read(nin,*) muoverrts ! write(nout,*)'Please give a seed (0<i<259200) for random numbers.' read(nin,*) seed call randominit(seed) ! call daytime(date) write(nout,100)date 100 format(// & ,' beowulf version 3.1 ',a,/,50('-')) ! !**************************************************** write(nout,*) revisiondate !**************************************************** ! call version ! 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 I''ll 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,'.') IF (gauge.EQ.'coulomb') THEN write(nout,103) 103 format('beowulf will use the Coulomb gauge and calculate') ELSE IF (gauge.EQ.'feynman') THEN write(nout,104) 104 format('beowulf will use the Feynman gauge and calculate') ELSE write(nout,*)'Not familiar with gauge choice ',gauge,'.' stop END IF ! IF (mode.EQ.'born ') THEN write(nout,105) 105 format('the Born contribution to observable/sigma_0.') write(nout,106)alphasofmz,mz,externalrts 106 format('alpha_s(mz) is',f7.4,' with mz = ',f7.4,'.',/, & 'sqrt(s) is ',1p g11.6,'.') ELSE IF (mode.EQ.'nlo ') THEN write(nout,107) 107 format('observable/sigma_0 ', & 'at next-to-leading order.') write(nout,106)alphasofmz,mz,externalrts ELSE IF (mode.EQ.'shower') THEN write(nout,108) 108 format('observable/sigma_0 ', & 'at next-to-leading order with showers.') write(nout,106)alphasofmz,mz,externalrts ELSE IF (mode.EQ.'hocoef') THEN write(nout,109) 109 format('the coefficient of (alpha_s/pi)^2 ', & 'for observable/sigma_0.') ELSE IF (mode.EQ.'shower') THEN write(nout,110) 110 format('the coefficient of (alpha_s/pi)^2 ', & 'for observable/sigma_0 with showers.') ELSE write(nout,*)'Not prepared for mode ',mode,'.' stop END IF ! write(nout,111)nc,nf 111 format('beowulf will use',1p g9.2,'colors and',1p g9.2,'flavors.') write(nout,112) 112 format('Renormalization parameter:') write(nout,113)muoverrts 113 format(' mu /sqrt(s):',1p g11.4) write(nout,1827)pi*alpi(muoverrts*externalrts) 1827 format(' alpha_s(mu):',1p g11.4) write(nout,114) 114 format('Cutoff parameters:') write(nout,115)badnesslimit,cancellimit 115 format(' badness limit:',1p g9.2,/, & ' cancellation limit:',1p g9.2,/) IF (mode.EQ.'shower') THEN IF (abs(showercut - 1.0d0).GT.0.00001) THEN write(nout,751)showercut 751 format('Born showering includes virtualities 0 < qbar^2 < ', & f4.2,'*\vec q^2.',/) END IF write(nout,752) showerendratio 752 format('Secondary showering stops at virtuality qbar^2/rts0^2 =', & 1p g9.2) IF (.NOT.dosoft) THEN write(nout,753) 753 format('Showering from higher order graphs and soft gluons OFF.') END IF IF (thrustcut.GT.0.0d0) THEN write(nout,754)thrustcut 754 format('Weights suppressed if 1 - thrust before showering is less than', & f8.4) END IF END IF ! allgraphs = .true. IF (mode.EQ.'born ') THEN DO n = 11,maxgraphs IF (.NOT.usegraph(n) ) THEN allgraphs = .false. END IF END DO ELSE IF (mode.EQ.'hocoef') THEN DO n = 1,10 IF (.NOT.usegraph(n) ) THEN allgraphs = .false. END IF END DO ELSE IF (mode.EQ.'nlo ') THEN DO n = 1,maxgraphs IF (.NOT.usegraph(n) ) THEN allgraphs = .false. END IF END DO ELSE IF (mode.EQ.'shower') THEN DO n = 1,maxgraphs IF (.NOT.usegraph(n) ) THEN allgraphs = .false. END IF END DO ELSE write(nout,*)'Not prepared for mode ',mode,'.' stop END IF IF (allgraphs) THEN write(nout,116) 116 format('All relevant graphs are used.') ELSE DO n = 1,maxgraphs write(nout,117) n,usegraph(n) 117 format('Use of graph ',i2,': ',l1) END DO END IF ! 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) END DO write(nout,125)groupsizegraph(graphnumber) 125 format(/,'Total points for this graph',i8,/) END IF END DO 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) END IF ! write(nout,*)' ' ! call reno( & sumr,errorr,sumi,errori, & sumbis,errorbis, & sumchkr,errorchkr,sumchki,errorchki,fluct, & included,extracount,omitted, & nvalpt,valptmax,badpointinfo, & nreno,cputime) ! 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.',/) ! 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,*)' ' ! 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) ! ! Check integral. ! IF (mode.NE.'shower') THEN write(nout,604) 604 format(/,'Check integral, equal to 1.0 plus a cutoff error.',/, & 'for badness limit = cancellation limit = 1.00e+04,',/, & 'and mode = hocoef, it 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) ! write(nout,*)' ' END IF ! ! Call Hrothgar to tell him that we are done with the calculation. ! call hrothgar(1,kf0,1.0d0,nreno,'calcdone ') ! write(nout,*)' ' write(nout,*)' ' write(nout,700) 700 format('***********************') write(nout,701) 701 format('diagnostic information: ') write(nout,*)' ' ! 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) END DO ! write(nout,*)' ' write(nout,703)valptmax 703 format('Biggest contribution was',1p g12.3) write(nout,704)badpointinfo%graphnumber,badpointinfo%mapnumber 704 format('from graph',i3,', point choice method',i3) ! ! Setting WRITEFLUCT to .TRUE. will result in the generation of ! information where the fluctions in SUMR come from, including ! a suggested set of COUNTFACTOR values that should make the ! fluctions in SUMR smaller. ! 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 END DO END DO 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) END DO END DO 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 END IF write(nout,706) graphnumber,mapnumber,itemp 706 format(' countfactor(',i2,',',i2,') = ',i4) END DO END DO write(nout,*)' ' END IF ! writediagnostic = .true. IF (writediagnostic) THEN call diagnostic(badpointinfo) END IF ! !---------------------------------------------- ! call daytime(date) write(nout,707)date 707 format(/ & ,'Done ',a,/) ! !--- ! stop END program beowulf ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine hrothgar(nparticles,kf,weight,nreno,whattodo) ! use beowulf_helpers implicit none ! In: integer :: nparticles real(kind=dbl) :: kf(maxparticles,0:3) real(kind=dbl) :: weight integer :: nreno character(len=10) :: whattodo ! ! Accumulates the results. ! ! The results are stored as they accumulate in arrays of size NBINSTOT. ! Each array element corresponds to a different observable -- for ! example a certain moment of the thrust distribution. ! The arrays are divided into subsets of size NBINS1, NBINS2, ... , ! with NBINS1 + NBINS2 + ... = NBINSTOT. Each subset corresponds ! to a different type of observable, such as "moments of the ! thrust distribution." The calculations of the ! observation function for each type of observable is done ! by a function named CALStype(nparticles,kf,index), where the index ! for type j takes values 1,...,NBINSj. Users can remove ! observable types and can add new observable types by writing new ! functions CALStype(nparticles,kf,index) and by modfifying the reporting ! of the results, which is controlled by this subroutine. ! ! Note that, in addition to sending sets of final state momenta along ! with weights to HROTHGAR, the other subroutines of beowulf compute ! a single integral, which is defined by the observable in the ! function CALS0(nparticles,kf) and is the average value of (1-thrust)^2 ! in the default version of the beowulf.f subroutines. If there is ! too much cancellation in the calculation of this observable among ! the contributions from different final states corresponding to ! a given point in the space of loop momenta, then HROTHGAR is sent ! a BADPOINT or XTRAPOINT signal. ! ! The action of HROTHGAR depends on value of WHATTODO. ! ! 'STARTCALC ': Initialize variables at start of run. ! 'STARTGROUP': Initialize variables at start of new group of points. ! 'STARTPOINT': Initialize variables at start of a new point in the ! space of loop momenta. ! 'BADPOINT ': The point is too near a singularity ! OR there is too much cancellation. ! Set switch BADPOINTQ to true. ! 'XTRAPOINT ': The point is too near a singularity ! OR there is too much cancellation. ! But not so bad as to trigger a BADPOINT signal. ! Set switch XTRAPOINTQ to true. ! 'NEWRESULT ': We have final state momentum variables nparticles, kf ! and corresponding WEIGHT for a given final state cut. ! We compute the corresponding INTEGRAND = WEIGHT*CALSVAL ! and add it to VALUE for the current point. In Hrothgar, ! WEIGHT and VALUE are real. ! 'POINTDONE ': We now have the accumulated VALUE for this point. ! *If we we earlier set BADPOINTQ to true, then we ! simply throw away this point. If BADPOINTQ is false ! we proceed according to whether there was too much ! cancellation or not: ! *If XTRAPOINTQ is false we are ok and ! we add VALUE to INTEGAL ! *If XTRAPOINTQ is true ! we add VALUE to INTEGALBIS ! 'GROUPDONE ': We increment the sums for this group. This includes ! sums for the observables being computed and also ! sums to use for computing the statistical errors. ! 'CALCDONE ': Compute the results and print them. ! ! The limits: ! Cancel<Lim Lim<Cancel<100*Lim 100*Lim<Cancel ! Bad < Lim OK XTRAPOINT BADPOINT ! Lim<Bad<100*Lim XTRAPOINT XTRAPOINT BADPOINT ! 100*Lim < Bad BADPOINT BADPOINT BADPOINT ! ! The BADPOINT signal overrides the XTRAPOINT signal in Hrothgar. ! ! 26 February 1998 ! Latest version: ! 3 July 2002 ! real(kind=dbl) :: badnesslimit,cancellimit,thrustcut common /limits/ badnesslimit,cancellimit,thrustcut ! save sum,sumbis,sum2,sum2bis,sum3,sum4,integral,integralbis, & value,maxpart,badpointq,xtrapointq ! integer, parameter :: nbins1 = 9 integer, parameter :: nbins2 = 10 integer, parameter :: nbinstot = 19 real(kind=dbl) :: sum(nbinstot),sumbis(nbinstot) real(kind=dbl) :: sum2(nbinstot),sum2bis(nbinstot) real(kind=dbl) :: sum3(nbinstot),sum4(nbinstot) real(kind=dbl) :: integral(nbinstot),integralbis(nbinstot) real(kind=dbl) :: value(nbinstot),maxpart(nbinstot) real(kind=dbl) :: error(nbinstot),error4(nbinstot),errorbis(nbinstot) ! logical :: badpointq,xtrapointq real(kind=dbl) :: calsvalue,calsthrust,calstmoments,calsymoments integer :: n,ntemp real(kind=dbl) :: thrustdist,thrustval,thrust ! 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 END DO ELSE IF (whattodo .EQ. 'startgroup') THEN DO n = 1,nbinstot integral(n) = 0.0d0 integralbis(n) = 0.0d0 END DO ELSE IF (whattodo .EQ. 'startpoint') THEN DO n = 1,nbinstot value(n) = 0.0d0 END DO badpointq = .false. xtrapointq = .false. ELSE IF (whattodo .EQ. 'badpoint ') THEN badpointq = .true. ELSE IF (whattodo .EQ. 'xtrapoint ') THEN xtrapointq = .true. !--------------------------------------------------------------- ! Here is where we calculate the functions CALS for each type of ! observable using functions CALStype(nparticles,kf,index). ! ELSE IF (whattodo .EQ. 'newresult ') THEN thrustval = thrust(kf,nparticles) n = 0 DO ntemp = 1,nbins1 n = n+1 calsvalue = calsthrust(thrustval,ntemp) value(n) = value(n) + weight*calsvalue END DO DO ntemp = 1,nbins2 n = n+1 calsvalue = calstmoments(thrustval,ntemp) value(n) = value(n) + weight*calsvalue END DO !--------------------------------------------------------------- ELSE IF (whattodo .EQ. 'pointdone ') THEN IF (.NOT.badpointq) THEN IF (xtrapointq) THEN DO n = 1,nbinstot integralbis(n) = integralbis(n) + value(n) END DO ELSE DO n = 1,nbinstot integral(n) = integral(n) + value(n) END DO END IF END IF ELSE IF (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 END DO ELSE IF (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 END DO !------------------------------------------------------- ! Write the labels for each subset of our results arrays ! and write the results themselves. ! n = 0 DO ntemp = 1,nbins1 n = n+1 write(nout,002)(0.68d0 + ntemp*0.03d0), & thrustdist(0.68d0 + ntemp*0.03d0) write(nout,003)sum(n),error(n) write(nout,004)sumbis(n),errorbis(n) write(nout,005)error4(n) END DO 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) END DO !------------------------------------------------------- END IF RETURN ! ! FORMAT statements for the output. ! 001 format(/, ' --- Hrothgar reports ---',/,/, & 'First, the average near the thrust t of the thrust ',/, & 'distribution ((1-t)/sigma_0) d sigma /d t divided by',/, & 'an approximation to the same quantity.'/) 002 format('t =',f6.3,' approximation function =',1p g10.3) 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) END subroutine hrothgar ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function calsthrust(t,n) ! use beowulf_helpers implicit none ! In: integer :: n ! bin number real(kind=dbl) :: t ! value of thrust ! Out: real(kind=dbl) :: calsthrust ! ! The thrust distribution in bins. ! The function THRUSTDIST(T) gives an approximation to the perturbative ! expansion of the thrust distribution (1/sigma_0) * (1-T) d sigma/dT. ! The result depends on Mode. ! ! 9 March 1998 ! 5 August 1999 ! 3 July 2002 ! real(kind=dbl) :: t1,t2 real(kind=dbl) :: rho,thrustdist real(kind=dbl), parameter :: deltathrust = 0.03d0 real(kind=dbl), parameter :: minthrust = 0.68d0 ! calsthrust = 0.0d0 ! t1 = minthrust + (n-1) * deltathrust t2 = minthrust + (n+1) * deltathrust 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/thrustdist(t) END IF ! RETURN END function calsthrust ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function calstmoments(t,n) ! use beowulf_helpers implicit none ! In: real(kind=dbl) :: t ! value of thrust integer :: n ! moment number ! Out: real(kind=dbl) :: calstmoments ! ! Moments of 1 - thrust. ! 20 September 1999 ! 3 July 2002 ! real(kind=dbl) :: power integer :: i real(kind=dbl), parameter :: tiny = 1.0d-12 ! calstmoments = 0.0d0 ! power = 1.0d0 + 0.5d0*n IF ( (1.0d0 - t) .GT. tiny ) THEN calstmoments = (1.0d0 - t)**power ELSE calstmoments = 0.0d0 END IF ! RETURN END function calstmoments ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function thrust(kf,nparticles) ! kf gives momenta for nparticles final state particles. ! Then thrust is the thrust of this final state. use beowulf_helpers implicit none real(kind=dbl) :: thrust real(kind=dbl) :: kf(maxparticles,0:3) integer :: nparticles ! ! One can increase imax for a more accurate answer at the expense of ! more time. ! For imax = 1, there are 3 starting thrust axes used. ! For imax = 2, there are 9 starting thrust axes used. ! For imax = 3, there are 20 starting thrust axes used. ! Approximately, there are somewhat fewer than 3*(imax)^2 starting ! thrust axes used. ! Here are some results for the largest error after 100 events. ! These are without the auxiliary trial axes choices. ! imax = 1, <thrust> ~ 0.80, maxerror = 0.0002 ! imax = 1, <thrust> ~ 0.65, maxerror = 0.004 ! imax = 1, <thrust> ~ 0.55, maxerror = 0.02 ! imax = 3, <thrust> ~ 0.55, maxerror = 0.003 ! imax = 10, <thrust> ~ 0.55, maxerror = 0.0002 ! later investigations with perturbative final states suggest that ! imax = 2 (with the auxiliary choices) gets the thrust right to 1E-3. ! ! 3 July 2002 ! real(kind=dbl) :: rts integer, parameter :: imax = 2 integer :: jmax integer :: i,j,n,m,mu real(kind=dbl), parameter :: tiny = 1.0d-3 real (kind=dbl), parameter :: pi = 3.141592653589793239d0 real(kind=dbl) :: sintheta,costheta,theta,phi real(kind=dbl), dimension(3) :: vin,vout real(kind=dbl) :: temp real(kind=dbl) :: thrustmod,thrusttest ! thrust = 0.0d0 ! Calculate rts. ! rts = 0.0d0 DO n = 1,nparticles rts = rts + kf(n,0) END DO ! ! We will loop over choices for the trial thrust axis. ! 1 - 0.391827 is (pi/2)*arccos(1/sqrt(3)), which is the angle ! from the center of one face of a regular solid with six vertices ! to one of the adjacent vertices. This optimizes the choice for ! imax = 1. ! DO i = 1,imax theta = (i - 0.391827d0)*0.5d0*pi/imax sintheta = sin(theta) costheta = cos(theta) jmax = nint(3.3d0*sintheta*imax) DO j = 1,jmax phi = (j - 0.5d0)*2.0*pi/jmax ! vin(1) = sintheta*cos(phi) vin(2) = sintheta*sin(phi) vin(3) = costheta ! ! Adjust thrust axis ! DO n = 1,100 call thrustaxis(kf,nparticles,vin,vout,thrustmod) temp = 0.0d0 DO mu = 1,3 temp = temp + (vin(mu) - vout(mu))**2 END DO IF (temp < tiny) EXIT vin = vout END DO IF (n > 30) print *,'In function thrust, n tries was ',n ! ! Now we have the thrust axis. Find trial value of thrust ! thrusttest = 2.0d0*thrustmod/rts ! ! If this is the biggest value so far, we save it ! IF (thrusttest > thrust) thrust = thrusttest ! END DO END DO ! !-- Auxilliary choices in the case of small thrust -- ! IF (thrust < 0.67d0) THEN ! ! We will loop over more choices for the trial thrust axis. ! DO i = 1,3*imax theta = (i - 0.391827d0)*0.5d0*pi/3/imax sintheta = sin(theta) costheta = cos(theta) jmax = nint(3.3d0*sintheta*3*imax) DO j = 1,jmax phi = (j - 0.5d0)*2.0*pi/jmax ! vin(1) = sintheta*cos(phi) vin(2) = sintheta*sin(phi) vin(3) = costheta ! ! Adjust thrust axis ! DO n = 1,100 call thrustaxis(kf,nparticles,vin,vout,thrustmod) temp = 0.0d0 DO mu = 1,3 temp = temp + (vin(mu) - vout(mu))**2 END DO IF (temp < tiny) EXIT vin = vout END DO IF (n > 30) print *,'In function thrust, n was ',n ! ! Now we have the thrust axis. Find trial value of thrust ! thrusttest = 2.0d0*thrustmod/rts ! ! If this is the biggest value so far, we save it ! IF (thrusttest > thrust) thrust = thrusttest ! END DO END DO ! END IF ! END function thrust ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine thrustaxis(kf,nparticles,vin,vout,thrustmod) ! use beowulf_helpers implicit none real(kind=dbl) :: kf(maxparticles,0:3) integer :: nparticles real(kind=dbl), dimension(3) :: vin,vout real(kind=dbl) :: thrustmod ! thrust*rts/2 ! ! kf gives momenta for nparticles final state particles. ! vin is the trial thrust axis. ! Then vout is the new thrust axis. ! Also returns thrustmod, the value of thrust*rts/2 if this ! is the final thrust axis. ! ! 3 July 2002 ! real(kind=dbl) :: temp integer n,mu ! thrustmod = 0.0d0 DO mu = 1,3 vout(mu) = 0.0d0 END DO DO n = 1, nparticles temp = 0.0d0 DO mu = 1,3 temp = temp + kf(n,mu) * vin(mu) END DO IF (temp>0.0d0) THEN DO mu = 1,3 vout(mu) = vout(mu) + kf(n,mu) END DO thrustmod = thrustmod + temp END IF END DO temp = 0.0d0 DO mu = 1,3 temp = temp + vout(mu)**2 END DO temp = sqrt(temp) DO mu = 1,3 vout(mu) = vout(mu)/temp END DO ! END subroutine thrustaxis ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function thrustdist(t) ! use beowulf_helpers implicit none ! In: real(kind=dbl) :: t ! Out: real(kind=dbl) :: thrustdist ! ! Approximation to the perturbative expansion of the thrust distribution ! (1/sigma_0) * (1-T) d sigma/dT. The result depends on Mode. Defining ! f(T) = (1/sigma_0) * (1-T) d sigma/dT ! = (alpha_s/Pi) KN0(T) + (alpha_s/Pi)^2 KN(T) + ... ! we calculate for mode = born, ! THRUSTDIST(T) = (alpha_s/Pi) KN0(T), ! for Mode = nlo, ! THRUSTDIST(T) = (alpha_s/Pi) KN0(T) + (alpha_s/Pi)^2 KN(T), ! for Mode = hocoef ! THRUSTDIST(T) = KN(T). ! ! 12 February 2002 ! ! What the program should do character(len=6) :: mode common /programmode/ mode ! Physics data real(kind=dbl) :: alphasofmz,mz,externalrts common /physicsdata/ alphasofmz,mz,externalrts real(kind=dbl) :: muoverrts common /renormalize/ muoverrts ! real(kind=dbl) :: alpi,kn0,kn ! IF (mode.EQ.'born ') THEN thrustdist = alpi(muoverrts*externalrts)*kn0(t) ELSE IF (mode.EQ.'nlo ') THEN thrustdist = alpi(muoverrts*externalrts)*kn0(t) & + alpi(muoverrts*externalrts)**2 * kn(t) ELSE IF (mode.EQ.'hocoef') THEN thrustdist = kn(t) ELSE IF (mode.EQ.'shower') THEN thrustdist = alpi(muoverrts*externalrts)*kn0(t) & + alpi(muoverrts*externalrts)**2 * kn(t) ELSE write(nout,*)'thrustdist misses this mode. ',mode stop END IF RETURN END function thrustdist ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function kn0(t) ! use beowulf_helpers implicit none ! In: real(kind=dbl) :: t ! Out: real(kind=dbl) :: kn0 ! ! The coefficient of (alpha_s/Pi)^1 in the thrust distribution, ! (1/sigma_0) * (1-T) d sigma/dT. ! The formula is copied from Ellis, Stirling and Webber eq. (3.44). ! Note that this function is only for NC = 3, NF = 5. ! ! 11 February 2002 ! kn0 = 2.0d0 * (3.0d0 * t**2 - 3.0d0 * t + 2.0d0) /t kn0 = kn0 * log( (2.0d0 * t - 1.0d0)/(1.0d0 - t) ) kn0 = kn0 - 3.0d0 * (3.0d0 * t - 2.0d0) * (2.0d0 - t) kn0 = kn0 * 2.0d0/3.0d0 ! RETURN END function kn0 ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function kn(t) ! use beowulf_helpers implicit none ! In: real(kind=dbl) :: t ! Out: real(kind=dbl) :: kn ! ! The coefficient of (alpha_s/Pi)^2 in the thrust distribution as ! given by Kunszt and Nason. ! Note that this function is only for NC = 3, NF = 5, ! and is for (1/sigma_0) * (1-T) d sigma/dT. ! ! 9 April 1998 ! real(kind=dbl) :: muoverrts common /renormalize/ muoverrts real(kind=dbl), parameter :: nf = 5.0d0 ! real(kind=dbl) :: l,nonsing,sing,lo l = log(1.0d0/(1.0d0 - t)) ! nonsing = -16100.550195d0 * (0.46789210556 + t) & * (1.1184765748d0 - 2.1058049743d0 * t + t**2) & * (0.4925129711d0 - 1.3620871313d0 * t + t**2) ! sing = -14.222222222d0 * l * (-5.098072232d0 + l ) & *( 0.691822232d0 + l) ! 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 ! kn = sing + nonsing & + (11.0d0 - 2.0d0*nf/3.0d0)*log(muoverrts)*lo ! ! We divide by 4 because KN report the coefficient of (alpha_s/2Pi)^2 ! while I want the coefficient of (alpha_s/Pi)^2. ! kn = kn/4.0d0 ! RETURN END function kn ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function cals0(nparticles,kf) ! use beowulf_helpers implicit none ! In: integer :: nparticles real(kind=dbl) :: kf(maxparticles,0:3) ! Out: real(kind=dbl) :: cals0 ! ! Average value of (1 - thrust)**2. ! ! 10 March 1998 ! 20 August 1999 ! real(kind=dbl) :: badnesslimit,cancellimit,thrustcut common /limits/ badnesslimit,cancellimit,thrustcut ! real(kind=dbl) :: rts real(kind=dbl) :: thrust,oneminusthrust integer :: i ! rts = 0.0d0 DO i = 1,nparticles rts = rts + kf(i,0) END DO ! oneminusthrust = 1.0d0 - thrust(kf,nparticles) ! cals0 = oneminusthrust**2 ! RETURN END function cals0 ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine daytime(date) ! use beowulf_helpers implicit none character(len=28) :: date ! ! Returns DATE = today's date and time as a character string. ! This version for Sun computers. Replace if necessary, for ! instance with ! date = 'today' ! character*10 thedate,thetime character(len=30) :: dt ! CALL FDATE(DT) ! DATE = DT(1:28) call date_and_time(thedate,thetime) date = thedate//' '//thetime RETURN END subroutine daytime ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine timing(deltatime) ! use beowulf_helpers implicit none real(kind=dbl) :: deltatime ! ! You call TIMING multiple times. Each time that you call it, ! DELTATIME gets set to the elapsed time since the previous time you ! called it. This version for Sun computers. Replace if necessary. ! ! REAL*4 DTIME,TIMEARRAY(2) ! DELTATIME = DTIME(TIMEARRAY) ! real :: t1,t2 data t1 /0.0/ call cpu_time(t2) deltatime = t2 - t1 t1 = t2 IF (deltatime .LT. 0.000d0) THEN write(nout,*)'deltatime =',deltatime,' replaced by 0.0d0.' deltatime = 0.0d0 END IF RETURN END subroutine timing ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine getcountscoul(countfactor) ! use beowulf_helpers implicit none ! Out: integer :: countfactor(maxgraphs,maxmaps) ! ! This subroutine contains information on how many points to devote ! to each of the ten graphs and, within each graph, how many to ! devote to each of the maps from the random numbers x to the loop ! momenta. These are for Coulomb gauge. ! ! Latest revision 15 February 2002. ! integer :: graphnumber,mapnumber ! DO graphnumber = 1,maxgraphs DO mapnumber = 1,maxmaps countfactor(graphnumber,mapnumber) = 20 END DO END DO ! !-------- countfactor( 1, 1) = 13 countfactor( 1, 2) = 39 countfactor( 2, 1) = 13 countfactor( 2, 2) = 18 countfactor( 3, 1) = 22 countfactor( 3, 2) = 12 countfactor( 3, 3) = 13 countfactor( 3, 4) = 15 countfactor( 3, 5) = 11 countfactor( 3, 6) = 23 countfactor( 4, 1) = 2 countfactor( 4, 2) = 2 countfactor( 4, 3) = 10 countfactor( 4, 4) = 9 countfactor( 5, 1) = 10 countfactor( 5, 2) = 3 countfactor( 5, 3) = 7 countfactor( 5, 4) = 4 countfactor( 5, 5) = 10 countfactor( 5, 6) = 5 countfactor( 6, 1) = 9 countfactor( 6, 2) = 4 countfactor( 6, 3) = 8 countfactor( 6, 4) = 3 countfactor( 6, 5) = 6 countfactor( 6, 6) = 8 countfactor( 7, 1) = 9 countfactor( 7, 2) = 5 countfactor( 7, 3) = 15 countfactor( 7, 4) = 3 countfactor( 7, 5) = 9 countfactor( 7, 6) = 1 countfactor( 7, 7) = 42 countfactor( 7, 8) = 12 countfactor( 7, 9) = 5 countfactor( 7,10) = 9 countfactor( 7,11) = 3 countfactor( 7,12) = 18 countfactor( 7,13) = 2 countfactor( 7,14) = 5 countfactor( 7,15) = 2 countfactor( 7,16) = 17 countfactor( 7,17) = 9 countfactor( 7,18) = 11 countfactor( 8, 1) = 1 countfactor( 8, 2) = 1 countfactor( 8, 3) = 3 countfactor( 8, 4) = 1 countfactor( 8, 5) = 1 countfactor( 8, 6) = 3 countfactor( 8, 7) = 2 countfactor( 8, 8) = 3 countfactor( 8, 9) = 3 countfactor( 8,10) = 2 countfactor( 8,11) = 1 countfactor( 8,12) = 3 countfactor( 9, 1) = 5 countfactor( 9, 2) = 8 countfactor( 9, 3) = 3 countfactor( 9, 4) = 5 countfactor(10, 1) = 1 countfactor(10, 2) = 1 countfactor(10, 3) = 1 countfactor(10, 4) = 1 countfactor(10, 5) = 1 countfactor(10, 6) = 1 countfactor(10, 7) = 1 countfactor(10, 8) = 1 countfactor(10, 9) = 1 countfactor(10,10) = 1 countfactor(10,11) = 1 countfactor(10,12) = 1 countfactor(10,13) = 1 countfactor(10,14) = 1 countfactor(10,15) = 1 countfactor(10,16) = 1 countfactor(10,17) = 1 countfactor(10,18) = 1 countfactor(10,19) = 1 countfactor(10,20) = 1 countfactor(10,21) = 1 countfactor(10,22) = 1 countfactor(10,23) = 1 countfactor(10,24) = 1 countfactor(11, 1) = 101 countfactor(12, 1) = 36 countfactor(12, 2) = 25 !-------- ! RETURN END subroutine getcountscoul ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine getcountsfeyn(countfactor) ! use beowulf_helpers implicit none ! Out: integer :: countfactor(maxgraphs,maxmaps) ! ! This subroutine contains information on how many points to devote ! to each of the ten graphs and, within each graph, how many to ! devote to each of the maps from the random numbers x to the loop ! momenta. These are for Feynman gauge. ! ! Latest revision 8 December 2002. ! integer :: graphnumber,mapnumber ! DO graphnumber = 1,maxgraphs DO mapnumber = 1,maxmaps countfactor(graphnumber,mapnumber) = 20 END DO END DO ! !-------- countfactor( 1, 1) = 7 countfactor( 1, 2) = 10 countfactor( 2, 1) = 2 countfactor( 2, 2) = 2 countfactor( 3, 1) = 13 countfactor( 3, 2) = 5 countfactor( 3, 3) = 7 countfactor( 3, 4) = 6 countfactor( 3, 5) = 7 countfactor( 3, 6) = 11 countfactor( 4, 1) = 4 countfactor( 4, 2) = 2 countfactor( 4, 3) = 3 countfactor( 4, 4) = 3 countfactor( 5, 1) = 5 countfactor( 5, 2) = 5 countfactor( 5, 3) = 9 countfactor( 5, 4) = 2 countfactor( 5, 5) = 4 countfactor( 5, 6) = 2 countfactor( 6, 1) = 5 countfactor( 6, 2) = 2 countfactor( 6, 3) = 4 countfactor( 6, 4) = 3 countfactor( 6, 5) = 7 countfactor( 6, 6) = 7 countfactor( 7, 1) = 17 countfactor( 7, 2) = 7 countfactor( 7, 3) = 21 countfactor( 7, 4) = 6 countfactor( 7, 5) = 11 countfactor( 7, 6) = 3 countfactor( 7, 7) = 30 countfactor( 7, 8) = 15 countfactor( 7, 9) = 11 countfactor( 7,10) = 17 countfactor( 7,11) = 11 countfactor( 7,12) = 21 countfactor( 7,13) = 4 countfactor( 7,14) = 11 countfactor( 7,15) = 3 countfactor( 7,16) = 47 countfactor( 7,17) = 29 countfactor( 7,18) = 9 countfactor( 8, 1) = 5 countfactor( 8, 2) = 5 countfactor( 8, 3) = 6 countfactor( 8, 4) = 2 countfactor( 8, 5) = 8 countfactor( 8, 6) = 5 countfactor( 8, 7) = 3 countfactor( 8, 8) = 5 countfactor( 8, 9) = 6 countfactor( 8,10) = 3 countfactor( 8,11) = 10 countfactor( 8,12) = 6 countfactor( 9, 1) = 1 countfactor( 9, 2) = 1 countfactor( 9, 3) = 1 countfactor( 9, 4) = 1 countfactor(10, 1) = 1 countfactor(10, 2) = 1 countfactor(10, 3) = 1 countfactor(10, 4) = 1 countfactor(10, 5) = 1 countfactor(10, 6) = 4 countfactor(10, 7) = 1 countfactor(10, 8) = 1 countfactor(10, 9) = 1 countfactor(10,10) = 1 countfactor(10,11) = 1 countfactor(10,12) = 1 countfactor(10,13) = 1 countfactor(10,14) = 1 countfactor(10,15) = 1 countfactor(10,16) = 2 countfactor(10,17) = 1 countfactor(10,18) = 4 countfactor(10,19) = 1 countfactor(10,20) = 1 countfactor(10,21) = 1 countfactor(10,22) = 1 countfactor(10,23) = 1 countfactor(10,24) = 2 countfactor(11, 1) = 47 countfactor(12, 1) = 105 countfactor(12, 2) = 77 !-------- ! RETURN END subroutine getcountsfeyn ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine getcountsshwr(countfactor) ! use beowulf_helpers implicit none ! Out: integer :: countfactor(maxgraphs,maxmaps) ! ! This subroutine contains information on how many points to devote ! to each of the ten graphs and, within each graph, how many to ! devote to each of the maps from the random numbers x to the loop ! momenta. These are for Coulomb gauge with showers. ! ! Latest revision 6 April 2002. ! integer :: graphnumber,mapnumber ! DO graphnumber = 1,maxgraphs DO mapnumber = 1,maxmaps countfactor(graphnumber,mapnumber) = 20 END DO END DO ! !-------- countfactor( 1, 1) = 3 countfactor( 1, 2) = 3 countfactor( 2, 1) = 7 countfactor( 2, 2) = 4 countfactor( 3, 1) = 3 countfactor( 3, 2) = 1 countfactor( 3, 3) = 3 countfactor( 3, 4) = 1 countfactor( 3, 5) = 1 countfactor( 3, 6) = 23 countfactor( 4, 1) = 1 countfactor( 4, 2) = 1 countfactor( 4, 3) = 1 countfactor( 4, 4) = 1 countfactor( 5, 1) = 8 countfactor( 5, 2) = 1 countfactor( 5, 3) = 4 countfactor( 5, 4) = 1 countfactor( 5, 5) = 2 countfactor( 5, 6) = 1 countfactor( 6, 1) = 4 countfactor( 6, 2) = 2 countfactor( 6, 3) = 1 countfactor( 6, 4) = 1 countfactor( 6, 5) = 4 countfactor( 6, 6) = 5 countfactor( 7, 1) = 2 countfactor( 7, 2) = 1 countfactor( 7, 3) = 9 countfactor( 7, 4) = 1 countfactor( 7, 5) = 1 countfactor( 7, 6) = 1 countfactor( 7, 7) = 26 countfactor( 7, 8) = 1 countfactor( 7, 9) = 1 countfactor( 7,10) = 2 countfactor( 7,11) = 1 countfactor( 7,12) = 5 countfactor( 7,13) = 1 countfactor( 7,14) = 1 countfactor( 7,15) = 1 countfactor( 7,16) = 3 countfactor( 7,17) = 2 countfactor( 7,18) = 2 countfactor( 8, 1) = 1 countfactor( 8, 2) = 1 countfactor( 8, 3) = 1 countfactor( 8, 4) = 1 countfactor( 8, 5) = 1 countfactor( 8, 6) = 1 countfactor( 8, 7) = 1 countfactor( 8, 8) = 1 countfactor( 8, 9) = 1 countfactor( 8,10) = 1 countfactor( 8,11) = 1 countfactor( 8,12) = 1 countfactor( 9, 1) = 1 countfactor( 9, 2) = 1 countfactor( 9, 3) = 1 countfactor( 9, 4) = 1 countfactor(10, 1) = 1 countfactor(10, 2) = 1 countfactor(10, 3) = 1 countfactor(10, 4) = 1 countfactor(10, 5) = 1 countfactor(10, 6) = 1 countfactor(10, 7) = 1 countfactor(10, 8) = 1 countfactor(10, 9) = 1 countfactor(10,10) = 1 countfactor(10,11) = 1 countfactor(10,12) = 1 countfactor(10,13) = 1 countfactor(10,14) = 1 countfactor(10,15) = 1 countfactor(10,16) = 1 countfactor(10,17) = 1 countfactor(10,18) = 1 countfactor(10,19) = 1 countfactor(10,20) = 1 countfactor(10,21) = 1 countfactor(10,22) = 1 countfactor(10,23) = 1 countfactor(10,24) = 1 countfactor(11, 1) = 77 countfactor(12, 1) = 4 countfactor(12, 2) = 39 !-------- ! RETURN END subroutine getcountsshwr ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890