ProcessContainer.cc

// ProcessContainer.cc is a part of the PYTHIA event generator.
// Copyright (C) 2018 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
// ProcessContainer and SetupContainers classes.

#include "Pythia8/ProcessContainer.h"

// Internal headers for special processes.
#include "Pythia8/SigmaCompositeness.h"
#include "Pythia8/SigmaEW.h"
#include "Pythia8/SigmaExtraDim.h"
#include "Pythia8/SigmaGeneric.h"
#include "Pythia8/SigmaHiggs.h"
#include "Pythia8/SigmaLeftRightSym.h"
#include "Pythia8/SigmaLeptoquark.h"
#include "Pythia8/SigmaNewGaugeBosons.h"
#include "Pythia8/SigmaQCD.h"
#include "Pythia8/SigmaSUSY.h"
#include "Pythia8/SigmaDM.h"

namespace Pythia8 {

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

// ProcessContainer class.
// Information allowing the generation of a specific process.

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

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

// Number of event tries to check maximization finding reliability.
const int ProcessContainer::N12SAMPLE = 100;

// Ditto, but increased for 2 -> 3 processes.
const int ProcessContainer::N3SAMPLE  = 1000;

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

// Initialize phase space and counters.
// Argument isFirst distinguishes two hard processes in same event.

bool ProcessContainer::init(bool isFirst, Info* infoPtrIn,
  Settings& settings, ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
  BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, Couplings* couplingsPtr,
  SigmaTotal* sigmaTotPtrIn, ResonanceDecays* resDecaysPtrIn,
  SLHAinterface* slhaInterfacePtr, UserHooks* userHooksPtrIn,
  GammaKinematics* gammaKinPtrIn) {

  // Extract info about current process from SigmaProcess object.
  isLHA       = sigmaProcessPtr->isLHA();
  isNonDiff   = sigmaProcessPtr->isNonDiff();
  isResolved  = sigmaProcessPtr->isResolved();
  isDiffA     = sigmaProcessPtr->isDiffA();
  isDiffB     = sigmaProcessPtr->isDiffB();
  isDiffC     = sigmaProcessPtr->isDiffC();
  isCentralDiff = isDiffC;
  isDoubleDiff  = isDiffA && isDiffB;
  isSingleDiff  = isDiff && !isDoubleDiff  && !isCentralDiff;
  isQCD3body  = sigmaProcessPtr->isQCD3body();
  int nFin    = sigmaProcessPtr->nFinal();
  lhaStrat    = (isLHA) ? lhaUpPtr->strategy() : 0;
  lhaStratAbs = abs(lhaStrat);
  allowNegSig = sigmaProcessPtr->allowNegativeSigma();

  // Switch to ensure that the scales set in a LH input event are
  // respected, even in resonance showers.
  useStrictLHEFscales = settings.flag("Beams:strictLHEFscale");

  // Maximal number of events for more in-depth control of
  // how many events are selected/tried/accepted.
  nTryRequested = settings.mode("Main:numberOfTriedEvents");
  nSelRequested = settings.mode("Main:numberOfSelectedEvents");
  nAccRequested = settings.mode("Main:numberOfAcceptedEvents");

  // Flag for maximum violation handling.
  increaseMaximum = settings.flag("PhaseSpace:increaseMaximum");

  // Store whether beam particle has a photon and save the mode.
  beamAPtr           = beamAPtrIn;
  beamBPtr           = beamBPtrIn;
  gammaKinPtr        = gammaKinPtrIn;
  beamHasGamma       = settings.flag("PDF:lepton2gamma");
  int gammaMode      = settings.mode("Photon:ProcessType");
  bool resolvedGamma = beamHasGamma && (gammaMode < 2);

  // Use external photon flux.
  externalFlux = (settings.mode("PDF:lepton2gammaSet") == 2);

  // Pick and create phase space generator. Send pointers where required.
  if (phaseSpacePtr != 0) ;
  else if (isLHA)       phaseSpacePtr = new PhaseSpaceLHA();
  else if (isNonDiff && !beamHasGamma)
                        phaseSpacePtr = new PhaseSpace2to2nondiffractive();
  else if (isNonDiff && resolvedGamma)
    phaseSpacePtr = new PhaseSpace2to2nondiffractiveGamma();
  else if (!isResolved && !isDiffA  && !isDiffB  && !isDiffC )
                        phaseSpacePtr = new PhaseSpace2to2elastic();
  else if (!isResolved && !isDiffA  && !isDiffB && isDiffC)
                        phaseSpacePtr = new PhaseSpace2to3diffractive();
  else if (!isResolved) phaseSpacePtr = new PhaseSpace2to2diffractive(
                                        isDiffA, isDiffB);
  else if (nFin == 1)   phaseSpacePtr = new PhaseSpace2to1tauy();
  else if (nFin == 2)   phaseSpacePtr = new PhaseSpace2to2tauyz();
  else if (isQCD3body)  phaseSpacePtr = new PhaseSpace2to3yyycyl();
  else                  phaseSpacePtr = new PhaseSpace2to3tauycyl();

  // Store pointers and perform simple initialization.
  infoPtr         = infoPtrIn;
  particleDataPtr = particleDataPtrIn;
  rndmPtr         = rndmPtrIn;
  resDecaysPtr    = resDecaysPtrIn;
  sigmaTotPtr     = sigmaTotPtrIn;
  userHooksPtr    = userHooksPtrIn;
  canVetoResDecay = (userHooksPtr != 0)
                  ? userHooksPtr->canVetoResonanceDecays() : false;
  if (isLHA) {
    sigmaProcessPtr->setLHAPtr(lhaUpPtr);
    phaseSpacePtr->setLHAPtr(lhaUpPtr);
  }
  sigmaProcessPtr->init(infoPtr, &settings, particleDataPtr, rndmPtr,
    beamAPtr, beamBPtr, couplingsPtr, sigmaTotPtr, slhaInterfacePtr);

  // Store the state of photon beams using inFlux: 0 = not a photon beam;
  // 1 = resolved photon; 2 = unresolved photon.
  string inState = sigmaProcessPtr->inFlux();
  beamAgammaMode = 0;
  beamBgammaMode = 0;
  gammaModeEvent = 0;
  if ( beamAPtr->isGamma() || (beamHasGamma && beamAPtr->isLepton() ) ) {
    if ( inState == "gmg" || inState == "gmq" || inState == "gmgm" )
      beamAgammaMode = 2;
    else if ( !(beamAPtr->isHadron()) ) beamAgammaMode = 1;
  }
  if ( beamBPtr->isGamma() || (beamHasGamma && beamBPtr->isLepton() ) ) {
    if ( inState == "ggm" || inState == "qgm" || inState == "gmgm" )
      beamBgammaMode = 2;
    else if ( !(beamBPtr->isHadron()) ) beamBgammaMode = 1;
  }

  // Save the photon modes and propagate further.
  if ( ( beamAPtr->isGamma() || beamBPtr->isGamma() ) || beamHasGamma )
    setBeamModes();

  // Save the state of photon beams.
  beamAhasResGamma = beamAPtr->hasResGamma();
  beamBhasResGamma = beamBPtr->hasResGamma();
  beamHasResGamma  = beamAhasResGamma || beamBhasResGamma;

  // Initialize also phaseSpace pointer.
  phaseSpacePtr->init( isFirst, sigmaProcessPtr, infoPtr, &settings,
    particleDataPtr, rndmPtr, beamAPtr,  beamBPtr, couplingsPtr, sigmaTotPtr,
    userHooksPtr);

  // Send the pointer to gammaKinematics for non-diffractive processes.
  if ( (beamAhasResGamma || beamBhasResGamma) && isNonDiffractive() )
    phaseSpacePtr->setGammaKinPtr( gammaKinPtr);

  // Reset cross section statistics.
  nTry      = 0;
  nSel      = 0;
  nAcc      = 0;
  nTryStat  = 0;
  sigmaMx   = 0.;
  sigmaSum  = 0.;
  sigma2Sum = 0.;
  sigmaNeg  = 0.;
  sigmaAvg  = 0.;
  sigmaFin  = 0.;
  deltaFin  = 0.;
  wtAccSum  = 0.;

  // Reset temporary cross section summand.
  sigmaTemp   = 0.;
  sigma2Temp  = 0.;

  // Initialize process and allowed incoming partons.
  sigmaProcessPtr->initProc();
  if (!sigmaProcessPtr->initFlux()) return false;

  // Find maximum of differential cross section * phasespace.
  bool physical       = phaseSpacePtr->setupSampling();
  sigmaMx             = phaseSpacePtr->sigmaMax();
  double sigmaHalfWay = sigmaMx;

  // Separate signed maximum needed for LHA with negative weight.
  sigmaSgn            = phaseSpacePtr->sigmaSumSigned();

  // Check maximum by a few events, and extrapolate a further increase.
  if (physical & !isLHA) {
    int nSample = (nFin < 3) ? N12SAMPLE : N3SAMPLE;
    for (int iSample = 0; iSample < nSample; ++iSample) {
      bool test = false;
      while (!test) test = phaseSpacePtr->trialKin(false);
      if (iSample == nSample/2) sigmaHalfWay = phaseSpacePtr->sigmaMax();
    }
    double sigmaFullWay = phaseSpacePtr->sigmaMax();
    sigmaMx = (sigmaHalfWay > 0.) ? pow2(sigmaFullWay) / sigmaHalfWay
                                  : sigmaFullWay;
    phaseSpacePtr->setSigmaMax(sigmaMx);
  }

  // Allow Pythia to overwrite incoming beams or parts of Les Houches input.
  idRenameBeams = settings.mode("LesHouches:idRenameBeams");
  setLifetime   = settings.mode("LesHouches:setLifetime");
  setQuarkMass  = settings.mode("LesHouches:setQuarkMass");
  setLeptonMass = settings.mode("LesHouches:setLeptonMass");
  mRecalculate  = settings.parm("LesHouches:mRecalculate");
  matchInOut    = settings.flag("LesHouches:matchInOut");
  for (int i = 0; i < 6; ++i) idNewM[i] = i;
  for (int i = 6; i < 9; ++i) idNewM[i] = 2 * i - 1;
  for (int i = 1; i < 9; ++i) mNewM[i]  = particleDataPtr->m0(idNewM[i]);

  // Done.
  return physical;
}

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

// Generate a trial event; selected or not.

bool ProcessContainer::trialProcess() {

  // Weights for photon flux oversampling with external flux.
  double wtPDF  = 1.;
  double wtFlux = 1.;

  // For photon beams set PDF pointer to resolved or unresolved.
  if ( beamAPtr->isGamma() || beamBPtr->isGamma() || beamHasGamma )
    setBeamModes();

  // Loop over tries only occurs for Les Houches strategy = +-2.
  for (int iTry = 0;  ; ++iTry) {

    // Generate a trial phase space point, if meaningful.
    if (sigmaMx == 0.) return false;
    infoPtr->setEndOfFile(false);
    bool repeatSame = (iTry > 0);
    bool physical = phaseSpacePtr->trialKin(true, repeatSame);

    // For acceptable kinematics sample also kT for photons from leptons.
    if ( physical && !isNonDiffractive() && beamHasGamma ) {

      // Save the x_gamma values for unresolved photons.
      if ( !beamAhasResGamma ) beamAPtr->xGamma( x1());
      if ( !beamBhasResGamma ) beamBPtr->xGamma( x2());

      // Sample the kinematics of virtual photons.
      if ( !gammaKinPtr->sampleKTgamma() ) physical = false;

      // Processes with direct photons rescale momenta and cross section.
      if ( physical && !(beamAhasResGamma && beamBhasResGamma) ) {
        double sHatNew = gammaKinPtr->calcNewSHat( phaseSpacePtr->sHat() );
        phaseSpacePtr->rescaleSigma( sHatNew);
        phaseSpacePtr->rescaleMomenta( sHatNew);
      }

      // With external photon flux calculate weight wrt. approximated fluxes.
      if ( beamHasGamma && externalFlux ) {
        wtPDF  = phaseSpacePtr->weightGammaPDFApprox();
        wtFlux = gammaKinPtr->fluxWeight();
      }

    }

    // Flag to check if more events should be generated.
    bool doTryNext = true;

    // Note if at end of Les Houches file, else do statistics.
    if (isLHA && !physical) infoPtr->setEndOfFile(true);
    else {

      // Check if the number of selected events should be limited.
      if (nSelRequested > 0) doTryNext = (nTry < nSelRequested);

      // Only add to counter if event should count towards cross section.
      if (doTryNext) ++nTry;

      // Statistics for separate Les Houches process codes. Insert new codes.
      if (isLHA) {
        int codeLHANow = lhaUpPtr->idProcess();
        int iFill = -1;
        for (int i = 0; i < int(codeLHA.size()); ++i)
          if (codeLHANow == codeLHA[i]) iFill = i;
        if (iFill >= 0) {
          // Only add to counter if event should count towards cross section.
          if (doTryNext) ++nTryLHA[iFill];
        } else {
          codeLHA.push_back(codeLHANow);
          nTryLHA.push_back(1);
          nSelLHA.push_back(0);
          nAccLHA.push_back(0);
          for (int i = int(codeLHA.size()) - 1; i > 0; --i) {
            if (codeLHA[i] < codeLHA[i - 1]) {
              swap(codeLHA[i], codeLHA[i - 1]);
              swap(nTryLHA[i], nTryLHA[i - 1]);
              swap(nSelLHA[i], nSelLHA[i - 1]);
              swap(nAccLHA[i], nAccLHA[i - 1]);
            }
            else break;
          }
        }
      }
    }

    // Incoming top beams will not be handled, although allowed by LHAPDF.
    if (isLHA && (abs(lhaUpPtr->id(1)) == 6 || abs(lhaUpPtr->id(2)) == 6)) {
      infoPtr->errorMsg("Error in ProcessContainer::trialProcess(): top not "
        "allowed incoming beam parton; event skipped");
      return false;
    }

    // Possibly fail, else cross section.
    if (!physical) return false;
    double sigmaNow = phaseSpacePtr->sigmaNow();

    // For photons with external flux correct the cross section.
    if (beamHasGamma && externalFlux) sigmaNow *= wtFlux * wtPDF;

    // Tell if this event comes with weight from cross section.
    double sigmaWeight = 1.;
    if (!isLHA && !increaseMaximum && sigmaNow > sigmaMx)
      sigmaWeight = sigmaNow / sigmaMx;
    if ( lhaStrat < 0 && sigmaNow < 0.) sigmaWeight = -1.;
    if ( lhaStratAbs == 4) sigmaWeight = sigmaNow;

    // Also compensating weight from biased phase-space selection.
    double biasWeight = phaseSpacePtr->biasSelectionWeight();
    weightNow = sigmaWeight * biasWeight;

    // Set event weight to zero if this event should not count.
    if (!doTryNext) weightNow = 0.0;

    infoPtr->setWeight( weightNow, lhaStrat);

    // Check that not negative cross section when not allowed.
    if (!allowNegSig) {
      if (sigmaNow < sigmaNeg) {
        infoPtr->errorMsg("Warning in ProcessContainer::trialProcess: neg"
          "ative cross section set 0", "for " +  sigmaProcessPtr->name() );
        sigmaNeg = sigmaNow;
      }
      if (sigmaNow < 0.) sigmaNow = 0.;
    }

    // Cross section of process may come with a weight. Update sum.
    double sigmaAdd = sigmaNow * biasWeight;
    if (lhaStratAbs == 2 || lhaStratAbs == 3) sigmaAdd = sigmaSgn;

    // Do not add event weight to the cross section if it should not count.
    if (!doTryNext) {
      sigmaAdd=0.;
      // Reset stored weight.
      sigmaTemp  = 0.;
      sigma2Temp = 0.;
    }

    // Only update once an event is accepted, as is needed if the event weight
    // can still change by a finite amount due to reweighting.
    if (lhaStratAbs >= 3 ) {
      sigmaTemp  = sigmaAdd;
      sigma2Temp = pow2(sigmaAdd);
    } else {
      sigmaTemp  += sigmaAdd;
      sigma2Temp += pow2(sigmaAdd);
    }

    // Check if maximum violated.
    newSigmaMx = phaseSpacePtr->newSigmaMax();
    if (newSigmaMx) sigmaMx = phaseSpacePtr->sigmaMax();

    // Select or reject trial point. Statistics.
    bool select = true;
    if (lhaStratAbs < 3) select
      = (newSigmaMx || rndmPtr->flat() * abs(sigmaMx) < abs(sigmaNow));
    if (select) {
      // Only add to counter if event should count towards cross section.
      if (doTryNext)               ++nSel;

     if (isLHA) {
        int codeLHANow = lhaUpPtr->idProcess();
        int iFill = -1;
        for (int i = 0; i < int(codeLHA.size()); ++i)
          if (codeLHANow == codeLHA[i]) iFill = i;
        // Only add to counter if event should count towards cross section.
        if (iFill >= 0 && doTryNext) ++nSelLHA[iFill];
      }
    }
    if (select || lhaStratAbs != 2) return select;

  }

}

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

// Accumulate statistics after user veto, including LHA code.

void ProcessContainer::accumulate() {

  // Only update for non-vanishing weight - vanishing weight means
  // event should not count as accepted.
  double wgtNow = infoPtr->weight();
  if (wgtNow != 0.0) {

    ++nAcc;

    // infoPtr->weight() performs coversion to pb if lhaStratAbs = 4
    if (lhaStratAbs == 4) wgtNow /= 1e9;

    wtAccSum += wgtNow;
    if (isLHA) {
      int codeLHANow = lhaUpPtr->idProcess();
      int iFill = -1;
      for (int i = 0; i < int(codeLHA.size()); ++i)
        if (codeLHANow == codeLHA[i]) iFill = i;
      if (iFill >= 0) ++nAccLHA[iFill];
    }
  }

}

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

// Pick flavours and colour flow of process.

bool ProcessContainer::constructState() {

  // Construct flavour and colours for selected event.
  if (isResolved && !isNonDiff) sigmaProcessPtr->pickInState();
  sigmaProcessPtr->setIdColAcol();

  // Set beam modes correctly for given process.
  if ( beamAPtr->isResolvedUnresolved() || beamBPtr->isResolvedUnresolved() )
    setBeamModes();

  // Done.
  return true;

}

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

// Set the photon modes according to the process.

void ProcessContainer::setBeamModes() {

  // Set the modes for the current beams.
  beamAPtr->setGammaMode(beamAgammaMode);
  beamBPtr->setGammaMode(beamBgammaMode);

  // Set the gammaMode for given process.
  if      (beamAgammaMode <= 1 && beamBgammaMode <= 1) gammaModeEvent = 1;
  else if (beamAgammaMode <= 1 && beamBgammaMode == 2) gammaModeEvent = 2;
  else if (beamAgammaMode == 2 && beamBgammaMode <= 1) gammaModeEvent = 3;
  else if (beamAgammaMode == 2 && beamBgammaMode == 2) gammaModeEvent = 4;
  else                                                 gammaModeEvent = 0;

  // Propagate gammaMode and VMD state to info pointer.
  infoPtr->setGammaMode(gammaModeEvent);
  if (gammaModeEvent != 4) {
    double scale[3] = {0.00331524, 0.00030905, 0.00039639};
    pair<int, int> idAB = chooseVMDstates(beamAPtr->id(), beamBPtr->id());
    if (idAB.first == 113 || idAB.first == 223 || idAB.first == 333) {
      double mA  = particleDataPtr->mSel(idAB.first);
      double scA = scale[idAB.first/100 - 1];
      beamAPtr->setVMDstate(true, idAB.first, mA, scA, false);
      infoPtr->setVMDstateA(true, idAB.first, mA, scA);
    }
    if (idAB.second == 113 || idAB.second == 223 || idAB.second == 333) {
      double mB  = particleDataPtr->mSel(idAB.second);
      double scB = scale[idAB.second/100 - 1];
      beamBPtr->setVMDstate(true, idAB.second, mB, scB, false);
      infoPtr->setVMDstateB(true, idAB.second, mB, scB);
    }
  }

}

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

// VMD processes: Choose particle.

pair<int, int> ProcessContainer::chooseVMDstates( int idA, int idB) {

  // Constants and initial values.
  double gammaFac[3] = {2.2, 23.6, 18.4};
  double alphaEM     = 0.00729353;
  double idVMD[3]    = {113, 223, 333};
  double pVP[3]      = {};
  double pVV[3][3]   = {};
  double pSum        = 0.;

  // gamma-gamma.
  if (idA == 22 && idB == 22) {
    for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j) {
      // Evaluate the cross sections individually.
      sigmaTotPtr->calc(idVMD[i], idVMD[j], infoPtr->eCM());
      pVV[i][j] = pow2(alphaEM) / (gammaFac[i] * gammaFac[j]);
      if (isSingleDiff && isDiffA)      pVV[i][j] *= sigmaTotPtr->sigmaXB();
      else if (isSingleDiff && isDiffB) pVV[i][j] *= sigmaTotPtr->sigmaAX();
      else if (isDoubleDiff)            pVV[i][j] *= sigmaTotPtr->sigmaXX();
      else return make_pair(idA, idB);
      pSum     += pVV[i][j];
    }
    // Choose VMD states based on relative fractions.
    double pickMode = rndmPtr->flat() * pSum;
    for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j) {
      pickMode -= pVV[i][j];
      if (pickMode < 0.) return make_pair(113 + 110 * i, 113 + 110 * j);
    }
    return make_pair(113, 113);

