src/JunctionSplitting.cc

// JunctionSplitting.cc 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.

// Function definitions (not found in the header) for the
// JunctionSplitting class.

// Setup the list of colours, this is needed later for finding colour chains.

#include "Pythia8/JunctionSplitting.h"

namespace Pythia8 {

//==========================================================================

// The JunctionSplitting class.

//--------------------------------------------------------------------------

// Constants: could be changed here if desired, but normally should not.
// These are of technical nature, as described for each.

// For breaking J-J string, pick a Gamma by taking a step with fictitious mass.
const double JunctionSplitting::JJSTRINGM2MAX  = 25.;
const double JunctionSplitting::JJSTRINGM2FRAC = 0.1;

// Iterate junction rest frame boost until convergence or too many tries.
const double JunctionSplitting::CONVJNREST     = 1e-5;
const int JunctionSplitting::NTRYJNREST        = 20;

// Typical average transvere primary hadron mass <mThad>.
const double JunctionSplitting::MTHAD          = 0.9;

// Minimum angle between two partons, to avoid problems with infinities.
const double JunctionSplitting::MINANGLE       = 1e-7;

//--------------------------------------------------------------------------

// Initialize the class and all the created classes.

void JunctionSplitting::init() {

  // Initialize
  colTrace.init(loggerPtr);
  stringLength.init(infoPtr, *settingsPtr);

  // Initialize auxiliary fragmentation classes.
  flavSel.init();
  pTSel.init();
  zSel.init();

  // Initialize string and ministring fragmentation.
  stringFrag.init(&flavSel, &pTSel, &zSel);

  // For junction processing.
  pNormJunction     = parm("StringFragmentation:pNormJunction");
  allowDoubleJunRem = flag("ColourReconnection:allowDoubleJunRem");
}

//--------------------------------------------------------------------------

// Check that all colours are connected in physical way. Also split
// junction pairs, such that the hadronization can handle the configuration.

bool JunctionSplitting::checkColours( Event& event) {

  // Not really a colour check, but a check all numbers are valid.
  for (int i = 0; i < event.size(); ++i)
    if (abs(event[i].px()) >= 0. && abs(event[i].py()) >= 0.
      && abs(event[i].pz()) >= 0.  && abs(event[i].e()) >= 0.
       && abs(event[i].m()) >= 0.);
    else {
       loggerPtr->WARNING_MSG("not-a-number energy/momentum/mass");
       return false;
    }

  // Check if any singlet gluons were made, and if so return false.
  for (int i  = 0; i < event.size(); ++i) {
    if (event[i].isFinal() && event[i].col() != 0 &&
        event[i].col() == event[i].acol()) {
      loggerPtr->WARNING_MSG("made a gluon colour singlet; redoing colours");
      return false;
    }
  }

  // Need to try and split junction structures.
  colTrace.setupColList(event);
  vector<int> iParton;
  vector<vector <int > > iPartonJun, iPartonAntiJun;
  getPartonLists(event, iPartonJun, iPartonAntiJun);

  // Try to split up the junction chains by splitting gluons
  if (!splitJunGluons(event, iPartonJun, iPartonAntiJun) ) {
    loggerPtr->WARNING_MSG(
      "not possible to split junctions; making new colours");
    return false;
  }

  // Remove junctions if more than 2 are connected.
  if (!splitJunChains(event) ) {
    loggerPtr->WARNING_MSG(
      "not possible to split junctions; making new colours");
    return false;
  }

  // Split up junction pairs.
  getPartonLists(event, iPartonJun, iPartonAntiJun);
  if (!splitJunPairs(event, iPartonJun, iPartonAntiJun) ) {
    loggerPtr->WARNING_MSG(
      "not possible to split junctions; making new colours");
    return false;
  }

  // Done checking.
  return true;
}

//--------------------------------------------------------------------------

// Split connected junction chains into separated, mainly by splitting gluons
// into q-qbar pairs. Other methods are applied if the junctions are directly
// connected. (Note: implementation below assumes any intermediate gluons
// are bookkept in iPartonJun, not in iPartonAntiJun.)

bool JunctionSplitting::splitJunGluons(Event& event,
  vector<vector< int > >& iPartonJun, vector<vector< int > >& iPartonAntiJun) {

  // Loop over all junctions and all junction legs.
  for (int iJun = 0; iJun < int(iPartonJun.size()); ++iJun) {

    // Fill in vector of the legs content.
    vector<vector <int> > iJunLegs;
    iJunLegs.resize(3);
    int leg = -1;
    for (int i = 0; i < int(iPartonJun[iJun].size()); ++i) {
      if ( iPartonJun[iJun][i]/10 == iPartonJun[iJun][0]/10) ++leg;
      iJunLegs[leg].push_back(iPartonJun[iJun][i]);
    }

    // Loop over legs.
    for (int i = 0;i < int(iJunLegs.size()); ++i) {

      // If it is not connected to another junction, no need to do anything.
      if (iJunLegs[i].back() > 0)
        continue;
      int identJun = iJunLegs[i][0];
      // If no gluons in between two junctions, not possible to do anything.
      if (iJunLegs[i].size() == 2)
        continue;

      int identAntiJun = 0, iAntiLeg = -1;

      // Pick a new quark at random; for simplicity no diquarks.
      int colQ = 0, acolQ = 0;
      int idQ = flavSel.pickLightQ();

      // If a single gluon in between the two junctions, change it to a
      // quark-anti quark system.
      if ( iJunLegs[i].size() == 3) {

        // Verify that intermediate particle is a gluon (not a gluino).
        if (event[ iJunLegs[i][1] ].idAbs() != 21) continue;

        // Store the new q qbar pair, sharing gluon colour and momentum.
        colQ = event[ iJunLegs[i][1] ].col();
        acolQ = event[ iJunLegs[i][1] ].acol();
        Vec4 pQ = 0.5 * event[ iJunLegs[i][1] ].p();
        double mQ = 0.5 * event[ iJunLegs[i][1] ].m();
        int iQ = event.append( idQ, 75, iJunLegs[i][1], 0, 0, 0, colQ, 0,
          pQ, mQ );
        int iQbar = event.append( -idQ, 75, iJunLegs[i][1], 0, 0, 0, 0, acolQ,
          pQ, mQ );

        // Mark split gluon.
        event[ iJunLegs[i][1] ].statusNeg();
        event[ iJunLegs[i][1] ].daughters( iQ, iQbar);

        // Update junction and antijunction list.
        identAntiJun = iJunLegs[i].back();
        int iOld = iJunLegs[i][1];
        bool erasing = false;
        for (int j = 0; j < int(iPartonJun[iJun].size()); ++j) {
          if (iPartonJun[iJun][j] == iOld)
            erasing = true;
          if (iPartonJun[iJun][j] == identAntiJun) {
            iPartonJun[iJun][j] = iQ;
            break;
          }
          if (erasing) {
            iPartonJun[iJun].erase(iPartonJun[iJun].begin() + j);
            --j;
          }
        }

        // Find the connected antijunction from the list of antijunctions.
        int iAntiJun = -1;
        for (int j = 0; j < int(iPartonAntiJun.size()); j++)
          if ( iPartonAntiJun[j][0]/10 == identAntiJun/10) {
            iAntiJun = j;
            break;
          }
        // If no antijunction found, something went wrong earlier.
        if (iAntiJun == -1) {
           loggerPtr->WARNING_MSG(
            "something went wrong in finding anti-junction");
           return false;
        }

        // Update the antijunction list.
        for (int j = 0; j < int(iPartonAntiJun[iAntiJun].size()); ++j) {
          if ( iPartonAntiJun[iAntiJun][j] / 10 == identAntiJun / 10)
            iAntiLeg++;
          if ( iPartonAntiJun[iAntiJun][j] == identJun) {
            iPartonAntiJun[iAntiJun][j] = iQbar;
            break;
          }
        }
      }
      // If more than a single gluon, decide depending on mass.
      else if (iJunLegs[i].size() > 3) {

        // Evaluate mass-squared for all adjacent gluon pairs.
        vector<double> m2Pair;
        double m2Sum = 0.;
        for (int j = 1; j < int(iJunLegs[i].size()) - 2; ++j) {
          double m2Now = 0.5 * event[ iJunLegs[i][j] ].p()
            * event[ iJunLegs[i][j + 1] ].p();
          m2Pair.push_back(m2Now);
          m2Sum += m2Now;
        }

        // Pick breakup region with probability proportional to mass-squared.
        double m2Reg = m2Sum * rndmPtr->flat();
        int iReg = -1;
        do m2Reg -= m2Pair[++iReg];
        while (m2Reg > 0. && iReg < int(m2Pair.size()) - 1);
        m2Reg = m2Pair[iReg];

        // increase iReg with one, since it should not point towards itself.
        iReg++;

        // Pick breaking point of string in chosen region (symmetrically).
        double m2Temp = min( JJSTRINGM2MAX, JJSTRINGM2FRAC * m2Reg);
        double xPos = 0.5;
        double xNeg = 0.5;
        do {
          double zTemp = zSel.zFrag( idQ, 0, m2Temp);
          xPos = 1. - zTemp;
          xNeg = m2Temp / (zTemp * m2Reg);
        } while (xNeg > 1.);
        if (rndmPtr->flat() > 0.5) swap(xPos, xNeg);

        // Verify that intermediate particles are gluons (not gluinos).
        if ( event[ iJunLegs[i][iReg] ].idAbs() != 21
          || event[ iJunLegs[i][iReg + 1] ].idAbs() != 21 ) continue;

        // Pick up two "mother" gluons of breakup. Mark them decayed.
        Particle& gJun = event[ iJunLegs[i][iReg] ];
        Particle& gAnti = event[ iJunLegs[i][iReg + 1] ];
        gJun.statusNeg();
        gAnti.statusNeg();
        int dau1 = event.size();
        gJun.daughters(dau1, dau1 + 3);
        gAnti.daughters(dau1, dau1 + 3);
        int mother1 = min( iJunLegs[i][iReg], iJunLegs[i][iReg + 1]);
        int mother2 = max( iJunLegs[i][iReg], iJunLegs[i][iReg + 1]);

        // Need to store variables, since it is not safe to use references
        // with append.
        int gJunCol   = gJun.col();
        int gJunAcol  = gJun.acol();
        int gAntiAcol = gAnti.acol();
        Vec4 gJunP    = gJun.p();
        Vec4 gAntiP   = gAnti.p();
        double gJunM  = gJun.m();
        double gAntiM = gAnti.m();

        // Can keep one of old colours but need one new so unambiguous.
        colQ          = gJunAcol;
        acolQ         = event.nextColTag();

         // Store copied gluons with reduced momenta.
        int iGjun = event.append( 21, 75, mother1, mother2, 0, 0,
          gJunCol, gJunAcol, (1. - 0.5 * xPos) * gJunP,
          (1. - 0.5 * xPos) * gJunM);
        event.append( 21, 75, mother1, mother2, 0, 0,
          acolQ, gAntiAcol, (1. - 0.5 * xNeg) * gAntiP,
          (1. - 0.5 * xNeg) * gAntiM);

        // Store the new q qbar pair with remaining momenta.
        int iQ = event.append( idQ, 75, mother1, mother2, 0, 0,
          colQ, 0, 0.5 * xNeg * gAntiP, 0.5 * xNeg * gAntiM );
        int iQbar = event.append( -idQ, 75, mother1, mother2, 0, 0,
          0, acolQ, 0.5 * xPos * gJunP, 0.5 * xPos * gJunM );

        // Update the list of junctions to reflect the splitting.
        identAntiJun = iJunLegs[i].back();
        bool erasing = false;
        for (int j = 0; j < int(iPartonJun[iJun].size()); ++j) {
          if (iPartonJun[iJun][j] == mother1 ||
              iPartonJun[iJun][j] == mother2)
            erasing = true;

          if ( iPartonJun[iJun][j] == identAntiJun) {
            iPartonJun[iJun][j] = iQ;
            iPartonJun[iJun].insert(iPartonJun[iJun].begin() + j, iGjun);
            break;
          }
          if (erasing) {
            iPartonJun[iJun].erase(iPartonJun[iJun].begin() + j);
            j--;
          }
        }

        // Find the connected antijunction from the list of antijunctions.
        int iAntiJun = -1;
        for (int j = 0; j < int(iPartonAntiJun.size());j++)
          if ( iPartonAntiJun[j][0]/10 == identAntiJun/10) {
            iAntiJun = j;
            break;
          }
        // If no antijunction found, something went wrong earlier.
        if (iAntiJun == -1) {
           loggerPtr->WARNING_MSG(
            "something went wrong in finding anti-junction");
           return false;
        }

        // Update the antijunction list to reflect the splitting
        for (int j = 0; j < int(iPartonAntiJun[iAntiJun].size()); ++j) {
          if ( iPartonAntiJun[iAntiJun][j] / 10 == identAntiJun / 10)
            iAntiLeg++;
          if (iPartonAntiJun[iAntiJun][j] == identJun) {
            iPartonAntiJun[iAntiJun][j] = iQbar;
            break;
          }
        }
      }

      // Update end colours for both g -> q qbar and g g -> g g q qbar.
      event.endColJunction((-identJun)/10 - 1, i, colQ);
      event.endColJunction((-identAntiJun)/10 - 1, iAntiLeg, acolQ);
    }
  }

  // Done.
  return true;
}

//--------------------------------------------------------------------------

// Fix chains that contain more than two junctions.
// This is done by removing the minimum needed amount of junctions.
// Might need to make choice based on String length, now randomly chosen.

bool JunctionSplitting::splitJunChains(Event& event) {

  // Get junction chains.
  event.saveJunctionSize();
  vector<vector<int> > junChains = colTrace.getJunChains(event);

  // Remove junctions.
  vector<int> junRem;
  for (int i = 0; i < int(junChains.size()); ++i) {
    if (junChains[i].size() < 3)
      continue;

    vector<int> cols, acols;
    for (int j = 0; j < int(junChains[i].size()); ++j) {

      junRem.push_back(junChains[i][j]);
      if (event.kindJunction(junChains[i][j]) % 2 == 0)
        for (int jLeg = 0; jLeg < 3; ++jLeg)
          acols.push_back(event.colJunction(junChains[i][j],jLeg));
      else
        for (int jLeg = 0; jLeg < 3; ++jLeg)
          cols.push_back(event.colJunction(junChains[i][j],jLeg));
    }

    for (int j = 0; j < int(cols.size()); ++j)
      for (int k = 0; k < int(acols.size()); ++k)
        if (cols[j] == acols[k]) {
          cols.erase(cols.begin() + j);
          acols.erase(acols.begin() + k);
          j--;
          break;
        }

    // Find junctions if we have more colours than anti colours
    while (cols.size() > acols.size()) {
      int i1 = int(rndmPtr->flat() *cols.size());
      int col1 = cols[i1];
      cols.erase(cols.begin() + i1);
      int i2 = int(rndmPtr->flat() *cols.size());
      int col2 = cols[i2];
      cols.erase(cols.begin() + i2);
      int i3 = int(rndmPtr->flat() *cols.size());
      int col3 = cols[i3];
      cols.erase(cols.begin() + i3);
      event.appendJunction(1, col1, col2, col3);
    }

    // Find junctions if we have more anti colours than colours
    while (acols.size() > cols.size()) {
      int i1 = int(rndmPtr->flat() *acols.size());
      int acol1 = acols[i1];
      acols.erase(acols.begin() + i1);
      int i2 = int(rndmPtr->flat() *acols.size());
      int acol2 = acols[i2];
      acols.erase(acols.begin() + i2);
      int i3 = int(rndmPtr->flat() *acols.size());
      int acol3 = acols[i3];
      acols.erase(acols.begin() + i3);
      event.appendJunction(2,acol1,acol2,acol3);
    }

    // If we have more than two colour anti colour pairs
    // form junction antijunction pair.
    while (int(acols.size()) > 1) {
      int i1 = int(rndmPtr->flat() *cols.size());
      int col1 = cols[i1];
      cols.erase(cols.begin() + i1);
      int i2 = int(rndmPtr->flat() *cols.size());
      int col2 = cols[i2];
      cols.erase(cols.begin() + i2);
      int i3 = int(rndmPtr->flat() *acols.size());
      int acol1 = acols[i3];
      acols.erase(acols.begin() + i3);
      int i4 = int(rndmPtr->flat() *acols.size());
      int acol2 = acols[i4];
      acols.erase(acols.begin() + i4);
      int newCol = event.nextColTag();
      event.appendJunction(1, col1, col2, newCol);
      event.appendJunction(2, acol1, acol2, newCol);
    }

    // If we have one colour and one anti colour, form normal string.
    if (int(acols.size()) == 1) {
      int iCol = -1;
      for (int iPar = 0; iPar < event.size(); ++iPar)
        if (event[iPar].isFinal() && event[iPar].col() == cols[0])
          iCol = iPar;
      if (iCol == -1) {
        loggerPtr->WARNING_MSG(
          "splitting multiple directly connected junctions failed");
        return false;
      }
      int iNew = event.copy(iCol, 76);
      event[iNew].col(acols[0]);
    }
  }

  // Delete the junctions from the event record.
  sort(junRem.begin(),junRem.end());
  reverse(junRem.begin(),junRem.end());
  for (int i = 0; i < int(junRem.size()); ++i)
    event.eraseJunction(junRem[i]);
  event.saveJunctionSize();

  return true;
}

//--------------------------------------------------------------------------

// Split junction pairs.
// If it has 3 connections just ignore junctions.
// If it has 2 connections colapse into single string.
// If it has 1 connection, depend on the string length.

bool JunctionSplitting::splitJunPairs(Event& event,
  vector<vector< int > >& iPartonJun, vector<vector< int > >& iPartonAntiJun) {

  // Clear old memory.
  event.saveJunctionSize();
  vector<int> junRem;

  // Get junction chains.
  vector<vector<int> > junChains = colTrace.getJunChains(event);

  for (int i = 0; i < int(junChains.size()); ++i) {
    if (junChains[i].size() == 2) {
      vector<pair<int,int> > matchedLegs;
      for (int j = 0; j < 3; ++j)
        for (int k = 0; k < 3; ++k)
          if (event.colJunction(junChains[i][0],j) ==
              event.colJunction(junChains[i][1],k))
            matchedLegs.push_back(make_pair(j,k));

      // For three connected legs, just remove the junctions,
      // since the pair contains no energy/momentum.
      if (matchedLegs.size() == 3) {
        junRem.push_back(junChains[i][0]);
        junRem.push_back(junChains[i][1]);

      }

      // Split into string if two legs are combined.
      else if (matchedLegs.size() == 2) {

        // Find first leg.
        int i1 = 0;
        if (matchedLegs[0].first != 1 && matchedLegs[1].first != 1) i1 = 1;
        if (matchedLegs[0].first != 2 && matchedLegs[1].first != 2) i1 = 2;

        // Find second leg.
        int j1 = 0;
        if (matchedLegs[0].second != 1 && matchedLegs[1].second != 1) j1 = 1;
        if (matchedLegs[0].second != 2 && matchedLegs[1].second != 2) j1 = 2;

        // Find corresponding colours.
        int col = event.colJunction(junChains[i][0],i1);
        int acol = event.colJunction(junChains[i][1],j1);
        if (event.kindJunction(junChains[i][1]) % 2 == 1)
          swap(col,acol);

        // Find index of anti particle.
        int iAcol = -1;
        for (int j = 0;j < event.size();++j)
          if (event[j].isFinal() && event[j].acol() == acol) {
            iAcol = j;
            break;
          }
        if (iAcol == -1) {
          loggerPtr->WARNING_MSG(
            "anti-colour not found when combining two junctions to a string");
          return false;
        }

        // Update anti colour of anti particle.
        int iNew = event.copy(iAcol,66);
        event[iNew].acol(col);

        // Remove the junctions from the event record.
        junRem.push_back(junChains[i][0]);
        junRem.push_back(junChains[i][1]);
      }

      // Split into string if two legs are combined.
      else if (matchedLegs.size() == 1) {

        // store junction numbers.
        int iJun = junChains[i][0];
        int iAnti = junChains[i][1];
        int iLeg = matchedLegs[0].first;
        int iAntiLeg = matchedLegs[0].second;
        if (event.kindJunction(iAnti) % 2 == 1)  {
          swap(iJun, iAnti);
          swap(iLeg, iAntiLeg);
        }

        // Find the junctions in the parton list.
        int iJunList = -1, iAntiList = -1;
        for (int l = 0;l < int(iPartonJun.size()); ++l)
          if (- iPartonJun[l][0]/10 - 1 == iJun) {
            iJunList = l;
            break;
        }
        for (int l = 0;l < int(iPartonAntiJun.size()); ++l)
          if (- iPartonAntiJun[l][0]/10 - 1 == iAnti) {
            iAntiList = l;
            break;
        }
        if (iJunList == -1 || iAntiList == -1) {
          loggerPtr->ERROR_MSG("failed to find junctions in the parton list");
          return false;
        }

        // Fill in vector of the legs content.
        vector<vector <int> > iJunLegs;
        iJunLegs.resize(3);
        int leg = -1;

        for (int l = 0; l < int(iPartonJun[iJunList].size()); ++l) {
          if ( iPartonJun[iJunList][l]/10 == iPartonJun[iJunList][0]/10) ++leg;
          iJunLegs[leg].push_back(iPartonJun[iJunList][l]);
        }

         // Fill in vector of the legs content.
        vector<vector <int> > iAntiLegs;
        iAntiLegs.resize(3);
        leg = -1;
        for (int l = 0; l < int(iPartonAntiJun[iAntiList].size()); ++l) {
          if ( iPartonAntiJun[iAntiList][l]/10
            == iPartonAntiJun[iAntiList][0]/10) ++leg;
          iAntiLegs[leg].push_back(iPartonAntiJun[iAntiList][l]);
        }

        // Identify the two external legs of either junction.
        vector<int>& iJunLeg0 = (iLeg == 0) ? iJunLegs[1] : iJunLegs[0];
        vector<int>& iJunLeg1 = (iLeg == 2) ? iJunLegs[1] : iJunLegs[2];
        vector<int>& iAntiLeg0 = (iAntiLeg == 0) ? iAntiLegs[1] : iAntiLegs[0];
        vector<int>& iAntiLeg1 = (iAntiLeg == 2) ? iAntiLegs[1] : iAntiLegs[2];

        // Check that the anti-leg is not stored as a junction.
        // This should only happen for gluinos, which are not split.
        if (iAntiLeg0[1] < 0 || iAntiLeg1[1] < 0)  continue;

        // Simplified procedure: mainly study first parton on each leg.
        Vec4 pJunLeg0 = event[ iJunLeg0[1] ].p();
        Vec4 pJunLeg1 = event[ iJunLeg1[1] ].p();
        Vec4 pAntiLeg0 = event[ iAntiLeg0[1] ].p();
        Vec4 pAntiLeg1 = event[ iAntiLeg1[1] ].p();

        // Check that no two legs are parallel.
        if ( theta(pJunLeg0,pJunLeg1)   < MINANGLE
          || theta(pAntiLeg0,pAntiLeg1) < MINANGLE
          || theta(pJunLeg0,pAntiLeg0)  < MINANGLE
          || theta(pJunLeg0,pAntiLeg1)  < MINANGLE
          || theta(pJunLeg1,pAntiLeg0)  < MINANGLE
          || theta(pJunLeg1,pAntiLeg1)  < MINANGLE) {
          loggerPtr->WARNING_MSG("parallel junction state not allowed");
          return false;
        }

        // For third side add both legs beyond other junction, weighted.
        double wt2 = 1. - exp( -pAntiLeg0.pAbs() / pNormJunction);
        double wt3 = 1. - exp( -pAntiLeg1.pAbs() / pNormJunction);
        Vec4 pCombined = wt2 * pAntiLeg0 + wt3 * pAntiLeg1;
        Vec4 vJun = stringFrag.junctionRestFrame(pJunLeg0, pJunLeg1,
                           pCombined);
        vJun /= sqrt( 1 - vJun.pAbs2() );
        Vec4 pInJRF[3] = { pJunLeg0, pJunLeg1, pCombined};
        for (int l = 0; l < 3; ++l) pInJRF[l].bstback(vJun);

        // For third side add both legs beyond other junction, weighted.
        wt2 = 1. - exp( -pJunLeg0.pAbs() / pNormJunction);
        wt3 = 1. - exp( -pJunLeg1.pAbs() / pNormJunction);
        pCombined = wt2 * pJunLeg0 + wt3 * pJunLeg1;
        Vec4 vAnti = stringFrag.junctionRestFrame(pAntiLeg0, pAntiLeg1,
                      pCombined);
        vAnti /= sqrt( 1 - vAnti.pAbs2() );
        Vec4 pInARF[3] = { pAntiLeg0, pAntiLeg1, pCombined};
        for (int l = 0; l < 3; ++l) pInARF[l].bstback(vAnti);

        // Calculate Lambda-measure length of three possible topologies.
        double vJvA = vJun * vAnti;
        double vJvAe2y = vJvA + sqrt(vJvA*vJvA - 1.);
        double lambdaJA = stringLength.getJuncLength(
          pInJRF[0], pInJRF[1], pInARF[0], pInARF[1]);

        double lambda00 = stringLength.getStringLength(pJunLeg0, pAntiLeg0)
          + stringLength.getStringLength(pJunLeg1, pAntiLeg1);

        double lambda01 = stringLength.getStringLength(pJunLeg0, pAntiLeg1)
          + stringLength.getStringLength(pJunLeg1, pAntiLeg0);

        // Case when either topology without junctions is the shorter one.
        if (lambdaJA > min( lambda00, lambda01) && allowDoubleJunRem) {

          // Find indices of particles.
          int iCol1 = iJunLeg0[1];
          int iCol2 = iJunLeg1[1];
          int iAcol1 = iAntiLeg0[1];
          int iAcol2 = iAntiLeg1[1];
          if (lambda00 > lambda01)
            swap(iAcol1, iAcol2);

          // Change the colour index and mark junctions to be removed.
          int iNew1 = event.copy(iAcol1, 66);
          event[iNew1].acol(event[iCol1].col());
          int iNew2 = event.copy(iAcol2, 66);
          event[iNew2].acol(event[iCol2].col());
          junRem.push_back(junChains[i][0]);
          junRem.push_back(junChains[i][1]);
          continue;
        }

        // Case where junction and antijunction to be separated.
        // Shuffle (p+/p-)  momentum of order <mThad> between systems,
        // times 2/3 for 120 degree in JRF, times 1/2 for two legs,
        // but not more than half of what nearest parton carries.

        // Only allow to take momentum from non-heavy resonances
        // (i.e. id 1-5 and 21 and diquarks). If none is available
        // return false.
        int idJ0Abs = event[ iJunLeg0[1] ].idAbs();
        bool acceptJ0 = idJ0Abs < 6 || idJ0Abs == 21 ||
          (idJ0Abs > 1000 && idJ0Abs < 6000 && (idJ0Abs / 10) % 10 == 0);
        int idJ1Abs = event[ iJunLeg1[1] ].idAbs();
        bool acceptJ1 = idJ1Abs < 6 || idJ1Abs == 21 ||
          (idJ1Abs > 1000 && idJ1Abs < 6000 && (idJ1Abs / 10) % 10 == 0);
        int idA0Abs = event[ iAntiLeg0[1] ].idAbs();
        bool acceptA0 = idA0Abs < 6 || idA0Abs == 21 ||
          (idA0Abs > 1000 && idA0Abs < 6000 && (idA0Abs / 10) % 10 == 0);
        int idA1Abs = event[ iAntiLeg1[1] ].idAbs();
        bool acceptA1 = idA1Abs < 6 || idA1Abs == 21 ||
          (idA1Abs > 1000 && idA1Abs < 6000 && (idA1Abs / 10) % 10 == 0);

        if ( !(acceptJ0 || acceptJ1)) {
          loggerPtr->WARNING_MSG(
            "no light quarks available in junction split");
          return false;
        }

        if ( !(acceptA0 || acceptA1)) {
          loggerPtr->WARNING_MSG(
            "no light quarks available in junction split");
          return false;
        }

        double eShift = MTHAD / (3. * sqrt(vJvAe2y));
        double fracJ0 = 0, fracJ1 = 0, fracA0 = 0, fracA1 = 0;
        if (acceptJ0)
          fracJ0 = min(0.5, eShift / pInJRF[0].e());
        if (acceptJ1)
          fracJ1 = min(0.5, eShift / pInJRF[1].e());
        Vec4 pFromJun = fracJ0 * pJunLeg0 + fracJ1 * pJunLeg1;
        if (acceptA0)
          fracA0 = min(0.5, eShift / pInARF[0].e());
        if (acceptA1)
          fracA1 = min(0.5, eShift / pInARF[1].e());
        Vec4 pFromAnti = fracA0 * pAntiLeg0 + fracA1 * pAntiLeg1;

        // Pick a new quark at random; for simplicity no diquarks.
        int idQ = flavSel.pickLightQ();

        int junMother1  = min(iJunLeg0[1], iJunLeg1[1]);
        int junMother2  = max(iJunLeg0[1], iJunLeg1[1]);
        int antiMother1 = min(iAntiLeg0[1], iAntiLeg1[1]);
        int antiMother2 = max(iAntiLeg0[1], iAntiLeg1[1]);

        // Copy junction partons with scaled-down momenta and update legs.
        int iJunNew1 = event.copy(iJunLeg0[1], 76);
        event[iJunNew1].rescale5(1. - fracJ0);
        iJunLeg0[1] = iJunNew1;
        event[iJunNew1].mothers(junMother2, junMother1);

        int iJunNew2 = event.copy(iJunLeg1[1], 76);
        event[iJunNew2].rescale5(1. - fracJ1);
        iJunLeg1[1] = iJunNew2;
        event[iJunNew2].mothers(junMother2, junMother1);

        // Update antijunction anticolour and store antiquark with junction
        // momentum.
        int acolQ = event.nextColTag();
        event.endColJunction(iAnti, iAntiLeg, acolQ);
        event.colJunction(iAnti, iAntiLeg, acolQ);
        int iNewA = event.append( -idQ, 76, junMother2, junMother1, 0, 0,
          0, acolQ, pFromJun, pFromJun.mCalc() );

        // Copy antijunction partons with scaled-down momenta and update legs.
        int iAntiNew1 = event.copy(iAntiLeg0[1], 76);
        event[iAntiNew1].rescale5(1. - fracA0);
        iAntiLeg0[1] = iAntiNew1;
        event[iAntiNew1].mothers(antiMother2, antiMother1);

        int iAntiNew2 = event.copy(iAntiLeg1[1], 76);
        event[iAntiNew2].rescale5(1. - fracA1);
        iAntiLeg1[1] = iAntiNew2;
        event[iAntiNew2].mothers(antiMother2, antiMother1);

        // Update junction colour and store quark with antijunction momentum.
        int colQ = event.nextColTag();
        event.endColJunction(iJun, iLeg, colQ);
        event.colJunction(iJun, iLeg, colQ);
        int iNewJ = event.append( idQ, 76, antiMother2, antiMother1, 0, 0,
          colQ, 0, pFromAnti, pFromAnti.mCalc() );

        // Set daughters.
        event[event[iJunNew1].mother1()].daughters( iJunNew1, iNewA);
        event[event[iJunNew1].mother2()].daughters( iJunNew1, iNewA);
        event[event[iAntiNew1].mother1()].daughters( iAntiNew1, iNewJ);
        event[event[iAntiNew1].mother2()].daughters( iAntiNew1, iNewJ);

        // Done with splitting junction from antijunction.
      }
    }
  }

  // Delete the junctions from the event record.
  sort(junRem.begin(),junRem.end());
  reverse(junRem.begin(),junRem.end());
  for (int i = 0; i < int(junRem.size()); ++i)
    event.eraseJunction(junRem[i]);
  event.saveJunctionSize();

  // Done.
  return true;
}

//--------------------------------------------------------------------------

// Create the lists of partons connected to junctions.
// Input: event
// Output: iPartonJun and iPartonAntiJun (for JJ connections).

bool JunctionSplitting::getPartonLists(Event& event,
  vector<vector< int > > & iPartonJun, vector<vector<int > >& iPartonAntiJun) {

  // Need to try and split junction structures.
  colTrace.setupColList(event);
  vector<int> iParton;
  iPartonJun.clear();
  iPartonAntiJun.clear();

  // Loop over junctions (first junctions, then antijunctions; this ensures
  // that any gluons between junction-antijunction pairs are bookkept in
  // iPartonJun rather than in iPartonAntiJun; assumed by subsequent
  // splitting methods).
  for (int iLoop = 0; iLoop < 2*event.sizeJunction(); ++iLoop) {

    int iJun = iLoop % event.sizeJunction();
    if ( !event.remainsJunction(iJun)) continue;

    // Do junctions first, then antijunctions
    int kindJun = event.kindJunction(iJun);
    if ( iLoop < event.sizeJunction() && kindJun%2 == 0) continue;
    else if ( iLoop >= event.sizeJunction() && kindJun%2 == 1) continue;

    iParton.resize(0);
    // Loop over junction legs
    // Afterwards (using iJun=0 as example), iParton should look something
    // like {-10, parton-chain-to-first-colour-end, -11, parton-chain-to-next-
    // colour-end, -12, parton-chain-to-last-colour-end}
    // (Note: -10*N indicates leg connected to another junction, with iJun=N-1)
    for (int iCol = 0; iCol < 3; ++iCol) {
      int indxCol = event.colJunction(iJun, iCol);
      iParton.push_back( -(10 + 10 * iJun + iCol) );
      // Junctions: check that we can find colour ends.
      if ( kindJun%2 == 1 && !colTrace.traceFromAcol(indxCol,event, iJun,
        iCol, iParton) ) return false;
      // Antijunctions: check that we can find anticolour ends.
      else if ( kindJun%2 == 0 && !colTrace.traceFromCol(indxCol,event, iJun,
        iCol, iParton) ) return false;
    }

    // Save lists for junction-junction systems in iPartonJun, iPartonAntiJun
    int nNeg = 0;
    for (int i = 0; i < int(iParton.size()); ++i) if (iParton[i] < 0) ++nNeg;
    // If 3
    if (nNeg > 3) {
      if ( kindJun%2 == 1 ) iPartonJun.push_back(iParton);
      else iPartonAntiJun.push_back(iParton);
    }

  }

  // Done.
  return true;

}

//--------------------------------------------------------------------------

// Change the anticolour of the particle that has acol to be col.

bool JunctionSplitting::setAcol(Event& event, int col, int acol) {

  // Update anticolour if it belongs to a particle.
  for (int j = 0;j < event.size(); ++j)
    if (event[j].isFinal() && event[j].acol() == acol) {
      int iNew = event.copy(j,66);
      event[iNew].acol(col);
      return true;
    }
  // Check if antijunction is connected to a junction.
  for (int j = 0;j < event.sizeJunction(); ++j)
    for (int jLeg = 0;jLeg < 3; ++jLeg)
      if (event.colJunction(j, jLeg) == acol) {
        event.colJunction(j, jLeg, col);
        return true;
      }

  // If no acol was found something went wrong.
  loggerPtr->WARNING_MSG(
     "anti-colour not found when combing two junctions to a string");
  return false;
}

//==========================================================================

} // end namespace Pythia8