Update master to aliroot
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / luclus_hijing.F
1 * $Id$
2     
3 C*********************************************************************  
4     
5       SUBROUTINE LUCLUS_HIJING(NJET)   
6     
7 C...Purpose: to subdivide the particle content of an event into 
8 C...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     
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)))   
20     
21 C...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     
38 C...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     
55 C...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     
74 C...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     
82 C...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     
94 C...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     
110 C...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     
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    
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     
136 C...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     
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    
160       NJET=NPRE 
161     
162 C...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     
190 C...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     
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)   
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     
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    
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     
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))    
256      &THEN  
257         TSAV=PSJT/PSS   
258         GOTO 300    
259       ENDIF 
260     
261 C...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     
281 C...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     
298 C...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