Coding violation fixes
[u/mrichter/AliRoot.git] / THydjet / hydjet1_1 / jetset_73.f
1 C*********************************************************************
2 CCPH   This file has enlarged event record, LUJETS size=30000
3 C*********************************************************************
4 C*********************************************************************
5 C*********************************************************************
6 C*                                                                  **
7 C*                                                     June 1991    **
8 C*                                                                  **
9 C*   The Lund Monte Carlo for Jet Fragmentation and e+e- Physics    **
10 C*                                                                  **
11 C*                        JETSET version 7.3                        **
12 C*                                                                  **
13 C*                        Torbjorn Sjostrand                        **
14 C*                                                                  **
15 C*                    CERN/TH, CH-1211 Geneva 23                    **
16 C*                BITNET/EARN address TORSJO@CERNVM                 **
17 C*                       Tel. +22 - 767 28 20                       **
18 C*                                                                  **
19 C*          LUSHOW is written together with Mats Bengtsson          **
20 C*                                                                  **
21 C*           A complete manual exists on a separate file            **
22 C*         Please report any program errors to the author!          **
23 C*                                                                  **
24 C*                   Copyright Torbjorn Sjostrand                   **
25 C*                                                                  **
26 C*********************************************************************
27 C*********************************************************************
28 C                                                                    *
29 C  List of subprograms in order of appearance, with main purpose     *
30 C  (S = subroutine, F = function, B = block data)                    *
31 C                                                                    *
32 C  S   LU1ENT   to fill one entry (= parton or particle)             *
33 C  S   LU2ENT   to fill two entries                                  *
34 C  S   LU3ENT   to fill three entries                                *
35 C  S   LU4ENT   to fill four entries                                 *
36 C  S   LUJOIN   to connect entries with colour flow information      *
37 C  S   LUGIVE   to fill (or query) commonblock variables             *
38 C  S   LUEXEC   to administrate fragmentation and decay chain        *
39 C  S   LUPREP   to rearrange showered partons along strings          *
40 C  S   LUSTRF   to do string fragmentation of jet system             *
41 C  S   LUINDF   to do independent fragmentation of one or many jets  *
42 C  S   LUDECY   to do the decay of a particle                        *
43 C  S   LUKFDI   to select parton and hadron flavours in fragm        *
44 C  S   LUPTDI   to select transverse momenta in fragm                *
45 C  S   LUZDIS   to select longitudinal scaling variable in fragm     *
46 C  S   LUSHOW   to do timelike parton shower evolution               *
47 C  S   LUBOEI   to include Bose-Einstein effects (crudely)           *
48 C  F   ULMASS   to give the mass of a particle or parton             *
49 C  S   LUNAME   to give the name of a particle or parton             *
50 C  F   LUCHGE   to give three times the electric charge              *
51 C  F   LUCOMP   to compress standard KF flavour code to internal KC  *
52 C  S   LUERRM   to write error messages and abort faulty run         *
53 C  F   ULALEM   to give the alpha_electromagnetic value              *
54 C  F   ULALPS   to give the alpha_strong value                       *
55 C  F   ULANGL   to give the angle from known x and y components      *
56 C  F   RLU      to provide a random number generator                 *
57 C  S   RLUGET   to save the state of the random number generator     *
58 C  S   RLUSET   to set the state of the random number generator      *
59 C  S   LUROBO   to rotate and/or boost an event                      *
60 C  S   LUEDIT   to remove unwanted entries from record               *
61 C  S   LULIST   to list event record or particle data                *
62 C  S   LUUPDA   to update particle data                              *
63 C  F   KLU      to provide integer-valued event information          *
64 C  F   PLU      to provide real-valued event information             *
65 C  S   LUSPHE   to perform sphericity analysis                       *
66 C  S   LUTHRU   to perform thrust analysis                           *
67 C  S   LUCLUS   to perform three-dimensional cluster analysis        *
68 C  S   LUCELL   to perform cluster analysis in (eta, phi, E_T)       *
69 C  S   LUJMAS   to give high and low jet mass of event               *
70 C  S   LUFOWO   to give Fox-Wolfram moments                          *
71 C  S   LUTABU   to analyze events, with tabular output               *
72 C                                                                    *
73 C  S   LUEEVT   to administrate the generation of an e+e- event      *
74 C  S   LUXTOT   to give the total cross-section at given CM energy   *
75 C  S   LURADK   to generate initial state photon radiation           *
76 C  S   LUXKFL   to select flavour of primary qqbar pair              *
77 C  S   LUXJET   to select (matrix element) jet multiplicity          *
78 C  S   LUX3JT   to select kinematics of three-jet event              *
79 C  S   LUX4JT   to select kinematics of four-jet event               *
80 C  S   LUXDIF   to select angular orientation of event               *
81 C  S   LUONIA   to perform generation of onium decay to gluons       *
82 C                                                                    *
83 C  S   LUHEPC   to convert between /LUJETS/ and /HEPEVT/ records     *
84 C  S   LUTEST   to test the proper functioning of the package        *
85 C  B   LUDATA   to contain default values and particle data          *
86 C                                                                    *
87 C*********************************************************************
88
89       SUBROUTINE LU1ENT(IP,KF,PE,THE,PHI)
90
91 C...Purpose: to store one parton/particle in commonblock LUJETS.
92       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
93       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
94       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
95       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
96
97 C...Standard checks.
98       MSTU(28)=0
99       IF(MSTU(12).GE.1) CALL LULIST(0)
100       IPA=MAX(1,IABS(IP))
101       IF(IPA.GT.MSTU(4)) CALL LUERRM(21,
102      &'(LU1ENT:) writing outside LUJETS memory')
103       KC=LUCOMP(KF)
104       IF(KC.EQ.0) CALL LUERRM(12,'(LU1ENT:) unknown flavour code')
105
106 C...Find mass. Reset K, P and V vectors.
107       PM=0.
108       IF(MSTU(10).EQ.1) PM=P(IPA,5)
109       IF(MSTU(10).GE.2) PM=ULMASS(KF)
110       DO 100 J=1,5
111       K(IPA,J)=0
112       P(IPA,J)=0.
113   100 V(IPA,J)=0.
114
115 C...Store parton/particle in K and P vectors.
116       K(IPA,1)=1
117       IF(IP.LT.0) K(IPA,1)=2
118       K(IPA,2)=KF
119       P(IPA,5)=PM
120       P(IPA,4)=MAX(PE,PM)
121       PA=SQRT(P(IPA,4)**2-P(IPA,5)**2)
122       P(IPA,1)=PA*SIN(THE)*COS(PHI)
123       P(IPA,2)=PA*SIN(THE)*SIN(PHI)
124       P(IPA,3)=PA*COS(THE)
125
126 C...Set N. Optionally fragment/decay.
127       N=IPA
128       IF(IP.EQ.0) CALL LUEXEC
129
130       RETURN
131       END
132
133 C*********************************************************************
134
135       SUBROUTINE LU2ENT(IP,KF1,KF2,PECM)
136
137 C...Purpose: to store two partons/particles in their CM frame,
138 C...with the first along the +z axis.
139       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
140       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
141       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
142       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
143
144 C...Standard checks.
145       MSTU(28)=0
146       IF(MSTU(12).GE.1) CALL LULIST(0)
147       IPA=MAX(1,IABS(IP))
148       IF(IPA.GT.MSTU(4)-1) CALL LUERRM(21,
149      &'(LU2ENT:) writing outside LUJETS memory')
150       KC1=LUCOMP(KF1)
151       KC2=LUCOMP(KF2)
152       IF(KC1.EQ.0.OR.KC2.EQ.0) CALL LUERRM(12,
153      &'(LU2ENT:) unknown flavour code')
154
155 C...Find masses. Reset K, P and V vectors.
156       PM1=0.
157       IF(MSTU(10).EQ.1) PM1=P(IPA,5)
158       IF(MSTU(10).GE.2) PM1=ULMASS(KF1)
159       PM2=0.
160       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
161       IF(MSTU(10).GE.2) PM2=ULMASS(KF2)
162       DO 100 I=IPA,IPA+1
163       DO 100 J=1,5
164       K(I,J)=0
165       P(I,J)=0.
166   100 V(I,J)=0.
167
168 C...Check flavours.
169       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
170       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
171       IF(MSTU(19).EQ.1) THEN
172         MSTU(19)=0
173       ELSE
174         IF(KQ1+KQ2.NE.0.AND.KQ1+KQ2.NE.4) CALL LUERRM(2,
175      &  '(LU2ENT:) unphysical flavour combination')
176       ENDIF
177       K(IPA,2)=KF1
178       K(IPA+1,2)=KF2
179
180 C...Store partons/particles in K vectors for normal case.
181       IF(IP.GE.0) THEN
182         K(IPA,1)=1
183         IF(KQ1.NE.0.AND.KQ2.NE.0) K(IPA,1)=2
184         K(IPA+1,1)=1
185
186 C...Store partons in K vectors for parton shower evolution.
187       ELSE
188         K(IPA,1)=3
189         K(IPA+1,1)=3
190         K(IPA,4)=MSTU(5)*(IPA+1)
191         K(IPA,5)=K(IPA,4)
192         K(IPA+1,4)=MSTU(5)*IPA
193         K(IPA+1,5)=K(IPA+1,4)
194       ENDIF
195
196 C...Check kinematics and store partons/particles in P vectors.
197       IF(PECM.LE.PM1+PM2) CALL LUERRM(13,
198      &'(LU2ENT:) energy smaller than sum of masses')
199       PA=SQRT(MAX(0.,(PECM**2-PM1**2-PM2**2)**2-(2.*PM1*PM2)**2))/
200      &(2.*PECM)
201       P(IPA,3)=PA
202       P(IPA,4)=SQRT(PM1**2+PA**2)
203       P(IPA,5)=PM1
204       P(IPA+1,3)=-PA
205       P(IPA+1,4)=SQRT(PM2**2+PA**2)
206       P(IPA+1,5)=PM2
207
208 C...Set N. Optionally fragment/decay.
209       N=IPA+1
210       IF(IP.EQ.0) CALL LUEXEC
211
212       RETURN
213       END
214
215 C*********************************************************************
216
217       SUBROUTINE LU3ENT(IP,KF1,KF2,KF3,PECM,X1,X3)
218
219 C...Purpose: to store three partons or particles in their CM frame,
220 C...with the first along the +z axis and the third in the (x,z)
221 C...plane with x > 0.
222       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
223       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
224       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
225       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
226
227 C...Standard checks.
228       MSTU(28)=0
229       IF(MSTU(12).GE.1) CALL LULIST(0)
230       IPA=MAX(1,IABS(IP))
231       IF(IPA.GT.MSTU(4)-2) CALL LUERRM(21,
232      &'(LU3ENT:) writing outside LUJETS memory')
233       KC1=LUCOMP(KF1)
234       KC2=LUCOMP(KF2)
235       KC3=LUCOMP(KF3)
236       IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0) CALL LUERRM(12,
237      &'(LU3ENT:) unknown flavour code')
238
239 C...Find masses. Reset K, P and V vectors.
240       PM1=0.
241       IF(MSTU(10).EQ.1) PM1=P(IPA,5)
242       IF(MSTU(10).GE.2) PM1=ULMASS(KF1)
243       PM2=0.
244       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
245       IF(MSTU(10).GE.2) PM2=ULMASS(KF2)
246       PM3=0.
247       IF(MSTU(10).EQ.1) PM3=P(IPA+2,5)
248       IF(MSTU(10).GE.2) PM3=ULMASS(KF3)
249       DO 100 I=IPA,IPA+2
250       DO 100 J=1,5
251       K(I,J)=0
252       P(I,J)=0.
253   100 V(I,J)=0.
254
255 C...Check flavours.
256       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
257       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
258       KQ3=KCHG(KC3,2)*ISIGN(1,KF3)
259       IF(MSTU(19).EQ.1) THEN
260         MSTU(19)=0
261       ELSEIF(KQ1.EQ.0.AND.KQ2.EQ.0.AND.KQ3.EQ.0) THEN
262       ELSEIF(KQ1.NE.0.AND.KQ2.EQ.2.AND.(KQ1+KQ3.EQ.0.OR.
263      &KQ1+KQ3.EQ.4)) THEN
264       ELSE
265         CALL LUERRM(2,'(LU3ENT:) unphysical flavour combination')
266       ENDIF
267       K(IPA,2)=KF1
268       K(IPA+1,2)=KF2
269       K(IPA+2,2)=KF3
270
271 C...Store partons/particles in K vectors for normal case.
272       IF(IP.GE.0) THEN
273         K(IPA,1)=1
274         IF(KQ1.NE.0.AND.(KQ2.NE.0.OR.KQ3.NE.0)) K(IPA,1)=2
275         K(IPA+1,1)=1
276         IF(KQ2.NE.0.AND.KQ3.NE.0) K(IPA+1,1)=2
277         K(IPA+2,1)=1
278
279 C...Store partons in K vectors for parton shower evolution.
280       ELSE
281         K(IPA,1)=3
282         K(IPA+1,1)=3
283         K(IPA+2,1)=3
284         KCS=4
285         IF(KQ1.EQ.-1) KCS=5
286         K(IPA,KCS)=MSTU(5)*(IPA+1)
287         K(IPA,9-KCS)=MSTU(5)*(IPA+2)
288         K(IPA+1,KCS)=MSTU(5)*(IPA+2)
289         K(IPA+1,9-KCS)=MSTU(5)*IPA
290         K(IPA+2,KCS)=MSTU(5)*IPA
291         K(IPA+2,9-KCS)=MSTU(5)*(IPA+1)
292       ENDIF
293
294 C...Check kinematics.
295       MKERR=0
296       IF(0.5*X1*PECM.LE.PM1.OR.0.5*(2.-X1-X3)*PECM.LE.PM2.OR.
297      &0.5*X3*PECM.LE.PM3) MKERR=1
298       PA1=SQRT(MAX(1E-10,(0.5*X1*PECM)**2-PM1**2))
299       PA2=SQRT(MAX(1E-10,(0.5*(2.-X1-X3)*PECM)**2-PM2**2))
300       PA3=SQRT(MAX(1E-10,(0.5*X3*PECM)**2-PM3**2))
301       CTHE2=(PA3**2-PA1**2-PA2**2)/(2.*PA1*PA2)
302       CTHE3=(PA2**2-PA1**2-PA3**2)/(2.*PA1*PA3)
303       IF(ABS(CTHE2).GE.1.001.OR.ABS(CTHE3).GE.1.001) MKERR=1
304       CTHE3=MAX(-1.,MIN(1.,CTHE3))
305       IF(MKERR.NE.0) CALL LUERRM(13,
306      &'(LU3ENT:) unphysical kinematical variable setup')
307
308 C...Store partons/particles in P vectors.
309       P(IPA,3)=PA1
310       P(IPA,4)=SQRT(PA1**2+PM1**2)
311       P(IPA,5)=PM1
312       P(IPA+2,1)=PA3*SQRT(1.-CTHE3**2)
313       P(IPA+2,3)=PA3*CTHE3
314       P(IPA+2,4)=SQRT(PA3**2+PM3**2)
315       P(IPA+2,5)=PM3
316       P(IPA+1,1)=-P(IPA+2,1)
317       P(IPA+1,3)=-P(IPA,3)-P(IPA+2,3)
318       P(IPA+1,4)=SQRT(P(IPA+1,1)**2+P(IPA+1,3)**2+PM2**2)
319       P(IPA+1,5)=PM2
320
321 C...Set N. Optionally fragment/decay.
322       N=IPA+2
323       IF(IP.EQ.0) CALL LUEXEC
324
325       RETURN
326       END
327
328 C*********************************************************************
329
330       SUBROUTINE LU4ENT(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)
331
332 C...Purpose: to store four partons or particles in their CM frame, with
333 C...the first along the +z axis, the last in the xz plane with x > 0
334 C...and the second having y < 0 and y > 0 with equal probability.
335       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
336       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
337       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
338       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
339
340 C...Standard checks.
341       MSTU(28)=0
342       IF(MSTU(12).GE.1) CALL LULIST(0)
343       IPA=MAX(1,IABS(IP))
344       IF(IPA.GT.MSTU(4)-3) CALL LUERRM(21,
345      &'(LU4ENT:) writing outside LUJETS momory')
346       KC1=LUCOMP(KF1)
347       KC2=LUCOMP(KF2)
348       KC3=LUCOMP(KF3)
349       KC4=LUCOMP(KF4)
350       IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0.OR.KC4.EQ.0) CALL LUERRM(12,
351      &'(LU4ENT:) unknown flavour code')
352
353 C...Find masses. Reset K, P and V vectors.
354       PM1=0.
355       IF(MSTU(10).EQ.1) PM1=P(IPA,5)
356       IF(MSTU(10).GE.2) PM1=ULMASS(KF1)
357       PM2=0.
358       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
359       IF(MSTU(10).GE.2) PM2=ULMASS(KF2)
360       PM3=0.
361       IF(MSTU(10).EQ.1) PM3=P(IPA+2,5)
362       IF(MSTU(10).GE.2) PM3=ULMASS(KF3)
363       PM4=0.
364       IF(MSTU(10).EQ.1) PM4=P(IPA+3,5)
365       IF(MSTU(10).GE.2) PM4=ULMASS(KF4)
366       DO 100 I=IPA,IPA+3
367       DO 100 J=1,5
368       K(I,J)=0
369       P(I,J)=0.
370   100 V(I,J)=0.
371
372 C...Check flavours.
373       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
374       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
375       KQ3=KCHG(KC3,2)*ISIGN(1,KF3)
376       KQ4=KCHG(KC4,2)*ISIGN(1,KF4)
377       IF(MSTU(19).EQ.1) THEN
378         MSTU(19)=0
379       ELSEIF(KQ1.EQ.0.AND.KQ2.EQ.0.AND.KQ3.EQ.0.AND.KQ4.EQ.0) THEN
380       ELSEIF(KQ1.NE.0.AND.KQ2.EQ.2.AND.KQ3.EQ.2.AND.(KQ1+KQ4.EQ.0.OR.
381      &KQ1+KQ4.EQ.4)) THEN
382       ELSEIF(KQ1.NE.0.AND.KQ1+KQ2.EQ.0.AND.KQ3.NE.0.AND.KQ3+KQ4.EQ.0.)
383      &THEN
384       ELSE
385         CALL LUERRM(2,'(LU4ENT:) unphysical flavour combination')
386       ENDIF
387       K(IPA,2)=KF1
388       K(IPA+1,2)=KF2
389       K(IPA+2,2)=KF3
390       K(IPA+3,2)=KF4
391
392 C...Store partons/particles in K vectors for normal case.
393       IF(IP.GE.0) THEN
394         K(IPA,1)=1
395         IF(KQ1.NE.0.AND.(KQ2.NE.0.OR.KQ3.NE.0.OR.KQ4.NE.0)) K(IPA,1)=2
396         K(IPA+1,1)=1
397         IF(KQ2.NE.0.AND.KQ1+KQ2.NE.0.AND.(KQ3.NE.0.OR.KQ4.NE.0))
398      &  K(IPA+1,1)=2
399         K(IPA+2,1)=1
400         IF(KQ3.NE.0.AND.KQ4.NE.0) K(IPA+2,1)=2
401         K(IPA+3,1)=1
402
403 C...Store partons for parton shower evolution from q-g-g-qbar or
404 C...g-g-g-g event.
405       ELSEIF(KQ1+KQ2.NE.0) THEN
406         K(IPA,1)=3
407         K(IPA+1,1)=3
408         K(IPA+2,1)=3
409         K(IPA+3,1)=3
410         KCS=4
411         IF(KQ1.EQ.-1) KCS=5
412         K(IPA,KCS)=MSTU(5)*(IPA+1)
413         K(IPA,9-KCS)=MSTU(5)*(IPA+3)
414         K(IPA+1,KCS)=MSTU(5)*(IPA+2)
415         K(IPA+1,9-KCS)=MSTU(5)*IPA
416         K(IPA+2,KCS)=MSTU(5)*(IPA+3)
417         K(IPA+2,9-KCS)=MSTU(5)*(IPA+1)
418         K(IPA+3,KCS)=MSTU(5)*IPA
419         K(IPA+3,9-KCS)=MSTU(5)*(IPA+2)
420
421 C...Store partons for parton shower evolution from q-qbar-q-qbar event.
422       ELSE
423         K(IPA,1)=3
424         K(IPA+1,1)=3
425         K(IPA+2,1)=3
426         K(IPA+3,1)=3
427         K(IPA,4)=MSTU(5)*(IPA+1)
428         K(IPA,5)=K(IPA,4)
429         K(IPA+1,4)=MSTU(5)*IPA
430         K(IPA+1,5)=K(IPA+1,4)
431         K(IPA+2,4)=MSTU(5)*(IPA+3)
432         K(IPA+2,5)=K(IPA+2,4)
433         K(IPA+3,4)=MSTU(5)*(IPA+2)
434         K(IPA+3,5)=K(IPA+3,4)
435       ENDIF
436
437 C...Check kinematics.
438       MKERR=0
439       IF(0.5*X1*PECM.LE.PM1.OR.0.5*X2*PECM.LE.PM2.OR.0.5*(2.-X1-X2-X4)*
440      &PECM.LE.PM3.OR.0.5*X4*PECM.LE.PM4) MKERR=1
441       PA1=SQRT(MAX(1E-10,(0.5*X1*PECM)**2-PM1**2))
442       PA2=SQRT(MAX(1E-10,(0.5*X2*PECM)**2-PM2**2))
443       PA4=SQRT(MAX(1E-10,(0.5*X4*PECM)**2-PM4**2))
444       X24=X1+X2+X4-1.-X12-X14+(PM3**2-PM1**2-PM2**2-PM4**2)/PECM**2
445       CTHE4=(X1*X4-2.*X14)*PECM**2/(4.*PA1*PA4)
446       IF(ABS(CTHE4).GE.1.002) MKERR=1
447       CTHE4=MAX(-1.,MIN(1.,CTHE4))
448       STHE4=SQRT(1.-CTHE4**2)
449       CTHE2=(X1*X2-2.*X12)*PECM**2/(4.*PA1*PA2)
450       IF(ABS(CTHE2).GE.1.002) MKERR=1
451       CTHE2=MAX(-1.,MIN(1.,CTHE2))
452       STHE2=SQRT(1.-CTHE2**2)
453       CPHI2=((X2*X4-2.*X24)*PECM**2-4.*PA2*CTHE2*PA4*CTHE4)/
454      &MAX(1E-8*PECM**2,4.*PA2*STHE2*PA4*STHE4)
455       IF(ABS(CPHI2).GE.1.05) MKERR=1
456       CPHI2=MAX(-1.,MIN(1.,CPHI2))
457       IF(MKERR.EQ.1) CALL LUERRM(13,
458      &'(LU4ENT:) unphysical kinematical variable setup')
459
460 C...Store partons/particles in P vectors.
461       P(IPA,3)=PA1
462       P(IPA,4)=SQRT(PA1**2+PM1**2)
463       P(IPA,5)=PM1
464       P(IPA+3,1)=PA4*STHE4
465       P(IPA+3,3)=PA4*CTHE4
466       P(IPA+3,4)=SQRT(PA4**2+PM4**2)
467       P(IPA+3,5)=PM4
468       P(IPA+1,1)=PA2*STHE2*CPHI2
469       P(IPA+1,2)=PA2*STHE2*SQRT(1.-CPHI2**2)*(-1.)**INT(RLU(0)+0.5)
470       P(IPA+1,3)=PA2*CTHE2
471       P(IPA+1,4)=SQRT(PA2**2+PM2**2)
472       P(IPA+1,5)=PM2
473       P(IPA+2,1)=-P(IPA+1,1)-P(IPA+3,1)
474       P(IPA+2,2)=-P(IPA+1,2)
475       P(IPA+2,3)=-P(IPA,3)-P(IPA+1,3)-P(IPA+3,3)
476       P(IPA+2,4)=SQRT(P(IPA+2,1)**2+P(IPA+2,2)**2+P(IPA+2,3)**2+PM3**2)
477       P(IPA+2,5)=PM3
478
479 C...Set N. Optionally fragment/decay.
480       N=IPA+3
481       IF(IP.EQ.0) CALL LUEXEC
482
483       RETURN
484       END
485
486 C*********************************************************************
487
488       SUBROUTINE LUJOIN(NJOIN,IJOIN)
489
490 C...Purpose: to connect a sequence of partons with colour flow indices,
491 C...as required for subsequent shower evolution (or other operations).
492       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
493       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
494       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
495       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
496       DIMENSION IJOIN(*)
497
498 C...Check that partons are of right types to be connected.
499       IF(NJOIN.LT.2) GOTO 120
500       KQSUM=0
501       DO 100 IJN=1,NJOIN
502       I=IJOIN(IJN)
503       IF(I.LE.0.OR.I.GT.N) GOTO 120
504       IF(K(I,1).LT.1.OR.K(I,1).GT.3) GOTO 120
505       KC=LUCOMP(K(I,2))
506       IF(KC.EQ.0) GOTO 120
507       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
508       IF(KQ.EQ.0) GOTO 120
509       IF(IJN.NE.1.AND.IJN.NE.NJOIN.AND.KQ.NE.2) GOTO 120
510       IF(KQ.NE.2) KQSUM=KQSUM+KQ
511   100 IF(IJN.EQ.1) KQS=KQ
512       IF(KQSUM.NE.0) GOTO 120
513
514 C...Connect the partons sequentially (closing for gluon loop).
515       KCS=(9-KQS)/2
516       IF(KQS.EQ.2) KCS=INT(4.5+RLU(0))
517       DO 110 IJN=1,NJOIN
518       I=IJOIN(IJN)
519       K(I,1)=3
520       IF(IJN.NE.1) IP=IJOIN(IJN-1)
521       IF(IJN.EQ.1) IP=IJOIN(NJOIN)
522       IF(IJN.NE.NJOIN) IN=IJOIN(IJN+1)
523       IF(IJN.EQ.NJOIN) IN=IJOIN(1)
524       K(I,KCS)=MSTU(5)*IN
525       K(I,9-KCS)=MSTU(5)*IP
526       IF(IJN.EQ.1.AND.KQS.NE.2) K(I,9-KCS)=0
527   110 IF(IJN.EQ.NJOIN.AND.KQS.NE.2) K(I,KCS)=0
528
529 C...Error exit: no action taken.
530       RETURN
531   120 CALL LUERRM(12,
532      &'(LUJOIN:) given entries can not be joined by one string')
533
534       RETURN
535       END
536
537 C*********************************************************************
538
539       SUBROUTINE LUGIVE(CHIN)
540
541 C...Purpose: to set values of commonblock variables (also in PYTHIA!).
542       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
543       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
544       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
545       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
546       COMMON/LUDAT4/CHAF(500)
547       CHARACTER CHAF*8
548       COMMON/LUDATR/MRLU(6),RRLU(100)
549       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
550       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
551       COMMON/PYINT1/MINT(400),VINT(400)
552       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
553       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
554       COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
555       COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
556       COMMON/PYINT6/PROC(0:200)
557       CHARACTER PROC*28
558       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/,/LUDAT4/,/LUDATR/
559       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,
560      &/PYINT5/,/PYINT6/
561       CHARACTER CHIN*(*),CHFIX*104,CHBIT*104,CHOLD*8,CHNEW*8,CHOLD2*28,
562      &CHNEW2*28,CHNAM*4,CHVAR(42)*4,CHALP(2)*26,CHIND*8,CHINI*10,
563      &CHINR*16
564       DIMENSION MSVAR(42,8)
565
566 C...For each variable to be translated give: name,
567 C...integer/real/character, no. of indices, lower&upper index bounds.
568       DATA CHVAR/'N','K','P','V','MSTU','PARU','MSTJ','PARJ','KCHG',
569      &'PMAS','PARF','VCKM','MDCY','MDME','BRAT','KFDP','CHAF','MRLU',
570      &'RRLU','MSEL','MSUB','KFIN','CKIN','MSTP','PARP','MSTI','PARI',
571      &'MINT','VINT','ISET','KFPR','COEF','ICOL','XSFX','ISIG','SIGH',
572      &'WIDP','WIDE','WIDS','NGEN','XSEC','PROC'/
573       DATA ((MSVAR(I,J),J=1,8),I=1,42)/ 1,7*0,  1,2,1,4000,1,5,2*0,
574      & 2,2,1,4000,1,5,2*0,  2,2,1,4000,1,5,2*0,  1,1,1,200,4*0,
575      & 2,1,1,200,4*0,  1,1,1,200,4*0,  2,1,1,200,4*0,
576      & 1,2,1,500,1,3,2*0,  2,2,1,500,1,4,2*0,  2,1,1,2000,4*0,
577      & 2,2,1,4,1,4,2*0,  1,2,1,500,1,3,2*0,  1,2,1,2000,1,2,2*0,
578      & 2,1,1,2000,4*0,  1,2,1,2000,1,5,2*0,  3,1,1,500,4*0,
579      & 1,1,1,6,4*0,  2,1,1,100,4*0,
580      & 1,7*0,  1,1,1,200,4*0,  1,2,1,2,-40,40,2*0,  2,1,1,200,4*0,
581      & 1,1,1,200,4*0,  2,1,1,200,4*0,  1,1,1,200,4*0,  2,1,1,200,4*0,
582      & 1,1,1,400,4*0,  2,1,1,400,4*0,  1,1,1,200,4*0,
583      & 1,2,1,200,1,2,2*0,  2,2,1,200,1,20,2*0,  1,3,1,40,1,4,1,2,
584      & 2,2,1,2,-40,40,2*0,  1,2,1,1000,1,3,2*0,  2,1,1,1000,4*0,
585      & 2,2,21,40,0,40,2*0,  2,2,21,40,0,40,2*0,  2,2,21,40,1,3,2*0,
586      & 1,2,0,200,1,3,2*0,  2,2,0,200,1,3,2*0,  4,1,0,200,4*0/
587       DATA CHALP/'abcdefghijklmnopqrstuvwxyz',
588      &'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
589
590 C...Length of character variable. Subdivide it into instructions.
591       IF(MSTU(12).GE.1) CALL LULIST(0)
592       CHBIT=CHIN//' '
593       LBIT=101
594   100 LBIT=LBIT-1
595       IF(CHBIT(LBIT:LBIT).EQ.' ') GOTO 100
596       LTOT=0
597       DO 110 LCOM=1,LBIT
598       IF(CHBIT(LCOM:LCOM).EQ.' ') GOTO 110
599       LTOT=LTOT+1
600       CHFIX(LTOT:LTOT)=CHBIT(LCOM:LCOM)
601   110 CONTINUE
602       LLOW=0
603   120 LHIG=LLOW+1
604   130 LHIG=LHIG+1
605       IF(LHIG.LE.LTOT.AND.CHFIX(LHIG:LHIG).NE.';') GOTO 130
606       LBIT=LHIG-LLOW-1
607       CHBIT(1:LBIT)=CHFIX(LLOW+1:LHIG-1)
608
609 C...Identify commonblock variable.
610       LNAM=1
611   140 LNAM=LNAM+1
612       IF(CHBIT(LNAM:LNAM).NE.'('.AND.CHBIT(LNAM:LNAM).NE.'='.AND.
613      &LNAM.LE.4) GOTO 140
614       CHNAM=CHBIT(1:LNAM-1)//' '
615       DO 150 LCOM=1,LNAM-1
616       DO 150 LALP=1,26
617   150 IF(CHNAM(LCOM:LCOM).EQ.CHALP(1)(LALP:LALP)) CHNAM(LCOM:LCOM)=
618      &CHALP(2)(LALP:LALP)
619       IVAR=0
620       DO 160 IV=1,42
621   160 IF(CHNAM.EQ.CHVAR(IV)) IVAR=IV
622       IF(IVAR.EQ.0) THEN
623         CALL LUERRM(18,'(LUGIVE:) do not recognize variable '//CHNAM)
624         LLOW=LHIG
625         IF(LLOW.LT.LTOT) GOTO 120
626         RETURN
627       ENDIF
628
629 C...Identify any indices.
630       I1=0
631       I2=0
632       I3=0
633       NINDX=0
634       IF(CHBIT(LNAM:LNAM).EQ.'(') THEN
635         LIND=LNAM
636   170   LIND=LIND+1
637         IF(CHBIT(LIND:LIND).NE.')'.AND.CHBIT(LIND:LIND).NE.',') GOTO 170
638         CHIND=' '
639         IF((CHBIT(LNAM+1:LNAM+1).EQ.'C'.OR.CHBIT(LNAM+1:LNAM+1).EQ.'c').
640      &  AND.(IVAR.EQ.9.OR.IVAR.EQ.10.OR.IVAR.EQ.13.OR.IVAR.EQ.17)) THEN
641           CHIND(LNAM-LIND+11:8)=CHBIT(LNAM+2:LIND-1)
642           READ(CHIND,'(I8)') KF
643           I1=LUCOMP(KF)
644         ELSEIF(CHBIT(LNAM+1:LNAM+1).EQ.'C'.OR.CHBIT(LNAM+1:LNAM+1).EQ.
645      &  'c') THEN
646           CALL LUERRM(18,'(LUGIVE:) not allowed to use C index for '//
647      &    CHNAM)
648           LLOW=LHIG
649           IF(LLOW.LT.LTOT) GOTO 120
650           RETURN
651         ELSE
652           CHIND(LNAM-LIND+10:8)=CHBIT(LNAM+1:LIND-1)
653           READ(CHIND,'(I8)') I1
654         ENDIF
655         LNAM=LIND
656         IF(CHBIT(LNAM:LNAM).EQ.')') LNAM=LNAM+1
657         NINDX=1
658       ENDIF
659       IF(CHBIT(LNAM:LNAM).EQ.',') THEN
660         LIND=LNAM
661   180   LIND=LIND+1
662         IF(CHBIT(LIND:LIND).NE.')'.AND.CHBIT(LIND:LIND).NE.',') GOTO 180
663         CHIND=' '
664         CHIND(LNAM-LIND+10:8)=CHBIT(LNAM+1:LIND-1)
665         READ(CHIND,'(I8)') I2
666         LNAM=LIND
667         IF(CHBIT(LNAM:LNAM).EQ.')') LNAM=LNAM+1
668         NINDX=2
669       ENDIF
670       IF(CHBIT(LNAM:LNAM).EQ.',') THEN
671         LIND=LNAM
672   190   LIND=LIND+1
673         IF(CHBIT(LIND:LIND).NE.')'.AND.CHBIT(LIND:LIND).NE.',') GOTO 190
674         CHIND=' '
675         CHIND(LNAM-LIND+10:8)=CHBIT(LNAM+1:LIND-1)
676         READ(CHIND,'(I8)') I3
677         LNAM=LIND+1
678         NINDX=3
679       ENDIF
680
681 C...Check that indices allowed.
682       IERR=0
683       IF(NINDX.NE.MSVAR(IVAR,2)) IERR=1
684       IF(NINDX.GE.1.AND.(I1.LT.MSVAR(IVAR,3).OR.I1.GT.MSVAR(IVAR,4)))
685      &IERR=2
686       IF(NINDX.GE.2.AND.(I2.LT.MSVAR(IVAR,5).OR.I2.GT.MSVAR(IVAR,6)))
687      &IERR=3
688       IF(NINDX.EQ.3.AND.(I3.LT.MSVAR(IVAR,7).OR.I3.GT.MSVAR(IVAR,8)))
689      &IERR=4
690       IF(CHBIT(LNAM:LNAM).NE.'=') IERR=5
691       IF(IERR.GE.1) THEN
692         CALL LUERRM(18,'(LUGIVE:) unallowed indices for '//
693      &  CHBIT(1:LNAM-1))
694         LLOW=LHIG
695         IF(LLOW.LT.LTOT) GOTO 120
696         RETURN
697       ENDIF
698
699 C...Save old value of variable.
700       IF(IVAR.EQ.1) THEN
701         IOLD=N
702       ELSEIF(IVAR.EQ.2) THEN
703         IOLD=K(I1,I2)
704       ELSEIF(IVAR.EQ.3) THEN
705         ROLD=P(I1,I2)
706       ELSEIF(IVAR.EQ.4) THEN
707         ROLD=V(I1,I2)
708       ELSEIF(IVAR.EQ.5) THEN
709         IOLD=MSTU(I1)
710       ELSEIF(IVAR.EQ.6) THEN
711         ROLD=PARU(I1)
712       ELSEIF(IVAR.EQ.7) THEN
713         IOLD=MSTJ(I1)
714       ELSEIF(IVAR.EQ.8) THEN
715         ROLD=PARJ(I1)
716       ELSEIF(IVAR.EQ.9) THEN
717         IOLD=KCHG(I1,I2)
718       ELSEIF(IVAR.EQ.10) THEN
719         ROLD=PMAS(I1,I2)
720       ELSEIF(IVAR.EQ.11) THEN
721         ROLD=PARF(I1)
722       ELSEIF(IVAR.EQ.12) THEN
723         ROLD=VCKM(I1,I2)
724       ELSEIF(IVAR.EQ.13) THEN
725         IOLD=MDCY(I1,I2)
726       ELSEIF(IVAR.EQ.14) THEN
727         IOLD=MDME(I1,I2)
728       ELSEIF(IVAR.EQ.15) THEN
729         ROLD=BRAT(I1)
730       ELSEIF(IVAR.EQ.16) THEN
731         IOLD=KFDP(I1,I2)
732       ELSEIF(IVAR.EQ.17) THEN
733         CHOLD=CHAF(I1)
734       ELSEIF(IVAR.EQ.18) THEN
735         IOLD=MRLU(I1)
736       ELSEIF(IVAR.EQ.19) THEN
737         ROLD=RRLU(I1)
738       ELSEIF(IVAR.EQ.20) THEN
739         IOLD=MSEL
740       ELSEIF(IVAR.EQ.21) THEN
741         IOLD=MSUB(I1)
742       ELSEIF(IVAR.EQ.22) THEN
743         IOLD=KFIN(I1,I2)
744       ELSEIF(IVAR.EQ.23) THEN
745         ROLD=CKIN(I1)
746       ELSEIF(IVAR.EQ.24) THEN
747         IOLD=MSTP(I1)
748       ELSEIF(IVAR.EQ.25) THEN
749         ROLD=PARP(I1)
750       ELSEIF(IVAR.EQ.26) THEN
751         IOLD=MSTI(I1)
752       ELSEIF(IVAR.EQ.27) THEN
753         ROLD=PARI(I1)
754       ELSEIF(IVAR.EQ.28) THEN
755         IOLD=MINT(I1)
756       ELSEIF(IVAR.EQ.29) THEN
757         ROLD=VINT(I1)
758       ELSEIF(IVAR.EQ.30) THEN
759         IOLD=ISET(I1)
760       ELSEIF(IVAR.EQ.31) THEN
761         IOLD=KFPR(I1,I2)
762       ELSEIF(IVAR.EQ.32) THEN
763         ROLD=COEF(I1,I2)
764       ELSEIF(IVAR.EQ.33) THEN
765         IOLD=ICOL(I1,I2,I3)
766       ELSEIF(IVAR.EQ.34) THEN
767         ROLD=XSFX(I1,I2)
768       ELSEIF(IVAR.EQ.35) THEN
769         IOLD=ISIG(I1,I2)
770       ELSEIF(IVAR.EQ.36) THEN
771         ROLD=SIGH(I1)
772       ELSEIF(IVAR.EQ.37) THEN
773         ROLD=WIDP(I1,I2)
774       ELSEIF(IVAR.EQ.38) THEN
775         ROLD=WIDE(I1,I2)
776       ELSEIF(IVAR.EQ.39) THEN
777         ROLD=WIDS(I1,I2)
778       ELSEIF(IVAR.EQ.40) THEN
779         IOLD=NGEN(I1,I2)
780       ELSEIF(IVAR.EQ.41) THEN
781         ROLD=XSEC(I1,I2)
782       ELSEIF(IVAR.EQ.42) THEN
783         CHOLD2=PROC(I1)
784       ENDIF
785
786 C...Print current value of variable. Loop back.
787       IF(LNAM.GE.LBIT) THEN
788         CHBIT(LNAM:14)=' '
789         CHBIT(15:60)=' has the value                                '
790         IF(MSVAR(IVAR,1).EQ.1) THEN
791           WRITE(CHBIT(51:60),'(I10)') IOLD
792         ELSEIF(MSVAR(IVAR,1).EQ.2) THEN
793           WRITE(CHBIT(47:60),'(F14.5)') ROLD
794         ELSEIF(MSVAR(IVAR,1).EQ.3) THEN
795           CHBIT(53:60)=CHOLD
796         ELSE
797           CHBIT(33:60)=CHOLD
798         ENDIF
799         IF(MSTU(13).GE.1) WRITE(MSTU(11),5000) CHBIT(1:60)
800         LLOW=LHIG
801         IF(LLOW.LT.LTOT) GOTO 120
802         RETURN
803       ENDIF
804
805 C...Read in new variable value.
806       IF(MSVAR(IVAR,1).EQ.1) THEN
807         CHINI=' '
808         CHINI(LNAM-LBIT+11:10)=CHBIT(LNAM+1:LBIT)
809         READ(CHINI,'(I10)') INEW
810       ELSEIF(MSVAR(IVAR,1).EQ.2) THEN
811         CHINR=' '
812         CHINR(LNAM-LBIT+17:16)=CHBIT(LNAM+1:LBIT)
813         READ(CHINR,'(F16.2)') RNEW
814       ELSEIF(MSVAR(IVAR,1).EQ.3) THEN
815         CHNEW=CHBIT(LNAM+1:LBIT)//' '
816       ELSE
817         CHNEW2=CHBIT(LNAM+1:LBIT)//' '
818       ENDIF
819
820 C...Store new variable value.
821       IF(IVAR.EQ.1) THEN
822         N=INEW
823       ELSEIF(IVAR.EQ.2) THEN
824         K(I1,I2)=INEW
825       ELSEIF(IVAR.EQ.3) THEN
826         P(I1,I2)=RNEW
827       ELSEIF(IVAR.EQ.4) THEN
828         V(I1,I2)=RNEW
829       ELSEIF(IVAR.EQ.5) THEN
830         MSTU(I1)=INEW
831       ELSEIF(IVAR.EQ.6) THEN
832         PARU(I1)=RNEW
833       ELSEIF(IVAR.EQ.7) THEN
834         MSTJ(I1)=INEW
835       ELSEIF(IVAR.EQ.8) THEN
836         PARJ(I1)=RNEW
837       ELSEIF(IVAR.EQ.9) THEN
838         KCHG(I1,I2)=INEW
839       ELSEIF(IVAR.EQ.10) THEN
840         PMAS(I1,I2)=RNEW
841       ELSEIF(IVAR.EQ.11) THEN
842         PARF(I1)=RNEW
843       ELSEIF(IVAR.EQ.12) THEN
844         VCKM(I1,I2)=RNEW
845       ELSEIF(IVAR.EQ.13) THEN
846         MDCY(I1,I2)=INEW
847       ELSEIF(IVAR.EQ.14) THEN
848         MDME(I1,I2)=INEW
849       ELSEIF(IVAR.EQ.15) THEN
850         BRAT(I1)=RNEW
851       ELSEIF(IVAR.EQ.16) THEN
852         KFDP(I1,I2)=INEW
853       ELSEIF(IVAR.EQ.17) THEN
854         CHAF(I1)=CHNEW
855       ELSEIF(IVAR.EQ.18) THEN
856         MRLU(I1)=INEW
857       ELSEIF(IVAR.EQ.19) THEN
858         RRLU(I1)=RNEW
859       ELSEIF(IVAR.EQ.20) THEN
860         MSEL=INEW
861       ELSEIF(IVAR.EQ.21) THEN
862         MSUB(I1)=INEW
863       ELSEIF(IVAR.EQ.22) THEN
864         KFIN(I1,I2)=INEW
865       ELSEIF(IVAR.EQ.23) THEN
866         CKIN(I1)=RNEW
867       ELSEIF(IVAR.EQ.24) THEN
868         MSTP(I1)=INEW
869       ELSEIF(IVAR.EQ.25) THEN
870         PARP(I1)=RNEW
871       ELSEIF(IVAR.EQ.26) THEN
872         MSTI(I1)=INEW
873       ELSEIF(IVAR.EQ.27) THEN
874         PARI(I1)=RNEW
875       ELSEIF(IVAR.EQ.28) THEN
876         MINT(I1)=INEW
877       ELSEIF(IVAR.EQ.29) THEN
878         VINT(I1)=RNEW
879       ELSEIF(IVAR.EQ.30) THEN
880         ISET(I1)=INEW
881       ELSEIF(IVAR.EQ.31) THEN
882         KFPR(I1,I2)=INEW
883       ELSEIF(IVAR.EQ.32) THEN
884         COEF(I1,I2)=RNEW
885       ELSEIF(IVAR.EQ.33) THEN
886         ICOL(I1,I2,I3)=INEW
887       ELSEIF(IVAR.EQ.34) THEN
888         XSFX(I1,I2)=RNEW
889       ELSEIF(IVAR.EQ.35) THEN
890         ISIG(I1,I2)=INEW
891       ELSEIF(IVAR.EQ.36) THEN
892         SIGH(I1)=RNEW
893       ELSEIF(IVAR.EQ.37) THEN
894         WIDP(I1,I2)=RNEW
895       ELSEIF(IVAR.EQ.38) THEN
896         WIDE(I1,I2)=RNEW
897       ELSEIF(IVAR.EQ.39) THEN
898         WIDS(I1,I2)=RNEW
899       ELSEIF(IVAR.EQ.40) THEN
900         NGEN(I1,I2)=INEW
901       ELSEIF(IVAR.EQ.41) THEN
902         XSEC(I1,I2)=RNEW
903       ELSEIF(IVAR.EQ.42) THEN
904         PROC(I1)=CHNEW2
905       ENDIF
906
907 C...Write old and new value. Loop back.
908       CHBIT(LNAM:14)=' '
909       CHBIT(15:60)=' changed from                to               '
910       IF(MSVAR(IVAR,1).EQ.1) THEN
911         WRITE(CHBIT(33:42),'(I10)') IOLD
912         WRITE(CHBIT(51:60),'(I10)') INEW
913         IF(MSTU(13).GE.1) WRITE(MSTU(11),5000) CHBIT(1:60)
914       ELSEIF(MSVAR(IVAR,1).EQ.2) THEN
915         WRITE(CHBIT(29:42),'(F14.5)') ROLD
916         WRITE(CHBIT(47:60),'(F14.5)') RNEW
917         IF(MSTU(13).GE.1) WRITE(MSTU(11),5000) CHBIT(1:60)
918       ELSEIF(MSVAR(IVAR,1).EQ.3) THEN
919         CHBIT(35:42)=CHOLD
920         CHBIT(53:60)=CHNEW
921         IF(MSTU(13).GE.1) WRITE(MSTU(11),5000) CHBIT(1:60)
922       ELSE
923         CHBIT(15:88)=' changed from '//CHOLD2//' to '//CHNEW2
924         IF(MSTU(13).GE.1) WRITE(MSTU(11),5100) CHBIT(1:88)
925       ENDIF
926       LLOW=LHIG
927       IF(LLOW.LT.LTOT) GOTO 120
928
929 C...Format statement for output on unit MSTU(11) (by default 6).
930  5000 FORMAT(5X,A60)
931  5100 FORMAT(5X,A88)
932
933       RETURN
934       END
935
936 C*********************************************************************
937
938       SUBROUTINE LUEXEC
939
940 C...Purpose: to administrate the fragmentation and decay chain.
941       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
942       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
943       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
944       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
945       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
946       DIMENSION PS(2,6)
947
948 C...Initialize and reset.
949       MSTU(24)=0
950       IF(MSTU(12).GE.1) CALL LULIST(0)
951       MSTU(31)=MSTU(31)+1
952       MSTU(1)=0
953       MSTU(2)=0
954       MSTU(3)=0
955       IF(MSTU(17).LE.0) MSTU(90)=0
956       MCONS=1
957
958 C...Sum up momentum, energy and charge for starting entries.
959       NSAV=N
960       DO 100 I=1,2
961       DO 100 J=1,6
962   100 PS(I,J)=0.
963       DO 120 I=1,N
964       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120
965       DO 110 J=1,4
966   110 PS(1,J)=PS(1,J)+P(I,J)
967       PS(1,6)=PS(1,6)+LUCHGE(K(I,2))
968   120 CONTINUE
969       PARU(21)=PS(1,4)
970
971 C...Prepare system for subsequent fragmentation/decay.
972       CALL LUPREP(0)
973
974 C...Loop through jet fragmentation and particle decays.
975       MBE=0
976   130 MBE=MBE+1
977       IP=0
978   140 IP=IP+1
979       KC=0
980       IF(K(IP,1).GT.0.AND.K(IP,1).LE.10) KC=LUCOMP(K(IP,2))
981       IF(KC.EQ.0) THEN
982
983 C...Particle decay if unstable and allowed. Save long-lived particle
984 C...decays until second pass after Bose-Einstein effects.
985       ELSEIF(KCHG(KC,2).EQ.0) THEN
986         IF(MSTJ(21).GE.1.AND.MDCY(KC,1).GE.1.AND.(MSTJ(51).LE.0.OR.MBE.
987      &  EQ.2.OR.PMAS(KC,2).GE.PARJ(91).OR.IABS(K(IP,2)).EQ.311))
988      &  CALL LUDECY(IP)
989
990 C...Decay products may develop a shower.
991         IF(MSTJ(92).GT.0) THEN
992           IP1=MSTJ(92)
993           QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,
994      &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))
995           CALL LUSHOW(IP1,IP1+1,QMAX)
996           CALL LUPREP(IP1)
997           MSTJ(92)=0
998         ELSEIF(MSTJ(92).LT.0) THEN
999           IP1=-MSTJ(92)
1000           CALL LUSHOW(IP1,-3,P(IP,5))
1001           CALL LUPREP(IP1)
1002           MSTJ(92)=0
1003         ENDIF
1004
1005 C...Jet fragmentation: string or independent fragmentation.
1006       ELSEIF(K(IP,1).EQ.1.OR.K(IP,1).EQ.2) THEN
1007         MFRAG=MSTJ(1)
1008         IF(MFRAG.GE.1.AND.K(IP,1).EQ.1) MFRAG=2
1009         IF(MSTJ(21).GE.2.AND.K(IP,1).EQ.2.AND.N.GT.IP) THEN
1010           IF(K(IP+1,1).EQ.1.AND.K(IP+1,3).EQ.K(IP,3).AND.
1011      &    K(IP,3).GT.0.AND.K(IP,3).LT.IP) THEN
1012             IF(KCHG(LUCOMP(K(K(IP,3),2)),2).EQ.0) MFRAG=MIN(1,MFRAG)
1013           ENDIF
1014         ENDIF
1015         IF(MFRAG.EQ.1) CALL LUSTRF(IP)
1016         IF(MFRAG.EQ.2) CALL LUINDF(IP)
1017         IF(MFRAG.EQ.2.AND.K(IP,1).EQ.1) MCONS=0
1018         IF(MFRAG.EQ.2.AND.(MSTJ(3).LE.0.OR.MOD(MSTJ(3),5).EQ.0)) MCONS=0
1019       ENDIF
1020
1021 C...Loop back if enough space left in LUJETS and no error abort.
1022       IF(MSTU(24).NE.0.AND.MSTU(21).GE.2) THEN
1023       ELSEIF(IP.LT.N.AND.N.LT.MSTU(4)-20-MSTU(32)) THEN
1024         GOTO 140
1025       ELSEIF(IP.LT.N) THEN
1026         CALL LUERRM(11,'(LUEXEC:) no more memory left in LUJETS')
1027       ENDIF
1028
1029 C...Include simple Bose-Einstein effect parametrization if desired.
1030       IF(MBE.EQ.1.AND.MSTJ(51).GE.1) THEN
1031         CALL LUBOEI(NSAV)
1032         GOTO 130
1033       ENDIF
1034
1035 C...Check that momentum, energy and charge were conserved.
1036       DO 160 I=1,N
1037       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 160
1038       DO 150 J=1,4
1039   150 PS(2,J)=PS(2,J)+P(I,J)
1040       PS(2,6)=PS(2,6)+LUCHGE(K(I,2))
1041   160 CONTINUE
1042       PDEV=(ABS(PS(2,1)-PS(1,1))+ABS(PS(2,2)-PS(1,2))+ABS(PS(2,3)-
1043      &PS(1,3))+ABS(PS(2,4)-PS(1,4)))/(1.+ABS(PS(2,4))+ABS(PS(1,4)))
1044       IF(MCONS.EQ.1.AND.PDEV.GT.PARU(11)) CALL LUERRM(15,
1045      &'(LUEXEC:) four-momentum was not conserved')
1046       IF(MCONS.EQ.1.AND.ABS(PS(2,6)-PS(1,6)).GT.0.1) CALL LUERRM(15,
1047      &'(LUEXEC:) charge was not conserved')
1048
1049       RETURN
1050       END
1051
1052 C*********************************************************************
1053
1054       SUBROUTINE LUPREP(IP)
1055
1056 C...Purpose: to rearrange partons along strings, to allow small systems
1057 C...to collapse into one or two particles and to check flavours.
1058       IMPLICIT DOUBLE PRECISION(D)
1059       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
1060       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1061       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
1062       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
1063       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
1064       DIMENSION DPS(5),DPC(5),UE(3)
1065
1066 C...Rearrange parton shower product listing along strings: begin loop.
1067       I1=N
1068       DO 130 MQGST=1,2
1069       DO 120 I=MAX(1,IP),N
1070       IF(K(I,1).NE.3) GOTO 120
1071       KC=LUCOMP(K(I,2))
1072       IF(KC.EQ.0) GOTO 120
1073       KQ=KCHG(KC,2)
1074       IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 120
1075
1076 C...Pick up loose string end.
1077       KCS=4
1078       IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
1079       IA=I
1080       NSTP=0
1081   100 NSTP=NSTP+1
1082       IF(NSTP.GT.4*N) THEN
1083         CALL LUERRM(14,'(LUPREP:) caught in infinite loop')
1084         RETURN
1085       ENDIF
1086
1087 C...Copy undecayed parton.
1088       IF(K(IA,1).EQ.3) THEN
1089         IF(I1.GE.MSTU(4)-MSTU(32)-5) THEN
1090           CALL LUERRM(11,'(LUPREP:) no more memory left in LUJETS')
1091           RETURN
1092         ENDIF
1093         I1=I1+1
1094         K(I1,1)=2
1095         IF(NSTP.GE.2.AND.IABS(K(IA,2)).NE.21) K(I1,1)=1
1096         K(I1,2)=K(IA,2)
1097         K(I1,3)=IA
1098         K(I1,4)=0
1099         K(I1,5)=0
1100         DO 110 J=1,5
1101         P(I1,J)=P(IA,J)
1102   110   V(I1,J)=V(IA,J)
1103         K(IA,1)=K(IA,1)+10
1104         IF(K(I1,1).EQ.1) GOTO 120
1105       ENDIF
1106
1107 C...Go to next parton in colour space.
1108       IB=IA
1109       IF(MOD(K(IB,KCS)/MSTU(5)**2,2).EQ.0.AND.MOD(K(IB,KCS),MSTU(5)).
1110      &NE.0) THEN
1111         IA=MOD(K(IB,KCS),MSTU(5))
1112         K(IB,KCS)=K(IB,KCS)+MSTU(5)**2
1113         MREV=0
1114       ELSE
1115         IF(K(IB,KCS).GE.2*MSTU(5)**2.OR.MOD(K(IB,KCS)/MSTU(5),MSTU(5)).
1116      &  EQ.0) KCS=9-KCS
1117         IA=MOD(K(IB,KCS)/MSTU(5),MSTU(5))
1118         K(IB,KCS)=K(IB,KCS)+2*MSTU(5)**2
1119         MREV=1
1120       ENDIF
1121       IF(IA.LE.0.OR.IA.GT.N) THEN
1122         CALL LUERRM(12,'(LUPREP:) colour rearrangement failed')
1123         RETURN
1124       ENDIF
1125       IF(MOD(K(IA,4)/MSTU(5),MSTU(5)).EQ.IB.OR.MOD(K(IA,5)/MSTU(5),
1126      &MSTU(5)).EQ.IB) THEN
1127         IF(MREV.EQ.1) KCS=9-KCS
1128         IF(MOD(K(IA,KCS)/MSTU(5),MSTU(5)).NE.IB) KCS=9-KCS
1129         K(IA,KCS)=K(IA,KCS)+2*MSTU(5)**2
1130       ELSE
1131         IF(MREV.EQ.0) KCS=9-KCS
1132         IF(MOD(K(IA,KCS),MSTU(5)).NE.IB) KCS=9-KCS
1133         K(IA,KCS)=K(IA,KCS)+MSTU(5)**2
1134       ENDIF
1135       IF(IA.NE.I) GOTO 100
1136       K(I1,1)=1
1137   120 CONTINUE
1138   130 CONTINUE
1139       N=I1
1140       IF(MSTJ(14).LT.0) RETURN
1141
1142 C...Find lowest-mass colour singlet jet system, OK if above threshold.
1143       IF(MSTJ(14).EQ.0) GOTO 320
1144       NS=N
1145   140 NSIN=N-NS
1146       PDM=1.+PARJ(32)
1147       IC=0
1148       DO 190 I=MAX(1,IP),NS
1149       IF(K(I,1).NE.1.AND.K(I,1).NE.2) THEN
1150       ELSEIF(K(I,1).EQ.2.AND.IC.EQ.0) THEN
1151         NSIN=NSIN+1
1152         IC=I
1153         DO 150 J=1,4
1154   150   DPS(J)=P(I,J)
1155         MSTJ(93)=1
1156         DPS(5)=ULMASS(K(I,2))
1157       ELSEIF(K(I,1).EQ.2) THEN
1158         DO 160 J=1,4
1159   160   DPS(J)=DPS(J)+P(I,J)
1160       ELSEIF(IC.NE.0.AND.KCHG(LUCOMP(K(I,2)),2).NE.0) THEN
1161         DO 170 J=1,4
1162   170   DPS(J)=DPS(J)+P(I,J)
1163         MSTJ(93)=1
1164         DPS(5)=DPS(5)+ULMASS(K(I,2))
1165         PD=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))-DPS(5)
1166         IF(PD.LT.PDM) THEN
1167           PDM=PD
1168           DO 180 J=1,5
1169   180     DPC(J)=DPS(J)
1170           IC1=IC
1171           IC2=I
1172         ENDIF
1173         IC=0
1174       ELSE
1175         NSIN=NSIN+1
1176       ENDIF
1177   190 CONTINUE
1178       IF(PDM.GE.PARJ(32)) GOTO 320
1179
1180 C...Fill small-mass system as cluster.
1181       NSAV=N
1182       PECM=SQRT(MAX(0D0,DPC(4)**2-DPC(1)**2-DPC(2)**2-DPC(3)**2))
1183       K(N+1,1)=11
1184       K(N+1,2)=91
1185       K(N+1,3)=IC1
1186       K(N+1,4)=N+2
1187       K(N+1,5)=N+3
1188       P(N+1,1)=DPC(1)
1189       P(N+1,2)=DPC(2)
1190       P(N+1,3)=DPC(3)
1191       P(N+1,4)=DPC(4)
1192       P(N+1,5)=PECM
1193
1194 C...Form two particles from flavours of lowest-mass system, if feasible.
1195       K(N+2,1)=1
1196       K(N+3,1)=1
1197       IF(MSTU(16).NE.2) THEN
1198         K(N+2,3)=N+1
1199         K(N+3,3)=N+1
1200       ELSE
1201         K(N+2,3)=IC1
1202         K(N+3,3)=IC2
1203       ENDIF
1204       K(N+2,4)=0
1205       K(N+3,4)=0
1206       K(N+2,5)=0
1207       K(N+3,5)=0
1208       IF(IABS(K(IC1,2)).NE.21) THEN
1209         KC1=LUCOMP(K(IC1,2))
1210         KC2=LUCOMP(K(IC2,2))
1211         IF(KC1.EQ.0.OR.KC2.EQ.0) GOTO 320
1212         KQ1=KCHG(KC1,2)*ISIGN(1,K(IC1,2))
1213         KQ2=KCHG(KC2,2)*ISIGN(1,K(IC2,2))
1214         IF(KQ1+KQ2.NE.0) GOTO 320
1215   200   CALL LUKFDI(K(IC1,2),0,KFLN,K(N+2,2))
1216         CALL LUKFDI(K(IC2,2),-KFLN,KFLDMP,K(N+3,2))
1217         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 200
1218       ELSE
1219         IF(IABS(K(IC2,2)).NE.21) GOTO 320
1220   210   CALL LUKFDI(1+INT((2.+PARJ(2))*RLU(0)),0,KFLN,KFDMP)
1221         CALL LUKFDI(KFLN,0,KFLM,K(N+2,2))
1222         CALL LUKFDI(-KFLN,-KFLM,KFLDMP,K(N+3,2))
1223         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 210
1224       ENDIF
1225       P(N+2,5)=ULMASS(K(N+2,2))
1226       P(N+3,5)=ULMASS(K(N+3,2))
1227       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM.AND.NSIN.EQ.1) GOTO 320
1228       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) GOTO 260
1229
1230 C...Perform two-particle decay of jet system, if possible.
1231       IF(PECM.GE.0.02*DPC(4)) THEN
1232         PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-
1233      &  (P(N+2,5)-P(N+3,5))**2))/(2.*PECM)
1234         UE(3)=2.*RLU(0)-1.
1235         PHI=PARU(2)*RLU(0)
1236         UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
1237         UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
1238         DO 220 J=1,3
1239         P(N+2,J)=PA*UE(J)
1240   220   P(N+3,J)=-PA*UE(J)
1241         P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)
1242         P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)
1243         MSTU(33)=1
1244         CALL LUDBRB(N+2,N+3,0.,0.,DPC(1)/DPC(4),DPC(2)/DPC(4),
1245      &  DPC(3)/DPC(4))
1246       ELSE
1247         NP=0
1248         DO 230 I=IC1,IC2
1249   230   IF(K(I,1).EQ.1.OR.K(I,1).EQ.2) NP=NP+1
1250         HA=P(IC1,4)*P(IC2,4)-P(IC1,1)*P(IC2,1)-P(IC1,2)*P(IC2,2)-
1251      &  P(IC1,3)*P(IC2,3)
1252         IF(NP.GE.3.OR.HA.LE.1.25*P(IC1,5)*P(IC2,5)) GOTO 260
1253         HD1=0.5*(P(N+2,5)**2-P(IC1,5)**2)
1254         HD2=0.5*(P(N+3,5)**2-P(IC2,5)**2)
1255         HR=SQRT(MAX(0.,((HA-HD1-HD2)**2-(P(N+2,5)*P(N+3,5))**2)/
1256      &  (HA**2-(P(IC1,5)*P(IC2,5))**2)))-1.
1257         HC=P(IC1,5)**2+2.*HA+P(IC2,5)**2
1258         HK1=((P(IC2,5)**2+HA)*HR+HD1-HD2)/HC
1259         HK2=((P(IC1,5)**2+HA)*HR+HD2-HD1)/HC
1260         DO 240 J=1,4
1261         P(N+2,J)=(1.+HK1)*P(IC1,J)-HK2*P(IC2,J)
1262   240   P(N+3,J)=(1.+HK2)*P(IC2,J)-HK1*P(IC1,J)
1263       ENDIF
1264       DO 250 J=1,4
1265       V(N+1,J)=V(IC1,J)
1266       V(N+2,J)=V(IC1,J)
1267   250 V(N+3,J)=V(IC2,J)
1268       V(N+1,5)=0.
1269       V(N+2,5)=0.
1270       V(N+3,5)=0.
1271       N=N+3
1272       GOTO 300
1273
1274 C...Else form one particle from the flavours available, if possible.
1275   260 K(N+1,5)=N+2
1276       IF(IABS(K(IC1,2)).GT.100.AND.IABS(K(IC2,2)).GT.100) THEN
1277         GOTO 320
1278       ELSEIF(IABS(K(IC1,2)).NE.21) THEN
1279         CALL LUKFDI(K(IC1,2),K(IC2,2),KFLDMP,K(N+2,2))
1280       ELSE
1281         KFLN=1+INT((2.+PARJ(2))*RLU(0))
1282         CALL LUKFDI(KFLN,-KFLN,KFLDMP,K(N+2,2))
1283       ENDIF
1284       IF(K(N+2,2).EQ.0) GOTO 260
1285       P(N+2,5)=ULMASS(K(N+2,2))
1286
1287 C...Find parton/particle which combines to largest extra mass.
1288       IR=0
1289       HA=0.
1290       HSM=0.
1291       DO 280 MCOMB=1,3
1292       IF(IR.NE.0) GOTO 280
1293       DO 270 I=MAX(1,IP),N
1294       IF(K(I,1).LE.0.OR.K(I,1).GT.10.OR.(I.GE.IC1.AND.I.LE.IC2.
1295      &AND.K(I,1).GE.1.AND.K(I,1).LE.2)) GOTO 270
1296       IF(MCOMB.EQ.1) KCI=LUCOMP(K(I,2))
1297       IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 270
1298       IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 270
1299       IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100)
1300      &GOTO 270
1301       HCR=DPC(4)*P(I,4)-DPC(1)*P(I,1)-DPC(2)*P(I,2)-DPC(3)*P(I,3)
1302       HSR=2.*HCR+PECM**2-P(N+2,5)**2-2.*P(N+2,5)*P(I,5)
1303       IF(HSR.GT.HSM) THEN
1304         IR=I
1305         HA=HCR
1306         HSM=HSR
1307       ENDIF
1308   270 CONTINUE
1309   280 CONTINUE
1310
1311 C...Shuffle energy and momentum to put new particle on mass shell.
1312       IF(IR.NE.0) THEN
1313         HB=PECM**2+HA
1314         HC=P(N+2,5)**2+HA
1315         HD=P(IR,5)**2+HA
1316         HK2=0.5*(HB*SQRT(MAX(0.,((HB+HC)**2-4.*(HB+HD)*P(N+2,5)**2)/
1317      &  (HA**2-(PECM*P(IR,5))**2)))-(HB+HC))/(HB+HD)
1318         HK1=(0.5*(P(N+2,5)**2-PECM**2)+HD*HK2)/HB
1319         DO 290 J=1,4
1320         P(N+2,J)=(1.+HK1)*DPC(J)-HK2*P(IR,J)
1321         P(IR,J)=(1.+HK2)*P(IR,J)-HK1*DPC(J)
1322         V(N+1,J)=V(IC1,J)
1323   290   V(N+2,J)=V(IC1,J)
1324         V(N+1,5)=0.
1325         V(N+2,5)=0.
1326         N=N+2
1327       ELSE
1328         CALL LUERRM(3,'(LUPREP:) no match for collapsing cluster')
1329         RETURN
1330       ENDIF
1331
1332 C...Mark collapsed system and store daughter pointers. Iterate.
1333   300 DO 310 I=IC1,IC2
1334       IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.KCHG(LUCOMP(K(I,2)),2).NE.0)
1335      &THEN
1336         K(I,1)=K(I,1)+10
1337         IF(MSTU(16).NE.2) THEN
1338           K(I,4)=NSAV+1
1339           K(I,5)=NSAV+1
1340         ELSE
1341           K(I,4)=NSAV+2
1342           K(I,5)=N
1343         ENDIF
1344       ENDIF
1345   310 CONTINUE
1346       IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 140
1347
1348 C...Check flavours and invariant masses in parton systems.
1349   320 NP=0
1350       KFN=0
1351       KQS=0
1352       DO 330 J=1,5
1353   330 DPS(J)=0.
1354       DO 360 I=MAX(1,IP),N
1355       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360
1356       KC=LUCOMP(K(I,2))
1357       IF(KC.EQ.0) GOTO 360
1358       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
1359       IF(KQ.EQ.0) GOTO 360
1360       NP=NP+1
1361       IF(KQ.NE.2) THEN
1362         KFN=KFN+1
1363         KQS=KQS+KQ
1364         MSTJ(93)=1
1365         DPS(5)=DPS(5)+ULMASS(K(I,2))
1366       ENDIF
1367       DO 340 J=1,4
1368   340 DPS(J)=DPS(J)+P(I,J)
1369       IF(K(I,1).EQ.1) THEN
1370         IF(NP.NE.1.AND.(KFN.EQ.1.OR.KFN.GE.3.OR.KQS.NE.0)) CALL
1371      &  LUERRM(2,'(LUPREP:) unphysical flavour combination')
1372         IF(NP.NE.1.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.
1373      &  (0.9*PARJ(32)+DPS(5))**2) CALL LUERRM(3,
1374      &  '(LUPREP:) too small mass in jet system')
1375         NP=0
1376         KFN=0
1377         KQS=0
1378         DO 350 J=1,5
1379   350   DPS(J)=0.
1380       ENDIF
1381   360 CONTINUE
1382
1383       RETURN
1384       END
1385
1386 C*********************************************************************
1387
1388       SUBROUTINE LUSTRF(IP)
1389 C...Purpose: to handle the fragmentation of an arbitrary colour singlet
1390 C...jet system according to the Lund string fragmentation model.
1391       IMPLICIT DOUBLE PRECISION(D)
1392       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
1393       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1394       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
1395       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
1396       DIMENSION DPS(5),KFL(3),PMQ(3),PX(3),PY(3),GAM(3),IE(2),PR(2),
1397      &IN(9),DHM(4),DHG(4),DP(5,5),IRANK(2),MJU(4),IJU(3),PJU(5,5),
1398      &TJU(5),KFJH(2),NJS(2),KFJS(2),PJS(4,5),MSTU9T(8),PARU9T(8)
1399
1400 C...Function: four-product of two vectors.
1401       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)
1402       DFOUR(I,J)=DP(I,4)*DP(J,4)-DP(I,1)*DP(J,1)-DP(I,2)*DP(J,2)-
1403      &DP(I,3)*DP(J,3)
1404
1405 C...Reset counters. Identify parton system.
1406       MSTJ(91)=0
1407       NSAV=N
1408       MSTU90=MSTU(90)
1409       NP=0
1410       KQSUM=0
1411       DO 100 J=1,5
1412   100 DPS(J)=0D0
1413       MJU(1)=0
1414       MJU(2)=0
1415       I=IP-1
1416   110 I=I+1
1417       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
1418         CALL LUERRM(12,'(LUSTRF:) failed to reconstruct jet system')
1419         IF(MSTU(21).GE.1) RETURN
1420       ENDIF
1421       IF(K(I,1).NE.1.AND.K(I,1).NE.2.AND.K(I,1).NE.41) GOTO 110
1422       KC=LUCOMP(K(I,2))
1423       IF(KC.EQ.0) GOTO 110
1424       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
1425       IF(KQ.EQ.0) GOTO 110
1426       IF(N+5*NP+11.GT.MSTU(4)-MSTU(32)-5) THEN
1427         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS')
1428         IF(MSTU(21).GE.1) RETURN
1429       ENDIF
1430
1431 C...Take copy of partons to be considered. Check flavour sum.
1432       NP=NP+1
1433       DO 120 J=1,5
1434       K(N+NP,J)=K(I,J)
1435       P(N+NP,J)=P(I,J)
1436   120 IF(J.NE.4) DPS(J)=DPS(J)+P(I,J)
1437       DPS(4)=DPS(4)+SQRT(DBLE(P(I,1))**2+DBLE(P(I,2))**2+
1438      &DBLE(P(I,3))**2+DBLE(P(I,5))**2)
1439       K(N+NP,3)=I
1440       IF(KQ.NE.2) KQSUM=KQSUM+KQ
1441       IF(K(I,1).EQ.41) THEN
1442         KQSUM=KQSUM+2*KQ
1443         IF(KQSUM.EQ.KQ) MJU(1)=N+NP
1444         IF(KQSUM.NE.KQ) MJU(2)=N+NP
1445       ENDIF
1446       IF(K(I,1).EQ.2.OR.K(I,1).EQ.41) GOTO 110
1447       IF(KQSUM.NE.0) THEN
1448         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination')
1449         IF(MSTU(21).GE.1) RETURN
1450       ENDIF
1451
1452 C...Boost copied system to CM frame (for better numerical precision).
1453       IF(ABS(DPS(3)).LT.0.99D0*DPS(4)) THEN
1454         MBST=0
1455         MSTU(33)=1
1456         CALL LUDBRB(N+1,N+NP,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4),
1457      &  -DPS(3)/DPS(4))
1458       ELSE
1459         MBST=1
1460         HHBZ=SQRT(MAX(1D-6,DPS(4)+DPS(3))/MAX(1D-6,DPS(4)-DPS(3)))
1461         DO 130 I=N+1,N+NP
1462         HHPMT=P(I,1)**2+P(I,2)**2+P(I,5)**2
1463         IF(P(I,3).GT.0.) THEN
1464           HHPEZ=(P(I,4)+P(I,3))/HHBZ
1465           P(I,3)=0.5*(HHPEZ-HHPMT/HHPEZ)
1466           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ)
1467         ELSE
1468           HHPEZ=(P(I,4)-P(I,3))*HHBZ
1469           P(I,3)=-0.5*(HHPEZ-HHPMT/HHPEZ)
1470           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ)
1471         ENDIF
1472   130   CONTINUE
1473       ENDIF
1474
1475 C...Search for very nearby partons that may be recombined.
1476       NTRYR=0
1477       PARU12=PARU(12)
1478       PARU13=PARU(13)
1479       MJU(3)=MJU(1)
1480       MJU(4)=MJU(2)
1481       NR=NP
1482   140 IF(NR.GE.3) THEN
1483         PDRMIN=2.*PARU12
1484         DO 150 I=N+1,N+NR
1485         IF(I.EQ.N+NR.AND.IABS(K(N+1,2)).NE.21) GOTO 150
1486         I1=I+1
1487         IF(I.EQ.N+NR) I1=N+1
1488         IF(K(I,1).EQ.41.OR.K(I1,1).EQ.41) GOTO 150
1489         IF(MJU(1).NE.0.AND.I1.LT.MJU(1).AND.IABS(K(I1,2)).NE.21)
1490      &  GOTO 150
1491         IF(MJU(2).NE.0.AND.I.GT.MJU(2).AND.IABS(K(I,2)).NE.21) GOTO 150
1492         PAP=SQRT((P(I,1)**2+P(I,2)**2+P(I,3)**2)*(P(I1,1)**2+
1493      &  P(I1,2)**2+P(I1,3)**2))
1494         PVP=P(I,1)*P(I1,1)+P(I,2)*P(I1,2)+P(I,3)*P(I1,3)
1495         PDR=4.*(PAP-PVP)**2/MAX(1E-6,PARU13**2*PAP+2.*(PAP-PVP))
1496         IF(PDR.LT.PDRMIN) THEN
1497           IR=I
1498           PDRMIN=PDR
1499         ENDIF
1500   150   CONTINUE
1501
1502 C...Recombine very nearby partons to avoid machine precision problems.
1503         IF(PDRMIN.LT.PARU12.AND.IR.EQ.N+NR) THEN
1504           DO 160 J=1,4
1505   160     P(N+1,J)=P(N+1,J)+P(N+NR,J)
1506           P(N+1,5)=SQRT(MAX(0.,P(N+1,4)**2-P(N+1,1)**2-P(N+1,2)**2-
1507      &    P(N+1,3)**2))
1508           NR=NR-1
1509           GOTO 140
1510         ELSEIF(PDRMIN.LT.PARU12) THEN
1511           DO 170 J=1,4
1512   170     P(IR,J)=P(IR,J)+P(IR+1,J)
1513           P(IR,5)=SQRT(MAX(0.,P(IR,4)**2-P(IR,1)**2-P(IR,2)**2-
1514      &    P(IR,3)**2))
1515           DO 180 I=IR+1,N+NR-1
1516           K(I,2)=K(I+1,2)
1517           DO 180 J=1,5
1518   180     P(I,J)=P(I+1,J)
1519           IF(IR.EQ.N+NR-1) K(IR,2)=K(N+NR,2)
1520           NR=NR-1
1521           IF(MJU(1).GT.IR) MJU(1)=MJU(1)-1
1522           IF(MJU(2).GT.IR) MJU(2)=MJU(2)-1
1523           GOTO 140
1524         ENDIF
1525       ENDIF
1526       NTRYR=NTRYR+1
1527
1528 C...Reset particle counter. Skip ahead if no junctions are present;
1529 C...this is usually the case!
1530       NRS=MAX(5*NR+11,NP)
1531       NTRY=0
1532   190 NTRY=NTRY+1
1533       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
1534         PARU12=4.*PARU12
1535         PARU13=2.*PARU13
1536         GOTO 140
1537       ELSEIF(NTRY.GT.100) THEN
1538         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop')
1539         IF(MSTU(21).GE.1) RETURN
1540       ENDIF
1541       I=N+NRS
1542       MSTU(90)=MSTU90
1543       IF(MJU(1).EQ.0.AND.MJU(2).EQ.0) GOTO 510
1544       DO 500 JT=1,2
1545       NJS(JT)=0
1546       IF(MJU(JT).EQ.0) GOTO 500
1547       JS=3-2*JT
1548
1549 C...Find and sum up momentum on three sides of junction. Check flavours.
1550       DO 200 IU=1,3
1551       IJU(IU)=0
1552       DO 200 J=1,5
1553   200 PJU(IU,J)=0.
1554       IU=0
1555       DO 210 I1=N+1+(JT-1)*(NR-1),N+NR+(JT-1)*(1-NR),JS
1556       IF(K(I1,2).NE.21.AND.IU.LE.2) THEN
1557         IU=IU+1
1558         IJU(IU)=I1
1559       ENDIF
1560       DO 210 J=1,4
1561   210 PJU(IU,J)=PJU(IU,J)+P(I1,J)
1562       DO 220 IU=1,3
1563   220 PJU(IU,5)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)
1564       IF(K(IJU(3),2)/100.NE.10*K(IJU(1),2)+K(IJU(2),2).AND.
1565      &K(IJU(3),2)/100.NE.10*K(IJU(2),2)+K(IJU(1),2)) THEN
1566         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination')
1567         IF(MSTU(21).GE.1) RETURN
1568       ENDIF
1569
1570 C...Calculate (approximate) boost to rest frame of junction.
1571       T12=(PJU(1,1)*PJU(2,1)+PJU(1,2)*PJU(2,2)+PJU(1,3)*PJU(2,3))/
1572      &(PJU(1,5)*PJU(2,5))
1573       T13=(PJU(1,1)*PJU(3,1)+PJU(1,2)*PJU(3,2)+PJU(1,3)*PJU(3,3))/
1574      &(PJU(1,5)*PJU(3,5))
1575       T23=(PJU(2,1)*PJU(3,1)+PJU(2,2)*PJU(3,2)+PJU(2,3)*PJU(3,3))/
1576      &(PJU(2,5)*PJU(3,5))
1577       T11=SQRT((2./3.)*(1.-T12)*(1.-T13)/(1.-T23))
1578       T22=SQRT((2./3.)*(1.-T12)*(1.-T23)/(1.-T13))
1579       TSQ=SQRT((2.*T11*T22+T12-1.)*(1.+T12))
1580       T1F=(TSQ-T22*(1.+T12))/(1.-T12**2)
1581       T2F=(TSQ-T11*(1.+T12))/(1.-T12**2)
1582       DO 230 J=1,3
1583   230 TJU(J)=-(T1F*PJU(1,J)/PJU(1,5)+T2F*PJU(2,J)/PJU(2,5))
1584       TJU(4)=SQRT(1.+TJU(1)**2+TJU(2)**2+TJU(3)**2)
1585       DO 240 IU=1,3
1586   240 PJU(IU,5)=TJU(4)*PJU(IU,4)-TJU(1)*PJU(IU,1)-TJU(2)*PJU(IU,2)-
1587      &TJU(3)*PJU(IU,3)
1588
1589 C...Put junction at rest if motion could give inconsistencies.
1590       IF(PJU(1,5)+PJU(2,5).GT.PJU(1,4)+PJU(2,4)) THEN
1591         DO 250 J=1,3
1592   250   TJU(J)=0.
1593         TJU(4)=1.
1594         PJU(1,5)=PJU(1,4)
1595         PJU(2,5)=PJU(2,4)
1596         PJU(3,5)=PJU(3,4)
1597       ENDIF
1598
1599 C...Start preparing for fragmentation of two strings from junction.
1600       ISTA=I
1601       DO 480 IU=1,2
1602       NS=IJU(IU+1)-IJU(IU)
1603
1604 C...Junction strings: find longitudinal string directions.
1605       DO 270 IS=1,NS
1606       IS1=IJU(IU)+IS-1
1607       IS2=IJU(IU)+IS
1608       DO 260 J=1,5
1609       DP(1,J)=0.5*P(IS1,J)
1610       IF(IS.EQ.1) DP(1,J)=P(IS1,J)
1611       DP(2,J)=0.5*P(IS2,J)
1612   260 IF(IS.EQ.NS) DP(2,J)=-PJU(IU,J)
1613       IF(IS.EQ.NS) DP(2,4)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2)
1614       IF(IS.EQ.NS) DP(2,5)=0.
1615       DP(3,5)=DFOUR(1,1)
1616       DP(4,5)=DFOUR(2,2)
1617       DHKC=DFOUR(1,2)
1618       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN
1619         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
1620         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
1621         DP(3,5)=0D0
1622         DP(4,5)=0D0
1623         DHKC=DFOUR(1,2)
1624       ENDIF
1625       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
1626       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.)
1627       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.)
1628       IN1=N+NR+4*IS-3
1629       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5))
1630       DO 270 J=1,4
1631       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J)
1632   270 P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J)
1633
1634 C...Junction strings: initialize flavour, momentum and starting pos.
1635       ISAV=I
1636       MSTU91=MSTU(90)
1637   280 NTRY=NTRY+1
1638       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
1639         PARU12=4.*PARU12
1640         PARU13=2.*PARU13
1641         GOTO 140
1642       ELSEIF(NTRY.GT.100) THEN
1643         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop')
1644         IF(MSTU(21).GE.1) RETURN
1645       ENDIF
1646       I=ISAV
1647       MSTU(90)=MSTU91
1648       IRANKJ=0
1649       IE(1)=K(N+1+(JT/2)*(NP-1),3)
1650       IN(4)=N+NR+1
1651       IN(5)=IN(4)+1
1652       IN(6)=N+NR+4*NS+1
1653       DO 290 JQ=1,2
1654       DO 290 IN1=N+NR+2+JQ,N+NR+4*NS-2+JQ,4
1655       P(IN1,1)=2-JQ
1656       P(IN1,2)=JQ-1
1657   290 P(IN1,3)=1.
1658       KFL(1)=K(IJU(IU),2)
1659       PX(1)=0.
1660       PY(1)=0.
1661       GAM(1)=0.
1662       DO 300 J=1,5
1663   300 PJU(IU+3,J)=0.
1664
1665 C...Junction strings: find initial transverse directions.
1666       DO 310 J=1,4
1667       DP(1,J)=P(IN(4),J)
1668       DP(2,J)=P(IN(4)+1,J)
1669       DP(3,J)=0.
1670   310 DP(4,J)=0.
1671       DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
1672       DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
1673       DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
1674       DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
1675       DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
1676       IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
1677       IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
1678       IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
1679       IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
1680       DHC12=DFOUR(1,2)
1681       DHCX1=DFOUR(3,1)/DHC12
1682       DHCX2=DFOUR(3,2)/DHC12
1683       DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
1684       DHCY1=DFOUR(4,1)/DHC12
1685       DHCY2=DFOUR(4,2)/DHC12
1686       DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
1687       DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
1688       DO 320 J=1,4
1689       DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
1690       P(IN(6),J)=DP(3,J)
1691   320 P(IN(6)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
1692      &DHCYX*DP(3,J))
1693
1694 C...Junction strings: produce new particle, origin.
1695   330 I=I+1
1696       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
1697         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS')
1698         IF(MSTU(21).GE.1) RETURN
1699       ENDIF
1700       IRANKJ=IRANKJ+1
1701       K(I,1)=1
1702       K(I,3)=IE(1)
1703       K(I,4)=0
1704       K(I,5)=0
1705
1706 C...Junction strings: generate flavour, hadron, pT, z and Gamma.
1707   340 CALL LUKFDI(KFL(1),0,KFL(3),K(I,2))
1708       IF(K(I,2).EQ.0) GOTO 280
1709       IF(MSTJ(12).GE.3.AND.IRANKJ.EQ.1.AND.IABS(KFL(1)).LE.10.AND.
1710      &IABS(KFL(3)).GT.10) THEN
1711         IF(RLU(0).GT.PARJ(19)) GOTO 340
1712       ENDIF
1713       P(I,5)=ULMASS(K(I,2))
1714       CALL LUPTDI(KFL(1),PX(3),PY(3))
1715       PR(1)=P(I,5)**2+(PX(1)+PX(3))**2+(PY(1)+PY(3))**2
1716       CALL LUZDIS(KFL(1),KFL(3),PR(1),Z)
1717       IF(IABS(KFL(1)).GE.4.AND.IABS(KFL(1)).LE.8.AND.
1718      &MSTU(90).LT.8) THEN
1719         MSTU(90)=MSTU(90)+1
1720         MSTU(90+MSTU(90))=I
1721         PARU(90+MSTU(90))=Z
1722       ENDIF
1723       GAM(3)=(1.-Z)*(GAM(1)+PR(1)/Z)
1724       DO 350 J=1,3
1725   350 IN(J)=IN(3+J)
1726
1727 C...Junction strings: stepping within or from 'low' string region easy.
1728       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
1729      &P(IN(1),5)**2.GE.PR(1)) THEN
1730         P(IN(1)+2,4)=Z*P(IN(1)+2,3)
1731         P(IN(2)+2,4)=PR(1)/(P(IN(1)+2,4)*P(IN(1),5)**2)
1732         DO 360 J=1,4
1733   360   P(I,J)=(PX(1)+PX(3))*P(IN(3),J)+(PY(1)+PY(3))*P(IN(3)+1,J)
1734         GOTO 430
1735       ELSEIF(IN(1)+1.EQ.IN(2)) THEN
1736         P(IN(2)+2,4)=P(IN(2)+2,3)
1737         P(IN(2)+2,1)=1.
1738         IN(2)=IN(2)+4
1739         IF(IN(2).GT.N+NR+4*NS) GOTO 280
1740         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
1741           P(IN(1)+2,4)=P(IN(1)+2,3)
1742           P(IN(1)+2,1)=0.
1743           IN(1)=IN(1)+4
1744         ENDIF
1745       ENDIF
1746
1747 C...Junction strings: find new transverse directions.
1748   370 IF(IN(1).GT.N+NR+4*NS.OR.IN(2).GT.N+NR+4*NS.OR.
1749      &IN(1).GT.IN(2)) GOTO 280
1750       IF(IN(1).NE.IN(4).OR.IN(2).NE.IN(5)) THEN
1751         DO 380 J=1,4
1752         DP(1,J)=P(IN(1),J)
1753         DP(2,J)=P(IN(2),J)
1754         DP(3,J)=0.
1755   380   DP(4,J)=0.
1756         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
1757         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
1758         DHC12=DFOUR(1,2)
1759         IF(DHC12.LE.1E-2) THEN
1760           P(IN(1)+2,4)=P(IN(1)+2,3)
1761           P(IN(1)+2,1)=0.
1762           IN(1)=IN(1)+4
1763           GOTO 370
1764         ENDIF
1765         IN(3)=N+NR+4*NS+5
1766         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
1767         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
1768         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
1769         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
1770         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
1771         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
1772         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
1773         DHCX1=DFOUR(3,1)/DHC12
1774         DHCX2=DFOUR(3,2)/DHC12
1775         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
1776         DHCY1=DFOUR(4,1)/DHC12
1777         DHCY2=DFOUR(4,2)/DHC12
1778         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
1779         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
1780         DO 390 J=1,4
1781         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
1782         P(IN(3),J)=DP(3,J)
1783   390   P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
1784      &  DHCYX*DP(3,J))
1785 C...Express pT with respect to new axes, if sensible.
1786         PXP=-(PX(3)*FOUR(IN(6),IN(3))+PY(3)*FOUR(IN(6)+1,IN(3)))
1787         PYP=-(PX(3)*FOUR(IN(6),IN(3)+1)+PY(3)*FOUR(IN(6)+1,IN(3)+1))
1788         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN
1789           PX(3)=PXP
1790           PY(3)=PYP
1791         ENDIF
1792       ENDIF
1793
1794 C...Junction strings: sum up known four-momentum, coefficients for m2.
1795       DO 410 J=1,4
1796       DHG(J)=0.
1797       P(I,J)=PX(1)*P(IN(6),J)+PY(1)*P(IN(6)+1,J)+PX(3)*P(IN(3),J)+
1798      &PY(3)*P(IN(3)+1,J)
1799       DO 400 IN1=IN(4),IN(1)-4,4
1800   400 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
1801       DO 410 IN2=IN(5),IN(2)-4,4
1802   410 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
1803       DHM(1)=FOUR(I,I)
1804       DHM(2)=2.*FOUR(I,IN(1))
1805       DHM(3)=2.*FOUR(I,IN(2))
1806       DHM(4)=2.*FOUR(IN(1),IN(2))
1807
1808 C...Junction strings: find coefficients for Gamma expression.
1809       DO 420 IN2=IN(1)+1,IN(2),4
1810       DO 420 IN1=IN(1),IN2-1,4
1811       DHC=2.*FOUR(IN1,IN2)
1812       DHG(1)=DHG(1)+P(IN1+2,1)*P(IN2+2,1)*DHC
1813       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-P(IN2+2,1)*DHC
1814       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+P(IN1+2,1)*DHC
1815   420 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
1816
1817 C...Junction strings: solve (m2, Gamma) equation system for energies.
1818       DHS1=DHM(3)*DHG(4)-DHM(4)*DHG(3)
1819       IF(ABS(DHS1).LT.1E-4) GOTO 280
1820       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(2)*DHG(3)-DHG(4)*
1821      &(P(I,5)**2-DHM(1))+DHG(2)*DHM(3)
1822       DHS3=DHM(2)*(GAM(3)-DHG(1))-DHG(2)*(P(I,5)**2-DHM(1))
1823       P(IN(2)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)-
1824      &DHS2/DHS1)
1825       IF(DHM(2)+DHM(4)*P(IN(2)+2,4).LE.0.) GOTO 280
1826       P(IN(1)+2,4)=(P(I,5)**2-DHM(1)-DHM(3)*P(IN(2)+2,4))/
1827      &(DHM(2)+DHM(4)*P(IN(2)+2,4))
1828
1829 C...Junction strings: step to new region if necessary.
1830       IF(P(IN(2)+2,4).GT.P(IN(2)+2,3)) THEN
1831         P(IN(2)+2,4)=P(IN(2)+2,3)
1832         P(IN(2)+2,1)=1.
1833         IN(2)=IN(2)+4
1834         IF(IN(2).GT.N+NR+4*NS) GOTO 280
1835         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
1836           P(IN(1)+2,4)=P(IN(1)+2,3)
1837           P(IN(1)+2,1)=0.
1838           IN(1)=IN(1)+4
1839         ENDIF
1840         GOTO 370
1841       ELSEIF(P(IN(1)+2,4).GT.P(IN(1)+2,3)) THEN
1842         P(IN(1)+2,4)=P(IN(1)+2,3)
1843         P(IN(1)+2,1)=0.
1844         IN(1)=IN(1)+JS
1845         GOTO 720
1846       ENDIF
1847
1848 C...Junction strings: particle four-momentum, remainder, loop back.
1849   430 DO 440 J=1,4
1850       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J)
1851   440 PJU(IU+3,J)=PJU(IU+3,J)+P(I,J)
1852       IF(P(I,4).LT.P(I,5)) GOTO 280
1853       PJU(IU+3,5)=TJU(4)*PJU(IU+3,4)-TJU(1)*PJU(IU+3,1)-
1854      &TJU(2)*PJU(IU+3,2)-TJU(3)*PJU(IU+3,3)
1855       IF(PJU(IU+3,5).LT.PJU(IU,5)) THEN
1856         KFL(1)=-KFL(3)
1857         PX(1)=-PX(3)
1858         PY(1)=-PY(3)
1859         GAM(1)=GAM(3)
1860         IF(IN(3).NE.IN(6)) THEN
1861           DO 450 J=1,4
1862           P(IN(6),J)=P(IN(3),J)
1863   450     P(IN(6)+1,J)=P(IN(3)+1,J)
1864         ENDIF
1865         DO 460 JQ=1,2
1866         IN(3+JQ)=IN(JQ)
1867         P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
1868   460   P(IN(JQ)+2,1)=P(IN(JQ)+2,1)-(3-2*JQ)*P(IN(JQ)+2,4)
1869         GOTO 330
1870       ENDIF
1871
1872 C...Junction strings: save quantities left after each string.
1873       IF(IABS(KFL(1)).GT.10) GOTO 280
1874       I=I-1
1875       KFJH(IU)=KFL(1)
1876       DO 470 J=1,4
1877   470 PJU(IU+3,J)=PJU(IU+3,J)-P(I+1,J)
1878   480 CONTINUE
1879
1880 C...Junction strings: put together to new effective string endpoint.
1881       NJS(JT)=I-ISTA
1882       KFJS(JT)=K(K(MJU(JT+2),3),2)
1883       KFLS=2*INT(RLU(0)+3.*PARJ(4)/(1.+3.*PARJ(4)))+1
1884       IF(KFJH(1).EQ.KFJH(2)) KFLS=3
1885       IF(ISTA.NE.I) KFJS(JT)=ISIGN(1000*MAX(IABS(KFJH(1)),
1886      &IABS(KFJH(2)))+100*MIN(IABS(KFJH(1)),IABS(KFJH(2)))+
1887      &KFLS,KFJH(1))
1888       DO 490 J=1,4
1889       PJS(JT,J)=PJU(1,J)+PJU(2,J)+P(MJU(JT),J)
1890   490 PJS(JT+2,J)=PJU(4,J)+PJU(5,J)
1891       PJS(JT,5)=SQRT(MAX(0.,PJS(JT,4)**2-PJS(JT,1)**2-PJS(JT,2)**2-
1892      &PJS(JT,3)**2))
1893   500 CONTINUE
1894
1895 C...Open versus closed strings. Choose breakup region for latter.
1896   510 IF(MJU(1).NE.0.AND.MJU(2).NE.0) THEN
1897         NS=MJU(2)-MJU(1)
1898         NB=MJU(1)-N
1899       ELSEIF(MJU(1).NE.0) THEN
1900         NS=N+NR-MJU(1)
1901         NB=MJU(1)-N
1902       ELSEIF(MJU(2).NE.0) THEN
1903         NS=MJU(2)-N
1904         NB=1
1905       ELSEIF(IABS(K(N+1,2)).NE.21) THEN
1906         NS=NR-1
1907         NB=1
1908       ELSE
1909         NS=NR+1
1910         W2SUM=0.
1911         DO 520 IS=1,NR
1912         P(N+NR+IS,1)=0.5*FOUR(N+IS,N+IS+1-NR*(IS/NR))
1913   520   W2SUM=W2SUM+P(N+NR+IS,1)
1914         W2RAN=RLU(0)*W2SUM
1915         NB=0
1916   530   NB=NB+1
1917         W2SUM=W2SUM-P(N+NR+NB,1)
1918         IF(W2SUM.GT.W2RAN.AND.NB.LT.NR) GOTO 530
1919       ENDIF
1920
1921 C...Find longitudinal string directions (i.e. lightlike four-vectors).
1922       DO 550 IS=1,NS
1923       IS1=N+IS+NB-1-NR*((IS+NB-2)/NR)
1924       IS2=N+IS+NB-NR*((IS+NB-1)/NR)
1925       DO 540 J=1,5
1926       DP(1,J)=P(IS1,J)
1927       IF(IABS(K(IS1,2)).EQ.21) DP(1,J)=0.5*DP(1,J)
1928       IF(IS1.EQ.MJU(1)) DP(1,J)=PJS(1,J)-PJS(3,J)
1929       DP(2,J)=P(IS2,J)
1930       IF(IABS(K(IS2,2)).EQ.21) DP(2,J)=0.5*DP(2,J)
1931   540 IF(IS2.EQ.MJU(2)) DP(2,J)=PJS(2,J)-PJS(4,J)
1932       DP(3,5)=DFOUR(1,1)
1933       DP(4,5)=DFOUR(2,2)
1934       DHKC=DFOUR(1,2)
1935       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN
1936         DP(3,5)=DP(1,5)**2
1937         DP(4,5)=DP(2,5)**2
1938         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2+DP(1,5)**2)
1939         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2+DP(2,5)**2)
1940         DHKC=DFOUR(1,2)
1941       ENDIF
1942       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5))
1943       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.)
1944       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.)
1945       IN1=N+NR+4*IS-3
1946       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5))
1947       DO 550 J=1,4
1948       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J)
1949   550 P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J)
1950
1951 C...Begin initialization: sum up energy, set starting position.
1952       ISAV=I
1953       MSTU91=MSTU(90)
1954   560 NTRY=NTRY+1
1955       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN
1956         PARU12=4.*PARU12
1957         PARU13=2.*PARU13
1958         GOTO 140
1959       ELSEIF(NTRY.GT.100) THEN
1960         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop')
1961         IF(MSTU(21).GE.1) RETURN
1962       ENDIF
1963       I=ISAV
1964       MSTU(90)=MSTU91
1965       DO 570 J=1,4
1966       P(N+NRS,J)=0.
1967       DO 570 IS=1,NR
1968   570 P(N+NRS,J)=P(N+NRS,J)+P(N+IS,J)
1969       DO 580 JT=1,2
1970       IRANK(JT)=0
1971       IF(MJU(JT).NE.0) IRANK(JT)=NJS(JT)
1972       IF(NS.GT.NR) IRANK(JT)=1
1973       IE(JT)=K(N+1+(JT/2)*(NP-1),3)
1974       IN(3*JT+1)=N+NR+1+4*(JT/2)*(NS-1)
1975       IN(3*JT+2)=IN(3*JT+1)+1
1976       IN(3*JT+3)=N+NR+4*NS+2*JT-1
1977       DO 580 IN1=N+NR+2+JT,N+NR+4*NS-2+JT,4
1978       P(IN1,1)=2-JT
1979       P(IN1,2)=JT-1
1980   580 P(IN1,3)=1.
1981
1982 C...Initialize flavour and pT variables for open string.
1983       IF(NS.LT.NR) THEN
1984         PX(1)=0.
1985         PY(1)=0.
1986         IF(NS.EQ.1.AND.MJU(1)+MJU(2).EQ.0) CALL LUPTDI(0,PX(1),PY(1))
1987         PX(2)=-PX(1)
1988         PY(2)=-PY(1)
1989         DO 590 JT=1,2
1990         KFL(JT)=K(IE(JT),2)
1991         IF(MJU(JT).NE.0) KFL(JT)=KFJS(JT)
1992         MSTJ(93)=1
1993         PMQ(JT)=ULMASS(KFL(JT))
1994   590   GAM(JT)=0.
1995
1996 C...Closed string: random initial breakup flavour, pT and vertex.
1997       ELSE
1998         KFL(3)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
1999         CALL LUKFDI(KFL(3),0,KFL(1),KDUMP)
2000         KFL(2)=-KFL(1)
2001         IF(IABS(KFL(1)).GT.10.AND.RLU(0).GT.0.5) THEN
2002           KFL(2)=-(KFL(1)+ISIGN(10000,KFL(1)))
2003         ELSEIF(IABS(KFL(1)).GT.10) THEN
2004           KFL(1)=-(KFL(2)+ISIGN(10000,KFL(2)))
2005         ENDIF
2006         CALL LUPTDI(KFL(1),PX(1),PY(1))
2007         PX(2)=-PX(1)
2008         PY(2)=-PY(1)
2009         PR3=MIN(25.,0.1*P(N+NR+1,5)**2)
2010   600   CALL LUZDIS(KFL(1),KFL(2),PR3,Z)
2011         ZR=PR3/(Z*P(N+NR+1,5)**2)
2012         IF(ZR.GE.1.) GOTO 600
2013         DO 610 JT=1,2
2014         MSTJ(93)=1
2015         PMQ(JT)=ULMASS(KFL(JT))
2016         GAM(JT)=PR3*(1.-Z)/Z
2017         IN1=N+NR+3+4*(JT/2)*(NS-1)
2018         P(IN1,JT)=1.-Z
2019         P(IN1,3-JT)=JT-1
2020         P(IN1,3)=(2-JT)*(1.-Z)+(JT-1)*Z
2021         P(IN1+1,JT)=ZR
2022         P(IN1+1,3-JT)=2-JT
2023   610   P(IN1+1,3)=(2-JT)*(1.-ZR)+(JT-1)*ZR
2024       ENDIF
2025
2026 C...Find initial transverse directions (i.e. spacelike four-vectors).
2027       DO 650 JT=1,2
2028       IF(JT.EQ.1.OR.NS.EQ.NR-1) THEN
2029         IN1=IN(3*JT+1)
2030         IN3=IN(3*JT+3)
2031         DO 620 J=1,4
2032         DP(1,J)=P(IN1,J)
2033         DP(2,J)=P(IN1+1,J)
2034         DP(3,J)=0.
2035   620   DP(4,J)=0.
2036         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
2037         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
2038         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
2039         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
2040         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
2041         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
2042         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
2043         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
2044         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
2045         DHC12=DFOUR(1,2)
2046         DHCX1=DFOUR(3,1)/DHC12
2047         DHCX2=DFOUR(3,2)/DHC12
2048         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
2049         DHCY1=DFOUR(4,1)/DHC12
2050         DHCY2=DFOUR(4,2)/DHC12
2051         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
2052         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
2053         DO 630 J=1,4
2054         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
2055         P(IN3,J)=DP(3,J)
2056   630   P(IN3+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
2057      &  DHCYX*DP(3,J))
2058       ELSE
2059         DO 640 J=1,4
2060         P(IN3+2,J)=P(IN3,J)
2061   640   P(IN3+3,J)=P(IN3+1,J)
2062       ENDIF
2063   650 CONTINUE
2064
2065 C...Remove energy used up in junction string fragmentation.
2066       IF(MJU(1)+MJU(2).GT.0) THEN
2067         DO 670 JT=1,2
2068         IF(NJS(JT).EQ.0) GOTO 670
2069         DO 660 J=1,4
2070   660   P(N+NRS,J)=P(N+NRS,J)-PJS(JT+2,J)
2071   670   CONTINUE
2072       ENDIF
2073
2074 C...Produce new particle: side, origin.
2075   680 I=I+1
2076       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN
2077         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS')
2078         IF(MSTU(21).GE.1) RETURN
2079       ENDIF
2080       JT=1.5+RLU(0)
2081       IF(IABS(KFL(3-JT)).GT.10) JT=3-JT
2082       IF(IABS(KFL(3-JT)).GE.4.AND.IABS(KFL(3-JT)).LE.8) JT=3-JT
2083       JR=3-JT
2084       JS=3-2*JT
2085       IRANK(JT)=IRANK(JT)+1
2086       K(I,1)=1
2087       K(I,3)=IE(JT)
2088       K(I,4)=0
2089       K(I,5)=0
2090
2091 C...Generate flavour, hadron and pT.
2092   690 CALL LUKFDI(KFL(JT),0,KFL(3),K(I,2))
2093       IF(K(I,2).EQ.0) GOTO 560
2094       IF(MSTJ(12).GE.3.AND.IRANK(JT).EQ.1.AND.IABS(KFL(JT)).LE.10.AND.
2095      &IABS(KFL(3)).GT.10) THEN
2096         IF(RLU(0).GT.PARJ(19)) GOTO 690
2097       ENDIF
2098       P(I,5)=ULMASS(K(I,2))
2099       CALL LUPTDI(KFL(JT),PX(3),PY(3))
2100       PR(JT)=P(I,5)**2+(PX(JT)+PX(3))**2+(PY(JT)+PY(3))**2
2101
2102 C...Final hadrons for small invariant mass.
2103       MSTJ(93)=1
2104       PMQ(3)=ULMASS(KFL(3))
2105       PARJST=PARJ(33)
2106       IF(MSTJ(11).EQ.2) PARJST=PARJ(34)
2107       WMIN=PARJST+PMQ(1)+PMQ(2)+PARJ(36)*PMQ(3)
2108       IF(IABS(KFL(JT)).GT.10.AND.IABS(KFL(3)).GT.10) WMIN=
2109      &WMIN-0.5*PARJ(36)*PMQ(3)
2110       WREM2=FOUR(N+NRS,N+NRS)
2111       IF(WREM2.LT.0.10) GOTO 560
2112       IF(WREM2.LT.MAX(WMIN*(1.+(2.*RLU(0)-1.)*PARJ(37)),
2113      &PARJ(32)+PMQ(1)+PMQ(2))**2) GOTO 820
2114
2115 C...Choose z, which gives Gamma. Shift z for heavy flavours.
2116       CALL LUZDIS(KFL(JT),KFL(3),PR(JT),Z)
2117       IF(IABS(KFL(JT)).GE.4.AND.IABS(KFL(JT)).LE.8.AND.
2118      &MSTU(90).LT.8) THEN
2119         MSTU(90)=MSTU(90)+1
2120         MSTU(90+MSTU(90))=I
2121         PARU(90+MSTU(90))=Z
2122       ENDIF
2123       KFL1A=IABS(KFL(1))
2124       KFL2A=IABS(KFL(2))
2125       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
2126      &MOD(KFL2A/1000,10)).GE.4) THEN
2127         PR(JR)=(PMQ(JR)+PMQ(3))**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
2128         PW12=SQRT(MAX(0.,(WREM2-PR(1)-PR(2))**2-4.*PR(1)*PR(2)))
2129         Z=(WREM2+PR(JT)-PR(JR)+PW12*(2.*Z-1.))/(2.*WREM2)
2130         PR(JR)=(PMQ(JR)+PARJST)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
2131         IF((1.-Z)*(WREM2-PR(JT)/Z).LT.PR(JR)) GOTO 820
2132       ENDIF
2133       GAM(3)=(1.-Z)*(GAM(JT)+PR(JT)/Z)
2134       DO 700 J=1,3
2135   700 IN(J)=IN(3*JT+J)
2136
2137 C...Stepping within or from 'low' string region easy.
2138       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)*
2139      &P(IN(1),5)**2.GE.PR(JT)) THEN
2140         P(IN(JT)+2,4)=Z*P(IN(JT)+2,3)
2141         P(IN(JR)+2,4)=PR(JT)/(P(IN(JT)+2,4)*P(IN(1),5)**2)
2142         DO 710 J=1,4
2143   710   P(I,J)=(PX(JT)+PX(3))*P(IN(3),J)+(PY(JT)+PY(3))*P(IN(3)+1,J)
2144         GOTO 780
2145       ELSEIF(IN(1)+1.EQ.IN(2)) THEN
2146         P(IN(JR)+2,4)=P(IN(JR)+2,3)
2147         P(IN(JR)+2,JT)=1.
2148         IN(JR)=IN(JR)+4*JS
2149         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 560
2150         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
2151           P(IN(JT)+2,4)=P(IN(JT)+2,3)
2152           P(IN(JT)+2,JT)=0.
2153           IN(JT)=IN(JT)+4*JS
2154         ENDIF
2155       ENDIF
2156
2157 C...Find new transverse directions (i.e. spacelike string vectors).
2158   720 IF(JS*IN(1).GT.JS*IN(3*JR+1).OR.JS*IN(2).GT.JS*IN(3*JR+2).OR.
2159      &IN(1).GT.IN(2)) GOTO 560
2160       IF(IN(1).NE.IN(3*JT+1).OR.IN(2).NE.IN(3*JT+2)) THEN
2161         DO 730 J=1,4
2162         DP(1,J)=P(IN(1),J)
2163         DP(2,J)=P(IN(2),J)
2164         DP(3,J)=0.
2165   730   DP(4,J)=0.
2166         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2)
2167         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2)
2168         DHC12=DFOUR(1,2)
2169         IF(DHC12.LE.1E-2) THEN
2170           P(IN(JT)+2,4)=P(IN(JT)+2,3)
2171           P(IN(JT)+2,JT)=0.
2172           IN(JT)=IN(JT)+4*JS
2173           GOTO 720
2174         ENDIF
2175         IN(3)=N+NR+4*NS+5
2176         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4)
2177         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4)
2178         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4)
2179         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1.
2180         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1.
2181         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1.
2182         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1.
2183         DHCX1=DFOUR(3,1)/DHC12
2184         DHCX2=DFOUR(3,2)/DHC12
2185         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12)
2186         DHCY1=DFOUR(4,1)/DHC12
2187         DHCY2=DFOUR(4,2)/DHC12
2188         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12
2189         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2)
2190         DO 740 J=1,4
2191         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J))
2192         P(IN(3),J)=DP(3,J)
2193   740   P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)-
2194      &  DHCYX*DP(3,J))
2195 C...Express pT with respect to new axes, if sensible.
2196         PXP=-(PX(3)*FOUR(IN(3*JT+3),IN(3))+PY(3)*
2197      &  FOUR(IN(3*JT+3)+1,IN(3)))
2198         PYP=-(PX(3)*FOUR(IN(3*JT+3),IN(3)+1)+PY(3)*
2199      &  FOUR(IN(3*JT+3)+1,IN(3)+1))
2200         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN
2201           PX(3)=PXP
2202           PY(3)=PYP
2203         ENDIF
2204       ENDIF
2205
2206 C...Sum up known four-momentum. Gives coefficients for m2 expression.
2207       DO 760 J=1,4
2208       DHG(J)=0.
2209       P(I,J)=PX(JT)*P(IN(3*JT+3),J)+PY(JT)*P(IN(3*JT+3)+1,J)+
2210      &PX(3)*P(IN(3),J)+PY(3)*P(IN(3)+1,J)
2211       DO 750 IN1=IN(3*JT+1),IN(1)-4*JS,4*JS
2212   750 P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J)
2213       DO 760 IN2=IN(3*JT+2),IN(2)-4*JS,4*JS
2214   760 P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J)
2215       DHM(1)=FOUR(I,I)
2216       DHM(2)=2.*FOUR(I,IN(1))
2217       DHM(3)=2.*FOUR(I,IN(2))
2218       DHM(4)=2.*FOUR(IN(1),IN(2))
2219
2220 C...Find coefficients for Gamma expression.
2221       DO 770 IN2=IN(1)+1,IN(2),4
2222       DO 770 IN1=IN(1),IN2-1,4
2223       DHC=2.*FOUR(IN1,IN2)
2224       DHG(1)=DHG(1)+P(IN1+2,JT)*P(IN2+2,JT)*DHC
2225       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-JS*P(IN2+2,JT)*DHC
2226       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+JS*P(IN1+2,JT)*DHC
2227   770 IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC
2228
2229 C...Solve (m2, Gamma) equation system for energies taken.
2230       DHS1=DHM(JR+1)*DHG(4)-DHM(4)*DHG(JR+1)
2231       IF(ABS(DHS1).LT.1E-4) GOTO 560
2232       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(JT+1)*DHG(JR+1)-DHG(4)*
2233      &(P(I,5)**2-DHM(1))+DHG(JT+1)*DHM(JR+1)
2234       DHS3=DHM(JT+1)*(GAM(3)-DHG(1))-DHG(JT+1)*(P(I,5)**2-DHM(1))
2235       P(IN(JR)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)-
2236      &DHS2/DHS1)
2237       IF(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4).LE.0.) GOTO 560
2238       P(IN(JT)+2,4)=(P(I,5)**2-DHM(1)-DHM(JR+1)*P(IN(JR)+2,4))/
2239      &(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4))
2240
2241 C...Step to new region if necessary.
2242       IF(P(IN(JR)+2,4).GT.P(IN(JR)+2,3)) THEN
2243         P(IN(JR)+2,4)=P(IN(JR)+2,3)
2244         P(IN(JR)+2,JT)=1.
2245         IN(JR)=IN(JR)+4*JS
2246         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 560
2247         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN
2248           P(IN(JT)+2,4)=P(IN(JT)+2,3)
2249           P(IN(JT)+2,JT)=0.
2250           IN(JT)=IN(JT)+4*JS
2251         ENDIF
2252         GOTO 720
2253       ELSEIF(P(IN(JT)+2,4).GT.P(IN(JT)+2,3)) THEN
2254         P(IN(JT)+2,4)=P(IN(JT)+2,3)
2255         P(IN(JT)+2,JT)=0.
2256         IN(JT)=IN(JT)+4*JS
2257         GOTO 720
2258       ENDIF
2259
2260 C...Four-momentum of particle. Remaining quantities. Loop back.
2261   780 DO 790 J=1,4
2262       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J)
2263   790 P(N+NRS,J)=P(N+NRS,J)-P(I,J)
2264       IF(P(I,4).LT.P(I,5)) GOTO 560
2265       KFL(JT)=-KFL(3)
2266       PMQ(JT)=PMQ(3)
2267       PX(JT)=-PX(3)
2268       PY(JT)=-PY(3)
2269       GAM(JT)=GAM(3)    
2270       IF(IN(3).NE.IN(3*JT+3)) THEN
2271         DO 800 J=1,4
2272         P(IN(3*JT+3),J)=P(IN(3),J)
2273   800   P(IN(3*JT+3)+1,J)=P(IN(3)+1,J)
2274       ENDIF
2275       DO 810 JQ=1,2
2276       IN(3*JT+JQ)=IN(JQ)
2277       P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4)
2278   810 P(IN(JQ)+2,JT)=P(IN(JQ)+2,JT)-JS*(3-2*JQ)*P(IN(JQ)+2,4)
2279       GOTO 680
2280
2281 C...Final hadron: side, flavour, hadron, mass.
2282   820 I=I+1
2283       K(I,1)=1
2284       K(I,3)=IE(JR)
2285       K(I,4)=0
2286       K(I,5)=0
2287       CALL LUKFDI(KFL(JR),-KFL(3),KFLDMP,K(I,2))
2288       IF(K(I,2).EQ.0) GOTO 560
2289       P(I,5)=ULMASS(K(I,2))
2290       PR(JR)=P(I,5)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2
2291
2292 C...Final two hadrons: find common setup of four-vectors.
2293       JQ=1
2294       IF(P(IN(4)+2,3)*P(IN(5)+2,3)*FOUR(IN(4),IN(5)).LT.P(IN(7),3)*
2295      &P(IN(8),3)*FOUR(IN(7),IN(8))) JQ=2
2296       DHC12=FOUR(IN(3*JQ+1),IN(3*JQ+2))
2297       DHR1=FOUR(N+NRS,IN(3*JQ+2))/DHC12
2298       DHR2=FOUR(N+NRS,IN(3*JQ+1))/DHC12
2299       IF(IN(4).NE.IN(7).OR.IN(5).NE.IN(8)) THEN
2300         PX(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3))-PX(JQ)
2301         PY(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3)+1)-PY(JQ)
2302         PR(3-JQ)=P(I+(JT+JQ-3)**2-1,5)**2+(PX(3-JQ)+(2*JQ-3)*JS*
2303      &  PX(3))**2+(PY(3-JQ)+(2*JQ-3)*JS*PY(3))**2
2304       ENDIF
2305
2306 C...Solve kinematics for final two hadrons, if possible.
2307       WREM2=WREM2+(PX(1)+PX(2))**2+(PY(1)+PY(2))**2
2308       FD=(SQRT(PR(1))+SQRT(PR(2)))/SQRT(WREM2)
2309       IF(MJU(1)+MJU(2).NE.0.AND.I.EQ.ISAV+2.AND.FD.GE.1.) GOTO 190
2310       IF(FD.GE.1.) GOTO 560
2311       FA=WREM2+PR(JT)-PR(JR)
2312       IF(MSTJ(11).NE.2) PREV=0.5*EXP(MAX(-80.,LOG(FD)*PARJ(38)*
2313      &(PR(1)+PR(2))**2))
2314       IF(MSTJ(11).EQ.2) PREV=0.5*FD**PARJ(39)
2315       FB=SIGN(SQRT(MAX(0.,FA**2-4.*WREM2*PR(JT))),JS*(RLU(0)-PREV))
2316       KFL1A=IABS(KFL(1))
2317       KFL2A=IABS(KFL(2))
2318       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10),
2319      &MOD(KFL2A/1000,10)).GE.6) FB=SIGN(SQRT(MAX(0.,FA**2-
2320      &4.*WREM2*PR(JT))),FLOAT(JS))
2321       DO 830 J=1,4
2322       P(I-1,J)=(PX(JT)+PX(3))*P(IN(3*JQ+3),J)+(PY(JT)+PY(3))*
2323      &P(IN(3*JQ+3)+1,J)+0.5*(DHR1*(FA+FB)*P(IN(3*JQ+1),J)+
2324      &DHR2*(FA-FB)*P(IN(3*JQ+2),J))/WREM2
2325   830 P(I,J)=P(N+NRS,J)-P(I-1,J)
2326       IF(P(I-1,4).LT.P(I-1,5).OR.P(I,4).LT.P(I,5)) GOTO 560
2327 C...Mark jets as fragmented and give daughter pointers.
2328       N=I-NRS+1
2329       DO 840 I=NSAV+1,NSAV+NP
2330       IM=K(I,3)
2331       K(IM,1)=K(IM,1)+10
2332       IF(MSTU(16).NE.2) THEN
2333         K(IM,4)=NSAV+1
2334         K(IM,5)=NSAV+1
2335       ELSE
2336         K(IM,4)=NSAV+2
2337         K(IM,5)=N
2338       ENDIF
2339   840 CONTINUE
2340
2341 C...Document string system. Move up particles.
2342       NSAV=NSAV+1
2343       K(NSAV,1)=11
2344       K(NSAV,2)=92
2345       K(NSAV,3)=IP
2346       K(NSAV,4)=NSAV+1
2347       K(NSAV,5)=N
2348       DO 850 J=1,4
2349       P(NSAV,J)=DPS(J)
2350   850 V(NSAV,J)=V(IP,J)
2351       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
2352       V(NSAV,5)=0.
2353       DO 860 I=NSAV+1,N
2354       DO 860 J=1,5
2355       K(I,J)=K(I+NRS-1,J)
2356       P(I,J)=P(I+NRS-1,J)
2357   860 V(I,J)=0.
2358       MSTU91=MSTU(90)
2359       DO 870 IZ=MSTU90+1,MSTU91
2360       MSTU9T(IZ)=MSTU(90+IZ)-NRS+1-NSAV+N
2361   870 PARU9T(IZ)=PARU(90+IZ)
2362       MSTU(90)=MSTU90
2363
2364 C...Order particles in rank along the chain. Update mother pointer.
2365       DO 880 I=NSAV+1,N
2366       DO 880 J=1,5
2367       K(I-NSAV+N,J)=K(I,J)
2368   880 P(I-NSAV+N,J)=P(I,J)
2369       I1=NSAV
2370       DO 910 I=N+1,2*N-NSAV
2371       IF(K(I,3).NE.IE(1)) GOTO 910
2372       I1=I1+1
2373       DO 890 J=1,5
2374       K(I1,J)=K(I,J)
2375   890 P(I1,J)=P(I,J)
2376       IF(MSTU(16).NE.2) K(I1,3)=NSAV
2377       DO 900 IZ=MSTU90+1,MSTU91
2378       IF(MSTU9T(IZ).EQ.I) THEN
2379         MSTU(90)=MSTU(90)+1
2380         MSTU(90+MSTU(90))=I1
2381         PARU(90+MSTU(90))=PARU9T(IZ)
2382       ENDIF
2383   900 CONTINUE
2384   910 CONTINUE
2385       DO 940 I=2*N-NSAV,N+1,-1
2386       IF(K(I,3).EQ.IE(1)) GOTO 940
2387       I1=I1+1
2388       DO 920 J=1,5
2389       K(I1,J)=K(I,J)
2390   920 P(I1,J)=P(I,J)
2391       IF(MSTU(16).NE.2) K(I1,3)=NSAV
2392       DO 930 IZ=MSTU90+1,MSTU91
2393       IF(MSTU9T(IZ).EQ.I) THEN
2394         MSTU(90)=MSTU(90)+1
2395         MSTU(90+MSTU(90))=I1
2396         PARU(90+MSTU(90))=PARU9T(IZ)
2397       ENDIF
2398   930 CONTINUE
2399   940 CONTINUE
2400
2401 C...Boost back particle system. Set production vertices.
2402       IF(MBST.EQ.0) THEN
2403         MSTU(33)=1
2404         CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4),
2405      &  DPS(3)/DPS(4))
2406       ELSE
2407         DO 950 I=NSAV+1,N
2408         HHPMT=P(I,1)**2+P(I,2)**2+P(I,5)**2
2409         IF(P(I,3).GT.0.) THEN
2410           HHPEZ=(P(I,4)+P(I,3))*HHBZ
2411           P(I,3)=0.5*(HHPEZ-HHPMT/HHPEZ)
2412           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ)
2413         ELSE
2414           HHPEZ=(P(I,4)-P(I,3))/HHBZ
2415           P(I,3)=-0.5*(HHPEZ-HHPMT/HHPEZ)
2416           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ)
2417         ENDIF
2418   950   CONTINUE
2419       ENDIF
2420       DO 960 I=NSAV+1,N 
2421       DO 960 J=1,4
2422   960 V(I,J)=V(IP,J)
2423
2424       RETURN
2425       END
2426
2427 C*********************************************************************
2428
2429       SUBROUTINE LUINDF(IP)
2430
2431 C...Purpose: to handle the fragmentation of a jet system (or a single
2432 C...jet) according to independent fragmentation models.
2433       IMPLICIT DOUBLE PRECISION(D)
2434       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
2435       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2436       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
2437       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
2438       DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
2439      &KFLO(2),PXO(2),PYO(2),WO(2)
2440
2441 C...Reset counters. Identify parton system and take copy. Check flavour.
2442       NSAV=N
2443       MSTU90=MSTU(90)
2444       NJET=0
2445       KQSUM=0
2446       DO 100 J=1,5
2447   100 DPS(J)=0.
2448       I=IP-1
2449   110 I=I+1
2450       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
2451         CALL LUERRM(12,'(LUINDF:) failed to reconstruct jet system')
2452         IF(MSTU(21).GE.1) RETURN
2453       ENDIF
2454       IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
2455       KC=LUCOMP(K(I,2))
2456       IF(KC.EQ.0) GOTO 110
2457       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
2458       IF(KQ.EQ.0) GOTO 110
2459       NJET=NJET+1
2460       IF(KQ.NE.2) KQSUM=KQSUM+KQ
2461       DO 120 J=1,5
2462       K(NSAV+NJET,J)=K(I,J)
2463       P(NSAV+NJET,J)=P(I,J)
2464   120 DPS(J)=DPS(J)+P(I,J)
2465       K(NSAV+NJET,3)=I
2466       IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
2467      &K(I+1,1).EQ.2)) GOTO 110
2468       IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
2469         CALL LUERRM(12,'(LUINDF:) unphysical flavour combination')
2470         IF(MSTU(21).GE.1) RETURN
2471       ENDIF
2472
2473 C...Boost copied system to CM frame. Find CM energy and sum flavours.
2474       IF(NJET.NE.1) THEN
2475         MSTU(33)=1
2476         CALL LUDBRB(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)/DPS(4),
2477      &  -DPS(2)/DPS(4),-DPS(3)/DPS(4))
2478       ENDIF
2479       PECM=0.
2480       DO 130 J=1,3
2481   130 NFI(J)=0
2482       DO 140 I=NSAV+1,NSAV+NJET
2483       PECM=PECM+P(I,4)
2484       KFA=IABS(K(I,2))
2485       IF(KFA.LE.3) THEN
2486         NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
2487       ELSEIF(KFA.GT.1000) THEN
2488         KFLA=MOD(KFA/1000,10)
2489         KFLB=MOD(KFA/100,10)
2490         IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
2491         IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
2492       ENDIF
2493   140 CONTINUE
2494
2495 C...Loop over attempts made. Reset counters.
2496       NTRY=0
2497   150 NTRY=NTRY+1
2498       IF(NTRY.GT.200) THEN
2499         CALL LUERRM(14,'(LUINDF:) caught in infinite loop')
2500         IF(MSTU(21).GE.1) RETURN
2501       ENDIF
2502       N=NSAV+NJET
2503       MSTU(90)=MSTU90
2504       DO 160 J=1,3
2505       NFL(J)=NFI(J)
2506       IFET(J)=0
2507   160 KFLF(J)=0
2508
2509 C...Loop over jets to be fragmented.
2510       DO 230 IP1=NSAV+1,NSAV+NJET
2511       MSTJ(91)=0
2512       NSAV1=N
2513       MSTU91=MSTU(90)
2514
2515 C...Initial flavour and momentum values. Jet along +z axis.
2516       KFLH=IABS(K(IP1,2))
2517       IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
2518       KFLO(2)=0
2519       WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
2520
2521 C...Initial values for quark or diquark jet.
2522   170 IF(IABS(K(IP1,2)).NE.21) THEN
2523         NSTR=1
2524         KFLO(1)=K(IP1,2)
2525         CALL LUPTDI(0,PXO(1),PYO(1))
2526         WO(1)=WF
2527
2528 C...Initial values for gluon treated like random quark jet.
2529       ELSEIF(MSTJ(2).LE.2) THEN
2530         NSTR=1
2531         IF(MSTJ(2).EQ.2) MSTJ(91)=1
2532         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
2533         CALL LUPTDI(0,PXO(1),PYO(1))
2534         WO(1)=WF
2535
2536 C...Initial values for gluon treated like quark-antiquark jet pair,
2537 C...sharing energy according to Altarelli-Parisi splitting function.
2538       ELSE
2539         NSTR=2
2540         IF(MSTJ(2).EQ.4) MSTJ(91)=1
2541         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
2542         KFLO(2)=-KFLO(1)
2543         CALL LUPTDI(0,PXO(1),PYO(1))
2544         PXO(2)=-PXO(1)
2545         PYO(2)=-PYO(1)
2546         WO(1)=WF*RLU(0)**(1./3.)
2547         WO(2)=WF-WO(1)
2548       ENDIF
2549
2550 C...Initial values for rank, flavour, pT and W+.
2551       DO 220 ISTR=1,NSTR
2552   180 I=N
2553       MSTU(90)=MSTU91
2554       IRANK=0
2555       KFL1=KFLO(ISTR)
2556       PX1=PXO(ISTR)
2557       PY1=PYO(ISTR)
2558       W=WO(ISTR)
2559
2560 C...New hadron. Generate flavour and hadron species.
2561   190 I=I+1
2562       IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
2563         CALL LUERRM(11,'(LUINDF:) no more memory left in LUJETS')
2564         IF(MSTU(21).GE.1) RETURN
2565       ENDIF
2566       IRANK=IRANK+1
2567       K(I,1)=1
2568       K(I,3)=IP1
2569       K(I,4)=0
2570       K(I,5)=0
2571   200 CALL LUKFDI(KFL1,0,KFL2,K(I,2))
2572       IF(K(I,2).EQ.0) GOTO 180
2573       IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.
2574      &IABS(KFL2).GT.10) THEN
2575         IF(RLU(0).GT.PARJ(19)) GOTO 200
2576       ENDIF
2577
2578 C...Find hadron mass. Generate four-momentum.
2579       P(I,5)=ULMASS(K(I,2))
2580       CALL LUPTDI(KFL1,PX2,PY2)
2581       P(I,1)=PX1+PX2
2582       P(I,2)=PY1+PY2
2583       PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
2584       CALL LUZDIS(KFL1,KFL2,PR,Z)
2585       MZSAV=0
2586       IF(IABS(KFL1).GE.4.AND.IABS(KFL1).LE.8.AND.MSTU(90).LT.8) THEN
2587         MZSAV=1
2588         MSTU(90)=MSTU(90)+1
2589         MSTU(90+MSTU(90))=I
2590         PARU(90+MSTU(90))=Z
2591       ENDIF
2592       P(I,3)=0.5*(Z*W-PR/(Z*W))
2593       P(I,4)=0.5*(Z*W+PR/(Z*W))
2594       IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
2595      &P(I,3).LE.0.001) THEN
2596         IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180
2597         P(I,3)=0.0001
2598         P(I,4)=SQRT(PR)
2599         Z=P(I,4)/W
2600       ENDIF
2601
2602 C...Remaining flavour and momentum.
2603       KFL1=-KFL2
2604       PX1=-PX2
2605       PY1=-PY2
2606       W=(1.-Z)*W
2607       DO 210 J=1,5
2608   210 V(I,J)=0.
2609
2610 C...Check if pL acceptable. Go back for new hadron if enough energy.
2611       IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) THEN
2612         I=I-1
2613         IF(MZSAV.EQ.1) MSTU(90)=MSTU(90)-1
2614       ENDIF
2615       IF(W.GT.PARJ(31)) GOTO 190
2616   220 N=I
2617       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32)
2618       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
2619
2620 C...Rotate jet to new direction.
2621       THE=ULANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
2622       PHI=ULANGL(P(IP1,1),P(IP1,2))
2623       MSTU(33)=1
2624       CALL LUDBRB(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
2625       K(K(IP1,3),4)=NSAV1+1
2626       K(K(IP1,3),5)=N
2627
2628 C...End of jet generation loop. Skip conservation in some cases.
2629   230 CONTINUE
2630       IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 470
2631       IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
2632
2633 C...Subtract off produced hadron flavours, finished if zero.
2634       DO 240 I=NSAV+NJET+1,N
2635       KFA=IABS(K(I,2))
2636       KFLA=MOD(KFA/1000,10)
2637       KFLB=MOD(KFA/100,10)
2638       KFLC=MOD(KFA/10,10)
2639       IF(KFLA.EQ.0) THEN
2640         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
2641         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
2642       ELSE
2643         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
2644         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
2645         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
2646       ENDIF
2647   240 CONTINUE
2648       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
2649      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
2650       IF(NREQ.EQ.0) GOTO 320
2651
2652 C...Take away flavour of low-momentum particles until enough freedom.
2653       NREM=0
2654   250 IREM=0
2655       P2MIN=PECM**2
2656       DO 260 I=NSAV+NJET+1,N
2657       P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
2658       IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
2659   260 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
2660       IF(IREM.EQ.0) GOTO 150
2661       K(IREM,1)=7
2662       KFA=IABS(K(IREM,2))
2663       KFLA=MOD(KFA/1000,10)
2664       KFLB=MOD(KFA/100,10)
2665       KFLC=MOD(KFA/10,10)
2666       IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
2667       IF(K(IREM,1).EQ.8) GOTO 250
2668       IF(KFLA.EQ.0) THEN
2669         ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
2670         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
2671         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
2672       ELSE
2673         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
2674         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
2675         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
2676       ENDIF
2677       NREM=NREM+1
2678       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
2679      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
2680       IF(NREQ.GT.NREM) GOTO 250
2681       DO 270 I=NSAV+NJET+1,N
2682   270 IF(K(I,1).EQ.8) K(I,1)=1
2683
2684 C...Find combination of existing and new flavours for hadron.
2685   280 NFET=2
2686       IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
2687       IF(NREQ.LT.NREM) NFET=1
2688       IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
2689       DO 290 J=1,NFET
2690       IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU(0)
2691       KFLF(J)=ISIGN(1,NFL(1))
2692       IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
2693   290 IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
2694       IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
2695      &GOTO 280
2696       IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
2697      &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3).
2698      &LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
2699       IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU(0))
2700       IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
2701       IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU(0)),-KFLF(1))
2702       IF(NFET.LE.2) KFLF(3)=0
2703       IF(KFLF(3).NE.0) THEN
2704         KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
2705      &  100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
2706         IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU(0).GT.1.)
2707      &  KFLFC=KFLFC+ISIGN(2,KFLFC)
2708       ELSE
2709         KFLFC=KFLF(1)
2710       ENDIF
2711       CALL LUKFDI(KFLFC,KFLF(2),KFLDMP,KF)
2712       IF(KF.EQ.0) GOTO 280
2713       DO 300 J=1,MAX(2,NFET)
2714   300 NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
2715
2716 C...Store hadron at random among free positions.
2717       NPOS=MIN(1+INT(RLU(0)*NREM),NREM)
2718       DO 310 I=NSAV+NJET+1,N
2719       IF(K(I,1).EQ.7) NPOS=NPOS-1
2720       IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
2721       K(I,1)=1
2722       K(I,2)=KF
2723       P(I,5)=ULMASS(K(I,2))
2724       P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
2725   310 CONTINUE
2726       NREM=NREM-1
2727       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
2728      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
2729       IF(NREM.GT.0) GOTO 280
2730
2731 C...Compensate for missing momentum in global scheme (3 options).
2732   320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
2733         DO 330 J=1,3
2734         PSI(J)=0.
2735         DO 330 I=NSAV+NJET+1,N
2736   330   PSI(J)=PSI(J)+P(I,J)
2737         PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
2738         PWS=0.
2739         DO 340 I=NSAV+NJET+1,N
2740         IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
2741         IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
2742      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
2743   340   IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1.
2744         DO 360 I=NSAV+NJET+1,N
2745         IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
2746         IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
2747      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
2748         IF(MOD(MSTJ(3),5).EQ.3) PW=1.
2749         DO 350 J=1,3
2750   350   P(I,J)=P(I,J)-PSI(J)*PW/PWS
2751   360   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
2752
2753 C...Compensate for missing momentum withing each jet separately.
2754       ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
2755         DO 370 I=N+1,N+NJET
2756         K(I,1)=0
2757         DO 370 J=1,5
2758   370   P(I,J)=0.
2759         DO 390 I=NSAV+NJET+1,N
2760         IR1=K(I,3)
2761         IR2=N+IR1-NSAV
2762         K(IR2,1)=K(IR2,1)+1
2763         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
2764      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
2765         DO 380 J=1,3
2766   380   P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
2767         P(IR2,4)=P(IR2,4)+P(I,4)
2768   390   P(IR2,5)=P(IR2,5)+PLS
2769         PSS=0.
2770         DO 400 I=N+1,N+NJET
2771   400   IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2))
2772         DO 420 I=NSAV+NJET+1,N
2773         IR1=K(I,3)
2774         IR2=N+IR1-NSAV
2775         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
2776      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
2777         DO 410 J=1,3
2778   410   P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS*
2779      &  P(IR1,J)
2780   420   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
2781       ENDIF
2782
2783 C...Scale momenta for energy conservation.
2784       IF(MOD(MSTJ(3),5).NE.0) THEN
2785         PMS=0.
2786         PES=0.
2787         PQS=0.
2788         DO 430 I=NSAV+NJET+1,N
2789         PMS=PMS+P(I,5)
2790         PES=PES+P(I,4)
2791   430   PQS=PQS+P(I,5)**2/P(I,4)
2792         IF(PMS.GE.PECM) GOTO 150
2793         NECO=0
2794   440   NECO=NECO+1
2795         PFAC=(PECM-PQS)/(PES-PQS)
2796         PES=0.
2797         PQS=0.
2798         DO 460 I=NSAV+NJET+1,N
2799         DO 450 J=1,3
2800   450   P(I,J)=PFAC*P(I,J)
2801         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
2802         PES=PES+P(I,4)
2803   460   PQS=PQS+P(I,5)**2/P(I,4)
2804         IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 440
2805       ENDIF
2806
2807 C...Origin of produced particles and parton daughter pointers.
2808   470 DO 480 I=NSAV+NJET+1,N
2809       IF(MSTU(16).NE.2) K(I,3)=NSAV+1
2810   480 IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
2811       DO 490 I=NSAV+1,NSAV+NJET
2812       I1=K(I,3)
2813       K(I1,1)=K(I1,1)+10
2814       IF(MSTU(16).NE.2) THEN
2815         K(I1,4)=NSAV+1
2816         K(I1,5)=NSAV+1
2817       ELSE
2818         K(I1,4)=K(I1,4)-NJET+1
2819         K(I1,5)=K(I1,5)-NJET+1
2820         IF(K(I1,5).LT.K(I1,4)) THEN
2821           K(I1,4)=0
2822           K(I1,5)=0
2823         ENDIF
2824       ENDIF
2825   490 CONTINUE
2826
2827 C...Document independent fragmentation system. Remove copy of jets.
2828       NSAV=NSAV+1
2829       K(NSAV,1)=11
2830       K(NSAV,2)=93
2831       K(NSAV,3)=IP
2832       K(NSAV,4)=NSAV+1
2833       K(NSAV,5)=N-NJET+1
2834       DO 500 J=1,4
2835       P(NSAV,J)=DPS(J)
2836   500 V(NSAV,J)=V(IP,J)
2837       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
2838       V(NSAV,5)=0.
2839       DO 510 I=NSAV+NJET,N
2840       DO 510 J=1,5
2841       K(I-NJET+1,J)=K(I,J)
2842       P(I-NJET+1,J)=P(I,J)
2843   510 V(I-NJET+1,J)=V(I,J)
2844       N=N-NJET+1
2845       DO 520 IZ=MSTU90+1,MSTU(90)
2846   520 MSTU(90+IZ)=MSTU(90+IZ)-NJET+1
2847
2848 C...Boost back particle system. Set production vertices.
2849       IF(NJET.NE.1) CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),
2850      &DPS(2)/DPS(4),DPS(3)/DPS(4))
2851       DO 530 I=NSAV+1,N
2852       DO 530 J=1,4
2853   530 V(I,J)=V(IP,J)
2854
2855       RETURN
2856       END
2857
2858 C*********************************************************************
2859
2860       SUBROUTINE LUDECY(IP)
2861
2862 C...Purpose: to handle the decay of unstable particles.
2863       COMMON/LUJETS/N,K(150000,5),P(150000,5),V(150000,5)
2864       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
2865       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
2866       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
2867       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
2868       DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3),
2869      &WTCOR(10)
2870       DATA WTCOR/2.,5.,15.,60.,250.,1500.,1.2E4,1.2E5,150.,16./
2871
2872 C...Functions: momentum in two-particle decays, four-product and
2873 C...matrix element times phase space in weak decays.
2874       PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A)
2875       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)
2876       HMEPS(HA)=((1.-HRQ-HA)**2+3.*HA*(1.+HRQ-HA))*
2877      &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA)
2878
2879 C...Initial values.
2880       NTRY=0
2881       NSAV=N
2882       KFA=IABS(K(IP,2))
2883       KFS=ISIGN(1,K(IP,2))
2884       KC=LUCOMP(KFA)
2885       MSTJ(92)=0
2886
2887 C...Choose lifetime and determine decay vertex.
2888       IF(K(IP,1).EQ.5) THEN
2889         V(IP,5)=0.
2890       ELSEIF(K(IP,1).NE.4) THEN
2891         V(IP,5)=-PMAS(KC,4)*LOG(RLU(0))
2892       ENDIF
2893       DO 100 J=1,4
2894   100 VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
2895
2896 C...Determine whether decay allowed or not.
2897       MOUT=0
2898       IF(MSTJ(22).EQ.2) THEN
2899         IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1
2900       ELSEIF(MSTJ(22).EQ.3) THEN
2901         IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1
2902       ELSEIF(MSTJ(22).EQ.4) THEN
2903         IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1
2904         IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1
2905       ENDIF
2906       IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN
2907         K(IP,1)=4
2908         RETURN
2909       ENDIF
2910
2911 C...B-B~ mixing: flip sign of meson appropriately.
2912       MMIX=0
2913       IF((KFA.EQ.511.OR.KFA.EQ.531).AND.MSTJ(26).GE.1) THEN
2914         XBBMIX=PARJ(76)
2915         IF(KFA.EQ.531) XBBMIX=PARJ(77)
2916         IF(SIN(0.5*XBBMIX*V(IP,5)/PMAS(KC,4))**2.GT.RLU(0)) MMIX=1
2917         IF(MMIX.EQ.1) KFS=-KFS
2918       ENDIF
2919
2920 C...Check existence of decay channels. Particle/antiparticle rules.
2921       KCA=KC
2922       IF(MDCY(KC,2).GT.0) THEN
2923         MDMDCY=MDME(MDCY(KC,2),2)
2924         IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY
2925       ENDIF
2926       IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN
2927         CALL LUERRM(9,'(LUDECY:) no decay channel defined')
2928         RETURN
2929       ENDIF
2930       IF(MOD(KFA/1000,10).EQ.0.AND.(KCA.EQ.85.OR.KCA.EQ.87)) KFS=-KFS
2931       IF(KCHG(KC,3).EQ.0) THEN
2932         KFSP=1
2933         KFSN=0
2934         IF(RLU(0).GT.0.5) KFS=-KFS
2935       ELSEIF(KFS.GT.0) THEN
2936         KFSP=1
2937         KFSN=0
2938       ELSE
2939         KFSP=0
2940         KFSN=1
2941       ENDIF
2942
2943 C...Sum branching ratios of allowed decay channels.
2944   110 NOPE=0
2945       BRSU=0.
2946       DO 120 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1
2947       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
2948      &KFSN*MDME(IDL,1).NE.3) GOTO 120
2949       IF(MDME(IDL,2).GT.100) GOTO 120
2950       NOPE=NOPE+1
2951       BRSU=BRSU+BRAT(IDL)
2952   120 CONTINUE
2953       IF(NOPE.EQ.0) THEN
2954         CALL LUERRM(2,'(LUDECY:) all decay channels closed by user')
2955         RETURN
2956       ENDIF
2957
2958 C...Select decay channel among allowed ones.
2959   130 RBR=BRSU*RLU(0)
2960       IDL=MDCY(KCA,2)-1
2961   140 IDL=IDL+1
2962       IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
2963      &KFSN*MDME(IDL,1).NE.3) THEN
2964         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140
2965       ELSEIF(MDME(IDL,2).GT.100) THEN
2966         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140
2967       ELSE
2968         IDC=IDL
2969         RBR=RBR-BRAT(IDL)
2970         IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0.) GOTO 140
2971       ENDIF
2972
2973 C...Start readout of decay channel: matrix element, reset counters.
2974       MMAT=MDME(IDC,2)
2975   150 NTRY=NTRY+1
2976       IF(NTRY.GT.1000) THEN
2977         CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
2978         IF(MSTU(21).GE.1) RETURN
2979       ENDIF
2980       I=N
2981       NP=0
2982       NQ=0
2983       MBST=0
2984       IF(MMAT.GE.11.AND.MMAT.NE.46.AND.P(IP,4).GT.20.*P(IP,5)) MBST=1
2985       DO 160 J=1,4
2986       PV(1,J)=0.
2987   160 IF(MBST.EQ.0) PV(1,J)=P(IP,J)
2988       IF(MBST.EQ.1) PV(1,4)=P(IP,5)
2989       PV(1,5)=P(IP,5)
2990       PS=0.
2991       PSQ=0.
2992       MREM=0
2993
2994 C...Read out decay products. Convert to standard flavour code.
2995       JTMAX=5
2996       IF(MDME(IDC+1,2).EQ.101) JTMAX=10
2997       DO 170 JT=1,JTMAX
2998       IF(JT.LE.5) KP=KFDP(IDC,JT)
2999       IF(JT.GE.6) KP=KFDP(IDC+1,JT-5)
3000       IF(KP.EQ.0) GOTO 170
3001       KPA=IABS(KP)
3002       KCP=LUCOMP(KPA)
3003       IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN
3004         KFP=KP
3005       ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN
3006         KFP=KFS*KP
3007       ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN
3008         KFP=-KFS*MOD(KFA/10,10)
3009       ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN
3010         KFP=KFS*(100*MOD(KFA/10,100)+3)
3011       ELSEIF(KPA.EQ.81) THEN
3012         KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1)
3013       ELSEIF(KP.EQ.82) THEN
3014         CALL LUKFDI(-KFS*INT(1.+(2.+PARJ(2))*RLU(0)),0,KFP,KDUMP)
3015         IF(KFP.EQ.0) GOTO 150
3016         MSTJ(93)=1
3017         IF(PV(1,5).LT.PARJ(32)+2.*ULMASS(KFP)) GOTO 150
3018       ELSEIF(KP.EQ.-82) THEN
3019         KFP=-KFP
3020         IF(IABS(KFP).GT.10) KFP=KFP+ISIGN(10000,KFP)
3021       ENDIF
3022       IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=LUCOMP(KFP)
3023
3024 C...Add decay product to event record or to quark flavour list.
3025       KFPA=IABS(KFP)
3026       KQP=KCHG(KCP,2)
3027       IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN
3028         NQ=NQ+1
3029         KFLO(NQ)=KFP
3030         MSTJ(93)=2
3031         PSQ=PSQ+ULMASS(KFLO(NQ))
3032       ELSEIF(MMAT.GE.42.AND.MMAT.LE.43.AND.NP.EQ.3.AND.MOD(NQ,2).EQ.1)
3033      &THEN
3034         NQ=NQ-1
3035         PS=PS-P(I,5)
3036         K(I,1)=1
3037         KFI=K(I,2)
3038         CALL LUKFDI(KFP,KFI,KFLDMP,K(I,2))
3039         IF(K(I,2).EQ.0) GOTO 150
3040         MSTJ(93)=1
3041         P(I,5)=ULMASS(K(I,2))
3042         PS=PS+P(I,5)
3043       ELSE
3044         I=I+1
3045         NP=NP+1
3046         IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1
3047         IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1
3048         K(I,1)=1+MOD(NQ,2)
3049         IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2
3050         IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1
3051         K(I,2)=KFP
3052         K(I,3)=IP
3053         K(I,4)=0
3054         K(I,5)=0
3055         P(I,5)=ULMASS(KFP)
3056         IF(MMAT.EQ.45.AND.KFPA.EQ.89) P(I,5)=PARJ(32)
3057         PS=PS+P(I,5)
3058       ENDIF
3059   170 CONTINUE
3060
3061 C...Choose decay multiplicity in phase space model.
3062   180 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN
3063         PSP=PS
3064         CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1))
3065         IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63)
3066   190   NTRY=NTRY+1
3067         IF(NTRY.GT.1000) THEN
3068           CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
3069           IF(MSTU(21).GE.1) RETURN
3070         ENDIF
3071         IF(MMAT.LE.20) THEN
3072           GAUSS=SQRT(-2.*CNDE*LOG(MAX(1E-10,RLU(0))))*
3073      &    SIN(PARU(2)*RLU(0))
3074           ND=0.5+0.5*NP+0.25*NQ+CNDE+GAUSS
3075           IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 190
3076           IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 190
3077           IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 190
3078           IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 190
3079         ELSE
3080           ND=MMAT-20
3081         ENDIF
3082
3083 C...Form hadrons from flavour content.
3084         DO 200 JT=1,4
3085   200   KFL1(JT)=KFLO(JT)
3086         IF(ND.EQ.NP+NQ/2) GOTO 220
3087         DO 210 I=N+NP+1,N+ND-NQ/2
3088         JT=1+INT((NQ-1)*RLU(0))
3089         CALL LUKFDI(KFL1(JT),0,KFL2,K(I,2))
3090         IF(K(I,2).EQ.0) GOTO 190
3091   210   KFL1(JT)=-KFL2
3092   220   JT=2
3093         JT2=3
3094         JT3=4
3095         IF(NQ.EQ.4.AND.RLU(0).LT.PARJ(66)) JT=4
3096         IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))*
3097      &  ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3
3098         IF(JT.EQ.3) JT2=2
3099         IF(JT.EQ.4) JT3=2
3100         CALL LUKFDI(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2))
3101         IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 190
3102         IF(NQ.EQ.4) CALL LUKFDI(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND,2))
3103         IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 190
3104
3105 C...Check that sum of decay product masses not too large.
3106         PS=PSP
3107         DO 230 I=N+NP+1,N+ND
3108         K(I,1)=1
3109         K(I,3)=IP
3110         K(I,4)=0
3111         K(I,5)=0
3112         P(I,5)=ULMASS(K(I,2))
3113   230   PS=PS+P(I,5)
3114         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 190
3115
3116 C...Rescale energy to subtract off spectator quark mass.
3117       ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44.OR.MMAT.EQ.45).
3118      &AND.NP.GE.3) THEN
3119         PS=PS-P(N+NP,5)
3120         PQT=(P(N+NP,5)+PARJ(65))/PV(1,5)
3121         DO 240 J=1,5
3122         P(N+NP,J)=PQT*PV(1,J)
3123   240   PV(1,J)=(1.-PQT)*PV(1,J)
3124         IF(PS+PARJ(64).GT.PV(1,5)) GOTO 150
3125         ND=NP-1
3126         MREM=1
3127
3128 C...Phase space factors imposed in W decay.
3129       ELSEIF(MMAT.EQ.46) THEN
3130         MSTJ(93)=1
3131         PSMC=ULMASS(K(N+1,2))
3132         MSTJ(93)=1
3133         PSMC=PSMC+ULMASS(K(N+2,2))
3134         IF(MAX(PS,PSMC)+PARJ(32).GT.PV(1,5)) GOTO 130
3135         HR1=(P(N+1,5)/PV(1,5))**2
3136         HR2=(P(N+2,5)/PV(1,5))**2
3137         IF((1.-HR1-HR2)*(2.+HR1+HR2)*SQRT((1.-HR1-HR2)**2-4.*HR1*HR2).
3138      &  LT.2.*RLU(0)) GOTO 130
3139         ND=NP
3140
3141 C...Fully specified final state: check mass broadening effects.
3142       ELSE
3143         IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 150
3144         ND=NP
3145       ENDIF
3146
3147 C...Select W mass in decay Q -> W + q, without W propagator.
3148       IF(MMAT.EQ.45.AND.MSTJ(25).LE.0) THEN
3149         HLQ=(PARJ(32)/PV(1,5))**2
3150         HUQ=(1.-(P(N+2,5)+PARJ(64))/PV(1,5))**2
3151         HRQ=(P(N+2,5)/PV(1,5))**2
3152   250   HW=HLQ+RLU(0)*(HUQ-HLQ)
3153         IF(HMEPS(HW).LT.RLU(0)) GOTO 250
3154         P(N+1,5)=PV(1,5)*SQRT(HW)
3155
3156 C...Ditto, including W propagator. Divide mass range into three regions.
3157       ELSEIF(MMAT.EQ.45) THEN
3158         HQW=(PV(1,5)/PMAS(24,1))**2
3159         HLW=(PARJ(32)/PMAS(24,1))**2
3160         HUW=((PV(1,5)-P(N+2,5)-PARJ(64))/PMAS(24,1))**2
3161         HRQ=(P(N+2,5)/PV(1,5))**2
3162         HG=PMAS(24,2)/PMAS(24,1)
3163         HATL=ATAN((HLW-1.)/HG)
3164         HM=MIN(1.,HUW-0.001)
3165         HMV1=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
3166   260   HM=HM-HG
3167         HMV2=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
3168         IF(HMV2.GT.HMV1.AND.HM-HG.GT.HLW) THEN
3169           HMV1=HMV2
3170           GOTO 260
3171         ENDIF
3172         HMV=MIN(2.*HMV1,HMEPS(HM/HQW)/HG**2)
3173         HM1=1.-SQRT(1./HMV-HG**2)
3174         IF(HM1.GT.HLW.AND.HM1.LT.HM) THEN
3175           HM=HM1
3176         ELSEIF(HMV2.LE.HMV1) THEN
3177           HM=MAX(HLW,HM-MIN(0.1,1.-HM))
3178         ENDIF
3179         HATM=ATAN((HM-1.)/HG)
3180         HWT1=(HATM-HATL)/HG
3181         HWT2=HMV*(MIN(1.,HUW)-HM)
3182         HWT3=0.
3183         IF(HUW.GT.1.) THEN
3184           HATU=ATAN((HUW-1.)/HG)
3185           HMP1=HMEPS(1./HQW)
3186           HWT3=HMP1*HATU/HG
3187         ENDIF
3188
3189 C...Select mass region and W mass there. Accept according to weight.
3190   270   HREG=RLU(0)*(HWT1+HWT2+HWT3)
3191         IF(HREG.LE.HWT1) THEN
3192           HW=1.+HG*TAN(HATL+RLU(0)*(HATM-HATL))
3193           HACC=HMEPS(HW/HQW)
3194         ELSEIF(HREG.LE.HWT1+HWT2) THEN
3195           HW=HM+RLU(0)*(MIN(1.,HUW)-HM)
3196           HACC=HMEPS(HW/HQW)/((HW-1.)**2+HG**2)/HMV
3197         ELSE
3198           HW=1.+HG*TAN(RLU(0)*HATU)
3199           HACC=HMEPS(HW/HQW)/HMP1
3200         ENDIF
3201         IF(HACC.LT.RLU(0)) GOTO 270
3202         P(N+1,5)=PMAS(24,1)*SQRT(HW)
3203       ENDIF
3204
3205 C...Determine position of grandmother, number of sisters, Q -> W sign.
3206       NM=0
3207       KFAS=0
3208       MSGN=0
3209       IF(MMAT.EQ.3.OR.MMAT.EQ.46) THEN
3210         IM=K(IP,3)
3211         IF(IM.LT.0.OR.IM.GE.IP) IM=0
3212         IF(IM.NE.0) KFAM=IABS(K(IM,2))
3213         IF(IM.NE.0.AND.MMAT.EQ.3) THEN
3214           DO 280 IL=MAX(IP-2,IM+1),MIN(IP+2,N)
3215           IF(K(IL,3).EQ.IM) NM=NM+1
3216   280     IF(K(IL,3).EQ.IM.AND.IL.NE.IP) ISIS=IL
3217           IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR.
3218      &    MOD(KFAM/1000,10).NE.0) NM=0
3219           IF(NM.EQ.2) THEN
3220             KFAS=IABS(K(ISIS,2))
3221             IF((KFAS.LE.100.OR.MOD(KFAS,10).NE.1.OR.
3222      &      MOD(KFAS/1000,10).NE.0).AND.KFAS.NE.22) NM=0
3223           ENDIF
3224         ELSEIF(IM.NE.0.AND.MMAT.EQ.46) THEN
3225           MSGN=ISIGN(1,K(IM,2)*K(IP,2))
3226           IF(KFAM.GT.100.AND.MOD(KFAM/1000,10).EQ.0) MSGN=
3227      &    MSGN*(-1)**MOD(KFAM/100,10)
3228         ENDIF
3229       ENDIF
3230
3231 C...Kinematics of one-particle decays.
3232       IF(ND.EQ.1) THEN
3233         DO 290 J=1,4
3234   290   P(N+1,J)=P(IP,J)
3235         GOTO 520
3236       ENDIF
3237
3238 C...Calculate maximum weight ND-particle decay.
3239       PV(ND,5)=P(N+ND,5)
3240       IF(ND.GE.3) THEN
3241         WTMAX=1./WTCOR(ND-2)
3242         PMAX=PV(1,5)-PS+P(N+ND,5)
3243         PMIN=0.
3244         DO 300 IL=ND-1,1,-1
3245         PMAX=PMAX+P(N+IL,5)
3246         PMIN=PMIN+P(N+IL+1,5)
3247   300   WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5))
3248       ENDIF
3249
3250 C...Find virtual gamma mass in Dalitz decay.
3251   310 IF(ND.EQ.2) THEN
3252       ELSEIF(MMAT.EQ.2) THEN
3253         PMES=4.*PMAS(11,1)**2
3254         PMRHO2=PMAS(131,1)**2
3255         PGRHO2=PMAS(131,2)**2
3256   320   PMST=PMES*(P(IP,5)**2/PMES)**RLU(0)
3257         WT=(1+0.5*PMES/PMST)*SQRT(MAX(0.,1.-PMES/PMST))*
3258      &  (1.-PMST/P(IP,5)**2)**3*(1.+PGRHO2/PMRHO2)/
3259      &  ((1.-PMST/PMRHO2)**2+PGRHO2/PMRHO2)
3260         IF(WT.LT.RLU(0)) GOTO 320
3261         PV(2,5)=MAX(2.00001*PMAS(11,1),SQRT(PMST))
3262
3263 C...M-generator gives weight. If rejected, try again.
3264       ELSE
3265   330   RORD(1)=1.
3266         DO 350 IL1=2,ND-1
3267         RSAV=RLU(0)
3268         DO 340 IL2=IL1-1,1,-1
3269         IF(RSAV.LE.RORD(IL2)) GOTO 350
3270   340   RORD(IL2+1)=RORD(IL2)
3271   350   RORD(IL2+1)=RSAV
3272         RORD(ND)=0.
3273         WT=1.
3274         DO 360 IL=ND-1,1,-1
3275         PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
3276   360   WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
3277         IF(WT.LT.RLU(0)*WTMAX) GOTO 330
3278       ENDIF
3279
3280 C...Perform two-particle decays in respective CM frame.
3281   370 DO 390 IL=1,ND-1
3282       PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
3283       UE(3)=2.*RLU(0)-1.
3284       PHI=PARU(2)*RLU(0)
3285       UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
3286       UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
3287       DO 380 J=1,3
3288       P(N+IL,J)=PA*UE(J)
3289   380 PV(IL+1,J)=-PA*UE(J)
3290       P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2)
3291   390 PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
3292
3293 C...Lorentz transform decay products to lab frame.
3294       DO 400 J=1,4
3295   400 P(N+ND,J)=PV(ND,J)
3296       DO 430 IL=ND-1,1,-1
3297       DO 410 J=1,3
3298   410 BE(J)=PV(IL,J)/PV(IL,4)
3299       GA=PV(IL,4)/PV(IL,5)
3300       DO 430 I=N+IL,N+ND
3301       BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
3302       DO 420 J=1,3
3303   420 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
3304   430 P(I,4)=GA*(P(I,4)+BEP)
3305
3306 C...Check that no infinite loop in matrix element weight.
3307       NTRY=NTRY+1
3308       IF(NTRY.GT.800) GOTO 450
3309
3310 C...Matrix elements for omega and phi decays.
3311       IF(MMAT.EQ.1) THEN
3312         WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2
3313      &  -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2
3314      &  +2.*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3)
3315         IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001).LT.RLU(0)) GOTO 310
3316
3317 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
3318       ELSEIF(MMAT.EQ.2) THEN
3319         FOUR12=FOUR(N+1,N+2)
3320         FOUR13=FOUR(N+1,N+3)
3321         WT=(PMST-0.5*PMES)*(FOUR12**2+FOUR13**2)+
3322      &  PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2)
3323         IF(WT.LT.RLU(0)*0.25*PMST*(P(IP,5)**2-PMST)**2) GOTO 370
3324
3325 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
3326 C...V vector), of form cos**2(theta02) in V1 rest frame, and for
3327 C...S0 -> gamma + V1 -> gamma + S2 + S3, of form sin**2(theta02).
3328       ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN
3329         FOUR10=FOUR(IP,IM)
3330         FOUR12=FOUR(IP,N+1)
3331         FOUR02=FOUR(IM,N+1)
3332         PMS1=P(IP,5)**2
3333         PMS0=P(IM,5)**2
3334         PMS2=P(N+1,5)**2
3335         IF(KFAS.NE.22) HNUM=(FOUR10*FOUR12-PMS1*FOUR02)**2
3336         IF(KFAS.EQ.22) HNUM=PMS1*(2.*FOUR10*FOUR12*FOUR02-
3337      &  PMS1*FOUR02**2-PMS0*FOUR12**2-PMS2*FOUR10**2+PMS1*PMS0*PMS2)
3338         HNUM=MAX(1E-6*PMS1**2*PMS0*PMS2,HNUM)
3339         HDEN=(FOUR10**2-PMS1*PMS0)*(FOUR12**2-PMS1*PMS2)
3340         IF(HNUM.LT.RLU(0)*HDEN) GOTO 370
3341
3342 C...Matrix element for "onium" -> g + g + g or gamma + g + g.
3343       ELSEIF(MMAT.EQ.4) THEN
3344         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
3345         HX2=2.*FOUR(IP,N+2)/P(IP,5)**2
3346         HX3=2.*FOUR(IP,N+3)/P(IP,5)**2
3347         WT=((1.-HX1)/(HX2*HX3))**2+((1.-HX2)/(HX1*HX3))**2+
3348      &  ((1.-HX3)/(HX1*HX2))**2
3349         IF(WT.LT.2.*RLU(0)) GOTO 310
3350         IF(K(IP+1,2).EQ.22.AND.(1.-HX1)*P(IP,5)**2.LT.4.*PARJ(32)**2)
3351      &  GOTO 310
3352
3353 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
3354       ELSEIF(MMAT.EQ.41) THEN
3355         HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
3356         IF(8.*HX1*(3.-2.*HX1)/9..LT.RLU(0)) GOTO 310
3357
3358 C...Matrix elements for weak decays (only semileptonic for c and b)
3359       ELSEIF(MMAT.GE.42.AND.MMAT.LE.44.AND.ND.EQ.3) THEN
3360         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3)
3361         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3)
3362         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310
3363       ELSEIF(MMAT.GE.42.AND.MMAT.LE.44) THEN
3364         DO 440 J=1,4
3365         P(N+NP+1,J)=0.
3366         DO 440 IS=N+3,N+NP
3367   440   P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J)
3368         IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1)
3369         IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1)
3370         IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310
3371
3372 C...Angular distribution in W decay.
3373       ELSEIF(MMAT.EQ.46.AND.MSGN.NE.0) THEN
3374         IF(MSGN.GT.0) WT=FOUR(IM,N+1)*FOUR(N+2,IP+1)
3375         IF(MSGN.LT.0) WT=FOUR(IM,N+2)*FOUR(N+1,IP+1)
3376         IF(WT.LT.RLU(0)*P(IM,5)**4/WTCOR(10)) GOTO 370
3377       ENDIF
3378
3379 C...Scale back energy and reattach spectator.
3380   450 IF(MREM.EQ.1) THEN
3381         DO 460 J=1,5
3382   460   PV(1,J)=PV(1,J)/(1.-PQT)
3383         ND=ND+1
3384         MREM=0
3385       ENDIF
3386
3387 C...Low invariant mass for system with spectator quark gives particle,
3388 C...not two jets. Readjust momenta accordingly.
3389       IF((MMAT.EQ.31.OR.MMAT.EQ.45).AND.ND.EQ.3) THEN
3390         MSTJ(93)=1
3391         PM2=ULMASS(K(N+2,2))
3392         MSTJ(93)=1
3393         PM3=ULMASS(K(N+3,2))
3394         IF(P(N+2,5)**2+P(N+3,5)**2+2.*FOUR(N+2,N+3).GE.
3395      &  (PARJ(32)+PM2+PM3)**2) GOTO 520
3396         K(N+2,1)=1
3397         KFTEMP=K(N+2,2)
3398         CALL LUKFDI(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2))
3399         IF(K(N+2,2).EQ.0) GOTO 150
3400         P(N+2,5)=ULMASS(K(N+2,2))
3401         PS=P(N+1,5)+P(N+2,5)
3402         PV(2,5)=P(N+2,5)
3403         MMAT=0
3404         ND=2
3405         GOTO 370
3406       ELSEIF(MMAT.EQ.44) THEN
3407         MSTJ(93)=1
3408         PM3=ULMASS(K(N+3,2))
3409         MSTJ(93)=1
3410         PM4=ULMASS(K(N+4,2))
3411         IF(P(N+3,5)**2+P(N+4,5)**2+2.*FOUR(N+3,N+4).GE.
3412      &  (PARJ(32)+PM3+PM4)**2) GOTO 490
3413         K(N+3,1)=1
3414         KFTEMP=K(N+3,2)
3415         CALL LUKFDI(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2))
3416         IF(K(N+3,2).EQ.0) GOTO 150
3417         P(N+3,5)=ULMASS(K(N+3,2))
3418         DO 470 J=1,3
3419   470   P(N+3,J)=P(N+3,J)+P(N+4,J)
3420         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)
3421         HA=P(N+1,4)**2-P(N+2,4)**2
3422         HB=HA-(P(N+1,5)**2-P(N+2,5)**2)
3423         HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+
3424      &  (P(N+1,3)-P(N+2,3))**2
3425         HD=(PV(1,4)-P(N+3,4))**2
3426         HE=HA**2-2.*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2
3427         HF=HD*HC-HB**2
3428         HG=HD*HC-HA*HB
3429         HH=(SQRT(HG**2+HE*HF)-HG)/(2.*HF)
3430         DO 480 J=1,3
3431         PCOR=HH*(P(N+1,J)-P(N+2,J))
3432         P(N+1,J)=P(N+1,J)+PCOR
3433   480   P(N+2,J)=P(N+2,J)-PCOR
3434         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)
3435         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)
3436         ND=ND-1
3437       ENDIF
3438
3439 C...Check invariant mass of W jets. May give one particle or start over.
3440   490 IF(MMAT.GE.42.AND.MMAT.LE.44.AND.IABS(K(N+1,2)).LT.10) THEN
3441         PMR=SQRT(MAX(0.,P(N+1,5)**2+P(N+2,5)**2+2.*FOUR(N+1,N+2)))
3442         MSTJ(93)=1
3443         PM1=ULMASS(K(N+1,2))
3444         MSTJ(93)=1
3445         PM2=ULMASS(K(N+2,2))
3446         IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 500
3447         KFLDUM=INT(1.5+RLU(0))
3448         CALL LUKFDI(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1)
3449         CALL LUKFDI(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2)
3450         IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 150
3451         PSM=ULMASS(KF1)+ULMASS(KF2)
3452         IF(MMAT.EQ.42.AND.PMR.GT.PARJ(64)+PSM) GOTO 500
3453         IF(MMAT.GE.43.AND.PMR.GT.0.2*PARJ(32)+PSM) GOTO 500
3454         IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 150
3455         K(N+1,1)=1
3456         KFTEMP=K(N+1,2)
3457         CALL LUKFDI(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2))
3458         IF(K(N+1,2).EQ.0) GOTO 150
3459         P(N+1,5)=ULMASS(K(N+1,2))
3460         K(N+2,2)=K(N+3,2)
3461         P(N+2,5)=P(N+3,5)
3462         PS=P(N+1,5)+P(N+2,5)
3463         PV(2,5)=P(N+3,5)
3464         MMAT=0
3465         ND=2
3466         GOTO 370
3467       ENDIF
3468
3469 C...Phase space decay of partons from W decay.
3470   500 IF(MMAT.EQ.42.AND.IABS(K(N+1,2)).LT.10) THEN
3471         KFLO(1)=K(N+1,2)
3472         KFLO(2)=K(N+2,2)
3473         K(N+1,1)=K(N+3,1)
3474         K(N+1,2)=K(N+3,2)