VinciaCommon.h

// VinciaCommon.h is a part of the PYTHIA event generator.
// Copyright (C) 2020 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 clashes.
namespace Vincia {

// Global numerical precision target (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 TINY     = 1.0e-9;
const double TINYMASS = 1.0e-6;
const double SMALL    = 1.0e-3;

// Color 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;

// Indices of FF antenna functions.
const int iQQemitFF  = 0;
const int iQGemitFF  = 1;
const int iGQemitFF  = 2;
const int iGGemitFF  = 3;
const int iGXsplitFF = 4;

// Indices of RF antenna functions.
const int iQQemitRF  = 5;
const int iQGemitRF  = 6;
const int iXGsplitRF = 7;

// Indices of II antenna functions.
const int iQQemitII  = 0;
const int iGQemitII  = 1;
const int iGGemitII  = 2;
const int iQXsplitII = 3;
const int iGXconvII  = 4;

// Indices of IF antenna functions.
const int iQQemitIF  = 5;
const int iQGemitIF  = 6;
const int iGQemitIF  = 7;
const int iGGemitIF  = 8;
const int iQXsplitIF = 9;
const int iGXconvIF  = 10;
const int iXGsplitIF = 11;

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

// Verbosity levels.
// Suppressed verbosity.
const int silent        = 0;
const int quiet         = 1;
// Normal verbosity for warnings.
const int normal        = 2;
// Extra verbosity for warnings.
const int quiteloud     = 3;
const int loud          = 4;
const int veryloud      = 5;
// Debug verbosity levels.
const int debug         = 6;
const int louddebug     = 7;
const int verylouddebug = 8;
const int superdebug    = 9;

}

// Everything else lives in the Pythia 8 namespace.

namespace Pythia8 {

// Include namespace with global Vincia constants.
using namespace Vincia;

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

// 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) {
  size_t colons = prettyFunction.find("::");
  size_t 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__

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

// Functors are required to use zbrent for the rescaling variable
// needed in the massive Rambo version.

// Pure virtual function overloading operator().
class TFunctor {

public:
  virtual double operator()(double) = 0;
  virtual ~TFunctor() {}

};

// Functor for xi-solution in Rambo.
class TXiFunctor : public TFunctor {

public:

  // Constructor takes a vector with the masses and one with the energies.
  TXiFunctor(vector<double> mIn, vector<double> energiesIn);
  double operator() (double);

private:
  vector<double> m;
  vector<double> energies;

};

// Derive a class which can wrap function pointers.
class TPtrFunctor : public TFunctor {

public:
  TPtrFunctor(double (*fPtrIn)(double) ){ fPtr = fPtrIn; };
  double operator() (double arg) { return fPtr(arg); };

private:

  // Pointer to a function taking one double and returning one double.
  double (*fPtr)(double);

};

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

// Global functions accessible in the Vincia namespace.

// Utilities for printing out VINCIA Messages.
void printOut(string,string);

// 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 few useful auxiliary functions like extra invariants and dot products.
double m          (const Vec4&);
double m2         (const Vec4&);
double m2         (const Vec4&, const Vec4&, const Vec4&);
double m2         (const Vec4&, const Vec4&, const Vec4&, const Vec4&);
double m2         (const Particle&, const Particle&, const Particle&);
double dot4       (const Particle&, const Particle&);
double getCosTheta(double E1, double E2, double m1, double m2, double s12);

// Gram determinant, invariants used in the argument = 2*pi*pj.
double gramDet(double s01tilde, double s12tilde, double s02tilde,
  double m0, double m1, double m2);
double gramDet(Vec4 p0, Vec4 p1, Vec4 p2);

// Math support functions.
double Li2      (const double,const double kmax = 100.0,
                 const double xerr = Vincia::TINY);
double factorial(const int);
int    binomial (const int,int);
double LambertW (const double x);

// Zero-finder zbrent with Functor (used by massive Rambo). Note,
// does not work if the reference to TFunctor is declared const.
double zbrent(TFunctor&, double, double, double, double);

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

// Rambo flat phase-space generator.

// This is an implementation of the Rambo phase-space generator as
// presented in A New Monte Carlo Treatment Of Multiparticle Phase
// Space At High-Energies, R. Kleiss, W.J. Stirling, S.D. Ellis, CPC40
// (1986) 359.

class Rambo {

 public:

  // Deafult constructor.
  Rambo() { rndmPtr=nullptr; isInitPtr=false;}

  // Initializing constructor.
  Rambo(Rndm* rndmPtrIn) { initPtr(rndmPtrIn); }

  // Destructor.
  virtual ~Rambo() {}

  // Initialize pointers.
  void initPtr(Rndm* rndmPtrIn) {rndmPtr = rndmPtrIn; isInitPtr = true;}

