1 #include "isajet/pilot.h"
4 C Finish generation of wpair 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.
8 C Also generate massless decay vectors PZERO for matrix
9 C element -- double precision for 32-bit machines.
11 C Ver 6.26: Check kinematics for W -> ff decay, since Z0 from
12 C Higgs decay can be virtual.
13 C Ver. 6.30: Added check in loop 201.
14 C Ver. 7.14: Add MSSM Higgs hooks
16 #if defined(CERNLIB_IMPNONE)
19 #include "isajet/itapes.inc"
20 #include "isajet/qcdpar.inc"
21 #include "isajet/jetpar.inc"
22 #include "isajet/primar.inc"
23 #include "isajet/q1q2.inc"
24 #include "isajet/jetsig.inc"
25 #include "isajet/wwsig.inc"
26 #include "isajet/wwpar.inc"
27 #include "isajet/const.inc"
28 #include "isajet/qsave.inc"
29 #include "isajet/wcon.inc"
30 #include "isajet/pjets.inc"
31 #include "isajet/pinits.inc"
32 #include "isajet/keys.inc"
33 #include "isajet/wsig.inc"
34 #include "isajet/hcon.inc"
35 #include "isajet/xmssm.inc"
37 DIMENSION X(2),LIST(25),P1WCM(4),P2WCM(4),P1LAB(4),P2LAB(4)
38 $,P1CM0(4),P2CM0(4),P1LAB0(4),P2LAB0(4)
42 EQUIVALENCE (PWW(1,1),P3WW(1))
43 DIMENSION JWWTYP(2),THWFF(2),PHIWFF(2)
44 #if defined(CERNLIB_SINGLE)
45 REAL P1CM0,P2CM0,DPHI,DCTH,DSTH,DAM0,PWW,BOOST
47 #if defined(CERNLIB_DOUBLE)
48 DOUBLE PRECISION P1CM0,P2CM0,DPHI,DCTH,DSTH,DAM0,PWW,BOOST
50 REAL AMWW1,AMWW2,X,STRUC,STRUCW,RND,RANF,CBRWW,AMASS,AM0,AM1,AM2,
51 $E1CM,E2CM,P12CM,CTHCM,STHCM,PHICM,CPHICM,SPHICM,P1WCM,P2WCM,
52 $PBOOST,P1LAB,P2LAB,AFX,SGWWMX,P1LAB0,P2LAB0,THWFF,PHIWFF
53 INTEGER IH,IQ,JWWTYP,JET,JWT,JQ,IQ1,IQ2,LIST,NREJ,NJ0,K
55 INTEGER IDABS,IDABS1,IDABS2
57 DATA LIST/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
58 $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16/
60 C Initialize for given W type.
63 CALL WWKIN(AMWW1,AMWW2)
66 C Calculate and save structure functions.
69 121 QSAVE(IQ,IH)=STRUC(X(IH),QSQ,IQ,IDIN(IH))/X(IH)
72 IF(KEYS(7).OR.KEYS(9)) THEN
74 123 QSAVE(IQ,IH)=STRUCW(X(IH),IQ-25,IDIN(IH))/X(IH)
77 C JWWTYP points to W types 1,2,3,4
81 ELSEIF((KEYS(7).AND..NOT.GOMSSM).OR.KEYS(9)) THEN
82 JWWTYP(1)=JETTYP(1)-25
83 JWWTYP(2)=JETTYP(2)-25
84 ELSEIF(KEYS(7).AND.GOMSSM) THEN
85 JWWTYP(1)=JETTYP(1)-76
86 JWWTYP(2)=JETTYP(2)-76
89 C Select W decay modes and put in /JETSET/. First particle
90 C is always the fermion.
93 IDABS=IABS(IDJETS(JET))
94 IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 200
97 C Must only consider allowed decays for this mass
108 IF(AM1+AM2.LT.PJETS(5,JET)) THEN
109 BRANCH(JQ)=RBRWW(JQ,JWT,JET)
110 SUMBR=SUMBR+BRANCH(JQ)
115 IF(SUMBR.LE.0.) GO TO 998
117 202 BRANCH(JQ)=BRANCH(JQ)/SUMBR
121 CBRWW=CBRWW+BRANCH(JQ)
122 IF(RND.GT.CBRWW) GO TO 210
125 IDPAIR(NPAIR+1)=LIST(IQ1)
126 IDPAIR(NPAIR+2)=LIST(IQ2)
127 PPAIR(5,NPAIR+1)=AMASS(LIST(IQ1))
128 PPAIR(5,NPAIR+2)=AMASS(LIST(IQ2))
137 C Generate decay uniformly in angle and put in PPAIR.
138 C Will check cross section later.
143 IDABS=IABS(IDJETS(JET))
144 IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 310
145 C Construct W com momenta.
149 E1CM=(AM0**2+AM1**2-AM2**2)/(2.*AM0)
150 E2CM=(AM0**2+AM2**2-AM1**2)/(2.*AM0)
151 P12CM=(AM0**2-AM1**2-AM2**2)**2-4.*(AM1*AM2)**2
152 P12CM=SQRT(P12CM)/(2.*AM0)
154 STHCM=SQRT(1.-CTHCM**2)
158 P1WCM(1)=P12CM*STHCM*CPHICM
160 P1WCM(2)=P12CM*STHCM*SPHICM
166 C Also construct zero mass vectors at same angle
167 #if defined(CERNLIB_SINGLE)
169 P1CM0(1)=.5*AM0*STHCM*CPHICM
171 P1CM0(2)=.5*AM0*STHCM*SPHICM
173 P1CM0(3)=.5*AM0*CTHCM
178 #if defined(CERNLIB_DOUBLE)
182 DSTH=DSQRT(1.D0-DCTH**2)
184 P1CM0(1)=.5*AM0*DSTH*DCOS(DPHI)
186 P1CM0(2)=.5*AM0*DSTH*DSIN(DPHI)
193 C Boost to lab frame.
195 320 PBOOST(K)=-PJETS(K,JET)
196 PBOOST(4)=PJETS(4,JET)
197 CALL LBOOST(PBOOST,1,P1WCM,P1LAB)
198 CALL LBOOST(PBOOST,1,P2WCM,P2LAB)
200 PPAIR(K,NJ0-1)=P1LAB(K)
201 PPAIR(K,NJ0)=P2LAB(K)
203 C Boost zero mass vectors -- double precision for 32 bits.
204 PZERO(4,NJ0-1)=(P1CM0(4)*PWW(4,JET)+P1CM0(1)*PWW(1,JET)
205 $ +P1CM0(2)*PWW(2,JET)+P1CM0(3)*PWW(3,JET))/PWW(5,JET)
206 BOOST=(P1CM0(4)+PZERO(4,NJ0-1))/(PWW(4,JET)+PWW(5,JET))
208 340 PZERO(K,NJ0-1)=P1CM0(K)+BOOST*PWW(K,JET)
209 PZERO(4,NJ0)=(P2CM0(4)*PWW(4,JET)+P2CM0(1)*PWW(1,JET)
210 $ +P2CM0(2)*PWW(2,JET)+P2CM0(3)*PWW(3,JET))/PWW(5,JET)
211 BOOST=(P2CM0(4)+PZERO(4,NJ0))/(PWW(4,JET)+PWW(5,JET))
213 350 PZERO(K,NJ0)=P2CM0(K)+BOOST*PWW(K,JET)
217 C Calculate cross section SIGWW2 containing TBRWW*RBRWW.
218 C Compare with WW cross section containing TBRWW. Ratio
219 C must be bounded by 3/(4*PI) for each W.
225 IF(IDJETS(1).NE.10) SGWWMX=SGWWMX*RBRWW(JQWW(1),JWWTYP(1),1)*AFX
226 IF(IDJETS(2).NE.10) SGWWMX=SGWWMX*RBRWW(JQWW(2),JWWTYP(2),2)*AFX
228 C Note that except for WW -> WW processes, SIGH3 just computes
229 C the decay angular distribution, so it can be used for both
230 C for SM and SUSY HL0/HH0 decays; HA0 -> WW is forbidden.
231 C For Z + HL0 decays, we just return, ie use phase space.
232 IDABS1=IABS(IDJETS(1))
233 IDABS2=IABS(IDJETS(2))
234 IF(.NOT.(IDABS1.EQ.10.OR.IDABS1.EQ.80.OR.IDABS1.EQ.90)) RETURN
235 IF(.NOT.(IDABS2.EQ.10.OR.IDABS2.EQ.80.OR.IDABS2.EQ.90)) RETURN
242 IF(WWSIG.GT.SGWWMX*RANF()) GO TO 400
244 IF(NREJ.LT.NTRIES) GO TO 300
253 WRITE(ITLIS,9991) NREJ
254 9991 FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE ',
255 1'A GOOD WPAIR EVENT.'/' CHECK LIMITS OR INCREASE NTRIES.')
258 WRITE(ITLIS,9981) JET
259 9981 FORMAT(//' ERROR IN WPAIR ... NO DECAY POSSIBLE FOR JET',I3)