! ------------------------ ! 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