1 #include "isajet/pilot.h"
4 C Finish generation of whiggs events started bY TWOJET.
5 C Select W decay modes as allowed by WMODE1, WMODE2.
6 C Generate W decay angles and put vectors in PPAIR.
9 #if defined(CERNLIB_IMPNONE)
12 #include "isajet/itapes.inc"
13 #include "isajet/qcdpar.inc"
14 #include "isajet/jetpar.inc"
15 #include "isajet/primar.inc"
16 #include "isajet/q1q2.inc"
17 #include "isajet/jetsig.inc"
18 #include "isajet/const.inc"
19 #include "isajet/qsave.inc"
20 #include "isajet/wcon.inc"
21 #include "isajet/pjets.inc"
22 #include "isajet/pinits.inc"
23 #include "isajet/keys.inc"
24 #include "isajet/hcon.inc"
25 #include "isajet/wwpar.inc"
26 #include "isajet/xmssm.inc"
28 DIMENSION X(2),LIST(25),P1WCM(4),P2WCM(4),P1LAB(4),P2LAB(4)
32 REAL GVQ(2),GAQ(2),GVL(2),GAL(2)
33 REAL X,RND,RANF,CBRWW,AMASS,AM0,AM1,AM2,
34 $E1CM,E2CM,P12CM,CTHCM,STHCM,PHICM,CPHICM,SPHICM,P1WCM,P2WCM,
35 $PBOOST,P1LAB,P2LAB,ZHSIG,ZHMAX
36 INTEGER JWWTYP,JET,JWT,JQ,IQ1,IQ2,LIST,NREJ,NJ0,K
37 REAL BRANCH(12),SUMBR,BETAWH,GAMWH,PZWHCM,CTHD,WHSIG
38 INTEGER IDABS,IDABS1,IDABSJ,IDIABS
40 DATA LIST/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
41 $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16/
52 IF(KEYS(10).AND..NOT.GOMSSM) THEN
53 JWWTYP(1)=JETTYP(1)-25
54 JWWTYP(2)=JETTYP(2)-25
55 ELSEIF(KEYS(10).AND.GOMSSM) THEN
56 JWWTYP(1)=JETTYP(1)-76
57 JWWTYP(2)=JETTYP(2)-76
60 C Select W decay modes and put in /JETSET/. First particle
61 C is always the fermion.
64 IDABS=IABS(IDJETS(JET))
65 IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 200
68 C Must only consider allowed decays for this mass
79 IF(AM1+AM2.LT.PJETS(5,JET)) THEN
80 BRANCH(JQ)=RBRWW(JQ,JWT,JET)
81 SUMBR=SUMBR+BRANCH(JQ)
86 IF(SUMBR.LE.0.) GO TO 998
88 202 BRANCH(JQ)=BRANCH(JQ)/SUMBR
92 CBRWW=CBRWW+BRANCH(JQ)
93 IF(RND.GT.CBRWW) GO TO 210
96 IDPAIR(NPAIR+1)=LIST(IQ1)
97 IDPAIR(NPAIR+2)=LIST(IQ2)
98 PPAIR(5,NPAIR+1)=AMASS(LIST(IQ1))
99 PPAIR(5,NPAIR+2)=AMASS(LIST(IQ2))
108 C Generate decay uniformly in angle and put in PPAIR.
109 C Will check cross section later.
114 IDABS=IABS(IDJETS(JET))
115 IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 310
116 C Construct W com momenta.
121 E1CM=(AM0**2+AM1**2-AM2**2)/(2.*AM0)
122 E2CM=(AM0**2+AM2**2-AM1**2)/(2.*AM0)
123 P12CM=(AM0**2-AM1**2-AM2**2)**2-4.*(AM1*AM2)**2
124 P12CM=SQRT(P12CM)/(2.*AM0)
126 STHCM=SQRT(1.-CTHCM**2)
130 P1WCM(1)=P12CM*STHCM*CPHICM
132 P1WCM(2)=P12CM*STHCM*SPHICM
138 C Boost to lab frame.
140 320 PBOOST(K)=-PJETS(K,JET)
141 PBOOST(4)=PJETS(4,JET)
142 CALL LBOOST(PBOOST,1,P1WCM,P1LAB)
143 CALL LBOOST(PBOOST,1,P2WCM,P2LAB)
145 PPAIR(K,NJ0-1)=P1LAB(K)
146 PPAIR(K,NJ0)=P2LAB(K)
151 C Impose simple (1+-cos(theta))**2 decay distribution for WH
152 C Must use P1 in WH CoM frame
153 IF (IDABSJ.NE.80.AND.IDABSJ.NE.90) GO TO 400
154 BETAWH=(PJETS(3,1)+PJETS(3,2))/(PJETS(4,1)+PJETS(4,2))
155 GAMWH=1./SQRT(1.-BETAWH**2)
156 PZWHCM=GAMWH*(P1LAB(3)-BETAWH*P1LAB(4))
157 CTHD=PZWHCM/SQRT(P1LAB(1)**2+P1LAB(2)**2+PZWHCM**2)
158 IF (IDINIT(1).LT.0) CTHD=-CTHD
159 IDIABS=IABS(IDINIT(1))
160 IDABS1=IABS(IDPAIR(1))
161 IF (IDABSJ.EQ.80) THEN
163 IF(WHSIG.GT.4*RANF()) GO TO 400
165 IF (IDABSJ.EQ.90) THEN
166 IF (IDIABS.EQ.1.OR.IDIABS.EQ.4) THEN
167 IF (IDABS1.EQ.1.OR.IDABS1.EQ.4) THEN
168 ZHSIG=(GVQ(1)**2+GAQ(1)**2)**2*(1.+CTHD**2)
169 $ +8*GVQ(1)*GAQ(1)*GVQ(1)*GAQ(1)*CTHD
170 ZHMAX=2*(GVQ(1)**2+GAQ(1)**2)**2
171 $ +8*GVQ(1)*GAQ(1)*GVQ(1)*GAQ(1)
172 ELSEIF (IDABS1.EQ.2.OR.IDABS1.EQ.3.OR.IDABS1.EQ.5) THEN
173 ZHSIG=(GVQ(1)**2+GAQ(1))*(GVQ(2)**2+GAQ(2)**2)*(1.+CTHD**2)
174 $ +8*GVQ(1)*GAQ(1)*GVQ(2)*GAQ(2)*CTHD
175 ZHMAX=(GVQ(1)**2+GAQ(1))*(GVQ(2)**2+GAQ(2)**2)*2
176 $ +8*GVQ(1)*GAQ(1)*GVQ(2)*GAQ(2)
177 ELSEIF (IDABS1.EQ.11.OR.IDABS1.EQ.13.OR.IDABS1.EQ.15) THEN
178 ZHSIG=(GVQ(1)**2+GAQ(1))*(GVL(1)**2+GAL(1)**2)*(1.+CTHD**2)
179 $ +8*GVQ(1)*GAQ(1)*GVL(1)*GAL(1)*CTHD
180 ZHMAX=(GVQ(1)**2+GAQ(1))*(GVL(1)**2+GAL(1)**2)*2
181 $ +8*GVQ(1)*GAQ(1)*GVL(1)*GAL(1)
182 ELSEIF (IDABS1.EQ.12.OR.IDABS1.EQ.14.OR.IDABS1.EQ.16) THEN
183 ZHSIG=(GVQ(1)**2+GAQ(1))*(GVL(2)**2+GAL(2)**2)*(1.+CTHD**2)
184 $ +8*GVQ(1)*GAQ(1)*GVL(2)*GAL(2)*CTHD
185 ZHMAX=(GVQ(1)**2+GAQ(1))*(GVL(2)**2+GAL(2)**2)*2
186 $ +8*GVQ(1)*GAQ(1)*GVL(2)*GAL(2)
188 ELSE IF (IDIABS.EQ.2.OR.IDIABS.EQ.3.OR.IDIABS.EQ.5) THEN
189 IF (IDABS1.EQ.1.OR.IDABS1.EQ.4) THEN
190 ZHSIG=(GVQ(2)**2+GAQ(2)**2)**2*(1.+CTHD**2)
191 $ +8*GVQ(2)*GAQ(2)*GVQ(1)*GAQ(1)*CTHD
192 ZHMAX=(GVQ(2)**2+GAQ(2)**2)**2*2
193 $ +8*GVQ(2)*GAQ(2)*GVQ(1)*GAQ(1)
194 ELSEIF (IDABS1.EQ.2.OR.IDABS1.EQ.3.OR.IDABS1.EQ.5) THEN
195 ZHSIG=(GVQ(2)**2+GAQ(2))*(GVQ(2)**2+GAQ(2)**2)*(1.+CTHD**2)
196 $ +8*GVQ(2)*GAQ(2)*GVQ(2)*GAQ(2)*CTHD
197 ZHMAX=(GVQ(2)**2+GAQ(2))*(GVQ(2)**2+GAQ(2)**2)*2
198 $ +8*GVQ(2)*GAQ(2)*GVQ(2)*GAQ(2)
199 ELSEIF (IDABS1.EQ.11.OR.IDABS1.EQ.13.OR.IDABS1.EQ.15) THEN
200 ZHSIG=(GVQ(2)**2+GAQ(2))*(GVL(1)**2+GAL(1)**2)*(1.+CTHD**2)
201 $ +8*GVQ(2)*GAQ(2)*GVL(1)*GAL(1)*CTHD
202 ZHMAX=(GVQ(2)**2+GAQ(2))*(GVL(1)**2+GAL(1)**2)*2
203 $ +8*GVQ(2)*GAQ(2)*GVL(1)*GAL(1)
204 ELSEIF (IDABS1.EQ.12.OR.IDABS1.EQ.14.OR.IDABS1.EQ.16) THEN
205 ZHSIG=(GVQ(2)**2+GAQ(2))*(GVL(2)**2+GAL(2)**2)*(1.+CTHD**2)
206 $ +8*GVQ(2)*GAQ(2)*GVL(2)*GAL(2)*CTHD
207 ZHMAX=(GVQ(2)**2+GAQ(2))*(GVL(2)**2+GAL(2)**2)*2
208 $ +8*GVQ(2)*GAQ(2)*GVL(2)*GAL(2)
211 IF(ZHSIG.GT.RANF()*ZHMAX) GO TO 400
214 IF(NREJ.LT.NTRIES) GO TO 300
223 WRITE(ITLIS,9991) NREJ
224 9991 FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE ',
225 1'A GOOD WHIGGS EVENT.'/' CHECK LIMITS OR INCREASE NTRIES.')
228 WRITE(ITLIS,9981) JET
229 9981 FORMAT(//' ERROR IN WHIGGS ... NO DECAY POSSIBLE FOR JET',I3)