PYTHIA
8.303

#include <History.h>
Public Member Functions  
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)  
~History ()  
The destructor deletes each child.  
bool  projectOntoDesiredHistories () 
Function to project paths onto desired paths. More...  
vector< double >  weightCKKWL (PartonLevel *trial, AlphaStrong *asFSR, AlphaStrong *asISR, AlphaEM *aemFSR, AlphaEM *aemISR, double RN) 
vector< double >  weightNL3Loop (PartonLevel *trial, double RN) 
vector< double >  weightNL3First (PartonLevel *trial, AlphaStrong *asFSR, AlphaStrong *asISR, AlphaEM *aemFSR, AlphaEM *aemISR, double RN, Rndm *rndmPtr) 
Return O()term of CKKWLweight for NL3 merging. More...  
vector< double >  weightNL3Tree (PartonLevel *trial, AlphaStrong *asFSR, AlphaStrong *asISR, AlphaEM *aemFSR, AlphaEM *aemISR, double RN) 
vector< double >  weightUMEPSTree (PartonLevel *trial, AlphaStrong *asFSR, AlphaStrong *asISR, AlphaEM *aemFSR, AlphaEM *aemISR, double RN) 
For UMEPS: More...  
vector< double >  weightUMEPSSubt (PartonLevel *trial, AlphaStrong *asFSR, AlphaStrong *asISR, AlphaEM *aemFSR, AlphaEM *aemISR, double RN) 
Function to return weight of virtual correction events for NLO merging. More...  
vector< double >  weightUNLOPSTree (PartonLevel *trial, AlphaStrong *asFSR, AlphaStrong *asISR, AlphaEM *aemFSR, AlphaEM *aemISR, double RN, int depthIn=1) 
For unitary NL3: More...  
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 calculate O()term of CKKWLweight for NLO merging. More...  
bool  foundAllowedHistories () 
Function to check if any allowed histories were found.  
bool  foundOrderedHistories () 
Function to check if any ordered histories were found.  
bool  foundCompleteHistories () 
Function to check if any ordered histories were found.  
void  getStartingConditions (const double RN, Event &outState) 
Function to set the state with complete scales for evolution. More...  
bool  getClusteredEvent (const double RN, int nSteps, Event &outState) 
Function to get the state with complete scales for evolution. More...  
bool  getFirstClusteredEventAboveTMS (const double RN, int nDesired, Event &process, int &nPerformed, bool updateProcess=true) 
Function to get the first reclustered state above the merging scale. More...  
int  nClusterings () 
Event  lowestMultProc (const double RN) 
Function to get the lowest multiplicity reclustered event. More...  
double  getPDFratio (int side, bool forSudakov, bool useHardPDF, int flavNum, double xNum, double muNum, int flavDen, double xDen, double muDen) 
Calculate and return pdf ratio. More...  
double  getWeakProb () 
Envelope function that calls the recursive getWeakProb. More...  
double  getWeakProb (vector< int > &mode, vector< Vec4 > &mom, vector< int > fermionLines) 
double  getSingleWeakProb (vector< int > &mode, vector< Vec4 > &mom, vector< int > fermionLines) 
void  findStateTransfer (map< int, int > &transfer) 
void  printHistory (const double RN) 
Function to print the history that would be chosen from the number RN. More...  
void  printStates () 
Function to print the states in a history, starting from the hard process. More...  
Friends  
class  Pythia 
Make Pythia class friend.  
class  Merging 
Make Merging class friend.  
Declaration of History class
A History object represents an event in a given step in the CKKWL clustering procedure. It defines a treelike 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 njets has been clustered down to the Bornlevel state. Also states which cannot be clustered down to the Bornlevel are possible  these will be called incomplete. The leaves are characterized by the vector of children being empty.
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 = true , 

bool  isStronglyOrdered = true , 

bool  isAllowed = true , 

bool  isNextInInput = true , 

double  probin = 1.0 , 

