pyevnw.f

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

 
C*********************************************************************
 
C...PYEVNW
C...Administers the generation of a high-pT event via calls to
C...a number of subroutines for the new multiple interactions and
C...showering framework.
 
      SUBROUTINE pyevnw
 
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)
      common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
     &     xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
     &     xmi(2,240),pt2mi(240),imisep(0:240)
      SAVE /pyjets/,/pyctag/,/pydat1/,/pydat2/,/pydat3/,
     &     /pypars/,/pyint1/,/pyint2/,/pyint4/,/pyint5/,/pyintm/
C...Local arrays.
      dimension vtx(4)
 
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 /PYCT/.
      nct=0
      mint(33)=0
 
C...Let called routines know call is from PYEVNW (not PYEVNT).
      mint(35)=3
 
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 300 ipile=1,npile
        IF(mint(84)+100.GE.mstu(4)) THEN
          CALL pyerrm(11,
     &    '(PYEVNW:) no more space in PYJETS for pileup events')
          IF(mstu(21).GE.1) goto 310
        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(36)=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 290
 
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
          npartd=n
          nfin=n
 
C...Intertwined initial state showers and multiple interactions.
C...Force no IS showers if no pdfs defined: MSTP(61) -> 0 for PYEVOL.
C...Force no MI if cross section not known: MSTP(81) -> 0 for PYEVOL.
          mstp61=mstp(61)
          IF (mint(47).LT.2) mstp(61)=0
          mstp81=mstp(81)
          IF (mint(50).EQ.0) mstp(81)=0
          IF ((mstp(61).GE.1.OR.mod(mstp(81),10).GE.0).AND.
     &    mint(111).NE.12) THEN
C...Absolute max pT2 scale for evolution: phase space limit.
            pt2mxs=0.25d0*vint(2)
C...Check if more constrained by ISR and MI max scales:
            pt2mxs=min(pt2mxs,max(vint(56),vint(62)))
C...Loopback point in case of failure in evolution.
            loop=0
  130       loop=loop+1
            mint(51)=0
            IF(loop.GT.100) THEN
              CALL pyerrm(9,'(PYEVNW:) failed to evolve shower or '
     &             //'multiple interactions.')
              mint(51)=1
              RETURN
            ENDIF
 
C...Pre-initialization of interleaved MI/ISR/JI evolution, only done
C...once per event. (E.g. compute constants and save variables to be
C...restored later in case of failure.)
            IF (loop.EQ.1) CALL pyevol(-1,dummy1,dummy2)
 
C...Initialize interleaved MI/ISR/JI evolution.
C...PT2MAX: absolute upper limit for evolution - Initialization may
C...        return a PT2MAX which is lower than this.
C...PT2MIN: absolute lower limit for evolution - Initialization may
C...        return a PT2MIN which is larger than this (e.g. Lambda_QCD).
            pt2max=pt2mxs
            pt2min=0d0
            CALL pyevol(0,pt2max,pt2min)
            IF (mint(51).EQ.1) goto 130
 
C...Perform interleaved MI/ISR/JI evolution from PT2MAX to PT2MIN.
C...In principle factorized, so can be stopped and restarted.
C...Example: stop/start at pT=10 GeV. (Commented out for now.)
C            PT2MED=MAX(10D0**2,PT2MIN)
C            CALL PYEVOL(1,PT2MAX,PT2MED)
C            IF (MINT(51).EQ.1) GOTO 160
C            PT2MAX=PT2MED
            CALL pyevol(1,pt2max,pt2min)
            IF (mint(51).EQ.1) goto 130
 
C...Finalize interleaved MI/ISR/JI evolution.
            CALL pyevol(2,pt2max,pt2min)
            IF (mint(51).EQ.1) goto 130
 
          ENDIF
          mstp(61)=mstp61
          mstp(81)=mstp81
          IF(mint(51).EQ.1) goto 100
C...(MINT(52) is actually obsolete in this routine. Set anyway
C...to ensure PYDOCU stable.)
          mint(52)=n
          mint(53)=n
 
C...Beam remnants - new scheme.
  140     IF(mint(50).EQ.1) THEN
            IF (isub.EQ.95) mint(31)=1
 
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)=0
          ELSEIF(mint(111).NE.12) THEN
