10 #include "Pythia8/History.h"
11 #include "Pythia8/SharedPointers.h"
27 void Clustering::list()
const {
28 cout <<
" emt " << emitted
30 <<
" rec " << recoiler
31 <<
" partner " << partner
32 <<
" pTscale " << pTscale << endl;
54 const int History::NTRIAL = 1;
71 History::History(
int depthIn,
75 MergingHooksPtr mergingHooksPtrIn,
78 ParticleData* particleDataPtrIn,
80 PartonLevel* showersIn,
82 bool isOrdered =
true,
83 bool isStronglyOrdered =
true,
84 bool isAllowed =
true,
85 bool isNextInInput =
true,
94 foundOrderedPath(false),
95 foundStronglyOrderedPath(false),
96 foundAllowedPath(false),
97 foundCompletePath(false),
99 nextInInput(isNextInInput),
105 mergingHooksPtr(mergingHooksPtrIn),
108 particleDataPtr(particleDataPtrIn),
111 coupSMPtr(coupSMPtrIn),
121 if (mother && mergingHooksPtr->includeRedundant()) prob *= pdfForSudakov();
126 double acoll = (mother->state[clusterIn.emittor].isFinal())
127 ? mergingHooksPtr->herwigAcollFSR()
128 : mergingHooksPtr->herwigAcollISR();
129 sumScalarPT = mother->sumScalarPT + acoll*scale;
134 if ( mother ) iReclusteredOld = mother->iReclusteredNew;
137 int nFinalP = 0, nFinalW = 0, nFinalZ = 0;
138 int nL = 0, nA= 0, nH = 0;
139 for (
int i = 0; i < int(state.size()); ++i )
140 if ( state[i].isFinal() ) {
141 if ( state[i].colType() != 0 )
143 if ( state[i].idAbs() == 23 )
145 if ( state[i].idAbs() == 24 )
147 if ( state[i].idAbs() < 20 && state[i].idAbs() > 10)
149 if ( state[i].idAbs() == 22)
151 if ( state[i].idAbs() == 23
152 || state[i].idAbs() == 24
153 || state[i].idAbs() == 25)
156 if ( mergingHooksPtr->doWeakClustering()
157 && nFinalP == 2 && nFinalW == 0 && nFinalZ == 0) depth = 0;
162 bool qcd = ( nFinalP > mergingHooksPtr->hardProcess->nQuarksOut() );
166 vector<Clustering> clusterings;
167 if ( qcd && depth > 0 ) clusterings = getAllQCDClusterings();
169 bool dow = ( mergingHooksPtr->doWeakClustering()
170 && nFinalP > 1 && nFinalW+nFinalZ > 0 );
173 vector<Clustering> clusteringsEW;
175 if ( depth > 0 && dow )
176 clusteringsEW = getAllEWClusterings();
177 if ( !clusteringsEW.empty() ) {
178 clusterings.insert( clusterings.end(), clusteringsEW.begin(),
179 clusteringsEW.end() );
183 vector<Clustering> clusteringsSQCD;
184 if ( depth > 0 && mergingHooksPtr->doSQCDClustering() )
185 clusteringsSQCD = getAllSQCDClusterings();
186 if ( !clusteringsSQCD.empty() )
187 clusterings.insert( clusterings.end(), clusteringsSQCD.begin(),
188 clusteringsSQCD.end() );
192 if ( clusterings.empty() ) {
194 prob *= hardProcessME(state);
195 if (registerPath( *
this, isOrdered, isStronglyOrdered, isAllowed,
196 depth == 0 )) updateMinDepth(depth);
203 multimap<double, Clustering *> sort;
204 for (
unsigned int i = 0; i < clusterings.size(); ++i) {
205 double t = clusterings[i].pT();
217 sort.insert(make_pair(index, &clusterings[i]));
220 for ( multimap<double, Clustering *>::iterator it = sort.begin();
221 it != sort.end(); ++it ) {
223 double t = it->second->pT();
226 bool stronglyOrdered = isStronglyOrdered;
227 if ( mergingHooksPtr->enforceStrongOrdering()
228 && ( !stronglyOrdered
230 mergingHooksPtr->scaleSeparationFactor()*scale ) ))) {
231 if ( onlyStronglyOrderedPaths() )
continue;
232 stronglyOrdered =
false;
236 bool ordered = isOrdered;
237 if ( mergingHooksPtr->orderInRapidity()
238 && mergingHooksPtr->orderHistories() ) {
240 double z = getCurrentZ((*it->second).emittor,
241 (*it->second).recoiler,(*it->second).emitted,
242 (*it->second).flavRadBef);
244 double zOld = (!mother) ? 0. : mother->getCurrentZ(clusterIn.emittor,
245 clusterIn.recoiler,clusterIn.emitted,
246 clusterIn.flavRadBef);
249 if ( !ordered || ( mother && (t < scale
250 || t < pow(1. - z,2) / (z * (1. - zOld ))*scale ))) {
251 if ( onlyOrderedPaths() )
continue;
254 }
else if ( mergingHooksPtr->orderHistories() ) {
258 if ( !ordered || ( mother && (t < scale) ) ) {
259 if ( depth >= minDepth() && onlyOrderedPaths() && onlyAllowedPaths() )
266 bool doCut = mergingHooksPtr->canCutOnRecState()
267 || mergingHooksPtr->allowCutOnRecState();
268 bool allowed = isAllowed;
270 && mergingHooksPtr->doCutOnRecState(cluster(*it->second)) ) {
271 if ( onlyAllowedPaths() )
continue;
276 double p = getProb(*it->second);
277 if (abs(p)*prob < 1e-10*probMax())
continue;
278 updateProbMax(abs(p)*prob,depth==0);
282 children.push_back(
new History(depth - 1, t, cluster(*it->second),
283 *it->second, mergingHooksPtr, beamA, beamB, particleDataPtr,
284 infoPtr, showers, coupSMPtr, ordered, stronglyOrdered, allowed,
285 true, prob*p,
this ));
294 bool History::projectOntoDesiredHistories() {
296 return trimHistories();
311 vector<double> History::weightCKKWL(PartonLevel* trial, AlphaStrong * asFSR,
312 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN) {
314 if ( mergingHooksPtr->canCutOnRecState() && !foundAllowedPath ) {
315 string message=
"Warning in History::weightCKKWL: No allowed history";
316 message+=
" found. Using disallowed history.";
317 infoPtr->errorMsg(message);
319 if ( mergingHooksPtr->orderHistories() && !foundOrderedPath ) {
320 string message=
"Warning in History::weightCKKWL: No ordered history";
321 message+=
" found. Using unordered history.";
322 infoPtr->errorMsg(message);
324 if ( mergingHooksPtr->canCutOnRecState()
325 && mergingHooksPtr->orderHistories()
326 && !foundAllowedPath && !foundOrderedPath ) {
327 string message=
"Warning in History::weightCKKWL: No allowed or ordered";
328 message+=
" history found.";
329 infoPtr->errorMsg(message);
333 double asME = infoPtr->alphaS();
334 double aemME = infoPtr->alphaEM();
335 double maxScale = (foundCompletePath) ? infoPtr->eCM()
336 : mergingHooksPtr->muFinME();
339 History * selected = select(RN);
342 selected->setScalesInHistory();
344 int nWgts = mergingHooksPtr->nWgts;
346 vector<double> sudakov( nWgts, 1. );
347 vector<double> asWeight( nWgts, 1. );
348 vector<double> aemWeight( nWgts, 1. );
349 vector<double> pdfWeight( nWgts, 1. );
352 sudakov = selected->weightTree( trial, asME, aemME, maxScale,
353 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
354 aemWeight, pdfWeight );
357 int njetsMaxMPI = mergingHooksPtr->nMinMPI();
358 vector<double> mpiwt = selected->weightTreeEmissions( trial, -1, 0,
359 njetsMaxMPI, maxScale );
362 bool resetScales = mergingHooksPtr->resetHardQRen();
368 && mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0) {
370 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
371 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
372 for (
double& asW: asWeight) asW *= pow2(runningCoupling);
373 }
else if (mergingHooksPtr->doWeakClustering()
374 && isQCD2to2(selected->state)) {
376 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
377 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
378 for (
double& asW: asWeight) asW *= pow2(runningCoupling);
382 if (mergingHooksPtr->doWeakClustering() && isEW2to1(selected->state)) {
384 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
385 double runningCoupling = (*aemFSR).alphaEM(newQ2Ren) / aemME;
386 for (
double& aemW: aemWeight) aemW *= runningCoupling;
393 && mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) {
395 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
396 double runningCoupling =
397 (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME;
398 for (
double& asW: asWeight) asW *= runningCoupling;
403 for (
int iVar = 0; iVar < nWgts; ++iVar)
404 ret.push_back(sudakov[iVar]*asWeight[iVar]*aemWeight[iVar]*pdfWeight[iVar]*
415 vector<double> History::weightNL3Loop(PartonLevel* trial,
double RN ) {
417 if ( mergingHooksPtr->canCutOnRecState() && !foundAllowedPath ) {
418 string message=
"Warning in History::weightNL3Loop: No allowed history";
419 message+=
" found. Using disallowed history.";
420 infoPtr->errorMsg(message);
424 History * selected = select(RN);
426 selected->setScalesInHistory();
429 int nWgts = mergingHooksPtr->nWgts;
430 vector<double> wt( nWgts, 1. );
433 double maxScale = (foundCompletePath) ? infoPtr->eCM()
434 : mergingHooksPtr->muFinME();
435 int njetsMaxMPI = mergingHooksPtr->nMinMPI();
436 vector<double> mpiwt = selected->weightTreeEmissions( trial, -1, 0,
437 njetsMaxMPI, maxScale );
447 vector<double> History::weightNL3First(PartonLevel* trial, AlphaStrong* asFSR,
448 AlphaStrong* asISR, AlphaEM* , AlphaEM* ,
double RN,
452 double asME = infoPtr->alphaS();
453 double muR = mergingHooksPtr->muRinME();
454 double maxScale = (foundCompletePath)
456 : mergingHooksPtr->muFinME();
459 History * selected = select(RN);
461 selected->setScalesInHistory();
463 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
466 double kFactor = asME * mergingHooksPtr->k1Factor(nSteps);
470 double wt = 1. + kFactor;
473 double wtFirst = selected->weightFirst(trial,asME, muR, maxScale, asFSR,
477 double startingScale = (selected->mother) ? state.scale()
484 double nWeight1 = 0.;
485 for(
int i=0; i < NTRIAL; ++i) {
487 vector<double> unresolvedEmissionTerm = countEmissions( trial,
488 startingScale, mergingHooksPtr->tms(), 2, asME, asFSR, asISR, 3,
490 nWeight1 += unresolvedEmissionTerm[1];
493 wtFirst += nWeight1/double(NTRIAL);
496 int nWgts = mergingHooksPtr->nWgts;
497 vector<double> wtVec({wt+wtFirst});
500 for (
int iVar = 1; iVar < nWgts; ++iVar) {
501 double asFix = asFSR->alphaS(pow2(muR*mergingHooksPtr->
502 muRVarFactors[iVar-1])) / asME;
503 wtVec.push_back( wt + asFix*wtFirst );
507 for (
int iVar = 1; iVar < nWgts; ++iVar ) {
508 double corrFac = std::pow(asFSR->alphaS(pow2(muR*mergingHooksPtr->
509 muRVarFactors[iVar-1])) / asME, nSteps);
510 wtVec[iVar] *= corrFac;
520 vector<double> History::weightNL3Tree(PartonLevel* trial, AlphaStrong*
521 asFSR, AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
524 return weightCKKWL( trial, asFSR, asISR, aemFSR, aemISR, RN);
529 vector<double> History::weightUMEPSTree(PartonLevel* trial, AlphaStrong*
530 asFSR, AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
533 return weightCKKWL( trial, asFSR, asISR, aemFSR, aemISR, RN);
540 vector<double> History::weightUMEPSSubt(PartonLevel* trial, AlphaStrong*
541 asFSR, AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
545 double asME = infoPtr->alphaS();
546 double aemME = infoPtr->alphaEM();
547 double maxScale = (foundCompletePath) ? infoPtr->eCM()
548 : mergingHooksPtr->muFinME();
550 History * selected = select(RN);
552 selected->setScalesInHistory();
555 int nWgts = mergingHooksPtr->nWgts;
556 vector<double> sudakov( nWgts, 1.);
557 vector<double> asWeight( nWgts, 1.);
558 vector<double> aemWeight( nWgts, 1.);
559 vector<double> pdfWeight( nWgts, 1.);
562 sudakov = selected->weightTree(trial, asME, aemME, maxScale,
563 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
564 aemWeight, pdfWeight);
567 int njetsMaxMPI = mergingHooksPtr->nMinMPI()+1;
568 vector<double> mpiwt = selected->weightTreeEmissions( trial, -1, 0,
569 njetsMaxMPI, maxScale );
572 bool resetScales = mergingHooksPtr->resetHardQRen();
577 && mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0) {
579 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
580 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
581 for (
double& asW: asWeight) asW *= pow(runningCoupling,2);
588 && mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) {
590 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
591 double runningCoupling =
592 (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME;
593 for (
double& asW: asWeight) asW *= runningCoupling;
598 for (
int iVar = 0; iVar < nWgts; ++iVar)
599 ret.push_back(sudakov[iVar]*asWeight[iVar]*aemWeight[iVar]*pdfWeight[iVar]*
607 vector<double> History::weightUNLOPSTree(PartonLevel* trial, AlphaStrong*
608 asFSR, AlphaStrong* asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
611 if ( mergingHooksPtr->canCutOnRecState() && !foundAllowedPath ) {
612 string message=
"Warning in History::weightUNLOPSTree: No allowed";
613 message+=
" history found. Using disallowed history.";
614 infoPtr->errorMsg(message);
616 if ( mergingHooksPtr->orderHistories() && !foundOrderedPath ) {
617 string message=
"Warning in History::weightUNLOPSTree: No ordered";
618 message+=
" history found. Using unordered history.";
619 infoPtr->errorMsg(message);
621 if ( mergingHooksPtr->canCutOnRecState()
622 && mergingHooksPtr->orderHistories()
623 && !foundAllowedPath && !foundOrderedPath ) {
624 string message=
"Warning in History::weightUNLOPSTree: No allowed or";
625 message+=
" ordered history found.";
626 infoPtr->errorMsg(message);
630 double asME = infoPtr->alphaS();
631 double aemME = infoPtr->alphaEM();
632 double maxScale = (foundCompletePath) ? infoPtr->eCM()
633 : mergingHooksPtr->muFinME();
635 History * selected = select(RN);
637 selected->setScalesInHistory();
640 int nWgts = mergingHooksPtr->nWgts;
641 vector<double> asWeight( nWgts, 1. );
642 vector<double> aemWeight( nWgts, 1. );
643 vector<double> pdfWeight( nWgts, 1. );
646 vector<double> wt( nWgts, 1.);
647 if (depthIn < 0) wt = selected->weightTree(trial, asME, aemME, maxScale,
648 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
649 aemWeight, pdfWeight);
651 wt = selected->weightTreeEmissions( trial, 1, 0, depthIn, maxScale );
653 asWeight = selected->weightTreeAlphaS( asME, asFSR, asISR, depthIn);
654 aemWeight = selected->weightTreeAlphaEM( aemME, aemFSR, aemISR, depthIn);
655 pdfWeight = selected->weightTreePDFs( maxScale,
656 selected->clusterIn.pT(), depthIn);
661 int njetsMaxMPI = mergingHooksPtr->nMinMPI();
662 vector<double> mpiwt = selected->weightTreeEmissions( trial, -1, 0,
663 njetsMaxMPI, maxScale );
666 bool resetScales = mergingHooksPtr->resetHardQRen();
671 && mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0) {
673 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
674 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
675 for (
double& asW: asWeight) asW *= pow(runningCoupling,2);
682 && mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) {
684 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
685 double runningCoupling =
686 (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME;
687 for (
double& asW: asWeight) asW *= runningCoupling;
692 for (
int iVar = 0; iVar < nWgts; ++iVar)
693 ret.push_back(wt[iVar]*asWeight[iVar]*aemWeight[iVar]*pdfWeight[iVar]*
698 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
699 double muR = mergingHooksPtr->muRinME();
700 for (
int iVar = 1; iVar < nWgts; ++iVar)
701 asWeight[iVar] *= std::pow((*asFSR).alphaS(muR*muR) /
702 (*asFSR).alphaS(pow2(muR*mergingHooksPtr->muRVarFactors[iVar-1])),
706 mergingHooksPtr->individualWeights.wtSave = wt;
707 mergingHooksPtr->individualWeights.asWeightSave = asWeight;
708 mergingHooksPtr->individualWeights.aemWeightSave = aemWeight;
709 mergingHooksPtr->individualWeights.pdfWeightSave = pdfWeight;
710 mergingHooksPtr->individualWeights.mpiWeightSave = mpiwt;
718 vector<double> History::weightUNLOPSLoop(PartonLevel* trial, AlphaStrong*
719 asFSR, AlphaStrong* asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
722 if (depthIn < 0)
return weightNL3Loop(trial, RN);
725 double asME = infoPtr->alphaS();
726 double aemME = infoPtr->alphaEM();
727 double maxScale = (foundCompletePath) ? infoPtr->eCM() :
728 mergingHooksPtr->muFinME();
731 History * selected = select(RN);
733 selected->setScalesInHistory();
736 int nWgts = mergingHooksPtr->nWgts;
737 vector<double> wt( nWgts, 1.);
738 vector<double> asWeight( nWgts, 1.);
739 vector<double> aemWeight( nWgts, 1.);
740 vector<double> pdfWeight( nWgts, 1.);
743 if (depthIn < 0) wt = selected->weightTree(trial, asME, aemME, maxScale,
744 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
745 aemWeight, pdfWeight);
747 wt = selected->weightTreeEmissions( trial, 1, 0, depthIn, maxScale );
749 asWeight = selected->weightTreeAlphaS( asME, asFSR, asISR, depthIn,
751 aemWeight = selected->weightTreeAlphaEM( aemME, aemFSR, aemISR, depthIn);
752 pdfWeight = selected->weightTreePDFs( maxScale,
753 selected->clusterIn.pT(), depthIn);
758 int njetsMaxMPI = mergingHooksPtr->nMinMPI();
759 vector<double> mpiwt = selected->weightTreeEmissions( trial, -1, 0,
760 njetsMaxMPI, maxScale );
763 bool resetScales = mergingHooksPtr->resetHardQRen();
768 && mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0) {
770 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
771 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
772 for (
double& asW: asWeight) asW *= pow(runningCoupling,2);
779 && mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) {
781 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
782 double runningCoupling =
783 (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME;
784 for (
double& asW: asWeight) asW *= runningCoupling;
789 for (
int iVar = 0; iVar < nWgts; ++iVar)
790 ret.push_back(wt[iVar]*asWeight[iVar]*aemWeight[iVar]*pdfWeight[iVar]*
794 mergingHooksPtr->individualWeights.wtSave = wt;
795 mergingHooksPtr->individualWeights.asWeightSave = asWeight;
796 mergingHooksPtr->individualWeights.aemWeightSave = aemWeight;
797 mergingHooksPtr->individualWeights.pdfWeightSave = pdfWeight;
798 mergingHooksPtr->individualWeights.mpiWeightSave = mpiwt;
806 vector<double> History::weightUNLOPSSubt(PartonLevel* trial, AlphaStrong*
807 asFSR, AlphaStrong* asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
811 History* selected = select(RN);
813 selected->setScalesInHistory();
816 double asME = infoPtr->alphaS();
817 double aemME = infoPtr->alphaEM();
818 double maxScale = (foundCompletePath)
820 : mergingHooksPtr->muFinME();
822 int nWgts = mergingHooksPtr->nWgts;
826 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
827 if ( nSteps == 2 && mergingHooksPtr->nRecluster() == 2
828 && ( !foundCompletePath
829 || !selected->allIntermediateAboveRhoMS( mergingHooksPtr->tms() )) )
830 return vector<double>( nWgts, 0. );
833 vector<double> asWeight( nWgts, 1.);
834 vector<double> aemWeight( nWgts, 1.);
835 vector<double> pdfWeight( nWgts, 1.);
837 vector<double> wt( nWgts, 1.);
839 wt = selected->weightTree(trial, asME, aemME, maxScale,
840 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
841 aemWeight, pdfWeight);
843 wt = selected->weightTreeEmissions( trial, 1, 0, depthIn, maxScale );
845 asWeight = selected->weightTreeAlphaS( asME, asFSR, asISR, depthIn);
846 aemWeight = selected->weightTreeAlphaEM( aemME, aemFSR, aemISR, depthIn);
847 pdfWeight = selected->weightTreePDFs( maxScale,
848 selected->clusterIn.pT(), depthIn);
853 int njetsMaxMPI = mergingHooksPtr->nMinMPI()+1;
854 vector<double> mpiwt = selected->weightTreeEmissions( trial, -1, 0,
855 njetsMaxMPI, maxScale );
859 if (mergingHooksPtr->nRecluster() == 2 )
860 ret = wt = asWeight = aemWeight = pdfWeight = mpiwt =
861 vector<double>( nWgts, 1. );
863 for (
int iVar = 0; iVar < nWgts; ++iVar)
864 ret.push_back(asWeight[iVar]*aemWeight[iVar]*pdfWeight[iVar]*
865 wt[iVar]*mpiwt[iVar]);
870 double muR = mergingHooksPtr->muRinME();
871 for (
int iVar = 1; iVar < nWgts; ++iVar)
872 asWeight[iVar] *= std::pow((*asFSR).alphaS(muR*muR) /
873 (*asFSR).alphaS(pow2(muR*mergingHooksPtr->muRVarFactors[iVar-1])),
877 mergingHooksPtr->individualWeights.wtSave = wt;
878 mergingHooksPtr->individualWeights.asWeightSave = asWeight;
879 mergingHooksPtr->individualWeights.aemWeightSave = aemWeight;
880 mergingHooksPtr->individualWeights.pdfWeightSave = pdfWeight;
881 mergingHooksPtr->individualWeights.mpiWeightSave = mpiwt;
890 vector<double> History::weightUNLOPSSubtNLO(PartonLevel* trial, AlphaStrong*
891 asFSR, AlphaStrong* asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
895 int nWgts = mergingHooksPtr->nWgts;
896 vector<double> wt( nWgts, 1. );
901 History * selected = select(RN);
903 selected->setScalesInHistory();
905 double maxScale = (foundCompletePath) ? infoPtr->eCM()
906 : mergingHooksPtr->muFinME();
907 int njetsMaxMPI = mergingHooksPtr->nMinMPI()+1;
908 vector<double> mpiwt = selected->weightTreeEmissions
909 ( trial, -1, 0, njetsMaxMPI, maxScale );
916 History* selected = select(RN);
918 selected->setScalesInHistory();
921 double asME = infoPtr->alphaS();
922 double aemME = infoPtr->alphaEM();
923 double maxScale = (foundCompletePath)
925 : mergingHooksPtr->muFinME();
929 double nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
930 if ( nSteps == 2 && mergingHooksPtr->nRecluster() == 2
931 && ( !foundCompletePath
932 || !selected->allIntermediateAboveRhoMS( mergingHooksPtr->tms() )) )
933 return vector<double>( nWgts, 0. );
936 vector<double> asWeight( nWgts, 1.);
937 vector<double> aemWeight( nWgts, 1.);
938 vector<double> pdfWeight( nWgts, 1.);
941 wt = selected->weightTree(trial, asME, aemME, maxScale,
942 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
943 aemWeight, pdfWeight);
945 wt = selected->weightTreeEmissions( trial, 1, 0, depthIn, maxScale );
947 asWeight = selected->weightTreeAlphaS( asME, asFSR, asISR, depthIn,
949 aemWeight = selected->weightTreeAlphaEM( aemME, aemFSR, aemISR, depthIn);
950 pdfWeight = selected->weightTreePDFs( maxScale,
951 selected->clusterIn.pT(), depthIn);
956 int njetsMaxMPI = mergingHooksPtr->nMinMPI()+1;
957 vector<double> mpiwt = selected->weightTreeEmissions( trial, -1, 0,
958 njetsMaxMPI, maxScale );
962 if (mergingHooksPtr->nRecluster() == 2 )
963 ret = wt = asWeight = aemWeight = pdfWeight = mpiwt =
964 vector<double>( nWgts, 1. );
966 for (
int iVar = 0; iVar < nWgts; ++iVar)
967 ret.push_back(asWeight[iVar]*aemWeight[iVar]*pdfWeight[iVar]*
968 wt[iVar]*mpiwt[iVar]);
972 mergingHooksPtr->individualWeights.wtSave = wt;
973 mergingHooksPtr->individualWeights.asWeightSave = asWeight;
974 mergingHooksPtr->individualWeights.aemWeightSave = aemWeight;
975 mergingHooksPtr->individualWeights.pdfWeightSave = pdfWeight;
976 mergingHooksPtr->individualWeights.mpiWeightSave = mpiwt;
986 vector<double> History::weightUNLOPSFirst(
int order, PartonLevel*
987 trial, AlphaStrong* asFSR, AlphaStrong* asISR, AlphaEM*, AlphaEM*,
988 double RN, Rndm* rndmPtr ) {
990 int nWgts = mergingHooksPtr->nWgts;
993 if ( order < 0 )
return vector<double>( nWgts, 0. );
996 double asME = infoPtr->alphaS();
998 double muR = mergingHooksPtr->muRinME();
999 double maxScale = (foundCompletePath)
1001 : mergingHooksPtr->muFinME();
1004 History * selected = select(RN);
1006 selected->setScalesInHistory();
1008 double nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
1011 double kFactor = asME * mergingHooksPtr->k1Factor(nSteps);
1017 vector<double> wtVec(nWgts, wt);
1021 if ( mergingHooksPtr->orderHistories() && foundOrderedPath )
1027 double wA = selected->weightFirstAlphaS( asME, muR, asFSR, asISR );
1029 double nWeight = 0.;
1030 for (
int i=0; i < NTRIAL; ++i ) {
1032 double wE = selected->weightFirstEmissions(trial,asME, maxScale,
1033 asFSR, asISR,
true,
true );
1037 double pscale = selected->clusterIn.pT();
1038 double wP = selected->weightFirstPDFs(asME, maxScale, pscale, rndmPtr);
1045 wtVec = vector<double>({wt + wA + nWeight/double(NTRIAL)});
1048 for (
int iVar = 1; iVar < nWgts; ++iVar) {
1049 double asFix = asFSR->alphaS(pow2(muR*mergingHooksPtr->
1050 muRVarFactors[iVar-1])) / asME;
1051 wtVec.push_back( wt + asFix*(wA + nWeight /
double(NTRIAL)) );
1056 mergingHooksPtr->individualWeights.bornAsVarFac =
1057 vector<double>( nWgts, 1. );
1058 for (
int iVar = 1; iVar < nWgts; ++iVar ) {
1059 double corrFac = std::pow(asFSR->alphaS(pow2(muR*mergingHooksPtr->
1060 muRVarFactors[iVar-1])) / asME, nSteps);
1061 wtVec[iVar] *= corrFac;
1062 mergingHooksPtr->individualWeights.bornAsVarFac[iVar] =
1067 if ( order <= 1 )
return wtVec;
1070 return vector<double>( nWgts, 0. );
1078 void History::getStartingConditions(
const double RN,
Event& outState ) {
1081 History * selected = select(RN);
1084 selected->setScalesInHistory();
1087 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
1090 if (!selected->mother) {
1092 for(
int i=0; i < int(state.size()); ++i)
1093 if ( state[i].isFinal()) nFinal++;
1095 state.scale(mergingHooksPtr->muF());
1101 if (mergingHooksPtr->getNumberOfClusteringSteps(state) == 0) {
1102 infoPtr->zNowISR(0.5);
1103 infoPtr->pT2NowISR(pow(state[0].e(),2));
1104 infoPtr->hasHistory(
true);
1107 infoPtr->zNowISR(selected->zISR());
1108 infoPtr->pT2NowISR(pow(selected->pTISR(),2));
1109 infoPtr->hasHistory(
true);
1115 double muf = state[0].e();
1116 for (
int i=0; i < state.size(); ++i )
1117 if ( state[i].isFinal()
1118 && ( state[i].colType() != 0 || state[i].id() == 22 ) ) {
1120 muf = min( muf, abs(state[i].mT()) );
1124 if ( nSteps == 0 && nFinalCol == 2
1125 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
1126 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) ) {
1128 for (
int i = 3;i < state.size();++i)
1129 state[i].scale(muf);
1133 if (nSteps == 0 && nFinalCol == 2 &&
1134 mergingHooksPtr->getProcessString().find(
"inc") != string::npos) {
1136 for (
int i = 3;i < state.size();++i)
1137 state[i].scale(muf);
1138 for (
int i=0; i < min(state.size(),outState.size()); ++i )
1139 state[i].pol(outState[i].pol());
1147 infoPtr->zNowISR(selected->zISR());
1148 infoPtr->pT2NowISR(pow(selected->pTISR(),2));
1149 infoPtr->hasHistory(
true);
1158 mergingHooksPtr->muMI(infoPtr->eCM());
1160 mergingHooksPtr->muMI(outState.scale());
1163 if (mergingHooksPtr->doWeakClustering()) setupSimpleWeakShower(0);
1165 mergingHooksPtr->setShowerStoppingScale( 0.0 );
1172 void History::printHistory(
const double RN ) {
1173 History * selected = select(RN);
1174 selected->printStates();
1182 void History::printStates() {
1184 cout << scientific << setprecision(6) <<
"Probability=" << prob << endl;
1190 double p = prob/mother->prob;
1191 cout << scientific << setprecision(6) <<
"Probability=" << p
1192 <<
" scale=" << clusterIn.pT() << endl;
1195 mother->printStates();
1204 bool History::getClusteredEvent(
const double RN,
int nSteps,
1208 History * selected = select(RN);
1212 selected->setScalesInHistory();
1215 if (nSteps > selected->nClusterings())
return false;
1218 outState = selected->clusteredState(nSteps-1);
1226 bool History::getFirstClusteredEventAboveTMS(
const double RN,
int nDesired,
1227 Event& process,
int& nPerformed,
bool doUpdate ) {
1230 int nTried = nDesired - 1;
1232 int nSteps = select(RN)->nClusterings();
1234 select(RN)->setScalesInHistory();
1241 dummy.init(
"(hard process-modified)", particleDataPtr );
1246 if ( !getClusteredEvent( RN, nSteps-nTried+1, dummy ) )
return false;
1247 if ( nTried >= nSteps )
break;
1250 }
while ( mergingHooksPtr->getNumberOfClusteringSteps(dummy) > 0
1251 && mergingHooksPtr->tmsNow( dummy) < mergingHooksPtr->tms() );
1254 if ( doUpdate ) process = dummy;
1257 if ( nTried > nSteps )
return false;
1259 nPerformed = nTried;
1262 mergingHooksPtr->nReclusterSave = nPerformed;
1264 if (mergingHooksPtr->getNumberOfClusteringSteps(state) == 0)
1265 mergingHooksPtr->muMI(infoPtr->eCM());
1267 mergingHooksPtr->muMI(state.scale());
1279 double History::getPDFratio(
int side,
bool forSudakov,
bool useHardPDFs,
1280 int flavNum,
double xNum,
double muNum,
1281 int flavDen,
double xDen,
double muDen) {
1284 if ( abs(flavNum) > 10 && flavNum != 21 )
return 1.0;
1285 if ( abs(flavDen) > 10 && flavDen != 21 )
return 1.0;
1288 double pdfRatio = 1.0;
1291 double pdfNum = 0.0;
1292 double pdfDen = 0.0;
1295 if ( useHardPDFs ) {
1298 pdfNum = mother->beamA.xfHard( flavNum, xNum, muNum*muNum);
1300 pdfNum = beamA.xfHard( flavNum, xNum, muNum*muNum);
1301 pdfDen = max(1e-10, beamA.xfHard( flavDen, xDen, muDen*muDen));
1304 pdfNum = mother->beamB.xfHard( flavNum, xNum, muNum*muNum);
1306 pdfNum = beamB.xfHard( flavNum, xNum, muNum*muNum);
1307 pdfDen = max(1e-10,beamB.xfHard( flavDen, xDen, muDen*muDen));
1314 pdfNum = mother->beamA.xfISR(0, flavNum, xNum, muNum*muNum);
1316 pdfNum = beamA.xfISR(0, flavNum, xNum, muNum*muNum);
1317 pdfDen = max(1e-10, beamA.xfISR(0, flavDen, xDen, muDen*muDen));
1321 pdfNum = mother->beamB.xfISR(0, flavNum, xNum, muNum*muNum);
1323 pdfNum = beamB.xfISR(0, flavNum, xNum, muNum*muNum);
1324 pdfDen = max(1e-10,beamB.xfISR(0, flavDen, xDen, muDen*muDen));
1329 if ( forSudakov && abs(flavNum) ==4 && abs(flavDen) == 4 && muDen == muNum
1330 && muNum < particleDataPtr->m0(4))
1331 pdfDen = pdfNum = 1.0;
1334 if ( pdfNum > 1e-15 && pdfDen > 1e-10 ) {
1335 pdfRatio *= pdfNum / pdfDen;
1336 }
else if ( pdfNum < pdfDen ) {
1338 }
else if ( pdfNum > pdfDen ) {
1354 void History::setScalesInHistory() {
1362 setScales(ident,
true);
1373 double History::getWeakProb() {
1374 vector<int> modes, fermionLines;
1376 return getWeakProb(modes, mom, fermionLines);
1385 double History::getWeakProb(vector<int> &mode, vector<Vec4> &mom,
1386 vector<int> fermionLines) {
1389 if (!mother)
return 1.;
1392 map<int,int> stateTransfer;
1393 findStateTransfer(stateTransfer);
1396 if (mode.empty()) setupWeakHard(mode,fermionLines,mom);
1399 vector<int> modeNew = updateWeakModes(mode, stateTransfer);
1400 vector<int> fermionLinesNew = updateWeakFermionLines(fermionLines,
1404 if (mother->state[clusterIn.emitted].idAbs() == 24 ||
1405 mother->state[clusterIn.emitted].idAbs() == 23)
1406 return getSingleWeakProb(modeNew, mom, fermionLinesNew) *
1407 mother->getWeakProb(modeNew, mom, fermionLinesNew);
1408 else return mother->getWeakProb(modeNew, mom, fermionLinesNew);
1413 double History::getSingleWeakProb(vector<int> &mode, vector<Vec4> &mom,
1414 vector<int> fermionLines) {
1417 double weakCoupling = 0.0;
1418 if (mother->state[clusterIn.emitted].idAbs() == 24) {
1420 if (clusterIn.spinRadBef == 1)
return 0.0;
1421 else if (clusterIn.spinRadBef == -1)
1422 weakCoupling = 4.*M_PI/ coupSMPtr->sin2thetaW()
1423 * coupSMPtr->V2CKMid(abs(clusterIn.flavRadBef),
1424 mother->state[clusterIn.emittor].idAbs());
1426 infoPtr->errorMsg(
"Warning in History::getSingleWeakProb: "
1427 "Spin not properly configurated. Skipping history");
1430 }
else if (mother->state[clusterIn.emitted].idAbs() == 23) {
1432 if (clusterIn.spinRadBef == 1)
1433 weakCoupling = 4.*M_PI*pow2(coupSMPtr->rf( abs(clusterIn.flavRadBef)))
1434 / (coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW()) ;
1435 else if (clusterIn.spinRadBef == -1)
1436 weakCoupling = 4.*M_PI*pow2(coupSMPtr->lf( abs(clusterIn.flavRadBef)))
1437 / (coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW()) ;
1439 infoPtr->errorMsg(
"Warning in History::getSingleWeakProb: "
1440 "Spin not properly configurated. Skipping history");
1444 infoPtr->errorMsg(
"Warning in History::getSingleWeakProb: "
1445 "Did not emit W/Z. Skipping history.");
1452 Vec4 pRadAft = mother->state[clusterIn.emittor].p();
1453 Vec4 pEmtAft = mother->state[clusterIn.emitted].p();
1454 Vec4 pRecAft = mother->state[clusterIn.recoiler].p();
1455 Vec4 pSum = pRadAft + pEmtAft + pRecAft;
1456 double m2sum = pSum.m2Calc();
1457 double Qsq = (pRadAft + pEmtAft).m2Calc();
1459 double m2Rad0 = pRadAft.m2Calc();
1460 double m2Emt0 = pEmtAft.m2Calc();
1461 double lambda13 = sqrt( pow2(Qsq - m2Rad0 - m2Emt0 ) - 4. * m2Rad0*m2Emt0 );
1462 double k1 = ( Qsq - lambda13 + (m2Emt0 - m2Rad0 ) ) / ( 2. * Qsq );
1463 double k3 = ( Qsq - lambda13 - (m2Emt0 - m2Rad0 ) ) / ( 2. * Qsq );
1465 double z = mother->getCurrentZ(clusterIn.emittor, clusterIn.recoiler,
1466 clusterIn.emitted, clusterIn.flavRadBef);
1467 double pT2 = pow2(clusterIn.pTscale);
1469 double x1 = 2. * pRadAft * pSum / m2sum;
1470 double x2 = 2. * pRecAft * pSum / m2sum;
1471 double x3 = 2. * pEmtAft * pSum / m2sum;
1474 if ( state[clusterIn.radBef].status() > 0) {
1476 if (mode[clusterIn.emittor] == 1) {
1478 double eCMME = pSum.mCalc();
1479 double r1 = mother->state[clusterIn.emittor].m() / eCMME;
1480 double r2 = mother->state[clusterIn.recoiler].m() / eCMME;
1481 double r3 = mother->state[clusterIn.emitted].m() / eCMME;
1482 double x1s = x1 * x1;
1483 double x2s = x2 * x2;
1484 double r1s = r1 * r1;
1485 double r2s = r2 * r2;
1486 double r3s = r3 * r3;
1487 double prop1 = 1. + r1s - r2s - x1;
1488 double prop2 = 1. + r2s - r1s - x2;
1489 double prop1s = prop1 * prop1;
1490 double prop2s = prop2 * prop2;
1491 double prop12 = prop1 * prop2;
1494 double jac = 1./(1.-z) * 1./pT2 * (1-x2+r2-r1)*(x3 - k1*(x1+x3))
1495 * (1.-x1+r1-r2) / x3;
1496 return jac * weakCoupling * ((2. * r3s * r3s + 2. * r3s *
1497 (x1 + x2) + x1s + x2s) / prop12 - r3s / prop1s - r3s / prop2s);
1502 Vec4 p1 = mother->state[clusterIn.emittor].p();
1503 Vec4 p2 = mother->state[clusterIn.recoiler].p();
1504 Vec4 p3 = mother->state[clusterIn.emitted].p();
1505 Vec4 radBef = state[clusterIn.radBef].p();
1506 Vec4 recBef = state[clusterIn.recBef].p();
1511 if (fermionLines[2] == clusterIn.emittor);
1512 else if (fermionLines[3] == clusterIn.emittor) swap(pIn1, pIn2);
1515 double scaleFactor2 = (p1 + p2 + p3).m2Calc() / (pIn1 + pIn2).m2Calc();
1516 double scaleFactor = sqrt(scaleFactor2);
1517 pIn1 *= scaleFactor;
1518 pIn2 *= scaleFactor;
1522 RotBstMatrix rot2to2frame;
1523 rot2to2frame.bstback(pIn1 + pIn2);
1524 pIn1.rotbst(rot2to2frame);
1525 pIn2.rotbst(rot2to2frame);
1526 p1.rotbst(rot2to2frame);
1527 p2.rotbst(rot2to2frame);
1528 p3.rotbst(rot2to2frame);
1529 recBef.rotbst(rot2to2frame);
1530 radBef.rotbst(rot2to2frame);
1533 RotBstMatrix rot2to3frame;
1534 rot2to3frame.bstback(p1 + p2 + p3);
1535 p1.rotbst(rot2to3frame);
1536 p2.rotbst(rot2to3frame);
1537 p3.rotbst(rot2to3frame);
1538 recBef.rotbst(rot2to3frame);
1539 radBef.rotbst(rot2to3frame);
1542 double sHat = (pIn1 + pIn2).m2Calc();
1543 double tHat = (radBef - pIn1).m2Calc();
1544 double uHat = (recBef - pIn1).m2Calc();
1545 double localProb = 0;
1546 double Q2 = pT2 / (z*(1.-z));
1547 double jac = 1./(4. * M_PI) * 1./( (1.-z) * z ) * sHat / (sHat - Q2)
1551 if (mode[clusterIn.emittor] == 2)
1552 localProb = simpleWeakShowerMEs.getMEqg2qgZ( pIn1, pIn2, p2, p3, p1)
1553 / simpleWeakShowerMEs.getMEqg2qg( sHat, tHat, uHat);
1554 else if (mode[clusterIn.emittor] == 3)
1555 localProb = simpleWeakShowerMEs.getMEqq2qqZ( pIn1, pIn2, p3, p2, p1)
1556 / simpleWeakShowerMEs.getMEqq2qq( sHat, tHat, uHat,
false);
1557 else if (mode[clusterIn.emittor] == 4)
1558 localProb = simpleWeakShowerMEs.getMEqq2qqZ( pIn1, pIn2, p3, p2, p1)
1559 / simpleWeakShowerMEs.getMEqq2qq( sHat, tHat, uHat,
true);
1561 string message=
"Warning in History::getSingleWeakProb: Wrong";
1562 message+=
" mode setup. Setting probability for path to zero.";
1563 infoPtr->errorMsg(message);
1567 localProb *= abs((-p3 + pIn1).m2Calc())
1568 / ((p3 + p1).m2Calc() + abs((-pIn1 + p3).m2Calc()));
1570 return jac * weakCoupling * localProb;
1576 if (mode[clusterIn.emittor] == 1) {
1577 Vec4 pIn1 = mother->state[clusterIn.emittor].p();
1578 Vec4 pIn2 = mother->state[clusterIn.recoiler].p();
1579 Vec4 p1 = mother->state[clusterIn.emitted].p();
1580 Vec4 p2 = pIn1 + pIn2 -p1;
1582 double sH = (pIn1 + pIn2).m2Calc();
1583 double tH = (p1 - pIn1).m2Calc();
1584 double uH = (p1 - pIn2).m2Calc();
1585 double m3s = p1.m2Calc();
1586 double m4s = p2.m2Calc();
1588 double jac = 1./sH * tH*uH / ( tH * (tH + uH) );
1589 return jac * weakCoupling * ((uH*uH + tH*tH + 2 * sH * (m3s + m4s))
1590 / (uH*tH) - m3s * m4s * (1/(tH*tH) + 1/(uH*uH)));
1595 Vec4 pIn1 = mother->state[clusterIn.emittor].p();
1596 Vec4 pIn2 = mother->state[clusterIn.recoiler].p();
1597 Vec4 p3 = mother->state[clusterIn.emitted].p();
1602 if (fermionLines[0] == clusterIn.emittor);
1603 else if (fermionLines[1] == clusterIn.emittor)
1607 Vec4 pDaughter, pRecoiler;
1608 int signBeam = (pIn1.pz() > 0.) ? 1 : -1;
1609 double eCM = state[0].e(), phi = 0;
1612 reverseBoostISR(pIn1, p3, pIn2, pDaughter, pRecoiler, signBeam,
1617 double scaleFactor2 = (pIn1 + pIn2 - p3).m2Calc() / (p1 + p2).m2Calc();
1618 double scaleFactor = sqrt(scaleFactor2);
1619 RotBstMatrix rot2to2frame;
1620 rot2to2frame.bstback(p1 + p2);
1621 p1.rotbst(rot2to2frame);
1622 p2.rotbst(rot2to2frame);
1628 Vec4 radBef = state[clusterIn.radBef].p();
1629 Vec4 recBef = state[clusterIn.recBef].p();
1631 RotBstMatrix rot2to2frameInc;
1632 rot2to2frameInc.bstback(radBef + recBef);
1633 radBef.rotbst(rot2to2frameInc);
1634 recBef.rotbst(rot2to2frameInc);
1635 double sHat = (p1 + p2).m2Calc();
1636 double tHat = (p1 - radBef).m2Calc();
1637 double uHat = (p1 - recBef).m2Calc();
1638 double localProb = 0;
1643 double jac = z / (4.*M_PI);
1645 if (mode[clusterIn.emittor] == 2)
1646 localProb = simpleWeakShowerMEs.getMEqg2qgZ(pIn1, pIn2, p2, p3, p1)
1647 / simpleWeakShowerMEs.getMEqg2qg(sHat, tHat, uHat);
1648 else if (mode[clusterIn.emittor] == 4)
1649 localProb = simpleWeakShowerMEs.getMEqq2qqZ(pIn1, pIn2, p3, p2, p1)
1650 / simpleWeakShowerMEs.getMEqq2qq(sHat, tHat, uHat,
true);
1651 else if (mode[clusterIn.emittor] == 3)
1652 localProb = simpleWeakShowerMEs.getMEqq2qqZ(pIn1, pIn2, p3, p2, p1)
1653 / simpleWeakShowerMEs.getMEqq2qq(sHat, tHat, uHat,
false);
1655 string message=
"Warning in History::getSingleWeakProb: Wrong";
1656 message+=
" mode setup. Setting probability for path to zero.";
1657 infoPtr->errorMsg(message);
1661 localProb *= (p3 + p1).m2Calc() / ( (p3 + p1).m2Calc()
1662 + abs((-pIn1 + p3).m2Calc()) );
1664 return jac * weakCoupling * localProb;
1673 bool History::checkWeakRecoils(map<int,int> &allowedRecoils,
bool isFirst) {
1674 if (!mother)
return true;
1679 if (state.size() != 8) {
1680 if (state[3].isQuark() || state[3].isLepton())
1681 allowedRecoils.insert(pair<int,int>(3,4));
1682 if (state[4].isQuark() || state[4].isLepton())
1683 allowedRecoils.insert(pair<int,int>(4,3));
1686 if (state[3].isQuark() || state[3].isLepton())
1687 allowedRecoils.insert(pair<int,int>(3,4));
1688 if (state[4].isQuark() || state[4].isLepton())
1689 allowedRecoils.insert(pair<int,int>(4,3));
1690 if (state[5].isQuark() || state[5].isLepton())
1691 allowedRecoils.insert(pair<int,int>(5,6));
1692 if (state[6].isQuark() || state[6].isLepton())
1693 allowedRecoils.insert(pair<int,int>(6,5));
1698 map<int,int> transfer;
1699 findStateTransfer(transfer);
1702 map<int,int> allowedRecoilsNew;
1703 for (map<int,int>::iterator it = allowedRecoils.begin();
1704 it != allowedRecoils.end(); ++it) {
1707 if (state[clusterIn.radBef].status() > 0) {
1709 if (it->first != clusterIn.radBef &&
1710 it->second != clusterIn.radBef)
1711 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1712 transfer[it->second]));
1714 else if (it->second == clusterIn.radBef) {
1716 if (state[clusterIn.recBef].isQuark() ||
1717 state[clusterIn.recBef].isLepton()) {
1718 if (mother->state[clusterIn.emittor].isQuark() ||
1719 mother->state[clusterIn.emittor].isLepton())
1720 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1721 clusterIn.emittor));
1723 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1724 clusterIn.emitted));
1728 double mEmittor = (mother->state[clusterIn.emittor].p() +
1729 mother->state[transfer[it->first]].p()).mCalc();
1730 double mEmitted = (mother->state[clusterIn.emitted].p() +
1731 mother->state[transfer[it->first]].p()).mCalc();
1732 if (mEmitted > mEmittor)
1733 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1734 clusterIn.emitted));
1736 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1737 clusterIn.emittor));
1741 if (mother->state[clusterIn.emittor].isQuark() ||
1742 mother->state[clusterIn.emittor].isLepton())
1743 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1744 transfer[it->second]));
1746 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emitted,
1747 transfer[it->second]));
1753 if (it->first != clusterIn.radBef &&
1754 it->second != clusterIn.radBef)
1755 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1756 transfer[it->second]));
1759 else if (it->second == clusterIn.radBef)
1760 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1761 clusterIn.emittor));
1766 if (mother->state[clusterIn.emittor].isQuark() ||
1767 mother->state[clusterIn.emittor].isLepton())
1768 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1769 clusterIn.recoiler));
1773 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1774 findISRRecoiler()));
1781 if ( (state[clusterIn.radBef].idAbs() == 22 ||
1782 state[clusterIn.radBef].idAbs() == 21) &&
1783 (mother->state[clusterIn.emittor].isQuark() ||
1784 mother->state[clusterIn.emittor].isLepton() ) ) {
1786 if (state[clusterIn.radBef].status() > 0) {
1787 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1788 clusterIn.emitted));
1789 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emitted,
1790 clusterIn.emittor));
1795 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1796 clusterIn.recoiler));
1797 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emitted,
1798 findISRRecoiler()));
1806 if (mother->state[clusterIn.emitted].idAbs() == 24 ||
1807 mother->state[clusterIn.emitted].idAbs() == 23)
1808 if ( clusterIn.recoiler != allowedRecoilsNew[clusterIn.emittor])
1812 return mother->checkWeakRecoils(allowedRecoilsNew);
1821 int History::findISRRecoiler() {
1823 int flavRad = mother->state[clusterIn.emitted].id();
1824 Vec4 pRad = mother->state[clusterIn.emitted].p();
1825 double mRad = mother->state[clusterIn.emitted].m();
1826 int iRad = clusterIn.emitted;
1828 double ppMin = 1E20;
1829 for (
int i = 0;i < mother->state.size(); ++i) {
1830 if (i == iRad)
continue;
1831 if (mother->state[i].isFinal() && mother->state[i].id() == - flavRad) {
1832 double ppNow = mother->state[i].p() * pRad
1833 - mother->state[i].m() - mRad;
1834 if (ppNow < ppMin) {
1840 if (iRec)
return iRec;
1843 for (
int i = 0;i < mother->state.size(); ++i) {
1844 if (i == iRad)
continue;
1845 if (mother->state[i].isFinal() && mother->state[i].idAbs() < 20) {
1846 double weakCoupNow = 1.;
1847 double ppNow = (mother->state[i].p() * pRad
1848 - mother->state[i].m() - mRad) / weakCoupNow;
1849 if (ppNow < ppMin) {
1855 if (iRec)
return iRec;
1858 for (
int i = 0;i < mother->state.size(); ++i) {
1859 if (i == iRad)
continue;
1860 if (mother->state[i].isFinal()) {
1861 double ppNow = mother->state[i].p() * pRad
1862 - mother->state[i].m() - mRad;
1863 if (ppNow < ppMin) {
1869 if (iRec)
return iRec;
1879 void History::findStateTransfer(map<int,int> &transfer) {
1881 if (!mother)
return;
1885 for(
int i = 0;i < 3; ++i)
1886 transfer.insert(pair<int,int>(i,i));
1888 transfer.insert(pair<int,int>(clusterIn.radBef, clusterIn.emittor));
1889 transfer.insert(pair<int,int>(clusterIn.recBef, clusterIn.recoiler));
1892 for (
int i = 0; i < int(mother->state.size()); ++i) {
1893 if (clusterIn.emitted == i ||
1894 clusterIn.emittor == i ||
1895 clusterIn.recoiler == i)
1898 for (
int j = 0;j < int(state.size()); ++j) {
1899 if (mother->state[i].id() == state[j].id()
1900 && mother->state[i].colType() == state[j].colType()
1901 && mother->state[i].chargeType() == state[j].chargeType()
1902 && mother->state[i].col() == state[j].col()
1903 && mother->state[i].acol() == state[j].acol()
1904 && mother->state[i].status() == state[j].status()) {
1905 transfer.insert(pair<int,int>(j,i));
1922 void History::findPath(vector<int>& out) {
1925 if (!mother &&
int(children.size()) < 1)
return;
1930 int size = int(mother->children.size());
1932 for (
int i=0; i < size; ++i) {
1933 if ( mother->children[i]->scale == scale
1934 && mother->children[i]->prob == prob
1935 && equalClustering(mother->children[i]->clusterIn,clusterIn)) {
1942 out.push_back(iChild);
1943 mother->findPath(out);
1968 void History::setScales( vector<int> index,
bool forward) {
1972 if ( children.empty() && forward ) {
1975 double scaleNew = 1.;
1976 if (mergingHooksPtr->incompleteScalePrescip()==0) {
1977 scaleNew = mergingHooksPtr->muF();
1978 }
else if (mergingHooksPtr->incompleteScalePrescip()==1) {
1980 pOut.p(0.,0.,0.,0.);
1981 for(
int i=0; i<int(state.size()); ++i)
1982 if (state[i].isFinal())
1983 pOut += state[i].p();
1984 scaleNew = pOut.mCalc();
1985 }
else if (mergingHooksPtr->incompleteScalePrescip()==2) {
1986 scaleNew = state[0].e();
1989 scaleNew = max( mergingHooksPtr->pTcut(), scaleNew);
1991 state.scale(scaleNew);
1992 for(
int i=3; i < int(state.size());++i)
1993 if (state[i].colType() != 0)
1994 state[i].scale(scaleNew);
1997 state.scale( state[0].e() );
1999 bool isLEP = ( state[3].isLepton() && state[4].isLepton() );
2001 int nFinalPartons = 0;
2002 int nFinalPhotons = 0;
2003 for (
int i=0; i < int(state.size()); ++i ) {
2004 if ( state[i].isFinal() ) {
2006 if ( state[i].colType() != 0 ) nFinalPartons++;
2007 if ( state[i].
id() == 22 ) nFinalPhotons++;
2010 bool isQCD = ( nFinal == 2 && nFinal == nFinalPartons );
2011 bool isPPh = ( nFinal == 2 && nFinalPartons == 1 && nFinalPhotons == 1);
2013 if ( !isLEP && ( isQCD || isPPh ) ) {
2014 double scaleNew = hardFacScale(state);
2015 state.scale( scaleNew );
2021 if (mother && forward) {
2023 double scaleNew = 1.;
2024 if (mergingHooksPtr->unorderedScalePrescip() == 0) {
2026 scaleNew = max( mergingHooksPtr->pTcut(), max(scale,mother->scale));
2027 }
else if (mergingHooksPtr->unorderedScalePrescip() == 1) {
2029 if (scale < mother->scale)
2030 scaleNew *= max( mergingHooksPtr->pTcut(), min(scale,mother->scale));
2032 scaleNew *= max( mergingHooksPtr->pTcut(), max(scale,mother->scale));
2037 mother->state[clusterIn.emitted].scale(scaleNew);
2038 mother->state[clusterIn.emittor].scale(scaleNew);
2039 mother->state[clusterIn.recoiler].scale(scaleNew);
2043 mother->scaleCopies(clusterIn.emitted, mother->state, scaleNew);
2044 mother->scaleCopies(clusterIn.emittor, mother->state, scaleNew);
2045 mother->scaleCopies(clusterIn.recoiler, mother->state, scaleNew);
2048 mother->setScales(index,
true);
2053 if (!mother || !forward) {
2056 if (
int(index.size()) > 0 ) {
2057 iChild = index.back();
2063 scale = max(mergingHooksPtr->pTcut(), scale);
2066 if (iChild != -1 && !children.empty()) {
2067 if (scale > children[iChild]->scale ) {
2068 if (mergingHooksPtr->unorderedScalePrescip() == 0) {
2070 double scaleNew = max( mergingHooksPtr->pTcut(),
2071 max(scale,children[iChild]->scale));
2073 for(
int i = 0; i < int(children[iChild]->state.size()); ++i)
2074 if (children[iChild]->state[i].scale() == children[iChild]->scale)
2075 children[iChild]->state[i].scale(scaleNew);
2077 children[iChild]->scale = scaleNew;
2079 }
else if ( mergingHooksPtr->unorderedScalePrescip() == 1) {
2081 double scaleNew = max(mergingHooksPtr->pTcut(),
2082 min(scale,children[iChild]->scale));
2084 for(
int i = 0; i < int(state.size()); ++i)
2085 if (state[i].scale() == scale)
2086 state[i].scale(scaleNew);
2092 double scalemin = state[0].e();
2093 for(
int i = 0; i < int(state.size()); ++i)
2094 if (state[i].colType() != 0)
2095 scalemin = max(mergingHooksPtr->pTcut(),
2096 min(scalemin,state[i].scale()));
2097 state.scale(scalemin);
2098 scale = max(mergingHooksPtr->pTcut(), scale);
2101 children[iChild]->setScales(index,
false);
2117 void History::scaleCopies(
int iPart,
const Event& refEvent,
double rho) {
2122 for(
int i=0; i < mother->state.size(); ++i) {
2123 if ( ( mother->state[i].id() == refEvent[iPart].id()
2124 && mother->state[i].colType() == refEvent[iPart].colType()
2125 && mother->state[i].chargeType() == refEvent[iPart].chargeType()
2126 && mother->state[i].col() == refEvent[iPart].col()
2127 && mother->state[i].acol() == refEvent[iPart].acol() )
2130 mother->state[i].scale(rho);
2133 mother->scaleCopies( iPart, refEvent, rho );
2146 void History::setEventScales() {
2150 mother->state.scale(scale);
2152 mother->setEventScales();
2162 double History::zISR() {
2165 if (!mother)
return 0.0;
2167 if (mother->state[clusterIn.emittor].isFinal())
return mother->zISR();
2169 int rad = clusterIn.emittor;
2170 int rec = clusterIn.recoiler;
2171 int emt = clusterIn.emitted;
2172 double z = (mother->state[rad].p() + mother->state[rec].p()
2173 - mother->state[emt].p()).m2Calc()
2174 / (mother->state[rad].p() + mother->state[rec].p()).m2Calc();
2176 double znew = mother->zISR();
2178 if (znew > 0.) z = znew;
2189 double History::zFSR() {
2192 if (!mother)
return 0.0;
2194 if (!mother->state[clusterIn.emittor].isFinal())
return mother->zFSR();
2196 int rad = clusterIn.emittor;
2197 int rec = clusterIn.recoiler;
2198 int emt = clusterIn.emitted;
2200 Vec4 sum = mother->state[rad].p() + mother->state[rec].p()
2201 + mother->state[emt].p();
2202 double m2Dip = sum.m2Calc();
2203 double x1 = 2. * (sum * mother->state[rad].p()) / m2Dip;
2204 double x3 = 2. * (sum * mother->state[emt].p()) / m2Dip;
2206 double z = x1/(x1+x3);
2208 double znew = mother->zFSR();
2210 if (znew > 0.) z = znew;
2221 double History::pTISR() {
2223 if (!mother)
return 0.0;
2225 if (mother->state[clusterIn.emittor].isFinal())
return mother->pTISR();
2226 double pT = mother->state.scale();
2228 double pTnew = mother->pTISR();
2230 if (pTnew > 0.) pT = pTnew;
2241 double History::pTFSR() {
2244 if (!mother)
return 0.0;
2246 if (!mother->state[clusterIn.emittor].isFinal())
return mother->pTFSR();
2247 double pT = mother->state.scale();
2249 double pTnew = mother->pTFSR();
2251 if (pTnew > 0.) pT = pTnew;
2262 int History::nClusterings() {
2263 if (!mother)
return 0;
2264 int w = mother->nClusterings();
2277 Event History::clusteredState(
int nSteps) {
2280 Event outState = state;
2282 if (mother && nSteps > 0)
2283 outState = mother->clusteredState(nSteps - 1);
2296 History * History::select(
double rnd) {
2299 if ( goodBranches.empty() && badBranches.empty() )
return this;
2303 map<double, History*> selectFrom;
2304 if ( !goodBranches.empty() ) {
2305 selectFrom = goodBranches;
2306 sum = sumGoodBranches;
2308 selectFrom = badBranches;
2309 sum = sumBadBranches;
2312 if (mergingHooksPtr->pickBySumPT()) {
2315 for (
int i=0; i < state.size(); ++i)
2316 if (state[i].isFinal())
2319 double sumMin = (nFinal-2)*state[0].e();
2320 for ( map<double, History*>::iterator it = selectFrom.begin();
2321 it != selectFrom.end(); ++it ) {
2323 if (it->second->sumScalarPT < sumMin) {
2324 sumMin = it->second->sumScalarPT;
2329 return selectFrom.lower_bound(iMin)->second;
2333 return selectFrom.upper_bound(sum*rnd)->second;
2335 return selectFrom.lower_bound(sum*rnd)->second;
2345 bool History::trimHistories() {
2347 if ( paths.empty() )
return false;
2349 for ( map<double, History*>::iterator it = paths.begin();
2350 it != paths.end(); ++it ) {
2352 if ( it->second->keep() && !it->second->keepHistory() )
2353 it->second->remove();
2356 double sumold, sumnew, sumprob, mismatch;
2357 sumold = sumnew = sumprob = mismatch = 0.;
2360 for ( map<double, History*>::iterator it = paths.begin();
2361 it != paths.end(); ++it ) {
2364 if ( it->second->keep() ) {
2366 goodBranches.insert( make_pair( sumnew - mismatch, it->second) );
2368 sumGoodBranches = sumnew - mismatch;
2372 double mismatchOld = mismatch;
2373 mismatch += sumnew - sumold;
2375 badBranches.insert( make_pair( mismatchOld + sumnew - sumold,
2378 sumBadBranches = mismatchOld + sumnew - sumold;
2386 return !goodBranches.empty();
2393 bool History::keepHistory() {
2394 bool keepPath =
true;
2397 if ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
2398 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
2399 || isQCD2to2(state) ) {
2402 double maxScale = hardFacScale(state);
2403 return keepPath = isOrderedPath( maxScale );
2407 if (isEW2to1(state)) {
2409 for (
int i = 0;i < state.size(); ++i)
2410 if (state[i].isFinal()) pSum += state[i].p();
2411 return isOrderedPath( pSum.mCalc());
2414 keepPath = isOrderedPath( infoPtr->eCM() );
2420 if (probMax() > 0. && abs(prob) < 1e-10*probMax()) keepPath=
false;
2431 bool History::isOrderedPath(
double maxscale ) {
2432 double newscale = clusterIn.pT();
2433 if ( !mother )
return true;
2434 if ( mother->state[clusterIn.emittor].idAbs() == 21
2435 && mother->state[clusterIn.emitted].idAbs() == 5
2436 && !mother->state[clusterIn.emittor].isFinal())
2438 bool ordered = mother->isOrderedPath(newscale);
2439 if ( !ordered || maxscale < newscale)
return false;
2448 bool History::allIntermediateAboveRhoMS(
double rhoms,
bool good ) {
2451 if ( !good )
return false;
2454 for (
int i = 0; i < state.size(); ++i )
2455 if ( state[i].isFinal() && state[i].colType() != 0 )
2457 double rhoNew = (nFinal > 0 ) ? mergingHooksPtr->tmsNow( state )
2460 if ( !mother )
return good;
2462 return mother->allIntermediateAboveRhoMS( rhoms, (rhoNew > rhoms) );
2469 bool History::foundAnyOrderedPaths() {
2471 if ( paths.empty() )
return false;
2472 double maxscale = infoPtr->eCM();
2474 for ( map<double, History*>::iterator it = paths.begin();
2475 it != paths.end(); ++it )
2476 if ( it->second->isOrderedPath(maxscale) )
2497 vector<double> History::weightTree(PartonLevel* trial,
double as0,
double aem0,
2498 double maxscale,
double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,
2499 AlphaEM * aemFSR, AlphaEM * aemISR, vector<double>& asWeight,
2500 vector<double>& aemWeight, vector<double>& pdfWeight) {
2503 double newScale = scale;
2505 int nWgts = mergingHooksPtr->nWgts;
2510 int sideRad = (state[3].pz() > 0) ? 1 :-1;
2511 int sideRec = (state[4].pz() > 0) ? 1 :-1;
2514 if (state[3].colType() != 0) {
2516 double x = 2.*state[3].e() / state[0].e();
2517 int flav = state[3].id();
2519 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2520 double scaleDen = mergingHooksPtr->muFinME();
2522 double ratio = getPDFratio(sideRad,
false,
false, flav, x, scaleNum,
2524 for (
double& pdfW: pdfWeight) pdfW *= ratio;
2528 if (state[4].colType() != 0) {
2530 double x = 2.*state[4].e() / state[0].e();
2531 int flav = state[4].id();
2533 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2534 double scaleDen = mergingHooksPtr->muFinME();
2536 double ratio = getPDFratio(sideRec,
false,
false, flav, x, scaleNum,
2538 for (
double& pdfW: pdfWeight) pdfW *= ratio;
2540 return vector<double>(nWgts,1.);
2545 double newPDFscale = newScale;
2546 if (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2547 newPDFscale = clusterIn.pT();
2550 vector<double> w = mother->weightTree(trial, as0, aem0, newScale,
2551 newPDFscale, asFSR, asISR, aemFSR, aemISR, asWeight, aemWeight,
2555 if (state.size() < 3)
return vector<double>(nWgts,1.);
2557 if ( w[0] < 1e-12 )
return vector<double>(nWgts,0.);
2559 vector<double> wTrial = doTrialShower(trial, 1, maxscale);
2560 for (
int iVar = 0; iVar < nWgts; ++iVar)
2561 w[iVar] *= wTrial[iVar];
2562 if ( w[0] < 1e-12 )
return vector<double>( nWgts, 0. );
2564 int emtType = mother->state[clusterIn.emitted].colType();
2566 if ( asFSR && asISR && emtType != 0) {
2567 double asScale = pow2( newScale );
2568 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2569 asScale = pow2( clusterIn.pT() );
2572 bool FSR = mother->state[clusterIn.emittor].isFinal();
2573 if (!FSR) asScale += pow2(mergingHooksPtr->pT0ISR());
2576 if (mergingHooksPtr->useShowerPlugin() )
2577 asScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2578 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale);
2580 double alphaSinPS = (FSR) ? (*asFSR).alphaS(asScale)
2581 : (*asISR).alphaS(asScale);
2582 asWeight[0] *= alphaSinPS / as0;
2585 for (
int iVar = 1; iVar < nWgts; ++iVar) {
2586 alphaSinPS = (FSR) ?
2587 asFSR->alphaS(pow2(mergingHooksPtr->muRVarFactors[iVar-1])*asScale) :
2588 asISR->alphaS(pow2(mergingHooksPtr->muRVarFactors[iVar-1])*asScale);
2589 asWeight[iVar] *= alphaSinPS / as0;
2594 if ( aemFSR && aemISR && emtType == 0 ) {
2595 double aemScale = pow2( newScale );
2596 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2597 aemScale = pow2( clusterIn.pT() );
2600 bool FSR = mother->state[clusterIn.emittor].isFinal();
2601 if (!FSR) aemScale += pow2(mergingHooksPtr->pT0ISR());
2604 if (mergingHooksPtr->useShowerPlugin() )
2605 aemScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2606 clusterIn.emitted, clusterIn.recoiler,
"scaleEM", aemScale);
2608 double alphaEMinPS = (FSR) ? (*aemFSR).alphaEM(aemScale)
2609 : (*aemISR).alphaEM(aemScale);
2610 for (
double& aemW: aemWeight) aemW *= alphaEMinPS / aem0;
2616 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
2617 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
2619 if ( mother->state[inP].colType() != 0 ) {
2621 double x = getCurrentX(sideP);
2622 int flav = getCurrentFlav(sideP);
2624 double scaleNum = (children.empty())
2625 ? hardFacScale(state)
2626 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2627 ? pdfScale : maxscale );
2628 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2629 ? clusterIn.pT() : newScale;
2631 double ratio = getPDFratio(sideP,
false,
false, flav, x, scaleNum,
2633 for (
double& pdfW: pdfWeight) pdfW *= ratio;
2636 if ( mother->state[inM].colType() != 0 ) {
2638 double x = getCurrentX(sideM);
2639 int flav = getCurrentFlav(sideM);
2641 double scaleNum = (children.empty())
2642 ? hardFacScale(state)
2643 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2644 ? pdfScale : maxscale );
2645 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2646 ? clusterIn.pT() : newScale;
2648 double ratio = getPDFratio(sideM,
false,
false, flav, x, scaleNum,
2650 for (
double& pdfW: pdfWeight) pdfW *= ratio;
2661 vector<double> History::weightTreeAlphaS(
double as0, AlphaStrong * asFSR,
2662 AlphaStrong * asISR,
int njetMax,
bool asVarInME ) {
2664 int nWgts = mergingHooksPtr->nWgts;
2667 if ( !mother )
return vector<double>(nWgts,1.);
2669 vector<double> w = mother->weightTreeAlphaS( as0, asFSR, asISR, njetMax,
2672 if (state.size() < 3)
return w;
2675 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2676 if (njetNow >= njetMax)
return vector<double>( nWgts, 1.0 );
2679 bool FSR = mother->state[clusterIn.emittor].isFinal();
2680 int emtID = mother->state[clusterIn.emitted].id();
2683 if (abs(emtID) == 22 || abs(emtID) == 23 || abs(emtID) == 24)
return w;
2686 if ( asFSR && asISR ) {
2687 double asScale = pow2( scale );
2688 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2689 asScale = pow2( clusterIn.pT() );
2692 if (!FSR) asScale += pow2(mergingHooksPtr->pT0ISR());
2695 if (mergingHooksPtr->useShowerPlugin() )
2696 asScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2697 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale);
2699 double alphaSinPS = (FSR) ? (*asFSR).alphaS(asScale)
2700 : (*asISR).alphaS(asScale);
2701 w[0] *= alphaSinPS / as0;
2704 for (
int iVar = 1; iVar < nWgts; ++iVar) {
2705 alphaSinPS = (FSR) ?
2706 asFSR->alphaS(pow2(mergingHooksPtr->muRVarFactors[iVar-1])*asScale) :
2707 asISR->alphaS(pow2(mergingHooksPtr->muRVarFactors[iVar-1])*asScale);
2709 double muR2 = pow2(mergingHooksPtr->muRinMESave);
2710 double alphaSinME = (!asVarInME) ? as0 : ((FSR) ?
2711 asFSR->alphaS(muR2*pow2(mergingHooksPtr->muRVarFactors[iVar-1])) :
2712 asISR->alphaS(muR2*pow2(mergingHooksPtr->muRVarFactors[iVar-1])));
2713 w[iVar] *= alphaSinPS / alphaSinME;
2725 vector<double> History::weightTreeAlphaEM(
double aem0, AlphaEM * aemFSR,
2726 AlphaEM * aemISR,
int njetMax ) {
2728 int nWgts = mergingHooksPtr->nWgts;
2731 if ( !mother )
return vector<double>( nWgts, 1. );
2733 vector<double> w = mother->weightTreeAlphaEM( aem0, aemFSR, aemISR,
2736 if (state.size() < 3)
return w;
2739 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2740 if (njetNow >= njetMax)
return vector<double>( nWgts, 1.0 );
2743 bool FSR = mother->state[clusterIn.emittor].isFinal();
2744 int emtID = mother->state[clusterIn.emitted].id();
2747 if (!(abs(emtID) == 22 || abs(emtID) == 23 || abs(emtID) == 24))
return w;
2750 if ( aemFSR && aemISR ) {
2751 double aemScale = pow2( scale );
2752 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2753 aemScale = pow2( clusterIn.pT() );
2756 if (!FSR) aemScale += pow2(mergingHooksPtr->pT0ISR());
2759 if (mergingHooksPtr->useShowerPlugin() )
2760 aemScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2761 clusterIn.emitted, clusterIn.recoiler,
"scaleEM", aemScale);
2763 double alphaEMinPS = (FSR) ? (*aemFSR).alphaEM(aemScale)
2764 : (*aemISR).alphaEM(aemScale);
2765 for (
double& wi: w) wi *= alphaEMinPS / aem0;
2776 vector<double> History::weightTreePDFs(
double maxscale,
double pdfScale,
2780 double newScale = scale;
2782 double nWgts = mergingHooksPtr->nWgts;
2788 int njet = mergingHooksPtr->getNumberOfClusteringSteps( state);
2789 if (njet > njetMax)
return vector<double>( nWgts, 1.0 );
2791 vector<double> wt( nWgts, 1. );
2792 int sideRad = (state[3].pz() > 0) ? 1 :-1;
2793 int sideRec = (state[4].pz() > 0) ? 1 :-1;
2796 if (state[3].colType() != 0) {
2798 double x = 2.*state[3].e() / state[0].e();
2799 int flav = state[3].id();
2801 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2802 double scaleDen = mergingHooksPtr->muFinME();
2804 double ratio = getPDFratio(sideRad,
false,
false, flav, x, scaleNum,
2806 for (
double& w: wt) w *= ratio;
2810 if (state[4].colType() != 0) {
2812 double x = 2.*state[4].e() / state[0].e();
2813 int flav = state[4].id();
2815 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2816 double scaleDen = mergingHooksPtr->muFinME();
2818 double ratio = getPDFratio(sideRec,
false,
false, flav, x, scaleNum,
2820 for (
double& w: wt) w *= ratio;
2828 double newPDFscale = newScale;
2829 if ( mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2830 newPDFscale = clusterIn.pT();
2833 vector<double> w = mother->weightTreePDFs( newScale, newPDFscale, njetMax );
2836 if (state.size() < 3)
return w;
2839 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2840 if (njetNow >= njetMax)
return vector<double>( nWgts, 1. );
2845 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
2846 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
2848 if ( mother->state[inP].colType() != 0 ) {
2850 double x = getCurrentX(sideP);
2851 int flav = getCurrentFlav(sideP);
2853 double scaleNum = (children.empty())
2854 ? hardFacScale(state)
2855 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2856 ? pdfScale : maxscale );
2857 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2858 ? clusterIn.pT() : newScale;
2860 double xDen = (njetNow == njetMax) ? mother->getCurrentX(sideP) : x;
2861 int flavDen = (njetNow == njetMax) ? mother->getCurrentFlav(sideP) : flav;
2862 double sDen = (njetNow == njetMax) ? mergingHooksPtr->muFinME() : scaleDen;
2863 double ratio = getPDFratio(sideP,
false,
false, flav, x, scaleNum,
2864 flavDen, xDen, sDen);
2865 for (
double& wi: w) wi *= ratio;
2868 if ( mother->state[inM].colType() != 0 ) {
2870 double x = getCurrentX(sideM);
2871 int flav = getCurrentFlav(sideM);
2873 double scaleNum = (children.empty())
2874 ? hardFacScale(state)
2875 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2876 ? pdfScale : maxscale );
2877 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2878 ? clusterIn.pT() : newScale;
2880 double xDen = (njetNow == njetMax) ? mother->getCurrentX(sideM) : x;
2881 int flavDen = (njetNow == njetMax) ? mother->getCurrentFlav(sideM) : flav;
2882 double sDen = (njetNow == njetMax) ? mergingHooksPtr->muFinME() : scaleDen;
2883 double ratio = getPDFratio(sideM,
false,
false, flav, x, scaleNum,
2884 flavDen, xDen, sDen);
2885 for (
double& wi: w) wi *= ratio;
2896 vector<double> History::weightTreeEmissions( PartonLevel* trial,
int type,
2897 int njetMin,
int njetMax,
double maxscale ) {
2899 int nWgts = mergingHooksPtr->nWgts;
2901 if (type == -1 && !mergingHooksPtr->settingsPtr->flag(
"PartonLevel:MPI"))
2902 return vector<double>( nWgts, 1. );
2905 double newScale = scale;
2908 if ( !mother )
return vector<double>( nWgts, 1.0 );
2910 vector<double> w = mother->weightTreeEmissions( trial, type, njetMin,
2911 njetMax, newScale );
2913 if (state.size() < 3)
return vector<double>( nWgts, 1.0 );
2915 if ( w[0] < 1e-12 )
return vector<double>( nWgts, 0.0 );
2917 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2918 if (njetNow >= njetMax)
return vector<double>( nWgts, 1.0);
2921 vector<double> wTrial = doTrialShower( trial, type, maxscale );
2922 for (
int iVar = 0; iVar < nWgts; ++iVar) w[iVar] *= wTrial[iVar];
2924 if ( w[0] < 1e-12 )
return vector<double>( nWgts, 0.0 );
2934 double History::weightFirst(PartonLevel* trial,
double as0,
double muR,
2935 double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR, Rndm* rndmPtr ) {
2938 double newScale = scale;
2945 if (state[3].colType() != 0) {
2947 double x = 2.*state[3].e() / state[0].e();
2948 int flav = state[3].id();
2950 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2951 double scaleDen = mergingHooksPtr->muFinME();
2953 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2954 mergingHooksPtr->muFinME(), as0, rndmPtr);
2959 if (state[4].colType() != 0) {
2961 double x = 2.*state[4].e() / state[0].e();
2962 int flav = state[4].id();
2964 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2965 double scaleDen = mergingHooksPtr->muFinME();
2967 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2968 mergingHooksPtr->muFinME(), as0, rndmPtr);
2976 double w = mother->weightFirst(trial, as0, muR, newScale, asFSR, asISR,
2980 if (state.size() < 3)
return 0.0;
2984 double asScale2 = newScale*newScale;
2985 int showerType = (mother->state[clusterIn.emittor].isFinal() ) ? 1 : -1;
2986 if (showerType == -1) {
2987 asScale2 += pow(mergingHooksPtr->pT0ISR(),2);
2992 if (mergingHooksPtr->useShowerPlugin() ){
2993 asScale2 = getShowerPluginScale(mother->state, clusterIn.emittor,
2994 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale2);
3000 double BETA0 = 11. - 2./3.* NF;
3002 w += as0 / (2.*M_PI) * 0.5 * BETA0 * log( (muR*muR) / (b*asScale2) );
3008 double nWeight1 = 0.;
3009 double nWeight2 = 0.;
3011 for(
int i=0; i < NTRIAL; ++i) {
3013 vector<double> unresolvedEmissionTerm = countEmissions(trial, maxscale,
3014 newScale, 2, as0, asFSR, asISR, 3, fixpdf, fixas);
3015 nWeight1 += unresolvedEmissionTerm[1];
3017 w += nWeight1/double(NTRIAL) + nWeight2/double(NTRIAL);
3022 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
3023 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
3025 if ( mother->state[inP].colType() != 0 ) {
3027 double x = getCurrentX(sideP);
3028 int flav = getCurrentFlav(sideP);
3030 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
3032 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, newScale,
3033 mergingHooksPtr->muFinME(), as0, rndmPtr);
3038 if ( mother->state[inM].colType() != 0 ) {
3040 double x = getCurrentX(sideM);
3041 int flav = getCurrentFlav(sideM);
3043 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
3045 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, newScale,
3046 mergingHooksPtr->muFinME(), as0, rndmPtr);
3061 double History::weightFirstAlphaS(
double as0,
double muR,
3062 AlphaStrong * asFSR, AlphaStrong * asISR ) {
3065 double newScale = scale;
3067 if ( !mother )
return 0.;
3069 double w = mother->weightFirstAlphaS( as0, muR, asFSR, asISR );
3071 int showerType = (mother->state[clusterIn.emittor].isFinal() ) ? 1 : -1;
3073 double asScale = pow2( newScale );
3074 if ( mergingHooksPtr->unorderedASscalePrescip() == 1 )
3075 asScale = pow2( clusterIn.pT() );
3076 if (showerType == -1) {
3077 asScale += pow2( mergingHooksPtr->pT0ISR() );
3082 if (mergingHooksPtr->useShowerPlugin() ) {
3083 asScale = getShowerPluginScale(mother->state, clusterIn.emittor,
3084 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale);
3090 double BETA0 = 11. - 2./3.* NF;
3092 w += as0 / (2.*M_PI) * 0.5 * BETA0 * log( (muR*muR) / (b*asScale) );
3104 double History::weightFirstPDFs(
double as0,
double maxscale,
double pdfScale,
3108 double newScale = scale;
3115 if (state[3].colType() != 0) {
3117 double x = 2.*state[3].e() / state[0].e();
3118 int flav = state[3].id();
3120 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
3121 double scaleDen = mergingHooksPtr->muFinME();
3123 wt += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
3124 mergingHooksPtr->muFinME(), as0, rndmPtr);
3127 if (state[4].colType() != 0) {
3129 double x = 2.*state[4].e() / state[0].e();
3130 int flav = state[4].id();
3132 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
3133 double scaleDen = mergingHooksPtr->muFinME();
3135 wt += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
3136 mergingHooksPtr->muFinME(), as0, rndmPtr);
3145 double newPDFscale = newScale;
3146 if (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
3147 newPDFscale = clusterIn.pT();
3150 double w = mother->weightFirstPDFs( as0, newScale, newPDFscale, rndmPtr);
3155 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
3156 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
3158 if ( mother->state[inP].colType() != 0 ) {
3160 double x = getCurrentX(sideP);
3161 int flav = getCurrentFlav(sideP);
3163 double scaleNum = (children.empty())
3164 ? hardFacScale(state)
3165 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
3166 ? pdfScale : maxscale );
3167 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
3168 ? clusterIn.pT() : newScale;
3170 w += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
3171 mergingHooksPtr->muFinME(), as0, rndmPtr);
3174 if ( mother->state[inM].colType() != 0 ) {
3176 double x = getCurrentX(sideM);
3177 int flav = getCurrentFlav(sideM);
3179 double scaleNum = (children.empty())
3180 ? hardFacScale(state)
3181 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
3182 ? pdfScale : maxscale );
3183 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
3184 ? clusterIn.pT() : newScale;
3186 w += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
3187 mergingHooksPtr->muFinME(), as0, rndmPtr);
3201 double History::weightFirstEmissions(PartonLevel* trial,
double as0,
3202 double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR,
3203 bool fixpdf,
bool fixas ) {
3206 double newScale = scale;
3207 if ( !mother )
return 0.0;
3209 double w = mother->weightFirstEmissions(trial, as0, newScale, asFSR, asISR,
3212 if (state.size() < 3)
return 0.0;
3214 double nWeight1 = 0.;
3215 double nWeight2 = 0.;
3216 for(
int i=0; i < NTRIAL; ++i) {
3218 vector<double> unresolvedEmissionTerm = countEmissions(trial, maxscale,
3219 newScale, 2, as0, asFSR, asISR, 3, fixpdf, fixas);
3220 nWeight1 += unresolvedEmissionTerm[1];
3223 w += nWeight1/double(NTRIAL) + nWeight2/double(NTRIAL);
3234 double History::hardFacScale(
const Event& event) {
3236 double hardscale = 0.;
3238 if ( !mergingHooksPtr->resetHardQFac() )
return mergingHooksPtr->muF();
3242 if ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
3243 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
3244 || isQCD2to2(event)) {
3247 for (
int i=0; i <
event.size(); ++i)
3248 if ( event[i].isFinal() &&
event[i].colType() != 0 )
3249 mT.push_back( abs(event[i].mT2()) );
3250 if (
int(mT.size()) != 2 )
3251 hardscale = infoPtr->QFac();
3253 hardscale = sqrt( min( mT[0], mT[1] ) );
3255 hardscale = mergingHooksPtr->muF();
3265 double History::hardRenScale(
const Event& event) {
3267 double hardscale = 0.;
3269 if ( !mergingHooksPtr->resetHardQRen() )
return mergingHooksPtr->muR();
3273 if ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
3274 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
3275 || isQCD2to2(event)) {
3278 for (
int i=0; i <
event.size(); ++i)
3279 if ( event[i].isFinal()
3280 && (
event[i].colType() != 0 ||
event[i].id() == 22 ) )
3281 mT.push_back( abs(event[i].mT()) );
3282 if (
int(mT.size()) != 2 )
3283 hardscale = infoPtr->QRen();
3285 hardscale = sqrt( mT[0]*mT[1] );
3287 hardscale = mergingHooksPtr->muR();
3304 vector<double> History::doTrialShower( PartonLevel* trial,
int type,
3305 double maxscaleIn,
double minscaleIn ) {
3308 Event process = state;
3310 double startingScale = maxscaleIn;
3313 if ( mergingHooksPtr->getNumberOfClusteringSteps(process) == 0
3314 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
3315 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
3316 || isQCD2to2(state) ) )
3317 startingScale = min( startingScale, hardFacScale(process) );
3319 int nWgts = mergingHooksPtr->nWgts;
3322 bool doVeto =
false;
3324 bool canEnhanceTrial = trial->canEnhanceTrial();
3327 vector<double> showerWeightVecSave = infoPtr->
3328 weightContainerPtr->weightsPS.weightValues;
3333 trial->resetTrial();
3336 event.init(
"(hard process-modified)", particleDataPtr);
3340 process.scale(startingScale);
3344 double minScale = (minscaleIn > 0.) ? minscaleIn : scale;
3349 if (minScale >= startingScale)
break;
3355 double z = ( mergingHooksPtr->getNumberOfClusteringSteps(state) == 0
3358 : mother->getCurrentZ(clusterIn.emittor,clusterIn.recoiler,
3359 clusterIn.emitted, clusterIn.flavRadBef);
3361 infoPtr->zNowISR(z);
3362 infoPtr->pT2NowISR(pow(startingScale,2));
3363 infoPtr->hasHistory(
true);
3366 if (mergingHooksPtr->doWeakClustering()) setupSimpleWeakShower(0);
3369 mergingHooksPtr->setShowerStoppingScale(minScale);
3371 trial->next(process,event);
3373 double pTtrial = trial->pTLastInShower();
3374 int typeTrial = trial->typeLastInShower();
3377 trial->resetTrial();
3380 double pTEnhanced = trial->getEnhancedTrialPT();
3381 double wtEnhanced = trial->getEnhancedTrialWeight();
3382 if ( canEnhanceTrial && pTEnhanced > 0.) pTtrial = pTEnhanced;
3385 double vetoScale = (mother) ? 0. : mergingHooksPtr->tms();
3387 double tnow = mergingHooksPtr->tmsNow( event );
3390 if ( pTtrial < minScale )
break;
3392 startingScale = pTtrial;
3395 if ( tnow < vetoScale && vetoScale > 0. )
continue;
3398 if ( mergingHooksPtr->canVetoTrialEmission()
3399 && mergingHooksPtr->doVetoTrialEmission( process, event) )
continue;
3401 int iRecAft =
event.size() - 1;
3402 int iEmt =
event.size() - 2;
3403 int iRadAft =
event.size() - 3;
3404 if ( (event[iRecAft].status() != 52 && event[iRecAft].status() != -53) ||
3405 event[iEmt].status() != 51 || event[iRadAft].status() != 51)
3406 iRecAft = iEmt = iRadAft = -1;
3407 for (
int i = event.size() - 1; i > 0; i--) {
3408 if (iRadAft == -1 && event[i].status() == -41) iRadAft = i;
3409 else if (iEmt == -1 && event[i].status() == 43) iEmt = i;
3410 else if (iRecAft == -1 && event[i].status() == -42) iRecAft = i;
3411 if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
break;
3416 if ( type == -1 && typeTrial != 1 ) {
3418 infoPtr->weightContainerPtr->weightsPS.weightValues =
3419 showerWeightVecSave;
3423 if ( type == 1 && !(typeTrial == 2 || typeTrial >= 3) ) {
3425 infoPtr->weightContainerPtr->weightsPS.weightValues =
3426 showerWeightVecSave;
3431 if (canEnhanceTrial && pTtrial > minScale) wt *= (1. - 1./wtEnhanced);
3433 if ( canEnhanceTrial && wt == 0.)
break;
3435 if ( canEnhanceTrial && pTtrial > minScale)
continue;
3439 if ( pTtrial > minScale ) doVeto =
true;
3449 && mergingHooksPtr->getNumberOfClusteringSteps(process) == 0
3450 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
3451 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
3452 || isQCD2to2(state))
3453 && pTtrial > hardFacScale(process) )
3454 return vector<double>( nWgts, 0.0 );
3458 if ( pTtrial < minScale ) doVeto =
false;
3467 vector<double> trialShowerWt = infoPtr->weightContainerPtr->weightsPS.
3468 getMuRWeightVector();
3470 trialShowerWt.insert(trialShowerWt.begin(),1.);
3474 if (type == -1) trialShowerWt = vector<double>(trialShowerWt.size(),1.);
3476 infoPtr->weightContainerPtr->weightsPS.weightValues = showerWeightVecSave;
3479 for (
double& swt: trialShowerWt)
3480 swt *= (canEnhanceTrial) ? wt : ( (doVeto) ? 0. : 1. );
3481 return trialShowerWt;
3492 bool History::updateind(vector<int> & ind,
int i,
int N) {
3493 if ( i < 0 )
return false;
3494 if ( ++ind[i] < N )
return true;
3495 if ( !updateind(ind, i - 1, N - 1) )
return false;
3496 ind[i] = ind[i - 1] + 1;
3506 vector<double> History::countEmissions(PartonLevel* trial,
double maxscale,
3507 double minscale,
int showerType,
double as0,
3508 AlphaStrong * asFSR, AlphaStrong * asISR,
int N = 1,
3509 bool fixpdf =
true,
bool fixas =
true) {
3511 if ( N < 0 )
return vector<double>();
3512 vector<double> result(N+1);
3514 if ( N < 1 )
return result;
3517 Event process = state;
3519 double startingScale = maxscale;
3522 if ( mergingHooksPtr->getNumberOfClusteringSteps(process) == 0
3523 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
3524 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
3525 || isQCD2to2(state) ) )
3526 startingScale = min( startingScale, hardFacScale(process) );
3529 bool canEnhanceTrial = trial->canEnhanceTrial();
3532 vector<double> showerWeightVecSave = infoPtr->
3533 weightContainerPtr->weightsPS.weightValues;
3537 trial->resetTrial();
3540 event.init(
"(hard process-modified)", particleDataPtr);
3544 process.scale(startingScale);
3549 if (minscale >= startingScale)
return result;
3555 double z = ( mergingHooksPtr->getNumberOfClusteringSteps(state) == 0)
3557 : mother->getCurrentZ(clusterIn.emittor,clusterIn.recoiler,
3560 infoPtr->zNowISR(z);
3561 infoPtr->pT2NowISR(pow(startingScale,2));
3562 infoPtr->hasHistory(
true);
3566 if (mergingHooksPtr->doWeakClustering()) setupSimpleWeakShower(0);
3569 trial->next(process,event);
3572 infoPtr->weightContainerPtr->weightsPS.weightValues = showerWeightVecSave;
3575 double pTtrial = trial->pTLastInShower();
3576 int typeTrial = trial->typeLastInShower();
3579 trial->resetTrial();
3582 double pTEnhanced = trial->getEnhancedTrialPT();
3583 double wtEnhanced = trial->getEnhancedTrialWeight();
3584 if ( canEnhanceTrial && pTEnhanced > 0.) pTtrial = pTEnhanced;
3587 double vetoScale = (mother) ? 0. : mergingHooksPtr->tms();
3589 double tnow = mergingHooksPtr->tmsNow( event );
3592 startingScale = pTtrial;
3594 if ( pTtrial < minscale )
break;
3596 if ( tnow < vetoScale && vetoScale > 0. )
continue;
3598 if ( mergingHooksPtr->canVetoTrialEmission()
3599 && mergingHooksPtr->doVetoTrialEmission( process, event) )
continue;
3602 double enhance = (canEnhanceTrial && pTtrial > minscale) ? wtEnhanced : 1.;
3607 double alphaSinPS = as0;
3610 double asScale2 = pTtrial*pTtrial;
3612 if (mergingHooksPtr->useShowerPlugin() )
3613 asScale2 = getShowerPluginScale(mother->state, clusterIn.emittor,
3614 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale2);
3617 if ( (showerType == -1 || showerType == 2) && typeTrial == 2 ) {
3619 if ( fixas ) alphaSinPS = (*asISR).alphaS(asScale2);
3622 pdfs = pdfFactor( event, typeTrial, pTtrial,
3623 mergingHooksPtr->muFinME() );
3625 }
else if ( (showerType == 1 || showerType == 2) && typeTrial >= 3 ) {
3627 if ( fixas ) alphaSinPS = (*asFSR).alphaS(asScale2);
3631 pdfs = pdfFactor( event, typeTrial, pTtrial,
3632 mergingHooksPtr->muFinME() );
3636 if ( typeTrial == 2 || typeTrial >= 3 )
3637 wts.push_back(as0/alphaSinPS * pdfs * 1./enhance);
3641 for (
int n = 1; n <= min(N,
int(wts.size())); ++n ) {
3643 for (
int i = 0; i < N; ++i ) ind[i] = i;
3646 for (
int j = 0; j < n; ++j ) x *= wts[ind[j]];
3648 }
while ( updateind(ind, n - 1, wts.size()) );
3649 if ( n%2 ) result[n] *= -1.0;
3661 double History::monteCarloPDFratios(
int flav,
double x,
double maxScale,
3662 double minScale,
double pdfScale,
double asME, Rndm* rndmPtr) {
3666 double factor = asME / (2.*M_PI);
3668 factor *= log(maxScale/minScale);
3671 if (factor == 0.)
return 0.;
3679 double integral = 0.;
3680 double RN = rndmPtr->flat();
3683 double zTrial = pow(x,RN);
3684 integral = -log(x) * zTrial *
3685 integrand(flav, x, pdfScale, zTrial);
3686 integral += 1./6.*(11.*CA - 4.*NF*TR)
3689 double zTrial = x + RN*(1. - x);
3691 integrand(flav, x, pdfScale, zTrial);
3692 integral += 3./2.*CF
3697 return (factor*integral);
3706 bool History::onlyOrderedPaths() {
3707 if ( !mother || foundOrderedPath )
return foundOrderedPath;
3708 return foundOrderedPath = mother->onlyOrderedPaths();
3717 bool History::onlyStronglyOrderedPaths() {
3718 if ( !mother || foundStronglyOrderedPath )
return foundStronglyOrderedPath;
3719 return foundStronglyOrderedPath = mother->onlyStronglyOrderedPaths();
3728 bool History::onlyAllowedPaths() {
3729 if ( !mother || foundAllowedPath )
return foundAllowedPath;
3730 return foundAllowedPath = mother->onlyAllowedPaths();
3742 bool History::registerPath(History & l,
bool isOrdered,
3743 bool isStronglyOrdered,
bool isAllowed,
bool isComplete) {
3749 if ( mother )
return mother->registerPath(l, isOrdered,
3750 isStronglyOrdered, isAllowed, isComplete);
3753 if ( sumpath == sumpath + l.prob )
3755 if ( mergingHooksPtr->canCutOnRecState()
3756 && foundAllowedPath && !isAllowed )
3758 if ( mergingHooksPtr->enforceStrongOrdering()
3759 && foundStronglyOrderedPath && !isStronglyOrdered )
3761 if ( mergingHooksPtr->orderHistories()
3762 && foundOrderedPath && !isOrdered ) {
3764 if ( (!foundCompletePath && isComplete)
3765 || (!foundAllowedPath && isAllowed) ) ;
3769 if ( foundCompletePath && !isComplete)
3771 if ( !mergingHooksPtr->canCutOnRecState()
3772 && !mergingHooksPtr->allowCutOnRecState() )
3773 foundAllowedPath =
true;
3775 if ( mergingHooksPtr->canCutOnRecState() && isAllowed && isComplete) {
3776 if ( !foundAllowedPath || !foundCompletePath ) {
3782 foundAllowedPath =
true;
3786 if ( mergingHooksPtr->enforceStrongOrdering() && isStronglyOrdered
3788 if ( !foundStronglyOrderedPath || !foundCompletePath ) {
3794 foundStronglyOrderedPath =
true;
3795 foundCompletePath =
true;
3799 if ( mergingHooksPtr->orderHistories() && isOrdered && isComplete ) {
3800 if ( !foundOrderedPath || !foundCompletePath ) {
3806 foundOrderedPath =
true;
3807 foundCompletePath =
true;
3812 if ( !foundCompletePath ) {
3818 foundCompletePath =
true;
3823 foundOrderedPath =
true;
3827 double weakProb = 1.;
3828 if (mergingHooksPtr->doWeakClustering())
3829 weakProb = l.getWeakProb();
3832 sumpath += l.prob * weakProb;
3833 paths[sumpath] = &l;
3835 updateProbMax(l.prob * weakProb, isComplete);
3847 vector<Clustering> History::getAllQCDClusterings() {
3848 vector<Clustering> ret;
3851 vector <int> posFinalPartn;
3852 vector <int> posInitPartn;
3853 vector <int> posFinalGluon;
3854 vector <int> posFinalQuark;
3855 vector <int> posFinalAntiq;
3856 vector <int> posInitGluon;
3857 vector <int> posInitQuark;
3858 vector <int> posInitAntiq;
3862 for (
int i=0; i < state.size(); ++i )
3863 if ( state[i].isFinal() && state[i].colType() !=0 ) {
3865 if ( state[i].
id() == 21 ) posFinalGluon.push_back(i);
3866 else if ( state[i].idAbs() < 10 && state[i].
id() > 0)
3867 posFinalQuark.push_back(i);
3868 else if ( state[i].idAbs() < 10 && state[i].
id() < 0)
3869 posFinalAntiq.push_back(i);
3870 }
else if (state[i].status() == -21 && state[i].colType() != 0 ) {
3872 if ( state[i].
id() == 21 ) posInitGluon.push_back(i);
3873 else if ( state[i].idAbs() < 10 && state[i].
id() > 0)
3874 posInitQuark.push_back(i);
3875 else if ( state[i].idAbs() < 10 && state[i].
id() < 0)
3876 posInitAntiq.push_back(i);
3880 vector<Clustering> systems;
3881 systems = getQCDClusterings(state);
3882 ret.insert(ret.end(), systems.begin(), systems.end());
3886 if ( !ret.empty() )
return ret;
3891 else if ( mergingHooksPtr->allowColourShuffling() ) {
3894 for(
int i = 0; i < int(posFinalQuark.size()); ++i) {
3896 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(posFinalQuark[i],
3899 int col = NewState[posFinalQuark[i]].col();
3900 for(
int j = 0; j < int(posInitAntiq.size()); ++j) {
3902 int acl = NewState[posInitAntiq[j]].acol();
3903 if ( col == acl )
continue;
3904 NewState[posFinalQuark[i]].col(acl);
3905 NewState[posInitAntiq[j]].acol(col);
3906 systems = getQCDClusterings(NewState);
3907 if (!systems.empty()) {
3910 ret.insert(ret.end(), systems.begin(), systems.end());
3917 for(
int i = 0; i < int(posFinalAntiq.size()); ++i) {
3919 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(posFinalAntiq[i],
3922 int acl = NewState[posFinalAntiq[i]].acol();
3923 for(
int j = 0; j < int(posInitQuark.size()); ++j) {
3925 int col = NewState[posInitQuark[j]].col();
3926 if ( col == acl )
continue;
3927 NewState[posFinalAntiq[i]].acol(col);
3928 NewState[posInitQuark[j]].col(acl);
3929 systems = getQCDClusterings(NewState);
3930 if (!systems.empty()) {
3933 ret.insert(ret.end(), systems.begin(), systems.end());
3952 vector<Clustering> History::getQCDClusterings(
const Event& event) {
3953 vector<Clustering> ret;
3957 vector <int> posFinalPartn;
3958 vector <int> posInitPartn;
3960 vector <int> posFinalGluon;
3961 vector <int> posFinalQuark;
3962 vector <int> posFinalAntiq;
3963 vector <int> posInitGluon;
3964 vector <int> posInitQuark;
3965 vector <int> posInitAntiq;
3969 for (
int i=0; i <
event.size(); ++i)
3970 if ( event[i].isFinal() &&
event[i].colType() !=0 ) {
3972 posFinalPartn.push_back(i);
3973 if ( event[i].
id() == 21 ) posFinalGluon.push_back(i);
3974 else if ( event[i].idAbs() < 10 && event[i].
id() > 0)
3975 posFinalQuark.push_back(i);
3976 else if ( event[i].idAbs() < 10 && event[i].
id() < 0)
3977 posFinalAntiq.push_back(i);
3978 }
else if ( event[i].status() == -21 && event[i].colType() != 0 ) {
3980 posInitPartn.push_back(i);
3981 if ( event[i].
id() == 21 ) posInitGluon.push_back(i);
3982 else if ( event[i].idAbs() < 10 && event[i].
id() > 0)
3983 posInitQuark.push_back(i);
3984 else if ( event[i].idAbs() < 10 && event[i].
id() < 0)
3985 posInitAntiq.push_back(i);
3988 int nFiGluon = int(posFinalGluon.size());
3989 int nFiQuark = int(posFinalQuark.size());
3990 int nFiAntiq = int(posFinalAntiq.size());
3991 int nInGluon = int(posInitGluon.size());
3992 int nInQuark = int(posInitQuark.size());
3993 int nInAntiq = int(posInitAntiq.size());
3995 vector<Clustering> systems;
3999 for (
int i = 0; i < nFiGluon; ++i) {
4000 int EmtGluon = posFinalGluon[i];
4001 systems = findQCDTriple( EmtGluon, 2, event, posFinalPartn, posInitPartn);
4002 ret.insert(ret.end(), systems.begin(), systems.end());
4008 bool check_g2qq =
true;
4009 if ( ( ( nInQuark + nInAntiq == 0 )
4011 && (nFiQuark == 1) && (nFiAntiq == 1) )
4012 || ( ( nFiQuark + nFiAntiq == 0)
4013 && (nInQuark == 1) && (nInAntiq == 1) ) )
4020 for(
int i=0; i < nFiQuark; ++i) {
4021 int EmtQuark = posFinalQuark[i];
4022 systems = findQCDTriple( EmtQuark,1,event, posFinalPartn, posInitPartn);
4023 ret.insert(ret.end(), systems.begin(), systems.end());
4029 for(
int i=0; i < nFiAntiq; ++i) {
4030 int EmtAntiq = posFinalAntiq[i];
4031 systems = findQCDTriple( EmtAntiq,1,event, posFinalPartn, posInitPartn);
4032 ret.insert(ret.end(), systems.begin(), systems.end());
4044 void History::attachClusterings (vector<Clustering>& clus,
int iEmt,
int iRad,
4045 int iRec,
int iPartner,
double pT,
const Event& event) {
4047 if ( !mergingHooksPtr->doWeakClustering() ) {
4050 if (pT <= 0.)
return;
4051 clus.push_back( Clustering(iEmt, iRad, iRec, iPartner,
4052 pT, 0, 0, 0, 0, 9));
4057 int radSpin =
event[iRad].intPol();
4058 int emtSpin =
event[iEmt].intPol();
4059 int recSpin =
event[iRec].intPol();
4060 bool hasRadSpin = (radSpin != 9);
4061 bool hasEmtSpin = (emtSpin != 9);
4062 bool hasRecSpin = (recSpin != 9);
4065 bool radQuark =
event[iRad].idAbs() < 10;
4066 bool emtQuark =
event[iEmt].idAbs() < 10;
4067 bool recQuark =
event[iRec].idAbs() < 10;
4070 vector < vector<int> > structs;
4072 for (
int i = 0; i < 3; ++i){
4073 int sRad = (i==0) ? -1 : (i==1) ? 1 : 9;
4074 for (
int j = 0; j < 3; ++j){
4075 int sEmt = (j==0) ? -1 : (j==1) ? 1 : 9;
4076 for (
int k = 0; k < 3; ++k){
4077 int sRec = (k==0) ? -1 : (k==1) ? 1 : 9;
4079 s.push_back(sRad); s.push_back(sEmt); s.push_back(sRec);
4080 structs.push_back(s);
4085 vector < vector<int> > allStructs;
4086 for (
int i = 0; i< int(structs.size()); ++i) {
4088 if (hasRadSpin && radQuark && structs[i][0] != radSpin)
continue;
4089 if (hasEmtSpin && emtQuark && structs[i][1] != emtSpin)
continue;
4090 if (hasRecSpin && recQuark && structs[i][2] != recSpin)
continue;
4092 if (!hasRadSpin && radQuark && structs[i][0] == 9)
continue;
4093 if (!hasEmtSpin && emtQuark && structs[i][1] == 9)
continue;
4094 if (!hasRecSpin && recQuark && structs[i][2] == 9)
continue;
4096 if (!radQuark && structs[i][0] != radSpin)
continue;
4097 if (!emtQuark && structs[i][1] != emtSpin)
continue;
4098 if (!recQuark && structs[i][2] != recSpin)
continue;
4100 if (radQuark && emtQuark && structs[i][0] != structs[i][1])
continue;
4102 allStructs.push_back(structs[i]);
4106 int flavRadBef = getRadBeforeFlav(iRad, iEmt, event);
4108 for (
int i = 0; i< int(allStructs.size()); ++i) {
4110 int spinRadBef = getRadBeforeSpin(iRad, iEmt, allStructs[i][0],
4111 allStructs[i][1], event);
4112 clus.push_back( Clustering(iEmt, iRad, iRec, iPartner, pT, flavRadBef,
4113 allStructs[i][0], allStructs[i][1], allStructs[i][2], spinRadBef) );
4134 vector<Clustering> History::findQCDTriple (
int EmtTagIn,
int colTopIn,
4136 vector<int> posFinalPartn,
4137 vector <int> posInitPartn ) {
4140 int EmtTag = EmtTagIn;
4143 int colTop = colTopIn;
4146 int finalSize = int(posFinalPartn.size());
4147 int initSize = int(posInitPartn.size());
4148 int size = initSize + finalSize;
4150 vector<Clustering> clus;
4154 for (
int a = 0; a < size; ++a ) {
4155 int i = (a < finalSize)? a : (a - finalSize) ;
4156 int iRad = (a < finalSize)? posFinalPartn[i] : posInitPartn[i];
4158 if ( event[iRad].col() ==
event[EmtTag].col()
4159 &&
event[iRad].acol() ==
event[EmtTag].acol() )
4162 if (iRad != EmtTag ) {
4163 int pTdef =
event[iRad].isFinal() ? 1 : -1;
4164 int sign = (a < finalSize)? 1 : -1 ;
4170 if ( event[iRad].
id() == -sign*
event[EmtTag].id() ) {
4174 if (event[iRad].isFinal() ) {
4175 if (event[iRad].
id() < 0) {
4176 col =
event[EmtTag].col();
4179 acl =
event[EmtTag].acol();
4183 if (event[iRad].
id() < 0) {
4184 acl =
event[EmtTag].acol();
4187 col =
event[EmtTag].col();
4199 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
4202 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4206 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4207 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4213 if ( iRec != 0 && iPartner != 0
4214 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4215 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4216 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4223 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
4226 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4230 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4231 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4237 if ( iRec != 0 && iPartner != 0
4238 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4239 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4240 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4251 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
4254 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4258 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4259 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4265 if ( iRec != 0 && iPartner != 0
4266 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4267 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4268 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4275 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
4278 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4282 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4283 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4289 if ( iRec != 0 && iPartner != 0
4290 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4291 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4292 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4297 }
else if ( event[iRad].
id() == 21
4298 &&( event[iRad].col() == event[EmtTag].col()
4299 || event[iRad].acol() == event[EmtTag].acol() )) {
4305 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4306 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
4310 int col = getRadBeforeCol(iRad, EmtTag, event);
4311 int acl = getRadBeforeAcol(iRad, EmtTag, event);
4321 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
4322 ? event[iRad].col() :
event[iRad].acol();
4325 if (colRemove > 0 && col > 0 && col != colRemove)
4326 iPartner = FindCol(col,iRad,EmtTag,event,1,
true)
4327 + FindCol(col,iRad,EmtTag,event,2,
true);
4328 else if (colRemove > 0 && acl > 0 && acl != colRemove)
4329 iPartner = FindCol(acl,iRad,EmtTag,event,1,
true)
4330 + FindCol(acl,iRad,EmtTag,event,2,
true);
4332 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event ) ) {
4333 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
4334 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
4342 if ( (event[iRad].col() == event[EmtTag].acol())
4343 || (event[iRad].acol() == event[EmtTag].col())
4344 || (event[iRad].col() == event[EmtTag].col())
4345 || (event[iRad].acol() == event[EmtTag].acol()) ) {
4352 if (event[iRad].isFinal() ) {
4354 if ( event[iRad].
id() < 0) {
4355 acl =
event[EmtTag].acol();
4356 col =
event[iRad].col();
4357 }
else if ( event[iRad].
id() > 0 && event[iRad].
id() < 10) {
4358 col =
event[EmtTag].col();
4359 acl =
event[iRad].acol();
4361 col =
event[EmtTag].col();
4362 acl =
event[EmtTag].acol();
4367 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
4368 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4370 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4371 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4372 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4376 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
4377 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4379 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4380 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4381 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4387 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
4388 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4390 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4391 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4392 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4396 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
4397 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4399 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4400 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4401 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4415 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4416 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
4420 col = getRadBeforeCol(iRad, EmtTag, event);
4421 acl = getRadBeforeAcol(iRad, EmtTag, event);
4431 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
4432 ? event[iRad].col() : 0;
4433 iPartner = (colRemove > 0)
4434 ? FindCol(col,iRad,EmtTag,event,1,
true)
4435 + FindCol(col,iRad,EmtTag,event,2,
true)
4436 : FindCol(acl,iRad,EmtTag,event,1,true)
4437 + FindCol(acl,iRad,EmtTag,event,2,true);
4439 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event)) {
4440 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
4441 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
4461 vector<Clustering> History::getAllEWClusterings() {
4462 vector<Clustering> ret;
4465 vector<Clustering> systems;
4466 systems = getEWClusterings(state);
4467 ret.insert(ret.end(), systems.begin(), systems.end());
4478 vector<Clustering> History::getEWClusterings(
const Event& event) {
4479 vector<Clustering> ret;
4483 vector <int> posFinalPartn;
4484 vector <int> posInitPartn;
4485 vector <int> posFinalW;
4486 vector <int> posFinalZ;
4490 for (
int i=3; i <
event.size(); ++i )
4491 if ( event[i].isFinal() ) {
4493 posFinalPartn.push_back(i);
4496 posInitPartn.push_back(i);
4499 for (
int i=0; i <
event.size(); ++i )
4500 if ( event[i].isFinal() &&
event[i].idAbs() == 24 )
4501 posFinalW.push_back( i );
4504 for (
int i=0; i <
event.size(); ++i )
4505 if ( event[i].isFinal() &&
event[i].idAbs() == 23 )
4506 posFinalZ.push_back( i );
4509 vector<Clustering> systems;
4512 for (
int i = 0; i < int(posFinalW.size()); ++i ) {
4513 int emtW = posFinalW[i];
4514 systems = findEWTripleW( emtW, event, posFinalPartn, posInitPartn);
4515 ret.insert(ret.end(), systems.begin(), systems.end());
4520 for (
int i = 0; i < int(posFinalZ.size()); ++i ) {
4521 int emtZ = posFinalZ[i];
4523 systems = findEWTripleZ( emtZ, event, posFinalPartn, posInitPartn);
4524 ret.insert(ret.end(), systems.begin(), systems.end());
4543 vector<Clustering> History::findEWTripleW (
int emtTagIn,
const Event& event,
4544 vector<int> posFinalPartn, vector<int> posInitPartn ) {
4546 int emtTag = emtTagIn;
4547 int flavEmt =
event[emtTag].id();
4553 int finalSize = int(posFinalPartn.size());
4554 int initSize = int(posInitPartn.size());
4557 vector<int> flavCounts(30,0);
4559 for (
int a = 0; a < finalSize; ++a ) {
4560 if (event[posFinalPartn[a]].idAbs() < 20) {
4562 if (event[posFinalPartn[a]].
id() < 0)
4564 flavCounts[
event[posFinalPartn[a]].idAbs()] += sign;
4566 if (event[posFinalPartn[a]].idAbs() == 24)
4570 for (
int a = 0; a < initSize; ++a ) {
4571 if (event[posInitPartn[a]].idAbs() < 20) {
4573 if (event[posInitPartn[a]].
id() < 0)
4575 flavCounts[
event[posInitPartn[a]].idAbs()] -= sign;
4579 vector<Clustering> clus;
4583 for (
int a = 0; a < finalSize; ++a ) {
4585 int iRad = posFinalPartn[a];
4586 if (iRad != emtTag) {
4589 int spinRad =
event[iRad].intPol();
4590 if (spinRad == -1 || spinRad == 9 || spinRad == 0) {
4594 int flavRad =
event[iRad].id();
4596 if (event[iRad].isQuark() || event[iRad].isLepton()) {
4599 int flavExp = (flavRad > 0) ? 24 : -24;
4600 if (abs(flavRad) % 2 == 0) flavExp = -flavExp;
4601 if (flavExp == flavEmt) {
4604 vector<int> flavRadBefs = posFlavCKM(flavRad);
4608 for (
int i = 0;i < int(flavRadBefs.size()); ++i)
4609 flavRadBefs[i] = - flavRadBefs[i];
4612 for (
int i = 0; i < finalSize; ++i ) {
4613 int iRec = posFinalPartn[i];
4616 if ( iRec != iRad && iRec != emtTag ) {
4618 for (
int j = 0;j < int(flavRadBefs.size()); ++j) {
4621 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4622 flavRad, flavRadBefs[j], 1))
4625 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4626 pTLund(event, iRad, emtTag, iRec, pTdef, flavRadBefs[j]),
4627 flavRadBefs[j], -1 ) );
4638 for (
int a = 0;a < int(posInitPartn.size()); ++a) {
4639 int iRad = posInitPartn[a];
4640 int flavRad =
event[iRad].id();
4642 if (event[iRad].isQuark() || event[iRad].isLepton()) {
4645 int spinRad =
event[iRad].intPol();
4646 if (spinRad == -1 || spinRad == 9 || spinRad == 0) {
4649 int flavExp = (flavRad > 0) ? 24 : -24;
4650 if (abs(flavRad) % 2 == 1) flavExp = -flavExp;
4651 if (flavExp == flavEmt) {
4653 vector<int> flavRadBefs = posFlavCKM(flavRad);
4657 for (
int i = 0;i < int(flavRadBefs.size()); ++i)
4658 flavRadBefs[i] = - flavRadBefs[i];
4661 for (
int i = 0; i < int(posInitPartn.size()); ++i ) {
4662 int iRec = posInitPartn[i];
4665 if ( i != a && iRec != emtTag) {
4666 for (
int j = 0;j < int(flavRadBefs.size()); ++j) {
4670 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4671 flavRad, flavRadBefs[j], -1))
4673 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4674 pTLund(event, iRad, emtTag, iRec, -1, flavRadBefs[j]),
4675 flavRadBefs[j], -1 ) );
4700 vector<Clustering> History::findEWTripleZ (
int emtTagIn,
const Event& event,
4701 vector<int> posFinalPartn, vector<int> posInitPartn ) {
4703 int emtTag = emtTagIn;
4709 int finalSize = int(posFinalPartn.size());
4710 int initSize = int(posInitPartn.size());
4713 vector<int> flavCounts(30,0);
4715 for (
int a = 0; a < finalSize; ++a ) {
4716 if (event[posFinalPartn[a]].idAbs() < 20) {
4718 if (event[posFinalPartn[a]].
id() < 0)
4720 flavCounts[
event[posFinalPartn[a]].idAbs()] += sign;
4722 if (event[posFinalPartn[a]].idAbs() == 24)
4726 for (
int a = 0; a < initSize; ++a ) {
4727 if (event[posInitPartn[a]].idAbs() < 20) {
4729 if (event[posInitPartn[a]].
id() < 0)
4731 flavCounts[
event[posInitPartn[a]].idAbs()] -= sign;
4735 vector<Clustering> clus;
4740 for (
int a = 0; a < finalSize; ++a ) {
4742 int iRad = posFinalPartn[a];
4743 if (iRad != emtTag) {
4746 int flavRad =
event[iRad].id();
4749 if (event[iRad].isQuark() || event[iRad].isLepton())
4751 for (
int i = 0; i < finalSize; ++i ) {
4752 int iRec = posFinalPartn[i];
4755 if ( iRec != iRad && iRec != emtTag ) {
4758 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4759 flavRad, flavRad, 1))
4762 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4763 pTLund(event, iRad, emtTag, iRec, pTdef, flavRad),
4771 for (
int a = 0;a < int(posInitPartn.size()); ++a) {
4772 int iRad = posInitPartn[a];
4773 int flavRad =
event[iRad].id();
4775 if (event[iRad].isQuark() || event[iRad].isLepton()) {
4778 for (
int i = 0; i < int(posInitPartn.size()); ++i ) {
4779 int iRec = posInitPartn[i];
4782 if ( i != a && iRec != emtTag) {
4785 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4786 flavRad, flavRad, -1))
4788 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4789 pTLund(event, iRad, emtTag, iRec, -1, flavRad),
4807 vector<Clustering> History::getAllSQCDClusterings() {
4808 vector<Clustering> ret;
4811 vector<Clustering> systems;
4812 systems = getSQCDClusterings(state);
4813 ret.insert(ret.end(), systems.begin(), systems.end());
4824 vector<Clustering> History::getSQCDClusterings(
const Event& event) {
4825 vector<Clustering> ret;
4829 vector <int> posFinalPartn;
4830 vector <int> posInitPartn;
4832 vector <int> posFinalGluon;
4833 vector <int> posFinalQuark;
4834 vector <int> posFinalAntiq;
4835 vector <int> posInitGluon;
4836 vector <int> posInitQuark;
4837 vector <int> posInitAntiq;
4841 for (
int i=0; i <
event.size(); ++i)
4842 if ( event[i].isFinal() &&
event[i].colType() !=0 ) {
4844 posFinalPartn.push_back(i);
4845 if ( event[i].
id() == 21 || event[i].
id() == 1000021)
4846 posFinalGluon.push_back(i);
4847 else if ( (event[i].idAbs() < 10 && event[i].
id() > 0)
4848 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4849 && event[i].
id() > 0)
4850 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4851 && event[i].
id() > 0))
4852 posFinalQuark.push_back(i);
4853 else if ( (event[i].idAbs() < 10 && event[i].
id() < 0)
4854 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4855 && event[i].
id() < 0)
4856 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4857 && event[i].
id() < 0))
4858 posFinalAntiq.push_back(i);
4859 }
else if ( event[i].status() == -21 && event[i].colType() != 0 ) {
4861 posInitPartn.push_back(i);
4862 if ( event[i].
id() == 21 || event[i].
id() == 1000021)
4863 posInitGluon.push_back(i);
4864 else if ( (event[i].idAbs() < 10 && event[i].
id() > 0)
4865 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4866 && event[i].
id() > 0)
4867 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4868 && event[i].
id() > 0))
4869 posInitQuark.push_back(i);
4870 else if ( (event[i].idAbs() < 10 && event[i].
id() < 0)
4871 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4872 && event[i].
id() < 0)
4873 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4874 && event[i].
id() < 0))
4875 posInitAntiq.push_back(i);
4878 int nFiGluon = int(posFinalGluon.size());
4879 int nFiQuark = int(posFinalQuark.size());
4880 int nFiAntiq = int(posFinalAntiq.size());
4881 int nInGluon = int(posInitGluon.size());
4882 int nInQuark = int(posInitQuark.size());
4883 int nInAntiq = int(posInitAntiq.size());
4884 vector<Clustering> systems;
4888 for (
int i = 0; i < nFiGluon; ++i) {
4889 int EmtGluon = posFinalGluon[i];
4890 systems = findSQCDTriple( EmtGluon, 2, event, posFinalPartn, posInitPartn);
4891 ret.insert(ret.end(), systems.begin(), systems.end());
4897 bool check_g2qq =
true;
4898 if ( ( ( nInQuark + nInAntiq == 0 )
4900 && (nFiQuark == 1) && (nFiAntiq == 1) )
4901 || ( ( nFiQuark + nFiAntiq == 0)
4902 && (nInQuark == 1) && (nInAntiq == 1) ) )
4909 for(
int i=0; i < nFiQuark; ++i) {
4910 int EmtQuark = posFinalQuark[i];
4911 systems = findSQCDTriple( EmtQuark,1,event, posFinalPartn, posInitPartn);
4912 ret.insert(ret.end(), systems.begin(), systems.end());
4918 for(
int i=0; i < nFiAntiq; ++i) {
4919 int EmtAntiq = posFinalAntiq[i];
4920 systems = findSQCDTriple( EmtAntiq,1,event, posFinalPartn, posInitPartn);
4921 ret.insert(ret.end(), systems.begin(), systems.end());
4942 vector<Clustering> History::findSQCDTriple (
int EmtTagIn,
int colTopIn,
4944 vector<int> posFinalPartn,
4945 vector <int> posInitPartn ) {
4948 int EmtTag = EmtTagIn;
4951 int colTop = colTopIn;
4954 int offsetL = 1000000;
4955 int offsetR = 2000000;
4958 int finalSize = int(posFinalPartn.size());
4959 int initSize = int(posInitPartn.size());
4960 int size = initSize + finalSize;
4962 vector<Clustering> clus;
4966 for (
int a = 0; a < size; ++a ) {
4967 int i = (a < finalSize)? a : (a - finalSize) ;
4968 int iRad = (a < finalSize)? posFinalPartn[i] : posInitPartn[i];
4970 if ( event[iRad].col() ==
event[EmtTag].col()
4971 &&
event[iRad].acol() ==
event[EmtTag].acol() )
4975 int radID =
event[iRad].id();
4977 bool isSQCDrad = (abs(radID) > offsetL);
4979 bool isSQCDemt = (
event[EmtTag].idAbs() > offsetL );
4981 if (iRad != EmtTag ) {
4982 int pTdef =
event[iRad].isFinal() ? 1 : -1;
4983 int sign = (a < finalSize)? 1 : -1 ;
4986 int radBefID = getRadBeforeFlav(iRad,EmtTag,event);
4987 if ( pTdef == -1 && abs(radBefID) > offsetL )
continue;
4993 int radSign = (
event[iRad].id() < 0) ? -1 : 1;
4994 int emtSign = (
event[EmtTag].id() < 0) ? -1 : 1;
4997 bool finalSplitting =
false;
4998 if ( abs(radID) < 10
4999 && radSign*(abs(radID)+offsetL) == -sign*event[EmtTag].
id() )
5000 finalSplitting =
true;
5001 if ( abs(radID) < 10
5002 && radSign*(abs(radID)+offsetR) == -sign*event[EmtTag].
id() )
5003 finalSplitting =
true;
5004 if ( abs(radID) > offsetL && abs(radID) < offsetL+10
5005 && radID == -sign*emtSign*(
event[EmtTag].idAbs() + offsetL) )
5006 finalSplitting =
true;
5007 if ( abs(radID) > offsetR && abs(radID) < offsetR+10
5008 && radID == -sign*emtSign*(
event[EmtTag].idAbs() + offsetR) )
5009 finalSplitting =
true;
5012 bool initialSplitting =
false;
5013 if ( radID == 21 && ( ( event[EmtTag].idAbs() > offsetL
5014 && event[EmtTag].idAbs() < offsetL+10)
5015 || (
event[EmtTag].idAbs() > offsetR
5016 &&
event[EmtTag].idAbs() < offsetR+10) )
5017 && (
event[iRad].col() ==
event[EmtTag].col()
5018 ||
event[iRad].acol() ==
event[EmtTag].acol() ) )
5019 initialSplitting =
true;
5021 if ( finalSplitting ) {
5025 if ( radID < 0 && event[iRad].colType() == -1) {
5026 acl =
event[EmtTag].acol();
5027 col =
event[iRad].acol();
5028 }
else if ( event[iRad].colType() == 1 ) {
5029 col =
event[EmtTag].col();
5030 acl =
event[iRad].col();
5040 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
5043 if ( (sign < 0) && (event[iRec].isFinal()) ) {
5047 for(
int l = 0; l < int(posInitPartn.size()); ++l)
5048 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
5056 if ( !isSQCDrad && !isSQCDemt
5057 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
5060 if ( iRec != 0 && iPartner != 0
5061 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
5062 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
5063 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5070 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
5073 if ( (sign < 0) && (event[iRec].isFinal()) ) {
5077 for(
int l = 0; l < int(posInitPartn.size()); ++l)
5078 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
5086 if ( !isSQCDrad && !isSQCDemt
5087 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
5090 if ( iRec != 0 && iPartner != 0
5091 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
5092 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
5093 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5103 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
5106 if ( (sign < 0) && (event[iRec].isFinal()) ) {
5110 for(
int l = 0; l < int(posInitPartn.size()); ++l)
5111 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
5119 if ( !isSQCDrad && !isSQCDemt
5120 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
5123 if ( iRec != 0 && iPartner != 0
5124 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
5125 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
5126 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5133 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
5136 if ( (sign < 0) && (event[iRec].isFinal()) ) {
5140 for(
int l = 0; l < int(posInitPartn.size()); ++l)
5141 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
5149 if ( !isSQCDrad && !isSQCDemt
5150 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
5153 if ( iRec != 0 && iPartner != 0
5154 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
5155 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
5156 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5161 }
else if ( initialSplitting ) {
5164 if ( !isSQCDrad && !isSQCDemt )
continue;
5171 for(
int l = 0; l < int(posInitPartn.size()); ++l)
5172 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
5176 int col = getRadBeforeCol(iRad, EmtTag, event);
5177 int acl = getRadBeforeAcol(iRad, EmtTag, event);
5187 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
5188 ? event[iRad].col() : 0;
5191 if (colRemove > 0 && col > 0)
5192 iPartner = FindCol(col,iRad,EmtTag,event,1,
true)
5193 + FindCol(col,iRad,EmtTag,event,2,
true);
5194 else if (colRemove > 0 && acl > 0)
5195 iPartner = FindCol(acl,iRad,EmtTag,event,1,
true)
5196 + FindCol(acl,iRad,EmtTag,event,2,
true);
5198 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event ) ) {
5199 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
5200 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
5209 if ( (event[iRad].col() == event[EmtTag].acol())
5210 || (event[iRad].acol() == event[EmtTag].col())
5211 || (event[iRad].col() == event[EmtTag].col())
5212 || (event[iRad].acol() == event[EmtTag].acol()) ) {
5219 if (event[iRad].isFinal() ) {
5221 if ( radID < 0 && event[iRad].colType() == -1) {
5222 acl =
event[EmtTag].acol();
5223 col =
event[iRad].col();
5224 }
else if ( radID > 0 && event[iRad].colType() == 1 ) {
5225 col =
event[EmtTag].col();
5226 acl =
event[iRad].acol();
5228 col =
event[iRad].col();
5229 acl =
event[iRad].acol();
5234 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
5235 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
5237 if ( !isSQCDrad && !isSQCDemt
5238 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
5241 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
5242 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
5243 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5247 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
5248 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
5250 if ( !isSQCDrad && !isSQCDemt
5251 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
5254 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
5255 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
5256 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5262 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
5263 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
5265 if ( !isSQCDrad && !isSQCDemt
5266 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
5269 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
5270 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
5271 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5275 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
5276 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
5278 if ( !isSQCDrad && !isSQCDemt
5279 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
5282 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
5283 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
5284 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
5296 if ( !isSQCDrad || !isSQCDemt )
continue;
5304 for(
int l = 0; l < int(posInitPartn.size()); ++l)
5305 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
5309 col = getRadBeforeCol(iRad, EmtTag, event);
5310 acl = getRadBeforeAcol(iRad, EmtTag, event);
5320 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
5321 ? event[iRad].col() : 0;
5322 iPartner = (colRemove > 0)
5323 ? FindCol(col,iRad,EmtTag,event,1,
true)
5324 + FindCol(col,iRad,EmtTag,event,2,
true)
5325 : FindCol(acl,iRad,EmtTag,event,1,true)
5326 + FindCol(acl,iRad,EmtTag,event,2,true);
5328 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event)) {
5329 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
5330 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
5351 double History::getProb(
const Clustering & SystemIn) {
5354 int Rad = SystemIn.emittor;
5355 int Rec = SystemIn.recoiler;
5356 int Emt = SystemIn.emitted;
5359 double showerProb = 0.0;
5363 if (SystemIn.pT() <= 0.)
return 0.;
5372 bool isFSR = (state[Rad].isFinal() && state[Rec].isFinal());
5373 bool isFSRinREC = (state[Rad].isFinal() && !state[Rec].isFinal());
5374 bool isISR = !state[Rad].isFinal();
5377 if ( mergingHooksPtr->useShowerPlugin() ) {
5378 int iPartner = (isISR && SystemIn.partner > 0) ? SystemIn.partner : Rec;
5381 bool isFSR2 = showers->timesPtr->isTimelike(state, Rad, Emt, iPartner,
"");
5383 string name = showers->timesPtr->getSplittingName( state, Rad, Emt,
5385 pr = showers->timesPtr->getSplittingProb( state, Rad, Emt,
5388 string name = showers->spacePtr->getSplittingName(state, Rad, Emt,
5390 pr = showers->spacePtr->getSplittingProb(state, Rad, Emt,
5399 for(
int i=0; i < state.size(); ++i)
5400 if (state[i].isFinal()) nFinal++;
5401 bool isLast = (nFinal == (mergingHooksPtr->hardProcess->nQuarksOut()
5402 +mergingHooksPtr->hardProcess->nLeptonOut()+1));
5405 bool isElectroWeak = (state[Emt].idAbs() == 23 || state[Emt].idAbs() == 24);
5411 for(
int i=0;i< int(state.size()); ++i) {
5412 if (state[i].mother1() == 1) inP = i;
5413 if (state[i].mother1() == 2) inM = i;
5416 Vec4 sum = state[Rad].p() + state[Rec].p() - state[Emt].p();
5417 double m2Dip = sum.m2Calc();
5418 double sHat = (state[inM].p() + state[inP].p()).m2Calc();
5420 double z1 = m2Dip / sHat;
5422 Vec4 Q1( state[Rad].p() - state[Emt].p() );
5423 Vec4 Q2( state[Rec].p() - state[Emt].p() );
5425 double Q1sq = -Q1.m2Calc();
5427 double pT1sq = pow2(pTLund(state, Rad, Emt, Rec, -1, SystemIn.flavRadBef));
5430 bool g2QQmassive = mergingHooksPtr->includeMassive()
5431 && state[Rad].id() == 21
5432 && ( (state[Emt].idAbs() >= 4 && state[Emt].idAbs() < 7)
5433 || (state[Emt].idAbs() > 1000000 && state[Emt].idAbs() < 1000010 )
5434 || (state[Emt].idAbs() > 2000000 && state[Emt].idAbs() < 2000010 ));
5435 bool Q2Qgmassive = mergingHooksPtr->includeMassive()
5436 && state[Emt].id() == 21
5437 && ( (state[Rad].idAbs() >= 4 && state[Rad].idAbs() < 7)
5438 || (state[Rad].idAbs() > 1000000 && state[Rad].idAbs() < 1000010 )
5439 || (state[Rad].idAbs() > 2000000 && state[Rad].idAbs() < 2000010 ));
5440 bool isMassive = mergingHooksPtr->includeMassive()
5441 && ( g2QQmassive || Q2Qgmassive
5442 || state[Emt].id() == 1000021);
5443 double m2Emt0 = state[Emt].p().m2Calc();
5444 double m2Rad0 = pow(particleDataPtr->m0(state[Rad].id()),2);
5447 if ( g2QQmassive) Q1sq += m2Emt0;
5448 else if (Q2Qgmassive) Q1sq += m2Rad0;
5451 double pT0sq = pow(mergingHooksPtr->pT0ISR(),2);
5452 double Q2sq = -Q2.m2Calc();
5455 bool g2QQmassiveRec = mergingHooksPtr->includeMassive()
5456 && state[Rec].id() == 21
5457 && ( state[Emt].idAbs() >= 4 && state[Emt].idAbs() < 7);
5458 bool Q2QgmassiveRec = mergingHooksPtr->includeMassive()
5459 && state[Emt].id() == 21
5460 && ( state[Rec].idAbs() >= 4 && state[Rec].idAbs() < 7);
5461 double m2Rec0 = pow(particleDataPtr->m0(state[Rec].id()),2);
5462 if ( g2QQmassiveRec) Q2sq += m2Emt0;
5463 else if (Q2QgmassiveRec) Q2sq += m2Rec0;
5465 bool hasJoinedEvol = (state[Emt].id() == 21
5466 || state[Rad].id() == state[Rec].id());
5471 if ( mergingHooksPtr->pickByFull() || mergingHooksPtr->pickBySumPT()) {
5472 double facJoined = ( Q2sq + pT0sq/(1.-z1) )
5473 * 1./(Q1sq*Q2sq + pT0sq*sHat + pow(pT0sq/(1.-z1),2));
5474 double facSingle = mergingHooksPtr->nonJoinedNorm()*1./( pT1sq + pT0sq);
5476 fac = (hasJoinedEvol && isLast) ? facJoined : facSingle;
5478 }
else if (mergingHooksPtr->pickByPoPT2()) {
5479 fac = 1./(pT1sq + pT0sq);
5481 string message=
"Error in History::getProb: Scheme for calculating";
5482 message+=
" shower splitting probability is undefined.";
5483 infoPtr->errorMsg(message);
5489 if ( isElectroWeak ) {
5498 }
else if ( state[Emt].
id() == 21 && state[Rad].
id() != 21) {
5500 double num = CF*(1. + pow(z1,2)) / (1.-z1);
5501 if (isMassive) num -= CF * z1 * (1.-z1) * (m2Rad0/pT1sq);
5504 double meReweighting = 1.;
5508 for(
int i=0; i < state.size(); ++i)
5509 if (state[i].isFinal() && state[i].colType() != 0
5510 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5515 &&
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1) {
5516 double sH = m2Dip / z1;
5518 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5519 double misMatch = (uH*tH - (uH + tH)*pT0sq/(1.-z1)
5520 + pow(pT0sq/(1.-z1),2) ) / (uH*tH);
5521 meReweighting *= (tH*tH + uH*uH + 2. * m2Dip * sH)
5522 / (sH*sH + m2Dip*m2Dip);
5523 meReweighting *= misMatch;
5526 showerProb = num*fac*meReweighting;
5529 }
else if ( state[Emt].
id() == 21 && state[Rad].id() == 21) {
5531 double num = 2.*CA*pow2(1. - z1*(1.-z1)) / (z1*(1.-z1));
5535 double meReweighting = 1.;
5539 for(
int i=0; i < state.size(); ++i)
5540 if (state[i].isFinal() && state[i].colType() != 0
5541 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5546 && mergingHooksPtr->getProcessString().compare(
"pp>h") == 0
5547 && int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1 ) {
5548 double sH = m2Dip / z1;
5550 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5551 meReweighting *= 0.5
5552 * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(m2Dip))
5553 / pow2(sH*sH - m2Dip * (sH - m2Dip));
5557 showerProb = num*fac*meReweighting;
5560 }
else if ( state[Emt].
id() != 21 && state[Rad].id() != 21 ) {
5562 double num = CF*(1. + pow2(1.-z1)) / z1;
5566 double meReweighting = 1.;
5570 for (
int i=0; i < state.size(); ++i )
5571 if ( state[i].isFinal() && state[i].colType() != 0
5572 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5577 && mergingHooksPtr->getProcessString().compare(
"pp>h") == 0
5578 && int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1) {
5579 double sH = m2Dip / z1;
5580 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5581 meReweighting *= (sH*sH + uH*uH)
5582 / (sH*sH + pow2(sH -m2Dip));
5586 showerProb = num*fac*meReweighting;
5589 }
else if ( state[Emt].
id() != 21 && state[Rad].id() == 21 ) {
5592 double num = TR * ( pow(z1,2) + pow(1.-z1,2) );
5593 if (isMassive) num += TR * 2.*z1*(1.-z1)*(m2Emt0/pT1sq);
5595 double meReweighting = 1.;
5599 for (
int i=0; i < state.size(); ++i )
5600 if ( state[i].isFinal() && state[i].colType() != 0
5601 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5606 &&
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1) {
5607 double sH = m2Dip / z1;
5609 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5611 double misMatch = ( uH - pT0sq/(1.-z1) ) / uH;
5612 double me = (sH*sH + uH*uH + 2. * m2Dip * tH)
5613 / (pow2(sH - m2Dip) + m2Dip*m2Dip);
5615 meReweighting *= me;
5616 meReweighting *= misMatch;
5619 showerProb = num*fac*meReweighting;
5623 string message =
"Error in History::getProb: Splitting kernel"
5624 " undefined in ISR in clustering.";
5625 infoPtr->errorMsg(message);
5629 double m2Sister = pow(state[Emt].m(),2);
5630 double pT2corr = (Q1sq - z1*(m2Dip + Q1sq)*(Q1sq + m2Sister)/m2Dip);
5631 if (pT2corr < 0.) showerProb = 0.0;
5635 if ( state[Emt].
id() == state[Rad].id()
5636 && ( state[Rad].idAbs() == 4 || state[Rad].idAbs() == 5 )) {
5637 double m2QQsister = 2.*4.*m2Sister;
5638 double pT2QQcorr = Q1sq - z1*(m2Dip + Q1sq)*(Q1sq + m2QQsister)
5640 if (pT2QQcorr < 0.0) showerProb = 0.0;
5645 = pow2(mergingHooksPtr->settingsPtr->parm(
"SpaceShower:pTmin"));
5647 double zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
5648 ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
5649 zMaxAbs = min(1.,zMaxAbs);
5651 double zMinAbs = 2. * state[Rad].e() / state[0].e() * z1;
5654 int radBefID = getRadBeforeFlav(Rad, Emt, state);
5655 if ( abs(radBefID) == 4 || abs(radBefID) == 5 ) {
5656 double m2Massive = pow2(particleDataPtr->m0(radBefID));
5657 double mRatio = sqrt( m2Massive / m2Dip );
5658 double zMaxMassive = (1. - mRatio) / ( 1. + mRatio * (1. - mRatio) );
5659 zMaxAbs = min(zMaxAbs, zMaxMassive);
5662 if (z1 < zMinAbs || z1 > zMaxAbs) showerProb = 0.0;
5664 if (mergingHooksPtr->includeRedundant()) {
5666 AlphaStrong* asISR = mergingHooksPtr->AlphaS_ISR();
5667 double as = (*asISR).alphaS(pT1sq + pT0sq) / (2.*M_PI);
5674 }
else if (isFSR || isFSRinREC) {
5677 int recSign = (state[Rec].isFinal()) ? 1 : -1;
5678 Vec4 sum = state[Rad].p() + recSign*state[Rec].p() + state[Emt].p();
5679 double m2Dip = abs(sum.m2Calc());
5682 Vec4 Q1( state[Rad].p() + state[Emt].p() );
5683 Vec4 Q2( state[Rec].p() + state[Emt].p() );
5686 double z1 = getCurrentZ( Rad, Rec, Emt, SystemIn.flavRadBef);
5689 double Q1sq = Q1.m2Calc();
5691 double pT1sq = pow(SystemIn.pT(),2);
5693 double Q2sq = Q2.m2Calc();
5696 bool isMassiveRad = ( state[Rad].idAbs() >= 4
5697 && state[Rad].id() != 21 );
5698 bool isMassiveRec = ( state[Rec].idAbs() >= 4
5699 && state[Rec].id() != 21 );
5702 double m2Rad0 = pow(particleDataPtr->m0(state[Rad].id()),2);
5703 double m2Rec0 = pow(particleDataPtr->m0(state[Rec].id()),2);
5704 if ( mergingHooksPtr->includeMassive() && isMassiveRad ) Q1sq -= m2Rad0;
5705 if ( mergingHooksPtr->includeMassive() && isMassiveRec ) Q2sq -= m2Rec0;
5710 if ( mergingHooksPtr->pickByFull() || mergingHooksPtr->pickBySumPT()) {
5711 double facJoined = (1.-z1)/Q1sq * m2Dip/( Q1sq + Q2sq );
5712 double facSingle = mergingHooksPtr->fsrInRecNorm() * 1./ pT1sq;
5713 fac = (!isFSRinREC && isLast) ? facJoined : facSingle;
5715 }
else if (mergingHooksPtr->pickByPoPT2()) {
5718 string message=
"Error in History::getProb: Scheme for calculating";
5719 message+=
" shower splitting probability is undefined.";
5720 infoPtr->errorMsg(message);
5725 if ( isElectroWeak ) {
5734 }
else if ( state[Emt].
id() == 21 && state[Rad].colType() == 2) {
5737 double num = 0.5* CA * (1. + pow3(z1)) / (1.-z1);
5739 showerProb = num*fac;
5742 }
else if ( state[Emt].
id() == 21
5743 && state[Rad].colType() != 2 && state[Rec].colType() != 2) {
5746 double num = CF * 2./(1.-z1);
5750 for(
int i=0; i < state.size(); ++i)
5751 if (state[i].isFinal() && state[i].colType() != 0
5752 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5756 ||
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) > 1)
5757 num = CF * (1. + pow2(z1)) /(1.-z1);
5759 double meReweighting = 1.;
5763 &&
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1 ) {
5765 double x1 = 2. * (sum * state[Rad].p()) / m2Dip;
5766 double x2 = 2. * (sum * state[Rec].p()) / m2Dip;
5767 double prop1 = max(1e-12, 1. - x1);
5768 double prop2 = max(1e-12, 1. - x2);
5769 double x3 = max(1e-12, 2. - x1 - x2);
5771 double ShowerRate1 = 2./( x3 * prop2 );
5772 double meDividingFactor1 = prop1 / x3;
5773 double me = (pow(x1,2) + pow(x2,2))/(prop1*prop2);
5774 meReweighting = meDividingFactor1 * me / ShowerRate1;
5777 showerProb = num*fac*meReweighting;
5780 }
else if ( state[Emt].
id() == 21 && state[Rad].colType() != 2
5781 && state[Rec].colType() == 2 ) {
5785 double num = CF * (1. + pow2(z1)) / (1.-z1);
5786 showerProb = num*fac;
5789 }
else if ( state[Emt].
id() != 21 ) {
5791 int flavour = state[Emt].id();
5794 double mFlavour = particleDataPtr->m0(flavour);
5796 double mDipole = m(state[Rad].p(), state[Emt].p());
5798 double beta = sqrtpos( 1. - 4.*pow2(mFlavour)/pow2(mDipole) );
5800 double num = 0.5*TR * ( z1*z1 + (1.-z1)*(1.-z1) );
5801 if (beta <= 0.) beta = 0.;
5803 showerProb = num*fac*beta;
5807 string message=
"Error in History::getProb: Splitting kernel undefined";
5808 message+=
" in FSR clustering.";
5809 infoPtr->errorMsg(message);
5812 if (mergingHooksPtr->includeRedundant()) {
5814 AlphaStrong* asFSR = mergingHooksPtr->AlphaS_FSR();
5815 double as = (*asFSR).alphaS(pT1sq) / (2.*M_PI);
5820 double m2DipCorr = pow2(sqrt(m2Dip) - sqrt(m2Rec0)) - m2Rad0;
5821 double zMin = 0.5 - sqrtpos( 0.25 - pT1sq / m2DipCorr );
5822 double m2 = m2Rad0 + 2. * state[Rad].p()*state[Emt].p();
5823 bool keepMassive = (z1 > zMin && z1 < 1. - zMin
5824 && m2 * m2Dip < z1 * (1. - z1) * pow2(m2Dip + m2 - m2Rec0) );
5826 if (!keepMassive) showerProb *= 0.0;
5830 string message=
"Error in History::getProb: Radiation could not be";
5831 message+=
" interpreted as FSR or ISR.";
5832 infoPtr->errorMsg(message);
5835 if (mergingHooksPtr->doWeakClustering()) {
5842 if (state[Rad].idAbs() < 10 && state[Rad].intPol() == 9
5843 && SystemIn.spinRad != 9) factor *= 0.5;
5844 if (state[Emt].idAbs() < 10 && state[Emt].intPol() == 9
5845 && SystemIn.spinEmt != 9) factor *= 0.5;
5846 if (state[Rec].idAbs() < 10 && state[Rec].intPol() == 9
5847 && SystemIn.spinRec != 9) factor *= 0.5;
5848 if ( state[Emt].colType() != 0 ) {
5849 showerProb *= factor;
5855 if ( state[Emt].idAbs() < 10 && state[Rad].colType() != 0) {
5878 void History::setupBeams() {
5883 if (state.size() < 4)
return;
5885 if ( state[3].colType() == 0 )
return;
5886 if ( state[4].colType() == 0 )
return;
5892 for(
int i=0;i< int(state.size()); ++i) {
5893 if (state[i].mother1() == 1) inP = i;
5894 if (state[i].mother1() == 2) inM = i;
5899 int motherPcompRes = -1;
5900 int motherMcompRes = -1;
5902 bool sameFlavP =
false;
5903 bool sameFlavM =
false;
5908 for(
int i=0;i< int(mother->state.size()); ++i) {
5909 if (mother->state[i].mother1() == 1) inMotherP = i;
5910 if (mother->state[i].mother1() == 2) inMotherM = i;
5912 sameFlavP = (state[inP].id() == mother->state[inMotherP].id());
5913 sameFlavM = (state[inM].id() == mother->state[inMotherM].id());
5915 motherPcompRes = (sameFlavP) ? beamA[0].companion() : -2;
5916 motherMcompRes = (sameFlavM) ? beamB[0].companion() : -2;
5924 double Ep = 2. * state[inP].e();
5925 double Em = 2. * state[inM].e();
5928 if (state[inP].m() != 0. || state[inM].m() != 0.) {
5929 Ep = state[inP].pPos() + state[inM].pPos();
5930 Em = state[inP].pNeg() + state[inM].pNeg();
5934 double x1 = Ep / state[inS].m();
5935 beamA.append( inP, state[inP].
id(), x1);
5936 double x2 = Em / state[inS].m();
5937 beamB.append( inM, state[inM].
id(), x2);
5941 double scalePDF = (mother) ? scale : infoPtr->QFac();
5945 beamA.xfISR( 0, state[inP].
id(), x1, scalePDF*scalePDF);
5947 beamA.pickValSeaComp();
5949 beamA[0].companion(motherPcompRes);
5951 beamB.xfISR( 0, state[inM].
id(), x2, scalePDF*scalePDF);
5953 beamB.pickValSeaComp();
5955 beamB[0].companion(motherMcompRes);
5965 double History::pdfForSudakov() {
5968 if ( state[3].colType() == 0 )
return 1.0;
5969 if ( state[4].colType() == 0 )
return 1.0;
5972 bool FSR = ( mother->state[clusterIn.emittor].isFinal()
5973 && mother->state[clusterIn.recoiler].isFinal());
5974 bool FSRinRec = ( mother->state[clusterIn.emittor].isFinal()
5975 && !mother->state[clusterIn.recoiler].isFinal());
5978 if (FSR)
return 1.0;
5980 int iInMother = (FSRinRec)? clusterIn.recoiler : clusterIn.emittor;
5982 int side = ( mother->state[iInMother].pz() > 0 ) ? 1 : -1;
5986 for(
int i=0;i< int(state.size()); ++i) {
5987 if (state[i].mother1() == 1) inP = i;
5988 if (state[i].mother1() == 2) inM = i;
5992 int idMother = mother->state[iInMother].id();
5994 int iDau = (side == 1) ? inP : inM;
5995 int idDaughter = state[iDau].id();
5997 double xMother = 2. * mother->state[iInMother].e() / mother->state[0].e();
5999 double xDaughter = 2.*state[iDau].e() / state[0].e();
6002 double ratio = getPDFratio(side,
true,
false, idMother, xMother, scale,
6003 idDaughter, xDaughter, scale);
6008 return ( (FSRinRec)? min(1.,ratio) : ratio);
6016 double History::hardProcessME(
const Event& event ) {
6019 if (isEW2to1(event)) {
6022 if (event[5].idAbs() == 24) {
6023 int idIn1 =
event[3].id();
6024 int idIn2 =
event[4].id();
6025 double mW = particleDataPtr->m0(24);
6026 double gW = particleDataPtr->mWidth(24) / mW;
6027 double sH = (
event[3].p()+
event[4].p()).m2Calc();
6029 double thetaWRat = 1. / (12. * coupSMPtr->sin2thetaW());
6030 double ckmW = coupSMPtr->V2CKMid(abs(idIn1), abs(idIn2));
6032 double bwW = 12. * M_PI / ( pow2(sH - pow2(mW)) + pow2(sH * gW) );
6033 double preFac = thetaWRat * sqrt(sH) * particleDataPtr->mWidth(24);
6034 return ckmW * preFac * bwW;
6037 else if (event[5].idAbs() == 23) {
6038 double mZ = particleDataPtr->m0(23);
6039 double gZ = particleDataPtr->mWidth(23) / mZ;
6040 double sH = (
event[3].p()+
event[4].p()).m2Calc();
6042 double thetaZRat = (pow2(coupSMPtr->rf( abs(clusterIn.flavRadBef))) +
6043 pow2(coupSMPtr->lf( abs(clusterIn.flavRadBef)))) /
6044 (24. * coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW());
6045 double bwW = 12. * M_PI / ( pow2(sH - pow2(mZ)) + pow2(sH * gZ) );
6046 double preFac = thetaZRat * sqrt(sH) * particleDataPtr->mWidth(23);
6047 return preFac * bwW;
6050 string message=
"Warning in History::hardProcessME: Only Z/W are";
6051 message+=
" supported as 2->1 processes. Skipping history.";
6052 infoPtr->errorMsg(message);
6057 else if (isQCD2to2(event)) {
6058 int idIn1 =
event[3].id();
6059 int idIn2 =
event[4].id();
6060 int idOut1 =
event[5].id();
6061 int idOut2 =
event[6].id();
6063 double sH = (
event[3].p()+
event[4].p()).m2Calc();
6064 double tH = (
event[3].p()-
event[5].p()).m2Calc();
6065 double uH = (
event[3].p()-
event[6].p()).m2Calc();
6069 if (!(abs(idIn1) < 10 || abs(idIn1) == 21) ) isQCD =
false;
6070 if (!(abs(idIn2) < 10 || abs(idIn2) == 21) ) isQCD =
false;
6071 if (!(abs(idOut1) < 10 || abs(idOut1) == 21) ) isQCD =
false;
6072 if (!(abs(idOut2) < 10 || abs(idOut2) == 21) ) isQCD =
false;
6075 double cor = M_PI / (9. * pow2(sH));
6082 if (abs(idIn1) == 21 && abs(idIn2) == 21) {
6083 if (abs(idOut1) == 21 && abs(idOut2) == 21)
6084 return cor * simpleWeakShowerMEs.getMEgg2gg(sH, tH, uH);
6085 else return cor * simpleWeakShowerMEs.getMEgg2qqbar(sH, tH, uH);
6088 }
else if (abs(idIn1) == 21 || abs(idIn2) == 21) {
6089 if (idIn1 != idOut1) swap(uH, tH);
6090 return cor * simpleWeakShowerMEs.getMEqg2qg(sH, tH, uH);
6095 if (abs(idOut1) == 21 && abs(idOut2) == 21)
6096 return cor * simpleWeakShowerMEs.getMEqqbar2gg(sH, tH, uH);
6097 if (idIn1 == -idIn2) {
6098 if (abs(idIn1) == abs(idOut1)) {
6099 if (idIn1 != idOut1) swap(uH, tH);
6100 return cor * simpleWeakShowerMEs.getMEqqbar2qqbar(sH,tH,uH,
true);
6103 return cor * simpleWeakShowerMEs.getMEqqbar2qqbar(sH,tH,uH,
false);
6106 else if (idIn1 == idIn2)
6107 return cor * simpleWeakShowerMEs.getMEqq2qq(sH, tH, uH,
true);
6109 if (idIn1 == idOut1) swap(uH,tH);
6110 return cor * simpleWeakShowerMEs.getMEqq2qq(sH, tH, uH,
false);
6118 string process = mergingHooksPtr->getProcessString();
6121 if ( process.compare(
"pp>e+ve") == 0
6122 || process.compare(
"pp>e-ve~") == 0
6123 || process.compare(
"pp>LEPTONS,NEUTRINOS") == 0 ) {
6126 for (
int i=0; i < int(event.size()); ++i )
6127 if ( event[i].isFinal() ) nFinal++;
6128 if ( nFinal != 2 )
return 1.;
6130 double mW = particleDataPtr->m0(24);
6131 double gW = particleDataPtr->mWidth(24) / mW;
6133 int inP = (
event[3].pz() > 0) ? 3 : 4;
6134 int inM = (
event[3].pz() > 0) ? 4 : 3;
6137 for (
int i=0; i < int(event.size()); ++i ) {
6138 if ( event[i].isFinal() &&
event[i].px() > 0 ) outP = i;
6141 double sH = (
event[inP].p() +
event[inM].p()).m2Calc();
6142 double tH = (
event[inP].p() -
event[outP].p()).m2Calc();
6143 double uH = - sH - tH;
6146 result = ( 1. + (tH - uH)/sH ) / ( pow2(sH - mW*mW) + pow2(sH*gW) );
6148 result = mergingHooksPtr->hardProcessME(event);
6161 Event History::cluster( Clustering & inSystem ) {
6164 int Rad = inSystem.emittor;
6165 int Rec = inSystem.recoiler;
6166 int Emt = inSystem.emitted;
6168 double eCM = state[0].e();
6170 int radType = state[Rad].isFinal() ? 1 : -1;
6171 int recType = state[Rec].isFinal() ? 1 : -1;
6175 NewEvent.init(
"(hard process-modified)", particleDataPtr);
6179 if ( mergingHooksPtr->useShowerPlugin() ) {
6180 int iPartner = (radType == -1 && inSystem.partner > 0)
6181 ? inSystem.partner : Rec;
6182 bool isFSR = showers->timesPtr->isTimelike(state, Rad, Emt, iPartner,
"");
6184 string name = showers->timesPtr->getSplittingName(state,Rad,Emt,
6186 return showers->timesPtr->clustered(state,Rad,Emt,iPartner,name);
6188 string name = showers->spacePtr->getSplittingName(state,Rad,Emt,
6190 return showers->spacePtr->clustered(state,Rad,Emt,iPartner,name);
6195 for (
int i = 0; i < state.size(); ++i)
6196 if ( i != Rad && i != Rec && i != Emt )
6197 NewEvent.append( state[i] );
6200 for (
int i = 0; i < state.sizeJunction(); ++i)
6201 NewEvent.appendJunction( state.getJunction(i) );
6203 double mu = choseHardScale(state);
6205 NewEvent.saveSize();
6206 NewEvent.saveJunctionSize();
6208 NewEvent.scaleSecond(mu);
6212 Particle RecBefore = Particle( state[Rec] );
6213 RecBefore.setEvtPtr(&NewEvent);
6214 RecBefore.daughters(0,0);
6215 RecBefore.pol(inSystem.spinRec);
6217 int radID = inSystem.flavRadBef;
6218 if (radID == 0) radID = getRadBeforeFlav(Rad, Emt, state);
6219 int recID = state[Rec].id();
6220 Particle RadBefore = Particle( state[Rad] );
6221 RadBefore.setEvtPtr(&NewEvent);
6222 RadBefore.id(radID);
6223 RadBefore.pol(inSystem.spinRadBef);
6224 RadBefore.daughters(0,0);
6226 RadBefore.cols(RecBefore.acol(),RecBefore.col());
6229 if ( particleDataPtr->isResonance(radID) && radType == 1)
6230 RadBefore.status(state[Rad].status());
6233 double radMass = particleDataPtr->m0(radID);
6234 double recMass = particleDataPtr->m0(recID);
6235 if (radType == 1 ) RadBefore.m(radMass);
6236 else RadBefore.m(0.0);
6237 if (recType == 1 ) RecBefore.m(recMass);
6238 else RecBefore.m(0.0);
6242 if ( radType + recType == 2 && state[Emt].idAbs() != 23 &&
6243 state[Emt].idAbs() != 24) {
6246 Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p();
6247 double eCMME = sum.mCalc();
6253 double mDsq = pow2(eCMME);
6255 double mRsq = (radID == state[Rad].id() )
6256 ? abs(state[Rad].p().m2Calc())
6257 : pow2(particleDataPtr->m0(radID));
6258 double mSsq = abs(state[Rec].p().m2Calc());
6259 double a = 0.5*sqrt(mDsq);
6260 double b = 0.25*(mRsq - mSsq) / a;
6261 double c = sqrt(pow2(a) + pow2(b) - 2.*a*b - mSsq);
6263 Rad4mom.p( 0., 0., c, a+b);
6264 Rec4mom.p( 0., 0.,-c, a-b);
6267 Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p());
6268 Vec4 old2 = Vec4(state[Rec].p());
6269 RotBstMatrix fromCM;
6270 fromCM.fromCMframe(old1, old2);
6272 Rad4mom.rotbst(fromCM);
6273 Rec4mom.rotbst(fromCM);
6275 RadBefore.p(Rad4mom);
6276 RecBefore.p(Rec4mom);
6277 RadBefore.m(sqrt(mRsq));
6278 RecBefore.m(sqrt(mSsq));
6280 }
else if ( radType + recType == 2 && (state[Emt].idAbs() == 23 ||
6281 state[Emt].idAbs() == 24) ) {
6284 Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p();
6285 double eCMME = sum.mCalc();
6291 double mDsq = pow2(eCMME);
6293 double mRsq = (radID == state[Rad].id() )
6294 ? abs(state[Rad].p().m2Calc())
6295 : pow2(particleDataPtr->m0(radID));
6296 double mSsq = abs(state[Rec].p().m2Calc());
6297 double a = 0.5*sqrt(mDsq);
6298 double b = 0.25*(mRsq - mSsq) / a;
6299 double c = sqrt(pow2(a) + pow2(b) - 2.*a*b - mSsq);
6301 Rad4mom.p( 0., 0., c, a+b);
6302 Rec4mom.p( 0., 0.,-c, a-b);
6305 Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p());
6306 Vec4 old2 = Vec4(state[Rec].p());
6307 RotBstMatrix fromCM;
6308 fromCM.fromCMframe(old1, old2);
6310 Rad4mom.rotbst(fromCM);
6311 Rec4mom.rotbst(fromCM);
6313 RadBefore.p(Rad4mom);
6314 RecBefore.p(Rec4mom);
6315 RadBefore.m(sqrt(mRsq));
6316 RecBefore.m(sqrt(mSsq));
6318 }
else if ( radType + recType == 0 ) {
6322 Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p();
6323 double eCMME = sum.mCalc();
6328 double mDsq = pow2(eCMME);
6330 double mRsq = (radID == state[Rad].id() )
6331 ? abs(state[Rad].p().m2Calc())
6332 : pow2(particleDataPtr->m0(radID));
6333 double mSsq = abs(state[Rec].p().m2Calc());
6334 double a = 0.5*sqrt(mDsq);
6335 double b = 0.25*(mRsq - mSsq) / a;
6336 double c = sqrt(pow2(a) + pow2(b) - 2.*a*b - mSsq);
6338 Rad4mom.p( 0., 0., c, a+b);
6339 Rec4mom.p( 0., 0.,-c, a-b);
6342 Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p());
6343 Vec4 old2 = Vec4(state[Rec].p());
6344 RotBstMatrix fromCM;
6345 fromCM.fromCMframe(old1, old2);
6347 Rad4mom.rotbst(fromCM);
6348 Rec4mom.rotbst(fromCM);
6351 Rec4mom = 2.*state[Rec].p() - Rec4mom;
6355 if ( abs(Rec4mom.mCalc()) > 1e-7 ) {
6356 double pzSign = (Rec4mom.pz() > 0.) ? 1. : -1.;
6357 double eRec = Rec4mom.e();
6358 Rec4mom.p(0., 0., pzSign*eRec, eRec);
6361 RadBefore.p(Rad4mom);
6362 RecBefore.p(Rec4mom);
6363 RadBefore.m(sqrt(mRsq));
6376 Vec4 pMother( state[Rad].p() );
6377 Vec4 pSister( state[Emt].p() );
6378 Vec4 pPartner( state[Rec].p() );
6379 Vec4 pDaughterBef( 0.,0.,0.,0. );
6380 Vec4 pRecoilerBef( 0.,0.,0.,0. );
6381 Vec4 pDaughter( 0.,0.,0.,0. );
6382 Vec4 pRecoiler( 0.,0.,0.,0. );
6385 int sign = (state[Rad].pz() > 0.) ? 1 : -1;
6389 double phi = pSister.phi();
6391 RotBstMatrix rot_by_mphi;
6392 rot_by_mphi.rot(0.,-phi);
6394 RotBstMatrix rot_by_pphi;
6395 rot_by_pphi.rot(0.,phi);
6399 double x1 = 2. * pMother.e() / eCM;
6401 double x2 = 2. * pPartner.e() / eCM;
6404 Vec4 qDip( pMother - pSister);
6405 Vec4 qAfter(pMother + pPartner);
6406 Vec4 qBefore(qDip + pPartner);
6407 double z = qBefore.m2Calc() / qAfter.m2Calc();
6410 double x1New = z*x1;
6412 double sHat = x1New*x2New*eCM*eCM;
6417 pDaughterBef.p( 0., 0., sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
6418 pRecoilerBef.p( 0., 0., -sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
6421 pMother.rotbst( rot_by_mphi );
6422 pSister.rotbst( rot_by_mphi );
6423 pPartner.rotbst( rot_by_mphi );
6424 for(
int i=3; i< NewEvent.size(); ++i)
6425 NewEvent[i].rotbst( rot_by_mphi );
6429 pDaughter.p( pMother - pSister);
6430 pRecoiler.p( pPartner );
6431 RotBstMatrix from_CM_to_DRoff;
6433 from_CM_to_DRoff.toCMframe(pDaughter, pRecoiler);
6435 from_CM_to_DRoff.toCMframe(pRecoiler, pDaughter);
6439 pMother.rotbst( from_CM_to_DRoff );
6440 pPartner.rotbst( from_CM_to_DRoff );
6441 pSister.rotbst( from_CM_to_DRoff );
6442 for(
int i=3; i< NewEvent.size(); ++i)
6443 NewEvent[i].rotbst( from_CM_to_DRoff );
6448 RotBstMatrix from_DR_to_CM;
6449 from_DR_to_CM.bst( 0., 0., sign*( x1New - x2New ) / ( x1New + x2New ) );
6454 pDaughterBef.rotbst( from_DR_to_CM );
6455 pRecoilerBef.rotbst( from_DR_to_CM );
6456 for(
int i=3; i< NewEvent.size(); ++i)
6457 NewEvent[i].rotbst( from_DR_to_CM );
6460 for(
int i=3; i< NewEvent.size(); ++i)
6461 NewEvent[i].rotbst( rot_by_pphi );
6465 if ( abs(pRecoilerBef.mCalc()) > 1e-7 ) {
6466 double pzSign = (pRecoilerBef.pz() > 0.) ? 1. : -1.;
6467 double eRec = pRecoilerBef.e();
6468 pRecoilerBef.p(0., 0., pzSign*eRec, eRec);
6470 if ( abs(pDaughterBef.mCalc()) > 1e-7 ) {
6471 double pzSign = (pDaughterBef.pz() > 0.) ? 1. : -1.;
6472 double eDau = pDaughterBef.e();
6473 pDaughterBef.p(0., 0., pzSign*eDau, eDau);
6477 RecBefore.p( pRecoilerBef );
6478 RadBefore.p( pDaughterBef );
6479 if (RecBefore.pz() > 0.) RecBefore.mother1(1);
6480 else RecBefore.mother1(2);
6481 if (RadBefore.pz() > 0.) RadBefore.mother1(1);
6482 else RadBefore.mother1(2);
6487 RecBefore.scale(mu);
6488 RadBefore.scale(mu);
6491 NewEvent.append(RecBefore);
6495 int emtID = state[Emt].id();
6496 if ( emtID == 22 || emtID == 23 || abs(emtID) == 24 )
6497 RadBefore.cols( state[Rad].col(), state[Rad].acol() );
6499 else if ( !connectRadiator( RadBefore, radType, RecBefore, recType,
6509 outState.init(
"(hard process-modified)", particleDataPtr);
6513 for (
int i = 0; i < 3; ++i)
6514 outState.append( NewEvent[i] );
6516 for (
int i = 0; i < state.sizeJunction(); ++i)
6517 outState.appendJunction( state.getJunction(i) );
6519 outState.saveSize();
6520 outState.saveJunctionSize();
6522 outState.scaleSecond(mu);
6523 bool radAppended =
false;
6524 bool recAppended =
false;
6525 int size = int(outState.size());
6527 int radPos = 0, recPos = 0;
6530 if ( RecBefore.mother1() == 1) {
6531 recPos = outState.append( RecBefore );
6533 }
else if ( RadBefore.mother1() == 1 ) {
6534 radPos = outState.append( RadBefore );
6539 for(
int i=0; i < int(state.size()); ++i)
6540 if (state[i].mother1() == 1) in1 =i;
6541 outState.append( state[in1] );
6545 if ( RecBefore.mother1() == 2) {
6546 recPos = outState.append( RecBefore );
6548 }
else if ( RadBefore.mother1() == 2 ) {
6549 radPos = outState.append( RadBefore );
6554 for(
int i=0; i < int(state.size()); ++i)
6555 if (state[i].mother1() == 2) in2 =i;
6557 outState.append( state[in2] );
6562 if (!recAppended && !RecBefore.isFinal()) {
6564 recPos = outState.append( RecBefore);
6567 if (!radAppended && !RadBefore.isFinal()) {
6569 radPos = outState.append( RadBefore);
6576 for (
int i = 0; i < int(NewEvent.size()-1); ++i)
6577 if (NewEvent[i].status() == -22) outState.append( NewEvent[i] );
6579 for (
int i = 0; i < int(NewEvent.size()-1); ++i)
6580 if (NewEvent[i].status() == 22) outState.append( NewEvent[i] );
6582 if (!radAppended && RadBefore.statusAbs() == 22)
6583 radPos = outState.append(RadBefore);
6585 recPos = outState.append(RecBefore);
6586 if (!radAppended && RadBefore.statusAbs() != 22)
6587 radPos = outState.append(RadBefore);
6589 for(
int i = 0; i < int(NewEvent.size()-1); ++i)
6590 if ( NewEvent[i].status() != 22
6591 && NewEvent[i].colType() != 0
6592 && NewEvent[i].isFinal())
6593 outState.append( NewEvent[i] );
6595 for(
int i = 0; i < int(NewEvent.size()-1); ++i)
6596 if ( NewEvent[i].status() != 22
6597 && NewEvent[i].colType() == 0
6598 && NewEvent[i].isFinal() )
6599 outState.append( NewEvent[i]);
6602 vector<int> posIntermediate;
6603 vector<int> posDaughter1;
6604 vector<int> posDaughter2;
6605 for(
int i=0; i < int(outState.size()); ++i)
6606 if (outState[i].status() == -22) {
6607 posIntermediate.push_back(i);
6608 int d1 = outState[i].daughter1();
6609 int d2 = outState[i].daughter2();
6611 int daughter1 = FindParticle( state[d1], outState);
6612 int daughter2 = FindParticle( state[d2], outState);
6617 posDaughter1.push_back( daughter1);
6620 while(!outState[daughter1].isFinal() ) daughter1++;
6621 posDaughter1.push_back( daughter1);
6624 posDaughter2.push_back( daughter2);
6626 daughter2 = outState.size()-1;
6627 while(!outState[daughter2].isFinal() ) daughter2--;
6628 posDaughter2.push_back( daughter2);
6632 for(
int i=0; i < int(posIntermediate.size()); ++i) {
6633 outState[posIntermediate[i]].daughters(posDaughter1[i],posDaughter2[i]);
6634 outState[posDaughter1[i]].mother1(posIntermediate[i]);
6635 outState[posDaughter2[i]].mother1(posIntermediate[i]);
6639 int minParFinal = int(outState.size());
6640 int maxParFinal = 0;
6641 for(
int i=0; i < int(outState.size()); ++i)
6642 if (outState[i].mother1() == 3 && outState[i].mother2() == 4) {
6643 minParFinal = min(i,minParFinal);
6644 maxParFinal = max(i,maxParFinal);
6647 if (minParFinal == maxParFinal) maxParFinal = 0;
6648 outState[3].daughters(minParFinal,maxParFinal);
6649 outState[4].daughters(minParFinal,maxParFinal);
6652 outState.saveSize();
6653 outState.saveJunctionSize();
6656 inSystem.recBef = recPos;
6657 inSystem.radBef = radPos;
6671 if ( radType == -1 && state[Rad].colType() == 1) {
6673 for(
int i=0; i < int(state.size()); ++i)
6674 if ( i != Rad && i != Emt && i != Rec
6675 && state[i].status() == -22
6676 && state[i].col() == state[Rad].col() )
6678 }
else if ( radType == -1 && state[Rad].colType() == -1) {
6680 for(
int i=0; i < int(state.size()); ++i)
6681 if ( i != Rad && i != Emt && i != Rec
6682 && state[i].status() == -22
6683 && state[i].acol() == state[Rad].acol() )
6685 }
else if ( radType == 1 && state[Rad].colType() == 1) {
6687 for(
int i=0; i < int(state.size()); ++i)
6688 if ( i != Rad && i != Emt && i != Rec
6689 && state[i].status() == -22
6690 && state[i].col() == state[Rad].col() )
6692 }
else if ( radType == 1 && state[Rad].colType() == -1) {
6694 for(
int i=0; i < int(state.size()); ++i)
6695 if ( i != Rad && i != Emt && i != Rec
6696 && state[i].status() == -22
6697 && state[i].acol() == state[Rad].acol() )
6703 int iColResNow = FindParticle( state[iColRes], outState);
6706 int radCol = outState[radPos].col();
6707 int radAcl = outState[radPos].acol();
6709 int resCol = outState[iColResNow].col();
6710 int resAcl = outState[iColResNow].acol();
6712 bool matchesRes = (radCol > 0
6713 && ( radCol == resCol || radCol == resAcl))
6715 && ( radAcl == resCol || radAcl == resAcl));
6719 if (!matchesRes && iColResNow > 0) {
6720 if ( radType == -1 && outState[radPos].colType() == 1)
6721 outState[iColResNow].col(radCol);
6722 else if ( radType ==-1 && outState[radPos].colType() ==-1)
6723 outState[iColResNow].acol(radAcl);
6724 else if ( radType == 1 && outState[radPos].colType() == 1)
6725 outState[iColResNow].col(radCol);
6726 else if ( radType == 1 && outState[radPos].colType() ==-1)
6727 outState[iColResNow].acol(radAcl);
6734 if (!matchesRes && iColResNow > 0 && iColRes != iColResNow)
6735 outState[radPos].mother1(iColResNow);
6740 if ( !validEvent(outState) ) {
6747 iReclusteredNew = radPos;
6761 int History::getRadBeforeFlav(
const int RadAfter,
const int EmtAfter,
6762 const Event& event) {
6764 int type =
event[RadAfter].isFinal() ? 1 :-1;
6765 int emtID =
event[EmtAfter].id();
6766 int radID =
event[RadAfter].id();
6767 int emtCOL =
event[EmtAfter].col();
6768 int radCOL =
event[RadAfter].col();
6769 int emtACL =
event[EmtAfter].acol();
6770 int radACL =
event[RadAfter].acol();
6772 bool colConnected = ((type == 1) && ( (emtCOL !=0 && (emtCOL ==radACL))
6773 || (emtACL !=0 && (emtACL ==radCOL)) ))
6774 ||((type ==-1) && ( (emtCOL !=0 && (emtCOL ==radCOL))
6775 || (emtACL !=0 && (emtACL ==radACL)) ));
6781 if ( type == 1 && emtID == -radID && !colConnected )
6784 if ( type ==-1 && radID == 21 )
6787 if ( type ==-1 && !colConnected
6788 && radID != 21 && abs(emtID) < 10 && abs(radID) < 10)
6792 int radSign = (radID < 0) ? -1 : 1;
6793 int offsetL = 1000000;
6794 int offsetR = 2000000;
6796 if ( emtID == 1000021 ) {
6798 if (abs(radID) < 10 ) {
6799 int offset = offsetL;
6802 for (
int i=0; i < int(event.size()); ++i)
6803 if ( event[i].isFinal()
6804 &&
event[i].idAbs() < offsetR+10 &&
event[i].idAbs() > offsetR)
6806 return radSign*(abs(radID)+offset);
6809 if (abs(radID) > offsetL && abs(radID) < offsetL+10 )
6810 return radSign*(abs(radID)-offsetL);
6811 if (abs(radID) > offsetR && abs(radID) < offsetR+10 )
6812 return radSign*(abs(radID)-offsetR);
6814 if (radID == 21 )
return emtID;
6817 int emtSign = (emtID < 0) ? -1 : 1;
6820 if ( abs(emtID) > offsetL && abs(emtID) < offsetL+10 )
6821 emtOffset = offsetL;
6822 if ( abs(emtID) > offsetR && abs(emtID) < offsetR+10 )
6823 emtOffset = offsetR;
6825 if ( abs(radID) > offsetL && abs(radID) < offsetL+10 )
6826 radOffset = offsetL;
6827 if ( abs(radID) > offsetR && abs(radID) < offsetR+10 )
6828 radOffset = offsetR;
6831 if ( type == 1 && !colConnected ) {
6833 if ( emtOffset > 0 && radOffset == 0
6834 && emtSign*(abs(emtID) - emtOffset) == -radID )
6837 if ( emtOffset == 0 && radOffset > 0
6838 && emtID == -radSign*(abs(radID) - radOffset) )
6843 if ( type ==-1 && radID == 1000021 ) {
6845 if ( emtOffset > 0 )
return -emtSign*(abs(emtID) - emtOffset);
6847 else return -emtSign*(abs(emtID) + emtOffset);
6852 && ( (abs(emtID) > offsetL && abs(emtID) < offsetL+10)
6853 || (abs(emtID) > offsetR && abs(emtID) < offsetR+10))
6854 && ( (abs(radID) > offsetL && abs(radID) < offsetL+10)
6855 || (abs(radID) > offsetR && abs(radID) < offsetR+10))
6856 && emtSign*(abs(emtID)+emtOffset) == radSign*(abs(radID) - radOffset)
6857 && !colConnected ) {
6863 double m2final = (
event[RadAfter].p()+
event[EmtAfter].p()).m2Calc();
6865 if ( emtID == 22 || emtID == 23 )
return radID;
6867 if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) <= 10. )
6870 if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) > 10. )
6873 if ( type ==-1 && (radID == 22 || radID == 23) )
6876 if ( type ==-1 && abs(emtID) < 10 && abs(radID) < 10 && colConnected )
6881 if ( emtID == 24 && radID < 0 )
return radID + 1;
6882 if ( emtID == 24 && radID > 0 )
return radID + 1;
6886 if ( emtID ==-24 && radID < 0 )
return radID - 1;
6887 if ( emtID ==-24 && radID > 0 )
return radID - 1;
6901 int History::getRadBeforeSpin(
const int radAfter,
const int emtAfter,
6902 const int spinRadAfter,
const int spinEmtAfter,
6903 const Event& event) {
6906 int radBeforeFlav = getRadBeforeFlav(radAfter, emtAfter, event);
6909 if ( event[radAfter].isFinal()
6910 && event[radAfter].
id() == -event[emtAfter].
id())
6911 return (spinRadAfter == 9) ? spinEmtAfter : spinRadAfter;
6914 if ( event[radAfter].isFinal() && abs(radBeforeFlav) < 10
6915 && event[radAfter].idAbs() < 10)
6917 return spinRadAfter;
6920 if ( event[radAfter].isFinal() && abs(radBeforeFlav) < 10
6921 && event[emtAfter].idAbs() < 10)
6923 return spinEmtAfter;
6926 if ( event[radAfter].isFinal() && radBeforeFlav == 21
6927 && event[radAfter].
id() == 21)
6929 return (spinRadAfter == 9) ? spinEmtAfter : spinRadAfter;
6932 if ( !event[radAfter].isFinal()
6933 && radBeforeFlav == -event[emtAfter].
id())
6934 return (spinRadAfter == 9) ? spinEmtAfter : spinRadAfter;
6937 if ( !event[radAfter].isFinal() && abs(radBeforeFlav) < 10
6938 && event[radAfter].idAbs() < 10)
6940 return spinRadAfter;
6943 if ( !event[radAfter].isFinal() && radBeforeFlav == 21
6944 && event[emtAfter].idAbs() < 10)
6946 return spinEmtAfter;
6965 bool History::connectRadiator( Particle& Radiator,
const int RadType,
6966 const Particle& Recoiler,
const int RecType,
6967 const Event& event ) {
6970 Radiator.cols( -1, -1 );
6974 if ( Radiator.colType() == -1 ) {
6977 if ( RadType + RecType == 2 )
6978 Radiator.cols( 0, Recoiler.col());
6979 else if ( RadType + RecType == 0 )
6980 Radiator.cols( 0, Recoiler.acol());
6988 for (
int i = 0; i <
event.size(); ++i) {
6989 int col =
event[i].col();
6990 int acl =
event[i].acol();
6992 if ( event[i].isFinal()) {
6994 if ( acl > 0 && FindCol(acl,i,0,event,1,
true) == 0
6995 && FindCol(acl,i,0,event,2,
true) == 0 )
6996 Radiator.acol(event[i].acol());
6999 if ( col > 0 && FindCol(col,i,0,event,1,
true) == 0
7000 && FindCol(col,i,0,event,2,
true) == 0 )
7001 Radiator.acol(event[i].col());
7006 }
else if ( Radiator.colType() == 1 ) {
7009 if ( RadType + RecType == 2 )
7010 Radiator.cols( Recoiler.acol(), 0);
7012 else if ( RadType + RecType == 0 )
7013 Radiator.cols( Recoiler.col(), 0);
7022 for (
int i = 0; i <
event.size(); ++i) {
7023 int col =
event[i].col();
7024 int acl =
event[i].acol();
7026 if ( event[i].isFinal()) {
7028 if ( col > 0 && FindCol(col,i,0,event,1,
true) == 0
7029 && FindCol(col,i,0,event,2,
true) == 0)
7030 Radiator.col(event[i].col());
7033 if ( acl > 0 && FindCol(acl,i,0,event,1,
true) == 0
7034 && FindCol(acl,i,0,event,2,
true) == 0)
7035 Radiator.col(event[i].acol());
7041 }
else if ( Radiator.colType() == 2 ) {
7047 for (
int i = 0; i <
event.size(); ++i) {
7048 int col =
event[i].col();
7049 int acl =
event[i].acol();
7052 if ( event[i].isFinal()) {
7053 if ( col > 0 && FindCol(col,iEx,0,event,1,
true) == 0
7054 && FindCol(col,iEx,0,event,2,
true) == 0) {
7055 if (Radiator.status() < 0 ) Radiator.col(event[i].col());
7056 else Radiator.acol(event[i].col());
7058 if ( acl > 0 && FindCol(acl,iEx,0,event,2,
true) == 0
7059 && FindCol(acl,iEx,0,event,1,
true) == 0 ) {
7060 if (Radiator.status() < 0 ) Radiator.acol(event[i].acol());
7061 else Radiator.col(event[i].acol());
7064 if ( col > 0 && FindCol(col,iEx,0,event,1,
true) == 0
7065 && FindCol(col,iEx,0,event,2,
true) == 0) {
7066 if (Radiator.status() < 0 ) Radiator.acol(event[i].col());
7067 else Radiator.col(event[i].col());
7069 if ( acl > 0 && (FindCol(acl,iEx,0,event,2,
true) == 0
7070 && FindCol(acl,iEx,0,event,1,
true) == 0)) {
7071 if (Radiator.status() < 0 ) Radiator.col(event[i].acol());
7072 else Radiator.acol(event[i].acol());
7079 if (Radiator.col() < 0 || Radiator.acol() < 0)
return false;
7101 int History::FindCol(
int col,
int iExclude1,
int iExclude2,
7102 const Event& event,
int type,
bool isHardIn) {
7104 bool isHard = isHardIn;
7109 for(
int n = 0; n <
event.size(); ++n) {
7110 if ( n != iExclude1 && n != iExclude2
7111 && event[n].colType() != 0
7112 &&(
event[n].status() > 0
7113 ||
event[n].status() == -21) ) {
7114 if ( event[n].acol() == col ) {
7118 if ( event[n].col() == col ) {
7127 for(
int n = 0; n <
event.size(); ++n) {
7128 if ( n != iExclude1 && n != iExclude2
7129 && event[n].colType() != 0
7130 &&(
event[n].status() == 43
7131 ||
event[n].status() == 51
7132 ||
event[n].status() == -41
7133 ||
event[n].status() == -42) ) {
7134 if ( event[n].acol() == col ) {
7138 if ( event[n].col() == col ) {
7146 if ( type == 1 && index < 0)
return abs(index);
7147 if ( type == 2 && index > 0)
return abs(index);
7161 int History::FindParticle(
const Particle& particle,
const Event& event,
7162 bool checkStatus ) {
7166 for (
int i =
int(event.size()) - 1; i > 0; --i )
7167 if ( event[i].
id() == particle.id()
7168 &&
event[i].colType() == particle.colType()
7169 &&
event[i].chargeType() == particle.chargeType()
7170 &&
event[i].col() == particle.col()
7171 &&
event[i].acol() == particle.acol()
7172 &&
event[i].charge() == particle.charge() ) {
7177 if (index >= 0 && checkStatus && event[index].status() != particle.status())
7193 int History::getRadBeforeCol(
const int rad,
const int emt,
7194 const Event& event) {
7197 int type = (
event[rad].isFinal()) ? 1 :-1;
7199 int radBeforeFlav = getRadBeforeFlav(rad,emt,event);
7201 int radBeforeCol = -1;
7203 if (radBeforeFlav == 21) {
7206 if (type == 1 && event[emt].
id() != 21) {
7207 radBeforeCol = (
event[rad].col() > 0)
7208 ? event[rad].col() :
event[emt].col();
7210 }
else if (type == -1 && event[emt].
id() != 21) {
7211 radBeforeCol = (
event[rad].col() > 0)
7212 ? event[rad].col() :
event[emt].acol();
7214 }
else if (type == 1 && event[emt].
id() == 21) {
7217 int colRemove = (
event[rad].col() ==
event[emt].acol())
7218 ? event[rad].col() :
event[rad].acol();
7219 radBeforeCol = (
event[rad].col() == colRemove)
7220 ? event[emt].col() :
event[rad].col();
7222 }
else if (type == -1 && event[emt].
id() == 21) {
7225 int colRemove = (
event[rad].col() ==
event[emt].col())
7226 ? event[rad].col() :
event[rad].acol();
7227 radBeforeCol = (
event[rad].col() == colRemove)
7228 ? event[emt].acol() :
event[rad].col();
7232 }
else if (radBeforeFlav > 0) {
7235 if (type == 1 && event[emt].
id() != 21) {
7238 int colRemove = (
event[rad].col() ==
event[emt].acol())
7239 ? event[rad].acol() : 0;
7240 radBeforeCol = (
event[rad].col() == colRemove)
7241 ? event[emt].col() :
event[rad].col();
7243 }
else if (type == 1 && event[emt].
id() == 21) {
7246 int colRemove = (
event[rad].col() ==
event[emt].acol())
7247 ? event[rad].col() : 0;
7248 radBeforeCol = (
event[rad].col() == colRemove)
7249 ? event[emt].col() :
event[rad].col();
7251 }
else if (type == -1 && event[emt].
id() != 21) {
7254 int colRemove = (
event[rad].col() ==
event[emt].col())
7255 ? event[rad].col() : 0;
7256 radBeforeCol = (
event[rad].col() == colRemove)
7257 ? event[emt].acol() :
event[rad].col();
7259 }
else if (type == -1 && event[emt].
id() == 21) {
7262 int colRemove = (
event[rad].col() ==
event[emt].col())
7263 ? event[rad].col() : 0;
7264 radBeforeCol = (
event[rad].col() == colRemove)
7265 ? event[emt].acol() :
event[rad].col();
7272 return radBeforeCol;
7285 int History::getRadBeforeAcol(
const int rad,
const int emt,
7286 const Event& event) {
7289 int type = (
event[rad].isFinal()) ? 1 :-1;
7292 int radBeforeFlav = getRadBeforeFlav(rad,emt,event);
7294 int radBeforeAcl = -1;
7296 if (radBeforeFlav == 21) {
7299 if (type == 1 && event[emt].
id() != 21) {
7300 radBeforeAcl = (
event[rad].acol() > 0)
7301 ? event[rad].acol() :
event[emt].acol();
7303 }
else if (type == -1 && event[emt].
id() != 21) {
7304 radBeforeAcl = (
event[rad].acol() > 0)
7305 ? event[rad].acol() :
event[emt].col();
7307 }
else if (type == 1 && event[emt].
id() == 21) {
7310 int colRemove = (
event[rad].col() ==
event[emt].acol())
7311 ? event[rad].col() :
event[rad].acol();
7312 radBeforeAcl = (
event[rad].acol() == colRemove)
7313 ? event[emt].acol() :
event[rad].acol();
7315 }
else if (type == -1 && event[emt].
id() == 21) {
7318 int colRemove = (
event[rad].col() ==
event[emt].col())
7319 ? event[rad].col() :
event[rad].acol();
7320 radBeforeAcl = (
event[rad].acol() == colRemove)
7321 ? event[emt].col() :
event[rad].acol();
7325 }
else if (radBeforeFlav < 0) {
7328 if (type == 1 && event[emt].
id() != 21) {
7331 int colRemove = (
event[rad].col() ==
event[emt].acol())
7332 ? event[rad].acol() : 0;
7333 radBeforeAcl = (
event[rad].acol() == colRemove)
7334 ? event[emt].acol() :
event[rad].acol();
7336 }
else if (type == 1 && event[emt].
id() == 21) {
7339 int colRemove = (
event[rad].acol() ==
event[emt].col())
7340 ? event[rad].acol() : 0;
7341 radBeforeAcl = (
event[rad].acol() == colRemove)
7342 ? event[emt].acol() :
event[rad].acol();
7344 }
else if (type == -1 && event[emt].
id() != 21) {
7347 int colRemove = (
event[rad].acol() ==
event[emt].acol())
7348 ? event[rad].acol() : 0;
7349 radBeforeAcl = (
event[rad].acol() == colRemove)
7350 ? event[emt].col() :
event[rad].acol();
7352 }
else if (type == -1 && event[emt].
id() == 21) {
7355 int colRemove = (
event[rad].acol() ==
event[emt].acol())
7356 ? event[rad].acol() : 0;
7357 radBeforeAcl = (
event[rad].acol() == colRemove)
7358 ? event[emt].col() :
event[rad].acol();
7365 return radBeforeAcl;
7377 int History::getColPartner(
const int in,
const Event& event) {
7379 if (event[in].col() == 0)
return 0;
7383 partner = FindCol(event[in].col(),in,0,event,1,
true);
7386 partner = FindCol(event[in].col(),in,0,event,2,
true);
7401 int History::getAcolPartner(
const int in,
const Event& event) {
7403 if (event[in].acol() == 0)
return 0;
7407 partner = FindCol(event[in].acol(),in,0,event,2,
true);
7410 partner = FindCol(event[in].acol(),in,0,event,1,
true);
7427 vector<int> History::getReclusteredPartners(
const int rad,
const int emt,
7428 const Event& event) {
7431 int type =
event[rad].isFinal() ? 1 : -1;
7433 int radBeforeCol = getRadBeforeCol(rad, emt, event);
7434 int radBeforeAcl = getRadBeforeAcol(rad, emt, event);
7436 vector<int> partners;
7442 for(
int i=0; i < int(event.size()); ++i) {
7444 if ( i != emt && i != rad
7445 && event[i].status() == -21
7446 &&
event[i].col() > 0
7447 &&
event[i].col() == radBeforeCol)
7448 partners.push_back(i);
7450 if ( i != emt && i != rad
7451 && event[i].isFinal()
7452 && event[i].acol() > 0
7453 && event[i].acol() == radBeforeCol)
7454 partners.push_back(i);
7456 if ( i != emt && i != rad
7457 && event[i].status() == -21
7458 && event[i].acol() > 0
7459 && event[i].acol() == radBeforeAcl)
7460 partners.push_back(i);
7462 if ( i != emt && i != rad
7463 && event[i].isFinal()
7464 && event[i].col() > 0
7465 && event[i].col() == radBeforeAcl)
7466 partners.push_back(i);
7471 for(
int i=0; i < int(event.size()); ++i) {
7473 if ( i != emt && i != rad
7474 && event[i].status() == -21
7475 &&
event[i].acol() > 0
7476 &&
event[i].acol() == radBeforeCol)
7477 partners.push_back(i);
7479 if ( i != emt && i != rad
7480 && event[i].isFinal()
7481 && event[i].col() > 0
7482 && event[i].col() == radBeforeCol)
7483 partners.push_back(i);
7485 if ( i != emt && i != rad
7486 && event[i].status() == -21
7487 && event[i].col() > 0
7488 && event[i].col() == radBeforeAcl)
7489 partners.push_back(i);
7491 if ( i != emt && i != rad
7492 && event[i].isFinal()
7493 && event[i].acol() > 0
7494 && event[i].acol() == radBeforeAcl)
7495 partners.push_back(i);
7522 bool History::getColSinglet(
const int flavType,
const int iParton,
7523 const Event& event, vector<int>& exclude, vector<int>& colSinglet) {
7526 if (iParton < 0)
return false;
7534 for(
int i=0; i < int(event.size()); ++i)
7535 if ( event[i].isFinal() &&
event[i].colType() != 0)
7540 int nExclude = int(exclude.size());
7541 int nInitExclude = 0;
7542 if (!event[exclude[2]].isFinal())
7544 if (!event[exclude[3]].isFinal())
7548 if (nFinal == nExclude - nInitExclude)
7558 colSinglet.push_back(iParton);
7560 exclude.push_back(iParton);
7563 colP = getColPartner(iParton,event);
7566 colP = getAcolPartner(iParton,event);
7569 for(
int i = 0; i < int(exclude.size()); ++i)
7570 if (colP == exclude[i])
7574 return getColSinglet(flavType,colP,event,exclude,colSinglet);
7585 bool History::isColSinglet(
const Event& event,
7586 vector<int> system ) {
7589 for(
int i=0; i < int(system.size()); ++i ) {
7592 && (event[system[i]].colType() == 1
7593 ||
event[system[i]].colType() == 2) ) {
7594 for(
int j=0; j < int(system.size()); ++j)
7597 && event[system[i]].col() ==
event[system[j]].acol()) {
7606 && (event[system[i]].colType() == -1
7607 || event[system[i]].colType() == 2) ) {
7608 for(
int j=0; j < int(system.size()); ++j)
7611 && event[system[i]].acol() ==
event[system[j]].col()) {
7623 bool isColSing =
true;
7624 for(
int i=0; i < int(system.size()); ++i)
7625 if ( system[i] != 0 )
7643 bool History::isFlavSinglet(
const Event& event,
7644 vector<int> system,
int flav) {
7649 for(
int i=0; i < int(system.size()); ++i)
7650 if ( system[i] > 0 ) {
7651 for(
int j=0; j < int(system.size()); ++j) {
7655 if ( event[i].idAbs() != 21
7656 &&
event[i].idAbs() != 22
7657 &&
event[i].idAbs() != 23
7658 &&
event[i].idAbs() != 24
7660 &&
event[system[i]].isFinal()
7661 &&
event[system[j]].isFinal()
7662 &&
event[system[i]].id() == -1*
event[system[j]].id()) {
7665 if (abs(flav) > 0 &&
event[system[i]].idAbs() != flav)
7675 if ( event[i].idAbs() != 21
7676 &&
event[i].idAbs() != 22
7677 &&
event[i].idAbs() != 23
7678 &&
event[i].idAbs() != 24
7680 &&
event[system[i]].isFinal() !=
event[system[j]].isFinal()
7681 &&
event[system[i]].id() ==
event[system[j]].id()) {
7684 if (abs(flav) > 0 &&
event[system[i]].idAbs() != flav)
7697 bool isFlavSing =
true;
7698 for(
int i=0; i < int(system.size()); ++i)
7699 if ( system[i] != 0 )
7714 bool History::allowedClustering(
int rad,
int emt,
int rec,
int partner,
7715 const Event& event ) {
7718 bool allowed =
true;
7723 bool isSing = isSinglett(rad,emt,partner,event);
7724 int type = (
event[rad].isFinal()) ? 1 :-1;
7726 int radBeforeFlav = getRadBeforeFlav(rad,emt,event);
7728 int radBeforeCol = getRadBeforeCol(rad,emt,event);
7729 int radBeforeAcl = getRadBeforeAcol(rad,emt,event);
7731 vector<int> radBeforeColP = getReclusteredPartners(rad, emt, event);
7734 int nPartonInHard = 0;
7735 for(
int i=0; i < int(event.size()); ++i)
7737 if ( event[i].isFinal()
7738 &&
event[i].colType() != 0
7739 && mergingHooksPtr->hardProcess->matchesAnyOutgoing(i, event) )
7745 for(
int i=0; i < int(event.size()); ++i)
7747 if ( i!=emt && i!=rad && i!=rec
7748 && event[i].isFinal()
7749 &&
event[i].colType() != 0
7750 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i, event) )
7754 int nInitialPartons = 0;
7755 for(
int i=0; i < int(event.size()); ++i)
7756 if ( event[i].status() == -21
7757 &&
event[i].colType() != 0 )
7762 for(
int i=0; i < int(event.size()); ++i)
7763 if ( event[i].isFinal()
7764 &&(
event[i].id() == 22
7765 ||
event[i].id() == 23
7766 ||
event[i].id() == 24
7767 ||(
event[i].idAbs() > 10 &&
event[i].idAbs() < 20)
7768 ||(event[i].idAbs() > 1000010 &&
event[i].idAbs() < 1000020)
7769 ||(event[i].idAbs() > 2000010 &&
event[i].idAbs() < 2000020) ))
7776 int nFinalQuark = 0;
7778 int nFinalQuarkExc = 0;
7779 for(
int i=0; i < int(event.size()); ++i) {
7780 if (i !=rad && i != emt && i != rec) {
7781 if (event[i].isFinal() && abs(event[i].colType()) == 1 ) {
7782 if ( !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,event) )
7791 if (event[rec].isFinal() &&
event[rec].isQuark()) nFinalQuark++;
7793 if ( event[rad].isFinal()
7794 && abs(particleDataPtr->colType(radBeforeFlav)) == 1) nFinalQuark++;
7797 int nInitialQuark = 0;
7799 if (event[rec].isFinal()) {
7800 if (event[3].isQuark()) nInitialQuark++;
7801 if (event[4].isQuark()) nInitialQuark++;
7803 int iOtherIn = (rec == 3) ? 4 : 3;
7804 if (event[rec].isQuark()) nInitialQuark++;
7805 if (event[iOtherIn].isQuark()) nInitialQuark++;
7809 if (event[rec].isQuark()) nInitialQuark++;
7811 if (abs(radBeforeFlav) < 10) nInitialQuark++;
7817 int proton[] = {1,2,3,4,5,21,22,23,24};
7818 bool isInProton =
false;
7819 for(
int i=0; i < 9; ++i)
7820 if (abs(radBeforeFlav) == proton[i]) isInProton =
true;
7821 if (type == -1 && !isInProton)
return false;
7824 vector<int> unmatchedCol;
7825 vector<int> unmatchedAcl;
7827 for (
int i = 0; i <
event.size(); ++i)
7828 if ( i != emt && i != rad
7829 && (event[i].isFinal() || event[i].status() == -21)
7830 &&
event[i].colType() != 0 ) {
7832 int colP = getColPartner(i,event);
7833 int aclP = getAcolPartner(i,event);
7835 if (event[i].col() > 0
7836 && (colP == emt || colP == rad || colP == 0) )
7837 unmatchedCol.push_back(i);
7838 if (event[i].acol() > 0
7839 && (aclP == emt || aclP == rad || aclP == 0) )
7840 unmatchedAcl.push_back(i);
7846 if (
int(unmatchedCol.size()) +
int(unmatchedAcl.size()) > 2)
7853 if ( isSing && event[rec].isQuark()
7854 && abs(particleDataPtr->colType(radBeforeFlav)) == 1)
7858 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(emt,event) ) {
7862 bool canReplace = mergingHooksPtr->hardProcess->findOtherCandidates(emt,
7864 if (canReplace) allowed =
true;
7865 else allowed =
false;
7870 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(rad,event)
7871 &&
event[rad].id() != radBeforeFlav )
7876 if (nFinalEW != 0 && nInitialQuark == 0
7877 && nFinalQuark == 0 && nFinalQuarkExc == 0)
7880 if ( (nInitialQuark + nFinalQuark + nFinalQuarkExc)%2 != 0 )
7890 for(
int i=0; i < int(event.size()); ++i)
7891 if ( i!=emt && i!=rad && i!=rec
7892 && (event[i].mother1() == 1 ||
event[i].mother1() == 2))
7893 in.push_back(event[i].
id());
7894 if (!event[rad].isFinal()) in.push_back(radBeforeFlav);
7895 if (!event[rec].isFinal()) in.push_back(event[rec].
id());
7897 for(
int i=0; i < int(event.size()); ++i)
7898 if ( i!=emt && i!=rad && i!=rec && event[i].isFinal())
7899 out.push_back(event[i].id());
7900 if (event[rad].isFinal()) out.push_back(radBeforeFlav);
7901 if (event[rec].isFinal()) out.push_back(event[rec].
id());
7902 if (event[3].col() == event[4].acol()
7903 && event[3].acol() == event[4].col()
7904 && !mergingHooksPtr->allowEffectiveVertex( in, out)
7905 && nFinalQuark == 0){
7908 int nTripletts = abs(event[rec].colType())
7909 + abs(particleDataPtr->colType(radBeforeFlav));
7910 if (event[3].isGluon()) allowed =
false;
7911 else if (nTripletts != 2 && nFinalQuarkExc%2 == 0) allowed =
false;
7915 if ( abs((event[rad].p()+type*event[emt].p()+event[rec].p()).pz())
7916 > (
event[rad].p()+type*
event[emt].p()+
event[rec].p()).e()
7918 && (
event[rad].p()-
event[emt].p()+
event[rec].p()).m2Calc() < 0.) ){
7923 if (event[emt].
id() == 21)
return allowed;
7926 if (event[emt].
id() == 1000021)
return allowed;
7929 vector<int> outgoingParticles;
7930 int nOut1 = int(mergingHooksPtr->hardProcess->PosOutgoing1.size());
7931 for (
int i=0; i < nOut1; ++i ) {
7932 int iPos = mergingHooksPtr->hardProcess->PosOutgoing1[i];
7933 outgoingParticles.push_back(
7934 mergingHooksPtr->hardProcess->state[iPos].id() );
7936 int nOut2 = int(mergingHooksPtr->hardProcess->PosOutgoing2.size());
7937 for (
int i=0; i < nOut2; ++i ) {
7938 int iPos = mergingHooksPtr->hardProcess->PosOutgoing2[i];
7939 outgoingParticles.push_back(
7940 mergingHooksPtr->hardProcess->state[iPos].id() );
7947 vector<int> iInQuarkFlav;
7948 for(
int i=0; i < int(event.size()); ++i)
7950 if ( i != emt && i != rad
7951 && event[i].status() == -21
7952 &&
event[i].idAbs() ==
event[emt].idAbs() )
7953 iInQuarkFlav.push_back(i);
7956 vector<int> iOutQuarkFlav;
7957 for(
int i=0; i < int(event.size()); ++i)
7959 if ( i != emt && i != rad
7960 && event[i].isFinal()
7961 &&
event[i].idAbs() ==
event[emt].idAbs() ) {
7965 bool matchOut =
false;
7966 for (
int j = 0; j < int(outgoingParticles.size()); ++j)
7967 if ( event[i].idAbs() == abs(outgoingParticles[j])) {
7969 outgoingParticles[j] = 99;
7971 if (!matchOut) iOutQuarkFlav.push_back(i);
7976 int nInQuarkFlav = int(iInQuarkFlav.size());
7977 int nOutQuarkFlav = int(iOutQuarkFlav.size());
7982 if ( event[partner].isFinal()
7983 &&
event[partner].id() == 21
7984 && radBeforeFlav == 21
7985 &&
event[partner].col() == radBeforeCol
7986 &&
event[partner].acol() == radBeforeAcl)
7990 if (nInQuarkFlav + nOutQuarkFlav == 0)
7997 vector<int> partons;
7998 for(
int i=0; i < int(event.size()); ++i)
8000 if ( i!=emt && i!=rad
8001 && event[i].colType() != 0
8002 && (
event[i].isFinal() ||
event[i].status() == -21) ) {
8004 partons.push_back(i);
8006 if (event[i].colType() == 2)
8008 else if (event[i].colType() == 1)
8010 else if (event[i].colType() == -1)
8016 bool isFSRg2qq = ((type == 1) && (event[rad].
id() == -1*
event[emt].id()) );
8017 bool isISRg2qq = ((type ==-1) && (event[rad].
id() ==
event[emt].id()) );
8022 if ( (isFSRg2qq || isISRg2qq)
8023 && int(quark.size()) +
int(antiq.size())
8024 +
int(gluon.size()) > nPartonInHard ) {
8026 vector<int> colours;
8027 vector<int> anticolours;
8032 colours.push_back(radBeforeCol);
8033 anticolours.push_back(radBeforeAcl);
8035 colours.push_back(radBeforeAcl);
8036 anticolours.push_back(radBeforeCol);
8039 for(
int i=0; i < int(gluon.size()); ++i)
8040 if (event[gluon[i]].isFinal()) {
8041 colours.push_back(event[gluon[i]].col());
8042 anticolours.push_back(event[gluon[i]].acol());
8044 colours.push_back(event[gluon[i]].acol());
8045 anticolours.push_back(event[gluon[i]].col());
8050 for(
int i=0; i < int(colours.size()); ++i)
8051 for(
int j=0; j < int(anticolours.size()); ++j)
8052 if (colours[i] > 0 && anticolours[j] > 0
8053 && colours[i] == anticolours[j]) {
8061 bool allMatched =
true;
8062 for(
int i=0; i < int(colours.size()); ++i)
8063 if (colours[i] != 0)
8065 for(
int i=0; i < int(anticolours.size()); ++i)
8066 if (anticolours[i] != 0)
8074 for(
int i=0; i < int(quark.size()); ++i)
8075 if ( event[quark[i]].isFinal()
8076 && mergingHooksPtr->hardProcess->matchesAnyOutgoing(quark[i], event) )
8077 colours.push_back(event[quark[i]].col());
8079 for(
int i=0; i < int(antiq.size()); ++i)
8080 if ( event[antiq[i]].isFinal()
8081 && mergingHooksPtr->hardProcess->matchesAnyOutgoing(antiq[i], event) )
8082 anticolours.push_back(event[antiq[i]].acol());
8086 for(
int i=0; i < int(colours.size()); ++i)
8088 for(
int j=0; j < int(anticolours.size()); ++j)
8089 if (colours[i] > 0 && anticolours[j] > 0
8090 && colours[i] == anticolours[j]) {
8097 for (
int i=0; i < int(quark.size()); ++i )
8098 if ( !mergingHooksPtr->hardProcess->matchesAnyOutgoing( quark[i],
8101 for (
int i=0; i < int(antiq.size()); ++i )
8102 if ( !mergingHooksPtr->hardProcess->matchesAnyOutgoing( antiq[i],
8105 for(
int i=0; i < int(gluon.size()); ++i)
8106 if ( event[gluon[i]].isFinal() )
8114 for(
int i=0; i < int(colours.size()); ++i)
8115 if (colours[i] != 0)
8117 for(
int i=0; i < int(anticolours.size()); ++i)
8118 if (anticolours[i] != 0)
8121 if (allMatched && nNotInHard > 0)
8128 if (isFSRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) {
8132 for(
int i=0; i < int(gluon.size()); ++i) {
8133 if (!event[gluon[i]].isFinal()
8134 &&
event[gluon[i]].col() == radBeforeCol
8135 &&
event[gluon[i]].acol() == radBeforeAcl)
8141 for(
int i=0; i < int(gluon.size()); ++i) {
8142 if (event[gluon[i]].isFinal()
8143 &&
event[gluon[i]].col() == radBeforeAcl
8144 &&
event[gluon[i]].acol() == radBeforeCol)
8150 if (
int(radBeforeColP.size()) == 2
8151 && event[radBeforeColP[0]].isFinal()
8152 &&
event[radBeforeColP[1]].isFinal()
8153 &&
event[radBeforeColP[0]].id() == -1*
event[radBeforeColP[1]].id() ) {
8157 if (nInitialPartons > 0)
8164 if (
int(radBeforeColP.size()) == 2
8165 && (( event[radBeforeColP[0]].status() == -21
8166 && event[radBeforeColP[1]].isFinal())
8167 ||(
event[radBeforeColP[0]].isFinal()
8168 &&
event[radBeforeColP[1]].status() == -21))
8169 &&
event[radBeforeColP[0]].id() ==
event[radBeforeColP[1]].id() ) {
8176 int incoming = (
event[radBeforeColP[0]].isFinal())
8177 ? radBeforeColP[1] : radBeforeColP[0];
8178 int outgoing = (
event[radBeforeColP[0]].isFinal())
8179 ? radBeforeColP[0] : radBeforeColP[1];
8182 bool clusPossible =
false;
8183 for(
int i=0; i < int(event.size()); ++i)
8184 if ( i != emt && i != rad
8185 && i != incoming && i != outgoing
8186 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,event) ) {
8188 if ( event[i].status() == -21
8189 && (
event[i].id() ==
event[outgoing].id()
8190 ||
event[i].id() == -1*
event[incoming].id()) )
8191 clusPossible =
true;
8193 if ( event[i].isFinal()
8194 && (
event[i].id() == -1*
event[outgoing].id()
8195 ||
event[i].id() ==
event[incoming].id()) )
8196 clusPossible =
true;
8207 vector<int> excludeIn1;
8208 for(
int i=0; i < 4; ++i)
8209 excludeIn1.push_back(0);
8210 vector<int> colSingletIn1;
8211 int flavIn1Type = (
event[incoming].id() > 0) ? 1 : -1;
8213 bool isColSingIn1 = getColSinglet(flavIn1Type,incoming,event,
8214 excludeIn1,colSingletIn1);
8216 bool isFlavSingIn1 = isFlavSinglet(event,colSingletIn1);
8221 int incoming2 = (incoming == 3) ? 4 : 3;
8222 vector<int> excludeIn2;
8223 for(
int i=0; i < 4; ++i)
8224 excludeIn2.push_back(0);
8225 vector<int> colSingletIn2;
8226 int flavIn2Type = (
event[incoming2].id() > 0) ? 1 : -1;
8228 bool isColSingIn2 = getColSinglet(flavIn2Type,incoming2,event,
8229 excludeIn2,colSingletIn2);
8231 bool isFlavSingIn2 = isFlavSinglet(event,colSingletIn2);
8235 && (!isColSingIn1 || !isFlavSingIn1
8236 || !isColSingIn2 || !isFlavSingIn2))
8248 int flav =
event[emt].id();
8249 vector<int> exclude;
8250 exclude.push_back(emt);
8251 exclude.push_back(rad);
8252 exclude.push_back(radBeforeColP[0]);
8253 exclude.push_back(radBeforeColP[1]);
8254 vector<int> colSinglet;
8258 for(
int i=0; i < int(event.size()); ++i)
8263 && i != radBeforeColP[0]
8264 && i != radBeforeColP[1]
8265 && event[i].isFinal() ) {
8267 if (event[i].
id() == flav) {
8273 int flavType = (iOther > 0 &&
event[iOther].id() > 0) ? 1
8274 : (iOther > 0) ? -1 : 0;
8276 bool isColSing = getColSinglet(flavType,iOther,event,exclude,colSinglet);
8278 bool isFlavSing = isFlavSinglet(event,colSinglet);
8282 bool isHardSys =
true;
8283 for(
int i=0; i < int(colSinglet.size()); ++i)
8284 isHardSys = mergingHooksPtr->hardProcess->matchesAnyOutgoing(
8285 colSinglet[i], event);
8290 if (isColSing && isFlavSing && !isHardSys) {
8296 vector<int> allFinal;
8297 for(
int i=0; i < int(event.size()); ++i)
8298 if ( event[i].isFinal() )
8299 allFinal.push_back(i);
8302 bool isFullColSing = isColSinglet(event,allFinal);
8304 bool isFullFlavSing = isFlavSinglet(event,allFinal,flav);
8309 if (!isFullColSing || !isFullFlavSing)
8316 if (isISRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) {
8320 for(
int i=0; i < int(gluon.size()); ++i) {
8321 if (event[gluon[i]].isFinal()
8322 &&
event[gluon[i]].col() == radBeforeCol
8323 &&
event[gluon[i]].acol() == radBeforeAcl)
8329 for(
int i=0; i < int(gluon.size()); ++i) {
8330 if (event[gluon[i]].status() == -21
8331 &&
event[gluon[i]].acol() == radBeforeCol
8332 &&
event[gluon[i]].col() == radBeforeAcl)
8338 if (
int(radBeforeColP.size()) == 2
8339 && event[radBeforeColP[0]].isFinal()
8340 &&
event[radBeforeColP[1]].isFinal()
8341 &&
event[radBeforeColP[0]].id() == -1*
event[radBeforeColP[1]].id() ) {
8348 bool clusPossible =
false;
8349 for(
int i=0; i < int(event.size()); ++i)
8350 if ( i != emt && i != rad
8351 && i != radBeforeColP[0]
8352 && i != radBeforeColP[1]
8353 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,event) ) {
8354 if (event[i].status() == -21
8355 && (
event[radBeforeColP[0]].id() ==
event[i].id()
8356 ||
event[radBeforeColP[1]].id() ==
event[i].id() ))
8358 clusPossible =
true;
8359 if (event[i].isFinal()
8360 && (
event[radBeforeColP[0]].id() == -1*
event[i].id()
8361 ||
event[radBeforeColP[1]].id() == -1*
event[i].id() ))
8362 clusPossible =
true;
8374 vector<int> excludeIn1;
8375 for(
int i=0; i < 4; ++i)
8376 excludeIn1.push_back(0);
8377 vector<int> colSingletIn1;
8378 int flavIn1Type = (
event[incoming1].id() > 0) ? 1 : -1;
8380 bool isColSingIn1 = getColSinglet(flavIn1Type,incoming1,event,
8381 excludeIn1,colSingletIn1);
8383 bool isFlavSingIn1 = isFlavSinglet(event,colSingletIn1);
8389 vector<int> excludeIn2;
8390 for(
int i=0; i < 4; ++i)
8391 excludeIn2.push_back(0);
8392 vector<int> colSingletIn2;
8393 int flavIn2Type = (
event[incoming2].id() > 0) ? 1 : -1;
8395 bool isColSingIn2 = getColSinglet(flavIn2Type,incoming2,event,
8396 excludeIn2,colSingletIn2);
8398 bool isFlavSingIn2 = isFlavSinglet(event,colSingletIn2);
8402 && (!isColSingIn1 || !isFlavSingIn1
8403 || !isColSingIn2 || !isFlavSingIn2))
8422 bool History::isSinglett(
int rad,
int emt,
int rec,
const Event& event ) {
8424 int radCol =
event[rad].col();
8425 int emtCol =
event[emt].col();
8426 int recCol =
event[rec].col();
8427 int radAcl =
event[rad].acol();
8428 int emtAcl =
event[emt].acol();
8429 int recAcl =
event[rec].acol();
8430 int recType =
event[rec].isFinal() ? 1 : -1;
8432 bool isSing =
false;
8434 if ( ( recType == -1
8435 && radCol + emtCol == recCol && radAcl + emtAcl == recAcl)
8437 && radCol + emtCol == recAcl && radAcl + emtAcl == recCol) )
8453 bool History::validEvent(
const Event& event ) {
8456 bool validColour =
true;
8457 for (
int i = 0; i <
event.size(); ++i)
8459 if ( event[i].isFinal() &&
event[i].colType() == 1
8461 && ( FindCol(event[i].col(),i,0,event,1,
true) == 0
8463 && FindCol(event[i].col(),i,0,event,2,
true) == 0 )) {
8464 validColour =
false;
8467 }
else if ( event[i].isFinal() &&
event[i].colType() == -1
8469 && ( FindCol(event[i].acol(),i,0,event,2,
true) == 0
8471 && FindCol(event[i].acol(),i,0,event,1,
true) == 0 )) {
8472 validColour =
false;
8475 }
else if ( event[i].isFinal() &&
event[i].colType() == 2
8477 && ( FindCol(event[i].col(),i,0,event,1,
true) == 0
8479 && FindCol(event[i].col(),i,0,event,2,
true) == 0 )
8481 && ( FindCol(event[i].acol(),i,0,event,2,
true) == 0
8483 && FindCol(event[i].acol(),i,0,event,1,
true) == 0 )) {
8484 validColour =
false;
8489 bool validCharge =
true;
8490 double initCharge =
event[3].charge() +
event[4].charge();
8491 double finalCharge = 0.0;
8492 for(
int i = 0; i <
event.size(); ++i)
8493 if (event[i].isFinal()) finalCharge += event[i].charge();
8494 if (abs(initCharge-finalCharge) > 1e-12) validCharge =
false;
8496 return (validColour && validCharge);
8505 bool History::equalClustering( Clustering clus1 , Clustering clus2 ) {
8506 return ( (clus1.emittor == clus2.emittor)
8507 && (clus1.emitted == clus2.emitted)
8508 && (clus1.recoiler == clus2.recoiler)
8509 && (clus1.partner == clus2.partner)
8510 && (clus1.pT() == clus2.pT())
8511 && (clus1.spinRadBef == clus2.spinRadBef)
8512 && (clus1.spinRad == clus2.spinRad)
8513 && (clus1.spinEmt == clus2.spinEmt)
8514 && (clus1.spinRec == clus2.spinRec)
8515 && (clus1.flavRadBef == clus2.flavRadBef));
8524 double History::choseHardScale(
const Event& event )
const {
8527 double mHat = (
event[3].p() +
event[4].p()).mCalc();
8534 for(
int i = 0; i <
event.size(); ++i)
8535 if ( event[i].isFinal() ) {
8538 if ( event[i].idAbs() == 23
8539 ||
event[i].idAbs() == 24 ) {
8542 mBos +=
event[i].m();
8544 }
else if ( abs(event[i].status()) == 22
8545 && ( event[i].idAbs() == 23
8546 || event[i].idAbs() == 24 )) {
8548 mBos +=
event[i].m();
8552 if ( nBosons > 0 && (nFinal + nFinBos*2) <= 3)
8553 return (mBos /
double(nBosons));
8564 int History::getCurrentFlav(
const int side)
const {
8565 int in = (side == 1) ? 3 : 4;
8566 return state[in].id();
8571 double History::getCurrentX(
const int side)
const {
8572 int in = (side == 1) ? 3 : 4;
8573 return ( 2.*state[in].e()/state[0].e() );
8578 double History::getCurrentZ(
const int rad,
8579 const int rec,
const int emt,
int idRadBef)
const {
8581 int type = state[rad].isFinal() ? 1 : -1;
8586 Vec4 radAfterBranch(state[rad].p());
8587 Vec4 recAfterBranch(state[rec].p());
8588 Vec4 emtAfterBranch(state[emt].p());
8591 double m2RadAft = radAfterBranch.m2Calc();
8592 double m2EmtAft = emtAfterBranch.m2Calc();
8593 double m2RadBef = 0.;
8594 if ( state[rad].idAbs() != 21 && state[rad].idAbs() != 22
8595 && state[emt].idAbs() != 24 && state[rad].idAbs() != state[emt].idAbs())
8596 m2RadBef = m2RadAft;
8597 else if ( state[emt].idAbs() == 24) {
8599 m2RadBef = pow2(particleDataPtr->m0(abs(idRadBef)));
8602 double Qsq = (radAfterBranch + emtAfterBranch).m2Calc();
8606 = (radAfterBranch + recAfterBranch + emtAfterBranch).m2Calc();
8608 if ( !state[rec].isFinal() ){
8609 double mar2 = m2final - 2. * Qsq + 2. * m2RadBef;
8610 recAfterBranch *= (1. - (Qsq - m2RadBef)/(mar2 - m2RadBef))
8611 /(1. + (Qsq - m2RadBef)/(mar2 - m2RadBef));
8614 if (Qsq > mar2)
return 0.5;
8617 Vec4 sum = radAfterBranch + recAfterBranch + emtAfterBranch;
8618 double m2Dip = sum.m2Calc();
8620 double x1 = 2. * (sum * radAfterBranch) / m2Dip;
8621 double x2 = 2. * (sum * recAfterBranch) / m2Dip;
8624 double lambda13 = sqrt( pow2(Qsq - m2RadAft - m2EmtAft )
8625 - 4.*m2RadAft*m2EmtAft);
8626 double k1 = ( Qsq - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8627 double k3 = ( Qsq - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8629 z = 1./ ( 1- k1 -k3) * ( x1 / (2.-x2) - k3);
8633 Vec4 qBR(state[rad].p() - state[emt].p() + state[rec].p());
8634 Vec4 qAR(state[rad].p() + state[rec].p());
8636 z = (qBR.m2Calc())/( qAR.m2Calc());
8647 double History::pTLund(
const Event& event,
int rad,
int emt,
int rec,
8648 int ShowerType,
int idRadBef) {
8651 Particle RadAfterBranch =
event[rad];
8652 Particle EmtAfterBranch =
event[emt];
8653 Particle RecAfterBranch =
event[rec];
8656 if ( mergingHooksPtr->useShowerPlugin() ) {
8657 map<string,double> stateVars;
8658 bool isFSR = showers->timesPtr->isTimelike(event, rad, emt, rec,
"");
8660 string name = showers->timesPtr->getSplittingName(event, rad, emt,
8662 stateVars = showers->timesPtr->getStateVariables(event, rad, emt, rec,
8665 string name = showers->spacePtr->getSplittingName(event, rad, emt,
8667 stateVars = showers->spacePtr->getStateVariables(event, rad, emt, rec,
8671 return ( (stateVars.size() > 0 && stateVars.find(
"t") != stateVars.end())
8672 ? sqrt(stateVars[
"t"]) : -1.0 );
8676 int Type = ShowerType;
8678 int sign = (Type==1) ? 1 : -1;
8679 Vec4 Q(RadAfterBranch.p() + sign*EmtAfterBranch.p());
8680 double Qsq = sign * Q.m2Calc();
8683 Vec4 radAft(RadAfterBranch.p());
8684 Vec4 recAft(RecAfterBranch.p());
8685 Vec4 emtAft(EmtAfterBranch.p());
8688 double m2RadAft = radAft.m2Calc();
8689 double m2EmtAft = emtAft.m2Calc();
8691 double m2RadBef = 0.;
8692 if ( RadAfterBranch.idAbs() != 21 && RadAfterBranch.idAbs() != 22
8693 && EmtAfterBranch.idAbs() != 24
8694 && RadAfterBranch.idAbs() != EmtAfterBranch.idAbs() )
8695 m2RadBef = m2RadAft;
8696 else if (EmtAfterBranch.idAbs() == 24) {
8698 m2RadBef = pow2(particleDataPtr->m0(abs(idRadBef)));
8699 }
else if (!RadAfterBranch.isFinal()) {
8700 if (RadAfterBranch.idAbs() == 21 && EmtAfterBranch.idAbs() != 21)
8701 m2RadBef = m2EmtAft;
8705 double m2final = (radAft + recAft + emtAft).m2Calc();
8707 if ( !RecAfterBranch.isFinal() && RadAfterBranch.isFinal() ){
8708 double mar2 = m2final - 2. * Qsq + 2. * m2RadBef;
8709 recAft *= (1. - (Qsq - m2RadBef)/(mar2 - m2RadBef))
8710 /(1. + (Qsq - m2RadBef)/(mar2 - m2RadBef));
8712 if (Qsq > mar2)
return 0.;
8715 Vec4 sum = radAft + recAft + emtAft;
8716 double m2Dip = sum.m2Calc();
8717 double x1 = 2. * (sum * radAft) / m2Dip;
8718 double x2 = 2. * (sum * recAft) / m2Dip;
8721 double q2BR = (RadAfterBranch.p() - EmtAfterBranch.p()
8722 + RecAfterBranch.p()).m2Calc();
8723 double q2AR = (RadAfterBranch.p() + RecAfterBranch.p()).m2Calc();
8726 double lambda13 = sqrt( pow2(Qsq - m2RadAft - m2EmtAft )
8727 - 4. * m2RadAft*m2EmtAft );
8728 double k1 = ( Qsq - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8729 double k3 = ( Qsq - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8732 double z = (Type==1) ? 1./ ( 1- k1 -k3) * ( x1 / (2.-x2) - k3)
8736 double pTpyth = (Type==1) ? z*(1.-z) : (1.-z);
8739 if (Type == 1) pTpyth *= (Qsq - m2RadBef);
8745 if ( (RadAfterBranch.idAbs() == 4 || EmtAfterBranch.idAbs() == 4)
8746 && RadAfterBranch.idAbs() != EmtAfterBranch.idAbs()) {
8747 if (pTpyth < 2 * pow2(particleDataPtr->m0(4)))
8748 pTpyth = (Qsq + pow2(particleDataPtr->m0(4)) ) * (1. - q2BR/q2AR);
8749 }
else if ( (RadAfterBranch.idAbs() == 5 || EmtAfterBranch.idAbs() == 5)
8750 && RadAfterBranch.idAbs() != EmtAfterBranch.idAbs()) {
8751 if (pTpyth < 2 * pow2(particleDataPtr->m0(5)))
8752 pTpyth = (Qsq + pow2(particleDataPtr->m0(5)) ) * (1. - q2BR/q2AR);
8756 if ( pTpyth < 0. ) pTpyth = 0.;
8759 return sqrt(pTpyth);
8768 int History::posChangedIncoming(
const Event& event,
bool before) {
8774 for (
int i =0; i <
event.size(); ++i)
8775 if (event[i].status() == 43) {
8781 if (iSister > 0) iMother =
event[iSister].mother1();
8784 if (iSister > 0 && iMother > 0) {
8787 int flavSister =
event[iSister].id();
8788 int flavMother =
event[iMother].id();
8791 int flavDaughter = 0;
8792 if ( abs(flavMother) < 21 && flavSister == 21)
8793 flavDaughter = flavMother;
8794 else if ( flavMother == 21 && flavSister == 21)
8795 flavDaughter = flavMother;
8796 else if ( flavMother == 21 && abs(flavSister) < 21)
8797 flavDaughter = -1*flavSister;
8798 else if ( abs(flavMother) < 21 && abs(flavSister) < 21)
8803 for (
int i =0; i <
event.size(); ++i)
8804 if ( !event[i].isFinal()
8805 &&
event[i].mother1() == iMother
8806 &&
event[i].id() == flavDaughter )
8810 if ( !before )
return iMother;
8811 else return iDaughter;
8819 for (
int i =0; i <
event.size(); ++i)
8820 if ( abs(event[i].status()) == 53 || abs(event[i].status()) == 54) {
8826 if (iMother > 0) iDaughter =
event[iMother].daughter1();
8829 if (iDaughter > 0 && iMother > 0) {
8832 if ( !before )
return iMother;
8833 else return iDaughter;
8849 double History::pdfFactor(
const Event& event,
const int type,
8850 double pdfScale,
double mu ) {
8859 for (
int i =0; i <
event.size(); ++i)
8860 if ( abs(event[i].status()) == 53 || abs(event[i].status()) == 54) {
8864 int flavMother =
event[iMother].id();
8867 if ( iMother == 0 )
return 1.;
8870 int iDaughter =
event[iMother].daughter1();
8871 int flavDaughter =
event[iDaughter].id();
8874 double xMother = 2.*
event[iMother].e() /
event[0].e();
8875 double xDaughter = 2.*
event[iDaughter].e() /
event[0].e();
8879 int sideSplit = (
event[iMother].pz() > 0.) ? 1 : -1;
8880 double pdfDen1, pdfDen2, pdfNum1, pdfNum2;
8881 pdfDen1 = pdfDen2 = pdfNum1 = pdfNum2 = 1.;
8882 if ( sideSplit == 1 ) {
8884 pdfDen1 = max(1e-15,beamA.xfISR(0, flavDaughter, xDaughter, pow2(mu)) );
8885 pdfNum1 = beamA.xfISR(0, flavDaughter, xDaughter, pow2(pdfScale) );
8886 pdfNum2 = beamA.xfISR(0, flavMother, xMother, pow2(mu) );
8887 pdfDen2 = max(1e-15,beamA.xfISR(0,flavMother, xMother, pow2(pdfScale)) );
8890 pdfDen1 = max(1e-15,beamB.xfISR(0, flavDaughter, xDaughter, pow2(mu)) );
8891 pdfNum1 = beamB.xfISR(0, flavDaughter, xDaughter, pow2(pdfScale) );
8892 pdfNum2 = beamB.xfISR(0, flavMother, xMother, pow2(mu) );
8893 pdfDen2 = max(1e-15,beamB.xfISR(0,flavMother, xMother, pow2(pdfScale)) );
8898 if ( pdfDen2/pdfNum1 > 1. )
return 1.;
8903 weight = (pdfNum1/pdfDen1) * (pdfNum2)/(pdfDen2);
8906 }
else if (type == 2) {
8910 for (
int i =0; i <
event.size(); ++i)
8911 if (event[i].status() == 43) {
8915 int flavSister =
event[iSister].id();
8918 int iMother =
event[iSister].mother1();
8919 int flavMother =
event[iMother].id();
8922 int flavDaughter = 0;
8923 if ( abs(flavMother) < 21 && flavSister == 21)
8924 flavDaughter = flavMother;
8925 else if ( flavMother == 21 && flavSister == 21)
8926 flavDaughter = flavMother;
8927 else if ( flavMother == 21 && abs(flavSister) < 21)
8928 flavDaughter = -1*flavSister;
8929 else if ( abs(flavMother) < 21 && abs(flavSister) < 21)
8933 double xMother = 2.*
event[iMother].e() /
event[0].e();
8937 for (
int i =0; i <
event.size(); ++i)
8938 if ( !event[i].isFinal()
8939 &&
event[i].mother1() == iMother
8940 &&
event[i].id() == flavDaughter )
8942 double xDaughter = 2.*
event[iDaughter].e() /
event[0].e();
8947 int sideSplit = (
event[iMother].pz() > 0.) ? 1 : -1;
8948 double ratio1 = getPDFratio( sideSplit,
false,
false, flavDaughter,
8949 xDaughter, pdfScale, flavDaughter, xDaughter, mu );
8950 double ratio2 = getPDFratio( sideSplit,
false,
false, flavMother,
8951 xMother, mu, flavMother, xMother, pdfScale );
8953 weight = ratio1*ratio2;
8970 double History::integrand(
int flav,
double x,
double scaleInt,
double z) {
8982 AlphaStrong* as = mergingHooksPtr->AlphaS_ISR();
8983 double asNow = (*as).alphaS(z);
8984 result = 1./z *asNow*asNow* ( log(scaleInt/z) -3./2. );
8988 }
else if (flav==21) {
8990 double measure1 = 1./(1. - z);
8991 double measure2 = 1.;
8995 * z * getPDFratio( 2,
false,
true, 21, x/z, scaleInt, 21, x, scaleInt )
9000 2.*CA *((1. -z)/z + z*(1.-z))
9001 * getPDFratio( 2,
false,
true, 21, x/z, scaleInt, 21, x, scaleInt )
9003 + CF * ((1+pow(1-z,2))/z)
9004 *( getPDFratio(2,
false,
true,1,x/z,scaleInt,21,x,scaleInt)
9005 + getPDFratio(2,
false,
true,-1,x/z,scaleInt,21,x,scaleInt)
9006 + getPDFratio(2,
false,
true,2,x/z,scaleInt,21,x,scaleInt)
9007 + getPDFratio(2,
false,
true,-2,x/z,scaleInt,21,x,scaleInt)
9008 + getPDFratio(2,
false,
true,3,x/z,scaleInt,21,x,scaleInt)
9009 + getPDFratio(2,
false,
true,-3,x/z,scaleInt,21,x,scaleInt)
9010 + getPDFratio(2,
false,
true,4,x/z,scaleInt,21,x,scaleInt)
9011 + getPDFratio(2,
false,
true,-4,x/z,scaleInt,21,x,scaleInt) );
9014 result = integrand1*measure1 + integrand2*measure2;
9018 double measure1 = 1./(1. -z);
9019 double measure2 = 1.;
9024 * getPDFratio(2,
false,
true,flav,x/z,scaleInt,flav,x,scaleInt)
9029 + TR * (pow(z,2) + pow(1-z,2))
9030 * getPDFratio(2,
false,
true,21,x/z,scaleInt,flav,x,scaleInt);
9033 result = measure1*integrand1 + measure2*integrand2;
9045 vector<int> History::posFlavCKM(
int flav) {
9048 int flavAbs = abs(flav);
9050 vector<int> flavRadBefs;
9052 if (flavAbs > 10 && flavAbs % 2 == 1)
9053 flavRadBefs.push_back(flavAbs + 1);
9055 else if (flavAbs > 10 && flavAbs % 2 == 0)
9056 flavRadBefs.push_back(flavAbs - 1);
9058 else if (flavAbs < 10 && flavAbs % 2 == 1) {
9059 flavRadBefs.push_back(2);
9060 flavRadBefs.push_back(4);
9061 flavRadBefs.push_back(6);
9063 else if (flavAbs < 10 && flavAbs % 2 == 0) {
9064 flavRadBefs.push_back(1);
9065 flavRadBefs.push_back(3);
9066 flavRadBefs.push_back(5);
9079 bool History::checkFlavour(vector<int>& flavCounts,
int flavRad,
9080 int flavRadBef,
int clusType) {
9083 for(
int k = 0; k < 20; ++k) {
9086 if (abs(flavRad) == k) {
9088 if (flavRad < 0) cor = 1;
9091 if (abs(flavRadBef) == k) {
9093 if (flavRadBef < 0) cor = -1;
9097 if (flavRadBef == flavRad) cor = 0;
9100 if (clusType == 1) {
9101 if (flavCounts[k] + cor != 0)
return false;
9104 if (flavCounts[k] - cor != 0)
return false;
9117 void History::reverseBoostISR(Vec4& pMother, Vec4& pSister, Vec4& pPartner,
9118 Vec4& pDaughter, Vec4& pRecoiler,
int sign,
double eCM,
double& phi ) {
9122 phi = pSister.phi();
9124 RotBstMatrix rot_by_mphi;
9125 rot_by_mphi.rot(0.,-phi);
9127 RotBstMatrix rot_by_pphi;
9128 rot_by_pphi.rot(0.,phi);
9132 double x1 = 2. * pMother.e() / eCM;
9134 double x2 = 2. * pPartner.e() / eCM;
9137 Vec4 qDip( pMother - pSister);
9138 Vec4 qAfter(pMother + pPartner);
9139 Vec4 qBefore(qDip + pPartner);
9140 double z = qBefore.m2Calc() / qAfter.m2Calc();
9143 double x1New = z*x1;
9145 double sHat = x1New*x2New*eCM*eCM;
9150 Vec4 pDaughterBef( 0., 0., sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
9151 Vec4 pRecoilerBef( 0., 0., -sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
9154 pMother.rotbst( rot_by_mphi );
9155 pSister.rotbst( rot_by_mphi );
9156 pPartner.rotbst( rot_by_mphi );
9160 pDaughter.p( pMother - pSister);
9161 pRecoiler.p( pPartner );
9162 RotBstMatrix from_CM_to_DRoff;
9164 from_CM_to_DRoff.toCMframe(pDaughter, pRecoiler);
9166 from_CM_to_DRoff.toCMframe(pRecoiler, pDaughter);
9170 pMother.rotbst( from_CM_to_DRoff );
9171 pPartner.rotbst( from_CM_to_DRoff );
9172 pSister.rotbst( from_CM_to_DRoff );
9177 RotBstMatrix from_DR_to_CM;
9178 from_DR_to_CM.bst( 0., 0., sign*( x1New - x2New ) / ( x1New + x2New ) );
9183 pDaughterBef.rotbst( from_DR_to_CM );
9184 pRecoilerBef.rotbst( from_DR_to_CM );
9188 if ( abs(pRecoilerBef.mCalc()) > 1e-7 ) {
9189 double pzSign = (pRecoilerBef.pz() > 0.) ? 1. : -1.;
9190 double eRec = pRecoilerBef.e();
9191 pRecoilerBef.p(0., 0., pzSign*eRec, eRec);
9193 if ( abs(pDaughterBef.mCalc()) > 1e-7 ) {
9194 double pzSign = (pDaughterBef.pz() > 0.) ? 1. : -1.;
9195 double eDau = pDaughterBef.e();
9196 pDaughterBef.p(0., 0., pzSign*eDau, eDau);
9208 bool History::isQCD2to2(
const Event& event) {
9210 if (!mergingHooksPtr->doWeakClustering())
return false;
9213 int nFinalPartons = 0, nFinal = 0;;
9214 for (
int i = 0;i <
event.size();++i)
9215 if (event[i].isFinal()) {
9217 if ( event[i].idAbs() < 10 ||
event[i].idAbs() == 21)
9220 if (nFinalPartons == 2 && nFinal == 2)
return true;
9230 bool History::isEW2to1(
const Event& event) {
9232 if (!mergingHooksPtr->doWeakClustering())
return false;
9235 for (
int i = 0;i <
event.size();++i) {
9236 if (event[i].isFinal()) {
9237 if (event[i].idAbs() == 23 ||
9238 event[i].idAbs() == 24 ||
9239 event[i].idAbs() == 22) nVector++;
9245 if (nVector == 1)
return true;
9255 void History::setSelectedChild() {
9256 if (mother == 0)
return;
9257 for (
int i = 0;i < int(mother->children.size());++i)
9258 if (mother->children[i] ==
this) mother->selectedChild = i;
9259 mother->setSelectedChild();
9264 void History::setupSimpleWeakShower(
int nSteps) {
9267 if (selectedChild != -1) {
9268 children[selectedChild]->setupSimpleWeakShower(nSteps+1);
9273 vector<int> mode, fermionLines;
9275 vector<pair<int,int> > dipoles;
9278 setupWeakHard(mode,fermionLines, mom);
9281 if (isQCD2to2(state)) {
9283 if (state[3].idAbs() < 10) dipoles.push_back(make_pair(3,4));
9284 if (state[4].idAbs() < 10) dipoles.push_back(make_pair(4,3));
9285 if (state[5].idAbs() < 10) dipoles.push_back(make_pair(5,6));
9286 if (state[6].idAbs() < 10) dipoles.push_back(make_pair(6,5));
9287 }
else if (isEW2to1(state)) {
9288 if (state[3].idAbs() < 10) dipoles.push_back(make_pair(3,4));
9289 if (state[4].idAbs() < 10) dipoles.push_back(make_pair(4,3));
9293 transferSimpleWeakShower(mode, mom, fermionLines, dipoles, nSteps);
9299 void History::transferSimpleWeakShower(vector<int> &mode, vector<Vec4> &mom,
9300 vector<int> fermionLines, vector<pair<int,int> > &dipoles,
9305 infoPtr->setWeakModes(mode);
9306 infoPtr->setWeakDipoles(dipoles);
9307 infoPtr->setWeakMomenta(mom);
9308 infoPtr->setWeak2to2lines(fermionLines);
9313 map<int,int> stateTransfer;
9314 findStateTransfer(stateTransfer);
9317 vector<int> modeNew = updateWeakModes(mode, stateTransfer);
9318 vector<int> fermionLinesNew = updateWeakFermionLines(fermionLines,
9320 vector<pair<int, int> > dipolesNew = updateWeakDipoles(dipoles,
9324 mother->transferSimpleWeakShower(modeNew, mom, fermionLinesNew, dipolesNew,
9331 vector<int> History::updateWeakModes(vector<int>& mode,
9332 map<int,int>& stateTransfer) {
9334 vector<int> modeNew(mode.size() + 1,0);
9337 for (map<int,int>::iterator it = stateTransfer.begin();
9338 it != stateTransfer.end(); ++it)
9339 modeNew[it->second] = mode[it->first];
9341 modeNew[clusterIn.emitted] = mode[clusterIn.radBef];
9345 if (state[clusterIn.radBef].idAbs() == 21 &&
9346 mother->state[clusterIn.emittor].idAbs() != 21) {
9348 if (state[clusterIn.radBef].status() > 0) modeNew[clusterIn.emittor] = 1;
9351 if (modeNew[clusterIn.emittor] == 1);
9352 else if ( mother->state[clusterIn.recoiler].id() == 21)
9353 modeNew[clusterIn.emittor] = 2;
9354 else if ( mother->state[clusterIn.recoiler].id()
9355 == mother->state[clusterIn.emittor].id() )
9356 modeNew[clusterIn.emittor] = 4;
9357 else modeNew[clusterIn.emittor] = 3;
9360 modeNew[clusterIn.emitted] = 1;
9364 if (state[clusterIn.radBef].idAbs() < 10 &&
9365 mother->state[clusterIn.emittor].idAbs() == 21) {
9366 if (state[clusterIn.radBef].status() < 0) {
9367 modeNew[clusterIn.emitted] = 1;
9372 if (state[clusterIn.radBef].idAbs() == 22) {
9374 if (state[clusterIn.radBef].status() > 0) modeNew[clusterIn.emittor] = 1;
9377 if (modeNew[clusterIn.emittor] == 1);
9378 else if ( mother->state[clusterIn.recoiler].id() == 21)
9379 modeNew[clusterIn.emittor] = 2;
9380 else if ( mother->state[clusterIn.recoiler].id()
9381 == mother->state[clusterIn.emittor].id() )
9382 modeNew[clusterIn.emittor] = 4;
9383 else modeNew[clusterIn.emittor] = 3;
9386 modeNew[clusterIn.emitted] = 1;
9395 vector<int> History::updateWeakFermionLines(vector<int> fermionLines,
9396 map<int,int>& stateTransfer) {
9399 if (!fermionLines.empty()) {
9401 fermionLines[0] = stateTransfer[fermionLines[0]];
9402 fermionLines[1] = stateTransfer[fermionLines[1]];
9405 bool lines[2] = {
false,
false};
9406 if (fermionLines[2] != clusterIn.radBef)
9407 fermionLines[2] = stateTransfer[fermionLines[2]];
9408 else lines[0] =
true;
9409 if (fermionLines[3] != clusterIn.radBef)
9410 fermionLines[3] = stateTransfer[fermionLines[3]];
9411 else lines[1] =
true;
9414 for (
int i = 0;i < 2; ++i) {
9416 if (state[fermionLines[2 + i]].isQuark() ||
9417 state[fermionLines[2 + i]].isLepton()) {
9418 if (mother->state[clusterIn.emittor].isQuark() ||
9419 mother->state[clusterIn.emittor].isLepton())
9420 fermionLines[2 + i] = clusterIn.emittor;
9421 else fermionLines[2 + i] = clusterIn.emitted;
9424 else fermionLines[2 + i] = 0;
9428 return fermionLines;
9434 vector<pair<int,int> > History::updateWeakDipoles(
9435 vector<pair<int,int> > &dipoles, map<int,int>& stateTransfer) {
9437 vector<pair<int,int> > dipolesNew;
9438 for (
int i = 0;i < int(dipoles.size());++i) {
9439 int iRecNew = -1, iRadNew = -1;
9442 if (dipoles[i].first != clusterIn.radBef)
9443 iRadNew = stateTransfer[dipoles[i].first];
9445 else if (state[clusterIn.radBef].status() > 0) {
9446 if (mother->state[clusterIn.emitted].id() ==
9447 state[clusterIn.radBef].id())
9448 iRadNew = clusterIn.emitted;
9449 else iRadNew = clusterIn.emittor;
9451 }
else if (mother->state[clusterIn.emittor].idAbs() < 10)
9452 iRadNew = clusterIn.emittor;
9455 if (iRadNew == -1)
continue;
9458 if (dipoles[i].second != clusterIn.radBef)
9459 iRecNew = stateTransfer[dipoles[i].second];
9461 else if (state[clusterIn.radBef].status() > 0) {
9463 if (mother->state[clusterIn.emitted].id() == 21 &&
9464 mother->state[clusterIn.emittor].id() == 21) {
9465 double m1 = (mother->state[clusterIn.emitted].p()
9466 + mother->state[iRadNew].p()).m2Calc();
9467 double m2 = (mother->state[clusterIn.emittor].p()
9468 + mother->state[iRadNew].p()).m2Calc();
9469 iRecNew = (m1 > m2) ? clusterIn.emitted : clusterIn.emittor;
9472 else if (mother->state[clusterIn.emitted].id() ==
9473 state[clusterIn.radBef].id())
9474 iRecNew = clusterIn.emitted;
9475 else iRecNew = clusterIn.emittor;
9477 }
else iRecNew = clusterIn.emittor;
9479 dipolesNew.push_back(make_pair(iRadNew,iRecNew));
9483 if (state[clusterIn.radBef].idAbs() == 21 &&
9484 mother->state[clusterIn.emittor].idAbs() != 21) {
9486 if (state[clusterIn.radBef].status() > 0) {
9487 dipolesNew.push_back(make_pair(clusterIn.emittor,clusterIn.emitted));
9488 dipolesNew.push_back(make_pair(clusterIn.emitted,clusterIn.emittor));
9491 int iRad = clusterIn.emittor;
9492 int iRec = (iRad == 3) ? 4 : 3;
9493 dipolesNew.push_back(make_pair(iRad,iRec));
9494 dipolesNew.push_back(make_pair(clusterIn.emitted,findISRRecoiler()));
9499 if (state[clusterIn.radBef].idAbs() < 10 &&
9500 mother->state[clusterIn.emittor].idAbs() == 21 &&
9501 state[clusterIn.radBef].status() < 0)
9502 dipolesNew.push_back(make_pair(clusterIn.emitted,findISRRecoiler()));
9511 void History::setupWeakHard(vector<int>& mode, vector<int>& fermionLines,
9512 vector<Vec4>& mom) {
9514 if (!isQCD2to2(state)) {
9516 mode.resize(state.size(), 1);
9520 for (
int i = 0;i < 4; ++i) {
9521 mom.push_back(state[3 + i].p());
9522 fermionLines.push_back(3 + i);
9525 if ( state[3].idAbs() == 21 && state[4].idAbs() == 21 &&
9526 state[5].idAbs() == 21 && state[6].idAbs() == 21)
9527 mode.resize(state.size(), 1);
9530 else if (state[5].
id() == -state[6].id() ||
9531 (state[5].idAbs() == 21 && state[6].idAbs() == 21))
9532 mode.resize(state.size(), 1);
9535 else if (state[5].idAbs() == 21 || state[6].idAbs() == 21) {
9536 mode.resize(state.size(), 2);
9537 if (state[3].
id() != state[5].id()) {
9538 swap(mom[0], mom[1]);
9539 swap(mom[2], mom[3]);
9544 else if (state[5].
id() != state[6].
id()) {
9545 mode.resize(state.size(), 3);
9546 if (state[3].
id() != state[5].id()) {
9547 swap(mom[0], mom[1]);
9548 swap(mom[2], mom[3]);
9554 else if (state[5].
id() == state[6].
id()) {
9555 mode.resize(state.size(), 4);
9564 double History::getShowerPluginScale(
const Event& event,
int rad,
int emt,
9565 int rec,
string key,
double scalePythia) {
9568 if ( !mergingHooksPtr->useShowerPlugin() )
return scalePythia;
9571 map<string,double> stateVars;
9572 bool isFSR = showers->timesPtr->isTimelike(event, rad, emt, rec,
"");
9574 string name = showers->timesPtr->getSplittingName(event, rad, emt,
9576 stateVars = showers->timesPtr->getStateVariables(event, rad, emt, rec,
9579 string name = showers->spacePtr->getSplittingName(event, rad, emt,
9581 stateVars = showers->spacePtr->getStateVariables(event, rad, emt, rec,
9585 return ( (stateVars.size() > 0 && stateVars.find(key) != stateVars.end())
9586 ? stateVars[key] : -1.0 );