VinciaCommon.h

// VinciaCommon.h is a part of the PYTHIA event generator.
// Copyright (C) 2024 Peter Skands, Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This file contains global constants for the Vincia antenna shower as
// well as various auxiliary classes used by the Vincia shower model.

#ifndef Pythia8_VinciaCommon_H
#define Pythia8_VinciaCommon_H

// Maths headers.
#include <limits>

// Include Pythia 8 headers.
#include "Pythia8/Event.h"
#include "Pythia8/Info.h"
#include "Pythia8/ParticleData.h"
#include "Pythia8/PartonSystems.h"
#include "Pythia8/PythiaStdlib.h"
#include "Pythia8/StandardModel.h"

// Global Vincia constants, defined to live in a Vincia-specific
// namespace to avoid potential clashes with anything else defined in Pythia8.
namespace VinciaConstants {

// Global numerical precision targets (should be kept within a
// reasonable margin of machine precision).  Large values may result
// in strange results due to non-zero terms being artificially forced
// to zero.
const double UNITY = 1.0;
const double DECI  = 1.0e-1;
const double CENTI = 1.0e-2;
const double MILLI = 1.0e-3;
const double MICRO = 1.0e-6;
const double NANO  = 1.0e-9;
const double PICO  = 1.0e-12;

// Colour factors in Vincia normalization.
const double CA = 3.0;
const double CF = 8.0/3.0;
const double TR = 1.0;
const double NC = 3.0;

// Mathematical constants (Euler–Mascheroni constant).
const double gammaE = 0.577215664901532860606512090082402431042;

// Verbosity levels. Vincia has one more level (debug) beyond report.
const int DEBUG   = 4;

// Padding length for dashes in standardised Vincia verbose output.
const int DASHLEN = 80;

}

// Everything else lives in the Pythia 8 namespace.

namespace Pythia8 {

// Forward declaration.
class VinciaCommon;

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

// Enumerator for antenna function types, with "void" member NoFun.
enum AntFunType { NoFun,
                  QQEmitFF, QGEmitFF, GQEmitFF, GGEmitFF, GXSplitFF,
                  QQEmitRF, QGEmitRF, XGSplitRF,
                  QQEmitII, GQEmitII, GGEmitII, QXConvII, GXConvII,
                  QQEmitIF, QGEmitIF, GQEmitIF, GGEmitIF, QXConvIF,
                  GXConvIF, XGSplitIF };

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

// Convenient typedef for unsigned integers.

typedef unsigned int uint;

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

// Print a method name using the appropritae pre-processor macro.

//  The following method was modified from
//  boost/current_function.hpp - BOOST_CURRENT_FUNCTION
//
//  Copyright (c) 2002 Peter Dimov and Multi Media Ltd.
//
//  Distributed under the Boost Software License, Version 1.0.
//  Boost Software License - Version 1.0 - August 17th, 2003
//
//  Permission is hereby granted, free of charge, to any person or
//  organization obtaining a copy of the software and accompanying
//  documentation covered by this license (the "Software") to use,
//  reproduce, display, distribute, execute, and transmit the
//  Software, and to prepare derivative works of the Software, and to
//  permit third-parties to whom the Software is furnished to do so,
//  all subject to the following:
//
//  The copyright notices in the Software and this entire statement,
//  including the above license grant, this restriction and the
//  following disclaimer, must be included in all copies of the
//  Software, in whole or in part, and all derivative works of the
//  Software, unless such copies or derivative works are solely in the
//  form of machine-executable object code generated by a source
//  language processor.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
//  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
//  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
//  NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
//  ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
//  OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
//  OTHER DEALINGS IN THE SOFTWARE.
//
//  http://www.boost.org/libs/utility/current_function.html
//
//  Note that Boost Software License - Version 1.0 is fully compatible
//  with GPLV2
//  For more information see https://www.gnu.org/licenses/license-list.en.html

#ifndef __METHOD_NAME__

#ifndef VINCIA_FUNCTION
#if ( defined(__GNUC__) || (defined(__MWERKS__) && (__MWERKS__ >= 0x3000)) \
|| (defined(__ICC) && (__ICC >= 600)) )
# define VINCIA_FUNCTION __PRETTY_FUNCTION__
#elif defined(__DMC__) && (__DMC__ >= 0x810)
# define VINCIA_FUNCTION __PRETTY_FUNCTION__
#elif defined(__FUNCSIG__)
# define VINCIA_FUNCTION __FUNCSIG__
#elif ( (defined(__INTEL_COMPILER) && (__INTEL_COMPILER >= 600)) \
|| (defined(__IBMCPP__) && (__IBMCPP__ >= 500)) )
# define VINCIA_FUNCTION __FUNCTION__
#elif defined(__BORLANDC__) && (__BORLANDC__ >= 0x550)
# define VINCIA_FUNCTION __FUNC__
#elif defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901)
# define VINCIA_FUNCTION __func__
#else
# define VINCIA_FUNCTION "unknown"
#endif
#endif // end VINCIA_FUNCTION

