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) 2018 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 ColState class.
16 
17 //--------------------------------------------------------------------------
18 
19 // Small storage class needed to find the full colour state
20 // of all the scattered partons.
21 
22 class ColState {
23 public:
24  ColState() : nTotal(0.) {}
25  vector< pair<int,int> > lastSteps;
26  // nTotal is integer but can become extremely big.
27  double nTotal;
28 };
29 
30 //==========================================================================
31 
32 // The BeamParticle class.
33 
34 //--------------------------------------------------------------------------
35 
36 // Constants: could be changed here if desired, but normally should not.
37 // These are of technical nature, as described for each.
38 
39 // A lepton that takes (almost) the full beam energy does not leave a remnant.
40 const double BeamParticle::XMINUNRESOLVED = 1. - 1e-10;
41 
42 // Fictitious Pomeron mass to leave some room for beam remnant.
43 const double BeamParticle::POMERONMASS = 1.;
44 
45 // Avoid numerical instability in the x -> 1 limit for companion quark.
46 const double BeamParticle::XMAXCOMPANION = 0.99;
47 
48 // Avoid too extremely uneven momentum sharing.
49 const double BeamParticle::TINYZREL = 1e-8;
50 
51 // Maximum number of tries to find a suitable colour.
52 const int BeamParticle::NMAX = 1000;
53 
54 // Maximum number of tries to randomly assign colours to beam remnant.
55 // After this number is reached, a systematic approach is used.
56 const int BeamParticle::NRANDOMTRIES = 1000;
57 
58 //--------------------------------------------------------------------------
59 
60 // Initialize data on a beam particle and save pointers.
61 
62 void BeamParticle::init( int idIn, double pzIn, double eIn, double mIn,
63  Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtrIn,
64  Rndm* rndmPtrIn, PDF* pdfInPtr, PDF* pdfHardInPtr, bool isUnresolvedIn,
65  StringFlav* flavSelPtrIn) {
66 
67  // Store input pointers (and one bool) for future use.
68  infoPtr = infoPtrIn;
69  particleDataPtr = particleDataPtrIn;
70  rndmPtr = rndmPtrIn;
71  pdfBeamPtr = pdfInPtr;
72  pdfHardBeamPtr = pdfHardInPtr;
73  isUnresolvedBeam = isUnresolvedIn;
74  flavSelPtr = flavSelPtrIn;
75 
76  // Save the usual PDF pointers as the normal ones may be overwritten
77  // with unresolved PDFs when mixing different photoproduction modes.
78  pdfBeamPtrSave = pdfInPtr;
79  pdfHardBeamPtrSave = pdfHardInPtr;
80 
81  // Check whether beam has a supplementary photon beam.
82  bool beamHasGamma = settings.flag("PDF:lepton2gamma");
83 
84  // To be set process by process but start with this.
85  hasResGammaInBeam = beamHasGamma;
86 
87  // Maximum quark kind in allowed incoming beam hadrons.
88  maxValQuark = settings.mode("BeamRemnants:maxValQuark");
89 
90  // Power of (1-x)^power/sqrt(x) for remnant valence quark distribution.
91  valencePowerMeson = settings.parm("BeamRemnants:valencePowerMeson");
92  valencePowerUinP = settings.parm("BeamRemnants:valencePowerUinP");
93  valencePowerDinP = settings.parm("BeamRemnants:valencePowerDinP");
94 
95  // Enhancement factor of x of diquark.
96  valenceDiqEnhance = settings.parm("BeamRemnants:valenceDiqEnhance");
97 
98  // Assume g(x) ~ (1-x)^power/x to constrain companion to sea quark.
99  companionPower = settings.mode("BeamRemnants:companionPower");
100 
101  // Assume g(x) ~ (1-x)^power/x with a cut-off for low x.
102  gluonPower = settings.parm("BeamRemnants:gluonPower");
103  xGluonCutoff = settings.parm("BeamRemnants:xGluonCutoff");
104 
105  // Allow or not more than one valence quark to be kicked out.
106  allowJunction = settings.flag("BeamRemnants:allowJunction");
107 
108  // Choose whether to form a di-quark or
109  // a junction with new colur reconnection scheme.
110  beamJunction = settings.flag("beamRemnants:beamJunction");
111 
112  // Allow junctions in the outgoing colour state.
113  allowBeamJunctions = settings.flag("beamRemnants:allowBeamJunction");
114 
115  // For low-mass diffractive system kick out q/g = norm / mass^power.
116  pickQuarkNorm = settings.parm("Diffraction:pickQuarkNorm");
117  pickQuarkPower = settings.parm("Diffraction:pickQuarkPower");
118 
119  // Controls the amount of saturation in the new model.
120  beamSat = settings.parm("BeamRemnants:saturation");
121 
122  // Width of primordial kT distribution in low-mass diffractive systems.
123  diffPrimKTwidth = settings.parm("Diffraction:primKTwidth");
124 
125  // Suppress large masses of beam remnant in low-mass diffractive systems.
126  diffLargeMassSuppress = settings.parm("Diffraction:largeMassSuppress");
127 
128  // Check if ISR for photon collisions is applied and set pTmin.
129  doND = settings.flag("SoftQCD:nonDiffractive");
130  doISR = settings.flag("PartonLevel:ISR");
131  doMPI = settings.flag("PartonLevel:MPI");
132  pTminISR = settings.parm("SpaceShower:pTmin");
133 
134  // Store info on the incoming beam.
135  idBeam = idIn;
136  initBeamKind();
137  pBeam = Vec4( 0., 0., pzIn, eIn);
138  mBeam = mIn;
139 
140  // Initialize parameters related to photon beams.
141  resetGamma();
142  resetGammaInLepton();
143 
144  clear();
145 
146 }
147 
148 //--------------------------------------------------------------------------
149 
150 // Initialize kind and valence flavour content of incoming beam.
151 // For recognized hadrons one can generate multiparton interactions.
152 // Dynamic choice of meson valence flavours in newValenceContent below.
153 
154 void BeamParticle::initBeamKind() {
155 
156  // Reset.
157  idBeamAbs = abs(idBeam);
158  isLeptonBeam = false;
159  isHadronBeam = false;
160  isMesonBeam = false;
161  isBaryonBeam = false;
162  isGammaBeam = false;
163  nValKinds = 0;
164 
165  // To be modified according to the process.
166  gammaMode = 0;
167  isResUnres = false;
168 
169  // Check for leptons or DM beams.
170  if ( (idBeamAbs > 10 && idBeamAbs < 17)
171  || (idBeamAbs > 50 && idBeamAbs < 60) ) {
172  nValKinds = 1;
173  nVal[0] = 1;
174  idVal[0] = idBeam;
175  isLeptonBeam = true;
176  }
177 
178  // Valence content for photons.
179  if (idBeamAbs == 22) {
180  isGammaBeam = true;
181  nValKinds = 2;
182  nVal[0] = 1;
183  nVal[1] = 1;
184  newValenceContent();
185  iPosVal = -1;
186  }
187 
188  // Done if cannot be lowest-lying hadron state.
189  if (idBeamAbs < 101 || idBeamAbs > 9999) return;
190 
191  // Resolve valence content for assumed Pomeron.
192  if (idBeamAbs == 990) {
193  isMesonBeam = true;
194  nValKinds = 2;
195  nVal[0] = 1 ;
196  nVal[1] = 1 ;
197  newValenceContent();
198 
199  // Resolve valence content for assumed meson. Flunk unallowed codes.
200  } else if (idBeamAbs < 1000) {
201  int id1 = idBeamAbs/100;
202  int id2 = (idBeamAbs/10)%10;
203  if ( id1 < 1 || id1 > maxValQuark
204  || id2 < 1 || id2 > maxValQuark ) return;
205  isMesonBeam = true;
206 
207  // Store valence content of a confirmed meson.
208  nValKinds = 2;
209  nVal[0] = 1 ;
210  nVal[1] = 1;
211  if (id1%2 == 0) {
212  idVal[0] = id1;
213  idVal[1] = -id2;
214  } else {
215  idVal[0] = id2;
216  idVal[1] = -id1;
217  }
218  newValenceContent();
219 
220  // Resolve valence content for assumed baryon. Flunk unallowed codes.
221  } else {
222  int id1 = idBeamAbs/1000;
223  int id2 = (idBeamAbs/100)%10;
224  int id3 = (idBeamAbs/10)%10;
225  if ( id1 < 1 || id1 > maxValQuark || id2 < 1 || id2 > maxValQuark
226  || id3 < 1 || id3 > maxValQuark) return;
227  if (id2 > id1 || id3 > id1) return;
228  isBaryonBeam = true;
229 
230  // Store valence content of a confirmed baryon.
231  nValKinds = 1; idVal[0] = id1; nVal[0] = 1;
232  if (id2 == id1) ++nVal[0];
233  else {
234  nValKinds = 2;
235  idVal[1] = id2;
236  nVal[1] = 1;
237  }
238  if (id3 == id1) ++nVal[0];
239  else if (id3 == id2) ++nVal[1];
240  else {
241  idVal[nValKinds] = id3;
242  nVal[nValKinds] = 1;
243  ++nValKinds;
244  }
245  }
246 
247  // Flip flavours for antimeson or antibaryon, and then done.
248  if (idBeam < 0) for (int i = 0; i < nValKinds; ++i) idVal[i] = -idVal[i];
249  isHadronBeam = true;
250  Q2ValFracSav = -1.;
251 
252 }
253 
254 //--------------------------------------------------------------------------
255 
256 // Initialize the photon beam with additional unresolved PDF pointer.
257 
258 void BeamParticle::initUnres(PDF* pdfUnresInPtr) {
259 
260  // Set the pointer and check that pointer exists.
261  pdfUnresBeamPtr = pdfUnresInPtr;
262  isResUnres = (pdfUnresBeamPtr != 0 ) ? true : false;
263 
264 }
265 
266 //--------------------------------------------------------------------------
267 
268 // Dynamic choice of meson valence flavours for pi0, K0S, K0L, Pomeron.
269 
270 void BeamParticle::newValenceContent() {
271 
272  // A pi0, rho and omega oscillates between d dbar and u ubar.
273  if (idBeam == 111 || idBeam == 113 || idBeam == 223) {
274  idVal[0] = (rndmPtr->flat() < 0.5) ? 1 : 2;
275  idVal[1] = -idVal[0];
276 
277  // A K0S or K0L oscillates between d sbar and s dbar.
278  } else if (idBeam == 130 || idBeam == 310) {
279  idVal[0] = (rndmPtr->flat() < 0.5) ? 1 : 3;
280  idVal[1] = (idVal[0] == 1) ? -3 : -1;
281 
282  // For a Pomeron split gluon remnant into d dbar or u ubar.
283  } else if (idBeam == 990) {
284  idVal[0] = (rndmPtr->flat() < 0.5) ? 1 : 2;
285  idVal[1] = -idVal[0];
286 
287  // For photons a VMD content may be chosen, or the choice may be delayed.
288  } else if (idBeam == 22) {
289  if (hasVMDstate()) {
290  int idTmp = idVMD();
291  // A rho and omega oscillates between a uubar and ddbar.
292  if (idTmp == 113 || idTmp == 223) {
293  idVal[0] = (rndmPtr->flat() < 0.5) ? 1 : 2;
294  idVal[1] = -idVal[0];
295  // A phi is an ssbar system.
296  } else if (idTmp == 333) {
297  idVal[0] = 3;
298  idVal[1] = -3;
299  } else return;
300  // For non-VMD photons set an usused code to indicate that the flavour
301  // is not chosen yet but will be done later.
302  } else {
303  idVal[0] = 10;
304  idVal[1] = -idVal[0];
305  }
306 
307  // If phi meson set content to s sbar.
308  } else if (idBeam == 333) {
309  idVal[0] = 3;
310  idVal[1] = -idVal[0];
311 
312  // Other hadrons so far do not require any event-by-event change.
313  } else return;
314 
315  // Propagate change to PDF routine(s).
316  pdfBeamPtr->newValenceContent( idVal[0], idVal[1]);
317  if (pdfHardBeamPtr != pdfBeamPtr && pdfHardBeamPtr != 0)
318  pdfHardBeamPtr->newValenceContent( idVal[0], idVal[1]);
319 
320 }
321 
322 //--------------------------------------------------------------------------
323 
324 double BeamParticle::xMax(int iSkip) {
325 
326  // Minimum requirement on remaining energy > nominal mass for hadron.
327  double xLeft = 1.;
328  if (idBeam == 990) xLeft -= POMERONMASS / e();
329  else if (isHadron()) xLeft -= m() / e();
330  if (size() == 0) return xLeft;
331 
332  // Subtract what was carried away by initiators (to date).
333  for (int i = 0; i < size(); ++i)
334  if (i != iSkip && resolved[i].isFromBeam()) xLeft -= resolved[i].x();
335  return xLeft;
336 
337 }
338 
339 //--------------------------------------------------------------------------
340 
341 // Parton distributions, reshaped to take into account previous
342 // multiparton interactions. By picking a non-negative iSkip value,
343 // one particular interaction is skipped, as needed for ISR
344 
345 double BeamParticle::xfModified(int iSkip, int idIn, double x, double Q2) {
346 
347  // Initial values.
348  idSave = idIn;
349  iSkipSave = iSkip;
350  xqVal = 0.;
351  xqgSea = 0.;
352  xqCompSum = 0.;
353 
354  // Fast procedure for first interaction.
355  if (size() == 0) {
356  if (x >= 1.) return 0.;
357  bool canBeVal = false;
358  for (int i = 0; i < nValKinds; ++i)
359  if (idIn == idVal[i]) canBeVal = true;
360  if (canBeVal) {
361  xqVal = xfVal( idIn, x, Q2);
362  xqgSea = xfSea( idIn, x, Q2);
363  }
364  else xqgSea = xf( idIn, x, Q2);
365 
366  // More complicated procedure for non-first interaction.
367  } else {
368 
369  // Sum up the x already removed, and check that remaining x is enough.
370  double xUsed = 0.;
371  for (int i = 0; i < size(); ++i)
372  if (i != iSkip) xUsed += resolved[i].x();
373  double xLeft = 1. - xUsed;
374  if (x >= xLeft) return 0.;
375  double xRescaled = x / xLeft;
376 
377  // Calculate total and remaining amount of x carried by valence quarks.
378  double xValTot = 0.;
379  double xValLeft = 0.;
380  for (int i = 0; i < nValKinds; ++i) {
381  nValLeft[i] = nVal[i];
382  for (int j = 0; j < size(); ++j)
383  if (j != iSkip && resolved[j].isValence()
384  && resolved[j].id() == idVal[i]) --nValLeft[i];
385  double xValNow = xValFrac(i, Q2);
386  xValTot += nVal[i] * xValNow;
387  xValLeft += nValLeft[i] * xValNow;
388  }
389 
390  // Calculate total amount of x carried by unmatched companion quarks.
391  double xCompAdded = 0.;
392  for (int i = 0; i < size(); ++i)
393  if (i != iSkip && resolved[i].isUnmatched()) xCompAdded
394  += xCompFrac( resolved[i].x() / (xLeft + resolved[i].x()) )
395  // Typo warning: extrafactor missing in Skands&Sjostrand article;
396  // <x> for companion refers to fraction of x left INCLUDING sea quark.
397  // To be modified further??
398  * (1. + resolved[i].x() / xLeft);
399 
400  // Calculate total rescaling factor and pdf for sea and gluon.
401  double rescaleGS = max( 0., (1. - xValLeft - xCompAdded)
402  / (1. - xValTot) );
403  xqgSea = rescaleGS * xfSea( idIn, xRescaled, Q2);
404 
405  // Find valence part and rescale it to remaining number of quarks.
406  for (int i = 0; i < nValKinds; ++i)
407  if (idIn == idVal[i] && nValLeft[i] > 0)
408  xqVal = xfVal( idIn, xRescaled, Q2)
409  * double(nValLeft[i]) / double(nVal[i]);
410 
411  // Find companion part, by adding all companion contributions.
412  for (int i = 0; i < size(); ++i)
413  if (i != iSkip && resolved[i].id() == -idIn
414  && resolved[i].isUnmatched()) {
415  double xsRescaled = resolved[i].x() / (xLeft + resolved[i].x());
416  double xcRescaled = x / (xLeft + resolved[i].x());
417  double xqCompNow = xCompDist( xcRescaled, xsRescaled);
418  // Normalize the companion quark PDF to the total momentum carried
419  // by the partons in case of photon beam at given scale Q^2.
420  if ( isGamma() ) xqCompNow *= xIntegratedPDFs(Q2);
421  resolved[i].xqCompanion( xqCompNow);
422  xqCompSum += xqCompNow;
423  }
424  }
425 
426  // Add total, but only return relevant part for ISR. More cases??
427  // Watch out, e.g. g can come from either kind of quark.??
428  xqgTot = xqVal + xqgSea + xqCompSum;
429 
430  // If ISR with photon beams no distinction between valence and sea.
431  if (isGammaBeam && doISR) return xqgTot;
432 
433  if (iSkip >= 0) {
434  if (resolved[iSkip].isValence()) return xqVal;
435  if (resolved[iSkip].isUnmatched()) return xqgSea + xqCompSum;
436  }
437  return xqgTot;
438 
439 }
440 
441 //--------------------------------------------------------------------------
442 
443 // Decide whether a quark extracted from the beam is of valence, sea or
444 // companion kind; in the latter case also pick its companion.
445 // Assumes xfModified has already been called.
446 
447 int BeamParticle::pickValSeaComp() {
448 
449  // If parton already has a companion than reset code for this.
450  int oldCompanion = resolved[iSkipSave].companion();
451  if (oldCompanion >= 0) resolved[oldCompanion].companion(-2);
452 
453  // Default assignment is sea.
454  int vsc = -2;
455 
456  // For gluons or photons no sense of valence or sea.
457  if (idSave == 21 || idSave == 22) vsc = -1;
458 
459  // For lepton beam assume same-kind lepton inside is valence.
460  else if (isLeptonBeam && idSave == idBeam) vsc = -3;
461 
462  // Decide if valence or sea quark.
463  // For photons, consider all partons as sea until valence content fixed.
464  else {
465  double xqRndm = xqgTot * rndmPtr->flat();
466  if (xqRndm < xqVal && !isGammaBeam ) vsc = -3;
467  else if (xqRndm < xqVal + xqgSea) vsc = -2;
468 
469  // If not either, loop over all possible companion quarks.
470  else {
471  xqRndm -= xqVal + xqgSea;
472  for (int i = 0; i < size(); ++i)
473  if (i != iSkipSave && resolved[i].id() == -idSave
474  && resolved[i].isUnmatched()) {
475  xqRndm -= resolved[i].xqCompanion();
476  if (xqRndm < 0.) vsc = i;
477  break;
478  }
479  }
480  }
481 
482  // Bookkeep assignment; for sea--companion pair both ways.
483  resolved[iSkipSave].companion(vsc);
484  if (vsc >= 0) resolved[vsc].companion(iSkipSave);
485 
486  // Done; return code for choice (to distinguish valence/sea in Info).
487  return vsc;
488 
489 }
490 
491 //--------------------------------------------------------------------------
492 
493 // Fraction of hadron momentum sitting in a valence quark distribution.
494 // Based on hardcoded parametrizations of CTEQ 5L numbers.
495 
496 double BeamParticle::xValFrac(int j, double Q2) {
497 
498  // Only recalculate when required.
499  if (Q2 != Q2ValFracSav) {
500  Q2ValFracSav = Q2;
501 
502  // Q2-dependence of log-log form; assume fixed Lambda = 0.2.
503  double llQ2 = log( log( max( 1., Q2) / 0.04 ));
504 
505  // Fractions carried by u and d in proton.
506  uValInt = 0.48 / (1. + 1.56 * llQ2);
507  dValInt = 0.385 / (1. + 1.60 * llQ2);
508  }
509 
510  // Baryon with three different quark kinds: (2 * u + d) / 3 of proton.
511  if (isBaryonBeam && nValKinds == 3) return (2. * uValInt + dValInt) / 3.;
512 
513  // Baryon with one or two identical: like d or u of proton.
514  if (isBaryonBeam && nVal[j] == 1) return dValInt;
515  if (isBaryonBeam && nVal[j] == 2) return uValInt;
516 
517  // Meson: (2 * u + d) / 2 of proton so same total valence quark fraction.
518  return 0.5 * (2. * uValInt + dValInt);
519 
520 }
521 
522 //--------------------------------------------------------------------------
523 
524 // The momentum integral of a companion quark, with its partner at x_s,
525 // using an approximate gluon density like (1 - x_g)^power / x_g.
526 // The value corresponds to an unrescaled range between 0 and 1 - x_s.
527 
528 double BeamParticle::xCompFrac(double xs) {
529 
530  // Tiny answer for xs -> 1 is numerically unstable, so set = 0.
531  if (xs > XMAXCOMPANION) return 0.;
532 
533  // Select case by power of gluon (1-x_g) shape.
534  switch (companionPower) {
535 
536  case 0:
537  return xs * ( 5. + xs * (-9. - 2. * xs * (-3. + xs)) + 3. * log(xs) )
538  / ( (-1. + xs) * (2. + xs * (-1. + 2. * xs)) );
539 
540  case 1:
541  return -1. -3. * xs + ( 2. * pow2(-1. + xs) * (1. + xs + xs*xs))
542  / ( 2. + xs*xs * (xs - 3.) + 3. * xs * log(xs) );
543 
544  case 2:
545  return xs * ( (1. - xs) * (19. + xs * (43. + 4. * xs))
546  + 6. * log(xs) * (1. + 6. * xs + 4.*xs*xs) ) /
547  ( 4. * ( (xs - 1.) * (1. + xs * (4. + xs) )
548  - 3. * xs * log(xs) * (1 + xs) ) );
549 
550  case 3:
551  return 3. * xs * ( (xs - 1.) * (7. + xs * (28. + 13. * xs))
552  - 2. * log(xs) * (1. + xs * (9. + 2. * xs * (6. + xs))) )
553  / ( 4. + 27. * xs - 31. * pow3(xs)
554  + 6. * xs * log(xs) * (3. + 2. * xs * (3.+xs)) );
555 
556  default:
557  return ( -9. * xs * (xs*xs - 1.) * (5. + xs * (24. + xs)) + 12. * xs
558  * log(xs) * (1. + 2. * xs) * (1. + 2. * xs * (5. + 2. * xs)) )
559  / ( 8. * (1. + 2. * xs) * ((xs - 1.) * (1. + xs * (10. + xs))
560  - 6. * xs * log(xs) * (1. + xs)) );
561 
562  }
563 }
564 
565 //--------------------------------------------------------------------------
566 
567 // The x*f pdf of a companion quark at x_c, with its sea partner at x_s,
568 // using an approximate gluon density like (1 - x_g)^power / x_g.
569 // The value corresponds to an unrescaled range between 0 and 1 - x_s.
570 
571 double BeamParticle::xCompDist(double xc, double xs) {
572 
573  // Tiny answer for xs -> 1 is numerically unstable, so set = 0.
574  if (xs > XMAXCOMPANION) return 0.;
575 
576  // Mother gluon momentum fraction. Check physical limit.
577  double xg = xc + xs;
578  if (xg > 1.) return 0.;
579 
580  // Common factor, including splitting kernel and part of gluon density
581  // (and that it is x_c * f that is coded).
582  double fac = 3. * xc * xs * (xc*xc + xs*xs) / pow4(xg);
583 
584  // Select case by power of gluon (1-x_g) shape.
585  switch (companionPower) {
586 
587  case 0:
588  return fac / ( 2. - xs * (3. - xs * (3. - 2. * xs)) );
589 
590  case 1:
591  return fac * (1. - xg) / ( 2. + xs*xs * (-3. + xs) + 3. * xs * log(xs) );
592 
593  case 2:
594  return fac * pow2(1. - xg) / ( 2. * ((1. - xs) * (1. + xs * (4. + xs))
595  + 3. * xs * (1. + xs) * log(xs)) );
596 
597  case 3:
598  return fac * pow3(1. - xg) * 2. / ( 4. + 27. * xs - 31. * pow3(xs)
599  + 6. * xs * log(xs) * (3. + 2. * xs * (3. + xs)) );
600 
601  default:
602  return fac * pow4(1. - xg) / ( 2. * (1. + 2. * xs) * ((1. - xs)
603  * (1. + xs * (10. + xs)) + 6. * xs * log(xs) * (1. + xs)) );
604 
605  }
606 }
607 
608 //--------------------------------------------------------------------------
609 
610 // Set the photon mode (none (0), resolved (1), unresolved (2)) of the beam.
611 
612 void BeamParticle::setGammaMode(int gammaModeIn) {
613 
614  // For hadrons mode always 0.
615  if (isHadron()) {
616  gammaMode = 0;
617  pdfBeamPtr = pdfBeamPtrSave;
618  pdfHardBeamPtr = pdfHardBeamPtrSave;
619  hasResGammaInBeam = false;
620  isResolvedGamma = false;
621  return;
622  }
623 
624  // Save the mode of the photon beam.
625  gammaMode = gammaModeIn;
626 
627  // Set the beam and PDF pointers to unresolved mode.
628  if (gammaMode == 2 && isResUnres) {
629  pdfBeamPtr = pdfUnresBeamPtr;
630  pdfHardBeamPtr = pdfUnresBeamPtr;
631  isResolvedGamma = false;
632  hasResGammaInBeam = false;
633 
634  // Only a photon beam can be unresolved with gammaMode == 2.
635  if ( isGamma()) isUnresolvedBeam = true;
636 
637  // Set the beam and PDF pointers to resolved mode.
638  } else {
639  pdfBeamPtr = pdfBeamPtrSave;
640  pdfHardBeamPtr = pdfHardBeamPtrSave;
641  isUnresolvedBeam = false;
642  if ( isGamma()) isResolvedGamma = true;
643  else isResolvedGamma = false;
644  if ( isLepton() && gammaMode == 1 ) hasResGammaInBeam = true;
645  else hasResGammaInBeam = false;
646  }
647 
648 }
649 
650 //--------------------------------------------------------------------------
651 
652 // Check whether parton iResolved with given Q^2 is a valence quark.
653 
654 bool BeamParticle::gammaInitiatorIsVal(int iResolved, double Q2) {
655  return gammaInitiatorIsVal( iResolved, resolved[iResolved].id(),
656  resolved[iResolved].x(), Q2 );
657 }
658 
659 //--------------------------------------------------------------------------
660 
661 // Check whether initiator parton is a valence quark using the PDFs.
662 // Set the position of the valence quark to iGamVal.
663 
664 bool BeamParticle::gammaInitiatorIsVal(int iResolved, int idInit,
665  double x, double Q2) {
666 
667  // Reset the valence quark position.
668  iPosVal = -1;
669 
670  // Gluon is not a valence parton. Sample content accordingly.
671  if ( idInit == 0 || abs(idInit) == 21 ) {
672  idVal[0] = pdfBeamPtr->sampleGammaValFlavor(Q2);
673  idVal[1] = -idVal[0];
674  return false;
675  } else {
676 
677  // Set the valence content to match with the hard process to get the
678  // correct PDFs and to store the choice. Changed by sampleGammaValFlavor.
679  idVal[0] = idInit;
680  idVal[1] = -idInit;
681  pdfBeamPtr->newValenceContent( idVal[0], idVal[1]);
682 
683  // If initiator from gamma->qqbar splitting then it is a valence quark.
684  if ( iResolved == iGamVal ) {
685  iPosVal = iResolved;
686  return true;
687 
688  // If Q^2 is smaller than mass of quark set to valence.
689  } else if ( Q2 < pdfBeamPtr->gammaPDFRefScale(idInit) ) {
690  iPosVal = iResolved;
691  return true;
692 
693  // Use PDFs to decide if valence parton.
694  } else {
695  double xVal = xfVal( idInit, x, Q2);
696  double xSea = xfSea( idInit, x, Q2);
697  if ( rndmPtr->flat() < xVal/( xVal + xSea ) ) {
698  iPosVal = iResolved;
699  return true;
700 
701  // If the initiator not valence sample the flavour.
702  } else {
703  idVal[0] = pdfBeamPtr->sampleGammaValFlavor(Q2);
704  idVal[1] = -idVal[0];
705  return false;
706  }
707  }
708  }
709 }
710 
711 //--------------------------------------------------------------------------
712 
713 // Return the type of the hard parton from a photon beam.
714 
715 int BeamParticle::gammaValSeaComp(int iResolved) {
716 
717  // Default choice a sea quark.
718  int vsc = -2;
719 
720  // Gluons and photons -1.
721  if ( resolved[iResolved].id() == 21 ||
722  resolved[iResolved].id() == 22 ) vsc = -1;
723 
724  // Quarks are valence partons if decided so earlier.
725  else if (iResolved == iPosVal) vsc = -3;
726  resolved[iResolved].companion(vsc);
727 
728  return vsc;
729 }
730 
731 //--------------------------------------------------------------------------
732 
733 // Add required extra remnant flavour content. Also initial colours.
734 
735 bool BeamParticle::remnantFlavours(Event& event, bool isDIS) {
736 
737  // Elastically scattered beam, e.g. for coherent photon from proton.
738  if (isHadronBeam && isUnresolvedBeam) {
739  append( 0, idBeam, 0., -1);
740  resolved[1].m( particleDataPtr->m0( idBeamAbs ) );
741  return true;
742  }
743 
744  // A baryon will have a junction, unless a diquark is formed later.
745  hasJunctionBeam = (isBaryon());
746 
747  // Store how many hard-scattering partons were removed from beam.
748  nInit = size();
749  if (isDIS && nInit != 1) return false;
750 
751  // Decide the valence content of photon beam here if ISR is applied.
752  if ( isGammaBeam ) {
753 
754  // For direct-resolved processes no need for remnants on direct side.
755  if (resolved[0].id() == 22) return true;
756 
757  // If ISR but no MPIs check only for the one initiator.
758  if ( doISR && !doMPI ) {
759 
760  // If remnants are constructed fix the valence content using the
761  // scale where ISR have stopped.
762  if ( isResolvedGamma ) {
763  gammaInitiatorIsVal(0, pow2(pTminISR));
764  }
765 
766  // Set the initiator companion code after the valence content is fixed.
767  gammaValSeaComp(0);
768 
769  } else if ( doMPI || doND ) {
770 
771  // If ISR is applied, use the min. scale of evolution for the valence
772  // parton decision. Otherwise use the pT of latest MPI (set in scatter).
773  double pTmin = doISR ? pTminISR : pTminMPI;
774 
775  // If a quark from gamma->qqbar exists, this must be the valence so
776  // set it to valence and the possible companion to other valence.
777  if ( iGamVal >= 0 ) {
778  gammaInitiatorIsVal(iGamVal, pow2(pTmin));
779  int valComp = resolved[iGamVal].companion();
780  if ( valComp >= 0 ) resolved[valComp].companion(-3);
781  gammaValSeaComp(iGamVal);
782  } else {
783 
784  // Loop through initiators starting from parton from the hardest
785  // interaction and check whether valence. When found, no need to
786  // check further.
787  for ( int i = 0; i < size(); ++i) {
788  bool isValence = gammaInitiatorIsVal(i, pow2(pTminISR));
789  int origComp = resolved[i].companion();
790  gammaValSeaComp(i);
791  if (isValence) {
792  // If the chosen valence parton has a companion, set this to be
793  // the other valence parton.
794  if ( origComp >= 0 ) resolved[origComp].companion(-3);
795  break;
796  }
797  }
798  }
799  }
800  }
801 
802  // No need for remnants with unresolved photon from leptons.
803  else if ( isLeptonBeam && gammaMode == 2) return true;
804 
805  // Find remaining valence quarks.
806  for (int i = 0; i < nValKinds; ++i) {
807  nValLeft[i] = nVal[i];
808  for (int j = 0; j < nInit; ++j) if (resolved[j].isValence()
809  && resolved[j].id() == idVal[i]) --nValLeft[i];
810  // No valence quarks if ISR find the original beam photon.
811  if ( isGammaBeam && doISR && !isResolvedGamma && !doMPI ) nValLeft[i] = 0;
812  // Add remaining valence quarks to record. Partly temporary values.
813  for (int k = 0; k < nValLeft[i]; ++k) append(0, idVal[i], 0., -3);
814  }
815 
816  // If at least two valence quarks left in baryon remnant then form diquark.
817  int nInitPlusVal = size();
818  if (isBaryon() && nInitPlusVal - nInit >= 2) {
819 
820  // If three, pick two at random to form diquark, else trivial.
821  int iQ1 = nInit;
822  int iQ2 = nInit + 1;
823  if (nInitPlusVal - nInit == 3) {
824  double pickDq = 3. * rndmPtr->flat();
825  if (pickDq > 1.) iQ2 = nInit + 2;
826  if (pickDq > 2.) iQ1 = nInit + 1;
827  }
828 
829  // Pick spin 0 or 1 according to SU(6) wave function factors.
830  int idDq = flavSelPtr->makeDiquark( resolved[iQ1].id(),
831  resolved[iQ2].id(), idBeam);
832 
833  // Overwrite with diquark flavour and remove one slot. No more junction.
834  resolved[iQ1].id(idDq);
835  if (nInitPlusVal - nInit == 3 && iQ2 == nInit + 1)
836  resolved[nInit + 1].id( resolved[nInit + 2].id() );
837  resolved.pop_back();
838  hasJunctionBeam = false;
839  }
840 
841  // Find companion quarks to unmatched sea quarks.
842  for (int i = 0; i < nInit; ++i)
843  if (resolved[i].isUnmatched()) {
844 
845  // Add companion quark to record; and bookkeep both ways.
846  append(0, -resolved[i].id(), 0., i);
847  resolved[i].companion(size() - 1);
848  }
849 
850  // If no other remnants found, add a gluon or photon to carry momentum.
851  // Add partons for photons only if remnants needed.
852  if ( size() == nInit && !isUnresolvedBeam &&
853  (!isGammaBeam || isResolvedGamma || doMPI) ) {
854  int idRemnant = (isHadronBeam || isGammaBeam) ? 21 : 22;
855  append(0, idRemnant, 0., -1);
856  }
857 
858  // For DIS allow collapse to one colour singlet hadron.
859  if (isHadronBeam && isDIS && size() > 2 && resolved[0].id() != 21) {
860  if (size() != 4) {
861  infoPtr->errorMsg("Error in BeamParticle::remnantFlavours: "
862  "unexpected number of beam remnants for DIS");
863  return false;
864  }
865 
866  // Companion last; find parton with matching colour.
867  int colTypeComp = particleDataPtr->colType( resolved[3].id() );
868  int colType1 = particleDataPtr->colType( resolved[1].id() );
869  int i12 = (colType1 == -colTypeComp) ? 1 : 2;
870 
871  // Combine to new hadron flavour.
872  int idHad = flavSelPtr->combineId( resolved[i12].id(), resolved[3].id() );
873  if (idHad == 0) {
874  infoPtr->errorMsg("Error in BeamParticle::remnantFlavours: "
875  "failed to combine hadron for DIS");
876  return false;
877  }
878 
879  // Overwrite with hadron flavour and remove companion.
880  resolved[i12].id(idHad);
881  resolved.pop_back();
882  resolved[0].companion(-3);
883  }
884 
885  // Set initiator and remnant masses.
886  for (int i = 0; i < size(); ++i) {
887  if (i < nInit) resolved[i].m(0.);
888  else resolved[i].m( particleDataPtr->m0( resolved[i].id() ) );
889  }
890 
891  // For debug purposes: reject beams with resolved junction topology.
892  if (hasJunctionBeam && !allowJunction) return false;
893 
894  // Pick initial colours for remnants.
895  for (int i = nInit; i < size(); ++i) {
896  int colType = particleDataPtr->colType( resolved[i].id() );
897  int col = (colType == 1 || colType == 2) ? event.nextColTag() : 0;
898  int acol = (colType == -1 || colType == 2) ? event.nextColTag() : 0;
899  resolved[i].cols( col, acol);
900  }
901 
902  // Done.
903  return true;
904 
905 }
906 
907 //--------------------------------------------------------------------------
908 
909 // Correlate all initiators and remnants to make a colour singlet.
910 
911 bool BeamParticle::remnantColours(Event& event, vector<int>& colFrom,
912  vector<int>& colTo) {
913 
914  // No colours in lepton beams so no need to do anything.
915  if (isLeptonBeam) return true;
916 
917  // Copy initiator colour info from the event record to the beam.
918  for (int i = 0; i < size(); ++i) {
919  int j = resolved[i].iPos();
920  resolved[i].cols( event[j].col(), event[j].acol());
921  }
922 
923  // Find number and position of valence quarks, of gluons, and
924  // of sea-companion pairs (counted as gluons) in the beam remnants.
925  // Skip gluons with same colour as anticolour and rescattering partons.
926  vector<int> iVal;
927  vector<int> iGlu;
928  for (int i = 0; i < size(); ++i)
929  if (resolved[i].isFromBeam()) {
930  if ( resolved[i].isValence() ) iVal.push_back(i);
931  else if ( resolved[i].isCompanion() && resolved[i].companion() > i )
932  iGlu.push_back(i);
933  else if ( resolved[i].id() == 21
934  && resolved[i].col() != resolved[i].acol() ) iGlu.push_back(i);
935  }
936 
937  // Pick a valence quark to which gluons are attached.
938  // Do not resolve quarks in diquark. (More sophisticated??)
939  int iValSel = 0;
940  if (iVal.size() != 0) {
941  iValSel = iVal[0];
942  if (iVal.size() == 2) {
943  if ( abs(resolved[iValSel].id()) > 10 ) iValSel = iVal[1];
944  } else {
945  double rndmValSel = 3. * rndmPtr->flat();
946  if (rndmValSel > 1.) iValSel= iVal[1];
947  if (rndmValSel > 2.) iValSel= iVal[2];
948  }
949  }
950 
951  // This valence quark defines initial (anti)colour.
952  int iBeg = iValSel;
953  bool hasCol = (resolved[iBeg].col() > 0);
954  int begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol();
955 
956  // Do random stepping through gluon/(sea+companion) list.
957  vector<int> iGluRndm;
958  for (int i = 0; i < int(iGlu.size()); ++i)
959  iGluRndm.push_back( iGlu[i] );
960  for (int iOrder = 0; iOrder < int(iGlu.size()); ++iOrder) {
961  int iRndm = int( double(iGluRndm.size()) * rndmPtr->flat());
962  int iGluSel = iGluRndm[iRndm];
963  iGluRndm[iRndm] = iGluRndm[iGluRndm.size() - 1];
964  iGluRndm.pop_back();
965 
966  // Find matching anticolour/colour to current colour/anticolour.
967  int iEnd = iGluSel;
968  int endCol = (hasCol) ? resolved[iEnd].acol() : resolved[iEnd].col();
969  // Not gluon but sea+companion pair: go to other.
970  if (endCol == 0) {
971  iEnd = resolved[iEnd].companion();
972  endCol = (hasCol) ? resolved[iEnd].acol() : resolved[iEnd].col();
973  }
974 
975  // Collapse this colour-anticolour pair to the lowest one.
976  if (begCol < endCol) {
977  if (hasCol) resolved[iEnd].acol(begCol);
978  else resolved[iEnd].col(begCol);
979  colFrom.push_back(endCol);
980  colTo.push_back(begCol);
981  } else {
982  if (hasCol) resolved[iBeg].col(endCol);
983  else resolved[iBeg].acol(endCol);
984  colFrom.push_back(begCol);
985  colTo.push_back(endCol);
986  }
987 
988  // Pick up the other colour of the recent gluon and repeat.
989  iBeg = iEnd;
990  begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol();
991  // Not gluon but sea+companion pair: go to other.
992  if (begCol == 0) {
993  iBeg = resolved[iBeg].companion();
994  begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol();
995  }
996 
997  // At end of gluon/(sea+companion) list.
998  }
999 
1000  // Now begin checks, and also finding junction information.
1001  // Loop through remnant partons; isolate all colours and anticolours.
1002  vector<int> colList;
1003  vector<int> acolList;
1004  for (int i = 0; i < size(); ++i)
1005  if ( resolved[i].isFromBeam() )
1006  if ( resolved[i].col() != resolved[i].acol() ) {
1007  if (resolved[i].col() > 0) colList.push_back( resolved[i].col() );
1008  if (resolved[i].acol() > 0) acolList.push_back( resolved[i].acol() );
1009  }
1010 
1011  // Remove all matching colour-anticolour pairs.
1012  bool foundPair = true;
1013  while (foundPair && colList.size() > 0 && acolList.size() > 0) {
1014  foundPair = false;
1015  for (int iCol = 0; iCol < int(colList.size()); ++iCol) {
1016  for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) {
1017  if (acolList[iAcol] == colList[iCol]) {
1018  colList[iCol] = colList.back();
1019  colList.pop_back();
1020  acolList[iAcol] = acolList.back();
1021  acolList.pop_back();
1022  foundPair = true;
1023  break;
1024  }
1025  } if (foundPair) break;
1026  }
1027  }
1028 
1029  // Usually one unmatched pair left to collapse.
1030  if (colList.size() == 1 && acolList.size() == 1) {
1031  int finalFrom = max( colList[0], acolList[0]);
1032  int finalTo = min( colList[0], acolList[0]);
1033  for (int i = 0; i < size(); ++i)
1034  if ( resolved[i].isFromBeam() ) {
1035  if (resolved[i].col() == finalFrom) resolved[i].col(finalTo);
1036  if (resolved[i].acol() == finalFrom) resolved[i].acol(finalTo);
1037  }
1038  colFrom.push_back(finalFrom);
1039  colTo.push_back(finalTo);
1040 
1041  // Store an (anti)junction when three (anti)coloured daughters.
1042  } else if (hasJunctionBeam && colList.size() == 3
1043  && acolList.size() == 0) {
1044  event.appendJunction( 1, colList[0], colList[1], colList[2]);
1045  junCol[0] = colList[0];
1046  junCol[1] = colList[1];
1047  junCol[2] = colList[2];
1048  } else if (hasJunctionBeam && acolList.size() == 3
1049  && colList.size() == 0) {
1050  event.appendJunction( 2, acolList[0], acolList[1], acolList[2]);
1051  junCol[0] = acolList[0];
1052  junCol[1] = acolList[1];
1053  junCol[2] = acolList[2];
1054 
1055  // Any other nonvanishing values indicate failure.
1056  } else if (colList.size() > 0 || acolList.size() > 0) {
1057  infoPtr->errorMsg("Error in BeamParticle::remnantColours: "
1058  "leftover unmatched colours");
1059  return false;
1060  }
1061 
1062  // Store colour assignment of beam particles.
1063  for (int i = nInit; i < size(); ++i)
1064  event[resolved[i].iPos()].cols( resolved[i].col(), resolved[i].acol() );
1065 
1066  // Done.
1067  return true;
1068 }
1069 
1070 
1071 //--------------------------------------------------------------------------
1072 
1073 // Pick unrescaled x values for beam remnant sharing.
1074 
1075 double BeamParticle::xRemnant( int i) {
1076 
1077  double x = 0.;
1078 
1079  // Hadrons (only used for DIS) rather primitive for now (probably OK).
1080  int idAbs = abs(resolved[i].id());
1081  if (idAbs > 100 && (idAbs/10)%10 != 0) {
1082  x = 1.;
1083 
1084  // Calculation of x of valence quark or diquark, for latter as sum.
1085  } else if (resolved[i].isValence()) {
1086 
1087  // Resolve diquark into sum of two quarks.
1088  int id1 = resolved[i].id();
1089  int id2 = 0;
1090  if (abs(id1) > 1000) {
1091  id2 = (id1 > 0) ? (id1/100)%10 : -(((-id1)/100)%10);
1092  id1 = (id1 > 0) ? id1/1000 : -((-id1)/1000);
1093  }
1094 
1095  // Loop over (up to) two quarks; add their contributions.
1096  for (int iId = 0; iId < 2; ++iId) {
1097  int idNow = (iId == 0) ? id1 : id2;
1098  if (idNow == 0) break;
1099  double xPart = 0.;
1100 
1101  // Assume form (1-x)^a / sqrt(x).
1102  double xPow = valencePowerMeson;
1103  if (isBaryonBeam) {
1104  if (nValKinds == 3 || nValKinds == 1)
1105  xPow = (3. * rndmPtr->flat() < 2.)
1106  ? valencePowerUinP : valencePowerDinP ;
1107  else if (nValence(idNow) == 2) xPow = valencePowerUinP;
1108  else xPow = valencePowerDinP;
1109  }
1110  do xPart = pow2( rndmPtr->flat() );
1111  while ( pow(1. - xPart, xPow) < rndmPtr->flat() );
1112 
1113  // End loop over (up to) two quarks. Possibly enhancement for diquarks.
1114  x += xPart;
1115  }
1116  if (id2 != 0) x *= valenceDiqEnhance;
1117 
1118  // Calculation of x of sea quark, based on companion association.
1119  } else if (resolved[i].isCompanion()) {
1120 
1121  // Find rescaled x value of companion.
1122  double xLeft = 1.;
1123  for (int iInit = 0; iInit < nInit; ++iInit)
1124  if (resolved[iInit].isFromBeam()) xLeft -= resolved[iInit].x();
1125  double xCompanion = resolved[ resolved[i].companion() ].x();
1126  xCompanion /= (xLeft + xCompanion);
1127 
1128  // Now use ansatz q(x; x_c) < N/(x +x_c) to pick x.
1129  do x = pow( xCompanion, rndmPtr->flat()) - xCompanion;
1130  while ( pow( (1. - x - xCompanion) / (1. - xCompanion), companionPower)
1131  * (pow2(x) + pow2(xCompanion)) / pow2(x + xCompanion)
1132  < rndmPtr->flat() );
1133 
1134  // Else a gluon remnant.
1135  // Rarely it is a single gluon remnant, for that case value does not matter.
1136  } else {
1137  do x = pow(xGluonCutoff, 1 - rndmPtr->flat());
1138  while ( pow(1. - x, gluonPower) < rndmPtr->flat() );
1139  }
1140 
1141  return x;
1142 
1143 }
1144 
1145 //--------------------------------------------------------------------------
1146 
1147 // Approximate the remnant mass according to the initiator.
1148 
1149 double BeamParticle::remnantMass(int idIn) {
1150 
1151  int idLight = 2;
1152 
1153  // Hadrons: remove valence flavour masses from the hadron mass, add others.
1154  if ( isHadron() ) {
1155  double mRem = particleDataPtr->m0( id());
1156  int valSign1 = (nValence(idIn) > 0) ? -1 : 1;
1157  mRem += valSign1 * particleDataPtr->m0(idIn);
1158  return mRem;
1159 
1160  // Photons: For gluons, add two light partons to act as valence, otherwise
1161  // add mass of companion. No remnants for unresolved photons.
1162  } else if ( isGamma() ) {
1163  if ( isUnresolved() ) return 0.;
1164  double mRem = (idIn == 21) ? 2*( particleDataPtr->m0( idLight ) ) :
1165  particleDataPtr->m0( idIn );
1166  return mRem;
1167 
1168  } else return 0.;
1169 
1170 }
1171 
1172 //--------------------------------------------------------------------------
1173 
1174 // Check whether room for a one remnant system.
1175 
1176 bool BeamParticle::roomFor1Remnant(double eCM) {
1177 
1178  // If no remnants for the beam return true.
1179  if (!isResolvedGamma) return true;
1180 
1181  // Else check whether room with given kinematics.
1182  double x1 = resolved[0].x();
1183  int id1 = resolved[0].id();
1184  return roomFor1Remnant(id1, x1, eCM);
1185 }
1186 
1187 //--------------------------------------------------------------------------
1188 
1189 // Check whether room for a one remnant system.
1190 
1191 bool BeamParticle::roomFor1Remnant(int id1, double x1, double eCM) {
1192 
1193  // Use u-quark mass as a lower limit for the remnant mass.
1194  int idLight = 2;
1195 
1196  // For gluons minimum requirement two light quarks.
1197  // For quarks need room for one quark of same flavor.
1198  double mRemnant = (id1 == 21) ? 2*( particleDataPtr->m0( idLight ) ) :
1199  particleDataPtr->m0( id1 );
1200  return ( mRemnant < eCM*( 1 - sqrt(x1) ) );
1201 }
1202 
1203 //--------------------------------------------------------------------------
1204 
1205 // Check whether room for two remnants in the event.
1206 
1207 bool BeamParticle::roomFor2Remnants(int id1, double x1, double eCM) {
1208 
1209  // Use u-quark mass as a lower limit for the remnant mass.
1210  int idLight = 2;
1211  double id2 = resolved[0].id();
1212  double x2 = resolved[0].x();
1213 
1214  // For gluons minimum requirement two light quarks.
1215  // For quarks need room for one quark of same flavor.
1216  double m1 = (id1 == 21) ? 2*( particleDataPtr->m0( idLight ) ) :
1217  particleDataPtr->m0( id1 );
1218  double m2 = (id2 == 21) ? 2*( particleDataPtr->m0( idLight ) ) :
1219  particleDataPtr->m0( id2 );
1220  return ( (m1 + m2) < eCM*sqrt( (1.0 - x1)*(1.0 - x2) ) );
1221 }
1222 
1223 
1224 //--------------------------------------------------------------------------
1225 
1226 // Check whether room for two remnants in the event. This used by MPI.
1227 
1228 bool BeamParticle::roomForRemnants(BeamParticle beamOther) {
1229 
1230  // Calculate the invariant mass remaining after MPIs.
1231  double xLeftA = this->xMax(-1);
1232  double xLeftB = beamOther.xMax(-1);
1233  double eCM = infoPtr->eCM();
1234  double Wleft = eCM * sqrt(xLeftA * xLeftB);
1235  double mRemA = 0;
1236  double mRemB = 0;
1237  bool allGluonsA = true;
1238  bool allGluonsB = true;
1239 
1240  // Calculate the total mass of each beam remnant.
1241  for (int i = 0; i < this->size(); ++i)
1242  if ( resolved[i].id() != 21 ) {
1243  allGluonsA = false;
1244  // If initiator a valence, no need for a companion remnant.
1245  if ( resolved[i].companion() < 0 && resolved[i].companion() != -3 )
1246  mRemA += particleDataPtr->m0( resolved[i].id() );
1247  }
1248  for (int i = 0; i < beamOther.size(); ++i)
1249  if ( beamOther[i].id() != 21 ) {
1250  allGluonsB = false;
1251  if ( beamOther[i].companion() < 0 && beamOther[i].companion() != -3 )
1252  mRemB += particleDataPtr->m0( beamOther[i].id() );
1253  }
1254 
1255  // If all initiators are gluons leave room for two light quarks.
1256  // In case of hadrons the mass is taken into account already in xMax().
1257  if ( allGluonsA) mRemA = this->isGamma() ? 2*particleDataPtr->m0(2) : 0.;
1258  if ( allGluonsB) mRemB = beamOther.isGamma() ? 2*particleDataPtr->m0(2) : 0.;
1259 
1260  // If not enough invariant mass left for remnants reject scattering.
1261  if ( Wleft < mRemA + mRemB ) return false;
1262  else return true;
1263 }
1264 
1265 //--------------------------------------------------------------------------
1266 
1267 // Print the list of resolved partons in a beam.
1268 
1269 void BeamParticle::list() const {
1270 
1271  // Header.
1272  cout << "\n -------- PYTHIA Partons resolved in beam -----------------"
1273  << "-------------------------------------------------------------\n"
1274  << "\n i iPos id x comp xqcomp pTfact "
1275  << "colours p_x p_y p_z e m \n";
1276 
1277  // Loop over list of removed partons and print it.
1278  double xSum = 0.;
1279  Vec4 pSum;
1280  for (int i = 0; i < size(); ++i) {
1281  ResolvedParton res = resolved[i];
1282  cout << fixed << setprecision(6) << setw(5) << i << setw(6) << res.iPos()
1283  << setw(8) << res.id() << setw(10) << res.x() << setw(6)
1284  << res.companion() << setw(10) << res.xqCompanion() << setw(10)
1285  << res.pTfactor() << setprecision(3) << setw(6) << res.col()
1286  << setw(6) << res.acol() << setw(11) << res.px() << setw(11)
1287  << res.py() << setw(11) << res.pz() << setw(11) << res.e()
1288  << setw(11) << res.m() << "\n";
1289 
1290  // Also find sum of x and p values.
1291  if (res.companion() != -10) {
1292  xSum += res.x();
1293  pSum += res.p();
1294  }
1295  }
1296 
1297  // Print sum and endline.
1298  cout << setprecision(6) << " x sum:" << setw(10) << xSum
1299  << setprecision(3) << " p sum:"
1300  << setw(11) << pSum.px() << setw(11) << pSum.py() << setw(11)
1301  << pSum.pz() << setw(11) << pSum.e()
1302  << "\n\n -------- End PYTHIA Partons resolved in beam -----------"
1303  << "---------------------------------------------------------------"
1304  << endl;
1305 }
1306 
1307 //--------------------------------------------------------------------------
1308 
1309 // Test whether a lepton is to be considered as unresolved.
1310 
1311 bool BeamParticle::isUnresolvedLepton() {
1312 
1313  // Require record to consist of lepton with full energy plus a photon.
1314  if (!isLeptonBeam || resolved.size() > 2 || resolved[1].id() != 22
1315  || resolved[0].x() < XMINUNRESOLVED) return false;
1316  return true;
1317 
1318 }
1319 
1320 //--------------------------------------------------------------------------
1321 
1322 // For a diffractive system, decide whether to kick out gluon or quark.
1323 
1324 bool BeamParticle::pickGluon(double mDiff) {
1325 
1326  // Relative weight to pick a quark, assumed falling with energy.
1327  double probPickQuark = pickQuarkNorm / pow( mDiff, pickQuarkPower);
1328  return ( (1. + probPickQuark) * rndmPtr->flat() < 1. );
1329 
1330 }
1331 
1332 //--------------------------------------------------------------------------
1333 
1334 // Pick a valence quark at random. (Used for diffractive systems.)
1335 
1336 int BeamParticle::pickValence() {
1337 
1338  // Pick one valence quark at random.
1339  int nTotVal = (isBaryonBeam) ? 3 : 2;
1340  double rnVal = rndmPtr->flat() * nTotVal;
1341  int iVal = (rnVal < 1.) ? 1 : ( (rnVal < 2.) ? 2 : 3 );
1342 
1343  // This valence in slot 1, the rest thereafter.
1344  idVal1 = 0;
1345  idVal2 = 0;
1346  idVal3 = 0;
1347  int iNow = 0;
1348  for (int i = 0; i < nValKinds; ++i)
1349  for (int j = 0; j < nVal[i]; ++j) {
1350  ++iNow;
1351  if (iNow == iVal) idVal1 = idVal[i];
1352  else if ( idVal2 == 0) idVal2 = idVal[i];
1353  else idVal3 = idVal[i];
1354  }
1355 
1356  // Construct diquark if baryon.
1357  if (idVal3 != 0) idVal2 = flavSelPtr->makeDiquark( idVal2, idVal3, idBeam);
1358 
1359  // Done.
1360  return idVal1;
1361 
1362 }
1363 
1364 //--------------------------------------------------------------------------
1365 
1366 // Share lightcone momentum between two remnants in a diffractive system.
1367 
1368 double BeamParticle::zShare( double mDiff, double m1, double m2) {
1369 
1370  // Set up as valence in normal beam so can use xRemnant code.
1371  append(0, idVal1, 0., -3);
1372  append(0, idVal2, 0., -3);
1373  double m2Diff = mDiff*mDiff;
1374 
1375  // Begin to generate z and pT until acceptable solution.
1376  double wtAcc = 0.;
1377  do {
1378  double x1 = xRemnant(0);
1379  double x2 = xRemnant(0);
1380  zRel = max( TINYZREL, min( 1. - TINYZREL, x1 / (x1 + x2) ) );
1381  pair<double, double> gauss2 = rndmPtr->gauss2();
1382  pxRel = diffPrimKTwidth * gauss2.first;
1383  pyRel = diffPrimKTwidth * gauss2.second;
1384 
1385  // Suppress large invariant masses of remnant system.
1386  double mTS1 = m1*m1 + pxRel*pxRel + pyRel*pyRel;
1387  double mTS2 = m2*m2 + pxRel*pxRel + pyRel*pyRel;
1388  double m2Sys = mTS1 / zRel + mTS2 / (1. - zRel);
1389  wtAcc = (m2Sys < m2Diff)
1390  ? pow( 1. - m2Sys / m2Diff, diffLargeMassSuppress) : 0.;
1391  } while (wtAcc < rndmPtr->flat());
1392 
1393  // Done.
1394  return zRel;
1395 
1396 }
1397 
1398 // -------------------------------------------------------------------------
1399 
1400 // Find the colour configuration of the scattered partons
1401 // and update the colour tags to match this configuration.
1402 
1403 void BeamParticle::findColSetup(Event& event) {
1404 
1405  // Reset current colour setup;
1406  colSetup = make_pair(0,0);
1407  cols.clear();
1408  acols.clear();
1409  colUpdates.clear();
1410  nJuncs = 0;
1411  nAjuncs = 0;
1412 
1413  // Setup colour states.
1414  vector<vector <vector <ColState> > > colStates;
1415  colStates.resize(size() + 1);
1416  for (int i = 0; i < size() + 1; ++i) {
1417  colStates[i].resize(2*(i + 1));
1418  for (int j = 0; j < 2*(i+1); ++j)
1419  colStates[i][j].resize(2*(i+1));
1420  }
1421  colStates[0][0][0].nTotal = 1.;
1422 
1423  bool noColouredParticles = true;
1424  // Find all possible multiplets and their degeneracies.
1425  for (int i = 0; i < size(); ++i) {
1426  for (int j = 0; j < int(colStates[i].size()); ++j) {
1427  for (int k = 0; k < int(colStates[i][j].size()); ++k) {
1428  if (colStates[i][j][k].nTotal < 0.5) continue;
1429  int idParton = resolved[i].id();
1430 
1431  // If particle is a quark.
1432  if (idParton > 0 && idParton < 9) {
1433  colStates[i+1][j+1][k].lastSteps.push_back(make_pair(j,k));
1434  colStates[i+1][j+1][k].nTotal += colStates[i][j][k].nTotal;
1435  if (k > 0) {
1436  colStates[i+1][j][k -1].lastSteps.push_back(make_pair(j,k));
1437  colStates[i+1][j][k -1].nTotal += colStates[i][j][k].nTotal;
1438  }
1439 
1440  // Junction combination.
1441  if (j > 0 && allowBeamJunctions) {
1442  colStates[i+1][j - 1][k + 1].lastSteps.push_back(make_pair(j,k));
1443  colStates[i+1][j - 1][k + 1].nTotal += colStates[i][j][k].nTotal;
1444  }
1445  }
1446 
1447  // If particle is an anti quark.
1448  if (idParton < 0 && idParton > -9) {
1449  colStates[i+1][j][k + 1].lastSteps.push_back(make_pair(j,k));
1450  colStates[i+1][j][k + 1].nTotal += colStates[i][j][k].nTotal;
1451  if (j > 0) {
1452  colStates[i+1][j - 1][k].lastSteps.push_back(make_pair(j,k));
1453  colStates[i+1][j - 1][k].nTotal += colStates[i][j][k].nTotal;
1454  }
1455 
1456  // Junction combination.
1457  if (k > 0 && allowBeamJunctions) {
1458  colStates[i+1][j + 1][k - 1].lastSteps.push_back(make_pair(j,k));
1459  colStates[i+1][j + 1][k - 1].nTotal += colStates[i][j][k].nTotal;
1460  }
1461  }
1462 
1463  // If particle is a gluon.
1464  if (idParton == 21) {
1465  colStates[i+1][j + 1][k + 1].lastSteps.push_back(make_pair(j,k));
1466  colStates[i+1][j + 1][k + 1].nTotal += colStates[i][j][k].nTotal;
1467  if (j > 0) {
1468  colStates[i+1][j][k].lastSteps.push_back(make_pair(j,k));
1469  colStates[i+1][j][k].nTotal += colStates[i][j][k].nTotal;
1470  }
1471  if (k > 0) {
1472  colStates[i+1][j][k].lastSteps.push_back(make_pair(j,k));
1473  colStates[i+1][j][k].nTotal += colStates[i][j][k].nTotal;
1474  }
1475  if (j > 0 && k > 0) {
1476  colStates[i+1][j - 1][k - 1].lastSteps.push_back(make_pair(j,k));
1477  colStates[i+1][j - 1][k - 1].nTotal += colStates[i][j][k].nTotal;
1478  }
1479 
1480  // Junction combinations.
1481  if (k > 0 && allowBeamJunctions) {
1482  colStates[i+1][j + 2][k - 1].lastSteps.push_back(make_pair(j,k));
1483  colStates[i+1][j + 2][k - 1].nTotal += colStates[i][j][k].nTotal;
1484  }
1485  if (j > 0 && allowBeamJunctions) {
1486  colStates[i+1][j - 1][k + 2].lastSteps.push_back(make_pair(j,k));
1487  colStates[i+1][j - 1][k + 2].nTotal += colStates[i][j][k].nTotal;
1488  }
1489  if (j > 1 && allowBeamJunctions) {
1490  colStates[i+1][j - 2][k + 1].lastSteps.push_back(make_pair(j,k));
1491  colStates[i+1][j - 2][k + 1].nTotal += colStates[i][j][k].nTotal;
1492  }
1493  if (k > 1 && allowBeamJunctions) {
1494  colStates[i+1][j + 1][k - 2].lastSteps.push_back(make_pair(j,k));
1495  colStates[i+1][j + 1][k - 2].nTotal += colStates[i][j][k].nTotal;
1496  }
1497  }
1498  // If the parton is not a quark or a gluon.
1499  if (idParton != 21 && abs(idParton) > 9) {
1500  colStates[i+1][j][k].lastSteps.push_back(make_pair(j,k));
1501  colStates[i+1][j][k].nTotal += colStates[i][j][k].nTotal;
1502  } else
1503  noColouredParticles = false;
1504  }
1505  }
1506  }
1507 
1508  // Pick a specific multiplet depending on colour weight and saturation.
1509  // Start by calculating the sum of all weights.
1510  double totalSize = 0;
1511  for (int i = 0;i < int(colStates[size()].size()); ++i) {
1512  for (int j = 0;j < int(colStates[size()][i].size()); ++j) {
1513  // Do not allow colour singlet states, since this will overlap
1514  // with diffractive events described elsewhere in PYTHIA.
1515  if (i == 0 && j == 0 && !noColouredParticles) continue;
1516 
1517  double multipletSize = (i + 1) * (j + 1) * (i + j + 2) / 2.;
1518  totalSize += colStates[size()][i][j].nTotal *
1519  multipletSize * exp(-multipletSize / beamSat);
1520  }
1521  }
1522 
1523  // Choose one colour configuration.
1524  double curSize = 0;
1525  double chosenSize = rndmPtr->flat() * totalSize;
1526  for (int i = 0;i < int(colStates[size()].size()); ++i) {
1527  for (int j = 0;j < int(colStates[size()][i].size()); ++j) {
1528 
1529  // Do not allow singlets.
1530  if (i == 0 && j == 0 && !noColouredParticles) continue;
1531 
1532  // Reweight according to multiplet size.
1533  double multipletSize = (i + 1) * (j + 1) * (i + j + 2) / 2.;
1534  curSize += colStates[size()][i][j].nTotal *
1535  multipletSize * exp(-multipletSize / beamSat);
1536  if (curSize > chosenSize) {
1537  colSetup.first = i;
1538  colSetup.second = j;
1539  break;
1540  }
1541  }
1542  if (curSize > chosenSize) break;
1543  }
1544 
1545  // Find the whole colour flow chain.
1546  vector<pair<int, int> > colSetupChain;
1547  colSetupChain.resize(size() + 1);
1548  pair<int,int> curColSetup = make_pair(colSetup.first, colSetup.second);
1549  for (int i = size(); i > 0; --i) {
1550  colSetupChain[i] = curColSetup;
1551  int curColSize = colStates[i][curColSetup.first][curColSetup.second].
1552  lastSteps.size();
1553  int iCurCol = int(rndmPtr->flat() * curColSize);
1554  curColSetup = colStates[i][curColSetup.first][curColSetup.second].
1555  lastSteps[iCurCol];
1556  }
1557  colSetupChain[0] = curColSetup;
1558 
1559  // Update scattered partons to reflect new colour configuration and
1560  // store still unconnected colours and anti-colours.
1561  for (int i = 0; i < size() ; ++i) {
1562 
1563  // If particle is a quark.
1564  if (resolved[i].id() > 0 && resolved[i].id() < 9) {
1565  // Add quark to list of colours.
1566  if (colSetupChain[i].first + 1 == colSetupChain[i + 1].first)
1567  cols.push_back(resolved[i].col());
1568 
1569  // Find anti colour that quark connects to and update event record.
1570  else if (colSetupChain[i].second - 1 == colSetupChain[i + 1].second) {
1571  int iAcol = int(acols.size() * rndmPtr->flat());
1572  int acol = acols[iAcol];
1573  acols.erase(acols.begin() + iAcol);
1574  updateSingleCol(acol, resolved[i].col());
1575  }
1576 
1577  // Else must be junction:
1578  else {
1579  int iCol = int(cols.size() * rndmPtr->flat());
1580  int juncCol = event.nextColTag();
1581  event.appendJunction(1, cols[iCol], resolved[i].col(), juncCol);
1582  event.saveJunctionSize();
1583  acols.push_back(juncCol);
1584  cols.erase(cols.begin() + iCol);
1585  nJuncs++;
1586  }
1587  }
1588 
1589  // If particle is an anti quark.
1590  if (resolved[i].id() < 0 && resolved[i].id() > -9) {
1591  // Add anti quark to list of anti colours
1592  if (colSetupChain[i].second + 1 == colSetupChain[i + 1].second)
1593  acols.push_back(resolved[i].acol());
1594 
1595  // Find colour that anti quark connects to and update event record.
1596  else if (colSetupChain[i].first - 1 == colSetupChain[i + 1].first) {
1597  int iCol = int(cols.size() * rndmPtr->flat());
1598  int col = cols[iCol];
1599  cols.erase(cols.begin() + iCol);
1600  updateSingleCol(col, resolved[i].acol());
1601  }
1602 
1603  // Else it has to be a junction configuration:
1604  else {
1605  int iAcol = int(acols.size() * rndmPtr->flat());
1606  int juncCol = event.nextColTag();
1607  event.appendJunction(2, acols[iAcol], resolved[i].acol(), juncCol);
1608  event.saveJunctionSize();
1609  cols.push_back(juncCol);
1610  acols.erase(acols.begin() + iAcol);
1611  nAjuncs++;
1612  }
1613  }
1614 
1615  // If particle is a gluon.
1616  if (resolved[i].id() == 21) {
1617  // Add gluon to list of colours.
1618  if (colSetupChain[i].first + 1 == colSetupChain[i + 1].first &&
1619  colSetupChain[i].second + 1 == colSetupChain[i + 1].second ) {
1620  acols.push_back(resolved[i].acol());
1621  cols.push_back(resolved[i].col());
1622  }
1623 
1624  // Remove colour and anti colour.
1625  if (colSetupChain[i].first - 1 == colSetupChain[i + 1].first &&
1626  colSetupChain[i].second - 1 == colSetupChain[i + 1].second ) {
1627  // First remove colour.
1628  int iCol = int(cols.size() * rndmPtr->flat());
1629  int col = cols[iCol];
1630  cols.erase(cols.begin() + iCol);
1631  updateSingleCol(col, resolved[i].acol());
1632  // Then remove anti colour.
1633  int iAcol = int(acols.size() * rndmPtr->flat());
1634  int acol = acols[iAcol];
1635  acols.erase(acols.begin() + iAcol);
1636  updateSingleCol(acol, resolved[i].col());
1637  }
1638 
1639  // If the gluon connects to a single end.
1640  // If it is possible to both go to a colour or anti colour pick randomly.
1641  if (colSetupChain[i].first == colSetupChain[i + 1].first &&
1642  colSetupChain[i].second == colSetupChain[i + 1].second ) {
1643  bool removeColour = true;
1644  if (cols.size() > 0 && acols.size() > 0)
1645  removeColour = (rndmPtr->flat() > 0.5);
1646  else if (acols.size() > 0)
1647  removeColour = false;
1648  // Remove colour and add new colour.
1649  if (removeColour) {
1650  int iCol = int(cols.size() * rndmPtr->flat());
1651  int col = cols[iCol];
1652  cols.erase(cols.begin() + iCol);
1653  cols.push_back(resolved[i].col());
1654  updateSingleCol(col, resolved[i].acol());
1655  }
1656  // Else remove anti colour and add new anti colour.
1657  else {
1658  int iAcol = int(acols.size() * rndmPtr->flat());
1659  int acol = acols[iAcol];
1660  acols.erase(acols.begin() + iAcol);
1661  acols.push_back(resolved[i].acol());
1662  updateSingleCol(acol, resolved[i].col());
1663  }
1664  }
1665 
1666  // Junction configuratons.
1667 
1668  // Gluon anti-junction case.
1669  if (colSetupChain[i].first + 2 == colSetupChain[i + 1].first &&
1670  colSetupChain[i].second - 1 == colSetupChain[i + 1].second ) {
1671  int iAcol = int(acols.size() * rndmPtr->flat());
1672  int acol = acols[iAcol];
1673  acols.erase(acols.begin() + iAcol);
1674  int acol3 = event.nextColTag();
1675  event.appendJunction(2,acol,resolved[i].acol(),acol3);
1676  cols.push_back(resolved[i].col());
1677  cols.push_back(acol3);
1678  nAjuncs++;
1679  }
1680 
1681  // Gluon junction case.
1682  if (colSetupChain[i].first - 1 == colSetupChain[i + 1].first &&
1683  colSetupChain[i].second + 2 == colSetupChain[i + 1].second ) {
1684  int iCol = int(cols.size() * rndmPtr->flat());
1685  int col = cols[iCol];
1686  cols.erase(cols.begin() + iCol);
1687  int col3 = event.nextColTag();
1688  event.appendJunction(1,col,resolved[i].col(),col3);
1689  acols.push_back(resolved[i].acol());
1690  acols.push_back(col3);
1691  nJuncs++;
1692  }
1693 
1694  // Gluon anti-junction case.
1695  if (colSetupChain[i].first + 1 == colSetupChain[i + 1].first &&
1696  colSetupChain[i].second - 2 == colSetupChain[i + 1].second ) {
1697  int iAcol = int(acols.size() * rndmPtr->flat());
1698  int acol = acols[iAcol];
1699  acols.erase(acols.begin() + iAcol);
1700  int iAcol2 = int(acols.size() * rndmPtr->flat());
1701  int acol2 = acols[iAcol2];
1702  acols.erase(acols.begin() + iAcol2);
1703  event.appendJunction(2,acol,resolved[i].acol(),acol2);
1704  cols.push_back(resolved[i].col());
1705  nAjuncs++;
1706  }
1707 
1708  // Gluon junction case.
1709  if (colSetupChain[i].first - 2 == colSetupChain[i + 1].first &&
1710  colSetupChain[i].second + 1 == colSetupChain[i + 1].second ) {
1711  int iCol = int(cols.size() * rndmPtr->flat());
1712  int col = cols[iCol];
1713  cols.erase(cols.begin() + iCol);
1714  int iCol2 = int(cols.size() * rndmPtr->flat());
1715  int col2 = cols[iCol2];
1716  cols.erase(cols.begin() + iCol2);
1717  event.appendJunction(1,col,resolved[i].col(),col2);
1718  acols.push_back(resolved[i].acol());
1719  nJuncs++;
1720  }
1721  }
1722  }
1723 
1724  // Done updating event.
1725 }
1726 
1727 // -------------------------------------------------------------------------
1728 
1729 // Add required extra remnant flavour content. Also initial colours.
1730 
1731 bool BeamParticle::remnantFlavoursNew(Event& event) {
1732 
1733  // A baryon will have a junction, unless a diquark is formed later.
1734  hasJunctionBeam = (isBaryon());
1735 
1736  // Store how many hard-scattering partons were removed from beam.
1737  nInit = size();
1738 
1739  // Find remaining valence quarks.
1740  for (int i = 0; i < nValKinds; ++i) {
1741  nValLeft[i] = nVal[i];
1742  for (int j = 0; j < nInit; ++j) if (resolved[j].isValence()
1743  && resolved[j].id() == idVal[i]) --nValLeft[i];
1744  // Add remaining valence quarks to record. Partly temporary values.
1745  for (int k = 0; k < nValLeft[i]; ++k) append(0, idVal[i], 0., -3);
1746  }
1747  int nInitPlusVal = size();
1748 
1749  // Find companion quarks to unmatched sea quarks.
1750  for (int i = 0; i < nInit; ++i)
1751  if (resolved[i].isUnmatched()) {
1752 
1753  // Add companion quark to record; and bookkeep both ways.
1754  append(0, -resolved[i].id(), 0., i);
1755  resolved[i].companion(size() - 1);
1756  }
1757  int beamJunc = 0;
1758  if (isBaryon()) beamJunc = 1;
1759  if (id() < 0) beamJunc = -beamJunc;
1760 
1761  // Count the number of gluons that needs to be added.
1762  int nGluons = (colSetup.first + colSetup.second - (size() - nInit)
1763  +abs( (nJuncs - nAjuncs) - beamJunc)) / 2;
1764  for (int i = 0; i < nGluons; i++) append(0,21,0.,-1);
1765 
1766  // If no other remnants found, add a light q-qbar pair or a photon
1767  // to carry momentum.
1768  if (size() == nInit) {
1769  if (!isHadron())
1770  append(0, 22, 0., -1);
1771  else {
1772  int idRemnant = int(3*rndmPtr->flat())+1;
1773  append(0, -idRemnant, 0., -1);
1774  append(0, idRemnant, 0., -1);
1775  resolved[size()-2].companion(size() - 1);
1776  resolved[size()-1].companion(size() - 2);
1777  }
1778  }
1779 
1780  usedCol = vector<bool>(size(),false);
1781  usedAcol = vector<bool>(size(),false);
1782 
1783  // If at least two valence quarks left in baryon and no junction formed.
1784  // First check if junction already was moved into beam.
1785  nDiffJuncs = nJuncs - nAjuncs - beamJunc;
1786  if (isBaryon() && nInitPlusVal - nInit >= 2 && (
1787  (nDiffJuncs > 0 && beamJunc < 0) ||
1788  (nDiffJuncs < 0 && beamJunc > 0)) ) {
1789 
1790  // If three, pick two at random to form junction, else trivial.
1791  int iQ1 = nInit;
1792  int iQ2 = nInit + 1;
1793  if (nInitPlusVal - nInit == 3) {
1794  double pickDq = 3. * rndmPtr->flat();
1795  if (pickDq > 1.) iQ2 = nInit + 2;
1796  if (pickDq > 2.) iQ1 = nInit + 1;
1797  }
1798 
1799  // Either form di-quark or (anti-)junction.
1800  if (beamJunction) {
1801  // Form antijunction.
1802  if (resolved[iQ1].id() < 0) {
1803 
1804  // Start by finding last colour in the out going particles.
1805  usedAcol[iQ1] = true;
1806  usedAcol[iQ2] = true;
1807 
1808  // Find matching anti-colour.
1809  int acol = findSingleCol(event, true, true);
1810  if ( acol == 0) return false;
1811 
1812  // Make the antijunction.
1813  int newCol1 = event.nextColTag();
1814  int newCol2 = event.nextColTag();
1815  resolved[iQ1].acol(newCol1);
1816  resolved[iQ2].acol(newCol2);
1817  event.appendJunction(2, resolved[iQ1].acol(), resolved[iQ2].acol(),
1818  acol);
1819  nDiffJuncs--;
1820  }
1821 
1822  // Form Junction.
1823  else {
1824  // Start by finding last colour in the out going particles.
1825  usedCol[iQ1] = true;
1826  usedCol[iQ2] = true;
1827 
1828  // Find matching colour.
1829  int col = findSingleCol(event, false, true);
1830  if (col == 0) return false;
1831 
1832  // Make the junction.
1833  int newCol1 = event.nextColTag();
1834  int newCol2 = event.nextColTag();
1835  resolved[iQ1].col(newCol1);
1836  resolved[iQ2].col(newCol2);
1837  event.appendJunction(1,resolved[iQ1].col(),resolved[iQ2].col(),col);
1838  nDiffJuncs++;
1839  }
1840 
1841  // Form diquark.
1842  } else {
1843 
1844  // Pick spin 0 or 1 according to SU(6) wave function factors.
1845  int idDq = flavSelPtr->makeDiquark( resolved[iQ1].id(),
1846  resolved[iQ2].id(), idBeam);
1847 
1848  // Overwrite with diquark flavour and remove one slot. No more junction.
1849  if (nInitPlusVal - nInit == 3)
1850  resolved[nInit + 2].id( resolved[3 * nInit + 3 - iQ2 - iQ1].id() );
1851  resolved[nInit].id(idDq);
1852  resolved.erase(resolved.begin() + nInit + 1);
1853  hasJunctionBeam = false;
1854 
1855  // Di-quark changes the baryon number.
1856  if (idDq > 0) nDiffJuncs++;
1857  else nDiffJuncs--;
1858  }
1859  }
1860 
1861  // Form anti-junction out of any beam remnants if needed.
1862  while (nDiffJuncs > 0) {
1863  int acol1 = findSingleCol(event, true, false);
1864  int acol2 = findSingleCol(event, true, false);
1865  int acol3 = findSingleCol(event, true, true);
1866  event.appendJunction(2,acol1,acol2,acol3);
1867  nDiffJuncs--;
1868  }
1869  // Form junction out of any beam remnants if needed.
1870  while (nDiffJuncs < 0) {
1871  int col1 = findSingleCol(event, false, false);
1872  int col2 = findSingleCol(event, false, false);
1873  int col3 = findSingleCol(event, false, true);
1874  event.appendJunction(1,col1,col2,col3);
1875  nDiffJuncs++;
1876  }
1877 
1878  // Set remaining colours first in random order.
1879  for (int j = 0;j < NRANDOMTRIES; ++j) {
1880  int i = int(rndmPtr->flat() * (size() - nInit) + nInit );
1881  // Check if resolved has colour.
1882  if ( resolved[i].hasCol() && !usedCol[i]) {
1883  usedCol[i] = true;
1884  int acol = findSingleCol(event,true,true);
1885  if ( acol == 0) return false;
1886  resolved[i].col(acol);
1887  }
1888  // Check if resolved has anti colour.
1889  if ( resolved[i].hasAcol() && !usedAcol[i]) {
1890  usedAcol[i] = true;
1891  int col = findSingleCol(event, false, true);
1892  if (col == 0) return false;
1893  resolved[i].acol(col);
1894  }
1895  }
1896 
1897  // Add all missed colours from the random assignment.
1898  for (int i = nInit;i < size();i++) {
1899  // Check if resolved has colour.
1900  if ( resolved[i].hasCol() && !usedCol[i]) {
1901  usedCol[i] = true;
1902  int acol = findSingleCol(event,true,true);
1903  if ( acol == 0) return false;
1904  resolved[i].col(acol);
1905  }
1906  // Check if resolved has anti colour.
1907  if ( resolved[i].hasAcol() && !usedAcol[i]) {
1908  usedAcol[i] = true;
1909  int col = findSingleCol(event, false, true);
1910  if (col == 0) return false;
1911  resolved[i].acol(col);
1912  }
1913  }
1914 
1915  // Need to end in a colour singlet.
1916  if (cols.size() != 0 || acols.size() != 0) {
1917  infoPtr->errorMsg("Error in BeamParticle::RemnantFlavours: "
1918  "Colour not conserved in beamRemnants");
1919  return false;
1920  }
1921 
1922  // Set initiator and remnant masses.
1923  for (int i = 0; i < size(); ++i) {
1924  if (i < nInit) resolved[i].m(0.);
1925  else resolved[i].m( particleDataPtr->m0( resolved[i].id() ) );
1926  }
1927 
1928  // For debug purposes: reject beams with resolved junction topology.
1929  if (hasJunctionBeam && !allowJunction) return false;
1930 
1931  // Done.
1932  return true;
1933 
1934 }
1935 
1936 //--------------------------------------------------------------------------
1937 
1938 // Set initial colours.
1939 
1940 void BeamParticle::setInitialCol(Event& event) {
1941 
1942  // Set beam colours equal to those in the event record.
1943  for (int i = 0;i < size(); ++i) {
1944  if (event[resolved[i].iPos()].col() != 0)
1945  resolved[i].col(event[resolved[i].iPos()].col());
1946  if (event[resolved[i].iPos()].acol() != 0)
1947  resolved[i].acol(event[resolved[i].iPos()].acol());
1948  }
1949 }
1950 
1951 //--------------------------------------------------------------------------
1952 
1953 // Find a single (anti-) colour in the beam, if it is a
1954 // beam remnant set the new colour.
1955 
1956 int BeamParticle::findSingleCol(Event& event, bool isAcol,
1957  bool useHardScatters) {
1958 
1959  // Look in the already scattered parton list.
1960  if (useHardScatters) {
1961  if (isAcol) {
1962  if (acols.size() > 0) {
1963  int iAcol = int(acols.size() * rndmPtr->flat());
1964  int acol = acols[iAcol];
1965  acols.erase(acols.begin() + iAcol);
1966  return acol;
1967  }
1968  } else {
1969  if (int(cols.size()) > 0) {
1970  int iCol = int(cols.size() * rndmPtr->flat());
1971  int col = cols[iCol];
1972  cols.erase(cols.begin() + iCol);
1973  return col;
1974  }
1975  }
1976  }
1977 
1978  // Look inside the beam remnants.
1979  if (isAcol) {
1980  for (int i = 0;i < NMAX; ++i) {
1981  int iBeam = int((size() - nInit) * rndmPtr->flat()) + nInit;
1982  if (resolved[iBeam].hasAcol() && !usedAcol[iBeam]) {
1983  int acol = event.nextColTag();
1984  resolved[iBeam].acol(acol);
1985  usedAcol[iBeam] = true;
1986  return acol;
1987  }
1988  }
1989  } else {
1990  for (int i = 0; i < NMAX; ++i) {
1991  int iBeam = int((size() - nInit) * rndmPtr->flat()) + nInit;
1992  if (resolved[iBeam].hasCol() && !usedCol[iBeam]) {
1993  int col = event.nextColTag();
1994  resolved[iBeam].col(col);
1995  usedCol[iBeam] = true;
1996  return col;
1997  }
1998  }
1999  }
2000 
2001  // Return 0 if no particle was found.
2002  infoPtr->errorMsg("Error in BeamParticle::findSingleCol: "
2003  "could not find matching anti colour");
2004  return 0;
2005 }
2006 
2007 //--------------------------------------------------------------------------
2008 
2009 // Update list of all colours in beam.
2010 
2011 void BeamParticle::updateCol(vector<pair<int,int> > colourChanges) {
2012 
2013  for (int iCol = 0;iCol < int(colourChanges.size()); ++iCol) {
2014  int oldCol = colourChanges[iCol].first;
2015  int newCol = colourChanges[iCol].second;
2016 
2017  // Update acols and cols.
2018  for (int i = 0;i < int(cols.size()); ++i)
2019  if (cols[i] == oldCol) cols[i] = newCol;
2020  for (int i = 0;i < int(acols.size()); ++i)
2021  if (acols[i] == oldCol) acols[i] = newCol;
2022 
2023  // Update resolved partons colours.
2024  for (int i = 0;i < int(resolved.size()); ++i) {
2025  if (resolved[i].acol() == oldCol) resolved[i].acol(newCol);
2026  if (resolved[i].col() == oldCol) resolved[i].col(newCol);
2027  }
2028  }
2029  return;
2030 }
2031 
2032 //--------------------------------------------------------------------------
2033 
2034 void BeamParticle::updateSingleCol(int oldCol, int newCol) {
2035 
2036  // Update acols and cols.
2037  for (int i = 0; i < int(cols.size()); ++i)
2038  if (cols[i] == oldCol) cols[i] = newCol;
2039 
2040  for ( int i = 0; i < int(acols.size()); ++i)
2041  if (acols[i] == oldCol) acols[i] = newCol;
2042 
2043  // Update resolved partons colours.
2044  for (int i = 0; i < int(resolved.size()); ++i) {
2045  if (resolved[i].acol() == oldCol) resolved[i].acol(newCol);
2046  if (resolved[i].col() == oldCol) resolved[i].col(newCol);
2047  }
2048 
2049  colUpdates.push_back(make_pair(oldCol,newCol));
2050 }
2051 
2052 //==========================================================================
2053 
2054 } // end namespace Pythia8
Definition: AgUStep.h:26