2 C*********************************************************************
6 C...Allows resonances to decay (including parton showers for hadronic
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)
25 C...The F, Xi and Xj functions of Gunion and Kunszt
26 C...(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))
35 C...Some general constants.
38 IF(MSTP(8).GE.2) XW=1.-(PMAS(24,1)/PMAS(23,1))**2
44 C...Define initial one, two or three objects.
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
64 ELSEIF(ISET(ISUB).EQ.6) THEN
73 C...Check if initial resonance has been moved (in resonance + jet).
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
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
90 C...Loop over decay history.
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
101 C...Start treatment of one or two resonances in parallel.
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))
123 C...Select decay channel.
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
132 IPM=(5-ISIGN(1,K(ID,2)))/2
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.
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
157 C...Read out and classify decay channel chosen.
158 KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
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))
164 IF(KFC2A.GT.100) KFC2A=LUCOMP(KFC2A)
165 IF(KCHG(KFC2A,3).EQ.0) KFL2(JT)=IABS(KFL2(JT))
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)
175 C...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')
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
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
193 C...Fill decay products, prepared for parton showers for quarks.
194 C...Special cases, done by hand, for techni-eta, t, leptoquark and q*.
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
199 CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
201 IF(K(ID,2).LT.0) ISID=5
203 IF(KFC1A.EQ.21.AND.RLU(0).GT.0.5) ISID=9-ISID
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)))
216 K(ID,ISID)=K(ID,ISID)+N
218 ELSEIF(KFA.EQ.39) THEN
221 K(ID,ISID)=K(ID,ISID)+(N-1)
222 K(N-1,ISID)=MSTU(5)*ID
223 ELSEIF(KFL1(JT).NE.21) THEN
226 K(ID,ISID)=K(ID,ISID)+N
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)
236 ELSEIF(KDCY(JT).EQ.1) THEN
237 CALL LU2ENT(-(N+1),KFL1(JT),KFL2(JT),P(ID,5))
239 CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
242 160 IF(KFA.GE.23.AND.KFA.LE.40.AND.KFL1(JT).EQ.0) NINH=NINH+1
245 C...Check for allowed combinations. Skip if no decays.
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
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
255 C...Order incoming partons and outgoing resonances.
256 IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
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)
262 IF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.IREF(IP,7)
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
272 C...Order decay products of resonances.
273 DO 180 JT=IORD,3-IORD,3-2*IORD
274 IF(KDCY(JT).EQ.0) THEN
277 ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
278 ILIN(IMAX+1)=N+2*JT-1
281 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
282 K(N+2*JT,2)=K(NSD(JT)+2,2)
285 ILIN(IMAX+2)=N+2*JT-1
287 K(N+2*JT-1,2)=K(NSD(JT)+1,2)
288 K(N+2*JT,2)=K(NSD(JT)+2,2)
292 C...Find charge, isospin, left- and righthanded couplings.
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)
305 C...Full propagator dependence and flavour correlations for 2 gamma*/Z.
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
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
326 C...Select angular orientation type - Z'/W' only.
329 IF(RLU(0).LT.PARU(130)) MZPWP=1
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
336 ELSEIF(ISUB.EQ.142) THEN
337 IF(RLU(0).LT.PARU(136)) MZPWP=1
339 IAKIR=IABS(K(IREF(2,2),2))
340 IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
345 C...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
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)
354 CTHE(JT)=2.*RLU(0)-1.
355 PHI(JT)=PARU(2)*RLU(0)
359 IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
360 C...Construct massless four-vectors.
369 IF(KDCY(JT).EQ.0) GOTO 280
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)))
379 C...Store incoming and outgoing momenta, with random rotation to
380 C...avoid accidental zeroes in HA expressions.
383 P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
385 P(N+4+I,5)=P(ILIN(I),5)
387 P(N+4+I,J)=P(ILIN(I),J)
390 310 THERR=ACOS(2.*RLU(0)-1.)
392 CALL LUDBRB(N+5,N+4+IMAX,THERR,PHIRR,0D0,0D0,0D0)
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
400 C...Calculate internal products.
401 IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25.OR.ISUB.EQ.141.OR.
403 DO 350 I1=IMIN,IMAX-1
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)
422 DO 390 I1=IMIN,IMAX-1
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)
431 KFAGM=IABS(IREF(IP,7))
432 IF(MSTP(47).LE.0.OR.NINH.NE.0) THEN
433 C...Isotropic decay selected by user.
437 ELSEIF(ITLH.GE.1) THEN
438 C... Isotropic decay t -> b + W etc for 4th generation q and l.
442 ELSEIF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.
443 &IREF(IP,7).EQ.36) THEN
444 C...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
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
454 C...Angular correlation in f -> f' + W -> f' + 2 quarks/leptons.
456 IF(MOD(KFAGM,2).EQ.0) THEN
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
471 ELSEIF(ISUB.EQ.1) THEN
472 C...Angular weight for gamma*/Z0 -> 2 quarks/leptons.
473 EI=KCHG(IABS(MINT(15)),1)/3.
476 EF=KCHG(IABS(KFL1(1)),1)/3.
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
485 ELSEIF(ISUB.EQ.2) THEN
486 C...Angular weight for W+/- -> 2 quarks/leptons.
487 WT=(1.+CTHE(1)*ISIGN(1,MINT(15)*KFL1(1)))**2
490 ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
491 C...Angular weight for f + f~ -> gluon/gamma + (gamma*/Z0) ->
492 C...-> 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)
510 ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
511 C...Angular weight for f + f~' -> gluon/gamma + W+/- ->
512 C...-> 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
516 ELSEIF(ISUB.EQ.22) THEN
517 C...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
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+
540 ELSEIF(ISUB.EQ.23) THEN
541 C...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))
555 ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.171.OR.ISUB.EQ.176) THEN
556 C...Angular weight for f + f~ -> Z0 + H0 -> 2 quarks/leptons + H0
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))
564 ELSEIF(ISUB.EQ.25) THEN
565 C...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)))
581 ELSEIF(ISUB.EQ.26.OR.ISUB.EQ.172.OR.ISUB.EQ.177) THEN
582 C...Angular weight for f + f~' -> W+/- + H0 -> 2 quarks/leptons + H0
585 WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
587 ELSEIF(ISUB.EQ.30.OR.ISUB.EQ.35) THEN
588 C...Angular weight for f + g/gamma -> f + (gamma*/Z0)
589 C...-> 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)
609 ELSEIF(ISUB.EQ.31) THEN
610 C...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
615 ELSEIF(ISUB.EQ.71.OR.ISUB.EQ.72.OR.ISUB.EQ.73.OR.ISUB.EQ.76.OR.
617 C...Angular weight for V_L1 + V_L2 -> V_L3 + V_L4 (V = Z/W).
618 WT=16.*PKK(3,5)*PKK(4,6)
621 ELSEIF(ISUB.EQ.110) THEN
622 C...Angular weight for f + f~ -> gamma + H0 -> gamma + X is isotropic.
626 ELSEIF(ISUB.EQ.141) THEN
627 IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
628 C...Angular weight for f + f~ -> gamma*/Z0/Z'0 -> 2 quarks/leptons.
629 C...Couplings of incoming flavour.
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)
640 C...Couplings of final flavour.
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)
651 C...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
661 ELSEIF(IP.EQ.1.AND.IABS(KFL1(1)).EQ.24) THEN
662 C...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
672 C...Angular weight for f + f~ -> Z' -> H0 + A0, H'0 + A0, H+ + H-.
675 ELSEIF(IP.EQ.1.AND.KDCY(1).EQ.3) THEN
676 C...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
683 C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
684 C...(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
695 C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
696 C...(W:s approximately longitudinal, like if intermediate H).
697 WT=16.*PKK(3,5)*PKK(4,6)
700 C...Angular weight for f + f~ -> Z' -> H+ + H-, Z0 + H0, H0 + A0,
701 C...H'0 + A0 -> 4 quarks/leptons.
706 ELSEIF(ISUB.EQ.142) THEN
707 IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
708 C...Angular weight for f + f~' -> W'+/- -> 2 quarks/leptons.
711 IF(KFAI.GT.10) KFAIC=2
716 IF(KFAF.GT.10) KFAFC=2
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
722 ELSEIF(IP.EQ.1.AND.IABS(KFL2(1)).EQ.23) THEN
723 C...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
732 C...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
739 C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
740 C...(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
751 C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
752 C...(W/Z approximately longitudinal, like if intermediate H).
753 WT=16.*PKK(3,5)*PKK(4,6)
756 C...Angular weight for f + f~ -> W' -> W + H0 -> whatever.
761 ELSEIF(ISUB.EQ.145.OR.ISUB.EQ.162.OR.ISUB.EQ.163.OR.ISUB.EQ.164)
763 C...Isotropic decay of leptoquarks (assumed spin 0).
767 ELSEIF(ISUB.EQ.147.OR.ISUB.EQ.148) THEN
768 C...Decays of (spin 1/2) q* -> q + (g,gamma) or (Z0,W+-).
770 IF(MINT(16).EQ.21) SIDE=-1.
771 IF(IP.EQ.1.AND.(KFL1(1).EQ.21.OR.KFL1(1).EQ.22)) 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)
779 C...W/Z decay assumed isotropic, since not known.
784 ELSEIF(ISUB.EQ.149) THEN
785 C...Isotropic decay of techni-eta.
789 C...Obtain correct angular distribution by rejection techniques.
794 IF(WT.LT.RLU(0)*WTMAX) GOTO 240
796 C...Construct massive four-vectors using angles chosen. Mark decayed
797 C...resonances, add documentation lines. Shower evolution.
798 400 DO 470 JT=1,JTMAX
799 IF(KDCY(JT).EQ.0) GOTO 470
801 IF(ISET(ISUB).NE.6.OR.JT.NE.3) THEN
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))
809 C...Z + q + q~ : angles already known, in rest frame of system.
811 DBEZQQ(J)=(P(ID,J)+P(ID+1,J)+P(ID+2,J))/(P(ID,4)+P(ID+1,4)+
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)/
824 CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,THEZQQ,PHIZQQ,DBEZQQ(1),
825 & DBEZQQ(2),DBEZQQ(3))
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
831 C...Do not kill colour flow through techni-eta, t, leptoquark or q*!
836 IDOC=MINT(83)+MINT(4)
837 DO 450 I=NSD(JT)+1,NSD(JT)+2
838 I1=MINT(83)+MINT(4)+1
840 IF(MSTP(128).GE.1) K(I,3)=ID
841 IF(MSTP(128).LE.1.AND.MINT(4).LT.MSTP(126)) THEN
845 K(I1,3)=IREF(IP,JT+3)
851 C...Shower - top currently special case.
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))
857 C...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
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
868 IF(NSD1.NE.0.AND.NSD2.NE.0) THEN
876 IREF(NP,7)=K(IREF(IP,JT),2)
877 IREF(NP,8)=IREF(IP,JT)
879 ELSEIF(KDCY(JT).EQ.3.OR.KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8) THEN
883 IF(NSHAFT-NSHBEF.GT.0) THEN
891 IREF(NP,7)=K(IREF(IP,JT),2)
892 IREF(NP,8)=IREF(IP,JT)
896 C...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)
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))
911 480 IF(IP.LT.NP) GOTO 130