pyk.f

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

 
C*********************************************************************
 
C...PYK
C...Provides various integer-valued event related data.
 
      FUNCTION pyk(I,J)
 
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...Default value. For I=0 number of entries, number of stable entries
C...or 3 times total charge.
      pyk=0
      IF(i.LT.0.OR.i.GT.mstu(4).OR.j.LE.0) THEN
      ELSEIF(i.EQ.0.AND.j.EQ.1) THEN
        pyk=n
      ELSEIF(i.EQ.0.AND.(j.EQ.2.OR.j.EQ.6)) THEN
        DO 100 i1=1,n
          IF(j.EQ.2.AND.k(i1,1).GE.1.AND.k(i1,1).LE.10) pyk=pyk+1
          IF(j.EQ.6.AND.k(i1,1).GE.1.AND.k(i1,1).LE.10) pyk=pyk+
     &    pychge(k(i1,2))
  100   CONTINUE
      ELSEIF(i.EQ.0) THEN
 
C...For I > 0 direct readout of K matrix or charge.
      ELSEIF(j.LE.5) THEN
        pyk=k(i,j)
      ELSEIF(j.EQ.6) THEN
        pyk=pychge(k(i,2))
 
C...Status (existing/fragmented/decayed), parton/hadron separation.
      ELSEIF(j.LE.8) THEN
        IF(k(i,1).GE.1.AND.k(i,1).LE.10) pyk=1
        IF(j.EQ.8) pyk=pyk*k(i,2)
      ELSEIF(j.LE.12) THEN
        kfa=iabs(k(i,2))
        kc=pycomp(kfa)
        kq=0
        IF(kc.NE.0) kq=kchg(kc,2)
        IF(j.EQ.9.AND.kc.NE.0.AND.kq.NE.0) pyk=k(i,2)
        IF(j.EQ.10.AND.kc.NE.0.AND.kq.EQ.0) pyk=k(i,2)
        IF(j.EQ.11) pyk=kc
        IF(j.EQ.12) pyk=kq*isign(1,k(i,2))
 
C...Heaviest flavour in hadron/diquark.
      ELSEIF(j.EQ.13) THEN
        kfa=iabs(k(i,2))
        pyk=mod(kfa/100,10)*(-1)**mod(kfa/100,10)
        IF(kfa.LT.10) pyk=kfa
        IF(mod(kfa/1000,10).NE.0) pyk=mod(kfa/1000,10)
        pyk=pyk*isign(1,k(i,2))
 
C...Particle history: generation, ancestor, rank.
      ELSEIF(j.LE.15) THEN
        i2=i
        i1=i
  110   pyk=pyk+1
        i2=i1
        i1=k(i1,3)
        IF(i1.GT.0) THEN
          IF(k(i1,1).GT.0.AND.k(i1,1).LE.20) goto 110
        ENDIF
        IF(j.EQ.15) pyk=i2
      ELSEIF(j.EQ.16) THEN
        kfa=iabs(k(i,2))
        IF(k(i,1).LE.20.AND.((kfa.GE.11.AND.kfa.LE.20).OR.kfa.EQ.22.OR.
     &  (kfa.GT.100.AND.mod(kfa/10,10).NE.0))) THEN
          i1=i
  120     i2=i1
          i1=k(i1,3)
          IF(i1.GT.0) THEN
            kfam=iabs(k(i1,2))
            ilp=1
            IF(kfam.NE.0.AND.kfam.LE.10) ilp=0
            IF(kfam.EQ.21.OR.kfam.EQ.91.OR.kfam.EQ.92.OR.kfam.EQ.93)
     &      ilp=0
            IF(kfam.GT.100.AND.mod(kfam/10,10).EQ.0) ilp=0
            IF(ilp.EQ.1) goto 120
          ENDIF
          IF(k(i1,1).EQ.12) THEN
            DO 130 i3=i1+1,i2
              IF(k(i3,3).EQ.k(i2,3).AND.k(i3,2).NE.91.AND.k(i3,2).NE.92
     &        .AND.k(i3,2).NE.93) pyk=pyk+1
  130       CONTINUE
          ELSE
            i3=i2
  140       pyk=pyk+1
            i3=i3+1
            IF(i3.LT.n.AND.k(i3,3).EQ.k(i2,3)) goto 140
          ENDIF
        ENDIF
 
C...Particle coming from collapsing jet system or not.
      ELSEIF(j.EQ.17) THEN
        i1=i
  150   pyk=pyk+1
        i3=i1
        i1=k(i1,3)
        i0=max(1,i1)
        kc=pycomp(k(i0,2))
        IF(i1.EQ.0.OR.k(i0,1).LE.0.OR.k(i0,1).GT.20.OR.kc.EQ.0) THEN
          IF(pyk.EQ.1) pyk=-1
          IF(pyk.GT.1) pyk=0
          RETURN
        ENDIF
        IF(kchg(kc,2).EQ.0) goto 150
        IF(k(i1,1).NE.12) pyk=0
        IF(k(i1,1).NE.12) RETURN
        i2=i1
  160   i2=i2+1
        IF(i2.LT.n.AND.k(i2,1).NE.11) goto 160
        k3m=k(i3-1,3)
        IF(k3m.GE.i1.AND.k3m.LE.i2) pyk=0
        k3p=k(i3+1,3)
        IF(i3.LT.n.AND.k3p.GE.i1.AND.k3p.LE.i2) pyk=0
 
C...Number of decay products. Colour flow.
      ELSEIF(j.EQ.18) THEN
        IF(k(i,1).EQ.11.OR.k(i,1).EQ.12) pyk=max(0,k(i,5)-k(i,4)+1)
        IF(k(i,4).EQ.0.OR.k(i,5).EQ.0) pyk=0
      ELSEIF(j.LE.22) THEN
        IF(k(i,1).NE.3.AND.k(i,1).NE.13.AND.k(i,1).NE.14) RETURN
        IF(j.EQ.19) pyk=mod(k(i,4)/mstu(5),mstu(5))
        IF(j.EQ.20) pyk=mod(k(i,5)/mstu(5),mstu(5))
        IF(j.EQ.21) pyk=mod(k(i,4),mstu(5))
        IF(j.EQ.22) pyk=mod(k(i,5),mstu(5))
      ELSE
      ENDIF
 
      RETURN
      END