  // gamma + p.
  } else if (idA == 22 && idB == 2212) {
    for (int i = 0; i < 3; ++i) {
      // Evaluate the cross sections individually.
      sigmaTotPtr->calc(idVMD[i], 2212, infoPtr->eCM());
      pVP[i] = alphaEM / gammaFac[i];
      if (isSingleDiff && isDiffA)      pVP[i] *= sigmaTotPtr->sigmaXB();
      else if (isSingleDiff && isDiffB) pVP[i] *= sigmaTotPtr->sigmaAX();
      else if (isDoubleDiff)            pVP[i] *= sigmaTotPtr->sigmaXX();
      else return make_pair(idA, idB);
      pSum     += pVP[i];
    }
    // Choose VMD state based on relative fractions.
    double pickMode = rndmPtr->flat() * pSum;
    if (pickMode > pVP[1] + pVP[2]) return make_pair(113, 2212);
    else if (pickMode > pVP[2])     return make_pair(223, 2212);
    else                            return make_pair(333, 2212);

  // p + gamma.
  } else if (idA == 2212 && idB == 22) {
    for (int i = 0; i < 3; ++i) {
      // Evaluate the cross sections individually.
      sigmaTotPtr->calc(2212, idVMD[i], infoPtr->eCM());
      pVP[i] = alphaEM / gammaFac[i];
      if (isSingleDiff && isDiffA)      pVP[i] *= sigmaTotPtr->sigmaXB();
      else if (isSingleDiff && isDiffB) pVP[i] *= sigmaTotPtr->sigmaAX();
      else if (isDoubleDiff)            pVP[i] *= sigmaTotPtr->sigmaXX();
      else return make_pair(idA, idB);
      pSum     += pVP[i];
    }
    // Choose VMD state based on relative fractions.
    double pickMode = rndmPtr->flat() * pSum;
    if (pickMode > pVP[1] + pVP[2]) return make_pair(2212, 113);
    else if (pickMode > pVP[2])     return make_pair(2212, 223);
    else                            return make_pair(2212, 333);

  // Do nothing for other beam combinations.
  } else return make_pair(idA, idB);

}

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

// Give the hard subprocess with kinematics.

