Line data Source code
1 : // ColosurReconnection.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 : // ColourReconnection class.
8 :
9 : #include "Pythia8/ColourReconnection.h"
10 : namespace Pythia8 {
11 :
12 : //==========================================================================
13 :
14 : // The BeamDipole class is purely internal to reconnectMPIs.
15 :
16 : class BeamDipole {
17 :
18 : public:
19 :
20 : // Constructor.
21 : BeamDipole( int colIn = 0, int iColIn = 0, int iAcolIn = 0)
22 0 : : col(colIn), iCol(iColIn), iAcol(iAcolIn) {}
23 :
24 : // Members.
25 : int col, iCol, iAcol;
26 : double p1p2;
27 :
28 : };
29 :
30 : //==========================================================================
31 :
32 : // The ColourDipole class.
33 :
34 : //--------------------------------------------------------------------------
35 :
36 : // Printing function, inteded for debugging.
37 :
38 : void ColourDipole::print() {
39 :
40 0 : cout << setw(10) << this << setw(6) << col << setw(3) << colReconnection
41 0 : << setw(6) << iCol << setw(5) << iAcol << setw(6) << iColLeg << setw(5)
42 0 : << iAcolLeg << setw(6) << isJun << setw(5) << isAntiJun << setw(10)
43 0 : << p1p2 << " colDips: ";
44 0 : for (int i = 0;i < int(colDips.size());++i)
45 0 : cout << setw(10) << colDips[i];
46 0 : cout << " acolDips: ";
47 0 : for (int i = 0;i < int(acolDips.size());++i)
48 0 : cout << setw(10) << acolDips[i];
49 0 : cout << setw(3) << isActive << endl;
50 :
51 0 : }
52 :
53 : //==========================================================================
54 :
55 : // The InfoGluonMove class is purely internal to reconnectMove.
56 :
57 : class InfoGluonMove{
58 :
59 : public:
60 :
61 : // Constructors.
62 : InfoGluonMove(int i1in, int col1in, int acol1in, int iCol1in, int iAcol1in,
63 : int col2in, int iCol2in, int iAcol2in, double lambdaRefIn,
64 0 : double dLambdaIn) : i1(i1in), i2(0), col1(col1in), acol1(acol1in),
65 0 : iCol1(iCol1in), iAcol1(iAcol1in), col2(col2in), iCol2(iCol2in),
66 0 : iAcol2(iAcol2in), lambdaRef(lambdaRefIn), dLambda(dLambdaIn) {}
67 : InfoGluonMove(int i1in, int i2in, int iCol1in, int iAcol1in, int iCol2in,
68 0 : int iAcol2in, int dLambdaIn) : i1(i1in), i2(i2in), col1(0), acol1(0),
69 0 : iCol1(iCol1in), iAcol1(iAcol1in), col2(0), iCol2(iCol2in),
70 0 : iAcol2(iAcol2in), lambdaRef(0.), dLambda(dLambdaIn) {}
71 :
72 : // Members.
73 : int i1, i2, col1, acol1, iCol1, iAcol1, col2, iCol2, iAcol2;
74 : double lambdaRef, dLambda;
75 :
76 : };
77 :
78 : //==========================================================================
79 :
80 : // The ColourJunction class.
81 :
82 : //--------------------------------------------------------------------------
83 :
84 : // Printing function, inteded for debugging.
85 :
86 : void ColourJunction::print() {
87 :
88 0 : cout << setw(6) << kind() << setw(6)
89 0 : << col(0) << setw(6) << col(1) << setw(6) << col(2) << setw(6)
90 0 : << endCol(0) << setw(6) << endCol(1) << setw(6) << endCol(2) << setw(6)
91 0 : << status(0) << setw(6) << status(1) << setw(6) << status(2) << setw(10)
92 0 : << dips[0] << setw(10) << dips[1] << setw(10) << dips[2] << setw(10)
93 0 : << "\n";
94 0 : cout << " " << setw(10) << dipsOrig[0] << setw(10) << dipsOrig[1]
95 0 : << setw(10) << dipsOrig[2] << endl;
96 :
97 0 : }
98 :
99 : //==========================================================================
100 :
101 : // The ColourParticle class.
102 :
103 : //--------------------------------------------------------------------------
104 :
105 : // Printing function, inteded for debugging.
106 :
107 : void ColourParticle::list() {
108 :
109 0 : const Particle& pt = (*this);
110 :
111 : // Basic line for a particle, always printed.
112 0 : cout << setw(10) << pt.id() << " " << left
113 0 : << setw(18) << pt.nameWithStatus(18) << right << setw(4)
114 0 : << pt.status() << setw(6) << pt.mother1() << setw(6)
115 0 : << pt.mother2() << setw(6) << pt.daughter1() << setw(6)
116 0 : << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol()
117 0 : << setprecision(3)
118 0 : << setw(11) << pt.px() << setw(11) << pt.py() << setw(11)
119 0 : << pt.pz() << setw(11) << pt.e() << setw(11) << pt.m() << "\n";
120 :
121 0 : }
122 :
123 : //--------------------------------------------------------------------------
124 :
125 : // Printing function, inteded for debugging.
126 :
127 : void ColourParticle::listActiveDips() {
128 :
129 0 : cout << "active dips: " << endl;
130 0 : for (int i = 0; i < int(activeDips.size()); ++i)
131 0 : activeDips[i]->print();
132 :
133 0 : }
134 :
135 : //--------------------------------------------------------------------------
136 :
137 : // Printing function, inteded for debugging.
138 :
139 : void ColourParticle::print() {
140 :
141 0 : cout << "--- Particle ---" << endl;
142 0 : for (int i = 0; i < int(dips.size()); ++i) {
143 0 : cout << "(" <<colEndIncluded[i] << ") ";
144 0 : for (int j = 0; j < int(dips[i].size()); ++j) {
145 0 : cout << dips[i][j]->iCol << " (" << dips[i][j]->col << ") ";
146 0 : if (j == int(dips[i].size() - 1))
147 0 : cout << dips[i][j]->iAcol << " (" << acolEndIncluded[i] << ")" << endl;
148 : }
149 : }
150 :
151 0 : }
152 :
153 : //==========================================================================
154 :
155 : // The ColourReconnection class.
156 :
157 : // Minimum needed gain in lambda for a reconnection (to avoid infinity loops).
158 : const double ColourReconnection::MINIMUMGAIN = 1E-10;
159 :
160 : // Minimum needed gain in lambda for a junction.
161 : const double ColourReconnection::MINIMUMGAINJUN = 1E-10;
162 :
163 : // Conversion of GeV^{-1} to fm for time calculations.
164 : const double ColourReconnection::HBAR = 0.197327;
165 :
166 : // Maximum number of reconnection per trial.
167 : // For very large number of outgoing partons, ie. if multiple pp collisions
168 : // are stacked on top of each other, this number needs to be raised.
169 : const int ColourReconnection::MAXRECONNECTIONS = 1000;
170 :
171 : //--------------------------------------------------------------------------
172 :
173 : // Simple comparison function for sort.
174 :
175 : bool cmpTrials(TrialReconnection j1, TrialReconnection j2) {
176 0 : return (j1.lambdaDiff < j2.lambdaDiff);
177 : }
178 :
179 : //--------------------------------------------------------------------------
180 :
181 : // Initialization.
182 :
183 : bool ColourReconnection::init( Info* infoPtrIn, Settings& settings,
184 : Rndm* rndmPtrIn, ParticleData* particleDataPtrIn, BeamParticle* beamAPtrIn,
185 : BeamParticle* beamBPtrIn, PartonSystems* partonSystemsPtrIn) {
186 :
187 : // Save pointers.
188 0 : infoPtr = infoPtrIn;
189 0 : rndmPtr = rndmPtrIn;
190 0 : particleDataPtr = particleDataPtrIn;
191 0 : beamAPtr = beamAPtrIn;
192 0 : beamBPtr = beamBPtrIn;
193 0 : partonSystemsPtr = partonSystemsPtrIn;
194 :
195 : // Total and squared CM energy at nominal energy.
196 0 : eCM = infoPtr->eCM();
197 0 : sCM = eCM * eCM;
198 :
199 : // Choice of reconnection model.
200 0 : reconnectMode = settings.mode("ColourReconnection:mode");
201 :
202 : // pT0 scale of MPI; used in the MPI-based reconnection model.
203 0 : pT0Ref = settings.parm("MultipartonInteractions:pT0Ref");
204 0 : ecmRef = settings.parm("MultipartonInteractions:ecmRef");
205 0 : ecmPow = settings.parm("MultipartonInteractions:ecmPow");
206 0 : pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow);
207 :
208 : // Parameter of the MPI-based reconnection model.
209 0 : reconnectRange = settings.parm("ColourReconnection:range");
210 0 : pT20Rec = pow2(reconnectRange * pT0);
211 :
212 : // Parameters of the new reconnection model.
213 0 : m0 = settings.parm("ColourReconnection:m0");
214 0 : m0sqr = pow2(m0);
215 0 : allowJunctions = settings.flag("ColourReconnection:allowJunctions");
216 0 : nReconCols = settings.mode("ColourReconnection:nColours");
217 0 : sameNeighbourCol = settings.flag("ColourReconnection:sameNeighbourColours");
218 :
219 0 : timeDilationMode = settings.mode("ColourReconnection:timeDilationMode");
220 0 : timeDilationPar = settings.parm("ColourReconnection:timeDilationPar");
221 0 : timeDilationParGeV = timeDilationPar / HBAR;
222 :
223 : // Parameters of gluon-move model.
224 0 : m2Lambda = settings.parm("ColourReconnection:m2Lambda");
225 0 : fracGluon = settings.parm("ColourReconnection:fracGluon");
226 0 : dLambdaCut = settings.parm("ColourReconnection:dLambdaCut");
227 0 : flipMode = settings.mode("ColourReconnection:flipMode");
228 :
229 : // Parameters of the e+e- CR models.
230 0 : singleReconOnly = settings.flag("ColourReconnection:singleReconnection");
231 0 : lowerLambdaOnly = settings.flag("ColourReconnection:lowerLambdaOnly");
232 0 : tfrag = settings.parm("ColourReconnection:fragmentationTime");
233 0 : blowR = settings.parm("ColourReconnection:blowR");
234 0 : blowT = settings.parm("ColourReconnection:blowT");
235 0 : rHadron = settings.parm("ColourReconnection:rHadron");
236 0 : kI = settings.parm("ColourReconnection:kI");
237 :
238 : // Initialize StringLength class.
239 0 : stringLength.init(infoPtr, settings);
240 :
241 : // Done.
242 0 : return true;
243 :
244 0 : }
245 :
246 : //--------------------------------------------------------------------------
247 :
248 : // Do colour reconnection for current event.
249 :
250 : bool ColourReconnection::next( Event& event, int iFirst) {
251 :
252 : // MPI-based reconnection model.
253 0 : if (reconnectMode == 0) return reconnectMPIs(event, iFirst);
254 :
255 : // New reconnection model.
256 0 : else if (reconnectMode == 1) return nextNew(event, iFirst);
257 :
258 : // Gluon-move model.
259 0 : else if (reconnectMode == 2) return reconnectMove(event, iFirst);
260 :
261 : // e+e- Type I CR model.
262 0 : else if (reconnectMode == 3 || reconnectMode == 4)
263 0 : return reconnectTypeCommon(event, iFirst);
264 :
265 : // Undefined.
266 : else {
267 0 : infoPtr->errorMsg("Warning in ColourReconnection::next: "
268 : "Colour reconnecion mode not found");
269 0 : return true;
270 : }
271 :
272 0 : }
273 :
274 : //--------------------------------------------------------------------------
275 :
276 : // Do new colour reconnection for current event.
277 :
278 : bool ColourReconnection::nextNew( Event& event, int iFirst) {
279 :
280 : // Clear old records.
281 0 : while (!dipoles.empty()) {
282 0 : delete dipoles.back();
283 0 : dipoles.pop_back();
284 : }
285 0 : particles.clear();
286 0 : junctions.clear();
287 0 : junTrials.clear();
288 0 : dipTrials.clear();
289 0 : formationTimes.clear();
290 :
291 : // Setup dipoles and make pseudo particles.
292 0 : setupDipoles(event, iFirst);
293 0 : makeAllPseudoParticles(event, iFirst);
294 :
295 : // Setup all dipole reconnections.
296 : // Split dipoles into the 9 different "colours".
297 0 : vector<vector<int> > iDips;
298 0 : iDips.resize(nReconCols);
299 0 : for (int i = 0; i < int(iDips.size()); ++i)
300 0 : iDips[i] = vector<int>();
301 :
302 0 : for (int i = 0; i < int(dipoles.size()); ++i)
303 0 : if (dipoles[i]->isActive)
304 0 : iDips[dipoles[i]->colReconnection].push_back(i);
305 :
306 : // Loop over each colour individually.
307 0 : for (int i = 0;i < int(iDips.size()); ++i)
308 0 : for (int j = 0; j < int(iDips[i].size()); ++j)
309 0 : for (int k = j + 1; k < int(iDips[i].size()); ++k)
310 0 : singleReconnection(dipoles[iDips[i][j]], dipoles[iDips[i][k]]);
311 :
312 : // Only do warning once per event.
313 : bool alreadyWarned = false;
314 :
315 : // Start outer loop over reconnections.
316 0 : for (int iOuterLoop = 0; iOuterLoop < 20; ++iOuterLoop) {
317 : bool finished = true;
318 :
319 : // Do inner loop for string reconnections
320 0 : for (int iInnerLoop = 0;dipTrials.size() > 0; ++iInnerLoop) {
321 :
322 : // Break if too many reonnections are carried out.
323 0 : if (iInnerLoop > MAXRECONNECTIONS) {
324 0 : if (!alreadyWarned)
325 0 : infoPtr->errorMsg("Warning in ColourReconnection::nextNew:"
326 : "Too many reconnections, stopping before minimum reached.");
327 : alreadyWarned = true;
328 0 : break;
329 : }
330 :
331 : // Store all dipoles connected to the chosen dipole.
332 0 : usedDipoles.clear();
333 0 : storeUsedDips(dipTrials.back());
334 :
335 : // Do the reconnection.
336 0 : doDipoleTrial(dipTrials.back());
337 :
338 : // Sort the used dipoles and remove copies of the same.
339 0 : sort(usedDipoles.begin(), usedDipoles.end());
340 0 : for (int i = 0;i < int(usedDipoles.size() - 1); ++i)
341 0 : if (usedDipoles[i] == usedDipoles[i + 1]) {
342 0 : usedDipoles.erase(usedDipoles.begin() + i);
343 0 : i--;
344 0 : }
345 :
346 : // Updating the dipole trials.
347 0 : updateDipoleTrials();
348 : }
349 :
350 : // Loop over list of dipoles to try and form junction structures.
351 0 : if (allowJunctions) {
352 :
353 : // Split dipoles into three categories.
354 0 : iDips.clear();
355 0 : iDips.resize(3);
356 0 : for (int i = 0; i < int(iDips.size()); ++i)
357 0 : iDips[i] = vector<int>();
358 :
359 0 : for (int i = 0; i < int(dipoles.size()); ++i)
360 0 : if (dipoles[i]->isActive)
361 0 : iDips[dipoles[i]->colReconnection % 3].push_back(i);
362 :
363 : // Loop over different "colours" (now only three different groups).
364 0 : for (int i = 0;i < int(iDips.size()); ++i)
365 0 : for (int j = 0; j < int(iDips[i].size()); ++j)
366 0 : for (int k = j + 1; k < int(iDips[i].size()); ++k)
367 0 : singleJunction(dipoles[iDips[i][j]], dipoles[iDips[i][k]]);
368 :
369 :
370 : // Loop over different "colours" (now only three different groups).
371 0 : for (int i = 0;i < int(iDips.size()); ++i)
372 0 : for (int j = 0; j < int(iDips[i].size()); ++j)
373 0 : for (int k = j + 1; k < int(iDips[i].size()); ++k)
374 0 : for (int l = k + 1; l < int(iDips[i].size()); ++l)
375 0 : singleJunction(dipoles[iDips[i][j]], dipoles[iDips[i][k]],
376 0 : dipoles[iDips[i][l]]);
377 :
378 : // Do inner loop for junction reconnections
379 0 : for (int iInnerLoop = 0;junTrials.size() > 0; ++iInnerLoop) {
380 :
381 : // Break if too many reonnections are carried out.
382 0 : if (iInnerLoop > MAXRECONNECTIONS) {
383 0 : if (!alreadyWarned)
384 0 : infoPtr->errorMsg("Warning in ColourReconnection::nextNew:"
385 : "Too many reconnections, stopping before minimum reached.");
386 : alreadyWarned = true;
387 0 : break;
388 : }
389 :
390 : // Find all dipoles connected to the reconnection.
391 0 : usedDipoles.clear();
392 0 : storeUsedDips(junTrials.back());
393 :
394 : // Do the reconnection.
395 0 : doJunctionTrial(event, junTrials.back());
396 :
397 : // Sort the used dipoles and remove copies of the same.
398 0 : sort(usedDipoles.begin(), usedDipoles.end());
399 0 : for (int i = 0;i < int(usedDipoles.size() - 1); ++i)
400 0 : if (usedDipoles[i] == usedDipoles[i + 1]) {
401 0 : usedDipoles.erase(usedDipoles.begin() + i);
402 0 : i--;
403 0 : }
404 :
405 : // Update lists.
406 0 : updateJunctionTrials();
407 0 : updateDipoleTrials();
408 :
409 : finished = false;
410 : }
411 0 : }
412 :
413 : // If no junctions were made, the overall loop is finished.
414 0 : if (finished)
415 0 : break;
416 0 : }
417 :
418 0 : updateEvent(event, iFirst);
419 :
420 : // Done.
421 : return true;
422 0 : }
423 :
424 :
425 : //--------------------------------------------------------------------------
426 :
427 : // Setup initial guess on dipoles, here all colours are assumed
428 : // to be found in the final state.
429 :
430 : void ColourReconnection::setupDipoles( Event& event, int iFirst) {
431 :
432 : // Make vectors needed for storage of chains.
433 0 : vector< vector<int > > chains;
434 0 : vector<bool> isJun;
435 0 : vector<bool> isAntiJun;
436 0 : vector<bool> isGluonLoop;
437 0 : vector<bool> inChain(event.size(),false);
438 :
439 : // Find all quarks and anti-diquarks and follow untill no more colour.
440 0 : for (int i = iFirst; i < event.size(); ++i) {
441 0 : if ( (event[i].isFinal() && !inChain[i]
442 0 : && event[i].isQuark() && event[i].id() > 0)
443 0 : || (event[i].isFinal() && !inChain[i]
444 0 : && event[i].isDiquark() && event[i].id() < 0) ) {
445 0 : int curCol = event[i].col();
446 0 : inChain[i] = true;
447 0 : vector<int> chain;
448 0 : chain.push_back(i);
449 0 : isAntiJun.push_back(false);
450 0 : isJun.push_back(false);
451 0 : isGluonLoop.push_back(false);
452 0 : for (int iSteps = 0; curCol != 0 && iSteps < 1000; ++iSteps) {
453 :
454 : // Check for particles with correct anti colour.
455 0 : for (int j = iFirst; j < event.size(); j++) {
456 0 : if (event[j].isFinal() && !inChain[j] && event[j].acol() == curCol) {
457 0 : chain.push_back(j);
458 0 : inChain[j] = true;
459 0 : curCol = event[j].col();
460 0 : break;
461 : }
462 : }
463 :
464 : // Check for junction with correct colour.
465 0 : for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
466 0 : for (int j = 0; j < 3; ++j) {
467 0 : if (event.colJunction(iJun,j) == curCol) {
468 0 : isJun[isJun.size() -1] = true;
469 : curCol = 0;
470 0 : chain.push_back( -(10 + 10 * iJun + j) );
471 : }
472 : }
473 : }
474 : }
475 0 : chains.push_back(chain);
476 0 : }
477 : }
478 :
479 : // Start from anti-junction and make chains.
480 0 : for (int i = 0; i < event.sizeJunction(); ++i) {
481 :
482 : // First check if junction belongs to the right diffractive system.
483 0 : int checkCol = event.colJunction(i,0);
484 : bool wrongSystem = false;
485 0 : for (int j = 0; j < iFirst; ++j)
486 0 : if (event[j].isFinal() && event[j].acol() == checkCol)
487 0 : wrongSystem = true;
488 0 : if (wrongSystem)
489 0 : continue;
490 :
491 : // Loop over legs of anti junction.
492 0 : if (event.kindJunction(i) == 2)
493 0 : for (int jCol = 0; jCol < 3; ++jCol) {
494 0 : int curCol = event.colJunction(i,jCol);
495 0 : vector<int> chain;
496 0 : chain.push_back( -(10 + 10 * i + jCol));
497 0 : isAntiJun.push_back(true);
498 0 : isJun.push_back(false);
499 0 : isGluonLoop.push_back(false);
500 0 : for (int iSteps = 0; curCol != 0 && iSteps < 1000; ++iSteps) {
501 :
502 : // Check for particles with correct anti colour.
503 0 : for (int j = iFirst; j < event.size(); ++j)
504 0 : if (event[j].isFinal() && !inChain[j] &&
505 0 : event[j].acol() == curCol) {
506 0 : chain.push_back(j);
507 0 : inChain[j] = true;
508 0 : curCol = event[j].col();
509 0 : break;
510 : }
511 :
512 : // Check for junction with correct colour.
513 0 : for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
514 0 : if (event.kindJunction(iJun) == 1)
515 0 : for (int j = 0; j < 3; ++j)
516 0 : if (event.colJunction(iJun,j) == curCol) {
517 0 : isJun[isJun.size() - 1] = true;
518 : curCol = 0;
519 0 : chain.push_back( -(10 + 10 * iJun + j));
520 0 : }
521 : }
522 0 : chains.push_back(chain);
523 0 : }
524 0 : }
525 :
526 : // Find all gluon loops.
527 0 : for (int i = iFirst; i < event.size(); ++i)
528 0 : if (event[i].isFinal() && !inChain[i] && event[i].col() != 0) {
529 0 : int curCol = event[i].col();
530 0 : inChain[i] = true;
531 0 : vector<int> chain;
532 0 : chain.push_back(i);
533 0 : isAntiJun.push_back(false);
534 0 : isJun.push_back(false);
535 0 : isGluonLoop.push_back(true);
536 0 : for (int iSteps = 0; curCol != 0 && iSteps < 1000; ++iSteps) {
537 : bool foundNext = false;
538 0 : for (int j = iFirst; j < event.size(); ++j)
539 0 : if (event[j].isFinal() && !inChain[j] && event[j].acol() == curCol) {
540 0 : chain.push_back(j);
541 0 : inChain[j] = true;
542 0 : curCol = event[j].col();
543 : foundNext = true;
544 0 : break;
545 : }
546 :
547 0 : if (!foundNext)
548 0 : break;
549 0 : }
550 0 : chains.push_back(chain);
551 0 : }
552 :
553 : // Form dipoles from chains.
554 0 : for (int i = 0; i < int(chains.size()); ++i) {
555 0 : if (chains[i].size() == 1) continue;
556 : int lastCol = -1;
557 : int firstCol = 0;
558 :
559 : // Start from the first and form the dipoles.
560 : // Two consececutive dipoles can not share the same colour.
561 0 : for (int j = 0; j < int(chains[i].size()); ++j) {
562 0 : if (j != int(chains[i].size() - 1)) {
563 :
564 : // Start by picking new colour.
565 0 : int newCol = int(rndmPtr->flat() * nReconCols);
566 0 : while (newCol == lastCol && !sameNeighbourCol) {
567 0 : newCol = int(rndmPtr->flat() * nReconCols);
568 : }
569 :
570 : // Need to check whether the quark comes from a g->qqbar split.
571 : // If that is the case, it can not have the same as qbar.
572 0 : if (j == 0 && !isAntiJun[i] && !isGluonLoop[i]) {
573 :
574 0 : int iMother = event[event[ chains[i][j] ].iTopCopy()].mother1();
575 0 : if ( event[iMother].idAbs() == 21) {
576 0 : vector<int> sisters = event[ chains[i][j] ].sisterList(true);
577 : // Need to have only one sister and need to be the anti particle.
578 0 : if (sisters.size() == 1 && event[ chains[i][j] ].id()
579 0 : == - event[ sisters[0] ].id()) {
580 :
581 : // Try to find dipole with sister.
582 : int colSis = -1;
583 0 : for (int k = 0; k < int(dipoles.size()); ++k)
584 0 : if (dipoles[k]->iAcol == sisters[0]) {
585 0 : colSis = dipoles[k]->colReconnection;
586 0 : break;
587 : }
588 :
589 : // If the two colours are the same, pick a new.
590 0 : while (colSis == newCol && !sameNeighbourCol)
591 0 : newCol = int(rndmPtr->flat() * nReconCols);
592 0 : }
593 0 : }
594 0 : }
595 :
596 : // Check if quark end comes from g->qqbar split.
597 : // If so check that the two quarks get different colours.
598 0 : if (j == int(chains[i].size() - 2) && !isJun[i] && !isGluonLoop[i]) {
599 :
600 0 : int iMother = event[event[chains[i][j + 1]].iTopCopy()].mother1();
601 0 : if (event[iMother].idAbs() == 21) {
602 0 : vector<int> sisters = event[ chains[i][j + 1] ].sisterList(true);
603 : // Need to have only one sister and need to be the anti particle.
604 0 : if (sisters.size() == 1 && event[ chains[i][j + 1] ].id()
605 0 : == - event[ sisters[0] ].id()) {
606 :
607 : // Try to find dipole with sister.
608 : int colSis = -1;
609 0 : for (int k = 0; k < int(dipoles.size()); ++k)
610 0 : if (dipoles[k]->iCol == sisters[0]) {
611 0 : colSis = dipoles[k]->colReconnection;
612 0 : break;
613 : }
614 :
615 : // If the two colours are the same, pick a new.
616 0 : while ((colSis == newCol || newCol == lastCol)
617 0 : && !sameNeighbourCol)
618 0 : newCol = int(rndmPtr->flat() * nReconCols);
619 0 : }
620 0 : }
621 0 : }
622 :
623 : // Special case for junction splitting if multiple gluons
624 : // between the junctions.
625 0 : if ((chains[i][j] > 0 && event[chains[i][j]].status() == 75) ||
626 0 : (chains[i][j + 1] > 0 &&
627 0 : event[ chains[i][j + 1] ].status() == 75) ) {
628 :
629 : // Find sisters.
630 0 : vector<int> sisters;
631 0 : if (chains[i][j] > 0 && event[ chains[i][j] ].status() == 75)
632 0 : sisters = event[ chains[i][j] ].sisterList();
633 : else
634 0 : sisters = event[ chains[i][j + 1] ].sisterList();
635 :
636 0 : if (sisters.size() == 3 ) {
637 :
638 : // Find colour of sisters.
639 : int acolSis1 = -1, acolSis2 = -1, acolSis3 = -1;
640 : int colSis1 = -1, colSis2 = -1, colSis3 = -1;
641 0 : for (int k = 0;k < int(dipoles.size()); ++k) {
642 0 : if (dipoles[k]->iAcol == sisters[0])
643 0 : acolSis1 = dipoles[k]->colReconnection;
644 :
645 0 : if (dipoles[k]->iAcol == sisters[1])
646 0 : acolSis2 = dipoles[k]->colReconnection;
647 :
648 0 : if (dipoles[k]->iAcol == sisters[2])
649 0 : acolSis3 = dipoles[k]->colReconnection;
650 :
651 0 : if (dipoles[k]->iCol == sisters[0])
652 0 : colSis1 = dipoles[k]->colReconnection;
653 :
654 0 : if (dipoles[k]->iCol == sisters[1])
655 0 : colSis2 = dipoles[k]->colReconnection;
656 :
657 0 : if (dipoles[k]->iCol == sisters[2])
658 0 : colSis3 = dipoles[k]->colReconnection;
659 : }
660 :
661 : // If the two colours are the same, pick a new.
662 0 : while ((colSis1 == newCol || colSis2 == newCol ||
663 0 : colSis3 == newCol || acolSis1 == newCol ||
664 0 : acolSis2 == newCol || acolSis3 == newCol)
665 0 : && !sameNeighbourCol)
666 0 : newCol = int(rndmPtr->flat() * nReconCols);
667 0 : }
668 0 : }
669 :
670 : // Update stored colours.
671 0 : if (j == 0) firstCol = newCol;
672 : lastCol = newCol;
673 :
674 : // Check if it is anti junction need special dipole.
675 0 : if (j == 0 && isAntiJun[i]) {
676 0 : int col = event.colJunction( - int(chains[i][j]/10) - 1,
677 0 : -chains[i][j] % 10);
678 0 : dipoles.push_back(new ColourDipole(col, chains[i][j],
679 0 : chains[i][j+1], newCol));
680 0 : dipoles.back()->isAntiJun = true;
681 0 : }
682 :
683 : // Otherwise just make the dipole.
684 0 : else dipoles.push_back(new ColourDipole(event[ chains[i][j] ].col(),
685 0 : chains[i][j], chains[i][j+1], newCol));
686 :
687 : // If the chain in end a junction mark it.
688 0 : if (j == int(chains[i].size() - 2) && isJun[i])
689 0 : dipoles.back()->isJun = true;
690 :
691 : // Update the links between dipoles.
692 0 : if (j > 0) {
693 0 : dipoles[dipoles.size() - 1]->leftDip = dipoles[dipoles.size() - 2];
694 0 : dipoles[dipoles.size() - 2]->rightDip = dipoles[dipoles.size() - 1];
695 0 : }
696 0 : }
697 :
698 : // If last particle has anti colour it should be possible to connect it
699 : // to the first particle in the chain. (e.g. gluon loop)
700 : else
701 0 : if (isGluonLoop[i])
702 0 : if (event[ chains[i][j] ].col() == event[ chains[i][0] ].acol()) {
703 0 : int newCol = int(rndmPtr->flat() * nReconCols);
704 0 : while ( (newCol == lastCol || newCol == firstCol)
705 0 : && !sameNeighbourCol) {
706 0 : newCol = int(rndmPtr->flat() * nReconCols);
707 : }
708 0 : dipoles.push_back(new ColourDipole(event[ chains[i][j] ].col(),
709 0 : chains[i][j], chains[i][0], newCol));
710 :
711 : // Update links between dipoles.
712 0 : dipoles[dipoles.size() - 1]->leftDip = dipoles[dipoles.size() - 2];
713 0 : dipoles[dipoles.size() - 2]->rightDip = dipoles[dipoles.size() - 1];
714 0 : dipoles[dipoles.size() - chains[i].size()]->leftDip =
715 0 : dipoles[dipoles.size() -1];
716 0 : dipoles[dipoles.size() - 1]->rightDip =
717 0 : dipoles[dipoles.size() - chains[i].size()];
718 :
719 0 : }
720 : }
721 0 : }
722 :
723 : // Setup junction list.
724 0 : iColJun.clear();
725 0 : iColJun.resize(event.sizeJunction());
726 0 : for (int i = 0; i < int(iColJun.size()); ++i) iColJun[i] = vector<int>(3,0);
727 :
728 : // Loop over event and store indices.
729 0 : for (int i = iFirst; i < event.size(); ++i)
730 0 : if (event[i].isFinal())
731 0 : for (int j = 0; j < event.sizeJunction(); ++j)
732 0 : for (int jLeg = 0; jLeg < 3; ++jLeg)
733 0 : if (event[i].col() == event.colJunction(j,jLeg) ||
734 0 : event[i].acol() == event.colJunction(j,jLeg))
735 0 : iColJun[j][jLeg] = i;
736 :
737 : // Loop over junction and store indices.
738 0 : for (int i = 0;i < event.sizeJunction(); ++i)
739 0 : for (int iLeg = 0; iLeg < 3; ++iLeg)
740 0 : for (int j = i + 1;j < event.sizeJunction(); ++j)
741 0 : for (int jLeg = 0; jLeg < 3; ++jLeg)
742 0 : if (event.colJunction(i, iLeg) == event.colJunction(j, jLeg)) {
743 0 : iColJun[i][iLeg] = -(10*j + 10 + jLeg);
744 0 : iColJun[j][jLeg] = -(10*i + 10 + iLeg);
745 0 : }
746 :
747 : // Setup formation times.
748 0 : setupFormationTimes(event);
749 :
750 : // Done.
751 0 : }
752 :
753 : //--------------------------------------------------------------------------
754 :
755 : // Calculate the string length of a dipole.
756 :
757 : double ColourReconnection::calculateStringLength(ColourDipole * dip,
758 : vector<ColourDipole *> &dips) {
759 :
760 : // Check if dipole was already included.
761 0 : for (int i = 0; i < int(dips.size()); ++i)
762 0 : if (dips[i] == dip) return 0.;
763 :
764 : // Ordinay dipole.
765 0 : if (!dip->isJun && !dip->isAntiJun) {
766 0 : return calculateStringLength(dip->iCol, dip->iAcol);
767 : }
768 : else {
769 :
770 : // Start by finding all particles connected to the junction system.
771 0 : vector<int> iParticles;
772 0 : vector<bool> usedJuns(junctions.size(),false);
773 0 : int nJuns = 0;
774 0 : if (dip->isJun) {
775 0 : if (!findJunctionParticles( -int(dip->iAcol/10) - 1, iParticles,
776 0 : usedJuns, nJuns, dips)) return 1e9;
777 : } else
778 0 : if (!findJunctionParticles(-int(dip->iCol/10) - 1, iParticles,
779 0 : usedJuns, nJuns, dips)) return 1e9;
780 :
781 : // If it is a single junction.
782 0 : if (int(iParticles.size()) == 3)
783 0 : return calculateJunctionLength(iParticles[0], iParticles[1],
784 0 : iParticles[2]);
785 :
786 : // If it is a junction pair.
787 0 : else if (int(iParticles.size()) == 4) {
788 0 : return calculateDoubleJunctionLength(iParticles[0], iParticles[1],
789 0 : iParticles[2], iParticles[3]);
790 : }
791 : // If any other number of junction legs return high number.
792 0 : else return 1e9;
793 0 : }
794 :
795 0 : }
796 : //--------------------------------------------------------------------------
797 :
798 : // Update all colours in the event.
799 :
800 : void ColourReconnection::updateEvent( Event& event, int iFirst) {
801 :
802 : // Start by making a new copy of particles.
803 0 : int oldSize = event.size();
804 0 : for (int i = iFirst; i < oldSize;++i)
805 0 : if (event[i].isFinal()) event.copy(i, 79);
806 :
807 : // Copy over junctions.
808 0 : event.clearJunctions();
809 0 : for (int i = 0; i < int(junctions.size()); ++i) {
810 0 : for (int j = 0; j < 3; ++j) {
811 0 : if ( junctions[i].dipsOrig[j] != 0) {
812 0 : junctions[i].col(j, junctions[i].dipsOrig[j]->col);
813 0 : }
814 : }
815 0 : event.appendJunction(Junction(junctions[i]));
816 : }
817 :
818 : // Assign colour according to the real dipoles.
819 0 : for (int i = 0; i < int(dipoles.size()); ++i)
820 0 : if (dipoles[i]->isReal) {
821 0 : if (dipoles[i]->iCol >= 0)
822 0 : event[ event[ dipoles[i]->iCol ].daughter1() ].col(dipoles[i]->col);
823 : else
824 0 : event.colJunction(-(dipoles[i]->iCol/10 + 1), -dipoles[i]->iCol % 10,
825 0 : dipoles[i]->col);
826 0 : if (dipoles[i]->iAcol >= 0)
827 0 : event[ event[ dipoles[i]->iAcol ].daughter1() ].acol(dipoles[i]->col);
828 : else
829 0 : event.colJunction(-(dipoles[i]->iAcol/10 + 1), -dipoles[i]->iAcol % 10,
830 0 : dipoles[i]->col);
831 : }
832 0 : }
833 :
834 : //--------------------------------------------------------------------------
835 :
836 : // Find all the particles connected in the junction.
837 : // If a single junction, the size of iParticles should be 3.
838 : // For multiple junction structures, the size will increase.
839 :
840 : bool ColourReconnection::findJunctionParticles(int iJun,
841 : vector<int>& iParticles, vector<bool> &usedJuns, int &nJuns,
842 : vector<ColourDipole*> &dips ) {
843 :
844 : // Mark current junction as used.
845 0 : usedJuns[iJun] = true;
846 0 : nJuns++;
847 :
848 : // It is not possible to handle junction structures larger than two.
849 0 : if (nJuns > 2)
850 0 : return false;
851 :
852 : // Find particles connected to the
853 0 : if (junctions[iJun].kind() % 2 == 1)
854 0 : for (int i = 0; i < 3; ++i)
855 0 : iParticles.push_back(junctions[iJun].dips[i]->iCol);
856 : else
857 0 : for (int i = 0; i < 3; ++i)
858 0 : iParticles.push_back(junctions[iJun].dips[i]->iAcol);
859 :
860 : // Add dipoles if not already included.
861 0 : for (int i = 0; i < 3; ++i) {
862 : bool addDip = true;
863 0 : for (int j = 0; j < int(dips.size()); ++j) {
864 0 : if (dips[j] == junctions[iJun].dips[i]) {
865 : addDip = false;
866 0 : break;
867 : }
868 : }
869 0 : if (addDip) dips.push_back(junctions[iJun].dips[i]);
870 : }
871 :
872 : // Check whether it connects to any other junctions.
873 0 : for (int i = 0; i < int(iParticles.size()); ++i)
874 0 : if (iParticles[i] < 0) {
875 0 : int iNewJun = - int(iParticles[i] / 10) -1;
876 0 : iParticles.erase(iParticles.begin() + i);
877 0 : i--;
878 0 : if (!usedJuns[iNewJun] && !findJunctionParticles( iNewJun, iParticles,
879 : usedJuns, nJuns, dips) )
880 0 : return false;
881 0 : }
882 :
883 : // Done.
884 0 : return true;
885 0 : }
886 :
887 : //--------------------------------------------------------------------------
888 :
889 : // Calculate string length for two indices in the particles record.
890 :
891 : double ColourReconnection::calculateStringLength(int i, int j) {
892 0 : return stringLength.getStringLength(particles[i].p(), particles[j].p());
893 : }
894 :
895 : //--------------------------------------------------------------------------
896 :
897 : // Calculate the length of a single junction given the 3 entries in the event.
898 :
899 : double ColourReconnection::calculateJunctionLength(int i,
900 : int j, int k) {
901 :
902 : // Need to be separate indices.
903 0 : if ( i == j || i == k || j == k) return 1e9;
904 :
905 0 : Vec4 p1 = particles[i].p();
906 0 : Vec4 p2 = particles[j].p();
907 0 : Vec4 p3 = particles[k].p();
908 :
909 0 : return stringLength.getJuncLength(p1, p2, p3);
910 :
911 0 : }
912 :
913 : //--------------------------------------------------------------------------
914 :
915 : // Calculate the length of a double junction given the 4 particle entries.
916 : // The first two are expected to be quarks, the second two to be antiquarks.
917 :
918 : double ColourReconnection::calculateDoubleJunctionLength( int i, int j,
919 : int k, int l) {
920 :
921 : // Need to be separate indices.
922 0 : if (i == j || i == k || i == l || j == k || j == l || k == l) return 1e9;
923 :
924 0 : Vec4 p1 = particles[i].p();
925 0 : Vec4 p2 = particles[j].p();
926 0 : Vec4 p3 = particles[k].p();
927 0 : Vec4 p4 = particles[l].p();
928 :
929 0 : return stringLength.getJuncLength(p1, p2, p3, p4);
930 :
931 0 : }
932 :
933 : //--------------------------------------------------------------------------
934 :
935 : // Do a single trial emission.
936 :
937 : void ColourReconnection::singleReconnection(ColourDipole* dip1,
938 : ColourDipole* dip2) {
939 :
940 : // Do nothing if it is the same dipole.
941 0 : if (dip1 == dip2) return;
942 :
943 : // No colour reconnection if colours do not match.
944 0 : if (dip1->colReconnection != dip2->colReconnection) return;
945 :
946 : // If it is not active return
947 0 : if (!dip1->isActive || !dip2->isActive) return;
948 :
949 : // Not possible to connect a gluon with itself.
950 0 : if (dip1->iCol == dip2->iAcol || dip1->iAcol == dip2->iCol) return;
951 :
952 : // Check that reconnection is allowed by time dilation.
953 0 : if (!checkTimeDilation(dip1, dip2)) return;
954 :
955 : // Calculate the difference in lambda.
956 0 : double lambdaDiff = getLambdaDiff(dip1, dip2);
957 :
958 : // Insert into trial reconnection if anything is gained.
959 0 : if (lambdaDiff > MINIMUMGAIN) {
960 0 : TrialReconnection dipTrial(dip1, dip2, 0, 0, 5, lambdaDiff);
961 0 : dipTrials.insert(lower_bound(dipTrials.begin(), dipTrials.end(),
962 : dipTrial, cmpTrials), dipTrial);
963 0 : }
964 :
965 0 : }
966 :
967 : //--------------------------------------------------------------------------
968 :
969 : // Simple test swap between two dipoles.
970 :
971 : void ColourReconnection::swapDipoles(ColourDipole* dip1,
972 : ColourDipole* dip2, bool back) {
973 :
974 : // Swap the anti colour of the dipoles.
975 0 : swap(dip1->iAcol, dip2->iAcol);
976 0 : swap(dip1->isJun, dip2->isJun);
977 0 : swap(dip1->iAcolLeg, dip2->iAcolLeg);
978 :
979 : // Update the active dipoles. Only change 1 active dipole;
980 : // this is to avoid problems when switching back.
981 0 : if (dip1->iAcol != dip2->iAcol) {
982 0 : if (!back) {
983 0 : if (dip1->iAcol >= 0)
984 0 : for (int i = 0; i < int(particles[dip1->iAcol].activeDips.size()); ++i)
985 0 : if (particles[dip1->iAcol].activeDips[i] == dip2) {
986 0 : particles[dip1->iAcol].activeDips[i] = dip1;
987 0 : swap1 = i;
988 0 : break;
989 0 : }
990 0 : if (dip2->iAcol >= 0)
991 0 : for (int i = 0; i < int(particles[dip2->iAcol].activeDips.size()); ++i)
992 0 : if (particles[dip2->iAcol].activeDips[i] == dip1) {
993 0 : particles[dip2->iAcol].activeDips[i] = dip2;
994 0 : swap2 = i;
995 0 : break;
996 0 : }
997 : } else {
998 0 : if (dip1->iAcol >= 0) particles[dip1->iAcol].activeDips[swap2] = dip1;
999 0 : if (dip2->iAcol >= 0) particles[dip2->iAcol].activeDips[swap1] = dip2;
1000 : }
1001 : }
1002 :
1003 : // Update list of junctions (only junctions, anti junctions stay the same).
1004 0 : for (int i = 0; i < int(junctions.size()); ++i)
1005 0 : if (junctions[i].kind() % 2 == 1)
1006 0 : for (int iLeg = 0; iLeg < 3; ++iLeg) {
1007 0 : if (junctions[i].dips[iLeg] == dip1) {
1008 0 : junctions[i].dips[iLeg] = dip2;
1009 0 : continue;
1010 : }
1011 0 : if (junctions[i].dips[iLeg] == dip2) {
1012 0 : junctions[i].dips[iLeg] = dip1;
1013 0 : continue;
1014 : }
1015 0 : }
1016 :
1017 : // Done.
1018 0 : }
1019 :
1020 : //--------------------------------------------------------------------------
1021 :
1022 : // Do a single trial reconnection to form a junction.
1023 :
1024 : void ColourReconnection::singleJunction(ColourDipole* dip1,
1025 : ColourDipole* dip2) {
1026 :
1027 : // Do nothing if it is the same dipole.
1028 0 : if (dip1 == dip2)
1029 : return;
1030 :
1031 0 : int iCol1 = dip1->iCol;
1032 0 : int iCol2 = dip2->iCol;
1033 0 : int iAcol1 = dip1->iAcol;
1034 0 : int iAcol2 = dip2->iAcol;
1035 :
1036 : // Not possible to connect a gluon with itself.
1037 0 : if (iCol1 == iCol2) return;
1038 0 : if (iAcol1 == iAcol2) return;
1039 :
1040 : // Check that all dipoles are active.
1041 0 : if (!dip1->isActive || !dip2->isActive) return;
1042 :
1043 : // Do nothing if one of the dipoles is a junction or anti junction.
1044 0 : if (dip1->isJun || dip1->isAntiJun) return;
1045 0 : if (dip2->isJun || dip2->isAntiJun) return;
1046 :
1047 : // Do nothing if it is a pseudo particle that already contains a
1048 : // baryon number inside of it.
1049 0 : if (int(particles[iCol1].dips.size()) != 1 ||
1050 0 : int(particles[iAcol1].dips.size()) != 1 ||
1051 0 : int(particles[iCol2].dips.size()) != 1 ||
1052 0 : int(particles[iAcol2].dips.size()) != 1)
1053 0 : return;
1054 :
1055 : // Only accept 2/9 of the colour configurations.
1056 0 : if ( (dip1->colReconnection) % 3 !=
1057 0 : dip2->colReconnection % 3) return;
1058 :
1059 0 : if ( (dip1->colReconnection) ==
1060 0 : dip2->colReconnection) return;
1061 :
1062 : // Check that reconnection is allowed by time dilation.
1063 0 : if (!checkTimeDilation(dip1, dip2)) return;
1064 :
1065 : // Find the colour of the last junction leg.
1066 0 : int junCol = (3 - (dip1->colReconnection / 3)
1067 0 : - (dip2->colReconnection / 3) ) * 3
1068 0 : + (dip1->colReconnection % 3);
1069 :
1070 : // if other than 9 colours.
1071 0 : if (nReconCols != 9) {
1072 0 : while (junCol < 0 || junCol % 3 != dip1->colReconnection % 3 ||
1073 0 : junCol == dip1->colReconnection || junCol == dip2->colReconnection)
1074 0 : junCol = int(nReconCols * rndmPtr->flat());
1075 : }
1076 :
1077 : // Store two new dipoles, these will form the anti-junction.
1078 0 : ColourDipole* dip3 = dip1;
1079 0 : ColourDipole* dip4 = dip2;
1080 :
1081 0 : double lambdaDiff = getLambdaDiff(dip1, dip2, dip3, dip4, 0);
1082 0 : if (lambdaDiff > MINIMUMGAINJUN) {
1083 0 : TrialReconnection junTrial(dip1, dip2, dip3, dip4, 0, lambdaDiff);
1084 0 : junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(),
1085 : junTrial, cmpTrials), junTrial);
1086 0 : }
1087 : // Outer loop
1088 : while (true) {
1089 :
1090 : // Reset dip4.
1091 0 : dip4 = dip2;
1092 :
1093 : // If the colour matches that of the junction.
1094 0 : if (dip3->colReconnection == junCol)
1095 : while (true) {
1096 : // Check if the new colour matches.
1097 0 : if (dip4->colReconnection == dip2->colReconnection)
1098 0 : if (checkTimeDilation(dip1, dip2, dip3, dip4)) {
1099 :
1100 : // Calculate lambda measure and store new dipole if anything is gained.
1101 0 : lambdaDiff = getLambdaDiff(dip1, dip2, dip3, dip4, 1);
1102 :
1103 0 : if (lambdaDiff > MINIMUMGAINJUN) {
1104 :
1105 0 : TrialReconnection junTrial(dip1, dip2, dip3, dip4, 1, lambdaDiff);
1106 0 : junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(),
1107 : junTrial, cmpTrials), junTrial);
1108 0 : }
1109 : }
1110 :
1111 :
1112 : // Find the next neighbour.
1113 0 : if (!findAntiNeighbour(dip4))
1114 : break;
1115 :
1116 : // Check for gluon loop.
1117 0 : if (dip4 == dip2 || dip4 == dip1)
1118 : break;
1119 :
1120 : } // Done with inner loop.
1121 :
1122 : // Reset dip4.
1123 0 : dip4 = dip2;
1124 :
1125 : // If the colour matches that of the other dipole.
1126 0 : if (dip3->colReconnection == dip1->colReconnection)
1127 : while (true) {
1128 : // Check if the new colour matches.
1129 0 : if (dip4->colReconnection == junCol)
1130 0 : if (checkTimeDilation(dip1, dip2, dip3, dip4)) {
1131 :
1132 : // Calculate lambda measure and store new dipole if anything is gained.
1133 0 : lambdaDiff = getLambdaDiff(dip1, dip2, dip3, dip4, 2);
1134 :
1135 0 : if (lambdaDiff > MINIMUMGAINJUN) {
1136 :
1137 0 : TrialReconnection junTrial(dip1, dip2, dip3, dip4, 2, lambdaDiff);
1138 0 : junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(),
1139 : junTrial, cmpTrials), junTrial);
1140 0 : }
1141 : }
1142 :
1143 : // Find the next neighbour.
1144 0 : if (!findAntiNeighbour(dip4))
1145 : break;
1146 :
1147 : // Check for gluon loop.
1148 0 : if (dip4 == dip2 || dip4 == dip1)
1149 : break;
1150 :
1151 : } // Done with inner loop.
1152 :
1153 : // Find the next neighbour.
1154 0 : if (!findAntiNeighbour(dip3))
1155 : break;
1156 :
1157 : // Check for gluon loop.
1158 0 : if (dip3 == dip1 || dip3 == dip2)
1159 : break;
1160 : }
1161 :
1162 : // Done.
1163 :
1164 0 : }
1165 :
1166 : //--------------------------------------------------------------------------
1167 :
1168 : // Do a single trial reconnection to form a junction.
1169 :
1170 : void ColourReconnection::singleJunction(ColourDipole* dip1,
1171 : ColourDipole* dip2, ColourDipole* dip3) {
1172 :
1173 : // Do nothing if one of the dipoles is a junction or anti junction.
1174 0 : if (dip1->isJun || dip1->isAntiJun) return;
1175 0 : if (dip2->isJun || dip2->isAntiJun) return;
1176 0 : if (dip3->isJun || dip3->isAntiJun) return;
1177 :
1178 :
1179 : // Check that all dipoles are active.
1180 0 : if (!dip1->isActive || !dip2->isActive || !dip3->isActive) return;
1181 :
1182 : // Only allow 0-3-6, 1-4-7 or 2-5-8.
1183 0 : if ( dip1->colReconnection % 3 != dip2->colReconnection % 3
1184 0 : || dip1->colReconnection % 3 != dip3->colReconnection % 3) return;
1185 :
1186 0 : if ( !(dip1->colReconnection != dip2->colReconnection
1187 0 : && dip1->colReconnection != dip3->colReconnection
1188 0 : && dip2->colReconnection != dip3->colReconnection) )
1189 : return;
1190 :
1191 :
1192 0 : if (int(particles[dip1->iCol].dips.size()) != 1 ||
1193 0 : int(particles[dip1->iAcol].dips.size()) != 1 ||
1194 0 : int(particles[dip2->iCol].dips.size()) != 1 ||
1195 0 : int(particles[dip2->iAcol].dips.size()) != 1 ||
1196 0 : int(particles[dip3->iCol].dips.size()) != 1 ||
1197 0 : int(particles[dip3->iAcol].dips.size()) != 1 )
1198 : return;
1199 :
1200 : // Check that reconnection is allowed by time dilation.
1201 0 : if (!checkTimeDilation(dip1, dip2, dip3)) return;
1202 :
1203 0 : double lambdaDiff = getLambdaDiff(dip1, dip2, dip3, 0, 3);
1204 :
1205 0 : if (lambdaDiff > MINIMUMGAINJUN) {
1206 0 : TrialReconnection junTrial(dip1, dip2, dip3, 0, 3, lambdaDiff);
1207 0 : junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(), junTrial,
1208 : cmpTrials), junTrial);
1209 0 : }
1210 :
1211 : // Done.
1212 : return;
1213 0 : }
1214 :
1215 : // ------------------------------------------------------------------
1216 :
1217 : // Form pseuparticle of a given dipole (or junction system).
1218 :
1219 : void ColourReconnection::makePseudoParticle(ColourDipole* dip , int status,
1220 : bool setupDone) {
1221 :
1222 : // If it is a normal dipole that needs to be combined.
1223 0 : if (!dip->isJun && !dip->isAntiJun) {
1224 :
1225 : // Start by storing variables for easier use.
1226 0 : int iCol = dip->iCol;
1227 0 : int iAcol = dip->iAcol;
1228 0 : int iColLeg = dip->iColLeg;
1229 0 : int iAcolLeg = dip->iAcolLeg;
1230 :
1231 : // Make new pseudo particle.
1232 0 : int iNew = particles.size();
1233 0 : particles.push_back(particles[iCol]);
1234 0 : particles[iNew].acol(particles[iCol].acol());
1235 0 : particles[iNew].col(particles[iAcol].col());
1236 0 : particles[iNew].mother1(iCol);
1237 0 : particles[iNew].mother2(iAcol);
1238 0 : particles[iNew].id(99);
1239 0 : particles[iNew].status(status);
1240 0 : particles[iNew].isJun = false;
1241 0 : particles[iAcol].statusNeg();
1242 0 : particles[iAcol].daughter1(iNew);
1243 0 : particles[iCol].statusNeg();
1244 0 : particles[iCol].daughter1(iNew);
1245 0 : if (iCol != iAcol)
1246 0 : particles[iNew].p(particles[iCol].p() + particles[iAcol].p());
1247 : else
1248 0 : particles[iNew].p(particles[iCol].p());
1249 :
1250 : // Add all the dipoles from the old pseudo particle.
1251 : // First from particle 1.
1252 0 : particles[iNew].dips = particles[dip->iCol].dips;
1253 0 : particles[iNew].colEndIncluded = particles[dip->iCol].colEndIncluded;
1254 0 : particles[iNew].acolEndIncluded = particles[dip->iCol].acolEndIncluded;
1255 :
1256 : // Then particle 2.
1257 0 : if (iCol != iAcol) {
1258 0 : for (int i = 0; i < int(particles[dip->iAcol].dips.size()); ++i) {
1259 0 : if (i != dip->iAcolLeg) {
1260 : // If it is not the same leg, add as separate vector.
1261 0 : particles[iNew].dips.push_back(particles[dip->iAcol].dips[i]);
1262 0 : particles[iNew].colEndIncluded.push_back(
1263 0 : particles[dip->iAcol].colEndIncluded[i]);
1264 0 : particles[iNew].acolEndIncluded.push_back(
1265 0 : particles[dip->iAcol].acolEndIncluded[i]);
1266 0 : } // If it is the same leg, at at the end of the vector.
1267 : else {
1268 0 : particles[iNew].acolEndIncluded[iColLeg] =
1269 0 : particles[iAcol].acolEndIncluded[i];
1270 0 : particles[iNew].dips[iColLeg].pop_back();
1271 0 : particles[iNew].dips[iColLeg].insert(
1272 0 : particles[iNew].dips[iColLeg].end(),
1273 0 : particles[iAcol].dips[i].begin(), particles[iAcol].dips[i].end() );
1274 : }
1275 : }
1276 0 : }
1277 0 : if (iCol != iAcol) {
1278 : // Update the dipole legs to the new particle.
1279 0 : for (int i = 0; i < int(particles[iAcol].activeDips.size()); ++i) {
1280 0 : if ( particles[iAcol].activeDips[i]->iAcol == iAcol) {
1281 0 : if (particles[iAcol].activeDips[i]->iAcolLeg < iAcolLeg)
1282 0 : particles[iAcol].activeDips[i]->iAcolLeg +=
1283 0 : particles[iCol].dips.size();
1284 0 : else if (particles[iAcol].activeDips[i]->iAcolLeg == iAcolLeg)
1285 0 : particles[iAcol].activeDips[i]->iAcolLeg = iColLeg;
1286 0 : else if (particles[iAcol].activeDips[i]->iAcolLeg > iAcolLeg)
1287 0 : particles[iAcol].activeDips[i]->iAcolLeg +=
1288 0 : particles[iCol].dips.size() - 1;
1289 : }
1290 0 : if (particles[iAcol].activeDips[i]->iCol == iAcol) {
1291 0 : if (particles[iAcol].activeDips[i]->iColLeg < iAcolLeg)
1292 0 : particles[iAcol].activeDips[i]->iColLeg +=
1293 0 : particles[iCol].dips.size();
1294 0 : else if (particles[iAcol].activeDips[i]->iColLeg == iAcolLeg)
1295 0 : particles[iAcol].activeDips[i]->iColLeg = iColLeg;
1296 0 : else if (particles[iAcol].activeDips[i]->iColLeg > iAcolLeg)
1297 0 : particles[iAcol].activeDips[i]->iColLeg +=
1298 0 : particles[iCol].dips.size() - 1;
1299 : }
1300 : }
1301 0 : }
1302 :
1303 : // Update list of active dipoles.
1304 0 : particles[iNew].activeDips.clear();
1305 0 : particles[iNew].activeDips.insert(particles[iNew].activeDips.end(),
1306 0 : particles[iCol].activeDips.begin(), particles[iCol].activeDips.end());
1307 0 : if (iCol != iAcol)
1308 0 : particles[iNew].activeDips.insert(particles[iNew].activeDips.end(),
1309 0 : particles[iAcol].activeDips.begin(), particles[iAcol].activeDips.end());
1310 :
1311 : // Remove the now inactive dipole.
1312 0 : for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i)
1313 0 : if (particles[iNew].activeDips[i] == dip) {
1314 0 : particles[iNew].activeDips.erase(
1315 0 : particles[iNew].activeDips.begin() + i);
1316 0 : i--;
1317 0 : }
1318 :
1319 : // Update the indices in the active dipoles.
1320 0 : for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) {
1321 0 : if (particles[iNew].activeDips[i]->iCol == iAcol)
1322 0 : particles[iNew].activeDips[i]->iCol = iNew;
1323 0 : if (particles[iNew].activeDips[i]->iCol == iCol)
1324 0 : particles[iNew].activeDips[i]->iCol = iNew;
1325 0 : if (particles[iNew].activeDips[i]->iAcol == iAcol)
1326 0 : particles[iNew].activeDips[i]->iAcol = iNew;
1327 0 : if (particles[iNew].activeDips[i]->iAcol == iCol)
1328 0 : particles[iNew].activeDips[i]->iAcol = iNew;
1329 0 : particles[iNew].activeDips[i]->p1p2
1330 0 : = mDip(particles[iNew].activeDips[i]);
1331 : }
1332 :
1333 : // If it is a combination of the same particle,
1334 : // check if any double active dipoles
1335 0 : if (iCol == iAcol)
1336 0 : for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i)
1337 0 : for (int j = i + 1; j < int(particles[iNew].activeDips.size()); ++j)
1338 0 : if (particles[iNew].activeDips[i] == particles[iNew].activeDips[j]) {
1339 0 : particles[iNew].activeDips.erase(particles[iNew].activeDips.begin() + j);
1340 0 : j--;
1341 0 : }
1342 :
1343 : // Add dips changed to used dips.
1344 0 : for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) {
1345 0 : if (particles[iNew].activeDips[i]->iCol >= 0)
1346 0 : usedDipoles.push_back(particles[iNew].activeDips[i]);
1347 : else
1348 0 : for (int j = 0;j < 3; ++j)
1349 0 : usedDipoles.push_back(junctions[-(particles[iNew].
1350 0 : activeDips[i]->iCol / 10 + 1)].getColDip(j));
1351 :
1352 0 : if (particles[iNew].activeDips[i]->iAcol >= 0)
1353 0 : usedDipoles.push_back(particles[iNew].activeDips[i]);
1354 : else
1355 0 : for (int j = 0;j < 3; ++j)
1356 0 : usedDipoles.push_back(junctions[-(particles[iNew].
1357 0 : activeDips[i]->iAcol / 10 + 1)].getColDip(j));
1358 : }
1359 :
1360 : // mark the internal dipole as not active.
1361 0 : dip->isActive = false;
1362 :
1363 : // Done.
1364 : return;
1365 : }
1366 :
1367 : // If both ends are connected to a junction something went wrong!
1368 0 : else if (dip->isJun && dip->isAntiJun) {
1369 : return;
1370 : }
1371 : else {
1372 :
1373 : // Find junction index and first leg to combine.
1374 0 : int iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2;
1375 0 : getJunctionIndices(dip, iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2);
1376 0 : ColourDipole* dip2 = junctions[iJun].dips[junLeg1];
1377 0 : ColourDipole* dip3 = junctions[iJun].dips[junLeg2];
1378 :
1379 : // Add new particle.
1380 0 : int iNew = particles.size();
1381 0 : particles.push_back(ColourParticle( Particle( 99, status, i0, i1, 0, 0, 0,
1382 0 : 0, particles[i0].p() + particles[i1].p() ) ) );
1383 0 : particles[iNew].isJun = true;
1384 0 : particles[iNew].junKind = junctions[iJun].kind();
1385 0 : if (i0 == i1) particles[iNew].p(particles[i0].p());
1386 :
1387 : // Update old particles.
1388 0 : particles[i0].statusNeg();
1389 0 : particles[i0].daughter1(iNew);
1390 0 : particles[i1].statusNeg();
1391 0 : particles[i1].daughter1(iNew);
1392 :
1393 : // Update list of internal dipoles.
1394 0 : particles[iNew].dips.clear();
1395 0 : particles[iNew].dips.insert(particles[iNew].dips.end(),
1396 0 : particles[i0].dips.begin(),particles[i0].dips.end());
1397 0 : if (i0 != i1)
1398 0 : particles[iNew].dips.insert(particles[iNew].dips.end(),
1399 0 : particles[i1].dips.begin(),particles[i1].dips.end());
1400 :
1401 : // Update list of whether colour ending is included.
1402 0 : particles[iNew].colEndIncluded.clear();
1403 0 : particles[iNew].colEndIncluded.insert(
1404 0 : particles[iNew].colEndIncluded.end(),
1405 0 : particles[i0].colEndIncluded.begin(),
1406 0 : particles[i0].colEndIncluded.end() );
1407 0 : if (i0 != i1)
1408 0 : particles[iNew].colEndIncluded.insert(
1409 0 : particles[iNew].colEndIncluded.end(),
1410 0 : particles[i1].colEndIncluded.begin(),
1411 0 : particles[i1].colEndIncluded.end() );
1412 :
1413 : // Update list of whether anti colour ending is included.
1414 0 : particles[iNew].acolEndIncluded.clear();
1415 0 : particles[iNew].acolEndIncluded.insert(
1416 0 : particles[iNew].acolEndIncluded.end(),
1417 0 : particles[i0].acolEndIncluded.begin(),
1418 0 : particles[i0].acolEndIncluded.end() );
1419 0 : if (i0 != i1)
1420 0 : particles[iNew].acolEndIncluded.insert(
1421 0 : particles[iNew].acolEndIncluded.end(),
1422 0 : particles[i1].acolEndIncluded.begin(),
1423 0 : particles[i1].acolEndIncluded.end() );
1424 :
1425 : // Third particle just need to add one to list of dipoles.
1426 0 : if (dip->isJun && i2 >= 0 && i2 != i0 && i2 != i1) {
1427 0 : particles[iNew].dips.push_back(particles[i2].dips[dip3->iColLeg]);
1428 0 : particles[iNew].dips.back().erase(particles[iNew].dips.back().begin(),
1429 0 : particles[iNew].dips.back().end() - 1);
1430 :
1431 0 : }
1432 0 : if (dip->isAntiJun && i2 >= 0 && i2 != i0 && i2 != i1) {
1433 0 : particles[iNew].dips.push_back(particles[i2].dips[dip3->iAcolLeg]);
1434 0 : particles[iNew].dips.back().erase(
1435 0 : particles[iNew].dips.back().begin() + 1,
1436 0 : particles[iNew].dips.back().end() );
1437 0 : }
1438 :
1439 : // Add endings for the third particle.
1440 0 : if (i2 != i0 && i2 != i1) {
1441 0 : particles[iNew].acolEndIncluded.push_back(false);
1442 0 : particles[iNew].colEndIncluded.push_back(false);
1443 0 : }
1444 :
1445 : // Special case if it is J-J connection.
1446 0 : if (i2 < 0) {
1447 0 : particles[iNew].dips.push_back(vector<ColourDipole *>());
1448 :
1449 : // Find the real dipole to add to dipole list.
1450 0 : for (int i = 0; i < int(dipoles.size()); ++i)
1451 0 : if (dipoles[i]->isReal && dipoles[i]->iCol == dip3->iCol &&
1452 0 : dipoles[i]->iAcol == dip3->iAcol)
1453 0 : particles[iNew].dips.back().push_back(dipoles[i]);
1454 :
1455 : // Change ending.
1456 0 : particles[iNew].acolEndIncluded.back() = true;
1457 0 : particles[iNew].colEndIncluded.back() = true;
1458 0 : }
1459 :
1460 : // The endings need to reflect the new junction structure.
1461 0 : if (dip->isJun)
1462 0 : for (int i = 0; i < int(particles[iNew].acolEndIncluded.size()); ++i)
1463 0 : particles[iNew].acolEndIncluded[i] = true;
1464 : else
1465 0 : for (int i = 0; i < int(particles[iNew].colEndIncluded.size()); ++i)
1466 0 : particles[iNew].colEndIncluded[i] = true;
1467 :
1468 : // Update active dipoles, first junction case.
1469 : // Set the now internal dipoles as inactive.
1470 0 : dip->isActive = false;
1471 0 : dip2->isActive = false;
1472 0 : dip3->isActive = true;
1473 :
1474 : // Update the dipole legs to the new particle.
1475 : // Only need to do it for the iAcol particle,
1476 : // since nothing changes for the iCol particle.
1477 0 : if (i0 != i1)
1478 0 : for (int i = 0; i < int(particles[i1].activeDips.size()); ++i) {
1479 0 : if (particles[i1].activeDips[i]->iAcol == i1)
1480 0 : particles[i1].activeDips[i]->iAcolLeg += particles[i0].dips.size();
1481 0 : if (particles[i1].activeDips[i]->iCol == i1)
1482 0 : particles[i1].activeDips[i]->iColLeg += particles[i0].dips.size();
1483 0 : }
1484 :
1485 : // Update list of active dipoles.
1486 0 : particles[iNew].activeDips.clear();
1487 0 : particles[iNew].activeDips.insert(particles[iNew].activeDips.end(),
1488 0 : particles[i0].activeDips.begin(), particles[i0].activeDips.end());
1489 0 : if (i0 != i1)
1490 0 : particles[iNew].activeDips.insert(particles[iNew].activeDips.end(),
1491 0 : particles[i1].activeDips.begin(), particles[i1].activeDips.end());
1492 0 : if (i2 != i0 && i2 != i1)
1493 0 : particles[iNew].activeDips.push_back(dip3);
1494 :
1495 : // Remove the now inactive dipoles.
1496 0 : for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) {
1497 0 : if (particles[iNew].activeDips[i] == dip) {
1498 0 : particles[iNew].activeDips.erase(
1499 0 : particles[iNew].activeDips.begin() + i);
1500 0 : i--;
1501 0 : continue;
1502 : }
1503 0 : if (particles[iNew].activeDips[i] == dip2) {
1504 0 : particles[iNew].activeDips.erase(
1505 0 : particles[iNew].activeDips.begin() + i);
1506 0 : i--;
1507 0 : continue;
1508 : }
1509 : }
1510 :
1511 : // Update the indices in the active dipoles.
1512 0 : for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) {
1513 0 : if (particles[iNew].activeDips[i]->iCol == i1)
1514 0 : particles[iNew].activeDips[i]->iCol = iNew;
1515 0 : if (particles[iNew].activeDips[i]->iCol == i0)
1516 0 : particles[iNew].activeDips[i]->iCol = iNew;
1517 0 : if (particles[iNew].activeDips[i]->iAcol == i1)
1518 0 : particles[iNew].activeDips[i]->iAcol = iNew;
1519 0 : if (particles[iNew].activeDips[i]->iAcol == i0)
1520 0 : particles[iNew].activeDips[i]->iAcol = iNew;
1521 0 : particles[iNew].activeDips[i]->p1p2
1522 0 : = mDip(particles[iNew].activeDips[i]);
1523 : }
1524 :
1525 : // The third dip is no longer connected to a junction.
1526 0 : if (dip->isJun) {
1527 0 : dip3->isJun = false;
1528 0 : dip3->iAcol = iNew;
1529 0 : if (i2 != i0 && i2 != i1)
1530 0 : dip3->iAcolLeg = particles[iNew].dips.size() - 1;
1531 : }
1532 : else {
1533 0 : dip3->isAntiJun = false;
1534 0 : dip3->iCol = iNew;
1535 0 : if (i2 != i0 && i2 != i1)
1536 0 : dip3->iColLeg = particles[iNew].dips.size() - 1;
1537 : }
1538 :
1539 : // Add dips changed to used dips.
1540 0 : for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) {
1541 : bool added = false;
1542 0 : for (int j = 0;j < int(usedDipoles.size()); ++j)
1543 0 : if (particles[iNew].activeDips[i] == usedDipoles[j]) {
1544 : added = true;
1545 0 : break;
1546 : }
1547 0 : if (!added) usedDipoles.push_back(particles[iNew].activeDips[i]);
1548 : }
1549 0 : usedDipoles.push_back(dip);
1550 0 : usedDipoles.push_back(dip2);
1551 0 : usedDipoles.push_back(dip3);
1552 :
1553 :
1554 : // Possible for the new dip to have a low m0.
1555 0 : if (setupDone && mDip(dip3) < m0)
1556 0 : makePseudoParticle(dip3, status, true);
1557 0 : }
1558 :
1559 : // Done.
1560 :
1561 0 : }
1562 :
1563 : // ------------------------------------------------------------------
1564 :
1565 : // Help function to sort dipoles in right order.
1566 :
1567 : bool sortFunc(ColourDipole* a, ColourDipole* b) {
1568 0 : return (a->p1p2 < b->p1p2);
1569 : }
1570 :
1571 : // ------------------------------------------------------------------
1572 :
1573 : // Form all pseudoparticles below m0.
1574 :
1575 : void ColourReconnection::makeAllPseudoParticles( Event & event, int iFirst) {
1576 :
1577 : // Make junctions.
1578 0 : for (int i = 0; i < event.sizeJunction(); ++i)
1579 0 : junctions.push_back(event.getJunction(i));
1580 :
1581 : // Make new copy of all the dipoles.
1582 0 : int oldSize = int(dipoles.size());
1583 0 : for (int i = 0; i < oldSize; ++i) {
1584 0 : dipoles.push_back(new ColourDipole(*dipoles[i]));
1585 0 : dipoles[i + oldSize]->iColLeg = 0;
1586 0 : dipoles[i + oldSize]->iAcolLeg = 0;
1587 0 : dipoles[i]->iColLeg = 0;
1588 0 : dipoles[i]->iAcolLeg = 0;
1589 0 : dipoles[i]->isActive = false;
1590 0 : dipoles[i]->isReal = true;
1591 0 : dipoles[i + oldSize]->isReal = false;
1592 :
1593 : // Store original dipoles connected to junctions.
1594 0 : if (dipoles[i]->iCol < 0) {
1595 0 : junctions[-(dipoles[i]->iCol / 10 + 1)].dipsOrig[(-dipoles[i]->iCol)
1596 0 : % 10] = dipoles[i];
1597 0 : }
1598 0 : if (dipoles[i]->iAcol < 0) {
1599 0 : junctions[-(dipoles[i]->iAcol / 10 + 1)].dipsOrig[-(dipoles[i]->iAcol
1600 0 : % 10)] = dipoles[i];
1601 0 : }
1602 : }
1603 :
1604 : // Set up the coldDips and acolDips.
1605 0 : for (int i = 0; i < oldSize; ++i) {
1606 0 : if (dipoles[i]->leftDip != 0)
1607 0 : for (int j = 0; j < oldSize; ++j)
1608 0 : if (dipoles[i]->leftDip == dipoles[j]) {
1609 0 : dipoles[i + oldSize]->colDips.push_back(dipoles[j + oldSize]);
1610 0 : break;
1611 0 : }
1612 :
1613 0 : if (dipoles[i]->rightDip != 0)
1614 0 : for (int j = 0; j < oldSize; ++j)
1615 0 : if (dipoles[i]->rightDip == dipoles[j]) {
1616 0 : dipoles[i + oldSize]->acolDips.push_back(dipoles[j + oldSize]);
1617 0 : break;
1618 0 : }
1619 : }
1620 :
1621 : // Start by copying event record to make pseudoparticles.
1622 : // The pseudoparticles also need to gain
1623 0 : for (int i = iFirst; i < event.size(); ++i)
1624 0 : if (event[i].isFinal()) {
1625 0 : particles.push_back(ColourParticle(event[i]));
1626 0 : particles.back().dips.resize(1,vector<ColourDipole *>());
1627 :
1628 : // Set up dipoles.
1629 0 : for (int j = 0; j < int(dipoles.size()); ++j) {
1630 0 : if (dipoles[j]->iCol == i) {
1631 0 : if (dipoles[j]->isActive) {
1632 0 : dipoles[j]->iCol = particles.size() - 1;
1633 0 : particles.back().activeDips.push_back(dipoles[j]);
1634 0 : }
1635 0 : else particles.back().dips[0].push_back(dipoles[j]);
1636 : }
1637 :
1638 0 : if (dipoles[j]->iAcol == i) {
1639 0 : if (dipoles[j]->isActive) {
1640 0 : dipoles[j]->iAcol = particles.size() - 1;
1641 0 : particles.back().activeDips.push_back(dipoles[j]);
1642 0 : }
1643 0 : else particles.back().dips[0].insert(particles.back().dips[0].begin(),
1644 0 : dipoles[j]);
1645 : }
1646 : }
1647 :
1648 : // Tell whether dipoles are connected to other dipoles.
1649 0 : if (event[i].isQuark() && event[i].id() > 0)
1650 0 : particles.back().colEndIncluded.push_back(true);
1651 0 : else particles.back().colEndIncluded.push_back(false);
1652 :
1653 0 : if (event[i].isQuark() && event[i].id() < 0)
1654 0 : particles.back().acolEndIncluded.push_back(true);
1655 0 : else particles.back().acolEndIncluded.push_back(false);
1656 : }
1657 :
1658 : // Inserting a copy of the event record, but now with full
1659 : // pseudo particle setup.
1660 : // This is mainly to avoid having to distinguish between combining
1661 : // original particles and pseudoparticles.
1662 :
1663 : // Set right dipole connections in junctions.
1664 0 : for (int i = 0; i < int(dipoles.size()); ++i) {
1665 0 : if (dipoles[i]->iCol < 0) {
1666 0 : int j = (- dipoles[i]->iCol / 10) - 1;
1667 0 : int jLeg = - dipoles[i]->iCol % 10;
1668 0 : junctions[j].setColDip(jLeg, dipoles[i]);
1669 0 : }
1670 0 : if (dipoles[i]->iAcol < 0) {
1671 0 : int j = (- dipoles[i]->iAcol / 10) - 1;
1672 0 : int jLeg = - dipoles[i]->iAcol % 10;
1673 0 : junctions[j].setColDip(jLeg, dipoles[i]);
1674 0 : }
1675 : }
1676 :
1677 : // Make sure all dipoles masses are set correctly.
1678 0 : for (int i = 0; i < int(dipoles.size()); ++i) {
1679 0 : if (dipoles[i]->isActive)
1680 0 : dipoles[i]->p1p2 = mDip(dipoles[i]);
1681 : else
1682 0 : dipoles[i]->p1p2 = 1e9;
1683 : }
1684 :
1685 : // Keep making pseudo particles until they are above the threshold.
1686 0 : while (true) {
1687 0 : sort(dipoles.begin(), dipoles.end(), sortFunc);
1688 : bool finished = true;
1689 0 : for (int i = 0; i < int(dipoles.size()); ++i) {
1690 0 : if (!dipoles[i]->isActive) continue;
1691 0 : if (dipoles[i]->p1p2 < m0) {
1692 0 : makePseudoParticle( dipoles[i], 110);
1693 : finished = false;
1694 0 : break;
1695 : }
1696 0 : else break;
1697 : }
1698 0 : if (finished) break;
1699 0 : }
1700 :
1701 : // Sort the dipoles.
1702 0 : sort(dipoles.begin(), dipoles.end(), sortFunc);
1703 :
1704 : // Done.
1705 : return;
1706 :
1707 0 : }
1708 :
1709 : // ------------------------------------------------------------------
1710 :
1711 : // Print statements if something is wrong in dipole setup.
1712 : // Does not have a return statement -- DEBUG PURPOSE ONLY --.
1713 :
1714 : void ColourReconnection::checkRealDipoles(Event& event, int iFirst) {
1715 0 : vector<int> dipConnections(event.size(),0);
1716 0 : for (int i = 0;i < int(dipoles.size()); ++i)
1717 0 : if (dipoles[i]->isReal) {
1718 0 : if (dipoles[i]->iCol >= 0)
1719 0 : dipConnections[dipoles[i]->iCol]++;
1720 0 : if (dipoles[i]->iAcol >= 0)
1721 0 : dipConnections[dipoles[i]->iAcol]++;
1722 : }
1723 : bool working = true;
1724 0 : for (int i = iFirst ;i < event.size(); ++i) {
1725 0 : if (event[i].isFinal()) {
1726 0 : if (event[i].isQuark() && dipConnections[i] != 1) {
1727 0 : cout << "quark " << i << " is wrong!!" << endl;
1728 : working = false;
1729 0 : }
1730 0 : else if (event[i].idAbs() == 21 && dipConnections[i] != 2) {
1731 0 : cout << "gluon " << i << " is wrong!!" << endl;
1732 : working = false;
1733 0 : }
1734 : }
1735 : }
1736 0 : if (!working) {
1737 0 : infoPtr->errorMsg("Error in ColourReconnection::checkRealDipoles:"
1738 : "Real dipoles not setup properply");
1739 :
1740 0 : }
1741 :
1742 0 : }
1743 : // ------------------------------------------------------------------
1744 :
1745 : // Print statements if something is wrong in dipole setup.
1746 : // Does not have a return statement -- DEBUG PURPOSE ONLY --.
1747 :
1748 : void ColourReconnection::checkDipoles() {
1749 :
1750 0 : for (int i = 0;i < int(dipoles.size()); ++i) {
1751 0 : if (dipoles[i] == 0) { cout << "dipole empty" << endl;}
1752 0 : if (dipoles[i]->isActive) {
1753 0 : if (dipoles[i]->iCol >= 0) {
1754 : bool foundMyself = false;
1755 0 : for (int j = 0; j < int(particles[ dipoles[i]->iCol ].
1756 0 : activeDips.size()); ++j) {
1757 0 : if (!particles[dipoles[i]->iCol].activeDips[j]->isActive) {
1758 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1759 : "Found inactive dipole, where only active was expected");
1760 0 : }
1761 0 : if (particles[dipoles[i]->iCol].activeDips[j] == dipoles[i])
1762 0 : foundMyself = true;
1763 : }
1764 :
1765 0 : if (!foundMyself) {
1766 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1767 : "Linking between active dipoles and particles is wrong");
1768 0 : }
1769 0 : if (dipoles[i]->iColLeg
1770 0 : >= int(particles[dipoles[i]->iCol].dips.size())) {
1771 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1772 : "Original dipoles not stored correct");
1773 0 : }
1774 :
1775 : // Check that linking to old dipoles work.
1776 0 : if (dipoles[i]->col !=
1777 0 : particles[dipoles[i]->iCol].dips[dipoles[i]->iColLeg].back()->col) {
1778 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1779 : "Original dipoles do not match in");
1780 0 : }
1781 0 : }
1782 :
1783 0 : if (dipoles[i]->iAcol >= 0) {
1784 : bool foundMyself = false;
1785 0 : for (int j = 0;j < int(particles[ dipoles[i]->iAcol ].
1786 0 : activeDips.size()); ++j) {
1787 :
1788 0 : if (!particles[dipoles[i]->iAcol].activeDips[j]->isActive) {
1789 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1790 : "Found inactive dipole, where only active was expected");
1791 0 : }
1792 0 : if (particles[dipoles[i]->iAcol].activeDips[j] == dipoles[i])
1793 0 : foundMyself = true;
1794 : }
1795 :
1796 0 : if (!foundMyself) {
1797 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1798 : "Linking between active dipoles and particles is wrong");
1799 0 : }
1800 0 : if (dipoles[i]->iAcolLeg >= int(particles[dipoles[i]->iAcol].
1801 0 : dips.size() )) {
1802 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1803 : "Original dipoles not stored correct");
1804 0 : }
1805 :
1806 : // Check that linking to old dipoles work
1807 0 : if (dipoles[i]->col != particles[dipoles[i]->iAcol].
1808 0 : dips[dipoles[i]->iAcolLeg].front()->col) {
1809 0 : infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:"
1810 : "Original dipoles do not match in");
1811 0 : }
1812 0 : }
1813 : }
1814 : }
1815 0 : }
1816 :
1817 : // ------------------------------------------------------------------
1818 :
1819 : // Print all the chains.
1820 :
1821 : void ColourReconnection::listAllChains() {
1822 :
1823 0 : cout << " ----- PRINTING CHAINS ----- " << dipoles.size() << endl;
1824 0 : for (int i = 0; i < int(dipoles.size()); ++i)
1825 0 : dipoles[i]->printed = false;
1826 :
1827 0 : for (int i = 0;i < int(dipoles.size()); ++i)
1828 0 : if (!dipoles[i]->printed)
1829 0 : listChain(dipoles[i]);
1830 0 : cout << " ----- PRINTED CHAINS ----- " << endl;
1831 :
1832 0 : }
1833 :
1834 : // ------------------------------------------------------------------
1835 :
1836 : // Print the chain containing the dipole.
1837 :
1838 : void ColourReconnection::listChain(ColourDipole *dip) {
1839 :
1840 : // Make sure not an empty pointer.
1841 0 : if (dip == 0) return;
1842 :
1843 : // If chain is not active, just print it.
1844 0 : if (!dip->isActive) {
1845 : return;
1846 : }
1847 :
1848 0 : ColourDipole * colDip = dip;
1849 :
1850 : // Try to reach one end of the chain.
1851 0 : while (particles[colDip->iCol].dips.size() == 1 && findColNeighbour(colDip))
1852 0 : if (dip == colDip)
1853 : break;
1854 :
1855 0 : ColourDipole * endDip = colDip;
1856 0 : do {
1857 0 : cout << colDip->iCol << " (" << colDip->p1p2 << ", " << colDip->col
1858 0 : << ") (" << colDip->isActive << ") ";
1859 0 : colDip->printed = true;
1860 0 : }
1861 : // Start the printing.
1862 0 : while (particles[colDip->iAcol].dips.size() == 1 && findAntiNeighbour(colDip)
1863 0 : && colDip != endDip);
1864 :
1865 : // Print the last part.
1866 0 : cout << colDip->iAcol<< endl;
1867 :
1868 : // Done.
1869 0 : }
1870 :
1871 : // ------------------------------------------------------------------
1872 :
1873 : // Return relevant indices for the junction.
1874 :
1875 : bool ColourReconnection::getJunctionIndices(ColourDipole * dip, int &iJun,
1876 : int &i0, int &i1, int &i2, int &junLeg0, int &junLeg1, int &junLeg2) {
1877 :
1878 : // Find junction index and first leg to combine.
1879 0 : int indxJun = dip->iCol;
1880 0 : if (dip->iAcol < 0)
1881 0 : indxJun = dip->iAcol;
1882 0 : iJun = (- indxJun / 10) - 1;
1883 0 : junLeg0 = -(indxJun % 10);
1884 0 : junLeg1 = 1;
1885 0 : junLeg2 = 2;
1886 0 : if (junLeg0 == 1) junLeg1 = 0;
1887 0 : else if (junLeg0 == 2) junLeg2 = 0;
1888 :
1889 0 : if (dip->iCol < 0) {
1890 0 : i0 = dip->iAcol;
1891 0 : i1 = junctions[iJun].dips[junLeg1]->iAcol;
1892 0 : i2 = junctions[iJun].dips[junLeg2]->iAcol;
1893 0 : }
1894 : else {
1895 0 : i0 = dip->iCol;
1896 0 : i1 = junctions[iJun].dips[junLeg1]->iCol;
1897 0 : i2 = junctions[iJun].dips[junLeg2]->iCol;
1898 : }
1899 :
1900 : // It is not possible to form a pseudoparticle if only a single particle is
1901 : // connected to the junction.
1902 0 : if (i1 < 0 && i2 < 0) return false;
1903 :
1904 : // Check which two particle should form the pseudoparticle.
1905 : double m1 = 1e9, m2 = 1e9;
1906 0 : if (i1 >= 0)
1907 0 : m1 = m(particles[i0].p(),particles[i1].p());
1908 0 : if (i2 >= 0)
1909 0 : m2 = m(particles[i0].p(),particles[i2].p());
1910 :
1911 0 : if (m1 > m2) {
1912 0 : swap(i1,i2);
1913 0 : swap(junLeg1,junLeg2);
1914 0 : }
1915 : // Force switch if i0 == i2
1916 0 : if (i0 == i2) {
1917 0 : swap(i1,i2);
1918 0 : swap(junLeg1,junLeg2);
1919 0 : }
1920 :
1921 : return true;
1922 0 : }
1923 :
1924 :
1925 : // ------------------------------------------------------------------
1926 :
1927 : // Check whether up to four dipoles are 'causally' connected.
1928 :
1929 : bool ColourReconnection::checkTimeDilation(ColourDipole * dip1,
1930 : ColourDipole * dip2, ColourDipole * dip3, ColourDipole * dip4) {
1931 :
1932 0 : if (timeDilationMode == 0) return true;
1933 :
1934 : // 2 dipole case.
1935 0 : if (dip3 == 0) {
1936 0 : Vec4 p1 = getDipoleMomentum(dip1);
1937 0 : Vec4 p2 = getDipoleMomentum(dip2);
1938 0 : double t1 = formationTimes[dip1->col];
1939 0 : double t2 = formationTimes[dip2->col];
1940 0 : if (dip1 == dip2) return true;
1941 0 : else return checkTimeDilation(p1, p2, t1, t2);
1942 :
1943 : // 3 dipole case.
1944 0 : } else if (dip4 == 0) {
1945 0 : Vec4 p1 = getDipoleMomentum(dip1);
1946 0 : Vec4 p2 = getDipoleMomentum(dip2);
1947 0 : Vec4 p3 = getDipoleMomentum(dip3);
1948 0 : double t1 = formationTimes[dip1->col];
1949 0 : double t2 = formationTimes[dip2->col];
1950 0 : double t3 = formationTimes[dip3->col];
1951 : // Modes that require all dipoles to be causally connected.
1952 0 : if (timeDilationMode == 1 || timeDilationMode == 2 ||
1953 0 : timeDilationMode == 4) {
1954 0 : if (dip1 != dip2 && !checkTimeDilation(p1, p2, t1, t2)) return false;
1955 0 : if (dip1 != dip3 && !checkTimeDilation(p1, p3, t1, t3)) return false;
1956 0 : if (dip2 != dip3 && !checkTimeDilation(p2, p3, t2, t3)) return false;
1957 0 : return true;
1958 : // Modes that require a single pair of dipoles to be causally connected.
1959 : } else {
1960 0 : if (dip1 != dip2 && checkTimeDilation(p1, p2, t1, t2)) return true;
1961 0 : if (dip1 != dip3 && checkTimeDilation(p1, p3, t1, t3)) return true;
1962 0 : if (dip2 != dip3 && checkTimeDilation(p2, p3, t2, t3)) return true;
1963 0 : return false;
1964 : }
1965 :
1966 : // 4 dipole case.
1967 0 : } else {
1968 0 : Vec4 p1 = getDipoleMomentum(dip1);
1969 0 : Vec4 p2 = getDipoleMomentum(dip2);
1970 0 : Vec4 p3 = getDipoleMomentum(dip3);
1971 0 : Vec4 p4 = getDipoleMomentum(dip4);
1972 0 : double t1 = formationTimes[dip1->col];
1973 0 : double t2 = formationTimes[dip2->col];
1974 0 : double t3 = formationTimes[dip3->col];
1975 0 : double t4 = formationTimes[dip4->col];
1976 : // Modes that require all dipoles to be causally connected.
1977 0 : if (timeDilationMode == 1 || timeDilationMode == 2 ||
1978 0 : timeDilationMode == 4) {
1979 0 : if (dip1 != dip2 && !checkTimeDilation(p1, p2, t1, t2)) return false;
1980 0 : if (dip1 != dip3 && !checkTimeDilation(p1, p3, t1, t3)) return false;
1981 0 : if (dip1 != dip4 && !checkTimeDilation(p1, p4, t1, t4)) return false;
1982 0 : if (dip2 != dip3 && !checkTimeDilation(p2, p3, t2, t3)) return false;
1983 0 : if (dip2 != dip4 && !checkTimeDilation(p2, p4, t2, t4)) return false;
1984 0 : if (dip3 != dip4 && !checkTimeDilation(p3, p4, t3, t4)) return false;
1985 0 : return true;
1986 : // Modes that require a single pair of dipoles to be causally connected.
1987 : } else {
1988 0 : if (dip1 != dip2 && checkTimeDilation(p1, p2, t1, t2)) return true;
1989 0 : if (dip1 != dip3 && checkTimeDilation(p1, p3, t1, t3)) return true;
1990 0 : if (dip1 != dip4 && checkTimeDilation(p1, p4, t1, t4)) return true;
1991 0 : if (dip2 != dip3 && checkTimeDilation(p2, p3, t2, t3)) return true;
1992 0 : if (dip2 != dip4 && checkTimeDilation(p2, p4, t2, t4)) return true;
1993 0 : if (dip3 != dip4 && checkTimeDilation(p3, p4, t3, t4)) return true;
1994 0 : return false;
1995 : }
1996 0 : }
1997 :
1998 : // Done.
1999 0 : }
2000 :
2001 : // ------------------------------------------------------------------
2002 :
2003 : // Find the momentum of the dipole.
2004 :
2005 : Vec4 ColourReconnection::getDipoleMomentum(ColourDipole * dip) {
2006 0 : vector<int> iPar, usedJuncs;
2007 0 : if (!dip->isJun) iPar.push_back(dip->iAcol);
2008 0 : else addJunctionIndices(dip->iAcol, iPar, usedJuncs);
2009 0 : if (!dip->isAntiJun) iPar.push_back(dip->iCol);
2010 0 : else addJunctionIndices(dip->iCol, iPar, usedJuncs);
2011 :
2012 : // Remove any duplicates.
2013 0 : sort(iPar.begin(),iPar.end());
2014 0 : for (int i = 0;i < int(iPar.size()) - 1; ++i)
2015 0 : if (iPar[i] == iPar[i+1]) {
2016 0 : iPar.erase(iPar.begin()+i);
2017 0 : i--;
2018 0 : }
2019 :
2020 0 : if (iPar.size() == 0) {
2021 0 : infoPtr->errorMsg("Error in ColourReconnection::getDipoleMomentum: "
2022 : "No particles connected to junction.");
2023 0 : return Vec4(0,0,0,0);
2024 : }
2025 :
2026 0 : Vec4 p = particles[iPar[0]].p();
2027 0 : for (int i = 1;i < int(iPar.size()); ++i)
2028 0 : p += particles[iPar[i]].p();
2029 :
2030 0 : return p;
2031 0 : }
2032 :
2033 : // ------------------------------------------------------------------
2034 :
2035 : // Check whether two four momenta are 'causally' connected.
2036 :
2037 : bool ColourReconnection::checkTimeDilation(Vec4 p1,
2038 : Vec4 p2, double t1, double t2) {
2039 : // No time dilation check.
2040 0 : if (timeDilationMode == 0) return true;
2041 :
2042 : // Check if gamma is below parameter.
2043 0 : else if (timeDilationMode == 1) {
2044 0 : p2.bstback(p1);
2045 0 : if (p2.e() / p2.mCalc() > timeDilationPar) return false;
2046 0 : else return true;
2047 :
2048 : // Check if gamma * mDip is below parameter for both dipoles.
2049 0 : } else if (timeDilationMode == 2) {
2050 : bool part1, part2;
2051 0 : p2.bstback(p1);
2052 0 : if (p2.e() / p2.mCalc() > timeDilationParGeV * p2.mCalc()) part1 = false;
2053 : else part1 = true;
2054 0 : p2.bst(p1);
2055 0 : p1.bstback(p2);
2056 0 : if (p1.e() / p1.mCalc() > timeDilationParGeV * p1.mCalc()) part2 = false;
2057 : else part2 = true;
2058 0 : if (part1 && part2) return true;
2059 0 : else return false;
2060 :
2061 : // Check if gamma * mDip is below parameter for a single dipole.
2062 0 : } else if (timeDilationMode == 3) {
2063 : bool part1, part2;
2064 0 : p2.bstback(p1);
2065 0 : if (p2.e() / p2.mCalc() > timeDilationParGeV * p2.mCalc()) part1 = false;
2066 : else part1 = true;
2067 0 : p2.bst(p1);
2068 0 : p1.bstback(p2);
2069 0 : if (p1.e() / p1.mCalc() > timeDilationParGeV * p1.mCalc()) part2 = false;
2070 : else part2 = true;
2071 0 : if (part1 || part2) return true;
2072 0 : else return false;
2073 :
2074 : // Check if gamma * mDip' is below parameter for both dipoles.
2075 0 : } else if (timeDilationMode == 4) {
2076 0 : p2.bstback(p1);
2077 0 : if (p2.e() / p2.mCalc() < timeDilationParGeV * min(t1,t2)) return true;
2078 0 : else return false;
2079 :
2080 : // Check if gamma * mDip' is below parameter for a single dipole.
2081 0 : } else if (timeDilationMode == 5) {
2082 0 : p2.bstback(p1);
2083 0 : if (p2.e() / p2.mCalc() < timeDilationParGeV * max(t1,t2)) return true;
2084 0 : else return false;
2085 :
2086 : // If mode is set wrong, should never happen.
2087 0 : } else return true;
2088 0 : }
2089 :
2090 : // ------------------------------------------------------------------
2091 :
2092 : // Store the formation times.
2093 :
2094 : void ColourReconnection::setupFormationTimes(Event & event) {
2095 :
2096 0 : for (int i = 0;i < event.size(); ++i) {
2097 : // First check the colour.
2098 0 : if (event[i].col() != 0 && formationTimes.count(event[i].col()) == 0) {
2099 0 : int col = event[i].col();
2100 : // Find first time the colour appears as an anticolour.
2101 : bool foundCol = false;
2102 : int iAcol = 0;
2103 0 : for (int j = i;j < event.size(); ++j) {
2104 0 : if (event[j].acol() == col) {
2105 : foundCol = true;
2106 : iAcol = j;
2107 0 : break;
2108 : }
2109 : }
2110 :
2111 : // If it was found add it to the list.
2112 0 : if (foundCol) {
2113 0 : formationTimes[col] = max( m0,
2114 0 : (event[i].p() + event[iAcol].p()).mCalc() );
2115 : // Otherwise it must be stored in a junction.
2116 0 : } else {
2117 0 : formationTimes[col] = max(m0, getJunctionMass(event, col));
2118 : }
2119 0 : }
2120 :
2121 : // Next check the anti colour.
2122 0 : if (event[i].acol() != 0 && formationTimes.count(event[i].acol()) == 0) {
2123 0 : int acol = event[i].acol();
2124 : // Find first time the colour appears as an anticolour.
2125 : bool foundCol = false;
2126 : int iCol = 0;
2127 0 : for (int j = i;j < event.size(); ++j) {
2128 0 : if (event[j].col() == acol) {
2129 : foundCol = true;
2130 : iCol = j;
2131 0 : break;
2132 : }
2133 : }
2134 :
2135 : // If it was found add it to the list.
2136 0 : if (foundCol) {
2137 0 : formationTimes[acol] = max(m0, (event[i].p() + event[iCol].p())
2138 : .mCalc());
2139 : // Otherwise it must be stored in a junction.
2140 0 : } else {
2141 0 : formationTimes[acol] = max(m0, getJunctionMass(event, acol));
2142 : }
2143 0 : }
2144 : }
2145 :
2146 : // Finally check if junction colours are stored.
2147 0 : for (int i = 0; i < event.sizeJunction(); ++i)
2148 0 : for (int j = 0; j < 3; ++j)
2149 0 : if (formationTimes.count(event.colJunction(i,j)) == 0)
2150 0 : formationTimes[event.colJunction(i,j)] = max(m0,
2151 0 : getJunctionMass(event, event.colJunction(i,j)));
2152 :
2153 : // Done.
2154 0 : }
2155 :
2156 : // ------------------------------------------------------------------
2157 :
2158 : // Find the invariant mass of all the partons connected to a junction system.
2159 :
2160 : double ColourReconnection::getJunctionMass(Event & event, int col) {
2161 :
2162 : // Find the partons connected to the junction system.
2163 0 : vector<int> iPar, usedJuncs;
2164 0 : addJunctionIndices(event, col, iPar, usedJuncs);
2165 :
2166 : // Check for doubles.
2167 0 : sort(iPar.begin(), iPar.end());
2168 0 : for (int i = 0;i < int(iPar.size() -1); ++i) {
2169 0 : if (iPar[i] == iPar[i + 1]) {
2170 0 : iPar.erase(iPar.begin() + i);
2171 0 : i--;
2172 0 : }
2173 : }
2174 :
2175 : // If no partons are connected to the junction system
2176 : // (or it was not a junction system).
2177 0 : if (int(iPar.size()) == 0) return 0;
2178 :
2179 0 : Vec4 p = event[iPar[0]].p();
2180 0 : for (int i = 1; i < int(iPar.size()); ++i)
2181 0 : p += event[iPar[i]].p();
2182 :
2183 0 : return p.mCalc();
2184 :
2185 0 : }
2186 :
2187 : // ------------------------------------------------------------------
2188 :
2189 : // Find all particles connected to a junction system for junctions stored in
2190 : // the event record.
2191 :
2192 : void ColourReconnection::addJunctionIndices(Event & event, int col,
2193 : vector<int> &iPar, vector<int> &usedJuncs) {
2194 :
2195 : // Find the junction.
2196 0 : vector<int> juncs;
2197 0 : for (int j = 0;j < event.sizeJunction(); ++j) {
2198 0 : for (int k = 0; k < 3; ++k) {
2199 0 : if (event.colJunction(j,k) == col) {
2200 0 : juncs.push_back(j);
2201 0 : break;
2202 : }
2203 : }
2204 : }
2205 :
2206 : // Check if they were already used.
2207 0 : for (int i = 0;i < int(juncs.size()); ++i) {
2208 0 : for (int j = 0;j < int(usedJuncs.size()); ++j) {
2209 0 : if (juncs[i] == usedJuncs[j]) {
2210 0 : juncs.erase(juncs.begin() + i);
2211 0 : i--;
2212 0 : break;
2213 : }
2214 : }
2215 : }
2216 :
2217 : // If list of junctions is empty return.
2218 0 : if (juncs.size() == 0) return;
2219 :
2220 : // Store the juncstions as used.
2221 0 : for (int i = 0;i < int(juncs.size()); ++i)
2222 0 : usedJuncs.push_back(juncs[i]);
2223 :
2224 : // Find the partons connected to it.
2225 0 : for (int iJunc = 0; iJunc < int(juncs.size()); ++iJunc) {
2226 0 : int iTempPars[3] = {-1,-1,-1};
2227 0 : int cols[3] = {event.colJunction(juncs[iJunc],0),
2228 0 : event.colJunction(juncs[iJunc],1), event.colJunction(juncs[iJunc],2)};
2229 :
2230 : // Store the first time the colour appear.
2231 0 : for (int i = 0;i < event.size(); ++i) {
2232 0 : for (int j = 0;j < 3; ++j) {
2233 0 : if (iTempPars[j] == -1 && event.kindJunction(juncs[iJunc]) % 2 == 1 &&
2234 0 : event[i].col() == cols[j]) iTempPars[j] = i;
2235 0 : if (iTempPars[j] == -1 && event.kindJunction(juncs[iJunc]) % 2 == 0 &&
2236 0 : event[i].acol() == cols[j]) iTempPars[j] = i;
2237 : }
2238 : }
2239 :
2240 0 : for (int i = 0;i < 3;++i) {
2241 0 : if (iTempPars[i] >= 0) iPar.push_back(iTempPars[i]);
2242 0 : else addJunctionIndices(event, cols[i], iPar, usedJuncs);
2243 : }
2244 0 : }
2245 :
2246 : // Done.
2247 0 : }
2248 :
2249 : // ------------------------------------------------------------------
2250 :
2251 : // Find all particles connected to a junction system.
2252 :
2253 : void ColourReconnection::addJunctionIndices(int iSinglePar, vector<int> &iPar,
2254 : vector<int> &usedJuncs) {
2255 :
2256 : // Check if junction was already considered.
2257 0 : int iJun = -(1 + iSinglePar/10);
2258 0 : for (int i = 0;i < int(usedJuncs.size()); ++i)
2259 0 : if (iJun == usedJuncs[i]) return;
2260 :
2261 : // Add particles connected to the junction.
2262 0 : usedJuncs.push_back(iJun);
2263 0 : for (int i = 0; i < 3; ++i)
2264 0 : if (junctions[iJun].kind() % 2 == 1) {
2265 0 : int iCol = junctions[iJun].dips[i]->iCol;
2266 0 : if (iCol >= 0) iPar.push_back(iCol);
2267 0 : else addJunctionIndices(iCol, iPar, usedJuncs);
2268 0 : } else {
2269 0 : int iAcol = junctions[iJun].dips[i]->iAcol;
2270 0 : if (iAcol >= 0) iPar.push_back(iAcol);
2271 0 : else addJunctionIndices(iAcol, iPar, usedJuncs);
2272 0 : }
2273 0 : }
2274 :
2275 : // ------------------------------------------------------------------
2276 :
2277 : // Calculate the invariant mass of a dipole.
2278 :
2279 : double ColourReconnection::mDip(ColourDipole* dip) {
2280 :
2281 : // If double junction no invariant mass is given.
2282 0 : if (dip->isJun && dip->isAntiJun) return 1e9;
2283 : // If it has a single junction end.
2284 0 : else if (dip->isJun || dip->isAntiJun) {
2285 0 : int iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2;
2286 0 : getJunctionIndices(dip, iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2);
2287 0 : if (i0 == i1)
2288 0 : return particles[i0].m();
2289 0 : if (i1 < 0)
2290 0 : return 1e9;
2291 0 : return m(particles[i0].p(),particles[i1].p());
2292 0 : } // No junction ends.
2293 : else {
2294 0 : if (dip->iCol == dip->iAcol)
2295 0 : return particles[dip->iCol].m();
2296 : else
2297 0 : return m(particles[dip->iCol].p(),particles[dip->iAcol].p());
2298 : }
2299 :
2300 0 : }
2301 :
2302 : // ------------------------------------------------------------------
2303 :
2304 : // Print dipoles, intended for debuggning purposes.
2305 :
2306 : void ColourReconnection::listDipoles(bool onlyActive, bool onlyReal) {
2307 :
2308 0 : cout << " --- listing dipoles ---" << endl;
2309 0 : for (int i = 0; i < int(dipoles.size()); ++i) {
2310 0 : if (onlyActive && !dipoles[i]->isActive)
2311 : continue;
2312 0 : if (onlyReal && !dipoles[i]->isReal)
2313 : continue;
2314 0 : dipoles[i]->print();
2315 0 : }
2316 0 : cout << " --- finished listing ---" << endl;
2317 :
2318 0 : }
2319 :
2320 : // ------------------------------------------------------------------
2321 :
2322 : // Print particles, intended for debugging purposes.
2323 :
2324 : void ColourReconnection::listParticles() {
2325 :
2326 0 : for (int i = 0; i < int(particles.size()); ++i) {
2327 0 : const ColourParticle& pt = particles[i];
2328 :
2329 : // Basic line for a particle, always printed.
2330 0 : cout << setw(6) << i << setw(10) << pt.id() << " " << left
2331 0 : << setw(18) << pt.nameWithStatus(18) << right << setw(4)
2332 0 : << pt.status() << setw(6) << pt.mother1() << setw(6)
2333 0 : << pt.mother2() << setw(6) << pt.daughter1() << setw(6)
2334 0 : << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol()
2335 0 : << setprecision(3)
2336 0 : << setw(11) << pt.px() << setw(11) << pt.py() << setw(11)
2337 0 : << pt.pz() << setw(11) << pt.e() << setw(11) << pt.m();
2338 0 : for (int j = 0;j < int(pt.activeDips.size());++j)
2339 0 : cout << setw(10) << pt.activeDips[j];
2340 0 : cout << "\n";
2341 : }
2342 :
2343 0 : }
2344 :
2345 : // ------------------------------------------------------------------
2346 :
2347 : // Print junctions, intended for debugging purposes.
2348 :
2349 : void ColourReconnection::listJunctions() {
2350 :
2351 0 : cout << " --- listing junctions ---" << endl;
2352 0 : for (int i = 0; i < int(junctions.size()); ++i)
2353 0 : junctions[i].print();
2354 0 : cout << " --- finished listing ---" << endl;
2355 :
2356 0 : }
2357 :
2358 : // ------------------------------------------------------------------
2359 :
2360 : // Allow colour reconnections by mergings of MPI collision subsystems.
2361 : // iRec is system that may be reconnected, by moving its gluons to iSys,
2362 : // where minimal pT (or equivalently Lambda) is used to pick location.
2363 : // Therefore all dipoles in iSys have to be found, and all gluons in iRec.
2364 : // Matching q-qbar pairs are treated by analogy with gluons.
2365 : // Note: owing to rescatterings some outgoing partons must be skipped.
2366 :
2367 : bool ColourReconnection::reconnectMPIs( Event& event, int oldSize) {
2368 :
2369 : // References to beams to simplify indexing.
2370 0 : BeamParticle& beamA = *beamAPtr;
2371 0 : BeamParticle& beamB = *beamBPtr;
2372 :
2373 : // Prepare record of which systems should be merged onto another.
2374 : // The iSys system must have colour in final state to attach to it.
2375 0 : nSys = partonSystemsPtr->sizeSys();
2376 0 : vector<int> iMerge(nSys);
2377 0 : vector<bool> hasColour(nSys);
2378 0 : for (int iSys = 0; iSys < nSys; ++iSys) {
2379 0 : iMerge[iSys] = iSys;
2380 : bool hasCol = false;
2381 0 : for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
2382 0 : int iNow = partonSystemsPtr->getOut( iSys, iMem);
2383 0 : if (event[iNow].isFinal() && (event[iNow].col() > 0
2384 0 : || event[iNow].acol() > 0) ) {
2385 : hasCol = true;
2386 0 : break;
2387 : }
2388 0 : }
2389 0 : hasColour[iSys] = hasCol;
2390 : }
2391 :
2392 : // Loop over systems to decide which should be reconnected.
2393 0 : for (int iRec = nSys - 1; iRec > 0; --iRec) {
2394 :
2395 : // Determine reconnection strength from pT scale of system.
2396 0 : double pT2Rec = pow2( partonSystemsPtr->getPTHat(iRec) );
2397 0 : double probRec = pT20Rec / (pT20Rec + pT2Rec);
2398 :
2399 : // Loop over other systems iSys at higher pT scale and
2400 : // decide whether to reconnect the iRec gluons onto one of them.
2401 0 : for (int iSys = iRec - 1; iSys >= 0; --iSys)
2402 0 : if (hasColour[iSys] && probRec > rndmPtr->flat()) {
2403 :
2404 : // The iRec system and all merged with it to be merged with iSys.
2405 0 : iMerge[iRec] = iSys;
2406 0 : for (int iRec2 = iRec + 1; iRec2 < nSys; ++iRec2)
2407 0 : if (iMerge[iRec2] == iRec) iMerge[iRec2] = iSys;
2408 :
2409 : // Once a system has been merged do not test it anymore.
2410 0 : break;
2411 : }
2412 : }
2413 :
2414 : // Loop over systems. Check whether other systems to be merged with it.
2415 0 : for (int iSys = 0; iSys < nSys; ++iSys) {
2416 : int nMerge = 0;
2417 0 : for (int iRec = iSys + 1; iRec < nSys; ++iRec)
2418 0 : if (iMerge[iRec] == iSys) ++nMerge;
2419 0 : if (nMerge == 0) continue;
2420 :
2421 : // Incoming partons not counted if rescattered.
2422 0 : int iInASys = partonSystemsPtr->getInA(iSys);
2423 0 : bool hasInA = (beamA[iSys].isFromBeam());
2424 0 : int iInBSys = partonSystemsPtr->getInB(iSys);
2425 0 : bool hasInB = (beamB[iSys].isFromBeam());
2426 :
2427 : // Begin find dipoles in iSys system.
2428 0 : vector<BeamDipole> bmdipoles;
2429 0 : int sizeOut = partonSystemsPtr->sizeOut(iSys);
2430 0 : for (int iMem = 0; iMem < sizeOut; ++iMem) {
2431 :
2432 : // Find colour dipoles to beam remnant.
2433 0 : int iNow = partonSystemsPtr->getOut( iSys, iMem);
2434 0 : if (!event[iNow].isFinal()) continue;
2435 0 : int col = event[iNow].col();
2436 0 : if (col > 0) {
2437 0 : if (hasInA && event[iInASys].col() == col)
2438 0 : bmdipoles.push_back( BeamDipole( col, iNow, iInASys ) );
2439 0 : else if (hasInB && event[iInBSys].col() == col)
2440 0 : bmdipoles.push_back( BeamDipole( col, iNow, iInBSys ) );
2441 :
2442 : // Find colour dipole between final-state partons.
2443 0 : else for (int iMem2 = 0; iMem2 < sizeOut; ++iMem2)
2444 0 : if (iMem2 != iMem) {
2445 0 : int iNow2 = partonSystemsPtr->getOut( iSys, iMem2);
2446 0 : if (!event[iNow2].isFinal()) continue;
2447 0 : if (event[iNow2].acol() == col) {
2448 0 : bmdipoles.push_back( BeamDipole( col, iNow, iNow2) );
2449 0 : break;
2450 : }
2451 0 : }
2452 : }
2453 :
2454 : // Find anticolour dipoles to beam remnant.
2455 0 : int acol = event[iNow].acol();
2456 0 : if (acol > 0) {
2457 0 : if (hasInA && event[iInASys].acol() == acol)
2458 0 : bmdipoles.push_back( BeamDipole( acol, iInASys, iNow ) );
2459 0 : else if (hasInB && event[iInBSys].acol() == acol)
2460 0 : bmdipoles.push_back( BeamDipole( acol, iInBSys, iNow ) );
2461 : }
2462 0 : }
2463 :
2464 : // Skip mergings if no dipoles found.
2465 0 : if (bmdipoles.size() == 0) continue;
2466 :
2467 : // Find dipole sizes.
2468 0 : for (int iDip = 0; iDip < int(bmdipoles.size()); ++iDip)
2469 0 : bmdipoles[iDip].p1p2 = event[bmdipoles[iDip].iCol].p()
2470 0 : * event[bmdipoles[iDip].iAcol].p();
2471 :
2472 : // Loop over systems iRec to be merged with iSys.
2473 0 : for (int iRec = iSys + 1; iRec < nSys; ++iRec) {
2474 0 : if (iMerge[iRec] != iSys) continue;
2475 :
2476 : // Information on iRec. Vectors for gluons and anything else.
2477 0 : int sizeRec = partonSystemsPtr->sizeOut(iRec);
2478 0 : int iInARec = partonSystemsPtr->getInA(iRec);
2479 0 : int iInBRec = partonSystemsPtr->getInB(iRec);
2480 : int nGluRec = 0;
2481 0 : vector<int> iGluRec;
2482 0 : vector<double> pT2GluRec;
2483 : int nAnyRec = 0;
2484 0 : vector<int> iAnyRec;
2485 0 : vector<bool> freeAnyRec;
2486 :
2487 : // Copy of gluon positions in descending order.
2488 0 : for (int iMem = 0; iMem < sizeRec; ++iMem) {
2489 0 : int iNow = partonSystemsPtr->getOut( iRec, iMem);
2490 0 : if (!event[iNow].isFinal()) continue;
2491 0 : if (event[iNow].isGluon()) {
2492 0 : ++nGluRec;
2493 0 : iGluRec.push_back( iNow );
2494 0 : pT2GluRec.push_back( event[iNow].pT2() );
2495 0 : for (int i = nGluRec - 1; i > 1; --i) {
2496 0 : if (pT2GluRec[i - 1] > pT2GluRec[i]) break;
2497 0 : swap( iGluRec[i - 1], iGluRec[i] );
2498 0 : swap( pT2GluRec[i - 1], pT2GluRec[i] );
2499 : }
2500 : // Copy of anything else, mainly quarks, in no particular order.
2501 0 : } else {
2502 0 : ++nAnyRec;
2503 0 : iAnyRec.push_back( iNow );
2504 0 : freeAnyRec.push_back( true );
2505 : }
2506 0 : }
2507 :
2508 : // For each gluon in iRec now find the dipole that gives the smallest
2509 : // (pGlu * pI) (pGlu * pJ) / (pI * pJ), i.e. minimal pT (and Lambda).
2510 0 : for (int iGRec = 0; iGRec < nGluRec; ++iGRec) {
2511 0 : int iGlu = iGluRec[iGRec];
2512 0 : Vec4 pGlu = event[iGlu].p();
2513 : int iDipMin = 0;
2514 0 : double pT2DipMin = sCM;
2515 0 : for (int iDip = 0; iDip < int(bmdipoles.size()); ++iDip) {
2516 0 : double pT2Dip = (pGlu * event[bmdipoles[iDip].iCol].p())
2517 0 : * (pGlu * event[bmdipoles[iDip].iAcol].p()) / bmdipoles[iDip].p1p2;
2518 0 : if (pT2Dip < pT2DipMin) {
2519 : iDipMin = iDip;
2520 : pT2DipMin = pT2Dip;
2521 0 : }
2522 : }
2523 :
2524 : // Attach the gluon to the dipole, i.e. split the dipole in two.
2525 0 : int colGlu = event[iGlu].col();
2526 0 : int acolGlu = event[iGlu].acol();
2527 0 : int colDip = bmdipoles[iDipMin].col;
2528 0 : int iColDip = bmdipoles[iDipMin].iCol;
2529 0 : int iAcolDip = bmdipoles[iDipMin].iAcol;
2530 0 : event[iGlu].acol( colDip );
2531 0 : if (event[iAcolDip].acol() == colDip)
2532 0 : event[iAcolDip].acol( colGlu );
2533 0 : else event[iAcolDip].col( colGlu );
2534 0 : bmdipoles[iDipMin].iAcol = iGlu;
2535 0 : bmdipoles[iDipMin].p1p2 = event[iColDip].p() * pGlu;
2536 0 : bmdipoles.push_back( BeamDipole( colGlu, iGlu, iAcolDip ) );
2537 0 : bmdipoles.back().p1p2 = pGlu * event[iAcolDip].p();
2538 :
2539 : // Remove gluon from old system: reconnect colours.
2540 0 : for (int i = oldSize; i < event.size(); ++i)
2541 0 : if (i != iGlu && i != iAcolDip) {
2542 0 : if (event[i].isFinal()) {
2543 0 : if (event[i].acol() == colGlu) event[i].acol( acolGlu );
2544 : } else {
2545 0 : if (event[i].col() == colGlu) event[i].col( acolGlu );
2546 : }
2547 : }
2548 :
2549 : // Update any junction legs that match reconnected dipole.
2550 0 : for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
2551 :
2552 : // Only junctions need to be updated, not antijunctions.
2553 0 : if (event.kindJunction(iJun) % 2 == 0) continue;
2554 0 : for (int leg = 0; leg < 3; ++leg) {
2555 0 : int col = event.colJunction(iJun, leg);
2556 0 : if (col == colDip)
2557 0 : event.colJunction(iJun, leg, colGlu);
2558 : }
2559 0 : }
2560 :
2561 0 : }
2562 :
2563 : // See if any matching quark-antiquark pairs among the rest.
2564 0 : for (int iQRec = 0; iQRec < nAnyRec; ++iQRec) {
2565 0 : int iQ = iAnyRec[iQRec];
2566 0 : int idQ = event[iQ].id();
2567 0 : if (freeAnyRec[iQRec] && idQ > 0 && idQ < 6)
2568 0 : for (int iQbarRec = 0; iQbarRec < nAnyRec; ++iQbarRec) {
2569 0 : int iQbar = iAnyRec[iQbarRec];
2570 0 : if (freeAnyRec[iQbarRec] && event[iQbar].id() == -idQ) {
2571 :
2572 : // Check that these can be traced back to same gluon splitting.
2573 : // For now also avoid qqbar pairs produced in rescatterings.??
2574 0 : int iTopQ = event[iQ].iTopCopyId();
2575 0 : int iTopQbar = event[iQbar].iTopCopyId();
2576 0 : int iMother = event[iTopQ].mother1();
2577 0 : if (event[iTopQbar].mother1() == iMother
2578 0 : && event[iMother].isGluon() && event[iMother].status() != -34
2579 0 : && event[iMother + 1].status() != -34 ) {
2580 :
2581 : // Now find the dipole that gives the smallest
2582 : // ((pQ + pQbar) * pI) ((pQ + pQbar) * pJ) / (pI * pJ).
2583 0 : Vec4 pGlu = event[iQ].p() + event[iQbar].p();
2584 : int iDipMin = 0;
2585 0 : double pT2DipMin = sCM;
2586 0 : for (int iDip = 0; iDip < int(bmdipoles.size()); ++iDip) {
2587 0 : double pT2Dip = (pGlu * event[bmdipoles[iDip].iCol].p())
2588 0 : * (pGlu * event[bmdipoles[iDip].iAcol].p())
2589 0 : / bmdipoles[iDip].p1p2;
2590 0 : if (pT2Dip < pT2DipMin) {
2591 : iDipMin = iDip;
2592 : pT2DipMin = pT2Dip;
2593 0 : }
2594 : }
2595 :
2596 : // Attach the q-qbar pair to the dipole, i.e. split the dipole.
2597 0 : int colGlu = event[iQ].col();
2598 0 : int acolGlu = event[iQbar].acol();
2599 0 : int colDip = bmdipoles[iDipMin].col;
2600 0 : int iColDip = bmdipoles[iDipMin].iCol;
2601 0 : int iAcolDip = bmdipoles[iDipMin].iAcol;
2602 0 : event[iQbar].acol( colDip );
2603 0 : if (event[iAcolDip].acol() == colDip)
2604 0 : event[iAcolDip].acol( colGlu );
2605 0 : else event[iAcolDip].col( colGlu );
2606 0 : bmdipoles[iDipMin].iAcol = iQbar;
2607 0 : bmdipoles[iDipMin].p1p2 = event[iColDip].p() * event[iQbar].p();
2608 0 : bmdipoles.push_back( BeamDipole( colGlu, iQ, iAcolDip ) );
2609 0 : bmdipoles.back().p1p2 = event[iQ].p() * event[iAcolDip].p();
2610 :
2611 : // Remove q-qbar pair from old system: reconnect colours.
2612 0 : freeAnyRec[iQRec] = false;
2613 0 : freeAnyRec[iQbarRec] = false;
2614 0 : for (int i = oldSize; i < event.size(); ++i)
2615 0 : if (i != iQRec && i != iQbarRec && i != iColDip
2616 0 : && i != iAcolDip) {
2617 0 : if (event[i].isFinal()) {
2618 0 : if (event[i].acol() == colGlu) event[i].acol( acolGlu );
2619 : } else {
2620 0 : if (event[i].col() == colGlu) event[i].col( acolGlu );
2621 : }
2622 : }
2623 :
2624 : // Update any junction legs that match reconnected dipole.
2625 0 : for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
2626 :
2627 : // Only junctions need to be updated, not antijunctions.
2628 0 : if (event.kindJunction(iJun) % 2 == 0) continue;
2629 0 : for (int leg = 0; leg < 3; ++leg) {
2630 0 : int col = event.colJunction(iJun, leg);
2631 0 : if (col == colDip)
2632 0 : event.colJunction(iJun, leg, colGlu);
2633 : }
2634 0 : }
2635 :
2636 : // Done with processing of q-qbar pairs.
2637 0 : }
2638 0 : }
2639 0 : }
2640 : }
2641 :
2642 : // If only two beam gluons left of system, set their colour = anticolour.
2643 : // Used by BeamParticle::remnantColours to skip irrelevant gluons.
2644 0 : if ( event[iInARec].isGluon() && !event[iInARec].isRescatteredIncoming()
2645 0 : && event[iInBRec].isGluon() && !event[iInBRec].isRescatteredIncoming()
2646 0 : && event[iInARec].col() == event[iInBRec].acol()
2647 0 : && event[iInARec].acol() == event[iInBRec].col() ) {
2648 0 : event[iInARec].acol( event[iInARec].col() );
2649 0 : event[iInBRec].acol( event[iInBRec].col() );
2650 0 : }
2651 :
2652 : // End of loops over iRec and iSys systems.
2653 0 : }
2654 0 : }
2655 :
2656 : // Done.
2657 : return true;
2658 0 : }
2659 :
2660 : //--------------------------------------------------------------------------
2661 :
2662 : // Find the neighbour to the anticolour side. Return false if the dipole
2663 : // is connected to a junction or the new particle has a junction inside of it.
2664 :
2665 : bool ColourReconnection::findAntiNeighbour(ColourDipole*& dip) {
2666 : // If only one active dipole, it has to be an antiquark.
2667 0 : if (int(particles[dip->iAcol].activeDips.size()) == 1)
2668 0 : return false;
2669 :
2670 : // Has to have to active dipoles, otherwise something went wrong.
2671 0 : if (int(particles[dip->iAcol].activeDips.size()) != 2) {
2672 0 : infoPtr->errorMsg("Warning in ColourReconnection::findAntiNeighbour: "
2673 : "Wrong number of active dipoles");
2674 0 : return false;
2675 : }
2676 :
2677 : // Otherwise find new dipole.
2678 0 : if (dip == particles[dip->iAcol].activeDips[0])
2679 0 : dip = particles[dip->iAcol].activeDips[1];
2680 0 : else dip = particles[dip->iAcol].activeDips[0];
2681 :
2682 : // Do not allow the new dipole to be connected to an antijunction.
2683 0 : if (dip->isAntiJun || dip->isJun)
2684 0 : return false;
2685 :
2686 : // Do not allow new dipole to have a pseudoparticle with
2687 : // a baryon number inside.
2688 0 : if (int(particles[dip->iAcol].dips.size()) != 1)
2689 0 : return false;
2690 :
2691 0 : return true;
2692 0 : }
2693 :
2694 : //--------------------------------------------------------------------------
2695 :
2696 : // Check that trials do not contain junctions/ unusable pseudoparticles.
2697 :
2698 : bool ColourReconnection::checkJunctionTrials() {
2699 0 : for (int i = 0;i < int(junTrials.size());++i) {
2700 : int minus = 0;
2701 0 : if (junTrials[i].mode == 3)
2702 0 : minus = 1;
2703 0 : for (int j = 0;j < int(junTrials[i].dips.size()) - minus; ++j) {
2704 0 : ColourDipole* dip = junTrials[i].dips[j];
2705 0 : if (dip->isJun || dip->isAntiJun) {
2706 0 : junTrials[i].list();
2707 0 : return false;
2708 : }
2709 0 : if (particles[dip->iCol].dips.size() != 1 ||
2710 0 : particles[dip->iAcol].dips.size() != 1) {
2711 0 : junTrials[i].list();
2712 0 : return false;
2713 : }
2714 0 : }
2715 0 : }
2716 0 : return true;
2717 0 : }
2718 :
2719 :
2720 : //--------------------------------------------------------------------------
2721 :
2722 : // Find the neighbour to the colour side. Return false if the dipole
2723 : // is connected to a junction or the new particle has a junction inside of it.
2724 :
2725 : bool ColourReconnection::findColNeighbour(ColourDipole*& dip) {
2726 : // If only one active dipole, it has to be an antiquark.
2727 0 : if (int(particles[dip->iCol].activeDips.size()) == 1)
2728 0 : return false;
2729 :
2730 : // Has to have to active dipoles, otherwise something went wrong.
2731 0 : if (int(particles[dip->iCol].activeDips.size()) != 2) {
2732 0 : infoPtr->errorMsg("Warning in ColourReconnection::findAntiNeighbour: "
2733 : "Wrong number of active dipoles");
2734 0 : return false;
2735 : }
2736 : // Otherwise find new dipole.
2737 0 : if (dip == particles[dip->iCol].activeDips[0])
2738 0 : dip = particles[dip->iCol].activeDips[1];
2739 0 : else dip = particles[dip->iCol].activeDips[0];
2740 :
2741 : // Do not allow the new dipole to be connected to an antijunction.
2742 0 : if (dip->isJun || dip->isAntiJun)
2743 0 : return false;
2744 :
2745 : // Do not allow new dipole to have a pseudoparticle with
2746 : // a baryon number inside.
2747 0 : if (int(particles[dip->iCol].dips.size()) != 1)
2748 0 : return false;
2749 :
2750 0 : return true;
2751 0 : }
2752 :
2753 : //--------------------------------------------------------------------------
2754 :
2755 : // Store used dipoles for a junction formation.
2756 :
2757 : void ColourReconnection::storeUsedDips(TrialReconnection& trial) {
2758 : // Normal dipole swap.
2759 0 : if (trial.mode == 5) {
2760 :
2761 0 : for (int i = 0;i < 2; ++i) {
2762 0 : ColourDipole* dip = trial.dips[i];
2763 0 : if (dip->iCol < 0)
2764 0 : for (int j = 0;j < 3; ++j)
2765 0 : usedDipoles.push_back(junctions[-(dip->iCol / 10 + 1)].getColDip(j));
2766 0 : if (dip->iAcol < 0)
2767 0 : for (int j = 0;j < 3; ++j)
2768 0 : usedDipoles.push_back(junctions[-(dip->iAcol / 10 + 1)].getColDip(j));
2769 :
2770 0 : usedDipoles.push_back(dip);
2771 0 : }
2772 :
2773 0 : } else {
2774 :
2775 0 : for (int i = 0;i < 4; ++i) {
2776 0 : if (trial.mode == 3 && i == 3)
2777 : continue;
2778 0 : usedDipoles.push_back(trial.dips[i]);
2779 0 : ColourDipole* dip = trial.dips[i];
2780 :
2781 :
2782 0 : while (findAntiNeighbour(dip) && dip != trial.dips[i])
2783 0 : usedDipoles.push_back(dip);
2784 :
2785 0 : dip = trial.dips[i];
2786 0 : while (findColNeighbour(dip) && dip != trial.dips[i])
2787 0 : usedDipoles.push_back(dip);
2788 0 : }
2789 : }
2790 0 : }
2791 :
2792 : //--------------------------------------------------------------------------
2793 :
2794 : // Calculate the difference between the old and new lambda for dipole swap.
2795 :
2796 : double ColourReconnection::getLambdaDiff(ColourDipole* dip1,
2797 : ColourDipole* dip2) {
2798 :
2799 : // Needed to make sure the same dipoles are compared.
2800 0 : vector<ColourDipole*> oldDips, newDips;
2801 :
2802 : // Calculate old string length.
2803 0 : double oldLambda = calculateStringLength(dip1, oldDips)
2804 0 : + calculateStringLength( dip2, oldDips);
2805 :
2806 : // Make test configuration.
2807 0 : swapDipoles(dip1,dip2);
2808 :
2809 : // Calculate new string lengths
2810 0 : double newLambda = calculateStringLength(dip1, newDips)
2811 0 : + calculateStringLength(dip2, newDips);
2812 :
2813 : // Swap back.
2814 0 : swapDipoles(dip1, dip2, true);
2815 :
2816 : // First check if new combination was not useable.
2817 0 : if (newLambda >= 0.5E9 || newLambda < 1.) return -1e9;
2818 :
2819 : // Return the difference.
2820 0 : return oldLambda - newLambda;
2821 0 : }
2822 :
2823 : //--------------------------------------------------------------------------
2824 :
2825 : // Calculate the difference between the old and new lambda.
2826 :
2827 : double ColourReconnection::getLambdaDiff(ColourDipole* dip1,
2828 : ColourDipole* dip2, ColourDipole* dip3, ColourDipole* dip4, int mode) {
2829 :
2830 : // Calculate old lambda measure.
2831 :
2832 0 : double oldLambda = calculateStringLength(dip1->iCol, dip1->iAcol)
2833 0 : + calculateStringLength(dip2->iCol, dip2->iAcol);
2834 0 : if (dip3 != dip1)
2835 0 : oldLambda += calculateStringLength(dip3->iCol, dip3->iAcol);
2836 0 : if (dip4 != 0 && dip4 != dip2)
2837 0 : oldLambda += calculateStringLength(dip4->iCol, dip4->iAcol);
2838 :
2839 : // Calculate new lambda.
2840 : double newLambda = 0;
2841 :
2842 0 : if (mode == 0)
2843 0 : newLambda = calculateDoubleJunctionLength(dip1->iCol, dip2->iCol,
2844 0 : dip1->iAcol, dip2->iAcol);
2845 0 : else if (mode == 1) {
2846 0 : if (dip2 == dip4)
2847 0 : newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip3->iCol)
2848 0 : + calculateJunctionLength(dip1->iAcol, dip2->iAcol, dip3->iAcol);
2849 : else
2850 : newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip3->iCol)
2851 0 : + calculateJunctionLength(dip2->iAcol, dip3->iAcol, dip4->iAcol)
2852 0 : + calculateStringLength(dip4->iCol, dip1->iAcol);
2853 : }
2854 :
2855 0 : else if (mode == 2) {
2856 0 : if (dip1 == dip3)
2857 0 : newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip4->iCol)
2858 0 : + calculateJunctionLength(dip1->iAcol, dip2->iAcol, dip4->iAcol);
2859 : else
2860 : newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip4->iCol)
2861 0 : + calculateJunctionLength(dip1->iAcol, dip3->iAcol, dip4->iAcol)
2862 0 : + calculateStringLength(dip3->iCol, dip2->iAcol);
2863 : }
2864 :
2865 : // Triple junction connection.
2866 0 : else if (mode == 3)
2867 0 : newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip3->iCol)
2868 0 : + calculateJunctionLength(dip1->iAcol, dip2->iAcol, dip3->iAcol);
2869 :
2870 : // First check if new combination was not useable.
2871 0 : if (newLambda >= 0.5E9 || newLambda < 1.) return -1e9;
2872 :
2873 : // Returning result.
2874 0 : return oldLambda - newLambda;
2875 :
2876 0 : }
2877 :
2878 : //--------------------------------------------------------------------------
2879 :
2880 : // Change colour structure to describe the reconnection in juncTrial.
2881 :
2882 : void ColourReconnection::doDipoleTrial(TrialReconnection& trial) {
2883 :
2884 : // Store for easier use.
2885 0 : ColourDipole* dip1 = trial.dips[0];
2886 0 : ColourDipole* dip2 = trial.dips[1];
2887 :
2888 : // If both acols ends are normal particles.
2889 0 : if (dip1->iAcol >= 0 && dip2->iAcol >= 0) {
2890 0 : swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front()->iAcol,
2891 0 : particles[dip2->iAcol].dips[dip2->iAcolLeg].front()->iAcol);
2892 0 : swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front(),
2893 0 : particles[dip2->iAcol].dips[dip2->iAcolLeg].front());
2894 : // If only dip1 has normal acol end.
2895 0 : } else if (dip1->iAcol >= 0) {
2896 0 : swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front()->iAcol,
2897 0 : junctions[-(dip2->iAcol / 10 + 1)].dipsOrig[-dip2->iAcol % 10]->iAcol);
2898 0 : swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front(),
2899 0 : junctions[-(dip2->iAcol / 10 + 1)].dipsOrig[-dip2->iAcol % 10]);
2900 : // If only dip2 has normal acol end.
2901 0 : } else if (dip2->iAcol >= 0) {
2902 0 : swap(particles[dip2->iAcol].dips[dip2->iAcolLeg].front()->iAcol,
2903 0 : junctions[-(dip1->iAcol / 10 + 1)].dipsOrig[-dip1->iAcol % 10]->iAcol);
2904 0 : swap(particles[dip2->iAcol].dips[dip2->iAcolLeg].front(),
2905 0 : junctions[-(dip1->iAcol / 10 + 1)].dipsOrig[-dip1->iAcol % 10]);
2906 : // If both ends are junctions.
2907 0 : } else {
2908 0 : swap(junctions[ -(dip1->iAcol / 10 + 1) ].dipsOrig[
2909 0 : -dip1->iAcol % 10 ]->iAcol,
2910 0 : junctions[ -(dip2->iAcol / 10 + 1) ].dipsOrig[
2911 0 : -dip2->iAcol % 10 ]->iAcol);
2912 0 : swap(junctions[ -(dip1->iAcol / 10 + 1) ].dipsOrig[ -dip1->iAcol % 10],
2913 0 : junctions[ -(dip2->iAcol / 10 + 1) ].dipsOrig[ -dip2->iAcol % 10] );
2914 : }
2915 :
2916 : // Swap the dipoles.
2917 0 : swapDipoles(dip1, dip2);
2918 :
2919 : // If new particles are below treshhold, form pseudoParticles.
2920 0 : if (mDip(dip1) < m0) makePseudoParticle(dip1, 110, true);
2921 0 : if (mDip(dip2) < m0) makePseudoParticle(dip2, 110, true);
2922 :
2923 : // Done.
2924 :
2925 0 : }
2926 :
2927 : //--------------------------------------------------------------------------
2928 :
2929 : // Update the list of dipole trial swaps to account for latest swap.
2930 :
2931 : void ColourReconnection::updateDipoleTrials() {
2932 :
2933 : // Remove any dipTrials that contains an used dipole.
2934 0 : for (int i = 0; i < int(dipTrials.size()); ++i)
2935 0 : for (int j = 0;j < 2; ++j) {
2936 0 : if (binary_search(usedDipoles.begin(), usedDipoles.end(),
2937 0 : dipTrials[i].dips[j]) ) {
2938 0 : dipTrials.erase(dipTrials.begin() + i);
2939 0 : i--;
2940 0 : break;
2941 : }
2942 : }
2943 :
2944 : // Make list of active dipoles.
2945 0 : vector<ColourDipole*> activeDipoles;
2946 0 : for (int i = 0;i < int(dipoles.size()); ++i)
2947 0 : if (dipoles[i]->isActive)
2948 0 : activeDipoles.push_back(dipoles[i]);
2949 :
2950 : // Loop over list of used dipoles and create new trial reconnections.
2951 0 : for (int i = 0;i < int(usedDipoles.size()); ++i)
2952 0 : if (usedDipoles[i]->isActive)
2953 0 : for (int j = 0; j < int(activeDipoles.size()); ++j)
2954 0 : singleReconnection(usedDipoles[i], activeDipoles[j]);
2955 :
2956 0 : }
2957 :
2958 : //--------------------------------------------------------------------------
2959 :
2960 : // Update the list of dipole trial swaps to account for latest swap.
2961 :
2962 : void ColourReconnection::updateJunctionTrials() {
2963 :
2964 : // Remove any junTrials that contains an used dipole.
2965 0 : for (int i = 0; i < int(junTrials.size()); ++i)
2966 0 : for (int j = 0; j < 4; ++j) {
2967 0 : if (binary_search(usedDipoles.begin(), usedDipoles.end(),
2968 0 : junTrials[i].dips[j]) ) {
2969 0 : junTrials.erase(junTrials.begin() + i);
2970 0 : i--;
2971 0 : break;
2972 : }
2973 : }
2974 :
2975 : // Make list of active dipoles.
2976 0 : vector<ColourDipole*> activeDipoles;
2977 0 : for (int i = 0;i < int(dipoles.size()); ++i)
2978 0 : if (dipoles[i]->isActive)
2979 0 : activeDipoles.push_back(dipoles[i]);
2980 :
2981 : // Loop over used dipoles and form new junction trials.
2982 0 : for (int i = 0;i < int(usedDipoles.size()); ++i)
2983 0 : if (usedDipoles[i]->isActive)
2984 0 : for (int j = 0; j < int(activeDipoles.size()); ++j)
2985 0 : singleJunction(usedDipoles[i], activeDipoles[j]);
2986 :
2987 : // Loop over used dipoles and form new junction trials.
2988 0 : for (int i = 0;i < int(usedDipoles.size()); ++i)
2989 0 : if (usedDipoles[i]->isActive)
2990 0 : for (int j = 0; j < int(activeDipoles.size()); ++j)
2991 0 : for (int k = j + 1; k < int(activeDipoles.size()); ++k)
2992 0 : singleJunction(usedDipoles[i], activeDipoles[j], activeDipoles[k]);
2993 :
2994 0 : }
2995 :
2996 : //--------------------------------------------------------------------------
2997 :
2998 : // Change colour structure to describe the reconnection in juncTrial.
2999 :
3000 : void ColourReconnection::doJunctionTrial(Event& event,
3001 : TrialReconnection& juncTrial) {
3002 :
3003 0 : int mode = juncTrial.mode;
3004 : // If trial mode is 3 (three dipoles -> 2 junctions) use its own update.
3005 0 : if (mode == 3) {
3006 0 : doTripleJunctionTrial(event, juncTrial);
3007 0 : return;
3008 : }
3009 :
3010 : // Store dipoles and numbers for easier acces.
3011 0 : ColourDipole* dip1 = juncTrial.dips[0];
3012 0 : ColourDipole* dip2 = juncTrial.dips[1];
3013 0 : ColourDipole* dip3 = juncTrial.dips[2];
3014 0 : ColourDipole* dip4 = juncTrial.dips[3];
3015 :
3016 0 : int iCol1 = dip1->iCol;
3017 0 : int iCol2 = dip2->iCol;
3018 0 : int iCol3 = dip3->iCol;
3019 0 : int iCol4 = dip4->iCol;
3020 0 : int iAcol1 = dip1->iAcol;
3021 0 : int iAcol2 = dip2->iAcol;
3022 0 : int iAcol3 = dip3->iAcol;
3023 0 : int iAcol4 = dip4->iAcol;
3024 :
3025 0 : int oldCol1 = dip1->col;
3026 0 : int oldCol2 = dip2->col;
3027 0 : int oldCol3 = dip3->col;
3028 0 : int oldCol4 = dip4->col;
3029 :
3030 : // New colour tags needed, since three more dipoles will be made.
3031 0 : int newCol1 = event.nextColTag();
3032 0 : int newCol2 = event.nextColTag();
3033 0 : int newCol3 = event.nextColTag();
3034 :
3035 : // Add new formation times.
3036 0 : double mCalc = (particles[iCol1].p() + particles[iAcol1].p() +
3037 0 : particles[iCol2].p() + particles[iAcol2].p() +
3038 0 : particles[iCol3].p() + particles[iAcol3].p() +
3039 0 : particles[iCol4].p() + particles[iAcol4].p()).mCalc();
3040 0 : formationTimes[newCol1] = mCalc;
3041 0 : formationTimes[newCol2] = mCalc;
3042 0 : formationTimes[newCol3] = mCalc;
3043 :
3044 : // Need to make 3 new real dipoles and 3 active dipoles.
3045 :
3046 : // First make dipoles between junctions.
3047 :
3048 : // Find the junction colour.
3049 0 : int junCol = 3 * (3 - (dip1->colReconnection / 3)
3050 0 : - (dip2->colReconnection / 3) ) + dip1->colReconnection % 3;
3051 :
3052 : // if other than 9 colours.
3053 0 : if (nReconCols != 9) {
3054 0 : while (junCol < 0 || junCol % 3 != dip1->colReconnection % 3 ||
3055 0 : junCol == dip1->colReconnection || junCol == dip2->colReconnection)
3056 0 : junCol = int(nReconCols * rndmPtr->flat());
3057 : }
3058 : // Need one active and one real dipole.
3059 0 : int iJun = junctions.size();
3060 0 : int iAntiJun = junctions.size() + 1;
3061 :
3062 : // Store real dipoles.
3063 0 : ColourDipole* dip3real = particles[iCol3].dips[dip3->iColLeg].back();
3064 0 : ColourDipole* dip4real = particles[iCol4].dips[dip4->iColLeg].back();
3065 :
3066 : // If the junction and antijunction are directly connected.
3067 : int iActive1 = 0, iReal1 = 0;
3068 0 : if (mode == 0) {
3069 0 : dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10 + 2) ,
3070 0 : -( iJun * 10 + 10 + 2), junCol, true, true, false, true));
3071 0 : iReal1 = dipoles.size() - 1;
3072 0 : dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10 + 2) ,
3073 : -( iJun * 10 + 10 + 2), junCol, true, true));
3074 0 : iActive1 = dipoles.size() - 1;
3075 0 : } else if (mode == 1) {
3076 0 : int iCol3real = particles[iCol3].dips[dip3->iColLeg].back()->iCol;
3077 0 : dipoles.push_back(new ColourDipole(newCol1, iCol3real ,
3078 0 : -( iJun * 10 + 10 + 2), junCol, true, false, false, true));
3079 0 : iReal1 = dipoles.size() - 1;
3080 0 : particles[iCol3].dips[dip3->iColLeg].back() = dipoles.back();
3081 0 : dipoles.push_back(new ColourDipole(newCol1, dip3->iCol,
3082 : -( iJun * 10 + 10 + 2), junCol, true, false));
3083 0 : iActive1 = dipoles.size() - 1;
3084 0 : } else if (mode == 2) {
3085 0 : int iCol4real = particles[iCol4].dips[dip4->iColLeg].back()->iCol;
3086 0 : dipoles.push_back(new ColourDipole(newCol1, iCol4real,
3087 0 : -( iJun * 10 + 10 + 2), junCol, true, false, false, true));
3088 0 : iReal1 = dipoles.size() - 1;
3089 0 : particles[iCol4].dips[dip4->iColLeg].back() = dipoles.back();
3090 0 : dipoles.push_back(new ColourDipole(newCol1, dip4->iCol,
3091 : -( iJun * 10 + 10 + 2), junCol, true, false));
3092 0 : iActive1 = dipoles.size() - 1;
3093 0 : }
3094 :
3095 : // Now make dipole between anti junction and iAcol1.
3096 : // Start by finding real iAcol.
3097 0 : int iAcol3real = particles[iAcol3].dips[dip3->iAcolLeg].front()->iAcol;
3098 0 : dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10),
3099 0 : iAcol3real, dip3->colReconnection, false, true, false, true));
3100 0 : int iReal2 = dipoles.size() - 1;
3101 0 : particles[iAcol3].dips[dip3->iAcolLeg].front() = dipoles.back();
3102 :
3103 0 : dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10),
3104 0 : iAcol3, dip3->colReconnection, false, true));
3105 0 : dipoles.back()->iAcolLeg = dip3->iAcolLeg;
3106 0 : int iActive2 = dipoles.size() - 1;
3107 :
3108 : // Now make dipole between anti junction and iAcol1.
3109 : // Start by finding real iAcol.
3110 0 : int iAcol4real = particles[iAcol4].dips[dip4->iAcolLeg].front()->iAcol;
3111 0 : dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 1),
3112 0 : iAcol4real, dip4->colReconnection, false, true, false, true));
3113 0 : int iReal3 = dipoles.size() - 1;
3114 0 : particles[iAcol4].dips[dip4->iAcolLeg].front() = dipoles.back();
3115 :
3116 0 : dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 1),
3117 0 : iAcol4, dip4->colReconnection, false, true));
3118 0 : dipoles.back()->iAcolLeg = dip4->iAcolLeg;
3119 0 : int iActive3 = dipoles.size() - 1;
3120 :
3121 : // Update already existing dipoles, start by internal dipoles.
3122 : // Now take dipoles connected to the anti junction
3123 : // and a possible gluon-gluon connection.
3124 0 : if (mode == 1) {
3125 0 : if (dip2 == dip4) {
3126 :
3127 : // Update real dipole.
3128 0 : dip3real->iAcol = particles[dip1->iAcol].dips[dip1->iAcolLeg].
3129 0 : front()->iAcol;
3130 0 : dip3real->iCol = -( iAntiJun * 10 + 10 + 2);
3131 0 : dip3real->isAntiJun = true;
3132 :
3133 : // Update active dipoles.
3134 0 : dip3->iAcol = dip1->iAcol;
3135 0 : dip3->iAcolLeg = dip1->iAcolLeg;
3136 0 : dip3->isAntiJun = true;
3137 0 : dip3->iCol = -( iAntiJun * 10 + 10 + 2);
3138 0 : dip3->iColLeg = 0;
3139 :
3140 : // Store real dipole
3141 0 : particles[dip3->iAcol].dips[dip3->iAcolLeg].front() = dip3real;
3142 :
3143 0 : } else {
3144 :
3145 : // Update real dipole.
3146 0 : dip3real->iAcol = particles[dip2->iAcol].dips[dip2->iAcolLeg].
3147 0 : front()->iAcol;
3148 0 : dip3real->iCol = -( iAntiJun * 10 + 10 + 2);
3149 0 : dip3real->isAntiJun = true;
3150 0 : dip4real->iAcol = particles[dip1->iAcol].dips[dip1->iAcolLeg].
3151 0 : front()->iAcol;
3152 :
3153 : // Change the dipole connected to the antijunction.
3154 0 : dip3->iAcol = dip2->iAcol;
3155 0 : dip3->iAcolLeg = dip2->iAcolLeg;
3156 0 : dip3->isAntiJun = true;
3157 0 : dip3->iCol = -( iAntiJun * 10 + 10 + 2);
3158 0 : dip3->iColLeg = 0;
3159 :
3160 : // Change the dipole between the two gluons.
3161 0 : dip4->iAcol = dip1->iAcol;
3162 0 : dip4->iAcolLeg = dip1->iAcolLeg;
3163 :
3164 : // Store real dipole
3165 0 : particles[dip3->iAcol].dips[dip3->iAcolLeg].front() = dip3real;
3166 0 : particles[dip4->iAcol].dips[dip4->iAcolLeg].front() = dip4real;
3167 :
3168 : }
3169 0 : } else if (mode == 2) {
3170 0 : if (dip1 == dip3) {
3171 :
3172 : // Update real dipole.
3173 0 : dip4real->iAcol = particles[dip2->iAcol].dips[dip2->iAcolLeg].
3174 0 : front()->iAcol;
3175 0 : dip4real->iCol = -( iAntiJun * 10 + 10 + 2);
3176 0 : dip4real->isAntiJun = true;
3177 :
3178 : // Update active dipoles.
3179 0 : dip4->iAcol = dip2->iAcol;
3180 0 : dip4->iAcolLeg = dip2->iAcolLeg;
3181 0 : dip4->isAntiJun = true;
3182 0 : dip4->iCol = -( iAntiJun * 10 + 10 + 2);
3183 0 : dip4->iColLeg = 0;
3184 :
3185 : // Store real dipole
3186 0 : particles[dip4->iAcol].dips[dip4->iAcolLeg].front() = dip4real;
3187 :
3188 0 : } else {
3189 :
3190 : // Update real dipole.
3191 0 : dip4real->iAcol = particles[dip1->iAcol].dips[dip1->iAcolLeg].
3192 0 : front()->iAcol;
3193 0 : dip4real->iCol = -( iAntiJun * 10 + 10 + 2);
3194 0 : dip4real->isAntiJun = true;
3195 0 : dip3real->iAcol = particles[dip2->iAcol].dips[dip2->iAcolLeg].
3196 0 : front()->iAcol;
3197 :
3198 : // Change the dipole connected to the antijunction.
3199 0 : dip4->iAcol = dip1->iAcol;
3200 0 : dip4->iAcolLeg = dip1->iAcolLeg;
3201 0 : dip4->isAntiJun = true;
3202 0 : dip4->iCol = -( iAntiJun * 10 + 10 + 2);
3203 0 : dip4->iColLeg = 0;
3204 :
3205 : // Change the dipole between the two gluons.
3206 0 : dip3->iAcol = dip2->iAcol;
3207 0 : dip3->iAcolLeg = dip2->iAcolLeg;
3208 :
3209 : // Store real dipole
3210 0 : particles[dip3->iAcol].dips[dip3->iAcolLeg].front() = dip3real;
3211 0 : particles[dip4->iAcol].dips[dip4->iAcolLeg].front() = dip4real;
3212 : }
3213 : }
3214 :
3215 : // Dipoles connected to the junction.
3216 : // Update real dipoles.
3217 0 : particles[iCol1].dips[dip1->iColLeg].back()->iAcol = - (iJun * 10 + 10);
3218 0 : particles[iCol2].dips[dip2->iColLeg].back()->iAcol = - (iJun * 10 + 10 + 1);
3219 0 : particles[iCol1].dips[dip1->iColLeg].back()->isJun = true;
3220 0 : particles[iCol2].dips[dip2->iColLeg].back()->isJun = true;
3221 :
3222 : // Update active dipoles.
3223 0 : dip1->isJun = true;
3224 0 : dip2->isJun = true;
3225 0 : dip1->iAcol = - (iJun * 10 + 10);
3226 0 : dip2->iAcol = - (iJun * 10 + 10 + 1);
3227 0 : dip1->iAcolLeg = 0;
3228 0 : dip2->iAcolLeg = 0;
3229 :
3230 : // Update active dipoles for anti particles.
3231 : // Normally should only contain active dipoles once,
3232 : // only problem is if the two dipole ends are the same particle.
3233 : // Start by settings common dipoles.
3234 0 : for (int i = 0; i < int(particles[iAcol3].activeDips.size()); ++i)
3235 0 : if (particles[iAcol3].activeDips[i] == dip3) {
3236 0 : particles[iAcol3].activeDips[i] = dipoles[iActive2];
3237 0 : break;
3238 : }
3239 :
3240 0 : for (int i = 0; i < int(particles[iAcol4].activeDips.size()); ++i)
3241 0 : if (particles[iAcol4].activeDips[i] == dip4) {
3242 0 : particles[iAcol4].activeDips[i] = dipoles[iActive3];
3243 0 : break;
3244 : }
3245 :
3246 : // Depending on how the new string is connected, the active dipoles vary.
3247 0 : if (mode == 1) {
3248 0 : for (int i = 0; i < int(particles[iCol3].activeDips.size()); ++i)
3249 0 : if (particles[iCol3].activeDips[i] == dip3) {
3250 0 : particles[iCol3].activeDips[i] = dipoles[iActive1];
3251 0 : break;
3252 : }
3253 :
3254 0 : if (dip2 == dip4) {
3255 0 : for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i)
3256 0 : if (particles[iAcol1].activeDips[i] == dip1) {
3257 0 : particles[iAcol1].activeDips[i] = dip3;
3258 0 : break;
3259 : }
3260 0 : } else {
3261 0 : for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i)
3262 0 : if (particles[iAcol2].activeDips[i] == dip2) {
3263 0 : particles[iAcol2].activeDips[i] = dip3;
3264 0 : break;
3265 : }
3266 :
3267 0 : for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i)
3268 0 : if (particles[iAcol1].activeDips[i] == dip1) {
3269 0 : particles[iAcol1].activeDips[i] = dip4;
3270 0 : break;
3271 : }
3272 : }
3273 0 : } else if (mode == 2) {
3274 0 : for (int i = 0; i < int(particles[iCol4].activeDips.size()); ++i)
3275 0 : if (particles[iCol4].activeDips[i] == dip4) {
3276 0 : particles[iCol4].activeDips[i] = dipoles[iActive1];
3277 0 : break;
3278 : }
3279 :
3280 0 : if (dip1 == dip3) {
3281 0 : for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i)
3282 0 : if (particles[iAcol2].activeDips[i] == dip2) {
3283 0 : particles[iAcol2].activeDips[i] = dip4;
3284 0 : break;
3285 : }
3286 0 : } else {
3287 0 : for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i)
3288 0 : if (particles[iAcol1].activeDips[i] == dip1) {
3289 0 : particles[iAcol1].activeDips[i] = dip4;
3290 0 : break;
3291 : }
3292 :
3293 0 : for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i)
3294 0 : if (particles[iAcol2].activeDips[i] == dip2) {
3295 0 : particles[iAcol2].activeDips[i] = dip3;
3296 0 : break;
3297 : }
3298 : }
3299 : }
3300 :
3301 : // Add the junctions to the event.
3302 0 : junctions.push_back(Junction(1, oldCol1, oldCol2, newCol1));
3303 0 : if (mode == 0) junctions.push_back(Junction(2, newCol2, newCol3, newCol1));
3304 0 : else if (mode == 1)
3305 0 : junctions.push_back(Junction(2, newCol2, newCol3, oldCol3));
3306 0 : else if (mode == 2)
3307 0 : junctions.push_back(Junction(2, newCol2, newCol3, oldCol4));
3308 :
3309 : // Set junction information.
3310 0 : junctions[iJun].dipsOrig[0] =
3311 0 : particles[iCol1].dips[dip1->iColLeg].back();
3312 0 : junctions[iJun].dipsOrig[1] =
3313 0 : particles[iCol2].dips[dip2->iColLeg].back();
3314 0 : junctions[iJun].dipsOrig[2] = dipoles[iReal1];
3315 0 : junctions[iJun].dips[0] = dip1;
3316 0 : junctions[iJun].dips[1] = dip2;
3317 0 : junctions[iJun].dips[2] = dipoles[iActive1];
3318 :
3319 : // Set anti junction information.
3320 0 : junctions[iAntiJun].dips[0] = dipoles[iActive2];
3321 0 : junctions[iAntiJun].dips[1] = dipoles[iActive3];
3322 0 : junctions[iAntiJun].dipsOrig[0] = dipoles[iReal2];
3323 0 : junctions[iAntiJun].dipsOrig[1] = dipoles[iReal3];
3324 :
3325 0 : if (mode == 0) {
3326 0 : junctions[iAntiJun].dips[2] = dipoles[iActive1];
3327 0 : junctions[iAntiJun].dipsOrig[2] = dipoles[iReal1];
3328 0 : } else if (mode == 1) {
3329 0 : junctions[iAntiJun].dips[2] = dip3;
3330 0 : junctions[iAntiJun].dipsOrig[2] =
3331 0 : particles[dip3->iAcol].dips[dip3->iAcolLeg].front();
3332 0 : } else if (mode == 2) {
3333 0 : junctions[iAntiJun].dips[2] = dip4;
3334 0 : junctions[iAntiJun].dipsOrig[2] =
3335 0 : particles[dip4->iAcol].dips[dip4->iAcolLeg].front();
3336 0 : }
3337 :
3338 : // Make pseudo particles.
3339 0 : if (dip1->isActive && mDip(dip1) < m0)
3340 0 : makePseudoParticle(dip1, 110, true);
3341 0 : if (dip2->isActive && mDip(dip2) < m0)
3342 0 : makePseudoParticle(dip2, 110, true);
3343 0 : if (dip3->isActive && mDip(dip3) < m0)
3344 0 : makePseudoParticle(dip3, 110, true);
3345 0 : if (dip4->isActive && mDip(dip4) < m0)
3346 0 : makePseudoParticle(dip4, 110, true);
3347 :
3348 0 : if (dipoles[iActive1]->isActive && mDip(dipoles[iActive1]) < m0)
3349 0 : makePseudoParticle(dipoles[iActive1], 110, true);
3350 0 : if (dipoles[iActive2]->isActive && mDip(dipoles[iActive2]) < m0)
3351 0 : makePseudoParticle(dipoles[iActive2], 110, true);
3352 0 : if (dipoles[iActive3]->isActive && mDip(dipoles[iActive3]) < m0)
3353 0 : makePseudoParticle(dipoles[iActive3], 110, true);
3354 :
3355 : // Add new dipoles to usedDipoles.
3356 0 : usedDipoles.push_back(dipoles[iActive1]);
3357 0 : usedDipoles.push_back(dipoles[iActive2]);
3358 0 : usedDipoles.push_back(dipoles[iActive3]);
3359 :
3360 : // Done.
3361 0 : }
3362 :
3363 : //--------------------------------------------------------------------------
3364 :
3365 : // Change colour structure if it is three dipoles forming a junction system.
3366 :
3367 : void ColourReconnection::doTripleJunctionTrial(Event& event,
3368 : TrialReconnection& juncTrial) {
3369 :
3370 : // store information for easier acces.
3371 0 : ColourDipole* dip1 = juncTrial.dips[0];
3372 0 : ColourDipole* dip2 = juncTrial.dips[1];
3373 0 : ColourDipole* dip3 = juncTrial.dips[2];
3374 :
3375 : // Store indices.
3376 0 : int iCol1 = dip1->iCol;
3377 0 : int iCol2 = dip2->iCol;
3378 0 : int iCol3 = dip3->iCol;
3379 0 : int iAcol1 = dip1->iAcol;
3380 0 : int iAcol2 = dip2->iAcol;
3381 0 : int iAcol3 = dip3->iAcol;
3382 :
3383 : // Store colours
3384 0 : int oldCol1 = dip1->col;
3385 0 : int oldCol2 = dip2->col;
3386 0 : int oldCol3 = dip3->col;
3387 0 : int newCol1 = event.nextColTag();
3388 0 : int newCol2 = event.nextColTag();
3389 0 : int newCol3 = event.nextColTag();
3390 :
3391 : // Store new junction indices.
3392 0 : int iJun = junctions.size();
3393 0 : int iAntiJun = junctions.size() + 1;
3394 :
3395 : // Now make dipole between anti junction and iAcol1.
3396 : // Start by finding real iAcol.
3397 : int iAcol1real
3398 0 : = particles[iAcol1].dips[dip1->iAcolLeg].front()->iAcol;
3399 0 : dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10),
3400 0 : iAcol1real, dip1->colReconnection, false, true, false, true));
3401 0 : int iReal1 = dipoles.size() - 1;
3402 0 : particles[iAcol1].dips[dip1->iAcolLeg].front() = dipoles.back();
3403 :
3404 0 : dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10),
3405 0 : iAcol1, dip1->colReconnection, false, true));
3406 0 : dipoles.back()->iAcolLeg = dip1->iAcolLeg;
3407 0 : int iActive1 = dipoles.size() - 1;
3408 :
3409 : // Now make dipole between anti junction and iAcol2.
3410 : // Start by finding real iAcol2.
3411 : int iAcol2real
3412 0 : = particles[iAcol2].dips[dip2->iAcolLeg].front()->iAcol;
3413 0 : dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10 + 1),
3414 0 : iAcol2real, dip2->colReconnection, false, true, false, true));
3415 0 : int iReal2 = dipoles.size() - 1;
3416 0 : particles[iAcol2].dips[dip2->iAcolLeg].front() = dipoles.back();
3417 :
3418 0 : dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10 + 1),
3419 0 : iAcol2, dip2->colReconnection, false, true));
3420 0 : dipoles.back()->iAcolLeg = dip2->iAcolLeg;
3421 0 : int iActive2 = dipoles.size() - 1;
3422 :
3423 : // Now make dipole between anti junction and iAcol3.
3424 : // Start by finding real iAcol3.
3425 : int iAcol3real
3426 0 : = particles[iAcol3].dips[dip3->iAcolLeg].front()->iAcol;
3427 0 : dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 2),
3428 0 : iAcol3real, dip3->colReconnection, false, true, false, true));
3429 0 : int iReal3 = dipoles.size() - 1;
3430 0 : particles[iAcol3].dips[dip3->iAcolLeg].front() = dipoles.back();
3431 :
3432 0 : dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 2),
3433 0 : iAcol3, dip3->colReconnection, false, true));
3434 0 : dipoles.back()->iAcolLeg = dip3->iAcolLeg;
3435 0 : int iActive3 = dipoles.size() - 1;
3436 :
3437 : // Update already existing dipoles.
3438 :
3439 : // Update real dipoles.
3440 0 : particles[iCol1].dips[dip1->iColLeg].back()->iAcol = - (iJun * 10 + 10);
3441 0 : particles[iCol2].dips[dip2->iColLeg].back()->iAcol = - (iJun * 10 + 10 + 1);
3442 0 : particles[iCol3].dips[dip3->iColLeg].back()->iAcol = - (iJun * 10 + 10 + 2);
3443 0 : particles[iCol1].dips[dip1->iColLeg].back()->isJun = true;
3444 0 : particles[iCol2].dips[dip2->iColLeg].back()->isJun = true;
3445 0 : particles[iCol3].dips[dip3->iColLeg].back()->isJun = true;
3446 :
3447 : // Update active dipoles.
3448 0 : dip1->isJun = true;
3449 0 : dip2->isJun = true;
3450 0 : dip3->isJun = true;
3451 0 : dip1->iAcol = - (iJun * 10 + 10);
3452 0 : dip2->iAcol = - (iJun * 10 + 10 + 1);
3453 0 : dip3->iAcol = - (iJun * 10 + 10 + 2);
3454 0 : dip1->iAcolLeg = 0;
3455 0 : dip2->iAcolLeg = 0;
3456 0 : dip3->iAcolLeg = 0;
3457 :
3458 : // Update active dipoles for anti particles.
3459 0 : for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i)
3460 0 : if (particles[iAcol1].activeDips[i] == dip1)
3461 0 : particles[iAcol1].activeDips[i] = dipoles[iActive1];
3462 0 : for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i)
3463 0 : if (particles[iAcol2].activeDips[i] == dip2)
3464 0 : particles[iAcol2].activeDips[i] = dipoles[iActive2];
3465 0 : for (int i = 0; i < int(particles[iAcol3].activeDips.size()); ++i)
3466 0 : if (particles[iAcol3].activeDips[i] == dip3)
3467 0 : particles[iAcol3].activeDips[i] = dipoles[iActive3];
3468 :
3469 : // Add the junctions to the event.
3470 0 : junctions.push_back(Junction(1, oldCol1, oldCol2, oldCol3));
3471 0 : junctions.push_back(Junction(2, newCol1, newCol3, newCol3));
3472 :
3473 : // Update junction ends.
3474 0 : junctions[iJun].dipsOrig[0] =
3475 0 : particles[iCol1].dips[dip1->iColLeg].back();
3476 0 : junctions[iJun].dipsOrig[1] =
3477 0 : particles[iCol2].dips[dip2->iColLeg].back();
3478 0 : junctions[iJun].dipsOrig[2] =
3479 0 : particles[iCol3].dips[dip3->iColLeg].back();
3480 0 : junctions[iJun].dips[0] = dip1;
3481 0 : junctions[iJun].dips[1] = dip2;
3482 0 : junctions[iJun].dips[2] = dip3;
3483 :
3484 : // Update the anti junction.
3485 0 : junctions[iAntiJun].dips[0] = dipoles[iActive1];
3486 0 : junctions[iAntiJun].dips[1] = dipoles[iActive2];
3487 0 : junctions[iAntiJun].dips[2] = dipoles[iActive3];
3488 0 : junctions[iAntiJun].dipsOrig[0] = dipoles[iReal1];
3489 0 : junctions[iAntiJun].dipsOrig[1] = dipoles[iReal2];
3490 0 : junctions[iAntiJun].dipsOrig[2] = dipoles[iReal3];
3491 :
3492 : // Make pseudo particles if needed.
3493 0 : if (dip1->isActive && mDip(dip1) < m0)
3494 0 : makePseudoParticle(dip1, 110, true);
3495 0 : if (dip2->isActive && mDip(dip2) < m0)
3496 0 : makePseudoParticle(dip2, 110, true);
3497 0 : if (dip3->isActive && mDip(dip3) < m0)
3498 0 : makePseudoParticle(dip3, 110, true);
3499 :
3500 0 : if (dipoles[iActive1]->isActive && mDip(dipoles[iActive1]) < m0)
3501 0 : makePseudoParticle(dipoles[iActive1], 110, true);
3502 0 : if (dipoles[iActive2]->isActive && mDip(dipoles[iActive2]) < m0)
3503 0 : makePseudoParticle(dipoles[iActive2], 110, true);
3504 0 : if (dipoles[iActive3]->isActive && mDip(dipoles[iActive3]) < m0)
3505 0 : makePseudoParticle(dipoles[iActive3], 110, true);
3506 :
3507 : // Add to newly created dipoles to used dipoles.
3508 0 : usedDipoles.push_back(dipoles[iActive1]);
3509 0 : usedDipoles.push_back(dipoles[iActive2]);
3510 0 : usedDipoles.push_back(dipoles[iActive3]);
3511 :
3512 : // Done.
3513 0 : }
3514 :
3515 : //--------------------------------------------------------------------------
3516 :
3517 : // Allow colour reconnections by moving gluons from their current location
3518 : // to another colour line. Also optionally flip two colour chains.
3519 :
3520 : bool ColourReconnection::reconnectMove( Event& event, int oldSize) {
3521 :
3522 : // Create or reset arrays to prepare for the new event analysis.
3523 0 : vector<int> iGlu;
3524 0 : iReduceCol.resize( event.size() );
3525 0 : iExpandCol.clear();
3526 0 : map<int, int> colMap, acolMap;
3527 0 : map<int, int>::iterator colM, acolM;
3528 0 : vector<InfoGluonMove> infoGM;
3529 :
3530 : // Temporary variables.
3531 0 : int iNow = 0;
3532 0 : int colNow = 0;
3533 0 : int acolNow = 0;
3534 : int iColNow = 0;
3535 : int iAcolNow = 0;
3536 0 : int col2Now = 0;
3537 : int iCol2Now = 0;
3538 : int iAcol2Now = 0;
3539 : double lambdaRefNow = 0.;
3540 : double dLambdaNow = 0.;
3541 :
3542 : // Loop over all final particles. Store (fraction of) gluons to move.
3543 0 : for (int i = oldSize; i < event.size(); ++i) if (event[i].isFinal()) {
3544 0 : if (flipMode < 3 && event[i].id() == 21 && rndmPtr->flat() < fracGluon)
3545 0 : iGlu.push_back(i);
3546 :
3547 : // Store location of all colour and anticolour particles and indices.
3548 0 : if (event[i].col() > 0 || event[i].acol() > 0) {
3549 0 : iReduceCol[i] = iExpandCol.size();
3550 0 : iExpandCol.push_back(i);
3551 0 : if (event[i].col() > 0) colMap[event[i].col()] = i;
3552 0 : if (event[i].acol() > 0) acolMap[event[i].acol()] = i;
3553 : }
3554 : }
3555 :
3556 : // Erase (anti)colours for (anti)junctions and skip adjacent gluons.
3557 0 : for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
3558 0 : if (event.kindJunction(iJun) == 1) {
3559 0 : for (int j = 0; j < 3; ++j) {
3560 0 : int jCol = event.colJunction( iJun, j);
3561 0 : for (colM = colMap.begin(); colM != colMap.end(); ++colM)
3562 0 : if (colM->first == jCol) {
3563 0 : colMap.erase( colM);
3564 0 : break;
3565 : }
3566 0 : for (int iG = 0; iG < int(iGlu.size()); ++iG)
3567 0 : if (event[iGlu[iG]].col() == jCol) {
3568 0 : iGlu.erase(iGlu.begin() + iG);
3569 0 : break;
3570 : }
3571 : }
3572 0 : } else if (event.kindJunction(iJun) == 2) {
3573 0 : for (int j = 0; j < 3; ++j) {
3574 0 : int jCol = event.colJunction( iJun, j);
3575 0 : for (acolM = acolMap.begin(); acolM != acolMap.end(); ++acolM)
3576 0 : if (acolM->first == jCol) {
3577 0 : acolMap.erase( acolM);
3578 0 : break;
3579 : }
3580 0 : for (int iG = 0; iG < int(iGlu.size()); ++iG)
3581 0 : if (event[iGlu[iG]].acol() == jCol) {
3582 0 : iGlu.erase(iGlu.begin() + iG);
3583 0 : break;
3584 : }
3585 : }
3586 0 : }
3587 : }
3588 :
3589 : // Error checks.
3590 0 : int nGlu = iGlu.size();
3591 0 : int nCol = colMap.size();
3592 0 : if (int(acolMap.size()) != nCol) {
3593 0 : infoPtr->errorMsg("Error in MBReconUserHooks: map sizes do not match");
3594 0 : return false;
3595 : }
3596 0 : colM = colMap.begin();
3597 0 : acolM = acolMap.begin();
3598 0 : for (int iCol = 0; iCol < nCol; ++iCol) {
3599 0 : if (colM->first != acolM->first) {
3600 0 : infoPtr->errorMsg("Error in MBReconUserHooks: map elements"
3601 : " do not match");
3602 0 : return false;
3603 : }
3604 0 : ++colM;
3605 0 : ++acolM;
3606 : }
3607 :
3608 : // Calculate and tabulate lambda between any pair of coloured partons.
3609 0 : nColMove = iExpandCol.size();
3610 0 : lambdaijMove.resize( pow2(nColMove) );
3611 0 : for (int iAC = 0; iAC < nColMove - 1; ++iAC) {
3612 0 : int i = iExpandCol[iAC];
3613 0 : for (int jAC = iAC + 1; jAC < nColMove; ++jAC) {
3614 0 : int j = iExpandCol[jAC];
3615 0 : lambdaijMove[nColMove * iAC + jAC]
3616 0 : = log(1. + m2( event[i], event[j]) / m2Lambda);
3617 : }
3618 : }
3619 :
3620 : // Set up initial possible gluon moves with lambda gains/losses.
3621 0 : for (int iG = 0; iG < nGlu; ++iG) {
3622 :
3623 : // Gluon and its current neighbours.
3624 0 : iNow = iGlu[iG];
3625 0 : colNow = event[iNow].col();
3626 0 : acolNow = event[iNow].acol();
3627 0 : iColNow = acolMap[colNow];
3628 0 : iAcolNow = colMap[acolNow];
3629 :
3630 : // Addition to Lambda of gluon in current position.
3631 0 : lambdaRefNow = lambda123Move( iNow, iColNow, iAcolNow);
3632 :
3633 : // Loop over all colour lines where gluon could be inserted.
3634 0 : for (colM = colMap.begin(); colM != colMap.end(); ++colM) {
3635 0 : col2Now = colM->first;
3636 0 : iCol2Now = colMap[col2Now];
3637 0 : iAcol2Now = acolMap[col2Now];
3638 :
3639 : // Addition to total Lambda if gluon moved to be inserted on line.
3640 0 : dLambdaNow = (iCol2Now == iNow || iAcol2Now == iNow
3641 0 : || iColNow == iAcolNow) ? 2e4
3642 0 : : lambda123Move( iNow, iCol2Now, iAcol2Now) - lambdaRefNow;
3643 :
3644 : // Add new container for gluon and colour line information.
3645 0 : infoGM.push_back( InfoGluonMove( iNow, colNow, acolNow, iColNow,
3646 0 : iAcolNow, col2Now, iCol2Now, iAcol2Now, lambdaRefNow, dLambdaNow ));
3647 : }
3648 : }
3649 0 : int nPair = infoGM.size();
3650 :
3651 : // Keep on looping over moves until no further negative dLambda.
3652 0 : for ( int iMove = 0; iMove < nGlu; ++iMove) {
3653 : int iPairMin = -1;
3654 : double dLambdaMin = 1e4;
3655 :
3656 : // Find lowest dLambda.
3657 0 : for (int iPair = 0; iPair < nPair; ++iPair)
3658 0 : if (infoGM[iPair].dLambda < dLambdaMin) {
3659 : iPairMin = iPair;
3660 0 : dLambdaMin = infoGM[iPair].dLambda;
3661 0 : }
3662 :
3663 : // Break if no shift below upper limit found.
3664 0 : if (dLambdaMin > -dLambdaCut) break;
3665 :
3666 : // Partons and colours involved in move.
3667 0 : InfoGluonMove& selSM = infoGM[iPairMin];
3668 0 : int i1Sel = selSM.i1;
3669 0 : int iCol1Sel = selSM.iCol1;
3670 0 : int iAcol1Sel = selSM.iAcol1;
3671 0 : int iCol2Sel = selSM.iCol2;
3672 0 : int iAcol2Sel = selSM.iAcol2;
3673 0 : int iCol2Mod[3] = { iAcol1Sel , i1Sel , iCol2Sel };
3674 0 : int col2Mod[3] = { selSM.col1, selSM.col2, selSM.acol1};
3675 :
3676 : // Remove gluon from old colour line and insert on new colour line.
3677 0 : for (int i = 0; i < 3; ++i) {
3678 0 : event[ iCol2Mod[i] ].col( col2Mod[i] );
3679 0 : colMap[ col2Mod[i] ] = iCol2Mod[i];
3680 : }
3681 :
3682 : // Update info for partons with new colors.
3683 : int i1Now = 0;
3684 : bool doUpdate = false;
3685 0 : for (int iPair = 0; iPair < nPair; ++iPair) {
3686 0 : InfoGluonMove& tmpSM = infoGM[iPair];
3687 0 : if (tmpSM.i1 != i1Now) {
3688 : i1Now = tmpSM.i1;
3689 : doUpdate = false;
3690 0 : if (i1Now == i1Sel || i1Now == iCol1Sel || i1Now == iAcol1Sel
3691 0 : || i1Now == iCol2Sel || i1Now == iAcol2Sel) {
3692 0 : colNow = event[i1Now].col();
3693 0 : acolNow = event[i1Now].acol();
3694 0 : iColNow = acolMap[colNow];
3695 0 : iAcolNow = colMap[acolNow];
3696 0 : lambdaRefNow = lambda123Move( i1Now, iColNow, iAcolNow);
3697 : doUpdate = true;
3698 0 : }
3699 : }
3700 0 : if (doUpdate) {
3701 0 : tmpSM.col1 = colNow;
3702 0 : tmpSM.acol1 = acolNow;
3703 0 : tmpSM.iCol1 = iColNow;
3704 0 : tmpSM.iAcol1 = iAcolNow;
3705 0 : tmpSM.lambdaRef = lambdaRefNow;
3706 0 : }
3707 : }
3708 :
3709 : // Update info on dLambda for affected particles and colour lines.
3710 0 : for (int iPair = 0; iPair < nPair; ++iPair) {
3711 0 : InfoGluonMove& tmpSM = infoGM[iPair];
3712 : int iMod = -1;
3713 0 : for (int i = 0; i < 3; ++i) if (tmpSM.col2 == col2Mod[i]) iMod = i;
3714 0 : if (iMod > -1) tmpSM.iCol2 = iCol2Mod[iMod];
3715 0 : if (tmpSM.i1 == i1Sel || tmpSM.i1 == iCol1Sel || tmpSM.i1 == iAcol1Sel
3716 0 : || tmpSM.i1 == iCol2Sel || tmpSM.i1 == iAcol2Sel || iMod > -1)
3717 0 : tmpSM.dLambda = (tmpSM.iCol2 == tmpSM.i1 || tmpSM.iAcol2 == tmpSM.i1
3718 0 : || tmpSM.iCol1 == tmpSM.iAcol1) ? 2e4
3719 0 : : lambda123Move( tmpSM.i1, tmpSM.iCol2, tmpSM.iAcol2)
3720 0 : - tmpSM.lambdaRef;
3721 : }
3722 :
3723 : // End of loop over gluon shifting.
3724 0 : }
3725 :
3726 : // Done if no flip.
3727 0 : if (flipMode == 0) return true;
3728 :
3729 : // Array with colour lines, and where each line begins and ends.
3730 0 : vector<int> iTmpFlip, iVecFlip, iBegFlip, iEndFlip;
3731 :
3732 : // Variables for minimum search.
3733 : int i1c, i1a, i2c, i2a, i1cMin, i1aMin, i2cMin, i2aMin, iSMin;
3734 : double dLambdaFlip, dLambdaFlipMin;
3735 0 : vector<InfoGluonMove> flipMin;
3736 :
3737 : // Grab all colour ends.
3738 0 : for (int i = oldSize; i < event.size(); ++i)
3739 0 : if (event[i].isFinal() && event[i].col() > 0 && event[i].acol() == 0) {
3740 0 : iTmpFlip.clear();
3741 0 : iTmpFlip.push_back( i);
3742 :
3743 : // Step through colour neighbours to catch system.
3744 0 : iNow = i;
3745 0 : acolM = acolMap.find( event[iNow].col() );
3746 : bool foundEnd = false;
3747 0 : while (acolM != acolMap.end()) {
3748 0 : iNow = acolM->second;
3749 0 : iTmpFlip.push_back( iNow);
3750 0 : if (event[iNow].col() == 0) {
3751 : foundEnd = true;
3752 0 : break;
3753 : }
3754 0 : acolM = acolMap.find( event[iNow].col() );
3755 : }
3756 :
3757 : // Store acceptable system, optionally including junction legs.
3758 0 : if (foundEnd || flipMode == 2 || flipMode == 4) {
3759 0 : iBegFlip.push_back( iVecFlip.size());
3760 0 : for (int j = 0; j < int(iTmpFlip.size()); ++j)
3761 0 : iVecFlip.push_back( iTmpFlip[j]);
3762 0 : iEndFlip.push_back( iVecFlip.size());
3763 : }
3764 0 : }
3765 :
3766 : // Optionally search for antijunction legs: grab all anticolour ends.
3767 0 : if (flipMode == 2 || flipMode == 4)
3768 0 : for (int i = oldSize; i < event.size(); ++i)
3769 0 : if (event[i].isFinal() && event[i].acol() > 0 && event[i].col() == 0) {
3770 0 : iTmpFlip.clear();
3771 0 : iTmpFlip.push_back( i);
3772 :
3773 : // Step through anticolour neighbours to catch system.
3774 0 : iNow = i;
3775 0 : colM = colMap.find( event[iNow].acol() );
3776 : bool foundEnd = false;
3777 0 : while (colM != colMap.end()) {
3778 0 : iNow = colM->second;
3779 0 : iTmpFlip.push_back( iNow);
3780 0 : if (event[iNow].acol() == 0) {
3781 : foundEnd = true;
3782 0 : break;
3783 : }
3784 0 : colM = colMap.find( event[iNow].acol() );
3785 : }
3786 :
3787 : // Store acceptable system, but do not doublecount q - (n g) - qbar.
3788 0 : if (!foundEnd) {
3789 0 : iBegFlip.push_back( iVecFlip.size());
3790 0 : for (int j = 0; j < int(iTmpFlip.size()); ++j)
3791 0 : iVecFlip.push_back( iTmpFlip[j]);
3792 0 : iEndFlip.push_back( iVecFlip.size());
3793 : }
3794 0 : }
3795 :
3796 : // Loop through all system pairs.
3797 0 : int nSysFlip = iBegFlip.size();
3798 0 : for (int iSys1 = 0; iSys1 < nSysFlip - 1; ++iSys1)
3799 0 : if (iBegFlip[iSys1] >= 0)
3800 0 : for (int iSys2 = iSys1 + 1; iSys2 < nSysFlip; ++iSys2)
3801 0 : if (iBegFlip[iSys2] >= 0) {
3802 : i1cMin = 0;
3803 : i1aMin = 0;
3804 : i2cMin = 0;
3805 : i2aMin = 0;
3806 : dLambdaFlipMin = 1e4;
3807 :
3808 : // Loop through all possible flip locations for a pair.
3809 0 : for (int j1 = iBegFlip[iSys1]; j1 < iEndFlip[iSys1] - 1; ++j1)
3810 0 : for (int j2 = iBegFlip[iSys2]; j2 < iEndFlip[iSys2] - 1; ++j2) {
3811 0 : i1c = iVecFlip[j1];
3812 0 : i1a = iVecFlip[j1 + 1];
3813 0 : i2c = iVecFlip[j2];
3814 0 : i2a = iVecFlip[j2 + 1];
3815 0 : dLambdaFlip = lambda12Move( i1c, i2a) + lambda12Move( i2c, i1a)
3816 0 : - lambda12Move( i1c, i1a) - lambda12Move( i2c, i2a);
3817 0 : if (dLambdaFlip < dLambdaFlipMin) {
3818 : i1cMin = i1c;
3819 : i1aMin = i1a;
3820 : i2cMin = i2c;
3821 : i2aMin = i2a;
3822 : dLambdaFlipMin = dLambdaFlip;
3823 0 : }
3824 : }
3825 :
3826 : // Store possible flips if low enough dLambdaMin.
3827 0 : if (dLambdaFlipMin < -dLambdaCut) flipMin.push_back( InfoGluonMove(
3828 0 : iSys1, iSys2, i1cMin, i1aMin, i2cMin, i2aMin, dLambdaFlipMin) );
3829 0 : }
3830 0 : int flipSize = flipMin.size();
3831 :
3832 : // Search for lowest possible flip among unused systems.
3833 0 : for (int iFlip = 0; iFlip < min( nSysFlip / 2, flipSize); ++iFlip) {
3834 : iSMin = -1;
3835 : dLambdaFlipMin = 1e4;
3836 0 : for (int iSys12 = 0; iSys12 < flipSize; ++iSys12)
3837 0 : if (flipMin[iSys12].i1 >= 0 && flipMin[iSys12].dLambda < dLambdaFlipMin) {
3838 : iSMin = iSys12;
3839 0 : dLambdaFlipMin = flipMin[iSys12].dLambda;
3840 0 : }
3841 :
3842 : // Do flip. Mark flipped systems.
3843 0 : if (iSMin >= 0) {
3844 0 : InfoGluonMove& flipNow = flipMin[iSMin];
3845 0 : int iS1 = flipNow.i1;
3846 0 : int iS2 = flipNow.i2;
3847 0 : event[ flipNow.iAcol1 ].acol( event[flipNow.iCol2].col() );
3848 0 : event[ flipNow.iAcol2 ].acol( event[flipNow.iCol1].col() );
3849 0 : for (int iSys12 = 0; iSys12 < flipSize; ++iSys12)
3850 0 : if ( flipMin[iSys12].i1 == iS1 || flipMin[iSys12].i1 == iS2
3851 0 : || flipMin[iSys12].i2 == iS1 || flipMin[iSys12].i2 == iS2)
3852 0 : flipMin[iSys12].i1 = -1;
3853 : }
3854 0 : else break;
3855 : }
3856 :
3857 : // Done.
3858 : return true;
3859 :
3860 0 : }
3861 :
3862 : //--------------------------------------------------------------------------
3863 :
3864 : // Common code for the SK I and SK II models for WW/ZZ/WZ systems.
3865 :
3866 : bool ColourReconnection::reconnectTypeCommon( Event& event, int ) {
3867 :
3868 : // Make storage containers. Check that at least two parton systems.
3869 0 : vector<vector< ColourDipole> > dips;
3870 0 : int iBosons[2];
3871 0 : Vec4 decays[2];
3872 0 : if (partonSystemsPtr->sizeSys() < 2) {
3873 0 : infoPtr->errorMsg("Error in ColourReconnection::reconnectTypeCommon: "
3874 : "expect at least two parton systems");
3875 0 : return false;
3876 : }
3877 :
3878 : // Find the dipoles connected to their respective resonance decays.
3879 0 : for (int i = 0; i < 2; ++i) {
3880 0 : dips.push_back(vector<ColourDipole>());
3881 0 : int iSys = partonSystemsPtr->sizeSys() - i - 1;
3882 0 : for (int j = 0; j < partonSystemsPtr->sizeOut(iSys); ++j) {
3883 0 : int iPar = partonSystemsPtr->getOut(iSys, j);
3884 :
3885 : // Find decayed boson. Only need to do it once.
3886 0 : if (j == 0) {
3887 0 : int iMot = event[iPar].mother1();
3888 0 : while (iMot != 0 && event[iMot].idAbs() != 23
3889 0 : && event[iMot].idAbs() != 24) iMot = event[iMot].mother1();
3890 0 : if (iMot == 0) {
3891 0 : infoPtr->errorMsg("Error in ColourReconnection::reconnectTypeCommon:"
3892 : " Not a resonance decay of a W/Z");
3893 0 : return false;
3894 : }
3895 0 : iBosons[i] = iMot;
3896 0 : }
3897 :
3898 : // Pick up all dipoles by starting from colour end.
3899 0 : int col = event[iPar].col();
3900 0 : if (col == 0) continue;
3901 0 : for (int k = 0; k < partonSystemsPtr->sizeOut(iSys); ++k) {
3902 0 : int iAntiPar = partonSystemsPtr->getOut(iSys, k);
3903 0 : if (event[iAntiPar].acol() == col) {
3904 0 : dips.back().push_back(ColourDipole(col, iPar, iAntiPar));
3905 0 : break;
3906 : }
3907 0 : }
3908 0 : }
3909 0 : }
3910 :
3911 : // Boost system to W+W- rest frame.
3912 0 : Vec4 boost = event[iBosons[0]].p() + event[iBosons[1]].p();
3913 0 : for (int i = 1; i < event.size(); ++i) event[i].bstback(boost);
3914 :
3915 : // Find the time and position of the decay vertices.
3916 0 : for (int i = 0; i < 2; ++i) {
3917 0 : int iBoson = iBosons[i];
3918 0 : double mBoson = particleDataPtr->m0(event[iBoson].idAbs());
3919 0 : double gammaBoson = particleDataPtr->mWidth(event[iBoson].idAbs());
3920 0 : double mReal = event[iBoson].mCalc();
3921 0 : decays[i][0] = -HBAR * log(rndmPtr->flat()) * event[iBoson].e() /
3922 0 : sqrt(pow2(pow2(mBoson) - pow2(mReal)) + pow2(gammaBoson * pow2(mReal)
3923 0 : / mBoson));
3924 0 : for (int j = 1; j < 4; ++j)
3925 0 : decays[i][j] = event[iBoson].p()[j]/ event[iBoson].e() * decays[i][0];
3926 0 : }
3927 :
3928 : // Find the possible reconnections, depeding on choice of model.
3929 0 : map<double,pair<int,int> > reconnections;
3930 0 : if (reconnectMode == 3)
3931 0 : reconnections = reconnectTypeI(event, dips, decays);
3932 0 : else if (reconnectMode == 4)
3933 0 : reconnections = reconnectTypeII(event, dips, decays);
3934 :
3935 : // Carry out the reconnections.
3936 0 : vector<bool> used1(dips[0].size(), false), used2(dips[1].size(), false);
3937 0 : for (map<double,pair<int,int> >::iterator it=reconnections.begin();
3938 0 : it != reconnections.end(); ++it) {
3939 :
3940 : // Check if any of the dipoles already reconnected.
3941 0 : if (used1[it->second.first] || used2[it->second.second]) continue;
3942 :
3943 : // Do the reconnection.
3944 0 : int iAcol1 = dips[0][it->second.first].iAcol;
3945 0 : int iAcol2 = dips[1][it->second.second].iAcol;
3946 0 : event.copy(iAcol1, 79);
3947 0 : event.back().acol(event[iAcol2].acol());
3948 0 : event.copy(iAcol2, 79);
3949 0 : event.back().acol(event[iAcol1].acol());
3950 :
3951 : // Mark dipoles as used.
3952 0 : used1[it->second.first] = true;
3953 0 : used2[it->second.second] = true;
3954 :
3955 : // If only a single reconnection is wanted, break the loop.
3956 0 : if (singleReconOnly) break;
3957 0 : }
3958 :
3959 : // Boost system back to origianl rest frame.
3960 0 : for (int i = 1; i < event.size(); ++i) event[i].bst(boost);
3961 :
3962 : // Done.
3963 : return true;
3964 0 : }
3965 :
3966 : //--------------------------------------------------------------------------
3967 :
3968 : // The SK I model for WW/ZZ/WZ systems.
3969 :
3970 : map<double,pair<int,int> > ColourReconnection::reconnectTypeI( Event& event,
3971 : vector<vector<ColourDipole> > &dips, Vec4 decays[2]) {
3972 :
3973 : // Make storage containers.
3974 0 : multimap<double,pair<int,int> > reconnections;
3975 :
3976 : // Velocity stored as: vx vy, vz , gamma (uses Vec4 for easy storage).
3977 0 : vector<vector<Vec4> > vel;
3978 : // string direction.
3979 0 : vector<vector<Vec4> > dir;
3980 :
3981 : // Calculate velocities.
3982 0 : for (int i = 0; i < 2; ++i) {
3983 0 : vel.push_back(vector<Vec4>(dips[i].size(),Vec4()));
3984 0 : dir.push_back(vector<Vec4>(dips[i].size(),Vec4()));
3985 0 : for (int j = 0; j < int(dips[i].size()); ++j) {
3986 :
3987 : // Calculate sum of momenta.
3988 0 : double pSumCol = event[dips[i][j].iCol].pAbs();
3989 0 : double pSumAcol = event[dips[i][j].iAcol].pAbs();
3990 :
3991 : // Velocities and directions.
3992 0 : for (int k = 1; k < 4; ++k) {
3993 0 : vel[i][j][k] = 0.5 *(event[dips[i][j].iCol].p()[k] / pSumCol
3994 0 : + event[dips[i][j].iAcol].p()[k] / pSumAcol);
3995 0 : dir[i][j][k] = event[dips[i][j].iCol].p()[k] / pSumCol
3996 0 : - event[dips[i][j].iAcol].p()[k] / pSumAcol;
3997 : }
3998 0 : vel[i][j][0] = 1/sqrt(1 - vel[i][j].pAbs2());
3999 0 : dir[i][j] /= dir[i][j].pAbs();
4000 :
4001 : }
4002 : }
4003 :
4004 : // Select random point.
4005 : int nPoints = 100;
4006 : double sumW = 0;
4007 0 : for (int i = 0; i < nPoints; ++i) {
4008 0 : Vec4 x;
4009 0 : double r = sqrt(- log(rndmPtr->flat()));
4010 0 : double phi = 2 * M_PI * rndmPtr->flat();
4011 0 : x[1] = blowR * rHadron * r * cos(phi);
4012 0 : x[2] = blowR * rHadron * r * sin(phi);
4013 0 : r = sqrt(- log(rndmPtr->flat()));
4014 0 : phi = 2 * M_PI * rndmPtr->flat();
4015 0 : x[3] = blowR * rHadron * r * cos(phi);
4016 0 : x[0] = max(decays[0][0], decays[1][0])
4017 0 : + blowT * sqrt(0.5) * tfrag * r * abs(sin(phi));
4018 0 : if (x.m2Calc() < 0) continue;
4019 :
4020 : // Find weight of trial point.
4021 0 : double weightTrial = exp(-x.pAbs2() / pow2(blowR*rHadron))
4022 0 : * exp( -2. * pow2(x[0] - max(decays[0][0],decays[1][0]))
4023 0 : / pow2(blowT * tfrag) );
4024 :
4025 : // Find max weight and their indices.
4026 0 : double maxWeights[2] = {0,0}; // {1E-10,1E-10};
4027 0 : int maxIndices[2] = {-1,-1};
4028 :
4029 : // Loop over W decays.
4030 0 : for (int j = 0;j < 2;++j) {
4031 :
4032 : // Calculate the difference between decay point and random point.
4033 0 : Vec4 xd = x - decays[j];
4034 :
4035 : // Loop over strings.
4036 0 : for (int k = 0; k < int(dips[j].size()); ++k) {
4037 :
4038 : // Boost to rest frame of string.
4039 0 : double gam = vel[j][k][0];
4040 0 : double dotProd = dot3(xd,vel[j][k]);
4041 0 : double xBoost = gam * (gam * dotProd/ (1. + gam) - xd[0]);
4042 0 : Vec4 xb = xd + xBoost * vel[j][k];
4043 0 : xb[0] = gam * (xd[0] - dotProd);
4044 :
4045 : // Check that position is reachable.
4046 0 : if (xb.m2Calc() < 0) continue;
4047 :
4048 : // Find and store largest weight.
4049 0 : double weight = exp( -(xb.pAbs2() - pow2(dot3(xb,dir[j][k])))
4050 0 : / (2 * pow2(rHadron)) )
4051 0 : * exp( -(pow2(xb[0]) - pow2(dot3(xb, dir[j][k]))) / pow2(tfrag) );
4052 0 : if (weight > maxWeights[j]) {
4053 0 : maxWeights[j] = weight;
4054 0 : maxIndices[j] = k;
4055 0 : }
4056 0 : }
4057 0 : }
4058 :
4059 : // Store weight.
4060 0 : if (maxIndices[0] != -1 && maxIndices[1] != -1) {
4061 :
4062 : // Check if the new configuration lowers the lambda measure.
4063 0 : if (lowerLambdaOnly) {
4064 0 : double oldLambda = stringLength.getStringLength(event,
4065 0 : dips[0][maxIndices[0]].iCol, dips[0][maxIndices[0]].iAcol) +
4066 0 : stringLength.getStringLength(event, dips[1][maxIndices[1]].iCol,
4067 0 : dips[1][maxIndices[1]].iAcol);
4068 0 : double newLambda = stringLength.getStringLength(event,
4069 0 : dips[0][maxIndices[0]].iCol, dips[1][maxIndices[1]].iAcol) +
4070 0 : stringLength.getStringLength(event, dips[1][maxIndices[1]].iCol,
4071 0 : dips[0][maxIndices[0]].iAcol);
4072 0 : if (oldLambda < newLambda) continue;
4073 0 : }
4074 :
4075 : // Save weights.
4076 0 : double weight = maxWeights[0] * maxWeights[1] / weightTrial;
4077 0 : reconnections.insert(make_pair(weight,
4078 0 : make_pair(maxIndices[0], maxIndices[1])));
4079 0 : sumW += weight;
4080 0 : }
4081 0 : }
4082 :
4083 :
4084 : // check whether to do any reconnections.
4085 0 : map<double,pair<int,int> > singleRecon;
4086 0 : double result = pow3(blowR) * blowT * sumW/nPoints;
4087 : // No reconections.
4088 0 : if (1 - exp(-kI * result) < rndmPtr->flat()) {
4089 0 : singleRecon.clear();
4090 0 : return singleRecon;
4091 : // Find reconnection.
4092 : } else {
4093 0 : double rSum = sumW * rndmPtr->flat();
4094 0 : for (map<double,pair<int,int> >::iterator it = reconnections.begin();
4095 0 : it != reconnections.end(); ++it) {
4096 0 : rSum -= it->first;
4097 0 : if (rSum < 0.) {
4098 0 : singleRecon.insert(make_pair(1., it->second));
4099 0 : return singleRecon;
4100 : }
4101 : }
4102 :
4103 : // If no solution was found (shoult not happen) do no reconnections.
4104 0 : singleRecon.clear();
4105 0 : return singleRecon;
4106 : }
4107 0 : }
4108 :
4109 : //--------------------------------------------------------------------------
4110 :
4111 : // The SK II model for WW/ZZ/WZ systems.
4112 :
4113 : map<double,pair<int,int> > ColourReconnection::reconnectTypeII( Event& event,
4114 : vector<vector<ColourDipole> > &dips, Vec4 decays[2]) {
4115 :
4116 : // Make storage containers.
4117 0 : map<double,pair<int,int> > reconnections;
4118 :
4119 : // Find dipole velocities.
4120 0 : for (int i = 0;i < int(dips[0].size()); ++i) {
4121 0 : for (int j = 0;j < int(dips[1].size()); ++j) {
4122 0 : Vec4 v1,v2,u1,u2;
4123 0 : v1 = event[dips[0][i].iCol].p() / event[dips[0][i].iCol].e();
4124 0 : v2 = event[dips[0][i].iAcol].p() / event[dips[0][i].iAcol].e();
4125 0 : u1 = event[dips[1][j].iCol].p() / event[dips[1][j].iCol].e();
4126 0 : u2 = event[dips[1][j].iAcol].p() / event[dips[1][j].iAcol].e();
4127 :
4128 : // Begin setup to solve system of equations.
4129 0 : vector<vector<double> > matUpper, matLower;
4130 0 : for (int k = 0; k < 3; ++k) {
4131 0 : matUpper.push_back(vector<double>(3,0));
4132 0 : matLower.push_back(vector<double>(3,0));
4133 : }
4134 :
4135 : // Insert in matrix, beware of different convention in index.
4136 0 : for (int k = 0;k < 3; ++k) {
4137 0 : matUpper[0][k] = matLower[0][k] = v2[k+1]-v1[k+1];
4138 0 : matUpper[1][k] = matLower[1][k] = -(u2[k+1]-u1[k+1]);
4139 0 : matUpper[2][k] = decays[0][k+1] - decays[1][k+1]
4140 0 : - decays[0][0] * v1[k+1] + decays[1][0] * u1[k+1];
4141 0 : matLower[2][k] = v1[k+1] - u1[k+1];
4142 : }
4143 0 : double t = -determinant3(matUpper) / determinant3(matLower);
4144 :
4145 : // Find alpha and beta of string crossing.
4146 0 : double s11 = matUpper[0][0]*(t - decays[0][0]);
4147 0 : double s12 = matUpper[1][0]*(t - decays[1][0]);
4148 0 : double s13 = matUpper[2][0] + matLower[2][0]*t;
4149 0 : double s21 = matUpper[0][1]*(t - decays[0][0]);
4150 0 : double s22 = matUpper[1][1]*(t - decays[1][0]);
4151 0 : double s23 = matUpper[2][1] + matLower[2][1]*t;
4152 0 : double den = s11*s22 - s12*s21;
4153 0 : double alpha = (s12*s23 - s22*s13)/den;
4154 0 : double beta = (s21*s13 - s11*s23)/den;
4155 :
4156 : // Check if solution is physical.
4157 0 : if (alpha < 0 || alpha > 1) continue;
4158 0 : if (beta < 0 || beta > 1) continue;
4159 0 : if (t < max(decays[0][0], decays[1][0])) continue;
4160 :
4161 : // Find the crossing points.
4162 0 : Vec4 x1,x2;
4163 0 : x1 = decays[0] + (alpha * v2 + (1 - alpha) * v1) * (t - decays[0][0]);
4164 0 : x2 = decays[1] + (beta * u2 + (1 - beta) * u1) * (t - decays[1][0]);
4165 0 : x1[0] = x2[0] = t;
4166 :
4167 : // Check that both points are the same.
4168 0 : if (dot3(x1-x2,x1-x2) > 1E-4 * (x1.pAbs2() + x2.pAbs2())) continue;
4169 :
4170 : // Find string eigentimes.
4171 0 : double tauPlus = (x1 - decays[0]).mCalc();
4172 0 : double tauMinus = (x2 - decays[1]).mCalc();
4173 :
4174 : // Check if strings already decayed.
4175 0 : if (rndmPtr->flat() > exp( -(pow2(tauPlus) + pow2(tauMinus))
4176 0 : / pow2(tfrag))) continue;
4177 :
4178 : // Optionally only accept if reconnection reduces string length.
4179 0 : if (lowerLambdaOnly) {
4180 0 : double oldLambda = stringLength.getStringLength(event,
4181 0 : dips[0][i].iCol, dips[0][i].iAcol) +
4182 0 : stringLength.getStringLength(event, dips[1][j].iCol,
4183 0 : dips[1][j].iAcol);
4184 0 : double newLambda = stringLength.getStringLength(event,
4185 0 : dips[0][i].iCol, dips[1][j].iAcol) +
4186 0 : stringLength.getStringLength(event, dips[1][j].iCol,
4187 0 : dips[0][i].iAcol);
4188 0 : if (oldLambda < newLambda) continue;
4189 0 : }
4190 :
4191 : // Store dipole pair.
4192 0 : reconnections.insert(make_pair(t,make_pair(i,j)));
4193 0 : }
4194 : }
4195 :
4196 : return reconnections;
4197 0 : }
4198 :
4199 :
4200 : //--------------------------------------------------------------------------
4201 :
4202 : // Calculate the determinant of a 3 * 3 matrix.
4203 :
4204 : double ColourReconnection::determinant3(vector<vector< double> >& vec) {
4205 0 : return vec[0][0]*vec[1][1]*vec[2][2] + vec[0][1]*vec[1][2]*vec[2][0]
4206 0 : + vec[0][2]*vec[1][0]*vec[2][1] - vec[0][0]*vec[2][1]*vec[1][2]
4207 0 : - vec[0][1]*vec[1][0]*vec[2][2] - vec[0][2]*vec[1][1]*vec[2][0];
4208 : }
4209 :
4210 : //==========================================================================
4211 :
4212 : } // end namespace Pythia8
|
