pxluc5.f

Go to the documentation of this file.

CDECK  ID>, PXLUC5.
      SUBROUTINE pxluc5 (N,NRLUDM,P,MSTU46,MSTU47,
     +                   paru44,paru45,njet,ijmul,ipass)
*.*********************************************************
*. ------
*. PXLUC5
*. ------
*. An "in-house" version of the Jetset jet-finding algorithm
*. which works entirely through an argument list rather than
*. through e.g. the Jetset common blocks. Its operation is
*. therefore entirely decoupled from the the operation of Jetset
*. (i.e. the values of MST or MSTJ etc. in Jetset common do not
*. affect this routine whatsoever).
*. The main purpose of an in-house version of the
*. Jetset jetfinding algorithm is to have a version
*. which is compatible with both Jetset6.3 and Jetset7.1 etc.
*. (because of the change in the Jetset common blocks between
*. these two versions, the version of this algorithm in the
*. Jetset library is version specific).
*. The input arguments MSTU46, MSTU47, PARU44, PARU45 correspond
*. to the parameters MSTU(46), MSTU(47), PARU(44), PARU(45) of
*. Jetset7.1, see "A manual to ... Jetset7.1," T.Sjostrand
*. (filename "Jetset7.1 MANUAL A" on the Opal generator disk).
*.
*.      INTEGER  NRLUDM,MXJET
*.      PARAMETER (NRLUDM=1000.or.so,MXJET=10.or.so)
*.      INTEGER  IPASS (NRLUDM),IJMUL (MXJET)
*.      INTEGER  NTRAK,NJET
*.      REAL PLUND (NRLUDM,5)
*.      REAL  MSTU46,MSTU47,PARU44,PARU45
*.
*.      (define NTRAK, fill PLUND)
*.      CALL PXLUC5 (NTRAK,NRLUDM,PLUND,MSTU46,MSTU47,
*.     +             PARU44,PARU45,NJET,IJMUL,IPASS)
*.
*. INPUT     : NTRAK    Number of tracks
*. INPUT     : NRLUDM   First dimension of PLUND
*. IN/OUTPUT : PLUND    4-momenta in Jetset format
*. INPUT     : MSTU46   same as MSTU(46) in Jetset7.1:jet-finder mode
*. INPUT     : MSTU47   same as MSTU(47) in Jetset7.1
*. INPUT     : PARU44   same as PARU(44) in Jetset7.1
*. INPUT     : PARU45   same as PARU(45) in Jetset7.1
*. OUTPUT    : NJET      Number of jets found
*. OUTPUT    : IJMUL(i)  Jet i contains IJMUL(i) particles
*. OUTPUT    : IPASS(k)  Particle k belongs to jet number IPASS(k)
*.
*. CALLS     : PXANXY,PXPLU3,PXRMX3,PXROF3,PXROB3
*. CALLED    : PXLTH4
*.
*. AUTHOR    : Modified from LUCLUS (T.Sjostrand) by J.W.Gary
*. CREATED   : 31-Jan-89
*. LAST MOD  : 31-Jan-89
*.
*. Modification Log.
*. 05-May-97 D. Chrisman, remove declaration of unused variables
*.             I1, I2 and PMAS.
*.
*.*********************************************************
      IMPLICIT NONE
      INTEGER  locdim
      parameter(locdim=2000)
      REAL paru48,pimas,paru43
      parameter(paru48=0.0001,pimas=0.13957,paru43=0.25)
      INTEGER  mstu42,mstu48,mstu43,mstu46,mstu47,itry1,itry2,
     +         i,j,idel,irec,iemp,ijet,imin,imax,njet,nrem,
     +         inew,iori,npre,ispl,nsav,itry,np,imin1,imin2,nrludm,
     +         n,nloop,ip,ij
      INTEGER  k (locdim,5),ijmul (*),ipass (*)
      REAL  paru44,paru45,r2max,pemax,rinit,pxrr2m,pxrr2t,pss,
     +      pmax,r2,tsav,psjt,r2acc,r2min,pxmas
      REAL  p(nrludm,*),v(locdim,5),ps(5)
      DATA  mstu42 / 2 /,mstu48 / 0 /, mstu43 / 1 /