bool ProcessContainer::constructProcess( Event& process, bool isHardest) {

  // Construct kinematics from selected phase space point.
  if (!phaseSpacePtr->finalKin()) return false;
  int nFin = sigmaProcessPtr->nFinal();

  // Basic info on process.
  if (isHardest) infoPtr->setType( name(), code(), nFin, isNonDiff,
    isResolved, isDiffA, isDiffB, isDiffC, isLHA);

  // Save sampled values for further use, requires info stored by setType.
  if ( beamHasGamma && !isNonDiffractive() ) gammaKinPtr->finalize();

  // Rescale the momenta again when unresolved photons after finalKin.
  if ( beamHasGamma && !(beamAhasResGamma && beamBhasResGamma) ) {
    double sHatNew = infoPtr->sHatNew();
    phaseSpacePtr->rescaleMomenta( sHatNew);
  }

  // Let hard process record begin with the event as a whole and
  // the two incoming beam particles.
  int idA = infoPtr->idA();
  if (abs(idA) == idRenameBeams) idA = 16;
  int idB = infoPtr->idB();
  if (abs(idB) == idRenameBeams) idB = -16;
  process.append( 90, -11, 0, 0, 0, 0, 0, 0,
    Vec4(0., 0., 0., infoPtr->eCM()), infoPtr->eCM(), 0. );
  process.append( idA, -12, 0, 0, 0, 0, 0, 0,
    Vec4(0., 0., infoPtr->pzA(), infoPtr->eA()), infoPtr->mA(), 0. );
  process.append( idB, -12, 0, 0, 0, 0, 0, 0,
    Vec4(0., 0., infoPtr->pzB(), infoPtr->eB()), infoPtr->mB(), 0. );

  // Add intermediate gammas for lepton -> gamma -> parton processes
  // for both non-diffractive and hard processes, including direct-resolved.
  // Add a copy of hadron beam when e->gamma + p.
  int nOffsetGamma = 0;
  bool isGammaHadronDir = (beamAgammaMode == 2 && beamBgammaMode == 0)
                       || (beamAgammaMode == 0 && beamBgammaMode == 2);
  if ( beamHasResGamma || (isGammaHadronDir && beamHasGamma) ) {
    double xGm1 = beamAPtr->xGamma();
    if ( beamAPtr->isHadron()) {
      process.append( beamAPtr->id(), -13, 1, 0, 0, 0, 0, 0,
        Vec4(0., 0., infoPtr->pzA(), infoPtr->eA()), beamAPtr->m(), 0. );
    } else {
      process.append( 22, -13, 1, 0, 0, 0, 0, 0,
        Vec4(0., 0., xGm1*infoPtr->pzA(), xGm1*infoPtr->eA()), 0, 0. );
    }
    process[1].daughter1(3);
    ++nOffsetGamma;
    double xGm2 = beamBPtr->xGamma();
    if ( beamBPtr->isHadron()) {
      process.append( beamBPtr->id(), -13, 2, 0, 0, 0, 0, 0,
        Vec4(0., 0., infoPtr->pzB(), infoPtr->eB()), beamBPtr->m(), 0. );
    } else {
      process.append( 22, -13, 2, 0, 0, 0, 0, 0,
        Vec4(0., 0., xGm2*infoPtr->pzB(), xGm2*infoPtr->eB()), 0, 0. );
    }
    process[1 + nOffsetGamma].daughter1(3 + nOffsetGamma);
    ++nOffsetGamma;
  }

  // For nondiffractive process no interaction selected so far, so done.
  if (isNonDiff) return true;

  // Entries 3 and 4, now to be added, come from 1 and 2.
  // Offset from normal locations possible due to intermediate photons.
  process[1 + nOffsetGamma].daughter1(3 + nOffsetGamma);
  process[2 + nOffsetGamma].daughter1(4 + nOffsetGamma);
  double scale  = 0.;
  double scalup = 0.;

  // For DiffC entries 3 - 5 come jointly from 1 and 2 (to keep HepMC happy).
  if (isDiffC) {
    process[1].daughters(3, 5);
    process[2].daughters(3, 5);
  }

  // Insert the subprocess partons - resolved processes.
  int idRes = sigmaProcessPtr->idSChannel();
  if (isResolved && !isLHA) {

    // NOAM: Mothers and daughters without/with intermediate state.
    int m_M1 = 3;
    int m_M2 = 4;
    int m_D1 = 5;
    int m_D2 = (nFin == 1) ? 0 : 4 + nFin;
    if (idRes != 0) {
      m_M1   = 5;
      m_M2   = 0;
      m_D1   = 5;
      m_D2   = 0;
    }

    // Find scale from which to begin MPI/ISR/FSR evolution.
    scale = sqrt(Q2Fac());
    process.scale( scale );

    // Loop over incoming and outgoing partons.
    int colOffset = process.lastColTag();
    for (int i = 1; i <= 2 + nFin; ++i) {

      // Read out particle info from SigmaProcess object.
      int id        = sigmaProcessPtr->id(i);
      int status    = (i <= 2) ? -21 : 23;
      int mother1   = (i <= 2) ? i : m_M1 ;
      int mother2   = (i <= 2) ? 0 : m_M2 ;
      int daughter1 = (i <= 2) ? m_D1 : 0;
      int daughter2 = (i <= 2) ? m_D2 : 0;
      int col       = sigmaProcessPtr->col(i);
      if      (col > 0) col += colOffset;
      else if (col < 0) col -= colOffset;
      int acol      = sigmaProcessPtr->acol(i);
      if      (acol > 0) acol += colOffset;
      else if (acol < 0) acol -= colOffset;

      // If extra photons in event record, offset the mother/daughter list.
      if ( beamAhasResGamma || beamBhasResGamma
         || (beamAgammaMode == 2 && beamBgammaMode == 0)
         || (beamAgammaMode == 0 && beamBgammaMode == 2) ) {
        if (mother1 > 0)   mother1   += nOffsetGamma;
        if (mother2 > 0)   mother2   += nOffsetGamma;
        if (daughter1 > 0) daughter1 += nOffsetGamma;
        if (daughter2 > 0) daughter2 += nOffsetGamma;
      }

      // Append to process record.
      int iNow = process.append( id, status, mother1, mother2,
        daughter1, daughter2, col, acol, phaseSpacePtr->p(i),
        phaseSpacePtr->m(i), scale);

      // NOAM: If there is an intermediate state, insert the it in
      // the process record after the two incoming particles.
      if (i == 2 && idRes != 0) {

        // Sign of intermediate state: go by charge.
        if (particleDataPtr->hasAnti(idRes)
          && process[3].chargeType() + process[4].chargeType() < 0)
          idRes *= -1;

        // The colour configuration of the intermediate state has to be
        // resolved separately.
        col         = 0;
        acol        = 0;
        int m_col1  = sigmaProcessPtr->col(1);
        int m_acol1 = sigmaProcessPtr->acol(1);
        int m_col2  = sigmaProcessPtr->col(2);
        int m_acol2 = sigmaProcessPtr->acol(2);
        if (m_col1 == m_acol2 && m_col2 != m_acol1) {
          col       = m_col2;
          acol       = m_acol1;
        } else if (m_col2 == m_acol1 && m_col1 != m_acol2) {
          col        = m_col1;
          acol       = m_acol2;
        }
        if      ( col > 0)  col += colOffset;
        else if ( col < 0)  col -= colOffset;
        if      (acol > 0) acol += colOffset;
        else if (acol < 0) acol -= colOffset;

        // Insert the intermediate state into the event record.
        Vec4 pIntMed = phaseSpacePtr->p(1) + phaseSpacePtr->p(2);
        process.append( idRes, -22, 3, 4,  6, 5 + nFin, col, acol,
          pIntMed, pIntMed.mCalc(), scale);
      }

      // Pick lifetime where relevant, else not.
      if (process[iNow].tau0() > 0.) process[iNow].tau(
        process[iNow].tau0() * rndmPtr->exp() );
    }
  }

  // Insert the outgoing particles - unresolved processes.
  else if (!isLHA) {
    int id3     = sigmaProcessPtr->id(3);
    int status3 = (id3 == process[1].id()) ? 14 : 15;
    process.append( id3, status3, 1, 0, 0, 0, 0, 0,
      phaseSpacePtr->p(3), phaseSpacePtr->m(3));
    int id4     = sigmaProcessPtr->id(4);
    int status4 = (id4 == process[2].id()) ? 14 : 15;
    process.append( id4, status4, 2, 0, 0, 0, 0, 0,
      phaseSpacePtr->p(4), phaseSpacePtr->m(4));

    // For central diffraction: two scattered protons inserted so far,
    // but modify mothers, add also centrally-produced hadronic system.
    if (isDiffC) {
      process[3].mothers( 1, 2);
      process[4].mothers( 1, 2);
      int id5     = sigmaProcessPtr->id(5);
      int status5 = 15;
      process.append( id5, status5, 1, 2, 0, 0, 0, 0,
        phaseSpacePtr->p(5), phaseSpacePtr->m(5));
    }
  }

  // Insert the outgoing particles - Les Houches Accord processes.
  else {

    // Since LHA partons may be out of order, determine correct one.
    // (Recall that zeroth particle is empty.)
    vector<int> newPos;
    newPos.reserve(lhaUpPtr->sizePart());
    newPos.push_back(0);
    for (int iNew = 0; iNew < lhaUpPtr->sizePart(); ++iNew) {
      // For iNew == 0 look for the two incoming partons, then for
      // partons having them as mothers, and so on layer by layer.
      for (int i = 1; i < lhaUpPtr->sizePart(); ++i)
        if (lhaUpPtr->mother1(i) == newPos[iNew]) newPos.push_back(i);
      if (int(newPos.size()) <= iNew) break;
    }

    // Find scale from which to begin MPI/ISR/FSR evolution.
    scalup = lhaUpPtr->scale();
    scale  = scalup;
    double scalePr = (scale < 0.) ? sqrt(Q2Fac()) : scale;
    process.scale( scalePr);

    // Copy over info from LHA event to process, in proper order.
    vector<int> iFinal;
    for (int i = 1; i < lhaUpPtr->sizePart(); ++i) {
      int iOld = newPos[i];
      int id = lhaUpPtr->id(iOld);
      if (i == 1 && abs(id) == idRenameBeams) id = 16;
      if (i == 2 && abs(id) == idRenameBeams) id = -16;
      int idAbs = abs(id);

      // Translate from LHA status codes.
      int lhaStatus =  lhaUpPtr->status(iOld);
      int status = -21;
      if (lhaStatus == 2 || lhaStatus == 3) status = -22;
      if (lhaStatus == 1) status = 23;

      // Find where mothers have been moved by reordering.
      int mother1Old = lhaUpPtr->mother1(iOld);
      int mother2Old = lhaUpPtr->mother2(iOld);
      int mother1 = 0;
      int mother2 = 0;
      for (int im = 1; im < i; ++im) {
        if (mother1Old == newPos[im]) mother1 = im + 2;
        if (mother2Old == newPos[im]) mother2 = im + 2;
      }
      if (i <= 2) mother1 = i;

      // Ensure that second mother = 0 except for bona fide carbon copies.
      if (mother1 > 0 && mother2 == mother1) {
        int sister1 = process[mother1].daughter1();
        int sister2 = process[mother1].daughter2();
        if (sister2 != sister1 && sister2 != 0) mother2 = 0;
      }

      // Find daughters and where they have been moved by reordering.
      // (Values shifted two steps to account for inserted beams.)
      int daughter1 = 0;
      int daughter2 = 0;
      for (int im = i + 1; im < lhaUpPtr->sizePart(); ++im) {
        if (lhaUpPtr->mother1(newPos[im]) == iOld
          || lhaUpPtr->mother2(newPos[im]) == iOld) {
          if (daughter1 == 0 || im + 2 < daughter1) daughter1 = im + 2;
          if (daughter2 == 0 || im + 2 > daughter2) daughter2 = im + 2;
        }
      }
      // For 2 -> 1 hard scatterings reset second daughter to 0.
      if (daughter2 == daughter1) daughter2 = 0;

      // Colour trivial, except reset irrelevant colour indices.
      int colType = particleDataPtr->colType(id);
      int col1   = (colType == 1 || colType == 2 || abs(colType) == 3)
                 ? lhaUpPtr->col1(iOld) : 0;
      int col2   = (colType == -1 || colType == 2 || abs(colType) == 3)
                 ?  lhaUpPtr->col2(iOld) : 0;

      // Copy momentum and ensure consistent (E, p, m) set.
      double px  = lhaUpPtr->px(iOld);
      double py  = lhaUpPtr->py(iOld);
      double pz  = lhaUpPtr->pz(iOld);
      double e   = lhaUpPtr->e(iOld);
      double m   = lhaUpPtr->m(iOld);
      if (mRecalculate > 0. && m > mRecalculate)
        m = sqrtpos( e*e - px*px - py*py - pz*pz);
      else e = sqrt( m*m + px*px + py*py + pz*pz);

      // Polarization.
      double pol = lhaUpPtr->spin(iOld);

      // Allow scale setting for generic partons.
      double scaleShow = lhaUpPtr->scale(iOld);

      // For resonance decay products use resonance mass as scale.
      double scaleNow = scalePr;
      if (mother1 > 4 && !useStrictLHEFscales) scaleNow = process[mother1].m();
      if (scaleShow >= 0.0) scaleNow = scaleShow;

      // Store Les Houches Accord partons.
      int iNow = process.append( id, status, mother1, mother2, daughter1,
        daughter2, col1, col2, Vec4(px, py, pz, e), m, scaleNow, pol);

      // Check if need to store lifetime.
      double tau = lhaUpPtr->tau(iOld);
      if ( (setLifetime == 1 && idAbs == 15) || setLifetime == 2)
         tau = process[iNow].tau0() * rndmPtr->exp();
      if (tau > 0.) process[iNow].tau(tau);

      // Track outgoing final particles.
      if (status == 23) iFinal.push_back( iNow);
    }

    // Second pass to catch vanishing final lepton and quark masses.
    int iFinalSz = iFinal.size();
    for (int iF = 0; iF < iFinalSz; ++iF) {
      int iMod = iFinal[iF];
      int iQLmod  = 0;
      double mOld = process[iMod].m();
      for (int iQL = 1; iQL < 9; ++iQL)
      if (process[iMod].idAbs() == idNewM[iQL]) {
        if ( iQL < 6 && setQuarkMass > 0 && (iQL == 4 || iQL == 5
          || setQuarkMass == 2) && (mOld < 0.5 * mNewM[iQL]
          || mOld > 1.5 * mNewM[iQL]) ) iQLmod = iQL;
        if ( iQL > 5 && setLeptonMass > 0 && (setLeptonMass == 2
          || mOld < 0.9 * mNewM[iQL] || mOld > 1.1 * mNewM[iQL]) )
          iQLmod = iQL;
      }
      if (iQLmod == 0) continue;
      double mNew = mNewM[iQLmod];

      // Find partner to exchange energy and momentum with: general.
      int iRec = 0;
      if (iFinalSz == 2) iRec = iFinal[1 - iF];
      else if (process[iMod].mother1() > 2) {
        int iMother = process[iMod].mother1();
        int iMDau1  = process[iMother].daughter1();
        int iMDau2  = process[iMother].daughter2();
        if (iMDau2 == iMDau1 + 1 && iMod == iMDau1) iRec = iMDau2;
        if (iMDau2 == iMDau1 + 1 && iMod == iMDau2) iRec = iMDau1;
      }

      // Find partner to exchange energy and momentum with: lepton.
      if (iRec == 0 && iQLmod > 5) {
        for (int iR = 0; iR < iFinalSz; ++iR) if (iR != iF) {
          int iRtmp = iFinal[iR];
          if (process[iRtmp].idAbs() == idNewM[iQLmod] + 1
            && process[iRtmp].id() * process[iMod].id() < 0) iRec = iRtmp;
        }
      }
      if (iRec == 0 && iQLmod > 5) {
        for (int iR = 0; iR < iFinalSz; ++iR) if (iR != iF) {
          int iRtmp = iFinal[iR];
          if (process[iRtmp].id() == -process[iMod].id()) iRec = iRtmp;
        }
      }

      // Find partner to exchange energy and momentum with: quark.
      if (iRec == 0 && iQLmod < 6 && process[iMod].col() != 0) {
        for (int iR = 0; iR < iFinalSz; ++iR) if (iR != iF) {
          int iRtmp = iFinal[iR];
          if (process[iRtmp].acol() == process[iMod].col()) iRec = iRtmp;
        }
      }
      if (iRec == 0 && iQLmod < 6 && process[iMod].acol() != 0) {
        for (int iR = 0; iR < iFinalSz; ++iR) if (iR != iF) {
          int iRtmp = iFinal[iR];
          if (process[iRtmp].col() == process[iMod].acol()) iRec = iRtmp;
        }
      }

      // Find partner to exchange energy and momentum with: largest mass.
      if (iRec == 0) {
        double mMax = 0.;
        for (int iR = 0; iR < iFinalSz; ++iR) if (iR != iF) {
          int iRtmp   = iFinal[iR];
          double mTmp = m( process[iMod], process[iRtmp]) - process[iRtmp].m();
          if (mTmp > mMax) { iRec = iRtmp; mMax = mTmp;}
        }
      }

      // Carry out exchange of energy and momentum, if possible.
      bool doneShift = false;
      Vec4 pMod   = process[iMod].p();
      if (iRec != 0) {
        Vec4 pRecOld = process[iRec].p();
        Vec4 pRecNew = pRecOld;
        double mRec  = process[iRec].m();
        if ( pShift( pMod, pRecNew, mNew, mRec) ) {
          process[iMod].p( pMod);
          process[iMod].m( mNew);
          process[iRec].p( pRecNew);
          // Recursive boost of recoiler decay products, if necessary.
          if (!process[iRec].isFinal()) {
            vector<int> iDauRec = process[iRec].daughterListRecursive();
            RotBstMatrix Mdau;
            Mdau.bst( pRecOld, pRecNew);
            for (int iD = 0; iD < int(iDauRec.size()); ++iD)
              process[iDauRec[iD]].rotbst( Mdau);
          }
          doneShift = true;
        }
      }

      // Emergency solution: just add energy as needed, => new incoming.
      if (!doneShift) {
        pMod.e( sqrtpos( pMod.pAbs2() + mNew * mNew) );
        process[iMod].p( pMod);
        process[iMod].m( mNew);
        infoPtr->errorMsg("Warning in ProcessContainer::constructProcess: "
          "unsuitable recoiler found");
      }
    }

    // Reassign momenta and masses for incoming partons.
    if (matchInOut) {
      Vec4 pSumOut;
      for (int iF = 0; iF < iFinalSz; ++iF)
        pSumOut += process[iFinal[iF]].p();
      double e1 = 0.5 * (pSumOut.e() + pSumOut.pz());
      double e2 = 0.5 * (pSumOut.e() - pSumOut.pz());
      process[3].pz( e1);
      process[3].e(  e1);
      process[3].m(  0.);
      process[4].pz(-e2);
      process[4].e(  e2);
      process[4].m(  0.);
      if (max(e1, e2) > 0.500001 * process[0].e()) {
        infoPtr->errorMsg("Error in ProcessContainer::constructProcess: "
          "setting mass failed");
        return false;
      }
    }
  }

  // Loop through decay chains and set secondary vertices when needed.
  for (int i = 3; i < process.size(); ++i) {
    int iMother  = process[i].mother1();

    // If sister to already assigned vertex then assign same.
    if ( process[i - 1].mother1() == iMother && process[i - 1].hasVertex() )
      process[i].vProd( process[i - 1].vProd() );

    // Else if mother already has vertex and/or lifetime then assign.
    else if ( process[iMother].hasVertex() || process[iMother].tau() > 0.)
      process[i].vProd( process[iMother].vDec() );
  }

  // Further info on process. Reset quantities that may or may not be known.
  int    id1Now  = process[3 + nOffsetGamma].id();
  int    id2Now  = process[4 + nOffsetGamma].id();
  int    id1pdf  = 0;
  int    id2pdf  = 0;
  double x1pdf   = 0.;
  double x2pdf   = 0.;
  double pdf1    = 0.;
  double pdf2    = 0.;
  double tHat    = 0.;
  double uHat    = 0.;
  double pTHatL  = 0.;
  double m3      = 0.;
  double m4      = 0.;
  double theta   = 0.;
  double phi     = 0.;
  double x1Now, x2Now, Q2FacNow, alphaEM, alphaS, Q2Ren, sHat;

  // Internally generated and stored information.
  if (!isLHA) {
    id1pdf       = id1Now;
    id2pdf       = id2Now;
    x1Now        = phaseSpacePtr->x1();
    x2Now        = phaseSpacePtr->x2();
    x1pdf        = x1Now;
    x2pdf        = x2Now;
    pdf1         = sigmaProcessPtr->pdf1();
    pdf2         = sigmaProcessPtr->pdf2();
    Q2FacNow     = sigmaProcessPtr->Q2Fac();
    alphaEM      = sigmaProcessPtr->alphaEMRen();
    alphaS       = sigmaProcessPtr->alphaSRen();
    Q2Ren        = sigmaProcessPtr->Q2Ren();
    sHat         = phaseSpacePtr->sHat();
    tHat         = phaseSpacePtr->tHat();
    uHat         = phaseSpacePtr->uHat();
    pTHatL       = phaseSpacePtr->pTHat();
    m3           = phaseSpacePtr->m(3);
    m4           = phaseSpacePtr->m(4);
    theta        = phaseSpacePtr->thetaHat();
    phi          = phaseSpacePtr->phiHat();
  }

  // Les Houches Accord process partly available, partly to be constructed.
  else {
    x1Now        = 2. * process[3].e() / infoPtr->eCM();
    x2Now        = 2. * process[4].e() / infoPtr->eCM();
    Q2FacNow     = (scale < 0.) ? sigmaProcessPtr->Q2Fac() : pow2(scale);
    alphaEM      = lhaUpPtr->alphaQED();
    if (alphaEM < 0.001) alphaEM = sigmaProcessPtr->alphaEMRen();
    alphaS       = lhaUpPtr->alphaQCD();
    if (alphaS  < 0.001) alphaS  = sigmaProcessPtr->alphaSRen();
    Q2Ren        = (scale < 0.) ? sigmaProcessPtr->Q2Ren() : pow2(scale);
    Vec4 pSum    = process[3].p() + process[4].p();
    sHat         = pSum * pSum;
    int nFinLH   = 0;
    double pTLH  = 0.;
    for (int i = 5; i < process.size(); ++i)
    if (process[i].mother1() == 3 && process[i].mother2() == 4) {
      ++nFinLH;
      pTLH      += process[i].pT();
    }
    if (nFinLH > 0) pTHatL = pTLH / nFinLH;

    // Read info on parton densities if provided.
    id1pdf       = lhaUpPtr->id1pdf();
    id2pdf       = lhaUpPtr->id2pdf();
    x1pdf        = lhaUpPtr->x1pdf();
    x2pdf        = lhaUpPtr->x2pdf();
    if (lhaUpPtr->pdfIsSet()) {
      pdf1       = lhaUpPtr->pdf1();
      pdf2       = lhaUpPtr->pdf2();
      Q2FacNow   = pow2(lhaUpPtr->scalePDF());
    }

    // Reconstruct kinematics of 2 -> 2 processes from momenta.
    if (nFin == 2) {
      Vec4 pDifT = process[3].p() - process[5].p();
      tHat       = pDifT * pDifT;
      Vec4 pDifU = process[3].p() - process[6].p();
      uHat       = pDifU * pDifU;
      pTHatL     = process[5].pT();
      m3         = process[5].m();
      m4         = process[6].m();
      Vec4 p5    = process[5].p();
      p5.bstback(pSum);
      theta      = p5.theta();
      phi        = process[5].phi();
    }
  }

  // Store information.
  if (isHardest) {
    infoPtr->setPDFalpha( 0, id1pdf, id2pdf, x1pdf, x2pdf, pdf1, pdf2,
      Q2FacNow, alphaEM, alphaS, Q2Ren, scalup);
    infoPtr->setKin( 0, id1Now, id2Now, x1Now, x2Now, sHat, tHat, uHat,
      pTHatL, m3, m4, theta, phi);
  }
  infoPtr->setTypeMPI( code(), pTHatL);

  // For Les Houches event store subprocess classification.
  if (isLHA) {
    int codeSub  = lhaUpPtr->idProcess();
    ostringstream nameSub;
    nameSub << "user process " << codeSub;
    infoPtr->setSubType( 0, nameSub.str(), codeSub, nFin);
  }

  // Done.
  return true;

}

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

