3 C*********************************************************************
5 SUBROUTINE LUCLUS_HIJING(NJET)
7 C...Purpose: to subdivide the particle content of an event into
9 #include "lujets_hijing.inc"
10 #include "ludat1_hijing.inc"
11 #include "ludat2_hijing.inc"
13 SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM
15 C...Functions: distance measure in pT or (pseudo)mass.
16 R2T(I1,I2)=(P(I1,5)*P(I2,5)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
17 &P(I1,3)*P(I2,3))*2.*P(I1,5)*P(I2,5)/(0.0001+P(I1,5)+P(I2,5))**2
18 R2M(I1,I2)=2.*P(I1,4)*P(I2,4)*(1.-(P(I1,1)*P(I2,1)+P(I1,2)*
19 &P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5)))
21 C...If first time, reset. If reentering, skip preliminaries.
22 IF(MSTU(48).LE.0) THEN
29 IF(MSTU(43).GE.2) N=N-NJET
31 110 P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
32 IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
33 IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
38 C...Find which particles are to be considered in cluster search.
40 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 140
41 IF(MSTU(41).GE.2) THEN
42 KC=LUCOMP_HIJING(K(I,2))
43 IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
45 IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
48 IF(N+2*NP.GE.MSTU(4)-MSTU(32)-5) THEN
50 $ ,'(LUCLUS_HIJING:) no more memory left in LUJETS_HIJING')
55 C...Take copy of these particles, with space left for jets later on.
60 IF(MSTU(42).EQ.0) P(N+NP,5)=0.
61 IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
62 P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
63 P(N+NP,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
65 130 PS(J)=PS(J)+P(N+NP,J)
72 PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
74 C...Very low multiplicities not considered.
75 IF(NP.LT.MSTU(47)) THEN
77 $ ,'(LUCLUS_HIJING:) too few particles for analysis')
82 C...Find precluster configuration. If too few jets, make harder cuts.
84 IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
85 IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
87 IF(NP.LE.MSTU(47)+2) RINIT=0.
91 DO 170 I=N+NP+1,N+2*NP
94 C...Sum up small momentum region. Jet if enough absolute momentum.
95 IF(MSTU(46).LE.2) THEN
98 DO 200 I=N+NP+1,N+2*NP
99 IF(P(I,5).GT.2.*RINIT) GOTO 200
103 190 P(N+1,J)=P(N+1,J)+P(I,J)
105 P(N+1,5)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)
106 IF(P(N+1,5).GT.2.*RINIT) NPRE=1
107 IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 160
110 C...Find fastest remaining particle.
113 DO 220 I=N+NP+1,N+2*NP
114 IF(K(I,4).NE.0.OR.P(I,5).LE.PMAX) GOTO 220
119 230 P(N+NPRE,J)=P(IMAX,J)
123 C...Sum up precluster around it according to pT separation.
124 IF(MSTU(46).LE.2) THEN
125 DO 250 I=N+NP+1,N+2*NP
126 IF(K(I,4).NE.0) GOTO 250
128 IF(R2.GT.RINIT**2) GOTO 250
132 240 P(N+NPRE,J)=P(N+NPRE,J)+P(I,J)
134 P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
136 C...Sum up precluster around it according to mass separation.
140 DO 270 I=N+NP+1,N+2*NP
141 IF(K(I,4).NE.0) GOTO 270
143 IF(R2.GE.R2MIN) GOTO 270
149 280 P(N+NPRE,J)=P(N+NPRE,J)+P(IMIN,J)
150 P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
157 C...Check if more preclusters to be found. Start over if too few.
158 IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 160
159 IF(NREM.GT.0) GOTO 210
162 C...Reassign all particles to nearest jet. Sum up new jet momenta.
165 300 IF(MSTU(46).LE.1) THEN
169 DO 340 I=N+NP+1,N+2*NP
171 DO 320 IJET=N+1,N+NJET
172 IF(P(IJET,5).LT.RINIT) GOTO 320
174 IF(R2.GE.R2MIN) GOTO 320
180 330 V(IMIN,J)=V(IMIN,J)+P(I,J)
186 P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
190 C...Find two closest jets.
192 DO 370 ITRY1=N+1,N+NJET-1
193 DO 370 ITRY2=ITRY1+1,N+NJET
194 IF(MSTU(46).LE.2) R2=R2T(ITRY1,ITRY2)
195 IF(MSTU(46).GE.3) R2=R2M(ITRY1,ITRY2)
196 IF(R2.GE.R2MIN) GOTO 370
202 C...If allowed, join two closest jets and start over.
203 IF(NJET.GT.MSTU(47).AND.R2MIN.LT.R2ACC) THEN
204 IREC=MIN(IMIN1,IMIN2)
205 IDEL=MAX(IMIN1,IMIN2)
207 380 P(IREC,J)=P(IMIN1,J)+P(IMIN2,J)
208 P(IREC,5)=SQRT(P(IREC,1)**2+P(IREC,2)**2+P(IREC,3)**2)
209 DO 390 I=IDEL+1,N+NJET
212 IF(MSTU(46).GE.2) THEN
213 DO 400 I=N+NP+1,N+2*NP
215 IF(IORI.EQ.IDEL) K(I,4)=IREC-N
216 400 IF(IORI.GT.IDEL) K(I,4)=K(I,4)-1
221 C...Divide up broad jet if empty cluster in list of final ones.
222 ELSEIF(NJET.EQ.MSTU(47).AND.MSTU(46).LE.1.AND.NLOOP.LE.2) THEN
225 DO 420 I=N+NP+1,N+2*NP
226 420 K(N+K(I,4),5)=K(N+K(I,4),5)+1
229 430 IF(K(I,5).EQ.0) IEMP=I
234 DO 440 I=N+NP+1,N+2*NP
235 IF(K(N+K(I,4),5).LE.1.OR.P(I,5).LT.RINIT) GOTO 440
238 IF(R2.LE.R2MAX) GOTO 440
246 450 P(IJET,J)=P(IJET,J)-P(ISPL,J)
248 P(IJET,5)=SQRT(P(IJET,1)**2+P(IJET,2)**2+P(IJET,3)**2)
249 IF(NLOOP.LE.2) GOTO 290
254 C...If generalized thrust has not yet converged, continue iteration.
255 IF(MSTU(46).LE.1.AND.NLOOP.LE.2.AND.PSJT/PSS.GT.TSAV+PARU(48))
261 C...Reorder jets according to energy.
265 DO 490 INEW=N+1,N+NJET
267 DO 470 ITRY=N+1,N+NJET
268 IF(V(ITRY,4).LE.PEMAX) GOTO 470
277 480 P(INEW,J)=V(IMAX,J)
281 C...Clean up particle-jet assignments and jet information.
282 DO 500 I=N+NP+1,N+2*NP
285 IF(K(K(I,3),1).NE.3) K(K(I,3),4)=IORI-N
286 K(IORI,4)=K(IORI,4)+1
293 P(I,5)=SQRT(MAX(P(I,4)**2-P(I,5)**2,0.))
296 520 IF(K(I,4).EQ.0) IEMP=I
298 C...Select storing option. Output variables. Check for failure.
305 IF(NJET.LE.1) PARU(63)=0.
308 $ ,'(LUCLUS_HIJING:) failed to reconstruct as requested')
311 IF(MSTU(43).LE.1) MSTU(3)=NJET
312 IF(MSTU(43).GE.2) N=N+NJET