]>
Commit | Line | Data |
---|---|---|
7fac8669 | 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 |