// Give the hard resonance decay chain from Les Houches input.
// Note: mildly modified copy of some constructProcess code; to streamline.

bool ProcessContainer::constructDecays( Event& process) {

  // Let hard process record begin with the event as a whole.
  process.clear();
  process.append( 90, -11, 0, 0, 0, 0, 0, 0, Vec4(0., 0., 0., 0.), 0., 0. );

  // Since LHA partons may be out of order, determine correct one.
  // (Recall that zeroth particle is empty.)
  vector<int> newPos;
  newPos.reserve(lhaUpPtr->sizePart());
  newPos.push_back(0);
  for (int iNew = 0; iNew < lhaUpPtr->sizePart(); ++iNew) {
    // For iNew == 0 look for the two incoming partons, then for
    // partons having them as mothers, and so on layer by layer.
    for (int i = 1; i < lhaUpPtr->sizePart(); ++i)
      if (lhaUpPtr->mother1(i) == newPos[iNew]) newPos.push_back(i);
    if (int(newPos.size()) <= iNew) break;
  }

  // Find scale from which to begin MPI/ISR/FSR evolution.
  double scale = lhaUpPtr->scale();
  process.scale( scale);
  Vec4 pSum;

  // Copy over info from LHA event to process, in proper order.
  for (int i = 1; i < lhaUpPtr->sizePart(); ++i) {
    int iOld = newPos[i];
    int id = lhaUpPtr->id(iOld);

    // Translate from LHA status codes.
    int lhaStatus =  lhaUpPtr->status(iOld);
    int status = -21;
    if (lhaStatus == 2 || lhaStatus == 3) status = -22;
    if (lhaStatus == 1) status = 23;

    // Find where mothers have been moved by reordering.
    int mother1Old = lhaUpPtr->mother1(iOld);
    int mother2Old = lhaUpPtr->mother2(iOld);
    int mother1 = 0;
    int mother2 = 0;
    for (int im = 1; im < i; ++im) {
      if (mother1Old == newPos[im]) mother1 = im;
      if (mother2Old == newPos[im]) mother2 = im;
    }

    // Ensure that second mother = 0 except for bona fide carbon copies.
    if (mother1 > 0 && mother2 == mother1) {
      int sister1 = process[mother1].daughter1();
      int sister2 = process[mother1].daughter2();
      if (sister2 != sister1 && sister2 != 0) mother2 = 0;
    }

    // Find daughters and where they have been moved by reordering.
    int daughter1 = 0;
    int daughter2 = 0;
    for (int im = i + 1; im < lhaUpPtr->sizePart(); ++im) {
      if (lhaUpPtr->mother1(newPos[im]) == iOld
        || lhaUpPtr->mother2(newPos[im]) == iOld) {
        if (daughter1 == 0 || im < daughter1) daughter1 = im;
        if (daughter2 == 0 || im > daughter2) daughter2 = im;
      }
    }
    // For 2 -> 1 hard scatterings reset second daughter to 0.
    if (daughter2 == daughter1) daughter2 = 0;

    // Colour trivial, except reset irrelevant colour indices.
    int colType = particleDataPtr->colType(id);
    int col1   = (colType == 1 || colType == 2 || abs(colType) == 3)
               ? lhaUpPtr->col1(iOld) : 0;
    int col2   = (colType == -1 || colType == 2 || abs(colType) == 3)
               ?  lhaUpPtr->col2(iOld) : 0;

    // Momentum trivial.
    double px  = lhaUpPtr->px(iOld);
    double py  = lhaUpPtr->py(iOld);
    double pz  = lhaUpPtr->pz(iOld);
    double e   = lhaUpPtr->e(iOld);
    double m   = lhaUpPtr->m(iOld);
    if (status > 0) pSum += Vec4( px, py, pz, e);

    // Polarization.
    double pol = lhaUpPtr->spin(iOld);

    // For resonance decay products use resonance mass as scale.
    double scaleNow = scale;
    if (mother1 > 0) scaleNow = process[mother1].m();

    // Store Les Houches Accord partons.
    int iNow = process.append( id, status, mother1, mother2, daughter1,
      daughter2, col1, col2, Vec4(px, py, pz, e), m, scaleNow, pol);

    // Check if need to store lifetime.
    double tau = lhaUpPtr->tau(iOld);
    if ( (setLifetime == 1 && abs(id) == 15) || setLifetime == 2)
       tau = process[iNow].tau0() * rndmPtr->exp();
    if (tau > 0.) process[iNow].tau(tau);
  }

  // Update four-momentum of system as a whole.
  process[0].p( pSum);
  process[0].m( pSum.mCalc());

  // Loop through decay chains and set secondary vertices when needed.
  for (int i = 1; i < process.size(); ++i) {
    int iMother  = process[i].mother1();

    // If sister to already assigned vertex then assign same.
    if ( process[i - 1].mother1() == iMother && process[i - 1].hasVertex()
      && i > 1) process[i].vProd( process[i - 1].vProd() );

    // Else if mother already has vertex and/or lifetime then assign.
    else if ( process[iMother].hasVertex() || process[iMother].tau() > 0.)
      process[i].vProd( process[iMother].vDec() );
  }

  // Done.
  return true;

}

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

// Handle resonance decays.

bool ProcessContainer::decayResonances( Event& process) {

  // Save current event-record size and status codes.
  process.saveSize();
  vector<int> statusSave( process.size());
  for (int i = 0; i < process.size(); ++i)
    statusSave[i] = process[i].status();
  bool physical    = true;
  bool newChain    = false;
  bool newFlavours = false;

  // Do loop over user veto.
  do {

    // Do sequential chain of uncorrelated isotropic decays.
    do {
      physical = resDecaysPtr->next( process);
      if (!physical) return false;

      // Check whether flavours should be correlated.
      // (Currently only relevant for f fbar -> gamma*/Z0 gamma*/Z0.)
      newFlavours = ( sigmaProcessPtr->weightDecayFlav( process)
                    < rndmPtr->flat() );

      // Reset the decay chains if have to redo.
      if (newFlavours) {
        process.restoreSize();
        for (int i = 0; i < process.size(); ++i)
          process[i].status( statusSave[i]);
      }

    // Loop back where required to generate new decays with new flavours.
    } while (newFlavours);

    // Correct to nonisotropic decays.
    phaseSpacePtr->decayKinematics( process);

    // Optionally user hooks check/veto on decay chain.
    if (canVetoResDecay)
      newChain = userHooksPtr->doVetoResonanceDecays( process);

    // Reset the decay chains if have to redo.
    if (newChain) {
      process.restoreSize();
      for (int i = 0; i < process.size(); ++i)
        process[i].status( statusSave[i]);
    }

  // Loop back where required to generate new decay chain.
  } while(newChain);

  // Done.
  return true;

}

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

// Reset event generation statistics; but NOT maximum of cross section.

void ProcessContainer::reset() {

  nTry      = 0;
  nSel      = 0;
  nAcc      = 0;
  nTryStat  = 0;
  sigmaSum  = 0.;
  sigma2Sum = 0.;
  sigmaNeg  = 0.;
  sigmaAvg  = 0.;
  sigmaFin  = 0.;
  deltaFin  = 0.;
  wtAccSum  = 0.;
  sigmaTemp = 0.;
  sigma2Temp= 0.;

}

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

// Estimate integrated cross section and its uncertainty.

void ProcessContainer::sigmaDelta() {

  // Initial values. No analysis meaningful unless accepted events.
  nTryStat = nTry;
  sigmaAvg = 0.;
  sigmaFin = 0.;
  deltaFin = 0.;
  if (nAcc == 0) return;

  // Only update once an event is accepted, as is needed if the event weight
  // can still change by a finite amount due to reweighting.
  double wgtNow = infoPtr->weight();
  // infoPtr->weight() performs coversion to pb if lhaStratAbs = 4
  if (lhaStrat    == 0) wgtNow  = sigmaTemp;
  if (lhaStratAbs == 3) wgtNow *= sigmaTemp;
  if (lhaStratAbs == 4) wgtNow /= 1e9;
  sigmaSum += wgtNow;

  double wgtNow2 = 1.0;
  if (lhaStrat    == 0) wgtNow2 = sigma2Temp;
  if (lhaStratAbs == 3) wgtNow2 = pow2(wgtNow)*sigma2Temp;
  if (lhaStratAbs == 4) wgtNow2 = pow2(wgtNow/1e9);
  sigma2Sum += wgtNow2;

  sigmaTemp  = 0.0;
  sigma2Temp = 0.0;

  // Average value. No error analysis unless at least two events.
  double nTryInv  = 1. / nTry;
  double nSelInv  = 1. / nSel;
  double nAccInv  = 1. / nAcc;
  sigmaAvg        = sigmaSum * nTryInv;
  double fracAcc  = nAcc * nSelInv;

  // If Pythia should not perform the unweighting, always simply add accepted
  // event weights.
  if (lhaStratAbs >= 3 ) fracAcc = 1.;

  sigmaFin        = sigmaAvg * fracAcc;
  deltaFin        = sigmaFin;

  if (nAcc == 1) return;

  // Estimated error. Quadratic sum of cross section term and
  // binomial from accept/reject step.
  double delta2Sig   = (lhaStratAbs != 3)
    ? (sigma2Sum * nTryInv - pow2(sigmaAvg)) * nTryInv / pow2(sigmaAvg)
    : pow2( lhaUpPtr->xErrSum() / lhaUpPtr->xSecSum());
  double delta2Veto  = (nSel - nAcc) * nAccInv * nSelInv;
  double delta2Sum   = delta2Sig + delta2Veto;
  deltaFin           = sqrtpos(delta2Sum) * sigmaFin;

}

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

// SetupContainer class.
// Turns list of user-desired processes into a vector of containers.

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

// Main routine to initialize list of processes.

