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