]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hipyset1_35/luindf_hijing.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Include Lambda(1520) decay
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / luindf_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE LUINDF_HIJING(IP)
6
7C...Purpose: to handle the fragmentation of a jet system (or a single
8C...jet) according to independent fragmentation models.
9 IMPLICIT DOUBLE PRECISION(D)
10#include "lujets_hijing.inc"
11#include "ludat1_hijing.inc"
12#include "ludat2_hijing.inc"
13 DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
14 &KFLO(2),PXO(2),PYO(2),WO(2)
15
16C...Reset counters. Identify parton system and take copy. Check flavour.
17 NSAV=N
18 NJET=0
19 KQSUM=0
20 DO 100 J=1,5
21 100 DPS(J)=0.
22 I=IP-1
23 110 I=I+1
24 IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
25 CALL LUERRM_HIJING(12
26 $ ,'(LUINDF_HIJING:) failed to reconstruct jet system')
27 IF(MSTU(21).GE.1) RETURN
28 ENDIF
29 IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
30 KC=LUCOMP_HIJING(K(I,2))
31 IF(KC.EQ.0) GOTO 110
32 KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
33 IF(KQ.EQ.0) GOTO 110
34 NJET=NJET+1
35 IF(KQ.NE.2) KQSUM=KQSUM+KQ
36 DO 120 J=1,5
37 K(NSAV+NJET,J)=K(I,J)
38 P(NSAV+NJET,J)=P(I,J)
39 120 DPS(J)=DPS(J)+P(I,J)
40 K(NSAV+NJET,3)=I
41 IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
42 &K(I+1,1).EQ.2)) GOTO 110
43 IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
44 CALL LUERRM_HIJING(12
45 $ ,'(LUINDF_HIJING:) unphysical flavour combination')
46 IF(MSTU(21).GE.1) RETURN
47 ENDIF
48
49C...Boost copied system to CM frame. Find CM energy and sum flavours.
50 IF(NJET.NE.1) CALL LUDBRB_HIJING(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)
51 $ /DPS(4),-DPS(2)/DPS(4),-DPS(3)/DPS(4))
52 PECM=0.
53 DO 130 J=1,3
54 130 NFI(J)=0
55 DO 140 I=NSAV+1,NSAV+NJET
56 PECM=PECM+P(I,4)
57 KFA=IABS(K(I,2))
58 IF(KFA.LE.3) THEN
59 NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
60 ELSEIF(KFA.GT.1000) THEN
61 KFLA=MOD(KFA/1000,10)
62 KFLB=MOD(KFA/100,10)
63 IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
64 IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
65 ENDIF
66 140 CONTINUE
67
68C...Loop over attempts made. Reset counters.
69 NTRY=0
70 150 NTRY=NTRY+1
71 N=NSAV+NJET
72 IF(NTRY.GT.200) THEN
73 CALL LUERRM_HIJING(14
74 $ ,'(LUINDF_HIJING:) caught in infinite loop')
75 IF(MSTU(21).GE.1) RETURN
76 ENDIF
77 DO 160 J=1,3
78 NFL(J)=NFI(J)
79 IFET(J)=0
80 160 KFLF(J)=0
81
82C...Loop over jets to be fragmented.
83 DO 230 IP1=NSAV+1,NSAV+NJET
84 MSTJ(91)=0
85 NSAV1=N
86
87C...Initial flavour and momentum values. Jet along +z axis.
88 KFLH=IABS(K(IP1,2))
89 IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
90 KFLO(2)=0
91 WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
92
93C...Initial values for quark or diquark jet.
94 170 IF(IABS(K(IP1,2)).NE.21) THEN
95 NSTR=1
96 KFLO(1)=K(IP1,2)
97 CALL LUPTDI_HIJING(0,PXO(1),PYO(1))
98 WO(1)=WF
99
100C...Initial values for gluon treated like random quark jet.
101 ELSEIF(MSTJ(2).LE.2) THEN
102 NSTR=1
103 IF(MSTJ(2).EQ.2) MSTJ(91)=1
104 KFLO(1)=INT(1.+(2.+PARJ(2))*RLU_HIJING(0))*(-1)
105 $ **INT(RLU_HIJING(0)+0.5)
106 CALL LUPTDI_HIJING(0,PXO(1),PYO(1))
107 WO(1)=WF
108
109C...Initial values for gluon treated like quark-antiquark jet pair,
110C...sharing energy according to Altarelli-Parisi splitting function.
111 ELSE
112 NSTR=2
113 IF(MSTJ(2).EQ.4) MSTJ(91)=1
114 KFLO(1)=INT(1.+(2.+PARJ(2))*RLU_HIJING(0))*(-1)
115 $ **INT(RLU_HIJING(0)+0.5)
116 KFLO(2)=-KFLO(1)
117 CALL LUPTDI_HIJING(0,PXO(1),PYO(1))
118 PXO(2)=-PXO(1)
119 PYO(2)=-PYO(1)
66f2039c 120 111 WO(1)=WF*RLU_HIJING(0)**(1./3.)
121 WO(2)=WF-WO(1)
122 IF (ABS(WO(2)).LE.0..OR.ABS(WO(1)).LE.0.) GOTO 111
e74335a4 123 ENDIF
124
125C...Initial values for rank, flavour, pT and W+.
126 DO 220 ISTR=1,NSTR
127 180 I=N
128 IRANK=0
129 KFL1=KFLO(ISTR)
130 PX1=PXO(ISTR)
131 PY1=PYO(ISTR)
132 W=WO(ISTR)
133
134C...New hadron. Generate flavour and hadron species.
135 190 I=I+1
136 IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
137 CALL LUERRM_HIJING(11
138 $ ,'(LUINDF_HIJING:) no more memory left in LUJETS_HIJING')
139 IF(MSTU(21).GE.1) RETURN
140 ENDIF
141 IRANK=IRANK+1
142 K(I,1)=1
143 K(I,3)=IP1
144 K(I,4)=0
145 K(I,5)=0
146 200 CALL LUKFDI_HIJING(KFL1,0,KFL2,K(I,2))
147 IF(K(I,2).EQ.0) GOTO 180
148 IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.
149 &IABS(KFL2).GT.10) THEN
150 IF(RLU_HIJING(0).GT.PARJ(19)) GOTO 200
151 ENDIF
152
153C...Find hadron mass. Generate four-momentum.
154 P(I,5)=ULMASS_HIJING(K(I,2))
155 CALL LUPTDI_HIJING(KFL1,PX2,PY2)
156 P(I,1)=PX1+PX2
157 P(I,2)=PY1+PY2
158 PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
159 CALL LUZDIS_HIJING(KFL1,KFL2,PR,Z)
160 P(I,3)=0.5*(Z*W-PR/(Z*W))
161 P(I,4)=0.5*(Z*W+PR/(Z*W))
162 IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
163 &P(I,3).LE.0.001) THEN
164 IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180
165 P(I,3)=0.0001
166 P(I,4)=SQRT(PR)
167 Z=P(I,4)/W
168 ENDIF
169
170C...Remaining flavour and momentum.
171 KFL1=-KFL2
172 PX1=-PX2
173 PY1=-PY2
174 W=(1.-Z)*W
175 DO 210 J=1,5
176 210 V(I,J)=0.
177
178C...Check if pL acceptable. Go back for new hadron if enough energy.
179 IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) I=I-1
180 IF(W.GT.PARJ(31)) GOTO 190
181 220 N=I
182 IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32)
183 IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
184
185C...Rotate jet to new direction.
186 THE=ULANGL_HIJING(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
187 PHI=ULANGL_HIJING(P(IP1,1),P(IP1,2))
188 CALL LUDBRB_HIJING(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
189 K(K(IP1,3),4)=NSAV1+1
190 K(K(IP1,3),5)=N
191
192C...End of jet generation loop. Skip conservation in some cases.
193 230 CONTINUE
194 IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 470
195 IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
196
197C...Subtract off produced hadron flavours, finished if zero.
198 DO 240 I=NSAV+NJET+1,N
199 KFA=IABS(K(I,2))
200 KFLA=MOD(KFA/1000,10)
201 KFLB=MOD(KFA/100,10)
202 KFLC=MOD(KFA/10,10)
203 IF(KFLA.EQ.0) THEN
204 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
205 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
206 ELSE
207 IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
208 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
209 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
210 ENDIF
211 240 CONTINUE
212 NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
213 &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
214 IF(NREQ.EQ.0) GOTO 320
215
216C...Take away flavour of low-momentum particles until enough freedom.
217 NREM=0
218 250 IREM=0
219 P2MIN=PECM**2
220 DO 260 I=NSAV+NJET+1,N
221 P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
222 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
223 260 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
224 IF(IREM.EQ.0) GOTO 150
225 K(IREM,1)=7
226 KFA=IABS(K(IREM,2))
227 KFLA=MOD(KFA/1000,10)
228 KFLB=MOD(KFA/100,10)
229 KFLC=MOD(KFA/10,10)
230 IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
231 IF(K(IREM,1).EQ.8) GOTO 250
232 IF(KFLA.EQ.0) THEN
233 ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
234 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
235 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
236 ELSE
237 IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
238 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
239 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
240 ENDIF
241 NREM=NREM+1
242 NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
243 &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
244 IF(NREQ.GT.NREM) GOTO 250
245 DO 270 I=NSAV+NJET+1,N
246 270 IF(K(I,1).EQ.8) K(I,1)=1
247
248C...Find combination of existing and new flavours for hadron.
249 280 NFET=2
250 IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
251 IF(NREQ.LT.NREM) NFET=1
252 IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
253 DO 290 J=1,NFET
254 IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU_HIJING(0)
255 KFLF(J)=ISIGN(1,NFL(1))
256 IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
257 290 IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
258 IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
259 &GOTO 280
260 IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
261 &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3).
262 &LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
263 IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU_HIJING(0))
264 IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
265 IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU_HIJING(0)),
266 $ -KFLF(1))
267 IF(NFET.LE.2) KFLF(3)=0
268 IF(KFLF(3).NE.0) THEN
269 KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
270 & 100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
271 IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU_HIJING(0).GT.1.)
272 & KFLFC=KFLFC+ISIGN(2,KFLFC)
273 ELSE
274 KFLFC=KFLF(1)
275 ENDIF
276 CALL LUKFDI_HIJING(KFLFC,KFLF(2),KFLDMP,KF)
277 IF(KF.EQ.0) GOTO 280
278 DO 300 J=1,MAX(2,NFET)
279 300 NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
280
281C...Store hadron at random among free positions.
282 NPOS=MIN(1+INT(RLU_HIJING(0)*NREM),NREM)
283 DO 310 I=NSAV+NJET+1,N
284 IF(K(I,1).EQ.7) NPOS=NPOS-1
285 IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
286 K(I,1)=1
287 K(I,2)=KF
288 P(I,5)=ULMASS_HIJING(K(I,2))
289 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
290 310 CONTINUE
291 NREM=NREM-1
292 NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
293 &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
294 IF(NREM.GT.0) GOTO 280
295
296C...Compensate for missing momentum in global scheme (3 options).
297 320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
298 DO 330 J=1,3
299 PSI(J)=0.
300 DO 330 I=NSAV+NJET+1,N
301 330 PSI(J)=PSI(J)+P(I,J)
302 PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
303 PWS=0.
304 DO 340 I=NSAV+NJET+1,N
305 IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
306 IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
307 & PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
308 340 IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1.
309 DO 360 I=NSAV+NJET+1,N
310 IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
311 IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
312 & PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
313 IF(MOD(MSTJ(3),5).EQ.3) PW=1.
314 DO 350 J=1,3
315 350 P(I,J)=P(I,J)-PSI(J)*PW/PWS
316 360 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
317
318C...Compensate for missing momentum withing each jet separately.
319 ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
320 DO 370 I=N+1,N+NJET
321 K(I,1)=0
322 DO 370 J=1,5
323 370 P(I,J)=0.
324 DO 390 I=NSAV+NJET+1,N
325 IR1=K(I,3)
326 IR2=N+IR1-NSAV
327 K(IR2,1)=K(IR2,1)+1
328 PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
329 & (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
330 DO 380 J=1,3
331 380 P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
332 P(IR2,4)=P(IR2,4)+P(I,4)
333 390 P(IR2,5)=P(IR2,5)+PLS
334 PSS=0.
335 DO 400 I=N+1,N+NJET
336 400 IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2))
337 DO 420 I=NSAV+NJET+1,N
338 IR1=K(I,3)
339 IR2=N+IR1-NSAV
340 PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
341 & (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
342 DO 410 J=1,3
343 410 P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS*
344 & P(IR1,J)
345 420 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
346 ENDIF
347
348C...Scale momenta for energy conservation.
349 IF(MOD(MSTJ(3),5).NE.0) THEN
350 PMS=0.
351 PES=0.
352 PQS=0.
353 DO 430 I=NSAV+NJET+1,N
354 PMS=PMS+P(I,5)
355 PES=PES+P(I,4)
356 430 PQS=PQS+P(I,5)**2/P(I,4)
357 IF(PMS.GE.PECM) GOTO 150
358 NECO=0
359 440 NECO=NECO+1
360 PFAC=(PECM-PQS)/(PES-PQS)
361 PES=0.
362 PQS=0.
363 DO 460 I=NSAV+NJET+1,N
364 DO 450 J=1,3
365 450 P(I,J)=PFAC*P(I,J)
366 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
367 PES=PES+P(I,4)
368 460 PQS=PQS+P(I,5)**2/P(I,4)
369 IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 440
370 ENDIF
371
372C...Origin of produced particles and parton daughter pointers.
373 470 DO 480 I=NSAV+NJET+1,N
374 IF(MSTU(16).NE.2) K(I,3)=NSAV+1
375 480 IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
376 DO 490 I=NSAV+1,NSAV+NJET
377 I1=K(I,3)
378 K(I1,1)=K(I1,1)+10
379 IF(MSTU(16).NE.2) THEN
380 K(I1,4)=NSAV+1
381 K(I1,5)=NSAV+1
382 ELSE
383 K(I1,4)=K(I1,4)-NJET+1
384 K(I1,5)=K(I1,5)-NJET+1
385 IF(K(I1,5).LT.K(I1,4)) THEN
386 K(I1,4)=0
387 K(I1,5)=0
388 ENDIF
389 ENDIF
390 490 CONTINUE
391
392C...Document independent fragmentation system. Remove copy of jets.
393 NSAV=NSAV+1
394 K(NSAV,1)=11
395 K(NSAV,2)=93
396 K(NSAV,3)=IP
397 K(NSAV,4)=NSAV+1
398 K(NSAV,5)=N-NJET+1
399 DO 500 J=1,4
400 P(NSAV,J)=DPS(J)
401 500 V(NSAV,J)=V(IP,J)
402 P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
403 V(NSAV,5)=0.
404 DO 510 I=NSAV+NJET,N
405 DO 510 J=1,5
406 K(I-NJET+1,J)=K(I,J)
407 P(I-NJET+1,J)=P(I,J)
408 510 V(I-NJET+1,J)=V(I,J)
409 N=N-NJET+1
410
411C...Boost back particle system. Set production vertices.
412 IF(NJET.NE.1) CALL LUDBRB_HIJING(NSAV+1,N,0.,0.,DPS(1)/DPS(4),
413 &DPS(2)/DPS(4),DPS(3)/DPS(4))
414 DO 520 I=NSAV+1,N
415 DO 520 J=1,4
416 520 V(I,J)=V(IP,J)
417
418 RETURN
419 END