pyevnt.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyevnt.f

 
C*********************************************************************
 
C...PYEVNT
C...Administers the generation of a high-pT event via calls to
C...a number of subroutines.
 
      SUBROUTINE pyevnt
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pyctag/nct,mct(4000,2)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
      common/pyint4/mwid(500),wids(500,5)
      common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
      SAVE /pyjets/,/pydat1/,/pyctag/,/pydat2/,/pydat3/,/pypars/,
     &/pyint1/,/pyint2/,/pyint4/,/pyint5/
C...Local array.
      dimension vtx(4)
 
C...Optionally let PYEVNW do the whole job.
      IF(mstp(81).GE.20) THEN
        CALL pyevnw
        RETURN
      ENDIF
 
C...Stop if no subprocesses on.
      IF(mint(121).EQ.1.AND.msti(53).EQ.1) THEN
        WRITE(mstu(11),5100)
        CALL pystop(1)
      ENDIF
 
C...Initial values for some counters.
      mstu(1)=0
      mstu(2)=0
      n=0
      mint(5)=mint(5)+1
      mint(7)=0
      mint(8)=0
      mint(30)=0
      mint(83)=0
      mint(84)=mstp(126)
      mstu(24)=0
      mstu70=0
      mstj14=mstj(14)
C...Normally, use K(I,4:5) colour info rather than /PYCTAG/.
      nct=0
      mint(33)=0
 
C...Let called routines know call is from PYEVNT (not PYEVNW).
      mint(35)=1
      IF (mstp(81).GE.10) mint(35)=2
 
C...If variable energies: redo incoming kinematics and cross-section.
      msti(61)=0
      IF(mstp(171).EQ.1) THEN
        CALL pyinki(1)
        IF(msti(61).EQ.1) THEN
          mint(5)=mint(5)-1
          RETURN
        ENDIF
        IF(mint(121).GT.1) CALL pysave(3,1)
        CALL pyxtot
      ENDIF
 
C...Loop over number of pileup events; check space left.
      IF(mstp(131).LE.0) THEN
        npile=1
      ELSE
        CALL pypile(2)
        npile=mint(81)
      ENDIF
      DO 270 ipile=1,npile
        IF(mint(84)+100.GE.mstu(4)) THEN
          CALL pyerrm(11,
     &    '(PYEVNT:) no more space in PYJETS for pileup events')
          IF(mstu(21).GE.1) goto 280
        ENDIF
        mint(82)=ipile
 
C...Generate variables of hard scattering.
        mint(51)=0
        msti(52)=0
  100   CONTINUE
        IF(mint(51).NE.0.OR.mstu(24).NE.0) msti(52)=msti(52)+1
        mint(31)=0
        mint(39)=0
        mint(51)=0
        mint(57)=0
        CALL pyrand
        IF(msti(61).EQ.1) THEN
          mint(5)=mint(5)-1
          RETURN
        ENDIF
        IF(mint(51).EQ.2) RETURN
        isub=mint(1)
        IF(mstp(111).EQ.-1) goto 260
 
C...Loopback point if PYPREP fails, especially for junction topologies.
        nprep=0
        mnt31s=mint(31)
  110   nprep=nprep+1
        mint(31)=mnt31s
 
        IF((isub.LE.90.OR.isub.GE.95).AND.isub.NE.99) THEN
C...Hard scattering (including low-pT):
C...reconstruct kinematics and colour flow of hard scattering.
          mint31=mint(31)
  120     mint(31)=mint31
          mint(51)=0
          CALL pyscat
          IF(mint(51).EQ.1) goto 100
          ipu1=mint(84)+1
          ipu2=mint(84)+2
          IF(isub.EQ.95) goto 140
 
C...Reset statistics on activity in event.
        DO 130 j=351,359
          mint(j)=0
          vint(j)=0d0
  130   CONTINUE
 
C...Showering of initial state partons (optional).
          nfin=n
          alamsv=parj(81)
          parj(81)=parp(72)
          IF(mstp(61).GE.1.AND.mint(47).GE.2.AND.mint(111).NE.12)
     &    CALL pysspa(ipu1,ipu2)
          parj(81)=alamsv
          IF(mint(51).EQ.1) goto 100
 
C...Showering of final state partons (optional).
          alamsv=parj(81)
          parj(81)=parp(72)
          IF(mstp(71).GE.1.AND.iset(isub).GE.2.AND.iset(isub).LE.10)
     &    THEN
            ipu3=mint(84)+3
            ipu4=mint(84)+4
            IF(iset(isub).EQ.5) ipu4=-3
            qmax=vint(55)
            IF(iset(isub).EQ.2) qmax=sqrt(parp(71))*vint(55)
            CALL pyshow(ipu3,ipu4,qmax)
          ELSEIF(iset(isub).EQ.11) THEN
            CALL pyadsh(nfin)
          ENDIF
          parj(81)=alamsv
 
