]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hipyset1_35/luclus_hijing.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
-> Changes to triggering scheme for p-Pb
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / luclus_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE LUCLUS_HIJING(NJET)
6
7C...Purpose: to subdivide the particle content of an event into
8C...jets/clusters.
9#include "lujets_hijing.inc"
10#include "ludat1_hijing.inc"
11#include "ludat2_hijing.inc"
12 DIMENSION PS(5)
13 SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM
14
15C...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)))
20
21C...If first time, reset. If reentering, skip preliminaries.
22 IF(MSTU(48).LE.0) THEN
23 NP=0
24 DO 100 J=1,5
25 100 PS(J)=0.
26 PSS=0.
27 ELSE
28 NJET=NSAV
29 IF(MSTU(43).GE.2) N=N-NJET
30 DO 110 I=N+1,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
34 NLOOP=0
35 GOTO 290
36 ENDIF
37
38C...Find which particles are to be considered in cluster search.
39 DO 140 I=1,N
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.
44 & KC.EQ.18) GOTO 140
45 IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
46 $ .EQ.0)GOTO 140
47 ENDIF
48 IF(N+2*NP.GE.MSTU(4)-MSTU(32)-5) THEN
49 CALL LUERRM_HIJING(11
50 $ ,'(LUCLUS_HIJING:) no more memory left in LUJETS_HIJING')
51 NJET=-1
52 RETURN
53 ENDIF
54
55C...Take copy of these particles, with space left for jets later on.
56 NP=NP+1
57 K(N+NP,3)=I
58 DO 120 J=1,5
59 120 P(N+NP,J)=P(I,J)
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)
64 DO 130 J=1,4
65 130 PS(J)=PS(J)+P(N+NP,J)
66 PSS=PSS+P(N+NP,5)
67 140 CONTINUE
68 DO 150 I=N+1,N+NP
69 K(I+NP,3)=K(I,3)
70 DO 150 J=1,5
71 150 P(I+NP,J)=P(I,J)
72 PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
73
74C...Very low multiplicities not considered.
75 IF(NP.LT.MSTU(47)) THEN
76 CALL LUERRM_HIJING(8
77 $ ,'(LUCLUS_HIJING:) too few particles for analysis')
78 NJET=-1
79 RETURN
80 ENDIF
81
82C...Find precluster configuration. If too few jets, make harder cuts.
83 NLOOP=0
84 IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
85 IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
86 RINIT=1.25*PARU(43)
87 IF(NP.LE.MSTU(47)+2) RINIT=0.
88 160 RINIT=0.8*RINIT
89 NPRE=0
90 NREM=NP
91 DO 170 I=N+NP+1,N+2*NP
92 170 K(I,4)=0
93
94C...Sum up small momentum region. Jet if enough absolute momentum.
95 IF(MSTU(46).LE.2) THEN
96 DO 180 J=1,4
97 180 P(N+1,J)=0.
98 DO 200 I=N+NP+1,N+2*NP
99 IF(P(I,5).GT.2.*RINIT) GOTO 200
100 NREM=NREM-1
101 K(I,4)=1
102 DO 190 J=1,4
103 190 P(N+1,J)=P(N+1,J)+P(I,J)
104 200 CONTINUE
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
108 ENDIF
109
110C...Find fastest remaining particle.
111 210 NPRE=NPRE+1
112 PMAX=0.
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
115 IMAX=I
116 PMAX=P(I,5)
117 220 CONTINUE
118 DO 230 J=1,5
119 230 P(N+NPRE,J)=P(IMAX,J)
120 NREM=NREM-1
121 K(IMAX,4)=NPRE
122
123C...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
127 R2=R2T(I,IMAX)
128 IF(R2.GT.RINIT**2) GOTO 250
129 NREM=NREM-1
130 K(I,4)=NPRE
131 DO 240 J=1,4
132 240 P(N+NPRE,J)=P(N+NPRE,J)+P(I,J)
133 250 CONTINUE
134 P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
135
136C...Sum up precluster around it according to mass separation.
137 ELSE
138 260 IMIN=0
139 R2MIN=RINIT**2
140 DO 270 I=N+NP+1,N+2*NP
141 IF(K(I,4).NE.0) GOTO 270
142 R2=R2M(I,N+NPRE)
143 IF(R2.GE.R2MIN) GOTO 270
144 IMIN=I
145 R2MIN=R2
146 270 CONTINUE
147 IF(IMIN.NE.0) THEN
148 DO 280 J=1,4
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)
151 NREM=NREM-1
152 K(IMIN,4)=NPRE
153 GOTO 260
154 ENDIF
155 ENDIF
156
157C...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
160 NJET=NPRE
161
162C...Reassign all particles to nearest jet. Sum up new jet momenta.
163 290 TSAV=0.
164 PSJT=0.
165 300 IF(MSTU(46).LE.1) THEN
166 DO 310 I=N+1,N+NJET
167 DO 310 J=1,4
168 310 V(I,J)=0.
169 DO 340 I=N+NP+1,N+2*NP
170 R2MIN=PSS**2
171 DO 320 IJET=N+1,N+NJET
172 IF(P(IJET,5).LT.RINIT) GOTO 320
173 R2=R2T(I,IJET)
174 IF(R2.GE.R2MIN) GOTO 320
175 IMIN=IJET
176 R2MIN=R2
177 320 CONTINUE
178 K(I,4)=IMIN-N
179 DO 330 J=1,4
180 330 V(IMIN,J)=V(IMIN,J)+P(I,J)
181 340 CONTINUE
182 PSJT=0.
183 DO 360 I=N+1,N+NJET
184 DO 350 J=1,4
185 350 P(I,J)=V(I,J)
186 P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
187 360 PSJT=PSJT+P(I,5)
188 ENDIF
189
190C...Find two closest jets.
191 R2MIN=2.*R2ACC
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
197 IMIN1=ITRY1
198 IMIN2=ITRY2
199 R2MIN=R2
200 370 CONTINUE
201
202C...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)
206 DO 380 J=1,4
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
210 DO 390 J=1,5
211 390 P(I-1,J)=P(I,J)
212 IF(MSTU(46).GE.2) THEN
213 DO 400 I=N+NP+1,N+2*NP
214 IORI=N+K(I,4)
215 IF(IORI.EQ.IDEL) K(I,4)=IREC-N
216 400 IF(IORI.GT.IDEL) K(I,4)=K(I,4)-1
217 ENDIF
218 NJET=NJET-1
219 GOTO 290
220
221C...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
223 DO 410 I=N+1,N+NJET
224 410 K(I,5)=0
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
227 IEMP=0
228 DO 430 I=N+1,N+NJET
229 430 IF(K(I,5).EQ.0) IEMP=I
230 IF(IEMP.NE.0) THEN
231 NLOOP=NLOOP+1
232 ISPL=0
233 R2MAX=0.
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
236 IJET=N+K(I,4)
237 R2=R2T(I,IJET)
238 IF(R2.LE.R2MAX) GOTO 440
239 ISPL=I
240 R2MAX=R2
241 440 CONTINUE
242 IF(ISPL.NE.0) THEN
243 IJET=N+K(ISPL,4)
244 DO 450 J=1,4
245 P(IEMP,J)=P(ISPL,J)
246 450 P(IJET,J)=P(IJET,J)-P(ISPL,J)
247 P(IEMP,5)=P(ISPL,5)
248 P(IJET,5)=SQRT(P(IJET,1)**2+P(IJET,2)**2+P(IJET,3)**2)
249 IF(NLOOP.LE.2) GOTO 290
250 ENDIF
251 ENDIF
252 ENDIF
253
254C...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))
256 &THEN
257 TSAV=PSJT/PSS
258 GOTO 300
259 ENDIF
260
261C...Reorder jets according to energy.
262 DO 460 I=N+1,N+NJET
263 DO 460 J=1,5
264 460 V(I,J)=P(I,J)
265 DO 490 INEW=N+1,N+NJET
266 PEMAX=0.
267 DO 470 ITRY=N+1,N+NJET
268 IF(V(ITRY,4).LE.PEMAX) GOTO 470
269 IMAX=ITRY
270 PEMAX=V(ITRY,4)
271 470 CONTINUE
272 K(INEW,1)=31
273 K(INEW,2)=97
274 K(INEW,3)=INEW-N
275 K(INEW,4)=0
276 DO 480 J=1,5
277 480 P(INEW,J)=V(IMAX,J)
278 V(IMAX,4)=-1.
279 490 K(IMAX,5)=INEW
280
281C...Clean up particle-jet assignments and jet information.
282 DO 500 I=N+NP+1,N+2*NP
283 IORI=K(N+K(I,4),5)
284 K(I,4)=IORI-N
285 IF(K(K(I,3),1).NE.3) K(K(I,3),4)=IORI-N
286 K(IORI,4)=K(IORI,4)+1
287 500 CONTINUE
288 IEMP=0
289 PSJT=0.
290 DO 520 I=N+1,N+NJET
291 K(I,5)=0
292 PSJT=PSJT+P(I,5)
293 P(I,5)=SQRT(MAX(P(I,4)**2-P(I,5)**2,0.))
294 DO 510 J=1,5
295 510 V(I,J)=0.
296 520 IF(K(I,4).EQ.0) IEMP=I
297
298C...Select storing option. Output variables. Check for failure.
299 MSTU(61)=N+1
300 MSTU(62)=NP
301 MSTU(63)=NPRE
302 PARU(61)=PS(5)
303 PARU(62)=PSJT/PSS
304 PARU(63)=SQRT(R2MIN)
305 IF(NJET.LE.1) PARU(63)=0.
306 IF(IEMP.NE.0) THEN
307 CALL LUERRM_HIJING(8
308 $ ,'(LUCLUS_HIJING:) failed to reconstruct as requested')
309 NJET=-1
310 ENDIF
311 IF(MSTU(43).LE.1) MSTU(3)=NJET
312 IF(MSTU(43).GE.2) N=N+NJET
313 NSAV=NJET
314
315 RETURN
316 END