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