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