]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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(4000,5),P(4000,5),V(4000,5) | |
9 | COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) | |
10 | COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) | |
11 | SAVE /LUJETS/,/LUDAT1/,/LUDAT2/ | |
12 | DIMENSION PS(5) | |
13 | SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM | |
14 | ||
15 | C...Functions: distance measure in pT or (pseudo)mass. | |
16 | R2T(I1,I2)=(P(I1,5)*P(I2,5)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)- | |
17 | &P(I1,3)*P(I2,3))*2.*P(I1,5)*P(I2,5)/(0.0001+P(I1,5)+P(I2,5))**2 | |
18 | R2M(I1,I2)=2.*P(I1,4)*P(I2,4)*(1.-(P(I1,1)*P(I2,1)+P(I1,2)* | |
19 | &P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5))) | |
20 | ||
21 | C...If first time, reset. If reentering, skip preliminaries. | |
22 | IF(MSTU(48).LE.0) THEN | |
23 | NP=0 | |
24 | DO 100 J=1,5 | |
25 | PS(J)=0. | |
26 | 100 CONTINUE | |
27 | PSS=0. | |
28 | ELSE | |
29 | NJET=NSAV | |
30 | IF(MSTU(43).GE.2) N=N-NJET | |
31 | DO 110 I=N+1,N+NJET | |
32 | P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) | |
33 | 110 CONTINUE | |
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 300 | |
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 | P(N+NP,J)=P(I,J) | |
61 | 120 CONTINUE | |
62 | IF(MSTU(42).EQ.0) P(N+NP,5)=0. | |
63 | IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1) | |
64 | P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2) | |
65 | P(N+NP,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) | |
66 | DO 130 J=1,4 | |
67 | PS(J)=PS(J)+P(N+NP,J) | |
68 | 130 CONTINUE | |
69 | PSS=PSS+P(N+NP,5) | |
70 | 140 CONTINUE | |
71 | DO 160 I=N+1,N+NP | |
72 | K(I+NP,3)=K(I,3) | |
73 | DO 150 J=1,5 | |
74 | P(I+NP,J)=P(I,J) | |
75 | 150 CONTINUE | |
76 | 160 CONTINUE | |
77 | PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2)) | |
78 | ||
79 | C...Very low multiplicities not considered. | |
80 | IF(NP.LT.MSTU(47)) THEN | |
81 | CALL LUERRM(8,'(LUCLUS:) too few particles for analysis') | |
82 | NJET=-1 | |
83 | RETURN | |
84 | ENDIF | |
85 | ||
86 | C...Find precluster configuration. If too few jets, make harder cuts. | |
87 | NLOOP=0 | |
88 | IF(MSTU(46).LE.3) R2ACC=PARU(44)**2 | |
89 | IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2 | |
90 | RINIT=1.25*PARU(43) | |
91 | IF(NP.LE.MSTU(47)+2) RINIT=0. | |
92 | 170 RINIT=0.8*RINIT | |
93 | NPRE=0 | |
94 | NREM=NP | |
95 | DO 180 I=N+NP+1,N+2*NP | |
96 | K(I,4)=0 | |
97 | 180 CONTINUE | |
98 | ||
99 | C...Sum up small momentum region. Jet if enough absolute momentum. | |
100 | IF(MSTU(46).LE.2) THEN | |
101 | DO 190 J=1,4 | |
102 | P(N+1,J)=0. | |
103 | 190 CONTINUE | |
104 | DO 210 I=N+NP+1,N+2*NP | |
105 | IF(P(I,5).GT.2.*RINIT) GOTO 210 | |
106 | NREM=NREM-1 | |
107 | K(I,4)=1 | |
108 | DO 200 J=1,4 | |
109 | P(N+1,J)=P(N+1,J)+P(I,J) | |
110 | 200 CONTINUE | |
111 | 210 CONTINUE | |
112 | P(N+1,5)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2) | |
113 | IF(P(N+1,5).GT.2.*RINIT) NPRE=1 | |
114 | IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 170 | |
115 | IF(NREM.EQ.0) GOTO 170 | |
116 | ENDIF | |
117 | ||
118 | C...Find fastest remaining particle. | |
119 | 220 NPRE=NPRE+1 | |
120 | PMAX=0. | |
121 | DO 230 I=N+NP+1,N+2*NP | |
122 | IF(K(I,4).NE.0.OR.P(I,5).LE.PMAX) GOTO 230 | |
123 | IMAX=I | |
124 | PMAX=P(I,5) | |
125 | 230 CONTINUE | |
126 | DO 240 J=1,5 | |
127 | P(N+NPRE,J)=P(IMAX,J) | |
128 | 240 CONTINUE | |
129 | NREM=NREM-1 | |
130 | K(IMAX,4)=NPRE | |
131 | ||
132 | C...Sum up precluster around it according to pT separation. | |
133 | IF(MSTU(46).LE.2) THEN | |
134 | DO 260 I=N+NP+1,N+2*NP | |
135 | IF(K(I,4).NE.0) GOTO 260 | |
136 | R2=R2T(I,IMAX) | |
137 | IF(R2.GT.RINIT**2) GOTO 260 | |
138 | NREM=NREM-1 | |
139 | K(I,4)=NPRE | |
140 | DO 250 J=1,4 | |
141 | P(N+NPRE,J)=P(N+NPRE,J)+P(I,J) | |
142 | 250 CONTINUE | |
143 | 260 CONTINUE | |
144 | P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2) | |
145 | ||
146 | C...Sum up precluster around it according to mass separation. | |
147 | ELSE | |
148 | 270 IMIN=0 | |
149 | R2MIN=RINIT**2 | |
150 | DO 280 I=N+NP+1,N+2*NP | |
151 | IF(K(I,4).NE.0) GOTO 280 | |
152 | R2=R2M(I,N+NPRE) | |
153 | IF(R2.GE.R2MIN) GOTO 280 | |
154 | IMIN=I | |
155 | R2MIN=R2 | |
156 | 280 CONTINUE | |
157 | IF(IMIN.NE.0) THEN | |
158 | DO 290 J=1,4 | |
159 | P(N+NPRE,J)=P(N+NPRE,J)+P(IMIN,J) | |
160 | 290 CONTINUE | |
161 | P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2) | |
162 | NREM=NREM-1 | |
163 | K(IMIN,4)=NPRE | |
164 | GOTO 270 | |
165 | ENDIF | |
166 | ENDIF | |
167 | ||
168 | C...Check if more preclusters to be found. Start over if too few. | |
169 | IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 170 | |
170 | IF(NREM.GT.0) GOTO 220 | |
171 | NJET=NPRE | |
172 | ||
173 | C...Reassign all particles to nearest jet. Sum up new jet momenta. | |
174 | 300 TSAV=0. | |
175 | PSJT=0. | |
176 | 310 IF(MSTU(46).LE.1) THEN | |
177 | DO 330 I=N+1,N+NJET | |
178 | DO 320 J=1,4 | |
179 | V(I,J)=0. | |
180 | 320 CONTINUE | |
181 | 330 CONTINUE | |
182 | DO 360 I=N+NP+1,N+2*NP | |
183 | R2MIN=PSS**2 | |
184 | DO 340 IJET=N+1,N+NJET | |
185 | IF(P(IJET,5).LT.RINIT) GOTO 340 | |
186 | R2=R2T(I,IJET) | |
187 | IF(R2.GE.R2MIN) GOTO 340 | |
188 | IMIN=IJET | |
189 | R2MIN=R2 | |
190 | 340 CONTINUE | |
191 | K(I,4)=IMIN-N | |
192 | DO 350 J=1,4 | |
193 | V(IMIN,J)=V(IMIN,J)+P(I,J) | |
194 | 350 CONTINUE | |
195 | 360 CONTINUE | |
196 | PSJT=0. | |
197 | DO 380 I=N+1,N+NJET | |
198 | DO 370 J=1,4 | |
199 | P(I,J)=V(I,J) | |
200 | 370 CONTINUE | |
201 | P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) | |
202 | PSJT=PSJT+P(I,5) | |
203 | 380 CONTINUE | |
204 | ENDIF | |
205 | ||
206 | C...Find two closest jets. | |
207 | R2MIN=2.*MAX(R2ACC,PS(5)**2) | |
208 | DO 400 ITRY1=N+1,N+NJET-1 | |
209 | DO 390 ITRY2=ITRY1+1,N+NJET | |
210 | IF(MSTU(46).LE.2) R2=R2T(ITRY1,ITRY2) | |
211 | IF(MSTU(46).GE.3) R2=R2M(ITRY1,ITRY2) | |
212 | IF(R2.GE.R2MIN) GOTO 390 | |
213 | IMIN1=ITRY1 | |
214 | IMIN2=ITRY2 | |
215 | R2MIN=R2 | |
216 | 390 CONTINUE | |
217 | 400 CONTINUE | |
218 | ||
219 | C...If allowed, join two closest jets and start over. | |
220 | IF(NJET.GT.MSTU(47).AND.R2MIN.LT.R2ACC) THEN | |
221 | IREC=MIN(IMIN1,IMIN2) | |
222 | IDEL=MAX(IMIN1,IMIN2) | |
223 | DO 410 J=1,4 | |
224 | P(IREC,J)=P(IMIN1,J)+P(IMIN2,J) | |
225 | 410 CONTINUE | |
226 | P(IREC,5)=SQRT(P(IREC,1)**2+P(IREC,2)**2+P(IREC,3)**2) | |
227 | DO 430 I=IDEL+1,N+NJET | |
228 | DO 420 J=1,5 | |
229 | P(I-1,J)=P(I,J) | |
230 | 420 CONTINUE | |
231 | 430 CONTINUE | |
232 | IF(MSTU(46).GE.2) THEN | |
233 | DO 440 I=N+NP+1,N+2*NP | |
234 | IORI=N+K(I,4) | |
235 | IF(IORI.EQ.IDEL) K(I,4)=IREC-N | |
236 | IF(IORI.GT.IDEL) K(I,4)=K(I,4)-1 | |
237 | 440 CONTINUE | |
238 | ENDIF | |
239 | NJET=NJET-1 | |
240 | GOTO 300 | |
241 | ||
242 | C...Divide up broad jet if empty cluster in list of final ones. | |
243 | ELSEIF(NJET.EQ.MSTU(47).AND.MSTU(46).LE.1.AND.NLOOP.LE.2) THEN | |
244 | DO 450 I=N+1,N+NJET | |
245 | K(I,5)=0 | |
246 | 450 CONTINUE | |
247 | DO 460 I=N+NP+1,N+2*NP | |
248 | K(N+K(I,4),5)=K(N+K(I,4),5)+1 | |
249 | 460 CONTINUE | |
250 | IEMP=0 | |
251 | DO 470 I=N+1,N+NJET | |
252 | IF(K(I,5).EQ.0) IEMP=I | |
253 | 470 CONTINUE | |
254 | IF(IEMP.NE.0) THEN | |
255 | NLOOP=NLOOP+1 | |
256 | ISPL=0 | |
257 | R2MAX=0. | |
258 | DO 480 I=N+NP+1,N+2*NP | |
259 | IF(K(N+K(I,4),5).LE.1.OR.P(I,5).LT.RINIT) GOTO 480 | |
260 | IJET=N+K(I,4) | |
261 | R2=R2T(I,IJET) | |
262 | IF(R2.LE.R2MAX) GOTO 480 | |
263 | ISPL=I | |
264 | R2MAX=R2 | |
265 | 480 CONTINUE | |
266 | IF(ISPL.NE.0) THEN | |
267 | IJET=N+K(ISPL,4) | |
268 | DO 490 J=1,4 | |
269 | P(IEMP,J)=P(ISPL,J) | |
270 | P(IJET,J)=P(IJET,J)-P(ISPL,J) | |
271 | 490 CONTINUE | |
272 | P(IEMP,5)=P(ISPL,5) | |
273 | P(IJET,5)=SQRT(P(IJET,1)**2+P(IJET,2)**2+P(IJET,3)**2) | |
274 | IF(NLOOP.LE.2) GOTO 300 | |
275 | ENDIF | |
276 | ENDIF | |
277 | ENDIF | |
278 | ||
279 | C...If generalized thrust has not yet converged, continue iteration. | |
280 | IF(MSTU(46).LE.1.AND.NLOOP.LE.2.AND.PSJT/PSS.GT.TSAV+PARU(48)) | |
281 | &THEN | |
282 | TSAV=PSJT/PSS | |
283 | GOTO 310 | |
284 | ENDIF | |
285 | ||
286 | C...Reorder jets according to energy. | |
287 | DO 510 I=N+1,N+NJET | |
288 | DO 500 J=1,5 | |
289 | V(I,J)=P(I,J) | |
290 | 500 CONTINUE | |
291 | 510 CONTINUE | |
292 | DO 540 INEW=N+1,N+NJET | |
293 | PEMAX=0. | |
294 | DO 520 ITRY=N+1,N+NJET | |
295 | IF(V(ITRY,4).LE.PEMAX) GOTO 520 | |
296 | IMAX=ITRY | |
297 | PEMAX=V(ITRY,4) | |
298 | 520 CONTINUE | |
299 | K(INEW,1)=31 | |
300 | K(INEW,2)=97 | |
301 | K(INEW,3)=INEW-N | |
302 | K(INEW,4)=0 | |
303 | DO 530 J=1,5 | |
304 | P(INEW,J)=V(IMAX,J) | |
305 | 530 CONTINUE | |
306 | V(IMAX,4)=-1. | |
307 | K(IMAX,5)=INEW | |
308 | 540 CONTINUE | |
309 | ||
310 | C...Clean up particle-jet assignments and jet information. | |
311 | DO 550 I=N+NP+1,N+2*NP | |
312 | IORI=K(N+K(I,4),5) | |
313 | K(I,4)=IORI-N | |
314 | IF(K(K(I,3),1).NE.3) K(K(I,3),4)=IORI-N | |
315 | K(IORI,4)=K(IORI,4)+1 | |
316 | 550 CONTINUE | |
317 | IEMP=0 | |
318 | PSJT=0. | |
319 | DO 570 I=N+1,N+NJET | |
320 | K(I,5)=0 | |
321 | PSJT=PSJT+P(I,5) | |
322 | P(I,5)=SQRT(MAX(P(I,4)**2-P(I,5)**2,0.)) | |
323 | DO 560 J=1,5 | |
324 | V(I,J)=0. | |
325 | 560 CONTINUE | |
326 | IF(K(I,4).EQ.0) IEMP=I | |
327 | 570 CONTINUE | |
328 | ||
329 | C...Select storing option. Output variables. Check for failure. | |
330 | MSTU(61)=N+1 | |
331 | MSTU(62)=NP | |
332 | MSTU(63)=NPRE | |
333 | PARU(61)=PS(5) | |
334 | PARU(62)=PSJT/PSS | |
335 | PARU(63)=SQRT(R2MIN) | |
336 | IF(NJET.LE.1) PARU(63)=0. | |
337 | IF(IEMP.NE.0) THEN | |
338 | CALL LUERRM(8,'(LUCLUS:) failed to reconstruct as requested') | |
339 | NJET=-1 | |
340 | ENDIF | |
341 | IF(MSTU(43).LE.1) MSTU(3)=NJET | |
342 | IF(MSTU(43).GE.2) N=N+NJET | |
343 | NSAV=NJET | |
344 | ||
345 | RETURN | |
346 | END |