pyeevt.f

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

 
C*********************************************************************
 
C...PYEEVT
C...Handles the generation of an e+e- annihilation jet event.
 
      SUBROUTINE pyeevt(KFL,ECM)
 
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/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pyjets/,/pydat1/,/pydat2/
 
C...Check input parameters.
      IF(mstu(12).NE.12345) CALL pylist(0)
      IF(kfl.LT.0.OR.kfl.GT.8) THEN
        CALL pyerrm(16,'(PYEEVT:) called with unknown flavour code')
        IF(mstu(21).GE.1) RETURN
      ENDIF
      IF(kfl.LE.5) ecmmin=parj(127)+2.02d0*parf(100+max(1,kfl))
      IF(kfl.GE.6) ecmmin=parj(127)+2.02d0*pmas(kfl,1)
      IF(ecm.LT.ecmmin) THEN
        CALL pyerrm(16,'(PYEEVT:) called with too small CM energy')
        IF(mstu(21).GE.1) RETURN
      ENDIF
 
C...Check consistency of MSTJ options set.
      IF(mstj(109).EQ.2.AND.mstj(110).NE.1) THEN
        CALL pyerrm(6,
     &  '(PYEEVT:) MSTJ(109) value requires MSTJ(110) = 1')
        mstj(110)=1
      ENDIF
      IF(mstj(109).EQ.2.AND.mstj(111).NE.0) THEN
        CALL pyerrm(6,
     &  '(PYEEVT:) MSTJ(109) value requires MSTJ(111) = 0')
        mstj(111)=0
      ENDIF
 
C...Initialize alpha_strong and total cross-section.
      mstu(111)=mstj(108)
      IF(mstj(108).EQ.2.AND.(mstj(101).EQ.0.OR.mstj(101).EQ.1))
     &mstu(111)=1
      paru(112)=parj(121)
      IF(mstu(111).EQ.2) paru(112)=parj(122)
      IF(mstj(116).GT.0.AND.(mstj(116).GE.2.OR.abs(ecm-parj(151)).GE.
     &parj(139).OR.10*mstj(102)+kfl.NE.mstj(119))) CALL pyxtee(kfl,ecm,
     &xtot)
      IF(mstj(116).GE.3) mstj(116)=1
      parj(171)=0d0
 
C...Add initial e+e- to event record (documentation only).
      ntry=0
  100 ntry=ntry+1
      IF(ntry.GT.100) THEN
        CALL pyerrm(14,'(PYEEVT:) caught in an infinite loop')
        RETURN
      ENDIF
      mstu(24)=0
      nc=0
      IF(mstj(115).GE.2) THEN
        nc=nc+2
        CALL py1ent(nc-1,11,0.5d0*ecm,0d0,0d0)
        k(nc-1,1)=21
        CALL py1ent(nc,-11,0.5d0*ecm,paru(1),0d0)
        k(nc,1)=21
      ENDIF
 
C...Radiative photon (in initial state).
      mk=0
      ecmc=ecm
      IF(mstj(107).GE.1.AND.mstj(116).GE.1) CALL pyradk(ecm,mk,pak,
     &thek,phik,alpk)
      IF(mk.EQ.1) ecmc=sqrt(ecm*(ecm-2d0*pak))
      IF(mstj(115).GE.1.AND.mk.EQ.1) THEN
        nc=nc+1
        CALL py1ent(nc,22,pak,thek,phik)
        k(nc,3)=min(mstj(115)/2,1)
      ENDIF
 
C...Virtual exchange boson (gamma or Z0).
      IF(mstj(115).GE.3) THEN
        nc=nc+1
        kf=22
        IF(mstj(102).EQ.2) kf=23
        mstu10=mstu(10)
        mstu(10)=1
        p(nc,5)=ecmc
        CALL py1ent(nc,kf,ecmc,0d0,0d0)
        k(nc,1)=21
        k(nc,3)=1
        mstu(10)=mstu10
      ENDIF
 
C...Choice of flavour and jet configuration.
      CALL pyxkfl(kfl,ecm,ecmc,kflc)
      IF(kflc.EQ.0) goto 100
      CALL pyxjet(ecmc,njet,cut)
      kfln=21
      IF(njet.EQ.4) CALL pyx4jt(njet,cut,kflc,ecmc,kfln,x1,x2,x4,
     &x12,x14)
      IF(njet.EQ.3) CALL pyx3jt(njet,cut,kflc,ecmc,x1,x3)
      IF(njet.EQ.2) mstj(120)=1
 
C...Fill jet configuration and origin.
      IF(njet.EQ.2.AND.mstj(101).NE.5) CALL py2ent(nc+1,kflc,-kflc,ecmc)
      IF(njet.EQ.2.AND.mstj(101).EQ.5) CALL py2ent(-(nc+1),kflc,-kflc,
     &ecmc)
      IF(njet.EQ.3) CALL py3ent(nc+1,kflc,21,-kflc,ecmc,x1,x3)
      IF(njet.EQ.4.AND.kfln.EQ.21) CALL py4ent(nc+1,kflc,kfln,kfln,
     &-kflc,ecmc,x1,x2,x4,x12,x14)
      IF(njet.EQ.4.AND.kfln.NE.21) CALL py4ent(nc+1,kflc,-kfln,kfln,
     &-kflc,ecmc,x1,x2,x4,x12,x14)
      IF(mstu(24).NE.0) goto 100
      DO 110 ip=nc+1,n
        k(ip,3)=k(ip,3)+min(mstj(115)/2,1)+(mstj(115)/3)*(nc-1)
  110 CONTINUE
 
C...Angular orientation according to matrix element.
      IF(mstj(106).EQ.1) THEN
        CALL pyxdif(nc,njet,kflc,ecmc,chi,the,phi)
        CALL pyrobo(nc+1,n,0d0,chi,0d0,0d0,0d0)
        CALL pyrobo(nc+1,n,the,phi,0d0,0d0,0d0)
      ENDIF
 
C...Rotation and boost from radiative photon.
      IF(mk.EQ.1) THEN
        dbek=-pak/(ecm-pak)
        nmin=nc+1-mstj(115)/3
        CALL pyrobo(nmin,n,0d0,-phik,0d0,0d0,0d0)
        CALL pyrobo(nmin,n,alpk,0d0,dbek*sin(thek),0d0,dbek*cos(thek))
        CALL pyrobo(nmin,n,0d0,phik,0d0,0d0,0d0)
      ENDIF
 
C...Generate parton shower. Rearrange along strings and check.
      IF(mstj(101).EQ.5) THEN
        CALL pyshow(n-1,n,ecmc)
        mstj14=mstj(14)
        IF(mstj(105).EQ.-1) mstj(14)=-1
        IF(mstj(105).GE.0) mstu(28)=0
        CALL pyprep(0)
        mstj(14)=mstj14
        IF(mstj(105).GE.0.AND.mstu(28).NE.0) goto 100
      ENDIF
 
C...Fragmentation/decay generation. Information for PYTABU.
      IF(mstj(105).EQ.1) CALL pyexec
      mstu(161)=kflc
      mstu(162)=-kflc
 
      RETURN
      END