// MergingHooks.cc is a part of the PYTHIA event generator.
// Copyright (C) 2012 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This file is written by Stefan Prestel.
// Function definitions (not found in the header) for the HardProcess and
// MergingHooks classes.

#include "MergingHooks.h"

namespace Pythia8 {
 
//==========================================================================

// The HardProcess class.

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

// Declaration of hard process class
// This class holds information on the desired hard 2->2 process to be merged
// This class is a container class for History class use.

// Initialisation on the process string

void HardProcess::initOnProcess( string process, ParticleData* particleData) {
  state.init("(hard process)", particleData);
  translateProcessString(process);
}

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

// Initialisation on the path to LHE file

void HardProcess::initOnLHEF( string LHEfile, ParticleData* particleData) {
  state.init("(hard process)", particleData);
  translateLHEFString(LHEfile);
}

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

// Function to access the LHE file and read relevant information.
// The merging scale will be read from the +1 jet sample, called
//   LHEpath_1.lhe
// while the hard process will be read from
//   LHEpath_0.lhe
// Currently, only read from MadEvent- or Sherpa-generated LHE files
// is automatic, else the user is asked to supply the necessary
// information.

void HardProcess::translateLHEFString( string LHEpath){

  // Open path to LHEF and extract merging scale
  ifstream infile;
  infile.open( (char*)( LHEpath +"_0.lhe").c_str());

  // Check with ME generator has been used
  int iLine =0;
  int nLinesMax = 200;
  string lineGenerator;
  while( iLine < nLinesMax
        && lineGenerator.find("SHERPA", 0) == string::npos
        && lineGenerator.find("POWHEG-BOX", 0) == string::npos
        && lineGenerator.find("Pythia8", 0) == string::npos
        && lineGenerator.find("MadGraph", 0) == string::npos){
    iLine++;
    lineGenerator = " ";
    getline(infile,lineGenerator);
  }
  infile.close();

  vector <int> incom;
  vector <int> inter;
  vector <int> outgo;
  // Particle identifiers, ordered in such a way that e.g. the "u"
  // in a mu is not mistaken for an u quark
  int inParticleNumbers[] = {
        // Leptons
        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
        // Jet container
        2212,2212,0,0,0,0,
        // Quarks
        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6};

  string inParticleNamesSH[] = {
        // Leptons
        "-11","11","-12","12","-13","13","-14","14","-15","15","-16","16",
        // Jet container
        "-93","93","-90","90","-91","91",
        // Quarks
        "-1","1","-2","2","-3","3","-4","4","-5","5","-6","6"};
  string inParticleNamesMG[] =  {
        // Leptons
        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
        // Jet container
        "p~","p","l+","l-","vl~","vl",
        // Quarks
        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t"};

  // Declare intermediate particle identifiers
  int interParticleNumbers[] = {
         // Electroweak gauge bosons
         22,23,-24,24,25,2400,
         // Top quarks
         -6,6,
         // Dummy index as back-up
         0,
         // All squarks
        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
  // Declare names of intermediate particles
  string interParticleNamesMG[] = {
        // Electroweak gauge bosons
        "a","z","w-","w+","h","W",
         // Top quarks
         "t~","t",
        // Dummy index as back-up
        "xx",
         // All squarks
         "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
         "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"}; 

  // Declare final state particle identifiers
  int outParticleNumbers[] = {
        // Leptons
        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
        // Jet container and lepton containers
        2212,2212,0,0,0,0,1200,1100,5000,
        // Quarks
        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6,
        // SM uncoloured bosons
        22,23,-24,24,25,2400,
        // Neutralino in SUSY
        1000022,
        // All squarks
        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
  // Declare names of final state particles
  string outParticleNamesMG[] =  {
        // Leptons
        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
        // Jet container and lepton containers
        "j~","j","l+","l-","vl~","vl","NEUTRINOS","LEPTONS","BQUARKS",
        // Quarks
        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t",
        // SM uncoloured bosons
        "a","z","w-","w+","h","W",
        // Neutralino in SUSY
        "n1",
        // All squarks
        "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
        "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"};

  string outParticleNamesSH[] = {
        // Leptons
        "-11","11","-12","12","-13","13","-14","14","-15","15","-16","16",
        // Jet container and lepton containers
        "-93","93","-90","90","-91","91","0","0","0",
        // Quarks
        "-1","1","-2","2","-3","3","-4","4","-5","5","-6","6",
        // SM uncoloured bosons
        "22","23","-24","24","25","0",
        // Neutralino in SUSY
        "1000022",
        // All squarks
        "-1000001","1000001","-1000002","1000002","-1000003","1000003",
                   "-1000004","1000004",
        "-1000005","1000005","-1000006","1000006","-2000001","2000001",
                   "-2000002","2000002",
        "-2000003","2000003","-2000004","2000004","-2000005","2000005",
                   "-2000006","2000006"};

  // Declare size of particle name arrays
  int nIn   = 30;
  int nInt  = 33;
  int nOut  = 64;

  // Save type of the generator, in order to be able to extract
  // the tms definition
  meGenType = (lineGenerator.find("MadGraph", 0) != string::npos)     ? -1
              : (lineGenerator.find("SHERPA", 0) != string::npos)     ? -2
              : (lineGenerator.find("POWHEG-BOX", 0) != string::npos) ? -3
              : (lineGenerator.find("Pythia8", 0) != string::npos)    ? -4
              : 0;

  if (meGenType == -2){
    // Now read merging scale
    // Open path to LHEF and extract merging scale
    infile.open( (char*)( LHEpath +"_1.lhe").c_str());
    string lineTMS;
    while(lineTMS.find("NJetFinder ", 0) == string::npos){
      lineTMS = " ";
      getline(infile,lineTMS);
    }
    infile.close();
    lineTMS = lineTMS.substr(0,lineTMS.find(" 0.0 ",0));
    lineTMS = lineTMS.substr(lineTMS.find(" ", 0)+3,lineTMS.size());
    // Remove whitespaces
    while(lineTMS.find(" ", 0) != string::npos)
      lineTMS.erase(lineTMS.begin()+lineTMS.find(" ",0));
    // Replace d with e
    if ( lineTMS.find("d", 0) != string::npos)
      lineTMS.replace(lineTMS.find("d", 0),1,1,'e');
    tms = atof((char*)lineTMS.c_str());

    // Now read hard process
    // Open path to LHEF and extract hard process
    infile.open( (char*)( LHEpath +"_0.lhe").c_str());
    string line;
    while(line.find("Process", 0) == string::npos){
      line = " ";
      getline(infile,line);
    }
    infile.close();
    line = line.substr(line.find(" ",0),line.size());

    // Cut string into incoming and outgoing pieces 
    vector <string> pieces;
    pieces.push_back( line.substr(0,line.find("->", 0)) );
    // Do not count additional final jets
    int end = (line.find("{", 0) != string::npos) ? line.find("{", 0)-2
            : line.size();
    pieces.push_back( line.substr(line.find("->", 0)+2, end) );

    // Get incoming particles
    for(int i=0; i < nIn; ++i) {
      for(int n = pieces[0].find(inParticleNamesSH[i], 0);
             n != int(string::npos);
             n = pieces[0].find(inParticleNamesSH[i], n)) {
        incom.push_back(inParticleNumbers[i]);
        pieces[0].erase(pieces[0].begin()+n,
                        pieces[0].begin()+n+inParticleNamesSH[i].size());
        n=0;
      }
    }
    // Get intermediate particles
    // If intermediates are still empty, fill intermediate with default value
    inter.push_back(0);
    // Get final particles
    for(int i=0; i < nOut; ++i) {
      for(int n = pieces[1].find(outParticleNamesSH[i], 0);
             n != int(string::npos);
             n = pieces[1].find(outParticleNamesSH[i], n)) {
        outgo.push_back(outParticleNumbers[i]);
        pieces[1].erase(pieces[1].begin()+n,
                        pieces[1].begin()+n+outParticleNamesSH[i].size());
        n=0;
      }
    }

  } else if (meGenType == -1 || meGenType == -3 || meGenType == -4){

    // Now read merging scale
    string lineTMS;

    if (meGenType == -1) {
      // Open path to LHEF and extract merging scale
      infile.open( (char*)( LHEpath +"_1.lhe").c_str());
      while(lineTMS.find("ktdurham", 0) == string::npos){
        lineTMS = " ";
        getline(infile,lineTMS);
      }
      lineTMS = lineTMS.substr(0,lineTMS.find("=",0));
      infile.close();
    } else {
      lineTMS = "30.";
    }

    // Remove whitespaces
    while(lineTMS.find(" ", 0) != string::npos)
      lineTMS.erase(lineTMS.begin()+lineTMS.find(" ",0));
    // Replace d with e
    if ( lineTMS.find("d", 0) != string::npos)
      lineTMS.replace(lineTMS.find("d", 0),1,1,'e');
    tms = atof((char*)lineTMS.c_str());

    // Now read hard process
    // Open path to LHEF and extract hard process
    infile.open( (char*)( LHEpath +"_0.lhe").c_str());
    string line;
    while(line.find("@1", 0) == string::npos){
      line = " ";
      getline(infile,line);
    }
    infile.close();
    line = line.substr(0,line.find("@",0));

    // Count number of resonances
    int appearances = 0;
    for(int n = line.find("(", 0); n != int(string::npos);
            n = line.find("(", n)) {
      appearances++;
      n++;
    }

    // Cut string in incoming, resonance+decay and outgoing pieces 
    vector <string> pieces;
    for(int i =0; i < appearances;++i) {
      int n = line.find("(", 0);
      pieces.push_back(line.substr(0,n));
      line = line.substr(n+1,line.size());
    }
    // Cut last resonance from rest
    if (appearances > 0) {
      pieces.push_back( line.substr(0,line.find(")",0)) );
      pieces.push_back( line.substr(line.find(")",0)+1,line.size()) );
    }

    // If the string was not cut into pieces, i.e. no resonance was
    // required, cut string using '>' as delimiter
    if (pieces.empty() ){
      appearances = 0;
      for(int n = line.find(">", 0); n != int(string::npos);
              n = line.find(">", n)) {
        appearances++;
        n++;
      }

      // Cut string in incoming and outgoing pieces 
      for(int i =0; i < appearances;++i) {
        int n = line.find(">", 0);
        pieces.push_back(line.substr(0,n));
        line = line.substr(n+1,line.size());
      }

      if (appearances == 1) pieces.push_back(line);
      if (appearances > 1) {
        pieces.push_back( line.substr(0,line.find(">",0)) );
        pieces.push_back( line.substr(line.find(">",0)+1,line.size()) );
      }
    }

    // Get incoming particles
    for(int i=0; i < nIn; ++i) {
      for(int n = pieces[0].find(inParticleNamesMG[i], 0);
             n != int(string::npos);
             n = pieces[0].find(inParticleNamesMG[i], n)) {
        incom.push_back(inParticleNumbers[i]);
        pieces[0].erase(pieces[0].begin()+n,
                        pieces[0].begin()+n+inParticleNamesMG[i].size());
        n=0;
      }
    }

    // Check intermediate resonances and decay products
    for(int i =1; i < int(pieces.size()); ++i){
      // Seperate strings into intermediate and outgoing, if not already done
      int k = pieces[i].find(">", 0);

      string intermediate = (pieces[i].find(">", 0) != string::npos) ? 
                             pieces[i].substr(0,k) : "";
      string outgoing = (pieces[i].find(">", 0) != string::npos) ? 
                         pieces[i].substr(k+1,pieces[i].size()) : pieces[i];

      // Get intermediate particles
      for(int j=0; j < nInt; ++j) {
        for(int n = intermediate.find(interParticleNamesMG[j], 0);
               n != int(string::npos);
               n = intermediate.find(interParticleNamesMG[j], n)) {
          inter.push_back(interParticleNumbers[j]);
          intermediate.erase(intermediate.begin()+n,
                      intermediate.begin()+n+interParticleNamesMG[j].size());
          n=0;
        }
      }

      // Get outgoing particles
      for(int j=0; j < nOut; ++j) {
        for(int n = outgoing.find(outParticleNamesMG[j], 0);
               n != int(string::npos);
               n = outgoing.find(outParticleNamesMG[j], n)) {
          outgo.push_back(outParticleNumbers[j]);
          outgoing.erase(outgoing.begin()+n,
                         outgoing.begin()+n+outParticleNamesMG[j].size());
          n=0;
        }
      }

      // For arbitrary or non-existing intermediate, remember zero for each 
      // two outgoing particles, without bosons.
      if (inter.empty()) {

        // For final state bosons, bookkeep the final state boson as
        // intermediate as well.
        int nBosons = 0;
        for(int l=0; l < int(outgo.size()); ++l)
          if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
            nBosons++;

        int nZeros = (outgo.size() - nBosons)/2;
        for(int l=0; l < nZeros; ++l)
          inter.push_back(0);
      }

      // For final state bosons, bookkeep the final state boson as
      // intermediate as well.
      for(int l=0; l < int(outgo.size()); ++l)
        if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
          inter.push_back(outgo[l]);

    }

  } else {

    cout << "Reading of tms and hard process information from LHEF currently"
         << " only automated for MadEvent- or SHERPA-produced LHEF" << endl;
    int tempInt      = 0;
    cout << "Use default process pp -> e+ve + jets? (0:no / 1:yes): ";
    cin >> tempInt;
    cout << endl;

    if (tempInt == 0){
      tempInt = 0;
      double tempDouble  = 0.0;
      cout << "Please specify merging scale (kT Durham, in GeV): ";
      cin >> tempDouble;
      tms = tempDouble;
      meGenType = -1;
      cout << endl;
      cout << "Please specify first incoming particle ";
      cout << "(p+/p- = 2212, e- = 11, e+ = -11): ";
      cin >> tempInt;
      incom.push_back(tempInt);
      tempInt = 0;
      cout << endl;
      cout << "Please specify second incoming particle ";
      cout << "(p+/p- = 2212, e- = 11, e+ = -11): ";
      cin >> tempInt;
      incom.push_back(tempInt);
      tempInt = 0;
      cout << endl;
      cout << "Please specify intermediate particle, if any ";
      cout << "(0 for none, else PDG code): ";
      cin >> tempInt;
      inter.push_back(tempInt);
      cout << endl;
      do {
        tempInt = 0;
        cout << "Please specify outgoing particle ";
        cout << "(jet=2212, else PDG code, exit with 99): ";
        cin >> tempInt;
        if (tempInt != 99) outgo.push_back(tempInt);
      } while(tempInt != 99);
      cout << endl;
    } else {
      cout << "LHE file not produced by SHERPA or MG/ME - ";
      cout << "Using default process and tms" << endl;
      incom.push_back(2212);
      incom.push_back(2212);
      inter.push_back(24);
      outgo.push_back(-11);
      outgo.push_back(12);
      tms = 10.;
      meGenType = -1;
    }
  }

  // Now store incoming, intermediate and outgoing
  // Set intermediate tags
  for(int i=0; i < int(inter.size()); ++i)
    hardIntermediate.push_back(inter[i]);

  // Set the incoming particle tags
  if (incom.size() != 2)
    cout << "Only two incoming particles allowed" << endl;
  else {
    hardIncoming1 = incom[0];
    hardIncoming2 = incom[1];
  }

  // Remember final state bosons
  int nBosons = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( (abs(outgo[i]) > 20 && abs(outgo[i]) <= 25) || outgo[i] == 2400)
      nBosons++;
  // Remember b-quark container
  int nBQuarks = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 5000)
      nBQuarks++;
  // Remember jet container
  int nJets = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 2212)
      nJets++;
  // Remember lepton container
  int nLeptons = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 1100)
      nLeptons++;
  // Remember lepton container
  int nNeutrinos = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 1200)
      nNeutrinos++;
  int nContainers = nLeptons + nNeutrinos + nJets + nBQuarks;

