]>
Commit | Line | Data |
---|---|---|
99e5fe42 | 1 | SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM, |
2 | + MXJET,NJET,PJET,IPASS,IJMUL,IERR) | |
3 | *.********************************************************* | |
4 | *. ------ | |
5 | *. PXCONE | |
6 | *. ------ | |
7 | *. | |
8 | *.********** Pre Release Version 26.2.93 | |
9 | *. | |
10 | *. Driver for the Cone Jet finding algorithm of L.A. del Pozo. | |
11 | *. Based on algorithm from D.E. Soper. | |
12 | *. Finds jets inside cone of half angle CONER with energy > EPSLON. | |
13 | *. Jets which receive more than a fraction OVLIM of their energy from | |
14 | *. overlaps with other jets are excluded. | |
15 | *. Output jets are ordered in energy. | |
16 | *. If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt) | |
17 | *. Usage : | |
18 | *. | |
19 | *. INTEGER ITKDM,MXTRK | |
20 | *. PARAMETER (ITKDM=4.or.more,MXTRK=1.or.more) | |
21 | *. INTEGER MXJET, MXTRAK, MXPROT | |
22 | *. PARAMETER (MXJET=10,MXTRAK=500,MXPROT=500) | |
23 | *. INTEGER IPASS (MXTRAK),IJMUL (MXJET) | |
24 | *. INTEGER NTRAK,NJET,IERR,MODE | |
25 | *. DOUBLE PRECISION PTRAK (ITKDM,MXTRK),PJET (5,MXJET) | |
26 | *. DOUBLE PRECISION CONER, EPSLON, OVLIM | |
27 | *. NTRAK = 1.to.MXTRAK | |
28 | *. CONER = ... | |
29 | *. EPSLON = ... | |
30 | *. OVLIM = ... | |
31 | *. CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET, | |
32 | *. + NJET,PJET,IPASS,IJMUL,IERR) | |
33 | *. | |
34 | *. INPUT : MODE 1=>e+e-, 2=>hadron-hadron | |
35 | *. INPUT : NTRAK Number of particles | |
36 | *. INPUT : ITKDM First dimension of PTRAK array | |
37 | *. INPUT : PTRAK Array of particle 4-momenta (Px,Py,Pz,E) | |
38 | *. INPUT : CONER Cone size (half angle) in radians | |
39 | *. INPUT : EPSLON Minimum Jet energy (GeV) | |
40 | *. INPUT : OVLIM Maximum fraction of overlap energy in a jet | |
41 | *. INPUT : MXJET Maximum possible number of jets | |
42 | *. OUTPUT : NJET Number of jets found | |
43 | *. OUTPUT : PJET 5-vectors of jets | |
44 | *. OUTPUT : IPASS(k) Particle k belongs to jet number IPASS(k) | |
45 | *. IPASS = -1 if not assosciated to a jet | |
46 | *. OUTPUT : IJMUL(i) Jet i contains IJMUL(i) particles | |
47 | *. OUTPUT : IERR = 0 if all is OK ; = -1 otherwise | |
48 | *. | |
49 | *. CALLS : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP | |
50 | *. CALLED : User | |
51 | *. | |
52 | *. AUTHOR : L.A. del Pozo | |
53 | *. CREATED : 26-Feb-93 | |
54 | *. LAST MOD : 2-Mar-93 | |
55 | *. | |
56 | *. Modification Log. | |
57 | *. 2-Jan-97: M Wobisch - fix bug concerning COS2R in eta phi mode | |
58 | *. 4-Apr-93: M H Seymour - Change 2d arrays to 1d in PXTRY & PXNEW | |
59 | *. 2-Apr-93: M H Seymour - Major changes to add boost-invariant mode | |
60 | *. 1-Apr-93: M H Seymour - Increase all array sizes | |
61 | *. 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION | |
62 | *. 30-Mar-93: M H Seymour - Change OVLIM into an input parameter | |
63 | *. 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP | |
64 | *. 1-Mar-93: L A del Pozo - Remove Cern library routine calls | |
65 | *. 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon | |
66 | *. | |
67 | *.********************************************************* | |
68 | C+SEQ,DECLARE. | |
69 | *** External Arrays | |
70 | INTEGER ITKDM,MXJET,NTRAK,NJET,IERR,MODE | |
71 | INTEGER IPASS (NTRAK),IJMUL (MXJET) | |
72 | DOUBLE PRECISION PTRAK (ITKDM,NTRAK),PJET (5,MXJET), | |
73 | + CONER, EPSLON, OVLIM | |
74 | *** Internal Arrays | |
75 | INTEGER MXPROT, MXTRAK | |
76 | PARAMETER (MXPROT=5000, MXTRAK=5000) | |
77 | DOUBLE PRECISION PP(4,MXTRAK), PU(3,MXTRAK), PJ(4,MXPROT) | |
78 | LOGICAL JETLIS(MXPROT,MXTRAK) | |
79 | *** Used in the routine. | |
80 | DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ, | |
81 | + COSVAL,PXMDPI | |
82 | cMWobisch | |
83 | DOUBLE PRECISION RSEP | |
84 | cMWobisch | |
85 | LOGICAL UNSTBL | |
86 | INTEGER I,J,N,MU,N1,N2, ITERR | |
87 | INTEGER NCALL, NPRINT | |
88 | DOUBLE PRECISION ROLD, EPSOLD, OVOLD | |
89 | SAVE NCALL,NPRINT,ROLD, EPSOLD, OVOLD | |
90 | DATA NCALL,NPRINT /0,0/ | |
91 | DATA ROLD,EPSOLD,OVOLD/0.,0.,0./ | |
92 | ||
93 | cMWobisch | |
94 | c*************************************** | |
95 | RSEP = 2D0 | |
96 | c*************************************** | |
97 | cMWobisch | |
98 | IERR=0 | |
99 | * | |
100 | *** INITIALIZE | |
101 | IF(NCALL.LE.0) THEN | |
102 | ROLD = 0. | |
103 | EPSOLD = 0. | |
104 | OVOLD = 0. | |
105 | ENDIF | |
106 | NCALL = NCALL + 1 | |
107 | * | |
108 | *** Print welcome and Jetfinder parameters | |
109 | IF((CONER.NE.ROLD .OR. EPSLON.NE.EPSOLD .OR. OVLIM.NE.OVOLD) | |
110 | + .AND. NPRINT.LE.10) THEN | |
111 | WRITE (6,*) | |
112 | WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********' | |
113 | WRITE (6,*) ' Written by Luis Del Pozo of OPAL' | |
114 | WRITE (6,*) ' Modified for eta-phi by Mike Seymour' | |
115 | WRITE(6,1000)' Cone Size R = ',CONER,' Radians' | |
116 | WRITE(6,1001)' Min Jet energy Epsilon = ',EPSLON,' GeV' | |
117 | WRITE(6,1002)' Overlap fraction parameter = ',OVLIM | |
118 | WRITE (6,*) ' ***********************************************' | |
119 | cMWobisch | |
120 | IF (RSEP .lt. 1.999) THEN | |
121 | WRITE(6,*) ' ' | |
122 | WRITE (6,*) ' ******************************************' | |
123 | WRITE (6,*) ' ******************************************' | |
124 | WRITE(6,*) ' M Wobisch: private change !!!!!!!!!!!! ' | |
125 | WRITE(6,*) ' Rsep is set to ',RSEP | |
126 | WRITE(6,*) ' this is ONLY meaningful in a NLO calculation' | |
127 | WRITE(6,*) ' ------------------------ ' | |
128 | WRITE(6,*) ' please check what you''re doing!!' | |
129 | WRITE(6,*) ' or ask: Markus.Wobisch@desy.de --' | |
130 | WRITE (6,*) ' ******************************************' | |
131 | WRITE (6,*) ' ******************************************' | |
132 | WRITE (6,*) ' ******************************************' | |
133 | WRITE(6,*) ' ' | |
134 | WRITE(6,*) ' ' | |
135 | ENDIF | |
136 | cMWobisch | |
137 | ||
138 | WRITE (6,*) | |
139 | 1000 FORMAT(A18,F5.2,A10) | |
140 | 1001 FORMAT(A29,F5.2,A5) | |
141 | 1002 FORMAT(A33,F5.2) | |
142 | NPRINT = NPRINT + 1 | |
143 | ROLD=CONER | |
144 | EPSOLD=EPSLON | |
145 | OVOLD=OVLIM | |
146 | ENDIF | |
147 | * | |
148 | *** Copy calling array PTRAK to internal array PP(4,NTRAK) | |
149 | * | |
150 | IF (NTRAK .GT. MXTRAK) THEN | |
151 | WRITE (6,*) ' PXCONE: Ntrak too large: ',NTRAK | |
152 | IERR=-1 | |
153 | RETURN | |
154 | ENDIF | |
155 | IF (MODE.NE.2) THEN | |
156 | DO 100 I=1, NTRAK | |
157 | DO 101 J=1,4 | |
158 | PP(J,I)=PTRAK(J,I) | |
159 | 101 CONTINUE | |
160 | 100 CONTINUE | |
161 | ELSE | |
162 | *** Converting to eta,phi,pt if necessary | |
163 | DO 104 I=1,NTRAK | |
164 | PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2 | |
165 | PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2 | |
166 | IF (PTSQ.LE.4.25E-18*PPSQ) THEN | |
167 | PP(1,I)=20 | |
168 | ELSE | |
169 | PP(1,I)=0.5*LOG(PPSQ/PTSQ) | |
170 | ENDIF | |
171 | PP(1,I)=SIGN(PP(1,I),PTRAK(3,I)) | |
172 | IF (PTSQ.EQ.0) THEN | |
173 | PP(2,I)=0 | |
174 | ELSE | |
175 | PP(2,I)=ATAN2(PTRAK(2,I),PTRAK(1,I)) | |
176 | ENDIF | |
177 | PP(3,I)=0 | |
178 | PP(4,I)=SQRT(PTSQ) | |
179 | PU(1,I)=PP(1,I) | |
180 | PU(2,I)=PP(2,I) | |
181 | PU(3,I)=PP(3,I) | |
182 | 104 CONTINUE | |
183 | ENDIF | |
184 | * | |
185 | *** Zero output variables | |
186 | * | |
187 | NJET=0 | |
188 | DO 102 I = 1, NTRAK | |
189 | DO 103 J = 1, MXPROT | |
190 | JETLIS(J,I) = .FALSE. | |
191 | 103 CONTINUE | |
192 | 102 CONTINUE | |
193 | CALL PXZERV(4*MXPROT,PJ) | |
194 | CALL PXZERI(MXJET,IJMUL) | |
195 | * | |
196 | IF (MODE.NE.2) THEN | |
197 | COSR = COS(CONER) | |
198 | COS2R = COS(CONER) | |
199 | ELSE | |
200 | *** Purely for convenience, work in terms of 1-R**2 | |
201 | COSR = 1-CONER**2 | |
202 | cMW -- select Rsep: 1-(Rsep*CONER)**2 | |
203 | COS2R = 1-(RSEP*CONER)**2 | |
204 | cORIGINAL COS2R = 1-(2*CONER)**2 | |
205 | ENDIF | |
206 | UNSTBL = .FALSE. | |
207 | IF (MODE.NE.2) THEN | |
208 | CALL PXUVEC(NTRAK,PP,PU,IERR) | |
209 | IF (IERR .NE. 0) RETURN | |
210 | ENDIF | |
211 | *** Look for jets using particle diretions as seed axes | |
212 | * | |
213 | DO 110 N = 1,NTRAK | |
214 | DO 120 MU = 1,3 | |
215 | VSEED(MU) = PU(MU,N) | |
216 | 120 CONTINUE | |
217 | CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED, | |
218 | & NJET,JETLIS,PJ,UNSTBL,IERR) | |
219 | IF (IERR .NE. 0) RETURN | |
220 | 110 CONTINUE | |
221 | ||
222 | cMW - for Rsep=1 goto 145 | |
223 | c GOTO 145 | |
224 | ||
225 | *** Now look between all pairs of jets as seed axes. | |
226 | DO 140 N1 = 1,NJET-1 | |
227 | VEC1(1)=PJ(1,N1) | |
228 | VEC1(2)=PJ(2,N1) | |
229 | VEC1(3)=PJ(3,N1) | |
230 | IF (MODE.NE.2) CALL PXNORV(3,VEC1,VEC1,ITERR) | |
231 | DO 150 N2 = N1+1,NJET | |
232 | VEC2(1)=PJ(1,N2) | |
233 | VEC2(2)=PJ(2,N2) | |
234 | VEC2(3)=PJ(3,N2) | |
235 | IF (MODE.NE.2) CALL PXNORV(3,VEC2,VEC2,ITERR) | |
236 | CALL PXADDV(3,VEC1,VEC2,VSEED,ITERR) | |
237 | IF (MODE.NE.2) THEN | |
238 | CALL PXNORV(3,VSEED,VSEED,ITERR) | |
239 | ELSE | |
240 | VSEED(1)=VSEED(1)/2 | |
241 | VSEED(2)=VSEED(2)/2 | |
242 | ENDIF | |
243 | C---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART | |
244 | IF (MODE.NE.2) THEN | |
245 | COSVAL=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3) | |
246 | ELSE | |
247 | IF (ABS(VEC1(1)).GE.20.OR.ABS(VEC2(1)).GE.20) THEN | |
248 | COSVAL=-1000 | |
249 | ELSE | |
250 | COSVAL=1- | |
251 | + ((VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2) | |
252 | ENDIF | |
253 | ENDIF | |
254 | ||
255 | IF (COSVAL.LE.COSR.AND.COSVAL.GE.COS2R) | |
256 | + CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, | |
257 | + JETLIS,PJ,UNSTBL,IERR) | |
258 | c CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, | |
259 | c + JETLIS,PJ,UNSTBL,IERR) | |
260 | IF (IERR .NE. 0) RETURN | |
261 | 150 CONTINUE | |
262 | 140 CONTINUE | |
263 | IF (UNSTBL) THEN | |
264 | IERR=-1 | |
265 | WRITE (6,*) ' PXCONE: Too many iterations to find a proto-jet' | |
266 | RETURN | |
267 | ENDIF | |
268 | ||
269 | 145 CONTINUE | |
270 | *** Now put the jet list into order by jet energy, eliminating jets | |
271 | *** with energy less than EPSLON. | |
272 | CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ) | |
273 | * | |
274 | *** Take care of jet overlaps | |
275 | CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM) | |
276 | * | |
277 | *** Order jets again as some have been eliminated, or lost energy. | |
278 | CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ) | |
279 | * | |
280 | *** All done!, Copy output into output arrays | |
281 | IF (NJET .GT. MXJET) THEN | |
282 | WRITE (6,*) ' PXCONE: Found more than MXJET jets' | |
283 | IERR=-1 | |
284 | GOTO 99 | |
285 | ENDIF | |
286 | IF (MODE.NE.2) THEN | |
287 | DO 300 I=1, NJET | |
288 | DO 310 J=1,4 | |
289 | PJET(J,I)=PJ(J,I) | |
290 | 310 CONTINUE | |
291 | 300 CONTINUE | |
292 | ELSE | |
293 | DO 315 I=1, NJET | |
294 | PJET(1,I)=PJ(4,I)*COS(PJ(2,I)) | |
295 | PJET(2,I)=PJ(4,I)*SIN(PJ(2,I)) | |
296 | PJET(3,I)=PJ(4,I)*SINH(PJ(1,I)) | |
297 | PJET(4,I)=PJ(4,I)*COSH(PJ(1,I)) | |
298 | 315 CONTINUE | |
299 | ENDIF | |
300 | DO 320 I=1, NTRAK | |
301 | IPASS(I)=-1 | |
302 | DO 330 J=1, NJET | |
303 | IF (JETLIS(J,I)) THEN | |
304 | IJMUL(J)=IJMUL(J)+1 | |
305 | IPASS(I)=J | |
306 | ENDIF | |
307 | 330 CONTINUE | |
308 | 320 CONTINUE | |
309 | 99 RETURN | |
310 | END | |
311 | *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper | |
312 | *-- Author : | |
313 | C----------------------------------------------------------------------- | |
314 | SUBROUTINE PXNORV(N,A,B,ITERR) | |
315 | INTEGER I,N,ITERR | |
316 | DOUBLE PRECISION A(N),B(N),C | |
317 | C=0 | |
318 | DO 10 I=1,N | |
319 | C=C+A(I)**2 | |
320 | 10 CONTINUE | |
321 | IF (C.LE.0) RETURN | |
322 | C=1/SQRT(C) | |
323 | DO 20 I=1,N | |
324 | B(I)=A(I)*C | |
325 | 20 CONTINUE | |
326 | END | |
327 | *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper | |
328 | *CMZ : 1.06/00 15/03/94 12.17.46 by P. Schleper | |
329 | *-- Author : | |
330 | * | |
331 | C+DECK,PXOLAP. | |
332 | SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM) | |
333 | * | |
334 | *** Looks for particles assigned to more than 1 jet, and reassigns them | |
335 | *** If more than a fraction OVLIM of a jet's energy is contained in | |
336 | *** higher energy jets, that jet is neglected. | |
337 | *** Particles assigned to the jet closest in angle (a la CDF, Snowmass). | |
338 | C+SEQ,DECLARE. | |
339 | INTEGER MXTRAK, MXPROT | |
340 | PARAMETER (MXTRAK=5000,MXPROT=5000) | |
341 | INTEGER NJET, NTRAK, MODE | |
342 | LOGICAL JETLIS(MXPROT,MXTRAK) | |
343 | DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI | |
344 | INTEGER I,J,N,MU | |
345 | LOGICAL OVELAP | |
346 | DOUBLE PRECISION EOVER | |
347 | DOUBLE PRECISION OVLIM | |
348 | INTEGER ITERR, IJMIN, IJET(MXPROT), NJ | |
349 | DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN | |
350 | DATA IJMIN/0/ | |
351 | * | |
352 | IF (NJET.LE.1) RETURN | |
353 | *** Look for jets with large overlaps with higher energy jets. | |
354 | DO 100 I = 2,NJET | |
355 | *** Find overlap energy between jets I and all higher energy jets. | |
356 | EOVER = 0.0 | |
357 | DO 110 N = 1,NTRAK | |
358 | OVELAP = .FALSE. | |
359 | DO 120 J= 1,I-1 | |
360 | IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN | |
361 | OVELAP = .TRUE. | |
362 | ENDIF | |
363 | 120 CONTINUE | |
364 | IF (OVELAP) THEN | |
365 | EOVER = EOVER + PP(4,N) | |
366 | ENDIF | |
367 | 110 CONTINUE | |
368 | *** Is the fraction of energy shared larger than OVLIM? | |
369 | IF (EOVER.GT.OVLIM*PJ(4,I)) THEN | |
370 | *** De-assign all particles from Jet I | |
371 | DO 130 N = 1,NTRAK | |
372 | JETLIS(I,N) = .FALSE. | |
373 | 130 CONTINUE | |
374 | ENDIF | |
375 | 100 CONTINUE | |
376 | *** Now there are no big overlaps, assign every particle in | |
377 | *** more than 1 jet to the closet jet. | |
378 | *** Any particles now in more than 1 jet are assigned to the CLOSET | |
379 | *** jet (in angle). | |
380 | DO 140 I=1,NTRAK | |
381 | NJ=0 | |
382 | DO 150 J=1, NJET | |
383 | IF(JETLIS(J,I)) THEN | |
384 | NJ=NJ+1 | |
385 | IJET(NJ)=J | |
386 | ENDIF | |
387 | 150 CONTINUE | |
388 | IF (NJ .GT. 1) THEN | |
389 | *** Particle in > 1 jet - calc angles... | |
390 | VEC1(1)=PP(1,I) | |
391 | VEC1(2)=PP(2,I) | |
392 | VEC1(3)=PP(3,I) | |
393 | THMIN=0. | |
394 | DO 160 J=1,NJ | |
395 | VEC2(1)=PJ(1,IJET(J)) | |
396 | VEC2(2)=PJ(2,IJET(J)) | |
397 | VEC2(3)=PJ(3,IJET(J)) | |
398 | IF (MODE.NE.2) THEN | |
399 | CALL PXANG3(VEC1,VEC2,COST,THET,ITERR) | |
400 | ELSE | |
401 | THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2 | |
402 | ENDIF | |
403 | IF (J .EQ. 1) THEN | |
404 | THMIN=THET | |
405 | IJMIN=IJET(J) | |
406 | ELSEIF (THET .LT. THMIN) THEN | |
407 | THMIN=THET | |
408 | IJMIN=IJET(J) | |
409 | ENDIF | |
410 | 160 CONTINUE | |
411 | *** Assign track to IJMIN | |
412 | DO 170 J=1,NJET | |
413 | JETLIS(J,I) = .FALSE. | |
414 | 170 CONTINUE | |
415 | JETLIS(IJMIN,I)=.TRUE. | |
416 | ENDIF | |
417 | 140 CONTINUE | |
418 | *** Recompute PJ | |
419 | DO 200 I = 1,NJET | |
420 | DO 210 MU = 1,4 | |
421 | PJ(MU,I) = 0.0 | |
422 | 210 CONTINUE | |
423 | DO 220 N = 1,NTRAK | |
424 | IF( JETLIS(I,N) ) THEN | |
425 | IF (MODE.NE.2) THEN | |
426 | DO 230 MU = 1,4 | |
427 | PJ(MU,I) = PJ(MU,I) + PP(MU,N) | |
428 | 230 CONTINUE | |
429 | ELSE | |
430 | PJ(1,I)=PJ(1,I) | |
431 | + + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I)) | |
432 | PJ(2,I)=PJ(2,I) | |
433 | + + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I)) | |
434 | PJ(4,I)=PJ(4,I)+PP(4,N) | |
435 | ENDIF | |
436 | ENDIF | |
437 | 220 CONTINUE | |
438 | 200 CONTINUE | |
439 | RETURN | |
440 | END | |
441 | *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper | |
442 | *CMZ : 1.06/00 14/03/94 15.37.45 by P. Schleper | |
443 | *-- Author : | |
444 | * | |
445 | C+DECK,PXORD. | |
446 | SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ) | |
447 | * | |
448 | *** Routine to put jets into order and eliminate tose less than EPSLON | |
449 | C+SEQ,DECLARE. | |
450 | INTEGER MXTRAK,MXPROT | |
451 | PARAMETER (MXTRAK=5000,MXPROT=5000) | |
452 | INTEGER I, J, INDEX(MXPROT) | |
453 | DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT) | |
454 | INTEGER NJET,NTRAK | |
455 | LOGICAL JETLIS(MXPROT,MXTRAK) | |
456 | LOGICAL LOGTMP(MXPROT,MXTRAK) | |
457 | DOUBLE PRECISION EPSLON,PJ(4,MXPROT) | |
458 | *** Puts jets in order of energy: 1 = highest energy etc. | |
459 | *** Then Eliminate jets with energy below EPSLON | |
460 | * | |
461 | *** Copy input arrays. | |
462 | DO 100 I=1,NJET | |
463 | DO 110 J=1,4 | |
464 | PTEMP(J,I)=PJ(J,I) | |
465 | 110 CONTINUE | |
466 | DO 120 J=1,NTRAK | |
467 | LOGTMP(I,J)=JETLIS(I,J) | |
468 | 120 CONTINUE | |
469 | 100 CONTINUE | |
470 | DO 150 I=1,NJET | |
471 | ELIST(I)=PJ(4,I) | |
472 | 150 CONTINUE | |
473 | *** Sort the energies... | |
474 | CALL PXSORV(NJET,ELIST,INDEX,'I') | |
475 | *** Fill PJ and JETLIS according to sort ( sort is in ascending order!!) | |
476 | DO 200 I=1, NJET | |
477 | DO 210 J=1,4 | |
478 | PJ(J,I)=PTEMP(J,INDEX(NJET+1-I)) | |
479 | 210 CONTINUE | |
480 | DO 220 J=1,NTRAK | |
481 | JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J) | |
482 | 220 CONTINUE | |
483 | 200 CONTINUE | |
484 | ** Jets are now in order | |
485 | *** Now eliminate jets with less than Epsilon energy | |
486 | DO 300, I=1, NJET | |
487 | IF (PJ(4,I) .LT. EPSLON) THEN | |
488 | NJET=NJET-1 | |
489 | PJ(4,I)=0. | |
490 | ENDIF | |
491 | 300 CONTINUE | |
492 | RETURN | |
493 | END | |
494 | ||
495 | ******************************************************************** | |
496 | *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper | |
497 | *CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper | |
498 | *-- Author : | |
499 | C+DECK,PXSEAR. | |
500 | SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, | |
501 | + JETLIS,PJ,UNSTBL,IERR) | |
502 | * | |
503 | C+SEQ,DECLARE. | |
504 | INTEGER MXTRAK, MXPROT | |
505 | PARAMETER (MXTRAK=5000,MXPROT=5000) | |
506 | INTEGER NTRAK, IERR, MODE | |
507 | DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3) | |
508 | LOGICAL UNSTBL | |
509 | LOGICAL JETLIS(MXPROT,MXTRAK) | |
510 | INTEGER NJET | |
511 | DOUBLE PRECISION PJ(4,MXPROT) | |
512 | *** Using VSEED as a trial axis , look for a stable jet. | |
513 | *** Check stable jets against those already found and add to PJ. | |
514 | *** Will try up to MXITER iterations to get a stable set of particles | |
515 | *** in the cone. | |
516 | INTEGER MU,N,ITER | |
517 | LOGICAL PXSAME,PXNEW,OK | |
518 | LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK) | |
519 | DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4) | |
520 | INTEGER MXITER | |
521 | PARAMETER(MXITER = 30) | |
522 | * | |
523 | DO 100 MU=1,3 | |
524 | OAXIS(MU) = VSEED(MU) | |
525 | 100 CONTINUE | |
526 | DO 110 N = 1,NTRAK | |
527 | OLDLIS(N) = .FALSE. | |
528 | 110 CONTINUE | |
529 | DO 120 ITER = 1,MXITER | |
530 | CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK) | |
531 | *** Return immediately if there were no particles in the cone. | |
532 | IF (.NOT.OK) THEN | |
533 | RETURN | |
534 | ENDIF | |
535 | IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN | |
536 | *** We have a stable jet. | |
537 | IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN | |
538 | *** And the jet is a new one. So add it to our arrays. | |
539 | *** Check arrays are big anough... | |
540 | IF (NJET .EQ. MXPROT) THEN | |
541 | WRITE (6,*) ' PXCONE: Found more than MXPROT proto-jets' | |
542 | IERR = -1 | |
543 | RETURN | |
544 | ENDIF | |
545 | NJET = NJET + 1 | |
546 | DO 130 N = 1,NTRAK | |
547 | JETLIS(NJET,N) = NEWLIS(N) | |
548 | 130 CONTINUE | |
549 | DO 140 MU=1,4 | |
550 | PJ(MU,NJET)=PNEW(MU) | |
551 | 140 CONTINUE | |
552 | ENDIF | |
553 | RETURN | |
554 | ENDIF | |
555 | *** The jet was not stable, so we iterate again | |
556 | DO 150 N=1,NTRAK | |
557 | OLDLIS(N)=NEWLIS(N) | |
558 | 150 CONTINUE | |
559 | DO 160 MU=1,3 | |
560 | OAXIS(MU)=NAXIS(MU) | |
561 | 160 CONTINUE | |
562 | 120 CONTINUE | |
563 | UNSTBL = .TRUE. | |
564 | RETURN | |
565 | END | |
566 | *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper | |
567 | *-- Author : | |
568 | C----------------------------------------------------------------------- | |
569 | SUBROUTINE PXSORV(N,A,K,OPT) | |
570 | C Sort A(N) into ascending order | |
571 | C OPT = 'I' : return index array K only | |
572 | C OTHERWISE : return sorted A and index array K | |
573 | C----------------------------------------------------------------------- | |
574 | INTEGER NMAX | |
575 | PARAMETER (NMAX=5000) | |
576 | * | |
577 | * INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX) | |
578 | *LUND | |
579 | INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX) | |
580 | CHARACTER OPT | |
581 | * | |
582 | * DOUBLE PRECISION A(N),B(NMAX) | |
583 | DOUBLE PRECISION A(NMAX),B(NMAX) | |
584 | *LUND | |
585 | IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV' | |
586 | IL(1)=0 | |
587 | IR(1)=0 | |
588 | DO 10 I=2,N | |
589 | IL(I)=0 | |
590 | IR(I)=0 | |
591 | J=1 | |
592 | 2 IF(A(I).GT.A(J)) GO TO 5 | |
593 | 3 IF(IL(J).EQ.0) GO TO 4 | |
594 | J=IL(J) | |
595 | GO TO 2 | |
596 | 4 IR(I)=-J | |
597 | IL(J)=I | |
598 | GO TO 10 | |
599 | 5 IF(IR(J).LE.0) GO TO 6 | |
600 | J=IR(J) | |
601 | GO TO 2 | |
602 | 6 IR(I)=IR(J) | |
603 | IR(J)=I | |
604 | 10 CONTINUE | |
605 | I=1 | |
606 | J=1 | |
607 | GO TO 8 | |
608 | 20 J=IL(J) | |
609 | 8 IF(IL(J).GT.0) GO TO 20 | |
610 | 9 K(I)=J | |
611 | B(I)=A(J) | |
612 | I=I+1 | |
613 | IF(IR(J)) 12,30,13 | |
614 | 13 J=IR(J) | |
615 | GO TO 8 | |
616 | 12 J=-IR(J) | |
617 | GO TO 9 | |
618 | 30 IF(OPT.EQ.'I') RETURN | |
619 | DO 31 I=1,N | |
620 | 31 A(I)=B(I) | |
621 | 999 END | |
622 | ||
623 | ********************************************************************* | |
624 | *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper | |
625 | *CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper | |
626 | *-- Author : | |
627 | * | |
628 | C+DECK,PXTRY. | |
629 | SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS, | |
630 | + PNEW,NEWLIS,OK) | |
631 | * | |
632 | C+SEQ,DECLARE. | |
633 | INTEGER MXTRAK | |
634 | PARAMETER (MXTRAK=5000) | |
635 | INTEGER NTRAK,MODE | |
636 | *** Note that although PU and PP are assumed to be 2d arrays, they | |
637 | *** are used as 1d in this routine for efficiency | |
638 | DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI | |
639 | LOGICAL OK | |
640 | LOGICAL NEWLIS(MXTRAK) | |
641 | DOUBLE PRECISION NAXIS(3),PNEW(4) | |
642 | *** Finds all particles in cone of size COSR about OAXIS direction. | |
643 | *** Calculates 4-momentum sum of all particles in cone (PNEW) , and | |
644 | *** returns this as new jet axis NAXIS (Both unit Vectors) | |
645 | INTEGER N,MU,NPU,NPP | |
646 | DOUBLE PRECISION COSVAL,NORMSQ,NORM | |
647 | * | |
648 | OK = .FALSE. | |
649 | DO 100 MU=1,4 | |
650 | PNEW(MU)=0.0 | |
651 | 100 CONTINUE | |
652 | NPU=-3 | |
653 | NPP=-4 | |
654 | DO 110 N=1,NTRAK | |
655 | NPU=NPU+3 | |
656 | NPP=NPP+4 | |
657 | IF (MODE.NE.2) THEN | |
658 | COSVAL=0.0 | |
659 | DO 120 MU=1,3 | |
660 | COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU) | |
661 | 120 CONTINUE | |
662 | ELSE | |
663 | IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN | |
664 | COSVAL=-1000 | |
665 | ELSE | |
666 | COSVAL=1- | |
667 | + ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2) | |
668 | ENDIF | |
669 | ENDIF | |
670 | IF (COSVAL.GE.COSR)THEN | |
671 | NEWLIS(N) = .TRUE. | |
672 | OK = .TRUE. | |
673 | IF (MODE.NE.2) THEN | |
674 | DO 130 MU=1,4 | |
675 | PNEW(MU) = PNEW(MU) + PP(MU+NPP) | |
676 | 130 CONTINUE | |
677 | ELSE | |
678 | PNEW(1)=PNEW(1) | |
679 | + + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1)) | |
680 | PNEW(2)=PNEW(2) | |
681 | + + PP(4+NPP)/(PP(4+NPP)+PNEW(4)) | |
682 | + *PXMDPI(PP(2+NPP)-PNEW(2)) | |
683 | PNEW(4)=PNEW(4)+PP(4+NPP) | |
684 | ENDIF | |
685 | ELSE | |
686 | NEWLIS(N)=.FALSE. | |
687 | ENDIF | |
688 | 110 CONTINUE | |
689 | *** If there are particles in the cone, calc new jet axis | |
690 | IF (OK) THEN | |
691 | IF (MODE.NE.2) THEN | |
692 | NORMSQ = 0.0 | |
693 | DO 140 MU = 1,3 | |
694 | NORMSQ = NORMSQ + PNEW(MU)**2 | |
695 | 140 CONTINUE | |
696 | NORM = SQRT(NORMSQ) | |
697 | ELSE | |
698 | NORM = 1 | |
699 | ENDIF | |
700 | DO 150 MU=1,3 | |
701 | NAXIS(MU) = PNEW(MU)/NORM | |
702 | 150 CONTINUE | |
703 | ENDIF | |
704 | RETURN | |
705 | END | |
706 | ||
707 | ********************************************************************* | |
708 | *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper | |
709 | *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper | |
710 | *-- Author : | |
711 | C+DECK,PXUVEC. | |
712 | * | |
713 | SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR) | |
714 | * | |
715 | *** Routine to calculate unit vectors PU of all particles PP | |
716 | C+SEQ,DECLARE. | |
717 | INTEGER MXTRAK | |
718 | PARAMETER (MXTRAK=5000) | |
719 | INTEGER NTRAK, IERR | |
720 | DOUBLE PRECISION PP(4,MXTRAK) | |
721 | DOUBLE PRECISION PU(3,MXTRAK) | |
722 | INTEGER N,MU | |
723 | DOUBLE PRECISION MAG | |
724 | DO 100 N=1,NTRAK | |
725 | MAG=0.0 | |
726 | DO 110 MU=1,3 | |
727 | MAG=MAG+PP(MU,N)**2 | |
728 | 110 CONTINUE | |
729 | MAG=SQRT(MAG) | |
730 | IF (MAG.EQ.0.0) THEN | |
731 | WRITE(6,*)' PXCONE: An input particle has zero mod(p)' | |
732 | IERR=-1 | |
733 | RETURN | |
734 | ENDIF | |
735 | DO 120 MU=1,3 | |
736 | PU(MU,N)=PP(MU,N)/MAG | |
737 | 120 CONTINUE | |
738 | 100 CONTINUE | |
739 | RETURN | |
740 | END | |
741 | *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper | |
742 | *-- Author : | |
743 | C----------------------------------------------------------------------- | |
744 | SUBROUTINE PXZERI(N,A) | |
745 | INTEGER I,N,A(N) | |
746 | DO 10 I=1,N | |
747 | A(I)=0 | |
748 | 10 CONTINUE | |
749 | END | |
750 | *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper | |
751 | *-- Author : | |
752 | C----------------------------------------------------------------------- | |
753 | C This is a set of routines written by Mike Seymour to provide the | |
754 | C services presumably normally provided by standard OPAL routines | |
755 | C PXZERV zeroes a vector | |
756 | C PXZERI zeroes a vector of integers | |
757 | C PXNORV normalizes a vector | |
758 | C PXADDV adds two vectors | |
759 | C PXSORV sorts a vector (copied from HERWIG) | |
760 | C PXANG3 finds the angle (and its cosine) between two vectors | |
761 | C PXMDPI moves its argument onto the range [-pi,pi) | |
762 | C----------------------------------------------------------------------- | |
763 | SUBROUTINE PXZERV(N,A) | |
764 | INTEGER I,N | |
765 | DOUBLE PRECISION A(N) | |
766 | DO 10 I=1,N | |
767 | A(I)=0 | |
768 | 10 CONTINUE | |
769 | END | |
770 | *-- Author : | |
771 | C----------------------------------------------------------------------- | |
772 | SUBROUTINE PXADDV(N,A,B,C,ITERR) | |
773 | INTEGER I,N,ITERR | |
774 | DOUBLE PRECISION A(N),B(N),C(N) | |
775 | DO 10 I=1,N | |
776 | C(I)=A(I)+B(I) | |
777 | 10 CONTINUE | |
778 | END | |
779 | *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper | |
780 | *-- Author : | |
781 | C----------------------------------------------------------------------- | |
782 | SUBROUTINE PXANG3(A,B,COST,THET,ITERR) | |
783 | INTEGER ITERR | |
784 | DOUBLE PRECISION A(3),B(3),C,COST,THET | |
785 | C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2) | |
786 | IF (C.LE.0) RETURN | |
787 | C=1/SQRT(C) | |
788 | COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C | |
789 | THET=ACOS(COST) | |
790 | END | |
791 | *CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper | |
792 | *-- Author : P. Schleper 28/02/94 | |
793 | LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET) | |
794 | * | |
795 | INTEGER MXTRAK,MXPROT | |
796 | PARAMETER (MXTRAK=5000,MXPROT=5000) | |
797 | INTEGER NTRAK,NJET | |
798 | *** Note that although JETLIS is assumed to be a 2d array, it | |
799 | *** it is used as 1d in this routine for efficiency | |
800 | LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK) | |
801 | *** Checks to see if TSTLIS entries correspond to a jet already found | |
802 | *** and entered in JETLIS | |
803 | INTEGER N, I, IN | |
804 | LOGICAL MATCH | |
805 | * | |
806 | PXNEW = .TRUE. | |
807 | DO 100 I = 1,NJET | |
808 | MATCH = .TRUE. | |
809 | IN=I-MXPROT | |
810 | DO 110 N = 1,NTRAK | |
811 | IN=IN+MXPROT | |
812 | IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN | |
813 | MATCH = .FALSE. | |
814 | GO TO 100 | |
815 | ENDIF | |
816 | 110 CONTINUE | |
817 | IF (MATCH) THEN | |
818 | PXNEW = .FALSE. | |
819 | RETURN | |
820 | ENDIF | |
821 | 100 CONTINUE | |
822 | RETURN | |
823 | END | |
824 | *CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper | |
825 | *-- Author : P. Schleper 28/02/94 | |
826 | LOGICAL FUNCTION PXSAME(LIST1,LIST2,N) | |
827 | * | |
828 | LOGICAL LIST1(*),LIST2(*) | |
829 | INTEGER N | |
830 | *** Returns T if the first N elements of LIST1 are the same as the | |
831 | *** first N elements of LIST2. | |
832 | INTEGER I | |
833 | * | |
834 | PXSAME = .TRUE. | |
835 | DO 100 I = 1,N | |
836 | IF ( LIST1(I).NEQV.LIST2(I) ) THEN | |
837 | PXSAME = .FALSE. | |
838 | RETURN | |
839 | ENDIF | |
840 | 100 CONTINUE | |
841 | RETURN | |
842 | END | |
843 | *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper | |
844 | *-- Author : | |
845 | C----------------------------------------------------------------------- | |
846 | FUNCTION PXMDPI(PHI) | |
847 | IMPLICIT NONE | |
848 | C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) | |
849 | DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS | |
850 | PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961) | |
851 | PARAMETER (EPS=1E-15) | |
852 | PXMDPI=PHI | |
853 | IF (PXMDPI.LE.PI) THEN | |
854 | IF (PXMDPI.GT.-PI) THEN | |
855 | GOTO 100 | |
856 | ELSEIF (PXMDPI.GT.-THRPI) THEN | |
857 | PXMDPI=PXMDPI+TWOPI | |
858 | ELSE | |
859 | PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI | |
860 | ENDIF | |
861 | ELSEIF (PXMDPI.LE.THRPI) THEN | |
862 | PXMDPI=PXMDPI-TWOPI | |
863 | ELSE | |
864 | PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI | |
865 | ENDIF | |
866 | 100 IF (ABS(PXMDPI).LT.EPS) PXMDPI=0 | |
867 | END |