o First Version of TRDnSigma implementation (Xianguo) o still requires some catching...
[u/mrichter/AliRoot.git] / JETAN / pxcone.F
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