  // Set final particle identifiers
  if ( (outgo.size() - nBosons - nContainers)%2 == 1) {
    cout << "Only even number of outgoing particles allowed" << endl;  
    for(int i=0; i < int(outgo.size()); ++i)
      cout << outgo[i] << endl;
  } else {

    // Push back particles / antiparticles
    for(int i=0; i < int(outgo.size()); ++i)
      if (outgo[i] > 0
        && outgo[i] != 2212
        && outgo[i] != 5000
        && outgo[i] != 1100
        && outgo[i] != 1200
        && outgo[i] != 2400
        && outgo[i] != 1000022)
        hardOutgoing2.push_back( outgo[i]);
      else if (outgo[i] < 0)
        hardOutgoing1.push_back( outgo[i]);

    // Save final state W-boson container as particle
    for(int i=0; i < int(outgo.size()); ++i)
      if ( outgo[i] == 2400)
        hardOutgoing2.push_back( outgo[i]);

    // Push back jets, distribute evenly amongst particles / antiparticles
    // Push back majorana particles, distribute evenly
    int iNow = 0;
    for(int i=0; i < int(outgo.size()); ++i)
      if ( (outgo[i] == 2212
        || outgo[i] == 5000
        || outgo[i] == 1200
        || outgo[i] == 1000022)
        && iNow%2 == 0 ){
        hardOutgoing2.push_back( outgo[i]);
        iNow++;
      } else if ( (outgo[i] == 2212
               || outgo[i] == 5000
               || outgo[i] == 1100
               || outgo[i] == 1000022)
               && iNow%2 == 1 ){
        hardOutgoing1.push_back( outgo[i]);
        iNow++;
      }
  }

  // Done
}

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

// Function to translate a string specitying the core process into the
// internal notation
// Currently, the input string has to be in MadEvent notation

