]>
Commit | Line | Data |
---|---|---|
0795afa3 | 1 | #include "isajet/pilot.h" |
2 | SUBROUTINE WPAIR | |
3 | C | |
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. | |
7 | C | |
8 | C Also generate massless decay vectors PZERO for matrix | |
9 | C element -- double precision for 32-bit machines. | |
10 | C | |
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 | |
15 | C | |
16 | #if defined(CERNLIB_IMPNONE) | |
17 | IMPLICIT NONE | |
18 | #endif | |
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" | |
36 | C | |
37 | DIMENSION X(2),LIST(25),P1WCM(4),P2WCM(4),P1LAB(4),P2LAB(4) | |
38 | $,P1CM0(4),P2CM0(4),P1LAB0(4),P2LAB0(4) | |
39 | 1,PBOOST(4) | |
40 | EQUIVALENCE (X(1),X1) | |
41 | DIMENSION PWW(5,2) | |
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 | |
46 | #endif | |
47 | #if defined(CERNLIB_DOUBLE) | |
48 | DOUBLE PRECISION P1CM0,P2CM0,DPHI,DCTH,DSTH,DAM0,PWW,BOOST | |
49 | #endif | |
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 | |
54 | REAL BRANCH(12),SUMBR | |
55 | INTEGER IDABS,IDABS1,IDABS2 | |
56 | C | |
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/ | |
59 | C | |
60 | C Initialize for given W type. | |
61 | AMWW1=PJETS(5,1) | |
62 | AMWW2=PJETS(5,2) | |
63 | CALL WWKIN(AMWW1,AMWW2) | |
64 | NPAIR=0 | |
65 | C | |
66 | C Calculate and save structure functions. | |
67 | DO 120 IH=1,2 | |
68 | DO 121 IQ=1,13 | |
69 | 121 QSAVE(IQ,IH)=STRUC(X(IH),QSQ,IQ,IDIN(IH))/X(IH) | |
70 | DO 122 IQ=14,25 | |
71 | 122 QSAVE(IQ,IH)=0. | |
72 | IF(KEYS(7).OR.KEYS(9)) THEN | |
73 | DO 123 IQ=26,29 | |
74 | 123 QSAVE(IQ,IH)=STRUCW(X(IH),IQ-25,IDIN(IH))/X(IH) | |
75 | ENDIF | |
76 | 120 CONTINUE | |
77 | C JWWTYP points to W types 1,2,3,4 | |
78 | IF(KEYS(6)) THEN | |
79 | JWWTYP(1)=JETTYP(1) | |
80 | JWWTYP(2)=JETTYP(2) | |
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 | |
87 | ENDIF | |
88 | C | |
89 | C Select W decay modes and put in /JETSET/. First particle | |
90 | C is always the fermion. | |
91 | ||
92 | DO 200 JET=1,2 | |
93 | IDABS=IABS(IDJETS(JET)) | |
94 | IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 200 | |
95 | RND=RANF() | |
96 | JWT=JWWTYP(JET) | |
97 | C Must only consider allowed decays for this mass | |
98 | SUMBR=0. | |
99 | DO 201 JQ=1,12 | |
100 | IQ1=2*JQ | |
101 | IQ2=MATCH(IQ1,JWT) | |
102 | IF(IQ2.EQ.0) THEN | |
103 | BRANCH(JQ)=0. | |
104 | GO TO 201 | |
105 | ENDIF | |
106 | AM1=AMASS(LIST(IQ1)) | |
107 | AM2=AMASS(LIST(IQ2)) | |
108 | IF(AM1+AM2.LT.PJETS(5,JET)) THEN | |
109 | BRANCH(JQ)=RBRWW(JQ,JWT,JET) | |
110 | SUMBR=SUMBR+BRANCH(JQ) | |
111 | ELSE | |
112 | BRANCH(JQ)=0. | |
113 | ENDIF | |
114 | 201 CONTINUE | |
115 | IF(SUMBR.LE.0.) GO TO 998 | |
116 | DO 202 JQ=1,12 | |
117 | 202 BRANCH(JQ)=BRANCH(JQ)/SUMBR | |
118 | C | |
119 | CBRWW=0. | |
120 | DO 210 JQ=1,12 | |
121 | CBRWW=CBRWW+BRANCH(JQ) | |
122 | IF(RND.GT.CBRWW) GO TO 210 | |
123 | IQ1=2*JQ | |
124 | IQ2=MATCH(IQ1,JWT) | |
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)) | |
129 | JPAIR(NPAIR+1)=JET | |
130 | JPAIR(NPAIR+2)=JET | |
131 | NPAIR=NPAIR+2 | |
132 | JQWW(JET)=JQ | |
133 | GO TO 200 | |
134 | 210 CONTINUE | |
135 | 200 CONTINUE | |
136 | C | |
137 | C Generate decay uniformly in angle and put in PPAIR. | |
138 | C Will check cross section later. | |
139 | C | |
140 | NREJ=0 | |
141 | 300 NJ0=2 | |
142 | DO 310 JET=1,2 | |
143 | IDABS=IABS(IDJETS(JET)) | |
144 | IF(IDABS.NE.80.AND.IDABS.NE.90) GO TO 310 | |
145 | C Construct W com momenta. | |
146 | AM0=PJETS(5,JET) | |
147 | AM1=PPAIR(5,NJ0-1) | |
148 | AM2=PPAIR(5,NJ0) | |
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) | |
153 | CTHCM=2.*RANF()-1. | |
154 | STHCM=SQRT(1.-CTHCM**2) | |
155 | PHICM=2.*PI*RANF() | |
156 | CPHICM=COS(PHICM) | |
157 | SPHICM=SIN(PHICM) | |
158 | P1WCM(1)=P12CM*STHCM*CPHICM | |
159 | P2WCM(1)=-P1WCM(1) | |
160 | P1WCM(2)=P12CM*STHCM*SPHICM | |
161 | P2WCM(2)=-P1WCM(2) | |
162 | P1WCM(3)=P12CM*CTHCM | |
163 | P2WCM(3)=-P1WCM(3) | |
164 | P1WCM(4)=E1CM | |
165 | P2WCM(4)=E2CM | |
166 | C Also construct zero mass vectors at same angle | |
167 | #if defined(CERNLIB_SINGLE) | |
168 | C Single precision. | |
169 | P1CM0(1)=.5*AM0*STHCM*CPHICM | |
170 | P2CM0(1)=-P1CM0(1) | |
171 | P1CM0(2)=.5*AM0*STHCM*SPHICM | |
172 | P2CM0(2)=-P1CM0(2) | |
173 | P1CM0(3)=.5*AM0*CTHCM | |
174 | P2CM0(3)=-P1CM0(3) | |
175 | P1CM0(4)=.5*AM0 | |
176 | P2CM0(4)=P1CM0(4) | |
177 | #endif | |
178 | #if defined(CERNLIB_DOUBLE) | |
179 | C Double precision. | |
180 | DAM0=AM0 | |
181 | DCTH=CTHCM | |
182 | DSTH=DSQRT(1.D0-DCTH**2) | |
183 | DPHI=PHICM | |
184 | P1CM0(1)=.5*AM0*DSTH*DCOS(DPHI) | |
185 | P2CM0(1)=-P1CM0(1) | |
186 | P1CM0(2)=.5*AM0*DSTH*DSIN(DPHI) | |
187 | P2CM0(2)=-P1CM0(2) | |
188 | P1CM0(3)=.5*AM0*DCTH | |
189 | P2CM0(3)=-P1CM0(3) | |
190 | P1CM0(4)=.5*AM0 | |
191 | P2CM0(4)=P1CM0(4) | |
192 | #endif | |
193 | C Boost to lab frame. | |
194 | DO 320 K=1,3 | |
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) | |
199 | DO 330 K=1,4 | |
200 | PPAIR(K,NJ0-1)=P1LAB(K) | |
201 | PPAIR(K,NJ0)=P2LAB(K) | |
202 | 330 CONTINUE | |
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)) | |
207 | DO 340 K=1,3 | |
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)) | |
212 | DO 350 K=1,3 | |
213 | 350 PZERO(K,NJ0)=P2CM0(K)+BOOST*PWW(K,JET) | |
214 | NJ0=NJ0+2 | |
215 | 310 CONTINUE | |
216 | C | |
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. | |
220 | C | |
221 | AFX=3./(2.*PI) | |
222 | IF(KEYS(6)) THEN | |
223 | CALL SIGWW2 | |
224 | SGWWMX=SIGEVT | |
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 | |
227 | ELSEIF(KEYS(7)) THEN | |
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 | |
236 | CALL SIGH3 | |
237 | SGWWMX=SIGLLQ*AFX**2 | |
238 | ELSEIF(KEYS(9)) THEN | |
239 | CALL SIGTC3 | |
240 | SGWWMX=SIGLLQ*AFX**2 | |
241 | ENDIF | |
242 | IF(WWSIG.GT.SGWWMX*RANF()) GO TO 400 | |
243 | NREJ=NREJ+1 | |
244 | IF(NREJ.LT.NTRIES) GO TO 300 | |
245 | GO TO 999 | |
246 | C | |
247 | C Good event | |
248 | C | |
249 | 400 CONTINUE | |
250 | RETURN | |
251 | C | |
252 | 999 CALL PRTEVT(0) | |
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.') | |
256 | STOP 99 | |
257 | 998 CALL PRTEVT(0) | |
258 | WRITE(ITLIS,9981) JET | |
259 | 9981 FORMAT(//' ERROR IN WPAIR ... NO DECAY POSSIBLE FOR JET',I3) | |
260 | STOP 99 | |
261 | END |