--- /dev/null
+* $Id$
+
+C*********************************************************************
+
+ SUBROUTINE LUCLUS_HIJING(NJET)
+
+C...Purpose: to subdivide the particle content of an event into
+C...jets/clusters.
+#include "lujets_hijing.inc"
+#include "ludat1_hijing.inc"
+#include "ludat2_hijing.inc"
+ DIMENSION PS(5)
+ SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM
+
+C...Functions: distance measure in pT or (pseudo)mass.
+ 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))*2.*P(I1,5)*P(I2,5)/(0.0001+P(I1,5)+P(I2,5))**2
+ R2M(I1,I2)=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)))
+
+C...If first time, reset. If reentering, skip preliminaries.
+ IF(MSTU(48).LE.0) THEN
+ NP=0
+ DO 100 J=1,5
+ 100 PS(J)=0.
+ PSS=0.
+ ELSE
+ NJET=NSAV
+ IF(MSTU(43).GE.2) N=N-NJET
+ DO 110 I=N+1,N+NJET
+ 110 P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
+ IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
+ IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
+ NLOOP=0
+ GOTO 290
+ 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=LUCOMP_HIJING(K(I,2))
+ IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
+ & KC.EQ.18) GOTO 140
+ IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
+ $ .EQ.0)GOTO 140
+ ENDIF
+ IF(N+2*NP.GE.MSTU(4)-MSTU(32)-5) THEN
+ CALL LUERRM_HIJING(11
+ $ ,'(LUCLUS_HIJING:) no more memory left in LUJETS_HIJING')
+ 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
+ 120 P(N+NP,J)=P(I,J)
+ IF(MSTU(42).EQ.0) P(N+NP,5)=0.
+ IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
+ 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.MSTU(47)) THEN
+ CALL LUERRM_HIJING(8
+ $ ,'(LUCLUS_HIJING:) 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) R2ACC=PARU(44)**2
+ IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
+ RINIT=1.25*PARU(43)
+ IF(NP.LE.MSTU(47)+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(MSTU(46).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*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) 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(MSTU(46).LE.2) THEN
+ DO 250 I=N+NP+1,N+2*NP
+ IF(K(I,4).NE.0) GOTO 250
+ R2=R2T(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=R2M(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*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) 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(MSTU(46).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=R2T(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(MSTU(46).LE.2) R2=R2T(ITRY1,ITRY2)
+ IF(MSTU(46).GE.3) R2=R2M(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.MSTU(47).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(MSTU(46).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.MSTU(47).AND.MSTU(46).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=R2T(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(MSTU(46).LE.1.AND.NLOOP.LE.2.AND.PSJT/PSS.GT.TSAV+PARU(48))
+ &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.
+ 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)=0.
+ IF(IEMP.NE.0) THEN
+ CALL LUERRM_HIJING(8
+ $ ,'(LUCLUS_HIJING:) failed to reconstruct as requested')
+ NJET=-1
+ ENDIF
+ IF(MSTU(43).LE.1) MSTU(3)=NJET
+ IF(MSTU(43).GE.2) N=N+NJET
+ NSAV=NJET
+
+ RETURN
+ END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff