]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/pythia/pyresd.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New version of MUON from M.Bondila & A.Morsch
[u/mrichter/AliRoot.git] / PYTHIA / pythia / pyresd.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE PYRESD
5
6C...Allows resonances to decay (including parton showers for hadronic
7C...channels).
8 IMPLICIT DOUBLE PRECISION(D)
9 COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
10 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
11 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
12 COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
13 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
14 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
15 COMMON/PYINT1/MINT(400),VINT(400)
16 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
17 COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
18 SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
19 SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT4/
20 DIMENSION IREF(20,8),KDCY(3),KFL1(3),KFL2(3),KEQL(3),NSD(3),
21 &ILIN(6),HGZ(3,3),COUP(6,4),CORL(2,2,2),PK(6,4),PKK(6,6),
22 &CTHE(3),PHI(3),WDTP(0:40),WDTE(0:40,0:5),DBEZQQ(3),DPMO(5)
23 COMPLEX FGK,HA(6,6),HC(6,6)
24
25C...The F, Xi and Xj functions of Gunion and Kunszt
26C...(Phys. Rev. D33, 665, plus errata from the authors).
27 FGK(I1,I2,I3,I4,I5,I6)=4.*HA(I1,I3)*HC(I2,I6)*(HA(I1,I5)*
28 &HC(I1,I4)+HA(I3,I5)*HC(I3,I4))
29 DIGK(DT,DU)=-4.*D34*D56+DT*(3.*DT+4.*DU)+DT**2*(DT*DU/(D34*D56)-
30 &2.*(1./D34+1./D56)*(DT+DU)+2.*(D34/D56+D56/D34))
31 DJGK(DT,DU)=8.*(D34+D56)**2-8.*(D34+D56)*(DT+DU)-6.*DT*DU-
32 &2.*DT*DU*(DT*DU/(D34*D56)-2.*(1./D34+1./D56)*(DT+DU)+
33 &2.*(D34/D56+D56/D34))
34
35C...Some general constants.
36 XW=PARU(102)
37 XWV=XW
38 IF(MSTP(8).GE.2) XW=1.-(PMAS(24,1)/PMAS(23,1))**2
39 XW1=1.-XW
40 SQMZ=PMAS(23,1)**2
41 SQMW=PMAS(24,1)**2
42 SH=VINT(44)
43
44C...Define initial one, two or three objects.
45 ISUB=MINT(1)
46 DO 100 JT=1,8
47 IREF(1,JT)=0
48 100 CONTINUE
49 IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
50 IREF(1,1)=MINT(84)+2+ISET(ISUB)
51 IREF(1,4)=MINT(83)+6+ISET(ISUB)
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,4)=MINT(83)+5+ISET(ISUB)
56 IREF(1,5)=MINT(83)+6+ISET(ISUB)
57 ELSEIF(ISET(ISUB).EQ.5) THEN
58 IREF(1,1)=MINT(84)+3
59 IREF(1,2)=MINT(84)+4
60 IREF(1,3)=MINT(84)+5
61 IREF(1,4)=MINT(83)+7
62 IREF(1,5)=MINT(83)+8
63 IREF(1,6)=MINT(83)+9
64 ELSEIF(ISET(ISUB).EQ.6) THEN
65 IREF(1,1)=MINT(84)+4
66 IREF(1,2)=MINT(84)+5
67 IREF(1,3)=MINT(84)+3
68 IREF(1,4)=MINT(83)+8
69 IREF(1,5)=MINT(83)+9
70 IREF(1,6)=MINT(83)+7
71 ENDIF
72
73C...Check if initial resonance has been moved (in resonance + jet).
74 DO 120 JT=1,3
75 IF(IREF(1,JT).GT.0) THEN
76 IF(K(IREF(1,JT),1).GT.10) THEN
77 KFA=IABS(K(IREF(1,JT),2))
78 IF(KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8.OR.KFA.EQ.39) THEN
79 DO 110 I=IREF(1,JT)+1,N
80 IF(K(I,1).LE.10.AND.K(I,2).EQ.K(IREF(1,JT),2)) IREF(1,JT)=I
81 110 CONTINUE
82 ELSE
83 KDA=MOD(K(IREF(1,JT),4),MSTU(4))
84 IF(KFA.GE.23.AND.KFA.LE.40.AND.KDA.GT.1) IREF(1,JT)=KDA
85 ENDIF
86 ENDIF
87 ENDIF
88 120 CONTINUE
89
90C...Loop over decay history.
91 NP=1
92 IP=0
93 130 IP=IP+1
94 NINH=0
95 JTMAX=2
96 IF(IP.EQ.1.AND.IREF(1,2).EQ.0) JTMAX=1
97 IF(IP.EQ.1.AND.IREF(1,3).NE.0) JTMAX=3
98 ITLH=0
99 NSAV=N
100
101C...Start treatment of one or two resonances in parallel.
102 140 N=NSAV
103 DO 170 JT=1,JTMAX
104 ID=IREF(IP,JT)
105 KDCY(JT)=0
106 KFL1(JT)=0
107 KFL2(JT)=0
108 KEQL(JT)=0
109 NSD(JT)=ID
110 IF(ID.EQ.0) GOTO 160
111 KFA=IABS(K(ID,2))
112 IF((KFA.LT.23.OR.KFA.GT.40).AND.KFA.NE.6.AND.KFA.NE.7.AND.
113 &KFA.NE.8.AND.KFA.NE.17.AND.KFA.NE.18) GOTO 160
114 IF(MSTP(48).LE.0.AND.KFA.EQ.6) GOTO 160
115 IF(MSTP(6).NE.1.AND.MSTP(49).LE.0.AND.(KFA.EQ.7.OR.KFA.EQ.8.OR.
116 &KFA.EQ.17.OR.KFA.EQ.18)) GOTO 160
117 IF(K(ID,1).GT.10.OR.MDCY(KFA,1).EQ.0) GOTO 160
118 IF(KFA.EQ.6.OR.(MSTP(6).NE.1.AND.(KFA.EQ.7.OR.KFA.EQ.8.OR.
119 &KFA.EQ.17.OR.KFA.EQ.18))) ITLH=ITLH+1
120 K(ID,4)=MSTU(5)*(K(ID,4)/MSTU(5))
121 K(ID,5)=MSTU(5)*(K(ID,5)/MSTU(5))
122
123C...Select decay channel.
124 KFB=0
125 IF(ISET(ISUB).NE.6.OR.JT.NE.3) THEN
126 IF(ISUB.EQ.1.OR.ISUB.EQ.15.OR.ISUB.EQ.19.OR.ISUB.EQ.22.OR.
127 & ISUB.EQ.30.OR.ISUB.EQ.35.OR.ISUB.EQ.141) MINT(61)=1
128 CALL PYWIDT(KFA,P(ID,5)**2,WDTP,WDTE)
129 IF(KCHG(KFA,3).EQ.0) THEN
130 IPM=2
131 ELSE
132 IPM=(5-ISIGN(1,K(ID,2)))/2
133 ENDIF
134 IF(JTMAX.GE.2.AND.JT.LE.2) KFB=IABS(K(IREF(IP,3-JT),2))
135 WDTE0S=WDTE(0,1)+WDTE(0,IPM)+WDTE(0,4)
136 IF(KFB.EQ.KFA) WDTE0S=WDTE0S+WDTE(0,5)
137 IF(WDTE0S.LE.0.) THEN
138 IF(KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8.OR.KFA.EQ.17.OR.
139 & KFA.EQ.18) THEN
140 MINT(51)=1
141 RETURN
142 ELSE
143 GOTO 160
144 ENDIF
145 ENDIF
146 RKFL=WDTE0S*RLU(0)
147 IDL=0
148 150 IDL=IDL+1
149 IDC=IDL+MDCY(KFA,2)-1
150 RKFL=RKFL-(WDTE(IDL,1)+WDTE(IDL,IPM)+WDTE(IDL,4))
151 IF(KFB.EQ.KFA) RKFL=RKFL-WDTE(IDL,5)
152 IF(IDL.LT.MDCY(KFA,3).AND.RKFL.GT.0.) GOTO 150
153 ELSE
154 IDC=MINT(35)
155 ENDIF
156
157C...Read out and classify decay channel chosen.
158 KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
159 KFC1A=IABS(KFL1(JT))
160 IF(KFC1A.GT.100) KFC1A=LUCOMP(KFC1A)
161 IF(KCHG(KFC1A,3).EQ.0) KFL1(JT)=IABS(KFL1(JT))
162 KFL2(JT)=KFDP(IDC,2)*ISIGN(1,K(ID,2))
163 KFC2A=IABS(KFL2(JT))
164 IF(KFC2A.GT.100) KFC2A=LUCOMP(KFC2A)
165 IF(KCHG(KFC2A,3).EQ.0) KFL2(JT)=IABS(KFL2(JT))
166 KDCY(JT)=2
167 IF(IABS(KFL1(JT)).LE.10.OR.KFL1(JT).EQ.21) KDCY(JT)=1
168 IF(IABS(KFL1(JT)).GE.23.AND.IABS(KFL1(JT)).LE.40) KDCY(JT)=3
169 IF(KFB.EQ.KFA) KEQL(JT)=MDME(IDC,1)
170 NSD(JT)=N
171 HGZ(JT,1)=VINT(111)
172 HGZ(JT,2)=VINT(112)
173 HGZ(JT,3)=VINT(114)
174
175C...Select masses and check that mass sum not too large.
176 IF(MSTP(42).LE.0.OR.(PMAS(KFC1A,2).LT.PARP(41).AND.
177 &PMAS(KFC2A,2).LT.PARP(41))) THEN
178 P(N+1,5)=PMAS(KFC1A,1)
179 P(N+2,5)=PMAS(KFC2A,1)
180 IF(P(N+1,5)+P(N+2,5)+PARJ(64).GT.P(ID,5)) THEN
181 CALL LUERRM(13,'(PYRESD:) daughter masses too large')
182 MINT(51)=1
183 RETURN
184 ENDIF
185 ELSEIF(IP.EQ.1) THEN
186 CALL PYOFSH(2,KFA,KFL1(JT),KFL2(JT),P(ID,5),P(N+1,5),P(N+2,5))
187 IF(MINT(51).EQ.1) RETURN
188 ELSE
189 CALL PYOFSH(7,KFA,KFL1(JT),KFL2(JT),P(ID,5),P(N+1,5),P(N+2,5))
190 IF(MINT(51).EQ.1) RETURN
191 ENDIF
192
193C...Fill decay products, prepared for parton showers for quarks.
194C...Special cases, done by hand, for techni-eta, t, leptoquark and q*.
195 MSTU(10)=1
196 IF(KFA.EQ.38.OR.KFA.EQ.39.OR.((MSTP(6).EQ.1.OR.MSTP(49).GE.1).AND.
197 &(KFA.EQ.7.OR.KFA.EQ.8)).OR.KFA.EQ.6) THEN
198 MSTU(19)=1
199 CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
200 ISID=4
201 IF(K(ID,2).LT.0) ISID=5
202 IF(KFA.EQ.38) THEN
203 IF(KFC1A.EQ.21.AND.RLU(0).GT.0.5) ISID=9-ISID
204 K(N-1,1)=3
205 K(N,1)=3
206 K(ID,ISID)=K(ID,ISID)+(N-1)
207 K(ID,9-ISID)=K(ID,9-ISID)+N
208 K(N-1,ISID)=MSTU(5)*ID
209 K(N-1,9-ISID)=MSTU(5)*N
210 K(N,ISID)=MSTU(5)*(N-1)
211 K(N,9-ISID)=MSTU(5)*ID
212 ELSEIF(KFA.EQ.6.OR.(MSTP(6).NE.1.AND.(KFA.EQ.7.OR.KFA.EQ.8)))
213 & THEN
214 K(N-1,1)=1
215 K(N,1)=3
216 K(ID,ISID)=K(ID,ISID)+N
217 K(N,ISID)=MSTU(5)*ID
218 ELSEIF(KFA.EQ.39) THEN
219 K(N-1,1)=3
220 K(N,1)=1
221 K(ID,ISID)=K(ID,ISID)+(N-1)
222 K(N-1,ISID)=MSTU(5)*ID
223 ELSEIF(KFL1(JT).NE.21) THEN
224 K(N-1,1)=1
225 K(N,1)=3
226 K(ID,ISID)=K(ID,ISID)+N
227 K(N,ISID)=MSTU(5)*ID
228 ELSE
229 K(N-1,1)=3
230 K(N,1)=3
231 K(ID,ISID)=K(ID,ISID)+(N-1)
232 K(N-1,ISID)=MSTU(5)*ID
233 K(N-1,9-ISID)=MSTU(5)*N
234 K(N,ISID)=MSTU(5)*(N-1)
235 ENDIF
236 ELSEIF(KDCY(JT).EQ.1) THEN
237 CALL LU2ENT(-(N+1),KFL1(JT),KFL2(JT),P(ID,5))
238 ELSE
239 CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
240 ENDIF
241 MSTU(10)=2
242 160 IF(KFA.GE.23.AND.KFA.LE.40.AND.KFL1(JT).EQ.0) NINH=NINH+1
243 170 CONTINUE
244
245C...Check for allowed combinations. Skip if no decays.
246 IF(JTMAX.GE.2) THEN
247 IF(KEQL(1).EQ.4.AND.KEQL(2).EQ.4) GOTO 140
248 IF(KEQL(1).EQ.5.AND.KEQL(2).EQ.5) GOTO 140
249 ENDIF
250 IF(JTMAX.EQ.1.AND.KDCY(1).EQ.0) GOTO 480
251 IF(JTMAX.EQ.2.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0) GOTO 480
252 IF(JTMAX.EQ.3.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0.AND.
253 &KDCY(3).EQ.0) GOTO 480
254
255C...Order incoming partons and outgoing resonances.
256 IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
257 ILIN(1)=MINT(84)+1
258 IF(K(MINT(84)+1,2).GT.0) ILIN(1)=MINT(84)+2
259 IF(K(ILIN(1),2).EQ.21) ILIN(1)=2*MINT(84)+3-ILIN(1)
260 ILIN(2)=2*MINT(84)+3-ILIN(1)
261 IMIN=1
262 IF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.IREF(IP,7)
263 & .EQ.36) IMIN=3
264 IMAX=2
265 IORD=1
266 IF(K(IREF(IP,1),2).EQ.23) IORD=2
267 IF(K(IREF(IP,1),2).EQ.24.AND.K(IREF(IP,2),2).EQ.-24) IORD=2
268 IAKIPD=IABS(K(IREF(IP,IORD),2))
269 IF(IAKIPD.EQ.25.OR.IAKIPD.EQ.35.OR.IAKIPD.EQ.36) IORD=3-IORD
270 IF(KDCY(IORD).EQ.0) IORD=3-IORD
271
272C...Order decay products of resonances.
273 DO 180 JT=IORD,3-IORD,3-2*IORD
274 IF(KDCY(JT).EQ.0) THEN
275 ILIN(IMAX+1)=NSD(JT)
276 IMAX=IMAX+1
277 ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
278 ILIN(IMAX+1)=N+2*JT-1
279 ILIN(IMAX+2)=N+2*JT
280 IMAX=IMAX+2
281 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
282 K(N+2*JT,2)=K(NSD(JT)+2,2)
283 ELSE
284 ILIN(IMAX+1)=N+2*JT
285 ILIN(IMAX+2)=N+2*JT-1
286 IMAX=IMAX+2
287 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
288 K(N+2*JT,2)=K(NSD(JT)+2,2)
289 ENDIF
290 180 CONTINUE
291
292C...Find charge, isospin, left- and righthanded couplings.
293 DO 200 I=IMIN,IMAX
294 DO 190 J=1,4
295 COUP(I,J)=0.
296 190 CONTINUE
297 KFA=IABS(K(ILIN(I),2))
298 IF(KFA.EQ.0.OR.KFA.GT.20) GOTO 200
299 COUP(I,1)=KCHG(KFA,1)/3.
300 COUP(I,2)=(-1)**MOD(KFA,2)
301 COUP(I,4)=-2.*COUP(I,1)*XWV
302 COUP(I,3)=COUP(I,2)+COUP(I,4)
303 200 CONTINUE
304
305C...Full propagator dependence and flavour correlations for 2 gamma*/Z.
306 IF(ISUB.EQ.22) THEN
307 DO 230 I=3,5,2
308 I1=IORD
309 IF(I.EQ.5) I1=3-IORD
310 DO 220 J1=1,2
311 DO 210 J2=1,2
312 CORL(I/2,J1,J2)=COUP(1,1)**2*HGZ(I1,1)*COUP(I,1)**2/16.+
313 & COUP(1,1)*COUP(1,J1+2)*HGZ(I1,2)*COUP(I,1)*COUP(I,J2+2)/4.+
314 & COUP(1,J1+2)**2*HGZ(I1,3)*COUP(I,J2+2)**2
315 210 CONTINUE
316 220 CONTINUE
317 230 CONTINUE
318 COWT12=(CORL(1,1,1)+CORL(1,1,2))*(CORL(2,1,1)+CORL(2,1,2))+
319 & (CORL(1,2,1)+CORL(1,2,2))*(CORL(2,2,1)+CORL(2,2,2))
320 COMX12=(CORL(1,1,1)+CORL(1,1,2)+CORL(1,2,1)+CORL(1,2,2))*
321 & (CORL(2,1,1)+CORL(2,1,2)+CORL(2,2,1)+CORL(2,2,2))
322 IF(COWT12.LT.RLU(0)*COMX12) GOTO 140
323 ENDIF
324 ENDIF
325
326C...Select angular orientation type - Z'/W' only.
327 MZPWP=0
328 IF(ISUB.EQ.141) THEN
329 IF(RLU(0).LT.PARU(130)) MZPWP=1
330 IF(IP.EQ.2) THEN
331 IF(IABS(K(IREF(2,1),2)).EQ.37) MZPWP=2
332 IAKIR=IABS(K(IREF(2,2),2))
333 IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
334 ENDIF
335 IF(IP.GE.3) MZPWP=2
336 ELSEIF(ISUB.EQ.142) THEN
337 IF(RLU(0).LT.PARU(136)) MZPWP=1
338 IF(IP.EQ.2) THEN
339 IAKIR=IABS(K(IREF(2,2),2))
340 IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
341 ENDIF
342 IF(IP.GE.3) MZPWP=2
343 ENDIF
344
345C...Select random angles (begin of weighting procedure).
346 240 DO 250 JT=1,JTMAX
347 IF(KDCY(JT).EQ.0) GOTO 250
348 IF(ISET(ISUB).EQ.6.AND.JT.EQ.3) GOTO 250
349 IF(JTMAX.EQ.1) THEN
350 CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))*RLU(0)
351 IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
352 PHI(JT)=VINT(24)
353 ELSE
354 CTHE(JT)=2.*RLU(0)-1.
355 PHI(JT)=PARU(2)*RLU(0)
356 ENDIF
357 250 CONTINUE
358
359 IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
360C...Construct massless four-vectors.
361 DO 270 I=N+1,N+4
362 K(I,1)=1
363 DO 260 J=1,5
364 P(I,J)=0.
365 V(I,J)=0.
366 260 CONTINUE
367 270 CONTINUE
368 DO 280 JT=1,JTMAX
369 IF(KDCY(JT).EQ.0) GOTO 280
370 ID=IREF(IP,JT)
371 P(N+2*JT-1,3)=0.5*P(ID,5)
372 P(N+2*JT-1,4)=0.5*P(ID,5)
373 P(N+2*JT,3)=-0.5*P(ID,5)
374 P(N+2*JT,4)=0.5*P(ID,5)
375 CALL LUDBRB(N+2*JT-1,N+2*JT,ACOS(CTHE(JT)),PHI(JT),DBLE(P(ID,1)/
376 & P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
377 280 CONTINUE
378
379C...Store incoming and outgoing momenta, with random rotation to
380C...avoid accidental zeroes in HA expressions.
381 DO 300 I=1,IMAX
382 K(N+4+I,1)=1
383 P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
384 & P(ILIN(I),5)**2)
385 P(N+4+I,5)=P(ILIN(I),5)
386 DO 290 J=1,3
387 P(N+4+I,J)=P(ILIN(I),J)
388 290 CONTINUE
389 300 CONTINUE
390 310 THERR=ACOS(2.*RLU(0)-1.)
391 PHIRR=PARU(2)*RLU(0)
392 CALL LUDBRB(N+5,N+4+IMAX,THERR,PHIRR,0D0,0D0,0D0)
393 DO 330 I=1,IMAX
394 IF(P(N+4+I,1)**2+P(N+4+I,2)**2.LT.1E-4*P(N+4+I,4)**2) GOTO 310
395 DO 320 J=1,4
396 PK(I,J)=P(N+4+I,J)
397 320 CONTINUE
398 330 CONTINUE
399
400C...Calculate internal products.
401 IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25.OR.ISUB.EQ.141.OR.
402 & ISUB.EQ.142) THEN
403 DO 350 I1=IMIN,IMAX-1
404 DO 340 I2=I1+1,IMAX
405 HA(I1,I2)=SQRT((PK(I1,4)-PK(I1,3))*(PK(I2,4)+PK(I2,3))/
406 & (1E-20+PK(I1,1)**2+PK(I1,2)**2))*CMPLX(PK(I1,1),PK(I1,2))-
407 & SQRT((PK(I1,4)+PK(I1,3))*(PK(I2,4)-PK(I2,3))/
408 & (1E-20+PK(I2,1)**2+PK(I2,2)**2))*CMPLX(PK(I2,1),PK(I2,2))
409 HC(I1,I2)=CONJG(HA(I1,I2))
410 IF(I1.LE.2) HA(I1,I2)=CMPLX(0.,1.)*HA(I1,I2)
411 IF(I1.LE.2) HC(I1,I2)=CMPLX(0.,1.)*HC(I1,I2)
412 HA(I2,I1)=-HA(I1,I2)
413 HC(I2,I1)=-HC(I1,I2)
414 340 CONTINUE
415 350 CONTINUE
416 ENDIF
417 DO 370 I=1,2
418 DO 360 J=1,4
419 PK(I,J)=-PK(I,J)
420 360 CONTINUE
421 370 CONTINUE
422 DO 390 I1=IMIN,IMAX-1
423 DO 380 I2=I1+1,IMAX
424 PKK(I1,I2)=2.*(PK(I1,4)*PK(I2,4)-PK(I1,1)*PK(I2,1)-
425 & PK(I1,2)*PK(I2,2)-PK(I1,3)*PK(I2,3))
426 PKK(I2,I1)=PKK(I1,I2)
427 380 CONTINUE
428 390 CONTINUE
429 ENDIF
430
431 KFAGM=IABS(IREF(IP,7))
432 IF(MSTP(47).LE.0.OR.NINH.NE.0) THEN
433C...Isotropic decay selected by user.
434 WT=1.
435 WTMAX=1.
436
437 ELSEIF(ITLH.GE.1) THEN
438C... Isotropic decay t -> b + W etc for 4th generation q and l.
439 WT=1.
440 WTMAX=1.
441
442 ELSEIF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.
443 &IREF(IP,7).EQ.36) THEN
444C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons.
445 WT=16.*PKK(3,5)*PKK(4,6)
446 IF(IP.EQ.1) WTMAX=SH**2
447 IF(IP.GE.2) WTMAX=P(IREF(IP,8),5)**4
448 KFA=IABS(K(IREF(IP,1),2))
449 IF(KFA.NE.23.AND.KFA.NE.24) WT=WTMAX
450
451 ELSEIF((KFAGM.EQ.6.OR.(MSTP(6).NE.1.AND.(KFAGM.EQ.7.OR.
452 &KFAGM.EQ.8.OR.KFAGM.EQ.17.OR.KFAGM.EQ.18))).AND.
453 &IABS(K(IREF(IP,1),2)).EQ.24) THEN
454C...Angular correlation in f -> f' + W -> f' + 2 quarks/leptons.
455 I1=IREF(IP,8)
456 IF(MOD(KFAGM,2).EQ.0) THEN
457 I2=N+1
458 I3=N+2
459 ELSE
460 I2=N+2
461 I3=N+1
462 ENDIF
463 I4=IREF(IP,2)
464 WT=(P(I1,4)*P(I2,4)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
465 & P(I1,3)*P(I2,3))*(P(I3,4)*P(I4,4)-P(I3,1)*P(I4,1)-
466 & P(I3,2)*P(I4,2)-P(I3,3)*P(I4,3))
467 WTMAX=(P(I1,5)**4-P(IREF(IP,1),5)**4)/8.
468 IF(KFAGM.EQ.6.AND.MSTP(48).LE.1) WT=WTMAX
469 IF(KFAGM.NE.6.AND.MSTP(49).LE.1) WT=WTMAX
470
471 ELSEIF(ISUB.EQ.1) THEN
472C...Angular weight for gamma*/Z0 -> 2 quarks/leptons.
473 EI=KCHG(IABS(MINT(15)),1)/3.
474 AI=SIGN(1.,EI+0.1)
475 VI=AI-4.*EI*XWV
476 EF=KCHG(IABS(KFL1(1)),1)/3.
477 AF=SIGN(1.,EF+0.1)
478 VF=AF-4.*EF*XWV
479 ASYM=2.*(EI*AI*VINT(112)*EF*AF+4.*VI*AI*VINT(114)*VF*AF)/
480 & (EI**2*VINT(111)*EF**2+EI*VI*VINT(112)*EF*VF+
481 & (VI**2+AI**2)*VINT(114)*(VF**2+AF**2))
482 WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
483 WTMAX=2.+ABS(ASYM)
484
485 ELSEIF(ISUB.EQ.2) THEN
486C...Angular weight for W+/- -> 2 quarks/leptons.
487 WT=(1.+CTHE(1)*ISIGN(1,MINT(15)*KFL1(1)))**2
488 WTMAX=4.
489
490 ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
491C...Angular weight for f + f~ -> gluon/gamma + (gamma*/Z0) ->
492C...-> gluon/gamma + 2 quarks/leptons.
493 CLILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
494 & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
495 & COUP(1,3)**2*HGZ(2,3)*COUP(3,3)**2
496 CLIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
497 & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
498 & COUP(1,3)**2*HGZ(2,3)*COUP(3,4)**2
499 CRILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
500 & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
501 & COUP(1,4)**2*HGZ(2,3)*COUP(3,3)**2
502 CRIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
503 & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
504 & COUP(1,4)**2*HGZ(2,3)*COUP(3,4)**2
505 WT=(CLILF+CRIRF)*(PKK(1,3)**2+PKK(2,4)**2)+
506 & (CLIRF+CRILF)*(PKK(1,4)**2+PKK(2,3)**2)
507 WTMAX=(CLILF+CLIRF+CRILF+CRIRF)*
508 & ((PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2)
509
510 ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
511C...Angular weight for f + f~' -> gluon/gamma + W+/- ->
512C...-> gluon/gamma + 2 quarks/leptons.
513 WT=PKK(1,3)**2+PKK(2,4)**2
514 WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2
515
516 ELSEIF(ISUB.EQ.22) THEN
517C...Angular weight for f + f~ -> Z0 + Z0 -> 4 quarks/leptons.
518 S34=P(IREF(IP,IORD),5)**2
519 S56=P(IREF(IP,3-IORD),5)**2
520 TI=PKK(1,3)+PKK(1,4)+S34
521 UI=PKK(1,5)+PKK(1,6)+S56
522 FGK135=ABS(FGK(1,2,3,4,5,6)/TI+FGK(1,2,5,6,3,4)/UI)**2
523 FGK145=ABS(FGK(1,2,4,3,5,6)/TI+FGK(1,2,5,6,4,3)/UI)**2
524 FGK136=ABS(FGK(1,2,3,4,6,5)/TI+FGK(1,2,6,5,3,4)/UI)**2
525 FGK146=ABS(FGK(1,2,4,3,6,5)/TI+FGK(1,2,6,5,4,3)/UI)**2
526 FGK253=ABS(FGK(2,1,5,6,3,4)/TI+FGK(2,1,3,4,5,6)/UI)**2
527 FGK263=ABS(FGK(2,1,6,5,3,4)/TI+FGK(2,1,3,4,6,5)/UI)**2
528 FGK254=ABS(FGK(2,1,5,6,4,3)/TI+FGK(2,1,4,3,5,6)/UI)**2
529 FGK264=ABS(FGK(2,1,6,5,4,3)/TI+FGK(2,1,4,3,6,5)/UI)**2
530 WT=
531 & CORL(1,1,1)*CORL(2,1,1)*FGK135+CORL(1,1,2)*CORL(2,1,1)*FGK145+
532 & CORL(1,1,1)*CORL(2,1,2)*FGK136+CORL(1,1,2)*CORL(2,1,2)*FGK146+
533 & CORL(1,2,1)*CORL(2,2,1)*FGK253+CORL(1,2,2)*CORL(2,2,1)*FGK263+
534 & CORL(1,2,1)*CORL(2,2,2)*FGK254+CORL(1,2,2)*CORL(2,2,2)*FGK264
535 WTMAX=16.*((CORL(1,1,1)+CORL(1,1,2))*(CORL(2,1,1)+CORL(2,1,2))+
536 & (CORL(1,2,1)+CORL(1,2,2))*(CORL(2,2,1)+CORL(2,2,2)))*S34*S56*
537 & ((TI**2+UI**2+2.*SH*(S34+S56))/(TI*UI)-S34*S56*(1./TI**2+
538 & 1./UI**2))
539
540 ELSEIF(ISUB.EQ.23) THEN
541C...Angular weight for f + f~' -> Z0 + W+/- -> 4 quarks/leptons.
542 D34=P(IREF(IP,IORD),5)**2
543 D56=P(IREF(IP,3-IORD),5)**2
544 DT=PKK(1,3)+PKK(1,4)+D34
545 DU=PKK(1,5)+PKK(1,6)+D56
546 FACBW=1./((SH-SQMW)**2+SQMW*PMAS(24,2)**2)
547 CAWZ=COUP(2,3)/SNGL(DT)-2.*XW1*COUP(1,2)*(SH-SQMW)*FACBW
548 CBWZ=COUP(1,3)/SNGL(DU)+2.*XW1*COUP(1,2)*(SH-SQMW)*FACBW
549 FGK135=ABS(CAWZ*FGK(1,2,3,4,5,6)+CBWZ*FGK(1,2,5,6,3,4))
550 FGK136=ABS(CAWZ*FGK(1,2,3,4,6,5)+CBWZ*FGK(1,2,6,5,3,4))
551 WT=(COUP(5,3)*FGK135)**2+(COUP(5,4)*FGK136)**2
552 WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*(CAWZ**2*
553 & DIGK(DT,DU)+CBWZ**2*DIGK(DU,DT)+CAWZ*CBWZ*DJGK(DT,DU))
554
555 ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.171.OR.ISUB.EQ.176) THEN
556C...Angular weight for f + f~ -> Z0 + H0 -> 2 quarks/leptons + H0
557C...(or H'0, or A0).
558 WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
559 & PKK(1,3)*PKK(2,4)+((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*
560 & COUP(3,3))**2)*PKK(1,4)*PKK(2,3)
561 WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
562 & (PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
563
564 ELSEIF(ISUB.EQ.25) THEN
565C...Angular weight for f + f~ -> W+ + W- -> 4 quarks/leptons.
566 D34=P(IREF(IP,IORD),5)**2
567 D56=P(IREF(IP,3-IORD),5)**2
568 DT=PKK(1,3)+PKK(1,4)+D34
569 DU=PKK(1,5)+PKK(1,6)+D56
570 FACBW=1./((SH-SQMZ)**2+SQMZ*PMAS(23,2)**2)
571 CDWW=(COUP(1,3)*SQMZ*(SH-SQMZ)*FACBW+COUP(1,2))/SH
572 CAWW=CDWW+0.5*(COUP(1,2)+1.)/SNGL(DT)
573 CBWW=CDWW+0.5*(COUP(1,2)-1.)/SNGL(DU)
574 CCWW=COUP(1,4)*SQMZ*(SH-SQMZ)*FACBW/SH
575 FGK135=ABS(CAWW*FGK(1,2,3,4,5,6)-CBWW*FGK(1,2,5,6,3,4))
576 FGK253=ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))
577 WT=FGK135**2+(CCWW*FGK253)**2
578 WTMAX=4.*D34*D56*(CAWW**2*DIGK(DT,DU)+CBWW**2*DIGK(DU,DT)-CAWW*
579 & CBWW*DJGK(DT,DU)+CCWW**2*(DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU)))
580
581 ELSEIF(ISUB.EQ.26.OR.ISUB.EQ.172.OR.ISUB.EQ.177) THEN
582C...Angular weight for f + f~' -> W+/- + H0 -> 2 quarks/leptons + H0
583C...(or H'0, or A0).
584 WT=PKK(1,3)*PKK(2,4)
585 WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
586
587 ELSEIF(ISUB.EQ.30.OR.ISUB.EQ.35) THEN
588C...Angular weight for f + g/gamma -> f + (gamma*/Z0)
589C...-> f + 2 quarks/leptons.
590 CLILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
591 & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
592 & COUP(1,3)**2*HGZ(2,3)*COUP(3,3)**2
593 CLIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
594 & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
595 & COUP(1,3)**2*HGZ(2,3)*COUP(3,4)**2
596 CRILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
597 & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
598 & COUP(1,4)**2*HGZ(2,3)*COUP(3,3)**2
599 CRIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
600 & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
601 & COUP(1,4)**2*HGZ(2,3)*COUP(3,4)**2
602 IF(K(ILIN(1),2).GT.0) WT=(CLILF+CRIRF)*(PKK(1,4)**2+
603 & PKK(3,5)**2)+(CLIRF+CRILF)*(PKK(1,3)**2+PKK(4,5)**2)
604 IF(K(ILIN(1),2).LT.0) WT=(CLILF+CRIRF)*(PKK(1,3)**2+
605 & PKK(4,5)**2)+(CLIRF+CRILF)*(PKK(1,4)**2+PKK(3,5)**2)
606 WTMAX=(CLILF+CLIRF+CRILF+CRIRF)*
607 & ((PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2)
608
609 ELSEIF(ISUB.EQ.31) THEN
610C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons.
611 IF(K(ILIN(1),2).GT.0) WT=PKK(1,4)**2+PKK(3,5)**2
612 IF(K(ILIN(1),2).LT.0) WT=PKK(1,3)**2+PKK(4,5)**2
613 WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2
614
615 ELSEIF(ISUB.EQ.71.OR.ISUB.EQ.72.OR.ISUB.EQ.73.OR.ISUB.EQ.76.OR.
616 &ISUB.EQ.77) THEN
617C...Angular weight for V_L1 + V_L2 -> V_L3 + V_L4 (V = Z/W).
618 WT=16.*PKK(3,5)*PKK(4,6)
619 WTMAX=SH**2
620
621 ELSEIF(ISUB.EQ.110) THEN
622C...Angular weight for f + f~ -> gamma + H0 -> gamma + X is isotropic.
623 WT=1.
624 WTMAX=1.
625
626 ELSEIF(ISUB.EQ.141) THEN
627 IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
628C...Angular weight for f + f~ -> gamma*/Z0/Z'0 -> 2 quarks/leptons.
629C...Couplings of incoming flavour.
630 KFAI=IABS(MINT(15))
631 EI=KCHG(KFAI,1)/3.
632 AI=SIGN(1.,EI+0.1)
633 VI=AI-4.*EI*XWV
634 KFAIC=1
635 IF(KFAI.LE.10.AND.MOD(KFAI,2).EQ.0) KFAIC=2
636 IF(KFAI.GT.10.AND.MOD(KFAI,2).NE.0) KFAIC=3
637 IF(KFAI.GT.10.AND.MOD(KFAI,2).EQ.0) KFAIC=4
638 VPI=PARU(119+2*KFAIC)
639 API=PARU(120+2*KFAIC)
640C...Couplings of final flavour.
641 KFAF=IABS(KFL1(1))
642 EF=KCHG(KFAF,1)/3.
643 AF=SIGN(1.,EF+0.1)
644 VF=AF-4.*EF*XWV
645 KFAFC=1
646 IF(KFAF.LE.10.AND.MOD(KFAF,2).EQ.0) KFAFC=2
647 IF(KFAF.GT.10.AND.MOD(KFAF,2).NE.0) KFAFC=3
648 IF(KFAF.GT.10.AND.MOD(KFAF,2).EQ.0) KFAFC=4
649 VPF=PARU(119+2*KFAFC)
650 APF=PARU(120+2*KFAFC)
651C...Asymmetry and weight.
652 ASYM=2.*(EI*AI*VINT(112)*EF*AF+EI*API*VINT(113)*EF*APF+
653 & 4.*VI*AI*VINT(114)*VF*AF+(VI*API+VPI*AI)*VINT(115)*
654 & (VF*APF+VPF*AF)+4.*VPI*API*VINT(116)*VPF*APF)/
655 & (EI**2*VINT(111)*EF**2+EI*VI*VINT(112)*EF*VF+
656 & EI*VPI*VINT(113)*EF*VPF+(VI**2+AI**2)*VINT(114)*
657 & (VF**2+AF**2)+(VI*VPI+AI*API)*VINT(115)*(VF*VPF+AF*APF)+
658 & (VPI**2+API**2)*VINT(116)*(VPF**2+APF**2))
659 WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
660 WTMAX=2.+ABS(ASYM)
661 ELSEIF(IP.EQ.1.AND.IABS(KFL1(1)).EQ.24) THEN
662C...Angular weight for f + f~ -> Z' -> W+ + W-.
663 RM1=P(NSD(1)+1,5)**2/SH
664 RM2=P(NSD(1)+2,5)**2/SH
665 CCOS2=-(1./16.)*((1.-RM1-RM2)**2-4.*RM1*RM2)*
666 & (1.-2.*RM1-2.*RM2+RM1**2+RM2**2+10.*RM1*RM2)
667 CFLAT=-CCOS2+0.5*(RM1+RM2)*(1.-2.*RM1-2.*RM2+(RM2-RM1)**2)
668 WT=CFLAT+CCOS2*CTHE(1)**2
669 WTMAX=CFLAT+MAX(0.,CCOS2)
670 ELSEIF(IP.EQ.1.AND.(KFL1(1).EQ.25.OR.KFL1(1).EQ.35.OR.
671 & IABS(KFL1(1)).EQ.37)) THEN
672C...Angular weight for f + f~ -> Z' -> H0 + A0, H'0 + A0, H+ + H-.
673 WT=1.-CTHE(1)**2
674 WTMAX=1.
675 ELSEIF(IP.EQ.1.AND.KDCY(1).EQ.3) THEN
676C...Angular weight for f + f~ -> Z' -> Z0 + H0.
677 RM1=P(NSD(1)+1,5)**2/SH
678 RM2=P(NSD(1)+2,5)**2/SH
679 FLAM2=MAX(0.,(1.-RM1-RM2)**2-4.*RM1*RM2)
680 WT=1.+FLAM2*(1.-CTHE(1)**2)/(8.*RM1)
681 WTMAX=1.+FLAM2/(8.*RM1)
682 ELSEIF(MZPWP.EQ.0) THEN
683C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
684C...(W:s like if intermediate Z).
685 D34=P(IREF(IP,IORD),5)**2
686 D56=P(IREF(IP,3-IORD),5)**2
687 DT=PKK(1,3)+PKK(1,4)+D34
688 DU=PKK(1,5)+PKK(1,6)+D56
689 FGK135=ABS(FGK(1,2,3,4,5,6)-FGK(1,2,5,6,3,4))
690 FGK253=ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))
691 WT=(COUP(1,3)*FGK135)**2+(COUP(1,4)*FGK253)**2
692 WTMAX=4.*D34*D56*(COUP(1,3)**2+COUP(1,4)**2)*
693 & (DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU))
694 ELSEIF(MZPWP.EQ.1) THEN
695C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
696C...(W:s approximately longitudinal, like if intermediate H).
697 WT=16.*PKK(3,5)*PKK(4,6)
698 WTMAX=SH**2
699 ELSE
700C...Angular weight for f + f~ -> Z' -> H+ + H-, Z0 + H0, H0 + A0,
701C...H'0 + A0 -> 4 quarks/leptons.
702 WT=1.
703 WTMAX=1.
704 ENDIF
705
706 ELSEIF(ISUB.EQ.142) THEN
707 IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
708C...Angular weight for f + f~' -> W'+/- -> 2 quarks/leptons.
709 KFAI=IABS(MINT(15))
710 KFAIC=1
711 IF(KFAI.GT.10) KFAIC=2
712 VI=PARU(129+2*KFAIC)
713 AI=PARU(130+2*KFAIC)
714 KFAF=IABS(KFL1(1))
715 KFAFC=1
716 IF(KFAF.GT.10) KFAFC=2
717 VF=PARU(129+2*KFAFC)
718 AF=PARU(130+2*KFAFC)
719 ASYM=8.*VI*AI*VF*AF/((VI**2+AI**2)*(VF**2+AF**2))
720 WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
721 WTMAX=2.+ABS(ASYM)
722 ELSEIF(IP.EQ.1.AND.IABS(KFL2(1)).EQ.23) THEN
723C...Angular weight for f + f~' -> W'+/- -> W+/- + Z0.
724 RM1=P(NSD(1)+1,5)**2/SH
725 RM2=P(NSD(1)+2,5)**2/SH
726 CCOS2=-(1./16.)*((1.-RM1-RM2)**2-4.*RM1*RM2)*
727 & (1.-2.*RM1-2.*RM2+RM1**2+RM2**2+10.*RM1*RM2)
728 CFLAT=-CCOS2+0.5*(RM1+RM2)*(1.-2.*RM1-2.*RM2+(RM2-RM1)**2)
729 WT=CFLAT+CCOS2*CTHE(1)**2
730 WTMAX=CFLAT+MAX(0.,CCOS2)
731 ELSEIF(IP.EQ.1.AND.KDCY(1).EQ.3) THEN
732C...Angular weight for f + f~ -> W'+/- -> W+/- + H0.
733 RM1=P(NSD(1)+1,5)**2/SH
734 RM2=P(NSD(1)+2,5)**2/SH
735 FLAM2=MAX(0.,(1.-RM1-RM2)**2-4.*RM1*RM2)
736 WT=1.+FLAM2*(1.-CTHE(1)**2)/(8.*RM1)
737 WTMAX=1.+FLAM2/(8.*RM1)
738 ELSEIF(MZPWP.EQ.0) THEN
739C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
740C...(W/Z like if intermediate W).
741 D34=P(IREF(IP,IORD),5)**2
742 D56=P(IREF(IP,3-IORD),5)**2
743 DT=PKK(1,3)+PKK(1,4)+D34
744 DU=PKK(1,5)+PKK(1,6)+D56
745 FGK135=ABS(FGK(1,2,3,4,5,6)-FGK(1,2,5,6,3,4))
746 FGK136=ABS(FGK(1,2,3,4,6,5)-FGK(1,2,6,5,3,4))
747 WT=(COUP(5,3)*FGK135)**2+(COUP(5,4)*FGK136)**2
748 WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*
749 & (DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU))
750 ELSEIF(MZPWP.EQ.1) THEN
751C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
752C...(W/Z approximately longitudinal, like if intermediate H).
753 WT=16.*PKK(3,5)*PKK(4,6)
754 WTMAX=SH**2
755 ELSE
756C...Angular weight for f + f~ -> W' -> W + H0 -> whatever.
757 WT=1.
758 WTMAX=1.
759 ENDIF
760
761 ELSEIF(ISUB.EQ.145.OR.ISUB.EQ.162.OR.ISUB.EQ.163.OR.ISUB.EQ.164)
762 &THEN
763C...Isotropic decay of leptoquarks (assumed spin 0).
764 WT=1.
765 WTMAX=1.
766
767 ELSEIF(ISUB.EQ.147.OR.ISUB.EQ.148) THEN
768C...Decays of (spin 1/2) q* -> q + (g,gamma) or (Z0,W+-).
769 SIDE=1.
770 IF(MINT(16).EQ.21) SIDE=-1.
771 IF(IP.EQ.1.AND.(KFL1(1).EQ.21.OR.KFL1(1).EQ.22)) THEN
772 WT=1.+SIDE*CTHE(1)
773 WTMAX=2.
774 ELSEIF(IP.EQ.1) THEN
775 RM1=P(NSD(1)+1,5)**2/SH
776 WT=1.+SIDE*CTHE(1)*(1.-0.5*RM1)/(1.+0.5*RM1)
777 WTMAX=1.+(1.-0.5*RM1)/(1.+0.5*RM1)
778 ELSE
779C...W/Z decay assumed isotropic, since not known.
780 WT=1.
781 WTMAX=1.
782 ENDIF
783
784 ELSEIF(ISUB.EQ.149) THEN
785C...Isotropic decay of techni-eta.
786 WT=1.
787 WTMAX=1.
788
789C...Obtain correct angular distribution by rejection techniques.
790 ELSE
791 WT=1.
792 WTMAX=1.
793 ENDIF
794 IF(WT.LT.RLU(0)*WTMAX) GOTO 240
795
796C...Construct massive four-vectors using angles chosen. Mark decayed
797C...resonances, add documentation lines. Shower evolution.
798 400 DO 470 JT=1,JTMAX
799 IF(KDCY(JT).EQ.0) GOTO 470
800 ID=IREF(IP,JT)
801 IF(ISET(ISUB).NE.6.OR.JT.NE.3) THEN
802 DO 410 J=1,5
803 DPMO(J)=P(ID,J)
804 410 CONTINUE
805 DPMO(4)=SQRT(DPMO(1)**2+DPMO(2)**2+DPMO(3)**2+DPMO(5)**2)
806 CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,ACOS(CTHE(JT)),PHI(JT),
807 & DPMO(1)/DPMO(4),DPMO(2)/DPMO(4),DPMO(3)/DPMO(4))
808 ELSE
809C...Z + q + q~ : angles already known, in rest frame of system.
810 DO 420 J=1,3
811 DBEZQQ(J)=(P(ID,J)+P(ID+1,J)+P(ID+2,J))/(P(ID,4)+P(ID+1,4)+
812 & P(ID+2,4))
813 420 CONTINUE
814 K(N+1,1)=1
815 DO 430 J=1,5
816 P(N+1,J)=P(ID,J)
817 430 CONTINUE
818 CALL LUDBRB(N+1,N+1,0.,0.,-DBEZQQ(1),-DBEZQQ(2),-DBEZQQ(3))
819 PHIZQQ=ULANGL(P(N+1,1),P(N+1,2))
820 THEZQQ=ULANGL(P(N+1,3),SQRT(P(N+1,1)**2+P(N+1,2)**2))
821 CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,ACOS(VINT(81)),VINT(82),
822 & 0D0,0D0,DBLE(SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)/
823 & P(N+1,4)))
824 CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,THEZQQ,PHIZQQ,DBEZQQ(1),
825 & DBEZQQ(2),DBEZQQ(3))
826 ENDIF
827 K(ID,1)=K(ID,1)+10
828 KFA=IABS(K(ID,2))
829 IF(KFA.EQ.38.OR.KFA.EQ.39.OR.((MSTP(6).EQ.1.OR.MSTP(49).GE.1).AND.
830 &(KFA.EQ.7.OR.KFA.EQ.8)).OR.(MSTP(48).GE.1.AND.KFA.EQ.6)) THEN
831C...Do not kill colour flow through techni-eta, t, leptoquark or q*!
832 ELSE
833 K(ID,4)=NSD(JT)+1
834 K(ID,5)=NSD(JT)+2
835 ENDIF
836 IDOC=MINT(83)+MINT(4)
837 DO 450 I=NSD(JT)+1,NSD(JT)+2
838 I1=MINT(83)+MINT(4)+1
839 K(I,3)=I1
840 IF(MSTP(128).GE.1) K(I,3)=ID
841 IF(MSTP(128).LE.1.AND.MINT(4).LT.MSTP(126)) THEN
842 MINT(4)=MINT(4)+1
843 K(I1,1)=21
844 K(I1,2)=K(I,2)
845 K(I1,3)=IREF(IP,JT+3)
846 DO 440 J=1,5
847 P(I1,J)=P(I,J)
848 440 CONTINUE
849 ENDIF
850 450 CONTINUE
851C...Shower - top currently special case.
852 NSHBEF=N
853 IF(MSTP(71).GE.1.AND.(KDCY(JT).LE.2.OR.KFA.EQ.6.OR.KFA.EQ.7.OR.
854 &KFA.EQ.8)) CALL LUSHOW(NSD(JT)+1,NSD(JT)+2,P(ID,5))
855 NSHAFT=N
856
857C...Check if new resonances were produced.
858 KNSDA1=IABS(K(NSD(JT)+1,2))
859 KNSDA2=IABS(K(NSD(JT)+2,2))
860 IF(KNSDA1.EQ.6.OR.KNSDA2.EQ.6.OR.KNSDA1.EQ.7.OR.KNSDA2.EQ.7.OR.
861 &KNSDA1.EQ.8.OR.KNSDA2.EQ.8) THEN
862 NSD1=0
863 NSD2=0
864 DO 460 I=NSD(JT)+1,N
865 IF(K(I,1).LT.10.AND.K(I,2).EQ.K(NSD(JT)+1,2)) NSD1=I
866 IF(K(I,1).LT.10.AND.K(I,2).EQ.K(NSD(JT)+2,2)) NSD2=I
867 460 CONTINUE
868 IF(NSD1.NE.0.AND.NSD2.NE.0) THEN
869 NP=NP+1
870 IREF(NP,1)=NSD1
871 IREF(NP,2)=NSD2
872 IREF(NP,3)=0
873 IREF(NP,4)=IDOC+1
874 IREF(NP,5)=IDOC+2
875 IREF(NP,6)=0
876 IREF(NP,7)=K(IREF(IP,JT),2)
877 IREF(NP,8)=IREF(IP,JT)
878 ENDIF
879 ELSEIF(KDCY(JT).EQ.3.OR.KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8) THEN
880 NP=NP+1
881 IREF(NP,1)=NSD(JT)+1
882 IREF(NP,2)=NSD(JT)+2
883 IF(NSHAFT-NSHBEF.GT.0) THEN
884 IREF(NP,1)=NSHBEF+2
885 IREF(NP,2)=NSHBEF+3
886 ENDIF
887 IREF(NP,3)=0
888 IREF(NP,4)=IDOC+1
889 IREF(NP,5)=IDOC+2
890 IREF(NP,6)=0
891 IREF(NP,7)=K(IREF(IP,JT),2)
892 IREF(NP,8)=IREF(IP,JT)
893 ENDIF
894 470 CONTINUE
895
896C...Fill information for 2 -> 1 -> 2. Loop back if needed.
897 IF(JTMAX.EQ.1.AND.KDCY(1).NE.0) THEN
898 MINT(7)=MINT(83)+6+2*ISET(ISUB)
899 MINT(8)=MINT(83)+7+2*ISET(ISUB)
900 MINT(25)=KFL1(1)
901 MINT(26)=KFL2(1)
902 VINT(23)=CTHE(1)
903 RM3=P(N-1,5)**2/SH
904 RM4=P(N,5)**2/SH
905 BE34=SQRT(MAX(0.,(1.-RM3-RM4)**2-4.*RM3*RM4))
906 VINT(45)=-0.5*SH*(1.-RM3-RM4-BE34*CTHE(1))
907 VINT(46)=-0.5*SH*(1.-RM3-RM4+BE34*CTHE(1))
908 VINT(48)=0.25*SH*BE34**2*MAX(0.,1.-CTHE(1)**2)
909 VINT(47)=SQRT(VINT(48))
910 ENDIF
911 480 IF(IP.LT.NP) GOTO 130
912
913 RETURN
914 END