ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhiscat.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhiscat.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhiscat
5 
6 C...Finds outgoing flavours and event type; sets up the kinematics
7 C...and colour flow of the hard scattering.
8  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
9  SAVE /lujets/
10  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11  SAVE /ludat1/
12  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13  SAVE /ludat2/
14  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
15  SAVE /ludat3/
16  common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
17  SAVE /pyhisubs/
18  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
19  SAVE /pyhipars/
20  common/pyhiint1/mint(400),vint(400)
21  SAVE /pyhiint1/
22  common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
23  SAVE /pyhiint2/
24  common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
25  SAVE /pyhiint3/
26  common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3)
27  SAVE /pyhiint4/
28  common/pyhiint5/ngen(0:200,3),xsec(0:200,3)
29  SAVE /pyhiint5/
30  dimension wdtp(0:40),wdte(0:40,0:5),pmq(2),z(2),cthe(2),phi(2)
31 
32 C...Choice of subprocess, number of documentation lines.
33  isub=mint(1)
34  idoc=6+iset(isub)
35  IF(isub.EQ.95) idoc=8
36  mint(3)=idoc-6
37  IF(idoc.GE.9) idoc=idoc+2
38  mint(4)=idoc
39  ipu1=mint(84)+1
40  ipu2=mint(84)+2
41  ipu3=mint(84)+3
42  ipu4=mint(84)+4
43  ipu5=mint(84)+5
44  ipu6=mint(84)+6
45 
46 C...Reset K, P and V vectors. Store incoming particles.
47  DO 100 jt=1,mstp(126)+10
48  i=mint(83)+jt
49  DO 100 j=1,5
50  k(i,j)=0
51  p(i,j)=0.
52  100 v(i,j)=0.
53  DO 110 jt=1,2
54  i=mint(83)+jt
55  k(i,1)=21
56  k(i,2)=mint(10+jt)
57  p(i,1)=0.
58  p(i,2)=0.
59  p(i,5)=vint(2+jt)
60  p(i,3)=vint(5)*(-1)**(jt+1)
61  110 p(i,4)=sqrt(p(i,3)**2+p(i,5)**2)
62  mint(6)=2
63  kfres=0
64 
65 C...Store incoming partons in their CM-frame.
66  sh=vint(44)
67  shr=sqrt(sh)
68  shp=vint(26)*vint(2)
69  shpr=sqrt(shp)
70  shuser=shr
71  IF(iset(isub).GE.3) shuser=shpr
72  DO 120 jt=1,2
73  i=mint(84)+jt
74  k(i,1)=14
75  k(i,2)=mint(14+jt)
76  k(i,3)=mint(83)+2+jt
77  120 p(i,5)=ulmass(k(i,2))
78  IF(p(ipu1,5)+p(ipu2,5).GE.shuser) THEN
79  p(ipu1,5)=0.
80  p(ipu2,5)=0.
81  ENDIF
82  p(ipu1,4)=0.5*(shuser+(p(ipu1,5)**2-p(ipu2,5)**2)/shuser)
83  p(ipu1,3)=sqrt(max(0.,p(ipu1,4)**2-p(ipu1,5)**2))
84  p(ipu2,4)=shuser-p(ipu1,4)
85  p(ipu2,3)=-p(ipu1,3)
86 
87 C...Copy incoming partons to documentation lines.
88  DO 130 jt=1,2
89  i1=mint(83)+4+jt
90  i2=mint(84)+jt
91  k(i1,1)=21
92  k(i1,2)=k(i2,2)
93  k(i1,3)=i1-2
94  DO 130 j=1,5
95  130 p(i1,j)=p(i2,j)
96 
97 C...Choose new quark flavour for relevant annihilation graphs.
98  IF(isub.EQ.12.OR.isub.EQ.53) THEN
99  CALL pyhiwidt(21,shr,wdtp,wdte)
100  rkfl=(wdte(0,1)+wdte(0,2)+wdte(0,4))*rlu(0)
101  DO 140 i=1,2*mstp(1)
102  kflq=i
103  rkfl=rkfl-(wdte(i,1)+wdte(i,2)+wdte(i,4))
104  IF(rkfl.LE.0.) goto 150
105  140 CONTINUE
106  150 CONTINUE
107  ENDIF
108 
109 C...Final state flavours and colour flow: default values.
110  js=1
111  mint(21)=mint(15)
112  mint(22)=mint(16)
113  mint(23)=0
114  mint(24)=0
115  kcc=20
116  kcs=isign(1,mint(15))
117 
118  IF(isub.LE.10) THEN
119  IF(isub.EQ.1) THEN
120 C...f + fb -> gamma*/Z0.
121  kfres=23
122 
123  ELSEIF(isub.EQ.2) THEN
124 C...f + fb' -> W+/- .
125  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
126  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
127  kfres=isign(24,kch1+kch2)
128 
129  ELSEIF(isub.EQ.3) THEN
130 C...f + fb -> H0.
131  kfres=25
132 
133  ELSEIF(isub.EQ.4) THEN
134 C...gamma + W+/- -> W+/-.
135 
136  ELSEIF(isub.EQ.5) THEN
137 C...Z0 + Z0 -> H0.
138  xh=sh/shp
139  mint(21)=mint(15)
140  mint(22)=mint(16)
141  pmq(1)=ulmass(mint(21))
142  pmq(2)=ulmass(mint(22))
143  240 jt=int(1.5+rlu(0))
144  zmin=2.*pmq(jt)/shpr
145  zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))
146  zmax=min(1.-xh,zmax)
147  z(jt)=zmin+(zmax-zmin)*rlu(0)
148  IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.
149  & (1.-xh)**2/(4.*xh)*rlu(0)) goto 240
150  sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)
151  IF(sqc1.LT.1.e-8) goto 240
152  c1=sqrt(sqc1)
153  c2=1.+2.*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp)
154  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
155  cthe(jt)=min(1.,max(-1.,cthe(jt)))
156  z(3-jt)=1.-xh/(1.-z(jt))
157  sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)
158  IF(sqc1.LT.1.e-8) goto 240
159  c1=sqrt(sqc1)
160  c2=1.+2.*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
161  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
162  cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))
163  phir=paru(2)*rlu(0)
164  cphi=cos(phir)
165  ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi
166  z1=2.-z(jt)
167  z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp)
168  z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
169  z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
170  & pmq(3-jt)**2/shp))
171  zmin=2.*pmq(3-jt)/shpr
172  zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
173  zmax=min(1.-xh,zmax)
174  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 240
175  kcc=22
176  kfres=25
177 
178  ELSEIF(isub.EQ.6) THEN
179 C...Z0 + W+/- -> W+/-.
180 
181  ELSEIF(isub.EQ.7) THEN
182 C...W+ + W- -> Z0.
183 
184  ELSEIF(isub.EQ.8) THEN
185 C...W+ + W- -> H0.
186  xh=sh/shp
187  250 DO 280 jt=1,2
188  i=mint(14+jt)
189  ia=iabs(i)
190  IF(ia.LE.10) THEN
191  rvckm=vint(180+i)*rlu(0)
192  DO 270 j=1,mstp(1)
193  ib=2*j-1+mod(ia,2)
194  ipm=(5-isign(1,i))/2
195  idc=j+mdcy(ia,2)+2
196  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 270
197  mint(20+jt)=isign(ib,i)
198  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
199  IF(rvckm.LE.0.) goto 280
200  270 CONTINUE
201  ELSE
202  ib=2*((ia+1)/2)-1+mod(ia,2)
203  mint(20+jt)=isign(ib,i)
204  ENDIF
205  280 pmq(jt)=ulmass(mint(20+jt))
206  jt=int(1.5+rlu(0))
207  zmin=2.*pmq(jt)/shpr
208  zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))
209  zmax=min(1.-xh,zmax)
210  z(jt)=zmin+(zmax-zmin)*rlu(0)
211  IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.
212  & (1.-xh)**2/(4.*xh)*rlu(0)) goto 250
213  sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)
214  IF(sqc1.LT.1.e-8) goto 250
215  c1=sqrt(sqc1)
216  c2=1.+2.*(pmas(24,1)**2-pmq(jt)**2)/(z(jt)*shp)
217  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
218  cthe(jt)=min(1.,max(-1.,cthe(jt)))
219  z(3-jt)=1.-xh/(1.-z(jt))
220  sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)
221  IF(sqc1.LT.1.e-8) goto 250
222  c1=sqrt(sqc1)
223  c2=1.+2.*(pmas(24,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
224  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
225  cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))
226  phir=paru(2)*rlu(0)
227  cphi=cos(phir)
228  ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi
229  z1=2.-z(jt)
230  z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp)
231  z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
232  z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
233  & pmq(3-jt)**2/shp))
234  zmin=2.*pmq(3-jt)/shpr
235  zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
236  zmax=min(1.-xh,zmax)
237  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 250
238  kcc=22
239  kfres=25
240  ENDIF
241 
242  ELSEIF(isub.LE.20) THEN
243  IF(isub.EQ.11) THEN
244 C...f + f' -> f + f'; th = (p(f)-p(f))**2.
245  kcc=mint(2)
246  IF(mint(15)*mint(16).LT.0) kcc=kcc+2
247 
248  ELSEIF(isub.EQ.12) THEN
249 C...f + fb -> f' + fb'; th = (p(f)-p(f'))**2.
250  mint(21)=isign(kflq,mint(15))
251  mint(22)=-mint(21)
252  kcc=4
253 
254  ELSEIF(isub.EQ.13) THEN
255 C...f + fb -> g + g; th arbitrary.
256  mint(21)=21
257  mint(22)=21
258  kcc=mint(2)+4
259 
260  ELSEIF(isub.EQ.14) THEN
261 C...f + fb -> g + gam; th arbitrary.
262  IF(rlu(0).GT.0.5) js=2
263  mint(20+js)=21
264  mint(23-js)=22
265  kcc=17+js
266 
267  ELSEIF(isub.EQ.15) THEN
268 C...f + fb -> g + Z0; th arbitrary.
269  IF(rlu(0).GT.0.5) js=2
270  mint(20+js)=21
271  mint(23-js)=23
272  kcc=17+js
273 
274  ELSEIF(isub.EQ.16) THEN
275 C...f + fb' -> g + W+/-; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2.
276  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
277  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
278  IF(mint(15)*(kch1+kch2).LT.0) js=2
279  mint(20+js)=21
280  mint(23-js)=isign(24,kch1+kch2)
281  kcc=17+js
282 
283  ELSEIF(isub.EQ.17) THEN
284 C...f + fb -> g + H0; th arbitrary.
285  IF(rlu(0).GT.0.5) js=2
286  mint(20+js)=21
287  mint(23-js)=25
288  kcc=17+js
289 
290  ELSEIF(isub.EQ.18) THEN
291 C...f + fb -> gamma + gamma; th arbitrary.
292  mint(21)=22
293  mint(22)=22
294 
295  ELSEIF(isub.EQ.19) THEN
296 C...f + fb -> gamma + Z0; th arbitrary.
297  IF(rlu(0).GT.0.5) js=2
298  mint(20+js)=22
299  mint(23-js)=23
300 
301  ELSEIF(isub.EQ.20) THEN
302 C...f + fb' -> gamma + W+/-; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2.
303  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
304  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
305  IF(mint(15)*(kch1+kch2).LT.0) js=2
306  mint(20+js)=22
307  mint(23-js)=isign(24,kch1+kch2)
308  ENDIF
309 
310  ELSEIF(isub.LE.30) THEN
311  IF(isub.EQ.21) THEN
312 C...f + fb -> gamma + H0; th arbitrary.
313  IF(rlu(0).GT.0.5) js=2
314  mint(20+js)=22
315  mint(23-js)=25
316 
317  ELSEIF(isub.EQ.22) THEN
318 C...f + fb -> Z0 + Z0; th arbitrary.
319  mint(21)=23
320  mint(22)=23
321 
322  ELSEIF(isub.EQ.23) THEN
323 C...f + fb' -> Z0 + W+/-; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2.
324  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
325  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
326  IF(mint(15)*(kch1+kch2).LT.0) js=2
327  mint(20+js)=23
328  mint(23-js)=isign(24,kch1+kch2)
329 
330  ELSEIF(isub.EQ.24) THEN
331 C...f + fb -> Z0 + H0; th arbitrary.
332  IF(rlu(0).GT.0.5) js=2
333  mint(20+js)=23
334  mint(23-js)=25
335 
336  ELSEIF(isub.EQ.25) THEN
337 C...f + fb -> W+ + W-; th = (p(f)-p(W-))**2.
338  mint(21)=-isign(24,mint(15))
339  mint(22)=-mint(21)
340 
341  ELSEIF(isub.EQ.26) THEN
342 C...f + fb' -> W+/- + H0; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2.
343  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
344  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
345  IF(mint(15)*(kch1+kch2).GT.0) js=2
346  mint(20+js)=isign(24,kch1+kch2)
347  mint(23-js)=25
348 
349  ELSEIF(isub.EQ.27) THEN
350 C...f + fb -> H0 + H0.
351 
352  ELSEIF(isub.EQ.28) THEN
353 C...f + g -> f + g; th = (p(f)-p(f))**2.
354  kcc=mint(2)+6
355  IF(mint(15).EQ.21) kcc=kcc+2
356  IF(mint(15).NE.21) kcs=isign(1,mint(15))
357  IF(mint(16).NE.21) kcs=isign(1,mint(16))
358 
359  ELSEIF(isub.EQ.29) THEN
360 C...f + g -> f + gamma; th = (p(f)-p(f))**2.
361  IF(mint(15).EQ.21) js=2
362  mint(23-js)=22
363  kcc=15+js
364  kcs=isign(1,mint(14+js))
365 
366  ELSEIF(isub.EQ.30) THEN
367 C...f + g -> f + Z0; th = (p(f)-p(f))**2.
368  IF(mint(15).EQ.21) js=2
369  mint(23-js)=23
370  kcc=15+js
371  kcs=isign(1,mint(14+js))
372  ENDIF
373 
374  ELSEIF(isub.LE.40) THEN
375  IF(isub.EQ.31) THEN
376 C...f + g -> f' + W+/-; th = (p(f)-p(f'))**2; choose flavour f'.
377  IF(mint(15).EQ.21) js=2
378  i=mint(14+js)
379  ia=iabs(i)
380  mint(23-js)=isign(24,kchg(ia,1)*i)
381  rvckm=vint(180+i)*rlu(0)
382  DO 220 j=1,mstp(1)
383  ib=2*j-1+mod(ia,2)
384  ipm=(5-isign(1,i))/2
385  idc=j+mdcy(ia,2)+2
386  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 220
387  mint(20+js)=isign(ib,i)
388  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
389  IF(rvckm.LE.0.) goto 230
390  220 CONTINUE
391  230 kcc=15+js
392  kcs=isign(1,mint(14+js))
393 
394  ELSEIF(isub.EQ.32) THEN
395 C...f + g -> f + H0; th = (p(f)-p(f))**2.
396  IF(mint(15).EQ.21) js=2
397  mint(23-js)=25
398  kcc=15+js
399  kcs=isign(1,mint(14+js))
400 
401  ELSEIF(isub.EQ.33) THEN
402 C...f + gamma -> f + g.
403 
404  ELSEIF(isub.EQ.34) THEN
405 C...f + gamma -> f + gamma.
406 
407  ELSEIF(isub.EQ.35) THEN
408 C...f + gamma -> f + Z0.
409 
410  ELSEIF(isub.EQ.36) THEN
411 C...f + gamma -> f' + W+/-.
412 
413  ELSEIF(isub.EQ.37) THEN
414 C...f + gamma -> f + H0.
415 
416  ELSEIF(isub.EQ.38) THEN
417 C...f + Z0 -> f + g.
418 
419  ELSEIF(isub.EQ.39) THEN
420 C...f + Z0 -> f + gamma.
421 
422  ELSEIF(isub.EQ.40) THEN
423 C...f + Z0 -> f + Z0.
424  ENDIF
425 
426  ELSEIF(isub.LE.50) THEN
427  IF(isub.EQ.41) THEN
428 C...f + Z0 -> f' + W+/-.
429 
430  ELSEIF(isub.EQ.42) THEN
431 C...f + Z0 -> f + H0.
432 
433  ELSEIF(isub.EQ.43) THEN
434 C...f + W+/- -> f' + g.
435 
436  ELSEIF(isub.EQ.44) THEN
437 C...f + W+/- -> f' + gamma.
438 
439  ELSEIF(isub.EQ.45) THEN
440 C...f + W+/- -> f' + Z0.
441 
442  ELSEIF(isub.EQ.46) THEN
443 C...f + W+/- -> f' + W+/-.
444 
445  ELSEIF(isub.EQ.47) THEN
446 C...f + W+/- -> f' + H0.
447 
448  ELSEIF(isub.EQ.48) THEN
449 C...f + H0 -> f + g.
450 
451  ELSEIF(isub.EQ.49) THEN
452 C...f + H0 -> f + gamma.
453 
454  ELSEIF(isub.EQ.50) THEN
455 C...f + H0 -> f + Z0.
456  ENDIF
457 
458  ELSEIF(isub.LE.60) THEN
459  IF(isub.EQ.51) THEN
460 C...f + H0 -> f' + W+/-.
461 
462  ELSEIF(isub.EQ.52) THEN
463 C...f + H0 -> f + H0.
464 
465  ELSEIF(isub.EQ.53) THEN
466 C...g + g -> f + fb; th arbitrary.
467  kcs=(-1)**int(1.5+rlu(0))
468  mint(21)=isign(kflq,kcs)
469  mint(22)=-mint(21)
470  kcc=mint(2)+10
471 
472  ELSEIF(isub.EQ.54) THEN
473 C...g + gamma -> f + fb.
474 
475  ELSEIF(isub.EQ.55) THEN
476 C...g + Z0 -> f + fb.
477 
478  ELSEIF(isub.EQ.56) THEN
479 C...g + W+/- -> f + fb'.
480 
481  ELSEIF(isub.EQ.57) THEN
482 C...g + H0 -> f + fb.
483 
484  ELSEIF(isub.EQ.58) THEN
485 C...gamma + gamma -> f + fb.
486 
487  ELSEIF(isub.EQ.59) THEN
488 C...gamma + Z0 -> f + fb.
489 
490  ELSEIF(isub.EQ.60) THEN
491 C...gamma + W+/- -> f + fb'.
492  ENDIF
493 
494  ELSEIF(isub.LE.70) THEN
495  IF(isub.EQ.61) THEN
496 C...gamma + H0 -> f + fb.
497 
498  ELSEIF(isub.EQ.62) THEN
499 C...Z0 + Z0 -> f + fb.
500 
501  ELSEIF(isub.EQ.63) THEN
502 C...Z0 + W+/- -> f + fb'.
503 
504  ELSEIF(isub.EQ.64) THEN
505 C...Z0 + H0 -> f + fb.
506 
507  ELSEIF(isub.EQ.65) THEN
508 C...W+ + W- -> f + fb.
509 
510  ELSEIF(isub.EQ.66) THEN
511 C...W+/- + H0 -> f + fb'.
512 
513  ELSEIF(isub.EQ.67) THEN
514 C...H0 + H0 -> f + fb.
515 
516  ELSEIF(isub.EQ.68) THEN
517 C...g + g -> g + g; th arbitrary.
518  kcc=mint(2)+12
519  kcs=(-1)**int(1.5+rlu(0))
520 
521  ELSEIF(isub.EQ.69) THEN
522 C...gamma + gamma -> W+ + W-.
523 
524  ELSEIF(isub.EQ.70) THEN
525 C...gamma + W+/- -> gamma + W+/-
526  ENDIF
527 
528  ELSEIF(isub.LE.80) THEN
529  IF(isub.EQ.71.OR.isub.EQ.72) THEN
530 C...Z0 + Z0 -> Z0 + Z0; Z0 + Z0 -> W+ + W-.
531  xh=sh/shp
532  mint(21)=mint(15)
533  mint(22)=mint(16)
534  pmq(1)=ulmass(mint(21))
535  pmq(2)=ulmass(mint(22))
536  290 jt=int(1.5+rlu(0))
537  zmin=2.*pmq(jt)/shpr
538  zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))
539  zmax=min(1.-xh,zmax)
540  z(jt)=zmin+(zmax-zmin)*rlu(0)
541  IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.
542  & (1.-xh)**2/(4.*xh)*rlu(0)) goto 290
543  sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)
544  IF(sqc1.LT.1.e-8) goto 290
545  c1=sqrt(sqc1)
546  c2=1.+2.*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp)
547  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
548  cthe(jt)=min(1.,max(-1.,cthe(jt)))
549  z(3-jt)=1.-xh/(1.-z(jt))
550  sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)
551  IF(sqc1.LT.1.e-8) goto 290
552  c1=sqrt(sqc1)
553  c2=1.+2.*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
554  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
555  cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))
556  phir=paru(2)*rlu(0)
557  cphi=cos(phir)
558  ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi
559  z1=2.-z(jt)
560  z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp)
561  z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
562  z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
563  & pmq(3-jt)**2/shp))
564  zmin=2.*pmq(3-jt)/shpr
565  zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
566  zmax=min(1.-xh,zmax)
567  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 290
568  kcc=22
569 
570  ELSEIF(isub.EQ.73) THEN
571 C...Z0 + W+/- -> Z0 + W+/-.
572  xh=sh/shp
573  300 jt=int(1.5+rlu(0))
574  i=mint(14+jt)
575  ia=iabs(i)
576  IF(ia.LE.10) THEN
577  rvckm=vint(180+i)*rlu(0)
578  DO 320 j=1,mstp(1)
579  ib=2*j-1+mod(ia,2)
580  ipm=(5-isign(1,i))/2
581  idc=j+mdcy(ia,2)+2
582  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 320
583  mint(20+jt)=isign(ib,i)
584  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
585  IF(rvckm.LE.0.) goto 330
586  320 CONTINUE
587  ELSE
588  ib=2*((ia+1)/2)-1+mod(ia,2)
589  mint(20+jt)=isign(ib,i)
590  ENDIF
591  330 pmq(jt)=ulmass(mint(20+jt))
592  mint(23-jt)=mint(17-jt)
593  pmq(3-jt)=ulmass(mint(23-jt))
594  jt=int(1.5+rlu(0))
595  zmin=2.*pmq(jt)/shpr
596  zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))
597  zmax=min(1.-xh,zmax)
598  z(jt)=zmin+(zmax-zmin)*rlu(0)
599  IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.
600  & (1.-xh)**2/(4.*xh)*rlu(0)) goto 300
601  sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)
602  IF(sqc1.LT.1.e-8) goto 300
603  c1=sqrt(sqc1)
604  c2=1.+2.*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp)
605  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
606  cthe(jt)=min(1.,max(-1.,cthe(jt)))
607  z(3-jt)=1.-xh/(1.-z(jt))
608  sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)
609  IF(sqc1.LT.1.e-8) goto 300
610  c1=sqrt(sqc1)
611  c2=1.+2.*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
612  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
613  cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))
614  phir=paru(2)*rlu(0)
615  cphi=cos(phir)
616  ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi
617  z1=2.-z(jt)
618  z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp)
619  z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
620  z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
621  & pmq(3-jt)**2/shp))
622  zmin=2.*pmq(3-jt)/shpr
623  zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
624  zmax=min(1.-xh,zmax)
625  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 300
626  kcc=22
627 
628  ELSEIF(isub.EQ.74) THEN
629 C...Z0 + H0 -> Z0 + H0.
630 
631  ELSEIF(isub.EQ.75) THEN
632 C...W+ + W- -> gamma + gamma.
633 
634  ELSEIF(isub.EQ.76.OR.isub.EQ.77) THEN
635 C...W+ + W- -> Z0 + Z0; W+ + W- -> W+ + W-.
636  xh=sh/shp
637  340 DO 370 jt=1,2
638  i=mint(14+jt)
639  ia=iabs(i)
640  IF(ia.LE.10) THEN
641  rvckm=vint(180+i)*rlu(0)
642  DO 360 j=1,mstp(1)
643  ib=2*j-1+mod(ia,2)
644  ipm=(5-isign(1,i))/2
645  idc=j+mdcy(ia,2)+2
646  IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 360
647  mint(20+jt)=isign(ib,i)
648  rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)
649  IF(rvckm.LE.0.) goto 370
650  360 CONTINUE
651  ELSE
652  ib=2*((ia+1)/2)-1+mod(ia,2)
653  mint(20+jt)=isign(ib,i)
654  ENDIF
655  370 pmq(jt)=ulmass(mint(20+jt))
656  jt=int(1.5+rlu(0))
657  zmin=2.*pmq(jt)/shpr
658  zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))
659  zmax=min(1.-xh,zmax)
660  z(jt)=zmin+(zmax-zmin)*rlu(0)
661  IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.
662  & (1.-xh)**2/(4.*xh)*rlu(0)) goto 340
663  sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)
664  IF(sqc1.LT.1.e-8) goto 340
665  c1=sqrt(sqc1)
666  c2=1.+2.*(pmas(24,1)**2-pmq(jt)**2)/(z(jt)*shp)
667  cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
668  cthe(jt)=min(1.,max(-1.,cthe(jt)))
669  z(3-jt)=1.-xh/(1.-z(jt))
670  sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)
671  IF(sqc1.LT.1.e-8) goto 340
672  c1=sqrt(sqc1)
673  c2=1.+2.*(pmas(24,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp)
674  cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1
675  cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))
676  phir=paru(2)*rlu(0)
677  cphi=cos(phir)
678  ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi
679  z1=2.-z(jt)
680  z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp)
681  z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp
682  z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*
683  & pmq(3-jt)**2/shp))
684  zmin=2.*pmq(3-jt)/shpr
685  zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))
686  zmax=min(1.-xh,zmax)
687  IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 340
688  kcc=22
689 
690  ELSEIF(isub.EQ.78) THEN
691 C...W+/- + H0 -> W+/- + H0.
692 
693  ELSEIF(isub.EQ.79) THEN
694 C...H0 + H0 -> H0 + H0.
695  ENDIF
696 
697  ELSEIF(isub.LE.90) THEN
698  IF(isub.EQ.81) THEN
699 C...q + qb -> Q' + Qb'; th = (p(q)-p(q'))**2.
700  mint(21)=isign(mint(46),mint(15))
701  mint(22)=-mint(21)
702  kcc=4
703 
704  ELSEIF(isub.EQ.82) THEN
705 C...g + g -> Q + Qb; th arbitrary.
706  kcs=(-1)**int(1.5+rlu(0))
707  mint(21)=isign(mint(46),kcs)
708  mint(22)=-mint(21)
709  kcc=mint(2)+10
710  ENDIF
711 
712  ELSEIF(isub.LE.100) THEN
713  IF(isub.EQ.95) THEN
714 C...Low-pT ( = energyless g + g -> g + g).
715  kcc=mint(2)+12
716  kcs=(-1)**int(1.5+rlu(0))
717 
718  ELSEIF(isub.EQ.96) THEN
719 C...Multiple interactions (should be reassigned to QCD process).
720  ENDIF
721 
722  ELSEIF(isub.LE.110) THEN
723  IF(isub.EQ.101) THEN
724 C...g + g -> gamma*/Z0.
725  kcc=21
726  kfres=22
727 
728  ELSEIF(isub.EQ.102) THEN
729 C...g + g -> H0.
730  kcc=21
731  kfres=25
732  ENDIF
733 
734  ELSEIF(isub.LE.120) THEN
735  IF(isub.EQ.111) THEN
736 C...f + fb -> g + H0; th arbitrary.
737  IF(rlu(0).GT.0.5) js=2
738  mint(20+js)=21
739  mint(23-js)=25
740  kcc=17+js
741 
742  ELSEIF(isub.EQ.112) THEN
743 C...f + g -> f + H0; th = (p(f) - p(f))**2.
744  IF(mint(15).EQ.21) js=2
745  mint(23-js)=25
746  kcc=15+js
747  kcs=isign(1,mint(14+js))
748 
749  ELSEIF(isub.EQ.113) THEN
750 C...g + g -> g + H0; th arbitrary.
751  IF(rlu(0).GT.0.5) js=2
752  mint(23-js)=25
753  kcc=22+js
754  kcs=(-1)**int(1.5+rlu(0))
755 
756  ELSEIF(isub.EQ.114) THEN
757 C...g + g -> gamma + gamma; th arbitrary.
758  IF(rlu(0).GT.0.5) js=2
759  mint(21)=22
760  mint(22)=22
761  kcc=21
762 
763  ELSEIF(isub.EQ.115) THEN
764 C...g + g -> gamma + Z0.
765 
766  ELSEIF(isub.EQ.116) THEN
767 C...g + g -> Z0 + Z0.
768 
769  ELSEIF(isub.EQ.117) THEN
770 C...g + g -> W+ + W-.
771  ENDIF
772 
773  ELSEIF(isub.LE.140) THEN
774  IF(isub.EQ.121) THEN
775 C...g + g -> f + fb + H0.
776  ENDIF
777 
778  ELSEIF(isub.LE.160) THEN
779  IF(isub.EQ.141) THEN
780 C...f + fb -> gamma*/Z0/Z'0.
781  kfres=32
782 
783  ELSEIF(isub.EQ.142) THEN
784 C...f + fb' -> H+/-.
785  kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))
786  kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))
787  kfres=isign(37,kch1+kch2)
788 
789  ELSEIF(isub.EQ.143) THEN
790 C...f + fb' -> R.
791  kfres=isign(40,mint(15)+mint(16))
792  ENDIF
793 
794  ELSE
795  IF(isub.EQ.161) THEN
796 C...g + f -> H+/- + f'; th = (p(f)-p(f))**2.
797  IF(mint(16).EQ.21) js=2
798  ia=iabs(mint(17-js))
799  mint(20+js)=isign(37,kchg(ia,1)*mint(17-js))
800  ja=ia+mod(ia,2)-mod(ia+1,2)
801  mint(23-js)=isign(ja,mint(17-js))
802  kcc=18-js
803  IF(mint(15).NE.21) kcs=isign(1,mint(15))
804  IF(mint(16).NE.21) kcs=isign(1,mint(16))
805  ENDIF
806  ENDIF
807 
808  IF(idoc.EQ.7) THEN
809 C...Resonance not decaying: store colour connection indices.
810  i=mint(83)+7
811  k(ipu3,1)=1
812  k(ipu3,2)=kfres
813  k(ipu3,3)=i
814  p(ipu3,4)=shuser
815  p(ipu3,5)=shuser
816  k(ipu1,4)=ipu2
817  k(ipu1,5)=ipu2
818  k(ipu2,4)=ipu1
819  k(ipu2,5)=ipu1
820  k(i,1)=21
821  k(i,2)=kfres
822  p(i,4)=shuser
823  p(i,5)=shuser
824  n=ipu3
825  mint(21)=kfres
826  mint(22)=0
827 
828  ELSEIF(idoc.EQ.8) THEN
829 C...2 -> 2 processes: store outgoing partons in their CM-frame.
830  DO 390 jt=1,2
831  i=mint(84)+2+jt
832  k(i,1)=1
833  IF(iabs(mint(20+jt)).LE.10.OR.mint(20+jt).EQ.21) k(i,1)=3
834  k(i,2)=mint(20+jt)
835  k(i,3)=mint(83)+idoc+jt-2
836  IF(iabs(k(i,2)).LE.10.OR.k(i,2).EQ.21) THEN
837  p(i,5)=ulmass(k(i,2))
838  ELSE
839  p(i,5)=sqrt(vint(63+mod(js+jt,2)))
840  ENDIF
841  390 CONTINUE
842  IF(p(ipu3,5)+p(ipu4,5).GE.shr) THEN
843  kfa1=iabs(mint(21))
844  kfa2=iabs(mint(22))
845  IF((kfa1.GT.3.AND.kfa1.NE.21).OR.(kfa2.GT.3.AND.kfa2.NE.21))
846  & THEN
847  mint(51)=1
848  RETURN
849  ENDIF
850  p(ipu3,5)=0.
851  p(ipu4,5)=0.
852  ENDIF
853  p(ipu3,4)=0.5*(shr+(p(ipu3,5)**2-p(ipu4,5)**2)/shr)
854  p(ipu3,3)=sqrt(max(0.,p(ipu3,4)**2-p(ipu3,5)**2))
855  p(ipu4,4)=shr-p(ipu3,4)
856  p(ipu4,3)=-p(ipu3,3)
857  n=ipu4
858  mint(7)=mint(83)+7
859  mint(8)=mint(83)+8
860 
861 C...Rotate outgoing partons using cos(theta)=(th-uh)/lam(sh,sqm3,sqm4).
862  CALL ludbrb(ipu3,ipu4,acos(vint(23)),vint(24),0d0,0d0,0d0)
863 
864  ELSEIF(idoc.EQ.9) THEN
865 C'''2 -> 3 processes:
866 
867  ELSEIF(idoc.EQ.11) THEN
868 C...Z0 + Z0 -> H0, W+ + W- -> H0: store Higgs and outgoing partons.
869  phi(1)=paru(2)*rlu(0)
870  phi(2)=phi(1)-phir
871  DO 400 jt=1,2
872  i=mint(84)+2+jt
873  k(i,1)=1
874  IF(iabs(mint(20+jt)).LE.10.OR.mint(20+jt).EQ.21) k(i,1)=3
875  k(i,2)=mint(20+jt)
876  k(i,3)=mint(83)+idoc+jt-2
877  p(i,5)=ulmass(k(i,2))
878  IF(0.5*shpr*z(jt).LE.p(i,5)) p(i,5)=0.
879  pabs=sqrt(max(0.,(0.5*shpr*z(jt))**2-p(i,5)**2))
880  ptabs=pabs*sqrt(max(0.,1.-cthe(jt)**2))
881  p(i,1)=ptabs*cos(phi(jt))
882  p(i,2)=ptabs*sin(phi(jt))
883  p(i,3)=pabs*cthe(jt)*(-1)**(jt+1)
884  p(i,4)=0.5*shpr*z(jt)
885  izw=mint(83)+6+jt
886  k(izw,1)=21
887  k(izw,2)=23
888  IF(isub.EQ.8) k(izw,2)=isign(24,luchge(mint(14+jt)))
889  k(izw,3)=izw-2
890  p(izw,1)=-p(i,1)
891  p(izw,2)=-p(i,2)
892  p(izw,3)=(0.5*shpr-pabs*cthe(jt))*(-1)**(jt+1)
893  p(izw,4)=0.5*shpr*(1.-z(jt))
894  400 p(izw,5)=-sqrt(max(0.,p(izw,3)**2+ptabs**2-p(izw,4)**2))
895  i=mint(83)+9
896  k(ipu5,1)=1
897  k(ipu5,2)=kfres
898  k(ipu5,3)=i
899  p(ipu5,5)=shr
900  p(ipu5,1)=-p(ipu3,1)-p(ipu4,1)
901  p(ipu5,2)=-p(ipu3,2)-p(ipu4,2)
902  p(ipu5,3)=-p(ipu3,3)-p(ipu4,3)
903  p(ipu5,4)=shpr-p(ipu3,4)-p(ipu4,4)
904  k(i,1)=21
905  k(i,2)=kfres
906  DO 410 j=1,5
907  410 p(i,j)=p(ipu5,j)
908  n=ipu5
909  mint(23)=kfres
910 
911  ELSEIF(idoc.EQ.12) THEN
912 C...Z0 and W+/- scattering: store bosons and outgoing partons.
913  phi(1)=paru(2)*rlu(0)
914  phi(2)=phi(1)-phir
915  DO 420 jt=1,2
916  i=mint(84)+2+jt
917  k(i,1)=1
918  IF(iabs(mint(20+jt)).LE.10.OR.mint(20+jt).EQ.21) k(i,1)=3
919  k(i,2)=mint(20+jt)
920  k(i,3)=mint(83)+idoc+jt-2
921  p(i,5)=ulmass(k(i,2))
922  IF(0.5*shpr*z(jt).LE.p(i,5)) p(i,5)=0.
923  pabs=sqrt(max(0.,(0.5*shpr*z(jt))**2-p(i,5)**2))
924  ptabs=pabs*sqrt(max(0.,1.-cthe(jt)**2))
925  p(i,1)=ptabs*cos(phi(jt))
926  p(i,2)=ptabs*sin(phi(jt))
927  p(i,3)=pabs*cthe(jt)*(-1)**(jt+1)
928  p(i,4)=0.5*shpr*z(jt)
929  izw=mint(83)+6+jt
930  k(izw,1)=21
931  IF(mint(14+jt).EQ.mint(20+jt)) THEN
932  k(izw,2)=23
933  ELSE
934  k(izw,2)=isign(24,luchge(mint(14+jt))-luchge(mint(20+jt)))
935  ENDIF
936  k(izw,3)=izw-2
937  p(izw,1)=-p(i,1)
938  p(izw,2)=-p(i,2)
939  p(izw,3)=(0.5*shpr-pabs*cthe(jt))*(-1)**(jt+1)
940  p(izw,4)=0.5*shpr*(1.-z(jt))
941  p(izw,5)=-sqrt(max(0.,p(izw,3)**2+ptabs**2-p(izw,4)**2))
942  ipu=mint(84)+4+jt
943  k(ipu,1)=3
944  k(ipu,2)=kfpr(isub,jt)
945  k(ipu,3)=mint(83)+8+jt
946  IF(iabs(k(ipu,2)).LE.10.OR.k(ipu,2).EQ.21) THEN
947  p(ipu,5)=ulmass(k(ipu,2))
948  ELSE
949  p(ipu,5)=sqrt(vint(63+mod(js+jt,2)))
950  ENDIF
951  mint(22+jt)=k(izw,2)
952  420 CONTINUE
953  IF(isub.EQ.72) k(mint(84)+4+int(1.5+rlu(0)),2)=-24
954 C...Find rotation and boost for hard scattering subsystem.
955  i1=mint(83)+7
956  i2=mint(83)+8
957  bexcm=(p(i1,1)+p(i2,1))/(p(i1,4)+p(i2,4))
958  beycm=(p(i1,2)+p(i2,2))/(p(i1,4)+p(i2,4))
959  bezcm=(p(i1,3)+p(i2,3))/(p(i1,4)+p(i2,4))
960  gamcm=(p(i1,4)+p(i2,4))/shr
961  bepcm=bexcm*p(i1,1)+beycm*p(i1,2)+bezcm*p(i1,3)
962  px=p(i1,1)+gamcm*(gamcm/(1.+gamcm)*bepcm-p(i1,4))*bexcm
963  py=p(i1,2)+gamcm*(gamcm/(1.+gamcm)*bepcm-p(i1,4))*beycm
964  pz=p(i1,3)+gamcm*(gamcm/(1.+gamcm)*bepcm-p(i1,4))*bezcm
965  thecm=ulangl(pz,sqrt(px**2+py**2))
966  phicm=ulangl(px,py)
967 C...Store hard scattering subsystem. Rotate and boost it.
968  sqlam=(sh-p(ipu5,5)**2-p(ipu6,5)**2)**2-4.*p(ipu5,5)**2*
969  & p(ipu6,5)**2
970  pabs=sqrt(max(0.,sqlam/(4.*sh)))
971  cthwz=vint(23)
972  sthwz=sqrt(max(0.,1.-cthwz**2))
973  phiwz=vint(24)-phicm
974  p(ipu5,1)=pabs*sthwz*cos(phiwz)
975  p(ipu5,2)=pabs*sthwz*sin(phiwz)
976  p(ipu5,3)=pabs*cthwz
977  p(ipu5,4)=sqrt(pabs**2+p(ipu5,5)**2)
978  p(ipu6,1)=-p(ipu5,1)
979  p(ipu6,2)=-p(ipu5,2)
980  p(ipu6,3)=-p(ipu5,3)
981  p(ipu6,4)=sqrt(pabs**2+p(ipu6,5)**2)
982  CALL ludbrb(ipu5,ipu6,thecm,phicm,dble(bexcm),dble(beycm),
983  & dble(bezcm))
984  DO 430 jt=1,2
985  i1=mint(83)+8+jt
986  i2=mint(84)+4+jt
987  k(i1,1)=21
988  k(i1,2)=k(i2,2)
989  DO 430 j=1,5
990  430 p(i1,j)=p(i2,j)
991  n=ipu6
992  mint(7)=mint(83)+9
993  mint(8)=mint(83)+10
994  ENDIF
995 
996  IF(idoc.GE.8) THEN
997 C...Store colour connection indices.
998  DO 440 j=1,2
999  jc=j
1000  IF(kcs.EQ.-1) jc=3-j
1001  IF(icol(kcc,1,jc).NE.0.AND.k(ipu1,1).EQ.14) k(ipu1,j+3)=
1002  & k(ipu1,j+3)+mint(84)+icol(kcc,1,jc)
1003  IF(icol(kcc,2,jc).NE.0.AND.k(ipu2,1).EQ.14) k(ipu2,j+3)=
1004  & k(ipu2,j+3)+mint(84)+icol(kcc,2,jc)
1005  IF(icol(kcc,3,jc).NE.0.AND.k(ipu3,1).EQ.3) k(ipu3,j+3)=
1006  & mstu(5)*(mint(84)+icol(kcc,3,jc))
1007  440 IF(icol(kcc,4,jc).NE.0.AND.k(ipu4,1).EQ.3) k(ipu4,j+3)=
1008  & mstu(5)*(mint(84)+icol(kcc,4,jc))
1009 
1010 C...Copy outgoing partons to documentation lines.
1011  DO 450 i=1,2
1012  i1=mint(83)+idoc-2+i
1013  i2=mint(84)+2+i
1014  k(i1,1)=21
1015  k(i1,2)=k(i2,2)
1016  IF(idoc.LE.9) k(i1,3)=0
1017  IF(idoc.GE.11) k(i1,3)=mint(83)+2+i
1018  DO 450 j=1,5
1019  450 p(i1,j)=p(i2,j)
1020  ENDIF
1021  mint(52)=n
1022 
1023 C...Low-pT events: remove gluons used for string drawing purposes.
1024  IF(isub.EQ.95) THEN
1025  k(ipu3,1)=k(ipu3,1)+10
1026  k(ipu4,1)=k(ipu4,1)+10
1027  DO 460 j=41,66
1028  460 vint(j)=0.
1029  DO 470 i=mint(83)+5,mint(83)+8
1030  DO 470 j=1,5
1031  470 p(i,j)=0.
1032  ENDIF
1033 
1034  RETURN
1035  END