inline std::string methodName(const std::string& prettyFunction, bool
  withPythia=false) {
  size_t colons = prettyFunction.find("::");
  // Include Pythia8:: or not.
  size_t begin = colons + 2;
  if (withPythia) begin = prettyFunction.substr(0,colons).rfind(" ") + 1;
  size_t end = prettyFunction.rfind("(") - begin;
  return prettyFunction.substr(begin,end) + "()";
}

#define __METHOD_NAME__ methodName(VINCIA_FUNCTION)
#endif // end __METHOD_NAME__

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

// Global functions accessible in the Vincia namespace.

// Utilities for printing out VINCIA Messages.
void printOut(string,string,int nPad=0,char padChar='-');

// String utilities.
string num2str(int,int width=4) ;
string num2str(double,int width=9) ;
string bool2str(bool, int width=3) ;

// Search and replace a string.
inline string replaceString(string subject, const string& search,
  const string& replace) {
  string::size_type pos = 0;
  while ((pos = subject.find(search, pos)) != string::npos) {
    subject.replace(pos, search.length(), replace);
    pos += replace.length();
  }
  return subject;

}

// Remove ":" and "/" from a file name.
inline string sanitizeFileName(string fileName) {
  map<string, string> rep;
  rep["/"] = "_div_";
  rep[":"] = "_colon_";
  string retVal = fileName;
  for (map<string, string>::const_iterator it = rep.begin(); it != rep.end();
       ++it) {
    retVal = replaceString(retVal, it->first, it->second);
  }
  return retVal;

}

// Utility for checking if a file exists.
inline bool fileExists(const std::string& name) {
  if (FILE *file = fopen(name.c_str(), "r")) {
    fclose(file);
    return true;
  } else {
    return false;
  }

}

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

// A class to store and process colour information, e.g. colour maps,
// reconnections, etc.

class VinciaColour {

public:

  // Initialize pointers (must be done before init).
  void initPtr(Info* infoPtrIn) {
    infoPtr       = infoPtrIn;
    particleDataPtr  = infoPtr->particleDataPtr;
    settingsPtr      = infoPtr->settingsPtr;
    partonSystemsPtr = infoPtr->partonSystemsPtr;
    rndmPtr          = infoPtr->rndmPtr;
    isInitPtr=true;
  }

  // Initialize.
  bool init();

  // Get translation map from ordinary Les Houches tags to antenna
  // indices (tags always ending on [1-9] with gluons having
  // non-identical colour and anticolour indices).
  bool colourise(int iSys, Event& event);

  // Method to sort list of partons in Vincia colour order. Returns
  // vector<int> with indices in following order:
  //   (1) colourless incoming particles
  //   (2) triplet-ngluon-antitriplet contractions
  //   (3) ngluon rings
  //   (4) colourless outgoing particles
  // No specific order is imposed on the relative ordering inside each
  // of the four classes.
  vector<int> colourSort(vector<Particle*>);

  // Method to create LC colour map and list of LC antennae.
  void makeColourMaps(const int iSysIn, const Event& event,
    map<int,int>& indexOfAcol, map<int,int>& indexOfCol,
    vector< pair<int,int> >& antLC, const bool findFF, const bool findIX);

  // Determine whether 01 or 12 inherit the colour tag from a parent.
  bool inherit01(double s01, double s12);

  // Set verbose level.
  void setVerbose(int verboseIn) {verbose = verboseIn;};

private:

  // Internal parameters.
  int inheritMode{};

  // Is initialized.
  bool isInitPtr{false}, isInit{false};

  // Pointers to PYTHIA 8 objects.
  Info*          infoPtr;
  ParticleData*  particleDataPtr;
  Settings*      settingsPtr;
  PartonSystems* partonSystemsPtr;
  Rndm*          rndmPtr;