void HardProcess::translateProcessString( string process){

  vector <int> incom;
  vector <int> inter;
  vector <int> outgo;
  // Particle identifiers, ordered in such a way that e.g. the "u"
  // in a mu is not mistaken for an u quark
  int inParticleNumbers[] = {
        // Leptons
        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
        // Jet container
        2212,2212,0,0,0,0,
        // Quarks
        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6};
  string inParticleNamesMG[] =  {
        // Leptons
        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
        // Jet container
        "p~","p","l+","l-","vl~","vl",
        // Quarks
        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t"};

  // Declare intermediate particle identifiers
  int interParticleNumbers[] = {
         // Electroweak gauge bosons
         22,23,-24,24,25,2400,
         // Top quarks
         -6,6,
         // Dummy index as back-up
         0,
         // All squarks
        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
  // Declare names of intermediate particles
  string interParticleNamesMG[] = {
        // Electroweak gauge bosons
        "a","z","w-","w+","h","W",
         // Top quarks
         "t~","t",
        // Dummy index as back-up
        "xx",
         // All squarks
         "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
         "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"}; 

  // Declare final state particle identifiers
  int outParticleNumbers[] = {
        // Leptons
        -11,11,-12,12,-13,13,-14,14,-15,15,-16,16,
        // Jet container and lepton containers
        2212,2212,0,0,0,0,1200,1100,5000,
        // Quarks
        -1,1,-2,2,-3,3,-4,4,-5,5,-6,6,
        // SM uncoloured bosons
        22,23,-24,24,25,2400,
        // Neutralino in SUSY
        1000022,
        // All squarks
        -1000001,1000001,-1000002,1000002,-1000003,1000003,-1000004,1000004,
        -1000005,1000005,-1000006,1000006,-2000001,2000001,-2000002,2000002,
        -2000003,2000003,-2000004,2000004,-2000005,2000005,-2000006,2000006};
  // Declare names of final state particles
  string outParticleNamesMG[] =  {
        // Leptons
        "e+","e-","ve~","ve","mu+","mu-","vm~","vm","ta+","ta-","vt~","vt",
        // Jet container and lepton containers
        "j~","j","l+","l-","vl~","vl","NEUTRINOS","LEPTONS","BQUARKS",
        // Quarks
        "d~","d","u~","u","s~","s","c~","c","b~","b","t~","t",
        // SM uncoloured bosons
        "a","z","w-","w+","h","W",
        // Neutralino in SUSY
        "n1",
        // All squarks
        "dl~","dl","ul~","ul","sl~","sl","cl~","cl","b1~","b1","t1~","t1",
        "dr~","dr","ur~","ur","sr~","sr","cr~","cr","b2~","b2","t2~","t2"};

  // Declare size of particle name arrays
  int nIn   = 30;
  int nInt  = 33;
  int nOut  = 64;

  // Start mapping user-defined particles onto particle ids.
  //string fullProc = "pp>{blaa,124}LEPTONS,NEUTRINOS";
  string fullProc = process;

  // Find user-defined hard process content
  // Count number of user particles
  int nUserParticles = 0;
  for(int n = fullProc.find("{", 0); n != int(string::npos);
          n = fullProc.find("{", n)) {
    nUserParticles++;
    n++;
  }
  // Cut user-defined particles from remaining process 
  vector <string> userParticleStrings;
  for(int i =0; i < nUserParticles;++i) {
    int n = fullProc.find("{", 0);
    userParticleStrings.push_back(fullProc.substr(0,n));
    fullProc = fullProc.substr(n+1,fullProc.size());
  }
  // Cut remaining process string from rest
  if (nUserParticles > 0)
    userParticleStrings.push_back(
      fullProc.substr( 0, fullProc.find("}",0) ) );
  // Remove curly brackets and whitespace
  for(int i =0; i < int(userParticleStrings.size());++i) {
    while(userParticleStrings[i].find("{", 0) != string::npos)
      userParticleStrings[i].erase(userParticleStrings[i].begin()
                                  +userParticleStrings[i].find("{", 0));
    while(userParticleStrings[i].find("}", 0) != string::npos)
      userParticleStrings[i].erase(userParticleStrings[i].begin()
                                  +userParticleStrings[i].find("}", 0));
    while(userParticleStrings[i].find(" ", 0) != string::npos)
      userParticleStrings[i].erase(userParticleStrings[i].begin()
                                  +userParticleStrings[i].find(" ", 0));
  }
  // Convert particle numbers in user particle to integers
  vector<int>userParticleNumbers;
  if ( int(userParticleStrings.size()) > 1) {
    for( int i = 1; i < int(userParticleStrings.size()); ++i) {
      userParticleNumbers.push_back( 
        atoi((char*)userParticleStrings[i].substr(
          userParticleStrings[i].find(",",0)+1, 
          userParticleStrings[i].size()).c_str() ) );
    }
  }

  // Save remaining process string
  if (nUserParticles > 0)
    userParticleStrings.push_back(
      fullProc.substr(
        fullProc.find("}",0)+1, fullProc.size() ) );
  // Remove curly brackets and whitespace
  for( int i = 0; i < int(userParticleStrings.size()); ++i) {
    while(userParticleStrings[i].find("{", 0) != string::npos)
      userParticleStrings[i].erase(userParticleStrings[i].begin()
                                  +userParticleStrings[i].find("{", 0));
    while(userParticleStrings[i].find("}", 0) != string::npos)
      userParticleStrings[i].erase(userParticleStrings[i].begin()
                                  +userParticleStrings[i].find("}", 0));
    while(userParticleStrings[i].find(" ", 0) != string::npos)
      userParticleStrings[i].erase(userParticleStrings[i].begin()
                                  +userParticleStrings[i].find(" ", 0));
  }

  // Start mapping residual process string onto particle IDs.
  // Declare leftover process after user-defined particles have been converted
  string residualProc;
  if ( int(userParticleStrings.size()) > 1 )
    residualProc = userParticleStrings.front() + userParticleStrings.back();
  else
    residualProc = fullProc;

  // Remove comma separation
  while(residualProc.find(",", 0) != string::npos)
    residualProc.erase(residualProc.begin()+residualProc.find(",",0));

  // Count number of resonances
  int appearances = 0;
  for(int n = residualProc.find("(", 0); n != int(string::npos);
          n = residualProc.find("(", n)) {
    appearances++;
    n++;
  }

  // Cut string in incoming, resonance+decay and outgoing pieces 
  vector <string> pieces;
  for(int i =0; i < appearances;++i) {
    int n = residualProc.find("(", 0);
    pieces.push_back(residualProc.substr(0,n));
    residualProc = residualProc.substr(n+1,residualProc.size());
  }
  // Cut last resonance from rest
  if (appearances > 0) {
    pieces.push_back( residualProc.substr(0,residualProc.find(")",0)) );
    pieces.push_back( residualProc.substr(
      residualProc.find(")",0)+1, residualProc.size()) );
  }

  // If the string was not cut into pieces, i.e. no resonance was
  // required, cut string using '>' as delimiter
  if (pieces.empty() ){
    appearances = 0;
    for(int n = residualProc.find(">", 0); n != int(string::npos);
            n = residualProc.find(">", n)) {
      appearances++;
      n++;
    }

    // Cut string in incoming and outgoing pieces 
    for(int i =0; i < appearances;++i) {
      int n = residualProc.find(">", 0);
      pieces.push_back(residualProc.substr(0,n));
      residualProc = residualProc.substr(n+1,residualProc.size());
    }

    if (appearances == 1) pieces.push_back(residualProc);
    if (appearances > 1) {
      pieces.push_back( residualProc.substr(0,residualProc.find(">",0)) );
      pieces.push_back( residualProc.substr(
        residualProc.find(">",0)+1, residualProc.size()) );
    }
  }

  // Get incoming particles
  for(int i=0; i < nIn; ++i) {
    for(int n = pieces[0].find(inParticleNamesMG[i], 0);
           n != int(string::npos);
           n = pieces[0].find(inParticleNamesMG[i], n)) {
      incom.push_back(inParticleNumbers[i]);
      pieces[0].erase(pieces[0].begin()+n,
                      pieces[0].begin()+n+inParticleNamesMG[i].size());
      n=0;
    }
  }

  // Check intermediate resonances and decay products
  for(int i =1; i < int(pieces.size()); ++i){
    // Seperate strings into intermediate and outgoing, if not already done
    int k = pieces[i].find(">", 0);

    string intermediate = (pieces[i].find(">", 0) != string::npos) ? 
                           pieces[i].substr(0,k) : "";
    string outgoing = (pieces[i].find(">", 0) != string::npos) ? 
                       pieces[i].substr(k+1,pieces[i].size()) : pieces[i];

    // Get intermediate particles
    for(int j=0; j < nInt; ++j) {
      for(int n = intermediate.find(interParticleNamesMG[j], 0);
             n != int(string::npos);
             n = intermediate.find(interParticleNamesMG[j], n)) {
        inter.push_back(interParticleNumbers[j]);
        intermediate.erase(intermediate.begin()+n,
                    intermediate.begin()+n+interParticleNamesMG[j].size());
        n=0;
      }
    }

    // Get outgoing particles
    for(int j=0; j < nOut; ++j) {
      for(int n = outgoing.find(outParticleNamesMG[j], 0);
             n != int(string::npos);
             n = outgoing.find(outParticleNamesMG[j], n)) {
        outgo.push_back(outParticleNumbers[j]);
        outgoing.erase(outgoing.begin()+n,
                       outgoing.begin()+n+outParticleNamesMG[j].size());
        n=0;
      }
    }

    // For arbitrary or non-existing intermediate, remember zero for each 
    // two outgoing particles, without bosons.
    if (inter.empty()) {

      // For final state bosons, bookkeep the final state boson as
      // intermediate as well.
      int nBosons = 0;
      for(int l=0; l < int(outgo.size()); ++l)
        if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
          nBosons++;

      int nZeros = (outgo.size() - nBosons)/2;
      for(int l=0; l < nZeros; ++l)
        inter.push_back(0);
    }

    // For final state bosons, bookkeep the final state boson as
    // intermediate as well.
    for(int l=0; l < int(outgo.size()); ++l)
      if ( (abs(outgo[l]) > 20 && abs(outgo[l]) <= 25) || outgo[l] == 2400)
        inter.push_back(outgo[l]);

  }

  // Now store incoming, intermediate and outgoing
  // Set intermediate tags
  for(int i=0; i < int(inter.size()); ++i)
    hardIntermediate.push_back(inter[i]);

  // Set the incoming particle tags
  if (incom.size() != 2)
    cout << "Only two incoming particles allowed" << endl;
  else {
    hardIncoming1 = incom[0];
    hardIncoming2 = incom[1];
  }

  // Remember final state bosons
  int nBosons = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( (abs(outgo[i]) > 20 && abs(outgo[i]) <= 25) || outgo[i] == 2400)
      nBosons++;
  // Remember b-quark container
  int nBQuarks = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 5000)
      nBQuarks++;
  // Remember jet container
  int nJets = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 2212)
      nJets++;
  // Remember lepton container
  int nLeptons = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 1100)
      nLeptons++;
  // Remember lepton container
  int nNeutrinos = 0;
  for(int i=0; i < int(outgo.size()); ++i)
    if ( outgo[i] == 1200)
      nNeutrinos++;
  int nContainers = nLeptons + nNeutrinos + nJets + nBQuarks;

  // Set final particle identifiers
  if ( (outgo.size() - nBosons - nContainers)%2 == 1) {
    cout << "Only even number of outgoing particles allowed" << endl;  
    for(int i=0; i < int(outgo.size()); ++i)
      cout << outgo[i] << endl;
  } else {

    // Start with user-defined particles.
    for( int i = 0; i < int(userParticleNumbers.size()); ++i)
      if (userParticleNumbers[i] > 0) {
        hardOutgoing2.push_back( userParticleNumbers[i]);
        hardIntermediate.push_back(0);
        // For non-existing intermediate, remember zero.
      } else if (userParticleNumbers[i] < 0) {
        hardOutgoing1.push_back( userParticleNumbers[i]);
        // For non-existing intermediate, remember zero.
        hardIntermediate.push_back(0);
      }

    // Push back particles / antiparticles
    for(int i=0; i < int(outgo.size()); ++i)
      if (outgo[i] > 0
        && outgo[i] != 2212
        && outgo[i] != 5000
        && outgo[i] != 1100
        && outgo[i] != 1200
        && outgo[i] != 2400
        && outgo[i] != 1000022)
        hardOutgoing2.push_back( outgo[i]);
      else if (outgo[i] < 0)
        hardOutgoing1.push_back( outgo[i]);

    // Save final state W-boson container as particle
    for(int i=0; i < int(outgo.size()); ++i)
      if ( outgo[i] == 2400)
        hardOutgoing2.push_back( outgo[i]);

    // Push back jets, distribute evenly among particles / antiparticles
    // Push back majorana particles, distribute evenly
    int iNow = 0;
    for(int i=0; i < int(outgo.size()); ++i)
      if ( (outgo[i] == 2212
        || outgo[i] == 5000
        || outgo[i] == 1200
        || outgo[i] == 1000022)
        && iNow%2 == 0 ){
        hardOutgoing2.push_back( outgo[i]);
        iNow++;
      } else if ( (outgo[i] == 2212
               || outgo[i] == 5000
               || outgo[i] == 1100
               || outgo[i] == 1000022)
               && iNow%2 == 1 ){
        hardOutgoing1.push_back( outgo[i]);
        iNow++;
      }
  }

  // Done
}

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

// Function to check if the candidates stored in Pos1 and Pos2, together with
// a proposed candidate iPos are allowed.

bool HardProcess::allowCandidates(int iPos, vector<int> Pos1, 
  vector<int> Pos2, const Event& event){

  bool allowed = true;

  // Find colour-partner of new candidate
  int type = (event[iPos].col() > 0) ? 1 : (event[iPos].acol() > 0) ? -1 : 0;

  if (type == 0) return true;

  if (type == 1){
    int col = event[iPos].col();
    int iPartner = 0;
    for(int i=0; i < int(event.size()); ++i)
      if ( i != iPos
        && (( event[i].isFinal() && event[i].acol() == col)
          ||( event[i].status() == -21 && event[i].col() == col) ))
      iPartner = i;

    vector<int> partners;
    for(int i=0; i < int(event.size()); ++i)
      for(int j=0; j < int(Pos1.size()); ++j)
        if ( Pos1[j] != 0 && i != Pos1[j] && event[Pos1[j]].colType() != 0
        && (( event[i].isFinal() && event[i].col() == event[Pos1[j]].acol())
          ||( event[i].status() == -21
           && event[i].acol() == event[Pos1[j]].acol()) ))
         partners.push_back(i);

    // Never allow equal initial partners!
    if (event[iPartner].status() == -21){
      for(int i=0; i < int(partners.size()); ++i)
        if ( partners[i] == iPartner)
          allowed = false;
    }

  } else {
    int col = event[iPos].acol();
    int iPartner = 0;
    for(int i=0; i < int(event.size()); ++i)
      if ( i != iPos
        && (( event[i].isFinal() && event[i].col()  == col)
          ||(!event[i].isFinal() && event[i].acol() == col) ))
      iPartner = i;

    vector<int> partners;
    for(int i=0; i < int(event.size()); ++i)
      for(int j=0; j < int(Pos2.size()); ++j)
        if ( Pos2[j] != 0 && i != Pos2[j] && event[Pos2[j]].colType() != 0
        && (( event[i].isFinal() && event[i].acol() == event[Pos2[j]].col())
          ||( event[i].status() == -21
           && event[i].col() == event[Pos2[j]].col()) ))
         partners.push_back(i);


    // Never allow equal initial partners!
    if (event[iPartner].status() == -21){
      for(int i=0; i < int(partners.size()); ++i){
        if ( partners[i] == iPartner)
          allowed = false;
      }
    }

  }

  return allowed;

}

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

// Function to identify the hard subprocess in the current event

