-// 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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff