EIC code displayed by LXR master/​include/​Pythia8Plugins/​ColourReconnectionHooks.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-18 10:06:36

 // ColourReconnectionHooks.h is a part of the PYTHIA event generator.
 // Copyright (C) 2024 Torbjorn Sjostrand.
 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
 // Please respect the MCnet Guidelines, see GUIDELINES for details.
 
 // This file contains two UserHooks that, along with the internal models,
 // implement all the models used for the top mass study in
 // S. Argyropoulos and T. Sjostrand,
 // arXiv:1407.6653 [hep-ph] (LU TP 14-23, DESY 14-134, MCnet-14-15)
 
 // MBReconUserHooks: can be used for all kinds of events, not only top ones.
 // TopReconUserHooks: models intended specifically for top decay products,
 // whereas the underlying event is handled by the default model.
 
 // Warning: some small modifications have been made when collecting
 // the models, but nothing intended to change the behaviour.
 // Note: the move model is also available with ColourReconnection:mode = 2,
 // while the ColourReconnection:mode = 1 model has not been used here.
 // Note: the new models tend to be slower than the default CR scenario,
 // since they have to probe many more reconnection possibilities.
 
 #ifndef Pythia8_ColourReconnectionHooks_H
 #define Pythia8_ColourReconnectionHooks_H
 
 // Includes
 #include "Pythia8/Pythia.h"
 namespace Pythia8 {
 
 //==========================================================================
 
 // Class for colour reconnection models of general validity.
 
 class MBReconUserHooks : public UserHooks {
 
 public:
 
   // Constructor and destructor.
   // mode = 0: no reconnection (dummy option, does nothing);
   //      = 1: swap gluons to minimize lambda.
   //      = 2: move gluons to minimize lambda.
   // flip = 0: no flip between quark-antiquark ends.
   //      = 1: flip between quark-antiquark ends, excluding junction systems.
   //      = 2: flip between quark-antiquark ends, including junction systems.
   // dLamCut: smallest -delta-lambda value for which to swap/mode (positive).
   // fracGluon: the fraction of gluons that will be studied for reconnection.
   // m2Ref   : squared reference mass scale for lambda measure calculation.
   MBReconUserHooks(int modeIn = 0, int flipIn = 0, double dLamCutIn = 0.,
     double fracGluonIn = 1.) : mode(modeIn), flip(flipIn), dLamCut(dLamCutIn),
     fracGluon(fracGluonIn) { m2Ref = 1.; dLamCut = max(0., dLamCut); }
   ~MBReconUserHooks() {}
 
   // Allow colour reconnection after resonance decays (early or late)...
   virtual bool canReconnectResonanceSystems() {return true;}
 
   // ...which gives access to the event, for modification.
   virtual bool doReconnectResonanceSystems( int, Event& event) {
 
     // Return without action for relevant mode numbers.
     if (mode <= 0 || mode > 2) return true;
 
     // Double diffraction not yet implemented, so return without action.
     // (But works for internal move implementation.)
     if (infoPtr->isDiffractiveA() && infoPtr->isDiffractiveB())
       return true;
 
     // Initial setup: relevant gluons and coloured partons.
     if (!setupConfig( event)) return false;
 
     // Done if not enough gluons.
     if ( (mode == 1 && nGlu < 2) || (mode == 2 && nGlu < 1) ) return true;
 
     // Colour reconnect. Return if failed.
     bool hasRec = (mode == 1) ? doReconnectSwap( event)
                               : doReconnectMove( event);
     if (!hasRec) return false;
 
     // Colour flip afterburner.
     if (flip > 0) return doReconnectFlip( event);
     return true;
 
   }
 
   // Return number of reconnections for current event.
   //int numberReconnections() {return nRec;}
   //double dLambdaReconnections() {return -dLamTot;}
 
 private:
 
   // Mode. Number of reconnections. lambda measure reference scale.
   int    mode, flip, nRec, nGlu, nAllCol, nCol;
   double dLamCut, fracGluon, m2Ref, dLamTot;
 
   // Array of (indices of) final gluons.
   vector<int> iGlu;
 
   // Array of (indices of) all final coloured particles.
   vector<int> iToAllCol, iAllCol;
 
   // Maps telling where all colours and anticolours are stored.
   map<int, int> colMap, acolMap;
 
   // Array of all lambda distances between coloured partons.
   vector<double> lambdaij;
 
   // Function to return lambda value from array.
   double lambda12( int i, int j) {
     int iAC = iToAllCol[i]; int jAC = iToAllCol[j];
     return lambdaij[nAllCol * min( iAC, jAC) + max( iAC, jAC)];
   }
 
   // Function to return lambda(i,j) + lambda(i,k) - lambda(j,k).
   double lambda123( int i, int j, int k) {
     int iAC = iToAllCol[i]; int jAC = iToAllCol[j]; int kAC = iToAllCol[k];
     return lambdaij[nAllCol * min( iAC, jAC) + max( iAC, jAC)]
          + lambdaij[nAllCol * min( iAC, kAC) + max( iAC, kAC)]
          - lambdaij[nAllCol * min( jAC, kAC) + max( jAC, kAC)];
   }
 
   // Small nested class for the effect of a potential gluon swap or move.
   class InfoSwapMove{
     public:
     InfoSwapMove(int i1in = 0, int i2in = 0) : i1(i1in), i2(i2in) {}
     InfoSwapMove(int i1in, int i2in, int iCol1in, int iAcol1in, int iCol2in,
       int iAcol2in, int dLamIn) : i1(i1in), i2(i2in), iCol1(iCol1in),
       iAcol1(iAcol1in), iCol2(iCol2in), iAcol2(iAcol2in), dLam(dLamIn) {}
     ~InfoSwapMove() {}
     int i1, i2, col1, acol1, iCol1, iAcol1, col2, acol2, iCol2, iAcol2;
     double lamNow, dLam;
   };
 
   // Vector of (1) gluon-pair swap or (2) gluon move properties.
   vector<InfoSwapMove> infoSM;
 
   //----------------------------------------------------------------------
 
   // Initial setup: relevant gluons and coloured partons.
 
   inline bool setupConfig(Event& event) {
 
     // Reset arrays to prepare for the new event analysis.
     iGlu.clear();
     iToAllCol.resize( event.size() );
     iAllCol.clear();
     colMap.clear();
     acolMap.clear();
     infoSM.clear();
     nRec = 0;
     dLamTot = 0.;
 
     // Loop over all final particles. Store (fraction of) gluons.
     for (int i = 3; i < event.size(); ++i) if (event[i].isFinal()) {
       if (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) {
         iToAllCol[i] = iAllCol.size();
         iAllCol.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 (map<int, int>::iterator 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 (map<int, int>::iterator 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.
     nGlu = iGlu.size();
     nCol = colMap.size();
     if (int(acolMap.size()) != nCol) {
       loggerPtr->ERROR_MSG("map sizes do not match");
       return false;
     }
     map<int, int>::iterator colM = colMap.begin();
     map<int, int>::iterator acolM = acolMap.begin();
     for (int iCol = 0; iCol < nCol; ++iCol) {
       if (colM->first != acolM->first) {
         loggerPtr->ERROR_MSG("map elements do not match");
         return false;
       }
       ++colM;
       ++acolM;
     }
 
     // Calculate and tabulate lambda between any pair of coloured partons.
     nAllCol = iAllCol.size();
     lambdaij.resize( pow2(nAllCol) );
     int i, j;
     for (int iAC = 0; iAC < nAllCol - 1; ++iAC) {
       i = iAllCol[iAC];
       for (int jAC = iAC + 1; jAC < nAllCol; ++jAC) {
         j = iAllCol[jAC];
         lambdaij[nAllCol * iAC + jAC]
           = log1p(m2( event[i], event[j]) / m2Ref);
       }
     }
 
     // Done.
     return true;
 
   }
 
   //----------------------------------------------------------------------
 
   // Swap gluons by lambda measure.
 
   inline bool doReconnectSwap(Event& event) {
 
     // Set up initial possible gluon swap pairs with lambda gains/losses.
     for (int iG1 = 0; iG1 < nGlu - 1; ++iG1) {
       int i1 = iGlu[iG1];
       for (int iG2 = iG1 + 1; iG2 < nGlu; ++iG2) {
         int i2 = iGlu[iG2];
         InfoSwapMove tmpSM( i1, i2);
         calcLamSwap( tmpSM, event);
         infoSM.push_back( tmpSM);
       }
     }
     int nPair = infoSM.size();
 
     // Keep on looping over swaps until no further negative dLambda.
     for ( int iSwap = 0; iSwap < nGlu; ++iSwap) {
       int    iPairMin = -1;
       double dLamMin  = 1e4;
 
       // Find lowest dLambda.
       for (int iPair = 0; iPair < nPair; ++iPair)
       if (infoSM[iPair].dLam < dLamMin) {
         iPairMin = iPair;
         dLamMin  = infoSM[iPair].dLam;
       }
 
       // Break if no shift below upper limit found.
       if (dLamMin > -dLamCut) break;
       ++nRec;
       dLamTot += dLamMin;
       int i1min = infoSM[iPairMin].i1;
       int i2min = infoSM[iPairMin].i2;
 
       // Swap the colours in the event record.
       int col1  = event[i1min].col();
       int acol1 = event[i1min].acol();
       int col2  = event[i2min].col();
       int acol2 = event[i2min].acol();
       event[i1min].cols( col2, acol2);
       event[i2min].cols( col1, acol1);
 
       // Swap the indices in the colour maps.
       colMap[col1]   = i2min;
       acolMap[acol1] = i2min;
       colMap[col2]   = i1min;
       acolMap[acol2] = i1min;
 
       // Remove already swapped pair from further consideration.
       infoSM[iPairMin] = infoSM.back();
       infoSM.pop_back();
       --nPair;
 
       // Update all pairs that have been affected.
       for (int iPair = 0; iPair < nPair; ++iPair) {
         InfoSwapMove& tmpSM = infoSM[iPair];
         if ( tmpSM.i1     == i1min || tmpSM.i1     == i2min
           || tmpSM.i2     == i1min || tmpSM.i2     == i2min
           || tmpSM.iCol1  == i1min || tmpSM.iCol1  == i2min
           || tmpSM.iAcol1 == i1min || tmpSM.iAcol1 == i2min
           || tmpSM.iCol2  == i1min || tmpSM.iCol2  == i2min
           || tmpSM.iAcol2 == i1min || tmpSM.iAcol2 == i2min)
           calcLamSwap( tmpSM, event);
       }
     }
 
     // Done.
     return true;
 
   }
 
   //----------------------------------------------------------------------
 
   // Calculate pair swap properties.
 
   inline void calcLamSwap( InfoSwapMove& tmpSM, Event& event) {
 
     // Colour line tracing to neighbours.
     tmpSM.col1   = event[tmpSM.i1].col();
     tmpSM.acol1  = event[tmpSM.i1].acol();
     tmpSM.iCol1  = acolMap[tmpSM.col1];
     tmpSM.iAcol1 = colMap[tmpSM.acol1];
     tmpSM.col2   = event[tmpSM.i2].col();
     tmpSM.acol2  = event[tmpSM.i2].acol();
     tmpSM.iCol2  = acolMap[tmpSM.col2];
     tmpSM.iAcol2 = colMap[tmpSM.acol2];
 
     // Lambda swap properties.
     double lam1c = lambda12( tmpSM.i1, tmpSM.iCol1);
     double lam1a = lambda12( tmpSM.i1, tmpSM.iAcol1);
     double lam2c = lambda12( tmpSM.i2, tmpSM.iCol2);
     double lam2a = lambda12( tmpSM.i2, tmpSM.iAcol2);
     double lam3c = lambda12( tmpSM.i1, tmpSM.iCol2);
     double lam3a = lambda12( tmpSM.i1, tmpSM.iAcol2);
     double lam4c = lambda12( tmpSM.i2, tmpSM.iCol1);
     double lam4a = lambda12( tmpSM.i2, tmpSM.iAcol1);
     if (tmpSM.col1 == tmpSM.acol2 && tmpSM.acol1 == tmpSM.col2)
        tmpSM.dLam = 2e4;
     else if (tmpSM.col1 == tmpSM.acol2)
        tmpSM.dLam = (lam3c + lam4a) - (lam1a + lam2c);
     else if (tmpSM.acol1 == tmpSM.col2)
        tmpSM.dLam = (lam3a + lam4c) - (lam1c + lam2a);
     else tmpSM.dLam = (lam3c + lam3a + lam4c + lam4a)
                    - (lam1c + lam1a + lam2c + lam2a);
 
   // Done.
   }
 
   //----------------------------------------------------------------------
 
   // Move gluons by lambda measure.
 
   inline bool doReconnectMove(Event& event) {
 
     // Temporary variables.
     int    iNow, colNow, acolNow, iColNow, iAcolNow, col2Now;
     double lamNow;
 
     // Set up initial possible gluon moves with lambda gains/losses.
     for (int iG = 0; iG < nGlu; ++iG) {
 
       // Gluon and its 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.
       lamNow   = lambda123( iNow, iColNow, iAcolNow);
 
       // Loop over all colour lines where gluon could be inserted.
       for (map<int, int>::iterator colM = colMap.begin();
       colM != colMap.end(); ++colM) {
         col2Now = colM->first;
 
         // New container for gluon and colour line information.
         InfoSwapMove tmpSM( iNow);
         tmpSM.col1   = colNow;
         tmpSM.acol1  = acolNow;
         tmpSM.iCol1  = iColNow;
         tmpSM.iAcol1 = iAcolNow;
         tmpSM.lamNow = lamNow;
         tmpSM.col2   = col2Now;
         tmpSM.iCol2  = colMap[col2Now];
         tmpSM.iAcol2 = acolMap[col2Now];
 
         // Addition to Lambda if gluon inserted on line.
         tmpSM.dLam = (tmpSM.iCol2 == tmpSM.i1 || tmpSM.iAcol2 == tmpSM.i1
           || tmpSM.iCol1 == tmpSM.iAcol1) ? 2e4
           : lambda123( tmpSM.i1, tmpSM.iCol2, tmpSM.iAcol2) - tmpSM.lamNow;
         infoSM.push_back( tmpSM);
       }
     }
     int nPair = infoSM.size();
 
     // Keep on looping over moves until no further negative dLambda.
     for ( int iMove = 0; iMove < nGlu; ++iMove) {
       int    iPairMin = -1;
       double dLamMin  = 1e4;
 
       // Find lowest dLambda.
       for (int iPair = 0; iPair < nPair; ++iPair)
       if (infoSM[iPair].dLam < dLamMin) {
         iPairMin = iPair;
         dLamMin  = infoSM[iPair].dLam;
       }
 
       // Break if no shift below upper limit found.
       if (dLamMin > -dLamCut) break;
       ++nRec;
       dLamTot += dLamMin;
 
       // Partons and colours involved in move.
       InfoSwapMove& selSM = infoSM[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) {
         InfoSwapMove& tmpSM = infoSM[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];
             lamNow   = lambda123( i1Now, iColNow, iAcolNow);
             doUpdate = true;
           }
         }
         if (doUpdate) {
           tmpSM.col1   = colNow;
           tmpSM.acol1  = acolNow;
           tmpSM.iCol1  = iColNow;
           tmpSM.iAcol1 = iAcolNow;
           tmpSM.lamNow = lamNow;
         }
       }
 
       // Update info on dLambda for affected particles and colour lines.
       for (int iPair = 0; iPair < nPair; ++iPair) {
         InfoSwapMove& tmpSM = infoSM[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.dLam = (tmpSM.iCol2 == tmpSM.i1 || tmpSM.iAcol2 == tmpSM.i1
             || tmpSM.iCol1 == tmpSM.iAcol1) ? 2e4
             : lambda123( tmpSM.i1, tmpSM.iCol2, tmpSM.iAcol2) - tmpSM.lamNow;
       }
 
     // End of loop over gluon shifting.
     }
 
     // Done.
     return true;
 
   }
 
   //----------------------------------------------------------------------
 
   // Flip colour chains by lambda measure.
 
   inline bool doReconnectFlip(Event& event) {
 
     // Array with colour lines, and where each line begins and ends.
     vector<int> iTmp, iVec, iBeg, iEnd;
 
     // Grab all colour ends.
     for (int i = 3; i < event.size(); ++i)
     if (event[i].isFinal() && event[i].col() > 0 && event[i].acol() == 0) {
       iTmp.clear();
       iTmp.push_back( i);
 
       // Step through colour neighbours to catch system.
       int iNow = i;
       map<int, int>::iterator acolM = acolMap.find( event[iNow].col() );
       bool foundEnd = false;
       while (acolM != acolMap.end()) {
         iNow = acolM->second;
         iTmp.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 || flip == 2) {
         iBeg.push_back( iVec.size());
         for (int j = 0; j < int(iTmp.size()); ++j) iVec.push_back( iTmp[j]);
         iEnd.push_back( iVec.size());
       }
     }
 
     // Optionally search for antijunction legs: grab all anticolour ends.
     if (flip == 2) for (int i = 3; i < event.size(); ++i)
     if (event[i].isFinal() && event[i].acol() > 0 && event[i].col() == 0) {
       iTmp.clear();
       iTmp.push_back( i);
 
       // Step through anticolour neighbours to catch system.
       int iNow = i;
       map<int, int>::iterator colM = colMap.find( event[iNow].acol() );
       bool foundEnd = false;
       while (colM != colMap.end()) {
         iNow = colM->second;
         iTmp.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) {
         iBeg.push_back( iVec.size());
         for (int j = 0; j < int(iTmp.size()); ++j) iVec.push_back( iTmp[j]);
         iEnd.push_back( iVec.size());
       }
     }
 
 
     // Variables for minimum search.
     int nSys = iBeg.size();
     int i1c, i1a, i2c, i2a, i1cMin, i1aMin, i2cMin, i2aMin, iSMin;
     double dLam, dLamMin;
     vector<InfoSwapMove> flipMin;
 
     // Loop through all system pairs.
     for (int iSys1 = 0; iSys1 < nSys - 1; ++iSys1) if (iBeg[iSys1] >= 0)
     for (int iSys2 = iSys1 + 1; iSys2 < nSys; ++iSys2) if (iBeg[iSys2] >= 0) {
       i1cMin = 0;
       i1aMin = 0;
       i2cMin = 0;
       i2aMin = 0;
       dLamMin = 1e4;
 
       // Loop through all possible flip locations for a pair.
       for (int j1 = iBeg[iSys1]; j1 < iEnd[iSys1] - 1; ++j1)
       for (int j2 = iBeg[iSys2]; j2 < iEnd[iSys2] - 1; ++j2) {
         i1c = iVec[j1];
         i1a = iVec[j1 + 1];
         i2c = iVec[j2];
         i2a = iVec[j2 + 1];
         dLam = lambda12( i1c, i2a) + lambda12( i2c, i1a)
              - lambda12( i1c, i1a) - lambda12( i2c, i2a);
         if (dLam < dLamMin) {
           i1cMin = i1c;
           i1aMin = i1a;
           i2cMin = i2c;
           i2aMin = i2a;
           dLamMin = dLam;
         }
       }
 
       // Store possible flips if low enough dLamMin.
       if (dLamMin < -dLamCut) flipMin.push_back( InfoSwapMove(
         iSys1, iSys2, i1cMin, i1aMin, i2cMin, i2aMin, dLamMin) );
     }
     int flipSize = flipMin.size();
 
     // Search for lowest possible flip among unused systems.
     for (int iFlip = 0; iFlip < min( nSys / 2, flipSize); ++iFlip) {
       iSMin = -1;
       dLamMin  = 1e4;
       for (int iSys12 = 0; iSys12 < flipSize; ++iSys12)
       if (flipMin[iSys12].i1 >= 0 && flipMin[iSys12].dLam < dLamMin) {
         iSMin   = iSys12;
         dLamMin = flipMin[iSys12].dLam;
       }
 
       // Do flip. Mark flipped systems.
       if (iSMin >= 0) {
         InfoSwapMove& 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() );
         ++nRec;
         dLamTot += flipNow.dLam;
         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;
 
   }
 
 };
 
 //==========================================================================
 
 
 // Class for colour reconnection models specifically aimed at top decays.
 
 class TopReconUserHooks : public UserHooks {
 
 public:
 
   // Constructor and destructor.
   // mode = 0: no reconnection of tops (dummy option, does nothing);
   //      = 1: reconnect with random background gluon;
   //      = 2: reconnect with nearest (smallest-mass) background gluon;
   //      = 3: reconnect with furthest (largest-mass) background gluon;
   //      = 4: reconnect with smallest (with sign) lambda measure shift;
   //      = 5: reconnect only if reduced lamda, and then to most reduction.
   // strength: fraction of top gluons that is to be colour reconnected.
   // nList: list first nList parton classifications.
   // pTolerance: acceptable total momentum error (in reconstruction check).
   // m2Ref: squared reference mass scale for lambda measure calculation.
   // Possible variants for the future: swap with nearest in angle, not mass,
   // and/or only allow a background gluon to swap colours once.
 
   TopReconUserHooks(int modeIn = 0, double strengthIn = 1.) : mode(modeIn),
     strength(strengthIn) { iList = 0; nList = 0; pTolerance = 0.01;
     m2Ref = 1.;}
   ~TopReconUserHooks() {}
 
   // Allow colour reconnection after resonance decays (early or late)...
   virtual bool canReconnectResonanceSystems() {return true;}
 
   // ...which gives access to the event, for modification.
   virtual bool doReconnectResonanceSystems( int, Event& event) {
 
     // Return without action for relevant mode numbers.
     if (mode <= 0 || mode > 5) return true;
 
     // Classify coloured final partons.
     classifyFinalPartons(event);
 
     // Check that classification worked as expected.
     if (!checkClassification(event)) return false;
 
     // List first few classifications, along with the event.
     if (iList++ < nList) {
       listClassification();
       event.list();
     }
 
     // Perform reconnection for t and tbar in random order.
     bool tqrkFirst = (rndmPtr->flat() < 0.5);
     doReconnect( tqrkFirst, event);
     doReconnect(!tqrkFirst, event);
 
     // Done.
     return true;
   }
 
   // Return number of reconnections for current event.
   //int numberReconnections() {return nRec;}
 
 private:
 
   // Mode. Counters for how many events to list. Allowed momentum error.
   int    mode, iList, nList, nRec;
   double strength, pTolerance, m2Ref;
 
   // Arrays of (indices of) final partons from different sources.
   // So far geared towards t -> b W decays only.
   vector<int> iBqrk, iWpos, iTqrk, iBbar, iWneg, iTbar, iRest;
 
   // Maps telling where all colours and anticolours are stored.
   map<int, int> colMap, acolMap;
 
   //----------------------------------------------------------------------
 
   // Classify all coloured partons at the end of showers by origin.
   // Note: for now only t -> b W is fully classified.
 
   inline bool classifyFinalPartons(Event& event) {
 
     // Reset arrays to prepare for the new event analysis.
     iBqrk.clear();
     iWpos.clear();
     iTqrk.clear();
     iBbar.clear();
     iWneg.clear();
     iTbar.clear();
     iRest.clear();
     colMap.clear();
     acolMap.clear();
     nRec = 0;
 
     // Loop over all final particles. Tag coloured ones.
     for (int i = 3; i < event.size(); ++i) if (event[i].isFinal()) {
       bool hasCol = (event[i].colType() != 0);
 
       // Set up to find where each parton comes from.
       bool fsrFromT = false;
       bool fromTqrk = false;
       bool fromBqrk = false;
       bool fromWpos = false;
       bool fromTbar = false;
       bool fromBbar = false;
       bool fromWneg = false;
 
       // Identify current particle.
       int iNow  = i;
       int idOld = 0;
       do {
         int idNow = event[iNow].id();
 
         // Exclude FSR of gluons/photons from the t quark proper.
         if (abs(idNow) == 6 && (idOld == 21 || idOld == 22)) fsrFromT = true;
 
         // Check if current particle matches any of the categories.
         else if (idNow ==   6) fromTqrk = true;
         else if (idNow ==   5) fromBqrk = true;
         else if (idNow ==  24) fromWpos = true;
         else if (idNow ==  -6) fromTbar = true;
         else if (idNow ==  -5) fromBbar = true;
         else if (idNow == -24) fromWneg = true;
 
         // Step up through the history to the very top.
         iNow  = event[iNow].mother1();
         idOld = idNow;
       } while (iNow > 2 && !fsrFromT);
 
       // Bookkeep where the parton comes from. Note that b quarks also
       // can appear in W decays, so order of checks is relevant.
       if      (fromTqrk && fromWpos && hasCol) iWpos.push_back(i);
       else if (fromTqrk && fromBqrk && hasCol) iBqrk.push_back(i);
       else if (fromTqrk)                       iTqrk.push_back(i);
       else if (fromTbar && fromWneg && hasCol) iWneg.push_back(i);
       else if (fromTbar && fromBbar && hasCol) iBbar.push_back(i);
       else if (fromTbar)                       iTbar.push_back(i);
       else if (hasCol)                         iRest.push_back(i);
 
       // Store location of all colour and anticolour indices.
       if (hasCol && (mode == 4 || mode == 5)) {
         if (event[i].col() > 0) colMap[event[i].col()] = i;
         if (event[i].acol() > 0) acolMap[event[i].acol()] = i;
       }
     }
 
     // So far method always returns true.
     return true;
 
   }
 
   //----------------------------------------------------------------------
 
   // Check that classification worked by summing up partons to t/tbar.
 
   inline bool checkClassification(Event& event) {
 
     // Find final copy of t and tbar quarks.
     int iTqrkLoc = 0;
     int iTbarLoc = 0;
     for (int i = 3; i < event.size(); ++i) {
       if(event[i].id() ==  6) iTqrkLoc = i;
       if(event[i].id() == -6) iTbarLoc = i;
     }
 
     // Four-momentum of t minus all its decay products.
     Vec4 tqrkDiff = event[iTqrkLoc].p();
     for (int i = 0; i < int(iBqrk.size()); ++i)
       tqrkDiff -= event[iBqrk[i]].p();
     for (int i = 0; i < int(iWpos.size()); ++i)
       tqrkDiff -= event[iWpos[i]].p();
     for (int i = 0; i < int(iTqrk.size()); ++i)
       tqrkDiff -= event[iTqrk[i]].p();
 
     // Four-momentum of tbar minus all its decay products.
     Vec4 tbarDiff = event[iTbarLoc].p();
     for (int i = 0; i < int(iBbar.size()); ++i)
       tbarDiff -= event[iBbar[i]].p();
     for (int i = 0; i < int(iWneg.size()); ++i)
       tbarDiff -= event[iWneg[i]].p();
     for (int i = 0; i < int(iTbar.size()); ++i)
       tbarDiff -= event[iTbar[i]].p();
 
     // Print difference vectors and event if sum deviation is too big.
     double totErr = abs(tqrkDiff.px()) + abs(tqrkDiff.py())
       + abs(tqrkDiff.pz()) + abs(tqrkDiff.e()) + abs(tbarDiff.px())
       + abs(tbarDiff.py()) + abs(tbarDiff.pz()) + abs(tqrkDiff.e());
     if (totErr > pTolerance) {
       loggerPtr->ERROR_MSG("Error in t/tbar daughter search");
       cout << "\n Error in t/tbar daughter search: \n t    difference "
            << tqrkDiff << " tbar difference "<< tbarDiff;
       listClassification();
       event.list();
     }
 
     // Done.
     return (totErr < pTolerance);
   }
 
   //----------------------------------------------------------------------
 
   // Print how final-state (mainly coloured) particles were classified.
 
   inline void listClassification() {
 
     cout << "\n Final-state coloured partons classified by source: ";
     cout << "\n From Bqrk:";
     for (int i = 0; i < int(iBqrk.size()); ++i) cout << "  " << iBqrk[i];
     cout << "\n From Wpos:";
     for (int i = 0; i < int(iWpos.size()); ++i) cout << "  " << iWpos[i];
     cout << "\n From Tqrk:";
     for (int i = 0; i < int(iTqrk.size()); ++i) cout << "  " << iTqrk[i];
     cout << "\n From Bbar:";
     for (int i = 0; i < int(iBbar.size()); ++i) cout << "  " << iBbar[i];
     cout << "\n From Wneg:";
     for (int i = 0; i < int(iWneg.size()); ++i) cout << "  " << iWneg[i];
     cout << "\n From Tbar:";
     for (int i = 0; i < int(iTbar.size()); ++i) cout << "  " << iTbar[i];
     cout << "\n From Rest:";
     for (int i = 0; i < int(iRest.size()); ++i) {
       cout << "  " << iRest[i];
       if (i%20 == 19 && i + 1 != int(iRest.size())) cout << "\n           ";
     }
     cout << endl;
   }
 
   //----------------------------------------------------------------------
 
   // Reconnect gluons either from t or from tbar quark.
 
   inline bool doReconnect(bool doTqrk, Event& event) {
 
     // Gather coloured decay products either of t or of tbar.
     vector<int> iTdec;
     if (doTqrk) {
       for (int i = 0; i < int(iBqrk.size()); ++i) iTdec.push_back(iBqrk[i]);
       for (int i = 0; i < int(iWpos.size()); ++i) iTdec.push_back(iWpos[i]);
     } else {
       for (int i = 0; i < int(iBbar.size()); ++i) iTdec.push_back(iBbar[i]);
       for (int i = 0; i < int(iWneg.size()); ++i) iTdec.push_back(iWneg[i]);
     }
 
     // Extract the gluons from the t quark decay.
     vector<int> iGT;
     for (int i = 0; i < int(iTdec.size()); ++i) {
       int colNow  = event[iTdec[i]].col();
       int acolNow = event[iTdec[i]].acol();
       if (colNow > 0 && acolNow > 0) iGT.push_back(iTdec[i]);
     }
     int nGT    = iGT.size();
 
     // Randomize their stored order.
     if (nGT > 1) for (int i = 0; i < nGT; ++i) {
       int j = min( int(nGT * rndmPtr->flat()), nGT - 1 );
       swap( iGT[i], iGT[j]);
     }
 
     // Also extract the rest of the gluons in the event.
     vector<int> iGR;
     for (int i = 0; i < int(iRest.size()); ++i) {
       int colNow = event[iRest[i]].col();
       int acolNow = event[iRest[i]].acol();
       if (colNow > 0 && acolNow > 0) iGR.push_back(iRest[i]);
     }
     int nGR    = iGR.size();
     int iR, colT, acolT, iColT, iAcolT, colR, acolR, iColR, iAcolR;
     double mTR2, mTR2now, dLam = 0.0, lamT, lamNow, lamRec;
 
     // Loop through all top gluons; study fraction given by strength.
     if (nGT > 0 && nGR > 0)
     for (int iT = 0; iT < nGT; ++iT) {
       if (strength < rndmPtr->flat()) continue;
 
       // Pick random gluon from rest of event.
       if (mode == 1) iR = min( int(nGR * rndmPtr->flat()), nGR - 1 );
 
       // Find gluon from rest with lowest or highest invariant mass.
       else if (mode < 4) {
         iR   = 0;
         mTR2 = m2( event[iGT[iT]], event[iGR[iR]]);
         for (int ii = 1; ii < nGR; ++ii) {
           mTR2now = m2( event[iGT[iT]], event[iGR[ii]]);
           if (mode == 2 && mTR2now < mTR2) {iR = ii; mTR2 = mTR2now;}
           if (mode == 3 && mTR2now > mTR2) {iR = ii; mTR2 = mTR2now;}
         }
 
       // Find gluon from rest with smallest lambda value shift.
       } else {
         iR     = -1;
         dLam   = 1e10;
         colT   = event[iGT[iT]].col();
         acolT  = event[iGT[iT]].acol();
         iColT  = acolMap[colT];
         iAcolT = colMap[acolT];
         lamT   = log1p(m2( event[iGT[iT]], event[iColT]) / m2Ref)
                + log1p(m2( event[iGT[iT]], event[iAcolT]) / m2Ref);
         for (int ii = 0; ii < nGR; ++ii) {
           colR   = event[iGR[ii]].col();
           acolR  = event[iGR[ii]].acol();
           iColR  = acolMap[colR];
           iAcolR = colMap[acolR];
           lamNow = lamT
                  + log1p(m2( event[iGR[ii]], event[iColR]) / m2Ref)
                  + log1p(m2( event[iGR[ii]], event[iAcolR]) / m2Ref);
           lamRec = log1p(m2( event[iGT[iT]], event[iColR]) / m2Ref)
                  + log1p(m2( event[iGT[iT]], event[iAcolR]) / m2Ref)
                  + log1p(m2( event[iGR[ii]], event[iColT]) / m2Ref)
                  + log1p(m2( event[iGR[ii]], event[iAcolT]) / m2Ref);
           if (lamRec - lamNow < dLam) {iR = ii; dLam = lamRec - lamNow;}
         }
         if (mode == 5 && dLam > 0.) continue;
       }
 
       // Swap top and rest gluon colour and anticolour.
       ++nRec;
       swapCols( iGT[iT], iGR[iR], event);
     }
 
     // Done.
     return true;
 
   }
 
   //----------------------------------------------------------------------
 
   // Swap colours and/or anticolours in the event listing.
 
   inline void swapCols( int i, int j, Event& event) {
 
     // Swap the colours in the event record.
     int coli  = event[i].col();
     int acoli = event[i].acol();
     int colj  = event[j].col();
     int acolj = event[j].acol();
     event[i].cols( colj, acolj);
     event[j].cols( coli, acoli);
 
     // Swap the indices in the colour maps.
     colMap[coli]   = j;
     acolMap[acoli] = j;
     colMap[colj]   = i;
     acolMap[acolj] = i;
 
   }
 
 };
 
 //==========================================================================
 
 } // end namespace Pythia8
 
 #endif // end Pythia8_ColourReconnectionHooks_H

This page was automatically generated by the 2.3.7 LXR engine. The LXR team