void HardProcess::storeCandidates( const Event& event, string process){

  // Store the reference event
  state.clear();
  state = event;

  // Local copy of intermediate bosons
  vector<int> intermediates;
  for(int i =0; i < int(hardIntermediate.size());++i)
    intermediates.push_back( hardIntermediate[i]);

  // Local copy of outpoing partons
  vector<int> outgoing1;
  for(int i =0; i < int(hardOutgoing1.size());++i)
    outgoing1.push_back( hardOutgoing1[i]);
  vector<int> outgoing2;
  for(int i =0; i < int(hardOutgoing2.size());++i)
    outgoing2.push_back( hardOutgoing2[i]);

  // Clear positions of intermediate and outgoing particles
  PosIntermediate.resize(0);
  PosOutgoing1.resize(0);
  PosOutgoing2.resize(0);
  for(int i =0; i < int(hardIntermediate.size());++i)
    PosIntermediate.push_back(0);
  for(int i =0; i < int(hardOutgoing1.size());++i)
    PosOutgoing1.push_back(0);
  for(int i =0; i < int(hardOutgoing2.size());++i)
    PosOutgoing2.push_back(0);

  // For QCD dijet or e+e- > jets hard process, do not store any candidates,
  // as to not discrimintate clusterings
  if (  process.compare("pp>jj") == 0
    || process.compare("e+e->jj") == 0
    || process.compare("e+e->(z>jj)") == 0 ){
    for(int i =0; i < int(hardOutgoing1.size());++i)
      PosOutgoing1[i] = 0;
    for(int i =0; i < int(hardOutgoing2.size());++i)
      PosOutgoing2[i] = 0;
    // Done
    return;
  }

  // Initialise vector of particles that were already identified as
  // hard process particles
  vector<int> iPosChecked;

  // If the hard process is specified only by containers, then add all
  // particles matching with the containers to the hard process.
  bool hasOnlyContainers = true;
  for(int i =0; i < int(hardOutgoing1.size());++i)
    if (  hardOutgoing1[i] != 1100
      && hardOutgoing1[i] != 1200
      && hardOutgoing1[i] != 5000)
      hasOnlyContainers = false;
  for(int i =0; i < int(hardOutgoing2.size());++i)
    if (  hardOutgoing2[i] != 1100
      && hardOutgoing2[i] != 1200
      && hardOutgoing2[i] != 5000)
      hasOnlyContainers = false;

  if (hasOnlyContainers){

    PosOutgoing1.resize(0);
    PosOutgoing2.resize(0);

    // Try to find all unmatched hard process leptons.
    // Loop through event to find outgoing lepton
    for(int i=0; i < int(event.size()); ++i){

      // Skip non-final particles
      if ( !event[i].isFinal() ) continue;

      // Skip all particles that have already been identified
      bool skip = false;
      for(int k=0; k < int(iPosChecked.size()); ++k){
        if (i == iPosChecked[k])
          skip = true;
      }
      if (skip) continue;

      for(int j=0; j < int(outgoing2.size()); ++j){

        // If the particle matches an outgoing neutrino, save it
        if ( outgoing2[j] == 1100
          && ( event[i].idAbs() == 11
            || event[i].idAbs() == 13
            || event[i].idAbs() == 15) ){
          PosOutgoing2.push_back(i);
          iPosChecked.push_back(i);
        }

        // If the particle matches an outgoing lepton, save it
        if ( outgoing2[j] == 1200
          && ( event[i].idAbs() == 12
            || event[i].idAbs() == 14
            || event[i].idAbs() == 16) ){
          PosOutgoing2.push_back(i);
          iPosChecked.push_back(i);
        }

        // If the particle matches an outgoing b-quark, save it
        if ( outgoing2[j] == 5000 && event[i].idAbs() == 5 ){
          PosOutgoing2.push_back(i);
          iPosChecked.push_back(i);
        }

      }

      // Skip all particles that have already been identified
      skip = false;
      for(int k=0; k < int(iPosChecked.size()); ++k){
        if (i == iPosChecked[k])
          skip = true;
      }
      if (skip) continue;

      for(int j=0; j < int(outgoing1.size()); ++j){

        // If the particle matches an outgoing neutrino, save it
        if ( outgoing1[j] == 1100
          && ( event[i].idAbs() == 11
            || event[i].idAbs() == 13
            || event[i].idAbs() == 15) ){
          PosOutgoing1.push_back(i);
          iPosChecked.push_back(i);
        }

        // If the particle matches an outgoing lepton, save it
        if ( outgoing1[j] == 1200
          && ( event[i].idAbs() == 12
            || event[i].idAbs() == 14
            || event[i].idAbs() == 16) ){
          PosOutgoing1.push_back(i);
          iPosChecked.push_back(i);
        }

        // If the particle matches an outgoing b-quark, save it
        if ( outgoing1[j] == 5000 && event[i].idAbs() == 5 ){
          PosOutgoing1.push_back(i);
          iPosChecked.push_back(i);
        }

      }
    }

    // Done
    return;
  }

  // Now begin finding candidates when not only containers are used.

  // First try to find final state bosons
  for(int i=0; i < int(intermediates.size()); ++i){

    // Do nothing if the intermediate boson is absent
    if (intermediates[i] == 0) continue;

    // Do nothing if this boson does not match any final state boson
    bool matchesFinalBoson = false;
    for(int j =0; j< int(outgoing1.size()); ++j){
      if ( intermediates[i] == outgoing1[j] )
        matchesFinalBoson = true;
    }
    for(int j =0; j< int(outgoing2.size()); ++j){
      if ( intermediates[i] == outgoing2[j] )
        matchesFinalBoson = true;
    }
    if (!matchesFinalBoson) continue;

    // Loop through event
    for(int j=0; j < int(event.size()); ++j) {
      // If the particle has a requested intermediate id, check if
      // if is a final state boson
      if ( (event[j].id() == intermediates[i])
        ||(event[j].idAbs() == 24 && intermediates[i] == 2400) ) {
        PosIntermediate[i] = j;
        intermediates[i] = 0;

        for(int k=0; k < int(outgoing1.size()); ++k) {
          if (event[j].id() == outgoing1[k]){
            PosOutgoing1[k] = j;
            outgoing1[k] = 99;
          }
        }

        for(int k=0; k < int(outgoing2.size()); ++k) {
          if (event[j].id() == outgoing2[k]){
            PosOutgoing2[k] = j;
            outgoing2[k] = 99;
          }
        }

        // Check for W-boson container
        for(int k=0; k < int(outgoing2.size()); ++k) {
          if (event[j].idAbs() == 24 && outgoing2[k] == 2400){
            PosOutgoing2[k] = j;
            outgoing2[k] = 99;
          }
        }

        iPosChecked.push_back(j);

      }
    }
  }

  // Second try to find particles coupled to intermediate bosons
  for(int i=0; i < int(intermediates.size()); ++i){

    // Do nothing if the intermediate boson is absent
    if (intermediates[i] == 0) continue;

    // Loop through event
    for(int j=0; j < int(event.size()); ++j) {
      // If the particle has a requested intermediate id, check if
      // daughters are hard process particles
      if ( (event[j].id() == intermediates[i])
        ||(event[j].idAbs() == 24 && intermediates[i] == 2400) ) {
        // If this particle is a potential intermediate
        PosIntermediate[i] = j;
        intermediates[i] = 0;
        // If id's of daughters are good, store position
        int iPos1 = event[j].daughter1();
        int iPos2 = event[j].daughter2();

        // Loop through daughters to check if these contain some hard
        // outgoing particles
        for( int k=iPos1; k <= iPos2; ++k){
          int id = event[k].id();

          // Check if daughter is hard outgoing particle
          for(int l=0; l < int(outgoing2.size()); ++l)
            if ( outgoing2[l] != 99 ){
                // Found particle id
              if (id == outgoing2[l]
                // Found jet
                || (id > 0 && abs(id) < 10 && outgoing2[l] == 2212) ){
                // Store position
                PosOutgoing2[l] = k;
                // Remove the matched particle from the list
                outgoing2[l] = 99;
                iPosChecked.push_back(k);
                break;
              }

            }
 
          // Check if daughter is hard outgoing antiparticle
          for(int l=0; l < int(outgoing1.size()); ++l)
            if ( outgoing1[l] != 99 ){
                // Found particle id
              if (id == outgoing1[l]
                // Found jet
                || (id < 0 && abs(id) < 10 && outgoing1[l] == 2212) ){
                // Store position
                PosOutgoing1[l] = k;
                // Remove the matched particle from the list
                outgoing1[l] = 99;
                iPosChecked.push_back(k);
                break;
            }

          }

        } // End loop through daughters
      } // End if ids match
    } // End loop through event
  } // End loop though requested intermediates 

  // If all outgoing particles were found, done
  bool done = true;
  for(int i=0; i < int(outgoing1.size()); ++i)
    if (outgoing1[i] != 99)
      done = false;
  for(int i=0; i < int(outgoing2.size()); ++i)
    if (outgoing2[i] != 99)
      done = false;
  // Return 
  if (done) return;

  // Leptons not associated with resonance are allowed.
  // Try to find all unmatched hard process leptons.
  // Loop through event to find outgoing lepton
  for(int i=0; i < int(event.size()); ++i){
    // Skip non-final particles and final partons
    if ( !event[i].isFinal() || event[i].colType() != 0)
      continue;
    // Skip all particles that have already been identified
    bool skip = false;
    for(int k=0; k < int(iPosChecked.size()); ++k){
      if (i == iPosChecked[k])
        skip = true;
    }
    if (skip) continue;

    // Check if any hard outgoing leptons remain
    for(int j=0; j < int(outgoing2.size()); ++j){
      // Do nothing if this particle has already be found,
      // or if this particle is a jet or quark
      if (  outgoing2[j] == 99
        || outgoing2[j] == 2212
        || abs(outgoing2[j]) < 10)
        continue;

      // If the particle matches an outgoing lepton, save it
      if (  event[i].id() == outgoing2[j] ){
        PosOutgoing2[j] = i;
        outgoing2[j] = 99;
        iPosChecked.push_back(i);
      }
    }

    // Check if any hard outgoing antileptons remain
    for(int j=0; j < int(outgoing1.size()); ++j){
      // Do nothing if this particle has already be found,
      // or if this particle is a jet or quark
      if (  outgoing1[j] == 99
        || outgoing1[j] == 2212
        || abs(outgoing1[j]) < 10)
        continue;

      // If the particle matches an outgoing lepton, save it
      if (event[i].id() == outgoing1[j] ){
        PosOutgoing1[j] = i;
        outgoing1[j] = 99;
        iPosChecked.push_back(i);
      }
    }
  }

  multimap<int,int> out2copy;
  for(int i=0; i < int(event.size()); ++i)
    for(int j=0; j < int(outgoing2.size()); ++j)
      // Do nothing if this particle has already be found,
      // or if this particle is a jet, lepton container or lepton
      if (  outgoing2[j] != 99
        && outgoing2[j] != 2212
        && abs(outgoing2[j]) < 10
        && event[i].isFinal()
        && event[i].id() == outgoing2[j] ){
        out2copy.insert(make_pair(j, i));
      }

  multimap<int,int> out1copy;
  for(int i=0; i < int(event.size()); ++i)
    for(int j=0; j < int(outgoing1.size()); ++j)
      // Do nothing if this particle has already be found,
      // or if this particle is a jet, lepton container or lepton
      if (  outgoing1[j] != 99
        && outgoing1[j] != 2212
        && abs(outgoing1[j]) < 10
        && event[i].isFinal()
        && event[i].id() == outgoing1[j] ){
        out1copy.insert(make_pair(j, i));
      }

  if ( out1copy.size() >  out2copy.size()){

    for ( multimap<int, int>::iterator it = out2copy.begin(); 
      it != out2copy.end(); ++it ) {
      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){

        // Save parton
        PosOutgoing2[it->first] = it->second;
        // remove entry form lists
        outgoing2[it->first] = 99;
        iPosChecked.push_back(it->second);

      }
    }

    for ( multimap<int, int>::iterator it = out1copy.begin(); 
      it != out1copy.end(); ++it ) {
      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){

        // Save parton
        PosOutgoing1[it->first] = it->second;
        // remove entry form lists
        outgoing1[it->first] = 99;
        iPosChecked.push_back(it->second);

      }
    }

  } else {

    for ( multimap<int, int>::iterator it = out1copy.begin(); 
      it != out1copy.end(); ++it ) {
      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){

        // Save parton
        PosOutgoing1[it->first] = it->second;
        // remove entry form lists
        outgoing1[it->first] = 99;
        iPosChecked.push_back(it->second);

      }
    }

    for ( multimap<int, int>::iterator it = out2copy.begin(); 
      it != out2copy.end(); ++it ) {
      if ( allowCandidates(it->second, PosOutgoing1, PosOutgoing2, event) ){

        // Save parton
        PosOutgoing2[it->first] = it->second;
        // remove entry form lists
        outgoing2[it->first] = 99;
        iPosChecked.push_back(it->second);

      }
    }
  }

  // It sometimes happens that MadEvent does not put a
  // heavy coloured resonance into the LHE file, even if requested.
  // This means that the decay products of this resonance need to be
  // found separately.
  // Loop through event to find hard process (anti)quarks
  for(int i=0; i < int(event.size()); ++i){

    // Skip non-final particles and final partons
    if ( !event[i].isFinal() || event[i].colType() == 0)
      continue;

    // Skip all particles that have already been identified
    bool skip = false;
    for(int k=0; k < int(iPosChecked.size()); ++k){
      if (i == iPosChecked[k])
        skip = true;
    }
    if (skip) continue;

    // Check if any hard outgoing quarks remain
    for(int j=0; j < int(outgoing2.size()); ++j){
      // Do nothing if this particle has already be found,
      // or if this particle is a jet, lepton container or lepton
      if (  outgoing2[j] == 99
        || outgoing2[j] == 2212
        || abs(outgoing2[j]) > 10)
        continue;
      // If the particle matches an outgoing quark, save it
      if (event[i].id() == outgoing2[j]){
        // Save parton
        PosOutgoing2[j] = i;
        // remove entry form lists
        outgoing2[j] = 99;
        iPosChecked.push_back(i);
      }
    }

    // Check if any hard outgoing antiquarks remain
    for(int j=0; j < int(outgoing1.size()); ++j){
      // Do nothing if this particle has already be found,
      // or if this particle is a jet, lepton container or lepton
      if (  outgoing1[j] == 99
        || outgoing1[j] == 2212
        || abs(outgoing1[j]) > 10)
        continue;
      // If the particle matches an outgoing antiquark, save it
      if (event[i].id() == outgoing1[j]){
        // Save parton
        PosOutgoing1[j] = i;
        // Remove parton from list
        outgoing1[j] = 99;
        iPosChecked.push_back(i);
      }
    }
  }

  // Done
}

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