C...Allow possibility for user to abort event generation.
          iveto=0
          IF(ipile.EQ.1.AND.mstp(143).EQ.1) CALL pyveto(iveto)
          IF(iveto.EQ.1) goto 100
 
C...Decay of final state resonances.
          mint(32)=0
          IF(mstp(41).GE.1.AND.iset(isub).LE.10) CALL pyresd(0)
          IF(mint(51).EQ.1) goto 100
          mint(52)=n
 
 
C...Multiple interactions - PYTHIA 6.3 intermediate style.
  140     IF(mstp(81).GE.10.AND.mint(50).EQ.1) THEN
            IF(isub.EQ.95) mint(31)=mint(31)+1
            CALL pymign(6)
            IF(mint(51).EQ.1) goto 100
            mint(53)=n
 
C...Beam remnant flavour and colour assignments - new scheme.
            CALL pymihk
            IF(mint(51).EQ.1.AND.mint(57).GE.1.AND.mint(57).LE.5)
     &      goto 120
            IF(mint(51).EQ.1) goto 100
 
C...Primordial kT and beam remnant momentum sharing - new scheme.
            CALL pymirm
            IF(mint(51).EQ.1.AND.mint(57).GE.1.AND.mint(57).LE.5)
     &      goto 120
            IF(mint(51).EQ.1) goto 100
            IF(isub.EQ.95) mint(31)=mint(31)-1
 
C...Multiple interactions - PYTHIA 6.2 style.
          ELSEIF(mint(111).NE.12) THEN
            IF (mstp(81).GE.1.AND.mint(50).EQ.1.AND.isub.NE.95) THEN
              CALL pymult(6)
              mint(53)=n
            ENDIF
 
C...Hadron remnants and primordial kT.
            CALL pyremn(ipu1,ipu2)
            IF(mint(51).EQ.1.AND.mint(57).GE.1.AND.mint(57).LE.5) goto
     &           110
            IF(mint(51).EQ.1) goto 100
          ENDIF
 
        ELSEIF(isub.NE.99) THEN
C...Diffractive and elastic scattering.
          CALL pydiff
 
        ELSE
C...DIS scattering (photon flux external).
          CALL pydisg
          IF(mint(51).EQ.1) goto 100
        ENDIF
 
C...Check that no odd resonance left undecayed.
        mint(54)=n
        IF(mstp(111).GE.1) THEN
          nfix=n
          DO 150 i=mint(84)+1,nfix
            IF(k(i,1).GE.1.AND.k(i,1).LE.10.AND.k(i,2).NE.21.AND.
     &      k(i,2).NE.22) THEN
              kca=pycomp(k(i,2))
              IF(mwid(kca).NE.0.AND.mdcy(kca,1).GE.1) THEN
                CALL pyresd(i)
                IF(mint(51).EQ.1) goto 100
              ENDIF
            ENDIF
  150     CONTINUE
        ENDIF
 
C...Boost hadronic subsystem to overall rest frame.
C..(Only relevant when photon inside lepton beam.)
        IF(mint(141).NE.0.OR.mint(142).NE.0) CALL pygaga(4,wtgaga)
 
C...Recalculate energies from momenta and masses (if desired).
        IF(mstp(113).GE.1) THEN
          DO 160 i=mint(83)+1,n
            IF(k(i,1).GT.0.AND.k(i,1).LE.10) p(i,4)=sqrt(p(i,1)**2+
     &      p(i,2)**2+p(i,3)**2+p(i,5)**2)
  160     CONTINUE
          nrecal=n
        ENDIF
 
C...Colour reconnection before string formation
        IF (mstp(95).GE.2) CALL pyfscr(mint(84)+1)

