!                     ------------------------
!                     beowulf.f90  Version 4.0
!                     ------------------------
!
! 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.f90. In particular, a
! user may very well want to modify parameter settings in subroutines
! beowulfinit and beowulfconclude 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.f90.
!
! These routines are designed to be used with a suitably modified version of
! Pythia, which supplies showering beyond the "primary" level and supplies
! hadronization. For this purpose, you need a main program with calls to
! Pythia routines. The main program calls the Pythia routine PYINIT, which will
! call the user the process subroutine upinit given below. The main program also
! calls the Pythia routine PYEVNT, which calls the user process subroutine 
! upevnt given below. At the user's discretion, a standard set of event analysis 
! can be performed with a call to hrothgarII(0,'newresult '), given below, for 
! each event.
!
! To use the beowulf routines without Pythia, you need
!
! program main
! call beowulf
! end
!
! subroutine PREP_EVENT  ! Not needed here but called in subroutine upevnt.
! return
! end
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!                        Modules for program beowulf
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
module beowulf_revisiondate
  character(len=36), parameter ::   &         ! date for beowulf.vnn.f90
          revisiondate = 'Latest revision 4 September 2005    '
!                        '------------------------------------'
end module beowulf_revisiondate
!
!----------------
!
module beowulf_parameters
!
  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 = 512
  integer, parameter :: eventrecordsize = 1000
!
end module beowulf_parameters
!
!----------------
!
module beowulf_structures
!
use beowulf_parameters
!
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 graphstructure
  integer                          :: graphnumber
  integer                          :: order
  integer, dimension(0:3*size-1,2) :: vrtx
  integer, dimension(2*size,3)     :: prop
end type graphstructure
!
type flavorstructure
  integer                               :: setnumber 
  character(len=5), dimension(3*size-1) :: type    ! parton type of propagator
  integer, dimension(3*size-1)          :: flavor  ! 0, +1, -1, +2, ...
end type flavorstructure
!
type cutstructure
  integer                        :: cutnumber ! identifies the cut
  integer                        :: ncut      ! number of cut propagators
  integer, dimension(size+1)     :: cutindex  ! index of cut propagator
  integer, dimension(size+1)     :: cutsign   ! direction of cut propagator
  logical                        :: leftloop  ! there is a loop to left
  logical                        :: rightloop ! there is a loop to right
  integer                        :: ninloop   ! number of props in loop
  integer, dimension(size+1)     :: loopindex ! indices around loop
  integer, dimension(size+1)     :: loopsign  ! propagator directions
end type cutstructure
!
type softinformation
  character(len=5) :: type0   ! quark, qbar, gluon
  character(len=5) :: type1   ! quark, qbar, gluon
  character(len=5) :: type2   ! quark, qbar, gluon
  real(kind=dbl) :: cos1
  real(kind=dbl) :: cos2
  real(kind=dbl) :: msoftsq
end type softinformation
!
type parton
  character(len=5) :: type   ! quark, qbar, gluon
  integer :: flavor          ! 0, +1, -1, +2, -2, etc
  integer :: self            ! index of this parton
  integer :: parent          ! index of parton's parent
  integer :: ancestor        ! starting parton
  integer :: stringquark     ! index of string attached to quark index
  integer :: stringqbar      ! index of string attached to to qbar index
  integer :: child1          ! index of parton's first daughter
  integer :: child2          ! index of parton's second daughter
  logical :: childless       ! no children (yet)
  logical :: done            ! .true. means don't split this parton
  real(kind=dbl), dimension(3) :: momentum 
  real(kind=dbl) :: kappasq  ! scale for subsequent branchings
end type parton
! The parton types are quark, qbar, gluon. For type quark, the flavors
! are 1, 2, ..., nf in the order up, down, strange, charm, bottom, top.
! For type qbar, the flavors are -1, -2, ..., -nf. For type gluon the
! flavor is 0.
!
type showerlist
  integer :: length                    ! length of the shower list
  integer :: nstrings                  ! number of strings joining partons
  real(kind=dbl) :: rts0               ! sqrt(s) of the starting graph
  real(kind=dbl) :: onemthrust         ! 1 - thrust of the starting graph
  real(kind=dbl) :: msoftsq            ! the soft scale
  real(kind=dbl) :: multfactor         ! factor to multiply |matrix element|^2
  integer :: pii                       ! parton emitting soft gluon
  integer :: pjj                       ! partons absorbing soft gluon
  type(parton), dimension(maxparticles) :: ptn
end type showerlist
!
type renoresults
  real(kind=dbl) :: r                 ! sum Re(integral)
  real(kind=dbl) :: i                 ! sum Im(integral)
  real(kind=dbl) :: bis               ! sum Re(integral) over extra points
  real(kind=dbl) :: chkr              ! sum Re(check_integral)
  real(kind=dbl) :: chki              ! sum Im(check_integral)
  real(kind=dbl) :: rsq               ! sum (Re(integral))**2
  real(kind=dbl) :: isq               ! sum (Im(integral))**2
  real(kind=dbl) :: bissq             ! sum (Re(integral))**2 over extra points
  real(kind=dbl) :: chkrsq            ! sum (Re(check_integral))**2
  real(kind=dbl) :: chkisq            ! sum (Im(check_integral))**2
  integer        :: included          ! number of points included in sums
  integer        :: extracount        ! number of extra points
  integer        :: omitted           ! number of points omitted entirely
  integer,dimension(-9:6) :: nvalpt   ! number pts with values in range n  
  real(kind=dbl) :: valptmax          ! max contribution from a point
  real(kind=dbl),dimension(maxgraphs,maxmaps) :: fluct
                                      ! sum integrand**2 from map,graph   
  type(badpointstate) :: badpointinfo ! random number state for worst point
end type renoresults
!
end module beowulf_structures
!
!----------------
!
module pythiawulf
use beowulf_parameters
!
type partonmod
  integer :: flavor          ! 0, +1, -1, +2, -2, etc
  integer :: self            ! index of this parton
  integer :: parent          ! index of parton's parent
  integer :: stringquark     ! index of string attached to quark index
  integer :: stringqbar      ! index of string attached to to qbar index
  integer :: child1          ! index of parton's first daughter
  integer :: child2          ! index of parton's second daughter
  logical :: childless       ! no children
  logical :: done            ! .true. means don't split this parton
  real(kind=dbl), dimension(0:3) :: momentum 
  real(kind=dbl) :: masssq   ! p**2 
  real(kind=dbl) :: qsqscale ! virtuality of suitable ancestor
  real(kind=dbl) :: costheta ! cosine of angle between daughters
end type partonmod
!
type eventdata
  integer        :: length        ! number of entries in parton list
  real(kind=dbl) :: weight        ! revised weight after event rejections
  real(kind=dbl) :: ktmax         ! largest kt for combining two jets
  integer        :: nparticles    ! number of final state particles
  type(partonmod), dimension(maxparticles) :: ptn
end type eventdata
!
type(eventdata), dimension(eventrecordsize) :: event
integer :: eventindex
!
end module pythiawulf
!
! Note that for many compilers, the modules have to come before the 
! programs in which they are used.
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine beowulf            
!
! To use this as the main program, you need
!
! program main
! call beowulf
! end
!
! subroutine PREP_EVENT
! return
! end
!
! 27 February 2004
!
use beowulf_parameters
use beowulf_structures
implicit none
!
! For communication with Pythia
logical :: pythia
character(len=10):: filename
common /pythiainfo/ pythia,filename
! Timing and counting variables for reno
real(kind=dbl) :: timelimit
common /maxtime/ timelimit
real(kind=dbl) :: deltatime, cputime
integer :: nreno
! reno results variables
type(renoresults) :: sums
!
!
pythia = .false.
call beowulfinit        ! Reads user inputs, writes run information.
call renoinit(sums)     ! Sets all of the sums to zero.
cputime = 0.0
nreno = 0
call timing(deltatime)
!
! -----------------------Loop for calls to reno--------------------------
!
DO WHILE (cputime.LT.timelimit)
nreno = nreno + 1
!                 ! Generates one group of points.
call reno(sums)   ! Nominal standard results and diagnostics are in "sums."
                  ! Events sent to hrothgar for accumulation and analysis.
!
call timing(deltatime)
cputime = cputime + deltatime
!
END DO   ! while (cputime.LT.timelimit)
!
! -----------------------------------------------------------------------
!
call beowulfconclude(sums,nreno,cputime)  ! Writes results and diagnostics.
!
return
END subroutine beowulf
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine beowulfinit
!
use beowulf_parameters
use beowulf_revisiondate
use beowulf_structures
implicit none
!
! 3 December 2003
! 2 October 2004 latest modification
!
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=8) :: mode
common /programmode/ mode
! For communication with Pythia
logical :: pythia
character(len=10):: filename
common /pythiainfo/ pythia,filename
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,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
! Factor in msoftsq =(msoftfactor*rts0*onemthrust)**2
real(kind=dbl) :: msoftfactor
common /softcutoff/ msoftfactor
! Factor for relating kappasq to onemthrust for small onemthrust
real(kind=dbl) :: onemthrustcrit
common /thrustfactor/ onemthrustcrit
! 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
! A blank shower for hrothgar:
type(showerlist) :: blankshower
! Seed for random number generator.
integer :: seed
! Today's date:
character(len=28) :: date
! Loop control variables:
integer :: n,graphnumber,mapnumber
! Logical variables that control aspects of the output:
logical :: allgraphs,writeprograminfo
! A temporary real number:
real(kind=dbl) :: temp 
! A temporary integer:
integer :: itemp
! alpha_s/pi
real(kind=dbl) :: alpi
real (kind=dbl), parameter :: pi = 3.141592653589793239d0
! Function for calculating thrustcut
real(kind=dbl) :: thrustlimit
!
!-----------------------------------------------------------------------
!
! 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 and number of flavors.
nc = 3.0d0
nf = 5.0d0
cf = (nc**2 - 1)/2.0d0/nc
! Physics data.
mz = 91.1876d0
widthz = 2.4952d0
alphasofmz = 0.118d0
sinsqthetaw = 0.23113
externalrts = mz
! Initialize common /GRAPHCOUNTS/.
call graphcountinit
! Parameters for subroutine CALCULATE for throwing away bad points.
badnesslimit = 1.0d4
cancellimit = 1.0d4
! Factor in msoftsq =(msoftfactor*rts0*onemthrust)**2
msoftfactor = 0.333333333333333333d0
! Factor for relating kappasq to onemthrust for small onemthrust
onemthrustcrit = 0.05d0
! 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 = 1.0d0/externalrts**2
! 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.
allgraphs = .true.
!
IF (pythia) THEN
!  mode = 'born    '
!  showerendratio = 0.0d0/externalrts**2  ! Not used.
  mode = 'showerI '
  showerendratio = 1.0d0/externalrts**2  ! Applied to showerI showers.
!   mode = 'showerII'        ! An alternative.
!   showerendratio = 1.0d0/externalrts**2  ! For showerII, use Pythia default.
ELSE
  write(nout,*) 'Please give program mode, '
  write(nout,*) '(born, nlo, hocoef, showr0I, showr0II, showerI, showerII).'
  read(nin,98)mode
  98    format(a8)
END IF
!
! 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.
  allgraphs = .false.
ELSE IF (mode.EQ.'hocoef  ') THEN
  usegraph(11) = .false.
  usegraph(12) = .false.
  allgraphs = .false.
ELSE IF (mode.EQ.'nlo     ') THEN
  continue
ELSE IF (mode.EQ.'showerI ') THEN
  continue
ELSE IF (mode.EQ.'showerII') THEN
  continue
ELSE IF (mode.EQ.'showr0I ') 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.
  allgraphs = .false.
  mode = 'showerI '
ELSE IF (mode.EQ.'showr0II') 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.
  allgraphs = .false.
  mode = 'showerII'
ELSE
  write(nout,*)'Not prepared for mode ',mode,'.'
  stop
END IF
! So from here on, the possible modes are
! born, nlo, hocoef, showerI, showerII
!
! Use this if you want to write a little file in the form
! Number of events written was      188
! Number nreno of main loop executions was        7
! IF (pythia) THEN
!     write(nout,*)'Please give file name. (Up to 10 characters, no suffix.)'
!     read(nin,*)filename
! END IF
filename = 'nofile    '  ! With this file name, no file is written.
!
!
IF (pythia) THEN
  temp = 1.0d0
ELSE
  temp = 1.0d0
!  write(nout,*)'Please give alpha_s/0.118'
!  read(nin,*)temp
END IF
alphasofmz = 0.118d0*temp
!
IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN
  gauge = 'coulomb'
ELSE
  write(nout,*) &
   'Please give gauge choice (coulomb or feynman).'
  read(nin,99)gauge
