git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
First commit of new jet reconstruction and analysis code to be used for the
[u/mrichter/AliRoot.git] / JETAN / pxcone.F
CommitLineData
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*.*********************************************************
68C+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
82cMWobisch
83 DOUBLE PRECISION RSEP
84cMWobisch
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
93cMWobisch
94c***************************************
95 RSEP = 2D0
96c***************************************
97cMWobisch
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,*) ' ***********************************************'
119cMWobisch
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
136cMWobisch
137
138 WRITE (6,*)
1391000 FORMAT(A18,F5.2,A10)
1401001 FORMAT(A29,F5.2,A5)
1411002 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)
159101 CONTINUE
160100 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)
182104 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.
191103 CONTINUE
192102 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
202cMW -- select Rsep: 1-(Rsep*CONER)**2
203 COS2R = 1-(RSEP*CONER)**2
204cORIGINAL 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)
216120 CONTINUE
217 CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,
218 & NJET,JETLIS,PJ,UNSTBL,IERR)
219 IF (IERR .NE. 0) RETURN
220110 CONTINUE
221
222cMW - for Rsep=1 goto 145
223c 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
243C---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)
258c CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
259c + JETLIS,PJ,UNSTBL,IERR)
260 IF (IERR .NE. 0) RETURN
261150 CONTINUE
262140 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)
290310 CONTINUE
291300 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
307330 CONTINUE
308320 CONTINUE
30999 RETURN
310 END
311*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
312*-- Author :
313C-----------------------------------------------------------------------
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*
331C+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).
338C+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
363120 CONTINUE
364 IF (OVELAP) THEN
365 EOVER = EOVER + PP(4,N)
366 ENDIF
367110 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.
373130 CONTINUE
374 ENDIF
375100 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
387150 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
410160 CONTINUE
411*** Assign track to IJMIN
412 DO 170 J=1,NJET
413 JETLIS(J,I) = .FALSE.
414170 CONTINUE
415 JETLIS(IJMIN,I)=.TRUE.
416 ENDIF
417140 CONTINUE
418*** Recompute PJ
419 DO 200 I = 1,NJET
420 DO 210 MU = 1,4
421 PJ(MU,I) = 0.0
422210 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)
428230 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
437220 CONTINUE
438200 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*
445C+DECK,PXORD.
446 SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
447*
448*** Routine to put jets into order and eliminate tose less than EPSLON
449C+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)
465110 CONTINUE
466 DO 120 J=1,NTRAK
467 LOGTMP(I,J)=JETLIS(I,J)
468120 CONTINUE
469100 CONTINUE
470 DO 150 I=1,NJET
471 ELIST(I)=PJ(4,I)
472150 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))
479210 CONTINUE
480 DO 220 J=1,NTRAK
481 JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J)
482220 CONTINUE
483200 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
491300 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 :
499C+DECK,PXSEAR.
500 SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
501 + JETLIS,PJ,UNSTBL,IERR)
502*
503C+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)
525100 CONTINUE
526 DO 110 N = 1,NTRAK
527 OLDLIS(N) = .FALSE.
528110 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)
548130 CONTINUE
549 DO 140 MU=1,4
550 PJ(MU,NJET)=PNEW(MU)
551140 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)
558150 CONTINUE
559 DO 160 MU=1,3
560 OAXIS(MU)=NAXIS(MU)
561160 CONTINUE
562120 CONTINUE
563 UNSTBL = .TRUE.
564 RETURN
565 END
566*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
567*-- Author :
568C-----------------------------------------------------------------------
569 SUBROUTINE PXSORV(N,A,K,OPT)
570C Sort A(N) into ascending order
571C OPT = 'I' : return index array K only
572C OTHERWISE : return sorted A and index array K
573C-----------------------------------------------------------------------
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*
628C+DECK,PXTRY.
629 SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,
630 + PNEW,NEWLIS,OK)
631*
632C+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
651100 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)
661120 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)
676130 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
688110 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
695140 CONTINUE
696 NORM = SQRT(NORMSQ)
697 ELSE
698 NORM = 1
699 ENDIF
700 DO 150 MU=1,3
701 NAXIS(MU) = PNEW(MU)/NORM
702150 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 :
711C+DECK,PXUVEC.
712*
713 SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR)
714*
715*** Routine to calculate unit vectors PU of all particles PP
716C+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
728110 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
737120 CONTINUE
738100 CONTINUE
739 RETURN
740 END
741*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
742*-- Author :
743C-----------------------------------------------------------------------
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 :
752C-----------------------------------------------------------------------
753C This is a set of routines written by Mike Seymour to provide the
754C services presumably normally provided by standard OPAL routines
755C PXZERV zeroes a vector
756C PXZERI zeroes a vector of integers
757C PXNORV normalizes a vector
758C PXADDV adds two vectors
759C PXSORV sorts a vector (copied from HERWIG)
760C PXANG3 finds the angle (and its cosine) between two vectors
761C PXMDPI moves its argument onto the range [-pi,pi)
762C-----------------------------------------------------------------------
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 :
771C-----------------------------------------------------------------------
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 :
781C-----------------------------------------------------------------------
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
816110 CONTINUE
817 IF (MATCH) THEN
818 PXNEW = .FALSE.
819 RETURN
820 ENDIF
821100 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
840100 CONTINUE
841 RETURN
842 END
843*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
844*-- Author :
845C-----------------------------------------------------------------------
846 FUNCTION PXMDPI(PHI)
847 IMPLICIT NONE
848C---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