PYTHIA8/pythia8170/src/SigmaExtraDim.cxx

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