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