ykern.f

Go to the documentation of this file.

CDECK  ID>, YKERN.
      SUBROUTINE ykern(IMODE,NT,ITKDM,PP,IERR)
C
C  JET FINDING ROUTINE A LA THE ORIGINAL JADE (E0) SCHEME, INCLUDING
C  THE E-,P- AND P0-VARIANTS. ALSO FEATURES THE NEW DURHAM (D)
C  AND GENEVA (G) SCHEMES.
C
C  INPUTS: IMODE (JET FINDING SCHEME); PP-ARRAY CONTAINING THE
C  FOUR-MOMENTA OF THE SELECTED PARTICLES, NT SPECIFYING HOW
C  MANY LOCATIONS IN PP ARE FILLED.
C
C  OUTPUTS:  YREC(I) CONTAINS VALUE OF Y WHERE EVENT FLIPS FROM
C  (I+1)-JET TO I-JET; PJET(K,I,J) CONTAINS THE JET AXES FOUR VECTOR
C  (K=1-4) OF JET NUMBER I WHEN THE EVENT IS CLASSIFIED AS J-JET.
C  LOCATION K=7 GIVES THE NUMBER OF PARTICLES ASSIGNED TO THIS JET.
C  NOTE THAT EACH EVENT IS ALWAYS FULLY RECONSTRUCTED DOWN TO
C  1-JET CONFIGURATION.
C  PROGRAM IS SELF-CONTAINED.
C
C  IMODE=1: E0 (=JADE) 2: P  3: P0  4: E  5: DURHAM (KT)  6: GENEVA
C
C  LAST MOD  :  21-Jun-98
C
C  Modification Log.
C  08-Oct-97 D. Chrisman, Increase NYCLMX from 250 to 500.
C  21-Jun-98 D. Chrisman, Want to be able to handle more than 10 jets.
C                         Introduce parameter NJETMX. Use this
C                         in declaration of PJET and YREC in order to remove
C                         hard coded numbers. Remove variable JCHECK and
C                         replace with NJETMX.
C
      IMPLICIT NONE
C     IMPLICIT NONE
      INTEGER nt,itkdm,nyclmx,njetmx,ncall,nprint,i,j,k,ierr,imodeo
      INTEGER imode
      parameter(nyclmx = 500, njetmx = 20)
      INTEGER nto,njeto,njj,im,jm,nold,kk,ii,histor(2,nyclmx)
      REAL pl(7,nyclmx),y(nyclmx,nyclmx),pint(10,nyclmx)
      REAL yrec(njetmx),eviso
      REAL pp(itkdm,*),jade,d,g,e,evis,pvis,pjet(10,njetmx,njetmx),ymini
      COMMON /ycl/yrec,pjet,histor
      COMMON /yint/ imodeo,nto,njeto,pint
      CHARACTER*7 cm
      SAVE ncall,nprint,y
      DATA ncall,nprint /0,0/
C
C  JET RESOLUTION FUNCTIONS
C
      jade(i,j) = 2.*pl(4,i)*pl(4,j)*max(0.,(1.-
     + (pl(1,i)*pl(1,j)+pl(2,i)*pl(2,j)+pl(3,i)*pl(3,j))/
     + (pl(6,i)*pl(6,j))))
      d(i,j) = 2.*min(pl(4,i)*pl(4,i),pl(4,j)*pl(4,j))*max(0.,(1.-
     + (pl(1,i)*pl(1,j)+pl(2,i)*pl(2,j)+pl(3,i)*pl(3,j))/
     + (pl(6,i)*pl(6,j))))
      e(i,j) = max(0.,(pl(4,i)+pl(4,j))**2-(pl(1,i)+pl(1,j))**2-
     + (pl(2,i)+pl(2,j))**2-(pl(3,i)+pl(3,j))**2)
      g(i,j) = 8.*pl(4,i)*pl(4,j)*max(0.,(1.-
     + (pl(1,i)*pl(1,j)+pl(2,i)*pl(2,j)+pl(3,i)*pl(3,j))/
     + (pl(6,i)*pl(6,j))))/(9.*(pl(4,i)+pl(4,j))**2)
