end-of-line normalization
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8175 / include / History.h
index 1ce1c0f..37b5de1 100644 (file)
-// History.h is a part of the PYTHIA event generator.\r
-// Copyright (C) 2013 Torbjorn Sjostrand.\r
-// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.\r
-// Please respect the MCnet Guidelines, see GUIDELINES for details.\r
-\r
-// This file is written by Stefan Prestel.\r
-// It contains the main class for matrix element merging.\r
-// Header file for the Clustering and History classes.\r
-\r
-#ifndef Pythia8_History_H\r
-#define Pythia8_History_H\r
-\r
-#include "Basics.h"\r
-#include "BeamParticle.h"\r
-#include "Event.h"\r
-#include "Info.h"\r
-#include "ParticleData.h"\r
-#include "PythiaStdlib.h"\r
-#include "Settings.h"\r
-#include "PartonLevel.h"\r
-\r
-namespace Pythia8 {\r
-\r
-//==========================================================================\r
-\r
-// Declaration of Clustering class.\r
-// This class holds information about one radiator, recoiler, emitted system.\r
-// This class is a container class for History class use.\r
-\r
-class Clustering {\r
-\r
-public:\r
-\r
-   // The emitted parton location.\r
-  int emitted;\r
-  // The emittor parton\r
-  int emittor;\r
-  // The recoiler parton.\r
-  int recoiler;\r
-  // The colour connected recoiler (Can be different for ISR)\r
-  int partner;\r
-  // The scale associated with this clustering.\r
-  double pTscale;\r
-\r
-  // Default constructor\r
-  Clustering(){\r
-    emitted   = 0;\r
-    emittor   = 0;\r
-    recoiler  = 0;\r
-    partner   = 0;\r
-    pTscale   = 0.0;\r
-  }\r
-\r
-  // Default destructor\r
-  ~Clustering(){}\r
-\r
-  // Copy constructor\r
-  Clustering( const Clustering& inSystem ){\r
-    emitted  = inSystem.emitted;\r
-    emittor  = inSystem.emittor;\r
-    recoiler = inSystem.recoiler;\r
-    partner  = inSystem.partner;\r
-    pTscale  = inSystem.pTscale;\r
-  }\r
-\r
-  // Constructor with input\r
-  Clustering( int emtIn, int radIn, int recIn, int partnerIn,\r
-    double pTscaleIn ){\r
-    emitted  = emtIn;\r
-    emittor  = radIn;\r
-    recoiler = recIn;\r
-    partner  = partnerIn;\r
-    pTscale  = pTscaleIn;\r
-  }\r
-\r
-  // Function to return pythia pT scale of clustering\r
-  double pT() const { return pTscale; }\r
-\r
-  // print for debug\r
-  void list() const;\r
-\r
-};\r
-\r
-//==========================================================================\r
-\r
-// Declaration of History class\r
-//\r
-// A History object represents an event in a given step in the CKKW-L\r
-// clustering procedure. It defines a tree-like recursive structure,\r
-// where the root node represents the state with n jets as given by\r
-// the matrix element generator, and is characterized by the member\r
-// variable mother being null. The leaves on the tree corresponds to a\r
-// fully clustered paths where the original n-jets has been clustered\r
-// down to the Born-level state. Also states which cannot be clustered\r
-// down to the Born-level are possible - these will be called\r
-// incomplete. The leaves are characterized by the vector of children\r
-// being empty.\r
-\r
-class History {\r
-\r
-public:\r
-\r
-  // The only constructor. Default arguments are used when creating\r
-  // the initial history node. The \a depth is the maximum number of\r
-  // clusterings requested. \a scalein is the scale at which the \a\r
-  // statein was clustered (should be set to the merging scale for the\r
-  // initial history node. \a beamAIn and beamBIn are needed to\r
-  // calcutate PDF ratios, \a particleDataIn to have access to the\r
-  // correct masses of particles. If \a isOrdered is true, the previous\r
-  // clusterings has been ordered. \a is the PDF ratio for this \r
-  // clustering (=1 for FSR clusterings). \a probin is the accumulated\r
-  // probabilities for the previous clusterings, and \ mothin is the\r
-  // previous history node (null for the initial node).\r
-  History( int depth,\r
-           double scalein,\r
-           Event statein,\r
-           Clustering c,\r
-           MergingHooks* mergingHooksPtrIn,\r
-           BeamParticle beamAIn,\r
-           BeamParticle beamBIn,\r
-           ParticleData* particleDataPtrIn,\r
-           Info* infoPtrIn,\r
-           bool isOrdered,\r
-           bool isStronglyOrdered,\r
-           bool isAllowed,\r
-           bool isNextInInput,\r
-           double probin,\r
-           History * mothin);\r
-\r
-  // The destructor deletes each child.\r
-  ~History() {\r
-    for ( int i = 0, N = children.size(); i < N; ++i ) delete children[i];\r
-  }\r
-\r
-  // Function to project paths onto desired paths.\r
-  bool projectOntoDesiredHistories();\r
-\r
-  // For CKKW-L, NL3 and UMEPS:\r
-  // In the initial history node, select one of the paths according to\r
-  // the probabilities. This function should be called for the initial\r
-  // history node.\r
-  // IN  trialShower*    : Previously initialised trialShower object,\r
-  //                       to perform trial showering and as\r
-  //                       repository of pointers to initialise alphaS\r
-  //     PartonSystems* : PartonSystems object needed to initialise\r
-  //                      shower objects\r
-  // OUT double         : (Sukadov) , (alpha_S ratios) , (PDF ratios)\r
-  double weightTREE(PartonLevel* trial, AlphaStrong * asFSR,\r
-                    AlphaStrong * asISR, double RN);\r
-\r
-  // For default NL3:\r
-  // Return weight of virtual correction and subtractive for NL3 merging\r
-  double weightLOOP(PartonLevel* trial, double RN);\r
-  // Return O(\alpha_s)-term of CKKWL-weight for NL3 merging\r
-  double weightFIRST(PartonLevel* trial, AlphaStrong* asFSR,\r
-                  AlphaStrong* asISR, double RN, Rndm* rndmPtr );\r
-\r
-  // For UMEPS:\r
-  double weight_UMEPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,\r
-                    AlphaStrong * asISR, double RN);\r
-  double weight_UMEPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,\r
-                    AlphaStrong * asISR, double RN);\r
-\r
-  // For unitary NL3:\r
-  double weight_UNLOPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,\r
-                    AlphaStrong * asISR, double RN);\r
-  double weight_UNLOPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,\r
-                    AlphaStrong * asISR, double RN);\r
-  double weight_UNLOPS_LOOP(PartonLevel* trial, double RN);\r
-  double weight_UNLOPS_SUBTNLO(PartonLevel* trial, double RN);\r
-  double weight_UNLOPS_CORRECTION( int order, PartonLevel* trial,\r
-                  AlphaStrong* asFSR, AlphaStrong* asISR, \r
-                  double RN, Rndm* rndmPtr );\r
-\r
-  // Function to check if any allowed histories were found\r
-  bool foundAllowedHistories() {\r
-    return (children.size() > 0 && foundAllowedPath); }\r
-  // Function to check if any ordered histories were found\r
-  bool foundOrderedHistories() {\r
-    return (children.size() > 0 && foundOrderedPath); }\r
-  // Function to check if any ordered histories were found\r
-  bool foundCompleteHistories() {\r
-    return (children.size() > 0 && foundCompletePath); }\r
-\r
-  // Function to set the state with complete scales for evolution \r
-  void getStartingConditions( const double RN, Event& outState );\r
-  // Function to get the state with complete scales for evolution \r
-  bool getClusteredEvent( const double RN, int nSteps, Event& outState);\r
-  // Function to get the first reclustered state above the merging scale. \r
-  bool getFirstClusteredEventAboveTMS( const double RN, int nDesired,\r
-    Event& process, int & nPerformed, bool updateProcess = true );\r
-  // Function to return the depth of the history (i.e. the number of\r
-  // reclustered splittings)\r
-  int nClusterings();\r
-\r
-  // Function to get the lowest multiplicity reclustered event\r
-  Event lowestMultProc( const double RN) {\r
-    // Return lowest multiplicity state\r
-    return (select(RN)->state);\r
-  }\r
-\r
-  // Calculate and return pdf ratio\r
-  double getPDFratio( int side, bool forSudakov, bool useHardPDF,\r
-                      int flavNum, double xNum, double muNum,\r
-                      int flavDen, double xDen, double muDen);\r
-\r
-  // Make Pythia class friend\r
-  friend class Pythia;\r
-  // Make Merging class friend\r
-  friend class Merging;\r
-\r
-private:\r
-\r
-  // Number of trial emission to use for calculating the average number of \r
-  // emissions\r
-  static const int NTRIAL;\r
-\r
-  // Function to set all scales in the sequence of states. This is a\r
-  // wrapper routine for setScales and setEventScales methods\r
-  void setScalesInHistory();\r
-\r
-  // Function to find the index (in the mother histories) of the\r
-  // child history, thus providing a way access the path from both\r
-  // initial history (mother == 0) and final history (all children == 0)\r
-  // IN vector<int> : The index of each child in the children vector\r
-  //                  of the current history node will be saved in\r
-  //                  this vector\r
-  // NO OUTPUT\r
-  void findPath(vector<int>& out);\r
-\r
-  // Functions to set the  parton production scales and enforce\r
-  // ordering on the scales of the respective clusterings stored in\r
-  // the History node:\r
-  // Method will start from lowest multiplicity state and move to\r
-  // higher states, setting the production scales the shower would\r
-  // have used.\r
-  // When arriving at the highest multiplicity, the method will switch\r
-  // and go back in direction of lower states to check and enforce\r
-  // ordering for unordered histories.\r
-  // IN vector<int> : Vector of positions of the chosen child\r
-  //                  in the mother history to allow to move\r
-  //                  in direction initial->final along path\r
-  //    bool        : True: Move in direction low->high\r
-  //                       multiplicity and set production scales\r
-  //                  False: Move in direction high->low\r
-  //                       multiplicity and check and enforce\r
-  //                       ordering\r
-  // NO OUTPUT\r
-  void setScales( vector<int> index, bool forward);\r
-\r
-  // Function to find a particle in all higher multiplicity events\r
-  // along the history path and set its production scale to the input\r
-  // scale\r
-  // IN  int iPart       : Parton in refEvent to be checked / rescaled\r
-  //     Event& refEvent : Reference event for iPart\r
-  //     double scale    : Scale to be set as production scale for\r
-  //                       unchanged copies of iPart in subsequent steps\r
-  void scaleCopies(int iPart, const Event& refEvent, double rho);\r
-\r
-  // Function to set the OVERALL EVENT SCALES [=state.scale()] to\r
-  // the scale of the last clustering\r
-  // NO INPUT\r
-  // NO OUTPUT\r
-  void setEventScales();\r
-\r
-  // Function to print information on the reconstructed scales in one path.\r
-  // For debug only.\r
-  void printScales() { if ( mother ) mother->printScales();\r
-    cout << " size " << state.size() << " scale " << scale << " clusterIn "\r
-      << clusterIn.pT() << " state.scale() " << state.scale() << endl; }\r
-\r
-  // Function to project paths onto desired paths.\r
-  bool trimHistories();\r
-  // Function to tag history for removal.\r
-  void remove(){ doInclude = false; }\r
-  // Function to return flag of allowed histories to choose from.\r
-  bool keep(){ return doInclude; }\r
-  // Function implementing checks on a paths, for deciding if the path should\r
-  // be considered valid or not.\r
-  bool keepHistory();\r
-  // Function to check if a path is ordered in evolution pT.\r
-  bool isOrderedPath( double maxscale );\r
-\r
-  bool followsInputPath();\r
-\r
-  // Function to check if all reconstucted states in a path pass the merging\r
-  // scale cut.\r
-  bool allIntermediateAboveRhoMS( double rhoms, bool good = true );\r
-  // Function to check if any ordered paths were found (and kept).\r
-  bool foundAnyOrderedPaths();\r
-\r
-  // Functions to return the z value of the last ISR splitting\r
-  // NO INPUT\r
-  // OUTPUT double : z value of last ISR splitting in history \r
-  double zISR();\r
-\r
-  // Functions to return the z value of the last FSR splitting\r
-  // NO INPUT\r
-  // OUTPUT double : z value of last FSR splitting in history \r
-  double zFSR();\r
-\r
-  // Functions to return the pT scale of the last ISR splitting\r
-  // NO INPUT\r
-  // OUTPUT double : pT scale of last ISR splitting in history \r
-  double pTISR();\r
-\r
-  // Functions to return the pT scale of the last FSR splitting\r
-  // NO INPUT\r
-  // OUTPUT double : pT scale of last FSR splitting in history \r
-  double pTFSR();\r
-\r
-  // Functions to return the event with nSteps additional partons\r
-  // INPUT  int   : Number of splittings in the event,\r
-  //                as counted from core 2->2 process\r
-  // OUTPUT Event : event with nSteps additional partons \r
-  Event clusteredState( int nSteps);\r
-\r
-  // Function to choose a path from all paths in the tree\r
-  // according to their splitting probabilities\r
-  // IN double    : Random number\r
-  // OUT History* : Leaf of history path chosen\r
-  History * select(double rnd);\r
-\r
-  // For a full path, find the weight calculated from the ratio of\r
-  // couplings, the no-emission probabilities, and possible PDF\r
-  // ratios. This function should only be called for the last history\r
-  // node of a full path.\r
-  // IN  TimeShower : Already initialised shower object to be used as\r
-  //                  trial shower\r
-  //     double     : alpha_s value used in ME calculation\r
-  //     double     : Maximal mass scale of the problem (e.g. E_CM)\r
-  //     AlphaStrong: Initialised shower alpha_s object for FSR alpha_s\r
-  //                  ratio calculation\r
-  //     AlphaStrong: Initialised shower alpha_s object for ISR alpha_s\r
-  //                  ratio calculation (can be different from previous)\r
-  double weightTree(PartonLevel* trial, double as0, double maxscale,\r
-    double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,\r
-    double& asWeight, double& pdfWeight);\r
-\r
-  // Function to return the \alpha_s-ratio part of the CKKWL weight.\r
-  double weightTreeALPHAS( double as0, AlphaStrong * asFSR, \r
-    AlphaStrong * asISR );\r
-  // Function to return the PDF-ratio part of the CKKWL weight.\r
-  double weightTreePDFs( double maxscale, double pdfScale );\r
-  // Function to return the no-emission probability part of the CKKWL weight.\r
-  double weightTreeEmissions( PartonLevel* trial, int type, int njetMax, \r
-    double maxscale );\r
-\r
-  // Function to generate the O(\alpha_s)-term of the CKKWL-weight\r
-  double weightFirst(PartonLevel* trial, double as0, double muR,\r
-    double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR, Rndm* rndmPtr );\r
-\r
-  // Function to generate the O(\alpha_s)-term of the \alpha_s-ratios\r
-  // appearing in the CKKWL-weight.\r
-  double weightFirstALPHAS( double as0, double muR, AlphaStrong * asFSR,\r
-    AlphaStrong * asISR);\r
-  // Function to generate the O(\alpha_s)-term of the PDF-ratios\r
-  // appearing in the CKKWL-weight.\r
-  double weightFirstPDFs( double as0, double maxscale, double pdfScale,\r
-    Rndm* rndmPtr );\r
-  // Function to generate the O(\alpha_s)-term of the no-emission\r
-  // probabilities appearing in the CKKWL-weight.\r
-  double weightFirstEmissions(PartonLevel* trial, double as0, double maxscale,\r
-    AlphaStrong * asFSR, AlphaStrong * asISR, bool fixpdf, bool fixas );\r
-\r
-  // Function to return the default factorisation scale of the hard process.\r
-  double hardFacScale(const Event& event);\r
-  // Function to return the default renormalisation scale of the hard process.\r
-  double hardRenScale(const Event& event);\r
-\r
-  // Perform a trial shower using the \a pythia object between\r
-  // maxscale down to this scale and return the corresponding Sudakov\r
-  // form factor.\r
-  // IN  trialShower : Shower object used as trial shower\r
-  //     double     : Maximum scale for trial shower branching\r
-  // OUT  0.0       : trial shower emission outside allowed pT range\r
-  //      1.0       : trial shower successful (any emission was below\r
-  //                  the minimal scale )\r
-  double doTrialShower(PartonLevel* trial, int type, double maxscale, \r
-    double minscale = 0.);\r
-\r
-  // Function to bookkeep the indices of weights generated in countEmissions\r
-  bool updateind(vector<int> & ind, int i, int N);\r
-\r
-  // Function to count number of emissions between two scales for NLO merging\r
-  vector<double> countEmissions(PartonLevel* trial, double maxscale,\r
-    double minscale, int showerType, double as0, AlphaStrong * asFSR,\r
-    AlphaStrong * asISR, int N, bool fixpdf, bool fixas);\r
-\r
-  // Function to integrate PDF ratios between two scales over x and t,\r
-  // where the PDFs are always evaluated at the lower t-integration limit\r
-  double monteCarloPDFratios(int flav, double x, double maxScale,\r
-           double minScale, double pdfScale, double asME, Rndm* rndmPtr);\r
-\r
-  // Default: Check if a ordered (and complete) path has been found in\r
-  // the initial node, in which case we will no longer be interested in\r
-  // any unordered paths.\r
-  bool onlyOrderedPaths();\r
-\r
-  // Check if a strongly ordered (and complete) path has been found in the\r
-  // initial node, in which case we will no longer be interested in\r
-  // any unordered paths.\r
-  bool onlyStronglyOrderedPaths();\r
-\r
-  // Check if an allowed (according to user-criterion) path has been found in \r
-  // the initial node, in which case we will no longer be interested in\r
-  // any forbidden paths.\r
-  bool onlyAllowedPaths();\r
-\r
-  // When a full path has been found, register it with the initial\r
-  // history node.\r
-  // IN  History : History to be registered as path\r
-  //     bool    : Specifying if clusterings so far were ordered\r
-  //     bool    : Specifying if path is complete down to 2->2 process\r
-  // OUT true if History object forms a plausible path (eg prob>0 ...)\r
-  bool registerPath(History & l, bool isOrdered, bool isStronglyOrdered,\r
-         bool isAllowed, bool isComplete);\r
-\r
-  // For the history-defining state (and if necessary interfering\r
-  // states), find all possible clusterings.\r
-  // NO INPUT\r
-  // OUT vector of all (rad,rec,emt) systems\r
-  vector<Clustering> getAllQCDClusterings();\r
-\r
-  // For one given state, find all possible clusterings.\r
-  // IN  Event : state to be investigated\r
-  // OUT vector of all (rad,rec,emt) systems in the state\r
-  vector<Clustering> getQCDClusterings( const Event& event);\r
-\r
-  // Function to construct (rad,rec,emt) triples from the event\r
-  // IN  int   : Position of Emitted in event record for which\r
-  //             dipoles should be constructed\r
-  //     int   : Colour topogy to be tested\r
-  //             1= g -> qqbar, causing 2 -> 2 dipole splitting\r
-  //             2= q(bar) -> q(bar) g && g -> gg,\r
-  //              causing a 2 -> 3 dipole splitting\r
-  //     Event : event record to be checked for ptential partners\r
-  // OUT vector of all allowed radiator+recoiler+emitted triples\r
-  vector<Clustering> findQCDTriple (int EmtTagIn, int colTopIn, \r
-                       const Event& event, vector<int> PosFinalPartn,\r
-                       vector <int> PosInitPartn );\r
-\r
-  vector<Clustering> getAllEWClusterings();\r
-  vector<Clustering> getEWClusterings( const Event& event);\r
-  vector<Clustering> findEWTriple( int EmtTagIn, const Event& event,\r
-                       vector<int> PosFinalPartn );\r
-\r
-  vector<Clustering> getAllSQCDClusterings();\r
-  vector<Clustering> getSQCDClusterings( const Event& event);\r
-  vector<Clustering> findSQCDTriple (int EmtTagIn, int colTopIn, \r
-                       const Event& event, vector<int> PosFinalPartn,\r
-                       vector <int> PosInitPartn );\r
-\r
-  // Calculate and return the probability of a clustering.\r
-  // IN  Clustering : rad,rec,emt - System for which the splitting\r
-  //                  probability should be calcuated\r
-  // OUT splitting probability\r
-  double getProb(const Clustering & SystemIn);\r
-\r
-  // Set up the beams (fill the beam particles with the correct\r
-  // current incoming particles) to allow calculation of splitting\r
-  // probability.\r
-  // For interleaved evolution, set assignments dividing PDFs into\r
-  // sea and valence content. This assignment is, until a history path\r
-  // is chosen and a first trial shower performed, not fully correct\r
-  // (since content is chosen form too high x and too low scale). The\r
-  // assignment used for reweighting will be corrected after trial\r
-  // showering\r
-  void setupBeams();\r
-\r
-  // Calculate the PDF ratio used in the argument of the no-emission\r
-  // probability.\r
-  double pdfForSudakov();\r
-\r
-  // Calculate the hard process matrix element to include in the selection\r
-  // probability.\r
-  double hardProcessME( const Event& event);\r
-\r
-  // Perform the clustering of the current state and return the\r
-  // clustered state.\r
-  // IN Clustering : rad,rec,emt triple to be clustered to two partons\r
-  // OUT clustered state\r
-  Event cluster(const Clustering & inSystem);\r
-\r
-  // Function to get the flavour of the radiator before the splitting\r
-  // for clustering\r
-  // IN  int   : Position of the radiator after the splitting, in the event\r
-  //     int   : Position of the emitted after the splitting, in the event\r
-  //     Event : Reference event   \r
-  // OUT int   : Flavour of the radiator before the splitting\r
-  int getRadBeforeFlav(const int RadAfter, const int EmtAfter,\r
-        const Event& event);\r
-\r
-  // Function to get the colour of the radiator before the splitting\r
-  // for clustering\r
-  // IN  int   : Position of the radiator after the splitting, in the event\r
-  //     int   : Position of the emitted after the splitting, in the event\r
-  //     Event : Reference event   \r
-  // OUT int   : Colour of the radiator before the splitting\r
-  int getRadBeforeCol(const int RadAfter, const int EmtAfter,\r
-        const Event& event);\r
-\r
-  // Function to get the anticolour of the radiator before the splitting\r
-  // for clustering\r
-  // IN  int   : Position of the radiator after the splitting, in the event\r
-  //     int   : Position of the emitted after the splitting, in the event\r
-  //     Event : Reference event   \r
-  // OUT int   : Anticolour of the radiator before the splitting\r
-  int getRadBeforeAcol(const int RadAfter, const int EmtAfter,\r
-        const Event& event);\r
-\r
-  // Function to get the parton connected to in by a colour line\r
-  // IN  int   : Position of parton for which partner should be found\r
-  //     Event : Reference event   \r
-  // OUT int   : If a colour line connects the "in" parton with another\r
-  //             parton, return the Position of the partner, else return 0\r
-  int getColPartner(const int in, const Event& event);\r
-\r
-  // Function to get the parton connected to in by an anticolour line\r
-  // IN  int   : Position of parton for which partner should be found\r
-  //     Event : Reference event   \r
-  // OUT int   : If an anticolour line connects the "in" parton with another\r
-  //             parton, return the Position of the partner, else return 0\r
-  int getAcolPartner(const int in, const Event& event);\r
-\r
-  // Function to get the list of partons connected to the particle\r
-  // formed by reclusterinf emt and rad by colour and anticolour lines\r
-  // IN  int          : Position of radiator in the clustering\r
-  // IN  int          : Position of emitted in the clustering\r
-  //     Event        : Reference event   \r
-  // OUT vector<int>  : List of positions of all partons that are connected\r
-  //                    to the parton that will be formed\r
-  //                    by clustering emt and rad.\r
-  vector<int> getReclusteredPartners(const int rad, const int emt,\r
-    const Event& event);\r
-\r
-  // Function to extract a chain of colour-connected partons in\r
-  // the event \r
-  // IN     int          : Type of parton from which to start extracting a\r
-  //                       parton chain. If the starting point is a quark\r
-  //                       i.e. flavType = 1, a chain of partons that are\r
-  //                       consecutively connected by colour-lines will be\r
-  //                       extracted. If the starting point is an antiquark\r
-  //                       i.e. flavType =-1, a chain of partons that are\r
-  //                       consecutively connected by anticolour-lines\r
-  //                       will be extracted.\r
-  // IN      int         : Position of the parton from which a\r
-  //                       colour-connected chain should be derived\r
-  // IN      Event       : Refernence event\r
-  // IN/OUT  vector<int> : Partons that should be excluded from the search.\r
-  // OUT     vector<int> : Positions of partons along the chain\r
-  // OUT     bool        : Found singlet / did not find singlet\r
-  bool getColSinglet(const int flavType, const int iParton,\r
-    const Event& event, vector<int>& exclude,\r
-    vector<int>& colSinglet);\r
-\r
-  // Function to check that a set of partons forms a colour singlet\r
-  // IN  Event       : Reference event\r
-  // IN  vector<int> : Positions of the partons in the set\r
-  // OUT bool        : Is a colour singlet / is not \r
-  bool isColSinglet( const Event& event, vector<int> system);\r
-  // Function to check that a set of partons forms a flavour singlet\r
-  // IN  Event       : Reference event\r
-  // IN  vector<int> : Positions of the partons in the set\r
-  // IN  int         : Flavour of all the quarks in the set, if\r
-  //                   all quarks in a set should have a fixed flavour\r
-  // OUT bool        : Is a flavour singlet / is not \r
-  bool isFlavSinglet( const Event& event,\r
-    vector<int> system, int flav=0);\r
-\r
-  // Function to properly colour-connect the radiator to the rest of\r
-  // the event, as needed during clustering\r
-  // IN  Particle& : Particle to be connected\r
-  //     Particle  : Recoiler forming a dipole with Radiator\r
-  //     Event     : event to which Radiator shall be appended\r
-  // OUT true               : Radiator could be connected to the event\r
-  //     false              : Radiator could not be connected to the\r
-  //                          event or the resulting event was\r
-  //                          non-valid\r
-  bool connectRadiator( Particle& Radiator, const int RadType,\r
-                        const Particle& Recoiler, const int RecType, \r
-                        const Event& event );\r
-\r
-  // Function to find a colour (anticolour) index in the input event\r
-  // IN  int col       : Colour tag to be investigated\r
-  //     int iExclude1 : Identifier of first particle to be excluded\r
-  //                     from search\r
-  //     int iExclude2 : Identifier of second particle to be excluded\r
-  //                     from  search\r
-  //     Event event   : event to be searched for colour tag\r
-  //     int type      : Tag to define if col should be counted as\r
-  //                      colour (type = 1) [->find anti-colour index\r
-  //                                         contracted with col]\r
-  //                      anticolour (type = 2) [->find colour index\r
-  //                                         contracted with col]\r
-  // OUT int           : Position of particle in event record\r
-  //                     contraced with col [0 if col is free tag]\r
-  int FindCol(int col, int iExclude1, int iExclude2,\r
-              const Event& event, int type, bool isHardIn);\r
-\r
-  // Function to in the input event find a particle with quantum\r
-  // numbers matching those of the input particle\r
-  // IN  Particle : Particle to be searched for\r
-  //     Event    : Event to be searched in\r
-  // OUT int      : > 0 : Position of matching particle in event\r
-  //                < 0 : No match in event\r
-  int FindParticle( const Particle& particle, const Event& event,\r
-    bool checkStatus = true );\r
-\r
-  // Function to check if rad,emt,rec triple is allowed for clustering\r
-  // IN int rad,emt,rec,partner : Positions (in event record) of the three\r
-  //                      particles considered for clustering, and the\r
-  //                      correct colour-connected recoiler (=partner)\r
-  //    Event event     : Reference event                  \r
-  bool allowedClustering( int rad, int emt, int rec, int partner,\r
-    const Event& event );\r
-\r
-  // Function to check if rad,emt,rec triple is results in\r
-  // colour singlet radBefore+recBefore\r
-  // IN int rad,emt,rec : Positions (in event record) of the three\r
-  //                      particles considered for clustering\r
-  //    Event event     : Reference event                  \r
-  bool isSinglett( int rad, int emt, int rec, const Event& event );\r
-\r
-  // Function to check if event is sensibly constructed: Meaning\r
-  // that all colour indices are contracted and that the charge in\r
-  // initial and final states matches\r
-  // IN  event : event to be checked\r
-  // OUT TRUE  : event is properly construced\r
-  //     FALSE : event not valid\r
-  bool validEvent( const Event& event );\r
-\r
-  // Function to check whether two clusterings are identical, used\r
-  // for finding the history path in the mother -> children direction\r
-  bool equalClustering( Clustering clus1 , Clustering clus2 );\r
-\r
-  // Chose dummy scale for event construction. By default, choose\r
-  //     sHat     for 2->Boson(->2)+ n partons processes and\r
-  //     M_Boson  for 2->Boson(->)             processes \r
-  double choseHardScale( const Event& event ) const;\r
-\r
-  // If the state has an incoming hadron return the flavour of the\r
-  // parton entering the hard interaction. Otherwise return 0\r
-  int getCurrentFlav(const int) const;\r
-\r
-   // If the state has an incoming hadron return the x-value for the\r
-   // parton entering the hard interaction. Otherwise return 0.\r
-  double getCurrentX(const int) const;\r
-\r
-  double getCurrentZ(const int rad, const int rec, const int emt) const;\r
-\r
-  // Function to compute "pythia pT separation" from Particle input\r
-  double pTLund(const Particle& RadAfterBranch,const Particle& EmtAfterBranch,\r
-                const Particle& RecAfterBranch, int ShowerType);\r
-\r
-  // Function to return the position of the initial line before (or after)\r
-  // a single (!) splitting.\r
-  int posChangedIncoming(const Event& event, bool before);\r
-\r
-  // Function to give back the ratio of PDFs and PDF * splitting kernels\r
-  // needed to convert a splitting at scale pdfScale, chosen with running\r
-  // PDFs, to a splitting chosen with PDFs at a fixed scale mu. As needed to\r
-  // properly count emissions.\r
-  double pdfFactor( const Event& event, const int type, double pdfScale,\r
-    double mu );\r
-\r
-  // Function giving the product of splitting kernels and PDFs so that the\r
-  // resulting flavour is given by flav. This is used as a helper routine \r
-  // to dgauss\r
-  double integrand(int flav, double x, double scaleInt, double z);\r
-\r
-  //----------------------------------------------------------------------//\r
-  // Class members.\r
-  //----------------------------------------------------------------------//\r
-\r
-  // The state of the event correponding to this step in the\r
-  // reconstruction.\r
-  Event state;\r
-\r
-  // The previous step from which this step has been clustered. If\r
-  // null, this is the initial step with the n-jet state generated by\r
-  // the matrix element.\r
-  History * mother;\r
-\r
-  // The different steps corresponding to possible clusterings of this\r
-  // state.\r
-  vector<History *> children;\r
-\r
-  // The different paths which have been reconstructed indexed with\r
-  // the (incremental) corresponding probability. This map is empty\r
-  // unless this is the initial step (mother == 0).\r
-  map<double,History *> paths;\r
-\r
-  // The sum of the probabilities of the full paths. This is zero\r
-  // unless this is the initial step (mother == 0).\r
-  double sumpath;\r
-\r
-  // The different allowed paths after projection, indexed with\r
-  // the (incremental) corresponding probability. This map is empty\r
-  // unless this is the initial step (mother == 0).\r
-  map<double,History *> goodBranches, badBranches;\r
-  // The sum of the probabilities of allowed paths after projection. This is\r
-  // zero unless this is the initial step (mother == 0).\r
-  double sumGoodBranches, sumBadBranches;\r
-\r
-  // This is set true if an ordered (and complete) path has been found\r
-  // and inserted in paths.\r
-  bool foundOrderedPath;\r
-\r
-  // This is set true if a strongly ordered (and complete) path has been found\r
-  // and inserted in paths.\r
-  bool foundStronglyOrderedPath;\r
-\r
-  // This is set true if an allowed (according to a user criterion) path has \r
-  // been found and inserted in paths.\r
-  bool foundAllowedPath;\r
-\r
-  // This is set true if a complete (with the required number of\r
-  // clusterings) path has been found and inserted in paths.\r
-  bool foundCompletePath;\r
-\r
-  // The scale of this step, corresponding to clustering which\r
-  // constructed the corresponding state (or the merging scale in case\r
-  // mother == 0).\r
-  double scale;\r
-\r
-  // Flag indicating if a clustering in the construction of all histories is\r
-  // the next clustering demanded by inout clusterings in LesHouches 2.0\r
-  // accord.\r
-  bool nextInInput;\r
-\r
-  // The probability associated with this step and the previous steps.\r
-  double prob;\r
-\r
-  // The partons and scale of the last clustering.\r
-  Clustering clusterIn;\r
-  int iReclusteredOld, iReclusteredNew;\r
-\r
-  // Flag to include the path amongst allowed paths.\r
-  bool doInclude;\r
-\r
-  // Pointer to MergingHooks object to get all the settings.\r
-  MergingHooks* mergingHooksPtr;\r
-\r
-   // The default constructor is private.\r
-  History() {}\r
-\r
-  // The copy-constructor is private.\r
-  History(const History &) {}\r
-\r
-  // The assignment operator is private.\r
-  History & operator=(const History &) {\r
-    return *this;\r
-  }\r
-\r
-  // BeamParticle to get access to PDFs\r
-  BeamParticle beamA;\r
-  BeamParticle beamB;\r
-  // ParticleData needed to initialise the shower AND to get the\r
-  // correct masses of partons needed in calculation of probability\r
-  ParticleData* particleDataPtr;\r
-\r
-  // Info object to have access to all information read from LHE file\r
-  Info* infoPtr;\r
-\r
-  // Minimal scalar sum of pT used in Herwig to choose history\r
-  double sumScalarPT;\r
-\r
-};\r
-\r
-//==========================================================================\r
-\r
-} // end namespace Pythia8\r
-\r
-#endif // end Pythia8_History_H \r
+// History.h is a part of the PYTHIA event generator.
+// Copyright (C) 2013 Torbjorn Sjostrand.
+// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
+// Please respect the MCnet Guidelines, see GUIDELINES for details.
+
+// This file is written by Stefan Prestel.
+// It contains the main class for matrix element merging.
+// Header file for the Clustering and History classes.
+
+#ifndef Pythia8_History_H
+#define Pythia8_History_H
+
+#include "Basics.h"
+#include "BeamParticle.h"
+#include "Event.h"
+#include "Info.h"
+#include "ParticleData.h"
+#include "PythiaStdlib.h"
+#include "Settings.h"
+#include "PartonLevel.h"
+
+namespace Pythia8 {
+
+//==========================================================================
+
+// Declaration of Clustering class.
+// This class holds information about one radiator, recoiler, emitted system.
+// This class is a container class for History class use.
+
+class Clustering {
+
+public:
+
+   // The emitted parton location.
+  int emitted;
+  // The emittor parton
+  int emittor;
+  // The recoiler parton.
+  int recoiler;
+  // The colour connected recoiler (Can be different for ISR)
+  int partner;
+  // The scale associated with this clustering.
+  double pTscale;
+
+  // Default constructor
+  Clustering(){
+    emitted   = 0;
+    emittor   = 0;
+    recoiler  = 0;
+    partner   = 0;
+    pTscale   = 0.0;
+  }
+
+  // Default destructor
+  ~Clustering(){}
+
+  // Copy constructor
+  Clustering( const Clustering& inSystem ){
+    emitted  = inSystem.emitted;
+    emittor  = inSystem.emittor;
+    recoiler = inSystem.recoiler;
+    partner  = inSystem.partner;
+    pTscale  = inSystem.pTscale;
+  }
+
+  // Constructor with input
+  Clustering( int emtIn, int radIn, int recIn, int partnerIn,
+    double pTscaleIn ){
+    emitted  = emtIn;
+    emittor  = radIn;
+    recoiler = recIn;
+    partner  = partnerIn;
+    pTscale  = pTscaleIn;
+  }
+
+  // Function to return pythia pT scale of clustering
+  double pT() const { return pTscale; }
+
+  // print for debug
+  void list() const;
+
+};
+
+//==========================================================================
+
+// Declaration of History class
+//
+// A History object represents an event in a given step in the CKKW-L
+// clustering procedure. It defines a tree-like recursive structure,
+// where the root node represents the state with n jets as given by
+// the matrix element generator, and is characterized by the member
+// variable mother being null. The leaves on the tree corresponds to a
+// fully clustered paths where the original n-jets has been clustered
+// down to the Born-level state. Also states which cannot be clustered
+// down to the Born-level are possible - these will be called
+// incomplete. The leaves are characterized by the vector of children
+// being empty.
+
+class History {
+
+public:
+
+  // The only constructor. Default arguments are used when creating
+  // the initial history node. The \a depth is the maximum number of
+  // clusterings requested. \a scalein is the scale at which the \a
+  // statein was clustered (should be set to the merging scale for the
+  // initial history node. \a beamAIn and beamBIn are needed to
+  // calcutate PDF ratios, \a particleDataIn to have access to the
+  // correct masses of particles. If \a isOrdered is true, the previous
+  // clusterings has been ordered. \a is the PDF ratio for this 
+  // clustering (=1 for FSR clusterings). \a probin is the accumulated
+  // probabilities for the previous clusterings, and \ mothin is the
+  // previous history node (null for the initial node).
+  History( int depth,
+           double scalein,
+           Event statein,
+           Clustering c,
+           MergingHooks* mergingHooksPtrIn,
+           BeamParticle beamAIn,
+           BeamParticle beamBIn,
+           ParticleData* particleDataPtrIn,
+           Info* infoPtrIn,
+           bool isOrdered,
+           bool isStronglyOrdered,
+           bool isAllowed,
+           bool isNextInInput,
+           double probin,
+           History * mothin);
+
+  // The destructor deletes each child.
+  ~History() {
+    for ( int i = 0, N = children.size(); i < N; ++i ) delete children[i];
+  }
+
+  // Function to project paths onto desired paths.
+  bool projectOntoDesiredHistories();
+
+  // For CKKW-L, NL3 and UMEPS:
+  // In the initial history node, select one of the paths according to
+  // the probabilities. This function should be called for the initial
+  // history node.
+  // IN  trialShower*    : Previously initialised trialShower object,
+  //                       to perform trial showering and as
+  //                       repository of pointers to initialise alphaS
+  //     PartonSystems* : PartonSystems object needed to initialise
+  //                      shower objects
+  // OUT double         : (Sukadov) , (alpha_S ratios) , (PDF ratios)
+  double weightTREE(PartonLevel* trial, AlphaStrong * asFSR,
+                    AlphaStrong * asISR, double RN);
+
+  // For default NL3:
+  // Return weight of virtual correction and subtractive for NL3 merging
+  double weightLOOP(PartonLevel* trial, double RN);
+  // Return O(\alpha_s)-term of CKKWL-weight for NL3 merging
+  double weightFIRST(PartonLevel* trial, AlphaStrong* asFSR,
+                  AlphaStrong* asISR, double RN, Rndm* rndmPtr );
+
+  // For UMEPS:
+  double weight_UMEPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,
+                    AlphaStrong * asISR, double RN);
+  double weight_UMEPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,
+                    AlphaStrong * asISR, double RN);
+
+  // For unitary NL3:
+  double weight_UNLOPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,
+                    AlphaStrong * asISR, double RN);
+  double weight_UNLOPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,
+                    AlphaStrong * asISR, double RN);
+  double weight_UNLOPS_LOOP(PartonLevel* trial, double RN);
+  double weight_UNLOPS_SUBTNLO(PartonLevel* trial, double RN);
+  double weight_UNLOPS_CORRECTION( int order, PartonLevel* trial,
+                  AlphaStrong* asFSR, AlphaStrong* asISR, 
+                  double RN, Rndm* rndmPtr );
+
+  // Function to check if any allowed histories were found
+  bool foundAllowedHistories() {
+    return (children.size() > 0 && foundAllowedPath); }
+  // Function to check if any ordered histories were found
+  bool foundOrderedHistories() {
+    return (children.size() > 0 && foundOrderedPath); }
+  // Function to check if any ordered histories were found
+  bool foundCompleteHistories() {
+    return (children.size() > 0 && foundCompletePath); }
+
+  // Function to set the state with complete scales for evolution 
+  void getStartingConditions( const double RN, Event& outState );
+  // Function to get the state with complete scales for evolution 
+  bool getClusteredEvent( const double RN, int nSteps, Event& outState);
+  // Function to get the first reclustered state above the merging scale. 
+  bool getFirstClusteredEventAboveTMS( const double RN, int nDesired,
+    Event& process, int & nPerformed, bool updateProcess = true );
+  // Function to return the depth of the history (i.e. the number of
+  // reclustered splittings)
+  int nClusterings();
+
+  // Function to get the lowest multiplicity reclustered event
+  Event lowestMultProc( const double RN) {
+    // Return lowest multiplicity state
+    return (select(RN)->state);
+  }
+
+  // Calculate and return pdf ratio
+  double getPDFratio( int side, bool forSudakov, bool useHardPDF,
+                      int flavNum, double xNum, double muNum,
+                      int flavDen, double xDen, double muDen);
+
+  // Make Pythia class friend
+  friend class Pythia;
+  // Make Merging class friend
+  friend class Merging;
+
+private:
+
+  // Number of trial emission to use for calculating the average number of 
+  // emissions
+  static const int NTRIAL;
+
+  // Function to set all scales in the sequence of states. This is a
+  // wrapper routine for setScales and setEventScales methods
+  void setScalesInHistory();
+
+  // Function to find the index (in the mother histories) of the
+  // child history, thus providing a way access the path from both
+  // initial history (mother == 0) and final history (all children == 0)
+  // IN vector<int> : The index of each child in the children vector
+  //                  of the current history node will be saved in
+  //                  this vector
+  // NO OUTPUT
+  void findPath(vector<int>& out);
+
+  // Functions to set the  parton production scales and enforce
+  // ordering on the scales of the respective clusterings stored in
+  // the History node:
+  // Method will start from lowest multiplicity state and move to
+  // higher states, setting the production scales the shower would
+  // have used.
+  // When arriving at the highest multiplicity, the method will switch
+  // and go back in direction of lower states to check and enforce
+  // ordering for unordered histories.
+  // IN vector<int> : Vector of positions of the chosen child
+  //                  in the mother history to allow to move
+  //                  in direction initial->final along path
+  //    bool        : True: Move in direction low->high
+  //                       multiplicity and set production scales
+  //                  False: Move in direction high->low
+  //                       multiplicity and check and enforce
+  //                       ordering
+  // NO OUTPUT
+  void setScales( vector<int> index, bool forward);
+
+  // Function to find a particle in all higher multiplicity events
+  // along the history path and set its production scale to the input
+  // scale
+  // IN  int iPart       : Parton in refEvent to be checked / rescaled
+  //     Event& refEvent : Reference event for iPart
+  //     double scale    : Scale to be set as production scale for
+  //                       unchanged copies of iPart in subsequent steps
+  void scaleCopies(int iPart, const Event& refEvent, double rho);
+
+  // Function to set the OVERALL EVENT SCALES [=state.scale()] to
+  // the scale of the last clustering
+  // NO INPUT
+  // NO OUTPUT
+  void setEventScales();
+
+  // Function to print information on the reconstructed scales in one path.
+  // For debug only.
+  void printScales() { if ( mother ) mother->printScales();
+    cout << " size " << state.size() << " scale " << scale << " clusterIn "
+      << clusterIn.pT() << " state.scale() " << state.scale() << endl; }
+
+  // Function to project paths onto desired paths.
+  bool trimHistories();
+  // Function to tag history for removal.
+  void remove(){ doInclude = false; }
+  // Function to return flag of allowed histories to choose from.
+  bool keep(){ return doInclude; }
+  // Function implementing checks on a paths, for deciding if the path should
+  // be considered valid or not.
+  bool keepHistory();
+  // Function to check if a path is ordered in evolution pT.
+  bool isOrderedPath( double maxscale );
+
+  bool followsInputPath();
+
+  // Function to check if all reconstucted states in a path pass the merging
+  // scale cut.
+  bool allIntermediateAboveRhoMS( double rhoms, bool good = true );
+  // Function to check if any ordered paths were found (and kept).
+  bool foundAnyOrderedPaths();
+
+  // Functions to return the z value of the last ISR splitting
+  // NO INPUT
+  // OUTPUT double : z value of last ISR splitting in history 
+  double zISR();
+
+  // Functions to return the z value of the last FSR splitting
+  // NO INPUT
+  // OUTPUT double : z value of last FSR splitting in history 
+  double zFSR();
+
+  // Functions to return the pT scale of the last ISR splitting
+  // NO INPUT
+  // OUTPUT double : pT scale of last ISR splitting in history 
+  double pTISR();
+
+  // Functions to return the pT scale of the last FSR splitting
+  // NO INPUT
+  // OUTPUT double : pT scale of last FSR splitting in history 
+  double pTFSR();
+
+  // Functions to return the event with nSteps additional partons
+  // INPUT  int   : Number of splittings in the event,
+  //                as counted from core 2->2 process
+  // OUTPUT Event : event with nSteps additional partons 
+  Event clusteredState( int nSteps);
+
+  // Function to choose a path from all paths in the tree
+  // according to their splitting probabilities
+  // IN double    : Random number
+  // OUT History* : Leaf of history path chosen
+  History * select(double rnd);
+
+  // For a full path, find the weight calculated from the ratio of
+  // couplings, the no-emission probabilities, and possible PDF
+  // ratios. This function should only be called for the last history
+  // node of a full path.
+  // IN  TimeShower : Already initialised shower object to be used as
+  //                  trial shower
+  //     double     : alpha_s value used in ME calculation
+  //     double     : Maximal mass scale of the problem (e.g. E_CM)
+  //     AlphaStrong: Initialised shower alpha_s object for FSR alpha_s
+  //                  ratio calculation
+  //     AlphaStrong: Initialised shower alpha_s object for ISR alpha_s
+  //                  ratio calculation (can be different from previous)
+  double weightTree(PartonLevel* trial, double as0, double maxscale,
+    double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,
+    double& asWeight, double& pdfWeight);
+
+  // Function to return the \alpha_s-ratio part of the CKKWL weight.
+  double weightTreeALPHAS( double as0, AlphaStrong * asFSR, 
+    AlphaStrong * asISR );
+  // Function to return the PDF-ratio part of the CKKWL weight.
+  double weightTreePDFs( double maxscale, double pdfScale );
+  // Function to return the no-emission probability part of the CKKWL weight.
+  double weightTreeEmissions( PartonLevel* trial, int type, int njetMax, 
+    double maxscale );
+
+  // Function to generate the O(\alpha_s)-term of the CKKWL-weight
+  double weightFirst(PartonLevel* trial, double as0, double muR,
+    double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR, Rndm* rndmPtr );
+
+  // Function to generate the O(\alpha_s)-term of the \alpha_s-ratios
+  // appearing in the CKKWL-weight.
+  double weightFirstALPHAS( double as0, double muR, AlphaStrong * asFSR,
+    AlphaStrong * asISR);
+  // Function to generate the O(\alpha_s)-term of the PDF-ratios
+  // appearing in the CKKWL-weight.
+  double weightFirstPDFs( double as0, double maxscale, double pdfScale,
+    Rndm* rndmPtr );
+  // Function to generate the O(\alpha_s)-term of the no-emission
+  // probabilities appearing in the CKKWL-weight.
+  double weightFirstEmissions(PartonLevel* trial, double as0, double maxscale,
+    AlphaStrong * asFSR, AlphaStrong * asISR, bool fixpdf, bool fixas );
+
+  // Function to return the default factorisation scale of the hard process.
+  double hardFacScale(const Event& event);
+  // Function to return the default renormalisation scale of the hard process.
+  double hardRenScale(const Event& event);
+
+  // Perform a trial shower using the \a pythia object between
+  // maxscale down to this scale and return the corresponding Sudakov
+  // form factor.
+  // IN  trialShower : Shower object used as trial shower
+  //     double     : Maximum scale for trial shower branching
+  // OUT  0.0       : trial shower emission outside allowed pT range
+  //      1.0       : trial shower successful (any emission was below
+  //                  the minimal scale )
+  double doTrialShower(PartonLevel* trial, int type, double maxscale, 
+    double minscale = 0.);
+
+  // Function to bookkeep the indices of weights generated in countEmissions
+  bool updateind(vector<int> & ind, int i, int N);
+
+  // Function to count number of emissions between two scales for NLO merging
+  vector<double> countEmissions(PartonLevel* trial, double maxscale,
+    double minscale, int showerType, double as0, AlphaStrong * asFSR,
+    AlphaStrong * asISR, int N, bool fixpdf, bool fixas);
+
+  // Function to integrate PDF ratios between two scales over x and t,
+  // where the PDFs are always evaluated at the lower t-integration limit
+  double monteCarloPDFratios(int flav, double x, double maxScale,
+           double minScale, double pdfScale, double asME, Rndm* rndmPtr);
+
+  // Default: Check if a ordered (and complete) path has been found in
+  // the initial node, in which case we will no longer be interested in
+  // any unordered paths.
+  bool onlyOrderedPaths();
+
+  // Check if a strongly ordered (and complete) path has been found in the
+  // initial node, in which case we will no longer be interested in
+  // any unordered paths.
+  bool onlyStronglyOrderedPaths();
+
+  // Check if an allowed (according to user-criterion) path has been found in 
+  // the initial node, in which case we will no longer be interested in
+  // any forbidden paths.
+  bool onlyAllowedPaths();
+
+  // When a full path has been found, register it with the initial
+  // history node.
+  // IN  History : History to be registered as path
+  //     bool    : Specifying if clusterings so far were ordered
+  //     bool    : Specifying if path is complete down to 2->2 process
+  // OUT true if History object forms a plausible path (eg prob>0 ...)
+  bool registerPath(History & l, bool isOrdered, bool isStronglyOrdered,
+         bool isAllowed, bool isComplete);
+
+  // For the history-defining state (and if necessary interfering
+  // states), find all possible clusterings.
+  // NO INPUT
+  // OUT vector of all (rad,rec,emt) systems
+  vector<Clustering> getAllQCDClusterings();
+
+  // For one given state, find all possible clusterings.
+  // IN  Event : state to be investigated
+  // OUT vector of all (rad,rec,emt) systems in the state
+  vector<Clustering> getQCDClusterings( const Event& event);
+
+  // Function to construct (rad,rec,emt) triples from the event
+  // IN  int   : Position of Emitted in event record for which
+  //             dipoles should be constructed
+  //     int   : Colour topogy to be tested
+  //             1= g -> qqbar, causing 2 -> 2 dipole splitting
+  //             2= q(bar) -> q(bar) g && g -> gg,
+  //              causing a 2 -> 3 dipole splitting
+  //     Event : event record to be checked for ptential partners
+  // OUT vector of all allowed radiator+recoiler+emitted triples
+  vector<Clustering> findQCDTriple (int EmtTagIn, int colTopIn, 
+                       const Event& event, vector<int> PosFinalPartn,
+                       vector <int> PosInitPartn );
+
+  vector<Clustering> getAllEWClusterings();
+  vector<Clustering> getEWClusterings( const Event& event);
+  vector<Clustering> findEWTriple( int EmtTagIn, const Event& event,
+                       vector<int> PosFinalPartn );
+
+  vector<Clustering> getAllSQCDClusterings();
+  vector<Clustering> getSQCDClusterings( const Event& event);
+  vector<Clustering> findSQCDTriple (int EmtTagIn, int colTopIn, 
+                       const Event& event, vector<int> PosFinalPartn,
+                       vector <int> PosInitPartn );
+
+  // Calculate and return the probability of a clustering.
+  // IN  Clustering : rad,rec,emt - System for which the splitting
+  //                  probability should be calcuated
+  // OUT splitting probability
+  double getProb(const Clustering & SystemIn);
+
+  // Set up the beams (fill the beam particles with the correct
+  // current incoming particles) to allow calculation of splitting
+  // probability.
+  // For interleaved evolution, set assignments dividing PDFs into
+  // sea and valence content. This assignment is, until a history path
+  // is chosen and a first trial shower performed, not fully correct
+  // (since content is chosen form too high x and too low scale). The
+  // assignment used for reweighting will be corrected after trial
+  // showering
+  void setupBeams();
+
+  // Calculate the PDF ratio used in the argument of the no-emission
+  // probability.
+  double pdfForSudakov();
+
+  // Calculate the hard process matrix element to include in the selection
+  // probability.
+  double hardProcessME( const Event& event);
+
+  // Perform the clustering of the current state and return the
+  // clustered state.
+  // IN Clustering : rad,rec,emt triple to be clustered to two partons
+  // OUT clustered state
+  Event cluster(const Clustering & inSystem);
+
+  // Function to get the flavour of the radiator before the splitting
+  // for clustering
+  // IN  int   : Position of the radiator after the splitting, in the event
+  //     int   : Position of the emitted after the splitting, in the event
+  //     Event : Reference event   
+  // OUT int   : Flavour of the radiator before the splitting
+  int getRadBeforeFlav(const int RadAfter, const int EmtAfter,
+        const Event& event);
+
+  // Function to get the colour of the radiator before the splitting
+  // for clustering
+  // IN  int   : Position of the radiator after the splitting, in the event
+  //     int   : Position of the emitted after the splitting, in the event
+  //     Event : Reference event   
+  // OUT int   : Colour of the radiator before the splitting
+  int getRadBeforeCol(const int RadAfter, const int EmtAfter,
+        const Event& event);
+
+  // Function to get the anticolour of the radiator before the splitting
+  // for clustering
+  // IN  int   : Position of the radiator after the splitting, in the event
+  //     int   : Position of the emitted after the splitting, in the event
+  //     Event : Reference event   
+  // OUT int   : Anticolour of the radiator before the splitting
+  int getRadBeforeAcol(const int RadAfter, const int EmtAfter,
+        const Event& event);
+
+  // Function to get the parton connected to in by a colour line
+  // IN  int   : Position of parton for which partner should be found
+  //     Event : Reference event   
+  // OUT int   : If a colour line connects the "in" parton with another
+  //             parton, return the Position of the partner, else return 0
+  int getColPartner(const int in, const Event& event);
+
+  // Function to get the parton connected to in by an anticolour line
+  // IN  int   : Position of parton for which partner should be found
+  //     Event : Reference event   
+  // OUT int   : If an anticolour line connects the "in" parton with another
+  //             parton, return the Position of the partner, else return 0
+  int getAcolPartner(const int in, const Event& event);
+
+  // Function to get the list of partons connected to the particle
+  // formed by reclusterinf emt and rad by colour and anticolour lines
+  // IN  int          : Position of radiator in the clustering
+  // IN  int          : Position of emitted in the clustering
+  //     Event        : Reference event   
+  // OUT vector<int>  : List of positions of all partons that are connected
+  //                    to the parton that will be formed
+  //                    by clustering emt and rad.
+  vector<int> getReclusteredPartners(const int rad, const int emt,
+    const Event& event);
+
+  // Function to extract a chain of colour-connected partons in
+  // the event 
+  // IN     int          : Type of parton from which to start extracting a
+  //                       parton chain. If the starting point is a quark
+  //                       i.e. flavType = 1, a chain of partons that are
+  //                       consecutively connected by colour-lines will be
+  //                       extracted. If the starting point is an antiquark
+  //                       i.e. flavType =-1, a chain of partons that are
+  //                       consecutively connected by anticolour-lines
+  //                       will be extracted.
+  // IN      int         : Position of the parton from which a
+  //                       colour-connected chain should be derived
+  // IN      Event       : Refernence event
+  // IN/OUT  vector<int> : Partons that should be excluded from the search.
+  // OUT     vector<int> : Positions of partons along the chain
+  // OUT     bool        : Found singlet / did not find singlet
+  bool getColSinglet(const int flavType, const int iParton,
+    const Event& event, vector<int>& exclude,
+    vector<int>& colSinglet);
+
+  // Function to check that a set of partons forms a colour singlet
+  // IN  Event       : Reference event
+  // IN  vector<int> : Positions of the partons in the set
+  // OUT bool        : Is a colour singlet / is not 
+  bool isColSinglet( const Event& event, vector<int> system);
+  // Function to check that a set of partons forms a flavour singlet
+  // IN  Event       : Reference event
+  // IN  vector<int> : Positions of the partons in the set
+  // IN  int         : Flavour of all the quarks in the set, if
+  //                   all quarks in a set should have a fixed flavour
+  // OUT bool        : Is a flavour singlet / is not 
+  bool isFlavSinglet( const Event& event,
+    vector<int> system, int flav=0);
+
+  // Function to properly colour-connect the radiator to the rest of
+  // the event, as needed during clustering
+  // IN  Particle& : Particle to be connected
+  //     Particle  : Recoiler forming a dipole with Radiator
+  //     Event     : event to which Radiator shall be appended
+  // OUT true               : Radiator could be connected to the event
+  //     false              : Radiator could not be connected to the
+  //                          event or the resulting event was
+  //                          non-valid
+  bool connectRadiator( Particle& Radiator, const int RadType,
+                        const Particle& Recoiler, const int RecType, 
+                        const Event& event );
+
+  // Function to find a colour (anticolour) index in the input event
+  // IN  int col       : Colour tag to be investigated
+  //     int iExclude1 : Identifier of first particle to be excluded
+  //                     from search
+  //     int iExclude2 : Identifier of second particle to be excluded
+  //                     from  search
+  //     Event event   : event to be searched for colour tag
+  //     int type      : Tag to define if col should be counted as
+  //                      colour (type = 1) [->find anti-colour index
+  //                                         contracted with col]
+  //                      anticolour (type = 2) [->find colour index
+  //                                         contracted with col]
+  // OUT int           : Position of particle in event record
+  //                     contraced with col [0 if col is free tag]
+  int FindCol(int col, int iExclude1, int iExclude2,
+              const Event& event, int type, bool isHardIn);
+
+  // Function to in the input event find a particle with quantum
+  // numbers matching those of the input particle
+  // IN  Particle : Particle to be searched for
+  //     Event    : Event to be searched in
+  // OUT int      : > 0 : Position of matching particle in event
+  //                < 0 : No match in event
+  int FindParticle( const Particle& particle, const Event& event,
+    bool checkStatus = true );
+
+  // Function to check if rad,emt,rec triple is allowed for clustering
+  // IN int rad,emt,rec,partner : Positions (in event record) of the three
+  //                      particles considered for clustering, and the
+  //                      correct colour-connected recoiler (=partner)
+  //    Event event     : Reference event                  
+  bool allowedClustering( int rad, int emt, int rec, int partner,
+    const Event& event );
+
+  // Function to check if rad,emt,rec triple is results in
+  // colour singlet radBefore+recBefore
+  // IN int rad,emt,rec : Positions (in event record) of the three
+  //                      particles considered for clustering
+  //    Event event     : Reference event                  
+  bool isSinglett( int rad, int emt, int rec, const Event& event );
+
+  // Function to check if event is sensibly constructed: Meaning
+  // that all colour indices are contracted and that the charge in
+  // initial and final states matches
+  // IN  event : event to be checked
+  // OUT TRUE  : event is properly construced
+  //     FALSE : event not valid
+  bool validEvent( const Event& event );
+
+  // Function to check whether two clusterings are identical, used
+  // for finding the history path in the mother -> children direction
+  bool equalClustering( Clustering clus1 , Clustering clus2 );
+
+  // Chose dummy scale for event construction. By default, choose
+  //     sHat     for 2->Boson(->2)+ n partons processes and
+  //     M_Boson  for 2->Boson(->)             processes 
+  double choseHardScale( const Event& event ) const;
+
+  // If the state has an incoming hadron return the flavour of the
+  // parton entering the hard interaction. Otherwise return 0
+  int getCurrentFlav(const int) const;
+
+   // If the state has an incoming hadron return the x-value for the
+   // parton entering the hard interaction. Otherwise return 0.
+  double getCurrentX(const int) const;
+
+  double getCurrentZ(const int rad, const int rec, const int emt) const;
+
+  // Function to compute "pythia pT separation" from Particle input
+  double pTLund(const Particle& RadAfterBranch,const Particle& EmtAfterBranch,
+                const Particle& RecAfterBranch, int ShowerType);
+
+  // Function to return the position of the initial line before (or after)
+  // a single (!) splitting.
+  int posChangedIncoming(const Event& event, bool before);
+
+  // Function to give back the ratio of PDFs and PDF * splitting kernels
+  // needed to convert a splitting at scale pdfScale, chosen with running
+  // PDFs, to a splitting chosen with PDFs at a fixed scale mu. As needed to
+  // properly count emissions.
+  double pdfFactor( const Event& event, const int type, double pdfScale,
+    double mu );
+
+  // Function giving the product of splitting kernels and PDFs so that the
+  // resulting flavour is given by flav. This is used as a helper routine 
+  // to dgauss
+  double integrand(int flav, double x, double scaleInt, double z);
+
+  //----------------------------------------------------------------------//
+  // Class members.
+  //----------------------------------------------------------------------//
+
+  // The state of the event correponding to this step in the
+  // reconstruction.
+  Event state;
+
+  // The previous step from which this step has been clustered. If
+  // null, this is the initial step with the n-jet state generated by
+  // the matrix element.
+  History * mother;
+
+  // The different steps corresponding to possible clusterings of this
+  // state.
+  vector<History *> children;
+
+  // The different paths which have been reconstructed indexed with
+  // the (incremental) corresponding probability. This map is empty
+  // unless this is the initial step (mother == 0).
+  map<double,History *> paths;
+
+  // The sum of the probabilities of the full paths. This is zero
+  // unless this is the initial step (mother == 0).
+  double sumpath;
+
+  // The different allowed paths after projection, indexed with
+  // the (incremental) corresponding probability. This map is empty
+  // unless this is the initial step (mother == 0).
+  map<double,History *> goodBranches, badBranches;
+  // The sum of the probabilities of allowed paths after projection. This is
+  // zero unless this is the initial step (mother == 0).
+  double sumGoodBranches, sumBadBranches;
+
+  // This is set true if an ordered (and complete) path has been found
+  // and inserted in paths.
+  bool foundOrderedPath;
+
+  // This is set true if a strongly ordered (and complete) path has been found
+  // and inserted in paths.
+  bool foundStronglyOrderedPath;
+
+  // This is set true if an allowed (according to a user criterion) path has 
+  // been found and inserted in paths.
+  bool foundAllowedPath;
+
+  // This is set true if a complete (with the required number of
+  // clusterings) path has been found and inserted in paths.
+  bool foundCompletePath;
+
+  // The scale of this step, corresponding to clustering which
+  // constructed the corresponding state (or the merging scale in case
+  // mother == 0).
+  double scale;
+
+  // Flag indicating if a clustering in the construction of all histories is
+  // the next clustering demanded by inout clusterings in LesHouches 2.0
+  // accord.
+  bool nextInInput;
+
+  // The probability associated with this step and the previous steps.
+  double prob;
+
+  // The partons and scale of the last clustering.
+  Clustering clusterIn;
+  int iReclusteredOld, iReclusteredNew;
+
+  // Flag to include the path amongst allowed paths.
+  bool doInclude;
+
+  // Pointer to MergingHooks object to get all the settings.
+  MergingHooks* mergingHooksPtr;
+
+   // The default constructor is private.
+  History() {}
+
+  // The copy-constructor is private.
+  History(const History &) {}
+
+  // The assignment operator is private.
+  History & operator=(const History &) {
+    return *this;
+  }
+
+  // BeamParticle to get access to PDFs
+  BeamParticle beamA;
+  BeamParticle beamB;
+  // ParticleData needed to initialise the shower AND to get the
+  // correct masses of partons needed in calculation of probability
+  ParticleData* particleDataPtr;
+
+  // Info object to have access to all information read from LHE file
+  Info* infoPtr;
+
+  // Minimal scalar sum of pT used in Herwig to choose history
+  double sumScalarPT;
+
+};
+
+//==========================================================================
+
+} // end namespace Pythia8
+
+#endif // end Pythia8_History_H