99    format(a7)
END IF
!
! Parameter for hrothgar for not sending events to Pythia. 
! NLO cross section for thrust < thrustcut will be 0.92 * sigma_T. 
! This makes the total cross section as generated about right.
thrustcut = thrustlimit(0.92d0)
write(*,*)'For Pythia, thrustcut = ',thrustcut
!
! Here we set the count factors COUNTFACTOR(GRAPHNUMBER,MAPNUMBER).
! These, multiplied by NSETS, give GROUPSIZE(GRAPHNUMBER,MAPNUMBER).
!
IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) 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.
!
call hrothgar(blankshower,1.0d0,1,'startcalc ')
!
write(nout,*)'Please give the approximate cpu time limit (hours).'
read(nin,*) hourlimit
timelimit = hourlimit * 3600
IF (hourlimit.LT.0.9d0) THEN
  nsets = 1 + nint(10.0d0*hourlimit)
ELSE
  nsets = 10
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.16667d0 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 4.0 ',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,sinsqthetaw,mz,widthz,externalrts
106     format('alpha_s(mz) = ',f7.4,' and sin^2(theta_W) = ',f7.4,/, &
               'with mz = ',f7.4,' and Zwidth = ',f7.4,'.',/, &
               'sqrt(s) = ',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,sinsqthetaw,mz,widthz,externalrts
ELSE IF ((mode.EQ.'showerI ').AND.allgraphs) THEN
  write(nout,108)
108     format('observable/sigma_0 ', &
         'at next-to-leading order with primary showers.')
  write(nout,106)alphasofmz,sinsqthetaw,mz,widthz,externalrts
ELSE IF ((mode.EQ.'showerII').AND.allgraphs) THEN
  write(nout,109)
109     format('observable/sigma_0 ', &
         'at next-to-leading order',/,'with primary and secondary showers.')
  write(nout,106)alphasofmz,sinsqthetaw,mz,widthz,externalrts
ELSE IF ((mode.EQ.'showerI ').AND.(.NOT.allgraphs)) THEN
  write(nout,110)
110     format('observable/sigma_0 with primary showers.')
  write(nout,106)alphasofmz,sinsqthetaw,mz,widthz,externalrts
ELSE IF (mode.EQ.'showerII'.AND.(.NOT.allgraphs)) THEN
  write(nout,111)
111     format('observable/sigma_0 with primary and secondary showers.')
  write(nout,106)alphasofmz,sinsqthetaw,mz,widthz,externalrts
ELSE IF (mode.EQ.'hocoef  ') THEN
  write(nout,112)
112     format('the coefficient of (alpha_s/pi)^2 ', &
         'for observable/sigma_0.')
ELSE
  write(nout,*)'Not prepared for mode ',mode,'.'
  stop
END IF
!
IF ((pythia).AND.(filename.NE.'nofile    ')) THEN
  write(nout,114)adjustr(filename)//'.info'
114     format('Event count information will be in file ',A15,'.')
END IF
!
IF (pythia) THEN
  write(nout,115)
115    format('Events will be sent to Pythia; report by HrothgarII.',/)
END IF
!
write(nout,121)nc,nf
121   format('beowulf will use',1p g9.2,'colors and',1p g9.2,'flavors.')
write(nout,122)
122   format('Renormalization parameter:')
write(nout,123)muoverrts
123   format('                   mu /sqrt(s):',1p g11.4)
write(nout,124)pi*alpi(muoverrts*externalrts)
124   format('                   alpha_s(mu):',1p g11.4)
write(nout,125)
125   format('Cutoff parameters:')
write(nout,126)badnesslimit,cancellimit
126   format('                 badness limit:',1p g9.2,/, &
       '            cancellation limit:',1p g9.2,/)
IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) 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
  IF ((pythia).AND.(thrustcut.LT.1.0d0)) THEN
    write(nout,754)thrustcut
754 format('Events not sent to Pythia when partonic thrust is greater than', &
            f8.4)
  END IF
END IF
IF (mode.EQ.'showerII') THEN
    write(nout,752) sqrt(showerendratio*externalrts**2)
752 format('Secondary showering stops at sqrt(virtuality) =', &
            f7.2,' GeV approximately.')
END IF
!
IF (allgraphs) THEN
  write(nout,116)
116     format('All graphs are used.')
  ELSE
  write(nout,117)
117     format('The following graphs are used.',/, &
               '(1-10 are alpha_s^2, 11 and 12 are alpha_s^1.)')
  DO n = 1,maxgraphs
    write(nout,118) n,usegraph(n)
118       format('Use of graph ',i2,': ',l1)
  END DO
END IF
!
writeprograminfo = .false.
IF (writeprograminfo) THEN
 write(nout,133)
133    format(/,'Points per group:',/)
 DO graphnumber = 1,numberofgraphs
   IF (usegraph(graphnumber)) THEN
     DO mapnumber = 1,numberofmaps(graphnumber)
       write(nout,134) graphnumber,mapnumber, &
                    groupsize(graphnumber,mapnumber)
134          format('Graph',i3,' map',i3,' points used',i6)
     END DO
   write(nout,135)groupsizegraph(graphnumber)
135      format(/,'Total points for this graph',i8,/)
   END IF
 END DO
 write(nout,136)groupsizetotal
136    format('Total points in each group:',i8,/)
 write(nout,137)
137    format('Parameters for contour deformation:')
 write(nout,138)deformalpha,deformbeta,deformgamma
138    format('                         alpha:',1p g9.2,/, &
        '                         beta :',1p g9.2,/, &
        '                         gamma:',1p g9.2)
 write(nout,139)
139    format('Parameters for smearing in rts:')
 write(nout,140)smearfctr,lowpwr,highpwr
140    format('       scale factor for rts/e0:',1p g9.2,/, &
        '               small rts power:  ',i3,/, &
        '               large rts power:  ',i3)
END IF
!
write(nout,*)' '
!
end subroutine beowulfinit
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine beowulfconclude(sums,nreno,cputime)
use beowulf_parameters
use beowulf_structures
implicit none
! In:
type(renoresults) :: sums
integer :: nreno
real(kind=dbl) :: cputime
!
! Write results and diagnostic information.
! 16 November 2003
!
integer :: countfactor(maxgraphs,maxmaps)
integer :: groupsize(maxgraphs,maxmaps)
integer :: groupsizegraph(maxgraphs)
integer :: groupsizetotal
common /montecarlo/groupsize,groupsizegraph,groupsizetotal
integer :: graphnumber,mapnumber
integer :: itemp,n
real(kind=dbl) :: temp
logical :: writefluct, writediagnostic
! Today's date:
character(len=28) :: date
! Errors for results variables:
real(kind=dbl) :: errorr,errori,errorbis,errorchkr,errorchki
! Average of FLUCT variables:
real(kind=dbl) :: avfluct
! What the program should do
character(len=8) :: mode
common /programmode/ mode
! How many graphs and how many maps for each:
integer :: numberofgraphs
integer :: numberofmaps(maxgraphs)   
common /graphcounts/ numberofgraphs,numberofmaps
! For communication with Pythia
logical :: pythia
character(len=10):: filename
common /pythiainfo/ pythia,filename
! A blank shower for hrothgar:
type(showerlist) :: blankshower
!
! Calculate the results for the standard integrals and their errors.
!
sums%r = sums%r/nreno
sums%i = sums%i/nreno
sums%bis = sums%bis/nreno
sums%chkr = sums%chkr/nreno
sums%chki = sums%chki/nreno
!
sums%rsq = sums%rsq/nreno
sums%isq = sums%isq/nreno
sums%bissq = sums%bissq/nreno
sums%chkrsq = sums%chkrsq/nreno
sums%chkisq = sums%chkisq/nreno
!
IF (nreno.EQ.1) THEN
  write(nout,*)'nreno = 1 changed to 2 to avoid 1/0.'
  write(nout,*)'Results will be finite but wrong.'
  write(nout,*)' '
  nreno = 2
END IF
errorr = sqrt((sums%rsq - sums%r**2)/(nreno - 1))
errori = sqrt((sums%isq - sums%i**2)/(nreno - 1))
errorbis = sqrt((sums%bissq - sums%bis**2)/(nreno - 1))
errorchkr = sqrt((sums%chkrsq - sums%chkr**2)/(nreno-1))
errorchki = sqrt((sums%chkisq - sums%chki**2)/(nreno-1))
!
DO graphnumber = 1,numberofgraphs
DO mapnumber = 1,numberofmaps(graphnumber)
  sums%fluct(graphnumber,mapnumber) = &
         sums%fluct(graphnumber,mapnumber)/nreno
END DO
END DO
! -----------------------------------------------------------------------
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)sums%included
write(nout,801)sums%extracount
write(nout,802)sums%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)sums%r,errorr
601   format('         re(result) =',1p g13.5,' +/-',1p g10.2)
write(nout,602)sums%i,errori
602   format('         im(result) =',1p g13.5,' +/-',1p g10.2)
write(nout,603)sums%bis,errorbis
603   format('  cutoff correction =',1p g13.5,' +/-',1p g10.2)
!
! Check integral.
!
IF ((mode.NE.'showerI ').AND.(mode.NE.'showerII')) 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)sums%chkr,errorchkr
  605   format('          re(check) =',1p g13.5,' +/-',1p g10.2)
  write(nout,606)sums%chki,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(blankshower,1.0d0,nreno,'calcdone  ')
IF (pythia) THEN
   call hrothgarII(nreno,'calcdone  ')
END IF
!
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,sums%nvalpt(n)
702    format('Number of points with ',i2,' < log_10(|v|) <', i2, &
        ' is ',i10)
END DO
!
write(nout,*)' '
write(nout,703)sums%valptmax
703   format('Biggest contribution was',1p g12.3)
write(nout,704)sums%badpointinfo%graphnumber,sums%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  &
       + sums%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), &
                   sums%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 = sums%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(sums%badpointinfo)
END IF
!
!----------------------------------------------
! 
call daytime(date)
write(nout,707)date
707   format(/ &
    ,'Done ',a,/)
!
!---
!
END subroutine beowulfconclude
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine upinit
use beowulf_parameters
use beowulf_structures
implicit none
!
! Puts initialization information in Pythia common block. Also gets beowulf
! started.
!
! 29 November 2003
!  7 April 2005
!
!------------------Pythia common blocks-------------------
!
integer, parameter :: MAXPUP = 100
integer :: IDBMUP(2), PDFGUP(2), PDFSUP(2), IDWTUP, NPRUP, LPRUP(MAXPUP)
real(kind=dbl) :: EBMUP(2), XSECUP(MAXPUP), XERRUP(MAXPUP), XMAXUP(MAXPUP)
common /HEPRUP/ IDBMUP, EBMUP, PDFGUP, PDFSUP, IDWTUP, NPRUP, XSECUP,     &
                XERRUP, XMAXUP, LPRUP
!
INTEGER :: ISHWVAR
COMMON/SHWPSS/ISHWVAR
!
INTEGER :: TYPE,ANGL,MONO,RECO
COMMON/CLSTR/TYPE,ANGL,MONO,RECO
!
save /HEPRUP/
save /SHWPSS/
save /CLSTR/
!
!--------------------------------------------------------

logical :: pythia
character(len=10):: filename
common /pythiainfo/ pythia,filename ! For pythia information
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts
!
integer j,im
!
pythia = .true.      ! This changes the behavior of some beowulf routines.
call beowulfinit
!
! The information for Pythia.
!
IDBMUP(1) = 11    ! e-
IDBMUP(2) = -11   ! e+
EBMUP(1) = 0.5d0*externalrts  ! beam1 energy
EBMUP(2) = 0.5d0*externalrts  ! beam2 energy
PDFGUP(1) = -1  ! No PDFs used
PDFGUP(2) = -1  ! No PDFs used
PDFSUP(1) = -1  ! No PDFs used
PDFSUP(2) = -1  ! No PDFs used
IDWTUP = -4     ! input events +/- weighted, output events +/- weighted
NPRUP = 1       ! only one user process
DO j = 1,MAXPUP
 XSECUP(j) = 0.0d0  ! we don't need a cross section
 XERRUP(j) = 0.0d0  ! nor do we need an error on the cross section
 XMAXUP(j) = 0.0d0  ! max weight is not needed
 LPRUP(j) = 0       ! user process number for the unused user processes
