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