pythiaeRHIC/pythia/pysphe.f

0001
0002 C*********************************************************************
0003
0004 C...PYSPHE
0005 C...Performs sphericity tensor analysis to give sphericity,
0006 C...aplanarity and the related event axes.
0007
0008       SUBROUTINE PYSPHE(SPH,APL)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Parameter statement to help give large particle numbers.
0015       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0016      &KEXCIT=4000000,KDIMEN=5000000)
0017 C...Commonblocks.
0018       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0022 C...Local arrays.
0023       DIMENSION SM(3,3),SV(3,3)
0024
0025 C...Calculate matrix to be diagonalized.
0026       NP=0
0027       DO 110 J1=1,3
0028         DO 100 J2=J1,3
0029           SM(J1,J2)=0D0
0030   100   CONTINUE
0031   110 CONTINUE
0032       PS=0D0
0033       DO 140 I=1,N
0034         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 140
0035         IF(MSTU(41).GE.2) THEN
0036           KC=PYCOMP(K(I,2))
0037           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0038      &    KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0039      &    K(I,2).EQ.KSUSY1+39) GOTO 140
0040           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0)
0041      &    GOTO 140
0042         ENDIF
0043         NP=NP+1
0044         PA=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0045         PWT=1D0
0046         IF(ABS(PARU(41)-2D0).GT.0.001D0) PWT=
0047      &  MAX(1D-10,PA)**(PARU(41)-2D0)
0048         DO 130 J1=1,3
0049           DO 120 J2=J1,3
0050             SM(J1,J2)=SM(J1,J2)+PWT*P(I,J1)*P(I,J2)
0051   120     CONTINUE
0052   130   CONTINUE
0053         PS=PS+PWT*PA**2
0054   140 CONTINUE
0055
0056 C...Very low multiplicities (0 or 1) not considered.
0057       IF(NP.LE.1) THEN
0058         CALL PYERRM(8,'(PYSPHE:) too few particles for analysis')
0059         SPH=-1D0
0060         APL=-1D0
0061         RETURN
0062       ENDIF
0063       DO 160 J1=1,3
0064         DO 150 J2=J1,3
0065           SM(J1,J2)=SM(J1,J2)/PS
0066   150   CONTINUE
0067   160 CONTINUE
0068
0069 C...Find eigenvalues to matrix (third degree equation).
0070       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-
0071      &SM(1,2)**2-SM(1,3)**2-SM(2,3)**2)/3D0-1D0/9D0
0072       SR=-0.5D0*(SQ+1D0/9D0+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+
0073      &SM(3,3)*SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+
0074      &SM(1,2)*SM(1,3)*SM(2,3)+1D0/27D0
0075       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1D0),-1D0))/3D0)
0076       P(N+1,4)=1D0/3D0+SQRT(-SQ)*MAX(2D0*SP,SQRT(3D0*(1D0-SP**2))-SP)
0077       P(N+3,4)=1D0/3D0+SQRT(-SQ)*MIN(2D0*SP,-SQRT(3D0*(1D0-SP**2))-SP)
0078       P(N+2,4)=1D0-P(N+1,4)-P(N+3,4)
0079       IF(P(N+2,4).LT.1D-5) THEN
0080         CALL PYERRM(8,'(PYSPHE:) all particles back-to-back')
0081         SPH=-1D0
0082         APL=-1D0
0083         RETURN
0084       ENDIF
0085
0086 C...Find first and last eigenvector by solving equation system.
0087       DO 240 I=1,3,2
0088         DO 180 J1=1,3
0089           SV(J1,J1)=SM(J1,J1)-P(N+I,4)
0090           DO 170 J2=J1+1,3
0091             SV(J1,J2)=SM(J1,J2)
0092             SV(J2,J1)=SM(J1,J2)
0093   170     CONTINUE
0094   180   CONTINUE
0095         SMAX=0D0
0096         DO 200 J1=1,3
0097           DO 190 J2=1,3
0098             IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 190
0099             JA=J1
0100             JB=J2
0101             SMAX=ABS(SV(J1,J2))
0102   190     CONTINUE
0103   200   CONTINUE
0104         SMAX=0D0
0105         DO 220 J3=JA+1,JA+2
0106           J1=J3-3*((J3-1)/3)
0107           RL=SV(J1,JB)/SV(JA,JB)
0108           DO 210 J2=1,3
0109             SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
0110             IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 210
0111             JC=J1
0112             SMAX=ABS(SV(J1,J2))
0113   210     CONTINUE
0114   220   CONTINUE
0115         JB1=JB+1-3*(JB/3)
0116         JB2=JB+2-3*((JB+1)/3)
0117         P(N+I,JB1)=-SV(JC,JB2)
0118         P(N+I,JB2)=SV(JC,JB1)
0119         P(N+I,JB)=-(SV(JA,JB1)*P(N+I,JB1)+SV(JA,JB2)*P(N+I,JB2))/
0120      &  SV(JA,JB)
0121         PA=SQRT(P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2)
0122         SGN=(-1D0)**INT(PYR(0)+0.5D0)
0123         DO 230 J=1,3
0124           P(N+I,J)=SGN*P(N+I,J)/PA
0125   230   CONTINUE
0126   240 CONTINUE
0127
0128 C...Middle axis orthogonal to other two. Fill other codes.
0129       SGN=(-1D0)**INT(PYR(0)+0.5D0)
0130       P(N+2,1)=SGN*(P(N+1,2)*P(N+3,3)-P(N+1,3)*P(N+3,2))
0131       P(N+2,2)=SGN*(P(N+1,3)*P(N+3,1)-P(N+1,1)*P(N+3,3))
0132       P(N+2,3)=SGN*(P(N+1,1)*P(N+3,2)-P(N+1,2)*P(N+3,1))
0133       DO 260 I=1,3
0134         K(N+I,1)=31
0135         K(N+I,2)=95
0136         K(N+I,3)=I
0137         K(N+I,4)=0
0138         K(N+I,5)=0
0139         P(N+I,5)=0D0
0140         DO 250 J=1,5
0141           V(I,J)=0D0
0142   250   CONTINUE
0143   260 CONTINUE
0144
0145 C...Calculate sphericity and aplanarity. Select storing option.
0146       SPH=1.5D0*(P(N+2,4)+P(N+3,4))
0147       APL=1.5D0*P(N+3,4)
0148       MSTU(61)=N+1
0149       MSTU(62)=NP
0150       IF(MSTU(43).LE.1) MSTU(3)=3
0151       IF(MSTU(43).GE.2) N=N+3
0152
0153       RETURN
0154       END