pyp.f

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

 
C*********************************************************************
 
C...PYP
C...Provides various real-valued event related data.
 
      FUNCTION pyp(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...Local array.
      dimension psum(4)
 
C...Set default value. For I = 0 sum of momenta or charges,
C...or invariant mass of system.
      pyp=0d0
      IF(i.LT.0.OR.i.GT.mstu(4).OR.j.LE.0) THEN
      ELSEIF(i.EQ.0.AND.j.LE.4) THEN
        DO 100 i1=1,n
          IF(k(i1,1).GT.0.AND.k(i1,1).LE.10) pyp=pyp+p(i1,j)
  100   CONTINUE
      ELSEIF(i.EQ.0.AND.j.EQ.5) THEN
        DO 120 j1=1,4
          psum(j1)=0d0
          DO 110 i1=1,n
            IF(k(i1,1).GT.0.AND.k(i1,1).LE.10) psum(j1)=psum(j1)+
     &      p(i1,j1)
  110     CONTINUE
  120   CONTINUE
        pyp=sqrt(max(0d0,psum(4)**2-psum(1)**2-psum(2)**2-psum(3)**2))
      ELSEIF(i.EQ.0.AND.j.EQ.6) THEN
        DO 130 i1=1,n
          IF(k(i1,1).GT.0.AND.k(i1,1).LE.10) pyp=pyp+pychge(k(i1,2))/3d0
  130   CONTINUE
      ELSEIF(i.EQ.0) THEN
 
C...Direct readout of P matrix.
      ELSEIF(j.LE.5) THEN
        pyp=p(i,j)
 
C...Charge, total momentum, transverse momentum, transverse mass.
      ELSEIF(j.LE.12) THEN
        IF(j.EQ.6) pyp=pychge(k(i,2))/3d0
        IF(j.EQ.7.OR.j.EQ.8) pyp=p(i,1)**2+p(i,2)**2+p(i,3)**2
        IF(j.EQ.9.OR.j.EQ.10) pyp=p(i,1)**2+p(i,2)**2
        IF(j.EQ.11.OR.j.EQ.12) pyp=p(i,5)**2+p(i,1)**2+p(i,2)**2
        IF(j.EQ.8.OR.j.EQ.10.OR.j.EQ.12) pyp=sqrt(pyp)
 
C...Theta and phi angle in radians or degrees.
      ELSEIF(j.LE.16) THEN
        IF(j.LE.14) pyp=pyangl(p(i,3),sqrt(p(i,1)**2+p(i,2)**2))
        IF(j.GE.15) pyp=pyangl(p(i,1),p(i,2))
        IF(j.EQ.14.OR.j.EQ.16) pyp=pyp*180d0/paru(1)
 
C...True rapidity, rapidity with pion mass, pseudorapidity.
      ELSEIF(j.LE.19) THEN
        pmr=0d0
        IF(j.EQ.17) pmr=p(i,5)
        IF(j.EQ.18) pmr=pymass(211)
        pr=max(1d-20,pmr**2+p(i,1)**2+p(i,2)**2)
        pyp=sign(log(min((sqrt(pr+p(i,3)**2)+abs(p(i,3)))/sqrt(pr),
     &  1d20)),p(i,3))
 
C...Energy and momentum fractions (only to be used in CM frame).
      ELSEIF(j.LE.25) THEN
        IF(j.EQ.20) pyp=2d0*sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)/paru(21)
        IF(j.EQ.21) pyp=2d0*p(i,3)/paru(21)
        IF(j.EQ.22) pyp=2d0*sqrt(p(i,1)**2+p(i,2)**2)/paru(21)
        IF(j.EQ.23) pyp=2d0*p(i,4)/paru(21)
        IF(j.EQ.24) pyp=(p(i,4)+p(i,3))/paru(21)
        IF(j.EQ.25) pyp=(p(i,4)-p(i,3))/paru(21)
      ENDIF
 
      RETURN
      END