// ColourReconnection.cc is a part of the PYTHIA event generator. // Copyright (C) 2017 Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. // Function definitions (not found in the header) for the // ColourReconnection class. #include "Pythia8/ColourReconnection.h" namespace Pythia8 { //========================================================================== // The BeamDipole class is purely internal to reconnectMPIs. class BeamDipole { public: // Constructor. BeamDipole( int colIn = 0, int iColIn = 0, int iAcolIn = 0) : col(colIn), iCol(iColIn), iAcol(iAcolIn) {} // Members. int col, iCol, iAcol; double p1p2; }; //========================================================================== // The ColourDipole class. //-------------------------------------------------------------------------- // Printing function, intended for debugging. void ColourDipole::list() { cout << setw(10) << this << setw(6) << col << setw(3) << colReconnection << setw(6) << iCol << setw(5) << iAcol << setw(6) << iColLeg << setw(5) << iAcolLeg << setw(6) << isJun << setw(5) << isAntiJun << setw(10) << p1p2 << " colDips: "; for (int i = 0;i < int(colDips.size());++i) cout << setw(10) << colDips[i]; cout << " acolDips: "; for (int i = 0;i < int(acolDips.size());++i) cout << setw(10) << acolDips[i]; cout << setw(3) << isActive << endl; } //========================================================================== // The InfoGluonMove class is purely internal to reconnectMove. class InfoGluonMove{ public: // Constructors. InfoGluonMove(int i1in, int col1in, int acol1in, int iCol1in, int iAcol1in, int col2in, int iCol2in, int iAcol2in, double lambdaRefIn, double dLambdaIn) : i1(i1in), i2(0), col1(col1in), acol1(acol1in), iCol1(iCol1in), iAcol1(iAcol1in), col2(col2in), iCol2(iCol2in), iAcol2(iAcol2in), lambdaRef(lambdaRefIn), dLambda(dLambdaIn) {} InfoGluonMove(int i1in, int i2in, int iCol1in, int iAcol1in, int iCol2in, int iAcol2in, int dLambdaIn) : i1(i1in), i2(i2in), col1(0), acol1(0), iCol1(iCol1in), iAcol1(iAcol1in), col2(0), iCol2(iCol2in), iAcol2(iAcol2in), lambdaRef(0.), dLambda(dLambdaIn) {} // Members. int i1, i2, col1, acol1, iCol1, iAcol1, col2, iCol2, iAcol2; double lambdaRef, dLambda; }; //========================================================================== // The ColourJunction class. //-------------------------------------------------------------------------- // Printing function, intended for debugging. void ColourJunction::list() { cout << setw(6) << kind() << setw(6) << col(0) << setw(6) << col(1) << setw(6) << col(2) << setw(6) << endCol(0) << setw(6) << endCol(1) << setw(6) << endCol(2) << setw(6) << status(0) << setw(6) << status(1) << setw(6) << status(2) << setw(10) << dips[0] << setw(10) << dips[1] << setw(10) << dips[2] << setw(10) << "\n"; cout << " " << setw(10) << dipsOrig[0] << setw(10) << dipsOrig[1] << setw(10) << dipsOrig[2] << endl; } //========================================================================== // The ColourParticle class. //-------------------------------------------------------------------------- // Printing function, intended for debugging. void ColourParticle::listParticle() { const Particle& pt = (*this); // Basic line for a particle, always printed. cout << setw(10) << pt.id() << " " << left << setw(18) << pt.nameWithStatus(18) << right << setw(4) << pt.status() << setw(6) << pt.mother1() << setw(6) << pt.mother2() << setw(6) << pt.daughter1() << setw(6) << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol() << setprecision(3) << setw(11) << pt.px() << setw(11) << pt.py() << setw(11) << pt.pz() << setw(11) << pt.e() << setw(11) << pt.m() << "\n"; } //-------------------------------------------------------------------------- // Printing function, intended for debugging. void ColourParticle::listActiveDips() { cout << "active dips: " << endl; for (int i = 0; i < int(activeDips.size()); ++i) activeDips[i]->list(); } //-------------------------------------------------------------------------- // Printing function, intended for debugging. void ColourParticle::listDips() { cout << "--- Particle ---" << endl; for (int i = 0; i < int(dips.size()); ++i) { cout << "(" <iCol << " (" << dips[i][j]->col << ") "; if (j == int(dips[i].size() - 1)) cout << dips[i][j]->iAcol << " (" << acolEndIncluded[i] << ")" << endl; } } } //========================================================================== // The ColourReconnection class. // Minimum needed gain in lambda for a reconnection (to avoid infinity loops). const double ColourReconnection::MINIMUMGAIN = 1E-10; // Minimum needed gain in lambda for a junction. const double ColourReconnection::MINIMUMGAINJUN = 1E-10; // Conversion of GeV^{-1} to fm for time calculations. const double ColourReconnection::HBAR = 0.197327; // Maximum number of reconnection per trial. // For very large number of outgoing partons, ie. if multiple pp collisions // are stacked on top of each other, this number needs to be raised. const int ColourReconnection::MAXRECONNECTIONS = 1000; //-------------------------------------------------------------------------- // Simple comparison function for sort. bool cmpTrials(TrialReconnection j1, TrialReconnection j2) { return (j1.lambdaDiff < j2.lambdaDiff); } //-------------------------------------------------------------------------- // Initialization. bool ColourReconnection::init( Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn, ParticleData* particleDataPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, PartonSystems* partonSystemsPtrIn) { // Save pointers. infoPtr = infoPtrIn; rndmPtr = rndmPtrIn; particleDataPtr = particleDataPtrIn; beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn; partonSystemsPtr = partonSystemsPtrIn; // Total and squared CM energy at nominal energy. eCM = infoPtr->eCM(); sCM = eCM * eCM; // Choice of reconnection model. reconnectMode = settings.mode("ColourReconnection:mode"); // pT0 scale of MPI; used in the MPI-based reconnection model. pT0Ref = settings.parm("MultipartonInteractions:pT0Ref"); ecmRef = settings.parm("MultipartonInteractions:ecmRef"); ecmPow = settings.parm("MultipartonInteractions:ecmPow"); pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow); // Parameter of the MPI-based reconnection model. reconnectRange = settings.parm("ColourReconnection:range"); pT20Rec = pow2(reconnectRange * pT0); // Parameters of the new reconnection model. m0 = settings.parm("ColourReconnection:m0"); m0sqr = pow2(m0); allowJunctions = settings.flag("ColourReconnection:allowJunctions"); nReconCols = settings.mode("ColourReconnection:nColours"); sameNeighbourCol = settings.flag("ColourReconnection:sameNeighbourColours"); timeDilationMode = settings.mode("ColourReconnection:timeDilationMode"); timeDilationPar = settings.parm("ColourReconnection:timeDilationPar"); timeDilationParGeV = timeDilationPar / HBAR; // Parameters of gluon-move model. m2Lambda = settings.parm("ColourReconnection:m2Lambda"); fracGluon = settings.parm("ColourReconnection:fracGluon"); dLambdaCut = settings.parm("ColourReconnection:dLambdaCut"); flipMode = settings.mode("ColourReconnection:flipMode"); // Parameters of the e+e- CR models. singleReconOnly = settings.flag("ColourReconnection:singleReconnection"); lowerLambdaOnly = settings.flag("ColourReconnection:lowerLambdaOnly"); tfrag = settings.parm("ColourReconnection:fragmentationTime"); blowR = settings.parm("ColourReconnection:blowR"); blowT = settings.parm("ColourReconnection:blowT"); rHadron = settings.parm("ColourReconnection:rHadron"); kI = settings.parm("ColourReconnection:kI"); // Initialize StringLength class. stringLength.init(infoPtr, settings); // Done. return true; } //-------------------------------------------------------------------------- // Do colour reconnection for current event. bool ColourReconnection::next( Event& event, int iFirst) { // MPI-based reconnection model. if (reconnectMode == 0) return reconnectMPIs(event, iFirst); // New reconnection model. else if (reconnectMode == 1) return nextNew(event, iFirst); // Gluon-move model. else if (reconnectMode == 2) return reconnectMove(event, iFirst); // e+e- Type I CR model. else if (reconnectMode == 3 || reconnectMode == 4) return reconnectTypeCommon(event, iFirst); // Undefined. else { infoPtr->errorMsg("Warning in ColourReconnection::next: " "Colour reconnecion mode not found"); return true; } } //-------------------------------------------------------------------------- // Do new colour reconnection for current event. bool ColourReconnection::nextNew( Event& event, int iFirst) { // Clear old records. while (!dipoles.empty()) { delete dipoles.back(); dipoles.pop_back(); } particles.clear(); junctions.clear(); junTrials.clear(); dipTrials.clear(); formationTimes.clear(); // Setup dipoles and make pseudo particles. setupDipoles(event, iFirst); makeAllPseudoParticles(event, iFirst); // Setup all dipole reconnections. // Split dipoles into the 9 different "colours". vector > iDips; iDips.resize(nReconCols); for (int i = 0; i < int(iDips.size()); ++i) iDips[i] = vector(); for (int i = 0; i < int(dipoles.size()); ++i) if (dipoles[i]->isActive) iDips[dipoles[i]->colReconnection].push_back(i); // Loop over each colour individually. for (int i = 0;i < int(iDips.size()); ++i) for (int j = 0; j < int(iDips[i].size()); ++j) for (int k = j + 1; k < int(iDips[i].size()); ++k) singleReconnection(dipoles[iDips[i][j]], dipoles[iDips[i][k]]); // Only do warning once per event. bool alreadyWarned = false; // Start outer loop over reconnections. for (int iOuterLoop = 0; iOuterLoop < 20; ++iOuterLoop) { bool finished = true; // Do inner loop for string reconnections for (int iInnerLoop = 0;dipTrials.size() > 0; ++iInnerLoop) { // Break if too many reonnections are carried out. if (iInnerLoop > MAXRECONNECTIONS) { if (!alreadyWarned) infoPtr->errorMsg("Warning in ColourReconnection::nextNew:" "Too many reconnections, stopping before minimum reached."); alreadyWarned = true; break; } // Store all dipoles connected to the chosen dipole. usedDipoles.clear(); storeUsedDips(dipTrials.back()); // Do the reconnection. doDipoleTrial(dipTrials.back()); // Sort the used dipoles and remove copies of the same. sort(usedDipoles.begin(), usedDipoles.end()); for (int i = 0;i < int(usedDipoles.size() - 1); ++i) if (usedDipoles[i] == usedDipoles[i + 1]) { usedDipoles.erase(usedDipoles.begin() + i); i--; } // Updating the dipole trials. updateDipoleTrials(); } // Loop over list of dipoles to try and form junction structures. if (allowJunctions) { // Split dipoles into three categories. iDips.clear(); iDips.resize(3); for (int i = 0; i < int(iDips.size()); ++i) iDips[i] = vector(); for (int i = 0; i < int(dipoles.size()); ++i) if (dipoles[i]->isActive) iDips[dipoles[i]->colReconnection % 3].push_back(i); // Loop over different "colours" (now only three different groups). for (int i = 0;i < int(iDips.size()); ++i) for (int j = 0; j < int(iDips[i].size()); ++j) for (int k = j + 1; k < int(iDips[i].size()); ++k) singleJunction(dipoles[iDips[i][j]], dipoles[iDips[i][k]]); // Loop over different "colours" (now only three different groups). for (int i = 0;i < int(iDips.size()); ++i) for (int j = 0; j < int(iDips[i].size()); ++j) for (int k = j + 1; k < int(iDips[i].size()); ++k) for (int l = k + 1; l < int(iDips[i].size()); ++l) singleJunction(dipoles[iDips[i][j]], dipoles[iDips[i][k]], dipoles[iDips[i][l]]); // Do inner loop for junction reconnections for (int iInnerLoop = 0;junTrials.size() > 0; ++iInnerLoop) { // Break if too many reonnections are carried out. if (iInnerLoop > MAXRECONNECTIONS) { if (!alreadyWarned) infoPtr->errorMsg("Warning in ColourReconnection::nextNew:" "Too many reconnections, stopping before minimum reached."); alreadyWarned = true; break; } // Find all dipoles connected to the reconnection. usedDipoles.clear(); storeUsedDips(junTrials.back()); // Do the reconnection. doJunctionTrial(event, junTrials.back()); // Sort the used dipoles and remove copies of the same. sort(usedDipoles.begin(), usedDipoles.end()); for (int i = 0;i < int(usedDipoles.size() - 1); ++i) if (usedDipoles[i] == usedDipoles[i + 1]) { usedDipoles.erase(usedDipoles.begin() + i); i--; } // Update lists. updateJunctionTrials(); updateDipoleTrials(); finished = false; } } // If no junctions were made, the overall loop is finished. if (finished) break; } updateEvent(event, iFirst); // Done. return true; } //-------------------------------------------------------------------------- // Setup initial guess on dipoles, here all colours are assumed // to be found in the final state. void ColourReconnection::setupDipoles( Event& event, int iFirst) { // Make vectors needed for storage of chains. vector< vector > chains; vector isJun; vector isAntiJun; vector isGluonLoop; vector inChain(event.size(),false); // Find all quarks and anti-diquarks and follow untill no more colour. for (int i = iFirst; i < event.size(); ++i) { if ( (event[i].isFinal() && !inChain[i] && event[i].isQuark() && event[i].id() > 0) || (event[i].isFinal() && !inChain[i] && event[i].isDiquark() && event[i].id() < 0) ) { int curCol = event[i].col(); inChain[i] = true; vector chain; chain.push_back(i); isAntiJun.push_back(false); isJun.push_back(false); isGluonLoop.push_back(false); for (int iSteps = 0; curCol != 0 && iSteps < 1000; ++iSteps) { // Check for particles with correct anti colour. for (int j = iFirst; j < event.size(); j++) { if (event[j].isFinal() && !inChain[j] && event[j].acol() == curCol) { chain.push_back(j); inChain[j] = true; curCol = event[j].col(); break; } } // Check for junction with correct colour. for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) { for (int j = 0; j < 3; ++j) { if (event.colJunction(iJun,j) == curCol) { isJun[isJun.size() -1] = true; curCol = 0; chain.push_back( -(10 + 10 * iJun + j) ); } } } } chains.push_back(chain); } } // Start from anti-junction and make chains. for (int i = 0; i < event.sizeJunction(); ++i) { // First check if junction belongs to the right diffractive system. int checkCol = event.colJunction(i,0); bool wrongSystem = false; for (int j = 0; j < iFirst; ++j) if (event[j].isFinal() && event[j].acol() == checkCol) wrongSystem = true; if (wrongSystem) continue; // Loop over legs of antijunction. if (event.kindJunction(i) == 2) for (int jCol = 0; jCol < 3; ++jCol) { int curCol = event.colJunction(i,jCol); vector chain; chain.push_back( -(10 + 10 * i + jCol)); isAntiJun.push_back(true); isJun.push_back(false); isGluonLoop.push_back(false); for (int iSteps = 0; curCol != 0 && iSteps < 1000; ++iSteps) { // Check for particles with correct anti colour. for (int j = iFirst; j < event.size(); ++j) if (event[j].isFinal() && !inChain[j] && event[j].acol() == curCol) { chain.push_back(j); inChain[j] = true; curCol = event[j].col(); break; } // Check for junction with correct colour. for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) if (event.kindJunction(iJun) == 1) for (int j = 0; j < 3; ++j) if (event.colJunction(iJun,j) == curCol) { isJun[isJun.size() - 1] = true; curCol = 0; chain.push_back( -(10 + 10 * iJun + j)); } } chains.push_back(chain); } } // Find all gluon loops. for (int i = iFirst; i < event.size(); ++i) if (event[i].isFinal() && !inChain[i] && event[i].col() != 0) { int curCol = event[i].col(); inChain[i] = true; vector chain; chain.push_back(i); isAntiJun.push_back(false); isJun.push_back(false); isGluonLoop.push_back(true); for (int iSteps = 0; curCol != 0 && iSteps < 1000; ++iSteps) { bool foundNext = false; for (int j = iFirst; j < event.size(); ++j) if (event[j].isFinal() && !inChain[j] && event[j].acol() == curCol) { chain.push_back(j); inChain[j] = true; curCol = event[j].col(); foundNext = true; break; } if (!foundNext) break; } chains.push_back(chain); } // Form dipoles from chains. for (int i = 0; i < int(chains.size()); ++i) { if (chains[i].size() == 1) continue; int lastCol = -1; int firstCol = 0; // Start from the first and form the dipoles. // Two consececutive dipoles can not share the same colour. for (int j = 0; j < int(chains[i].size()); ++j) { if (j != int(chains[i].size() - 1)) { // Start by picking new colour. int newCol = int(rndmPtr->flat() * nReconCols); while (newCol == lastCol && !sameNeighbourCol) { newCol = int(rndmPtr->flat() * nReconCols); } // Need to check whether the quark comes from a g->qqbar split. // If that is the case, it can not have the same as qbar. if (j == 0 && !isAntiJun[i] && !isGluonLoop[i]) { int iMother = event[event[ chains[i][j] ].iTopCopy()].mother1(); if ( event[iMother].idAbs() == 21) { vector sisters = event[ chains[i][j] ].sisterList(true); // Need to have only one sister and need to be the anti particle. if (sisters.size() == 1 && event[ chains[i][j] ].id() == - event[ sisters[0] ].id()) { // Try to find dipole with sister. int colSis = -1; for (int k = 0; k < int(dipoles.size()); ++k) if (dipoles[k]->iAcol == sisters[0]) { colSis = dipoles[k]->colReconnection; break; } // If the two colours are the same, pick a new. while (colSis == newCol && !sameNeighbourCol) newCol = int(rndmPtr->flat() * nReconCols); } } } // Check if quark end comes from g->qqbar split. // If so check that the two quarks get different colours. if (j == int(chains[i].size() - 2) && !isJun[i] && !isGluonLoop[i]) { int iMother = event[event[chains[i][j + 1]].iTopCopy()].mother1(); if (event[iMother].idAbs() == 21) { vector sisters = event[ chains[i][j + 1] ].sisterList(true); // Need to have only one sister and need to be the anti particle. if (sisters.size() == 1 && event[ chains[i][j + 1] ].id() == - event[ sisters[0] ].id()) { // Try to find dipole with sister. int colSis = -1; for (int k = 0; k < int(dipoles.size()); ++k) if (dipoles[k]->iCol == sisters[0]) { colSis = dipoles[k]->colReconnection; break; } // If the two colours are the same, pick a new. while ((colSis == newCol || newCol == lastCol) && !sameNeighbourCol) newCol = int(rndmPtr->flat() * nReconCols); } } } // Special case for junction splitting if multiple gluons // between the junctions. if ((chains[i][j] > 0 && event[chains[i][j]].status() == 75) || (chains[i][j + 1] > 0 && event[ chains[i][j + 1] ].status() == 75) ) { // Find sisters. vector sisters; if (chains[i][j] > 0 && event[ chains[i][j] ].status() == 75) sisters = event[ chains[i][j] ].sisterList(); else sisters = event[ chains[i][j + 1] ].sisterList(); if (sisters.size() == 3 ) { // Find colour of sisters. int acolSis1 = -1, acolSis2 = -1, acolSis3 = -1; int colSis1 = -1, colSis2 = -1, colSis3 = -1; for (int k = 0;k < int(dipoles.size()); ++k) { if (dipoles[k]->iAcol == sisters[0]) acolSis1 = dipoles[k]->colReconnection; if (dipoles[k]->iAcol == sisters[1]) acolSis2 = dipoles[k]->colReconnection; if (dipoles[k]->iAcol == sisters[2]) acolSis3 = dipoles[k]->colReconnection; if (dipoles[k]->iCol == sisters[0]) colSis1 = dipoles[k]->colReconnection; if (dipoles[k]->iCol == sisters[1]) colSis2 = dipoles[k]->colReconnection; if (dipoles[k]->iCol == sisters[2]) colSis3 = dipoles[k]->colReconnection; } // If the two colours are the same, pick a new. while ((colSis1 == newCol || colSis2 == newCol || colSis3 == newCol || acolSis1 == newCol || acolSis2 == newCol || acolSis3 == newCol) && !sameNeighbourCol) newCol = int(rndmPtr->flat() * nReconCols); } } // Update stored colours. if (j == 0) firstCol = newCol; lastCol = newCol; // Check if it is antijunction need special dipole. if (j == 0 && isAntiJun[i]) { int col = event.colJunction( - int(chains[i][j]/10) - 1, -chains[i][j] % 10); dipoles.push_back(new ColourDipole(col, chains[i][j], chains[i][j+1], newCol)); dipoles.back()->isAntiJun = true; } // Otherwise just make the dipole. else dipoles.push_back(new ColourDipole(event[ chains[i][j] ].col(), chains[i][j], chains[i][j+1], newCol)); // If the chain in end a junction mark it. if (j == int(chains[i].size() - 2) && isJun[i]) dipoles.back()->isJun = true; // Update the links between dipoles. if (j > 0) { dipoles[dipoles.size() - 1]->leftDip = dipoles[dipoles.size() - 2]; dipoles[dipoles.size() - 2]->rightDip = dipoles[dipoles.size() - 1]; } } // If last particle has anti colour it should be possible to connect it // to the first particle in the chain. (e.g. gluon loop) else if (isGluonLoop[i]) if (event[ chains[i][j] ].col() == event[ chains[i][0] ].acol()) { int newCol = int(rndmPtr->flat() * nReconCols); while ( (newCol == lastCol || newCol == firstCol) && !sameNeighbourCol) { newCol = int(rndmPtr->flat() * nReconCols); } dipoles.push_back(new ColourDipole(event[ chains[i][j] ].col(), chains[i][j], chains[i][0], newCol)); // Update links between dipoles. dipoles[dipoles.size() - 1]->leftDip = dipoles[dipoles.size() - 2]; dipoles[dipoles.size() - 2]->rightDip = dipoles[dipoles.size() - 1]; dipoles[dipoles.size() - chains[i].size()]->leftDip = dipoles[dipoles.size() -1]; dipoles[dipoles.size() - 1]->rightDip = dipoles[dipoles.size() - chains[i].size()]; } } } // Setup junction list. iColJun.clear(); iColJun.resize(event.sizeJunction()); for (int i = 0; i < int(iColJun.size()); ++i) iColJun[i] = vector(3,0); // Loop over event and store indices. for (int i = iFirst; i < event.size(); ++i) if (event[i].isFinal()) for (int j = 0; j < event.sizeJunction(); ++j) for (int jLeg = 0; jLeg < 3; ++jLeg) if (event[i].col() == event.colJunction(j,jLeg) || event[i].acol() == event.colJunction(j,jLeg)) iColJun[j][jLeg] = i; // Loop over junction and store indices. for (int i = 0;i < event.sizeJunction(); ++i) for (int iLeg = 0; iLeg < 3; ++iLeg) for (int j = i + 1;j < event.sizeJunction(); ++j) for (int jLeg = 0; jLeg < 3; ++jLeg) if (event.colJunction(i, iLeg) == event.colJunction(j, jLeg)) { iColJun[i][iLeg] = -(10*j + 10 + jLeg); iColJun[j][jLeg] = -(10*i + 10 + iLeg); } // Setup formation times. setupFormationTimes(event); // Done. } //-------------------------------------------------------------------------- // Calculate the string length of a dipole. double ColourReconnection::calculateStringLength(ColourDipole * dip, vector &dips) { // Check if dipole was already included. for (int i = 0; i < int(dips.size()); ++i) if (dips[i] == dip) return 0.; // Ordinay dipole. if (!dip->isJun && !dip->isAntiJun) { return calculateStringLength(dip->iCol, dip->iAcol); } else { // Start by finding all particles connected to the junction system. vector iParticles; vector usedJuns(junctions.size(),false); int nJuns = 0; if (dip->isJun) { if (!findJunctionParticles( -int(dip->iAcol/10) - 1, iParticles, usedJuns, nJuns, dips)) return 1e9; } else if (!findJunctionParticles(-int(dip->iCol/10) - 1, iParticles, usedJuns, nJuns, dips)) return 1e9; // If it is a single junction. if (int(iParticles.size()) == 3) return calculateJunctionLength(iParticles[0], iParticles[1], iParticles[2]); // If it is a junction pair. else if (int(iParticles.size()) == 4) { return calculateDoubleJunctionLength(iParticles[0], iParticles[1], iParticles[2], iParticles[3]); } // If any other number of junction legs return high number. else return 1e9; } } //-------------------------------------------------------------------------- // Update all colours in the event. void ColourReconnection::updateEvent( Event& event, int iFirst) { // Start by making a new copy of particles. int oldSize = event.size(); for (int i = iFirst; i < oldSize;++i) if (event[i].isFinal()) event.copy(i, 79); // Copy over junctions. event.clearJunctions(); for (int i = 0; i < int(junctions.size()); ++i) { for (int j = 0; j < 3; ++j) { if ( junctions[i].dipsOrig[j] != 0) { junctions[i].col(j, junctions[i].dipsOrig[j]->col); } } event.appendJunction(Junction(junctions[i])); } // Assign colour according to the real dipoles. for (int i = 0; i < int(dipoles.size()); ++i) if (dipoles[i]->isReal) { if (dipoles[i]->iCol >= 0) event[ event[ dipoles[i]->iCol ].daughter1() ].col(dipoles[i]->col); else event.colJunction(-(dipoles[i]->iCol/10 + 1), -dipoles[i]->iCol % 10, dipoles[i]->col); if (dipoles[i]->iAcol >= 0) event[ event[ dipoles[i]->iAcol ].daughter1() ].acol(dipoles[i]->col); else event.colJunction(-(dipoles[i]->iAcol/10 + 1), -dipoles[i]->iAcol % 10, dipoles[i]->col); } } //-------------------------------------------------------------------------- // Find all the particles connected in the junction. // If a single junction, the size of iParticles should be 3. // For multiple junction structures, the size will increase. bool ColourReconnection::findJunctionParticles(int iJun, vector& iParticles, vector &usedJuns, int &nJuns, vector &dips ) { // Mark current junction as used. usedJuns[iJun] = true; nJuns++; // It is not possible to handle junction structures larger than two. if (nJuns > 2) return false; // Find particles connected to the if (junctions[iJun].kind() % 2 == 1) for (int i = 0; i < 3; ++i) iParticles.push_back(junctions[iJun].dips[i]->iCol); else for (int i = 0; i < 3; ++i) iParticles.push_back(junctions[iJun].dips[i]->iAcol); // Add dipoles if not already included. for (int i = 0; i < 3; ++i) { bool addDip = true; for (int j = 0; j < int(dips.size()); ++j) { if (dips[j] == junctions[iJun].dips[i]) { addDip = false; break; } } if (addDip) dips.push_back(junctions[iJun].dips[i]); } // Check whether it connects to any other junctions. for (int i = 0; i < int(iParticles.size()); ++i) if (iParticles[i] < 0) { int iNewJun = - int(iParticles[i] / 10) -1; iParticles.erase(iParticles.begin() + i); i--; if (!usedJuns[iNewJun] && !findJunctionParticles( iNewJun, iParticles, usedJuns, nJuns, dips) ) return false; } // Done. return true; } //-------------------------------------------------------------------------- // Calculate string length for two indices in the particles record. double ColourReconnection::calculateStringLength(int i, int j) { return stringLength.getStringLength(particles[i].p(), particles[j].p()); } //-------------------------------------------------------------------------- // Calculate the length of a single junction given the 3 entries in the event. double ColourReconnection::calculateJunctionLength(int i, int j, int k) { // Need to be separate indices. if ( i == j || i == k || j == k) return 1e9; Vec4 p1 = particles[i].p(); Vec4 p2 = particles[j].p(); Vec4 p3 = particles[k].p(); return stringLength.getJuncLength(p1, p2, p3); } //-------------------------------------------------------------------------- // Calculate the length of a double junction given the 4 particle entries. // The first two are expected to be quarks, the second two to be antiquarks. double ColourReconnection::calculateDoubleJunctionLength( int i, int j, int k, int l) { // Need to be separate indices. if (i == j || i == k || i == l || j == k || j == l || k == l) return 1e9; Vec4 p1 = particles[i].p(); Vec4 p2 = particles[j].p(); Vec4 p3 = particles[k].p(); Vec4 p4 = particles[l].p(); return stringLength.getJuncLength(p1, p2, p3, p4); } //-------------------------------------------------------------------------- // Do a single trial emission. void ColourReconnection::singleReconnection(ColourDipole* dip1, ColourDipole* dip2) { // Do nothing if it is the same dipole. if (dip1 == dip2) return; // No colour reconnection if colours do not match. if (dip1->colReconnection != dip2->colReconnection) return; // If it is not active return if (!dip1->isActive || !dip2->isActive) return; // Not possible to connect a gluon with itself. if (dip1->iCol == dip2->iAcol || dip1->iAcol == dip2->iCol) return; // Check that reconnection is allowed by time dilation. if (!checkTimeDilation(dip1, dip2)) return; // Calculate the difference in lambda. double lambdaDiff = getLambdaDiff(dip1, dip2); // Insert into trial reconnection if anything is gained. if (lambdaDiff > MINIMUMGAIN) { TrialReconnection dipTrial(dip1, dip2, 0, 0, 5, lambdaDiff); dipTrials.insert(lower_bound(dipTrials.begin(), dipTrials.end(), dipTrial, cmpTrials), dipTrial); } } //-------------------------------------------------------------------------- // Simple test swap between two dipoles. void ColourReconnection::swapDipoles(ColourDipole* dip1, ColourDipole* dip2, bool back) { // Swap the anti colour of the dipoles. swap(dip1->iAcol, dip2->iAcol); swap(dip1->isJun, dip2->isJun); swap(dip1->iAcolLeg, dip2->iAcolLeg); // Update the active dipoles. Only change 1 active dipole; // this is to avoid problems when switching back. if (dip1->iAcol != dip2->iAcol) { if (!back) { if (dip1->iAcol >= 0) for (int i = 0; i < int(particles[dip1->iAcol].activeDips.size()); ++i) if (particles[dip1->iAcol].activeDips[i] == dip2) { particles[dip1->iAcol].activeDips[i] = dip1; swap1 = i; break; } if (dip2->iAcol >= 0) for (int i = 0; i < int(particles[dip2->iAcol].activeDips.size()); ++i) if (particles[dip2->iAcol].activeDips[i] == dip1) { particles[dip2->iAcol].activeDips[i] = dip2; swap2 = i; break; } } else { if (dip1->iAcol >= 0) particles[dip1->iAcol].activeDips[swap2] = dip1; if (dip2->iAcol >= 0) particles[dip2->iAcol].activeDips[swap1] = dip2; } } // Update list of junctions (only junctions, antijunctions stay the same). for (int i = 0; i < int(junctions.size()); ++i) if (junctions[i].kind() % 2 == 1) for (int iLeg = 0; iLeg < 3; ++iLeg) { if (junctions[i].dips[iLeg] == dip1) { junctions[i].dips[iLeg] = dip2; continue; } if (junctions[i].dips[iLeg] == dip2) { junctions[i].dips[iLeg] = dip1; continue; } } // Done. } //-------------------------------------------------------------------------- // Do a single trial reconnection to form a junction. void ColourReconnection::singleJunction(ColourDipole* dip1, ColourDipole* dip2) { // Do nothing if it is the same dipole. if (dip1 == dip2) return; int iCol1 = dip1->iCol; int iCol2 = dip2->iCol; int iAcol1 = dip1->iAcol; int iAcol2 = dip2->iAcol; // Not possible to connect a gluon with itself. if (iCol1 == iCol2) return; if (iAcol1 == iAcol2) return; // Check that all dipoles are active. if (!dip1->isActive || !dip2->isActive) return; // Do nothing if one of the dipoles is a junction or antijunction. if (dip1->isJun || dip1->isAntiJun) return; if (dip2->isJun || dip2->isAntiJun) return; // Do nothing if it is a pseudo particle that already contains a // baryon number inside of it. if (int(particles[iCol1].dips.size()) != 1 || int(particles[iAcol1].dips.size()) != 1 || int(particles[iCol2].dips.size()) != 1 || int(particles[iAcol2].dips.size()) != 1) return; // Only accept 2/9 of the colour configurations. if ( (dip1->colReconnection) % 3 != dip2->colReconnection % 3) return; if ( (dip1->colReconnection) == dip2->colReconnection) return; // Check that reconnection is allowed by time dilation. if (!checkTimeDilation(dip1, dip2)) return; // Find the colour of the last junction leg. int junCol = (3 - (dip1->colReconnection / 3) - (dip2->colReconnection / 3) ) * 3 + (dip1->colReconnection % 3); // if other than 9 colours. if (nReconCols != 9) { while (junCol < 0 || junCol % 3 != dip1->colReconnection % 3 || junCol == dip1->colReconnection || junCol == dip2->colReconnection) junCol = int(nReconCols * rndmPtr->flat()); } // Store two new dipoles, these will form the anti-junction. ColourDipole* dip3 = dip1; ColourDipole* dip4 = dip2; double lambdaDiff = getLambdaDiff(dip1, dip2, dip3, dip4, 0); if (lambdaDiff > MINIMUMGAINJUN) { TrialReconnection junTrial(dip1, dip2, dip3, dip4, 0, lambdaDiff); junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(), junTrial, cmpTrials), junTrial); } // Outer loop while (true) { // Reset dip4. dip4 = dip2; // If the colour matches that of the junction. if (dip3->colReconnection == junCol) while (true) { // Check if the new colour matches. if (dip4->colReconnection == dip2->colReconnection) if (checkTimeDilation(dip1, dip2, dip3, dip4)) { // Calculate lambda measure and store new dipole if anything is gained. lambdaDiff = getLambdaDiff(dip1, dip2, dip3, dip4, 1); if (lambdaDiff > MINIMUMGAINJUN) { TrialReconnection junTrial(dip1, dip2, dip3, dip4, 1, lambdaDiff); junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(), junTrial, cmpTrials), junTrial); } } // Find the next neighbour. if (!findAntiNeighbour(dip4)) break; // Check for gluon loop. if (dip4 == dip2 || dip4 == dip1) break; } // Done with inner loop. // Reset dip4. dip4 = dip2; // If the colour matches that of the other dipole. if (dip3->colReconnection == dip1->colReconnection) while (true) { // Check if the new colour matches. if (dip4->colReconnection == junCol) if (checkTimeDilation(dip1, dip2, dip3, dip4)) { // Calculate lambda measure and store new dipole if anything is gained. lambdaDiff = getLambdaDiff(dip1, dip2, dip3, dip4, 2); if (lambdaDiff > MINIMUMGAINJUN) { TrialReconnection junTrial(dip1, dip2, dip3, dip4, 2, lambdaDiff); junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(), junTrial, cmpTrials), junTrial); } } // Find the next neighbour. if (!findAntiNeighbour(dip4)) break; // Check for gluon loop. if (dip4 == dip2 || dip4 == dip1) break; } // Done with inner loop. // Find the next neighbour. if (!findAntiNeighbour(dip3)) break; // Check for gluon loop. if (dip3 == dip1 || dip3 == dip2) break; } // Done. } //-------------------------------------------------------------------------- // Do a single trial reconnection to form a junction. void ColourReconnection::singleJunction(ColourDipole* dip1, ColourDipole* dip2, ColourDipole* dip3) { // Do nothing if one of the dipoles is a junction or antijunction. if (dip1->isJun || dip1->isAntiJun) return; if (dip2->isJun || dip2->isAntiJun) return; if (dip3->isJun || dip3->isAntiJun) return; // Check that all dipoles are active. if (!dip1->isActive || !dip2->isActive || !dip3->isActive) return; // Only allow 0-3-6, 1-4-7 or 2-5-8. if ( dip1->colReconnection % 3 != dip2->colReconnection % 3 || dip1->colReconnection % 3 != dip3->colReconnection % 3) return; if ( !(dip1->colReconnection != dip2->colReconnection && dip1->colReconnection != dip3->colReconnection && dip2->colReconnection != dip3->colReconnection) ) return; if (int(particles[dip1->iCol].dips.size()) != 1 || int(particles[dip1->iAcol].dips.size()) != 1 || int(particles[dip2->iCol].dips.size()) != 1 || int(particles[dip2->iAcol].dips.size()) != 1 || int(particles[dip3->iCol].dips.size()) != 1 || int(particles[dip3->iAcol].dips.size()) != 1 ) return; // Check that reconnection is allowed by time dilation. if (!checkTimeDilation(dip1, dip2, dip3)) return; double lambdaDiff = getLambdaDiff(dip1, dip2, dip3, 0, 3); if (lambdaDiff > MINIMUMGAINJUN) { TrialReconnection junTrial(dip1, dip2, dip3, 0, 3, lambdaDiff); junTrials.insert(lower_bound(junTrials.begin(), junTrials.end(), junTrial, cmpTrials), junTrial); } // Done. return; } // ------------------------------------------------------------------ // Form pseuparticle of a given dipole (or junction system). void ColourReconnection::makePseudoParticle(ColourDipole* dip , int status, bool setupDone) { // If it is a normal dipole that needs to be combined. if (!dip->isJun && !dip->isAntiJun) { // Start by storing variables for easier use. int iCol = dip->iCol; int iAcol = dip->iAcol; int iColLeg = dip->iColLeg; int iAcolLeg = dip->iAcolLeg; // Make new pseudo particle. int iNew = particles.size(); particles.push_back(particles[iCol]); particles[iNew].acol(particles[iCol].acol()); particles[iNew].col(particles[iAcol].col()); particles[iNew].mother1(iCol); particles[iNew].mother2(iAcol); particles[iNew].id(99); particles[iNew].status(status); particles[iNew].isJun = false; particles[iAcol].statusNeg(); particles[iAcol].daughter1(iNew); particles[iCol].statusNeg(); particles[iCol].daughter1(iNew); if (iCol != iAcol) particles[iNew].p(particles[iCol].p() + particles[iAcol].p()); else particles[iNew].p(particles[iCol].p()); // Add all the dipoles from the old pseudo particle. // First from particle 1. particles[iNew].dips = particles[dip->iCol].dips; particles[iNew].colEndIncluded = particles[dip->iCol].colEndIncluded; particles[iNew].acolEndIncluded = particles[dip->iCol].acolEndIncluded; // Then particle 2. if (iCol != iAcol) { for (int i = 0; i < int(particles[dip->iAcol].dips.size()); ++i) { if (i != dip->iAcolLeg) { // If it is not the same leg, add as separate vector. particles[iNew].dips.push_back(particles[dip->iAcol].dips[i]); particles[iNew].colEndIncluded.push_back( particles[dip->iAcol].colEndIncluded[i]); particles[iNew].acolEndIncluded.push_back( particles[dip->iAcol].acolEndIncluded[i]); } // If it is the same leg, at at the end of the vector. else { particles[iNew].acolEndIncluded[iColLeg] = particles[iAcol].acolEndIncluded[i]; particles[iNew].dips[iColLeg].pop_back(); particles[iNew].dips[iColLeg].insert( particles[iNew].dips[iColLeg].end(), particles[iAcol].dips[i].begin(), particles[iAcol].dips[i].end() ); } } } if (iCol != iAcol) { // Update the dipole legs to the new particle. for (int i = 0; i < int(particles[iAcol].activeDips.size()); ++i) { if ( particles[iAcol].activeDips[i]->iAcol == iAcol) { if (particles[iAcol].activeDips[i]->iAcolLeg < iAcolLeg) particles[iAcol].activeDips[i]->iAcolLeg += particles[iCol].dips.size(); else if (particles[iAcol].activeDips[i]->iAcolLeg == iAcolLeg) particles[iAcol].activeDips[i]->iAcolLeg = iColLeg; else if (particles[iAcol].activeDips[i]->iAcolLeg > iAcolLeg) particles[iAcol].activeDips[i]->iAcolLeg += particles[iCol].dips.size() - 1; } if (particles[iAcol].activeDips[i]->iCol == iAcol) { if (particles[iAcol].activeDips[i]->iColLeg < iAcolLeg) particles[iAcol].activeDips[i]->iColLeg += particles[iCol].dips.size(); else if (particles[iAcol].activeDips[i]->iColLeg == iAcolLeg) particles[iAcol].activeDips[i]->iColLeg = iColLeg; else if (particles[iAcol].activeDips[i]->iColLeg > iAcolLeg) particles[iAcol].activeDips[i]->iColLeg += particles[iCol].dips.size() - 1; } } } // Update list of active dipoles. particles[iNew].activeDips.clear(); particles[iNew].activeDips.insert(particles[iNew].activeDips.end(), particles[iCol].activeDips.begin(), particles[iCol].activeDips.end()); if (iCol != iAcol) particles[iNew].activeDips.insert(particles[iNew].activeDips.end(), particles[iAcol].activeDips.begin(), particles[iAcol].activeDips.end()); // Remove the now inactive dipole. for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) if (particles[iNew].activeDips[i] == dip) { particles[iNew].activeDips.erase( particles[iNew].activeDips.begin() + i); i--; } // Update the indices in the active dipoles. for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) { if (particles[iNew].activeDips[i]->iCol == iAcol) particles[iNew].activeDips[i]->iCol = iNew; if (particles[iNew].activeDips[i]->iCol == iCol) particles[iNew].activeDips[i]->iCol = iNew; if (particles[iNew].activeDips[i]->iAcol == iAcol) particles[iNew].activeDips[i]->iAcol = iNew; if (particles[iNew].activeDips[i]->iAcol == iCol) particles[iNew].activeDips[i]->iAcol = iNew; particles[iNew].activeDips[i]->p1p2 = mDip(particles[iNew].activeDips[i]); } // If it is a combination of the same particle, // check if any double active dipoles if (iCol == iAcol) for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) for (int j = i + 1; j < int(particles[iNew].activeDips.size()); ++j) if (particles[iNew].activeDips[i] == particles[iNew].activeDips[j]) { particles[iNew].activeDips.erase(particles[iNew].activeDips.begin() + j); j--; } // Add dips changed to used dips. for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) { if (particles[iNew].activeDips[i]->iCol >= 0) usedDipoles.push_back(particles[iNew].activeDips[i]); else for (int j = 0;j < 3; ++j) usedDipoles.push_back(junctions[-(particles[iNew]. activeDips[i]->iCol / 10 + 1)].getColDip(j)); if (particles[iNew].activeDips[i]->iAcol >= 0) usedDipoles.push_back(particles[iNew].activeDips[i]); else for (int j = 0;j < 3; ++j) usedDipoles.push_back(junctions[-(particles[iNew]. activeDips[i]->iAcol / 10 + 1)].getColDip(j)); } // mark the internal dipole as not active. dip->isActive = false; // Done. return; } // If both ends are connected to a junction something went wrong! else if (dip->isJun && dip->isAntiJun) { return; } else { // Find junction index and first leg to combine. int iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2; getJunctionIndices(dip, iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2); ColourDipole* dip2 = junctions[iJun].dips[junLeg1]; ColourDipole* dip3 = junctions[iJun].dips[junLeg2]; // Add new particle. int iNew = particles.size(); particles.push_back(ColourParticle( Particle( 99, status, i0, i1, 0, 0, 0, 0, particles[i0].p() + particles[i1].p() ) ) ); particles[iNew].isJun = true; particles[iNew].junKind = junctions[iJun].kind(); if (i0 == i1) particles[iNew].p(particles[i0].p()); // Update old particles. particles[i0].statusNeg(); particles[i0].daughter1(iNew); particles[i1].statusNeg(); particles[i1].daughter1(iNew); // Update list of internal dipoles. particles[iNew].dips.clear(); particles[iNew].dips.insert(particles[iNew].dips.end(), particles[i0].dips.begin(),particles[i0].dips.end()); if (i0 != i1) particles[iNew].dips.insert(particles[iNew].dips.end(), particles[i1].dips.begin(),particles[i1].dips.end()); // Update list of whether colour ending is included. particles[iNew].colEndIncluded.clear(); particles[iNew].colEndIncluded.insert( particles[iNew].colEndIncluded.end(), particles[i0].colEndIncluded.begin(), particles[i0].colEndIncluded.end() ); if (i0 != i1) particles[iNew].colEndIncluded.insert( particles[iNew].colEndIncluded.end(), particles[i1].colEndIncluded.begin(), particles[i1].colEndIncluded.end() ); // Update list of whether anti colour ending is included. particles[iNew].acolEndIncluded.clear(); particles[iNew].acolEndIncluded.insert( particles[iNew].acolEndIncluded.end(), particles[i0].acolEndIncluded.begin(), particles[i0].acolEndIncluded.end() ); if (i0 != i1) particles[iNew].acolEndIncluded.insert( particles[iNew].acolEndIncluded.end(), particles[i1].acolEndIncluded.begin(), particles[i1].acolEndIncluded.end() ); // Third particle just need to add one to list of dipoles. if (dip->isJun && i2 >= 0 && i2 != i0 && i2 != i1) { particles[iNew].dips.push_back(particles[i2].dips[dip3->iColLeg]); particles[iNew].dips.back().erase(particles[iNew].dips.back().begin(), particles[iNew].dips.back().end() - 1); } if (dip->isAntiJun && i2 >= 0 && i2 != i0 && i2 != i1) { particles[iNew].dips.push_back(particles[i2].dips[dip3->iAcolLeg]); particles[iNew].dips.back().erase( particles[iNew].dips.back().begin() + 1, particles[iNew].dips.back().end() ); } // Add endings for the third particle. if (i2 != i0 && i2 != i1) { particles[iNew].acolEndIncluded.push_back(false); particles[iNew].colEndIncluded.push_back(false); } // Special case if it is J-J connection. if (i2 < 0) { particles[iNew].dips.push_back(vector()); // Find the real dipole to add to dipole list. for (int i = 0; i < int(dipoles.size()); ++i) if (dipoles[i]->isReal && dipoles[i]->iCol == dip3->iCol && dipoles[i]->iAcol == dip3->iAcol) particles[iNew].dips.back().push_back(dipoles[i]); // Change ending. particles[iNew].acolEndIncluded.back() = true; particles[iNew].colEndIncluded.back() = true; } // The endings need to reflect the new junction structure. if (dip->isJun) for (int i = 0; i < int(particles[iNew].acolEndIncluded.size()); ++i) particles[iNew].acolEndIncluded[i] = true; else for (int i = 0; i < int(particles[iNew].colEndIncluded.size()); ++i) particles[iNew].colEndIncluded[i] = true; // Update active dipoles, first junction case. // Set the now internal dipoles as inactive. dip->isActive = false; dip2->isActive = false; dip3->isActive = true; // Update the dipole legs to the new particle. // Only need to do it for the iAcol particle, // since nothing changes for the iCol particle. if (i0 != i1) for (int i = 0; i < int(particles[i1].activeDips.size()); ++i) { if (particles[i1].activeDips[i]->iAcol == i1) particles[i1].activeDips[i]->iAcolLeg += particles[i0].dips.size(); if (particles[i1].activeDips[i]->iCol == i1) particles[i1].activeDips[i]->iColLeg += particles[i0].dips.size(); } // Update list of active dipoles. particles[iNew].activeDips.clear(); particles[iNew].activeDips.insert(particles[iNew].activeDips.end(), particles[i0].activeDips.begin(), particles[i0].activeDips.end()); if (i0 != i1) particles[iNew].activeDips.insert(particles[iNew].activeDips.end(), particles[i1].activeDips.begin(), particles[i1].activeDips.end()); if (i2 != i0 && i2 != i1) particles[iNew].activeDips.push_back(dip3); // Remove the now inactive dipoles. for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) { if (particles[iNew].activeDips[i] == dip) { particles[iNew].activeDips.erase( particles[iNew].activeDips.begin() + i); i--; continue; } if (particles[iNew].activeDips[i] == dip2) { particles[iNew].activeDips.erase( particles[iNew].activeDips.begin() + i); i--; continue; } } // Update the indices in the active dipoles. for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) { if (particles[iNew].activeDips[i]->iCol == i1) particles[iNew].activeDips[i]->iCol = iNew; if (particles[iNew].activeDips[i]->iCol == i0) particles[iNew].activeDips[i]->iCol = iNew; if (particles[iNew].activeDips[i]->iAcol == i1) particles[iNew].activeDips[i]->iAcol = iNew; if (particles[iNew].activeDips[i]->iAcol == i0) particles[iNew].activeDips[i]->iAcol = iNew; particles[iNew].activeDips[i]->p1p2 = mDip(particles[iNew].activeDips[i]); } // The third dip is no longer connected to a junction. if (dip->isJun) { dip3->isJun = false; dip3->iAcol = iNew; if (i2 != i0 && i2 != i1) dip3->iAcolLeg = particles[iNew].dips.size() - 1; } else { dip3->isAntiJun = false; dip3->iCol = iNew; if (i2 != i0 && i2 != i1) dip3->iColLeg = particles[iNew].dips.size() - 1; } // Add dips changed to used dips. for (int i = 0; i < int(particles[iNew].activeDips.size()); ++i) { bool added = false; for (int j = 0;j < int(usedDipoles.size()); ++j) if (particles[iNew].activeDips[i] == usedDipoles[j]) { added = true; break; } if (!added) usedDipoles.push_back(particles[iNew].activeDips[i]); } usedDipoles.push_back(dip); usedDipoles.push_back(dip2); usedDipoles.push_back(dip3); // Possible for the new dip to have a low m0. if (setupDone && mDip(dip3) < m0) makePseudoParticle(dip3, status, true); } // Done. } // ------------------------------------------------------------------ // Help function to sort dipoles in right order. bool sortFunc(ColourDipole* a, ColourDipole* b) { return (a->p1p2 < b->p1p2); } // ------------------------------------------------------------------ // Form all pseudoparticles below m0. void ColourReconnection::makeAllPseudoParticles( Event & event, int iFirst) { // Make junctions. for (int i = 0; i < event.sizeJunction(); ++i) junctions.push_back(event.getJunction(i)); // Make new copy of all the dipoles. int oldSize = int(dipoles.size()); for (int i = 0; i < oldSize; ++i) { dipoles.push_back(new ColourDipole(*dipoles[i])); dipoles[i + oldSize]->iColLeg = 0; dipoles[i + oldSize]->iAcolLeg = 0; dipoles[i]->iColLeg = 0; dipoles[i]->iAcolLeg = 0; dipoles[i]->isActive = false; dipoles[i]->isReal = true; dipoles[i + oldSize]->isReal = false; // Store original dipoles connected to junctions. if (dipoles[i]->iCol < 0) { junctions[-(dipoles[i]->iCol / 10 + 1)].dipsOrig[(-dipoles[i]->iCol) % 10] = dipoles[i]; } if (dipoles[i]->iAcol < 0) { junctions[-(dipoles[i]->iAcol / 10 + 1)].dipsOrig[-(dipoles[i]->iAcol % 10)] = dipoles[i]; } } // Set up the coldDips and acolDips. for (int i = 0; i < oldSize; ++i) { if (dipoles[i]->leftDip != 0) for (int j = 0; j < oldSize; ++j) if (dipoles[i]->leftDip == dipoles[j]) { dipoles[i + oldSize]->colDips.push_back(dipoles[j + oldSize]); break; } if (dipoles[i]->rightDip != 0) for (int j = 0; j < oldSize; ++j) if (dipoles[i]->rightDip == dipoles[j]) { dipoles[i + oldSize]->acolDips.push_back(dipoles[j + oldSize]); break; } } // Start by copying event record to make pseudoparticles. // The pseudoparticles also need to gain for (int i = iFirst; i < event.size(); ++i) if (event[i].isFinal()) { particles.push_back(ColourParticle(event[i])); particles.back().dips.resize(1,vector()); // Set up dipoles. for (int j = 0; j < int(dipoles.size()); ++j) { if (dipoles[j]->iCol == i) { if (dipoles[j]->isActive) { dipoles[j]->iCol = particles.size() - 1; particles.back().activeDips.push_back(dipoles[j]); } else particles.back().dips[0].push_back(dipoles[j]); } if (dipoles[j]->iAcol == i) { if (dipoles[j]->isActive) { dipoles[j]->iAcol = particles.size() - 1; particles.back().activeDips.push_back(dipoles[j]); } else particles.back().dips[0].insert(particles.back().dips[0].begin(), dipoles[j]); } } // Tell whether dipoles are connected to other dipoles. if (event[i].isQuark() && event[i].id() > 0) particles.back().colEndIncluded.push_back(true); else particles.back().colEndIncluded.push_back(false); if (event[i].isQuark() && event[i].id() < 0) particles.back().acolEndIncluded.push_back(true); else particles.back().acolEndIncluded.push_back(false); } // Inserting a copy of the event record, but now with full // pseudo particle setup. // This is mainly to avoid having to distinguish between combining // original particles and pseudoparticles. // Set right dipole connections in junctions. for (int i = 0; i < int(dipoles.size()); ++i) { if (dipoles[i]->iCol < 0) { int j = (- dipoles[i]->iCol / 10) - 1; int jLeg = - dipoles[i]->iCol % 10; junctions[j].setColDip(jLeg, dipoles[i]); } if (dipoles[i]->iAcol < 0) { int j = (- dipoles[i]->iAcol / 10) - 1; int jLeg = - dipoles[i]->iAcol % 10; junctions[j].setColDip(jLeg, dipoles[i]); } } // Make sure all dipoles masses are set correctly. for (int i = 0; i < int(dipoles.size()); ++i) { if (dipoles[i]->isActive) dipoles[i]->p1p2 = mDip(dipoles[i]); else dipoles[i]->p1p2 = 1e9; } // Keep making pseudo particles until they are above the threshold. while (true) { sort(dipoles.begin(), dipoles.end(), sortFunc); bool finished = true; for (int i = 0; i < int(dipoles.size()); ++i) { if (!dipoles[i]->isActive) continue; if (dipoles[i]->p1p2 < m0) { makePseudoParticle( dipoles[i], 110); finished = false; break; } else break; } if (finished) break; } // Sort the dipoles. sort(dipoles.begin(), dipoles.end(), sortFunc); // Done. return; } // ------------------------------------------------------------------ // Print statements if something is wrong in dipole setup. // Does not have a return statement -- DEBUG PURPOSE ONLY --. void ColourReconnection::checkRealDipoles(Event& event, int iFirst) { vector dipConnections(event.size(),0); for (int i = 0;i < int(dipoles.size()); ++i) if (dipoles[i]->isReal) { if (dipoles[i]->iCol >= 0) dipConnections[dipoles[i]->iCol]++; if (dipoles[i]->iAcol >= 0) dipConnections[dipoles[i]->iAcol]++; } bool working = true; for (int i = iFirst ;i < event.size(); ++i) { if (event[i].isFinal()) { if (event[i].isQuark() && dipConnections[i] != 1) { cout << "quark " << i << " is wrong!!" << endl; working = false; } else if (event[i].idAbs() == 21 && dipConnections[i] != 2) { cout << "gluon " << i << " is wrong!!" << endl; working = false; } } } if (!working) { infoPtr->errorMsg("Error in ColourReconnection::checkRealDipoles:" "Real dipoles not setup properply"); } } // ------------------------------------------------------------------ // Print statements if something is wrong in dipole setup. // Does not have a return statement -- DEBUG PURPOSE ONLY --. void ColourReconnection::checkDipoles() { for (int i = 0;i < int(dipoles.size()); ++i) { if (dipoles[i] == 0) { cout << "dipole empty" << endl;} if (dipoles[i]->isActive) { if (dipoles[i]->iCol >= 0) { bool foundMyself = false; for (int j = 0; j < int(particles[ dipoles[i]->iCol ]. activeDips.size()); ++j) { if (!particles[dipoles[i]->iCol].activeDips[j]->isActive) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Found inactive dipole, where only active was expected"); } if (particles[dipoles[i]->iCol].activeDips[j] == dipoles[i]) foundMyself = true; } if (!foundMyself) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Linking between active dipoles and particles is wrong"); } if (dipoles[i]->iColLeg >= int(particles[dipoles[i]->iCol].dips.size())) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Original dipoles not stored correct"); } // Check that linking to old dipoles work. if (dipoles[i]->col != particles[dipoles[i]->iCol].dips[dipoles[i]->iColLeg].back()->col) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Original dipoles do not match in"); } } if (dipoles[i]->iAcol >= 0) { bool foundMyself = false; for (int j = 0;j < int(particles[ dipoles[i]->iAcol ]. activeDips.size()); ++j) { if (!particles[dipoles[i]->iAcol].activeDips[j]->isActive) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Found inactive dipole, where only active was expected"); } if (particles[dipoles[i]->iAcol].activeDips[j] == dipoles[i]) foundMyself = true; } if (!foundMyself) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Linking between active dipoles and particles is wrong"); } if (dipoles[i]->iAcolLeg >= int(particles[dipoles[i]->iAcol]. dips.size() )) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Original dipoles not stored correct"); } // Check that linking to old dipoles work if (dipoles[i]->col != particles[dipoles[i]->iAcol]. dips[dipoles[i]->iAcolLeg].front()->col) { infoPtr->errorMsg("Error in ColourReconnection::checkDipoles:" "Original dipoles do not match in"); } } } } } // ------------------------------------------------------------------ // Print all the chains. void ColourReconnection::listAllChains() { cout << " ----- PRINTING CHAINS ----- " << dipoles.size() << endl; for (int i = 0; i < int(dipoles.size()); ++i) dipoles[i]->printed = false; for (int i = 0;i < int(dipoles.size()); ++i) if (!dipoles[i]->printed) listChain(dipoles[i]); cout << " ----- PRINTED CHAINS ----- " << endl; } // ------------------------------------------------------------------ // Print the chain containing the dipole. void ColourReconnection::listChain(ColourDipole *dip) { // Make sure not an empty pointer. if (dip == 0) return; // If chain is not active, just print it. if (!dip->isActive) { return; } ColourDipole * colDip = dip; // Try to reach one end of the chain. while (particles[colDip->iCol].dips.size() == 1 && findColNeighbour(colDip)) if (dip == colDip) break; ColourDipole * endDip = colDip; do { cout << colDip->iCol << " (" << colDip->p1p2 << ", " << colDip->col << ") (" << colDip->isActive << ") "; colDip->printed = true; } // Start the printing. while (particles[colDip->iAcol].dips.size() == 1 && findAntiNeighbour(colDip) && colDip != endDip); // Print the last part. cout << colDip->iAcol<< endl; // Done. } // ------------------------------------------------------------------ // Return relevant indices for the junction. bool ColourReconnection::getJunctionIndices(ColourDipole * dip, int &iJun, int &i0, int &i1, int &i2, int &junLeg0, int &junLeg1, int &junLeg2) { // Find junction index and first leg to combine. int indxJun = dip->iCol; if (dip->iAcol < 0) indxJun = dip->iAcol; iJun = (- indxJun / 10) - 1; junLeg0 = -(indxJun % 10); junLeg1 = 1; junLeg2 = 2; if (junLeg0 == 1) junLeg1 = 0; else if (junLeg0 == 2) junLeg2 = 0; if (dip->iCol < 0) { i0 = dip->iAcol; i1 = junctions[iJun].dips[junLeg1]->iAcol; i2 = junctions[iJun].dips[junLeg2]->iAcol; } else { i0 = dip->iCol; i1 = junctions[iJun].dips[junLeg1]->iCol; i2 = junctions[iJun].dips[junLeg2]->iCol; } // It is not possible to form a pseudoparticle if only a single particle is // connected to the junction. if (i1 < 0 && i2 < 0) return false; // Check which two particle should form the pseudoparticle. double m1 = 1e9, m2 = 1e9; if (i1 >= 0) m1 = m(particles[i0].p(),particles[i1].p()); if (i2 >= 0) m2 = m(particles[i0].p(),particles[i2].p()); if (m1 > m2) { swap(i1,i2); swap(junLeg1,junLeg2); } // Force switch if i0 == i2 if (i0 == i2) { swap(i1,i2); swap(junLeg1,junLeg2); } return true; } // ------------------------------------------------------------------ // Check whether up to four dipoles are 'causally' connected. bool ColourReconnection::checkTimeDilation(ColourDipole * dip1, ColourDipole * dip2, ColourDipole * dip3, ColourDipole * dip4) { if (timeDilationMode == 0) return true; // 2 dipole case. if (dip3 == 0) { Vec4 p1 = getDipoleMomentum(dip1); Vec4 p2 = getDipoleMomentum(dip2); double t1 = formationTimes[dip1->col]; double t2 = formationTimes[dip2->col]; if (dip1 == dip2) return true; else return checkTimeDilation(p1, p2, t1, t2); // 3 dipole case. } else if (dip4 == 0) { Vec4 p1 = getDipoleMomentum(dip1); Vec4 p2 = getDipoleMomentum(dip2); Vec4 p3 = getDipoleMomentum(dip3); double t1 = formationTimes[dip1->col]; double t2 = formationTimes[dip2->col]; double t3 = formationTimes[dip3->col]; // Modes that require all dipoles to be causally connected. if (timeDilationMode == 1 || timeDilationMode == 2 || timeDilationMode == 4) { if (dip1 != dip2 && !checkTimeDilation(p1, p2, t1, t2)) return false; if (dip1 != dip3 && !checkTimeDilation(p1, p3, t1, t3)) return false; if (dip2 != dip3 && !checkTimeDilation(p2, p3, t2, t3)) return false; return true; // Modes that require a single pair of dipoles to be causally connected. } else { if (dip1 != dip2 && checkTimeDilation(p1, p2, t1, t2)) return true; if (dip1 != dip3 && checkTimeDilation(p1, p3, t1, t3)) return true; if (dip2 != dip3 && checkTimeDilation(p2, p3, t2, t3)) return true; return false; } // 4 dipole case. } else { Vec4 p1 = getDipoleMomentum(dip1); Vec4 p2 = getDipoleMomentum(dip2); Vec4 p3 = getDipoleMomentum(dip3); Vec4 p4 = getDipoleMomentum(dip4); double t1 = formationTimes[dip1->col]; double t2 = formationTimes[dip2->col]; double t3 = formationTimes[dip3->col]; double t4 = formationTimes[dip4->col]; // Modes that require all dipoles to be causally connected. if (timeDilationMode == 1 || timeDilationMode == 2 || timeDilationMode == 4) { if (dip1 != dip2 && !checkTimeDilation(p1, p2, t1, t2)) return false; if (dip1 != dip3 && !checkTimeDilation(p1, p3, t1, t3)) return false; if (dip1 != dip4 && !checkTimeDilation(p1, p4, t1, t4)) return false; if (dip2 != dip3 && !checkTimeDilation(p2, p3, t2, t3)) return false; if (dip2 != dip4 && !checkTimeDilation(p2, p4, t2, t4)) return false; if (dip3 != dip4 && !checkTimeDilation(p3, p4, t3, t4)) return false; return true; // Modes that require a single pair of dipoles to be causally connected. } else { if (dip1 != dip2 && checkTimeDilation(p1, p2, t1, t2)) return true; if (dip1 != dip3 && checkTimeDilation(p1, p3, t1, t3)) return true; if (dip1 != dip4 && checkTimeDilation(p1, p4, t1, t4)) return true; if (dip2 != dip3 && checkTimeDilation(p2, p3, t2, t3)) return true; if (dip2 != dip4 && checkTimeDilation(p2, p4, t2, t4)) return true; if (dip3 != dip4 && checkTimeDilation(p3, p4, t3, t4)) return true; return false; } } // Done. } // ------------------------------------------------------------------ // Find the momentum of the dipole. Vec4 ColourReconnection::getDipoleMomentum(ColourDipole * dip) { vector iPar, usedJuncs; if (!dip->isJun) iPar.push_back(dip->iAcol); else addJunctionIndices(dip->iAcol, iPar, usedJuncs); if (!dip->isAntiJun) iPar.push_back(dip->iCol); else addJunctionIndices(dip->iCol, iPar, usedJuncs); // Remove any duplicates. sort(iPar.begin(),iPar.end()); for (int i = 0;i < int(iPar.size()) - 1; ++i) if (iPar[i] == iPar[i+1]) { iPar.erase(iPar.begin()+i); i--; } if (iPar.size() == 0) { infoPtr->errorMsg("Error in ColourReconnection::getDipoleMomentum: " "No particles connected to junction."); return Vec4(0,0,0,0); } Vec4 p = particles[iPar[0]].p(); for (int i = 1;i < int(iPar.size()); ++i) p += particles[iPar[i]].p(); return p; } // ------------------------------------------------------------------ // Check whether two four momenta are 'causally' connected. bool ColourReconnection::checkTimeDilation(Vec4 p1, Vec4 p2, double t1, double t2) { // No time dilation check. if (timeDilationMode == 0) return true; // Check if gamma is below parameter. else if (timeDilationMode == 1) { p2.bstback(p1); if (p2.e() / p2.mCalc() > timeDilationPar) return false; else return true; // Check if gamma * mDip is below parameter for both dipoles. } else if (timeDilationMode == 2) { bool part1, part2; p2.bstback(p1); if (p2.e() / p2.mCalc() > timeDilationParGeV * p2.mCalc()) part1 = false; else part1 = true; p2.bst(p1); p1.bstback(p2); if (p1.e() / p1.mCalc() > timeDilationParGeV * p1.mCalc()) part2 = false; else part2 = true; if (part1 && part2) return true; else return false; // Check if gamma * mDip is below parameter for a single dipole. } else if (timeDilationMode == 3) { bool part1, part2; p2.bstback(p1); if (p2.e() / p2.mCalc() > timeDilationParGeV * p2.mCalc()) part1 = false; else part1 = true; p2.bst(p1); p1.bstback(p2); if (p1.e() / p1.mCalc() > timeDilationParGeV * p1.mCalc()) part2 = false; else part2 = true; if (part1 || part2) return true; else return false; // Check if gamma * mDip' is below parameter for both dipoles. } else if (timeDilationMode == 4) { p2.bstback(p1); if (p2.e() / p2.mCalc() < timeDilationParGeV * min(t1,t2)) return true; else return false; // Check if gamma * mDip' is below parameter for a single dipole. } else if (timeDilationMode == 5) { p2.bstback(p1); if (p2.e() / p2.mCalc() < timeDilationParGeV * max(t1,t2)) return true; else return false; // If mode is set wrong, should never happen. } else return true; } // ------------------------------------------------------------------ // Store the formation times. void ColourReconnection::setupFormationTimes(Event & event) { for (int i = 0;i < event.size(); ++i) { // First check the colour. if (event[i].col() != 0 && formationTimes.count(event[i].col()) == 0) { int col = event[i].col(); // Find first time the colour appears as an anticolour. bool foundCol = false; int iAcol = 0; for (int j = i;j < event.size(); ++j) { if (event[j].acol() == col) { foundCol = true; iAcol = j; break; } } // If it was found add it to the list. if (foundCol) { formationTimes[col] = max( m0, (event[i].p() + event[iAcol].p()).mCalc() ); // Otherwise it must be stored in a junction. } else { formationTimes[col] = max(m0, getJunctionMass(event, col)); } } // Next check the anti colour. if (event[i].acol() != 0 && formationTimes.count(event[i].acol()) == 0) { int acol = event[i].acol(); // Find first time the colour appears as an anticolour. bool foundCol = false; int iCol = 0; for (int j = i;j < event.size(); ++j) { if (event[j].col() == acol) { foundCol = true; iCol = j; break; } } // If it was found add it to the list. if (foundCol) { formationTimes[acol] = max(m0, (event[i].p() + event[iCol].p()) .mCalc()); // Otherwise it must be stored in a junction. } else { formationTimes[acol] = max(m0, getJunctionMass(event, acol)); } } } // Finally check if junction colours are stored. for (int i = 0; i < event.sizeJunction(); ++i) for (int j = 0; j < 3; ++j) if (formationTimes.count(event.colJunction(i,j)) == 0) formationTimes[event.colJunction(i,j)] = max(m0, getJunctionMass(event, event.colJunction(i,j))); // Done. } // ------------------------------------------------------------------ // Find the invariant mass of all the partons connected to a junction system. double ColourReconnection::getJunctionMass(Event & event, int col) { // Find the partons connected to the junction system. vector iPar, usedJuncs; addJunctionIndices(event, col, iPar, usedJuncs); // Check for doubles. sort(iPar.begin(), iPar.end()); for (int i = 0;i < int(iPar.size() -1); ++i) { if (iPar[i] == iPar[i + 1]) { iPar.erase(iPar.begin() + i); i--; } } // If no partons are connected to the junction system // (or it was not a junction system). if (int(iPar.size()) == 0) return 0; Vec4 p = event[iPar[0]].p(); for (int i = 1; i < int(iPar.size()); ++i) p += event[iPar[i]].p(); return p.mCalc(); } // ------------------------------------------------------------------ // Find all particles connected to a junction system for junctions stored in // the event record. void ColourReconnection::addJunctionIndices(Event & event, int col, vector &iPar, vector &usedJuncs) { // Find the junction. vector juncs; for (int j = 0;j < event.sizeJunction(); ++j) { for (int k = 0; k < 3; ++k) { if (event.colJunction(j,k) == col) { juncs.push_back(j); break; } } } // Check if they were already used. for (int i = 0;i < int(juncs.size()); ++i) { for (int j = 0;j < int(usedJuncs.size()); ++j) { if (juncs[i] == usedJuncs[j]) { juncs.erase(juncs.begin() + i); i--; break; } } } // If list of junctions is empty return. if (juncs.size() == 0) return; // Store the juncstions as used. for (int i = 0;i < int(juncs.size()); ++i) usedJuncs.push_back(juncs[i]); // Find the partons connected to it. for (int iJunc = 0; iJunc < int(juncs.size()); ++iJunc) { int iTempPars[3] = {-1,-1,-1}; int cols[3] = {event.colJunction(juncs[iJunc],0), event.colJunction(juncs[iJunc],1), event.colJunction(juncs[iJunc],2)}; // Store the first time the colour appear. for (int i = 0;i < event.size(); ++i) { for (int j = 0;j < 3; ++j) { if (iTempPars[j] == -1 && event.kindJunction(juncs[iJunc]) % 2 == 1 && event[i].col() == cols[j]) iTempPars[j] = i; if (iTempPars[j] == -1 && event.kindJunction(juncs[iJunc]) % 2 == 0 && event[i].acol() == cols[j]) iTempPars[j] = i; } } for (int i = 0;i < 3;++i) { if (iTempPars[i] >= 0) iPar.push_back(iTempPars[i]); else addJunctionIndices(event, cols[i], iPar, usedJuncs); } } // Done. } // ------------------------------------------------------------------ // Find all particles connected to a junction system. void ColourReconnection::addJunctionIndices(int iSinglePar, vector &iPar, vector &usedJuncs) { // Check if junction was already considered. int iJun = -(1 + iSinglePar/10); for (int i = 0;i < int(usedJuncs.size()); ++i) if (iJun == usedJuncs[i]) return; // Add particles connected to the junction. usedJuncs.push_back(iJun); for (int i = 0; i < 3; ++i) if (junctions[iJun].kind() % 2 == 1) { int iCol = junctions[iJun].dips[i]->iCol; if (iCol >= 0) iPar.push_back(iCol); else addJunctionIndices(iCol, iPar, usedJuncs); } else { int iAcol = junctions[iJun].dips[i]->iAcol; if (iAcol >= 0) iPar.push_back(iAcol); else addJunctionIndices(iAcol, iPar, usedJuncs); } } // ------------------------------------------------------------------ // Calculate the invariant mass of a dipole. double ColourReconnection::mDip(ColourDipole* dip) { // If double junction no invariant mass is given. if (dip->isJun && dip->isAntiJun) return 1e9; // If it has a single junction end. else if (dip->isJun || dip->isAntiJun) { int iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2; getJunctionIndices(dip, iJun, i0, i1, i2, junLeg0, junLeg1, junLeg2); if (i0 == i1) return particles[i0].m(); if (i1 < 0) return 1e9; return m(particles[i0].p(),particles[i1].p()); } // No junction ends. else { if (dip->iCol == dip->iAcol) return particles[dip->iCol].m(); else return m(particles[dip->iCol].p(),particles[dip->iAcol].p()); } } // ------------------------------------------------------------------ // Print dipoles, intended for debuggning purposes. void ColourReconnection::listDipoles(bool onlyActive, bool onlyReal) { cout << " --- listing dipoles ---" << endl; for (int i = 0; i < int(dipoles.size()); ++i) { if (onlyActive && !dipoles[i]->isActive) continue; if (onlyReal && !dipoles[i]->isReal) continue; dipoles[i]->list(); } cout << " --- finished listing ---" << endl; } // ------------------------------------------------------------------ // Print particles, intended for debugging purposes. void ColourReconnection::listParticles() { for (int i = 0; i < int(particles.size()); ++i) { const ColourParticle& pt = particles[i]; // Basic line for a particle, always printed. cout << setw(6) << i << setw(10) << pt.id() << " " << left << setw(18) << pt.nameWithStatus(18) << right << setw(4) << pt.status() << setw(6) << pt.mother1() << setw(6) << pt.mother2() << setw(6) << pt.daughter1() << setw(6) << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol() << setprecision(3) << setw(11) << pt.px() << setw(11) << pt.py() << setw(11) << pt.pz() << setw(11) << pt.e() << setw(11) << pt.m(); for (int j = 0;j < int(pt.activeDips.size());++j) cout << setw(10) << pt.activeDips[j]; cout << "\n"; } } // ------------------------------------------------------------------ // Print junctions, intended for debugging purposes. void ColourReconnection::listJunctions() { cout << " --- listing junctions ---" << endl; for (int i = 0; i < int(junctions.size()); ++i) junctions[i].list(); cout << " --- finished listing ---" << endl; } // ------------------------------------------------------------------ // Allow colour reconnections by mergings of MPI collision subsystems. // iRec is system that may be reconnected, by moving its gluons to iSys, // where minimal pT (or equivalently Lambda) is used to pick location. // Therefore all dipoles in iSys have to be found, and all gluons in iRec. // Matching q-qbar pairs are treated by analogy with gluons. // Note: owing to rescatterings some outgoing partons must be skipped. bool ColourReconnection::reconnectMPIs( Event& event, int oldSize) { // References to beams to simplify indexing. BeamParticle& beamA = *beamAPtr; BeamParticle& beamB = *beamBPtr; // Prepare record of which systems should be merged onto another. // The iSys system must have colour in final state to attach to it. nSys = partonSystemsPtr->sizeSys(); vector iMerge(nSys); vector hasColour(nSys); for (int iSys = 0; iSys < nSys; ++iSys) { iMerge[iSys] = iSys; bool hasCol = false; for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) { int iNow = partonSystemsPtr->getOut( iSys, iMem); if (event[iNow].isFinal() && (event[iNow].col() > 0 || event[iNow].acol() > 0) ) { hasCol = true; break; } } hasColour[iSys] = hasCol; } // Loop over systems to decide which should be reconnected. for (int iRec = nSys - 1; iRec > 0; --iRec) { // Determine reconnection strength from pT scale of system. double pT2Rec = pow2( partonSystemsPtr->getPTHat(iRec) ); double probRec = pT20Rec / (pT20Rec + pT2Rec); // Loop over other systems iSys at higher pT scale and // decide whether to reconnect the iRec gluons onto one of them. for (int iSys = iRec - 1; iSys >= 0; --iSys) if (hasColour[iSys] && probRec > rndmPtr->flat()) { // The iRec system and all merged with it to be merged with iSys. iMerge[iRec] = iSys; for (int iRec2 = iRec + 1; iRec2 < nSys; ++iRec2) if (iMerge[iRec2] == iRec) iMerge[iRec2] = iSys; // Once a system has been merged do not test it anymore. break; } } // Loop over systems. Check whether other systems to be merged with it. for (int iSys = 0; iSys < nSys; ++iSys) { int nMerge = 0; for (int iRec = iSys + 1; iRec < nSys; ++iRec) if (iMerge[iRec] == iSys) ++nMerge; if (nMerge == 0) continue; // Incoming partons not counted if rescattered. int iInASys = partonSystemsPtr->getInA(iSys); bool hasInA = (beamA[iSys].isFromBeam()); int iInBSys = partonSystemsPtr->getInB(iSys); bool hasInB = (beamB[iSys].isFromBeam()); // Begin find dipoles in iSys system. vector bmdipoles; int sizeOut = partonSystemsPtr->sizeOut(iSys); for (int iMem = 0; iMem < sizeOut; ++iMem) { // Find colour dipoles to beam remnant. int iNow = partonSystemsPtr->getOut( iSys, iMem); if (!event[iNow].isFinal()) continue; int col = event[iNow].col(); if (col > 0) { if (hasInA && event[iInASys].col() == col) bmdipoles.push_back( BeamDipole( col, iNow, iInASys ) ); else if (hasInB && event[iInBSys].col() == col) bmdipoles.push_back( BeamDipole( col, iNow, iInBSys ) ); // Find colour dipole between final-state partons. else for (int iMem2 = 0; iMem2 < sizeOut; ++iMem2) if (iMem2 != iMem) { int iNow2 = partonSystemsPtr->getOut( iSys, iMem2); if (!event[iNow2].isFinal()) continue; if (event[iNow2].acol() == col) { bmdipoles.push_back( BeamDipole( col, iNow, iNow2) ); break; } } } // Find anticolour dipoles to beam remnant. int acol = event[iNow].acol(); if (acol > 0) { if (hasInA && event[iInASys].acol() == acol) bmdipoles.push_back( BeamDipole( acol, iInASys, iNow ) ); else if (hasInB && event[iInBSys].acol() == acol) bmdipoles.push_back( BeamDipole( acol, iInBSys, iNow ) ); } } // Skip mergings if no dipoles found. if (bmdipoles.size() == 0) continue; // Find dipole sizes. for (int iDip = 0; iDip < int(bmdipoles.size()); ++iDip) bmdipoles[iDip].p1p2 = event[bmdipoles[iDip].iCol].p() * event[bmdipoles[iDip].iAcol].p(); // Loop over systems iRec to be merged with iSys. for (int iRec = iSys + 1; iRec < nSys; ++iRec) { if (iMerge[iRec] != iSys) continue; // Information on iRec. Vectors for gluons and anything else. int sizeRec = partonSystemsPtr->sizeOut(iRec); int iInARec = partonSystemsPtr->getInA(iRec); int iInBRec = partonSystemsPtr->getInB(iRec); int nGluRec = 0; vector iGluRec; vector pT2GluRec; int nAnyRec = 0; vector iAnyRec; vector freeAnyRec; // Copy of gluon positions in descending order. for (int iMem = 0; iMem < sizeRec; ++iMem) { int iNow = partonSystemsPtr->getOut( iRec, iMem); if (!event[iNow].isFinal()) continue; if (event[iNow].isGluon()) { ++nGluRec; iGluRec.push_back( iNow ); pT2GluRec.push_back( event[iNow].pT2() ); for (int i = nGluRec - 1; i > 1; --i) { if (pT2GluRec[i - 1] > pT2GluRec[i]) break; swap( iGluRec[i - 1], iGluRec[i] ); swap( pT2GluRec[i - 1], pT2GluRec[i] ); } // Copy of anything else, mainly quarks, in no particular order. } else { ++nAnyRec; iAnyRec.push_back( iNow ); freeAnyRec.push_back( true ); } } // For each gluon in iRec now find the dipole that gives the smallest // (pGlu * pI) (pGlu * pJ) / (pI * pJ), i.e. minimal pT (and Lambda). for (int iGRec = 0; iGRec < nGluRec; ++iGRec) { int iGlu = iGluRec[iGRec]; Vec4 pGlu = event[iGlu].p(); int iDipMin = 0; double pT2DipMin = sCM; for (int iDip = 0; iDip < int(bmdipoles.size()); ++iDip) { double pT2Dip = (pGlu * event[bmdipoles[iDip].iCol].p()) * (pGlu * event[bmdipoles[iDip].iAcol].p()) / bmdipoles[iDip].p1p2; if (pT2Dip < pT2DipMin) { iDipMin = iDip; pT2DipMin = pT2Dip; } } // Attach the gluon to the dipole, i.e. split the dipole in two. int colGlu = event[iGlu].col(); int acolGlu = event[iGlu].acol(); int colDip = bmdipoles[iDipMin].col; int iColDip = bmdipoles[iDipMin].iCol; int iAcolDip = bmdipoles[iDipMin].iAcol; event[iGlu].acol( colDip ); if (event[iAcolDip].acol() == colDip) event[iAcolDip].acol( colGlu ); else event[iAcolDip].col( colGlu ); bmdipoles[iDipMin].iAcol = iGlu; bmdipoles[iDipMin].p1p2 = event[iColDip].p() * pGlu; bmdipoles.push_back( BeamDipole( colGlu, iGlu, iAcolDip ) ); bmdipoles.back().p1p2 = pGlu * event[iAcolDip].p(); // Remove gluon from old system: reconnect colours. for (int i = oldSize; i < event.size(); ++i) if (i != iGlu && i != iAcolDip) { if (event[i].isFinal()) { if (event[i].acol() == colGlu) event[i].acol( acolGlu ); } else { if (event[i].col() == colGlu) event[i].col( acolGlu ); } } // Update any junction legs that match reconnected dipole. for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) { // Only junctions need to be updated, not antijunctions. if (event.kindJunction(iJun) % 2 == 0) continue; for (int leg = 0; leg < 3; ++leg) { int col = event.colJunction(iJun, leg); if (col == colDip) event.colJunction(iJun, leg, colGlu); } } } // See if any matching quark-antiquark pairs among the rest. for (int iQRec = 0; iQRec < nAnyRec; ++iQRec) { int iQ = iAnyRec[iQRec]; int idQ = event[iQ].id(); if (freeAnyRec[iQRec] && idQ > 0 && idQ < 6) for (int iQbarRec = 0; iQbarRec < nAnyRec; ++iQbarRec) { int iQbar = iAnyRec[iQbarRec]; if (freeAnyRec[iQbarRec] && event[iQbar].id() == -idQ) { // Check that these can be traced back to same gluon splitting. // For now also avoid qqbar pairs produced in rescatterings.?? int iTopQ = event[iQ].iTopCopyId(); int iTopQbar = event[iQbar].iTopCopyId(); int iMother = event[iTopQ].mother1(); if (event[iTopQbar].mother1() == iMother && event[iMother].isGluon() && event[iMother].status() != -34 && event[iMother + 1].status() != -34 ) { // Now find the dipole that gives the smallest // ((pQ + pQbar) * pI) ((pQ + pQbar) * pJ) / (pI * pJ). Vec4 pGlu = event[iQ].p() + event[iQbar].p(); int iDipMin = 0; double pT2DipMin = sCM; for (int iDip = 0; iDip < int(bmdipoles.size()); ++iDip) { double pT2Dip = (pGlu * event[bmdipoles[iDip].iCol].p()) * (pGlu * event[bmdipoles[iDip].iAcol].p()) / bmdipoles[iDip].p1p2; if (pT2Dip < pT2DipMin) { iDipMin = iDip; pT2DipMin = pT2Dip; } } // Attach the q-qbar pair to the dipole, i.e. split the dipole. int colGlu = event[iQ].col(); int acolGlu = event[iQbar].acol(); int colDip = bmdipoles[iDipMin].col; int iColDip = bmdipoles[iDipMin].iCol; int iAcolDip = bmdipoles[iDipMin].iAcol; event[iQbar].acol( colDip ); if (event[iAcolDip].acol() == colDip) event[iAcolDip].acol( colGlu ); else event[iAcolDip].col( colGlu ); bmdipoles[iDipMin].iAcol = iQbar; bmdipoles[iDipMin].p1p2 = event[iColDip].p() * event[iQbar].p(); bmdipoles.push_back( BeamDipole( colGlu, iQ, iAcolDip ) ); bmdipoles.back().p1p2 = event[iQ].p() * event[iAcolDip].p(); // Remove q-qbar pair from old system: reconnect colours. freeAnyRec[iQRec] = false; freeAnyRec[iQbarRec] = false; for (int i = oldSize; i < event.size(); ++i) if (i != iQRec && i != iQbarRec && i != iColDip && i != iAcolDip) { if (event[i].isFinal()) { if (event[i].acol() == colGlu) event[i].acol( acolGlu ); } else { if (event[i].col() == colGlu) event[i].col( acolGlu ); } } // Update any junction legs that match reconnected dipole. for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) { // Only junctions need to be updated, not antijunctions. if (event.kindJunction(iJun) % 2 == 0) continue; for (int leg = 0; leg < 3; ++leg) { int col = event.colJunction(iJun, leg); if (col == colDip) event.colJunction(iJun, leg, colGlu); } } // Done with processing of q-qbar pairs. } } } } // If only two beam gluons left of system, set their colour = anticolour. // Used by BeamParticle::remnantColours to skip irrelevant gluons. if ( event[iInARec].isGluon() && !event[iInARec].isRescatteredIncoming() && event[iInBRec].isGluon() && !event[iInBRec].isRescatteredIncoming() && event[iInARec].col() == event[iInBRec].acol() && event[iInARec].acol() == event[iInBRec].col() ) { event[iInARec].acol( event[iInARec].col() ); event[iInBRec].acol( event[iInBRec].col() ); } // End of loops over iRec and iSys systems. } } // Done. return true; } //-------------------------------------------------------------------------- // Find the neighbour to the anticolour side. Return false if the dipole // is connected to a junction or the new particle has a junction inside of it. bool ColourReconnection::findAntiNeighbour(ColourDipole*& dip) { // If only one active dipole, it has to be an antiquark. if (int(particles[dip->iAcol].activeDips.size()) == 1) return false; // Has to have to active dipoles, otherwise something went wrong. if (int(particles[dip->iAcol].activeDips.size()) != 2) { infoPtr->errorMsg("Warning in ColourReconnection::findAntiNeighbour: " "Wrong number of active dipoles"); return false; } // Otherwise find new dipole. if (dip == particles[dip->iAcol].activeDips[0]) dip = particles[dip->iAcol].activeDips[1]; else dip = particles[dip->iAcol].activeDips[0]; // Do not allow the new dipole to be connected to an antijunction. if (dip->isAntiJun || dip->isJun) return false; // Do not allow new dipole to have a pseudoparticle with // a baryon number inside. if (int(particles[dip->iAcol].dips.size()) != 1) return false; return true; } //-------------------------------------------------------------------------- // Check that trials do not contain junctions/ unusable pseudoparticles. bool ColourReconnection::checkJunctionTrials() { for (int i = 0;i < int(junTrials.size());++i) { int minus = 0; if (junTrials[i].mode == 3) minus = 1; for (int j = 0;j < int(junTrials[i].dips.size()) - minus; ++j) { ColourDipole* dip = junTrials[i].dips[j]; if (dip->isJun || dip->isAntiJun) { junTrials[i].list(); return false; } if (particles[dip->iCol].dips.size() != 1 || particles[dip->iAcol].dips.size() != 1) { junTrials[i].list(); return false; } } } return true; } //-------------------------------------------------------------------------- // Find the neighbour to the colour side. Return false if the dipole // is connected to a junction or the new particle has a junction inside of it. bool ColourReconnection::findColNeighbour(ColourDipole*& dip) { // If only one active dipole, it has to be an antiquark. if (int(particles[dip->iCol].activeDips.size()) == 1) return false; // Has to have to active dipoles, otherwise something went wrong. if (int(particles[dip->iCol].activeDips.size()) != 2) { infoPtr->errorMsg("Warning in ColourReconnection::findAntiNeighbour: " "Wrong number of active dipoles"); return false; } // Otherwise find new dipole. if (dip == particles[dip->iCol].activeDips[0]) dip = particles[dip->iCol].activeDips[1]; else dip = particles[dip->iCol].activeDips[0]; // Do not allow the new dipole to be connected to an antijunction. if (dip->isJun || dip->isAntiJun) return false; // Do not allow new dipole to have a pseudoparticle with // a baryon number inside. if (int(particles[dip->iCol].dips.size()) != 1) return false; return true; } //-------------------------------------------------------------------------- // Store used dipoles for a junction formation. void ColourReconnection::storeUsedDips(TrialReconnection& trial) { // Normal dipole swap. if (trial.mode == 5) { for (int i = 0;i < 2; ++i) { ColourDipole* dip = trial.dips[i]; if (dip->iCol < 0) for (int j = 0;j < 3; ++j) usedDipoles.push_back(junctions[-(dip->iCol / 10 + 1)].getColDip(j)); if (dip->iAcol < 0) for (int j = 0;j < 3; ++j) usedDipoles.push_back(junctions[-(dip->iAcol / 10 + 1)].getColDip(j)); usedDipoles.push_back(dip); } } else { for (int i = 0;i < 4; ++i) { if (trial.mode == 3 && i == 3) continue; usedDipoles.push_back(trial.dips[i]); ColourDipole* dip = trial.dips[i]; while (findAntiNeighbour(dip) && dip != trial.dips[i]) usedDipoles.push_back(dip); dip = trial.dips[i]; while (findColNeighbour(dip) && dip != trial.dips[i]) usedDipoles.push_back(dip); } } } //-------------------------------------------------------------------------- // Calculate the difference between the old and new lambda for dipole swap. double ColourReconnection::getLambdaDiff(ColourDipole* dip1, ColourDipole* dip2) { // Needed to make sure the same dipoles are compared. vector oldDips, newDips; // Calculate old string length. double oldLambda = calculateStringLength(dip1, oldDips) + calculateStringLength( dip2, oldDips); // Make test configuration. swapDipoles(dip1,dip2); // Calculate new string lengths double newLambda = calculateStringLength(dip1, newDips) + calculateStringLength(dip2, newDips); // Swap back. swapDipoles(dip1, dip2, true); // First check if new combination was not useable. if (newLambda >= 0.5E9) return -1e9; // Return the difference. return oldLambda - newLambda; } //-------------------------------------------------------------------------- // Calculate the difference between the old and new lambda. double ColourReconnection::getLambdaDiff(ColourDipole* dip1, ColourDipole* dip2, ColourDipole* dip3, ColourDipole* dip4, int mode) { // Calculate old lambda measure. double oldLambda = calculateStringLength(dip1->iCol, dip1->iAcol) + calculateStringLength(dip2->iCol, dip2->iAcol); if (dip3 != dip1) oldLambda += calculateStringLength(dip3->iCol, dip3->iAcol); if (dip4 != 0 && dip4 != dip2) oldLambda += calculateStringLength(dip4->iCol, dip4->iAcol); // Calculate new lambda. double newLambda = 0; if (mode == 0) newLambda = calculateDoubleJunctionLength(dip1->iCol, dip2->iCol, dip1->iAcol, dip2->iAcol); else if (mode == 1) { if (dip2 == dip4) newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip3->iCol) + calculateJunctionLength(dip1->iAcol, dip2->iAcol, dip3->iAcol); else newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip3->iCol) + calculateJunctionLength(dip2->iAcol, dip3->iAcol, dip4->iAcol) + calculateStringLength(dip4->iCol, dip1->iAcol); } else if (mode == 2) { if (dip1 == dip3) newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip4->iCol) + calculateJunctionLength(dip1->iAcol, dip2->iAcol, dip4->iAcol); else newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip4->iCol) + calculateJunctionLength(dip1->iAcol, dip3->iAcol, dip4->iAcol) + calculateStringLength(dip3->iCol, dip2->iAcol); } // Triple junction connection. else if (mode == 3) newLambda = calculateJunctionLength(dip1->iCol, dip2->iCol, dip3->iCol) + calculateJunctionLength(dip1->iAcol, dip2->iAcol, dip3->iAcol); // First check if new combination was not useable. if (newLambda >= 0.5E9) return -1e9; // Returning result. return oldLambda - newLambda; } //-------------------------------------------------------------------------- // Change colour structure to describe the reconnection in juncTrial. void ColourReconnection::doDipoleTrial(TrialReconnection& trial) { // Store for easier use. ColourDipole* dip1 = trial.dips[0]; ColourDipole* dip2 = trial.dips[1]; // If both acols ends are normal particles. if (dip1->iAcol >= 0 && dip2->iAcol >= 0) { swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front()->iAcol, particles[dip2->iAcol].dips[dip2->iAcolLeg].front()->iAcol); swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front(), particles[dip2->iAcol].dips[dip2->iAcolLeg].front()); // If only dip1 has normal acol end. } else if (dip1->iAcol >= 0) { swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front()->iAcol, junctions[-(dip2->iAcol / 10 + 1)].dipsOrig[-dip2->iAcol % 10]->iAcol); swap(particles[dip1->iAcol].dips[dip1->iAcolLeg].front(), junctions[-(dip2->iAcol / 10 + 1)].dipsOrig[-dip2->iAcol % 10]); // If only dip2 has normal acol end. } else if (dip2->iAcol >= 0) { swap(particles[dip2->iAcol].dips[dip2->iAcolLeg].front()->iAcol, junctions[-(dip1->iAcol / 10 + 1)].dipsOrig[-dip1->iAcol % 10]->iAcol); swap(particles[dip2->iAcol].dips[dip2->iAcolLeg].front(), junctions[-(dip1->iAcol / 10 + 1)].dipsOrig[-dip1->iAcol % 10]); // If both ends are junctions. } else { swap(junctions[ -(dip1->iAcol / 10 + 1) ].dipsOrig[ -dip1->iAcol % 10 ]->iAcol, junctions[ -(dip2->iAcol / 10 + 1) ].dipsOrig[ -dip2->iAcol % 10 ]->iAcol); swap(junctions[ -(dip1->iAcol / 10 + 1) ].dipsOrig[ -dip1->iAcol % 10], junctions[ -(dip2->iAcol / 10 + 1) ].dipsOrig[ -dip2->iAcol % 10] ); } // Swap the dipoles. swapDipoles(dip1, dip2); // If new particles are below treshhold, form pseudoParticles. if (mDip(dip1) < m0) makePseudoParticle(dip1, 110, true); if (mDip(dip2) < m0) makePseudoParticle(dip2, 110, true); // Done. } //-------------------------------------------------------------------------- // Update the list of dipole trial swaps to account for latest swap. void ColourReconnection::updateDipoleTrials() { // Remove any dipTrials that contains an used dipole. for (int i = 0; i < int(dipTrials.size()); ++i) for (int j = 0;j < 2; ++j) { if (binary_search(usedDipoles.begin(), usedDipoles.end(), dipTrials[i].dips[j]) ) { dipTrials.erase(dipTrials.begin() + i); i--; break; } } // Make list of active dipoles. vector activeDipoles; for (int i = 0;i < int(dipoles.size()); ++i) if (dipoles[i]->isActive) activeDipoles.push_back(dipoles[i]); // Loop over list of used dipoles and create new trial reconnections. for (int i = 0;i < int(usedDipoles.size()); ++i) if (usedDipoles[i]->isActive) for (int j = 0; j < int(activeDipoles.size()); ++j) singleReconnection(usedDipoles[i], activeDipoles[j]); } //-------------------------------------------------------------------------- // Update the list of dipole trial swaps to account for latest swap. void ColourReconnection::updateJunctionTrials() { // Remove any junTrials that contains an used dipole. for (int i = 0; i < int(junTrials.size()); ++i) for (int j = 0; j < 4; ++j) { if (binary_search(usedDipoles.begin(), usedDipoles.end(), junTrials[i].dips[j]) ) { junTrials.erase(junTrials.begin() + i); i--; break; } } // Make list of active dipoles. vector activeDipoles; for (int i = 0;i < int(dipoles.size()); ++i) if (dipoles[i]->isActive) activeDipoles.push_back(dipoles[i]); // Loop over used dipoles and form new junction trials. for (int i = 0;i < int(usedDipoles.size()); ++i) if (usedDipoles[i]->isActive) for (int j = 0; j < int(activeDipoles.size()); ++j) singleJunction(usedDipoles[i], activeDipoles[j]); // Loop over used dipoles and form new junction trials. for (int i = 0;i < int(usedDipoles.size()); ++i) if (usedDipoles[i]->isActive) for (int j = 0; j < int(activeDipoles.size()); ++j) for (int k = j + 1; k < int(activeDipoles.size()); ++k) singleJunction(usedDipoles[i], activeDipoles[j], activeDipoles[k]); } //-------------------------------------------------------------------------- // Change colour structure to describe the reconnection in juncTrial. void ColourReconnection::doJunctionTrial(Event& event, TrialReconnection& juncTrial) { int mode = juncTrial.mode; // If trial mode is 3 (three dipoles -> 2 junctions) use its own update. if (mode == 3) { doTripleJunctionTrial(event, juncTrial); return; } // Store dipoles and numbers for easier acces. ColourDipole* dip1 = juncTrial.dips[0]; ColourDipole* dip2 = juncTrial.dips[1]; ColourDipole* dip3 = juncTrial.dips[2]; ColourDipole* dip4 = juncTrial.dips[3]; int iCol1 = dip1->iCol; int iCol2 = dip2->iCol; int iCol3 = dip3->iCol; int iCol4 = dip4->iCol; int iAcol1 = dip1->iAcol; int iAcol2 = dip2->iAcol; int iAcol3 = dip3->iAcol; int iAcol4 = dip4->iAcol; int oldCol1 = dip1->col; int oldCol2 = dip2->col; int oldCol3 = dip3->col; int oldCol4 = dip4->col; // New colour tags needed, since three more dipoles will be made. int newCol1 = event.nextColTag(); int newCol2 = event.nextColTag(); int newCol3 = event.nextColTag(); // Add new formation times. double mCalc = (particles[iCol1].p() + particles[iAcol1].p() + particles[iCol2].p() + particles[iAcol2].p() + particles[iCol3].p() + particles[iAcol3].p() + particles[iCol4].p() + particles[iAcol4].p()).mCalc(); formationTimes[newCol1] = mCalc; formationTimes[newCol2] = mCalc; formationTimes[newCol3] = mCalc; // Need to make 3 new real dipoles and 3 active dipoles. // First make dipoles between junctions. // Find the junction colour. int junCol = 3 * (3 - (dip1->colReconnection / 3) - (dip2->colReconnection / 3) ) + dip1->colReconnection % 3; // if other than 9 colours. if (nReconCols != 9) { while (junCol < 0 || junCol % 3 != dip1->colReconnection % 3 || junCol == dip1->colReconnection || junCol == dip2->colReconnection) junCol = int(nReconCols * rndmPtr->flat()); } // Need one active and one real dipole. int iJun = junctions.size(); int iAntiJun = junctions.size() + 1; // Store real dipoles. ColourDipole* dip3real = particles[iCol3].dips[dip3->iColLeg].back(); ColourDipole* dip4real = particles[iCol4].dips[dip4->iColLeg].back(); // If the junction and antijunction are directly connected. int iActive1 = 0, iReal1 = 0; if (mode == 0) { dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10 + 2) , -( iJun * 10 + 10 + 2), junCol, true, true, false, true)); iReal1 = dipoles.size() - 1; dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10 + 2) , -( iJun * 10 + 10 + 2), junCol, true, true)); iActive1 = dipoles.size() - 1; } else if (mode == 1) { int iCol3real = particles[iCol3].dips[dip3->iColLeg].back()->iCol; dipoles.push_back(new ColourDipole(newCol1, iCol3real , -( iJun * 10 + 10 + 2), junCol, true, false, false, true)); iReal1 = dipoles.size() - 1; particles[iCol3].dips[dip3->iColLeg].back() = dipoles.back(); dipoles.push_back(new ColourDipole(newCol1, dip3->iCol, -( iJun * 10 + 10 + 2), junCol, true, false)); iActive1 = dipoles.size() - 1; } else if (mode == 2) { int iCol4real = particles[iCol4].dips[dip4->iColLeg].back()->iCol; dipoles.push_back(new ColourDipole(newCol1, iCol4real, -( iJun * 10 + 10 + 2), junCol, true, false, false, true)); iReal1 = dipoles.size() - 1; particles[iCol4].dips[dip4->iColLeg].back() = dipoles.back(); dipoles.push_back(new ColourDipole(newCol1, dip4->iCol, -( iJun * 10 + 10 + 2), junCol, true, false)); iActive1 = dipoles.size() - 1; } // Now make dipole between antijunction and iAcol1. // Start by finding real iAcol. int iAcol3real = particles[iAcol3].dips[dip3->iAcolLeg].front()->iAcol; dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10), iAcol3real, dip3->colReconnection, false, true, false, true)); int iReal2 = dipoles.size() - 1; particles[iAcol3].dips[dip3->iAcolLeg].front() = dipoles.back(); dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10), iAcol3, dip3->colReconnection, false, true)); dipoles.back()->iAcolLeg = dip3->iAcolLeg; int iActive2 = dipoles.size() - 1; // Now make dipole between antijunction and iAcol1. // Start by finding real iAcol. int iAcol4real = particles[iAcol4].dips[dip4->iAcolLeg].front()->iAcol; dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 1), iAcol4real, dip4->colReconnection, false, true, false, true)); int iReal3 = dipoles.size() - 1; particles[iAcol4].dips[dip4->iAcolLeg].front() = dipoles.back(); dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 1), iAcol4, dip4->colReconnection, false, true)); dipoles.back()->iAcolLeg = dip4->iAcolLeg; int iActive3 = dipoles.size() - 1; // Update already existing dipoles, start by internal dipoles. // Now take dipoles connected to the antijunction // and a possible gluon-gluon connection. if (mode == 1) { if (dip2 == dip4) { // Update real dipole. dip3real->iAcol = particles[dip1->iAcol].dips[dip1->iAcolLeg]. front()->iAcol; dip3real->iCol = -( iAntiJun * 10 + 10 + 2); dip3real->isAntiJun = true; // Update active dipoles. dip3->iAcol = dip1->iAcol; dip3->iAcolLeg = dip1->iAcolLeg; dip3->isAntiJun = true; dip3->iCol = -( iAntiJun * 10 + 10 + 2); dip3->iColLeg = 0; // Store real dipole particles[dip3->iAcol].dips[dip3->iAcolLeg].front() = dip3real; } else { // Update real dipole. dip3real->iAcol = particles[dip2->iAcol].dips[dip2->iAcolLeg]. front()->iAcol; dip3real->iCol = -( iAntiJun * 10 + 10 + 2); dip3real->isAntiJun = true; dip4real->iAcol = particles[dip1->iAcol].dips[dip1->iAcolLeg]. front()->iAcol; // Change the dipole connected to the antijunction. dip3->iAcol = dip2->iAcol; dip3->iAcolLeg = dip2->iAcolLeg; dip3->isAntiJun = true; dip3->iCol = -( iAntiJun * 10 + 10 + 2); dip3->iColLeg = 0; // Change the dipole between the two gluons. dip4->iAcol = dip1->iAcol; dip4->iAcolLeg = dip1->iAcolLeg; // Store real dipole particles[dip3->iAcol].dips[dip3->iAcolLeg].front() = dip3real; particles[dip4->iAcol].dips[dip4->iAcolLeg].front() = dip4real; } } else if (mode == 2) { if (dip1 == dip3) { // Update real dipole. dip4real->iAcol = particles[dip2->iAcol].dips[dip2->iAcolLeg]. front()->iAcol; dip4real->iCol = -( iAntiJun * 10 + 10 + 2); dip4real->isAntiJun = true; // Update active dipoles. dip4->iAcol = dip2->iAcol; dip4->iAcolLeg = dip2->iAcolLeg; dip4->isAntiJun = true; dip4->iCol = -( iAntiJun * 10 + 10 + 2); dip4->iColLeg = 0; // Store real dipole particles[dip4->iAcol].dips[dip4->iAcolLeg].front() = dip4real; } else { // Update real dipole. dip4real->iAcol = particles[dip1->iAcol].dips[dip1->iAcolLeg]. front()->iAcol; dip4real->iCol = -( iAntiJun * 10 + 10 + 2); dip4real->isAntiJun = true; dip3real->iAcol = particles[dip2->iAcol].dips[dip2->iAcolLeg]. front()->iAcol; // Change the dipole connected to the antijunction. dip4->iAcol = dip1->iAcol; dip4->iAcolLeg = dip1->iAcolLeg; dip4->isAntiJun = true; dip4->iCol = -( iAntiJun * 10 + 10 + 2); dip4->iColLeg = 0; // Change the dipole between the two gluons. dip3->iAcol = dip2->iAcol; dip3->iAcolLeg = dip2->iAcolLeg; // Store real dipole particles[dip3->iAcol].dips[dip3->iAcolLeg].front() = dip3real; particles[dip4->iAcol].dips[dip4->iAcolLeg].front() = dip4real; } } // Dipoles connected to the junction. // Update real dipoles. particles[iCol1].dips[dip1->iColLeg].back()->iAcol = - (iJun * 10 + 10); particles[iCol2].dips[dip2->iColLeg].back()->iAcol = - (iJun * 10 + 10 + 1); particles[iCol1].dips[dip1->iColLeg].back()->isJun = true; particles[iCol2].dips[dip2->iColLeg].back()->isJun = true; // Update active dipoles. dip1->isJun = true; dip2->isJun = true; dip1->iAcol = - (iJun * 10 + 10); dip2->iAcol = - (iJun * 10 + 10 + 1); dip1->iAcolLeg = 0; dip2->iAcolLeg = 0; // Update active dipoles for anti particles. // Normally should only contain active dipoles once, // only problem is if the two dipole ends are the same particle. // Start by settings common dipoles. for (int i = 0; i < int(particles[iAcol3].activeDips.size()); ++i) if (particles[iAcol3].activeDips[i] == dip3) { particles[iAcol3].activeDips[i] = dipoles[iActive2]; break; } for (int i = 0; i < int(particles[iAcol4].activeDips.size()); ++i) if (particles[iAcol4].activeDips[i] == dip4) { particles[iAcol4].activeDips[i] = dipoles[iActive3]; break; } // Depending on how the new string is connected, the active dipoles vary. if (mode == 1) { for (int i = 0; i < int(particles[iCol3].activeDips.size()); ++i) if (particles[iCol3].activeDips[i] == dip3) { particles[iCol3].activeDips[i] = dipoles[iActive1]; break; } if (dip2 == dip4) { for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i) if (particles[iAcol1].activeDips[i] == dip1) { particles[iAcol1].activeDips[i] = dip3; break; } } else { for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i) if (particles[iAcol2].activeDips[i] == dip2) { particles[iAcol2].activeDips[i] = dip3; break; } for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i) if (particles[iAcol1].activeDips[i] == dip1) { particles[iAcol1].activeDips[i] = dip4; break; } } } else if (mode == 2) { for (int i = 0; i < int(particles[iCol4].activeDips.size()); ++i) if (particles[iCol4].activeDips[i] == dip4) { particles[iCol4].activeDips[i] = dipoles[iActive1]; break; } if (dip1 == dip3) { for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i) if (particles[iAcol2].activeDips[i] == dip2) { particles[iAcol2].activeDips[i] = dip4; break; } } else { for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i) if (particles[iAcol1].activeDips[i] == dip1) { particles[iAcol1].activeDips[i] = dip4; break; } for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i) if (particles[iAcol2].activeDips[i] == dip2) { particles[iAcol2].activeDips[i] = dip3; break; } } } // Add the junctions to the event. junctions.push_back(Junction(1, oldCol1, oldCol2, newCol1)); if (mode == 0) junctions.push_back(Junction(2, newCol2, newCol3, newCol1)); else if (mode == 1) junctions.push_back(Junction(2, newCol2, newCol3, oldCol3)); else if (mode == 2) junctions.push_back(Junction(2, newCol2, newCol3, oldCol4)); // Set junction information. junctions[iJun].dipsOrig[0] = particles[iCol1].dips[dip1->iColLeg].back(); junctions[iJun].dipsOrig[1] = particles[iCol2].dips[dip2->iColLeg].back(); junctions[iJun].dipsOrig[2] = dipoles[iReal1]; junctions[iJun].dips[0] = dip1; junctions[iJun].dips[1] = dip2; junctions[iJun].dips[2] = dipoles[iActive1]; // Set antijunction information. junctions[iAntiJun].dips[0] = dipoles[iActive2]; junctions[iAntiJun].dips[1] = dipoles[iActive3]; junctions[iAntiJun].dipsOrig[0] = dipoles[iReal2]; junctions[iAntiJun].dipsOrig[1] = dipoles[iReal3]; if (mode == 0) { junctions[iAntiJun].dips[2] = dipoles[iActive1]; junctions[iAntiJun].dipsOrig[2] = dipoles[iReal1]; } else if (mode == 1) { junctions[iAntiJun].dips[2] = dip3; junctions[iAntiJun].dipsOrig[2] = particles[dip3->iAcol].dips[dip3->iAcolLeg].front(); } else if (mode == 2) { junctions[iAntiJun].dips[2] = dip4; junctions[iAntiJun].dipsOrig[2] = particles[dip4->iAcol].dips[dip4->iAcolLeg].front(); } // Make pseudo particles. if (dip1->isActive && mDip(dip1) < m0) makePseudoParticle(dip1, 110, true); if (dip2->isActive && mDip(dip2) < m0) makePseudoParticle(dip2, 110, true); if (dip3->isActive && mDip(dip3) < m0) makePseudoParticle(dip3, 110, true); if (dip4->isActive && mDip(dip4) < m0) makePseudoParticle(dip4, 110, true); if (dipoles[iActive1]->isActive && mDip(dipoles[iActive1]) < m0) makePseudoParticle(dipoles[iActive1], 110, true); if (dipoles[iActive2]->isActive && mDip(dipoles[iActive2]) < m0) makePseudoParticle(dipoles[iActive2], 110, true); if (dipoles[iActive3]->isActive && mDip(dipoles[iActive3]) < m0) makePseudoParticle(dipoles[iActive3], 110, true); // Add new dipoles to usedDipoles. usedDipoles.push_back(dipoles[iActive1]); usedDipoles.push_back(dipoles[iActive2]); usedDipoles.push_back(dipoles[iActive3]); // Done. } //-------------------------------------------------------------------------- // Change colour structure if it is three dipoles forming a junction system. void ColourReconnection::doTripleJunctionTrial(Event& event, TrialReconnection& juncTrial) { // store information for easier acces. ColourDipole* dip1 = juncTrial.dips[0]; ColourDipole* dip2 = juncTrial.dips[1]; ColourDipole* dip3 = juncTrial.dips[2]; // Store indices. int iCol1 = dip1->iCol; int iCol2 = dip2->iCol; int iCol3 = dip3->iCol; int iAcol1 = dip1->iAcol; int iAcol2 = dip2->iAcol; int iAcol3 = dip3->iAcol; // Store colours int oldCol1 = dip1->col; int oldCol2 = dip2->col; int oldCol3 = dip3->col; int newCol1 = event.nextColTag(); int newCol2 = event.nextColTag(); int newCol3 = event.nextColTag(); // Store new junction indices. int iJun = junctions.size(); int iAntiJun = junctions.size() + 1; // Now make dipole between antijunction and iAcol1. // Start by finding real iAcol. int iAcol1real = particles[iAcol1].dips[dip1->iAcolLeg].front()->iAcol; dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10), iAcol1real, dip1->colReconnection, false, true, false, true)); int iReal1 = dipoles.size() - 1; particles[iAcol1].dips[dip1->iAcolLeg].front() = dipoles.back(); dipoles.push_back(new ColourDipole(newCol1, -( iAntiJun * 10 + 10), iAcol1, dip1->colReconnection, false, true)); dipoles.back()->iAcolLeg = dip1->iAcolLeg; int iActive1 = dipoles.size() - 1; // Now make dipole between antijunction and iAcol2. // Start by finding real iAcol2. int iAcol2real = particles[iAcol2].dips[dip2->iAcolLeg].front()->iAcol; dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10 + 1), iAcol2real, dip2->colReconnection, false, true, false, true)); int iReal2 = dipoles.size() - 1; particles[iAcol2].dips[dip2->iAcolLeg].front() = dipoles.back(); dipoles.push_back(new ColourDipole(newCol2, -( iAntiJun * 10 + 10 + 1), iAcol2, dip2->colReconnection, false, true)); dipoles.back()->iAcolLeg = dip2->iAcolLeg; int iActive2 = dipoles.size() - 1; // Now make dipole between antijunction and iAcol3. // Start by finding real iAcol3. int iAcol3real = particles[iAcol3].dips[dip3->iAcolLeg].front()->iAcol; dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 2), iAcol3real, dip3->colReconnection, false, true, false, true)); int iReal3 = dipoles.size() - 1; particles[iAcol3].dips[dip3->iAcolLeg].front() = dipoles.back(); dipoles.push_back(new ColourDipole(newCol3, -( iAntiJun * 10 + 10 + 2), iAcol3, dip3->colReconnection, false, true)); dipoles.back()->iAcolLeg = dip3->iAcolLeg; int iActive3 = dipoles.size() - 1; // Update already existing dipoles. // Update real dipoles. particles[iCol1].dips[dip1->iColLeg].back()->iAcol = - (iJun * 10 + 10); particles[iCol2].dips[dip2->iColLeg].back()->iAcol = - (iJun * 10 + 10 + 1); particles[iCol3].dips[dip3->iColLeg].back()->iAcol = - (iJun * 10 + 10 + 2); particles[iCol1].dips[dip1->iColLeg].back()->isJun = true; particles[iCol2].dips[dip2->iColLeg].back()->isJun = true; particles[iCol3].dips[dip3->iColLeg].back()->isJun = true; // Update active dipoles. dip1->isJun = true; dip2->isJun = true; dip3->isJun = true; dip1->iAcol = - (iJun * 10 + 10); dip2->iAcol = - (iJun * 10 + 10 + 1); dip3->iAcol = - (iJun * 10 + 10 + 2); dip1->iAcolLeg = 0; dip2->iAcolLeg = 0; dip3->iAcolLeg = 0; // Update active dipoles for anti particles. for (int i = 0; i < int(particles[iAcol1].activeDips.size()); ++i) if (particles[iAcol1].activeDips[i] == dip1) particles[iAcol1].activeDips[i] = dipoles[iActive1]; for (int i = 0; i < int(particles[iAcol2].activeDips.size()); ++i) if (particles[iAcol2].activeDips[i] == dip2) particles[iAcol2].activeDips[i] = dipoles[iActive2]; for (int i = 0; i < int(particles[iAcol3].activeDips.size()); ++i) if (particles[iAcol3].activeDips[i] == dip3) particles[iAcol3].activeDips[i] = dipoles[iActive3]; // Add the junctions to the event. junctions.push_back(Junction(1, oldCol1, oldCol2, oldCol3)); junctions.push_back(Junction(2, newCol1, newCol3, newCol3)); // Update junction ends. junctions[iJun].dipsOrig[0] = particles[iCol1].dips[dip1->iColLeg].back(); junctions[iJun].dipsOrig[1] = particles[iCol2].dips[dip2->iColLeg].back(); junctions[iJun].dipsOrig[2] = particles[iCol3].dips[dip3->iColLeg].back(); junctions[iJun].dips[0] = dip1; junctions[iJun].dips[1] = dip2; junctions[iJun].dips[2] = dip3; // Update the antijunction. junctions[iAntiJun].dips[0] = dipoles[iActive1]; junctions[iAntiJun].dips[1] = dipoles[iActive2]; junctions[iAntiJun].dips[2] = dipoles[iActive3]; junctions[iAntiJun].dipsOrig[0] = dipoles[iReal1]; junctions[iAntiJun].dipsOrig[1] = dipoles[iReal2]; junctions[iAntiJun].dipsOrig[2] = dipoles[iReal3]; // Make pseudo particles if needed. if (dip1->isActive && mDip(dip1) < m0) makePseudoParticle(dip1, 110, true); if (dip2->isActive && mDip(dip2) < m0) makePseudoParticle(dip2, 110, true); if (dip3->isActive && mDip(dip3) < m0) makePseudoParticle(dip3, 110, true); if (dipoles[iActive1]->isActive && mDip(dipoles[iActive1]) < m0) makePseudoParticle(dipoles[iActive1], 110, true); if (dipoles[iActive2]->isActive && mDip(dipoles[iActive2]) < m0) makePseudoParticle(dipoles[iActive2], 110, true); if (dipoles[iActive3]->isActive && mDip(dipoles[iActive3]) < m0) makePseudoParticle(dipoles[iActive3], 110, true); // Add to newly created dipoles to used dipoles. usedDipoles.push_back(dipoles[iActive1]); usedDipoles.push_back(dipoles[iActive2]); usedDipoles.push_back(dipoles[iActive3]); // Done. } //-------------------------------------------------------------------------- // Allow colour reconnections by moving gluons from their current location // to another colour line. Also optionally flip two colour chains. bool ColourReconnection::reconnectMove( Event& event, int oldSize) { // Create or reset arrays to prepare for the new event analysis. vector iGlu; iReduceCol.resize( event.size() ); iExpandCol.clear(); map colMap, acolMap; map::iterator colM, acolM; vector infoGM; // Temporary variables. int iNow = 0; int colNow = 0; int acolNow = 0; int iColNow = 0; int iAcolNow = 0; int col2Now = 0; int iCol2Now = 0; int iAcol2Now = 0; double lambdaRefNow = 0.; double dLambdaNow = 0.; // Loop over all final particles. Store (fraction of) gluons to move. for (int i = oldSize; i < event.size(); ++i) if (event[i].isFinal()) { if (flipMode < 3 && event[i].id() == 21 && rndmPtr->flat() < fracGluon) iGlu.push_back(i); // Store location of all colour and anticolour particles and indices. if (event[i].col() > 0 || event[i].acol() > 0) { iReduceCol[i] = iExpandCol.size(); iExpandCol.push_back(i); if (event[i].col() > 0) colMap[event[i].col()] = i; if (event[i].acol() > 0) acolMap[event[i].acol()] = i; } } // Erase (anti)colours for (anti)junctions and skip adjacent gluons. for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) { if (event.kindJunction(iJun) == 1) { for (int j = 0; j < 3; ++j) { int jCol = event.colJunction( iJun, j); for (colM = colMap.begin(); colM != colMap.end(); ++colM) if (colM->first == jCol) { colMap.erase( colM); break; } for (int iG = 0; iG < int(iGlu.size()); ++iG) if (event[iGlu[iG]].col() == jCol) { iGlu.erase(iGlu.begin() + iG); break; } } } else if (event.kindJunction(iJun) == 2) { for (int j = 0; j < 3; ++j) { int jCol = event.colJunction( iJun, j); for (acolM = acolMap.begin(); acolM != acolMap.end(); ++acolM) if (acolM->first == jCol) { acolMap.erase( acolM); break; } for (int iG = 0; iG < int(iGlu.size()); ++iG) if (event[iGlu[iG]].acol() == jCol) { iGlu.erase(iGlu.begin() + iG); break; } } } } // Error checks. int nGlu = iGlu.size(); int nCol = colMap.size(); if (int(acolMap.size()) != nCol) { infoPtr->errorMsg("Error in MBReconUserHooks: map sizes do not match"); return false; } colM = colMap.begin(); acolM = acolMap.begin(); for (int iCol = 0; iCol < nCol; ++iCol) { if (colM->first != acolM->first) { infoPtr->errorMsg("Error in MBReconUserHooks: map elements" " do not match"); return false; } ++colM; ++acolM; } // Calculate and tabulate lambda between any pair of coloured partons. nColMove = iExpandCol.size(); lambdaijMove.resize( pow2(nColMove) ); for (int iAC = 0; iAC < nColMove - 1; ++iAC) { int i = iExpandCol[iAC]; for (int jAC = iAC + 1; jAC < nColMove; ++jAC) { int j = iExpandCol[jAC]; lambdaijMove[nColMove * iAC + jAC] = log(1. + m2( event[i], event[j]) / m2Lambda); } } // Set up initial possible gluon moves with lambda gains/losses. for (int iG = 0; iG < nGlu; ++iG) { // Gluon and its current neighbours. iNow = iGlu[iG]; colNow = event[iNow].col(); acolNow = event[iNow].acol(); iColNow = acolMap[colNow]; iAcolNow = colMap[acolNow]; // Addition to Lambda of gluon in current position. lambdaRefNow = lambda123Move( iNow, iColNow, iAcolNow); // Loop over all colour lines where gluon could be inserted. for (colM = colMap.begin(); colM != colMap.end(); ++colM) { col2Now = colM->first; iCol2Now = colMap[col2Now]; iAcol2Now = acolMap[col2Now]; // Addition to total Lambda if gluon moved to be inserted on line. dLambdaNow = (iCol2Now == iNow || iAcol2Now == iNow || iColNow == iAcolNow) ? 2e4 : lambda123Move( iNow, iCol2Now, iAcol2Now) - lambdaRefNow; // Add new container for gluon and colour line information. infoGM.push_back( InfoGluonMove( iNow, colNow, acolNow, iColNow, iAcolNow, col2Now, iCol2Now, iAcol2Now, lambdaRefNow, dLambdaNow )); } } int nPair = infoGM.size(); // Keep on looping over moves until no further negative dLambda. for ( int iMove = 0; iMove < nGlu; ++iMove) { int iPairMin = -1; double dLambdaMin = 1e4; // Find lowest dLambda. for (int iPair = 0; iPair < nPair; ++iPair) if (infoGM[iPair].dLambda < dLambdaMin) { iPairMin = iPair; dLambdaMin = infoGM[iPair].dLambda; } // Break if no shift below upper limit found. if (dLambdaMin > -dLambdaCut) break; // Partons and colours involved in move. InfoGluonMove& selSM = infoGM[iPairMin]; int i1Sel = selSM.i1; int iCol1Sel = selSM.iCol1; int iAcol1Sel = selSM.iAcol1; int iCol2Sel = selSM.iCol2; int iAcol2Sel = selSM.iAcol2; int iCol2Mod[3] = { iAcol1Sel , i1Sel , iCol2Sel }; int col2Mod[3] = { selSM.col1, selSM.col2, selSM.acol1}; // Remove gluon from old colour line and insert on new colour line. for (int i = 0; i < 3; ++i) { event[ iCol2Mod[i] ].col( col2Mod[i] ); colMap[ col2Mod[i] ] = iCol2Mod[i]; } // Update info for partons with new colors. int i1Now = 0; bool doUpdate = false; for (int iPair = 0; iPair < nPair; ++iPair) { InfoGluonMove& tmpSM = infoGM[iPair]; if (tmpSM.i1 != i1Now) { i1Now = tmpSM.i1; doUpdate = false; if (i1Now == i1Sel || i1Now == iCol1Sel || i1Now == iAcol1Sel || i1Now == iCol2Sel || i1Now == iAcol2Sel) { colNow = event[i1Now].col(); acolNow = event[i1Now].acol(); iColNow = acolMap[colNow]; iAcolNow = colMap[acolNow]; lambdaRefNow = lambda123Move( i1Now, iColNow, iAcolNow); doUpdate = true; } } if (doUpdate) { tmpSM.col1 = colNow; tmpSM.acol1 = acolNow; tmpSM.iCol1 = iColNow; tmpSM.iAcol1 = iAcolNow; tmpSM.lambdaRef = lambdaRefNow; } } // Update info on dLambda for affected particles and colour lines. for (int iPair = 0; iPair < nPair; ++iPair) { InfoGluonMove& tmpSM = infoGM[iPair]; int iMod = -1; for (int i = 0; i < 3; ++i) if (tmpSM.col2 == col2Mod[i]) iMod = i; if (iMod > -1) tmpSM.iCol2 = iCol2Mod[iMod]; if (tmpSM.i1 == i1Sel || tmpSM.i1 == iCol1Sel || tmpSM.i1 == iAcol1Sel || tmpSM.i1 == iCol2Sel || tmpSM.i1 == iAcol2Sel || iMod > -1) tmpSM.dLambda = (tmpSM.iCol2 == tmpSM.i1 || tmpSM.iAcol2 == tmpSM.i1 || tmpSM.iCol1 == tmpSM.iAcol1) ? 2e4 : lambda123Move( tmpSM.i1, tmpSM.iCol2, tmpSM.iAcol2) - tmpSM.lambdaRef; } // End of loop over gluon shifting. } // Done if no flip. if (flipMode == 0) return true; // Array with colour lines, and where each line begins and ends. vector iTmpFlip, iVecFlip, iBegFlip, iEndFlip; // Variables for minimum search. int i1c, i1a, i2c, i2a, i1cMin, i1aMin, i2cMin, i2aMin, iSMin; double dLambdaFlip, dLambdaFlipMin; vector flipMin; // Grab all colour ends. for (int i = oldSize; i < event.size(); ++i) if (event[i].isFinal() && event[i].col() > 0 && event[i].acol() == 0) { iTmpFlip.clear(); iTmpFlip.push_back( i); // Step through colour neighbours to catch system. iNow = i; acolM = acolMap.find( event[iNow].col() ); bool foundEnd = false; while (acolM != acolMap.end()) { iNow = acolM->second; iTmpFlip.push_back( iNow); if (event[iNow].col() == 0) { foundEnd = true; break; } acolM = acolMap.find( event[iNow].col() ); } // Store acceptable system, optionally including junction legs. if (foundEnd || flipMode == 2 || flipMode == 4) { iBegFlip.push_back( iVecFlip.size()); for (int j = 0; j < int(iTmpFlip.size()); ++j) iVecFlip.push_back( iTmpFlip[j]); iEndFlip.push_back( iVecFlip.size()); } } // Optionally search for antijunction legs: grab all anticolour ends. if (flipMode == 2 || flipMode == 4) for (int i = oldSize; i < event.size(); ++i) if (event[i].isFinal() && event[i].acol() > 0 && event[i].col() == 0) { iTmpFlip.clear(); iTmpFlip.push_back( i); // Step through anticolour neighbours to catch system. iNow = i; colM = colMap.find( event[iNow].acol() ); bool foundEnd = false; while (colM != colMap.end()) { iNow = colM->second; iTmpFlip.push_back( iNow); if (event[iNow].acol() == 0) { foundEnd = true; break; } colM = colMap.find( event[iNow].acol() ); } // Store acceptable system, but do not doublecount q - (n g) - qbar. if (!foundEnd) { iBegFlip.push_back( iVecFlip.size()); for (int j = 0; j < int(iTmpFlip.size()); ++j) iVecFlip.push_back( iTmpFlip[j]); iEndFlip.push_back( iVecFlip.size()); } } // Loop through all system pairs. int nSysFlip = iBegFlip.size(); for (int iSys1 = 0; iSys1 < nSysFlip - 1; ++iSys1) if (iBegFlip[iSys1] >= 0) for (int iSys2 = iSys1 + 1; iSys2 < nSysFlip; ++iSys2) if (iBegFlip[iSys2] >= 0) { i1cMin = 0; i1aMin = 0; i2cMin = 0; i2aMin = 0; dLambdaFlipMin = 1e4; // Loop through all possible flip locations for a pair. for (int j1 = iBegFlip[iSys1]; j1 < iEndFlip[iSys1] - 1; ++j1) for (int j2 = iBegFlip[iSys2]; j2 < iEndFlip[iSys2] - 1; ++j2) { i1c = iVecFlip[j1]; i1a = iVecFlip[j1 + 1]; i2c = iVecFlip[j2]; i2a = iVecFlip[j2 + 1]; dLambdaFlip = lambda12Move( i1c, i2a) + lambda12Move( i2c, i1a) - lambda12Move( i1c, i1a) - lambda12Move( i2c, i2a); if (dLambdaFlip < dLambdaFlipMin) { i1cMin = i1c; i1aMin = i1a; i2cMin = i2c; i2aMin = i2a; dLambdaFlipMin = dLambdaFlip; } } // Store possible flips if low enough dLambdaMin. if (dLambdaFlipMin < -dLambdaCut) flipMin.push_back( InfoGluonMove( iSys1, iSys2, i1cMin, i1aMin, i2cMin, i2aMin, dLambdaFlipMin) ); } int flipSize = flipMin.size(); // Search for lowest possible flip among unused systems. for (int iFlip = 0; iFlip < min( nSysFlip / 2, flipSize); ++iFlip) { iSMin = -1; dLambdaFlipMin = 1e4; for (int iSys12 = 0; iSys12 < flipSize; ++iSys12) if (flipMin[iSys12].i1 >= 0 && flipMin[iSys12].dLambda < dLambdaFlipMin) { iSMin = iSys12; dLambdaFlipMin = flipMin[iSys12].dLambda; } // Do flip. Mark flipped systems. if (iSMin >= 0) { InfoGluonMove& flipNow = flipMin[iSMin]; int iS1 = flipNow.i1; int iS2 = flipNow.i2; event[ flipNow.iAcol1 ].acol( event[flipNow.iCol2].col() ); event[ flipNow.iAcol2 ].acol( event[flipNow.iCol1].col() ); for (int iSys12 = 0; iSys12 < flipSize; ++iSys12) if ( flipMin[iSys12].i1 == iS1 || flipMin[iSys12].i1 == iS2 || flipMin[iSys12].i2 == iS1 || flipMin[iSys12].i2 == iS2) flipMin[iSys12].i1 = -1; } else break; } // Done. return true; } //-------------------------------------------------------------------------- // Common code for the SK I and SK II models for WW/ZZ/WZ systems. bool ColourReconnection::reconnectTypeCommon( Event& event, int ) { // Make storage containers. Check that at least two parton systems. vector > dips; int iBosons[2]; Vec4 decays[2]; if (partonSystemsPtr->sizeSys() < 2) { infoPtr->errorMsg("Error in ColourReconnection::reconnectTypeCommon: " "expect at least two parton systems"); return false; } // Find the dipoles connected to their respective resonance decays. for (int i = 0; i < 2; ++i) { dips.push_back(vector()); int iSys = partonSystemsPtr->sizeSys() - i - 1; for (int j = 0; j < partonSystemsPtr->sizeOut(iSys); ++j) { int iPar = partonSystemsPtr->getOut(iSys, j); // Find decayed boson. Only need to do it once. if (j == 0) { int iMot = event[iPar].mother1(); while (iMot != 0 && event[iMot].idAbs() != 23 && event[iMot].idAbs() != 24) iMot = event[iMot].mother1(); if (iMot == 0) { infoPtr->errorMsg("Error in ColourReconnection::reconnectTypeCommon:" " Not a resonance decay of a W/Z"); return false; } iBosons[i] = iMot; } // Pick up all dipoles by starting from colour end. int col = event[iPar].col(); if (col == 0) continue; for (int k = 0; k < partonSystemsPtr->sizeOut(iSys); ++k) { int iAntiPar = partonSystemsPtr->getOut(iSys, k); if (event[iAntiPar].acol() == col) { dips.back().push_back(ColourDipole(col, iPar, iAntiPar)); break; } } } } // Boost system to W+W- rest frame. Vec4 boost = event[iBosons[0]].p() + event[iBosons[1]].p(); for (int i = 1; i < event.size(); ++i) event[i].bstback(boost); // Find the time and position of the decay vertices. for (int i = 0; i < 2; ++i) { int iBoson = iBosons[i]; double mBoson = particleDataPtr->m0(event[iBoson].idAbs()); double gammaBoson = particleDataPtr->mWidth(event[iBoson].idAbs()); double mReal = event[iBoson].mCalc(); decays[i][0] = -HBAR * log(rndmPtr->flat()) * event[iBoson].e() / sqrt(pow2(pow2(mBoson) - pow2(mReal)) + pow2(gammaBoson * pow2(mReal) / mBoson)); for (int j = 1; j < 4; ++j) decays[i][j] = event[iBoson].p()[j]/ event[iBoson].e() * decays[i][0]; } // Find the possible reconnections, depeding on choice of model. map > reconnections; if (reconnectMode == 3) reconnections = reconnectTypeI(event, dips, decays); else if (reconnectMode == 4) reconnections = reconnectTypeII(event, dips, decays); // Carry out the reconnections. vector used1(dips[0].size(), false), used2(dips[1].size(), false); for (map >::iterator it=reconnections.begin(); it != reconnections.end(); ++it) { // Check if any of the dipoles already reconnected. if (used1[it->second.first] || used2[it->second.second]) continue; // Do the reconnection. int iAcol1 = dips[0][it->second.first].iAcol; int iAcol2 = dips[1][it->second.second].iAcol; event.copy(iAcol1, 79); event.back().acol(event[iAcol2].acol()); event.copy(iAcol2, 79); event.back().acol(event[iAcol1].acol()); // Mark dipoles as used. used1[it->second.first] = true; used2[it->second.second] = true; // If only a single reconnection is wanted, break the loop. if (singleReconOnly) break; } // Boost system back to origianl rest frame. for (int i = 1; i < event.size(); ++i) event[i].bst(boost); // Done. return true; } //-------------------------------------------------------------------------- // The SK I model for WW/ZZ/WZ systems. map > ColourReconnection::reconnectTypeI( Event& event, vector > &dips, Vec4 decays[2]) { // Make storage containers. multimap > reconnections; // Velocity stored as: vx vy, vz , gamma (uses Vec4 for easy storage). vector > vel; // string direction. vector > dir; // Calculate velocities. for (int i = 0; i < 2; ++i) { vel.push_back(vector(dips[i].size(),Vec4())); dir.push_back(vector(dips[i].size(),Vec4())); for (int j = 0; j < int(dips[i].size()); ++j) { // Calculate sum of momenta. double pSumCol = event[dips[i][j].iCol].pAbs(); double pSumAcol = event[dips[i][j].iAcol].pAbs(); // Velocities and directions. for (int k = 1; k < 4; ++k) { vel[i][j][k] = 0.5 *(event[dips[i][j].iCol].p()[k] / pSumCol + event[dips[i][j].iAcol].p()[k] / pSumAcol); dir[i][j][k] = event[dips[i][j].iCol].p()[k] / pSumCol - event[dips[i][j].iAcol].p()[k] / pSumAcol; } vel[i][j][0] = 1/sqrt(1 - vel[i][j].pAbs2()); dir[i][j] /= dir[i][j].pAbs(); } } // Select random point. int nPoints = 100; double sumW = 0; for (int i = 0; i < nPoints; ++i) { Vec4 x; double r = sqrt(- log(rndmPtr->flat())); double phi = 2 * M_PI * rndmPtr->flat(); x[1] = blowR * rHadron * r * cos(phi); x[2] = blowR * rHadron * r * sin(phi); r = sqrt(- log(rndmPtr->flat())); phi = 2 * M_PI * rndmPtr->flat(); x[3] = blowR * rHadron * r * cos(phi); x[0] = max(decays[0][0], decays[1][0]) + blowT * sqrt(0.5) * tfrag * r * abs(sin(phi)); if (x.m2Calc() < 0) continue; // Find weight of trial point. double weightTrial = exp(-x.pAbs2() / pow2(blowR*rHadron)) * exp( -2. * pow2(x[0] - max(decays[0][0],decays[1][0])) / pow2(blowT * tfrag) ); // Find max weight and their indices. double maxWeights[2] = {0,0}; // {1E-10,1E-10}; int maxIndices[2] = {-1,-1}; // Loop over W decays. for (int j = 0;j < 2;++j) { // Calculate the difference between decay point and random point. Vec4 xd = x - decays[j]; // Loop over strings. for (int k = 0; k < int(dips[j].size()); ++k) { // Boost to rest frame of string. double gam = vel[j][k][0]; double dotProd = dot3(xd,vel[j][k]); double xBoost = gam * (gam * dotProd/ (1. + gam) - xd[0]); Vec4 xb = xd + xBoost * vel[j][k]; xb[0] = gam * (xd[0] - dotProd); // Check that position is reachable. if (xb.m2Calc() < 0) continue; // Find and store largest weight. double weight = exp( -(xb.pAbs2() - pow2(dot3(xb,dir[j][k]))) / (2 * pow2(rHadron)) ) * exp( -(pow2(xb[0]) - pow2(dot3(xb, dir[j][k]))) / pow2(tfrag) ); if (weight > maxWeights[j]) { maxWeights[j] = weight; maxIndices[j] = k; } } } // Store weight. if (maxIndices[0] != -1 && maxIndices[1] != -1) { // Check if the new configuration lowers the lambda measure. if (lowerLambdaOnly) { double oldLambda = stringLength.getStringLength(event, dips[0][maxIndices[0]].iCol, dips[0][maxIndices[0]].iAcol) + stringLength.getStringLength(event, dips[1][maxIndices[1]].iCol, dips[1][maxIndices[1]].iAcol); double newLambda = stringLength.getStringLength(event, dips[0][maxIndices[0]].iCol, dips[1][maxIndices[1]].iAcol) + stringLength.getStringLength(event, dips[1][maxIndices[1]].iCol, dips[0][maxIndices[0]].iAcol); if (oldLambda < newLambda) continue; } // Save weights. double weight = maxWeights[0] * maxWeights[1] / weightTrial; reconnections.insert(make_pair(weight, make_pair(maxIndices[0], maxIndices[1]))); sumW += weight; } } // check whether to do any reconnections. map > singleRecon; double result = pow3(blowR) * blowT * sumW/nPoints; // No reconections. if (1 - exp(-kI * result) < rndmPtr->flat()) { singleRecon.clear(); return singleRecon; // Find reconnection. } else { double rSum = sumW * rndmPtr->flat(); for (map >::iterator it = reconnections.begin(); it != reconnections.end(); ++it) { rSum -= it->first; if (rSum < 0.) { singleRecon.insert(make_pair(1., it->second)); return singleRecon; } } // If no solution was found (shoult not happen) do no reconnections. singleRecon.clear(); return singleRecon; } } //-------------------------------------------------------------------------- // The SK II model for WW/ZZ/WZ systems. map > ColourReconnection::reconnectTypeII( Event& event, vector > &dips, Vec4 decays[2]) { // Make storage containers. map > reconnections; // Find dipole velocities. for (int i = 0;i < int(dips[0].size()); ++i) { for (int j = 0;j < int(dips[1].size()); ++j) { Vec4 v1,v2,u1,u2; v1 = event[dips[0][i].iCol].p() / event[dips[0][i].iCol].e(); v2 = event[dips[0][i].iAcol].p() / event[dips[0][i].iAcol].e(); u1 = event[dips[1][j].iCol].p() / event[dips[1][j].iCol].e(); u2 = event[dips[1][j].iAcol].p() / event[dips[1][j].iAcol].e(); // Begin setup to solve system of equations. vector > matUpper, matLower; for (int k = 0; k < 3; ++k) { matUpper.push_back(vector(3,0)); matLower.push_back(vector(3,0)); } // Insert in matrix, beware of different convention in index. for (int k = 0;k < 3; ++k) { matUpper[0][k] = matLower[0][k] = v2[k+1]-v1[k+1]; matUpper[1][k] = matLower[1][k] = -(u2[k+1]-u1[k+1]); matUpper[2][k] = decays[0][k+1] - decays[1][k+1] - decays[0][0] * v1[k+1] + decays[1][0] * u1[k+1]; matLower[2][k] = v1[k+1] - u1[k+1]; } double t = -determinant3(matUpper) / determinant3(matLower); // Find alpha and beta of string crossing. double s11 = matUpper[0][0]*(t - decays[0][0]); double s12 = matUpper[1][0]*(t - decays[1][0]); double s13 = matUpper[2][0] + matLower[2][0]*t; double s21 = matUpper[0][1]*(t - decays[0][0]); double s22 = matUpper[1][1]*(t - decays[1][0]); double s23 = matUpper[2][1] + matLower[2][1]*t; double den = s11*s22 - s12*s21; double alpha = (s12*s23 - s22*s13)/den; double beta = (s21*s13 - s11*s23)/den; // Check if solution is physical. if (alpha < 0 || alpha > 1) continue; if (beta < 0 || beta > 1) continue; if (t < max(decays[0][0], decays[1][0])) continue; // Find the crossing points. Vec4 x1,x2; x1 = decays[0] + (alpha * v2 + (1 - alpha) * v1) * (t - decays[0][0]); x2 = decays[1] + (beta * u2 + (1 - beta) * u1) * (t - decays[1][0]); x1[0] = x2[0] = t; // Check that both points are the same. if (dot3(x1-x2,x1-x2) > 1E-4 * (x1.pAbs2() + x2.pAbs2())) continue; // Find string eigentimes. double tauPlus = (x1 - decays[0]).mCalc(); double tauMinus = (x2 - decays[1]).mCalc(); // Check if strings already decayed. if (rndmPtr->flat() > exp( -(pow2(tauPlus) + pow2(tauMinus)) / pow2(tfrag))) continue; // Optionally only accept if reconnection reduces string length. if (lowerLambdaOnly) { double oldLambda = stringLength.getStringLength(event, dips[0][i].iCol, dips[0][i].iAcol) + stringLength.getStringLength(event, dips[1][j].iCol, dips[1][j].iAcol); double newLambda = stringLength.getStringLength(event, dips[0][i].iCol, dips[1][j].iAcol) + stringLength.getStringLength(event, dips[1][j].iCol, dips[0][i].iAcol); if (oldLambda < newLambda) continue; } // Store dipole pair. reconnections.insert(make_pair(t,make_pair(i,j))); } } return reconnections; } //-------------------------------------------------------------------------- // Calculate the determinant of a 3 * 3 matrix. double ColourReconnection::determinant3(vector >& vec) { return vec[0][0]*vec[1][1]*vec[2][2] + vec[0][1]*vec[1][2]*vec[2][0] + vec[0][2]*vec[1][0]*vec[2][1] - vec[0][0]*vec[2][1]*vec[1][2] - vec[0][1]*vec[1][0]*vec[2][2] - vec[0][2]*vec[1][1]*vec[2][0]; } //========================================================================== } // end namespace Pythia8