  // Verbose level.
  int verbose{};

};

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

// Simple struct to store information about a 3 -> 2 clustering.

struct VinciaClustering {

  // Set information about daughters from current event and daughters indices.
  void setDaughters(const Event& state, int dau1In, int dau2In, int dau3In);
  void setDaughters(const vector<Particle>& state, int dau1In, int dau2In,
    int dau3In);

  // Set mother particle ids.
  void setMothers(int idMot1In, int idMot2In) {
    idMot1 = idMot1In;
    idMot2 = idMot2In;
  }

  // Set antenna information.
  void setAntenna(bool isFSRin, enum AntFunType antFunTypeIn) {
    isFSR = isFSRin;
    antFunType = antFunTypeIn;
  }

  // Initialise vectors of invariants and masses.
  bool initInvariantAndMassVecs();

  // Set invariants and masses.
  void setInvariantsAndMasses(const Event& state);
  void setInvariantsAndMasses(const vector<Particle>& state);

  // Swap 1 <-> 3, including information about parents.
  void swap13() {
    swap(dau1,dau3);
    swap(idMot1,idMot2);
    swap(saj,sjb);
    if (mDau.size() == 3)
      swap(mDau[0],mDau[2]);
    if (mMot.size() == 2)
      swap(mMot[0],mMot[1]);
    if (invariants.size() == 3) {
      swap(invariants[1],invariants[2]);
    }
  }

  // Methods to get antenna type.
  bool isFF() const {
    if (!isFSR) return false;
    else if (antFunType >= QQEmitFF && antFunType < QQEmitRF) return true;
    else return false;
  }
  bool isRF() const {
    if (!isFSR) return false;
    else if (antFunType >= QQEmitRF && antFunType < QQEmitII) return true;
    else return false;
  }
  bool isII() const {
    if (isFSR) return false;
    else if (antFunType >= QQEmitII && antFunType < QQEmitIF) return true;
    return false;
  }
  bool isIF() const {
    if (isFSR) return false;
    else if (antFunType >= QQEmitIF) return true;
    else return false;
  }
  string getAntName() const;

  // Methods to get branching type (currently only 2 -> 3).
  bool is2to3() const { return true; }

  // Branching daughter information (indices in event record).
  int dau1{}, dau2{}, dau3{};

  // Antenna information.
  bool isFSR{true};
  AntFunType antFunType{NoFun};

  // Mother ids.
  int idMot1{}, idMot2{};

  // Helicities.
  vector<int> helDau = {9, 9, 9};
  vector<int> helMot = {9, 9};

  // Masses.
  vector<double> mDau;
  vector<double> mMot;

  // Invariants.
  double saj{}, sjb{}, sab{};
  // Vector of invariants (stored as sAB, saj, sjb, sab).
  vector<double> invariants;

  // Value of sector resolution variable that this clustering produces.
  double q2res{};

  // Value of evolution variable that this clustering produces.
  double q2evol{};

  // Kinematic map (only used for FF).
  int kMapType{};

};

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

// A simple class for containing evolution variable definitions.

class Resolution {

public:

  // Initialize pointers (must be done before init).
  void initPtr(Settings* settingsPtrIn, Info* infoPtrIn,
    VinciaCommon* vinComPtrIn) {
    settingsPtr = settingsPtrIn;
    infoPtr     = infoPtrIn;
    loggerPtr   = infoPtrIn->loggerPtr;
    vinComPtr   = vinComPtrIn;
    isInitPtr   = true;
  }

  // Initialize.
  bool init();

  // Method to calculate (and set) evolution variable.
  double q2evol(VinciaClustering& clus);

  // Method to calculate dimensionless evolution variable.
  // Note: calls q2evol(), so will set dimensionful evolution
  // variable in the VinciaClustering object.
  double xTevol(VinciaClustering& clus);

  // Top-level function to calculate (and set) sector resolution.
  double q2sector(VinciaClustering& clus);

  // Find sector with minimal resolution.
  // Optionally resolve Born: avoid clusterings that would not lead
  // to a specified (Born) configuration.
  VinciaClustering findSector(vector<Particle>& state,
    map<int,int> flavsBorn);
  // Find sector with minimal resolution.
  // Optionally avoid sectors that cluster beyond a minimal number
  // of quark pairs or gluons.
  VinciaClustering findSector(vector<Particle>& state,
    int nqpMin = 0, int ngMin = 0);

