pyonia.f

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

 
C*********************************************************************
 
C...PYONIA
C...Generates Upsilon and toponium decays into three gluons
C...or two gluons and a photon.
 
      SUBROUTINE pyonia(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...Printout. Check input parameters.
      IF(mstu(12).NE.12345) CALL pylist(0)
      IF(kfl.LT.0.OR.kfl.GT.8) THEN
        CALL pyerrm(16,'(PYONIA:) called with unknown flavour code')
        IF(mstu(21).GE.1) RETURN
      ENDIF
      IF(ecm.LT.parj(127)+2.02d0*parf(101)) THEN
        CALL pyerrm(16,'(PYONIA:) called with too small CM energy')
        IF(mstu(21).GE.1) RETURN
      ENDIF
 
C...Initial e+e- and onium state (optional).
      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
      kflc=iabs(kfl)
      IF(mstj(115).GE.3.AND.kflc.GE.5) THEN
        nc=nc+1
        kf=110*kflc+3
        mstu10=mstu(10)
        mstu(10)=1
        p(nc,5)=ecm
        CALL py1ent(nc,kf,ecm,0d0,0d0)
        k(nc,1)=21
        k(nc,3)=1
        mstu(10)=mstu10
      ENDIF
 
C...Choose x1 and x2 according to matrix element.
      ntry=0
  100 x1=pyr(0)
      x2=pyr(0)
      x3=2d0-x1-x2
      IF(x3.GE.1d0.OR.((1d0-x1)/(x2*x3))**2+((1d0-x2)/(x1*x3))**2+
     &((1d0-x3)/(x1*x2))**2.LE.2d0*pyr(0)) goto 100
      ntry=ntry+1
      njet=3
      IF(mstj(101).LE.4) CALL py3ent(nc+1,21,21,21,ecm,x1,x3)
      IF(mstj(101).GE.5) CALL py3ent(-(nc+1),21,21,21,ecm,x1,x3)
 
C...Photon-gluon-gluon events. Small system modifications. Jet origin.
      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)
      qf=0d0
      IF(kflc.NE.0) qf=kchg(kflc,1)/3d0
      rgam=7.2d0*qf**2*paru(101)/pyalps(ecm**2)
      mk=0
      ecmc=ecm
      IF(pyr(0).GT.rgam/(1d0+rgam)) THEN
        IF(1d0-max(x1,x2,x3).LE.max((parj(126)/ecm)**2,parj(125)))
     &  njet=2
        IF(njet.EQ.2.AND.mstj(101).LE.4) CALL py2ent(nc+1,21,21,ecm)
        IF(njet.EQ.2.AND.mstj(101).GE.5) CALL py2ent(-(nc+1),21,21,ecm)
      ELSE
        mk=1
        ecmc=sqrt(1d0-x1)*ecm
        IF(ecmc.LT.2d0*parj(127)) goto 100
        k(nc+1,1)=1
        k(nc+1,2)=22
        k(nc+1,4)=0
        k(nc+1,5)=0
        IF(mstj(101).GE.5) k(nc+2,4)=mstu(5)*(nc+3)
        IF(mstj(101).GE.5) k(nc+2,5)=mstu(5)*(nc+3)
        IF(mstj(101).GE.5) k(nc+3,4)=mstu(5)*(nc+2)
        IF(mstj(101).GE.5) k(nc+3,5)=mstu(5)*(nc+2)
        njet=2
        IF(ecmc.LT.4d0*parj(127)) THEN
          mstu10=mstu(10)
          mstu(10)=1
          p(nc+2,5)=ecmc
          CALL py1ent(nc+2,83,0.5d0*(x2+x3)*ecm,paru(1),0d0)
          mstu(10)=mstu10
          njet=0
        ENDIF
      ENDIF
      DO 110 ip=nc+1,n
        k(ip,3)=k(ip,3)+(mstj(115)/2)+(kflc/5)*(mstj(115)/3)*(nc-1)
  110 CONTINUE
 
C...Differential cross-sections. Upper limit for cross-section.
      IF(mstj(106).EQ.1) THEN
        sq2=sqrt(2d0)
        hf1=1d0-parj(131)*parj(132)
        hf3=parj(133)**2
        ct13=(x1*x3-2d0*x1-2d0*x3+2d0)/(x1*x3)
        st13=sqrt(1d0-ct13**2)
        sigl=0.5d0*x3**2*((1d0-x2)**2+(1d0-x3)**2)*st13**2
        sigu=(x1*(1d0-x1))**2+(x2*(1d0-x2))**2+(x3*(1d0-x3))**2-sigl
        sigt=0.5d0*sigl
        sigi=(sigl*ct13/st13+0.5d0*x1*x3*(1d0-x2)**2*st13)/sq2
        sigmax=(2d0*hf1+hf3)*abs(sigu)+2d0*(hf1+hf3)*abs(sigl)+2d0*(hf1+
     &  2d0*hf3)*abs(sigt)+2d0*sq2*(hf1+2d0*hf3)*abs(sigi)
 
C...Angular orientation of event.
  120   chi=paru(2)*pyr(0)
        cthe=2d0*pyr(0)-1d0
        phi=paru(2)*pyr(0)
        cchi=cos(chi)
        schi=sin(chi)
        c2chi=cos(2d0*chi)
        s2chi=sin(2d0*chi)
        the=acos(cthe)
        sthe=sin(the)
        c2phi=cos(2d0*(phi-parj(134)))
        s2phi=sin(2d0*(phi-parj(134)))
        sig=((1d0+cthe**2)*hf1+sthe**2*c2phi*hf3)*sigu+2d0*(sthe**2*hf1-
     &  sthe**2*c2phi*hf3)*sigl+2d0*(sthe**2*c2chi*hf1+((1d0+cthe**2)*
     &  c2chi*c2phi-2d0*cthe*s2chi*s2phi)*hf3)*sigt-
     &  2d0*sq2*(2d0*sthe*cthe*cchi*hf1-2d0*sthe*
     &  (cthe*cchi*c2phi-schi*s2phi)*hf3)*sigi
        IF(sig.LT.sigmax*pyr(0)) goto 120
        CALL pyrobo(nc+1,n,0d0,chi,0d0,0d0,0d0)
        CALL pyrobo(nc+1,n,the,phi,0d0,0d0,0d0)
      ENDIF
 
C...Generate parton shower. Rearrange along strings and check.
      IF(mstj(101).GE.5.AND.njet.GE.2) THEN
        CALL pyshow(nc+mk+1,-njet,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...Generate fragmentation. Information for PYTABU:
      IF(mstj(105).EQ.1) CALL pyexec
      mstu(161)=110*kflc+3
      mstu(162)=0
 
      RETURN
      END