ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhiresd.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhiresd.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhiresd
5 
6 C...Allows resonances to decay (including parton showers for hadronic
7 C...channels).
8  IMPLICIT DOUBLE PRECISION(d)
9  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
10  SAVE /lujets/
11  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
12  SAVE /ludat1/
13  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
14  SAVE /ludat2/
15  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
16  SAVE /ludat3/
17  common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
18  SAVE /pyhisubs/
19  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
20  SAVE /pyhipars/
21  common/pyhiint1/mint(400),vint(400)
22  SAVE /pyhiint1/
23  common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
24  SAVE /pyhiint2/
25  common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3)
26  SAVE /pyhiint4/
27  dimension iref(10,6),kdcy(2),kfl1(2),kfl2(2),nsd(2),ilin(6),
28  &coup(6,4),pk(6,4),pkk(6,6),cthe(2),phi(2),wdtp(0:40),
29  &wdte(0:40,0:5)
30  COMPLEX fgk,ha(6,6),hc(6,6)
31 
32 C...The F, Xi and Xj functions of Gunion and Kunszt
33 C...(Phys. Rev. D33, 665, plus errata from the authors).
34  fgk(i1,i2,i3,i4,i5,i6)=4.*ha(i1,i3)*hc(i2,i6)*(ha(i1,i5)*
35  &hc(i1,i4)+ha(i3,i5)*hc(i3,i4))
36  digk(dt,du)=-4.*d34*d56+dt*(3.*dt+4.*du)+dt**2*(dt*du/(d34*d56)-
37  &2.*(1./d34+1./d56)*(dt+du)+2.*(d34/d56+d56/d34))
38  djgk(dt,du)=8.*(d34+d56)**2-8.*(d34+d56)*(dt+du)-6.*dt*du-
39  &2.*dt*du*(dt*du/(d34*d56)-2.*(1./d34+1./d56)*(dt+du)+
40  &2.*(d34/d56+d56/d34))
41 
42 C...Define initial two objects, initialize loop.
43  isub=mint(1)
44  sh=vint(44)
45  iref(1,5)=0
46  iref(1,6)=0
47  IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN
48  iref(1,1)=mint(84)+2+iset(isub)
49  iref(1,2)=0
50  iref(1,3)=mint(83)+6+iset(isub)
51  iref(1,4)=0
52  ELSEIF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN
53  iref(1,1)=mint(84)+1+iset(isub)
54  iref(1,2)=mint(84)+2+iset(isub)
55  iref(1,3)=mint(83)+5+iset(isub)
56  iref(1,4)=mint(83)+6+iset(isub)
57  ENDIF
58  np=1
59  ip=0
60  100 ip=ip+1
61  ninh=0
62 
63 C...Loop over one/two resonances; reset decay rates.
64  jtmax=2
65  IF(ip.EQ.1.AND.(iset(isub).EQ.1.OR.iset(isub).EQ.3)) jtmax=1
66  DO 140 jt=1,jtmax
67  kdcy(jt)=0
68  kfl1(jt)=0
69  kfl2(jt)=0
70  nsd(jt)=iref(ip,jt)
71  id=iref(ip,jt)
72  IF(id.EQ.0) goto 140
73  kfa=iabs(k(id,2))
74  IF(kfa.LT.23.OR.kfa.GT.40) goto 140
75  IF(mdcy(kfa,1).NE.0) THEN
76  IF(isub.EQ.1.OR.isub.EQ.141) mint(61)=1
77  CALL pyhiwidt(kfa,p(id,5),wdtp,wdte)
78  IF(kchg(kfa,3).EQ.0) THEN
79  ipm=2
80  ELSE
81  ipm=(5+isign(1,k(id,2)))/2
82  ENDIF
83  IF(jtmax.EQ.1.OR.iabs(k(iref(ip,1),2)).NE.iabs(k(iref(ip,2),2)))
84  & THEN
85  i12=4
86  ELSE
87  IF(jt.EQ.1) i12=int(4.5+rlu(0))
88  i12=9-i12
89  ENDIF
90  rkfl=(wdte(0,1)+wdte(0,ipm)+wdte(0,i12))*rlu(0)
91  DO 120 i=1,mdcy(kfa,3)
92  idc=i+mdcy(kfa,2)-1
93  kfl1(jt)=kfdp(idc,1)*isign(1,k(id,2))
94  kfl2(jt)=kfdp(idc,2)*isign(1,k(id,2))
95  rkfl=rkfl-(wdte(i,1)+wdte(i,ipm)+wdte(i,i12))
96  IF(rkfl.LE.0.) goto 130
97  120 CONTINUE
98  130 CONTINUE
99  ENDIF
100 
101 C...Summarize result on decay channel chosen.
102  IF((kfa.EQ.23.OR.kfa.EQ.24).AND.kfl1(jt).EQ.0) ninh=ninh+1
103  IF(kfl1(jt).EQ.0) goto 140
104  kdcy(jt)=2
105  IF(iabs(kfl1(jt)).LE.10.OR.kfl1(jt).EQ.21) kdcy(jt)=1
106  IF((iabs(kfl1(jt)).GE.23.AND.iabs(kfl1(jt)).LE.25).OR.
107  &(iabs(kfl1(jt)).EQ.37)) kdcy(jt)=3
108  nsd(jt)=n
109 
110 C...Fill decay products, prepared for parton showers for quarks.
111  IF(kdcy(jt).EQ.1) THEN
112  CALL lu2ent(-(n+1),kfl1(jt),kfl2(jt),p(id,5))
113  ELSE
114  CALL lu2ent(n+1,kfl1(jt),kfl2(jt),p(id,5))
115  ENDIF
116  IF(jtmax.EQ.1) THEN
117  cthe(jt)=vint(13)+(vint(33)-vint(13)+vint(34)-vint(14))*rlu(0)
118  IF(cthe(jt).GT.vint(33)) cthe(jt)=cthe(jt)+vint(14)-vint(33)
119  phi(jt)=vint(24)
120  ELSE
121  cthe(jt)=2.*rlu(0)-1.
122  phi(jt)=paru(2)*rlu(0)
123  ENDIF
124  140 CONTINUE
125  IF(mint(3).EQ.1.AND.ip.EQ.1) THEN
126  mint(25)=kfl1(1)
127  mint(26)=kfl2(1)
128  ENDIF
129  IF(jtmax.EQ.1.AND.kdcy(1).EQ.0) goto 530
130  IF(jtmax.EQ.2.AND.kdcy(1).EQ.0.AND.kdcy(2).EQ.0) goto 530
131  IF(mstp(45).LE.0.OR.iref(ip,2).EQ.0.OR.ninh.GE.1) goto 500
132  IF(k(iref(1,1),2).EQ.25.AND.ip.EQ.1) goto 500
133  IF(k(iref(1,1),2).EQ.25.AND.kdcy(1)*kdcy(2).EQ.0) goto 500
134 
135 C...Order incoming partons and outgoing resonances.
136  ilin(1)=mint(84)+1
137  IF(k(mint(84)+1,2).GT.0) ilin(1)=mint(84)+2
138  IF(k(ilin(1),2).EQ.21) ilin(1)=2*mint(84)+3-ilin(1)
139  ilin(2)=2*mint(84)+3-ilin(1)
140  imin=1
141  IF(iref(ip,5).EQ.25) imin=3
142  imax=2
143  iord=1
144  IF(k(iref(ip,1),2).EQ.23) iord=2
145  IF(k(iref(ip,1),2).EQ.24.AND.k(iref(ip,2),2).EQ.-24) iord=2
146  IF(iabs(k(iref(ip,iord),2)).EQ.25) iord=3-iord
147  IF(kdcy(iord).EQ.0) iord=3-iord
148 
149 C...Order decay products of resonances.
150  DO 390 jt=iord,3-iord,3-2*iord
151  IF(kdcy(jt).EQ.0) THEN
152  ilin(imax+1)=nsd(jt)
153  imax=imax+1
154  ELSEIF(k(nsd(jt)+1,2).GT.0) THEN
155  ilin(imax+1)=n+2*jt-1
156  ilin(imax+2)=n+2*jt
157  imax=imax+2
158  k(n+2*jt-1,2)=k(nsd(jt)+1,2)
159  k(n+2*jt,2)=k(nsd(jt)+2,2)
160  ELSE
161  ilin(imax+1)=n+2*jt
162  ilin(imax+2)=n+2*jt-1
163  imax=imax+2
164  k(n+2*jt-1,2)=k(nsd(jt)+1,2)
165  k(n+2*jt,2)=k(nsd(jt)+2,2)
166  ENDIF
167  390 CONTINUE
168 
169 C...Find charge, isospin, left- and righthanded couplings.
170  xw=paru(102)
171  DO 410 i=imin,imax
172  DO 400 j=1,4
173  400 coup(i,j)=0.
174  kfa=iabs(k(ilin(i),2))
175  IF(kfa.GT.20) goto 410
176  coup(i,1)=luchge(kfa)/3.
177  coup(i,2)=(-1)**mod(kfa,2)
178  coup(i,4)=-2.*coup(i,1)*xw
179  coup(i,3)=coup(i,2)+coup(i,4)
180  410 CONTINUE
181  sqmz=pmas(23,1)**2
182  gzmz=pmas(23,1)*pmas(23,2)
183  sqmw=pmas(24,1)**2
184  gzmw=pmas(24,1)*pmas(24,2)
185  sqmzp=pmas(32,1)**2
186  gzmzp=pmas(32,1)*pmas(32,2)
187 
188 C...Select random angles; construct massless four-vectors.
189  420 DO 430 i=n+1,n+4
190  k(i,1)=1
191  DO 430 j=1,5
192  430 p(i,j)=0.
193  DO 440 jt=1,jtmax
194  IF(kdcy(jt).EQ.0) goto 440
195  id=iref(ip,jt)
196  p(n+2*jt-1,3)=0.5*p(id,5)
197  p(n+2*jt-1,4)=0.5*p(id,5)
198  p(n+2*jt,3)=-0.5*p(id,5)
199  p(n+2*jt,4)=0.5*p(id,5)
200  cthe(jt)=2.*rlu(0)-1.
201  phi(jt)=paru(2)*rlu(0)
202  CALL ludbrb(n+2*jt-1,n+2*jt,acos(cthe(jt)),phi(jt),
203  &dble(p(id,1)/p(id,4)),dble(p(id,2)/p(id,4)),dble(p(id,3)/p(id,4)))
204  440 CONTINUE
205 
206 C...Store incoming and outgoing momenta, with random rotation to
207 C...avoid accidental zeroes in HA expressions.
208  DO 450 i=1,imax
209  k(n+4+i,1)=1
210  p(n+4+i,4)=sqrt(p(ilin(i),1)**2+p(ilin(i),2)**2+p(ilin(i),3)**2+
211  &p(ilin(i),5)**2)
212  p(n+4+i,5)=p(ilin(i),5)
213  DO 450 j=1,3
214  450 p(n+4+i,j)=p(ilin(i),j)
215  therr=acos(2.*rlu(0)-1.)
216  phirr=paru(2)*rlu(0)
217  CALL ludbrb(n+5,n+4+imax,therr,phirr,0d0,0d0,0d0)
218  DO 460 i=1,imax
219  DO 460 j=1,4
220  460 pk(i,j)=p(n+4+i,j)
221 
222 C...Calculate internal products.
223  IF(isub.EQ.22.OR.isub.EQ.23.OR.isub.EQ.25) THEN
224  DO 470 i1=imin,imax-1
225  DO 470 i2=i1+1,imax
226  ha(i1,i2)=sqrt((pk(i1,4)-pk(i1,3))*(pk(i2,4)+pk(i2,3))/
227  & (1e-20+pk(i1,1)**2+pk(i1,2)**2))*cmplx(pk(i1,1),pk(i1,2))-
228  & sqrt((pk(i1,4)+pk(i1,3))*(pk(i2,4)-pk(i2,3))/
229  & (1e-20+pk(i2,1)**2+pk(i2,2)**2))*cmplx(pk(i2,1),pk(i2,2))
230  hc(i1,i2)=conjg(ha(i1,i2))
231  IF(i1.LE.2) ha(i1,i2)=cmplx(0.,1.)*ha(i1,i2)
232  IF(i1.LE.2) hc(i1,i2)=cmplx(0.,1.)*hc(i1,i2)
233  ha(i2,i1)=-ha(i1,i2)
234  470 hc(i2,i1)=-hc(i1,i2)
235  ENDIF
236  DO 480 i=1,2
237  DO 480 j=1,4
238  480 pk(i,j)=-pk(i,j)
239  DO 490 i1=imin,imax-1
240  DO 490 i2=i1+1,imax
241  pkk(i1,i2)=2.*(pk(i1,4)*pk(i2,4)-pk(i1,1)*pk(i2,1)-
242  &pk(i1,2)*pk(i2,2)-pk(i1,3)*pk(i2,3))
243  490 pkk(i2,i1)=pkk(i1,i2)
244 
245  IF(iref(ip,5).EQ.25) THEN
246 C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons
247  wt=16.*pkk(3,5)*pkk(4,6)
248  IF(ip.EQ.1) wtmax=sh**2
249  IF(ip.GE.2) wtmax=p(iref(ip,6),5)**4
250 
251  ELSEIF(isub.EQ.1) THEN
252  IF(kfa.NE.37) THEN
253 C...Angular weight for gamma*/Z0 -> 2 quarks/leptons
254  ei=kchg(iabs(mint(15)),1)/3.
255  ai=sign(1.,ei+0.1)
256  vi=ai-4.*ei*xw
257  ef=kchg(kfa,1)/3.
258  af=sign(1.,ef+0.1)
259  vf=af-4.*ef*xw
260  gg=1.
261  gz=1./(8.*xw*(1.-xw))*sh*(sh-sqmz)/((sh-sqmz)**2+gzmz**2)
262  zz=1./(16.*xw*(1.-xw))**2*sh**2/((sh-sqmz)**2+gzmz**2)
263  IF(mstp(43).EQ.1) THEN
264 C...Only gamma* production included
265  gz=0.
266  zz=0.
267  ELSEIF(mstp(43).EQ.2) THEN
268 C...Only Z0 production included
269  gg=0.
270  gz=0.
271  ENDIF
272  asym=2.*(ei*ai*gz*ef*af+4.*vi*ai*zz*vf*af)/(ei**2*gg*ef**2+
273  & ei*vi*gz*ef*vf+(vi**2+ai**2)*zz*(vf**2+af**2))
274  wt=1.+asym*cthe(jt)+cthe(jt)**2
275  wtmax=2.+abs(asym)
276  ELSE
277 C...Angular weight for gamma*/Z0 -> H+ + H-
278  wt=1.-cthe(jt)**2
279  wtmax=1.
280  ENDIF
281 
282  ELSEIF(isub.EQ.2) THEN
283 C...Angular weight for W+/- -> 2 quarks/leptons
284  wt=(1.+cthe(jt))**2
285  wtmax=4.
286 
287  ELSEIF(isub.EQ.15.OR.isub.EQ.19) THEN
288 C...Angular weight for f + fb -> gluon/gamma + Z0 ->
289 C...-> gluon/gamma + 2 quarks/leptons
290  wt=((coup(1,3)*coup(3,3))**2+(coup(1,4)*coup(3,4))**2)*
291  & (pkk(1,3)**2+pkk(2,4)**2)+((coup(1,3)*coup(3,4))**2+
292  & (coup(1,4)*coup(3,3))**2)*(pkk(1,4)**2+pkk(2,3)**2)
293  wtmax=(coup(1,3)**2+coup(1,4)**2)*(coup(3,3)**2+coup(3,4)**2)*
294  & ((pkk(1,3)+pkk(1,4))**2+(pkk(2,3)+pkk(2,4))**2)
295 
296  ELSEIF(isub.EQ.16.OR.isub.EQ.20) THEN
297 C...Angular weight for f + fb' -> gluon/gamma + W+/- ->
298 C...-> gluon/gamma + 2 quarks/leptons
299  wt=pkk(1,3)**2+pkk(2,4)**2
300  wtmax=(pkk(1,3)+pkk(1,4))**2+(pkk(2,3)+pkk(2,4))**2
301 
302  ELSEIF(isub.EQ.22) THEN
303 C...Angular weight for f + fb -> Z0 + Z0 -> 4 quarks/leptons
304  s34=p(iref(ip,iord),5)**2
305  s56=p(iref(ip,3-iord),5)**2
306  ti=pkk(1,3)+pkk(1,4)+s34
307  ui=pkk(1,5)+pkk(1,6)+s56
308  wt=coup(1,3)**4*((coup(3,3)*coup(5,3)*abs(fgk(1,2,3,4,5,6)/
309  & ti+fgk(1,2,5,6,3,4)/ui))**2+(coup(3,4)*coup(5,3)*abs(
310  & fgk(1,2,4,3,5,6)/ti+fgk(1,2,5,6,4,3)/ui))**2+(coup(3,3)*
311  & coup(5,4)*abs(fgk(1,2,3,4,6,5)/ti+fgk(1,2,6,5,3,4)/ui))**2+
312  & (coup(3,4)*coup(5,4)*abs(fgk(1,2,4,3,6,5)/ti+fgk(1,2,6,5,4,3)/
313  & ui))**2)+coup(1,4)**4*((coup(3,3)*coup(5,3)*abs(
314  & fgk(2,1,5,6,3,4)/ti+fgk(2,1,3,4,5,6)/ui))**2+(coup(3,4)*
315  & coup(5,3)*abs(fgk(2,1,6,5,3,4)/ti+fgk(2,1,3,4,6,5)/ui))**2+
316  & (coup(3,3)*coup(5,4)*abs(fgk(2,1,5,6,4,3)/ti+fgk(2,1,4,3,5,6)/
317  & ui))**2+(coup(3,4)*coup(5,4)*abs(fgk(2,1,6,5,4,3)/ti+
318  & fgk(2,1,4,3,6,5)/ui))**2)
319  wtmax=4.*s34*s56*(coup(1,3)**4+coup(1,4)**4)*(coup(3,3)**2+
320  & coup(3,4)**2)*(coup(5,3)**2+coup(5,4)**2)*4.*(ti/ui+ui/ti+
321  & 2.*sh*(s34+s56)/(ti*ui)-s34*s56*(1./ti**2+1./ui**2))
322 
323  ELSEIF(isub.EQ.23) THEN
324 C...Angular weight for f + fb' -> Z0 + W +/- -> 4 quarks/leptons
325  d34=p(iref(ip,iord),5)**2
326  d56=p(iref(ip,3-iord),5)**2
327  dt=pkk(1,3)+pkk(1,4)+d34
328  du=pkk(1,5)+pkk(1,6)+d56
329  cawz=coup(2,3)/sngl(dt)-2.*(1.-xw)*coup(1,2)/(sh-sqmw)
330  cbwz=coup(1,3)/sngl(du)+2.*(1.-xw)*coup(1,2)/(sh-sqmw)
331  wt=coup(5,3)**2*abs(cawz*fgk(1,2,3,4,5,6)+cbwz*
332  & fgk(1,2,5,6,3,4))**2+coup(5,4)**2*abs(cawz*
333  & fgk(1,2,3,4,6,5)+cbwz*fgk(1,2,6,5,3,4))**2
334  wtmax=4.*d34*d56*(coup(5,3)**2+coup(5,4)**2)*(cawz**2*
335  & digk(dt,du)+cbwz**2*digk(du,dt)+cawz*cbwz*djgk(dt,du))
336 
337  ELSEIF(isub.EQ.24) THEN
338 C...Angular weight for f + fb -> Z0 + H0 -> 2 quarks/leptons + H0
339  wt=((coup(1,3)*coup(3,3))**2+(coup(1,4)*coup(3,4))**2)*
340  & pkk(1,3)*pkk(2,4)+((coup(1,3)*coup(3,4))**2+(coup(1,4)*
341  & coup(3,3))**2)*pkk(1,4)*pkk(2,3)
342  wtmax=(coup(1,3)**2+coup(1,4)**2)*(coup(3,3)**2+coup(3,4)**2)*
343  & (pkk(1,3)+pkk(1,4))*(pkk(2,3)+pkk(2,4))
344 
345  ELSEIF(isub.EQ.25) THEN
346 C...Angular weight for f + fb -> W+ + W- -> 4 quarks/leptons
347  d34=p(iref(ip,iord),5)**2
348  d56=p(iref(ip,3-iord),5)**2
349  dt=pkk(1,3)+pkk(1,4)+d34
350  du=pkk(1,5)+pkk(1,6)+d56
351  cdww=(coup(1,3)*sqmz/(sh-sqmz)+coup(1,2))/sh
352  caww=cdww+0.5*(coup(1,2)+1.)/sngl(dt)
353  cbww=cdww+0.5*(coup(1,2)-1.)/sngl(du)
354  ccww=coup(1,4)*sqmz/(sh-sqmz)/sh
355  wt=abs(caww*fgk(1,2,3,4,5,6)-cbww*fgk(1,2,5,6,3,4))**2+
356  & ccww**2*abs(fgk(2,1,5,6,3,4)-fgk(2,1,3,4,5,6))**2
357  wtmax=4.*d34*d56*(caww**2*digk(dt,du)+cbww**2*digk(du,dt)-caww*
358  & cbww*djgk(dt,du)+ccww**2*(digk(dt,du)+digk(du,dt)-djgk(dt,du)))
359 
360  ELSEIF(isub.EQ.26) THEN
361 C...Angular weight for f + fb' -> W+/- + H0 -> 2 quarks/leptons + H0
362  wt=pkk(1,3)*pkk(2,4)
363  wtmax=(pkk(1,3)+pkk(1,4))*(pkk(2,3)+pkk(2,4))
364 
365  ELSEIF(isub.EQ.30) THEN
366 C...Angular weight for f + g -> f + Z0 -> f + 2 quarks/leptons
367  IF(k(ilin(1),2).GT.0) wt=((coup(1,3)*coup(3,3))**2+
368  & (coup(1,4)*coup(3,4))**2)*(pkk(1,4)**2+pkk(3,5)**2)+
369  & ((coup(1,3)*coup(3,4))**2+(coup(1,4)*coup(3,3))**2)*
370  & (pkk(1,3)**2+pkk(4,5)**2)
371  IF(k(ilin(1),2).LT.0) wt=((coup(1,3)*coup(3,3))**2+
372  & (coup(1,4)*coup(3,4))**2)*(pkk(1,3)**2+pkk(4,5)**2)+
373  & ((coup(1,3)*coup(3,4))**2+(coup(1,4)*coup(3,3))**2)*
374  & (pkk(1,4)**2+pkk(3,5)**2)
375  wtmax=(coup(1,3)**2+coup(1,4)**2)*(coup(3,3)**2+coup(3,4)**2)*
376  & ((pkk(1,3)+pkk(1,4))**2+(pkk(3,5)+pkk(4,5))**2)
377 
378  ELSEIF(isub.EQ.31) THEN
379 C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons
380  IF(k(ilin(1),2).GT.0) wt=pkk(1,4)**2+pkk(3,5)**2
381  IF(k(ilin(1),2).LT.0) wt=pkk(1,3)**2+pkk(4,5)**2
382  wtmax=(pkk(1,3)+pkk(1,4))**2+(pkk(3,5)+pkk(4,5))**2
383 
384  ELSEIF(isub.EQ.141) THEN
385 C...Angular weight for gamma*/Z0/Z'0 -> 2 quarks/leptons
386  ei=kchg(iabs(mint(15)),1)/3.
387  ai=sign(1.,ei+0.1)
388  vi=ai-4.*ei*xw
389  api=sign(1.,ei+0.1)
390  vpi=api-4.*ei*xw
391  ef=kchg(kfa,1)/3.
392  af=sign(1.,ef+0.1)
393  vf=af-4.*ef*xw
394  apf=sign(1.,ef+0.1)
395  vpf=apf-4.*ef*xw
396  gg=1.
397  gz=1./(8.*xw*(1.-xw))*sh*(sh-sqmz)/((sh-sqmz)**2+gzmz**2)
398  gzp=1./(8.*xw*(1.-xw))*sh*(sh-sqmzp)/((sh-sqmzp)**2+gzmzp**2)
399  zz=1./(16.*xw*(1.-xw))**2*sh**2/((sh-sqmz)**2+gzmz**2)
400  zzp=2./(16.*xw*(1.-xw))**2*
401  & sh**2*((sh-sqmz)*(sh-sqmzp)+gzmz*gzmzp)/
402  & (((sh-sqmz)**2+gzmz**2)*((sh-sqmzp)**2+gzmzp**2))
403  zpzp=1./(16.*xw*(1.-xw))**2*sh**2/((sh-sqmzp)**2+gzmzp**2)
404  IF(mstp(44).EQ.1) THEN
405 C...Only gamma* production included
406  gz=0.
407  gzp=0.
408  zz=0.
409  zzp=0.
410  zpzp=0.
411  ELSEIF(mstp(44).EQ.2) THEN
412 C...Only Z0 production included
413  gg=0.
414  gz=0.
415  gzp=0.
416  zzp=0.
417  zpzp=0.
418  ELSEIF(mstp(44).EQ.3) THEN
419 C...Only Z'0 production included
420  gg=0.
421  gz=0.
422  gzp=0.
423  zz=0.
424  zzp=0.
425  ELSEIF(mstp(44).EQ.4) THEN
426 C...Only gamma*/Z0 production included
427  gzp=0.
428  zzp=0.
429  zpzp=0.
430  ELSEIF(mstp(44).EQ.5) THEN
431 C...Only gamma*/Z'0 production included
432  gz=0.
433  zz=0.
434  zzp=0.
435  ELSEIF(mstp(44).EQ.6) THEN
436 C...Only Z0/Z'0 production included
437  gg=0.
438  gz=0.
439  gzp=0.
440  ENDIF
441  asym=2.*(ei*ai*gz*ef*af+ei*api*gzp*ef*apf+4.*vi*ai*zz*vf*af+
442  & (vi*api+vpi*ai)*zzp*(vf*apf+vpf*af)+4.*vpi*api*zpzp*vpf*apf)/
443  & (ei**2*gg*ef**2+ei*vi*gz*ef*vf+ei*vpi*gzp*ef*vpf+
444  & (vi**2+ai**2)*zz*(vf**2+af**2)+(vi*vpi+ai*api)*zzp*
445  & (vf*vpf+af*apf)+(vpi**2+api**2)*zpzp*(vpf**2+apf**2))
446  wt=1.+asym*cthe(jt)+cthe(jt)**2
447  wtmax=2.+abs(asym)
448 
449  ELSE
450  wt=1.
451  wtmax=1.
452  ENDIF
453 C...Obtain correct angular distribution by rejection techniques.
454  IF(wt.LT.rlu(0)*wtmax) goto 420
455 
456 C...Construct massive four-vectors using angles chosen. Mark decayed
457 C...resonances, add documentation lines. Shower evolution.
458  500 DO 520 jt=1,jtmax
459  IF(kdcy(jt).EQ.0) goto 520
460  id=iref(ip,jt)
461  CALL ludbrb(nsd(jt)+1,nsd(jt)+2,acos(cthe(jt)),phi(jt),
462  &dble(p(id,1)/p(id,4)),dble(p(id,2)/p(id,4)),dble(p(id,3)/p(id,4)))
463  k(id,1)=k(id,1)+10
464  k(id,4)=nsd(jt)+1
465  k(id,5)=nsd(jt)+2
466  idoc=mint(83)+mint(4)
467  DO 510 i=nsd(jt)+1,nsd(jt)+2
468  mint(4)=mint(4)+1
469  i1=mint(83)+mint(4)
470  k(i,3)=i1
471  k(i1,1)=21
472  k(i1,2)=k(i,2)
473  k(i1,3)=iref(ip,jt+2)
474  DO 510 j=1,5
475  510 p(i1,j)=p(i,j)
476  IF(jtmax.EQ.1) THEN
477  mint(7)=mint(83)+6+2*iset(isub)
478  mint(8)=mint(83)+7+2*iset(isub)
479  ENDIF
480  IF(mstp(71).GE.1.AND.kdcy(jt).EQ.1) CALL lushow(nsd(jt)+1,
481  &nsd(jt)+2,p(id,5))
482 
483 C...Check if new resonances were produced, loop back if needed.
484  IF(kdcy(jt).NE.3) goto 520
485  np=np+1
486  iref(np,1)=nsd(jt)+1
487  iref(np,2)=nsd(jt)+2
488  iref(np,3)=idoc+1
489  iref(np,4)=idoc+2
490  iref(np,5)=k(iref(ip,jt),2)
491  iref(np,6)=iref(ip,jt)
492  520 CONTINUE
493  530 IF(ip.LT.np) goto 100
494 
495  RETURN
496  END