// Function to check if the particle event[iPos] matches any of
// the stored outgoing particles of the hard subprocess

bool HardProcess::matchesAnyOutgoing(int iPos, const Event& event){

  // Match quantum numbers of any first outgoing candidate
  bool matchQN1 = false;
  // Match quantum numbers of any second outgoing candidate
  bool matchQN2 = false;
  // Match parton in the hard process,
  // or parton from decay of electroweak boson in hard process,
  // or parton from decay of electroweak boson from decay of top
  bool matchHP = false;

  // Check outgoing candidates
  for(int i=0; i < int(PosOutgoing1.size()); ++i)
    // Compare particle properties
    if ( event[iPos].id()         == state[PosOutgoing1[i]].id()
     && event[iPos].colType()    == state[PosOutgoing1[i]].colType() 
     && event[iPos].chargeType() == state[PosOutgoing1[i]].chargeType() 
     && event[iPos].col()        == state[PosOutgoing1[i]].col() 
     && event[iPos].acol()       == state[PosOutgoing1[i]].acol()
     && event[iPos].charge()     == state[PosOutgoing1[i]].charge() )
      matchQN1 = true;

  // Check outgoing candidates
  for(int i=0; i < int(PosOutgoing2.size()); ++i)
    // Compare particle properties
    if ( event[iPos].id()         == state[PosOutgoing2[i]].id()
     && event[iPos].colType()    == state[PosOutgoing2[i]].colType() 
     && event[iPos].chargeType() == state[PosOutgoing2[i]].chargeType() 
     && event[iPos].col()        == state[PosOutgoing2[i]].col() 
     && event[iPos].acol()       == state[PosOutgoing2[i]].acol()
     && event[iPos].charge()     == state[PosOutgoing2[i]].charge() )
      matchQN2 = true;

  // Check if maps to hard process:
  // Check that particle is in hard process
  if ( event[iPos].mother1()*event[iPos].mother2() == 12
      // Or particle has taken recoil from first splitting
      || (  event[iPos].status() == 44
         && event[event[iPos].mother1()].mother1()
           *event[event[iPos].mother1()].mother2() == 12 )
      // Or particle has on-shell resonace as mother
      || (  event[iPos].status() == 23
         && event[event[iPos].mother1()].mother1()
           *event[event[iPos].mother1()].mother2() == 12 )
      // Or particle has on-shell resonace as mother,
      // which again has and on-shell resonance as mother
      || (  event[iPos].status() == 23
         && event[event[iPos].mother1()].status() == -22
         && event[event[event[iPos].mother1()].mother1()].status() == -22
         && event[event[event[iPos].mother1()].mother1()].mother1()
           *event[event[event[iPos].mother1()].mother1()].mother2() == 12 ) )
      matchHP = true;

  // Done
  return ( matchHP && (matchQN1 || matchQN2) );

}


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

// Function to return the type of the ME generator

int HardProcess::MEgenType(){ return meGenType;}

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

// Function to get the number of coloured final state partons in the
// hard process

int HardProcess::nQuarksOut(){
  int nFin =0;
  for(int i =0; i< int(hardOutgoing1.size()); ++i){
    if (hardOutgoing1[i] == 2212 || abs(hardOutgoing1[i]) < 10) nFin++;
  }
  for(int i =0; i< int(hardOutgoing2.size()); ++i){
    if (hardOutgoing2[i] == 2212 || abs(hardOutgoing2[i]) < 10) nFin++;
  }
  // For very loose hard process definition, check number of hard process
  // b-quarks explicitly.
  for(int i =0; i< int(hardOutgoing1.size()); ++i)
    if (hardOutgoing1[i] == 5000)
      for(int j =0; j< int(PosOutgoing1.size()); ++j)
        if (state[PosOutgoing1[j]].idAbs() == 5)
          nFin++;
  for(int i =0; i< int(hardOutgoing2.size()); ++i)
    if (hardOutgoing2[i] == 5000)
      for(int j =0; j< int(PosOutgoing2.size()); ++j)
        if (state[PosOutgoing2[j]].idAbs() == 5)
          nFin++;
  return nFin;
}

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

// Function to get the number of uncoloured final state particles in the
// hard process

int HardProcess::nLeptonOut(){
  int nFin =0;
  for(int i =0; i< int(hardOutgoing1.size()); ++i){
    if (abs(hardOutgoing1[i]) > 10 && abs(hardOutgoing1[i]) < 20) nFin++;
    // Bookkeep MSSM neutralinos as leptons
    if (abs(hardOutgoing1[i]) == 1000022) nFin++;
  }
  for(int i =0; i< int(hardOutgoing2.size()); ++i){
    if (abs(hardOutgoing2[i]) > 10 && abs(hardOutgoing2[i]) < 20) nFin++;
    // Bookkeep MSSM neutralinos as leptons
    if (abs(hardOutgoing2[i]) == 1000022) nFin++;
  }
  // For very loose hard process definition, check number of hard process
  // lepton explicitly.
  // Check lepton / neutrino containers as leptons
  for(int i =0; i< int(hardOutgoing1.size()); ++i)
    if (hardOutgoing1[i] == 1100)
      for(int j =0; j< int(PosOutgoing1.size()); ++j)
        if (  abs(state[PosOutgoing1[j]].id()) == 11
          || abs(state[PosOutgoing1[j]].id()) == 13
          || abs(state[PosOutgoing1[j]].id()) == 15 )
          nFin++;
  for(int i =0; i< int(hardOutgoing2.size()); ++i)
    if (hardOutgoing2[i] == 1200)
      for(int j =0; j< int(PosOutgoing2.size()); ++j)
        if (  abs(state[PosOutgoing2[j]].id()) == 12
          || abs(state[PosOutgoing2[j]].id()) == 14
          || abs(state[PosOutgoing2[j]].id()) == 16 )
          nFin++;
  return nFin;
}

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

// Function to get the number of uncoloured final state particles in the
// hard process

int HardProcess::nBosonsOut(){
  int nFin =0;
  for(int i =0; i< int(hardOutgoing1.size()); ++i){
    if (abs(hardOutgoing1[i]) > 20 && abs(hardOutgoing1[i]) <= 25) nFin++;
  }
  for(int i =0; i< int(hardOutgoing2.size()); ++i){
    if (abs(hardOutgoing2[i]) > 20 && abs(hardOutgoing2[i]) <= 25) nFin++;
    if ( hardOutgoing2[i] == 2400) nFin++;
  }
  return nFin;
}

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

// Function to get the number of coloured initial state partons in the
// hard process

int HardProcess::nQuarksIn(){
  int nIn =0;
  if (hardIncoming1 == 2212 || abs(hardIncoming1) < 10) nIn++;
  if (hardIncoming2 == 2212 || abs(hardIncoming2) < 10) nIn++;
  return nIn;
}

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

// Function to get the number of uncoloured initial state particles in the
// hard process

int HardProcess::nLeptonIn(){
  int nIn =0;
  if (abs(hardIncoming1) > 10 && abs(hardIncoming1) < 20) nIn++;
  if (abs(hardIncoming2) > 10 && abs(hardIncoming2) < 20) nIn++;
  return nIn;
}

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

// Function to report if a resonace decay was found in the
// 2->2 hard sub-process in the current state

int HardProcess::hasResInCurrent(){
  for(int i =0; i< int(PosIntermediate.size()); ++i)
    if (PosIntermediate[i] == 0) return 0;
  // Do not count final state bosons as resonaces
  for(int i =0; i< int(PosIntermediate.size()); ++i){
    for(int j =0; j< int(PosOutgoing1.size()); ++j){
      if ( PosIntermediate[i] == PosOutgoing1[j] )
        return 0;
    }
    for(int j =0; j< int(PosOutgoing2.size()); ++j){
      if ( PosIntermediate[i] == PosOutgoing2[j] )
        return 0;
    }
  }
  return 1;
}

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

// Function to report the number of resonace decays in the 2->2 sub-process 
// of the  current state

int HardProcess::nResInCurrent(){
  int nRes = 0;
  for(int i =0; i< int(PosIntermediate.size()); ++i){
    if (PosIntermediate[i] != 0) {
      bool matchesFinalBoson = false;
      for(int j =0; j< int(PosOutgoing1.size()); ++j){
        if ( PosIntermediate[i] == PosOutgoing1[j] )
          matchesFinalBoson = true;
      }
      for(int j =0; j< int(PosOutgoing2.size()); ++j){
        if ( PosIntermediate[i] == PosOutgoing2[j] )
          matchesFinalBoson = true;
      }
      if (!matchesFinalBoson) nRes++;
    }
  }
  return nRes;
}

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

