]>
Commit | Line | Data |
---|---|---|
63ba5337 | 1 | // History.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2012 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 | // Function definitions (not found in the header) for the | |
8 | // Clustering and History classes. | |
9 | ||
10 | #include "History.h" | |
11 | ||
12 | namespace Pythia8 { | |
13 | ||
14 | //========================================================================== | |
15 | ||
16 | // The Clustering class. | |
17 | ||
18 | //-------------------------------------------------------------------------- | |
19 | ||
20 | // Declaration of Clustering class | |
21 | // This class holds information about one radiator, recoiler, | |
22 | // emitted system. | |
23 | // This class is a container class for History class use. | |
24 | ||
25 | // print for debug | |
26 | void Clustering::list() const { | |
27 | cout << " emt " << emitted | |
28 | << " rad " << emittor | |
29 | << " rec " << recoiler | |
30 | << " partner " << partner | |
31 | << " pTscale " << pTscale << endl; | |
32 | } | |
33 | ||
34 | //========================================================================== | |
35 | ||
36 | // The History class. | |
37 | ||
38 | // A History object represents an event in a given step in the CKKW-L | |
39 | // clustering procedure. It defines a tree-like recursive structure, | |
40 | // where the root node represents the state with n jets as given by | |
41 | // the matrix element generator, and is characterized by the member | |
42 | // variable mother being null. The leaves on the tree corresponds to a | |
43 | // fully clustered paths where the original n-jets has been clustered | |
44 | // down to the Born-level state. Also states which cannot be clustered | |
45 | // down to the Born-level are possible - these will be called | |
46 | // incomplete. The leaves are characterized by the vector of children | |
47 | // being empty. | |
48 | ||
49 | //-------------------------------------------------------------------------- | |
50 | ||
51 | // Declaration of History class | |
52 | // The only constructor. Default arguments are used when creating | |
53 | // the initial history node. The \a depth is the maximum number of | |
54 | // clusterings requested. \a scalein is the scale at which the \a | |
55 | // statein was clustered (should be set to the merging scale for the | |
56 | // initial history node. \a beamAIn and beamBIn are needed to | |
57 | // calcutate PDF ratios, \a particleDataIn to have access to the | |
58 | // correct masses of particles. If \a isOrdered is true, the previous | |
59 | // clusterings has been ordered. \a is the PDF ratio for this | |
60 | // clustering (=1 for FSR clusterings). \a probin is the accumulated | |
61 | // probabilities for the previous clusterings, and \ mothin is the | |
62 | // previous history node (null for the initial node). | |
63 | ||
64 | History::History( int depth, | |
65 | double scalein, | |
66 | Event statein, | |
67 | Clustering c, | |
68 | MergingHooks* mergingHooksPtrIn, | |
69 | BeamParticle beamAIn, | |
70 | BeamParticle beamBIn, | |
71 | ParticleData* particleDataPtrIn, | |
72 | Info* infoPtrIn, | |
73 | bool isOrdered = true, | |
74 | bool isUnordered = true, | |
75 | bool isStronglyOrdered = true, | |
76 | bool isAllowed = true, | |
77 | double probin = 1.0, | |
78 | History * mothin = 0) | |
79 | : state(statein), | |
80 | mother(mothin), | |
81 | sumpath(0.0), | |
82 | sumGoodBranches(0.0), | |
83 | sumBadBranches(0.0), | |
84 | foundOrderedPath(false), | |
85 | foundUnorderedPath(false), | |
86 | foundStronglyOrderedPath(false), | |
87 | foundAllowedPath(false), | |
88 | foundCompletePath(false), | |
89 | scale(scalein), | |
90 | prob(probin), | |
91 | clusterIn(c), | |
92 | iReclusteredOld(0), | |
93 | doInclude(true), | |
94 | mergingHooksPtr(mergingHooksPtrIn), | |
95 | beamA(beamAIn), | |
96 | beamB(beamBIn), | |
97 | particleDataPtr(particleDataPtrIn), | |
98 | infoPtr(infoPtrIn) | |
99 | { | |
100 | ||
101 | // Initialise beam particles | |
102 | setupBeams(); | |
103 | // Update probability with PDF ratio | |
104 | if (mother && mergingHooksPtr->includeRedundant()) prob *= pdfForSudakov(); | |
105 | ||
106 | // Minimal scalar sum of pT used in Herwig to choose history | |
107 | // Keep track of scalar PT | |
108 | if (mother) { | |
109 | double acoll = (mother->state[clusterIn.emittor].isFinal()) | |
110 | ? mergingHooksPtr->herwigAcollFSR() | |
111 | : mergingHooksPtr->herwigAcollISR(); | |
112 | sumScalarPT = mother->sumScalarPT + acoll*scale; | |
113 | } else | |
114 | sumScalarPT = 0.0; | |
115 | ||
116 | // Remember reclustered radiator in lower multiplicity state | |
117 | if ( mother ) iReclusteredOld = mother->iReclusteredNew; | |
118 | ||
119 | // Check if more steps should be taken. | |
120 | int nFinalP = 0; | |
121 | int nFinalW = 0; | |
122 | for ( int i = 0; i < int(state.size()); ++i ) | |
123 | if ( state[i].isFinal() ) { | |
124 | if ( state[i].colType() != 0 ) | |
125 | nFinalP++; | |
126 | if ( state[i].idAbs() == 24 ) | |
127 | nFinalW++; | |
128 | } | |
129 | if ( mergingHooksPtr->doWClustering() | |
130 | && nFinalP == 2 && nFinalW == 0 ) depth = 0; | |
131 | ||
132 | // If this is not the fully clustered state, try to find possible | |
133 | // QCD clusterings. | |
134 | vector<Clustering> clusterings; | |
135 | if ( depth > 0 ) clusterings = getAllQCDClusterings(); | |
136 | ||
137 | // If necessary, try to find possible EW clusterings. | |
138 | vector<Clustering> clusteringsEW; | |
139 | if ( depth > 0 && mergingHooksPtr->doWClustering() ) | |
140 | clusteringsEW = getAllEWClusterings(); | |
141 | if ( !clusteringsEW.empty() ) { | |
142 | clusterings.resize(0); | |
143 | clusterings.insert( clusterings.end(), clusteringsEW.begin(), | |
144 | clusteringsEW.end() ); | |
145 | } | |
146 | // If no clusterings were found, the recursion is done and we | |
147 | // register this node. | |
148 | if ( clusterings.empty() ) { | |
149 | registerPath( *this, isOrdered, isUnordered, isStronglyOrdered, | |
150 | isAllowed, depth == 0 ); | |
151 | return; | |
152 | } | |
153 | ||
154 | // Now we sort the possible clusterings so that we try the | |
155 | // smallest scale first. | |
156 | multimap<double, Clustering *> sorted; | |
157 | for ( int i = 0, N = clusterings.size(); i < N; ++i ) { | |
158 | sorted.insert(make_pair(clusterings[i].pT(), &clusterings[i])); | |
159 | } | |
160 | ||
161 | for ( multimap<double, Clustering *>::iterator it = sorted.begin(); | |
162 | it != sorted.end(); ++it ) { | |
163 | ||
164 | // If this path is not strongly ordered and we already have found an | |
165 | // ordered path, then we don't need to continue along this path. | |
166 | bool stronglyOrdered = isStronglyOrdered; | |
167 | if ( mergingHooksPtr->enforceStrongOrdering() | |
168 | && ( !stronglyOrdered | |
169 | || ( mother && ( it->first < | |
170 | mergingHooksPtr->scaleSeparationFactor()*scale ) ))) { | |
171 | if ( onlyStronglyOrderedPaths() ) continue; | |
172 | stronglyOrdered = false; | |
173 | } | |
174 | ||
175 | bool ordered = isOrdered; | |
176 | bool unordered = isUnordered; | |
177 | if ( mergingHooksPtr->orderInRapidity() | |
178 | && mergingHooksPtr->orderHistories() ) { | |
179 | // Get new z value | |
180 | double z = getCurrentZ((*it->second).emittor, | |
181 | (*it->second).recoiler,(*it->second).emitted); | |
182 | // Get z value of splitting that produced this state | |
183 | double zOld = (!mother) ? 0. : mother->getCurrentZ(clusterIn.emittor, | |
184 | clusterIn.recoiler,clusterIn.emitted); | |
185 | // If this path is not ordered in pT and y, and we already have found | |
186 | // an ordered path, then we don't need to continue along this path. | |
187 | if ( !ordered || ( mother && (it->first < scale | |
188 | || it->first < pow(1. - z,2) / (z * (1. - zOld ))*scale ))) { | |
189 | if ( onlyOrderedPaths() ) continue; | |
190 | ordered = false; | |
191 | } | |
192 | ||
193 | } else if ( mergingHooksPtr->orderHistories() ) { | |
194 | // If this path is not ordered in pT and we already have found an | |
195 | // ordered path, then we don't need to continue along this path. | |
196 | if ( !ordered || ( mother && (it->first < scale) ) ) { | |
197 | if ( onlyOrderedPaths() ) continue; | |
198 | ordered = false; | |
199 | } | |
200 | } else if ( mergingHooksPtr->forceUnorderedHistories() ) { | |
201 | // If this path is not unordered in pT and we already have found an | |
202 | // unordered path, then we don't need to continue along this path. | |
203 | if ( !unordered || ( mother && (it->first > scale) ) ) { | |
204 | if ( onlyUnorderedPaths() ) continue; | |
205 | unordered = false; | |
206 | } | |
207 | } | |
208 | ||
209 | // Check if reclustered state should be disallowed | |
210 | bool doCut = mergingHooksPtr->canCutOnRecState() | |
211 | || mergingHooksPtr->allowCutOnRecState(); | |
212 | bool allowed = isAllowed; | |
213 | ||
214 | if ( doCut | |
215 | && mergingHooksPtr->doCutOnRecState(cluster(*it->second)) ) { | |
216 | if ( onlyAllowedPaths() ) continue; | |
217 | allowed = false; | |
218 | } | |
219 | ||
220 | // Perform the clustering and recurse and construct the next | |
221 | // history node. | |
222 | children.push_back(new History(depth - 1,it->first,cluster(*it->second), | |
223 | *it->second, mergingHooksPtr, beamA, beamB, particleDataPtr, | |
224 | infoPtr, ordered, unordered, stronglyOrdered, allowed, | |
225 | prob*getProb(*it->second), this )); | |
226 | } | |
227 | } | |
228 | ||
229 | //-------------------------------------------------------------------------- | |
230 | ||
231 | // Function to project all possible paths onto only the desired paths. | |
232 | ||
233 | bool History::projectOntoDesiredHistories() { | |
234 | // At the moment, only trim histories. | |
235 | return trimHistories(); | |
236 | } | |
237 | ||
238 | //-------------------------------------------------------------------------- | |
239 | ||
240 | // In the initial history node, select one of the paths according to | |
241 | // the probabilities. This function should be called for the initial | |
242 | // history node. | |
243 | // IN trialShower* : Previously initialised trialShower object, | |
244 | // to perform trial showering and as | |
245 | // repository of pointers to initialise alphaS | |
246 | // PartonSystems* : PartonSystems object needed to initialise | |
247 | // shower objects | |
248 | // OUT double : (Sukadov) , (alpha_S ratios) , (PDF ratios) | |
249 | ||
250 | double History::weightTREE(PartonLevel* trial, AlphaStrong * asFSR, | |
251 | AlphaStrong * asISR, double RN) { | |
252 | ||
253 | // Read alpha_S in ME calculation and maximal scale (eCM) | |
254 | double asME = infoPtr->alphaS(); | |
255 | double maxScale = (foundCompletePath) ? infoPtr->eCM() : infoPtr->QFac(); | |
256 | // Select a path of clusterings | |
257 | History * selected = select(RN); | |
258 | // Set scales in the states to the scales pythia would have set | |
259 | selected->setScalesInHistory(); | |
260 | // Get weight. | |
261 | double asWeight = 1.; | |
262 | double pdfWeight = 1.; | |
263 | ||
264 | // Do trial shower, calculation of alpha_S ratios, PDF ratios | |
265 | double wt = selected->weightTree(trial, asME, maxScale, | |
266 | selected->clusterIn.pT(), asFSR, asISR, asWeight, pdfWeight); | |
267 | ||
268 | // Set hard process renormalisation scale to default Pythia value. | |
269 | bool resetScales = mergingHooksPtr->resetHardQRen(); | |
270 | // For pure QCD dijet events, evaluate the coupling of the hard process at | |
271 | // a more reasonable pT, rather than evaluation \alpha_s at a fixed | |
272 | // arbitrary scale. | |
273 | if ( resetScales | |
274 | && mergingHooksPtr->getProcessString().compare("pp>jj") == 0) { | |
275 | // Reset to a running coupling. Here we choose FSR for simplicity. | |
276 | double newQ2Ren = pow2( selected->hardRenScale() ); | |
277 | double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME; | |
278 | asWeight *= pow(runningCoupling,2); | |
279 | } | |
280 | ||
281 | // For prompt photon events, evaluate the coupling of the hard process at | |
282 | // a more reasonable pT, rather than evaluation \alpha_s at a fixed | |
283 | // arbitrary scale. | |
284 | if ( resetScales | |
285 | && mergingHooksPtr->getProcessString().compare("pp>aj") == 0) { | |
286 | // Reset to a running coupling. In prompt photon always ISR. | |
287 | double newQ2Ren = pow2( selected->hardRenScale() ); | |
288 | double runningCoupling = | |
289 | (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME; | |
290 | asWeight *= runningCoupling; | |
291 | } | |
292 | ||
293 | // Done | |
294 | return (wt*asWeight*pdfWeight); | |
295 | } | |
296 | ||
297 | //-------------------------------------------------------------------------- | |
298 | ||
299 | // Function to set the state with complete scales for evolution. | |
300 | ||
301 | void History::getStartingConditions( const double RN, Event& outState ) { | |
302 | ||
303 | // Select the history | |
304 | History * selected = select(RN); | |
305 | ||
306 | ||
307 | // Copy the output state | |
308 | outState = state; | |
309 | ||
310 | // Set the scale of the lowest order process | |
311 | if (!selected->mother) { | |
312 | int nFinal = 0; | |
313 | for(int i=0; i < int(outState.size()); ++i) | |
314 | if (outState[i].isFinal()) nFinal++; | |
315 | if (nFinal <=2) | |
316 | outState.scale(infoPtr->QFac()); | |
317 | // // If the hard process has a resonance decay which is not | |
318 | // // corrected (e.g. for pp -> (V->jj) + jets merging), set | |
319 | // // factorization scale as starting scale | |
320 | // if (mergingHooksPtr->hardProcess.hasResInProc()) | |
321 | // outState.scale(infoPtr->QFac()); | |
322 | // // If the hard process has a resonance decay which is | |
323 | // // corrected (e.g. for e+e- -> 2 + n jets merging), set | |
324 | // // half the intermediate mass as starting scale | |
325 | // else | |
326 | // outState.scale(0.5*outState[5].m()); | |
327 | ||
328 | // Save information on last splitting, to allow the next | |
329 | // emission in the shower to have smaller rapidity with | |
330 | // respect to the last ME splitting | |
331 | // For hard process, use dummy values | |
332 | if (mergingHooksPtr->getNumberOfClusteringSteps(state) == 0) { | |
333 | infoPtr->zNowISR(0.5); | |
334 | infoPtr->pT2NowISR(pow(state[0].e(),2)); | |
335 | infoPtr->hasHistory(true); | |
336 | // For incomplete process, try to use real values | |
337 | } else { | |
338 | infoPtr->zNowISR(selected->zISR()); | |
339 | infoPtr->pT2NowISR(pow(selected->pTISR(),2)); | |
340 | infoPtr->hasHistory(true); | |
341 | } | |
342 | ||
343 | } else { | |
344 | // Save information on last splitting, to allow the next | |
345 | // emission in the shower to have smaller rapidity with | |
346 | // respect to the last ME splitting | |
347 | infoPtr->zNowISR(selected->zISR()); | |
348 | infoPtr->pT2NowISR(pow(selected->pTISR(),2)); | |
349 | infoPtr->hasHistory(true); | |
350 | } | |
351 | ||
352 | } | |
353 | ||
354 | //-------------------------------------------------------------------------- | |
355 | ||
356 | // Function to set the state with complete scales for evolution. | |
357 | ||
358 | bool History::getClusteredEvent( const double RN, Event& outState, | |
359 | int nSteps) { | |
360 | ||
361 | // Select history | |
362 | History * selected = select(RN); | |
363 | // Set scales in the states to the scales pythia would have set | |
364 | // (Only needed if not done before in calculation of weights or | |
365 | // setting of starting conditions) | |
366 | selected->setScalesInHistory(); | |
367 | // If the history does not allow for nSteps clusterings (e.g. because the | |
368 | // history is incomplete), return false | |
369 | if (nSteps > selected->nClusterings()) return false; | |
370 | // Return event with nSteps-1 additional partons (i.e. recluster the last | |
371 | // splitting) and copy the output state | |
372 | outState = selected->clusteredState(nSteps-1); | |
373 | // Done | |
374 | return true; | |
375 | ||
376 | } | |
377 | ||
378 | //-------------------------------------------------------------------------- | |
379 | ||
380 | // Calculate and return pdf ratio. | |
381 | ||
382 | double History::getPDFratio( int side, bool forSudakov, | |
383 | int flavNum, double xNum, double muNum, | |
384 | int flavDen, double xDen, double muDen) { | |
385 | ||
386 | // Do nothing for e+e- beams | |
387 | if ( abs(flavNum) > 10 && flavNum != 21 ) return 1.0; | |
388 | if ( abs(flavDen) > 10 && flavDen != 21 ) return 1.0; | |
389 | ||
390 | // Now calculate PDF ratio if necessary | |
391 | double pdfRatio = 1.0; | |
392 | ||
393 | // Get mother and daughter pdfs | |
394 | double pdfNum = 0.0; | |
395 | double pdfDen = 0.0; | |
396 | ||
397 | // Use rescaled PDFs in the presence of multiparton interactions | |
398 | if (side == 1) { | |
399 | if (forSudakov) | |
400 | pdfNum = mother->beamA.xfISR(0, flavNum, xNum, muNum*muNum); | |
401 | else | |
402 | pdfNum = beamA.xfISR(0, flavNum, xNum, muNum*muNum); | |
403 | if (forSudakov) | |
404 | pdfDen = max(1e-10, beamA.xfISR(0, flavDen, xDen, muDen*muDen)); | |
405 | else | |
406 | pdfDen = max(1e-10, beamA.xfISR(0, flavDen, xDen, muDen*muDen)); | |
407 | ||
408 | } else { | |
409 | if (forSudakov) | |
410 | pdfNum = mother->beamB.xfISR(0, flavNum, xNum, muNum*muNum); | |
411 | else | |
412 | pdfNum = beamB.xfISR(0, flavNum, xNum, muNum*muNum); | |
413 | ||
414 | if (forSudakov) | |
415 | pdfDen = max(1e-10,beamB.xfISR(0, flavDen, xDen, muDen*muDen)); | |
416 | else | |
417 | pdfDen = max(1e-10,beamB.xfISR(0, flavDen, xDen, muDen*muDen)); | |
418 | } | |
419 | ||
420 | // Return ratio of pdfs | |
421 | if ( pdfNum > 1e-15 && pdfDen > 1e-10 ) { | |
422 | pdfRatio *= pdfNum / pdfDen; | |
423 | ||
424 | } else if ( pdfNum < pdfDen ) { | |
425 | pdfRatio = 0.; | |
426 | } else if ( pdfNum > pdfDen ) { | |
427 | pdfRatio = 1.; | |
428 | } | |
429 | ||
430 | // Done | |
431 | return pdfRatio; | |
432 | ||
433 | } | |
434 | ||
435 | //-------------------------------------------------------------------------- | |
436 | ||
437 | /*--------------- METHODS USED FOR ONLY ONE PATH OF HISTORY NODES ------- */ | |
438 | ||
439 | // Function to set all scales in the sequence of states. This is a | |
440 | // wrapper routine for setScales and setEventScales methods | |
441 | ||
442 | void History::setScalesInHistory() { | |
443 | // Find correct links from n+1 to n states (mother --> child), as | |
444 | // needed for enforcing ordered scale sequences | |
445 | vector<int> ident; | |
446 | findPath(ident); | |
447 | // Set production scales in the states to the scales pythia would | |
448 | // have set and enforce ordering | |
449 | setScales(ident,true); | |
450 | // Set the overall event scales to the scale of the last branching | |
451 | setEventScales(); | |
452 | } | |
453 | ||
454 | //-------------------------------------------------------------------------- | |
455 | ||
456 | // Function to find the index (in the mother histories) of the | |
457 | // child history, thus providing a way access the path from both | |
458 | // initial history (mother == 0) and final history (all children == 0) | |
459 | // IN vector<int> : The index of each child in the children vector | |
460 | // of the current history node will be saved in | |
461 | // this vector | |
462 | // NO OUTPUT | |
463 | ||
464 | void History::findPath(vector<int>& out) { | |
465 | ||
466 | // If the initial and final nodes are identical, return | |
467 | if (!mother && int(children.size()) < 1) return; | |
468 | // Find the child by checking the children vector for the perfomed | |
469 | // clustering | |
470 | int iChild=-1; | |
471 | if ( mother ) { | |
472 | int size = int(mother->children.size()); | |
473 | // Loop through children and identify child chosen | |
474 | for ( int i=0; i < size; ++i) { | |
475 | if ( mother->children[i]->scale == scale | |
476 | && mother->children[i]->prob == prob | |
477 | && equalClustering(mother->children[i]->clusterIn,clusterIn)) { | |
478 | iChild = i; | |
479 | break; | |
480 | } | |
481 | } | |
482 | // Save the index of the child in the children vector and recurse | |
483 | if (iChild >-1) | |
484 | out.push_back(iChild); | |
485 | mother->findPath(out); | |
486 | } | |
487 | } | |
488 | ||
489 | //-------------------------------------------------------------------------- | |
490 | ||
491 | // Functions to set the parton production scales and enforce | |
492 | // ordering on the scales of the respective clusterings stored in | |
493 | // the History node: | |
494 | // Method will start from lowest multiplicity state and move to | |
495 | // higher states, setting the production scales the shower would | |
496 | // have used. | |
497 | // When arriving at the highest multiplicity, the method will switch | |
498 | // and go back in direction of lower states to check and enforce | |
499 | // ordering for unordered histories. | |
500 | // IN vector<int> : Vector of positions of the chosen child | |
501 | // in the mother history to allow to move | |
502 | // in direction initial->final along path | |
503 | // bool : True: Move in direction low->high | |
504 | // multiplicity and set production scales | |
505 | // False: Move in direction high->low | |
506 | // multiplicity and check and enforce | |
507 | // ordering | |
508 | // NO OUTPUT | |
509 | ||
510 | void History::setScales( vector<int> index, bool forward) { | |
511 | ||
512 | // First, set the scales of the hard process to the kinematial | |
513 | // limit (=s) | |
514 | if ( children.empty() && forward ) { | |
515 | // New "incomplete" configurations showered from mu | |
516 | if (!mother) { | |
517 | double scaleNew = 1.; | |
518 | ||
519 | if (mergingHooksPtr->incompleteScalePrescip()==0) { | |
520 | scaleNew = infoPtr->QFac(); | |
521 | } else if (mergingHooksPtr->incompleteScalePrescip()==1) { | |
522 | Vec4 pOut; | |
523 | pOut.p(0.,0.,0.,0.); | |
524 | for(int i=0; i<int(state.size()); ++i) | |
525 | if (state[i].isFinal()) | |
526 | pOut += state[i].p(); | |
527 | scaleNew = pOut.mCalc(); | |
528 | } else if (mergingHooksPtr->incompleteScalePrescip()==2) { | |
529 | scaleNew = state[0].e(); | |
530 | } | |
531 | ||
532 | scaleNew = max( mergingHooksPtr->pTcut(), scaleNew); | |
533 | ||
534 | state.scale(scaleNew); | |
535 | for(int i=3; i < int(state.size());++i) | |
536 | if (state[i].colType() != 0) | |
537 | state[i].scale(scaleNew); | |
538 | } else { | |
539 | // 2->2 with non-parton particles showered from eCM | |
540 | state.scale( state[0].e() ); | |
541 | // Count final partons | |
542 | bool isLEP = ( state[3].isLepton() && state[4].isLepton() ); | |
543 | bool isQCD = true; | |
544 | vector<int> finalPartons; | |
545 | for(int i=0; i < int(state.size());++i) { | |
546 | if (state[i].isFinal() && state[i].colType() != 0) { | |
547 | finalPartons.push_back(i); | |
548 | } | |
549 | if (state[i].isFinal() && state[i].colType() == 0) { | |
550 | isQCD = false; | |
551 | break; | |
552 | } | |
553 | } | |
554 | // If 2->2, purely partonic, set event scale to kinematic pT | |
555 | if (!isLEP && isQCD && int(finalPartons.size()) == 2) | |
556 | state.scale(state[finalPartons[0]].pT()); | |
557 | ||
558 | } | |
559 | } | |
560 | // Set all particle production scales, starting from lowest | |
561 | // multiplicity (final) state | |
562 | if (mother && forward) { | |
563 | // When choosing splitting scale, beware of unordered splittings: | |
564 | double scaleNew = 1.; | |
565 | if (mergingHooksPtr->unorderedScalePrescip() == 0) { | |
566 | // Use larger scale as common splitting scale for mother and child | |
567 | scaleNew = max( mergingHooksPtr->pTcut(), max(scale,mother->scale)); | |
568 | } else if (mergingHooksPtr->unorderedScalePrescip() == 1) { | |
569 | // Use smaller scale as common splitting scale for mother and child | |
570 | if (scale < mother->scale) | |
571 | scaleNew *= max( mergingHooksPtr->pTcut(), min(scale,mother->scale)); | |
572 | else | |
573 | scaleNew *= max( mergingHooksPtr->pTcut(), max(scale,mother->scale)); | |
574 | } | |
575 | ||
576 | // Rescale the mother state partons to the clustering scales | |
577 | // that have been found along the path | |
578 | mother->state[clusterIn.emitted].scale(scaleNew); | |
579 | mother->state[clusterIn.emittor].scale(scaleNew); | |
580 | mother->state[clusterIn.recoiler].scale(scaleNew); | |
581 | ||
582 | // Find unchanged copies of partons in higher multiplicity states | |
583 | // and rescale those | |
584 | mother->scaleCopies(clusterIn.emitted, mother->state, scaleNew); | |
585 | mother->scaleCopies(clusterIn.emittor, mother->state, scaleNew); | |
586 | mother->scaleCopies(clusterIn.recoiler, mother->state, scaleNew); | |
587 | ||
588 | // Recurse | |
589 | mother->setScales(index,true); | |
590 | } | |
591 | ||
592 | // Now, check and correct ordering from the highest multiplicity | |
593 | // state backwards to all the clustered states | |
594 | if (!mother || !forward) { | |
595 | // Get index of child along the path | |
596 | int iChild = -1; | |
597 | if ( int(index.size()) > 0 ) { | |
598 | iChild = index.back(); | |
599 | index.pop_back(); | |
600 | } | |
601 | ||
602 | // Check that the reclustered scale is above the shower cut | |
603 | if (mother) { | |
604 | scale = max(mergingHooksPtr->pTcut(), scale); | |
605 | } | |
606 | // If this is NOT the 2->2 process, check and enforce ordering | |
607 | if (iChild != -1 && !children.empty()) { | |
608 | ||
609 | if (scale > children[iChild]->scale ) { | |
610 | if (mergingHooksPtr->unorderedScalePrescip() == 0) { | |
611 | // Use larger scale as common splitting scale for mother and child | |
612 | double scaleNew = max( mergingHooksPtr->pTcut(), | |
613 | max(scale,children[iChild]->scale)); | |
614 | // Enforce ordering in particle production scales | |
615 | for( int i = 0; i < int(children[iChild]->state.size()); ++i) | |
616 | if (children[iChild]->state[i].scale() == children[iChild]->scale) | |
617 | children[iChild]->state[i].scale(scaleNew); | |
618 | // Enforce ordering in saved clustering scale | |
619 | children[iChild]->scale = scaleNew; | |
620 | ||
621 | } else if ( mergingHooksPtr->unorderedScalePrescip() == 1) { | |
622 | // Use smaller scale as common splitting scale for mother & child | |
623 | double scaleNew = max(mergingHooksPtr->pTcut(), | |
624 | min(scale,children[iChild]->scale)); | |
625 | // Enforce ordering in particle production scales | |
626 | for( int i = 0; i < int(state.size()); ++i) | |
627 | if (state[i].scale() == scale) | |
628 | state[i].scale(scaleNew); | |
629 | // Enforce ordering in saved clustering scale | |
630 | scale = scaleNew; | |
631 | } | |
632 | // Just set the overall event scale to the minimal scale | |
633 | } else { | |
634 | double scalemin = state[0].e(); | |
635 | for( int i = 0; i < int(state.size()); ++i) | |
636 | if (state[i].colType() != 0) | |
637 | scalemin = max(mergingHooksPtr->pTcut(), | |
638 | min(scalemin,state[i].scale())); | |
639 | state.scale(scalemin); | |
640 | scale = max(mergingHooksPtr->pTcut(), scale); | |
641 | } | |
642 | //Recurse | |
643 | children[iChild]->setScales(index, false); | |
644 | } | |
645 | } | |
646 | ||
647 | } | |
648 | ||
649 | //-------------------------------------------------------------------------- | |
650 | ||
651 | // Function to find a particle in all higher multiplicity events | |
652 | // along the history path and set its production scale to the input | |
653 | // scale | |
654 | // IN int iPart : Parton in refEvent to be checked / rescaled | |
655 | // Event& refEvent : Reference event for iPart | |
656 | // double scale : Scale to be set as production scale for | |
657 | // unchanged copies of iPart in subsequent steps | |
658 | ||
659 | void History::scaleCopies(int iPart, const Event& refEvent, double rho) { | |
660 | ||
661 | // Check if any parton recently rescaled is found unchanged: | |
662 | // Same charge, colours in mother->state | |
663 | if ( mother ) { | |
664 | for( int i=0; i < mother->state.size(); ++i) { | |
665 | if ( ( mother->state[i].id() == refEvent[iPart].id() | |
666 | && mother->state[i].colType() == refEvent[iPart].colType() | |
667 | && mother->state[i].chargeType() == refEvent[iPart].chargeType() | |
668 | && mother->state[i].col() == refEvent[iPart].col() | |
669 | && mother->state[i].acol() == refEvent[iPart].acol() ) | |
670 | ) { | |
671 | // Rescale the unchanged parton | |
672 | mother->state[i].scale(rho); | |
673 | // Recurse | |
674 | if (mother->mother) | |
675 | mother->scaleCopies( iPart, refEvent, rho ); | |
676 | } // end if found unchanged parton case | |
677 | } // end loop over particle entries in event | |
678 | } | |
679 | } | |
680 | ||
681 | //-------------------------------------------------------------------------- | |
682 | ||
683 | // Functions to set the OVERALL EVENT SCALES [=state.scale()] to | |
684 | // the scale of the last clustering | |
685 | // NO INPUT | |
686 | // NO OUTPUT | |
687 | ||
688 | void History::setEventScales() { | |
689 | // Set the event scale to the scale of the last clustering, | |
690 | // except for the very lowest multiplicity state | |
691 | if (mother) { | |
692 | mother->state.scale(scale); | |
693 | // Recurse | |
694 | mother->setEventScales(); | |
695 | } | |
696 | } | |
697 | ||
698 | //-------------------------------------------------------------------------- | |
699 | ||
700 | // Functions to return the z value of the last ISR splitting | |
701 | // NO INPUT | |
702 | // OUTPUT double : z value of last ISR splitting in history | |
703 | ||
704 | double History::zISR() { | |
705 | ||
706 | // Do nothing for ME level state | |
707 | if (!mother) return 0.0; | |
708 | // Skip FSR splitting | |
709 | if (mother->state[clusterIn.emittor].isFinal()) return mother->zISR(); | |
710 | // Calculate z | |
711 | int rad = clusterIn.emittor; | |
712 | int rec = clusterIn.recoiler; | |
713 | int emt = clusterIn.emitted; | |
714 | double z = (mother->state[rad].p() + mother->state[rec].p() | |
715 | - mother->state[emt].p()).m2Calc() | |
716 | / (mother->state[rad].p() + mother->state[rec].p()).m2Calc(); | |
717 | // Recurse | |
718 | double znew = mother->zISR(); | |
719 | // Update z | |
720 | if (znew > 0.) z = znew; | |
721 | ||
722 | return z; | |
723 | } | |
724 | ||
725 | //-------------------------------------------------------------------------- | |
726 | ||
727 | // Functions to return the z value of the last FSR splitting | |
728 | // NO INPUT | |
729 | // OUTPUT double : z value of last FSR splitting in history | |
730 | ||
731 | double History::zFSR() { | |
732 | ||
733 | // Do nothing for ME level state | |
734 | if (!mother) return 0.0; | |
735 | // Skip ISR splitting | |
736 | if (!mother->state[clusterIn.emittor].isFinal()) return mother->zFSR(); | |
737 | // Calculate z | |
738 | int rad = clusterIn.emittor; | |
739 | int rec = clusterIn.recoiler; | |
740 | int emt = clusterIn.emitted; | |
741 | // Construct 2->3 variables for FSR | |
742 | Vec4 sum = mother->state[rad].p() + mother->state[rec].p() | |
743 | + mother->state[emt].p(); | |
744 | double m2Dip = sum.m2Calc(); | |
745 | double x1 = 2. * (sum * mother->state[rad].p()) / m2Dip; | |
746 | double x3 = 2. * (sum * mother->state[emt].p()) / m2Dip; | |
747 | // Calculate z of splitting for FSR | |
748 | double z = x1/(x1+x3); | |
749 | // Recurse | |
750 | double znew = mother->zFSR(); | |
751 | // Update z | |
752 | if (znew > 0.) z = znew; | |
753 | ||
754 | return z; | |
755 | } | |
756 | ||
757 | //-------------------------------------------------------------------------- | |
758 | ||
759 | // Functions to return the pT scale of the last FSR splitting | |
760 | // NO INPUT | |
761 | // OUTPUT double : pT scale of last FSR splitting in history | |
762 | ||
763 | double History::pTISR() { | |
764 | // Do nothing for ME level state | |
765 | if (!mother) return 0.0; | |
766 | // Skip FSR splitting | |
767 | if (mother->state[clusterIn.emittor].isFinal()) return mother->pTISR(); | |
768 | double pT = mother->state.scale(); | |
769 | // Recurse | |
770 | double pTnew = mother->pTISR(); | |
771 | // Update pT | |
772 | if (pTnew > 0.) pT = pTnew; | |
773 | ||
774 | return pT; | |
775 | } | |
776 | ||
777 | //-------------------------------------------------------------------------- | |
778 | ||
779 | // Functions to return the pT scale of the last FSR splitting | |
780 | // NO INPUT | |
781 | // OUTPUT double : pT scale of last FSR splitting in history | |
782 | ||
783 | double History::pTFSR() { | |
784 | // Do nothing for ME level state | |
785 | if (!mother) return 0.0; | |
786 | // Skip ISR splitting | |
787 | if (!mother->state[clusterIn.emittor].isFinal()) return mother->pTFSR(); | |
788 | double pT = mother->state.scale(); | |
789 | // Recurse | |
790 | double pTnew = mother->pTFSR(); | |
791 | // Update pT | |
792 | if (pTnew > 0.) pT = pTnew; | |
793 | return pT; | |
794 | } | |
795 | ||
796 | //-------------------------------------------------------------------------- | |
797 | ||
798 | // Functions to return the depth of the history (i.e. the number of | |
799 | // reclustered splittings) | |
800 | // NO INPUT | |
801 | // OUTPUT int : Depth of history | |
802 | ||
803 | int History::nClusterings() { | |
804 | ||
805 | if (!mother) return 0; | |
806 | int w = mother->nClusterings(); | |
807 | w += 1; | |
808 | return w; | |
809 | } | |
810 | ||
811 | //-------------------------------------------------------------------------- | |
812 | ||
813 | // Functions to return the event after nSteps splittings of the 2->2 process | |
814 | // Example: nSteps = 1 -> return event with one additional parton | |
815 | // INPUT int : Number of splittings in the event, | |
816 | // as counted from core 2->2 process | |
817 | // OUTPUT Event : event with nSteps additional partons | |
818 | ||
819 | Event History::clusteredState(int nSteps) { | |
820 | ||
821 | // Save state | |
822 | Event outState = state; | |
823 | // As long as there are steps to do, recursively save state | |
824 | if (mother && nSteps > 0) | |
825 | outState = mother->clusteredState(nSteps - 1); | |
826 | // Done | |
827 | return outState; | |
828 | ||
829 | } | |
830 | ||
831 | //-------------------------------------------------------------------------- | |
832 | ||
833 | // Function to choose a path from all paths in the tree | |
834 | // according to their splitting probabilities | |
835 | // IN double : Random number | |
836 | // OUT History* : Leaf of history path chosen | |
837 | ||
838 | History * History::select(double rnd) { | |
839 | ||
840 | // No need to choose if no paths have been constructed. | |
841 | if ( goodBranches.empty() && badBranches.empty() ) return this; | |
842 | ||
843 | // Choose amongst paths allowed by projections. | |
844 | double sum = 0.; | |
845 | map<double, History*> selectFrom; | |
846 | if ( !goodBranches.empty() ) { | |
847 | selectFrom = goodBranches; | |
848 | sum = sumGoodBranches; | |
849 | } else { | |
850 | selectFrom = badBranches; | |
851 | sum = sumBadBranches; | |
852 | } | |
853 | ||
854 | if (mergingHooksPtr->pickBySumPT()) { | |
855 | // Find index of history with minimal sum of scalar pT | |
856 | int nFinal = 0; | |
857 | for (int i=0; i < state.size(); ++i) | |
858 | if (state[i].isFinal()) | |
859 | nFinal++; | |
860 | double iMin = 0.; | |
861 | double sumMin = (nFinal-2)*state[0].e(); | |
862 | for ( map<double, History*>::iterator it = selectFrom.begin(); | |
863 | it != selectFrom.end(); ++it ) { | |
864 | ||
865 | if (it->second->sumScalarPT < sumMin) { | |
866 | sumMin = it->second->sumScalarPT; | |
867 | iMin = it->first; | |
868 | } | |
869 | } | |
870 | // Choose history with smallest sum of scalar pT | |
871 | return selectFrom.lower_bound(iMin)->second; | |
872 | } else { | |
873 | // Choose history according to probability, be careful about upper bound | |
874 | if ( rnd != 1. ) { | |
875 | return selectFrom.upper_bound(sum*rnd)->second; | |
876 | } else { | |
877 | return selectFrom.lower_bound(sum*rnd)->second; | |
878 | } | |
879 | } | |
880 | // Done | |
881 | } | |
882 | ||
883 | //-------------------------------------------------------------------------- | |
884 | ||
885 | // Function to project paths onto desired paths. | |
886 | ||
887 | bool History::trimHistories() { | |
888 | // Do nothing if no paths have been constructed. | |
889 | if ( paths.empty() ) return false; | |
890 | // Loop through all constructed paths. Check all removal conditions. | |
891 | for ( map<double, History*>::iterator it = paths.begin(); | |
892 | it != paths.end(); ++it ) { | |
893 | // Check if history is allowed. | |
894 | if ( it->second->keep() && !it->second->keepHistory() ) | |
895 | it->second->remove(); | |
896 | } | |
897 | // Project onto desired / undesired branches. | |
898 | double sumold, sumnew, sumprob, mismatch; | |
899 | sumold = sumnew = sumprob = mismatch = 0.; | |
900 | // Loop through all constructed paths and store allowed paths. | |
901 | // Skip undesired paths. | |
902 | for ( map<double, History*>::iterator it = paths.begin(); | |
903 | it != paths.end(); ++it ) { | |
904 | // Update index | |
905 | sumnew = it->first; | |
906 | if ( it->second->keep() ) { | |
907 | // Fill branches with allowed paths. | |
908 | goodBranches.insert( make_pair( sumnew - mismatch, it->second) ); | |
909 | // Add probability of this path. | |
910 | sumGoodBranches = sumnew - mismatch; | |
911 | } else { | |
912 | // Update mismatch in probabilities resulting from not including this | |
913 | // path | |
914 | double mismatchOld = mismatch; | |
915 | mismatch += sumnew - sumold; | |
916 | // Fill branches with allowed paths. | |
917 | badBranches.insert( make_pair( mismatchOld + sumnew - sumold, | |
918 | it->second ) ); | |
919 | // Add probability of this path. | |
920 | sumBadBranches = mismatchOld + sumnew - sumold; | |
921 | } | |
922 | // remember index of this path in order to caclulate probability of | |
923 | // subsequent path. | |
924 | sumold = it->first; | |
925 | } | |
926 | // Done | |
927 | return !goodBranches.empty(); | |
928 | } | |
929 | ||
930 | //-------------------------------------------------------------------------- | |
931 | ||
932 | // Function implementing checks on a paths, for deciding if the path should | |
933 | // be considered valid or not. | |
934 | ||
935 | bool History::keepHistory() { | |
936 | bool keepPath = true; | |
937 | //// Tag unordered paths for removal. | |
938 | //double maxScale = (foundCompletePath) | |
939 | // ? infoPtr->eCM() | |
940 | // : infoPtr->QFac(); | |
941 | //keepPath = isOrderedPath( maxScale ); | |
942 | //Done | |
943 | return keepPath; | |
944 | } | |
945 | ||
946 | //-------------------------------------------------------------------------- | |
947 | ||
948 | // Function to check if a path is ordered in evolution pT. | |
949 | ||
950 | bool History::isOrderedPath( double maxscale ) { | |
951 | double newscale = clusterIn.pT(); | |
952 | if ( !mother ) return true; | |
953 | bool ordered = mother->isOrderedPath(newscale); | |
954 | if ( !ordered || maxscale < newscale) return false; | |
955 | return ordered; | |
956 | } | |
957 | ||
958 | //-------------------------------------------------------------------------- | |
959 | ||
960 | // Function to check if all reconstucted states in a path pass the merging | |
961 | // scale cut. | |
962 | ||
963 | bool History::allIntermediateAboveRhoMS( double rhoms ) { | |
964 | ||
965 | int nFinal = 0; | |
966 | for ( int i = 0; i < state.size(); ++i ) | |
967 | if ( state[i].isFinal() && state[i].colType() != 0 ) | |
968 | nFinal++; | |
969 | ||
970 | double rhoNew = mergingHooksPtr->rhoms( state, false ); | |
971 | ||
972 | if ( !mother ) return true; | |
973 | ||
974 | bool passRHOMS = mother->allIntermediateAboveRhoMS( rhoms ); | |
975 | ||
976 | if ( !passRHOMS || ( nFinal > 0 && rhoNew < rhoms ) ) return false; | |
977 | ||
978 | // Done | |
979 | return passRHOMS; | |
980 | } | |
981 | ||
982 | //-------------------------------------------------------------------------- | |
983 | ||
984 | // Function to check if reconstucted states in a path are ordered in the | |
985 | // merging scale variable. | |
986 | ||
987 | bool History::intermediateRhoMSOrdered( double maxscale ) { | |
988 | // Count number of final state partons. | |
989 | int nFinal = 0; | |
990 | for ( int i = 0; i < state.size(); ++i ) | |
991 | if ( state[i].isFinal() && state[i].colType() != 0 ) | |
992 | nFinal++; | |
993 | // Get new min(rho) scale. | |
994 | double newscale = (nFinal == 0 ) | |
995 | ? maxscale | |
996 | : mergingHooksPtr->rhoms( state, false ); | |
997 | // For no final state particles or available mother, check mother. | |
998 | bool isOrdered = ( nFinal == 0 || mother ) | |
999 | ? mother->intermediateRhoMSOrdered( newscale ) | |
1000 | : true; | |
1001 | // Done. | |
1002 | return isOrdered && ( maxscale >= newscale ); | |
1003 | } | |
1004 | ||
1005 | //-------------------------------------------------------------------------- | |
1006 | ||
1007 | // Function to check if any ordered paths were found (and kept). | |
1008 | ||
1009 | bool History::foundAnyOrderedPaths() { | |
1010 | //Do nothing if no paths were found | |
1011 | if ( paths.empty() ) return false; | |
1012 | double maxscale = infoPtr->eCM(); | |
1013 | // Loop through paths. Divide probability into ordered and unordered pieces. | |
1014 | for ( map<double, History*>::iterator it = paths.begin(); | |
1015 | it != paths.end(); ++it ) | |
1016 | if ( it->second->isOrderedPath(maxscale) ) | |
1017 | return true; | |
1018 | // Done | |
1019 | return false; | |
1020 | } | |
1021 | ||
1022 | //-------------------------------------------------------------------------- | |
1023 | ||
1024 | // Function to check if any unordered paths were found (and kept). | |
1025 | ||
1026 | bool History::foundAnyUnorderedPaths() { | |
1027 | // Do nothing if no paths were found | |
1028 | if ( paths.empty() ) return false; | |
1029 | double maxscale = infoPtr->eCM(); | |
1030 | // Loop through paths. Return false if an ordered path was found. | |
1031 | for ( map<double, History*>::iterator it = paths.begin(); | |
1032 | it != paths.end(); ++it ) | |
1033 | if ( !it->second->isOrderedPath(maxscale) ) | |
1034 | return true; | |
1035 | // Done | |
1036 | return false; | |
1037 | } | |
1038 | ||
1039 | //-------------------------------------------------------------------------- | |
1040 | ||
1041 | // For a full path, find the weight calculated from the ratio of | |
1042 | // couplings, the no-emission probabilities, and possible PDF | |
1043 | // ratios. This function should only be called for the last history | |
1044 | // node of a full path. | |
1045 | // IN TimeShower : Already initialised shower object to be used as | |
1046 | // trial shower | |
1047 | // double : alpha_s value used in ME calculation | |
1048 | // double : Maximal mass scale of the problem (e.g. E_CM) | |
1049 | // AlphaStrong: Initialised shower alpha_s object for FSR | |
1050 | // alpha_s ratio calculation | |
1051 | // AlphaStrong: Initialised shower alpha_s object for ISR | |
1052 | // alpha_s ratio calculation (can be different from previous) | |
1053 | ||
1054 | double History::weightTree(PartonLevel* trial, double as0, double maxscale, | |
1055 | double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR, | |
1056 | double& asWeight, double& pdfWeight) { | |
1057 | ||
1058 | // Use correct scale | |
1059 | double newScale = scale; | |
1060 | ||
1061 | // For ME state, just multiply by PDF ratios | |
1062 | if ( !mother ) { | |
1063 | ||
1064 | int sideRad = (state[3].pz() > 0) ? 1 :-1; | |
1065 | int sideRec = (state[4].pz() > 0) ? 1 :-1; | |
1066 | // Calculate PDF first leg | |
1067 | if (state[3].colType() != 0) { | |
1068 | // Find x value and flavour | |
1069 | double x = 2.*state[3].e() / state[0].e(); | |
1070 | int flav = state[3].id(); | |
1071 | ||
1072 | // Find numerator/denominator scale | |
1073 | double scaleNum = (children.empty()) ? hardFacScale() : maxscale; | |
1074 | double scaleDen = infoPtr->QFac(); | |
1075 | // For initial parton, multiply by PDF ratio | |
1076 | double ratio = getPDFratio(sideRad, false, flav, x, scaleNum, | |
1077 | flav, x, scaleDen); | |
1078 | pdfWeight *= ratio; | |
1079 | } | |
1080 | ||
1081 | // Calculate PDF ratio for second leg | |
1082 | if (state[4].colType() != 0) { | |
1083 | // Find x value and flavour | |
1084 | double x = 2.*state[4].e() / state[0].e(); | |
1085 | int flav = state[4].id(); | |
1086 | // Find numerator/denominator scale | |
1087 | double scaleNum = (children.empty()) ? hardFacScale() : maxscale; | |
1088 | double scaleDen = infoPtr->QFac(); | |
1089 | // For initial parton, multiply with PDF ratio | |
1090 | double ratio = getPDFratio(sideRec, false, flav, x, scaleNum, | |
1091 | flav, x, scaleDen); | |
1092 | pdfWeight *= ratio; | |
1093 | } | |
1094 | ||
1095 | return 1.0; | |
1096 | } | |
1097 | ||
1098 | // Remember new PDF scale n case true sclae should be used for un-ordered | |
1099 | // splittings. | |
1100 | double newPDFscale = newScale; | |
1101 | if (mergingHooksPtr->unorderedPDFscalePrescip() == 1) | |
1102 | newPDFscale = clusterIn.pT(); | |
1103 | ||
1104 | // Recurse | |
1105 | double w = mother->weightTree(trial, as0, newScale, newPDFscale, | |
1106 | asFSR, asISR, asWeight, pdfWeight); | |
1107 | ||
1108 | // Do nothing for empty state | |
1109 | if (state.size() < 3) return 1.0; | |
1110 | // If up to now, trial shower was not successful, return zero | |
1111 | if ( w < 1e-12 ) return 0.0; | |
1112 | // Do trial shower on current state, return zero if not successful | |
1113 | w *= doTrialShower(trial, maxscale); | |
1114 | if ( w < 1e-12 ) return 0.0; | |
1115 | ||
1116 | // Calculate alpha_s ratio for current state | |
1117 | if ( asFSR && asISR ) { | |
1118 | double asScale = pow2( newScale ); | |
1119 | if (mergingHooksPtr->unorderedASscalePrescip() == 1) | |
1120 | asScale = pow2( clusterIn.pT() ); | |
1121 | bool FSR = mother->state[clusterIn.emittor].isFinal(); | |
1122 | double alphaSinPS = (FSR) ? (*asFSR).alphaS(asScale) | |
1123 | : (*asISR).alphaS(asScale | |
1124 | + pow2(mergingHooksPtr->pT0ISR()) ); | |
1125 | asWeight *= alphaSinPS / as0; | |
1126 | } | |
1127 | ||
1128 | // Calculate pdf ratios: Get both sides of event | |
1129 | int inP = 3; | |
1130 | int inM = 4; | |
1131 | int sideP = (mother->state[inP].pz() > 0) ? 1 :-1; | |
1132 | int sideM = (mother->state[inM].pz() > 0) ? 1 :-1; | |
1133 | ||
1134 | if ( mother->state[inP].colType() != 0 ) { | |
1135 | // Find x value and flavour | |
1136 | double x = getCurrentX(sideP); | |
1137 | int flav = getCurrentFlav(sideP); | |
1138 | // Find numerator scale | |
1139 | double scaleNum = (children.empty()) | |
1140 | ? hardFacScale() | |
1141 | : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1) | |
1142 | ? pdfScale : maxscale ); | |
1143 | double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1) | |
1144 | ? clusterIn.pT() : newScale; | |
1145 | // Multiply PDF ratio | |
1146 | double ratio = getPDFratio(sideP, false, flav, x, scaleNum, | |
1147 | flav, x, scaleDen); | |
1148 | pdfWeight *= ratio; | |
1149 | } | |
1150 | ||
1151 | if ( mother->state[inM].colType() != 0 ) { | |
1152 | // Find x value and flavour | |
1153 | double x = getCurrentX(sideM); | |
1154 | int flav = getCurrentFlav(sideM); | |
1155 | // Find numerator scale | |
1156 | double scaleNum = (children.empty()) | |
1157 | ? hardFacScale() | |
1158 | : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1) | |
1159 | ? pdfScale : maxscale ); | |
1160 | double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1) | |
1161 | ? clusterIn.pT() : newScale; | |
1162 | // Multiply PDF ratio | |
1163 | double ratio = getPDFratio(sideM, false, flav, x, scaleNum, | |
1164 | flav, x, scaleDen); | |
1165 | pdfWeight *= ratio; | |
1166 | } | |
1167 | ||
1168 | // Done | |
1169 | return w; | |
1170 | } | |
1171 | ||
1172 | //-------------------------------------------------------------------------- | |
1173 | ||
1174 | // Function to return the factorisation scale of the hard process in Pythia. | |
1175 | ||
1176 | double History::hardFacScale() { | |
1177 | // Declare output scale. | |
1178 | double hardscale = 0.; | |
1179 | // If scale should not be reset, done. | |
1180 | if ( !mergingHooksPtr->resetHardQFac() ) return infoPtr->QFac(); | |
1181 | // For pure QCD dijet events, calculate the hadronic cross section | |
1182 | // of the hard process at the pT of the dijet system, rather than at fixed | |
1183 | // arbitrary scale. | |
1184 | if ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0 | |
1185 | || mergingHooksPtr->getProcessString().compare("pp>aj") == 0 ) { | |
1186 | // Find the mT in the hard sub-process. | |
1187 | vector <double> mT; | |
1188 | for ( int i=0; i < state.size(); ++i) | |
1189 | if ( state[i].isFinal() && state[i].colType() != 0 ) | |
1190 | mT.push_back( abs(state[i].mT2()) ); | |
1191 | if ( int(mT.size()) != 2 ) | |
1192 | hardscale = infoPtr->QFac(); | |
1193 | else | |
1194 | hardscale = sqrt( min( mT[0], mT[1] ) ); | |
1195 | } else { | |
1196 | hardscale = infoPtr->QFac(); | |
1197 | } | |
1198 | // Done | |
1199 | return hardscale; | |
1200 | } | |
1201 | ||
1202 | //-------------------------------------------------------------------------- | |
1203 | ||
1204 | // Function to return the factorisation scale of the hard process in Pythia. | |
1205 | ||
1206 | double History::hardRenScale() { | |
1207 | // Declare output scale. | |
1208 | double hardscale = 0.; | |
1209 | // For pure QCD dijet events, calculate the hadronic cross section | |
1210 | // of the hard process at the pT of the dijet system, rather than at fixed | |
1211 | // arbitrary scale. | |
1212 | if ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0 | |
1213 | || mergingHooksPtr->getProcessString().compare("pp>aj") == 0 ) { | |
1214 | // Find the mT in the hard sub-process. | |
1215 | vector <double> mT; | |
1216 | for ( int i=0; i < state.size(); ++i) | |
1217 | if ( state[i].isFinal() | |
1218 | && ( state[i].colType() != 0 || state[i].id() == 22 ) ) | |
1219 | mT.push_back( abs(state[i].mT()) ); | |
1220 | if ( int(mT.size()) != 2 ) | |
1221 | hardscale = infoPtr->QRen(); | |
1222 | else | |
1223 | hardscale = sqrt( mT[0]*mT[1] ); | |
1224 | } else { | |
1225 | hardscale = infoPtr->QRen(); | |
1226 | } | |
1227 | // Done | |
1228 | return hardscale; | |
1229 | } | |
1230 | ||
1231 | //-------------------------------------------------------------------------- | |
1232 | ||
1233 | // Perform a trial shower using the \a pythia object between | |
1234 | // maxscale down to this scale and return the corresponding Sudakov | |
1235 | // form factor. | |
1236 | // IN trialShower : Shower object used as trial shower | |
1237 | // double : Maximum scale for trial shower branching | |
1238 | // OUT 0.0 : trial shower emission outside allowed pT range | |
1239 | // 1.0 : trial shower successful (any emission was below | |
1240 | // the minimal scale ) | |
1241 | ||
1242 | double History::doTrialShower( PartonLevel* trial, double maxscale, | |
1243 | double minscale ) { | |
1244 | ||
1245 | // Copy state to local process | |
1246 | Event process = state; | |
1247 | // Set starting scale. | |
1248 | double startingScale = maxscale; | |
1249 | // Set output. | |
1250 | bool doVeto = false; | |
1251 | ||
1252 | while ( true ) { | |
1253 | ||
1254 | // Reset trialShower object | |
1255 | trial->resetTrial(); | |
1256 | // Construct event to be showered | |
1257 | Event event = Event(); | |
1258 | event.init("(hard process-modified)", particleDataPtr); | |
1259 | event.clear(); | |
1260 | // Reset process scale so thatshower starting scale is correctly set. | |
1261 | process.scale(startingScale); | |
1262 | ||
1263 | doVeto = false; | |
1264 | ||
1265 | // Get pT before reclustering | |
1266 | double minScale = (minscale > 0.) ? minscale : scale; | |
1267 | // If the maximal scale and the minimal scale coincide (as would | |
1268 | // be the case for the corrected scales of unordered histories), | |
1269 | // do not generate Sudakov | |
1270 | if (minScale >= startingScale) break; | |
1271 | ||
1272 | // Find z and pT values at which the current state was formed, to | |
1273 | // ensure that the showers can order the next emission correctly in | |
1274 | // rapidity, if required. | |
1275 | // NOT CORRECTLY SET FOR HIGHEST MULTIPLICITY STATE! | |
1276 | double z = ( mergingHooksPtr->getNumberOfClusteringSteps(state) == 0 | |
1277 | || !mother ) | |
1278 | ? 0.5 | |
1279 | : mother->getCurrentZ(clusterIn.emittor,clusterIn.recoiler, | |
1280 | clusterIn.emitted); | |
1281 | // Store z and pT values at which the current state was formed | |
1282 | infoPtr->zNowISR(z); | |
1283 | infoPtr->pT2NowISR(pow(startingScale,2)); | |
1284 | infoPtr->hasHistory(true); | |
1285 | ||
1286 | // Perform trial shower emission | |
1287 | trial->next(process,event); | |
1288 | // Get trial shower pT. | |
1289 | double pTtrial = trial->pTLastInShower(); | |
1290 | double typeTrial = trial->typeLastInShower(); | |
1291 | ||
1292 | // Get veto (merging) scale value | |
1293 | double vetoScale = (mother) ? 0. : mergingHooksPtr->tms(); | |
1294 | // Get merging scale in current event | |
1295 | double tnow = 0.; | |
1296 | // Use KT/Durham merging scale definition. | |
1297 | if (mergingHooksPtr->doKTMerging() || mergingHooksPtr->doMGMerging()) | |
1298 | tnow = mergingHooksPtr->kTms(event); | |
1299 | // Use Lund PT merging scale definition. | |
1300 | else if (mergingHooksPtr->doPTLundMerging()) | |
1301 | tnow = mergingHooksPtr->rhoms(event, false); | |
1302 | // Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale definition. | |
1303 | else if (mergingHooksPtr->doCutBasedMerging()) | |
1304 | tnow = mergingHooksPtr->cutbasedms(event); | |
1305 | // Use user-defined merging scale. | |
1306 | else | |
1307 | tnow = mergingHooksPtr->tmsDefinition(event); | |
1308 | ||
1309 | // Done if evolution scale has fallen below minimum | |
1310 | if ( pTtrial < minScale ) break; | |
1311 | // Reset starting scale. | |
1312 | startingScale = pTtrial; | |
1313 | ||
1314 | // Continue if this state is below the veto scale | |
1315 | if ( tnow < vetoScale && vetoScale > 0. ) continue; | |
1316 | ||
1317 | // Retry if the trial emission was not allowed. | |
1318 | if ( mergingHooksPtr->canVetoTrialEmission() | |
1319 | && mergingHooksPtr->doVetoTrialEmission( process, event) ) continue; | |
1320 | ||
1321 | // Veto event if trial pT was above the next nodal scale. | |
1322 | if ( pTtrial > minScale ) doVeto = true; | |
1323 | ||
1324 | //// For last no-emission probability, veto event if the trial state | |
1325 | //// is above the merging scale, i.e. in the matrix element region. | |
1326 | //if ( !mother && tnow > vetoScale && vetoScale > 0. ) doVeto = true; | |
1327 | ||
1328 | // For 2 -> 2 pure QCD state, do not allow multiparton interactions | |
1329 | // above the kinematical pT of the 2 -> 2 state | |
1330 | if (typeTrial == 1) { | |
1331 | // Count number of final state particles and remember partons | |
1332 | int nFinal = 0; | |
1333 | vector<int> finalPartons; | |
1334 | for(int i=0; i < state.size(); ++i) { | |
1335 | if (state[i].isFinal()) { | |
1336 | nFinal++; | |
1337 | if ( state[i].colType() != 0) | |
1338 | finalPartons.push_back(i); | |
1339 | } | |
1340 | } | |
1341 | // Veto if MPI was above 2 -> 2 pT | |
1342 | if ( nFinal == 2 && int(finalPartons.size()) == 2 | |
1343 | && pTtrial > event[finalPartons[0]].pT() ) { | |
1344 | return 0.0; | |
1345 | } | |
1346 | } | |
1347 | ||
1348 | // If pT of trial emission was in suitable range (trial shower | |
1349 | // successful), return false | |
1350 | if ( pTtrial < minScale ) doVeto = false; | |
1351 | ||
1352 | // Done | |
1353 | break; | |
1354 | ||
1355 | } | |
1356 | ||
1357 | // Done | |
1358 | return ( (doVeto) ? 0. : 1. ); | |
1359 | } | |
1360 | ||
1361 | /*--------------- METHODS USED FOR CONTRUCTION OF ALL HISTORIES --------- */ | |
1362 | ||
1363 | // Check if a ordered (and complete) path has been found in the | |
1364 | // initial node, in which case we will no longer be interested in | |
1365 | // any unordered paths. | |
1366 | ||
1367 | bool History::onlyOrderedPaths() { | |
1368 | if ( !mother || foundOrderedPath ) return foundOrderedPath; | |
1369 | return foundOrderedPath = mother->onlyOrderedPaths(); | |
1370 | } | |
1371 | ||
1372 | bool History::onlyUnorderedPaths() { | |
1373 | if ( !mother || foundUnorderedPath ) return foundUnorderedPath; | |
1374 | return foundUnorderedPath = mother->onlyUnorderedPaths(); | |
1375 | } | |
1376 | ||
1377 | //-------------------------------------------------------------------------- | |
1378 | ||
1379 | // Check if a STRONGLY ordered (and complete) path has been found in the | |
1380 | // initial node, in which case we will no longer be interested in | |
1381 | // any unordered paths. | |
1382 | ||
1383 | bool History::onlyStronglyOrderedPaths() { | |
1384 | if ( !mother || foundStronglyOrderedPath ) return foundStronglyOrderedPath; | |
1385 | return foundStronglyOrderedPath = mother->onlyStronglyOrderedPaths(); | |
1386 | } | |
1387 | ||
1388 | //-------------------------------------------------------------------------- | |
1389 | ||
1390 | // Check if an allowed (according to user-criterion) path has been found in | |
1391 | // the initial node, in which case we will no longer be interested in | |
1392 | // any forbidden paths. | |
1393 | ||
1394 | bool History::onlyAllowedPaths() { | |
1395 | if ( !mother || foundAllowedPath ) return foundAllowedPath; | |
1396 | return foundAllowedPath = mother->onlyAllowedPaths(); | |
1397 | } | |
1398 | ||
1399 | //-------------------------------------------------------------------------- | |
1400 | ||
1401 | // When a full path has been found, register it with the initial | |
1402 | // history node. | |
1403 | // IN History : History to be registered as path | |
1404 | // bool : Specifying if clusterings so far were ordered | |
1405 | // bool : Specifying if path is complete down to 2->2 process | |
1406 | // OUT true if History object forms a plausible path (eg prob>0 ...) | |
1407 | ||
1408 | bool History::registerPath(History & l, bool isOrdered, bool isUnordered, | |
1409 | bool isStronglyOrdered, bool isAllowed, bool isComplete) { | |
1410 | ||
1411 | // We are not interested in improbable paths. | |
1412 | if ( l.prob <= 0.0) | |
1413 | return false; | |
1414 | // We only register paths in the initial node. | |
1415 | if ( mother ) return mother->registerPath(l, isOrdered, isUnordered, | |
1416 | isStronglyOrdered, isAllowed, isComplete); | |
1417 | // Again, we are not interested in improbable paths. | |
1418 | if ( sumpath == sumpath + l.prob ) | |
1419 | return false; | |
1420 | if ( mergingHooksPtr->canCutOnRecState() | |
1421 | && foundAllowedPath && !isAllowed ) | |
1422 | return false; | |
1423 | if ( mergingHooksPtr->enforceStrongOrdering() | |
1424 | && foundStronglyOrderedPath && !isStronglyOrdered ) | |
1425 | return false; | |
1426 | if ( mergingHooksPtr->orderHistories() | |
1427 | && foundOrderedPath && !isOrdered ) | |
1428 | return false; | |
1429 | ||
1430 | if ( mergingHooksPtr->forceUnorderedHistories() | |
1431 | && foundUnorderedPath && !isUnordered ) | |
1432 | return false; | |
1433 | ||
1434 | if ( foundCompletePath && !isComplete ) | |
1435 | return false; | |
1436 | if ( !mergingHooksPtr->canCutOnRecState() | |
1437 | && !mergingHooksPtr->allowCutOnRecState() ) | |
1438 | foundAllowedPath = true; | |
1439 | ||
1440 | if ( mergingHooksPtr->canCutOnRecState() && isAllowed && isComplete) { | |
1441 | if ( !foundAllowedPath || !foundCompletePath ) { | |
1442 | // If this is the first complete, allowed path, discard the | |
1443 | // old, disallowed or incomplete ones. | |
1444 | paths.clear(); | |
1445 | sumpath = 0.0; | |
1446 | } | |
1447 | foundAllowedPath = true; | |
1448 | ||
1449 | } | |
1450 | ||
1451 | if ( mergingHooksPtr->enforceStrongOrdering() && isStronglyOrdered | |
1452 | && isComplete ) { | |
1453 | if ( !foundStronglyOrderedPath || !foundCompletePath ) { | |
1454 | // If this is the first complete, ordered path, discard the | |
1455 | // old, non-ordered or incomplete ones. | |
1456 | paths.clear(); | |
1457 | sumpath = 0.0; | |
1458 | } | |
1459 | foundStronglyOrderedPath = true; | |
1460 | foundCompletePath = true; | |
1461 | ||
1462 | } | |
1463 | ||
1464 | if ( mergingHooksPtr->orderHistories() && isOrdered && isComplete ) { | |
1465 | if ( !foundOrderedPath || !foundCompletePath ) { | |
1466 | // If this is the first complete, ordered path, discard the | |
1467 | // old, non-ordered or incomplete ones. | |
1468 | paths.clear(); | |
1469 | sumpath = 0.0; | |
1470 | } | |
1471 | foundOrderedPath = true; | |
1472 | foundCompletePath = true; | |
1473 | ||
1474 | } | |
1475 | ||
1476 | if ( mergingHooksPtr->forceUnorderedHistories() && isUnordered | |
1477 | && isComplete ) { | |
1478 | if ( !foundUnorderedPath || !foundCompletePath ) { | |
1479 | // If this is the first complete, unordered path, discard the | |
1480 | // old, ordered or incomplete ones. | |
1481 | paths.clear(); | |
1482 | sumpath = 0.0; | |
1483 | } | |
1484 | foundUnorderedPath = true; | |
1485 | foundCompletePath = true; | |
1486 | } | |
1487 | ||
1488 | if ( isComplete ) { | |
1489 | if ( !foundCompletePath ) { | |
1490 | // If this is the first complete path, discard the old, | |
1491 | // incomplete ones. | |
1492 | paths.clear(); | |
1493 | sumpath = 0.0; | |
1494 | } | |
1495 | foundCompletePath = true; | |
1496 | } | |
1497 | ||
1498 | // Remember, if this path is ordered, even if no ordering is required | |
1499 | if ( isOrdered ) { | |
1500 | foundOrderedPath = true; | |
1501 | } | |
1502 | ||
1503 | // Index path by probability | |
1504 | sumpath += l.prob; | |
1505 | paths[sumpath] = &l; | |
1506 | ||
1507 | return true; | |
1508 | } | |
1509 | ||
1510 | //-------------------------------------------------------------------------- | |
1511 | ||
1512 | // For the history-defining state (and if necessary interfering | |
1513 | // states), find all possible clusterings. | |
1514 | // NO INPUT | |
1515 | // OUT vector of all (rad,rec,emt) systems | |
1516 | ||
1517 | vector<Clustering> History::getAllQCDClusterings() { | |
1518 | vector<Clustering> ret; | |
1519 | // Initialise vectors to keep track of position of partons in the | |
1520 | // history-defining state | |
1521 | vector <int> PosFinalPartn; | |
1522 | vector <int> PosInitPartn; | |
1523 | vector <int> PosFinalGluon; | |
1524 | vector <int> PosFinalQuark; | |
1525 | vector <int> PosFinalAntiq; | |
1526 | vector <int> PosInitGluon; | |
1527 | vector <int> PosInitQuark; | |
1528 | vector <int> PosInitAntiq; | |
1529 | ||
1530 | // Search event record for final state particles and store these in | |
1531 | // quark, anti-quark and gluon vectors | |
1532 | for ( int i=0; i < state.size(); ++i ) | |
1533 | if ( state[i].isFinal() && state[i].colType() !=0 ) { | |
1534 | // Store final partons | |
1535 | if ( state[i].id() == 21 ) PosFinalGluon.push_back(i); | |
1536 | else if ( state[i].idAbs() < 10 && state[i].id() > 0) | |
1537 | PosFinalQuark.push_back(i); | |
1538 | else if ( state[i].idAbs() < 10 && state[i].id() < 0) | |
1539 | PosFinalAntiq.push_back(i); | |
1540 | } else if (state[i].status() == -21 && state[i].colType() != 0 ) { | |
1541 | // Store initial partons | |
1542 | if ( state[i].id() == 21 ) PosInitGluon.push_back(i); | |
1543 | else if ( state[i].idAbs() < 10 && state[i].id() > 0) | |
1544 | PosInitQuark.push_back(i); | |
1545 | else if ( state[i].idAbs() < 10 && state[i].id() < 0) | |
1546 | PosInitAntiq.push_back(i); | |
1547 | } | |
1548 | ||
1549 | // Get all clusterings for input state | |
1550 | vector<Clustering> systems; | |
1551 | systems = getQCDClusterings(state); | |
1552 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
1553 | systems.resize(0); | |
1554 | ||
1555 | // If valid clusterings were found, return | |
1556 | if ( !ret.empty() ) return ret; | |
1557 | // If no clusterings have been found until now, try to find | |
1558 | // clusterings of diagrams that interfere with the current one | |
1559 | // (i.e. change the colours of the current event slightly and run | |
1560 | // search again) | |
1561 | else if ( ret.empty() | |
1562 | && mergingHooksPtr->allowColourShuffling() ) { | |
1563 | Event NewState = Event(state); | |
1564 | // Start with changing final state quark colour | |
1565 | for(int i = 0; i < int(PosFinalQuark.size()); ++i) { | |
1566 | // Never change the hard process candidates | |
1567 | if ( mergingHooksPtr->hardProcess.matchesAnyOutgoing(PosFinalQuark[i], | |
1568 | NewState) ) | |
1569 | continue; | |
1570 | int col = NewState[PosFinalQuark[i]].col(); | |
1571 | for(int j = 0; j < int(PosInitAntiq.size()); ++j) { | |
1572 | // Now swap colours | |
1573 | int acl = NewState[PosInitAntiq[j]].acol(); | |
1574 | if ( col == acl ) continue; | |
1575 | NewState[PosFinalQuark[i]].col(acl); | |
1576 | NewState[PosInitAntiq[j]].acol(col); | |
1577 | systems = getQCDClusterings(NewState); | |
1578 | if (!systems.empty()) { | |
1579 | state = NewState; | |
1580 | NewState.clear(); | |
1581 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
1582 | systems.resize(0); | |
1583 | return ret; | |
1584 | } | |
1585 | } | |
1586 | } | |
1587 | // Now change final state antiquark anticolour | |
1588 | for(int i = 0; i < int(PosFinalAntiq.size()); ++i) { | |
1589 | // Never change the hard process candidates | |
1590 | if ( mergingHooksPtr->hardProcess.matchesAnyOutgoing(PosFinalAntiq[i], | |
1591 | NewState) ) | |
1592 | continue; | |
1593 | int acl = NewState[PosFinalAntiq[i]].acol(); | |
1594 | for(int j = 0; j < int(PosInitQuark.size()); ++j) { | |
1595 | // Now swap colours | |
1596 | int col = NewState[PosInitQuark[j]].col(); | |
1597 | if ( col == acl ) continue; | |
1598 | NewState[PosFinalAntiq[i]].acol(col); | |
1599 | NewState[PosInitQuark[j]].col(acl); | |
1600 | systems = getQCDClusterings(NewState); | |
1601 | if (!systems.empty()) { | |
1602 | state = NewState; | |
1603 | NewState.clear(); | |
1604 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
1605 | systems.resize(0); | |
1606 | return ret; | |
1607 | } | |
1608 | } | |
1609 | } | |
1610 | ||
1611 | if ( !ret.empty() ) { | |
1612 | string message="Warning in History::getAllQCDClusterings: Changed"; | |
1613 | message+=" colour structure to allow at least one clustering"; | |
1614 | infoPtr->errorMsg(message); | |
1615 | } | |
1616 | ||
1617 | // If no colour rearrangements should be done, print warning and return | |
1618 | } else { | |
1619 | string message="Warning in History::getAllQCDClusterings: No clusterings"; | |
1620 | message+=" found. History incomplete"; | |
1621 | infoPtr->errorMsg(message); | |
1622 | } | |
1623 | // Done | |
1624 | return ret; | |
1625 | } | |
1626 | ||
1627 | //-------------------------------------------------------------------------- | |
1628 | ||
1629 | // For one given state, find all possible clusterings. | |
1630 | // IN Event : state to be investigated | |
1631 | // OUT vector of all (rad,rec,emt) systems in the state | |
1632 | ||
1633 | vector<Clustering> History::getQCDClusterings( const Event& event) { | |
1634 | vector<Clustering> ret; | |
1635 | ||
1636 | // Initialise vectors to keep track of position of partons in the | |
1637 | // input event | |
1638 | vector <int> PosFinalPartn; | |
1639 | vector <int> PosInitPartn; | |
1640 | ||
1641 | vector <int> PosFinalGluon; | |
1642 | vector <int> PosFinalQuark; | |
1643 | vector <int> PosFinalAntiq; | |
1644 | vector <int> PosInitGluon; | |
1645 | vector <int> PosInitQuark; | |
1646 | vector <int> PosInitAntiq; | |
1647 | ||
1648 | // Search event record for final state particles and store these in | |
1649 | // quark, anti-quark and gluon vectors | |
1650 | for (int i=0; i < event.size(); ++i) | |
1651 | if ( event[i].isFinal() && event[i].colType() !=0 ) { | |
1652 | // Store final partons | |
1653 | PosFinalPartn.push_back(i); | |
1654 | if ( event[i].id() == 21 ) PosFinalGluon.push_back(i); | |
1655 | else if ( event[i].idAbs() < 10 && event[i].id() > 0) | |
1656 | PosFinalQuark.push_back(i); | |
1657 | else if ( event[i].idAbs() < 10 && event[i].id() < 0) | |
1658 | PosFinalAntiq.push_back(i); | |
1659 | } else if ( event[i].status() == -21 && event[i].colType() != 0 ) { | |
1660 | // Store initial partons | |
1661 | PosInitPartn.push_back(i); | |
1662 | if ( event[i].id() == 21 ) PosInitGluon.push_back(i); | |
1663 | else if ( event[i].idAbs() < 10 && event[i].id() > 0) | |
1664 | PosInitQuark.push_back(i); | |
1665 | else if ( event[i].idAbs() < 10 && event[i].id() < 0) | |
1666 | PosInitAntiq.push_back(i); | |
1667 | } | |
1668 | ||
1669 | int nFiGluon = int(PosFinalGluon.size()); | |
1670 | int nFiQuark = int(PosFinalQuark.size()); | |
1671 | int nFiAntiq = int(PosFinalAntiq.size()); | |
1672 | int nInGluon = int(PosInitGluon.size()); | |
1673 | int nInQuark = int(PosInitQuark.size()); | |
1674 | int nInAntiq = int(PosInitAntiq.size()); | |
1675 | ||
1676 | vector<Clustering> systems; | |
1677 | ||
1678 | // Find rad + emt + rec systems: | |
1679 | // (1) Start from gluon and find all (rad,rec,emt=gluon) triples | |
1680 | for (int i = 0; i < nFiGluon; ++i) { | |
1681 | //if ( mergingHooksPtr->nRecluster() > 0 | |
1682 | // && PosFinalGluon[i] == iReclusteredOld ) | |
1683 | // continue; | |
1684 | int EmtGluon = PosFinalGluon[i]; | |
1685 | systems = findQCDTriple( EmtGluon, 2, event, PosFinalPartn, PosInitPartn); | |
1686 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
1687 | systems.resize(0); | |
1688 | } | |
1689 | ||
1690 | // For more than one quark-antiquark pair in final state, check for | |
1691 | // g -> qqbar splittings | |
1692 | bool check_g2qq = true; | |
1693 | if ( ( ( nInQuark + nInAntiq == 0 ) | |
1694 | && (nInGluon == 0) | |
1695 | && (nFiQuark == 1) && (nFiAntiq == 1) ) | |
1696 | || ( ( nFiQuark + nFiAntiq == 0) | |
1697 | && (nInQuark == 1) && (nInAntiq == 1) ) ) | |
1698 | check_g2qq = false; | |
1699 | ||
1700 | if ( check_g2qq ) { | |
1701 | ||
1702 | // (2) Start from quark and find all (rad,rec,emt=quark) triples | |
1703 | // ( when g -> q qbar occured ) | |
1704 | for( int i=0; i < nFiQuark; ++i) { | |
1705 | //if ( mergingHooksPtr->nRecluster() > 0 | |
1706 | // && PosFinalQuark[i] == iReclusteredOld ) | |
1707 | // continue; | |
1708 | int EmtQuark = PosFinalQuark[i]; | |
1709 | systems = findQCDTriple( EmtQuark,1,event, PosFinalPartn, PosInitPartn); | |
1710 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
1711 | systems.resize(0); | |
1712 | } | |
1713 | ||
1714 | // (3) Start from anti-quark and find all (rad,rec,emt=anti-quark) | |
1715 | // triples ( when g -> q qbar occured ) | |
1716 | for( int i=0; i < nFiAntiq; ++i) { | |
1717 | //if ( mergingHooksPtr->nRecluster() > 0 | |
1718 | // && PosFinalAntiq[i] == iReclusteredOld ) | |
1719 | // continue; | |
1720 | int EmtAntiq = PosFinalAntiq[i]; | |
1721 | systems = findQCDTriple( EmtAntiq,1,event, PosFinalPartn, PosInitPartn); | |
1722 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
1723 | systems.resize(0); | |
1724 | } | |
1725 | } | |
1726 | ||
1727 | return ret; | |
1728 | } | |
1729 | ||
1730 | //-------------------------------------------------------------------------- | |
1731 | ||
1732 | // Function to construct (rad,rec,emt) triples from the event | |
1733 | // IN int : Position of Emitted in event record for which | |
1734 | // dipoles should be constructed | |
1735 | // int : Colour topogy to be tested | |
1736 | // 1= g -> qqbar, causing 2 -> 2 dipole splitting | |
1737 | // 2= q(bar) -> q(bar) g && g -> gg, | |
1738 | // causing a 2 -> 3 dipole splitting | |
1739 | // Event : event record to be checked for ptential partners | |
1740 | // OUT vector of all allowed radiator+recoiler+emitted triples | |
1741 | ||
1742 | vector<Clustering> History::findQCDTriple (int EmtTagIn, int colTopIn, | |
1743 | const Event& event, | |
1744 | vector<int> PosFinalPartn, | |
1745 | vector <int> PosInitPartn ) { | |
1746 | ||
1747 | // Copy input parton tag | |
1748 | int EmtTag = EmtTagIn; | |
1749 | // Copy input colour topology tag | |
1750 | // (1: g --> qqbar splitting present, 2:rest) | |
1751 | int colTop = colTopIn; | |
1752 | ||
1753 | // Initialise FinalSize | |
1754 | int FinalSize = int(PosFinalPartn.size()); | |
1755 | int InitSize = int(PosInitPartn.size()); | |
1756 | int Size = InitSize + FinalSize; | |
1757 | ||
1758 | vector<Clustering> clus; | |
1759 | ||
1760 | // Search final partons to find partons colour-connected to | |
1761 | // event[EmtTag], choose radiator, then choose recoiler | |
1762 | for ( int a = 0; a < Size; ++a ) { | |
1763 | int i = (a < FinalSize)? a : (a - FinalSize) ; | |
1764 | int iRad = (a < FinalSize)? PosFinalPartn[i] : PosInitPartn[i]; | |
1765 | ||
1766 | if ( event[iRad].col() == event[EmtTag].col() | |
1767 | && event[iRad].acol() == event[EmtTag].acol() ) | |
1768 | continue; | |
1769 | ||
1770 | //if ( mergingHooksPtr->nRecluster() > 0 | |
1771 | // && iRad == iReclusteredOld ) continue; | |
1772 | ||
1773 | if (iRad != EmtTag ) { | |
1774 | int pTdef = event[iRad].isFinal() ? 1 : -1; | |
1775 | int sign = (a < FinalSize)? 1 : -1 ; | |
1776 | ||
1777 | // First colour topology: g --> qqbar. Here, emt & rad should | |
1778 | // have same flavour (causes problems for gamma->qqbar). | |
1779 | if (colTop == 1) { | |
1780 | ||
1781 | if ( event[iRad].id() == -sign*event[EmtTag].id() ) { | |
1782 | int col = -1; | |
1783 | int acl = -1; | |
1784 | if (event[iRad].id() < 0) { | |
1785 | col = event[EmtTag].acol(); | |
1786 | acl = event[iRad].acol(); | |
1787 | } else { | |
1788 | col = event[EmtTag].col(); | |
1789 | acl = event[iRad].col(); | |
1790 | } | |
1791 | // Recoiler | |
1792 | int iRec = 0; | |
1793 | // Colour partner | |
1794 | int iPartner = 0; | |
1795 | ||
1796 | if (col > 0) { | |
1797 | // Find recoiler by colour | |
1798 | iRec = FindCol(col,iRad,EmtTag,event,1,true); | |
1799 | // In initial state splitting has final state colour partner, | |
1800 | // Save both partner and recoiler | |
1801 | if ( (sign < 0) && (event[iRec].isFinal()) ) { | |
1802 | // Save colour recoiler | |
1803 | iPartner = iRec; | |
1804 | // Reset kinematic recoiler to initial state parton | |
1805 | for(int l = 0; l < int(PosInitPartn.size()); ++l) | |
1806 | if (PosInitPartn[l] != iRad) iRec = PosInitPartn[l]; | |
1807 | // For final state splittings, colour partner and recoiler are | |
1808 | // identical | |
1809 | } else { | |
1810 | iPartner = iRec; | |
1811 | } | |
1812 | if ( iRec != 0 && iPartner != 0 | |
1813 | && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) { | |
1814 | clus.push_back( Clustering(EmtTag, iRad, iRec, iPartner, | |
1815 | pTLund(event[iRad], event[EmtTag], event[iRec], pTdef) )); | |
1816 | continue; | |
1817 | } | |
1818 | ||
1819 | // Reset partner | |
1820 | iPartner = 0; | |
1821 | // Find recoiler by colour | |
1822 | iRec = FindCol(col,iRad,EmtTag,event,2,true); | |
1823 | // In initial state splitting has final state colour partner, | |
1824 | // Save both partner and recoiler | |
1825 | if ( (sign < 0) && (event[iRec].isFinal()) ) { | |
1826 | // Save colour recoiler | |
1827 | iPartner = iRec; | |
1828 | // Reset kinematic recoiler to initial state parton | |
1829 | for(int l = 0; l < int(PosInitPartn.size()); ++l) | |
1830 | if (PosInitPartn[l] != iRad) iRec = PosInitPartn[l]; | |
1831 | // For final state splittings, colour partner and recoiler are | |
1832 | // identical | |
1833 | } else { | |
1834 | iPartner = iRec; | |
1835 | } | |
1836 | if ( iRec != 0 && iPartner != 0 | |
1837 | && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) { | |
1838 | clus.push_back( Clustering(EmtTag, iRad, iRec, iPartner, | |
1839 | pTLund(event[iRad], event[EmtTag], event[iRec], pTdef) )); | |
1840 | continue; | |
1841 | } | |
1842 | } | |
1843 | ||
1844 | if (acl > 0) { | |
1845 | ||
1846 | // Reset partner | |
1847 | iPartner = 0; | |
1848 | // Find recoiler by colour | |
1849 | iRec = FindCol(acl,iRad,EmtTag,event,1,true); | |
1850 | // In initial state splitting has final state colour partner, | |
1851 | // Save both partner and recoiler | |
1852 | if ( (sign < 0) && (event[iRec].isFinal()) ) { | |
1853 | // Save colour recoiler | |
1854 | iPartner = iRec; | |
1855 | // Reset kinematic recoiler to initial state parton | |
1856 | for(int l = 0; l < int(PosInitPartn.size()); ++l) | |
1857 | if (PosInitPartn[l] != iRad) iRec = PosInitPartn[l]; | |
1858 | // For final state splittings, colour partner and recoiler are | |
1859 | // identical | |
1860 | } else { | |
1861 | iPartner = iRec; | |
1862 | } | |
1863 | if ( iRec != 0 && iPartner != 0 | |
1864 | && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) { | |
1865 | clus.push_back( Clustering(EmtTag, iRad, iRec, iPartner, | |
1866 | pTLund(event[iRad], event[EmtTag], event[iRec], pTdef) )); | |
1867 | continue; | |
1868 | } | |
1869 | ||
1870 | // Reset partner | |
1871 | iPartner = 0; | |
1872 | // Find recoiler by colour | |
1873 | iRec = FindCol(acl,iRad,EmtTag,event,2,true); | |
1874 | // In initial state splitting has final state colour partner, | |
1875 | // Save both partner and recoiler | |
1876 | if ( (sign < 0) && (event[iRec].isFinal()) ) { | |
1877 | // Save colour recoiler | |
1878 | iPartner = iRec; | |
1879 | // Reset kinematic recoiler to initial state parton | |
1880 | for(int l = 0; l < int(PosInitPartn.size()); ++l) | |
1881 | if (PosInitPartn[l] != iRad) iRec = PosInitPartn[l]; | |
1882 | // For final state splittings, colour partner and recoiler are | |
1883 | // identical | |
1884 | } else { | |
1885 | iPartner = iRec; | |
1886 | } | |
1887 | if ( iRec != 0 && iPartner != 0 | |
1888 | && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) { | |
1889 | clus.push_back( Clustering(EmtTag, iRad, iRec, iPartner, | |
1890 | pTLund(event[iRad], event[EmtTag], event[iRec], pTdef) )); | |
1891 | continue; | |
1892 | } | |
1893 | } | |
1894 | // Initial gluon splitting | |
1895 | } else if ( event[iRad].id() == 21 | |
1896 | &&( event[iRad].col() == event[EmtTag].col() | |
1897 | || event[iRad].acol() == event[EmtTag].acol() )) { | |
1898 | // For an initial state radiator, always set recoiler | |
1899 | // to the other initial state parton (recoil is taken | |
1900 | // by full remaining system, so this is just a | |
1901 | // labelling for such a process) | |
1902 | int RecInit = 0; | |
1903 | for(int l = 0; l < int(PosInitPartn.size()); ++l) | |
1904 | if (PosInitPartn[l] != iRad) RecInit = PosInitPartn[l]; | |
1905 | ||
1906 | // Find the colour connected partner | |
1907 | // Find colour index of radiator before splitting | |
1908 | int col = getRadBeforeCol(iRad, EmtTag, event); | |
1909 | int acl = getRadBeforeAcol(iRad, EmtTag, event); | |
1910 | ||
1911 | // Find the correct partner: If a colour line has split, | |
1912 | // the partner is connected to the radiator before the splitting | |
1913 | // by a colour line (same reasoning for anticolour). The colour | |
1914 | // that split is the colour appearing twice in the | |
1915 | // radiator + emitted pair. | |
1916 | // Thus, if we remove a colour index with the clustering, | |
1917 | // we should look for a colour partner, else look for | |
1918 | // an anticolour partner | |
1919 | int colRemove = (event[iRad].col() == event[EmtTag].col()) | |
1920 | ? event[iRad].col() : 0; | |
1921 | ||
1922 | int iPartner = 0; | |
1923 | if (colRemove > 0 && col > 0) | |
1924 | iPartner = FindCol(col,iRad,EmtTag,event,1,true) | |
1925 | + FindCol(col,iRad,EmtTag,event,2,true); | |
1926 | else if (colRemove > 0 && acl > 0) | |
1927 | iPartner = FindCol(acl,iRad,EmtTag,event,1,true) | |
1928 | + FindCol(acl,iRad,EmtTag,event,2,true); | |
1929 | ||
1930 | if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event ) ) { | |
1931 | clus.push_back( Clustering(EmtTag, iRad, RecInit, iPartner, | |
1932 | pTLund(event[iRad],event[EmtTag],event[RecInit], pTdef) )); | |
1933 | continue; | |
1934 | } | |
1935 | } | |
1936 | ||
1937 | // Second colour topology: Gluon emission | |
1938 | ||
1939 | } else { | |
1940 | if ( (event[iRad].col() == event[EmtTag].acol()) | |
1941 | || (event[iRad].acol() == event[EmtTag].col()) | |
1942 | || (event[iRad].col() == event[EmtTag].col()) | |
1943 | || (event[iRad].acol() == event[EmtTag].acol()) ) { | |
1944 | // For the rest, choose recoiler to have a common colour | |
1945 | // tag with radiator, while not being the "Emitted" | |
1946 | ||
1947 | int col = -1; | |
1948 | int acl = -1; | |
1949 | ||
1950 | if (event[iRad].isFinal() ) { | |
1951 | ||
1952 | if ( event[iRad].id() < 0) { | |
1953 | acl = event[EmtTag].acol(); | |
1954 | col = event[iRad].col(); | |
1955 | } else if ( event[iRad].id() > 0 && event[iRad].id() < 10) { | |
1956 | col = event[EmtTag].col(); | |
1957 | acl = event[iRad].acol(); | |
1958 | } else { | |
1959 | col = event[iRad].col(); | |
1960 | acl = event[iRad].acol(); | |
1961 | } | |
1962 | ||
1963 | int iRec = 0; | |
1964 | if (col > 0) { | |
1965 | iRec = FindCol(col,iRad,EmtTag,event,1,true); | |
1966 | if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0; | |
1967 | if (iRec != 0 | |
1968 | && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) { | |
1969 | clus.push_back( Clustering(EmtTag, iRad, iRec, iRec, | |
1970 | pTLund(event[iRad],event[EmtTag],event[iRec], pTdef) )); | |
1971 | continue; | |
1972 | } | |
1973 | ||
1974 | iRec = FindCol(col,iRad,EmtTag,event,2,true); | |
1975 | if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0; | |
1976 | if (iRec != 0 | |
1977 | && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) { | |
1978 | clus.push_back( Clustering(EmtTag, iRad, iRec, iRec, | |
1979 | pTLund(event[iRad],event[EmtTag],event[iRec], pTdef) )); | |
1980 | continue; | |
1981 | } | |
1982 | } | |
1983 | ||
1984 | if (acl > 0) { | |
1985 | iRec = FindCol(acl,iRad,EmtTag,event,1,true); | |
1986 | if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0; | |
1987 | if (iRec != 0 | |
1988 | && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) { | |
1989 | clus.push_back( Clustering(EmtTag, iRad, iRec, iRec, | |
1990 | pTLund(event[iRad],event[EmtTag],event[iRec], pTdef) )); | |
1991 | continue; | |
1992 | } | |
1993 | ||
1994 | iRec = FindCol(acl,iRad,EmtTag,event,2,true); | |
1995 | if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0; | |
1996 | if (iRec != 0 | |
1997 | && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) { | |
1998 | clus.push_back( Clustering(EmtTag, iRad, iRec, iRec, | |
1999 | pTLund(event[iRad],event[EmtTag],event[iRec], pTdef) )); | |
2000 | continue; | |
2001 | } | |
2002 | } | |
2003 | ||
2004 | } else { | |
2005 | ||
2006 | // For an initial state radiator, always set recoiler | |
2007 | // to the other initial state parton (recoil is taken | |
2008 | ||
2009 | // by full remaining system, so this is just a | |
2010 | // labelling for such a process) | |
2011 | int RecInit = 0; | |
2012 | int iPartner = 0; | |
2013 | for(int l = 0; l < int(PosInitPartn.size()); ++l) | |
2014 | if (PosInitPartn[l] != iRad) RecInit = PosInitPartn[l]; | |
2015 | ||
2016 | // Find the colour connected partner | |
2017 | // Find colour index of radiator before splitting | |
2018 | col = getRadBeforeCol(iRad, EmtTag, event); | |
2019 | acl = getRadBeforeAcol(iRad, EmtTag, event); | |
2020 | ||
2021 | // Find the correct partner: If a colour line has split, | |
2022 | // the partner is connected to the radiator before the splitting | |
2023 | // by a colour line (same reasoning for anticolour). The colour | |
2024 | // that split is the colour appearing twice in the | |
2025 | // radiator + emitted pair. | |
2026 | // Thus, if we remove a colour index with the clustering, | |
2027 | // we should look for a colour partner, else look for | |
2028 | // an anticolour partner | |
2029 | int colRemove = (event[iRad].col() == event[EmtTag].col()) | |
2030 | ? event[iRad].col() : 0; | |
2031 | iPartner = (colRemove > 0) | |
2032 | ? FindCol(col,iRad,EmtTag,event,1,true) | |
2033 | + FindCol(col,iRad,EmtTag,event,2,true) | |
2034 | : FindCol(acl,iRad,EmtTag,event,1,true) | |
2035 | + FindCol(acl,iRad,EmtTag,event,2,true); | |
2036 | ||
2037 | if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event)) { | |
2038 | clus.push_back( Clustering(EmtTag, iRad, RecInit, iPartner, | |
2039 | pTLund(event[iRad],event[EmtTag],event[RecInit], pTdef))); | |
2040 | ||
2041 | continue; | |
2042 | } | |
2043 | } | |
2044 | } | |
2045 | } | |
2046 | } | |
2047 | } | |
2048 | ||
2049 | // Done | |
2050 | return clus; | |
2051 | } | |
2052 | ||
2053 | //-------------------------------------------------------------------------- | |
2054 | ||
2055 | // For the history-defining state (and if necessary interfering | |
2056 | // states), find all possible clusterings. | |
2057 | // NO INPUT | |
2058 | // OUT vector of all (rad,rec,emt) systems | |
2059 | ||
2060 | vector<Clustering> History::getAllEWClusterings() { | |
2061 | vector<Clustering> ret; | |
2062 | ||
2063 | // Get all clusterings for input state | |
2064 | vector<Clustering> systems; | |
2065 | systems = getEWClusterings(state); | |
2066 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
2067 | // Done | |
2068 | return ret; | |
2069 | } | |
2070 | ||
2071 | //-------------------------------------------------------------------------- | |
2072 | ||
2073 | // For one given state, find all possible clusterings. | |
2074 | // IN Event : state to be investigated | |
2075 | // OUT vector of all (rad,rec,emt) systems in the state | |
2076 | ||
2077 | vector<Clustering> History::getEWClusterings( const Event& event) { | |
2078 | vector<Clustering> ret; | |
2079 | ||
2080 | // Initialise vectors to keep track of position of partons in the | |
2081 | // input event | |
2082 | vector <int> PosFinalPartn; | |
2083 | vector <int> PosInitPartn; | |
2084 | vector <int> PosFinalW; | |
2085 | ||
2086 | // Search event record for final state particles and store these in | |
2087 | // quark, anti-quark and gluon vectors | |
2088 | for ( int i=0; i < event.size(); ++i ) | |
2089 | if ( event[i].isFinal() && abs(event[i].colType()) == 1 ) { | |
2090 | // Store final partons | |
2091 | PosFinalPartn.push_back(i); | |
2092 | } else if ( event[i].status() == -21 && abs(event[i].colType()) == 1 ) { | |
2093 | // Store initial partons | |
2094 | PosInitPartn.push_back(i); | |
2095 | } | |
2096 | // Search event record for final W | |
2097 | for ( int i=0; i < event.size(); ++i ) | |
2098 | if ( event[i].isFinal() && event[i].idAbs() == 24 ) | |
2099 | PosFinalW.push_back( i ); | |
2100 | ||
2101 | vector<Clustering> systems; | |
2102 | // Find rad + emt + rec systems: | |
2103 | // (1) Start from W boson and find all (rad,rec,emt=W) triples | |
2104 | for ( int i = 0; i < int(PosFinalW.size()); ++i ) { | |
2105 | int EmtW = PosFinalW[i]; | |
2106 | systems = findEWTriple( EmtW, event, PosFinalPartn); | |
2107 | ret.insert(ret.end(), systems.begin(), systems.end()); | |
2108 | systems.resize(0); | |
2109 | } | |
2110 | ||
2111 | return ret; | |
2112 | } | |
2113 | ||
2114 | //-------------------------------------------------------------------------- | |
2115 | ||
2116 | // Function to construct (rad,rec,emt) triples from the event | |
2117 | // IN int : Position of Emitted in event record for which | |
2118 | // dipoles should be constructed | |
2119 | // int : Colour topogy to be tested | |
2120 | // 1= g -> qqbar, causing 2 -> 2 dipole splitting | |
2121 | // 2= q(bar) -> q(bar) g && g -> gg, | |
2122 | // causing a 2 -> 3 dipole splitting | |
2123 | // Event : event record to be checked for ptential partners | |
2124 | // OUT vector of all allowed radiator+recoiler+emitted triples | |
2125 | ||
2126 | vector<Clustering> History::findEWTriple ( int EmtTagIn, const Event& event, | |
2127 | vector<int> PosFinalPartn ) { | |
2128 | // Copy input parton tag | |
2129 | int EmtTag = EmtTagIn; | |
2130 | // Copy input colour topology tag | |
2131 | // (1: g --> qqbar splitting present, 2:rest) | |
2132 | ||
2133 | // Initialise FinalSize | |
2134 | int FinalSize = int(PosFinalPartn.size()); | |
2135 | ||
2136 | vector<Clustering> clus; | |
2137 | ||
2138 | // Search final partons to find partons colour-connected to | |
2139 | // event[EmtTag], choose radiator, then choose recoiler | |
2140 | for ( int a = 0; a < FinalSize; ++a ) { | |
2141 | int iRad = PosFinalPartn[a]; | |
2142 | if (iRad != EmtTag ) { | |
2143 | int pTdef = 1; | |
2144 | // Find recoiler by flavour. | |
2145 | int flavRad = event[iRad].id(); | |
2146 | int flavEmt = event[EmtTag].id(); | |
2147 | ||
2148 | // Loop through final partons and try to find matching flavours. | |
2149 | for ( int i = 0; i < FinalSize; ++i ) { | |
2150 | int iRec = PosFinalPartn[i]; | |
2151 | if ( i != a && flavEmt > 0 | |
2152 | && event[iRec].id() == -flavRad - 1 ) | |
2153 | clus.push_back( Clustering(EmtTag, iRad, iRec, iRec, | |
2154 | pTLund(event[iRad],event[EmtTag],event[iRec], pTdef) ) ); | |
2155 | } | |
2156 | } | |
2157 | } | |
2158 | ||
2159 | // Done | |
2160 | return clus; | |
2161 | } | |
2162 | ||
2163 | //-------------------------------------------------------------------------- | |
2164 | ||
2165 | // Calculate and return the probability of a clustering. | |
2166 | // IN Clustering : rad,rec,emt - System for which the splitting | |
2167 | // probability should be calcuated | |
2168 | // OUT splitting probability | |
2169 | ||
2170 | double History::getProb(const Clustering & SystemIn) { | |
2171 | ||
2172 | // Get local copies of input system | |
2173 | int Rad = SystemIn.emittor; | |
2174 | int Rec = SystemIn.recoiler; | |
2175 | int Emt = SystemIn.emitted; | |
2176 | ||
2177 | // Initialise shower probability | |
2178 | double showerProb = 0.0; | |
2179 | ||
2180 | // If the splitting resulted in disallowed evolution variable, | |
2181 | // disallow the splitting | |
2182 | if (SystemIn.pT() <= 0.) return 0.; | |
2183 | ||
2184 | // Initialise all combinatorical factors | |
2185 | double CF = 4./3.; | |
2186 | double NC = 3.; | |
2187 | // Flavour is known when reclustring, thus n_f=1 | |
2188 | double TR = 1./2.; | |
2189 | ||
2190 | // Split up in FSR and ISR | |
2191 | bool isFSR = (state[Rad].isFinal() && state[Rec].isFinal()); | |
2192 | bool isFSRinREC = (state[Rad].isFinal() && !state[Rec].isFinal()); | |
2193 | bool isISR = !state[Rad].isFinal(); | |
2194 | ||
2195 | // Check if this is the clustering 2->3 to 2->2. | |
2196 | // If so, use weight for joined evolution | |
2197 | int nFinal = 0; | |
2198 | for(int i=0; i < state.size(); ++i) | |
2199 | if (state[i].isFinal()) nFinal++; | |
2200 | bool isLast = (nFinal == (mergingHooksPtr->hardProcess.nQuarksOut() | |
2201 | +mergingHooksPtr->hardProcess.nLeptonOut()+1)); | |
2202 | ||
2203 | if (isISR) { | |
2204 | // Find incoming particles | |
2205 | ||
2206 | int inP = 0; | |
2207 | int inM = 0; | |
2208 | for(int i=0;i< int(state.size()); ++i) { | |
2209 | if (state[i].mother1() == 1) inP = i; | |
2210 | if (state[i].mother1() == 2) inM = i; | |
2211 | } | |
2212 | // Construct dipole mass, eCM and sHat = x1*x2*s | |
2213 | Vec4 sum = state[Rad].p() + state[Rec].p() - state[Emt].p(); | |
2214 | double m2Dip = sum.m2Calc(); | |
2215 | double sHat = (state[inM].p() + state[inP].p()).m2Calc(); | |
2216 | // Energy fraction z=E_q1/E_qi in branch q(i)q(2) -> q(1)g(3)q(2) | |
2217 | double z1 = m2Dip / sHat; | |
2218 | // Virtuality of the splittings | |
2219 | Vec4 Q1( state[Rad].p() - state[Emt].p() ); | |
2220 | Vec4 Q2( state[Rec].p() - state[Emt].p() ); | |
2221 | // Q^2 for emission off radiator line | |
2222 | double Q1sq = -Q1.m2Calc(); | |
2223 | // pT^2 for emission off radiator line | |
2224 | double pT1sq = pow(SystemIn.pT(),2); | |
2225 | // Remember if massive particles involved: Mass corrections for | |
2226 | // to g->QQ and Q->Qg splittings | |
2227 | bool g2QQmassive = mergingHooksPtr->includeMassive() | |
2228 | && state[Rad].id() == 21 | |
2229 | && ( state[Emt].idAbs() >= 4 && state[Emt].idAbs() < 7); | |
2230 | bool Q2Qgmassive = mergingHooksPtr->includeMassive() | |
2231 | && state[Emt].id() == 21 | |
2232 | && ( state[Rad].idAbs() >= 4 && state[Rad].idAbs() < 7); | |
2233 | bool isMassive = mergingHooksPtr->includeMassive() | |
2234 | && (g2QQmassive || Q2Qgmassive); | |
2235 | double m2Emt0 = pow(particleDataPtr->m0(state[Emt].id()),2); | |
2236 | double m2Rad0 = pow(particleDataPtr->m0(state[Rad].id()),2); | |
2237 | ||
2238 | // Correction of virtuality for massive splittings | |
2239 | if ( g2QQmassive) Q1sq += m2Emt0; | |
2240 | else if (Q2Qgmassive) Q1sq += m2Rad0; | |
2241 | ||
2242 | // pT0 dependence!!! | |
2243 | double pT0sq = pow(mergingHooksPtr->pT0ISR(),2); | |
2244 | double Q2sq = -Q2.m2Calc(); | |
2245 | ||
2246 | // Correction of virtuality of other splitting | |
2247 | bool g2QQmassiveRec = mergingHooksPtr->includeMassive() | |
2248 | && state[Rec].id() == 21 | |
2249 | && ( state[Emt].idAbs() >= 4 && state[Emt].idAbs() < 7); | |
2250 | bool Q2QgmassiveRec = mergingHooksPtr->includeMassive() | |
2251 | && state[Emt].id() == 21 | |
2252 | && ( state[Rec].idAbs() >= 4 && state[Rec].idAbs() < 7); | |
2253 | double m2Rec0 = pow(particleDataPtr->m0(state[Rad].id()),2); | |
2254 | if ( g2QQmassiveRec) Q2sq += m2Emt0; | |
2255 | else if (Q2QgmassiveRec) Q2sq += m2Rec0; | |
2256 | ||
2257 | bool hasJoinedEvol = (state[Emt].id() == 21 | |
2258 | || state[Rad].id() == state[Rec].id()); | |
2259 | ||
2260 | // Initialise normalization factor multiplying the splitting | |
2261 | // function numerator | |
2262 | double fac = 1.; | |
2263 | if ( mergingHooksPtr->pickByFull() || mergingHooksPtr->pickBySumPT()) { | |
2264 | double facJoined = ( Q2sq + pT0sq/(1.-z1) ) | |
2265 | * 1./(Q1sq*Q2sq + pT0sq*sHat + pow(pT0sq/(1.-z1),2)); | |
2266 | double facSingle = mergingHooksPtr->nonJoinedNorm()*1./( pT1sq + pT0sq); | |
2267 | ||
2268 | fac = (hasJoinedEvol && isLast) ? facJoined : facSingle; | |
2269 | ||
2270 | } else if (mergingHooksPtr->pickByPoPT2()) { | |
2271 | fac = 1./(pT1sq + pT0sq); | |
2272 | } | |
2273 | ||
2274 | // Calculate shower splitting probability: | |
2275 | // Splitting functions*normalization*ME reweighting factors | |
2276 | ||
2277 | // Calculate branching probability for q -> q g | |
2278 | if ( state[Emt].id() == 21 && state[Rad].id() != 21) { | |
2279 | // Find splitting kernel | |
2280 | double num = CF*(1. + pow(z1,2)) / (1.-z1); | |
2281 | if (isMassive) num -= CF * z1 * (1.-z1) * (m2Rad0/pT1sq); | |
2282 | ||
2283 | // Find ME reweighting factor | |
2284 | double meReweighting = 1.; | |
2285 | // Find the number of final state coloured particles, apart | |
2286 | // from those coming from the hard process | |
2287 | int nCol = 0; | |
2288 | for(int i=0; i < state.size(); ++i) | |
2289 | if (state[i].isFinal() && state[i].colType() != 0 | |
2290 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) ) | |
2291 | nCol++; | |
2292 | // For first splitting of single vector boson production, | |
2293 | // apply ME corrections | |
2294 | if (nCol == 1 | |
2295 | && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1) { | |
2296 | double sH = m2Dip / z1; | |
2297 | double tH = -Q1sq; | |
2298 | double uH = Q1sq - m2Dip * (1. - z1) / z1; | |
2299 | double misMatch = (uH*tH - (uH + tH)*pT0sq/(1.-z1) | |
2300 | + pow(pT0sq/(1.-z1),2) ) / (uH*tH); | |
2301 | meReweighting *= (tH*tH + uH*uH + 2. * m2Dip * sH) | |
2302 | / (sH*sH + m2Dip*m2Dip); | |
2303 | meReweighting *= misMatch; | |
2304 | } | |
2305 | // Multiply factors | |
2306 | showerProb = num*fac*meReweighting; | |
2307 | ||
2308 | // Calculate branching probability for g -> g g | |
2309 | } else if ( state[Emt].id() == 21 && state[Rad].id() == 21) { | |
2310 | // Calculate splitting kernel | |
2311 | double num = 2.*NC*pow2(1. - z1*(1.-z1)) / (z1*(1.-z1)); | |
2312 | ||
2313 | // Include ME reweighting for higgs!! | |
2314 | // Find ME reweighting factor | |
2315 | double meReweighting = 1.; | |
2316 | // Find the number of final state coloured particles, apart | |
2317 | // from those coming from the hard process | |
2318 | int nCol = 0; | |
2319 | for(int i=0; i < state.size(); ++i) | |
2320 | if (state[i].isFinal() && state[i].colType() != 0 | |
2321 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) ) | |
2322 | nCol++; | |
2323 | // For first splitting of single vector boson production, | |
2324 | // apply ME corrections | |
2325 | if ( nCol == 1 | |
2326 | && mergingHooksPtr->getProcessString().compare("pp>h") == 0 | |
2327 | && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1 ) { | |
2328 | double sH = m2Dip / z1; | |
2329 | double tH = -Q1sq; | |
2330 | double uH = Q1sq - m2Dip * (1. - z1) / z1; | |
2331 | meReweighting *= 0.5 | |
2332 | * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(m2Dip)) | |
2333 | / pow2(sH*sH - m2Dip * (sH - m2Dip)); | |
2334 | } | |
2335 | ||
2336 | // Multiply factors | |
2337 | showerProb = num*fac*meReweighting; | |
2338 | ||
2339 | // Calculate branching probability for q -> g q | |
2340 | } else if ( state[Emt].id() != 21 && state[Rad].id() != 21 ) { | |
2341 | // Calculate splitting kernel | |
2342 | double num = CF*(1. + pow2(1.-z1)) / z1; | |
2343 | ||
2344 | // Include ME reweighting for higgs!! | |
2345 | // Find ME reweighting factor | |
2346 | double meReweighting = 1.; | |
2347 | // Find the number of final state coloured particles, apart | |
2348 | // from those coming from the hard process | |
2349 | int nCol = 0; | |
2350 | for ( int i=0; i < state.size(); ++i ) | |
2351 | if ( state[i].isFinal() && state[i].colType() != 0 | |
2352 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) ) | |
2353 | nCol++; | |
2354 | // For first splitting of single vector boson production, | |
2355 | // apply ME corrections | |
2356 | if (nCol == 1 | |
2357 | && mergingHooksPtr->getProcessString().compare("pp>h") == 0 | |
2358 | && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1) { | |
2359 | double sH = m2Dip / z1; | |
2360 | double uH = Q1sq - m2Dip * (1. - z1) / z1; | |
2361 | meReweighting *= (sH*sH + uH*uH) | |
2362 | / (sH*sH + pow2(sH -m2Dip)); | |
2363 | } | |
2364 | ||
2365 | // Multiply factors | |
2366 | showerProb = num*fac*meReweighting; | |
2367 | ||
2368 | // Calculate branching probability for g -> q qbar | |
2369 | } else if ( state[Emt].id() != 21 && state[Rad].id() == 21 ) { | |
2370 | // Calculate splitting kernel | |
2371 | double num = TR * ( pow(z1,2) + pow(1.-z1,2) ); | |
2372 | if (isMassive) num += TR * 2.*z1*(1.-z1)*(m2Emt0/pT1sq); | |
2373 | // Calculate ME reweighting factor | |
2374 | double meReweighting = 1.; | |
2375 | // Find the number of final state coloured particles, apart | |
2376 | // from those coming from the hard process | |
2377 | int nCol = 0; | |
2378 | for ( int i=0; i < state.size(); ++i ) | |
2379 | if ( state[i].isFinal() && state[i].colType() != 0 | |
2380 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) ) | |
2381 | nCol++; | |
2382 | // For first splitting of single vector boson production, | |
2383 | // apply ME corrections | |
2384 | if (nCol == 1 | |
2385 | && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1) { | |
2386 | double sH = m2Dip / z1; | |
2387 | double tH = -Q1sq; | |
2388 | double uH = Q1sq - m2Dip * (1. - z1) / z1; | |
2389 | swap( tH, uH); | |
2390 | double misMatch = ( uH - pT0sq/(1.-z1) ) / uH; | |
2391 | double me = (sH*sH + uH*uH + 2. * m2Dip * tH) | |
2392 | / (pow2(sH - m2Dip) + m2Dip*m2Dip); | |
2393 | // Weight with me/overestimate | |
2394 | meReweighting *= me; | |
2395 | meReweighting *= misMatch; | |
2396 | } | |
2397 | // Multiply factors | |
2398 | showerProb = num*fac*meReweighting; | |
2399 | ||
2400 | // Print error if no kernel calculated | |
2401 | } else { | |
2402 | string message="Error in History::getProb: Splitting kernel undefined"; | |
2403 | message+=" in ISR clustering."; | |
2404 | infoPtr->errorMsg(message); | |
2405 | } | |
2406 | ||
2407 | // If corrected pT below zero in ISR, put probability to zero | |
2408 | double m2Sister0 = pow(state[Emt].m0(),2); | |
2409 | double pT2corr = (Q1sq - z1*(m2Dip + Q1sq)*(Q1sq + m2Sister0)/m2Dip); | |
2410 | if (pT2corr < 0.) showerProb *= 1e-9; | |
2411 | ||
2412 | // If creating heavy quark by Q -> gQ then next need g -> Q + Qbar. | |
2413 | // So minimum total mass2 is 4 * m2Sister, but use more to be safe. | |
2414 | if ( state[Emt].id() == state[Rad].id() | |
2415 | && ( state[Rad].idAbs() == 4 || state[Rad].idAbs() == 5 )) { | |
2416 | double m2QQsister = 2.*4.*m2Sister0; | |
2417 | double pT2QQcorr = Q1sq - z1*(m2Dip + Q1sq)*(Q1sq + m2QQsister) | |
2418 | / m2Dip; | |
2419 | if (pT2QQcorr < 0.0) showerProb *= 1e-9; | |
2420 | } | |
2421 | ||
2422 | if (mergingHooksPtr->includeRedundant()) { | |
2423 | // Initialise the spacelike shower alpha_S | |
2424 | AlphaStrong* asISR = mergingHooksPtr->AlphaS_ISR(); | |
2425 | double as = (*asISR).alphaS(pT1sq + pT0sq) / (2.*M_PI); | |
2426 | // Multiply with alpha_S | |
2427 | showerProb *= as; | |
2428 | } | |
2429 | ||
2430 | // Done for ISR case, begin FSR case | |
2431 | ||
2432 | } else if (isFSR || isFSRinREC) { | |
2433 | ||
2434 | // Construct dipole mass | |
2435 | Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p(); | |
2436 | double m2Dip = sum.m2Calc(); | |
2437 | // Construct 2->3 variables for FSR | |
2438 | double x1 = 2. * (sum * state[Rad].p()) / m2Dip; | |
2439 | double x2 = 2. * (sum * state[Rec].p()) / m2Dip; | |
2440 | double prop1 = max(1e-12, 1. - x1); | |
2441 | double prop2 = max(1e-12, 1. - x2); | |
2442 | double x3 = max(1e-12, 2. - x1 - x2); | |
2443 | // Energy fraction z=E_q1/E_qi in branch q(i)q(2) -> q(1)g(3)q(2) | |
2444 | double z1 = x1/(x1 + x3); | |
2445 | ||
2446 | // Virtuality of the splittings | |
2447 | Vec4 Q1( state[Rad].p() + state[Emt].p() ); | |
2448 | Vec4 Q2( state[Rec].p() + state[Emt].p() ); | |
2449 | ||
2450 | // Q^2 for emission off radiator line | |
2451 | double Q1sq = Q1.m2Calc(); | |
2452 | // pT^2 for emission off radiator line | |
2453 | double pT1sq = pow(SystemIn.pT(),2); | |
2454 | // Q^2 for emission off recoiler line | |
2455 | double Q2sq = Q2.m2Calc(); | |
2456 | ||
2457 | // Correction of virtuality for massive splittings | |
2458 | double m2Rad0 = pow(particleDataPtr->m0(state[Rad].id()),2); | |
2459 | double m2Rec0 = pow(particleDataPtr->m0(state[Rad].id()),2); | |
2460 | if ( state[Rad].idAbs() >= 4 && state[Rad].idAbs() < 7) | |
2461 | Q1sq -= m2Rad0; | |
2462 | if ( state[Rec].idAbs() >= 4 && state[Rec].idAbs() < 7) | |
2463 | Q2sq -= m2Rec0; | |
2464 | ||
2465 | // Initialise normalization factor multiplying the splitting | |
2466 | // function numerator | |
2467 | double fac = 1.; | |
2468 | if ( mergingHooksPtr->pickByFull() || mergingHooksPtr->pickBySumPT()) { | |
2469 | double facJoined = (1.-z1)/Q1sq * m2Dip/( Q1sq + Q2sq ); | |
2470 | double facSingle = mergingHooksPtr->fsrInRecNorm() * 1./ pT1sq; | |
2471 | fac = (!isFSRinREC && isLast) ? facJoined : facSingle; | |
2472 | } else if (mergingHooksPtr->pickByPoPT2()) { | |
2473 | fac = 1. / pT1sq; | |
2474 | } | |
2475 | // Calculate shower splitting probability: | |
2476 | // Splitting functions*normalization*ME reweighting factors | |
2477 | ||
2478 | // Calculate branching probability for g -> g_1 g_2 | |
2479 | if ( state[Emt].id() == 21 && state[Rad].id() == 21) { | |
2480 | // Calculate splitting kernel | |
2481 | double num = 0.5* NC * (1. + pow3(z1)) / (1.-z1); | |
2482 | // Multiply factors | |
2483 | showerProb = num*fac; | |
2484 | ||
2485 | // Calculate branching probability for q -> q g with quark recoiler | |
2486 | } else if ( state[Emt].id() == 21 | |
2487 | && state[Rad].id() != 21 && state[Rec].id() != 21) { | |
2488 | // For a qqbar dipole in FSR, ME corrections exist and the | |
2489 | // splitting function "z-weight" is set to 1.0 (only for 2->2 ??) | |
2490 | double num = CF * 2./(1.-z1); | |
2491 | // Find the number of final state coloured particles, apart | |
2492 | // from those coming from the hard process | |
2493 | int nCol = 0; | |
2494 | for(int i=0; i < state.size(); ++i) | |
2495 | if (state[i].isFinal() && state[i].colType() != 0 | |
2496 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) ) | |
2497 | nCol++; | |
2498 | // Calculate splitting kernel | |
2499 | if ( nCol > 3 | |
2500 | || int(mergingHooksPtr->hardProcess.hardIntermediate.size()) > 1) | |
2501 | num = CF * (1. + pow2(z1)) /(1.-z1); | |
2502 | // Calculate ME reweighting factor | |
2503 | double meReweighting = 1.; | |
2504 | // Correct if this is the process created by the first | |
2505 | // FSR splitting of a 2->2 process | |
2506 | if ( nCol == 3 | |
2507 | && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1 ) { | |
2508 | // Calculate the ME reweighting factor | |
2509 | double ShowerRate1 = 2./( x3 * prop2 ); | |
2510 | double meDividingFactor1 = prop1 / x3; | |
2511 | double me = (pow(x1,2) + pow(x2,2))/(prop1*prop2); | |
2512 | meReweighting = meDividingFactor1 * me / ShowerRate1; | |
2513 | } | |
2514 | // Multiply factors | |
2515 | showerProb = num*fac*meReweighting; | |
2516 | ||
2517 | // Calculate branching probability for q1 -> q2 w+- with quark recoiler | |
2518 | } else if ( state[Emt].idAbs() == 24 | |
2519 | && state[Rad].id() != 21 && state[Rec].id() != 21 ) { | |
2520 | double m2W = state[Emt].p().m2Calc(); | |
2521 | double num = ( 3.*pow2(m2W / m2Dip) | |
2522 | + 2.* (m2W/m2Dip)*(x1 + x2) | |
2523 | + pow2(x1) + pow2(x2) ) / ( prop1*prop2 ) | |
2524 | - (m2W/m2Dip) / pow2(prop1) | |
2525 | - (m2W/m2Dip) / pow2(prop2); | |
2526 | // Multiply factors | |
2527 | showerProb = num*fac; | |
2528 | ||
2529 | // Calculate branching probability for q -> q g with gluon recoiler | |
2530 | } else if ( state[Emt].id() == 21 && state[Rad].id() != 21 | |
2531 | && state[Rec].id() == 21 ) { | |
2532 | // For qg /qbarg dipoles, the splitting function is | |
2533 | // calculated and not weighted by a ME correction factor | |
2534 | // Shower splitting function | |
2535 | double num = CF * (1. + pow2(z1)) / (1.-z1); | |
2536 | showerProb = num*fac; | |
2537 | ||
2538 | // Calculate branching probability for g -> q qbar | |
2539 | } else if ( state[Emt].id() != 21 ) { | |
2540 | // Get flavour of quark / antiquark | |
2541 | int flavour = state[Emt].id(); | |
2542 | // Get correct masses for the quarks | |
2543 | // (needed to calculate splitting function?) | |
2544 | double mFlavour = particleDataPtr->m0(flavour); | |
2545 | // Get mass of quark/antiquark pair | |
2546 | double mDipole = m(state[Rad].p(), state[Emt].p()); | |
2547 | // Factor determining if gluon decay was allowed | |
2548 | double beta = sqrtpos( 1. - 4.*pow2(mFlavour)/pow2(mDipole) ); | |
2549 | // Shower splitting function | |
2550 | double num = 0.5*TR * ( z1*z1 + (1.-z1)*(1.-z1) ); | |
2551 | ||
2552 | showerProb = num*fac*beta; | |
2553 | ||
2554 | // Print error if no kernel calculated | |
2555 | } else { | |
2556 | string message="Error in History::getProb: Splitting kernel undefined"; | |
2557 | message+=" in FSR clustering."; | |
2558 | infoPtr->errorMsg(message); | |
2559 | } | |
2560 | ||
2561 | if (mergingHooksPtr->includeRedundant()) { | |
2562 | // Initialise the spacelike shower alpha_S | |
2563 | AlphaStrong* asFSR = mergingHooksPtr->AlphaS_FSR(); | |
2564 | double as = (*asFSR).alphaS(pT1sq) / (2.*M_PI); | |
2565 | // Multiply with alpha_S | |
2566 | showerProb *= as; | |
2567 | } | |
2568 | ||
2569 | // Done for FSR | |
2570 | } else { | |
2571 | string message="Error in History::getProb: Radiation could not be"; | |
2572 | message+=" interpreted as FSR or ISR."; | |
2573 | infoPtr->errorMsg(message); | |
2574 | } | |
2575 | ||
2576 | if (showerProb <= 0.) showerProb = 0.; | |
2577 | ||
2578 | //// Different coupling constants for qcd and ew splittings | |
2579 | //if ( state[Emt].colType() != 0 ) { | |
2580 | // AlphaStrong* asFSR = mergingHooksPtr->AlphaS_FSR(); | |
2581 | // double as = (*asFSR).alphaS(91.188*91.188) / (2.*M_PI); | |
2582 | // showerProb *= as; | |
2583 | //} else { | |
2584 | // AlphaEM* aEMFSR = mergingHooksPtr->AlphaEM_FSR(); | |
2585 | // double aEM = (*aEMFSR).alphaEM(91.188*91.188) / (2.*M_PI); | |
2586 | // showerProb *= aEM; | |
2587 | //} | |
2588 | ||
2589 | // Done | |
2590 | return showerProb; | |
2591 | } | |
2592 | ||
2593 | ||
2594 | //-------------------------------------------------------------------------- | |
2595 | ||
2596 | // Set up the beams (fill the beam particles with the correct | |
2597 | // current incoming particles) to allow calculation of splitting | |
2598 | // probability. | |
2599 | // For interleaved evolution, set assignments dividing PDFs into | |
2600 | // sea and valence content. This assignment is, until a history path | |
2601 | // is chosen and a first trial shower performed, not fully correct | |
2602 | // (since content is chosen form too high x and too low scale). The | |
2603 | // assignment used for reweighting will be corrected after trial | |
2604 | // showering | |
2605 | ||
2606 | void History::setupBeams() { | |
2607 | ||
2608 | // Do nothing for empty event, possible if sequence of | |
2609 | // clusterings was ill-advised in that it results in | |
2610 | // colour-disconnected states | |
2611 | if (state.size() < 4) return; | |
2612 | // Do nothing for e+e- beams | |
2613 | if ( state[3].colType() == 0 ) return; | |
2614 | if ( state[4].colType() == 0 ) return; | |
2615 | ||
2616 | // Incoming partons to hard process are stored in slots 3 and 4. | |
2617 | int inS = 0; | |
2618 | int inP = 0; | |
2619 | int inM = 0; | |
2620 | for(int i=0;i< int(state.size()); ++i) { | |
2621 | if (state[i].mother1() == 1) inP = i; | |
2622 | if (state[i].mother1() == 2) inM = i; | |
2623 | } | |
2624 | ||
2625 | // Save some info before clearing beams | |
2626 | // Mothers of incoming partons companion code | |
2627 | int motherPcompRes = -1; | |
2628 | int motherMcompRes = -1; | |
2629 | ||
2630 | bool sameFlavP = false; | |
2631 | bool sameFlavM = false; | |
2632 | ||
2633 | if (mother) { | |
2634 | int inMotherP = 0; | |
2635 | int inMotherM = 0; | |
2636 | for(int i=0;i< int(mother->state.size()); ++i) { | |
2637 | if (mother->state[i].mother1() == 1) inMotherP = i; | |
2638 | if (mother->state[i].mother1() == 2) inMotherM = i; | |
2639 | } | |
2640 | sameFlavP = (state[inP].id() == mother->state[inMotherP].id()); | |
2641 | sameFlavM = (state[inM].id() == mother->state[inMotherM].id()); | |
2642 | ||
2643 | motherPcompRes = (sameFlavP) ? beamA[0].companion() : -2; | |
2644 | motherMcompRes = (sameFlavM) ? beamB[0].companion() : -2; | |
2645 | } | |
2646 | ||
2647 | // Append the current incoming particles to the beam | |
2648 | beamA.clear(); | |
2649 | beamB.clear(); | |
2650 | ||
2651 | // Get energy of incoming particles | |
2652 | double Ep = 2. * state[inP].e(); | |
2653 | double Em = 2. * state[inM].e(); | |
2654 | ||
2655 | // If incoming partons are massive then recalculate to put them massless. | |
2656 | if (state[inP].m() != 0. || state[inM].m() != 0.) { | |
2657 | Ep = state[inP].pPos() + state[inM].pPos(); | |
2658 | Em = state[inP].pNeg() + state[inM].pNeg(); | |
2659 | } | |
2660 | ||
2661 | // Add incoming hard-scattering partons to list in beam remnants. | |
2662 | double x1 = Ep / state[inS].m(); | |
2663 | beamA.append( inP, state[inP].id(), x1); | |
2664 | double x2 = Em / state[inS].m(); | |
2665 | beamB.append( inM, state[inM].id(), x2); | |
2666 | ||
2667 | // Scale. For ME multiplicity history, put scale to mu_F | |
2668 | // (since sea/valence quark content is chosen from this scale) | |
2669 | double scalePDF = (mother) ? scale : infoPtr->QFac(); | |
2670 | // Find whether incoming partons are valence or sea. Store. | |
2671 | // Can I do better, e.g. by setting the scale to the hard process | |
2672 | // scale (= M_W) or by replacing one of the x values by some x/z?? | |
2673 | beamA.xfISR( 0, state[inP].id(), x1, scalePDF*scalePDF); | |
2674 | if (!mother) { | |
2675 | beamA.pickValSeaComp(); | |
2676 | } else { | |
2677 | beamA[0].companion(motherPcompRes); | |
2678 | } | |
2679 | beamB.xfISR( 0, state[inM].id(), x2, scalePDF*scalePDF); | |
2680 | if (!mother) { | |
2681 | beamB.pickValSeaComp(); | |
2682 | } else { | |
2683 | beamB[0].companion(motherMcompRes); | |
2684 | } | |
2685 | ||
2686 | } | |
2687 | ||
2688 | //-------------------------------------------------------------------------- | |
2689 | ||
2690 | // Calculate the PDF ratio used in the argument of the no-emission | |
2691 | // probability | |
2692 | ||
2693 | double History::pdfForSudakov() { | |
2694 | ||
2695 | // Do nothing for e+e- beams | |
2696 | if ( state[3].colType() == 0 ) return 1.0; | |
2697 | if ( state[4].colType() == 0 ) return 1.0; | |
2698 | ||
2699 | // Check if splittings was ISR or FSR | |
2700 | bool FSR = ( mother->state[clusterIn.emittor].isFinal() | |
2701 | && mother->state[clusterIn.recoiler].isFinal()); | |
2702 | bool FSRinRec = ( mother->state[clusterIn.emittor].isFinal() | |
2703 | && !mother->state[clusterIn.recoiler].isFinal()); | |
2704 | ||
2705 | // Done for pure FSR | |
2706 | if (FSR) return 1.0; | |
2707 | ||
2708 | int iInMother = (FSRinRec)? clusterIn.recoiler : clusterIn.emittor; | |
2709 | // Find side of event that was reclustered | |
2710 | int side = ( mother->state[iInMother].pz() > 0 ) ? 1 : -1; | |
2711 | ||
2712 | int inP = 0; | |
2713 | int inM = 0; | |
2714 | for(int i=0;i< int(state.size()); ++i) { | |
2715 | if (state[i].mother1() == 1) inP = i; | |
2716 | if (state[i].mother1() == 2) inM = i; | |
2717 | } | |
2718 | ||
2719 | // Save mother id | |
2720 | int idMother = mother->state[iInMother].id(); | |
2721 | // Find daughter position and id | |
2722 | int iDau = (side == 1) ? inP : inM; | |
2723 | int idDaughter = state[iDau].id(); | |
2724 | // Get mother x value | |
2725 | double xMother = 2. * mother->state[iInMother].e() / mother->state[0].e(); | |
2726 | // Get daughter x value of daughter | |
2727 | double xDaughter = 2.*state[iDau].e() / state[0].e(); // x1 before isr | |
2728 | ||
2729 | // Calculate pdf ratio | |
2730 | double ratio = getPDFratio(side, true, idMother, xMother, scale, | |
2731 | idDaughter, xDaughter, scale); | |
2732 | ||
2733 | // For FSR with incoming recoiler, maximally return 1.0, as | |
2734 | // is done in Pythia::TimeShower. | |
2735 | // For ISR, return ratio | |
2736 | return ( (FSRinRec)? min(1.,ratio) : ratio); | |
2737 | } | |
2738 | ||
2739 | //-------------------------------------------------------------------------- | |
2740 | ||
2741 | // Perform the clustering of the current state and return the | |
2742 | // clustered state. | |
2743 | // IN Clustering : rad,rec,emt triple to be clustered to two partons | |
2744 | // OUT clustered state | |
2745 | ||
2746 | Event History::cluster( const Clustering & inSystem ) { | |
2747 | ||
2748 | // Initialise tags of particles to be changed | |
2749 | int Rad = inSystem.emittor; | |
2750 | int Rec = inSystem.recoiler; | |
2751 | int Emt = inSystem.emitted; | |
2752 | // Initialise eCM,mHat | |
2753 | double eCM = state[0].e(); | |
2754 | // Flags for type of radiation | |
2755 | int radType = state[Rad].isFinal() ? 1 : -1; | |
2756 | int recType = state[Rec].isFinal() ? 1 : -1; | |
2757 | ||
2758 | // Construct the clustered event | |
2759 | Event NewEvent = Event(); | |
2760 | NewEvent.init("(hard process-modified)", particleDataPtr); | |
2761 | NewEvent.clear(); | |
2762 | // Copy all unchanged particles to NewEvent | |
2763 | for (int i = 0; i < state.size(); ++i) | |
2764 | if ( i != Rad && i != Rec && i != Emt ) | |
2765 | NewEvent.append( state[i] ); | |
2766 | ||
2767 | // Copy all the junctions one by one | |
2768 | for (int i = 0; i < state.sizeJunction(); ++i) | |
2769 | NewEvent.appendJunction( state.getJunction(i) ); | |
2770 | // Set particle data table pointer | |
2771 | NewEvent.setPDTPtr(); | |
2772 | // Find an appropriate scale for the hard process | |
2773 | double mu = choseHardScale(state); | |
2774 | // Initialise scales for new event | |
2775 | NewEvent.saveSize(); | |
2776 | NewEvent.saveJunctionSize(); | |
2777 | NewEvent.scale(mu); | |
2778 | NewEvent.scaleSecond(mu); | |
2779 | ||
2780 | // Set properties of radiator/recoiler after the clustering | |
2781 | // Recoiler properties will be unchanged | |
2782 | Particle RecBefore = Particle( state[Rec] ); | |
2783 | RecBefore.daughters(0,0); | |
2784 | // Find flavour of radiator before splitting | |
2785 | int radID = getRadBeforeFlav(Rad, Emt, state); | |
2786 | Particle RadBefore = Particle( state[Rad] ); | |
2787 | RadBefore.id(radID); | |
2788 | RadBefore.daughters(0,0); | |
2789 | // Put dummy values for colours | |
2790 | RadBefore.cols(RecBefore.acol(),RecBefore.col()); | |
2791 | // Put mass for radiator and recoiler | |
2792 | RecBefore.m(particleDataPtr->m0(state[Rec].id())); | |
2793 | RadBefore.m(particleDataPtr->m0(radID)); | |
2794 | ||
2795 | // Construct momenta and colours of clustered particles | |
2796 | // ISR/FSR splittings are treated differently | |
2797 | if ( radType + recType == 2 && state[Emt].idAbs() != 24 ) { | |
2798 | // Clustering of final(rad)/final(rec) dipole splitting | |
2799 | // Get eCM of (rad,rec,emt) triple | |
2800 | Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p(); | |
2801 | double eCMME = sum.mCalc(); | |
2802 | ||
2803 | // Define radiator and recoiler back-to-back in the dipole | |
2804 | // rest frame [=(state[rad]+state[emt])+state[rec] rest frame] | |
2805 | Vec4 Rad4mom; | |
2806 | Vec4 Rec4mom; | |
2807 | Rad4mom.p( 0., 0., 0.5*eCMME, 0.5*eCMME); | |
2808 | Rec4mom.p( 0., 0.,-0.5*eCMME, 0.5*eCMME); | |
2809 | ||
2810 | // Find boost from Rad4mom+Rec4mom rest frame to event cm frame | |
2811 | Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p()); | |
2812 | Vec4 old2 = Vec4(state[Rec].p()); | |
2813 | RotBstMatrix fromCM; | |
2814 | fromCM.fromCMframe(old1, old2); | |
2815 | // Transform momenta | |
2816 | Rad4mom.rotbst(fromCM); | |
2817 | Rec4mom.rotbst(fromCM); | |
2818 | ||
2819 | RadBefore.p(Rad4mom); | |
2820 | RecBefore.p(Rec4mom); | |
2821 | ||
2822 | } else if ( radType + recType == 2 && state[Emt].idAbs() == 24 ) { | |
2823 | // Clustering of final(rad)/final(rec) dipole splitting | |
2824 | // Get eCM of (rad,rec,emt) triple | |
2825 | Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p(); | |
2826 | double eCMME = sum.mCalc(); | |
2827 | ||
2828 | // Define radiator and recoiler back-to-back in the dipole | |
2829 | // rest frame [=(state[rad]+state[emt])+state[rec] rest frame] | |
2830 | Vec4 Rad4mom; | |
2831 | Vec4 Rec4mom; | |
2832 | Rad4mom.p( 0., 0., 0.5*eCMME, 0.5*eCMME); | |
2833 | Rec4mom.p( 0., 0.,-0.5*eCMME, 0.5*eCMME); | |
2834 | ||
2835 | // Find boost from Rad4mom+Rec4mom rest frame to event cm frame | |
2836 | Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p()); | |
2837 | Vec4 old2 = Vec4(state[Rec].p()); | |
2838 | RotBstMatrix fromCM; | |
2839 | fromCM.fromCMframe(old1, old2); | |
2840 | // Transform momenta | |
2841 | Rad4mom.rotbst(fromCM); | |
2842 | Rec4mom.rotbst(fromCM); | |
2843 | ||
2844 | RadBefore.p(Rad4mom); | |
2845 | RecBefore.p(Rec4mom); | |
2846 | ||
2847 | } else if ( radType + recType == 0 ) { | |
2848 | // Clustering of final(rad)/initial(rec) dipole splitting | |
2849 | // Get eCM of (rad,rec,emt) triple | |
2850 | Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p(); | |
2851 | double eCMME = sum.mCalc(); | |
2852 | // Define radiator and recoiler back-to-back in the dipole | |
2853 | // rest frame [=(state[rad]+state[emt])+state[rec] rest frame] | |
2854 | Vec4 Rad4mom; | |
2855 | Vec4 Rec4mom; | |
2856 | Rad4mom.p( 0., 0., 0.5*eCMME, 0.5*eCMME); | |
2857 | Rec4mom.p( 0., 0.,-0.5*eCMME, 0.5*eCMME); | |
2858 | ||
2859 | // Find boost from Rad4mom+Rec4mom rest frame to event cm frame | |
2860 | Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p()); | |
2861 | Vec4 old2 = Vec4(state[Rec].p()); | |
2862 | RotBstMatrix fromCM; | |
2863 | fromCM.fromCMframe(old1, old2); | |
2864 | // Transform momenta | |
2865 | Rad4mom.rotbst(fromCM); | |
2866 | Rec4mom.rotbst(fromCM); | |
2867 | ||
2868 | // Rescale recoiler momentum | |
2869 | Rec4mom = 2.*state[Rec].p() - Rec4mom; | |
2870 | ||
2871 | RadBefore.p(Rad4mom); | |
2872 | RecBefore.p(Rec4mom); | |
2873 | ||
2874 | // Set mass of initial recoiler to zero | |
2875 | RecBefore.m( 0.0 ); | |
2876 | ||
2877 | } else { | |
2878 | // Clustering of initial(rad)/initial(rec) dipole splitting | |
2879 | // We want to cluster: Meaning doing the inverse of a process | |
2880 | // ( pDaughter + pRecoiler -> pOut ) | |
2881 | // ==> ( pMother + pPartner -> pOut' + pSister ) | |
2882 | // produced by an initial state splitting. The matrix element | |
2883 | // provides us with pMother, pPartner, pSister and pOut' | |
2884 | Vec4 pMother( state[Rad].p() ); | |
2885 | Vec4 pSister( state[Emt].p() ); | |
2886 | Vec4 pPartner( state[Rec].p() ); | |
2887 | Vec4 pDaughter( 0.,0.,0.,0. ); | |
2888 | Vec4 pRecoiler( 0.,0.,0.,0. ); | |
2889 | ||
2890 | // Find side that radiates event (mother moving in | |
2891 | // sign * p_z direction) | |
2892 | ||
2893 | int sign = state[Rad].pz() > 0 ? 1 : -1; | |
2894 | ||
2895 | // Find rotation by phi that would have been done for a | |
2896 | // splitting daughter -> mother + sister | |
2897 | double phi = pSister.phi(); | |
2898 | // Find rotation with -phi | |
2899 | RotBstMatrix rot_by_mphi; | |
2900 | rot_by_mphi.rot(0.,-phi); | |
2901 | // Find rotation with +phi | |
2902 | RotBstMatrix rot_by_pphi; | |
2903 | rot_by_pphi.rot(0.,phi); | |
2904 | ||
2905 | // Transform pMother and outgoing momenta | |
2906 | pMother.rotbst( rot_by_mphi ); | |
2907 | pSister.rotbst( rot_by_mphi ); | |
2908 | pPartner.rotbst( rot_by_mphi ); | |
2909 | for(int i=3; i< NewEvent.size(); ++i) | |
2910 | NewEvent[i].rotbst( rot_by_mphi ); | |
2911 | ||
2912 | // Get mother and partner x values | |
2913 | // x1 after isr | |
2914 | double x1 = 2. * pMother.e() / eCM; | |
2915 | // x2 after isr | |
2916 | double x2 = 2. * pPartner.e() / eCM; | |
2917 | ||
2918 | pDaughter.p( pMother - pSister); | |
2919 | pRecoiler.p( pPartner ); | |
2920 | ||
2921 | // Find boost from event cm frame to rest frame of | |
2922 | // of-shell daughter + on-shell recoiler | |
2923 | RotBstMatrix from_CM_to_DR; | |
2924 | if (sign == 1) | |
2925 | from_CM_to_DR.toCMframe(pDaughter, pRecoiler); | |
2926 | else | |
2927 | from_CM_to_DR.toCMframe(pRecoiler, pDaughter); | |
2928 | ||
2929 | // Transform all momenta | |
2930 | pMother.rotbst( from_CM_to_DR ); | |
2931 | pPartner.rotbst( from_CM_to_DR ); | |
2932 | pSister.rotbst( from_CM_to_DR ); | |
2933 | for(int i=3; i< NewEvent.size(); ++i) | |
2934 | NewEvent[i].rotbst( from_CM_to_DR ); | |
2935 | ||
2936 | // Find theta angle between pMother and z-axis and undo | |
2937 | // rotation that would have been done by shower | |
2938 | double theta = pMother.theta(); | |
2939 | if ( pMother.px() < 0. ) theta *= -1.; | |
2940 | if (sign == -1) theta += M_PI; | |
2941 | // Find rotation by +theta | |
2942 | RotBstMatrix rot_by_ptheta; | |
2943 | rot_by_ptheta.rot(theta, 0.); | |
2944 | ||
2945 | // Transform all momenta | |
2946 | pMother.rotbst( rot_by_ptheta ); | |
2947 | pPartner.rotbst( rot_by_ptheta ); | |
2948 | pSister.rotbst( rot_by_ptheta ); | |
2949 | for(int i=3; i< NewEvent.size(); ++i) | |
2950 | NewEvent[i].rotbst( rot_by_ptheta ); | |
2951 | ||
2952 | // Find z of the splitting | |
2953 | Vec4 qDip( pMother - pSister); | |
2954 | Vec4 qAfter(pMother + pPartner); | |
2955 | Vec4 qBefore(qDip + pPartner); | |
2956 | double z = qBefore.m2Calc() / qAfter.m2Calc(); | |
2957 | ||
2958 | // Calculate new e_CM^2 | |
2959 | double x1New = z*x1; // x1 before isr | |
2960 | double x2New = x2; // x2 before isr | |
2961 | double sHat = x1New*x2New*eCM*eCM; | |
2962 | ||
2963 | // Construct daughter and recoiler momenta | |
2964 | pDaughter.p( 0., 0., sign*0.5*sqrt(sHat), 0.5*sqrt(sHat)); | |
2965 | pRecoiler.p( 0., 0., -sign*0.5*sqrt(sHat), 0.5*sqrt(sHat)); | |
2966 | ||
2967 | // Find boost from current (daughter+recoiler rest frame) | |
2968 | // frame to rest frame of daughter+unchanged recoiler to | |
2969 | // recover the old x2 value | |
2970 | RotBstMatrix from_DR_to_CM; | |
2971 | from_DR_to_CM.bst( 0., 0., sign*( x1New - x2New ) / ( x1New + x2New ) ); | |
2972 | ||
2973 | // Correct for momentum mismatch by transforming all momenta | |
2974 | pMother.rotbst( from_DR_to_CM ); | |
2975 | pPartner.rotbst( from_DR_to_CM ); | |
2976 | pSister.rotbst( from_DR_to_CM ); | |
2977 | pDaughter.rotbst( from_DR_to_CM ); | |
2978 | pRecoiler.rotbst( from_DR_to_CM ); | |
2979 | for(int i=3; i< NewEvent.size(); ++i) | |
2980 | NewEvent[i].rotbst( from_DR_to_CM ); | |
2981 | ||
2982 | // Transform pMother and outgoing momenta | |
2983 | pMother.rotbst( rot_by_pphi ); | |
2984 | pPartner.rotbst( rot_by_pphi ); | |
2985 | pSister.rotbst( rot_by_pphi ); | |
2986 | pDaughter.rotbst( rot_by_pphi ); | |
2987 | pRecoiler.rotbst( rot_by_pphi ); | |
2988 | for(int i=3; i< NewEvent.size(); ++i) | |
2989 | NewEvent[i].rotbst( rot_by_pphi ); | |
2990 | ||
2991 | // Set momenta of particles to be attached to new event record | |
2992 | RecBefore.p( pRecoiler ); | |
2993 | RadBefore.p( pDaughter ); | |
2994 | ||
2995 | } | |
2996 | ||
2997 | // Put some dummy production scales for RecBefore, RadBefore | |
2998 | RecBefore.scale(mu); | |
2999 | RadBefore.scale(mu); | |
3000 | RecBefore.setPDTPtr(particleDataPtr); | |
3001 | RadBefore.setPDTPtr(particleDataPtr); | |
3002 | ||
3003 | // Append new recoiler and find new radiator colour | |
3004 | NewEvent.append(RecBefore); | |
3005 | ||
3006 | // Assign the correct colour to re-clustered radiator. | |
3007 | // Keep old radiator colours for electroweak emission. | |
3008 | int emtID = state[Emt].id(); | |
3009 | if ( emtID == 22 || emtID == 23 || abs(emtID) == 24 ) | |
3010 | RadBefore.cols( state[Rad].col(), state[Rad].acol() ); | |
3011 | // For QCD, carefully construct colour. | |
3012 | else if ( !connectRadiator( RadBefore, radType, RecBefore, recType, | |
3013 | NewEvent ) ) { | |
3014 | // Could happen if previous clustering produced several colour | |
3015 | // singlett subsystems in the event | |
3016 | NewEvent.reset(); | |
3017 | return NewEvent; | |
3018 | } | |
3019 | ||
3020 | // Build the clustered event | |
3021 | Event outState = Event(); | |
3022 | outState.init("(hard process-modified)", particleDataPtr); | |
3023 | outState.clear(); | |
3024 | ||
3025 | // Copy system and incoming beam particles to outState | |
3026 | for (int i = 0; i < 3; ++i) | |
3027 | outState.append( NewEvent[i] ); | |
3028 | // Copy all the junctions one by one | |
3029 | for (int i = 0; i < state.sizeJunction(); ++i) | |
3030 | outState.appendJunction( state.getJunction(i) ); | |
3031 | // Set particle data table pointer | |
3032 | outState.setPDTPtr(); | |
3033 | // Initialise scales for new event | |
3034 | outState.saveSize(); | |
3035 | outState.saveJunctionSize(); | |
3036 | outState.scale(mu); | |
3037 | outState.scaleSecond(mu); | |
3038 | bool radAppended = false; | |
3039 | bool recAppended = false; | |
3040 | int size = int(outState.size()); | |
3041 | // Save position of radiator in new event record | |
3042 | int radPos = 0; | |
3043 | // Append first incoming particle | |
3044 | if ( RecBefore.mother1() == 1) { | |
3045 | outState.append( RecBefore ); | |
3046 | recAppended = true; | |
3047 | } else if ( RadBefore.mother1() == 1 ) { | |
3048 | radPos = outState.append( RadBefore ); | |
3049 | radAppended = true; | |
3050 | } else { | |
3051 | // Find second incoming in input event | |
3052 | int in1 = 0; | |
3053 | for(int i=0; i < int(state.size()); ++i) | |
3054 | if (state[i].mother1() == 1) in1 =i; | |
3055 | outState.append( state[in1] ); | |
3056 | size++; | |
3057 | } | |
3058 | // Append second incoming particle | |
3059 | if ( RecBefore.mother1() == 2) { | |
3060 | outState.append( RecBefore ); | |
3061 | recAppended = true; | |
3062 | } else if ( RadBefore.mother1() == 2 ) { | |
3063 | radPos = outState.append( RadBefore ); | |
3064 | radAppended = true; | |
3065 | } else { | |
3066 | // Find second incoming in input event | |
3067 | int in2 = 0; | |
3068 | for(int i=0; i < int(state.size()); ++i) | |
3069 | if (state[i].mother1() == 2) in2 =i; | |
3070 | ||
3071 | outState.append( state[in2] ); | |
3072 | size++; | |
3073 | } | |
3074 | ||
3075 | // Append new recoiler if not done already | |
3076 | if (!recAppended && !RecBefore.isFinal()) { | |
3077 | recAppended = true; | |
3078 | outState.append( RecBefore); | |
3079 | } | |
3080 | // Append new radiator if not done already | |
3081 | if (!radAppended && !RadBefore.isFinal()) { | |
3082 | radAppended = true; | |
3083 | radPos = outState.append( RadBefore); | |
3084 | } | |
3085 | ||
3086 | // Append intermediate particle | |
3087 | // (careful not to append reclustered recoiler) | |
3088 | for (int i = 0; i < int(NewEvent.size()-1); ++i) | |
3089 | if (NewEvent[i].status() == -22) outState.append( NewEvent[i] ); | |
3090 | ||
3091 | if (!recAppended) outState.append(RecBefore); | |
3092 | if (!radAppended) radPos = outState.append(RadBefore); | |
3093 | ||
3094 | // Append final state particles, partons first (not reclustered recoiler) | |
3095 | for(int i = 0; i < int(NewEvent.size()-1); ++i) | |
3096 | if (NewEvent[i].colType() != 0 && NewEvent[i].isFinal()) | |
3097 | outState.append( NewEvent[i] ); | |
3098 | ||
3099 | for(int i = 0; i < int(NewEvent.size()-1); ++i) | |
3100 | if (NewEvent[i].colType() == 0 && NewEvent[i].isFinal()) | |
3101 | outState.append( NewEvent[i]); | |
3102 | ||
3103 | // Find intermediate and respective daughters | |
3104 | vector<int> PosIntermediate; | |
3105 | vector<int> PosDaughter1; | |
3106 | vector<int> PosDaughter2; | |
3107 | for(int i=0; i < int(outState.size()); ++i) | |
3108 | if (outState[i].status() == -22) { | |
3109 | PosIntermediate.push_back(i); | |
3110 | int d1 = outState[i].daughter1(); | |
3111 | int d2 = outState[i].daughter2(); | |
3112 | // Find daughters in output state | |
3113 | int daughter1 = FindParticle( state[d1], outState); | |
3114 | int daughter2 = FindParticle( state[d2], outState); | |
3115 | // If both daughters found, done | |
3116 | // Else put first final particle as first daughter | |
3117 | // and last final particle as second daughter | |
3118 | if (daughter1 > 0) | |
3119 | PosDaughter1.push_back( daughter1); | |
3120 | else { | |
3121 | daughter1 = 0; | |
3122 | while(!outState[daughter1].isFinal() ) daughter1++; | |
3123 | PosDaughter1.push_back( daughter1); | |
3124 | } | |
3125 | if (daughter2 > 0) | |
3126 | PosDaughter2.push_back( daughter2); | |
3127 | else { | |
3128 | daughter2 = outState.size()-1; | |
3129 | while(!outState[daughter2].isFinal() ) daughter2--; | |
3130 | PosDaughter2.push_back( daughter2); | |
3131 | } | |
3132 | } | |
3133 | // Set daughters and mothers | |
3134 | for(int i=0; i < int(PosIntermediate.size()); ++i) { | |
3135 | outState[PosIntermediate[i]].daughters(PosDaughter1[i],PosDaughter2[i]); | |
3136 | outState[PosDaughter1[i]].mother1(PosIntermediate[i]); | |
3137 | outState[PosDaughter2[i]].mother1(PosIntermediate[i]); | |
3138 | } | |
3139 | ||
3140 | // Find range of final state partons | |
3141 | int minParFinal = int(outState.size()); | |
3142 | int maxParFinal = 0; | |
3143 | for(int i=0; i < int(outState.size()); ++i) | |
3144 | if (outState[i].mother1() == 3 && outState[i].mother2() == 4) { | |
3145 | minParFinal = min(i,minParFinal); | |
3146 | maxParFinal = max(i,maxParFinal); | |
3147 | } | |
3148 | ||
3149 | if (minParFinal == maxParFinal) maxParFinal = 0; | |
3150 | outState[3].daughters(minParFinal,maxParFinal); | |
3151 | outState[4].daughters(minParFinal,maxParFinal); | |
3152 | ||
3153 | // Update event properties | |
3154 | outState.saveSize(); | |
3155 | outState.saveJunctionSize(); | |
3156 | ||
3157 | // Almost there... | |
3158 | // If an intermediate coloured parton exists which was directly | |
3159 | // colour connected to the radiator before the splitting, and the | |
3160 | // radiator before and after the splitting had only one colour, problems | |
3161 | // will arise since the colour of the radiator will be changed, whereas | |
3162 | // the intermediate parton still has the old colour. In effect, this | |
3163 | // means that when setting up a event for trial showering, one colour will | |
3164 | // be free. | |
3165 | // Hence, check for an intermediate coloured triplet resonance has been | |
3166 | // colour-connected to the "old" radiator. | |
3167 | // Find resonance | |
3168 | int iColRes = 0; | |
3169 | if ( radType == -1 && state[Rad].colType() == 1) { | |
3170 | // Find resonance connected to initial colour | |
3171 | for(int i=0; i < int(state.size()); ++i) | |
3172 | if ( i != Rad && i != Emt && i != Rec | |
3173 | && state[i].status() == -22 | |
3174 | && state[i].col() == state[Rad].col() ) | |
3175 | iColRes = i; | |
3176 | } else if ( radType == -1 && state[Rad].colType() == -1) { | |
3177 | // Find resonance connected to initial anticolour | |
3178 | for(int i=0; i < int(state.size()); ++i) | |
3179 | if ( i != Rad && i != Emt && i != Rec | |
3180 | && state[i].status() == -22 | |
3181 | && state[i].acol() == state[Rad].acol() ) | |
3182 | iColRes = i; | |
3183 | } else if ( radType == 1 && state[Rad].colType() == 1) { | |
3184 | // Find resonance connected to final state colour | |
3185 | for(int i=0; i < int(state.size()); ++i) | |
3186 | if ( i != Rad && i != Emt && i != Rec | |
3187 | && state[i].status() == -22 | |
3188 | && state[i].col() == state[Rad].col() ) | |
3189 | iColRes = i; | |
3190 | } else if ( radType == 1 && state[Rad].colType() == -1) { | |
3191 | // Find resonance connected to final state anticolour | |
3192 | for(int i=0; i < int(state.size()); ++i) | |
3193 | if ( i != Rad && i != Emt && i != Rec | |
3194 | && state[i].status() == -22 | |
3195 | && state[i].acol() == state[Rad].acol() ) | |
3196 | iColRes = i; | |
3197 | } | |
3198 | ||
3199 | if (iColRes > 0) { | |
3200 | // Now find this resonance in the reclustered state | |
3201 | int iColResNow = FindParticle( state[iColRes], outState); | |
3202 | // Find reclustered radiator colours | |
3203 | int radCol = outState[radPos].col(); | |
3204 | int radAcl = outState[radPos].acol(); | |
3205 | // Find resonance radiator colours | |
3206 | int resCol = outState[iColResNow].col(); | |
3207 | int resAcl = outState[iColResNow].acol(); | |
3208 | // Check if any of the reclustered radiators colours match the resonance | |
3209 | bool matchesRes = (radCol > 0 | |
3210 | && ( radCol == resCol || radCol == resAcl)) | |
3211 | || (radAcl > 0 | |
3212 | && ( radAcl == resCol || radAcl == resAcl)); | |
3213 | ||
3214 | // If a resonance has been found, but no colours match, change | |
3215 | // the colour of the resonance | |
3216 | if (!matchesRes && iColResNow > 0) { | |
3217 | if ( radType == -1 && outState[radPos].colType() == 1) | |
3218 | outState[iColResNow].col(radCol); | |
3219 | else if ( radType ==-1 && outState[radPos].colType() ==-1) | |
3220 | outState[iColResNow].acol(radAcl); | |
3221 | else if ( radType == 1 && outState[radPos].colType() == 1) | |
3222 | outState[iColResNow].col(radCol); | |
3223 | else if ( radType == 1 && outState[radPos].colType() ==-1) | |
3224 | outState[iColResNow].acol(radAcl); | |
3225 | } | |
3226 | ||
3227 | } | |
3228 | ||
3229 | // If event is not constructed properly, return false | |
3230 | if ( !validEvent(outState) ) { | |
3231 | outState.reset(); | |
3232 | return outState; | |
3233 | } | |
3234 | ||
3235 | // Remember position of reclustered radiator in state | |
3236 | iReclusteredNew = radPos; | |
3237 | ||
3238 | // Done | |
3239 | return outState; | |
3240 | } | |
3241 | ||
3242 | //-------------------------------------------------------------------------- | |
3243 | ||
3244 | // Function to get the flavour of the radiator before the splitting | |
3245 | // for clustering | |
3246 | // IN int : Flavour of the radiator after the splitting | |
3247 | // int : Flavour of the emitted after the splitting | |
3248 | // OUT int : Flavour of the radiator before the splitting | |
3249 | ||
3250 | int History::getRadBeforeFlav(const int RadAfter, const int EmtAfter, | |
3251 | const Event& event) { | |
3252 | ||
3253 | int type = event[RadAfter].isFinal() ? 1 :-1; | |
3254 | int emtID = event[EmtAfter].id(); | |
3255 | int radID = event[RadAfter].id(); | |
3256 | int emtCOL = event[EmtAfter].col(); | |
3257 | int radCOL = event[RadAfter].col(); | |
3258 | int emtACL = event[EmtAfter].acol(); | |
3259 | int radACL = event[RadAfter].acol(); | |
3260 | ||
3261 | bool colConnected = ((type == 1) && ( (emtCOL !=0 && (emtCOL ==radACL)) | |
3262 | || (emtACL !=0 && (emtACL ==radCOL)) )) | |
3263 | ||((type ==-1) && ( (emtCOL !=0 && (emtCOL ==radCOL)) | |
3264 | || (emtACL !=0 && (emtACL ==radACL)) )); | |
3265 | // QCD splittings | |
3266 | // Gluon radiation | |
3267 | if ( emtID == 21 ) | |
3268 | return radID; | |
3269 | // Final state gluon splitting | |
3270 | if ( type == 1 && emtID == -radID && !colConnected ) | |
3271 | return 21; | |
3272 | // Initial state s-channel gluon splitting | |
3273 | if ( type ==-1 && radID == 21 ) | |
3274 | return -emtID; | |
3275 | // Initial state t-channel gluon splitting | |
3276 | if ( type ==-1 && emtID != 21 && radID != 21 && !colConnected ) | |
3277 | return 21; | |
3278 | ||
3279 | // Electroweak splittings splittings | |
3280 | // Photon / Z radiation: Calculate invariant mass of system | |
3281 | double m2final = (event[RadAfter].p()+ event[EmtAfter].p()).m2Calc(); | |
3282 | ||
3283 | if ( emtID == 22 || emtID == 23 ) return radID; | |
3284 | // Final state Photon splitting | |
3285 | if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) <= 10. ) | |
3286 | return 22; | |
3287 | // Final state Photon splitting | |
3288 | if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) > 10. ) | |
3289 | return 23; | |
3290 | // Initial state s-channel photon/ Z splitting | |
3291 | if ( type ==-1 && (radID == 22 || radID == 23) ) | |
3292 | return -emtID; | |
3293 | // Initial state t-channel photon / Z splitting: Always bookkeep as photon | |
3294 | if ( type ==-1 && abs(emtID) < 10 && abs(radID) < 10 && colConnected ) | |
3295 | return 22; | |
3296 | ||
3297 | // W+ radiation | |
3298 | // Final state W+ splitting | |
3299 | if ( emtID == 24 && radID < 0 ) return radID + 1; | |
3300 | if ( emtID == 24 && radID > 0 ) return radID + 1; | |
3301 | ||
3302 | // W- radiation | |
3303 | // Final state W- splitting | |
3304 | if ( emtID ==-24 && radID < 0 ) return radID - 1; | |
3305 | if ( emtID ==-24 && radID > 0 ) return radID - 1; | |
3306 | ||
3307 | return 0; | |
3308 | ||
3309 | } | |
3310 | ||
3311 | //-------------------------------------------------------------------------- | |
3312 | ||
3313 | // Function to properly colour-connect the radiator to the rest of | |
3314 | // the event, as needed during clustering | |
3315 | // IN Particle& : Particle to be connected | |
3316 | // Particle : Recoiler forming a dipole with Radiator | |
3317 | // Event : event to which Radiator shall be appended | |
3318 | // OUT true : Radiator could be connected to the event | |
3319 | // false : Radiator could not be connected to the | |
3320 | // event or the resulting event was | |
3321 | // non-valid | |
3322 | ||
3323 | bool History::connectRadiator( Particle& Radiator, const int RadType, | |
3324 | const Particle& Recoiler, const int RecType, | |
3325 | const Event& event ) { | |
3326 | ||
3327 | // Start filling radiator colour indices with dummy values | |
3328 | Radiator.cols( -1, -1 ); | |
3329 | ||
3330 | // Radiator should always be colour-connected to recoiler. | |
3331 | // Three cases (rad = Anti-Quark, Quark, Gluon) to be considered | |
3332 | if ( Radiator.colType() == -1 ) { | |
3333 | // For final state antiquark radiator, the anticolour is fixed | |
3334 | // by the final / initial state recoiler colour / anticolour | |
3335 | if ( RadType + RecType == 2 ) | |
3336 | Radiator.cols( 0, Recoiler.col()); | |
3337 | else if ( RadType + RecType == 0 ) | |
3338 | Radiator.cols( 0, Recoiler.acol()); | |
3339 | // For initial state antiquark radiator, the anticolour is fixed | |
3340 | // by the colour of the emitted gluon (which will be the | |
3341 | // leftover anticolour of a final state particle or the leftover | |
3342 | // colour of an initial state particle ( = the recoiler)) | |
3343 | else { | |
3344 | // Set colour of antiquark radiator to zero | |
3345 | Radiator.col( 0 ); | |
3346 | for (int i = 0; i < event.size(); ++i) { | |
3347 | int col = event[i].col(); | |
3348 | int acl = event[i].acol(); | |
3349 | ||
3350 | if ( event[i].isFinal()) { | |
3351 | // Search for leftover anticolour in final / initial state | |
3352 | if ( acl > 0 && FindCol(acl,i,0,event,1,true) == 0 | |
3353 | && FindCol(acl,i,0,event,2,true) == 0 ) | |
3354 | Radiator.acol(event[i].acol()); | |
3355 | } else { | |
3356 | // Search for leftover colour in initial / final state | |
3357 | if ( col > 0 && FindCol(col,i,0,event,1,true) == 0 | |
3358 | && FindCol(col,i,0,event,2,true) == 0 ) | |
3359 | Radiator.acol(event[i].col()); | |
3360 | } | |
3361 | } // end loop over particles in event record | |
3362 | } | |
3363 | ||
3364 | } else if ( Radiator.colType() == 1 ) { | |
3365 | // For final state quark radiator, the colour is fixed | |
3366 | // by the final / initial state recoiler anticolour / colour | |
3367 | if ( RadType + RecType == 2 ) | |
3368 | Radiator.cols( Recoiler.acol(), 0); | |
3369 | ||
3370 | else if ( RadType + RecType == 0 ) | |
3371 | Radiator.cols( Recoiler.col(), 0); | |
3372 | // For initial state quark radiator, the colour is fixed | |
3373 | // by the anticolour of the emitted gluon (which will be the | |
3374 | // leftover colour of a final state particle or the leftover | |
3375 | // anticolour of an initial state particle ( = the recoiler)) | |
3376 | ||
3377 | else { | |
3378 | // Set anticolour of quark radiator to zero | |
3379 | Radiator.acol( 0 ); | |
3380 | for (int i = 0; i < event.size(); ++i) { | |
3381 | int col = event[i].col(); | |
3382 | int acl = event[i].acol(); | |
3383 | ||
3384 | if ( event[i].isFinal()) { | |
3385 | // Search for leftover colour in final / initial state | |
3386 | if ( col > 0 && FindCol(col,i,0,event,1,true) == 0 | |
3387 | && FindCol(col,i,0,event,2,true) == 0) | |
3388 | Radiator.col(event[i].col()); | |
3389 | } else { | |
3390 | // Search for leftover anticolour in initial / final state | |
3391 | if ( acl > 0 && FindCol(acl,i,0,event,1,true) == 0 | |
3392 | && FindCol(acl,i,0,event,2,true) == 0) | |
3393 | Radiator.col(event[i].acol()); | |
3394 | } | |
3395 | } // end loop over particles in event record | |
3396 | ||
3397 | } // end distinction between fsr / fsr+initial recoiler / isr | |
3398 | ||
3399 | } else if ( Radiator.colType() == 2 ) { | |
3400 | // For a gluon radiator, one (anticolour) colour index is defined | |
3401 | // by the recoiler colour (anticolour). | |
3402 | // The remaining index is chosen to match the free index in the | |
3403 | // event | |
3404 | // Search for leftover colour (anticolour) in the final state | |
3405 | for (int i = 0; i < event.size(); ++i) { | |
3406 | int col = event[i].col(); | |
3407 | int acl = event[i].acol(); | |
3408 | int iEx = i; | |
3409 | ||
3410 | if ( event[i].isFinal()) { | |
3411 | if ( col > 0 && FindCol(col,iEx,0,event,1,true) == 0 | |
3412 | && FindCol(col,iEx,0,event,2,true) == 0) { | |
3413 | if (Radiator.status() < 0 ) Radiator.col(event[i].col()); | |
3414 | else Radiator.acol(event[i].col()); | |
3415 | } | |
3416 | if ( acl > 0 && FindCol(acl,iEx,0,event,2,true) == 0 | |
3417 | && FindCol(acl,iEx,0,event,1,true) == 0 ) { | |
3418 | if (Radiator.status() < 0 ) Radiator.acol(event[i].acol()); | |
3419 | else Radiator.col(event[i].acol()); | |
3420 | } | |
3421 | } else { | |
3422 | if ( col > 0 && FindCol(col,iEx,0,event,1,true) == 0 | |
3423 | && FindCol(col,iEx,0,event,2,true) == 0) { | |
3424 | if (Radiator.status() < 0 ) Radiator.acol(event[i].col()); | |
3425 | else Radiator.col(event[i].col()); | |
3426 | } | |
3427 | if ( acl > 0 && (FindCol(acl,iEx,0,event,2,true) == 0 | |
3428 | && FindCol(acl,iEx,0,event,1,true) == 0)) { | |
3429 | if (Radiator.status() < 0 ) Radiator.col(event[i].acol()); | |
3430 | else Radiator.acol(event[i].acol()); | |
3431 | } | |
3432 | } | |
3433 | } // end loop over particles in event record | |
3434 | } // end cases of different radiator colour type | |
3435 | ||
3436 | // If either colour or anticolour has not been set, return false | |
3437 | if (Radiator.col() < 0 || Radiator.acol() < 0) return false; | |
3438 | // Done | |
3439 | return true; | |
3440 | } | |
3441 | ||
3442 | //-------------------------------------------------------------------------- | |
3443 | ||
3444 | // Function to find a colour (anticolour) index in the input event | |
3445 | // IN int col : Colour tag to be investigated | |
3446 | // int iExclude1 : Identifier of first particle to be excluded | |
3447 | // from search | |
3448 | // int iExclude2 : Identifier of second particle to be excluded | |
3449 | // from search | |
3450 | // Event event : event to be searched for colour tag | |
3451 | // int type : Tag to define if col should be counted as | |
3452 | // colour (type = 1) [->find anti-colour index | |
3453 | // contracted with col] | |
3454 | // anticolour (type = 2) [->find colour index | |
3455 | // contracted with col] | |
3456 | // OUT int : Position of particle in event record | |
3457 | // contraced with col [0 if col is free tag] | |
3458 | ||
3459 | int History::FindCol(int col, int iExclude1, int iExclude2, | |
3460 | const Event& event, int type, bool isHardIn) { | |
3461 | ||
3462 | bool isHard = isHardIn; | |
3463 | int index = 0; | |
3464 | ||
3465 | if (isHard) { | |
3466 | // Search event record for matching colour & anticolour | |
3467 | for(int n = 0; n < event.size(); ++n) { | |
3468 | if ( n != iExclude1 && n != iExclude2 | |
3469 | && event[n].colType() != 0 | |
3470 | &&( event[n].status() > 0 // Check outgoing | |
3471 | || event[n].status() == -21) ) { // Check incoming | |
3472 | if ( event[n].acol() == col ) { | |
3473 | index = -n; | |
3474 | break; | |
3475 | } | |
3476 | if ( event[n].col() == col ) { | |
3477 | index = n; | |
3478 | break; | |
3479 | } | |
3480 | } | |
3481 | } | |
3482 | } else { | |
3483 | ||
3484 | // Search event record for matching colour & anticolour | |
3485 | for(int n = 0; n < event.size(); ++n) { | |
3486 | if ( n != iExclude1 && n != iExclude2 | |
3487 | && event[n].colType() != 0 | |
3488 | &&( event[n].status() == 43 // Check outgoing from ISR | |
3489 | || event[n].status() == 51 // Check outgoing from FSR | |
3490 | || event[n].status() == -41 // first initial | |
3491 | || event[n].status() == -42) ) { // second initial | |
3492 | if ( event[n].acol() == col ) { | |
3493 | index = -n; | |
3494 | break; | |
3495 | } | |
3496 | if ( event[n].col() == col ) { | |
3497 | index = n; | |
3498 | break; | |
3499 | } | |
3500 | } | |
3501 | } | |
3502 | } | |
3503 | // if no matching colour / anticolour has been found, return false | |
3504 | if ( type == 1 && index < 0) return abs(index); | |
3505 | if ( type == 2 && index > 0) return abs(index); | |
3506 | ||
3507 | return 0; | |
3508 | } | |
3509 | ||
3510 | //-------------------------------------------------------------------------- | |
3511 | ||
3512 | // Function to in the input event find a particle with quantum | |
3513 | // numbers matching those of the input particle | |
3514 | // IN Particle : Particle to be searched for | |
3515 | // Event : Event to be searched in | |
3516 | // OUT int : > 0 : Position of matching particle in event | |
3517 | // < 0 : No match in event | |
3518 | ||
3519 | int History::FindParticle( const Particle& particle, const Event& event, | |
3520 | bool checkStatus ) { | |
3521 | ||
3522 | int index = -1; | |
3523 | ||
3524 | for ( int i = int(event.size()) - 1; i > 0; --i ) | |
3525 | if ( event[i].id() == particle.id() | |
3526 | && event[i].colType() == particle.colType() | |
3527 | && event[i].chargeType() == particle.chargeType() | |
3528 | && event[i].col() == particle.col() | |
3529 | && event[i].acol() == particle.acol() | |
3530 | && event[i].charge() == particle.charge() ) { | |
3531 | index = i; | |
3532 | break; | |
3533 | } | |
3534 | ||
3535 | if ( checkStatus && event[index].status() != particle.status() ) | |
3536 | index = -1; | |
3537 | ||
3538 | return index; | |
3539 | } | |
3540 | ||
3541 | //-------------------------------------------------------------------------- | |
3542 | ||
3543 | // Function to get the colour of the radiator before the splitting | |
3544 | // for clustering | |
3545 | // IN int : Position of the radiator after the splitting, in the event | |
3546 | // int : Position of the emitted after the splitting, in the event | |
3547 | // Event : Reference event | |
3548 | // OUT int : Colour of the radiator before the splitting | |
3549 | ||
3550 | int History::getRadBeforeCol(const int rad, const int emt, | |
3551 | const Event& event) { | |
3552 | ||
3553 | // Save type of splitting | |
3554 | int type = (event[rad].isFinal()) ? 1 :-1; | |
3555 | // Get flavour of radiator after potential clustering | |
3556 | int radBeforeFlav = getRadBeforeFlav(rad,emt,event); | |
3557 | // Get colours of the radiator before the potential clustering | |
3558 | int radBeforeCol = -1; | |
3559 | // Get reconstructed gluon colours | |
3560 | if (radBeforeFlav == 21) { | |
3561 | ||
3562 | // Start with quark emissions in FSR | |
3563 | if (type == 1 && event[emt].id() != 21) { | |
3564 | radBeforeCol = (event[rad].col() > 0) | |
3565 | ? event[rad].col() : event[emt].col(); | |
3566 | // Quark emissions in ISR | |
3567 | } else if (type == -1 && event[emt].id() != 21) { | |
3568 | radBeforeCol = (event[rad].col() > 0) | |
3569 | ? event[rad].col() : event[emt].acol(); | |
3570 | //Gluon emissions in FSR | |
3571 | } else if (type == 1 && event[emt].id() == 21) { | |
3572 | // If emitted is a gluon, remove the repeated index, and take | |
3573 | // the remaining indices as colour and anticolour | |
3574 | int colRemove = (event[rad].col() == event[emt].acol()) | |
3575 | ? event[rad].acol() : event[rad].col(); | |
3576 | radBeforeCol = (event[rad].col() == colRemove) | |
3577 | ? event[emt].col() : event[rad].col(); | |
3578 | //Gluon emissions in ISR | |
3579 | } else if (type == -1 && event[emt].id() == 21) { | |
3580 | // If emitted is a gluon, remove the repeated index, and take | |
3581 | // the remaining indices as colour and anticolour | |
3582 | int colRemove = (event[rad].col() == event[emt].col()) | |
3583 | ? event[rad].col() : event[rad].acol(); | |
3584 | radBeforeCol = (event[rad].col() == colRemove) | |
3585 | ? event[emt].acol() : event[rad].col(); | |
3586 | } | |
3587 | ||
3588 | // Get reconstructed quark colours | |
3589 | } else if ( radBeforeFlav != 21 && radBeforeFlav > 0) { | |
3590 | ||
3591 | // Quark emission in FSR | |
3592 | if (type == 1 && event[emt].id() != 21) { | |
3593 | // If radiating is a quark, remove the repeated index, and take | |
3594 | // the remaining indices as colour and anticolour | |
3595 | int colRemove = (event[rad].col() == event[emt].acol()) | |
3596 | ? event[rad].acol() : 0; | |
3597 | radBeforeCol = (event[rad].col() == colRemove) | |
3598 | ? event[emt].col() : event[rad].col(); | |
3599 | //Gluon emissions in FSR | |
3600 | } else if (type == 1 && event[emt].id() == 21) { | |
3601 | // If emitted is a gluon, remove the repeated index, and take | |
3602 | // the remaining indices as colour and anticolour | |
3603 | int colRemove = (event[rad].col() == event[emt].acol()) | |
3604 | ? event[rad].col() : 0; | |
3605 | radBeforeCol = (event[rad].col() == colRemove) | |
3606 | ? event[emt].col() : event[rad].col(); | |
3607 | //Quark emissions in ISR | |
3608 | } else if (type == -1 && event[emt].id() != 21) { | |
3609 | // If emitted is a quark, remove the repeated index, and take | |
3610 | // the remaining indices as colour and anticolour | |
3611 | int colRemove = (event[rad].col() == event[emt].col()) | |
3612 | ? event[rad].col() : 0; | |
3613 | radBeforeCol = (event[rad].col() == colRemove) | |
3614 | ? event[emt].acol() : event[rad].col(); | |
3615 | //Gluon emissions in ISR | |
3616 | } else if (type == -1 && event[emt].id() == 21) { | |
3617 | // If emitted is a gluon, remove the repeated index, and take | |
3618 | // the remaining indices as colour and anticolour | |
3619 | int colRemove = (event[rad].col() == event[emt].col()) | |
3620 | ? event[rad].col() : 0; | |
3621 | radBeforeCol = (event[rad].col() == colRemove) | |
3622 | ? event[emt].acol() : event[rad].col(); | |
3623 | } | |
3624 | // Other particles are assumed uncoloured | |
3625 | } else { | |
3626 | radBeforeCol = 0; | |
3627 | } | |
3628 | ||
3629 | return radBeforeCol; | |
3630 | ||
3631 | } | |
3632 | ||
3633 | //-------------------------------------------------------------------------- | |
3634 | ||
3635 | // Function to get the anticolour of the radiator before the splitting | |
3636 | // for clustering | |
3637 | // IN int : Position of the radiator after the splitting, in the event | |
3638 | // int : Position of the emitted after the splitting, in the event | |
3639 | // Event : Reference event | |
3640 | // OUT int : Anticolour of the radiator before the splitting | |
3641 | ||
3642 | int History::getRadBeforeAcol(const int rad, const int emt, | |
3643 | const Event& event) { | |
3644 | ||
3645 | // Save type of splitting | |
3646 | int type = (event[rad].isFinal()) ? 1 :-1; | |
3647 | // Get flavour of radiator after potential clustering | |
3648 | int radBeforeFlav = getRadBeforeFlav(rad,emt,event); | |
3649 | // Get colours of the radiator before the potential clustering | |
3650 | int radBeforeAcl = -1; | |
3651 | // Get reconstructed gluon colours | |
3652 | if (radBeforeFlav == 21) { | |
3653 | ||
3654 | // Start with quark emissions in FSR | |
3655 | if (type == 1 && event[emt].id() != 21) { | |
3656 | radBeforeAcl = (event[rad].acol() > 0) | |
3657 | ? event[rad].acol() : event[emt].acol(); | |
3658 | // Quark emissions in ISR | |
3659 | } else if (type == -1 && event[emt].id() != 21) { | |
3660 | radBeforeAcl = (event[rad].acol() > 0) | |
3661 | ? event[rad].acol() : event[emt].col(); | |
3662 | //Gluon emissions in FSR | |
3663 | } else if (type == 1 && event[emt].id() == 21) { | |
3664 | // If emitted is a gluon, remove the repeated index, and take | |
3665 | // the remaining indices as colour and anticolour | |
3666 | int colRemove = (event[rad].col() == event[emt].acol()) | |
3667 | ? event[rad].acol() : event[rad].col(); | |
3668 | radBeforeAcl = (event[rad].acol() == colRemove) | |
3669 | ? event[emt].acol() : event[rad].acol(); | |
3670 | //Gluon emissions in ISR | |
3671 | } else if (type == -1 && event[emt].id() == 21) { | |
3672 | // If emitted is a gluon, remove the repeated index, and take | |
3673 | // the remaining indices as colour and anticolour | |
3674 | int colRemove = (event[rad].col() == event[emt].col()) | |
3675 | ? event[rad].col() : event[rad].acol(); | |
3676 | radBeforeAcl = (event[rad].acol() == colRemove) | |
3677 | ? event[emt].col() : event[rad].acol(); | |
3678 | } | |
3679 | ||
3680 | // Get reconstructed anti-quark colours | |
3681 | } else if ( radBeforeFlav != 21 && radBeforeFlav < 0) { | |
3682 | ||
3683 | // Antiquark emission in FSR | |
3684 | if (type == 1 && event[emt].id() != 21) { | |
3685 | // If radiating is a antiquark, remove the repeated index, and take | |
3686 | // the remaining indices as colour and anticolour | |
3687 | int colRemove = (event[rad].col() == event[emt].acol()) | |
3688 | ? event[rad].acol() : 0; | |
3689 | radBeforeAcl = (event[rad].acol() == colRemove) | |
3690 | ? event[emt].acol() : event[rad].acol(); | |
3691 | //Gluon emissions in FSR | |
3692 | } else if (type == 1 && event[emt].id() == 21) { | |
3693 | // If emitted is a gluon, remove the repeated index, and take | |
3694 | // the remaining indices as colour and anticolour | |
3695 | int colRemove = (event[rad].acol() == event[emt].col()) | |
3696 | ? event[rad].acol() : 0; | |
3697 | radBeforeAcl = (event[rad].acol() == colRemove) | |
3698 | ? event[emt].acol() : event[rad].acol(); | |
3699 | //Antiquark emissions in ISR | |
3700 | } else if (type == -1 && event[emt].id() != 21) { | |
3701 | // If emitted is an antiquark, remove the repeated index, and take | |
3702 | // the remaining indices as colour and anticolour | |
3703 | int colRemove = (event[rad].acol() == event[emt].acol()) | |
3704 | ? event[rad].acol() : 0; | |
3705 | radBeforeAcl = (event[rad].acol() == colRemove) | |
3706 | ? event[emt].col() : event[rad].acol(); | |
3707 | //Gluon emissions in ISR | |
3708 | } else if (type == -1 && event[emt].id() == 21) { | |
3709 | // If emitted is a gluon, remove the repeated index, and take | |
3710 | // the remaining indices as colour and anticolour | |
3711 | int colRemove = (event[rad].acol() == event[emt].acol()) | |
3712 | ? event[rad].acol() : 0; | |
3713 | radBeforeAcl = (event[rad].acol() == colRemove) | |
3714 | ? event[emt].col() : event[rad].acol(); | |
3715 | } | |
3716 | // Other particles are considered uncoloured | |
3717 | } else { | |
3718 | radBeforeAcl = 0; | |
3719 | } | |
3720 | ||
3721 | return radBeforeAcl; | |
3722 | ||
3723 | } | |
3724 | ||
3725 | //-------------------------------------------------------------------------- | |
3726 | ||
3727 | // Function to get the parton connected to in by a colour line | |
3728 | // IN int : Position of parton for which partner should be found | |
3729 | // Event : Reference event | |
3730 | // OUT int : If a colour line connects the "in" parton with another | |
3731 | // parton, return the Position of the partner, else return 0 | |
3732 | ||
3733 | int History::getColPartner(const int in, const Event& event) { | |
3734 | ||
3735 | if (event[in].col() == 0) return 0; | |
3736 | ||
3737 | int partner = 0; | |
3738 | // Try to find anticolour index first | |
3739 | partner = FindCol(event[in].col(),in,0,event,1,true); | |
3740 | // If no anticolour index has been found, try colour | |
3741 | if (partner == 0) | |
3742 | partner = FindCol(event[in].col(),in,0,event,2,true); | |
3743 | ||
3744 | return partner; | |
3745 | ||
3746 | } | |
3747 | ||
3748 | //-------------------------------------------------------------------------- | |
3749 | ||
3750 | // Function to get the parton connected to in by an anticolour line | |
3751 | // IN int : Position of parton for which partner should be found | |
3752 | // Event : Reference event | |
3753 | // OUT int : If an anticolour line connects the "in" parton with another | |
3754 | // parton, return the Position of the partner, else return 0 | |
3755 | ||
3756 | int History::getAcolPartner(const int in, const Event& event) { | |
3757 | ||
3758 | if (event[in].acol() == 0) return 0; | |
3759 | ||
3760 | int partner = 0; | |
3761 | // Try to find colour index first | |
3762 | partner = FindCol(event[in].acol(),in,0,event,2,true); | |
3763 | // If no colour index has been found, try anticolour | |
3764 | if (partner == 0) | |
3765 | partner = FindCol(event[in].acol(),in,0,event,1,true); | |
3766 | ||
3767 | return partner; | |
3768 | ||
3769 | } | |
3770 | ||
3771 | //-------------------------------------------------------------------------- | |
3772 | ||
3773 | // Function to get the list of partons connected to the particle | |
3774 | // formed by reclusterinf emt and rad by colour and anticolour lines | |
3775 | // IN int : Position of radiator in the clustering | |
3776 | // IN int : Position of emitted in the clustering | |
3777 | // Event : Reference event | |
3778 | // OUT vector<int> : List of positions of all partons that are connected | |
3779 | // to the parton that will be formed | |
3780 | // by clustering emt and rad. | |
3781 | ||
3782 | vector<int> History::getReclusteredPartners(const int rad, const int emt, | |
3783 | const Event& event) { | |
3784 | ||
3785 | // Save type | |
3786 | int type = event[rad].isFinal() ? 1 : -1; | |
3787 | // Get reclustered colours | |
3788 | int radBeforeCol = getRadBeforeCol(rad, emt, event); | |
3789 | int radBeforeAcl = getRadBeforeAcol(rad, emt, event); | |
3790 | // Declare output | |
3791 | vector<int> partners; | |
3792 | ||
3793 | // Start with FSR clusterings | |
3794 | if (type == 1) { | |
3795 | ||
3796 | for(int i=0; i < int(event.size()); ++i) { | |
3797 | // Check all initial state partons | |
3798 | if ( i != emt && i != rad | |
3799 | && event[i].status() == -21 | |
3800 | && event[i].col() > 0 | |
3801 | && event[i].col() == radBeforeCol) | |
3802 | partners.push_back(i); | |
3803 | // Check all final state partons | |
3804 | if ( i != emt && i != rad | |
3805 | && event[i].isFinal() | |
3806 | && event[i].acol() > 0 | |
3807 | && event[i].acol() == radBeforeCol) | |
3808 | partners.push_back(i); | |
3809 | // Check all initial state partons | |
3810 | if ( i != emt && i != rad | |
3811 | && event[i].status() == -21 | |
3812 | && event[i].acol() > 0 | |
3813 | && event[i].acol() == radBeforeAcl) | |
3814 | partners.push_back(i); | |
3815 | // Check all final state partons | |
3816 | if ( i != emt && i != rad | |
3817 | && event[i].isFinal() | |
3818 | && event[i].col() > 0 | |
3819 | && event[i].col() == radBeforeAcl) | |
3820 | partners.push_back(i); | |
3821 | } | |
3822 | // Start with ISR clusterings | |
3823 | } else { | |
3824 | ||
3825 | for(int i=0; i < int(event.size()); ++i) { | |
3826 | // Check all initial state partons | |
3827 | if ( i != emt && i != rad | |
3828 | && event[i].status() == -21 | |
3829 | && event[i].acol() > 0 | |
3830 | && event[i].acol() == radBeforeCol) | |
3831 | partners.push_back(i); | |
3832 | // Check all final state partons | |
3833 | if ( i != emt && i != rad | |
3834 | && event[i].isFinal() | |
3835 | && event[i].col() > 0 | |
3836 | && event[i].col() == radBeforeCol) | |
3837 | partners.push_back(i); | |
3838 | // Check all initial state partons | |
3839 | if ( i != emt && i != rad | |
3840 | && event[i].status() == -21 | |
3841 | && event[i].col() > 0 | |
3842 | && event[i].col() == radBeforeAcl) | |
3843 | partners.push_back(i); | |
3844 | // Check all final state partons | |
3845 | if ( i != emt && i != rad | |
3846 | && event[i].isFinal() | |
3847 | && event[i].acol() > 0 | |
3848 | && event[i].acol() == radBeforeAcl) | |
3849 | partners.push_back(i); | |
3850 | } | |
3851 | ||
3852 | } | |
3853 | // Done | |
3854 | return partners; | |
3855 | } | |
3856 | ||
3857 | //-------------------------------------------------------------------------- | |
3858 | ||
3859 | // Function to extract a chain of colour-connected partons in | |
3860 | // the event | |
3861 | // IN int : Type of parton from which to start extracting a | |
3862 | // parton chain. If the starting point is a quark | |
3863 | // i.e. flavType = 1, a chain of partons that are | |
3864 | // consecutively connected by colour-lines will be | |
3865 | // extracted. If the starting point is an antiquark | |
3866 | // i.e. flavType =-1, a chain of partons that are | |
3867 | // consecutively connected by anticolour-lines | |
3868 | // will be extracted. | |
3869 | // IN int : Position of the parton from which a | |
3870 | // colour-connected chain should be derived | |
3871 | // IN Event : Refernence event | |
3872 | // IN/OUT vector<int> : Partons that should be excluded from the search. | |
3873 | // OUT vector<int> : Positions of partons along the chain | |
3874 | // OUT bool : Found singlet / did not find singlet | |
3875 | ||
3876 | bool History::getColSinglet( const int flavType, const int iParton, | |
3877 | const Event& event, vector<int>& exclude, vector<int>& colSinglet) { | |
3878 | ||
3879 | // If no possible flavour to start from has been found | |
3880 | if (iParton < 0) return false; | |
3881 | ||
3882 | // If no further partner has been found in a previous iteration, | |
3883 | // and the whole final state has been excluded, we're done | |
3884 | if (iParton == 0) { | |
3885 | ||
3886 | // Count number of final state partons | |
3887 | int nFinal = 0; | |
3888 | for(int i=0; i < int(event.size()); ++i) | |
3889 | if ( event[i].isFinal() && event[i].colType() != 0) | |
3890 | nFinal++; | |
3891 | ||
3892 | // Get number of initial state partons in the list of | |
3893 | // excluded partons | |
3894 | int nExclude = int(exclude.size()); | |
3895 | int nInitExclude = 0; | |
3896 | if (!event[exclude[2]].isFinal()) | |
3897 | nInitExclude++; | |
3898 | if (!event[exclude[3]].isFinal()) | |
3899 | nInitExclude++; | |
3900 | ||
3901 | // If the whole final state has been considered, return | |
3902 | if (nFinal == nExclude - nInitExclude) | |
3903 | return true; | |
3904 | else | |
3905 | return false; | |
3906 | ||
3907 | } | |
3908 | ||
3909 | // Declare colour partner | |
3910 | int colP = 0; | |
3911 | // Save the colour partner | |
3912 | colSinglet.push_back(iParton); | |
3913 | // Remove the partner from the list | |
3914 | exclude.push_back(iParton); | |
3915 | // When starting out from a quark line, follow the colour lines | |
3916 | if (flavType == 1) | |
3917 | colP = getColPartner(iParton,event); | |
3918 | // When starting out from an antiquark line, follow the anticolour lines | |
3919 | else | |
3920 | colP = getAcolPartner(iParton,event); | |
3921 | ||
3922 | // Do not count excluded partons twice | |
3923 | for(int i = 0; i < int(exclude.size()); ++i) | |
3924 | if (colP == exclude[i]) | |
3925 | return true; | |
3926 | ||
3927 | // Recurse | |
3928 | return getColSinglet(flavType,colP,event,exclude,colSinglet); | |
3929 | ||
3930 | } | |
3931 | ||
3932 | //-------------------------------------------------------------------------- | |
3933 | ||
3934 | // Function to check that a set of partons forms a colour singlet | |
3935 | // IN Event : Reference event | |
3936 | // IN vector<int> : Positions of the partons in the set | |
3937 | // OUT bool : Is a colour singlet / is not | |
3938 | ||
3939 | bool History::isColSinglet( const Event& event, | |
3940 | vector<int> system ) { | |
3941 | ||
3942 | // Check if system forms a colour singlet | |
3943 | for(int i=0; i < int(system.size()); ++i ) { | |
3944 | // Match quark and gluon colours | |
3945 | if ( system[i] > 0 | |
3946 | && (event[system[i]].colType() == 1 | |
3947 | || event[system[i]].colType() == 2) ) { | |
3948 | for(int j=0; j < int(system.size()); ++j) | |
3949 | // If flavour matches, remove both partons and continue | |
3950 | if ( system[i] > 0 | |
3951 | && system[j] > 0 | |
3952 | && event[system[i]].col() == event[system[j]].acol()) { | |
3953 | // Remove index and break | |
3954 | system[i] = 0; | |
3955 | system[j] = 0; | |
3956 | break; | |
3957 | } | |
3958 | } | |
3959 | // Match antiquark and gluon anticolours | |
3960 | if ( system[i] > 0 | |
3961 | && (event[system[i]].colType() == -1 | |
3962 | || event[system[i]].colType() == 2) ) { | |
3963 | for(int j=0; j < int(system.size()); ++j) | |
3964 | // If flavour matches, remove both partons and continue | |
3965 | if ( system[i] > 0 | |
3966 | && system[j] > 0 | |
3967 | && event[system[i]].acol() == event[system[j]].col()) { | |
3968 | // Remove index and break | |
3969 | system[i] = 0; | |
3970 | system[j] = 0; | |
3971 | break; | |
3972 | } | |
3973 | } | |
3974 | ||
3975 | } | |
3976 | ||
3977 | // The system is a colour singlet if for all colours, | |
3978 | // an anticolour was found | |
3979 | bool isColSing = true; | |
3980 | for(int i=0; i < int(system.size()); ++i) | |
3981 | if ( system[i] != 0 ) | |
3982 | isColSing = false; | |
3983 | ||
3984 | // Return | |
3985 | return isColSing; | |
3986 | ||
3987 | } | |
3988 | ||
3989 | //-------------------------------------------------------------------------- | |
3990 | ||
3991 | // Function to check that a set of partons forms a flavour singlet | |
3992 | // IN Event : Reference event | |
3993 | // IN vector<int> : Positions of the partons in the set | |
3994 | // IN int : Flavour of all the quarks in the set, if | |
3995 | // all quarks in a set should have a fixed flavour | |
3996 | // OUT bool : Is a flavour singlet / is not | |
3997 | ||
3998 | bool History::isFlavSinglet( const Event& event, | |
3999 | vector<int> system, int flav) { | |
4000 | ||
4001 | // If a decoupled colour singlet has been found, check if this is also | |
4002 | // a flavour singlet | |
4003 | // Check that each quark matches an antiquark | |
4004 | for(int i=0; i < int(system.size()); ++i) | |
4005 | if ( system[i] > 0 ) { | |
4006 | for(int j=0; j < int(system.size()); ++j) { | |
4007 | // If flavour of outgoing partons matches, | |
4008 | // remove both partons and continue. | |
4009 | // Skip all bosons | |
4010 | if ( event[i].idAbs() != 21 | |
4011 | && event[i].idAbs() != 22 | |
4012 | && event[i].idAbs() != 23 | |
4013 | && event[i].idAbs() != 24 | |
4014 | && system[i] > 0 | |
4015 | && system[j] > 0 | |
4016 | && event[system[i]].isFinal() | |
4017 | && event[system[j]].isFinal() | |
4018 | && event[system[i]].id() == -1*event[system[j]].id()) { | |
4019 | // If we want to check if only one flavour of quarks | |
4020 | // exists | |
4021 | if (abs(flav) > 0 && event[system[i]].idAbs() != flav) | |
4022 | return false; | |
4023 | // Remove index and break | |
4024 | system[i] = 0; | |
4025 | system[j] = 0; | |
4026 | break; | |
4027 | } | |
4028 | // If flavour of outgoing and incoming partons match, | |
4029 | // remove both partons and continue. | |
4030 | // Skip all bosons | |
4031 | if ( event[i].idAbs() != 21 | |
4032 | && event[i].idAbs() != 22 | |
4033 | && event[i].idAbs() != 23 | |
4034 | && event[i].idAbs() != 24 | |
4035 | && system[i] > 0 | |
4036 | && system[j] > 0 | |
4037 | && ( ( !event[system[i]].isFinal() && event[system[j]].isFinal()) | |
4038 | ||( !event[system[j]].isFinal() && event[system[i]].isFinal()) ) | |
4039 | && event[system[i]].id() == event[system[j]].id()) { | |
4040 | // If we want to check if only one flavour of quarks | |
4041 | // exists | |
4042 | if (abs(flav) > 0 && event[system[i]].idAbs() != flav) | |
4043 | return false; | |
4044 | // Remove index and break | |
4045 | system[i] = 0; | |
4046 | system[j] = 0; | |
4047 | break; | |
4048 | } | |
4049 | ||
4050 | } | |
4051 | } | |
4052 | ||
4053 | // The colour singlet is a flavour singlet if for all quarks, | |
4054 | // an antiquark was found | |
4055 | bool isFlavSing = true; | |
4056 | for(int i=0; i < int(system.size()); ++i) | |
4057 | if ( system[i] != 0 ) | |
4058 | isFlavSing = false; | |
4059 | ||
4060 | // Return | |
4061 | return isFlavSing; | |
4062 | ||
4063 | } | |
4064 | ||
4065 | ||
4066 | //-------------------------------------------------------------------------- | |
4067 | ||
4068 | // Function to check if rad,emt,rec triple is allowed for clustering | |
4069 | // IN int rad,emt,rec : Positions (in event record) of the three | |
4070 | // particles considered for clustering | |
4071 | // Event event : Reference event | |
4072 | ||
4073 | bool History::allowedClustering( int rad, int emt, int rec, int partner, | |
4074 | const Event& event ) { | |
4075 | ||
4076 | // Declare output | |
4077 | bool allowed = true; | |
4078 | ||
4079 | // CONSTRUCT SOME PROPERTIES FOR LATER INVESTIGATION | |
4080 | ||
4081 | // Check if the triple forms a colour singlett | |
4082 | bool isSing = isSinglett(rad,emt,partner,event); | |
4083 | int type = (event[rad].isFinal()) ? 1 :-1; | |
4084 | // Get flavour of radiator after potential clustering | |
4085 | int radBeforeFlav = getRadBeforeFlav(rad,emt,event); | |
4086 | // Get colours of the radiator before the potential clustering | |
4087 | int radBeforeCol = getRadBeforeCol(rad,emt,event); | |
4088 | int radBeforeAcl = getRadBeforeAcol(rad,emt,event); | |
4089 | // Get colour partner of reclustered parton | |
4090 | vector<int> radBeforeColP = getReclusteredPartners(rad, emt, event); | |
4091 | ||
4092 | // Count coloured partons in hard process | |
4093 | int nPartonInHard = 0; | |
4094 | for(int i=0; i < int(event.size()); ++i) | |
4095 | // Check all final state partons | |
4096 | if ( event[i].isFinal() | |
4097 | && event[i].colType() != 0 | |
4098 | && mergingHooksPtr->hardProcess.matchesAnyOutgoing(i, event) ) | |
4099 | nPartonInHard++; | |
4100 | ||
4101 | ||
4102 | // Count coloured final state partons in event, excluding | |
4103 | // rad, rec, emt and hard process | |
4104 | int nPartons = 0; | |
4105 | for(int i=0; i < int(event.size()); ++i) | |
4106 | // Check all final state partons | |
4107 | if ( i!=emt && i!=rad && i!=rec | |
4108 | && event[i].isFinal() | |
4109 | && event[i].colType() != 0 | |
4110 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i, event) ) | |
4111 | nPartons++; | |
4112 | ||
4113 | // Count number of initial state partons | |
4114 | int nInitialPartons = 0; | |
4115 | for(int i=0; i < int(event.size()); ++i) | |
4116 | if ( event[i].status() == -21 | |
4117 | && event[i].colType() != 0 ) | |
4118 | nInitialPartons++; | |
4119 | ||
4120 | // Get number of non-charged final state particles | |
4121 | int nFinalEW = 0; | |
4122 | for(int i=0; i < int(event.size()); ++i) | |
4123 | if ( event[i].isFinal() | |
4124 | &&( event[i].id() == 22 | |
4125 | || event[i].id() == 23 | |
4126 | || event[i].id() == 24 | |
4127 | ||(event[i].idAbs() > 10 && event[i].idAbs() < 20))) | |
4128 | nFinalEW++; | |
4129 | ||
4130 | // Check if event after potential clustering contains an even | |
4131 | // number of quarks and/or antiquarks | |
4132 | // (otherwise no electroweak vertex could be formed!) | |
4133 | // Get number of final quarks | |
4134 | int nFinalQuark = 0; | |
4135 | // Get number of excluded final state quarks as well | |
4136 | int nFinalQuarkExc = 0; | |
4137 | for(int i=0; i < int(event.size()); ++i) { | |
4138 | if (i !=rad && i != emt && i != rec) { | |
4139 | if (event[i].isFinal() && event[i].isQuark() ) { | |
4140 | if ( !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,event) ) | |
4141 | nFinalQuark++; | |
4142 | else | |
4143 | nFinalQuarkExc++; | |
4144 | } | |
4145 | } | |
4146 | } | |
4147 | ||
4148 | // Add recoiler to number of final quarks | |
4149 | if (event[rec].isFinal() && event[rec].isQuark()) nFinalQuark++; | |
4150 | // Add radiator after clustering to number of final quarks | |
4151 | if (event[rad].isFinal() && abs(radBeforeFlav) < 10) nFinalQuark++; | |
4152 | ||
4153 | // Get number of initial quarks | |
4154 | int nInitialQuark = 0; | |
4155 | if (type == 1) { | |
4156 | if (event[rec].isFinal()) { | |
4157 | if (event[3].isQuark()) nInitialQuark++; | |
4158 | if (event[4].isQuark()) nInitialQuark++; | |
4159 | } else { | |
4160 | int iOtherIn = (rec == 3) ? 4 : 3; | |
4161 | if (event[rec].isQuark()) nInitialQuark++; | |
4162 | if (event[iOtherIn].isQuark()) nInitialQuark++; | |
4163 | } | |
4164 | } else { | |
4165 | // Add recoiler to number of initial quarks | |
4166 | if (event[rec].isQuark()) nInitialQuark++; | |
4167 | // Add radiator after clustering to number of initial quarks | |
4168 | if (abs(radBeforeFlav) < 10) nInitialQuark++; | |
4169 | } | |
4170 | ||
4171 | // If only odd number of quarks in state, | |
4172 | // reject (no electroweak vertex can be formed). | |
4173 | // This is only true of no primary partonic resonance could have produced | |
4174 | // electroweak bosons. | |
4175 | // Check for tops | |
4176 | int nTop = 0; | |
4177 | for(int i=0; i < int(event.size()); ++i) | |
4178 | if (event[i].idAbs() == 6) | |
4179 | nTop++; | |
4180 | ||
4181 | // BEGIN CHECKING THE CLUSTERING | |
4182 | ||
4183 | // Check if colour is conserved | |
4184 | vector<int> unmatchedCol; | |
4185 | vector<int> unmatchedAcl; | |
4186 | // Check all unmatched colours | |
4187 | for ( int i = 0; i < event.size(); ++i) | |
4188 | if ( i != emt && i != rad | |
4189 | && (event[i].isFinal() || event[i].status() == -21) | |
4190 | && event[i].colType() != 0 ) { | |
4191 | ||
4192 | int colP = getColPartner(i,event); | |
4193 | int aclP = getAcolPartner(i,event); | |
4194 | ||
4195 | if (event[i].col() > 0 | |
4196 | && (colP == emt || colP == rad || colP == 0) ) | |
4197 | unmatchedCol.push_back(i); | |
4198 | if (event[i].acol() > 0 | |
4199 | && (aclP == emt || aclP == rad || aclP == 0) ) | |
4200 | unmatchedAcl.push_back(i); | |
4201 | ||
4202 | } | |
4203 | ||
4204 | // If more than one colour or more than one anticolour are unmatched, | |
4205 | // there is no way to make this clustering work | |
4206 | if (int(unmatchedCol.size()) + int(unmatchedAcl.size()) > 2) | |
4207 | return false; | |
4208 | ||
4209 | // If triple forms colour singlett, check that resulting state | |
4210 | // matches hard core process | |
4211 | if (isSing) | |
4212 | allowed = false; | |
4213 | if (isSing && (abs(radBeforeFlav)<10 && event[rec].isQuark()) ) | |
4214 | allowed = true; | |
4215 | ||
4216 | // Never recluster any outgoing partons of the core V -> qqbar' splitting! | |
4217 | if ( mergingHooksPtr->hardProcess.matchesAnyOutgoing(emt,event) ) | |
4218 | allowed = false; | |
4219 | ||
4220 | // Never allow clustering of any outgoing partons of the hard process | |
4221 | // which would change the flavour of one of the hard process partons! | |
4222 | if ( mergingHooksPtr->hardProcess.matchesAnyOutgoing(rad,event) | |
4223 | && event[rad].id() != radBeforeFlav ) | |
4224 | allowed = false; | |
4225 | ||
4226 | // If only gluons in initial state and no quarks in final state, | |
4227 | // reject (no electroweak vertex can be formed) | |
4228 | if (nFinalEW != 0 && nInitialQuark == 0 | |
4229 | && nFinalQuark == 0 && nFinalQuarkExc == 0) | |
4230 | allowed = false; | |
4231 | ||
4232 | if ( (nInitialQuark + nFinalQuark + nFinalQuarkExc)%2 != 0 ) | |
4233 | allowed = false; | |
4234 | ||
4235 | // Do not allow final state splitting that produces only | |
4236 | // allowed final state gluons, and has a colour-connected initial state | |
4237 | // This means forbidding clusterings that do not allow for a | |
4238 | // t-channel gluon, which is needed to have a quark-antiquark initial state. | |
4239 | // Here, partons excluded from clustering are not counted as possible | |
4240 | // partners to form a t-channel gluon | |
4241 | if (event[3].col() == event[4].acol() | |
4242 | && event[3].acol() == event[4].col() | |
4243 | && nFinalQuark == 0) | |
4244 | allowed = false; | |
4245 | ||
4246 | // No problems with gluon radiation | |
4247 | if (event[emt].id() == 21) | |
4248 | return allowed; | |
4249 | ||
4250 | // Save all hard process candidates | |
4251 | vector<int> outgoingParticles; | |
4252 | int nOut1 = int(mergingHooksPtr->hardProcess.PosOutgoing1.size()); | |
4253 | for ( int i=0; i < nOut1; ++i ) { | |
4254 | int iPos = mergingHooksPtr->hardProcess.PosOutgoing1[i]; | |
4255 | outgoingParticles.push_back( | |
4256 | mergingHooksPtr->hardProcess.state[iPos].id() ); | |
4257 | } | |
4258 | int nOut2 = int(mergingHooksPtr->hardProcess.PosOutgoing2.size()); | |
4259 | for ( int i=0; i < nOut2; ++i ) { | |
4260 | int iPos = mergingHooksPtr->hardProcess.PosOutgoing2[i]; | |
4261 | outgoingParticles.push_back( | |
4262 | mergingHooksPtr->hardProcess.state[iPos].id() ); | |
4263 | } | |
4264 | ||
4265 | // Start more involved checks. g -> q_1 qbar_1 splittings are | |
4266 | // particularly problematic if more than one quark of the emitted | |
4267 | // flavour is present. | |
4268 | // Count number of initial quarks of radiator or emitted flavour | |
4269 | vector<int> iInQuarkFlav; | |
4270 | for(int i=0; i < int(event.size()); ++i) | |
4271 | // Check all initial state partons | |
4272 | if ( i != emt && i != rad | |
4273 | && event[i].status() == -21 | |
4274 | && event[i].idAbs() == event[emt].idAbs() ) | |
4275 | iInQuarkFlav.push_back(i); | |
4276 | ||
4277 | // Count number of final quarks of radiator or emitted flavour | |
4278 | vector<int> iOutQuarkFlav; | |
4279 | for(int i=0; i < int(event.size()); ++i) | |
4280 | // Check all final state partons | |
4281 | if ( i != emt && i != rad | |
4282 | && event[i].isFinal() | |
4283 | && event[i].idAbs() == event[emt].idAbs() ) { | |
4284 | ||
4285 | // Loop through final state hard particles. If one matches, remove the | |
4286 | // matching one, and do not count. | |
4287 | bool matchOut = false; | |
4288 | for (int j = 0; j < int(outgoingParticles.size()); ++j) | |
4289 | if ( event[i].idAbs() == abs(outgoingParticles[j])) { | |
4290 | matchOut = true; | |
4291 | outgoingParticles[j] = 99; | |
4292 | } | |
4293 | if (!matchOut) iOutQuarkFlav.push_back(i); | |
4294 | ||
4295 | } | |
4296 | ||
4297 | // Save number of potentially dangerous quarks | |
4298 | int nInQuarkFlav = int(iInQuarkFlav.size()); | |
4299 | int nOutQuarkFlav = int(iOutQuarkFlav.size()); | |
4300 | ||
4301 | // Easiest problem 0: | |
4302 | // Radiator before splitting exactly matches the partner | |
4303 | // after the splitting | |
4304 | if ( event[partner].isFinal() | |
4305 | && event[partner].id() == 21 | |
4306 | && radBeforeFlav == 21 | |
4307 | && event[partner].col() == radBeforeCol | |
4308 | && event[partner].acol() == radBeforeAcl) | |
4309 | return false; | |
4310 | ||
4311 | // If there are no ambiguities in qqbar pairs, return | |
4312 | if (nInQuarkFlav + nOutQuarkFlav == 0) | |
4313 | return allowed; | |
4314 | ||
4315 | ||
4316 | // Save all quarks and gluons that will not change colour | |
4317 | vector<int> gluon; | |
4318 | vector<int> quark; | |
4319 | vector<int> antiq; | |
4320 | vector<int> partons; | |
4321 | for(int i=0; i < int(event.size()); ++i) | |
4322 | // Check initial and final state partons | |
4323 | if ( i!=emt && i!=rad | |
4324 | && event[i].colType() != 0 | |
4325 | && (event[i].isFinal() || event[i].status() == -21) ) { | |
4326 | // Save index | |
4327 | partons.push_back(i); | |
4328 | // Split into components | |
4329 | if (event[i].colType() == 2) | |
4330 | gluon.push_back(i); | |
4331 | else if (event[i].colType() == 1) | |
4332 | quark.push_back(i); | |
4333 | else if (event[i].colType() == -1) | |
4334 | antiq.push_back(i); | |
4335 | } | |
4336 | ||
4337 | // We split up the test of the g->qq splitting into final state | |
4338 | // and initial state problems | |
4339 | bool isFSRg2qq = ((type == 1) && (event[rad].id() == -1*event[emt].id()) ); | |
4340 | bool isISRg2qq = ((type ==-1) && (event[rad].id() == event[emt].id()) ); | |
4341 | ||
4342 | // First check general things about colour connections | |
4343 | // Check that clustering does not produce a gluon that is exactly | |
4344 | // matched in the final state, or does not have any colour connections | |
4345 | if ( (isFSRg2qq || isISRg2qq) | |
4346 | && int(quark.size()) + int(antiq.size()) | |
4347 | + int(gluon.size()) > nPartonInHard ) { | |
4348 | ||
4349 | vector<int> colours; | |
4350 | vector<int> anticolours; | |
4351 | // Add the colour and anticolour of the gluon before the emission | |
4352 | // to the list, bookkeep initial colour as final anticolour, and | |
4353 | // initial anticolour as final colour | |
4354 | if (type == 1) { | |
4355 | colours.push_back(radBeforeCol); | |
4356 | anticolours.push_back(radBeforeAcl); | |
4357 | } else { | |
4358 | colours.push_back(radBeforeAcl); | |
4359 | anticolours.push_back(radBeforeCol); | |
4360 | } | |
4361 | // Now store gluon colours and anticolours. | |
4362 | for(int i=0; i < int(gluon.size()); ++i) | |
4363 | if (event[gluon[i]].isFinal()) { | |
4364 | colours.push_back(event[gluon[i]].col()); | |
4365 | anticolours.push_back(event[gluon[i]].acol()); | |
4366 | } else { | |
4367 | colours.push_back(event[gluon[i]].acol()); | |
4368 | anticolours.push_back(event[gluon[i]].col()); | |
4369 | } | |
4370 | ||
4371 | // Loop through colours and check if any match with | |
4372 | // anticolours. If colour matches, remove from list | |
4373 | for(int i=0; i < int(colours.size()); ++i) | |
4374 | for(int j=0; j < int(anticolours.size()); ++j) | |
4375 | if (colours[i] > 0 && anticolours[j] > 0 | |
4376 | && colours[i] == anticolours[j]) { | |
4377 | colours[i] = 0; | |
4378 | anticolours[j] = 0; | |
4379 | } | |
4380 | ||
4381 | ||
4382 | // If all gluon anticolours and all colours matched, disallow | |
4383 | // the clustering | |
4384 | bool allMatched = true; | |
4385 | for(int i=0; i < int(colours.size()); ++i) | |
4386 | if (colours[i] != 0) | |
4387 | allMatched = false; | |
4388 | for(int i=0; i < int(anticolours.size()); ++i) | |
4389 | if (anticolours[i] != 0) | |
4390 | allMatched = false; | |
4391 | ||
4392 | if (allMatched) | |
4393 | return false; | |
4394 | ||
4395 | // Now add the colours of the hard process, and check if all | |
4396 | // colours match. | |
4397 | for(int i=0; i < int(quark.size()); ++i) | |
4398 | if ( event[quark[i]].isFinal() | |
4399 | && mergingHooksPtr->hardProcess.matchesAnyOutgoing(quark[i], event) ) | |
4400 | colours.push_back(event[quark[i]].col()); | |
4401 | ||
4402 | for(int i=0; i < int(antiq.size()); ++i) | |
4403 | if ( event[antiq[i]].isFinal() | |
4404 | && mergingHooksPtr->hardProcess.matchesAnyOutgoing(antiq[i], event) ) | |
4405 | anticolours.push_back(event[antiq[i]].acol()); | |
4406 | ||
4407 | // Loop through colours again and check if any match with | |
4408 | // anticolours. If colour matches, remove from list | |
4409 | for(int i=0; i < int(colours.size()); ++i) | |
4410 | ||
4411 | for(int j=0; j < int(anticolours.size()); ++j) | |
4412 | if (colours[i] > 0 && anticolours[j] > 0 | |
4413 | && colours[i] == anticolours[j]) { | |
4414 | colours[i] = 0; | |
4415 | anticolours[j] = 0; | |
4416 | } | |
4417 | ||
4418 | // Check if clustering would produce the hard process | |
4419 | int nNotInHard = 0; | |
4420 | for ( int i=0; i < int(quark.size()); ++i ) | |
4421 | if ( !mergingHooksPtr->hardProcess.matchesAnyOutgoing( quark[i], | |
4422 | event) ) | |
4423 | nNotInHard++; | |
4424 | for ( int i=0; i < int(antiq.size()); ++i ) | |
4425 | if ( !mergingHooksPtr->hardProcess.matchesAnyOutgoing( antiq[i], | |
4426 | event) ) | |
4427 | nNotInHard++; | |
4428 | for(int i=0; i < int(gluon.size()); ++i) | |
4429 | if ( event[gluon[i]].isFinal() ) | |
4430 | nNotInHard++; | |
4431 | if ( type == 1 ) | |
4432 | nNotInHard++; | |
4433 | ||
4434 | // If all colours are matched now, and since we have more quarks than | |
4435 | // present in the hard process, disallow the clustering | |
4436 | allMatched = true; | |
4437 | for(int i=0; i < int(colours.size()); ++i) | |
4438 | if (colours[i] != 0) | |
4439 | allMatched = false; | |
4440 | for(int i=0; i < int(anticolours.size()); ++i) | |
4441 | if (anticolours[i] != 0) | |
4442 | allMatched = false; | |
4443 | ||
4444 | if (allMatched && nNotInHard > 0) | |
4445 | return false; | |
4446 | ||
4447 | } | |
4448 | ||
4449 | // FSR PROBLEMS | |
4450 | ||
4451 | if (isFSRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) { | |
4452 | ||
4453 | // Easiest problem 1: | |
4454 | // RECLUSTERED FINAL STATE GLUON MATCHES INITIAL STATE GLUON | |
4455 | for(int i=0; i < int(gluon.size()); ++i) { | |
4456 | if (!event[gluon[i]].isFinal() | |
4457 | && event[gluon[i]].col() == radBeforeCol | |
4458 | && event[gluon[i]].acol() == radBeforeAcl) | |
4459 | return false; | |
4460 | } | |
4461 | ||
4462 | // Easiest problem 2: | |
4463 | // RECLUSTERED FINAL STATE GLUON MATCHES FINAL STATE GLUON | |
4464 | for(int i=0; i < int(gluon.size()); ++i) { | |
4465 | if (event[gluon[i]].isFinal() | |
4466 | && event[gluon[i]].col() == radBeforeAcl | |
4467 | && event[gluon[i]].acol() == radBeforeCol) | |
4468 | return false; | |
4469 | } | |
4470 | ||
4471 | // Easiest problem 3: | |
4472 | // RECLUSTERED FINAL STATE GLUON MATCHES FINAL STATE Q-QBAR PAIR | |
4473 | if ( int(radBeforeColP.size()) == 2 | |
4474 | && event[radBeforeColP[0]].isFinal() | |
4475 | && event[radBeforeColP[1]].isFinal() | |
4476 | && event[radBeforeColP[0]].id() == -1*event[radBeforeColP[1]].id() ) { | |
4477 | ||
4478 | // This clustering is allowed if there is no colour in the | |
4479 | // initial state | |
4480 | if (nInitialPartons > 0) | |
4481 | return false; | |
4482 | } | |
4483 | ||
4484 | // Next-to-easiest problem 1: | |
4485 | // RECLUSTERED FINAL STATE GLUON MATCHES ONE FINAL STARE Q_1 | |
4486 | // AND ONE INITIAL STATE Q_1 | |
4487 | if ( int(radBeforeColP.size()) == 2 | |
4488 | && (( event[radBeforeColP[0]].status() == -21 | |
4489 | && event[radBeforeColP[1]].isFinal()) | |
4490 | ||( event[radBeforeColP[0]].isFinal() | |
4491 | && event[radBeforeColP[1]].status() == -21)) | |
4492 | && event[radBeforeColP[0]].id() == event[radBeforeColP[1]].id() ) { | |
4493 | ||
4494 | // In principle, clustering this splitting can disconnect | |
4495 | // the colour lines of a graph. However, the colours can be connected | |
4496 | // again if a final or initial partons of the correct flavour exists. | |
4497 | ||
4498 | // Check which of the partners are final / initial | |
4499 | int incoming = (event[radBeforeColP[0]].isFinal()) | |
4500 | ? radBeforeColP[1] : radBeforeColP[0]; | |
4501 | int outgoing = (event[radBeforeColP[0]].isFinal()) | |
4502 | ? radBeforeColP[0] : radBeforeColP[1]; | |
4503 | ||
4504 | // Loop through event to find "recovery partons" | |
4505 | bool clusPossible = false; | |
4506 | for(int i=0; i < int(event.size()); ++i) | |
4507 | if ( i != emt && i != rad | |
4508 | && i != incoming && i != outgoing | |
4509 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,event) ) { | |
4510 | // Check if an incoming parton matches | |
4511 | if ( event[i].status() == -21 | |
4512 | && (event[i].id() == event[outgoing].id() | |
4513 | ||event[i].id() == -1*event[incoming].id()) ) | |
4514 | clusPossible = true; | |
4515 | // Check if a final parton matches | |
4516 | if ( event[i].isFinal() | |
4517 | && (event[i].id() == -1*event[outgoing].id() | |
4518 | ||event[i].id() == event[incoming].id()) ) | |
4519 | clusPossible = true; | |
4520 | } | |
4521 | ||
4522 | // There can be a further complication: If e.g. in | |
4523 | // t-channel photon exchange topologies, both incoming | |
4524 | // partons are quarks, and form colour singlets with any | |
4525 | // number of final state partons, at least try to | |
4526 | // recluster as much as possible. | |
4527 | // For this, check if the incoming parton | |
4528 | // connected to the radiator is connected to a | |
4529 | // colour and flavour singlet | |
4530 | vector<int> excludeIn1; | |
4531 | for(int i=0; i < 4; ++i) | |
4532 | excludeIn1.push_back(0); | |
4533 | vector<int> colSingletIn1; | |
4534 | int flavIn1Type = (event[incoming].id() > 0) ? 1 : -1; | |
4535 | // Try finding colour singlets | |
4536 | bool isColSingIn1 = getColSinglet(flavIn1Type,incoming,event, | |
4537 | excludeIn1,colSingletIn1); | |
4538 | // Check if colour singlet also is a flavour singlet | |
4539 | bool isFlavSingIn1 = isFlavSinglet(event,colSingletIn1); | |
4540 | ||
4541 | // Check if the incoming parton not | |
4542 | // connected to the radiator is connected to a | |
4543 | // colour and flavour singlet | |
4544 | int incoming2 = (incoming == 3) ? 4 : 3; | |
4545 | vector<int> excludeIn2; | |
4546 | for(int i=0; i < 4; ++i) | |
4547 | excludeIn2.push_back(0); | |
4548 | vector<int> colSingletIn2; | |
4549 | int flavIn2Type = (event[incoming2].id() > 0) ? 1 : -1; | |
4550 | // Try finding colour singlets | |
4551 | bool isColSingIn2 = getColSinglet(flavIn2Type,incoming2,event, | |
4552 | excludeIn2,colSingletIn2); | |
4553 | // Check if colour singlet also is a flavour singlet | |
4554 | bool isFlavSingIn2 = isFlavSinglet(event,colSingletIn2); | |
4555 | ||
4556 | // If no "recovery clustering" is possible, reject clustering | |
4557 | if (!clusPossible | |
4558 | && (!isColSingIn1 || !isFlavSingIn1 | |
4559 | || !isColSingIn2 || !isFlavSingIn2)) | |
4560 | return false; | |
4561 | ||
4562 | } | |
4563 | ||
4564 | // Next-to-easiest problem 2: | |
4565 | // FINAL STATE Q-QBAR CLUSTERING DISCONNECTS SINGLETT SUBSYSTEM WITH | |
4566 | // FINAL STATE Q-QBAR PAIR FROM GRAPH | |
4567 | ||
4568 | // Prepare to check for colour singlet combinations of final state quarks | |
4569 | // Start by building a list of partons to exclude when checking for | |
4570 | // colour singlet combinations | |
4571 | int flav = event[emt].id(); | |
4572 | vector<int> exclude; | |
4573 | exclude.push_back(emt); | |
4574 | exclude.push_back(rad); | |
4575 | exclude.push_back(radBeforeColP[0]); | |
4576 | exclude.push_back(radBeforeColP[1]); | |
4577 | vector<int> colSinglet; | |
4578 | // Now find parton from which to start checking colour singlets | |
4579 | int iOther = -1; | |
4580 | // Loop through event to find a parton of correct flavour | |
4581 | for(int i=0; i < int(event.size()); ++i) | |
4582 | // Check final state for parton equalling emitted flavour. | |
4583 | // Exclude the colour system coupled to the clustering | |
4584 | if ( i != emt | |
4585 | && i != rad | |
4586 | && i != radBeforeColP[0] | |
4587 | && i != radBeforeColP[1] | |
4588 | && event[i].isFinal() ) { | |
4589 | // Stop if one parton of the correct flavour is found | |
4590 | if (event[i].id() == flav) { | |
4591 | iOther = i; | |
4592 | break; | |
4593 | } | |
4594 | } | |
4595 | // Save the type of flavour | |
4596 | int flavType = (event[iOther].id() > 0) ? 1 : -1; | |
4597 | // Try finding colour singlets | |
4598 | bool isColSing = getColSinglet(flavType,iOther,event,exclude,colSinglet); | |
4599 | // Check if colour singlet also is a flavour singlet | |
4600 | bool isFlavSing = isFlavSinglet(event,colSinglet); | |
4601 | ||
4602 | // Nearly there... | |
4603 | if (isColSing && isFlavSing) { | |
4604 | ||
4605 | // In a final check, ensure that the final state does not only | |
4606 | // consist of colour singlets that are also flavour singlets | |
4607 | // of the identical (!) flavours | |
4608 | // Loop through event and save all final state partons | |
4609 | vector<int> allFinal; | |
4610 | for(int i=0; i < int(event.size()); ++i) | |
4611 | if ( event[i].isFinal() ) | |
4612 | allFinal.push_back(i); | |
4613 | ||
4614 | // Check if all final partons form a colour singlet | |
4615 | bool isFullColSing = isColSinglet(event,allFinal); | |
4616 | // Check if all final partons form a flavour singlet | |
4617 | bool isFullFlavSing = isFlavSinglet(event,allFinal,flav); | |
4618 | ||
4619 | // If all final quarks are of identical flavour, | |
4620 | // no possible clustering should be discriminated. | |
4621 | // Otherwise, disallow | |
4622 | if (!isFullColSing || !isFullFlavSing) | |
4623 | return false; | |
4624 | } | |
4625 | } | |
4626 | ||
4627 | // ISR PROBLEMS | |
4628 | ||
4629 | if (isISRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) { | |
4630 | ||
4631 | // Easiest problem 1: | |
4632 | // RECLUSTERED INITIAL STATE GLUON MATCHES FINAL STATE GLUON | |
4633 | for(int i=0; i < int(gluon.size()); ++i) { | |
4634 | if (event[gluon[i]].isFinal() | |
4635 | && event[gluon[i]].col() == radBeforeCol | |
4636 | && event[gluon[i]].acol() == radBeforeAcl) | |
4637 | return false; | |
4638 | } | |
4639 | ||
4640 | // Easiest problem 2: | |
4641 | // RECLUSTERED INITIAL STATE GLUON MATCHES INITIAL STATE GLUON | |
4642 | for(int i=0; i < int(gluon.size()); ++i) { | |
4643 | if (event[gluon[i]].status() == -21 | |
4644 | && event[gluon[i]].acol() == radBeforeCol | |
4645 | && event[gluon[i]].col() == radBeforeAcl) | |
4646 | return false; | |
4647 | } | |
4648 | ||
4649 | // Next-to-easiest problem 1: | |
4650 | // RECLUSTERED INITIAL STATE GLUON MATCHES FINAL STATE Q-QBAR PAIR | |
4651 | if ( int(radBeforeColP.size()) == 2 | |
4652 | && event[radBeforeColP[0]].isFinal() | |
4653 | && event[radBeforeColP[1]].isFinal() | |
4654 | && event[radBeforeColP[0]].id() == -1*event[radBeforeColP[1]].id() ) { | |
4655 | ||
4656 | // In principle, clustering this splitting can disconnect | |
4657 | // the colour lines of a graph. However, the colours can be connected | |
4658 | // again if final state partons of the correct (anti)flavour, or | |
4659 | // initial state partons of the correct flavour exist | |
4660 | // Loop through event to check | |
4661 | bool clusPossible = false; | |
4662 | for(int i=0; i < int(event.size()); ++i) | |
4663 | if ( i != emt && i != rad | |
4664 | && i != radBeforeColP[0] | |
4665 | && i != radBeforeColP[1] | |
4666 | && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,event) ) { | |
4667 | if (event[i].status() == -21 | |
4668 | && ( event[radBeforeColP[0]].id() == event[i].id() | |
4669 | || event[radBeforeColP[1]].id() == event[i].id() )) | |
4670 | ||
4671 | clusPossible = true; | |
4672 | if (event[i].isFinal() | |
4673 | && ( event[radBeforeColP[0]].id() == -1*event[i].id() | |
4674 | || event[radBeforeColP[1]].id() == -1*event[i].id() )) | |
4675 | clusPossible = true; | |
4676 | } | |
4677 | ||
4678 | // There can be a further complication: If e.g. in | |
4679 | // t-channel photon exchange topologies, both incoming | |
4680 | // partons are quarks, and form colour singlets with any | |
4681 | // number of final state partons, at least try to | |
4682 | // recluster as much as possible. | |
4683 | // For this, check if the incoming parton | |
4684 | // connected to the radiator is connected to a | |
4685 | // colour and flavour singlet | |
4686 | int incoming1 = 3; | |
4687 | vector<int> excludeIn1; | |
4688 | for(int i=0; i < 4; ++i) | |
4689 | excludeIn1.push_back(0); | |
4690 | vector<int> colSingletIn1; | |
4691 | int flavIn1Type = (event[incoming1].id() > 0) ? 1 : -1; | |
4692 | // Try finding colour singlets | |
4693 | bool isColSingIn1 = getColSinglet(flavIn1Type,incoming1,event, | |
4694 | excludeIn1,colSingletIn1); | |
4695 | // Check if colour singlet also is a flavour singlet | |
4696 | bool isFlavSingIn1 = isFlavSinglet(event,colSingletIn1); | |
4697 | ||
4698 | // Check if the incoming parton not | |
4699 | // connected to the radiator is connected to a | |
4700 | // colour and flavour singlet | |
4701 | int incoming2 = 4; | |
4702 | vector<int> excludeIn2; | |
4703 | for(int i=0; i < 4; ++i) | |
4704 | excludeIn2.push_back(0); | |
4705 | vector<int> colSingletIn2; | |
4706 | int flavIn2Type = (event[incoming2].id() > 0) ? 1 : -1; | |
4707 | // Try finding colour singlets | |
4708 | bool isColSingIn2 = getColSinglet(flavIn2Type,incoming2,event, | |
4709 | excludeIn2,colSingletIn2); | |
4710 | // Check if colour singlet also is a flavour singlet | |
4711 | bool isFlavSingIn2 = isFlavSinglet(event,colSingletIn2); | |
4712 | ||
4713 | // If no "recovery clustering" is possible, reject clustering | |
4714 | if (!clusPossible | |
4715 | && (!isColSingIn1 || !isFlavSingIn1 | |
4716 | || !isColSingIn2 || !isFlavSingIn2)) | |
4717 | return false; | |
4718 | ||
4719 | } | |
4720 | ||
4721 | } | |
4722 | ||
4723 | // Done | |
4724 | return allowed; | |
4725 | } | |
4726 | ||
4727 | //-------------------------------------------------------------------------- | |
4728 | ||
4729 | // Function to check if rad,emt,rec triple is results in | |
4730 | // colour singlet radBefore+recBefore | |
4731 | // IN int rad,emt,rec : Positions (in event record) of the three | |
4732 | ||
4733 | // particles considered for clustering | |
4734 | // Event event : Reference event | |
4735 | ||
4736 | bool History::isSinglett( int rad, int emt, int rec, const Event& event ) { | |
4737 | ||
4738 | int radCol = event[rad].col(); | |
4739 | int emtCol = event[emt].col(); | |
4740 | int recCol = event[rec].col(); | |
4741 | int radAcl = event[rad].acol(); | |
4742 | int emtAcl = event[emt].acol(); | |
4743 | int recAcl = event[rec].acol(); | |
4744 | int recType = event[rec].isFinal() ? 1 : -1; | |
4745 | ||
4746 | bool isSing = false; | |
4747 | ||
4748 | if ( ( recType == -1 | |
4749 | && radCol + emtCol == recCol && radAcl + emtAcl == recAcl) | |
4750 | ||( recType == 1 | |
4751 | && radCol + emtCol == recAcl && radAcl + emtAcl == recCol) ) | |
4752 | isSing = true; | |
4753 | ||
4754 | return isSing; | |
4755 | ||
4756 | } | |
4757 | ||
4758 | //-------------------------------------------------------------------------- | |
4759 | ||
4760 | // Function to check if event is sensibly constructed: Meaning | |
4761 | // that all colour indices are contracted and that the charge in | |
4762 | // initial and final states matches | |
4763 | // IN event : event to be checked | |
4764 | // OUT TRUE : event is properly construced | |
4765 | // FALSE : event not valid | |
4766 | ||
4767 | bool History::validEvent( const Event& event ) { | |
4768 | ||
4769 | // Check if event is coloured | |
4770 | bool validColour = true; | |
4771 | for ( int i = 0; i < event.size(); ++i) | |
4772 | // Check colour of quarks | |
4773 | if ( event[i].isFinal() && event[i].colType() == 1 | |
4774 | // No corresponding anticolour in final state | |
4775 | && ( FindCol(event[i].col(),i,0,event,1,true) == 0 | |
4776 | // No corresponding colour in initial state | |
4777 | && FindCol(event[i].col(),i,0,event,2,true) == 0 )) { | |
4778 | validColour = false; | |
4779 | break; | |
4780 | // Check anticolour of antiquarks | |
4781 | } else if ( event[i].isFinal() && event[i].colType() == -1 | |
4782 | // No corresponding colour in final state | |
4783 | && ( FindCol(event[i].acol(),i,0,event,2,true) == 0 | |
4784 | // No corresponding anticolour in initial state | |
4785 | && FindCol(event[i].acol(),i,0,event,1,true) == 0 )) { | |
4786 | validColour = false; | |
4787 | break; | |
4788 | // No uncontracted colour (anticolour) charge of gluons | |
4789 | } else if ( event[i].isFinal() && event[i].colType() == 2 | |
4790 | // No corresponding anticolour in final state | |
4791 | && ( FindCol(event[i].col(),i,0,event,1,true) == 0 | |
4792 | // No corresponding colour in initial state | |
4793 | && FindCol(event[i].col(),i,0,event,2,true) == 0 ) | |
4794 | // No corresponding colour in final state | |
4795 | && ( FindCol(event[i].acol(),i,0,event,2,true) == 0 | |
4796 | // No corresponding anticolour in initial state | |
4797 | && FindCol(event[i].acol(),i,0,event,1,true) == 0 )) { | |
4798 | validColour = false; | |
4799 | break; | |
4800 | } | |
4801 | ||
4802 | // Check charge sum in initial and final state | |
4803 | bool validCharge = true; | |
4804 | double initCharge = event[3].charge() + event[4].charge(); | |
4805 | double finalCharge = 0.0; | |
4806 | for(int i = 0; i < event.size(); ++i) | |
4807 | if (event[i].isFinal()) finalCharge += event[i].charge(); | |
4808 | if (abs(initCharge-finalCharge) > 1e-12) validCharge = false; | |
4809 | ||
4810 | return (validColour && validCharge); | |
4811 | ||
4812 | } | |
4813 | ||
4814 | //-------------------------------------------------------------------------- | |
4815 | ||
4816 | // Function to check whether two clusterings are identical, used | |
4817 | // for finding the history path in the mother -> children direction | |
4818 | ||
4819 | bool History::equalClustering( Clustering clus1 , Clustering clus2 ) { | |
4820 | return ( (clus1.emittor == clus2.emittor) | |
4821 | && (clus1.emitted == clus2.emitted) | |
4822 | && (clus1.recoiler == clus2.recoiler) | |
4823 | && (clus1.partner == clus2.partner) | |
4824 | && (clus1.pT() == clus2.pT()) ); | |
4825 | } | |
4826 | ||
4827 | //-------------------------------------------------------------------------- | |
4828 | ||
4829 | // Chose dummy scale for event construction. By default, choose | |
4830 | // sHat for 2->Boson(->2)+ n partons processes and | |
4831 | // M_Boson for 2->Boson(->) processes | |
4832 | ||
4833 | double History::choseHardScale( const Event& event ) const { | |
4834 | ||
4835 | // Get sHat | |
4836 | double mHat = (event[3].p() + event[4].p()).mCalc(); | |
4837 | ||
4838 | // Find number of final state particles and bosons | |
4839 | int nFinal = 0; | |
4840 | int nFinBos= 0; | |
4841 | int nBosons= 0; | |
4842 | double mBos = 0.0; | |
4843 | for(int i = 0; i < event.size(); ++i) | |
4844 | if ( event[i].isFinal() ) { | |
4845 | nFinal++; | |
4846 | // Remember final state unstable bosons | |
4847 | if ( event[i].idAbs() == 23 | |
4848 | || event[i].idAbs() == 24 ) { | |
4849 | nFinBos++; | |
4850 | nBosons++; | |
4851 | mBos += event[i].m(); | |
4852 | } | |
4853 | } else if ( abs(event[i].status()) == 22 | |
4854 | && ( event[i].idAbs() == 23 | |
4855 | || event[i].idAbs() == 24 )) { | |
4856 | nBosons++; | |
4857 | mBos += event[i].m(); // Real mass | |
4858 | } | |
4859 | ||
4860 | // Return averaged boson masses | |
4861 | if ( nBosons > 0 && (nFinal + nFinBos*2) <= 3) | |
4862 | return (mBos / double(nBosons)); | |
4863 | else return | |
4864 | mHat; | |
4865 | } | |
4866 | ||
4867 | //-------------------------------------------------------------------------- | |
4868 | ||
4869 | // If the state has an incoming hadron return the flavour of the | |
4870 | // parton entering the hard interaction. Otherwise return 0 | |
4871 | ||
4872 | int History::getCurrentFlav(const int side) const { | |
4873 | int in = (side == 1) ? 3 : 4; | |
4874 | return state[in].id(); | |
4875 | } | |
4876 | ||
4877 | //-------------------------------------------------------------------------- | |
4878 | ||
4879 | double History::getCurrentX(const int side) const { | |
4880 | int in = (side == 1) ? 3 : 4; | |
4881 | return ( 2.*state[in].e()/state[0].e() ); | |
4882 | } | |
4883 | ||
4884 | //-------------------------------------------------------------------------- | |
4885 | ||
4886 | double History::getCurrentZ(const int rad, | |
4887 | const int rec, const int emt) const { | |
4888 | ||
4889 | int type = state[rad].isFinal() ? 1 : -1; | |
4890 | double z = 0.; | |
4891 | ||
4892 | if (type == 1) { | |
4893 | // Construct 2->3 variables for FSR | |
4894 | Vec4 sum = state[rad].p() + state[rec].p() | |
4895 | + state[emt].p(); | |
4896 | double m2Dip = sum.m2Calc(); | |
4897 | double x1 = 2. * (sum * state[rad].p()) / m2Dip; | |
4898 | double x3 = 2. * (sum * state[emt].p()) / m2Dip; | |
4899 | // Calculate z of splitting, different for FSR | |
4900 | z = x1/(x1+x3); | |
4901 | } else { | |
4902 | // Construct momenta of dipole before/after splitting for ISR | |
4903 | Vec4 qBR(state[rad].p() - state[emt].p() + state[rec].p()); | |
4904 | Vec4 qAR(state[rad].p() + state[rec].p()); | |
4905 | // Calculate z of splitting, different for ISR | |
4906 | z = (qBR.m2Calc())/( qAR.m2Calc()); | |
4907 | } | |
4908 | ||
4909 | return z; | |
4910 | ||
4911 | } | |
4912 | ||
4913 | //-------------------------------------------------------------------------- | |
4914 | ||
4915 | // Function to compute "pythia pT separation" from Particle input | |
4916 | ||
4917 | double History::pTLund(const Particle& RadAfterBranch, | |
4918 | const Particle& EmtAfterBranch, | |
4919 | const Particle& RecAfterBranch, int ShowerType) { | |
4920 | ||
4921 | // Save type: 1 = FSR pT definition, else ISR definition | |
4922 | int Type = ShowerType; | |
4923 | // Calculate virtuality of splitting | |
4924 | int sign = (Type==1) ? 1 : -1; | |
4925 | Vec4 Q(RadAfterBranch.p() + sign*EmtAfterBranch.p()); | |
4926 | double Qsq = sign * Q.m2Calc(); | |
4927 | // Mass term of radiator | |
4928 | double m2Rad = (mergingHooksPtr->includeMassive() | |
4929 | && RadAfterBranch.idAbs() >= 4 | |
4930 | && RadAfterBranch.idAbs() < 7) | |
4931 | ? pow(particleDataPtr->m0(RadAfterBranch.id()), 2) | |
4932 | : 0.; | |
4933 | // Construct 2->3 variables for FSR | |
4934 | Vec4 sum = RadAfterBranch.p() + RecAfterBranch.p() | |
4935 | + EmtAfterBranch.p(); | |
4936 | double m2Dip = sum.m2Calc(); | |
4937 | double x1 = 2. * (sum * RadAfterBranch.p()) / m2Dip; | |
4938 | double x3 = 2. * (sum * EmtAfterBranch.p()) / m2Dip; | |
4939 | // Construct momenta of dipole before/after splitting for ISR | |
4940 | Vec4 qBR(RadAfterBranch.p() - EmtAfterBranch.p() + RecAfterBranch.p()); | |
4941 | Vec4 qAR(RadAfterBranch.p() + RecAfterBranch.p()); | |
4942 | // Calculate z of splitting, different for FSR and ISR | |
4943 | double z = (Type==1) ? x1/(x1+x3) | |
4944 | : (qBR.m2Calc())/( qAR.m2Calc()); | |
4945 | // Separation of splitting, different for FSR and ISR | |
4946 | double pTpyth = (Type==1) ? z*(1.-z) : (1.-z); | |
4947 | // pT^2 = separation*virtuality | |
4948 | pTpyth *= (Qsq - sign*m2Rad); | |
4949 | if ( pTpyth < 0. ) pTpyth = 0.; | |
4950 | ||
4951 | // Return pT | |
4952 | return sqrt(pTpyth); | |
4953 | } | |
4954 | ||
4955 | //-------------------------------------------------------------------------- | |
4956 | ||
4957 | // Function to return the position of the initial line before (or after) | |
4958 | // a single (!) splitting. | |
4959 | ||
4960 | int History::posChangedIncoming(const Event& event, bool before) { | |
4961 | ||
4962 | // Check for initial state splittings. | |
4963 | // Consider a splitting to exist if both mother and sister were found. | |
4964 | // Find sister | |
4965 | int iSister = 0; | |
4966 | for (int i =0; i < event.size(); ++i) | |
4967 | if (event[i].status() == 43) { | |
4968 | iSister = i; | |
4969 | break; | |
4970 | } | |
4971 | // Find mother | |
4972 | int iMother = 0; | |
4973 | if (iSister > 0) iMother = event[iSister].mother1(); | |
4974 | ||
4975 | // Initial state splitting has been found. | |
4976 | if (iSister > 0 && iMother > 0) { | |
4977 | ||
4978 | // Find flavour, mother flavour | |
4979 | int flavSister = event[iSister].id(); | |
4980 | int flavMother = event[iMother].id(); | |
4981 | ||
4982 | // Find splitting flavour | |
4983 | int flavDaughter = 0; | |
4984 | if ( abs(flavMother) < 21 && flavSister == 21) | |
4985 | flavDaughter = flavMother; | |
4986 | else if ( flavMother == 21 && flavSister == 21) | |
4987 | flavDaughter = flavMother; | |
4988 | else if ( flavMother == 21 && abs(flavSister) < 21) | |
4989 | flavDaughter = -1*flavSister; | |
4990 | else if ( abs(flavMother) < 21 && abs(flavSister) < 21) | |
4991 | flavDaughter = 21; | |
4992 | ||
4993 | // Find initial state (!) daughter | |
4994 | int iDaughter = 0; | |
4995 | for (int i =0; i < event.size(); ++i) | |
4996 | if ( !event[i].isFinal() | |
4997 | && event[i].mother1() == iMother | |
4998 | && event[i].id() == flavDaughter ) | |
4999 | iDaughter = i; | |
5000 | ||
5001 | // Done for initial state splitting. | |
5002 | if ( !before ) return iMother; | |
5003 | else return iDaughter; | |
5004 | ||
5005 | } | |
5006 | ||
5007 | // Check for final state splittings with initial state recoiler. | |
5008 | // Consider a splitting to exist if both mother and daughter were found. | |
5009 | // Find new mother | |
5010 | iMother = 0; | |
5011 | for (int i =0; i < event.size(); ++i) | |
5012 | if ( abs(event[i].status()) == 53 || abs(event[i].status()) == 54) { | |
5013 | iMother = i; | |
5014 | break; | |
5015 | ||
5016 | } | |
5017 | // Find daughter | |
5018 | int iDaughter = 0; | |
5019 | if (iMother > 0) iDaughter = event[iMother].daughter1(); | |
5020 | ||
5021 | // Done if final state splitting has been found. | |
5022 | if (iDaughter > 0 && iMother > 0) { | |
5023 | ||
5024 | // Done for final state splitting. | |
5025 | if ( !before ) return iMother; | |
5026 | else return iDaughter; | |
5027 | ||
5028 | } | |
5029 | ||
5030 | // If no splitting has been found, return zero. | |
5031 | return 0; | |
5032 | ||
5033 | } | |
5034 | ||
5035 | //========================================================================== | |
5036 | ||
5037 | } // end namespace Pythia8 |