END DO
LPRUP(1) = 661      ! user process number for the one process is 661.
!
ISHWVAR=12
im=1111
TYPE=MOD(im/1000,10)
ANGL=MOD(im/100 ,10)
MONO=MOD(im/10  ,10)
RECO=MOD(im     ,10)
!
RETURN
END subroutine upinit
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine upevnt
!
use beowulf_parameters
use beowulf_revisiondate
use beowulf_structures
use pythiawulf
! Defines type(eventdata), dimension(eventrecordsize) :: event
!         integer :: eventindex
! Data for event and eventindex are created in hrothgar.
implicit none
!
! Puts events in Pythia common block. When out of events, gets more using
! subroutine reno. A synthetic shower history is produced by subroutine
! translate and is available to upevnt, which passes the information on
! to Pythia by means of the common blocks /IPASS/ and /HEPEUP/.
!
! For common /IPASS/, upevnt supplies
! KTLIMIT = max kt of any splitting in synthetic shower history
! SCALEQ(n) = sqrt(virtuality) at which parton n was ``created''
! THORD(n) = angle between parton n and its sister in splitting of mother
!
! For common /HEPEUP/, upevnt supplies
! NUP = number of entries in event record = event(j)%length + 3 
! IDPRUP = 661                            ! Our ID of the process.
! XWGTUP = event(j)%weight                ! Weight for this event.
! SCALUP = muoverrts*externalrts          ! Scale, maybe used in alpha_s.
! AQEDUP = 0.0078125d0                    ! alpha_QED, = 1/128
! AQCDUP = pi*alpi(muoverrts*externalrts) ! alpha_QCD
! IDUP(n) = flavor code for parton n, eg. 21 for gluon
! ISTUP(n) = 1     ! outgoing final state particle, splitting allowed.
!          = 7     ! outgoing final state particle, no splitting.
!          = 2     ! intermediate resonance, preserve mass
! MOTHUP(1,n) = MOTHUP(2,n) = parent of parton n
! ICOLUP(1,n) = index of quark string for parton n
! ICOLUP(2,n) = index of qbar string for parton n
! PUP(mu,n) = momentum of parton n, mu = 1,...,4
! PUP(5,n) = p**2 of parton n
! VTIMUP(n) = 0.0d0             ! lifetime c*tau
! SPINUP(n) = 9.0d0             ! spin (code for unknown spin)
!
! 19 December 2003
!  6 June 2005
!
!------------------Pythia common block-------------------
!
integer, parameter :: MAXNUP = 500
integer :: NUP, IDPRUP, IDUP(MAXNUP), ISTUP(MAXNUP),  &
           MOTHUP(2,MAXNUP), ICOLUP(2,MAXNUP)
real(kind=dbl) :: XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP(5,MAXNUP), &
                  VTIMUP(MAXNUP), SPINUP(MAXNUP)
common /HEPEUP/ NUP, IDPRUP, XWGTUP, SCALUP, AQEDUP, AQCDUP, IDUP,   &
                ISTUP, MOTHUP, ICOLUP, PUP, VTIMUP, SPINUP
save /HEPEUP/
!
!---Common block for extra shower information for Pythia---
!
integer :: NMATCH
real(kind=dbl) :: Q0RES,QXRES,QHARD,PASSKT(20),SCALEQ(20),THORD(20)
real(kind=dbl) :: KTLIMIT
common /IPASS/ Q0RES,QXRES,QHARD,PASSKT,SCALEQ,THORD,KTLIMIT,NMATCH
save /IPASS/
!
! ------- beowulf counting variable nreno -------
!
integer :: beowulfnreno
common /beowulfcounting/ beowulfnreno
save /beowulfcounting/
!
! -------
real(kind=dbl) :: muoverrts
common /renormalize/ muoverrts
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts
!
logical :: pythia
character(len=10):: filename
common /pythiainfo/ pythia,filename ! For pythia information
! Timing and counting variables for reno
real(kind=dbl) :: timelimit
common /maxtime/ timelimit
real(kind=dbl),save :: deltatime, cputime
integer,save :: nreno
integer,save :: eventnumber,j,n,nn       ! j = eventindex
! reno results variables
type(renoresults),save :: sums
!
! ---- for testing and reporting sample events
!
integer :: nparticles
real(kind=dbl) :: kf(maxparticles,0:3)  ! final state momenta
real(kind=dbl) :: y43,y32
real(kind=dbl) :: jetmass(3) 
!
real(kind=dbl) :: alpi
real (kind=dbl), parameter :: pi = 3.141592653589793239d0
!
logical,save :: init = .true.
integer :: ntemp
!
!-----------------------------------------------------------------------
!
IF (init) THEN
init = .false.
cputime = 0.0
nreno = 0
eventnumber = 0
eventindex = 0
call timing(deltatime)
call hrothgarII(0,'startcalc ')
END IF
!
IF (eventindex.EQ.0) THEN           ! Out of events.
  IF (.NOT.init) THEN               ! We were out of events because Pythia
    call hrothgarII(0,'groupdone ') ! finished a whole group of evenets.
  END IF
IF (cputime.LT.timelimit) THEN   ! Still have time.
! Out of events but we have time, so get some more events.
!
ntemp = 0
DO
  nreno = nreno + 1
  ntemp = ntemp + 1
  call hrothgarII(0,'startgroup')
                    ! Generates one group of points.
                    ! Nominal standard results and diagnostics are in "sums."
   call reno(sums)  ! Events sent to hrothgar for accumulation and analysis.
                    ! hrothgar puts event data in array event
                    ! and records eventindex = number of events in group.
  IF (eventindex.GT.0) EXIT
  IF (ntemp.GT.1) THEN
    write(nout,*)'That is strange but not impossible:'
    write(nout,*)'Subroutine reno failed to return any events.'
    IF (ntemp.GT.100) THEN
      write(nout,*)'Too many tries in upevnt'
      STOP
    END IF
  END IF
