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