C
C INITIALIZE
C
      ierr = -1
      IF(ncall.LE.0) THEN
        imodeo = 0
        nto = 0
        njeto = 0
      ENDIF
      DO 5001 i=1,njetmx
        yrec(i) = 0.
 5001 CONTINUE
      ncall = ncall + 1
C
C PRINT JET SCHEME
C
      IF(imode.NE.imodeo .AND. nprint.LE.7) THEN
        IF(imode.EQ.1) THEN
          cm = 'JADE E0'
        ELSEIF(imode.EQ.2) THEN
          cm = 'JADE P '
        ELSEIF(imode.EQ.3) THEN
          cm = 'JADE P0'
        ELSEIF(imode.EQ.4) THEN
          cm = 'JADE E '
        ELSEIF(imode.EQ.5) THEN
          cm = 'DURHAM '
        ELSEIF(imode.EQ.6) THEN
          cm = 'GENEVA '
        ELSE
          WRITE(6,281) imode
 281      FORMAT(/,' ### YKERN: IMODE =',i3,' INVALID; SET TO 1 ###')
          imode = 1
          cm = 'JADE E0'
        ENDIF
        print 789, cm
 789    FORMAT(/,8x,54('#'),/,8x,
     +  '# YCLUS JET FINDER WITH ',a7,' RECOMBINATION SCHEME #',/,
     +  8x,54('#'),/)
        nprint = nprint + 1
      ENDIF
C
C CHECK INPUT PARAMETERS
C
      IF(nt.LE.1) THEN
        print 701, nt
 701    FORMAT(/,' ### YKERN: NUMBER OF INPUT PARTICLES ',i4,' < 2;',
     +  ' NO JET RECONSTRUCTION DONE. ###')
        RETURN
      ENDIF
      IF(nt.GT.nyclmx) THEN
        print 700, nt, nyclmx,nyclmx
  700   FORMAT(/,' #### YKERN: NUMBER OF INPUT PARTICLES ',i4,' > ',i4,
     +  '; SET TO ',i4,/,12x,' INCREASE NYCLMX IF THIS',
     +  ' WARNING OCCURS MORE OFTEN. #### ')
        nt = nyclmx
      ENDIF
C
C COPY INPUT VECTORS INTO INTERNAL MOMENTUM ARRAY (PL)
C  (POSITION  7 OF EACH VECTOR IN PL WILL BE OVERWRITTEN BY THIS ROUTINE
C   IN ORDER TO RECORD NUMBER OF INITIAL PARTICLES BELONGING TO THIS
C   POSITION [= 1. INITIALLY])
C
      evis = 0.
      pvis = 0.
      DO 5002 i=1,nt
        pl(6,i)=sqrt(pp(1,i)**2 + pp(2,i)**2 + pp(3,i)**2)
        pl(1,i)=pp(1,i)
        pl(2,i)=pp(2,i)
        pl(3,i)=pp(3,i)
        IF(imode.EQ.2 .OR. imode.EQ.3) THEN
          pl(4,i) = pl(6,i)
        ELSE
          pl(4,i) = pp(4,i)
        ENDIF
        pl(7,i) = 1.
        DO 5003 k=1,7
          pint(k,i) = pl(k,i)
 5003   CONTINUE
        evis=evis+pl(4,i)
        pvis=pvis+pl(6,i)
 5002 CONTINUE
      njj = nt
C
      IF(evis.LE.0. .OR. pvis.GT.1.001*evis) THEN
        WRITE(6,283) evis,pvis
 283    FORMAT(' #### YKERN: INCOMPATIBLE SUMS OF ENERGIES',
     +  ' AND/OR MOMENTA:',
     +  /,12x,' SUM(E) = ',f11.4,' SUM(P) = ',f11.4,/,12x,' CHECK',
     +  ' INPUT VECTORS AND/OR ARRAY DIMENSIONS. ####')
        RETURN
      ENDIF
C
      IF(njj.LE.njetmx) THEN
        DO 5004 j=1,njj
          DO 5005 i=1,7
            pjet(i,j,njj) =  pl(i,j)
 5005     CONTINUE
 5004   CONTINUE
      ENDIF
