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