pxcone.f

      SUBROUTINE pxcone(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,
     +                   mxjet,njet,pjet,ipass,ijmul,ierr)
*.*********************************************************
*. ------
*. PXCONE
*. ------
*.
*. Code downloaded from the following web page 
*.
*.   http://aliceinfo.cern.ch/alicvs/viewvc/JETAN/pxcone.F?view=markup&pathrev=v4-05-04
*. 
*. on 17/10/2006 by G. Salam. Permission subsequently granted by Michael
*. H. Seymour (on behalf of the PxCone authors) for this code to be
*. distributed together with FastJet under the terms of the GNU Public
*. License v2 (see the file COPYING in the main FastJet directory).
*.
*.********** 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.
*. 25-Feb-07: G P Salam   - fix bugs concerning 2pi periodicity in eta phi mode
*.                        - added commented code to get consistent behaviour
*.                          regardless of particle order (replaces n-way 
*.                          midpoints with 2-way midpoints however...)
*. 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, njtorg
      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,*) '    Includes bug fixes by Wobisch, Salam'
         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,*) '    PXCONE is not a supported product and is'
         WRITE (6,*) '    is provided for comparative purposes only'
         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.
c      NJTORG = NJET           ! GPS -- to get consistent behaviour (2-way midpnts)
c      DO 140 N1 = 1,NJTORG-1  ! GPS -- to get consistent behaviour (2-way midpnts)
      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)
C         DO 150 N2 = N1+1,NJTORG ! GPS -- to get consistent behaviour
         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
               ! GPS 25/02/07
               vseed(2)=pxmdpi(vec1(2)+0.5d0*pxmdpi(vec2(2)-vec1(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))
c GPS 25/02/07
                pj(2,i)=pxmdpi(pj(2,i)
     +             + pp(4,n)/(pp(4,n)+pj(4,i))*pxmdpi(pp(2,n)-pj(2,i)))
c                PJ(2,I)=PJ(2,I)
c     +               + 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))
c                PNEW(2)=PNEW(2)
c     +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))
c     +               *PXMDPI(PP(2+NPP)-PNEW(2))
! GPS 25/02/07
                pnew(2)=pxmdpi(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