END DO
call timing(deltatime)
cputime = cputime + deltatime
ELSE  
! Out of events and out of time. So quit.
! Set NUP to zero. Write beowulf results to screen.
! Write nreno and eventnumber = number of accepted events to file:
!
!  Note that in beowulf's calculation one sums over events and divides by nreno.
!
!   <S> = (1/nreno) sum_j w_j S_j
!
!  We now accept events with probability p_j, so that
!
!   <S> = (1/nreno) sum_j (w_j/p_j) \theta(acceptevent_j) S_j
!
!  Thus we record a new weight equal to oldweight/probability. We are
!  still dividing by nreno. So if standard weights are defined by
!
!   <S> = (1/naccepted) sum_j \theta(acceptevent_j) * standardweight_j S_j
!
!  then the standard weights are the beowulf weights times naccepted/nreno.
!
!
  NUP = 0
  beowulfnreno = nreno ! Program main gets this information
  IF (filename.NE.'nofile    ') THEN  ! By default, this is not used.
    open(888,file=trim(filename)//'.info',status='new')
    write(888,1201)eventnumber
1201    format('Number of events written was',i9)
    write(888,1202)nreno
1202    format('Number nreno of main loop executions was',i9)
    close(888)
  END IF
  call beowulfconclude(sums,nreno,cputime)  ! Writes results and diagnostics.
  RETURN
!
END IF  ! (cputime.LT.timelimit)
END IF  ! (eventindex.EQ.0)
!
! There are events left. Report the next one to Pythia.
!
eventnumber = eventnumber + 1
j = eventindex
!
! Put event(eventindex) into Pythia common blocks. Note that the main program
! will call hrothgarII with the event information after Pythia processes it.
!
NUP = event(j)%length + 3               ! Number of particle entries.
IDPRUP = 661                            ! Our ID of the process.
XWGTUP = event(j)%weight                ! Weight for this event.
SCALUP = muoverrts*externalrts          ! Scale, maybe used in alpha_s.
AQEDUP = 0.0078125d0                    ! alpha_QED, = 1/128
AQCDUP = pi*alpi(muoverrts*externalrts) ! alpha_QCD
!
IDUP(1) = 11    ! e-
ISTUP(1) = -1   ! incoming
MOTHUP(1,1) = 0 ! no mother1
MOTHUP(2,1) = 0 ! no mother2
ICOLUP(1,1) = 0 ! no color-q
ICOLUP(2,1) = 0 ! no color-qbar
PUP(1,1) = 0.0d0              ! p_x
PUP(2,1) = 0.0d0              ! p_y
PUP(3,1) = 0.5d0*externalrts  ! p_z
PUP(4,1) = 0.5d0*externalrts  ! E
PUP(5,1) = 0.0d0              ! mass
VTIMUP(1) = 0.0d0             ! lifetime c*tau
SPINUP(1) = 9.0d0             ! spin (code for unknown spin)
!
IDUP(2) = -11   ! e+
ISTUP(2) = -1   ! incoming
MOTHUP(1,2) = 0 ! no mother1
MOTHUP(2,2) = 0 ! no mother2
ICOLUP(1,2) = 0 ! no color-q
ICOLUP(2,2) = 0 ! no color-qbar
PUP(1,2) = 0.0d0              ! p_x
PUP(2,2) = 0.0d0              ! p_y
PUP(3,2) = -0.5d0*externalrts ! p_z
PUP(4,2) = 0.5d0*externalrts  ! E
PUP(5,2) = 0.0d0              ! mass
VTIMUP(2) = 0.0d0             ! lifetime c*tau
SPINUP(2) = 9.0d0             ! spin (code for unknown spin)
!
IDUP(3) = 22    ! gamma
ISTUP(3) = 2    ! intermediate resonance, preserve mass
MOTHUP(1,3) = 1 ! mother1 is the e-
MOTHUP(2,3) = 2 ! mother2 is the e+
ICOLUP(1,3) = 0 ! no color-q
ICOLUP(2,3) = 0 ! no color-qbar
PUP(1,3) = 0.0d0              ! p_x
PUP(2,3) = 0.0d0              ! p_y
PUP(3,3) = 0.0d0              ! p_z
PUP(4,3) = externalrts        ! E
PUP(5,3) = externalrts        ! mass
VTIMUP(3) = 0.0d0             ! lifetime c*tau
SPINUP(3) = 9.0d0             ! spin (code for unknown spin)
!
DO n=1,event(j)%length
  nn = n + 3           ! index for Pythia, including the three particles above.
  IF (event(j)%ptn(n)%flavor.EQ.0) THEN
    IDUP(nn) = 21      ! gluon
  ELSE IF (event(j)%ptn(n)%flavor.EQ.1) THEN
    IDUP(nn) = 2       ! up quark
  ELSE IF (event(j)%ptn(n)%flavor.EQ.-1) THEN
    IDUP(nn) = -2       ! up qbar
  ELSE IF (event(j)%ptn(n)%flavor.EQ.2) THEN
    IDUP(nn) = 1       ! down quark
  ELSE IF (event(j)%ptn(n)%flavor.EQ.-2) THEN
    IDUP(nn) = -1       ! down qbar
  ELSE
    IDUP(nn) = event(j)%ptn(n)%flavor  ! for s (3), c (4), b (5).
  END IF
  IF (event(j)%ptn(n)%childless) THEN
    IF (event(j)%ptn(n)%done) THEN
     ISTUP(nn) = 7     ! outgoing final state particle, no splitting.
!    ISTUP(nn) = 1     ! Needed for pythia version 6221
    ELSE
      ISTUP(nn) = 1     ! outgoing final state particle, splitting allowed.
    END IF
  ELSE
    ISTUP(nn) = 2     ! intermediate resonance, preserve mass
  END IF
  MOTHUP(1,nn) = event(j)%ptn(n)%parent + 3
  MOTHUP(2,nn) = event(j)%ptn(n)%parent + 3
  ICOLUP(1,nn) = event(j)%ptn(n)%stringquark
  ICOLUP(2,nn) = event(j)%ptn(n)%stringqbar
  PUP(1,nn) = event(j)%ptn(n)%momentum(1)  ! p_x
  PUP(2,nn) = event(j)%ptn(n)%momentum(2)  ! p_y
  PUP(3,nn) = event(j)%ptn(n)%momentum(3)  ! p_z
  PUP(4,nn) = event(j)%ptn(n)%momentum(0)  ! E
  PUP(5,nn) = sqrt(event(j)%ptn(n)%masssq)  ! mass
  VTIMUP(nn) = 0.0d0             ! lifetime c*tau
  SPINUP(nn) = 9.0d0             ! spin (code for unknown spin)
  SCALEQ(nn) = sqrt(event(j)%ptn(n)%qsqscale) ! Q scale for splitting
  THORD(nn) = acos(event(j)%ptn(n)%costheta)  ! angle for angular ordering
  PASSKT(nn) = event(j)%ktmax                 ! kT limit for splitting
END DO
!
! Extra information for showers.
!
KTLIMIT = event(j)%ktmax   ! Largest kt value for event, not to be increased.
call prep_event
!
! Write it out so we can see what has happened: change .LE.0 to .LE.1 here
!
IF (eventnumber.LE.0) THEN
      write(nout,870)eventnumber
870     format(/,'Event number       ',i10)
      write(nout,871)event(j)%nparticles
871     format('N final state prtns',i10)
      write(nout,872)event(j)%ktmax
872     format('Value of ktmax     ',f10.4)
      write(nout,873)event(j)%weight
873     format('Weight             ',f10.4)
      DO n=1,event(j)%length
        IF (event(j)%ptn(n)%done) THEN
          IF (event(j)%ptn(n)%childless) THEN
            write(nout,875)n,        &
                        event(j)%ptn(n)%flavor,      &
                        event(j)%ptn(n)%stringquark, &
                        event(j)%ptn(n)%stringqbar,  &
                        event(j)%ptn(n)%parent
875         format(i3,':   flavor',i3, &
                 ' strings ',i3,' ',i3,' mother ',i3,&
                 ' final state, done')
          ELSE
            write(nout,*)'This cannot happen.'
          END IF
        ELSE ! parton is not marked done
          IF (event(j)%ptn(n)%childless) THEN
            write(nout,876)n,        &
                        event(j)%ptn(n)%flavor,      &
                        event(j)%ptn(n)%stringquark, &
                        event(j)%ptn(n)%stringqbar,  &
                        event(j)%ptn(n)%parent
876         format(i3,':   flavor',i3, &
                 ' strings ',i3,' ',i3,' mother ',i3,&
                 ' final state')
          ELSE ! parton is not childless
            write(nout,877)n,        &
                        event(j)%ptn(n)%flavor,      &
                        event(j)%ptn(n)%stringquark, &
                        event(j)%ptn(n)%stringqbar,  &
                        event(j)%ptn(n)%parent
877         format(i3,':   flavor',i3, &
                 ' strings ',i3,' ',i3,' mother ',i3,&
                 ' intermediate')
          END IF
        END IF
        write(nout,878) event(j)%ptn(n)%momentum(0), &
                       event(j)%ptn(n)%momentum(1), &
                       event(j)%ptn(n)%momentum(2), &
                       event(j)%ptn(n)%momentum(3), &
                       sqrt(event(j)%ptn(n)%masssq)
878       format(4 f10.4,' m =',f10.4)
        write(nout,879) sqrt(event(j)%ptn(n)%qsqscale),      &
                        acos(event(j)%ptn(n)%costheta)
879       format(' scaleq =',f10.4,' theta ='f10.4)
      END DO
      call makekfmod(event(j),nparticles,kf) 
      call findjets(kf,nparticles,'durham',y43,y32,jetmass)
      write(*,*) 'Before Pythia, kt3(Durham) =',sqrt(y32)*externalrts
      call findjets(kf,nparticles,'pythia',y43,y32,jetmass)
      write(*,*) 'Before Pythia, kt3(Pythia) =',sqrt(y32)*externalrts
END IF
!
eventindex = eventindex - 1 ! When eventindex = 0, the next time upevnt
                            ! is called, it will ask reno for more events.
RETURN
!
END subroutine upevnt
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
function thrustlimit(v)
!
use beowulf_parameters
implicit none
! In:
real(kind=dbl) :: v
! Out:
real(kind=dbl) :: thrustlimit
!
! This function gives the thrust value such that the integral of 
! d sigma/dt, calculated at order alpha_s^2, from 2/3 to thrustlimit
! is v*sigma_T, the total cross section for electron-positron annihilation
! calculated at order alpha_s^1. Here alpha_s is evaluated at scale
! muoverrts*externalrts.
!
! 6 September 2004
! 18 September 2004 latest change
!
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts
! Global parameters:
real(kind=dbl) :: nc,nf,cf
common /colorfactors/ nc,nf,cf
real(kind=dbl) :: muoverrts
common /renormalize/ muoverrts
!
real(kind=dbl) :: kn0,kn,alpi,alpival,factor
real(kind=dbl) :: integral,x,dx,t
!
alpival = alpi(muoverrts*externalrts)  ! alpha_s/pi
! This factor corrects 1/sigma_0 to 1/sigma_T.
factor = 1.0d0/(1.0d0 + alpival)
!
! kn0(t) and kn(t) are defined by
!
! (1/sigma_0) * (1-t) d sigma/dt
!      = (alpha_s/Pi) KN0(t) + (alpha_s/Pi)^2 KN(t) + ...
!
! Note that these functions are only for NC = 3, NF = 5.
! Note that there is a factor (1-t) to change the integration variable
! from t to x = log(1/(1-t)).
!
IF ( (nc-3.0d0)**2 + (nf-5.0d0)**2 .LT. 0.01d0 ) THEN
  integral = 0.00d0
  x = log(3.0d0)
  dx = 0.001d0
  DO ! Until integral > 1.
    t = 1.0d0 - exp(-x)
    integral = integral + dx*alpival*factor*(kn0(t) + alpival*kn(t))
    IF (integral.GT.v) THEN
      thrustlimit = t
      RETURN
    END IF
    x = x + dx
  END DO
ELSE
  write(nout,*)'(nc,nf) .NE. (3,5) so we calculate function thrustlimit'
  write(nout,*)'at lowest order only.' 
  factor = 3.0d0*cf/4.0d0  ! Here factor corrects the lowest order color factor
  integral = 0.00d0
  x = log(3.0d0)
  dx = 0.001d0
  DO ! Until integral > 1.
    t = 1.0d0 - exp(-x)
    integral = integral + dx*alpival*factor*kn0(t)
    IF (integral.GT.v) THEN
      thrustlimit = t
      RETURN
    END IF
    x = x + dx
  END DO          
END IF
END function thrustlimit
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine hrothgar(theshower,weight,nreno,whattodo)
!
use beowulf_parameters
use beowulf_structures
use pythiawulf   ! For passing list of events, event, to upevnt
implicit none
! In:
type(showerlist) :: theshower
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:
! 04 February 2005
!
character(len=8) :: mode
common /programmode/ mode
logical :: pythia
character(len=10):: filename
common /pythiainfo/ pythia,filename ! For pythia information
!
real(kind=dbl) :: badnesslimit,cancellimit,thrustcut
common /limits/ badnesslimit,cancellimit,thrustcut
integer :: groupsize(maxgraphs,maxmaps)
integer :: groupsizegraph(maxgraphs)
integer :: groupsizetotal
common /montecarlo/groupsize,groupsizegraph,groupsizetotal
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts
!
real(kind=dbl) :: muoverrts
common /renormalize/ muoverrts
real(kind=dbl) :: nc,nf,cf
common /colorfactors/ nc,nf,cf
real(kind=dbl) :: alpi
!
integer, parameter :: nbins1 =  9
integer, parameter :: nbins2 = 10
integer, parameter :: nbins3 =  6
integer, parameter :: nbins4 =  24
integer, parameter :: nbins5 =  6
integer, parameter :: nbins6 =  1
integer, parameter :: nbinstot = 56
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)
real(kind=dbl) :: error(nbinstot),error4(nbinstot),errorbis(nbinstot)
!
real(kind=dbl), parameter :: metric(0:3) = (/1.0d0,-1.0d0,-1.0d0,-1.0d0/)
real(kind=dbl) :: kf(maxparticles,0:3)  ! final state momenta
integer :: nparticles
real(kind=dbl) :: temp
integer :: i,j
!
real(kind=dbl), parameter :: ycut(6) = &
              (/0.005d0,0.01d0,0.02d0,0.05d0,0.1d0,0.2d0/)
real(kind=dbl), parameter :: y3cut(7) = &               ! 3.16228 is the square root of 10
              (/0.000316228d0,0.001d0,0.00316228d0,0.01d0,0.0316228d0,0.1d0,0.316228d0/)
real(kind=dbl) :: y43,y32
real(kind=dbl) :: jetmass(3), mass0,mass1
real(kind=dbl), parameter :: jetmassbinsize = 1.0d0
!
type(eventdata) :: theevent  ! event data for Pythia
! Module pythiawulf has predefined the following
!   type(eventdata), dimension(eventrecordsize) :: event
!   integer :: eventindex 
integer :: evt      ! stored event number
integer :: evt0     ! stored event number at start of new point
real(kind=dbl) :: probability
real(kind=dbl) :: oldweight,newweight,rts
logical :: acceptevent
!
real(kind=dbl) :: random
logical :: badpointq,xtrapointq
real(kind=dbl) :: calsvalue,calsthrust,calstmoments
integer :: n,mu
real(kind=dbl) :: thrustdist,thrustval,thrust
!
save sum,sumbis,sum2,sum2bis,sum3,sum4,integral,integralbis, &
     value,badpointq,xtrapointq,evt,evt0
!
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
  IF (pythia) THEN 
    evt = 0      ! stored event number
    evt0 = 0     ! stored event number at start of new point
  END IF
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.
!---------------------------------------------------------------
!
ELSE IF (whattodo .EQ. 'newresult ') THEN
!
! First, it is convenient to have a list of the momenta of the final state
! momenta, kf(n,mu). We also extract the number of final state particles.
!
  j = 0
  rts = 0.0d0
  DO n = 1,theshower%length
  IF (theshower%ptn(n)%childless) THEN
    j = j + 1
    temp = 0.0d0
    DO mu = 1,3
      kf(j,mu) = theshower%ptn(n)%momentum(mu)
      temp = temp +  kf(j,mu)**2
    END DO
    temp = sqrt(temp)
    kf(j,0) = temp
    rts = rts + temp
  END IF
  END DO
  nparticles = j
!
! Calculate the functions calS for each type of
! observable using functions calstype(nparticles,kf,index).
!
  thrustval = thrust(kf,nparticles)
  n = 0
  DO j = 1,nbins1
     n = n+1
     calsvalue = calsthrust(thrustval,j)
     value(n) = value(n) + weight*calsvalue
  END DO
  DO j = 1,nbins2
     n = n+1
     calsvalue = calstmoments(thrustval,j)
     value(n) = value(n) + weight*calsvalue
  END DO
!
! Jets with Durham (KT) algorithm.
!
  call findjets(kf,nparticles,'durham',y43,y32,jetmass)
  ! We divide by the ratio of sigma_T to sigma_0 at order alpha_s^1.
  calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts))
  DO j = 1,nbins3
    n = n+1
    IF ((y43.LT.ycut(j)).AND.(ycut(j).LE.y32)) THEN
      ! We have three jets with ycut = ycut(j).
      value(n) = value(n) + weight*calsvalue
    END IF
  END DO
!
! Jet mass distribution for three jets. Note that kf vectors given to findjets 
! were of the dimensionless variety, so we have to translate to get distribution ! in bins with width given in GeV. 
!
  IF ((y43.LT.ycut(4)).AND.(ycut(4).LE.y32)) THEN
    ! We have three jets with ycut = ycut(4).
    ! We divide by the ratio of sigma_T to sigma_0 at order alpha_s^1
    ! and by the bin size.
    calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts)) &
                /jetmassbinsize
    DO j = 1,nbins4
      n = n+1
      mass0 = (j-1)*jetmassbinsize*rts/externalrts
      IF (j.EQ.1) THEN
         mass0 = -1.0d0 ! Just to make sure we catch zero mass jets.
      END IF
      mass1 = j*jetmassbinsize*rts/externalrts
      DO i = 1,3
        IF ((jetmass(i).GE.mass0).AND.(jetmass(i).LT.mass1)) THEN
          value(n) = value(n) + weight*calsvalue
        END IF
      END DO
    END DO
  ELSE
    n = n + nbins4
  END IF
!
! Distribution of y32 in bins.
!
  calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts))
  DO j = 1,nbins5
    n = n+1
    IF ((y3cut(j).LT.y32).AND.(y32.LT.y3cut(j+1))) THEN
      ! y32 is in the jth bin.
      value(n) = value(n) + weight*calsvalue
    END IF
  END DO
!
! Jets with Jade algorithm.
!
!  call findjets(kf,nparticles,'jade  ',y43,y32,jetmass)
!
! For now, we don't do Jade jets.
!
!
! Fraction with thrust < thrustcut.
!
n = n + nbins6  ! nbins6 is 1
IF (thrustval.LT.thrustcut) THEN
  calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts))
  value(n) = value(n) + weight*calsvalue
END IF
!
! In addition, we may want to produce an event record for Pythia.
!
  IF (pythia) THEN
