StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BeamParticle.cc
1 // BeamParticle.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 "Pythia8/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, idBeam);
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
Definition: AgUStep.h:26