C...Rearrange partons along strings, check invariant mass cuts.
        mstu(28)=0
        IF(mstp(111).LE.0) mstj(14)=-1
        CALL pyprep(mint(84)+1)
        mstj(14)=mstj14
        IF(mint(51).EQ.1.AND.mstu(24).EQ.1) THEN
          mstu(24)=0
          goto 100
        ENDIF
        IF (mint(51).EQ.1.AND.nprep.LE.5) goto 110
        IF (mint(51).EQ.1) goto 100
        IF(mstp(112).EQ.1.AND.mstu(28).EQ.3) goto 100
        IF(mstp(125).EQ.0.OR.mstp(125).EQ.1) THEN
          DO 190 i=mint(84)+1,n
            IF(k(i,2).EQ.94) THEN
              DO 180 i1=i+1,min(n,i+10)
                IF(k(i1,3).EQ.i) THEN
                  k(i1,3)=mod(k(i1,4)/mstu(5),mstu(5))
                  IF(k(i1,3).EQ.0) THEN
                    DO 170 ii=mint(84)+1,i-1
                        IF(k(ii,2).EQ.k(i1,2)) THEN
                          IF(mod(k(ii,4),mstu(5)).EQ.i1.OR.
     &                    mod(k(ii,5),mstu(5)).EQ.i1) k(i1,3)=ii
                        ENDIF
  170               CONTINUE
                    IF(k(i+1,3).EQ.0) k(i+1,3)=k(i,3)
                  ENDIF
                ENDIF
  180         CONTINUE
            ENDIF
  190     CONTINUE
          CALL pyedit(12)
          CALL pyedit(14)
          IF(mstp(125).EQ.0) CALL pyedit(15)
          IF(mstp(125).EQ.0) mint(4)=0
          DO 210 i=mint(83)+1,n
            IF(k(i,1).EQ.11.AND.k(i,4).EQ.0.AND.k(i,5).EQ.0) THEN
              DO 200 i1=i+1,n
                IF(k(i1,3).EQ.i.AND.k(i,4).EQ.0) k(i,4)=i1
                IF(k(i1,3).EQ.i) k(i,5)=i1
  200         CONTINUE
            ENDIF
  210     CONTINUE
        ENDIF
 
C...Introduce separators between sections in PYLIST event listing.
        IF(ipile.EQ.1.AND.mstp(125).LE.0) THEN
          mstu70=1
          mstu(71)=n
        ELSEIF(ipile.EQ.1) THEN
          mstu70=3
          mstu(71)=2
          mstu(72)=mint(4)
          mstu(73)=n
        ENDIF
 
C...Go back to lab frame (needed for vertices, also in fragmentation).
        CALL pyfram(1)
 
C...Set nonvanishing production vertex (optional).
        IF(mstp(151).EQ.1) THEN
          DO 220 j=1,4
            vtx(j)=parp(150+j)*sqrt(-2d0*log(max(1d-10,pyr(0))))*
     &      sin(paru(2)*pyr(0))
  220     CONTINUE
          DO 240 i=mint(83)+1,n
            DO 230 j=1,4
              v(i,j)=v(i,j)+vtx(j)
  230       CONTINUE
  240     CONTINUE
        ENDIF
 
C...Perform hadronization (if desired).
        IF(mstp(111).GE.1) THEN
          CALL pyexec
          IF(mstu(24).NE.0) goto 100
        ENDIF
        IF(mstp(113).GE.1) THEN
          DO 250 i=nrecal,n
            IF(p(i,5).GT.0d0) p(i,4)=sqrt(p(i,1)**2+
     &      p(i,2)**2+p(i,3)**2+p(i,5)**2)
  250     CONTINUE
        ENDIF
        IF(mstp(125).EQ.0.OR.mstp(125).EQ.1) CALL pyedit(14)
 
C...Store event information and calculate Monte Carlo estimates of
C...subprocess cross-sections.
  260   IF(ipile.EQ.1) CALL pydocu
 
C...Set counters for current pileup event and loop to next one.
        msti(41)=ipile
        IF(ipile.GE.2.AND.ipile.LE.10) msti(40+ipile)=isub
        IF(mstu70.LT.10) THEN
          mstu70=mstu70+1
          mstu(70+mstu70)=n
        ENDIF
        mint(83)=n
        mint(84)=n+mstp(126)
        IF(ipile.LT.npile) CALL pyfram(2)
  270 CONTINUE
 
C...Generic information on pileup events. Reconstruct missing history.
      IF(mstp(131).EQ.1.AND.mstp(133).GE.1) THEN
        pari(91)=vint(132)
        pari(92)=vint(133)
        pari(93)=vint(134)
        IF(mstp(133).GE.2) pari(93)=pari(93)*xsec(0,3)/vint(131)
      ENDIF
      CALL pyedit(16)
 
C...Transform to the desired coordinate frame.
  280 CALL pyfram(mstp(124))
      mstu(70)=mstu70
      paru(21)=vint(1)
 
C...Error messages
 5100 FORMAT(1x,'Error: no subprocess switched on.'/
     &1x,'Execution stopped.')
 
      RETURN
      END