ckern.f

Go to the documentation of this file.

CDECK  ID>, CKERN.
      SUBROUTINE ckern(ITKDM,NT,PT,NJREQ,IERR)
*.-----------------------------------------------------------------------
*.
*.    CKERN: Resolves the clustering, and stores information
*.           in common block CKCOM.
*.           INPUT: ITKDM      (integer) First dimension of PT array
*.                                 (>=4 required)
*.                  NT         (integer) Number of tracks
*.                  PT(ITKDM,*)(real)    Four-momenta of tracks
*.                  NJREQ      (integer) Clustering limit. Note that the CPU
*.                                       consumption is roughly proportional to
*.                                       calling JADE/Durham algorithm NJREQ
*.                                       times!!
*.                                       1) If interested in ALL possible N-jet
*.                                          configurations, set NJREQ equal to
*.                                          NT. Note the remark on CPU time!
*.                                       2) If interested in first N-jet
*.                                          configuration (with largest
*.                                          ycut values), set NJREQ equal
*.                                          to number of jets N-jet.
*.           OUTPUT: IERR      (integer) Error flag, 0=OK.
*.           CALLS:  CKSORD
*.
*.
*.    Once CKERN is called, all information on the event clustering,
*.    final state jets and particle association can be accessed
*.    without additional CPU time using the utility routines.
*.
*.---CAMBRIDGE JET CLUSTERING ALGORITHM
*.   BASED ON YCLUS BY S BETHKE
*.   REF: YU L DOKSHITZER, G D LEDER, S MORETTI, B R WEBBER
*.   CAVENDISH-HEP-97/06 (JUNE 1997)
*.   07/07/97 FIRST RELEASED BY BRW
*.   23/08/97 COMMENTS REVISED BY BRW
*.   23/09/97 IMPLEMENT FOR PX LIBRARY BY STK AND SB
*.
*.
*.  CREATED :  11-12-1997, STAN BENTVELSEN
*.  LAST MOD:  11-02-1998, Z. Troscyani, Add factor epsilon
*.-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER nmxy , nmxp
      parameter(nmxy = 300)
      parameter(nmxp = 30)
      DOUBLE PRECISION ytrans(nmxy)
      DOUBLE PRECISION pcmj(nmxp,4,nmxp)
      INTEGER          ntrans(nmxy), njiter, njmax, ntrack
      INTEGER          icmj(nmxp,nmxy)
      COMMON / ckcom / ytrans, pcmj, ntrans, icmj, njiter, njmax, ntrack
      INTEGER itkdm,nt,nj,njreq,ierr
      REAL pt(itkdm,*)
      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,ij(ntrk),jj(ntrk)
      DOUBLE PRECISION pp(5),pl(5,ntrk),v(nv),pm,vmini,ysca,evis
C
      DOUBLE PRECISION epsilon
      DOUBLE PRECISION ycur(ntrk)
      INTEGER          icur, kix(ntrk)
C

      SAVE lcall
      DATA lcall / .false. /
      DATA epsilon / 1d-12 /
      SAVE epsilon

C  Welcome message:
      IF( .NOT.lcall ) THEN
c        PRINT *, ' '
c        PRINT *, 'Cambridge jet finding algorithm, please refer to:'
c        PRINT *, 'Yu.L. Dokshitzer, G.D. Leder, S. Moretti, B.R. Webber'
c        PRINT *, 'CAVENDISH-HEP-97/06 (June 1997)'
        print *, ' '
        print *, ' #################################################'
        print *, ' ###     C K E R N    V E R S I O N    104     ###'
c        PRINT *, ' ###  S. Bentvelsen and I. Meyer, EPXXX, 1998  ###'
        print *, ' #################################################'
        WRITE(*,'(A,I2,A)')
     +      '  ##### CAMBRIDGE JET CLUSTERING 1 - ',njreq,' JETS  #####'
        print *, ' #################################################'
        lcall= .true.
      ENDIF
C---WARNINGS
      ierr = 0
      evis = 0d0
      icur = 0
      ntrack = nt
      njreq = min(njreq,nt)
      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
C..   RESET ALL VALUES
C
      njiter = 0
      DO i=1,nmxy
         ytrans(i) = 0d0
         ntrans(i) = 0
         DO j=1,nmxp
            icmj(j,i)   = 0
         ENDDO
      ENDDO

      DO i=1,nt
         evis = evis + dble(pt(4,i))
      ENDDO
      ysca   = evis*evis
C---COPY MOMENTA INTO PL-ARRAY
 2    DO i=1,nt
         ycur(i) = 0d0
         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
         IF (pm.GT.0d0) THEN
            pl(5,i)=1d0/sqrt(pm)
         ELSE
            pl(5,i)=1d0
         ENDIF
      ENDDO
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
         icur       = icur + 1
         ycur(icur) = vmini*min(pl(4,imini),pl(4,jmini))**2
         IF (ycur(icur)+epsilon.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
C..   ITERATION IN JET-FINDING FINISHED
C
         njiter = njiter + 1
         IF(njiter.LE.nmxy) THEN
C
C..   STORE THE VALUE OF YFLIP AND THE CORRESPONDING NUMBER OF JETS
C
            ytrans(njiter) = ysca
            ntrans(njiter) = nj

            IF(njiter.LE.nmxp) THEN
C
C..   STORE THE JET-DIRECTIONS AND PARTICLE ASSOCIATION
C
               j=0
               DO i=1,nt
                  IF (ij(i).EQ.i) THEN
                     j=j+1
                     jj(i)=j
                     DO ii=1,4
                        pcmj(njiter,ii,j)= pl(ii,i)
                     ENDDO
                  ENDIF
               ENDDO
               DO i=1,nt
                  icmj(njiter,i)=jj(ij(i))
               ENDDO
            ENDIF
            IF(nj.LE.njreq) THEN
C
C..   ANOTHER ITERATION?
C
               CALL cksord(nt,ycur,kix,'I')
               DO i=1,nt
C
C..   DETERMINE THE NEXT VALUE FOR YSCA, RESET VARIABLES
C
                  IF(ycur(kix(i)).LT.ysca) THEN
                     ysca = ycur(kix(i))
                     icur = 0
                     DO j=1,nt
                        ycur(j) = 0d0
                     ENDDO
C
C..   YUP, GO FOR THE NEXT ROUND
C
                     IF(ysca.NE.0) goto 2
                  ENDIF
               ENDDO
            ENDIF
         ENDIF
      ENDIF
C
C..   DEVIDE BY ECM**2
C
      njmax = -1
      DO i=1,njiter
         ytrans(i) = ytrans(i)/evis**2
         IF(ntrans(i).GT.njmax) njmax = ntrans(i)
      ENDDO

      END