source: HiSusy/trunk/Pythia8/pythia8170/src/SigmaExtraDim.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 105.6 KB
Line 
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
12namespace 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
19complex 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
79double 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 
94void 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
131void 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
153void 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 
169double 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 
243void 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
280void 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
298double 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
316void 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
334double 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 
404void 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
438void 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
489double 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
504void 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 
519double 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 
576void 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
597void 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
612void 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 
627double 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 
651void 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
672void 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
687void 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 
707double 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 
731void 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
752void 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
767void 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 
782double 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
807void 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
904int Sigma2ffbar2TEVffbar::resonanceB() {
905
906  return 23;
907
908}
909
910//--------------------------------------------------------------------------
911
912// For improving the phase-space sampling (there can be 2 resonances)
913
914int 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
932void 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
955double 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
1107void 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
1126double 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
1144void 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
1198void 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
1257double 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
1288void 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
1307void 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
1368void 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
1431double 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
1464void 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
1488void 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
1548void 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
1609double 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
1642void 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.
1669const double Sigma2ffbar2LEDUnparticleZ::FIXRATIO = 1.;
1670
1671//--------------------------------------------------------------------------
1672
1673void 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
1747void 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
1844double 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
1885void 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.
1913const double Sigma2ffbar2LEDUnparticlegamma::FIXRATIO = 1.;
1914
1915//--------------------------------------------------------------------------
1916
1917void 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
1993void 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
2092double 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
2130void 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
2150void 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
2197void 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
2237double 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
2267void 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
2287void 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
2331void 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
2366double 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
2386void 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
2405void 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
2459void 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
2507double 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
2621void 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
2647void 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
2691void 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
2719void 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 
2741void 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
2758void 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
2825void 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 
2849void 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
2868void 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
2928void 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 
2949void 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
2966void 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
3017void 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 
3040void 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
3057void 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
3110double 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
3133void 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 
3158void 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
3175void 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
3227void 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 
3249void 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
3267void 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
3312void 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
Note: See TracBrowser for help on using the repository browser.