This commit was generated by cvs2svn to compensate for changes in r18048,
[u/mrichter/AliRoot.git] / THydjet / hydjet1_1 / jetset_73.f
CommitLineData
cb220f83 1C*********************************************************************
2CCPH This file has enlarged event record, LUJETS size=30000
3C*********************************************************************
4C*********************************************************************
5C*********************************************************************
6C* **
7C* June 1991 **
8C* **
9C* The Lund Monte Carlo for Jet Fragmentation and e+e- Physics **
10C* **
11C* JETSET version 7.3 **
12C* **
13C* Torbjorn Sjostrand **
14C* **
15C* CERN/TH, CH-1211 Geneva 23 **
16C* BITNET/EARN address TORSJO@CERNVM **
17C* Tel. +22 - 767 28 20 **
18C* **
19C* LUSHOW is written together with Mats Bengtsson **
20C* **
21C* A complete manual exists on a separate file **
22C* Please report any program errors to the author! **
23C* **
24C* Copyright Torbjorn Sjostrand **
25C* **
26C*********************************************************************
27C*********************************************************************
28C *
29C List of subprograms in order of appearance, with main purpose *
30C (S = subroutine, F = function, B = block data) *
31C *
32C S LU1ENT to fill one entry (= parton or particle) *
33C S LU2ENT to fill two entries *
34C S LU3ENT to fill three entries *
35C S LU4ENT to fill four entries *
36C S LUJOIN to connect entries with colour flow information *
37C S LUGIVE to fill (or query) commonblock variables *
38C S LUEXEC to administrate fragmentation and decay chain *
39C S LUPREP to rearrange showered partons along strings *
40C S LUSTRF to do string fragmentation of jet system *
41C S LUINDF to do independent fragmentation of one or many jets *
42C S LUDECY to do the decay of a particle *
43C S LUKFDI to select parton and hadron flavours in fragm *
44C S LUPTDI to select transverse momenta in fragm *
45C S LUZDIS to select longitudinal scaling variable in fragm *
46C S LUSHOW to do timelike parton shower evolution *
47C S LUBOEI to include Bose-Einstein effects (crudely) *
48C F ULMASS to give the mass of a particle or parton *
49C S LUNAME to give the name of a particle or parton *
50C F LUCHGE to give three times the electric charge *
51C F LUCOMP to compress standard KF flavour code to internal KC *
52C S LUERRM to write error messages and abort faulty run *
53C F ULALEM to give the alpha_electromagnetic value *
54C F ULALPS to give the alpha_strong value *
55C F ULANGL to give the angle from known x and y components *
56C F RLU to provide a random number generator *
57C S RLUGET to save the state of the random number generator *
58C S RLUSET to set the state of the random number generator *
59C S LUROBO to rotate and/or boost an event *
60C S LUEDIT to remove unwanted entries from record *
61C S LULIST to list event record or particle data *
62C S LUUPDA to update particle data *
63C F KLU to provide integer-valued event information *
64C F PLU to provide real-valued event information *
65C S LUSPHE to perform sphericity analysis *
66C S LUTHRU to perform thrust analysis *
67C S LUCLUS to perform three-dimensional cluster analysis *
68C S LUCELL to perform cluster analysis in (eta, phi, E_T) *
69C S LUJMAS to give high and low jet mass of event *
70C S LUFOWO to give Fox-Wolfram moments *
71C S LUTABU to analyze events, with tabular output *
72C *
73C S LUEEVT to administrate the generation of an e+e- event *
74C S LUXTOT to give the total cross-section at given CM energy *
75C S LURADK to generate initial state photon radiation *
76C S LUXKFL to select flavour of primary qqbar pair *
77C S LUXJET to select (matrix element) jet multiplicity *
78C S LUX3JT to select kinematics of three-jet event *
79C S LUX4JT to select kinematics of four-jet event *
80C S LUXDIF to select angular orientation of event *
81C S LUONIA to perform generation of onium decay to gluons *
82C *
83C S LUHEPC to convert between /LUJETS/ and /HEPEVT/ records *
84C S LUTEST to test the proper functioning of the package *
85C B LUDATA to contain default values and particle data *
86C *
87C*********************************************************************
88
89 SUBROUTINE LU1ENT(IP,KF,PE,THE,PHI)
90
91C...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
97C...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
106C...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
115C...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
126C...Set N. Optionally fragment/decay.
127 N=IPA
128 IF(IP.EQ.0) CALL LUEXEC
129
130 RETURN
131 END
132
133C*********************************************************************
134
135 SUBROUTINE LU2ENT(IP,KF1,KF2,PECM)
136
137C...Purpose: to store two partons/particles in their CM frame,
138C...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
144C...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
155C...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
168C...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
180C...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
186C...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
196C...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
208C...Set N. Optionally fragment/decay.
209 N=IPA+1
210 IF(IP.EQ.0) CALL LUEXEC
211
212 RETURN
213 END
214
215C*********************************************************************
216
217 SUBROUTINE LU3ENT(IP,KF1,KF2,KF3,PECM,X1,X3)
218
219C...Purpose: to store three partons or particles in their CM frame,
220C...with the first along the +z axis and the third in the (x,z)
221C...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
227C...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
239C...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
255C...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
271C...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
279C...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
294C...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
308C...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
321C...Set N. Optionally fragment/decay.
322 N=IPA+2
323 IF(IP.EQ.0) CALL LUEXEC
324
325 RETURN
326 END
327
328C*********************************************************************
329
330 SUBROUTINE LU4ENT(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)
331
332C...Purpose: to store four partons or particles in their CM frame, with
333C...the first along the +z axis, the last in the xz plane with x > 0
334C...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
340C...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
353C...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
372C...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
392C...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
403C...Store partons for parton shower evolution from q-g-g-qbar or
404C...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
421C...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
437C...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
460C...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
479C...Set N. Optionally fragment/decay.
480 N=IPA+3
481 IF(IP.EQ.0) CALL LUEXEC
482
483 RETURN
484 END
485
486C*********************************************************************
487
488 SUBROUTINE LUJOIN(NJOIN,IJOIN)
489
490C...Purpose: to connect a sequence of partons with colour flow indices,
491C...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
498C...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
514C...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
529C...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
537C*********************************************************************
538
539 SUBROUTINE LUGIVE(CHIN)
540
541C...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
566C...For each variable to be translated give: name,
567C...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
590C...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
609C...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
629C...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
681C...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
699C...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
786C...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
805C...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
820C...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
907C...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
929C...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
936C*********************************************************************
937
938 SUBROUTINE LUEXEC
939
940C...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
948C...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
958C...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
971C...Prepare system for subsequent fragmentation/decay.
972 CALL LUPREP(0)
973
974C...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
983C...Particle decay if unstable and allowed. Save long-lived particle
984C...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
990C...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
1005C...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
1021C...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
1029C...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
1035C...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
1052C*********************************************************************
1053
1054 SUBROUTINE LUPREP(IP)
1055
1056C...Purpose: to rearrange partons along strings, to allow small systems
1057C...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
1066C...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
1076C...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
1087C...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
1107C...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
1142C...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
1180C...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
1194C...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
1230C...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
1274C...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
1287C...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
1311C...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
1332C...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
1348C...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
1386C*********************************************************************
1387
1388 SUBROUTINE LUSTRF(IP)
1389C...Purpose: to handle the fragmentation of an arbitrary colour singlet
1390C...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
1400C...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
1405C...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
1431C...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
1452C...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
1475C...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
1502C...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
1528C...Reset particle counter. Skip ahead if no junctions are present;
1529C...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
1549C...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
1570C...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
1589C...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
1599C...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
1604C...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
1634C...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
1665C...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
1694C...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
1706C...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
1727C...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
1747C...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))
1785C...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
1794C...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
1808C...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
1817C...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
1829C...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
1848C...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
1872C...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
1880C...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
1895C...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
1921C...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
1951C...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
1982C...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
1996C...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
2026C...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
2065C...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
2074C...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
2091C...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
2102C...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
2115C...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
2137C...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
2157C...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))
2195C...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
2206C...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
2220C...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
2229C...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
2241C...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
2260C...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
2281C...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
2292C...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
2306C...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
2327C...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
2341C...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
2364C...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
2401C...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
2427C*********************************************************************
2428
2429 SUBROUTINE LUINDF(IP)
2430
2431C...Purpose: to handle the fragmentation of a jet system (or a single
2432C...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
2441C...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
2473C...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
2495C...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
2509C...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
2515C...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
2521C...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
2528C...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
2536C...Initial values for gluon treated like quark-antiquark jet pair,
2537C...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
2550C...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
2560C...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
2578C...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
2602C...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
2610C...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
2620C...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
2628C...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
2633C...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
2652C...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
2684C...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
2716C...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
2731C...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
2753C...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
2783C...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
2807C...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
2827C...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
2848C...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
2858C*********************************************************************
2859
2860 SUBROUTINE LUDECY(IP)
2861
2862C...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
2872C...Functions: momentum in two-particle decays, four-product and
2873C...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
2879C...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
2887C...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
2896C...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
2911C...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
2920C...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
2943C...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
2958C...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
2973C...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
2994C...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
3024C...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
3061C...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
3083C...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
3105C...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
3116C...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
3128C...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
3141C...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
3147C...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
3156C...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
3189C...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
3205C...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
3231C...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
3238C...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
3250C...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
3263C...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
3280C...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
3293C...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
3306C...Check that no infinite loop in matrix element weight.
3307 NTRY=NTRY+1
3308 IF(NTRY.GT.800) GOTO 450
3309
3310C...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
3317C...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
3325C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
3326C...V vector), of form cos**2(theta02) in V1 rest frame, and for
3327C...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
3342C...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
3353C...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
3358C...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
3372C...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
3379C...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
3387C...Low invariant mass for system with spectator quark gives particle,
3388C...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
3439C...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
3469C...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)
3475 DO 510 J=1,5
3476 PV(1,J)=P(N+1,J)+P(N+2,J)
3477 510 P(N+1,J)=P(N+3,J)
3478 PV(1,5)=PMR
3479 N=N+1
3480 NP=0
3481 NQ=2
3482 PS=0.
3483 MSTJ(93)=2
3484 PSQ=ULMASS(KFLO(1))
3485 MSTJ(93)=2
3486 PSQ=PSQ+ULMASS(KFLO(2))
3487 MMAT=11
3488 GOTO 180
3489 ENDIF
3490
3491C...Boost back for rapidly moving particle.
3492 520 N=N+ND
3493 IF(MBST.EQ.1) THEN
3494 DO 530 J=1,3
3495 530 BE(J)=P(IP,J)/P(IP,4)
3496 GA=P(IP,4)/P(IP,5)
3497 DO 550 I=NSAV+1,N
3498 BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
3499 DO 540 J=1,3
3500 540 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
3501 550 P(I,4)=GA*(P(I,4)+BEP)
3502 ENDIF
3503
3504C...Fill in position of decay vertex.
3505 DO 570 I=NSAV+1,N
3506 DO 560 J=1,4
3507 560 V(I,J)=VDCY(J)
3508 570 V(I,5)=0.
3509
3510C...Set up for parton shower evolution from jets.
3511 IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN
3512 K(NSAV+1,1)=3
3513 K(NSAV+2,1)=3
3514 K(NSAV+3,1)=3
3515 K(NSAV+1,4)=MSTU(5)*(NSAV+2)
3516 K(NSAV+1,5)=MSTU(5)*(NSAV+3)
3517 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
3518 K(NSAV+2,5)=MSTU(5)*(NSAV+1)
3519 K(NSAV+3,4)=MSTU(5)*(NSAV+1)
3520 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
3521 MSTJ(92)=-(NSAV+1)
3522 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN
3523 K(NSAV+2,1)=3
3524 K(NSAV+3,1)=3
3525 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
3526 K(NSAV+2,5)=MSTU(5)*(NSAV+3)
3527 K(NSAV+3,4)=MSTU(5)*(NSAV+2)
3528 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
3529 MSTJ(92)=NSAV+2
3530 ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46).
3531 &AND.IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN
3532 K(NSAV+1,1)=3
3533 K(NSAV+2,1)=3
3534 K(NSAV+1,4)=MSTU(5)*(NSAV+2)
3535 K(NSAV+1,5)=MSTU(5)*(NSAV+2)
3536 K(NSAV+2,4)=MSTU(5)*(NSAV+1)
3537 K(NSAV+2,5)=MSTU(5)*(NSAV+1)
3538 MSTJ(92)=NSAV+1
3539 ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46).
3540 &AND.IABS(K(NSAV+1,2)).LE.20.AND.IABS(K(NSAV+2,2)).LE.20) THEN
3541 MSTJ(92)=NSAV+1
3542 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21)
3543 &THEN
3544 K(NSAV+1,1)=3
3545 K(NSAV+2,1)=3
3546 K(NSAV+3,1)=3
3547 KCP=LUCOMP(K(NSAV+1,2))
3548 KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2))
3549 JCON=4
3550 IF(KQP.LT.0) JCON=5
3551 K(NSAV+1,JCON)=MSTU(5)*(NSAV+2)
3552 K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1)
3553 K(NSAV+2,JCON)=MSTU(5)*(NSAV+3)
3554 K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2)
3555 MSTJ(92)=NSAV+1
3556 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN
3557 K(NSAV+1,1)=3
3558 K(NSAV+3,1)=3
3559 K(NSAV+1,4)=MSTU(5)*(NSAV+3)
3560 K(NSAV+1,5)=MSTU(5)*(NSAV+3)
3561 K(NSAV+3,4)=MSTU(5)*(NSAV+1)
3562 K(NSAV+3,5)=MSTU(5)*(NSAV+1)
3563 MSTJ(92)=NSAV+1
3564 ENDIF
3565
3566C...Mark decayed particle; special option for B-B~ mixing.
3567 IF(K(IP,1).EQ.5) K(IP,1)=15
3568 IF(K(IP,1).LE.10) K(IP,1)=11
3569 IF(MMIX.EQ.1.AND.MSTJ(26).EQ.2.AND.K(IP,1).EQ.11) K(IP,1)=12
3570 K(IP,4)=NSAV+1
3571 K(IP,5)=N
3572
3573 RETURN
3574 END
3575
3576C*********************************************************************
3577
3578 SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF)
3579
3580C...Purpose: to generate a new flavour pair and combine off a hadron.
3581 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
3582 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
3583 SAVE /LUDAT1/,/LUDAT2/
3584
3585C...Default flavour values. Input consistency checks.
3586 KF1A=IABS(KFL1)
3587 KF2A=IABS(KFL2)
3588 KFL3=0
3589 KF=0
3590 IF(KF1A.EQ.0) RETURN
3591 IF(KF2A.NE.0) THEN
3592 IF(KF1A.LE.10.AND.KF2A.LE.10.AND.KFL1*KFL2.GT.0) RETURN
3593 IF(KF1A.GT.10.AND.KF2A.GT.10) RETURN
3594 IF((KF1A.GT.10.OR.KF2A.GT.10).AND.KFL1*KFL2.LT.0) RETURN
3595 ENDIF
3596
3597C...Check if tabulated flavour probabilities are to be used.
3598 IF(MSTJ(15).EQ.1) THEN
3599 KTAB1=-1
3600 IF(KF1A.GE.1.AND.KF1A.LE.6) KTAB1=KF1A
3601 KFL1A=MOD(KF1A/1000,10)
3602 KFL1B=MOD(KF1A/100,10)
3603 KFL1S=MOD(KF1A,10)
3604 IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1B.GE.1.AND.KFL1B.LE.4)
3605 & KTAB1=6+KFL1A*(KFL1A-2)+2*KFL1B+(KFL1S-1)/2
3606 IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1A.EQ.KFL1B) KTAB1=KTAB1-1
3607 IF(KF1A.GE.1.AND.KF1A.LE.6) KFL1A=KF1A
3608 KTAB2=0
3609 IF(KF2A.NE.0) THEN
3610 KTAB2=-1
3611 IF(KF2A.GE.1.AND.KF2A.LE.6) KTAB2=KF2A
3612 KFL2A=MOD(KF2A/1000,10)
3613 KFL2B=MOD(KF2A/100,10)
3614 KFL2S=MOD(KF2A,10)
3615 IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2B.GE.1.AND.KFL2B.LE.4)
3616 & KTAB2=6+KFL2A*(KFL2A-2)+2*KFL2B+(KFL2S-1)/2
3617 IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2A.EQ.KFL2B) KTAB2=KTAB2-1
3618 ENDIF
3619 IF(KTAB1.GE.0.AND.KTAB2.GE.0) GOTO 140
3620 ENDIF
3621
3622C...Parameters and breaking diquark parameter combinations.
3623 100 PAR2=PARJ(2)
3624 PAR3=PARJ(3)
3625 PAR4=3.*PARJ(4)
3626 IF(MSTJ(12).GE.2) THEN
3627 PAR3M=SQRT(PARJ(3))
3628 PAR4M=1./(3.*SQRT(PARJ(4)))
3629 PARDM=PARJ(7)/(PARJ(7)+PAR3M*PARJ(6))
3630 PARS0=PARJ(5)*(2.+(1.+PAR2*PAR3M*PARJ(7))*(1.+PAR4M))
3631 PARS1=PARJ(7)*PARS0/(2.*PAR3M)+PARJ(5)*(PARJ(6)*(1.+PAR4M)+
3632 & PAR2*PAR3M*PARJ(6)*PARJ(7))
3633 PARS2=PARJ(5)*2.*PARJ(6)*PARJ(7)*(PAR2*PARJ(7)+(1.+PAR4M)/PAR3M)
3634 PARSM=MAX(PARS0,PARS1,PARS2)
3635 PAR4=PAR4*(1.+PARSM)/(1.+PARSM/(3.*PAR4M))
3636 ENDIF
3637
3638C...Choice of whether to generate meson or baryon.
3639 MBARY=0
3640 KFDA=0
3641 IF(KF1A.LE.10) THEN
3642 IF(KF2A.EQ.0.AND.MSTJ(12).GE.1.AND.(1.+PARJ(1))*RLU(0).GT.1.)
3643 & MBARY=1
3644 IF(KF2A.GT.10) MBARY=2
3645 IF(KF2A.GT.10.AND.KF2A.LE.10000) KFDA=KF2A
3646 ELSE
3647 MBARY=2
3648 IF(KF1A.LE.10000) KFDA=KF1A
3649 ENDIF
3650
3651C...Possibility of process diquark -> meson + new diquark.
3652 IF(KFDA.NE.0.AND.MSTJ(12).GE.2) THEN
3653 KFLDA=MOD(KFDA/1000,10)
3654 KFLDB=MOD(KFDA/100,10)
3655 KFLDS=MOD(KFDA,10)
3656 WTDQ=PARS0
3657 IF(MAX(KFLDA,KFLDB).EQ.3) WTDQ=PARS1
3658 IF(MIN(KFLDA,KFLDB).EQ.3) WTDQ=PARS2
3659 IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
3660 IF((1.+WTDQ)*RLU(0).GT.1.) MBARY=-1
3661 IF(MBARY.EQ.-1.AND.KF2A.NE.0) RETURN
3662 ENDIF
3663
3664C...Flavour for meson, possibly with new flavour.
3665 IF(MBARY.LE.0) THEN
3666 KFS=ISIGN(1,KFL1)
3667 IF(MBARY.EQ.0) THEN
3668 IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),-KFL1)
3669 KFLA=MAX(KF1A,KF2A+IABS(KFL3))
3670 KFLB=MIN(KF1A,KF2A+IABS(KFL3))
3671 IF(KFLA.NE.KF1A) KFS=-KFS
3672
3673C...Splitting of diquark into meson plus new diquark.
3674 ELSE
3675 KFL1A=MOD(KF1A/1000,10)
3676 KFL1B=MOD(KF1A/100,10)
3677 110 KFL1D=KFL1A+INT(RLU(0)+0.5)*(KFL1B-KFL1A)
3678 KFL1E=KFL1A+KFL1B-KFL1D
3679 IF((KFL1D.EQ.3.AND.RLU(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.
3680 & RLU(0).LT.PARDM)) THEN
3681 KFL1D=KFL1A+KFL1B-KFL1D
3682 KFL1E=KFL1A+KFL1B-KFL1E
3683 ENDIF
3684 KFL3A=1+INT((2.+PAR2*PAR3M*PARJ(7))*RLU(0))
3685 IF((KFL1E.NE.KFL3A.AND.RLU(0).GT.(1.+PAR4M)/MAX(2.,1.+PAR4M)).
3686 & OR.(KFL1E.EQ.KFL3A.AND.RLU(0).GT.2./MAX(2.,1.+PAR4M)))
3687 & GOTO 110
3688 KFLDS=3
3689 IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU(0)+1./(1.+PAR4M))+1
3690 KFL3=ISIGN(10000+1000*MAX(KFL1E,KFL3A)+100*MIN(KFL1E,KFL3A)+
3691 & KFLDS,-KFL1)
3692 KFLA=MAX(KFL1D,KFL3A)
3693 KFLB=MIN(KFL1D,KFL3A)
3694 IF(KFLA.NE.KFL1D) KFS=-KFS
3695 ENDIF
3696
3697C...Form meson, with spin and flavour mixing for diagonal states.
3698 IF(KFLA.LE.2) KMUL=INT(PARJ(11)+RLU(0))
3699 IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU(0))
3700 IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU(0))
3701 IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN
3702 IF(RLU(0).LT.PARJ(14)) KMUL=2
3703 ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN
3704 RMUL=RLU(0)
3705 IF(RMUL.LT.PARJ(15)) KMUL=3
3706 IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)) KMUL=4
3707 IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)+PARJ(17)) KMUL=5
3708 ENDIF
3709 KFLS=3
3710 IF(KMUL.EQ.0.OR.KMUL.EQ.3) KFLS=1
3711 IF(KMUL.EQ.5) KFLS=5
3712 IF(KFLA.NE.KFLB) THEN
3713 KF=(100*KFLA+10*KFLB+KFLS)*KFS*(-1)**KFLA
3714 ELSE
3715 RMIX=RLU(0)
3716 IMIX=2*KFLA+10*KMUL
3717 IF(KFLA.LE.3) KF=110*(1+INT(RMIX+PARF(IMIX-1))+
3718 & INT(RMIX+PARF(IMIX)))+KFLS
3719 IF(KFLA.GE.4) KF=110*KFLA+KFLS
3720 ENDIF
3721 IF(KMUL.EQ.2.OR.KMUL.EQ.3) KF=KF+ISIGN(10000,KF)
3722 IF(KMUL.EQ.4) KF=KF+ISIGN(20000,KF)
3723
3724C...Generate diquark flavour.
3725 ELSE
3726 120 IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN
3727 KFLA=KF1A
3728 130 KFLB=1+INT((2.+PAR2*PAR3)*RLU(0))
3729 KFLC=1+INT((2.+PAR2*PAR3)*RLU(0))
3730 KFLDS=1
3731 IF(KFLB.GE.KFLC) KFLDS=3
3732 IF(KFLDS.EQ.1.AND.PAR4*RLU(0).GT.1.) GOTO 130
3733 IF(KFLDS.EQ.3.AND.PAR4.LT.RLU(0)) GOTO 130
3734 KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1)
3735
3736C...Take diquark flavour from input.
3737 ELSEIF(KF1A.LE.10) THEN
3738 KFLA=KF1A
3739 KFLB=MOD(KF2A/1000,10)
3740 KFLC=MOD(KF2A/100,10)
3741 KFLDS=MOD(KF2A,10)
3742
3743C...Generate (or take from input) quark to go with diquark.
3744 ELSE
3745 IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),KFL1)
3746 KFLA=KF2A+IABS(KFL3)
3747 KFLB=MOD(KF1A/1000,10)
3748 KFLC=MOD(KF1A/100,10)
3749 KFLDS=MOD(KF1A,10)
3750 ENDIF
3751