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 |
---|