This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lusphe.F
1  
2 C********************************************************************* 
3  
4       SUBROUTINE LUSPHE(SPH,APL) 
5  
6 C...Purpose: to perform sphericity tensor analysis to give sphericity, 
7 C...aplanarity and the related event axes. 
8       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) 
9       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
10       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) 
11       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/ 
12       DIMENSION SM(3,3),SV(3,3) 
13  
14 C...Calculate matrix to be diagonalized. 
15       NP=0 
16       DO 110 J1=1,3 
17       DO 100 J2=J1,3 
18       SM(J1,J2)=0. 
19   100 CONTINUE 
20   110 CONTINUE 
21       PS=0. 
22       DO 140 I=1,N 
23       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 140 
24       IF(MSTU(41).GE.2) THEN 
25         KC=LUCOMP(K(I,2)) 
26         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR. 
27      &  KC.EQ.18) GOTO 140 
28         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0) 
29      &  GOTO 140 
30       ENDIF 
31       NP=NP+1 
32       PA=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) 
33       PWT=1. 
34       IF(ABS(PARU(41)-2.).GT.0.001) PWT=MAX(1E-10,PA)**(PARU(41)-2.) 
35       DO 130 J1=1,3 
36       DO 120 J2=J1,3 
37       SM(J1,J2)=SM(J1,J2)+PWT*P(I,J1)*P(I,J2) 
38   120 CONTINUE 
39   130 CONTINUE 
40       PS=PS+PWT*PA**2 
41   140 CONTINUE 
42  
43 C...Very low multiplicities (0 or 1) not considered. 
44       IF(NP.LE.1) THEN 
45         CALL LUERRM(8,'(LUSPHE:) too few particles for analysis') 
46         SPH=-1. 
47         APL=-1. 
48         RETURN 
49       ENDIF 
50       DO 160 J1=1,3 
51       DO 150 J2=J1,3 
52       SM(J1,J2)=SM(J1,J2)/PS 
53   150 CONTINUE 
54   160 CONTINUE 
55  
56 C...Find eigenvalues to matrix (third degree equation). 
57       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2- 
58      &SM(1,3)**2-SM(2,3)**2)/3.-1./9. 
59       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)* 
60      &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27. 
61       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.) 
62       P(N+1,4)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP) 
63       P(N+3,4)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP) 
64       P(N+2,4)=1.-P(N+1,4)-P(N+3,4) 
65       IF(P(N+2,4).LT.1E-5) THEN 
66         CALL LUERRM(8,'(LUSPHE:) all particles back-to-back') 
67         SPH=-1. 
68         APL=-1. 
69         RETURN 
70       ENDIF 
71  
72 C...Find first and last eigenvector by solving equation system. 
73       DO 240 I=1,3,2 
74       DO 180 J1=1,3 
75       SV(J1,J1)=SM(J1,J1)-P(N+I,4) 
76       DO 170 J2=J1+1,3 
77       SV(J1,J2)=SM(J1,J2) 
78       SV(J2,J1)=SM(J1,J2) 
79   170 CONTINUE 
80   180 CONTINUE 
81       SMAX=0. 
82       DO 200 J1=1,3 
83       DO 190 J2=1,3 
84       IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 190 
85       JA=J1 
86       JB=J2 
87       SMAX=ABS(SV(J1,J2)) 
88   190 CONTINUE 
89   200 CONTINUE 
90       SMAX=0. 
91       DO 220 J3=JA+1,JA+2 
92       J1=J3-3*((J3-1)/3) 
93       RL=SV(J1,JB)/SV(JA,JB) 
94       DO 210 J2=1,3 
95       SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2) 
96       IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 210 
97       JC=J1 
98       SMAX=ABS(SV(J1,J2)) 
99   210 CONTINUE 
100   220 CONTINUE 
101       JB1=JB+1-3*(JB/3) 
102       JB2=JB+2-3*((JB+1)/3) 
103       P(N+I,JB1)=-SV(JC,JB2) 
104       P(N+I,JB2)=SV(JC,JB1) 
105       P(N+I,JB)=-(SV(JA,JB1)*P(N+I,JB1)+SV(JA,JB2)*P(N+I,JB2))/ 
106      &SV(JA,JB) 
107       PA=SQRT(P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2) 
108       SGN=(-1.)**INT(RLU(0)+0.5) 
109       DO 230 J=1,3 
110       P(N+I,J)=SGN*P(N+I,J)/PA 
111   230 CONTINUE 
112   240 CONTINUE 
113  
114 C...Middle axis orthogonal to other two. Fill other codes. 
115       SGN=(-1.)**INT(RLU(0)+0.5) 
116       P(N+2,1)=SGN*(P(N+1,2)*P(N+3,3)-P(N+1,3)*P(N+3,2)) 
117       P(N+2,2)=SGN*(P(N+1,3)*P(N+3,1)-P(N+1,1)*P(N+3,3)) 
118       P(N+2,3)=SGN*(P(N+1,1)*P(N+3,2)-P(N+1,2)*P(N+3,1)) 
119       DO 260 I=1,3 
120       K(N+I,1)=31 
121       K(N+I,2)=95 
122       K(N+I,3)=I 
123       K(N+I,4)=0 
124       K(N+I,5)=0 
125       P(N+I,5)=0. 
126       DO 250 J=1,5 
127       V(I,J)=0. 
128   250 CONTINUE 
129   260 CONTINUE 
130  
131 C...Calculate sphericity and aplanarity. Select storing option. 
132       SPH=1.5*(P(N+2,4)+P(N+3,4)) 
133       APL=1.5*P(N+3,4) 
134       MSTU(61)=N+1 
135       MSTU(62)=NP 
136       IF(MSTU(43).LE.1) MSTU(3)=3 
137       IF(MSTU(43).GE.2) N=N+3 
138  
139       RETURN 
140       END