!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