StRoot  1
1 // SigmaExtraDim.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 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 // Author: Stefan Ask for the *LED* routines.
8 // extra-dimensional simulation classes.
9
11
12 namespace Pythia8 {
13
14 //==========================================================================
15
16 // ampLedS (amplitude) method for LED graviton tree level exchange.
17 // Based on Eq. (8) in JHEP 1105 (2011) 092, arXiv:1101.4919.
18
19 complex ampLedS(double x, double n, double L, double M) {
20
21  complex cS(0., 0.);
22  if (n <= 0) return cS;
23
24  // Constants.
25  double exp1 = n - 2;
26  double exp2 = n + 2;
27  double rC = sqrt(pow(M_PI,n)) * pow(L,exp1)
28  / (GammaReal(n/2.) * pow(M,exp2));
29
30  // Base functions, F1 and F2.
31  complex I(0., 1.);
32  if (x < 0) {
33  double sqrX = sqrt(-x);
34  if (int(n) % 2 == 0) {
35  cS = -log(fabs(1 - 1/x));
36  } else {
37  cS = (2.*atan(sqrX) - M_PI)/sqrX;
38  }
39  } else if ((x > 0) && (x < 1)) {
40  double sqrX = sqrt(x);
41  if (int(n) % 2 == 0) {
42  cS = -log(fabs(1 - 1/x)) - M_PI*I;
43  } else {
44  double rat = (sqrX + 1)/(sqrX - 1);
45  cS = log(fabs(rat))/sqrX - M_PI*I/sqrX;
46  }
47  } else if (x > 1){
48  double sqrX = sqrt(x);
49  if (int(n) % 2 == 0) {
50  cS = -log(fabs(1 - 1/x));
51  } else {
52  double rat = (sqrX + 1)/(sqrX - 1);
53  cS = log(fabs(rat))/sqrX;
54  }
55  }
56
57  // Recursive part.
58  int nL;
59  int nD;
60  if (int(n) % 2 == 0) {
61  nL = int(n/2.);
62  nD = 2;
63  } else {
64  nL = int((n + 1)/2.);
65  nD = 1;
66  }
67  for (int i=1; i<nL; ++i) {
68  cS = x*cS - 2./nD;
69  nD += 2;
70  }
71
72  return rC*cS;
73 }
74
75 //--------------------------------------------------------------------------
76
77 // Common method, "Mandelstam polynomial", for LED dijet processes.
78
79 double funLedG(double x, double y) {
80  double ret = pow(x,4) + 10. * pow(x,3) * y + 42. * pow2(x) * pow2(y)
81  + 64. * x * pow(y,3) + 32. * pow(y,4);
82  return ret;
83 }
84
85 //==========================================================================
86
87 // Sigma1gg2GravitonStar class.
88 // Cross section for g g -> G* (excited graviton state).
89
90 //--------------------------------------------------------------------------
91
92 // Initialize process.
93
94 void Sigma1gg2GravitonStar::initProc() {
95
96  // Store G* mass and width for propagator.
97  idGstar = 5100039;
98  mRes = particleDataPtr->m0(idGstar);
99  GammaRes = particleDataPtr->mWidth(idGstar);
100  m2Res = mRes*mRes;
101  GamMRat = GammaRes / mRes;
102
103  // SMinBulk = off/on, use universal coupling (kappaMG)
104  // or individual (Gxx) between graviton and SM particles.
106  eDvlvl = false;
107  if (eDsmbulk) eDvlvl = settingsPtr->flag("ExtraDimensionsG*:VLVL");
109  for (int i = 0; i < 27; ++i) eDcoupling[i] = 0.;
111  for (int i = 1; i <= 4; ++i) eDcoupling[i] = tmPcoup;
115  for (int i = 11; i <= 16; ++i) eDcoupling[i] = tmPcoup;
121
122  // Set pointer to particle properties and decay table.
123  gStarPtr = particleDataPtr->particleDataEntryPtr(idGstar);
124
125 }
126
127 //--------------------------------------------------------------------------
128
129 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
130
131 void Sigma1gg2GravitonStar::sigmaKin() {
132
133  // Incoming width for gluons.
134  double widthIn = mH / (160. * M_PI);
135
136  // RS graviton coupling
137  if (eDsmbulk) widthIn *= 2. * pow2(eDcoupling[21] * mH);
138  else widthIn *= pow2(kappaMG * mH / mRes);
139
140  // Set up Breit-Wigner. Width out only includes open channels.
141  double sigBW = 5. * M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
142  double widthOut = gStarPtr->resWidthOpen(idGstar, mH);
143
144  // Modify cross section in wings of peak. Done.
145  sigma = widthIn * sigBW * widthOut;
146
147 }
148
149 //--------------------------------------------------------------------------
150
151 // Select identity, colour and anticolour.
152
153 void Sigma1gg2GravitonStar::setIdColAcol() {
154
155  // Flavours trivial.
156  setId( 21, 21, idGstar);
157
158  // Colour flow topology.
159  setColAcol( 1, 2, 2, 1, 0, 0);
160
161 }
162
163 //--------------------------------------------------------------------------
164
165 // Evaluate weight for G* decay angle.
166 // SA: Angle dist. for decay G* -> W/Z/h, based on
167 // Phys.Rev. D65 (2002) 075008, [arXiv:hep-ph/0103308v3]
168
169 double Sigma1gg2GravitonStar::weightDecay( Event& process, int iResBeg,
170  int iResEnd) {
171
172  // Identity of mother of decaying reseonance(s).
173  int idMother = process[process[iResBeg].mother1()].idAbs();
174
175  // For top decay hand over to standard routine.
176  if (idMother == 6)
177  return weightTopDecay( process, iResBeg, iResEnd);
178
179  // G* should sit in entry 5.
180  if (iResBeg != 5 || iResEnd != 5) return 1.;
181
182  // Phase space factors. Reconstruct decay angle.
183  double mr1 = pow2(process[6].m()) / sH;
184  double mr2 = pow2(process[7].m()) / sH;
185  double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
186  double cosThe = (process[3].p() - process[4].p())
187  * (process[7].p() - process[6].p()) / (sH * betaf);
188
189  // Default is isotropic decay.
190  double wt = 1.;
191
192  // Angular weight for g + g -> G* -> f + fbar.
193  if (process[6].idAbs() < 19) {
194  wt = 1. - pow4(cosThe);
195
196  // Angular weight for g + g -> G* -> g + g or gamma + gamma.
197  } else if (process[6].id() == 21 || process[6].id() == 22) {
198  wt = (1. + 6. * pow2(cosThe) + pow4(cosThe)) / 8.;
199
200  // Angular weight for g + g -> G* -> Z + Z or W + W.
201  } else if (process[6].id() == 23 || process[6].id() == 24) {
202  double beta2 = pow2(betaf);
203  double cost2 = pow2(cosThe);
204  double cost4 = pow2(cost2);
205  wt = pow2(beta2 - 2.)*(1. - 2.*cost2 + cost4);
206  // Longitudinal W/Z only.
207  if(eDvlvl) {
208  wt /= 4.;
209  // Transverse W/Z contributions as well.
210  } else {
211  double beta4 = pow2(beta2);
212  double beta8 = pow2(beta4);
213  wt += 2.*pow2(beta4 - 1.)*beta4*cost4;
214  wt += 2.*pow2(beta2 - 1.)*(1. - 2.*beta4*cost2 + beta8*cost4);
215  wt += 2.*(1. + 6.*beta4*cost2 + beta8*cost4);
216  wt += 8.*(1. - beta2)*(1. - cost4);
217  wt /= 18.;
218  }
219
220  // Angular weight for g + g -> G* -> h + h
221  } else if (process[6].id() == 25) {
222  double beta2 = pow2(betaf);
223  double cost2 = pow2(cosThe);
224  double cost4 = pow2(cost2);
225  wt = pow2(beta2 - 2.)*(1. - 2.*cost2 + cost4);
226  wt /= 4.;
227  }
228
229  // Done.
230  return wt;
231
232 }
233
234 //==========================================================================
235
236 // Sigma1ffbar2GravitonStar class.
237 // Cross section for f fbar -> G* (excited graviton state).
238
239 //--------------------------------------------------------------------------
240
241 // Initialize process.
242
243 void Sigma1ffbar2GravitonStar::initProc() {
244
245  // Store G* mass and width for propagator.
246  idGstar = 5100039;
247  mRes = particleDataPtr->m0(idGstar);
248  GammaRes = particleDataPtr->mWidth(idGstar);
249  m2Res = mRes*mRes;
250  GamMRat = GammaRes / mRes;
251
252  // SMinBulk = off/on, use universal coupling (kappaMG)
253  // or individual (Gxx) between graviton and SM particles.
255  eDvlvl = false;
256  if (eDsmbulk) eDvlvl = settingsPtr->flag("ExtraDimensionsG*:VLVL");
258  for (int i = 0; i < 27; ++i) eDcoupling[i] = 0.;
260  for (int i = 1; i <= 4; ++i) eDcoupling[i] = tmPcoup;
264  for (int i = 11; i <= 16; ++i) eDcoupling[i] = tmPcoup;
270
271  // Set pointer to particle properties and decay table.
272  gStarPtr = particleDataPtr->particleDataEntryPtr(idGstar);
273
274 }
275
276 //--------------------------------------------------------------------------
277
278 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
279
280 void Sigma1ffbar2GravitonStar::sigmaKin() {
281
282  // Incoming width for fermions, disregarding colour factor.
283  double widthIn = mH / (80. * M_PI);
284
285  // Set up Breit-Wigner. Width out only includes open channels.
286  double sigBW = 5. * M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
287  double widthOut = gStarPtr->resWidthOpen(idGstar, mH);
288
289  // Modify cross section in wings of peak. Done.
290  sigma0 = widthIn * sigBW * widthOut;
291 }
292
293 //--------------------------------------------------------------------------
294
295 // Evaluate sigmaHat(sHat), part dependent of incoming flavour.
296
297 double Sigma1ffbar2GravitonStar::sigmaHat() {
298
299  double sigma = sigma0;
300
301  // RS graviton coupling
302  if (eDsmbulk) sigma *= 2. * pow2(eDcoupling[min( abs(id1), 26)] * mH);
303  else sigma *= pow2(kappaMG * mH / mRes);
304
305  // If initial quarks, 1/N_C
306  if (abs(id1) < 9) sigma /= 3.;
307
308  return sigma;
309 }
310
311 //--------------------------------------------------------------------------
312
313 // Select identity, colour and anticolour.
314
315 void Sigma1ffbar2GravitonStar::setIdColAcol() {
316
317  // Flavours trivial.
318  setId( id1, id2, idGstar);
319
320  // Colour flow topologies. Swap when antiquarks.
321  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
322  else setColAcol( 0, 0, 0, 0, 0, 0);
323  if (id1 < 0) swapColAcol();
324
325 }
326
327 //--------------------------------------------------------------------------
328
329 // Evaluate weight for G* decay angle.
330 // SA: Angle dist. for decay G* -> W/Z/h, based on
331 // Phys.Rev. D65 (2002) 075008, [arXiv:hep-ph/0103308v3]
332
333 double Sigma1ffbar2GravitonStar::weightDecay( Event& process, int iResBeg,
334  int iResEnd) {
335
336  // Identity of mother of decaying reseonance(s).
337  int idMother = process[process[iResBeg].mother1()].idAbs();
338
339  // For top decay hand over to standard routine.
340  if (idMother == 6)
341  return weightTopDecay( process, iResBeg, iResEnd);
342
343  // G* should sit in entry 5.
344  if (iResBeg != 5 || iResEnd != 5) return 1.;
345
346  // Phase space factors. Reconstruct decay angle.
347  double mr1 = pow2(process[6].m()) / sH;
348  double mr2 = pow2(process[7].m()) / sH;
349  double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
350  double cosThe = (process[3].p() - process[4].p())
351  * (process[7].p() - process[6].p()) / (sH * betaf);
352
353  // Default is isotropic decay.
354  double wt = 1.;
355
356  // Angular weight for f + fbar -> G* -> f + fbar.
357  if (process[6].idAbs() < 19) {
358  wt = (1. - 3. * pow2(cosThe) + 4. * pow4(cosThe)) / 2.;
359
360  // Angular weight for f + fbar -> G* -> g + g or gamma + gamma.
361  } else if (process[6].id() == 21 || process[6].id() == 22) {
362  wt = 1. - pow4(cosThe);
363
364  // Angular weight for f + fbar -> G* -> Z + Z or W + W.
365  } else if (process[6].id() == 23 || process[6].id() == 24) {
366  double beta2 = pow2(betaf);
367  double cost2 = pow2(cosThe);
368  double cost4 = pow2(cost2);
369  wt = pow2(beta2 - 2.)*cost2*(1. - cost2);
370  // Longitudinal W/Z only.
371  if (eDvlvl) {
372  wt /= 4.;
373  // Transverse W/Z contributions as well.
374  } else {
375  wt += pow2(beta2 - 1.)*cost2*(1. - cost2);
376  wt += 2.*(1. - cost4);
377  wt += (1. - beta2)*(1. - 3.*cost2 + 4.*cost4);
378  wt /= 8.;
379  }
380
381  // Angular weight for f + fbar -> G* -> h + h
382  } else if (process[6].id() == 25) {
383  double beta2 = pow2(betaf);
384  double cost2 = pow2(cosThe);
385  wt = pow2(beta2 - 2.)*cost2*(1. - cost2);
386  wt /= 4.;
387  }
388
389  // Done.
390  return wt;
391
392 }
393
394 //==========================================================================
395
396 // Sigma1qqbar2KKgluonStar class.
397 // Cross section for q qbar -> g^*/KK-gluon^* (excited KK-gluon state).
398
399 //--------------------------------------------------------------------------
400
401 // Initialize process.
402
403 void Sigma1qqbar2KKgluonStar::initProc() {
404
405  // Store kk-gluon* mass and width for propagator.
406  idKKgluon = 5100021;
407  mRes = particleDataPtr->m0(idKKgluon);
408  GammaRes = particleDataPtr->mWidth(idKKgluon);
409  m2Res = mRes*mRes;
410  GamMRat = GammaRes / mRes;
411
412  // KK-gluon gv/ga couplings and interference.
413  for (int i = 0; i < 10; ++i) { eDgv[i] = 0.; eDga[i] = 0.; }
416  for (int i = 1; i <= 4; ++i) {
417  eDgv[i] = 0.5 * (tmPgL + tmPgR);
418  eDga[i] = 0.5 * (tmPgL - tmPgR);
419  }
422  eDgv[5] = 0.5 * (tmPgL + tmPgR); eDga[5] = 0.5 * (tmPgL - tmPgR);
425  eDgv[6] = 0.5 * (tmPgL + tmPgR); eDga[6] = 0.5 * (tmPgL - tmPgR);
427
428  // Set pointer to particle properties and decay table.
429  gStarPtr = particleDataPtr->particleDataEntryPtr(idKKgluon);
430
431 }
432
433 //--------------------------------------------------------------------------
434
435 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
436
437 void Sigma1qqbar2KKgluonStar::sigmaKin() {
438
439  // Incoming width for fermions.
440  double widthIn = alpS * mH * 4 / 27;
441  double widthOut = alpS * mH / 6;
442
443  // Loop over all decay channels.
444  sumSM = 0.;
445  sumInt = 0.;
446  sumKK = 0.;
447
448  for (int i = 0; i < gStarPtr->sizeChannels(); ++i) {
449  int idAbs = abs( gStarPtr->channel(i).product(0) );
450
451  // Only contributions quarks.
452  if ( idAbs > 0 && idAbs <= 6 ) {
453  double mf = particleDataPtr->m0(idAbs);
454
455  // Check that above threshold. Phase space.
456  if (mH > 2. * mf + MASSMARGIN) {
457  double mr = pow2(mf / mH);
458  double beta = sqrtpos(1. - 4. * mr);
459
460  // Store sum of combinations. For outstate only open channels.
461  int onMode = gStarPtr->channel(i).onMode();
462  if (onMode == 1 || onMode == 2) {
463  sumSM += beta * (1. + 2. * mr);
464  sumInt += beta * eDgv[min(idAbs, 9)] * (1. + 2. * mr);
465  sumKK += beta * (pow2(eDgv[min(idAbs, 9)]) * (1. + 2.*mr)
466  + pow2(eDga[min(idAbs, 9)]) * (1. - 4.*mr));
467  }
468  }
469  }
470  }
471
472  // Set up Breit-Wigner. Width out only includes open channels.
473  sigSM = widthIn * 12. * M_PI * widthOut / sH2;
474  sigInt = 2. * sigSM * sH * (sH - m2Res)
475  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
476  sigKK = sigSM * sH2 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
477
478  // Optionally only keep g* or gKK term.
479  if (interfMode == 1) {sigInt = 0.; sigKK = 0.;}
480  if (interfMode == 2) {sigSM = 0.; sigInt = 0.;}
481
482 }
483
484 //--------------------------------------------------------------------------
485
486 // Evaluate sigmaHat(sHat), part dependent of incoming flavour.
487
488 double Sigma1qqbar2KKgluonStar::sigmaHat() {
489
490  // RS graviton coupling.
491  double sigma = sigSM * sumSM
492  + eDgv[min(abs(id1), 9)] * sigInt * sumInt
493  + ( pow2(eDgv[min(abs(id1), 9)])
494  + pow2(eDga[min(abs(id1), 9)]) ) * sigKK * sumKK;
495
496  return sigma;
497 }
498
499 //--------------------------------------------------------------------------
500
501 // Select identity, colour and anticolour.
502
503 void Sigma1qqbar2KKgluonStar::setIdColAcol() {
504
505  // Flavours trivial.
506  setId( id1, id2, idKKgluon);
507
508  // Colour flow topologies. Swap when antiquarks.
509  setColAcol( 1, 0, 0, 2, 1, 2);
510  if (id1 < 0) swapColAcol();
511
512 }
513
514 //--------------------------------------------------------------------------
515
516 // Evaluate weight for KK-gluon* decay angle (based on ffbar2gmZ).
517
518 double Sigma1qqbar2KKgluonStar::weightDecay( Event& process, int iResBeg,
519  int iResEnd) {
520
521  // Identity of mother of decaying reseonance(s).
522  int idMother = process[process[iResBeg].mother1()].idAbs();
523
524  // For top decay hand over to standard routine.
525  if (idMother == 6)
526  return weightTopDecay( process, iResBeg, iResEnd);
527
528  // g* should sit in entry 5.
529  if (iResBeg != 5 || iResEnd != 5) return 1.;
530
531  // Couplings for in- and out-flavours (alpS already included).
532  int idInAbs = process[3].idAbs();
533  double vi = eDgv[min(idInAbs, 9)];
534  double ai = eDga[min(idInAbs, 9)];
535  int idOutAbs = process[6].idAbs();
536  double vf = eDgv[min(idOutAbs, 9)];
537  double af = eDga[min(idOutAbs, 9)];
538
539  // Phase space factors. (One power of beta left out in formulae.)
540  double mf = process[6].m();
541  double mr = mf*mf / sH;
542  double betaf = sqrtpos(1. - 4. * mr);
543
544  // Coefficients of angular expression.
545  double coefTran = sigSM + vi * sigInt * vf
546  + (vi*vi + ai*ai) * sigKK * (vf*vf + pow2(betaf) * af*af);
547  double coefLong = 4. * mr * ( sigSM + vi * sigInt * vf
548  + (vi*vi + ai*ai) * sigKK * vf*vf );
549  double coefAsym = betaf * ( ai * sigInt * af
550  + 4. * vi * ai * sigKK * vf * af );
551
552  // Flip asymmetry for in-fermion + out-antifermion.
553  if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
554
555  // Reconstruct decay angle and weight for it.
556  double cosThe = (process[3].p() - process[4].p())
557  * (process[7].p() - process[6].p()) / (sH * betaf);
558  double wtMax = 2. * (coefTran + abs(coefAsym));
559  double wt = coefTran * (1. + pow2(cosThe))
560  + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
561
562  // Done.
563  return (wt / wtMax);
564 }
565
566 //==========================================================================
567
568 // Sigma2gg2GravitonStarg class.
569 // Cross section for g g -> G* g (excited graviton state).
570
571 //--------------------------------------------------------------------------
572
573 // Initialize process.
574
575 void Sigma2gg2GravitonStarg::initProc() {
576
577  // Store G* mass and width for propagator.
578  idGstar = 5100039;
579  mRes = particleDataPtr->m0(idGstar);
580  GammaRes = particleDataPtr->mWidth(idGstar);
581  m2Res = mRes*mRes;
582  GamMRat = GammaRes / mRes;
583
584  // Overall coupling strength kappa * m_G*.
586
587  // Secondary open width fraction.
588  openFrac = particleDataPtr->resOpenFrac(idGstar);
589
590 }
591
592 //--------------------------------------------------------------------------
593
594 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
595
596 void Sigma2gg2GravitonStarg::sigmaKin() {
597
598  // Evaluate cross section. Secondary width for G*.
599  sigma = (3. * pow2(kappaMG) * alpS) / (32. * sH * m2Res)
600  * ( pow2(tH2 + tH * uH + uH2) / (sH2 * tH * uH)
601  + 2. * (tH2 / uH + uH2 / tH) / sH + 3. * (tH / uH + uH / tH)
602  + 2. * (sH / uH + sH/tH) + sH2 / (tH * uH) );
603  sigma *= openFrac;
604
605 }
606
607 //--------------------------------------------------------------------------
608
609 // Select identity, colour and anticolour.
610
611 void Sigma2gg2GravitonStarg::setIdColAcol() {
612
613  // Flavours trivial.
614  setId( 21, 21, idGstar, 21);
615
616  // Colour flow topologies: random choice between two mirrors.
617  if (rndmPtr->flat() < 0.5) setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
618  else setColAcol( 1, 2, 3, 1, 0, 0, 3, 2);
619
620 }
621
622 //--------------------------------------------------------------------------
623
624 // Evaluate weight for decay angles: currently G* assumed isotropic.
625
626 double Sigma2gg2GravitonStarg::weightDecay( Event& process, int iResBeg,
627  int iResEnd) {
628
629  // Identity of mother of decaying reseonance(s).
630  int idMother = process[process[iResBeg].mother1()].idAbs();
631
632  // For top decay hand over to standard routine.
633  if (idMother == 6)
634  return weightTopDecay( process, iResBeg, iResEnd);
635
636  // No equations for G* decay so assume isotropic.
637  return 1.;
638
639 }
640
641 //==========================================================================
642
643 // Sigma2qg2GravitonStarq class.
644 // Cross section for q g -> G* q (excited graviton state).
645
646 //--------------------------------------------------------------------------
647
648 // Initialize process.
649
650 void Sigma2qg2GravitonStarq::initProc() {
651
652  // Store G* mass and width for propagator.
653  idGstar = 5100039;
654  mRes = particleDataPtr->m0(idGstar);
655  GammaRes = particleDataPtr->mWidth(idGstar);
656  m2Res = mRes*mRes;
657  GamMRat = GammaRes / mRes;
658
659  // Overall coupling strength kappa * m_G*.
661
662  // Secondary open width fraction.
663  openFrac = particleDataPtr->resOpenFrac(idGstar);
664
665 }
666
667 //--------------------------------------------------------------------------
668
669 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
670
671 void Sigma2qg2GravitonStarq::sigmaKin() {
672
673  // Evaluate cross section. Secondary width for G*.
674  sigma = -(pow2(kappaMG) * alpS) / (192. * sH * m2Res )
675  * ( 4. * (sH2 + uH2) / (tH * sH) + 9. * (sH + uH) / sH + sH / uH
676  + uH2 / sH2 + 3. * tH * (4. + sH / uH + uH / sH) / sH
677  + 4. * tH2 * (1. / uH + 1. / sH) / sH + 2. * tH2 * tH / (uH * sH2) );
678  sigma *= openFrac;
679
680 }
681
682 //--------------------------------------------------------------------------
683
684 // Select identity, colour and anticolour.
685
686 void Sigma2qg2GravitonStarq::setIdColAcol() {
687
688  // Flavour set up for q g -> H q.
689  int idq = (id2 == 21) ? id1 : id2;
690  setId( id1, id2, idGstar, idq);
691
692  // tH defined between f and f': must swap tHat <-> uHat if q g in.
693  swapTU = (id2 == 21);
694
695  // Colour flow topologies. Swap when antiquarks.
696  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
697  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
698  if (idq < 0) swapColAcol();
699
700 }
701
702 //--------------------------------------------------------------------------
703
704 // Evaluate weight for decay angles: currently G* assumed isotropic.
705
706 double Sigma2qg2GravitonStarq::weightDecay( Event& process, int iResBeg,
707  int iResEnd) {
708
709  // Identity of mother of decaying reseonance(s).
710  int idMother = process[process[iResBeg].mother1()].idAbs();
711
712  // For top decay hand over to standard routine.
713  if (idMother == 6)
714  return weightTopDecay( process, iResBeg, iResEnd);
715
716  // No equations for G* decay so assume isotropic.
717  return 1.;
718
719 }
720
721 //==========================================================================
722
723 // Sigma2qqbar2GravitonStarg class.
724 // Cross section for q qbar -> G* g (excited graviton state).
725
726 //--------------------------------------------------------------------------
727
728 // Initialize process.
729
730 void Sigma2qqbar2GravitonStarg::initProc() {
731
732  // Store G* mass and width for propagator.
733  idGstar = 5100039;
734  mRes = particleDataPtr->m0(idGstar);
735  GammaRes = particleDataPtr->mWidth(idGstar);
736  m2Res = mRes*mRes;
737  GamMRat = GammaRes / mRes;
738
739  // Overall coupling strength kappa * m_G*.
741
742  // Secondary open width fraction.
743  openFrac = particleDataPtr->resOpenFrac(idGstar);
744
745 }
746
747 //--------------------------------------------------------------------------
748
749 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
750
751 void Sigma2qqbar2GravitonStarg::sigmaKin() {
752
753  // Evaluate cross section. Secondary width for G*.
754  sigma = (pow2(kappaMG) * alpS) / (72. * sH * m2Res)
755  * ( 4. * (tH2 + uH2) / sH2 + 9. * (tH + uH) / sH
756  + (tH2 / uH + uH2 / tH) / sH + 3. * (4. + tH / uH + uH/ tH)
757  + 4. * (sH / uH + sH / tH) + 2. * sH2 / (tH * uH) );
758  sigma *= openFrac;
759
760 }
761
762 //--------------------------------------------------------------------------
763
764 // Select identity, colour and anticolour.
765
766 void Sigma2qqbar2GravitonStarg::setIdColAcol() {
767
768  // Flavours trivial.
769  setId( id1, id2, idGstar, 21);
770
771  // Colour flow topologies. Swap when antiquarks.
772  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
773  if (id1 < 0) swapColAcol();
774
775 }
776
777 //--------------------------------------------------------------------------
778
779 // Evaluate weight for decay angles: currently G* assumed isotropic.
780
781 double Sigma2qqbar2GravitonStarg::weightDecay( Event& process, int iResBeg,
782  int iResEnd) {
783
784  // Identity of mother of decaying reseonance(s).
785  int idMother = process[process[iResBeg].mother1()].idAbs();
786
787  // For top decay hand over to standard routine.
788  if (idMother == 6)
789  return weightTopDecay( process, iResBeg, iResEnd);
790
791  // No equations for G* decay so assume isotropic.
792  return 1.;
793
794 }
795
796 //==========================================================================
797
798 // NOAM: Sigma2ffbar2TEVffbar class.
799 // Cross section for, f fbar -> gammaKK/ZKK -> F Fbar.
800 // Process provided by N. Hod et al. and is described in arXiv:XXXX.YYYY
801
802 //--------------------------------------------------------------------------
803
804 // Initialize process.
805
806 void Sigma2ffbar2TEVffbar::initProc() {
807
808  // Process name.
809  if (idNew == 1) nameSave = "f fbar -> d dbar (s-channel gamma_KK/Z_KK)";
810  if (idNew == 2) nameSave = "f fbar -> u ubar (s-channel gamma_KK/Z_KK)";
811  if (idNew == 3) nameSave = "f fbar -> s sbar (s-channel gamma_KK/Z_KK)";
812  if (idNew == 4) nameSave = "f fbar -> c cbar (s-channel gamma_KK/Z_KK)";
813  if (idNew == 5) nameSave = "f fbar -> b bbar (s-channel gamma_KK/Z_KK)";
814  if (idNew == 6) nameSave = "f fbar -> t tbar (s-channel gamma_KK/Z_KK)";
815  if (idNew == 11) nameSave = "f fbar -> e+ e- (s-channel gamma_KK/Z_KK)";
816  if (idNew == 12) nameSave = "f fbar -> nue nuebar (s-channel gamma_KK/Z_KK)";
817  if (idNew == 13) nameSave = "f fbar -> mu+ mu- (s-channel gamma_KK/Z_KK)";
818  if (idNew == 14) nameSave
819  = "f fbar -> numu numubar (s-channel gamma_KK/Z_KK)";
820  if (idNew == 15) nameSave = "f fbar -> tau+ tau- (s-channel gamma_KK/Z_KK)";
821  if (idNew == 16) nameSave
822  = "f fbar -> nutau nutaubar (s-channel gamma_KK/Z_KK)";
823
824  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
826
827  // Pick number of KK excitations
829
830  // Initialize the widths of the KK propogators.
831  // partial width of the KK photon
832  wgmKKFactor = 0.;
833  // total width of the KK photon
834  wgmKKn = 0.;
835  // will be proportional to "wZ0" + ttbar addition
836  wZKKn = 0.;
837
838  // Store Z0 mass and width for propagator.
839  wZ0 = particleDataPtr->mWidth(23);
840  mRes = particleDataPtr->m0(23);
841  m2Res = mRes*mRes;
842
843  // Store the top mass for the ttbar width calculations
844  mTop = particleDataPtr->m0(6);
845  m2Top = mTop*mTop;
846
847  // Store the KK mass parameter, equivalent to the mass of the first KK
848  // excitation: particleDataPtr->m0(5000023);
850
851  // Get alphaEM - relevant for the calculation of the widths
852  alphaemfixed = settingsPtr->parm("StandardModel:alphaEM0");
853
854  // initialize imaginari number
855  mI = complex(0.,1.);
856
857  // Sum all partial widths of the KK photon except for the ttbar channel
858  // which is handeled afterwards seperately
859  if (gmZmode>=0 && gmZmode<=5) {
860  for (int i=1 ; i<17 ; i++) {
861  if (i==7) { i=11; }
862  // skip the ttbar decay and add its contribution later
863  if (i==6) { continue; }
864  if (i<9) {
865  wgmKKFactor += ( (alphaemfixed / 6.) * 4.
866  * couplingsPtr->ef(i) * couplingsPtr->ef(i) * 3. );
867  }
868  else {
869  wgmKKFactor += (alphaemfixed / 6.) * 4.
870  * couplingsPtr->ef(i) * couplingsPtr->ef(i);
871  }
872  }
873  }
874
875  // Get the helicity-couplings of the Z0 to all the fermions except top
876  gMinusF = ( couplingsPtr->t3f(idNew) - couplingsPtr->ef(idNew)
877  * couplingsPtr->sin2thetaW() )
878  / sqrt( couplingsPtr->sin2thetaW()*couplingsPtr->cos2thetaW() );
879  gPlusF = -1. * couplingsPtr->ef(idNew) * couplingsPtr->sin2thetaW()
880  / sqrt( couplingsPtr->sin2thetaW() * couplingsPtr->cos2thetaW() );
881  // Get the helicity-couplings of the Z0 to the top quark
882  gMinusTop = ( couplingsPtr->t3f(6) - couplingsPtr->ef(6)
883  * couplingsPtr->sin2thetaW() )
884  / sqrt( couplingsPtr->sin2thetaW()*couplingsPtr->cos2thetaW() );
885
886  gPlusTop = -1. * couplingsPtr->ef(6) * couplingsPtr->sin2thetaW()
887  / sqrt( couplingsPtr->sin2thetaW() * couplingsPtr->cos2thetaW() );
888  // calculate the constant factor of the unique ttbar decay width
889  ttbarwFactorA = pow2(gMinusTop) + pow2(gPlusTop);
890  ttbarwFactorB = 6.*gMinusTop*gPlusTop - pow2(gMinusTop) - pow2(gPlusTop);
891
892  // Secondary open width fraction, relevant for top (or heavier).
893  openFracPair = 1.;
894  if ((idNew >=6 && idNew <=8) || idNew == 17 || idNew == 18)
895  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
896
897 }
898
899 //--------------------------------------------------------------------------
900
901 // For improving the phase-space sampling (there can be 2 resonances)
902
903 int Sigma2ffbar2TEVffbar::resonanceB() {
904
905  return 23;
906
907 }
908
909 //--------------------------------------------------------------------------
910
911 // For improving the phase-space sampling (there can be 2 resonances)
912
913 int Sigma2ffbar2TEVffbar::resonanceA() {
914
915  if (gmZmode>=3) {
916  phaseSpacemHatMin = settingsPtr->parm("PhaseSpace:mHatMin");
917  phaseSpacemHatMax = settingsPtr->parm("PhaseSpace:mHatMax");
918  double mResFirstKKMode = sqrt(pow2(particleDataPtr->m0(23)) + pow2(mStar));
919  if (mResFirstKKMode/2. <= phaseSpacemHatMax
920  || 3*mResFirstKKMode/2. >= phaseSpacemHatMin) { return 5000023; }
921  else { return 23; }
922  // no KK terms at all
923  } else { return 23; }
924
925 }
926
927 //--------------------------------------------------------------------------
928
929 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
930
931 void Sigma2ffbar2TEVffbar::sigmaKin() {
932
933  // Check that above threshold.
934  isPhysical = true;
935  if (mH < m3 + m4 + MASSMARGIN) {
936  isPhysical = false;
937  return;
938  }
939
940  // Define average F, Fbar mass so same beta. Phase space.
941  double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
942  mr = s34Avg / sH;
943  betaf = sqrtpos(1. - 4. * mr);
944
945  // Reconstruct decay angle so can reuse 2 -> 1 cross section.
946  cosThe = (tH - uH) / (betaf * sH);
947
948 }
949
950 //--------------------------------------------------------------------------
951
952 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
953
954 double Sigma2ffbar2TEVffbar::sigmaHat() {
955
956  // Fail if below threshold.
957  if (!isPhysical) return 0.;
958
959  // Couplings for in/out-flavours.
960  int idAbs = abs(id1);
961
962  // The couplings of the Z0 to the fermions for in/out flavors
963  gMinusf = ( couplingsPtr->t3f(idAbs) - couplingsPtr->ef(idAbs)
964  * couplingsPtr->sin2thetaW() )
965  / sqrt( couplingsPtr->sin2thetaW()*couplingsPtr->cos2thetaW() );
966  gPlusf = -1. * couplingsPtr->ef(idAbs)*couplingsPtr->sin2thetaW()
967  / sqrt( couplingsPtr->sin2thetaW()*couplingsPtr->cos2thetaW() );
968
969  // Initialize the some values
970  helicityME2 = 0.;
971  coefAngular = 0.;
972  gf=0.;
973  gF=0.;
974  gammaProp = complex(0.,0.);
975  resProp = complex(0.,0.);
976  gmPropKK = complex(0.,0.);
977  ZPropKK = complex(0.,0.);
978  totalProp = complex(0.,0.);
979
980  // Sum all initial and final helicity states this corresponds to an
981  // unpolarized beams and unmeasured polarization final-state
982  for (double helicityf=-0.5 ; helicityf<=0.5 ; helicityf++) {
983  for (double helicityF=-0.5 ; helicityF<=0.5 ; helicityF++) {
984  // the couplings for the initial-final helicity configuration
985  gF = (helicityF == +0.5) ? gMinusF : gPlusF;
986  gf = (helicityf == +0.5) ? gMinusf : gPlusf;
987  // 0=SM gmZ, 1=SM gm, 2=SM Z, 3=SM+KK gmZ, 4=KK gm, 5=KK Z
988  switch(gmZmode) {
989  // SM photon and Z0 only
990  case 0:
991  gammaProp = couplingsPtr->ef(idAbs)*couplingsPtr->ef(idNew)/sH;
992  resProp = gf*gF/( sH - m2Res + mI*sH*(wZ0/mRes) );
993  break;
994  // SM photon only
995  case 1:
996  gammaProp = couplingsPtr->ef(idAbs)*couplingsPtr->ef(idNew)/sH;
997  break;
998  // SM Z0 only
999  case 2:
1000  resProp = gf*gF/( sH - m2Res + mI*sH*(wZ0/mRes) );
1001  break;
1002  // KK photon and Z
1003  case 3:
1004  gammaProp = couplingsPtr->ef(idAbs)*couplingsPtr->ef(idNew)/sH;
1005  resProp = gf*gF/( sH - m2Res + mI*sH*(wZ0/mRes) );
1006  ZPropKK = complex(0.,0.);
1007  gmPropKK = complex(0.,0.);
1008  // Sum all KK excitations contributions
1009  for(int nexcitation = 1; nexcitation <= nexcitationmax;
1010  nexcitation++) {
1011  mZKKn = sqrt(m2Res + pow2(mStar * nexcitation));
1012  m2ZKKn = m2Res + pow2(mStar * nexcitation);
1013  mgmKKn = mStar * nexcitation;
1014  m2gmKKn = (mStar*nexcitation)*(mStar*nexcitation);
1015  // calculate the width of the n'th excitation of the KK Z
1016  // (proportional to the Z0 width + ttbar partial width)
1017  ttbarwZKKn = 2.*(alphaemfixed*3./6.)*mZKKn
1018  * sqrt(1.-4.*m2Top/m2ZKKn)
1019  * (ttbarwFactorA+(m2Top/m2ZKKn)*ttbarwFactorB);
1020  wZKKn = 2.*wZ0*mZKKn/mRes+ttbarwZKKn;
1021  // calculate the width of the n'th excitation of the
1022  // KK photon
1023  ttbarwgmKKn = 2.*(alphaemfixed*3./6.)*mgmKKn
1024  * sqrt(1.-4.*m2Top/m2gmKKn)
1025  * 2.*pow2(couplingsPtr->ef(6))*(1.+2.*(m2Top/m2gmKKn));
1026  wgmKKn = wgmKKFactor*mgmKKn+ttbarwgmKKn;
1027  // the propogators
1028  gmPropKK += (2.*couplingsPtr->ef(idAbs)*couplingsPtr->ef(idNew))
1029  / (sH-m2gmKKn+mI*sH*wgmKKn/mgmKKn);
1030  ZPropKK += (2.*gf*gF)/(sH-m2ZKKn+mI*sH*wZKKn/mZKKn );
1031  }
1032  break;
1033  // SM photon and Z0 with KK photon only
1034  case 4:
1035  gammaProp = couplingsPtr->ef(idAbs)*couplingsPtr->ef(idNew)/sH;
1036  resProp = gf*gF/( sH - m2Res + mI*sH*(wZ0/mRes) );
1037  gmPropKK = complex(0.,0.);
1038  for (int nexcitation = 1; nexcitation <= nexcitationmax;
1039  nexcitation++ ) {
1040  mgmKKn = mStar * nexcitation;
1041  m2gmKKn = (mStar*nexcitation)*(mStar*nexcitation);
1042
1043  ttbarwgmKKn = 2.*(alphaemfixed*3./6.)*mgmKKn
1044  * sqrt(1.-4.*m2Top/m2gmKKn)
1045  * 2.*pow2(couplingsPtr->ef(6))
1046  * (1.+2.*(m2Top/m2gmKKn));
1047  wgmKKn = wgmKKFactor*mgmKKn+ttbarwgmKKn;
1048  gmPropKK += (2.*couplingsPtr->ef(idAbs)*couplingsPtr->ef(idNew))
1049  / (sH-m2gmKKn+mI*sH*wgmKKn/mgmKKn);
1050  }
1051  break;
1052  // SM photon and Z0 with KK Z only
1053  case 5:
1054  gammaProp = couplingsPtr->ef(idAbs)*couplingsPtr->ef(idNew)/sH;
1055  resProp = gf*gF/( sH - m2Res + mI*sH*(wZ0/mRes) );
1056  ZPropKK = complex(0.,0.);
1057  for (int nexcitation = 1; nexcitation <= nexcitationmax;
1058  nexcitation++ ) {
1059  mZKKn = sqrt(m2Res + pow2(mStar * nexcitation));
1060  m2ZKKn = m2Res + pow2(mStar * nexcitation);
1061
1062  ttbarwZKKn = 2.*(alphaemfixed*3./6.)*mZKKn
1063  * sqrt(1.-4.*m2Top/m2ZKKn)
1064  * (ttbarwFactorA+(m2Top/m2ZKKn)*ttbarwFactorB);
1065  wZKKn = 2.*wZ0*mZKKn/mRes+ttbarwZKKn;
1066  ZPropKK += (2.*gf*gF)/(sH-m2ZKKn+mI*sH*wZKKn/mZKKn );
1067  }
1068  break;
1069  default: break;
1070  // end run over initial and final helicity states
1071  }
1072
1073  // sum all contributing amplitudes
1074  totalProp = gammaProp + resProp + ZPropKK + gmPropKK;
1075
1076  // angular distribution for the helicity configuration
1077  coefAngular = 1. + 4. * helicityF * helicityf * cosThe;
1078
1079  // the squared helicity matrix element
1080  helicityME2 += real(totalProp*conj(totalProp))*pow2(coefAngular);
1081  }
1082  }
1083
1084  // calculate the coefficient of the squared helicity matrix element.
1085  coefTot = (2./sH) * 2*M_PI * pow2(alpEM)/(4.*sH) * pow2(sH)/4.;
1086
1087  // the full squared helicity matrix element.
1088  double sigma = helicityME2 * coefTot;
1089
1090  // Top: corrections for closed decay channels.
1091  sigma *= openFracPair;
1092
1093  // Initial-state colour factor. Answer.
1094  if (idAbs < 9) sigma /= 3.;
1095
1096  // Final-state colour factor. Answer.
1097  if (idNew < 9) sigma *= 3.*(1.+alpS/M_PI);
1098
1099  return sigma;
1100 }
1101
1102 //--------------------------------------------------------------------------
1103
1104 // Select identity, colour and anticolour.
1105
1106 void Sigma2ffbar2TEVffbar::setIdColAcol() {
1107
1108  // Set outgoing flavours.
1109  id3 = (id1 > 0) ? idNew : -idNew;
1110  setId( id1, id2, id3, -id3);
1111
1112  // Colour flow topologies. Swap when antiquarks.
1113  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1114  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1115  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1116  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1117  if (id1 < 0) swapColAcol();
1118
1119 }
1120
1121 //--------------------------------------------------------------------------
1122
1123 // Evaluate weight for decay angles of W in top decay.
1124
1125 double Sigma2ffbar2TEVffbar::weightDecay( Event& process, int iResBeg,
1126  int iResEnd) {
1127
1128  // For top decay hand over to standard routine, else done.
1129  if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1130  return weightTopDecay( process, iResBeg, iResEnd);
1131  else return 1.;
1132
1133 }
1134
1135 //==========================================================================
1136
1137 // Sigma2gg2LEDUnparticleg class.
1138 // Cross section for g g -> U/G g (real graviton emission in
1139 // large extra dimensions or unparticle emission).
1140
1141 //--------------------------------------------------------------------------
1142
1143 void Sigma2gg2LEDUnparticleg::initProc() {
1144
1145  // Init model parameters.
1146  eDidG = 5000039;
1147  if (eDgraviton) {
1148  eDspin = (settingsPtr->flag("ExtraDimensionsLED:GravScalar")) ? 0 : 2;
1150  eDdU = 0.5 * eDnGrav + 1;
1152  eDlambda = 1;
1156  } else {
1162  }
1163
1164  // The A(dU) or S'(n) value.
1166  if (eDgraviton) {
1167  tmpAdU = 2 * M_PI * sqrt( pow(M_PI, double(eDnGrav)) )
1168  / GammaReal(0.5 * eDnGrav);
1169  if (eDspin == 0) { // Scalar graviton
1170  tmpAdU *= sqrt( pow(2., double(eDnGrav)) );
1171  eDcf *= eDcf;
1172  }
1173  } else {
1174  tmpAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
1175  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
1176  }
1177
1178  // Cross section related constants
1179  // and ME dependent powers of lambda / LambdaU.
1180  double tmpExp = eDdU - 2;
1181  double tmpLS = pow2(eDLambdaU);
1182  eDconstantTerm = tmpAdU / (2 * 16 * pow2(M_PI) * tmpLS * pow(tmpLS,tmpExp));
1183  if (eDgraviton) {
1184  eDconstantTerm /= tmpLS;
1185  } else if (eDspin == 0) {
1186  eDconstantTerm *= pow2(eDlambda) / tmpLS;
1187  } else {
1188  eDconstantTerm = 0;
1189  infoPtr->errorMsg("Error in Sigma2gg2LEDUnparticleg::initProc: "
1190  "Incorrect spin value (turn process off)!");
1191  }
1192
1193 }
1194
1195 //--------------------------------------------------------------------------
1196
1197 void Sigma2gg2LEDUnparticleg::sigmaKin() {
1198
1199  // Set graviton mass.
1200  mG = m3;
1201  mGS = mG*mG;
1202
1203  // Set mandelstam variables and ME expressions.
1204  if (eDgraviton) {
1205
1206  double A0 = 1/sH;
1207  if (eDspin == 0) { // Scalar graviton
1208  double tmpTerm1 = uH + tH;
1209  double tmpTerm2 = uH + sH;
1210  double tmpTerm3 = tH + sH;
1211  double T0 = pow(tmpTerm1,4) + pow(tmpTerm2,4) + pow(tmpTerm3,4)
1212  + 12. * sH * tH * uH * mGS;
1213  eDsigma0 = eDcf * A0 * T0 / (sH2 * tH * uH);
1214  } else {
1215  double xH = tH/sH;
1216  double yH = mGS/sH;
1217  double xHS = pow2(xH);
1218  double yHS = pow2(yH);
1219  double xHC = pow(xH,3);
1220  double yHC = pow(yH,3);
1221  double xHQ = pow(xH,4);
1222  double yHQ = pow(yH,4);
1223
1224  double T0 = 1/(xH*(yH-1-xH));
1225  double T1 = 1 + 2*xH + 3*xHS + 2*xHC + xHQ;
1226  double T2 = -2*yH*(1 + xHC);
1227  double T3 = 3*yHS*(1 + xHS);
1228  double T4 = -2*yHC*(1 + xH);
1229  double T5 = yHQ;
1230
1231  eDsigma0 = A0 * T0 *( T1 + T2 + T3 + T4 + T5 );
1232  }
1233
1234  } else if (eDspin == 0) {
1235
1236  double A0 = 1/pow2(sH);
1237  double sHQ = pow(sH,4);
1238  double tHQ = pow(tH,4);
1239  double uHQ = pow(uH,4);
1240
1241  eDsigma0 = A0 * (pow(mGS,4) + sHQ + tHQ + uHQ) / (sH * tH * uH);
1242
1243  }
1244
1245  // Mass measure, (m^2)^(d-2).
1246  double tmpExp = eDdU - 2;
1247  eDsigma0 *= pow(mGS, tmpExp);
1248
1249  // Constants.
1250  eDsigma0 *= eDconstantTerm;
1251
1252 }
1253
1254 //--------------------------------------------------------------------------
1255
1256 double Sigma2gg2LEDUnparticleg::sigmaHat() {
1257
1258  // Mass spectrum weighting.
1259  double sigma = eDsigma0 / runBW3;
1260
1261  // SM couplings...
1262  if (eDgraviton) {
1263  sigma *= 16 * M_PI * alpS * 3 / 16;
1264  } else if (eDspin == 0) {
1265  sigma *= 6 * M_PI * alpS;
1266  }
1267
1268  // Truncate sH region or use form factor.
1269  // Form factor uses either pythia8 renormScale2
1270  // or E_jet in cms.
1271  if (eDcutoff == 1) {
1272  if (sH > pow2(eDLambdaU) ) { sigma *= pow(eDLambdaU,4)/pow2(sH); }
1273  } else if ( (eDgraviton && (eDspin == 2))
1274  && ((eDcutoff == 2) || (eDcutoff == 3)) ) {
1275  double tmPmu = sqrt(Q2RenSave);
1276  if (eDcutoff == 3) tmPmu = (sH + s4 - s3) / (2 * mH);
1277  double tmPformfact = tmPmu / (eDtff * eDLambdaU);
1278  double tmPexp = double(eDnGrav) + 2;
1279  sigma *= 1 / (1 + pow(tmPformfact, tmPexp));
1280  }
1281
1282  return sigma;
1283 }
1284
1285 //--------------------------------------------------------------------------
1286
1287 void Sigma2gg2LEDUnparticleg::setIdColAcol() {
1288
1289  // Flavours trivial.
1290  setId( 21, 21, eDidG, 21);
1291
1292  // Colour flow topologies: random choice between two mirrors.
1293  if (rndmPtr->flat() < 0.5) setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
1294  else setColAcol( 1, 2, 3, 1, 0, 0, 3, 2);
1295
1296 }
1297
1298 //==========================================================================
1299
1300 // Sigma2qg2LEDUnparticleq class.
1301 // Cross section for q g -> U/G q (real graviton emission in
1302 // large extra dimensions or unparticle emission).
1303
1304 //--------------------------------------------------------------------------
1305
1306 void Sigma2qg2LEDUnparticleq::initProc() {
1307
1308  // Init model parameters.
1309  eDidG = 5000039;
1310  if (eDgraviton) {
1311  eDspin = (settingsPtr->flag("ExtraDimensionsLED:GravScalar")) ? 0 : 2;
1313  eDdU = 0.5 * eDnGrav + 1;
1315  eDlambda = 1;
1320  } else {
1326  }
1327
1328  // The A(dU) or S'(n) value.
1330  if (eDgraviton) {
1331  tmpAdU = 2 * M_PI * sqrt( pow(M_PI, double(eDnGrav)) )
1332  / GammaReal(0.5 * eDnGrav);
1333  // Scalar graviton
1334  if (eDspin == 0) {
1335  tmpAdU *= 2. * sqrt( pow(2., double(eDnGrav)) );
1336  eDcf *= 4. * eDcf / pow2(eDLambdaU);
1337  double tmpExp = 2. * double(eDnGrav) / (double(eDnGrav) + 2.);
1338  eDgf *= eDgf / pow(2. * M_PI, tmpExp);
1339  }
1340  } else {
1341  tmpAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
1342  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
1343  }
1344
1345  // Cross section related constants
1346  // and ME dependent powers of lambda / LambdaU.
1347  double tmpExp = eDdU - 2;
1348  double tmpLS = pow2(eDLambdaU);
1349  eDconstantTerm = tmpAdU / (2 * 16 * pow2(M_PI) * tmpLS * pow(tmpLS,tmpExp));
1350  if (eDgraviton && (eDspin == 2)) {
1351  eDconstantTerm /= tmpLS;
1352  } else if (eDspin == 1) {
1353  eDconstantTerm *= pow2(eDlambda);
1354  } else if (eDspin == 0) {
1355  eDconstantTerm *= pow2(eDlambda);
1356  } else {
1357  eDconstantTerm = 0;
1358  infoPtr->errorMsg("Error in Sigma2qg2LEDUnparticleq::initProc: "
1359  "Incorrect spin value (turn process off)!");
1360  }
1361
1362
1363 }
1364
1365 //--------------------------------------------------------------------------
1366
1367 void Sigma2qg2LEDUnparticleq::sigmaKin() {
1368
1369  // Set graviton mass.
1370  mG = m3;
1371  mGS = mG*mG;
1372
1373  // Set mandelstam variables and ME expressions.
1374  if (eDgraviton) {
1375
1376  double A0 = 1/sH;
1377  // Scalar graviton
1378  if (eDspin == 0) {
1379  A0 /= sH;
1380  double T0 = -(uH2 + pow2(mGS)) / (sH * tH);
1381  double T1 = -(tH2 + sH2)/ uH;
1382  eDsigma0 = A0 * (eDgf * T0 + eDcf * T1);
1383  } else {
1384  double xH = tH/sH;
1385  double yH = mGS/sH;
1386  double x2H = xH/(yH - 1 - xH);
1387  double y2H = yH/(yH - 1 - xH);
1388  double x2HS = pow2(x2H);
1389  double y2HS = pow2(y2H);
1390  double x2HC = pow(x2H,3);
1391  double y2HC = pow(y2H,3);
1392
1393  double T0 = -(yH - 1 - xH);
1394  double T20 = 1/(x2H*(y2H-1-x2H));
1395  double T21 = -4*x2H*(1 + x2H)*(1 + 2*x2H + 2*x2HS);
1396  double T22 = y2H*(1 + 6*x2H + 18*x2HS + 16*x2HC);
1397  double T23 = -6*y2HS*x2H*(1+2*x2H);
1398  double T24 = y2HC*(1 + 4*x2H);
1399
1400  eDsigma0 = A0 * T0 * T20 * ( T21 + T22 + T23 + T24 );
1401  }
1402
1403  } else if (eDspin == 1) {
1404
1405  double A0 = 1/pow2(sH);
1406  double tmpTerm1 = tH - mGS;
1407  double tmpTerm2 = sH - mGS;
1408
1409  eDsigma0 = A0 * (pow2(tmpTerm1) + pow2(tmpTerm2)) / (sH*tH);
1410
1411  } else if (eDspin == 0) {
1412
1413  double A0 = 1/pow2(sH);
1414  // Sign correction by Tom
1415  eDsigma0 = A0 * (pow2(tH) + pow2(mGS)) / (sH*uH);
1416
1417  }
1418
1419  // Mass measure, (m^2)^(d-2).
1420  double tmpExp = eDdU - 2;
1421  eDsigma0 *= pow(mGS, tmpExp);
1422
1423  // Constants.
1424  eDsigma0 *= eDconstantTerm;
1425
1426 }
1427
1428 //--------------------------------------------------------------------------
1429
1430 double Sigma2qg2LEDUnparticleq::sigmaHat() {
1431
1432  // Mass spactrum weighting.
1433  double sigma = eDsigma0 /runBW3;
1434
1435  // SM couplings...
1436  if (eDgraviton) {
1437  sigma *= 16 * M_PI * alpS / 96;
1438  } else if (eDspin == 1) {
1439  sigma *= - 4 * M_PI * alpS / 3;
1440  } else if (eDspin == 0) {
1441  sigma *= - 2 * M_PI * alpS / 3;
1442  }
1443
1444  // Truncate sH region or use form factor.
1445  // Form factor uses either pythia8 renormScale2
1446  // or E_jet in cms.
1447  if (eDcutoff == 1) {
1448  if (sH > pow2(eDLambdaU) ) { sigma *= pow(eDLambdaU,4)/pow2(sH); }
1449  } else if ( (eDgraviton && (eDspin == 2))
1450  && ((eDcutoff == 2) || (eDcutoff == 3)) ) {
1451  double tmPmu = sqrt(Q2RenSave);
1452  if (eDcutoff == 3) tmPmu = (sH + s4 - s3) / (2 * mH);
1453  double tmPformfact = tmPmu / (eDtff * eDLambdaU);
1454  double tmPexp = double(eDnGrav) + 2;
1455  sigma *= 1 / (1 + pow(tmPformfact, tmPexp));
1456  }
1457
1458  return sigma;
1459 }
1460
1461 //--------------------------------------------------------------------------
1462
1463 void Sigma2qg2LEDUnparticleq::setIdColAcol() {
1464
1465  // Flavour set up for q g -> G* q.
1466  int idq = (id2 == 21) ? id1 : id2;
1467  setId( id1, id2, eDidG, idq);
1468
1469  // tH defined between f and f': must swap tHat <-> uHat if q g in.
1470  swapTU = (id2 == 21);
1471
1472  // Colour flow topologies. Swap when antiquarks.
1473  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
1474  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
1475  if (idq < 0) swapColAcol();
1476
1477 }
1478
1479 //==========================================================================
1480
1481 // Sigma2qqbar2LEDUnparticleg class.
1482 // Cross section for q qbar -> U/G g (real graviton emission in
1483 // large extra dimensions or unparticle emission).
1484
1485 //--------------------------------------------------------------------------
1486
1487 void Sigma2qqbar2LEDUnparticleg::initProc() {
1488
1489  // Init model parameters.
1490  eDidG = 5000039;
1491  if (eDgraviton) {
1492  eDspin = (settingsPtr->flag("ExtraDimensionsLED:GravScalar")) ? 0 : 2;
1494  eDdU = 0.5 * eDnGrav + 1;
1496  eDlambda = 1;
1501  } else {
1507  }
1508
1509  // The A(dU) or S'(n) value.
1511  if (eDgraviton) {
1512  tmpAdU = 2 * M_PI * sqrt( pow(M_PI, double(eDnGrav)) )
1513  / GammaReal(0.5 * eDnGrav);
1514  // Scalar graviton
1515  if (eDspin == 0) {
1516  tmpAdU *= 2. * sqrt( pow(2., double(eDnGrav)) );
1517  eDcf *= 4. * eDcf / pow2(eDLambdaU);
1518  double tmpExp = 2. * double(eDnGrav) / (double(eDnGrav) + 2.);
1519  eDgf *= eDgf / pow(2. * M_PI, tmpExp);
1520  }
1521  } else {
1522  tmpAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
1523  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
1524  }
1525
1526  // Cross section related constants
1527  // and ME dependent powers of lambda / LambdaU.
1528  double tmpExp = eDdU - 2;
1529  double tmpLS = pow2(eDLambdaU);
1530  eDconstantTerm = tmpAdU / (2 * 16 * pow2(M_PI) * tmpLS * pow(tmpLS,tmpExp));
1531  if (eDgraviton && (eDspin == 2)) {
1532  eDconstantTerm /= tmpLS;
1533  } else if (eDspin == 1) {
1534  eDconstantTerm *= pow2(eDlambda);
1535  } else if (eDspin == 0) {
1536  eDconstantTerm *= pow2(eDlambda);
1537  } else {
1538  eDconstantTerm = 0;
1539  infoPtr->errorMsg("Error in Sigma2qqbar2LEDUnparticleg::initProc: "
1540  "Incorrect spin value (turn process off)!");
1541  }
1542
1543 }
1544
1545 //--------------------------------------------------------------------------
1546
1547 void Sigma2qqbar2LEDUnparticleg::sigmaKin() {
1548
1549  // Set graviton mass.
1550  mG = m3;
1551  mGS = mG*mG;
1552
1553  // Set mandelstam variables and ME expressions.
1554  if (eDgraviton) {
1555
1556  double A0 = 1/sH;
1557  // Scalar graviton
1558  if (eDspin == 0) {
1559  A0 /= sH;
1560  double tmpTerm1 = uH + tH;
1561  double T0 = (2. * mGS * sH + pow2(tmpTerm1)) / (uH * tH);
1562  double T1 = (tH2 + uH2) / sH;
1563  eDsigma0 = A0 * (eDgf * T0 + eDcf * T1);
1564  } else {
1565  double xH = tH/sH;
1566  double yH = mGS/sH;
1567  double xHS = pow2(xH);
1568  double yHS = pow2(yH);
1569  double xHC = pow(xH,3);
1570  double yHC = pow(yH,3);
1571
1572  double T0 = 1/(xH*(yH-1-xH));
1573  double T1 = -4*xH*(1 + xH)*(1 + 2*xH + 2*xHS);
1574  double T2 = yH*(1 + 6*xH + 18*xHS + 16*xHC);
1575  double T3 = -6*yHS*xH*(1+2*xH);
1576  double T4 = yHC*(1 + 4*xH);
1577
1578  eDsigma0 = A0 * T0 *( T1 + T2 + T3 + T4 );
1579  }
1580
1581  } else if (eDspin == 1) {
1582
1583  double A0 = 1/pow2(sH);
1584  double tmpTerm1 = tH - mGS;
1585  double tmpTerm2 = uH - mGS;
1586
1587  eDsigma0 = A0 * (pow2(tmpTerm1) + pow2(tmpTerm2)) / (tH * uH);
1588
1589  } else if (eDspin == 0) {
1590
1591  double A0 = 1/pow2(sH);
1592
1593  eDsigma0 = A0 * (pow2(sH) - pow2(mGS)) / (tH * uH);
1594
1595  }
1596
1597  // Mass measure, (m^2)^(d-2).
1598  double tmpExp = eDdU - 2;
1599  eDsigma0 *= pow(mGS, tmpExp);
1600
1601  // Constants.
1602  eDsigma0 *= eDconstantTerm;
1603
1604 }
1605
1606 //--------------------------------------------------------------------------
1607
1608 double Sigma2qqbar2LEDUnparticleg::sigmaHat() {
1609
1610  // Mass spactrum weighting.
1611  double sigma = eDsigma0 /runBW3;
1612
1613  // SM couplings...
1614  if (eDgraviton) {
1615  sigma *= 16 * M_PI * alpS / 36;
1616  } else if (eDspin == 1) {
1617  sigma *= 4 * M_PI * 8 * alpS / 9;
1618  } else if (eDspin == 0) {
1619  sigma *= 4 * M_PI * 4 * alpS / 9;
1620  }
1621
1622  // Truncate sH region or use form factor.
1623  // Form factor uses either pythia8 renormScale2
1624  // or E_jet in cms.
1625  if (eDcutoff == 1) {
1626  if (sH > pow2(eDLambdaU) ) { sigma *= pow(eDLambdaU,4)/pow2(sH); }
1627  } else if ( (eDgraviton && (eDspin == 2))
1628  && ((eDcutoff == 2) || (eDcutoff == 3)) ) {
1629  double tmPmu = sqrt(Q2RenSave);
1630  if (eDcutoff == 3) tmPmu = (sH + s4 - s3) / (2 * mH);
1631  double tmPformfact = tmPmu / (eDtff * eDLambdaU);
1632  double tmPexp = double(eDnGrav) + 2;
1633  sigma *= 1 / (1 + pow(tmPformfact, tmPexp));
1634  }
1635
1636  return sigma;
1637 }
1638
1639 //--------------------------------------------------------------------------
1640
1641 void Sigma2qqbar2LEDUnparticleg::setIdColAcol() {
1642
1643  // Flavours trivial.
1644  setId( id1, id2, eDidG, 21);
1645
1646  // Colour flow topologies. Swap when antiquarks.
1647  if (abs(id1) < 9) setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
1648  if (id1 < 0) swapColAcol();
1649
1650 }
1651
1652 //==========================================================================
1653
1654 // Sigma2ffbar2LEDUnparticleZ class.
1655 // Cross section for f fbar -> U/G Z (real LED graviton or unparticle
1656 // emission).
1657
1658 //--------------------------------------------------------------------------
1659
1660 // Constants: could be changed here if desired, but normally should not.
1661 // These are of technical nature, as described for each.
1662
1663 // FIXRATIO:
1664 // Ratio between the two possible coupling constants of the spin-2 ME.
1665 // A value different from one give rise to an IR divergence which makes
1666 // the event generation very slow, so this values is fixed to 1 until
1667 // investigated further.
1668 const double Sigma2ffbar2LEDUnparticleZ::FIXRATIO = 1.;
1669
1670 //--------------------------------------------------------------------------
1671
1672 void Sigma2ffbar2LEDUnparticleZ::initProc() {
1673
1674  // Init model parameters.
1675  eDidG = 5000039;
1676  if (eDgraviton) {
1677  eDspin = 2;
1679  eDdU = 0.5 * eDnGrav + 1;
1681  eDlambda = 1;
1684  } else {
1689  eDratio = FIXRATIO;
1692  }
1693
1694  // Store Z0 mass and width for propagator.
1695  mZ = particleDataPtr->m0(23);
1696  widZ = particleDataPtr->mWidth(23);
1697  mZS = mZ*mZ;
1698  mwZS = pow2(mZ * widZ);
1699
1700  // Init spin-2 parameters
1701  if ( eDspin != 2 ){
1702  eDgraviton = false;
1703  eDlambdaPrime = 0;
1704  } else if (eDgraviton) {
1705  eDlambda = 1;
1706  eDratio = 1;
1707  eDlambdaPrime = eDlambda;
1708  } else {
1709  eDlambdaPrime = eDratio * eDlambda;
1710  }
1711
1712  // The A(dU) or S'(n) value
1713  double tmpAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
1714  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
1715
1716  if (eDgraviton) {
1717  tmpAdU = 2 * M_PI * sqrt( pow(M_PI, double(eDnGrav)) )
1718  / GammaReal(0.5 * eDnGrav);
1719  }
1720
1721  // Standard 2 to 2 cross section related constants
1722  double tmpTerm1 = 1/(2 * 16 * pow2(M_PI));
1723  double tmpLS = pow2(eDLambdaU);
1724
1725  // Spin dependent constants from ME.
1726  double tmpTerm2 = 0;
1727  if ( eDspin == 0 ) {
1728  tmpTerm2 = 2 * pow2(eDlambda);
1729  } else if (eDspin == 1) {
1730  tmpTerm2 = 4 * pow2(eDlambda);
1731  } else if (eDspin == 2) {
1732  tmpTerm2 = pow2(eDlambda)/(4 * 3 * tmpLS);
1733  }
1734
1735  // Unparticle phase space related
1736  double tmpExp2 = eDdU - 2;
1737  double tmpTerm3 = tmpAdU / (tmpLS * pow(tmpLS, tmpExp2));
1738
1739  // All in total
1740  eDconstantTerm = tmpTerm1 * tmpTerm2 * tmpTerm3;
1741
1742  // Secondary open width fraction.
1743  openFrac = particleDataPtr->resOpenFrac(23);
1744
1745 }
1746
1747 //--------------------------------------------------------------------------
1748
1749 void Sigma2ffbar2LEDUnparticleZ::sigmaKin() {
1750
1751  // Set graviton mass and some powers of mandelstam variables
1752  mU = m3;
1753  mUS = mU*mU;
1754
1755  sHS = pow2(sH);
1756  tHS = pow2(tH);
1757  uHS = pow2(uH);
1758  tHC = pow(tH,3);
1759  uHC = pow(uH,3);
1760  tHQ = pow(tH,4);
1761  uHQ = pow(uH,4);
1762  tHuH = tH+uH;
1763
1764  // Evaluate (m**2, t, u) part of differential cross section.
1765  // Extra 1/sHS comes from standard 2 to 2 cross section
1766  // phase space factors.
1767
1768  if ( eDspin == 0 ) {
1769
1770  double A0 = 1/sHS;
1771  double T1 = - sH/tH - sH/uH;
1772  double T2 = - (1 - mZS/tH)*(1 - mUS/tH);
1773  double T3 = - (1 - mZS/uH)*(1 - mUS/uH);
1774  double T4 = 2*(1 - mUS/tH)*(1 - mUS/uH);
1775
1776  eDsigma0 = A0 * ( T1 + T2 + T3 + T4);
1777
1778  } else if ( eDspin == 1 ) {
1779
1780  double A0 = 1/sHS;
1781  double T1 = 0.5 * (tH/uH + uH/tH);
1782  double T2 = pow2(mZS + mUS)/(tH * uH);
1783  double T3 = - 0.5 * mUS * (mZS/tHS + mZS/uHS) ;
1784  double T4 = - (mZS+mUS)*(1/tH + 1/uH);
1785
1786  eDsigma0 = A0 * ( T1 + T2 + T3 + T4 );
1787
1788  } else if ( eDspin == 2 ) {
1789
1790  double A0 = 1 / ( sHS * uHS * tHS * pow2(sH-mZS) );
1791  double F0 = 2*tHS*uHS*( 16*pow(mZS,3) + mUS*(7*tHS + 12*tH*uH + 7*uHS)
1792  - 3*(3*tHC + 11*tHS*uH + 11*tH*uHS + 3*uHC)
1793  + 6*pow(mZS,2)*(7*mUS - 2*tHuH) + mZS*(14*pow(mUS,2)
1794  - 15*tHS - 44*tH*uH - 15*uHS + 2*mUS*tHuH) );
1795  double F2 = 2*tHS*uHS*tHuH*( -8*pow(mZS,2)*tHuH
1796  + 4*mZS*(tHS + 3*tH*uH + uHS)
1797  + 3*(tHC + 5*tHS*uH + 5*tH*uHS + uHC) );
1798  double F4 = -2*tHS*uHS*pow(tHuH,3)*(tHS + uHS - mZS*tHuH);
1799
1800  double G0 = 4*tH*uH*( 6*pow(mZS,3)*(mUS - tH - uH)*tHuH
1801  + pow(mZS,2)*( 9*tHC + 7*tHS*uH + 7*tH*uHS + 9*uHC
1802  + 15*pow2(mUS)*tHuH - 2*mUS*(12*tHS + 19*tH*uH + 12*uHS) )
1803  + tH*uH*( 6*pow(mUS,3) - 9*pow(mUS,2)*tHuH - mUS*(tHS
1804  + 12*tH*uH + uHS) + 6*(tHC + 6*tHS*uH + 6*tH*uHS + uHC) )
1805  + mZS*(-3*tHQ + 25*tHC*uH + 58*tHS*uHS + 25*tH*uHC
1806  - 3*uHQ + 6*pow(mUS,3)*tHuH
1807  - pow(mUS,2)*(15*tHS + 2*tH*uH + 15*uHS) + 2*mUS*(6*tHC
1808  - 11*tHS*uH - 11*tH*uHS + 6*uHC)) );
1809  double G2 = -4*tHS*uHS*tHuH*( -10*pow2(mZS)*tHuH + 2*mZS*(3*tHS
1810  + 7*tH*uH + 3*uHS) + 3*(tHC + 5*tHS*uH + 5*tH*uHS + uHC) );
1811  double G4 = -2*F4;
1812
1813  double H0 = 24*pow(mZS,3)*tH*uH*pow2(-mUS + tHuH)
1814  - 6*pow(mZS,2)*tH*uH*( -9*pow(mUS,3) + 24*pow(mUS,2)*tHuH
1815  - mUS*(21*tHS + 38*tH*uH + 21*uHS)
1816  + 2*(3*tHC + 5*tHS*uH + 5*tH*uHS + 3*uHC) )
1817  - mZS*( 3*pow(mUS,4)*(tHS - 12*tH*uH + uHS)
1818  - 2*tH*uH*pow2(tHuH)*(6*tHS - 29*tH*uH + 6*uHS)
1819  - 6*pow(mUS,3)*(tHC - 16*tHS*uH - 16*tH*uHS + uHC)
1820  + 54*mUS*tH*uH*(tHC + tHS*uH + tH*uHS + uHC)
1821  + pow2(mUS)*(3*tHQ - 102*tHC*uH - 166*tHS*uHS
1822  - 102*tH*uHC + 3*uHQ) )
1823  + tH*uH*( 6*pow(mUS,5) - 18*pow(mUS,4)*tHuH
1824  - 12*pow(mUS,2)*pow(tHuH,3)
1825  + 3*pow(mUS,3)*(7*tHS + 12*tH*uH + 7*uHS)
1826  - 18*tH*uH*(tHC + 5*tHS*uH + 5*tH*uHS + uHC)
1827  + mUS*(3*tHQ + 32*tHC*uH + 78*tHS*uHS + 32*tH*uHC + 3*uHQ) );
1828  double H2 = 2*tHS*uHS*pow2(tHuH)*( -12*pow2(mZS) + 8*mZS*tHuH
1829  + 3*(tHS + 4*tH*uH + uHS) );
1830  double H4 = F4;
1831
1832  eDsigma0 = A0*( F0 + 1/mUS*F2 + 1/pow2(mUS)*F4
1833  + eDratio*(G0 + 1/mUS*G2 + 1/pow2(mUS)*G4)
1834  + pow2(eDratio)*(H0 + 1/mUS*H2 + 1/pow2(mUS)*H4) );
1835
1836  } else {
1837
1838  eDsigma0 = 0;
1839
1840  }
1841
1842 }
1843
1844 //--------------------------------------------------------------------------
1845
1846 double Sigma2ffbar2LEDUnparticleZ::sigmaHat() {
1847
1848  // Electroweak couplings.
1849  int idAbs = abs(id1);
1850  // Note: 1/2 * (g_L^2 + g_R^2) = (g_v^2 + g_a^2)
1851  double facEWS = 4 * M_PI * alpEM
1852  / (couplingsPtr->sin2thetaW() * couplingsPtr->cos2thetaW())
1853  * ( 0.25 * 0.25 * couplingsPtr->vf2af2(idAbs) );
1854
1855  // Mass Spectrum, (m^2)^(d-2)
1856  double tmpExp = eDdU - 2;
1857  double facSpect = pow(mUS, tmpExp);
1858
1859  // Total cross section
1860  double sigma = eDconstantTerm * facEWS * facSpect * eDsigma0;
1861
1862  // Secondary width for Z0
1863  sigma *= openFrac;
1864
1865  // If f fbar are quarks (1/N_c)
1866  if (idAbs < 9) sigma /= 3.;
1867
1868  // Related to mass spactrum weighting.
1869  sigma /= runBW3;
1870
1871  // Truncate sH region or use form factor.
1872  // Form factor uses either pythia8 renormScale2
1873  // or E_jet in cms.
1874  if (eDcutoff == 1) {
1875  if (sH > pow2(eDLambdaU) ) { sigma *= pow(eDLambdaU,4)/pow2(sH); }
1876  } else if (eDgraviton && ((eDcutoff == 2) || (eDcutoff == 3))) {
1877  double tmPmu = sqrt(Q2RenSave);
1878  if (eDcutoff == 3) tmPmu = (sH + s4 - s3) / (2 * mH);
1879  double tmPformfact = tmPmu / (eDtff * eDLambdaU);
1880  double tmPexp = double(eDnGrav) + 2;
1881  sigma *= 1 / (1 + pow(tmPformfact, tmPexp));
1882  }
1883
1884  return sigma;
1885
1886 }
1887
1888 //--------------------------------------------------------------------------
1889
1890 void Sigma2ffbar2LEDUnparticleZ::setIdColAcol() {
1891
1892  // Flavours trivial.
1893  setId( id1, id2, eDidG, 23);
1894
1895  // Colour flow topologies. Swap when antiquarks.
1896  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
1897  else setColAcol( 0, 0, 0, 0, 0, 0);
1898  if (id1 < 0) swapColAcol();
1899
1900 }
1901
1902 //==========================================================================
1903
1904 // Sigma2ffbar2LEDUnparticlegamma class.
1905 // Cross section for f fbar -> U/G gamma (real LED graviton or unparticle
1906 // emission).
1907
1908 //--------------------------------------------------------------------------
1909
1910 // Constants: could be changed here if desired, but normally should not.
1911 // These are of technical nature, as described for each.
1912
1913 // FIXRATIO:
1914 // Ratio between the two possible coupling constants of the spin-2 ME.
1915 // A value different from one give rise to an IR divergence which makes
1916 // the event generation very slow, so this values is fixed to 1 until
1917 // investigated further.
1918 const double Sigma2ffbar2LEDUnparticlegamma::FIXRATIO = 1.;
1919
1920 //--------------------------------------------------------------------------
1921
1922 void Sigma2ffbar2LEDUnparticlegamma::initProc() {
1923
1924  // WARNING: Keep in mind that this class uses the photon limit
1925  // of the Z+G/U ME code. This might give rise to some
1926  // confusing things, e.g. mZ = particleDataPtr->m0(22);
1927
1928  // Init model parameters.
1929  eDidG = 5000039;
1930  if (eDgraviton) {
1931  eDspin = 2;
1933  eDdU = 0.5 * eDnGrav + 1;
1935  eDlambda = 1;
1938  } else {
1943  eDratio = FIXRATIO;
1946  }
1947
1948  // Store Z0 mass.
1949  mZ = particleDataPtr->m0(22);
1950  mZS = mZ*mZ;
1951
1952  // Init spin-2 parameters
1953  if ( eDspin != 2 ){
1954  eDgraviton = false;
1955  eDlambdaPrime = 0;
1956  } else if (eDgraviton) {
1957  eDlambda = 1;
1958  eDratio = 1;
1959  eDlambdaPrime = eDlambda;
1960  } else {
1961  eDlambdaPrime = eDratio * eDlambda;
1962  }
1963
1964  // The A(dU) or S'(n) value
1965  double tmpAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
1966  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
1967
1968  if (eDgraviton) {
1969  tmpAdU = 2 * M_PI * sqrt( pow(M_PI, double(eDnGrav)) )
1970  / GammaReal(0.5 * eDnGrav);
1971  }
1972
1973  // Standard 2 to 2 cross section related constants
1974  double tmpTerm1 = 1/(2 * 16 * pow2(M_PI));
1975  double tmpLS = pow2(eDLambdaU);
1976
1977  // Spin dependent constants from ME.
1978  double tmpTerm2 = 0;
1979  if ( eDspin == 0 ) {
1980  tmpTerm2 = 2 * pow2(eDlambda);
1981  } else if (eDspin == 1) {
1982  tmpTerm2 = 4 * pow2(eDlambda);
1983  } else if (eDspin == 2) {
1984  tmpTerm2 = pow2(eDlambda)/(4 * 3 * tmpLS);
1985  }
1986
1987  // Unparticle phase space related
1988  double tmpExp2 = eDdU - 2;
1989  double tmpTerm3 = tmpAdU / (tmpLS * pow(tmpLS, tmpExp2));
1990
1991  // All in total
1992  eDconstantTerm = tmpTerm1 * tmpTerm2 * tmpTerm3;
1993
1994 }
1995
1996 //--------------------------------------------------------------------------
1997
1998 void Sigma2ffbar2LEDUnparticlegamma::sigmaKin() {
1999
2000  // Set graviton mass and some powers of mandelstam variables
2001  mU = m3;
2002  mUS = mU*mU;
2003
2004  sHS = pow2(sH);
2005  tHS = pow2(tH);
2006  uHS = pow2(uH);
2007  tHC = pow(tH,3);
2008  uHC = pow(uH,3);
2009  tHQ = pow(tH,4);
2010  uHQ = pow(uH,4);
2011  tHuH = tH+uH;
2012
2013  // Evaluate (m**2, t, u) part of differential cross section.
2014  // Extra 1/sHS comes from standard 2 to 2 cross section
2015  // phase space factors.
2016
2017  if ( eDspin == 0 ) {
2018
2019  double A0 = 1/sHS;
2020  double T1 = - sH/tH - sH/uH;
2021  double T2 = - (1 - mZS/tH)*(1 - mUS/tH);
2022  double T3 = - (1 - mZS/uH)*(1 - mUS/uH);
2023  double T4 = 2*(1 - mUS/tH)*(1 - mUS/uH);
2024
2025  eDsigma0 = A0 * ( T1 + T2 + T3 + T4);
2026
2027  } else if ( eDspin == 1 ) {
2028
2029  double A0 = 1/sHS;
2030  double T1 = 0.5 * (tH/uH + uH/tH);
2031  double T2 = pow2(mZS + mUS)/(tH * uH);
2032  double T3 = - 0.5 * mUS * (mZS/tHS + mZS/uHS) ;
2033  double T4 = - (mZS+mUS)*(1/tH + 1/uH);
2034
2035  eDsigma0 = A0 * ( T1 + T2 + T3 + T4 );
2036
2037  } else if ( eDspin == 2 ) {
2038
2039  double A0 = 1 / ( sHS * uHS * tHS * pow2(sH-mZS) );
2040  double F0 = 2*tHS*uHS*( 16*pow(mZS,3) + mUS*(7*tHS + 12*tH*uH + 7*uHS)
2041  - 3*(3*tHC + 11*tHS*uH + 11*tH*uHS + 3*uHC)
2042  + 6*pow(mZS,2)*(7*mUS - 2*tHuH) + mZS*(14*pow(mUS,2)
2043  - 15*tHS - 44*tH*uH - 15*uHS + 2*mUS*tHuH) );
2044  double F2 = 2*tHS*uHS*tHuH*( -8*pow(mZS,2)*tHuH
2045  + 4*mZS*(tHS + 3*tH*uH + uHS)
2046  + 3*(tHC + 5*tHS*uH + 5*tH*uHS + uHC) );
2047  double F4 = -2*tHS*uHS*pow(tHuH,3)*(tHS + uHS - mZS*tHuH);
2048
2049  double G0 = 4*tH*uH*( 6*pow(mZS,3)*(mUS - tH - uH)*tHuH
2050  + pow(mZS,2)*( 9*tHC + 7*tHS*uH + 7*tH*uHS + 9*uHC
2051  + 15*pow2(mUS)*tHuH - 2*mUS*(12*tHS + 19*tH*uH + 12*uHS) )
2052  + tH*uH*( 6*pow(mUS,3) - 9*pow(mUS,2)*tHuH
2053  - mUS*(tHS + 12*tH*uH + uHS)
2054  + 6*(tHC + 6*tHS*uH + 6*tH*uHS + uHC) )
2055  + mZS*(-3*tHQ + 25*tHC*uH + 58*tHS*uHS + 25*tH*uHC
2056  - 3*uHQ + 6*pow(mUS,3)*tHuH
2057  - pow(mUS,2)*(15*tHS + 2*tH*uH + 15*uHS)
2058  + 2*mUS*(6*tHC - 11*tHS*uH - 11*tH*uHS + 6*uHC)) );
2059  double G2 = -4*tHS*uHS*tHuH*( -10*pow2(mZS)*tHuH
2060  + 2*mZS*(3*tHS + 7*tH*uH + 3*uHS)
2061  + 3*(tHC + 5*tHS*uH + 5*tH*uHS + uHC) );
2062  double G4 = -2*F4;
2063
2064  double H0 = 24*pow(mZS,3)*tH*uH*pow2(-mUS + tHuH)
2065  - 6*pow(mZS,2)*tH*uH*( -9*pow(mUS,3) + 24*pow(mUS,2)*tHuH
2066  - mUS*(21*tHS + 38*tH*uH + 21*uHS)
2067  + 2*(3*tHC + 5*tHS*uH + 5*tH*uHS + 3*uHC) )
2068  - mZS*( 3*pow(mUS,4)*(tHS - 12*tH*uH + uHS)
2069  - 2*tH*uH*pow2(tHuH)*(6*tHS - 29*tH*uH + 6*uHS)
2070  - 6*pow(mUS,3)*(tHC - 16*tHS*uH - 16*tH*uHS + uHC)
2071  + 54*mUS*tH*uH*(tHC + tHS*uH + tH*uHS + uHC)
2072  + pow2(mUS)*(3*tHQ - 102*tHC*uH - 166*tHS*uHS
2073  - 102*tH*uHC + 3*uHQ) )
2074  + tH*uH*( 6*pow(mUS,5) - 18*pow(mUS,4)*tHuH
2075  - 12*pow(mUS,2)*pow(tHuH,3)
2076  + 3*pow(mUS,3)*(7*tHS + 12*tH*uH + 7*uHS)
2077  - 18*tH*uH*(tHC + 5*tHS*uH + 5*tH*uHS + uHC)
2078  + mUS*(3*tHQ + 32*tHC*uH + 78*tHS*uHS + 32*tH*uHC + 3*uHQ) );
2079  double H2 = 2*tHS*uHS*pow2(tHuH)*( -12*pow2(mZS) + 8*mZS*tHuH
2080  + 3*(tHS + 4*tH*uH + uHS) );
2081  double H4 = F4;
2082
2083  eDsigma0 = A0*( F0 + 1/mUS*F2 + 1/pow2(mUS)*F4
2084  + eDratio*(G0 + 1/mUS*G2 + 1/pow2(mUS)*G4)
2085  + pow2(eDratio)*(H0 + 1/mUS*H2 + 1/pow2(mUS)*H4) );
2086
2087  } else {
2088
2089  eDsigma0 = 0;
2090
2091  }
2092
2093 }
2094
2095 //--------------------------------------------------------------------------
2096
2097 double Sigma2ffbar2LEDUnparticlegamma::sigmaHat() {
2098
2099  // Electroweak couplings..
2100  int idAbs = abs(id1);
2101  double facEWS = 4 * M_PI * alpEM * couplingsPtr->ef2(idAbs);
2102
2103  // Mass Spectrum, (m^2)^(d-2)
2104  double tmpExp = eDdU - 2;
2105  double facSpect = pow(mUS, tmpExp);
2106
2107  // Total cross section
2108  double sigma = eDconstantTerm * facEWS * facSpect * eDsigma0;
2109
2110  // If f fbar are quarks
2111  if (idAbs < 9) sigma /= 3.;
2112
2113  // Related to mass spactrum weighting.
2114  sigma /= runBW3;
2115
2116  // Truncate sH region or use form factor.
2117  // Form factor uses either pythia8 renormScale2
2118  // or E_jet in cms.
2119  if (eDcutoff == 1) {
2120  if (sH > pow2(eDLambdaU) ) { sigma *= pow(eDLambdaU,4)/pow2(sH); }
2121  } else if (eDgraviton && ((eDcutoff == 2) || (eDcutoff == 3))) {
2122  double tmPmu = sqrt(Q2RenSave);
2123  if (eDcutoff == 3) tmPmu = (sH + s4 - s3) / (2 * mH);
2124  double tmPformfact = tmPmu / (eDtff * eDLambdaU);
2125  double tmPexp = double(eDnGrav) + 2;
2126  sigma *= 1 / (1 + pow(tmPformfact, tmPexp));
2127  }
2128
2129  return sigma;
2130
2131 }
2132
2133 //--------------------------------------------------------------------------
2134
2135 void Sigma2ffbar2LEDUnparticlegamma::setIdColAcol() {
2136
2137  // Flavours trivial.
2138  setId( id1, id2, eDidG, 22);
2139
2140  // Colour flow topologies. Swap when antiquarks.
2141  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2142  else setColAcol( 0, 0, 0, 0, 0, 0);
2143  if (id1 < 0) swapColAcol();
2144
2145 }
2146
2147 //==========================================================================
2148
2149 // Sigma2ffbar2LEDgammagamma class.
2150 // Cross section for f fbar -> (LED G*/U*) -> gamma gamma
2151 // (virtual graviton/unparticle exchange).
2152
2153 //--------------------------------------------------------------------------
2154
2155 void Sigma2ffbar2LEDgammagamma::initProc() {
2156
2157  // Init model parameters.
2158  if (eDgraviton) {
2159  eDspin = 2;
2161  eDdU = 2;
2163  eDlambda = 1;
2167  } else {
2172  eDnegInt = 0;
2173  }
2174
2175  // Model dependent constants.
2176  if (eDgraviton) {
2177  eDlambda2chi = 4*M_PI;
2178  if (eDnegInt == 1) eDlambda2chi *= -1.;
2179  } else {
2180  double tmPAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
2181  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
2182  double tmPdUpi = eDdU * M_PI;
2183  eDlambda2chi = pow2(eDlambda) * tmPAdU / (2 * sin(tmPdUpi));
2184  }
2185
2186  // Model parameter check (if not applicable, sigma = 0).
2187  // Note: SM contribution still generated.
2188  if ( !(eDspin==0 || eDspin==2) ) {
2189  eDlambda2chi = 0;
2190  infoPtr->errorMsg("Error in Sigma2ffbar2LEDgammagamma::initProc: "
2191  "Incorrect spin value (turn process off)!");
2192  } else if ( !eDgraviton && (eDdU >= 2)) {
2193  eDlambda2chi = 0;
2194  infoPtr->errorMsg("Error in Sigma2ffbar2LEDgammagamma::initProc: "
2195  "This process requires dU < 2 (turn process off)!");
2196  }
2197
2198 }
2199
2200 //--------------------------------------------------------------------------
2201
2202 void Sigma2ffbar2LEDgammagamma::sigmaKin() {
2203
2204  // Mandelstam variables.
2205  double sHS = pow2(sH);
2206  double sHQ = pow(sH, 4);
2207  double tHS = pow2(tH);
2208  double uHS = pow2(uH);
2209
2210  // Form factor.
2211  double tmPeffLambdaU = eDLambdaU;
2212  if (eDgraviton && ((eDcutoff == 2) || (eDcutoff == 3))) {
2213  double tmPffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaU);
2214  double tmPexp = double(eDnGrav) + 2;
2215  double tmPformfact = 1 + pow(tmPffterm, tmPexp);
2216  tmPeffLambdaU *= pow(tmPformfact,0.25);
2217  }
2218
2219  // ME from spin-0 and spin-2 unparticles
2220  // including extra 1/sHS from 2-to-2 phase space.
2221  if (eDspin == 0) {
2222  double tmPsLambda2 = sH / pow2(tmPeffLambdaU);
2223  double tmPexp = 2 * eDdU - 1;
2224  eDterm1 = pow(tmPsLambda2,tmPexp);
2225  eDterm1 /= sHS;
2226  } else {
2227  eDterm1 = (uH / tH + tH / uH);
2228  eDterm1 /= sHS;
2229  double tmPsLambda2 = sH / pow2(tmPeffLambdaU);
2230  double tmPexp = eDdU;
2231  eDterm2 = pow(tmPsLambda2,tmPexp) * (uHS + tHS) / sHS;
2232  eDterm2 /= sHS;
2233  tmPexp = 2 * eDdU;
2234  eDterm3 = pow(tmPsLambda2,tmPexp) * tH * uH * (uHS + tHS) / sHQ;
2235  eDterm3 /= sHS;
2236  }
2237
2238 }
2239
2240 //--------------------------------------------------------------------------
2241
2242 double Sigma2ffbar2LEDgammagamma::sigmaHat() {
2243
2244  // Incoming fermion flavor.
2245  int idAbs = abs(id1);
2246
2247  // Couplings and constants.
2248  // Note: ME already contain 1/2 for identical
2249  // particles in the final state.
2250  double sigma = 0;
2251  if (eDspin == 0) {
2252  sigma = pow2(eDlambda2chi) * eDterm1 / 8;
2253  } else {
2254  double tmPe2Q2 = 4 * M_PI * alpEM * couplingsPtr->ef2(idAbs);
2255  double tmPdUpi = eDdU * M_PI;
2256  sigma = pow2(tmPe2Q2) * eDterm1
2257  - tmPe2Q2 * eDlambda2chi * cos(tmPdUpi) * eDterm2
2258  + pow2(eDlambda2chi) * eDterm3 / 4;
2259  }
2260
2261  // dsigma/dt, 2-to-2 phase space factors.
2262  sigma /= 16 * M_PI;
2263
2264  // If f fbar are quarks.
2265  if (idAbs < 9) sigma /= 3.;
2266
2267  return sigma;
2268 }
2269
2270 //--------------------------------------------------------------------------
2271
2272 void Sigma2ffbar2LEDgammagamma::setIdColAcol() {
2273
2274  // Flavours trivial.
2275  setId( id1, id2, 22, 22);
2276
2277  // Colour flow topologies. Swap when antiquarks.
2278  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2279  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2280  if (id1 < 0) swapColAcol();
2281
2282 }
2283
2284 //==========================================================================
2285
2286 // Sigma2gg2LEDgammagamma class.
2287 // Cross section for g g -> (LED G*/U*) -> gamma gamma
2288 // (virtual graviton/unparticle exchange).
2289
2290 //--------------------------------------------------------------------------
2291
2292 void Sigma2gg2LEDgammagamma::initProc() {
2293
2294  // Init model parameters.
2295  if (eDgraviton) {
2296  eDspin = 2;
2298  eDdU = 2;
2300  eDlambda = 1;
2303  } else {
2308  }
2309
2310  // Model dependent constants.
2311  if (eDgraviton) {
2312  eDlambda2chi = 4 * M_PI;
2313
2314  } else {
2315  double tmPAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
2316  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
2317  double tmPdUpi = eDdU * M_PI;
2318  eDlambda2chi = pow2(eDlambda) * tmPAdU / (2 * sin(tmPdUpi));
2319  }
2320
2321  // Model parameter check (if not applicable, sigma = 0).
2322  if ( !(eDspin==0 || eDspin==2) ) {
2323  eDlambda2chi = 0;
2324  infoPtr->errorMsg("Error in Sigma2gg2LEDgammagamma::initProc: "
2325  "Incorrect spin value (turn process off)!");
2326  } else if ( !eDgraviton && (eDdU >= 2)) {
2327  eDlambda2chi = 0;
2328  infoPtr->errorMsg("Error in Sigma2gg2LEDgammagamma::initProc: "
2329  "This process requires dU < 2 (turn process off)!");
2330  }
2331
2332 }
2333
2334 //--------------------------------------------------------------------------
2335
2336 void Sigma2gg2LEDgammagamma::sigmaKin() {
2337
2338  // Mandelstam variables.
2339  double sHS = pow2(sH);
2340  double sHQ = pow(sH, 4);
2341  double tHQ = pow(tH, 4);
2342  double uHQ = pow(uH, 4);
2343
2344  // Form factor.
2345  double tmPeffLambdaU = eDLambdaU;
2346  if (eDgraviton && ((eDcutoff == 2) || (eDcutoff == 3))) {
2347  double tmPffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaU);
2348  double tmPexp = double(eDnGrav) + 2;
2349  double tmPformfact = 1 + pow(tmPffterm, tmPexp);
2350  tmPeffLambdaU *= pow(tmPformfact,0.25);
2351  }
2352
2353  // ME from spin-0 and spin-2 unparticles.
2354  if (eDspin == 0) {
2355  double tmPsLambda2 = sH / pow2(tmPeffLambdaU);
2356  double tmPexp = 2 * eDdU;
2357  eDsigma0 = pow(tmPsLambda2,tmPexp);
2358  } else {
2359  double tmPsLambda2 = sH / pow2(tmPeffLambdaU);
2360  double tmPexp = 2 * eDdU;
2361  eDsigma0 = pow(tmPsLambda2,tmPexp) * (uHQ + tHQ) / sHQ;
2362  }
2363
2364  // extra 1/sHS from 2-to-2 phase space.
2365  eDsigma0 /= sHS;
2366
2367 }
2368
2369 //--------------------------------------------------------------------------
2370
2371 double Sigma2gg2LEDgammagamma::sigmaHat() {
2372
2373  // Couplings and constants.
2374  // Note: ME already contain 1/2 for identical
2375  // particles in the final state.
2376  double sigma = eDsigma0;
2377  if (eDspin == 0) {
2378  sigma *= pow2(eDlambda2chi) / 256;
2379  } else {
2380  sigma *= pow2(eDlambda2chi) / 32;
2381  }
2382
2383  // dsigma/dt, 2-to-2 phase space factors.
2384  sigma /= 16 * M_PI;
2385
2386  return sigma;
2387 }
2388
2389 //--------------------------------------------------------------------------
2390
2391 void Sigma2gg2LEDgammagamma::setIdColAcol() {
2392
2393  // Flavours trivial.
2394  setId( 21, 21, 22, 22);
2395
2396  // Colour flow topologies.
2397  setColAcol( 1, 2, 2, 1, 0, 0, 0, 0);
2398
2399 }
2400
2401 //==========================================================================
2402
2403 // Sigma2ffbar2LEDllbar class.
2404 // Cross section for f fbar -> (LED G*/U*) -> l lbar
2405 // (virtual graviton/unparticle exchange).
2406 // Does not include t-channel contributions relevant for e^+e^- to e^+e^-
2407
2408 //--------------------------------------------------------------------------
2409
2410 void Sigma2ffbar2LEDllbar::initProc() {
2411
2412  // Init model parameters.
2413  if (eDgraviton) {
2414  eDspin = 2;
2416  eDdU = 2;
2418  eDlambda = 1;
2422  } else {
2429  eDnegInt = 0;
2430  }
2431
2432  eDmZ = particleDataPtr->m0(23);
2433  eDmZS = eDmZ * eDmZ;
2434  eDGZ = particleDataPtr->mWidth(23);
2435  eDGZS = eDGZ * eDGZ;
2436
2437  // Model dependent constants.
2438  if (eDgraviton) {
2439  eDlambda2chi = 4*M_PI;
2440  if (eDnegInt == 1) eDlambda2chi *= -1.;
2441  } else {
2442  double tmPAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
2443  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
2444  double tmPdUpi = eDdU * M_PI;
2445  eDlambda2chi = pow2(eDlambda) * tmPAdU / (2 * sin(tmPdUpi));
2446  }
2447
2448  // Model parameter check (if not applicable, sigma = 0).
2449  // Note: SM contribution still generated.
2450  if ( !(eDspin==1 || eDspin==2) ) {
2451  eDlambda2chi = 0;
2452  infoPtr->errorMsg("Error in Sigma2ffbar2LEDllbar::initProc: "
2453  "Incorrect spin value (turn process off)!");
2454  } else if ( !eDgraviton && (eDdU >= 2)) {
2455  eDlambda2chi = 0;
2456  infoPtr->errorMsg("Error in Sigma2ffbar2LEDllbar::initProc: "
2457  "This process requires dU < 2 (turn process off)!");
2458  }
2459
2460 }
2461
2462 //--------------------------------------------------------------------------
2463
2464 void Sigma2ffbar2LEDllbar::sigmaKin() {
2465
2466  // Mandelstam variables.
2467  double tHS = pow2(tH);
2468  double uHS = pow2(uH);
2469  double tHC = pow(tH,3);
2470  double uHC = pow(uH,3);
2471  double tHQ = pow(tH,4);
2472  double uHQ = pow(uH,4);
2473
2474  // Form factor.
2475  double tmPeffLambdaU = eDLambdaU;
2476  if (eDgraviton && ((eDcutoff == 2) || (eDcutoff == 3))) {
2477  double tmPffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaU);
2478  double tmPexp = double(eDnGrav) + 2;
2479  double tmPformfact = 1 + pow(tmPffterm, tmPexp);
2480  tmPeffLambdaU *= pow(tmPformfact,0.25);
2481  }
2482
2483  // ME from spin-1 and spin-2 unparticles
2484  eDdenomPropZ = pow2(sH - eDmZS) + eDmZS * eDGZS;
2485  eDrePropZ = (sH - eDmZS) / eDdenomPropZ;
2486  eDimPropZ = -eDmZ * eDGZ / eDdenomPropZ;
2487  eDrePropGamma = 1 / sH;
2488  if (eDspin == 1) {
2489  double tmPsLambda2 = sH / pow2(tmPeffLambdaU);
2490  double tmPexp = eDdU - 2;
2491  eDabsMeU = eDlambda2chi * pow(tmPsLambda2,tmPexp)
2492  / pow2(tmPeffLambdaU);
2493  } else {
2494  double tmPsLambda2 = sH / pow2(tmPeffLambdaU);
2495  double tmPexp = eDdU - 2;
2496  double tmPA = -eDlambda2chi * pow(tmPsLambda2,tmPexp)
2497  / (8 * pow(tmPeffLambdaU,4));
2498  eDabsAS = pow2(tmPA);
2499  eDreA = tmPA * cos(M_PI * eDdU);
2500  eDreABW = tmPA * ((sH - eDmZS) * cos(M_PI * eDdU) + eDmZ * eDGZ
2501  * sin(M_PI * eDdU)) / eDdenomPropZ;
2502  eDpoly1 = tHQ + uHQ - 6*tHC*uH - 6*tH*uHC + 18*tHS*uHS;
2503  double tmPdiffUT = uH - tH;
2504  eDpoly2 = pow(tmPdiffUT,3);
2505  eDpoly3 = tHC - 3*tHS*uH - 3*tH*uHS + uHC;
2506  }
2507
2508 }
2509
2510 //--------------------------------------------------------------------------
2511
2512 double Sigma2ffbar2LEDllbar::sigmaHat() {
2513
2514  // Incoming fermion flavor.
2515  int idAbs = abs(id1);
2516
2517  // Couplings and constants.
2518  // Qq = couplingsPtr->ef(idAbs), quark, i.e. id > 0.
2519  // Ql = couplingsPtr->ef(11), electron.
2520  double tmPe2QfQl = 4 * M_PI * alpEM * couplingsPtr->ef(idAbs)
2521  * couplingsPtr->ef(11);
2522  double tmPgvq = 0.25 * couplingsPtr->vf(idAbs);
2523  double tmPgaq = 0.25 * couplingsPtr->af(idAbs);
2524  double tmPgLq = tmPgvq + tmPgaq;
2525  double tmPgRq = tmPgvq - tmPgaq;
2526  double tmPgvl = 0.25 * couplingsPtr->vf(11);
2527  double tmPgal = 0.25 * couplingsPtr->af(11);
2528  double tmPgLl = tmPgvl + tmPgal;
2529  double tmPgRl = tmPgvl - tmPgal;
2530  double tmPe2s2c2 = 4 * M_PI * alpEM
2531  / (couplingsPtr->sin2thetaW() * couplingsPtr->cos2thetaW());
2532
2533  // LL, RR, LR, RL couplings.
2534  vector<double> tmPcoupZ;
2535  tmPcoupZ.push_back(tmPe2s2c2 * tmPgLq * tmPgLl);
2536  tmPcoupZ.push_back(tmPe2s2c2 * tmPgRq * tmPgRl);
2537  tmPcoupZ.push_back(tmPe2s2c2 * tmPgRq * tmPgLl);
2538  tmPcoupZ.push_back(tmPe2s2c2 * tmPgLq * tmPgRl);
2539  vector<double> tmPcoupU;
2540  if (eDnxx == 1) {
2541  // LL
2542  tmPcoupU.push_back(-1);
2543  // RR
2544  tmPcoupU.push_back(-1);
2545  } else if (eDnxx == 2) {
2546  // LL
2547  tmPcoupU.push_back(0);
2548  // RR
2549  tmPcoupU.push_back(0);
2550  } else {
2551  // LL
2552  tmPcoupU.push_back(1);
2553  // RR
2554  tmPcoupU.push_back(1);
2555  }
2556  if (eDnxy == 1) {
2557  // RL
2558  tmPcoupU.push_back(-1);
2559  // LR
2560  tmPcoupU.push_back(-1);
2561  } else if (eDnxy == 2) {
2562  // RL
2563  tmPcoupU.push_back(0);
2564  // LR
2565  tmPcoupU.push_back(0);
2566  } else {
2567  // RL
2568  tmPcoupU.push_back(1);
2569  // LR
2570  tmPcoupU.push_back(1);
2571  }
2572
2573  // Matrix elements
2574  double tmPMES = 0;
2575  if (eDspin == 1) {
2576
2577  for (unsigned int i = 0; i<tmPcoupZ.size(); ++i) {
2578  double tmPMS = pow2(tmPcoupU[i] * eDabsMeU)
2579  + pow2(tmPe2QfQl * eDrePropGamma)
2580  + pow2(tmPcoupZ[i]) / eDdenomPropZ
2581  + 2 * cos(M_PI * eDdU) * tmPcoupU[i] * eDabsMeU
2582  * tmPe2QfQl * eDrePropGamma
2583  + 2 * cos(M_PI * eDdU) * tmPcoupU[i] * eDabsMeU
2584  * tmPcoupZ[i] * eDrePropZ
2585  + 2 * tmPe2QfQl * eDrePropGamma
2586  * tmPcoupZ[i] * eDrePropZ
2587  - 2 * sin(M_PI * eDdU) * tmPcoupU[i] * eDabsMeU
2588  * tmPcoupZ[i] * eDimPropZ;
2589
2590  if (i<2) { tmPMES += 4 * pow2(uH) * tmPMS; }
2591  else if (i<4) { tmPMES += 4 * pow2(tH) * tmPMS; }
2592  }
2593
2594  } else {
2595
2596  for (unsigned int i = 0; i<tmPcoupZ.size(); ++i) {
2597  double tmPMS = pow2(tmPe2QfQl * eDrePropGamma)
2598  + pow2(tmPcoupZ[i]) / eDdenomPropZ
2599  + 2 * tmPe2QfQl * eDrePropGamma * tmPcoupZ[i] * eDrePropZ;
2600
2601  if (i<2) { tmPMES += 4 * pow2(uH) * tmPMS; }
2602  else if (i<4) { tmPMES += 4 * pow2(tH) * tmPMS; }
2603  }
2604  tmPMES += 8 * eDabsAS * eDpoly1;
2605  tmPMES += 16 * tmPe2QfQl * eDrePropGamma * eDreA * eDpoly2;
2606  tmPMES += 16 * tmPe2s2c2 * eDreABW * (tmPgaq * tmPgal * eDpoly3
2607  + tmPgvq * tmPgvl * eDpoly2);
2608
2609  }
2610
2611  // dsigma/dt, 2-to-2 phase space factors.
2612  double sigma = 0.25 * tmPMES; // 0.25, is the spin average
2613  sigma /= 16 * M_PI * pow2(sH);
2614
2615  // If f fbar are quarks.
2616  if (idAbs < 9) sigma /= 3.;
2617
2618  // sigma(ffbar->llbar) = 3 * sigma(ffbar->eebar)
2619  sigma *= 3.;
2620
2621  return sigma;
2622 }
2623
2624 //--------------------------------------------------------------------------
2625
2626 void Sigma2ffbar2LEDllbar::setIdColAcol() {
2627
2628  double tmPrand = rndmPtr->flat();
2629  // Flavours trivial.
2630  if (tmPrand < 0.33333333) { setId( id1, id2, 11, -11); }
2631  else if (tmPrand < 0.66666667) { setId( id1, id2, 13, -13); }
2632  else { setId( id1, id2, 15, -15); }
2633
2634  // tH defined between f and f': must swap tHat <-> uHat if id1 is fbar.
2635  swapTU = (id2 > 0);
2636
2637  // Colour flow topologies. Swap when antiquarks.
2638  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2639  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2640  if (id1 < 0) swapColAcol();
2641
2642 }
2643
2644 //==========================================================================
2645
2646 // Sigma2gg2LEDllbar class.
2647 // Cross section for g g -> (LED G*/U*) -> l lbar
2648 // (virtual graviton/unparticle exchange).
2649
2650 //--------------------------------------------------------------------------
2651
2652 void Sigma2gg2LEDllbar::initProc() {
2653
2654  // Init model parameters.
2655  if (eDgraviton) {
2656  eDspin = 2;
2658  eDdU = 2;
2660  eDlambda = 1;
2663  } else {
2668  }
2669
2670  // Model dependent constants.
2671  if (eDgraviton) {
2672  eDlambda2chi = 4 * M_PI;
2673
2674  } else {
2675  double tmPAdU = 16 * pow2(M_PI) * sqrt(M_PI) / pow(2. * M_PI, 2. * eDdU)
2676  * GammaReal(eDdU + 0.5) / (GammaReal(eDdU - 1.) * GammaReal(2. * eDdU));
2677  double tmPdUpi = eDdU * M_PI;
2678  eDlambda2chi = pow2(eDlambda) * tmPAdU / (2 * sin(tmPdUpi));
2679  }
2680
2681  // Model parameter check (if not applicable, sigma = 0).
2682  if ( !(eDspin==2) ) {
2683  eDlambda2chi = 0;
2684  infoPtr->errorMsg("Error in Sigma2gg2LEDllbar::initProc: "
2685  "Incorrect spin value (turn process off)!");
2686  } else if ( !eDgraviton && (eDdU >= 2)) {
2687  eDlambda2chi = 0;
2688  infoPtr->errorMsg("Error in Sigma2gg2LEDllbar::initProc: "
2689  "This process requires dU < 2 (turn process off)!");
2690  }
2691
2692 }
2693
2694 //--------------------------------------------------------------------------
2695
2696 void Sigma2gg2LEDllbar::sigmaKin() {
2697
2698  // Form factor.
2699  double tmPeffLambdaU = eDLambdaU;
2700  if (eDgraviton && ((eDcutoff == 2) || (eDcutoff == 3))) {
2701  double tmPffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaU);
2702  double tmPexp = double(eDnGrav) + 2;
2703  double tmPformfact = 1 + pow(tmPffterm, tmPexp);
2704  tmPeffLambdaU *= pow(tmPformfact,0.25);
2705  }
2706
2707  // ME from spin-2 unparticle.
2708  double tmPsLambda2 = sH / pow2(tmPeffLambdaU);
2709  double tmPexp = eDdU - 2;
2710  double tmPA = -eDlambda2chi * pow(tmPsLambda2,tmPexp)
2711  / (8 * pow(tmPeffLambdaU,4));
2712  eDsigma0 = 4 * pow2(tmPA) * uH * tH * (pow2(uH) + pow2(tH));
2713
2714  // extra 1/sHS from 2-to-2 phase space.
2715  eDsigma0 /= 16 * M_PI * pow2(sH);
2716
2717  // sigma(ffbar->llbar) = 3 * sigma(ffbar->eebar)
2718  eDsigma0 *= 3.;
2719
2720 }
2721
2722 //--------------------------------------------------------------------------
2723
2724 void Sigma2gg2LEDllbar::setIdColAcol() {
2725
2726  double tmPrand = rndmPtr->flat();
2727  // Flavours trivial.
2728  if (tmPrand < 0.33333333) { setId( 21, 21, 11, -11); }
2729  else if (tmPrand < 0.66666667) { setId( 21, 21, 13, -13); }
2730  else { setId( 21, 21, 15, -15); }
2731
2732  // Colour flow topologies.
2733  setColAcol( 1, 2, 2, 1, 0, 0, 0, 0);
2734
2735 }
2736
2737 //==========================================================================
2738
2739 // Sigma2gg2LEDgg class.
2740 // Cross section for g g -> (LED G*) -> g g.
2741
2742 //--------------------------------------------------------------------------
2743
2744 // Initialize process.
2745
2746 void Sigma2gg2LEDgg::initProc() {
2747
2756
2757 }
2758
2759 //--------------------------------------------------------------------------
2760
2761 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
2762
2763 void Sigma2gg2LEDgg::sigmaKin() {
2764
2765  // Get S(x) values for G amplitude.
2766  complex sS(0., 0.);
2767  complex sT(0., 0.);
2768  complex sU(0., 0.);
2769  if (eDopMode == 0) {
2770  sS = ampLedS( sH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2771  sT = ampLedS( tH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2772  sU = ampLedS( uH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2773  } else {
2774  // Form factor.
2775  double effLambda = eDLambdaT;
2776  if ((eDcutoff == 2) || (eDcutoff == 3)) {
2777  double ffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaT);
2778  double exp = double(eDnGrav) + 2.;
2779  double formfa = 1. + pow(ffterm, exp);
2780  effLambda *= pow(formfa,0.25);
2781  }
2782  sS = 4.*M_PI/pow(effLambda,4);
2783  sT = 4.*M_PI/pow(effLambda,4);
2784  sU = 4.*M_PI/pow(effLambda,4);
2785  if (eDnegInt == 1) {
2786  sS *= -1.;
2787  sT *= -1.;
2788  sU *= -1.;
2789  }
2790  }
2791
2792  // Calculate kinematics dependence.
2793  double sH3 = sH*sH2;
2794  double tH3 = tH*tH2;
2795  double uH3 = uH*uH2;
2796
2797  sigTS = (128. * pow2(M_PI) * pow2(alpS)) * (9./4.)
2798  * (tH2 / sH2 + 2. * tH / sH + 3. + 2. * sH / tH + sH2 / tH2)
2799  + 24.*M_PI*alpS*( (sH3/tH + tH2 + 3.*(sH*tH + sH2))*sS.real()
2800  + (tH3/sH + sH2 + 3.*(tH*sH + tH2))*sT.real())
2801  + pow2(uH2)*( 4.*real(sS*conj(sS)) + sS.real()*sT.real()
2802  + sS.imag()*sT.imag() + 4.*real(sT*conj(sT)));
2803
2804
2805  sigUS = (128. * pow2(M_PI) * pow2(alpS)) * (9./4.)
2806  * (uH2 / sH2 + 2. * uH / sH + 3. + 2. * sH / uH + sH2 / uH2)
2807  + 24.*M_PI*alpS*( (sH3/uH + uH2 + 3.*(sH*uH + sH2))*sS.real()
2808  + (uH3/sH + sH2 + 3.*(uH*sH + uH2))*sU.real())
2809  + pow2(tH2)*( 4.*real(sS*conj(sS)) + sS.real()*sU.real()
2810  + sS.imag()*sU.imag() + 4.*real(sU*conj(sU)));
2811
2812  sigTU = (128. * pow2(M_PI) * pow2(alpS)) * (9./4.)
2813  * (tH2 / uH2 + 2. * tH / uH + 3. + 2. * uH / tH + uH2 / tH2)
2814  + 24.*M_PI*alpS*( (tH3/uH + uH2 + 3.*(tH*uH + tH2))*sT.real()
2815  + (uH3/tH + tH2 + 3.*(uH*tH + uH2))*sU.real())
2816  + pow2(sH2)*( 4.*real(sT*conj(sT)) + sT.real()*sU.real()
2817  + sT.imag()*sU.imag() + 4.*real(sU*conj(sU)));
2818
2819  sigSum = sigTS + sigUS + sigTU;
2820
2821  // Answer contains factor 1/2 from identical gluons.
2822  sigma = 0.5 * sigSum / (128. * M_PI * sH2);
2823
2824 }
2825
2826 //--------------------------------------------------------------------------
2827
2828 // Select identity, colour and anticolour.
2829
2830 void Sigma2gg2LEDgg::setIdColAcol() {
2831
2832  // Flavours are trivial.
2833  setId( id1, id2, 21, 21);
2834
2835  // Three colour flow topologies, each with two orientations.
2836  double sigRand = sigSum * rndmPtr->flat();
2837  if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
2838  else if (sigRand < sigTS + sigUS)
2839  setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
2840  else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
2841  if (rndmPtr->flat() > 0.5) swapColAcol();
2842
2843 }
2844
2845 //==========================================================================
2846
2847 // Sigma2gg2LEDqqbar class.
2848 // Cross section for g g -> (LED G*) -> q qbar.
2849
2850 //--------------------------------------------------------------------------
2851
2852 // Initialize process.
2853
2854 void Sigma2gg2LEDqqbar::initProc() {
2855
2856  // Read number of quarks to be considered in massless approximation
2857  // as well as model parameters.
2866
2867 }
2868
2869 //--------------------------------------------------------------------------
2870
2871 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
2872
2873 void Sigma2gg2LEDqqbar::sigmaKin() {
2874
2875  // Get S(x) values for G amplitude.
2876  complex sS(0., 0.);
2877  complex sT(0., 0.);
2878  complex sU(0., 0.);
2879  if (eDopMode == 0) {
2880  sS = ampLedS( sH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2881  sT = ampLedS( tH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2882  sU = ampLedS( uH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2883  } else {
2884  // Form factor.
2885  double effLambda = eDLambdaT;
2886  if ((eDcutoff == 2) || (eDcutoff == 3)) {
2887  double ffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaT);
2888  double exp = double(eDnGrav) + 2.;
2889  double formfa = 1. + pow(ffterm, exp);
2890  effLambda *= pow(formfa,0.25);
2891  }
2892  sS = 4.*M_PI/pow(effLambda,4);
2893  sT = 4.*M_PI/pow(effLambda,4);
2894  sU = 4.*M_PI/pow(effLambda,4);
2895  if (eDnegInt == 1) {
2896  sS *= -1.;
2897  sT *= -1.;
2898  sU *= -1.;
2899  }
2900  }
2901
2902  // Pick new flavour.
2903  idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
2904  mNew = particleDataPtr->m0(idNew);
2905  m2New = mNew*mNew;
2906
2907  // Calculate kinematics dependence.
2908  sigTS = 0.;
2909  sigUS = 0.;
2910  if (sH > 4. * m2New) {
2911  double tH3 = tH*tH2;
2912  double uH3 = uH*uH2;
2913  sigTS = (16. * pow2(M_PI) * pow2(alpS))
2914  * ((1./6.) * uH / tH - (3./8.) * uH2 / sH2)
2915  - 0.5 * M_PI * alpS * uH2 * sS.real()
2916  + (3./16.) * uH3 * tH * real(sS*conj(sS));
2917  sigUS = (16. * pow2(M_PI) * pow2(alpS))
2918  * ((1./6.) * tH / uH - (3./8.) * tH2 / sH2)
2919  - 0.5 * M_PI * alpS * tH2 * sS.real()
2920  + (3./16.) * tH3 * uH * real(sS*conj(sS));
2921  }
2922  sigSum = sigTS + sigUS;
2923
2924  // Answer is proportional to number of outgoing flavours.
2925  sigma = nQuarkNew * sigSum / (16. * M_PI * sH2);
2926
2927 }
2928
2929 //--------------------------------------------------------------------------
2930
2931 // Select identity, colour and anticolour.
2932
2933 void Sigma2gg2LEDqqbar::setIdColAcol() {
2934
2935  // Flavours are trivial.
2936  setId( id1, id2, idNew, -idNew);
2937
2938  // Two colour flow topologies.
2939  double sigRand = sigSum * rndmPtr->flat();
2940  if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
2941  else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
2942
2943 }
2944
2945 //==========================================================================
2946
2947 // Sigma2qg2LEDqg class.
2948 // Cross section for q g -> (LED G*) -> q g.
2949
2950 //--------------------------------------------------------------------------
2951
2952 // Initialize process.
2953
2954 void Sigma2qg2LEDqg::initProc() {
2955
2964
2965 }
2966
2967 //--------------------------------------------------------------------------
2968
2969 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
2970
2971 void Sigma2qg2LEDqg::sigmaKin() {
2972
2973  // Get S(x) values for G amplitude.
2974  complex sS(0., 0.);
2975  complex sT(0., 0.);
2976  complex sU(0., 0.);
2977  if (eDopMode == 0) {
2978  sS = ampLedS( sH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2979  sT = ampLedS( tH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2980  sU = ampLedS( uH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
2981  } else {
2982  // Form factor.
2983  double effLambda = eDLambdaT;
2984  if ((eDcutoff == 2) || (eDcutoff == 3)) {
2985  double ffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaT);
2986  double exp = double(eDnGrav) + 2.;
2987  double formfa = 1. + pow(ffterm, exp);
2988  effLambda *= pow(formfa,0.25);
2989  }
2990  sS = 4.*M_PI/pow(effLambda,4);
2991  sT = 4.*M_PI/pow(effLambda,4);
2992  sU = 4.*M_PI/pow(effLambda,4);
2993  if (eDnegInt == 1) {
2994  sS *= -1.;
2995  sT *= -1.;
2996  sU *= -1.;
2997  }
2998  }
2999
3000  // Calculate kinematics dependence.
3001  double sH3 = sH*sH2;
3002  double uH3 = uH*uH2;
3003  sigTS = (16. * pow2(M_PI) * pow2(alpS))
3004  * (uH2 / tH2 - (4./9.) * uH / sH)
3005  + (4./3.) * M_PI * alpS * uH2 * sT.real()
3006  - 0.5 * uH3 * sH * real(sT*conj(sT));
3007  sigTU = (16. * pow2(M_PI) * pow2(alpS))
3008  * (sH2 / tH2 - (4./9.) * sH / uH)
3009  + (4./3.) * M_PI * alpS * sH2 * sT.real()
3010  - 0.5 * sH3 * uH * real(sT*conj(sT));
3011  sigSum = sigTS + sigTU;
3012
3014  sigma = sigSum / (16. * M_PI * sH2);
3015
3016 }
3017
3018 //--------------------------------------------------------------------------
3019
3020 // Select identity, colour and anticolour.
3021
3022 void Sigma2qg2LEDqg::setIdColAcol() {
3023
3024  // Outgoing = incoming flavours.
3025  setId( id1, id2, id1, id2);
3026
3027  // Two colour flow topologies. Swap if first is gluon, or when antiquark.
3028  double sigRand = sigSum * rndmPtr->flat();
3029  if (sigRand < sigTS) setColAcol( 1, 0, 2, 1, 3, 0, 2, 3);
3030  else setColAcol( 1, 0, 2, 3, 2, 0, 1, 3);
3031  if (id1 == 21) swapCol1234();
3032  if (id1 < 0 || id2 < 0) swapColAcol();
3033
3034 }
3035
3036 //==========================================================================
3037
3038 // Sigma2qq2LEDqq class.
3039 // Cross section for q q(bar)' -> (LED G*) -> q q(bar)'
3040
3041 //--------------------------------------------------------------------------
3042
3043 // Initialize process.
3044
3045 void Sigma2qq2LEDqq::initProc() {
3046
3055
3056 }
3057
3058 //--------------------------------------------------------------------------
3059
3060 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
3061
3062 void Sigma2qq2LEDqq::sigmaKin() {
3063
3064  // Get S(x) values for G amplitude.
3065  complex sS(0., 0.);
3066  complex sT(0., 0.);
3067  complex sU(0., 0.);
3068  if (eDopMode == 0) {
3069  sS = ampLedS( sH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3070  sT = ampLedS( tH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3071  sU = ampLedS( uH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3072  } else {
3073  // Form factor.
3074  double effLambda = eDLambdaT;
3075  if ((eDcutoff == 2) || (eDcutoff == 3)) {
3076  double ffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaT);
3077  double exp = double(eDnGrav) + 2.;
3078  double formfa = 1. + pow(ffterm, exp);
3079  effLambda *= pow(formfa,0.25);
3080  }
3081  sS = 4.*M_PI/pow(effLambda,4);
3082  sT = 4.*M_PI/pow(effLambda,4);
3083  sU = 4.*M_PI/pow(effLambda,4);
3084  if (eDnegInt == 1) {
3085  sS *= -1.;
3086  sT *= -1.;
3087  sU *= -1.;
3088  }
3089  }
3090
3091  // Calculate kinematics dependence for different terms.
3092  sigT = (4./9.) * (sH2 + uH2) / tH2;
3093  sigU = (4./9.) * (sH2 + tH2) / uH2;
3094  sigTU = - (8./27.) * sH2 / (tH * uH);
3095  sigST = - (8./27.) * uH2 / (sH * tH);
3096  // Graviton terms.
3097  sigGrT1 = funLedG(tH, uH) * real(sT*conj(sT)) / 8.;
3098  sigGrT2 = funLedG(tH, sH) * real(sT*conj(sT)) / 8.;
3099  sigGrU = funLedG(uH, tH) * real(sU*conj(sU)) / 8.;
3100  sigGrTU = (8./9.) * M_PI * alpS * sH2
3101  * ((4.*uH + tH)*sT.real()/uH + (4.*tH + uH)*sU.real()/tH)
3102  + (sT.real()*sU.real() + sT.imag()*sU.imag())
3103  * (4.*tH + uH)*(4.*uH + tH) * sH2 / 48.;
3104  sigGrST = (8./9.) * M_PI * alpS * uH2
3105  * ((4.*tH + sH)*sS.real()/tH + (4.*sH + tH)*sT.real()/sH)
3106  + (sS.real()*sT.real() + sS.imag()*sT.imag())
3107  * (4.*sH + tH)*(4.*tH + sH) * uH2 / 48.;
3108
3109 }
3110
3111 //--------------------------------------------------------------------------
3112
3113 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
3114
3115 double Sigma2qq2LEDqq::sigmaHat() {
3116
3117  // Combine cross section terms; factor 1/2 when identical quarks.
3118  if (id2 == id1) {
3119  sigSum = (16. * pow2(M_PI) * pow2(alpS)) * (sigT + sigU + sigTU)
3120  + sigGrT1 + sigGrU + sigGrTU;
3121  sigSum *= 0.5;
3122  } else if (id2 == -id1) {
3123  sigSum = (16. * pow2(M_PI) * pow2(alpS)) * (sigT + sigST)
3124  + sigGrT2 + sigGrST;
3125  } else {
3126  sigSum = 16. * pow2(M_PI) * pow2(alpS) * sigT + sigGrT1;
3127  }
3128
3130  return sigSum / (16. * M_PI * sH2);
3131
3132 }
3133
3134 //--------------------------------------------------------------------------
3135
3136 // Select identity, colour and anticolour.
3137
3138 void Sigma2qq2LEDqq::setIdColAcol() {
3139
3140  // Outgoing = incoming flavours.
3141  setId( id1, id2, id1, id2);
3142
3143  // Colour flow topologies. Swap when antiquarks.
3144  double sigTtot = sigT + sigGrT2;
3145  double sigUtot = sigU + sigGrU;
3146  if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
3147  else setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
3148  if (id2 == id1 && (sigTtot + sigUtot) * rndmPtr->flat() > sigTtot)
3149  setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
3150  if (id1 < 0) swapColAcol();
3151
3152 }
3153
3154 //==========================================================================
3155
3156 // Sigma2qqbar2LEDgg class.
3157 // Cross section for q qbar -> (LED G*) -> g g.
3158
3159 //--------------------------------------------------------------------------
3160
3161 // Initialize process.
3162
3163 void Sigma2qqbar2LEDgg::initProc() {
3164
3173
3174 }
3175
3176 //--------------------------------------------------------------------------
3177
3178 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
3179
3180 void Sigma2qqbar2LEDgg::sigmaKin() {
3181
3182  // Get S(x) values for G amplitude.
3183  complex sS(0., 0.);
3184  complex sT(0., 0.);
3185  complex sU(0., 0.);
3186  if (eDopMode == 0) {
3187  sS = ampLedS( sH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3188  sT = ampLedS( tH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3189  sU = ampLedS( uH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3190  } else {
3191  // Form factor.
3192  double effLambda = eDLambdaT;
3193  if ((eDcutoff == 2) || (eDcutoff == 3)) {
3194  double ffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaT);
3195  double exp = double(eDnGrav) + 2.;
3196  double formfa = 1. + pow(ffterm, exp);
3197  effLambda *= pow(formfa,0.25);
3198  }
3199  sS = 4.*M_PI/pow(effLambda,4);
3200  sT = 4.*M_PI/pow(effLambda,4);
3201  sU = 4.*M_PI/pow(effLambda,4);
3202  if (eDnegInt == 1) {
3203  sS *= -1.;
3204  sT *= -1.;
3205  sU *= -1.;
3206  }
3207  }
3208
3209  // Calculate kinematics dependence.
3210  double tH3 = tH*tH2;
3211  double uH3 = uH*uH2;
3212  sigTS = (16. * pow2(M_PI) * pow2(alpS))
3213  * ((1./6.) * uH / tH - (3./8.) * uH2 / sH2)
3214  - 0.5 * M_PI * alpS * uH2 * sS.real()
3215  + (3./16.) * uH3 * tH * real(sS*conj(sS));
3216  sigUS = (16. * pow2(M_PI) * pow2(alpS))
3217  * ((1./6.) * tH / uH - (3./8.) * tH2 / sH2)
3218  - 0.5 * M_PI * alpS * tH2 * sS.real()
3219  + (3./16.) * tH3 * uH * real(sS*conj(sS));
3220
3221  sigSum = sigTS + sigUS;
3222
3223  // Answer contains factor 1/2 from identical gluons.
3224  sigma = (64./9.) * 0.5 * sigSum / (16. * M_PI * sH2);
3225
3226 }
3227
3228 //--------------------------------------------------------------------------
3229
3230 // Select identity, colour and anticolour.
3231
3232 void Sigma2qqbar2LEDgg::setIdColAcol() {
3233
3234  // Outgoing flavours trivial.
3235  setId( id1, id2, 21, 21);
3236
3237  // Two colour flow topologies. Swap if first is antiquark.
3238  double sigRand = sigSum * rndmPtr->flat();
3239  if (sigRand < sigTS) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
3240  else setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
3241  if (id1 < 0) swapColAcol();
3242
3243 }
3244
3245 //==========================================================================
3246
3247 // Sigma2qqbar2LEDqqbarNew class.
3248 // Cross section q qbar -> (LED G*) -> q' qbar'.
3249
3250 //--------------------------------------------------------------------------
3251
3252 // Initialize process.
3253
3254 void Sigma2qqbar2LEDqqbarNew::initProc() {
3255
3256  // Read number of quarks to be considered in massless approximation
3257  // as well as model parameters.
3265
3266 }
3267
3268 //--------------------------------------------------------------------------
3269
3270 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
3271
3272 void Sigma2qqbar2LEDqqbarNew::sigmaKin() {
3273
3274  // Get S(x) values for G amplitude.
3275  complex sS(0., 0.);
3276  complex sT(0., 0.);
3277  complex sU(0., 0.);
3278  if (eDopMode == 0) {
3279  sS = ampLedS( sH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3280  sT = ampLedS( tH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3281  sU = ampLedS( uH/pow2(eDLambdaT), eDnGrav, eDLambdaT, eDMD);
3282  } else {
3283  // Form factor.
3284  double effLambda = eDLambdaT;
3285  if ((eDcutoff == 2) || (eDcutoff == 3)) {
3286  double ffterm = sqrt(Q2RenSave) / (eDtff * eDLambdaT);
3287  double exp = double(eDnGrav) + 2.;
3288  double formfa = 1. + pow(ffterm, exp);
3289  effLambda *= pow(formfa,0.25);
3290  }
3291  sS = 4.*M_PI/pow(effLambda,4);
3292  sT = 4.*M_PI/pow(effLambda,4);
3293  sU = 4.*M_PI/pow(effLambda,4);
3294  }
3295
3296  // Pick new flavour.
3297  idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
3298  mNew = particleDataPtr->m0(idNew);
3299  m2New = mNew*mNew;
3300
3301  // Calculate kinematics dependence.
3302  sigS = 0.;
3303  if (sH > 4. * m2New) {
3304  sigS = (16. * pow2(M_PI) * pow2(alpS))
3305  * (4./9.) * (tH2 + uH2) / sH2
3306  + funLedG(sH, tH) * real(sS*conj(sS)) / 8.;
3307  }
3308  // Answer is proportional to number of outgoing flavours.
3309  sigma = nQuarkNew * sigS / (16. * M_PI * sH2);
3310
3311 }
3312
3313 //--------------------------------------------------------------------------
3314
3315 // Select identity, colour and anticolour.
3316
3317 void Sigma2qqbar2LEDqqbarNew::setIdColAcol() {
3318
3319  // Set outgoing flavours ones.
3320  id3 = (id1 > 0) ? idNew : -idNew;
3321  setId( id1, id2, id3, -id3);
3322
3323  // Colour flow topologies. Swap when antiquarks.
3324  setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
3325  if (id1 < 0) swapColAcol();
3326
3327 }
3328
3329 //==========================================================================
3330
3331 } // end namespace Pythia8
Definition: AgUStep.h:26