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 |
---|