10 #include "Pythia8/History.h"
26 void Clustering::list()
const {
27 cout <<
" emt " << emitted
29 <<
" rec " << recoiler
30 <<
" partner " << partner
31 <<
" pTscale " << pTscale << endl;
53 const int History::NTRIAL = 1;
70 History::History(
int depth,
74 MergingHooks* mergingHooksPtrIn,
77 ParticleData* particleDataPtrIn,
79 PartonLevel* showersIn,
81 bool isOrdered =
true,
82 bool isStronglyOrdered =
true,
83 bool isAllowed =
true,
84 bool isNextInInput =
true,
93 foundOrderedPath(false),
94 foundStronglyOrderedPath(false),
95 foundAllowedPath(false),
96 foundCompletePath(false),
98 nextInInput(isNextInInput),
103 mergingHooksPtr(mergingHooksPtrIn),
106 particleDataPtr(particleDataPtrIn),
109 coupSMPtr(coupSMPtrIn)
117 if (mother && mergingHooksPtr->includeRedundant()) prob *= pdfForSudakov();
122 double acoll = (mother->state[clusterIn.emittor].isFinal())
123 ? mergingHooksPtr->herwigAcollFSR()
124 : mergingHooksPtr->herwigAcollISR();
125 sumScalarPT = mother->sumScalarPT + acoll*scale;
130 if ( mother ) iReclusteredOld = mother->iReclusteredNew;
133 int nFinalP = 0, nFinalW = 0, nFinalZ = 0;
134 int nL = 0, nA= 0, nH = 0;
135 for (
int i = 0; i < int(state.size()); ++i )
136 if ( state[i].isFinal() ) {
137 if ( state[i].colType() != 0 )
139 if ( state[i].idAbs() == 23 )
141 if ( state[i].idAbs() == 24 )
143 if ( state[i].idAbs() < 20 && state[i].idAbs() > 10)
145 if ( state[i].idAbs() == 22)
147 if ( state[i].idAbs() == 23
148 || state[i].idAbs() == 24
149 || state[i].idAbs() == 25)
152 if ( mergingHooksPtr->doWeakClustering()
153 && nFinalP == 2 && nFinalW == 0 && nFinalZ == 0) depth = 0;
158 bool qcd = ( nFinalP > mergingHooksPtr->hardProcess->nQuarksOut() );
162 vector<Clustering> clusterings;
163 if ( qcd && depth > 0 ) clusterings = getAllQCDClusterings();
165 bool dow = ( mergingHooksPtr->doWeakClustering()
166 && nFinalP > 1 && nFinalW+nFinalZ > 0 );
169 vector<Clustering> clusteringsEW;
171 if ( depth > 0 && dow )
172 clusteringsEW = getAllEWClusterings();
173 if ( !clusteringsEW.empty() ) {
174 clusterings.insert( clusterings.end(), clusteringsEW.begin(),
175 clusteringsEW.end() );
179 vector<Clustering> clusteringsSQCD;
180 if ( depth > 0 && mergingHooksPtr->doSQCDClustering() )
181 clusteringsSQCD = getAllSQCDClusterings();
182 if ( !clusteringsSQCD.empty() )
183 clusterings.insert( clusterings.end(), clusteringsSQCD.begin(),
184 clusteringsSQCD.end() );
188 if ( clusterings.empty() ) {
190 prob *= hardProcessME(state);
191 registerPath( *
this, isOrdered, isStronglyOrdered, isAllowed, depth == 0 );
197 multimap<double, Clustering *> sorted;
198 for (
int i = 0, N = clusterings.size(); i < N; ++i ) {
199 sorted.insert(make_pair(clusterings[i].pT(), &clusterings[i]));
202 for ( multimap<double, Clustering *>::iterator it = sorted.begin();
203 it != sorted.end(); ++it ) {
207 bool stronglyOrdered = isStronglyOrdered;
208 if ( mergingHooksPtr->enforceStrongOrdering()
209 && ( !stronglyOrdered
210 || ( mother && ( it->first <
211 mergingHooksPtr->scaleSeparationFactor()*scale ) ))) {
212 if ( onlyStronglyOrderedPaths() )
continue;
213 stronglyOrdered =
false;
217 bool ordered = isOrdered;
218 if ( mergingHooksPtr->orderInRapidity()
219 && mergingHooksPtr->orderHistories() ) {
221 double z = getCurrentZ((*it->second).emittor,
222 (*it->second).recoiler,(*it->second).emitted,
223 (*it->second).flavRadBef);
225 double zOld = (!mother) ? 0. : mother->getCurrentZ(clusterIn.emittor,
226 clusterIn.recoiler,clusterIn.emitted,
227 clusterIn.flavRadBef);
230 if ( !ordered || ( mother && (it->first < scale
231 || it->first < pow(1. - z,2) / (z * (1. - zOld ))*scale ))) {
232 if ( onlyOrderedPaths() )
continue;
236 }
else if ( mergingHooksPtr->orderHistories() ) {
240 if ( !ordered || ( mother && (it->first < scale) ) ) {
241 if ( onlyOrderedPaths() && onlyAllowedPaths() )
continue;
247 bool doCut = mergingHooksPtr->canCutOnRecState()
248 || mergingHooksPtr->allowCutOnRecState();
249 bool allowed = isAllowed;
251 && mergingHooksPtr->doCutOnRecState(cluster(*it->second)) ) {
252 if ( onlyAllowedPaths() )
continue;
258 children.push_back(
new History(depth - 1,it->first,cluster(*it->second),
259 *it->second, mergingHooksPtr, beamA, beamB, particleDataPtr,
260 infoPtr, showers, coupSMPtr, ordered, stronglyOrdered, allowed,
261 true, prob*getProb(*it->second),
this ));
269 bool History::projectOntoDesiredHistories() {
271 return trimHistories();
286 double History::weightTREE(PartonLevel* trial, AlphaStrong * asFSR,
287 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN) {
289 if ( mergingHooksPtr->canCutOnRecState() && !foundAllowedPath ) {
290 string message=
"Warning in History::weightTREE: No allowed history";
291 message+=
" found. Using disallowed history.";
292 infoPtr->errorMsg(message);
294 if ( mergingHooksPtr->orderHistories() && !foundOrderedPath ) {
295 string message=
"Warning in History::weightTREE: No ordered history";
296 message+=
" found. Using unordered history.";
297 infoPtr->errorMsg(message);
299 if ( mergingHooksPtr->canCutOnRecState()
300 && mergingHooksPtr->orderHistories()
301 && !foundAllowedPath && !foundOrderedPath ) {
302 string message=
"Warning in History::weightTREE: No allowed or ordered";
303 message+=
" history found.";
304 infoPtr->errorMsg(message);
308 double asME = infoPtr->alphaS();
309 double aemME = infoPtr->alphaEM();
310 double maxScale = (foundCompletePath) ? infoPtr->eCM()
311 : mergingHooksPtr->muFinME();
314 History * selected = select(RN);
317 selected->setScalesInHistory();
321 double asWeight = 1.;
322 double aemWeight = 1.;
323 double pdfWeight = 1.;
326 sudakov = selected->weightTree( trial, asME, aemME, maxScale,
327 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
328 aemWeight, pdfWeight );
331 int njetsMaxMPI = mergingHooksPtr->nMinMPI();
332 double mpiwt = selected->weightTreeEmissions( trial, -1, 0, njetsMaxMPI,
336 bool resetScales = mergingHooksPtr->resetHardQRen();
342 && mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0) {
344 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
345 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
346 asWeight *= pow2(runningCoupling);
347 }
else if (mergingHooksPtr->doWeakClustering()
348 && isQCD2to2(selected->state)) {
350 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
351 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
352 asWeight *= pow2(runningCoupling);
356 if (mergingHooksPtr->doWeakClustering() && isEW2to1(selected->state)) {
358 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
359 double runningCoupling = (*aemFSR).alphaEM(newQ2Ren) / aemME;
360 aemWeight *= runningCoupling;
367 && mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) {
369 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
370 double runningCoupling =
371 (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME;
372 asWeight *= runningCoupling;
376 return (sudakov*asWeight*aemWeight*pdfWeight*mpiwt);
385 double History::weightLOOP(PartonLevel* trial,
double RN ) {
387 if ( mergingHooksPtr->canCutOnRecState() && !foundAllowedPath ) {
388 string message=
"Warning in History::weightLOOP: No allowed history";
389 message+=
" found. Using disallowed history.";
390 infoPtr->errorMsg(message);
394 History * selected = select(RN);
396 selected->setScalesInHistory();
402 double maxScale = (foundCompletePath) ? infoPtr->eCM()
403 : mergingHooksPtr->muFinME();
404 int njetsMaxMPI = mergingHooksPtr->nMinMPI();
405 double mpiwt = selected->weightTreeEmissions( trial, -1, 0, njetsMaxMPI,
416 double History::weightFIRST(PartonLevel* trial, AlphaStrong* asFSR,
417 AlphaStrong* asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
421 if (
false) cout << aemFSR << aemISR;
424 double asME = infoPtr->alphaS();
425 double muR = mergingHooksPtr->muRinME();
426 double maxScale = (foundCompletePath)
428 : mergingHooksPtr->muFinME();
431 History * selected = select(RN);
433 selected->setScalesInHistory();
435 double nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
438 double kFactor = asME * mergingHooksPtr->k1Factor(nSteps);
442 double wt = 1. + kFactor;
445 wt += selected->weightFirst(trial,asME, muR, maxScale, asFSR, asISR,
449 double startingScale = (selected->mother) ? state.scale() : infoPtr->eCM();
455 double nWeight1 = 0.;
456 for(
int i=0; i < NTRIAL; ++i) {
458 vector<double> unresolvedEmissionTerm = countEmissions( trial,
459 startingScale, mergingHooksPtr->tms(), 2, asME, asFSR, asISR, 3,
461 nWeight1 += unresolvedEmissionTerm[1];
464 wt += nWeight1/double(NTRIAL);
473 double History::weight_UMEPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,
474 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN) {
476 return weightTREE( trial, asFSR, asISR, aemFSR, aemISR, RN);
483 double History::weight_UMEPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,
484 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN ) {
487 double asME = infoPtr->alphaS();
488 double aemME = infoPtr->alphaEM();
489 double maxScale = (foundCompletePath) ? infoPtr->eCM()
490 : mergingHooksPtr->muFinME();
492 History * selected = select(RN);
494 selected->setScalesInHistory();
498 double asWeight = 1.;
499 double aemWeight = 1.;
500 double pdfWeight = 1.;
503 sudakov = selected->weightTree(trial, asME, aemME, maxScale,
504 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
505 aemWeight, pdfWeight);
508 int njetsMaxMPI = mergingHooksPtr->nMinMPI()+1;
509 double mpiwt = selected->weightTreeEmissions( trial, -1, 0, njetsMaxMPI,
513 bool resetScales = mergingHooksPtr->resetHardQRen();
518 && mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0) {
520 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
521 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
522 asWeight *= pow(runningCoupling,2);
529 && mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) {
531 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
532 double runningCoupling =
533 (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME;
534 asWeight *= runningCoupling;
538 return (sudakov*asWeight*aemWeight*pdfWeight*mpiwt);
544 double History::weight_UNLOPS_TREE(PartonLevel* trial, AlphaStrong * asFSR,
545 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
549 double asME = infoPtr->alphaS();
550 double aemME = infoPtr->alphaEM();
551 double maxScale = (foundCompletePath) ? infoPtr->eCM()
552 : mergingHooksPtr->muFinME();
554 History * selected = select(RN);
556 selected->setScalesInHistory();
559 double asWeight = 1.;
560 double aemWeight = 1.;
561 double pdfWeight = 1.;
565 if (depth < 0) wt = selected->weightTree(trial, asME, aemME, maxScale,
566 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
567 aemWeight, pdfWeight);
569 wt = selected->weightTreeEmissions( trial, 1, 0, depth, maxScale );
570 if (wt != 0.) asWeight = selected->weightTreeALPHAS( asME, asFSR, asISR,
572 if (wt != 0.) aemWeight = selected->weightTreeALPHAEM( aemME, aemFSR,
574 if (wt != 0.) pdfWeight = selected->weightTreePDFs( maxScale,
575 selected->clusterIn.pT(), depth);
579 int njetsMaxMPI = mergingHooksPtr->nMinMPI();
580 double mpiwt = selected->weightTreeEmissions( trial, -1, 0, njetsMaxMPI,
584 bool resetScales = mergingHooksPtr->resetHardQRen();
589 && mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0) {
591 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
592 double runningCoupling = (*asFSR).alphaS(newQ2Ren) / asME;
593 asWeight *= pow(runningCoupling,2);
600 && mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) {
602 double newQ2Ren = pow2( selected->hardRenScale(selected->state) );
603 double runningCoupling =
604 (*asISR).alphaS( newQ2Ren + pow(mergingHooksPtr->pT0ISR(),2) ) / asME;
605 asWeight *= runningCoupling;
609 return (wt*asWeight*aemWeight*pdfWeight*mpiwt);
615 double History::weight_UNLOPS_LOOP(PartonLevel* trial, AlphaStrong * asFSR,
616 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
619 if (depth < 0)
return weightLOOP(trial, RN);
620 else return weight_UNLOPS_TREE(trial, asFSR,asISR, aemFSR,aemISR, RN,depth);
625 double History::weight_UNLOPS_SUBT(PartonLevel* trial, AlphaStrong * asFSR,
626 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
630 History * selected = select(RN);
632 selected->setScalesInHistory();
637 double asME = infoPtr->alphaS();
638 double aemME = infoPtr->alphaEM();
639 double maxScale = (foundCompletePath)
641 : mergingHooksPtr->muFinME();
645 double nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
646 if ( nSteps == 2 && mergingHooksPtr->nRecluster() == 2
647 && ( !foundCompletePath
648 || !selected->allIntermediateAboveRhoMS( mergingHooksPtr->tms() )) )
652 double asWeight = 1.;
653 double aemWeight = 1.;
654 double pdfWeight = 1.;
658 sudakov = selected->weightTree(trial, asME, aemME, maxScale,
659 selected->clusterIn.pT(), asFSR, asISR, aemFSR, aemISR, asWeight,
660 aemWeight, pdfWeight);
662 sudakov = selected->weightTreeEmissions( trial, 1, 0, depth, maxScale );
663 if (sudakov > 0.) asWeight = selected->weightTreeALPHAS( asME, asFSR,
665 if (sudakov > 0.) aemWeight = selected->weightTreeALPHAEM( aemME, aemFSR,
667 if (sudakov > 0.) pdfWeight = selected->weightTreePDFs( maxScale,
668 selected->clusterIn.pT(), depth);
672 int njetsMaxMPI = mergingHooksPtr->nMinMPI()+1;
673 double mpiwt = selected->weightTreeEmissions( trial, -1, 0, njetsMaxMPI,
677 wt = ( mergingHooksPtr->nRecluster() == 2 ) ? 1.
678 : asWeight*aemWeight*pdfWeight*sudakov*mpiwt;
687 double History::weight_UNLOPS_SUBTNLO(PartonLevel* trial, AlphaStrong * asFSR,
688 AlphaStrong * asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
double RN,
694 History * selected = select(RN);
696 selected->setScalesInHistory();
700 double maxScale = (foundCompletePath) ? infoPtr->eCM()
701 : mergingHooksPtr->muFinME();
702 int njetsMaxMPI = mergingHooksPtr->nMinMPI()+1;
703 double mpiwt = selected->weightTreeEmissions( trial, -1, 0, njetsMaxMPI,
709 }
else return weight_UNLOPS_SUBT(trial, asFSR, asISR, aemFSR, aemISR, RN,
718 double History::weight_UNLOPS_CORRECTION(
int order, PartonLevel* trial,
719 AlphaStrong* asFSR, AlphaStrong* asISR, AlphaEM * aemFSR, AlphaEM * aemISR,
720 double RN, Rndm* rndmPtr ) {
723 if (
false) cout << aemFSR << aemISR;
726 if ( order < 0 )
return 0.;
729 double asME = infoPtr->alphaS();
731 double muR = mergingHooksPtr->muRinME();
732 double maxScale = (foundCompletePath)
734 : mergingHooksPtr->muFinME();
737 History * selected = select(RN);
739 selected->setScalesInHistory();
741 double nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
744 double kFactor = asME * mergingHooksPtr->k1Factor(nSteps);
751 if ( order == 0 )
return wt;
761 double wA = selected->weightFirstALPHAS( asME, muR, asFSR, asISR );
763 wt += (fixas) ? wA : 0.;
766 for (
int i=0; i < NTRIAL; ++i ) {
768 double wE = selected->weightFirstEmissions(trial,asME, maxScale,
769 asFSR, asISR, fixpdf, fixas );
773 double pscale = selected->clusterIn.pT();
774 double wP = selected->weightFirstPDFs(asME, maxScale, pscale, rndmPtr);
776 nWeight += (fixpdf) ? wP : 0.;
778 wt += nWeight/double(NTRIAL);
781 if ( order == 1 )
return wt;
792 void History::getStartingConditions(
const double RN,
Event& outState ) {
795 History * selected = select(RN);
798 selected->setScalesInHistory();
801 int nSteps = mergingHooksPtr->getNumberOfClusteringSteps(state);
804 if (!selected->mother) {
806 for(
int i=0; i < int(state.size()); ++i)
807 if ( state[i].isFinal()) nFinal++;
809 state.scale(mergingHooksPtr->muF());
815 if (mergingHooksPtr->getNumberOfClusteringSteps(state) == 0) {
816 infoPtr->zNowISR(0.5);
817 infoPtr->pT2NowISR(pow(state[0].e(),2));
818 infoPtr->hasHistory(
true);
821 infoPtr->zNowISR(selected->zISR());
822 infoPtr->pT2NowISR(pow(selected->pTISR(),2));
823 infoPtr->hasHistory(
true);
829 double muf = state[0].e();
830 for (
int i=0; i < state.size(); ++i )
831 if ( state[i].isFinal()
832 && ( state[i].colType() != 0 || state[i].id() == 22 ) ) {
834 muf = min( muf, abs(state[i].mT()) );
838 if ( nSteps == 0 && nFinalCol == 2
839 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
840 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0) ) {
842 for (
int i = 3;i < state.size();++i)
847 if (nSteps == 0 && nFinalCol == 2 &&
848 mergingHooksPtr->getProcessString().find(
"inc") != string::npos) {
850 for (
int i = 3;i < state.size();++i)
852 for (
int i=0; i < min(state.size(),outState.size()); ++i )
853 state[i].pol(outState[i].pol());
861 infoPtr->zNowISR(selected->zISR());
862 infoPtr->pT2NowISR(pow(selected->pTISR(),2));
863 infoPtr->hasHistory(
true);
872 mergingHooksPtr->muMI(infoPtr->eCM());
874 mergingHooksPtr->muMI(outState.scale());
877 if (mergingHooksPtr->doWeakClustering()) setupWeakShower(0);
885 void History::printHistory(
const double RN ) {
886 History * selected = select(RN);
887 selected->printStates();
895 void History::printStates() {
897 cout << scientific << setprecision(6) <<
"Probability=" << prob << endl;
903 double p = (mother) ? prob/mother->prob : prob;
904 cout << scientific << setprecision(6) <<
"Probability=" << p
905 <<
" scale=" << clusterIn.pT() << endl;
908 mother->printStates();
917 bool History::getClusteredEvent(
const double RN,
int nSteps,
921 History * selected = select(RN);
925 selected->setScalesInHistory();
928 if (nSteps > selected->nClusterings())
return false;
931 outState = selected->clusteredState(nSteps-1);
939 bool History::getFirstClusteredEventAboveTMS(
const double RN,
int nDesired,
940 Event& process,
int& nPerformed,
bool doUpdate ) {
943 int nTried = nDesired - 1;
945 int nSteps = select(RN)->nClusterings();
947 select(RN)->setScalesInHistory();
954 dummy.init(
"(hard process-modified)", particleDataPtr );
959 if ( !getClusteredEvent( RN, nSteps-nTried+1, dummy ) )
return false;
960 if ( nTried >= nSteps )
break;
963 }
while ( mergingHooksPtr->getNumberOfClusteringSteps(dummy) > 0
964 && mergingHooksPtr->tmsNow( dummy) < mergingHooksPtr->tms() );
967 if ( doUpdate ) process = dummy;
970 if ( nTried > nSteps )
return false;
975 mergingHooksPtr->nReclusterSave = nPerformed;
977 if (mergingHooksPtr->getNumberOfClusteringSteps(state) == 0)
978 mergingHooksPtr->muMI(infoPtr->eCM());
980 mergingHooksPtr->muMI(state.scale());
992 double History::getPDFratio(
int side,
bool forSudakov,
bool useHardPDFs,
993 int flavNum,
double xNum,
double muNum,
994 int flavDen,
double xDen,
double muDen) {
997 if ( abs(flavNum) > 10 && flavNum != 21 )
return 1.0;
998 if ( abs(flavDen) > 10 && flavDen != 21 )
return 1.0;
1001 double pdfRatio = 1.0;
1004 double pdfNum = 0.0;
1005 double pdfDen = 0.0;
1008 if ( useHardPDFs ) {
1011 pdfNum = mother->beamA.xfHard( flavNum, xNum, muNum*muNum);
1013 pdfNum = beamA.xfHard( flavNum, xNum, muNum*muNum);
1015 pdfDen = max(1e-10, beamA.xfHard( flavDen, xDen, muDen*muDen));
1017 pdfDen = max(1e-10, beamA.xfHard( flavDen, xDen, muDen*muDen));
1020 pdfNum = mother->beamB.xfHard( flavNum, xNum, muNum*muNum);
1022 pdfNum = beamB.xfHard( flavNum, xNum, muNum*muNum);
1025 pdfDen = max(1e-10,beamB.xfHard( flavDen, xDen, muDen*muDen));
1027 pdfDen = max(1e-10,beamB.xfHard( flavDen, xDen, muDen*muDen));
1034 pdfNum = mother->beamA.xfISR(0, flavNum, xNum, muNum*muNum);
1036 pdfNum = beamA.xfISR(0, flavNum, xNum, muNum*muNum);
1038 pdfDen = max(1e-10, beamA.xfISR(0, flavDen, xDen, muDen*muDen));
1040 pdfDen = max(1e-10, beamA.xfISR(0, flavDen, xDen, muDen*muDen));
1044 pdfNum = mother->beamB.xfISR(0, flavNum, xNum, muNum*muNum);
1046 pdfNum = beamB.xfISR(0, flavNum, xNum, muNum*muNum);
1049 pdfDen = max(1e-10,beamB.xfISR(0, flavDen, xDen, muDen*muDen));
1051 pdfDen = max(1e-10,beamB.xfISR(0, flavDen, xDen, muDen*muDen));
1056 if ( forSudakov && abs(flavNum) ==4 && abs(flavDen) == 4 && muDen == muNum
1057 && muNum < particleDataPtr->m0(4))
1058 pdfDen = pdfNum = 1.0;
1061 if ( pdfNum > 1e-15 && pdfDen > 1e-10 ) {
1062 pdfRatio *= pdfNum / pdfDen;
1063 }
else if ( pdfNum < pdfDen ) {
1065 }
else if ( pdfNum > pdfDen ) {
1081 void History::setScalesInHistory() {
1089 setScales(ident,
true);
1100 double History::getWeakProb() {
1101 vector<int> modes, fermionLines;
1103 return getWeakProb(modes, mom, fermionLines);
1112 double History::getWeakProb(vector<int> &mode, vector<Vec4> &mom,
1113 vector<int> fermionLines) {
1116 if (!mother)
return 1.;
1119 map<int,int> stateTransfer;
1120 findStateTransfer(stateTransfer);
1123 if (mode.empty()) setupWeakHard(mode,fermionLines,mom);
1126 vector<int> modeNew = updateWeakModes(mode, stateTransfer);
1127 vector<int> fermionLinesNew = updateWeakFermionLines(fermionLines,
1131 if (mother->state[clusterIn.emitted].idAbs() == 24 ||
1132 mother->state[clusterIn.emitted].idAbs() == 23)
1133 return getSingleWeakProb(modeNew, mom, fermionLinesNew) *
1134 mother->getWeakProb(modeNew, mom, fermionLinesNew);
1135 else return mother->getWeakProb(modeNew, mom, fermionLinesNew);
1140 double History::getSingleWeakProb(vector<int> &mode, vector<Vec4> &mom,
1141 vector<int> fermionLines) {
1144 double weakCoupling = 0.0;
1145 if (mother->state[clusterIn.emitted].idAbs() == 24) {
1147 if (clusterIn.spinRadBef == 1)
return 0.0;
1148 else if (clusterIn.spinRadBef == -1)
1149 weakCoupling = 4.*M_PI/ coupSMPtr->sin2thetaW()
1150 * coupSMPtr->V2CKMid(abs(clusterIn.flavRadBef),
1151 mother->state[clusterIn.emittor].idAbs());
1153 infoPtr->errorMsg(
"Warning in History::getSingleWeakProb: "
1154 "Spin not properly configurated. Skipping history");
1157 }
else if (mother->state[clusterIn.emitted].idAbs() == 23) {
1159 if (clusterIn.spinRadBef == 1)
1160 weakCoupling = 4.*M_PI*pow2(coupSMPtr->rf( abs(clusterIn.flavRadBef)))
1161 / (coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW()) ;
1162 else if (clusterIn.spinRadBef == -1)
1163 weakCoupling = 4.*M_PI*pow2(coupSMPtr->lf( abs(clusterIn.flavRadBef)))
1164 / (coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW()) ;
1166 infoPtr->errorMsg(
"Warning in History::getSingleWeakProb: "
1167 "Spin not properly configurated. Skipping history");
1171 infoPtr->errorMsg(
"Warning in History::getSingleWeakProb: "
1172 "Did not emit W/Z. Skipping history.");
1179 Vec4 pRadAft = mother->state[clusterIn.emittor].p();
1180 Vec4 pEmtAft = mother->state[clusterIn.emitted].p();
1181 Vec4 pRecAft = mother->state[clusterIn.recoiler].p();
1182 Vec4 pSum = pRadAft + pEmtAft + pRecAft;
1183 double m2sum = pSum.m2Calc();
1184 double Qsq = (pRadAft + pEmtAft).m2Calc();
1186 double m2Rad0 = pRadAft.m2Calc();
1187 double m2Emt0 = pEmtAft.m2Calc();
1188 double lambda13 = sqrt( pow2(Qsq - m2Rad0 - m2Emt0 ) - 4. * m2Rad0*m2Emt0 );
1189 double k1 = ( Qsq - lambda13 + (m2Emt0 - m2Rad0 ) ) / ( 2. * Qsq );
1190 double k3 = ( Qsq - lambda13 - (m2Emt0 - m2Rad0 ) ) / ( 2. * Qsq );
1192 double z = mother->getCurrentZ(clusterIn.emittor, clusterIn.recoiler,
1193 clusterIn.emitted, clusterIn.flavRadBef);
1194 double pT2 = pow2(clusterIn.pTscale);
1196 double x1 = 2. * pRadAft * pSum / m2sum;
1197 double x2 = 2. * pRecAft * pSum / m2sum;
1198 double x3 = 2. * pEmtAft * pSum / m2sum;
1201 if ( state[clusterIn.radBef].status() > 0) {
1203 if (mode[clusterIn.emittor] == 1) {
1205 double eCMME = pSum.mCalc();
1206 double r1 = mother->state[clusterIn.emittor].m() / eCMME;
1207 double r2 = mother->state[clusterIn.recoiler].m() / eCMME;
1208 double r3 = mother->state[clusterIn.emitted].m() / eCMME;
1209 double x1s = x1 * x1;
1210 double x2s = x2 * x2;
1211 double r1s = r1 * r1;
1212 double r2s = r2 * r2;
1213 double r3s = r3 * r3;
1214 double prop1 = 1. + r1s - r2s - x1;
1215 double prop2 = 1. + r2s - r1s - x2;
1216 double prop1s = prop1 * prop1;
1217 double prop2s = prop2 * prop2;
1218 double prop12 = prop1 * prop2;
1221 double jac = 1./(1.-z) * 1./pT2 * (1-x2+r2-r1)*(x3 - k1*(x1+x3))
1222 * (1.-x1+r1-r2) / x3;
1223 return jac * weakCoupling * ((2. * r3s * r3s + 2. * r3s *
1224 (x1 + x2) + x1s + x2s) / prop12 - r3s / prop1s - r3s / prop2s);
1229 Vec4 p1 = mother->state[clusterIn.emittor].p();
1230 Vec4 p2 = mother->state[clusterIn.recoiler].p();
1231 Vec4 p3 = mother->state[clusterIn.emitted].p();
1232 Vec4 radBef = state[clusterIn.radBef].p();
1233 Vec4 recBef = state[clusterIn.recBef].p();
1238 if (fermionLines[2] == clusterIn.emittor);
1239 else if (fermionLines[3] == clusterIn.emittor) swap(pIn1, pIn2);
1242 double scaleFactor2 = (p1 + p2 + p3).m2Calc() / (pIn1 + pIn2).m2Calc();
1243 double scaleFactor = sqrt(scaleFactor2);
1244 pIn1 *= scaleFactor;
1245 pIn2 *= scaleFactor;
1249 RotBstMatrix rot2to2frame;
1250 rot2to2frame.bstback(pIn1 + pIn2);
1251 pIn1.rotbst(rot2to2frame);
1252 pIn2.rotbst(rot2to2frame);
1253 p1.rotbst(rot2to2frame);
1254 p2.rotbst(rot2to2frame);
1255 p3.rotbst(rot2to2frame);
1256 recBef.rotbst(rot2to2frame);
1257 radBef.rotbst(rot2to2frame);
1260 RotBstMatrix rot2to3frame;
1261 rot2to3frame.bstback(p1 + p2 + p3);
1262 p1.rotbst(rot2to3frame);
1263 p2.rotbst(rot2to3frame);
1264 p3.rotbst(rot2to3frame);
1265 recBef.rotbst(rot2to3frame);
1266 radBef.rotbst(rot2to3frame);
1269 double sHat = (pIn1 + pIn2).m2Calc();
1270 double tHat = (radBef - pIn1).m2Calc();
1271 double uHat = (recBef - pIn1).m2Calc();
1272 double localProb = 0;
1273 double Q2 = pT2 / (z*(1.-z));
1274 double jac = 1./(4. * M_PI) * 1./( (1.-z) * z ) * sHat / (sHat - Q2)
1278 if (mode[clusterIn.emittor] == 2)
1279 localProb = weakShowerMEs.getMEqg2qgZ( pIn1, pIn2, p2, p3, p1)
1280 / weakShowerMEs.getMEqg2qg( sHat, tHat, uHat);
1281 else if (mode[clusterIn.emittor] == 3)
1282 localProb = weakShowerMEs.getMEqq2qqZ( pIn1, pIn2, p3, p2, p1)
1283 / weakShowerMEs.getMEqq2qq( sHat, tHat, uHat,
false);
1284 else if (mode[clusterIn.emittor] == 4)
1285 localProb = weakShowerMEs.getMEqq2qqZ( pIn1, pIn2, p3, p2, p1)
1286 / weakShowerMEs.getMEqq2qq( sHat, tHat, uHat,
true);
1288 string message=
"Warning in History::getSingleWeakProb: Wrong";
1289 message+=
" mode setup. Setting probability for path to zero.";
1290 infoPtr->errorMsg(message);
1294 localProb *= abs((-p3 + pIn1).m2Calc())
1295 / ((p3 + p1).m2Calc() + abs((-pIn1 + p3).m2Calc()));
1297 return jac * weakCoupling * localProb;
1303 if (mode[clusterIn.emittor] == 1) {
1304 Vec4 pIn1 = mother->state[clusterIn.emittor].p();
1305 Vec4 pIn2 = mother->state[clusterIn.recoiler].p();
1306 Vec4 p1 = mother->state[clusterIn.emitted].p();
1307 Vec4 p2 = pIn1 + pIn2 -p1;
1309 double sH = (pIn1 + pIn2).m2Calc();
1310 double tH = (p1 - pIn1).m2Calc();
1311 double uH = (p1 - pIn2).m2Calc();
1312 double m3s = p1.m2Calc();
1313 double m4s = p2.m2Calc();
1315 double jac = 1./sH * tH*uH / ( tH * (tH + uH) );
1316 return jac * weakCoupling * ((uH*uH + tH*tH + 2 * sH * (m3s + m4s))
1317 / (uH*tH) - m3s * m4s * (1/(tH*tH) + 1/(uH*uH)));
1322 Vec4 pIn1 = mother->state[clusterIn.emittor].p();
1323 Vec4 pIn2 = mother->state[clusterIn.recoiler].p();
1324 Vec4 p3 = mother->state[clusterIn.emitted].p();
1329 if (fermionLines[0] == clusterIn.emittor);
1330 else if (fermionLines[1] == clusterIn.emittor)
1334 Vec4 pDaughter, pRecoiler;
1335 int signBeam = (pIn1.pz() > 0.) ? 1 : -1;
1336 double eCM = state[0].e(), phi = 0;
1339 reverseBoostISR(pIn1, p3, pIn2, pDaughter, pRecoiler, signBeam,
1344 double scaleFactor2 = (pIn1 + pIn2 - p3).m2Calc() / (p1 + p2).m2Calc();
1345 double scaleFactor = sqrt(scaleFactor2);
1346 RotBstMatrix rot2to2frame;
1347 rot2to2frame.bstback(p1 + p2);
1348 p1.rotbst(rot2to2frame);
1349 p2.rotbst(rot2to2frame);
1355 Vec4 radBef = state[clusterIn.radBef].p();
1356 Vec4 recBef = state[clusterIn.recBef].p();
1358 RotBstMatrix rot2to2frameInc;
1359 rot2to2frameInc.bstback(radBef + recBef);
1360 radBef.rotbst(rot2to2frameInc);
1361 recBef.rotbst(rot2to2frameInc);
1362 double sHat = (p1 + p2).m2Calc();
1363 double tHat = (p1 - radBef).m2Calc();
1364 double uHat = (p1 - recBef).m2Calc();
1365 double localProb = 0;
1370 double jac = z / (4.*M_PI);
1372 if (mode[clusterIn.emittor] == 2)
1373 localProb = weakShowerMEs.getMEqg2qgZ(pIn1, pIn2, p2, p3, p1)
1374 / weakShowerMEs.getMEqg2qg(sHat, tHat, uHat);
1375 else if (mode[clusterIn.emittor] == 4)
1376 localProb = weakShowerMEs.getMEqq2qqZ(pIn1, pIn2, p3, p2, p1)
1377 / weakShowerMEs.getMEqq2qq(sHat, tHat, uHat,
true);
1378 else if (mode[clusterIn.emittor] == 3)
1379 localProb = weakShowerMEs.getMEqq2qqZ(pIn1, pIn2, p3, p2, p1)
1380 / weakShowerMEs.getMEqq2qq(sHat, tHat, uHat,
false);
1382 string message=
"Warning in History::getSingleWeakProb: Wrong";
1383 message+=
" mode setup. Setting probability for path to zero.";
1384 infoPtr->errorMsg(message);
1388 localProb *= (p3 + p1).m2Calc() / ( (p3 + p1).m2Calc()
1389 + abs((-pIn1 + p3).m2Calc()) );
1391 return jac * weakCoupling * localProb;
1400 bool History::checkWeakRecoils(map<int,int> &allowedRecoils,
bool isFirst) {
1401 if (!mother)
return true;
1406 if (state.size() != 8) {
1407 if (state[3].isQuark() || state[3].isLepton())
1408 allowedRecoils.insert(pair<int,int>(3,4));
1409 if (state[4].isQuark() || state[4].isLepton())
1410 allowedRecoils.insert(pair<int,int>(4,3));
1413 if (state[3].isQuark() || state[3].isLepton())
1414 allowedRecoils.insert(pair<int,int>(3,4));
1415 if (state[4].isQuark() || state[4].isLepton())
1416 allowedRecoils.insert(pair<int,int>(4,3));
1417 if (state[5].isQuark() || state[5].isLepton())
1418 allowedRecoils.insert(pair<int,int>(5,6));
1419 if (state[6].isQuark() || state[6].isLepton())
1420 allowedRecoils.insert(pair<int,int>(6,5));
1425 map<int,int> transfer;
1426 findStateTransfer(transfer);
1429 map<int,int> allowedRecoilsNew;
1430 for (map<int,int>::iterator it = allowedRecoils.begin();
1431 it != allowedRecoils.end(); ++it) {
1434 if (state[clusterIn.radBef].status() > 0) {
1436 if (it->first != clusterIn.radBef &&
1437 it->second != clusterIn.radBef)
1438 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1439 transfer[it->second]));
1441 else if (it->second == clusterIn.radBef) {
1443 if (state[clusterIn.recBef].isQuark() ||
1444 state[clusterIn.recBef].isLepton()) {
1445 if (mother->state[clusterIn.emittor].isQuark() ||
1446 mother->state[clusterIn.emittor].isLepton())
1447 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1448 clusterIn.emittor));
1450 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1451 clusterIn.emitted));
1455 double mEmittor = (mother->state[clusterIn.emittor].p() +
1456 mother->state[transfer[it->first]].p()).mCalc();
1457 double mEmitted = (mother->state[clusterIn.emitted].p() +
1458 mother->state[transfer[it->first]].p()).mCalc();
1459 if (mEmitted > mEmittor)
1460 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1461 clusterIn.emitted));
1463 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1464 clusterIn.emittor));
1468 if (mother->state[clusterIn.emittor].isQuark() ||
1469 mother->state[clusterIn.emittor].isLepton())
1470 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1471 transfer[it->second]));
1473 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emitted,
1474 transfer[it->second]));
1480 if (it->first != clusterIn.radBef &&
1481 it->second != clusterIn.radBef)
1482 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1483 transfer[it->second]));
1486 else if (it->second == clusterIn.radBef)
1487 allowedRecoilsNew.insert(pair<int,int>(transfer[it->first],
1488 clusterIn.emittor));
1493 if (mother->state[clusterIn.emittor].isQuark() ||
1494 mother->state[clusterIn.emittor].isLepton())
1495 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1496 clusterIn.recoiler));
1500 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1501 findISRRecoiler()));
1508 if ( (state[clusterIn.radBef].idAbs() == 22 ||
1509 state[clusterIn.radBef].idAbs() == 21) &&
1510 (mother->state[clusterIn.emittor].isQuark() ||
1511 mother->state[clusterIn.emittor].isLepton() ) ) {
1513 if (state[clusterIn.radBef].status() > 0) {
1514 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1515 clusterIn.emitted));
1516 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emitted,
1517 clusterIn.emittor));
1522 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emittor,
1523 clusterIn.recoiler));
1524 allowedRecoilsNew.insert(pair<int,int>(clusterIn.emitted,
1525 findISRRecoiler()));
1533 if (mother->state[clusterIn.emitted].idAbs() == 24 ||
1534 mother->state[clusterIn.emitted].idAbs() == 23)
1535 if ( clusterIn.recoiler != allowedRecoilsNew[clusterIn.emittor])
1539 return mother->checkWeakRecoils(allowedRecoilsNew);
1548 int History::findISRRecoiler() {
1550 int flavRad = mother->state[clusterIn.emitted].id();
1551 Vec4 pRad = mother->state[clusterIn.emitted].p();
1552 double mRad = mother->state[clusterIn.emitted].m();
1553 int iRad = clusterIn.emitted;
1555 double ppMin = 1E20;
1556 for (
int i = 0;i < mother->state.size(); ++i) {
1557 if (i == iRad)
continue;
1558 if (mother->state[i].isFinal() && mother->state[i].id() == - flavRad) {
1559 double ppNow = mother->state[i].p() * pRad
1560 - mother->state[i].m() - mRad;
1561 if (ppNow < ppMin) {
1567 if (iRec)
return iRec;
1570 for (
int i = 0;i < mother->state.size(); ++i) {
1571 if (i == iRad)
continue;
1572 if (mother->state[i].isFinal() && mother->state[i].idAbs() < 20) {
1573 double weakCoupNow = 1.;
1574 double ppNow = (mother->state[i].p() * pRad
1575 - mother->state[i].m() - mRad) / weakCoupNow;
1576 if (ppNow < ppMin) {
1582 if (iRec)
return iRec;
1585 for (
int i = 0;i < mother->state.size(); ++i) {
1586 if (i == iRad)
continue;
1587 if (mother->state[i].isFinal()) {
1588 double ppNow = mother->state[i].p() * pRad
1589 - mother->state[i].m() - mRad;
1590 if (ppNow < ppMin) {
1596 if (iRec)
return iRec;
1606 void History::findStateTransfer(map<int,int> &transfer) {
1608 if (!mother)
return;
1612 for(
int i = 0;i < 3; ++i)
1613 transfer.insert(pair<int,int>(i,i));
1615 transfer.insert(pair<int,int>(clusterIn.radBef, clusterIn.emittor));
1616 transfer.insert(pair<int,int>(clusterIn.recBef, clusterIn.recoiler));
1619 for (
int i = 0; i < int(mother->state.size()); ++i) {
1620 if (clusterIn.emitted == i ||
1621 clusterIn.emittor == i ||
1622 clusterIn.recoiler == i)
1625 for (
int j = 0;j < int(state.size()); ++j) {
1626 if (mother->state[i].id() == state[j].id()
1627 && mother->state[i].colType() == state[j].colType()
1628 && mother->state[i].chargeType() == state[j].chargeType()
1629 && mother->state[i].col() == state[j].col()
1630 && mother->state[i].acol() == state[j].acol()
1631 && mother->state[i].status() == state[j].status()) {
1632 transfer.insert(pair<int,int>(j,i));
1649 void History::findPath(vector<int>& out) {
1652 if (!mother &&
int(children.size()) < 1)
return;
1657 int size = int(mother->children.size());
1659 for (
int i=0; i < size; ++i) {
1660 if ( mother->children[i]->scale == scale
1661 && mother->children[i]->prob == prob
1662 && equalClustering(mother->children[i]->clusterIn,clusterIn)) {
1669 out.push_back(iChild);
1670 mother->findPath(out);
1695 void History::setScales( vector<int> index,
bool forward) {
1699 if ( children.empty() && forward ) {
1702 double scaleNew = 1.;
1703 if (mergingHooksPtr->incompleteScalePrescip()==0) {
1704 scaleNew = mergingHooksPtr->muF();
1705 }
else if (mergingHooksPtr->incompleteScalePrescip()==1) {
1707 pOut.p(0.,0.,0.,0.);
1708 for(
int i=0; i<int(state.size()); ++i)
1709 if (state[i].isFinal())
1710 pOut += state[i].p();
1711 scaleNew = pOut.mCalc();
1712 }
else if (mergingHooksPtr->incompleteScalePrescip()==2) {
1713 scaleNew = state[0].e();
1716 scaleNew = max( mergingHooksPtr->pTcut(), scaleNew);
1718 state.scale(scaleNew);
1719 for(
int i=3; i < int(state.size());++i)
1720 if (state[i].colType() != 0)
1721 state[i].scale(scaleNew);
1724 state.scale( state[0].e() );
1726 bool isLEP = ( state[3].isLepton() && state[4].isLepton() );
1728 int nFinalPartons = 0;
1729 int nFinalPhotons = 0;
1730 for (
int i=0; i < int(state.size()); ++i ) {
1731 if ( state[i].isFinal() ) {
1733 if ( state[i].colType() != 0 ) nFinalPartons++;
1734 if ( state[i].
id() == 22 ) nFinalPhotons++;
1737 bool isQCD = ( nFinal == 2 && nFinal == nFinalPartons );
1738 bool isPPh = ( nFinal == 2 && nFinalPartons == 1 && nFinalPhotons == 1);
1740 if ( !isLEP && ( isQCD || isPPh ) ) {
1741 double scaleNew = hardFacScale(state);
1742 state.scale( scaleNew );
1748 if (mother && forward) {
1750 double scaleNew = 1.;
1751 if (mergingHooksPtr->unorderedScalePrescip() == 0) {
1753 scaleNew = max( mergingHooksPtr->pTcut(), max(scale,mother->scale));
1754 }
else if (mergingHooksPtr->unorderedScalePrescip() == 1) {
1756 if (scale < mother->scale)
1757 scaleNew *= max( mergingHooksPtr->pTcut(), min(scale,mother->scale));
1759 scaleNew *= max( mergingHooksPtr->pTcut(), max(scale,mother->scale));
1764 mother->state[clusterIn.emitted].scale(scaleNew);
1765 mother->state[clusterIn.emittor].scale(scaleNew);
1766 mother->state[clusterIn.recoiler].scale(scaleNew);
1770 mother->scaleCopies(clusterIn.emitted, mother->state, scaleNew);
1771 mother->scaleCopies(clusterIn.emittor, mother->state, scaleNew);
1772 mother->scaleCopies(clusterIn.recoiler, mother->state, scaleNew);
1775 mother->setScales(index,
true);
1780 if (!mother || !forward) {
1783 if (
int(index.size()) > 0 ) {
1784 iChild = index.back();
1790 scale = max(mergingHooksPtr->pTcut(), scale);
1793 if (iChild != -1 && !children.empty()) {
1794 if (scale > children[iChild]->scale ) {
1795 if (mergingHooksPtr->unorderedScalePrescip() == 0) {
1797 double scaleNew = max( mergingHooksPtr->pTcut(),
1798 max(scale,children[iChild]->scale));
1800 for(
int i = 0; i < int(children[iChild]->state.size()); ++i)
1801 if (children[iChild]->state[i].scale() == children[iChild]->scale)
1802 children[iChild]->state[i].scale(scaleNew);
1804 children[iChild]->scale = scaleNew;
1806 }
else if ( mergingHooksPtr->unorderedScalePrescip() == 1) {
1808 double scaleNew = max(mergingHooksPtr->pTcut(),
1809 min(scale,children[iChild]->scale));
1811 for(
int i = 0; i < int(state.size()); ++i)
1812 if (state[i].scale() == scale)
1813 state[i].scale(scaleNew);
1819 double scalemin = state[0].e();
1820 for(
int i = 0; i < int(state.size()); ++i)
1821 if (state[i].colType() != 0)
1822 scalemin = max(mergingHooksPtr->pTcut(),
1823 min(scalemin,state[i].scale()));
1824 state.scale(scalemin);
1825 scale = max(mergingHooksPtr->pTcut(), scale);
1828 children[iChild]->setScales(index,
false);
1844 void History::scaleCopies(
int iPart,
const Event& refEvent,
double rho) {
1849 for(
int i=0; i < mother->state.size(); ++i) {
1850 if ( ( mother->state[i].id() == refEvent[iPart].id()
1851 && mother->state[i].colType() == refEvent[iPart].colType()
1852 && mother->state[i].chargeType() == refEvent[iPart].chargeType()
1853 && mother->state[i].col() == refEvent[iPart].col()
1854 && mother->state[i].acol() == refEvent[iPart].acol() )
1857 mother->state[i].scale(rho);
1860 mother->scaleCopies( iPart, refEvent, rho );
1873 void History::setEventScales() {
1877 mother->state.scale(scale);
1879 mother->setEventScales();
1889 double History::zISR() {
1892 if (!mother)
return 0.0;
1894 if (mother->state[clusterIn.emittor].isFinal())
return mother->zISR();
1896 int rad = clusterIn.emittor;
1897 int rec = clusterIn.recoiler;
1898 int emt = clusterIn.emitted;
1899 double z = (mother->state[rad].p() + mother->state[rec].p()
1900 - mother->state[emt].p()).m2Calc()
1901 / (mother->state[rad].p() + mother->state[rec].p()).m2Calc();
1903 double znew = mother->zISR();
1905 if (znew > 0.) z = znew;
1916 double History::zFSR() {
1919 if (!mother)
return 0.0;
1921 if (!mother->state[clusterIn.emittor].isFinal())
return mother->zFSR();
1923 int rad = clusterIn.emittor;
1924 int rec = clusterIn.recoiler;
1925 int emt = clusterIn.emitted;
1927 Vec4 sum = mother->state[rad].p() + mother->state[rec].p()
1928 + mother->state[emt].p();
1929 double m2Dip = sum.m2Calc();
1930 double x1 = 2. * (sum * mother->state[rad].p()) / m2Dip;
1931 double x3 = 2. * (sum * mother->state[emt].p()) / m2Dip;
1933 double z = x1/(x1+x3);
1935 double znew = mother->zFSR();
1937 if (znew > 0.) z = znew;
1948 double History::pTISR() {
1950 if (!mother)
return 0.0;
1952 if (mother->state[clusterIn.emittor].isFinal())
return mother->pTISR();
1953 double pT = mother->state.scale();
1955 double pTnew = mother->pTISR();
1957 if (pTnew > 0.) pT = pTnew;
1968 double History::pTFSR() {
1971 if (!mother)
return 0.0;
1973 if (!mother->state[clusterIn.emittor].isFinal())
return mother->pTFSR();
1974 double pT = mother->state.scale();
1976 double pTnew = mother->pTFSR();
1978 if (pTnew > 0.) pT = pTnew;
1989 int History::nClusterings() {
1990 if (!mother)
return 0;
1991 int w = mother->nClusterings();
2004 Event History::clusteredState(
int nSteps) {
2007 Event outState = state;
2009 if (mother && nSteps > 0)
2010 outState = mother->clusteredState(nSteps - 1);
2023 History * History::select(
double rnd) {
2026 if ( goodBranches.empty() && badBranches.empty() )
return this;
2030 map<double, History*> selectFrom;
2031 if ( !goodBranches.empty() ) {
2032 selectFrom = goodBranches;
2033 sum = sumGoodBranches;
2035 selectFrom = badBranches;
2036 sum = sumBadBranches;
2039 if (mergingHooksPtr->pickBySumPT()) {
2042 for (
int i=0; i < state.size(); ++i)
2043 if (state[i].isFinal())
2046 double sumMin = (nFinal-2)*state[0].e();
2047 for ( map<double, History*>::iterator it = selectFrom.begin();
2048 it != selectFrom.end(); ++it ) {
2050 if (it->second->sumScalarPT < sumMin) {
2051 sumMin = it->second->sumScalarPT;
2056 return selectFrom.lower_bound(iMin)->second;
2060 return selectFrom.upper_bound(sum*rnd)->second;
2062 return selectFrom.lower_bound(sum*rnd)->second;
2072 bool History::trimHistories() {
2074 if ( paths.empty() )
return false;
2076 for ( map<double, History*>::iterator it = paths.begin();
2077 it != paths.end(); ++it ) {
2079 if ( it->second->keep() && !it->second->keepHistory() )
2080 it->second->remove();
2083 double sumold, sumnew, sumprob, mismatch;
2084 sumold = sumnew = sumprob = mismatch = 0.;
2087 for ( map<double, History*>::iterator it = paths.begin();
2088 it != paths.end(); ++it ) {
2091 if ( it->second->keep() ) {
2093 goodBranches.insert( make_pair( sumnew - mismatch, it->second) );
2095 sumGoodBranches = sumnew - mismatch;
2099 double mismatchOld = mismatch;
2100 mismatch += sumnew - sumold;
2102 badBranches.insert( make_pair( mismatchOld + sumnew - sumold,
2105 sumBadBranches = mismatchOld + sumnew - sumold;
2113 return !goodBranches.empty();
2120 bool History::keepHistory() {
2121 bool keepPath =
true;
2124 if ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
2125 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
2126 || isQCD2to2(state) ) {
2129 double maxScale = hardFacScale(state);
2130 return keepPath = isOrderedPath( maxScale );
2134 if (isEW2to1(state)) {
2136 for (
int i = 0;i < state.size(); ++i)
2137 if (state[i].isFinal()) pSum += state[i].p();
2138 return isOrderedPath( pSum.mCalc());
2141 keepPath = isOrderedPath( infoPtr->eCM() );
2154 bool History::isOrderedPath(
double maxscale ) {
2155 double newscale = clusterIn.pT();
2156 if ( !mother )
return true;
2157 if ( mother->state[clusterIn.emittor].idAbs() == 21
2158 && mother->state[clusterIn.emitted].idAbs() == 5
2159 && !mother->state[clusterIn.emittor].isFinal())
2161 bool ordered = mother->isOrderedPath(newscale);
2162 if ( !ordered || maxscale < newscale)
return false;
2171 bool History::allIntermediateAboveRhoMS(
double rhoms,
bool good ) {
2174 if ( !good )
return false;
2177 for (
int i = 0; i < state.size(); ++i )
2178 if ( state[i].isFinal() && state[i].colType() != 0 )
2180 double rhoNew = (nFinal > 0 ) ? mergingHooksPtr->tmsNow( state )
2183 if ( !mother )
return good;
2185 return good && mother->allIntermediateAboveRhoMS( rhoms, (rhoNew > rhoms) );
2192 bool History::foundAnyOrderedPaths() {
2194 if ( paths.empty() )
return false;
2195 double maxscale = infoPtr->eCM();
2197 for ( map<double, History*>::iterator it = paths.begin();
2198 it != paths.end(); ++it )
2199 if ( it->second->isOrderedPath(maxscale) )
2220 double History::weightTree(PartonLevel* trial,
double as0,
double aem0,
2221 double maxscale,
double pdfScale, AlphaStrong * asFSR, AlphaStrong * asISR,
2222 AlphaEM * aemFSR, AlphaEM * aemISR,
double& asWeight,
double& aemWeight,
2223 double& pdfWeight) {
2226 double newScale = scale;
2231 int sideRad = (state[3].pz() > 0) ? 1 :-1;
2232 int sideRec = (state[4].pz() > 0) ? 1 :-1;
2235 if (state[3].colType() != 0) {
2237 double x = 2.*state[3].e() / state[0].e();
2238 int flav = state[3].id();
2240 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2241 double scaleDen = mergingHooksPtr->muFinME();
2243 double ratio = getPDFratio(sideRad,
false,
false, flav, x, scaleNum,
2249 if (state[4].colType() != 0) {
2251 double x = 2.*state[4].e() / state[0].e();
2252 int flav = state[4].id();
2254 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2255 double scaleDen = mergingHooksPtr->muFinME();
2257 double ratio = getPDFratio(sideRec,
false,
false, flav, x, scaleNum,
2267 double newPDFscale = newScale;
2268 if (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2269 newPDFscale = clusterIn.pT();
2272 double w = mother->weightTree(trial, as0, aem0, newScale, newPDFscale,
2273 asFSR, asISR, aemFSR, aemISR, asWeight, aemWeight, pdfWeight);
2276 if (state.size() < 3)
return 1.0;
2278 if ( w < 1e-12 )
return 0.0;
2280 w *= doTrialShower(trial, 1, maxscale);
2281 if ( w < 1e-12 )
return 0.0;
2283 int emtType = mother->state[clusterIn.emitted].colType();
2285 if ( asFSR && asISR && emtType != 0) {
2286 double asScale = pow2( newScale );
2287 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2288 asScale = pow2( clusterIn.pT() );
2291 bool FSR = mother->state[clusterIn.emittor].isFinal();
2292 if (!FSR) asScale += pow2(mergingHooksPtr->pT0ISR());
2295 if (mergingHooksPtr->useShowerPlugin() )
2296 asScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2297 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale);
2299 double alphaSinPS = (FSR) ? (*asFSR).alphaS(asScale)
2300 : (*asISR).alphaS(asScale);
2301 asWeight *= alphaSinPS / as0;
2305 if ( aemFSR && aemISR && emtType == 0 ) {
2306 double aemScale = pow2( newScale );
2307 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2308 aemScale = pow2( clusterIn.pT() );
2311 bool FSR = mother->state[clusterIn.emittor].isFinal();
2312 if (!FSR) aemScale += pow2(mergingHooksPtr->pT0ISR());
2315 if (mergingHooksPtr->useShowerPlugin() )
2316 aemScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2317 clusterIn.emitted, clusterIn.recoiler,
"scaleEM", aemScale);
2319 double alphaEMinPS = (FSR) ? (*aemFSR).alphaEM(aemScale)
2320 : (*aemISR).alphaEM(aemScale);
2321 aemWeight *= alphaEMinPS / aem0;
2327 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
2328 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
2330 if ( mother->state[inP].colType() != 0 ) {
2332 double x = getCurrentX(sideP);
2333 int flav = getCurrentFlav(sideP);
2335 double scaleNum = (children.empty())
2336 ? hardFacScale(state)
2337 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2338 ? pdfScale : maxscale );
2339 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2340 ? clusterIn.pT() : newScale;
2342 double ratio = getPDFratio(sideP,
false,
false, flav, x, scaleNum,
2347 if ( mother->state[inM].colType() != 0 ) {
2349 double x = getCurrentX(sideM);
2350 int flav = getCurrentFlav(sideM);
2352 double scaleNum = (children.empty())
2353 ? hardFacScale(state)
2354 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2355 ? pdfScale : maxscale );
2356 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2357 ? clusterIn.pT() : newScale;
2359 double ratio = getPDFratio(sideM,
false,
false, flav, x, scaleNum,
2372 double History::weightTreeALPHAS(
double as0, AlphaStrong * asFSR,
2373 AlphaStrong * asISR,
int njetMax ) {
2376 if ( !mother )
return 1.;
2378 double w = mother->weightTreeALPHAS( as0, asFSR, asISR, njetMax );
2380 if (state.size() < 3)
return w;
2383 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2384 if (njetNow >= njetMax)
return 1.0;
2387 bool FSR = mother->state[clusterIn.emittor].isFinal();
2388 int emtID = mother->state[clusterIn.emitted].id();
2391 if (abs(emtID) == 22 || abs(emtID) == 23 || abs(emtID) == 24)
return w;
2394 if ( asFSR && asISR ) {
2395 double asScale = pow2( scale );
2396 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2397 asScale = pow2( clusterIn.pT() );
2400 if (!FSR) asScale += pow2(mergingHooksPtr->pT0ISR());
2403 if (mergingHooksPtr->useShowerPlugin() )
2404 asScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2405 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale);
2407 double alphaSinPS = (FSR) ? (*asFSR).alphaS(asScale)
2408 : (*asISR).alphaS(asScale);
2409 w *= alphaSinPS / as0;
2420 double History::weightTreeALPHAEM(
double aem0, AlphaEM * aemFSR,
2421 AlphaEM * aemISR,
int njetMax ) {
2424 if ( !mother )
return 1.;
2426 double w = mother->weightTreeALPHAEM( aem0, aemFSR, aemISR, njetMax );
2428 if (state.size() < 3)
return w;
2431 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2432 if (njetNow >= njetMax)
return 1.0;
2435 bool FSR = mother->state[clusterIn.emittor].isFinal();
2436 int emtID = mother->state[clusterIn.emitted].id();
2439 if (!(abs(emtID) == 22 || abs(emtID) == 23 || abs(emtID) == 24))
return w;
2442 if ( aemFSR && aemISR ) {
2443 double aemScale = pow2( scale );
2444 if (mergingHooksPtr->unorderedASscalePrescip() == 1)
2445 aemScale = pow2( clusterIn.pT() );
2448 if (!FSR) aemScale += pow2(mergingHooksPtr->pT0ISR());
2451 if (mergingHooksPtr->useShowerPlugin() )
2452 aemScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2453 clusterIn.emitted, clusterIn.recoiler,
"scaleEM", aemScale);
2455 double alphaEMinPS = (FSR) ? (*aemFSR).alphaEM(aemScale)
2456 : (*aemISR).alphaEM(aemScale);
2457 w *= alphaEMinPS / aem0;
2468 double History::weightTreePDFs(
double maxscale,
double pdfScale,
2472 double newScale = scale;
2478 int njet = mergingHooksPtr->getNumberOfClusteringSteps( state);
2479 if (njet > njetMax)
return 1.0;
2482 int sideRad = (state[3].pz() > 0) ? 1 :-1;
2483 int sideRec = (state[4].pz() > 0) ? 1 :-1;
2486 if (state[3].colType() != 0) {
2488 double x = 2.*state[3].e() / state[0].e();
2489 int flav = state[3].id();
2491 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2492 double scaleDen = mergingHooksPtr->muFinME();
2494 wt *= getPDFratio(sideRad,
false,
false, flav, x, scaleNum, flav, x,
2499 if (state[4].colType() != 0) {
2501 double x = 2.*state[4].e() / state[0].e();
2502 int flav = state[4].id();
2504 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2505 double scaleDen = mergingHooksPtr->muFinME();
2507 wt *= getPDFratio(sideRec,
false,
false, flav, x, scaleNum, flav, x,
2516 double newPDFscale = newScale;
2517 if ( mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2518 newPDFscale = clusterIn.pT();
2521 double w = mother->weightTreePDFs( newScale, newPDFscale, njetMax );
2524 if (state.size() < 3)
return w;
2527 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2532 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
2533 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
2535 if ( mother->state[inP].colType() != 0 ) {
2537 double x = getCurrentX(sideP);
2538 int flav = getCurrentFlav(sideP);
2540 double scaleNum = (children.empty())
2541 ? hardFacScale(state)
2542 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2543 ? pdfScale : maxscale );
2544 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2545 ? clusterIn.pT() : newScale;
2547 double xDen = (njetNow == njetMax) ? mother->getCurrentX(sideP) : x;
2548 int flavDen = (njetNow == njetMax) ? mother->getCurrentFlav(sideP) : flav;
2549 double sDen = (njetNow == njetMax) ? mergingHooksPtr->muFinME() : scaleDen;
2550 double ratio = getPDFratio(sideP,
false,
false, flav, x, scaleNum,
2551 flavDen, xDen, sDen);
2556 if ( mother->state[inM].colType() != 0 ) {
2558 double x = getCurrentX(sideM);
2559 int flav = getCurrentFlav(sideM);
2561 double scaleNum = (children.empty())
2562 ? hardFacScale(state)
2563 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2564 ? pdfScale : maxscale );
2565 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2566 ? clusterIn.pT() : newScale;
2568 double xDen = (njetNow == njetMax) ? mother->getCurrentX(sideM) : x;
2569 int flavDen = (njetNow == njetMax) ? mother->getCurrentFlav(sideM) : flav;
2570 double sDen = (njetNow == njetMax) ? mergingHooksPtr->muFinME() : scaleDen;
2571 double ratio = getPDFratio(sideM,
false,
false, flav, x, scaleNum,
2572 flavDen, xDen, sDen);
2584 double History::weightTreeEmissions( PartonLevel* trial,
int type,
2585 int njetMin,
int njetMax,
double maxscale ) {
2588 double newScale = scale;
2591 if ( !mother )
return 1.0;
2593 double w = mother->weightTreeEmissions(trial,type,njetMin,njetMax,newScale);
2595 if (state.size() < 3)
return 1.0;
2597 if ( w < 1e-12 )
return 0.0;
2599 int njetNow = mergingHooksPtr->getNumberOfClusteringSteps( state) ;
2600 if (njetNow >= njetMax)
return 1.0;
2601 if (njetNow < njetMin ) w *= 1.0;
2603 else w *= doTrialShower(trial, type, maxscale);
2605 if ( w < 1e-12 )
return 0.0;
2615 double History::weightFirst(PartonLevel* trial,
double as0,
double muR,
2616 double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR, Rndm* rndmPtr ) {
2619 double newScale = scale;
2626 if (state[3].colType() != 0) {
2628 double x = 2.*state[3].e() / state[0].e();
2629 int flav = state[3].id();
2631 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2632 double scaleDen = mergingHooksPtr->muFinME();
2634 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2635 mergingHooksPtr->muFinME(), as0, rndmPtr);
2640 if (state[4].colType() != 0) {
2642 double x = 2.*state[4].e() / state[0].e();
2643 int flav = state[4].id();
2645 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2646 double scaleDen = mergingHooksPtr->muFinME();
2648 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2649 mergingHooksPtr->muFinME(), as0, rndmPtr);
2657 double w = mother->weightFirst(trial, as0, muR, newScale, asFSR, asISR,
2661 if (state.size() < 3)
return 0.0;
2665 double asScale2 = newScale*newScale;
2666 int showerType = (mother->state[clusterIn.emittor].isFinal() ) ? 1 : -1;
2667 if (showerType == -1) {
2668 asScale2 += pow(mergingHooksPtr->pT0ISR(),2);
2673 if (mergingHooksPtr->useShowerPlugin() ){
2674 asScale2 = getShowerPluginScale(mother->state, clusterIn.emittor,
2675 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale2);
2681 double BETA0 = 11. - 2./3.* NF;
2683 w += as0 / (2.*M_PI) * 0.5 * BETA0 * log( (muR*muR) / (b*asScale2) );
2689 double nWeight1 = 0.;
2690 double nWeight2 = 0.;
2692 for(
int i=0; i < NTRIAL; ++i) {
2694 vector<double> unresolvedEmissionTerm = countEmissions(trial, maxscale,
2695 newScale, 2, as0, asFSR, asISR, 3, fixpdf, fixas);
2696 nWeight1 += unresolvedEmissionTerm[1];
2698 w += nWeight1/double(NTRIAL) + nWeight2/double(NTRIAL);
2703 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
2704 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
2706 if ( mother->state[inP].colType() != 0 ) {
2708 double x = getCurrentX(sideP);
2709 int flav = getCurrentFlav(sideP);
2711 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2713 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, newScale,
2714 mergingHooksPtr->muFinME(), as0, rndmPtr);
2719 if ( mother->state[inM].colType() != 0 ) {
2721 double x = getCurrentX(sideM);
2722 int flav = getCurrentFlav(sideM);
2724 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2726 double intPDF4 = monteCarloPDFratios(flav, x, scaleNum, newScale,
2727 mergingHooksPtr->muFinME(), as0, rndmPtr);
2742 double History::weightFirstALPHAS(
double as0,
double muR,
2743 AlphaStrong * asFSR, AlphaStrong * asISR ) {
2746 double newScale = scale;
2748 if ( !mother )
return 0.;
2750 double w = mother->weightFirstALPHAS( as0, muR, asFSR, asISR );
2752 int showerType = (mother->state[clusterIn.emittor].isFinal() ) ? 1 : -1;
2754 double asScale = pow2( newScale );
2755 if ( mergingHooksPtr->unorderedASscalePrescip() == 1 )
2756 asScale = pow2( clusterIn.pT() );
2757 if (showerType == -1) {
2758 asScale += pow2( mergingHooksPtr->pT0ISR() );
2763 if (mergingHooksPtr->useShowerPlugin() ) {
2764 asScale = getShowerPluginScale(mother->state, clusterIn.emittor,
2765 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale);
2771 double BETA0 = 11. - 2./3.* NF;
2773 w += as0 / (2.*M_PI) * 0.5 * BETA0 * log( (muR*muR) / (b*asScale) );
2785 double History::weightFirstPDFs(
double as0,
double maxscale,
double pdfScale,
2789 double newScale = scale;
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 wt += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2805 mergingHooksPtr->muFinME(), as0, rndmPtr);
2808 if (state[4].colType() != 0) {
2810 double x = 2.*state[4].e() / state[0].e();
2811 int flav = state[4].id();
2813 double scaleNum = (children.empty()) ? hardFacScale(state) : maxscale;
2814 double scaleDen = mergingHooksPtr->muFinME();
2816 wt += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2817 mergingHooksPtr->muFinME(), as0, rndmPtr);
2826 double newPDFscale = newScale;
2827 if (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2828 newPDFscale = clusterIn.pT();
2831 double w = mother->weightFirstPDFs( as0, newScale, newPDFscale, rndmPtr);
2836 int sideP = (mother->state[inP].pz() > 0) ? 1 :-1;
2837 int sideM = (mother->state[inM].pz() > 0) ? 1 :-1;
2839 if ( mother->state[inP].colType() != 0 ) {
2841 double x = getCurrentX(sideP);
2842 int flav = getCurrentFlav(sideP);
2844 double scaleNum = (children.empty())
2845 ? hardFacScale(state)
2846 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2847 ? pdfScale : maxscale );
2848 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2849 ? clusterIn.pT() : newScale;
2851 w += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2852 mergingHooksPtr->muFinME(), as0, rndmPtr);
2855 if ( mother->state[inM].colType() != 0 ) {
2857 double x = getCurrentX(sideM);
2858 int flav = getCurrentFlav(sideM);
2860 double scaleNum = (children.empty())
2861 ? hardFacScale(state)
2862 : ( (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2863 ? pdfScale : maxscale );
2864 double scaleDen = (mergingHooksPtr->unorderedPDFscalePrescip() == 1)
2865 ? clusterIn.pT() : newScale;
2867 w += monteCarloPDFratios(flav, x, scaleNum, scaleDen,
2868 mergingHooksPtr->muFinME(), as0, rndmPtr);
2882 double History::weightFirstEmissions(PartonLevel* trial,
double as0,
2883 double maxscale, AlphaStrong * asFSR, AlphaStrong * asISR,
2884 bool fixpdf,
bool fixas ) {
2887 double newScale = scale;
2888 if ( !mother )
return 0.0;
2890 double w = mother->weightFirstEmissions(trial, as0, newScale, asFSR, asISR,
2893 if (state.size() < 3)
return 0.0;
2895 double nWeight1 = 0.;
2896 double nWeight2 = 0.;
2897 for(
int i=0; i < NTRIAL; ++i) {
2899 vector<double> unresolvedEmissionTerm = countEmissions(trial, maxscale,
2900 newScale, 2, as0, asFSR, asISR, 3, fixpdf, fixas);
2901 nWeight1 += unresolvedEmissionTerm[1];
2904 w += nWeight1/double(NTRIAL) + nWeight2/double(NTRIAL);
2915 double History::hardFacScale(
const Event& event) {
2917 double hardscale = 0.;
2919 if ( !mergingHooksPtr->resetHardQFac() )
return mergingHooksPtr->muF();
2923 if ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
2924 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
2925 || isQCD2to2(event)) {
2928 for (
int i=0; i <
event.size(); ++i)
2929 if ( event[i].isFinal() &&
event[i].colType() != 0 )
2930 mT.push_back( abs(event[i].mT2()) );
2931 if (
int(mT.size()) != 2 )
2932 hardscale = infoPtr->QFac();
2934 hardscale = sqrt( min( mT[0], mT[1] ) );
2936 hardscale = mergingHooksPtr->muF();
2946 double History::hardRenScale(
const Event& event) {
2948 double hardscale = 0.;
2950 if ( !mergingHooksPtr->resetHardQRen() )
return mergingHooksPtr->muR();
2954 if ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
2955 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
2956 || isQCD2to2(event)) {
2959 for (
int i=0; i <
event.size(); ++i)
2960 if ( event[i].isFinal()
2961 && (
event[i].colType() != 0 ||
event[i].id() == 22 ) )
2962 mT.push_back( abs(event[i].mT()) );
2963 if (
int(mT.size()) != 2 )
2964 hardscale = infoPtr->QRen();
2966 hardscale = sqrt( mT[0]*mT[1] );
2968 hardscale = mergingHooksPtr->muR();
2985 double History::doTrialShower( PartonLevel* trial,
int type,
2986 double maxscaleIn,
double minscaleIn ) {
2989 Event process = state;
2991 double startingScale = maxscaleIn;
2994 if ( mergingHooksPtr->getNumberOfClusteringSteps(process) == 0
2995 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
2996 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
2997 || isQCD2to2(state) ) )
2998 startingScale = min( startingScale, hardFacScale(process) );
3001 bool doVeto =
false;
3003 bool canEnhanceTrial = (trial->userHooksPtr!=0)
3004 && trial->userHooksPtr->canEnhanceTrial();
3009 trial->resetTrial();
3012 event.init(
"(hard process-modified)", particleDataPtr);
3016 process.scale(startingScale);
3020 double minScale = (minscaleIn > 0.) ? minscaleIn : scale;
3025 if (minScale >= startingScale)
break;
3031 double z = ( mergingHooksPtr->getNumberOfClusteringSteps(state) == 0
3034 : mother->getCurrentZ(clusterIn.emittor,clusterIn.recoiler,
3035 clusterIn.emitted, clusterIn.flavRadBef);
3037 infoPtr->zNowISR(z);
3038 infoPtr->pT2NowISR(pow(startingScale,2));
3039 infoPtr->hasHistory(
true);
3042 if (mergingHooksPtr->doWeakClustering()) setupWeakShower(0);
3045 trial->next(process,event);
3047 double pTtrial = trial->pTLastInShower();
3048 int typeTrial = trial->typeLastInShower();
3051 trial->resetTrial();
3054 double pTEnhanced = (canEnhanceTrial)
3055 ? trial->userHooksPtr->getEnhancedTrialPT() : 0.;
3056 double wtEnhanced = (canEnhanceTrial)
3057 ? trial->userHooksPtr->getEnhancedTrialWeight() : 1.;
3058 if ( canEnhanceTrial && pTEnhanced > 0.) pTtrial = pTEnhanced;
3061 double vetoScale = (mother) ? 0. : mergingHooksPtr->tms();
3063 double tnow = mergingHooksPtr->tmsNow( event );
3066 if ( pTtrial < minScale )
break;
3068 startingScale = pTtrial;
3071 if ( tnow < vetoScale && vetoScale > 0. )
continue;
3074 if ( mergingHooksPtr->canVetoTrialEmission()
3075 && mergingHooksPtr->doVetoTrialEmission( process, event) )
continue;
3077 int iRecAft =
event.size() - 1;
3078 int iEmt =
event.size() - 2;
3079 int iRadAft =
event.size() - 3;
3080 if ( (event[iRecAft].status() != 52 && event[iRecAft].status() != -53) ||
3081 event[iEmt].status() != 51 || event[iRadAft].status() != 51)
3082 iRecAft = iEmt = iRadAft = -1;
3083 for (
int i = event.size() - 1; i > 0; i--) {
3084 if (iRadAft == -1 && event[i].status() == -41) iRadAft = i;
3085 else if (iEmt == -1 && event[i].status() == 43) iEmt = i;
3086 else if (iRecAft == -1 && event[i].status() == -42) iRecAft = i;
3087 if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
break;
3092 if ( type == -1 && typeTrial != 1 )
continue;
3094 if ( type == 1 && !(typeTrial == 2 || typeTrial >= 3) )
continue;
3097 if (canEnhanceTrial && pTtrial > minScale) wt *= (1. - 1./wtEnhanced);
3099 if ( canEnhanceTrial && wt == 0.)
break;
3101 if ( canEnhanceTrial && pTtrial > minScale)
continue;
3105 if ( pTtrial > minScale ) doVeto =
true;
3115 && mergingHooksPtr->getNumberOfClusteringSteps(process) == 0
3116 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
3117 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
3118 || isQCD2to2(state))
3119 && pTtrial > hardFacScale(process) )
3124 if ( pTtrial < minScale ) doVeto =
false;
3132 double res = (canEnhanceTrial) ? wt : ( (doVeto) ? 0. : 1. );
3144 bool History::updateind(vector<int> & ind,
int i,
int N) {
3145 if ( i < 0 )
return false;
3146 if ( ++ind[i] < N )
return true;
3147 if ( !updateind(ind, i - 1, N - 1) )
return false;
3148 ind[i] = ind[i - 1] + 1;
3159 History::countEmissions(PartonLevel* trial,
double maxscale,
3160 double minscale,
int showerType,
double as0,
3161 AlphaStrong * asFSR, AlphaStrong * asISR,
int N = 1,
3162 bool fixpdf =
true,
bool fixas =
true) {
3164 if ( N < 0 )
return vector<double>();
3165 vector<double> result(N+1);
3167 if ( N < 1 )
return result;
3170 Event process = state;
3172 double startingScale = maxscale;
3175 if ( mergingHooksPtr->getNumberOfClusteringSteps(process) == 0
3176 && ( mergingHooksPtr->getProcessString().compare(
"pp>jj") == 0
3177 || mergingHooksPtr->getProcessString().compare(
"pp>aj") == 0
3178 || isQCD2to2(state) ) )
3179 startingScale = min( startingScale, hardFacScale(process) );
3182 bool canEnhanceTrial = (trial->userHooksPtr!=0)
3183 && trial->userHooksPtr->canEnhanceTrial();
3187 trial->resetTrial();
3190 event.init(
"(hard process-modified)", particleDataPtr);
3194 process.scale(startingScale);
3199 if (minscale >= startingScale)
return result;
3205 double z = ( mergingHooksPtr->getNumberOfClusteringSteps(state) == 0)
3207 : mother->getCurrentZ(clusterIn.emittor,clusterIn.recoiler,
3210 infoPtr->zNowISR(z);
3211 infoPtr->pT2NowISR(pow(startingScale,2));
3212 infoPtr->hasHistory(
true);
3216 if (mergingHooksPtr->doWeakClustering()) setupWeakShower(0);
3219 trial->next(process,event);
3222 double pTtrial = trial->pTLastInShower();
3223 int typeTrial = trial->typeLastInShower();
3226 trial->resetTrial();
3229 double pTEnhanced = (canEnhanceTrial)
3230 ? trial->userHooksPtr->getEnhancedTrialPT() : 0.;
3231 double wtEnhanced = (canEnhanceTrial)
3232 ? trial->userHooksPtr->getEnhancedTrialWeight() : 1.;
3233 if ( canEnhanceTrial && pTEnhanced > 0.) pTtrial = pTEnhanced;
3236 double vetoScale = (mother) ? 0. : mergingHooksPtr->tms();
3238 double tnow = mergingHooksPtr->tmsNow( event );
3241 startingScale = pTtrial;
3243 if ( pTtrial < minscale )
break;
3245 if ( tnow < vetoScale && vetoScale > 0. )
continue;
3247 if ( mergingHooksPtr->canVetoTrialEmission()
3248 && mergingHooksPtr->doVetoTrialEmission( process, event) )
continue;
3251 double enhance = (canEnhanceTrial && pTtrial > minscale) ? wtEnhanced : 1.;
3256 double alphaSinPS = as0;
3259 double asScale2 = pTtrial*pTtrial;
3261 if (mergingHooksPtr->useShowerPlugin() )
3262 asScale2 = getShowerPluginScale(mother->state, clusterIn.emittor,
3263 clusterIn.emitted, clusterIn.recoiler,
"scaleAS", asScale2);
3266 if ( (showerType == -1 || showerType == 2) && typeTrial == 2 ) {
3268 if ( fixas ) alphaSinPS = (*asISR).alphaS(asScale2);
3271 pdfs = pdfFactor( event, typeTrial, pTtrial,
3272 mergingHooksPtr->muFinME() );
3274 }
else if ( (showerType == 1 || showerType == 2) && typeTrial >= 3 ) {
3276 if ( fixas ) alphaSinPS = (*asFSR).alphaS(asScale2);
3280 pdfs = pdfFactor( event, typeTrial, pTtrial,
3281 mergingHooksPtr->muFinME() );
3285 if ( typeTrial == 2 || typeTrial >= 3 )
3286 wts.push_back(as0/alphaSinPS * pdfs * 1./enhance);
3290 for (
int n = 1; n <= min(N,
int(wts.size())); ++n ) {
3292 for (
int i = 0; i < N; ++i ) ind[i] = i;
3295 for (
int j = 0; j < n; ++j ) x *= wts[ind[j]];
3297 }
while ( updateind(ind, n - 1, wts.size()) );
3298 if ( n%2 ) result[n] *= -1.0;
3310 double History::monteCarloPDFratios(
int flav,
double x,
double maxScale,
3311 double minScale,
double pdfScale,
double asME, Rndm* rndmPtr) {
3315 double factor = asME / (2.*M_PI);
3317 factor *= log(maxScale/minScale);
3320 if (factor == 0.)
return 0.;
3328 double integral = 0.;
3329 double RN = rndmPtr->flat();
3332 double zTrial = pow(x,RN);
3333 integral = -log(x) * zTrial *
3334 integrand(flav, x, pdfScale, zTrial);
3335 integral += 1./6.*(11.*CA - 4.*NF*TR)
3338 double zTrial = x + RN*(1. - x);
3340 integrand(flav, x, pdfScale, zTrial);
3341 integral += 3./2.*CF
3346 return (factor*integral);
3355 bool History::onlyOrderedPaths() {
3356 if ( !mother || foundOrderedPath )
return foundOrderedPath;
3357 return foundOrderedPath = mother->onlyOrderedPaths();
3366 bool History::onlyStronglyOrderedPaths() {
3367 if ( !mother || foundStronglyOrderedPath )
return foundStronglyOrderedPath;
3368 return foundStronglyOrderedPath = mother->onlyStronglyOrderedPaths();
3377 bool History::onlyAllowedPaths() {
3378 if ( !mother || foundAllowedPath )
return foundAllowedPath;
3379 return foundAllowedPath = mother->onlyAllowedPaths();
3391 bool History::registerPath(History & l,
bool isOrdered,
3392 bool isStronglyOrdered,
bool isAllowed,
bool isComplete) {
3398 if ( mother )
return mother->registerPath(l, isOrdered,
3399 isStronglyOrdered, isAllowed, isComplete);
3402 if ( sumpath == sumpath + l.prob )
3404 if ( mergingHooksPtr->canCutOnRecState()
3405 && foundAllowedPath && !isAllowed )
3407 if ( mergingHooksPtr->enforceStrongOrdering()
3408 && foundStronglyOrderedPath && !isStronglyOrdered )
3410 if ( mergingHooksPtr->orderHistories()
3411 && foundOrderedPath && !isOrdered ) {
3413 if ( (!foundCompletePath && isComplete)
3414 || (!foundAllowedPath && isAllowed) ) ;
3418 if ( foundCompletePath && !isComplete)
3420 if ( !mergingHooksPtr->canCutOnRecState()
3421 && !mergingHooksPtr->allowCutOnRecState() )
3422 foundAllowedPath =
true;
3424 if ( mergingHooksPtr->canCutOnRecState() && isAllowed && isComplete) {
3425 if ( !foundAllowedPath || !foundCompletePath ) {
3431 foundAllowedPath =
true;
3435 if ( mergingHooksPtr->enforceStrongOrdering() && isStronglyOrdered
3437 if ( !foundStronglyOrderedPath || !foundCompletePath ) {
3443 foundStronglyOrderedPath =
true;
3444 foundCompletePath =
true;
3448 if ( mergingHooksPtr->orderHistories() && isOrdered && isComplete ) {
3449 if ( !foundOrderedPath || !foundCompletePath ) {
3455 foundOrderedPath =
true;
3456 foundCompletePath =
true;
3461 if ( !foundCompletePath ) {
3467 foundCompletePath =
true;
3472 foundOrderedPath =
true;
3476 double weakProb = 1.;
3477 if (mergingHooksPtr->doWeakClustering())
3478 weakProb = l.getWeakProb();
3481 sumpath += l.prob * weakProb;
3482 paths[sumpath] = &l;
3494 vector<Clustering> History::getAllQCDClusterings() {
3495 vector<Clustering> ret;
3498 vector <int> posFinalPartn;
3499 vector <int> posInitPartn;
3500 vector <int> posFinalGluon;
3501 vector <int> posFinalQuark;
3502 vector <int> posFinalAntiq;
3503 vector <int> posInitGluon;
3504 vector <int> posInitQuark;
3505 vector <int> posInitAntiq;
3509 for (
int i=0; i < state.size(); ++i )
3510 if ( state[i].isFinal() && state[i].colType() !=0 ) {
3512 if ( state[i].
id() == 21 ) posFinalGluon.push_back(i);
3513 else if ( state[i].idAbs() < 10 && state[i].
id() > 0)
3514 posFinalQuark.push_back(i);
3515 else if ( state[i].idAbs() < 10 && state[i].
id() < 0)
3516 posFinalAntiq.push_back(i);
3517 }
else if (state[i].status() == -21 && state[i].colType() != 0 ) {
3519 if ( state[i].
id() == 21 ) posInitGluon.push_back(i);
3520 else if ( state[i].idAbs() < 10 && state[i].
id() > 0)
3521 posInitQuark.push_back(i);
3522 else if ( state[i].idAbs() < 10 && state[i].
id() < 0)
3523 posInitAntiq.push_back(i);
3527 vector<Clustering> systems;
3528 systems = getQCDClusterings(state);
3529 ret.insert(ret.end(), systems.begin(), systems.end());
3533 if ( !ret.empty() )
return ret;
3538 else if ( ret.empty()
3539 && mergingHooksPtr->allowColourShuffling() ) {
3542 for(
int i = 0; i < int(posFinalQuark.size()); ++i) {
3544 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(posFinalQuark[i],
3547 int col = NewState[posFinalQuark[i]].col();
3548 for(
int j = 0; j < int(posInitAntiq.size()); ++j) {
3550 int acl = NewState[posInitAntiq[j]].acol();
3551 if ( col == acl )
continue;
3552 NewState[posFinalQuark[i]].col(acl);
3553 NewState[posInitAntiq[j]].acol(col);
3554 systems = getQCDClusterings(NewState);
3555 if (!systems.empty()) {
3558 ret.insert(ret.end(), systems.begin(), systems.end());
3565 for(
int i = 0; i < int(posFinalAntiq.size()); ++i) {
3567 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(posFinalAntiq[i],
3570 int acl = NewState[posFinalAntiq[i]].acol();
3571 for(
int j = 0; j < int(posInitQuark.size()); ++j) {
3573 int col = NewState[posInitQuark[j]].col();
3574 if ( col == acl )
continue;
3575 NewState[posFinalAntiq[i]].acol(col);
3576 NewState[posInitQuark[j]].col(acl);
3577 systems = getQCDClusterings(NewState);
3578 if (!systems.empty()) {
3581 ret.insert(ret.end(), systems.begin(), systems.end());
3588 if ( !ret.empty() ) {
3589 string message=
"Warning in History::getAllQCDClusterings: Changed";
3590 message+=
" colour structure to allow at least one clustering.";
3591 infoPtr->errorMsg(message);
3606 vector<Clustering> History::getQCDClusterings(
const Event& event) {
3607 vector<Clustering> ret;
3611 vector <int> posFinalPartn;
3612 vector <int> posInitPartn;
3614 vector <int> posFinalGluon;
3615 vector <int> posFinalQuark;
3616 vector <int> posFinalAntiq;
3617 vector <int> posInitGluon;
3618 vector <int> posInitQuark;
3619 vector <int> posInitAntiq;
3623 for (
int i=0; i <
event.size(); ++i)
3624 if ( event[i].isFinal() &&
event[i].colType() !=0 ) {
3626 posFinalPartn.push_back(i);
3627 if ( event[i].
id() == 21 ) posFinalGluon.push_back(i);
3628 else if ( event[i].idAbs() < 10 && event[i].
id() > 0)
3629 posFinalQuark.push_back(i);
3630 else if ( event[i].idAbs() < 10 && event[i].
id() < 0)
3631 posFinalAntiq.push_back(i);
3632 }
else if ( event[i].status() == -21 && event[i].colType() != 0 ) {
3634 posInitPartn.push_back(i);
3635 if ( event[i].
id() == 21 ) posInitGluon.push_back(i);
3636 else if ( event[i].idAbs() < 10 && event[i].
id() > 0)
3637 posInitQuark.push_back(i);
3638 else if ( event[i].idAbs() < 10 && event[i].
id() < 0)
3639 posInitAntiq.push_back(i);
3642 int nFiGluon = int(posFinalGluon.size());
3643 int nFiQuark = int(posFinalQuark.size());
3644 int nFiAntiq = int(posFinalAntiq.size());
3645 int nInGluon = int(posInitGluon.size());
3646 int nInQuark = int(posInitQuark.size());
3647 int nInAntiq = int(posInitAntiq.size());
3649 vector<Clustering> systems;
3653 for (
int i = 0; i < nFiGluon; ++i) {
3654 int EmtGluon = posFinalGluon[i];
3655 systems = findQCDTriple( EmtGluon, 2, event, posFinalPartn, posInitPartn);
3656 ret.insert(ret.end(), systems.begin(), systems.end());
3662 bool check_g2qq =
true;
3663 if ( ( ( nInQuark + nInAntiq == 0 )
3665 && (nFiQuark == 1) && (nFiAntiq == 1) )
3666 || ( ( nFiQuark + nFiAntiq == 0)
3667 && (nInQuark == 1) && (nInAntiq == 1) ) )
3674 for(
int i=0; i < nFiQuark; ++i) {
3675 int EmtQuark = posFinalQuark[i];
3676 systems = findQCDTriple( EmtQuark,1,event, posFinalPartn, posInitPartn);
3677 ret.insert(ret.end(), systems.begin(), systems.end());
3683 for(
int i=0; i < nFiAntiq; ++i) {
3684 int EmtAntiq = posFinalAntiq[i];
3685 systems = findQCDTriple( EmtAntiq,1,event, posFinalPartn, posInitPartn);
3686 ret.insert(ret.end(), systems.begin(), systems.end());
3698 void History::attachClusterings (vector<Clustering>& clus,
int iEmt,
int iRad,
3699 int iRec,
int iPartner,
double pT,
const Event& event) {
3701 if ( !mergingHooksPtr->doWeakClustering() ) {
3703 clus.push_back( Clustering(iEmt, iRad, iRec, iPartner,
3704 pT, 0, 0, 0, 0, 9));
3709 int radSpin =
event[iRad].intPol();
3710 int emtSpin =
event[iEmt].intPol();
3711 int recSpin =
event[iRec].intPol();
3712 bool hasRadSpin = (radSpin != 9);
3713 bool hasEmtSpin = (emtSpin != 9);
3714 bool hasRecSpin = (recSpin != 9);
3717 bool radQuark =
event[iRad].idAbs() < 10;
3718 bool emtQuark =
event[iEmt].idAbs() < 10;
3719 bool recQuark =
event[iRec].idAbs() < 10;
3722 vector < vector<int> > structs;
3724 for (
int i = 0; i < 3; ++i){
3725 int sRad = (i==0) ? -1 : (i==1) ? 1 : 9;
3726 for (
int j = 0; j < 3; ++j){
3727 int sEmt = (j==0) ? -1 : (j==1) ? 1 : 9;
3728 for (
int k = 0; k < 3; ++k){
3729 int sRec = (k==0) ? -1 : (k==1) ? 1 : 9;
3731 s.push_back(sRad); s.push_back(sEmt); s.push_back(sRec);
3732 structs.push_back(s);
3737 vector < vector<int> > allStructs;
3738 for (
int i = 0; i< int(structs.size()); ++i) {
3740 if (hasRadSpin && radQuark && structs[i][0] != radSpin)
continue;
3741 if (hasEmtSpin && emtQuark && structs[i][1] != emtSpin)
continue;
3742 if (hasRecSpin && recQuark && structs[i][2] != recSpin)
continue;
3744 if (!hasRadSpin && radQuark && structs[i][0] == 9)
continue;
3745 if (!hasEmtSpin && emtQuark && structs[i][1] == 9)
continue;
3746 if (!hasRecSpin && recQuark && structs[i][2] == 9)
continue;
3748 if (!radQuark && structs[i][0] != radSpin)
continue;
3749 if (!emtQuark && structs[i][1] != emtSpin)
continue;
3750 if (!recQuark && structs[i][2] != recSpin)
continue;
3752 if (radQuark && emtQuark && structs[i][0] != structs[i][1])
continue;
3754 allStructs.push_back(structs[i]);
3758 int flavRadBef = getRadBeforeFlav(iRad, iEmt, event);
3760 for (
int i = 0; i< int(allStructs.size()); ++i) {
3762 int spinRadBef = getRadBeforeSpin(iRad, iEmt, allStructs[i][0],
3763 allStructs[i][1], event);
3764 clus.push_back( Clustering(iEmt, iRad, iRec, iPartner, pT, flavRadBef,
3765 allStructs[i][0], allStructs[i][1], allStructs[i][2], spinRadBef) );
3786 vector<Clustering> History::findQCDTriple (
int EmtTagIn,
int colTopIn,
3788 vector<int> posFinalPartn,
3789 vector <int> posInitPartn ) {
3792 int EmtTag = EmtTagIn;
3795 int colTop = colTopIn;
3798 int finalSize = int(posFinalPartn.size());
3799 int initSize = int(posInitPartn.size());
3800 int size = initSize + finalSize;
3802 vector<Clustering> clus;
3806 for (
int a = 0; a < size; ++a ) {
3807 int i = (a < finalSize)? a : (a - finalSize) ;
3808 int iRad = (a < finalSize)? posFinalPartn[i] : posInitPartn[i];
3810 if ( event[iRad].col() ==
event[EmtTag].col()
3811 &&
event[iRad].acol() ==
event[EmtTag].acol() )
3814 if (iRad != EmtTag ) {
3815 int pTdef =
event[iRad].isFinal() ? 1 : -1;
3816 int sign = (a < finalSize)? 1 : -1 ;
3822 if ( event[iRad].
id() == -sign*
event[EmtTag].id() ) {
3826 if (event[iRad].isFinal() ) {
3827 if (event[iRad].
id() < 0) {
3828 col =
event[EmtTag].col();
3831 acl =
event[EmtTag].acol();
3835 if (event[iRad].
id() < 0) {
3836 acl =
event[EmtTag].acol();
3839 col =
event[EmtTag].col();
3851 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
3854 if ( (sign < 0) && (event[iRec].isFinal()) ) {
3858 for(
int l = 0; l < int(posInitPartn.size()); ++l)
3859 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
3865 if ( iRec != 0 && iPartner != 0
3866 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
3867 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
3868 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
3875 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
3878 if ( (sign < 0) && (event[iRec].isFinal()) ) {
3882 for(
int l = 0; l < int(posInitPartn.size()); ++l)
3883 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
3889 if ( iRec != 0 && iPartner != 0
3890 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
3891 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
3892 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
3903 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
3906 if ( (sign < 0) && (event[iRec].isFinal()) ) {
3910 for(
int l = 0; l < int(posInitPartn.size()); ++l)
3911 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
3917 if ( iRec != 0 && iPartner != 0
3918 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
3919 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
3920 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
3927 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
3930 if ( (sign < 0) && (event[iRec].isFinal()) ) {
3934 for(
int l = 0; l < int(posInitPartn.size()); ++l)
3935 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
3941 if ( iRec != 0 && iPartner != 0
3942 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
3943 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
3944 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
3949 }
else if ( event[iRad].
id() == 21
3950 &&( event[iRad].col() == event[EmtTag].col()
3951 || event[iRad].acol() == event[EmtTag].acol() )) {
3957 for(
int l = 0; l < int(posInitPartn.size()); ++l)
3958 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
3962 int col = getRadBeforeCol(iRad, EmtTag, event);
3963 int acl = getRadBeforeAcol(iRad, EmtTag, event);
3973 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
3974 ? event[iRad].col() :
event[iRad].acol();
3977 if (colRemove > 0 && col > 0 && col != colRemove)
3978 iPartner = FindCol(col,iRad,EmtTag,event,1,
true)
3979 + FindCol(col,iRad,EmtTag,event,2,
true);
3980 else if (colRemove > 0 && acl > 0 && acl != colRemove)
3981 iPartner = FindCol(acl,iRad,EmtTag,event,1,
true)
3982 + FindCol(acl,iRad,EmtTag,event,2,
true);
3984 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event ) ) {
3985 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
3986 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
3994 if ( (event[iRad].col() == event[EmtTag].acol())
3995 || (event[iRad].acol() == event[EmtTag].col())
3996 || (event[iRad].col() == event[EmtTag].col())
3997 || (event[iRad].acol() == event[EmtTag].acol()) ) {
4004 if (event[iRad].isFinal() ) {
4006 if ( event[iRad].
id() < 0) {
4007 acl =
event[EmtTag].acol();
4008 col =
event[iRad].col();
4009 }
else if ( event[iRad].
id() > 0 && event[iRad].
id() < 10) {
4010 col =
event[EmtTag].col();
4011 acl =
event[iRad].acol();
4013 col =
event[EmtTag].col();
4014 acl =
event[EmtTag].acol();
4019 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
4020 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4022 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4023 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4024 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4028 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
4029 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4031 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4032 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4033 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4039 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
4040 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4042 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4043 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4044 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4048 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
4049 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4051 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4052 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4053 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4067 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4068 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
4072 col = getRadBeforeCol(iRad, EmtTag, event);
4073 acl = getRadBeforeAcol(iRad, EmtTag, event);
4083 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
4084 ? event[iRad].col() : 0;
4085 iPartner = (colRemove > 0)
4086 ? FindCol(col,iRad,EmtTag,event,1,
true)
4087 + FindCol(col,iRad,EmtTag,event,2,
true)
4088 : FindCol(acl,iRad,EmtTag,event,1,true)
4089 + FindCol(acl,iRad,EmtTag,event,2,true);
4091 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event)) {
4092 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
4093 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
4113 vector<Clustering> History::getAllEWClusterings() {
4114 vector<Clustering> ret;
4117 vector<Clustering> systems;
4118 systems = getEWClusterings(state);
4119 ret.insert(ret.end(), systems.begin(), systems.end());
4130 vector<Clustering> History::getEWClusterings(
const Event& event) {
4131 vector<Clustering> ret;
4135 vector <int> posFinalPartn;
4136 vector <int> posInitPartn;
4137 vector <int> posFinalW;
4138 vector <int> posFinalZ;
4142 for (
int i=3; i <
event.size(); ++i )
4143 if ( event[i].isFinal() ) {
4145 posFinalPartn.push_back(i);
4148 posInitPartn.push_back(i);
4151 for (
int i=0; i <
event.size(); ++i )
4152 if ( event[i].isFinal() &&
event[i].idAbs() == 24 )
4153 posFinalW.push_back( i );
4156 for (
int i=0; i <
event.size(); ++i )
4157 if ( event[i].isFinal() &&
event[i].idAbs() == 23 )
4158 posFinalZ.push_back( i );
4161 vector<Clustering> systems;
4164 for (
int i = 0; i < int(posFinalW.size()); ++i ) {
4165 int emtW = posFinalW[i];
4166 systems = findEWTripleW( emtW, event, posFinalPartn, posInitPartn);
4167 ret.insert(ret.end(), systems.begin(), systems.end());
4172 for (
int i = 0; i < int(posFinalZ.size()); ++i ) {
4173 int emtZ = posFinalZ[i];
4175 systems = findEWTripleZ( emtZ, event, posFinalPartn, posInitPartn);
4176 ret.insert(ret.end(), systems.begin(), systems.end());
4195 vector<Clustering> History::findEWTripleW (
int emtTagIn,
const Event& event,
4196 vector<int> posFinalPartn, vector<int> posInitPartn ) {
4198 int emtTag = emtTagIn;
4199 int flavEmt =
event[emtTag].id();
4205 int finalSize = int(posFinalPartn.size());
4206 int initSize = int(posInitPartn.size());
4209 vector<int> flavCounts(30,0);
4211 for (
int a = 0; a < finalSize; ++a ) {
4212 if (event[posFinalPartn[a]].idAbs() < 20) {
4214 if (event[posFinalPartn[a]].
id() < 0)
4216 flavCounts[
event[posFinalPartn[a]].idAbs()] += sign;
4218 if (event[posFinalPartn[a]].idAbs() == 24)
4222 for (
int a = 0; a < initSize; ++a ) {
4223 if (event[posInitPartn[a]].idAbs() < 20) {
4225 if (event[posInitPartn[a]].
id() < 0)
4227 flavCounts[
event[posInitPartn[a]].idAbs()] -= sign;
4231 vector<Clustering> clus;
4235 for (
int a = 0; a < finalSize; ++a ) {
4237 int iRad = posFinalPartn[a];
4238 if (iRad != emtTag) {
4241 int spinRad =
event[iRad].intPol();
4242 if (spinRad == -1 || spinRad == 9 || spinRad == 0) {
4246 int flavRad =
event[iRad].id();
4248 if (event[iRad].isQuark() || event[iRad].isLepton()) {
4251 int flavExp = (flavRad > 0) ? 24 : -24;
4252 if (abs(flavRad) % 2 == 0) flavExp = -flavExp;
4253 if (flavExp == flavEmt) {
4256 vector<int> flavRadBefs = posFlavCKM(flavRad);
4260 for (
int i = 0;i < int(flavRadBefs.size()); ++i)
4261 flavRadBefs[i] = - flavRadBefs[i];
4264 for (
int i = 0; i < finalSize; ++i ) {
4265 int iRec = posFinalPartn[i];
4268 if ( iRec != iRad && iRec != emtTag ) {
4270 for (
int j = 0;j < int(flavRadBefs.size()); ++j) {
4273 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4274 flavRad, flavRadBefs[j], 1))
4277 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4278 pTLund(event, iRad, emtTag, iRec, pTdef, flavRadBefs[j]),
4279 flavRadBefs[j], -1 ) );
4290 for (
int a = 0;a < int(posInitPartn.size()); ++a) {
4291 int iRad = posInitPartn[a];
4292 int flavRad =
event[iRad].id();
4294 if (event[iRad].isQuark() || event[iRad].isLepton()) {
4297 int spinRad =
event[iRad].intPol();
4298 if (spinRad == -1 || spinRad == 9 || spinRad == 0) {
4301 int flavExp = (flavRad > 0) ? 24 : -24;
4302 if (abs(flavRad) % 2 == 1) flavExp = -flavExp;
4303 if (flavExp == flavEmt) {
4305 vector<int> flavRadBefs = posFlavCKM(flavRad);
4309 for (
int i = 0;i < int(flavRadBefs.size()); ++i)
4310 flavRadBefs[i] = - flavRadBefs[i];
4313 for (
int i = 0; i < int(posInitPartn.size()); ++i ) {
4314 int iRec = posInitPartn[i];
4317 if ( i != a && iRec != emtTag) {
4318 for (
int j = 0;j < int(flavRadBefs.size()); ++j) {
4322 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4323 flavRad, flavRadBefs[j], -1))
4325 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4326 pTLund(event, iRad, emtTag, iRec, -1, flavRadBefs[j]),
4327 flavRadBefs[j], -1 ) );
4352 vector<Clustering> History::findEWTripleZ (
int emtTagIn,
const Event& event,
4353 vector<int> posFinalPartn, vector<int> posInitPartn ) {
4355 int emtTag = emtTagIn;
4361 int finalSize = int(posFinalPartn.size());
4362 int initSize = int(posInitPartn.size());
4365 vector<int> flavCounts(30,0);
4367 for (
int a = 0; a < finalSize; ++a ) {
4368 if (event[posFinalPartn[a]].idAbs() < 20) {
4370 if (event[posFinalPartn[a]].
id() < 0)
4372 flavCounts[
event[posFinalPartn[a]].idAbs()] += sign;
4374 if (event[posFinalPartn[a]].idAbs() == 24)
4378 for (
int a = 0; a < initSize; ++a ) {
4379 if (event[posInitPartn[a]].idAbs() < 20) {
4381 if (event[posInitPartn[a]].
id() < 0)
4383 flavCounts[
event[posInitPartn[a]].idAbs()] -= sign;
4387 vector<Clustering> clus;
4392 for (
int a = 0; a < finalSize; ++a ) {
4394 int iRad = posFinalPartn[a];
4395 if (iRad != emtTag) {
4398 int flavRad =
event[iRad].id();
4401 if (event[iRad].isQuark() || event[iRad].isLepton())
4403 for (
int i = 0; i < finalSize; ++i ) {
4404 int iRec = posFinalPartn[i];
4407 if ( iRec != iRad && iRec != emtTag ) {
4410 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4411 flavRad, flavRad, 1))
4414 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4415 pTLund(event, iRad, emtTag, iRec, pTdef, flavRad),
4423 for (
int a = 0;a < int(posInitPartn.size()); ++a) {
4424 int iRad = posInitPartn[a];
4425 int flavRad =
event[iRad].id();
4427 if (event[iRad].isQuark() || event[iRad].isLepton()) {
4430 for (
int i = 0; i < int(posInitPartn.size()); ++i ) {
4431 int iRec = posInitPartn[i];
4434 if ( i != a && iRec != emtTag) {
4437 if (flavCounts[24] <= 1 && !checkFlavour(flavCounts,
4438 flavRad, flavRad, -1))
4440 clus.push_back( Clustering(emtTag, iRad, iRec, iRec,
4441 pTLund(event, iRad, emtTag, iRec, -1, flavRad),
4459 vector<Clustering> History::getAllSQCDClusterings() {
4460 vector<Clustering> ret;
4463 vector<Clustering> systems;
4464 systems = getSQCDClusterings(state);
4465 ret.insert(ret.end(), systems.begin(), systems.end());
4476 vector<Clustering> History::getSQCDClusterings(
const Event& event) {
4477 vector<Clustering> ret;
4481 vector <int> posFinalPartn;
4482 vector <int> posInitPartn;
4484 vector <int> posFinalGluon;
4485 vector <int> posFinalQuark;
4486 vector <int> posFinalAntiq;
4487 vector <int> posInitGluon;
4488 vector <int> posInitQuark;
4489 vector <int> posInitAntiq;
4493 for (
int i=0; i <
event.size(); ++i)
4494 if ( event[i].isFinal() &&
event[i].colType() !=0 ) {
4496 posFinalPartn.push_back(i);
4497 if ( event[i].
id() == 21 || event[i].
id() == 1000021)
4498 posFinalGluon.push_back(i);
4499 else if ( (event[i].idAbs() < 10 && event[i].
id() > 0)
4500 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4501 && event[i].
id() > 0)
4502 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4503 && event[i].
id() > 0))
4504 posFinalQuark.push_back(i);
4505 else if ( (event[i].idAbs() < 10 && event[i].
id() < 0)
4506 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4507 && event[i].
id() < 0)
4508 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4509 && event[i].
id() < 0))
4510 posFinalAntiq.push_back(i);
4511 }
else if ( event[i].status() == -21 && event[i].colType() != 0 ) {
4513 posInitPartn.push_back(i);
4514 if ( event[i].
id() == 21 || event[i].
id() == 1000021)
4515 posInitGluon.push_back(i);
4516 else if ( (event[i].idAbs() < 10 && event[i].
id() > 0)
4517 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4518 && event[i].
id() > 0)
4519 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4520 && event[i].
id() > 0))
4521 posInitQuark.push_back(i);
4522 else if ( (event[i].idAbs() < 10 && event[i].
id() < 0)
4523 || (event[i].idAbs() < 1000010 && event[i].idAbs() > 1000000
4524 && event[i].
id() < 0)
4525 || (event[i].idAbs() < 2000010 && event[i].idAbs() > 2000000
4526 && event[i].
id() < 0))
4527 posInitAntiq.push_back(i);
4530 int nFiGluon = int(posFinalGluon.size());
4531 int nFiQuark = int(posFinalQuark.size());
4532 int nFiAntiq = int(posFinalAntiq.size());
4533 int nInGluon = int(posInitGluon.size());
4534 int nInQuark = int(posInitQuark.size());
4535 int nInAntiq = int(posInitAntiq.size());
4536 vector<Clustering> systems;
4540 for (
int i = 0; i < nFiGluon; ++i) {
4541 int EmtGluon = posFinalGluon[i];
4542 systems = findSQCDTriple( EmtGluon, 2, event, posFinalPartn, posInitPartn);
4543 ret.insert(ret.end(), systems.begin(), systems.end());
4549 bool check_g2qq =
true;
4550 if ( ( ( nInQuark + nInAntiq == 0 )
4552 && (nFiQuark == 1) && (nFiAntiq == 1) )
4553 || ( ( nFiQuark + nFiAntiq == 0)
4554 && (nInQuark == 1) && (nInAntiq == 1) ) )
4561 for(
int i=0; i < nFiQuark; ++i) {
4562 int EmtQuark = posFinalQuark[i];
4563 systems = findSQCDTriple( EmtQuark,1,event, posFinalPartn, posInitPartn);
4564 ret.insert(ret.end(), systems.begin(), systems.end());
4570 for(
int i=0; i < nFiAntiq; ++i) {
4571 int EmtAntiq = posFinalAntiq[i];
4572 systems = findSQCDTriple( EmtAntiq,1,event, posFinalPartn, posInitPartn);
4573 ret.insert(ret.end(), systems.begin(), systems.end());
4594 vector<Clustering> History::findSQCDTriple (
int EmtTagIn,
int colTopIn,
4596 vector<int> posFinalPartn,
4597 vector <int> posInitPartn ) {
4600 int EmtTag = EmtTagIn;
4603 int colTop = colTopIn;
4606 int offsetL = 1000000;
4607 int offsetR = 2000000;
4610 int finalSize = int(posFinalPartn.size());
4611 int initSize = int(posInitPartn.size());
4612 int size = initSize + finalSize;
4614 vector<Clustering> clus;
4618 for (
int a = 0; a < size; ++a ) {
4619 int i = (a < finalSize)? a : (a - finalSize) ;
4620 int iRad = (a < finalSize)? posFinalPartn[i] : posInitPartn[i];
4622 if ( event[iRad].col() ==
event[EmtTag].col()
4623 &&
event[iRad].acol() ==
event[EmtTag].acol() )
4627 int radID =
event[iRad].id();
4629 bool isSQCDrad = (abs(radID) > offsetL);
4631 bool isSQCDemt = (
event[EmtTag].idAbs() > offsetL );
4633 if (iRad != EmtTag ) {
4634 int pTdef =
event[iRad].isFinal() ? 1 : -1;
4635 int sign = (a < finalSize)? 1 : -1 ;
4638 int radBefID = getRadBeforeFlav(iRad,EmtTag,event);
4639 if ( pTdef == -1 && abs(radBefID) > offsetL )
continue;
4645 int radSign = (
event[iRad].id() < 0) ? -1 : 1;
4646 int emtSign = (
event[EmtTag].id() < 0) ? -1 : 1;
4649 bool finalSplitting =
false;
4650 if ( abs(radID) < 10
4651 && radSign*(abs(radID)+offsetL) == -sign*event[EmtTag].
id() )
4652 finalSplitting =
true;
4653 if ( abs(radID) < 10
4654 && radSign*(abs(radID)+offsetR) == -sign*event[EmtTag].
id() )
4655 finalSplitting =
true;
4656 if ( abs(radID) > offsetL && abs(radID) < offsetL+10
4657 && radID == -sign*emtSign*(
event[EmtTag].idAbs() + offsetL) )
4658 finalSplitting =
true;
4659 if ( abs(radID) > offsetR && abs(radID) < offsetR+10
4660 && radID == -sign*emtSign*(
event[EmtTag].idAbs() + offsetR) )
4661 finalSplitting =
true;
4664 bool initialSplitting =
false;
4665 if ( radID == 21 && ( ( event[EmtTag].idAbs() > offsetL
4666 && event[EmtTag].idAbs() < offsetL+10)
4667 || (
event[EmtTag].idAbs() > offsetR
4668 &&
event[EmtTag].idAbs() < offsetR+10) )
4669 && (
event[iRad].col() ==
event[EmtTag].col()
4670 ||
event[iRad].acol() ==
event[EmtTag].acol() ) )
4671 initialSplitting =
true;
4673 if ( finalSplitting ) {
4677 if ( radID < 0 && event[iRad].colType() == -1) {
4678 acl =
event[EmtTag].acol();
4679 col =
event[iRad].acol();
4680 }
else if ( event[iRad].colType() == 1 ) {
4681 col =
event[EmtTag].col();
4682 acl =
event[iRad].col();
4692 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
4695 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4699 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4700 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4708 if ( !isSQCDrad && !isSQCDemt
4709 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
4712 if ( iRec != 0 && iPartner != 0
4713 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4714 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4715 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4722 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
4725 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4729 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4730 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4738 if ( !isSQCDrad && !isSQCDemt
4739 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
4742 if ( iRec != 0 && iPartner != 0
4743 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4744 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4745 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4755 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
4758 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4762 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4763 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4771 if ( !isSQCDrad && !isSQCDemt
4772 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
4775 if ( iRec != 0 && iPartner != 0
4776 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4777 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4778 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4785 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
4788 if ( (sign < 0) && (event[iRec].isFinal()) ) {
4792 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4793 if (posInitPartn[l] != iRad) iRec = posInitPartn[l];
4801 if ( !isSQCDrad && !isSQCDemt
4802 && (event[iRec].idAbs() < 10 ||
event[iRec].id() == 21) )
4805 if ( iRec != 0 && iPartner != 0
4806 && allowedClustering( iRad, EmtTag, iRec, iPartner, event) ) {
4807 attachClusterings (clus, EmtTag, iRad, iRec, iPartner,
4808 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4813 }
else if ( initialSplitting ) {
4816 if ( !isSQCDrad && !isSQCDemt )
continue;
4823 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4824 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
4828 int col = getRadBeforeCol(iRad, EmtTag, event);
4829 int acl = getRadBeforeAcol(iRad, EmtTag, event);
4839 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
4840 ? event[iRad].col() : 0;
4843 if (colRemove > 0 && col > 0)
4844 iPartner = FindCol(col,iRad,EmtTag,event,1,
true)
4845 + FindCol(col,iRad,EmtTag,event,2,
true);
4846 else if (colRemove > 0 && acl > 0)
4847 iPartner = FindCol(acl,iRad,EmtTag,event,1,
true)
4848 + FindCol(acl,iRad,EmtTag,event,2,
true);
4850 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event ) ) {
4851 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
4852 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
4861 if ( (event[iRad].col() == event[EmtTag].acol())
4862 || (event[iRad].acol() == event[EmtTag].col())
4863 || (event[iRad].col() == event[EmtTag].col())
4864 || (event[iRad].acol() == event[EmtTag].acol()) ) {
4871 if (event[iRad].isFinal() ) {
4873 if ( radID < 0 && event[iRad].colType() == -1) {
4874 acl =
event[EmtTag].acol();
4875 col =
event[iRad].col();
4876 }
else if ( radID > 0 && event[iRad].colType() == 1 ) {
4877 col =
event[EmtTag].col();
4878 acl =
event[iRad].acol();
4880 col =
event[iRad].col();
4881 acl =
event[iRad].acol();
4886 iRec = FindCol(col,iRad,EmtTag,event,1,
true);
4887 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4889 if ( !isSQCDrad && !isSQCDemt
4890 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
4893 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4894 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4895 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4899 iRec = FindCol(col,iRad,EmtTag,event,2,
true);
4900 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4902 if ( !isSQCDrad && !isSQCDemt
4903 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
4906 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4907 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4908 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4914 iRec = FindCol(acl,iRad,EmtTag,event,1,
true);
4915 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4917 if ( !isSQCDrad && !isSQCDemt
4918 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
4921 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4922 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4923 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4927 iRec = FindCol(acl,iRad,EmtTag,event,2,
true);
4928 if ( (sign < 0) && (event[iRec].isFinal()) ) iRec = 0;
4930 if ( !isSQCDrad && !isSQCDemt
4931 && (event[iRec].idAbs() < 10 || event[iRec].
id() == 21) )
4934 && allowedClustering( iRad, EmtTag, iRec, iRec, event) ) {
4935 attachClusterings (clus, EmtTag, iRad, iRec, iRec,
4936 pTLund(event, iRad, EmtTag, iRec, pTdef), event);
4948 if ( !isSQCDrad || !isSQCDemt )
continue;
4956 for(
int l = 0; l < int(posInitPartn.size()); ++l)
4957 if (posInitPartn[l] != iRad) RecInit = posInitPartn[l];
4961 col = getRadBeforeCol(iRad, EmtTag, event);
4962 acl = getRadBeforeAcol(iRad, EmtTag, event);
4972 int colRemove = (
event[iRad].col() ==
event[EmtTag].col())
4973 ? event[iRad].col() : 0;
4974 iPartner = (colRemove > 0)
4975 ? FindCol(col,iRad,EmtTag,event,1,
true)
4976 + FindCol(col,iRad,EmtTag,event,2,
true)
4977 : FindCol(acl,iRad,EmtTag,event,1,true)
4978 + FindCol(acl,iRad,EmtTag,event,2,true);
4980 if ( allowedClustering( iRad, EmtTag, RecInit, iPartner, event)) {
4981 attachClusterings (clus, EmtTag, iRad, RecInit, iPartner,
4982 pTLund(event, iRad, EmtTag, RecInit, pTdef), event);
5003 double History::getProb(
const Clustering & SystemIn) {
5006 int Rad = SystemIn.emittor;
5007 int Rec = SystemIn.recoiler;
5008 int Emt = SystemIn.emitted;
5011 double showerProb = 0.0;
5015 if (SystemIn.pT() <= 0.)
return 0.;
5024 bool isFSR = (state[Rad].isFinal() && state[Rec].isFinal());
5025 bool isFSRinREC = (state[Rad].isFinal() && !state[Rec].isFinal());
5026 bool isISR = !state[Rad].isFinal();
5029 if ( mergingHooksPtr->useShowerPlugin() ) {
5030 int iPartner = (isISR && SystemIn.partner > 0) ? SystemIn.partner : Rec;
5033 bool isFSR2 = showers->timesPtr->isTimelike(state, Rad, Emt, iPartner,
"");
5035 string name = showers->timesPtr->getSplittingName( state, Rad, Emt,
5037 pr = showers->timesPtr->getSplittingProb( state, Rad, Emt,
5040 string name = showers->spacePtr->getSplittingName(state, Rad, Emt,
5042 pr = showers->spacePtr->getSplittingProb(state, Rad, Emt,
5051 for(
int i=0; i < state.size(); ++i)
5052 if (state[i].isFinal()) nFinal++;
5053 bool isLast = (nFinal == (mergingHooksPtr->hardProcess->nQuarksOut()
5054 +mergingHooksPtr->hardProcess->nLeptonOut()+1));
5057 bool isElectroWeak = (state[Emt].idAbs() == 23 || state[Emt].idAbs() == 24);
5064 for(
int i=0;i< int(state.size()); ++i) {
5065 if (state[i].mother1() == 1) inP = i;
5066 if (state[i].mother1() == 2) inM = i;
5069 Vec4 sum = state[Rad].p() + state[Rec].p() - state[Emt].p();
5070 double m2Dip = sum.m2Calc();
5071 double sHat = (state[inM].p() + state[inP].p()).m2Calc();
5073 double z1 = m2Dip / sHat;
5075 Vec4 Q1( state[Rad].p() - state[Emt].p() );
5076 Vec4 Q2( state[Rec].p() - state[Emt].p() );
5078 double Q1sq = -Q1.m2Calc();
5080 double pT1sq = pow2(pTLund(state, Rad, Emt, Rec, -1, SystemIn.flavRadBef));
5083 bool g2QQmassive = mergingHooksPtr->includeMassive()
5084 && state[Rad].id() == 21
5085 && ( (state[Emt].idAbs() >= 4 && state[Emt].idAbs() < 7)
5086 || (state[Emt].idAbs() > 1000000 && state[Emt].idAbs() < 1000010 )
5087 || (state[Emt].idAbs() > 2000000 && state[Emt].idAbs() < 2000010 ));
5088 bool Q2Qgmassive = mergingHooksPtr->includeMassive()
5089 && state[Emt].id() == 21
5090 && ( (state[Rad].idAbs() >= 4 && state[Rad].idAbs() < 7)
5091 || (state[Rad].idAbs() > 1000000 && state[Rad].idAbs() < 1000010 )
5092 || (state[Rad].idAbs() > 2000000 && state[Rad].idAbs() < 2000010 ));
5093 bool isMassive = mergingHooksPtr->includeMassive()
5094 && ( g2QQmassive || Q2Qgmassive
5095 || state[Emt].id() == 1000021);
5096 double m2Emt0 = state[Emt].p().m2Calc();
5097 double m2Rad0 = pow(particleDataPtr->m0(state[Rad].id()),2);
5100 if ( g2QQmassive) Q1sq += m2Emt0;
5101 else if (Q2Qgmassive) Q1sq += m2Rad0;
5104 double pT0sq = pow(mergingHooksPtr->pT0ISR(),2);
5105 double Q2sq = -Q2.m2Calc();
5108 bool g2QQmassiveRec = mergingHooksPtr->includeMassive()
5109 && state[Rec].id() == 21
5110 && ( state[Emt].idAbs() >= 4 && state[Emt].idAbs() < 7);
5111 bool Q2QgmassiveRec = mergingHooksPtr->includeMassive()
5112 && state[Emt].id() == 21
5113 && ( state[Rec].idAbs() >= 4 && state[Rec].idAbs() < 7);
5114 double m2Rec0 = pow(particleDataPtr->m0(state[Rec].id()),2);
5115 if ( g2QQmassiveRec) Q2sq += m2Emt0;
5116 else if (Q2QgmassiveRec) Q2sq += m2Rec0;
5118 bool hasJoinedEvol = (state[Emt].id() == 21
5119 || state[Rad].id() == state[Rec].id());
5124 if ( mergingHooksPtr->pickByFull() || mergingHooksPtr->pickBySumPT()) {
5125 double facJoined = ( Q2sq + pT0sq/(1.-z1) )
5126 * 1./(Q1sq*Q2sq + pT0sq*sHat + pow(pT0sq/(1.-z1),2));
5127 double facSingle = mergingHooksPtr->nonJoinedNorm()*1./( pT1sq + pT0sq);
5129 fac = (hasJoinedEvol && isLast) ? facJoined : facSingle;
5131 }
else if (mergingHooksPtr->pickByPoPT2()) {
5132 fac = 1./(pT1sq + pT0sq);
5134 string message=
"Error in History::getProb: Scheme for calculating";
5135 message+=
" shower splitting probability is undefined.";
5136 infoPtr->errorMsg(message);
5142 if ( isElectroWeak ) {
5151 }
else if ( state[Emt].
id() == 21 && state[Rad].
id() != 21) {
5153 double num = CF*(1. + pow(z1,2)) / (1.-z1);
5154 if (isMassive) num -= CF * z1 * (1.-z1) * (m2Rad0/pT1sq);
5157 double meReweighting = 1.;
5161 for(
int i=0; i < state.size(); ++i)
5162 if (state[i].isFinal() && state[i].colType() != 0
5163 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5168 &&
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1) {
5169 double sH = m2Dip / z1;
5171 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5172 double misMatch = (uH*tH - (uH + tH)*pT0sq/(1.-z1)
5173 + pow(pT0sq/(1.-z1),2) ) / (uH*tH);
5174 meReweighting *= (tH*tH + uH*uH + 2. * m2Dip * sH)
5175 / (sH*sH + m2Dip*m2Dip);
5176 meReweighting *= misMatch;
5179 showerProb = num*fac*meReweighting;
5182 }
else if ( state[Emt].
id() == 21 && state[Rad].id() == 21) {
5184 double num = 2.*CA*pow2(1. - z1*(1.-z1)) / (z1*(1.-z1));
5188 double meReweighting = 1.;
5192 for(
int i=0; i < state.size(); ++i)
5193 if (state[i].isFinal() && state[i].colType() != 0
5194 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5199 && mergingHooksPtr->getProcessString().compare(
"pp>h") == 0
5200 && int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1 ) {
5201 double sH = m2Dip / z1;
5203 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5204 meReweighting *= 0.5
5205 * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(m2Dip))
5206 / pow2(sH*sH - m2Dip * (sH - m2Dip));
5210 showerProb = num*fac*meReweighting;
5213 }
else if ( state[Emt].
id() != 21 && state[Rad].id() != 21 ) {
5215 double num = CF*(1. + pow2(1.-z1)) / z1;
5219 double meReweighting = 1.;
5223 for (
int i=0; i < state.size(); ++i )
5224 if ( state[i].isFinal() && state[i].colType() != 0
5225 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5230 && mergingHooksPtr->getProcessString().compare(
"pp>h") == 0
5231 && int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1) {
5232 double sH = m2Dip / z1;
5233 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5234 meReweighting *= (sH*sH + uH*uH)
5235 / (sH*sH + pow2(sH -m2Dip));
5239 showerProb = num*fac*meReweighting;
5242 }
else if ( state[Emt].
id() != 21 && state[Rad].id() == 21 ) {
5244 double num = TR * ( pow(z1,2) + pow(1.-z1,2) );
5245 if (isMassive) num += TR * 2.*z1*(1.-z1)*(m2Emt0/pT1sq);
5247 double meReweighting = 1.;
5251 for (
int i=0; i < state.size(); ++i )
5252 if ( state[i].isFinal() && state[i].colType() != 0
5253 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5258 &&
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1) {
5259 double sH = m2Dip / z1;
5261 double uH = Q1sq - m2Dip * (1. - z1) / z1;
5263 double misMatch = ( uH - pT0sq/(1.-z1) ) / uH;
5264 double me = (sH*sH + uH*uH + 2. * m2Dip * tH)
5265 / (pow2(sH - m2Dip) + m2Dip*m2Dip);
5267 meReweighting *= me;
5268 meReweighting *= misMatch;
5271 showerProb = num*fac*meReweighting;
5275 string message =
"Error in History::getProb: Splitting kernel"
5276 " undefined in ISR in clustering.";
5277 infoPtr->errorMsg(message);
5281 double m2Sister = pow(state[Emt].m(),2);
5282 double pT2corr = (Q1sq - z1*(m2Dip + Q1sq)*(Q1sq + m2Sister)/m2Dip);
5283 if (pT2corr < 0.) showerProb = 0.0;
5287 if ( state[Emt].
id() == state[Rad].id()
5288 && ( state[Rad].idAbs() == 4 || state[Rad].idAbs() == 5 )) {
5289 double m2QQsister = 2.*4.*m2Sister;
5290 double pT2QQcorr = Q1sq - z1*(m2Dip + Q1sq)*(Q1sq + m2QQsister)
5293 if (pT2QQcorr < 0.0) showerProb = 0.0;
5297 double pT2minNow = mergingHooksPtr->pTcut();
5298 double zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
5299 ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
5300 zMaxAbs = min(1.,zMaxAbs);
5301 double zMinAbs = max(0.,1. - zMaxAbs);
5304 int radBefID = getRadBeforeFlav(Rad, Emt, state);
5305 if ( abs(radBefID) == 4 || abs(radBefID) == 5 ) {
5306 double m2Massive = pow2(particleDataPtr->m0(radBefID));
5307 double mRatio = sqrt( m2Massive / m2Dip );
5308 double zMaxMassive = (1. - mRatio) / ( 1. + mRatio * (1. - mRatio) );
5309 zMaxAbs = min(zMaxAbs, zMaxMassive);
5312 if (z1 < zMinAbs || z1 > zMaxAbs) showerProb = 0.0;
5314 if (mergingHooksPtr->includeRedundant()) {
5316 AlphaStrong* asISR = mergingHooksPtr->AlphaS_ISR();
5317 double as = (*asISR).alphaS(pT1sq + pT0sq) / (2.*M_PI);
5324 }
else if (isFSR || isFSRinREC) {
5327 int recSign = (state[Rec].isFinal()) ? 1 : -1;
5328 Vec4 sum = state[Rad].p() + recSign*state[Rec].p() + state[Emt].p();
5329 double m2Dip = abs(sum.m2Calc());
5332 Vec4 Q1( state[Rad].p() + state[Emt].p() );
5333 Vec4 Q2( state[Rec].p() + state[Emt].p() );
5336 double z1 = getCurrentZ( Rad, Rec, Emt, SystemIn.flavRadBef);
5339 double Q1sq = Q1.m2Calc();
5341 double pT1sq = pow(SystemIn.pT(),2);
5343 double Q2sq = Q2.m2Calc();
5346 bool isMassiveRad = ( state[Rad].idAbs() >= 4
5347 && state[Rad].id() != 21 );
5348 bool isMassiveRec = ( state[Rec].idAbs() >= 4
5349 && state[Rec].id() != 21 );
5352 double m2Rad0 = pow(particleDataPtr->m0(state[Rad].id()),2);
5353 double m2Rec0 = pow(particleDataPtr->m0(state[Rec].id()),2);
5354 if ( mergingHooksPtr->includeMassive() && isMassiveRad ) Q1sq -= m2Rad0;
5355 if ( mergingHooksPtr->includeMassive() && isMassiveRec ) Q2sq -= m2Rec0;
5360 if ( mergingHooksPtr->pickByFull() || mergingHooksPtr->pickBySumPT()) {
5361 double facJoined = (1.-z1)/Q1sq * m2Dip/( Q1sq + Q2sq );
5362 double facSingle = mergingHooksPtr->fsrInRecNorm() * 1./ pT1sq;
5363 fac = (!isFSRinREC && isLast) ? facJoined : facSingle;
5365 }
else if (mergingHooksPtr->pickByPoPT2()) {
5368 string message=
"Error in History::getProb: Scheme for calculating";
5369 message+=
" shower splitting probability is undefined.";
5370 infoPtr->errorMsg(message);
5375 if ( isElectroWeak ) {
5384 }
else if ( state[Emt].
id() == 21 && state[Rad].colType() == 2) {
5387 double num = 0.5* CA * (1. + pow3(z1)) / (1.-z1);
5389 showerProb = num*fac;
5392 }
else if ( state[Emt].
id() == 21
5393 && state[Rad].colType() != 2 && state[Rec].colType() != 2) {
5396 double num = CF * 2./(1.-z1);
5400 for(
int i=0; i < state.size(); ++i)
5401 if (state[i].isFinal() && state[i].colType() != 0
5402 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,state) )
5406 ||
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) > 1)
5407 num = CF * (1. + pow2(z1)) /(1.-z1);
5409 double meReweighting = 1.;
5413 &&
int(mergingHooksPtr->hardProcess->hardIntermediate.size()) == 1 ) {
5415 double x1 = 2. * (sum * state[Rad].p()) / m2Dip;
5416 double x2 = 2. * (sum * state[Rec].p()) / m2Dip;
5417 double prop1 = max(1e-12, 1. - x1);
5418 double prop2 = max(1e-12, 1. - x2);
5419 double x3 = max(1e-12, 2. - x1 - x2);
5421 double ShowerRate1 = 2./( x3 * prop2 );
5422 double meDividingFactor1 = prop1 / x3;
5423 double me = (pow(x1,2) + pow(x2,2))/(prop1*prop2);
5424 meReweighting = meDividingFactor1 * me / ShowerRate1;
5427 showerProb = num*fac*meReweighting;
5430 }
else if ( state[Emt].
id() == 21 && state[Rad].colType() != 2
5431 && state[Rec].colType() == 2 ) {
5435 double num = CF * (1. + pow2(z1)) / (1.-z1);
5436 showerProb = num*fac;
5439 }
else if ( state[Emt].
id() != 21 ) {
5441 int flavour = state[Emt].id();
5444 double mFlavour = particleDataPtr->m0(flavour);
5446 double mDipole = m(state[Rad].p(), state[Emt].p());
5448 double beta = sqrtpos( 1. - 4.*pow2(mFlavour)/pow2(mDipole) );
5450 double num = 0.5*TR * ( z1*z1 + (1.-z1)*(1.-z1) );
5451 if (beta <= 0.) beta = 0.;
5453 showerProb = num*fac*beta;
5457 string message=
"Error in History::getProb: Splitting kernel undefined";
5458 message+=
" in FSR clustering.";
5459 infoPtr->errorMsg(message);
5462 if (mergingHooksPtr->includeRedundant()) {
5464 AlphaStrong* asFSR = mergingHooksPtr->AlphaS_FSR();
5465 double as = (*asFSR).alphaS(pT1sq) / (2.*M_PI);
5470 double m2DipCorr = pow2(sqrt(m2Dip) - sqrt(m2Rec0)) - m2Rad0;
5471 double zMin = 0.5 - sqrtpos( 0.25 - pT1sq / m2DipCorr );
5472 double m2 = m2Rad0 + 2. * state[Rad].p()*state[Emt].p();
5473 bool keepMassive = (z1 > zMin && z1 < 1. - zMin
5474 && m2 * m2Dip < z1 * (1. - z1) * pow2(m2Dip + m2 - m2Rec0) );
5476 if (!keepMassive) showerProb *= 0.0;
5480 string message=
"Error in History::getProb: Radiation could not be";
5481 message+=
" interpreted as FSR or ISR.";
5482 infoPtr->errorMsg(message);
5485 if (showerProb <= 0.) showerProb = 0.;
5487 if (mergingHooksPtr->doWeakClustering()) {
5494 if (state[Rad].idAbs() < 10 && state[Rad].intPol() == 9
5495 && SystemIn.spinRad != 9) factor *= 0.5;
5496 if (state[Emt].idAbs() < 10 && state[Emt].intPol() == 9
5497 && SystemIn.spinEmt != 9) factor *= 0.5;
5498 if (state[Rec].idAbs() < 10 && state[Rec].intPol() == 9
5499 && SystemIn.spinRec != 9) factor *= 0.5;
5500 if ( state[Emt].colType() != 0 ) {
5501 showerProb *= factor;
5507 if ( state[Emt].idAbs() < 10 && state[Rad].colType() != 0) {
5530 void History::setupBeams() {
5535 if (state.size() < 4)
return;
5537 if ( state[3].colType() == 0 )
return;
5538 if ( state[4].colType() == 0 )
return;
5544 for(
int i=0;i< int(state.size()); ++i) {
5545 if (state[i].mother1() == 1) inP = i;
5546 if (state[i].mother1() == 2) inM = i;
5551 int motherPcompRes = -1;
5552 int motherMcompRes = -1;
5554 bool sameFlavP =
false;
5555 bool sameFlavM =
false;
5560 for(
int i=0;i< int(mother->state.size()); ++i) {
5561 if (mother->state[i].mother1() == 1) inMotherP = i;
5562 if (mother->state[i].mother1() == 2) inMotherM = i;
5564 sameFlavP = (state[inP].id() == mother->state[inMotherP].id());
5565 sameFlavM = (state[inM].id() == mother->state[inMotherM].id());
5567 motherPcompRes = (sameFlavP) ? beamA[0].companion() : -2;
5568 motherMcompRes = (sameFlavM) ? beamB[0].companion() : -2;
5576 double Ep = 2. * state[inP].e();
5577 double Em = 2. * state[inM].e();
5580 if (state[inP].m() != 0. || state[inM].m() != 0.) {
5581 Ep = state[inP].pPos() + state[inM].pPos();
5582 Em = state[inP].pNeg() + state[inM].pNeg();
5586 double x1 = Ep / state[inS].m();
5587 beamA.append( inP, state[inP].
id(), x1);
5588 double x2 = Em / state[inS].m();
5589 beamB.append( inM, state[inM].
id(), x2);
5593 double scalePDF = (mother) ? scale : infoPtr->QFac();
5597 beamA.xfISR( 0, state[inP].
id(), x1, scalePDF*scalePDF);
5599 beamA.pickValSeaComp();
5601 beamA[0].companion(motherPcompRes);
5603 beamB.xfISR( 0, state[inM].
id(), x2, scalePDF*scalePDF);
5605 beamB.pickValSeaComp();
5607 beamB[0].companion(motherMcompRes);
5617 double History::pdfForSudakov() {
5620 if ( state[3].colType() == 0 )
return 1.0;
5621 if ( state[4].colType() == 0 )
return 1.0;
5624 bool FSR = ( mother->state[clusterIn.emittor].isFinal()
5625 && mother->state[clusterIn.recoiler].isFinal());
5626 bool FSRinRec = ( mother->state[clusterIn.emittor].isFinal()
5627 && !mother->state[clusterIn.recoiler].isFinal());
5630 if (FSR)
return 1.0;
5632 int iInMother = (FSRinRec)? clusterIn.recoiler : clusterIn.emittor;
5634 int side = ( mother->state[iInMother].pz() > 0 ) ? 1 : -1;
5638 for(
int i=0;i< int(state.size()); ++i) {
5639 if (state[i].mother1() == 1) inP = i;
5640 if (state[i].mother1() == 2) inM = i;
5644 int idMother = mother->state[iInMother].id();
5646 int iDau = (side == 1) ? inP : inM;
5647 int idDaughter = state[iDau].id();
5649 double xMother = 2. * mother->state[iInMother].e() / mother->state[0].e();
5651 double xDaughter = 2.*state[iDau].e() / state[0].e();
5654 double ratio = getPDFratio(side,
true,
false, idMother, xMother, scale,
5655 idDaughter, xDaughter, scale);
5660 return ( (FSRinRec)? min(1.,ratio) : ratio);
5668 double History::hardProcessME(
const Event& event ) {
5671 if (isEW2to1(event)) {
5674 if (event[5].idAbs() == 24) {
5675 int idIn1 =
event[3].id();
5676 int idIn2 =
event[4].id();
5677 double mW = particleDataPtr->m0(24);
5678 double gW = particleDataPtr->mWidth(24) / mW;
5679 double sH = (
event[3].p()+
event[4].p()).m2Calc();
5681 double thetaWRat = 1. / (12. * coupSMPtr->sin2thetaW());
5682 double ckmW = coupSMPtr->V2CKMid(abs(idIn1), abs(idIn2));
5684 double bwW = 12. * M_PI / ( pow2(sH - pow2(mW)) + pow2(sH * gW) );
5685 double preFac = thetaWRat * sqrt(sH) * particleDataPtr->mWidth(24);
5686 return ckmW * preFac * bwW;
5689 else if (event[5].idAbs() == 23) {
5690 double mZ = particleDataPtr->m0(23);
5691 double gZ = particleDataPtr->mWidth(23) / mZ;
5692 double sH = (
event[3].p()+
event[4].p()).m2Calc();
5694 double thetaZRat = (pow2(coupSMPtr->rf( abs(clusterIn.flavRadBef))) +
5695 pow2(coupSMPtr->lf( abs(clusterIn.flavRadBef)))) /
5696 (24. * coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW());
5697 double bwW = 12. * M_PI / ( pow2(sH - pow2(mZ)) + pow2(sH * gZ) );
5698 double preFac = thetaZRat * sqrt(sH) * particleDataPtr->mWidth(23);
5699 return preFac * bwW;
5702 string message=
"Warning in History::hardProcessME: Only Z/W are";
5703 message+=
" supported as 2->1 processes. Skipping history.";
5704 infoPtr->errorMsg(message);
5709 else if (isQCD2to2(event)) {
5710 int idIn1 =
event[3].id();
5711 int idIn2 =
event[4].id();
5712 int idOut1 =
event[5].id();
5713 int idOut2 =
event[6].id();
5715 double sH = (
event[3].p()+
event[4].p()).m2Calc();
5716 double tH = (
event[3].p()-
event[5].p()).m2Calc();
5717 double uH = (
event[3].p()-
event[6].p()).m2Calc();
5721 if (!(abs(idIn1) < 10 || abs(idIn1) == 21) ) isQCD =
false;
5722 if (!(abs(idIn2) < 10 || abs(idIn2) == 21) ) isQCD =
false;
5723 if (!(abs(idOut1) < 10 || abs(idOut1) == 21) ) isQCD =
false;
5724 if (!(abs(idOut2) < 10 || abs(idOut2) == 21) ) isQCD =
false;
5727 double cor = M_PI / (9. * pow2(sH));
5734 if (abs(idIn1) == 21 && abs(idIn2) == 21) {
5735 if (abs(idOut1) == 21 && abs(idOut2) == 21)
5736 return cor * weakShowerMEs.getMEgg2gg(sH, tH, uH);
5737 else return cor * weakShowerMEs.getMEgg2qqbar(sH, tH, uH);
5740 }
else if (abs(idIn1) == 21 || abs(idIn2) == 21) {
5741 if (idIn1 != idOut1) swap(uH, tH);
5742 return cor * weakShowerMEs.getMEqg2qg(sH, tH, uH);
5747 if (abs(idOut1) == 21 && abs(idOut2) == 21)
5748 return cor * weakShowerMEs.getMEqqbar2gg(sH, tH, uH);
5749 if (idIn1 == -idIn2) {
5750 if (abs(idIn1) == abs(idOut1)) {
5751 if (idIn1 != idOut1) swap(uH, tH);
5752 return cor * weakShowerMEs.getMEqqbar2qqbar(sH, tH, uH,
true);
5755 return cor * weakShowerMEs.getMEqqbar2qqbar(sH, tH, uH,
false);
5758 else if (idIn1 == idIn2)
5759 return cor * weakShowerMEs.getMEqq2qq(sH, tH, uH,
true);
5761 if (idIn1 == idOut1) swap(uH,tH);
5762 return cor * weakShowerMEs.getMEqq2qq(sH, tH, uH,
false);
5770 string process = mergingHooksPtr->getProcessString();
5773 if ( process.compare(
"pp>e+ve") == 0
5774 || process.compare(
"pp>e-ve~") == 0
5775 || process.compare(
"pp>LEPTONS,NEUTRINOS") == 0 ) {
5778 for (
int i=0; i < int(event.size()); ++i )
5779 if ( event[i].isFinal() ) nFinal++;
5780 if ( nFinal != 2 )
return 1.;
5782 double mW = particleDataPtr->m0(24);
5783 double gW = particleDataPtr->mWidth(24) / mW;
5785 int inP = (
event[3].pz() > 0) ? 3 : 4;
5786 int inM = (
event[3].pz() > 0) ? 4 : 3;
5789 for (
int i=0; i < int(event.size()); ++i ) {
5790 if ( event[i].isFinal() &&
event[i].px() > 0 ) outP = i;
5793 double sH = (
event[inP].p() +
event[inM].p()).m2Calc();
5794 double tH = (
event[inP].p() -
event[outP].p()).m2Calc();
5795 double uH = - sH - tH;
5798 result = ( 1. + (tH - uH)/sH ) / ( pow2(sH - mW*mW) + pow2(sH*gW) );
5800 result = mergingHooksPtr->hardProcessME(event);
5813 Event History::cluster( Clustering & inSystem ) {
5816 int Rad = inSystem.emittor;
5817 int Rec = inSystem.recoiler;
5818 int Emt = inSystem.emitted;
5820 double eCM = state[0].e();
5822 int radType = state[Rad].isFinal() ? 1 : -1;
5823 int recType = state[Rec].isFinal() ? 1 : -1;
5827 NewEvent.init(
"(hard process-modified)", particleDataPtr);
5831 if ( mergingHooksPtr->useShowerPlugin() ) {
5832 int iPartner = (radType == -1 && inSystem.partner > 0)
5833 ? inSystem.partner : Rec;
5834 bool isFSR = showers->timesPtr->isTimelike(state, Rad, Emt, iPartner,
"");
5836 string name = showers->timesPtr->getSplittingName(state,Rad,Emt,
5838 return showers->timesPtr->clustered(state,Rad,Emt,iPartner,name);
5840 string name = showers->spacePtr->getSplittingName(state,Rad,Emt,
5842 return showers->spacePtr->clustered(state,Rad,Emt,iPartner,name);
5847 for (
int i = 0; i < state.size(); ++i)
5848 if ( i != Rad && i != Rec && i != Emt )
5849 NewEvent.append( state[i] );
5852 for (
int i = 0; i < state.sizeJunction(); ++i)
5853 NewEvent.appendJunction( state.getJunction(i) );
5855 double mu = choseHardScale(state);
5857 NewEvent.saveSize();
5858 NewEvent.saveJunctionSize();
5860 NewEvent.scaleSecond(mu);
5864 Particle RecBefore = Particle( state[Rec] );
5865 RecBefore.setEvtPtr(&NewEvent);
5866 RecBefore.daughters(0,0);
5867 RecBefore.pol(inSystem.spinRec);
5869 int radID = inSystem.flavRadBef;
5870 if (radID == 0) radID = getRadBeforeFlav(Rad, Emt, state);
5871 int recID = state[Rec].id();
5872 Particle RadBefore = Particle( state[Rad] );
5873 RadBefore.setEvtPtr(&NewEvent);
5874 RadBefore.id(radID);
5875 RadBefore.pol(inSystem.spinRadBef);
5876 RadBefore.daughters(0,0);
5878 RadBefore.cols(RecBefore.acol(),RecBefore.col());
5881 if ( particleDataPtr->isResonance(radID) && radType == 1)
5882 RadBefore.status(state[Rad].status());
5885 double radMass = particleDataPtr->m0(radID);
5886 double recMass = particleDataPtr->m0(recID);
5887 if (radType == 1 ) RadBefore.m(radMass);
5888 else RadBefore.m(0.0);
5889 if (recType == 1 ) RecBefore.m(recMass);
5890 else RecBefore.m(0.0);
5894 if ( radType + recType == 2 && state[Emt].idAbs() != 23 &&
5895 state[Emt].idAbs() != 24) {
5898 Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p();
5899 double eCMME = sum.mCalc();
5905 double mDsq = pow2(eCMME);
5907 double mRsq = (radID == state[Rad].id() )
5908 ? abs(state[Rad].p().m2Calc())
5909 : pow2(particleDataPtr->m0(radID));
5910 double mSsq = abs(state[Rec].p().m2Calc());
5911 double a = 0.5*sqrt(mDsq);
5912 double b = 0.25*(mRsq - mSsq) / a;
5913 double c = sqrt(pow2(a) + pow2(b) - 2.*a*b - mSsq);
5915 Rad4mom.p( 0., 0., c, a+b);
5916 Rec4mom.p( 0., 0.,-c, a-b);
5919 Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p());
5920 Vec4 old2 = Vec4(state[Rec].p());
5921 RotBstMatrix fromCM;
5922 fromCM.fromCMframe(old1, old2);
5924 Rad4mom.rotbst(fromCM);
5925 Rec4mom.rotbst(fromCM);
5927 RadBefore.p(Rad4mom);
5928 RecBefore.p(Rec4mom);
5929 RadBefore.m(sqrt(mRsq));
5930 RecBefore.m(sqrt(mSsq));
5932 }
else if ( radType + recType == 2 && (state[Emt].idAbs() == 23 ||
5933 state[Emt].idAbs() == 24) ) {
5936 Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p();
5937 double eCMME = sum.mCalc();
5943 double mDsq = pow2(eCMME);
5945 double mRsq = (radID == state[Rad].id() )
5946 ? abs(state[Rad].p().m2Calc())
5947 : pow2(particleDataPtr->m0(radID));
5948 double mSsq = abs(state[Rec].p().m2Calc());
5949 double a = 0.5*sqrt(mDsq);
5950 double b = 0.25*(mRsq - mSsq) / a;
5951 double c = sqrt(pow2(a) + pow2(b) - 2.*a*b - mSsq);
5953 Rad4mom.p( 0., 0., c, a+b);
5954 Rec4mom.p( 0., 0.,-c, a-b);
5957 Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p());
5958 Vec4 old2 = Vec4(state[Rec].p());
5959 RotBstMatrix fromCM;
5960 fromCM.fromCMframe(old1, old2);
5962 Rad4mom.rotbst(fromCM);
5963 Rec4mom.rotbst(fromCM);
5965 RadBefore.p(Rad4mom);
5966 RecBefore.p(Rec4mom);
5967 RadBefore.m(sqrt(mRsq));
5968 RecBefore.m(sqrt(mSsq));
5970 }
else if ( radType + recType == 0 ) {
5974 Vec4 sum = state[Rad].p() + state[Rec].p() + state[Emt].p();
5975 double eCMME = sum.mCalc();
5980 double mDsq = pow2(eCMME);
5982 double mRsq = (radID == state[Rad].id() )
5983 ? abs(state[Rad].p().m2Calc())
5984 : pow2(particleDataPtr->m0(radID));
5985 double mSsq = abs(state[Rec].p().m2Calc());
5986 double a = 0.5*sqrt(mDsq);
5987 double b = 0.25*(mRsq - mSsq) / a;
5988 double c = sqrt(pow2(a) + pow2(b) - 2.*a*b - mSsq);
5990 Rad4mom.p( 0., 0., c, a+b);
5991 Rec4mom.p( 0., 0.,-c, a-b);
5994 Vec4 old1 = Vec4(state[Rad].p() + state[Emt].p());
5995 Vec4 old2 = Vec4(state[Rec].p());
5996 RotBstMatrix fromCM;
5997 fromCM.fromCMframe(old1, old2);
5999 Rad4mom.rotbst(fromCM);
6000 Rec4mom.rotbst(fromCM);
6003 Rec4mom = 2.*state[Rec].p() - Rec4mom;
6007 if ( abs(Rec4mom.mCalc()) > 1e-7 ) {
6008 double pzSign = (Rec4mom.pz() > 0.) ? 1. : -1.;
6009 double eRec = Rec4mom.e();
6010 Rec4mom.p(0., 0., pzSign*eRec, eRec);
6013 RadBefore.p(Rad4mom);
6014 RecBefore.p(Rec4mom);
6015 RadBefore.m(sqrt(mRsq));
6028 Vec4 pMother( state[Rad].p() );
6029 Vec4 pSister( state[Emt].p() );
6030 Vec4 pPartner( state[Rec].p() );
6031 Vec4 pDaughterBef( 0.,0.,0.,0. );
6032 Vec4 pRecoilerBef( 0.,0.,0.,0. );
6033 Vec4 pDaughter( 0.,0.,0.,0. );
6034 Vec4 pRecoiler( 0.,0.,0.,0. );
6037 int sign = (state[Rad].pz() > 0.) ? 1 : -1;
6041 double phi = pSister.phi();
6043 RotBstMatrix rot_by_mphi;
6044 rot_by_mphi.rot(0.,-phi);
6046 RotBstMatrix rot_by_pphi;
6047 rot_by_pphi.rot(0.,phi);
6051 double x1 = 2. * pMother.e() / eCM;
6053 double x2 = 2. * pPartner.e() / eCM;
6056 Vec4 qDip( pMother - pSister);
6057 Vec4 qAfter(pMother + pPartner);
6058 Vec4 qBefore(qDip + pPartner);
6059 double z = qBefore.m2Calc() / qAfter.m2Calc();
6062 double x1New = z*x1;
6064 double sHat = x1New*x2New*eCM*eCM;
6069 pDaughterBef.p( 0., 0., sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
6070 pRecoilerBef.p( 0., 0., -sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
6073 pMother.rotbst( rot_by_mphi );
6074 pSister.rotbst( rot_by_mphi );
6075 pPartner.rotbst( rot_by_mphi );
6076 for(
int i=3; i< NewEvent.size(); ++i)
6077 NewEvent[i].rotbst( rot_by_mphi );
6081 pDaughter.p( pMother - pSister);
6082 pRecoiler.p( pPartner );
6083 RotBstMatrix from_CM_to_DRoff;
6085 from_CM_to_DRoff.toCMframe(pDaughter, pRecoiler);
6087 from_CM_to_DRoff.toCMframe(pRecoiler, pDaughter);
6091 pMother.rotbst( from_CM_to_DRoff );
6092 pPartner.rotbst( from_CM_to_DRoff );
6093 pSister.rotbst( from_CM_to_DRoff );
6094 for(
int i=3; i< NewEvent.size(); ++i)
6095 NewEvent[i].rotbst( from_CM_to_DRoff );
6100 RotBstMatrix from_DR_to_CM;
6101 from_DR_to_CM.bst( 0., 0., sign*( x1New - x2New ) / ( x1New + x2New ) );
6106 pDaughterBef.rotbst( from_DR_to_CM );
6107 pRecoilerBef.rotbst( from_DR_to_CM );
6108 for(
int i=3; i< NewEvent.size(); ++i)
6109 NewEvent[i].rotbst( from_DR_to_CM );
6112 for(
int i=3; i< NewEvent.size(); ++i)
6113 NewEvent[i].rotbst( rot_by_pphi );
6117 if ( abs(pRecoilerBef.mCalc()) > 1e-7 ) {
6118 double pzSign = (pRecoilerBef.pz() > 0.) ? 1. : -1.;
6119 double eRec = pRecoilerBef.e();
6120 pRecoilerBef.p(0., 0., pzSign*eRec, eRec);
6122 if ( abs(pDaughterBef.mCalc()) > 1e-7 ) {
6123 double pzSign = (pDaughterBef.pz() > 0.) ? 1. : -1.;
6124 double eDau = pDaughterBef.e();
6125 pDaughterBef.p(0., 0., pzSign*eDau, eDau);
6129 RecBefore.p( pRecoilerBef );
6130 RadBefore.p( pDaughterBef );
6131 if (RecBefore.pz() > 0.) RecBefore.mother1(1);
6132 else RecBefore.mother1(2);
6133 if (RadBefore.pz() > 0.) RadBefore.mother1(1);
6134 else RadBefore.mother1(2);
6139 RecBefore.scale(mu);
6140 RadBefore.scale(mu);
6143 NewEvent.append(RecBefore);
6147 int emtID = state[Emt].id();
6148 if ( emtID == 22 || emtID == 23 || abs(emtID) == 24 )
6149 RadBefore.cols( state[Rad].col(), state[Rad].acol() );
6151 else if ( !connectRadiator( RadBefore, radType, RecBefore, recType,
6161 outState.init(
"(hard process-modified)", particleDataPtr);
6165 for (
int i = 0; i < 3; ++i)
6166 outState.append( NewEvent[i] );
6168 for (
int i = 0; i < state.sizeJunction(); ++i)
6169 outState.appendJunction( state.getJunction(i) );
6171 outState.saveSize();
6172 outState.saveJunctionSize();
6174 outState.scaleSecond(mu);
6175 bool radAppended =
false;
6176 bool recAppended =
false;
6177 int size = int(outState.size());
6179 int radPos = 0, recPos = 0;
6182 if ( RecBefore.mother1() == 1) {
6183 recPos = outState.append( RecBefore );
6185 }
else if ( RadBefore.mother1() == 1 ) {
6186 radPos = outState.append( RadBefore );
6191 for(
int i=0; i < int(state.size()); ++i)
6192 if (state[i].mother1() == 1) in1 =i;
6193 outState.append( state[in1] );
6197 if ( RecBefore.mother1() == 2) {
6198 recPos = outState.append( RecBefore );
6200 }
else if ( RadBefore.mother1() == 2 ) {
6201 radPos = outState.append( RadBefore );
6206 for(
int i=0; i < int(state.size()); ++i)
6207 if (state[i].mother1() == 2) in2 =i;
6209 outState.append( state[in2] );
6214 if (!recAppended && !RecBefore.isFinal()) {
6216 recPos = outState.append( RecBefore);
6219 if (!radAppended && !RadBefore.isFinal()) {
6221 radPos = outState.append( RadBefore);
6228 for (
int i = 0; i < int(NewEvent.size()-1); ++i)
6229 if (NewEvent[i].status() == -22) outState.append( NewEvent[i] );
6231 for (
int i = 0; i < int(NewEvent.size()-1); ++i)
6232 if (NewEvent[i].status() == 22) outState.append( NewEvent[i] );
6234 if (!radAppended && RadBefore.statusAbs() == 22)
6235 radPos = outState.append(RadBefore);
6237 recPos = outState.append(RecBefore);
6238 if (!radAppended && RadBefore.statusAbs() != 22)
6239 radPos = outState.append(RadBefore);
6241 for(
int i = 0; i < int(NewEvent.size()-1); ++i)
6242 if ( NewEvent[i].status() != 22
6243 && NewEvent[i].colType() != 0
6244 && NewEvent[i].isFinal())
6245 outState.append( NewEvent[i] );
6247 for(
int i = 0; i < int(NewEvent.size()-1); ++i)
6248 if ( NewEvent[i].status() != 22
6249 && NewEvent[i].colType() == 0
6250 && NewEvent[i].isFinal() )
6251 outState.append( NewEvent[i]);
6254 vector<int> posIntermediate;
6255 vector<int> posDaughter1;
6256 vector<int> posDaughter2;
6257 for(
int i=0; i < int(outState.size()); ++i)
6258 if (outState[i].status() == -22) {
6259 posIntermediate.push_back(i);
6260 int d1 = outState[i].daughter1();
6261 int d2 = outState[i].daughter2();
6263 int daughter1 = FindParticle( state[d1], outState);
6264 int daughter2 = FindParticle( state[d2], outState);
6269 posDaughter1.push_back( daughter1);
6272 while(!outState[daughter1].isFinal() ) daughter1++;
6273 posDaughter1.push_back( daughter1);
6276 posDaughter2.push_back( daughter2);
6278 daughter2 = outState.size()-1;
6279 while(!outState[daughter2].isFinal() ) daughter2--;
6280 posDaughter2.push_back( daughter2);
6284 for(
int i=0; i < int(posIntermediate.size()); ++i) {
6285 outState[posIntermediate[i]].daughters(posDaughter1[i],posDaughter2[i]);
6286 outState[posDaughter1[i]].mother1(posIntermediate[i]);
6287 outState[posDaughter2[i]].mother1(posIntermediate[i]);
6291 int minParFinal = int(outState.size());
6292 int maxParFinal = 0;
6293 for(
int i=0; i < int(outState.size()); ++i)
6294 if (outState[i].mother1() == 3 && outState[i].mother2() == 4) {
6295 minParFinal = min(i,minParFinal);
6296 maxParFinal = max(i,maxParFinal);
6299 if (minParFinal == maxParFinal) maxParFinal = 0;
6300 outState[3].daughters(minParFinal,maxParFinal);
6301 outState[4].daughters(minParFinal,maxParFinal);
6304 outState.saveSize();
6305 outState.saveJunctionSize();
6308 inSystem.recBef = recPos;
6309 inSystem.radBef = radPos;
6323 if ( radType == -1 && state[Rad].colType() == 1) {
6325 for(
int i=0; i < int(state.size()); ++i)
6326 if ( i != Rad && i != Emt && i != Rec
6327 && state[i].status() == -22
6328 && state[i].col() == state[Rad].col() )
6330 }
else if ( radType == -1 && state[Rad].colType() == -1) {
6332 for(
int i=0; i < int(state.size()); ++i)
6333 if ( i != Rad && i != Emt && i != Rec
6334 && state[i].status() == -22
6335 && state[i].acol() == state[Rad].acol() )
6337 }
else if ( radType == 1 && state[Rad].colType() == 1) {
6339 for(
int i=0; i < int(state.size()); ++i)
6340 if ( i != Rad && i != Emt && i != Rec
6341 && state[i].status() == -22
6342 && state[i].col() == state[Rad].col() )
6344 }
else if ( radType == 1 && state[Rad].colType() == -1) {
6346 for(
int i=0; i < int(state.size()); ++i)
6347 if ( i != Rad && i != Emt && i != Rec
6348 && state[i].status() == -22
6349 && state[i].acol() == state[Rad].acol() )
6355 int iColResNow = FindParticle( state[iColRes], outState);
6358 int radCol = outState[radPos].col();
6359 int radAcl = outState[radPos].acol();
6361 int resCol = outState[iColResNow].col();
6362 int resAcl = outState[iColResNow].acol();
6364 bool matchesRes = (radCol > 0
6365 && ( radCol == resCol || radCol == resAcl))
6367 && ( radAcl == resCol || radAcl == resAcl));
6371 if (!matchesRes && iColResNow > 0) {
6372 if ( radType == -1 && outState[radPos].colType() == 1)
6373 outState[iColResNow].col(radCol);
6374 else if ( radType ==-1 && outState[radPos].colType() ==-1)
6375 outState[iColResNow].acol(radAcl);
6376 else if ( radType == 1 && outState[radPos].colType() == 1)
6377 outState[iColResNow].col(radCol);
6378 else if ( radType == 1 && outState[radPos].colType() ==-1)
6379 outState[iColResNow].acol(radAcl);
6386 if (!matchesRes && iColResNow > 0 && iColRes != iColResNow)
6387 outState[radPos].mother1(iColResNow);
6392 if ( !validEvent(outState) ) {
6399 iReclusteredNew = radPos;
6413 int History::getRadBeforeFlav(
const int RadAfter,
const int EmtAfter,
6414 const Event& event) {
6416 int type =
event[RadAfter].isFinal() ? 1 :-1;
6417 int emtID =
event[EmtAfter].id();
6418 int radID =
event[RadAfter].id();
6419 int emtCOL =
event[EmtAfter].col();
6420 int radCOL =
event[RadAfter].col();
6421 int emtACL =
event[EmtAfter].acol();
6422 int radACL =
event[RadAfter].acol();
6424 bool colConnected = ((type == 1) && ( (emtCOL !=0 && (emtCOL ==radACL))
6425 || (emtACL !=0 && (emtACL ==radCOL)) ))
6426 ||((type ==-1) && ( (emtCOL !=0 && (emtCOL ==radCOL))
6427 || (emtACL !=0 && (emtACL ==radACL)) ));
6433 if ( type == 1 && emtID == -radID && !colConnected )
6436 if ( type ==-1 && radID == 21 )
6439 if ( type ==-1 && !colConnected
6440 && emtID != 21 && radID != 21 && abs(emtID) < 10 && abs(radID) < 10)
6444 int radSign = (radID < 0) ? -1 : 1;
6445 int offsetL = 1000000;
6446 int offsetR = 2000000;
6448 if ( emtID == 1000021 ) {
6450 if (abs(radID) < 10 ) {
6451 int offset = offsetL;
6454 for (
int i=0; i < int(event.size()); ++i)
6455 if ( event[i].isFinal()
6456 &&
event[i].idAbs() < offsetR+10 &&
event[i].idAbs() > offsetR)
6458 return radSign*(abs(radID)+offset);
6461 if (abs(radID) > offsetL && abs(radID) < offsetL+10 )
6462 return radSign*(abs(radID)-offsetL);
6463 if (abs(radID) > offsetR && abs(radID) < offsetR+10 )
6464 return radSign*(abs(radID)-offsetR);
6466 if (radID == 21 )
return emtID;
6469 int emtSign = (emtID < 0) ? -1 : 1;
6472 if ( abs(emtID) > offsetL && abs(emtID) < offsetL+10 )
6473 emtOffset = offsetL;
6474 if ( abs(emtID) > offsetR && abs(emtID) < offsetR+10 )
6475 emtOffset = offsetR;
6477 if ( abs(radID) > offsetL && abs(radID) < offsetL+10 )
6478 radOffset = offsetL;
6479 if ( abs(radID) > offsetR && abs(radID) < offsetR+10 )
6480 radOffset = offsetR;
6483 if ( type == 1 && !colConnected ) {
6485 if ( emtOffset > 0 && radOffset == 0
6486 && emtSign*(abs(emtID) - emtOffset) == -radID )
6489 if ( emtOffset == 0 && radOffset > 0
6490 && emtID == -radSign*(abs(radID) - radOffset) )
6495 if ( type ==-1 && radID == 1000021 ) {
6497 if ( emtOffset > 0 )
return -emtSign*(abs(emtID) - emtOffset);
6499 else return -emtSign*(abs(emtID) + emtOffset);
6504 && ( (abs(emtID) > offsetL && abs(emtID) < offsetL+10)
6505 || (abs(emtID) > offsetR && abs(emtID) < offsetR+10))
6506 && ( (abs(radID) > offsetL && abs(radID) < offsetL+10)
6507 || (abs(radID) > offsetR && abs(radID) < offsetR+10))
6508 && emtSign*(abs(emtID)+emtOffset) == radSign*(abs(radID) - radOffset)
6509 && !colConnected ) {
6515 double m2final = (
event[RadAfter].p()+
event[EmtAfter].p()).m2Calc();
6517 if ( emtID == 22 || emtID == 23 )
return radID;
6519 if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) <= 10. )
6522 if ( type == 1 && emtID == -radID && colConnected && sqrt(m2final) > 10. )
6525 if ( type ==-1 && (radID == 22 || radID == 23) )
6528 if ( type ==-1 && abs(emtID) < 10 && abs(radID) < 10 && colConnected )
6533 if ( emtID == 24 && radID < 0 )
return radID + 1;
6534 if ( emtID == 24 && radID > 0 )
return radID + 1;
6538 if ( emtID ==-24 && radID < 0 )
return radID - 1;
6539 if ( emtID ==-24 && radID > 0 )
return radID - 1;
6553 int History::getRadBeforeSpin(
const int radAfter,
const int emtAfter,
6554 const int spinRadAfter,
const int spinEmtAfter,
6555 const Event& event) {
6558 int radBeforeFlav = getRadBeforeFlav(radAfter, emtAfter, event);
6561 if ( event[radAfter].isFinal()
6562 && event[radAfter].
id() == -event[emtAfter].
id())
6563 return (spinRadAfter == 9) ? spinEmtAfter : spinRadAfter;
6566 if ( event[radAfter].isFinal() && abs(radBeforeFlav) < 10
6567 && event[radAfter].idAbs() < 10)
6569 return spinRadAfter;
6572 if ( event[radAfter].isFinal() && abs(radBeforeFlav) < 10
6573 && event[emtAfter].idAbs() < 10)
6575 return spinEmtAfter;
6578 if ( event[radAfter].isFinal() && radBeforeFlav == 21
6579 && event[radAfter].
id() == 21)
6581 return (spinRadAfter == 9) ? spinEmtAfter : spinRadAfter;
6584 if ( !event[radAfter].isFinal()
6585 && radBeforeFlav == -event[emtAfter].
id())
6586 return (spinRadAfter == 9) ? spinEmtAfter : spinRadAfter;
6589 if ( !event[radAfter].isFinal() && abs(radBeforeFlav) < 10
6590 && event[radAfter].idAbs() < 10)
6592 return spinRadAfter;
6595 if ( !event[radAfter].isFinal() && radBeforeFlav == 21
6596 && event[emtAfter].idAbs() < 10)
6598 return spinEmtAfter;
6617 bool History::connectRadiator( Particle& Radiator,
const int RadType,
6618 const Particle& Recoiler,
const int RecType,
6619 const Event& event ) {
6622 Radiator.cols( -1, -1 );
6626 if ( Radiator.colType() == -1 ) {
6629 if ( RadType + RecType == 2 )
6630 Radiator.cols( 0, Recoiler.col());
6631 else if ( RadType + RecType == 0 )
6632 Radiator.cols( 0, Recoiler.acol());
6640 for (
int i = 0; i <
event.size(); ++i) {
6641 int col =
event[i].col();
6642 int acl =
event[i].acol();
6644 if ( event[i].isFinal()) {
6646 if ( acl > 0 && FindCol(acl,i,0,event,1,
true) == 0
6647 && FindCol(acl,i,0,event,2,
true) == 0 )
6648 Radiator.acol(event[i].acol());
6651 if ( col > 0 && FindCol(col,i,0,event,1,
true) == 0
6652 && FindCol(col,i,0,event,2,
true) == 0 )
6653 Radiator.acol(event[i].col());
6658 }
else if ( Radiator.colType() == 1 ) {
6661 if ( RadType + RecType == 2 )
6662 Radiator.cols( Recoiler.acol(), 0);
6664 else if ( RadType + RecType == 0 )
6665 Radiator.cols( Recoiler.col(), 0);
6674 for (
int i = 0; i <
event.size(); ++i) {
6675 int col =
event[i].col();
6676 int acl =
event[i].acol();
6678 if ( event[i].isFinal()) {
6680 if ( col > 0 && FindCol(col,i,0,event,1,
true) == 0
6681 && FindCol(col,i,0,event,2,
true) == 0)
6682 Radiator.col(event[i].col());
6685 if ( acl > 0 && FindCol(acl,i,0,event,1,
true) == 0
6686 && FindCol(acl,i,0,event,2,
true) == 0)
6687 Radiator.col(event[i].acol());
6693 }
else if ( Radiator.colType() == 2 ) {
6699 for (
int i = 0; i <
event.size(); ++i) {
6700 int col =
event[i].col();
6701 int acl =
event[i].acol();
6704 if ( event[i].isFinal()) {
6705 if ( col > 0 && FindCol(col,iEx,0,event,1,
true) == 0
6706 && FindCol(col,iEx,0,event,2,
true) == 0) {
6707 if (Radiator.status() < 0 ) Radiator.col(event[i].col());
6708 else Radiator.acol(event[i].col());
6710 if ( acl > 0 && FindCol(acl,iEx,0,event,2,
true) == 0
6711 && FindCol(acl,iEx,0,event,1,
true) == 0 ) {
6712 if (Radiator.status() < 0 ) Radiator.acol(event[i].acol());
6713 else Radiator.col(event[i].acol());
6716 if ( col > 0 && FindCol(col,iEx,0,event,1,
true) == 0
6717 && FindCol(col,iEx,0,event,2,
true) == 0) {
6718 if (Radiator.status() < 0 ) Radiator.acol(event[i].col());
6719 else Radiator.col(event[i].col());
6721 if ( acl > 0 && (FindCol(acl,iEx,0,event,2,
true) == 0
6722 && FindCol(acl,iEx,0,event,1,
true) == 0)) {
6723 if (Radiator.status() < 0 ) Radiator.col(event[i].acol());
6724 else Radiator.acol(event[i].acol());
6731 if (Radiator.col() < 0 || Radiator.acol() < 0)
return false;
6753 int History::FindCol(
int col,
int iExclude1,
int iExclude2,
6754 const Event& event,
int type,
bool isHardIn) {
6756 bool isHard = isHardIn;
6761 for(
int n = 0; n <
event.size(); ++n) {
6762 if ( n != iExclude1 && n != iExclude2
6763 && event[n].colType() != 0
6764 &&(
event[n].status() > 0
6765 ||
event[n].status() == -21) ) {
6766 if ( event[n].acol() == col ) {
6770 if ( event[n].col() == col ) {
6779 for(
int n = 0; n <
event.size(); ++n) {
6780 if ( n != iExclude1 && n != iExclude2
6781 && event[n].colType() != 0
6782 &&(
event[n].status() == 43
6783 ||
event[n].status() == 51
6784 ||
event[n].status() == -41
6785 ||
event[n].status() == -42) ) {
6786 if ( event[n].acol() == col ) {
6790 if ( event[n].col() == col ) {
6798 if ( type == 1 && index < 0)
return abs(index);
6799 if ( type == 2 && index > 0)
return abs(index);
6813 int History::FindParticle(
const Particle& particle,
const Event& event,
6814 bool checkStatus ) {
6818 for (
int i =
int(event.size()) - 1; i > 0; --i )
6819 if ( event[i].
id() == particle.id()
6820 &&
event[i].colType() == particle.colType()
6821 &&
event[i].chargeType() == particle.chargeType()
6822 &&
event[i].col() == particle.col()
6823 &&
event[i].acol() == particle.acol()
6824 &&
event[i].charge() == particle.charge() ) {
6829 if ( checkStatus && event[index].status() != particle.status() )
6844 int History::getRadBeforeCol(
const int rad,
const int emt,
6845 const Event& event) {
6848 int type = (
event[rad].isFinal()) ? 1 :-1;
6850 int radBeforeFlav = getRadBeforeFlav(rad,emt,event);
6852 int radBeforeCol = -1;
6854 if (radBeforeFlav == 21) {
6857 if (type == 1 && event[emt].
id() != 21) {
6858 radBeforeCol = (
event[rad].col() > 0)
6859 ? event[rad].col() :
event[emt].col();
6861 }
else if (type == -1 && event[emt].
id() != 21) {
6862 radBeforeCol = (
event[rad].col() > 0)
6863 ? event[rad].col() :
event[emt].acol();
6865 }
else if (type == 1 && event[emt].
id() == 21) {
6868 int colRemove = (
event[rad].col() ==
event[emt].acol())
6869 ? event[rad].col() :
event[rad].acol();
6870 radBeforeCol = (
event[rad].col() == colRemove)
6871 ? event[emt].col() :
event[rad].col();
6873 }
else if (type == -1 && event[emt].
id() == 21) {
6876 int colRemove = (
event[rad].col() ==
event[emt].col())
6877 ? event[rad].col() :
event[rad].acol();
6878 radBeforeCol = (
event[rad].col() == colRemove)
6879 ? event[emt].acol() :
event[rad].col();
6883 }
else if ( radBeforeFlav != 21 && radBeforeFlav > 0) {
6886 if (type == 1 && event[emt].
id() != 21) {
6889 int colRemove = (
event[rad].col() ==
event[emt].acol())
6890 ? event[rad].acol() : 0;
6891 radBeforeCol = (
event[rad].col() == colRemove)
6892 ? event[emt].col() :
event[rad].col();
6894 }
else if (type == 1 && event[emt].
id() == 21) {
6897 int colRemove = (
event[rad].col() ==
event[emt].acol())
6898 ? event[rad].col() : 0;
6899 radBeforeCol = (
event[rad].col() == colRemove)
6900 ? event[emt].col() :
event[rad].col();
6902 }
else if (type == -1 && event[emt].
id() != 21) {
6905 int colRemove = (
event[rad].col() ==
event[emt].col())
6906 ? event[rad].col() : 0;
6907 radBeforeCol = (
event[rad].col() == colRemove)
6908 ? event[emt].acol() :
event[rad].col();
6910 }
else if (type == -1 && event[emt].
id() == 21) {
6913 int colRemove = (
event[rad].col() ==
event[emt].col())
6914 ? event[rad].col() : 0;
6915 radBeforeCol = (
event[rad].col() == colRemove)
6916 ? event[emt].acol() :
event[rad].col();
6923 return radBeforeCol;
6936 int History::getRadBeforeAcol(
const int rad,
const int emt,
6937 const Event& event) {
6940 int type = (
event[rad].isFinal()) ? 1 :-1;
6943 int radBeforeFlav = getRadBeforeFlav(rad,emt,event);
6945 int radBeforeAcl = -1;
6947 if (radBeforeFlav == 21) {
6950 if (type == 1 && event[emt].
id() != 21) {
6951 radBeforeAcl = (
event[rad].acol() > 0)
6952 ? event[rad].acol() :
event[emt].acol();
6954 }
else if (type == -1 && event[emt].
id() != 21) {
6955 radBeforeAcl = (
event[rad].acol() > 0)
6956 ? event[rad].acol() :
event[emt].col();
6958 }
else if (type == 1 && event[emt].
id() == 21) {
6961 int colRemove = (
event[rad].col() ==
event[emt].acol())
6962 ? event[rad].col() :
event[rad].acol();
6963 radBeforeAcl = (
event[rad].acol() == colRemove)
6964 ? event[emt].acol() :
event[rad].acol();
6966 }
else if (type == -1 && event[emt].
id() == 21) {
6969 int colRemove = (
event[rad].col() ==
event[emt].col())
6970 ? event[rad].col() :
event[rad].acol();
6971 radBeforeAcl = (
event[rad].acol() == colRemove)
6972 ? event[emt].col() :
event[rad].acol();
6976 }
else if ( radBeforeFlav != 21 && radBeforeFlav < 0) {
6979 if (type == 1 && event[emt].
id() != 21) {
6982 int colRemove = (
event[rad].col() ==
event[emt].acol())
6983 ? event[rad].acol() : 0;
6984 radBeforeAcl = (
event[rad].acol() == colRemove)
6985 ? event[emt].acol() :
event[rad].acol();
6987 }
else if (type == 1 && event[emt].
id() == 21) {
6990 int colRemove = (
event[rad].acol() ==
event[emt].col())
6991 ? event[rad].acol() : 0;
6992 radBeforeAcl = (
event[rad].acol() == colRemove)
6993 ? event[emt].acol() :
event[rad].acol();
6995 }
else if (type == -1 && event[emt].
id() != 21) {
6998 int colRemove = (
event[rad].acol() ==
event[emt].acol())
6999 ? event[rad].acol() : 0;
7000 radBeforeAcl = (
event[rad].acol() == colRemove)
7001 ? event[emt].col() :
event[rad].acol();
7003 }
else if (type == -1 && event[emt].
id() == 21) {
7006 int colRemove = (
event[rad].acol() ==
event[emt].acol())
7007 ? event[rad].acol() : 0;
7008 radBeforeAcl = (
event[rad].acol() == colRemove)
7009 ? event[emt].col() :
event[rad].acol();
7016 return radBeforeAcl;
7028 int History::getColPartner(
const int in,
const Event& event) {
7030 if (event[in].col() == 0)
return 0;
7034 partner = FindCol(event[in].col(),in,0,event,1,
true);
7037 partner = FindCol(event[in].col(),in,0,event,2,
true);
7052 int History::getAcolPartner(
const int in,
const Event& event) {
7054 if (event[in].acol() == 0)
return 0;
7058 partner = FindCol(event[in].acol(),in,0,event,2,
true);
7061 partner = FindCol(event[in].acol(),in,0,event,1,
true);
7078 vector<int> History::getReclusteredPartners(
const int rad,
const int emt,
7079 const Event& event) {
7082 int type =
event[rad].isFinal() ? 1 : -1;
7084 int radBeforeCol = getRadBeforeCol(rad, emt, event);
7085 int radBeforeAcl = getRadBeforeAcol(rad, emt, event);
7087 vector<int> partners;
7093 for(
int i=0; i < int(event.size()); ++i) {
7095 if ( i != emt && i != rad
7096 && event[i].status() == -21
7097 &&
event[i].col() > 0
7098 &&
event[i].col() == radBeforeCol)
7099 partners.push_back(i);
7101 if ( i != emt && i != rad
7102 && event[i].isFinal()
7103 && event[i].acol() > 0
7104 && event[i].acol() == radBeforeCol)
7105 partners.push_back(i);
7107 if ( i != emt && i != rad
7108 && event[i].status() == -21
7109 && event[i].acol() > 0
7110 && event[i].acol() == radBeforeAcl)
7111 partners.push_back(i);
7113 if ( i != emt && i != rad
7114 && event[i].isFinal()
7115 && event[i].col() > 0
7116 && event[i].col() == radBeforeAcl)
7117 partners.push_back(i);
7122 for(
int i=0; i < int(event.size()); ++i) {
7124 if ( i != emt && i != rad
7125 && event[i].status() == -21
7126 &&
event[i].acol() > 0
7127 &&
event[i].acol() == radBeforeCol)
7128 partners.push_back(i);
7130 if ( i != emt && i != rad
7131 && event[i].isFinal()
7132 && event[i].col() > 0
7133 && event[i].col() == radBeforeCol)
7134 partners.push_back(i);
7136 if ( i != emt && i != rad
7137 && event[i].status() == -21
7138 && event[i].col() > 0
7139 && event[i].col() == radBeforeAcl)
7140 partners.push_back(i);
7142 if ( i != emt && i != rad
7143 && event[i].isFinal()
7144 && event[i].acol() > 0
7145 && event[i].acol() == radBeforeAcl)
7146 partners.push_back(i);
7173 bool History::getColSinglet(
const int flavType,
const int iParton,
7174 const Event& event, vector<int>& exclude, vector<int>& colSinglet) {
7177 if (iParton < 0)
return false;
7185 for(
int i=0; i < int(event.size()); ++i)
7186 if ( event[i].isFinal() &&
event[i].colType() != 0)
7191 int nExclude = int(exclude.size());
7192 int nInitExclude = 0;
7193 if (!event[exclude[2]].isFinal())
7195 if (!event[exclude[3]].isFinal())
7199 if (nFinal == nExclude - nInitExclude)
7209 colSinglet.push_back(iParton);
7211 exclude.push_back(iParton);
7214 colP = getColPartner(iParton,event);
7217 colP = getAcolPartner(iParton,event);
7220 for(
int i = 0; i < int(exclude.size()); ++i)
7221 if (colP == exclude[i])
7225 return getColSinglet(flavType,colP,event,exclude,colSinglet);
7236 bool History::isColSinglet(
const Event& event,
7237 vector<int> system ) {
7240 for(
int i=0; i < int(system.size()); ++i ) {
7243 && (event[system[i]].colType() == 1
7244 ||
event[system[i]].colType() == 2) ) {
7245 for(
int j=0; j < int(system.size()); ++j)
7249 && event[system[i]].col() ==
event[system[j]].acol()) {
7258 && (event[system[i]].colType() == -1
7259 || event[system[i]].colType() == 2) ) {
7260 for(
int j=0; j < int(system.size()); ++j)
7264 && event[system[i]].acol() ==
event[system[j]].col()) {
7276 bool isColSing =
true;
7277 for(
int i=0; i < int(system.size()); ++i)
7278 if ( system[i] != 0 )
7296 bool History::isFlavSinglet(
const Event& event,
7297 vector<int> system,
int flav) {
7302 for(
int i=0; i < int(system.size()); ++i)
7303 if ( system[i] > 0 ) {
7304 for(
int j=0; j < int(system.size()); ++j) {
7308 if ( event[i].idAbs() != 21
7309 &&
event[i].idAbs() != 22
7310 &&
event[i].idAbs() != 23
7311 &&
event[i].idAbs() != 24
7314 &&
event[system[i]].isFinal()
7315 &&
event[system[j]].isFinal()
7316 &&
event[system[i]].id() == -1*
event[system[j]].id()) {
7319 if (abs(flav) > 0 &&
event[system[i]].idAbs() != flav)
7329 if ( event[i].idAbs() != 21
7330 &&
event[i].idAbs() != 22
7331 &&
event[i].idAbs() != 23
7332 &&
event[i].idAbs() != 24
7335 && ( ( !
event[system[i]].isFinal() &&
event[system[j]].isFinal())
7336 ||( !event[system[j]].isFinal() &&
event[system[i]].isFinal()) )
7337 &&
event[system[i]].id() ==
event[system[j]].id()) {
7340 if (abs(flav) > 0 &&
event[system[i]].idAbs() != flav)
7353 bool isFlavSing =
true;
7354 for(
int i=0; i < int(system.size()); ++i)
7355 if ( system[i] != 0 )
7370 bool History::allowedClustering(
int rad,
int emt,
int rec,
int partner,
7371 const Event& event ) {
7374 bool allowed =
true;
7379 bool isSing = isSinglett(rad,emt,partner,event);
7380 int type = (
event[rad].isFinal()) ? 1 :-1;
7382 int radBeforeFlav = getRadBeforeFlav(rad,emt,event);
7384 int radBeforeCol = getRadBeforeCol(rad,emt,event);
7385 int radBeforeAcl = getRadBeforeAcol(rad,emt,event);
7387 vector<int> radBeforeColP = getReclusteredPartners(rad, emt, event);
7390 int nPartonInHard = 0;
7391 for(
int i=0; i < int(event.size()); ++i)
7393 if ( event[i].isFinal()
7394 &&
event[i].colType() != 0
7395 && mergingHooksPtr->hardProcess->matchesAnyOutgoing(i, event) )
7401 for(
int i=0; i < int(event.size()); ++i)
7403 if ( i!=emt && i!=rad && i!=rec
7404 && event[i].isFinal()
7405 &&
event[i].colType() != 0
7406 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i, event) )
7410 int nInitialPartons = 0;
7411 for(
int i=0; i < int(event.size()); ++i)
7412 if ( event[i].status() == -21
7413 &&
event[i].colType() != 0 )
7418 for(
int i=0; i < int(event.size()); ++i)
7419 if ( event[i].isFinal()
7420 &&(
event[i].id() == 22
7421 ||
event[i].id() == 23
7422 ||
event[i].id() == 24
7423 ||(
event[i].idAbs() > 10 &&
event[i].idAbs() < 20)
7424 ||(event[i].idAbs() > 10 &&
event[i].idAbs() < 20)
7425 ||(event[i].idAbs() > 1000010 &&
event[i].idAbs() < 1000020)
7426 ||(event[i].idAbs() > 2000010 &&
event[i].idAbs() < 2000020) ))
7433 int nFinalQuark = 0;
7435 int nFinalQuarkExc = 0;
7436 for(
int i=0; i < int(event.size()); ++i) {
7437 if (i !=rad && i != emt && i != rec) {
7438 if (event[i].isFinal() && abs(event[i].colType()) == 1 ) {
7439 if ( !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,event) )
7448 if (event[rec].isFinal() &&
event[rec].isQuark()) nFinalQuark++;
7450 if ( event[rad].isFinal()
7451 && abs(particleDataPtr->colType(radBeforeFlav)) == 1) nFinalQuark++;
7454 int nInitialQuark = 0;
7456 if (event[rec].isFinal()) {
7457 if (event[3].isQuark()) nInitialQuark++;
7458 if (event[4].isQuark()) nInitialQuark++;
7460 int iOtherIn = (rec == 3) ? 4 : 3;
7461 if (event[rec].isQuark()) nInitialQuark++;
7462 if (event[iOtherIn].isQuark()) nInitialQuark++;
7466 if (event[rec].isQuark()) nInitialQuark++;
7468 if (abs(radBeforeFlav) < 10) nInitialQuark++;
7474 int proton[] = {1,2,3,4,5,21,22,23,24};
7475 bool isInProton =
false;
7476 for(
int i=0; i < 9; ++i)
7477 if (abs(radBeforeFlav) == proton[i]) isInProton =
true;
7478 if (type == -1 && !isInProton)
return false;
7481 vector<int> unmatchedCol;
7482 vector<int> unmatchedAcl;
7484 for (
int i = 0; i <
event.size(); ++i)
7485 if ( i != emt && i != rad
7486 && (event[i].isFinal() || event[i].status() == -21)
7487 &&
event[i].colType() != 0 ) {
7489 int colP = getColPartner(i,event);
7490 int aclP = getAcolPartner(i,event);
7492 if (event[i].col() > 0
7493 && (colP == emt || colP == rad || colP == 0) )
7494 unmatchedCol.push_back(i);
7495 if (event[i].acol() > 0
7496 && (aclP == emt || aclP == rad || aclP == 0) )
7497 unmatchedAcl.push_back(i);
7503 if (
int(unmatchedCol.size()) +
int(unmatchedAcl.size()) > 2)
7510 if ( isSing && event[rec].isQuark()
7511 && abs(particleDataPtr->colType(radBeforeFlav)) == 1)
7515 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(emt,event) ) {
7519 bool canReplace = mergingHooksPtr->hardProcess->findOtherCandidates(emt,
7521 if (canReplace) allowed =
true;
7522 else allowed =
false;
7527 if ( mergingHooksPtr->hardProcess->matchesAnyOutgoing(rad,event)
7528 &&
event[rad].id() != radBeforeFlav )
7533 if (nFinalEW != 0 && nInitialQuark == 0
7534 && nFinalQuark == 0 && nFinalQuarkExc == 0)
7537 if ( (nInitialQuark + nFinalQuark + nFinalQuarkExc)%2 != 0 )
7547 for(
int i=0; i < int(event.size()); ++i)
7548 if ( i!=emt && i!=rad && i!=rec
7549 && (event[i].mother1() == 1 ||
event[i].mother1() == 2))
7550 in.push_back(event[i].
id());
7551 if (!event[rad].isFinal()) in.push_back(radBeforeFlav);
7552 if (!event[rec].isFinal()) in.push_back(event[rec].
id());
7554 for(
int i=0; i < int(event.size()); ++i)
7555 if ( i!=emt && i!=rad && i!=rec && event[i].isFinal())
7556 out.push_back(event[i].id());
7557 if (event[rad].isFinal()) out.push_back(radBeforeFlav);
7558 if (event[rec].isFinal()) out.push_back(event[rec].
id());
7559 if (event[3].col() == event[4].acol()
7560 && event[3].acol() == event[4].col()
7561 && !mergingHooksPtr->allowEffectiveVertex( in, out)
7562 && nFinalQuark == 0){
7565 int nTripletts = abs(event[rec].colType())
7566 + abs(particleDataPtr->colType(radBeforeFlav));
7567 if (event[3].isGluon()) allowed =
false;
7568 else if (nTripletts != 2 && nFinalQuarkExc%2 == 0) allowed =
false;
7572 if ( abs((event[rad].p()+type*event[emt].p()+event[rec].p()).pz())
7573 > (
event[rad].p()+type*
event[emt].p()+
event[rec].p()).e()
7575 && (
event[rad].p()-
event[emt].p()+
event[rec].p()).m2Calc() < 0.) ){
7580 if (event[emt].
id() == 21)
return allowed;
7583 if (event[emt].
id() == 1000021)
return allowed;
7586 vector<int> outgoingParticles;
7587 int nOut1 = int(mergingHooksPtr->hardProcess->PosOutgoing1.size());
7588 for (
int i=0; i < nOut1; ++i ) {
7589 int iPos = mergingHooksPtr->hardProcess->PosOutgoing1[i];
7590 outgoingParticles.push_back(
7591 mergingHooksPtr->hardProcess->state[iPos].id() );
7593 int nOut2 = int(mergingHooksPtr->hardProcess->PosOutgoing2.size());
7594 for (
int i=0; i < nOut2; ++i ) {
7595 int iPos = mergingHooksPtr->hardProcess->PosOutgoing2[i];
7596 outgoingParticles.push_back(
7597 mergingHooksPtr->hardProcess->state[iPos].id() );
7604 vector<int> iInQuarkFlav;
7605 for(
int i=0; i < int(event.size()); ++i)
7607 if ( i != emt && i != rad
7608 && event[i].status() == -21
7609 &&
event[i].idAbs() ==
event[emt].idAbs() )
7610 iInQuarkFlav.push_back(i);
7613 vector<int> iOutQuarkFlav;
7614 for(
int i=0; i < int(event.size()); ++i)
7616 if ( i != emt && i != rad
7617 && event[i].isFinal()
7618 &&
event[i].idAbs() ==
event[emt].idAbs() ) {
7622 bool matchOut =
false;
7623 for (
int j = 0; j < int(outgoingParticles.size()); ++j)
7624 if ( event[i].idAbs() == abs(outgoingParticles[j])) {
7626 outgoingParticles[j] = 99;
7628 if (!matchOut) iOutQuarkFlav.push_back(i);
7633 int nInQuarkFlav = int(iInQuarkFlav.size());
7634 int nOutQuarkFlav = int(iOutQuarkFlav.size());
7639 if ( event[partner].isFinal()
7640 &&
event[partner].id() == 21
7641 && radBeforeFlav == 21
7642 &&
event[partner].col() == radBeforeCol
7643 &&
event[partner].acol() == radBeforeAcl)
7647 if (nInQuarkFlav + nOutQuarkFlav == 0)
7654 vector<int> partons;
7655 for(
int i=0; i < int(event.size()); ++i)
7657 if ( i!=emt && i!=rad
7658 && event[i].colType() != 0
7659 && (
event[i].isFinal() ||
event[i].status() == -21) ) {
7661 partons.push_back(i);
7663 if (event[i].colType() == 2)
7665 else if (event[i].colType() == 1)
7667 else if (event[i].colType() == -1)
7673 bool isFSRg2qq = ((type == 1) && (event[rad].
id() == -1*
event[emt].id()) );
7674 bool isISRg2qq = ((type ==-1) && (event[rad].
id() ==
event[emt].id()) );
7679 if ( (isFSRg2qq || isISRg2qq)
7680 && int(quark.size()) +
int(antiq.size())
7681 +
int(gluon.size()) > nPartonInHard ) {
7683 vector<int> colours;
7684 vector<int> anticolours;
7689 colours.push_back(radBeforeCol);
7690 anticolours.push_back(radBeforeAcl);
7692 colours.push_back(radBeforeAcl);
7693 anticolours.push_back(radBeforeCol);
7696 for(
int i=0; i < int(gluon.size()); ++i)
7697 if (event[gluon[i]].isFinal()) {
7698 colours.push_back(event[gluon[i]].col());
7699 anticolours.push_back(event[gluon[i]].acol());
7701 colours.push_back(event[gluon[i]].acol());
7702 anticolours.push_back(event[gluon[i]].col());
7707 for(
int i=0; i < int(colours.size()); ++i)
7708 for(
int j=0; j < int(anticolours.size()); ++j)
7709 if (colours[i] > 0 && anticolours[j] > 0
7710 && colours[i] == anticolours[j]) {
7718 bool allMatched =
true;
7719 for(
int i=0; i < int(colours.size()); ++i)
7720 if (colours[i] != 0)
7722 for(
int i=0; i < int(anticolours.size()); ++i)
7723 if (anticolours[i] != 0)
7731 for(
int i=0; i < int(quark.size()); ++i)
7732 if ( event[quark[i]].isFinal()
7733 && mergingHooksPtr->hardProcess->matchesAnyOutgoing(quark[i], event) )
7734 colours.push_back(event[quark[i]].col());
7736 for(
int i=0; i < int(antiq.size()); ++i)
7737 if ( event[antiq[i]].isFinal()
7738 && mergingHooksPtr->hardProcess->matchesAnyOutgoing(antiq[i], event) )
7739 anticolours.push_back(event[antiq[i]].acol());
7743 for(
int i=0; i < int(colours.size()); ++i)
7745 for(
int j=0; j < int(anticolours.size()); ++j)
7746 if (colours[i] > 0 && anticolours[j] > 0
7747 && colours[i] == anticolours[j]) {
7754 for (
int i=0; i < int(quark.size()); ++i )
7755 if ( !mergingHooksPtr->hardProcess->matchesAnyOutgoing( quark[i],
7758 for (
int i=0; i < int(antiq.size()); ++i )
7759 if ( !mergingHooksPtr->hardProcess->matchesAnyOutgoing( antiq[i],
7762 for(
int i=0; i < int(gluon.size()); ++i)
7763 if ( event[gluon[i]].isFinal() )
7771 for(
int i=0; i < int(colours.size()); ++i)
7772 if (colours[i] != 0)
7774 for(
int i=0; i < int(anticolours.size()); ++i)
7775 if (anticolours[i] != 0)
7778 if (allMatched && nNotInHard > 0)
7785 if (isFSRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) {
7789 for(
int i=0; i < int(gluon.size()); ++i) {
7790 if (!event[gluon[i]].isFinal()
7791 &&
event[gluon[i]].col() == radBeforeCol
7792 &&
event[gluon[i]].acol() == radBeforeAcl)
7798 for(
int i=0; i < int(gluon.size()); ++i) {
7799 if (event[gluon[i]].isFinal()
7800 &&
event[gluon[i]].col() == radBeforeAcl
7801 &&
event[gluon[i]].acol() == radBeforeCol)
7807 if (
int(radBeforeColP.size()) == 2
7808 && event[radBeforeColP[0]].isFinal()
7809 &&
event[radBeforeColP[1]].isFinal()
7810 &&
event[radBeforeColP[0]].id() == -1*
event[radBeforeColP[1]].id() ) {
7814 if (nInitialPartons > 0)
7821 if (
int(radBeforeColP.size()) == 2
7822 && (( event[radBeforeColP[0]].status() == -21
7823 && event[radBeforeColP[1]].isFinal())
7824 ||(
event[radBeforeColP[0]].isFinal()
7825 &&
event[radBeforeColP[1]].status() == -21))
7826 &&
event[radBeforeColP[0]].id() ==
event[radBeforeColP[1]].id() ) {
7833 int incoming = (
event[radBeforeColP[0]].isFinal())
7834 ? radBeforeColP[1] : radBeforeColP[0];
7835 int outgoing = (
event[radBeforeColP[0]].isFinal())
7836 ? radBeforeColP[0] : radBeforeColP[1];
7839 bool clusPossible =
false;
7840 for(
int i=0; i < int(event.size()); ++i)
7841 if ( i != emt && i != rad
7842 && i != incoming && i != outgoing
7843 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,event) ) {
7845 if ( event[i].status() == -21
7846 && (
event[i].id() ==
event[outgoing].id()
7847 ||
event[i].id() == -1*
event[incoming].id()) )
7848 clusPossible =
true;
7850 if ( event[i].isFinal()
7851 && (
event[i].id() == -1*
event[outgoing].id()
7852 ||
event[i].id() ==
event[incoming].id()) )
7853 clusPossible =
true;
7864 vector<int> excludeIn1;
7865 for(
int i=0; i < 4; ++i)
7866 excludeIn1.push_back(0);
7867 vector<int> colSingletIn1;
7868 int flavIn1Type = (
event[incoming].id() > 0) ? 1 : -1;
7870 bool isColSingIn1 = getColSinglet(flavIn1Type,incoming,event,
7871 excludeIn1,colSingletIn1);
7873 bool isFlavSingIn1 = isFlavSinglet(event,colSingletIn1);
7878 int incoming2 = (incoming == 3) ? 4 : 3;
7879 vector<int> excludeIn2;
7880 for(
int i=0; i < 4; ++i)
7881 excludeIn2.push_back(0);
7882 vector<int> colSingletIn2;
7883 int flavIn2Type = (
event[incoming2].id() > 0) ? 1 : -1;
7885 bool isColSingIn2 = getColSinglet(flavIn2Type,incoming2,event,
7886 excludeIn2,colSingletIn2);
7888 bool isFlavSingIn2 = isFlavSinglet(event,colSingletIn2);
7892 && (!isColSingIn1 || !isFlavSingIn1
7893 || !isColSingIn2 || !isFlavSingIn2))
7905 int flav =
event[emt].id();
7906 vector<int> exclude;
7907 exclude.push_back(emt);
7908 exclude.push_back(rad);
7909 exclude.push_back(radBeforeColP[0]);
7910 exclude.push_back(radBeforeColP[1]);
7911 vector<int> colSinglet;
7915 for(
int i=0; i < int(event.size()); ++i)
7920 && i != radBeforeColP[0]
7921 && i != radBeforeColP[1]
7922 && event[i].isFinal() ) {
7924 if (event[i].
id() == flav) {
7930 int flavType = (iOther > 0 &&
event[iOther].id() > 0) ? 1
7931 : (iOther > 0) ? -1 : 0;
7933 bool isColSing = getColSinglet(flavType,iOther,event,exclude,colSinglet);
7935 bool isFlavSing = isFlavSinglet(event,colSinglet);
7939 bool isHardSys =
true;
7940 for(
int i=0; i < int(colSinglet.size()); ++i)
7941 isHardSys = mergingHooksPtr->hardProcess->matchesAnyOutgoing(
7942 colSinglet[i], event);
7947 if (isColSing && isFlavSing && !isHardSys) {
7953 vector<int> allFinal;
7954 for(
int i=0; i < int(event.size()); ++i)
7955 if ( event[i].isFinal() )
7956 allFinal.push_back(i);
7959 bool isFullColSing = isColSinglet(event,allFinal);
7961 bool isFullFlavSing = isFlavSinglet(event,allFinal,flav);
7966 if (!isFullColSing || !isFullFlavSing)
7973 if (isISRg2qq && nInQuarkFlav + nOutQuarkFlav > 0) {
7977 for(
int i=0; i < int(gluon.size()); ++i) {
7978 if (event[gluon[i]].isFinal()
7979 &&
event[gluon[i]].col() == radBeforeCol
7980 &&
event[gluon[i]].acol() == radBeforeAcl)
7986 for(
int i=0; i < int(gluon.size()); ++i) {
7987 if (event[gluon[i]].status() == -21
7988 &&
event[gluon[i]].acol() == radBeforeCol
7989 &&
event[gluon[i]].col() == radBeforeAcl)
7995 if (
int(radBeforeColP.size()) == 2
7996 && event[radBeforeColP[0]].isFinal()
7997 &&
event[radBeforeColP[1]].isFinal()
7998 &&
event[radBeforeColP[0]].id() == -1*
event[radBeforeColP[1]].id() ) {
8005 bool clusPossible =
false;
8006 for(
int i=0; i < int(event.size()); ++i)
8007 if ( i != emt && i != rad
8008 && i != radBeforeColP[0]
8009 && i != radBeforeColP[1]
8010 && !mergingHooksPtr->hardProcess->matchesAnyOutgoing(i,event) ) {
8011 if (event[i].status() == -21
8012 && (
event[radBeforeColP[0]].id() ==
event[i].id()
8013 ||
event[radBeforeColP[1]].id() ==
event[i].id() ))
8015 clusPossible =
true;
8016 if (event[i].isFinal()
8017 && (
event[radBeforeColP[0]].id() == -1*
event[i].id()
8018 ||
event[radBeforeColP[1]].id() == -1*
event[i].id() ))
8019 clusPossible =
true;
8031 vector<int> excludeIn1;
8032 for(
int i=0; i < 4; ++i)
8033 excludeIn1.push_back(0);
8034 vector<int> colSingletIn1;
8035 int flavIn1Type = (
event[incoming1].id() > 0) ? 1 : -1;
8037 bool isColSingIn1 = getColSinglet(flavIn1Type,incoming1,event,
8038 excludeIn1,colSingletIn1);
8040 bool isFlavSingIn1 = isFlavSinglet(event,colSingletIn1);
8046 vector<int> excludeIn2;
8047 for(
int i=0; i < 4; ++i)
8048 excludeIn2.push_back(0);
8049 vector<int> colSingletIn2;
8050 int flavIn2Type = (
event[incoming2].id() > 0) ? 1 : -1;
8052 bool isColSingIn2 = getColSinglet(flavIn2Type,incoming2,event,
8053 excludeIn2,colSingletIn2);
8055 bool isFlavSingIn2 = isFlavSinglet(event,colSingletIn2);
8059 && (!isColSingIn1 || !isFlavSingIn1
8060 || !isColSingIn2 || !isFlavSingIn2))
8079 bool History::isSinglett(
int rad,
int emt,
int rec,
const Event& event ) {
8081 int radCol =
event[rad].col();
8082 int emtCol =
event[emt].col();
8083 int recCol =
event[rec].col();
8084 int radAcl =
event[rad].acol();
8085 int emtAcl =
event[emt].acol();
8086 int recAcl =
event[rec].acol();
8087 int recType =
event[rec].isFinal() ? 1 : -1;
8089 bool isSing =
false;
8091 if ( ( recType == -1
8092 && radCol + emtCol == recCol && radAcl + emtAcl == recAcl)
8094 && radCol + emtCol == recAcl && radAcl + emtAcl == recCol) )
8110 bool History::validEvent(
const Event& event ) {
8113 bool validColour =
true;
8114 for (
int i = 0; i <
event.size(); ++i)
8116 if ( event[i].isFinal() &&
event[i].colType() == 1
8118 && ( FindCol(event[i].col(),i,0,event,1,
true) == 0
8120 && FindCol(event[i].col(),i,0,event,2,
true) == 0 )) {
8121 validColour =
false;
8124 }
else if ( event[i].isFinal() &&
event[i].colType() == -1
8126 && ( FindCol(event[i].acol(),i,0,event,2,
true) == 0
8128 && FindCol(event[i].acol(),i,0,event,1,
true) == 0 )) {
8129 validColour =
false;
8132 }
else if ( event[i].isFinal() &&
event[i].colType() == 2
8134 && ( FindCol(event[i].col(),i,0,event,1,
true) == 0
8136 && FindCol(event[i].col(),i,0,event,2,
true) == 0 )
8138 && ( FindCol(event[i].acol(),i,0,event,2,
true) == 0
8140 && FindCol(event[i].acol(),i,0,event,1,
true) == 0 )) {
8141 validColour =
false;
8146 bool validCharge =
true;
8147 double initCharge =
event[3].charge() +
event[4].charge();
8148 double finalCharge = 0.0;
8149 for(
int i = 0; i <
event.size(); ++i)
8150 if (event[i].isFinal()) finalCharge += event[i].charge();
8151 if (abs(initCharge-finalCharge) > 1e-12) validCharge =
false;
8153 return (validColour && validCharge);
8162 bool History::equalClustering( Clustering clus1 , Clustering clus2 ) {
8163 return ( (clus1.emittor == clus2.emittor)
8164 && (clus1.emitted == clus2.emitted)
8165 && (clus1.recoiler == clus2.recoiler)
8166 && (clus1.partner == clus2.partner)
8167 && (clus1.pT() == clus2.pT())
8168 && (clus1.spinRadBef == clus2.spinRadBef)
8169 && (clus1.spinRad == clus2.spinRad)
8170 && (clus1.spinEmt == clus2.spinEmt)
8171 && (clus1.spinRec == clus2.spinRec)
8172 && (clus1.flavRadBef == clus2.flavRadBef));
8181 double History::choseHardScale(
const Event& event )
const {
8184 double mHat = (
event[3].p() +
event[4].p()).mCalc();
8191 for(
int i = 0; i <
event.size(); ++i)
8192 if ( event[i].isFinal() ) {
8195 if ( event[i].idAbs() == 23
8196 ||
event[i].idAbs() == 24 ) {
8199 mBos +=
event[i].m();
8201 }
else if ( abs(event[i].status()) == 22
8202 && ( event[i].idAbs() == 23
8203 || event[i].idAbs() == 24 )) {
8205 mBos +=
event[i].m();
8209 if ( nBosons > 0 && (nFinal + nFinBos*2) <= 3)
8210 return (mBos /
double(nBosons));
8221 int History::getCurrentFlav(
const int side)
const {
8222 int in = (side == 1) ? 3 : 4;
8223 return state[in].id();
8228 double History::getCurrentX(
const int side)
const {
8229 int in = (side == 1) ? 3 : 4;
8230 return ( 2.*state[in].e()/state[0].e() );
8235 double History::getCurrentZ(
const int rad,
8236 const int rec,
const int emt,
int idRadBef)
const {
8238 int type = state[rad].isFinal() ? 1 : -1;
8243 Vec4 radAfterBranch(state[rad].p());
8244 Vec4 recAfterBranch(state[rec].p());
8245 Vec4 emtAfterBranch(state[emt].p());
8248 double m2RadAft = radAfterBranch.m2Calc();
8249 double m2EmtAft = emtAfterBranch.m2Calc();
8250 double m2RadBef = 0.;
8251 if ( state[rad].idAbs() != 21 && state[rad].idAbs() != 22
8252 && state[emt].idAbs() != 24 && state[rad].idAbs() != state[emt].idAbs())
8253 m2RadBef = m2RadAft;
8254 else if ( state[emt].idAbs() == 24) {
8256 m2RadBef = pow2(particleDataPtr->m0(abs(idRadBef)));
8259 double Qsq = (radAfterBranch + emtAfterBranch).m2Calc();
8263 = (radAfterBranch + recAfterBranch + emtAfterBranch).m2Calc();
8265 if ( !state[rec].isFinal() ){
8266 double mar2 = m2final - 2. * Qsq + 2. * m2RadBef;
8267 recAfterBranch *= (1. - (Qsq - m2RadBef)/(mar2 - m2RadBef))
8268 /(1. + (Qsq - m2RadBef)/(mar2 - m2RadBef));
8271 if (Qsq > mar2)
return 0.5;
8274 Vec4 sum = radAfterBranch + recAfterBranch + emtAfterBranch;
8275 double m2Dip = sum.m2Calc();
8277 double x1 = 2. * (sum * radAfterBranch) / m2Dip;
8278 double x2 = 2. * (sum * recAfterBranch) / m2Dip;
8281 double lambda13 = sqrt( pow2(Qsq - m2RadAft - m2EmtAft )
8282 - 4.*m2RadAft*m2EmtAft);
8283 double k1 = ( Qsq - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8284 double k3 = ( Qsq - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8286 z = 1./ ( 1- k1 -k3) * ( x1 / (2.-x2) - k3);
8290 Vec4 qBR(state[rad].p() - state[emt].p() + state[rec].p());
8291 Vec4 qAR(state[rad].p() + state[rec].p());
8293 z = (qBR.m2Calc())/( qAR.m2Calc());
8304 double History::pTLund(
const Event& event,
int rad,
int emt,
int rec,
8305 int ShowerType,
int idRadBef) {
8308 Particle RadAfterBranch =
event[rad];
8309 Particle EmtAfterBranch =
event[emt];
8310 Particle RecAfterBranch =
event[rec];
8313 if ( mergingHooksPtr->useShowerPlugin() ) {
8314 map<string,double> stateVars;
8315 bool isFSR = showers->timesPtr->isTimelike(event, rad, emt, rec,
"");
8317 string name = showers->timesPtr->getSplittingName(event, rad, emt,
8319 stateVars = showers->timesPtr->getStateVariables(event, rad, emt, rec,
8322 string name = showers->spacePtr->getSplittingName(event, rad, emt,
8324 stateVars = showers->spacePtr->getStateVariables(event, rad, emt, rec,
8328 return ( (stateVars.size() > 0 && stateVars.find(
"t") != stateVars.end())
8329 ? sqrt(stateVars[
"t"]) : -1.0 );
8333 int Type = ShowerType;
8335 int sign = (Type==1) ? 1 : -1;
8336 Vec4 Q(RadAfterBranch.p() + sign*EmtAfterBranch.p());
8337 double Qsq = sign * Q.m2Calc();
8340 Vec4 radAft(RadAfterBranch.p());
8341 Vec4 recAft(RecAfterBranch.p());
8342 Vec4 emtAft(EmtAfterBranch.p());
8345 double m2RadAft = radAft.m2Calc();
8346 double m2EmtAft = emtAft.m2Calc();
8348 double m2RadBef = 0.;
8349 if ( RadAfterBranch.idAbs() != 21 && RadAfterBranch.idAbs() != 22
8350 && EmtAfterBranch.idAbs() != 24
8351 && RadAfterBranch.idAbs() != EmtAfterBranch.idAbs() )
8352 m2RadBef = m2RadAft;
8353 else if (EmtAfterBranch.idAbs() == 24) {
8355 m2RadBef = pow2(particleDataPtr->m0(abs(idRadBef)));
8356 }
else if (!RadAfterBranch.isFinal()) {
8357 if (RadAfterBranch.idAbs() == 21 && EmtAfterBranch.idAbs() != 21)
8358 m2RadBef = m2EmtAft;
8362 double m2final = (radAft + recAft + emtAft).m2Calc();
8364 if ( !RecAfterBranch.isFinal() && RadAfterBranch.isFinal() ){
8365 double mar2 = m2final - 2. * Qsq + 2. * m2RadBef;
8366 recAft *= (1. - (Qsq - m2RadBef)/(mar2 - m2RadBef))
8367 /(1. + (Qsq - m2RadBef)/(mar2 - m2RadBef));
8369 if (Qsq > mar2)
return 0.;
8372 Vec4 sum = radAft + recAft + emtAft;
8373 double m2Dip = sum.m2Calc();
8374 double x1 = 2. * (sum * radAft) / m2Dip;
8375 double x2 = 2. * (sum * recAft) / m2Dip;
8378 double q2BR = (RadAfterBranch.p() - EmtAfterBranch.p()
8379 + RecAfterBranch.p()).m2Calc();
8380 double q2AR = (RadAfterBranch.p() + RecAfterBranch.p()).m2Calc();
8383 double lambda13 = sqrt( pow2(Qsq - m2RadAft - m2EmtAft )
8384 - 4. * m2RadAft*m2EmtAft );
8385 double k1 = ( Qsq - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8386 double k3 = ( Qsq - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2. * Qsq );
8389 double z = (Type==1) ? 1./ ( 1- k1 -k3) * ( x1 / (2.-x2) - k3)
8393 double pTpyth = (Type==1) ? z*(1.-z) : (1.-z);
8396 if (Type == 1) pTpyth *= (Qsq - m2RadBef);
8402 if ( (RadAfterBranch.idAbs() == 4 || EmtAfterBranch.idAbs() == 4)
8403 && RadAfterBranch.idAbs() != EmtAfterBranch.idAbs()) {
8404 if (pTpyth < 2 * pow2(particleDataPtr->m0(4)))
8405 pTpyth = (Qsq + pow2(particleDataPtr->m0(4)) ) * (1. - q2BR/q2AR);
8406 }
else if ( (RadAfterBranch.idAbs() == 5 || EmtAfterBranch.idAbs() == 5)
8407 && RadAfterBranch.idAbs() != EmtAfterBranch.idAbs()) {
8408 if (pTpyth < 2 * pow2(particleDataPtr->m0(5)))
8409 pTpyth = (Qsq + pow2(particleDataPtr->m0(5)) ) * (1. - q2BR/q2AR);
8413 if ( pTpyth < 0. ) pTpyth = 0.;
8416 return sqrt(pTpyth);
8425 int History::posChangedIncoming(
const Event& event,
bool before) {
8431 for (
int i =0; i <
event.size(); ++i)
8432 if (event[i].status() == 43) {
8438 if (iSister > 0) iMother =
event[iSister].mother1();
8441 if (iSister > 0 && iMother > 0) {
8444 int flavSister =
event[iSister].id();
8445 int flavMother =
event[iMother].id();
8448 int flavDaughter = 0;
8449 if ( abs(flavMother) < 21 && flavSister == 21)
8450 flavDaughter = flavMother;
8451 else if ( flavMother == 21 && flavSister == 21)
8452 flavDaughter = flavMother;
8453 else if ( flavMother == 21 && abs(flavSister) < 21)
8454 flavDaughter = -1*flavSister;
8455 else if ( abs(flavMother) < 21 && abs(flavSister) < 21)
8460 for (
int i =0; i <
event.size(); ++i)
8461 if ( !event[i].isFinal()
8462 &&
event[i].mother1() == iMother
8463 &&
event[i].id() == flavDaughter )
8467 if ( !before )
return iMother;
8468 else return iDaughter;
8476 for (
int i =0; i <
event.size(); ++i)
8477 if ( abs(event[i].status()) == 53 || abs(event[i].status()) == 54) {
8483 if (iMother > 0) iDaughter =
event[iMother].daughter1();
8486 if (iDaughter > 0 && iMother > 0) {
8489 if ( !before )
return iMother;
8490 else return iDaughter;
8506 double History::pdfFactor(
const Event& event,
const int type,
8507 double pdfScale,
double mu ) {
8516 for (
int i =0; i <
event.size(); ++i)
8517 if ( abs(event[i].status()) == 53 || abs(event[i].status()) == 54) {
8521 int flavMother =
event[iMother].id();
8524 if ( iMother == 0 )
return 1.;
8527 int iDaughter =
event[iMother].daughter1();
8528 int flavDaughter =
event[iDaughter].id();
8531 double xMother = 2.*
event[iMother].e() /
event[0].e();
8532 double xDaughter = 2.*
event[iDaughter].e() /
event[0].e();
8536 int sideSplit = (
event[iMother].pz() > 0.) ? 1 : -1;
8537 double pdfDen1, pdfDen2, pdfNum1, pdfNum2;
8538 pdfDen1 = pdfDen2 = pdfNum1 = pdfNum2 = 1.;
8539 if ( sideSplit == 1 ) {
8541 pdfDen1 = max(1e-15,beamA.xfISR(0, flavDaughter, xDaughter, pow2(mu)) );
8542 pdfNum1 = beamA.xfISR(0, flavDaughter, xDaughter, pow2(pdfScale) );
8543 pdfNum2 = beamA.xfISR(0, flavMother, xMother, pow2(mu) );
8544 pdfDen2 = max(1e-15,beamA.xfISR(0,flavMother, xMother, pow2(pdfScale)) );
8547 pdfDen1 = max(1e-15,beamB.xfISR(0, flavDaughter, xDaughter, pow2(mu)) );
8548 pdfNum1 = beamB.xfISR(0, flavDaughter, xDaughter, pow2(pdfScale) );
8549 pdfNum2 = beamB.xfISR(0, flavMother, xMother, pow2(mu) );
8550 pdfDen2 = max(1e-15,beamB.xfISR(0,flavMother, xMother, pow2(pdfScale)) );
8555 if ( pdfDen2/pdfNum1 > 1. )
return 1.;
8560 weight = (pdfNum1/pdfDen1) * (pdfNum2)/(pdfDen2);
8563 }
else if (type == 2) {
8567 for (
int i =0; i <
event.size(); ++i)
8568 if (event[i].status() == 43) {
8572 int flavSister =
event[iSister].id();
8575 int iMother =
event[iSister].mother1();
8576 int flavMother =
event[iMother].id();
8579 int flavDaughter = 0;
8580 if ( abs(flavMother) < 21 && flavSister == 21)
8581 flavDaughter = flavMother;
8582 else if ( flavMother == 21 && flavSister == 21)
8583 flavDaughter = flavMother;
8584 else if ( flavMother == 21 && abs(flavSister) < 21)
8585 flavDaughter = -1*flavSister;
8586 else if ( abs(flavMother) < 21 && abs(flavSister) < 21)
8590 double xMother = 2.*
event[iMother].e() /
event[0].e();
8594 for (
int i =0; i <
event.size(); ++i)
8595 if ( !event[i].isFinal()
8596 &&
event[i].mother1() == iMother
8597 &&
event[i].id() == flavDaughter )
8599 double xDaughter = 2.*
event[iDaughter].e() /
event[0].e();
8604 int sideSplit = (
event[iMother].pz() > 0.) ? 1 : -1;
8605 double ratio1 = getPDFratio( sideSplit,
false,
false, flavDaughter,
8606 xDaughter, pdfScale, flavDaughter, xDaughter, mu );
8607 double ratio2 = getPDFratio( sideSplit,
false,
false, flavMother,
8608 xMother, mu, flavMother, xMother, pdfScale );
8610 weight = ratio1*ratio2;
8627 double History::integrand(
int flav,
double x,
double scaleInt,
double z) {
8639 AlphaStrong* as = mergingHooksPtr->AlphaS_ISR();
8640 double asNow = (*as).alphaS(z);
8641 result = 1./z *asNow*asNow* ( log(scaleInt/z) -3./2. );
8645 }
else if (flav==21) {
8647 double measure1 = 1./(1. - z);
8648 double measure2 = 1.;
8652 * z * beamB.xf( 21,x/z,pow(scaleInt,2))
8653 / beamB.xf( 21,x, pow(scaleInt,2))
8658 2.*CA *((1. -z)/z + z*(1.-z))
8659 * beamB.xf( 21,x/z,pow(scaleInt,2))
8660 / beamB.xf( 21,x, pow(scaleInt,2))
8662 + CF * ((1+pow(1-z,2))/z)
8663 *( beamB.xf( 1, x/z,pow(scaleInt,2))
8664 / beamB.xf( 21, x, pow(scaleInt,2))
8665 + beamB.xf( -1, x/z,pow(scaleInt,2))
8666 / beamB.xf( 21, x, pow(scaleInt,2))
8667 + beamB.xf( 2, x/z,pow(scaleInt,2))
8668 / beamB.xf( 21, x, pow(scaleInt,2))
8669 + beamB.xf( -2, x/z,pow(scaleInt,2))
8670 / beamB.xf( 21, x, pow(scaleInt,2))
8671 + beamB.xf( 3, x/z,pow(scaleInt,2))
8672 / beamB.xf( 21, x, pow(scaleInt,2))
8673 + beamB.xf( -3, x/z,pow(scaleInt,2))
8674 / beamB.xf( 21, x, pow(scaleInt,2))
8675 + beamB.xf( 4, x/z,pow(scaleInt,2))
8676 / beamB.xf( 21, x, pow(scaleInt,2))
8677 + beamB.xf( -4, x/z,pow(scaleInt,2))
8678 / beamB.xf( 21, x, pow(scaleInt,2)) );
8681 result = integrand1*measure1 + integrand2*measure2;
8685 double measure1 = 1./(1. -z);
8686 double measure2 = 1.;
8691 * beamB.xf( flav, x/z, pow(scaleInt,2))
8692 / beamB.xf( flav, x, pow(scaleInt,2))
8697 + TR * (pow(z,2) + pow(1-z,2))
8698 * beamB.xf( 21, x/z, pow(scaleInt,2))
8699 / beamB.xf( flav, x, pow(scaleInt,2));
8702 result = measure1*integrand1 + measure2*integrand2;
8714 vector<int> History::posFlavCKM(
int flav) {
8717 int flavAbs = abs(flav);
8719 vector<int> flavRadBefs;
8721 if (flavAbs > 10 && flavAbs % 2 == 1)
8722 flavRadBefs.push_back(flavAbs + 1);
8724 else if (flavAbs > 10 && flavAbs % 2 == 0)
8725 flavRadBefs.push_back(flavAbs - 1);
8727 else if (flavAbs < 10 && flavAbs % 2 == 1) {
8728 flavRadBefs.push_back(2);
8729 flavRadBefs.push_back(4);
8730 flavRadBefs.push_back(6);
8732 else if (flavAbs < 10 && flavAbs % 2 == 0) {
8733 flavRadBefs.push_back(1);
8734 flavRadBefs.push_back(3);
8735 flavRadBefs.push_back(5);
8748 bool History::checkFlavour(vector<int>& flavCounts,
int flavRad,
8749 int flavRadBef,
int clusType) {
8752 for(
int k = 0; k < 20; ++k) {
8755 if (abs(flavRad) == k) {
8757 if (flavRad < 0) cor = 1;
8760 if (abs(flavRadBef) == k) {
8762 if (flavRadBef < 0) cor = -1;
8766 if (flavRadBef == flavRad) cor = 0;
8769 if (clusType == 1) {
8770 if (flavCounts[k] + cor != 0)
return false;
8773 if (flavCounts[k] - cor != 0)
return false;
8786 void History::reverseBoostISR(Vec4& pMother, Vec4& pSister, Vec4& pPartner,
8787 Vec4& pDaughter, Vec4& pRecoiler,
int sign,
double eCM,
double& phi ) {
8791 phi = pSister.phi();
8793 RotBstMatrix rot_by_mphi;
8794 rot_by_mphi.rot(0.,-phi);
8796 RotBstMatrix rot_by_pphi;
8797 rot_by_pphi.rot(0.,phi);
8801 double x1 = 2. * pMother.e() / eCM;
8803 double x2 = 2. * pPartner.e() / eCM;
8806 Vec4 qDip( pMother - pSister);
8807 Vec4 qAfter(pMother + pPartner);
8808 Vec4 qBefore(qDip + pPartner);
8809 double z = qBefore.m2Calc() / qAfter.m2Calc();
8812 double x1New = z*x1;
8814 double sHat = x1New*x2New*eCM*eCM;
8819 Vec4 pDaughterBef( 0., 0., sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
8820 Vec4 pRecoilerBef( 0., 0., -sign*0.5*sqrt(sHat), 0.5*sqrt(sHat));
8823 pMother.rotbst( rot_by_mphi );
8824 pSister.rotbst( rot_by_mphi );
8825 pPartner.rotbst( rot_by_mphi );
8829 pDaughter.p( pMother - pSister);
8830 pRecoiler.p( pPartner );
8831 RotBstMatrix from_CM_to_DRoff;
8833 from_CM_to_DRoff.toCMframe(pDaughter, pRecoiler);
8835 from_CM_to_DRoff.toCMframe(pRecoiler, pDaughter);
8839 pMother.rotbst( from_CM_to_DRoff );
8840 pPartner.rotbst( from_CM_to_DRoff );
8841 pSister.rotbst( from_CM_to_DRoff );
8846 RotBstMatrix from_DR_to_CM;
8847 from_DR_to_CM.bst( 0., 0., sign*( x1New - x2New ) / ( x1New + x2New ) );
8852 pDaughterBef.rotbst( from_DR_to_CM );
8853 pRecoilerBef.rotbst( from_DR_to_CM );
8857 if ( abs(pRecoilerBef.mCalc()) > 1e-7 ) {
8858 double pzSign = (pRecoilerBef.pz() > 0.) ? 1. : -1.;
8859 double eRec = pRecoilerBef.e();
8860 pRecoilerBef.p(0., 0., pzSign*eRec, eRec);
8862 if ( abs(pDaughterBef.mCalc()) > 1e-7 ) {
8863 double pzSign = (pDaughterBef.pz() > 0.) ? 1. : -1.;
8864 double eDau = pDaughterBef.e();
8865 pDaughterBef.p(0., 0., pzSign*eDau, eDau);
8877 bool History::isQCD2to2(
const Event& event) {
8879 if (!mergingHooksPtr->doWeakClustering())
return false;
8882 int nFinalPartons = 0, nFinal = 0;;
8883 for (
int i = 0;i <
event.size();++i)
8884 if (event[i].isFinal()) {
8886 if ( event[i].idAbs() < 10 ||
event[i].idAbs() == 21)
8889 if (nFinalPartons == 2 && nFinal == 2)
return true;
8899 bool History::isEW2to1(
const Event& event) {
8901 if (!mergingHooksPtr->doWeakClustering())
return false;
8904 for (
int i = 0;i <
event.size();++i) {
8905 if (event[i].isFinal()) {
8906 if (event[i].idAbs() == 23 ||
8907 event[i].idAbs() == 24 ||
8908 event[i].idAbs() == 22) nVector++;
8914 if (nVector == 1)
return true;
8924 void History::setSelectedChild() {
8925 if (mother == 0)
return;
8926 for (
int i = 0;i < int(mother->children.size());++i)
8927 if (mother->children[i] ==
this) mother->selectedChild = i;
8928 mother->setSelectedChild();
8933 void History::setupWeakShower(
int nSteps) {
8936 if (selectedChild != -1) {
8937 children[selectedChild]->setupWeakShower(nSteps+1);
8942 vector<int> mode, fermionLines;
8944 vector<pair<int,int> > dipoles;
8947 setupWeakHard(mode,fermionLines, mom);
8950 if (isQCD2to2(state)) {
8952 if (state[3].idAbs() < 10) dipoles.push_back(make_pair(3,4));
8953 if (state[4].idAbs() < 10) dipoles.push_back(make_pair(4,3));
8954 if (state[5].idAbs() < 10) dipoles.push_back(make_pair(5,6));
8955 if (state[6].idAbs() < 10) dipoles.push_back(make_pair(6,5));
8956 }
else if (isEW2to1(state)) {
8957 if (state[3].idAbs() < 10) dipoles.push_back(make_pair(3,4));
8958 if (state[4].idAbs() < 10) dipoles.push_back(make_pair(4,3));
8962 transferWeakShower(mode, mom, fermionLines, dipoles, nSteps);
8968 void History::transferWeakShower(vector<int> &mode, vector<Vec4> &mom,
8969 vector<int> fermionLines, vector<pair<int,int> > &dipoles,
8974 infoPtr->setWeakModes(mode);
8975 infoPtr->setWeakDipoles(dipoles);
8976 infoPtr->setWeakMomenta(mom);
8977 infoPtr->setWeak2to2lines(fermionLines);
8982 map<int,int> stateTransfer;
8983 findStateTransfer(stateTransfer);
8986 vector<int> modeNew = updateWeakModes(mode, stateTransfer);
8987 vector<int> fermionLinesNew = updateWeakFermionLines(fermionLines,
8989 vector<pair<int, int> > dipolesNew = updateWeakDipoles(dipoles,
8993 mother->transferWeakShower(modeNew, mom, fermionLinesNew, dipolesNew,
9000 vector<int> History::updateWeakModes(vector<int>& mode,
9001 map<int,int>& stateTransfer) {
9003 vector<int> modeNew(mode.size() + 1,0);
9006 for (map<int,int>::iterator it = stateTransfer.begin();
9007 it != stateTransfer.end(); ++it)
9008 modeNew[it->second] = mode[it->first];
9010 modeNew[clusterIn.emitted] = mode[clusterIn.radBef];
9014 if (state[clusterIn.radBef].idAbs() == 21 &&
9015 mother->state[clusterIn.emittor].idAbs() != 21) {
9017 if (state[clusterIn.radBef].status() > 0) modeNew[clusterIn.emittor] = 1;
9020 if (modeNew[clusterIn.emittor] == 1);
9021 else if ( mother->state[clusterIn.recoiler].id() == 21)
9022 modeNew[clusterIn.emittor] = 2;
9023 else if ( mother->state[clusterIn.recoiler].id()
9024 == mother->state[clusterIn.emittor].id() )
9025 modeNew[clusterIn.emittor] = 4;
9026 else modeNew[clusterIn.emittor] = 3;
9029 modeNew[clusterIn.emitted] = 1;
9033 if (state[clusterIn.radBef].idAbs() < 10 &&
9034 mother->state[clusterIn.emittor].idAbs() == 21) {
9035 if (state[clusterIn.radBef].status() < 0) {
9036 modeNew[clusterIn.emitted] = 1;
9041 if (state[clusterIn.radBef].idAbs() == 22) {
9043 if (state[clusterIn.radBef].status() > 0) modeNew[clusterIn.emittor] = 1;
9046 if (modeNew[clusterIn.emittor] == 1);
9047 else if ( mother->state[clusterIn.recoiler].id() == 21)
9048 modeNew[clusterIn.emittor] = 2;
9049 else if ( mother->state[clusterIn.recoiler].id()
9050 == mother->state[clusterIn.emittor].id() )
9051 modeNew[clusterIn.emittor] = 4;
9052 else modeNew[clusterIn.emittor] = 3;
9055 modeNew[clusterIn.emitted] = 1;
9064 vector<int> History::updateWeakFermionLines(vector<int> fermionLines,
9065 map<int,int>& stateTransfer) {
9068 if (!fermionLines.empty()) {
9070 fermionLines[0] = stateTransfer[fermionLines[0]];
9071 fermionLines[1] = stateTransfer[fermionLines[1]];
9074 bool lines[2] = {
false,
false};
9075 if (fermionLines[2] != clusterIn.radBef)
9076 fermionLines[2] = stateTransfer[fermionLines[2]];
9077 else lines[0] =
true;
9078 if (fermionLines[3] != clusterIn.radBef)
9079 fermionLines[3] = stateTransfer[fermionLines[3]];
9080 else lines[1] =
true;
9083 for (
int i = 0;i < 2; ++i) {
9085 if (state[fermionLines[2 + i]].isQuark() ||
9086 state[fermionLines[2 + i]].isLepton()) {
9087 if (mother->state[clusterIn.emittor].isQuark() ||
9088 mother->state[clusterIn.emittor].isLepton())
9089 fermionLines[2 + i] = clusterIn.emittor;
9090 else fermionLines[2 + i] = clusterIn.emitted;
9093 else fermionLines[2 + i] = 0;
9097 return fermionLines;
9103 vector<pair<int,int> > History::updateWeakDipoles(
9104 vector<pair<int,int> > &dipoles, map<int,int>& stateTransfer) {
9106 vector<pair<int,int> > dipolesNew;
9107 for (
int i = 0;i < int(dipoles.size());++i) {
9108 int iRecNew = -1, iRadNew = -1;
9111 if (dipoles[i].first != clusterIn.radBef)
9112 iRadNew = stateTransfer[dipoles[i].first];
9114 else if (state[clusterIn.radBef].status() > 0) {
9115 if (mother->state[clusterIn.emitted].id() ==
9116 state[clusterIn.radBef].id())
9117 iRadNew = clusterIn.emitted;
9118 else iRadNew = clusterIn.emittor;
9120 }
else if (mother->state[clusterIn.emittor].idAbs() < 10)
9121 iRadNew = clusterIn.emittor;
9124 if (iRadNew == -1)
continue;
9127 if (dipoles[i].second != clusterIn.radBef)
9128 iRecNew = stateTransfer[dipoles[i].second];
9130 else if (state[clusterIn.radBef].status() > 0) {
9132 if (mother->state[clusterIn.emitted].id() == 21 &&
9133 mother->state[clusterIn.emittor].id() == 21) {
9134 double m1 = (mother->state[clusterIn.emitted].p()
9135 + mother->state[iRadNew].p()).m2Calc();
9136 double m2 = (mother->state[clusterIn.emittor].p()
9137 + mother->state[iRadNew].p()).m2Calc();
9138 iRecNew = (m1 > m2) ? clusterIn.emitted : clusterIn.emittor;
9141 else if (mother->state[clusterIn.emitted].id() ==
9142 state[clusterIn.radBef].id())
9143 iRecNew = clusterIn.emitted;
9144 else iRecNew = clusterIn.emittor;
9146 }
else iRecNew = clusterIn.emittor;
9148 dipolesNew.push_back(make_pair(iRadNew,iRecNew));
9152 if (state[clusterIn.radBef].idAbs() == 21 &&
9153 mother->state[clusterIn.emittor].idAbs() != 21) {
9155 if (state[clusterIn.radBef].status() > 0) {
9156 dipolesNew.push_back(make_pair(clusterIn.emittor,clusterIn.emitted));
9157 dipolesNew.push_back(make_pair(clusterIn.emitted,clusterIn.emittor));
9160 int iRad = clusterIn.emittor;
9161 int iRec = (iRad == 3) ? 4 : 3;
9162 dipolesNew.push_back(make_pair(iRad,iRec));
9163 dipolesNew.push_back(make_pair(clusterIn.emitted,findISRRecoiler()));
9168 if (state[clusterIn.radBef].idAbs() < 10 &&
9169 mother->state[clusterIn.emittor].idAbs() == 21 &&
9170 state[clusterIn.radBef].status() < 0)
9171 dipolesNew.push_back(make_pair(clusterIn.emitted,findISRRecoiler()));
9180 void History::setupWeakHard(vector<int>& mode, vector<int>& fermionLines,
9181 vector<Vec4>& mom) {
9183 if (!isQCD2to2(state)) {
9185 mode.resize(state.size(), 1);
9189 for (
int i = 0;i < 4; ++i) {
9190 mom.push_back(state[3 + i].p());
9191 fermionLines.push_back(3 + i);
9194 if ( state[3].idAbs() == 21 && state[4].idAbs() == 21 &&
9195 state[5].idAbs() == 21 && state[6].idAbs() == 21)
9196 mode.resize(state.size(), 1);
9199 else if (state[5].
id() == -state[6].id() ||
9200 (state[5].idAbs() == 21 && state[6].idAbs() == 21))
9201 mode.resize(state.size(), 1);
9204 else if (state[5].idAbs() == 21 || state[6].idAbs() == 21) {
9205 mode.resize(state.size(), 2);
9206 if (state[3].
id() != state[5].id()) {
9207 swap(mom[0], mom[1]);
9208 swap(mom[2], mom[3]);
9213 else if (state[5].
id() != state[6].
id()) {
9214 mode.resize(state.size(), 3);
9215 if (state[3].
id() != state[5].id()) {
9216 swap(mom[0], mom[1]);
9217 swap(mom[2], mom[3]);
9223 else if (state[5].
id() == state[6].
id()) {
9224 mode.resize(state.size(), 4);
9233 double History::getShowerPluginScale(
const Event& event,
int rad,
int emt,
9234 int rec,
string key,
double scalePythia) {
9237 if ( !mergingHooksPtr->useShowerPlugin() )
return scalePythia;
9240 map<string,double> stateVars;
9241 bool isFSR = showers->timesPtr->isTimelike(event, rad, emt, rec,
"");
9243 string name = showers->timesPtr->getSplittingName(event, rad, emt,
9245 stateVars = showers->timesPtr->getStateVariables(event, rad, emt, rec,
9248 string name = showers->spacePtr->getSplittingName(event, rad, emt,
9250 stateVars = showers->spacePtr->getStateVariables(event, rad, emt, rec,
9254 return ( (stateVars.size() > 0 && stateVars.find(key) != stateVars.end())
9255 ? stateVars[key] : -1.0 );