pxcone.f

Go to the documentation of this file.

CDECK  ID>, PXCONE.
      SUBROUTINE pxcone (NTRAK,ITKDM,PTRAK,CONER,EPSLON,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.
*. Output jets are ordered in energy.
*. Usage     :
*.
*.      INTEGER  ITKDM,MXTRK
*.      PARAMETER  (ITKDM=4.or.more,MXTRK=1.or.more)
*.      INTEGER  MXJET, MXTRAK, MXPROT
*.      PARAMETER  (MXJET=10,MXTRAK=200,MXPROT=100)
*.      INTEGER  IPASS (MXTRAK),IJMUL (MXJET)
*.      INTEGER  NTRAK,NJET,IERR
*.      REAL  PTRAK (ITKDM,MXTRK),PJET (5,MXJET)
*.      REAL  CONER, EPSLON
*.      NTRAK = 1.to.MXTRAK
*.      CONER   = 0.7   (suggested value)
*.      EPSLON  = 7.0   (suggested value)
*.      CALL PXCONE (NTRAK,ITKDM,PTRAK,CONER,EPSLON,MXJET,
*.     +             NJET,PJET,IPASS,IJMUL,IERR)
*.
*. 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     :  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  :  08-Oct-97
*.
*. Modification Log.
*. 08-Oct-97: D. Chrisman - Call PXADDV with the correct number
*.                          of arguments.
*. 26-Jun-96: D. Chrisman - Save ROLD and EPSOLD
*. 3-Mar-93: L A del Pozo - Check Cone size is sensible.
*. 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
*.
*.*********************************************************
      IMPLICIT NONE
*** External Arrays
      INTEGER  itkdm,mxjet,ntrak,njet,ierr
      INTEGER  ipass (*),ijmul (mxjet)
      REAL  ptrak (itkdm,*),pjet (5,*), coner, epslon
*** Internal Arrays
      INTEGER mxprot, mxtrak
      parameter(mxprot=100, mxtrak=200)
      REAL pp(4,mxtrak), pu(3,mxtrak), pj(4,mxprot)
      LOGICAL jetlis(mxprot,mxtrak)
*** Used in the routine.
      REAL cosr,cos2r, vseed(3), vec1(3), vec2(3)
      LOGICAL unstbl
      INTEGER i,j,n,mu,n1,n2, iterr
      REAL rold, epsold
      SAVE rold, epsold
      INTEGER ncall, nprint
      SAVE ncall,nprint
      DATA ncall,nprint /0,0/
      ierr=0
*
*** INITIALIZE
      IF(ncall.LE.0)  THEN
         rold = 0.
         epsold = 0.
      ENDIF
      ncall = ncall + 1
*
*** Print welcome and Jetfinder parameters
      IF((coner.NE.rold .OR. epslon.NE.epsold) .AND. nprint.LE.10) THEN
         WRITE (6,*)
         WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********'
         WRITE(6,1000)'   Cone Size R = ',coner,' Radians'
         WRITE(6,1001)'   Min Jet energy Epsilon = ',epslon,' GeV'
         WRITE (6,*) ' ***********************************************'
         WRITE (6,*)
1000     FORMAT(a18,f5.2,a10)
1001     FORMAT(a29,f5.2,a5)
         nprint = nprint + 1
         rold=coner
         epsold=epslon
      ENDIF
*** Check input Value of R is sensible.
      IF (coner .GT. 1.5708) THEN
         WRITE (6,*) ' PXCONE: CONER > 1.57 rad (90 degrees)'
         ierr=-1
         RETURN
      ENDIF
*
*** Copy calling array PTRAK  to internal array PP(4,NTRAK)
*
      IF (ntrak .GT. mxtrak) THEN
         WRITE (6,*) ' PXCONE: Ntrak too large'
         ierr=-1
         RETURN
      ENDIF
      DO  100 i=1, ntrak
         DO  101 j=1,4
            pp(j,i)=ptrak(j,i)
101      CONTINUE
100   CONTINUE
*
*** 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 pxzerv(mxjet,ijmul)
*
      cosr = cos(coner)
      cos2r = cos(2*coner)
      unstbl = .false.
      CALL pxuvec(ntrak,pp,pu,ierr)
      IF (ierr .NE. 0) RETURN
*** 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(cosr,ntrak,pu,pp,vseed,
     +                   njet,jetlis,pj,unstbl,ierr)
         IF (ierr .NE. 0) RETURN
110   CONTINUE
*** 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)
         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)
            CALL pxnorv(3,vec2,vec2,iterr)
            CALL pxaddv(3,vec1,vec2,vseed)
            CALL pxnorv(3,vseed,vseed,iterr)
            CALL pxsear(cosr,ntrak,pu,pp,vseed,njet,
     +      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
*** 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(njet,ntrak,jetlis,pj,pp)
*
*** 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
      DO 300 i=1, njet
         DO 310 j=1,4
            pjet(j,i)=pj(j,i)
310      CONTINUE
300   CONTINUE
      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