+++ /dev/null
-
-C*********************************************************************
-
- SUBROUTINE PYRESD
-
-C...Allows resonances to decay (including parton showers for hadronic
-C...channels).
- IMPLICIT DOUBLE PRECISION(D)
- COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
- COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
- COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
- COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
- COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
- COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
- COMMON/PYINT1/MINT(400),VINT(400)
- COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
- COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
- SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
- SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT4/
- DIMENSION IREF(20,8),KDCY(3),KFL1(3),KFL2(3),KEQL(3),NSD(3),
- &ILIN(6),HGZ(3,3),COUP(6,4),CORL(2,2,2),PK(6,4),PKK(6,6),
- &CTHE(3),PHI(3),WDTP(0:40),WDTE(0:40,0:5),DBEZQQ(3),DPMO(5)
- COMPLEX FGK,HA(6,6),HC(6,6)
-
-C...The F, Xi and Xj functions of Gunion and Kunszt
-C...(Phys. Rev. D33, 665, plus errata from the authors).
- FGK(I1,I2,I3,I4,I5,I6)=4.*HA(I1,I3)*HC(I2,I6)*(HA(I1,I5)*
- &HC(I1,I4)+HA(I3,I5)*HC(I3,I4))
- DIGK(DT,DU)=-4.*D34*D56+DT*(3.*DT+4.*DU)+DT**2*(DT*DU/(D34*D56)-
- &2.*(1./D34+1./D56)*(DT+DU)+2.*(D34/D56+D56/D34))
- DJGK(DT,DU)=8.*(D34+D56)**2-8.*(D34+D56)*(DT+DU)-6.*DT*DU-
- &2.*DT*DU*(DT*DU/(D34*D56)-2.*(1./D34+1./D56)*(DT+DU)+
- &2.*(D34/D56+D56/D34))
-
-C...Some general constants.
- XW=PARU(102)
- XWV=XW
- IF(MSTP(8).GE.2) XW=1.-(PMAS(24,1)/PMAS(23,1))**2
- XW1=1.-XW
- SQMZ=PMAS(23,1)**2
- SQMW=PMAS(24,1)**2
- SH=VINT(44)
-
-C...Define initial one, two or three objects.
- ISUB=MINT(1)
- DO 100 JT=1,8
- IREF(1,JT)=0
- 100 CONTINUE
- IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
- IREF(1,1)=MINT(84)+2+ISET(ISUB)
- IREF(1,4)=MINT(83)+6+ISET(ISUB)
- ELSEIF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
- IREF(1,1)=MINT(84)+1+ISET(ISUB)
- IREF(1,2)=MINT(84)+2+ISET(ISUB)
- IREF(1,4)=MINT(83)+5+ISET(ISUB)
- IREF(1,5)=MINT(83)+6+ISET(ISUB)
- ELSEIF(ISET(ISUB).EQ.5) THEN
- IREF(1,1)=MINT(84)+3
- IREF(1,2)=MINT(84)+4
- IREF(1,3)=MINT(84)+5
- IREF(1,4)=MINT(83)+7
- IREF(1,5)=MINT(83)+8
- IREF(1,6)=MINT(83)+9
- ELSEIF(ISET(ISUB).EQ.6) THEN
- IREF(1,1)=MINT(84)+4
- IREF(1,2)=MINT(84)+5
- IREF(1,3)=MINT(84)+3
- IREF(1,4)=MINT(83)+8
- IREF(1,5)=MINT(83)+9
- IREF(1,6)=MINT(83)+7
- ENDIF
-
-C...Check if initial resonance has been moved (in resonance + jet).
- DO 120 JT=1,3
- IF(IREF(1,JT).GT.0) THEN
- IF(K(IREF(1,JT),1).GT.10) THEN
- KFA=IABS(K(IREF(1,JT),2))
- IF(KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8.OR.KFA.EQ.39) THEN
- DO 110 I=IREF(1,JT)+1,N
- IF(K(I,1).LE.10.AND.K(I,2).EQ.K(IREF(1,JT),2)) IREF(1,JT)=I
- 110 CONTINUE
- ELSE
- KDA=MOD(K(IREF(1,JT),4),MSTU(4))
- IF(KFA.GE.23.AND.KFA.LE.40.AND.KDA.GT.1) IREF(1,JT)=KDA
- ENDIF
- ENDIF
- ENDIF
- 120 CONTINUE
-
-C...Loop over decay history.
- NP=1
- IP=0
- 130 IP=IP+1
- NINH=0
- JTMAX=2
- IF(IP.EQ.1.AND.IREF(1,2).EQ.0) JTMAX=1
- IF(IP.EQ.1.AND.IREF(1,3).NE.0) JTMAX=3
- ITLH=0
- NSAV=N
-
-C...Start treatment of one or two resonances in parallel.
- 140 N=NSAV
- DO 170 JT=1,JTMAX
- ID=IREF(IP,JT)
- KDCY(JT)=0
- KFL1(JT)=0
- KFL2(JT)=0
- KEQL(JT)=0
- NSD(JT)=ID
- IF(ID.EQ.0) GOTO 160
- KFA=IABS(K(ID,2))
- IF((KFA.LT.23.OR.KFA.GT.40).AND.KFA.NE.6.AND.KFA.NE.7.AND.
- &KFA.NE.8.AND.KFA.NE.17.AND.KFA.NE.18) GOTO 160
- IF(MSTP(48).LE.0.AND.KFA.EQ.6) GOTO 160
- IF(MSTP(6).NE.1.AND.MSTP(49).LE.0.AND.(KFA.EQ.7.OR.KFA.EQ.8.OR.
- &KFA.EQ.17.OR.KFA.EQ.18)) GOTO 160
- IF(K(ID,1).GT.10.OR.MDCY(KFA,1).EQ.0) GOTO 160
- IF(KFA.EQ.6.OR.(MSTP(6).NE.1.AND.(KFA.EQ.7.OR.KFA.EQ.8.OR.
- &KFA.EQ.17.OR.KFA.EQ.18))) ITLH=ITLH+1
- K(ID,4)=MSTU(5)*(K(ID,4)/MSTU(5))
- K(ID,5)=MSTU(5)*(K(ID,5)/MSTU(5))
-
-C...Select decay channel.
- KFB=0
- IF(ISET(ISUB).NE.6.OR.JT.NE.3) THEN
- IF(ISUB.EQ.1.OR.ISUB.EQ.15.OR.ISUB.EQ.19.OR.ISUB.EQ.22.OR.
- & ISUB.EQ.30.OR.ISUB.EQ.35.OR.ISUB.EQ.141) MINT(61)=1
- CALL PYWIDT(KFA,P(ID,5)**2,WDTP,WDTE)
- IF(KCHG(KFA,3).EQ.0) THEN
- IPM=2
- ELSE
- IPM=(5-ISIGN(1,K(ID,2)))/2
- ENDIF
- IF(JTMAX.GE.2.AND.JT.LE.2) KFB=IABS(K(IREF(IP,3-JT),2))
- WDTE0S=WDTE(0,1)+WDTE(0,IPM)+WDTE(0,4)
- IF(KFB.EQ.KFA) WDTE0S=WDTE0S+WDTE(0,5)
- IF(WDTE0S.LE.0.) THEN
- IF(KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8.OR.KFA.EQ.17.OR.
- & KFA.EQ.18) THEN
- MINT(51)=1
- RETURN
- ELSE
- GOTO 160
- ENDIF
- ENDIF
- RKFL=WDTE0S*RLU(0)
- IDL=0
- 150 IDL=IDL+1
- IDC=IDL+MDCY(KFA,2)-1
- RKFL=RKFL-(WDTE(IDL,1)+WDTE(IDL,IPM)+WDTE(IDL,4))
- IF(KFB.EQ.KFA) RKFL=RKFL-WDTE(IDL,5)
- IF(IDL.LT.MDCY(KFA,3).AND.RKFL.GT.0.) GOTO 150
- ELSE
- IDC=MINT(35)
- ENDIF
-
-C...Read out and classify decay channel chosen.
- KFL1(JT)=KFDP(IDC,1)*ISIGN(1,K(ID,2))
- KFC1A=IABS(KFL1(JT))
- IF(KFC1A.GT.100) KFC1A=LUCOMP(KFC1A)
- IF(KCHG(KFC1A,3).EQ.0) KFL1(JT)=IABS(KFL1(JT))
- KFL2(JT)=KFDP(IDC,2)*ISIGN(1,K(ID,2))
- KFC2A=IABS(KFL2(JT))
- IF(KFC2A.GT.100) KFC2A=LUCOMP(KFC2A)
- IF(KCHG(KFC2A,3).EQ.0) KFL2(JT)=IABS(KFL2(JT))
- KDCY(JT)=2
- IF(IABS(KFL1(JT)).LE.10.OR.KFL1(JT).EQ.21) KDCY(JT)=1
- IF(IABS(KFL1(JT)).GE.23.AND.IABS(KFL1(JT)).LE.40) KDCY(JT)=3
- IF(KFB.EQ.KFA) KEQL(JT)=MDME(IDC,1)
- NSD(JT)=N
- HGZ(JT,1)=VINT(111)
- HGZ(JT,2)=VINT(112)
- HGZ(JT,3)=VINT(114)
-
-C...Select masses and check that mass sum not too large.
- IF(MSTP(42).LE.0.OR.(PMAS(KFC1A,2).LT.PARP(41).AND.
- &PMAS(KFC2A,2).LT.PARP(41))) THEN
- P(N+1,5)=PMAS(KFC1A,1)
- P(N+2,5)=PMAS(KFC2A,1)
- IF(P(N+1,5)+P(N+2,5)+PARJ(64).GT.P(ID,5)) THEN
- CALL LUERRM(13,'(PYRESD:) daughter masses too large')
- MINT(51)=1
- RETURN
- ENDIF
- ELSEIF(IP.EQ.1) THEN
- CALL PYOFSH(2,KFA,KFL1(JT),KFL2(JT),P(ID,5),P(N+1,5),P(N+2,5))
- IF(MINT(51).EQ.1) RETURN
- ELSE
- CALL PYOFSH(7,KFA,KFL1(JT),KFL2(JT),P(ID,5),P(N+1,5),P(N+2,5))
- IF(MINT(51).EQ.1) RETURN
- ENDIF
-
-C...Fill decay products, prepared for parton showers for quarks.
-C...Special cases, done by hand, for techni-eta, t, leptoquark and q*.
- MSTU(10)=1
- IF(KFA.EQ.38.OR.KFA.EQ.39.OR.((MSTP(6).EQ.1.OR.MSTP(49).GE.1).AND.
- &(KFA.EQ.7.OR.KFA.EQ.8)).OR.KFA.EQ.6) THEN
- MSTU(19)=1
- CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
- ISID=4
- IF(K(ID,2).LT.0) ISID=5
- IF(KFA.EQ.38) THEN
- IF(KFC1A.EQ.21.AND.RLU(0).GT.0.5) ISID=9-ISID
- K(N-1,1)=3
- K(N,1)=3
- K(ID,ISID)=K(ID,ISID)+(N-1)
- K(ID,9-ISID)=K(ID,9-ISID)+N
- K(N-1,ISID)=MSTU(5)*ID
- K(N-1,9-ISID)=MSTU(5)*N
- K(N,ISID)=MSTU(5)*(N-1)
- K(N,9-ISID)=MSTU(5)*ID
- ELSEIF(KFA.EQ.6.OR.(MSTP(6).NE.1.AND.(KFA.EQ.7.OR.KFA.EQ.8)))
- & THEN
- K(N-1,1)=1
- K(N,1)=3
- K(ID,ISID)=K(ID,ISID)+N
- K(N,ISID)=MSTU(5)*ID
- ELSEIF(KFA.EQ.39) THEN
- K(N-1,1)=3
- K(N,1)=1
- K(ID,ISID)=K(ID,ISID)+(N-1)
- K(N-1,ISID)=MSTU(5)*ID
- ELSEIF(KFL1(JT).NE.21) THEN
- K(N-1,1)=1
- K(N,1)=3
- K(ID,ISID)=K(ID,ISID)+N
- K(N,ISID)=MSTU(5)*ID
- ELSE
- K(N-1,1)=3
- K(N,1)=3
- K(ID,ISID)=K(ID,ISID)+(N-1)
- K(N-1,ISID)=MSTU(5)*ID
- K(N-1,9-ISID)=MSTU(5)*N
- K(N,ISID)=MSTU(5)*(N-1)
- ENDIF
- ELSEIF(KDCY(JT).EQ.1) THEN
- CALL LU2ENT(-(N+1),KFL1(JT),KFL2(JT),P(ID,5))
- ELSE
- CALL LU2ENT(N+1,KFL1(JT),KFL2(JT),P(ID,5))
- ENDIF
- MSTU(10)=2
- 160 IF(KFA.GE.23.AND.KFA.LE.40.AND.KFL1(JT).EQ.0) NINH=NINH+1
- 170 CONTINUE
-
-C...Check for allowed combinations. Skip if no decays.
- IF(JTMAX.GE.2) THEN
- IF(KEQL(1).EQ.4.AND.KEQL(2).EQ.4) GOTO 140
- IF(KEQL(1).EQ.5.AND.KEQL(2).EQ.5) GOTO 140
- ENDIF
- IF(JTMAX.EQ.1.AND.KDCY(1).EQ.0) GOTO 480
- IF(JTMAX.EQ.2.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0) GOTO 480
- IF(JTMAX.EQ.3.AND.KDCY(1).EQ.0.AND.KDCY(2).EQ.0.AND.
- &KDCY(3).EQ.0) GOTO 480
-
-C...Order incoming partons and outgoing resonances.
- IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
- ILIN(1)=MINT(84)+1
- IF(K(MINT(84)+1,2).GT.0) ILIN(1)=MINT(84)+2
- IF(K(ILIN(1),2).EQ.21) ILIN(1)=2*MINT(84)+3-ILIN(1)
- ILIN(2)=2*MINT(84)+3-ILIN(1)
- IMIN=1
- IF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.IREF(IP,7)
- & .EQ.36) IMIN=3
- IMAX=2
- IORD=1
- IF(K(IREF(IP,1),2).EQ.23) IORD=2
- IF(K(IREF(IP,1),2).EQ.24.AND.K(IREF(IP,2),2).EQ.-24) IORD=2
- IAKIPD=IABS(K(IREF(IP,IORD),2))
- IF(IAKIPD.EQ.25.OR.IAKIPD.EQ.35.OR.IAKIPD.EQ.36) IORD=3-IORD
- IF(KDCY(IORD).EQ.0) IORD=3-IORD
-
-C...Order decay products of resonances.
- DO 180 JT=IORD,3-IORD,3-2*IORD
- IF(KDCY(JT).EQ.0) THEN
- ILIN(IMAX+1)=NSD(JT)
- IMAX=IMAX+1
- ELSEIF(K(NSD(JT)+1,2).GT.0) THEN
- ILIN(IMAX+1)=N+2*JT-1
- ILIN(IMAX+2)=N+2*JT
- IMAX=IMAX+2
- K(N+2*JT-1,2)=K(NSD(JT)+1,2)
- K(N+2*JT,2)=K(NSD(JT)+2,2)
- ELSE
- ILIN(IMAX+1)=N+2*JT
- ILIN(IMAX+2)=N+2*JT-1
- IMAX=IMAX+2
- K(N+2*JT-1,2)=K(NSD(JT)+1,2)
- K(N+2*JT,2)=K(NSD(JT)+2,2)
- ENDIF
- 180 CONTINUE
-
-C...Find charge, isospin, left- and righthanded couplings.
- DO 200 I=IMIN,IMAX
- DO 190 J=1,4
- COUP(I,J)=0.
- 190 CONTINUE
- KFA=IABS(K(ILIN(I),2))
- IF(KFA.EQ.0.OR.KFA.GT.20) GOTO 200
- COUP(I,1)=KCHG(KFA,1)/3.
- COUP(I,2)=(-1)**MOD(KFA,2)
- COUP(I,4)=-2.*COUP(I,1)*XWV
- COUP(I,3)=COUP(I,2)+COUP(I,4)
- 200 CONTINUE
-
-C...Full propagator dependence and flavour correlations for 2 gamma*/Z.
- IF(ISUB.EQ.22) THEN
- DO 230 I=3,5,2
- I1=IORD
- IF(I.EQ.5) I1=3-IORD
- DO 220 J1=1,2
- DO 210 J2=1,2
- CORL(I/2,J1,J2)=COUP(1,1)**2*HGZ(I1,1)*COUP(I,1)**2/16.+
- & COUP(1,1)*COUP(1,J1+2)*HGZ(I1,2)*COUP(I,1)*COUP(I,J2+2)/4.+
- & COUP(1,J1+2)**2*HGZ(I1,3)*COUP(I,J2+2)**2
- 210 CONTINUE
- 220 CONTINUE
- 230 CONTINUE
- COWT12=(CORL(1,1,1)+CORL(1,1,2))*(CORL(2,1,1)+CORL(2,1,2))+
- & (CORL(1,2,1)+CORL(1,2,2))*(CORL(2,2,1)+CORL(2,2,2))
- COMX12=(CORL(1,1,1)+CORL(1,1,2)+CORL(1,2,1)+CORL(1,2,2))*
- & (CORL(2,1,1)+CORL(2,1,2)+CORL(2,2,1)+CORL(2,2,2))
- IF(COWT12.LT.RLU(0)*COMX12) GOTO 140
- ENDIF
- ENDIF
-
-C...Select angular orientation type - Z'/W' only.
- MZPWP=0
- IF(ISUB.EQ.141) THEN
- IF(RLU(0).LT.PARU(130)) MZPWP=1
- IF(IP.EQ.2) THEN
- IF(IABS(K(IREF(2,1),2)).EQ.37) MZPWP=2
- IAKIR=IABS(K(IREF(2,2),2))
- IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
- ENDIF
- IF(IP.GE.3) MZPWP=2
- ELSEIF(ISUB.EQ.142) THEN
- IF(RLU(0).LT.PARU(136)) MZPWP=1
- IF(IP.EQ.2) THEN
- IAKIR=IABS(K(IREF(2,2),2))
- IF(IAKIR.EQ.25.OR.IAKIR.EQ.35.OR.IAKIR.EQ.36) MZPWP=2
- ENDIF
- IF(IP.GE.3) MZPWP=2
- ENDIF
-
-C...Select random angles (begin of weighting procedure).
- 240 DO 250 JT=1,JTMAX
- IF(KDCY(JT).EQ.0) GOTO 250
- IF(ISET(ISUB).EQ.6.AND.JT.EQ.3) GOTO 250
- IF(JTMAX.EQ.1) THEN
- CTHE(JT)=VINT(13)+(VINT(33)-VINT(13)+VINT(34)-VINT(14))*RLU(0)
- IF(CTHE(JT).GT.VINT(33)) CTHE(JT)=CTHE(JT)+VINT(14)-VINT(33)
- PHI(JT)=VINT(24)
- ELSE
- CTHE(JT)=2.*RLU(0)-1.
- PHI(JT)=PARU(2)*RLU(0)
- ENDIF
- 250 CONTINUE
-
- IF(JTMAX.EQ.2.AND.MSTP(47).GE.1.AND.NINH.EQ.0) THEN
-C...Construct massless four-vectors.
- DO 270 I=N+1,N+4
- K(I,1)=1
- DO 260 J=1,5
- P(I,J)=0.
- V(I,J)=0.
- 260 CONTINUE
- 270 CONTINUE
- DO 280 JT=1,JTMAX
- IF(KDCY(JT).EQ.0) GOTO 280
- ID=IREF(IP,JT)
- P(N+2*JT-1,3)=0.5*P(ID,5)
- P(N+2*JT-1,4)=0.5*P(ID,5)
- P(N+2*JT,3)=-0.5*P(ID,5)
- P(N+2*JT,4)=0.5*P(ID,5)
- CALL LUDBRB(N+2*JT-1,N+2*JT,ACOS(CTHE(JT)),PHI(JT),DBLE(P(ID,1)/
- & P(ID,4)),DBLE(P(ID,2)/P(ID,4)),DBLE(P(ID,3)/P(ID,4)))
- 280 CONTINUE
-
-C...Store incoming and outgoing momenta, with random rotation to
-C...avoid accidental zeroes in HA expressions.
- DO 300 I=1,IMAX
- K(N+4+I,1)=1
- P(N+4+I,4)=SQRT(P(ILIN(I),1)**2+P(ILIN(I),2)**2+P(ILIN(I),3)**2+
- & P(ILIN(I),5)**2)
- P(N+4+I,5)=P(ILIN(I),5)
- DO 290 J=1,3
- P(N+4+I,J)=P(ILIN(I),J)
- 290 CONTINUE
- 300 CONTINUE
- 310 THERR=ACOS(2.*RLU(0)-1.)
- PHIRR=PARU(2)*RLU(0)
- CALL LUDBRB(N+5,N+4+IMAX,THERR,PHIRR,0D0,0D0,0D0)
- DO 330 I=1,IMAX
- IF(P(N+4+I,1)**2+P(N+4+I,2)**2.LT.1E-4*P(N+4+I,4)**2) GOTO 310
- DO 320 J=1,4
- PK(I,J)=P(N+4+I,J)
- 320 CONTINUE
- 330 CONTINUE
-
-C...Calculate internal products.
- IF(ISUB.EQ.22.OR.ISUB.EQ.23.OR.ISUB.EQ.25.OR.ISUB.EQ.141.OR.
- & ISUB.EQ.142) THEN
- DO 350 I1=IMIN,IMAX-1
- DO 340 I2=I1+1,IMAX
- HA(I1,I2)=SQRT((PK(I1,4)-PK(I1,3))*(PK(I2,4)+PK(I2,3))/
- & (1E-20+PK(I1,1)**2+PK(I1,2)**2))*CMPLX(PK(I1,1),PK(I1,2))-
- & SQRT((PK(I1,4)+PK(I1,3))*(PK(I2,4)-PK(I2,3))/
- & (1E-20+PK(I2,1)**2+PK(I2,2)**2))*CMPLX(PK(I2,1),PK(I2,2))
- HC(I1,I2)=CONJG(HA(I1,I2))
- IF(I1.LE.2) HA(I1,I2)=CMPLX(0.,1.)*HA(I1,I2)
- IF(I1.LE.2) HC(I1,I2)=CMPLX(0.,1.)*HC(I1,I2)
- HA(I2,I1)=-HA(I1,I2)
- HC(I2,I1)=-HC(I1,I2)
- 340 CONTINUE
- 350 CONTINUE
- ENDIF
- DO 370 I=1,2
- DO 360 J=1,4
- PK(I,J)=-PK(I,J)
- 360 CONTINUE
- 370 CONTINUE
- DO 390 I1=IMIN,IMAX-1
- DO 380 I2=I1+1,IMAX
- PKK(I1,I2)=2.*(PK(I1,4)*PK(I2,4)-PK(I1,1)*PK(I2,1)-
- & PK(I1,2)*PK(I2,2)-PK(I1,3)*PK(I2,3))
- PKK(I2,I1)=PKK(I1,I2)
- 380 CONTINUE
- 390 CONTINUE
- ENDIF
-
- KFAGM=IABS(IREF(IP,7))
- IF(MSTP(47).LE.0.OR.NINH.NE.0) THEN
-C...Isotropic decay selected by user.
- WT=1.
- WTMAX=1.
-
- ELSEIF(ITLH.GE.1) THEN
-C... Isotropic decay t -> b + W etc for 4th generation q and l.
- WT=1.
- WTMAX=1.
-
- ELSEIF(IREF(IP,7).EQ.25.OR.IREF(IP,7).EQ.35.OR.
- &IREF(IP,7).EQ.36) THEN
-C...Angular weight for H0 -> Z0 + Z0 or W+ + W- -> 4 quarks/leptons.
- WT=16.*PKK(3,5)*PKK(4,6)
- IF(IP.EQ.1) WTMAX=SH**2
- IF(IP.GE.2) WTMAX=P(IREF(IP,8),5)**4
- KFA=IABS(K(IREF(IP,1),2))
- IF(KFA.NE.23.AND.KFA.NE.24) WT=WTMAX
-
- ELSEIF((KFAGM.EQ.6.OR.(MSTP(6).NE.1.AND.(KFAGM.EQ.7.OR.
- &KFAGM.EQ.8.OR.KFAGM.EQ.17.OR.KFAGM.EQ.18))).AND.
- &IABS(K(IREF(IP,1),2)).EQ.24) THEN
-C...Angular correlation in f -> f' + W -> f' + 2 quarks/leptons.
- I1=IREF(IP,8)
- IF(MOD(KFAGM,2).EQ.0) THEN
- I2=N+1
- I3=N+2
- ELSE
- I2=N+2
- I3=N+1
- ENDIF
- I4=IREF(IP,2)
- WT=(P(I1,4)*P(I2,4)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
- & P(I1,3)*P(I2,3))*(P(I3,4)*P(I4,4)-P(I3,1)*P(I4,1)-
- & P(I3,2)*P(I4,2)-P(I3,3)*P(I4,3))
- WTMAX=(P(I1,5)**4-P(IREF(IP,1),5)**4)/8.
- IF(KFAGM.EQ.6.AND.MSTP(48).LE.1) WT=WTMAX
- IF(KFAGM.NE.6.AND.MSTP(49).LE.1) WT=WTMAX
-
- ELSEIF(ISUB.EQ.1) THEN
-C...Angular weight for gamma*/Z0 -> 2 quarks/leptons.
- EI=KCHG(IABS(MINT(15)),1)/3.
- AI=SIGN(1.,EI+0.1)
- VI=AI-4.*EI*XWV
- EF=KCHG(IABS(KFL1(1)),1)/3.
- AF=SIGN(1.,EF+0.1)
- VF=AF-4.*EF*XWV
- ASYM=2.*(EI*AI*VINT(112)*EF*AF+4.*VI*AI*VINT(114)*VF*AF)/
- & (EI**2*VINT(111)*EF**2+EI*VI*VINT(112)*EF*VF+
- & (VI**2+AI**2)*VINT(114)*(VF**2+AF**2))
- WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
- WTMAX=2.+ABS(ASYM)
-
- ELSEIF(ISUB.EQ.2) THEN
-C...Angular weight for W+/- -> 2 quarks/leptons.
- WT=(1.+CTHE(1)*ISIGN(1,MINT(15)*KFL1(1)))**2
- WTMAX=4.
-
- ELSEIF(ISUB.EQ.15.OR.ISUB.EQ.19) THEN
-C...Angular weight for f + f~ -> gluon/gamma + (gamma*/Z0) ->
-C...-> gluon/gamma + 2 quarks/leptons.
- CLILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
- & COUP(1,3)**2*HGZ(2,3)*COUP(3,3)**2
- CLIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
- & COUP(1,3)**2*HGZ(2,3)*COUP(3,4)**2
- CRILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
- & COUP(1,4)**2*HGZ(2,3)*COUP(3,3)**2
- CRIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
- & COUP(1,4)**2*HGZ(2,3)*COUP(3,4)**2
- WT=(CLILF+CRIRF)*(PKK(1,3)**2+PKK(2,4)**2)+
- & (CLIRF+CRILF)*(PKK(1,4)**2+PKK(2,3)**2)
- WTMAX=(CLILF+CLIRF+CRILF+CRIRF)*
- & ((PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2)
-
- ELSEIF(ISUB.EQ.16.OR.ISUB.EQ.20) THEN
-C...Angular weight for f + f~' -> gluon/gamma + W+/- ->
-C...-> gluon/gamma + 2 quarks/leptons.
- WT=PKK(1,3)**2+PKK(2,4)**2
- WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(2,3)+PKK(2,4))**2
-
- ELSEIF(ISUB.EQ.22) THEN
-C...Angular weight for f + f~ -> Z0 + Z0 -> 4 quarks/leptons.
- S34=P(IREF(IP,IORD),5)**2
- S56=P(IREF(IP,3-IORD),5)**2
- TI=PKK(1,3)+PKK(1,4)+S34
- UI=PKK(1,5)+PKK(1,6)+S56
- FGK135=ABS(FGK(1,2,3,4,5,6)/TI+FGK(1,2,5,6,3,4)/UI)**2
- FGK145=ABS(FGK(1,2,4,3,5,6)/TI+FGK(1,2,5,6,4,3)/UI)**2
- FGK136=ABS(FGK(1,2,3,4,6,5)/TI+FGK(1,2,6,5,3,4)/UI)**2
- FGK146=ABS(FGK(1,2,4,3,6,5)/TI+FGK(1,2,6,5,4,3)/UI)**2
- FGK253=ABS(FGK(2,1,5,6,3,4)/TI+FGK(2,1,3,4,5,6)/UI)**2
- FGK263=ABS(FGK(2,1,6,5,3,4)/TI+FGK(2,1,3,4,6,5)/UI)**2
- FGK254=ABS(FGK(2,1,5,6,4,3)/TI+FGK(2,1,4,3,5,6)/UI)**2
- FGK264=ABS(FGK(2,1,6,5,4,3)/TI+FGK(2,1,4,3,6,5)/UI)**2
- WT=
- & CORL(1,1,1)*CORL(2,1,1)*FGK135+CORL(1,1,2)*CORL(2,1,1)*FGK145+
- & CORL(1,1,1)*CORL(2,1,2)*FGK136+CORL(1,1,2)*CORL(2,1,2)*FGK146+
- & CORL(1,2,1)*CORL(2,2,1)*FGK253+CORL(1,2,2)*CORL(2,2,1)*FGK263+
- & CORL(1,2,1)*CORL(2,2,2)*FGK254+CORL(1,2,2)*CORL(2,2,2)*FGK264
- WTMAX=16.*((CORL(1,1,1)+CORL(1,1,2))*(CORL(2,1,1)+CORL(2,1,2))+
- & (CORL(1,2,1)+CORL(1,2,2))*(CORL(2,2,1)+CORL(2,2,2)))*S34*S56*
- & ((TI**2+UI**2+2.*SH*(S34+S56))/(TI*UI)-S34*S56*(1./TI**2+
- & 1./UI**2))
-
- ELSEIF(ISUB.EQ.23) THEN
-C...Angular weight for f + f~' -> Z0 + W+/- -> 4 quarks/leptons.
- D34=P(IREF(IP,IORD),5)**2
- D56=P(IREF(IP,3-IORD),5)**2
- DT=PKK(1,3)+PKK(1,4)+D34
- DU=PKK(1,5)+PKK(1,6)+D56
- FACBW=1./((SH-SQMW)**2+SQMW*PMAS(24,2)**2)
- CAWZ=COUP(2,3)/SNGL(DT)-2.*XW1*COUP(1,2)*(SH-SQMW)*FACBW
- CBWZ=COUP(1,3)/SNGL(DU)+2.*XW1*COUP(1,2)*(SH-SQMW)*FACBW
- FGK135=ABS(CAWZ*FGK(1,2,3,4,5,6)+CBWZ*FGK(1,2,5,6,3,4))
- FGK136=ABS(CAWZ*FGK(1,2,3,4,6,5)+CBWZ*FGK(1,2,6,5,3,4))
- WT=(COUP(5,3)*FGK135)**2+(COUP(5,4)*FGK136)**2
- WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*(CAWZ**2*
- & DIGK(DT,DU)+CBWZ**2*DIGK(DU,DT)+CAWZ*CBWZ*DJGK(DT,DU))
-
- ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.171.OR.ISUB.EQ.176) THEN
-C...Angular weight for f + f~ -> Z0 + H0 -> 2 quarks/leptons + H0
-C...(or H'0, or A0).
- WT=((COUP(1,3)*COUP(3,3))**2+(COUP(1,4)*COUP(3,4))**2)*
- & PKK(1,3)*PKK(2,4)+((COUP(1,3)*COUP(3,4))**2+(COUP(1,4)*
- & COUP(3,3))**2)*PKK(1,4)*PKK(2,3)
- WTMAX=(COUP(1,3)**2+COUP(1,4)**2)*(COUP(3,3)**2+COUP(3,4)**2)*
- & (PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
-
- ELSEIF(ISUB.EQ.25) THEN
-C...Angular weight for f + f~ -> W+ + W- -> 4 quarks/leptons.
- D34=P(IREF(IP,IORD),5)**2
- D56=P(IREF(IP,3-IORD),5)**2
- DT=PKK(1,3)+PKK(1,4)+D34
- DU=PKK(1,5)+PKK(1,6)+D56
- FACBW=1./((SH-SQMZ)**2+SQMZ*PMAS(23,2)**2)
- CDWW=(COUP(1,3)*SQMZ*(SH-SQMZ)*FACBW+COUP(1,2))/SH
- CAWW=CDWW+0.5*(COUP(1,2)+1.)/SNGL(DT)
- CBWW=CDWW+0.5*(COUP(1,2)-1.)/SNGL(DU)
- CCWW=COUP(1,4)*SQMZ*(SH-SQMZ)*FACBW/SH
- FGK135=ABS(CAWW*FGK(1,2,3,4,5,6)-CBWW*FGK(1,2,5,6,3,4))
- FGK253=ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))
- WT=FGK135**2+(CCWW*FGK253)**2
- WTMAX=4.*D34*D56*(CAWW**2*DIGK(DT,DU)+CBWW**2*DIGK(DU,DT)-CAWW*
- & CBWW*DJGK(DT,DU)+CCWW**2*(DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU)))
-
- ELSEIF(ISUB.EQ.26.OR.ISUB.EQ.172.OR.ISUB.EQ.177) THEN
-C...Angular weight for f + f~' -> W+/- + H0 -> 2 quarks/leptons + H0
-C...(or H'0, or A0).
- WT=PKK(1,3)*PKK(2,4)
- WTMAX=(PKK(1,3)+PKK(1,4))*(PKK(2,3)+PKK(2,4))
-
- ELSEIF(ISUB.EQ.30.OR.ISUB.EQ.35) THEN
-C...Angular weight for f + g/gamma -> f + (gamma*/Z0)
-C...-> f + 2 quarks/leptons.
- CLILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
- & COUP(1,3)**2*HGZ(2,3)*COUP(3,3)**2
- CLIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,3)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
- & COUP(1,3)**2*HGZ(2,3)*COUP(3,4)**2
- CRILF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,3)/4.+
- & COUP(1,4)**2*HGZ(2,3)*COUP(3,3)**2
- CRIRF=COUP(1,1)**2*HGZ(2,1)*COUP(3,1)**2/16.+
- & COUP(1,1)*COUP(1,4)*HGZ(2,2)*COUP(3,1)*COUP(3,4)/4.+
- & COUP(1,4)**2*HGZ(2,3)*COUP(3,4)**2
- IF(K(ILIN(1),2).GT.0) WT=(CLILF+CRIRF)*(PKK(1,4)**2+
- & PKK(3,5)**2)+(CLIRF+CRILF)*(PKK(1,3)**2+PKK(4,5)**2)
- IF(K(ILIN(1),2).LT.0) WT=(CLILF+CRIRF)*(PKK(1,3)**2+
- & PKK(4,5)**2)+(CLIRF+CRILF)*(PKK(1,4)**2+PKK(3,5)**2)
- WTMAX=(CLILF+CLIRF+CRILF+CRIRF)*
- & ((PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2)
-
- ELSEIF(ISUB.EQ.31) THEN
-C...Angular weight for f + g -> f' + W+/- -> f' + 2 quarks/leptons.
- IF(K(ILIN(1),2).GT.0) WT=PKK(1,4)**2+PKK(3,5)**2
- IF(K(ILIN(1),2).LT.0) WT=PKK(1,3)**2+PKK(4,5)**2
- WTMAX=(PKK(1,3)+PKK(1,4))**2+(PKK(3,5)+PKK(4,5))**2
-
- ELSEIF(ISUB.EQ.71.OR.ISUB.EQ.72.OR.ISUB.EQ.73.OR.ISUB.EQ.76.OR.
- &ISUB.EQ.77) THEN
-C...Angular weight for V_L1 + V_L2 -> V_L3 + V_L4 (V = Z/W).
- WT=16.*PKK(3,5)*PKK(4,6)
- WTMAX=SH**2
-
- ELSEIF(ISUB.EQ.110) THEN
-C...Angular weight for f + f~ -> gamma + H0 -> gamma + X is isotropic.
- WT=1.
- WTMAX=1.
-
- ELSEIF(ISUB.EQ.141) THEN
- IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
-C...Angular weight for f + f~ -> gamma*/Z0/Z'0 -> 2 quarks/leptons.
-C...Couplings of incoming flavour.
- KFAI=IABS(MINT(15))
- EI=KCHG(KFAI,1)/3.
- AI=SIGN(1.,EI+0.1)
- VI=AI-4.*EI*XWV
- KFAIC=1
- IF(KFAI.LE.10.AND.MOD(KFAI,2).EQ.0) KFAIC=2
- IF(KFAI.GT.10.AND.MOD(KFAI,2).NE.0) KFAIC=3
- IF(KFAI.GT.10.AND.MOD(KFAI,2).EQ.0) KFAIC=4
- VPI=PARU(119+2*KFAIC)
- API=PARU(120+2*KFAIC)
-C...Couplings of final flavour.
- KFAF=IABS(KFL1(1))
- EF=KCHG(KFAF,1)/3.
- AF=SIGN(1.,EF+0.1)
- VF=AF-4.*EF*XWV
- KFAFC=1
- IF(KFAF.LE.10.AND.MOD(KFAF,2).EQ.0) KFAFC=2
- IF(KFAF.GT.10.AND.MOD(KFAF,2).NE.0) KFAFC=3
- IF(KFAF.GT.10.AND.MOD(KFAF,2).EQ.0) KFAFC=4
- VPF=PARU(119+2*KFAFC)
- APF=PARU(120+2*KFAFC)
-C...Asymmetry and weight.
- ASYM=2.*(EI*AI*VINT(112)*EF*AF+EI*API*VINT(113)*EF*APF+
- & 4.*VI*AI*VINT(114)*VF*AF+(VI*API+VPI*AI)*VINT(115)*
- & (VF*APF+VPF*AF)+4.*VPI*API*VINT(116)*VPF*APF)/
- & (EI**2*VINT(111)*EF**2+EI*VI*VINT(112)*EF*VF+
- & EI*VPI*VINT(113)*EF*VPF+(VI**2+AI**2)*VINT(114)*
- & (VF**2+AF**2)+(VI*VPI+AI*API)*VINT(115)*(VF*VPF+AF*APF)+
- & (VPI**2+API**2)*VINT(116)*(VPF**2+APF**2))
- WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
- WTMAX=2.+ABS(ASYM)
- ELSEIF(IP.EQ.1.AND.IABS(KFL1(1)).EQ.24) THEN
-C...Angular weight for f + f~ -> Z' -> W+ + W-.
- RM1=P(NSD(1)+1,5)**2/SH
- RM2=P(NSD(1)+2,5)**2/SH
- CCOS2=-(1./16.)*((1.-RM1-RM2)**2-4.*RM1*RM2)*
- & (1.-2.*RM1-2.*RM2+RM1**2+RM2**2+10.*RM1*RM2)
- CFLAT=-CCOS2+0.5*(RM1+RM2)*(1.-2.*RM1-2.*RM2+(RM2-RM1)**2)
- WT=CFLAT+CCOS2*CTHE(1)**2
- WTMAX=CFLAT+MAX(0.,CCOS2)
- ELSEIF(IP.EQ.1.AND.(KFL1(1).EQ.25.OR.KFL1(1).EQ.35.OR.
- & IABS(KFL1(1)).EQ.37)) THEN
-C...Angular weight for f + f~ -> Z' -> H0 + A0, H'0 + A0, H+ + H-.
- WT=1.-CTHE(1)**2
- WTMAX=1.
- ELSEIF(IP.EQ.1.AND.KDCY(1).EQ.3) THEN
-C...Angular weight for f + f~ -> Z' -> Z0 + H0.
- RM1=P(NSD(1)+1,5)**2/SH
- RM2=P(NSD(1)+2,5)**2/SH
- FLAM2=MAX(0.,(1.-RM1-RM2)**2-4.*RM1*RM2)
- WT=1.+FLAM2*(1.-CTHE(1)**2)/(8.*RM1)
- WTMAX=1.+FLAM2/(8.*RM1)
- ELSEIF(MZPWP.EQ.0) THEN
-C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
-C...(W:s like if intermediate Z).
- D34=P(IREF(IP,IORD),5)**2
- D56=P(IREF(IP,3-IORD),5)**2
- DT=PKK(1,3)+PKK(1,4)+D34
- DU=PKK(1,5)+PKK(1,6)+D56
- FGK135=ABS(FGK(1,2,3,4,5,6)-FGK(1,2,5,6,3,4))
- FGK253=ABS(FGK(2,1,5,6,3,4)-FGK(2,1,3,4,5,6))
- WT=(COUP(1,3)*FGK135)**2+(COUP(1,4)*FGK253)**2
- WTMAX=4.*D34*D56*(COUP(1,3)**2+COUP(1,4)**2)*
- & (DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU))
- ELSEIF(MZPWP.EQ.1) THEN
-C...Angular weight for f + f~ -> Z' -> W+ + W- -> 4 quarks/leptons
-C...(W:s approximately longitudinal, like if intermediate H).
- WT=16.*PKK(3,5)*PKK(4,6)
- WTMAX=SH**2
- ELSE
-C...Angular weight for f + f~ -> Z' -> H+ + H-, Z0 + H0, H0 + A0,
-C...H'0 + A0 -> 4 quarks/leptons.
- WT=1.
- WTMAX=1.
- ENDIF
-
- ELSEIF(ISUB.EQ.142) THEN
- IF(IP.EQ.1.AND.(KDCY(1).EQ.1.OR.KDCY(1).EQ.2)) THEN
-C...Angular weight for f + f~' -> W'+/- -> 2 quarks/leptons.
- KFAI=IABS(MINT(15))
- KFAIC=1
- IF(KFAI.GT.10) KFAIC=2
- VI=PARU(129+2*KFAIC)
- AI=PARU(130+2*KFAIC)
- KFAF=IABS(KFL1(1))
- KFAFC=1
- IF(KFAF.GT.10) KFAFC=2
- VF=PARU(129+2*KFAFC)
- AF=PARU(130+2*KFAFC)
- ASYM=8.*VI*AI*VF*AF/((VI**2+AI**2)*(VF**2+AF**2))
- WT=1.+ASYM*CTHE(1)*ISIGN(1,MINT(15)*KFL1(1))+CTHE(1)**2
- WTMAX=2.+ABS(ASYM)
- ELSEIF(IP.EQ.1.AND.IABS(KFL2(1)).EQ.23) THEN
-C...Angular weight for f + f~' -> W'+/- -> W+/- + Z0.
- RM1=P(NSD(1)+1,5)**2/SH
- RM2=P(NSD(1)+2,5)**2/SH
- CCOS2=-(1./16.)*((1.-RM1-RM2)**2-4.*RM1*RM2)*
- & (1.-2.*RM1-2.*RM2+RM1**2+RM2**2+10.*RM1*RM2)
- CFLAT=-CCOS2+0.5*(RM1+RM2)*(1.-2.*RM1-2.*RM2+(RM2-RM1)**2)
- WT=CFLAT+CCOS2*CTHE(1)**2
- WTMAX=CFLAT+MAX(0.,CCOS2)
- ELSEIF(IP.EQ.1.AND.KDCY(1).EQ.3) THEN
-C...Angular weight for f + f~ -> W'+/- -> W+/- + H0.
- RM1=P(NSD(1)+1,5)**2/SH
- RM2=P(NSD(1)+2,5)**2/SH
- FLAM2=MAX(0.,(1.-RM1-RM2)**2-4.*RM1*RM2)
- WT=1.+FLAM2*(1.-CTHE(1)**2)/(8.*RM1)
- WTMAX=1.+FLAM2/(8.*RM1)
- ELSEIF(MZPWP.EQ.0) THEN
-C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
-C...(W/Z like if intermediate W).
- D34=P(IREF(IP,IORD),5)**2
- D56=P(IREF(IP,3-IORD),5)**2
- DT=PKK(1,3)+PKK(1,4)+D34
- DU=PKK(1,5)+PKK(1,6)+D56
- FGK135=ABS(FGK(1,2,3,4,5,6)-FGK(1,2,5,6,3,4))
- FGK136=ABS(FGK(1,2,3,4,6,5)-FGK(1,2,6,5,3,4))
- WT=(COUP(5,3)*FGK135)**2+(COUP(5,4)*FGK136)**2
- WTMAX=4.*D34*D56*(COUP(5,3)**2+COUP(5,4)**2)*
- & (DIGK(DT,DU)+DIGK(DU,DT)-DJGK(DT,DU))
- ELSEIF(MZPWP.EQ.1) THEN
-C...Angular weight for f + f~' -> W' -> W + Z0 -> 4 quarks/leptons
-C...(W/Z approximately longitudinal, like if intermediate H).
- WT=16.*PKK(3,5)*PKK(4,6)
- WTMAX=SH**2
- ELSE
-C...Angular weight for f + f~ -> W' -> W + H0 -> whatever.
- WT=1.
- WTMAX=1.
- ENDIF
-
- ELSEIF(ISUB.EQ.145.OR.ISUB.EQ.162.OR.ISUB.EQ.163.OR.ISUB.EQ.164)
- &THEN
-C...Isotropic decay of leptoquarks (assumed spin 0).
- WT=1.
- WTMAX=1.
-
- ELSEIF(ISUB.EQ.147.OR.ISUB.EQ.148) THEN
-C...Decays of (spin 1/2) q* -> q + (g,gamma) or (Z0,W+-).
- SIDE=1.
- IF(MINT(16).EQ.21) SIDE=-1.
- IF(IP.EQ.1.AND.(KFL1(1).EQ.21.OR.KFL1(1).EQ.22)) THEN
- WT=1.+SIDE*CTHE(1)
- WTMAX=2.
- ELSEIF(IP.EQ.1) THEN
- RM1=P(NSD(1)+1,5)**2/SH
- WT=1.+SIDE*CTHE(1)*(1.-0.5*RM1)/(1.+0.5*RM1)
- WTMAX=1.+(1.-0.5*RM1)/(1.+0.5*RM1)
- ELSE
-C...W/Z decay assumed isotropic, since not known.
- WT=1.
- WTMAX=1.
- ENDIF
-
- ELSEIF(ISUB.EQ.149) THEN
-C...Isotropic decay of techni-eta.
- WT=1.
- WTMAX=1.
-
-C...Obtain correct angular distribution by rejection techniques.
- ELSE
- WT=1.
- WTMAX=1.
- ENDIF
- IF(WT.LT.RLU(0)*WTMAX) GOTO 240
-
-C...Construct massive four-vectors using angles chosen. Mark decayed
-C...resonances, add documentation lines. Shower evolution.
- 400 DO 470 JT=1,JTMAX
- IF(KDCY(JT).EQ.0) GOTO 470
- ID=IREF(IP,JT)
- IF(ISET(ISUB).NE.6.OR.JT.NE.3) THEN
- DO 410 J=1,5
- DPMO(J)=P(ID,J)
- 410 CONTINUE
- DPMO(4)=SQRT(DPMO(1)**2+DPMO(2)**2+DPMO(3)**2+DPMO(5)**2)
- CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,ACOS(CTHE(JT)),PHI(JT),
- & DPMO(1)/DPMO(4),DPMO(2)/DPMO(4),DPMO(3)/DPMO(4))
- ELSE
-C...Z + q + q~ : angles already known, in rest frame of system.
- DO 420 J=1,3
- DBEZQQ(J)=(P(ID,J)+P(ID+1,J)+P(ID+2,J))/(P(ID,4)+P(ID+1,4)+
- & P(ID+2,4))
- 420 CONTINUE
- K(N+1,1)=1
- DO 430 J=1,5
- P(N+1,J)=P(ID,J)
- 430 CONTINUE
- CALL LUDBRB(N+1,N+1,0.,0.,-DBEZQQ(1),-DBEZQQ(2),-DBEZQQ(3))
- PHIZQQ=ULANGL(P(N+1,1),P(N+1,2))
- THEZQQ=ULANGL(P(N+1,3),SQRT(P(N+1,1)**2+P(N+1,2)**2))
- CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,ACOS(VINT(81)),VINT(82),
- & 0D0,0D0,DBLE(SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)/
- & P(N+1,4)))
- CALL LUDBRB(NSD(JT)+1,NSD(JT)+2,THEZQQ,PHIZQQ,DBEZQQ(1),
- & DBEZQQ(2),DBEZQQ(3))
- ENDIF
- K(ID,1)=K(ID,1)+10
- KFA=IABS(K(ID,2))
- IF(KFA.EQ.38.OR.KFA.EQ.39.OR.((MSTP(6).EQ.1.OR.MSTP(49).GE.1).AND.
- &(KFA.EQ.7.OR.KFA.EQ.8)).OR.(MSTP(48).GE.1.AND.KFA.EQ.6)) THEN
-C...Do not kill colour flow through techni-eta, t, leptoquark or q*!
- ELSE
- K(ID,4)=NSD(JT)+1
- K(ID,5)=NSD(JT)+2
- ENDIF
- IDOC=MINT(83)+MINT(4)
- DO 450 I=NSD(JT)+1,NSD(JT)+2
- I1=MINT(83)+MINT(4)+1
- K(I,3)=I1
- IF(MSTP(128).GE.1) K(I,3)=ID
- IF(MSTP(128).LE.1.AND.MINT(4).LT.MSTP(126)) THEN
- MINT(4)=MINT(4)+1
- K(I1,1)=21
- K(I1,2)=K(I,2)
- K(I1,3)=IREF(IP,JT+3)
- DO 440 J=1,5
- P(I1,J)=P(I,J)
- 440 CONTINUE
- ENDIF
- 450 CONTINUE
-C...Shower - top currently special case.
- NSHBEF=N
- IF(MSTP(71).GE.1.AND.(KDCY(JT).LE.2.OR.KFA.EQ.6.OR.KFA.EQ.7.OR.
- &KFA.EQ.8)) CALL LUSHOW(NSD(JT)+1,NSD(JT)+2,P(ID,5))
- NSHAFT=N
-
-C...Check if new resonances were produced.
- KNSDA1=IABS(K(NSD(JT)+1,2))
- KNSDA2=IABS(K(NSD(JT)+2,2))
- IF(KNSDA1.EQ.6.OR.KNSDA2.EQ.6.OR.KNSDA1.EQ.7.OR.KNSDA2.EQ.7.OR.
- &KNSDA1.EQ.8.OR.KNSDA2.EQ.8) THEN
- NSD1=0
- NSD2=0
- DO 460 I=NSD(JT)+1,N
- IF(K(I,1).LT.10.AND.K(I,2).EQ.K(NSD(JT)+1,2)) NSD1=I
- IF(K(I,1).LT.10.AND.K(I,2).EQ.K(NSD(JT)+2,2)) NSD2=I
- 460 CONTINUE
- IF(NSD1.NE.0.AND.NSD2.NE.0) THEN
- NP=NP+1
- IREF(NP,1)=NSD1
- IREF(NP,2)=NSD2
- IREF(NP,3)=0
- IREF(NP,4)=IDOC+1
- IREF(NP,5)=IDOC+2
- IREF(NP,6)=0
- IREF(NP,7)=K(IREF(IP,JT),2)
- IREF(NP,8)=IREF(IP,JT)
- ENDIF
- ELSEIF(KDCY(JT).EQ.3.OR.KFA.EQ.6.OR.KFA.EQ.7.OR.KFA.EQ.8) THEN
- NP=NP+1
- IREF(NP,1)=NSD(JT)+1
- IREF(NP,2)=NSD(JT)+2
- IF(NSHAFT-NSHBEF.GT.0) THEN
- IREF(NP,1)=NSHBEF+2
- IREF(NP,2)=NSHBEF+3
- ENDIF
- IREF(NP,3)=0
- IREF(NP,4)=IDOC+1
- IREF(NP,5)=IDOC+2
- IREF(NP,6)=0
- IREF(NP,7)=K(IREF(IP,JT),2)
- IREF(NP,8)=IREF(IP,JT)
- ENDIF
- 470 CONTINUE
-
-C...Fill information for 2 -> 1 -> 2. Loop back if needed.
- IF(JTMAX.EQ.1.AND.KDCY(1).NE.0) THEN
- MINT(7)=MINT(83)+6+2*ISET(ISUB)
- MINT(8)=MINT(83)+7+2*ISET(ISUB)
- MINT(25)=KFL1(1)
- MINT(26)=KFL2(1)
- VINT(23)=CTHE(1)
- RM3=P(N-1,5)**2/SH
- RM4=P(N,5)**2/SH
- BE34=SQRT(MAX(0.,(1.-RM3-RM4)**2-4.*RM3*RM4))
- VINT(45)=-0.5*SH*(1.-RM3-RM4-BE34*CTHE(1))
- VINT(46)=-0.5*SH*(1.-RM3-RM4+BE34*CTHE(1))
- VINT(48)=0.25*SH*BE34**2*MAX(0.,1.-CTHE(1)**2)
- VINT(47)=SQRT(VINT(48))
- ENDIF
- 480 IF(IP.LT.NP) GOTO 130
-
- RETURN
- END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff