EIC code displayed by LXR master/​include/​Pythia8/​History.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-18 10:06:25

 // History.h is a part of the PYTHIA event generator.
 // Copyright (C) 2024 Torbjorn Sjostrand.
 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
 // Please respect the MCnet Guidelines, see GUIDELINES for details.
 
 // This file 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 "Pythia8/Basics.h"
 #include "Pythia8/BeamParticle.h"
 #include "Pythia8/Event.h"
 #include "Pythia8/Info.h"
 #include "Pythia8/ParticleData.h"
 #include "Pythia8/PartonLevel.h"
 #include "Pythia8/PythiaStdlib.h"
 #include "Pythia8/Settings.h"
 #include "Pythia8/StandardModel.h"
 #include "Pythia8/SimpleWeakShowerMEs.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;
   // The flavour of the radiator prior to the emission.
   int flavRadBef;
   // Spin of the radiator (-1 is left handed, +1 is right handed).
   int spinRad;
   // Spin of the emitted  (-1 is left handed, +1 is right handed).
   int spinEmt;
   // Spin of the recoiler (-1 is left handed, +1 is right handed).
   int spinRec;
   // Spin of the radiator before the splitting.
   int spinRadBef;
   // The radiator before the splitting.
   int radBef;
   // The recoiler before the splitting.
   int recBef;
 
   // Map between particle positions in the clustered states -> particle
   // positions in unclustered real-emission state.
   map<int,int> iPosInMother;
 
   // Default constructor
   Clustering() : emitted(0), emittor(0), recoiler(0), partner(0), pTscale(),
     flavRadBef(0), spinRad(9), spinEmt(9), spinRec(9), spinRadBef(9),
     radBef(0), recBef(0), iPosInMother() {}
 
   // Default destructor
   ~Clustering(){}
 
   // Copy constructor
   Clustering( const Clustering& inSystem ) :
     emitted(inSystem.emitted), emittor(inSystem.emittor),
     recoiler(inSystem.recoiler), partner(inSystem.partner),
     pTscale(inSystem.pTscale), flavRadBef(inSystem.flavRadBef),
     spinRad(inSystem.spinRad), spinEmt(inSystem.spinEmt),
     spinRec(inSystem.spinRec), spinRadBef(inSystem.spinRad),
     radBef(inSystem.radBef), recBef(inSystem.recBef),
     iPosInMother(inSystem.iPosInMother) {}
 
   // Constructor with input
   Clustering( int emtIn, int radIn, int recIn, int partnerIn,
     double pTscaleIn, int flavRadBefIn = 0, int spinRadIn = 9,
     int spinEmtIn = 9, int spinRecIn = 9, int spinRadBefIn = 9,
     int radBefIn = 0, int recBefIn = 0, map<int,int> posIn = map<int,int>())
     : emitted(emtIn), emittor(radIn), recoiler(recIn), partner(partnerIn),
       pTscale(pTscaleIn), flavRadBef(flavRadBefIn), spinRad(spinRadIn),
       spinEmt(spinEmtIn), spinRec(spinRecIn), spinRadBef(spinRadBefIn),
       radBef(radBefIn), recBef(recBefIn), iPosInMother(posIn) {}
 
   // 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 depthIn,
            double scalein,
            Event statein,
            Clustering c,
            MergingHooksPtr mergingHooksPtrIn,
            BeamParticle beamAIn,
            BeamParticle beamBIn,
            ParticleData* particleDataPtrIn,
            Info* infoPtrIn,
            PartonLevel* showersIn,
            CoupSM* coupSMPtrIn,
            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 vector<double>  : (Sukadov) , (alpha_S ratios) , (PDF ratios)
   vector<double> weightCKKWL(PartonLevel* trial, AlphaStrong * asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
 
 
   // For default NL3:
   // Return weight of virtual correction and subtractive for NL3 merging
   vector<double> weightNL3Loop(PartonLevel* trial, double RN);
   // Return O(\alpha_s)-term of CKKWL-weight for NL3 merging
   vector<double> weightNL3First(PartonLevel* trial, AlphaStrong* asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
     Rndm* rndmPtr);
   vector<double> weightNL3Tree(PartonLevel* trial, AlphaStrong * asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
 
 
   // For UMEPS:
   vector<double> weightUMEPSTree(PartonLevel* trial, AlphaStrong * asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
   vector<double> weightUMEPSSubt(PartonLevel* trial, AlphaStrong * asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
 
 
   // For unitary NL3:
   vector<double> weightUNLOPSTree(PartonLevel* trial, AlphaStrong * asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
     int depthIn = -1);
   vector<double> weightUNLOPSSubt(PartonLevel* trial, AlphaStrong * asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
     int depthIn = -1);
   vector<double> weightUNLOPSLoop(PartonLevel* trial, AlphaStrong * asFSR,
      AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
      int depthIn = -1);
   vector<double> weightUNLOPSSubtNLO(PartonLevel* trial, AlphaStrong * asFSR,
     AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
     int depthIn = -1);
   vector<double> weightUNLOPSFirst( int order, PartonLevel* trial,
     AlphaStrong* asFSR, AlphaStrong * asISR, AlphaEM * aemFSR,
     AlphaEM * aemISR, 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);
 
   // Envelope function that calls the recursive getWeakProb.
   double getWeakProb();
 
   // Recursive function that returns the weak probability for the given path.
   // Mode refers to which ME correction to use, 1 = sChannel, 2 = gluon
   // channel, 3 = double quark t-channel, 4 is double quark u-channel.
   double getWeakProb(vector<int>& mode, vector<Vec4>& mom,
      vector<int> fermionLines);
 
   // return the weak probability of a single reclustering.
   // Mode refers to which ME correction to use, 1 = sChannel, 2 = gluon
   // channel, 3 = double quark t-channel, 4 is double quark u-channel.
   double getSingleWeakProb(vector<int> &mode, vector<Vec4> &mom,
     vector<int> fermionLines);
 
   // Find map between indecies in the current state and the state after
   // the splitting.
   // NOT IMPLEMENTED FOR MULTIPLE W/Z/GAMMA (NEED TO HAVE A WAY TO
   // IDENTIFY THEM).
   void findStateTransfer(map<int,int> &transfer);
 
   // Function to print the history that would be chosen from the random number
   // RN. Mainly for debugging.
   void printHistory( const double RN );
   // Function to print the states in a history, starting from the hard process.
   // Mainly for debugging.
   void printStates();
 
   // 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)
   vector<double> weightTree(PartonLevel* trial, double as0, double aem0,
     double maxscale, double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,
     AlphaEM * aemFSR, AlphaEM * aemISR, vector<double>& asWeight,
     vector<double>& aemWeight, vector<double>& pdfWeight);
 
   // Function to return the \alpha_s-ratio part of the CKKWL weight.
   vector<double> weightTreeAlphaS( double as0, AlphaStrong * asFSR,
     AlphaStrong * asISR, int njetMax = -1, bool asVarInME = false );
   // Function to return the \alpha_em-ratio part of the CKKWL weight.
   vector<double> weightTreeAlphaEM( double aem0, AlphaEM * aemFSR,
     AlphaEM * aemISR, int njetMax = -1 );
   // Function to return the PDF-ratio part of the CKKWL weight.
   vector<double> weightTreePDFs( double maxscale, double pdfScale,
     int njetMax = -1 );
   // Function to return the no-emission probability part of the CKKWL weight.
   vector<double> weightTreeEmissions( PartonLevel* trial, int type,
     int njetMin, 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_em)-term of the \alpha_em-ratios
   // appearing in the CKKWL-weight.
   double weightFirstAlphaEM( double aem0, double muR, AlphaEM * aemFSR,
     AlphaEM * aemISR);
   // 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 )
   vector<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> findEWTripleW( int emtTagIn, const Event& event,
                        vector<int> posFinalPartn, vector<int> posInitPartn );
   vector<Clustering> findEWTripleZ( int emtTagIn, const Event& event,
                        vector<int> posFinalPartn, vector<int> posInitPartn );
 
   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 );
 
   // Function to attach (spin-dependent duplicates of) a clustering.
   void attachClusterings (vector<Clustering>& clus, int iEmt, int iRad,
     int iRec, int iPartner, double pT, const Event& event);
 
   // 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( 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 spin of the radiator before the splitting
   // IN int  : Spin of the radiator after the splitting
   //    int  : Spin of the emitted after the splitting
   // OUT int : Spin of the radiator before the splitting
   int getRadBeforeSpin(const int radAfter, const int emtAfter,
         const int spinRadAfter, const int spinEmtAfter,
         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,
     int idRadBef = 0) const;
 
   // Function to compute "pythia pT separation" from Particle input
   double pTLund(const Event& event, int radAfterBranch, int emtAfterBranch,
     int recAfterBranch, int showerType, int idRadBef = 0);
 
   // 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);
 
   // Function providing a list of possible new flavours after a w emssion
   // from the input flavour.
   vector<int> posFlavCKM(int flav);
 
   // Check if the new flavour structure is possible.
   // If clusType is 1 final clustering is assumed, otherwise initial clustering
   // is assumed.
   bool checkFlavour(vector<int>& flavCounts, int flavRad, int flavRadBef,
     int clusType);
 
   // Check if the weak recoil structure is allowed.
   bool checkWeakRecoils(map<int,int> &allowedRecoils, bool isFirst = false);
 
   // Find the recoiler for an ISR scattered weak particle.
 
   int findISRRecoiler();
 
   // Reverse the boost carried out by the ISR.
   // The three first momenta are given by the ME,
   // the last two are filled in by this function.
   void reverseBoostISR(Vec4& pMother, Vec4& pSister, Vec4& pPartner,
     Vec4& pDaughter, Vec4& pRecoiler, int sign, double eCM, double& phi);
 
   // Check if an event reclustered into a 2 -> 2 dijet.
   // (Only enabled if W reclustering is used).
   bool isQCD2to2(const Event& event);
 
   // Check if an event reclustered into a 2 -> 1 Drell-Yan.
   // (Only enabled if W reclustering is used).
   bool isEW2to1(const Event& event);
 
   // Set selected child indices.
   void setSelectedChild();
 
   // Setup the weak dipole showers.
   void setupSimpleWeakShower(int nSteps);
 
   // Update weak dipoles after an emission.
   void transferSimpleWeakShower(vector<int> &mode, vector<Vec4> &mom,
     vector<int> fermionLines, vector<pair<int,int> > &dipoles, int nSteps);
 
   // Update the weak modes after an emission.
   vector<int> updateWeakModes(vector<int>& weakModes,
     map<int,int>& stateTransfer);
 
   // Update the fermion lines for the 2 -> 2 process. This is needed for
   // the weak probabilities.
   vector<int> updateWeakFermionLines(vector<int> fermionLines,
     map<int,int>& stateTransfer);
 
   // Update the list of weak dipoles. This is needed to setup the PS correctly.
   vector<pair<int,int> > updateWeakDipoles(vector<pair<int,int> > &dipoles,
     map<int,int>& stateTransfer);
 
   // Setup the hard process information needed for calculating weak
   // probabilities and setting up the shower.
   void setupWeakHard(vector<int>& mode, vector<int>& fermionLines,
     vector<Vec4>& mom);
 
   // Functions to retrieve scale information from external showers.
   double getShowerPluginScale(const Event& event, int rad, int emt, int rec,
     string key, double scalePythia);
 
   //----------------------------------------------------------------------//
   // 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;
 
   // After a path is selected, store the child index.
   int selectedChild;
 
   // 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.
   MergingHooksPtr mergingHooksPtr;
 
    // The default constructor is private.
   History() : mother(), selectedChild(), sumpath(), sumGoodBranches(),
     sumBadBranches(), foundOrderedPath(), foundStronglyOrderedPath(),
     foundAllowedPath(), foundCompletePath(), scale(), nextInInput(), prob(),
     iReclusteredOld(), iReclusteredNew(), doInclude(), mergingHooksPtr(),
     particleDataPtr(), infoPtr(), loggerPtr(), showers(), coupSMPtr(),
     sumScalarPT(), probMaxSave(), depth(), minDepthSave() {}
 
   // 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;
 
   // Logger object.
   Logger* loggerPtr;
 
   // Class for calculation weak shower ME.
   SimpleWeakShowerMEs simpleWeakShowerMEs;
 
   // Pointer to showers, to simplify external clusterings.
   PartonLevel* showers;
 
   // Pointer to standard model couplings.
   CoupSM*     coupSMPtr;
 
   // Minimal scalar sum of pT used in Herwig to choose history
   double sumScalarPT;
 
   double probMaxSave;
   double probMax() {
     if (mother) return mother->probMax();
     return probMaxSave;
   }
   void updateProbMax(double probIn, bool isComplete = false) {
     if (mother) mother->updateProbMax(probIn, isComplete);
     if (!isComplete && !foundCompletePath) return;
     if (abs(probIn) > probMaxSave) probMaxSave = probIn;
   }
 
   int depth, minDepthSave;
   int minDepth() {
     if ( mother ) return mother->minDepth();
     return minDepthSave;
   }
   void updateMinDepth(int depthIn) {
     if ( mother ) return mother->updateMinDepth(depthIn);
     minDepthSave = (minDepthSave>0) ? min(minDepthSave,depthIn) : depthIn;
   }
 
 };
 
 //==========================================================================
 
 } // end namespace Pythia8
 
 #endif // end Pythia8_History_H

This page was automatically generated by the 2.3.7 LXR engine. The LXR team