ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lusphe.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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(9000,5),p(9000,5),v(9000,5)
9  SAVE /lujets/
10  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11  SAVE /ludat1/
12  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13  SAVE /ludat2/
14  dimension sm(3,3),sv(3,3)
15 
16 C...Calculate matrix to be diagonalized.
17  np=0
18  DO 100 j1=1,3
19  DO 100 j2=j1,3
20  100 sm(j1,j2)=0.
21  ps=0.
22  DO 120 i=1,n
23  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 120
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 120
28  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)
29  & goto 120
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 110 j1=1,3
36  DO 110 j2=j1,3
37  110 sm(j1,j2)=sm(j1,j2)+pwt*p(i,j1)*p(i,j2)
38  ps=ps+pwt*pa**2
39  120 CONTINUE
40 
41 C...Very low multiplicities (0 or 1) not considered.
42  IF(np.LE.1) THEN
43  CALL luerrm(8,'(LUSPHE:) too few particles for analysis')
44  sph=-1.
45  apl=-1.
46  RETURN
47  ENDIF
48  DO 130 j1=1,3
49  DO 130 j2=j1,3
50  130 sm(j1,j2)=sm(j1,j2)/ps
51 
52 C...Find eigenvalues to matrix (third degree equation).
53  sq=(sm(1,1)*sm(2,2)+sm(1,1)*sm(3,3)+sm(2,2)*sm(3,3)-sm(1,2)**2-
54  &sm(1,3)**2-sm(2,3)**2)/3.-1./9.
55  sr=-0.5*(sq+1./9.+sm(1,1)*sm(2,3)**2+sm(2,2)*sm(1,3)**2+sm(3,3)*
56  &sm(1,2)**2-sm(1,1)*sm(2,2)*sm(3,3))+sm(1,2)*sm(1,3)*sm(2,3)+1./27.
57  sp=cos(acos(max(min(sr/sqrt(-sq**3),1.),-1.))/3.)
58  p(n+1,4)=1./3.+sqrt(-sq)*max(2.*sp,sqrt(3.*(1.-sp**2))-sp)
59  p(n+3,4)=1./3.+sqrt(-sq)*min(2.*sp,-sqrt(3.*(1.-sp**2))-sp)
60  p(n+2,4)=1.-p(n+1,4)-p(n+3,4)
61  IF(p(n+2,4).LT.1e-5) THEN
62  CALL luerrm(8,'(LUSPHE:) 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(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(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