! ------------------------ ! 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 = (1/nreno) sum_j w_j S_j ! ! We now accept events with probability p_j, so that ! ! = (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 ! ! = (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 = (1/nreno) sum_j w_j S_j ! ! We now accept events with probability p_j, so that ! ! = (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 ! ! = (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, ~ 0.80, maxerror = 0.0002 ! imax = 1, ~ 0.65, maxerror = 0.004 ! imax = 1, ~ 0.55, maxerror = 0.02 ! imax = 3, ~ 0.55, maxerror = 0.003 ! imax = 10, ~ 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