!
!  Note that in beowulf's calculation one sums over events and divides by nreno.
!
!   <S> = (1/nreno) sum_j w_j S_j
!
!  We now accept events with probability p_j, so that
!
!   <S> = (1/nreno) sum_j (w_j/p_j) \theta(acceptevent_j) S_j
!
!  Thus we record a new weight equal to oldweight/probability. We are
!  still dividing by nreno. So if standard weights are defined by
!
!   <S> = (1/naccepted) sum_j \theta(acceptevent_j) * standardweight_j S_j
!
!  then the standard weights are the beowulf weights times naccepted/nrenno.
!
    oldweight = weight
    ! The formula for probability could be changed.
    probability = abs(1.0d4*abs(oldweight)*(1.0d0 - thrustval)**2)
    IF (probability.GT.1.0d0 ) THEN
      acceptevent = .true.
      newweight = oldweight  ! We set probability to 1.
    ELSE IF ( probability.GT.random(1) ) THEN
      acceptevent = .true.
      newweight = oldweight/probability
    ELSE
      acceptevent = .false.
    END IF
    IF (thrustval.GT.thrustcut) THEN
      acceptevent = .false.
    END IF
!
    IF (acceptevent) THEN
      evt = evt + 1                      ! new stored event number
      IF (evt.GT.eventrecordsize) THEN
        write(*,*) 'Need to enlarge eventrecordsize ',evt
        STOP
      END IF
!
      IF (mode.EQ.'born    ') THEN
      ! We may want born level, which did not define color connections.
        DO n = 1,3
          IF (theshower%ptn(n)%type.EQ.'quark') THEN
            theshower%ptn(n)%stringquark = 501
            theshower%ptn(n)%stringqbar  = 0
          ELSE IF (theshower%ptn(n)%type.EQ.'qbar ') THEN
            theshower%ptn(n)%stringquark = 0
            theshower%ptn(n)%stringqbar  = 502
          ELSE                             ! it is the gluon
            theshower%ptn(n)%stringquark = 502
            theshower%ptn(n)%stringqbar  = 501
          END IF
        END DO
      END IF
!
      call prune(theshower,rts)  ! Remove small virtuality branches.
      call translate(theshower,newweight,theevent)
      ! Produces theevent, with information needed for Pythia.
      event(evt) = theevent
    END IF ! (acceptevent)
!
  END IF ! (pythia)
!---------------------------------------------------------------
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
  IF (pythia) THEN
    IF (badpointq.OR.xtrapointq) THEN
      ! Bad point. Restore event counter to value before this point.
      ! This effectively erases the information recorded for this point.
      evt = evt0 
    ELSE
      evt0 = evt    ! Good Point. Update evt0 to current value of event counter.
    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
  IF (pythia) THEN
    eventindex = evt          ! eventindex of last event, for upevnt
    evt = 0
    evt0 = 0
  END IF
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 j = 1,nbins1
    n = n+1
    write(nout,002)(0.68d0 + j*0.03d0), &
         thrustdist(0.68d0 + j*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 j = 1,nbins2
    n = n+1
    write(nout,007)1.0d0+j*0.5d0
    write(nout,003)sum(n),error(n)
    write(nout,004)sumbis(n),errorbis(n)
    write(nout,005)error4(n)
  END DO
  write(nout,008)
  DO j = 1,nbins3
    n = n+1
    write(nout,009)ycut(j)
    write(nout,003)sum(n),error(n)
    write(nout,004)sumbis(n),errorbis(n)
    write(nout,005)error4(n)
  END DO
  write(nout,010)
  write(nout,009)ycut(4)
  write(nout,*)' '
  DO j = 1,nbins4
    n = n+1
    mass0 = (j-1)*jetmassbinsize
    mass1 = j*jetmassbinsize
    write(nout,011)mass0,mass1
    write(nout,003)sum(n),error(n)
    write(nout,004)sumbis(n),errorbis(n)
    write(nout,005)error4(n)
  END DO
  write(nout,012)
  DO j = 1,nbins5
    n = n+1
    write(nout,013)y3cut(j),y3cut(j+1)
    write(nout,003)sum(n),error(n)
    write(nout,004)sumbis(n),errorbis(n)
    write(nout,005)error4(n)
  END DO
  write(nout,14)
  write(nout,15)thrustcut
  n=n+1
  write(nout,003)sum(n),error(n)
  write(nout,004)sumbis(n),errorbis(n)
  write(nout,005)error4(n)
!-------------------------------------------------------
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)
008   format(/,  &
          'Next, three jet cross section/sigma_T, Durham algorithm',/)
009   format('ycut =',f10.4)
010   format('Finally, mass distributions of the three jets.')
011   format(f10.3,' < jet mass <',f10.3)
012   format(/,  &
          'Next, cross section/sigma_T, in bins of y32 (Durham algorithm)',/)
013   format(f10.5,' < y32 <',f10.5)
014   format(/,  &
          'Next, (cross section(thrust < thrustcut))/sigma_T',/)
015   format('with thrustcut = ',f10.3)
!
END subroutine hrothgar
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine hrothgarII(nreno,whattodo)
!
use beowulf_parameters
use beowulf_structures
implicit none
!
! In:
integer :: nreno
character(len=10) :: whattodo
!
! Accumulates the results from Pythia.
! 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.
! '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.
! '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.
!
! This subroutine is pretty much identical to subroutine hrothgar. However,
! it does not prepare events for Pythia since it is getting the events from
! Pythia. In addition, it does not accumulate 'cutoff correction' results
! in integralbis(n), sumbis(n), sum2bis(n). Events corresponding to reno points
! that are beyond the cutoffs are simply not sent to Pythia.  Finally, there
! are one fewer levels of results accumulation. In instead of 
! value --> integral --> sum in hrothgar, we have value --> sum in hrothgarII.
!
!  3 January 2004
!  4 February 2005
!
real(kind=dbl) :: badnesslimit,cancellimit,thrustcut
common /limits/ badnesslimit,cancellimit,thrustcut
integer :: groupsize(maxgraphs,maxmaps)
integer :: groupsizegraph(maxgraphs)
integer :: groupsizetotal
common /montecarlo/groupsize,groupsizegraph,groupsizetotal
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts
!
real(kind=dbl) :: muoverrts
common /renormalize/ muoverrts
real(kind=dbl) :: nc,nf,cf
common /colorfactors/ nc,nf,cf
real(kind=dbl) :: alpi
!
integer, parameter :: nbins1 =  9
integer, parameter :: nbins2 = 10
integer, parameter :: nbins3 =  6
integer, parameter :: nbins4 = 24
integer, parameter :: nbins5 =  6
integer, parameter :: nbins6 =  1
integer, parameter :: nbinstot = 56
real(kind=dbl) :: sum(nbinstot)
real(kind=dbl) :: sum2(nbinstot)
real(kind=dbl) :: sum3(nbinstot),sum4(nbinstot)
real(kind=dbl) :: value(nbinstot)
real(kind=dbl) :: error(nbinstot),error4(nbinstot)
!
real(kind=dbl), parameter :: metric(0:3) = (/1.0d0,-1.0d0,-1.0d0,-1.0d0/)
real(kind=dbl) :: kf(maxparticles,0:3)  ! final state momenta
real(kind=dbl) :: rts
integer :: nparticles
integer :: i,j
!
real(kind=dbl), parameter :: ycut(6) = &
              (/0.005d0,0.01d0,0.02d0,0.05d0,0.1d0,0.2d0/)
real(kind=dbl), parameter :: y3cut(7) = &    ! 3.16228 is the square root of 10
     (/0.000316228d0,0.001d0,0.00316228d0,0.01d0,0.0316228d0,0.1d0,0.316228d0/)
real(kind=dbl) :: y43,y32
real(kind=dbl) :: jetmass(3), mass0,mass1
real(kind=dbl), parameter :: jetmassbinsize = 1.0d0
!
real(kind=dbl) :: weight
real(kind=dbl) :: calsvalue,calsthrust,calstmoments
integer :: n
real(kind=dbl) :: thrustdist,thrustval,thrust
!
!---------Pythia user process event common block----------------
!
integer, parameter :: MAXNUP = 500
integer :: NUP, IDPRUP, IDUP(MAXNUP), ISTUP(MAXNUP),  &
           MOTHUP(2,MAXNUP), ICOLUP(2,MAXNUP)
real(kind=dbl) :: XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP(5,MAXNUP), &
                  VTIMUP(MAXNUP), SPINUP(MAXNUP)
common /HEPEUP/ NUP, IDPRUP, XWGTUP, SCALUP, AQEDUP, AQCDUP, IDUP,   &
                ISTUP, MOTHUP, ICOLUP, PUP, VTIMUP, SPINUP
save /HEPEUP/
!
save sum,sum2,sum3,sum4,value
!
rts = externalrts
IF (whattodo .EQ. 'startcalc ') THEN
  DO n = 1,nbinstot
    sum(n)     = 0.0d0
    sum2(n)    = 0.0d0
    sum3(n)    = 0.0d0
    sum4(n)    = 0.0d0
  END DO
ELSE IF (whattodo .EQ. 'startgroup') THEN
  DO n = 1,nbinstot
    value(n) = 0.0d0
  END DO
!---------------------------------------------------------------
!
ELSE IF (whattodo .EQ. 'newresult ') THEN
!
! First, it is convenient to have a list of the momenta of the final state
! momenta, kf(n,mu). We also extract the number of final state particles.
!
  call makekf(nparticles,kf)  ! Gets info from Pythia common block/PYJETS/.
  weight = XWGTUP             ! From Pythia common block /HEPEUP/.
!
! Calculate the functions calS for each type of
! observable using functions calstype(nparticles,kf,index).
!
  thrustval = thrust(kf,nparticles)
  n = 0
  DO j = 1,nbins1
     n = n+1
     calsvalue = calsthrust(thrustval,j)
     value(n) = value(n) + weight*calsvalue
  END DO
  DO j = 1,nbins2
     n = n+1
     calsvalue = calstmoments(thrustval,j)
     value(n) = value(n) + weight*calsvalue
  END DO
!
! Jets with Durham (KT) algorithm.
!
  call findjets(kf,nparticles,'durham',y43,y32,jetmass)
  ! We divide by the ratio of sigma_T to sigma_0 at order alpha_s^1.
  calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts))
  DO j = 1,nbins3
    n = n+1
    IF ((y43.LT.ycut(j)).AND.(ycut(j).LE.y32)) THEN
      ! We have three jets with ycut = ycut(j).
      value(n) = value(n) + weight*calsvalue
    END IF
  END DO
!
! Jet mass distribution. This code is copied from hrothgar. But here
! rts/externalrts = 1.
!
  IF ((y43.LT.ycut(4)).AND.(ycut(4).LE.y32)) THEN
    ! We have three jets with ycut = ycut(4).
    ! We have three jets with ycut = ycut(4).
    ! We divide by the ratio of sigma_T to sigma_0 at order alpha_s^1
    ! and by the bin size.
    calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts)) &
                /jetmassbinsize
    DO j = 1,nbins4
      n = n+1
      mass0 = (j-1)*jetmassbinsize*rts/externalrts
      IF (j.EQ.1) THEN
         mass0 = -1.0d0 ! Just to make sure we catch zero mass jets.
      END IF
      mass1 = j*jetmassbinsize*rts/externalrts
      DO i = 1,3
        IF ((jetmass(i).GE.mass0).AND.(jetmass(i).LT.mass1)) THEN
          value(n) = value(n) + weight*calsvalue
        END IF
      END DO
    END DO
  ELSE
    n = n + nbins4
  END IF
!
! Distribution of y32 in bins.
!
  calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts))
  DO j = 1,nbins5
    n = n+1
    IF ((y3cut(j).LT.y32).AND.(y32.LT.y3cut(j+1))) THEN
      ! y32 is in the jth bin.
      value(n) = value(n) + weight*calsvalue
    END IF
  END DO
! Jets with Jade algorithm.
!
!  call findjets(kf,nparticles,'jade  ',y43,y32,jetmass)
!
! For now, we don't do Jade jets.
!
! Fraction with thrust_0 < thrustcut, so counted by Pythia.
!
n = n + nbins6  ! nbins6 is 1
calsvalue = 1.0d0/(1.0d0 + 0.75d0*cf*alpi(muoverrts*externalrts))
value(n) = value(n) + weight*calsvalue
!
!---------------------------------------------------------------
!
ELSE IF (whattodo .EQ. 'groupdone ') THEN
  DO n = 1,nbinstot
    sum(n)         = sum(n)  + value(n)
    sum2(n)        = sum2(n) + value(n)**2
    sum3(n)        = sum3(n) + value(n)**3
    sum4(n)        = sum4(n) + value(n)**4
  END DO