CC...If first time, reset. If reentering, skip preliminaries.
**JWG      IF(MSTU48.LE.0) THEN
        np=0
        DO 100 j=1,5
  100   ps(j)=0.
        pss=0.
**JWG      ELSE
**JWG        NJET=NSAV
**JWG        IF(MSTU43.GE.2) N=N-NJET
**JWG        DO 110 I=N+1,N+NJET
**JWG  110   P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
**JWG        IF(MSTU46.LE.3) R2ACC=PARU44**2
**JWG        IF(MSTU46.GE.4) R2ACC=PARU45*PS(5)**2
**JWG        NLOOP=0
**JWG        GOTO 290
**JWG      ENDIF
C...Find which particles are to be considered in cluster search.
      DO 140 i=1,n
      IF(n+2*np.GE.locdim-5) THEN
          WRITE (6,fmt='('' PXLUC5: Error, not enough buffer'',
     +           ''space for jet-finer calculation'')')
          njet = -1
          go to 990
      ENDIF
C...Take copy of these particles, with space left for jets later on.
      np=np+1
      k(n+np,3)=i
      DO 120 j=1,5
  120 p(n+np,j)=p(i,j)
**JWG      IF(MSTU42.EQ.0) P(N+NP,5)=0.
      pxmas = p(i,5)
