PYTHIA  8.312
History.h
1 // History.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // This file is written by Stefan Prestel.
7 // It contains the main class for matrix element merging.
8 // Header file for the Clustering and History classes.
9 
10 #ifndef Pythia8_History_H
11 #define Pythia8_History_H
12 
13 #include "Pythia8/Basics.h"
14 #include "Pythia8/BeamParticle.h"
15 #include "Pythia8/Event.h"
16 #include "Pythia8/Info.h"
17 #include "Pythia8/ParticleData.h"
18 #include "Pythia8/PartonLevel.h"
19 #include "Pythia8/PythiaStdlib.h"
20 #include "Pythia8/Settings.h"
21 #include "Pythia8/StandardModel.h"
22 #include "Pythia8/SimpleWeakShowerMEs.h"
23 
24 namespace Pythia8 {
25 
26 //==========================================================================
27 
28 // Declaration of Clustering class.
29 // This class holds information about one radiator, recoiler, emitted system.
30 // This class is a container class for History class use.
31 
32 class Clustering {
33 
34 public:
35 
36  // The emitted parton location.
37  int emitted;
38  // The emittor parton
39  int emittor;
40  // The recoiler parton.
41  int recoiler;
42  // The colour connected recoiler (Can be different for ISR)
43  int partner;
44  // The scale associated with this clustering.
45  double pTscale;
46  // The flavour of the radiator prior to the emission.
48  // Spin of the radiator (-1 is left handed, +1 is right handed).
49  int spinRad;
50  // Spin of the emitted (-1 is left handed, +1 is right handed).
51  int spinEmt;
52  // Spin of the recoiler (-1 is left handed, +1 is right handed).
53  int spinRec;
54  // Spin of the radiator before the splitting.
56  // The radiator before the splitting.
57  int radBef;
58  // The recoiler before the splitting.
59  int recBef;
60 
61  // Map between particle positions in the clustered states -> particle
62  // positions in unclustered real-emission state.
63  map<int,int> iPosInMother;
64 
65  // Default constructor
66  Clustering() : emitted(0), emittor(0), recoiler(0), partner(0), pTscale(),
67  flavRadBef(0), spinRad(9), spinEmt(9), spinRec(9), spinRadBef(9),
68  radBef(0), recBef(0), iPosInMother() {}
69 
70  // Default destructor
72 
73  // Copy constructor
74  Clustering( const Clustering& inSystem ) :
75  emitted(inSystem.emitted), emittor(inSystem.emittor),
76  recoiler(inSystem.recoiler), partner(inSystem.partner),
77  pTscale(inSystem.pTscale), flavRadBef(inSystem.flavRadBef),
78  spinRad(inSystem.spinRad), spinEmt(inSystem.spinEmt),
79  spinRec(inSystem.spinRec), spinRadBef(inSystem.spinRad),
80  radBef(inSystem.radBef), recBef(inSystem.recBef),
81  iPosInMother(inSystem.iPosInMother) {}
82 
83  // Constructor with input
84  Clustering( int emtIn, int radIn, int recIn, int partnerIn,
85  double pTscaleIn, int flavRadBefIn = 0, int spinRadIn = 9,
86  int spinEmtIn = 9, int spinRecIn = 9, int spinRadBefIn = 9,
87  int radBefIn = 0, int recBefIn = 0, map<int,int> posIn = map<int,int>())
88  : emitted(emtIn), emittor(radIn), recoiler(recIn), partner(partnerIn),
89  pTscale(pTscaleIn), flavRadBef(flavRadBefIn), spinRad(spinRadIn),
90  spinEmt(spinEmtIn), spinRec(spinRecIn), spinRadBef(spinRadBefIn),
91  radBef(radBefIn), recBef(recBefIn), iPosInMother(posIn) {}
92 
93  // Function to return pythia pT scale of clustering
94  double pT() const { return pTscale; }
95 
96  // print for debug
97  void list() const;
98 
99 };
100 
101 //==========================================================================
102 
103 // Declaration of History class
104 //
105 // A History object represents an event in a given step in the CKKW-L
106 // clustering procedure. It defines a tree-like recursive structure,
107 // where the root node represents the state with n jets as given by
108 // the matrix element generator, and is characterized by the member
109 // variable mother being null. The leaves on the tree corresponds to a
110 // fully clustered paths where the original n-jets has been clustered
111 // down to the Born-level state. Also states which cannot be clustered
112 // down to the Born-level are possible - these will be called
113 // incomplete. The leaves are characterized by the vector of children
114 // being empty.
115 
116 class History {
117 
118 public:
119 
120  // The only constructor. Default arguments are used when creating
121  // the initial history node. The \a depth is the maximum number of
122  // clusterings requested. \a scalein is the scale at which the \a
123  // statein was clustered (should be set to the merging scale for the
124  // initial history node. \a beamAIn and beamBIn are needed to
125  // calcutate PDF ratios, \a particleDataIn to have access to the
126  // correct masses of particles. If \a isOrdered is true, the previous
127  // clusterings has been ordered. \a is the PDF ratio for this
128  // clustering (=1 for FSR clusterings). \a probin is the accumulated
129  // probabilities for the previous clusterings, and \ mothin is the
130  // previous history node (null for the initial node).
131  History( int depthIn,
132  double scalein,
133  Event statein,
134  Clustering c,
135  MergingHooksPtr mergingHooksPtrIn,
136  BeamParticle beamAIn,
137  BeamParticle beamBIn,
138  ParticleData* particleDataPtrIn,
139  Info* infoPtrIn,
140  PartonLevel* showersIn,
141  CoupSM* coupSMPtrIn,
142  bool isOrdered,
143  bool isStronglyOrdered,
144  bool isAllowed,
145  bool isNextInInput,
146  double probin,
147  History * mothin);
148 
149  // The destructor deletes each child.
151  for ( int i = 0, N = children.size(); i < N; ++i ) delete children[i];
152  }
153 
154  // Function to project paths onto desired paths.
155  bool projectOntoDesiredHistories();
156 
157  // For CKKW-L, NL3 and UMEPS:
158  // In the initial history node, select one of the paths according to
159  // the probabilities. This function should be called for the initial
160  // history node.
161  // IN trialShower* : Previously initialised trialShower object,
162  // to perform trial showering and as
163  // repository of pointers to initialise alphaS
164  // PartonSystems* : PartonSystems object needed to initialise
165  // shower objects
166  // OUT vector<double> : (Sukadov) , (alpha_S ratios) , (PDF ratios)
167  vector<double> weightCKKWL(PartonLevel* trial, AlphaStrong * asFSR,
168  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
169 
170 
171  // For default NL3:
172  // Return weight of virtual correction and subtractive for NL3 merging
173  vector<double> weightNL3Loop(PartonLevel* trial, double RN);
174  // Return O(\alpha_s)-term of CKKWL-weight for NL3 merging
175  vector<double> weightNL3First(PartonLevel* trial, AlphaStrong* asFSR,
176  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
177  Rndm* rndmPtr);
178  vector<double> weightNL3Tree(PartonLevel* trial, AlphaStrong * asFSR,
179  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
180 
181 
182  // For UMEPS:
183  vector<double> weightUMEPSTree(PartonLevel* trial, AlphaStrong * asFSR,
184  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
185  vector<double> weightUMEPSSubt(PartonLevel* trial, AlphaStrong * asFSR,
186  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN);
187 
188 
189  // For unitary NL3:
190  vector<double> weightUNLOPSTree(PartonLevel* trial, AlphaStrong * asFSR,
191  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
192  int depthIn = -1);
193  vector<double> weightUNLOPSSubt(PartonLevel* trial, AlphaStrong * asFSR,
194  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
195  int depthIn = -1);
196  vector<double> weightUNLOPSLoop(PartonLevel* trial, AlphaStrong * asFSR,
197  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
198  int depthIn = -1);
199  vector<double> weightUNLOPSSubtNLO(PartonLevel* trial, AlphaStrong * asFSR,
200  AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR, double RN,
201  int depthIn = -1);
202  vector<double> weightUNLOPSFirst( int order, PartonLevel* trial,
203  AlphaStrong* asFSR, AlphaStrong * asISR, AlphaEM * aemFSR,
204  AlphaEM * aemISR, double RN, Rndm* rndmPtr );
205 
206  // Function to check if any allowed histories were found
208  return (children.size() > 0 && foundAllowedPath); }
209  // Function to check if any ordered histories were found
211  return (children.size() > 0 && foundOrderedPath); }
212  // Function to check if any ordered histories were found
214  return (children.size() > 0 && foundCompletePath); }
215 
216  // Function to set the state with complete scales for evolution
217  void getStartingConditions( const double RN, Event& outState );
218  // Function to get the state with complete scales for evolution
219  bool getClusteredEvent( const double RN, int nSteps, Event& outState);
220  // Function to get the first reclustered state above the merging scale.
221  bool getFirstClusteredEventAboveTMS( const double RN, int nDesired,
222  Event& process, int & nPerformed, bool updateProcess = true );
223  // Function to return the depth of the history (i.e. the number of
224  // reclustered splittings)
225  int nClusterings();
226 
227  // Function to get the lowest multiplicity reclustered event
228  Event lowestMultProc( const double RN) {
229  // Return lowest multiplicity state
230  return (select(RN)->state);
231  }
232 
233  // Calculate and return pdf ratio
234  double getPDFratio( int side, bool forSudakov, bool useHardPDF,
235  int flavNum, double xNum, double muNum,
236  int flavDen, double xDen, double muDen);
237 
238  // Envelope function that calls the recursive getWeakProb.
239  double getWeakProb();
240 
241  // Recursive function that returns the weak probability for the given path.
242  // Mode refers to which ME correction to use, 1 = sChannel, 2 = gluon
243  // channel, 3 = double quark t-channel, 4 is double quark u-channel.
244  double getWeakProb(vector<int>& mode, vector<Vec4>& mom,
245  vector<int> fermionLines);
246 
247  // return the weak probability of a single reclustering.
248  // Mode refers to which ME correction to use, 1 = sChannel, 2 = gluon
249  // channel, 3 = double quark t-channel, 4 is double quark u-channel.
250  double getSingleWeakProb(vector<int> &mode, vector<Vec4> &mom,
251  vector<int> fermionLines);
252 
253  // Find map between indecies in the current state and the state after
254  // the splitting.
255  // NOT IMPLEMENTED FOR MULTIPLE W/Z/GAMMA (NEED TO HAVE A WAY TO
256  // IDENTIFY THEM).
257  void findStateTransfer(map<int,int> &transfer);
258 
259  // Function to print the history that would be chosen from the random number
260  // RN. Mainly for debugging.
261  void printHistory( const double RN );
262  // Function to print the states in a history, starting from the hard process.
263  // Mainly for debugging.
264  void printStates();
265 
266  // Make Pythia class friend
267  friend class Pythia;
268  // Make Merging class friend
269  friend class Merging;
270 
271 private:
272 
273  // Number of trial emission to use for calculating the average number of
274  // emissions
275  static const int NTRIAL;
276 
277  // Function to set all scales in the sequence of states. This is a
278  // wrapper routine for setScales and setEventScales methods
279  void setScalesInHistory();
280 
281  // Function to find the index (in the mother histories) of the
282  // child history, thus providing a way access the path from both
283  // initial history (mother == 0) and final history (all children == 0)
284  // IN vector<int> : The index of each child in the children vector
285  // of the current history node will be saved in
286  // this vector
287  // NO OUTPUT
288  void findPath(vector<int>& out);
289 
290  // Functions to set the parton production scales and enforce
291  // ordering on the scales of the respective clusterings stored in
292  // the History node:
293  // Method will start from lowest multiplicity state and move to
294  // higher states, setting the production scales the shower would
295  // have used.
296  // When arriving at the highest multiplicity, the method will switch
297  // and go back in direction of lower states to check and enforce
298  // ordering for unordered histories.
299  // IN vector<int> : Vector of positions of the chosen child
300  // in the mother history to allow to move
301  // in direction initial->final along path
302  // bool : True: Move in direction low->high
303  // multiplicity and set production scales
304  // False: Move in direction high->low
305  // multiplicity and check and enforce
306  // ordering
307  // NO OUTPUT
308  void setScales( vector<int> index, bool forward);
309 
310  // Function to find a particle in all higher multiplicity events
311  // along the history path and set its production scale to the input
312  // scale
313  // IN int iPart : Parton in refEvent to be checked / rescaled
314  // Event& refEvent : Reference event for iPart
315  // double scale : Scale to be set as production scale for
316  // unchanged copies of iPart in subsequent steps
317  void scaleCopies(int iPart, const Event& refEvent, double rho);
318 
319  // Function to set the OVERALL EVENT SCALES [=state.scale()] to
320  // the scale of the last clustering
321  // NO INPUT
322  // NO OUTPUT
323  void setEventScales();
324 
325  // Function to print information on the reconstructed scales in one path.
326  // For debug only.
327  void printScales() { if ( mother ) mother->printScales();
328  cout << " size " << state.size() << " scale " << scale << " clusterIn "
329  << clusterIn.pT() << " state.scale() " << state.scale() << endl; }
330 
331  // Function to project paths onto desired paths.
332  bool trimHistories();
333  // Function to tag history for removal.
334  void remove(){ doInclude = false; }
335  // Function to return flag of allowed histories to choose from.
336  bool keep(){ return doInclude; }
337  // Function implementing checks on a paths, for deciding if the path should
338  // be considered valid or not.
339  bool keepHistory();
340  // Function to check if a path is ordered in evolution pT.
341  bool isOrderedPath( double maxscale );
342 
343  bool followsInputPath();
344 
345  // Function to check if all reconstucted states in a path pass the merging
346  // scale cut.
347  bool allIntermediateAboveRhoMS( double rhoms, bool good = true );
348  // Function to check if any ordered paths were found (and kept).
349  bool foundAnyOrderedPaths();
350 
351  // Functions to return the z value of the last ISR splitting
352  // NO INPUT
353  // OUTPUT double : z value of last ISR splitting in history
354  double zISR();
355 
356  // Functions to return the z value of the last FSR splitting
357  // NO INPUT
358  // OUTPUT double : z value of last FSR splitting in history
359  double zFSR();
360 
361  // Functions to return the pT scale of the last ISR splitting
362  // NO INPUT
363  // OUTPUT double : pT scale of last ISR splitting in history
364  double pTISR();
365 
366  // Functions to return the pT scale of the last FSR splitting
367  // NO INPUT
368  // OUTPUT double : pT scale of last FSR splitting in history
369  double pTFSR();
370 
371  // Functions to return the event with nSteps additional partons
372  // INPUT int : Number of splittings in the event,
373  // as counted from core 2->2 process
374  // OUTPUT Event : event with nSteps additional partons
375  Event clusteredState( int nSteps);
376 
377  // Function to choose a path from all paths in the tree
378  // according to their splitting probabilities
379  // IN double : Random number
380  // OUT History* : Leaf of history path chosen
381  History * select(double rnd);
382 
383  // For a full path, find the weight calculated from the ratio of
384  // couplings, the no-emission probabilities, and possible PDF
385  // ratios. This function should only be called for the last history
386  // node of a full path.
387  // IN TimeShower : Already initialised shower object to be used as
388  // trial shower
389  // double : alpha_s value used in ME calculation
390  // double : Maximal mass scale of the problem (e.g. E_CM)
391  // AlphaStrong: Initialised shower alpha_s object for FSR alpha_s
392  // ratio calculation
393  // AlphaStrong: Initialised shower alpha_s object for ISR alpha_s
394  // ratio calculation (can be different from previous)
395  vector<double> weightTree(PartonLevel* trial, double as0, double aem0,
396  double maxscale, double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,
397  AlphaEM * aemFSR, AlphaEM * aemISR, vector<double>& asWeight,
398  vector<double>& aemWeight, vector<double>& pdfWeight);
399 
400  // Function to return the \alpha_s-ratio part of the CKKWL weight.
401  vector<double> weightTreeAlphaS( double as0, AlphaStrong * asFSR,
402  AlphaStrong * asISR, int njetMax = -1, bool asVarInME = false );
403  // Function to return the \alpha_em-ratio part of the CKKWL weight.
404  vector<double> weightTreeAlphaEM( double aem0, AlphaEM * aemFSR,
405  AlphaEM * aemISR, int njetMax = -1 );
406  // Function to return the PDF-ratio part of the CKKWL weight.
407  vector<double> weightTreePDFs( double maxscale, double pdfScale,
408  int njetMax = -1 );
409  // Function to return the no-emission probability part of the CKKWL weight.
410  vector<double> weightTreeEmissions( PartonLevel* trial, int type,
411  int njetMin, int njetMax, double maxscale );
412 
413  // Function to generate the O(\alpha_s)-term of the CKKWL-weight
414  double weightFirst(PartonLevel* trial, double as0, double muR,
415  double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR, Rndm* rndmPtr );
416 
417  // Function to generate the O(\alpha_s)-term of the \alpha_s-ratios
418  // appearing in the CKKWL-weight.
419  double weightFirstAlphaS( double as0, double muR, AlphaStrong * asFSR,
420  AlphaStrong * asISR);
421  // Function to generate the O(\alpha_em)-term of the \alpha_em-ratios
422  // appearing in the CKKWL-weight.
423  double weightFirstAlphaEM( double aem0, double muR, AlphaEM * aemFSR,
424  AlphaEM * aemISR);
425  // Function to generate the O(\alpha_s)-term of the PDF-ratios
426  // appearing in the CKKWL-weight.
427  double weightFirstPDFs( double as0, double maxscale, double pdfScale,
428  Rndm* rndmPtr );
429  // Function to generate the O(\alpha_s)-term of the no-emission
430  // probabilities appearing in the CKKWL-weight.
431  double weightFirstEmissions(PartonLevel* trial, double as0, double maxscale,
432  AlphaStrong * asFSR, AlphaStrong * asISR, bool fixpdf, bool fixas );
433 
434  // Function to return the default factorisation scale of the hard process.
435  double hardFacScale(const Event& event);
436  // Function to return the default renormalisation scale of the hard process.
437  double hardRenScale(const Event& event);
438 
439  // Perform a trial shower using the \a pythia object between
440  // maxscale down to this scale and return the corresponding Sudakov
441  // form factor.
442  // IN trialShower : Shower object used as trial shower
443  // double : Maximum scale for trial shower branching
444  // OUT 0.0 : trial shower emission outside allowed pT range
445  // 1.0 : trial shower successful (any emission was below
446  // the minimal scale )
447  vector<double> doTrialShower(PartonLevel* trial, int type, double maxscale,
448  double minscale = 0.);
449 
450  // Function to bookkeep the indices of weights generated in countEmissions
451  bool updateind(vector<int> & ind, int i, int N);
452 
453  // Function to count number of emissions between two scales for NLO merging
454  vector<double> countEmissions(PartonLevel* trial, double maxscale,
455  double minscale, int showerType, double as0, AlphaStrong * asFSR,
456  AlphaStrong * asISR, int N, bool fixpdf, bool fixas);
457 
458  // Function to integrate PDF ratios between two scales over x and t,
459  // where the PDFs are always evaluated at the lower t-integration limit
460  double monteCarloPDFratios(int flav, double x, double maxScale,
461  double minScale, double pdfScale, double asME, Rndm* rndmPtr);
462 
463  // Default: Check if a ordered (and complete) path has been found in
464  // the initial node, in which case we will no longer be interested in
465  // any unordered paths.
466  bool onlyOrderedPaths();
467 
468  // Check if a strongly ordered (and complete) path has been found in the
469  // initial node, in which case we will no longer be interested in
470  // any unordered paths.
471  bool onlyStronglyOrderedPaths();
472 
473  // Check if an allowed (according to user-criterion) path has been found in
474  // the initial node, in which case we will no longer be interested in
475  // any forbidden paths.
476  bool onlyAllowedPaths();
477 
478  // When a full path has been found, register it with the initial
479  // history node.
480  // IN History : History to be registered as path
481  // bool : Specifying if clusterings so far were ordered
482  // bool : Specifying if path is complete down to 2->2 process
483  // OUT true if History object forms a plausible path (eg prob>0 ...)
484  bool registerPath(History & l, bool isOrdered, bool isStronglyOrdered,
485  bool isAllowed, bool isComplete);
486 
487  // For the history-defining state (and if necessary interfering
488  // states), find all possible clusterings.
489  // NO INPUT
490  // OUT vector of all (rad,rec,emt) systems
491  vector<Clustering> getAllQCDClusterings();
492 
493  // For one given state, find all possible clusterings.
494  // IN Event : state to be investigated
495  // OUT vector of all (rad,rec,emt) systems in the state
496  vector<Clustering> getQCDClusterings( const Event& event);
497 
498  // Function to construct (rad,rec,emt) triples from the event
499  // IN int : Position of Emitted in event record for which
500  // dipoles should be constructed
501  // int : Colour topogy to be tested
502  // 1= g -> qqbar, causing 2 -> 2 dipole splitting
503  // 2= q(bar) -> q(bar) g && g -> gg,
504  // causing a 2 -> 3 dipole splitting
505  // Event : event record to be checked for ptential partners
506  // OUT vector of all allowed radiator+recoiler+emitted triples
507  vector<Clustering> findQCDTriple (int emtTagIn, int colTopIn,
508  const Event& event, vector<int> posFinalPartn,
509  vector <int> posInitPartn );
510 
511  vector<Clustering> getAllEWClusterings();
512  vector<Clustering> getEWClusterings( const Event& event);
513  vector<Clustering> findEWTripleW( int emtTagIn, const Event& event,
514  vector<int> posFinalPartn, vector<int> posInitPartn );
515  vector<Clustering> findEWTripleZ( int emtTagIn, const Event& event,
516  vector<int> posFinalPartn, vector<int> posInitPartn );
517 
518  vector<Clustering> getAllSQCDClusterings();
519  vector<Clustering> getSQCDClusterings( const Event& event);
520  vector<Clustering> findSQCDTriple (int emtTagIn, int colTopIn,
521  const Event& event, vector<int> posFinalPartn,
522  vector <int> posInitPartn );
523 
524  // Function to attach (spin-dependent duplicates of) a clustering.
525  void attachClusterings (vector<Clustering>& clus, int iEmt, int iRad,
526  int iRec, int iPartner, double pT, const Event& event);
527 
528  // Calculate and return the probability of a clustering.
529  // IN Clustering : rad,rec,emt - System for which the splitting
530  // probability should be calcuated
531  // OUT splitting probability
532  double getProb(const Clustering & SystemIn);
533 
534  // Set up the beams (fill the beam particles with the correct
535  // current incoming particles) to allow calculation of splitting
536  // probability.
537  // For interleaved evolution, set assignments dividing PDFs into
538  // sea and valence content. This assignment is, until a history path
539  // is chosen and a first trial shower performed, not fully correct
540  // (since content is chosen form too high x and too low scale). The
541  // assignment used for reweighting will be corrected after trial
542  // showering
543  void setupBeams();
544 
545  // Calculate the PDF ratio used in the argument of the no-emission
546  // probability.
547  double pdfForSudakov();
548 
549  // Calculate the hard process matrix element to include in the selection
550  // probability.
551  double hardProcessME( const Event& event);
552 
553  // Perform the clustering of the current state and return the
554  // clustered state.
555  // IN Clustering : rad,rec,emt triple to be clustered to two partons
556  // OUT clustered state
557  Event cluster( Clustering & inSystem);
558 
559  // Function to get the flavour of the radiator before the splitting
560  // for clustering
561  // IN int : Position of the radiator after the splitting, in the event
562  // int : Position of the emitted after the splitting, in the event
563  // Event : Reference event
564  // OUT int : Flavour of the radiator before the splitting
565  int getRadBeforeFlav(const int radAfter, const int emtAfter,
566  const Event& event);
567 
568  // Function to get the spin of the radiator before the splitting
569  // IN int : Spin of the radiator after the splitting
570  // int : Spin of the emitted after the splitting
571  // OUT int : Spin of the radiator before the splitting
572  int getRadBeforeSpin(const int radAfter, const int emtAfter,
573  const int spinRadAfter, const int spinEmtAfter,
574  const Event& event);
575 
576  // Function to get the colour of the radiator before the splitting
577  // for clustering
578  // IN int : Position of the radiator after the splitting, in the event
579  // int : Position of the emitted after the splitting, in the event
580  // Event : Reference event
581  // OUT int : Colour of the radiator before the splitting
582  int getRadBeforeCol(const int radAfter, const int emtAfter,
583  const Event& event);
584 
585  // Function to get the anticolour of the radiator before the splitting
586  // for clustering
587  // IN int : Position of the radiator after the splitting, in the event
588  // int : Position of the emitted after the splitting, in the event
589  // Event : Reference event
590  // OUT int : Anticolour of the radiator before the splitting
591  int getRadBeforeAcol(const int radAfter, const int emtAfter,
592  const Event& event);
593 
594  // Function to get the parton connected to in by a colour line
595  // IN int : Position of parton for which partner should be found
596  // Event : Reference event
597  // OUT int : If a colour line connects the "in" parton with another
598  // parton, return the Position of the partner, else return 0
599  int getColPartner(const int in, const Event& event);
600 
601  // Function to get the parton connected to in by an anticolour line
602  // IN int : Position of parton for which partner should be found
603  // Event : Reference event
604  // OUT int : If an anticolour line connects the "in" parton with another
605  // parton, return the Position of the partner, else return 0
606  int getAcolPartner(const int in, const Event& event);
607 
608  // Function to get the list of partons connected to the particle
609  // formed by reclusterinf emt and rad by colour and anticolour lines
610  // IN int : Position of radiator in the clustering
611  // IN int : Position of emitted in the clustering
612  // Event : Reference event
613  // OUT vector<int> : List of positions of all partons that are connected
614  // to the parton that will be formed
615  // by clustering emt and rad.
616  vector<int> getReclusteredPartners(const int rad, const int emt,
617  const Event& event);
618 
619  // Function to extract a chain of colour-connected partons in
620  // the event
621  // IN int : Type of parton from which to start extracting a
622  // parton chain. If the starting point is a quark
623  // i.e. flavType = 1, a chain of partons that are
624  // consecutively connected by colour-lines will be
625  // extracted. If the starting point is an antiquark
626  // i.e. flavType =-1, a chain of partons that are
627  // consecutively connected by anticolour-lines
628  // will be extracted.
629  // IN int : Position of the parton from which a
630  // colour-connected chain should be derived
631  // IN Event : Refernence event
632  // IN/OUT vector<int> : Partons that should be excluded from the search.
633  // OUT vector<int> : Positions of partons along the chain
634  // OUT bool : Found singlet / did not find singlet
635  bool getColSinglet(const int flavType, const int iParton,
636  const Event& event, vector<int>& exclude,
637  vector<int>& colSinglet);
638 
639  // Function to check that a set of partons forms a colour singlet
640  // IN Event : Reference event
641  // IN vector<int> : Positions of the partons in the set
642  // OUT bool : Is a colour singlet / is not
643  bool isColSinglet( const Event& event, vector<int> system);
644  // Function to check that a set of partons forms a flavour singlet
645  // IN Event : Reference event
646  // IN vector<int> : Positions of the partons in the set
647  // IN int : Flavour of all the quarks in the set, if
648  // all quarks in a set should have a fixed flavour
649  // OUT bool : Is a flavour singlet / is not
650  bool isFlavSinglet( const Event& event,
651  vector<int> system, int flav=0);
652 
653  // Function to properly colour-connect the radiator to the rest of
654  // the event, as needed during clustering
655  // IN Particle& : Particle to be connected
656  // Particle : Recoiler forming a dipole with Radiator
657  // Event : event to which Radiator shall be appended
658  // OUT true : Radiator could be connected to the event
659  // false : Radiator could not be connected to the
660  // event or the resulting event was
661  // non-valid
662  bool connectRadiator( Particle& radiator, const int radType,
663  const Particle& recoiler, const int recType,
664  const Event& event );
665 
666  // Function to find a colour (anticolour) index in the input event
667  // IN int col : Colour tag to be investigated
668  // int iExclude1 : Identifier of first particle to be excluded
669  // from search
670  // int iExclude2 : Identifier of second particle to be excluded
671  // from search
672  // Event event : event to be searched for colour tag
673  // int type : Tag to define if col should be counted as
674  // colour (type = 1) [->find anti-colour index
675  // contracted with col]
676  // anticolour (type = 2) [->find colour index
677  // contracted with col]
678  // OUT int : Position of particle in event record
679  // contraced with col [0 if col is free tag]
680  int FindCol(int col, int iExclude1, int iExclude2,
681  const Event& event, int type, bool isHardIn);
682 
683  // Function to in the input event find a particle with quantum
684  // numbers matching those of the input particle
685  // IN Particle : Particle to be searched for
686  // Event : Event to be searched in
687  // OUT int : > 0 : Position of matching particle in event
688  // < 0 : No match in event
689  int FindParticle( const Particle& particle, const Event& event,
690  bool checkStatus = true );
691 
692  // Function to check if rad,emt,rec triple is allowed for clustering
693  // IN int rad,emt,rec,partner : Positions (in event record) of the three
694  // particles considered for clustering, and the
695  // correct colour-connected recoiler (=partner)
696  // Event event : Reference event
697  bool allowedClustering( int rad, int emt, int rec, int partner,
698  const Event& event );
699 
700  // Function to check if rad,emt,rec triple is results in
701  // colour singlet radBefore+recBefore
702  // IN int rad,emt,rec : Positions (in event record) of the three
703  // particles considered for clustering
704  // Event event : Reference event
705  bool isSinglett( int rad, int emt, int rec, const Event& event );
706 
707  // Function to check if event is sensibly constructed: Meaning
708  // that all colour indices are contracted and that the charge in
709  // initial and final states matches
710  // IN event : event to be checked
711  // OUT TRUE : event is properly construced
712  // FALSE : event not valid
713  bool validEvent( const Event& event );
714 
715  // Function to check whether two clusterings are identical, used
716  // for finding the history path in the mother -> children direction
717  bool equalClustering( Clustering clus1 , Clustering clus2 );
718 
719  // Chose dummy scale for event construction. By default, choose
720  // sHat for 2->Boson(->2)+ n partons processes and
721  // M_Boson for 2->Boson(->) processes
722  double choseHardScale( const Event& event ) const;
723 
724  // If the state has an incoming hadron return the flavour of the
725  // parton entering the hard interaction. Otherwise return 0
726  int getCurrentFlav(const int) const;
727 
728  // If the state has an incoming hadron return the x-value for the
729  // parton entering the hard interaction. Otherwise return 0.
730  double getCurrentX(const int) const;
731 
732  double getCurrentZ(const int rad, const int rec, const int emt,
733  int idRadBef = 0) const;
734 
735  // Function to compute "pythia pT separation" from Particle input
736  double pTLund(const Event& event, int radAfterBranch, int emtAfterBranch,
737  int recAfterBranch, int showerType, int idRadBef = 0);
738 
739  // Function to return the position of the initial line before (or after)
740  // a single (!) splitting.
741  int posChangedIncoming(const Event& event, bool before);
742 
743  // Function to give back the ratio of PDFs and PDF * splitting kernels
744  // needed to convert a splitting at scale pdfScale, chosen with running
745  // PDFs, to a splitting chosen with PDFs at a fixed scale mu. As needed to
746  // properly count emissions.
747  double pdfFactor( const Event& event, const int type, double pdfScale,
748  double mu );
749 
750  // Function giving the product of splitting kernels and PDFs so that the
751  // resulting flavour is given by flav. This is used as a helper routine
752  // to dgauss
753  double integrand(int flav, double x, double scaleInt, double z);
754 
755  // Function providing a list of possible new flavours after a w emssion
756  // from the input flavour.
757  vector<int> posFlavCKM(int flav);
758 
759  // Check if the new flavour structure is possible.
760  // If clusType is 1 final clustering is assumed, otherwise initial clustering
761  // is assumed.
762  bool checkFlavour(vector<int>& flavCounts, int flavRad, int flavRadBef,
763  int clusType);
764 
765  // Check if the weak recoil structure is allowed.
766  bool checkWeakRecoils(map<int,int> &allowedRecoils, bool isFirst = false);
767 
768  // Find the recoiler for an ISR scattered weak particle.
769 
770  int findISRRecoiler();
771 
772  // Reverse the boost carried out by the ISR.
773  // The three first momenta are given by the ME,
774  // the last two are filled in by this function.
775  void reverseBoostISR(Vec4& pMother, Vec4& pSister, Vec4& pPartner,
776  Vec4& pDaughter, Vec4& pRecoiler, int sign, double eCM, double& phi);
777 
778  // Check if an event reclustered into a 2 -> 2 dijet.
779  // (Only enabled if W reclustering is used).
780  bool isQCD2to2(const Event& event);
781 
782  // Check if an event reclustered into a 2 -> 1 Drell-Yan.
783  // (Only enabled if W reclustering is used).
784  bool isEW2to1(const Event& event);
785 
786  // Set selected child indices.
787  void setSelectedChild();
788 
789  // Setup the weak dipole showers.
790  void setupSimpleWeakShower(int nSteps);
791 
792  // Update weak dipoles after an emission.
793  void transferSimpleWeakShower(vector<int> &mode, vector<Vec4> &mom,
794  vector<int> fermionLines, vector<pair<int,int> > &dipoles, int nSteps);
795 
796  // Update the weak modes after an emission.
797  vector<int> updateWeakModes(vector<int>& weakModes,
798  map<int,int>& stateTransfer);
799 
800  // Update the fermion lines for the 2 -> 2 process. This is needed for
801  // the weak probabilities.
802  vector<int> updateWeakFermionLines(vector<int> fermionLines,
803  map<int,int>& stateTransfer);
804 
805  // Update the list of weak dipoles. This is needed to setup the PS correctly.
806  vector<pair<int,int> > updateWeakDipoles(vector<pair<int,int> > &dipoles,
807  map<int,int>& stateTransfer);
808 
809  // Setup the hard process information needed for calculating weak
810  // probabilities and setting up the shower.
811  void setupWeakHard(vector<int>& mode, vector<int>& fermionLines,
812  vector<Vec4>& mom);
813 
814  // Functions to retrieve scale information from external showers.
815  double getShowerPluginScale(const Event& event, int rad, int emt, int rec,
816  string key, double scalePythia);
817 
818  //----------------------------------------------------------------------//
819  // Class members.
820  //----------------------------------------------------------------------//
821 
822  // The state of the event correponding to this step in the
823  // reconstruction.
824  Event state;
825 
826  // The previous step from which this step has been clustered. If
827  // null, this is the initial step with the n-jet state generated by
828  // the matrix element.
829  History * mother;
830 
831  // The different steps corresponding to possible clusterings of this
832  // state.
833  vector<History *> children;
834 
835  // After a path is selected, store the child index.
836  int selectedChild;
837 
838  // The different paths which have been reconstructed indexed with
839  // the (incremental) corresponding probability. This map is empty
840  // unless this is the initial step (mother == 0).
841  map<double,History *> paths;
842 
843  // The sum of the probabilities of the full paths. This is zero
844  // unless this is the initial step (mother == 0).
845  double sumpath;
846 
847  // The different allowed paths after projection, indexed with
848  // the (incremental) corresponding probability. This map is empty
849  // unless this is the initial step (mother == 0).
850  map<double,History *> goodBranches, badBranches;
851  // The sum of the probabilities of allowed paths after projection. This is
852  // zero unless this is the initial step (mother == 0).
853  double sumGoodBranches, sumBadBranches;
854 
855  // This is set true if an ordered (and complete) path has been found
856  // and inserted in paths.
857  bool foundOrderedPath;
858 
859  // This is set true if a strongly ordered (and complete) path has been found
860  // and inserted in paths.
861  bool foundStronglyOrderedPath;
862 
863  // This is set true if an allowed (according to a user criterion) path has
864  // been found and inserted in paths.
865  bool foundAllowedPath;
866 
867  // This is set true if a complete (with the required number of
868  // clusterings) path has been found and inserted in paths.
869  bool foundCompletePath;
870 
871  // The scale of this step, corresponding to clustering which
872  // constructed the corresponding state (or the merging scale in case
873  // mother == 0).
874  double scale;
875 
876  // Flag indicating if a clustering in the construction of all histories is
877  // the next clustering demanded by inout clusterings in LesHouches 2.0
878  // accord.
879  bool nextInInput;
880 
881  // The probability associated with this step and the previous steps.
882  double prob;
883 
884  // The partons and scale of the last clustering.
885  Clustering clusterIn;
886  int iReclusteredOld, iReclusteredNew;
887 
888  // Flag to include the path amongst allowed paths.
889  bool doInclude;
890 
891  // Pointer to MergingHooks object to get all the settings.
892  MergingHooksPtr mergingHooksPtr;
893 
894  // The default constructor is private.
895  History() : mother(), selectedChild(), sumpath(), sumGoodBranches(),
896  sumBadBranches(), foundOrderedPath(), foundStronglyOrderedPath(),
897  foundAllowedPath(), foundCompletePath(), scale(), nextInInput(), prob(),
898  iReclusteredOld(), iReclusteredNew(), doInclude(), mergingHooksPtr(),
899  particleDataPtr(), infoPtr(), loggerPtr(), showers(), coupSMPtr(),
900  sumScalarPT(), probMaxSave(), depth(), minDepthSave() {}
901 
902  // The copy-constructor is private.
903  History(const History &) {}
904 
905  // The assignment operator is private.
906  History & operator=(const History &) {
907  return *this;
908  }
909 
910  // BeamParticle to get access to PDFs
911  BeamParticle beamA;
912  BeamParticle beamB;
913  // ParticleData needed to initialise the shower AND to get the
914  // correct masses of partons needed in calculation of probability
915  ParticleData* particleDataPtr;
916 
917  // Info object to have access to all information read from LHE file
918  Info* infoPtr;
919 
920  // Logger object.
921  Logger* loggerPtr;
922 
923  // Class for calculation weak shower ME.
924  SimpleWeakShowerMEs simpleWeakShowerMEs;
925 
926  // Pointer to showers, to simplify external clusterings.
927  PartonLevel* showers;
928 
929  // Pointer to standard model couplings.
930  CoupSM* coupSMPtr;
931 
932  // Minimal scalar sum of pT used in Herwig to choose history
933  double sumScalarPT;
934 
935  double probMaxSave;
936  double probMax() {
937  if (mother) return mother->probMax();
938  return probMaxSave;
939  }
940  void updateProbMax(double probIn, bool isComplete = false) {
941  if (mother) mother->updateProbMax(probIn, isComplete);
942  if (!isComplete && !foundCompletePath) return;
943  if (abs(probIn) > probMaxSave) probMaxSave = probIn;
944  }
945 
946  int depth, minDepthSave;
947  int minDepth() {
948  if ( mother ) return mother->minDepth();
949  return minDepthSave;
950  }
951  void updateMinDepth(int depthIn) {
952  if ( mother ) return mother->updateMinDepth(depthIn);
953  minDepthSave = (minDepthSave>0) ? min(minDepthSave,depthIn) : depthIn;
954  }
955 
956 };
957 
958 //==========================================================================
959 
960 } // end namespace Pythia8
961 
962 #endif // end Pythia8_History_H
~Clustering()
Default destructor.
Definition: History.h:71
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 >())
Constructor with input.
Definition: History.h:84
int radBef
The radiator before the splitting.
Definition: History.h:57
Clustering(const Clustering &inSystem)
Copy constructor.
Definition: History.h:74
Definition: StandardModel.h:23
int flavRadBef
The flavour of the radiator prior to the emission.
Definition: History.h:47
double pT() const
Function to return pythia pT scale of clustering.
Definition: History.h:94
double pTscale
The scale associated with this clustering.
Definition: History.h:45
Definition: Info.h:45
bool foundAllowedHistories()
Function to check if any allowed histories were found.
Definition: History.h:207
The Event class holds all info on the generated event.
Definition: Event.h:453
int recBef
The recoiler before the splitting.
Definition: History.h:59
Definition: BeamParticle.h:133
void list() const
print for debug
Definition: History.cc:27
Definition: StandardModel.h:106
Definition: History.h:32
Definition: Logger.h:23
int spinEmt
Spin of the emitted (-1 is left handed, +1 is right handed).
Definition: History.h:51
Definition: Merging.h:33
Definition: Basics.h:386
int partner
The colour connected recoiler (Can be different for ISR)
Definition: History.h:43
int spinRad
Spin of the radiator (-1 is left handed, +1 is right handed).
Definition: History.h:49
int emittor
The emittor parton.
Definition: History.h:39
Definition: PartonLevel.h:45
map< int, int > iPosInMother
Definition: History.h:63
Definition: StandardModel.h:135
Definition: Event.h:32
~History()
The destructor deletes each child.
Definition: History.h:150
Definition: SimpleWeakShowerMEs.h:26
int emitted
The emitted parton location.
Definition: History.h:37
int recoiler
The recoiler parton.
Definition: History.h:41
double phi(const Vec4 &v1, const Vec4 &v2)
phi is azimuthal angle between v1 and v2 around z axis.
Definition: Basics.cc:693
bool foundOrderedHistories()
Function to check if any ordered histories were found.
Definition: History.h:210
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The Pythia class contains the top-level routines to generate an event.
Definition: Pythia.h:71
Event lowestMultProc(const double RN)
Function to get the lowest multiplicity reclustered event.
Definition: History.h:228
int spinRadBef
Spin of the radiator before the splitting.
Definition: History.h:55
bool foundCompleteHistories()
Function to check if any ordered histories were found.
Definition: History.h:213
Clustering()
Default constructor.
Definition: History.h:66
bool validEvent(const Event &event)
The Merging class.
Definition: DireMerging.cc:24
This class holds a map of all ParticleDataEntries.
Definition: ParticleData.h:422
Definition: Basics.h:32
Definition: History.h:116
int spinRec
Spin of the recoiler (-1 is left handed, +1 is right handed).
Definition: History.h:53