ELSE IF (whattodo .EQ. 'calcdone  ') THEN
  write(nout,001)
  DO n = 1,nbinstot
    sum(n)       = sum(n)/nreno
    sum2(n)      = sum2(n)/nreno
    error(n)     = sqrt((sum2(n) - sum(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 j = 1,nbins1
    n = n+1
    write(nout,002)(0.68d0 + j*0.03d0), &
         thrustdist(0.68d0 + j*0.03d0)
    write(nout,003)sum(n),error(n)
    write(nout,005)error4(n)
  END DO
  write(nout,006)
  DO j = 1,nbins2
    n = n+1
    write(nout,007)1.0d0+j*0.5d0
    write(nout,003)sum(n),error(n)
    write(nout,005)error4(n)
  END DO
  write(nout,008)
  DO j = 1,nbins3
    n = n+1
    write(nout,009)ycut(j)
    write(nout,003)sum(n),error(n)
    write(nout,005)error4(n)
  END DO
  write(nout,010)
  write(nout,009)ycut(4)
  write(nout,*)' '
  DO j = 1,nbins4
    n = n+1
    mass0 = (j-1)*jetmassbinsize
    mass1 = j*jetmassbinsize
    write(nout,011)mass0,mass1
    write(nout,003)sum(n),error(n)
     write(nout,005)error4(n)
  END DO
  write(nout,012)
  DO j = 1,nbins5
    n = n+1
    write(nout,013)y3cut(j),y3cut(j+1)
    write(nout,003)sum(n),error(n)
    write(nout,005)error4(n)
  END DO
  write(nout,14)
  write(nout,15)thrustcut
  n=n+1
  write(nout,003)sum(n),error(n)
  write(nout,005)error4(n)
!-------------------------------------------------------
END IF
RETURN
!
! FORMAT statements for the output.
!
001   format(/, '       ---   Hrothgar II 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)
008   format(/,  &
          'Next, three jet cross section/sigma_T, Durham algorithm',/)
009   format('ycut =',f10.4)
010   format('Finally, mass distributions of the three jets.')
011   format(f10.3,' < jet mass <',f10.3)
012   format(/,  &
          'Next, cross section/sigma_T, in bins of y32 (Durham algorithm)',/)
013   format(f10.5,' < y32 <',f10.5)
014   format(/,  &
          'Next, (cross section from Pythia)/sigma_T',/)
015   format('with thrustcut = ',f10.3)
!
END subroutine hrothgarII
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
function calsthrust(t,n)
!
use beowulf_parameters
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_parameters
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
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)  
!
use beowulf_parameters
implicit none
! In:
real(kind=dbl) :: kf(maxparticles,0:3)
integer :: nparticles
! Out:
real(kind=dbl) :: thrust
!
! kf gives momenta for nparticles final state particles.
! Then thrust is the thrust of this final state.
!
! 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,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_parameters
implicit none
! In :
real(kind=dbl) :: kf(maxparticles,0:3)
integer :: nparticles
real(kind=dbl), dimension(3) :: vin
! Out:
real(kind=dbl), dimension(3) :: 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_parameters
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=8) :: mode
common /programmode/ mode
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,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.'showerI ').OR.(mode.EQ.'showerII')) 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
!
subroutine findjets(kf,nparticles,algorithm,y43,y32,jetmass)
!
use beowulf_parameters
use beowulf_structures
implicit none
! In:
real(kind=dbl) :: kf(maxparticles,0:3)
integer :: nparticles
character(len=6) :: algorithm
! Out:
integer :: njets
real(kind=dbl) :: jetmass(3)
real(kind=dbl) :: y43,y32
!
! Groups the particles into jets according to the algorithm specified. 
! Reports the masses of the jets when there are three jets
! as well as the values of ycut at which the number of jets changes from
! 4 to 3 and from 3 to 2. The algorithm codes are
! 'durham' and 'jade  ' and 'pythia'.
!
! 19 March 2004
! 24 May 2005
!
type jetstructure
  logical :: alive
  real(kind=dbl) :: k(0:3)
  real(kind=dbl) :: absk
  integer :: nparticles
  integer :: label(maxparticles)
end type jetstructure
type(jetstructure), dimension(maxparticles) :: jet
!
integer n,mu,i,j,istar,jstar
real(kind=dbl) :: ksq
real(kind=dbl) :: rts       ! Energy of event
real(kind=dbl), dimension(maxparticles,maxparticles) :: y
real(kind=dbl) :: ymin
real(kind=dbl) :: dotij,cosij,temp
!
! Calculate rts and the initial jet list.
!
rts = 0.0d0
njets = nparticles
DO n = 1,nparticles
  jet(n)%alive = .true.
  jet(n)%nparticles = 1
  jet(n)%label(1) = n
  ksq = 0.0d0
  DO mu = 1,3
    temp = kf(n,mu)
    ksq = ksq + temp**2
    jet(n)%k(mu) = temp
  END DO
  IF (ksq.LT.1d-10) THEN
    jet(n)%alive = .false.  ! Pythia seems to give some zero energy photons
    njets = njets - 1
  ELSE
    jet(n)%k(0) = kf(n,0)
    rts = rts + jet(n)%k(0)
    jet(n)%absk = sqrt(ksq)
    IF (jet(n)%k(0).LT.0.99999999d0*jet(n)%absk) THEN
      write(nout,*)'Tachyon discovered in findjets',jet(n)%k(0),jet(n)%absk
      STOP
    END IF
  END IF
END DO
!
IF (njets.LE.2) THEN
  ! This can happen if there were three particles and one was very soft.
  y43 = 0.0d0
  y32 = 0.0d0
  jetmass(1) = 0.0d0
  jetmass(2) = 0.0d0
  jetmass(3) = 0.0d0
  RETURN
END IF
!
! Calculate the initial y(i,j)
!
DO i = 1,nparticles
DO j = i+1,nparticles
  IF ((jet(i)%alive).AND.(jet(j)%alive)) THEN
    dotij = 0.0d0
    DO mu = 1,3
      dotij = dotij + jet(i)%k(mu) * jet(j)%k(mu)
    END DO
    cosij = dotij/jet(i)%absk/jet(j)%absk
    IF (algorithm.EQ.'durham') THEN
      temp = 1.0d0 - cosij
      temp = temp * 2.0d0 * min(jet(i)%k(0),jet(j)%k(0))**2
      y(i,j) = temp / rts**2
    ELSE IF (algorithm.EQ.'jade  ') THEN
      temp = 1.0d0 - cosij
      temp = temp * 2.0d0 * jet(i)%k(0) * jet(j)%k(0)
      y(i,j) = temp / rts**2
    ELSE IF (algorithm.EQ.'pythia') THEN
      temp = 1.0d0 - cosij
      temp = temp * 2.0d0 * &
       ( jet(i)%k(0) * jet(j)%k(0)/(jet(i)%k(0) + jet(j)%k(0)) )**2
      y(i,j) = temp / rts**2
    ELSE
     write(nout,*)'Not prepared for algorithm ',algorithm,' in findjets.'
     STOP
    END IF
  END IF
END DO
END DO
!
y43 = 777.77d0        ! This will be changed.
y32 = 777.77d0        ! This will be changed.
jetmass(1) = 777.77d0 ! This will be changed.
jetmass(2) = 777.77d0 ! This will be changed.
jetmass(3) = 777.77d0 ! This will be changed.
ymin = 0.0d0 ! We will set y43 to this if we start with just three jets.
!
! Main loop. Stops when we get to two jets.
!
DO
!
IF (njets.EQ.3) THEN
  y43 = ymin  ! This is the ycut for 4 to 3 jets, = 0 if start with three jets.
  ! Find the jet masses.
  n = 0
  DO j = 1,nparticles
    IF (jet(j)%alive) THEN
      n = n + 1
      IF (n.GT.3) THEN
        write(nout,*)'Oops, more than three jets for njets = 3.'
        STOP
      END IF
      ksq = jet(j)%k(0)**2 - jet(j)%absk**2
      IF (ksq.GT.0.0d0) THEN
        jetmass(n) = sqrt(ksq)
      ELSE
        jetmass(n) = 0.0d0       ! This allows for roundoff errors.
      END IF
    END IF
  END DO
ELSE IF (njets.LT.2) THEN ! Test for njets.EQ.2 is below.
  write(nout,*)'Jets lost in findjets',njets,nparticles
  STOP
END IF
!
! Find smallest y(i,j), ymin = y(i*,j*).
!
ymin = 1.0d30
DO i = 1,nparticles
DO j = i+1,nparticles
  IF ((jet(i)%alive).AND.(jet(j)%alive)) THEN
    IF ( y(i,j).LT.ymin ) THEN
      ymin = y(i,j)
      istar = i
      jstar = j
    END IF
  END IF
END DO
END DO
!
! Combine jets istar and jstar. Record by overwriting jet istar.
!
ksq = 0.0d0
DO mu = 1,3
  temp = jet(istar)%k(mu) + jet(jstar)%k(mu)
  ksq = ksq + temp**2
  jet(istar)%k(mu) = temp
END DO
IF (ksq.LE.0.0d0) THEN
 write(nout,*)'New jet ksq.LE.0 in findjets',ksq
 STOP
END IF
jet(istar)%k(0) = jet(istar)%k(0) + jet(jstar)%k(0)
jet(istar)%absk = sqrt(ksq)
IF (jet(istar)%k(0).LT.0.99999999d0*jet(istar)%absk) THEN
  write(nout,*)'Tachyon discovered in findjets',jet(istar)%k(0),jet(istar)%absk
  STOP
END IF
DO n = 1,jet(jstar)%nparticles
  jet(istar)%label(jet(istar)%nparticles+n) = jet(jstar)%label(n)
END DO
jet(istar)%nparticles = jet(istar)%nparticles + jet(jstar)%nparticles
jet(jstar)%alive = .false.  ! This effectively erases jet jstar.
njets = njets - 1
!
IF (njets.EQ.2) THEN
  y32 = ymin   ! This is the ycut for 3 to 2 jets.
  RETURN       ! We are done.