// Function to report if a resonace decay was found in the
// 2->2 hard core process

bool HardProcess::hasResInProc(){

  for(int i =0; i< int(hardIntermediate.size()); ++i)
    if (hardIntermediate[i] == 0) return false;
  // Do not count final state bosons as resonaces
  for(int i =0; i< int(hardIntermediate.size()); ++i){
    for(int j =0; j< int(hardOutgoing1.size()); ++j){
      if ( hardIntermediate[i] == hardOutgoing1[j] )
        return false;
    }
    for(int j =0; j< int(hardOutgoing2.size()); ++j){
      if ( hardIntermediate[i] == hardOutgoing2[j] )
        return false;
    }
  }
  return true;
}

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

// Function to print the hard process (for debug)

void HardProcess::list() const {
  cout << "   Hard Process: ";
  cout << " \t " << hardIncoming1 << " + " << hardIncoming2;
  cout << " \t -----> \t ";
  for(int i =0; i < int(hardIntermediate.size());++i)
    cout << hardIntermediate[i] << " "; 
  cout << " \t -----> \t ";
  for(int i =0; i < int(hardOutgoing1.size());++i)
    cout << hardOutgoing1[i] << " ";
  for(int i =0; i < int(hardOutgoing2.size());++i)
    cout << hardOutgoing2[i] << " ";
  cout << endl;
}

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

// Function to list the hard process candiates in the matrix element state
// (for debug)

void HardProcess::listCandidates() const {
  cout << "   Hard Process candidates: ";
  cout << " \t " << hardIncoming1 << " + " << hardIncoming2;
  cout << " \t -----> \t ";
  for(int i =0; i < int(PosIntermediate.size());++i)
    cout << PosIntermediate[i] << " "; 
  cout << " \t -----> \t ";
  for(int i =0; i < int(PosOutgoing1.size());++i)
    cout << PosOutgoing1[i] << " ";
  for(int i =0; i < int(PosOutgoing2.size());++i)
    cout << PosOutgoing2[i] << " ";
  cout << endl;
}

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

// Function to clear hard process information

void HardProcess::clear() {

  // Clear flavour of the first incoming particle
  hardIncoming1 = hardIncoming2 = 0;
  // Clear outgoing particles
  hardOutgoing1.resize(0);
  hardOutgoing2.resize(0);
  // Clear intermediate bosons in the hard 2->2
  hardIntermediate.resize(0);

  // Clear reference event
  state.clear();

  // Clear potential positions of outgoing particles in reference event
  PosOutgoing1.resize(0);
  PosOutgoing2.resize(0);
  // Clear potential positions of intermediate bosons in reference event
  PosIntermediate.resize(0);

  // Clear merging scale read from LHE file
  tms = 0.;
  // Clear type of ME generator
  meGenType = 0;


}

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

// The MergingHooks class.

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

// Initialise MergingHooks class

void MergingHooks::init( Settings& settings, Info* infoPtrIn, 
  ParticleData* particleDataPtrIn, ostream& os){

  // Save pointers
  infoPtr          = infoPtrIn;
  particleDataPtr  = particleDataPtrIn;

  // Initialise AlphaS objects for reweighting
  double alphaSvalueFSR = settings.parm("TimeShower:alphaSvalue");
  int    alphaSorderFSR = settings.mode("TimeShower:alphaSorder");
  AlphaS_FSRSave.init(alphaSvalueFSR,alphaSorderFSR);
  double alphaSvalueISR = settings.parm("SpaceShower:alphaSvalue");
  int    alphaSorderISR = settings.mode("SpaceShower:alphaSorder");
  AlphaS_ISRSave.init(alphaSvalueISR,alphaSorderISR);

  // Initialise AlphaS objects for reweighting
  int    alphaEMFSRorder = settings.mode("TimeShower:alphaEMorder");
  AlphaEM_FSRSave.init(alphaEMFSRorder, &settings);

  // Initialise merging switches
  doUserMergingSave     =  settings.flag("Merging:doUserMerging");

  // Initialise automated MadGraph kT merging
  doMGMergingSave       =  settings.flag("Merging:doMGMerging");

  // Initialise kT merging
  doKTMergingSave       =  settings.flag("Merging:doKTMerging");

  // Initialise evolution-pT merging
  doPTLundMergingSave   =  settings.flag("Merging:doPTLundMerging");

  // Initialise \Delta_R_{ij}, pT_i Q_{ij} merging
  doCutBasedMergingSave =  settings.flag("Merging:doCutBasedMerging");

  // Initialise exact definition of kT
  ktTypeSave            =  settings.mode("Merging:ktType");

  // Get core process from user input
  processSave = settings.word("Merging:Process");

  // Clear hard process
  hardProcess.clear();

  bool doStandardMerging = doUserMergingSave || doKTMergingSave 
         || doPTLundMergingSave || doCutBasedMergingSave;

  // Initialise the hard process
  if ( doStandardMerging )
    hardProcess.initOnProcess(processSave, particleDataPtr);
  else
    hardProcess.initOnLHEF(lheInputFile, particleDataPtr);

  // Parameters for reconstruction of evolution scales
  includeMassiveSave        = settings.flag("Merging:includeMassive");
  enforceStrongOrderingSave = settings.flag("Merging:enforceStrongOrdering");
  scaleSeparationFactorSave = settings.parm("Merging:scaleSeparationFactor");
  orderInRapiditySave       = settings.flag("Merging:orderInRapidity");

  // Parameters for choosing history probabilistically
  nonJoinedNormSave     = settings.parm("Merging:nonJoinedNorm");
  fsrInRecNormSave      = settings.parm("Merging:fsrInRecNorm");
  pickByFullPSave       = settings.flag("Merging:pickByFullP");
  pickByPoPT2Save       = settings.flag("Merging:pickByPoPT2");
  includeRedundantSave  = settings.flag("Merging:includeRedundant");

  // Parameters for scale choices
  unorderedScalePrescipSave   =
    settings.mode("Merging:unorderedScalePrescrip");
  unorderedASscalePrescipSave =
    settings.mode("Merging:unorderedASscalePrescrip");
  unorderedPDFscalePrescipSave =
    settings.mode("Merging:unorderedPDFscalePrescrip");
  incompleteScalePrescipSave  =
    settings.mode("Merging:incompleteScalePrescrip");

  // Parameter for allowing swapping of one colour index while reclustering
  allowColourShufflingSave  =
    settings.flag("Merging:allowColourShuffling");

  // Parameters to allow setting hard process scales to default (dynamical)
  // Pythia values.
  resetHardQRenSave     =  settings.flag("Merging:usePythiaQRenHard");
  resetHardQFacSave     =  settings.flag("Merging:usePythiaQFacHard");

  // Parameters for choosing history by sum(|pT|)
  pickBySumPTSave       = settings.flag("Merging:pickBySumPT");
  herwigAcollFSRSave    = settings.parm("Merging:aCollFSR");
  herwigAcollISRSave    = settings.parm("Merging:aCollISR");

  // Information on the shower cut-off scale
  pT0ISRSave            = settings.parm("SpaceShower:pT0Ref");
  pTcutSave             = settings.parm("SpaceShower:pTmin");
  pTcutSave             = max(pTcutSave,pT0ISRSave);

  // Initialise CKKWL weight
  weightSave = 1.;
  weightCKKWLSave = 1.;

  // Initialise merging scale
  tmsValueSave = 0.;
  tmsListSave.resize(0);

  // Save merging scale on maximal number of jets
  if (  doKTMergingSave || doUserMergingSave || doPTLundMergingSave ) {
    // Read merging scale (defined in kT) from input parameter.
    tmsValueSave    = settings.parm("Merging:TMS");
    nJetMaxSave     = settings.mode("Merging:nJetMax");
  } else if (doMGMergingSave) {
    // Read merging scale (defined in kT) from LHE file.
    tmsValueSave    = hardProcess.tms;
    nJetMaxSave     = settings.mode("Merging:nJetMax");
  } else if (doCutBasedMergingSave) { 
    // Save list of cuts defining the merging scale.
    nJetMaxSave     = settings.mode("Merging:nJetMax");
    // Write tms cut values to list of cut values,
    // ordered by DeltaR_{ij}, pT_{i}, Q_{ij}.
    tmsListSave.resize(0);
    double drms     = settings.parm("Merging:dRijMS");
    double ptms     = settings.parm("Merging:pTiMS");
    double qms      = settings.parm("Merging:QijMS");
    tmsListSave.push_back(drms);
    tmsListSave.push_back(ptms);
    tmsListSave.push_back(qms);
  }

  bool writeBanner =  doKTMergingSave || doMGMergingSave || doUserMergingSave
                   || doPTLundMergingSave || doCutBasedMergingSave;

  if (!writeBanner) return;

  // Write banner
  os << "\n"
     << " *---------- MEPS Merging Initialization  -----------------------*"
     << "\n";
  if ( doKTMergingSave || doMGMergingSave || doUserMergingSave
    || doPTLundMergingSave || doCutBasedMergingSave )
     os << " | CKKW-L tree-level merging:\n"
        << " | We merge  " << processSave << "  with up to  " << nJetMaxSave
        << " additional jets \n";
  if ( doKTMergingSave )
    os << " | Merging scale is defined in kT with the value ktMS = "
       << tmsValueSave << " GeV";
  else if ( doMGMergingSave )
    os << " | Perform automanted MG/ME merging \n"
       << " | Merging scale is defined in kT with the value ktMS = "
       << tmsValueSave << " GeV";
  else if ( doUserMergingSave )
    os << " | Merging scale is defined by the user, with the value tMS = "
       << tmsValueSave;
  else if ( doPTLundMergingSave )
    os << " | Merging scale is defined by Lund pT, with the value tMS = "
       << tmsValueSave;
  else if ( doCutBasedMergingSave )
    os << " | Merging scale is defined by combination of Delta R_{ij}, pT_i\n"
       << " | and Q_{ij} cut, with values \n"
       << " | Delta R_{ij,min} = " << tmsListSave[0] << "\n"
       << " | pT_{i,min} = " << tmsListSave[1] << "\n"
       << " | Q_{ij,min} = " << tmsListSave[2];
  os << "\n *---------- END MEPS Merging Initialization  -------------------*"
     << "\n\n";

}

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

// Function to return the number of clustering steps for the current event

int MergingHooks::getNumberOfClusteringSteps(const Event& event){

  // Count the number of final state partons
  int nFinalPartons = 0;
  for( int i=0; i < event.size(); ++i)
    if ( event[i].isFinal() && (event[i].isQuark() || event[i].isGluon()) )
      nFinalPartons++;

  // Count the number of final state leptons
  int nFinalLeptons = 0;
  for( int i=0; i < event.size(); ++i)
    if ( event[i].isFinal() && event[i].isLepton())
      nFinalLeptons++;

  // Add neutralinos to number of leptons
  for( int i=0; i < event.size(); ++i)
    if ( event[i].isFinal()
       && event[i].idAbs() == 1000022)
      nFinalLeptons++;

  // Count the number of final state electroweak bosons
  int nFinalBosons = 0;
  for( int i=0; i < event.size(); ++i)
    if ( event[i].isFinal()
      && ( event[i].idAbs() == 22
        || event[i].idAbs() == 23
        || event[i].idAbs() == 24
        || event[i].idAbs() == 25 ) )
      nFinalBosons++;

  // Save sum of all final state particles
  int nFinal = nFinalPartons + nFinalLeptons
             + 2*(nFinalBosons - nHardOutBosons() );

  // Return the difference to the core process outgoing particles
  return (nFinal - nHardOutPartons() - nHardOutLeptons() );
}

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