  // Find sector with minimal q2sector in list of clusterings.
  VinciaClustering getMinSector(vector<VinciaClustering>& clusterings);

  // Sector veto to check whether given value of resolution is minimal,
  // given we want to preserve a certain Born configuration.
  // Returns true = vetoed, and false = not vetoed.
  bool sectorVeto(double q2In, vector<Particle>& state,
    map<int,int> nFlavsBorn) {
    VinciaClustering clusMin = findSector(state, nFlavsBorn);
    if (q2In > clusMin.q2res) return true;
    else return false;
  }
  bool sectorVeto(const VinciaClustering& clusMin,
    const VinciaClustering& clus);

  // Set verbosity level.
  void setVerbose(int verboseIn) {verbose = verboseIn;}

private:

  // Member functions.

  // Sector resolution for 2 -> 3 branchings.
  double q2sector2to3FF(VinciaClustering& clus);
  double q2sector2to3RF(VinciaClustering& clus);
  double q2sector2to3IF(VinciaClustering& clus);
  double q2sector2to3II(VinciaClustering& clus);

  // Sector resolution function for 3->4 branchings (currently only
  // used for gluon splitting, with m2qq as the measure).
  //TODO: currently disabled.
  //   double q2sector3to4(const Particle*, const Particle*,
  //     const Particle* j1, const Particle* j2) {return -1.;}

  // Sector resolution function for 2->4 branchings (double emission).
  // Assume j1 and j2 are colour connected, with a and b hard
  // recoilers.
  //TODO: currently disabled.
  //   double q2sector2to4(const Particle* a, const Particle* b,
  //     const Particle* j1, const Particle* j2) {return -1.;}

  // Sector resolution function for 3->5 branchings (emission +
  // splitting).
  //TODO: currently disabled.
  //   double q2sector3to5(Particle* a, Particle* b,
  //     Particle* j1, Particle* j2, Particle* j3) {return -1;}

  // Members.

  // Initialized.
  bool isInitPtr{false}, isInit{false};

  // Pointer to PYTHIA 8 settings database.
  Settings* settingsPtr{};
  Info* infoPtr{};
  Logger* loggerPtr{};

  // Pointer to VinciaCommon.
  VinciaCommon* vinComPtr{};

  // Number of flavours to be treated as massless.
  int nFlavZeroMassSav{};

  // Verbosity level.
  int verbose{};

};

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

// Class which contains functions and variables shared by the
// VinciaShower and VinciaMatching classes.

class VinciaCommon {

public:

  // Initialize pointers.
  bool initPtr(Info* infoPtrIn) {
    infoPtr          = infoPtrIn;
    particleDataPtr  = infoPtr->particleDataPtr;
    settingsPtr      = infoPtr->settingsPtr;
    loggerPtr        = infoPtr->loggerPtr;
    rndmPtr          = infoPtr->rndmPtr;
    partonSystemsPtr = infoPtr->partonSystemsPtr;
    isInitPtr        = true;
    return true;
  }

  // Initialize data members.
  bool init();

  // Function to check for the hadronization cutoff for a colour
  // connected parton pair.
  double mHadMin(const int id1, const int id2);

  // Function to check the event after each branching. Added by NF to
  // see if certain Pythia warnings/error are caused by the shower.
  bool showerChecks(Event& event, bool ISR);

  // Function to reset counters (print once every event for now).
  void resetCounters() {
    nUnkownPDG    = 0;
    nIncorrectCol = 0;
    nNAN          = 0;
    nVertex       = 0;
    nChargeCons   = 0;
    nMotDau       = 0;
    for (int i=0; i<2; i++) {
      nUnmatchedMass[i] = 0;
      nEPcons[i]        = 0;
    }
  }

  // Get number of active flavors at a given Q scale.
  int getNf(double q) {
    if (q <= mc) return 3;
    else if (q <= mb) return 4;
    else if (q <= mt) return 5;
    else return 6;
  }

  // Get the shower starting scale.
  double getShowerStartingScale(int iSys, const Event& event,
    double sbbSav);

  // Method to find all possible clusterings for a given system,
  // given we want to resolve a certain Born configuration, i.e.,
  // not cluster more gluons or quark flavours as we had in the Born.
  vector<VinciaClustering> findClusterings(const vector<Particle>& state,
    map<int, int> nFlavsBorn);
  // Method to find all possible clusterings while retaining a certain
  // minimal number of quark pairs and gluons.
  vector<VinciaClustering> findClusterings(const vector<Particle>& state,
    int nqpMin = 0, int ngMin = 0);

