Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luclus.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luclus.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luclus(NJET)
5 
6 C...Purpose: to subdivide the particle content of an event into
7 C...jets/clusters.
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 ps(5)
15  SAVE nsav,np,ps,pss,rinit,npre,nrem
16 
17 C...Functions: distance measure in pT or (pseudo)mass.
18  r2t(i1,i2)=(p(i1,5)*p(i2,5)-p(i1,1)*p(i2,1)-p(i1,2)*p(i2,2)-
19  &p(i1,3)*p(i2,3))*2.*p(i1,5)*p(i2,5)/(0.0001+p(i1,5)+p(i2,5))**2
20  r2m(i1,i2)=2.*p(i1,4)*p(i2,4)*(1.-(p(i1,1)*p(i2,1)+p(i1,2)*
21  &p(i2,2)+p(i1,3)*p(i2,3))/(p(i1,5)*p(i2,5)))
22 
23 C...If first time, reset. If reentering, skip preliminaries.
24  IF(mstu(48).LE.0) THEN
25  np=0
26  DO 100 j=1,5
27  100 ps(j)=0.
28  pss=0.
29  ELSE
30  njet=nsav
31  IF(mstu(43).GE.2) n=n-njet
32  DO 110 i=n+1,n+njet
33  110 p(i,5)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
34  IF(mstu(46).LE.3) r2acc=paru(44)**2
35  IF(mstu(46).GE.4) r2acc=paru(45)*ps(5)**2
36  nloop=0
37  goto 290
38  ENDIF
39 
40 C...Find which particles are to be considered in cluster search.
41  DO 140 i=1,n
42  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 140
43  IF(mstu(41).GE.2) THEN
44  kc=lucomp(k(i,2))
45  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
46  & kc.EQ.18) goto 140
47  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)
48  & goto 140
49  ENDIF
50  IF(n+2*np.GE.mstu(4)-mstu(32)-5) THEN
51  CALL luerrm(11,'(LUCLUS:) no more memory left in LUJETS')
52  njet=-1
53  RETURN
54  ENDIF
55 
56 C...Take copy of these particles, with space left for jets later on.
57  np=np+1
58  k(n+np,3)=i
59  DO 120 j=1,5
60  120 p(n+np,j)=p(i,j)
61  IF(mstu(42).EQ.0) p(n+np,5)=0.
62  IF(mstu(42).EQ.1.AND.k(i,2).NE.22) p(n+np,5)=pmas(101,1)
63  p(n+np,4)=sqrt(p(n+np,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
64  p(n+np,5)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
65  DO 130 j=1,4
66  130 ps(j)=ps(j)+p(n+np,j)
67  pss=pss+p(n+np,5)
68  140 CONTINUE
69  DO 150 i=n+1,n+np
70  k(i+np,3)=k(i,3)
71  DO 150 j=1,5
72  150 p(i+np,j)=p(i,j)
73  ps(5)=sqrt(max(0.,ps(4)**2-ps(1)**2-ps(2)**2-ps(3)**2))
74 
75 C...Very low multiplicities not considered.
76  IF(np.LT.mstu(47)) THEN
77  CALL luerrm(8,'(LUCLUS:) too few particles for analysis')
78  njet=-1
79  RETURN
80  ENDIF
81 
82 C...Find precluster configuration. If too few jets, make harder cuts.
83  nloop=0
84  IF(mstu(46).LE.3) r2acc=paru(44)**2
85  IF(mstu(46).GE.4) r2acc=paru(45)*ps(5)**2
86  rinit=1.25*paru(43)
87  IF(np.LE.mstu(47)+2) rinit=0.
88  160 rinit=0.8*rinit
89  npre=0
90  nrem=np
91  DO 170 i=n+np+1,n+2*np
92  170 k(i,4)=0
93 
94 C...Sum up small momentum region. Jet if enough absolute momentum.
95  IF(mstu(46).LE.2) THEN
96  DO 180 j=1,4
97  180 p(n+1,j)=0.
98  DO 200 i=n+np+1,n+2*np
99  IF(p(i,5).GT.2.*rinit) goto 200
100  nrem=nrem-1
101  k(i,4)=1
102  DO 190 j=1,4
103  190 p(n+1,j)=p(n+1,j)+p(i,j)
104  200 CONTINUE
105  p(n+1,5)=sqrt(p(n+1,1)**2+p(n+1,2)**2+p(n+1,3)**2)
106  IF(p(n+1,5).GT.2.*rinit) npre=1
107  IF(rinit.GE.0.2*paru(43).AND.npre+nrem.LT.mstu(47)) goto 160
108  ENDIF
109 
110 C...Find fastest remaining particle.
111  210 npre=npre+1
112  pmax=0.
113  DO 220 i=n+np+1,n+2*np
114  IF(k(i,4).NE.0.OR.p(i,5).LE.pmax) goto 220
115  imax=i
116  pmax=p(i,5)
117  220 CONTINUE
118  DO 230 j=1,5
119  230 p(n+npre,j)=p(imax,j)
120  nrem=nrem-1
121  k(imax,4)=npre
122 
123 C...Sum up precluster around it according to pT separation.
124  IF(mstu(46).LE.2) THEN
125  DO 250 i=n+np+1,n+2*np
126  IF(k(i,4).NE.0) goto 250
127  r2=r2t(i,imax)
128  IF(r2.GT.rinit**2) goto 250
129  nrem=nrem-1
130  k(i,4)=npre
131  DO 240 j=1,4
132  240 p(n+npre,j)=p(n+npre,j)+p(i,j)
133  250 CONTINUE
134  p(n+npre,5)=sqrt(p(n+npre,1)**2+p(n+npre,2)**2+p(n+npre,3)**2)
135 
136 C...Sum up precluster around it according to mass separation.
137  ELSE
138  260 imin=0
139  r2min=rinit**2
140  DO 270 i=n+np+1,n+2*np
141  IF(k(i,4).NE.0) goto 270
142  r2=r2m(i,n+npre)
143  IF(r2.GE.r2min) goto 270
144  imin=i
145  r2min=r2
146  270 CONTINUE
147  IF(imin.NE.0) THEN
148  DO 280 j=1,4
149  280 p(n+npre,j)=p(n+npre,j)+p(imin,j)
150  p(n+npre,5)=sqrt(p(n+npre,1)**2+p(n+npre,2)**2+p(n+npre,3)**2)
151  nrem=nrem-1
152  k(imin,4)=npre
153  goto 260
154  ENDIF
155  ENDIF
156 
157 C...Check if more preclusters to be found. Start over if too few.
158  IF(rinit.GE.0.2*paru(43).AND.npre+nrem.LT.mstu(47)) goto 160
159  IF(nrem.GT.0) goto 210
160  njet=npre
161 
162 C...Reassign all particles to nearest jet. Sum up new jet momenta.
163  290 tsav=0.
164  psjt=0.
165  300 IF(mstu(46).LE.1) THEN
166  DO 310 i=n+1,n+njet
167  DO 310 j=1,4
168  310 v(i,j)=0.
169  DO 340 i=n+np+1,n+2*np
170  r2min=pss**2
171  DO 320 ijet=n+1,n+njet
172  IF(p(ijet,5).LT.rinit) goto 320
173  r2=r2t(i,ijet)
174  IF(r2.GE.r2min) goto 320
175  imin=ijet
176  r2min=r2
177  320 CONTINUE
178  k(i,4)=imin-n
179  DO 330 j=1,4
180  330 v(imin,j)=v(imin,j)+p(i,j)
181  340 CONTINUE
182  psjt=0.
183  DO 360 i=n+1,n+njet
184  DO 350 j=1,4
185  350 p(i,j)=v(i,j)
186  p(i,5)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
187  360 psjt=psjt+p(i,5)
188  ENDIF
189 
190 C...Find two closest jets.
191  r2min=2.*r2acc
192  DO 370 itry1=n+1,n+njet-1
193  DO 370 itry2=itry1+1,n+njet
194  IF(mstu(46).LE.2) r2=r2t(itry1,itry2)
195  IF(mstu(46).GE.3) r2=r2m(itry1,itry2)
196  IF(r2.GE.r2min) goto 370
197  imin1=itry1
198  imin2=itry2
199  r2min=r2
200  370 CONTINUE
201 
202 C...If allowed, join two closest jets and start over.
203  IF(njet.GT.mstu(47).AND.r2min.LT.r2acc) THEN
204  irec=min(imin1,imin2)
205  idel=max(imin1,imin2)
206  DO 380 j=1,4
207  380 p(irec,j)=p(imin1,j)+p(imin2,j)
208  p(irec,5)=sqrt(p(irec,1)**2+p(irec,2)**2+p(irec,3)**2)
209  DO 390 i=idel+1,n+njet
210  DO 390 j=1,5
211  390 p(i-1,j)=p(i,j)
212  IF(mstu(46).GE.2) THEN
213  DO 400 i=n+np+1,n+2*np
214  iori=n+k(i,4)
215  IF(iori.EQ.idel) k(i,4)=irec-n
216  400 IF(iori.GT.idel) k(i,4)=k(i,4)-1
217  ENDIF
218  njet=njet-1
219  goto 290
220 
221 C...Divide up broad jet if empty cluster in list of final ones.
222  ELSEIF(njet.EQ.mstu(47).AND.mstu(46).LE.1.AND.nloop.LE.2) THEN
223  DO 410 i=n+1,n+njet
224  410 k(i,5)=0
225  DO 420 i=n+np+1,n+2*np
226  420 k(n+k(i,4),5)=k(n+k(i,4),5)+1
227  iemp=0
228  DO 430 i=n+1,n+njet
229  430 IF(k(i,5).EQ.0) iemp=i
230  IF(iemp.NE.0) THEN
231  nloop=nloop+1
232  ispl=0
233  r2max=0.
234  DO 440 i=n+np+1,n+2*np
235  IF(k(n+k(i,4),5).LE.1.OR.p(i,5).LT.rinit) goto 440
236  ijet=n+k(i,4)
237  r2=r2t(i,ijet)
238  IF(r2.LE.r2max) goto 440
239  ispl=i
240  r2max=r2
241  440 CONTINUE
242  IF(ispl.NE.0) THEN
243  ijet=n+k(ispl,4)
244  DO 450 j=1,4
245  p(iemp,j)=p(ispl,j)
246  450 p(ijet,j)=p(ijet,j)-p(ispl,j)
247  p(iemp,5)=p(ispl,5)
248  p(ijet,5)=sqrt(p(ijet,1)**2+p(ijet,2)**2+p(ijet,3)**2)
249  IF(nloop.LE.2) goto 290
250  ENDIF
251  ENDIF
252  ENDIF
253 
254 C...If generalized thrust has not yet converged, continue iteration.
255  IF(mstu(46).LE.1.AND.nloop.LE.2.AND.psjt/pss.GT.tsav+paru(48))
256  &THEN
257  tsav=psjt/pss
258  goto 300
259  ENDIF
260 
261 C...Reorder jets according to energy.
262  DO 460 i=n+1,n+njet
263  DO 460 j=1,5
264  460 v(i,j)=p(i,j)
265  DO 490 inew=n+1,n+njet
266  pemax=0.
267  DO 470 itry=n+1,n+njet
268  IF(v(itry,4).LE.pemax) goto 470
269  imax=itry
270  pemax=v(itry,4)
271  470 CONTINUE
272  k(inew,1)=31
273  k(inew,2)=97
274  k(inew,3)=inew-n
275  k(inew,4)=0
276  DO 480 j=1,5
277  480 p(inew,j)=v(imax,j)
278  v(imax,4)=-1.
279  490 k(imax,5)=inew
280 
281 C...Clean up particle-jet assignments and jet information.
282  DO 500 i=n+np+1,n+2*np
283  iori=k(n+k(i,4),5)
284  k(i,4)=iori-n
285  IF(k(k(i,3),1).NE.3) k(k(i,3),4)=iori-n
286  k(iori,4)=k(iori,4)+1
287  500 CONTINUE
288  iemp=0
289  psjt=0.
290  DO 520 i=n+1,n+njet
291  k(i,5)=0
292  psjt=psjt+p(i,5)
293  p(i,5)=sqrt(max(p(i,4)**2-p(i,5)**2,0.))
294  DO 510 j=1,5
295  510 v(i,j)=0.
296  520 IF(k(i,4).EQ.0) iemp=i
297 
298 C...Select storing option. Output variables. Check for failure.
299  mstu(61)=n+1
300  mstu(62)=np
301  mstu(63)=npre
302  paru(61)=ps(5)
303  paru(62)=psjt/pss
304  paru(63)=sqrt(r2min)
305  IF(njet.LE.1) paru(63)=0.
306  IF(iemp.NE.0) THEN
307  CALL luerrm(8,'(LUCLUS:) failed to reconstruct as requested')
308  njet=-1
309  ENDIF
310  IF(mstu(43).LE.1) mstu(3)=njet
311  IF(mstu(43).GE.2) n=n+njet
312  nsav=njet
313 
314  RETURN
315  END