b2a7675aa04370947900898836eee6503ac39dc4
[u/mrichter/AliRoot.git] / TAmpt / AMPT / hipyset1.35.f
1 c.................... hipyset1.35.f
2 C
3 C
4 C
5 C     Modified for HIJING program
6 c
7 c    modification July 22, 1997  In pyremnn put an upper limit
8 c     on the total pt kick the parton can accumulate via multiple
9 C     scattering. Set the upper limit to be the sqrt(s)/2,
10 c     this is fix cronin bug for Pb+Pb events at SPS energy.
11 c
12 C
13 C Last modification Oct. 1993 to comply with non-vax
14 C machines' compiler 
15 C
16 C*********************************************************************  
17     
18 cms
19 cms   gsfs 8/2009 Renamed common block PYINT4A due to conflict with something in CMSSW
20 cms 
21       SUBROUTINE LU2ENT(IP,KF1,KF2,PECM)    
22     
23 C...Purpose: to store two partons/particles in their CM frame,  
24 C...with the first along the +z axis.   
25       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
26       SAVE /LUJETSA/ 
27       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
28       SAVE /LUDAT1A/ 
29       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
30       SAVE /LUDAT2A/ 
31     
32 C...Standard checks.    
33       MSTU(28)=0    
34       IF(MSTU(12).GE.1) CALL LULIST(0)  
35       IPA=MAX(1,IABS(IP))   
36       IF(IPA.GT.MSTU(4)-1) CALL LUERRM(21,  
37      &'(LU2ENT:) writing outside LUJETSA memory')    
38       KC1=LUCOMP(KF1)   
39       KC2=LUCOMP(KF2)   
40       IF(KC1.EQ.0.OR.KC2.EQ.0) CALL LUERRM(12,  
41      &'(LU2ENT:) unknown flavour code') 
42     
43 C...Find masses. Reset K, P and V vectors.  
44       PM1=0.    
45       IF(MSTU(10).EQ.1) PM1=P(IPA,5)    
46       IF(MSTU(10).GE.2) PM1=ULMASS(KF1) 
47       PM2=0.    
48       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)  
49       IF(MSTU(10).GE.2) PM2=ULMASS(KF2) 
50       DO 100 I=IPA,IPA+1    
51       DO 100 J=1,5  
52       K(I,J)=0  
53       P(I,J)=0. 
54   100 V(I,J)=0. 
55     
56 C...Check flavours. 
57       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)  
58       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)  
59       IF(KQ1+KQ2.NE.0.AND.KQ1+KQ2.NE.4) CALL LUERRM(2,  
60      &'(LU2ENT:) unphysical flavour combination')   
61       K(IPA,2)=KF1  
62       K(IPA+1,2)=KF2    
63     
64 C...Store partons/particles in K vectors for normal case.   
65       IF(IP.GE.0) THEN  
66         K(IPA,1)=1  
67         IF(KQ1.NE.0.AND.KQ2.NE.0) K(IPA,1)=2    
68         K(IPA+1,1)=1    
69     
70 C...Store partons in K vectors for parton shower evolution. 
71       ELSE  
72         IF(KQ1.EQ.0.OR.KQ2.EQ.0) CALL LUERRM(2, 
73      &  '(LU2ENT:) requested flavours can not develop parton shower')   
74         K(IPA,1)=3  
75         K(IPA+1,1)=3    
76         K(IPA,4)=MSTU(5)*(IPA+1)    
77         K(IPA,5)=K(IPA,4)   
78         K(IPA+1,4)=MSTU(5)*IPA  
79         K(IPA+1,5)=K(IPA+1,4)   
80       ENDIF 
81     
82 C...Check kinematics and store partons/particles in P vectors.  
83       IF(PECM.LE.PM1+PM2) CALL LUERRM(13,   
84      &'(LU2ENT:) energy smaller than sum of masses')    
85       PA=SQRT(MAX(0.,(PECM**2-PM1**2-PM2**2)**2-(2.*PM1*PM2)**2))/  
86      &(2.*PECM) 
87       P(IPA,3)=PA   
88       P(IPA,4)=SQRT(PM1**2+PA**2)   
89       P(IPA,5)=PM1  
90       P(IPA+1,3)=-PA    
91       P(IPA+1,4)=SQRT(PM2**2+PA**2) 
92       P(IPA+1,5)=PM2    
93     
94 C...Set N. Optionally fragment/decay.   
95       N=IPA+1   
96       IF(IP.EQ.0) CALL LUEXEC   
97     
98       RETURN    
99       END   
100     
101 C*********************************************************************  
102     
103       SUBROUTINE LUGIVE(CHIN)   
104     
105 C...Purpose: to set values of commonblock variables.    
106       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
107       SAVE /LUJETSA/ 
108       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
109       SAVE /LUDAT1A/ 
110       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
111       SAVE /LUDAT2A/ 
112       COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)    
113       SAVE /LUDAT3A/ 
114       COMMON/LUDAT4A/CHAF(500)   
115       CHARACTER CHAF*8  
116       SAVE /LUDAT4A/ 
117       CHARACTER CHIN*(*),CHFIX*104,CHBIT*104,CHOLD*8,CHNEW*8,   
118      &CHNAM*4,CHVAR(17)*4,CHALP(2)*26,CHIND*8,CHINI*10,CHINR*16 
119       DATA CHVAR/'N','K','P','V','MSTU','PARU','MSTJ','PARJ','KCHG',    
120      &'PMAS','PARF','VCKM','MDCY','MDME','BRAT','KFDP','CHAF'/  
121       DATA CHALP/'abcdefghijklmnopqrstuvwxyz',  
122      &'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ 
123     
124 C...Length of character variable. Subdivide it into instructions.   
125       IF(MSTU(12).GE.1) CALL LULIST(0)  
126       CHBIT=CHIN//' '   
127       LBIT=101  
128   100 LBIT=LBIT-1   
129       IF(CHBIT(LBIT:LBIT).EQ.' ') GOTO 100  
130       LTOT=0    
131       DO 110 LCOM=1,LBIT    
132       IF(CHBIT(LCOM:LCOM).EQ.' ') GOTO 110  
133       LTOT=LTOT+1   
134       CHFIX(LTOT:LTOT)=CHBIT(LCOM:LCOM) 
135   110 CONTINUE  
136       LLOW=0    
137   120 LHIG=LLOW+1   
138   130 LHIG=LHIG+1   
139       IF(LHIG.LE.LTOT.AND.CHFIX(LHIG:LHIG).NE.';') GOTO 130 
140       LBIT=LHIG-LLOW-1  
141       CHBIT(1:LBIT)=CHFIX(LLOW+1:LHIG-1)    
142     
143 C...Identify commonblock variable.  
144       LNAM=1    
145   140 LNAM=LNAM+1   
146       IF(CHBIT(LNAM:LNAM).NE.'('.AND.CHBIT(LNAM:LNAM).NE.'='.AND.   
147      &LNAM.LE.4) GOTO 140   
148       CHNAM=CHBIT(1:LNAM-1)//' '    
149       DO 150 LCOM=1,LNAM-1  
150       DO 150 LALP=1,26  
151   150 IF(CHNAM(LCOM:LCOM).EQ.CHALP(1)(LALP:LALP)) CHNAM(LCOM:LCOM)= 
152      &CHALP(2)(LALP:LALP)   
153       IVAR=0    
154       DO 160 IV=1,17    
155   160 IF(CHNAM.EQ.CHVAR(IV)) IVAR=IV    
156       IF(IVAR.EQ.0) THEN    
157         CALL LUERRM(18,'(LUGIVE:) do not recognize variable '//CHNAM)   
158         LLOW=LHIG   
159         IF(LLOW.LT.LTOT) GOTO 120   
160         RETURN  
161       ENDIF 
162     
163 C...Identify any indices.   
164       I=0   
165       J=0   
166       IF(CHBIT(LNAM:LNAM).EQ.'(') THEN  
167         LIND=LNAM   
168   170   LIND=LIND+1 
169         IF(CHBIT(LIND:LIND).NE.')'.AND.CHBIT(LIND:LIND).NE.',') GOTO 170    
170         CHIND=' '   
171         IF((CHBIT(LNAM+1:LNAM+1).EQ.'C'.OR.CHBIT(LNAM+1:LNAM+1).EQ.'c').    
172      &  AND.(IVAR.EQ.9.OR.IVAR.EQ.10.OR.IVAR.EQ.13.OR.IVAR.EQ.17)) THEN 
173           CHIND(LNAM-LIND+11:8)=CHBIT(LNAM+2:LIND-1)    
174           READ(CHIND,'(I8)') I1 
175           I=LUCOMP(I1)  
176         ELSE    
177           CHIND(LNAM-LIND+10:8)=CHBIT(LNAM+1:LIND-1)    
178           READ(CHIND,'(I8)') I  
179         ENDIF   
180         LNAM=LIND   
181         IF(CHBIT(LNAM:LNAM).EQ.')') LNAM=LNAM+1 
182       ENDIF 
183       IF(CHBIT(LNAM:LNAM).EQ.',') THEN  
184         LIND=LNAM   
185   180   LIND=LIND+1 
186         IF(CHBIT(LIND:LIND).NE.')'.AND.CHBIT(LIND:LIND).NE.',') GOTO 180    
187         CHIND=' '   
188         CHIND(LNAM-LIND+10:8)=CHBIT(LNAM+1:LIND-1)  
189         READ(CHIND,'(I8)') J    
190         LNAM=LIND+1 
191       ENDIF 
192     
193 C...Check that indices allowed and save old value.  
194       IERR=1    
195       IF(CHBIT(LNAM:LNAM).NE.'=') GOTO 190  
196       IF(IVAR.EQ.1) THEN    
197         IF(I.NE.0.OR.J.NE.0) GOTO 190   
198         IOLD=N  
199       ELSEIF(IVAR.EQ.2) THEN    
200         IF(I.LT.1.OR.I.GT.MSTU(4).OR.J.LT.1.OR.J.GT.5) GOTO 190 
201         IOLD=K(I,J) 
202       ELSEIF(IVAR.EQ.3) THEN    
203         IF(I.LT.1.OR.I.GT.MSTU(4).OR.J.LT.1.OR.J.GT.5) GOTO 190 
204         ROLD=P(I,J) 
205       ELSEIF(IVAR.EQ.4) THEN    
206         IF(I.LT.1.OR.I.GT.MSTU(4).OR.J.LT.1.OR.J.GT.5) GOTO 190 
207         ROLD=V(I,J) 
208       ELSEIF(IVAR.EQ.5) THEN    
209         IF(I.LT.1.OR.I.GT.200.OR.J.NE.0) GOTO 190   
210         IOLD=MSTU(I)    
211       ELSEIF(IVAR.EQ.6) THEN    
212         IF(I.LT.1.OR.I.GT.200.OR.J.NE.0) GOTO 190   
213         ROLD=PARU(I)    
214       ELSEIF(IVAR.EQ.7) THEN    
215         IF(I.LT.1.OR.I.GT.200.OR.J.NE.0) GOTO 190   
216         IOLD=MSTJ(I)    
217       ELSEIF(IVAR.EQ.8) THEN    
218         IF(I.LT.1.OR.I.GT.200.OR.J.NE.0) GOTO 190   
219         ROLD=PARJ(I)    
220       ELSEIF(IVAR.EQ.9) THEN    
221         IF(I.LT.1.OR.I.GT.MSTU(6).OR.J.LT.1.OR.J.GT.3) GOTO 190 
222         IOLD=KCHG(I,J)  
223       ELSEIF(IVAR.EQ.10) THEN   
224         IF(I.LT.1.OR.I.GT.MSTU(6).OR.J.LT.1.OR.J.GT.4) GOTO 190 
225         ROLD=PMAS(I,J)  
226       ELSEIF(IVAR.EQ.11) THEN   
227         IF(I.LT.1.OR.I.GT.2000.OR.J.NE.0) GOTO 190  
228         ROLD=PARF(I)    
229       ELSEIF(IVAR.EQ.12) THEN   
230         IF(I.LT.1.OR.I.GT.4.OR.J.LT.1.OR.J.GT.4) GOTO 190   
231         ROLD=VCKM(I,J)  
232       ELSEIF(IVAR.EQ.13) THEN   
233         IF(I.LT.1.OR.I.GT.MSTU(6).OR.J.LT.1.OR.J.GT.3) GOTO 190 
234         IOLD=MDCY(I,J)  
235       ELSEIF(IVAR.EQ.14) THEN   
236         IF(I.LT.1.OR.I.GT.MSTU(7).OR.J.LT.1.OR.J.GT.2) GOTO 190 
237         IOLD=MDME(I,J)  
238       ELSEIF(IVAR.EQ.15) THEN   
239         IF(I.LT.1.OR.I.GT.MSTU(7).OR.J.NE.0) GOTO 190   
240         ROLD=BRAT(I)    
241       ELSEIF(IVAR.EQ.16) THEN   
242         IF(I.LT.1.OR.I.GT.MSTU(7).OR.J.LT.1.OR.J.GT.5) GOTO 190 
243         IOLD=KFDP(I,J)  
244       ELSEIF(IVAR.EQ.17) THEN   
245         IF(I.LT.1.OR.I.GT.MSTU(6).OR.J.NE.0) GOTO 190   
246         CHOLD=CHAF(I)   
247       ENDIF 
248       IERR=0    
249   190 IF(IERR.EQ.1) THEN    
250         CALL LUERRM(18,'(LUGIVE:) unallowed indices for '// 
251      &  CHBIT(1:LNAM-1))    
252         LLOW=LHIG   
253         IF(LLOW.LT.LTOT) GOTO 120   
254         RETURN  
255       ENDIF 
256     
257 C...Print current value of variable. Loop back. 
258       IF(LNAM.GE.LBIT) THEN 
259         CHBIT(LNAM:14)=' '  
260         CHBIT(15:60)=' has the value                                '   
261         IF(IVAR.EQ.1.OR.IVAR.EQ.2.OR.IVAR.EQ.5.OR.IVAR.EQ.7.OR. 
262      &  IVAR.EQ.9.OR.IVAR.EQ.13.OR.IVAR.EQ.14.OR.IVAR.EQ.16) THEN   
263           WRITE(CHBIT(51:60),'(I10)') IOLD  
264         ELSEIF(IVAR.NE.17) THEN 
265           WRITE(CHBIT(47:60),'(F14.5)') ROLD    
266         ELSE    
267           CHBIT(53:60)=CHOLD    
268         ENDIF   
269         IF(MSTU(13).GE.1) WRITE(MSTU(11),1000) CHBIT(1:60)  
270         LLOW=LHIG   
271         IF(LLOW.LT.LTOT) GOTO 120   
272         RETURN  
273       ENDIF 
274     
275 C...Read in new variable value. 
276       IF(IVAR.EQ.1.OR.IVAR.EQ.2.OR.IVAR.EQ.5.OR.IVAR.EQ.7.OR.   
277      &IVAR.EQ.9.OR.IVAR.EQ.13.OR.IVAR.EQ.14.OR.IVAR.EQ.16) THEN 
278         CHINI=' '   
279         CHINI(LNAM-LBIT+11:10)=CHBIT(LNAM+1:LBIT)   
280         READ(CHINI,'(I10)') INEW    
281       ELSEIF(IVAR.NE.17) THEN   
282         CHINR=' '   
283         CHINR(LNAM-LBIT+17:16)=CHBIT(LNAM+1:LBIT)   
284         READ(CHINR,'(F16.2)') RNEW  
285       ELSE  
286         CHNEW=CHBIT(LNAM+1:LBIT)//' '   
287       ENDIF 
288     
289 C...Store new variable value.   
290       IF(IVAR.EQ.1) THEN    
291         N=INEW  
292       ELSEIF(IVAR.EQ.2) THEN    
293         K(I,J)=INEW 
294       ELSEIF(IVAR.EQ.3) THEN    
295         P(I,J)=RNEW 
296       ELSEIF(IVAR.EQ.4) THEN    
297         V(I,J)=RNEW 
298       ELSEIF(IVAR.EQ.5) THEN    
299         MSTU(I)=INEW    
300       ELSEIF(IVAR.EQ.6) THEN    
301         PARU(I)=RNEW    
302       ELSEIF(IVAR.EQ.7) THEN    
303         MSTJ(I)=INEW    
304       ELSEIF(IVAR.EQ.8) THEN    
305         PARJ(I)=RNEW    
306       ELSEIF(IVAR.EQ.9) THEN    
307         KCHG(I,J)=INEW  
308       ELSEIF(IVAR.EQ.10) THEN   
309         PMAS(I,J)=RNEW  
310       ELSEIF(IVAR.EQ.11) THEN   
311         PARF(I)=RNEW    
312       ELSEIF(IVAR.EQ.12) THEN   
313         VCKM(I,J)=RNEW  
314       ELSEIF(IVAR.EQ.13) THEN   
315         MDCY(I,J)=INEW  
316       ELSEIF(IVAR.EQ.14) THEN   
317         MDME(I,J)=INEW  
318       ELSEIF(IVAR.EQ.15) THEN   
319         BRAT(I)=RNEW    
320       ELSEIF(IVAR.EQ.16) THEN   
321         KFDP(I,J)=INEW  
322       ELSEIF(IVAR.EQ.17) THEN   
323         CHAF(I)=CHNEW   
324       ENDIF 
325     
326 C...Write old and new value. Loop back. 
327       CHBIT(LNAM:14)=' '    
328       CHBIT(15:60)=' changed from                to               ' 
329       IF(IVAR.EQ.1.OR.IVAR.EQ.2.OR.IVAR.EQ.5.OR.IVAR.EQ.7.OR.   
330      &IVAR.EQ.9.OR.IVAR.EQ.13.OR.IVAR.EQ.14.OR.IVAR.EQ.16) THEN 
331         WRITE(CHBIT(33:42),'(I10)') IOLD    
332         WRITE(CHBIT(51:60),'(I10)') INEW    
333       ELSEIF(IVAR.NE.17) THEN   
334         WRITE(CHBIT(29:42),'(F14.5)') ROLD  
335         WRITE(CHBIT(47:60),'(F14.5)') RNEW  
336       ELSE  
337         CHBIT(35:42)=CHOLD  
338         CHBIT(53:60)=CHNEW  
339       ENDIF 
340       IF(MSTU(13).GE.1) WRITE(MSTU(11),1000) CHBIT(1:60)    
341       LLOW=LHIG 
342       IF(LLOW.LT.LTOT) GOTO 120 
343     
344 C...Format statement for output on unit MSTU(11) (by default 6).    
345  1000 FORMAT(5X,A60)    
346     
347       RETURN    
348       END   
349     
350 C*********************************************************************  
351     
352       SUBROUTINE LUEXEC 
353     
354 C...Purpose: to administrate the fragmentation and decay chain. 
355       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
356       SAVE /LUJETSA/ 
357       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
358       SAVE /LUDAT1A/ 
359       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
360       SAVE /LUDAT2A/ 
361       COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)    
362       SAVE /LUDAT3A/ 
363       DIMENSION PS(2,6) 
364     
365 C...Initialize and reset.   
366       MSTU(24)=0    
367       IF(MSTU(12).GE.1) CALL LULIST(0)  
368       MSTU(31)=MSTU(31)+1   
369       MSTU(1)=0 
370       MSTU(2)=0 
371       MSTU(3)=0 
372       MCONS=1   
373     
374 C...Sum up momentum, energy and charge for starting entries.    
375       NSAV=N    
376       DO 100 I=1,2  
377       DO 100 J=1,6  
378   100 PS(I,J)=0.    
379       DO 120 I=1,N  
380       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120  
381       DO 110 J=1,4  
382   110 PS(1,J)=PS(1,J)+P(I,J)    
383       PS(1,6)=PS(1,6)+LUCHGE(K(I,2))    
384   120 CONTINUE  
385       PARU(21)=PS(1,4)  
386     
387 C...Prepare system for subsequent fragmentation/decay.  
388       CALL LUPREP(0)    
389     
390 C...Loop through jet fragmentation and particle decays. 
391       MBE=0 
392   130 MBE=MBE+1 
393       IP=0  
394   140 IP=IP+1   
395       KC=0  
396       IF(K(IP,1).GT.0.AND.K(IP,1).LE.10) KC=LUCOMP(K(IP,2)) 
397       IF(KC.EQ.0) THEN  
398     
399 C...Particle decay if unstable and allowed. Save long-lived particle    
400 C...decays until second pass after Bose-Einstein effects.   
401       ELSEIF(KCHG(KC,2).EQ.0) THEN  
402 clin-4/2008 break up compound IF statements:
403 c        IF(MSTJ(21).GE.1.AND.MDCY(KC,1).GE.1.AND.(MSTJ(51).LE.0.OR.MBE. 
404 c     &  EQ.2.OR.PMAS(KC,2).GE.PARJ(91).OR.IABS(K(IP,2)).EQ.311))    
405 c     &  CALL LUDECY(IP) 
406          if(MSTJ(21).GE.1.AND.MDCY(KC,1).GE.1) then
407             if(MSTJ(51).LE.0.OR.MBE.EQ.2.OR.PMAS(KC,2).GE.PARJ(91)
408      &           .OR.IABS(K(IP,2)).EQ.311)
409      &           CALL LUDECY(IP) 
410          endif
411 c    
412 C...Decay products may develop a shower.    
413         IF(MSTJ(92).GT.0) THEN  
414           IP1=MSTJ(92)  
415           QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,  
416      &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))    
417           CALL LUSHOW(IP1,IP1+1,QMAX)   
418           CALL LUPREP(IP1)  
419           MSTJ(92)=0    
420         ELSEIF(MSTJ(92).LT.0) THEN  
421           IP1=-MSTJ(92) 
422 clin-8/19/02 avoid actual argument in common blocks of LUSHOW:
423 c          CALL LUSHOW(IP1,-3,P(IP,5))   
424           pip5=P(IP,5)
425           CALL LUSHOW(IP1,-3,pip5)   
426           CALL LUPREP(IP1)  
427           MSTJ(92)=0    
428         ENDIF   
429     
430 C...Jet fragmentation: string or independent fragmentation. 
431       ELSEIF(K(IP,1).EQ.1.OR.K(IP,1).EQ.2) THEN 
432         MFRAG=MSTJ(1)   
433         IF(MFRAG.GE.1.AND.K(IP,1).EQ.1) MFRAG=2 
434         IF(MSTJ(21).GE.2.AND.K(IP,1).EQ.2.AND.N.GT.IP) THEN 
435           IF(K(IP+1,1).EQ.1.AND.K(IP+1,3).EQ.K(IP,3).AND.   
436      &    K(IP,3).GT.0.AND.K(IP,3).LT.IP) THEN  
437             IF(KCHG(LUCOMP(K(K(IP,3),2)),2).EQ.0) MFRAG=MIN(1,MFRAG)    
438           ENDIF 
439         ENDIF   
440         IF(MFRAG.EQ.1) then
441            CALL LUSTRF(IP)  
442         endif
443         IF(MFRAG.EQ.2) CALL LUINDF(IP)  
444         IF(MFRAG.EQ.2.AND.K(IP,1).EQ.1) MCONS=0 
445         IF(MFRAG.EQ.2.AND.(MSTJ(3).LE.0.OR.MOD(MSTJ(3),5).EQ.0)) MCONS=0    
446       ENDIF 
447     
448 C...Loop back if enough space left in LUJETSA and no error abort.    
449       IF(MSTU(24).NE.0.AND.MSTU(21).GE.2) THEN  
450       ELSEIF(IP.LT.N.AND.N.LT.MSTU(4)-20-MSTU(32)) THEN 
451         GOTO 140    
452       ELSEIF(IP.LT.N) THEN  
453         CALL LUERRM(11,'(LUEXEC:) no more memory left in LUJETSA')   
454       ENDIF 
455     
456 C...Include simple Bose-Einstein effect parametrization if desired. 
457       IF(MBE.EQ.1.AND.MSTJ(51).GE.1) THEN   
458         CALL LUBOEI(NSAV)   
459         GOTO 130    
460       ENDIF 
461     
462 C...Check that momentum, energy and charge were conserved.  
463       DO 160 I=1,N  
464       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 160  
465       DO 150 J=1,4  
466   150 PS(2,J)=PS(2,J)+P(I,J)    
467       PS(2,6)=PS(2,6)+LUCHGE(K(I,2))    
468   160 CONTINUE  
469       PDEV=(ABS(PS(2,1)-PS(1,1))+ABS(PS(2,2)-PS(1,2))+ABS(PS(2,3)-  
470      &PS(1,3))+ABS(PS(2,4)-PS(1,4)))/(1.+ABS(PS(2,4))+ABS(PS(1,4))) 
471       IF(MCONS.EQ.1.AND.PDEV.GT.PARU(11)) CALL LUERRM(15,   
472      &'(LUEXEC:) four-momentum was not conserved')  
473 c      IF(MCONS.EQ.1.AND.PDEV.GT.PARU(11)) then
474 c         CALL LUERRM(15,   
475 c     &'(LUEXEC:) four-momentum was not conserved')  
476 c         write(6,*) 'PS1,2=',PS(1,1),PS(1,2),PS(1,3),PS(1,4),
477 c     1        '*',PS(2,1),PS(2,2),PS(2,3),PS(2,4)
478 c      endif
479
480       IF(MCONS.EQ.1.AND.ABS(PS(2,6)-PS(1,6)).GT.0.1) CALL LUERRM(15,    
481      &'(LUEXEC:) charge was not conserved') 
482     
483       RETURN    
484       END   
485     
486 C*********************************************************************  
487     
488       SUBROUTINE LUPREP(IP) 
489     
490 C...Purpose: to rearrange partons along strings, to allow small systems 
491 C...to collapse into one or two particles and to check flavours.    
492       IMPLICIT DOUBLE PRECISION(D)  
493       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
494       SAVE /LUJETSA/ 
495       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
496       SAVE /LUDAT1A/ 
497       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
498       SAVE /LUDAT2A/ 
499       COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)    
500       SAVE /LUDAT3A/ 
501       DIMENSION DPS(5),DPC(5),UE(3) 
502     
503       ic1=0
504       ic2=0
505       kci=0
506 C...Rearrange parton shower product listing along strings: begin loop.  
507       I1=N  
508       DO 130 MQGST=1,2  
509       DO 120 I=MAX(1,IP),N  
510       IF(K(I,1).NE.3) GOTO 120  
511       KC=LUCOMP(K(I,2)) 
512       IF(KC.EQ.0) GOTO 120  
513       KQ=KCHG(KC,2) 
514       IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 120  
515     
516 C...Pick up loose string end.   
517       KCS=4 
518       IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5 
519       IA=I  
520       NSTP=0    
521   100 NSTP=NSTP+1   
522       IF(NSTP.GT.4*N) THEN  
523         CALL LUERRM(14,'(LUPREP:) caught in infinite loop') 
524         RETURN  
525       ENDIF 
526     
527 C...Copy undecayed parton.  
528       IF(K(IA,1).EQ.3) THEN 
529         IF(I1.GE.MSTU(4)-MSTU(32)-5) THEN   
530           CALL LUERRM(11,'(LUPREP:) no more memory left in LUJETSA') 
531           RETURN    
532         ENDIF   
533         I1=I1+1 
534         K(I1,1)=2   
535         IF(NSTP.GE.2.AND.IABS(K(IA,2)).NE.21) K(I1,1)=1 
536         K(I1,2)=K(IA,2) 
537         K(I1,3)=IA  
538         K(I1,4)=0   
539         K(I1,5)=0   
540         DO 110 J=1,5    
541         P(I1,J)=P(IA,J) 
542   110   V(I1,J)=V(IA,J) 
543         K(IA,1)=K(IA,1)+10  
544         IF(K(I1,1).EQ.1) GOTO 120   
545       ENDIF 
546     
547 C...Go to next parton in colour space.  
548       IB=IA 
549       IF(MOD(K(IB,KCS)/MSTU(5)**2,2).EQ.0.AND.MOD(K(IB,KCS),MSTU(5)).   
550      &NE.0) THEN    
551         IA=MOD(K(IB,KCS),MSTU(5))   
552         K(IB,KCS)=K(IB,KCS)+MSTU(5)**2  
553         MREV=0  
554       ELSE  
555         IF(K(IB,KCS).GE.2*MSTU(5)**2.OR.MOD(K(IB,KCS)/MSTU(5),MSTU(5)). 
556      &  EQ.0) KCS=9-KCS 
557         IA=MOD(K(IB,KCS)/MSTU(5),MSTU(5))   
558         K(IB,KCS)=K(IB,KCS)+2*MSTU(5)**2    
559         MREV=1  
560       ENDIF 
561       IF(IA.LE.0.OR.IA.GT.N) THEN   
562         CALL LUERRM(12,'(LUPREP:) colour rearrangement failed') 
563         RETURN  
564       ENDIF 
565       IF(MOD(K(IA,4)/MSTU(5),MSTU(5)).EQ.IB.OR.MOD(K(IA,5)/MSTU(5), 
566      &MSTU(5)).EQ.IB) THEN  
567         IF(MREV.EQ.1) KCS=9-KCS 
568         IF(MOD(K(IA,KCS)/MSTU(5),MSTU(5)).NE.IB) KCS=9-KCS  
569         K(IA,KCS)=K(IA,KCS)+2*MSTU(5)**2    
570       ELSE  
571         IF(MREV.EQ.0) KCS=9-KCS 
572         IF(MOD(K(IA,KCS),MSTU(5)).NE.IB) KCS=9-KCS  
573         K(IA,KCS)=K(IA,KCS)+MSTU(5)**2  
574       ENDIF 
575       IF(IA.NE.I) GOTO 100  
576       K(I1,1)=1 
577   120 CONTINUE  
578   130 CONTINUE  
579       N=I1  
580     
581 C...Find lowest-mass colour singlet jet system, OK if above thresh.  
582       IF(MSTJ(14).LE.0) GOTO 320    
583       NS=N  
584   140 NSIN=N-NS 
585       PDM=1.+PARJ(32)   
586       IC=0  
587       DO 190 I=MAX(1,IP),NS 
588       IF(K(I,1).NE.1.AND.K(I,1).NE.2) THEN  
589       ELSEIF(K(I,1).EQ.2.AND.IC.EQ.0) THEN  
590         NSIN=NSIN+1 
591         IC=I    
592         DO 150 J=1,4    
593   150   DPS(J)=dble(P(I,J))
594         MSTJ(93)=1  
595         DPS(5)=dble(ULMASS(K(I,2)))
596       ELSEIF(K(I,1).EQ.2) THEN  
597         DO 160 J=1,4    
598   160   DPS(J)=DPS(J)+dble(P(I,J))
599       ELSEIF(IC.NE.0.AND.KCHG(LUCOMP(K(I,2)),2).NE.0) THEN  
600         DO 170 J=1,4    
601   170   DPS(J)=DPS(J)+dble(P(I,J))
602         MSTJ(93)=1  
603         DPS(5)=DPS(5)+dble(ULMASS(K(I,2)))
604         PD=sngl(SQRT(MAX(0D0,DPS(4)**2
605      1       -DPS(1)**2-DPS(2)**2-DPS(3)**2))-DPS(5))    
606         IF(PD.LT.PDM) THEN  
607           PDM=PD    
608           DO 180 J=1,5  
609   180     DPC(J)=DPS(J) 
610           IC1=IC    
611           IC2=I 
612         ENDIF   
613         IC=0    
614       ELSE  
615         NSIN=NSIN+1 
616       ENDIF 
617   190 CONTINUE  
618       IF(PDM.GE.PARJ(32)) GOTO 320  
619     
620 C...Fill small-mass system as cluster.  
621       NSAV=N    
622       PECM=sngl(SQRT(MAX(0D0,DPC(4)**2-DPC(1)**2-DPC(2)**2-DPC(3)**2)))
623       K(N+1,1)=11   
624       K(N+1,2)=91   
625       K(N+1,3)=IC1  
626       K(N+1,4)=N+2  
627       K(N+1,5)=N+3  
628       P(N+1,1)=sngl(DPC(1))
629       P(N+1,2)=sngl(DPC(2))  
630       P(N+1,3)=sngl(DPC(3))  
631       P(N+1,4)=sngl(DPC(4))
632       P(N+1,5)=PECM 
633     
634 C...Form two particles from flavours of lowest-mass system, if feasible.    
635       K(N+2,1)=1    
636       K(N+3,1)=1    
637       IF(MSTU(16).NE.2) THEN    
638         K(N+2,3)=N+1    
639         K(N+3,3)=N+1    
640       ELSE  
641         K(N+2,3)=IC1    
642         K(N+3,3)=IC2    
643       ENDIF 
644       K(N+2,4)=0    
645       K(N+3,4)=0    
646       K(N+2,5)=0    
647       K(N+3,5)=0    
648       IF(IABS(K(IC1,2)).NE.21) THEN 
649         KC1=LUCOMP(K(IC1,2))    
650         KC2=LUCOMP(K(IC2,2))    
651         IF(KC1.EQ.0.OR.KC2.EQ.0) GOTO 320   
652         KQ1=KCHG(KC1,2)*ISIGN(1,K(IC1,2))   
653         KQ2=KCHG(KC2,2)*ISIGN(1,K(IC2,2))   
654         IF(KQ1+KQ2.NE.0) GOTO 320   
655   200   CALL LUKFDI(K(IC1,2),0,KFLN,K(N+2,2))   
656         CALL LUKFDI(K(IC2,2),-KFLN,KFLDMP,K(N+3,2)) 
657         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 200 
658       ELSE  
659         IF(IABS(K(IC2,2)).NE.21) GOTO 320   
660   210   CALL LUKFDI(1+INT((2.+PARJ(2))*RLU(0)),0,KFLN,KFDMP)    
661         CALL LUKFDI(KFLN,0,KFLM,K(N+2,2))   
662         CALL LUKFDI(-KFLN,-KFLM,KFLDMP,K(N+3,2))    
663         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 210 
664       ENDIF 
665       P(N+2,5)=ULMASS(K(N+2,2)) 
666       P(N+3,5)=ULMASS(K(N+3,2)) 
667       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM.AND.NSIN.EQ.1) GOTO 320 
668       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) GOTO 260   
669     
670 C...Perform two-particle decay of jet system, if possible.  
671       IF(PECM.GE.0.02d0*DPC(4)) THEN  
672         PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-  
673      &  (P(N+2,5)-P(N+3,5))**2))/(2.*PECM)  
674         UE(3)=2.*RLU(0)-1.  
675         PHI=PARU(2)*RLU(0)  
676         UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)    
677         UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)    
678         DO 220 J=1,3    
679         P(N+2,J)=PA*UE(J)   
680   220   P(N+3,J)=-PA*UE(J)  
681         P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)    
682         P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)    
683         CALL LUDBRB(N+2,N+3,0.,0.,DPC(1)/DPC(4),DPC(2)/DPC(4),  
684      &  DPC(3)/DPC(4))  
685       ELSE  
686         NP=0    
687         DO 230 I=IC1,IC2    
688   230   IF(K(I,1).EQ.1.OR.K(I,1).EQ.2) NP=NP+1  
689         HA=P(IC1,4)*P(IC2,4)-P(IC1,1)*P(IC2,1)-P(IC1,2)*P(IC2,2)-   
690      &  P(IC1,3)*P(IC2,3)   
691         IF(NP.GE.3.OR.HA.LE.1.25*P(IC1,5)*P(IC2,5)) GOTO 260    
692         HD1=0.5*(P(N+2,5)**2-P(IC1,5)**2)   
693         HD2=0.5*(P(N+3,5)**2-P(IC2,5)**2)   
694         HR=SQRT(MAX(0.,((HA-HD1-HD2)**2-(P(N+2,5)*P(N+3,5))**2)/    
695      &  (HA**2-(P(IC1,5)*P(IC2,5))**2)))-1. 
696         HC=P(IC1,5)**2+2.*HA+P(IC2,5)**2    
697         HK1=((P(IC2,5)**2+HA)*HR+HD1-HD2)/HC    
698         HK2=((P(IC1,5)**2+HA)*HR+HD2-HD1)/HC    
699         DO 240 J=1,4    
700         P(N+2,J)=(1.+HK1)*P(IC1,J)-HK2*P(IC2,J) 
701   240   P(N+3,J)=(1.+HK2)*P(IC2,J)-HK1*P(IC1,J) 
702       ENDIF 
703       DO 250 J=1,4  
704       V(N+1,J)=V(IC1,J) 
705       V(N+2,J)=V(IC1,J) 
706   250 V(N+3,J)=V(IC2,J) 
707       V(N+1,5)=0.   
708       V(N+2,5)=0.   
709       V(N+3,5)=0.   
710       N=N+3 
711       GOTO 300  
712     
713 C...Else form one particle from the flavours available, if possible.    
714   260 K(N+1,5)=N+2  
715       IF(IABS(K(IC1,2)).GT.100.AND.IABS(K(IC2,2)).GT.100) THEN  
716         GOTO 320    
717       ELSEIF(IABS(K(IC1,2)).NE.21) THEN 
718         CALL LUKFDI(K(IC1,2),K(IC2,2),KFLDMP,K(N+2,2))  
719       ELSE  
720         KFLN=1+INT((2.+PARJ(2))*RLU(0)) 
721         CALL LUKFDI(KFLN,-KFLN,KFLDMP,K(N+2,2)) 
722       ENDIF 
723       IF(K(N+2,2).EQ.0) GOTO 260    
724       P(N+2,5)=ULMASS(K(N+2,2)) 
725     
726 C...Find parton/particle which combines to largest extra mass.  
727       IR=0  
728       HA=0. 
729       DO 280 MCOMB=1,3  
730       IF(IR.NE.0) GOTO 280  
731       DO 270 I=MAX(1,IP),N  
732       IF(K(I,1).LE.0.OR.K(I,1).GT.10.OR.(I.GE.IC1.AND.I.LE.IC2. 
733      &AND.K(I,1).GE.1.AND.K(I,1).LE.2)) GOTO 270    
734       IF(MCOMB.EQ.1) KCI=LUCOMP(K(I,2)) 
735       IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 270  
736       IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 270  
737       IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100) 
738      &GOTO 270  
739       HCR=sngl(DPC(4))*P(I,4)-sngl(DPC(1))*P(I,1)
740      1     -sngl(DPC(2))*P(I,2)-sngl(DPC(3))*P(I,3)   
741       IF(HCR.GT.HA) THEN    
742         IR=I    
743         HA=HCR  
744       ENDIF 
745   270 CONTINUE  
746   280 CONTINUE  
747     
748 C...Shuffle energy and momentum to put new particle on mass shell.  
749       HB=PECM**2+HA 
750       HC=P(N+2,5)**2+HA 
751       HD=P(IR,5)**2+HA
752 C******************CHANGES BY HIJING************  
753       HK2=0.0
754       IF(HA**2-(PECM*P(IR,5))**2.EQ.0.0.OR.HB+HD.EQ.0.0) GO TO 285
755 C******************
756       HK2=0.5*(HB*SQRT(((HB+HC)**2-4.*(HB+HD)*P(N+2,5)**2)/ 
757      &(HA**2-(PECM*P(IR,5))**2))-(HB+HC))/(HB+HD)   
758   285 HK1=(0.5*(P(N+2,5)**2-PECM**2)+HD*HK2)/HB 
759       DO 290 J=1,4  
760       P(N+2,J)=(1.+HK1)*sngl(DPC(J))-HK2*P(IR,J)  
761       P(IR,J)=(1.+HK2)*P(IR,J)-HK1*sngl(DPC(J))
762       V(N+1,J)=V(IC1,J) 
763   290 V(N+2,J)=V(IC1,J) 
764       V(N+1,5)=0.   
765       V(N+2,5)=0.   
766       N=N+2 
767     
768 C...Mark collapsed system and store daughter pointers. Iterate. 
769   300 DO 310 I=IC1,IC2  
770       IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.KCHG(LUCOMP(K(I,2)),2).NE.0)  
771      &THEN  
772         K(I,1)=K(I,1)+10    
773         IF(MSTU(16).NE.2) THEN  
774           K(I,4)=NSAV+1 
775           K(I,5)=NSAV+1 
776         ELSE    
777           K(I,4)=NSAV+2 
778           K(I,5)=N  
779         ENDIF   
780       ENDIF 
781   310 CONTINUE  
782       IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 140  
783     
784 C...Check flavours and invariant masses in parton systems.  
785   320 NP=0  
786       KFN=0 
787       KQS=0 
788       DO 330 J=1,5  
789   330 DPS(J)=0d0
790       DO 360 I=MAX(1,IP),N  
791       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360  
792       KC=LUCOMP(K(I,2)) 
793       IF(KC.EQ.0) GOTO 360  
794       KQ=KCHG(KC,2)*ISIGN(1,K(I,2)) 
795       IF(KQ.EQ.0) GOTO 360  
796       NP=NP+1   
797       IF(KQ.NE.2) THEN  
798         KFN=KFN+1   
799         KQS=KQS+KQ  
800         MSTJ(93)=1  
801         DPS(5)=DPS(5)+dble(ULMASS(K(I,2)))
802       ENDIF 
803       DO 340 J=1,4  
804   340 DPS(J)=DPS(J)+dble(P(I,J))
805
806 clin-4/12/01:
807 c     np: # of partons, KFN: number of quarks and diquarks, 
808 c     KC=0 for color singlet system, -1 for quarks and anti-diquarks, 
809 c     1 for quarks and anti-diquarks, and 2 for gluons:
810       IF(K(I,1).EQ.1) THEN  
811 clin-4/12/01     end of color singlet system.
812         IF(NP.NE.1.AND.(KFN.EQ.1.OR.KFN.GE.3.OR.KQS.NE.0)) CALL 
813      &  LUERRM(2,'(LUPREP:) unphysical flavour combination')    
814
815 clin-4/16/01: 'jet system' should be defined as np.ne.2:
816 c        IF(NP.NE.1.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.  
817 c     &  (0.9*PARJ(32)+DPS(5))**2) CALL LUERRM(3,    
818 c     &  '(LUPREP:) too small mass in jet system')   
819         IF(NP.NE.2.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.  
820      &  (0.9d0*dble(PARJ(32))+DPS(5))**2) then 
821            CALL LUERRM(3,    
822      &  '(LUPREP:) too small mass in jet system')   
823            write (6,*) 'DPS(1-5),KI1-5=',DPS(1),DPS(2),DPS(3),DPS(4),
824      1 DPS(5),'*',K(I,1),K(I,2),K(I,3),K(I,4),K(I,5)
825         endif
826
827         NP=0    
828         KFN=0   
829         KQS=0   
830         DO 350 J=1,5    
831   350   DPS(J)=0d0
832       ENDIF 
833   360 CONTINUE  
834     
835       RETURN    
836       END   
837     
838 C*********************************************************************  
839     
840       SUBROUTINE LUSTRF(IP) 
841 C...Purpose: to handle the fragmentation of an arbitrary colour singlet 
842 C...jet system according to the Lund string fragmentation model.    
843       IMPLICIT DOUBLE PRECISION(D)  
844       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
845       SAVE /LUJETSA/ 
846       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
847       SAVE /LUDAT1A/ 
848       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
849       SAVE /LUDAT2A/ 
850       DIMENSION DPS(5),KFL(3),PMQ(3),PX(3),PY(3),GAM(3),IE(2),PR(2),    
851      &IN(9),DHM(4),DHG(4),DP(5,5),IRANK(2),MJU(4),IJU(3),PJU(5,5),  
852      &TJU(5),KFJH(2),NJS(2),KFJS(2),PJS(4,5)    
853     
854 C...Function: four-product of two vectors.  
855       FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3) 
856       DFOUR(I,J)=DP(I,4)*DP(J,4)-DP(I,1)*DP(J,1)-DP(I,2)*DP(J,2)-   
857      &DP(I,3)*DP(J,3)   
858
859       ir=0
860       in3=0
861       jr=0
862       prev=0
863     
864 C...Reset counters. Identify parton system. 
865       MSTJ(91)=0    
866       NSAV=N    
867       NP=0  
868       KQSUM=0   
869       DO 100 J=1,5  
870   100 DPS(J)=0d0 
871       MJU(1)=0  
872       MJU(2)=0  
873       I=IP-1    
874   110 I=I+1 
875       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN 
876         CALL LUERRM(12,'(LUSTRF:) failed to reconstruct jet system')    
877         IF(MSTU(21).GE.1) RETURN    
878       ENDIF 
879       IF(K(I,1).NE.1.AND.K(I,1).NE.2.AND.K(I,1).NE.41) GOTO 110 
880       KC=LUCOMP(K(I,2)) 
881       IF(KC.EQ.0) GOTO 110  
882       KQ=KCHG(KC,2)*ISIGN(1,K(I,2)) 
883       IF(KQ.EQ.0) GOTO 110  
884       IF(N+5*NP+11.GT.MSTU(4)-MSTU(32)-5) THEN  
885         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETSA')   
886         IF(MSTU(21).GE.1) RETURN    
887       ENDIF 
888     
889 C...Take copy of partons to be considered. Check flavour sum.   
890       NP=NP+1   
891       DO 120 J=1,5  
892       K(N+NP,J)=K(I,J)  
893       P(N+NP,J)=P(I,J)  
894   120 DPS(J)=DPS(J)+dble(P(I,J))
895       K(N+NP,3)=I   
896       IF(P(N+NP,4)**2.LT.P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2) THEN   
897         P(N+NP,4)=SQRT(P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2+  
898      &  P(N+NP,5)**2)   
899         DPS(4)=DPS(4)+dble(MAX(0.,P(N+NP,4)-P(I,4)))
900       ENDIF 
901       IF(KQ.NE.2) KQSUM=KQSUM+KQ    
902       IF(K(I,1).EQ.41) THEN 
903         KQSUM=KQSUM+2*KQ    
904         IF(KQSUM.EQ.KQ) MJU(1)=N+NP 
905         IF(KQSUM.NE.KQ) MJU(2)=N+NP 
906       ENDIF 
907       IF(K(I,1).EQ.2.OR.K(I,1).EQ.41) GOTO 110  
908       IF(KQSUM.NE.0) THEN   
909         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination')  
910         IF(MSTU(21).GE.1) RETURN    
911       ENDIF 
912
913 C...Boost copied system to CM frame (for better numerical precision).   
914       CALL LUDBRB(N+1,N+NP,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4), 
915      &-DPS(3)/DPS(4))   
916
917 C...Search for very nearby partons that may be recombined.  
918       NTRYR=0   
919       PARU12=PARU(12)   
920       PARU13=PARU(13)   
921       MJU(3)=MJU(1) 
922       MJU(4)=MJU(2) 
923       NR=NP 
924   130 IF(NR.GE.3) THEN  
925         PDRMIN=2.*PARU12    
926         DO 140 I=N+1,N+NR   
927         IF(I.EQ.N+NR.AND.IABS(K(N+1,2)).NE.21) GOTO 140 
928         I1=I+1  
929         IF(I.EQ.N+NR) I1=N+1    
930         IF(K(I,1).EQ.41.OR.K(I1,1).EQ.41) GOTO 140  
931         IF(MJU(1).NE.0.AND.I1.LT.MJU(1).AND.IABS(K(I1,2)).NE.21)    
932      &  GOTO 140    
933         IF(MJU(2).NE.0.AND.I.GT.MJU(2).AND.IABS(K(I,2)).NE.21) GOTO 140 
934         PAP=SQRT((P(I,1)**2+P(I,2)**2+P(I,3)**2)*(P(I1,1)**2+   
935      &  P(I1,2)**2+P(I1,3)**2)) 
936         PVP=P(I,1)*P(I1,1)+P(I,2)*P(I1,2)+P(I,3)*P(I1,3)    
937         PDR=4.*(PAP-PVP)**2/(PARU13**2*PAP+2.*(PAP-PVP))    
938         IF(PDR.LT.PDRMIN) THEN  
939           IR=I  
940           PDRMIN=PDR    
941         ENDIF   
942   140   CONTINUE    
943     
944 C...Recombine very nearby partons to avoid machine precision problems.  
945         IF(PDRMIN.LT.PARU12.AND.IR.EQ.N+NR) THEN    
946           DO 150 J=1,4  
947   150     P(N+1,J)=P(N+1,J)+P(N+NR,J)   
948           P(N+1,5)=SQRT(MAX(0.,P(N+1,4)**2-P(N+1,1)**2-P(N+1,2)**2- 
949      &    P(N+1,3)**2)) 
950           NR=NR-1   
951           GOTO 130  
952         ELSEIF(PDRMIN.LT.PARU12) THEN   
953           DO 160 J=1,4  
954   160     P(IR,J)=P(IR,J)+P(IR+1,J) 
955           P(IR,5)=SQRT(MAX(0.,P(IR,4)**2-P(IR,1)**2-P(IR,2)**2- 
956      &    P(IR,3)**2))  
957           DO 170 I=IR+1,N+NR-1  
958           K(I,2)=K(I+1,2)   
959           DO 170 J=1,5  
960   170     P(I,J)=P(I+1,J)   
961           IF(IR.EQ.N+NR-1) K(IR,2)=K(N+NR,2)    
962           NR=NR-1   
963           IF(MJU(1).GT.IR) MJU(1)=MJU(1)-1  
964           IF(MJU(2).GT.IR) MJU(2)=MJU(2)-1  
965           GOTO 130  
966         ENDIF   
967       ENDIF 
968       NTRYR=NTRYR+1 
969     
970 C...Reset particle counter. Skip ahead if no junctions are present; 
971 C...this is usually the case!   
972       NRS=MAX(5*NR+11,NP)   
973       NTRY=0    
974   180 NTRY=NTRY+1   
975       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN   
976         PARU12=4.*PARU12    
977         PARU13=2.*PARU13    
978         GOTO 130    
979       ELSEIF(NTRY.GT.100) THEN  
980         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop') 
981         IF(MSTU(21).GE.1) RETURN    
982       ENDIF 
983       I=N+NRS   
984       IF(MJU(1).EQ.0.AND.MJU(2).EQ.0) GOTO 500  
985       DO 490 JT=1,2 
986       NJS(JT)=0 
987       IF(MJU(JT).EQ.0) GOTO 490 
988       JS=3-2*JT 
989     
990 C...Find and sum up momentum on three sides of junction. Check flavours.    
991       DO 190 IU=1,3 
992       IJU(IU)=0 
993       DO 190 J=1,5  
994   190 PJU(IU,J)=0.  
995       IU=0  
996       DO 200 I1=N+1+(JT-1)*(NR-1),N+NR+(JT-1)*(1-NR),JS 
997       IF(K(I1,2).NE.21.AND.IU.LE.2) THEN    
998         IU=IU+1 
999         IJU(IU)=I1  
1000       ENDIF 
1001       DO 200 J=1,4  
1002   200 PJU(IU,J)=PJU(IU,J)+P(I1,J)   
1003       DO 210 IU=1,3 
1004   210 PJU(IU,5)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)    
1005       IF(K(IJU(3),2)/100.NE.10*K(IJU(1),2)+K(IJU(2),2).AND. 
1006      &K(IJU(3),2)/100.NE.10*K(IJU(2),2)+K(IJU(1),2)) THEN   
1007         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination')  
1008         IF(MSTU(21).GE.1) RETURN    
1009       ENDIF 
1010     
1011 C...Calculate (approximate) boost to rest frame of junction.    
1012       T12=(PJU(1,1)*PJU(2,1)+PJU(1,2)*PJU(2,2)+PJU(1,3)*PJU(2,3))/  
1013      &(PJU(1,5)*PJU(2,5))   
1014       T13=(PJU(1,1)*PJU(3,1)+PJU(1,2)*PJU(3,2)+PJU(1,3)*PJU(3,3))/  
1015      &(PJU(1,5)*PJU(3,5))   
1016       T23=(PJU(2,1)*PJU(3,1)+PJU(2,2)*PJU(3,2)+PJU(2,3)*PJU(3,3))/  
1017      &(PJU(2,5)*PJU(3,5))   
1018       T11=SQRT((2./3.)*(1.-T12)*(1.-T13)/(1.-T23))  
1019       T22=SQRT((2./3.)*(1.-T12)*(1.-T23)/(1.-T13))  
1020       TSQ=SQRT((2.*T11*T22+T12-1.)*(1.+T12))    
1021       T1F=(TSQ-T22*(1.+T12))/(1.-T12**2)    
1022       T2F=(TSQ-T11*(1.+T12))/(1.-T12**2)    
1023       DO 220 J=1,3  
1024   220 TJU(J)=-(T1F*PJU(1,J)/PJU(1,5)+T2F*PJU(2,J)/PJU(2,5)) 
1025       TJU(4)=SQRT(1.+TJU(1)**2+TJU(2)**2+TJU(3)**2) 
1026       DO 230 IU=1,3 
1027   230 PJU(IU,5)=TJU(4)*PJU(IU,4)-TJU(1)*PJU(IU,1)-TJU(2)*PJU(IU,2)- 
1028      &TJU(3)*PJU(IU,3)  
1029     
1030 C...Put junction at rest if motion could give inconsistencies.  
1031       IF(PJU(1,5)+PJU(2,5).GT.PJU(1,4)+PJU(2,4)) THEN   
1032         DO 240 J=1,3    
1033   240   TJU(J)=0.   
1034         TJU(4)=1.   
1035         PJU(1,5)=PJU(1,4)   
1036         PJU(2,5)=PJU(2,4)   
1037         PJU(3,5)=PJU(3,4)   
1038       ENDIF 
1039     
1040 C...Start preparing for fragmentation of two strings from junction. 
1041       ISTA=I    
1042       DO 470 IU=1,2 
1043       NS=IJU(IU+1)-IJU(IU)  
1044     
1045 C...Junction strings: find longitudinal string directions.  
1046       DO 260 IS=1,NS    
1047       IS1=IJU(IU)+IS-1  
1048       IS2=IJU(IU)+IS    
1049       DO 250 J=1,5  
1050       DP(1,J)=dble(0.5*P(IS1,J))
1051       IF(IS.EQ.1) DP(1,J)=dble(P(IS1,J))
1052       DP(2,J)=dble(0.5*P(IS2,J))
1053   250 IF(IS.EQ.NS) DP(2,J)=-dble(PJU(IU,J))
1054       IF(IS.EQ.NS) DP(2,4)=dble(
1055      1     SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2))
1056       IF(IS.EQ.NS) DP(2,5)=0d0   
1057       DP(3,5)=DFOUR(1,1)    
1058       DP(4,5)=DFOUR(2,2)    
1059       DHKC=DFOUR(1,2)   
1060       IF(DP(3,5)+2d0*DHKC+DP(4,5).LE.0d0) THEN    
1061         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)  
1062         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)  
1063         DP(3,5)=0D0 
1064         DP(4,5)=0D0 
1065         DHKC=DFOUR(1,2) 
1066       ENDIF 
1067       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))    
1068       DHK1=0.5d0*((DP(4,5)+DHKC)/DHKS-1d0) 
1069       DHK2=0.5d0*((DP(3,5)+DHKC)/DHKS-1d0) 
1070       IN1=N+NR+4*IS-3   
1071       P(IN1,5)=sngl(SQRT(DP(3,5)+2d0*DHKC+DP(4,5)))
1072       DO 260 J=1,4  
1073       P(IN1,J)=sngl((1d0+DHK1)*DP(1,J)-DHK2*DP(2,J))
1074   260 P(IN1+1,J)=sngl((1d0+DHK2)*DP(2,J)-DHK1*DP(1,J))
1075     
1076 C...Junction strings: initialize flavour, momentum and starting pos.    
1077       ISAV=I    
1078   270 NTRY=NTRY+1   
1079       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN   
1080         PARU12=4.*PARU12    
1081         PARU13=2.*PARU13    
1082         GOTO 130    
1083       ELSEIF(NTRY.GT.100) THEN  
1084         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop') 
1085         IF(MSTU(21).GE.1) RETURN    
1086       ENDIF 
1087       I=ISAV    
1088       IRANKJ=0  
1089       IE(1)=K(N+1+(JT/2)*(NP-1),3)  
1090       IN(4)=N+NR+1  
1091       IN(5)=IN(4)+1 
1092       IN(6)=N+NR+4*NS+1 
1093       DO 280 JQ=1,2 
1094       DO 280 IN1=N+NR+2+JQ,N+NR+4*NS-2+JQ,4 
1095       P(IN1,1)=2-JQ 
1096       P(IN1,2)=JQ-1 
1097   280 P(IN1,3)=1.   
1098       KFL(1)=K(IJU(IU),2)   
1099       PX(1)=0.  
1100       PY(1)=0.  
1101       GAM(1)=0. 
1102       DO 290 J=1,5  
1103   290 PJU(IU+3,J)=0.    
1104     
1105 C...Junction strings: find initial transverse directions.   
1106       DO 300 J=1,4  
1107       DP(1,J)=dble(P(IN(4),J))
1108       DP(2,J)=dble(P(IN(4)+1,J))
1109       DP(3,J)=0d0    
1110   300 DP(4,J)=0d0    
1111       DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)    
1112       DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)    
1113       DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)   
1114       DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)   
1115       DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)   
1116       IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1d0    
1117       IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1d0    
1118       IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1d0    
1119       IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1d0    
1120       DHC12=DFOUR(1,2)  
1121       DHCX1=DFOUR(3,1)/DHC12    
1122       DHCX2=DFOUR(3,2)/DHC12    
1123       DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12) 
1124       DHCY1=DFOUR(4,1)/DHC12    
1125       DHCY2=DFOUR(4,2)/DHC12    
1126       DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12   
1127       DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)    
1128       DO 310 J=1,4  
1129       DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))   
1130       P(IN(6),J)=sngl(DP(3,J))
1131   310 P(IN(6)+1,J)=sngl(DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-  
1132      &DHCYX*DP(3,J)))    
1133     
1134 C...Junction strings: produce new particle, origin. 
1135   320 I=I+1 
1136       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN   
1137         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETSA')   
1138         IF(MSTU(21).GE.1) RETURN    
1139       ENDIF 
1140       IRANKJ=IRANKJ+1   
1141       K(I,1)=1  
1142       K(I,3)=IE(1)  
1143       K(I,4)=0  
1144       K(I,5)=0  
1145     
1146 C...Junction strings: generate flavour, hadron, pT, z and Gamma.    
1147   330 CALL LUKFDI(KFL(1),0,KFL(3),K(I,2))   
1148       IF(K(I,2).EQ.0) GOTO 270  
1149       IF(MSTJ(12).GE.3.AND.IRANKJ.EQ.1.AND.IABS(KFL(1)).LE.10.AND.  
1150      &IABS(KFL(3)).GT.10) THEN  
1151         IF(RLU(0).GT.PARJ(19)) GOTO 330 
1152       ENDIF 
1153       P(I,5)=ULMASS(K(I,2)) 
1154       CALL LUPTDI(KFL(1),PX(3),PY(3))   
1155       PR(1)=P(I,5)**2+(PX(1)+PX(3))**2+(PY(1)+PY(3))**2 
1156       CALL LUZDIS(KFL(1),KFL(3),PR(1),Z)    
1157       GAM(3)=(1.-Z)*(GAM(1)+PR(1)/Z)    
1158       DO 340 J=1,3  
1159   340 IN(J)=IN(3+J) 
1160
1161 C...Junction strings: stepping within or from 'low' string region easy. 
1162       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*  
1163      &P(IN(1),5)**2.GE.PR(1)) THEN  
1164         P(IN(1)+2,4)=Z*P(IN(1)+2,3) 
1165         P(IN(2)+2,4)=PR(1)/(P(IN(1)+2,4)*P(IN(1),5)**2) 
1166         DO 350 J=1,4    
1167   350   P(I,J)=(PX(1)+PX(3))*P(IN(3),J)+(PY(1)+PY(3))*P(IN(3)+1,J)  
1168         GOTO 420    
1169       ELSEIF(IN(1)+1.EQ.IN(2)) THEN 
1170         P(IN(2)+2,4)=P(IN(2)+2,3)   
1171         P(IN(2)+2,1)=1. 
1172         IN(2)=IN(2)+4   
1173         IF(IN(2).GT.N+NR+4*NS) GOTO 270 
1174         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN  
1175           P(IN(1)+2,4)=P(IN(1)+2,3) 
1176           P(IN(1)+2,1)=0.   
1177           IN(1)=IN(1)+4 
1178         ENDIF   
1179       ENDIF 
1180     
1181 C...Junction strings: find new transverse directions.   
1182   360 IF(IN(1).GT.N+NR+4*NS.OR.IN(2).GT.N+NR+4*NS.OR.   
1183      &IN(1).GT.IN(2)) GOTO 270  
1184       IF(IN(1).NE.IN(4).OR.IN(2).NE.IN(5)) THEN 
1185         DO 370 J=1,4    
1186         DP(1,J)=dble(P(IN(1),J))
1187         DP(2,J)=dble(P(IN(2),J))
1188         DP(3,J)=0d0  
1189   370   DP(4,J)=0d0  
1190         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)  
1191         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)  
1192         DHC12=DFOUR(1,2)    
1193         IF(DHC12.LE.1E-2) THEN  
1194           P(IN(1)+2,4)=P(IN(1)+2,3) 
1195           P(IN(1)+2,1)=0.   
1196           IN(1)=IN(1)+4 
1197           GOTO 360  
1198         ENDIF   
1199         IN(3)=N+NR+4*NS+5   
1200         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4) 
1201         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4) 
1202         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4) 
1203         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1d0  
1204         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1d0  
1205         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1d0  
1206         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1d0  
1207         DHCX1=DFOUR(3,1)/DHC12  
1208         DHCX2=DFOUR(3,2)/DHC12  
1209         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)   
1210         DHCY1=DFOUR(4,1)/DHC12  
1211         DHCY2=DFOUR(4,2)/DHC12  
1212         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12 
1213         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)  
1214         DO 380 J=1,4    
1215         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J)) 
1216         P(IN(3),J)=sngl(DP(3,J))
1217   380   P(IN(3)+1,J)=sngl(DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-    
1218      &  DHCYX*DP(3,J)))  
1219 C...Express pT with respect to new axes, if sensible.   
1220         PXP=-(PX(3)*FOUR(IN(6),IN(3))+PY(3)*FOUR(IN(6)+1,IN(3)))    
1221         PYP=-(PX(3)*FOUR(IN(6),IN(3)+1)+PY(3)*FOUR(IN(6)+1,IN(3)+1))    
1222         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN   
1223           PX(3)=PXP 
1224           PY(3)=PYP 
1225         ENDIF   
1226       ENDIF 
1227     
1228 C...Junction strings: sum up known four-momentum, coefficients for m2.  
1229       DO 400 J=1,4  
1230       DHG(J)=0d0 
1231       P(I,J)=PX(1)*P(IN(6),J)+PY(1)*P(IN(6)+1,J)+PX(3)*P(IN(3),J)+  
1232      &PY(3)*P(IN(3)+1,J)    
1233       DO 390 IN1=IN(4),IN(1)-4,4    
1234   390 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J) 
1235       DO 400 IN2=IN(5),IN(2)-4,4    
1236   400 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J) 
1237       DHM(1)=dble(FOUR(I,I))
1238       DHM(2)=dble(2.*FOUR(I,IN(1)))   
1239       DHM(3)=dble(2.*FOUR(I,IN(2)))  
1240       DHM(4)=dble(2.*FOUR(IN(1),IN(2))) 
1241     
1242 C...Junction strings: find coefficients for Gamma expression.   
1243       DO 410 IN2=IN(1)+1,IN(2),4    
1244       DO 410 IN1=IN(1),IN2-1,4  
1245       DHC=dble(2.*FOUR(IN1,IN2))
1246       DHG(1)=DHG(1)+dble(P(IN1+2,1)*P(IN2+2,1))*DHC   
1247       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-dble(P(IN2+2,1))*DHC 
1248       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+dble(P(IN1+2,1))*DHC 
1249   410 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC   
1250     
1251 C...Junction strings: solve (m2, Gamma) equation system for energies.   
1252       DHS1=DHM(3)*DHG(4)-DHM(4)*DHG(3)  
1253       IF(ABS(DHS1).LT.1E-4) GOTO 270    
1254       DHS2=DHM(4)*(dble(GAM(3))-DHG(1))-DHM(2)*DHG(3)-DHG(4)* 
1255      &(dble(P(I,5))**2-DHM(1))+DHG(2)*DHM(3)  
1256       DHS3=DHM(2)*(dble(GAM(3))-DHG(1))
1257      1     -DHG(2)*(dble(P(I,5))**2-DHM(1)) 
1258       P(IN(2)+2,4)=0.5*sngl(SQRT(MAX(0D0,DHS2**2-4d0*DHS1*DHS3))
1259      &     /ABS(DHS1)-DHS2/DHS1)
1260       IF(DHM(2)+DHM(4)*dble(P(IN(2)+2,4)).LE.0d0) GOTO 270 
1261       P(IN(1)+2,4)=(P(I,5)**2-sngl(DHM(1))-sngl(DHM(3))*P(IN(2)+2,4))/  
1262      &(sngl(DHM(2))+sngl(DHM(4))*P(IN(2)+2,4))  
1263
1264 C...Junction strings: step to new region if necessary.  
1265       IF(P(IN(2)+2,4).GT.P(IN(2)+2,3)) THEN 
1266         P(IN(2)+2,4)=P(IN(2)+2,3)   
1267         P(IN(2)+2,1)=1. 
1268         IN(2)=IN(2)+4   
1269         IF(IN(2).GT.N+NR+4*NS) GOTO 270 
1270         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN  
1271           P(IN(1)+2,4)=P(IN(1)+2,3) 
1272           P(IN(1)+2,1)=0.   
1273           IN(1)=IN(1)+4 
1274         ENDIF   
1275         GOTO 360    
1276       ELSEIF(P(IN(1)+2,4).GT.P(IN(1)+2,3)) THEN 
1277         P(IN(1)+2,4)=P(IN(1)+2,3)   
1278         P(IN(1)+2,1)=0. 
1279         IN(1)=IN(1)+JS  
1280         GOTO 710    
1281       ENDIF 
1282     
1283 C...Junction strings: particle four-momentum, remainder, loop back. 
1284   420 DO 430 J=1,4  
1285       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J) 
1286   430 PJU(IU+3,J)=PJU(IU+3,J)+P(I,J)    
1287       IF(P(I,4).LE.0.) GOTO 270 
1288       PJU(IU+3,5)=TJU(4)*PJU(IU+3,4)-TJU(1)*PJU(IU+3,1)-    
1289      &TJU(2)*PJU(IU+3,2)-TJU(3)*PJU(IU+3,3) 
1290       IF(PJU(IU+3,5).LT.PJU(IU,5)) THEN 
1291         KFL(1)=-KFL(3)  
1292         PX(1)=-PX(3)    
1293         PY(1)=-PY(3)    
1294         GAM(1)=GAM(3)   
1295         IF(IN(3).NE.IN(6)) THEN 
1296           DO 440 J=1,4  
1297           P(IN(6),J)=P(IN(3),J) 
1298   440     P(IN(6)+1,J)=P(IN(3)+1,J) 
1299         ENDIF   
1300         DO 450 JQ=1,2   
1301         IN(3+JQ)=IN(JQ) 
1302         P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)   
1303   450   P(IN(JQ)+2,1)=P(IN(JQ)+2,1)-(3-2*JQ)*P(IN(JQ)+2,4)  
1304         GOTO 320    
1305       ENDIF 
1306     
1307 C...Junction strings: save quantities left after each string.   
1308       IF(IABS(KFL(1)).GT.10) GOTO 270   
1309       I=I-1 
1310       KFJH(IU)=KFL(1)   
1311       DO 460 J=1,4  
1312   460 PJU(IU+3,J)=PJU(IU+3,J)-P(I+1,J)  
1313   470 CONTINUE  
1314     
1315 C...Junction strings: put together to new effective string endpoint.    
1316       NJS(JT)=I-ISTA    
1317       KFJS(JT)=K(K(MJU(JT+2),3),2)  
1318       KFLS=2*INT(RLU(0)+3.*PARJ(4)/(1.+3.*PARJ(4)))+1   
1319       IF(KFJH(1).EQ.KFJH(2)) KFLS=3 
1320       IF(ISTA.NE.I) KFJS(JT)=ISIGN(1000*MAX(IABS(KFJH(1)),  
1321      &IABS(KFJH(2)))+100*MIN(IABS(KFJH(1)),IABS(KFJH(2)))+  
1322      &KFLS,KFJH(1)) 
1323       DO 480 J=1,4  
1324       PJS(JT,J)=PJU(1,J)+PJU(2,J)+P(MJU(JT),J)  
1325   480 PJS(JT+2,J)=PJU(4,J)+PJU(5,J) 
1326       PJS(JT,5)=SQRT(MAX(0.,PJS(JT,4)**2-PJS(JT,1)**2-PJS(JT,2)**2- 
1327      &PJS(JT,3)**2))    
1328   490 CONTINUE  
1329     
1330 C...Open versus closed strings. Choose breakup region for latter.   
1331   500 IF(MJU(1).NE.0.AND.MJU(2).NE.0) THEN  
1332         NS=MJU(2)-MJU(1)    
1333         NB=MJU(1)-N 
1334       ELSEIF(MJU(1).NE.0) THEN  
1335         NS=N+NR-MJU(1)  
1336         NB=MJU(1)-N 
1337       ELSEIF(MJU(2).NE.0) THEN  
1338         NS=MJU(2)-N 
1339         NB=1    
1340       ELSEIF(IABS(K(N+1,2)).NE.21) THEN 
1341         NS=NR-1 
1342         NB=1    
1343       ELSE  
1344         NS=NR+1 
1345         W2SUM=0.    
1346         DO 510 IS=1,NR  
1347         P(N+NR+IS,1)=0.5*FOUR(N+IS,N+IS+1-NR*(IS/NR))   
1348   510   W2SUM=W2SUM+P(N+NR+IS,1)    
1349         W2RAN=RLU(0)*W2SUM  
1350         NB=0    
1351   520   NB=NB+1 
1352         W2SUM=W2SUM-P(N+NR+NB,1)    
1353         IF(W2SUM.GT.W2RAN.AND.NB.LT.NR) GOTO 520    
1354       ENDIF 
1355     
1356 C...Find longitudinal string directions (i.e. lightlike four-vectors).  
1357       DO 540 IS=1,NS    
1358       IS1=N+IS+NB-1-NR*((IS+NB-2)/NR)   
1359       IS2=N+IS+NB-NR*((IS+NB-1)/NR) 
1360       DO 530 J=1,5  
1361       DP(1,J)=dble(P(IS1,J))
1362       IF(IABS(K(IS1,2)).EQ.21) DP(1,J)=0.5d0*DP(1,J)  
1363       IF(IS1.EQ.MJU(1)) DP(1,J)=dble(PJS(1,J)-PJS(3,J))
1364       DP(2,J)=dble(P(IS2,J))
1365       IF(IABS(K(IS2,2)).EQ.21) DP(2,J)=0.5d0*DP(2,J)  
1366   530 IF(IS2.EQ.MJU(2)) DP(2,J)=dble(PJS(2,J)-PJS(4,J))
1367       DP(3,5)=DFOUR(1,1)    
1368       DP(4,5)=DFOUR(2,2)    
1369       DHKC=DFOUR(1,2)   
1370       IF(DP(3,5)+2.d0*DHKC+DP(4,5).LE.0.d0) THEN    
1371         DP(3,5)=DP(1,5)**2  
1372         DP(4,5)=DP(2,5)**2  
1373         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2+DP(1,5)**2)   
1374         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2+DP(2,5)**2)   
1375         DHKC=DFOUR(1,2) 
1376       ENDIF 
1377       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))    
1378       DHK1=0.5d0*((DP(4,5)+DHKC)/DHKS-1.d0) 
1379       DHK2=0.5d0*((DP(3,5)+DHKC)/DHKS-1.d0) 
1380       IN1=N+NR+4*IS-3   
1381       P(IN1,5)=SQRT(sngl(DP(3,5)+2.d0*DHKC+DP(4,5)))
1382       DO 540 J=1,4  
1383       P(IN1,J)=sngl((1.d0+DHK1)*DP(1,J)-DHK2*DP(2,J))
1384   540 P(IN1+1,J)=sngl((1.d0+DHK2)*DP(2,J)-DHK1*DP(1,J))
1385     
1386 C...Begin initialization: sum up energy, set starting position. 
1387       ISAV=I    
1388   550 NTRY=NTRY+1   
1389       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN   
1390         PARU12=4.*PARU12    
1391         PARU13=2.*PARU13    
1392         GOTO 130    
1393       ELSEIF(NTRY.GT.100) THEN  
1394         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop') 
1395         IF(MSTU(21).GE.1) RETURN    
1396       ENDIF 
1397       I=ISAV    
1398       DO 560 J=1,4  
1399       P(N+NRS,J)=0. 
1400       DO 560 IS=1,NR    
1401   560 P(N+NRS,J)=P(N+NRS,J)+P(N+IS,J)   
1402       DO 570 JT=1,2 
1403       IRANK(JT)=0   
1404       IF(MJU(JT).NE.0) IRANK(JT)=NJS(JT)    
1405       IF(NS.GT.NR) IRANK(JT)=1  
1406       IE(JT)=K(N+1+(JT/2)*(NP-1),3) 
1407       IN(3*JT+1)=N+NR+1+4*(JT/2)*(NS-1) 
1408       IN(3*JT+2)=IN(3*JT+1)+1   
1409       IN(3*JT+3)=N+NR+4*NS+2*JT-1   
1410       DO 570 IN1=N+NR+2+JT,N+NR+4*NS-2+JT,4 
1411       P(IN1,1)=2-JT 
1412       P(IN1,2)=JT-1 
1413   570 P(IN1,3)=1.   
1414     
1415 C...Initialize flavour and pT variables for open string.    
1416       IF(NS.LT.NR) THEN 
1417         PX(1)=0.    
1418         PY(1)=0.    
1419         IF(NS.EQ.1.AND.MJU(1)+MJU(2).EQ.0) CALL LUPTDI(0,PX(1),PY(1))   
1420         PX(2)=-PX(1)    
1421         PY(2)=-PY(1)    
1422         DO 580 JT=1,2   
1423         KFL(JT)=K(IE(JT),2) 
1424         IF(MJU(JT).NE.0) KFL(JT)=KFJS(JT)   
1425         MSTJ(93)=1  
1426         PMQ(JT)=ULMASS(KFL(JT)) 
1427   580   GAM(JT)=0.  
1428     
1429 C...Closed string: random initial breakup flavour, pT and vertex.   
1430       ELSE  
1431         KFL(3)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)    
1432         CALL LUKFDI(KFL(3),0,KFL(1),KDUMP)  
1433         KFL(2)=-KFL(1)  
1434         IF(IABS(KFL(1)).GT.10.AND.RLU(0).GT.0.5) THEN   
1435           KFL(2)=-(KFL(1)+ISIGN(10000,KFL(1)))  
1436         ELSEIF(IABS(KFL(1)).GT.10) THEN 
1437           KFL(1)=-(KFL(2)+ISIGN(10000,KFL(2)))  
1438         ENDIF   
1439         CALL LUPTDI(KFL(1),PX(1),PY(1)) 
1440         PX(2)=-PX(1)    
1441         PY(2)=-PY(1)    
1442         PR3=MIN(25.,0.1*P(N+NR+1,5)**2) 
1443   590   CALL LUZDIS(KFL(1),KFL(2),PR3,Z)    
1444         ZR=PR3/(Z*P(N+NR+1,5)**2)   
1445         IF(ZR.GE.1.) GOTO 590   
1446
1447         DO 600 JT=1,2   
1448         MSTJ(93)=1  
1449         PMQ(JT)=ULMASS(KFL(JT)) 
1450         GAM(JT)=PR3*(1.-Z)/Z    
1451         IN1=N+NR+3+4*(JT/2)*(NS-1)  
1452         P(IN1,JT)=1.-Z  
1453         P(IN1,3-JT)=JT-1    
1454         P(IN1,3)=(2-JT)*(1.-Z)+(JT-1)*Z 
1455         P(IN1+1,JT)=ZR  
1456         P(IN1+1,3-JT)=2-JT  
1457   600   P(IN1+1,3)=(2-JT)*(1.-ZR)+(JT-1)*ZR 
1458       ENDIF 
1459     
1460 C...Find initial transverse directions (i.e. spacelike four-vectors).   
1461       DO 640 JT=1,2 
1462       IF(JT.EQ.1.OR.NS.EQ.NR-1) THEN    
1463         IN1=IN(3*JT+1)  
1464         IN3=IN(3*JT+3)  
1465         DO 610 J=1,4    
1466         DP(1,J)=dble(P(IN1,J))
1467         DP(2,J)=dble(P(IN1+1,J))
1468         DP(3,J)=0.d0
1469   610   DP(4,J)=0.d0
1470         DP(1,4)=DSQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)  
1471         DP(2,4)=DSQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)  
1472         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4) 
1473         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4) 
1474         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4) 
1475         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.d0
1476         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.d0
1477         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.d0
1478         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.d0
1479         DHC12=DFOUR(1,2)    
1480         DHCX1=DFOUR(3,1)/DHC12  
1481         DHCX2=DFOUR(3,2)/DHC12  
1482         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)   
1483         DHCY1=DFOUR(4,1)/DHC12  
1484         DHCY2=DFOUR(4,2)/DHC12  
1485         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12 
1486         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)  
1487         DO 620 J=1,4    
1488         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J)) 
1489         P(IN3,J)=sngl(DP(3,J))
1490   620   P(IN3+1,J)=sngl(DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-  
1491      &  DHCYX*DP(3,J)))
1492       ELSE  
1493         DO 630 J=1,4    
1494         P(IN3+2,J)=P(IN3,J) 
1495   630   P(IN3+3,J)=P(IN3+1,J)   
1496       ENDIF 
1497   640 CONTINUE  
1498     
1499 C...Remove energy used up in junction string fragmentation. 
1500       IF(MJU(1)+MJU(2).GT.0) THEN   
1501         DO 660 JT=1,2   
1502         IF(NJS(JT).EQ.0) GOTO 660   
1503         DO 650 J=1,4    
1504   650   P(N+NRS,J)=P(N+NRS,J)-PJS(JT+2,J)   
1505   660   CONTINUE    
1506       ENDIF 
1507     
1508 C...Produce new particle: side, origin. 
1509   670 I=I+1 
1510       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN   
1511         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETSA')   
1512         IF(MSTU(21).GE.1) RETURN    
1513       ENDIF 
1514       JT=int(1.5+RLU(0))
1515       IF(IABS(KFL(3-JT)).GT.10) JT=3-JT 
1516       JR=3-JT   
1517       JS=3-2*JT 
1518       IRANK(JT)=IRANK(JT)+1 
1519       K(I,1)=1  
1520       K(I,3)=IE(JT) 
1521       K(I,4)=0  
1522       K(I,5)=0  
1523     
1524 C...Generate flavour, hadron and pT.    
1525   680 CALL LUKFDI(KFL(JT),0,KFL(3),K(I,2))  
1526       IF(K(I,2).EQ.0) GOTO 550  
1527       IF(MSTJ(12).GE.3.AND.IRANK(JT).EQ.1.AND.IABS(KFL(JT)).LE.10.AND.  
1528      &IABS(KFL(3)).GT.10) THEN  
1529         IF(RLU(0).GT.PARJ(19)) GOTO 680 
1530       ENDIF 
1531       P(I,5)=ULMASS(K(I,2)) 
1532       CALL LUPTDI(KFL(JT),PX(3),PY(3))  
1533       PR(JT)=P(I,5)**2+(PX(JT)+PX(3))**2+(PY(JT)+PY(3))**2  
1534     
1535 C...Final hadrons for small invariant mass. 
1536       MSTJ(93)=1    
1537       PMQ(3)=ULMASS(KFL(3)) 
1538       WMIN=PARJ(32+MSTJ(11))+PMQ(1)+PMQ(2)+PARJ(36)*PMQ(3)  
1539       IF(IABS(KFL(JT)).GT.10.AND.IABS(KFL(3)).GT.10) WMIN=  
1540      &WMIN-0.5*PARJ(36)*PMQ(3)  
1541       WREM2=FOUR(N+NRS,N+NRS)   
1542       IF(WREM2.LT.0.10) GOTO 550    
1543       IF(WREM2.LT.MAX(WMIN*(1.+(2.*RLU(0)-1.)*PARJ(37)),    
1544      &PARJ(32)+PMQ(1)+PMQ(2))**2) GOTO 810  
1545     
1546 C...Choose z, which gives Gamma. Shift z for heavy flavours.    
1547       CALL LUZDIS(KFL(JT),KFL(3),PR(JT),Z)  
1548
1549       KFL1A=IABS(KFL(1))    
1550       KFL2A=IABS(KFL(2))    
1551       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),    
1552      &MOD(KFL2A/1000,10)).GE.4) THEN    
1553         PR(JR)=(PMQ(JR)+PMQ(3))**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2  
1554         PW12=SQRT(MAX(0.,(WREM2-PR(1)-PR(2))**2-4.*PR(1)*PR(2)))    
1555         Z=(WREM2+PR(JT)-PR(JR)+PW12*(2.*Z-1.))/(2.*WREM2)   
1556         PR(JR)=(PMQ(JR)+PARJ(32+MSTJ(11)))**2+(PX(JR)-PX(3))**2+    
1557      &  (PY(JR)-PY(3))**2   
1558         IF((1.-Z)*(WREM2-PR(JT)/Z).LT.PR(JR)) GOTO 810  
1559       ENDIF 
1560       GAM(3)=(1.-Z)*(GAM(JT)+PR(JT)/Z)  
1561       DO 690 J=1,3  
1562   690 IN(J)=IN(3*JT+J)  
1563     
1564 C...Stepping within or from 'low' string region easy.   
1565       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*  
1566      &P(IN(1),5)**2.GE.PR(JT)) THEN 
1567         P(IN(JT)+2,4)=Z*P(IN(JT)+2,3)   
1568         P(IN(JR)+2,4)=PR(JT)/(P(IN(JT)+2,4)*P(IN(1),5)**2)  
1569         DO 700 J=1,4    
1570   700   P(I,J)=(PX(JT)+PX(3))*P(IN(3),J)+(PY(JT)+PY(3))*P(IN(3)+1,J)    
1571         GOTO 770    
1572       ELSEIF(IN(1)+1.EQ.IN(2)) THEN 
1573         P(IN(JR)+2,4)=P(IN(JR)+2,3) 
1574         P(IN(JR)+2,JT)=1.   
1575         IN(JR)=IN(JR)+4*JS  
1576         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 550   
1577         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN  
1578           P(IN(JT)+2,4)=P(IN(JT)+2,3)   
1579           P(IN(JT)+2,JT)=0. 
1580           IN(JT)=IN(JT)+4*JS    
1581         ENDIF   
1582       ENDIF 
1583     
1584 C...Find new transverse directions (i.e. spacelike string vectors). 
1585   710 IF(JS*IN(1).GT.JS*IN(3*JR+1).OR.JS*IN(2).GT.JS*IN(3*JR+2).OR. 
1586      &IN(1).GT.IN(2)) GOTO 550  
1587       IF(IN(1).NE.IN(3*JT+1).OR.IN(2).NE.IN(3*JT+2)) THEN   
1588         DO 720 J=1,4    
1589         DP(1,J)=dble(P(IN(1),J))
1590         DP(2,J)=dble(P(IN(2),J))
1591         DP(3,J)=0.d0
1592   720   DP(4,J)=0.d0
1593         DP(1,4)=DSQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)  
1594         DP(2,4)=DSQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)  
1595         DHC12=DFOUR(1,2)    
1596         IF(DHC12.LE.1E-2) THEN  
1597           P(IN(JT)+2,4)=P(IN(JT)+2,3)   
1598           P(IN(JT)+2,JT)=0. 
1599           IN(JT)=IN(JT)+4*JS    
1600           GOTO 710  
1601         ENDIF   
1602         IN(3)=N+NR+4*NS+5   
1603         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4) 
1604         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4) 
1605         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4) 
1606         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.d0
1607         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.d0
1608         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.d0
1609         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.d0
1610         DHCX1=DFOUR(3,1)/DHC12  
1611         DHCX2=DFOUR(3,2)/DHC12  
1612         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)   
1613         DHCY1=DFOUR(4,1)/DHC12  
1614         DHCY2=DFOUR(4,2)/DHC12  
1615         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12 
1616         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)  
1617         DO 730 J=1,4    
1618         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J)) 
1619         P(IN(3),J)=sngl(DP(3,J))
1620   730   P(IN(3)+1,J)=sngl(DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-    
1621      &  DHCYX*DP(3,J))) 
1622 C...Express pT with respect to new axes, if sensible.   
1623         PXP=-(PX(3)*FOUR(IN(3*JT+3),IN(3))+PY(3)*   
1624      &  FOUR(IN(3*JT+3)+1,IN(3)))   
1625         PYP=-(PX(3)*FOUR(IN(3*JT+3),IN(3)+1)+PY(3)* 
1626      &  FOUR(IN(3*JT+3)+1,IN(3)+1)) 
1627         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN   
1628           PX(3)=PXP 
1629           PY(3)=PYP 
1630         ENDIF   
1631       ENDIF 
1632     
1633 C...Sum up known four-momentum. Gives coefficients for m2 expression.   
1634       DO 750 J=1,4  
1635       DHG(J)=0.d0
1636       P(I,J)=PX(JT)*P(IN(3*JT+3),J)+PY(JT)*P(IN(3*JT+3)+1,J)+   
1637      &PX(3)*P(IN(3),J)+PY(3)*P(IN(3)+1,J)   
1638       DO 740 IN1=IN(3*JT+1),IN(1)-4*JS,4*JS 
1639   740 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J) 
1640       DO 750 IN2=IN(3*JT+2),IN(2)-4*JS,4*JS 
1641   750 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J) 
1642       DHM(1)=dble(FOUR(I,I))
1643       DHM(2)=dble(2.*FOUR(I,IN(1)))  
1644       DHM(3)=dble(2.*FOUR(I,IN(2)))
1645       DHM(4)=dble(2.*FOUR(IN(1),IN(2)))
1646     
1647 C...Find coefficients for Gamma expression. 
1648       DO 760 IN2=IN(1)+1,IN(2),4    
1649       DO 760 IN1=IN(1),IN2-1,4  
1650       DHC=dble(2.*FOUR(IN1,IN2))
1651       DHG(1)=DHG(1)+dble(P(IN1+2,JT)*P(IN2+2,JT))*DHC 
1652       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-dble(float(JS)*P(IN2+2,JT))*DHC 
1653       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+dble(float(JS)*P(IN1+2,JT))*DHC 
1654   760 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC   
1655     
1656 C...Solve (m2, Gamma) equation system for energies taken.   
1657       DHS1=DHM(JR+1)*DHG(4)-DHM(4)*DHG(JR+1)    
1658       IF(ABS(DHS1).LT.1E-4) GOTO 550    
1659       DHS2=DHM(4)*(dble(GAM(3))-DHG(1))-DHM(JT+1)*DHG(JR+1)-DHG(4)*   
1660      &(dble(P(I,5))**2-DHM(1))+DHG(JT+1)*DHM(JR+1)    
1661       DHS3=DHM(JT+1)*(dble(GAM(3))-DHG(1))-DHG(JT+1)
1662      &     *(dble(P(I,5))**2-DHM(1))   
1663       P(IN(JR)+2,4)=0.5*sngl((SQRT(MAX(0D0,DHS2**2-4.d0*DHS1*DHS3)))
1664      &/ABS(DHS1)-DHS2/DHS1)
1665       IF(DHM(JT+1)+DHM(4)*dble(P(IN(JR)+2,4)).LE.0.d0) GOTO 550 
1666       P(IN(JT)+2,4)=(P(I,5)**2-sngl(DHM(1))-sngl(DHM(JR+1))
1667      &     *P(IN(JR)+2,4))/(sngl(DHM(JT+1))+sngl(DHM(4))*P(IN(JR)+2,4))
1668     
1669 C...Step to new region if necessary.    
1670       IF(P(IN(JR)+2,4).GT.P(IN(JR)+2,3)) THEN   
1671         P(IN(JR)+2,4)=P(IN(JR)+2,3) 
1672         P(IN(JR)+2,JT)=1.   
1673         IN(JR)=IN(JR)+4*JS  
1674         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 550   
1675         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN  
1676           P(IN(JT)+2,4)=P(IN(JT)+2,3)   
1677           P(IN(JT)+2,JT)=0. 
1678           IN(JT)=IN(JT)+4*JS    
1679         ENDIF   
1680         GOTO 710    
1681       ELSEIF(P(IN(JT)+2,4).GT.P(IN(JT)+2,3)) THEN   
1682         P(IN(JT)+2,4)=P(IN(JT)+2,3) 
1683         P(IN(JT)+2,JT)=0.   
1684         IN(JT)=IN(JT)+4*JS  
1685         GOTO 710    
1686       ENDIF 
1687     
1688 C...Four-momentum of particle. Remaining quantities. Loop back. 
1689   770 DO 780 J=1,4  
1690       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J) 
1691   780 P(N+NRS,J)=P(N+NRS,J)-P(I,J)  
1692       IF(P(I,4).LE.0.) GOTO 550 
1693       KFL(JT)=-KFL(3)   
1694       PMQ(JT)=PMQ(3)    
1695       PX(JT)=-PX(3) 
1696       PY(JT)=-PY(3) 
1697       GAM(JT)=GAM(3)    
1698       IF(IN(3).NE.IN(3*JT+3)) THEN  
1699         DO 790 J=1,4    
1700         P(IN(3*JT+3),J)=P(IN(3),J)  
1701   790   P(IN(3*JT+3)+1,J)=P(IN(3)+1,J)  
1702       ENDIF 
1703       DO 800 JQ=1,2 
1704       IN(3*JT+JQ)=IN(JQ)    
1705       P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4) 
1706   800 P(IN(JQ)+2,JT)=P(IN(JQ)+2,JT)-JS*(3-2*JQ)*P(IN(JQ)+2,4)   
1707       GOTO 670  
1708     
1709 C...Final hadron: side, flavour, hadron, mass.  
1710   810 I=I+1 
1711       K(I,1)=1  
1712       K(I,3)=IE(JR) 
1713       K(I,4)=0  
1714       K(I,5)=0  
1715       CALL LUKFDI(KFL(JR),-KFL(3),KFLDMP,K(I,2))    
1716       IF(K(I,2).EQ.0) GOTO 550  
1717       P(I,5)=ULMASS(K(I,2)) 
1718       PR(JR)=P(I,5)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2  
1719
1720 C...Final two hadrons: find common setup of four-vectors.   
1721       JQ=1  
1722       IF(P(IN(4)+2,3)*P(IN(5)+2,3)*FOUR(IN(4),IN(5)).LT.P(IN(7),3)* 
1723      &P(IN(8),3)*FOUR(IN(7),IN(8))) JQ=2    
1724       DHC12=dble(FOUR(IN(3*JQ+1),IN(3*JQ+2)))
1725       DHR1=dble(FOUR(N+NRS,IN(3*JQ+2)))/DHC12
1726       DHR2=dble(FOUR(N+NRS,IN(3*JQ+1)))/DHC12
1727       IF(IN(4).NE.IN(7).OR.IN(5).NE.IN(8)) THEN 
1728         PX(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3))-PX(JQ) 
1729         PY(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3)+1)-PY(JQ)   
1730         PR(3-JQ)=P(I+(JT+JQ-3)**2-1,5)**2+(PX(3-JQ)+(2*JQ-3)*JS*    
1731      &  PX(3))**2+(PY(3-JQ)+(2*JQ-3)*JS*PY(3))**2   
1732       ENDIF 
1733     
1734 C...Solve kinematics for final two hadrons, if possible.    
1735       WREM2=WREM2+(PX(1)+PX(2))**2+(PY(1)+PY(2))**2 
1736       FD=(SQRT(PR(1))+SQRT(PR(2)))/SQRT(WREM2)  
1737       IF(MJU(1)+MJU(2).NE.0.AND.I.EQ.ISAV+2.AND.FD.GE.1.) GOTO 180  
1738       IF(FD.GE.1.) GOTO 550 
1739       FA=WREM2+PR(JT)-PR(JR)    
1740       IF(MSTJ(11).EQ.2) PREV=0.5*FD**PARJ(37+MSTJ(11))  
1741       IF(MSTJ(11).NE.2) PREV=0.5*EXP(MAX(-100.,LOG(FD)* 
1742      &PARJ(37+MSTJ(11))*(PR(1)+PR(2))**2))  
1743       FB=SIGN(SQRT(MAX(0.,FA**2-4.*WREM2*PR(JT))),JS*(RLU(0)-PREV)) 
1744       KFL1A=IABS(KFL(1))    
1745       KFL2A=IABS(KFL(2))    
1746       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),    
1747      &MOD(KFL2A/1000,10)).GE.6) FB=SIGN(SQRT(MAX(0.,FA**2-  
1748      &4.*WREM2*PR(JT))),FLOAT(JS))  
1749       DO 820 J=1,4  
1750       P(I-1,J)=(PX(JT)+PX(3))*P(IN(3*JQ+3),J)+(PY(JT)+PY(3))*   
1751      &P(IN(3*JQ+3)+1,J)+0.5*(sngl(DHR1)*(FA+FB)*P(IN(3*JQ+1),J)+  
1752      &sngl(DHR2)*(FA-FB)*P(IN(3*JQ+2),J))/WREM2   
1753   820 P(I,J)=P(N+NRS,J)-P(I-1,J)    
1754
1755 C...Mark jets as fragmented and give daughter pointers. 
1756       N=I-NRS+1 
1757       DO 830 I=NSAV+1,NSAV+NP   
1758       IM=K(I,3) 
1759       K(IM,1)=K(IM,1)+10    
1760       IF(MSTU(16).NE.2) THEN    
1761         K(IM,4)=NSAV+1  
1762         K(IM,5)=NSAV+1  
1763       ELSE  
1764         K(IM,4)=NSAV+2  
1765         K(IM,5)=N   
1766       ENDIF 
1767   830 CONTINUE  
1768     
1769 C...Document string system. Move up particles.  
1770       NSAV=NSAV+1   
1771       K(NSAV,1)=11  
1772       K(NSAV,2)=92  
1773       K(NSAV,3)=IP  
1774       K(NSAV,4)=NSAV+1  
1775       K(NSAV,5)=N   
1776       DO 840 J=1,4  
1777       P(NSAV,J)=sngl(DPS(J))
1778   840 V(NSAV,J)=V(IP,J) 
1779       P(NSAV,5)=SQRT(sngl(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2
1780      &     -DPS(3)**2)))
1781       V(NSAV,5)=0.
1782       DO 850 I=NSAV+1,N 
1783
1784       DO 850 J=1,5  
1785       K(I,J)=K(I+NRS-1,J)   
1786       P(I,J)=P(I+NRS-1,J)   
1787   850 V(I,J)=0. 
1788     
1789 C...Order particles in rank along the chain. Update mother pointer. 
1790       DO 860 I=NSAV+1,N 
1791       DO 860 J=1,5  
1792       K(I-NSAV+N,J)=K(I,J)  
1793   860 P(I-NSAV+N,J)=P(I,J)  
1794       I1=NSAV   
1795       DO 880 I=N+1,2*N-NSAV 
1796       IF(K(I,3).NE.IE(1)) GOTO 880  
1797       I1=I1+1   
1798       DO 870 J=1,5  
1799       K(I1,J)=K(I,J)    
1800   870 P(I1,J)=P(I,J)    
1801       IF(MSTU(16).NE.2) K(I1,3)=NSAV    
1802   880 CONTINUE  
1803       DO 900 I=2*N-NSAV,N+1,-1  
1804       IF(K(I,3).EQ.IE(1)) GOTO 900  
1805       I1=I1+1   
1806       DO 890 J=1,5  
1807       K(I1,J)=K(I,J)    
1808   890 P(I1,J)=P(I,J)    
1809       IF(MSTU(16).NE.2) K(I1,3)=NSAV    
1810   900 CONTINUE  
1811     
1812 C...Boost back particle system. Set production vertices.    
1813       CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),   
1814      &DPS(3)/DPS(4))    
1815       DO 910 I=NSAV+1,N 
1816
1817       DO 910 J=1,4  
1818   910 V(I,J)=V(IP,J)    
1819     
1820       RETURN    
1821       END   
1822     
1823 C*********************************************************************  
1824     
1825       SUBROUTINE LUINDF(IP) 
1826     
1827 C...Purpose: to handle the fragmentation of a jet system (or a single   
1828 C...jet) according to independent fragmentation models. 
1829       IMPLICIT DOUBLE PRECISION(D)  
1830       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
1831       SAVE /LUJETSA/ 
1832       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
1833       SAVE /LUDAT1A/ 
1834       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
1835       SAVE /LUDAT2A/ 
1836       DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),    
1837      &KFLO(2),PXO(2),PYO(2),WO(2)   
1838
1839       pw=0.
1840 C...Reset counters. Identify parton system and take copy. Check flavour.    
1841       NSAV=N    
1842       NJET=0    
1843       KQSUM=0   
1844       DO 100 J=1,5  
1845   100 DPS(J)=0.d0
1846       I=IP-1    
1847   110 I=I+1 
1848       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN 
1849         CALL LUERRM(12,'(LUINDF:) failed to reconstruct jet system')    
1850         IF(MSTU(21).GE.1) RETURN    
1851       ENDIF 
1852       IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110  
1853       KC=LUCOMP(K(I,2)) 
1854       IF(KC.EQ.0) GOTO 110  
1855       KQ=KCHG(KC,2)*ISIGN(1,K(I,2)) 
1856       IF(KQ.EQ.0) GOTO 110  
1857       NJET=NJET+1   
1858       IF(KQ.NE.2) KQSUM=KQSUM+KQ    
1859       DO 120 J=1,5  
1860       K(NSAV+NJET,J)=K(I,J) 
1861       P(NSAV+NJET,J)=P(I,J) 
1862   120 DPS(J)=DPS(J)+dble(P(I,J))
1863       K(NSAV+NJET,3)=I  
1864       IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.   
1865      &K(I+1,1).EQ.2)) GOTO 110  
1866       IF(NJET.NE.1.AND.KQSUM.NE.0) THEN 
1867         CALL LUERRM(12,'(LUINDF:) unphysical flavour combination')  
1868         IF(MSTU(21).GE.1) RETURN    
1869       ENDIF 
1870     
1871 C...Boost copied system to CM frame. Find CM energy and sum flavours.   
1872       IF(NJET.NE.1) CALL LUDBRB(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)/DPS(4),  
1873      &-DPS(2)/DPS(4),-DPS(3)/DPS(4))    
1874       PECM=0.   
1875       DO 130 J=1,3  
1876   130 NFI(J)=0  
1877       DO 140 I=NSAV+1,NSAV+NJET 
1878       PECM=PECM+P(I,4)  
1879       KFA=IABS(K(I,2))  
1880       IF(KFA.LE.3) THEN 
1881         NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))   
1882       ELSEIF(KFA.GT.1000) THEN  
1883         KFLA=MOD(KFA/1000,10)   
1884         KFLB=MOD(KFA/100,10)    
1885         IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))   
1886         IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))   
1887       ENDIF 
1888   140 CONTINUE  
1889     
1890 C...Loop over attempts made. Reset counters.    
1891       NTRY=0    
1892   150 NTRY=NTRY+1   
1893       N=NSAV+NJET   
1894       IF(NTRY.GT.200) THEN  
1895         CALL LUERRM(14,'(LUINDF:) caught in infinite loop') 
1896         IF(MSTU(21).GE.1) RETURN    
1897       ENDIF 
1898       DO 160 J=1,3  
1899       NFL(J)=NFI(J) 
1900       IFET(J)=0 
1901   160 KFLF(J)=0 
1902     
1903 C...Loop over jets to be fragmented.    
1904       DO 230 IP1=NSAV+1,NSAV+NJET   
1905       MSTJ(91)=0    
1906       NSAV1=N   
1907     
1908 C...Initial flavour and momentum values. Jet along +z axis. 
1909       KFLH=IABS(K(IP1,2))   
1910       IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10) 
1911       KFLO(2)=0 
1912       WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2) 
1913     
1914 C...Initial values for quark or diquark jet.    
1915   170 IF(IABS(K(IP1,2)).NE.21) THEN 
1916         NSTR=1  
1917         KFLO(1)=K(IP1,2)    
1918         CALL LUPTDI(0,PXO(1),PYO(1))    
1919         WO(1)=WF    
1920     
1921 C...Initial values for gluon treated like random quark jet. 
1922       ELSEIF(MSTJ(2).LE.2) THEN 
1923         NSTR=1  
1924         IF(MSTJ(2).EQ.2) MSTJ(91)=1 
1925         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)   
1926         CALL LUPTDI(0,PXO(1),PYO(1))    
1927         WO(1)=WF    
1928     
1929 C...Initial values for gluon treated like quark-antiquark jet pair, 
1930 C...sharing energy according to Altarelli-Parisi splitting function.    
1931       ELSE  
1932         NSTR=2  
1933         IF(MSTJ(2).EQ.4) MSTJ(91)=1 
1934         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)   
1935         KFLO(2)=-KFLO(1)    
1936         CALL LUPTDI(0,PXO(1),PYO(1))    
1937         PXO(2)=-PXO(1)  
1938         PYO(2)=-PYO(1)  
1939         WO(1)=WF*RLU(0)**(1./3.)    
1940         WO(2)=WF-WO(1)  
1941       ENDIF 
1942     
1943 C...Initial values for rank, flavour, pT and W+.    
1944       DO 220 ISTR=1,NSTR    
1945   180 I=N   
1946       IRANK=0   
1947       KFL1=KFLO(ISTR)   
1948       PX1=PXO(ISTR) 
1949       PY1=PYO(ISTR) 
1950       W=WO(ISTR)    
1951     
1952 C...New hadron. Generate flavour and hadron species.    
1953   190 I=I+1 
1954       IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN 
1955         CALL LUERRM(11,'(LUINDF:) no more memory left in LUJETSA')   
1956         IF(MSTU(21).GE.1) RETURN    
1957       ENDIF 
1958       IRANK=IRANK+1 
1959       K(I,1)=1  
1960       K(I,3)=IP1    
1961       K(I,4)=0  
1962       K(I,5)=0  
1963   200 CALL LUKFDI(KFL1,0,KFL2,K(I,2))   
1964       IF(K(I,2).EQ.0) GOTO 180  
1965       IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND. 
1966      &IABS(KFL2).GT.10) THEN    
1967         IF(RLU(0).GT.PARJ(19)) GOTO 200 
1968       ENDIF 
1969     
1970 C...Find hadron mass. Generate four-momentum.   
1971       P(I,5)=ULMASS(K(I,2)) 
1972       CALL LUPTDI(KFL1,PX2,PY2) 
1973       P(I,1)=PX1+PX2    
1974       P(I,2)=PY1+PY2    
1975       PR=P(I,5)**2+P(I,1)**2+P(I,2)**2  
1976       CALL LUZDIS(KFL1,KFL2,PR,Z)   
1977       P(I,3)=0.5*(Z*W-PR/(Z*W)) 
1978       P(I,4)=0.5*(Z*W+PR/(Z*W)) 
1979       IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND. 
1980      &P(I,3).LE.0.001) THEN 
1981         IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180   
1982         P(I,3)=0.0001   
1983         P(I,4)=SQRT(PR) 
1984         Z=P(I,4)/W  
1985       ENDIF 
1986     
1987 C...Remaining flavour and momentum. 
1988       KFL1=-KFL2    
1989       PX1=-PX2  
1990       PY1=-PY2  
1991       W=(1.-Z)*W    
1992       DO 210 J=1,5  
1993   210 V(I,J)=0. 
1994     
1995 C...Check if pL acceptable. Go back for new hadron if enough energy.    
1996       IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) I=I-1   
1997       IF(W.GT.PARJ(31)) GOTO 190    
1998   220 N=I   
1999       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32) 
2000       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170   
2001     
2002 C...Rotate jet to new direction.    
2003       THE=ULANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))    
2004       PHI=ULANGL(P(IP1,1),P(IP1,2)) 
2005       CALL LUDBRB(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)    
2006       K(K(IP1,3),4)=NSAV1+1 
2007       K(K(IP1,3),5)=N   
2008     
2009 C...End of jet generation loop. Skip conservation in some cases.    
2010   230 CONTINUE  
2011       IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 470    
2012       IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150 
2013     
2014 C...Subtract off produced hadron flavours, finished if zero.    
2015       DO 240 I=NSAV+NJET+1,N    
2016       KFA=IABS(K(I,2))  
2017       KFLA=MOD(KFA/1000,10) 
2018       KFLB=MOD(KFA/100,10)  
2019       KFLC=MOD(KFA/10,10)   
2020       IF(KFLA.EQ.0) THEN    
2021         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB    
2022         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB    
2023       ELSE  
2024         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))   
2025         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))   
2026         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))   
2027       ENDIF 
2028   240 CONTINUE  
2029       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+ 
2030      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3    
2031       IF(NREQ.EQ.0) GOTO 320    
2032     
2033 C...Take away flavour of low-momentum particles until enough freedom.   
2034       NREM=0    
2035   250 IREM=0    
2036       P2MIN=PECM**2 
2037       DO 260 I=NSAV+NJET+1,N    
2038       P2=P(I,1)**2+P(I,2)**2+P(I,3)**2  
2039       IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I    
2040   260 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2  
2041       IF(IREM.EQ.0) GOTO 150    
2042       K(IREM,1)=7   
2043       KFA=IABS(K(IREM,2))   
2044       KFLA=MOD(KFA/1000,10) 
2045       KFLB=MOD(KFA/100,10)  
2046       KFLC=MOD(KFA/10,10)   
2047       IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8    
2048       IF(K(IREM,1).EQ.8) GOTO 250   
2049       IF(KFLA.EQ.0) THEN    
2050         ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB  
2051         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN  
2052         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN  
2053       ELSE  
2054         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))    
2055         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))    
2056         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))    
2057       ENDIF 
2058       NREM=NREM+1   
2059       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+ 
2060      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3    
2061       IF(NREQ.GT.NREM) GOTO 250 
2062       DO 270 I=NSAV+NJET+1,N    
2063   270 IF(K(I,1).EQ.8) K(I,1)=1  
2064     
2065 C...Find combination of existing and new flavours for hadron.   
2066   280 NFET=2    
2067       IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3  
2068       IF(NREQ.LT.NREM) NFET=1   
2069       IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0    
2070       DO 290 J=1,NFET   
2071       IFET(J)=1+int((IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU(0))
2072       KFLF(J)=ISIGN(1,NFL(1))   
2073       IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))   
2074   290 IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))  
2075       IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))    
2076      &GOTO 280  
2077       IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.    
2078      &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3).    
2079      &LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280    
2080       IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU(0))  
2081       IF(NFET.EQ.0) KFLF(2)=-KFLF(1)    
2082       IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU(0)),-KFLF(1))  
2083       IF(NFET.LE.2) KFLF(3)=0   
2084       IF(KFLF(3).NE.0) THEN 
2085         KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+  
2086      &  100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1)) 
2087         IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU(0).GT.1.)  
2088      &  KFLFC=KFLFC+ISIGN(2,KFLFC)  
2089       ELSE  
2090         KFLFC=KFLF(1)   
2091       ENDIF 
2092       CALL LUKFDI(KFLFC,KFLF(2),KFLDMP,KF)  
2093       IF(KF.EQ.0) GOTO 280  
2094       DO 300 J=1,MAX(2,NFET)    
2095   300 NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))    
2096     
2097 C...Store hadron at random among free positions.    
2098       NPOS=MIN(1+INT(RLU(0)*NREM),NREM) 
2099       DO 310 I=NSAV+NJET+1,N    
2100       IF(K(I,1).EQ.7) NPOS=NPOS-1   
2101       IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310 
2102       K(I,1)=1  
2103       K(I,2)=KF 
2104       P(I,5)=ULMASS(K(I,2)) 
2105       P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)  
2106   310 CONTINUE  
2107       NREM=NREM-1   
2108       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+ 
2109      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3    
2110       IF(NREM.GT.0) GOTO 280    
2111     
2112 C...Compensate for missing momentum in global scheme (3 options).   
2113   320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN  
2114         DO 330 J=1,3    
2115         PSI(J)=0.   
2116         DO 330 I=NSAV+NJET+1,N  
2117   330   PSI(J)=PSI(J)+P(I,J)    
2118         PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2    
2119         PWS=0.  
2120         DO 340 I=NSAV+NJET+1,N  
2121         IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)  
2122         IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+  
2123      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4)) 
2124   340   IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1.  
2125         DO 360 I=NSAV+NJET+1,N  
2126         IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)   
2127         IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+   
2128      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4)) 
2129         IF(MOD(MSTJ(3),5).EQ.3) PW=1.   
2130         DO 350 J=1,3    
2131   350   P(I,J)=P(I,J)-PSI(J)*PW/PWS 
2132   360   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)    
2133     
2134 C...Compensate for missing momentum withing each jet separately.    
2135       ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN  
2136         DO 370 I=N+1,N+NJET 
2137         K(I,1)=0    
2138         DO 370 J=1,5    
2139   370   P(I,J)=0.   
2140         DO 390 I=NSAV+NJET+1,N  
2141         IR1=K(I,3)  
2142         IR2=N+IR1-NSAV  
2143         K(IR2,1)=K(IR2,1)+1 
2144         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/  
2145      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)   
2146         DO 380 J=1,3    
2147   380   P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)   
2148         P(IR2,4)=P(IR2,4)+P(I,4)    
2149   390   P(IR2,5)=P(IR2,5)+PLS   
2150         PSS=0.  
2151         DO 400 I=N+1,N+NJET 
2152   400   IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2))  
2153         DO 420 I=NSAV+NJET+1,N  
2154         IR1=K(I,3)  
2155         IR2=N+IR1-NSAV  
2156         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/  
2157      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)   
2158         DO 410 J=1,3    
2159   410   P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS* 
2160      &  P(IR1,J)    
2161   420   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)    
2162       ENDIF 
2163     
2164 C...Scale momenta for energy conservation.  
2165       IF(MOD(MSTJ(3),5).NE.0) THEN  
2166         PMS=0.  
2167         PES=0.  
2168         PQS=0.  
2169         DO 430 I=NSAV+NJET+1,N  
2170         PMS=PMS+P(I,5)  
2171         PES=PES+P(I,4)  
2172   430   PQS=PQS+P(I,5)**2/P(I,4)    
2173         IF(PMS.GE.PECM) GOTO 150    
2174         NECO=0  
2175   440   NECO=NECO+1 
2176         PFAC=(PECM-PQS)/(PES-PQS)   
2177         PES=0.  
2178         PQS=0.  
2179         DO 460 I=NSAV+NJET+1,N  
2180         DO 450 J=1,3    
2181   450   P(I,J)=PFAC*P(I,J)  
2182         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)    
2183         PES=PES+P(I,4)  
2184   460   PQS=PQS+P(I,5)**2/P(I,4)    
2185         IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 440  
2186       ENDIF 
2187     
2188 C...Origin of produced particles and parton daughter pointers.  
2189   470 DO 480 I=NSAV+NJET+1,N    
2190       IF(MSTU(16).NE.2) K(I,3)=NSAV+1   
2191   480 IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)  
2192       DO 490 I=NSAV+1,NSAV+NJET 
2193       I1=K(I,3) 
2194       K(I1,1)=K(I1,1)+10    
2195       IF(MSTU(16).NE.2) THEN    
2196         K(I1,4)=NSAV+1  
2197         K(I1,5)=NSAV+1  
2198       ELSE  
2199         K(I1,4)=K(I1,4)-NJET+1  
2200         K(I1,5)=K(I1,5)-NJET+1  
2201         IF(K(I1,5).LT.K(I1,4)) THEN 
2202           K(I1,4)=0 
2203           K(I1,5)=0 
2204         ENDIF   
2205       ENDIF 
2206   490 CONTINUE  
2207     
2208 C...Document independent fragmentation system. Remove copy of jets. 
2209       NSAV=NSAV+1   
2210       K(NSAV,1)=11  
2211       K(NSAV,2)=93  
2212       K(NSAV,3)=IP  
2213       K(NSAV,4)=NSAV+1  
2214       K(NSAV,5)=N-NJET+1    
2215       DO 500 J=1,4  
2216       P(NSAV,J)=sngl(DPS(J))
2217   500 V(NSAV,J)=V(IP,J) 
2218       P(NSAV,5)=SQRT(sngl(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2
2219      &     -DPS(3)**2)))
2220       V(NSAV,5)=0.  
2221       DO 510 I=NSAV+NJET,N  
2222       DO 510 J=1,5  
2223       K(I-NJET+1,J)=K(I,J)  
2224       P(I-NJET+1,J)=P(I,J)  
2225   510 V(I-NJET+1,J)=V(I,J)  
2226       N=N-NJET+1    
2227     
2228 C...Boost back particle system. Set production vertices.    
2229       IF(NJET.NE.1) CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),   
2230      &DPS(2)/DPS(4),DPS(3)/DPS(4))  
2231       DO 520 I=NSAV+1,N 
2232       DO 520 J=1,4  
2233   520 V(I,J)=V(IP,J)    
2234     
2235       RETURN    
2236       END   
2237     
2238 C*********************************************************************  
2239     
2240       SUBROUTINE LUDECY(IP) 
2241     
2242 C...Purpose: to handle the decay of unstable particles. 
2243       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
2244       SAVE /LUJETSA/ 
2245       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
2246       SAVE /LUDAT1A/ 
2247       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
2248       SAVE /LUDAT2A/ 
2249       COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)    
2250       SAVE /LUDAT3A/ 
2251       DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3),  
2252      &WTCOR(10) 
2253 clin-2/18/03 for resonance decay in hadron cascade:
2254       common/resdcy/NSAV,iksdcy
2255       SAVE /resdcy/
2256       DATA WTCOR/2.,5.,15.,60.,250.,1500.,1.2E4,1.2E5,150.,16./ 
2257     
2258 C...Functions: momentum in two-particle decays, four-product and    
2259 C...matrix element times phase space in weak decays.    
2260       PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A)  
2261       FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3) 
2262       HMEPS(HA)=((1.-HRQ-HA)**2+3.*HA*(1.+HRQ-HA))* 
2263      &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA)    
2264     
2265 C...Initial values. 
2266
2267       idc=0
2268       pqt=0.
2269       hatu=0.
2270       hmp1=0.
2271       im=0
2272       kfam=0
2273       wtmax=0.
2274       pmes=0.
2275       pmst=0.
2276       wt=0.
2277       pmr=0.
2278
2279       NTRY=0    
2280       NSAV=N    
2281       KFA=IABS(K(IP,2)) 
2282       KFS=ISIGN(1,K(IP,2))  
2283       KC=LUCOMP(KFA)    
2284       MSTJ(92)=0    
2285     
2286 C...Choose lifetime and determine decay vertex. 
2287       IF(K(IP,1).EQ.5) THEN 
2288         V(IP,5)=0.  
2289       ELSEIF(K(IP,1).NE.4) THEN 
2290         V(IP,5)=-PMAS(KC,4)*LOG(RLU(0)) 
2291       ENDIF 
2292       DO 100 J=1,4  
2293   100 VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)   
2294     
2295 C...Determine whether decay allowed or not. 
2296       MOUT=0    
2297       IF(MSTJ(22).EQ.2) THEN    
2298         IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1   
2299       ELSEIF(MSTJ(22).EQ.3) THEN    
2300         IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1  
2301       ELSEIF(MSTJ(22).EQ.4) THEN    
2302         IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1 
2303         IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1 
2304       ENDIF 
2305       IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN   
2306         K(IP,1)=4   
2307         RETURN  
2308       ENDIF 
2309     
2310 C...Check existence of decay channels. Particle/antiparticle rules. 
2311       KCA=KC    
2312       IF(MDCY(KC,2).GT.0) THEN  
2313         MDMDCY=MDME(MDCY(KC,2),2)   
2314         IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY    
2315       ENDIF 
2316       IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN 
2317         CALL LUERRM(9,'(LUDECY:) no decay channel defined') 
2318         RETURN  
2319       ENDIF 
2320       IF(MOD(KFA/1000,10).EQ.0.AND.(KCA.EQ.85.OR.KCA.EQ.87)) KFS=-KFS   
2321       IF(KCHG(KC,3).EQ.0) THEN  
2322         KFSP=1  
2323         KFSN=0  
2324         IF(RLU(0).GT.0.5) KFS=-KFS  
2325       ELSEIF(KFS.GT.0) THEN 
2326         KFSP=1  
2327         KFSN=0  
2328       ELSE  
2329         KFSP=0  
2330         KFSN=1  
2331       ENDIF 
2332     
2333 C...Sum branching ratios of allowed decay channels. 
2334 clin  110 NOPE=0    
2335       NOPE=0    
2336       BRSU=0.   
2337       DO 120 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1  
2338       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.    
2339      &KFSN*MDME(IDL,1).NE.3) GOTO 120   
2340       IF(MDME(IDL,2).GT.100) GOTO 120   
2341       NOPE=NOPE+1   
2342       BRSU=BRSU+BRAT(IDL)   
2343   120 CONTINUE  
2344       IF(NOPE.EQ.0) THEN    
2345         CALL LUERRM(2,'(LUDECY:) all decay channels closed by user')    
2346         RETURN  
2347       ENDIF 
2348     
2349 C...Select decay channel among allowed ones.    
2350   130 RBR=BRSU*RLU(0)   
2351       IDL=MDCY(KCA,2)-1 
2352   140 IDL=IDL+1 
2353       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.    
2354      &KFSN*MDME(IDL,1).NE.3) THEN   
2355         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140   
2356       ELSEIF(MDME(IDL,2).GT.100) THEN   
2357         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140   
2358       ELSE  
2359         IDC=IDL 
2360         RBR=RBR-BRAT(IDL)   
2361         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0.) GOTO 140 
2362       ENDIF 
2363     
2364 C...Start readout of decay channel: matrix element, reset counters. 
2365       MMAT=MDME(IDC,2)  
2366   150 NTRY=NTRY+1   
2367       IF(NTRY.GT.1000) THEN 
2368         CALL LUERRM(14,'(LUDECY:) caught in infinite loop') 
2369         IF(MSTU(21).GE.1) RETURN    
2370       ENDIF 
2371       I=N   
2372       NP=0  
2373       NQ=0  
2374       MBST=0    
2375       IF(MMAT.GE.11.AND.MMAT.NE.46.AND.P(IP,4).GT.20.*P(IP,5)) MBST=1   
2376       DO 160 J=1,4  
2377       PV(1,J)=0.    
2378   160 IF(MBST.EQ.0) PV(1,J)=P(IP,J) 
2379       IF(MBST.EQ.1) PV(1,4)=P(IP,5) 
2380       PV(1,5)=P(IP,5)   
2381       PS=0. 
2382       PSQ=0.    
2383       MREM=0    
2384     
2385 C...Read out decay products. Convert to standard flavour code.  
2386       JTMAX=5   
2387       IF(MDME(IDC+1,2).EQ.101) JTMAX=10 
2388       DO 170 JT=1,JTMAX 
2389       IF(JT.LE.5) KP=KFDP(IDC,JT)   
2390       IF(JT.GE.6) KP=KFDP(IDC+1,JT-5)   
2391       IF(KP.EQ.0) GOTO 170  
2392       KPA=IABS(KP)  
2393       KCP=LUCOMP(KPA)   
2394       IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN 
2395         KFP=KP  
2396       ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN  
2397         KFP=KFS*KP  
2398       ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN  
2399         KFP=-KFS*MOD(KFA/10,10) 
2400       ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN  
2401         KFP=KFS*(100*MOD(KFA/10,100)+3) 
2402       ELSEIF(KPA.EQ.81) THEN    
2403         KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1) 
2404       ELSEIF(KP.EQ.82) THEN 
2405         CALL LUKFDI(-KFS*INT(1.+(2.+PARJ(2))*RLU(0)),0,KFP,KDUMP)   
2406         IF(KFP.EQ.0) GOTO 150   
2407         MSTJ(93)=1  
2408         IF(PV(1,5).LT.PARJ(32)+2.*ULMASS(KFP)) GOTO 150 
2409       ELSEIF(KP.EQ.-82) THEN    
2410         KFP=-KFP    
2411         IF(IABS(KFP).GT.10) KFP=KFP+ISIGN(10000,KFP)    
2412       ENDIF 
2413       IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=LUCOMP(KFP)    
2414     
2415 C...Add decay product to event record or to quark flavour list. 
2416       KFPA=IABS(KFP)    
2417       KQP=KCHG(KCP,2)   
2418       IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN   
2419         NQ=NQ+1 
2420         KFLO(NQ)=KFP    
2421         MSTJ(93)=2  
2422         PSQ=PSQ+ULMASS(KFLO(NQ))    
2423       ELSEIF(MMAT.GE.42.AND.MMAT.LE.43.AND.NP.EQ.3.AND.MOD(NQ,2).EQ.1)  
2424      &THEN  
2425         NQ=NQ-1 
2426         PS=PS-P(I,5)    
2427         K(I,1)=1    
2428         KFI=K(I,2)  
2429         CALL LUKFDI(KFP,KFI,KFLDMP,K(I,2))  
2430         IF(K(I,2).EQ.0) GOTO 150    
2431         MSTJ(93)=1  
2432         P(I,5)=ULMASS(K(I,2))   
2433         PS=PS+P(I,5)    
2434       ELSE  
2435         I=I+1   
2436         NP=NP+1 
2437         IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1 
2438         IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1    
2439         K(I,1)=1+MOD(NQ,2)  
2440         IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2    
2441         IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1  
2442         K(I,2)=KFP  
2443         K(I,3)=IP   
2444         K(I,4)=0    
2445         K(I,5)=0    
2446         P(I,5)=ULMASS(KFP)  
2447         IF(MMAT.EQ.45.AND.KFPA.EQ.89) P(I,5)=PARJ(32)   
2448         PS=PS+P(I,5)    
2449       ENDIF 
2450   170 CONTINUE  
2451     
2452 C...Choose decay multiplicity in phase space model. 
2453   180 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN    
2454         PSP=PS  
2455         CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1))   
2456         IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63)   
2457   190   NTRY=NTRY+1 
2458         IF(NTRY.GT.1000) THEN   
2459           CALL LUERRM(14,'(LUDECY:) caught in infinite loop')   
2460           IF(MSTU(21).GE.1) RETURN  
2461         ENDIF   
2462         IF(MMAT.LE.20) THEN 
2463           GAUSS=SQRT(-2.*CNDE*LOG(MAX(1E-10,RLU(0))))*  
2464      &    SIN(PARU(2)*RLU(0))   
2465           ND=int(0.5+0.5*NP+0.25*NQ+CNDE+GAUSS)
2466           IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 190 
2467           IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 190   
2468           IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 190   
2469           IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 190   
2470         ELSE    
2471           ND=MMAT-20    
2472         ENDIF   
2473     
2474 C...Form hadrons from flavour content.  
2475         DO 200 JT=1,4   
2476   200   KFL1(JT)=KFLO(JT)   
2477         IF(ND.EQ.NP+NQ/2) GOTO 220  
2478         DO 210 I=N+NP+1,N+ND-NQ/2   
2479         JT=1+INT((NQ-1)*RLU(0)) 
2480         CALL LUKFDI(KFL1(JT),0,KFL2,K(I,2)) 
2481         IF(K(I,2).EQ.0) GOTO 190    
2482   210   KFL1(JT)=-KFL2  
2483   220   JT=2    
2484         JT2=3   
2485         JT3=4   
2486         IF(NQ.EQ.4.AND.RLU(0).LT.PARJ(66)) JT=4 
2487         IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))* 
2488      &  ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3    
2489         IF(JT.EQ.3) JT2=2   
2490         IF(JT.EQ.4) JT3=2   
2491         CALL LUKFDI(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2))   
2492         IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 190  
2493         IF(NQ.EQ.4) CALL LUKFDI(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND,2))   
2494         IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 190 
2495     
2496 C...Check that sum of decay product masses not too large.   
2497         PS=PSP  
2498         DO 230 I=N+NP+1,N+ND    
2499         K(I,1)=1    
2500         K(I,3)=IP   
2501         K(I,4)=0    
2502         K(I,5)=0    
2503         P(I,5)=ULMASS(K(I,2))   
2504   230   PS=PS+P(I,5)    
2505         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 190 
2506     
2507 C...Rescale energy to subtract off spectator quark mass.    
2508       ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44.OR.MMAT.EQ.45).    
2509      &AND.NP.GE.3) THEN 
2510         PS=PS-P(N+NP,5) 
2511         PQT=(P(N+NP,5)+PARJ(65))/PV(1,5)    
2512         DO 240 J=1,5    
2513         P(N+NP,J)=PQT*PV(1,J)   
2514   240   PV(1,J)=(1.-PQT)*PV(1,J)    
2515         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 150 
2516         ND=NP-1 
2517         MREM=1  
2518     
2519 C...Phase space factors imposed in W decay. 
2520       ELSEIF(MMAT.EQ.46) THEN   
2521         MSTJ(93)=1  
2522         PSMC=ULMASS(K(N+1,2))   
2523         MSTJ(93)=1  
2524         PSMC=PSMC+ULMASS(K(N+2,2))  
2525         IF(MAX(PS,PSMC)+PARJ(32).GT.PV(1,5)) GOTO 130   
2526         HR1=(P(N+1,5)/PV(1,5))**2   
2527         HR2=(P(N+2,5)/PV(1,5))**2   
2528         IF((1.-HR1-HR2)*(2.+HR1+HR2)*SQRT((1.-HR1-HR2)**2-4.*HR1*HR2).  
2529      &  LT.2.*RLU(0)) GOTO 130  
2530         ND=NP   
2531     
2532 C...Fully specified final state: check mass broadening effects. 
2533       ELSE  
2534         IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 150 
2535         ND=NP   
2536       ENDIF 
2537     
2538 C...Select W mass in decay Q -> W + q, without W propagator.    
2539       IF(MMAT.EQ.45.AND.MSTJ(25).LE.0) THEN 
2540         HLQ=(PARJ(32)/PV(1,5))**2   
2541         HUQ=(1.-(P(N+2,5)+PARJ(64))/PV(1,5))**2 
2542         HRQ=(P(N+2,5)/PV(1,5))**2   
2543   250   HW=HLQ+RLU(0)*(HUQ-HLQ) 
2544         IF(HMEPS(HW).LT.RLU(0)) GOTO 250    
2545         P(N+1,5)=PV(1,5)*SQRT(HW)   
2546     
2547 C...Ditto, including W propagator. Divide mass range into three regions.    
2548       ELSEIF(MMAT.EQ.45) THEN   
2549         HQW=(PV(1,5)/PMAS(24,1))**2 
2550         HLW=(PARJ(32)/PMAS(24,1))**2    
2551         HUW=((PV(1,5)-P(N+2,5)-PARJ(64))/PMAS(24,1))**2 
2552         HRQ=(P(N+2,5)/PV(1,5))**2   
2553         HG=PMAS(24,2)/PMAS(24,1)    
2554         HATL=ATAN((HLW-1.)/HG)  
2555         HM=MIN(1.,HUW-0.001)    
2556         HMV1=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)   
2557   260   HM=HM-HG    
2558         HMV2=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)   
2559         HSAV1=HMEPS(HM/HQW) 
2560         HSAV2=1./((HM-1.)**2+HG**2) 
2561         IF(HMV2.GT.HMV1.AND.HM-HG.GT.HLW) THEN  
2562           HMV1=HMV2 
2563           GOTO 260  
2564         ENDIF   
2565         HMV=MIN(2.*HMV1,HMEPS(HM/HQW)/HG**2)    
2566         HM1=1.-SQRT(1./HMV-HG**2)   
2567         IF(HM1.GT.HLW.AND.HM1.LT.HM) THEN   
2568           HM=HM1    
2569         ELSEIF(HMV2.LE.HMV1) THEN   
2570           HM=MAX(HLW,HM-MIN(0.1,1.-HM)) 
2571         ENDIF   
2572         HATM=ATAN((HM-1.)/HG)   
2573         HWT1=(HATM-HATL)/HG 
2574         HWT2=HMV*(MIN(1.,HUW)-HM)   
2575         HWT3=0. 
2576         IF(HUW.GT.1.) THEN  
2577           HATU=ATAN((HUW-1.)/HG)    
2578           HMP1=HMEPS(1./HQW)    
2579           HWT3=HMP1*HATU/HG 
2580         ENDIF   
2581     
2582 C...Select mass region and W mass there. Accept according to weight.    
2583   270   HREG=RLU(0)*(HWT1+HWT2+HWT3)    
2584         IF(HREG.LE.HWT1) THEN   
2585           HW=1.+HG*TAN(HATL+RLU(0)*(HATM-HATL)) 
2586           HACC=HMEPS(HW/HQW)    
2587         ELSEIF(HREG.LE.HWT1+HWT2) THEN  
2588           HW=HM+RLU(0)*(MIN(1.,HUW)-HM) 
2589           HACC=HMEPS(HW/HQW)/((HW-1.)**2+HG**2)/HMV 
2590         ELSE    
2591           HW=1.+HG*TAN(RLU(0)*HATU) 
2592           HACC=HMEPS(HW/HQW)/HMP1   
2593         ENDIF   
2594         IF(HACC.LT.RLU(0)) GOTO 270 
2595         P(N+1,5)=PMAS(24,1)*SQRT(HW)    
2596       ENDIF 
2597     
2598 C...Determine position of grandmother, number of sisters, Q -> W sign.  
2599       NM=0  
2600       MSGN=0    
2601       IF(MMAT.EQ.3.OR.MMAT.EQ.46) THEN  
2602         IM=K(IP,3)  
2603         IF(IM.LT.0.OR.IM.GE.IP) IM=0    
2604         IF(IM.NE.0) KFAM=IABS(K(IM,2))  
2605         IF(IM.NE.0.AND.MMAT.EQ.3) THEN  
2606           DO 280 IL=MAX(IP-2,IM+1),MIN(IP+2,N)  
2607   280     IF(K(IL,3).EQ.IM) NM=NM+1 
2608           IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR.    
2609      &    MOD(KFAM/1000,10).NE.0) NM=0  
2610         ELSEIF(IM.NE.0.AND.MMAT.EQ.46) THEN 
2611           MSGN=ISIGN(1,K(IM,2)*K(IP,2)) 
2612           IF(KFAM.GT.100.AND.MOD(KFAM/1000,10).EQ.0) MSGN=  
2613      &    MSGN*(-1)**MOD(KFAM/100,10)   
2614         ENDIF   
2615       ENDIF 
2616     
2617 C...Kinematics of one-particle decays.  
2618       IF(ND.EQ.1) THEN  
2619         DO 290 J=1,4    
2620   290   P(N+1,J)=P(IP,J)    
2621         GOTO 510    
2622       ENDIF 
2623     
2624 C...Calculate maximum weight ND-particle decay. 
2625       PV(ND,5)=P(N+ND,5)    
2626       IF(ND.GE.3) THEN  
2627         WTMAX=1./WTCOR(ND-2)    
2628         PMAX=PV(1,5)-PS+P(N+ND,5)   
2629         PMIN=0. 
2630         DO 300 IL=ND-1,1,-1 
2631         PMAX=PMAX+P(N+IL,5) 
2632         PMIN=PMIN+P(N+IL+1,5)   
2633   300   WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5))   
2634       ENDIF 
2635     
2636 C...Find virtual gamma mass in Dalitz decay.    
2637   310 IF(ND.EQ.2) THEN  
2638       ELSEIF(MMAT.EQ.2) THEN    
2639         PMES=4.*PMAS(11,1)**2   
2640         PMRHO2=PMAS(131,1)**2   
2641         PGRHO2=PMAS(131,2)**2   
2642   320   PMST=PMES*(P(IP,5)**2/PMES)**RLU(0) 
2643         WT=(1+0.5*PMES/PMST)*SQRT(MAX(0.,1.-PMES/PMST))*    
2644      &  (1.-PMST/P(IP,5)**2)**3*(1.+PGRHO2/PMRHO2)/ 
2645      &  ((1.-PMST/PMRHO2)**2+PGRHO2/PMRHO2) 
2646         IF(WT.LT.RLU(0)) GOTO 320   
2647         PV(2,5)=MAX(2.00001*PMAS(11,1),SQRT(PMST))  
2648     
2649 C...M-generator gives weight. If rejected, try again.   
2650       ELSE  
2651   330   RORD(1)=1.  
2652         DO 350 IL1=2,ND-1   
2653         RSAV=RLU(0) 
2654         DO 340 IL2=IL1-1,1,-1   
2655         IF(RSAV.LE.RORD(IL2)) GOTO 350  
2656   340   RORD(IL2+1)=RORD(IL2)   
2657   350   RORD(IL2+1)=RSAV    
2658         RORD(ND)=0. 
2659         WT=1.   
2660         DO 360 IL=ND-1,1,-1 
2661         PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)    
2662   360   WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))   
2663         IF(WT.LT.RLU(0)*WTMAX) GOTO 330 
2664       ENDIF 
2665     
2666 C...Perform two-particle decays in respective CM frame. 
2667   370 DO 390 IL=1,ND-1  
2668       PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))    
2669       UE(3)=2.*RLU(0)-1.    
2670       PHI=PARU(2)*RLU(0)    
2671       UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)  
2672       UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)  
2673       DO 380 J=1,3  
2674       P(N+IL,J)=PA*UE(J)    
2675   380 PV(IL+1,J)=-PA*UE(J)  
2676       P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2)    
2677   390 PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)  
2678     
2679 C...Lorentz transform decay products to lab frame.  
2680       DO 400 J=1,4  
2681   400 P(N+ND,J)=PV(ND,J)    
2682       DO 430 IL=ND-1,1,-1   
2683       DO 410 J=1,3  
2684   410 BE(J)=PV(IL,J)/PV(IL,4)   
2685       GA=PV(IL,4)/PV(IL,5)  
2686       DO 430 I=N+IL,N+ND    
2687       BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)    
2688       DO 420 J=1,3  
2689   420 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)    
2690   430 P(I,4)=GA*(P(I,4)+BEP)    
2691     
2692 C...Matrix elements for omega and phi decays.   
2693       IF(MMAT.EQ.1) THEN    
2694         WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2  
2695      &  -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2    
2696      &  +2.*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3)   
2697         IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001).LT.RLU(0)) GOTO 310    
2698     
2699 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-. 
2700       ELSEIF(MMAT.EQ.2) THEN    
2701         FOUR12=FOUR(N+1,N+2)    
2702         FOUR13=FOUR(N+1,N+3)    
2703         FOUR23=0.5*PMST-0.25*PMES   
2704         WT=(PMST-0.5*PMES)*(FOUR12**2+FOUR13**2)+   
2705      &  PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2)    
2706         IF(WT.LT.RLU(0)*0.25*PMST*(P(IP,5)**2-PMST)**2) GOTO 370    
2707     
2708 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar, 
2709 C...V vector), of form cos**2(theta02) in V1 rest frame.    
2710       ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN    
2711         IF((P(IP,5)**2*FOUR(IM,N+1)-FOUR(IP,IM)*FOUR(IP,N+1))**2.LE.    
2712      &  RLU(0)*(FOUR(IP,IM)**2-(P(IP,5)*P(IM,5))**2)*(FOUR(IP,N+1)**2-  
2713      &  (P(IP,5)*P(N+1,5))**2)) GOTO 370    
2714     
2715 C...Matrix element for "onium" -> g + g + g or gamma + g + g.   
2716       ELSEIF(MMAT.EQ.4) THEN    
2717         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2  
2718         HX2=2.*FOUR(IP,N+2)/P(IP,5)**2  
2719         HX3=2.*FOUR(IP,N+3)/P(IP,5)**2  
2720         WT=((1.-HX1)/(HX2*HX3))**2+((1.-HX2)/(HX1*HX3))**2+ 
2721      &  ((1.-HX3)/(HX1*HX2))**2 
2722         IF(WT.LT.2.*RLU(0)) GOTO 310    
2723         IF(K(IP+1,2).EQ.22.AND.(1.-HX1)*P(IP,5)**2.LT.4.*PARJ(32)**2)   
2724      &  GOTO 310    
2725     
2726 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.    
2727       ELSEIF(MMAT.EQ.41) THEN   
2728         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2  
2729         IF(8.*HX1*(3.-2.*HX1)/9..LT.RLU(0)) GOTO 310    
2730     
2731 C...Matrix elements for weak decays (only semileptonic for c and b) 
2732       ELSEIF(MMAT.GE.42.AND.MMAT.LE.44.AND.ND.EQ.3) THEN    
2733         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3) 
2734         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3) 
2735         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310  
2736       ELSEIF(MMAT.GE.42.AND.MMAT.LE.44) THEN    
2737         DO 440 J=1,4    
2738         P(N+NP+1,J)=0.  
2739         DO 440 IS=N+3,N+NP  
2740   440   P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J) 
2741         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1)  
2742         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1)  
2743         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310  
2744     
2745 C...Angular distribution in W decay.    
2746       ELSEIF(MMAT.EQ.46.AND.MSGN.NE.0) THEN 
2747         IF(MSGN.GT.0) WT=FOUR(IM,N+1)*FOUR(N+2,IP+1)    
2748         IF(MSGN.LT.0) WT=FOUR(IM,N+2)*FOUR(N+1,IP+1)    
2749         IF(WT.LT.RLU(0)*P(IM,5)**4/WTCOR(10)) GOTO 370  
2750       ENDIF 
2751     
2752 C...Scale back energy and reattach spectator.   
2753       IF(MREM.EQ.1) THEN    
2754         DO 450 J=1,5    
2755   450   PV(1,J)=PV(1,J)/(1.-PQT)    
2756         ND=ND+1 
2757         MREM=0  
2758       ENDIF 
2759     
2760 C...Low invariant mass for system with spectator quark gives particle,  
2761 C...not two jets. Readjust momenta accordingly. 
2762       IF((MMAT.EQ.31.OR.MMAT.EQ.45).AND.ND.EQ.3) THEN   
2763         MSTJ(93)=1  
2764         PM2=ULMASS(K(N+2,2))    
2765         MSTJ(93)=1  
2766         PM3=ULMASS(K(N+3,2))    
2767         IF(P(N+2,5)**2+P(N+3,5)**2+2.*FOUR(N+2,N+3).GE. 
2768      &  (PARJ(32)+PM2+PM3)**2) GOTO 510 
2769         K(N+2,1)=1  
2770         KFTEMP=K(N+2,2) 
2771         CALL LUKFDI(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2))    
2772         IF(K(N+2,2).EQ.0) GOTO 150  
2773         P(N+2,5)=ULMASS(K(N+2,2))   
2774         PS=P(N+1,5)+P(N+2,5)    
2775         PV(2,5)=P(N+2,5)    
2776         MMAT=0  
2777         ND=2    
2778         GOTO 370    
2779       ELSEIF(MMAT.EQ.44) THEN   
2780         MSTJ(93)=1  
2781         PM3=ULMASS(K(N+3,2))    
2782         MSTJ(93)=1  
2783         PM4=ULMASS(K(N+4,2))    
2784         IF(P(N+3,5)**2+P(N+4,5)**2+2.*FOUR(N+3,N+4).GE. 
2785      &  (PARJ(32)+PM3+PM4)**2) GOTO 480 
2786         K(N+3,1)=1  
2787         KFTEMP=K(N+3,2) 
2788         CALL LUKFDI(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2))    
2789         IF(K(N+3,2).EQ.0) GOTO 150  
2790         P(N+3,5)=ULMASS(K(N+3,2))   
2791         DO 460 J=1,3    
2792   460   P(N+3,J)=P(N+3,J)+P(N+4,J)  
2793         P(N+3,4)=SQRT(P(N+3,1)**2+P(N+3,2)**2+P(N+3,3)**2+P(N+3,5)**2)  
2794         HA=P(N+1,4)**2-P(N+2,4)**2  
2795         HB=HA-(P(N+1,5)**2-P(N+2,5)**2) 
2796         HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+   
2797      &  (P(N+1,3)-P(N+2,3))**2  
2798         HD=(PV(1,4)-P(N+3,4))**2    
2799         HE=HA**2-2.*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2  
2800         HF=HD*HC-HB**2  
2801         HG=HD*HC-HA*HB  
2802         HH=(SQRT(HG**2+HE*HF)-HG)/(2.*HF)   
2803         DO 470 J=1,3    
2804         PCOR=HH*(P(N+1,J)-P(N+2,J)) 
2805         P(N+1,J)=P(N+1,J)+PCOR  
2806   470   P(N+2,J)=P(N+2,J)-PCOR  
2807         P(N+1,4)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2+P(N+1,5)**2)  
2808         P(N+2,4)=SQRT(P(N+2,1)**2+P(N+2,2)**2+P(N+2,3)**2+P(N+2,5)**2)  
2809         ND=ND-1 
2810       ENDIF 
2811     
2812 C...Check invariant mass of W jets. May give one particle or start over.    
2813   480 IF(MMAT.GE.42.AND.MMAT.LE.44.AND.IABS(K(N+1,2)).LT.10) THEN   
2814         PMR=SQRT(MAX(0.,P(N+1,5)**2+P(N+2,5)**2+2.*FOUR(N+1,N+2)))  
2815         MSTJ(93)=1  
2816         PM1=ULMASS(K(N+1,2))    
2817         MSTJ(93)=1  
2818         PM2=ULMASS(K(N+2,2))    
2819         IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 490    
2820         KFLDUM=INT(1.5+RLU(0))  
2821         CALL LUKFDI(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1)    
2822         CALL LUKFDI(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2)    
2823         IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 150   
2824         PSM=ULMASS(KF1)+ULMASS(KF2) 
2825         IF(MMAT.EQ.42.AND.PMR.GT.PARJ(64)+PSM) GOTO 490 
2826         IF(MMAT.GE.43.AND.PMR.GT.0.2*PARJ(32)+PSM) GOTO 490 
2827         IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 150   
2828         K(N+1,1)=1  
2829         KFTEMP=K(N+1,2) 
2830         CALL LUKFDI(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2))    
2831         IF(K(N+1,2).EQ.0) GOTO 150  
2832         P(N+1,5)=ULMASS(K(N+1,2))   
2833         K(N+2,2)=K(N+3,2)   
2834         P(N+2,5)=P(N+3,5)   
2835         PS=P(N+1,5)+P(N+2,5)    
2836         PV(2,5)=P(N+3,5)    
2837         MMAT=0  
2838         ND=2    
2839         GOTO 370    
2840       ENDIF 
2841     
2842 C...Phase space decay of partons from W decay.  
2843   490 IF(MMAT.EQ.42.AND.IABS(K(N+1,2)).LT.10) THEN  
2844         KFLO(1)=K(N+1,2)    
2845         KFLO(2)=K(N+2,2)    
2846         K(N+1,1)=K(N+3,1)   
2847         K(N+1,2)=K(N+3,2)   
2848         DO 500 J=1,5    
2849         PV(1,J)=P(N+1,J)+P(N+2,J)   
2850   500   P(N+1,J)=P(N+3,J)   
2851         PV(1,5)=PMR 
2852         N=N+1   
2853         NP=0    
2854         NQ=2    
2855         PS=0.   
2856         MSTJ(93)=2  
2857         PSQ=ULMASS(KFLO(1)) 
2858         MSTJ(93)=2  
2859         PSQ=PSQ+ULMASS(KFLO(2)) 
2860         MMAT=11 
2861         GOTO 180    
2862       ENDIF 
2863     
2864 C...Boost back for rapidly moving particle. 
2865   510 N=N+ND    
2866       IF(MBST.EQ.1) THEN    
2867         DO 520 J=1,3    
2868   520   BE(J)=P(IP,J)/P(IP,4)   
2869         GA=P(IP,4)/P(IP,5)  
2870         DO 540 I=NSAV+1,N   
2871         BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)  
2872         DO 530 J=1,3    
2873   530   P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)  
2874   540   P(I,4)=GA*(P(I,4)+BEP)  
2875       ENDIF 
2876     
2877 C...Fill in position of decay vertex.   
2878       DO 560 I=NSAV+1,N 
2879       DO 550 J=1,4  
2880   550 V(I,J)=VDCY(J)    
2881   560 V(I,5)=0. 
2882     
2883 C...Set up for parton shower evolution from jets.   
2884       IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN    
2885         K(NSAV+1,1)=3   
2886         K(NSAV+2,1)=3   
2887         K(NSAV+3,1)=3   
2888         K(NSAV+1,4)=MSTU(5)*(NSAV+2)    
2889         K(NSAV+1,5)=MSTU(5)*(NSAV+3)    
2890         K(NSAV+2,4)=MSTU(5)*(NSAV+3)    
2891         K(NSAV+2,5)=MSTU(5)*(NSAV+1)    
2892         K(NSAV+3,4)=MSTU(5)*(NSAV+1)    
2893         K(NSAV+3,5)=MSTU(5)*(NSAV+2)    
2894         MSTJ(92)=-(NSAV+1)  
2895       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN  
2896         K(NSAV+2,1)=3   
2897         K(NSAV+3,1)=3   
2898         K(NSAV+2,4)=MSTU(5)*(NSAV+3)    
2899         K(NSAV+2,5)=MSTU(5)*(NSAV+3)    
2900         K(NSAV+3,4)=MSTU(5)*(NSAV+2)    
2901         K(NSAV+3,5)=MSTU(5)*(NSAV+2)    
2902         MSTJ(92)=NSAV+2 
2903       ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46).    
2904      &AND.IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN 
2905         K(NSAV+1,1)=3   
2906         K(NSAV+2,1)=3   
2907         K(NSAV+1,4)=MSTU(5)*(NSAV+2)    
2908         K(NSAV+1,5)=MSTU(5)*(NSAV+2)    
2909         K(NSAV+2,4)=MSTU(5)*(NSAV+1)    
2910         K(NSAV+2,5)=MSTU(5)*(NSAV+1)    
2911         MSTJ(92)=NSAV+1 
2912       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21)  
2913      &THEN  
2914         K(NSAV+1,1)=3   
2915         K(NSAV+2,1)=3   
2916         K(NSAV+3,1)=3   
2917         KCP=LUCOMP(K(NSAV+1,2)) 
2918         KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2))    
2919         JCON=4  
2920         IF(KQP.LT.0) JCON=5 
2921         K(NSAV+1,JCON)=MSTU(5)*(NSAV+2) 
2922         K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1)   
2923         K(NSAV+2,JCON)=MSTU(5)*(NSAV+3) 
2924         K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2)   
2925         MSTJ(92)=NSAV+1 
2926       ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN 
2927         K(NSAV+1,1)=3   
2928         K(NSAV+3,1)=3   
2929         K(NSAV+1,4)=MSTU(5)*(NSAV+3)    
2930         K(NSAV+1,5)=MSTU(5)*(NSAV+3)    
2931         K(NSAV+3,4)=MSTU(5)*(NSAV+1)    
2932         K(NSAV+3,5)=MSTU(5)*(NSAV+1)    
2933         MSTJ(92)=NSAV+1 
2934       ENDIF 
2935     
2936 C...Mark decayed particle.  
2937       IF(K(IP,1).EQ.5) K(IP,1)=15   
2938       IF(K(IP,1).LE.10) K(IP,1)=11  
2939       K(IP,4)=NSAV+1    
2940       K(IP,5)=N 
2941     
2942       RETURN    
2943       END   
2944     
2945 C*********************************************************************  
2946     
2947       SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF)  
2948     
2949 C...Purpose: to generate a new flavour pair and combine off a hadron.   
2950       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
2951       SAVE /LUDAT1A/ 
2952       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
2953       SAVE /LUDAT2A/ 
2954     
2955       par3m=0.
2956       par4m=0.
2957       pardm=0.
2958       pars0=0.
2959       pars1=0.
2960       pars2=0.
2961       parsm=0.
2962       kmul=0
2963       ktab3=0
2964
2965 C...Default flavour values. Input consistency checks.   
2966       KF1A=IABS(KFL1)   
2967       KF2A=IABS(KFL2)   
2968       KFL3=0    
2969       KF=0  
2970       IF(KF1A.EQ.0) RETURN  
2971       IF(KF2A.NE.0) THEN    
2972         IF(KF1A.LE.10.AND.KF2A.LE.10.AND.KFL1*KFL2.GT.0) RETURN 
2973         IF(KF1A.GT.10.AND.KF2A.GT.10) RETURN    
2974         IF((KF1A.GT.10.OR.KF2A.GT.10).AND.KFL1*KFL2.LT.0) RETURN    
2975       ENDIF 
2976     
2977 C...Check if tabulated flavour probabilities are to be used.    
2978       IF(MSTJ(15).EQ.1) THEN    
2979         KTAB1=-1    
2980         IF(KF1A.GE.1.AND.KF1A.LE.6) KTAB1=KF1A  
2981         KFL1A=MOD(KF1A/1000,10) 
2982         KFL1B=MOD(KF1A/100,10)  
2983         KFL1S=MOD(KF1A,10)  
2984         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1B.GE.1.AND.KFL1B.LE.4) 
2985      &  KTAB1=6+KFL1A*(KFL1A-2)+2*KFL1B+(KFL1S-1)/2 
2986         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1A.EQ.KFL1B) KTAB1=KTAB1-1  
2987         IF(KF1A.GE.1.AND.KF1A.LE.6) KFL1A=KF1A  
2988         KTAB2=0 
2989         IF(KF2A.NE.0) THEN  
2990           KTAB2=-1  
2991           IF(KF2A.GE.1.AND.KF2A.LE.6) KTAB2=KF2A    
2992           KFL2A=MOD(KF2A/1000,10)   
2993           KFL2B=MOD(KF2A/100,10)    
2994           KFL2S=MOD(KF2A,10)    
2995           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2B.GE.1.AND.KFL2B.LE.4)   
2996      &    KTAB2=6+KFL2A*(KFL2A-2)+2*KFL2B+(KFL2S-1)/2   
2997           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2A.EQ.KFL2B) KTAB2=KTAB2-1    
2998         ENDIF   
2999         IF(KTAB1.GE.0.AND.KTAB2.GE.0) GOTO 140  
3000       ENDIF 
3001     
3002 C...Parameters and breaking diquark parameter combinations. 
3003   100 PAR2=PARJ(2)  
3004       PAR3=PARJ(3)  
3005       PAR4=3.*PARJ(4)   
3006       IF(MSTJ(12).GE.2) THEN    
3007         PAR3M=SQRT(PARJ(3)) 
3008         PAR4M=1./(3.*SQRT(PARJ(4))) 
3009         PARDM=PARJ(7)/(PARJ(7)+PAR3M*PARJ(6))   
3010         PARS0=PARJ(5)*(2.+(1.+PAR2*PAR3M*PARJ(7))*(1.+PAR4M))   
3011         PARS1=PARJ(7)*PARS0/(2.*PAR3M)+PARJ(5)*(PARJ(6)*(1.+PAR4M)+ 
3012      &  PAR2*PAR3M*PARJ(6)*PARJ(7)) 
3013         PARS2=PARJ(5)*2.*PARJ(6)*PARJ(7)*(PAR2*PARJ(7)+(1.+PAR4M)/PAR3M)    
3014         PARSM=MAX(PARS0,PARS1,PARS2)    
3015         PAR4=PAR4*(1.+PARSM)/(1.+PARSM/(3.*PAR4M))  
3016       ENDIF 
3017     
3018 C...Choice of whether to generate meson or baryon.  
3019       MBARY=0   
3020       KFDA=0    
3021       IF(KF1A.LE.10) THEN   
3022         IF(KF2A.EQ.0.AND.MSTJ(12).GE.1.AND.(1.+PARJ(1))*RLU(0).GT.1.)   
3023      &  MBARY=1 
3024         IF(KF2A.GT.10) MBARY=2  
3025         IF(KF2A.GT.10.AND.KF2A.LE.10000) KFDA=KF2A  
3026       ELSE  
3027         MBARY=2 
3028         IF(KF1A.LE.10000) KFDA=KF1A 
3029       ENDIF 
3030     
3031 C...Possibility of process diquark -> meson + new diquark.  
3032       IF(KFDA.NE.0.AND.MSTJ(12).GE.2) THEN  
3033         KFLDA=MOD(KFDA/1000,10) 
3034         KFLDB=MOD(KFDA/100,10)  
3035         KFLDS=MOD(KFDA,10)  
3036         WTDQ=PARS0  
3037         IF(MAX(KFLDA,KFLDB).EQ.3) WTDQ=PARS1    
3038         IF(MIN(KFLDA,KFLDB).EQ.3) WTDQ=PARS2    
3039         IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M) 
3040         IF((1.+WTDQ)*RLU(0).GT.1.) MBARY=-1 
3041         IF(MBARY.EQ.-1.AND.KF2A.NE.0) RETURN    
3042       ENDIF 
3043     
3044 C...Flavour for meson, possibly with new flavour.   
3045       IF(MBARY.LE.0) THEN   
3046         KFS=ISIGN(1,KFL1)   
3047         IF(MBARY.EQ.0) THEN 
3048           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),-KFL1)   
3049           KFLA=MAX(KF1A,KF2A+IABS(KFL3))    
3050           KFLB=MIN(KF1A,KF2A+IABS(KFL3))    
3051           IF(KFLA.NE.KF1A) KFS=-KFS 
3052     
3053 C...Splitting of diquark into meson plus new diquark.   
3054         ELSE    
3055           KFL1A=MOD(KF1A/1000,10)   
3056           KFL1B=MOD(KF1A/100,10)    
3057   110     KFL1D=KFL1A+INT(RLU(0)+0.5)*(KFL1B-KFL1A) 
3058           KFL1E=KFL1A+KFL1B-KFL1D   
3059           IF((KFL1D.EQ.3.AND.RLU(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.   
3060      &    RLU(0).LT.PARDM)) THEN    
3061             KFL1D=KFL1A+KFL1B-KFL1D 
3062             KFL1E=KFL1A+KFL1B-KFL1E 
3063           ENDIF 
3064           KFL3A=1+INT((2.+PAR2*PAR3M*PARJ(7))*RLU(0))   
3065           IF((KFL1E.NE.KFL3A.AND.RLU(0).GT.(1.+PAR4M)/MAX(2.,1.+PAR4M)).    
3066      &    OR.(KFL1E.EQ.KFL3A.AND.RLU(0).GT.2./MAX(2.,1.+PAR4M)))    
3067      &    GOTO 110  
3068           KFLDS=3   
3069           IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU(0)+1./(1.+PAR4M))+1    
3070           KFL3=ISIGN(10000+1000*MAX(KFL1E,KFL3A)+100*MIN(KFL1E,KFL3A)+  
3071      &    KFLDS,-KFL1)  
3072           KFLA=MAX(KFL1D,KFL3A) 
3073           KFLB=MIN(KFL1D,KFL3A) 
3074           IF(KFLA.NE.KFL1D) KFS=-KFS    
3075         ENDIF   
3076     
3077 C...Form meson, with spin and flavour mixing for diagonal states.   
3078         IF(KFLA.LE.2) KMUL=INT(PARJ(11)+RLU(0)) 
3079         IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU(0)) 
3080         IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU(0)) 
3081         IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN   
3082           IF(RLU(0).LT.PARJ(14)) KMUL=2 
3083         ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN 
3084           RMUL=RLU(0)   
3085           IF(RMUL.LT.PARJ(15)) KMUL=3   
3086           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)) KMUL=4    
3087           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)+PARJ(17)) KMUL=5   
3088         ENDIF   
3089         KFLS=3  
3090         IF(KMUL.EQ.0.OR.KMUL.EQ.3) KFLS=1   
3091         IF(KMUL.EQ.5) KFLS=5    
3092         IF(KFLA.NE.KFLB) THEN   
3093           KF=(100*KFLA+10*KFLB+KFLS)*KFS*(-1)**KFLA 
3094         ELSE    
3095           RMIX=RLU(0)   
3096           IMIX=2*KFLA+10*KMUL   
3097           IF(KFLA.LE.3) KF=110*(1+INT(RMIX+PARF(IMIX-1))+   
3098      &    INT(RMIX+PARF(IMIX)))+KFLS    
3099           IF(KFLA.GE.4) KF=110*KFLA+KFLS    
3100         ENDIF   
3101         IF(KMUL.EQ.2.OR.KMUL.EQ.3) KF=KF+ISIGN(10000,KF)    
3102         IF(KMUL.EQ.4) KF=KF+ISIGN(20000,KF) 
3103     
3104 C...Generate diquark flavour.   
3105       ELSE  
3106   120   IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN   
3107           KFLA=KF1A 
3108   130     KFLB=1+INT((2.+PAR2*PAR3)*RLU(0)) 
3109           KFLC=1+INT((2.+PAR2*PAR3)*RLU(0)) 
3110           KFLDS=1   
3111           IF(KFLB.GE.KFLC) KFLDS=3  
3112           IF(KFLDS.EQ.1.AND.PAR4*RLU(0).GT.1.) GOTO 130 
3113           IF(KFLDS.EQ.3.AND.PAR4.LT.RLU(0)) GOTO 130    
3114           KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1) 
3115     
3116 C...Take diquark flavour from input.    
3117         ELSEIF(KF1A.LE.10) THEN 
3118           KFLA=KF1A 
3119           KFLB=MOD(KF2A/1000,10)    
3120           KFLC=MOD(KF2A/100,10) 
3121           KFLDS=MOD(KF2A,10)    
3122     
3123 C...Generate (or take from input) quark to go with diquark. 
3124         ELSE    
3125           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),KFL1)    
3126           KFLA=KF2A+IABS(KFL3)  
3127           KFLB=MOD(KF1A/1000,10)    
3128           KFLC=MOD(KF1A/100,10) 
3129           KFLDS=MOD(KF1A,10)    
3130         ENDIF   
3131     
3132 C...SU(6) factors for formation of baryon. Try again if fails.  
3133         KBARY=KFLDS 
3134         IF(KFLDS.EQ.3.AND.KFLB.NE.KFLC) KBARY=5 
3135         IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC) KBARY=KBARY+1 
3136         WT=PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY)   
3137         IF(MBARY.EQ.1.AND.MSTJ(12).GE.2) THEN   
3138           WTDQ=PARS0    
3139           IF(MAX(KFLB,KFLC).EQ.3) WTDQ=PARS1    
3140           IF(MIN(KFLB,KFLC).EQ.3) WTDQ=PARS2    
3141           IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)   
3142           IF(KFLDS.EQ.1) WT=WT*(1.+WTDQ)/(1.+PARSM/(3.*PAR4M))  
3143           IF(KFLDS.EQ.3) WT=WT*(1.+WTDQ)/(1.+PARSM) 
3144         ENDIF   
3145         IF(KF2A.EQ.0.AND.WT.LT.RLU(0)) GOTO 120 
3146     
3147 C...Form baryon. Distinguish Lambda- and Sigmalike baryons. 
3148         KFLD=MAX(KFLA,KFLB,KFLC)    
3149         KFLF=MIN(KFLA,KFLB,KFLC)    
3150         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF   
3151         KFLS=2  
3152         IF((PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY))*RLU(0).GT.  
3153      &  PARF(60+KBARY)) KFLS=4  
3154         KFLL=0  
3155         IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF) THEN    
3156           IF(KFLDS.EQ.1.AND.KFLA.EQ.KFLD) KFLL=1    
3157           IF(KFLDS.EQ.1.AND.KFLA.NE.KFLD) KFLL=INT(0.25+RLU(0)) 
3158           IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU(0)) 
3159         ENDIF   
3160         IF(KFLL.EQ.0) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)    
3161         IF(KFLL.EQ.1) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)    
3162       ENDIF 
3163       RETURN    
3164     
3165 C...Use tabulated probabilities to select new flavour and hadron.   
3166   140 IF(KTAB2.EQ.0.AND.MSTJ(12).LE.0) THEN 
3167         KT3L=1  
3168         KT3U=6  
3169       ELSEIF(KTAB2.EQ.0.AND.KTAB1.GE.7.AND.MSTJ(12).LE.1) THEN  
3170         KT3L=1  
3171         KT3U=6  
3172       ELSEIF(KTAB2.EQ.0) THEN   
3173         KT3L=1  
3174         KT3U=22 
3175       ELSE  
3176         KT3L=KTAB2  
3177         KT3U=KTAB2  
3178       ENDIF 
3179       RFL=0.    
3180       DO 150 KTS=0,2    
3181       DO 150 KT3=KT3L,KT3U  
3182       RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3) 
3183   150 CONTINUE  
3184       RFL=RLU(0)*RFL    
3185       DO 160 KTS=0,2    
3186       KTABS=KTS 
3187       DO 160 KT3=KT3L,KT3U  
3188       KTAB3=KT3 
3189       RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3) 
3190   160 IF(RFL.LE.0.) GOTO 170    
3191   170 CONTINUE  
3192     
3193 C...Reconstruct flavour of produced quark/diquark.  
3194       IF(KTAB3.LE.6) THEN   
3195         KFL3A=KTAB3 
3196         KFL3B=0 
3197         KFL3=ISIGN(KFL3A,KFL1*(2*KTAB1-13)) 
3198       ELSE  
3199         KFL3A=1 
3200         IF(KTAB3.GE.8) KFL3A=2  
3201         IF(KTAB3.GE.11) KFL3A=3 
3202         IF(KTAB3.GE.16) KFL3A=4 
3203         KFL3B=(KTAB3-6-KFL3A*(KFL3A-2))/2   
3204         KFL3=1000*KFL3A+100*KFL3B+1 
3205         IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A*(KFL3A-2)+2*KFL3B) KFL3=  
3206      &  KFL3+2  
3207         KFL3=ISIGN(KFL3,KFL1*(13-2*KTAB1))  
3208       ENDIF 
3209     
3210 C...Reconstruct meson code. 
3211       IF(KFL3A.EQ.KFL1A.AND.KFL3B.EQ.KFL1B.AND.(KFL3A.LE.3.OR.  
3212      &KFL3B.NE.0)) THEN 
3213         RFL=RLU(0)*(PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+  
3214      &  25*KTABS)+PARF(145+80*KTAB1+25*KTABS))  
3215         KF=110+2*KTABS+1    
3216         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)) KF=220+2*KTABS+1 
3217         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+    
3218      &  25*KTABS)) KF=330+2*KTABS+1 
3219       ELSEIF(KTAB1.LE.6.AND.KTAB3.LE.6) THEN    
3220         KFLA=MAX(KTAB1,KTAB3)   
3221         KFLB=MIN(KTAB1,KTAB3)   
3222         KFS=ISIGN(1,KFL1)   
3223         IF(KFLA.NE.KF1A) KFS=-KFS   
3224         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA  
3225       ELSEIF(KTAB1.GE.7.AND.KTAB3.GE.7) THEN    
3226         KFS=ISIGN(1,KFL1)   
3227         IF(KFL1A.EQ.KFL3A) THEN 
3228           KFLA=MAX(KFL1B,KFL3B) 
3229           KFLB=MIN(KFL1B,KFL3B) 
3230           IF(KFLA.NE.KFL1B) KFS=-KFS    
3231         ELSEIF(KFL1A.EQ.KFL3B) THEN 
3232           KFLA=KFL3A    
3233           KFLB=KFL1B    
3234           KFS=-KFS  
3235         ELSEIF(KFL1B.EQ.KFL3A) THEN 
3236           KFLA=KFL1A    
3237           KFLB=KFL3B    
3238         ELSEIF(KFL1B.EQ.KFL3B) THEN 
3239           KFLA=MAX(KFL1A,KFL3A) 
3240           KFLB=MIN(KFL1A,KFL3A) 
3241           IF(KFLA.NE.KFL1A) KFS=-KFS    
3242         ELSE    
3243           CALL LUERRM(2,'(LUKFDI:) no matching flavours for qq -> qq')  
3244           GOTO 100  
3245         ENDIF   
3246         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA  
3247     
3248 C...Reconstruct baryon code.    
3249       ELSE  
3250         IF(KTAB1.GE.7) THEN 
3251           KFLA=KFL3A    
3252           KFLB=KFL1A    
3253           KFLC=KFL1B    
3254         ELSE    
3255           KFLA=KFL1A    
3256           KFLB=KFL3A    
3257           KFLC=KFL3B    
3258         ENDIF   
3259         KFLD=MAX(KFLA,KFLB,KFLC)    
3260         KFLF=MIN(KFLA,KFLB,KFLC)    
3261         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF   
3262         IF(KTABS.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+2,KFL1)  
3263         IF(KTABS.GE.1) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+2*KTABS,KFL1)    
3264       ENDIF 
3265     
3266 C...Check that constructed flavour code is an allowed one.  
3267       IF(KFL2.NE.0) KFL3=0  
3268       KC=LUCOMP(KF) 
3269       IF(KC.EQ.0) THEN  
3270         CALL LUERRM(2,'(LUKFDI:) user-defined flavour probabilities '// 
3271      &  'failed')   
3272         GOTO 100    
3273       ENDIF 
3274     
3275       RETURN    
3276       END   
3277     
3278 C*********************************************************************  
3279     
3280       SUBROUTINE LUPTDI(KFL,PX,PY)  
3281     
3282 C...Purpose: to generate transverse momentum according to a Gaussian.   
3283       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
3284       SAVE /LUDAT1A/ 
3285     
3286 C...Generate p_T and azimuthal angle, gives p_x and p_y.    
3287       KFLA=IABS(KFL)    
3288       PT=PARJ(21)*SQRT(-LOG(MAX(1E-10,RLU(0)))) 
3289       IF(MSTJ(91).EQ.1) PT=PARJ(22)*PT  
3290       IF(KFLA.EQ.0.AND.MSTJ(13).LE.0) PT=0. 
3291       PHI=PARU(2)*RLU(0)    
3292       PX=PT*COS(PHI)    
3293       PY=PT*SIN(PHI)    
3294     
3295       RETURN    
3296       END   
3297     
3298 C*********************************************************************  
3299     
3300       SUBROUTINE LUZDIS(KFL1,KFL2,PR,Z) 
3301     
3302 C...Purpose: to generate the longitudinal splitting variable z. 
3303       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
3304       SAVE /LUDAT1A/ 
3305     
3306       zdiv=0.
3307       fint=0.
3308       zdivc=0.
3309
3310 C...Check if heavy flavour fragmentation.   
3311       KFLA=IABS(KFL1)   
3312       KFLB=IABS(KFL2)   
3313       KFLH=KFLA 
3314       IF(KFLA.GE.10) KFLH=MOD(KFLA/1000,10) 
3315     
3316 C...Lund symmetric scaling function: determine parameters of shape. 
3317       IF(MSTJ(11).EQ.1.OR.(MSTJ(11).EQ.3.AND.KFLH.LE.3)) THEN   
3318         FA=PARJ(41) 
3319         IF(MSTJ(91).EQ.1) FA=PARJ(43)   
3320         IF(KFLB.GE.10) FA=FA+PARJ(45)   
3321         FB=PARJ(42)*PR  
3322         IF(MSTJ(91).EQ.1) FB=PARJ(44)*PR    
3323         FC=1.   
3324         IF(KFLA.GE.10) FC=FC-PARJ(45)   
3325         IF(KFLB.GE.10) FC=FC+PARJ(45)   
3326         MC=1    
3327         IF(ABS(FC-1.).GT.0.01) MC=2 
3328     
3329 C...Determine position of maximum. Special cases for a = 0 or a = c.    
3330         IF(FA.LT.0.02) THEN 
3331           MA=1  
3332           ZMAX=1.   
3333           IF(FC.GT.FB) ZMAX=FB/FC   
3334         ELSEIF(ABS(FC-FA).LT.0.01) THEN 
3335           MA=2  
3336           ZMAX=FB/(FB+FC)   
3337         ELSE    
3338           MA=3  
3339           ZMAX=0.5*(FB+FC-SQRT((FB-FC)**2+4.*FA*FB))/(FC-FA)    
3340           IF(ZMAX.GT.0.99.AND.FB.GT.100.) ZMAX=1.-FA/FB 
3341         ENDIF   
3342     
3343 C...Subdivide z range if distribution very peaked near endpoint.    
3344         MMAX=2  
3345         IF(ZMAX.LT.0.1) THEN    
3346           MMAX=1    
3347           ZDIV=2.75*ZMAX    
3348           IF(MC.EQ.1) THEN  
3349             FINT=1.-LOG(ZDIV)   
3350           ELSE  
3351             ZDIVC=ZDIV**(1.-FC) 
3352             FINT=1.+(1.-1./ZDIVC)/(FC-1.)   
3353           ENDIF 
3354         ELSEIF(ZMAX.GT.0.85.AND.FB.GT.1.) THEN  
3355           MMAX=3    
3356           FSCB=SQRT(4.+(FC/FB)**2)  
3357           ZDIV=FSCB-1./ZMAX-(FC/FB)*LOG(ZMAX*0.5*(FSCB+FC/FB))  
3358           IF(MA.GE.2) ZDIV=ZDIV+(FA/FB)*LOG(1.-ZMAX)    
3359           ZDIV=MIN(ZMAX,MAX(0.,ZDIV))   
3360           FINT=1.+FB*(1.-ZDIV)  
3361         ENDIF   
3362     
3363 C...Choice of z, preweighted for peaks at low or high z.    
3364   100   Z=RLU(0)    
3365         FPRE=1. 
3366         IF(MMAX.EQ.1) THEN  
3367           IF(FINT*RLU(0).LE.1.) THEN    
3368             Z=ZDIV*Z    
3369           ELSEIF(MC.EQ.1) THEN  
3370             Z=ZDIV**Z   
3371             FPRE=ZDIV/Z 
3372           ELSE  
3373             Z=1./(ZDIVC+Z*(1.-ZDIVC))**(1./(1.-FC)) 
3374             FPRE=(ZDIV/Z)**FC   
3375           ENDIF 
3376         ELSEIF(MMAX.EQ.3) THEN  
3377           IF(FINT*RLU(0).LE.1.) THEN    
3378             Z=ZDIV+LOG(Z)/FB    
3379             FPRE=EXP(FB*(Z-ZDIV))   
3380           ELSE  
3381             Z=ZDIV+Z*(1.-ZDIV)  
3382           ENDIF 
3383         ENDIF   
3384     
3385 C...Weighting according to correct formula. 
3386         IF(Z.LE.FB/(50.+FB).OR.Z.GE.1.) GOTO 100    
3387         FVAL=(ZMAX/Z)**FC*EXP(FB*(1./ZMAX-1./Z))    
3388         IF(MA.GE.2) FVAL=((1.-Z)/(1.-ZMAX))**FA*FVAL    
3389         IF(FVAL.LT.RLU(0)*FPRE) GOTO 100    
3390     
3391 C...Generate z according to Field-Feynman, SLAC, (1-z)**c OR z**c.  
3392       ELSE  
3393         FC=PARJ(50+MAX(1,KFLH)) 
3394         IF(MSTJ(91).EQ.1) FC=PARJ(59)   
3395   110   Z=RLU(0)    
3396         IF(FC.GE.0..AND.FC.LE.1.) THEN  
3397           IF(FC.GT.RLU(0)) Z=1.-Z**(1./3.)  
3398         ELSEIF(FC.GT.-1.) THEN  
3399           IF(-4.*FC*Z*(1.-Z)**2.LT.RLU(0)*((1.-Z)**2-FC*Z)**2) GOTO 110 
3400         ELSE    
3401           IF(FC.GT.0.) Z=1.-Z**(1./FC)  
3402           IF(FC.LT.0.) Z=Z**(-1./FC)    
3403         ENDIF   
3404       ENDIF 
3405     
3406       RETURN    
3407       END   
3408     
3409 C*********************************************************************  
3410     
3411       SUBROUTINE LUSHOW(IP1,IP2,QMAX)   
3412     
3413 C...Purpose: to generate timelike parton showers from given partons.    
3414       IMPLICIT DOUBLE PRECISION(D)  
3415       COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
3416       SAVE /LUJETSA/ 
3417       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
3418       SAVE /LUDAT1A/ 
3419       COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)    
3420       SAVE /LUDAT2A/ 
3421       DIMENSION PMTH(5,40),PS(5),PMA(4),PMSD(4),IEP(4),IPA(4),  
3422      &KFLA(4),KFLD(4),KFL(4),ITRY(4),ISI(4),ISL(4),DP(4),DPT(5,4)   
3423
3424       npa=0
3425       kflm=0
3426       pem=0.
3427       pmed=0.
3428       fbre=0.
3429       pm2=0.
3430       ped=0.
3431       zm=0.
3432       pa1s=0.
3433       pa2s=0.
3434       pa3s=0.
3435       pts=0.
3436       pzm=0.
3437       pmls=0.
3438       pt=0.
3439       hazip=0.
3440       hazic=0.
3441     
3442 C...Initialization of cutoff masses etc.    
3443       IF(MSTJ(41).LE.0.OR.(MSTJ(41).EQ.1.AND.QMAX.LE.PARJ(82)).OR.  
3444      &QMAX.LE.MIN(PARJ(82),PARJ(83)).OR.MSTJ(41).GE.3) RETURN   
3445       PMTH(1,21)=ULMASS(21) 
3446       PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25*PARJ(82)**2)   
3447       PMTH(3,21)=2.*PMTH(2,21)  
3448       PMTH(4,21)=PMTH(3,21) 
3449       PMTH(5,21)=PMTH(3,21) 
3450       PMTH(1,22)=ULMASS(22) 
3451       PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25*PARJ(83)**2)   
3452       PMTH(3,22)=2.*PMTH(2,22)  
3453       PMTH(4,22)=PMTH(3,22) 
3454       PMTH(5,22)=PMTH(3,22) 
3455       PMQTH1=PARJ(82)   
3456       IF(MSTJ(41).EQ.2) PMQTH1=MIN(PARJ(82),PARJ(83))   
3457       PMQTH2=PMTH(2,21) 
3458       IF(MSTJ(41).EQ.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))   
3459       DO 100 IF=1,8 
3460       PMTH(1,IF)=ULMASS(IF) 
3461       PMTH(2,IF)=SQRT(PMTH(1,IF)**2+0.25*PMQTH1**2) 
3462       PMTH(3,IF)=PMTH(2,IF)+PMQTH2  
3463       PMTH(4,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(82)**2)+PMTH(2,21)    
3464   100 PMTH(5,IF)=SQRT(PMTH(1,IF)**2+0.25*PARJ(83)**2)+PMTH(2,22)    
3465       PT2MIN=MAX(0.5*PARJ(82),1.1*PARJ(81))**2  
3466       ALAMS=PARJ(81)**2 
3467       ALFM=LOG(PT2MIN/ALAMS)    
3468     
3469 C...Store positions of shower initiating partons.   
3470       M3JC=0    
3471       IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN 
3472         NPA=1   
3473         IPA(1)=IP1  
3474       ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-   
3475      &MSTU(32))) THEN   
3476         NPA=2   
3477         IPA(1)=IP1  
3478         IPA(2)=IP2  
3479       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0.  
3480      &AND.IP2.GE.-3) THEN   
3481         NPA=IABS(IP2)   
3482         DO 110 I=1,NPA  
3483   110   IPA(I)=IP1+I-1  
3484       ELSE  
3485         CALL LUERRM(12, 
3486      &  '(LUSHOW:) failed to reconstruct showering system') 
3487         IF(MSTU(21).GE.1) RETURN    
3488       ENDIF 
3489     
3490 C...Check on phase space available for emission.    
3491       IREJ=0    
3492       DO 120 J=1,5  
3493   120 PS(J)=0.  
3494       PM=0. 
3495       DO 130 I=1,NPA    
3496       KFLA(I)=IABS(K(IPA(I),2)) 
3497       PMA(I)=P(IPA(I),5)    
3498       IF(KFLA(I).NE.0.AND.(KFLA(I).LE.8.OR.KFLA(I).EQ.21))  
3499      &PMA(I)=PMTH(3,KFLA(I))    
3500       PM=PM+PMA(I)  
3501       IF(KFLA(I).EQ.0.OR.(KFLA(I).GT.8.AND.KFLA(I).NE.21).OR.   
3502      &PMA(I).GT.QMAX) IREJ=IREJ+1   
3503       DO 130 J=1,4  
3504   130 PS(J)=PS(J)+P(IPA(I),J)   
3505       IF(IREJ.EQ.NPA) RETURN    
3506       PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))   
3507       IF(NPA.EQ.1) PS(5)=PS(4)  
3508       IF(PS(5).LE.PM+PMQTH1) RETURN 
3509       IF(NPA.EQ.2.AND.MSTJ(47).GE.1) THEN   
3510         IF(KFLA(1).GE.1.AND.KFLA(1).LE.8.AND.KFLA(2).GE.1.AND.  
3511      &  KFLA(2).LE.8) M3JC=1    
3512         IF(MSTJ(47).GE.2) M3JC=1    
3513       ENDIF 
3514     
3515 C...Define imagined single initiator of shower for parton system.   
3516       NS=N  
3517       IF(N.GT.MSTU(4)-MSTU(32)-5) THEN  
3518         CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETSA')   
3519         IF(MSTU(21).GE.1) RETURN    
3520       ENDIF 
3521       IF(NPA.GE.2) THEN 
3522         K(N+1,1)=11 
3523         K(N+1,2)=21 
3524         K(N+1,3)=0  
3525         K(N+1,4)=0  
3526         K(N+1,5)=0  
3527         P(N+1,1)=0. 
3528         P(N+1,2)=0. 
3529         P(N+1,3)=0. 
3530         P(N+1,4)=PS(5)  
3531         P(N+1,5)=PS(5)  
3532         V(N+1,5)=PS(5)**2   
3533         N=N+1   
3534       ENDIF 
3535     
3536 C...Loop over partons that may branch.  
3537       NEP=NPA   
3538       IM=NS 
3539       IF(NPA.EQ.1) IM=NS-1  
3540   140 IM=IM+1   
3541       IF(N.GT.NS) THEN  
3542         IF(IM.GT.N) GOTO 380    
3543         KFLM=IABS(K(IM,2))  
3544         IF(KFLM.EQ.0.OR.(KFLM.GT.8.AND.KFLM.NE.21)) GOTO 140    
3545         IF(P(IM,5).LT.PMTH(2,KFLM)) GOTO 140    
3546         IGM=K(IM,3) 
3547       ELSE  
3548         IGM=-1  
3549       ENDIF 
3550       IF(N+NEP.GT.MSTU(4)-MSTU(32)-5) THEN  
3551         CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETSA')   
3552         IF(MSTU(21).GE.1) RETURN    
3553       ENDIF 
3554     
3555 C...Position of aunt (sister to branching parton).  
3556 C...Origin and flavour of daughters.    
3557       IAU=0 
3558       IF(IGM.GT.0) THEN 
3559         IF(K(IM-1,3).EQ.IGM) IAU=IM-1   
3560         IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1 
3561       ENDIF 
3562       IF(IGM.GE.0) THEN 
3563         K(IM,4)=N+1 
3564         DO 150 I=1,NEP  
3565   150   K(N+I,3)=IM 
3566       ELSE  
3567         K(N+1,3)=IPA(1) 
3568       ENDIF 
3569       IF(IGM.LE.0) THEN 
3570         DO 160 I=1,NEP  
3571   160   K(N+I,2)=K(IPA(I),2)    
3572       ELSEIF(KFLM.NE.21) THEN   
3573         K(N+1,2)=K(IM,2)    
3574         K(N+2,2)=K(IM,5)    
3575       ELSEIF(K(IM,5).EQ.21) THEN    
3576         K(N+1,2)=21 
3577         K(N+2,2)=21 
3578       ELSE  
3579         K(N+1,2)=K(IM,5)    
3580         K(N+2,2)=-K(IM,5)   
3581       ENDIF 
3582     
3583 C...Reset flags on daughers and tries made. 
3584       DO 170 IP=1,NEP   
3585       K(N+IP,1)=3   
3586       K(N+IP,4)=0   
3587       K(N+IP,5)=0   
3588       KFLD(IP)=IABS(K(N+IP,2))  
3589       ITRY(IP)=0    
3590       ISL(IP)=0 
3591       ISI(IP)=0 
3592   170 IF(KFLD(IP).GT.0.AND.(KFLD(IP).LE.8.OR.KFLD(IP).EQ.21)) ISI(IP)=1 
3593       ISLM=0    
3594     
3595 C...Maximum virtuality of daughters.    
3596       IF(IGM.LE.0) THEN 
3597         DO 180 I=1,NPA  
3598         IF(NPA.GE.3) P(N+I,4)=(PS(4)*P(IPA(I),4)-PS(1)*P(IPA(I),1)- 
3599      &  PS(2)*P(IPA(I),2)-PS(3)*P(IPA(I),3))/PS(5)  
3600         P(N+I,5)=MIN(QMAX,PS(5))    
3601         IF(NPA.GE.3) P(N+I,5)=MIN(P(N+I,5),P(N+I,4))    
3602   180   IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)    
3603       ELSE  
3604         IF(MSTJ(43).LE.2) PEM=V(IM,2)   
3605         IF(MSTJ(43).GE.3) PEM=P(IM,4)   
3606         P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)   
3607         P(N+2,5)=MIN(P(IM,5),(1.-V(IM,1))*PEM)  
3608         IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)  
3609       ENDIF 
3610       DO 190 I=1,NEP    
3611       PMSD(I)=P(N+I,5)  
3612       IF(ISI(I).EQ.1) THEN  
3613         IF(P(N+I,5).LE.PMTH(3,KFLD(I))) P(N+I,5)=PMTH(1,KFLD(I))    
3614       ENDIF 
3615   190 V(N+I,5)=P(N+I,5)**2  
3616     
3617 C...Choose one of the daughters for evolution.  
3618   200 INUM=0    
3619       IF(NEP.EQ.1) INUM=1   
3620       DO 210 I=1,NEP    
3621   210 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I  
3622       DO 220 I=1,NEP    
3623       IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN   
3624         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) INUM=I  
3625       ENDIF 
3626   220 CONTINUE  
3627       IF(INUM.EQ.0) THEN    
3628         RMAX=0. 
3629         DO 230 I=1,NEP  
3630         IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQTH2) THEN  
3631           RPM=P(N+I,5)/PMSD(I)  
3632           IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,KFLD(I))) THEN  
3633             RMAX=RPM    
3634             INUM=I  
3635           ENDIF 
3636         ENDIF   
3637   230   CONTINUE    
3638       ENDIF 
3639     
3640 C...Store information on choice of evolving daughter.   
3641       INUM=MAX(1,INUM)  
3642       IEP(1)=N+INUM 
3643       DO 240 I=2,NEP    
3644       IEP(I)=IEP(I-1)+1 
3645   240 IF(IEP(I).GT.N+NEP) IEP(I)=N+1    
3646       DO 250 I=1,NEP    
3647   250 KFL(I)=IABS(K(IEP(I),2))  
3648       ITRY(INUM)=ITRY(INUM)+1   
3649       IF(ITRY(INUM).GT.200) THEN    
3650         CALL LUERRM(14,'(LUSHOW:) caught in infinite loop') 
3651         IF(MSTU(21).GE.1) RETURN    
3652       ENDIF 
3653       Z=0.5 
3654       IF(KFL(1).EQ.0.OR.(KFL(1).GT.8.AND.KFL(1).NE.21)) GOTO 300    
3655       IF(P(IEP(1),5).LT.PMTH(2,KFL(1))) GOTO 300    
3656     
3657 C...Calculate allowed z range.  
3658       IF(NEP.EQ.1) THEN 
3659         PMED=PS(4)  
3660       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN    
3661         PMED=P(IM,5)    
3662       ELSE  
3663         IF(INUM.EQ.1) PMED=V(IM,1)*PEM  
3664         IF(INUM.EQ.2) PMED=(1.-V(IM,1))*PEM 
3665       ENDIF 
3666       IF(MOD(MSTJ(43),2).EQ.1) THEN 
3667         ZC=PMTH(2,21)/PMED  
3668         ZCE=PMTH(2,22)/PMED 
3669       ELSE  
3670         ZC=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,21)/PMED)**2)))    
3671         IF(ZC.LT.1E-4) ZC=(PMTH(2,21)/PMED)**2  
3672         ZCE=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,22)/PMED)**2)))   
3673         IF(ZCE.LT.1E-4) ZCE=(PMTH(2,22)/PMED)**2    
3674       ENDIF 
3675       ZC=MIN(ZC,0.491)  
3676       ZCE=MIN(ZCE,0.491)    
3677       IF((MSTJ(41).EQ.1.AND.ZC.GT.0.49).OR.(MSTJ(41).EQ.2.AND.  
3678      &MIN(ZC,ZCE).GT.0.49)) THEN    
3679         P(IEP(1),5)=PMTH(1,KFL(1))  
3680         V(IEP(1),5)=P(IEP(1),5)**2  
3681         GOTO 300    
3682       ENDIF 
3683     
3684 C...Integral of Altarelli-Parisi z kernel for QCD.  
3685       IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN   
3686         FBR=6.*LOG((1.-ZC)/ZC)+MSTJ(45)*(0.5-ZC)    
3687       ELSEIF(MSTJ(49).EQ.0) THEN    
3688         FBR=(8./3.)*LOG((1.-ZC)/ZC) 
3689     
3690 C...Integral of Altarelli-Parisi z kernel for scalar gluon. 
3691       ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN   
3692         FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1.-2.*ZC) 
3693       ELSEIF(MSTJ(49).EQ.1) THEN    
3694         FBR=(1.-2.*ZC)/3.   
3695         IF(IGM.EQ.0.AND.M3JC.EQ.1) FBR=4.*FBR   
3696     
3697 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon. 
3698       ELSEIF(KFL(1).EQ.21) THEN 
3699         FBR=6.*MSTJ(45)*(0.5-ZC)    
3700       ELSE  
3701         FBR=2.*LOG((1.-ZC)/ZC)  
3702       ENDIF 
3703     
3704 C...Integral of Altarelli-Parisi kernel for photon emission.    
3705       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8) 
3706      &FBRE=(KCHG(KFL(1),1)/3.)**2*2.*LOG((1.-ZCE)/ZCE)  
3707     
3708 C...Inner veto algorithm starts. Find maximum mass for evolution.   
3709   260 PMS=V(IEP(1),5)   
3710       IF(IGM.GE.0) THEN 
3711         PM2=0.  
3712         DO 270 I=2,NEP  
3713         PM=P(IEP(I),5)  
3714         IF(KFL(I).GT.0.AND.(KFL(I).LE.8.OR.KFL(I).EQ.21)) PM=   
3715      &  PMTH(2,KFL(I))  
3716   270   PM2=PM2+PM  
3717         PMS=MIN(PMS,(P(IM,5)-PM2)**2)   
3718       ENDIF 
3719     
3720 C...Select mass for daughter in QCD evolution.  
3721       B0=27./6. 
3722       DO 280 IF=4,MSTJ(45)  
3723   280 IF(PMS.GT.4.*PMTH(2,IF)**2) B0=(33.-2.*IF)/6. 
3724       IF(MSTJ(44).LE.0) THEN    
3725         PMSQCD=PMS*EXP(MAX(-100.,LOG(RLU(0))*PARU(2)/(PARU(111)*FBR)))  
3726       ELSEIF(MSTJ(44).EQ.1) THEN    
3727         PMSQCD=4.*ALAMS*(0.25*PMS/ALAMS)**(RLU(0)**(B0/FBR))    
3728       ELSE  
3729         PMSQCD=PMS*RLU(0)**(ALFM*B0/FBR)    
3730       ENDIF 
3731       IF(ZC.GT.0.49.OR.PMSQCD.LE.PMTH(4,KFL(1))**2) PMSQCD= 
3732      &PMTH(2,KFL(1))**2 
3733       V(IEP(1),5)=PMSQCD    
3734       MCE=1 
3735     
3736 C...Select mass for daughter in QED evolution.  
3737       IF(MSTJ(41).EQ.2.AND.KFL(1).GE.1.AND.KFL(1).LE.8) THEN    
3738         PMSQED=PMS*EXP(MAX(-100.,LOG(RLU(0))*PARU(2)/(PARU(101)*FBRE))) 
3739         IF(ZCE.GT.0.49.OR.PMSQED.LE.PMTH(5,KFL(1))**2) PMSQED=  
3740      &  PMTH(2,KFL(1))**2   
3741         IF(PMSQED.GT.PMSQCD) THEN   
3742           V(IEP(1),5)=PMSQED    
3743           MCE=2 
3744         ENDIF   
3745       ENDIF 
3746     
3747 C...Check whether daughter mass below cutoff.   
3748       P(IEP(1),5)=SQRT(V(IEP(1),5)) 
3749       IF(P(IEP(1),5).LE.PMTH(3,KFL(1))) THEN    
3750         P(IEP(1),5)=PMTH(1,KFL(1))  
3751         V(IEP(1),5)=P(IEP(1),5)**2  
3752         GOTO 300    
3753       ENDIF 
3754     
3755 C...Select z value of branching: q -> qgamma.   
3756       IF(MCE.EQ.2) THEN 
3757         Z=1.-(1.-ZCE)*(ZCE/(1.-ZCE))**RLU(0)    
3758         IF(1.+Z**2.LT.2.*RLU(0)) GOTO 260   
3759         K(IEP(1),5)=22  
3760     
3761 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.  
3762       ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN   
3763         Z=1.-(1.-ZC)*(ZC/(1.-ZC))**RLU(0)   
3764         IF(1.+Z**2.LT.2.*RLU(0)) GOTO 260   
3765         K(IEP(1),5)=21  
3766       ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*(0.5-ZC).LT.RLU(0)*FBR) THEN    
3767         Z=(1.-ZC)*(ZC/(1.-ZC))**RLU(0)  
3768         IF(RLU(0).GT.0.5) Z=1.-Z    
3769         IF((1.-Z*(1.-Z))**2.LT.RLU(0)) GOTO 260 
3770         K(IEP(1),5)=21  
3771       ELSEIF(MSTJ(49).NE.1) THEN    
3772         Z=ZC+(1.-2.*ZC)*RLU(0)  
3773         IF(Z**2+(1.-Z)**2.LT.RLU(0)) GOTO 260   
3774         KFLB=1+INT(MSTJ(45)*RLU(0)) 
3775         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)  
3776         IF(PMQ.GE.1.) GOTO 260  
3777         PMQ0=4.*PMTH(2,21)**2/V(IEP(1),5)   
3778         IF(MOD(MSTJ(43),2).EQ.0.AND.(1.+0.5*PMQ)*SQRT(1.-PMQ).LT.   
3779      &  RLU(0)*(1.+0.5*PMQ0)*SQRT(1.-PMQ0)) GOTO 260    
3780         K(IEP(1),5)=KFLB    
3781     
3782 C...Ditto for scalar gluon model.   
3783       ELSEIF(KFL(1).NE.21) THEN 
3784         Z=1.-SQRT(ZC**2+RLU(0)*(1.-2.*ZC))  
3785         K(IEP(1),5)=21  
3786       ELSEIF(RLU(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN  
3787         Z=ZC+(1.-2.*ZC)*RLU(0)  
3788         K(IEP(1),5)=21  
3789       ELSE  
3790         Z=ZC+(1.-2.*ZC)*RLU(0)  
3791         KFLB=1+INT(MSTJ(45)*RLU(0)) 
3792         PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)  
3793         IF(PMQ.GE.1.) GOTO 260  
3794         K(IEP(1),5)=KFLB    
3795       ENDIF 
3796       IF(MCE.EQ.1.AND.MSTJ(44).GE.2) THEN   
3797         IF(Z*(1.-Z)*V(IEP(1),5).LT.PT2MIN) GOTO 260 
3798         IF(ALFM/LOG(V(IEP(1),5)*Z*(1.-Z)/ALAMS).LT.RLU(0)) GOTO 260 
3799       ENDIF 
3800     
3801 C...Check if z consistent with chosen m.    
3802       IF(KFL(1).EQ.21) THEN 
3803         KFLGD1=IABS(K(IEP(1),5))    
3804         KFLGD2=KFLGD1   
3805       ELSE  
3806         KFLGD1=KFL(1)   
3807         KFLGD2=IABS(K(IEP(1),5))    
3808       ENDIF 
3809       IF(NEP.EQ.1) THEN 
3810         PED=PS(4)   
3811       ELSEIF(NEP.GE.3) THEN 
3812         PED=P(IEP(1),4) 
3813       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN    
3814         PED=0.5*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)    
3815       ELSE  
3816         IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM   
3817         IF(IEP(1).EQ.N+2) PED=(1.-V(IM,1))*PEM  
3818       ENDIF 
3819       IF(MOD(MSTJ(43),2).EQ.1) THEN 
3820         PMQTH3=0.5*PARJ(82) 
3821         IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)    
3822         PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(IEP(1),5)  
3823         PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(IEP(1),5)  
3824         ZD=SQRT(MAX(0.,(1.-V(IEP(1),5)/PED**2)*((1.-PMQ1-PMQ2)**2-  
3825      &  4.*PMQ1*PMQ2))) 
3826         ZH=1.+PMQ1-PMQ2 
3827       ELSE  
3828         ZD=SQRT(MAX(0.,1.-V(IEP(1),5)/PED**2))  
3829         ZH=1.   
3830       ENDIF 
3831       ZL=0.5*(ZH-ZD)    
3832       ZU=0.5*(ZH+ZD)    
3833       IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 260   
3834       IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL* 
3835      &(1.-ZU))) 
3836       IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))    
3837     
3838 C...Three-jet matrix element correction.    
3839       IF(IGM.EQ.0.AND.M3JC.EQ.1) THEN   
3840         X1=Z*(1.+V(IEP(1),5)/V(NS+1,5)) 
3841         X2=1.-V(IEP(1),5)/V(NS+1,5) 
3842         X3=(1.-X1)+(1.-X2)  
3843         IF(MCE.EQ.2) THEN   
3844           KI1=K(IPA(INUM),2)    
3845           KI2=K(IPA(3-INUM),2)  
3846           QF1=KCHG(IABS(KI1),1)*ISIGN(1,KI1)/3. 
3847           QF2=KCHG(IABS(KI2),1)*ISIGN(1,KI2)/3. 
3848           WSHOW=QF1**2*(1.-X1)/X3*(1.+(X1/(2.-X2))**2)+ 
3849      &    QF2**2*(1.-X2)/X3*(1.+(X2/(2.-X1))**2)    
3850           WME=(QF1*(1.-X1)/X3-QF2*(1.-X2)/X3)**2*(X1**2+X2**2)  
3851         ELSEIF(MSTJ(49).NE.1) THEN  
3852           WSHOW=1.+(1.-X1)/X3*(X1/(2.-X2))**2+  
3853      &    (1.-X2)/X3*(X2/(2.-X1))**2    
3854           WME=X1**2+X2**2   
3855         ELSE    
3856           WSHOW=4.*X3*((1.-X1)/(2.-X2)**2+(1.-X2)/(2.-X1)**2)   
3857           WME=X3**2 
3858         ENDIF   
3859         IF(WME.LT.RLU(0)*WSHOW) GOTO 260    
3860     
3861 C...Impose angular ordering by rejection of nonordered emission.    
3862       ELSEIF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2) THEN  
3863         MAOM=1  
3864         ZM=V(IM,1)  
3865         IF(IEP(1).EQ.N+2) ZM=1.-V(IM,1) 
3866         THE2ID=Z*(1.-Z)*(ZM*P(IM,4))**2/V(IEP(1),5) 
3867         IAOM=IM 
3868   290   IF(K(IAOM,5).EQ.22) THEN    
3869           IAOM=K(IAOM,3)    
3870           IF(K(IAOM,3).LE.NS) MAOM=0    
3871           IF(MAOM.EQ.1) GOTO 290    
3872         ENDIF   
3873         IF(MAOM.EQ.1) THEN  
3874           THE2IM=V(IAOM,1)*(1.-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)    
3875           IF(THE2ID.LT.THE2IM) GOTO 260 
3876         ENDIF   
3877       ENDIF 
3878     
3879 C...Impose user-defined maximum angle at first branching.   
3880       IF(MSTJ(48).EQ.1) THEN    
3881         IF(NEP.EQ.1.AND.IM.EQ.NS) THEN  
3882           THE2ID=Z*(1.-Z)*PS(4)**2/V(IEP(1),5)  
3883           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260 
3884         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN    
3885           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)  
3886           IF(THE2ID.LT.1./PARJ(85)**2) GOTO 260 
3887         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN    
3888           THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)  
3889           IF(THE2ID.LT.1./PARJ(86)**2) GOTO 260 
3890         ENDIF   
3891       ENDIF 
3892     
3893 C...End of inner veto algorithm. Check if only one leg evolved so far.  
3894   300 V(IEP(1),1)=Z 
3895       ISL(1)=0  
3896       ISL(2)=0  
3897       IF(NEP.EQ.1) GOTO 330 
3898       IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 200  
3899       DO 310 I=1,NEP    
3900       IF(ITRY(I).EQ.0.AND.KFLD(I).GT.0.AND.(KFLD(I).LE.8.OR.KFLD(I).EQ. 
3901      &21)) THEN 
3902         IF(P(N+I,5).GE.PMTH(2,KFLD(I))) GOTO 200    
3903       ENDIF 
3904   310 CONTINUE  
3905     
3906 C...Check if chosen multiplet m1,m2,z1,z2 is physical.  
3907       IF(NEP.EQ.3) THEN 
3908         PA1S=(P(N+1,4)+P(N+1,5))*(P(N+1,4)-P(N+1,5))    
3909         PA2S=(P(N+2,4)+P(N+2,5))*(P(N+2,4)-P(N+2,5))    
3910         PA3S=(P(N+3,4)+P(N+3,5))*(P(N+3,4)-P(N+3,5))    
3911         PTS=0.25*(2.*PA1S*PA2S+2.*PA1S*PA3S+2.*PA2S*PA3S-   
3912      &  PA1S**2-PA2S**2-PA3S**2)/PA1S   
3913         IF(PTS.LE.0.) GOTO 200  
3914       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN    
3915         DO 320 I1=N+1,N+2   
3916         KFLDA=IABS(K(I1,2)) 
3917         IF(KFLDA.EQ.0.OR.(KFLDA.GT.8.AND.KFLDA.NE.21)) GOTO 320 
3918         IF(P(I1,5).LT.PMTH(2,KFLDA)) GOTO 320   
3919         IF(KFLDA.EQ.21) THEN    
3920           KFLGD1=IABS(K(I1,5))  
3921           KFLGD2=KFLGD1 
3922         ELSE    
3923           KFLGD1=KFLDA  
3924           KFLGD2=IABS(K(I1,5))  
3925         ENDIF   
3926         I2=2*N+3-I1 
3927         IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN  
3928           PED=0.5*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5) 
3929         ELSE    
3930           IF(I1.EQ.N+1) ZM=V(IM,1)  
3931           IF(I1.EQ.N+2) ZM=1.-V(IM,1)   
3932           PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-  
3933      &    4.*V(N+1,5)*V(N+2,5)) 
3934           PED=PEM*(0.5*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/V(IM,5)    
3935         ENDIF   
3936         IF(MOD(MSTJ(43),2).EQ.1) THEN   
3937           PMQTH3=0.5*PARJ(82)   
3938           IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)  
3939           PMQ1=(PMTH(1,KFLGD1)**2+PMQTH3**2)/V(I1,5)    
3940           PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(I1,5)    
3941           ZD=SQRT(MAX(0.,(1.-V(I1,5)/PED**2)*((1.-PMQ1-PMQ2)**2-    
3942      &    4.*PMQ1*PMQ2)))   
3943           ZH=1.+PMQ1-PMQ2   
3944         ELSE    
3945           ZD=SQRT(MAX(0.,1.-V(I1,5)/PED**2))    
3946           ZH=1. 
3947         ENDIF   
3948         ZL=0.5*(ZH-ZD)  
3949         ZU=0.5*(ZH+ZD)  
3950         IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(1)=1 
3951         IF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(2)=1 
3952         IF(KFLDA.EQ.21) V(I1,4)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*(1.-ZU)))   
3953         IF(KFLDA.NE.21) V(I1,4)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))   
3954   320   CONTINUE    
3955         IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN  
3956           ISL(3-ISLM)=0 
3957           ISLM=3-ISLM   
3958         ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN    
3959           ZDR1=MAX(0.,V(N+1,3)/V(N+1,4)-1.) 
3960           ZDR2=MAX(0.,V(N+2,3)/V(N+2,4)-1.) 
3961           IF(ZDR2.GT.RLU(0)*(ZDR1+ZDR2)) ISL(1)=0   
3962           IF(ISL(1).EQ.1) ISL(2)=0  
3963           IF(ISL(1).EQ.0) ISLM=1    
3964           IF(ISL(2).EQ.0) ISLM=2    
3965         ENDIF   
3966         IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 200 
3967       ENDIF 
3968       IF(IGM.GT.0.AND.MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.    
3969      &PMTH(2,KFLD(1)).OR.P(N+2,5).GE.PMTH(2,KFLD(2)))) THEN 
3970         PMQ1=V(N+1,5)/V(IM,5)   
3971         PMQ2=V(N+2,5)/V(IM,5)   
3972         ZD=SQRT(MAX(0.,(1.-V(IM,5)/PEM**2)*((1.-PMQ1-PMQ2)**2-  
3973      &  4.*PMQ1*PMQ2))) 
3974         ZH=1.+PMQ1-PMQ2 
3975         ZL=0.5*(ZH-ZD)  
3976         ZU=0.5*(ZH+ZD)  
3977         IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 200 
3978       ENDIF 
3979     
3980 C...Accepted branch. Construct four-momentum for initial partons.   
3981   330 MAZIP=0   
3982       MAZIC=0