pythiaeRHIC/pythia/pyp.f

0001
0002 C*********************************************************************
0003
0004 C...PYP
0005 C...Provides various real-valued event related data.
0006
0007       FUNCTION PYP(I,J)
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0017       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0018 C...Local array.
0019       DIMENSION PSUM(4)
0020
0021 C...Set default value. For I = 0 sum of momenta or charges,
0022 C...or invariant mass of system.
0023       PYP=0D0
0024       IF(I.LT.0.OR.I.GT.MSTU(4).OR.J.LE.0) THEN
0025       ELSEIF(I.EQ.0.AND.J.LE.4) THEN
0026         DO 100 I1=1,N
0027           IF(K(I1,1).GT.0.AND.K(I1,1).LE.10) PYP=PYP+P(I1,J)
0028   100   CONTINUE
0029       ELSEIF(I.EQ.0.AND.J.EQ.5) THEN
0030         DO 120 J1=1,4
0031           PSUM(J1)=0D0
0032           DO 110 I1=1,N
0033             IF(K(I1,1).GT.0.AND.K(I1,1).LE.10) PSUM(J1)=PSUM(J1)+
0034      &      P(I1,J1)
0035   110     CONTINUE
0036   120   CONTINUE
0037         PYP=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2))
0038       ELSEIF(I.EQ.0.AND.J.EQ.6) THEN
0039         DO 130 I1=1,N
0040           IF(K(I1,1).GT.0.AND.K(I1,1).LE.10) PYP=PYP+PYCHGE(K(I1,2))/3D0
0041   130   CONTINUE
0042       ELSEIF(I.EQ.0) THEN
0043
0044 C...Direct readout of P matrix.
0045       ELSEIF(J.LE.5) THEN
0046         PYP=P(I,J)
0047
0048 C...Charge, total momentum, transverse momentum, transverse mass.
0049       ELSEIF(J.LE.12) THEN
0050         IF(J.EQ.6) PYP=PYCHGE(K(I,2))/3D0
0051         IF(J.EQ.7.OR.J.EQ.8) PYP=P(I,1)**2+P(I,2)**2+P(I,3)**2
0052         IF(J.EQ.9.OR.J.EQ.10) PYP=P(I,1)**2+P(I,2)**2
0053         IF(J.EQ.11.OR.J.EQ.12) PYP=P(I,5)**2+P(I,1)**2+P(I,2)**2
0054         IF(J.EQ.8.OR.J.EQ.10.OR.J.EQ.12) PYP=SQRT(PYP)
0055
0056 C...Theta and phi angle in radians or degrees.
0057       ELSEIF(J.LE.16) THEN
0058         IF(J.LE.14) PYP=PYANGL(P(I,3),SQRT(P(I,1)**2+P(I,2)**2))
0059         IF(J.GE.15) PYP=PYANGL(P(I,1),P(I,2))
0060         IF(J.EQ.14.OR.J.EQ.16) PYP=PYP*180D0/PARU(1)
0061
0062 C...True rapidity, rapidity with pion mass, pseudorapidity.
0063       ELSEIF(J.LE.19) THEN
0064         PMR=0D0
0065         IF(J.EQ.17) PMR=P(I,5)
0066         IF(J.EQ.18) PMR=PYMASS(211)
0067         PR=MAX(1D-20,PMR**2+P(I,1)**2+P(I,2)**2)
0068         PYP=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/SQRT(PR),
0069      &  1D20)),P(I,3))
0070
0071 C...Energy and momentum fractions (only to be used in CM frame).
0072       ELSEIF(J.LE.25) THEN
0073         IF(J.EQ.20) PYP=2D0*SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)/PARU(21)
0074         IF(J.EQ.21) PYP=2D0*P(I,3)/PARU(21)
0075         IF(J.EQ.22) PYP=2D0*SQRT(P(I,1)**2+P(I,2)**2)/PARU(21)
0076         IF(J.EQ.23) PYP=2D0*P(I,4)/PARU(21)
0077         IF(J.EQ.24) PYP=(P(I,4)+P(I,3))/PARU(21)
0078         IF(J.EQ.25) PYP=(P(I,4)-P(I,3))/PARU(21)
0079       ENDIF
0080
0081       RETURN
0082       END