C********************************************************************* C...PYCLUS C...Subdivides the particle content of an event into jets/clusters. SUBROUTINE PYCLUS(NJET) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) INTEGER PYK,PYCHGE,PYCOMP C...Parameter statement to help give large particle numbers. PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000, &KEXCIT=4000000,KDIMEN=5000000) C...Commonblocks. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) SAVE /PYJETS/,/PYDAT1/,/PYDAT2/ C...Local arrays and saved variables. DIMENSION PS(5) SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM SAVE C...Functions: distance measure in pT, (pseudo)mass or Durham pT. R2T(I1,I2)=(P(I1,5)*P(I2,5)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)- &P(I1,3)*P(I2,3))*2D0*P(I1,5)*P(I2,5)/(0.0001D0+P(I1,5)+P(I2,5))**2 R2M(I1,I2)=2D0*P(I1,4)*P(I2,4)*(1D0-(P(I1,1)*P(I2,1)+P(I1,2)* &P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5))) R2D(I1,I2)=2D0*MIN(P(I1,4),P(I2,4))**2*(1D0-(P(I1,1)*P(I2,1)+ &P(I1,2)*P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5))) C...If first time, reset. If reentering, skip preliminaries. IF(MSTU(48).LE.0) THEN NP=0 DO 100 J=1,5 PS(J)=0D0 100 CONTINUE PSS=0D0 PIMASS=PMAS(PYCOMP(211),1) ELSE NJET=NSAV IF(MSTU(43).GE.2) N=N-NJET DO 110 I=N+1,N+NJET P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) 110 CONTINUE IF(MSTU(46).LE.3.OR.MSTU(46).EQ.5) THEN R2ACC=PARU(44)**2 ELSE R2ACC=PARU(45)*PS(5)**2 ENDIF NLOOP=0 GOTO 300 ENDIF C...Find which particles are to be considered in cluster search. DO 140 I=1,N IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 140 IF(MSTU(41).GE.2) THEN KC=PYCOMP(K(I,2)) IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR. & KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR. & K(I,2).EQ.KSUSY1+39) GOTO 140 IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0) & GOTO 140 ENDIF IF(N+2*NP.GE.MSTU(4)-MSTU(32)-5) THEN CALL PYERRM(11,'(PYCLUS:) no more memory left in PYJETS') NJET=-1 RETURN 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 P(N+NP,J)=P(I,J) 120 CONTINUE IF(MSTU(42).EQ.0) P(N+NP,5)=0D0 IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PIMASS 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 PS(J)=PS(J)+P(N+NP,J) 130 CONTINUE PSS=PSS+P(N+NP,5) 140 CONTINUE DO 160 I=N+1,N+NP K(I+NP,3)=K(I,3) DO 150 J=1,5 P(I+NP,J)=P(I,J) 150 CONTINUE 160 CONTINUE PS(5)=SQRT(MAX(0D0,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2)) C...Very low multiplicities not considered. IF(NP.LT.MSTU(47)) THEN CALL PYERRM(8,'(PYCLUS:) too few particles for analysis') NJET=-1 RETURN ENDIF C...Find precluster configuration. If too few jets, make harder cuts. NLOOP=0 IF(MSTU(46).LE.3.OR.MSTU(46).EQ.5) THEN R2ACC=PARU(44)**2 ELSE R2ACC=PARU(45)*PS(5)**2 ENDIF RINIT=1.25D0*PARU(43) IF(NP.LE.MSTU(47)+2) RINIT=0D0 170 RINIT=0.8D0*RINIT NPRE=0 NREM=NP DO 180 I=N+NP+1,N+2*NP K(I,4)=0 180 CONTINUE C...Sum up small momentum region. Jet if enough absolute momentum. IF(MSTU(46).LE.2) THEN DO 190 J=1,4 P(N+1,J)=0D0 190 CONTINUE DO 210 I=N+NP+1,N+2*NP IF(P(I,5).GT.2D0*RINIT) GOTO 210 NREM=NREM-1 K(I,4)=1 DO 200 J=1,4 P(N+1,J)=P(N+1,J)+P(I,J) 200 CONTINUE 210 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.2D0*RINIT) NPRE=1 IF(RINIT.GE.0.2D0*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 170 IF(NREM.EQ.0) GOTO 170 ENDIF C...Find fastest remaining particle. 220 NPRE=NPRE+1 PMAX=0D0 DO 230 I=N+NP+1,N+2*NP IF(K(I,4).NE.0.OR.P(I,5).LE.PMAX) GOTO 230 IMAX=I PMAX=P(I,5) 230 CONTINUE DO 240 J=1,5 P(N+NPRE,J)=P(IMAX,J) 240 CONTINUE NREM=NREM-1 K(IMAX,4)=NPRE C...Sum up precluster around it according to pT separation. IF(MSTU(46).LE.2) THEN DO 260 I=N+NP+1,N+2*NP IF(K(I,4).NE.0) GOTO 260 R2=R2T(I,IMAX) IF(R2.GT.RINIT**2) GOTO 260 NREM=NREM-1 K(I,4)=NPRE DO 250 J=1,4 P(N+NPRE,J)=P(N+NPRE,J)+P(I,J) 250 CONTINUE 260 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 or C...Durham pT separation. ELSE 270 IMIN=0 R2MIN=RINIT**2 DO 280 I=N+NP+1,N+2*NP IF(K(I,4).NE.0) GOTO 280 IF(MSTU(46).LE.4) THEN R2=R2M(I,N+NPRE) ELSE R2=R2D(I,N+NPRE) ENDIF IF(R2.GE.R2MIN) GOTO 280 IMIN=I R2MIN=R2 280 CONTINUE IF(IMIN.NE.0) THEN DO 290 J=1,4 P(N+NPRE,J)=P(N+NPRE,J)+P(IMIN,J) 290 CONTINUE 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 270 ENDIF ENDIF C...Check if more preclusters to be found. Start over if too few. IF(RINIT.GE.0.2D0*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 170 IF(NREM.GT.0) GOTO 220 NJET=NPRE C...Reassign all particles to nearest jet. Sum up new jet momenta. 300 TSAV=0D0 PSJT=0D0 310 IF(MSTU(46).LE.1) THEN DO 330 I=N+1,N+NJET DO 320 J=1,4 V(I,J)=0D0 320 CONTINUE 330 CONTINUE DO 360 I=N+NP+1,N+2*NP R2MIN=PSS**2 DO 340 IJET=N+1,N+NJET IF(P(IJET,5).LT.RINIT) GOTO 340 R2=R2T(I,IJET) IF(R2.GE.R2MIN) GOTO 340 IMIN=IJET R2MIN=R2 340 CONTINUE K(I,4)=IMIN-N DO 350 J=1,4 V(IMIN,J)=V(IMIN,J)+P(I,J) 350 CONTINUE 360 CONTINUE PSJT=0D0 DO 380 I=N+1,N+NJET DO 370 J=1,4 P(I,J)=V(I,J) 370 CONTINUE P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2) PSJT=PSJT+P(I,5) 380 CONTINUE ENDIF C...Find two closest jets. R2MIN=2D0*MAX(R2ACC,PS(5)**2) DO 400 ITRY1=N+1,N+NJET-1 DO 390 ITRY2=ITRY1+1,N+NJET IF(MSTU(46).LE.2) THEN R2=R2T(ITRY1,ITRY2) ELSEIF(MSTU(46).LE.4) THEN R2=R2M(ITRY1,ITRY2) ELSE R2=R2D(ITRY1,ITRY2) ENDIF IF(R2.GE.R2MIN) GOTO 390 IMIN1=ITRY1 IMIN2=ITRY2 R2MIN=R2 390 CONTINUE 400 CONTINUE C...If allowed, join two closest jets and start over. IF(NJET.GT.MSTU(47).AND.R2MIN.LT.R2ACC) THEN IREC=MIN(IMIN1,IMIN2) IDEL=MAX(IMIN1,IMIN2) DO 410 J=1,4 P(IREC,J)=P(IMIN1,J)+P(IMIN2,J) 410 CONTINUE P(IREC,5)=SQRT(P(IREC,1)**2+P(IREC,2)**2+P(IREC,3)**2) DO 430 I=IDEL+1,N+NJET DO 420 J=1,5 P(I-1,J)=P(I,J) 420 CONTINUE 430 CONTINUE IF(MSTU(46).GE.2) THEN DO 440 I=N+NP+1,N+2*NP IORI=N+K(I,4) IF(IORI.EQ.IDEL) K(I,4)=IREC-N IF(IORI.GT.IDEL) K(I,4)=K(I,4)-1 440 CONTINUE ENDIF NJET=NJET-1 GOTO 300 C...Divide up broad jet if empty cluster in list of final ones. ELSEIF(NJET.EQ.MSTU(47).AND.MSTU(46).LE.1.AND.NLOOP.LE.2) THEN DO 450 I=N+1,N+NJET K(I,5)=0 450 CONTINUE DO 460 I=N+NP+1,N+2*NP K(N+K(I,4),5)=K(N+K(I,4),5)+1 460 CONTINUE IEMP=0 DO 470 I=N+1,N+NJET IF(K(I,5).EQ.0) IEMP=I 470 CONTINUE IF(IEMP.NE.0) THEN NLOOP=NLOOP+1 ISPL=0 R2MAX=0D0 DO 480 I=N+NP+1,N+2*NP IF(K(N+K(I,4),5).LE.1.OR.P(I,5).LT.RINIT) GOTO 480 IJET=N+K(I,4) R2=R2T(I,IJET) IF(R2.LE.R2MAX) GOTO 480 ISPL=I R2MAX=R2 480 CONTINUE IF(ISPL.NE.0) THEN IJET=N+K(ISPL,4) DO 490 J=1,4 P(IEMP,J)=P(ISPL,J) P(IJET,J)=P(IJET,J)-P(ISPL,J) 490 CONTINUE 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 300 ENDIF ENDIF ENDIF C...If generalized thrust has not yet converged, continue iteration. IF(MSTU(46).LE.1.AND.NLOOP.LE.2.AND.PSJT/PSS.GT.TSAV+PARU(48)) &THEN TSAV=PSJT/PSS GOTO 310 ENDIF C...Reorder jets according to energy. DO 510 I=N+1,N+NJET DO 500 J=1,5 V(I,J)=P(I,J) 500 CONTINUE 510 CONTINUE DO 540 INEW=N+1,N+NJET PEMAX=0D0 DO 520 ITRY=N+1,N+NJET IF(V(ITRY,4).LE.PEMAX) GOTO 520 IMAX=ITRY PEMAX=V(ITRY,4) 520 CONTINUE K(INEW,1)=31 K(INEW,2)=97 K(INEW,3)=INEW-N K(INEW,4)=0 DO 530 J=1,5 P(INEW,J)=V(IMAX,J) 530 CONTINUE V(IMAX,4)=-1D0 K(IMAX,5)=INEW 540 CONTINUE C...Clean up particle-jet assignments and jet information. DO 550 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 550 CONTINUE IEMP=0 PSJT=0D0 DO 570 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,0D0)) DO 560 J=1,5 V(I,J)=0D0 560 CONTINUE IF(K(I,4).EQ.0) IEMP=I 570 CONTINUE C...Select storing option. Output variables. Check for failure. MSTU(61)=N+1 MSTU(62)=NP MSTU(63)=NPRE PARU(61)=PS(5) PARU(62)=PSJT/PSS PARU(63)=SQRT(R2MIN) IF(NJET.LE.1) PARU(63)=0D0 IF(IEMP.NE.0) THEN CALL PYERRM(8,'(PYCLUS:) failed to reconstruct as requested') NJET=-1 RETURN ENDIF IF(MSTU(43).LE.1) MSTU(3)=MAX(0,NJET) IF(MSTU(43).GE.2) N=N+MAX(0,NJET) NSAV=NJET RETURN END