C
C CALCULATE AND STORE PAIR MASSES
C
       DO 5006 i=1,njj-1
         DO 5007 j=i+1,njj
           IF(imode.EQ.1) THEN
             y(i,j) = jade(i,j)
           ELSEIF(imode.EQ.5) THEN
             y(i,j) = d(i,j)
           ELSEIF(imode.EQ.6) THEN
             y(i,j) = g(i,j)
           ELSE
             y(i,j) = e(i,j)
           ENDIF
 5007    CONTINUE
 5006  CONTINUE
C
C FIND LOCAL MINIMUM OF PAIR MASSES
C
      im = 0
      jm = 0
 1000 CONTINUE
      ymini = 2.*evis**2
      DO 5008 i=1,njj-1
        DO 5009 j=i+1,njj
          IF(y(i,j).LT.ymini) THEN
           ymini = y(i,j)
           im = i
           jm = j
          ENDIF
 5009   CONTINUE
 5008 CONTINUE
C
      IF(im.EQ.0 .OR. jm.EQ.0) THEN
        WRITE(6,284)
  284   FORMAT(' #### YKERN: ERROR; NO MINIMUM FOUND IN Y-ARRAY! ####')
        RETURN
      ENDIF
C
C RECOMBINE PARTICLES IM AND JM, STORE AT POSITION IM
C
      pl(1,im) = pl(1,im) + pl(1,jm)
      pl(2,im) = pl(2,im) + pl(2,jm)
      pl(3,im) = pl(3,im) + pl(3,jm)
      pl(6,im) = sqrt(pl(1,im)**2 + pl(2,im)**2
     +         + pl(3,im)**2)
      IF(imode.EQ.2) THEN
        pl(4,im) = pl(6,im)
      ELSEIF(imode.EQ.3) THEN
        eviso = evis
        evis = evis-pl(4,im)-pl(4,jm)+pl(6,im)
        pl(4,im) = pl(6,im)
      ELSE
        pl(4,im) = pl(4,im) + pl(4,jm)
      ENDIF
C                  KEEP TRACK OF # OF PARTICLES ASSIGNED TO NEW CLUSTER:
      pl(7,im) = pl(7,im) + pl(7,jm)
C                  MOVE LAST PARTICLE IN LIST (NJJ) TO EMPTY SLOT (JM)
      nold = 0
      IF(jm.NE.njj) THEN
        DO 5010 kk=1,7
          pl(kk,jm) = pl(kk,njj)
 5010   CONTINUE
        nold = njj
      ENDIF
      histor(1,njj) = jm
      histor(2,njj) = im
      njj = njj - 1
C                  REMEMBER JET AXES AND VALUE OF YIJ OF FLIP
      IF(njj.LE.njetmx) THEN
        IF(imode.EQ.3) THEN
          yrec(njj) = ymini/eviso**2
        ELSEIF(imode.EQ.6) THEN
          yrec(njj) = ymini
        ELSE
          yrec(njj) = ymini/evis**2
        ENDIF
        DO 5011 i=1,njj
          DO 5012 k=1,7
            pjet(k,i,njj) = pl(k,i)
 5012     CONTINUE
 5011   CONTINUE
      ENDIF
C
C END IF 1-JET CASE REACHED
C
      IF(njj.LE.1) goto 9000
C
C NOW CALCULATE RELEVANT NEW MASS-COMBINATIONS
C
      DO 5013 ii=1,njj
        IF(ii.NE.im) THEN
         i = min(ii,im)
         j = max(ii,im)
         IF(imode.EQ.1) THEN
           y(i,j) = jade(i,j)
         ELSEIF(imode.EQ.5) THEN
           y(i,j) = d(i,j)
         ELSEIF(imode.EQ.6) THEN
           y(i,j) = g(i,j)
         ELSE
           y(i,j) = e(i,j)
         ENDIF
        ENDIF
        IF(nold.NE.0) THEN
          i = min(ii,jm)
          j = max(ii,jm)
          y(i,j) = y(ii,nold)
        ENDIF
 5013 CONTINUE
C
C  BACK TO START OF LOOP
C
      goto 1000
C
 9000 CONTINUE
C
      imodeo = imode
      nto = nt
      ierr = 0
CCCCC
C     WRITE(6,825) (J,(HISTOR(I,J),I=1,2),J=1,NT)
C825  FORMAT(/,' HISTORY:',/,250(I3,4X,2I4,/))
CCCCC
      RETURN
      END