Line data Source code
1 : // SigmaCompositeness.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 : // compositeness simulation classes.
8 :
9 : #include "Pythia8/SigmaCompositeness.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // Sigma1qg2qStar class.
16 : // Cross section for q g -> q^* (excited quark state).
17 :
18 : //--------------------------------------------------------------------------
19 :
20 : // Initialize process.
21 :
22 : void Sigma1qg2qStar::initProc() {
23 :
24 : // Set up process properties from the chosen quark flavour.
25 0 : idRes = 4000000 + idq;
26 0 : codeSave = 4000 + idq;
27 0 : if (idq == 1) nameSave = "d g -> d^*";
28 0 : else if (idq == 2) nameSave = "u g -> u^*";
29 0 : else if (idq == 3) nameSave = "s g -> s^*";
30 0 : else if (idq == 4) nameSave = "c g -> c^*";
31 0 : else nameSave = "b g -> b^*";
32 :
33 : // Store q* mass and width for propagator.
34 0 : mRes = particleDataPtr->m0(idRes);
35 0 : GammaRes = particleDataPtr->mWidth(idRes);
36 0 : m2Res = mRes*mRes;
37 0 : GamMRat = GammaRes / mRes;
38 :
39 : // Locally stored properties and couplings.
40 0 : Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
41 0 : coupFcol = settingsPtr->parm("ExcitedFermion:coupFcol");
42 :
43 : // Set pointer to particle properties and decay table.
44 0 : qStarPtr = particleDataPtr->particleDataEntryPtr(idRes);
45 :
46 0 : }
47 :
48 : //--------------------------------------------------------------------------
49 :
50 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
51 :
52 : void Sigma1qg2qStar::sigmaKin() {
53 :
54 : // Incoming width for correct quark.
55 0 : widthIn = pow3(mH) * alpS * pow2(coupFcol) / (3. * pow2(Lambda));
56 :
57 : // Set up Breit-Wigner.
58 0 : sigBW = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
59 :
60 0 : }
61 :
62 : //--------------------------------------------------------------------------
63 :
64 : // Evaluate sigmaHat(sHat) for specific incoming flavours.
65 :
66 : double Sigma1qg2qStar::sigmaHat() {
67 :
68 : // Identify whether correct incoming flavours.
69 0 : int idqNow = (id2 == 21) ? id1 : id2;
70 0 : if (abs(idqNow) != idq) return 0.;
71 :
72 : // Outgoing width and total sigma. Done.
73 0 : return widthIn * sigBW * qStarPtr->resWidthOpen(idqNow, mH);
74 :
75 0 : }
76 :
77 : //--------------------------------------------------------------------------
78 :
79 : // Select identity, colour and anticolour.
80 :
81 : void Sigma1qg2qStar::setIdColAcol() {
82 :
83 : // Flavours.
84 0 : int idqNow = (id2 == 21) ? id1 : id2;
85 0 : int idqStar = (idqNow > 0) ? idRes : -idRes;
86 0 : setId( id1, id2, idqStar);
87 :
88 : // Colour flow topology.
89 0 : if (id1 == idqNow) setColAcol( 1, 0, 2, 1, 2, 0);
90 0 : else setColAcol( 2, 1, 1, 0, 2, 0);
91 0 : if (idqNow < 0) swapColAcol();
92 :
93 0 : }
94 :
95 : //--------------------------------------------------------------------------
96 :
97 : // Evaluate weight for q* decay angle.
98 :
99 : double Sigma1qg2qStar::weightDecay( Event& process, int iResBeg,
100 : int iResEnd) {
101 :
102 : // q* should sit in entry 5. Sequential Z/W decay assumed isotropic.
103 0 : if (iResBeg != 5 || iResEnd != 5) return 1.;
104 :
105 : // Sign of asymmetry.
106 0 : int sideIn = (process[3].idAbs() < 20) ? 1 : 2;
107 0 : int sideOut = (process[6].idAbs() < 20) ? 1 : 2;
108 0 : double eps = (sideIn == sideOut) ? 1. : -1.;
109 :
110 : // Phase space factors.
111 0 : double mr1 = pow2(process[6].m()) / sH;
112 0 : double mr2 = pow2(process[7].m()) / sH;
113 0 : double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
114 :
115 : // Reconstruct decay angle. Default isotropic decay.
116 0 : double cosThe = (process[3].p() - process[4].p())
117 0 : * (process[7].p() - process[6].p()) / (sH * betaf);
118 : double wt = 1.;
119 : double wtMax = 1.;
120 :
121 : // Decay q* -> q (g/gamma) or q (Z^0/W^+-).
122 0 : int idBoson = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs();
123 0 : if (idBoson == 21 || idBoson == 22) {
124 0 : wt = 1. + eps * cosThe;
125 : wtMax = 2.;
126 0 : } else if (idBoson == 23 || idBoson == 24) {
127 0 : double mrB = (sideOut == 1) ? mr2 : mr1;
128 0 : double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB);
129 0 : wt = 1. + eps * cosThe * ratB;
130 0 : wtMax = 1. + ratB;
131 0 : }
132 :
133 : // Done.
134 0 : return (wt / wtMax);
135 :
136 0 : }
137 :
138 : //==========================================================================
139 :
140 : // Sigma1lgm2lStar class.
141 : // Cross section for l gamma -> l^* (excited lepton state).
142 :
143 : //--------------------------------------------------------------------------
144 :
145 : // Initialize process.
146 :
147 : void Sigma1lgm2lStar::initProc() {
148 :
149 : // Set up process properties from the chosen lepton flavour.
150 0 : idRes = 4000000 + idl;
151 0 : codeSave = 4000 + idl;
152 0 : if (idl == 11) nameSave = "e gamma -> e^*";
153 0 : else if (idl == 13) nameSave = "mu gamma -> mu^*";
154 0 : else nameSave = "tau gamma -> tau^*";
155 :
156 : // Store l* mass and width for propagator.
157 0 : mRes = particleDataPtr->m0(idRes);
158 0 : GammaRes = particleDataPtr->mWidth(idRes);
159 0 : m2Res = mRes*mRes;
160 0 : GamMRat = GammaRes / mRes;
161 :
162 : // Locally stored properties and couplings.
163 0 : Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
164 0 : double coupF = settingsPtr->parm("ExcitedFermion:coupF");
165 0 : double coupFp = settingsPtr->parm("ExcitedFermion:coupFprime");
166 0 : coupChg = -0.5 * coupF - 0.5 * coupFp;
167 :
168 : // Set pointer to particle properties and decay table.
169 0 : qStarPtr = particleDataPtr->particleDataEntryPtr(idRes);
170 :
171 0 : }
172 :
173 : //--------------------------------------------------------------------------
174 :
175 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
176 :
177 : void Sigma1lgm2lStar::sigmaKin() {
178 :
179 : // Incoming width for correct lepton.
180 0 : widthIn = pow3(mH) * alpEM * pow2(coupChg) / pow2(Lambda);
181 :
182 : // Set up Breit-Wigner.
183 0 : sigBW = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
184 :
185 0 : }
186 :
187 : //--------------------------------------------------------------------------
188 :
189 : // Evaluate sigmaHat(sHat) for specific incoming flavours.
190 :
191 : double Sigma1lgm2lStar::sigmaHat() {
192 :
193 : // Identify whether correct incoming flavours.
194 0 : int idlNow = (id2 == 22) ? id1 : id2;
195 0 : if (abs(idlNow) != idl) return 0.;
196 :
197 : // Outgoing width and total sigma. Done.
198 0 : return widthIn * sigBW * qStarPtr->resWidthOpen(idlNow, mH);
199 :
200 0 : }
201 :
202 : //--------------------------------------------------------------------------
203 :
204 : // Select identity, colour and anticolour.
205 :
206 : void Sigma1lgm2lStar::setIdColAcol() {
207 :
208 : // Flavours.
209 0 : int idlNow = (id2 == 22) ? id1 : id2;
210 0 : int idlStar = (idlNow > 0) ? idRes : -idRes;
211 0 : setId( id1, id2, idlStar);
212 :
213 : // No colour flow.
214 0 : setColAcol( 0, 0, 0, 0, 0, 0);
215 :
216 0 : }
217 :
218 : //--------------------------------------------------------------------------
219 :
220 : // Evaluate weight for l* decay angle.
221 :
222 : double Sigma1lgm2lStar::weightDecay( Event& process, int iResBeg,
223 : int iResEnd) {
224 :
225 : // l* should sit in entry 5. Sequential Z/W decay assumed isotropic.
226 0 : if (iResBeg != 5 || iResEnd != 5) return 1.;
227 :
228 : // Sign of asymmetry.
229 0 : int sideIn = (process[3].idAbs() < 20) ? 1 : 2;
230 0 : int sideOut = (process[6].idAbs() < 20) ? 1 : 2;
231 0 : double eps = (sideIn == sideOut) ? 1. : -1.;
232 :
233 : // Phase space factors.
234 0 : double mr1 = pow2(process[6].m()) / sH;
235 0 : double mr2 = pow2(process[7].m()) / sH;
236 0 : double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
237 :
238 : // Reconstruct decay angle. Default isotropic decay.
239 0 : double cosThe = (process[3].p() - process[4].p())
240 0 : * (process[7].p() - process[6].p()) / (sH * betaf);
241 : double wt = 1.;
242 : double wtMax = 1.;
243 :
244 : // Decay l* -> l gamma or l (Z^0/W^+-).
245 0 : int idBoson = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs();
246 0 : if (idBoson == 22) {
247 0 : wt = 1. + eps * cosThe;
248 : wtMax = 2.;
249 0 : } else if (idBoson == 23 || idBoson == 24) {
250 0 : double mrB = (sideOut == 1) ? mr2 : mr1;
251 0 : double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB);
252 0 : wt = 1. + eps * cosThe * ratB;
253 0 : wtMax = 1. + ratB;
254 0 : }
255 :
256 : // Done.
257 0 : return (wt / wtMax);
258 :
259 0 : }
260 :
261 : //==========================================================================
262 :
263 : // Sigma2qq2qStarq class.
264 : // Cross section for q q' -> q^* q' (excited quark state).
265 :
266 : //--------------------------------------------------------------------------
267 :
268 : // Initialize process.
269 :
270 : void Sigma2qq2qStarq::initProc() {
271 :
272 : // Set up process properties from the chosen quark flavour.
273 0 : idRes = 4000000 + idq;
274 0 : codeSave = 4020 + idq;
275 0 : if (idq == 1) nameSave = "q q -> d^* q";
276 0 : else if (idq == 2) nameSave = "q q -> u^* q";
277 0 : else if (idq == 3) nameSave = "q q -> s^* q";
278 0 : else if (idq == 4) nameSave = "q q -> c^* q";
279 0 : else nameSave = "q q -> b^* q";
280 :
281 : // Locally stored properties and couplings.
282 0 : Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
283 0 : preFac = M_PI / pow4(Lambda);
284 :
285 : // Secondary open width fractions.
286 0 : openFracPos = particleDataPtr->resOpenFrac( idRes);
287 0 : openFracNeg = particleDataPtr->resOpenFrac(-idRes);
288 :
289 0 : }
290 :
291 : //--------------------------------------------------------------------------
292 :
293 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
294 :
295 : void Sigma2qq2qStarq::sigmaKin() {
296 :
297 : // Two possible expressions, for like or unlike sign.
298 0 : sigmaA = preFac * (1. - s3 / sH);
299 0 : sigmaB = preFac * (-uH) * (sH + tH) / sH2;
300 :
301 0 : }
302 :
303 : //--------------------------------------------------------------------------
304 :
305 : // Evaluate sigmaHat(sHat) for specific incoming flavours.
306 :
307 : double Sigma2qq2qStarq::sigmaHat() {
308 :
309 : // Identify different allowed incoming flavour combinations.
310 0 : int id1Abs = abs(id1);
311 0 : int id2Abs = abs(id2);
312 0 : double open1 = (id1 > 0) ? openFracPos : openFracNeg;
313 0 : double open2 = (id2 > 0) ? openFracPos : openFracNeg;
314 : double sigma = 0.;
315 0 : if (id1 * id2 > 0) {
316 0 : if (id1Abs == idq) sigma += (4./3.) * sigmaA * open1;
317 0 : if (id2Abs == idq) sigma += (4./3.) * sigmaA * open2;
318 0 : } else if (id1Abs == idq && id2 == -id1)
319 0 : sigma = (8./3.) * sigmaB * (open1 + open2);
320 0 : else if (id2 == -id1) sigma = sigmaB * (open1 + open2);
321 0 : else if (id1Abs == idq) sigma = sigmaB * open1;
322 0 : else if (id2Abs == idq) sigma = sigmaB * open2;
323 :
324 : // Done.
325 0 : return sigma;
326 :
327 : }
328 :
329 : //--------------------------------------------------------------------------
330 :
331 : // Select identity, colour and anticolour.
332 :
333 : void Sigma2qq2qStarq::setIdColAcol() {
334 :
335 : // Flavours: either side may have been excited.
336 0 : int idAbs1 = abs(id1);
337 0 : int idAbs2 = abs(id2);
338 : double open1 = 0.;
339 : double open2 = 0.;
340 0 : if (idAbs1 == idq) open1 = (id1 > 0) ? openFracPos : openFracNeg;
341 0 : if (idAbs2 == idq) open2 = (id2 > 0) ? openFracPos : openFracNeg;
342 0 : if (open1 == 0. && open2 == 0.) {
343 0 : open1 = (id1 > 0) ? openFracPos : openFracNeg;
344 0 : open2 = (id2 > 0) ? openFracPos : openFracNeg;
345 0 : }
346 0 : bool excite1 = (open1 > 0.);
347 0 : if (open1 > 0. && open2 > 0.) excite1
348 0 : = (rndmPtr->flat() * (open1 + open2) < open1);
349 :
350 : // Always excited quark in slot 3 so colour flow flipped or not.
351 0 : if (excite1) {
352 0 : id3 = (id1 > 0) ? idRes : -idRes;
353 0 : id4 = id2;
354 : // Special case for s-channel like production.
355 0 : if ((idAbs1 == idAbs2) && (id1 * id2 < 0)) {
356 0 : id4 = (id3 > 0) ? -idq : idq;
357 0 : }
358 0 : if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
359 0 : else setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
360 0 : if (id1 < 0) swapColAcol();
361 : } else {
362 0 : id3 = (id2 > 0) ? idRes : -idRes;
363 0 : id4 = id1;
364 : // Special case for s-channel like production.
365 0 : if ((idAbs1 == idAbs2) && (id1 * id2 < 0)) {
366 0 : id4 = (id3 > 0) ? -idq : idq;
367 0 : }
368 0 : swapTU = true;
369 0 : if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
370 0 : else setColAcol( 1, 0, 0, 2, 0, 2, 1, 0);
371 0 : if (id1 < 0) swapColAcol();
372 : }
373 0 : setId( id1, id2, id3, id4);
374 :
375 0 : }
376 :
377 : //--------------------------------------------------------------------------
378 :
379 : // Evaluate weight for q* decay angle.
380 : // SA: Angles dist. for decay q* -> q V, based on Eq. 1.7
381 : // in CERN Yellow Reports 90-10 vol.2, p. 1014 to 1021.
382 :
383 : double Sigma2qq2qStarq::weightDecay( Event& process, int iResBeg,
384 : int iResEnd) {
385 :
386 : // q* should sit in entry 5. Sequential Z/W decay assumed isotropic.
387 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
388 :
389 : // Phase space factors.
390 0 : double mr1 = pow2(process[7].m() / process[5].m());
391 0 : double mr2 = pow2(process[8].m() / process[5].m());
392 :
393 : // Reconstruct decay angle in q* CoM frame.
394 0 : int idAbs3 = process[7].idAbs();
395 : // Olga Igonkina: flip sign of helicity for agreement with Baur et al.
396 : // Vec4 pQStarCom = (idAbs3 < 20) ? process[7].p() : process[8].p();
397 0 : Vec4 pQStarCom = (idAbs3 < 20) ? process[8].p() : process[7].p();
398 0 : pQStarCom.bstback(process[5].p());
399 0 : double cosThe = costheta(pQStarCom, process[5].p());
400 : double wt = 1.;
401 :
402 : // Decay q* -> q (g/gamma) or q (Z^0/W^+-).
403 0 : int idBoson = (idAbs3 < 20) ? process[8].idAbs() : process[7].idAbs();
404 0 : if (idBoson == 21 || idBoson == 22) {
405 0 : wt = 0.5 * (1. + cosThe);
406 0 : } else if (idBoson == 23 || idBoson == 24) {
407 0 : double mrB = (idAbs3 < 20) ? mr2 : mr1;
408 0 : double kTrm = 0.5 * (mrB * (1. - cosThe));
409 0 : wt = (1. + cosThe + kTrm) / (2 + mrB);
410 0 : }
411 :
412 : // Done.
413 : return wt;
414 0 : }
415 :
416 : //==========================================================================
417 :
418 : // Sigma2qqbar2lStarlbar class.
419 : // Cross section for q qbar -> l^* lbar (excited lepton state).
420 :
421 : //--------------------------------------------------------------------------
422 :
423 : // Initialize process.
424 :
425 : void Sigma2qqbar2lStarlbar::initProc() {
426 :
427 : // Set up process properties from the chosen lepton flavour.
428 0 : idRes = 4000000 + idl;
429 0 : codeSave = 4020 + idl;
430 0 : if (idl == 11) nameSave = "q qbar -> e^*+- e^-+";
431 0 : else if (idl == 12) nameSave = "q qbar -> nu_e^* nu_ebar";
432 0 : else if (idl == 13) nameSave = "q qbar -> mu^*+- mu^-+";
433 0 : else if (idl == 14) nameSave = "q qbar -> nu_mu^* nu_mubar";
434 0 : else if (idl == 15) nameSave = "q qbar -> tau^*+- tau^-+";
435 0 : else nameSave = "q qbar -> nu_tau^* nu_taubar";
436 :
437 : // Secondary open width fractions.
438 0 : openFracPos = particleDataPtr->resOpenFrac( idRes);
439 0 : openFracNeg = particleDataPtr->resOpenFrac(-idRes);
440 :
441 : // Locally stored properties and couplings.
442 0 : Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
443 0 : preFac = (M_PI / pow4(Lambda)) * (openFracPos + openFracNeg) / 3.;
444 :
445 0 : }
446 :
447 : //--------------------------------------------------------------------------
448 :
449 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
450 :
451 : void Sigma2qqbar2lStarlbar::sigmaKin() {
452 :
453 : // Only one possible expression.
454 0 : sigma = preFac * (-uH) * (sH + tH) / sH2;
455 :
456 0 : }
457 :
458 : //--------------------------------------------------------------------------
459 :
460 : // Select identity, colour and anticolour.
461 :
462 : void Sigma2qqbar2lStarlbar::setIdColAcol() {
463 :
464 : // Flavours: either lepton or antilepton may be excited.
465 0 : if (rndmPtr->flat() * (openFracPos + openFracNeg) < openFracPos) {
466 0 : setId( id1, id2, idRes, -idl);
467 0 : if (id1 < 0) swapTU = true;
468 : } else {
469 0 : setId( id1, id2, -idRes, idl);
470 0 : if (id1 > 0) swapTU = true;
471 : }
472 :
473 : // Colour flow trivial.
474 0 : if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
475 0 : else setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
476 :
477 0 : }
478 :
479 : //--------------------------------------------------------------------------
480 :
481 : // Evaluate weight for l* decay angle.
482 : // SA: Angles dist. for decay l* -> l V, based on Eq. 1.7
483 : // in CERN Yellow Reports 90-10 vol.2, p. 1014 to 1021.
484 :
485 : double Sigma2qqbar2lStarlbar::weightDecay( Event& process, int iResBeg,
486 : int iResEnd) {
487 :
488 : // l* should sit in entry 5. Sequential Z/W decay assumed isotropic.
489 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
490 :
491 : // Phase space factors.
492 0 : double mr1 = pow2(process[7].m() / process[5].m());
493 0 : double mr2 = pow2(process[8].m() / process[5].m());
494 :
495 : // Reconstruct decay angle in l* CoM frame.
496 0 : int idAbs3 = process[7].idAbs();
497 : // Olga Igonkina: flip sign of helicity for agreement with Baur et al.
498 : // Vec4 pLStarCom = (idAbs3 < 20) ? process[7].p() : process[8].p();
499 0 : Vec4 pLStarCom = (idAbs3 < 20) ? process[8].p() : process[7].p();
500 0 : pLStarCom.bstback(process[5].p());
501 0 : double cosThe = costheta(pLStarCom, process[5].p());
502 : double wt = 1.;
503 :
504 : // Decay, l* -> l + gamma/Z^0/W^+-.
505 0 : int idBoson = (idAbs3 < 20) ? process[8].idAbs() : process[7].idAbs();
506 0 : if (idBoson == 22) {
507 0 : wt = 0.5 * (1. + cosThe);
508 0 : } else if (idBoson == 23 || idBoson == 24) {
509 0 : double mrB = (idAbs3 < 20) ? mr2 : mr1;
510 0 : double kTrm = 0.5 * (mrB * (1. - cosThe));
511 0 : wt = (1. + cosThe + kTrm) / (2 + mrB);
512 0 : }
513 :
514 : // Done.
515 : return wt;
516 0 : }
517 :
518 : //==========================================================================
519 :
520 : // Sigma2qqbar2lStarlStarBar class.
521 : // Cross section for q qbar -> l^* l^*bar (excited lepton state).
522 : // Code contributed by Olga Igonkina.
523 :
524 : //--------------------------------------------------------------------------
525 :
526 : // Initialize process.
527 :
528 : void Sigma2qqbar2lStarlStarBar::initProc() {
529 :
530 : // Set up process properties from the chosen lepton flavour.
531 0 : idRes = 4000000 + idl;
532 0 : codeSave = 4040 + idl;
533 0 : if (idl == 11) nameSave = "q qbar -> e^*+- e^*-+";
534 0 : else if (idl == 12) nameSave = "q qbar -> nu_e^* nu_e^*bar";
535 0 : else if (idl == 13) nameSave = "q qbar -> mu^*+- mu^*-+";
536 0 : else if (idl == 14) nameSave = "q qbar -> nu_mu^* nu_mu^*bar";
537 0 : else if (idl == 15) nameSave = "q qbar -> tau^*+- tau^*-+";
538 0 : else nameSave = "q qbar -> nu_tau^* nu_tau^*bar";
539 :
540 : // Secondary open width fractions.
541 0 : openFracPos = particleDataPtr->resOpenFrac( idRes);
542 0 : openFracNeg = particleDataPtr->resOpenFrac(-idRes);
543 :
544 : // Locally stored properties and couplings.
545 0 : Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
546 0 : preFac = (M_PI / pow4(Lambda)) * openFracPos * openFracNeg / 12.;
547 :
548 0 : }
549 :
550 : //--------------------------------------------------------------------------
551 :
552 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
553 :
554 : void Sigma2qqbar2lStarlStarBar::sigmaKin() {
555 :
556 : // Only one possible expression.
557 0 : sigma = preFac * 2. * (tH2 + uH2 + sH * (s3 + s4) - 2. * s3 * s4) / sH2;
558 :
559 0 : }
560 :
561 : //--------------------------------------------------------------------------
562 :
563 : // Select identity, colour and anticolour.
564 :
565 : void Sigma2qqbar2lStarlStarBar::setIdColAcol() {
566 :
567 : // Flavours: both lepton and antilepton are excited.
568 0 : setId( id1, id2, idRes, -idRes);
569 :
570 : // Colour flow trivial.
571 0 : if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
572 0 : else setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
573 :
574 0 : }
575 :
576 : //--------------------------------------------------------------------------
577 :
578 : // Evaluate weight for l* -> l V decay angles.
579 :
580 : double Sigma2qqbar2lStarlStarBar::weightDecay( Event& process, int iResBeg,
581 : int iResEnd) {
582 :
583 : // l* should sit in 5 and 6. Sequential Z/W decay assumed isotropic.
584 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
585 :
586 : // Loop over two l* decays; check that they are two-body.
587 : double wt = 1.;
588 0 : for (int iResNow = iResBeg; iResNow <= iResEnd; ++iResNow) {
589 0 : int iDau1 = process[iResNow].daughter1();
590 0 : int iDau2 = process[iResNow].daughter2();
591 0 : if (iDau2 != iDau1 + 1) continue;
592 :
593 : // Phase space factors.
594 0 : double mr1 = pow2(process[iDau1].m() / process[iResNow].m());
595 0 : double mr2 = pow2(process[iDau2].m() / process[iResNow].m());
596 :
597 : // Reconstruct decay angle in l* CoM frame.
598 0 : int idAbs3 = process[iDau1].idAbs();
599 0 : Vec4 pLStarCom = (idAbs3 < 20) ? process[iDau2].p() : process[iDau1].p();
600 0 : pLStarCom.bstback(process[iResNow].p());
601 0 : double cosThe = costheta(pLStarCom, process[iResNow].p());
602 :
603 : // Decay, l* -> l + gamma/Z^0/W^+-.
604 0 : int idBoson = (idAbs3 < 20) ? process[iDau2].idAbs()
605 0 : : process[iDau1].idAbs();
606 0 : if (idBoson == 22) {
607 0 : wt *= 0.5 * (1. + cosThe);
608 0 : } else if (idBoson == 23 || idBoson == 24) {
609 0 : double mrB = (idAbs3 < 20) ? mr2 : mr1;
610 0 : double kTrm = 0.5 * (mrB * (1. - cosThe));
611 0 : wt *= (1. + cosThe + kTrm) / (2 + mrB);
612 0 : }
613 0 : }
614 :
615 : // Done.
616 : return wt;
617 :
618 0 : }
619 :
620 : //==========================================================================
621 :
622 : // Sigma2QCqq2qq class.
623 : // Cross section for q q -> q q (quark contact interactions).
624 :
625 : //--------------------------------------------------------------------------
626 :
627 : // Initialize process.
628 :
629 : void Sigma2QCqq2qq::initProc() {
630 :
631 0 : qCLambda2 = settingsPtr->parm("ContactInteractions:Lambda");
632 0 : qCetaLL = settingsPtr->mode("ContactInteractions:etaLL");
633 0 : qCetaRR = settingsPtr->mode("ContactInteractions:etaRR");
634 0 : qCetaLR = settingsPtr->mode("ContactInteractions:etaLR");
635 0 : qCLambda2 *= qCLambda2;
636 :
637 0 : }
638 :
639 : //--------------------------------------------------------------------------
640 :
641 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
642 :
643 : void Sigma2QCqq2qq::sigmaKin() {
644 :
645 : // Calculate kinematics dependence for different terms.
646 0 : sigT = (4./9.) * (sH2 + uH2) / tH2;
647 0 : sigU = (4./9.) * (sH2 + tH2) / uH2;
648 0 : sigTU = - (8./27.) * sH2 / (tH * uH);
649 0 : sigST = - (8./27.) * uH2 / (sH * tH);
650 :
651 0 : sigQCSTU = sH2 * (1 / tH + 1 / uH);
652 0 : sigQCUTS = uH2 * (1 / tH + 1 / sH);
653 :
654 0 : }
655 :
656 : //--------------------------------------------------------------------------
657 :
658 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
659 :
660 : double Sigma2QCqq2qq::sigmaHat() {
661 :
662 : // Terms from QC contact interactions.
663 : double sigQCLL = 0;
664 : double sigQCRR = 0;
665 : double sigQCLR = 0;
666 :
667 : // Combine cross section terms; factor 1/2 when identical quarks.
668 : // q q -> q q
669 0 : if (id2 == id1) {
670 :
671 : // SM terms.
672 0 : sigSum = 0.5 * (sigT + sigU + sigTU);
673 :
674 : // Contact terms.
675 0 : sigQCLL = (8./9.) * alpS * (qCetaLL/qCLambda2) * sigQCSTU
676 0 : + (8./3.) * pow2(qCetaLL/qCLambda2) * sH2;
677 0 : sigQCRR = (8./9.) * alpS * (qCetaRR/qCLambda2) * sigQCSTU
678 0 : + (8./3.) * pow2(qCetaRR/qCLambda2) * sH2;
679 0 : sigQCLR = 2. * (uH2 + tH2) * pow2(qCetaLR/qCLambda2);
680 :
681 0 : sigQCLL /= 2;
682 0 : sigQCRR /= 2;
683 0 : sigQCLR /= 2;
684 :
685 : // q qbar -> q qbar, without pure s-channel term.
686 0 : } else if (id2 == -id1) {
687 :
688 : // SM terms.
689 0 : sigSum = sigT + sigST;
690 :
691 : // Contact terms, minus the terms included in qqbar2qqbar.
692 0 : sigQCLL = (8./9.) * alpS * (qCetaLL/qCLambda2) * sigQCUTS
693 0 : + (5./3.) * pow2(qCetaLL/qCLambda2) * uH2;
694 0 : sigQCRR = (8./9.) * alpS * (qCetaRR/qCLambda2) * sigQCUTS
695 0 : + (5./3.) * pow2(qCetaRR/qCLambda2) * uH2;
696 0 : sigQCLR = 2. * sH2 * pow2(qCetaLR/qCLambda2);
697 :
698 : // q q' -> q q' or q qbar' -> q qbar'
699 0 : } else {
700 :
701 : // SM terms.
702 0 : sigSum = sigT;
703 :
704 : // Contact terms.
705 0 : if (id1 * id2 > 0) {
706 0 : sigQCLL = pow2(qCetaLL/qCLambda2) * sH2;
707 0 : sigQCRR = pow2(qCetaRR/qCLambda2) * sH2;
708 0 : sigQCLR = 2 * pow2(qCetaLR/qCLambda2) * uH2;
709 0 : } else {
710 0 : sigQCLL = pow2(qCetaLL/qCLambda2) * uH2;
711 0 : sigQCRR = pow2(qCetaRR/qCLambda2) * uH2;
712 0 : sigQCLR = 2 * pow2(qCetaLR/qCLambda2) * sH2;
713 : }
714 : }
715 :
716 : // Answer.
717 0 : double sigma = (M_PI/sH2) * ( pow2(alpS) * sigSum
718 0 : + sigQCLL + sigQCRR + sigQCLR );
719 0 : return sigma;
720 :
721 : }
722 :
723 : //--------------------------------------------------------------------------
724 :
725 : // Select identity, colour and anticolour.
726 :
727 : void Sigma2QCqq2qq::setIdColAcol() {
728 :
729 : // Outgoing = incoming flavours.
730 0 : setId( id1, id2, id1, id2);
731 :
732 : // Colour flow topologies. Swap when antiquarks.
733 0 : if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
734 0 : else setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
735 0 : if (id2 == id1 && (sigT + sigU) * rndmPtr->flat() > sigT)
736 0 : setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
737 0 : if (id1 < 0) swapColAcol();
738 :
739 0 : }
740 :
741 : //==========================================================================
742 :
743 : // Sigma2QCqqbar2qqbar class.
744 : // Cross section for q qbar -> q' qbar' (quark contact interactions).
745 :
746 : //--------------------------------------------------------------------------
747 :
748 : // Initialize process.
749 :
750 : void Sigma2QCqqbar2qqbar::initProc() {
751 :
752 0 : qCnQuarkNew = settingsPtr->mode("ContactInteractions:nQuarkNew");
753 0 : qCLambda2 = settingsPtr->parm("ContactInteractions:Lambda");
754 0 : qCetaLL = settingsPtr->mode("ContactInteractions:etaLL");
755 0 : qCetaRR = settingsPtr->mode("ContactInteractions:etaRR");
756 0 : qCetaLR = settingsPtr->mode("ContactInteractions:etaLR");
757 0 : qCLambda2 *= qCLambda2;
758 :
759 0 : }
760 :
761 : //--------------------------------------------------------------------------
762 :
763 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
764 :
765 : void Sigma2QCqqbar2qqbar::sigmaKin() {
766 :
767 : // Pick new flavour.
768 0 : idNew = 1 + int( qCnQuarkNew * rndmPtr->flat() );
769 0 : mNew = particleDataPtr->m0(idNew);
770 0 : m2New = mNew*mNew;
771 :
772 : // Calculate kinematics dependence.
773 : double sigQC = 0.;
774 0 : sigS = 0.;
775 0 : if (sH > 4. * m2New) {
776 0 : sigS = (4./9.) * (tH2 + uH2) / sH2;
777 0 : sigQC = pow2(qCetaLL/qCLambda2) * uH2
778 0 : + pow2(qCetaRR/qCLambda2) * uH2
779 0 : + 2 * pow2(qCetaLR/qCLambda2) * tH2;
780 0 : }
781 :
782 : // Answer is proportional to number of outgoing flavours.
783 0 : sigma = (M_PI / sH2) * qCnQuarkNew * ( pow2(alpS) * sigS + sigQC);
784 :
785 0 : }
786 :
787 : //--------------------------------------------------------------------------
788 :
789 : // Select identity, colour and anticolour.
790 :
791 : void Sigma2QCqqbar2qqbar::setIdColAcol() {
792 :
793 : // Set outgoing flavours ones.
794 0 : id3 = (id1 > 0) ? idNew : -idNew;
795 0 : setId( id1, id2, id3, -id3);
796 :
797 : // Colour flow topologies. Swap when antiquarks.
798 0 : setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
799 0 : if (id1 < 0) swapColAcol();
800 :
801 0 : }
802 :
803 : //==========================================================================
804 :
805 : // Sigma2QCffbar2llbar class.
806 : // Cross section for f fbar -> l lbar (contact interactions).
807 :
808 : //--------------------------------------------------------------------------
809 :
810 : // Initialize process.
811 :
812 : void Sigma2QCffbar2llbar::initProc() {
813 :
814 0 : qCLambda2 = settingsPtr->parm("ContactInteractions:Lambda");
815 0 : qCetaLL = settingsPtr->mode("ContactInteractions:etaLL");
816 0 : qCetaRR = settingsPtr->mode("ContactInteractions:etaRR");
817 0 : qCetaLR = settingsPtr->mode("ContactInteractions:etaLR");
818 0 : qCLambda2 *= qCLambda2;
819 :
820 : // Process name.
821 0 : if (idNew == 11) nameNew = "f fbar -> (QC) -> e- e+";
822 0 : if (idNew == 13) nameNew = "f fbar -> (QC) -> mu- mu+";
823 0 : if (idNew == 15) nameNew = "f fbar -> (QC) -> tau- tau+";
824 :
825 : // Kinematics.
826 0 : qCmNew = particleDataPtr->m0(idNew);
827 0 : qCmNew2 = qCmNew * qCmNew;
828 0 : qCmZ = particleDataPtr->m0(23);
829 0 : qCmZ2 = qCmZ * qCmZ;
830 0 : qCGZ = particleDataPtr->mWidth(23);
831 0 : qCGZ2 = qCGZ * qCGZ;
832 :
833 0 : }
834 :
835 : //--------------------------------------------------------------------------
836 :
837 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
838 :
839 : void Sigma2QCffbar2llbar::sigmaKin() {
840 :
841 0 : qCPropGm = 1./sH;
842 0 : double denomPropZ = pow2(sH - qCmZ2) + qCmZ2 * qCGZ2;
843 0 : qCrePropZ = (sH - qCmZ2) / denomPropZ;
844 0 : qCimPropZ = -qCmZ * qCGZ / denomPropZ;
845 :
846 0 : sigma0 = 0.;
847 0 : if (sH > 4. * qCmNew2) sigma0 = 1./(16. * M_PI * sH2);
848 :
849 0 : }
850 :
851 : //--------------------------------------------------------------------------
852 :
853 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
854 :
855 : double Sigma2QCffbar2llbar::sigmaHat() {
856 :
857 : // Incoming fermion flavor.
858 0 : int idAbs = abs(id1);
859 :
860 : // Couplings and constants.
861 0 : double tmPe2QfQl = 4. * M_PI * alpEM * couplingsPtr->ef(idAbs)
862 0 : * couplingsPtr->ef(idNew);
863 0 : double tmPgvf = 0.25 * couplingsPtr->vf(idAbs);
864 0 : double tmPgaf = 0.25 * couplingsPtr->af(idAbs);
865 0 : double tmPgLf = tmPgvf + tmPgaf;
866 0 : double tmPgRf = tmPgvf - tmPgaf;
867 0 : double tmPgvl = 0.25 * couplingsPtr->vf(idNew);
868 0 : double tmPgal = 0.25 * couplingsPtr->af(idNew);
869 0 : double tmPgLl = tmPgvl + tmPgal;
870 0 : double tmPgRl = tmPgvl - tmPgal;
871 0 : double tmPe2s2c2 = 4. * M_PI * alpEM
872 0 : / (couplingsPtr->sin2thetaW() * couplingsPtr->cos2thetaW());
873 :
874 : // Complex amplitudes.
875 0 : complex I(0., 1.);
876 0 : complex meLL(0., 0.);
877 0 : complex meRR(0., 0.);
878 0 : complex meLR(0., 0.);
879 0 : complex meRL(0., 0.);
880 :
881 : // Amplitudes, M = gamma + Z + CI.
882 0 : meLL = tmPe2QfQl * qCPropGm
883 0 : + tmPe2s2c2 * tmPgLf * tmPgLl * (qCrePropZ + I * qCimPropZ)
884 0 : + 4. * M_PI * qCetaLL / qCLambda2;
885 0 : meRR = tmPe2QfQl * qCPropGm
886 0 : + tmPe2s2c2 * tmPgRf * tmPgRl * (qCrePropZ + I * qCimPropZ)
887 0 : + 4. * M_PI * qCetaRR / qCLambda2;
888 0 : meLR = tmPe2QfQl * qCPropGm
889 0 : + tmPe2s2c2 * tmPgLf * tmPgRl * (qCrePropZ + I * qCimPropZ)
890 0 : + 4. * M_PI * qCetaLR / qCLambda2;
891 0 : meRL = tmPe2QfQl * qCPropGm
892 0 : + tmPe2s2c2 * tmPgRf * tmPgLl * (qCrePropZ + I * qCimPropZ)
893 0 : + 4. * M_PI * qCetaLR / qCLambda2;
894 :
895 0 : double sigma = sigma0 * uH2 * real(meLL*conj(meLL));
896 0 : sigma += sigma0 * uH2 * real(meRR*conj(meRR));
897 0 : sigma += sigma0 * tH2 * real(meLR*conj(meLR));
898 0 : sigma += sigma0 * tH2 * real(meRL*conj(meRL));
899 :
900 : // If f fbar are quarks.
901 0 : if (idAbs < 9) sigma /= 3.;
902 :
903 0 : return sigma;
904 0 : }
905 :
906 : //--------------------------------------------------------------------------
907 :
908 : // Select identity, colour and anticolour.
909 :
910 : void Sigma2QCffbar2llbar::setIdColAcol() {
911 :
912 : // Flavours trivial.
913 0 : setId(id1, id2, idNew, -idNew);
914 :
915 : // tH defined between f and f': must swap tHat <-> uHat if id1 is fbar.
916 0 : swapTU = (id2 > 0);
917 :
918 : // Colour flow topologies. Swap when antiquarks.
919 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
920 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
921 0 : if (id1 < 0) swapColAcol();
922 :
923 0 : }
924 :
925 : //==========================================================================
926 :
927 : } // end namespace Pythia8
|
