! Main program for running Pythia with beowulf. DES
! 4 September 2005.
!
implicit none
integer, parameter :: dbl = kind(1.0d0)
!
!---------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/
!
! ------Pythia common block for initial data ------------------
!
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
!
! ------- 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/
!
!---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/
!
! ------- Pythia switches and parameters -------
!
integer :: MSTU(200), MSTJ(200)
real(kind=dbl) :: PARU(200), PARJ(200)
common /PYDAT1/ MSTU, PARU, MSTJ, PARJ
!
integer :: MSTP(200), MSTI(200)
real(kind=dbl) :: PARP(200), PARI(200)
common /PYPARS/ MSTP, PARP, MSTI, PARI
!
! ------- beowulf counting variable nreno -------
!
integer :: beowulfnreno
common /beowulfcounting/ beowulfnreno
!
! ---- for testing and reporting sample events
integer, parameter :: maxparticles = 512
integer :: nparticles
real(kind=dbl) :: kf(maxparticles,0:3) ! final state momenta
real(kind=dbl) :: y43,y32
real(kind=dbl) :: jetmass(3)
!
! Color factors:
real(kind=dbl) :: nc,nf,cf
common /colorfactors/ nc,nf,cf
! Physics data:
real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts
common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts
! Renormaliation parameters:
real(kind=dbl) :: muoverrts
common /renormalize/ muoverrts
real(kind=dbl) :: testsum
!
! Mrenna variables
integer :: IFSOPT
!
!--
integer :: naccepted,indexfactor
!
WRITE(*,*) ' This is the main program for Pythia and beowulf.'
!
! Set parameters for Pythia or ask that parameters be reported.
!
IFSOPT=0 ! 0 for ordinary shower, 1 for pt ordered shower.
! Turn off beam structure: e+ e- have full beam energy
call pygive('mstp(11)=0')
IF(IFSOPT.EQ.1) THEN ! Final state shower pt ordered.
call pygive('MSTP(3)=1')
call pygive('PARJ(81)=.14')
call pygive('PARP(72)=.14')
call pygive('PARP(61)=.192')
call pygive('MSTJ(41)=11')
ELSE ! Ordinary final state shower.
call pygive('MSTP(3)=')
! call pygive('PARJ(81)=.192')
! call pygive('PARP(72)=.192')
! call pygive('PARP(61)=.192')
call pygive('PARJ(81)=')
call pygive('PARP(72)=')
call pygive('PARP(61)=')
call pygive('MSTJ(41)=')
ENDIF
!
call pygive('MSTJ(1)=') ! Default is 1 for string fragmentation.
! Use 0 for no fragmentation.
call pygive('MSTP(81)=')
call pygive('MSTP(91)=')
call pygive('MSTP(61)=')
call pygive('MSTP(71)=')
call pygive('MSTP(111)=')
call pygive('MDCY(C111,1)=')
call pygive('MDCY(C211,1)=')
call pygive('MSTU(46)=')
call pygive('PARU(44)=')
call pygive('MSTU(22)=20')
! switch kind of shower
call pygive('MSTJ(47)=0')
call pygive('MSTP(123)=1')
call pygive('MSTU(21)=1')
! call pygive('MSTP(125)=2')
! options for clustering routine
call pygive('MSTU(41)=1') ! use all particles in pyclus
call pygive('PARU(44)=500.0') ! maximum distance for combining clusters
call pygive('MSTU(47)=3')
!
! Call Pythia initialization.
!
call PYINIT('USER',' ',' ',0.0d0)
!
! Ready to start.
!
indexfactor = 1000 ! This just controls progress notification
naccepted = 0
testsum = 0.0D0
call usrinit('compare.hb') ! Initialization of interface to beowulf.
!
! The main event loop.
!
DO ! Loop ends when UPEVNT sends NUP = 0 to Phthia.
call PYEVNT
IF (NUP.EQ.0) EXIT
call usrevnt(xwgtup) ! Pythia routine. Eventually the beowulf routine upevnt
! will supply the event.
naccepted = naccepted + 1
IF (mod(naccepted,indexfactor).EQ.0) THEN
write(*,*)'Event number',naccepted,' nup =',nup
IF (naccepted.GE.3000) THEN
indexfactor = 100000
END IF
END IF
IF(naccepted.LE.0) THEN ! Print out information on first events if you want.
write(*,*)' '
write(*,*) 'Pythia receives event',naccepted
write(*,*) 'With ',NUP,'total entries'
write(*,*) 'and ktlimit =',KTLIMIT
call makekf(nparticles,kf) ! Gets info from Pythia common block/PYJETS/.
call findjets(kf,nparticles,'durham',y43,y32,jetmass)
write(*,*) 'After Pythia, kt3(Durham) =',sqrt(y32*4.0d0*EBMUP(1)*EBMUP(2))
call findjets(kf,nparticles,'pythia',y43,y32,jetmass)
write(*,*) 'After Pythia, kt3(Pythia) =',sqrt(y32*4.0d0*EBMUP(1)*EBMUP(2))
call PYLIST(7) ! In case you want a sample event at input
call PYLIST(2) ! In case you want a sample event completed
END IF
!
! Event analysis could go here.
! The call to hrothgarII tells beowulf to do its own event analysis.
!
call hrothgarII(0,'newresult ')
!
END DO
call usrdone(dble(beowulfnreno))
! At this point, beowulfnreno is available from common /beowulfcounting/.
! This would be needed to complete any event analysis:
! <S> = (1/beowulfnreno) sum_(i=1)^Ntot XWGTUP_i S_i.
!
write(*,*)'Pythia done after',naccepted,' events.'
END
!
! ----------------Date and TIme -------------------
!
subroutine idate(idtemp)
integer :: idtemp(3)
integer :: ival(8)
!
CALL DATE_AND_TIME(VALUES=IVAL)
idtemp(3) = ival(1)
idtemp(2) = ival(2)
idtemp(1) = ival(3)
return
end
!
subroutine itime(idtemp)
integer :: idtemp(3)
integer :: ival(8)
!
CALL DATE_AND_TIME(VALUES=IVAL)
idtemp(1) = ival(5)
idtemp(2) = ival(6)
idtemp(3) = ival(7)
return
end
!
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!2345678901234567890123456789012345678901234567890123456789012345678901234567890