bool SetupContainers::init(vector<ProcessContainer*>& containerPtrs,
       Info *infoPtr, Settings& settings, ParticleData* particleDataPtr,
       Couplings* couplings) {

  // Reset process list, if filled in previous subrun.
  if (containerPtrs.size() > 0) {
    for (int i = 0; i < int(containerPtrs.size()); ++i)
      delete containerPtrs[i];
    containerPtrs.clear();
  }
  SigmaProcess* sigmaPtr;

  // Set up requested objects for soft QCD processes.
  bool softQCD = settings.flag("SoftQCD:all");
  bool inelastic = settings.flag("SoftQCD:inelastic");
  if (softQCD || inelastic || settings.flag("SoftQCD:nonDiffractive")) {
    sigmaPtr = new Sigma0nonDiffractive;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (softQCD || settings.flag("SoftQCD:elastic")) {
    sigmaPtr = new Sigma0AB2AB;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (softQCD || inelastic || settings.flag("SoftQCD:singleDiffractive")) {
    sigmaPtr = new Sigma0AB2XB;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma0AB2AX;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (softQCD || inelastic || settings.flag("SoftQCD:doubleDiffractive")) {
    sigmaPtr = new Sigma0AB2XX;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (softQCD || inelastic || settings.flag("SoftQCD:centralDiffractive")) {
    sigmaPtr = new Sigma0AB2AXB;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for hard QCD processes.
  bool hardQCD = settings.flag("HardQCD:all");
  if (hardQCD || settings.flag("HardQCD:gg2gg")) {
    sigmaPtr = new Sigma2gg2gg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD || settings.flag("HardQCD:gg2qqbar")) {
    sigmaPtr = new Sigma2gg2qqbar;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD || settings.flag("HardQCD:qg2qg")) {
    sigmaPtr = new Sigma2qg2qg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD || settings.flag("HardQCD:qq2qq")) {
    sigmaPtr = new Sigma2qq2qq;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD || settings.flag("HardQCD:qqbar2gg")) {
    sigmaPtr = new Sigma2qqbar2gg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD || settings.flag("HardQCD:qqbar2qqbarNew")) {
    sigmaPtr = new Sigma2qqbar2qqbarNew;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for c cbar and b bbar, also hard QCD.
  bool hardccbar = settings.flag("HardQCD:hardccbar");
  if (hardQCD || hardccbar || settings.flag("HardQCD:gg2ccbar")) {
    sigmaPtr = new Sigma2gg2QQbar(4, 121);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD || hardccbar || settings.flag("HardQCD:qqbar2ccbar")) {
    sigmaPtr = new Sigma2qqbar2QQbar(4, 122);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  bool hardbbbar = settings.flag("HardQCD:hardbbbar");
  if (hardQCD || hardbbbar || settings.flag("HardQCD:gg2bbbar")) {
    sigmaPtr = new Sigma2gg2QQbar(5, 123);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD || hardbbbar || settings.flag("HardQCD:qqbar2bbbar")) {
    sigmaPtr = new Sigma2qqbar2QQbar(5, 124);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for hard QCD 2 -> 3 processes.
  bool hardQCD3parton = settings.flag("HardQCD:3parton");
  if (hardQCD3parton || settings.flag("HardQCD:gg2ggg")) {
    sigmaPtr = new Sigma3gg2ggg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qqbar2ggg")) {
    sigmaPtr = new Sigma3qqbar2ggg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qg2qgg")) {
    sigmaPtr = new Sigma3qg2qgg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qq2qqgDiff")) {
    sigmaPtr = new Sigma3qq2qqgDiff;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qq2qqgSame")) {
    sigmaPtr = new Sigma3qq2qqgSame;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qqbar2qqbargDiff")) {
    sigmaPtr = new Sigma3qqbar2qqbargDiff;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qqbar2qqbargSame")) {
    sigmaPtr = new Sigma3qqbar2qqbargSame;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:gg2qqbarg")) {
    sigmaPtr = new Sigma3gg2qqbarg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qg2qqqbarDiff")) {
    sigmaPtr = new Sigma3qg2qqqbarDiff;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hardQCD3parton || settings.flag("HardQCD:qg2qqqbarSame")) {
    sigmaPtr = new Sigma3qg2qqqbarSame;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for prompt photon processes.
  bool promptPhotons = settings.flag("PromptPhoton:all");
  if (promptPhotons
    || settings.flag("PromptPhoton:qg2qgamma")) {
    sigmaPtr = new Sigma2qg2qgamma;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (promptPhotons
    || settings.flag("PromptPhoton:qqbar2ggamma")) {
    sigmaPtr = new Sigma2qqbar2ggamma;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (promptPhotons
    || settings.flag("PromptPhoton:gg2ggamma")) {
    sigmaPtr = new Sigma2gg2ggamma;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (promptPhotons
    || settings.flag("PromptPhoton:ffbar2gammagamma")) {
    sigmaPtr = new Sigma2ffbar2gammagamma;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (promptPhotons
    || settings.flag("PromptPhoton:gg2gammagamma")) {
    sigmaPtr = new Sigma2gg2gammagamma;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for weak gauge boson t-channel exchange.
  bool weakBosonExchanges = settings.flag("WeakBosonExchange:all");
  if (weakBosonExchanges
    || settings.flag("WeakBosonExchange:ff2ff(t:gmZ)")) {
    sigmaPtr = new Sigma2ff2fftgmZ;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonExchanges
    || settings.flag("WeakBosonExchange:ff2ff(t:W)")) {
    sigmaPtr = new Sigma2ff2fftW;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for weak gauge boson processes.
  bool weakSingleBosons = settings.flag("WeakSingleBoson:all");
  if (weakSingleBosons
    || settings.flag("WeakSingleBoson:ffbar2gmZ")) {
    sigmaPtr = new Sigma1ffbar2gmZ;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakSingleBosons
    || settings.flag("WeakSingleBoson:ffbar2W")) {
    sigmaPtr = new Sigma1ffbar2W;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested object for s-channel gamma exchange.
  // Subset of gamma*/Z0 above, intended for multiparton interactions.
  if (settings.flag("WeakSingleBoson:ffbar2ffbar(s:gm)")) {
    sigmaPtr = new Sigma2ffbar2ffbarsgm;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested object for s-channel gamma*/Z0 or W+- exchange.
  if (settings.flag("WeakSingleBoson:ffbar2ffbar(s:gmZ)")) {
    sigmaPtr = new Sigma2ffbar2ffbarsgmZ;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("WeakSingleBoson:ffbar2ffbar(s:W)")) {
    sigmaPtr = new Sigma2ffbar2ffbarsW;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for weak gauge boson pair processes.
  bool weakDoubleBosons = settings.flag("WeakDoubleBoson:all");
  if (weakDoubleBosons
    || settings.flag("WeakDoubleBoson:ffbar2gmZgmZ")) {
    sigmaPtr = new Sigma2ffbar2gmZgmZ;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakDoubleBosons
    || settings.flag("WeakDoubleBoson:ffbar2ZW")) {
    sigmaPtr = new Sigma2ffbar2ZW;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakDoubleBosons
    || settings.flag("WeakDoubleBoson:ffbar2WW")) {
    sigmaPtr = new Sigma2ffbar2WW;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for weak gauge boson + parton processes.
  bool weakBosonAndPartons = settings.flag("WeakBosonAndParton:all");
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:qqbar2gmZg")) {
    sigmaPtr = new Sigma2qqbar2gmZg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:qg2gmZq")) {
    sigmaPtr = new Sigma2qg2gmZq;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:ffbar2gmZgm")) {
    sigmaPtr = new Sigma2ffbar2gmZgm;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:fgm2gmZf")) {
    sigmaPtr = new Sigma2fgm2gmZf;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:qqbar2Wg")) {
    sigmaPtr = new Sigma2qqbar2Wg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:qg2Wq")) {
    sigmaPtr = new Sigma2qg2Wq;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:ffbar2Wgm")) {
    sigmaPtr = new Sigma2ffbar2Wgm;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (weakBosonAndPartons
    || settings.flag("WeakBosonAndParton:fgm2Wf")) {
    sigmaPtr = new Sigma2fgm2Wf;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for photon collision processes.
  bool photonCollisions = settings.flag("PhotonCollision:all");
  bool hasGamma      = settings.flag("PDF:lepton2gamma");
  int  photonMode    = settings.mode("Photon:ProcessType");
  bool initGmGm      = ( ( (photonMode == 4) || (photonMode == 0) )
                     && hasGamma ) || !hasGamma;
  if ( initGmGm && ( photonCollisions
    || settings.flag("PhotonCollision:gmgm2qqbar") ) ) {
    sigmaPtr = new Sigma2gmgm2ffbar(1, 261);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if ( initGmGm && ( photonCollisions
    || settings.flag("PhotonCollision:gmgm2ccbar") ) ) {
    sigmaPtr = new Sigma2gmgm2ffbar(4, 262);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if ( initGmGm && ( photonCollisions
    || settings.flag("PhotonCollision:gmgm2bbbar") ) ) {
    sigmaPtr = new Sigma2gmgm2ffbar(5, 263);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if ( initGmGm && ( photonCollisions
    || settings.flag("PhotonCollision:gmgm2ee") ) ) {
    sigmaPtr = new Sigma2gmgm2ffbar(11, 264);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if ( initGmGm && ( photonCollisions
    || settings.flag("PhotonCollision:gmgm2mumu") ) ) {
    sigmaPtr = new Sigma2gmgm2ffbar(13, 265);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if ( initGmGm && ( photonCollisions
    || settings.flag("PhotonCollision:gmgm2tautau") ) ) {
    sigmaPtr = new Sigma2gmgm2ffbar(15, 266);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Read in settings related to photon beams.
  bool photonParton  = settings.flag("PhotonParton:all");
  bool beamAisHadron = abs(infoPtr->idA()) > 100;
  bool beamBisHadron = abs(infoPtr->idB()) > 100;
  bool isGammaGamma  = infoPtr->idA() == 22 && infoPtr->idB() == 22;
  bool isGammaHadron = (infoPtr->idA() == 22 && beamBisHadron)
                    || (infoPtr->idB() == 22 && beamAisHadron);

  // Check if photon beams present.
  int idAbsA     = abs(infoPtr->idA());
  int idAbsB     = abs(infoPtr->idB());
  bool hasGammaA = hasGamma && (idAbsA == 11 || idAbsA == 13 || idAbsA == 15);
  bool hasGammaB = hasGamma && (idAbsB == 11 || idAbsB == 13 || idAbsB == 15);
  bool gammaA    = hasGammaA || idAbsA == 22;
  bool gammaB    = hasGammaB || idAbsB == 22;

  // Initialize necessary amount of processes with initiating photons.
  // When no photon beams initialize one instance (initGammaB == true).
  bool initGammaA = gammaA && (photonMode == 0 || photonMode == 3);
  bool initGammaB = !(isGammaGamma || hasGamma || isGammaHadron)
    || ( (photonMode == 0 || photonMode == 2) && gammaB);

  // Set up requested objects for photon-parton processes.
  // In case of unresolved photons insert process with correct initiator.
  if (photonParton || settings.flag("PhotonParton:ggm2qqbar")) {
    if ( initGammaB ) {
      sigmaPtr = new Sigma2ggm2qqbar(1, 271, "ggm");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if ( initGammaA ) {
      sigmaPtr = new Sigma2ggm2qqbar(1, 281, "gmg");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
  }
  if (photonParton || settings.flag("PhotonParton:ggm2ccbar")) {
    if ( initGammaB ) {
      sigmaPtr = new Sigma2ggm2qqbar(4, 272, "ggm");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if ( initGammaA) {
      sigmaPtr = new Sigma2ggm2qqbar(4, 282, "gmg");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
  }
  if (photonParton || settings.flag("PhotonParton:ggm2bbbar")) {
    if ( initGammaB ) {
      sigmaPtr = new Sigma2ggm2qqbar(5, 273, "ggm");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if ( initGammaA ) {
      sigmaPtr = new Sigma2ggm2qqbar(5, 283, "gmg");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
  }
  if (photonParton || settings.flag("PhotonParton:qgm2qg")) {
    if ( initGammaB ) {
      sigmaPtr = new Sigma2qgm2qg(274, "qgm");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if ( initGammaA ) {
      sigmaPtr = new Sigma2qgm2qg(284, "gmq");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
  }
  if (settings.flag("PhotonParton:qgm2qgm")) {
    if ( initGammaB ) {
      sigmaPtr = new Sigma2qgm2qgm(275, "qgm");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if ( initGammaA ) {
      sigmaPtr = new Sigma2qgm2qgm(285, "gmq");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
  }

  // Set up requested objects for onia production.
  charmonium = SigmaOniaSetup(infoPtr, &settings, particleDataPtr, 4);
  bottomonium = SigmaOniaSetup(infoPtr, &settings, particleDataPtr, 5);
  vector<SigmaProcess*> charmoniumSigmaPtrs, bottomoniumSigmaPtrs;
  charmonium.setupSigma2gg(charmoniumSigmaPtrs);
  charmonium.setupSigma2qg(charmoniumSigmaPtrs);
  charmonium.setupSigma2qq(charmoniumSigmaPtrs);
  charmonium.setupSigma2dbl(charmoniumSigmaPtrs);
  bottomonium.setupSigma2gg(bottomoniumSigmaPtrs);
  bottomonium.setupSigma2qg(bottomoniumSigmaPtrs);
  bottomonium.setupSigma2qq(bottomoniumSigmaPtrs);
  bottomonium.setupSigma2dbl(bottomoniumSigmaPtrs);
  for (unsigned int i = 0; i < charmoniumSigmaPtrs.size(); ++i)
    containerPtrs.push_back( new ProcessContainer(charmoniumSigmaPtrs[i]) );
  for (unsigned int i = 0; i < bottomoniumSigmaPtrs.size(); ++i)
    containerPtrs.push_back( new ProcessContainer(bottomoniumSigmaPtrs[i]) );

  // Set up requested objects for top production.
  bool tops = settings.flag("Top:all");
  if (tops || settings.flag("Top:gg2ttbar")) {
    sigmaPtr = new Sigma2gg2QQbar(6, 601);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tops || settings.flag("Top:qqbar2ttbar")) {
    sigmaPtr = new Sigma2qqbar2QQbar(6, 602);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tops || settings.flag("Top:qq2tq(t:W)")) {
    sigmaPtr = new Sigma2qq2QqtW(6, 603);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tops || settings.flag("Top:ffbar2ttbar(s:gmZ)")) {
    sigmaPtr = new Sigma2ffbar2FFbarsgmZ(6, 604);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tops || settings.flag("Top:ffbar2tqbar(s:W)")) {
    sigmaPtr = new Sigma2ffbar2FfbarsW(6, 0, 605);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tops || settings.flag("Top:gmgm2ttbar")) {
    sigmaPtr = new Sigma2gmgm2ffbar(6, 606);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for fourth-generation b' production
  bool bPrimes = settings.flag("FourthBottom:all");
  if (bPrimes || settings.flag("FourthBottom:gg2bPrimebPrimebar")) {
    sigmaPtr = new Sigma2gg2QQbar(7, 801);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (bPrimes || settings.flag("FourthBottom:qqbar2bPrimebPrimebar")) {
    sigmaPtr = new Sigma2qqbar2QQbar(7, 802);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (bPrimes || settings.flag("FourthBottom:qq2bPrimeq(t:W)")) {
    sigmaPtr = new Sigma2qq2QqtW(7, 803);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (bPrimes || settings.flag("FourthBottom:ffbar2bPrimebPrimebar(s:gmZ)")) {
    sigmaPtr = new Sigma2ffbar2FFbarsgmZ(7, 804);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (bPrimes || settings.flag("FourthBottom:ffbar2bPrimeqbar(s:W)")) {
    sigmaPtr = new Sigma2ffbar2FfbarsW(7, 0, 805);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (bPrimes || settings.flag("FourthBottom:ffbar2bPrimetbar(s:W)")) {
    sigmaPtr = new Sigma2ffbar2FfbarsW(7, 6, 806);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for fourth-generation t' production
  bool tPrimes = settings.flag("FourthTop:all");
  if (tPrimes || settings.flag("FourthTop:gg2tPrimetPrimebar")) {
    sigmaPtr = new Sigma2gg2QQbar(8, 821);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tPrimes || settings.flag("FourthTop:qqbar2tPrimetPrimebar")) {
    sigmaPtr = new Sigma2qqbar2QQbar(8, 822);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tPrimes || settings.flag("FourthTop:qq2tPrimeq(t:W)")) {
    sigmaPtr = new Sigma2qq2QqtW(8, 823);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tPrimes || settings.flag("FourthTop:ffbar2tPrimetPrimebar(s:gmZ)")) {
    sigmaPtr = new Sigma2ffbar2FFbarsgmZ(8, 824);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (tPrimes || settings.flag("FourthTop:ffbar2tPrimeqbar(s:W)")) {
    sigmaPtr = new Sigma2ffbar2FfbarsW(8, 0, 825);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for two different fourth-generation fermions.
  if (bPrimes || tPrimes
    || settings.flag("FourthPair:ffbar2tPrimebPrimebar(s:W)")) {
    sigmaPtr = new Sigma2ffbar2FfbarsW(8, 7, 841);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("FourthPair:ffbar2tauPrimenuPrimebar(s:W)")) {
    sigmaPtr = new Sigma2ffbar2FfbarsW(17, 18, 842);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Flag for global choice between SM and BSM Higgses.
  bool useBSMHiggses = settings.flag("Higgs:useBSM");

  // Set up requested objects for Standard-Model Higgs production.
  if (!useBSMHiggses) {
    bool HiggsesSM = settings.flag("HiggsSM:all");
    if (HiggsesSM || settings.flag("HiggsSM:ffbar2H")) {
      sigmaPtr = new Sigma1ffbar2H(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:gg2H")) {
     sigmaPtr = new Sigma1gg2H(0);
     containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:gmgm2H")) {
      sigmaPtr = new Sigma1gmgm2H(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:ffbar2HZ")) {
      sigmaPtr = new Sigma2ffbar2HZ(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:ffbar2HW")) {
      sigmaPtr = new Sigma2ffbar2HW(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:ff2Hff(t:ZZ)")) {
      sigmaPtr = new Sigma3ff2HfftZZ(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:ff2Hff(t:WW)")) {
      sigmaPtr = new Sigma3ff2HfftWW(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:gg2Httbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(6,0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesSM || settings.flag("HiggsSM:qqbar2Httbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(6,0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Further Standard-Model Higgs processes, not included in "all".
    if (settings.flag("HiggsSM:qg2Hq")) {
      sigmaPtr = new Sigma2qg2Hq(4,0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
      sigmaPtr = new Sigma2qg2Hq(5,0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsSM:gg2Hbbbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(5,0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsSM:qqbar2Hbbbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(5,0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsSM:gg2Hg(l:t)")) {
      sigmaPtr = new Sigma2gg2Hglt(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsSM:qg2Hq(l:t)")) {
      sigmaPtr = new Sigma2qg2Hqlt(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsSM:qqbar2Hg(l:t)")) {
      sigmaPtr = new Sigma2qqbar2Hglt(0);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
  }

  // Common switch for the group of Higgs production BSM.
  if (useBSMHiggses) {
    bool HiggsesBSM = settings.flag("HiggsBSM:all");

    // Set up requested objects for BSM H1 production.
    bool HiggsesH1 = settings.flag("HiggsBSM:allH1");
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:ffbar2H1")) {
      sigmaPtr = new Sigma1ffbar2H(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:gg2H1")) {
      sigmaPtr = new Sigma1gg2H(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:gmgm2H1")) {
      sigmaPtr = new Sigma1gmgm2H(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:ffbar2H1Z")) {
      sigmaPtr = new Sigma2ffbar2HZ(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:ffbar2H1W")) {
      sigmaPtr = new Sigma2ffbar2HW(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:ff2H1ff(t:ZZ)")) {
      sigmaPtr = new Sigma3ff2HfftZZ(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:ff2H1ff(t:WW)")) {
      sigmaPtr = new Sigma3ff2HfftWW(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:gg2H1ttbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(6,1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH1 || settings.flag("HiggsBSM:qqbar2H1ttbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(6,1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Further BSM H1 processes, not included in "all".
    if (settings.flag("HiggsBSM:qg2H1q")) {
      sigmaPtr = new Sigma2qg2Hq(4,1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
      sigmaPtr = new Sigma2qg2Hq(5,1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:gg2H1bbbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(5,1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qqbar2H1bbbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(5,1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:gg2H1g(l:t)")) {
      sigmaPtr = new Sigma2gg2Hglt(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qg2H1q(l:t)")) {
      sigmaPtr = new Sigma2qg2Hqlt(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qqbar2H1g(l:t)")) {
      sigmaPtr = new Sigma2qqbar2Hglt(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Set up requested objects for BSM H2 production.
    bool HiggsesH2 = settings.flag("HiggsBSM:allH2");
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:ffbar2H2")) {
      sigmaPtr = new Sigma1ffbar2H(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:gg2H2")) {
      sigmaPtr = new Sigma1gg2H(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:gmgm2H2")) {
      sigmaPtr = new Sigma1gmgm2H(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:ffbar2H2Z")) {
      sigmaPtr = new Sigma2ffbar2HZ(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:ffbar2H2W")) {
      sigmaPtr = new Sigma2ffbar2HW(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:ff2H2ff(t:ZZ)")) {
      sigmaPtr = new Sigma3ff2HfftZZ(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:ff2H2ff(t:WW)")) {
      sigmaPtr = new Sigma3ff2HfftWW(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:gg2H2ttbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(6,2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesH2 || settings.flag("HiggsBSM:qqbar2H2ttbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(6,2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Further BSM H2 processes, not included in "all".
   if (settings.flag("HiggsBSM:qg2H2q")) {
      sigmaPtr = new Sigma2qg2Hq(4,2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
      sigmaPtr = new Sigma2qg2Hq(5,2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:gg2H2bbbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(5,2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qqbar2H2bbbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(5,2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:gg2H2g(l:t)")) {
      sigmaPtr = new Sigma2gg2Hglt(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qg2H2q(l:t)")) {
      sigmaPtr = new Sigma2qg2Hqlt(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qqbar2H2g(l:t)")) {
      sigmaPtr = new Sigma2qqbar2Hglt(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Set up requested objects for BSM A3 production.
    bool HiggsesA3 = settings.flag("HiggsBSM:allA3");
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:ffbar2A3")) {
      sigmaPtr = new Sigma1ffbar2H(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:gg2A3")) {
      sigmaPtr = new Sigma1gg2H(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:gmgm2A3")) {
      sigmaPtr = new Sigma1gmgm2H(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:ffbar2A3Z")) {
      sigmaPtr = new Sigma2ffbar2HZ(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:ffbar2A3W")) {
      sigmaPtr = new Sigma2ffbar2HW(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:ff2A3ff(t:ZZ)")) {
      sigmaPtr = new Sigma3ff2HfftZZ(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:ff2A3ff(t:WW)")) {
      sigmaPtr = new Sigma3ff2HfftWW(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:gg2A3ttbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(6,3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesA3 || settings.flag("HiggsBSM:qqbar2A3ttbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(6,3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Further BSM A3 processes, not included in "all".
    if (settings.flag("HiggsBSM:qg2A3q")) {
      sigmaPtr = new Sigma2qg2Hq(4,3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
      sigmaPtr = new Sigma2qg2Hq(5,3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:gg2A3bbbar")) {
      sigmaPtr = new Sigma3gg2HQQbar(5,3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qqbar2A3bbbar")) {
      sigmaPtr = new Sigma3qqbar2HQQbar(5,3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:gg2A3g(l:t)")) {
      sigmaPtr = new Sigma2gg2Hglt(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qg2A3q(l:t)")) {
      sigmaPtr = new Sigma2qg2Hqlt(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (settings.flag("HiggsBSM:qqbar2A3g(l:t)")) {
      sigmaPtr = new Sigma2qqbar2Hglt(3);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Set up requested objects for Charged Higgs production
    bool HiggsesChg = settings.flag("HiggsBSM:allH+-");
    if (HiggsesBSM || HiggsesChg || settings.flag("HiggsBSM:ffbar2H+-")) {
      sigmaPtr = new Sigma1ffbar2Hchg;
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesChg || settings.flag("HiggsBSM:bg2H+-t")) {
      sigmaPtr = new Sigma2qg2Hchgq(6, 1062, "b g -> H+- t");
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }

    // Set up requested objects for Higgs pair-production
    bool HiggsesPairs = settings.flag("HiggsBSM:allHpair");
    if (HiggsesBSM || HiggsesPairs || settings.flag("HiggsBSM:ffbar2A3H1")) {
      sigmaPtr = new Sigma2ffbar2A3H12(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesPairs || settings.flag("HiggsBSM:ffbar2A3H2")) {
      sigmaPtr = new Sigma2ffbar2A3H12(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesPairs || settings.flag("HiggsBSM:ffbar2H+-H1")) {
      sigmaPtr = new Sigma2ffbar2HchgH12(1);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesPairs || settings.flag("HiggsBSM:ffbar2H+-H2")) {
      sigmaPtr = new Sigma2ffbar2HchgH12(2);
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
    if (HiggsesBSM || HiggsesPairs || settings.flag("HiggsBSM:ffbar2H+H-")) {
      sigmaPtr = new Sigma2ffbar2HposHneg();
      containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
    }
  }

  // Set up requested objects for SUSY pair processes.
  if (couplings->isSUSY) {
    CoupSUSY* coupSUSY = (CoupSUSY *) couplings;

    bool SUSYs = settings.flag("SUSY:all");
    bool nmssm = settings.flag("SLHA:NMSSM");

    // Preselected SUSY codes.
    setupIdVecs( settings);

    // MSSM: 4 neutralinos
    int nNeut = 4;
    if (nmssm) nNeut = 5;

    // Gluino-gluino
    if (SUSYs || settings.flag("SUSY:gg2gluinogluino")) {
      // Skip if outgoing codes not asked for
      if (allowIdVals( 1000021, 1000021)) {
        sigmaPtr = new Sigma2gg2gluinogluino();
        containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
      }
    }
    if (SUSYs || settings.flag("SUSY:qqbar2gluinogluino")) {
      // Skip if outgoing codes not asked for
      if (allowIdVals( 1000021, 1000021)) {
        sigmaPtr = new Sigma2qqbar2gluinogluino();
        containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
      }
    }

    // Gluino-squark
    if (SUSYs || settings.flag("SUSY:qg2squarkgluino")) {
      int iproc = 1202;
      for (int idx = 1; idx <= 6; ++idx) {
        for (int iso = 1; iso <= 2; ++iso) {
          iproc += 2;
          int id3 = iso + ((idx <= 3)
                  ? 1000000+2*(idx-1) : 2000000+2*(idx-4));
          int id4 = 1000021;
          // Skip if outgoing codes not asked for
          if (!allowIdVals( id3, id4)) continue;
          sigmaPtr = new Sigma2qg2squarkgluino(id3,iproc-1);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
          sigmaPtr = new Sigma2qg2squarkgluino(-id3,iproc);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
        }
      }
    }

    // Squark-antisquark (gg initiated)
    if (SUSYs || settings.flag("SUSY:gg2squarkantisquark")) {
      int iproc = 1214;
      for (int idx = 1; idx <= 6; ++idx) {
        for (int iso = 1; iso <= 2; ++iso) {
          iproc++;
          int id = iso + ((idx <= 3)
                 ? 1000000+2*(idx-1) : 2000000+2*(idx-4));
          // Skip if outgoing codes not asked for
          if (!allowIdVals( id, id)) continue;
          sigmaPtr = new Sigma2gg2squarkantisquark(id,iproc);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
        }
      }
    }

    // Squark-antisquark (qqbar initiated)
    if (SUSYs || settings.flag("SUSY:qqbar2squarkantisquark")) {
      int iproc = 1230;
      for (int idx = 1; idx <= 6; ++idx) {
        for (int iso = 1; iso <= 2; ++iso) {
          for (int jso = iso; jso >= 1; --jso) {
            for (int jdx = 1; jdx <= 6; ++jdx) {
              if (iso == jso && jdx < idx) continue;
              int id1 = iso + ((idx <= 3) ? 1000000+2*(idx-1)
                               : 2000000+2*(idx-4));
              int id2 = jso + ((jdx <= 3) ? 1000000+2*(jdx-1)
                               : 2000000+2*(jdx-4));
              // Update process number counter (for ~q~q, +2 if not self-conj)
              //if (iproc == 1302) iproc=1310;
              iproc++;
              if (iso == jso && id1 != id2) iproc++;
              // Skip if outgoing codes not asked for
              if (!allowIdVals( id1, id2)) continue;
              if (iso == jso && id1 != id2) {
                sigmaPtr = new Sigma2qqbar2squarkantisquark(id1,-id2,iproc-1);
                containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
                sigmaPtr = new Sigma2qqbar2squarkantisquark(id2,-id1,iproc);
                containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
              } else {
                sigmaPtr = new Sigma2qqbar2squarkantisquark(id1,-id2,iproc);
                containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
              }
            }
          }
        }
      }
    }

    // Squark-squark
    if (SUSYs || settings.flag("SUSY:qq2squarksquark")) {
      int iproc = 1350;
      for (int idx = 1; idx <= 6; ++idx) {
        for (int iso = 1; iso <= 2; ++iso) {
          for (int jso = iso; jso >= 1; jso--) {
            for (int jdx = 1; jdx <= 6; ++jdx) {
              if (iso == jso && jdx < idx) continue;
              iproc++;
              int id1 = iso + ((idx <= 3)
                      ? 1000000+2*(idx-1) : 2000000+2*(idx-4));
              int id2 = jso + ((jdx <= 3)
                      ? 1000000+2*(jdx-1) : 2000000+2*(jdx-4));
              // Skip if outgoing codes not asked for.
              if (!allowIdVals( id1, id2)) continue;
              sigmaPtr = new Sigma2qq2squarksquark(id1,id2,iproc);
              containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
            }
          }
        }
      }
    }

    // Neutralino + squark
    if (SUSYs || settings.flag("SUSY:qg2chi0squark")) {
      int iproc = 1430;
      for (int iNeut= 1; iNeut <= nNeut; iNeut++) {
        for (int idx = 1; idx <= 6; idx++) {
          bool isUp = false;
          for (int iso = 1; iso <= 2; iso++) {
            if (iso == 2) isUp = true;
            iproc++;
            int id3 = coupSUSY->idNeut(iNeut);
            int id4 = iso + ((idx <= 3)
                    ? 1000000+2*(idx-1) : 2000000+2*(idx-4));
            // Skip if outgoing codes not asked for
            if (!allowIdVals( id3, id4)) continue;
            sigmaPtr = new Sigma2qg2chi0squark(iNeut,idx,isUp,iproc);
            containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
          }
        }
      }
    }

    // Chargino + squark
    if (SUSYs || settings.flag("SUSY:qg2chi+-squark")) {
      int iproc = 1490;
      for (int iChar = 1; iChar <= 2; iChar++) {
        for (int idx = 1; idx <= 6; idx++) {
          bool isUp = false;
          for (int iso = 1; iso <= 2; iso++) {
            if (iso == 2) isUp = true;
            iproc++;
            int id3 = coupSUSY->idChar(iChar);
            int id4 = iso + ((idx <= 3)
                    ? 1000000+2*(idx-1) : 2000000+2*(idx-4));
            // Skip if outgoing codes not asked for
            if (!allowIdVals( id3, id4)) continue;
            sigmaPtr = new Sigma2qg2charsquark(iChar,idx,isUp,iproc);
            containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
          }
        }
      }
    }

    // Neutralino pairs
    if (SUSYs || settings.flag("SUSY:qqbar2chi0chi0")) {
      int iproc = 1550;
      for (int iNeut2 = 1; iNeut2 <= nNeut; iNeut2++) {
        for (int iNeut1 = 1; iNeut1 <= iNeut2; iNeut1++) {
          iproc++;
          // Skip if outgoing codes not asked for
          if (!allowIdVals( coupSUSY->idNeut(iNeut1),
            coupSUSY->idNeut(iNeut2) ) ) continue;
          sigmaPtr = new Sigma2qqbar2chi0chi0(iNeut1, iNeut2,iproc);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
        }
      }
    }

    // Neutralino-Chargino
    if (SUSYs || settings.flag("SUSY:qqbar2chi+-chi0")) {
      int iproc = 1570;
      for (int iNeut = 1; iNeut <= nNeut; iNeut++) {
        for (int iChar = 1; iChar <= 2; ++iChar) {
          iproc += 2;
          // Skip if outgoing codes not asked for
          if (!allowIdVals( coupSUSY->idNeut(iNeut),
            coupSUSY->idChar(iChar) ) ) continue;
          sigmaPtr = new Sigma2qqbar2charchi0( iChar, iNeut, iproc-1);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
          sigmaPtr = new Sigma2qqbar2charchi0(-iChar, iNeut, iproc);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
        }
      }
    }

    // Chargino-Chargino
    if (SUSYs || settings.flag("SUSY:qqbar2chi+chi-")) {
      int iproc = 1590;
      for (int i = 1; i <= 2; ++i) {
        for (int j = 1; j <= 2; ++j) {
          iproc++;
          // Skip if outgoing codes not asked for
          if (!allowIdVals( coupSUSY->idChar(i),
            coupSUSY->idChar(j) ) ) continue;
          sigmaPtr = new Sigma2qqbar2charchar( i,-j, iproc);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
        }
      }
    }

    // RPV squark production
    if(SUSYs || settings.flag("SUSY:qq2antisquark")) {
      for (int idx = 1; idx <= 6; ++idx) {
        for (int iso = 1; iso <= 2; ++iso) {
          int id1 = iso + ((idx <= 3) ? 1000000+2*(idx-1) : 2000000+2*(idx-4));
          // Skip if outgoing code not asked for
          if (!allowIdVals( id1, 0)) continue;
          sigmaPtr = new Sigma1qq2antisquark(id1);
          containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
        }
      }
    }

    // Neutralino-gluino
    if (SUSYs || settings.flag("SUSY:qqbar2chi0gluino")) {
      int iproc = 1600;
      for (int iNeut = 1; iNeut <= nNeut; iNeut++) {
        iproc++;
        // Skip if outgoing codes not asked for
        if (!allowIdVals( coupSUSY->idNeut(iNeut), 1000021)) continue;
        sigmaPtr = new Sigma2qqbar2chi0gluino(iNeut, iproc);
        containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
      }
    }

    // Chargino-Gluino
    if (SUSYs || settings.flag("SUSY:qqbar2chi+-gluino")) {
      int iproc = 1620;
      for (int iChar = 1; iChar <= 2; ++iChar) {
        iproc += 2;
        // Skip if outgoing codes not asked for
        if (!allowIdVals( coupSUSY->idChar(iChar), 1000021)) continue;
        sigmaPtr = new Sigma2qqbar2chargluino( iChar, iproc-1);
        containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
        sigmaPtr = new Sigma2qqbar2chargluino( -iChar, iproc);
        containerPtrs.push_back( new ProcessContainer(sigmaPtr) );

      }
    }

    // Slepton-antislepton (qqbar initiated); Currently no RH sneutrinos
    if (SUSYs || settings.flag("SUSY:qqbar2sleptonantislepton")) {
      int iproc = 1650;
      for (int idx = 1; idx <= 6; ++idx) {
        for (int iso = 1; iso <= 2; ++iso) {
          for (int jso = iso; jso >= 1; --jso) {
            for (int jdx = 1; jdx <= 6; ++jdx) {
              if (iso == jso && jdx < idx) continue;
              int id1 = iso + ((idx <= 3) ? 1000010+2*(idx-1)
                               : 2000010+2*(idx-4));
              int id2 = jso + ((jdx <= 3) ? 1000010+2*(jdx-1)
                               : 2000010+2*(jdx-4));
              // Update process number counter
              iproc++;
              if (iso == jso && id1 != id2) iproc++;
              // Exclude RH neutrinos from allowed final states
              if (abs(id1) >= 2000012 && id1 % 2 == 0) continue;
              if (abs(id2) >= 2000012 && id2 % 2 == 0) continue;
              // Skip if outgoing codes not asked for
              if (!allowIdVals( id1, id2)) continue;
              if (iso == jso && id1 != id2) {
                sigmaPtr
                  = new Sigma2qqbar2sleptonantislepton(id1,-id2,iproc-1);
                containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
                sigmaPtr = new Sigma2qqbar2sleptonantislepton(id2,-id1,iproc);
                containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
              } else {
                sigmaPtr = new Sigma2qqbar2sleptonantislepton(id1,-id2,iproc);
                containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
              }
            }
          }
        }
      }
    }

  } // End of SUSY processes.

  // Set up requested objects for New-Gauge-Boson processes.
  if (settings.flag("NewGaugeBoson:ffbar2gmZZprime")) {
    sigmaPtr = new Sigma1ffbar2gmZZprime();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("NewGaugeBoson:ffbar2Wprime")) {
    sigmaPtr = new Sigma1ffbar2Wprime();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("NewGaugeBoson:ffbar2R0")) {
    sigmaPtr = new Sigma1ffbar2Rhorizontal();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for Left-Right-Symmetry processes.
  bool leftrights = settings.flag("LeftRightSymmmetry:all");
  if (leftrights || settings.flag("LeftRightSymmmetry:ffbar2ZR")) {
    sigmaPtr = new Sigma1ffbar2ZRight();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:ffbar2WR")) {
    sigmaPtr = new Sigma1ffbar2WRight();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:ll2HL")) {
    sigmaPtr = new Sigma1ll2Hchgchg(1);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:lgm2HLe")) {
    sigmaPtr = new Sigma2lgm2Hchgchgl(1, 11);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:lgm2HLmu")) {
    sigmaPtr = new Sigma2lgm2Hchgchgl(1, 13);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:lgm2HLtau")) {
    sigmaPtr = new Sigma2lgm2Hchgchgl(1, 15);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:ff2HLff")) {
    sigmaPtr = new Sigma3ff2HchgchgfftWW(1);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:ffbar2HLHL")) {
    sigmaPtr = new Sigma2ffbar2HchgchgHchgchg(1);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:ll2HR")) {
    sigmaPtr = new Sigma1ll2Hchgchg(2);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:lgm2HRe")) {
    sigmaPtr = new Sigma2lgm2Hchgchgl(2, 11);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:lgm2HRmu")) {
    sigmaPtr = new Sigma2lgm2Hchgchgl(2, 13);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:lgm2HRtau")) {
    sigmaPtr = new Sigma2lgm2Hchgchgl(2, 15);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:ff2HRff")) {
    sigmaPtr = new Sigma3ff2HchgchgfftWW(2);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leftrights || settings.flag("LeftRightSymmmetry:ffbar2HRHR")) {
    sigmaPtr = new Sigma2ffbar2HchgchgHchgchg(2);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for leptoquark LQ processes.
  bool leptoquarks = settings.flag("LeptoQuark:all");
  if (leptoquarks || settings.flag("LeptoQuark:ql2LQ")) {
    sigmaPtr = new Sigma1ql2LeptoQuark;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leptoquarks || settings.flag("LeptoQuark:qg2LQl")) {
    sigmaPtr = new Sigma2qg2LeptoQuarkl;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leptoquarks || settings.flag("LeptoQuark:gg2LQLQbar")) {
    sigmaPtr = new Sigma2gg2LQLQbar;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (leptoquarks || settings.flag("LeptoQuark:qqbar2LQLQbar")) {
    sigmaPtr = new Sigma2qqbar2LQLQbar;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for excited-fermion processes.
  bool excitedfermions = settings.flag("ExcitedFermion:all");
  if (excitedfermions || settings.flag("ExcitedFermion:dg2dStar")) {
    sigmaPtr = new Sigma1qg2qStar(1);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:ug2uStar")) {
    sigmaPtr = new Sigma1qg2qStar(2);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:sg2sStar")) {
    sigmaPtr = new Sigma1qg2qStar(3);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:cg2cStar")) {
    sigmaPtr = new Sigma1qg2qStar(4);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:bg2bStar")) {
    sigmaPtr = new Sigma1qg2qStar(5);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:egm2eStar")) {
    sigmaPtr = new Sigma1lgm2lStar(11);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:mugm2muStar")) {
    sigmaPtr = new Sigma1lgm2lStar(13);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:taugm2tauStar")) {
    sigmaPtr = new Sigma1lgm2lStar(15);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qq2dStarq")) {
    sigmaPtr = new Sigma2qq2qStarq(1);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qq2uStarq")) {
    sigmaPtr = new Sigma2qq2qStarq(2);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qq2sStarq")) {
    sigmaPtr = new Sigma2qq2qStarq(3);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qq2cStarq")) {
    sigmaPtr = new Sigma2qq2qStarq(4);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qq2bStarq")) {
    sigmaPtr = new Sigma2qq2qStarq(5);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qqbar2eStare")) {
    sigmaPtr = new Sigma2qqbar2lStarlbar(11);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qqbar2nueStarnue")) {
    sigmaPtr = new Sigma2qqbar2lStarlbar(12);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qqbar2muStarmu")) {
    sigmaPtr = new Sigma2qqbar2lStarlbar(13);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qqbar2numuStarnumu")) {
    sigmaPtr = new Sigma2qqbar2lStarlbar(14);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qqbar2tauStartau")) {
    sigmaPtr = new Sigma2qqbar2lStarlbar(15);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions
    || settings.flag("ExcitedFermion:qqbar2nutauStarnutau")) {
    sigmaPtr = new Sigma2qqbar2lStarlbar(16);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qqbar2eStareStar")) {
    sigmaPtr = new Sigma2qqbar2lStarlStarBar(11);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions
    || settings.flag("ExcitedFermion:qqbar2nueStarnueStar")) {
    sigmaPtr = new Sigma2qqbar2lStarlStarBar(12);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions || settings.flag("ExcitedFermion:qqbar2muStarmuStar")) {
    sigmaPtr = new Sigma2qqbar2lStarlStarBar(13);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions
    || settings.flag("ExcitedFermion:qqbar2numuStarnumuStar")) {
    sigmaPtr = new Sigma2qqbar2lStarlStarBar(14);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions
    || settings.flag("ExcitedFermion:qqbar2tauStartauStar")) {
    sigmaPtr = new Sigma2qqbar2lStarlStarBar(15);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (excitedfermions
    || settings.flag("ExcitedFermion:qqbar2nutauStarnutauStar")) {
    sigmaPtr = new Sigma2qqbar2lStarlStarBar(16);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for contact interaction processes.
  if (settings.flag("ContactInteractions:QCqq2qq")) {
    sigmaPtr = new Sigma2QCqq2qq();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ContactInteractions:QCqqbar2qqbar")) {
    sigmaPtr = new Sigma2QCqqbar2qqbar();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ContactInteractions:QCffbar2eebar")) {
    sigmaPtr = new Sigma2QCffbar2llbar(11, 4203);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ContactInteractions:QCffbar2mumubar")) {
    sigmaPtr = new Sigma2QCffbar2llbar(13, 4204);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ContactInteractions:QCffbar2tautaubar")) {
    sigmaPtr = new Sigma2QCffbar2llbar(15, 4205);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set spin of particles in the Hidden Valley scenario.
  int spinFv = settings.mode("HiddenValley:spinFv");
  for (int i = 1; i < 7; ++i) {
    if (particleDataPtr->spinType( 4900000 + i) != spinFv + 1)
        particleDataPtr->spinType( 4900000 + i,    spinFv + 1);
    if (particleDataPtr->spinType( 4900010 + i) != spinFv + 1)
        particleDataPtr->spinType( 4900010 + i,    spinFv + 1);
  }
  if (spinFv != 1) {
    if (particleDataPtr->spinType( 4900101) != 2)
       particleDataPtr->spinType( 4900101, 2);
  } else {
    int spinqv = settings.mode("HiddenValley:spinqv");
    if (particleDataPtr->spinType( 4900101) != 2 * spinqv + 1)
        particleDataPtr->spinType( 4900101,    2 * spinqv + 1);
  }

  // Set up requested objects for HiddenValley processes.
  bool hiddenvalleys = settings.flag("HiddenValley:all");
  if (hiddenvalleys || settings.flag("HiddenValley:gg2DvDvbar")) {
    sigmaPtr = new Sigma2gg2qGqGbar( 4900001, 4901, spinFv,
      "g g -> Dv Dvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:gg2UvUvbar")) {
    sigmaPtr = new Sigma2gg2qGqGbar( 4900002, 4902, spinFv,
      "g g -> Uv Uvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:gg2SvSvbar")) {
    sigmaPtr = new Sigma2gg2qGqGbar( 4900003, 4903, spinFv,
      "g g -> Sv Svbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:gg2CvCvbar")) {
    sigmaPtr = new Sigma2gg2qGqGbar( 4900004, 4904, spinFv,
      "g g -> Cv Cvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:gg2BvBvbar")) {
    sigmaPtr = new Sigma2gg2qGqGbar( 4900005, 4905, spinFv,
      "g g -> Bv Bvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:gg2TvTvbar")) {
    sigmaPtr = new Sigma2gg2qGqGbar( 4900006, 4906, spinFv,
      "g g -> Tv Tvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:qqbar2DvDvbar")) {
    sigmaPtr = new Sigma2qqbar2qGqGbar( 4900001, 4911, spinFv,
      "q qbar -> Dv Dvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:qqbar2UvUvbar")) {
    sigmaPtr = new Sigma2qqbar2qGqGbar( 4900002, 4912, spinFv,
      "q qbar -> Uv Uvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:qqbar2SvSvbar")) {
    sigmaPtr = new Sigma2qqbar2qGqGbar( 4900003, 4913, spinFv,
      "q qbar -> Sv Svbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:qqbar2CvCvbar")) {
    sigmaPtr = new Sigma2qqbar2qGqGbar( 4900004, 4914, spinFv,
      "q qbar -> Cv Cvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:qqbar2BvBvbar")) {
    sigmaPtr = new Sigma2qqbar2qGqGbar( 4900005, 4915, spinFv,
      "q qbar -> Bv Bvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:qqbar2TvTvbar")) {
    sigmaPtr = new Sigma2qqbar2qGqGbar( 4900006, 4916, spinFv,
      "q qbar -> Tv Tvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2DvDvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900001, 4921, spinFv,
      "f fbar -> Dv Dvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2UvUvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900002, 4922, spinFv,
      "f fbar -> Uv Uvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2SvSvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900003, 4923, spinFv,
      "f fbar -> Sv Svbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2CvCvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900004, 4924, spinFv,
      "f fbar -> Cv Cvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2BvBvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900005, 4925, spinFv,
      "f fbar -> Bv Bvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2TvTvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900006, 4926, spinFv,
      "f fbar -> Tv Tvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2EvEvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900011, 4931, spinFv,
      "f fbar -> Ev Evbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2nuEvnuEvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900012, 4932, spinFv,
      "f fbar -> nuEv nuEvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2MUvMUvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900013, 4933, spinFv,
      "f fbar -> MUv MUvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2nuMUvnuMUvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900014, 4934, spinFv,
      "f fbar -> nuMUv nuMUvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2TAUvTAUvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900015, 4935, spinFv,
      "f fbar -> TAUv TAUvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2nuTAUvnuTAUvbar")) {
    sigmaPtr = new Sigma2ffbar2fGfGbar( 4900016, 4936, spinFv,
      "f fbar -> nuTAUv nuTAUvbar");
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (hiddenvalleys || settings.flag("HiddenValley:ffbar2Zv")) {
    sigmaPtr = new Sigma1ffbar2Zv();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for RS extra-dimensional G* processes.
  bool extraDimGstars = settings.flag("ExtraDimensionsG*:all");
  if (extraDimGstars || settings.flag("ExtraDimensionsG*:gg2G*")) {
    sigmaPtr = new Sigma1gg2GravitonStar;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimGstars || settings.flag("ExtraDimensionsG*:ffbar2G*")) {
    sigmaPtr = new Sigma1ffbar2GravitonStar;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsG*:gg2G*g")) {
    sigmaPtr = new Sigma2gg2GravitonStarg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsG*:qg2G*q")) {
    sigmaPtr = new Sigma2qg2GravitonStarq;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsG*:qqbar2G*g")) {
    sigmaPtr = new Sigma2qqbar2GravitonStarg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  //  Set up requested objects for RS extra-dimensional KKgluon processes.
  if (settings.flag("ExtraDimensionsG*:qqbar2KKgluon*")) {
    sigmaPtr = new Sigma1qqbar2KKgluonStar;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // NOAM: Set up requested objects for TEV extra-dimensional processes.
  if (settings.flag("ExtraDimensionsTEV:ffbar2ddbar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(1, 5061);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2uubar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(2, 5062);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2ssbar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(3, 5063);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2ccbar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(4, 5064);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2bbbar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(5, 5065);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2ttbar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(6, 5066);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2e+e-")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(11, 5071);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2nuenuebar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(12, 5072);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2mu+mu-")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(13, 5073);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2numunumubar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(14, 5074);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2tau+tau-")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(15, 5075);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsTEV:ffbar2nutaunutaubar")) {
    sigmaPtr = new Sigma2ffbar2TEVffbar(16, 5076);
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for large extra-dimensional G processes.
  bool extraDimLEDmono = settings.flag("ExtraDimensionsLED:monojet");
  if (extraDimLEDmono || settings.flag("ExtraDimensionsLED:gg2Gg")) {
    sigmaPtr = new Sigma2gg2LEDUnparticleg( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimLEDmono || settings.flag("ExtraDimensionsLED:qg2Gq")) {
    sigmaPtr = new Sigma2qg2LEDUnparticleq( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimLEDmono || settings.flag("ExtraDimensionsLED:qqbar2Gg")) {
    sigmaPtr = new Sigma2qqbar2LEDUnparticleg( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsLED:ffbar2GZ")) {
    sigmaPtr = new Sigma2ffbar2LEDUnparticleZ( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsLED:ffbar2Ggamma")) {
    sigmaPtr = new Sigma2ffbar2LEDUnparticlegamma( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsLED:ffbar2gammagamma")) {
    sigmaPtr = new Sigma2ffbar2LEDgammagamma( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsLED:gg2gammagamma")) {
    sigmaPtr = new Sigma2gg2LEDgammagamma( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsLED:ffbar2llbar")) {
    sigmaPtr = new Sigma2ffbar2LEDllbar( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsLED:gg2llbar")) {
    sigmaPtr = new Sigma2gg2LEDllbar( true );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  bool extraDimLEDdij = settings.flag("ExtraDimensionsLED:dijets");
  if (extraDimLEDdij || settings.flag("ExtraDimensionsLED:gg2DJgg")) {
    sigmaPtr = new Sigma2gg2LEDgg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimLEDdij || settings.flag("ExtraDimensionsLED:gg2DJqqbar")) {
    sigmaPtr = new Sigma2gg2LEDqqbar;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimLEDdij || settings.flag("ExtraDimensionsLED:qg2DJqg")) {
    sigmaPtr = new Sigma2qg2LEDqg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimLEDdij || settings.flag("ExtraDimensionsLED:qq2DJqq")) {
    sigmaPtr = new Sigma2qq2LEDqq;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimLEDdij || settings.flag("ExtraDimensionsLED:qqbar2DJgg")) {
    sigmaPtr = new Sigma2qqbar2LEDgg;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimLEDdij || settings.flag("ExtraDimensionsLED:qqbar2DJqqbarNew")) {
    sigmaPtr = new Sigma2qqbar2LEDqqbarNew;
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for unparticle processes.
  bool extraDimUnpartmono = settings.flag("ExtraDimensionsUnpart:monojet");
  if (extraDimUnpartmono || settings.flag("ExtraDimensionsUnpart:gg2Ug")) {
    sigmaPtr = new Sigma2gg2LEDUnparticleg( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimUnpartmono || settings.flag("ExtraDimensionsUnpart:qg2Uq")) {
    sigmaPtr = new Sigma2qg2LEDUnparticleq( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (extraDimUnpartmono || settings.flag("ExtraDimensionsUnpart:qqbar2Ug")) {
    sigmaPtr = new Sigma2qqbar2LEDUnparticleg( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsUnpart:ffbar2UZ")) {
    sigmaPtr = new Sigma2ffbar2LEDUnparticleZ( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsUnpart:ffbar2Ugamma")) {
    sigmaPtr = new Sigma2ffbar2LEDUnparticlegamma( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsUnpart:ffbar2gammagamma")) {
    sigmaPtr = new Sigma2ffbar2LEDgammagamma( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsUnpart:gg2gammagamma")) {
    sigmaPtr = new Sigma2gg2LEDgammagamma( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsUnpart:ffbar2llbar")) {
    sigmaPtr = new Sigma2ffbar2LEDllbar( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("ExtraDimensionsUnpart:gg2llbar")) {
    sigmaPtr = new Sigma2gg2LEDllbar( false );
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Set up requested objects for Dark Matter processes.
  if (settings.flag("DM:ffbar2Zp2XX")) {
    sigmaPtr = new Sigma1ffbar2Zp2XX();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("DM:ffbar2Zp2XXj")) {
    sigmaPtr = new Sigma2qqbar2Zpg2XXj();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("DM:ffbar2ZpH")) {
    sigmaPtr = new Sigma2ffbar2ZpH();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("DM:gg2S2XX")) {
    sigmaPtr = new Sigma1gg2S2XX();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }
  if (settings.flag("DM:gg2S2XXj")) {
    sigmaPtr = new Sigma2gg2Sg2XXj();
    containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Done.
  return true;

}

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

// Routine to initialize list of second hard processes.

bool SetupContainers::init2(vector<ProcessContainer*>& container2Ptrs,
  Settings& settings) {

  // Reset process list, if filled in previous subrun.
  if (container2Ptrs.size() > 0) {
    for (int i = 0; i < int(container2Ptrs.size()); ++i)
      delete container2Ptrs[i];
    container2Ptrs.clear();
  }
  SigmaProcess* sigmaPtr;

  // Two hard QCD jets.
  if (settings.flag("SecondHard:TwoJets")) {
    sigmaPtr = new Sigma2gg2gg;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2gg2qqbar;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qg2qg;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qq2qq;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qqbar2gg;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qqbar2qqbarNew;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2gg2QQbar(4, 121);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qqbar2QQbar(4, 122);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2gg2QQbar(5, 123);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qqbar2QQbar(5, 124);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // A prompt photon and a hard jet.
  if (settings.flag("SecondHard:PhotonAndJet")) {
    sigmaPtr = new Sigma2qg2qgamma;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qqbar2ggamma;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2gg2ggamma;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Two prompt photons.
  if (settings.flag("SecondHard:TwoPhotons")) {
    sigmaPtr = new Sigma2ffbar2gammagamma;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2gg2gammagamma;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Charmonium.
  if (settings.flag("SecondHard:Charmonium")) {
    vector<SigmaProcess*> charmoniumSigmaPtrs;
    charmonium.setupSigma2gg(charmoniumSigmaPtrs, true);
    charmonium.setupSigma2qg(charmoniumSigmaPtrs, true);
    charmonium.setupSigma2qq(charmoniumSigmaPtrs, true);
    for (unsigned int i = 0; i < charmoniumSigmaPtrs.size(); ++i)
      container2Ptrs.push_back( new ProcessContainer(charmoniumSigmaPtrs[i]));
  }

  // Bottomonium.
  if (settings.flag("SecondHard:Bottomonium")) {
    vector<SigmaProcess*> bottomoniumSigmaPtrs;
    bottomonium.setupSigma2gg(bottomoniumSigmaPtrs, true);
    bottomonium.setupSigma2qg(bottomoniumSigmaPtrs, true);
    bottomonium.setupSigma2qq(bottomoniumSigmaPtrs, true);
    for (unsigned int i = 0; i < bottomoniumSigmaPtrs.size(); ++i)
      container2Ptrs.push_back( new ProcessContainer(bottomoniumSigmaPtrs[i]));
  }

  // A single gamma*/Z0.
  if (settings.flag("SecondHard:SingleGmZ")) {
    sigmaPtr = new Sigma1ffbar2gmZ;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // A single W+-.
  if (settings.flag("SecondHard:SingleW")) {
    sigmaPtr = new Sigma1ffbar2W;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // A gamma*/Z0 and a hard jet.
  if (settings.flag("SecondHard:GmZAndJet")) {
    sigmaPtr = new Sigma2qqbar2gmZg;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qg2gmZq;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // A W+- and a hard jet.
  if (settings.flag("SecondHard:WAndJet")) {
    sigmaPtr = new Sigma2qqbar2Wg;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qg2Wq;
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Top pair production.
  if (settings.flag("SecondHard:TopPair")) {
    sigmaPtr = new Sigma2gg2QQbar(6, 601);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qqbar2QQbar(6, 602);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2ffbar2FFbarsgmZ(6, 604);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Single top production.
  if (settings.flag("SecondHard:SingleTop")) {
    sigmaPtr = new Sigma2qq2QqtW(6, 603);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2ffbar2FfbarsW(6, 0, 605);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Two b jets - already part of TwoJets sample above.
  if (settings.flag("SecondHard:TwoBJets")) {
    sigmaPtr = new Sigma2gg2QQbar(5, 123);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
    sigmaPtr = new Sigma2qqbar2QQbar(5, 124);
    container2Ptrs.push_back( new ProcessContainer(sigmaPtr) );
  }

  // Done.
  return true;

}

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

// Set up arrays of allowed outgoing SUSY particles.

void SetupContainers::setupIdVecs( Settings& settings) {

  // First array either none, one or many particles.
  idVecA.resize(0);
  if (settings.mode("SUSY:idA") != 0) {
    idVecA.push_back( abs(settings.mode("SUSY:idA")) );
  } else {
    vector<int> idTmpA = settings.mvec("SUSY:idVecA");
    for (int i = 0; i < int(idTmpA.size()); ++i)
      if (idTmpA[i] != 0) idVecA.push_back( abs(idTmpA[i]) );
  }
  nVecA = idVecA.size();

  // Second array either none, one or many particles.
  idVecB.resize(0);
  if (settings.mode("SUSY:idB") != 0) {
    idVecB.push_back( abs(settings.mode("SUSY:idB")) );
  } else {
    vector<int> idTmpB = settings.mvec("SUSY:idVecB");
    for (int i = 0; i < int(idTmpB.size()); ++i)
    if (idTmpB[i] != 0) idVecB.push_back( abs(idTmpB[i]) );
  }
  nVecB = idVecB.size();

}

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

// Check final state for allowed outgoing SUSY particles.
// Normally check two codes, but allow for only one.

bool SetupContainers::allowIdVals( int idCheck1, int idCheck2) {

  // If empty arrays or id's no need for checks. Else need absolute values.
  if (nVecA == 0 && nVecB == 0) return true;
  if (idCheck1 == 0 && idCheck2 == 0) return true;
  int idChk1 = abs(idCheck1);
  int idChk2 = abs(idCheck2);

  // If only one outgoing particle then check idVecA and idVecB.
  if (idChk1 == 0) swap(idChk1, idChk2);
  if (idChk2 == 0) {
    for (int i = 0; i < nVecA; ++i) if (idChk1 == idVecA[i]) return true;
    for (int i = 0; i < nVecB; ++i) if (idChk1 == idVecB[i]) return true;
    return false;
  }

  // If empty array idVecB then compare with idVecA.
  if (nVecB == 0) {
    for (int i = 0; i < nVecA; ++i)
      if (idChk1 == idVecA[i] || idChk2 == idVecA[i]) return true;
    return false;
  }

  // If empty array idVecA then compare with idVecB.
  if (nVecA == 0) {
    for (int i = 0; i < nVecB; ++i)
      if (idChk1 == idVecB[i] || idChk2 == idVecB[i]) return true;
    return false;
  }

  // Else check that pair matches allowed combinations.
  for (int i = 0; i < nVecA; ++i)
  for (int j = 0; j < nVecB; ++j)
    if ( (idChk1 == idVecA[i] && idChk2 == idVecB[j])
      || (idChk2 == idVecA[i] && idChk1 == idVecB[j]) ) return true;
  return false;

}

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

} // end namespace Pythia8