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