+ SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,
+ + MXJET,NJET,PJET,IPASS,IJMUL,IERR)
+*.*********************************************************
+*. ------
+*. PXCONE
+*. ------
+*.
+*.********** Pre Release Version 26.2.93
+*.
+*. Driver for the Cone Jet finding algorithm of L.A. del Pozo.
+*. Based on algorithm from D.E. Soper.
+*. Finds jets inside cone of half angle CONER with energy > EPSLON.
+*. Jets which receive more than a fraction OVLIM of their energy from
+*. overlaps with other jets are excluded.
+*. Output jets are ordered in energy.
+*. If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt)
+*. Usage :
+*.
+*. INTEGER ITKDM,MXTRK
+*. PARAMETER (ITKDM=4.or.more,MXTRK=1.or.more)
+*. INTEGER MXJET, MXTRAK, MXPROT
+*. PARAMETER (MXJET=10,MXTRAK=500,MXPROT=500)
+*. INTEGER IPASS (MXTRAK),IJMUL (MXJET)
+*. INTEGER NTRAK,NJET,IERR,MODE
+*. DOUBLE PRECISION PTRAK (ITKDM,MXTRK),PJET (5,MXJET)
+*. DOUBLE PRECISION CONER, EPSLON, OVLIM
+*. NTRAK = 1.to.MXTRAK
+*. CONER = ...
+*. EPSLON = ...
+*. OVLIM = ...
+*. CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET,
+*. + NJET,PJET,IPASS,IJMUL,IERR)
+*.
+*. INPUT : MODE 1=>e+e-, 2=>hadron-hadron
+*. INPUT : NTRAK Number of particles
+*. INPUT : ITKDM First dimension of PTRAK array
+*. INPUT : PTRAK Array of particle 4-momenta (Px,Py,Pz,E)
+*. INPUT : CONER Cone size (half angle) in radians
+*. INPUT : EPSLON Minimum Jet energy (GeV)
+*. INPUT : OVLIM Maximum fraction of overlap energy in a jet
+*. INPUT : MXJET Maximum possible number of jets
+*. OUTPUT : NJET Number of jets found
+*. OUTPUT : PJET 5-vectors of jets
+*. OUTPUT : IPASS(k) Particle k belongs to jet number IPASS(k)
+*. IPASS = -1 if not assosciated to a jet
+*. OUTPUT : IJMUL(i) Jet i contains IJMUL(i) particles
+*. OUTPUT : IERR = 0 if all is OK ; = -1 otherwise
+*.
+*. CALLS : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP
+*. CALLED : User
+*.
+*. AUTHOR : L.A. del Pozo
+*. CREATED : 26-Feb-93
+*. LAST MOD : 2-Mar-93
+*.
+*. Modification Log.
+*. 2-Jan-97: M Wobisch - fix bug concerning COS2R in eta phi mode
+*. 4-Apr-93: M H Seymour - Change 2d arrays to 1d in PXTRY & PXNEW
+*. 2-Apr-93: M H Seymour - Major changes to add boost-invariant mode
+*. 1-Apr-93: M H Seymour - Increase all array sizes
+*. 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION
+*. 30-Mar-93: M H Seymour - Change OVLIM into an input parameter
+*. 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP
+*. 1-Mar-93: L A del Pozo - Remove Cern library routine calls
+*. 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon
+*.
+*.*********************************************************
+C+SEQ,DECLARE.
+*** External Arrays
+ INTEGER ITKDM,MXJET,NTRAK,NJET,IERR,MODE
+ INTEGER IPASS (NTRAK),IJMUL (MXJET)
+ DOUBLE PRECISION PTRAK (ITKDM,NTRAK),PJET (5,MXJET),
+ + CONER, EPSLON, OVLIM
+*** Internal Arrays
+ INTEGER MXPROT, MXTRAK
+ PARAMETER (MXPROT=5000, MXTRAK=5000)
+ DOUBLE PRECISION PP(4,MXTRAK), PU(3,MXTRAK), PJ(4,MXPROT)
+ LOGICAL JETLIS(MXPROT,MXTRAK)
+*** Used in the routine.
+ DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ,
+ + COSVAL,PXMDPI
+cMWobisch
+ DOUBLE PRECISION RSEP
+cMWobisch
+ LOGICAL UNSTBL
+ INTEGER I,J,N,MU,N1,N2, ITERR
+ INTEGER NCALL, NPRINT
+ DOUBLE PRECISION ROLD, EPSOLD, OVOLD
+ SAVE NCALL,NPRINT,ROLD, EPSOLD, OVOLD
+ DATA NCALL,NPRINT /0,0/
+ DATA ROLD,EPSOLD,OVOLD/0.,0.,0./
+
+cMWobisch
+c***************************************
+ RSEP = 2D0
+c***************************************
+cMWobisch
+ IERR=0
+*
+*** INITIALIZE
+ IF(NCALL.LE.0) THEN
+ ROLD = 0.
+ EPSOLD = 0.
+ OVOLD = 0.
+ ENDIF
+ NCALL = NCALL + 1
+*
+*** Print welcome and Jetfinder parameters
+ IF((CONER.NE.ROLD .OR. EPSLON.NE.EPSOLD .OR. OVLIM.NE.OVOLD)
+ + .AND. NPRINT.LE.10) THEN
+ WRITE (6,*)
+ WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********'
+ WRITE (6,*) ' Written by Luis Del Pozo of OPAL'
+ WRITE (6,*) ' Modified for eta-phi by Mike Seymour'
+ WRITE(6,1000)' Cone Size R = ',CONER,' Radians'
+ WRITE(6,1001)' Min Jet energy Epsilon = ',EPSLON,' GeV'
+ WRITE(6,1002)' Overlap fraction parameter = ',OVLIM
+ WRITE (6,*) ' ***********************************************'
+cMWobisch
+ IF (RSEP .lt. 1.999) THEN
+ WRITE(6,*) ' '
+ WRITE (6,*) ' ******************************************'
+ WRITE (6,*) ' ******************************************'
+ WRITE(6,*) ' M Wobisch: private change !!!!!!!!!!!! '
+ WRITE(6,*) ' Rsep is set to ',RSEP
+ WRITE(6,*) ' this is ONLY meaningful in a NLO calculation'
+ WRITE(6,*) ' ------------------------ '
+ WRITE(6,*) ' please check what you''re doing!!'
+ WRITE(6,*) ' or ask: Markus.Wobisch@desy.de --'
+ WRITE (6,*) ' ******************************************'
+ WRITE (6,*) ' ******************************************'
+ WRITE (6,*) ' ******************************************'
+ WRITE(6,*) ' '
+ WRITE(6,*) ' '
+ ENDIF
+cMWobisch
+
+ WRITE (6,*)
+1000 FORMAT(A18,F5.2,A10)
+1001 FORMAT(A29,F5.2,A5)
+1002 FORMAT(A33,F5.2)
+ NPRINT = NPRINT + 1
+ ROLD=CONER
+ EPSOLD=EPSLON
+ OVOLD=OVLIM
+ ENDIF
+*
+*** Copy calling array PTRAK to internal array PP(4,NTRAK)
+*
+ IF (NTRAK .GT. MXTRAK) THEN
+ WRITE (6,*) ' PXCONE: Ntrak too large: ',NTRAK
+ IERR=-1
+ RETURN
+ ENDIF
+ IF (MODE.NE.2) THEN
+ DO 100 I=1, NTRAK
+ DO 101 J=1,4
+ PP(J,I)=PTRAK(J,I)
+101 CONTINUE
+100 CONTINUE
+ ELSE
+*** Converting to eta,phi,pt if necessary
+ DO 104 I=1,NTRAK
+ PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2
+ PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2
+ IF (PTSQ.LE.4.25E-18*PPSQ) THEN
+ PP(1,I)=20
+ ELSE
+ PP(1,I)=0.5*LOG(PPSQ/PTSQ)
+ ENDIF
+ PP(1,I)=SIGN(PP(1,I),PTRAK(3,I))
+ IF (PTSQ.EQ.0) THEN
+ PP(2,I)=0
+ ELSE
+ PP(2,I)=ATAN2(PTRAK(2,I),PTRAK(1,I))
+ ENDIF
+ PP(3,I)=0
+ PP(4,I)=SQRT(PTSQ)
+ PU(1,I)=PP(1,I)
+ PU(2,I)=PP(2,I)
+ PU(3,I)=PP(3,I)
+104 CONTINUE
+ ENDIF
+*
+*** Zero output variables
+*
+ NJET=0
+ DO 102 I = 1, NTRAK
+ DO 103 J = 1, MXPROT
+ JETLIS(J,I) = .FALSE.
+103 CONTINUE
+102 CONTINUE
+ CALL PXZERV(4*MXPROT,PJ)
+ CALL PXZERI(MXJET,IJMUL)
+*
+ IF (MODE.NE.2) THEN
+ COSR = COS(CONER)
+ COS2R = COS(CONER)
+ ELSE
+*** Purely for convenience, work in terms of 1-R**2
+ COSR = 1-CONER**2
+cMW -- select Rsep: 1-(Rsep*CONER)**2
+ COS2R = 1-(RSEP*CONER)**2
+cORIGINAL COS2R = 1-(2*CONER)**2
+ ENDIF
+ UNSTBL = .FALSE.
+ IF (MODE.NE.2) THEN
+ CALL PXUVEC(NTRAK,PP,PU,IERR)
+ IF (IERR .NE. 0) RETURN
+ ENDIF
+*** Look for jets using particle diretions as seed axes
+*
+ DO 110 N = 1,NTRAK
+ DO 120 MU = 1,3
+ VSEED(MU) = PU(MU,N)
+120 CONTINUE
+ CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,
+ & NJET,JETLIS,PJ,UNSTBL,IERR)
+ IF (IERR .NE. 0) RETURN
+110 CONTINUE
+
+cMW - for Rsep=1 goto 145
+c GOTO 145
+
+*** Now look between all pairs of jets as seed axes.
+ DO 140 N1 = 1,NJET-1
+ VEC1(1)=PJ(1,N1)
+ VEC1(2)=PJ(2,N1)
+ VEC1(3)=PJ(3,N1)
+ IF (MODE.NE.2) CALL PXNORV(3,VEC1,VEC1,ITERR)
+ DO 150 N2 = N1+1,NJET
+ VEC2(1)=PJ(1,N2)
+ VEC2(2)=PJ(2,N2)
+ VEC2(3)=PJ(3,N2)
+ IF (MODE.NE.2) CALL PXNORV(3,VEC2,VEC2,ITERR)
+ CALL PXADDV(3,VEC1,VEC2,VSEED,ITERR)
+ IF (MODE.NE.2) THEN
+ CALL PXNORV(3,VSEED,VSEED,ITERR)
+ ELSE
+ VSEED(1)=VSEED(1)/2
+ VSEED(2)=VSEED(2)/2
+ ENDIF
+C---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART
+ IF (MODE.NE.2) THEN
+ COSVAL=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
+ ELSE
+ IF (ABS(VEC1(1)).GE.20.OR.ABS(VEC2(1)).GE.20) THEN
+ COSVAL=-1000
+ ELSE
+ COSVAL=1-
+ + ((VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2)
+ ENDIF
+ ENDIF
+
+ IF (COSVAL.LE.COSR.AND.COSVAL.GE.COS2R)
+ + CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
+ + JETLIS,PJ,UNSTBL,IERR)
+c CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
+c + JETLIS,PJ,UNSTBL,IERR)
+ IF (IERR .NE. 0) RETURN
+150 CONTINUE
+140 CONTINUE
+ IF (UNSTBL) THEN
+ IERR=-1
+ WRITE (6,*) ' PXCONE: Too many iterations to find a proto-jet'
+ RETURN
+ ENDIF
+
+ 145 CONTINUE
+*** Now put the jet list into order by jet energy, eliminating jets
+*** with energy less than EPSLON.
+ CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
+*
+*** Take care of jet overlaps
+ CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
+*
+*** Order jets again as some have been eliminated, or lost energy.
+ CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
+*
+*** All done!, Copy output into output arrays
+ IF (NJET .GT. MXJET) THEN
+ WRITE (6,*) ' PXCONE: Found more than MXJET jets'
+ IERR=-1
+ GOTO 99
+ ENDIF
+ IF (MODE.NE.2) THEN
+ DO 300 I=1, NJET
+ DO 310 J=1,4
+ PJET(J,I)=PJ(J,I)
+310 CONTINUE
+300 CONTINUE
+ ELSE
+ DO 315 I=1, NJET
+ PJET(1,I)=PJ(4,I)*COS(PJ(2,I))
+ PJET(2,I)=PJ(4,I)*SIN(PJ(2,I))
+ PJET(3,I)=PJ(4,I)*SINH(PJ(1,I))
+ PJET(4,I)=PJ(4,I)*COSH(PJ(1,I))
+ 315 CONTINUE
+ ENDIF
+ DO 320 I=1, NTRAK
+ IPASS(I)=-1
+ DO 330 J=1, NJET
+ IF (JETLIS(J,I)) THEN
+ IJMUL(J)=IJMUL(J)+1
+ IPASS(I)=J
+ ENDIF
+330 CONTINUE
+320 CONTINUE
+99 RETURN
+ END
+*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+ SUBROUTINE PXNORV(N,A,B,ITERR)
+ INTEGER I,N,ITERR
+ DOUBLE PRECISION A(N),B(N),C
+ C=0
+ DO 10 I=1,N
+ C=C+A(I)**2
+ 10 CONTINUE
+ IF (C.LE.0) RETURN
+ C=1/SQRT(C)
+ DO 20 I=1,N
+ B(I)=A(I)*C
+ 20 CONTINUE
+ END
+*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
+*CMZ : 1.06/00 15/03/94 12.17.46 by P. Schleper
+*-- Author :
+*
+C+DECK,PXOLAP.
+ SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
+*
+*** Looks for particles assigned to more than 1 jet, and reassigns them
+*** If more than a fraction OVLIM of a jet's energy is contained in
+*** higher energy jets, that jet is neglected.
+*** Particles assigned to the jet closest in angle (a la CDF, Snowmass).
+C+SEQ,DECLARE.
+ INTEGER MXTRAK, MXPROT
+ PARAMETER (MXTRAK=5000,MXPROT=5000)
+ INTEGER NJET, NTRAK, MODE
+ LOGICAL JETLIS(MXPROT,MXTRAK)
+ DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI
+ INTEGER I,J,N,MU
+ LOGICAL OVELAP
+ DOUBLE PRECISION EOVER
+ DOUBLE PRECISION OVLIM
+ INTEGER ITERR, IJMIN, IJET(MXPROT), NJ
+ DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN
+ DATA IJMIN/0/
+*
+ IF (NJET.LE.1) RETURN
+*** Look for jets with large overlaps with higher energy jets.
+ DO 100 I = 2,NJET
+*** Find overlap energy between jets I and all higher energy jets.
+ EOVER = 0.0
+ DO 110 N = 1,NTRAK
+ OVELAP = .FALSE.
+ DO 120 J= 1,I-1
+ IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN
+ OVELAP = .TRUE.
+ ENDIF
+120 CONTINUE
+ IF (OVELAP) THEN
+ EOVER = EOVER + PP(4,N)
+ ENDIF
+110 CONTINUE
+*** Is the fraction of energy shared larger than OVLIM?
+ IF (EOVER.GT.OVLIM*PJ(4,I)) THEN
+*** De-assign all particles from Jet I
+ DO 130 N = 1,NTRAK
+ JETLIS(I,N) = .FALSE.
+130 CONTINUE
+ ENDIF
+100 CONTINUE
+*** Now there are no big overlaps, assign every particle in
+*** more than 1 jet to the closet jet.
+*** Any particles now in more than 1 jet are assigned to the CLOSET
+*** jet (in angle).
+ DO 140 I=1,NTRAK
+ NJ=0
+ DO 150 J=1, NJET
+ IF(JETLIS(J,I)) THEN
+ NJ=NJ+1
+ IJET(NJ)=J
+ ENDIF
+150 CONTINUE
+ IF (NJ .GT. 1) THEN
+*** Particle in > 1 jet - calc angles...
+ VEC1(1)=PP(1,I)
+ VEC1(2)=PP(2,I)
+ VEC1(3)=PP(3,I)
+ THMIN=0.
+ DO 160 J=1,NJ
+ VEC2(1)=PJ(1,IJET(J))
+ VEC2(2)=PJ(2,IJET(J))
+ VEC2(3)=PJ(3,IJET(J))
+ IF (MODE.NE.2) THEN
+ CALL PXANG3(VEC1,VEC2,COST,THET,ITERR)
+ ELSE
+ THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2
+ ENDIF
+ IF (J .EQ. 1) THEN
+ THMIN=THET
+ IJMIN=IJET(J)
+ ELSEIF (THET .LT. THMIN) THEN
+ THMIN=THET
+ IJMIN=IJET(J)
+ ENDIF
+160 CONTINUE
+*** Assign track to IJMIN
+ DO 170 J=1,NJET
+ JETLIS(J,I) = .FALSE.
+170 CONTINUE
+ JETLIS(IJMIN,I)=.TRUE.
+ ENDIF
+140 CONTINUE
+*** Recompute PJ
+ DO 200 I = 1,NJET
+ DO 210 MU = 1,4
+ PJ(MU,I) = 0.0
+210 CONTINUE
+ DO 220 N = 1,NTRAK
+ IF( JETLIS(I,N) ) THEN
+ IF (MODE.NE.2) THEN
+ DO 230 MU = 1,4
+ PJ(MU,I) = PJ(MU,I) + PP(MU,N)
+230 CONTINUE
+ ELSE
+ PJ(1,I)=PJ(1,I)
+ + + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I))
+ PJ(2,I)=PJ(2,I)
+ + + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I))
+ PJ(4,I)=PJ(4,I)+PP(4,N)
+ ENDIF
+ ENDIF
+220 CONTINUE
+200 CONTINUE
+ RETURN
+ END
+*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
+*CMZ : 1.06/00 14/03/94 15.37.45 by P. Schleper
+*-- Author :
+*
+C+DECK,PXORD.
+ SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
+*
+*** Routine to put jets into order and eliminate tose less than EPSLON
+C+SEQ,DECLARE.
+ INTEGER MXTRAK,MXPROT
+ PARAMETER (MXTRAK=5000,MXPROT=5000)
+ INTEGER I, J, INDEX(MXPROT)
+ DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT)
+ INTEGER NJET,NTRAK
+ LOGICAL JETLIS(MXPROT,MXTRAK)
+ LOGICAL LOGTMP(MXPROT,MXTRAK)
+ DOUBLE PRECISION EPSLON,PJ(4,MXPROT)
+*** Puts jets in order of energy: 1 = highest energy etc.
+*** Then Eliminate jets with energy below EPSLON
+*
+*** Copy input arrays.
+ DO 100 I=1,NJET
+ DO 110 J=1,4
+ PTEMP(J,I)=PJ(J,I)
+110 CONTINUE
+ DO 120 J=1,NTRAK
+ LOGTMP(I,J)=JETLIS(I,J)
+120 CONTINUE
+100 CONTINUE
+ DO 150 I=1,NJET
+ ELIST(I)=PJ(4,I)
+150 CONTINUE
+*** Sort the energies...
+ CALL PXSORV(NJET,ELIST,INDEX,'I')
+*** Fill PJ and JETLIS according to sort ( sort is in ascending order!!)
+ DO 200 I=1, NJET
+ DO 210 J=1,4
+ PJ(J,I)=PTEMP(J,INDEX(NJET+1-I))
+210 CONTINUE
+ DO 220 J=1,NTRAK
+ JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J)
+220 CONTINUE
+200 CONTINUE
+** Jets are now in order
+*** Now eliminate jets with less than Epsilon energy
+ DO 300, I=1, NJET
+ IF (PJ(4,I) .LT. EPSLON) THEN
+ NJET=NJET-1
+ PJ(4,I)=0.
+ ENDIF
+300 CONTINUE
+ RETURN
+ END
+
+********************************************************************
+*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
+*CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper
+*-- Author :
+C+DECK,PXSEAR.
+ SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
+ + JETLIS,PJ,UNSTBL,IERR)
+*
+C+SEQ,DECLARE.
+ INTEGER MXTRAK, MXPROT
+ PARAMETER (MXTRAK=5000,MXPROT=5000)
+ INTEGER NTRAK, IERR, MODE
+ DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3)
+ LOGICAL UNSTBL
+ LOGICAL JETLIS(MXPROT,MXTRAK)
+ INTEGER NJET
+ DOUBLE PRECISION PJ(4,MXPROT)
+*** Using VSEED as a trial axis , look for a stable jet.
+*** Check stable jets against those already found and add to PJ.
+*** Will try up to MXITER iterations to get a stable set of particles
+*** in the cone.
+ INTEGER MU,N,ITER
+ LOGICAL PXSAME,PXNEW,OK
+ LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK)
+ DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4)
+ INTEGER MXITER
+ PARAMETER(MXITER = 30)
+*
+ DO 100 MU=1,3
+ OAXIS(MU) = VSEED(MU)
+100 CONTINUE
+ DO 110 N = 1,NTRAK
+ OLDLIS(N) = .FALSE.
+110 CONTINUE
+ DO 120 ITER = 1,MXITER
+ CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK)
+*** Return immediately if there were no particles in the cone.
+ IF (.NOT.OK) THEN
+ RETURN
+ ENDIF
+ IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN
+*** We have a stable jet.
+ IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN
+*** And the jet is a new one. So add it to our arrays.
+*** Check arrays are big anough...
+ IF (NJET .EQ. MXPROT) THEN
+ WRITE (6,*) ' PXCONE: Found more than MXPROT proto-jets'
+ IERR = -1
+ RETURN
+ ENDIF
+ NJET = NJET + 1
+ DO 130 N = 1,NTRAK
+ JETLIS(NJET,N) = NEWLIS(N)
+130 CONTINUE
+ DO 140 MU=1,4
+ PJ(MU,NJET)=PNEW(MU)
+140 CONTINUE
+ ENDIF
+ RETURN
+ ENDIF
+*** The jet was not stable, so we iterate again
+ DO 150 N=1,NTRAK
+ OLDLIS(N)=NEWLIS(N)
+150 CONTINUE
+ DO 160 MU=1,3
+ OAXIS(MU)=NAXIS(MU)
+160 CONTINUE
+120 CONTINUE
+ UNSTBL = .TRUE.
+ RETURN
+ END
+*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+ SUBROUTINE PXSORV(N,A,K,OPT)
+C Sort A(N) into ascending order
+C OPT = 'I' : return index array K only
+C OTHERWISE : return sorted A and index array K
+C-----------------------------------------------------------------------
+ INTEGER NMAX
+ PARAMETER (NMAX=5000)
+*
+* INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX)
+*LUND
+ INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX)
+ CHARACTER OPT
+*
+* DOUBLE PRECISION A(N),B(NMAX)
+ DOUBLE PRECISION A(NMAX),B(NMAX)
+*LUND
+ IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV'
+ IL(1)=0
+ IR(1)=0
+ DO 10 I=2,N
+ IL(I)=0
+ IR(I)=0
+ J=1
+ 2 IF(A(I).GT.A(J)) GO TO 5
+ 3 IF(IL(J).EQ.0) GO TO 4
+ J=IL(J)
+ GO TO 2
+ 4 IR(I)=-J
+ IL(J)=I
+ GO TO 10
+ 5 IF(IR(J).LE.0) GO TO 6
+ J=IR(J)
+ GO TO 2
+ 6 IR(I)=IR(J)
+ IR(J)=I
+ 10 CONTINUE
+ I=1
+ J=1
+ GO TO 8
+ 20 J=IL(J)
+ 8 IF(IL(J).GT.0) GO TO 20
+ 9 K(I)=J
+ B(I)=A(J)
+ I=I+1
+ IF(IR(J)) 12,30,13
+ 13 J=IR(J)
+ GO TO 8
+ 12 J=-IR(J)
+ GO TO 9
+ 30 IF(OPT.EQ.'I') RETURN
+ DO 31 I=1,N
+ 31 A(I)=B(I)
+ 999 END
+
+*********************************************************************
+*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
+*CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper
+*-- Author :
+*
+C+DECK,PXTRY.
+ SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,
+ + PNEW,NEWLIS,OK)
+*
+C+SEQ,DECLARE.
+ INTEGER MXTRAK
+ PARAMETER (MXTRAK=5000)
+ INTEGER NTRAK,MODE
+*** Note that although PU and PP are assumed to be 2d arrays, they
+*** are used as 1d in this routine for efficiency
+ DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI
+ LOGICAL OK
+ LOGICAL NEWLIS(MXTRAK)
+ DOUBLE PRECISION NAXIS(3),PNEW(4)
+*** Finds all particles in cone of size COSR about OAXIS direction.
+*** Calculates 4-momentum sum of all particles in cone (PNEW) , and
+*** returns this as new jet axis NAXIS (Both unit Vectors)
+ INTEGER N,MU,NPU,NPP
+ DOUBLE PRECISION COSVAL,NORMSQ,NORM
+*
+ OK = .FALSE.
+ DO 100 MU=1,4
+ PNEW(MU)=0.0
+100 CONTINUE
+ NPU=-3
+ NPP=-4
+ DO 110 N=1,NTRAK
+ NPU=NPU+3
+ NPP=NPP+4
+ IF (MODE.NE.2) THEN
+ COSVAL=0.0
+ DO 120 MU=1,3
+ COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU)
+120 CONTINUE
+ ELSE
+ IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN
+ COSVAL=-1000
+ ELSE
+ COSVAL=1-
+ + ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2)
+ ENDIF
+ ENDIF
+ IF (COSVAL.GE.COSR)THEN
+ NEWLIS(N) = .TRUE.
+ OK = .TRUE.
+ IF (MODE.NE.2) THEN
+ DO 130 MU=1,4
+ PNEW(MU) = PNEW(MU) + PP(MU+NPP)
+130 CONTINUE
+ ELSE
+ PNEW(1)=PNEW(1)
+ + + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1))
+ PNEW(2)=PNEW(2)
+ + + PP(4+NPP)/(PP(4+NPP)+PNEW(4))
+ + *PXMDPI(PP(2+NPP)-PNEW(2))
+ PNEW(4)=PNEW(4)+PP(4+NPP)
+ ENDIF
+ ELSE
+ NEWLIS(N)=.FALSE.
+ ENDIF
+110 CONTINUE
+*** If there are particles in the cone, calc new jet axis
+ IF (OK) THEN
+ IF (MODE.NE.2) THEN
+ NORMSQ = 0.0
+ DO 140 MU = 1,3
+ NORMSQ = NORMSQ + PNEW(MU)**2
+140 CONTINUE
+ NORM = SQRT(NORMSQ)
+ ELSE
+ NORM = 1
+ ENDIF
+ DO 150 MU=1,3
+ NAXIS(MU) = PNEW(MU)/NORM
+150 CONTINUE
+ ENDIF
+ RETURN
+ END
+
+*********************************************************************
+*CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
+*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
+*-- Author :
+C+DECK,PXUVEC.
+*
+ SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR)
+*
+*** Routine to calculate unit vectors PU of all particles PP
+C+SEQ,DECLARE.
+ INTEGER MXTRAK
+ PARAMETER (MXTRAK=5000)
+ INTEGER NTRAK, IERR
+ DOUBLE PRECISION PP(4,MXTRAK)
+ DOUBLE PRECISION PU(3,MXTRAK)
+ INTEGER N,MU
+ DOUBLE PRECISION MAG
+ DO 100 N=1,NTRAK
+ MAG=0.0
+ DO 110 MU=1,3
+ MAG=MAG+PP(MU,N)**2
+110 CONTINUE
+ MAG=SQRT(MAG)
+ IF (MAG.EQ.0.0) THEN
+ WRITE(6,*)' PXCONE: An input particle has zero mod(p)'
+ IERR=-1
+ RETURN
+ ENDIF
+ DO 120 MU=1,3
+ PU(MU,N)=PP(MU,N)/MAG
+120 CONTINUE
+100 CONTINUE
+ RETURN
+ END
+*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+ SUBROUTINE PXZERI(N,A)
+ INTEGER I,N,A(N)
+ DO 10 I=1,N
+ A(I)=0
+ 10 CONTINUE
+ END
+*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+C This is a set of routines written by Mike Seymour to provide the
+C services presumably normally provided by standard OPAL routines
+C PXZERV zeroes a vector
+C PXZERI zeroes a vector of integers
+C PXNORV normalizes a vector
+C PXADDV adds two vectors
+C PXSORV sorts a vector (copied from HERWIG)
+C PXANG3 finds the angle (and its cosine) between two vectors
+C PXMDPI moves its argument onto the range [-pi,pi)
+C-----------------------------------------------------------------------
+ SUBROUTINE PXZERV(N,A)
+ INTEGER I,N
+ DOUBLE PRECISION A(N)
+ DO 10 I=1,N
+ A(I)=0
+ 10 CONTINUE
+ END
+*-- Author :
+C-----------------------------------------------------------------------
+ SUBROUTINE PXADDV(N,A,B,C,ITERR)
+ INTEGER I,N,ITERR
+ DOUBLE PRECISION A(N),B(N),C(N)
+ DO 10 I=1,N
+ C(I)=A(I)+B(I)
+ 10 CONTINUE
+ END
+*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+ SUBROUTINE PXANG3(A,B,COST,THET,ITERR)
+ INTEGER ITERR
+ DOUBLE PRECISION A(3),B(3),C,COST,THET
+ C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2)
+ IF (C.LE.0) RETURN
+ C=1/SQRT(C)
+ COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C
+ THET=ACOS(COST)
+ END
+*CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper
+*-- Author : P. Schleper 28/02/94
+ LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET)
+*
+ INTEGER MXTRAK,MXPROT
+ PARAMETER (MXTRAK=5000,MXPROT=5000)
+ INTEGER NTRAK,NJET
+*** Note that although JETLIS is assumed to be a 2d array, it
+*** it is used as 1d in this routine for efficiency
+ LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK)
+*** Checks to see if TSTLIS entries correspond to a jet already found
+*** and entered in JETLIS
+ INTEGER N, I, IN
+ LOGICAL MATCH
+*
+ PXNEW = .TRUE.
+ DO 100 I = 1,NJET
+ MATCH = .TRUE.
+ IN=I-MXPROT
+ DO 110 N = 1,NTRAK
+ IN=IN+MXPROT
+ IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN
+ MATCH = .FALSE.
+ GO TO 100
+ ENDIF
+110 CONTINUE
+ IF (MATCH) THEN
+ PXNEW = .FALSE.
+ RETURN
+ ENDIF
+100 CONTINUE
+ RETURN
+ END
+*CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper
+*-- Author : P. Schleper 28/02/94
+ LOGICAL FUNCTION PXSAME(LIST1,LIST2,N)
+*
+ LOGICAL LIST1(*),LIST2(*)
+ INTEGER N
+*** Returns T if the first N elements of LIST1 are the same as the
+*** first N elements of LIST2.
+ INTEGER I
+*
+ PXSAME = .TRUE.
+ DO 100 I = 1,N
+ IF ( LIST1(I).NEQV.LIST2(I) ) THEN
+ PXSAME = .FALSE.
+ RETURN
+ ENDIF
+100 CONTINUE
+ RETURN
+ END
+*CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
+*-- Author :
+C-----------------------------------------------------------------------
+ FUNCTION PXMDPI(PHI)
+ IMPLICIT NONE
+C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI)
+ DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS
+ PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961)
+ PARAMETER (EPS=1E-15)
+ PXMDPI=PHI
+ IF (PXMDPI.LE.PI) THEN
+ IF (PXMDPI.GT.-PI) THEN
+ GOTO 100
+ ELSEIF (PXMDPI.GT.-THRPI) THEN
+ PXMDPI=PXMDPI+TWOPI
+ ELSE
+ PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI
+ ENDIF
+ ELSEIF (PXMDPI.LE.THRPI) THEN
+ PXMDPI=PXMDPI-TWOPI
+ ELSE
+ PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI
+ ENDIF
+ 100 IF (ABS(PXMDPI).LT.EPS) PXMDPI=0
+ END