  // Rambo phase space generator. Generates nOut uniformly distributed
  // massless 4-vectors with sqrt(s) = eCM. Output in pOut.
  double genPoint(double eCM,int nOut,vector<Vec4>& pOut);

  // Massive generalisation, weights NOT 1 anymore - literal implementation
  // of original RAMBO paper by Ellis, Kleiss and Stirling. Number of particles
  // determined from size of mIn vector.
  double genPoint(double eCM,vector<double> mIn,vector<Vec4>& pOut);

 private:

  // Is initialized.
  bool isInitPtr;

  // Pointer to the random number generator.
  Rndm*  rndmPtr;

};

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

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

class Colour {

public:

  // Constructor.

  Colour() {isInitPtr=false; isInit=false;}

  // Destructor.
  ~Colour() {}

  // 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, isInit;

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

  // Verbose level.
  int verbose;

};

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

// A simple class for containing evolution variable definitions.

class Resolution {

public:

  // Constructor.
  Resolution() = default;

  // Destructor.
  virtual ~Resolution() {};

  // Initialize pointers (must be done before init).
  void initPtr(Settings* settingsPtrIn) {
    settingsPtr = settingsPtrIn;
    isInitPtr   = true;
  }

  // Initialize.
  bool init();

  // Sector resolution functions.
  double q2sector2to3(const Particle* a, const Particle* b, const Particle* j,
    bool = false);

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

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

  // Sector resolution function for 3->5 branchings (emission +
  // splitting).
  double q2sector3to5(Particle* a, Particle* b,
    Particle* j1, Particle* j2, Particle* j3);

  // Sector accept function. Optionally prevent g->qq clusterings if
  // that would reduce the number of fermion lines below some minimum
  // (cheap way to indicate that Z->qq/ll always has at least one
  // fermion pair).
  double findSector(vector<int>& iSec, vector<Particle> state, int nFmin = 1);

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

private:

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

  // Pointer to PYTHIA 8 settings database.
  Settings* settingsPtr{};

  // 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:

  // Constructor.
  VinciaCommon() {isInitPtr = false; isInit = false;}

  // Destructor.
  virtual ~VinciaCommon() {}

  // Initialize pointers.
  bool initPtr(Info* infoPtrIn) {
    infoPtr      = infoPtrIn;
    particleDataPtr = infoPtr->particleDataPtr;
    settingsPtr     = infoPtr->settingsPtr;
    rndmPtr         = infoPtr->rndmPtr;
    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);

  // More lightweight function to check conservation of momentum.
  bool checkCoM(int iSys, Event& event,PartonSystems* partonSystemsPtr);

  // 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, PartonSystems* partonSystemsPtr,
    const Event& event, double sbbSav);

  // 3->2 clustering maps.
  bool map3to2FFmassive(vector<Vec4>& pClu, vector<Vec4> pIn,
    int kMapType, int a=0, int r=1, int b=2, double mI=0.0, double mK=0.0);
  bool map3to2FFmassless(vector<Vec4>& pClu, vector<Vec4> pIn,
    int kMapType, int a=0, int r=1, int b=2);
  bool map3to2IFmassive(vector<Vec4>& pClu, vector<Vec4>& pIn,
    double saj, double sjk, double sak);
  bool map3to2IImassive(vector<Vec4>& pClu, vector<Vec4>& pIn,
    vector<Vec4>& pRec, double saj, double sjb, double sab, bool doBoost);

  // 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);
    }
  }
  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);

  // 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);
  bool map2to3IImassless(vector<Vec4>& pNew, vector<Vec4>& pRec,
    vector<Vec4>& pOld, double sAB, double saj, double sjb, double sab,
    double phi);

  // 2->3 kinematics maps for IF branchings. General massive case
  // implemented by Verheyen.
  bool map2to3IFlocal(vector<Vec4>& pNew, 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,
    vector<Vec4>& pOld, Vec4 &pB,
    double sAK, double saj, double sjk, double sak, double phi,
    double mK2, double mj2, double mk2);

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

  // 2->3 kinematics map for RF branchings. Original implementation by Brooks.
  bool map2to3RFmassive(vector<Vec4>& pNew, vector<Vec4>& pRec,
    vector<Vec4> pOld, double saj, double sjk, double phi,
    double m2oldA, double m2j, double m2oldK);

  // 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, PartonSystems* partonSystemsPtr,
                     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));
  }

  // 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:

  // Pointers.
  Info*          infoPtr;
  Settings*      settingsPtr;
  ParticleData*  particleDataPtr;
  Rndm*          rndmPtr;

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

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

};

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

} // end namespace Pythia8

#endif // Pythia8_VinciaCommon_H