C...Hadron remnants and primordial kT - old model.
C...Happens e.g. for direct photon on one side.
            ipu1=imi(1,1,1)
            ipu2=imi(2,1,1)
            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
C...PYREMN does not set colour tags for BRs, so needs to be done now.
            DO 160 i=mint(53)+1,n
              DO 150 kcs=4,5
                ida=mod(k(i,kcs),mstu(5))
                IF (ida.NE.0) THEN
                  mct(i,kcs-3)=mct(ida,6-kcs)
                ELSE
                  mct(i,kcs-3)=0
                ENDIF
  150         CONTINUE
  160       CONTINUE
C...Instruct PYPREP to use colour tags
            mint(33)=1

            DO 360 mqgst=1,2
              DO 350 i=mint(84)+1,n
  
C...Look for coloured string endpoint, or (later) leftover gluon.
                IF (k(i,1).NE.3) goto 350
                kc=pycomp(k(i,2))
                IF(kc.EQ.0) goto 350
                kq=kchg(kc,2)
                IF(kq.EQ.0.OR.(mqgst.EQ.1.AND.kq.EQ.2)) goto 350
  
C...  Pick up loose string end with no previous tag.
                kcs=4
                IF(kq*isign(1,k(i,2)).LT.0) kcs=5
                IF(mct(i,kcs-3).NE.0) goto 350
                  
                CALL pycttr(i,kcs,i)
                IF(mint(51).NE.0) RETURN
  
 350          CONTINUE
 360        CONTINUE
C...Now delete any colour processing information if set (since partons
C...otherwise not FS showered!)
            DO 170 i=mint(84)+1,n
              IF (i.LE.n) THEN
                k(i,4)=mod(k(i,4),mstu(5)**2)
                k(i,5)=mod(k(i,5),mstu(5)**2)
              ENDIF
  170       CONTINUE
          ENDIF
 
C...Showering of final state partons (optional).
          alamsv=parj(81)
          parj(81)=parp(72)
          IF(mstp(71).GE.1.AND.iset(isub).GE.1.AND.iset(isub).LE.10)
     &    THEN
            qmax=vint(55)
            IF(iset(isub).EQ.2) qmax=sqrt(parp(71))*vint(55)
            CALL pyptfs(1,qmax,0d0,ptgen)
C...External processes: handle successive showers.
          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) ! sm
          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) THEN
            CALL pyresd(0)
            IF(mint(51).NE.0) goto 100
          ENDIF
 
          IF(mint(51).EQ.1) goto 100
 
        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 180 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
  180     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 190 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)
  190     CONTINUE
          nrecal=n
        ENDIF
 
C...Colour reconnection before string formation
        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) goto 110
        IF(mstp(112).EQ.1.AND.mstu(28).EQ.3) goto 100
        IF(mstp(125).EQ.0.OR.mstp(125).EQ.1) THEN
          DO 220 i=mint(84)+1,n
            IF(k(i,2).EQ.94) THEN
              DO 210 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 200 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
  200               CONTINUE
                    IF(k(i+1,3).EQ.0) k(i+1,3)=k(i,3)
                  ENDIF
                ENDIF
  210         CONTINUE
            ENDIF
  220     CONTINUE
          CALL pyedit(12)
          CALL pyedit(14)
          IF(mstp(125).EQ.0) CALL pyedit(15)
          IF(mstp(125).EQ.0) mint(4)=0
          DO 240 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 230 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
  230         CONTINUE
            ENDIF
  240     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 250 j=1,4
            vtx(j)=parp(150+j)*sqrt(-2d0*log(max(1d-10,pyr(0))))*
     &      sin(paru(2)*pyr(0))
  250     CONTINUE
          DO 270 i=mint(83)+1,n
            DO 260 j=1,4
              v(i,j)=v(i,j)+vtx(j)
  260       CONTINUE
  270     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 280 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)
  280     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.
  290   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)
  300 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.
  310 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