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