// Function to check if event contains an emission not present in the hard
// process.

bool MergingHooks::isFirstEmission(const Event& event ) {

  // Check that only one additional parton has been produced.
  // If not, we're already in the PS region (e.g. in MI).
  // Then, do not veto.
  int nMPI = infoPtr->nMPI();
  if (nMPI > 1) return false;

  // Count particle types
  int nFinalQuarks   = 0;
  int nFinalGluons   = 0;
  int nFinalLeptons  = 0;
  int nFinalBosons   = 0;
  int nFinalPhotons  = 0;
  int nFinal         = 0;
  for( int i=0; i < event.size(); ++i) {
    if (event[i].isFinal()){
      if (  event[i].id() != 21
        && event[i].id() != 22
        && event[i].id() != 23
        && event[i].idAbs() != 24
        && event[i].id() != 25
        && event[i].colType() == 0)
        nFinalLeptons++;
      if (  event[i].id() == 23
        || event[i].idAbs() == 24
        || event[i].id() == 25)
        nFinalBosons++;
      if (  event[i].id() == 22)
        nFinalPhotons++;
      if ( event[i].isQuark())
        nFinalQuarks++;
      if ( event[i].isGluon())
        nFinalGluons++;
      if ( !event[i].isDiquark() )
        nFinal++;
    }
  }

  // Return highest value if the event does not contain any final state
  // coloured particles.
  if (nFinalQuarks + nFinalGluons == 0) return false;

  // Use MergingHooks functions to get information on the hard process.
  int nLeptons      = nHardOutLeptons();

  // The state is already in the PS region if the number of leptons had been
  // increased bt QED splittings.
  if (nFinalLeptons > nLeptons) return false;

  // If the mumber of photons if larger than in the hard process, QED
  // radiation has pushed the state into the PS region.
  int nPhotons = 0;
  for(int i =0; i< int(hardProcess.hardOutgoing1.size()); ++i)
    if (hardProcess.hardOutgoing1[i] == 22)
      nPhotons++;
  for(int i =0; i< int(hardProcess.hardOutgoing2.size()); ++i)
    if (hardProcess.hardOutgoing2[i] == 22)
      nPhotons++;
  if (nFinalPhotons > nPhotons) return false;

  // Done
  return true;
}

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

// Function to return the minimal kT in the event. If doKTMerging = true, this
// function will be used as a merging scale definition.

double MergingHooks::kTms(const Event& event) {

  // Only check first emission.
  if (!isFirstEmission(event)) return 0.;

  // Find all electroweak decayed bosons in the state.
  vector<int> ewResonancePos;
  ewResonancePos.clear();
  for (int i=0; i < event.size(); ++i)
    if ( abs(event[i].status()) == 22 )
      ewResonancePos.push_back(i);

  // Declare final parton vectors
  vector <int> FinalPartPos;
  FinalPartPos.clear();
  // Search inEvent record for final state partons.
  // Exclude decay products of ew resonance.
  for (int i=0; i < event.size(); ++i){
    if ( event[i].isFinal()
      && event[i].colType() != 0
      && event[i].idAbs()   != 6 ){
      bool isDecayProduct = false;
      for(int j=0; j < int(ewResonancePos.size()); ++j)
        if ( event.isAncestor(i, ewResonancePos[j])
          && event[i].mother2() == 0 )
          isDecayProduct = true;
      // Except for e+e- -> jets, do not check radiation in resonance decays.
      if ( !isDecayProduct
        || getProcessString().compare("e+e->jj") == 0 
        || getProcessString().compare("e+e->(z>jj)") == 0 )
        FinalPartPos.push_back(i);
    }
  }

  // Declare kT algorithm parameters
  double Dparam = 0.4;

  // Find minimal Durham kT in event, using own function: Check
  // definition of separation
  int type = (event[3].colType() == 0
           && event[4].colType() == 0) ? -1 : ktTypeSave;
  // Find minimal kT
  double ktmin = event[0].e();
  for(int i=0; i < int(FinalPartPos.size()); ++i){
    double kt12  = ktmin;
    // Compute separation to the beam axis for hadronic collisions
    if (type == 1 || type == 2) {
      double temp = event[FinalPartPos[i]].pT();
      kt12 = min(kt12, temp);
    }
    // Compute separation to other final state jets
    for(int j=i+1; j < int(FinalPartPos.size()); ++j) {
      double temp = kTdurham( event[FinalPartPos[i]], event[FinalPartPos[j]],
                      type, Dparam);
      kt12 = min(kt12, temp);
    }
    // Keep the minimal Durham separation
    ktmin = min(ktmin,kt12);
  }

  // Return minimal Durham kT
  return ktmin;

}

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

// Function to compute durham y separation from Particle input.

double MergingHooks::kTdurham(const Particle& RadAfterBranch,
  const Particle& EmtAfterBranch, int Type, double D ){

  // Declare return variable
  double ktdur;
  // Save 4-momenta of final state particles
  Vec4 jet1 = RadAfterBranch.p();
  Vec4 jet2 = EmtAfterBranch.p();

  if ( Type == -1) {
    // Get angle between jets for e+e- collisions, make sure that
    // -1 <= cos(theta) <= 1
    double costh;
    if (jet1.pAbs()*jet2.pAbs() <=0.) costh = 1.;
    else {
      costh = costheta(jet1,jet2);
    }
    // Calculate kt durham separation between jets for e+e- collisions
    ktdur = 2.0*min( pow(jet1.e(),2) , (pow(jet2.e(),2)) )*(1.0 - costh);
  } else if ( Type == 1 ){
    // Get delta_y for hadronic collisions:
    // Get mT of first jet
    double mT1sq = jet1.m2Calc() + jet1.pT2();
    double mT1 = 0.;
    if (mT1sq < 0) mT1 = - sqrt(-mT1sq);
    else mT1 = sqrt(mT1sq);
    // Get mT of second jet
    double mT2sq = jet2.m2Calc() + jet2.pT2();
    double mT2 = 0.;
    if (mT2sq < 0) mT2 = - sqrt(-mT2sq);
    else mT2 = sqrt(mT2sq);
    // Get rapidity of first jet
    double y1 = log( ( jet1.e() + abs(jet1.pz()) ) / mT1 );
    if (jet1.pz() < 0) y1 *= -1.;
    // Get rapidity of second jet
    double y2 = log( ( jet2.e() + abs(jet2.pz()) ) / mT2 );
    if (jet2.pz() < 0) y2 *= -1.;
    // Get delta_phi for hadronic collisions
    double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
    double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );  
    double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
    double dPhi = acos( cosdPhi );
    // Calculate kT durham like fastjet,
    // but with rapidity instead of pseudo-rapidity
    ktdur = min( pow(pt1,2),pow(pt2,2) )
          * ( pow(y1-y2,2) + pow(dPhi,2) ) / pow(D,2);
  } else if ( Type == 2 ){
    // Get delta_eta for hadronic collisions
    double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
    double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
    // Get delta_phi and cos(Delta_phi) for hadronic collisions
    double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
    double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );  
    double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
    double dPhi = acos( cosdPhi );
    // Calculate kT durham like fastjet
    ktdur = min( pow(pt1,2),pow(pt2,2) )
          * ( pow(eta1-eta2,2) + pow(dPhi,2) ) / pow(D,2);
  } else if ( Type == 3 ){
    // Get cosh(Delta_eta) for hadronic collisions
    double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
    double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
    double coshdEta = cosh( eta1 - eta2 );
    // Get delta_phi and cos(Delta_phi) for hadronic collisions
    double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
    double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );  
    double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
    // Calculate kT durham separation "SHERPA-like"
    ktdur = 2.0*min( pow(pt1,2),pow(pt2,2) )
          * ( coshdEta - cosdPhi ) / pow(D,2);
  } else {
    ktdur = 0.0;
  }
  // Return kT
  return sqrt(ktdur);
}

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

// Find the minimal Lund pT between coloured partons in the input
// event. If doPTLundMerging = true, this function will be used as a merging
// scale definition.

