Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyshow.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyshow.f
1 
2 C*********************************************************************
3 
4 C...PYSHOW
5 C...Generates timelike parton showers from given partons.
6 
7  SUBROUTINE pyshow(IP1,IP2,QMAX)
8 
9 C...Double precision and integer declarations.
10  IMPLICIT DOUBLE PRECISION(a-h, o-z)
11  IMPLICIT INTEGER(i-n)
12  INTEGER pyk,pychge,pycomp
13 C...Parameter statement to help give large particle numbers.
14  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
15  &kexcit=4000000,kdimen=5000000)
16  parameter(maxnur=1000)
17 C...Commonblocks.
18  common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
19  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
20  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
21  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
22  common/pypars/mstp(200),parp(200),msti(200),pari(200)
23  common/pyint1/mint(400),vint(400)
24  SAVE /pypart/,/pyjets/,/pydat1/,/pydat2/,/pypars/,/pyint1/
25 C...Local arrays.
26  dimension pmth(5,140),ps(5),pma(100),pmsd(100),iep(100),ipa(100),
27  &kfla(100),kfld(100),kfl(100),itry(100),isi(100),isl(100),dp(100),
28  &dpt(5,4),ksh(0:140),kcii(2),niis(2),iiis(2,2),theiis(2,2),
29  &phiiis(2,2),isii(2),isset(2),iscol(0:140),ischg(0:140),
30  &iref(1000)
31 
32 C...Check that QMAX not too low.
33  IF(mstj(41).LE.0) THEN
34  RETURN
35  ELSEIF(mstj(41).EQ.1.OR.mstj(41).EQ.11) THEN
36  IF(qmax.LE.parj(82).AND.ip2.GE.-80) RETURN
37  ELSE
38  IF(qmax.LE.min(parj(82),parj(83),parj(90)).AND.ip2.GE.-80)
39  & RETURN
40  ENDIF
41 
42 C...Store positions of shower initiating partons.
43  mpspd=0
44  IF(ip1.GT.0.AND.ip1.LE.min(n,mstu(4)-mstu(32)).AND.ip2.EQ.0) THEN
45  npa=1
46  ipa(1)=ip1
47  ELSEIF(min(ip1,ip2).GT.0.AND.max(ip1,ip2).LE.min(n,mstu(4)-
48  & mstu(32))) THEN
49  npa=2
50  ipa(1)=ip1
51  ipa(2)=ip2
52  ELSEIF(ip1.GT.0.AND.ip1.LE.min(n,mstu(4)-mstu(32)).AND.ip2.LT.0
53  & .AND.ip2.GE.-80) THEN
54  npa=iabs(ip2)
55  DO 100 i=1,npa
56  ipa(i)=ip1+i-1
57  100 CONTINUE
58  ELSEIF(ip1.GT.0.AND.ip1.LE.min(n,mstu(4)-mstu(32)).AND.
59  &ip2.EQ.-100) THEN
60  mpspd=1
61  npa=2
62  ipa(1)=ip1+6
63  ipa(2)=ip1+7
64  ELSE
65  CALL pyerrm(12,
66  & '(PYSHOW:) failed to reconstruct showering system')
67  IF(mstu(21).GE.1) RETURN
68  ENDIF
69 
70 C...Send off to PYPTFS for pT-ordered evolution if requested,
71 C...if at least 2 partons, and without predefined shower branchings.
72  IF((mstj(41).EQ.11.OR.mstj(41).EQ.12).AND.npa.GE.2.AND.
73  &mpspd.EQ.0) THEN
74  npart=npa
75  DO 110 ii=1,npart
76  ipart(ii)=ipa(ii)
77  ptpart(ii)=0.5d0*qmax
78  110 CONTINUE
79  CALL pyptfs(2,0.5d0*qmax,0d0,ptgen)
80  RETURN
81  ENDIF
82 
83 C...Initialization of cutoff masses etc.
84  DO 120 ifl=0,40
85  iscol(ifl)=0
86  ischg(ifl)=0
87  ksh(ifl)=0
88  120 CONTINUE
89  iscol(21)=1
90  ksh(21)=1
91  pmth(1,21)=pymass(21)
92  pmth(2,21)=sqrt(pmth(1,21)**2+0.25d0*parj(82)**2)
93  pmth(3,21)=2d0*pmth(2,21)
94  pmth(4,21)=pmth(3,21)
95  pmth(5,21)=pmth(3,21)
96  pmth(1,22)=pymass(22)
97  pmth(2,22)=sqrt(pmth(1,22)**2+0.25d0*parj(83)**2)
98  pmth(3,22)=2d0*pmth(2,22)
99  pmth(4,22)=pmth(3,22)
100  pmth(5,22)=pmth(3,22)
101  pmqth1=parj(82)
102  IF(mstj(41).GE.2) pmqth1=min(parj(82),parj(83))
103  pmqt1e=min(pmqth1,parj(90))
104  pmqth2=pmth(2,21)
105  IF(mstj(41).GE.2) pmqth2=min(pmth(2,21),pmth(2,22))
106  pmqt2e=min(pmqth2,0.5d0*parj(90))
107  DO 130 ifl=1,5
108  iscol(ifl)=1
109  IF(mstj(41).GE.2) ischg(ifl)=1
110  ksh(ifl)=1
111  pmth(1,ifl)=pymass(ifl)
112  pmth(2,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*pmqth1**2)
113  pmth(3,ifl)=pmth(2,ifl)+pmqth2
114  pmth(4,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*parj(82)**2)+pmth(2,21)
115  pmth(5,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*parj(83)**2)+pmth(2,22)
116  130 CONTINUE
117  DO 140 ifl=11,15,2
118  IF(mstj(41).EQ.2.OR.mstj(41).GE.4) ischg(ifl)=1
119  IF(mstj(41).EQ.2.OR.mstj(41).GE.4) ksh(ifl)=1
120  pmth(1,ifl)=pymass(ifl)
121  pmth(2,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*parj(90)**2)
122  pmth(3,ifl)=pmth(2,ifl)+0.5d0*parj(90)
123  pmth(4,ifl)=pmth(3,ifl)
124  pmth(5,ifl)=pmth(3,ifl)
125  140 CONTINUE
126  pt2min=max(0.5d0*parj(82),1.1d0*parj(81))**2
127  alams=parj(81)**2
128  alfm=log(pt2min/alams)
129 
130 C...Check on phase space available for emission.
131  irej=0
132  DO 150 j=1,5
133  ps(j)=0d0
134  150 CONTINUE
135  pm=0d0
136  kfla(2)=0
137  DO 170 i=1,npa
138  kfla(i)=iabs(k(ipa(i),2))
139  pma(i)=p(ipa(i),5)
140 C...Special cutoff masses for initial partons (may be a heavy quark,
141 C...squark, ..., and need not be on the mass shell).
142  ir=30+i
143  IF(npa.LE.1) iref(i)=ir
144  IF(npa.GE.2) iref(i+1)=ir
145  iscol(ir)=0
146  ischg(ir)=0
147  ksh(ir)=0
148  IF(kfla(i).LE.8) THEN
149  iscol(ir)=1
150  IF(mstj(41).GE.2) ischg(ir)=1
151  ELSEIF(kfla(i).EQ.11.OR.kfla(i).EQ.13.OR.kfla(i).EQ.15.OR.
152  & kfla(i).EQ.17) THEN
153  IF(mstj(41).EQ.2.OR.mstj(41).GE.4) ischg(ir)=1
154  ELSEIF(kfla(i).EQ.21) THEN
155  iscol(ir)=1
156  ELSEIF((kfla(i).GE.ksusy1+1.AND.kfla(i).LE.ksusy1+8).OR.
157  & (kfla(i).GE.ksusy2+1.AND.kfla(i).LE.ksusy2+8)) THEN
158  iscol(ir)=1
159  ELSEIF(kfla(i).EQ.ksusy1+21) THEN
160  iscol(ir)=1
161 C...QUARKONIA+++
162 C...same for QQ~[3S18]
163  ELSEIF(mstp(148).GE.1.AND.(kfla(i).EQ.9900443.OR.
164  & kfla(i).EQ.9900553)) THEN
165  iscol(ir)=1
166 C...QUARKONIA---
167  ENDIF
168  IF(iscol(ir).EQ.1.OR.ischg(ir).EQ.1) ksh(ir)=1
169  pmth(1,ir)=pma(i)
170  IF(iscol(ir).EQ.1.AND.ischg(ir).EQ.1) THEN
171  pmth(2,ir)=sqrt(pmth(1,ir)**2+0.25d0*pmqth1**2)
172  pmth(3,ir)=pmth(2,ir)+pmqth2
173  pmth(4,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(82)**2)+pmth(2,21)
174  pmth(5,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(83)**2)+pmth(2,22)
175  ELSEIF(iscol(ir).EQ.1) THEN
176  pmth(2,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(82)**2)
177  pmth(3,ir)=pmth(2,ir)+0.5d0*parj(82)
178  pmth(4,ir)=pmth(3,ir)
179  pmth(5,ir)=pmth(3,ir)
180  ELSEIF(ischg(ir).EQ.1) THEN
181  pmth(2,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(90)**2)
182  pmth(3,ir)=pmth(2,ir)+0.5d0*parj(90)
183  pmth(4,ir)=pmth(3,ir)
184  pmth(5,ir)=pmth(3,ir)
185  ENDIF
186  IF(ksh(ir).EQ.1) pma(i)=pmth(3,ir)
187  pm=pm+pma(i)
188  IF(ksh(ir).EQ.0.OR.pma(i).GT.10d0*qmax) irej=irej+1
189  DO 160 j=1,4
190  ps(j)=ps(j)+p(ipa(i),j)
191  160 CONTINUE
192  170 CONTINUE
193  IF(irej.EQ.npa.AND.ip2.GE.-7) RETURN
194  ps(5)=sqrt(max(0d0,ps(4)**2-ps(1)**2-ps(2)**2-ps(3)**2))
195  IF(npa.EQ.1) ps(5)=ps(4)
196  IF(ps(5).LE.pm+pmqt1e) RETURN
197 
198 C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0).
199  kfsrce=0
200  IF(ip2.LE.0) THEN
201  ELSEIF(k(ip1,3).EQ.k(ip2,3).AND.k(ip1,3).GT.0) THEN
202  kfsrce=iabs(k(k(ip1,3),2))
203  ELSE
204  ipar1=max(1,k(ip1,3))
205  ipar2=max(1,k(ip2,3))
206  IF(k(ipar1,3).EQ.k(ipar2,3).AND.k(ipar1,3).GT.0)
207  & kfsrce=iabs(k(k(ipar1,3),2))
208  ENDIF
209  itypes=0
210  IF(kfsrce.GE.1.AND.kfsrce.LE.8) itypes=1
211  IF(kfsrce.GE.ksusy1+1.AND.kfsrce.LE.ksusy1+8) itypes=2
212  IF(kfsrce.GE.ksusy2+1.AND.kfsrce.LE.ksusy2+8) itypes=2
213  IF(kfsrce.GE.21.AND.kfsrce.LE.24) itypes=3
214  IF(kfsrce.GE.32.AND.kfsrce.LE.34) itypes=3
215  IF(kfsrce.EQ.25.OR.(kfsrce.GE.35.AND.kfsrce.LE.37)) itypes=4
216  IF(kfsrce.GE.ksusy1+22.AND.kfsrce.LE.ksusy1+37) itypes=5
217  IF(kfsrce.EQ.ksusy1+21) itypes=6
218 
219 C...Identify two primary showerers.
220  itype1=0
221  IF(kfla(1).GE.1.AND.kfla(1).LE.8) itype1=1
222  IF(kfla(1).GE.ksusy1+1.AND.kfla(1).LE.ksusy1+8) itype1=2
223  IF(kfla(1).GE.ksusy2+1.AND.kfla(1).LE.ksusy2+8) itype1=2
224  IF(kfla(1).GE.21.AND.kfla(1).LE.24) itype1=3
225  IF(kfla(1).GE.32.AND.kfla(1).LE.34) itype1=3
226  IF(kfla(1).EQ.25.OR.(kfla(1).GE.35.AND.kfla(1).LE.37)) itype1=4
227  IF(kfla(1).GE.ksusy1+22.AND.kfla(1).LE.ksusy1+37) itype1=5
228  IF(kfla(1).EQ.ksusy1+21) itype1=6
229  itype2=0
230  IF(kfla(2).GE.1.AND.kfla(2).LE.8) itype2=1
231  IF(kfla(2).GE.ksusy1+1.AND.kfla(2).LE.ksusy1+8) itype2=2
232  IF(kfla(2).GE.ksusy2+1.AND.kfla(2).LE.ksusy2+8) itype2=2
233  IF(kfla(2).GE.21.AND.kfla(2).LE.24) itype2=3
234  IF(kfla(2).GE.32.AND.kfla(2).LE.34) itype2=3
235  IF(kfla(2).EQ.25.OR.(kfla(2).GE.35.AND.kfla(2).LE.37)) itype2=4
236  IF(kfla(2).GE.ksusy1+22.AND.kfla(2).LE.ksusy1+37) itype2=5
237  IF(kfla(2).EQ.ksusy1+21) itype2=6
238 
239 C...Order of showerers. Presence of gluino.
240  itypmn=min(itype1,itype2)
241  itypmx=max(itype1,itype2)
242  iord=1
243  IF(itype1.GT.itype2) iord=2
244  iglui=0
245  IF(itype1.EQ.6.OR.itype2.EQ.6) iglui=1
246 
247 C...Check if 3-jet matrix elements to be used.
248  m3jc=0
249  alpha=0.5d0
250  IF(npa.EQ.2.AND.mstj(47).GE.1.AND.mpspd.EQ.0) THEN
251  IF(mstj(38).NE.0) THEN
252  m3jc=mstj(38)
253  alpha=parj(80)
254  mstj(38)=0
255  ELSEIF(mstj(47).GE.6) THEN
256  m3jc=mstj(47)
257  ELSE
258  iclass=1
259  icombi=4
260 
261 C...Vector/axial vector -> q + qbar; q -> q + V.
262  IF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.(itypes.EQ.0.OR.
263  & itypes.EQ.3)) THEN
264  iclass=2
265  IF(kfsrce.EQ.21.OR.kfsrce.EQ.22) THEN
266  icombi=1
267  ELSEIF(kfsrce.EQ.23.OR.(kfsrce.EQ.0.AND.
268  & k(ipa(1),2)+k(ipa(2),2).EQ.0)) THEN
269 C...gamma*/Z0: assume e+e- initial state if unknown.
270  ei=-1d0
271  IF(kfsrce.EQ.23) THEN
272  iannfl=k(k(ip1,3),3)
273  IF(iannfl.NE.0) THEN
274  kannfl=iabs(k(iannfl,2))
275  IF(kannfl.GE.1.AND.kannfl.LE.18) ei=kchg(kannfl,1)/3d0
276  ENDIF
277  ENDIF
278  ai=sign(1d0,ei+0.1d0)
279  vi=ai-4d0*ei*paru(102)
280  ef=kchg(kfla(1),1)/3d0
281  af=sign(1d0,ef+0.1d0)
282  vf=af-4d0*ef*paru(102)
283  xwc=1d0/(16d0*paru(102)*(1d0-paru(102)))
284  sh=ps(5)**2
285  sqmz=pmas(23,1)**2
286  sqwz=ps(5)*pmas(23,2)
287  sbwz=1d0/((sh-sqmz)**2+sqwz**2)
288  vect=ei**2*ef**2+2d0*ei*vi*ef*vf*xwc*sh*(sh-sqmz)*sbwz+
289  & (vi**2+ai**2)*vf**2*xwc**2*sh**2*sbwz
290  axiv=(vi**2+ai**2)*af**2*xwc**2*sh**2*sbwz
291  icombi=3
292  alpha=vect/(vect+axiv)
293  ELSEIF(kfsrce.EQ.24.OR.kfsrce.EQ.0) THEN
294  icombi=4
295  ENDIF
296 C...For chi -> chi q qbar, use V/A -> q qbar as first approximation.
297  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.itypes.EQ.5) THEN
298  iclass=2
299  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.3.AND.(itypes.EQ.0.OR.
300  & itypes.EQ.1)) THEN
301  iclass=3
302 
303 C...Scalar/pseudoscalar -> q + qbar; q -> q + S.
304  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.itypes.EQ.4) THEN
305  iclass=4
306  IF(kfsrce.EQ.25.OR.kfsrce.EQ.35.OR.kfsrce.EQ.37) THEN
307  icombi=1
308  ELSEIF(kfsrce.EQ.36) THEN
309  icombi=2
310  ENDIF
311  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.4.AND.(itypes.EQ.0.OR.
312  & itypes.EQ.1)) THEN
313  iclass=5
314 
315 C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S.
316  ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.2.AND.(itypes.EQ.0.OR.
317  & itypes.EQ.3)) THEN
318  iclass=6
319  ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.3.AND.(itypes.EQ.0.OR.
320  & itypes.EQ.2)) THEN
321  iclass=7
322  ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.2.AND.itypes.EQ.4) THEN
323  iclass=8
324  ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.4.AND.(itypes.EQ.0.OR.
325  & itypes.EQ.2)) THEN
326  iclass=9
327 
328 C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi.
329  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.2.AND.(itypes.EQ.0.OR.
330  & itypes.EQ.5)) THEN
331  iclass=10
332  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.5.AND.(itypes.EQ.0.OR.
333  & itypes.EQ.2)) THEN
334  iclass=11
335  ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.5.AND.(itypes.EQ.0.OR.
336  & itypes.EQ.1)) THEN
337  iclass=12
338 
339 C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g.
340  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.2.AND.itypes.EQ.6) THEN
341  iclass=13
342  ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.6.AND.(itypes.EQ.0.OR.
343  & itypes.EQ.2)) THEN
344  iclass=14
345  ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.6.AND.(itypes.EQ.0.OR.
346  & itypes.EQ.1)) THEN
347  iclass=15
348 
349 C...g -> ~g + ~g (eikonal approximation).
350  ELSEIF(itypmn.EQ.6.AND.itypmx.EQ.6.AND.itypes.EQ.0) THEN
351  iclass=16
352  ENDIF
353  m3jc=5*iclass+icombi
354  ENDIF
355  ENDIF
356 
357 C...Find if interference with initial state partons.
358  miis=0
359  IF(mstj(50).GE.1.AND.mstj(50).LE.3.AND.npa.EQ.2.AND.kfsrce.EQ.0
360  &.AND.mpspd.EQ.0) miis=mstj(50)
361  IF(mstj(50).GE.4.AND.mstj(50).LE.6.AND.npa.EQ.2.AND.mpspd.EQ.0)
362  &miis=mstj(50)-3
363  IF(miis.NE.0) THEN
364  DO 190 i=1,2
365  kcii(i)=0
366  kca=pycomp(kfla(i))
367  IF(kca.NE.0) kcii(i)=kchg(kca,2)*isign(1,k(ipa(i),2))
368  niis(i)=0
369  IF(kcii(i).NE.0) THEN
370  DO 180 j=1,2
371  icsi=mod(k(ipa(i),3+j)/mstu(5),mstu(5))
372  IF(icsi.GT.0.AND.icsi.NE.ipa(1).AND.icsi.NE.ipa(2).AND.
373  & (kcii(i).EQ.(-1)**(j+1).OR.kcii(i).EQ.2)) THEN
374  niis(i)=niis(i)+1
375  iiis(i,niis(i))=icsi
376  ENDIF
377  180 CONTINUE
378  ENDIF
379  190 CONTINUE
380  IF(niis(1)+niis(2).EQ.0) miis=0
381  ENDIF
382 
383 C...Boost interfering initial partons to rest frame
384 C...and reconstruct their polar and azimuthal angles.
385  IF(miis.NE.0) THEN
386  DO 210 i=1,2
387  DO 200 j=1,5
388  k(n+i,j)=k(ipa(i),j)
389  p(n+i,j)=p(ipa(i),j)
390  v(n+i,j)=0d0
391  200 CONTINUE
392  210 CONTINUE
393  DO 230 i=3,2+niis(1)
394  DO 220 j=1,5
395  k(n+i,j)=k(iiis(1,i-2),j)
396  p(n+i,j)=p(iiis(1,i-2),j)
397  v(n+i,j)=0d0
398  220 CONTINUE
399  230 CONTINUE
400  DO 250 i=3+niis(1),2+niis(1)+niis(2)
401  DO 240 j=1,5
402  k(n+i,j)=k(iiis(2,i-2-niis(1)),j)
403  p(n+i,j)=p(iiis(2,i-2-niis(1)),j)
404  v(n+i,j)=0d0
405  240 CONTINUE
406  250 CONTINUE
407  CALL pyrobo(n+1,n+2+niis(1)+niis(2),0d0,0d0,-ps(1)/ps(4),
408  & -ps(2)/ps(4),-ps(3)/ps(4))
409  phi=pyangl(p(n+1,1),p(n+1,2))
410  CALL pyrobo(n+1,n+2+niis(1)+niis(2),0d0,-phi,0d0,0d0,0d0)
411  the=pyangl(p(n+1,3),p(n+1,1))
412  CALL pyrobo(n+1,n+2+niis(1)+niis(2),-the,0d0,0d0,0d0,0d0)
413  DO 260 i=3,2+niis(1)
414  theiis(1,i-2)=pyangl(p(n+i,3),sqrt(p(n+i,1)**2+p(n+i,2)**2))
415  phiiis(1,i-2)=pyangl(p(n+i,1),p(n+i,2))
416  260 CONTINUE
417  DO 270 i=3+niis(1),2+niis(1)+niis(2)
418  theiis(2,i-2-niis(1))=paru(1)-pyangl(p(n+i,3),
419  & sqrt(p(n+i,1)**2+p(n+i,2)**2))
420  phiiis(2,i-2-niis(1))=pyangl(p(n+i,1),p(n+i,2))
421  270 CONTINUE
422  ENDIF
423 
424 C...Boost 3 or more partons to their rest frame.
425  IF(npa.GE.3) CALL pyrobo(ipa(1),ipa(npa),0d0,0d0,-ps(1)/ps(4),
426  &-ps(2)/ps(4),-ps(3)/ps(4))
427 
428 C...Define imagined single initiator of shower for parton system.
429  ns=n
430  IF(n.GT.mstu(4)-mstu(32)-10) THEN
431  CALL pyerrm(11,'(PYSHOW:) no more memory left in PYJETS')
432  IF(mstu(21).GE.1) RETURN
433  ENDIF
434  280 n=ns
435  IF(npa.GE.2) THEN
436  k(n+1,1)=11
437  k(n+1,2)=21
438  k(n+1,3)=0
439  k(n+1,4)=0
440  k(n+1,5)=0
441  p(n+1,1)=0d0
442  p(n+1,2)=0d0
443  p(n+1,3)=0d0
444  p(n+1,4)=ps(5)
445  p(n+1,5)=ps(5)
446  v(n+1,5)=ps(5)**2
447  n=n+1
448  iref(1)=21
449  ENDIF
450 
451 C...Loop over partons that may branch.
452  nep=npa
453  im=ns
454  IF(npa.EQ.1) im=ns-1
455  290 im=im+1
456  IF(n.GT.ns) THEN
457  IF(im.GT.n) goto 600
458  kflm=iabs(k(im,2))
459  ir=iref(im-ns)
460  IF(ksh(ir).EQ.0) goto 290
461  IF(p(im,5).LT.pmth(2,ir)) goto 290
462  igm=k(im,3)
463  ELSE
464  igm=-1
465  ENDIF
466  IF(n+nep.GT.mstu(4)-mstu(32)-10) THEN
467  CALL pyerrm(11,'(PYSHOW:) no more memory left in PYJETS')
468  IF(mstu(21).GE.1) RETURN
469  ENDIF
470 
471 C...Position of aunt (sister to branching parton).
472 C...Origin and flavour of daughters.
473  iau=0
474  IF(igm.GT.0) THEN
475  IF(k(im-1,3).EQ.igm) iau=im-1
476  IF(n.GE.im+1.AND.k(im+1,3).EQ.igm) iau=im+1
477  ENDIF
478  IF(igm.GE.0) THEN
479  k(im,4)=n+1
480  DO 300 i=1,nep
481  k(n+i,3)=im
482  300 CONTINUE
483  ELSE
484  k(n+1,3)=ipa(1)
485  ENDIF
486  IF(igm.LE.0) THEN
487  DO 310 i=1,nep
488  k(n+i,2)=k(ipa(i),2)
489  310 CONTINUE
490  ELSEIF(kflm.NE.21) THEN
491  k(n+1,2)=k(im,2)
492  k(n+2,2)=k(im,5)
493  iref(n+1-ns)=iref(im-ns)
494  iref(n+2-ns)=iabs(k(n+2,2))
495  ELSEIF(k(im,5).EQ.21) THEN
496  k(n+1,2)=21
497  k(n+2,2)=21
498  iref(n+1-ns)=21
499  iref(n+2-ns)=21
500  ELSE
501  k(n+1,2)=k(im,5)
502  k(n+2,2)=-k(im,5)
503  iref(n+1-ns)=iabs(k(n+1,2))
504  iref(n+2-ns)=iabs(k(n+2,2))
505  ENDIF
506 
507 C...Reset flags on daughters and tries made.
508  DO 320 ip=1,nep
509  k(n+ip,1)=3
510  k(n+ip,4)=0
511  k(n+ip,5)=0
512  kfld(ip)=iabs(k(n+ip,2))
513  IF(kchg(pycomp(kfld(ip)),2).EQ.0) k(n+ip,1)=1
514  itry(ip)=0
515  isl(ip)=0
516  isi(ip)=0
517  IF(ksh(iref(n+ip-ns)).EQ.1) isi(ip)=1
518  320 CONTINUE
519  islm=0
520 
521 C...Maximum virtuality of daughters.
522  IF(igm.LE.0) THEN
523  DO 330 i=1,npa
524  IF(npa.GE.3) p(n+i,4)=p(ipa(i),4)
525  p(n+i,5)=min(qmax,ps(5))
526  ir=iref(n+i-ns)
527  IF(ip2.LE.-8) p(n+i,5)=max(p(n+i,5),2d0*pmth(3,ir))
528  IF(isi(i).EQ.0) p(n+i,5)=p(ipa(i),5)
529  330 CONTINUE
530  ELSE
531  IF(mstj(43).LE.2) pem=v(im,2)
532  IF(mstj(43).GE.3) pem=p(im,4)
533  p(n+1,5)=min(p(im,5),v(im,1)*pem)
534  p(n+2,5)=min(p(im,5),(1d0-v(im,1))*pem)
535  IF(k(n+2,2).EQ.22) p(n+2,5)=pmth(1,22)
536  ENDIF
537  DO 340 i=1,nep
538  pmsd(i)=p(n+i,5)
539  IF(isi(i).EQ.1) THEN
540  ir=iref(n+i-ns)
541  IF(p(n+i,5).LE.pmth(3,ir)) p(n+i,5)=pmth(1,ir)
542  ENDIF
543  v(n+i,5)=p(n+i,5)**2
544  340 CONTINUE
545 
546 C...Choose one of the daughters for evolution.
547  350 inum=0
548  IF(nep.EQ.1) inum=1
549  DO 360 i=1,nep
550  IF(inum.EQ.0.AND.isl(i).EQ.1) inum=i
551  360 CONTINUE
552  DO 370 i=1,nep
553  IF(inum.EQ.0.AND.itry(i).EQ.0.AND.isi(i).EQ.1) THEN
554  ir=iref(n+i-ns)
555  IF(p(n+i,5).GE.pmth(2,ir)) inum=i
556  ENDIF
557  370 CONTINUE
558  IF(inum.EQ.0) THEN
559  rmax=0d0
560  DO 380 i=1,nep
561  IF(isi(i).EQ.1.AND.pmsd(i).GE.pmqt2e) THEN
562  rpm=p(n+i,5)/pmsd(i)
563  ir=iref(n+i-ns)
564  IF(rpm.GT.rmax.AND.p(n+i,5).GE.pmth(2,ir)) THEN
565  rmax=rpm
566  inum=i
567  ENDIF
568  ENDIF
569  380 CONTINUE
570  ENDIF
571 
572 C...Cancel choice of predetermined daughter already treated.
573  inum=max(1,inum)
574  inumt=inum
575  IF(mpspd.EQ.1.AND.igm.EQ.0.AND.itry(inumt).GE.1) THEN
576  IF(k(ip1-1+inum,4).GT.0) inum=3-inum
577  ELSEIF(mpspd.EQ.1.AND.im.EQ.ns+2.AND.itry(inumt).GE.1) THEN
578  IF(kfld(inumt).NE.21.AND.k(ip1+2,4).GT.0) inum=3-inum
579  IF(kfld(inumt).EQ.21.AND.k(ip1+3,4).GT.0) inum=3-inum
580  ENDIF
581 
582 C...Store information on choice of evolving daughter.
583  iep(1)=n+inum
584  DO 390 i=2,nep
585  iep(i)=iep(i-1)+1
586  IF(iep(i).GT.n+nep) iep(i)=n+1
587  390 CONTINUE
588  DO 400 i=1,nep
589  kfl(i)=iabs(k(iep(i),2))
590  400 CONTINUE
591  itry(inum)=itry(inum)+1
592  IF(itry(inum).GT.200) THEN
593  CALL pyerrm(14,'(PYSHOW:) caught in infinite loop')
594  IF(mstu(21).GE.1) RETURN
595  ENDIF
596  z=0.5d0
597  ir=iref(iep(1)-ns)
598  IF(ksh(ir).EQ.0) goto 450
599  IF(p(iep(1),5).LT.pmth(2,ir)) goto 450
600 
601 C...Check if evolution already predetermined for daughter.
602  ipspd=0
603  IF(mpspd.EQ.1.AND.igm.EQ.0) THEN
604  IF(k(ip1-1+inum,4).GT.0) ipspd=ip1-1+inum
605  ELSEIF(mpspd.EQ.1.AND.im.EQ.ns+2) THEN
606  IF(kfl(1).NE.21.AND.k(ip1+2,4).GT.0) ipspd=ip1+2
607  IF(kfl(1).EQ.21.AND.k(ip1+3,4).GT.0) ipspd=ip1+3
608  ENDIF
609  IF(inum.EQ.1.OR.inum.EQ.2) THEN
610  isset(inum)=0
611  IF(ipspd.NE.0) isset(inum)=1
612  ENDIF
613 
614 C...Select side for interference with initial state partons.
615  IF(miis.GE.1.AND.iep(1).LE.ns+3) THEN
616  iii=iep(1)-ns-1
617  isii(iii)=0
618  IF(iabs(kcii(iii)).EQ.1.AND.niis(iii).EQ.1) THEN
619  isii(iii)=1
620  ELSEIF(kcii(iii).EQ.2.AND.niis(iii).EQ.1) THEN
621  IF(pyr(0).GT.0.5d0) isii(iii)=1
622  ELSEIF(kcii(iii).EQ.2.AND.niis(iii).EQ.2) THEN
623  isii(iii)=1
624  IF(pyr(0).GT.0.5d0) isii(iii)=2
625  ENDIF
626  ENDIF
627 
628 C...Calculate allowed z range.
629  IF(nep.EQ.1) THEN
630  pmed=ps(4)
631  ELSEIF(igm.EQ.0.OR.mstj(43).LE.2) THEN
632  pmed=p(im,5)
633  ELSE
634  IF(inum.EQ.1) pmed=v(im,1)*pem
635  IF(inum.EQ.2) pmed=(1d0-v(im,1))*pem
636  ENDIF
637  IF(mod(mstj(43),2).EQ.1) THEN
638  zc=pmth(2,21)/pmed
639  zce=pmth(2,22)/pmed
640  IF(iscol(ir).EQ.0) zce=0.5d0*parj(90)/pmed
641  ELSE
642  zc=0.5d0*(1d0-sqrt(max(0d0,1d0-(2d0*pmth(2,21)/pmed)**2)))
643  IF(zc.LT.1d-6) zc=(pmth(2,21)/pmed)**2
644  pmtmpe=pmth(2,22)
645  IF(iscol(ir).EQ.0) pmtmpe=0.5d0*parj(90)
646  zce=0.5d0*(1d0-sqrt(max(0d0,1d0-(2d0*pmtmpe/pmed)**2)))
647  IF(zce.LT.1d-6) zce=(pmtmpe/pmed)**2
648  ENDIF
649  zc=min(zc,0.491d0)
650  zce=min(zce,0.49991d0)
651  IF(((mstj(41).EQ.1.AND.zc.GT.0.49d0).OR.(mstj(41).GE.2.AND.
652  &min(zc,zce).GT.0.4999d0)).AND.ipspd.EQ.0) THEN
653  p(iep(1),5)=pmth(1,ir)
654  v(iep(1),5)=p(iep(1),5)**2
655  goto 450
656  ENDIF
657 
658 C...Integral of Altarelli-Parisi z kernel for QCD.
659 C...(Includes squark and gluino; with factor N_C/C_F extra for latter).
660  IF(mstj(49).EQ.0.AND.kfl(1).EQ.21) THEN
661  fbr=6d0*log((1d0-zc)/zc)+mstj(45)*0.5d0
662 C...QUARKONIA+++
663 C...Evolution of QQ~[3S18] state if MSTP(148)=1.
664  ELSEIF(mstj(49).EQ.0.AND.mstp(149).GE.0.AND.
665  & (kfl(1).EQ.9900443.OR.kfl(1).EQ.9900553)) THEN
666  fbr=6d0*log((1d0-zc)/zc)
667 C...QUARKONIA---
668  ELSEIF(mstj(49).EQ.0) THEN
669  fbr=(8d0/3d0)*log((1d0-zc)/zc)
670  IF(iglui.EQ.1.AND.ir.GE.31) fbr=fbr*(9d0/4d0)
671 
672 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
673  ELSEIF(mstj(49).EQ.1.AND.kfl(1).EQ.21) THEN
674  fbr=(parj(87)+mstj(45)*parj(88))*(1d0-2d0*zc)
675  ELSEIF(mstj(49).EQ.1) THEN
676  fbr=(1d0-2d0*zc)/3d0
677  IF(igm.EQ.0.AND.m3jc.GE.1) fbr=4d0*fbr
678 
679 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
680  ELSEIF(kfl(1).EQ.21) THEN
681  fbr=6d0*mstj(45)*(0.5d0-zc)
682  ELSE
683  fbr=2d0*log((1d0-zc)/zc)
684  ENDIF
685 
686 C...Reset QCD probability for colourless.
687  IF(iscol(ir).EQ.0) fbr=0d0
688 
689 C...Integral of Altarelli-Parisi kernel for photon emission.
690  fbre=0d0
691  IF(mstj(41).GE.2.AND.ischg(ir).EQ.1) THEN
692  IF(kfl(1).LE.18) THEN
693  fbre=(kchg(kfl(1),1)/3d0)**2*2d0*log((1d0-zce)/zce)
694  ENDIF
695  IF(mstj(41).EQ.10) fbre=parj(84)*fbre
696  ENDIF
697 
698 C...Inner veto algorithm starts. Find maximum mass for evolution.
699  410 pms=v(iep(1),5)
700  IF(igm.GE.0) THEN
701  pm2=0d0
702  DO 420 i=2,nep
703  pm=p(iep(i),5)
704  iri=iref(iep(i)-ns)
705  IF(ksh(iri).EQ.1) pm=pmth(2,iri)
706  pm2=pm2+pm
707  420 CONTINUE
708  pms=min(pms,(p(im,5)-pm2)**2)
709  ENDIF
710 
711 C...Select mass for daughter in QCD evolution.
712  b0=27d0/6d0
713  DO 430 iff=4,mstj(45)
714  IF(pms.GT.4d0*pmth(2,iff)**2) b0=(33d0-2d0*iff)/6d0
715  430 CONTINUE
716 C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
717  pmsc=max(0.5d0*parj(82),pms-pmth(1,ir)**2)
718 C...Already predetermined choice.
719  IF(ipspd.NE.0) THEN
720  pmsqcd=p(ipspd,5)**2
721  ELSEIF(fbr.LT.1d-3) THEN
722  pmsqcd=0d0
723  ELSEIF(mstj(44).LE.0) THEN
724  pmsqcd=pmsc*exp(max(-50d0,log(pyr(0))*paru(2)/(paru(111)*fbr)))
725  ELSEIF(mstj(44).EQ.1) THEN
726  pmsqcd=4d0*alams*(0.25d0*pmsc/alams)**(pyr(0)**(b0/fbr))
727  ELSE
728  pmsqcd=pmsc*exp(max(-50d0,alfm*b0*log(pyr(0))/fbr))
729  ENDIF
730 C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
731  IF(ipspd.EQ.0) pmsqcd=pmsqcd+pmth(1,ir)**2
732  IF(zc.GT.0.49d0.OR.pmsqcd.LE.pmth(4,ir)**2) pmsqcd=pmth(2,ir)**2
733  v(iep(1),5)=pmsqcd
734  mce=1
735 
736 C...Select mass for daughter in QED evolution.
737  IF(mstj(41).GE.2.AND.ischg(ir).EQ.1.AND.ipspd.EQ.0) THEN
738 C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
739  pmse=max(0.5d0*parj(83),pms-pmth(1,ir)**2)
740  IF(fbre.LT.1d-3) THEN
741  pmsqed=0d0
742  ELSE
743  pmsqed=pmse*exp(max(-50d0,log(pyr(0))*paru(2)/
744  & (paru(101)*fbre)))
745  ENDIF
746 C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
747  pmsqed=pmsqed+pmth(1,ir)**2
748  IF(zce.GT.0.4999d0.OR.pmsqed.LE.pmth(5,ir)**2) pmsqed=
749  & pmth(2,ir)**2
750  IF(pmsqed.GT.pmsqcd) THEN
751  v(iep(1),5)=pmsqed
752  mce=2
753  ENDIF
754  ENDIF
755 
756 C...Check whether daughter mass below cutoff.
757  p(iep(1),5)=sqrt(v(iep(1),5))
758  IF(p(iep(1),5).LE.pmth(3,ir)) THEN
759  p(iep(1),5)=pmth(1,ir)
760  v(iep(1),5)=p(iep(1),5)**2
761  goto 450
762  ENDIF
763 
764 C...Already predetermined choice of z, and flavour in g -> qqbar.
765  IF(ipspd.NE.0) THEN
766  ipsgd1=k(ipspd,4)
767  ipsgd2=k(ipspd,5)
768  pmsgd1=p(ipsgd1,5)**2
769  pmsgd2=p(ipsgd2,5)**2
770  alamps=sqrt(max(1d-10,(pmsqcd-pmsgd1-pmsgd2)**2-
771  & 4d0*pmsgd1*pmsgd2))
772  z=0.5d0*(pmsqcd*(2d0*p(ipsgd1,4)/p(ipspd,4)-1d0)+alamps-
773  & pmsgd1+pmsgd2)/alamps
774  z=max(0.00001d0,min(0.99999d0,z))
775  IF(kfl(1).NE.21) THEN
776  k(iep(1),5)=21
777  ELSE
778  k(iep(1),5)=iabs(k(ipsgd1,2))
779  ENDIF
780 
781 C...Select z value of branching: q -> qgamma.
782  ELSEIF(mce.EQ.2) THEN
783  z=1d0-(1d0-zce)*(zce/(1d0-zce))**pyr(0)
784  IF(1d0+z**2.LT.2d0*pyr(0)) goto 410
785  k(iep(1),5)=22
786 
787 C...QUARKONIA+++
788 C...Select z value of branching: QQ~[3S18] -> QQ~[3S18]g.
789  ELSEIF(mstj(49).EQ.0.AND.
790  & (kfl(1).EQ.9900443.OR.kfl(1).EQ.9900553)) THEN
791  z=(1d0-zc)*(zc/(1d0-zc))**pyr(0)
792 C...Select always the harder 'gluon' if the switch MSTP(149)<=0.
793  IF(mstp(149).LE.0.OR.pyr(0).GT.0.5d0) z=1d0-z
794  IF((1d0-z*(1d0-z))**2.LT.pyr(0)) goto 410
795  k(iep(1),5)=21
796 C...QUARKONIA---
797 
798 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
799  ELSEIF(mstj(49).NE.1.AND.kfl(1).NE.21) THEN
800  z=1d0-(1d0-zc)*(zc/(1d0-zc))**pyr(0)
801 C...Only do z weighting when no ME correction afterwards.
802  IF(m3jc.EQ.0.AND.1d0+z**2.LT.2d0*pyr(0)) goto 410
803  k(iep(1),5)=21
804  ELSEIF(mstj(49).EQ.0.AND.mstj(45)*0.5d0.LT.pyr(0)*fbr) THEN
805  z=(1d0-zc)*(zc/(1d0-zc))**pyr(0)
806  IF(pyr(0).GT.0.5d0) z=1d0-z
807  IF((1d0-z*(1d0-z))**2.LT.pyr(0)) goto 410
808  k(iep(1),5)=21
809  ELSEIF(mstj(49).NE.1) THEN
810  z=pyr(0)
811  IF(z**2+(1d0-z)**2.LT.pyr(0)) goto 410
812  kflb=1+int(mstj(45)*pyr(0))
813  pmq=4d0*pmth(2,kflb)**2/v(iep(1),5)
814  IF(pmq.GE.1d0) goto 410
815  IF(mstj(44).LE.2.OR.mstj(44).EQ.4) THEN
816  IF(z.LT.zc.OR.z.GT.1d0-zc) goto 410
817  pmq0=4d0*pmth(2,21)**2/v(iep(1),5)
818  IF(mod(mstj(43),2).EQ.0.AND.(1d0+0.5d0*pmq)*sqrt(1d0-pmq)
819  & .LT.pyr(0)*(1d0+0.5d0*pmq0)*sqrt(1d0-pmq0)) goto 410
820  ELSE
821  IF((1d0+0.5d0*pmq)*sqrt(1d0-pmq).LT.pyr(0)) goto 410
822  ENDIF
823  k(iep(1),5)=kflb
824 
825 C...Ditto for scalar gluon model.
826  ELSEIF(kfl(1).NE.21) THEN
827  z=1d0-sqrt(zc**2+pyr(0)*(1d0-2d0*zc))
828  k(iep(1),5)=21
829  ELSEIF(pyr(0)*(parj(87)+mstj(45)*parj(88)).LE.parj(87)) THEN
830  z=zc+(1d0-2d0*zc)*pyr(0)
831  k(iep(1),5)=21
832  ELSE
833  z=zc+(1d0-2d0*zc)*pyr(0)
834  kflb=1+int(mstj(45)*pyr(0))
835  pmq=4d0*pmth(2,kflb)**2/v(iep(1),5)
836  IF(pmq.GE.1d0) goto 410
837  k(iep(1),5)=kflb
838  ENDIF
839 
840 C...Correct to alpha_s(pT^2) (optionally m^2/4 for g -> q qbar).
841  IF(mce.EQ.1.AND.mstj(44).GE.2.AND.ipspd.EQ.0) THEN
842  IF(kfl(1).EQ.21.AND.k(iep(1),5).LT.10.AND.
843  & (mstj(44).EQ.3.OR.mstj(44).EQ.5)) THEN
844  IF(alfm/log(v(iep(1),5)*0.25d0/alams).LT.pyr(0)) goto 410
845  ELSE
846  pt2app=z*(1d0-z)*v(iep(1),5)
847  IF(mstj(44).GE.4) pt2app=pt2app*
848  & (1d0-pmth(1,ir)**2/v(iep(1),5))**2
849  IF(pt2app.LT.pt2min) goto 410
850  IF(alfm/log(pt2app/alams).LT.pyr(0)) goto 410
851  ENDIF
852  ENDIF
853 
854 C...Check if z consistent with chosen m.
855  IF(kfl(1).EQ.21) THEN
856  irgd1=iabs(k(iep(1),5))
857  irgd2=irgd1
858  ELSE
859  irgd1=ir
860  irgd2=iabs(k(iep(1),5))
861  ENDIF
862  IF(nep.EQ.1) THEN
863  ped=ps(4)
864  ELSEIF(nep.GE.3) THEN
865  ped=p(iep(1),4)
866  ELSEIF(igm.EQ.0.OR.mstj(43).LE.2) THEN
867  ped=0.5d0*(v(im,5)+v(iep(1),5)-pm2**2)/p(im,5)
868  ELSE
869  IF(iep(1).EQ.n+1) ped=v(im,1)*pem
870  IF(iep(1).EQ.n+2) ped=(1d0-v(im,1))*pem
871  ENDIF
872  IF(mod(mstj(43),2).EQ.1) THEN
873  pmqth3=0.5d0*parj(82)
874  IF(irgd2.EQ.22) pmqth3=0.5d0*parj(83)
875  IF(irgd2.EQ.22.AND.iscol(ir).EQ.0) pmqth3=0.5d0*parj(90)
876  pmq1=(pmth(1,irgd1)**2+pmqth3**2)/v(iep(1),5)
877  pmq2=(pmth(1,irgd2)**2+pmqth3**2)/v(iep(1),5)
878  zd=sqrt(max(0d0,(1d0-v(iep(1),5)/ped**2)*((1d0-pmq1-pmq2)**2-
879  & 4d0*pmq1*pmq2)))
880  zh=1d0+pmq1-pmq2
881  ELSE
882  zd=sqrt(max(0d0,1d0-v(iep(1),5)/ped**2))
883  zh=1d0
884  ENDIF
885  IF(kfl(1).EQ.21.AND.k(iep(1),5).LT.10.AND.
886  &(mstj(44).EQ.3.OR.mstj(44).EQ.5)) THEN
887  ELSEIF(ipspd.NE.0) THEN
888  ELSE
889  zl=0.5d0*(zh-zd)
890  zu=0.5d0*(zh+zd)
891  IF(z.LT.zl.OR.z.GT.zu) goto 410
892  ENDIF
893  IF(kfl(1).EQ.21) v(iep(1),3)=log(zu*(1d0-zl)/max(1d-20,zl*
894  &(1d0-zu)))
895  IF(kfl(1).NE.21) v(iep(1),3)=log((1d0-zl)/max(1d-10,1d0-zu))
896 
897 C...Width suppression for q -> q + g.
898  IF(mstj(40).NE.0.AND.kfl(1).NE.21.AND.ipspd.EQ.0) THEN
899  IF(igm.EQ.0) THEN
900  eglu=0.5d0*ps(5)*(1d0-z)*(1d0+v(iep(1),5)/v(ns+1,5))
901  ELSE
902  eglu=pmed*(1d0-z)
903  ENDIF
904  chi=parj(89)**2/(parj(89)**2+eglu**2)
905  IF(mstj(40).EQ.1) THEN
906  IF(chi.LT.pyr(0)) goto 410
907  ELSEIF(mstj(40).EQ.2) THEN
908  IF(1d0-chi.LT.pyr(0)) goto 410
909  ENDIF
910  ENDIF
911 
912 C...Three-jet matrix element correction.
913  IF(m3jc.GE.1) THEN
914  wme=1d0
915  wshow=1d0
916 
917 C...QED matrix elements: only for massless case so far.
918  IF(mce.EQ.2.AND.igm.EQ.0) THEN
919  x1=z*(1d0+v(iep(1),5)/v(ns+1,5))
920  x2=1d0-v(iep(1),5)/v(ns+1,5)
921  x3=(1d0-x1)+(1d0-x2)
922  ki1=k(ipa(inum),2)
923  ki2=k(ipa(3-inum),2)
924  qf1=kchg(pycomp(ki1),1)*isign(1,ki1)/3d0
925  qf2=kchg(pycomp(ki2),1)*isign(1,ki2)/3d0
926  wshow=qf1**2*(1d0-x1)/x3*(1d0+(x1/(2d0-x2))**2)+
927  & qf2**2*(1d0-x2)/x3*(1d0+(x2/(2d0-x1))**2)
928  wme=(qf1*(1d0-x1)/x3-qf2*(1d0-x2)/x3)**2*(x1**2+x2**2)
929  ELSEIF(mce.EQ.2) THEN
930 
931 C...QCD matrix elements, including mass effects.
932  ELSEIF(mstj(49).NE.1.AND.k(iep(1),2).NE.21) THEN
933  ps1me=v(iep(1),5)
934  pm1me=pmth(1,ir)
935  m3jcc=m3jc
936  IF(ir.GE.31.AND.igm.EQ.0) THEN
937 C...QCD ME: original parton, first branching.
938  pm2me=pmth(1,63-ir)
939  ecmme=ps(5)
940  ELSEIF(ir.GE.31) THEN
941 C...QCD ME: original parton, subsequent branchings.
942  pm2me=pmth(1,63-ir)
943  pedme=pem*(v(im,1)+(1d0-v(im,1))*ps1me/v(im,5))
944  ecmme=pedme+sqrt(max(0d0,pedme**2-ps1me+pm2me**2))
945  ELSEIF(k(im,2).EQ.21) THEN
946 C...QCD ME: secondary partons, first branching.
947  pm2me=pm1me
948  zmme=v(im,1)
949  IF(iep(1).GT.iep(2)) zmme=1d0-zmme
950  pmlme=sqrt(max(0d0,(v(im,5)-ps1me-pm2me**2)**2-
951  & 4d0*ps1me*pm2me**2))
952  pedme=pem*(0.5d0*(v(im,5)-pmlme+ps1me-pm2me**2)+pmlme*zmme)/
953  & v(im,5)
954  ecmme=pedme+sqrt(max(0d0,pedme**2-ps1me+pm2me**2))
955  m3jcc=66
956  ELSE
957 C...QCD ME: secondary partons, subsequent branchings.
958  pm2me=pm1me
959  pedme=pem*(v(im,1)+(1d0-v(im,1))*ps1me/v(im,5))
960  ecmme=pedme+sqrt(max(0d0,pedme**2-ps1me+pm2me**2))
961  m3jcc=66
962  ENDIF
963 C...Construct ME variables.
964  r1me=pm1me/ecmme
965  r2me=pm2me/ecmme
966  x1=(1d0+ps1me/ecmme**2-r2me**2)*(z+(1d0-z)*pm1me**2/ps1me)
967  x2=1d0+r2me**2-ps1me/ecmme**2
968 C...Call ME, with right order important for two inequivalent showerers.
969  IF(ir.EQ.iord+30) THEN
970  wme=pymael(m3jcc,x1,x2,r1me,r2me,alpha)
971  ELSE
972  wme=pymael(m3jcc,x2,x1,r2me,r1me,alpha)
973  ENDIF
974 C...Split up total ME when two radiating partons.
975  isprad=1
976  IF((m3jcc.GE.16.AND.m3jcc.LE.19).OR.
977  & (m3jcc.GE.26.AND.m3jcc.LE.29).OR.
978  & (m3jcc.GE.36.AND.m3jcc.LE.39).OR.
979  & (m3jcc.GE.46.AND.m3jcc.LE.49).OR.
980  & (m3jcc.GE.56.AND.m3jcc.LE.64)) isprad=0
981  IF(isprad.EQ.1) wme=wme*max(1d-10,1d0+r1me**2-r2me**2-x1)/
982  & max(1d-10,2d0-x1-x2)
983 C...Evaluate shower rate to be compared with.
984  wshow=2d0/(max(1d-10,2d0-x1-x2)*
985  & max(1d-10,1d0+r2me**2-r1me**2-x2))
986  IF(iglui.EQ.1.AND.ir.GE.31) wshow=(9d0/4d0)*wshow
987  ELSEIF(mstj(49).NE.1) THEN
988 
989 C...Toy model scalar theory matrix elements; no mass effects.
990  ELSE
991  x1=z*(1d0+v(iep(1),5)/v(ns+1,5))
992  x2=1d0-v(iep(1),5)/v(ns+1,5)
993  x3=(1d0-x1)+(1d0-x2)
994  wshow=4d0*x3*((1d0-x1)/(2d0-x2)**2+(1d0-x2)/(2d0-x1)**2)
995  wme=x3**2
996  IF(mstj(102).GE.2) wme=x3**2-2d0*(1d0+x3)*(1d0-x1)*(1d0-x2)*
997  & parj(171)
998  ENDIF
999 
1000  IF(wme.LT.pyr(0)*wshow) goto 410
1001  ENDIF
1002 
1003 C...Impose angular ordering by rejection of nonordered emission.
1004  IF(mce.EQ.1.AND.igm.GT.0.AND.mstj(42).GE.2.AND.ipspd.EQ.0) THEN
1005  pemao=v(im,1)*p(im,4)
1006  IF(iep(1).EQ.n+2) pemao=(1d0-v(im,1))*p(im,4)
1007  IF(ir.GE.31.AND.mstj(42).GE.5) THEN
1008  maod=0
1009  ELSEIF(kfl(1).EQ.21.AND.k(iep(1),5).LE.10.AND.(mstj(42).EQ.4
1010  & .OR.mstj(42).EQ.7)) THEN
1011  maod=0
1012  ELSEIF(kfl(1).EQ.21.AND.k(iep(1),5).LE.10.AND.(mstj(42).EQ.3
1013  & .OR.mstj(42).EQ.6)) THEN
1014  maod=1
1015  pmdao=pmth(2,k(iep(1),5))
1016  the2id=z*(1d0-z)*pemao**2/(v(iep(1),5)-4d0*pmdao**2)
1017  ELSE
1018  maod=1
1019  the2id=z*(1d0-z)*pemao**2/v(iep(1),5)
1020  IF(mstj(42).GE.3.AND.mstj(42).NE.5) the2id=the2id*
1021  & (1d0+pmth(1,ir)**2*(1d0-z)/(v(iep(1),5)*z))**2
1022  ENDIF
1023  maom=1
1024  iaom=im
1025  440 IF(k(iaom,5).EQ.22) THEN
1026  iaom=k(iaom,3)
1027  IF(k(iaom,3).LE.ns) maom=0
1028  IF(maom.EQ.1) goto 440
1029  ENDIF
1030  IF(maom.EQ.1.AND.maod.EQ.1) THEN
1031  the2im=v(iaom,1)*(1d0-v(iaom,1))*p(iaom,4)**2/v(iaom,5)
1032  IF(the2id.LT.the2im) goto 410
1033  ENDIF
1034  ENDIF
1035 
1036 C...Impose user-defined maximum angle at first branching.
1037  IF(mstj(48).EQ.1.AND.ipspd.EQ.0) THEN
1038  IF(nep.EQ.1.AND.im.EQ.ns) THEN
1039  the2id=z*(1d0-z)*ps(4)**2/v(iep(1),5)
1040  IF(parj(85)**2*the2id.LT.1d0) goto 410
1041  ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+2) THEN
1042  the2id=z*(1d0-z)*(0.5d0*p(im,4))**2/v(iep(1),5)
1043  IF(parj(85)**2*the2id.LT.1d0) goto 410
1044  ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+3) THEN
1045  the2id=z*(1d0-z)*(0.5d0*p(im,4))**2/v(iep(1),5)
1046  IF(parj(86)**2*the2id.LT.1d0) goto 410
1047  ENDIF
1048  ENDIF
1049 
1050 C...Impose angular constraint in first branching from interference
1051 C...with initial state partons.
1052  IF(miis.GE.2.AND.iep(1).LE.ns+3) THEN
1053  the2d=max((1d0-z)/z,z/(1d0-z))*v(iep(1),5)/(0.5d0*p(im,4))**2
1054  IF(iep(1).EQ.ns+2.AND.isii(1).GE.1) THEN
1055  IF(the2d.GT.theiis(1,isii(1))**2) goto 410
1056  ELSEIF(iep(1).EQ.ns+3.AND.isii(2).GE.1) THEN
1057  IF(the2d.GT.theiis(2,isii(2))**2) goto 410
1058  ENDIF
1059  ENDIF
1060 
1061 C...End of inner veto algorithm. Check if only one leg evolved so far.
1062  450 v(iep(1),1)=z
1063  isl(1)=0
1064  isl(2)=0
1065  IF(nep.EQ.1) goto 490
1066  IF(nep.EQ.2.AND.p(iep(1),5)+p(iep(2),5).GE.p(im,5)) goto 350
1067  DO 460 i=1,nep
1068  ir=iref(n+i-ns)
1069  IF(itry(i).EQ.0.AND.ksh(ir).EQ.1) THEN
1070  IF(p(n+i,5).GE.pmth(2,ir)) goto 350
1071  ENDIF
1072  460 CONTINUE
1073 
1074 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
1075  IF(nep.GE.3) THEN
1076  pmsum=0d0
1077  DO 470 i=1,nep
1078  pmsum=pmsum+p(n+i,5)
1079  470 CONTINUE
1080  IF(pmsum.GE.ps(5)) goto 350
1081  ELSEIF(igm.EQ.0.OR.mstj(43).LE.2.OR.mod(mstj(43),2).EQ.0) THEN
1082  DO 480 i1=n+1,n+2
1083  irda=iref(i1-ns)
1084  IF(ksh(irda).EQ.0) goto 480
1085  IF(p(i1,5).LT.pmth(2,irda)) goto 480
1086  IF(irda.EQ.21) THEN
1087  irgd1=iabs(k(i1,5))
1088  irgd2=irgd1
1089  ELSE
1090  irgd1=irda
1091  irgd2=iabs(k(i1,5))
1092  ENDIF
1093  i2=2*n+3-i1
1094  IF(igm.EQ.0.OR.mstj(43).LE.2) THEN
1095  ped=0.5d0*(v(im,5)+v(i1,5)-v(i2,5))/p(im,5)
1096  ELSE
1097  IF(i1.EQ.n+1) zm=v(im,1)
1098  IF(i1.EQ.n+2) zm=1d0-v(im,1)
1099  pml=sqrt((v(im,5)-v(n+1,5)-v(n+2,5))**2-
1100  & 4d0*v(n+1,5)*v(n+2,5))
1101  ped=pem*(0.5d0*(v(im,5)-pml+v(i1,5)-v(i2,5))+pml*zm)/
1102  & v(im,5)
1103  ENDIF
1104  IF(mod(mstj(43),2).EQ.1) THEN
1105  pmqth3=0.5d0*parj(82)
1106  IF(irgd2.EQ.22) pmqth3=0.5d0*parj(83)
1107  IF(irgd2.EQ.22.AND.iscol(irda).EQ.0) pmqth3=0.5d0*parj(90)
1108  pmq1=(pmth(1,irgd1)**2+pmqth3**2)/v(i1,5)
1109  pmq2=(pmth(1,irgd2)**2+pmqth3**2)/v(i1,5)
1110  zd=sqrt(max(0d0,(1d0-v(i1,5)/ped**2)*((1d0-pmq1-pmq2)**2-
1111  & 4d0*pmq1*pmq2)))
1112  zh=1d0+pmq1-pmq2
1113  ELSE
1114  zd=sqrt(max(0d0,1d0-v(i1,5)/ped**2))
1115  zh=1d0
1116  ENDIF
1117  IF(irda.EQ.21.AND.irgd1.LT.10.AND.
1118  & (mstj(44).EQ.3.OR.mstj(44).EQ.5)) THEN
1119  ELSE
1120  zl=0.5d0*(zh-zd)
1121  zu=0.5d0*(zh+zd)
1122  IF(i1.EQ.n+1.AND.(v(i1,1).LT.zl.OR.v(i1,1).GT.zu).AND.
1123  & isset(1).EQ.0) THEN
1124  isl(1)=1
1125  ELSEIF(i1.EQ.n+2.AND.(v(i1,1).LT.zl.OR.v(i1,1).GT.zu).AND.
1126  & isset(2).EQ.0) THEN
1127  isl(2)=1
1128  ENDIF
1129  ENDIF
1130  IF(irda.EQ.21) v(i1,4)=log(zu*(1d0-zl)/max(1d-20,
1131  & zl*(1d0-zu)))
1132  IF(irda.NE.21) v(i1,4)=log((1d0-zl)/max(1d-10,1d0-zu))
1133  480 CONTINUE
1134  IF(isl(1).EQ.1.AND.isl(2).EQ.1.AND.islm.NE.0) THEN
1135  isl(3-islm)=0
1136  islm=3-islm
1137  ELSEIF(isl(1).EQ.1.AND.isl(2).EQ.1) THEN
1138  zdr1=max(0d0,v(n+1,3)/max(1d-6,v(n+1,4))-1d0)
1139  zdr2=max(0d0,v(n+2,3)/max(1d-6,v(n+2,4))-1d0)
1140  IF(zdr2.GT.pyr(0)*(zdr1+zdr2)) isl(1)=0
1141  IF(isl(1).EQ.1) isl(2)=0
1142  IF(isl(1).EQ.0) islm=1
1143  IF(isl(2).EQ.0) islm=2
1144  ENDIF
1145  IF(isl(1).EQ.1.OR.isl(2).EQ.1) goto 350
1146  ENDIF
1147  ird1=iref(n+1-ns)
1148  ird2=iref(n+2-ns)
1149  IF(igm.GT.0) THEN
1150  IF(mod(mstj(43),2).EQ.1.AND.(p(n+1,5).GE.
1151  & pmth(2,ird1).OR.p(n+2,5).GE.pmth(2,ird2))) THEN
1152  pmq1=v(n+1,5)/v(im,5)
1153  pmq2=v(n+2,5)/v(im,5)
1154  zd=sqrt(max(0d0,(1d0-v(im,5)/pem**2)*((1d0-pmq1-pmq2)**2-
1155  & 4d0*pmq1*pmq2)))
1156  zh=1d0+pmq1-pmq2
1157  zl=0.5d0*(zh-zd)
1158  zu=0.5d0*(zh+zd)
1159  IF(v(im,1).LT.zl.OR.v(im,1).GT.zu) goto 350
1160  ENDIF
1161  ENDIF
1162 
1163 C...Accepted branch. Construct four-momentum for initial partons.
1164  490 mazip=0
1165  mazic=0
1166  IF(nep.EQ.1) THEN
1167  p(n+1,1)=0d0
1168  p(n+1,2)=0d0
1169  p(n+1,3)=sqrt(max(0d0,(p(ipa(1),4)+p(n+1,5))*(p(ipa(1),4)-
1170  & p(n+1,5))))
1171  p(n+1,4)=p(ipa(1),4)
1172  v(n+1,2)=p(n+1,4)
1173  ELSEIF(igm.EQ.0.AND.nep.EQ.2) THEN
1174  ped1=0.5d0*(v(im,5)+v(n+1,5)-v(n+2,5))/p(im,5)
1175  p(n+1,1)=0d0
1176  p(n+1,2)=0d0
1177  p(n+1,3)=sqrt(max(0d0,(ped1+p(n+1,5))*(ped1-p(n+1,5))))
1178  p(n+1,4)=ped1
1179  p(n+2,1)=0d0
1180  p(n+2,2)=0d0
1181  p(n+2,3)=-p(n+1,3)
1182  p(n+2,4)=p(im,5)-ped1
1183  v(n+1,2)=p(n+1,4)
1184  v(n+2,2)=p(n+2,4)
1185  ELSEIF(nep.GE.3) THEN
1186 C...Rescale all momenta for energy conservation.
1187  loop=0
1188  pes=0d0
1189  pqs=0d0
1190  DO 510 i=1,nep
1191  DO 500 j=1,4
1192  p(n+i,j)=p(ipa(i),j)
1193  500 CONTINUE
1194  pes=pes+p(n+i,4)
1195  pqs=pqs+p(n+i,5)**2/p(n+i,4)
1196  510 CONTINUE
1197  520 loop=loop+1
1198  fac=(ps(5)-pqs)/(pes-pqs)
1199  pes=0d0
1200  pqs=0d0
1201  DO 540 i=1,nep
1202  DO 530 j=1,3
1203  p(n+i,j)=fac*p(n+i,j)
1204  530 CONTINUE
1205  p(n+i,4)=sqrt(p(n+i,5)**2+p(n+i,1)**2+p(n+i,2)**2+p(n+i,3)**2)
1206  v(n+i,2)=p(n+i,4)
1207  pes=pes+p(n+i,4)
1208  pqs=pqs+p(n+i,5)**2/p(n+i,4)
1209  540 CONTINUE
1210  IF(loop.LT.10.AND.abs(pes-ps(5)).GT.1d-12*ps(5)) goto 520
1211 
1212 C...Construct transverse momentum for ordinary branching in shower.
1213  ELSE
1214  zm=v(im,1)
1215  looppt=0
1216  550 looppt=looppt+1
1217  pzm=sqrt(max(0d0,(pem+p(im,5))*(pem-p(im,5))))
1218  pmls=(v(im,5)-v(n+1,5)-v(n+2,5))**2-4d0*v(n+1,5)*v(n+2,5)
1219  IF(pzm.LE.0d0) THEN
1220  pts=0d0
1221  ELSEIF(k(im,2).EQ.21.AND.iabs(k(n+1,2)).LE.10.AND.
1222  & (mstj(44).EQ.3.OR.mstj(44).EQ.5)) THEN
1223  pts=pmls*zm*(1d0-zm)/v(im,5)
1224  ELSEIF(mod(mstj(43),2).EQ.1) THEN
1225  pts=(pem**2*(zm*(1d0-zm)*v(im,5)-(1d0-zm)*v(n+1,5)-
1226  & zm*v(n+2,5))-0.25d0*pmls)/pzm**2
1227  ELSE
1228  pts=pmls*(zm*(1d0-zm)*pem**2/v(im,5)-0.25d0)/pzm**2
1229  ENDIF
1230  IF(pts.LT.0d0.AND.looppt.LT.10) THEN
1231  zm=0.05d0+0.9d0*zm
1232  goto 550
1233  ELSEIF(pts.LT.0d0) THEN
1234  goto 280
1235  ENDIF
1236  pt=sqrt(max(0d0,pts))
1237 
1238 C...Global statistics.
1239  mint(353)=mint(353)+1
1240  vint(353)=vint(353)+pt
1241  IF (mint(353).EQ.1) vint(358)=pt
1242 
1243 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
1244  hazip=0d0
1245  IF(mstj(49).NE.1.AND.mod(mstj(46),2).EQ.1.AND.k(im,2).EQ.21
1246  & .AND.iau.NE.0) THEN
1247  IF(k(igm,3).NE.0) mazip=1
1248  zau=v(igm,1)
1249  IF(iau.EQ.im+1) zau=1d0-v(igm,1)
1250  IF(mazip.EQ.0) zau=0d0
1251  IF(k(igm,2).NE.21) THEN
1252  hazip=2d0*zau/(1d0+zau**2)
1253  ELSE
1254  hazip=(zau/(1d0-zau*(1d0-zau)))**2
1255  ENDIF
1256  IF(k(n+1,2).NE.21) THEN
1257  hazip=hazip*(-2d0*zm*(1d0-zm))/(1d0-2d0*zm*(1d0-zm))
1258  ELSE
1259  hazip=hazip*(zm*(1d0-zm)/(1d0-zm*(1d0-zm)))**2
1260  ENDIF
1261  ENDIF
1262 
1263 C...Find coefficient of azimuthal asymmetry due to soft gluon
1264 C...interference.
1265  hazic=0d0
1266  IF(mstj(49).NE.2.AND.mstj(46).GE.2.AND.(k(n+1,2).EQ.21.OR.
1267  & k(n+2,2).EQ.21).AND.iau.NE.0) THEN
1268  IF(k(igm,3).NE.0) mazic=n+1
1269  IF(k(igm,3).NE.0.AND.k(n+1,2).NE.21) mazic=n+2
1270  IF(k(igm,3).NE.0.AND.k(n+1,2).EQ.21.AND.k(n+2,2).EQ.21.AND.
1271  & zm.GT.0.5d0) mazic=n+2
1272  IF(k(iau,2).EQ.22) mazic=0
1273  zs=zm
1274  IF(mazic.EQ.n+2) zs=1d0-zm
1275  zgm=v(igm,1)
1276  IF(iau.EQ.im-1) zgm=1d0-v(igm,1)
1277  IF(mazic.EQ.0) zgm=1d0
1278  IF(mazic.NE.0) hazic=(p(im,5)/p(igm,5))*
1279  & sqrt((1d0-zs)*(1d0-zgm)/(zs*zgm))
1280  hazic=min(0.95d0,hazic)
1281  ENDIF
1282  ENDIF
1283 
1284 C...Construct energies for ordinary branching in shower.
1285  560 IF(nep.EQ.2.AND.igm.GT.0) THEN
1286  IF(k(im,2).EQ.21.AND.iabs(k(n+1,2)).LE.10.AND.
1287  & (mstj(44).EQ.3.OR.mstj(44).EQ.5)) THEN
1288  p(n+1,4)=0.5d0*(pem*(v(im,5)+v(n+1,5)-v(n+2,5))+
1289  & pzm*sqrt(max(0d0,pmls))*(2d0*zm-1d0))/v(im,5)
1290  ELSEIF(mod(mstj(43),2).EQ.1) THEN
1291  p(n+1,4)=pem*v(im,1)
1292  ELSE
1293  p(n+1,4)=pem*(0.5d0*(v(im,5)-sqrt(pmls)+v(n+1,5)-v(n+2,5))+
1294  & sqrt(pmls)*zm)/v(im,5)
1295  ENDIF
1296 
1297 C...Already predetermined choice of phi angle or not
1298  phi=paru(2)*pyr(0)
1299  IF(mpspd.EQ.1.AND.igm.EQ.ns+1) THEN
1300  ipspd=ip1+im-ns-2
1301  IF(k(ipspd,4).GT.0) THEN
1302  ipsgd1=k(ipspd,4)
1303  IF(im.EQ.ns+2) THEN
1304  phi=pyangl(p(ipsgd1,1),p(ipsgd1,2))
1305  ELSE
1306  phi=pyangl(-p(ipsgd1,1),p(ipsgd1,2))
1307  ENDIF
1308  ENDIF
1309  ELSEIF(mpspd.EQ.1.AND.igm.EQ.ns+2) THEN
1310  ipspd=ip1+im-ns-2
1311  IF(k(ipspd,4).GT.0) THEN
1312  ipsgd1=k(ipspd,4)
1313  phipsm=pyangl(p(ipspd,1),p(ipspd,2))
1314  thepsm=pyangl(p(ipspd,3),sqrt(p(ipspd,1)**2+p(ipspd,2)**2))
1315  CALL pyrobo(ipsgd1,ipsgd1,0d0,-phipsm,0d0,0d0,0d0)
1316  CALL pyrobo(ipsgd1,ipsgd1,-thepsm,0d0,0d0,0d0,0d0)
1317  phi=pyangl(p(ipsgd1,1),p(ipsgd1,2))
1318  CALL pyrobo(ipsgd1,ipsgd1,thepsm,phipsm,0d0,0d0,0d0)
1319  ENDIF
1320  ENDIF
1321 
1322 C...Construct momenta for ordinary branching in shower.
1323  p(n+1,1)=pt*cos(phi)
1324  p(n+1,2)=pt*sin(phi)
1325  IF(k(im,2).EQ.21.AND.iabs(k(n+1,2)).LE.10.AND.
1326  & (mstj(44).EQ.3.OR.mstj(44).EQ.5)) THEN
1327  p(n+1,3)=0.5d0*(pzm*(v(im,5)+v(n+1,5)-v(n+2,5))+
1328  & pem*sqrt(max(0d0,pmls))*(2d0*zm-1d0))/v(im,5)
1329  ELSEIF(pzm.GT.0d0) THEN
1330  p(n+1,3)=0.5d0*(v(n+2,5)-v(n+1,5)-v(im,5)+
1331  & 2d0*pem*p(n+1,4))/pzm
1332  ELSE
1333  p(n+1,3)=0d0
1334  ENDIF
1335  p(n+2,1)=-p(n+1,1)
1336  p(n+2,2)=-p(n+1,2)
1337  p(n+2,3)=pzm-p(n+1,3)
1338  p(n+2,4)=pem-p(n+1,4)
1339  IF(mstj(43).LE.2) THEN
1340  v(n+1,2)=(pem*p(n+1,4)-pzm*p(n+1,3))/p(im,5)
1341  v(n+2,2)=(pem*p(n+2,4)-pzm*p(n+2,3))/p(im,5)
1342  ENDIF
1343  ENDIF
1344 
1345 C...Rotate and boost daughters.
1346  IF(igm.GT.0) THEN
1347  IF(mstj(43).LE.2) THEN
1348  bex=p(igm,1)/p(igm,4)
1349  bey=p(igm,2)/p(igm,4)
1350  bez=p(igm,3)/p(igm,4)
1351  ga=p(igm,4)/p(igm,5)
1352  gabep=ga*(ga*(bex*p(im,1)+bey*p(im,2)+bez*p(im,3))/(1d0+ga)-
1353  & p(im,4))
1354  ELSE
1355  bex=0d0
1356  bey=0d0
1357  bez=0d0
1358  ga=1d0
1359  gabep=0d0
1360  ENDIF
1361  ptimb=sqrt((p(im,1)+gabep*bex)**2+(p(im,2)+gabep*bey)**2)
1362  the=pyangl(p(im,3)+gabep*bez,ptimb)
1363  IF(ptimb.GT.1d-4) THEN
1364  phi=pyangl(p(im,1)+gabep*bex,p(im,2)+gabep*bey)
1365  ELSE
1366  phi=0d0
1367  ENDIF
1368  DO 570 i=n+1,n+2
1369  dp(1)=cos(the)*cos(phi)*p(i,1)-sin(phi)*p(i,2)+
1370  & sin(the)*cos(phi)*p(i,3)
1371  dp(2)=cos(the)*sin(phi)*p(i,1)+cos(phi)*p(i,2)+
1372  & sin(the)*sin(phi)*p(i,3)
1373  dp(3)=-sin(the)*p(i,1)+cos(the)*p(i,3)
1374  dp(4)=p(i,4)
1375  dbp=bex*dp(1)+bey*dp(2)+bez*dp(3)
1376  dgabp=ga*(ga*dbp/(1d0+ga)+dp(4))
1377  p(i,1)=dp(1)+dgabp*bex
1378  p(i,2)=dp(2)+dgabp*bey
1379  p(i,3)=dp(3)+dgabp*bez
1380  p(i,4)=ga*(dp(4)+dbp)
1381  570 CONTINUE
1382  ENDIF
1383 
1384 C...Weight with azimuthal distribution, if required.
1385  IF(mazip.NE.0.OR.mazic.NE.0) THEN
1386  DO 580 j=1,3
1387  dpt(1,j)=p(im,j)
1388  dpt(2,j)=p(iau,j)
1389  dpt(3,j)=p(n+1,j)
1390  580 CONTINUE
1391  dpma=dpt(1,1)*dpt(2,1)+dpt(1,2)*dpt(2,2)+dpt(1,3)*dpt(2,3)
1392  dpmd=dpt(1,1)*dpt(3,1)+dpt(1,2)*dpt(3,2)+dpt(1,3)*dpt(3,3)
1393  dpmm=dpt(1,1)**2+dpt(1,2)**2+dpt(1,3)**2
1394  DO 590 j=1,3
1395  dpt(4,j)=dpt(2,j)-dpma*dpt(1,j)/max(1d-10,dpmm)
1396  dpt(5,j)=dpt(3,j)-dpmd*dpt(1,j)/max(1d-10,dpmm)
1397  590 CONTINUE
1398  dpt(4,4)=sqrt(dpt(4,1)**2+dpt(4,2)**2+dpt(4,3)**2)
1399  dpt(5,4)=sqrt(dpt(5,1)**2+dpt(5,2)**2+dpt(5,3)**2)
1400  IF(min(dpt(4,4),dpt(5,4)).GT.0.1d0*parj(82)) THEN
1401  cad=(dpt(4,1)*dpt(5,1)+dpt(4,2)*dpt(5,2)+
1402  & dpt(4,3)*dpt(5,3))/(dpt(4,4)*dpt(5,4))
1403  IF(mazip.NE.0) THEN
1404  IF(1d0+hazip*(2d0*cad**2-1d0).LT.pyr(0)*(1d0+abs(hazip)))
1405  & goto 560
1406  ENDIF
1407  IF(mazic.NE.0) THEN
1408  IF(mazic.EQ.n+2) cad=-cad
1409  IF((1d0-hazic)*(1d0-hazic*cad)/(1d0+hazic**2-2d0*hazic*cad)
1410  & .LT.pyr(0)) goto 560
1411  ENDIF
1412  ENDIF
1413  ENDIF
1414 
1415 C...Azimuthal anisotropy due to interference with initial state partons.
1416  IF(mod(miis,2).EQ.1.AND.igm.EQ.ns+1.AND.(k(n+1,2).EQ.21.OR.
1417  &k(n+2,2).EQ.21)) THEN
1418  iii=im-ns-1
1419  IF(isii(iii).GE.1) THEN
1420  iaziid=n+1
1421  IF(k(n+1,2).NE.21) iaziid=n+2
1422  IF(k(n+1,2).EQ.21.AND.k(n+2,2).EQ.21.AND.
1423  & p(n+1,4).GT.p(n+2,4)) iaziid=n+2
1424  theiid=pyangl(p(iaziid,3),sqrt(p(iaziid,1)**2+p(iaziid,2)**2))
1425  IF(iii.EQ.2) theiid=paru(1)-theiid
1426  phiiid=pyangl(p(iaziid,1),p(iaziid,2))
1427  hazii=min(0.95d0,theiid/theiis(iii,isii(iii)))
1428  cad=cos(phiiid-phiiis(iii,isii(iii)))
1429  phirel=abs(phiiid-phiiis(iii,isii(iii)))
1430  IF(phirel.GT.paru(1)) phirel=paru(2)-phirel
1431  IF((1d0-hazii)*(1d0-hazii*cad)/(1d0+hazii**2-2d0*hazii*cad)
1432  & .LT.pyr(0)) goto 560
1433  ENDIF
1434  ENDIF
1435 
1436 C...Continue loop over partons that may branch, until none left.
1437  IF(igm.GE.0) k(im,1)=14
1438  n=n+nep
1439  nep=2
1440  IF(n.GT.mstu(4)-mstu(32)-10) THEN
1441  CALL pyerrm(11,'(PYSHOW:) no more memory left in PYJETS')
1442  IF(mstu(21).GE.1) n=ns
1443  IF(mstu(21).GE.1) RETURN
1444  ENDIF
1445  goto 290
1446 
1447 C...Set information on imagined shower initiator.
1448  600 IF(npa.GE.2) THEN
1449  k(ns+1,1)=11
1450  k(ns+1,2)=94
1451  k(ns+1,3)=ip1
1452  IF(ip2.GT.0.AND.ip2.LT.ip1) k(ns+1,3)=ip2
1453  k(ns+1,4)=ns+2
1454  k(ns+1,5)=ns+1+npa
1455  iim=1
1456  ELSE
1457  iim=0
1458  ENDIF
1459 
1460 C...Reconstruct string drawing information.
1461  DO 610 i=ns+1+iim,n
1462  kq=kchg(pycomp(k(i,2)),2)
1463  IF(k(i,1).LE.10.AND.k(i,2).EQ.22) THEN
1464  k(i,1)=1
1465  ELSEIF(k(i,1).LE.10.AND.iabs(k(i,2)).GE.11.AND.
1466  & iabs(k(i,2)).LE.18) THEN
1467  k(i,1)=1
1468  ELSEIF(k(i,1).LE.10) THEN
1469  k(i,4)=mstu(5)*(k(i,4)/mstu(5))
1470  k(i,5)=mstu(5)*(k(i,5)/mstu(5))
1471  ELSEIF(k(mod(k(i,4),mstu(5))+1,2).NE.22) THEN
1472  id1=mod(k(i,4),mstu(5))
1473  IF(kq.EQ.1.AND.k(i,2).GT.0) id1=mod(k(i,4),mstu(5))+1
1474  IF(kq.EQ.2.AND.(k(id1,2).EQ.21.OR.k(id1+1,2).EQ.21).AND.
1475  & pyr(0).GT.0.5d0) id1=mod(k(i,4),mstu(5))+1
1476  id2=2*mod(k(i,4),mstu(5))+1-id1
1477  k(i,4)=mstu(5)*(k(i,4)/mstu(5))+id1
1478  k(i,5)=mstu(5)*(k(i,5)/mstu(5))+id2
1479  k(id1,4)=k(id1,4)+mstu(5)*i
1480  k(id1,5)=k(id1,5)+mstu(5)*id2
1481  k(id2,4)=k(id2,4)+mstu(5)*id1
1482  k(id2,5)=k(id2,5)+mstu(5)*i
1483  ELSE
1484  id1=mod(k(i,4),mstu(5))
1485  id2=id1+1
1486  k(i,4)=mstu(5)*(k(i,4)/mstu(5))+id1
1487  k(i,5)=mstu(5)*(k(i,5)/mstu(5))+id1
1488  IF(kq.EQ.1.OR.k(id1,1).GE.11) THEN
1489  k(id1,4)=k(id1,4)+mstu(5)*i
1490  k(id1,5)=k(id1,5)+mstu(5)*i
1491  ELSE
1492  k(id1,4)=0
1493  k(id1,5)=0
1494  ENDIF
1495  k(id2,4)=0
1496  k(id2,5)=0
1497  ENDIF
1498  610 CONTINUE
1499 
1500 C...Transformation from CM frame.
1501  IF(npa.EQ.1) THEN
1502  the=pyangl(p(ipa(1),3),sqrt(p(ipa(1),1)**2+p(ipa(1),2)**2))
1503  phi=pyangl(p(ipa(1),1),p(ipa(1),2))
1504  mstu(33)=1
1505  CALL pyrobo(ns+1,n,the,phi,0d0,0d0,0d0)
1506  ELSEIF(npa.EQ.2) THEN
1507  bex=ps(1)/ps(4)
1508  bey=ps(2)/ps(4)
1509  bez=ps(3)/ps(4)
1510  ga=ps(4)/ps(5)
1511  gabep=ga*(ga*(bex*p(ipa(1),1)+bey*p(ipa(1),2)+bez*p(ipa(1),3))
1512  & /(1d0+ga)-p(ipa(1),4))
1513  the=pyangl(p(ipa(1),3)+gabep*bez,sqrt((p(ipa(1),1)
1514  & +gabep*bex)**2+(p(ipa(1),2)+gabep*bey)**2))
1515  phi=pyangl(p(ipa(1),1)+gabep*bex,p(ipa(1),2)+gabep*bey)
1516  mstu(33)=1
1517  CALL pyrobo(ns+1,n,the,phi,bex,bey,bez)
1518  ELSE
1519  CALL pyrobo(ipa(1),ipa(npa),0d0,0d0,ps(1)/ps(4),ps(2)/ps(4),
1520  & ps(3)/ps(4))
1521  mstu(33)=1
1522  CALL pyrobo(ns+1,n,0d0,0d0,ps(1)/ps(4),ps(2)/ps(4),ps(3)/ps(4))
1523  ENDIF
1524 
1525 C...Decay vertex of shower.
1526  DO 630 i=ns+1,n
1527  DO 620 j=1,5
1528  v(i,j)=v(ip1,j)
1529  620 CONTINUE
1530  630 CONTINUE
1531 
1532 C...Delete trivial shower, else connect initiators.
1533  IF(n.LE.ns+npa+iim) THEN
1534  n=ns
1535  ELSE
1536  DO 640 ip=1,npa
1537  k(ipa(ip),1)=14
1538  k(ipa(ip),4)=k(ipa(ip),4)+ns+iim+ip
1539  k(ipa(ip),5)=k(ipa(ip),5)+ns+iim+ip
1540  k(ns+iim+ip,3)=ipa(ip)
1541  IF(iim.EQ.1.AND.mstu(16).NE.2) k(ns+iim+ip,3)=ns+1
1542  IF(k(ns+iim+ip,1).NE.1) THEN
1543  k(ns+iim+ip,4)=mstu(5)*ipa(ip)+k(ns+iim+ip,4)
1544  k(ns+iim+ip,5)=mstu(5)*ipa(ip)+k(ns+iim+ip,5)
1545  ENDIF
1546  640 CONTINUE
1547  ENDIF
1548 
1549  RETURN
1550  END