| 5ad4eb21 |
1 | // SigmaProcess.cc is a part of the PYTHIA event generator. |
| 2 | // Copyright (C) 2008 Torbjorn Sjostrand. |
| 3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
| 4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
| 5 | |
| 6 | // Function definitions (not found in the header) for the |
| 7 | // SigmaProcess class, and classes derived from it. |
| 8 | |
| 9 | #include "SigmaProcess.h" |
| 10 | |
| 11 | namespace Pythia8 { |
| 12 | |
| 13 | //************************************************************************** |
| 14 | |
| 15 | // The SigmaProcess class. |
| 16 | // Base class for cross sections. |
| 17 | |
| 18 | //********* |
| 19 | |
| 20 | // Constants: could be changed here if desired, but normally should not. |
| 21 | // These are of technical nature, as described for each. |
| 22 | |
| 23 | // Conversion of GeV^{-2} to mb for cross section. |
| 24 | const double SigmaProcess::CONVERT2MB = 0.389380; |
| 25 | |
| 26 | // The sum of outgoing masses must not be too close to the cm energy. |
| 27 | const double SigmaProcess::MASSMARGIN = 0.1; |
| 28 | |
| 29 | //********* |
| 30 | |
| 31 | // Perform simple initialization and store pointers. |
| 32 | |
| 33 | void SigmaProcess::init(Info* infoPtrIn, BeamParticle* beamAPtrIn, |
| 34 | BeamParticle* beamBPtrIn, AlphaStrong* alphaSPtrIn, |
| 35 | AlphaEM* alphaEMPtrIn,SigmaTotal* sigmaTotPtrIn, |
| 36 | SusyLesHouches* slhaPtrIn) { |
| 37 | |
| 38 | // Store pointers. |
| 39 | infoPtr = infoPtrIn; |
| 40 | beamAPtr = beamAPtrIn; |
| 41 | beamBPtr = beamBPtrIn; |
| 42 | alphaSPtr = alphaSPtrIn; |
| 43 | alphaEMPtr = alphaEMPtrIn; |
| 44 | sigmaTotPtr = sigmaTotPtrIn; |
| 45 | slhaPtr = slhaPtrIn; |
| 46 | |
| 47 | // Read out some properties of beams to allow shorthand. |
| 48 | idA = beamAPtr->id(); |
| 49 | idB = beamBPtr->id(); |
| 50 | mA = beamAPtr->m(); |
| 51 | mB = beamBPtr->m(); |
| 52 | isLeptonA = beamAPtr->isLepton(); |
| 53 | isLeptonB = beamBPtr->isLepton(); |
| 54 | hasLeptonBeams = isLeptonA || isLeptonB; |
| 55 | |
| 56 | // K factor, multiplying resolved processes. (But not here for MI.) |
| 57 | Kfactor = Settings::parm("SigmaProcess:Kfactor"); |
| 58 | |
| 59 | // Maximum incoming quark flavour. |
| 60 | nQuarkIn = Settings::mode("PDFinProcess:nQuarkIn"); |
| 61 | |
| 62 | // Renormalization scale choice. |
| 63 | renormScale1 = Settings::mode("SigmaProcess:renormScale1"); |
| 64 | renormScale2 = Settings::mode("SigmaProcess:renormScale2"); |
| 65 | renormScale3 = Settings::mode("SigmaProcess:renormScale3"); |
| 66 | renormScale3VV = Settings::mode("SigmaProcess:renormScale3VV"); |
| 67 | renormMultFac = Settings::parm("SigmaProcess:renormMultFac"); |
| 68 | renormFixScale = Settings::parm("SigmaProcess:renormFixScale"); |
| 69 | |
| 70 | // Factorization scale choice. |
| 71 | factorScale1 = Settings::mode("SigmaProcess:factorScale1"); |
| 72 | factorScale2 = Settings::mode("SigmaProcess:factorScale2"); |
| 73 | factorScale3 = Settings::mode("SigmaProcess:factorScale3"); |
| 74 | factorScale3VV = Settings::mode("SigmaProcess:factorScale3VV"); |
| 75 | factorMultFac = Settings::parm("SigmaProcess:factorMultFac"); |
| 76 | factorFixScale = Settings::parm("SigmaProcess:factorFixScale"); |
| 77 | |
| 78 | // CP violation parameters for the BSM Higgs sector. |
| 79 | higgsH1parity = Settings::mode("HiggsH1:parity"); |
| 80 | higgsH1eta = Settings::parm("HiggsH1:etaParity"); |
| 81 | higgsH2parity = Settings::mode("HiggsH2:parity"); |
| 82 | higgsH2eta = Settings::parm("HiggsH2:etaParity"); |
| 83 | higgsA3parity = Settings::mode("HiggsA3:parity"); |
| 84 | higgsA3eta = Settings::parm("HiggsA3:etaParity"); |
| 85 | |
| 86 | // If BSM not switched on then H1 should have SM properties. |
| 87 | if (!Settings::flag("Higgs:useBSM")){ |
| 88 | higgsH1parity = 1; |
| 89 | higgsH1eta = 0.; |
| 90 | } |
| 91 | |
| 92 | } |
| 93 | |
| 94 | //********* |
| 95 | |
| 96 | // Set up allowed flux of incoming partons. |
| 97 | // addBeam: set up PDF's that need to be evaluated for the two beams. |
| 98 | // addPair: set up pairs of incoming partons from the two beams. |
| 99 | |
| 100 | bool SigmaProcess::initFlux() { |
| 101 | |
| 102 | // Read in process-specific channel information. |
| 103 | string fluxType = inFlux(); |
| 104 | |
| 105 | // Case with g g incoming state. |
| 106 | if (fluxType == "gg") { |
| 107 | addBeamA(21); |
| 108 | addBeamB(21); |
| 109 | addPair(21, 21); |
| 110 | } |
| 111 | |
| 112 | // Case with q g incoming state. |
| 113 | else if (fluxType == "qg") { |
| 114 | for (int i = -nQuarkIn; i <= nQuarkIn; ++i) { |
| 115 | int idNow = (i == 0) ? 21 : i; |
| 116 | addBeamA(idNow); |
| 117 | addBeamB(idNow); |
| 118 | } |
| 119 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 120 | if (idNow != 0) { |
| 121 | addPair(idNow, 21); |
| 122 | addPair(21, idNow); |
| 123 | } |
| 124 | } |
| 125 | |
| 126 | // Case with q q', q qbar' or qbar qbar' incoming state. |
| 127 | else if (fluxType == "qq") { |
| 128 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 129 | if (idNow != 0) { |
| 130 | addBeamA(idNow); |
| 131 | addBeamB(idNow); |
| 132 | } |
| 133 | for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now) |
| 134 | if (id1Now != 0) |
| 135 | for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now) |
| 136 | if (id2Now != 0) |
| 137 | addPair(id1Now, id2Now); |
| 138 | } |
| 139 | |
| 140 | // Case with q qbar incoming state. |
| 141 | else if (fluxType == "qqbarSame") { |
| 142 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 143 | if (idNow != 0) { |
| 144 | addBeamA(idNow); |
| 145 | addBeamB(idNow); |
| 146 | } |
| 147 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 148 | if (idNow != 0) |
| 149 | addPair(idNow, -idNow); |
| 150 | } |
| 151 | |
| 152 | // Case with f f', f fbar', fbar fbar' incoming state. |
| 153 | else if (fluxType == "ff") { |
| 154 | // If beams are leptons then they are also the colliding partons. |
| 155 | if ( isLeptonA && isLeptonB ) { |
| 156 | addBeamA(idA); |
| 157 | addBeamB(idB); |
| 158 | addPair(idA, idB); |
| 159 | // First beam is lepton and second is hadron. |
| 160 | } else if ( isLeptonA ) { |
| 161 | addBeamA(idA); |
| 162 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 163 | if (idNow != 0) { |
| 164 | addBeamB(idNow); |
| 165 | addPair(idA, idNow); |
| 166 | } |
| 167 | // First beam is hadron and second is lepton. |
| 168 | } else if ( isLeptonB ) { |
| 169 | addBeamB(idB); |
| 170 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 171 | if (idNow != 0) { |
| 172 | addBeamA(idNow); |
| 173 | addPair(idNow, idB); |
| 174 | } |
| 175 | // Hadron beams gives quarks. |
| 176 | } else { |
| 177 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 178 | if (idNow != 0) { |
| 179 | addBeamA(idNow); |
| 180 | addBeamB(idNow); |
| 181 | } |
| 182 | for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now) |
| 183 | if (id1Now != 0) |
| 184 | for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now) |
| 185 | if (id2Now != 0) |
| 186 | addPair(id1Now, id2Now); |
| 187 | } |
| 188 | } |
| 189 | |
| 190 | // Case with f fbar incoming state. |
| 191 | else if (fluxType == "ffbarSame") { |
| 192 | // If beams are antiparticle pair and leptons then also colliding partons. |
| 193 | if ( idA + idB == 0 && isLeptonA ) { |
| 194 | addBeamA(idA); |
| 195 | addBeamB(idB); |
| 196 | addPair(idA, idB); |
| 197 | // Else assume both to be hadrons, for better or worse. |
| 198 | } else { |
| 199 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 200 | if (idNow != 0) { |
| 201 | addBeamA(idNow); |
| 202 | addBeamB(idNow); |
| 203 | } |
| 204 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 205 | if (idNow != 0) |
| 206 | addPair(idNow, -idNow); |
| 207 | } |
| 208 | } |
| 209 | |
| 210 | // Case with f fbar' charged(+-1) incoming state. |
| 211 | else if (fluxType == "ffbarChg") { |
| 212 | // If beams are leptons then also colliding partons. |
| 213 | if ( isLeptonA && isLeptonB && abs( ParticleDataTable::chargeType(idA) |
| 214 | + ParticleDataTable::chargeType(idB) ) == 3 ) { |
| 215 | addBeamA(idA); |
| 216 | addBeamB(idB); |
| 217 | addPair(idA, idB); |
| 218 | // Hadron beams gives quarks. |
| 219 | } else { |
| 220 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 221 | if (idNow != 0) { |
| 222 | addBeamA(idNow); |
| 223 | addBeamB(idNow); |
| 224 | } |
| 225 | for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now) |
| 226 | if (id1Now != 0) |
| 227 | for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now) |
| 228 | if (id2Now != 0 && id1Now * id2Now < 0 |
| 229 | && (abs(id1Now) + abs(id2Now))%2 == 1) addPair(id1Now, id2Now); |
| 230 | } |
| 231 | } |
| 232 | |
| 233 | // Case with f fbar' generic incoming state. |
| 234 | else if (fluxType == "ffbar") { |
| 235 | // If beams are leptons then also colliding partons. |
| 236 | if (isLeptonA && isLeptonB && idA * idB < 0) { |
| 237 | addBeamA(idA); |
| 238 | addBeamB(idB); |
| 239 | addPair(idA, idB); |
| 240 | // Hadron beams gives quarks. |
| 241 | } else { |
| 242 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 243 | if (idNow != 0) { |
| 244 | addBeamA(idNow); |
| 245 | addBeamB(idNow); |
| 246 | } |
| 247 | for (int id1Now = -nQuarkIn; id1Now <= nQuarkIn; ++id1Now) |
| 248 | if (id1Now != 0) |
| 249 | for (int id2Now = -nQuarkIn; id2Now <= nQuarkIn; ++id2Now) |
| 250 | if (id2Now != 0 && id1Now * id2Now < 0) |
| 251 | addPair(id1Now, id2Now); |
| 252 | } |
| 253 | } |
| 254 | |
| 255 | // Case with f gamma incoming state. |
| 256 | else if (fluxType == "fgm") { |
| 257 | // Fermion from incoming side A. |
| 258 | if ( isLeptonA ) { |
| 259 | addBeamA(idA); |
| 260 | addPair(idA, 22); |
| 261 | } else { |
| 262 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 263 | if (idNow != 0) { |
| 264 | addBeamA(idNow); |
| 265 | addPair(idNow, 22); |
| 266 | } |
| 267 | } |
| 268 | // Fermion from incoming side B. |
| 269 | if ( isLeptonB ) { |
| 270 | addBeamB( idB); |
| 271 | addPair(22, idB); |
| 272 | } else { |
| 273 | for (int idNow = -nQuarkIn; idNow <= nQuarkIn; ++idNow) |
| 274 | if (idNow != 0) { |
| 275 | addBeamB(idNow); |
| 276 | addPair(22, idNow); |
| 277 | } |
| 278 | } |
| 279 | // Photons in the beams. |
| 280 | addBeamA(22); |
| 281 | addBeamB(22); |
| 282 | } |
| 283 | |
| 284 | // Case with gamma gamma incoming state. |
| 285 | else if (fluxType == "ggm") { |
| 286 | addBeamA(21); |
| 287 | addBeamA(22); |
| 288 | addBeamB(21); |
| 289 | addBeamB(22); |
| 290 | addPair(21, 22); |
| 291 | addPair(22, 21); |
| 292 | } |
| 293 | |
| 294 | // Case with gamma gamma incoming state. |
| 295 | else if (fluxType == "gmgm") { |
| 296 | addBeamA(22); |
| 297 | addBeamB(22); |
| 298 | addPair(22, 22); |
| 299 | } |
| 300 | |
| 301 | // Unrecognized fluxType is bad sign. Else done. |
| 302 | else { |
| 303 | infoPtr->errorMsg("Error in SigmaProcess::initFlux: " |
| 304 | "unrecognized inFlux type", fluxType); |
| 305 | return false; |
| 306 | } |
| 307 | return true; |
| 308 | |
| 309 | } |
| 310 | |
| 311 | //********* |
| 312 | |
| 313 | // Convolute matrix-element expression(s) with parton flux and K factor. |
| 314 | |
| 315 | double SigmaProcess::sigmaPDF() { |
| 316 | |
| 317 | // Evaluate and store the required parton densities. |
| 318 | for (int j = 0; j < sizeBeamA(); ++j) |
| 319 | inBeamA[j].pdf = beamAPtr->xfHard( inBeamA[j].id, x1Save, Q2FacSave); |
| 320 | for (int j = 0; j < sizeBeamB(); ++j) |
| 321 | inBeamB[j].pdf = beamBPtr->xfHard( inBeamB[j].id, x2Save, Q2FacSave); |
| 322 | |
| 323 | // Loop over allowed incoming channels. |
| 324 | sigmaSumSave = 0.; |
| 325 | for (int i = 0; i < sizePair(); ++i) { |
| 326 | |
| 327 | // Evaluate hard-scattering cross section. Include K factor. |
| 328 | inPair[i].pdfSigma = Kfactor |
| 329 | * sigmaHatWrap(inPair[i].idA, inPair[i].idB); |
| 330 | |
| 331 | // Multiply by respective parton densities. |
| 332 | for (int j = 0; j < sizeBeamA(); ++j) |
| 333 | if (inPair[i].idA == inBeamA[j].id) { |
| 334 | inPair[i].pdfA = inBeamA[j].pdf; |
| 335 | inPair[i].pdfSigma *= inBeamA[j].pdf; |
| 336 | break; |
| 337 | } |
| 338 | for (int j = 0; j < sizeBeamB(); ++j) |
| 339 | if (inPair[i].idB == inBeamB[j].id) { |
| 340 | inPair[i].pdfB = inBeamB[j].pdf; |
| 341 | inPair[i].pdfSigma *= inBeamB[j].pdf; |
| 342 | break; |
| 343 | } |
| 344 | |
| 345 | // Sum for all channels. |
| 346 | sigmaSumSave += inPair[i].pdfSigma; |
| 347 | } |
| 348 | |
| 349 | // Done. |
| 350 | return sigmaSumSave; |
| 351 | |
| 352 | } |
| 353 | |
| 354 | //********* |
| 355 | |
| 356 | // Select incoming parton channel and extract parton densities (resolved). |
| 357 | |
| 358 | void SigmaProcess::pickInState(int id1in, int id2in) { |
| 359 | |
| 360 | // Multiple interactions: partons already selected. |
| 361 | if (id1in != 0 && id2in != 0) { |
| 362 | id1 = id1in; |
| 363 | id2 = id2in; |
| 364 | } |
| 365 | |
| 366 | // Pick channel. Extract channel flavours and pdf's. |
| 367 | double sigmaRand = sigmaSumSave * Rndm::flat(); |
| 368 | for (int i = 0; i < sizePair(); ++i) { |
| 369 | sigmaRand -= inPair[i].pdfSigma; |
| 370 | if (sigmaRand <= 0.) { |
| 371 | id1 = inPair[i].idA; |
| 372 | id2 = inPair[i].idB; |
| 373 | pdf1Save = inPair[i].pdfA; |
| 374 | pdf2Save = inPair[i].pdfB; |
| 375 | break; |
| 376 | } |
| 377 | } |
| 378 | |
| 379 | } |
| 380 | |
| 381 | //********* |
| 382 | |
| 383 | // Evaluate weight for W decay distribution in t -> W b -> f fbar b. |
| 384 | |
| 385 | double SigmaProcess::weightTopDecay( Event& process, int iResBeg, |
| 386 | int iResEnd) { |
| 387 | |
| 388 | // If not pair W d/s/b and mother t then return unit weight. |
| 389 | if (iResEnd - iResBeg != 1) return 1.; |
| 390 | int iW1 = iResBeg; |
| 391 | int iB2 = iResBeg + 1; |
| 392 | int idW1 = process[iW1].idAbs(); |
| 393 | int idB2 = process[iB2].idAbs(); |
| 394 | if (idW1 != 24) { |
| 395 | swap(iW1, iB2); |
| 396 | swap(idW1, idB2); |
| 397 | } |
| 398 | if (idW1 != 24 || (idB2 != 1 && idB2 != 3 && idB2 != 5)) return 1.; |
| 399 | int iT = process[iW1].mother1(); |
| 400 | if (iT <= 0 || process[iT].idAbs() != 6) return 1.; |
| 401 | |
| 402 | // Find sign-matched order of W decay products. |
| 403 | int iF = process[iW1].daughter1(); |
| 404 | int iFbar = process[iW1].daughter2(); |
| 405 | if (iFbar - iF != 1) return 1.; |
| 406 | if (process[iT].id() * process[iF].id() < 0) swap(iF, iFbar); |
| 407 | |
| 408 | // Weight and maximum weight. |
| 409 | double wt = (process[iT].p() * process[iFbar].p()) |
| 410 | * (process[iF].p() * process[iB2].p()); |
| 411 | double wtMax = ( pow4(process[iT].m()) - pow4(process[iW1].m()) ) / 8.; |
| 412 | |
| 413 | // Done. |
| 414 | return wt / wtMax; |
| 415 | |
| 416 | } |
| 417 | |
| 418 | |
| 419 | //********* |
| 420 | |
| 421 | // Evaluate weight for Z0/W+- decay distributions in H -> Z0/W+ Z0/W- -> 4f. |
| 422 | |
| 423 | double SigmaProcess::weightHiggsDecay( Event& process, int iResBeg, |
| 424 | int iResEnd) { |
| 425 | |
| 426 | // If not pair Z0 Z0 or W+ W- then return unit weight. |
| 427 | if (iResEnd - iResBeg != 1) return 1.; |
| 428 | int iZW1 = iResBeg; |
| 429 | int iZW2 = iResBeg + 1; |
| 430 | int idZW1 = process[iZW1].id(); |
| 431 | int idZW2 = process[iZW2].id(); |
| 432 | if (idZW1 < 0) { |
| 433 | swap(iZW1, iZW2); |
| 434 | swap(idZW1, idZW2); |
| 435 | } |
| 436 | if ( (idZW1 != 23 || idZW2 != 23) && (idZW1 != 24 || idZW2 != -24) ) |
| 437 | return 1.; |
| 438 | |
| 439 | // If mother is not Higgs then return unit weight. |
| 440 | int iH = process[iZW1].mother1(); |
| 441 | if (iH <= 0) return 1.; |
| 442 | int idH = process[iH].id(); |
| 443 | if (idH != 25 && idH != 35 && idH !=36) return 1.; |
| 444 | |
| 445 | // Parameters depend on Higgs type: H0(H_1), H^0(H_2) or A^0(H_3). |
| 446 | int higgsParity = higgsH1parity; |
| 447 | double higgsEta = higgsH1eta; |
| 448 | if (idH == 35) { |
| 449 | higgsParity = higgsH2parity; |
| 450 | higgsEta = higgsH2eta; |
| 451 | } else if (idH == 36) { |
| 452 | higgsParity = higgsA3parity; |
| 453 | higgsEta = higgsA3eta; |
| 454 | } |
| 455 | |
| 456 | // Option with isotropic decays. |
| 457 | if (higgsParity == 0) return 1.; |
| 458 | |
| 459 | // Maximum and initial weight. |
| 460 | double wtMax = pow4(process[iH].m()); |
| 461 | double wt = wtMax; |
| 462 | |
| 463 | // Find sign-matched order of Z0/W+- decay products. |
| 464 | int i3 = process[iZW1].daughter1(); |
| 465 | int i4 = process[iZW1].daughter2(); |
| 466 | if (process[i3].id() < 0) swap( i3, i4); |
| 467 | int i5 = process[iZW2].daughter1(); |
| 468 | int i6 = process[iZW2].daughter2(); |
| 469 | if (process[i5].id() < 0) swap( i5, i6); |
| 470 | |
| 471 | // Evaluate four-vector products and find masses.. |
| 472 | double p35 = 2. * process[i3].p() * process[i5].p(); |
| 473 | double p36 = 2. * process[i3].p() * process[i6].p(); |
| 474 | double p45 = 2. * process[i4].p() * process[i5].p(); |
| 475 | double p46 = 2. * process[i4].p() * process[i6].p(); |
| 476 | double p34 = 2. * process[i3].p() * process[i4].p(); |
| 477 | double p56 = 2. * process[i5].p() * process[i6].p(); |
| 478 | double mZW1 = process[iZW1].m(); |
| 479 | double mZW2 = process[iZW2].m(); |
| 480 | |
| 481 | // For mixed CP states need epsilon product and gauge boson masses. |
| 482 | double epsilonProd = 0.; |
| 483 | if (higgsParity == 3) { |
| 484 | double p[4][4]; |
| 485 | for (int i = 0; i < 4; ++i) { |
| 486 | int ii = i3; |
| 487 | if (i == 1) ii = i4; |
| 488 | if (i == 2) ii = i5; |
| 489 | if (i == 3) ii = i6; |
| 490 | p[i][0] = process[ii].e(); |
| 491 | p[i][1] = process[ii].px(); |
| 492 | p[i][2] = process[ii].py(); |
| 493 | p[i][3] = process[ii].pz(); |
| 494 | } |
| 495 | epsilonProd |
| 496 | = p[0][0]*p[1][1]*p[2][2]*p[3][3] - p[0][0]*p[1][1]*p[2][3]*p[3][2] |
| 497 | - p[0][0]*p[1][2]*p[2][1]*p[3][3] + p[0][0]*p[1][2]*p[2][3]*p[3][1] |
| 498 | + p[0][0]*p[1][3]*p[2][1]*p[3][2] - p[0][0]*p[1][3]*p[2][2]*p[3][1] |
| 499 | - p[0][1]*p[1][0]*p[2][2]*p[3][3] + p[0][1]*p[1][0]*p[2][3]*p[3][2] |
| 500 | + p[0][1]*p[1][2]*p[2][0]*p[3][3] - p[0][1]*p[1][2]*p[2][3]*p[3][0] |
| 501 | - p[0][1]*p[1][3]*p[2][0]*p[3][2] + p[0][1]*p[1][3]*p[2][2]*p[3][0] |
| 502 | + p[0][2]*p[1][0]*p[2][1]*p[3][3] - p[0][2]*p[1][0]*p[2][3]*p[3][1] |
| 503 | - p[0][2]*p[1][1]*p[2][0]*p[3][3] + p[0][2]*p[1][1]*p[2][3]*p[3][0] |
| 504 | + p[0][2]*p[1][3]*p[2][0]*p[3][1] - p[0][2]*p[1][3]*p[2][1]*p[3][0] |
| 505 | - p[0][3]*p[1][0]*p[2][1]*p[3][2] + p[0][3]*p[1][0]*p[2][2]*p[3][1] |
| 506 | + p[0][3]*p[1][1]*p[2][0]*p[3][2] - p[0][3]*p[1][1]*p[2][2]*p[3][0] |
| 507 | - p[0][3]*p[1][2]*p[2][0]*p[3][1] + p[0][3]*p[1][2]*p[2][1]*p[3][0]; |
| 508 | } |
| 509 | |
| 510 | // Z0 Z0 decay: vector and axial couplings of two fermion pairs. |
| 511 | if (idZW1 == 23) { |
| 512 | double vf1 = CoupEW::vf(process[i3].idAbs()); |
| 513 | double af1 = CoupEW::af(process[i3].idAbs()); |
| 514 | double vf2 = CoupEW::vf(process[i5].idAbs()); |
| 515 | double af2 = CoupEW::af(process[i5].idAbs()); |
| 516 | double va12asym = 4. * vf1 * af1 * vf2 * af2 |
| 517 | / ( (vf1*vf1 + af1*af1) * (vf2*vf2 + af2*af2) ); |
| 518 | double etaMod = higgsEta / pow2( ParticleDataTable::m0(23) ); |
| 519 | |
| 520 | // Normal CP-even decay. |
| 521 | if (higgsParity == 1) wt = 8. * (1. + va12asym) * p35 * p46 |
| 522 | + 8. * (1. - va12asym) * p36 * p45; |
| 523 | |
| 524 | // CP-odd decay (normal for A0(H_3)). |
| 525 | else if (higgsParity == 2) wt = ( pow2(p35 + p46) |
| 526 | + pow2(p36 + p45) - 2. * p34 * p56 |
| 527 | - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56) |
| 528 | + va12asym * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ) |
| 529 | / (1. + va12asym); |
| 530 | |
| 531 | // Mixed CP states. |
| 532 | else wt = 32. * ( 0.25 * ( (1. + va12asym) * p35 * p46 |
| 533 | + (1. - va12asym) * p36 * p45 ) - 0.5 * etaMod * epsilonProd |
| 534 | * ( (1. + va12asym) * (p35 + p46) - (1. - va12asym) * (p36 + p45) ) |
| 535 | + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56) |
| 536 | - 2. * pow2(p35 * p46 - p36 * p45) |
| 537 | + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45)) |
| 538 | + va12asym * p34 * p56 * (p35 + p36 - p45 - p46) |
| 539 | * (p35 + p45 - p36 - p46) ) ) / ( 1. + 2. * etaMod * mZW1 * mZW2 |
| 540 | + 2. * pow2(etaMod * mZW1 * mZW2) * (1. + va12asym) ); |
| 541 | |
| 542 | // W+ W- decay. |
| 543 | } else if (idZW1 == 24) { |
| 544 | double etaMod = higgsEta / pow2( ParticleDataTable::m0(24) ); |
| 545 | |
| 546 | // Normal CP-even decay. |
| 547 | if (higgsParity == 1) wt = 16. * p35 * p46; |
| 548 | |
| 549 | // CP-odd decay (normal for A0(H_3)). |
| 550 | else if (higgsParity == 2) wt = 0.5 * ( pow2(p35 + p46) |
| 551 | + pow2(p36 + p45) - 2. * p34 * p56 |
| 552 | - 2. * pow2(p35 * p46 - p36 * p45) / (p34 * p56) |
| 553 | + (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ); |
| 554 | |
| 555 | // Mixed CP states. |
| 556 | else wt = 32. * ( 0.25 * 2. * p35 * p46 |
| 557 | - 0.5 * etaMod * epsilonProd * 2. * (p35 + p46) |
| 558 | + 0.0625 * etaMod * etaMod * (-2. * pow2(p34 * p56) |
| 559 | - 2. * pow2(p35 * p46 - p36 * p45) |
| 560 | + p34 * p56 * (pow2(p35 + p46) + pow2(p36 + p45)) |
| 561 | + p34 * p56 * (p35 + p36 - p45 - p46) * (p35 + p45 - p36 - p46) ) ) |
| 562 | / ( 1. * 2. * etaMod * mZW1 * mZW2 + 2. * pow2(etaMod * mZW1 * mZW2) ); |
| 563 | } |
| 564 | |
| 565 | // Done. |
| 566 | return wt / wtMax; |
| 567 | |
| 568 | } |
| 569 | |
| 570 | //************************************************************************** |
| 571 | |
| 572 | // The Sigma1Process class. |
| 573 | // Base class for resolved 2 -> 1 cross sections; derived from SigmaProcess. |
| 574 | |
| 575 | //********* |
| 576 | |
| 577 | // Input and complement kinematics for resolved 2 -> 1 process. |
| 578 | |
| 579 | void Sigma1Process::store1Kin( double x1in, double x2in, double sHin) { |
| 580 | |
| 581 | // Default value only sensible for these processes. |
| 582 | swapTU = false; |
| 583 | |
| 584 | // Incoming parton momentum fractions and sHat. |
| 585 | x1Save = x1in; |
| 586 | x2Save = x2in; |
| 587 | sH = sHin; |
| 588 | mH = sqrt(sH); |
| 589 | sH2 = sH * sH; |
| 590 | |
| 591 | // Different options for renormalization scale, but normally sHat. |
| 592 | Q2RenSave = renormMultFac * sH; |
| 593 | if (renormScale1 == 2) Q2RenSave = renormFixScale; |
| 594 | |
| 595 | // Different options for factorization scale, but normally sHat. |
| 596 | Q2FacSave = factorMultFac * sH; |
| 597 | if (factorScale1 == 2) Q2RenSave = factorFixScale; |
| 598 | |
| 599 | // Evaluate alpha_strong and alpha_EM. |
| 600 | alpS = alphaSPtr->alphaS(Q2RenSave); |
| 601 | alpEM = alphaEMPtr->alphaEM(Q2RenSave); |
| 602 | |
| 603 | } |
| 604 | |
| 605 | //************************************************************************** |
| 606 | |
| 607 | // The Sigma2Process class. |
| 608 | // Base class for resolved 2 -> 2 cross sections; derived from SigmaProcess. |
| 609 | |
| 610 | //********* |
| 611 | |
| 612 | // Input and complement kinematics for resolved 2 -> 2 process. |
| 613 | |
| 614 | void Sigma2Process::store2Kin( double x1in, double x2in, double sHin, |
| 615 | double tHin, double m3in, double m4in, double runBW3in, double runBW4in) { |
| 616 | |
| 617 | // Default ordering of particles 3 and 4. |
| 618 | swapTU = false; |
| 619 | |
| 620 | // Incoming parton momentum fractions. |
| 621 | x1Save = x1in; |
| 622 | x2Save = x2in; |
| 623 | |
| 624 | // Incoming masses and their squares. |
| 625 | bool masslessKin = (id3Mass() == 0) && (id4Mass() == 0); |
| 626 | if (masslessKin) { |
| 627 | m3 = 0.; |
| 628 | m4 = 0.; |
| 629 | } else { |
| 630 | m3 = m3in; |
| 631 | m4 = m4in; |
| 632 | } |
| 633 | mSave[3] = m3; |
| 634 | mSave[4] = m4; |
| 635 | s3 = m3 * m3; |
| 636 | s4 = m4 * m4; |
| 637 | |
| 638 | // Standard Mandelstam variables and their squares. |
| 639 | sH = sHin; |
| 640 | tH = tHin; |
| 641 | uH = (masslessKin) ? -(sH + tH) : s3 + s4 - (sH + tH); |
| 642 | mH = sqrt(sH); |
| 643 | sH2 = sH * sH; |
| 644 | tH2 = tH * tH; |
| 645 | uH2 = uH * uH; |
| 646 | |
| 647 | // The nominal Breit-Wigner factors with running width. |
| 648 | runBW3 = runBW3in; |
| 649 | runBW4 = runBW4in; |
| 650 | |
| 651 | // Calculate squared transverse momentum. |
| 652 | pT2 = (masslessKin) ? tH * uH / sH : (tH * uH - s3 * s4) / sH; |
| 653 | |
| 654 | // Special case: pick scale as if 2 -> 1 process in disguise. |
| 655 | if (isSChannel()) { |
| 656 | |
| 657 | // Different options for renormalization scale, but normally sHat. |
| 658 | Q2RenSave = renormMultFac * sH; |
| 659 | if (renormScale1 == 2) Q2RenSave = renormFixScale; |
| 660 | |
| 661 | // Different options for factorization scale, but normally sHat. |
| 662 | Q2FacSave = factorMultFac * sH; |
| 663 | if (factorScale1 == 2) Q2RenSave = factorFixScale; |
| 664 | |
| 665 | // Normal case with "true" 2 -> 2. |
| 666 | } else { |
| 667 | |
| 668 | // Different options for renormalization scale. |
| 669 | if (masslessKin) Q2RenSave = (renormScale2 < 4) ? pT2 : sH; |
| 670 | else if (renormScale2 == 1) Q2RenSave = pT2 + min(s3, s4); |
| 671 | else if (renormScale2 == 2) Q2RenSave = sqrt((pT2 + s3) * (pT2 + s4)); |
| 672 | else if (renormScale2 == 3) Q2RenSave = pT2 + 0.5 * (s3 + s4); |
| 673 | else Q2RenSave = sH; |
| 674 | Q2RenSave *= renormMultFac; |
| 675 | if (renormScale2 == 5) Q2RenSave = renormFixScale; |
| 676 | |
| 677 | // Different options for factorization scale. |
| 678 | if (masslessKin) Q2FacSave = (factorScale2 < 4) ? pT2 : sH; |
| 679 | else if (factorScale2 == 1) Q2FacSave = pT2 + min(s3, s4); |
| 680 | else if (factorScale2 == 2) Q2FacSave = sqrt((pT2 + s3) * (pT2 + s4)); |
| 681 | else if (factorScale2 == 3) Q2FacSave = pT2 + 0.5 * (s3 + s4); |
| 682 | else Q2FacSave = sH; |
| 683 | Q2FacSave *= factorMultFac; |
| 684 | if (factorScale2 == 5) Q2FacSave = factorFixScale; |
| 685 | } |
| 686 | |
| 687 | // Evaluate alpha_strong and alpha_EM. |
| 688 | alpS = alphaSPtr->alphaS(Q2RenSave); |
| 689 | alpEM = alphaEMPtr->alphaEM(Q2RenSave); |
| 690 | |
| 691 | } |
| 692 | |
| 693 | //********* |
| 694 | |
| 695 | // As above, special kinematics for multiple interactions. |
| 696 | |
| 697 | void Sigma2Process::store2KinMI( double x1in, double x2in, |
| 698 | double sHin, double tHin, double uHin, double alpSin, double alpEMin, |
| 699 | bool needMasses, double m3in, double m4in) { |
| 700 | |
| 701 | // Default ordering of particles 3 and 4. |
| 702 | swapTU = false; |
| 703 | |
| 704 | // Incoming x values. |
| 705 | x1Save = x1in; |
| 706 | x2Save = x2in; |
| 707 | |
| 708 | // Standard Mandelstam variables and their squares. |
| 709 | sH = sHin; |
| 710 | tH = tHin; |
| 711 | uH = uHin; |
| 712 | mH = sqrt(sH); |
| 713 | sH2 = sH * sH; |
| 714 | tH2 = tH * tH; |
| 715 | uH2 = uH * uH; |
| 716 | |
| 717 | // Strong and electroweak couplings. |
| 718 | alpS = alpSin; |
| 719 | alpEM = alpEMin; |
| 720 | |
| 721 | // Assume vanishing masses. (Will be modified in final kinematics.) |
| 722 | m3 = 0.; |
| 723 | s3 = 0.; |
| 724 | m4 = 0.; |
| 725 | s4 = 0.; |
| 726 | sHBeta = sH; |
| 727 | |
| 728 | // Scattering angle. |
| 729 | cosTheta = (tH - uH) / sH; |
| 730 | sinTheta = 2. * sqrtpos( tH * uH ) / sH; |
| 731 | |
| 732 | // In some cases must use masses and redefine meaning of tHat and uHat. |
| 733 | if (needMasses) { |
| 734 | m3 = m3in; |
| 735 | s3 = m3 * m3; |
| 736 | m4 = m4in; |
| 737 | s4 = m4 * m4; |
| 738 | sHMass = sH - s3 - s4; |
| 739 | sHBeta = sqrtpos(sHMass*sHMass - 4. * s3 * s4); |
| 740 | tH = -0.5 * (sHMass - sHBeta * cosTheta); |
| 741 | uH = -0.5 * (sHMass + sHBeta * cosTheta); |
| 742 | tH2 = tH * tH; |
| 743 | uH2 = uH * uH; |
| 744 | } |
| 745 | |
| 746 | // pT2 with masses (at this stage) included. |
| 747 | pT2Mass = 0.25 * sHBeta * pow2(sinTheta); |
| 748 | |
| 749 | } |
| 750 | |
| 751 | //********* |
| 752 | |
| 753 | // Perform kinematics for a Multiple Interaction. |
| 754 | |
| 755 | bool Sigma2Process::final2KinMI() { |
| 756 | |
| 757 | // Have to set flavours and colours. |
| 758 | setIdColAcol(); |
| 759 | |
| 760 | // Check that masses of outgoing particles not too big. |
| 761 | m3 = ParticleDataTable::m0(idSave[3]); |
| 762 | m4 = ParticleDataTable::m0(idSave[4]); |
| 763 | mH = sqrt(sH); |
| 764 | if (m3 + m4 + MASSMARGIN > mH) return false; |
| 765 | s3 = m3 * m3; |
| 766 | s4 = m4 * m4; |
| 767 | |
| 768 | // Do kinematics of the decay. |
| 769 | double eIn = 0.5 * mH; |
| 770 | double e3 = 0.5 * (sH + s3 - s4) / mH; |
| 771 | double e4 = 0.5 * (sH + s4 - s3) / mH; |
| 772 | double pAbs = sqrtpos( e3*e3 - s3 ); |
| 773 | phi = 2. * M_PI * Rndm::flat(); |
| 774 | double pZ = pAbs * cosTheta; |
| 775 | double pX = pAbs * sinTheta * sin(phi); |
| 776 | double pY = pAbs * sinTheta * cos(phi); |
| 777 | double scale = eIn * sinTheta; |
| 778 | |
| 779 | // Fill particle info. |
| 780 | parton[1] = Particle( idSave[1], -31, 0, 0, 3, 4, colSave[1], acolSave[1], |
| 781 | 0., 0., eIn, eIn, 0., scale); |
| 782 | parton[2] = Particle( idSave[2], -31, 0, 0, 3, 4, colSave[2], acolSave[2], |
| 783 | 0., 0., -eIn, eIn, 0., scale); |
| 784 | parton[3] = Particle( idSave[3], 33, 1, 2, 0, 0, colSave[3], acolSave[3], |
| 785 | pX, pY, pZ, e3, m3, scale); |
| 786 | parton[4] = Particle( idSave[4], 33, 1, 2, 0, 0, colSave[4], acolSave[4], |
| 787 | -pX, -pY, -pZ, e4, m4, scale); |
| 788 | |
| 789 | // Boost particles from subprocess rest frame to event rest frame. |
| 790 | double betaZ = (x1Save - x2Save) / (x1Save + x2Save); |
| 791 | for (int i = 1; i <= 4; ++i) parton[i].bst(0., 0., betaZ); |
| 792 | |
| 793 | // Done. |
| 794 | return true; |
| 795 | |
| 796 | } |
| 797 | |
| 798 | //************************************************************************** |
| 799 | |
| 800 | // The Sigma3Process class. |
| 801 | // Base class for resolved 2 -> 3 cross sections; derived from SigmaProcess. |
| 802 | |
| 803 | //********* |
| 804 | |
| 805 | // Input and complement kinematics for resolved 2 -> 3 process. |
| 806 | |
| 807 | void Sigma3Process::store3Kin( double x1in, double x2in, double sHin, |
| 808 | Vec4 p3cmIn, Vec4 p4cmIn, Vec4 p5cmIn, double m3in, double m4in, |
| 809 | double m5in, double runBW3in, double runBW4in, double runBW5in) { |
| 810 | |
| 811 | // Default ordering of particles 3 and 4 - not relevant here. |
| 812 | swapTU = false; |
| 813 | |
| 814 | // Incoming parton momentum fractions. |
| 815 | x1Save = x1in; |
| 816 | x2Save = x2in; |
| 817 | |
| 818 | // Incoming masses and their squares. |
| 819 | if (id3Mass() == 0 && id4Mass() == 0 && id5Mass() == 0) { |
| 820 | m3 = 0.; |
| 821 | m4 = 0.; |
| 822 | m5 = 0.; |
| 823 | } else { |
| 824 | m3 = m3in; |
| 825 | m4 = m4in; |
| 826 | m5 = m5in; |
| 827 | } |
| 828 | mSave[3] = m3; |
| 829 | mSave[4] = m4; |
| 830 | mSave[5] = m5; |
| 831 | s3 = m3 * m3; |
| 832 | s4 = m4 * m4; |
| 833 | s5 = m5 * m5; |
| 834 | |
| 835 | // Standard Mandelstam variables and four-momenta in rest frame. |
| 836 | sH = sHin; |
| 837 | mH = sqrt(sH); |
| 838 | p3cm = p3cmIn; |
| 839 | p4cm = p4cmIn; |
| 840 | p5cm = p5cmIn; |
| 841 | |
| 842 | // The nominal Breit-Wigner factors with running width. |
| 843 | runBW3 = runBW3in; |
| 844 | runBW4 = runBW4in; |
| 845 | runBW5 = runBW5in; |
| 846 | |
| 847 | // Special case: pick scale as if 2 -> 1 process in disguise. |
| 848 | if (isSChannel()) { |
| 849 | |
| 850 | // Different options for renormalization scale, but normally sHat. |
| 851 | Q2RenSave = renormMultFac * sH; |
| 852 | if (renormScale1 == 2) Q2RenSave = renormFixScale; |
| 853 | |
| 854 | // Different options for factorization scale, but normally sHat. |
| 855 | Q2FacSave = factorMultFac * sH; |
| 856 | if (factorScale1 == 2) Q2RenSave = factorFixScale; |
| 857 | |
| 858 | // "Normal" 2 -> 3 processes, i.e. not vector boson fusion. |
| 859 | } else if ( idTchan1() != 23 && idTchan1() != 24 && idTchan2() != 23 |
| 860 | && idTchan1() != 24 ) { |
| 861 | double mT3S = s3 + p3cm.pT2(); |
| 862 | double mT4S = s4 + p4cm.pT2(); |
| 863 | double mT5S = s5 + p5cm.pT2(); |
| 864 | |
| 865 | // Different options for renormalization scale. |
| 866 | if (renormScale3 == 1) Q2RenSave = min( mT3S, min(mT4S, mT5S) ); |
| 867 | else if (renormScale3 == 2) Q2RenSave = sqrt( mT3S * mT4S * mT5S |
| 868 | / max( mT3S, max(mT4S, mT5S) ) ); |
| 869 | else if (renormScale3 == 3) Q2RenSave = pow( mT3S * mT4S * mT5S, |
| 870 | 1./3. ); |
| 871 | else if (renormScale3 == 4) Q2RenSave = (mT3S * mT4S * mT5S) / 3.; |
| 872 | else Q2RenSave = sH; |
| 873 | Q2RenSave *= renormMultFac; |
| 874 | if (renormScale3 == 6) Q2RenSave = renormFixScale; |
| 875 | |
| 876 | // Different options for factorization scale. |
| 877 | if (factorScale3 == 1) Q2FacSave = min( mT3S, min(mT4S, mT5S) ); |
| 878 | else if (factorScale3 == 2) Q2FacSave = sqrt( mT3S * mT4S * mT5S |
| 879 | / max( mT3S, max(mT4S, mT5S) ) ); |
| 880 | else if (factorScale3 == 3) Q2FacSave = pow( mT3S * mT4S * mT5S, |
| 881 | 1./3. ); |
| 882 | else if (factorScale3 == 4) Q2FacSave = (mT3S * mT4S * mT5S) / 3.; |
| 883 | else Q2FacSave = sH; |
| 884 | Q2RenSave *= factorMultFac; |
| 885 | if (factorScale3 == 6) Q2FacSave = factorFixScale; |
| 886 | |
| 887 | // Vector boson fusion 2 -> 3 processes; recoils in positions 4 and 5. |
| 888 | } else { |
| 889 | double sV4 = pow2( ParticleDataTable::m0(idTchan1()) ); |
| 890 | double sV5 = pow2( ParticleDataTable::m0(idTchan2()) ); |
| 891 | double mT3S = s3 + p3cm.pT2(); |
| 892 | double mTV4S = sV4 + p4cm.pT2(); |
| 893 | double mTV5S = sV5 + p5cm.pT2(); |
| 894 | |
| 895 | // Different options for renormalization scale. |
| 896 | if (renormScale3VV == 1) Q2RenSave = max( sV4, sV5); |
| 897 | else if (renormScale3VV == 2) Q2RenSave = sqrt( mTV4S * mTV5S ); |
| 898 | else if (renormScale3VV == 3) Q2RenSave = pow( mT3S * mTV4S * mTV5S, |
| 899 | 1./3. ); |
| 900 | else if (renormScale3VV == 4) Q2RenSave = (mT3S * mTV4S * mTV5S) / 3.; |
| 901 | else Q2RenSave = sH; |
| 902 | Q2RenSave *= renormMultFac; |
| 903 | if (renormScale3VV == 6) Q2RenSave = renormFixScale; |
| 904 | |
| 905 | // Different options for factorization scale. |
| 906 | if (factorScale3VV == 1) Q2FacSave = max( sV4, sV5); |
| 907 | else if (factorScale3VV == 2) Q2FacSave = sqrt( mTV4S * mTV5S ); |
| 908 | else if (factorScale3VV == 3) Q2FacSave = pow( mT3S * mTV4S * mTV5S, |
| 909 | 1./3. ); |
| 910 | else if (factorScale3VV == 4) Q2FacSave = (mT3S * mTV4S * mTV5S) / 3.; |
| 911 | else Q2FacSave = sH; |
| 912 | Q2RenSave *= factorMultFac; |
| 913 | if (factorScale3VV == 6) Q2FacSave = factorFixScale; |
| 914 | } |
| 915 | |
| 916 | // Evaluate alpha_strong and alpha_EM. |
| 917 | alpS = alphaSPtr->alphaS(Q2RenSave); |
| 918 | alpEM = alphaEMPtr->alphaEM(Q2RenSave); |
| 919 | |
| 920 | } |
| 921 | |
| 922 | //************************************************************************** |
| 923 | |
| 924 | // The SigmaLHAProcess class. |
| 925 | // Wrapper for Les Houches Accord external input; derived from SigmaProcess. |
| 926 | // Note: arbitrary subdivision into PhaseSpaceLHA and SigmaLHAProcess tasks. |
| 927 | |
| 928 | //********* |
| 929 | |
| 930 | // Obtain number of final-state partons from LHA object. |
| 931 | |
| 932 | int SigmaLHAProcess::nFinal() const { |
| 933 | |
| 934 | // At initialization size unknown, so return 0. |
| 935 | if (lhaUpPtr->sizePart() <= 0) return 0; |
| 936 | |
| 937 | // Sum up all particles that has first mother = 1. |
| 938 | int nFin = 0; |
| 939 | for (int i = 3; i < lhaUpPtr->sizePart(); ++i) |
| 940 | if (lhaUpPtr->mother1(i) == 1) ++nFin; |
| 941 | return nFin; |
| 942 | |
| 943 | } |
| 944 | |
| 945 | //************************************************************************** |
| 946 | |
| 947 | } // end namespace Pythia8 |
git://git.uio.no / u / mrichter / AliRoot.git / blame