double MergingHooks::rhoms( const Event& event, bool withColour){

  // Only check first emission.
  if (!isFirstEmission(event)) return 0.;

  // Find all electroweak decayed bosons in the state.
  vector<int> ewResonancePos;
  ewResonancePos.clear();
  for (int i=0; i < event.size(); ++i)
    if ( abs(event[i].status()) == 22 )
      ewResonancePos.push_back(i);

  // Declare final parton vectors
  vector <int> FinalPartPos;
  FinalPartPos.clear();
  // Search inEvent record for final state partons.
  // Exclude decay products of ew resonance.
  for (int i=0; i < event.size(); ++i){
    if ( event[i].isFinal()
      && event[i].colType() != 0
      && event[i].idAbs()   != 6 ){
      bool isDecayProduct = false;
      for(int j=0; j < int(ewResonancePos.size()); ++j)
        if ( event.isAncestor(i, ewResonancePos[j])
          && event[i].mother2() == 0 )
          isDecayProduct = true;
      // Except for e+e- -> jets, do not check radiation in resonance decays.
      if ( !isDecayProduct
        || getProcessString().compare("e+e->jj") == 0 
        || getProcessString().compare("e+e->(z>jj)") == 0 )
        FinalPartPos.push_back(i);
    }
  }

  // Get index of first incoming
  int in1 = 0;
  for (int i=0; i < event.size(); ++i)
    if (abs(event[i].status()) == 41 ){
      in1 = i;
      break;
    }

  // Get index of second incoming
  int in2 = 0;
  for (int i=0; i < event.size(); ++i)
    if (abs(event[i].status()) == 42 ){
      in2 = i;
      break;
    }

  // If no incoming of the cascade are found, try incoming
  if (in1 == 0 || in2 == 0){
    // Find current incoming partons
    for(int i=0; i < int(event.size()); ++i){
      if (event[i].mother1() == 1) in1 = i;
      if (event[i].mother1() == 2) in2 = i;
    }
  }

  // Find minimal pythia pt in event
  double ptmin = event[0].e();
  for(int i=0; i < int(FinalPartPos.size()); ++i){
    double pt12  = ptmin;
    // Compute pythia ISR separation i-jet and first incoming
    if (event[in1].colType() != 0) {
      double temp = rhoPythia( event[in1],
                      event[FinalPartPos[i]], event[in2], -1 );
      pt12 = min(pt12, temp);
    }
    // Compute pythia ISR separation i-jet and second incoming
    if ( event[in2].colType() != 0) {
      double temp = rhoPythia( event[in2],
                      event[FinalPartPos[i]], event[in1], -1 );
      pt12 = min(pt12, temp);
    }

    if (withColour) {
      // Compute pythia FSR separation between two jets,
      // with knowledge of colour connections
      for(int j=0; j < int(FinalPartPos.size()); ++j) {

        // Find recoiler in event
        if ( i!=j ){
          bool isHard = false;
          int radAcl = event[FinalPartPos[i]].acol();
          int radCol = event[FinalPartPos[i]].col();
          int emtAcl = event[FinalPartPos[j]].acol();
          int emtCol = event[FinalPartPos[j]].col();
          int iRec = -1;
          // Check in final state
          if (iRec <= 0 && radAcl > 0 && radAcl != emtCol)
            iRec = findColour(radAcl, FinalPartPos[i], FinalPartPos[j],
                     event,1, isHard);
          if (iRec <= 0 && radCol > 0 && radCol != emtAcl)
            iRec = findColour(radCol, FinalPartPos[i], FinalPartPos[j],
                     event,1, isHard);
          if (iRec <= 0 && emtAcl > 0 && emtAcl != radCol)
            iRec = findColour(emtAcl, FinalPartPos[i], FinalPartPos[j],
                     event,1, isHard);
          if (iRec <= 0 && emtCol > 0 && emtCol != radAcl)
            iRec = findColour(emtCol, FinalPartPos[i], FinalPartPos[j],
                     event,1, isHard);

          // Check in initial state
          if (iRec <= 0 && radAcl > 0 && radAcl != emtCol)
            iRec = findColour(radAcl, FinalPartPos[i], FinalPartPos[j],
                     event,2, isHard);
          if (iRec <= 0 && radCol > 0 && radCol != emtAcl)
            iRec = findColour(radCol, FinalPartPos[i], FinalPartPos[j],
                     event,2, isHard);
          if (iRec <= 0 && emtAcl > 0 && emtAcl != radCol)
            iRec = findColour(emtAcl, FinalPartPos[i], FinalPartPos[j],
                     event,2, isHard);
          if (iRec <= 0 && emtCol > 0 && emtCol != radAcl)
            iRec = findColour(emtCol, FinalPartPos[i], FinalPartPos[j],
                     event,2, isHard);

          if (iRec > 0) {
            double temp = rhoPythia( event[FinalPartPos[i]],
                                     event[FinalPartPos[j]],
                                     event[iRec], 1 );
            pt12 = min(pt12, temp);
          }
        }

        // If minimal pT below shower cut-off, return
        if (pt12 < 0.4) return pt12;

      }

    } else {
      // Compute pythia FSR separation between two jets,
      // without any knowledge of colour connections
      for(int j=0; j < int(FinalPartPos.size()); ++j) {
        for(int k=0; k < int(FinalPartPos.size()); ++k) {
          // Allow any parton as recoiler
          if ( (i != j) && (i != k) && (j != k) ){

            double temp = 0.;
            // Only check splittings allowed by flavour, e.g.
            // only q -> qg and g -> qqbar
            temp = rhoPythia( event[FinalPartPos[i]],
                              event[FinalPartPos[j]],
                              event[FinalPartPos[k]], 1 );
            pt12 = min(pt12, temp);
          }
        }
        // If minimal pT below shower cut-off, return
        if (pt12 < 0.4) return pt12;
      }

      // Compute pythia FSR separation between two jets, with initial recoiler
      // without any knowledge of colour connections
      if ( event[in1].colType() != 0 && event[in2].colType() != 0) {
        for(int j=0; j < int(FinalPartPos.size()); ++j) {
          // Allow both initial partons as recoiler
          if ( i != j ){
            // Check with first initial as recoiler
            double temp = rhoPythia( event[FinalPartPos[i]],
                                     event[FinalPartPos[j]],
                                     event[in1], 1 );
            pt12 = min(pt12, temp);
            // Check with second initial as recoiler
            temp        = rhoPythia( event[FinalPartPos[i]],
                                     event[FinalPartPos[j]],
                                     event[in2], 1 );
            pt12 = min(pt12, temp);
          }

          // If minimal pT below shower cut-off, return
          if (pt12 < 0.4) return pt12;
        }
      }

    }
    // Reset minimal y separation
    ptmin = min(ptmin,pt12);
  }

  return ptmin;

}

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

// Function to compute "pythia pT separation" from Particle input, as a helper
// for rhoms(...).

double MergingHooks::rhoPythia(const Particle& RadAfterBranch,
              const Particle& EmtAfterBranch,
              const Particle& RecAfterBranch, int ShowerType){

  // Save type: 1 = FSR pT definition, else ISR definition
  int Type   = ShowerType;
  // Calculate virtuality of splitting
  int sign = (Type==1) ? 1 : -1;
  Vec4 Q(RadAfterBranch.p() + sign*EmtAfterBranch.p());
  double Qsq = sign * Q.m2Calc();
  // Mass term of radiator
  double m2Rad = ( includeMassive()
               && abs(RadAfterBranch.id()) >= 4
               && abs(RadAfterBranch.id()) < 7)
               ? pow(particleDataPtr->m0(RadAfterBranch.id()), 2)
               : 0.;
  // Construct 2->3 variables for FSR
  Vec4   sum     = RadAfterBranch.p() + RecAfterBranch.p()
                 + EmtAfterBranch.p();
  double m2Dip = sum.m2Calc();
  double x1 = 2. * (sum * RadAfterBranch.p()) / m2Dip;
  double x3 = 2. * (sum * EmtAfterBranch.p()) / m2Dip;
  // Construct momenta of dipole before/after splitting for ISR 
  Vec4 qBR(RadAfterBranch.p() - EmtAfterBranch.p() + RecAfterBranch.p());
  Vec4 qAR(RadAfterBranch.p() + RecAfterBranch.p());
  // Calculate z of splitting, different for FSR and ISR
  double z = (Type==1) ? x1/(x1+x3)
                     : (qBR.m2Calc())/( qAR.m2Calc());
  // Separation of splitting, different for FSR and ISR
  double pTpyth = (Type==1) ? z*(1.-z) : (1.-z);
  // pT^2 = separation*virtuality
  pTpyth *= (Qsq - sign*m2Rad);
  if (pTpyth < 0.) pTpyth = 0.;
  // Return pT
  return sqrt(pTpyth);
}

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

// Function to find a colour (anticolour) index in the input event.
// Helper for rhoms
// IN  int col       : Colour tag to be investigated
//     int iExclude1 : Identifier of first particle to be excluded
//                     from search
//     int iExclude2 : Identifier of second particle to be excluded
//                     from  search
//     Event event   : event to be searched for colour tag
//     int type      : Tag to define if col should be counted as
//                      colour (type = 1) [->find anti-colour index
//                                         contracted with col]
//                      anticolour (type = 2) [->find colour index
//                                         contracted with col]
// OUT int           : Position of particle in event record
//                     contraced with col [0 if col is free tag]

int MergingHooks::findColour(int col, int iExclude1, int iExclude2,
      const Event& event, int type, bool isHardIn){

  bool isHard = isHardIn;
  int index = 0;

  if (isHard){
    // Search event record for matching colour & anticolour
    for(int n = 0; n < event.size(); ++n) {
      if ( n != iExclude1 && n != iExclude2
        && event[n].colType() != 0
        &&(   event[n].status() > 0          // Check outgoing
           || event[n].status() == -21) ) {  // Check incoming
         if ( event[n].acol() == col ) {
          index = -n;
          break;
        }
        if ( event[n].col()  == col ){
          index =  n;
          break;
        }
      }
    }
  } else {

    // Search event record for matching colour & anticolour
    for(int n = 0; n < event.size(); ++n) {
      if (  n != iExclude1 && n != iExclude2
        && event[n].colType() != 0
        &&(   event[n].status() == 43        // Check outgoing from ISR
           || event[n].status() == 51        // Check outgoing from FSR
           || event[n].status() == 52        // Check outgoing from FSR
           || event[n].status() == -41       // first initial
           || event[n].status() == -42) ) {  // second initial
        if ( event[n].acol() == col ) {
          index = -n;
          break;
        }
        if ( event[n].col()  == col ){
          index =  n;
          break;
        }
      }
    }
  }
  // if no matching colour / anticolour has been found, return false
  if ( type == 1 && index < 0) return abs(index);
  if ( type == 2 && index > 0) return abs(index);

  return 0;
}

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

// Function to check if the properties of the input particle should be
// checked against the cut-based merging scale defintion.

bool MergingHooks::checkAgainstCut( const Particle& particle){

  // Do not check uncoloured particles.
  if (particle.colType() == 0)
    return false;
  // Do not check tops and bottoms.
  if (particle.idAbs() == 5 || particle.idAbs() == 6)
    return false;
  // Done
  return true;

}

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

// Find the if the event passes the Delta R_{ij}, pT_{i} and Q_{ij} cuts on
// the matrix element, as needed for cut-based merging scale definition.

double MergingHooks::cutbasedms( const Event& event ){

  // Only check first emission.
  if (!isFirstEmission(event)) return -1.;

  // Save allowed final state particles
  vector<int> partons;
  for( int i=0; i < event.size(); ++i)
    if ( event[i].isFinal() && checkAgainstCut(event[i]) )
      partons.push_back(i);

  // Declare overall veto
  bool doVeto = false;
  // Declare vetoes
  bool vetoPT  = false;
  bool vetoRjj = false;
  bool vetoMjj = false;
  // Declare cuts used in matrix element
  double pTjmin = pTiMS();
  double mjjmin = QijMS();
  double rjjmin = dRijMS();

  // Declare minimum values
  double minPT  = event[0].e();
  double minRJJ = 10.;
  double minMJJ = event[0].e();

  // Check matrix element cuts
  for( int i=0; i < int(partons.size()); ++i){
    // Save pT value
    minPT = min(minPT,event[partons[i]].pT());

    // Check two-parton cuts
    for( int j=0; j < int(partons.size()); ++j){
      if (i!=j){

        // Save delta R value
        minRJJ = min(minRJJ, deltaRij( event[partons[i]].p(),
              event[partons[j]].p()) );
        // Save delta R value
        minMJJ = min(minMJJ, ( event[partons[i]].p()
                              +event[partons[j]].p() ).mCalc() );

      }
    }
  // Done with cut evaluation
  }

  // Check if all partons are in the matrix element region
  if (minPT  > pTjmin) vetoPT  = true;
  if (minRJJ > rjjmin) vetoRjj = true;
  if (minMJJ > mjjmin) vetoMjj = true;

  // In the matrix element calculation, we have rejected the events if any of 
  // the cuts had not been fulfilled,
  // i.e. we are in the matrix element domain if all cuts are fulfilled.
  // We veto any emission in the ME region.
  // Disregard the two-parton cuts if only one parton in the final state
  if (int(partons.size() == 1))
    doVeto = vetoPT;
  else
  // Veto if the combination of cuts would be in the ME region
    doVeto = vetoPT && vetoRjj && vetoMjj;

  // If event is above merging scale, veto
  if (doVeto) return 1.;

  // Else, do nothing
  return -1.;

}

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

// Function to compute Delta R separation from 4-vector input.

double MergingHooks::deltaRij(Vec4 jet1, Vec4 jet2){

  // Declare return variable
  double deltaR = 0.;
  // Get delta_eta and cosh(Delta_eta) for hadronic collisions
  double eta1 = 0.5*log( (jet1.e() + jet1.pz()) / (jet1.e() - jet1.pz()) );
  double eta2 = 0.5*log( (jet2.e() + jet2.pz()) / (jet2.e() - jet2.pz()) );
  // Get delta_phi and cos(Delta_phi) for hadronic collisions
  double pt1 = sqrt( pow(jet1.px(),2) + pow(jet1.py(),2) );
  double pt2 = sqrt( pow(jet2.px(),2) + pow(jet2.py(),2) );  
  double cosdPhi = ( jet1.px()*jet2.px() + jet1.py()*jet2.py() ) / (pt1*pt2);
  double dPhi = acos( cosdPhi );
  // Calculate kT durham like fastjet
  deltaR = sqrt(pow(eta1-eta2,2) + pow(dPhi,2));
  // Return kT
  return deltaR;
}

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

} // end namespace Pythia8