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.
6 // This file is written by Stefan Prestel.
7 // Function definitions (not found in the header) for the
8 // Clustering and History classes.
14 //==========================================================================
16 // The Clustering class.
18 //--------------------------------------------------------------------------
20 // Declaration of Clustering class
21 // This class holds information about one radiator, recoiler,
23 // This class is a container class for History class use.
26 void Clustering::list() const {
27 cout << " emt " << emitted
29 << " rec " << recoiler
30 << " partner " << partner
31 << " pTscale " << pTscale << endl;
34 //==========================================================================
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
49 //--------------------------------------------------------------------------
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).
64 History::History( int depth,
68 MergingHooks* mergingHooksPtrIn,
71 ParticleData* particleDataPtrIn,
73 bool isOrdered = true,
74 bool isUnordered = true,
75 bool isStronglyOrdered = true,
76 bool isAllowed = true,
84 foundOrderedPath(false),
85 foundUnorderedPath(false),
86 foundStronglyOrderedPath(false),
87 foundAllowedPath(false),
88 foundCompletePath(false),
94 mergingHooksPtr(mergingHooksPtrIn),
97 particleDataPtr(particleDataPtrIn),
101 // Initialise beam particles
103 // Update probability with PDF ratio
104 if (mother && mergingHooksPtr->includeRedundant()) prob *= pdfForSudakov();
106 // Minimal scalar sum of pT used in Herwig to choose history
107 // Keep track of scalar PT
109 double acoll = (mother->state[clusterIn.emittor].isFinal())
110 ? mergingHooksPtr->herwigAcollFSR()
111 : mergingHooksPtr->herwigAcollISR();
112 sumScalarPT = mother->sumScalarPT + acoll*scale;
116 // Remember reclustered radiator in lower multiplicity state
117 if ( mother ) iReclusteredOld = mother->iReclusteredNew;
119 // Check if more steps should be taken.
122 for ( int i = 0; i < int(state.size()); ++i )
123 if ( state[i].isFinal() ) {
124 if ( state[i].colType() != 0 )
126 if ( state[i].idAbs() == 24 )
129 if ( mergingHooksPtr->doWClustering()
130 && nFinalP == 2 && nFinalW == 0 ) depth = 0;
132 // If this is not the fully clustered state, try to find possible
134 vector<Clustering> clusterings;
135 if ( depth > 0 ) clusterings = getAllQCDClusterings();
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() );
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 );
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]));
161 for ( multimap<double, Clustering *>::iterator it = sorted.begin();
162 it != sorted.end(); ++it ) {
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;
175 bool ordered = isOrdered;
176 bool unordered = isUnordered;
177 if ( mergingHooksPtr->orderInRapidity()
178 && mergingHooksPtr->orderHistories() ) {
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;
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;
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;
209 // Check if reclustered state should be disallowed
210 bool doCut = mergingHooksPtr->canCutOnRecState()
211 || mergingHooksPtr->allowCutOnRecState();
212 bool allowed = isAllowed;
215 && mergingHooksPtr->doCutOnRecState(cluster(*it->second)) ) {
216 if ( onlyAllowedPaths() ) continue;
220 // Perform the clustering and recurse and construct the next
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 ));
229 //--------------------------------------------------------------------------
231 // Function to project all possible paths onto only the desired paths.
233 bool History::projectOntoDesiredHistories() {
234 // At the moment, only trim histories.
235 return trimHistories();
238 //--------------------------------------------------------------------------
240 // In the initial history node, select one of the paths according to
241 // the probabilities. This function should be called for the initial
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
248 // OUT double : (Sukadov) , (alpha_S ratios) , (PDF ratios)
250 double History::weightTREE(PartonLevel* trial, AlphaStrong * asFSR,
251 AlphaStrong * asISR, double RN) {
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();
261 double asWeight = 1.;
262 double pdfWeight = 1.;
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);
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
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);
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
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;
294 return (wt*asWeight*pdfWeight);
297 //--------------------------------------------------------------------------
299 // Function to set the state with complete scales for evolution.
301 void History::getStartingConditions( const double RN, Event& outState ) {
303 // Select the history
304 History * selected = select(RN);
307 // Copy the output state
310 // Set the scale of the lowest order process
311 if (!selected->mother) {
313 for(int i=0; i < int(outState.size()); ++i)
314 if (outState[i].isFinal()) nFinal++;
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
326 // outState.scale(0.5*outState[5].m());
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
338 infoPtr->zNowISR(selected->zISR());
339 infoPtr->pT2NowISR(pow(selected->pTISR(),2));
340 infoPtr->hasHistory(true);
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);
354 //--------------------------------------------------------------------------
356 // Function to set the state with complete scales for evolution.
358 bool History::getClusteredEvent( const double RN, Event& outState,
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);
378 //--------------------------------------------------------------------------
380 // Calculate and return pdf ratio.
382 double History::getPDFratio( int side, bool forSudakov,
383 int flavNum, double xNum, double muNum,
384 int flavDen, double xDen, double muDen) {
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;
390 // Now calculate PDF ratio if necessary
391 double pdfRatio = 1.0;
393 // Get mother and daughter pdfs
397 // Use rescaled PDFs in the presence of multiparton interactions
400 pdfNum = mother->beamA.xfISR(0, flavNum, xNum, muNum*muNum);
402 pdfNum = beamA.xfISR(0, flavNum, xNum, muNum*muNum);
404 pdfDen = max(1e-10, beamA.xfISR(0, flavDen, xDen, muDen*muDen));
406 pdfDen = max(1e-10, beamA.xfISR(0, flavDen, xDen, muDen*muDen));
410 pdfNum = mother->beamB.xfISR(0, flavNum, xNum, muNum*muNum);
412 pdfNum = beamB.xfISR(0, flavNum, xNum, muNum*muNum);
415 pdfDen = max(1e-10,beamB.xfISR(0, flavDen, xDen, muDen*muDen));
417 pdfDen = max(1e-10,beamB.xfISR(0, flavDen, xDen, muDen*muDen));
420 // Return ratio of pdfs
421 if ( pdfNum > 1e-15 && pdfDen > 1e-10 ) {
422 pdfRatio *= pdfNum / pdfDen;
424 } else if ( pdfNum < pdfDen ) {
426 } else if ( pdfNum > pdfDen ) {
435 //--------------------------------------------------------------------------
437 /*--------------- METHODS USED FOR ONLY ONE PATH OF HISTORY NODES ------- */
439 // Function to set all scales in the sequence of states. This is a
440 // wrapper routine for setScales and setEventScales methods
442 void History::setScalesInHistory() {
443 // Find correct links from n+1 to n states (mother --> child), as
444 // needed for enforcing ordered scale sequences
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
454 //--------------------------------------------------------------------------
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
464 void History::findPath(vector<int>& out) {
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
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)) {
482 // Save the index of the child in the children vector and recurse
484 out.push_back(iChild);
485 mother->findPath(out);
489 //--------------------------------------------------------------------------
491 // Functions to set the parton production scales and enforce
492 // ordering on the scales of the respective clusterings stored in
494 // Method will start from lowest multiplicity state and move to
495 // higher states, setting the production scales the shower would
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
510 void History::setScales( vector<int> index, bool forward) {
512 // First, set the scales of the hard process to the kinematial
514 if ( children.empty() && forward ) {
515 // New "incomplete" configurations showered from mu
517 double scaleNew = 1.;
519 if (mergingHooksPtr->incompleteScalePrescip()==0) {
520 scaleNew = infoPtr->QFac();
521 } else if (mergingHooksPtr->incompleteScalePrescip()==1) {
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();
532 scaleNew = max( mergingHooksPtr->pTcut(), scaleNew);
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);
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() );
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);
549 if (state[i].isFinal() && state[i].colType() == 0) {
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());
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));
573 scaleNew *= max( mergingHooksPtr->pTcut(), max(scale,mother->scale));
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);
582 // Find unchanged copies of partons in higher multiplicity states
584 mother->scaleCopies(clusterIn.emitted, mother->state, scaleNew);
585 mother->scaleCopies(clusterIn.emittor, mother->state, scaleNew);
586 mother->scaleCopies(clusterIn.recoiler, mother->state, scaleNew);
589 mother->setScales(index,true);
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
597 if ( int(index.size()) > 0 ) {
598 iChild = index.back();
602 // Check that the reclustered scale is above the shower cut
604 scale = max(mergingHooksPtr->pTcut(), scale);
606 // If this is NOT the 2->2 process, check and enforce ordering
607 if (iChild != -1 && !children.empty()) {
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;
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
632 // Just set the overall event scale to the minimal scale
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);
643 children[iChild]->setScales(index, false);
649 //--------------------------------------------------------------------------
651 // Function to find a particle in all higher multiplicity events
652 // along the history path and set its production scale to the input
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
659 void History::scaleCopies(int iPart, const Event& refEvent, double rho) {
661 // Check if any parton recently rescaled is found unchanged:
662 // Same charge, colours in mother->state
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() )
671 // Rescale the unchanged parton
672 mother->state[i].scale(rho);
675 mother->scaleCopies( iPart, refEvent, rho );
676 } // end if found unchanged parton case
677 } // end loop over particle entries in event
681 //--------------------------------------------------------------------------
683 // Functions to set the OVERALL EVENT SCALES [=state.scale()] to
684 // the scale of the last clustering
688 void History::setEventScales() {
689 // Set the event scale to the scale of the last clustering,
690 // except for the very lowest multiplicity state
692 mother->state.scale(scale);
694 mother->setEventScales();
698 //--------------------------------------------------------------------------
700 // Functions to return the z value of the last ISR splitting
702 // OUTPUT double : z value of last ISR splitting in history
704 double History::zISR() {
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();
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();
718 double znew = mother->zISR();
720 if (znew > 0.) z = znew;
725 //--------------------------------------------------------------------------
727 // Functions to return the z value of the last FSR splitting
729 // OUTPUT double : z value of last FSR splitting in history
731 double History::zFSR() {
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();
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);
750 double znew = mother->zFSR();
752 if (znew > 0.) z = znew;
757 //--------------------------------------------------------------------------
759 // Functions to return the pT scale of the last FSR splitting
761 // OUTPUT double : pT scale of last FSR splitting in history
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();
770 double pTnew = mother->pTISR();
772 if (pTnew > 0.) pT = pTnew;
777 //--------------------------------------------------------------------------
779 // Functions to return the pT scale of the last FSR splitting
781 // OUTPUT double : pT scale of last FSR splitting in history
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();
790 double pTnew = mother->pTFSR();
792 if (pTnew > 0.) pT = pTnew;
796 //--------------------------------------------------------------------------
798 // Functions to return the depth of the history (i.e. the number of
799 // reclustered splittings)
801 // OUTPUT int : Depth of history
803 int History::nClusterings() {
805 if (!mother) return 0;
806 int w = mother->nClusterings();
811 //--------------------------------------------------------------------------
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
819 Event History::clusteredState(int nSteps) {
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);
831 //--------------------------------------------------------------------------
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
838 History * History::select(double rnd) {
840 // No need to choose if no paths have been constructed.
841 if ( goodBranches.empty() && badBranches.empty() ) return this;
843 // Choose amongst paths allowed by projections.
845 map<double, History*> selectFrom;
846 if ( !goodBranches.empty() ) {
847 selectFrom = goodBranches;
848 sum = sumGoodBranches;
850 selectFrom = badBranches;
851 sum = sumBadBranches;
854 if (mergingHooksPtr->pickBySumPT()) {
855 // Find index of history with minimal sum of scalar pT
857 for (int i=0; i < state.size(); ++i)
858 if (state[i].isFinal())
861 double sumMin = (nFinal-2)*state[0].e();
862 for ( map<double, History*>::iterator it = selectFrom.begin();
863 it != selectFrom.end(); ++it ) {
865 if (it->second->sumScalarPT < sumMin) {
866 sumMin = it->second->sumScalarPT;
870 // Choose history with smallest sum of scalar pT
871 return selectFrom.lower_bound(iMin)->second;
873 // Choose history according to probability, be careful about upper bound
875 return selectFrom.upper_bound(sum*rnd)->second;
877 return selectFrom.lower_bound(sum*rnd)->second;
883 //--------------------------------------------------------------------------
885 // Function to project paths onto desired paths.
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();
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 ) {
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;
912 // Update mismatch in probabilities resulting from not including this
914 double mismatchOld = mismatch;
915 mismatch += sumnew - sumold;
916 // Fill branches with allowed paths.
917 badBranches.insert( make_pair( mismatchOld + sumnew - sumold,
919 // Add probability of this path.
920 sumBadBranches = mismatchOld + sumnew - sumold;
922 // remember index of this path in order to caclulate probability of
927 return !goodBranches.empty();
930 //--------------------------------------------------------------------------
932 // Function implementing checks on a paths, for deciding if the path should
933 // be considered valid or not.
935 bool History::keepHistory() {
936 bool keepPath = true;
937 //// Tag unordered paths for removal.
938 //double maxScale = (foundCompletePath)
940 // : infoPtr->QFac();
941 //keepPath = isOrderedPath( maxScale );
946 //--------------------------------------------------------------------------
948 // Function to check if a path is ordered in evolution pT.
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;
958 //--------------------------------------------------------------------------
960 // Function to check if all reconstucted states in a path pass the merging
963 bool History::allIntermediateAboveRhoMS( double rhoms ) {
966 for ( int i = 0; i < state.size(); ++i )
967 if ( state[i].isFinal() && state[i].colType() != 0 )
970 double rhoNew = mergingHooksPtr->rhoms( state, false );
972 if ( !mother ) return true;
974 bool passRHOMS = mother->allIntermediateAboveRhoMS( rhoms );
976 if ( !passRHOMS || ( nFinal > 0 && rhoNew < rhoms ) ) return false;
982 //--------------------------------------------------------------------------
984 // Function to check if reconstucted states in a path are ordered in the
985 // merging scale variable.
987 bool History::intermediateRhoMSOrdered( double maxscale ) {
988 // Count number of final state partons.
990 for ( int i = 0; i < state.size(); ++i )
991 if ( state[i].isFinal() && state[i].colType() != 0 )
993 // Get new min(rho) scale.
994 double newscale = (nFinal == 0 )
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 )
1002 return isOrdered && ( maxscale >= newscale );
1005 //--------------------------------------------------------------------------
1007 // Function to check if any ordered paths were found (and kept).
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) )
1022 //--------------------------------------------------------------------------
1024 // Function to check if any unordered paths were found (and kept).
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) )
1039 //--------------------------------------------------------------------------
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
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)
1054 double History::weightTree(PartonLevel* trial, double as0, double maxscale,
1055 double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,
1056 double& asWeight, double& pdfWeight) {
1058 // Use correct scale
1059 double newScale = scale;
1061 // For ME state, just multiply by PDF ratios
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();
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,
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,
1098 // Remember new PDF scale n case true sclae should be used for un-ordered
1100 double newPDFscale = newScale;
1101 if (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
1102 newPDFscale = clusterIn.pT();
1105 double w = mother->weightTree(trial, as0, newScale, newPDFscale,
1106 asFSR, asISR, asWeight, pdfWeight);
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;
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;
1128 // Calculate pdf ratios: Get both sides of event
1131 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
1132 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
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())
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,
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())
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,
1172 //--------------------------------------------------------------------------
1174 // Function to return the factorisation scale of the hard process in Pythia.
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
1184 if ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
1185 || mergingHooksPtr->getProcessString().compare("pp>aj") == 0 ) {
1186 // Find the mT in the hard sub-process.
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();
1194 hardscale = sqrt( min( mT[0], mT[1] ) );
1196 hardscale = infoPtr->QFac();
1202 //--------------------------------------------------------------------------
1204 // Function to return the factorisation scale of the hard process in Pythia.
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
1212 if ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
1213 || mergingHooksPtr->getProcessString().compare("pp>aj") == 0 ) {
1214 // Find the mT in the hard sub-process.
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();
1223 hardscale = sqrt( mT[0]*mT[1] );
1225 hardscale = infoPtr->QRen();
1231 //--------------------------------------------------------------------------
1233 // Perform a trial shower using the \a pythia object between
1234 // maxscale down to this scale and return the corresponding Sudakov
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 )
1242 double History::doTrialShower( PartonLevel* trial, double maxscale,
1245 // Copy state to local process
1246 Event process = state;
1247 // Set starting scale.
1248 double startingScale = maxscale;
1250 bool doVeto = false;
1254 // Reset trialShower object
1255 trial->resetTrial();
1256 // Construct event to be showered
1257 Event event = Event();
1258 event.init("(hard process-modified)", particleDataPtr);
1260 // Reset process scale so thatshower starting scale is correctly set.
1261 process.scale(startingScale);
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;
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
1279 : mother->getCurrentZ(clusterIn.emittor,clusterIn.recoiler,
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);
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();
1292 // Get veto (merging) scale value
1293 double vetoScale = (mother) ? 0. : mergingHooksPtr->tms();
1294 // Get merging scale in current event
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.
1307 tnow = mergingHooksPtr->tmsDefinition(event);
1309 // Done if evolution scale has fallen below minimum
1310 if ( pTtrial < minScale ) break;
1311 // Reset starting scale.
1312 startingScale = pTtrial;
1314 // Continue if this state is below the veto scale
1315 if ( tnow < vetoScale && vetoScale > 0. ) continue;
1317 // Retry if the trial emission was not allowed.
1318 if ( mergingHooksPtr->canVetoTrialEmission()
1319 && mergingHooksPtr->doVetoTrialEmission( process, event) ) continue;
1321 // Veto event if trial pT was above the next nodal scale.
1322 if ( pTtrial > minScale ) doVeto = true;
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;
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
1333 vector<int> finalPartons;
1334 for(int i=0; i < state.size(); ++i) {
1335 if (state[i].isFinal()) {
1337 if ( state[i].colType() != 0)
1338 finalPartons.push_back(i);
1341 // Veto if MPI was above 2 -> 2 pT
1342 if ( nFinal == 2 && int(finalPartons.size()) == 2
1343 && pTtrial > event[finalPartons[0]].pT() ) {
1348 // If pT of trial emission was in suitable range (trial shower
1349 // successful), return false
1350 if ( pTtrial < minScale ) doVeto = false;
1358 return ( (doVeto) ? 0. : 1. );
1361 /*--------------- METHODS USED FOR CONTRUCTION OF ALL HISTORIES --------- */
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.
1367 bool History::onlyOrderedPaths() {
1368 if ( !mother || foundOrderedPath ) return foundOrderedPath;
1369 return foundOrderedPath = mother->onlyOrderedPaths();
1372 bool History::onlyUnorderedPaths() {
1373 if ( !mother || foundUnorderedPath ) return foundUnorderedPath;
1374 return foundUnorderedPath = mother->onlyUnorderedPaths();
1377 //--------------------------------------------------------------------------
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.
1383 bool History::onlyStronglyOrderedPaths() {
1384 if ( !mother || foundStronglyOrderedPath ) return foundStronglyOrderedPath;
1385 return foundStronglyOrderedPath = mother->onlyStronglyOrderedPaths();
1388 //--------------------------------------------------------------------------
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.
1394 bool History::onlyAllowedPaths() {
1395 if ( !mother || foundAllowedPath ) return foundAllowedPath;
1396 return foundAllowedPath = mother->onlyAllowedPaths();
1399 //--------------------------------------------------------------------------
1401 // When a full path has been found, register it with the initial
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 ...)
1408 bool History::registerPath(History & l, bool isOrdered, bool isUnordered,
1409 bool isStronglyOrdered, bool isAllowed, bool isComplete) {
1411 // We are not interested in improbable paths.
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 )
1420 if ( mergingHooksPtr->canCutOnRecState()
1421 && foundAllowedPath && !isAllowed )
1423 if ( mergingHooksPtr->enforceStrongOrdering()
1424 && foundStronglyOrderedPath && !isStronglyOrdered )
1426 if ( mergingHooksPtr->orderHistories()
1427 && foundOrderedPath && !isOrdered )
1430 if ( mergingHooksPtr->forceUnorderedHistories()
1431 && foundUnorderedPath && !isUnordered )
1434 if ( foundCompletePath && !isComplete )
1436 if ( !mergingHooksPtr->canCutOnRecState()
1437 && !mergingHooksPtr->allowCutOnRecState() )
1438 foundAllowedPath = true;
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.
1447 foundAllowedPath = true;
1451 if ( mergingHooksPtr->enforceStrongOrdering() && isStronglyOrdered
1453 if ( !foundStronglyOrderedPath || !foundCompletePath ) {
1454 // If this is the first complete, ordered path, discard the
1455 // old, non-ordered or incomplete ones.
1459 foundStronglyOrderedPath = true;
1460 foundCompletePath = true;
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.
1471 foundOrderedPath = true;
1472 foundCompletePath = true;
1476 if ( mergingHooksPtr->forceUnorderedHistories() && isUnordered
1478 if ( !foundUnorderedPath || !foundCompletePath ) {
1479 // If this is the first complete, unordered path, discard the
1480 // old, ordered or incomplete ones.
1484 foundUnorderedPath = true;
1485 foundCompletePath = true;
1489 if ( !foundCompletePath ) {
1490 // If this is the first complete path, discard the old,
1495 foundCompletePath = true;
1498 // Remember, if this path is ordered, even if no ordering is required
1500 foundOrderedPath = true;
1503 // Index path by probability
1505 paths[sumpath] = &l;
1510 //--------------------------------------------------------------------------
1512 // For the history-defining state (and if necessary interfering
1513 // states), find all possible clusterings.
1515 // OUT vector of all (rad,rec,emt) systems
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;
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);
1549 // Get all clusterings for input state
1550 vector<Clustering> systems;
1551 systems = getQCDClusterings(state);
1552 ret.insert(ret.end(), systems.begin(), systems.end());
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
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],
1570 int col = NewState[PosFinalQuark[i]].col();
1571 for(int j = 0; j < int(PosInitAntiq.size()); ++j) {
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()) {
1581 ret.insert(ret.end(), systems.begin(), systems.end());
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],
1593 int acl = NewState[PosFinalAntiq[i]].acol();
1594 for(int j = 0; j < int(PosInitQuark.size()); ++j) {
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()) {
1604 ret.insert(ret.end(), systems.begin(), systems.end());
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);
1617 // If no colour rearrangements should be done, print warning and return
1619 string message="Warning in History::getAllQCDClusterings: No clusterings";
1620 message+=" found. History incomplete";
1621 infoPtr->errorMsg(message);
1627 //--------------------------------------------------------------------------
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
1633 vector<Clustering> History::getQCDClusterings( const Event& event) {
1634 vector<Clustering> ret;
1636 // Initialise vectors to keep track of position of partons in the
1638 vector <int> PosFinalPartn;
1639 vector <int> PosInitPartn;
1641 vector <int> PosFinalGluon;
1642 vector <int> PosFinalQuark;
1643 vector <int> PosFinalAntiq;
1644 vector <int> PosInitGluon;
1645 vector <int> PosInitQuark;
1646 vector <int> PosInitAntiq;
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);
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());
1676 vector<Clustering> systems;
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 )
1684 int EmtGluon = PosFinalGluon[i];
1685 systems = findQCDTriple( EmtGluon, 2, event, PosFinalPartn, PosInitPartn);
1686 ret.insert(ret.end(), systems.begin(), systems.end());
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 )
1695 && (nFiQuark == 1) && (nFiAntiq == 1) )
1696 || ( ( nFiQuark + nFiAntiq == 0)
1697 && (nInQuark == 1) && (nInAntiq == 1) ) )
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 )
1708 int EmtQuark = PosFinalQuark[i];
1709 systems = findQCDTriple( EmtQuark,1,event, PosFinalPartn, PosInitPartn);
1710 ret.insert(ret.end(), systems.begin(), systems.end());
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 )
1720 int EmtAntiq = PosFinalAntiq[i];
1721 systems = findQCDTriple( EmtAntiq,1,event, PosFinalPartn, PosInitPartn);
1722 ret.insert(ret.end(), systems.begin(), systems.end());
1730 //--------------------------------------------------------------------------
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
1742 vector<Clustering> History::findQCDTriple (int EmtTagIn, int colTopIn,
1744 vector<int> PosFinalPartn,
1745 vector <int> PosInitPartn ) {
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;
1753 // Initialise FinalSize
1754 int FinalSize = int(PosFinalPartn.size());
1755 int InitSize = int(PosInitPartn.size());
1756 int Size = InitSize + FinalSize;
1758 vector<Clustering> clus;
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];
1766 if ( event[iRad].col() == event[EmtTag].col()
1767 && event[iRad].acol() == event[EmtTag].acol() )
1770 //if ( mergingHooksPtr->nRecluster() > 0
1771 // && iRad == iReclusteredOld ) continue;
1773 if (iRad != EmtTag ) {
1774 int pTdef = event[iRad].isFinal() ? 1 : -1;
1775 int sign = (a < FinalSize)? 1 : -1 ;
1777 // First colour topology: g --> qqbar. Here, emt & rad should
1778 // have same flavour (causes problems for gamma->qqbar).
1781 if ( event[iRad].id() == -sign*event[EmtTag].id() ) {
1784 if (event[iRad].id() < 0) {
1785 col = event[EmtTag].acol();
1786 acl = event[iRad].acol();
1788 col = event[EmtTag].col();
1789 acl = event[iRad].col();
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
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
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) ));
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
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
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) ));
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
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
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) ));
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
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
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) ));
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)
1903 for(int l = 0; l < int(PosInitPartn.size()); ++l)
1904 if (PosInitPartn[l] != iRad) RecInit = PosInitPartn[l];
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);
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;
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);
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) ));
1937 // Second colour topology: Gluon emission
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"
1950 if (event[iRad].isFinal() ) {
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();
1959 col = event[iRad].col();
1960 acl = event[iRad].acol();
1965 iRec = FindCol(col,iRad,EmtTag,event,1,true);
1966 if ( (sign < 0) && (event[iRec].isFinal()) ) 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) ));
1974 iRec = FindCol(col,iRad,EmtTag,event,2,true);
1975 if ( (sign < 0) && (event[iRec].isFinal()) ) 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) ));
1985 iRec = FindCol(acl,iRad,EmtTag,event,1,true);
1986 if ( (sign < 0) && (event[iRec].isFinal()) ) 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) ));
1994 iRec = FindCol(acl,iRad,EmtTag,event,2,true);
1995 if ( (sign < 0) && (event[iRec].isFinal()) ) 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) ));
2006 // For an initial state radiator, always set recoiler
2007 // to the other initial state parton (recoil is taken
2009 // by full remaining system, so this is just a
2010 // labelling for such a process)
2013 for(int l = 0; l < int(PosInitPartn.size()); ++l)
2014 if (PosInitPartn[l] != iRad) RecInit = PosInitPartn[l];
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);
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);
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)));
2053 //--------------------------------------------------------------------------
2055 // For the history-defining state (and if necessary interfering
2056 // states), find all possible clusterings.
2058 // OUT vector of all (rad,rec,emt) systems
2060 vector<Clustering> History::getAllEWClusterings() {
2061 vector<Clustering> ret;
2063 // Get all clusterings for input state
2064 vector<Clustering> systems;
2065 systems = getEWClusterings(state);
2066 ret.insert(ret.end(), systems.begin(), systems.end());
2071 //--------------------------------------------------------------------------
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
2077 vector<Clustering> History::getEWClusterings( const Event& event) {
2078 vector<Clustering> ret;
2080 // Initialise vectors to keep track of position of partons in the
2082 vector <int> PosFinalPartn;
2083 vector <int> PosInitPartn;
2084 vector <int> PosFinalW;
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);
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 );
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());
2114 //--------------------------------------------------------------------------
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
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)
2133 // Initialise FinalSize
2134 int FinalSize = int(PosFinalPartn.size());
2136 vector<Clustering> clus;
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 ) {
2144 // Find recoiler by flavour.
2145 int flavRad = event[iRad].id();
2146 int flavEmt = event[EmtTag].id();
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) ) );
2163 //--------------------------------------------------------------------------
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
2170 double History::getProb(const Clustering & SystemIn) {
2172 // Get local copies of input system
2173 int Rad = SystemIn.emittor;
2174 int Rec = SystemIn.recoiler;
2175 int Emt = SystemIn.emitted;
2177 // Initialise shower probability
2178 double showerProb = 0.0;
2180 // If the splitting resulted in disallowed evolution variable,
2181 // disallow the splitting
2182 if (SystemIn.pT() <= 0.) return 0.;
2184 // Initialise all combinatorical factors
2187 // Flavour is known when reclustring, thus n_f=1
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();
2195 // Check if this is the clustering 2->3 to 2->2.
2196 // If so, use weight for joined evolution
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));
2204 // Find incoming particles
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;
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);
2238 // Correction of virtuality for massive splittings
2239 if ( g2QQmassive) Q1sq += m2Emt0;
2240 else if (Q2Qgmassive) Q1sq += m2Rad0;
2242 // pT0 dependence!!!
2243 double pT0sq = pow(mergingHooksPtr->pT0ISR(),2);
2244 double Q2sq = -Q2.m2Calc();
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;
2257 bool hasJoinedEvol = (state[Emt].id() == 21
2258 || state[Rad].id() == state[Rec].id());
2260 // Initialise normalization factor multiplying the splitting
2261 // function numerator
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);
2268 fac = (hasJoinedEvol && isLast) ? facJoined : facSingle;
2270 } else if (mergingHooksPtr->pickByPoPT2()) {
2271 fac = 1./(pT1sq + pT0sq);
2274 // Calculate shower splitting probability:
2275 // Splitting functions*normalization*ME reweighting factors
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);
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
2288 for(int i=0; i < state.size(); ++i)
2289 if (state[i].isFinal() && state[i].colType() != 0
2290 && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) )
2292 // For first splitting of single vector boson production,
2293 // apply ME corrections
2295 && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1) {
2296 double sH = m2Dip / z1;
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;
2306 showerProb = num*fac*meReweighting;
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));
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
2319 for(int i=0; i < state.size(); ++i)
2320 if (state[i].isFinal() && state[i].colType() != 0
2321 && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) )
2323 // For first splitting of single vector boson production,
2324 // apply ME corrections
2326 && mergingHooksPtr->getProcessString().compare("pp>h") == 0
2327 && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1 ) {
2328 double sH = m2Dip / z1;
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));
2337 showerProb = num*fac*meReweighting;
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;
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
2350 for ( int i=0; i < state.size(); ++i )
2351 if ( state[i].isFinal() && state[i].colType() != 0
2352 && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) )
2354 // For first splitting of single vector boson production,
2355 // apply ME corrections
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));
2366 showerProb = num*fac*meReweighting;
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
2378 for ( int i=0; i < state.size(); ++i )
2379 if ( state[i].isFinal() && state[i].colType() != 0
2380 && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) )
2382 // For first splitting of single vector boson production,
2383 // apply ME corrections
2385 && int(mergingHooksPtr->hardProcess.hardIntermediate.size()) == 1) {
2386 double sH = m2Dip / z1;
2388 double uH = Q1sq - m2Dip * (1. - z1) / z1;
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;
2398 showerProb = num*fac*meReweighting;
2400 // Print error if no kernel calculated
2402 string message="Error in History::getProb: Splitting kernel undefined";
2403 message+=" in ISR clustering.";
2404 infoPtr->errorMsg(message);
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;
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)
2419 if (pT2QQcorr < 0.0) showerProb *= 1e-9;
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
2430 // Done for ISR case, begin FSR case
2432 } else if (isFSR || isFSRinREC) {
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);
2446 // Virtuality of the splittings
2447 Vec4 Q1( state[Rad].p() + state[Emt].p() );
2448 Vec4 Q2( state[Rec].p() + state[Emt].p() );
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();
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)
2462 if ( state[Rec].idAbs() >= 4 && state[Rec].idAbs() < 7)
2465 // Initialise normalization factor multiplying the splitting
2466 // function numerator
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()) {
2475 // Calculate shower splitting probability:
2476 // Splitting functions*normalization*ME reweighting factors
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);
2483 showerProb = num*fac;
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
2494 for(int i=0; i < state.size(); ++i)
2495 if (state[i].isFinal() && state[i].colType() != 0
2496 && !mergingHooksPtr->hardProcess.matchesAnyOutgoing(i,state) )
2498 // Calculate splitting kernel
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
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;
2515 showerProb = num*fac*meReweighting;
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);
2527 showerProb = num*fac;
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;
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) );
2552 showerProb = num*fac*beta;
2554 // Print error if no kernel calculated
2556 string message="Error in History::getProb: Splitting kernel undefined";
2557 message+=" in FSR clustering.";
2558 infoPtr->errorMsg(message);
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
2571 string message="Error in History::getProb: Radiation could not be";
2572 message+=" interpreted as FSR or ISR.";
2573 infoPtr->errorMsg(message);
2576 if (showerProb <= 0.) showerProb = 0.;
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;
2584 // AlphaEM* aEMFSR = mergingHooksPtr->AlphaEM_FSR();
2585 // double aEM = (*aEMFSR).alphaEM(91.188*91.188) / (2.*M_PI);
2586 // showerProb *= aEM;
2594 //--------------------------------------------------------------------------
2596 // Set up the beams (fill the beam particles with the correct
2597 // current incoming particles) to allow calculation of splitting
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
2606 void History::setupBeams() {
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;
2616 // Incoming partons to hard process are stored in slots 3 and 4.
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;
2625 // Save some info before clearing beams
2626 // Mothers of incoming partons companion code
2627 int motherPcompRes = -1;
2628 int motherMcompRes = -1;
2630 bool sameFlavP = false;
2631 bool sameFlavM = false;
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;
2640 sameFlavP = (state[inP].id() == mother->state[inMotherP].id());
2641 sameFlavM = (state[inM].id() == mother->state[inMotherM].id());
2643 motherPcompRes = (sameFlavP) ? beamA[0].companion() : -2;
2644 motherMcompRes = (sameFlavM) ? beamB[0].companion() : -2;
2647 // Append the current incoming particles to the beam
2651 // Get energy of incoming particles
2652 double Ep = 2. * state[inP].e();
2653 double Em = 2. * state[inM].e();
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();
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);
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);
2675 beamA.pickValSeaComp();
2677 beamA[0].companion(motherPcompRes);
2679 beamB.xfISR( 0, state[inM].id(), x2, scalePDF*scalePDF);
2681 beamB.pickValSeaComp();
2683 beamB[0].companion(motherMcompRes);
2688 //--------------------------------------------------------------------------
2690 // Calculate the PDF ratio used in the argument of the no-emission
2693 double History::pdfForSudakov() {
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;
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());
2705 // Done for pure FSR
2706 if (FSR) return 1.0;
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;
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;
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
2729 // Calculate pdf ratio
2730 double ratio = getPDFratio(side, true, idMother, xMother, scale,
2731 idDaughter, xDaughter, scale);
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);
2739 //--------------------------------------------------------------------------
2741 // Perform the clustering of the current state and return the
2743 // IN Clustering : rad,rec,emt triple to be clustered to two partons
2744 // OUT clustered state
2746 Event History::cluster( const Clustering & inSystem ) {
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;
2758 // Construct the clustered event
2759 Event NewEvent = Event();
2760 NewEvent.init("(hard process-modified)", particleDataPtr);
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] );
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();
2778 NewEvent.scaleSecond(mu);
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));
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();
2803 // Define radiator and recoiler back-to-back in the dipole
2804 // rest frame [=(state[rad]+state[emt])+state[rec] rest frame]
2807 Rad4mom.p( 0., 0., 0.5*eCMME, 0.5*eCMME);
2808 Rec4mom.p( 0., 0.,-0.5*eCMME, 0.5*eCMME);
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);
2819 RadBefore.p(Rad4mom);
2820 RecBefore.p(Rec4mom);
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();
2828 // Define radiator and recoiler back-to-back in the dipole
2829 // rest frame [=(state[rad]+state[emt])+state[rec] rest frame]
2832 Rad4mom.p( 0., 0., 0.5*eCMME, 0.5*eCMME);
2833 Rec4mom.p( 0., 0.,-0.5*eCMME, 0.5*eCMME);
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);
2844 RadBefore.p(Rad4mom);
2845 RecBefore.p(Rec4mom);
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]
2856 Rad4mom.p( 0., 0., 0.5*eCMME, 0.5*eCMME);
2857 Rec4mom.p( 0., 0.,-0.5*eCMME, 0.5*eCMME);
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);
2868 // Rescale recoiler momentum
2869 Rec4mom = 2.*state[Rec].p() - Rec4mom;
2871 RadBefore.p(Rad4mom);
2872 RecBefore.p(Rec4mom);
2874 // Set mass of initial recoiler to zero
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. );
2890 // Find side that radiates event (mother moving in
2891 // sign * p_z direction)
2893 int sign = state[Rad].pz() > 0 ? 1 : -1;
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);
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 );
2912 // Get mother and partner x values
2914 double x1 = 2. * pMother.e() / eCM;
2916 double x2 = 2. * pPartner.e() / eCM;
2918 pDaughter.p( pMother - pSister);
2919 pRecoiler.p( pPartner );
2921 // Find boost from event cm frame to rest frame of
2922 // of-shell daughter + on-shell recoiler
2923 RotBstMatrix from_CM_to_DR;
2925 from_CM_to_DR.toCMframe(pDaughter, pRecoiler);
2927 from_CM_to_DR.toCMframe(pRecoiler, pDaughter);
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 );
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.);
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 );
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();
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;
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));
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 ) );
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 );
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 );
2991 // Set momenta of particles to be attached to new event record
2992 RecBefore.p( pRecoiler );
2993 RadBefore.p( pDaughter );
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);
3003 // Append new recoiler and find new radiator colour
3004 NewEvent.append(RecBefore);
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,
3014 // Could happen if previous clustering produced several colour
3015 // singlett subsystems in the event
3020 // Build the clustered event
3021 Event outState = Event();
3022 outState.init("(hard process-modified)", particleDataPtr);
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();
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
3043 // Append first incoming particle
3044 if ( RecBefore.mother1() == 1) {
3045 outState.append( RecBefore );
3047 } else if ( RadBefore.mother1() == 1 ) {
3048 radPos = outState.append( RadBefore );
3051 // Find second incoming in input event
3053 for(int i=0; i < int(state.size()); ++i)
3054 if (state[i].mother1() == 1) in1 =i;
3055 outState.append( state[in1] );
3058 // Append second incoming particle
3059 if ( RecBefore.mother1() == 2) {
3060 outState.append( RecBefore );
3062 } else if ( RadBefore.mother1() == 2 ) {
3063 radPos = outState.append( RadBefore );
3066 // Find second incoming in input event
3068 for(int i=0; i < int(state.size()); ++i)
3069 if (state[i].mother1() == 2) in2 =i;
3071 outState.append( state[in2] );
3075 // Append new recoiler if not done already
3076 if (!recAppended && !RecBefore.isFinal()) {
3078 outState.append( RecBefore);
3080 // Append new radiator if not done already
3081 if (!radAppended && !RadBefore.isFinal()) {
3083 radPos = outState.append( RadBefore);
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] );
3091 if (!recAppended) outState.append(RecBefore);
3092 if (!radAppended) radPos = outState.append(RadBefore);
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] );
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]);
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
3119 PosDaughter1.push_back( daughter1);
3122 while(!outState[daughter1].isFinal() ) daughter1++;
3123 PosDaughter1.push_back( daughter1);
3126 PosDaughter2.push_back( daughter2);
3128 daughter2 = outState.size()-1;
3129 while(!outState[daughter2].isFinal() ) daughter2--;
3130 PosDaughter2.push_back( daughter2);
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]);
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);
3149 if (minParFinal == maxParFinal) maxParFinal = 0;
3150 outState[3].daughters(minParFinal,maxParFinal);
3151 outState[4].daughters(minParFinal,maxParFinal);
3153 // Update event properties
3154 outState.saveSize();
3155 outState.saveJunctionSize();
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
3165 // Hence, check for an intermediate coloured triplet resonance has been
3166 // colour-connected to the "old" radiator.
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() )
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() )
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() )
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() )
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))
3212 && ( radAcl == resCol || radAcl == resAcl));
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);
3229 // If event is not constructed properly, return false
3230 if ( !validEvent(outState) ) {
3235 // Remember position of reclustered radiator in state
3236 iReclusteredNew = radPos;
3242 //--------------------------------------------------------------------------
3244 // Function to get the flavour of the radiator before the splitting
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
3250 int History::getRadBeforeFlav(const int RadAfter, const int EmtAfter,
3251 const Event& event) {
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();
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)) ));
3269 // Final state gluon splitting
3270 if ( type == 1 && emtID == -radID && !colConnected )
3272 // Initial state s-channel gluon splitting
3273 if ( type ==-1 && radID == 21 )
3275 // Initial state t-channel gluon splitting
3276 if ( type ==-1 && emtID != 21 && radID != 21 && !colConnected )
3279 // Electroweak splittings splittings
3280 // Photon / Z radiation: Calculate invariant mass of system
3281 double m2final = (event[RadAfter].p()+ event[EmtAfter].p()).m2Calc();
3283 if ( emtID == 22 || emtID == 23 ) return radID;
3284 // Final state Photon splitting
3285 if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) <= 10. )
3287 // Final state Photon splitting
3288 if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) > 10. )
3290 // Initial state s-channel photon/ Z splitting
3291 if ( type ==-1 && (radID == 22 || radID == 23) )
3293 // Initial state t-channel photon / Z splitting: Always bookkeep as photon
3294 if ( type ==-1 && abs(emtID) < 10 && abs(radID) < 10 && colConnected )
3298 // Final state W+ splitting
3299 if ( emtID == 24 && radID < 0 ) return radID + 1;
3300 if ( emtID == 24 && radID > 0 ) return radID + 1;
3303 // Final state W- splitting
3304 if ( emtID ==-24 && radID < 0 ) return radID - 1;
3305 if ( emtID ==-24 && radID > 0 ) return radID - 1;
3311 //--------------------------------------------------------------------------
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
3323 bool History::connectRadiator( Particle& Radiator, const int RadType,
3324 const Particle& Recoiler, const int RecType,
3325 const Event& event ) {
3327 // Start filling radiator colour indices with dummy values
3328 Radiator.cols( -1, -1 );
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))
3344 // Set colour of antiquark radiator to zero
3346 for (int i = 0; i < event.size(); ++i) {
3347 int col = event[i].col();
3348 int acl = event[i].acol();
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());
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());
3361 } // end loop over particles in event record
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);
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))
3378 // Set anticolour of quark radiator to zero
3380 for (int i = 0; i < event.size(); ++i) {
3381 int col = event[i].col();
3382 int acl = event[i].acol();
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());
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());
3395 } // end loop over particles in event record
3397 } // end distinction between fsr / fsr+initial recoiler / isr
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
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();
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());
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());
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());
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());
3433 } // end loop over particles in event record
3434 } // end cases of different radiator colour type
3436 // If either colour or anticolour has not been set, return false
3437 if (Radiator.col() < 0 || Radiator.acol() < 0) return false;
3442 //--------------------------------------------------------------------------
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
3448 // int iExclude2 : Identifier of second particle to be excluded
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]
3459 int History::FindCol(int col, int iExclude1, int iExclude2,
3460 const Event& event, int type, bool isHardIn) {
3462 bool isHard = isHardIn;
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 ) {
3476 if ( event[n].col() == col ) {
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 ) {
3496 if ( event[n].col() == col ) {
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);
3510 //--------------------------------------------------------------------------
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
3519 int History::FindParticle( const Particle& particle, const Event& event,
3520 bool checkStatus ) {
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() ) {
3535 if ( checkStatus && event[index].status() != particle.status() )
3541 //--------------------------------------------------------------------------
3543 // Function to get the colour of the radiator before the splitting
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
3550 int History::getRadBeforeCol(const int rad, const int emt,
3551 const Event& event) {
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) {
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();
3588 // Get reconstructed quark colours
3589 } else if ( radBeforeFlav != 21 && radBeforeFlav > 0) {
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();
3624 // Other particles are assumed uncoloured
3629 return radBeforeCol;
3633 //--------------------------------------------------------------------------
3635 // Function to get the anticolour of the radiator before the splitting
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
3642 int History::getRadBeforeAcol(const int rad, const int emt,
3643 const Event& event) {
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) {
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();
3680 // Get reconstructed anti-quark colours
3681 } else if ( radBeforeFlav != 21 && radBeforeFlav < 0) {
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();
3716 // Other particles are considered uncoloured
3721 return radBeforeAcl;
3725 //--------------------------------------------------------------------------
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
3733 int History::getColPartner(const int in, const Event& event) {
3735 if (event[in].col() == 0) return 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
3742 partner = FindCol(event[in].col(),in,0,event,2,true);
3748 //--------------------------------------------------------------------------
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
3756 int History::getAcolPartner(const int in, const Event& event) {
3758 if (event[in].acol() == 0) return 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
3765 partner = FindCol(event[in].acol(),in,0,event,1,true);
3771 //--------------------------------------------------------------------------
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.
3782 vector<int> History::getReclusteredPartners(const int rad, const int emt,
3783 const Event& event) {
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);
3791 vector<int> partners;
3793 // Start with FSR clusterings
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);
3822 // Start with ISR clusterings
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);
3857 //--------------------------------------------------------------------------
3859 // Function to extract a chain of colour-connected partons in
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
3876 bool History::getColSinglet( const int flavType, const int iParton,
3877 const Event& event, vector<int>& exclude, vector<int>& colSinglet) {
3879 // If no possible flavour to start from has been found
3880 if (iParton < 0) return false;
3882 // If no further partner has been found in a previous iteration,
3883 // and the whole final state has been excluded, we're done
3886 // Count number of final state partons
3888 for(int i=0; i < int(event.size()); ++i)
3889 if ( event[i].isFinal() && event[i].colType() != 0)
3892 // Get number of initial state partons in the list of
3894 int nExclude = int(exclude.size());
3895 int nInitExclude = 0;
3896 if (!event[exclude[2]].isFinal())
3898 if (!event[exclude[3]].isFinal())
3901 // If the whole final state has been considered, return
3902 if (nFinal == nExclude - nInitExclude)
3909 // Declare colour partner
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
3917 colP = getColPartner(iParton,event);
3918 // When starting out from an antiquark line, follow the anticolour lines
3920 colP = getAcolPartner(iParton,event);
3922 // Do not count excluded partons twice
3923 for(int i = 0; i < int(exclude.size()); ++i)
3924 if (colP == exclude[i])
3928 return getColSinglet(flavType,colP,event,exclude,colSinglet);
3932 //--------------------------------------------------------------------------
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
3939 bool History::isColSinglet( const Event& event,
3940 vector<int> system ) {
3942 // Check if system forms a colour singlet
3943 for(int i=0; i < int(system.size()); ++i ) {
3944 // Match quark and gluon colours
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
3952 && event[system[i]].col() == event[system[j]].acol()) {
3953 // Remove index and break
3959 // Match antiquark and gluon anticolours
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
3967 && event[system[i]].acol() == event[system[j]].col()) {
3968 // Remove index and break
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 )
3989 //--------------------------------------------------------------------------
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
3998 bool History::isFlavSinglet( const Event& event,
3999 vector<int> system, int flav) {
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.
4010 if ( event[i].idAbs() != 21
4011 && event[i].idAbs() != 22
4012 && event[i].idAbs() != 23
4013 && event[i].idAbs() != 24
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
4021 if (abs(flav) > 0 && event[system[i]].idAbs() != flav)
4023 // Remove index and break
4028 // If flavour of outgoing and incoming partons match,
4029 // remove both partons and continue.
4031 if ( event[i].idAbs() != 21
4032 && event[i].idAbs() != 22
4033 && event[i].idAbs() != 23
4034 && event[i].idAbs() != 24
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
4042 if (abs(flav) > 0 && event[system[i]].idAbs() != flav)
4044 // Remove index and break
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 )
4066 //--------------------------------------------------------------------------
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
4073 bool History::allowedClustering( int rad, int emt, int rec, int partner,
4074 const Event& event ) {
4077 bool allowed = true;
4079 // CONSTRUCT SOME PROPERTIES FOR LATER INVESTIGATION
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);
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) )
4102 // Count coloured final state partons in event, excluding
4103 // rad, rec, emt and hard process
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) )
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 )
4120 // Get number of non-charged final state particles
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)))
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) )
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++;
4153 // Get number of initial quarks
4154 int nInitialQuark = 0;
4156 if (event[rec].isFinal()) {
4157 if (event[3].isQuark()) nInitialQuark++;
4158 if (event[4].isQuark()) nInitialQuark++;
4160 int iOtherIn = (rec == 3) ? 4 : 3;
4161 if (event[rec].isQuark()) nInitialQuark++;
4162 if (event[iOtherIn].isQuark()) nInitialQuark++;
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++;
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.
4177 for(int i=0; i < int(event.size()); ++i)
4178 if (event[i].idAbs() == 6)
4181 // BEGIN CHECKING THE CLUSTERING
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 ) {
4192 int colP = getColPartner(i,event);
4193 int aclP = getAcolPartner(i,event);
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);
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)
4209 // If triple forms colour singlett, check that resulting state
4210 // matches hard core process
4213 if (isSing && (abs(radBeforeFlav)<10 && event[rec].isQuark()) )
4216 // Never recluster any outgoing partons of the core V -> qqbar' splitting!
4217 if ( mergingHooksPtr->hardProcess.matchesAnyOutgoing(emt,event) )
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 )
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)
4232 if ( (nInitialQuark + nFinalQuark + nFinalQuarkExc)%2 != 0 )
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)
4246 // No problems with gluon radiation
4247 if (event[emt].id() == 21)
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() );
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() );
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);
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() ) {
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])) {
4291 outgoingParticles[j] = 99;
4293 if (!matchOut) iOutQuarkFlav.push_back(i);
4297 // Save number of potentially dangerous quarks
4298 int nInQuarkFlav = int(iInQuarkFlav.size());
4299 int nOutQuarkFlav = int(iOutQuarkFlav.size());
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)
4311 // If there are no ambiguities in qqbar pairs, return
4312 if (nInQuarkFlav + nOutQuarkFlav == 0)
4316 // Save all quarks and gluons that will not change colour
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) ) {
4327 partons.push_back(i);
4328 // Split into components
4329 if (event[i].colType() == 2)
4331 else if (event[i].colType() == 1)
4333 else if (event[i].colType() == -1)
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()) );
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 ) {
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
4355 colours.push_back(radBeforeCol);
4356 anticolours.push_back(radBeforeAcl);
4358 colours.push_back(radBeforeAcl);
4359 anticolours.push_back(radBeforeCol);
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());
4367 colours.push_back(event[gluon[i]].acol());
4368 anticolours.push_back(event[gluon[i]].col());
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]) {
4382 // If all gluon anticolours and all colours matched, disallow
4384 bool allMatched = true;
4385 for(int i=0; i < int(colours.size()); ++i)
4386 if (colours[i] != 0)
4388 for(int i=0; i < int(anticolours.size()); ++i)
4389 if (anticolours[i] != 0)
4395 // Now add the colours of the hard process, and check if all
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());
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());
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)
4411 for(int j=0; j < int(anticolours.size()); ++j)
4412 if (colours[i] > 0 && anticolours[j] > 0
4413 && colours[i] == anticolours[j]) {
4418 // Check if clustering would produce the hard process
4420 for ( int i=0; i < int(quark.size()); ++i )
4421 if ( !mergingHooksPtr->hardProcess.matchesAnyOutgoing( quark[i],
4424 for ( int i=0; i < int(antiq.size()); ++i )
4425 if ( !mergingHooksPtr->hardProcess.matchesAnyOutgoing( antiq[i],
4428 for(int i=0; i < int(gluon.size()); ++i)
4429 if ( event[gluon[i]].isFinal() )
4434 // If all colours are matched now, and since we have more quarks than
4435 // present in the hard process, disallow the clustering
4437 for(int i=0; i < int(colours.size()); ++i)
4438 if (colours[i] != 0)
4440 for(int i=0; i < int(anticolours.size()); ++i)
4441 if (anticolours[i] != 0)
4444 if (allMatched && nNotInHard > 0)
4451 if (isFSRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) {
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)
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)
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() ) {
4478 // This clustering is allowed if there is no colour in the
4480 if (nInitialPartons > 0)
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() ) {
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.
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];
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;
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);
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);
4556 // If no "recovery clustering" is possible, reject clustering
4558 && (!isColSingIn1 || !isFlavSingIn1
4559 || !isColSingIn2 || !isFlavSingIn2))
4564 // Next-to-easiest problem 2:
4565 // FINAL STATE Q-QBAR CLUSTERING DISCONNECTS SINGLETT SUBSYSTEM WITH
4566 // FINAL STATE Q-QBAR PAIR FROM GRAPH
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
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
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) {
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);
4603 if (isColSing && isFlavSing) {
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);
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);
4619 // If all final quarks are of identical flavour,
4620 // no possible clustering should be discriminated.
4621 // Otherwise, disallow
4622 if (!isFullColSing || !isFullFlavSing)
4629 if (isISRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) {
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)
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)
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() ) {
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() ))
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;
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
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);
4698 // Check if the incoming parton not
4699 // connected to the radiator is connected to a
4700 // colour and flavour singlet
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);
4713 // If no "recovery clustering" is possible, reject clustering
4715 && (!isColSingIn1 || !isFlavSingIn1
4716 || !isColSingIn2 || !isFlavSingIn2))
4727 //--------------------------------------------------------------------------
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
4733 // particles considered for clustering
4734 // Event event : Reference event
4736 bool History::isSinglett( int rad, int emt, int rec, const Event& event ) {
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;
4746 bool isSing = false;
4748 if ( ( recType == -1
4749 && radCol + emtCol == recCol && radAcl + emtAcl == recAcl)
4751 && radCol + emtCol == recAcl && radAcl + emtAcl == recCol) )
4758 //--------------------------------------------------------------------------
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
4767 bool History::validEvent( const Event& event ) {
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;
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;
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;
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;
4810 return (validColour && validCharge);
4814 //--------------------------------------------------------------------------
4816 // Function to check whether two clusterings are identical, used
4817 // for finding the history path in the mother -> children direction
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()) );
4827 //--------------------------------------------------------------------------
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
4833 double History::choseHardScale( const Event& event ) const {
4836 double mHat = (event[3].p() + event[4].p()).mCalc();
4838 // Find number of final state particles and bosons
4843 for(int i = 0; i < event.size(); ++i)
4844 if ( event[i].isFinal() ) {
4846 // Remember final state unstable bosons
4847 if ( event[i].idAbs() == 23
4848 || event[i].idAbs() == 24 ) {
4851 mBos += event[i].m();
4853 } else if ( abs(event[i].status()) == 22
4854 && ( event[i].idAbs() == 23
4855 || event[i].idAbs() == 24 )) {
4857 mBos += event[i].m(); // Real mass
4860 // Return averaged boson masses
4861 if ( nBosons > 0 && (nFinal + nFinBos*2) <= 3)
4862 return (mBos / double(nBosons));
4867 //--------------------------------------------------------------------------
4869 // If the state has an incoming hadron return the flavour of the
4870 // parton entering the hard interaction. Otherwise return 0
4872 int History::getCurrentFlav(const int side) const {
4873 int in = (side == 1) ? 3 : 4;
4874 return state[in].id();
4877 //--------------------------------------------------------------------------
4879 double History::getCurrentX(const int side) const {
4880 int in = (side == 1) ? 3 : 4;
4881 return ( 2.*state[in].e()/state[0].e() );
4884 //--------------------------------------------------------------------------
4886 double History::getCurrentZ(const int rad,
4887 const int rec, const int emt) const {
4889 int type = state[rad].isFinal() ? 1 : -1;
4893 // Construct 2->3 variables for FSR
4894 Vec4 sum = state[rad].p() + state[rec].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
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());
4913 //--------------------------------------------------------------------------
4915 // Function to compute "pythia pT separation" from Particle input
4917 double History::pTLund(const Particle& RadAfterBranch,
4918 const Particle& EmtAfterBranch,
4919 const Particle& RecAfterBranch, int ShowerType) {
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)
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.;
4952 return sqrt(pTpyth);
4955 //--------------------------------------------------------------------------
4957 // Function to return the position of the initial line before (or after)
4958 // a single (!) splitting.
4960 int History::posChangedIncoming(const Event& event, bool before) {
4962 // Check for initial state splittings.
4963 // Consider a splitting to exist if both mother and sister were found.
4966 for (int i =0; i < event.size(); ++i)
4967 if (event[i].status() == 43) {
4973 if (iSister > 0) iMother = event[iSister].mother1();
4975 // Initial state splitting has been found.
4976 if (iSister > 0 && iMother > 0) {
4978 // Find flavour, mother flavour
4979 int flavSister = event[iSister].id();
4980 int flavMother = event[iMother].id();
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)
4993 // Find initial state (!) daughter
4995 for (int i =0; i < event.size(); ++i)
4996 if ( !event[i].isFinal()
4997 && event[i].mother1() == iMother
4998 && event[i].id() == flavDaughter )
5001 // Done for initial state splitting.
5002 if ( !before ) return iMother;
5003 else return iDaughter;
5007 // Check for final state splittings with initial state recoiler.
5008 // Consider a splitting to exist if both mother and daughter were found.
5011 for (int i =0; i < event.size(); ++i)
5012 if ( abs(event[i].status()) == 53 || abs(event[i].status()) == 54) {
5019 if (iMother > 0) iDaughter = event[iMother].daughter1();
5021 // Done if final state splitting has been found.
5022 if (iDaughter > 0 && iMother > 0) {
5024 // Done for final state splitting.
5025 if ( !before ) return iMother;
5026 else return iDaughter;
5030 // If no splitting has been found, return zero.
5035 //==========================================================================
5037 } // end namespace Pythia8