**JWG      IF(MSTU42.EQ.1.AND.PXMAS.NE.0) P(N+NP,5)=PIMAS
      p(n+np,4)=sqrt(p(n+np,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
      p(n+np,5)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
      DO 130 j=1,4
  130 ps(j)=ps(j)+p(n+np,j)
      pss=pss+p(n+np,5)
  140 CONTINUE
      DO 150 i=n+1,n+np
      k(i+np,3)=k(i,3)
      DO 150 j=1,5
  150 p(i+np,j)=p(i,j)
      ps(5)=sqrt(max(0.,ps(4)**2-ps(1)**2-ps(2)**2-ps(3)**2))
C...Very low multiplicities not considered.
      IF(np.LT.mstu47) THEN
        njet=-1
        RETURN
      ENDIF
C...Find precluster configuration. If too few jets, make harder cuts.
      nloop=0
      IF(mstu46.LE.3) r2acc=paru44**2
      IF(mstu46.GE.4) r2acc=paru45*ps(5)**2
      rinit=1.25*paru43
      IF(np.LE.mstu47+2) rinit=0.
  160 rinit=0.8*rinit
      npre=0
      nrem=np
      DO 170 i=n+np+1,n+2*np
  170 k(i,4)=0
C...Sum up small momentum region. Jet if enough absolute momentum.
      IF(mstu46.LE.2) THEN
        DO 180 j=1,4
  180   p(n+1,j)=0.
        DO 200 i=n+np+1,n+2*np
        IF(p(i,5).GT.2.*rinit) goto 200
        nrem=nrem-1
        k(i,4)=1
        DO 190 j=1,4
  190   p(n+1,j)=p(n+1,j)+p(i,j)
  200   CONTINUE
        p(n+1,5)=sqrt(p(n+1,1)**2+p(n+1,2)**2+p(n+1,3)**2)
        IF(p(n+1,5).GT.2.*rinit) npre=1
        IF(rinit.GE.0.2*paru43.AND.npre+nrem.LT.mstu47) goto 160
      ENDIF
C...Find fastest remaining particle.
  210 npre=npre+1
      pmax=0.
      DO 220 i=n+np+1,n+2*np
      IF(k(i,4).NE.0.OR.p(i,5).LE.pmax) goto 220
      imax=i
      pmax=p(i,5)
  220 CONTINUE
      DO 230 j=1,5
  230 p(n+npre,j)=p(imax,j)
      nrem=nrem-1
      k(imax,4)=npre
C...Sum up precluster around it according to pT separation.
      IF(mstu46.LE.2) THEN
        DO 250 i=n+np+1,n+2*np
        IF(k(i,4).NE.0) goto 250
        r2=pxrr2t(nrludm,p,i,imax)
        IF(r2.GT.rinit**2) goto 250
        nrem=nrem-1
        k(i,4)=npre
        DO 240 j=1,4
  240   p(n+npre,j)=p(n+npre,j)+p(i,j)
  250   CONTINUE
        p(n+npre,5)=sqrt(p(n+npre,1)**2+p(n+npre,2)**2+p(n+npre,3)**2)
C...Sum up precluster around it according to mass separation.
      ELSE
  260   imin=0
        r2min=rinit**2
        DO 270 i=n+np+1,n+2*np
        IF(k(i,4).NE.0) goto 270
        r2=pxrr2m(nrludm,p,i,n+npre)
        IF(r2.GE.r2min) goto 270
        imin=i
        r2min=r2
  270   CONTINUE
        IF(imin.NE.0) THEN
          DO 280 j=1,4
  280     p(n+npre,j)=p(n+npre,j)+p(imin,j)
          p(n+npre,5)=sqrt(p(n+npre,1)**2+p(n+npre,2)**2+p(n+npre,3)**2)
          nrem=nrem-1
          k(imin,4)=npre
          goto 260
        ENDIF
      ENDIF
C...Check if more preclusters to be found. Start over if too few.
      IF(rinit.GE.0.2*paru43.AND.npre+nrem.LT.mstu47) goto 160
      IF(nrem.GT.0) goto 210
      njet=npre
C...Reassign all particles to nearest jet. Sum up new jet momenta.
  290 tsav=0.
      psjt=0.
  300 IF(mstu46.LE.1) THEN
        DO 310 i=n+1,n+njet
        DO 310 j=1,4
  310   v(i,j)=0.
        DO 340 i=n+np+1,n+2*np
        r2min=pss**2
        DO 320 ijet=n+1,n+njet
        IF(p(ijet,5).LT.rinit) goto 320
        r2=pxrr2t(nrludm,p,i,ijet)
        IF(r2.GE.r2min) goto 320
        imin=ijet
        r2min=r2
  320   CONTINUE
        k(i,4)=imin-n
        DO 330 j=1,4
  330   v(imin,j)=v(imin,j)+p(i,j)
  340   CONTINUE
        psjt=0.
        DO 360 i=n+1,n+njet
        DO 350 j=1,4
  350   p(i,j)=v(i,j)
        p(i,5)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
  360   psjt=psjt+p(i,5)
      ENDIF
C...Find two closest jets.
      r2min=2.*r2acc
      DO 370 itry1=n+1,n+njet-1
      DO 370 itry2=itry1+1,n+njet
      IF(mstu46.LE.2) r2=pxrr2t(nrludm,p,itry1,itry2)
      IF(mstu46.GE.3) r2=pxrr2m(nrludm,p,itry1,itry2)
      IF(r2.GE.r2min) goto 370
      imin1=itry1
      imin2=itry2
      r2min=r2
  370 CONTINUE
C...If allowed, join two closest jets and start over.
      IF(njet.GT.mstu47.AND.r2min.LT.r2acc) THEN
        irec=min(imin1,imin2)
        idel=max(imin1,imin2)
        DO 380 j=1,4
  380   p(irec,j)=p(imin1,j)+p(imin2,j)
        p(irec,5)=sqrt(p(irec,1)**2+p(irec,2)**2+p(irec,3)**2)
        DO 390 i=idel+1,n+njet
        DO 390 j=1,5
  390   p(i-1,j)=p(i,j)
        IF(mstu46.GE.2) THEN
          DO 400 i=n+np+1,n+2*np
          iori=n+k(i,4)
          IF(iori.EQ.idel) k(i,4)=irec-n
  400     IF(iori.GT.idel) k(i,4)=k(i,4)-1
        ENDIF
        njet=njet-1
        goto 290
C...Divide up broad jet if empty cluster in list of final ones.
      ELSEIF(njet.EQ.mstu47.AND.mstu46.LE.1.AND.nloop.LE.2) THEN
        DO 410 i=n+1,n+njet
  410   k(i,5)=0
        DO 420 i=n+np+1,n+2*np
  420   k(n+k(i,4),5)=k(n+k(i,4),5)+1
        iemp=0
        DO 430 i=n+1,n+njet
  430   IF(k(i,5).EQ.0) iemp=i
        IF(iemp.NE.0) THEN
          nloop=nloop+1
          ispl=0
          r2max=0.
          DO 440 i=n+np+1,n+2*np
          IF(k(n+k(i,4),5).LE.1.OR.p(i,5).LT.rinit) goto 440
          ijet=n+k(i,4)
          r2=pxrr2t(nrludm,p,i,ijet)
          IF(r2.LE.r2max) goto 440
          ispl=i
          r2max=r2
  440     CONTINUE
          IF(ispl.NE.0) THEN
            ijet=n+k(ispl,4)
            DO 450 j=1,4
            p(iemp,j)=p(ispl,j)
  450       p(ijet,j)=p(ijet,j)-p(ispl,j)
            p(iemp,5)=p(ispl,5)
            p(ijet,5)=sqrt(p(ijet,1)**2+p(ijet,2)**2+p(ijet,3)**2)
            IF(nloop.LE.2) goto 290
          ENDIF
        ENDIF
      ENDIF
C...If generalized thrust has not yet converged, continue iteration.
      IF(mstu46.LE.1.AND.nloop.LE.2.AND.psjt/pss.GT.tsav+paru48)
     +THEN
        tsav=psjt/pss
        goto 300
      ENDIF
C...Reorder jets according to energy.
      DO 460 i=n+1,n+njet
      DO 460 j=1,5
  460 v(i,j)=p(i,j)
      DO 490 inew=n+1,n+njet
      pemax=0.
      DO 470 itry=n+1,n+njet
      IF(v(itry,4).LE.pemax) goto 470
      imax=itry
      pemax=v(itry,4)
  470 CONTINUE
      k(inew,1)=31
      k(inew,2)=97
      k(inew,3)=inew-n
      k(inew,4)=0
      DO 480 j=1,5
  480 p(inew,j)=v(imax,j)
      v(imax,4)=-1.
  490 k(imax,5)=inew
C...Clean up particle-jet assignments and jet information.
      DO 500 i=n+np+1,n+2*np
      iori=k(n+k(i,4),5)
      k(i,4)=iori-n
      IF(k(k(i,3),1).NE.3) k(k(i,3),4)=iori-n
      k(iori,4)=k(iori,4)+1
  500 CONTINUE
      iemp=0
      psjt=0.
      DO 520 i=n+1,n+njet
      k(i,5)=0
      psjt=psjt+p(i,5)
      p(i,5)=sqrt(max(p(i,4)**2-p(i,5)**2,0.))
      DO 510 j=1,5
  510 v(i,j)=0.
  520 IF(k(i,4).EQ.0) iemp=i
C...Select storing option. Output variables. Check for failure.
      IF(iemp.NE.0) THEN
        njet=-1
      ENDIF
      nsav=njet
      DO 560  ij = 1,njet
          ijmul(ij) = k(n+ij,4)
 560  CONTINUE
      DO 580  ip = 1,np
          ipass(ip) = k(ip,4)
 580  CONTINUE

 990  RETURN
      END

C...Functions: distance measure in pT or (pseudo)mass.
      FUNCTION pxrr2t (NRLUDM,P,I1,I2)
      IMPLICIT NONE
      INTEGER  nrludm,i1,i2
      REAL  p (nrludm,*),pxrr2t
      pxrr2t = (p(i1,5)*p(i2,5)-p(i1,1)*p(i2,1)-p(i1,2)*p(i2,2)-
     +p(i1,3)*p(i2,3))*2.*p(i1,5)*p(i2,5)/(0.0001+p(i1,5)+p(i2,5))**2
      RETURN
      END

      FUNCTION pxrr2m (NRLUDM,P,I1,I2)
      IMPLICIT NONE
      INTEGER  nrludm,i1,i2
      REAL  p (nrludm,*),pxrr2m
      pxrr2m = 2.*p(i1,4)*p(i2,4)*(1.-(p(i1,1)*p(i2,1)+p(i1,2)*
     +p(i2,2)+p(i1,3)*p(i2,3))/(p(i1,5)*p(i2,5)))
      RETURN
      END