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