Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luthru.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luthru.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luthru(THR,OBL)
5 
6 C...Purpose: to perform thrust analysis to give thrust, oblateness
7 C...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 tdi(3),tpr(3)
15 
16 C...Take copy of particles that are to be considered in thrust analysis.
17  np=0
18  ps=0.
19  DO 100 i=1,n
20  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 100
21  IF(mstu(41).GE.2) THEN
22  kc=lucomp(k(i,2))
23  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
24  & kc.EQ.18) goto 100
25  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)
26  & goto 100
27  ENDIF
28  IF(n+np+mstu(44)+15.GE.mstu(4)-mstu(32)-5) THEN
29  CALL luerrm(11,'(LUTHRU:) no more memory left in LUJETS')
30  thr=-2.
31  obl=-2.
32  RETURN
33  ENDIF
34  np=np+1
35  k(n+np,1)=23
36  p(n+np,1)=p(i,1)
37  p(n+np,2)=p(i,2)
38  p(n+np,3)=p(i,3)
39  p(n+np,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
40  p(n+np,5)=1.
41  IF(abs(paru(42)-1.).GT.0.001) p(n+np,5)=p(n+np,4)**(paru(42)-1.)
42  ps=ps+p(n+np,4)*p(n+np,5)
43  100 CONTINUE
44 
45 C...Very low multiplicities (0 or 1) not considered.
46  IF(np.LE.1) THEN
47  CALL luerrm(8,'(LUTHRU:) too few particles for analysis')
48  thr=-1.
49  obl=-1.
50  RETURN
51  ENDIF
52 
53 C...Loop over thrust and major. T axis along z direction in latter case.
54  DO 280 ild=1,2
55  IF(ild.EQ.2) THEN
56  k(n+np+1,1)=31
57  phi=ulangl(p(n+np+1,1),p(n+np+1,2))
58  CALL ludbrb(n+1,n+np+1,0.,-phi,0d0,0d0,0d0)
59  the=ulangl(p(n+np+1,3),p(n+np+1,1))
60  CALL ludbrb(n+1,n+np+1,-the,0.,0d0,0d0,0d0)
61  ENDIF
62 
63 C...Find and order particles with highest p (pT for major).
64  DO 110 ilf=n+np+4,n+np+mstu(44)+4
65  110 p(ilf,4)=0.
66  DO 150 i=n+1,n+np
67  IF(ild.EQ.2) p(i,4)=sqrt(p(i,1)**2+p(i,2)**2)
68  DO 120 ilf=n+np+mstu(44)+3,n+np+4,-1
69  IF(p(i,4).LE.p(ilf,4)) goto 130
70  DO 120 j=1,5
71  120 p(ilf+1,j)=p(ilf,j)
72  ilf=n+np+3
73  130 DO 140 j=1,5
74  140 p(ilf+1,j)=p(i,j)
75  150 CONTINUE
76 
77 C...Find and order initial axes with highest thrust (major).
78  DO 160 ilg=n+np+mstu(44)+5,n+np+mstu(44)+15
79  160 p(ilg,4)=0.
80  nc=2**(min(mstu(44),np)-1)
81  DO 220 ilc=1,nc
82  DO 170 j=1,3
83  170 tdi(j)=0.
84  DO 180 ilf=1,min(mstu(44),np)
85  sgn=p(n+np+ilf+3,5)
86  IF(2**ilf*((ilc+2**(ilf-1)-1)/2**ilf).GE.ilc) sgn=-sgn
87  DO 180 j=1,4-ild
88  180 tdi(j)=tdi(j)+sgn*p(n+np+ilf+3,j)
89  tds=tdi(1)**2+tdi(2)**2+tdi(3)**2
90  DO 190 ilg=n+np+mstu(44)+min(ilc,10)+4,n+np+mstu(44)+5,-1
91  IF(tds.LE.p(ilg,4)) goto 200
92  DO 190 j=1,4
93  190 p(ilg+1,j)=p(ilg,j)
94  ilg=n+np+mstu(44)+4
95  200 DO 210 j=1,3
96  210 p(ilg+1,j)=tdi(j)
97  p(ilg+1,4)=tds
98  220 CONTINUE
99 
100 C...Iterate direction of axis until stable maximum.
101  p(n+np+ild,4)=0.
102  ilg=0
103  230 ilg=ilg+1
104  thp=0.
105  240 thps=thp
106  DO 250 j=1,3
107  IF(thp.LE.1e-10) tdi(j)=p(n+np+mstu(44)+4+ilg,j)
108  IF(thp.GT.1e-10) tdi(j)=tpr(j)
109  250 tpr(j)=0.
110  DO 260 i=n+1,n+np
111  sgn=sign(p(i,5),tdi(1)*p(i,1)+tdi(2)*p(i,2)+tdi(3)*p(i,3))
112  DO 260 j=1,4-ild
113  260 tpr(j)=tpr(j)+sgn*p(i,j)
114  thp=sqrt(tpr(1)**2+tpr(2)**2+tpr(3)**2)/ps
115  IF(thp.GE.thps+paru(48)) goto 240
116 
117 C...Save good axis. Try new initial axis until a number of tries agree.
118  IF(thp.LT.p(n+np+ild,4)-paru(48).AND.ilg.LT.min(10,nc)) goto 230
119  IF(thp.GT.p(n+np+ild,4)+paru(48)) THEN
120  iagr=0
121  sgn=(-1.)**int(rlu(0)+0.5)
122  DO 270 j=1,3
123  270 p(n+np+ild,j)=sgn*tpr(j)/(ps*thp)
124  p(n+np+ild,4)=thp
125  p(n+np+ild,5)=0.
126  ENDIF
127  iagr=iagr+1
128  280 IF(iagr.LT.mstu(45).AND.ilg.LT.min(10,nc)) goto 230
129 
130 C...Find minor axis and value by orthogonality.
131  sgn=(-1.)**int(rlu(0)+0.5)
132  p(n+np+3,1)=-sgn*p(n+np+2,2)
133  p(n+np+3,2)=sgn*p(n+np+2,1)
134  p(n+np+3,3)=0.
135  thp=0.
136  DO 290 i=n+1,n+np
137  290 thp=thp+p(i,5)*abs(p(n+np+3,1)*p(i,1)+p(n+np+3,2)*p(i,2))
138  p(n+np+3,4)=thp/ps
139  p(n+np+3,5)=0.
140 
141 C...Fill axis information. Rotate back to original coordinate system.
142  DO 300 ild=1,3
143  k(n+ild,1)=31
144  k(n+ild,2)=96
145  k(n+ild,3)=ild
146  k(n+ild,4)=0
147  k(n+ild,5)=0
148  DO 300 j=1,5
149  p(n+ild,j)=p(n+np+ild,j)
150  300 v(n+ild,j)=0.
151  CALL ludbrb(n+1,n+3,the,phi,0d0,0d0,0d0)
152 
153 C...Select storing option. Calculate thurst and oblateness.
154  mstu(61)=n+1
155  mstu(62)=np
156  IF(mstu(43).LE.1) mstu(3)=3
157  IF(mstu(43).GE.2) n=n+3
158  thr=p(n+1,4)
159  obl=p(n+2,4)-p(n+3,4)
160 
161  RETURN
162  END