History *  mothin = 0 

) 
The only constructor. Default arguments are used when creating the initial history node. The depth is the maximum number of clusterings requested. scalein is the scale at which the statein was clustered (should be set to the merging scale for the initial history node. beamAIn and beamBIn are needed to calcutate PDF ratios, particleDataIn to have access to the correct masses of particles. If isOrdered is true, the previous clusterings has been ordered. is the PDF ratio for this clustering (=1 for FSR clusterings). probin is the accumulated probabilities for the previous clusterings, and \ mothin is the previous history node (null for the initial node).
Declaration of History class The only constructor. Default arguments are used when creating the initial history node. The depth is the maximum number of clusterings requested. scalein is the scale at which the statein was clustered (should be set to the merging scale for the initial history node. beamAIn and beamBIn are needed to calcutate PDF ratios, particleDataIn to have access to the correct masses of particles. If isOrdered is true, the previous clusterings has been ordered. is the PDF ratio for this clustering (=1 for FSR clusterings). probin is the accumulated probabilities for the previous clusterings, and \ mothin is the previous history node (null for the initial node).
Initialise beam particles
Update probability with PDF ratio
Minimal scalar sum of pT used in Herwig to choose history Keep track of scalar PT
Remember reclustered radiator in lower multiplicity state
Check if more steps should be taken.
Stop clustering at 2>1 massive. Stop clustering at 2>2 massless.
If this is not the fully clustered state, try to find possible QCD clusterings.
If necessary, try to find possible EW clusterings.
if ( depth > 0 && mergingHooksPtr>doWeakClustering() )
If necessary, try to find possible SQCD clusterings.
If no clusterings were found, the recursion is done and we register this node.
Multiply with hard process matrix element.
We'll now order the clusterings in such a way that an ordered history is found more rapidly. Following the branches with small pT is a good heuristic, as is following ISR clusterings.
This might be a marginally faster ordering. double z = getCurrentZ(clusterings[i].emittor, clusterings[i].recoiler, clusterings[i].emitted, clusterings[i].flavRadBef); double index = t/z; if (!state[clusterings[i].emittor].isFinal()) sort.insert(make_pair(1./index, &clusterings[i])); else sort.insert(make_pair(index, &clusterings[i]));
If this path is not strongly ordered and we already have found an ordered path, then we don't need to continue along this path.
Check if reclustering follows ordered sequence.
Get new z value
Get z value of splitting that produced this state
If this path is not ordered in pT and y, and we already have found an ordered path, then we don't need to continue along this path.
If this path is not ordered in pT and we already have found an ordered path, then we don't need to continue along this path, unless we have not yet found an allowed path.
Check if reclustered state should be disallowed.
Skip if this branch is already strongly suppressed.
Perform the clustering and recurse and construct the next history node.
void findStateTransfer  (  map< int, int > &  transfer  ) 
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).
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)
No need to transfer if already at highest multiplicity.
Directly assign the 3 first particles (system, beam1, beam2);
Handle all particles that are not part of the clustering.
bool getClusteredEvent  (  const double  RN, 
int  nSteps,  
Event &  outState  
) 
Function to get the state with complete scales for evolution.
Function to set the state with complete scales for evolution.
Select history
Set scales in the states to the scales pythia would have set (Only needed if not done before in calculation of weights or setting of starting conditions)
If the history does not allow for nSteps clusterings (e.g. because the history is incomplete), return false
Return event with nSteps1 additional partons (i.e. recluster the last splitting) and copy the output state
Done.
bool getFirstClusteredEventAboveTMS  (  const double  RN, 
int  nDesired,  
Event &  process,  
int &  nPerformed,  
bool  updateProcess = true 

) 
Function to get the first reclustered state above the merging scale.
Do reclustering (looping) steps. Remember process scale.
Get number of clustering steps.
Set scales in the states to the scales pythia would have set.
Recluster until reclustered event is above the merging scale.
Initialise temporary output of reclustering.
Recluster once more.
If reclustered event does not exist, exit.
Continue loop if reclustered event has unresolved partons.
Update the hard process.
Failed to produce output state.
Update to the actual number of steps.
Save MPI starting scale
Done
double getPDFratio  (  int  side, 
bool  forSudakov,  
bool  useHardPDF,  
int  flavNum,  
double  xNum,  
double  muNum,  
int  flavDen,  
double  xDen,  
double  muDen  
) 
Calculate and return pdf ratio.
Do nothing for e+e beams
Now calculate PDF ratio if necessary
Get mother and daughter pdfs
Use hard process PDFs (i.e. PDFs NOT used in ISR, FSR or MPI).
Use rescaled PDFs in the presence of multiparton interactions
Cut out charm threshold.
Return ratio of pdfs
Done
double getSingleWeakProb  (  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 tchannel, 4 is double quark uchannel.
Find the correct coupling coefficient.
No emissions from right handed particles.
No emissions from right handed particles.
Find and store kinematics (e.g. z, pT, k1, k3).
Store momenta.
Final state clustering.
schannel
Calculate variables.
Calculate Jacobian.
tchannel.
Store momentas needed.
Check if a swap is needed.
Rescaling of incoming partons p3 and p4.
Longitudinal boost to rest frame of incoming partons of hard interaction.
Further boost to rest frame of outgoing state.
Calculate variables;
Calculate the ME depending on the top of process.
Split matrix element according to propagaters.
Initial clustering.
schannel
Store momenta.
Check if radiator is from beam one or two.
Undo the ISR boost.
Scale outgoing vectors to conserve energy / momentum. double scaleFactor2 = (pIn1 + pIn2  p3).m2Calc() / (p1 + p2).m2Calc();
Find 2 to 2 rest frame for incoming particles. This is done before one of the two are made virtual (Q^2 mass).
Calculating the Jacobian
Split of ME into an ISR part and FSR part.
void getStartingConditions  (  const double  RN, 
Event &  outState  
) 
Function to set the state with complete scales for evolution.
Select the history
Set scales in the states to the scales pythia would have set
Get number of clustering steps.
Update the lowest order process.
Save information on last splitting, to allow the next emission in the shower to have smaller rapidity with respect to the last ME splitting. For hard process, use dummy values.
For incomplete process, try to use real values.
Set QCD 2>2 starting scale different from arbitrary scale in LHEF! –> Set to minimal mT of partons.
For pure QCD dijet events (only!), set the process scale to the transverse momentum of the outgoing partons.
For weak inclusive merging, follow QCD 2>2 starting scale for dijet events. Also, restore input input polarisations.
Save information on last splitting, to allow the next emission in the shower to have smaller rapidity with respect to the last ME splitting.
Copy the output state.
Save MPI starting scale.
Setup the weak shower if W clustering is enabled.
double getWeakProb  (  ) 
Envelope function that calls the recursive getWeakProb.
Setup function that call the real getWeakProb.
double getWeakProb  (  vector< int > &  mode, 
vector< Vec4 > &  mom,  
vector< int >  fermionLines  
) 
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 tchannel, 4 is double quark uchannel.
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 tchannel, 4 is double quark uchannel.
If at end, return 1.
Find the transfer map given the splitting.
Setup hard process.
Update modes and fermionLines.
Get the probability if it is a weak emission.

inline 
Function to get the lowest multiplicity reclustered event.
Return lowest multiplicity state
int nClusterings  (  ) 
Function to return the depth of the history (i.e. the number of reclustered splittings)
Function to return the depth of the history (i.e. the number of reclustered splittings) NO INPUT OUTPUT int : Depth of history
void printHistory  (  const double  RN  ) 
Function to print the history that would be chosen from the number RN.
Function to print the history that would be chosen from the random number RN. Mainly for debugging.
Done
void printStates  (  ) 
Function to print the states in a history, starting from the hard process.
Function to print the states in a history, starting from the hard process. Mainly for debugging.
Print.
Recurse
Done
bool projectOntoDesiredHistories  (  ) 
Function to project paths onto desired paths.
Function to project all possible paths onto only the desired paths.
At the moment, only trim histories.
vector< double > weightCKKWL  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN  
) 
For CKKWL, 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)
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)
Read alpha_S in ME calculation and maximal scale (eCM)
Select a path of clusterings
Set scales in the states to the scales pythia would have set
Get weight.
Do trial shower, calculation of alpha_S ratios, PDF ratios
MPI noemission probability
Set hard process renormalisation scale to default Pythia value.
For pure QCD dijet events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. Here we choose FSR for simplicity.
Reset to a running coupling. Here we choose FSR for simplicity.
For W clustering, correct the .
Reset to a running coupling. Here we choose FSR for simplicity.
For prompt photon events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. In prompt photon always ISR.
Done
vector< double > weightNL3First  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN,  
Rndm *  rndmPtr  
) 
Return O()term of CKKWLweight for NL3 merging.
Function to calculate O()term of CKKWLweight for NLO merging.
Read alpha_S in ME calculation and maximal scale (eCM)
Pick path of clusterings
Set scales in the states to the scales pythia would have set
Get the lowest order kfactor and add first two terms in expansion
If using Bbar, which includes a treelevel part, subtract an additional one, i.e. the O(^0) contribution as well
Calculate sum of O(alpha) terms
Get starting scale for trial showers.
Count emissions: New variant Generate true average, not only onepoint
Get number of emissions
Introduce vector to allow variation of coefficient
Use the varied scale in the coefficient around which we expand
Introduce variation of stong coupling that is not done in Born input
Done
vector< double > weightNL3Loop  (  PartonLevel *  trial, 
double  RN  
) 
For default NL3: Return weight of virtual correction and subtractive for NL3 merging
Function to return weight of virtual correction and subtractive events for NL3 merging
Select a path of clusterings
Set scales in the states to the scales pythia would have set
So far, no reweighting
Only reweighting with MPI noemission probability
Done
vector< double > weightNL3Tree  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN  
) 
No difference to CKKWL. Recycle CKKWL function.
vector< double > weightUMEPSSubt  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN  
) 
Function to return weight of virtual correction events for NLO merging.
Read alpha_S in ME calculation and maximal scale (eCM)
Select a path of clusterings
Set scales in the states to the scales pythia would have set
Get weight.
Do trial shower, calculation of alpha_S ratios, PDF ratios
MPI noemission probability.
Set hard process renormalisation scale to default Pythia value.
For pure QCD dijet events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. Here we choose FSR for simplicity.
For prompt photon events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. In prompt photon always ISR.
Done
vector< double > weightUMEPSTree  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN  
) 
For UMEPS:
No difference to CKKWL. Recycle CKKWL function.
vector< double > weightUNLOPSFirst  (  int  order, 
PartonLevel *  trial,  
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN,  
Rndm *  rndmPtr  
) 
Function to calculate O()term of CKKWLweight for NLO merging.
Already done if no correction should be calculated
Read alpha_S in ME calculation and maximal scale (eCM)
double aemME = infoPtr>alphaEM();
Pick path of clusterings
Set scales in the states to the scales pythia would have set
Get the lowest order kfactor and add first two terms in expansion
If using Bbar, which includes a treelevel part, subtract an additional one, i.e. the O(^0) contribution as well
Set up vector for order == 0 case.
Start by adding the O(^1)term of the kfactor.
Calculate sum of O(^1)terms of the ckkwl weight WITHOUT the O(^1)term of the last noemission probability. Get first term in expansion of alpha_s ratios.
Generate true average, not only onepoint.
Get average number of emissions.
Add average number of emissions off reconstructed states to weight.
Get first term in expansion of PDF ratios.
Add integral of DGLAP shifted PDF ratios from expansion to wt.
Introduce vector to allow variation of coefficient
Use the varied scale in the coefficient around which we expand
Introduce variation of stong coupling that is not done in Born input
If O(^1)term + O(^1)term is to be calculated, done.
So far, no calculation of O(^2)term
vector< double > weightUNLOPSLoop  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN,  
int  depthIn = 1 

) 
No difference to default NL3
Read alpha_S in ME calculation and maximal scale (eCM)
Select a path of clusterings
Set scales in the status to the scales pythia would have set
Get weight.
Do trial shower, calculation of alpha_S ratios, PDF ratios
MPI noemission probability.
Set hard process renormalisation scale to default Pythia value.
For pure QCD dijet events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. Here we choose FSR for simplicity.
For prompt photon events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. In prompt photon always ISR.
Done
Save weight vectors interally for UNLOPSP and PC
vector< double > weightUNLOPSSubt  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN,  
int  depthIn = 1 

) 
Select a path of clusterings
Set scales in the states to the scales pythia would have set
Read alpha_S in ME calculation and maximal scale (eCM)
Only allow two clusterings if all intermediate states above the merging scale.
Get weights: alpha_S ratios and PDF ratios
Do trial shower, calculation of alpha_S ratios, PDF ratios
MPI noemission probability.
Set weight
For tree level, undo as variation applied to ME component to avoid double ratios when combining later.
Save weight vectors interally for UNLOPSP and PC
Done
vector< double > weightUNLOPSSubtNLO  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN,  
int  depthIn = 1 

) 
So far, no reweighting
Select a path of clusterings
Set scales in the states to the scales pythia would have set
Only reweighting with MPI noemission probability
Done
Select a path of clusterings
Set scales in the states to the scales pythia would have set
Read alpha_S in ME calculation and maximal scale (eCM)
Only allow two clusterings if all intermediate states above the merging scale.
Get weights: alpha_S ratios and PDF ratios
Do trial shower, calculation of alpha_S ratios, PDF ratios
MPI noemission probability.
Set weight
Save weight vectors interally for UNLOPSP and PC
Done
vector< double > weightUNLOPSTree  (  PartonLevel *  trial, 
AlphaStrong *  asFSR,  
AlphaStrong *  asISR,  
AlphaEM *  aemFSR,  
AlphaEM *  aemISR,  
double  RN,  
int  depthIn = 1 

) 
For unitary NL3:
Read alpha_S in ME calculation and maximal scale (eCM)
Select a path of clusterings
Set scales in the states to the scales pythia would have set
Get weight.
Do trial shower, calculation of alpha_S ratios, PDF ratios
MPI noemission probability.
Set hard process renormalisation scale to default Pythia value.
For pure QCD dijet events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. Here we choose FSR for simplicity.
For prompt photon events, evaluate the coupling of the hard process at a more reasonable pT, rather than evaluation at a fixed arbitrary scale.
Reset to a running coupling. In prompt photon always ISR.
Done
For tree level, undo as variation applied to ME component to avoid double ratios when combining later.
Save weight vectors internally for UNLOPSP and PC