END IF
!
! Revise list of y(i,j). (We don't do this if we have two jets.)
!
DO j = 1,nparticles
  IF (jet(j)%alive.AND.(j.NE.istar)) THEN
    dotij = 0.0d0
    DO mu = 1,3
      dotij = dotij + jet(istar)%k(mu) * jet(j)%k(mu)
    END DO
    cosij = dotij/jet(istar)%absk/jet(j)%absk
    IF (algorithm.EQ.'durham') THEN
      temp = 1.0d0 - cosij
      temp = temp * 2.0d0 * min(jet(istar)%k(0),jet(j)%k(0))**2
     ELSE IF (algorithm.EQ.'jade  ') THEN
      temp = 1.0d0 - cosij
      temp = temp * 2.0d0 * jet(istar)%k(0) * jet(j)%k(0)
    ELSE IF (algorithm.EQ.'pythia') THEN
      temp = 1.0d0 - cosij
      temp = temp * 2.0d0 * &
       ( jet(istar)%k(0) * jet(j)%k(0)/(jet(istar)%k(0) + jet(j)%k(0)) )**2
    ELSE
     write(nout,*)'Not prepared for algorithm ',algorithm,' in findjets.'
     STOP
    END IF
    IF (istar.LT.j) THEN
      y(istar,j) = temp / rts**2
    ELSE
      y(j,istar) = temp / rts**2
    END IF
  END IF
END DO
!
END DO ! Main loop
!
END subroutine findjets
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
function kn0(t)
!
use beowulf_parameters
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_parameters
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_parameters
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
!   3 November 2003
!
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
!
subroutine prune(theshower,rts)
!
use beowulf_parameters
use beowulf_structures
implicit none
! In:
real(kind=dbl) :: rts
! In and Out:
type(showerlist) :: theshower
!
! Prunes branches with small virtuality from shower.
!
!  9 December 2003
!  1 March 2004
!
real(kind=dbl) :: showerendratio
common /showerend/ showerendratio
type(showerlist) :: shower
!
real(kind=dbl) :: virtuality,minvirtuality,newrts
integer :: n,nn,mu,c1,c2,s1,s2
logical :: change
real(kind=dbl) :: pnsq,p1sq,p2sq,e1,e2
real(kind=dbl) :: pgluonsq,pgluon(3)
logical, dimension(maxparticles)  :: erased
logical :: nothingchanged
integer :: newname(0:maxparticles),oldname(maxparticles)
!
! Find splittings at the end of the tree and prune them if the virtuality
! is too low or if gluon is too soft.
!
nothingchanged = .true.
shower = theshower
DO n = 1,shower%length
  erased(n) = .false.
END DO
newrts = rts
DO  ! Do as long as there is change.
  change = .false.
  minvirtuality = showerendratio*newrts**2
  DO n = 1,shower%length
  IF (.NOT.shower%ptn(n)%childless) THEN
    c1 = shower%ptn(n)%child1
    c2 = shower%ptn(n)%child2
    IF (shower%ptn(c1)%childless.AND.shower%ptn(c2)%childless) THEN
      ! parton n is parent of final state partons
      pnsq = 0.0d0
      p1sq = 0.0d0
      p2sq = 0.0d0
      DO mu = 1,3
        pnsq = pnsq + shower%ptn(n)%momentum(mu)**2
        p1sq = p1sq + shower%ptn(c1)%momentum(mu)**2
        p2sq = p2sq + shower%ptn(c2)%momentum(mu)**2
      END DO
      e1 = sqrt(p1sq)
      e2 = sqrt(p2sq)
      virtuality =   (e1 + e2)**2 - pnsq
      IF (virtuality.LT.minvirtuality) THEN  ! Prune these branches.
        erased(c1) = .true.
        erased(c2) = .true.
        shower%ptn(n)%childless = .true.
        shower%ptn(n)%done = .true.
        shower%ptn(n)%child1 = -1
        shower%ptn(n)%child2 = -1
        newrts = newrts - e1 - e2 + sqrt(pnsq) ! Adjust rts.
        change = .true.
        nothingchanged = .false.
      END IF
    END IF  ! shower%ptn(c1)%childless.AND.shower%ptn(c2)%childless
  END IF   ! .NOT.shower%ptn(n)%childless
  END DO   ! n = 1,shower%length
  IF (.NOT.change) exit
END DO    ! Do as long as there is change.
!
! If any gluon with parent 0 is too soft, erase it. For sqrt(showerendratio)
! * sqrts = 1 GeV, we take 'too soft' to be E = 100 MeV. We also give both
! strings connected to this gluon the same name. Half of the momentum of the 
! erased gluon is assigned to the partons connected to its quark string and
! half is assigned to the partons connected to its qbar string.
!
DO n = 1,shower%length
IF ( (shower%ptn(n)%childless)       & 
     .AND.(shower%ptn(n)%type.EQ.'gluon') &
     .AND.(shower%ptn(n)%parent.EQ.0)      ) THEN
  pgluonsq = 0.0d0
  DO mu = 1,3
    pgluon(mu) = shower%ptn(n)%momentum(mu)
    pgluonsq = pgluonsq + pgluon(mu)**2
  END DO
  IF (pgluonsq.LT.showerendratio*newrts**2/1.0d2) THEN
    erased(n) = .true.
    nothingchanged = .false.
    s1 = shower%ptn(n)%stringquark
    s2 = shower%ptn(n)%stringqbar
    DO nn = 1,shower%length
      IF (shower%ptn(nn)%stringquark.EQ.s2) THEN
        shower%ptn(nn)%stringquark = s1
        DO mu = 1,3
          shower%ptn(nn)%momentum(mu) = &
             shower%ptn(nn)%momentum(mu) + 0.5d0*pgluon(mu)
        END DO
      END IF
      IF (shower%ptn(nn)%stringqbar.EQ.s1) THEN
        DO mu = 1,3
          shower%ptn(nn)%momentum(mu) = &
             shower%ptn(nn)%momentum(mu) + 0.5d0*pgluon(mu)
        END DO
      END IF
    END DO
  END IF  ! p1sq.LT.showerendratio*newrts**2/1.0d2
END IF  ! shower%ptn(n)%childless et cetera
END DO  ! n = 1,shower%length
!
IF (nothingchanged) THEN
  RETURN
END IF
!
! Something has changed. Now we need a new shower list with the 
! specified lines erased.
!
nn =0
DO n = 1,shower%length
  IF(.NOT.erased(n)) THEN
    nn = nn + 1
    newname(n) = nn
    oldname(nn) = n
  END IF
END DO
newname(0) = 0           ! This may be used for parent names.
theshower%length = nn
!
DO nn = 1,theshower%length
  n = oldname(nn)
  theshower%ptn(nn) = shower%ptn(n)
  theshower%ptn(nn)%self = nn
  theshower%ptn(nn)%parent = newname(shower%ptn(n)%parent)
  ! Ancestor names cannot have been renamed.
  IF (.NOT.theshower%ptn(nn)%childless) THEN
    theshower%ptn(nn)%child1 = newname(shower%ptn(n)%child1)
    theshower%ptn(nn)%child2 = newname(shower%ptn(n)%child2)
  END IF
END DO ! nn = 1,theshower%length
!
RETURN
END subroutine prune
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine translate(theshower,weight,theevent)
use beowulf_parameters
use beowulf_structures
use pythiawulf   ! For passing list of events, event, to upevnt
implicit none
! In:
type(showerlist) :: theshower
real(kind=dbl)   :: weight
! Out:
type(eventdata) :: theevent
!
! Produces event record with information we will need to send to Pythia,
! including dimensionful 4-momenta for partons instead of dimensionless
! 3-momenta. The final state partons are given Pythia default massses.
!
! 7  December 2003
! 1  March 2004
! 11 June 2005
!
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts
!
real(kind=dbl), parameter :: metric(0:3) = (/1.0d0,-1.0d0,-1.0d0,-1.0d0/)
real(kind=dbl) :: rts
real(kind=dbl) :: kf(maxparticles,0:3)    ! final state momenta
real(kind=dbl) :: kfsq(maxparticles)      ! square of 3-momenta
real(kind=dbl) :: kfmasssq(maxparticles)  ! final state masses
integer :: nparticles
real(kind=dbl) :: temp,delta
integer :: j,n,nn,n1,n2,mu,mother,grandmother,numberofquarks
type(partonmod), dimension(maxparticles) :: theparton  ! For Pythia.
real(kind=dbl), parameter :: partonmass(0:6) =  &       ! for Pythia
      (/0.0d0, 0.33d0, 0.33d0, 0.5d0, 1.5d0, 4.8d0, 175.0d0/)
!        g       u       d       s      c      b      t     beowulf convention
real(kind=dbl) :: factor,approxrts
integer :: count
logical :: active(maxparticles),new(maxparticles,maxparticles)
real(kind=dbl) :: y12,ymin,yvalues(maxparticles,maxparticles)
real(kind=dbl) :: dot11,dot22,dot12,cos12,dot12star
real(kind=dbl) :: dot12values(maxparticles,maxparticles)
integer :: n1star,n2star
logical :: ok,found
real(kind=dbl) :: ymax
integer :: partner(maxparticles)  ! partner of each quark or qbar
integer :: tempstring
logical :: partnerfound,done,known(maxparticles)
!
! First, it is convenient to have a list of the momenta of the final state
! momenta, kf(n,mu). We also extract the number of final state particles.
!
j = 0
rts = 0
DO n = 1,theshower%length
IF (theshower%ptn(n)%childless) THEN
  j = j + 1
  temp = 0.0d0
  DO mu = 1,3
    kf(j,mu) = theshower%ptn(n)%momentum(mu)
    temp = temp +  kf(j,mu)**2
  END DO
  kfsq(j) = temp
  temp = sqrt(temp)
  kf(j,0) = temp
  rts = rts + temp
END IF
END DO
nparticles = j
!
! Collect mass information and get dimensionful kf approximately.
!
j = 0
DO n = 1,theshower%length
  IF (theshower%ptn(n)%childless) THEN
    j = j + 1
    kfmasssq(j) = partonmass(abs(theshower%ptn(n)%flavor))**2
    DO mu = 0,3
      kf(j,mu) = kf(j,mu) * externalrts/rts ! Approximately right, in GeV.
    END DO
    kfsq(j) = kfsq(j) * (externalrts/rts)**2
  END IF
END DO
!
! Revise kf(j,mu) to account for masses.
!
count = 0
DO
  count = count + 1
  IF (count.GT.3) THEN
    write(nout,*)'Too many loops to get masses right', count
  END IF
  IF (count.GT.10) STOP
  factor = 0.0d0
  approxrts = 0.0d0
  DO j = 1,nparticles
     kf(j,0) = sqrt(kfsq(j) + kfmasssq(j))
     approxrts = approxrts + kf(j,0)
     factor = factor + kfsq(j)/kf(j,0)
  END DO
  delta = (externalrts - approxrts)/factor
  IF (abs(delta).LT.1.0d-6) EXIT
  DO j = 1,nparticles
    DO mu = 1,3
      kf(j,mu) = kf(j,mu)*(1.0d0 + delta)
    END DO
    kfsq(j) = kfsq(j)*(1.0d0 + delta)**2
  END DO
END DO
!
theevent%length     = 2*nparticles - 2
theevent%weight     = weight
theevent%nparticles = nparticles
!
! Insert momentum information for final state partons.
! Allowing for mothers numbered 1,...,nparticles - 2,
! we insert momenta for final state particles numbered from 
! (nparticles - 2) + 1  to (nparticles - 2) + nparticles.
!
DO n = 1,theevent%length
  active(n) = .false.     ! The active jets in the kT algorithm; so far none.
END DO
numberofquarks = 0
j = 0
nn = nparticles - 2
DO n = 1,theshower%length
  IF (theshower%ptn(n)%childless) THEN
    j = j + 1
    nn = nn + 1
    DO mu = 0,3
      theparton(nn)%momentum(mu) = kf(j,mu)
    END DO
    theparton(nn)%masssq      = kfmasssq(j)
    theparton(nn)%qsqscale    = 777.0d7  ! dummy
    theparton(nn)%flavor      = theshower%ptn(n)%flavor
    theparton(nn)%self        = nn
    theparton(nn)%parent      = 0  ! Two partons will be left with parent = 0.
    theparton(nn)%stringquark = theshower%ptn(n)%stringquark
    theparton(nn)%stringqbar  = theshower%ptn(n)%stringqbar
    theparton(nn)%child1      = -1
    theparton(nn)%child2      = -1
    theparton(nn)%childless   = .true.
    theparton(nn)%done        = theshower%ptn(n)%done
    active(nn) = .true.
    IF (theparton(nn)%flavor.GT.0) THEN
      numberofquarks = numberofquarks + 1
    END IF
  END IF
END DO
!
! We will need the color partner of each quark (if numberofquarks > 1) in order to correctly
! do joining of quarks to antiquarks to make a gluon. The q and qbar to be joined cannot
! be color partners.
!
IF (numberofquarks.GT.1) THEN  ! find color connections
  DO n1 = nparticles-1, 2*nparticles-2
    IF (theparton(n1)%flavor.GT.0) THEN ! find partner(n1)
      tempstring = theparton(n1)%stringquark
      partnerfound = .false.
      count = 0
      DO  ! until partner found
        count = count + 1
        DO n2 = nparticles-1, 2*nparticles-2
          IF (theparton(n2)%stringqbar.EQ.tempstring) THEN
            IF (theparton(n2)%flavor.EQ.0) THEN
              tempstring = theparton(n2)%stringquark
            ELSE IF (theparton(n2)%flavor.lT.0) THEN
              partnerfound = .true.
              partner(n1) = n2
              partner(n2) = n1
            ELSE
              write(nout,*)'Trouble with qbar partner in translate'
              STOP
            END IF
          END IF
        END DO
        IF (partnerfound) EXIT
        IF (count.GT.100) THEN
           write(nout,*)'Real trouble with qbar partner in translate'
           STOP
        END IF
      END DO ! until partner found
    END IF   ! theparton(n1)%flavor.GT.0
  END DO     ! n1 = nparticles-1, 2*nparticles-2
END IF       ! (numberofquarks.GT.1)
!
! Calculate angle for angular ordering. We define this to for each final state parton n1
! to be the smallest angle to a color connected final state parton n2.
!
DO n1 = nparticles-1, 2*nparticles-2
  theparton(n1)%costheta = - 777.0d0 
  DO n2 = nparticles-1, 2*nparticles-2
    IF (n1.NE.n2) THEN
      ! Check for parton matching stringquark, compute if found.
      IF (theparton(n1)%stringquark.NE.0) THEN
        IF (theparton(n1)%stringquark.EQ.theparton(n2)%stringqbar) THEN
          dot12 = 0.0d0
          dot11 = 0.0d0
          dot22 = 0.0d0
          DO mu = 1,3
            dot12 = dot12 + &
               theparton(n1)%momentum(mu) * theparton(n2)%momentum(mu)
            dot11 = dot11 + theparton(n1)%momentum(mu)**2
            dot22 = dot22 + theparton(n2)%momentum(mu)**2
          END DO
          cos12 = dot12/sqrt(dot11*dot22)
          IF (cos12.GT.theparton(n1)%costheta) THEN
            theparton(n1)%costheta = cos12
          END IF
        END IF ! parton(n1)%stringquark.EQ.parton(n2)%stringqbar
      END IF   ! parton(n1)%stringquark.NE.0
      ! Check for parton matching stringqbar, compute if found.
      IF (theparton(n1)%stringqbar.NE.0) THEN
        IF (theparton(n1)%stringqbar.EQ.theparton(n2)%stringquark) THEN
          dot12 = 0.0d0
          dot11 = 0.0d0
          dot22 = 0.0d0
          DO mu = 1,3
            dot12 = dot12 + &
               theparton(n1)%momentum(mu) * theparton(n2)%momentum(mu)
            dot11 = dot11 + theparton(n1)%momentum(mu)**2
            dot22 = dot22 + theparton(n2)%momentum(mu)**2
          END DO
          cos12 = dot12/sqrt(dot11*dot22)
          IF (cos12.GT.theparton(n1)%costheta) THEN
            theparton(n1)%costheta = cos12
          END IF
        END IF ! parton(n1)%stringqbar.EQ.parton(n2)%stringquark
      END IF   ! parton(n1)%stringqbar.NE.0
    END IF  ! n1.NE.n2
  END DO    ! n2 loop
  IF (theparton(n1)%costheta.LT.(-1.0d0)) THEN
    write(nout,*) 'We failed to find the color connected parton',n1
    STOP
  END IF
END DO      ! n1 looop
!
! Now find the mothers according to the kt algorithm, Pythia version. 
! But two daughters can be joined only if the flavors are 
! (0,f2) or (f1,0) or (f,-f) with f.ne.0. and we do not combine the last quarks.
! In case of (0,f2) or (f1,0) the two partons must be color connected.
! N.B., the color string of the mother is the qbar string
! from n1 and the quark string from n2.
!
DO n1 = 1, theevent%length
DO n2 = 1, theevent%length
   yvalues(n1,n2) = 1.0d30     ! Initialize with big values.
   dot12values(n1,n2) = 1.0d30 ! Initialize with big values.
   new(n1,n2) = .true.         ! We need to calculate yvalues for this n1,n2.
END DO
END DO
!
mother = nparticles - 2 ! Index of the current mother parton to be made.
!
ymax = 0.0d0
!
DO ! while (mother.GT.0) This is the main loop for combining partons.
!
ymin = 1.0d30
found = .false.
DO n1 = 1,theevent%length
IF(active(n1)) THEN ! a parton in the jet list.
  DO n2 = 1,theevent%length
  IF (n2.NE.n1) THEN
  IF(active(n2)) THEN ! a second parton in the jet list.
    ok = .false.
    IF((theparton(n1)%flavor.EQ.0).OR.(theparton(n2)%flavor.EQ.0)) THEN 
      ! We have gg, quark-g, or qbar-g
      IF(theparton(n1)%stringquark.GT.0) THEN
        IF(theparton(n1)%stringquark.EQ.theparton(n2)%stringqbar) THEN 
          ! 1 and 2 are color connected
          ok = .true.      ! These two could be joined.
          found = .true.   ! We have found at least one joinable pair.
        END IF
      END IF
    ELSE IF (numberofquarks.GT.1) THEN ! There is still a q-qbar pairs.
      IF  ((theparton(n1)%flavor.LT.0).AND. &
           (theparton(n1)%flavor + theparton(n2)%flavor .EQ. 0))  THEN  
        ! n1 is a qbar  and n2 is matching quark
        IF (n1.NE.partner(n2)) THEN
        ! They are not color connected, even through gluon chain.
          ok = .true.    ! These two could be joined.
          found = .true.   ! We have found at least one joinable pair.
        END IF
      END IF
    END IF
    IF (ok) THEN
      IF (new(n1,n2)) THEN    ! we need to calculate for this combination
        new(n1,n2) = .false.  ! but we will not need to calculate again
        dot12 = 0.0d0
        dot11 = 0.0d0
        dot22 = 0.0d0
        DO mu = 1,3
          dot12 = dot12 + &
             theparton(n1)%momentum(mu) * theparton(n2)%momentum(mu)
          dot11 = dot11 + theparton(n1)%momentum(mu)**2
          dot22 = dot22 + theparton(n2)%momentum(mu)**2
        END DO
        cos12 = dot12/sqrt(dot11*dot22)
        temp = 1.0d0 - cos12
        temp = temp * 2.0d0 *                                            &
            (theparton(n1)%momentum(0)*theparton(n2)%momentum(0)         &
            /(theparton(n1)%momentum(0) + theparton(n2)%momentum(0)))**2
         ! This is the Pythia definition of kt^2
        y12 = temp / externalrts**2
        yvalues(n1,n2) = y12
        dot12values(n1,n2) = dot12
      ELSE ! we have already calculated for this combination
        y12 = yvalues(n1,n2)
      END IF
      IF ( y12.LT.ymin ) THEN
        ymin = y12
        n1star = n1
        n2star = n2
        dot12star = dot12values(n1,n2)
      END IF
    END IF
  END IF    ! active(n2)
  END IF    ! n2.NE.n1
  END DO    ! n2 = 1,theevent%length
END IF      ! active(n1)
END DO      ! n1 = 1,theevent%length
IF (.NOT.found) THEN
  write(nout,*)'No possible pairs to join found in translate'
  DO n1 = 1,theevent%length
    IF (active(n1)) THEN
      write(nout,*)n1,theparton(n1)%flavor,  &
           theparton(n1)%stringqbar,theparton(n1)%stringquark
    END IF
  END DO
  DO n1 = 1,theshower%length
  IF (theshower%ptn(n1)%childless) THEN
    write(nout,*)theshower%ptn(n1)%flavor,  &
      theshower%ptn(n1)%stringqbar,theshower%ptn(n1)%stringquark
  END IF
  END DO
  STOP
END IF
active(mother) = .true.   ! mother is an active jet
active(n1star) = .false.  ! erase n1star as an active jet
active(n2star) = .false.  ! erase n2star as an active jet
IF (numberofquarks.GT.1) THEN  ! Fix partner information.
  IF ( (theparton(n1star)%flavor.NE.0) .AND.    &
       (theparton(n2star)%flavor.EQ.0) ) THEN
    partner(mother) = partner(n1star)
    partner(partner(n1star)) = mother
  ELSE IF ( (theparton(n2star)%flavor.NE.0) .AND.    &
            (theparton(n1star)%flavor.EQ.0) ) THEN
    partner(mother) = partner(n2star)
    partner(partner(n2star)) = mother
  END IF
END IF
IF ((theparton(n1star)%flavor .NE. 0).AND.(theparton(n2star)%flavor .NE. 0)) THEN
  numberofquarks = numberofquarks - 1  ! We combined a q with a qbar.
END IF
!
! Create mother. Modify daughters.
!
DO mu = 0,3
  theparton(mother)%momentum(mu) =   &
        theparton(n1star)%momentum(mu) + theparton(n2star)%momentum(mu)
END DO
theparton(mother)%masssq  = theparton(n1star)%masssq                    &
                            + theparton(n2star)%masssq                  &
                            + 2.0d0*( theparton(n1star)%momentum(0) *   & 
                                      theparton(n2star)%momentum(0)     &
                                      - dot12star) 
theparton(mother)%flavor      = theparton(n1star)%flavor + & 
                                theparton(n2star)%flavor 
theparton(mother)%self        = mother
theparton(mother)%parent      = 0  ! Two partons will be left with parent = 0.
theparton(mother)%stringqbar  = theparton(n1star)%stringqbar
theparton(mother)%stringquark = theparton(n2star)%stringquark
theparton(mother)%child1      = n1star
theparton(mother)%child2      = n2star
theparton(mother)%childless   = .false.
theparton(mother)%done        = .false.
theparton(n1star)%parent      = mother
theparton(n2star)%parent      = mother
theparton(mother)%costheta    = - 1.0d0  ! will remain for intermediate partons
!
mother = mother - 1
IF (ymax.LT.ymin) THEN
  ymax = ymin        ! Largest of the y values at which partons are combined.
END IF
IF (mother.EQ.0) EXIT
END DO
!
! Find qsqscale for each parton.
!
count = 0
DO n = 1,theevent%length
  known(n) = .false.
  IF (theparton(n)%parent.EQ.0) THEN
    known(n) = .true.
    theparton(n)%qsqscale = externalrts**2  ! Start for recursion.
  END IF
END DO
DO  ! until done
count = count + 1
done = .true.
DO n = 1,theevent%length
  IF (.NOT.known(n)) THEN
    mother = theparton(n)%parent
    IF (known(mother)) THEN
      IF (theparton(mother)%flavor.NE.0 ) THEN  ! mother is not a gluon
        IF (theparton(n)%flavor.EQ.0) THEN
           theparton(n)%qsqscale = theparton(mother)%masssq    ! q --> g
        ELSE
           theparton(n)%qsqscale = theparton(mother)%qsqscale  ! q --> q
        END IF
      ELSE  ! mother is a gluon
        IF (theparton(n)%flavor.NE.0) THEN
          theparton(n)%qsqscale = theparton(mother)%masssq     ! g --> q
        ELSE                                                   ! g --> g
          IF (theparton(n)%momentum(0).GT. &
              0.5d0* theparton(mother)%momentum(0) ) THEN
             ! This is the daughter gluon with z > 1/2
             grandmother = theparton(mother)%parent
             theparton(n)%qsqscale = theparton(grandmother)%masssq
             ! theparton(n)%qsqscale = theparton(mother)%qsqscale !  ALTERNATIVE VERSION
          ELSE
             ! This is the daughter gluon with z < 1/2
             theparton(n)%qsqscale = theparton(mother)%masssq
          END IF
        END IF  ! for check on mother flavor
      END IF    ! mother is known
      known(n) = .true.
    ELSE
      ! mother was not known, so we will come back to this parton
      done = .false. 
    END IF  ! mother known
  END IF  ! parton n known
END DO ! n = 1,theevent%length
IF (done) EXIT
IF (count.GT.100) THEN
  ! Actually, mother < child1 and child2 so we never get count > 1.
  write(nout,*)'Things are not going well in translate.'
  STOP
END IF
END DO ! main loop until done
!
! Record the parton information
!
theevent%ktmax = sqrt(ymax)*externalrts ! Usually, the kt for 3 jets to 2 jets.
theevent%ptn = theparton                ! ptn and theparton are arrays
!
END subroutine translate
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine makekf(nparticles,kf)
!
use beowulf_parameters
implicit none
! Out:
integer :: nparticles
real(kind=dbl) :: kf(maxparticles,0:3)
!
! Produces list of final state particles kf from event in Pythia event
! common block.
!
! 12 October 2004
!
! ------- Pythia  event common block  -------
!
integer:: N, NPAD, K(4000,5)
real(kind=dbl) :: P(4000,5), V(4000,5)
common /PYJETS/ N, NPAD, K, P, V
save /PYJETS/
!
! -------------------------------------------
integer :: i,j
!
j = 0
DO i = 1,N
  IF (K(i,1).LE.2) THEN
    j = j+1
    IF (j.GT.maxparticles) THEN
      write(nout,*)'Subroutine makekf found too many particles',j
      write(nout,*)'Probably the solution is to increase maxparticles in beowulf_parameters'
      write(nout,*)'and recompile beowulf and beowulfsubs.'
      STOP
    END IF
    kf(j,1) = P(i,1)
    kf(j,2) = P(i,2)
    kf(j,3) = P(i,3)
    kf(j,0) = P(i,4)
  END IF
END DO
nparticles = j
RETURN
END subroutine makekf
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine makekfmod(theevent,nparticles,kf)
!
use beowulf_parameters
use pythiawulf
implicit none
! In:
type(eventdata) :: theevent
! Out:
integer :: nparticles
real(kind=dbl) :: kf(maxparticles,0:3)
!
! Produces list of final state particles kf from theevent input.
!
! 23 May 2005
!
! -------------------------------------------
integer :: i,j
!
j = 0
DO i = 1,theevent%length
  IF (theevent%ptn(i)%childless) THEN
    j = j+1
    IF (j.GT.maxparticles) THEN
      write(nout,*)'Subroutine makekfmod found too many particles',j
      write(nout,*)'Probably the solution is to increase maxparticles in beowulf_parameters'
      write(nout,*)'and recompile beowulf and beowulfsubs.'
      STOP
    END IF
    kf(j,1) = theevent%ptn(i)%momentum(1)
    kf(j,2) = theevent%ptn(i)%momentum(2)
    kf(j,3) = theevent%ptn(i)%momentum(3)
    kf(j,0) = theevent%ptn(i)%momentum(0)
  END IF
END DO
nparticles = j
RETURN
END subroutine makekfmod
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine daytime(date)
!
use beowulf_parameters
implicit none
! Out:
character(len=28) :: date
! 
! Returns DATE = today's date and time as a character string.
!
character*10 thedate,thetime
call date_and_time(thedate,thetime)
date = thedate//' '//thetime
RETURN
END subroutine daytime
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!
subroutine timing(deltatime)
!
use beowulf_parameters
implicit none
! Out:
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/
save t1
!
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_parameters
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_parameters
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_parameters
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