pyupev.f

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

 
C*********************************************************************
 
C...PYUPEV
C...Administers the hard-process generation required for output to the
C...Les Houches event record.
 
      SUBROUTINE pyupev
 
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)
      SAVE /pyjets/,/pyctag/,/pydat1/,/pydat2/,/pydat3/,/pypars/,
     &/pyint1/,/pyint2/,/pyint4/
 
C...HEPEUP for output.
      INTEGER maxnup
      parameter(maxnup=500)
      INTEGER nup,idprup,idup,istup,mothup,icolup
      DOUBLE PRECISION xwgtup,scalup,aqedup,aqcdup,pup,vtimup,spinup
      common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
     &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
     &vtimup(maxnup),spinup(maxnup)
      SAVE /hepeup/
 
C...Stop if no subprocesses on.
      IF(mint(121).EQ.1.AND.msti(53).EQ.1) THEN
        WRITE(mstu(11),5100)
        stop
      ENDIF
 
C...Special flags for hard-process generation only.
      mstp71=mstp(71)
      mstp(71)=0
      mst128=mstp(128)
      mstp(128)=1
 
C...Initial values for some counters.
      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/.
      mint(33)=0
 
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...Do not allow pileup events.
      mint(82)=1
 
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(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((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)
  110   mint(31)=mint31
        mint(51)=0
        CALL pyscat
        IF(mint(51).EQ.1) goto 100
        ipu1=mint(84)+1
        ipu2=mint(84)+2
 
C...Decay of final state resonances.
        mint(32)=0
        IF(mstp(41).GE.1.AND.iset(isub).LE.10.AND.isub.NE.95)
     &  CALL pyresd(0)
        IF(mint(51).EQ.1) goto 100
        mint(52)=n
 
C...Longitudinal boost of hard scattering.
        betaz=(vint(41)-vint(42))/(vint(41)+vint(42))
        CALL pyrobo(mint(84)+1,n,0d0,0d0,0d0,0d0,betaz)
 
      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
      nfix=n
      DO 120 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
  120 CONTINUE
 
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...Store event information and calculate Monte Carlo estimates of
C...subprocess cross-sections.
  130 CALL pydocu
 
C...Transform to the desired coordinate frame.
  140 CALL pyfram(mstp(124))
      mstu(70)=mstu70
      paru(21)=vint(1)
 
C...Restore special flags for hard-process generation only.
      mstp(71)=mstp71
      mstp(128)=mst128
 
C...Trace colour tags; convert to LHA style labels.
      nct=100
      DO 150 i=mint(84)+1,n
        mct(i,1)=0
        mct(i,2)=0
  150 CONTINUE
      DO 160 i=mint(84)+1,n
        kq=kchg(pycomp(k(i,2)),2)*isign(1,k(i,2))
        IF(k(i,1).EQ.3.OR.k(i,1).EQ.13.OR.k(i,1).EQ.14) THEN
          IF(k(i,4).NE.0.AND.(kq.EQ.1.OR.kq.EQ.2).AND.mct(i,1).EQ.0)
     &    THEN
            imo=mod(k(i,4)/mstu(5),mstu(5))
            ida=mod(k(i,4),mstu(5))
            IF(imo.NE.0.AND.mod(k(imo,5)/mstu(5),mstu(5)).EQ.i.AND.
     &      mct(imo,2).NE.0) THEN
              mct(i,1)=mct(imo,2)
            ELSEIF(imo.NE.0.AND.mod(k(imo,4),mstu(5)).EQ.i.AND.
     &      mct(imo,1).NE.0) THEN
              mct(i,1)=mct(imo,1)
            ELSEIF(ida.NE.0.AND.mod(k(ida,5),mstu(5)).EQ.i.AND.
     &      mct(ida,2).NE.0) THEN
              mct(i,1)=mct(ida,2)
            ELSE
              nct=nct+1
              mct(i,1)=nct
            ENDIF
          ENDIF
          IF(k(i,5).NE.0.AND.(kq.EQ.-1.OR.kq.EQ.2).AND.mct(i,2).EQ.0)
     &    THEN
            imo=mod(k(i,5)/mstu(5),mstu(5))
            ida=mod(k(i,5),mstu(5))
            IF(imo.NE.0.AND.mod(k(imo,4)/mstu(5),mstu(5)).EQ.i.AND.
     &      mct(imo,1).NE.0) THEN
              mct(i,2)=mct(imo,1)
            ELSEIF(imo.NE.0.AND.mod(k(imo,5),mstu(5)).EQ.i.AND.
     &      mct(imo,2).NE.0) THEN
              mct(i,2)=mct(imo,2)
            ELSEIF(ida.NE.0.AND.mod(k(ida,4),mstu(5)).EQ.i.AND.
     &      mct(ida,1).NE.0) THEN
              mct(i,2)=mct(ida,1)
            ELSE
              nct=nct+1
              mct(i,2)=nct
            ENDIF
          ENDIF
        ENDIF
  160 CONTINUE
 
C...Put event in HEPEUP commonblock.
      nup=n-mint(84)
      idprup=mint(1)
      xwgtup=1d0
      scalup=vint(53)
      aqedup=vint(57)
      aqcdup=vint(58)
      DO 180 i=1,nup
        idup(i)=k(i+mint(84),2)
        IF(i.LE.2) THEN
          istup(i)=-1
          mothup(1,i)=0
          mothup(2,i)=0
        ELSEIF(k(i+4,3).EQ.0) THEN
          istup(i)=1
          mothup(1,i)=1
          mothup(2,i)=2
        ELSE
          istup(i)=1
          mothup(1,i)=k(i+mint(84),3)-mint(84)
          mothup(2,i)=0
        ENDIF
        IF(i.GE.3.AND.k(i+mint(84),3).GT.0)
     &  istup(k(i+mint(84),3)-mint(84))=2
        icolup(1,i)=mct(i+mint(84),1)
        icolup(2,i)=mct(i+mint(84),2)
        DO 170 j=1,5
          pup(j,i)=p(i+mint(84),j)
  170   CONTINUE
        vtimup(i)=v(i,5)
        spinup(i)=9d0
  180 CONTINUE
 
C...Optionally write out event to disk. Minimal size for time/spin fields.
      IF(mstp(162).GT.0) THEN
        WRITE(mstp(162),5200) nup,idprup,xwgtup,scalup,aqedup,aqcdup
        DO 190 i=1,nup
          IF(vtimup(i).EQ.0d0) THEN
            WRITE(mstp(162),5300) idup(i),istup(i),mothup(1,i),
     &      mothup(2,i),icolup(1,i),icolup(2,i),(pup(j,i),j=1,5),
     &      ' 0. 9.'
          ELSE
            WRITE(mstp(162),5400) idup(i),istup(i),mothup(1,i),
     &      mothup(2,i),icolup(1,i),icolup(2,i),(pup(j,i),j=1,5),
     &      vtimup(i),' 9.'
          ENDIF
  190   CONTINUE

C...Optional extra line with parton-density information.
        IF(mstp(165).GE.1) WRITE(mstp(162),5500) msti(15),msti(16),
     &  pari(33),pari(34),pari(23),pari(29),pari(30) 
      ENDIF
 
C...Error messages and other print formats.
 5100 FORMAT(1x,'Error: no subprocess switched on.'/
     &1x,'Execution stopped.')
 5200 FORMAT(1p,2i6,4e14.6)
 5300 FORMAT(1p,i8,5i5,5e18.10,a6)
 5400 FORMAT(1p,i8,5i5,5e18.10,e12.4,a3)
 5500 FORMAT(1p,'#pdf ',2i5,5e18.10)
 
      RETURN
      END