hijing/src/lusphe.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUSPHE(SPH,APL)
0005
0006 C...Purpose: to perform sphericity tensor analysis to give sphericity,
0007 C...aplanarity and the related event axes.
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       DIMENSION SM(3,3),SV(3,3)
0015
0016 C...Calculate matrix to be diagonalized.
0017       NP=0
0018       DO 100 J1=1,3
0019       DO 100 J2=J1,3
0020   100 SM(J1,J2)=0.
0021       PS=0.
0022       DO 120 I=1,N
0023       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120
0024       IF(MSTU(41).GE.2) THEN
0025         KC=LUCOMP(K(I,2))
0026         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0027      &  KC.EQ.18) GOTO 120
0028         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
0029      &  GOTO 120
0030       ENDIF
0031       NP=NP+1
0032       PA=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0033       PWT=1.
0034       IF(ABS(PARU(41)-2.).GT.0.001) PWT=MAX(1E-10,PA)**(PARU(41)-2.)
0035       DO 110 J1=1,3
0036       DO 110 J2=J1,3
0037   110 SM(J1,J2)=SM(J1,J2)+PWT*P(I,J1)*P(I,J2)
0038       PS=PS+PWT*PA**2
0039   120 CONTINUE
0040
0041 C...Very low multiplicities (0 or 1) not considered.
0042       IF(NP.LE.1) THEN
0043         CALL LUERRM(8,'(LUSPHE:) too few particles for analysis')
0044         SPH=-1.
0045         APL=-1.
0046         RETURN
0047       ENDIF
0048       DO 130 J1=1,3
0049       DO 130 J2=J1,3
0050   130 SM(J1,J2)=SM(J1,J2)/PS
0051
0052 C...Find eigenvalues to matrix (third degree equation).
0053       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
0054      &SM(1,3)**2-SM(2,3)**2)/3.-1./9.
0055       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
0056      &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27.
0057       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
0058       P(N+1,4)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
0059       P(N+3,4)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP)
0060       P(N+2,4)=1.-P(N+1,4)-P(N+3,4)
0061       IF(P(N+2,4).LT.1E-5) THEN
0062         CALL LUERRM(8,'(LUSPHE:) all particles back-to-back')
0063         SPH=-1.
0064         APL=-1.
0065         RETURN
0066       ENDIF
0067
0068 C...Find first and last eigenvector by solving equation system.
0069       DO 170 I=1,3,2
0070       DO 140 J1=1,3
0071       SV(J1,J1)=SM(J1,J1)-P(N+I,4)
0072       DO 140 J2=J1+1,3
0073       SV(J1,J2)=SM(J1,J2)
0074   140 SV(J2,J1)=SM(J1,J2)
0075       SMAX=0.
0076       DO 150 J1=1,3
0077       DO 150 J2=1,3
0078       IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 150
0079       JA=J1
0080       JB=J2
0081       SMAX=ABS(SV(J1,J2))
0082   150 CONTINUE
0083       SMAX=0.
0084       DO 160 J3=JA+1,JA+2
0085       J1=J3-3*((J3-1)/3)
0086       RL=SV(J1,JB)/SV(JA,JB)
0087       DO 160 J2=1,3
0088       SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
0089       IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 160
0090       JC=J1
0091       SMAX=ABS(SV(J1,J2))
0092   160 CONTINUE
0093       JB1=JB+1-3*(JB/3)
0094       JB2=JB+2-3*((JB+1)/3)
0095       P(N+I,JB1)=-SV(JC,JB2)
0096       P(N+I,JB2)=SV(JC,JB1)
0097       P(N+I,JB)=-(SV(JA,JB1)*P(N+I,JB1)+SV(JA,JB2)*P(N+I,JB2))/
0098      &SV(JA,JB)
0099       PA=SQRT(P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2)
0100       SGN=(-1.)**INT(RLU(0)+0.5)
0101       DO 170 J=1,3
0102   170 P(N+I,J)=SGN*P(N+I,J)/PA
0103
0104 C...Middle axis orthogonal to other two. Fill other codes.
0105       SGN=(-1.)**INT(RLU(0)+0.5)
0106       P(N+2,1)=SGN*(P(N+1,2)*P(N+3,3)-P(N+1,3)*P(N+3,2))
0107       P(N+2,2)=SGN*(P(N+1,3)*P(N+3,1)-P(N+1,1)*P(N+3,3))
0108       P(N+2,3)=SGN*(P(N+1,1)*P(N+3,2)-P(N+1,2)*P(N+3,1))
0109       DO 180 I=1,3
0110       K(N+I,1)=31
0111       K(N+I,2)=95
0112       K(N+I,3)=I
0113       K(N+I,4)=0
0114       K(N+I,5)=0
0115       P(N+I,5)=0.
0116       DO 180 J=1,5
0117   180 V(I,J)=0.
0118
0119 C...Select storing option. Calculate sphericity and aplanarity.
0120       MSTU(61)=N+1
0121       MSTU(62)=NP
0122       IF(MSTU(43).LE.1) MSTU(3)=3
0123       IF(MSTU(43).GE.2) N=N+3
0124       SPH=1.5*(P(N+2,4)+P(N+3,4))
0125       APL=1.5*P(N+3,4)
0126
0127       RETURN
0128       END