[1] | 1 | // BeamParticle.cc is a part of the PYTHIA event generator. |
---|
| 2 | // Copyright (C) 2012 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 | // BeamParticle class. |
---|
| 8 | |
---|
| 9 | #include "BeamParticle.h" |
---|
| 10 | |
---|
| 11 | namespace Pythia8 { |
---|
| 12 | |
---|
| 13 | //========================================================================== |
---|
| 14 | |
---|
| 15 | // The BeamParticle class. |
---|
| 16 | |
---|
| 17 | //-------------------------------------------------------------------------- |
---|
| 18 | |
---|
| 19 | // Constants: could be changed here if desired, but normally should not. |
---|
| 20 | // These are of technical nature, as described for each. |
---|
| 21 | |
---|
| 22 | // A lepton that takes (almost) the full beam energy does not leave a remnant. |
---|
| 23 | const double BeamParticle::XMINUNRESOLVED = 1. - 1e-10; |
---|
| 24 | |
---|
| 25 | //-------------------------------------------------------------------------- |
---|
| 26 | |
---|
| 27 | // Initialize data on a beam particle and save pointers. |
---|
| 28 | |
---|
| 29 | void BeamParticle::init( int idIn, double pzIn, double eIn, double mIn, |
---|
| 30 | Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtrIn, |
---|
| 31 | Rndm* rndmPtrIn,PDF* pdfInPtr, PDF* pdfHardInPtr, bool isUnresolvedIn, |
---|
| 32 | StringFlav* flavSelPtrIn) { |
---|
| 33 | |
---|
| 34 | // Store input pointers (and one bool) for future use. |
---|
| 35 | infoPtr = infoPtrIn; |
---|
| 36 | particleDataPtr = particleDataPtrIn; |
---|
| 37 | rndmPtr = rndmPtrIn; |
---|
| 38 | pdfBeamPtr = pdfInPtr; |
---|
| 39 | pdfHardBeamPtr = pdfHardInPtr; |
---|
| 40 | isUnresolvedBeam = isUnresolvedIn; |
---|
| 41 | flavSelPtr = flavSelPtrIn; |
---|
| 42 | |
---|
| 43 | // Maximum quark kind in allowed incoming beam hadrons. |
---|
| 44 | maxValQuark = settings.mode("BeamRemnants:maxValQuark"); |
---|
| 45 | |
---|
| 46 | // Power of (1-x)^power/sqrt(x) for remnant valence quark distribution. |
---|
| 47 | valencePowerMeson = settings.parm("BeamRemnants:valencePowerMeson"); |
---|
| 48 | valencePowerUinP = settings.parm("BeamRemnants:valencePowerUinP"); |
---|
| 49 | valencePowerDinP = settings.parm("BeamRemnants:valencePowerDinP"); |
---|
| 50 | |
---|
| 51 | // Enhancement factor of x of diquark. |
---|
| 52 | valenceDiqEnhance = settings.parm("BeamRemnants:valenceDiqEnhance"); |
---|
| 53 | |
---|
| 54 | // Assume g(x) ~ (1-x)^power/x to constrain companion to sea quark. |
---|
| 55 | companionPower = settings.mode("BeamRemnants:companionPower"); |
---|
| 56 | |
---|
| 57 | // Assume g(x) ~ (1-x)^power/x to constrain companion to sea quark. |
---|
| 58 | companionPower = settings.mode("BeamRemnants:companionPower"); |
---|
| 59 | |
---|
| 60 | // Allow or not more than one valence quark to be kicked out. |
---|
| 61 | allowJunction = settings.flag("BeamRemnants:allowJunction"); |
---|
| 62 | |
---|
| 63 | // For low-mass diffractive system kick out q/g = norm / mass^power. |
---|
| 64 | pickQuarkNorm = settings.parm("Diffraction:pickQuarkNorm"); |
---|
| 65 | pickQuarkPower = settings.parm("Diffraction:pickQuarkPower"); |
---|
| 66 | |
---|
| 67 | // Width of primordial kT distribution in low-mass diffractive systems. |
---|
| 68 | diffPrimKTwidth = settings.parm("Diffraction:primKTwidth"); |
---|
| 69 | |
---|
| 70 | // Suppress large masses of beam remnant in low-mass diffractive systems. |
---|
| 71 | diffLargeMassSuppress = settings.parm("Diffraction:largeMassSuppress"); |
---|
| 72 | |
---|
| 73 | // Store info on the incoming beam. |
---|
| 74 | idBeam = idIn; |
---|
| 75 | initBeamKind(); |
---|
| 76 | pBeam = Vec4( 0., 0., pzIn, eIn); |
---|
| 77 | mBeam = mIn; |
---|
| 78 | |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | //-------------------------------------------------------------------------- |
---|
| 82 | |
---|
| 83 | // Initialize kind and valence flavour content of incoming beam. |
---|
| 84 | // For recognized hadrons one can generate multiparton interactions. |
---|
| 85 | // Dynamic choice of meson valence flavours in newValenceContent below. |
---|
| 86 | |
---|
| 87 | void BeamParticle::initBeamKind() { |
---|
| 88 | |
---|
| 89 | // Reset. |
---|
| 90 | idBeamAbs = abs(idBeam); |
---|
| 91 | isLeptonBeam = false; |
---|
| 92 | isHadronBeam = false; |
---|
| 93 | isMesonBeam = false; |
---|
| 94 | isBaryonBeam = false; |
---|
| 95 | nValKinds = 0; |
---|
| 96 | |
---|
| 97 | // Check for leptons. |
---|
| 98 | if (idBeamAbs > 10 && idBeamAbs < 17) { |
---|
| 99 | nValKinds = 1; |
---|
| 100 | nVal[0] = 1; |
---|
| 101 | idVal[0] = idBeam; |
---|
| 102 | isLeptonBeam = true; |
---|
| 103 | } |
---|
| 104 | |
---|
| 105 | // Done if cannot be lowest-lying hadron state. |
---|
| 106 | if (idBeamAbs < 101 || idBeamAbs > 9999) return; |
---|
| 107 | |
---|
| 108 | // Resolve valence content for assumed Pomeron. |
---|
| 109 | if (idBeamAbs == 990) { |
---|
| 110 | isMesonBeam = true; |
---|
| 111 | nValKinds = 2; |
---|
| 112 | nVal[0] = 1 ; |
---|
| 113 | nVal[1] = 1 ; |
---|
| 114 | newValenceContent(); |
---|
| 115 | |
---|
| 116 | // Resolve valence content for assumed meson. Flunk unallowed codes. |
---|
| 117 | } else if (idBeamAbs < 1000) { |
---|
| 118 | int id1 = idBeamAbs/100; |
---|
| 119 | int id2 = (idBeamAbs/10)%10; |
---|
| 120 | if ( id1 < 1 || id1 > maxValQuark |
---|
| 121 | || id2 < 1 || id2 > maxValQuark ) return; |
---|
| 122 | isMesonBeam = true; |
---|
| 123 | |
---|
| 124 | // Store valence content of a confirmed meson. |
---|
| 125 | nValKinds = 2; |
---|
| 126 | nVal[0] = 1 ; |
---|
| 127 | nVal[1] = 1; |
---|
| 128 | if (id1%2 == 0) { |
---|
| 129 | idVal[0] = id1; |
---|
| 130 | idVal[1] = -id2; |
---|
| 131 | } else { |
---|
| 132 | idVal[0] = id2; |
---|
| 133 | idVal[1] = -id1; |
---|
| 134 | } |
---|
| 135 | newValenceContent(); |
---|
| 136 | |
---|
| 137 | // Resolve valence content for assumed baryon. Flunk unallowed codes. |
---|
| 138 | } else { |
---|
| 139 | int id1 = idBeamAbs/1000; |
---|
| 140 | int id2 = (idBeamAbs/100)%10; |
---|
| 141 | int id3 = (idBeamAbs/10)%10; |
---|
| 142 | if ( id1 < 1 || id1 > maxValQuark || id2 < 1 || id2 > maxValQuark |
---|
| 143 | || id3 < 1 || id3 > maxValQuark) return; |
---|
| 144 | if (id2 > id1 || id3 > id1) return; |
---|
| 145 | isBaryonBeam = true; |
---|
| 146 | |
---|
| 147 | // Store valence content of a confirmed baryon. |
---|
| 148 | nValKinds = 1; idVal[0] = id1; nVal[0] = 1; |
---|
| 149 | if (id2 == id1) ++nVal[0]; |
---|
| 150 | else { |
---|
| 151 | nValKinds = 2; |
---|
| 152 | idVal[1] = id2; |
---|
| 153 | nVal[1] = 1; |
---|
| 154 | } |
---|
| 155 | if (id3 == id1) ++nVal[0]; |
---|
| 156 | else if (id3 == id2) ++nVal[1]; |
---|
| 157 | else { |
---|
| 158 | idVal[nValKinds] = id3; |
---|
| 159 | nVal[nValKinds] = 1; |
---|
| 160 | ++nValKinds; |
---|
| 161 | } |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | // Flip flavours for antimeson or antibaryon, and then done. |
---|
| 165 | if (idBeam < 0) for (int i = 0; i < nValKinds; ++i) idVal[i] = -idVal[i]; |
---|
| 166 | isHadronBeam = true; |
---|
| 167 | Q2ValFracSav = -1.; |
---|
| 168 | |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | //-------------------------------------------------------------------------- |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | // Dynamic choice of meson valence flavours for pi0, K0S, K0L, Pomeron. |
---|
| 175 | |
---|
| 176 | void BeamParticle::newValenceContent() { |
---|
| 177 | |
---|
| 178 | // A pi0 oscillates between d dbar and u ubar. |
---|
| 179 | if (idBeam == 111) { |
---|
| 180 | idVal[0] = (rndmPtr->flat() < 0.5) ? 1 : 2; |
---|
| 181 | idVal[1] = -idVal[0]; |
---|
| 182 | |
---|
| 183 | // A K0S or K0L oscillates between d sbar and s dbar. |
---|
| 184 | } else if (idBeam == 130 || idBeam == 310) { |
---|
| 185 | idVal[0] = (rndmPtr->flat() < 0.5) ? 1 : 3; |
---|
| 186 | idVal[1] = (idVal[0] == 1) ? -3 : -1; |
---|
| 187 | |
---|
| 188 | // For a Pomeron split gluon remnant into d dbar or u ubar. |
---|
| 189 | } else if (idBeam == 990) { |
---|
| 190 | idVal[0] = (rndmPtr->flat() < 0.5) ? 1 : 2; |
---|
| 191 | idVal[1] = -idVal[0]; |
---|
| 192 | |
---|
| 193 | // Other hadrons so far do not require any event-by-event change. |
---|
| 194 | } else return; |
---|
| 195 | |
---|
| 196 | // Propagate change to PDF routine(s). |
---|
| 197 | pdfBeamPtr->newValenceContent( idVal[0], idVal[1]); |
---|
| 198 | if (pdfHardBeamPtr != pdfBeamPtr && pdfHardBeamPtr != 0) |
---|
| 199 | pdfHardBeamPtr->newValenceContent( idVal[0], idVal[1]); |
---|
| 200 | |
---|
| 201 | } |
---|
| 202 | |
---|
| 203 | //-------------------------------------------------------------------------- |
---|
| 204 | |
---|
| 205 | double BeamParticle::xMax(int iSkip) { |
---|
| 206 | |
---|
| 207 | // Minimum requirement on remaining energy > nominal mass for hadron. |
---|
| 208 | double xLeft = 1.; |
---|
| 209 | if (isHadron()) xLeft -= m() / e(); |
---|
| 210 | if (size() == 0) return xLeft; |
---|
| 211 | |
---|
| 212 | // Subtract what was carried away by initiators (to date). |
---|
| 213 | for (int i = 0; i < size(); ++i) |
---|
| 214 | if (i != iSkip && resolved[i].isFromBeam()) xLeft -= resolved[i].x(); |
---|
| 215 | return xLeft; |
---|
| 216 | |
---|
| 217 | } |
---|
| 218 | |
---|
| 219 | //-------------------------------------------------------------------------- |
---|
| 220 | |
---|
| 221 | // Parton distributions, reshaped to take into account previous |
---|
| 222 | // multiparton interactions. By picking a non-negative iSkip value, |
---|
| 223 | // one particular interaction is skipped, as needed for ISR |
---|
| 224 | |
---|
| 225 | double BeamParticle::xfModified(int iSkip, int idIn, double x, double Q2) { |
---|
| 226 | |
---|
| 227 | // Initial values. |
---|
| 228 | idSave = idIn; |
---|
| 229 | iSkipSave = iSkip; |
---|
| 230 | xqVal = 0.; |
---|
| 231 | xqgSea = 0.; |
---|
| 232 | xqCompSum = 0.; |
---|
| 233 | |
---|
| 234 | // Fast procedure for first interaction. |
---|
| 235 | if (size() == 0) { |
---|
| 236 | if (x >= 1.) return 0.; |
---|
| 237 | bool canBeVal = false; |
---|
| 238 | for (int i = 0; i < nValKinds; ++i) |
---|
| 239 | if (idIn == idVal[i]) canBeVal = true; |
---|
| 240 | if (canBeVal) { |
---|
| 241 | xqVal = xfVal( idIn, x, Q2); |
---|
| 242 | xqgSea = xfSea( idIn, x, Q2); |
---|
| 243 | } |
---|
| 244 | else xqgSea = xf( idIn, x, Q2); |
---|
| 245 | |
---|
| 246 | // More complicated procedure for non-first interaction. |
---|
| 247 | } else { |
---|
| 248 | |
---|
| 249 | // Sum up the x already removed, and check that remaining x is enough. |
---|
| 250 | double xUsed = 0.; |
---|
| 251 | for (int i = 0; i < size(); ++i) |
---|
| 252 | if (i != iSkip) xUsed += resolved[i].x(); |
---|
| 253 | double xLeft = 1. - xUsed; |
---|
| 254 | if (x >= xLeft) return 0.; |
---|
| 255 | double xRescaled = x / xLeft; |
---|
| 256 | |
---|
| 257 | // Calculate total and remaining amount of x carried by valence quarks. |
---|
| 258 | double xValTot = 0.; |
---|
| 259 | double xValLeft = 0.; |
---|
| 260 | for (int i = 0; i < nValKinds; ++i) { |
---|
| 261 | nValLeft[i] = nVal[i]; |
---|
| 262 | for (int j = 0; j < size(); ++j) |
---|
| 263 | if (j != iSkip && resolved[j].isValence() |
---|
| 264 | && resolved[j].id() == idVal[i]) --nValLeft[i]; |
---|
| 265 | double xValNow = xValFrac(i, Q2); |
---|
| 266 | xValTot += nVal[i] * xValNow; |
---|
| 267 | xValLeft += nValLeft[i] * xValNow; |
---|
| 268 | } |
---|
| 269 | |
---|
| 270 | // Calculate total amount of x carried by unmatched companion quarks. |
---|
| 271 | double xCompAdded = 0.; |
---|
| 272 | for (int i = 0; i < size(); ++i) |
---|
| 273 | if (i != iSkip && resolved[i].isUnmatched()) xCompAdded |
---|
| 274 | += xCompFrac( resolved[i].x() / (xLeft + resolved[i].x()) ) |
---|
| 275 | // Typo warning: extrafactor missing in Skands&Sjostrand article; |
---|
| 276 | // <x> for companion refers to fraction of x left INCLUDING sea quark. |
---|
| 277 | // To be modified further?? |
---|
| 278 | * (1. + resolved[i].x() / xLeft); |
---|
| 279 | |
---|
| 280 | // Calculate total rescaling factor and pdf for sea and gluon. |
---|
| 281 | double rescaleGS = max( 0., (1. - xValLeft - xCompAdded) |
---|
| 282 | / (1. - xValTot) ); |
---|
| 283 | xqgSea = rescaleGS * xfSea( idIn, xRescaled, Q2); |
---|
| 284 | |
---|
| 285 | // Find valence part and rescale it to remaining number of quarks. |
---|
| 286 | for (int i = 0; i < nValKinds; ++i) |
---|
| 287 | if (idIn == idVal[i] && nValLeft[i] > 0) |
---|
| 288 | xqVal = xfVal( idIn, xRescaled, Q2) |
---|
| 289 | * double(nValLeft[i]) / double(nVal[i]); |
---|
| 290 | |
---|
| 291 | // Find companion part, by adding all companion contributions. |
---|
| 292 | for (int i = 0; i < size(); ++i) |
---|
| 293 | if (i != iSkip && resolved[i].id() == -idIn |
---|
| 294 | && resolved[i].isUnmatched()) { |
---|
| 295 | double xsRescaled = resolved[i].x() / (xLeft + resolved[i].x()); |
---|
| 296 | double xcRescaled = x / (xLeft + resolved[i].x()); |
---|
| 297 | double xqCompNow = xCompDist( xcRescaled, xsRescaled); |
---|
| 298 | resolved[i].xqCompanion( xqCompNow); |
---|
| 299 | xqCompSum += xqCompNow; |
---|
| 300 | } |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | // Add total, but only return relevant part for ISR. More cases?? |
---|
| 304 | // Watch out, e.g. g can come from either kind of quark.?? |
---|
| 305 | xqgTot = xqVal + xqgSea + xqCompSum; |
---|
| 306 | if (iSkip >= 0) { |
---|
| 307 | if (resolved[iSkip].isValence()) return xqVal; |
---|
| 308 | if (resolved[iSkip].isUnmatched()) return xqgSea + xqCompSum; |
---|
| 309 | } |
---|
| 310 | return xqgTot; |
---|
| 311 | |
---|
| 312 | } |
---|
| 313 | |
---|
| 314 | //-------------------------------------------------------------------------- |
---|
| 315 | |
---|
| 316 | // Decide whether a quark extracted from the beam is of valence, sea or |
---|
| 317 | // companion kind; in the latter case also pick its companion. |
---|
| 318 | // Assumes xfModified has already been called. |
---|
| 319 | |
---|
| 320 | int BeamParticle::pickValSeaComp() { |
---|
| 321 | |
---|
| 322 | // If parton already has a companion than reset code for this. |
---|
| 323 | int oldCompanion = resolved[iSkipSave].companion(); |
---|
| 324 | if (oldCompanion >= 0) resolved[oldCompanion].companion(-2); |
---|
| 325 | |
---|
| 326 | // Default assignment is sea. |
---|
| 327 | int vsc = -2; |
---|
| 328 | |
---|
| 329 | // For gluons or photons no sense of valence or sea. |
---|
| 330 | if (idSave == 21 || idSave == 22) vsc = -1; |
---|
| 331 | |
---|
| 332 | // For lepton beam assume same-kind lepton inside is valence. |
---|
| 333 | else if (isLeptonBeam && idSave == idBeam) vsc = -3; |
---|
| 334 | |
---|
| 335 | // Decide if valence or sea quark. |
---|
| 336 | else { |
---|
| 337 | double xqRndm = xqgTot * rndmPtr->flat(); |
---|
| 338 | if (xqRndm < xqVal) vsc = -3; |
---|
| 339 | else if (xqRndm < xqVal + xqgSea) vsc = -2; |
---|
| 340 | |
---|
| 341 | // If not either, loop over all possible companion quarks. |
---|
| 342 | else { |
---|
| 343 | xqRndm -= xqVal + xqgSea; |
---|
| 344 | for (int i = 0; i < size(); ++i) |
---|
| 345 | if (i != iSkipSave && resolved[i].id() == -idSave |
---|
| 346 | && resolved[i].isUnmatched()) { |
---|
| 347 | xqRndm -= resolved[i].xqCompanion(); |
---|
| 348 | if (xqRndm < 0.) vsc = i; |
---|
| 349 | break; |
---|
| 350 | } |
---|
| 351 | } |
---|
| 352 | } |
---|
| 353 | |
---|
| 354 | // Bookkeep assignment; for sea--companion pair both ways. |
---|
| 355 | resolved[iSkipSave].companion(vsc); |
---|
| 356 | if (vsc >= 0) resolved[vsc].companion(iSkipSave); |
---|
| 357 | |
---|
| 358 | // Done; return code for choice (to distinguish valence/sea in Info). |
---|
| 359 | return vsc; |
---|
| 360 | |
---|
| 361 | } |
---|
| 362 | |
---|
| 363 | //-------------------------------------------------------------------------- |
---|
| 364 | |
---|
| 365 | // Fraction of hadron momentum sitting in a valence quark distribution. |
---|
| 366 | // Based on hardcoded parametrizations of CTEQ 5L numbers. |
---|
| 367 | |
---|
| 368 | double BeamParticle::xValFrac(int j, double Q2) { |
---|
| 369 | |
---|
| 370 | // Only recalculate when required. |
---|
| 371 | if (Q2 != Q2ValFracSav) { |
---|
| 372 | Q2ValFracSav = Q2; |
---|
| 373 | |
---|
| 374 | // Q2-dependence of log-log form; assume fixed Lambda = 0.2. |
---|
| 375 | double llQ2 = log( log( max( 1., Q2) / 0.04 )); |
---|
| 376 | |
---|
| 377 | // Fractions carried by u and d in proton. |
---|
| 378 | uValInt = 0.48 / (1. + 1.56 * llQ2); |
---|
| 379 | dValInt = 0.385 / (1. + 1.60 * llQ2); |
---|
| 380 | } |
---|
| 381 | |
---|
| 382 | // Baryon with three different quark kinds: (2 * u + d) / 3 of proton. |
---|
| 383 | if (isBaryonBeam && nValKinds == 3) return (2. * uValInt + dValInt) / 3.; |
---|
| 384 | |
---|
| 385 | // Baryon with one or two identical: like d or u of proton. |
---|
| 386 | if (isBaryonBeam && nVal[j] == 1) return dValInt; |
---|
| 387 | if (isBaryonBeam && nVal[j] == 2) return uValInt; |
---|
| 388 | |
---|
| 389 | // Meson: (2 * u + d) / 2 of proton so same total valence quark fraction. |
---|
| 390 | return 0.5 * (2. * uValInt + dValInt); |
---|
| 391 | |
---|
| 392 | } |
---|
| 393 | |
---|
| 394 | //-------------------------------------------------------------------------- |
---|
| 395 | |
---|
| 396 | // The momentum integral of a companion quark, with its partner at x_s, |
---|
| 397 | // using an approximate gluon density like (1 - x_g)^power / x_g. |
---|
| 398 | // The value corresponds to an unrescaled range between 0 and 1 - x_s. |
---|
| 399 | |
---|
| 400 | double BeamParticle::xCompFrac(double xs) { |
---|
| 401 | |
---|
| 402 | // Select case by power of gluon (1-x_g) shape. |
---|
| 403 | switch (companionPower) { |
---|
| 404 | |
---|
| 405 | case 0: |
---|
| 406 | return xs * ( 5. + xs * (-9. - 2. * xs * (-3. + xs)) + 3. * log(xs) ) |
---|
| 407 | / ( (-1. + xs) * (2. + xs * (-1. + 2. * xs)) ); |
---|
| 408 | |
---|
| 409 | case 1: |
---|
| 410 | return -1. -3. * xs + ( 2. * pow2(-1. + xs) * (1. + xs + xs*xs)) |
---|
| 411 | / ( 2. + xs*xs * (xs - 3.) + 3. * xs * log(xs) ); |
---|
| 412 | |
---|
| 413 | case 2: |
---|
| 414 | return xs * ( (1. - xs) * (19. + xs * (43. + 4. * xs)) |
---|
| 415 | + 6. * log(xs) * (1. + 6. * xs + 4.*xs*xs) ) / |
---|
| 416 | ( 4. * ( (xs - 1.) * (1. + xs * (4. + xs) ) |
---|
| 417 | - 3. * xs * log(xs) * (1 + xs) ) ); |
---|
| 418 | |
---|
| 419 | case 3: |
---|
| 420 | return 3. * xs * ( (xs - 1.) * (7. + xs * (28. + 13. * xs)) |
---|
| 421 | - 2. * log(xs) * (1. + xs * (9. + 2. * xs * (6. + xs))) ) |
---|
| 422 | / ( 4. + 27. * xs - 31. * pow3(xs) |
---|
| 423 | + 6. * xs * log(xs) * (3. + 2. * xs * (3.+xs)) ); |
---|
| 424 | |
---|
| 425 | default: |
---|
| 426 | return ( -9. * xs * (xs*xs - 1.) * (5. + xs * (24. + xs)) + 12. * xs |
---|
| 427 | * log(xs) * (1. + 2. * xs) * (1. + 2. * xs * (5. + 2. * xs)) ) |
---|
| 428 | / ( 8. * (1. + 2. * xs) * ((xs - 1.) * (1. + xs * (10. + xs)) |
---|
| 429 | - 6. * xs * log(xs) * (1. + xs)) ); |
---|
| 430 | |
---|
| 431 | } |
---|
| 432 | } |
---|
| 433 | |
---|
| 434 | //-------------------------------------------------------------------------- |
---|
| 435 | |
---|
| 436 | // The x*f pdf of a companion quark at x_c, with its sea partner at x_s, |
---|
| 437 | // using an approximate gluon density like (1 - x_g)^power / x_g. |
---|
| 438 | // The value corresponds to an unrescaled range between 0 and 1 - x_s. |
---|
| 439 | |
---|
| 440 | double BeamParticle::xCompDist(double xc, double xs) { |
---|
| 441 | |
---|
| 442 | // Mother gluon momentum fraction. Check physical limit. |
---|
| 443 | double xg = xc + xs; |
---|
| 444 | if (xg > 1.) return 0.; |
---|
| 445 | |
---|
| 446 | // Common factor, including splitting kernel and part of gluon density |
---|
| 447 | // (and that it is x_c * f that is coded). |
---|
| 448 | double fac = 3. * xc * xs * (xc*xc + xs*xs) / pow4(xg); |
---|
| 449 | |
---|
| 450 | // Select case by power of gluon (1-x_g) shape. |
---|
| 451 | switch (companionPower) { |
---|
| 452 | |
---|
| 453 | case 0: |
---|
| 454 | return fac / ( 2. - xs * (3. - xs * (3. - 2. * xs)) ); |
---|
| 455 | |
---|
| 456 | case 1: |
---|
| 457 | return fac * (1. - xg) / ( 2. + xs*xs * (-3. + xs) + 3. * xs * log(xs) ); |
---|
| 458 | |
---|
| 459 | case 2: |
---|
| 460 | return fac * pow2(1. - xg) / ( 2. * ((1. - xs) * (1. + xs * (4. + xs)) |
---|
| 461 | + 3. * xs * (1. + xs) * log(xs)) ); |
---|
| 462 | |
---|
| 463 | case 3: |
---|
| 464 | return fac * pow3(1. - xg) * 2. / ( 4. + 27. * xs - 31. * pow3(xs) |
---|
| 465 | + 6. * xs * log(xs) * (3. + 2. * xs * (3. + xs)) ); |
---|
| 466 | |
---|
| 467 | default: |
---|
| 468 | return fac * pow4(1. - xg) / ( 2. * (1. + 2. * xs) * ((1. - xs) |
---|
| 469 | * (1. + xs * (10. + xs)) + 6. * xs * log(xs) * (1. + xs)) ); |
---|
| 470 | |
---|
| 471 | } |
---|
| 472 | } |
---|
| 473 | |
---|
| 474 | //-------------------------------------------------------------------------- |
---|
| 475 | |
---|
| 476 | // Add required extra remnant flavour content. Also initial colours. |
---|
| 477 | |
---|
| 478 | bool BeamParticle::remnantFlavours(Event& event) { |
---|
| 479 | |
---|
| 480 | // A baryon will have a junction, unless a diquark is formed later. |
---|
| 481 | hasJunctionBeam = (isBaryon()); |
---|
| 482 | |
---|
| 483 | // Store how many hard-scattering partons were removed from beam. |
---|
| 484 | nInit = size(); |
---|
| 485 | |
---|
| 486 | // Find remaining valence quarks. |
---|
| 487 | for (int i = 0; i < nValKinds; ++i) { |
---|
| 488 | nValLeft[i] = nVal[i]; |
---|
| 489 | for (int j = 0; j < nInit; ++j) if (resolved[j].isValence() |
---|
| 490 | && resolved[j].id() == idVal[i]) --nValLeft[i]; |
---|
| 491 | // Add remaining valence quarks to record. Partly temporary values. |
---|
| 492 | for (int k = 0; k < nValLeft[i]; ++k) append(0, idVal[i], 0., -3); |
---|
| 493 | } |
---|
| 494 | |
---|
| 495 | // If at least two valence quarks left in baryon remnant then form diquark. |
---|
| 496 | int nInitPlusVal = size(); |
---|
| 497 | if (isBaryon() && nInitPlusVal - nInit >= 2) { |
---|
| 498 | |
---|
| 499 | // If three, pick two at random to form diquark, else trivial. |
---|
| 500 | int iQ1 = nInit; |
---|
| 501 | int iQ2 = nInit + 1; |
---|
| 502 | if (nInitPlusVal - nInit == 3) { |
---|
| 503 | double pickDq = 3. * rndmPtr->flat(); |
---|
| 504 | if (pickDq > 1.) iQ2 = nInit + 2; |
---|
| 505 | if (pickDq > 2.) iQ1 = nInit + 1; |
---|
| 506 | } |
---|
| 507 | |
---|
| 508 | // Pick spin 0 or 1 according to SU(6) wave function factors. |
---|
| 509 | int idDq = flavSelPtr->makeDiquark( resolved[iQ1].id(), |
---|
| 510 | resolved[iQ2].id(), idBeam); |
---|
| 511 | |
---|
| 512 | // Overwrite with diquark flavour and remove one slot. No more junction. |
---|
| 513 | resolved[iQ1].id(idDq); |
---|
| 514 | if (nInitPlusVal - nInit == 3 && iQ2 == nInit + 1) |
---|
| 515 | resolved[nInit + 1].id( resolved[nInit + 2].id() ); |
---|
| 516 | resolved.pop_back(); |
---|
| 517 | hasJunctionBeam = false; |
---|
| 518 | } |
---|
| 519 | |
---|
| 520 | // Find companion quarks to unmatched sea quarks. |
---|
| 521 | for (int i = 0; i < nInit; ++i) |
---|
| 522 | if (resolved[i].isUnmatched()) { |
---|
| 523 | |
---|
| 524 | // Add companion quark to record; and bookkeep both ways. |
---|
| 525 | append(0, -resolved[i].id(), 0., i); |
---|
| 526 | resolved[i].companion(size() - 1); |
---|
| 527 | } |
---|
| 528 | |
---|
| 529 | // If no other remnants found, add a gluon or photon to carry momentum. |
---|
| 530 | if (size() == nInit) { |
---|
| 531 | int idRemnant = (isHadronBeam) ? 21 : 22; |
---|
| 532 | append(0, idRemnant, 0., -1); |
---|
| 533 | } |
---|
| 534 | |
---|
| 535 | // Set initiator and remnant masses. |
---|
| 536 | for (int i = 0; i < size(); ++i) { |
---|
| 537 | if (i < nInit) resolved[i].m(0.); |
---|
| 538 | else resolved[i].m( particleDataPtr->m0( resolved[i].id() ) ); |
---|
| 539 | } |
---|
| 540 | |
---|
| 541 | // For debug purposes: reject beams with resolved junction topology. |
---|
| 542 | if (hasJunctionBeam && !allowJunction) return false; |
---|
| 543 | |
---|
| 544 | // Pick initial colours for remnants. |
---|
| 545 | for (int i = nInit; i < size(); ++i) { |
---|
| 546 | int colType = particleDataPtr->colType( resolved[i].id() ); |
---|
| 547 | int col = (colType == 1 || colType == 2) ? event.nextColTag() : 0; |
---|
| 548 | int acol = (colType == -1 || colType == 2) ? event.nextColTag() : 0; |
---|
| 549 | resolved[i].cols( col, acol); |
---|
| 550 | } |
---|
| 551 | |
---|
| 552 | // Done. |
---|
| 553 | return true; |
---|
| 554 | |
---|
| 555 | } |
---|
| 556 | |
---|
| 557 | //-------------------------------------------------------------------------- |
---|
| 558 | |
---|
| 559 | // Correlate all initiators and remnants to make a colour singlet. |
---|
| 560 | |
---|
| 561 | bool BeamParticle::remnantColours(Event& event, vector<int>& colFrom, |
---|
| 562 | vector<int>& colTo) { |
---|
| 563 | |
---|
| 564 | // No colours in lepton beams so no need to do anything. |
---|
| 565 | if (isLeptonBeam) return true; |
---|
| 566 | |
---|
| 567 | // Copy initiator colour info from the event record to the beam. |
---|
| 568 | for (int i = 0; i < size(); ++i) { |
---|
| 569 | int j = resolved[i].iPos(); |
---|
| 570 | resolved[i].cols( event[j].col(), event[j].acol()); |
---|
| 571 | } |
---|
| 572 | |
---|
| 573 | // Find number and position of valence quarks, of gluons, and |
---|
| 574 | // of sea-companion pairs (counted as gluons) in the beam remnants. |
---|
| 575 | // Skip gluons with same colour as anticolour and rescattering partons. |
---|
| 576 | vector<int> iVal; |
---|
| 577 | vector<int> iGlu; |
---|
| 578 | for (int i = 0; i < size(); ++i) |
---|
| 579 | if (resolved[i].isFromBeam()) { |
---|
| 580 | if ( resolved[i].isValence() ) iVal.push_back(i); |
---|
| 581 | else if ( resolved[i].isCompanion() && resolved[i].companion() > i ) |
---|
| 582 | iGlu.push_back(i); |
---|
| 583 | else if ( resolved[i].id() == 21 |
---|
| 584 | && resolved[i].col() != resolved[i].acol() ) iGlu.push_back(i); |
---|
| 585 | } |
---|
| 586 | |
---|
| 587 | // Pick a valence quark to which gluons are attached. |
---|
| 588 | // Do not resolve quarks in diquark. (More sophisticated??) |
---|
| 589 | int iValSel= iVal[0]; |
---|
| 590 | if (iVal.size() == 2) { |
---|
| 591 | if ( abs(resolved[iValSel].id()) > 10 ) iValSel = iVal[1]; |
---|
| 592 | } else { |
---|
| 593 | double rndmValSel = 3. * rndmPtr->flat(); |
---|
| 594 | if (rndmValSel > 1.) iValSel= iVal[1]; |
---|
| 595 | if (rndmValSel > 2.) iValSel= iVal[2]; |
---|
| 596 | } |
---|
| 597 | |
---|
| 598 | // This valence quark defines initial (anti)colour. |
---|
| 599 | int iBeg = iValSel; |
---|
| 600 | bool hasCol = (resolved[iBeg].col() > 0); |
---|
| 601 | int begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol(); |
---|
| 602 | |
---|
| 603 | // Do random stepping through gluon/(sea+companion) list. |
---|
| 604 | vector<int> iGluRndm; |
---|
| 605 | for (int i = 0; i < int(iGlu.size()); ++i) |
---|
| 606 | iGluRndm.push_back( iGlu[i] ); |
---|
| 607 | for (int iOrder = 0; iOrder < int(iGlu.size()); ++iOrder) { |
---|
| 608 | int iRndm = int( double(iGluRndm.size()) * rndmPtr->flat()); |
---|
| 609 | int iGluSel = iGluRndm[iRndm]; |
---|
| 610 | iGluRndm[iRndm] = iGluRndm[iGluRndm.size() - 1]; |
---|
| 611 | iGluRndm.pop_back(); |
---|
| 612 | |
---|
| 613 | // Find matching anticolour/colour to current colour/anticolour. |
---|
| 614 | int iEnd = iGluSel; |
---|
| 615 | int endCol = (hasCol) ? resolved[iEnd].acol() : resolved[iEnd].col(); |
---|
| 616 | // Not gluon but sea+companion pair: go to other. |
---|
| 617 | if (endCol == 0) { |
---|
| 618 | iEnd = resolved[iEnd].companion(); |
---|
| 619 | endCol = (hasCol) ? resolved[iEnd].acol() : resolved[iEnd].col(); |
---|
| 620 | } |
---|
| 621 | |
---|
| 622 | // Collapse this colour-anticolour pair to the lowest one. |
---|
| 623 | if (begCol < endCol) { |
---|
| 624 | if (hasCol) resolved[iEnd].acol(begCol); |
---|
| 625 | else resolved[iEnd].col(begCol); |
---|
| 626 | colFrom.push_back(endCol); |
---|
| 627 | colTo.push_back(begCol); |
---|
| 628 | } else { |
---|
| 629 | if (hasCol) resolved[iBeg].col(endCol); |
---|
| 630 | else resolved[iBeg].acol(endCol); |
---|
| 631 | colFrom.push_back(begCol); |
---|
| 632 | colTo.push_back(endCol); |
---|
| 633 | } |
---|
| 634 | |
---|
| 635 | // Pick up the other colour of the recent gluon and repeat. |
---|
| 636 | iBeg = iEnd; |
---|
| 637 | begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol(); |
---|
| 638 | // Not gluon but sea+companion pair: go to other. |
---|
| 639 | if (begCol == 0) { |
---|
| 640 | iBeg = resolved[iBeg].companion(); |
---|
| 641 | begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol(); |
---|
| 642 | } |
---|
| 643 | |
---|
| 644 | // At end of gluon/(sea+companion) list. |
---|
| 645 | } |
---|
| 646 | |
---|
| 647 | // Now begin checks, and also finding junction information. |
---|
| 648 | // Loop through remnant partons; isolate all colours and anticolours. |
---|
| 649 | vector<int> colList; |
---|
| 650 | vector<int> acolList; |
---|
| 651 | for (int i = 0; i < size(); ++i) |
---|
| 652 | if ( resolved[i].isFromBeam() ) |
---|
| 653 | if ( resolved[i].col() != resolved[i].acol() ) { |
---|
| 654 | if (resolved[i].col() > 0) colList.push_back( resolved[i].col() ); |
---|
| 655 | if (resolved[i].acol() > 0) acolList.push_back( resolved[i].acol() ); |
---|
| 656 | } |
---|
| 657 | |
---|
| 658 | // Remove all matching colour-anticolour pairs. |
---|
| 659 | bool foundPair = true; |
---|
| 660 | while (foundPair && colList.size() > 0 && acolList.size() > 0) { |
---|
| 661 | foundPair = false; |
---|
| 662 | for (int iCol = 0; iCol < int(colList.size()); ++iCol) { |
---|
| 663 | for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) { |
---|
| 664 | if (acolList[iAcol] == colList[iCol]) { |
---|
| 665 | colList[iCol] = colList.back(); |
---|
| 666 | colList.pop_back(); |
---|
| 667 | acolList[iAcol] = acolList.back(); |
---|
| 668 | acolList.pop_back(); |
---|
| 669 | foundPair = true; |
---|
| 670 | break; |
---|
| 671 | } |
---|
| 672 | } if (foundPair) break; |
---|
| 673 | } |
---|
| 674 | } |
---|
| 675 | |
---|
| 676 | // Usually one unmatched pair left to collapse. |
---|
| 677 | if (colList.size() == 1 && acolList.size() == 1) { |
---|
| 678 | int finalFrom = max( colList[0], acolList[0]); |
---|
| 679 | int finalTo = min( colList[0], acolList[0]); |
---|
| 680 | for (int i = 0; i < size(); ++i) |
---|
| 681 | if ( resolved[i].isFromBeam() ) { |
---|
| 682 | if (resolved[i].col() == finalFrom) resolved[i].col(finalTo); |
---|
| 683 | if (resolved[i].acol() == finalFrom) resolved[i].acol(finalTo); |
---|
| 684 | } |
---|
| 685 | colFrom.push_back(finalFrom); |
---|
| 686 | colTo.push_back(finalTo); |
---|
| 687 | |
---|
| 688 | // Store an (anti)junction when three (anti)coloured daughters. |
---|
| 689 | } else if (hasJunctionBeam && colList.size() == 3 |
---|
| 690 | && acolList.size() == 0) { |
---|
| 691 | event.appendJunction( 1, colList[0], colList[1], colList[2]); |
---|
| 692 | junCol[0] = colList[0]; |
---|
| 693 | junCol[1] = colList[1]; |
---|
| 694 | junCol[2] = colList[2]; |
---|
| 695 | } else if (hasJunctionBeam && acolList.size() == 3 |
---|
| 696 | && colList.size() == 0) { |
---|
| 697 | event.appendJunction( 2, acolList[0], acolList[1], acolList[2]); |
---|
| 698 | junCol[0] = acolList[0]; |
---|
| 699 | junCol[1] = acolList[1]; |
---|
| 700 | junCol[2] = acolList[2]; |
---|
| 701 | |
---|
| 702 | // Any other nonvanishing values indicate failure. |
---|
| 703 | } else if (colList.size() > 0 || acolList.size() > 0) { |
---|
| 704 | infoPtr->errorMsg("Error in BeamParticle::remnantColours: " |
---|
| 705 | "leftover unmatched colours"); |
---|
| 706 | return false; |
---|
| 707 | } |
---|
| 708 | |
---|
| 709 | // Store colour assignment of beam particles. |
---|
| 710 | for (int i = nInit; i < size(); ++i) |
---|
| 711 | event[resolved[i].iPos()].cols( resolved[i].col(), resolved[i].acol() ); |
---|
| 712 | |
---|
| 713 | // Done. |
---|
| 714 | return true; |
---|
| 715 | } |
---|
| 716 | |
---|
| 717 | |
---|
| 718 | //-------------------------------------------------------------------------- |
---|
| 719 | |
---|
| 720 | // Pick unrescaled x values for beam remnant sharing. |
---|
| 721 | |
---|
| 722 | double BeamParticle::xRemnant( int i) { |
---|
| 723 | |
---|
| 724 | double x = 0.; |
---|
| 725 | |
---|
| 726 | // Calculation of x of valence quark or diquark, for latter as sum. |
---|
| 727 | if (resolved[i].isValence()) { |
---|
| 728 | |
---|
| 729 | // Resolve diquark into sum of two quarks. |
---|
| 730 | int id1 = resolved[i].id(); |
---|
| 731 | int id2 = 0; |
---|
| 732 | if (abs(id1) > 10) { |
---|
| 733 | id2 = (id1 > 0) ? (id1/100)%10 : -(((-id1)/100)%10); |
---|
| 734 | id1 = (id1 > 0) ? id1/1000 : -((-id1)/1000); |
---|
| 735 | } |
---|
| 736 | |
---|
| 737 | // Loop over (up to) two quarks; add their contributions. |
---|
| 738 | for (int iId = 0; iId < 2; ++iId) { |
---|
| 739 | int idNow = (iId == 0) ? id1 : id2; |
---|
| 740 | if (idNow == 0) break; |
---|
| 741 | double xPart = 0.; |
---|
| 742 | |
---|
| 743 | // Assume form (1-x)^a / sqrt(x). |
---|
| 744 | double xPow = valencePowerMeson; |
---|
| 745 | if (isBaryonBeam) { |
---|
| 746 | if (nValKinds == 3 || nValKinds == 1) |
---|
| 747 | xPow = (3. * rndmPtr->flat() < 2.) |
---|
| 748 | ? valencePowerUinP : valencePowerDinP ; |
---|
| 749 | else if (nValence(idNow) == 2) xPow = valencePowerUinP; |
---|
| 750 | else xPow = valencePowerDinP; |
---|
| 751 | } |
---|
| 752 | do xPart = pow2( rndmPtr->flat() ); |
---|
| 753 | while ( pow(1. - xPart, xPow) < rndmPtr->flat() ); |
---|
| 754 | |
---|
| 755 | // End loop over (up to) two quarks. Possibly enhancement for diquarks. |
---|
| 756 | x += xPart; |
---|
| 757 | } |
---|
| 758 | if (id2 != 0) x *= valenceDiqEnhance; |
---|
| 759 | |
---|
| 760 | // Calculation of x of sea quark, based on companion association. |
---|
| 761 | } else if (resolved[i].isCompanion()) { |
---|
| 762 | |
---|
| 763 | // Find rescaled x value of companion. |
---|
| 764 | double xLeft = 1.; |
---|
| 765 | for (int iInit = 0; iInit < nInit; ++iInit) |
---|
| 766 | if (resolved[iInit].isFromBeam()) xLeft -= resolved[iInit].x(); |
---|
| 767 | double xCompanion = resolved[ resolved[i].companion() ].x(); |
---|
| 768 | xCompanion /= (xLeft + xCompanion); |
---|
| 769 | |
---|
| 770 | // Now use ansatz q(x; x_c) < N/(x +x_c) to pick x. |
---|
| 771 | do x = pow( xCompanion, rndmPtr->flat()) - xCompanion; |
---|
| 772 | while ( pow( (1. - x - xCompanion) / (1. - xCompanion), companionPower) |
---|
| 773 | * (pow2(x) + pow2(xCompanion)) / pow2(x + xCompanion) |
---|
| 774 | < rndmPtr->flat() ); |
---|
| 775 | |
---|
| 776 | // Else, rarely, a single gluon remnant, so value does not matter. |
---|
| 777 | } else x = 1.; |
---|
| 778 | return x; |
---|
| 779 | |
---|
| 780 | } |
---|
| 781 | |
---|
| 782 | //-------------------------------------------------------------------------- |
---|
| 783 | |
---|
| 784 | // Print the list of resolved partons in a beam. |
---|
| 785 | |
---|
| 786 | void BeamParticle::list(ostream& os) const { |
---|
| 787 | |
---|
| 788 | // Header. |
---|
| 789 | os << "\n -------- PYTHIA Partons resolved in beam -----------------" |
---|
| 790 | << "-------------------------------------------------------------\n" |
---|
| 791 | << "\n i iPos id x comp xqcomp pTfact " |
---|
| 792 | << "colours p_x p_y p_z e m \n"; |
---|
| 793 | |
---|
| 794 | // Loop over list of removed partons and print it. |
---|
| 795 | double xSum = 0.; |
---|
| 796 | Vec4 pSum; |
---|
| 797 | for (int i = 0; i < size(); ++i) { |
---|
| 798 | ResolvedParton res = resolved[i]; |
---|
| 799 | os << fixed << setprecision(6) << setw(5) << i << setw(6) << res.iPos() |
---|
| 800 | << setw(8) << res.id() << setw(10) << res.x() << setw(6) |
---|
| 801 | << res.companion() << setw(10) << res.xqCompanion() << setw(10) |
---|
| 802 | << res.pTfactor() << setprecision(3) << setw(6) << res.col() |
---|
| 803 | << setw(6) << res.acol() << setw(11) << res.px() << setw(11) |
---|
| 804 | << res.py() << setw(11) << res.pz() << setw(11) << res.e() |
---|
| 805 | << setw(11) << res.m() << "\n"; |
---|
| 806 | |
---|
| 807 | // Also find sum of x and p values. |
---|
| 808 | if (res.companion() != -10) { |
---|
| 809 | xSum += res.x(); |
---|
| 810 | pSum += res.p(); |
---|
| 811 | } |
---|
| 812 | } |
---|
| 813 | |
---|
| 814 | // Print sum and endline. |
---|
| 815 | os << setprecision(6) << " x sum:" << setw(10) << xSum |
---|
| 816 | << setprecision(3) << " p sum:" |
---|
| 817 | << setw(11) << pSum.px() << setw(11) << pSum.py() << setw(11) |
---|
| 818 | << pSum.pz() << setw(11) << pSum.e() |
---|
| 819 | << "\n\n -------- End PYTHIA Partons resolved in beam -----------" |
---|
| 820 | << "---------------------------------------------------------------" |
---|
| 821 | << endl; |
---|
| 822 | } |
---|
| 823 | |
---|
| 824 | //-------------------------------------------------------------------------- |
---|
| 825 | |
---|
| 826 | // Test whether a lepton is to be considered as unresolved. |
---|
| 827 | |
---|
| 828 | bool BeamParticle::isUnresolvedLepton() { |
---|
| 829 | |
---|
| 830 | // Require record to consist of lepton with full energy plus a photon. |
---|
| 831 | if (!isLeptonBeam || resolved.size() > 2 || resolved[1].id() != 22 |
---|
| 832 | || resolved[0].x() < XMINUNRESOLVED) return false; |
---|
| 833 | return true; |
---|
| 834 | |
---|
| 835 | } |
---|
| 836 | |
---|
| 837 | //-------------------------------------------------------------------------- |
---|
| 838 | |
---|
| 839 | // For a diffractive system, decide whether to kick out gluon or quark. |
---|
| 840 | |
---|
| 841 | bool BeamParticle::pickGluon(double mDiff) { |
---|
| 842 | |
---|
| 843 | // Relative weight to pick a quark, assumed falling with energy. |
---|
| 844 | double probPickQuark = pickQuarkNorm / pow( mDiff, pickQuarkPower); |
---|
| 845 | return ( (1. + probPickQuark) * rndmPtr->flat() < 1. ); |
---|
| 846 | |
---|
| 847 | } |
---|
| 848 | |
---|
| 849 | //-------------------------------------------------------------------------- |
---|
| 850 | |
---|
| 851 | // Pick a valence quark at random. (Used for diffractive systems.) |
---|
| 852 | |
---|
| 853 | int BeamParticle::pickValence() { |
---|
| 854 | |
---|
| 855 | // Pick one valence quark at random. |
---|
| 856 | int nTotVal = (isBaryonBeam) ? 3 : 2; |
---|
| 857 | double rnVal = rndmPtr->flat() * nTotVal; |
---|
| 858 | int iVal = (rnVal < 1.) ? 1 : ( (rnVal < 2.) ? 2 : 3 ); |
---|
| 859 | |
---|
| 860 | // This valence in slot 1, the rest thereafter. |
---|
| 861 | idVal1 = 0; |
---|
| 862 | idVal2 = 0; |
---|
| 863 | idVal3 = 0; |
---|
| 864 | int iNow = 0; |
---|
| 865 | for (int i = 0; i < nValKinds; ++i) |
---|
| 866 | for (int j = 0; j < nVal[i]; ++j) { |
---|
| 867 | ++iNow; |
---|
| 868 | if (iNow == iVal) idVal1 = idVal[i]; |
---|
| 869 | else if ( idVal2 == 0) idVal2 = idVal[i]; |
---|
| 870 | else idVal3 = idVal[i]; |
---|
| 871 | } |
---|
| 872 | |
---|
| 873 | // Construct diquark if baryon. |
---|
| 874 | if (idVal3 != 0) idVal2 = flavSelPtr->makeDiquark( idVal2, idVal3); |
---|
| 875 | |
---|
| 876 | // Done. |
---|
| 877 | return idVal1; |
---|
| 878 | |
---|
| 879 | } |
---|
| 880 | |
---|
| 881 | //-------------------------------------------------------------------------- |
---|
| 882 | |
---|
| 883 | // Share lightcone momentum between two remnants in a diffractive system. |
---|
| 884 | |
---|
| 885 | double BeamParticle::zShare( double mDiff, double m1, double m2) { |
---|
| 886 | |
---|
| 887 | // Set up as valence in normal beam so can use xRemnant code. |
---|
| 888 | append(0, idVal1, 0., -3); |
---|
| 889 | append(0, idVal2, 0., -3); |
---|
| 890 | double m2Diff = mDiff*mDiff; |
---|
| 891 | |
---|
| 892 | // Begin to generate z and pT until acceptable solution. |
---|
| 893 | double wtAcc = 0.; |
---|
| 894 | do { |
---|
| 895 | double x1 = xRemnant(0); |
---|
| 896 | double x2 = xRemnant(0); |
---|
| 897 | zRel = x1 / (x1 + x2); |
---|
| 898 | pair<double, double> gauss2 = rndmPtr->gauss2(); |
---|
| 899 | pxRel = diffPrimKTwidth * gauss2.first; |
---|
| 900 | pyRel = diffPrimKTwidth * gauss2.second; |
---|
| 901 | |
---|
| 902 | // Suppress large invariant masses of remnant system. |
---|
| 903 | double mTS1 = m1*m1 + pxRel*pxRel + pyRel*pyRel; |
---|
| 904 | double mTS2 = m2*m2 + pxRel*pxRel + pyRel*pyRel; |
---|
| 905 | double m2Sys = mTS1 / zRel + mTS2 / (1. - zRel); |
---|
| 906 | wtAcc = (m2Sys < m2Diff) |
---|
| 907 | ? pow( 1. - m2Sys / m2Diff, diffLargeMassSuppress) : 0.; |
---|
| 908 | } while (wtAcc < rndmPtr->flat()); |
---|
| 909 | |
---|
| 910 | // Done. |
---|
| 911 | return zRel; |
---|
| 912 | |
---|
| 913 | } |
---|
| 914 | |
---|
| 915 | //========================================================================== |
---|
| 916 | |
---|
| 917 | } // end namespace Pythia8 |
---|