- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / src / BeamRemnants.cxx
1 // BeamRemnants.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 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 // BeamRemnants class.
8
9 #include "BeamRemnants.h"
10
11 namespace Pythia8 {
12
13 //==========================================================================
14
15 // The BeamDipole class is purely internal to reconnectColours.
16
17 class BeamDipole {
18
19 public:
20
21   // Constructor.
22   BeamDipole( int colIn = 0, int iColIn = 0, int iAcolIn = 0) 
23     : col(colIn), iCol(iColIn), iAcol(iAcolIn) {}  
24
25   // Members.
26   int    col, iCol, iAcol;
27   double p1p2;
28  
29 };
30
31 //==========================================================================
32
33 // The BeamRemnants class.
34
35 //--------------------------------------------------------------------------
36
37 // Constants: could be changed here if desired, but normally should not.
38 // These are of technical nature, as described for each.
39
40 // If same (anti)colour appears twice in final state, repair or reject.
41 const bool   BeamRemnants::ALLOWCOLOURTWICE = true; 
42
43 // Maximum number of tries to match colours and kinematics in the event.
44 const int    BeamRemnants::NTRYCOLMATCH     = 10; 
45 const int    BeamRemnants::NTRYKINMATCH     = 10;  
46
47 // Overall correction step for energy-momentum conservation; only
48 // becomes relevant in rescattering scenarios when FSR dipole emissions
49 // and primordial kT is added. Should hopefully not be needed.
50 const bool   BeamRemnants::CORRECTMISMATCH  = false; 
51
52 //--------------------------------------------------------------------------
53
54 // Initialization.
55
56 bool BeamRemnants::init( Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn, 
57   BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, 
58   PartonSystems* partonSystemsPtrIn) {
59
60   // Save pointers.
61   infoPtr             = infoPtrIn;
62   rndmPtr             = rndmPtrIn;
63   beamAPtr            = beamAPtrIn;
64   beamBPtr            = beamBPtrIn;
65   partonSystemsPtr    = partonSystemsPtrIn;
66
67   // Width of primordial kT distribution.
68   doPrimordialKT      = settings.flag("BeamRemnants:primordialKT");
69   primordialKTsoft    = settings.parm("BeamRemnants:primordialKTsoft");
70   primordialKThard    = settings.parm("BeamRemnants:primordialKThard");
71   primordialKTremnant = settings.parm("BeamRemnants:primordialKTremnant");
72   halfScaleForKT      = settings.parm("BeamRemnants:halfScaleForKT");
73   halfMassForKT       = settings.parm("BeamRemnants:halfMassForKT");
74
75   // Handling of rescattering kinematics uncertainties from primodial kT.
76   allowRescatter      = settings.flag("MultipleInteractions:allowRescatter");
77   doRescatterRestoreY = settings.flag("BeamRemnants:rescatterRestoreY");
78
79   // Parameters for colour reconnection scenario, partly borrowed from
80   // multiple interactions not to introduce too many new ones.
81   doReconnect         = settings.flag("BeamRemnants:reconnectColours");
82   reconnectRange      = settings.parm("BeamRemnants:reconnectRange");
83   pT0Ref              = settings.parm("MultipleInteractions:pT0Ref");
84   ecmRef              = settings.parm("MultipleInteractions:ecmRef");
85   ecmPow              = settings.parm("MultipleInteractions:ecmPow");
86
87   // Total and squared CM energy at nominal energy.
88   eCM                 = infoPtr->eCM();
89   sCM                 = eCM * eCM;
90
91   // The MI pT0 smoothening scale and its reconnection-strength combination.
92   pT0                 = pT0Ref * pow(eCM / ecmRef, ecmPow);
93   pT20Rec             = pow2(reconnectRange * pT0); 
94   
95   // Done.
96   return true;
97
98 }
99
100 //--------------------------------------------------------------------------
101
102 // Select the flavours/kinematics/colours of the two beam remnants. 
103 // Notation: iPar = all partons, iSys = matched systems of two beams,
104 //           iRem = additional partons in remnants.
105
106 bool BeamRemnants::add( Event& event) {
107
108   // Update to current CM energy.
109   eCM     = infoPtr->eCM();
110   sCM     = eCM * eCM;
111
112   // Check that flavour bookkept in event and in beam particles agree.
113   for (int i = 0; i < beamAPtr->size(); ++i) {
114     int j = (*beamAPtr)[i].iPos();
115     if ((*beamAPtr)[i].id() != event[j].id()) {
116       infoPtr->errorMsg("Error in BeamRemnants::add: "
117         "event and beam flavours do not match"); 
118       return false;
119     }
120   }
121   for (int i = 0; i < beamBPtr->size(); ++i) {
122     int j =  (*beamBPtr)[i].iPos();
123     if ((*beamBPtr)[i].id() != event[j].id()) {
124       infoPtr->errorMsg("Error in BeamRemnants::add: "
125         "event and beam flavours do not match"); 
126       return false;
127     }
128   }
129
130   // Number of scattering subsystems. Size of event record before treatment.
131   nSys    = partonSystemsPtr->sizeSys();
132   oldSize = event.size();
133
134   // Add required extra remnant flavour content. 
135   // Start all over if fails (in option where junctions not allowed).
136   if ( !beamAPtr->remnantFlavours(event) 
137     || !beamBPtr->remnantFlavours(event) ) {
138     infoPtr->errorMsg("Error in BeamRemnants::add:"
139       " remnant flavour setup failed"); 
140     return false;
141   }
142
143   // Do the kinematics of the collision subsystems and two beam remnants.
144   if (!setKinematics(event)) return false;
145
146   // Allow colour reconnections.
147   if (doReconnect) reconnectColours(event);
148
149   // Save current modifiable colour configuration for fast restoration. 
150   vector<int> colSave;
151   vector<int> acolSave;
152   for (int i = oldSize; i < event.size(); ++i) {
153     colSave.push_back( event[i].col() );
154     acolSave.push_back( event[i].acol() );
155   }
156   event.saveJunctionSize();
157   
158   // Allow several tries to match colours of initiators and remnants.
159   // Frequent "failures" since shortcutting colours separately on
160   // the two event sides may give "colour singlet gluons" etc.
161   bool physical = true;
162   for (int iTry = 0; iTry < NTRYCOLMATCH; ++iTry) {
163     physical = true;
164
165     // Reset list of colour "collapses" (transformations).
166     colFrom.resize(0);
167     colTo.resize(0);      
168
169     // First process each set of beam colours on its own.
170     if (!beamAPtr->remnantColours(event, colFrom, colTo)) 
171       physical = false;
172     if (!beamBPtr->remnantColours(event, colFrom, colTo)) 
173       physical = false;
174
175     // Then check that colours and anticolours are matched in whole event.
176     if ( physical && !checkColours(event) ) physical = false;     
177
178     // If no problems then done, else restore and loop.
179     if (physical) break;
180     for (int i = oldSize; i < event.size(); ++i) 
181       event[i].cols( colSave[i - oldSize], acolSave[i - oldSize] );
182     event.restoreJunctionSize();
183     infoPtr->errorMsg("Warning in BeamRemnants::add:"
184       " colour tracing failed; will try again"); 
185   }
186
187   // If no solution after several tries then failed.
188   if (!physical) {
189     infoPtr->errorMsg("Error in BeamRemnants::add:"
190       " colour tracing failed after several attempts"); 
191     return false;
192   }
193
194   // Done.
195   return true;
196
197 }
198
199 //--------------------------------------------------------------------------
200
201 // Set up trial transverse and longitudinal kinematics for each beam 
202 // separately. Final decisions involve comparing the two beams.
203
204 bool BeamRemnants::setKinematics( Event& event) {
205
206   // References to beams to simplify indexing.
207   BeamParticle& beamA = *beamAPtr;  
208   BeamParticle& beamB = *beamBPtr;  
209
210   // Nothing to do for lepton-lepton scattering with all energy already used.
211   if ( beamA.isUnresolvedLepton() && beamB.isUnresolvedLepton() ) 
212     return true; 
213
214   // Check that has not already used up beams.
215   if ( (!beamA.isLepton() && beamA.xMax(-1) <= 0.)
216     || (!beamB.isLepton() && beamB.xMax(-1) <= 0.) ) {
217     infoPtr->errorMsg("Error in BeamRemnants::setKinematics:"
218       " no momentum left for beam remnants"); 
219     return false;
220   }
221
222   // Last beam-status particles. Offset relative to normal beam locations.
223   int nBeams   = 3;
224   for (int i = 3; i < event.size(); ++i) 
225     if (event[i].statusAbs() < 20) nBeams = i + 1; 
226   int nOffset  = nBeams - 3; 
227
228   // Reserve space for extra information on the systems and beams.
229   int nMaxBeam = max( beamA.size(), beamB.size() );
230   vector<double> sHatSys(nMaxBeam);
231   vector<double> kTwidth(nMaxBeam);
232   vector<double> kTcomp(nMaxBeam);
233   vector<RotBstMatrix> Msys(nSys);
234
235   // Loop over all subsystems. Default values. Find invariant mass.
236   double kTcompSumA   = 0.;
237   double kTcompSumB   = 0.;
238   for (int iSys = 0; iSys < nSys; ++iSys) { 
239     double kTwidthNow = 0.;
240     double mHatDamp   = 1.;
241     int iInA          = partonSystemsPtr->getInA(iSys);  
242     int iInB          = partonSystemsPtr->getInB(iSys);  
243     double sHatNow    = (event[iInA].p() + event[iInB].p()).m2Calc();
244
245     // Allow primordial kT reduction for small-mass and small-pT systems
246     // (for hardest interaction pT -> renormalization scale so also 2 -> 1).
247     if (doPrimordialKT) {
248       double mHat     = sqrt(sHatNow);
249       mHatDamp        = mHat / (mHat + halfMassForKT);
250       double scale    = (iSys == 0) ? infoPtr->QRen() : infoPtr->pTMI(iSys);
251       kTwidthNow      = ( (halfScaleForKT * primordialKTsoft
252       + scale * primordialKThard) / (halfScaleForKT + scale) ) * mHatDamp;
253     }
254
255     // Store properties of compensation systems and total compensation power.
256     // Rescattered partons do not compensate, but may be massive.
257     sHatSys[iSys] = sHatNow;
258     kTwidth[iSys] = kTwidthNow ;
259     kTcomp[iSys]  = mHatDamp;
260     if (beamA[iSys].isFromBeam()) kTcompSumA += mHatDamp;
261     else beamA[iSys].m( event[iInA].m() );
262     if (beamB[iSys].isFromBeam()) kTcompSumB += mHatDamp;
263     else beamB[iSys].m( event[iInB].m() );
264   } 
265
266   // Primordial kT and compensation power among remnants.
267   double kTwidthNow = (doPrimordialKT) ? primordialKTremnant : 0.;    
268   for (int iRem = nSys; iRem < nMaxBeam; ++iRem) { 
269     sHatSys[iRem] = 0.;
270     kTwidth[iRem] = kTwidthNow ;
271     kTcomp[iRem]  = 1.;
272   }     
273   kTcompSumA += beamA.size() - nSys;
274   kTcompSumB += beamB.size() - nSys;
275
276   // Allow ten tries to construct kinematics (but normally works first).
277   bool physical;
278   double xSum[2], xInvM[2], w2Beam[2], wPosRem, wNegRem, w2Rem;
279   for (int iTry = 0; iTry < NTRYKINMATCH; ++iTry) {
280     physical = true;
281
282     // Loop over the two beams. 
283     for (int iBeam = 0; iBeam < 2; ++iBeam) {
284       BeamParticle& beam = (iBeam == 0) ? beamA : beamB;
285       int nPar = beam.size();
286  
287       // Generate Gaussian pT for initiator partons inside hadrons.
288       // Store/restore rescattered parton momenta before primordial kT. 
289       if (beam.isHadron() && doPrimordialKT) { 
290         double pxSum = 0.;
291         double pySum = 0.;
292         for (int iPar = 0; iPar < nPar; ++iPar) { 
293           if ( beam[iPar].isFromBeam() ) {
294             pair<double, double> gauss2 = rndmPtr->gauss2();
295             double px = kTwidth[iPar] * gauss2.first;
296             double py = kTwidth[iPar] * gauss2.second;
297             beam[iPar].px(px);
298             beam[iPar].py(py);
299             pxSum += px;
300             pySum += py;
301           } else {
302             int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar)    
303                                      : partonSystemsPtr->getInB(iPar);    
304             beam[iPar].p( event[iInAB].p() );
305           }
306         }  
307
308         // Share recoil between all initiator partons, rescatterers excluded.
309         double kTcompSum = (iBeam == 0) ? kTcompSumA : kTcompSumB;
310         for (int iPar = 0; iPar < nPar; ++iPar) 
311         if (beam[iPar].isFromBeam() ) {
312           beam[iPar].px( beam[iPar].px() - pxSum * kTcomp[iPar] / kTcompSum );
313           beam[iPar].py( beam[iPar].py() - pySum * kTcomp[iPar] / kTcompSum );
314         }
315
316       // Without primordial kT: still need to store rescattered partons.
317       } else if (beam.isHadron()) {
318         for (int iPar = 0; iPar < nPar; ++iPar)  
319         if ( !beam[iPar].isFromBeam() ) {
320           int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar)    
321                                    : partonSystemsPtr->getInB(iPar);       
322           beam[iPar].p( event[iInAB].p() );
323         }   
324       }
325
326       // Pick unrescaled x values for remnants. Sum up (unscaled) p+ and p-.
327       xSum[iBeam]  = 0.;
328       xInvM[iBeam] = 0.;
329       for (int iRem = nSys; iRem < nPar; ++iRem) {
330         double xPrel = beam.xRemnant( iRem);
331         beam[iRem].x(xPrel);
332         xSum[iBeam]  += xPrel;
333         xInvM[iBeam] += beam[iRem].mT2()/xPrel;  
334       }
335
336       // Squared transverse mass for each beam, using lightcone x.
337       w2Beam[iBeam] = xSum[iBeam] * xInvM[iBeam];
338   
339     // End separate treatment of the two beams. 
340     } 
341
342     // Recalculate kinematics of initiator systems with primordial kT.
343     wPosRem = eCM;
344     wNegRem = eCM;
345     for (int iSys = 0; iSys < nSys; ++iSys) { 
346       int iA          = beamA[iSys].iPos();
347       int iB          = beamB[iSys].iPos();
348       double sHat     = sHatSys[iSys];
349
350       // Classify system: rescattering on either or both sides?
351       bool normalA    = beamA[iSys].isFromBeam();
352       bool normalB    = beamB[iSys].isFromBeam();
353       bool normalSys  = normalA && normalB;
354       bool doubleRes  = !normalA && !normalB;
355
356       // Check whether final-state system momentum matches initial-state one.
357       if (allowRescatter && CORRECTMISMATCH) {
358         Vec4 pInitial = event[iA].p() + event[iB].p();
359         Vec4 pFinal;
360         for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
361           int iAB      = partonSystemsPtr->getOut(iSys, iMem);
362           if (event[iAB].isFinal()) pFinal += event[iAB].p();
363         }
364
365         // Scale down primordial kT from side A if p+ increased.
366         if (normalA && pFinal.pPos() > pInitial.pPos())  
367           beamA[iSys].scalePT( pInitial.pPos() / pFinal.pPos() ); 
368
369         // Scale down primordial kT from side B if p- increased.
370         if (normalB && pFinal.pNeg() > pInitial.pNeg()) 
371           beamB[iSys].scalePT( pInitial.pNeg() / pFinal.pNeg() ); 
372       }
373
374       // Rescatter: possible change in sign of lightcone momentum of a
375       //            rescattered parton. If this happens, try to pick
376       //            new primordial kT values
377       if (allowRescatter 
378          && (event[iA].pPos() / beamA[iSys].pPos() < 0 
379          ||  event[iB].pNeg() / beamB[iSys].pNeg() < 0) ) {
380         infoPtr->errorMsg("Warning in BeamRemnants::setKinematics:"
381           " change in lightcone momentum sign; retrying kinematics"); 
382         physical = false;
383         break;
384       }
385
386       // Begin kinematics of partons after primordial kT has been added.
387       double sHatTAft = sHat + pow2( beamA[iSys].px() + beamB[iSys].px()) 
388                              + pow2( beamA[iSys].py() + beamB[iSys].py()); 
389       double w2A      = beamA[iSys].mT2();
390       double w2B      = beamB[iSys].mT2();
391       double w2Diff   = sHatTAft - w2A - w2B;
392       double lambda   = pow2(w2Diff) - 4. * w2A * w2B;
393
394       // Too large transverse momenta means that kinematics will not work.
395       if (lambda <= 0.) {
396         physical      = false;
397         break;
398       }  
399       double lamRoot  = sqrtpos( lambda );
400
401       // Mirror solution if the two incoming have reverse rapidity ordering.
402       if (allowRescatter && doubleRes && (event[iA].pPos() * event[iB].pNeg() 
403         < event[iA].pNeg() * event[iB].pPos()) ) lamRoot = -lamRoot;
404
405       // Two procedures, which agree for normal scattering, separate here.
406       // First option keeps rapidity (and mass) of system unchanged by 
407       // primordial kT, by modification of rescattered parton. 
408       if (normalSys || doRescatterRestoreY || doubleRes) {
409
410         // Find kinematics of initial system: transverse mass, and 
411         // longitudinal momentum carried by all or rescattered partons.
412         double sHatTBef = sHat;
413         double wPosBef, wNegBef, wPosBefRes, wNegBefRes;
414         // Normal scattering.
415         if (normalSys) {
416           wPosBef       = beamA[iSys].x() * eCM;
417           wNegBef       = beamB[iSys].x() * eCM;
418           wPosBefRes    = 0.;
419           wNegBefRes    = 0.;      
420         // Rescattering on side A.
421         } else if (normalB) {
422           sHatTBef     += event[iA].pT2();
423           wPosBef       = event[iA].pPos();
424           wNegBef       = event[iA].pNeg() + beamB[iSys].x() * eCM;
425           wPosBefRes    = beamA[iSys].pPos();
426           wNegBefRes    = beamA[iSys].pNeg();
427         // Rescattering on side B.
428         } else if (normalA) {
429           sHatTBef     += event[iB].pT2(); 
430           wPosBef       = beamA[iSys].x() * eCM + event[iB].pPos();
431           wNegBef       = event[iB].pNeg();
432           wPosBefRes    = beamB[iSys].pPos();
433           wNegBefRes    = beamB[iSys].pNeg();
434         // Rescattering on both sides.
435         } else {
436           sHatTBef     += pow2( event[iA].px() + event[iB].px())
437                         + pow2( event[iA].py() + event[iB].py());
438           wPosBef       = event[iA].pPos() + event[iB].pPos();
439           wNegBef       = event[iA].pNeg() + event[iB].pNeg();
440           wPosBefRes    = beamA[iSys].pPos() + beamB[iSys].pPos();
441           wNegBefRes    = beamA[iSys].pNeg() + beamB[iSys].pNeg();
442         }
443
444         // Rescale outgoing momenta to keep same mass and rapidity of system
445         // as originally, and solve for kinematics.
446         double rescale  = sqrt(sHatTAft / sHatTBef);
447         double wPosAft  = rescale * wPosBef;
448         double wNegAft  = rescale * wNegBef;
449         wPosRem        -= wPosAft - wPosBefRes;
450         wNegRem        -= wNegAft - wNegBefRes; 
451         double wPosA    = 0.5 * (sHatTAft + w2A - w2B + lamRoot) / wNegAft;
452         double wNegB    = 0.5 * (sHatTAft + w2B - w2A + lamRoot) / wPosAft;
453   
454         // Store modified beam parton momenta.
455         beamA[iSys].e(  0.5 * (wPosA + w2A / wPosA) );
456         beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) );        
457         beamB[iSys].e(  0.5 * (w2B / wNegB + wNegB) );
458         beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) );
459
460       // Second option keeps rescattered parton (and system mass) unchanged
461       // by primordial kT, by modification of system rapidity. 
462       } else {
463
464         // Rescattering on side A: preserve already scattered parton.
465         if (normalB) {
466           double wPosA  = beamA[iSys].pPos(); 
467           double wNegB  = 0.5 * (w2Diff + lamRoot) / wPosA;
468           beamB[iSys].e(  0.5 * (w2B / wNegB + wNegB) );
469           beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) );
470           wPosRem      -= w2B / wNegB;
471           wNegRem      -= wNegB; 
472      
473
474         // Rescattering on side B: preserve already scattered parton.
475         } else if (normalA) {
476           double wNegB  = beamB[iSys].pNeg(); 
477           double wPosA  = 0.5 * (w2Diff + lamRoot) / wNegB;
478           beamA[iSys].e(  0.5 * (wPosA + w2A / wPosA) );
479           beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) );
480           wPosRem      -= wPosA;
481           wNegRem      -= w2A / wPosA; 
482
483         // Primordial kT in double rescattering does change the mass of 
484         // the system without any possible fix in this framework, which 
485         // leads to incorrect boosts. Current choice is therefore always
486         // to handle it with the first procedure, where mass is retained.
487         } else {
488         }
489       }
490
491       // Construct system rotation and boost caused by primordial kT.
492       Msys[iSys].reset();
493       Msys[iSys].toCMframe( event[iA].p(), event[iB].p() );
494       Msys[iSys].fromCMframe( beamA[iSys].p(), beamB[iSys].p() );
495
496       // Boost rescattered partons in subsequent beam A list.
497       for (int iSys2 = iSys + 1; iSys2 < nSys; ++iSys2) { 
498         if (!beamA[iSys2].isFromBeam()) {
499           int iBefResc = event[ beamA[iSys2].iPos() ].mother1();     
500           for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem)
501           if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) {
502             Vec4 pTemp = event[iBefResc].p();
503             pTemp.rotbst( Msys[iSys] );
504             beamA[iSys2].p( pTemp );
505           }
506         } 
507
508         // Boost rescattered partons in subsequent beam B list.
509         if (!beamB[iSys2].isFromBeam()) {
510           int iBefResc = event[ beamB[iSys2].iPos() ].mother1();
511           for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem)
512           if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) {
513             Vec4 pTemp = event[iBefResc].p();
514             pTemp.rotbst( Msys[iSys] );
515             beamB[iSys2].p( pTemp );
516           }
517         } 
518       }
519     }
520
521     // Check that remaining momentum is enough for remnants.
522     if (wPosRem < 0. || wNegRem < 0.) physical = false;
523     w2Rem = wPosRem * wNegRem;
524     if (sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1]))
525       physical = false;
526
527     // End of loop over ten tries. Do not loop when solution found.  
528     if (physical) break;
529   }
530
531   // If no solution after ten tries then failed.
532   if (!physical) {
533     infoPtr->errorMsg("Error in BeamRemnants::setKinematics:"
534       " kinematics construction failed"); 
535     return false;
536   }
537
538   // For successful initiator kinematics process whole systems.
539   Vec4 pSumOut;
540   for (int iSys = 0; iSys < nSys; ++iSys) {
541
542     // Copy initiators and their systems and boost them accordingly.
543     // Update subsystem and beams info on new positions of partons.
544     // Update daughter info of mothers, i.e. of beams, for hardest interaction.
545     if (beamA[iSys].isFromBeam()) {
546       int iA       = beamA[iSys].iPos();
547       int iAcopy   = event.copy(iA, -61);
548       event[iAcopy].rotbst(Msys[iSys]);
549       partonSystemsPtr->setInA(iSys, iAcopy);
550       beamA[iSys].iPos( iAcopy);
551       if (iSys == 0) { 
552         int mother = event[iAcopy].mother1();
553         event[mother].daughter1(iAcopy);      
554       }   
555     }
556     if (beamB[iSys].isFromBeam()) {
557       int iB       = beamB[iSys].iPos();
558       int iBcopy   = event.copy(iB, -61);
559       event[iBcopy].rotbst(Msys[iSys]);
560       partonSystemsPtr->setInB(iSys, iBcopy);
561       beamB[iSys].iPos( iBcopy);
562       if (iSys == 0) { 
563         int mother = event[iBcopy].mother1();
564         event[mother].daughter1(iBcopy);      
565       }    
566     }
567
568     for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
569       int iAB      = partonSystemsPtr->getOut(iSys, iMem);
570       if (event[iAB].isFinal()) {
571         int iABcopy = event.copy(iAB, 62);
572         event[iABcopy].rotbst(Msys[iSys]); 
573         partonSystemsPtr->setOut(iSys, iMem, iABcopy);
574         pSumOut   += event[iABcopy].p();
575       }
576     }
577
578   }
579
580   // Colour dipoles spanning systems gives mismatch between FSR recoils
581   // and primordial kT boosts.
582   if (allowRescatter && CORRECTMISMATCH) { 
583
584     // Find summed pT of beam remnants = - wanted pT of systems.
585     double pxBeams = 0.;
586     double pyBeams = 0.;
587     for (int iRem = nSys; iRem < beamA.size(); ++iRem) { 
588       pxBeams     += beamA[iRem].px();
589       pyBeams     += beamA[iRem].py();
590     }
591     for (int iRem = nSys; iRem < beamB.size(); ++iRem) { 
592       pxBeams     += beamB[iRem].px();
593       pyBeams     += beamB[iRem].py();
594     }
595
596     // Boost all final partons in systems transversely, and also their sum.
597     Vec4 pSumTo( -pxBeams, -pyBeams, pSumOut.pz(), sqrt( pow2(pxBeams) 
598       + pow2(pyBeams) + pow2(pSumOut.pz()) + pSumOut.m2Calc() ) ); 
599     RotBstMatrix Mmismatch;
600     Mmismatch.bst( pSumOut, pSumTo);  
601     for (int iSys = 0; iSys < nSys; ++iSys)  
602     for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
603       int iAB = partonSystemsPtr->getOut(iSys, iMem);
604       if (event[iAB].isFinal()) event[iAB].rotbst(Mmismatch);      
605     }
606     pSumOut.rotbst(Mmismatch);
607
608    // Reset energy and momentum sum, to be compensated by beam remnants.
609     wPosRem = eCM - (pSumOut.e() + pSumOut.pz());
610     wNegRem = eCM - (pSumOut.e() - pSumOut.pz());
611     w2Rem    = wPosRem * wNegRem;
612     if ( wPosRem < 0. || wNegRem < 0. 
613       || sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1])) {
614       infoPtr->errorMsg("Error in BeamRemnants::setKinematics:"
615         " kinematics construction failed owing to recoil mismatch"); 
616       return false;
617     }
618   }
619
620   // Construct x rescaling factors for the two remants.
621   double lambdaRoot = sqrtpos( pow2(w2Rem - w2Beam[0] - w2Beam[1])
622     - 4. * w2Beam[0] * w2Beam[1] );
623   double rescaleA   = (w2Rem + w2Beam[0] - w2Beam[1] + lambdaRoot)
624     / (2. * w2Rem * xSum[0]) ;
625   double rescaleB   = (w2Rem + w2Beam[1] - w2Beam[0] + lambdaRoot)
626     / (2. * w2Rem * xSum[1]) ;
627
628   // Construct energy and pz for remnants in first beam.
629   for (int iRem = nSys; iRem < beamA.size(); ++iRem) {
630     double pPos = rescaleA * beamA[iRem].x() * wPosRem;
631     double pNeg = beamA[iRem].mT2() / pPos;
632     beamA[iRem].e( 0.5 * (pPos + pNeg) );
633     beamA[iRem].pz( 0.5 * (pPos - pNeg) );  
634
635     // Add these partons to the normal event record.
636     int iNew = event.append( beamA[iRem].id(), 63, 1 + nOffset, 0, 0, 0, 
637       beamA[iRem].col(), beamA[iRem].acol(), beamA[iRem].p(), 
638       beamA[iRem].m() );  
639     beamA[iRem].iPos( iNew);
640   }
641
642   // Construct energy and pz for remnants in second beam.
643   for (int iRem = nSys; iRem < beamB.size(); ++iRem) {
644     double pNeg = rescaleB * beamB[iRem].x() * wNegRem;
645     double pPos = beamB[iRem].mT2() / pNeg;
646     beamB[iRem].e( 0.5 * (pPos + pNeg) );
647     beamB[iRem].pz( 0.5 * (pPos - pNeg) );  
648
649     // Add these partons to the normal event record. 
650     int iNew = event.append( beamB[iRem].id(), 63, 2 + nOffset, 0, 0, 0, 
651       beamB[iRem].col(), beamB[iRem].acol(), beamB[iRem].p(), 
652       beamB[iRem].m() );  
653     beamB[iRem].iPos( iNew);
654   }
655
656   // Done.
657   return true;
658
659 }
660
661 //--------------------------------------------------------------------------
662
663 // Allow colour reconnections by mergings of collision subsystems.
664 // iRec is system that may be reconnected, by moving its gluons to iSys,   
665 // where minimal pT (or equivalently Lambda) is used to pick location.
666 // Therefore all dipoles in iSys have to be found, and all gluons in iRec.
667 // Matching q-qbar pairs are treated by analogy with gluons.
668 // Note: owing to rescatterings some outgoing partons must be skipped.
669
670 bool BeamRemnants::reconnectColours( Event&  event) {
671
672   // References to beams to simplify indexing.
673   BeamParticle& beamA = *beamAPtr;  
674   BeamParticle& beamB = *beamBPtr;  
675
676   // Prepare record of which systems should be merged onto another.
677   // The iSys system must have colour in final state to attach to it.
678   vector<int>  iMerge(nSys);
679   vector<bool> hasColour(nSys);
680   for (int iSys = 0; iSys < nSys; ++iSys) {
681     iMerge[iSys] = iSys;
682     bool hasCol = false;
683     for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
684       int iNow = partonSystemsPtr->getOut( iSys, iMem);
685       if (event[iNow].isFinal() && (event[iNow].col() > 0 
686         || event[iNow].acol() > 0) ) {
687         hasCol = true;
688         break;
689       }    
690     }
691     hasColour[iSys] = hasCol;
692   }
693
694   // Loop over systems to decide which should be reconnected. 
695   for (int iRec = nSys - 1; iRec > 0; --iRec) {
696
697     // Determine reconnection strength from pT scale of system.
698     double pT2Rec  = pow2( infoPtr->pTMI(iRec) );
699     double probRec = pT20Rec / (pT20Rec + pT2Rec); 
700
701     // Loop over other systems iSys at higher pT scale and 
702     // decide whether to reconnect the iRec gluons onto one of them.
703     for (int iSys = iRec - 1; iSys >= 0; --iSys)
704     if (hasColour[iSys] && probRec > rndmPtr->flat()) { 
705
706       // The iRec system and all merged with it to be merged with iSys.
707       iMerge[iRec] = iSys;       
708       for (int iRec2 = iRec + 1; iRec2 < nSys; ++iRec2)
709       if (iMerge[iRec2] == iRec) iMerge[iRec2] = iSys;    
710
711       // Once a system has been merged do not test it anymore.
712       break;
713     }
714   }
715
716   // Loop over systems. Check whether other systems to be merged with it.
717   for (int iSys = 0; iSys < nSys; ++iSys) {
718     int nMerge = 0;
719     for (int iRec = iSys + 1; iRec < nSys; ++iRec)
720     if (iMerge[iRec] == iSys) ++nMerge;
721     if (nMerge == 0) continue; 
722
723     // Incoming partons not counted if rescattered.
724     int  iInASys = partonSystemsPtr->getInA(iSys);
725     bool hasInA  = (beamA[iSys].isFromBeam());   
726     int  iInBSys = partonSystemsPtr->getInB(iSys);    
727     bool hasInB  = (beamB[iSys].isFromBeam()); 
728
729     // Begin find dipoles in iSys system.
730     vector<BeamDipole> dipoles;
731     int sizeOut = partonSystemsPtr->sizeOut(iSys);
732     for (int iMem = 0; iMem < sizeOut; ++iMem) {
733
734       // Find colour dipoles to beam remnant.
735       int iNow = partonSystemsPtr->getOut( iSys, iMem);
736       if (!event[iNow].isFinal()) continue;
737       int col = event[iNow].col();  
738       if (col > 0) {
739         if      (hasInA && event[iInASys].col() == col)
740           dipoles.push_back( BeamDipole( col, iNow, iInASys ) );
741         else if (hasInB && event[iInBSys].col() == col)
742           dipoles.push_back( BeamDipole( col, iNow, iInBSys ) );
743  
744         // Find colour dipole between final-state partons.
745         else for (int iMem2 = 0; iMem2 < sizeOut; ++iMem2) 
746         if (iMem2 != iMem) {
747           int iNow2 = partonSystemsPtr->getOut( iSys, iMem2); 
748           if (!event[iNow2].isFinal()) continue;
749           if (event[iNow2].acol() == col) {
750             dipoles.push_back( BeamDipole( col, iNow, iNow2) );
751             break;
752           }
753         }
754       }
755
756       // Find anticolour dipoles to beam remnant.
757       int acol = event[iNow].acol();  
758       if (acol > 0) {
759         if      (hasInA && event[iInASys].acol() == acol)
760           dipoles.push_back( BeamDipole( acol, iInASys, iNow ) );
761         else if (hasInB && event[iInBSys].acol() == acol)
762           dipoles.push_back( BeamDipole( acol, iInBSys, iNow ) ); 
763       }
764     }
765    
766     // Skip mergings if no dipoles found.
767     if (dipoles.size() == 0) continue; 
768
769     // Find dipole sizes.
770     for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) 
771       dipoles[iDip].p1p2 = event[dipoles[iDip].iCol].p() 
772                          * event[dipoles[iDip].iAcol].p();
773     
774     // Loop over systems iRec to be merged with iSys.
775     for (int iRec = iSys + 1; iRec < nSys; ++iRec) {
776       if (iMerge[iRec] != iSys) continue;
777
778       // Information on iRec. Vectors for gluons and anything else.
779       int sizeRec = partonSystemsPtr->sizeOut(iRec);
780       int iInARec = partonSystemsPtr->getInA(iRec);
781       int iInBRec = partonSystemsPtr->getInB(iRec);   
782       int nGluRec = 0;
783       vector<int>    iGluRec;
784       vector<double> pT2GluRec;
785       int nAnyRec = 0;
786       vector<int>    iAnyRec;
787       vector<bool>   freeAnyRec;
788
789       // Copy of gluon positions in descending order. 
790       for (int iMem = 0; iMem < sizeRec; ++iMem) {
791         int iNow = partonSystemsPtr->getOut( iRec, iMem);
792         if (!event[iNow].isFinal()) continue;
793         if (event[iNow].isGluon()) {
794           ++nGluRec;
795           iGluRec.push_back( iNow );  
796           pT2GluRec.push_back( event[iNow].pT2() );
797           for (int i = nGluRec - 1; i > 1; --i) {
798             if (pT2GluRec[i - 1] > pT2GluRec[i]) break;
799             swap(   iGluRec[i - 1],   iGluRec[i] );    
800             swap( pT2GluRec[i - 1], pT2GluRec[i] ); 
801           }  
802         // Copy of anything else, mainly quarks, in no particular order. 
803         } else {
804           ++nAnyRec;
805           iAnyRec.push_back( iNow ); 
806           freeAnyRec.push_back( true ); 
807         }
808       }
809
810       // For each gluon in iRec now find the dipole that gives the smallest
811       // (pGlu * pI) (pGlu * pJ) / (pI * pJ), i.e. minimal pT (and Lambda). 
812       for (int iGRec = 0; iGRec < nGluRec; ++iGRec) {
813         int    iGlu      = iGluRec[iGRec];
814         Vec4   pGlu      = event[iGlu].p();
815         int    iDipMin   = 0;
816         double pT2DipMin = sCM;
817         for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
818           double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p())
819             * (pGlu * event[dipoles[iDip].iAcol].p()) / dipoles[iDip].p1p2;
820           if (pT2Dip < pT2DipMin) {
821             iDipMin   = iDip;
822             pT2DipMin = pT2Dip;
823           }
824         }  
825
826         // Attach the gluon to the dipole, i.e. split the dipole in two.
827         int colGlu   = event[iGlu].col();
828         int acolGlu  = event[iGlu].acol();
829         int colDip   = dipoles[iDipMin].col;
830         int iColDip  = dipoles[iDipMin].iCol;
831         int iAcolDip = dipoles[iDipMin].iAcol;
832         event[iGlu].acol( colDip );
833         if (event[iAcolDip].acol() == colDip) 
834              event[iAcolDip].acol( colGlu );
835         else event[iAcolDip].col(  colGlu ); 
836         dipoles[iDipMin].iAcol = iGlu;
837         dipoles[iDipMin].p1p2 = event[iColDip].p() * pGlu;
838         dipoles.push_back( BeamDipole( colGlu, iGlu, iAcolDip ) );
839         dipoles.back().p1p2 = pGlu * event[iAcolDip].p();
840      
841         // Remove gluon from old system: reconnect colours.
842         for (int i = oldSize; i < event.size(); ++i)
843         if (i != iGlu && i != iAcolDip) { 
844           if (event[i].isFinal()) {     
845             if (event[i].acol() == colGlu) event[i].acol( acolGlu ); 
846           } else {      
847               if (event[i].col()  == colGlu) event[i].col( acolGlu );  
848           }       
849         }
850       }
851
852       // See if any matching quark-antiquark pairs among the rest.
853       for (int iQRec = 0; iQRec < nAnyRec; ++iQRec) {
854         int iQ  = iAnyRec[iQRec];
855         int idQ = event[iQ].id();
856         if (freeAnyRec[iQRec] && idQ > 0 && idQ < 6) 
857         for (int iQbarRec = 0; iQbarRec < nAnyRec; ++iQbarRec) {
858           int iQbar  = iAnyRec[iQbarRec];
859           if (freeAnyRec[iQbarRec] && event[iQbar].id() == -idQ) {
860
861             // Check that these can be traced back to same gluon splitting.
862             // For now also avoid qqbar pairs produced in rescatterings.??
863             int iTopQ    = event.iTopCopyId(iQ);
864             int iTopQbar = event.iTopCopyId(iQbar);
865             int iMother  = event[iTopQ].mother1();
866             if (event[iTopQbar].mother1() == iMother
867               && event[iMother].isGluon() && event[iMother].status() != -34
868               && event[iMother + 1].status() != -34 ) {
869
870               // Now find the dipole that gives the smallest
871               // ((pQ + pQbar) * pI) ((pQ + pQbar) * pJ) / (pI * pJ). 
872               Vec4   pGlu      = event[iQ].p() + event[iQbar].p();
873               int    iDipMin   = 0;
874               double pT2DipMin = sCM;
875               for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
876                 double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p())
877                   * (pGlu * event[dipoles[iDip].iAcol].p()) 
878                   / dipoles[iDip].p1p2;
879                 if (pT2Dip < pT2DipMin) {
880                   iDipMin   = iDip;
881                   pT2DipMin = pT2Dip;
882                 }
883               }  
884
885               // Attach the q-qbar pair to the dipole, i.e. split the dipole.
886               int colGlu   = event[iQ].col();
887               int acolGlu  = event[iQbar].acol();
888               int colDip   = dipoles[iDipMin].col;
889               int iColDip  = dipoles[iDipMin].iCol;
890               int iAcolDip = dipoles[iDipMin].iAcol;
891               event[iQbar].acol( colDip );
892               if (event[iAcolDip].acol() == colDip) 
893                    event[iAcolDip].acol( colGlu );
894               else event[iAcolDip].col(  colGlu ); 
895               dipoles[iDipMin].iAcol = iQbar;
896               dipoles[iDipMin].p1p2 = event[iColDip].p() * event[iQbar].p();
897               dipoles.push_back( BeamDipole( colGlu, iQ, iAcolDip ) );
898               dipoles.back().p1p2 = event[iQ].p() * event[iAcolDip].p();
899      
900               // Remove q-qbar pair from old system: reconnect colours. 
901               freeAnyRec[iQRec]    = false;
902               freeAnyRec[iQbarRec] = false;
903               for (int i = oldSize; i < event.size(); ++i)
904               if (i != iQRec && i != iQbarRec && i != iColDip && i != iAcolDip) { 
905                  if (event[i].isFinal()) {     
906                   if (event[i].acol() == colGlu) event[i].acol( acolGlu ); 
907                 } else {      
908                     if (event[i].col()  == colGlu) event[i].col( acolGlu );  
909                 }       
910               }
911               
912             // Done with processing of q-qbar pairs.
913             }
914           }
915         }
916       }
917
918       // If only two beam gluons left of system, set their colour = anticolour.
919       // Used by BeamParticle::remnantColours to skip irrelevant gluons.
920       if ( event[iInARec].isGluon() && !event[iInARec].isRescatteredIncoming() 
921         && event[iInBRec].isGluon() && !event[iInBRec].isRescatteredIncoming() 
922         && event[iInARec].col() == event[iInBRec].acol() 
923         && event[iInARec].acol() == event[iInBRec].col() ) { 
924           event[iInARec].acol( event[iInARec].col() );
925           event[iInBRec].acol( event[iInBRec].col() );
926       }
927
928     // End of loops over iRec and iSys systems.
929     }
930   }
931
932   // Done.
933   return true;    
934
935 }
936
937 //--------------------------------------------------------------------------
938
939 // Collapse colours and check that they are consistent.
940
941 bool BeamRemnants::checkColours( Event& event) {
942
943   // No colours in lepton beams so no need to do anything.
944   if (beamAPtr->isLepton() && beamBPtr->isLepton()) return true;
945
946   // Remove ambiguities when one colour collapses two ways.
947   // Resolve chains where one colour is mapped to another.
948   for (int iCol = 1; iCol < int(colFrom.size()); ++iCol) 
949   for (int iColRef = 0; iColRef < iCol; ++iColRef) { 
950     if (colFrom[iCol] == colFrom[iColRef]) {
951       colFrom[iCol] = colTo[iCol];
952       colTo[iCol] = colTo[iColRef]; 
953     }
954     if (colTo[iCol] == colFrom[iColRef]) colTo[iCol] = colTo[iColRef];
955   }   
956   
957   // Transform event record colours from beam remnant colour collapses.
958   for (int i = oldSize; i < event.size(); ++i) { 
959     int col = event[i].col();
960     int acol = event[i].acol(); 
961     for (int iCol = 0; iCol < int(colFrom.size()); ++iCol) {
962       if (col == colFrom[iCol]) {col = colTo[iCol]; event[i].col(col);} 
963       if (acol == colFrom[iCol]) {acol = colTo[iCol]; event[i].acol(acol);} 
964     }
965   }
966
967   // Transform junction colours from beam remnant colour collapses.
968   for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
969   for (int leg = 0; leg < 3; ++leg) {
970     int col = event.colJunction(iJun, leg); 
971     for (int iCol = 0; iCol < int(colFrom.size()); ++iCol) 
972     if (col == colFrom[iCol]) {
973       col = colTo[iCol]; 
974       event.colJunction(iJun, leg, col);
975     } 
976   }
977
978   // Arrays for current colours and anticolours, and for singlet gluons.
979   vector<int> colList;
980   vector<int> acolList;
981   vector<int> iSingletGluon;
982
983   // Find current colours and anticolours in the event record.
984   for (int i = oldSize; i < event.size(); ++i) 
985   if (event[i].isFinal()) {
986     int id   = event[i].id();
987     int col  = event[i].col();
988     int acol = event[i].acol(); 
989
990     // Quarks must have colour set, antiquarks anticolour, gluons both.
991     if ( (id > 0 && id < 9 && (col <= 0 || acol != 0) )
992       || (id < 0 && id > -9 && (col != 0 || acol <= 0) )
993       || (id == 21 && (col <= 0 || acol <= 0) ) ) {
994       infoPtr->errorMsg("Error in BeamRemnants::checkColours: "
995         "q/qbar/g has wrong colour slots set");
996       return false;
997     }
998
999     // Save colours/anticolours, and position of colour singlet gluons.
1000     if ( col > 0)  colList.push_back(  col );
1001     if (acol > 0) acolList.push_back( acol );
1002     if (col > 0 && acol == col) iSingletGluon.push_back(i);
1003   }
1004
1005   // Run though list of singlet gluons and put them on final-state dipole
1006   // (i,j) that offers smallest (p_g p_i) * (p_g p_j) / (p_i p_j).
1007   for (int iS = 0; iS < int(iSingletGluon.size()); ++iS) {
1008     int    iGlu      = iSingletGluon[iS];
1009     int    iAcolDip  = -1;
1010     double pT2DipMin = sCM;
1011     for (int iC = oldSize; iC < event.size(); ++iC) 
1012       if (iC != iGlu && event[iC].isFinal()) {
1013       int colDip = event[iC].col();
1014       if (colDip > 0 && event[iC].acol() !=colDip)
1015       for (int iA = oldSize; iA < event.size(); ++iA)
1016         if (iA != iGlu && iA != iC && event[iA].isFinal() 
1017         && event[iA].acol() == colDip && event[iA].col() !=colDip) {
1018         double pT2Dip = (event[iGlu].p() * event[iC].p()) 
1019           * (event[iGlu].p() * event[iA].p())
1020           / (event[iC].p() * event[iA].p());
1021         if (pT2Dip < pT2DipMin) {
1022           iAcolDip  = iA;
1023           pT2DipMin = pT2Dip; 
1024         }
1025       }
1026     }
1027
1028     // Fail if no dipole. Else insert singlet gluon onto relevant dipole.
1029     if (iAcolDip == -1)  return false;
1030     event[iGlu].acol( event[iAcolDip].acol() );
1031     event[iAcolDip].acol( event[iGlu].col() );
1032   }
1033
1034   // Check that not the same colour or anticolour appears twice.
1035   for (int iCol = 0; iCol < int(colList.size()) - 1; ++iCol) {
1036     int col = colList[iCol];
1037     for (int iCol2 = iCol + 1; iCol2 < int(colList.size()); ++iCol2) 
1038     if (colList[iCol2] == col) { 
1039       infoPtr->errorMsg("Warning in BeamRemnants::checkColours:"
1040         " colour appears twice"); 
1041       if (!ALLOWCOLOURTWICE) return false;
1042     }
1043   }
1044   for (int iAcol = 0; iAcol < int(acolList.size()) - 1; ++iAcol) {
1045     int acol = acolList[iAcol];
1046     for (int iAcol2 = iAcol + 1; iAcol2 < int(acolList.size()); ++iAcol2) 
1047     if (acolList[iAcol2] == acol) {
1048       infoPtr->errorMsg("Warning in BeamRemnants::checkColours:"
1049         " anticolour appears twice"); 
1050       if (!ALLOWCOLOURTWICE) return false;
1051     }
1052   }
1053
1054   // Remove all matching colour-anticolour pairs.
1055   bool foundPair = true;
1056   while (foundPair && colList.size() > 0 && acolList.size() > 0) {
1057     foundPair = false;
1058     for (int iCol = 0; iCol < int(colList.size()); ++iCol) {
1059       for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) {
1060         if (acolList[iAcol] == colList[iCol]) { 
1061           colList[iCol] = colList.back(); colList.pop_back();     
1062           acolList[iAcol] = acolList.back(); acolList.pop_back();     
1063           foundPair = true; break;
1064         }
1065       } if (foundPair) break;
1066     }
1067   } 
1068
1069   // Check that remaining (anti)colours are accounted for by junctions.
1070   for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
1071     int kindJun = event.kindJunction(iJun);
1072     for (int leg = 0; leg < 3; ++leg) {
1073       int colEnd = event.colJunction(iJun, leg); 
1074
1075       // Junction connected to three colours.
1076       if (kindJun == 1) {
1077         bool foundCol = false;
1078         for (int iCol = 0; iCol < int(colList.size()); ++iCol) 
1079         if (colList[iCol] == colEnd) { 
1080           colList[iCol] = colList.back(); 
1081           colList.pop_back();     
1082           foundCol = true; 
1083           break;
1084         }  
1085       } 
1086
1087       // Junction connected to three anticolours.
1088       else if (kindJun == 2) {
1089         bool foundCol = false;
1090         for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) 
1091         if (acolList[iAcol] == colEnd) { 
1092           acolList[iAcol] = acolList.back(); 
1093           acolList.pop_back();     
1094           foundCol = true; 
1095           break;
1096         }
1097       }
1098
1099     // End junction check. More junction cases to come??
1100     }
1101   }
1102
1103
1104   // Repair step - sometimes needed when rescattering allowed.
1105   if (colList.size() > 0 || acolList.size() > 0)    
1106     infoPtr->errorMsg("Warning in BeamRemnants::checkColours:"
1107     " need to repair unmatched colours"); 
1108   while (colList.size() > 0 && acolList.size() > 0) {
1109
1110     // Replace one colour and one anticolour index by a new common one.
1111     int  colMatch =  colList.back();
1112     int acolMatch = acolList.back();
1113     int  colNew   = event.nextColTag();
1114     colList.pop_back();     
1115     acolList.pop_back();     
1116     for (int i = oldSize; i < event.size(); ++i) 
1117     if (event[i].isFinal() && event[i].col() == colMatch) {
1118       event[i].col( colNew);
1119       break;
1120     }
1121     for (int i = oldSize; i < event.size(); ++i) 
1122     if (event[i].isFinal() && event[i].acol() == acolMatch) {
1123       event[i].acol( colNew);
1124       break;
1125     }
1126   }
1127
1128   // Done.
1129   return (colList.size() == 0 && acolList.size() == 0);
1130
1131 }
1132
1133 //==========================================================================
1134
1135 } // end namespace Pythia8