]>
Commit | Line | Data |
---|---|---|
e74335a4 | 1 | * $Id$ |
2 | ||
3 | C********************************************************************* | |
4 | ||
5 | SUBROUTINE LUSPHE_HIJING(SPH,APL) | |
6 | ||
7 | C...Purpose: to perform sphericity tensor analysis to give sphericity, | |
8 | C...aplanarity and the related event axes. | |
9 | #include "lujets_hijing.inc" | |
10 | #include "ludat1_hijing.inc" | |
11 | #include "ludat2_hijing.inc" | |
12 | DIMENSION SM(3,3),SV(3,3) | |
13 | ||
14 | C...Calculate matrix to be diagonalized. | |
15 | NP=0 | |
16 | DO 100 J1=1,3 | |
17 | DO 100 J2=J1,3 | |
18 | 100 SM(J1,J2)=0. | |
19 | PS=0. | |
20 | DO 120 I=1,N | |
21 | IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120 | |
22 | IF(MSTU(41).GE.2) THEN | |
23 | KC=LUCOMP_HIJING(K(I,2)) | |
24 | IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR. | |
25 | & KC.EQ.18) GOTO 120 | |
26 | IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2)) | |
27 | $ .EQ.0)GOTO 120 | |
28 | ENDIF | |
29 | NP=NP+1 | |
30 | PA=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) | |
31 | PWT=1. | |
32 | IF(ABS(PARU(41)-2.).GT.0.001) PWT=MAX(1E-10,PA)**(PARU(41)-2.) | |
33 | DO 110 J1=1,3 | |
34 | DO 110 J2=J1,3 | |
35 | 110 SM(J1,J2)=SM(J1,J2)+PWT*P(I,J1)*P(I,J2) | |
36 | PS=PS+PWT*PA**2 | |
37 | 120 CONTINUE | |
38 | ||
39 | C...Very low multiplicities (0 or 1) not considered. | |
40 | IF(NP.LE.1) THEN | |
41 | CALL LUERRM_HIJING(8 | |
42 | $ ,'(LUSPHE_HIJING:) too few particles for analysis') | |
43 | SPH=-1. | |
44 | APL=-1. | |
45 | RETURN | |
46 | ENDIF | |
47 | DO 130 J1=1,3 | |
48 | DO 130 J2=J1,3 | |
49 | 130 SM(J1,J2)=SM(J1,J2)/PS | |
50 | ||
51 | C...Find eigenvalues to matrix (third degree equation). | |
52 | SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2- | |
53 | &SM(1,3)**2-SM(2,3)**2)/3.-1./9. | |
54 | SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)* | |
55 | &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27. | |
56 | SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.) | |
57 | P(N+1,4)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP) | |
58 | P(N+3,4)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP) | |
59 | P(N+2,4)=1.-P(N+1,4)-P(N+3,4) | |
60 | IF(P(N+2,4).LT.1E-5) THEN | |
61 | CALL LUERRM_HIJING(8 | |
62 | $ ,'(LUSPHE_HIJING:) all particles back-to-back') | |
63 | SPH=-1. | |
64 | APL=-1. | |
65 | RETURN | |
66 | ENDIF | |
67 | ||
68 | C...Find first and last eigenvector by solving equation system. | |
69 | DO 170 I=1,3,2 | |
70 | DO 140 J1=1,3 | |
71 | SV(J1,J1)=SM(J1,J1)-P(N+I,4) | |
72 | DO 140 J2=J1+1,3 | |
73 | SV(J1,J2)=SM(J1,J2) | |
74 | 140 SV(J2,J1)=SM(J1,J2) | |
75 | SMAX=0. | |
76 | DO 150 J1=1,3 | |
77 | DO 150 J2=1,3 | |
78 | IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 150 | |
79 | JA=J1 | |
80 | JB=J2 | |
81 | SMAX=ABS(SV(J1,J2)) | |
82 | 150 CONTINUE | |
83 | SMAX=0. | |
84 | DO 160 J3=JA+1,JA+2 | |
85 | J1=J3-3*((J3-1)/3) | |
86 | RL=SV(J1,JB)/SV(JA,JB) | |
87 | DO 160 J2=1,3 | |
88 | SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2) | |
89 | IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 160 | |
90 | JC=J1 | |
91 | SMAX=ABS(SV(J1,J2)) | |
92 | 160 CONTINUE | |
93 | JB1=JB+1-3*(JB/3) | |
94 | JB2=JB+2-3*((JB+1)/3) | |
95 | P(N+I,JB1)=-SV(JC,JB2) | |
96 | P(N+I,JB2)=SV(JC,JB1) | |
97 | P(N+I,JB)=-(SV(JA,JB1)*P(N+I,JB1)+SV(JA,JB2)*P(N+I,JB2))/ | |
98 | &SV(JA,JB) | |
99 | PA=SQRT(P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2) | |
100 | SGN=(-1.)**INT(RLU_HIJING(0)+0.5) | |
101 | DO 170 J=1,3 | |
102 | 170 P(N+I,J)=SGN*P(N+I,J)/PA | |
103 | ||
104 | C...Middle axis orthogonal to other two. Fill other codes. | |
105 | SGN=(-1.)**INT(RLU_HIJING(0)+0.5) | |
106 | P(N+2,1)=SGN*(P(N+1,2)*P(N+3,3)-P(N+1,3)*P(N+3,2)) | |
107 | P(N+2,2)=SGN*(P(N+1,3)*P(N+3,1)-P(N+1,1)*P(N+3,3)) | |
108 | P(N+2,3)=SGN*(P(N+1,1)*P(N+3,2)-P(N+1,2)*P(N+3,1)) | |
109 | DO 180 I=1,3 | |
110 | K(N+I,1)=31 | |
111 | K(N+I,2)=95 | |
112 | K(N+I,3)=I | |
113 | K(N+I,4)=0 | |
114 | K(N+I,5)=0 | |
115 | P(N+I,5)=0. | |
116 | DO 180 J=1,5 | |
117 | 180 V(I,J)=0. | |
118 | ||
119 | C...Select storing option. Calculate sphericity and aplanarity. | |
120 | MSTU(61)=N+1 | |
121 | MSTU(62)=NP | |
122 | IF(MSTU(43).LE.1) MSTU(3)=3 | |
123 | IF(MSTU(43).GE.2) N=N+3 | |
124 | SPH=1.5*(P(N+2,4)+P(N+3,4)) | |
125 | APL=1.5*P(N+3,4) | |
126 | ||
127 | RETURN | |
128 | END |