1 | // SigmaNewGaugeBosons.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 Torbjorn Sjostrand. |
---|
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
5 | |
---|
6 | // Function definitions (not found in the header) for the |
---|
7 | // leptoquark simulation classes. |
---|
8 | |
---|
9 | #include "SigmaNewGaugeBosons.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | |
---|
14 | //========================================================================== |
---|
15 | |
---|
16 | // Sigma1ffbarZprimeWprime class. |
---|
17 | // Collects common methods for f fbar -> Z'/W' -> WW/WZ -> 4 fermions. |
---|
18 | // Copied from SigmaEW for gauge-boson-pair production. |
---|
19 | |
---|
20 | //-------------------------------------------------------------------------- |
---|
21 | |
---|
22 | // Calculate and store internal products. |
---|
23 | |
---|
24 | void Sigma1ffbarZprimeWprime::setupProd( Event& process, int i1, int i2, |
---|
25 | int i3, int i4, int i5, int i6) { |
---|
26 | |
---|
27 | // Store incoming and outgoing momenta, |
---|
28 | pRot[1] = process[i1].p(); |
---|
29 | pRot[2] = process[i2].p(); |
---|
30 | pRot[3] = process[i3].p(); |
---|
31 | pRot[4] = process[i4].p(); |
---|
32 | pRot[5] = process[i5].p(); |
---|
33 | pRot[6] = process[i6].p(); |
---|
34 | |
---|
35 | // Do random rotation to avoid accidental zeroes in HA expressions. |
---|
36 | bool smallPT = false; |
---|
37 | do { |
---|
38 | smallPT = false; |
---|
39 | double thetaNow = acos(2. * rndmPtr->flat() - 1.); |
---|
40 | double phiNow = 2. * M_PI * rndmPtr->flat(); |
---|
41 | for (int i = 1; i <= 6; ++i) { |
---|
42 | pRot[i].rot( thetaNow, phiNow); |
---|
43 | if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true; |
---|
44 | } |
---|
45 | } while (smallPT); |
---|
46 | |
---|
47 | // Calculate internal products. |
---|
48 | for (int i = 1; i < 6; ++i) { |
---|
49 | for (int j = i + 1; j <= 6; ++j) { |
---|
50 | hA[i][j] = |
---|
51 | sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz()) |
---|
52 | / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() ) |
---|
53 | - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz()) |
---|
54 | / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() ); |
---|
55 | hC[i][j] = conj( hA[i][j] ); |
---|
56 | if (i <= 2) { |
---|
57 | hA[i][j] *= complex( 0., 1.); |
---|
58 | hC[i][j] *= complex( 0., 1.); |
---|
59 | } |
---|
60 | hA[j][i] = - hA[i][j]; |
---|
61 | hC[j][i] = - hC[i][j]; |
---|
62 | } |
---|
63 | } |
---|
64 | |
---|
65 | } |
---|
66 | |
---|
67 | //-------------------------------------------------------------------------- |
---|
68 | |
---|
69 | // Evaluate the F function of Gunion and Kunszt. |
---|
70 | |
---|
71 | complex Sigma1ffbarZprimeWprime::fGK(int j1, int j2, int j3, int j4, |
---|
72 | int j5, int j6) { |
---|
73 | |
---|
74 | return 4. * hA[j1][j3] * hC[j2][j6] |
---|
75 | * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] ); |
---|
76 | |
---|
77 | } |
---|
78 | |
---|
79 | //-------------------------------------------------------------------------- |
---|
80 | |
---|
81 | // Evaluate the Xi function of Gunion and Kunszt. |
---|
82 | |
---|
83 | double Sigma1ffbarZprimeWprime::xiGK( double tHnow, double uHnow, |
---|
84 | double s3now, double s4now) { |
---|
85 | |
---|
86 | return - 4. * s3now * s4now + tHnow * (3. * tHnow + 4. * uHnow) |
---|
87 | + tHnow * tHnow * ( tHnow * uHnow / (s3now * s4now) |
---|
88 | - 2. * (1. / s3now + 1./s4now) * (tHnow + uHnow) |
---|
89 | + 2. * (s3now / s4now + s4now / s3now) ); |
---|
90 | |
---|
91 | } |
---|
92 | |
---|
93 | //-------------------------------------------------------------------------- |
---|
94 | |
---|
95 | // Evaluate the Xj function of Gunion and Kunszt. |
---|
96 | |
---|
97 | double Sigma1ffbarZprimeWprime::xjGK( double tHnow, double uHnow, |
---|
98 | double s3now, double s4now) { |
---|
99 | |
---|
100 | return 8. * pow2(s3now + s4now) - 8. * (s3now + s4now) * (tHnow + uHnow) |
---|
101 | - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow |
---|
102 | / (s3now * s4now) - 2. * (1. / s3now + 1. / s4now) * (tHnow + uHnow) |
---|
103 | + 2. * (s3now / s4now + s4now / s3now) ); |
---|
104 | |
---|
105 | } |
---|
106 | |
---|
107 | //========================================================================== |
---|
108 | |
---|
109 | // Sigma1ffbar2gmZZprime class. |
---|
110 | // Cross section for f fbar -> gamma*/Z0/Z'0 (f is quark or lepton). |
---|
111 | |
---|
112 | //-------------------------------------------------------------------------- |
---|
113 | |
---|
114 | // Initialize process. |
---|
115 | |
---|
116 | void Sigma1ffbar2gmZZprime::initProc() { |
---|
117 | |
---|
118 | // Allow to pick only parts of full gamma*/Z0/Z'0 expression. |
---|
119 | gmZmode = settingsPtr->mode("Zprime:gmZmode"); |
---|
120 | |
---|
121 | // Store Z'0 mass and width for propagator. |
---|
122 | mRes = particleDataPtr->m0(32); |
---|
123 | GammaRes = particleDataPtr->mWidth(32); |
---|
124 | m2Res = mRes*mRes; |
---|
125 | GamMRat = GammaRes / mRes; |
---|
126 | sin2tW = couplingsPtr->sin2thetaW(); |
---|
127 | cos2tW = 1. - sin2tW; |
---|
128 | thetaWRat = 1. / (16. * sin2tW * cos2tW); |
---|
129 | |
---|
130 | // Properties of Z0 resonance also needed. |
---|
131 | mZ = particleDataPtr->m0(23); |
---|
132 | GammaZ = particleDataPtr->mWidth(23); |
---|
133 | m2Z = mZ*mZ; |
---|
134 | GamMRatZ = GammaZ / mZ; |
---|
135 | |
---|
136 | // Ensure that arrays initially are empty. |
---|
137 | for (int i = 0; i < 20; ++i) afZp[i] = 0.; |
---|
138 | for (int i = 0; i < 20; ++i) vfZp[i] = 0.; |
---|
139 | |
---|
140 | // Store first-generation axial and vector couplings. |
---|
141 | afZp[1] = settingsPtr->parm("Zprime:ad"); |
---|
142 | afZp[2] = settingsPtr->parm("Zprime:au"); |
---|
143 | afZp[11] = settingsPtr->parm("Zprime:ae"); |
---|
144 | afZp[12] = settingsPtr->parm("Zprime:anue"); |
---|
145 | vfZp[1] = settingsPtr->parm("Zprime:vd"); |
---|
146 | vfZp[2] = settingsPtr->parm("Zprime:vu"); |
---|
147 | vfZp[11] = settingsPtr->parm("Zprime:ve"); |
---|
148 | vfZp[12] = settingsPtr->parm("Zprime:vnue"); |
---|
149 | |
---|
150 | // Second and third generation could be carbon copy of this... |
---|
151 | if (settingsPtr->flag("Zprime:universality")) { |
---|
152 | for (int i = 3; i <= 6; ++i) { |
---|
153 | afZp[i] = afZp[i-2]; |
---|
154 | vfZp[i] = vfZp[i-2]; |
---|
155 | afZp[i+10] = afZp[i+8]; |
---|
156 | vfZp[i+10] = vfZp[i+8]; |
---|
157 | } |
---|
158 | |
---|
159 | // ... or could have different couplings. |
---|
160 | } else { |
---|
161 | afZp[3] = settingsPtr->parm("Zprime:as"); |
---|
162 | afZp[4] = settingsPtr->parm("Zprime:ac"); |
---|
163 | afZp[5] = settingsPtr->parm("Zprime:ab"); |
---|
164 | afZp[6] = settingsPtr->parm("Zprime:at"); |
---|
165 | afZp[13] = settingsPtr->parm("Zprime:amu"); |
---|
166 | afZp[14] = settingsPtr->parm("Zprime:anumu"); |
---|
167 | afZp[15] = settingsPtr->parm("Zprime:atau"); |
---|
168 | afZp[16] = settingsPtr->parm("Zprime:anutau"); |
---|
169 | vfZp[3] = settingsPtr->parm("Zprime:vs"); |
---|
170 | vfZp[4] = settingsPtr->parm("Zprime:vc"); |
---|
171 | vfZp[5] = settingsPtr->parm("Zprime:vb"); |
---|
172 | vfZp[6] = settingsPtr->parm("Zprime:vt"); |
---|
173 | vfZp[13] = settingsPtr->parm("Zprime:vmu"); |
---|
174 | vfZp[14] = settingsPtr->parm("Zprime:vnumu"); |
---|
175 | vfZp[15] = settingsPtr->parm("Zprime:vtau"); |
---|
176 | vfZp[16] = settingsPtr->parm("Zprime:vnutau"); |
---|
177 | } |
---|
178 | |
---|
179 | // Coupling for Z' -> W+ W- and decay angular admixture. |
---|
180 | coupZpWW = settingsPtr->parm("Zprime:coup2WW"); |
---|
181 | anglesZpWW = settingsPtr->parm("Zprime:anglesWW"); |
---|
182 | |
---|
183 | // Set pointer to particle properties and decay table. |
---|
184 | particlePtr = particleDataPtr->particleDataEntryPtr(32); |
---|
185 | |
---|
186 | } |
---|
187 | |
---|
188 | //-------------------------------------------------------------------------- |
---|
189 | |
---|
190 | // Evaluate sigmaHat(sHat), part independent of incoming flavour. |
---|
191 | |
---|
192 | void Sigma1ffbar2gmZZprime::sigmaKin() { |
---|
193 | |
---|
194 | // Common coupling factors. |
---|
195 | double colQ = 3. * (1. + alpS / M_PI); |
---|
196 | |
---|
197 | // Reset quantities to sum. Declare variables in loop. |
---|
198 | gamSum = 0.; |
---|
199 | gamZSum = 0.; |
---|
200 | ZSum = 0.; |
---|
201 | gamZpSum = 0.; |
---|
202 | ZZpSum = 0.; |
---|
203 | ZpSum = 0.; |
---|
204 | int idAbs, onMode; |
---|
205 | double mf, mr, ps, kinFacA, kinFacV, ef, af, vf, apf, vpf, |
---|
206 | ef2, efvf, vaf2, efvpf, vafvapf, vapf2, colf; |
---|
207 | |
---|
208 | // Loop over all open Z'0 decay channels. |
---|
209 | for (int i = 0; i < particlePtr->sizeChannels(); ++i) { |
---|
210 | onMode = particlePtr->channel(i).onMode(); |
---|
211 | if (onMode != 1 && onMode != 2) continue; |
---|
212 | idAbs = abs( particlePtr->channel(i).product(0) ); |
---|
213 | |
---|
214 | // Contributions from three fermion generations. |
---|
215 | if ( (idAbs > 0 && idAbs < 7) || ( idAbs > 10 && idAbs < 17) ) { |
---|
216 | mf = particleDataPtr->m0(idAbs); |
---|
217 | |
---|
218 | // Check that above threshold. |
---|
219 | if (mH > 2. * mf + MASSMARGIN) { |
---|
220 | mr = pow2(mf / mH); |
---|
221 | ps = sqrtpos(1. - 4. * mr); |
---|
222 | |
---|
223 | // Couplings of gamma^*/Z^0/Z'^0 to final flavour |
---|
224 | ef = couplingsPtr->ef(idAbs); |
---|
225 | af = couplingsPtr->af(idAbs); |
---|
226 | vf = couplingsPtr->vf(idAbs); |
---|
227 | apf = afZp[idAbs]; |
---|
228 | vpf = vfZp[idAbs]; |
---|
229 | |
---|
230 | // Combine couplings with kinematical factors. |
---|
231 | kinFacA = pow3(ps); |
---|
232 | kinFacV = ps * (1. + 2. * mr); |
---|
233 | ef2 = ef * ef * kinFacV; |
---|
234 | efvf = ef * vf * kinFacV; |
---|
235 | vaf2 = vf * vf * kinFacV + af * af * kinFacA; |
---|
236 | efvpf = ef * vpf * kinFacV; |
---|
237 | vafvapf = vf * vpf * kinFacV + af * apf * kinFacA; |
---|
238 | vapf2 = vpf * vpf * kinFacV + apf * apf * kinFacA; |
---|
239 | |
---|
240 | // Colour factor. Additionally secondary width for top. |
---|
241 | colf = (idAbs < 9) ? colQ : 1.; |
---|
242 | if (idAbs == 6) colf *= particleDataPtr->resOpenFrac(6, -6); |
---|
243 | |
---|
244 | // Store sum of combinations. |
---|
245 | gamSum += colf * ef2; |
---|
246 | gamZSum += colf * efvf; |
---|
247 | ZSum += colf * vaf2; |
---|
248 | gamZpSum += colf * efvpf; |
---|
249 | ZZpSum += colf * vafvapf; |
---|
250 | ZpSum += colf * vapf2; |
---|
251 | } |
---|
252 | |
---|
253 | // Optional contribution from W+ W-. |
---|
254 | } else if (idAbs == 24) { |
---|
255 | mf = particleDataPtr->m0(idAbs); |
---|
256 | if (mH > 2. * mf + MASSMARGIN) { |
---|
257 | mr = pow2(mf / mH); |
---|
258 | ps = sqrtpos(1. - 4. * mr); |
---|
259 | ZpSum += pow2(coupZpWW * cos2tW) * pow3(ps) |
---|
260 | * (1. + 20. * mr + 12. * mr*mr) |
---|
261 | * particleDataPtr->resOpenFrac(24, -24); |
---|
262 | } |
---|
263 | } |
---|
264 | } |
---|
265 | |
---|
266 | // Calculate prefactors for gamma/Z0/Z'0 cross section terms. |
---|
267 | double propZ = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) ); |
---|
268 | double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); |
---|
269 | gamNorm = 4. * M_PI * pow2(alpEM) / (3. * sH); |
---|
270 | gamZNorm = gamNorm * 2. * thetaWRat * (sH - m2Z) * propZ; |
---|
271 | ZNorm = gamNorm * pow2(thetaWRat) * sH * propZ; |
---|
272 | gamZpNorm = gamNorm * 2. * thetaWRat * (sH - m2Res) * propZp; |
---|
273 | ZZpNorm = gamNorm * 2. * pow2(thetaWRat) * ((sH - m2Z) * (sH - m2Res) |
---|
274 | + sH * GamMRatZ * sH * GamMRat) * propZ * propZp; |
---|
275 | ZpNorm = gamNorm * pow2(thetaWRat) * sH * propZp; |
---|
276 | |
---|
277 | // Optionally only keep some of gamma*, Z0 and Z' terms. |
---|
278 | if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.; |
---|
279 | ZZpNorm = 0.; ZpNorm = 0.;} |
---|
280 | if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.; |
---|
281 | ZZpNorm = 0.; ZpNorm = 0.;} |
---|
282 | if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.; |
---|
283 | gamZpNorm = 0.; ZZpNorm = 0.;} |
---|
284 | if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;} |
---|
285 | if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;} |
---|
286 | if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;} |
---|
287 | |
---|
288 | } |
---|
289 | |
---|
290 | //-------------------------------------------------------------------------- |
---|
291 | |
---|
292 | // Evaluate sigmaHat(sHat), including incoming flavour dependence. |
---|
293 | |
---|
294 | double Sigma1ffbar2gmZZprime::sigmaHat() { |
---|
295 | |
---|
296 | // Couplings to an incoming flavour. |
---|
297 | int idAbs = abs(id1); |
---|
298 | double ei = couplingsPtr->ef(idAbs); |
---|
299 | double ai = couplingsPtr->af(idAbs); |
---|
300 | double vi = couplingsPtr->vf(idAbs); |
---|
301 | double api = afZp[idAbs]; |
---|
302 | double vpi = vfZp[idAbs]; |
---|
303 | double ei2 = ei * ei; |
---|
304 | double eivi = ei * vi; |
---|
305 | double vai2 = vi * vi + ai * ai; |
---|
306 | double eivpi = ei * vpi; |
---|
307 | double vaivapi = vi * vpi + ai * api;; |
---|
308 | double vapi2 = vpi * vpi + api * api; |
---|
309 | |
---|
310 | // Combine gamma, interference and Z0 parts. |
---|
311 | double sigma = ei2 * gamNorm * gamSum + eivi * gamZNorm * gamZSum |
---|
312 | + vai2 * ZNorm * ZSum + eivpi * gamZpNorm * gamZpSum |
---|
313 | + vaivapi * ZZpNorm * ZZpSum + vapi2 * ZpNorm * ZpSum; |
---|
314 | |
---|
315 | // Colour factor. Answer. |
---|
316 | if (idAbs < 9) sigma /= 3.; |
---|
317 | return sigma; |
---|
318 | |
---|
319 | } |
---|
320 | |
---|
321 | //-------------------------------------------------------------------------- |
---|
322 | |
---|
323 | // Select identity, colour and anticolour. |
---|
324 | |
---|
325 | void Sigma1ffbar2gmZZprime::setIdColAcol() { |
---|
326 | |
---|
327 | // Flavours trivial. |
---|
328 | setId( id1, id2, 32); |
---|
329 | |
---|
330 | // Colour flow topologies. Swap when antiquarks. |
---|
331 | if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0); |
---|
332 | else setColAcol( 0, 0, 0, 0, 0, 0); |
---|
333 | if (id1 < 0) swapColAcol(); |
---|
334 | |
---|
335 | } |
---|
336 | |
---|
337 | //-------------------------------------------------------------------------- |
---|
338 | |
---|
339 | // Evaluate weight for gamma*/Z0/Z'0 decay angle. |
---|
340 | |
---|
341 | double Sigma1ffbar2gmZZprime::weightDecay( Event& process, int iResBeg, |
---|
342 | int iResEnd) { |
---|
343 | |
---|
344 | // Default values, in- and out-flavours in process. |
---|
345 | double wt = 1.; |
---|
346 | double wtMax = 1.; |
---|
347 | int idInAbs = process[3].idAbs(); |
---|
348 | int idOutAbs = process[6].idAbs(); |
---|
349 | |
---|
350 | // Angular weight for outgoing fermion pair. |
---|
351 | if (iResBeg == 5 && iResEnd == 5 && |
---|
352 | (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) { |
---|
353 | |
---|
354 | // Couplings for in- and out-flavours. |
---|
355 | double ei = couplingsPtr->ef(idInAbs); |
---|
356 | double vi = couplingsPtr->vf(idInAbs); |
---|
357 | double ai = couplingsPtr->af(idInAbs); |
---|
358 | double vpi = vfZp[idInAbs]; |
---|
359 | double api = afZp[idInAbs]; |
---|
360 | double ef = couplingsPtr->ef(idOutAbs); |
---|
361 | double vf = couplingsPtr->vf(idOutAbs); |
---|
362 | double af = couplingsPtr->af(idOutAbs); |
---|
363 | double vpf = vfZp[idOutAbs]; |
---|
364 | double apf = afZp[idOutAbs]; |
---|
365 | |
---|
366 | // Phase space factors. (One power of beta left out in formulae.) |
---|
367 | double mr1 = pow2(process[6].m()) / sH; |
---|
368 | double mr2 = pow2(process[7].m()) / sH; |
---|
369 | double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); |
---|
370 | double mrAvg = 0.5 * (mr1 + mr2) - 0.25 * pow2(mr1 - mr2); |
---|
371 | |
---|
372 | // Coefficients of angular expression. |
---|
373 | double coefTran = ei*ei * gamNorm * ef*ef + ei * vi * gamZNorm * ef * vf |
---|
374 | + (vi*vi + ai*ai) * ZNorm * (vf*vf + ps*ps * af*af) |
---|
375 | + ei * vpi * gamZpNorm * ef * vpf |
---|
376 | + (vi * vpi + ai * api) * ZZpNorm * (vf * vpf + ps*ps * af * apf) |
---|
377 | + (vpi*vpi + api*api) * ZpNorm * (vpf*vpf + ps*ps * apf*apf); |
---|
378 | double coefLong = 4. * mrAvg * ( ei*ei * gamNorm * ef*ef |
---|
379 | + ei * vi * gamZNorm * ef * vf + (vi*vi + ai*ai) * ZNorm * vf*vf |
---|
380 | + ei * vpi * gamZpNorm * ef * vpf |
---|
381 | + (vi * vpi + ai * api) * ZZpNorm * vf * vpf |
---|
382 | + (vpi*vpi + api*api) * ZpNorm * vpf*vpf ); |
---|
383 | double coefAsym = ps * ( ei * ai * gamZNorm * ef * af |
---|
384 | + 4. * vi * ai * ZNorm * vf * af + ei * api * gamZpNorm * ef * apf |
---|
385 | + (vi * api + vpi * ai) * ZZpNorm * (vf * apf + vpf * af) |
---|
386 | + 4. * vpi * api * ZpNorm * vpf * apf ); |
---|
387 | |
---|
388 | // Flip asymmetry for in-fermion + out-antifermion. |
---|
389 | if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym; |
---|
390 | |
---|
391 | // Reconstruct decay angle and weight for it. |
---|
392 | double cosThe = (process[3].p() - process[4].p()) |
---|
393 | * (process[7].p() - process[6].p()) / (sH * ps); |
---|
394 | wt = coefTran * (1. + pow2(cosThe)) |
---|
395 | + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe; |
---|
396 | wtMax = 2. * (coefTran + abs(coefAsym)); |
---|
397 | } |
---|
398 | |
---|
399 | // Angular weight for Z' -> W+ W-. |
---|
400 | else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) { |
---|
401 | double mr1 = pow2(process[6].m()) / sH; |
---|
402 | double mr2 = pow2(process[7].m()) / sH; |
---|
403 | double ps = sqrtpos(pow2(1. - mr1 -mr2) - 4. * mr1 * mr2); |
---|
404 | double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2 |
---|
405 | + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2); |
---|
406 | double cFlat = -cCos2 + 0.5 * (mr1 + mr2) |
---|
407 | * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2)); |
---|
408 | |
---|
409 | // Reconstruct decay angle and weight for it. |
---|
410 | double cosThe = (process[3].p() - process[4].p()) |
---|
411 | * (process[7].p() - process[6].p()) / (sH * ps); |
---|
412 | wt = cFlat + cCos2 * cosThe*cosThe; |
---|
413 | wtMax = cFlat + max(0., cCos2); |
---|
414 | } |
---|
415 | |
---|
416 | // Angular weight for f + fbar -> Z' -> W+ + W- -> 4 fermions. |
---|
417 | else if (iResBeg == 6 && iResEnd == 7 && idOutAbs == 24) { |
---|
418 | |
---|
419 | // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6). |
---|
420 | // with f' fbar' from W- and f" fbar" from W+. |
---|
421 | int i1 = (process[3].id() < 0) ? 3 : 4; |
---|
422 | int i2 = 7 - i1; |
---|
423 | int i3 = (process[8].id() > 0) ? 8 : 9; |
---|
424 | int i4 = 17 - i3; |
---|
425 | int i5 = (process[10].id() > 0) ? 10 : 11; |
---|
426 | int i6 = 21 - i5; |
---|
427 | if (process[6].id() > 0) {swap(i3, i5); swap(i4, i6);} |
---|
428 | |
---|
429 | // Decay distribution like in f fbar -> Z^* -> W+ W-. |
---|
430 | if (rndmPtr->flat() > anglesZpWW) { |
---|
431 | |
---|
432 | // Set up four-products and internal products. |
---|
433 | setupProd( process, i1, i2, i3, i4, i5, i6); |
---|
434 | |
---|
435 | // tHat and uHat of fbar f -> W- W+, and their squared masses. |
---|
436 | int iNeg = (process[6].id() < 0) ? 6 : 7; |
---|
437 | int iPos = 13 - iNeg; |
---|
438 | double tHres = (process[i1].p() - process[iNeg].p()).m2Calc(); |
---|
439 | double uHres = (process[i1].p() - process[iPos].p()).m2Calc(); |
---|
440 | double s3now = process[iNeg].m2(); |
---|
441 | double s4now = process[iPos].m2(); |
---|
442 | |
---|
443 | // Kinematics combinations (norm(x) = |x|^2). |
---|
444 | double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) ); |
---|
445 | double fGK253 = norm(fGK( 2, 1, 5, 6, 3, 4) - fGK( 2, 1, 3, 4, 5, 6) ); |
---|
446 | double xiT = xiGK( tHres, uHres, s3now, s4now); |
---|
447 | double xiU = xiGK( uHres, tHres, s3now, s4now); |
---|
448 | double xjTU = xjGK( tHres, uHres, s3now, s4now); |
---|
449 | |
---|
450 | // Couplings of incoming (anti)fermion. Combine with kinematics. |
---|
451 | int idAbs = process[i1].idAbs(); |
---|
452 | double li = 0.5 * (vfZp[idAbs] + afZp[idAbs]); |
---|
453 | double ri = 0.5 * (vfZp[idAbs] - afZp[idAbs]); |
---|
454 | wt = li*li * fGK135 + ri*ri * fGK253; |
---|
455 | wtMax = 4. * s3now * s4now * (li*li + ri*ri) |
---|
456 | * (xiT + xiU - xjTU); |
---|
457 | |
---|
458 | // Decay distribution like in f fbar -> h^0 -> W+ W-. |
---|
459 | } else { |
---|
460 | double p35 = 2. * process[i3].p() * process[i5].p(); |
---|
461 | double p46 = 2. * process[i4].p() * process[i6].p(); |
---|
462 | wt = 16. * p35 * p46; |
---|
463 | wtMax = sH2; |
---|
464 | } |
---|
465 | } |
---|
466 | |
---|
467 | // Angular weight in top decay by standard routine. |
---|
468 | else if (process[process[iResBeg].mother1()].idAbs() == 6) |
---|
469 | return weightTopDecay( process, iResBeg, iResEnd); |
---|
470 | |
---|
471 | |
---|
472 | // Done. |
---|
473 | return (wt / wtMax); |
---|
474 | |
---|
475 | } |
---|
476 | |
---|
477 | //========================================================================== |
---|
478 | |
---|
479 | // Sigma1ffbar2Wprime class. |
---|
480 | // Cross section for f fbar' -> W'+- (f is quark or lepton). |
---|
481 | |
---|
482 | //-------------------------------------------------------------------------- |
---|
483 | |
---|
484 | // Initialize process. |
---|
485 | |
---|
486 | void Sigma1ffbar2Wprime::initProc() { |
---|
487 | |
---|
488 | // Store W+- mass and width for propagator. |
---|
489 | mRes = particleDataPtr->m0(34); |
---|
490 | GammaRes = particleDataPtr->mWidth(34); |
---|
491 | m2Res = mRes*mRes; |
---|
492 | GamMRat = GammaRes / mRes; |
---|
493 | thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW()); |
---|
494 | |
---|
495 | // Axial and vector couplings of fermions. |
---|
496 | aqWp = settingsPtr->parm("Wprime:aq"); |
---|
497 | vqWp = settingsPtr->parm("Wprime:vq"); |
---|
498 | alWp = settingsPtr->parm("Wprime:al"); |
---|
499 | vlWp = settingsPtr->parm("Wprime:vl"); |
---|
500 | |
---|
501 | // Coupling for W' -> W Z and decay angular admixture. |
---|
502 | coupWpWZ = settingsPtr->parm("Wprime:coup2WZ"); |
---|
503 | anglesWpWZ = settingsPtr->parm("Wprime:anglesWZ"); |
---|
504 | |
---|
505 | // Set pointer to particle properties and decay table. |
---|
506 | particlePtr = particleDataPtr->particleDataEntryPtr(34); |
---|
507 | |
---|
508 | } |
---|
509 | |
---|
510 | //-------------------------------------------------------------------------- |
---|
511 | |
---|
512 | // Evaluate sigmaHat(sHat), part independent of incoming flavour. |
---|
513 | |
---|
514 | void Sigma1ffbar2Wprime::sigmaKin() { |
---|
515 | |
---|
516 | // Set up Breit-Wigner. Cross section for W+ and W- separately. |
---|
517 | double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); |
---|
518 | double preFac = alpEM * thetaWRat * mH; |
---|
519 | sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(34, mH); |
---|
520 | sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-34, mH); |
---|
521 | |
---|
522 | } |
---|
523 | |
---|
524 | //-------------------------------------------------------------------------- |
---|
525 | |
---|
526 | // Evaluate sigmaHat(sHat), including incoming flavour dependence. |
---|
527 | |
---|
528 | double Sigma1ffbar2Wprime::sigmaHat() { |
---|
529 | |
---|
530 | // Secondary width for W+ or W-. CKM and colour factors. |
---|
531 | int idUp = (abs(id1)%2 == 0) ? id1 : id2; |
---|
532 | double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg; |
---|
533 | if (abs(id1) < 7) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.; |
---|
534 | |
---|
535 | // Couplings. |
---|
536 | if (abs(id1) < 7) sigma *= 0.5 * (aqWp * aqWp + vqWp * vqWp); |
---|
537 | else sigma *= 0.5 * (alWp * alWp + vlWp * vlWp); |
---|
538 | |
---|
539 | // Answer. |
---|
540 | return sigma; |
---|
541 | |
---|
542 | } |
---|
543 | |
---|
544 | //-------------------------------------------------------------------------- |
---|
545 | |
---|
546 | // Select identity, colour and anticolour. |
---|
547 | |
---|
548 | void Sigma1ffbar2Wprime::setIdColAcol() { |
---|
549 | |
---|
550 | // Sign of outgoing W. |
---|
551 | int sign = 1 - 2 * (abs(id1)%2); |
---|
552 | if (id1 < 0) sign = -sign; |
---|
553 | setId( id1, id2, 34 * sign); |
---|
554 | |
---|
555 | // Colour flow topologies. Swap when antiquarks. |
---|
556 | if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0); |
---|
557 | else setColAcol( 0, 0, 0, 0, 0, 0); |
---|
558 | if (id1 < 0) swapColAcol(); |
---|
559 | |
---|
560 | } |
---|
561 | |
---|
562 | //-------------------------------------------------------------------------- |
---|
563 | |
---|
564 | // Evaluate weight for W decay angle. |
---|
565 | |
---|
566 | double Sigma1ffbar2Wprime::weightDecay( Event& process, int iResBeg, |
---|
567 | int iResEnd) { |
---|
568 | |
---|
569 | // Default values, in- and out-flavours in process. |
---|
570 | double wt = 1.; |
---|
571 | double wtMax = 1.; |
---|
572 | int idInAbs = process[3].idAbs(); |
---|
573 | int idOutAbs = process[6].idAbs(); |
---|
574 | |
---|
575 | // Angular weight for outgoing fermion pair. |
---|
576 | if (iResBeg == 5 && iResEnd == 5 && |
---|
577 | (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) { |
---|
578 | |
---|
579 | // Couplings for in- and out-flavours. |
---|
580 | double ai = (idInAbs < 9) ? aqWp : alWp; |
---|
581 | double vi = (idInAbs < 9) ? vqWp : vlWp; |
---|
582 | double af = (idOutAbs < 9) ? aqWp : alWp; |
---|
583 | double vf = (idOutAbs < 9) ? vqWp : vlWp; |
---|
584 | |
---|
585 | // Asymmetry expression. |
---|
586 | double coefAsym = 8. * vi * ai * vf * af |
---|
587 | / ((vi*vi + ai*ai) * (vf*vf + af*af)); |
---|
588 | |
---|
589 | // Flip asymmetry for in-fermion + out-antifermion. |
---|
590 | if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym; |
---|
591 | |
---|
592 | // Phase space factors. |
---|
593 | double mr1 = pow2(process[6].m()) / sH; |
---|
594 | double mr2 = pow2(process[7].m()) / sH; |
---|
595 | double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); |
---|
596 | |
---|
597 | // Reconstruct decay angle and weight for it. |
---|
598 | double cosThe = (process[3].p() - process[4].p()) |
---|
599 | * (process[7].p() - process[6].p()) / (sH * ps); |
---|
600 | wt = 1. + coefAsym * cosThe + cosThe * cosThe; |
---|
601 | wtMax = 2. + abs(coefAsym); |
---|
602 | } |
---|
603 | |
---|
604 | // Angular weight for W' -> W Z. |
---|
605 | else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) { |
---|
606 | double mr1 = pow2(process[6].m()) / sH; |
---|
607 | double mr2 = pow2(process[7].m()) / sH; |
---|
608 | double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); |
---|
609 | double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2 |
---|
610 | + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2); |
---|
611 | double cFlat = -cCos2 + 0.5 * (mr1 + mr2) |
---|
612 | * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2)); |
---|
613 | |
---|
614 | // Reconstruct decay angle and weight for it. |
---|
615 | double cosThe = (process[3].p() - process[4].p()) |
---|
616 | * (process[7].p() - process[6].p()) / (sH * ps); |
---|
617 | wt = cFlat + cCos2 * cosThe*cosThe; |
---|
618 | wtMax = cFlat + max(0., cCos2); |
---|
619 | } |
---|
620 | |
---|
621 | // Angular weight for f + fbar -> W' -> W + Z -> 4 fermions. |
---|
622 | else if (iResBeg == 6 && iResEnd == 7 |
---|
623 | && (idOutAbs == 24 || idOutAbs == 23)) { |
---|
624 | |
---|
625 | // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6). |
---|
626 | // with f' fbar' from W and f" fbar" from Z. |
---|
627 | int i1 = (process[3].id() < 0) ? 3 : 4; |
---|
628 | int i2 = 7 - i1; |
---|
629 | int i3 = (process[8].id() > 0) ? 8 : 9; |
---|
630 | int i4 = 17 - i3; |
---|
631 | int i5 = (process[10].id() > 0) ? 10 : 11; |
---|
632 | int i6 = 21 - i5; |
---|
633 | if (process[6].id() == 23) {swap(i3, i5); swap(i4, i6);} |
---|
634 | |
---|
635 | // Decay distribution like in f fbar -> Z^* -> W+ W-. |
---|
636 | if (rndmPtr->flat() > anglesWpWZ) { |
---|
637 | |
---|
638 | // Set up four-products and internal products. |
---|
639 | setupProd( process, i1, i2, i3, i4, i5, i6); |
---|
640 | |
---|
641 | // tHat and uHat of fbar f -> W Z, and their squared masses. |
---|
642 | int iW = (process[6].id() == 23) ? 7 : 6; |
---|
643 | int iZ = 13 - iW; |
---|
644 | double tHres = (process[i1].p() - process[iW].p()).m2Calc(); |
---|
645 | double uHres = (process[i1].p() - process[iZ].p()).m2Calc(); |
---|
646 | double s3now = process[iW].m2(); |
---|
647 | double s4now = process[iZ].m2(); |
---|
648 | |
---|
649 | // Kinematics combinations (norm(x) = |x|^2). |
---|
650 | double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) ); |
---|
651 | double fGK136 = norm(fGK( 1, 2, 3, 4, 6, 5) - fGK( 1, 2, 6, 5, 3, 4) ); |
---|
652 | double xiT = xiGK( tHres, uHres, s3now, s4now); |
---|
653 | double xiU = xiGK( uHres, tHres, s3now, s4now); |
---|
654 | double xjTU = xjGK( tHres, uHres, s3now, s4now); |
---|
655 | |
---|
656 | // Couplings of outgoing fermion from Z. Combine with kinematics. |
---|
657 | int idAbs = process[i5].idAbs(); |
---|
658 | double lfZ = couplingsPtr->lf(idAbs); |
---|
659 | double rfZ = couplingsPtr->rf(idAbs); |
---|
660 | wt = lfZ*lfZ * fGK135 + rfZ*rfZ * fGK136; |
---|
661 | wtMax = 4. * s3now * s4now * (lfZ*lfZ + rfZ*rfZ) |
---|
662 | * (xiT + xiU - xjTU); |
---|
663 | |
---|
664 | // Decay distribution like in f fbar -> H^+- -> W+- Z0. |
---|
665 | } else { |
---|
666 | double p35 = 2. * process[i3].p() * process[i5].p(); |
---|
667 | double p46 = 2. * process[i4].p() * process[i6].p(); |
---|
668 | wt = 16. * p35 * p46; |
---|
669 | wtMax = sH2; |
---|
670 | } |
---|
671 | } |
---|
672 | |
---|
673 | // Angular weight in top decay by standard routine. |
---|
674 | else if (process[process[iResBeg].mother1()].idAbs() == 6) |
---|
675 | return weightTopDecay( process, iResBeg, iResEnd); |
---|
676 | |
---|
677 | // Done. |
---|
678 | return (wt / wtMax); |
---|
679 | |
---|
680 | } |
---|
681 | |
---|
682 | |
---|
683 | //========================================================================== |
---|
684 | |
---|
685 | // Sigma1ffbar2Rhorizontal class. |
---|
686 | // Cross section for f fbar' -> R^0 (f is a quark or lepton). |
---|
687 | |
---|
688 | //-------------------------------------------------------------------------- |
---|
689 | |
---|
690 | // Initialize process. |
---|
691 | |
---|
692 | void Sigma1ffbar2Rhorizontal::initProc() { |
---|
693 | |
---|
694 | // Store R^0 mass and width for propagator. |
---|
695 | mRes = particleDataPtr->m0(41); |
---|
696 | GammaRes = particleDataPtr->mWidth(41); |
---|
697 | m2Res = mRes*mRes; |
---|
698 | GamMRat = GammaRes / mRes; |
---|
699 | thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW()); |
---|
700 | |
---|
701 | // Set pointer to particle properties and decay table. |
---|
702 | particlePtr = particleDataPtr->particleDataEntryPtr(41); |
---|
703 | |
---|
704 | } |
---|
705 | |
---|
706 | //-------------------------------------------------------------------------- |
---|
707 | |
---|
708 | // Evaluate sigmaHat(sHat), part independent of incoming flavour. |
---|
709 | |
---|
710 | void Sigma1ffbar2Rhorizontal::sigmaKin() { |
---|
711 | |
---|
712 | // Set up Breit-Wigner. Cross section for W+ and W- separately. |
---|
713 | double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); |
---|
714 | double preFac = alpEM * thetaWRat * mH; |
---|
715 | sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(41, mH); |
---|
716 | sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-41, mH); |
---|
717 | |
---|
718 | } |
---|
719 | |
---|
720 | //-------------------------------------------------------------------------- |
---|
721 | |
---|
722 | // Evaluate sigmaHat(sHat), including incoming flavour dependence. |
---|
723 | |
---|
724 | double Sigma1ffbar2Rhorizontal::sigmaHat() { |
---|
725 | |
---|
726 | // Check for allowed flavour combinations, one generation apart. |
---|
727 | if (id1 * id2 > 0 || abs(id1 + id2) != 2) return 0.; |
---|
728 | |
---|
729 | // Find whether R0 or R0bar. Colour factors. |
---|
730 | double sigma = (id1 + id2 > 0) ? sigma0Pos : sigma0Neg; |
---|
731 | if (abs(id1) < 7) sigma /= 3.; |
---|
732 | |
---|
733 | // Answer. |
---|
734 | return sigma; |
---|
735 | |
---|
736 | } |
---|
737 | |
---|
738 | //-------------------------------------------------------------------------- |
---|
739 | |
---|
740 | // Select identity, colour and anticolour. |
---|
741 | |
---|
742 | void Sigma1ffbar2Rhorizontal::setIdColAcol() { |
---|
743 | |
---|
744 | // Outgoing R0 or R0bar. |
---|
745 | id3 = (id1 +id2 > 0) ? 41 : -41; |
---|
746 | setId( id1, id2, id3); |
---|
747 | |
---|
748 | // Colour flow topologies. Swap when antiquarks. |
---|
749 | if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0); |
---|
750 | else setColAcol( 0, 0, 0, 0, 0, 0); |
---|
751 | if (id1 < 0) swapColAcol(); |
---|
752 | |
---|
753 | } |
---|
754 | |
---|
755 | //========================================================================== |
---|
756 | |
---|
757 | } // end namespace Pythia8 |
---|