end-of-line normalization
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8175 / include / History.h
1 // History.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2013 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, 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 "Basics.h"
14 #include "BeamParticle.h"
15 #include "Event.h"
16 #include "Info.h"
17 #include "ParticleData.h"
18 #include "PythiaStdlib.h"
19 #include "Settings.h"
20 #include "PartonLevel.h"
21
22 namespace Pythia8 {
23
24 //==========================================================================
25
26 // Declaration of Clustering class.
27 // This class holds information about one radiator, recoiler, emitted system.
28 // This class is a container class for History class use.
29
30 class Clustering {
31
32 public:
33
34    // The emitted parton location.
35   int emitted;
36   // The emittor parton
37   int emittor;
38   // The recoiler parton.
39   int recoiler;
40   // The colour connected recoiler (Can be different for ISR)
41   int partner;
42   // The scale associated with this clustering.
43   double pTscale;
44
45   // Default constructor
46   Clustering(){
47     emitted   = 0;
48     emittor   = 0;
49     recoiler  = 0;
50     partner   = 0;
51     pTscale   = 0.0;
52   }
53
54   // Default destructor
55   ~Clustering(){}
56
57   // Copy constructor
58   Clustering( const Clustering& inSystem ){
59     emitted  = inSystem.emitted;
60     emittor  = inSystem.emittor;
61     recoiler = inSystem.recoiler;
62     partner  = inSystem.partner;
63     pTscale  = inSystem.pTscale;
64   }
65
66   // Constructor with input
67   Clustering( int emtIn, int radIn, int recIn, int partnerIn,
68     double pTscaleIn ){
69     emitted  = emtIn;
70     emittor  = radIn;
71     recoiler = recIn;
72     partner  = partnerIn;
73     pTscale  = pTscaleIn;
74   }
75
76   // Function to return pythia pT scale of clustering
77   double pT() const { return pTscale; }
78
79   // print for debug
80   void list() const;
81
82 };
83
84 //==========================================================================
85
86 // Declaration of History class
87 //
88 // A History object represents an event in a given step in the CKKW-L
89 // clustering procedure. It defines a tree-like recursive structure,
90 // where the root node represents the state with n jets as given by
91 // the matrix element generator, and is characterized by the member
92 // variable mother being null. The leaves on the tree corresponds to a
93 // fully clustered paths where the original n-jets has been clustered
94 // down to the Born-level state. Also states which cannot be clustered
95 // down to the Born-level are possible - these will be called
96 // incomplete. The leaves are characterized by the vector of children
97 // being empty.
98
99 class History {
100
101 public:
102
103   // The only constructor. Default arguments are used when creating
104   // the initial history node. The \a depth is the maximum number of
105   // clusterings requested. \a scalein is the scale at which the \a
106   // statein was clustered (should be set to the merging scale for the
107   // initial history node. \a beamAIn and beamBIn are needed to
108   // calcutate PDF ratios, \a particleDataIn to have access to the
109   // correct masses of particles. If \a isOrdered is true, the previous
110   // clusterings has been ordered. \a is the PDF ratio for this 
111   // clustering (=1 for FSR clusterings). \a probin is the accumulated
112   // probabilities for the previous clusterings, and \ mothin is the
113   // previous history node (null for the initial node).
114   History( int depth,
115            double scalein,
116            Event statein,
117            Clustering c,
118            MergingHooks* mergingHooksPtrIn,
119            BeamParticle beamAIn,
120            BeamParticle beamBIn,
121            ParticleData* particleDataPtrIn,
122            Info* infoPtrIn,
123            bool isOrdered,
124            bool isStronglyOrdered,
125            bool isAllowed,
126            bool isNextInInput,
127            double probin,
128            History * mothin);
129
130   // The destructor deletes each child.
131   ~History() {
132     for ( int i = 0, N = children.size(); i < N; ++i ) delete children[i];
133   }
134
135   // Function to project paths onto desired paths.
136   bool projectOntoDesiredHistories();
137
138   // For CKKW-L, NL3 and UMEPS:
139   // In the initial history node, select one of the paths according to
140   // the probabilities. This function should be called for the initial
141   // history node.
142   // IN  trialShower*    : Previously initialised trialShower object,
143   //                       to perform trial showering and as
144   //                       repository of pointers to initialise alphaS
145   //     PartonSystems* : PartonSystems object needed to initialise
146   //                      shower objects
147   // OUT double         : (Sukadov) , (alpha_S ratios) , (PDF ratios)
148   double weightTREE(PartonLevel* trial, AlphaStrong * asFSR,
149                     AlphaStrong * asISR, double RN);
150
151   // For default NL3:
152   // Return weight of virtual correction and subtractive for NL3 merging
153   double weightLOOP(PartonLevel* trial, double RN);
154   // Return O(\alpha_s)-term of CKKWL-weight for NL3 merging
155   double weightFIRST(PartonLevel* trial, AlphaStrong* asFSR,
156                   AlphaStrong* asISR, double RN, Rndm* rndmPtr );
157
158   // For UMEPS:
159   double weight_UMEPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,
160                     AlphaStrong * asISR, double RN);
161   double weight_UMEPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,
162                     AlphaStrong * asISR, double RN);
163
164   // For unitary NL3:
165   double weight_UNLOPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,
166                     AlphaStrong * asISR, double RN);
167   double weight_UNLOPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,
168                     AlphaStrong * asISR, double RN);
169   double weight_UNLOPS_LOOP(PartonLevel* trial, double RN);
170   double weight_UNLOPS_SUBTNLO(PartonLevel* trial, double RN);
171   double weight_UNLOPS_CORRECTION( int order, PartonLevel* trial,
172                   AlphaStrong* asFSR, AlphaStrong* asISR, 
173                   double RN, Rndm* rndmPtr );
174
175   // Function to check if any allowed histories were found
176   bool foundAllowedHistories() {
177     return (children.size() > 0 && foundAllowedPath); }
178   // Function to check if any ordered histories were found
179   bool foundOrderedHistories() {
180     return (children.size() > 0 && foundOrderedPath); }
181   // Function to check if any ordered histories were found
182   bool foundCompleteHistories() {
183     return (children.size() > 0 && foundCompletePath); }
184
185   // Function to set the state with complete scales for evolution 
186   void getStartingConditions( const double RN, Event& outState );
187   // Function to get the state with complete scales for evolution 
188   bool getClusteredEvent( const double RN, int nSteps, Event& outState);
189   // Function to get the first reclustered state above the merging scale. 
190   bool getFirstClusteredEventAboveTMS( const double RN, int nDesired,
191     Event& process, int & nPerformed, bool updateProcess = true );
192   // Function to return the depth of the history (i.e. the number of
193   // reclustered splittings)
194   int nClusterings();
195
196   // Function to get the lowest multiplicity reclustered event
197   Event lowestMultProc( const double RN) {
198     // Return lowest multiplicity state
199     return (select(RN)->state);
200   }
201
202   // Calculate and return pdf ratio
203   double getPDFratio( int side, bool forSudakov, bool useHardPDF,
204                       int flavNum, double xNum, double muNum,
205                       int flavDen, double xDen, double muDen);
206
207   // Make Pythia class friend
208   friend class Pythia;
209   // Make Merging class friend
210   friend class Merging;
211
212 private:
213
214   // Number of trial emission to use for calculating the average number of 
215   // emissions
216   static const int NTRIAL;
217
218   // Function to set all scales in the sequence of states. This is a
219   // wrapper routine for setScales and setEventScales methods
220   void setScalesInHistory();
221
222   // Function to find the index (in the mother histories) of the
223   // child history, thus providing a way access the path from both
224   // initial history (mother == 0) and final history (all children == 0)
225   // IN vector<int> : The index of each child in the children vector
226   //                  of the current history node will be saved in
227   //                  this vector
228   // NO OUTPUT
229   void findPath(vector<int>& out);
230
231   // Functions to set the  parton production scales and enforce
232   // ordering on the scales of the respective clusterings stored in
233   // the History node:
234   // Method will start from lowest multiplicity state and move to
235   // higher states, setting the production scales the shower would
236   // have used.
237   // When arriving at the highest multiplicity, the method will switch
238   // and go back in direction of lower states to check and enforce
239   // ordering for unordered histories.
240   // IN vector<int> : Vector of positions of the chosen child
241   //                  in the mother history to allow to move
242   //                  in direction initial->final along path
243   //    bool        : True: Move in direction low->high
244   //                       multiplicity and set production scales
245   //                  False: Move in direction high->low
246   //                       multiplicity and check and enforce
247   //                       ordering
248   // NO OUTPUT
249   void setScales( vector<int> index, bool forward);
250
251   // Function to find a particle in all higher multiplicity events
252   // along the history path and set its production scale to the input
253   // scale
254   // IN  int iPart       : Parton in refEvent to be checked / rescaled
255   //     Event& refEvent : Reference event for iPart
256   //     double scale    : Scale to be set as production scale for
257   //                       unchanged copies of iPart in subsequent steps
258   void scaleCopies(int iPart, const Event& refEvent, double rho);
259
260   // Function to set the OVERALL EVENT SCALES [=state.scale()] to
261   // the scale of the last clustering
262   // NO INPUT
263   // NO OUTPUT
264   void setEventScales();
265
266   // Function to print information on the reconstructed scales in one path.
267   // For debug only.
268   void printScales() { if ( mother ) mother->printScales();
269     cout << " size " << state.size() << " scale " << scale << " clusterIn "
270       << clusterIn.pT() << " state.scale() " << state.scale() << endl; }
271
272   // Function to project paths onto desired paths.
273   bool trimHistories();
274   // Function to tag history for removal.
275   void remove(){ doInclude = false; }
276   // Function to return flag of allowed histories to choose from.
277   bool keep(){ return doInclude; }
278   // Function implementing checks on a paths, for deciding if the path should
279   // be considered valid or not.
280   bool keepHistory();
281   // Function to check if a path is ordered in evolution pT.
282   bool isOrderedPath( double maxscale );
283
284   bool followsInputPath();
285
286   // Function to check if all reconstucted states in a path pass the merging
287   // scale cut.
288   bool allIntermediateAboveRhoMS( double rhoms, bool good = true );
289   // Function to check if any ordered paths were found (and kept).
290   bool foundAnyOrderedPaths();
291
292   // Functions to return the z value of the last ISR splitting
293   // NO INPUT
294   // OUTPUT double : z value of last ISR splitting in history 
295   double zISR();
296
297   // Functions to return the z value of the last FSR splitting
298   // NO INPUT
299   // OUTPUT double : z value of last FSR splitting in history 
300   double zFSR();
301
302   // Functions to return the pT scale of the last ISR splitting
303   // NO INPUT
304   // OUTPUT double : pT scale of last ISR splitting in history 
305   double pTISR();
306
307   // Functions to return the pT scale of the last FSR splitting
308   // NO INPUT
309   // OUTPUT double : pT scale of last FSR splitting in history 
310   double pTFSR();
311
312   // Functions to return the event with nSteps additional partons
313   // INPUT  int   : Number of splittings in the event,
314   //                as counted from core 2->2 process
315   // OUTPUT Event : event with nSteps additional partons 
316   Event clusteredState( int nSteps);
317
318   // Function to choose a path from all paths in the tree
319   // according to their splitting probabilities
320   // IN double    : Random number
321   // OUT History* : Leaf of history path chosen
322   History * select(double rnd);
323
324   // For a full path, find the weight calculated from the ratio of
325   // couplings, the no-emission probabilities, and possible PDF
326   // ratios. This function should only be called for the last history
327   // node of a full path.
328   // IN  TimeShower : Already initialised shower object to be used as
329   //                  trial shower
330   //     double     : alpha_s value used in ME calculation
331   //     double     : Maximal mass scale of the problem (e.g. E_CM)
332   //     AlphaStrong: Initialised shower alpha_s object for FSR alpha_s
333   //                  ratio calculation
334   //     AlphaStrong: Initialised shower alpha_s object for ISR alpha_s
335   //                  ratio calculation (can be different from previous)
336   double weightTree(PartonLevel* trial, double as0, double maxscale,
337     double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,
338     double& asWeight, double& pdfWeight);
339
340   // Function to return the \alpha_s-ratio part of the CKKWL weight.
341   double weightTreeALPHAS( double as0, AlphaStrong * asFSR, 
342     AlphaStrong * asISR );
343   // Function to return the PDF-ratio part of the CKKWL weight.
344   double weightTreePDFs( double maxscale, double pdfScale );
345   // Function to return the no-emission probability part of the CKKWL weight.
346   double weightTreeEmissions( PartonLevel* trial, int type, int njetMax, 
347     double maxscale );
348
349   // Function to generate the O(\alpha_s)-term of the CKKWL-weight
350   double weightFirst(PartonLevel* trial, double as0, double muR,
351     double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR, Rndm* rndmPtr );
352
353   // Function to generate the O(\alpha_s)-term of the \alpha_s-ratios
354   // appearing in the CKKWL-weight.
355   double weightFirstALPHAS( double as0, double muR, AlphaStrong * asFSR,
356     AlphaStrong * asISR);
357   // Function to generate the O(\alpha_s)-term of the PDF-ratios
358   // appearing in the CKKWL-weight.
359   double weightFirstPDFs( double as0, double maxscale, double pdfScale,
360     Rndm* rndmPtr );
361   // Function to generate the O(\alpha_s)-term of the no-emission
362   // probabilities appearing in the CKKWL-weight.
363   double weightFirstEmissions(PartonLevel* trial, double as0, double maxscale,
364     AlphaStrong * asFSR, AlphaStrong * asISR, bool fixpdf, bool fixas );
365
366   // Function to return the default factorisation scale of the hard process.
367   double hardFacScale(const Event& event);
368   // Function to return the default renormalisation scale of the hard process.
369   double hardRenScale(const Event& event);
370
371   // Perform a trial shower using the \a pythia object between
372   // maxscale down to this scale and return the corresponding Sudakov
373   // form factor.
374   // IN  trialShower : Shower object used as trial shower
375   //     double     : Maximum scale for trial shower branching
376   // OUT  0.0       : trial shower emission outside allowed pT range
377   //      1.0       : trial shower successful (any emission was below
378   //                  the minimal scale )
379   double doTrialShower(PartonLevel* trial, int type, double maxscale, 
380     double minscale = 0.);
381
382   // Function to bookkeep the indices of weights generated in countEmissions
383   bool updateind(vector<int> & ind, int i, int N);
384
385   // Function to count number of emissions between two scales for NLO merging
386   vector<double> countEmissions(PartonLevel* trial, double maxscale,
387     double minscale, int showerType, double as0, AlphaStrong * asFSR,
388     AlphaStrong * asISR, int N, bool fixpdf, bool fixas);
389
390   // Function to integrate PDF ratios between two scales over x and t,
391   // where the PDFs are always evaluated at the lower t-integration limit
392   double monteCarloPDFratios(int flav, double x, double maxScale,
393            double minScale, double pdfScale, double asME, Rndm* rndmPtr);
394
395   // Default: Check if a ordered (and complete) path has been found in
396   // the initial node, in which case we will no longer be interested in
397   // any unordered paths.
398   bool onlyOrderedPaths();
399
400   // Check if a strongly ordered (and complete) path has been found in the
401   // initial node, in which case we will no longer be interested in
402   // any unordered paths.
403   bool onlyStronglyOrderedPaths();
404
405   // Check if an allowed (according to user-criterion) path has been found in 
406   // the initial node, in which case we will no longer be interested in
407   // any forbidden paths.
408   bool onlyAllowedPaths();
409
410   // When a full path has been found, register it with the initial
411   // history node.
412   // IN  History : History to be registered as path
413   //     bool    : Specifying if clusterings so far were ordered
414   //     bool    : Specifying if path is complete down to 2->2 process
415   // OUT true if History object forms a plausible path (eg prob>0 ...)
416   bool registerPath(History & l, bool isOrdered, bool isStronglyOrdered,
417          bool isAllowed, bool isComplete);
418
419   // For the history-defining state (and if necessary interfering
420   // states), find all possible clusterings.
421   // NO INPUT
422   // OUT vector of all (rad,rec,emt) systems
423   vector<Clustering> getAllQCDClusterings();
424
425   // For one given state, find all possible clusterings.
426   // IN  Event : state to be investigated
427   // OUT vector of all (rad,rec,emt) systems in the state
428   vector<Clustering> getQCDClusterings( const Event& event);
429
430   // Function to construct (rad,rec,emt) triples from the event
431   // IN  int   : Position of Emitted in event record for which
432   //             dipoles should be constructed
433   //     int   : Colour topogy to be tested
434   //             1= g -> qqbar, causing 2 -> 2 dipole splitting
435   //             2= q(bar) -> q(bar) g && g -> gg,
436   //              causing a 2 -> 3 dipole splitting
437   //     Event : event record to be checked for ptential partners
438   // OUT vector of all allowed radiator+recoiler+emitted triples
439   vector<Clustering> findQCDTriple (int EmtTagIn, int colTopIn, 
440                        const Event& event, vector<int> PosFinalPartn,
441                        vector <int> PosInitPartn );
442
443   vector<Clustering> getAllEWClusterings();
444   vector<Clustering> getEWClusterings( const Event& event);
445   vector<Clustering> findEWTriple( int EmtTagIn, const Event& event,
446                        vector<int> PosFinalPartn );
447
448   vector<Clustering> getAllSQCDClusterings();
449   vector<Clustering> getSQCDClusterings( const Event& event);
450   vector<Clustering> findSQCDTriple (int EmtTagIn, int colTopIn, 
451                        const Event& event, vector<int> PosFinalPartn,
452                        vector <int> PosInitPartn );
453
454   // Calculate and return the probability of a clustering.
455   // IN  Clustering : rad,rec,emt - System for which the splitting
456   //                  probability should be calcuated
457   // OUT splitting probability
458   double getProb(const Clustering & SystemIn);
459
460   // Set up the beams (fill the beam particles with the correct
461   // current incoming particles) to allow calculation of splitting
462   // probability.
463   // For interleaved evolution, set assignments dividing PDFs into
464   // sea and valence content. This assignment is, until a history path
465   // is chosen and a first trial shower performed, not fully correct
466   // (since content is chosen form too high x and too low scale). The
467   // assignment used for reweighting will be corrected after trial
468   // showering
469   void setupBeams();
470
471   // Calculate the PDF ratio used in the argument of the no-emission
472   // probability.
473   double pdfForSudakov();
474
475   // Calculate the hard process matrix element to include in the selection
476   // probability.
477   double hardProcessME( const Event& event);
478
479   // Perform the clustering of the current state and return the
480   // clustered state.
481   // IN Clustering : rad,rec,emt triple to be clustered to two partons
482   // OUT clustered state
483   Event cluster(const Clustering & inSystem);
484
485   // Function to get the flavour of the radiator before the splitting
486   // for clustering
487   // IN  int   : Position of the radiator after the splitting, in the event
488   //     int   : Position of the emitted after the splitting, in the event
489   //     Event : Reference event   
490   // OUT int   : Flavour of the radiator before the splitting
491   int getRadBeforeFlav(const int RadAfter, const int EmtAfter,
492         const Event& event);
493
494   // Function to get the colour of the radiator before the splitting
495   // for clustering
496   // IN  int   : Position of the radiator after the splitting, in the event
497   //     int   : Position of the emitted after the splitting, in the event
498   //     Event : Reference event   
499   // OUT int   : Colour of the radiator before the splitting
500   int getRadBeforeCol(const int RadAfter, const int EmtAfter,
501         const Event& event);
502
503   // Function to get the anticolour of the radiator before the splitting
504   // for clustering
505   // IN  int   : Position of the radiator after the splitting, in the event
506   //     int   : Position of the emitted after the splitting, in the event
507   //     Event : Reference event   
508   // OUT int   : Anticolour of the radiator before the splitting
509   int getRadBeforeAcol(const int RadAfter, const int EmtAfter,
510         const Event& event);
511
512   // Function to get the parton connected to in by a colour line
513   // IN  int   : Position of parton for which partner should be found
514   //     Event : Reference event   
515   // OUT int   : If a colour line connects the "in" parton with another
516   //             parton, return the Position of the partner, else return 0
517   int getColPartner(const int in, const Event& event);
518
519   // Function to get the parton connected to in by an anticolour line
520   // IN  int   : Position of parton for which partner should be found
521   //     Event : Reference event   
522   // OUT int   : If an anticolour line connects the "in" parton with another
523   //             parton, return the Position of the partner, else return 0
524   int getAcolPartner(const int in, const Event& event);
525
526   // Function to get the list of partons connected to the particle
527   // formed by reclusterinf emt and rad by colour and anticolour lines
528   // IN  int          : Position of radiator in the clustering
529   // IN  int          : Position of emitted in the clustering
530   //     Event        : Reference event   
531   // OUT vector<int>  : List of positions of all partons that are connected
532   //                    to the parton that will be formed
533   //                    by clustering emt and rad.
534   vector<int> getReclusteredPartners(const int rad, const int emt,
535     const Event& event);
536
537   // Function to extract a chain of colour-connected partons in
538   // the event 
539   // IN     int          : Type of parton from which to start extracting a
540   //                       parton chain. If the starting point is a quark
541   //                       i.e. flavType = 1, a chain of partons that are
542   //                       consecutively connected by colour-lines will be
543   //                       extracted. If the starting point is an antiquark
544   //                       i.e. flavType =-1, a chain of partons that are
545   //                       consecutively connected by anticolour-lines
546   //                       will be extracted.
547   // IN      int         : Position of the parton from which a
548   //                       colour-connected chain should be derived
549   // IN      Event       : Refernence event
550   // IN/OUT  vector<int> : Partons that should be excluded from the search.
551   // OUT     vector<int> : Positions of partons along the chain
552   // OUT     bool        : Found singlet / did not find singlet
553   bool getColSinglet(const int flavType, const int iParton,
554     const Event& event, vector<int>& exclude,
555     vector<int>& colSinglet);
556
557   // Function to check that a set of partons forms a colour singlet
558   // IN  Event       : Reference event
559   // IN  vector<int> : Positions of the partons in the set
560   // OUT bool        : Is a colour singlet / is not 
561   bool isColSinglet( const Event& event, vector<int> system);
562   // Function to check that a set of partons forms a flavour singlet
563   // IN  Event       : Reference event
564   // IN  vector<int> : Positions of the partons in the set
565   // IN  int         : Flavour of all the quarks in the set, if
566   //                   all quarks in a set should have a fixed flavour
567   // OUT bool        : Is a flavour singlet / is not 
568   bool isFlavSinglet( const Event& event,
569     vector<int> system, int flav=0);
570
571   // Function to properly colour-connect the radiator to the rest of
572   // the event, as needed during clustering
573   // IN  Particle& : Particle to be connected
574   //     Particle  : Recoiler forming a dipole with Radiator
575   //     Event     : event to which Radiator shall be appended
576   // OUT true               : Radiator could be connected to the event
577   //     false              : Radiator could not be connected to the
578   //                          event or the resulting event was
579   //                          non-valid
580   bool connectRadiator( Particle& Radiator, const int RadType,
581                         const Particle& Recoiler, const int RecType, 
582                         const Event& event );
583
584   // Function to find a colour (anticolour) index in the input event
585   // IN  int col       : Colour tag to be investigated
586   //     int iExclude1 : Identifier of first particle to be excluded
587   //                     from search
588   //     int iExclude2 : Identifier of second particle to be excluded
589   //                     from  search
590   //     Event event   : event to be searched for colour tag
591   //     int type      : Tag to define if col should be counted as
592   //                      colour (type = 1) [->find anti-colour index
593   //                                         contracted with col]
594   //                      anticolour (type = 2) [->find colour index
595   //                                         contracted with col]
596   // OUT int           : Position of particle in event record
597   //                     contraced with col [0 if col is free tag]
598   int FindCol(int col, int iExclude1, int iExclude2,
599               const Event& event, int type, bool isHardIn);
600
601   // Function to in the input event find a particle with quantum
602   // numbers matching those of the input particle
603   // IN  Particle : Particle to be searched for
604   //     Event    : Event to be searched in
605   // OUT int      : > 0 : Position of matching particle in event
606   //                < 0 : No match in event
607   int FindParticle( const Particle& particle, const Event& event,
608     bool checkStatus = true );
609
610   // Function to check if rad,emt,rec triple is allowed for clustering
611   // IN int rad,emt,rec,partner : Positions (in event record) of the three
612   //                      particles considered for clustering, and the
613   //                      correct colour-connected recoiler (=partner)
614   //    Event event     : Reference event                  
615   bool allowedClustering( int rad, int emt, int rec, int partner,
616     const Event& event );
617
618   // Function to check if rad,emt,rec triple is results in
619   // colour singlet radBefore+recBefore
620   // IN int rad,emt,rec : Positions (in event record) of the three
621   //                      particles considered for clustering
622   //    Event event     : Reference event                  
623   bool isSinglett( int rad, int emt, int rec, const Event& event );
624
625   // Function to check if event is sensibly constructed: Meaning
626   // that all colour indices are contracted and that the charge in
627   // initial and final states matches
628   // IN  event : event to be checked
629   // OUT TRUE  : event is properly construced
630   //     FALSE : event not valid
631   bool validEvent( const Event& event );
632
633   // Function to check whether two clusterings are identical, used
634   // for finding the history path in the mother -> children direction
635   bool equalClustering( Clustering clus1 , Clustering clus2 );
636
637   // Chose dummy scale for event construction. By default, choose
638   //     sHat     for 2->Boson(->2)+ n partons processes and
639   //     M_Boson  for 2->Boson(->)             processes 
640   double choseHardScale( const Event& event ) const;
641
642   // If the state has an incoming hadron return the flavour of the
643   // parton entering the hard interaction. Otherwise return 0
644   int getCurrentFlav(const int) const;
645
646    // If the state has an incoming hadron return the x-value for the
647    // parton entering the hard interaction. Otherwise return 0.
648   double getCurrentX(const int) const;
649
650   double getCurrentZ(const int rad, const int rec, const int emt) const;
651
652   // Function to compute "pythia pT separation" from Particle input
653   double pTLund(const Particle& RadAfterBranch,const Particle& EmtAfterBranch,
654                 const Particle& RecAfterBranch, int ShowerType);
655
656   // Function to return the position of the initial line before (or after)
657   // a single (!) splitting.
658   int posChangedIncoming(const Event& event, bool before);
659
660   // Function to give back the ratio of PDFs and PDF * splitting kernels
661   // needed to convert a splitting at scale pdfScale, chosen with running
662   // PDFs, to a splitting chosen with PDFs at a fixed scale mu. As needed to
663   // properly count emissions.
664   double pdfFactor( const Event& event, const int type, double pdfScale,
665     double mu );
666
667   // Function giving the product of splitting kernels and PDFs so that the
668   // resulting flavour is given by flav. This is used as a helper routine 
669   // to dgauss
670   double integrand(int flav, double x, double scaleInt, double z);
671
672   //----------------------------------------------------------------------//
673   // Class members.
674   //----------------------------------------------------------------------//
675
676   // The state of the event correponding to this step in the
677   // reconstruction.
678   Event state;
679
680   // The previous step from which this step has been clustered. If
681   // null, this is the initial step with the n-jet state generated by
682   // the matrix element.
683   History * mother;
684
685   // The different steps corresponding to possible clusterings of this
686   // state.
687   vector<History *> children;
688
689   // The different paths which have been reconstructed indexed with
690   // the (incremental) corresponding probability. This map is empty
691   // unless this is the initial step (mother == 0).
692   map<double,History *> paths;
693
694   // The sum of the probabilities of the full paths. This is zero
695   // unless this is the initial step (mother == 0).
696   double sumpath;
697
698   // The different allowed paths after projection, indexed with
699   // the (incremental) corresponding probability. This map is empty
700   // unless this is the initial step (mother == 0).
701   map<double,History *> goodBranches, badBranches;
702   // The sum of the probabilities of allowed paths after projection. This is
703   // zero unless this is the initial step (mother == 0).
704   double sumGoodBranches, sumBadBranches;
705
706   // This is set true if an ordered (and complete) path has been found
707   // and inserted in paths.
708   bool foundOrderedPath;
709
710   // This is set true if a strongly ordered (and complete) path has been found
711   // and inserted in paths.
712   bool foundStronglyOrderedPath;
713
714   // This is set true if an allowed (according to a user criterion) path has 
715   // been found and inserted in paths.
716   bool foundAllowedPath;
717
718   // This is set true if a complete (with the required number of
719   // clusterings) path has been found and inserted in paths.
720   bool foundCompletePath;
721
722   // The scale of this step, corresponding to clustering which
723   // constructed the corresponding state (or the merging scale in case
724   // mother == 0).
725   double scale;
726
727   // Flag indicating if a clustering in the construction of all histories is
728   // the next clustering demanded by inout clusterings in LesHouches 2.0
729   // accord.
730   bool nextInInput;
731
732   // The probability associated with this step and the previous steps.
733   double prob;
734
735   // The partons and scale of the last clustering.
736   Clustering clusterIn;
737   int iReclusteredOld, iReclusteredNew;
738
739   // Flag to include the path amongst allowed paths.
740   bool doInclude;
741
742   // Pointer to MergingHooks object to get all the settings.
743   MergingHooks* mergingHooksPtr;
744
745    // The default constructor is private.
746   History() {}
747
748   // The copy-constructor is private.
749   History(const History &) {}
750
751   // The assignment operator is private.
752   History & operator=(const History &) {
753     return *this;
754   }
755
756   // BeamParticle to get access to PDFs
757   BeamParticle beamA;
758   BeamParticle beamB;
759   // ParticleData needed to initialise the shower AND to get the
760   // correct masses of partons needed in calculation of probability
761   ParticleData* particleDataPtr;
762
763   // Info object to have access to all information read from LHE file
764   Info* infoPtr;
765
766   // Minimal scalar sum of pT used in Herwig to choose history
767   double sumScalarPT;
768
769 };
770
771 //==========================================================================
772
773 } // end namespace Pythia8
774
775 #endif // end Pythia8_History_H