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