  // Check if clustering is sensible, i.e., corresponds to an existing antenna.
  bool isValidClustering(const VinciaClustering& clus,
    const Event& event, int verboseIn);

  // Perform a clustering.
  bool clus3to2(const VinciaClustering& clus, const Event& event,
    vector<Particle>& pClustered);
  bool clus3to2(const VinciaClustering& clus, const vector<Particle>& state,
    vector<Particle>& pClustered);
  // Helper functions to perform clustering.
  bool getCols3to2(const Particle* a, const Particle* j, const Particle* b,
    const VinciaClustering& clus, pair<int,int>& colsA, pair<int,int>& colsB);
  bool getMomenta3to2(vector<Vec4>& momNow, vector<Vec4>& momClus,
    const VinciaClustering& clus, int iOffset = 0);

  // 3->2 clustering maps.
  bool map3to2FF(vector<Vec4>& pClu, const vector<Vec4> pIn,
    int kMapType, int a=0, int r=1, int b=2, double mI=0.0, double mK=0.0) {
    if (mI == 0. && mK == 0.)
      return map3to2FFmassless(pClu, pIn, kMapType, a, r, b);
    else
      return map3to2FFmassive(pClu, pIn, kMapType, mI, mK, a, r, b);
  }
  bool map3to2RF(vector<Vec4>& pClu, const vector<Vec4>& pIn, int a=0,
    int r=1, int b=2, double mK = 0.);
  bool map3to2IF(vector<Vec4>& pClu, const vector<Vec4>& pIn,
    int a = 0, int r = 1, int b = 2,
    double mj = 0., double mk = 0., double mK = 0.);
  bool map3to2II(vector<Vec4>& pClu, const vector<Vec4>& pIn, bool doBoost,
    int a = 0, int r = 2, int b = 1, double mj = 0.);

  // 2->3 kinematics maps for FF branchings. Original implementations;
  // massless by Skands, massive by Ritzmann.
  bool map2to3FF(vector<Vec4>& pNew, const vector<Vec4>& pOld, int kMapType,
    const vector<double>& invariants, double phi, vector<double> masses) {
    if ( masses.size() <= 2 || ( masses[0] == 0.0 && masses[1] == 0.0
        && masses[2] == 0.0 )) {
      return map2to3FFmassless(pNew, pOld, kMapType, invariants, phi);
    } else {
      return map2to3FFmassive(pNew, pOld, kMapType, invariants, phi, masses);
    }
  }

  // 2->3 kinematics maps for II branchings. Original implementations:
  // massless by Fischer, massive by Verheyen.
  bool map2to3II(vector<Vec4>& pNew, vector<Vec4>& pRec,
    vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
    double phi, double m2j = 0.0) {
    if (m2j == 0.0)
      return map2to3IImassless(pNew, pRec, pOld, sAB, saj, sjb, sab, phi);
    else
      return map2to3IImassive(pNew, pRec, pOld, sAB, saj, sjb, sab, phi, m2j);
  }

  // 2->3 kinematics maps for IF branchings. General massive case
  // implemented by Verheyen.
  bool map2to3IFlocal(vector<Vec4>& pNew, const vector<Vec4>& pOld,
    double sOldAK, double saj, double sjk, double sak, double phi,
    double m2oldK, double m2j, double m2k);
  bool map2to3IFglobal(vector<Vec4>& pNew, vector<Vec4>& pRec,
    const vector<Vec4>& pOld, const Vec4 &pB,
    double sAK, double saj, double sjk, double sak, double phi,
    double mK2, double mj2, double mk2);

  // Resonance decay kinematic maps.
  bool map2toNRF(vector<Vec4>& pAfter, const vector<Vec4> pBefore,
    unsigned int posR, unsigned int posF,
    const vector<double> invariants, double phi,
    const vector<double> masses);

  // 1->2 decay map for (already offshell) resonance decay
  bool map1to2RF(vector<Vec4>& pNew, const Vec4 pRes, double m1,
    double m2, double theta, double phi);

  // Check if 2-particle system is on-shell and rescale if not.
  bool onShellCM(Vec4& p1, Vec4& p2, double m1, double m2, double tol = 1e-6);

  // Force initial-state and light-flavour partons to be massless.
  bool mapToMassless(int iSys, Event& event, bool makeNewCopies);

