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