pxcamj.f

Go to the documentation of this file.

CDECK  ID>, PXCAMJ.
      SUBROUTINE pxcamj(ITKDM,NT,PT,YCUT,NJ,IJ,PJ,IERR)
C
C---CAMBRIDGE JET CLUSTERING ALGORITHM
C   BASED ON YCLUS BY S BETHKE
C   REF: YU L DOKSHITZER, G D LEDER, S MORETTI, B R WEBBER
C   CAVENDISH-HEP-97/06 (JUNE 1997)
C   07/07/97 FIRST RELEASED BY BRW
C   23/08/97 COMMENTS REVISED BY BRW
C
C   INPUT:
C    ITKDM = 1st dimension of array PT, ITKDM >= 4 required
C       NT = NUMBER OF TRACKS
C   PT(,I) = 4-MOMENTUM OF TRACK I (I=1,NT)
C     YCUT = (DURHAM) JET RESOLUTION
C
C  OUTPUT:
C       NJ = NUMBER OF JETS
C    IJ(I) = J IF TRACK I BELONGS TO JET J (I=1,NT)
C  PJ(4,J) = 4-MOMENTUM OF JET J (J=1,NJ)
C
C      NB:   CLUSTERING SEQUENCE DEPENDS ON VALUE OF YCUT
C
C  Modifications:
C  23.09.97, STK: Variable 1st dim for PT, REAL call args, handle
C                 errors with error flag IERR
C  24.09.97, STK: Improved error handling, sorted declarations,
C                 added welcome message
C  29.09.97, SB : Calculate EVIS from PT array; introduce to PX library
C-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER itkdm,nt,nj,ij(*),ierr
      REAL pt(itkdm,*),evis,ycut,pj(4,*)
      INTEGER ntrk,nv
      parameter( ntrk=300, nv=ntrk*(ntrk-2)+ntrk-(ntrk-2)*(ntrk-1)/2 )
      LOGICAL ip(ntrk),lcall
      INTEGER i,ii,j,k,l,imini,jmini,iad,jj(ntrk)
      DOUBLE PRECISION pp(5),pl(5,ntrk),v(nv),pm,vmini,ysca
      SAVE lcall
      DATA lcall / .false. /
C  Welcome message:
      IF( .NOT.lcall ) THEN
        print *, ' '
        print *, 'Cambridge jet finding algorithm, please refer to:'
        print *, 'Yu.L. Dokshitzer, G.D. Leder, S. Moretti, B.R. Webber'
        print *, 'CAVENDISH-HEP-97/06 (June 1997)'
        print *, ' '
        lcall= .true.
      ENDIF
C---WARNINGS
      ierr = 0
      evis = 0.
      IF( nt.GT.ntrk ) THEN
        WRITE(*,'(''CAMJET: More than '',I3,'' input particles: '',I5)')
     &        ntrk,nt
        ierr= 1
        RETURN
      ELSEIF( nt.LT.2 ) THEN
        WRITE(*,'(''CAMJET: Less than 2 input particles: '',I5)') nt
        ierr= 2
        RETURN
      ENDIF
C---COPY MOMENTA INTO PL-ARRAY
      DO i=1,nt
        ip(i)=.true.
        ij(i)=i
        DO ii=1,4
          pl(ii,i)= dble(pt(ii,i))
        ENDDO
        pm=pl(1,i)**2+pl(2,i)**2+pl(3,i)**2
        evis = evis + pt(4,i)
        IF (pm.GT.0d0) THEN
          pl(5,i)=1d0/sqrt(pm)
        ELSE
          pl(5,i)=1d0
        ENDIF
      ENDDO
      ysca= dble(ycut)*dble(evis)**2
C---FILL V-ARRAY: V(I,J) IS V(NT*(I-1)+J-I(I+1)/2)
      iad = 0
      DO i=1,nt-1
        DO ii=1,5
          pp(ii)=pl(ii,i)
        ENDDO
        DO j=i+1,nt
          iad = iad + 1
          v(iad) = 2d0*(1d0-(pp(1)*pl(1,j) +pp(2)*pl(2,j)
     &                      +pp(3)*pl(3,j))*pp(5)*pl(5,j))
        ENDDO
      ENDDO
      nj=nt
C---START MAIN LOOP.  FIRST LOOK FOR MINIMUM V
    1 vmini = 1d10
      imini = 0
      iad = 0
      DO i=1,nt-1
        IF (ip(i)) THEN
          DO j=i+1,nt
            iad = iad + 1
            IF (ip(j).AND.v(iad).LT.vmini) THEN
              vmini = v(iad)
              imini = i
              jmini = j
            ENDIF
          ENDDO
        ELSE
          iad=iad+nt-i
        ENDIF
      ENDDO
C---END OF CLUSTER SEARCH FOR VMINI
      IF (imini.NE.0) THEN
C---NOT FINISHED YET
        IF (vmini*min(pl(4,imini),pl(4,jmini))**2.GE.ysca) THEN
C---SOFT FREEZING HERE
          IF (pl(4,imini).LT.pl(4,jmini)) THEN
            ip(imini)=.false.
          ELSE
            ip(jmini)=.false.
          ENDIF
        ELSE
C---COMBINE PARTICLES IMINI AND JMINI
          DO ii=1,4
            pl(ii,imini)=pl(ii,imini)+pl(ii,jmini)
          ENDDO
          pm=pl(1,imini)**2+pl(2,imini)**2+pl(3,imini)**2
          IF (pm.GT.0d0) THEN
            pl(5,imini)=1d0/sqrt(pm)
          ELSE
            pl(5,imini)=1d0
          ENDIF
C---FLAG PARTICLE JMINI AS COMBINED
          ip(jmini)=.false.
          ij(jmini)=imini
          nj=nj-1
C---CALCULATE RELEVANT NEW V VALUES
          DO i=1,nt
            IF (i.NE.imini) THEN
              IF (ij(i).EQ.jmini) ij(i)=imini
              IF (ip(i)) THEN
                k = min(i,imini)
                l = max(i,imini)
                iad = nt*(k-1) + l - (k*(k+1))/2
                v(iad) = 2d0*(1d0-(pl(1,k)*pl(1,l) +pl(2,k)*pl(2,l)
     &                            +pl(3,k)*pl(3,l))*pl(5,k)*pl(5,l))
              ENDIF
            ENDIF
          ENDDO
        ENDIF
C---BACK TO START OF LOOP
        go to 1
      ELSE
C---FINISHED: CONSTRUCT JETS
        j=0
        DO i=1,nt
          IF (ij(i).EQ.i) THEN
            j=j+1
            jj(i)=j
            DO ii=1,4
              pj(ii,j)= sngl(pl(ii,i))
            ENDDO
          ENDIF
        ENDDO
        DO i=1,nt
          ij(i)=jj(ij(i))
        ENDDO
      ENDIF
      END