  // Map a massless antenna to equivalent massive one. Boolean
  // returns true if a modification was made.
  bool mapToMassive(Vec4& p1, Vec4& p2, double m1, double m2) {
    return (!onShellCM(p1,p2,m1,m2,1e-9));
  }

  // Make list of particles as vector<Particle>.
  //   First 1 or 2 entries : incoming particle(s).
  //   Subseqent entries    : outgoing particles.
  // The two last arguments are optional and allow to specify a list
  // of indices to be ignored, and a set of particles to be added, e.g.
  // in the context of setting up a trial state after a branching.
  // The newly added particles are then at the end of the respective
  // lists, i.e. a newly added incoming particle is the last incoming
  // one and newly added outgoing ones are the last among the outgoing
  // ones.
  vector<Particle> makeParticleList(const int iSys, const Event& event,
    const vector<Particle> &pNew = vector<Particle>(),
    const vector<int> &iOld = vector<int>());

  // Method to find all antennae that can produce a branching.
  //   IN: indices of clustering in event, where i2 is the emission.
  //  OUT: vector of VinciaClusterings.
  // Also swap daughters to match antenna function convention if needed
  // (e.g. for GXSplitFF, when dau2 and dau3 form the quark pair).
  vector<VinciaClustering> findAntennae(Event& state, int i1, int i2, int i3);

  // Check whether two particles are colour connected.
  bool colourConnected(const Particle& ptcl1, const Particle& ptcl2);

  // Print a list of Particles.
  void list(const vector<Particle>& state, string title = "",
    bool footer = true);

  // Print a list of VinciaClusterings.
  void list(const vector<VinciaClustering>& clusterings, string title = "",
    bool footer = true);

  // Get/set verbose parameter.
  int getVerbose() {return verbose; };
  void setVerbose(int verboseIn) { verbose = verboseIn;};

  // Public data members: strong coupling in MSbar and CMW schemes,
  // user and default choices,
  AlphaStrong alphaStrong{}, alphaStrongCMW{}, alphaStrongDef{},
    alphaStrongDefCMW{};

  // Couplings for use in merging.
  AlphaStrong alphaS{};
  AlphaEM     alphaEM{};
  double mu2freeze{}, mu2min{}, alphaSmax{};

  // Quark masses.
  double ms{}, mc{}, mb{}, mt{};
  int nFlavZeroMass{};

  // Checks.
  double epTolErr{}, epTolWarn{}, mTolErr{}, mTolWarn{};

private:

  // Functions.

  // Special cases of 3 -> 2 maps.
  bool map3to2FFmassive(vector<Vec4>& pClu, const vector<Vec4> pIn,
    int kMapType, double mI, double mK, int a=0, int r=1, int b=2);
  bool map3to2FFmassless(vector<Vec4>& pClu, const vector<Vec4> pIn,
    int kMapType, int a=0, int r=1, int b=2);

  // Special cases of 2 -> 3 maps.
  bool map2to3FFmassive(vector<Vec4>& pNew, const vector<Vec4>& pOld,
    int kMapType, const vector<double>& invariants, double phi,
    vector<double> masses);
  bool map2to3FFmassless(vector<Vec4>& pNew, const vector<Vec4>& pOld,
    int kMapType, const vector<double>& invariants, double phi);
  bool map2to3IImassive(vector<Vec4>& pNew, vector<Vec4>& pRec,
    vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
    double phi, double m2j = 0.0);
  bool map2to3IImassless(vector<Vec4>& pNew, vector<Vec4>& pRec,
    vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
    double phi);
  bool map2to3RF(vector<Vec4>& pThree, const vector<Vec4> pTwo,
    const vector<double> invariants, double phi,
    const vector<double> masses);

  // Members.

  // Pointers.
  Info*          infoPtr{};
  Settings*      settingsPtr{};
  ParticleData*  particleDataPtr{};
  Rndm*          rndmPtr{};
  Logger*        loggerPtr{};
  PartonSystems* partonSystemsPtr{};

  // Counter for output control.
  int nUnkownPDG{}, nIncorrectCol{}, nNAN{}, nVertex{}, nChargeCons{},
    nMotDau{};
  vector<int> nUnmatchedMass, nEPcons;

  // Internal flags and settings.
  bool isInitPtr{false}, isInit{false};
  int verbose{};

};

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

} // end namespace Pythia8

#endif // Pythia8_VinciaCommon_H