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