8 #include "Pythia8/DireSpace.h"
9 #include "Pythia8/DireTimes.h"
10 #include "Pythia8/DireHistory.h"
25 const int DireSpace::MAXLOOPTINYPDF = 10;
28 const double DireSpace::MCMIN = 1.2;
29 const double DireSpace::MBMIN = 4.0;
33 const double DireSpace::CTHRESHOLD = 2.0;
34 const double DireSpace::BTHRESHOLD = 2.0;
38 const double DireSpace::EVALPDFSTEP = 0.1;
42 const double DireSpace::TINYPDF = 1e-5;
45 const double DireSpace::TINYKERNELPDF = 1e-15;
48 const double DireSpace::TINYPT2 = 0.25e-6;
52 const double DireSpace::HEAVYPT2EVOL = 1.1;
57 const double DireSpace::HEAVYXEVOL = 0.9;
62 const double DireSpace::EXTRASPACEQ = 2.0;
65 const double DireSpace::LAMBDA3MARGIN = 1.1;
68 const double DireSpace::PT2MINWARN = 1.;
72 const double DireSpace::LEPTONXMIN = 1e-10;
73 const double DireSpace::LEPTONXMAX = 1. - 1e-10;
76 const double DireSpace::LEPTONPT2MIN = 1.2;
79 const double DireSpace::LEPTONFUDGE = 10.;
82 const double DireSpace::HEADROOMQ2G = 1.35;
83 const double DireSpace::HEADROOMQ2Q = 1.15;
84 const double DireSpace::HEADROOMG2Q = 1.35;
85 const double DireSpace::HEADROOMG2G = 1.35;
87 const double DireSpace::TINYMASS = 1e-3;
89 const double DireSpace::DPHI_II = 2./3.;
90 const double DireSpace::DPHI_IF = 2./3.;
93 const double DireSpace::PT2_INCREASE_OVERESTIMATE = 2.;
95 const double DireSpace::KERNEL_HEADROOM = 1.;
101 void DireSpace::init( BeamParticle* beamAPtrIn,
102 BeamParticle* beamBPtrIn) {
107 CA = settingsPtr->parm(
"DireColorQCD:CA") > 0.0
108 ? settingsPtr->parm(
"DireColorQCD:CA") : 3.0;
109 CF = settingsPtr->parm(
"DireColorQCD:CF") > 0.0
110 ? settingsPtr->parm(
"DireColorQCD:CF") : 4./3.;
111 TR = settingsPtr->parm(
"DireColorQCD:TR") > 0.
112 ? settingsPtr->parm(
"DireColorQCD:TR") : 0.5;
113 NC = settingsPtr->parm(
"DireColorQCD:NC") > 0.
114 ? settingsPtr->parm(
"DireColorQCD:NC") : 3.0;
117 beamAPtr = beamAPtrIn;
118 beamBPtr = beamBPtrIn;
121 doQCDshower = settingsPtr->flag(
"SpaceShower:QCDshower");
122 doQEDshowerByQ = settingsPtr->flag(
"SpaceShower:QEDshowerByQ");
123 doQEDshowerByL = settingsPtr->flag(
"SpaceShower:QEDshowerByL");
124 doDecaysAsShower = settingsPtr->flag(
"DireSpace:DecaysAsShower");
127 pTmaxMatch = settingsPtr->mode(
"SpaceShower:pTmaxMatch");
128 pTdampMatch = settingsPtr->mode(
"SpaceShower:pTdampMatch");
129 pTmaxFudge = settingsPtr->parm(
"SpaceShower:pTmaxFudge");
130 pTmaxFudgeMPI = settingsPtr->parm(
"SpaceShower:pTmaxFudgeMPI");
131 pTdampFudge = settingsPtr->parm(
"SpaceShower:pTdampFudge");
132 pT2minVariations= pow2(max(0.,settingsPtr->parm(
"Variations:pTmin")));
133 pT2minEnhance = pow2(max(0.,settingsPtr->parm(
"Enhance:pTmin")));
134 pT2minMECs = pow2(max(0.,settingsPtr->parm(
"Dire:pTminMECs")));
135 Q2minMECs = pow2(max(0.,settingsPtr->parm(
"Dire:QminMECs")));
136 nFinalMaxMECs = settingsPtr->mode(
"Dire:nFinalMaxMECs");
137 suppressLargeMECs = settingsPtr->flag(
"Dire:suppressLargeMECs");
140 doRapidityOrder = settingsPtr->flag(
"SpaceShower:rapidityOrder");
143 mc = max( MCMIN, particleDataPtr->m0(4));
144 mb = max( MBMIN, particleDataPtr->m0(5));
149 renormMultFac = settingsPtr->parm(
"SpaceShower:renormMultFac");
150 factorMultFac = settingsPtr->parm(
"SpaceShower:factorMultFac");
151 useFixedFacScale = settingsPtr->flag(
"SpaceShower:useFixedFacScale");
152 fixedFacScale2 = pow2(settingsPtr->parm(
"SpaceShower:fixedFacScale"));
155 alphaSvalue = settingsPtr->parm(
"SpaceShower:alphaSvalue");
156 alphaSorder = settingsPtr->mode(
"SpaceShower:alphaSorder");
157 alphaSnfmax = settingsPtr->mode(
"StandardModel:alphaSnfmax");
158 alphaSuseCMW = settingsPtr->flag(
"SpaceShower:alphaSuseCMW");
159 alphaS2pi = 0.5 * alphaSvalue / M_PI;
160 asScheme = settingsPtr->mode(
"DireSpace:alphasScheme");
163 double mcpy = particleDataPtr->m0(4);
164 double mbpy = particleDataPtr->m0(5);
165 double mtpy = particleDataPtr->m0(6);
166 if (mcpy > 0.0 && mbpy > 0.0 && mtpy > 0.0)
167 alphaS.setThresholds(mcpy, mbpy, mtpy);
170 alphaS.init( alphaSvalue, alphaSorder, alphaSnfmax, alphaSuseCMW);
173 Lambda5flav = alphaS.Lambda5();
174 Lambda4flav = alphaS.Lambda4();
175 Lambda3flav = alphaS.Lambda3();
176 Lambda5flav2 = pow2(Lambda5flav);
177 Lambda4flav2 = pow2(Lambda4flav);
178 Lambda3flav2 = pow2(Lambda3flav);
182 useSamePTasMPI = settingsPtr->flag(
"SpaceShower:samePTasMPI");
183 if (useSamePTasMPI) {
184 pT0Ref = settingsPtr->parm(
"MultipartonInteractions:pT0Ref");
185 ecmRef = settingsPtr->parm(
"MultipartonInteractions:ecmRef");
186 ecmPow = settingsPtr->parm(
"MultipartonInteractions:ecmPow");
187 pTmin = settingsPtr->parm(
"MultipartonInteractions:pTmin");
189 pT0Ref = settingsPtr->parm(
"SpaceShower:pT0Ref");
190 ecmRef = settingsPtr->parm(
"SpaceShower:ecmRef");
191 ecmPow = settingsPtr->parm(
"SpaceShower:ecmPow");
192 pTmin = settingsPtr->parm(
"SpaceShower:pTmin");
196 sCM = m2( beamAPtr->p(), beamBPtr->p());
201 double pTminAbs = sqrtpos(pow2(LAMBDA3MARGIN) * Lambda3flav2 / renormMultFac
203 if (pTmin < pTminAbs) {
205 ostringstream newPTmin;
206 newPTmin << fixed << setprecision(3) << pTmin;
207 infoPtr->errorMsg(
"Warning in DireSpace::init: pTmin too low",
208 ", raised to " + newPTmin.str() );
209 infoPtr->setTooLowPTmin(
true);
214 pT2min = pow2(pTmin);
217 mTolErr = settingsPtr->parm(
"Check:mTolErr");
219 double pT2minQED = pow2(settingsPtr->parm(
"SpaceShower:pTminChgQ"));
220 pT2minQED = max(pT2minQED, pow2(settingsPtr->parm(
"SpaceShower:pTminChgL")));
222 pT2cutSave = create_unordered_map<int,double>
224 (1,pT2min)(-1,pT2min)(2,pT2min)(-2,pT2min)
225 (3,pT2min)(-3,pT2min)(4,pT2min)(-4,pT2min)
226 (5,pT2min)(-5,pT2min)(6,pT2min)(-6,pT2min)
228 (11,pT2minQED)(-11,pT2minQED)(13,pT2minQED)(-13,pT2minQED)
229 (15,pT2minQED)(-15,pT2minQED)
230 (900032,pT2minQED)(900012,pT2minQED)
233 bool_settings = create_unordered_map<string,bool>
234 (
"doQEDshowerByL",doQEDshowerByL)
235 (
"doQEDshowerByQ",doQEDshowerByQ);
237 usePDFalphas = settingsPtr->flag(
"ShowerPDF:usePDFalphas");
238 useSummedPDF = settingsPtr->flag(
"ShowerPDF:useSummedPDF");
239 usePDF = settingsPtr->flag(
"ShowerPDF:usePDF");
240 BeamParticle&
beam = (particleDataPtr->isHadron(beamAPtr->id())) ?
241 *beamAPtr : *beamBPtr;
242 alphaS2piOverestimate = (usePDFalphas) ? beam.alphaS(pT2min) * 0.5/M_PI
243 : (alphaSorder > 0) ? alphaS.alphaS(pT2min) * 0.5/M_PI
245 usePDFmasses = settingsPtr->flag(
"ShowerPDF:usePDFmasses");
246 BeamParticle* bb = ( particleDataPtr->isHadron(beamAPtr->id())) ? beamAPtr
247 : ( particleDataPtr->isHadron(beamBPtr->id())) ?
249 m2cPhys = (usePDFalphas && bb != NULL)
250 ? pow2(max(0.,bb->mQuarkPDF(4))) : alphaS.muThres2(4);
251 m2bPhys = (usePDFalphas && bb != NULL)
252 ? pow2(max(0.,bb->mQuarkPDF(5))) : alphaS.muThres2(5);
256 useMassiveBeams =
false;
260 rejectProbability.insert( make_pair(key, multimap<double,double>() ));
261 acceptProbability.insert( make_pair(key, map<double,double>() ));
262 doVariations = settingsPtr->flag(
"Variations:doVariations");
267 if (splittingsPtr) splits = splittingsPtr->getSplittings();
269 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
270 it != splits.end(); ++it ) overhead.insert(make_pair(it->first,1.));
272 nFinalMax = settingsPtr->mode(
"DireSpace:nFinalMax");
273 useGlobalMapIF = settingsPtr->flag(
"DireSpace:useGlobalMapIF");
275 forceMassiveMap = settingsPtr->flag(
"DireSpace:forceMassiveMap");
278 kernelOrder = settingsPtr->mode(
"DireSpace:kernelOrder");
279 kernelOrderMPI = settingsPtr->mode(
"DireSpace:kernelOrderMPI");
282 doMEcorrections = settingsPtr->flag(
"Dire:doMECs")
283 || settingsPtr->flag(
"Dire:doMOPS");
284 doMEafterFirst = settingsPtr->flag(
"SpaceShower:MEafterFirst");
285 doPhiPolAsym = settingsPtr->flag(
"SpaceShower:phiPolAsym");
286 doPhiIntAsym = settingsPtr->flag(
"SpaceShower:phiIntAsym");
287 strengthIntAsym = settingsPtr->parm(
"SpaceShower:strengthIntAsym");
288 nQuarkIn = settingsPtr->mode(
"SpaceShower:nQuarkIn");
291 doSecondHard = settingsPtr->flag(
"SecondHard:generate");
295 = settingsPtr->mode(
"MultipartonInteractions:enhanceScreening");
296 if (!useSamePTasMPI) enhanceScreening = 0;
299 hasUserHooks = (userHooksPtr != 0);
300 canVetoEmission = (userHooksPtr != 0)
301 ? userHooksPtr->canVetoISREmission() :
false;
312 void DireSpace::initVariations() {
315 for (
int i=0; i < weights->sizeWeights(); ++i) {
316 string key = weights->weightName(i);
317 if ( key.compare(
"base") == 0)
continue;
318 if ( key.find(
"fsr") != string::npos)
continue;
319 rejectProbability.insert( make_pair(key, multimap<double,double>() ));
320 acceptProbability.insert( make_pair(key, map<double,double>() ));
323 for ( unordered_map<
string, multimap<double,double> >::iterator
324 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
326 for ( unordered_map<
string, map<double,double> >::iterator
327 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
339 bool DireSpace::limitPTmax(
Event& event,
double,
double) {
342 bool dopTlimit =
false;
343 dopTlimit1 = dopTlimit2 =
false;
345 if (pTmaxMatch == 1) dopTlimit = dopTlimit1 = dopTlimit2 =
true;
346 else if (pTmaxMatch == 2) dopTlimit = dopTlimit1 = dopTlimit2 =
false;
349 else if (infoPtr->isNonDiffractive() || infoPtr->isDiffractiveA()
350 || infoPtr->isDiffractiveB() || infoPtr->isDiffractiveC() )
351 dopTlimit = dopTlimit1 = dopTlimit2 =
true;
357 for (
int i = 5; i <
event.size(); ++i) {
358 if (event[i].status() == -21) ++n21;
360 int idAbs =
event[i].idAbs();
361 if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit1 =
true;
362 if ( (event[i].col() != 0 || event[i].acol() != 0)
363 && idAbs > 5 && idAbs != 21 ) ++nHeavyCol;
364 }
else if (n21 == 2) {
365 int idAbs =
event[i].idAbs();
366 if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit2 =
true;
369 dopTlimit = (doSecondHard) ? (dopTlimit1 && dopTlimit2) : dopTlimit1;
385 void DireSpace::resetWeights() {
389 for ( unordered_map<
string, multimap<double,double> >::iterator
390 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
392 for ( unordered_map<
string, map<double,double> >::iterator
393 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
404 void DireSpace::prepare(
int iSys,
Event& event,
bool limitPTmaxIn) {
407 if (nMPI < infoPtr->getCounter(23) && iSys == infoPtr->getCounter(23)) {
408 weights->calcWeight(pow2(infoPtr->pTnow()));
411 for ( unordered_map<
string, multimap<double,double> >::iterator
412 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
414 for ( unordered_map<
string, map<double,double> >::iterator
415 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
420 nMPI = infoPtr->getCounter(23);
423 int in1 = getInA(iSys);
424 int in2 = getInB(iSys);
427 bool canRadiate1 = !(
event[in1].isRescatteredIncoming());
428 bool canRadiate2 = !(
event[in2].isRescatteredIncoming());
431 if (iSys == 0) {dipEnd.resize(0); idResFirst = 0;}
432 else if (iSys == 1) idResSecond = 0;
435 splits = splittingsPtr->getSplittings();
437 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
438 it != splits.end(); ++it ) overhead.insert(make_pair(it->first,1.));
444 if (doSecondHard && iSys == 0) limitPTmaxIn = dopTlimit1;
445 if (doSecondHard && iSys == 1) limitPTmaxIn = dopTlimit2;
452 int colTag =
event[in1].col();
453 if (canRadiate1 && colTag > 0) setupQCDdip( iSys, 1, colTag, 1, event,
454 MEtype, limitPTmaxIn);
456 int acolTag =
event[in1].acol();
457 if (canRadiate1 && acolTag > 0) setupQCDdip( iSys, 1, acolTag, -1, event,
458 MEtype, limitPTmaxIn);
460 colTag =
event[in2].col();
461 if (canRadiate2 && colTag > 0) setupQCDdip( iSys, 2, colTag, 1, event,
462 MEtype, limitPTmaxIn);
464 acolTag =
event[in2].acol();
465 if (canRadiate2 && acolTag > 0) setupQCDdip( iSys, 2, acolTag, -1, event,
466 MEtype, limitPTmaxIn);
470 getGenDip( iSys, 1, event, limitPTmaxIn, dipEnd);
471 getGenDip( iSys, 2, event, limitPTmaxIn, dipEnd);
475 if (iSys == 0 && infoPtr->hasHistory()) {
476 double zNow = infoPtr->zNowISR();
477 double pT2Now = infoPtr->pT2NowISR();
478 for (
int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
479 dipEnd[iDipEnd].zOld = zNow;
480 dipEnd[iDipEnd].pT2Old = pT2Now;
481 ++dipEnd[iDipEnd].nBranch;
486 updateDipoles(event, iSys);
490 if ( nProposedPT.find(iSys) == nProposedPT.end() )
491 nProposedPT.insert(make_pair(iSys,0));
498 for ( unordered_map<
string, multimap<double,double> >::iterator
499 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
501 for ( unordered_map<
string, map<double,double> >::iterator
502 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
509 void DireSpace::clear() {
519 for ( unordered_map<
string, multimap<double,double> >::iterator
520 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
522 for ( unordered_map<
string, map<double,double> >::iterator
523 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
532 void DireSpace::setupQCDdip(
int iSys,
int side,
int colTag,
int colSign,
533 const Event& event,
int MEtype,
bool limitPTmaxIn) {
536 int iRad = (side == 1) ? getInA(iSys) : getInB(iSys);
538 int sizeAllA = partonSystemsPtr->sizeAll(iSys);
539 int sizeOut = partonSystemsPtr->sizeOut(iSys);
540 int sizeAll = sizeAllA;
541 int sizeIn = sizeAll - sizeOut;
542 int sizeInA = sizeAllA - sizeIn - sizeOut;
546 for (
int j = 0; j < sizeAll; ++j) {
547 int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA);
548 if (iRecNow == iRad)
continue;
549 if ( ( j >= sizeIn && event[iRecNow].col() == colTag
550 && event[iRecNow].isFinal() )
551 || ( j < sizeIn && event[iRecNow].acol() == colTag
552 && !event[iRecNow].isRescatteredIncoming() ) ) {
560 for (
int j = 0; j < sizeAll; ++j) {
561 int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA);
562 if (iRecNow == iRad)
continue;
563 if ( ( j >= sizeIn && event[iRecNow].acol() == colTag
564 && event[iRecNow].isFinal() )
565 || ( j < sizeIn && event[iRecNow].col() == colTag
566 && !event[iRecNow].isRescatteredIncoming() ) ) {
573 if ( iPartner == 0 ) {
574 infoPtr->errorMsg(
"Error in DireSpace::setupQCDdip: "
575 "failed to locate any recoiling partner");
581 double pTmax =
event[iRad].scale();
583 if (iSys == 0 || (iSys == 1 && doSecondHard)) pTmax *= pTmaxFudge;
584 else if (sizeIn > 0) pTmax *= pTmaxFudgeMPI;
586 }
else pTmax = m( event[iRad], event[iPartner]);
589 if ( abs(event[iRad].status()) > 20 && abs(event[iRad].status()) < 24
590 && settingsPtr->flag(
"Beams:setProductionScalesFromLHEF")
591 &&
event[iRad].scale() > 0.)
592 pTmax =
event[iRad].scale();
595 double mups = infoPtr->getScalesAttribute(
"mups");
596 if ( abs(event[iRad].status()) > 20
597 && abs(event[iRad].status()) < 24
598 && settingsPtr->flag(
"Beams:setProductionScalesFromLHEF")
602 int colType = (
event[iRad].id() == 21) ? 2 * colSign : colSign;
603 dipEnd.push_back( DireSpaceEnd( iSys, side, iRad, iPartner, pTmax, colType,
605 dipEnd.back().init(event);
613 void DireSpace::getGenDip(
int iSys,
int side,
const Event& event,
614 bool limitPTmaxIn, vector<DireSpaceEnd>& dipEnds) {
617 int iRad = (iSys > -1) ? ((side == 1) ? getInA(iSys) : getInB(iSys))
619 int sizeAllA = (iSys > -1) ? partonSystemsPtr->sizeAll(iSys) :
event.size();
620 int sizeOut = (iSys > -1) ? partonSystemsPtr->sizeOut(iSys) :
event.size();
621 int sizeAll = (iSys > -1) ? sizeAllA : event.size();
622 int sizeIn = (iSys > -1) ? sizeAll - sizeOut : 0;
623 int sizeInA = (iSys > -1) ? sizeAllA - sizeIn - sizeOut : 0;
625 for (
int i = 0; i < sizeAll; ++i) {
626 int iRecNow = (iSys > -1) ?
627 partonSystemsPtr->getAll(iSys, i + sizeInA) : i;
628 if ( !event[iRecNow].isFinal()
629 && event[iRecNow].mother1() != 1
630 && event[iRecNow].mother1() != 2)
continue;
632 if ( iRecNow == iRad)
continue;
636 for (
int j = 0; j < int(dipEnds.size()); ++j)
637 if ( dipEnds[j].iRadiator == iRad && dipEnds[j].iRecoiler == iRecNow )
639 if (
int(iDip.size()) > 0) {
640 for (
int j = 0; j < int(iDip.size()); ++j)
641 updateAllowedEmissions(event, &dipEnds[iDip[j]]);
645 double pTmax = abs(2.*event[iRad].p()*event[iRecNow].p());
647 if (iSys == 0 || (iSys == 1 && doSecondHard)) pTmax *= pTmaxFudge;
648 else if (sizeIn > 0) pTmax *= pTmaxFudgeMPI;
649 }
else pTmax = m( event[iRad], event[iRecNow]);
651 appendDipole( event, iSys, side, iRad, iRecNow, pTmax, 0, 0, 0, 0,
true, 0,
652 vector<int>(), vector<double>(), dipEnds);
664 void DireSpace::getQCDdip(
int iRad,
int colTag,
int colSign,
665 const Event& event, vector<DireSpaceEnd>& dipEnds) {
671 for (
int iRecNow = 0; iRecNow <
event.size(); ++iRecNow) {
672 if (iRecNow == iRad)
continue;
673 if ( ( event[iRecNow].col() == colTag && event[iRecNow].isFinal() )
674 || (
event[iRecNow].acol() == colTag && !
event[iRecNow].isFinal() ) ) {
682 for (
int iRecNow = 0; iRecNow <
event.size(); ++iRecNow) {
683 if (iRecNow == iRad)
continue;
684 if ( ( event[iRecNow].acol() == colTag && event[iRecNow].isFinal() )
685 || (
event[iRecNow].col() == colTag && !
event[iRecNow].isFinal() ) ) {
693 double pTmax = abs(2.*event[iRad].p()*event[iPartner].p());
694 int side = (
event[iRad].pz() > 0.) ? 1 : 2;
695 int colType = (
event[iRad].id() == 21) ? 2 * colSign : colSign;
698 dipEnds.push_back( DireSpaceEnd( 0, side, iRad, iPartner, pTmax, colType));
699 dipEnds.back().init(event);
707 bool DireSpace::appendDipole(
const Event& state,
int sys,
int side,
708 int iRad,
int iRecNow,
double pTmax,
int colType,
int chgType,
int weakType,
709 int MEtype,
bool normalRecoil,
int weakPol, vector<int> iSpectator,
710 vector<double> mass, vector<DireSpaceEnd>& dipEnds) {
713 if (colType == 0 && state[iRad].colType() != 0) {
714 vector<int> shared = sharedColor(state[iRad], state[iRecNow]);
718 for (
int i=0; i < int(shared.size()); ++i) {
719 if ( state[iRad].colType() == 2 && state[iRad].col() == shared[i])
721 if ( state[iRad].colType() == 2 && state[iRad].acol() == shared[i])
723 if ( state[iRad].colType() == 1 && state[iRad].id() > 0
724 && state[iRad].col() == shared[i])
726 if ( state[iRad].colType() ==-1 && state[iRad].id() < 0
727 && state[iRad].acol() == shared[i])
730 for (
int j=0; j < int(dipEnds.size()); ++j) {
731 if ( dipEnds[j].iRadiator == iRad && dipEnds[j].iRecoiler == iRecNow
732 && dipEnds[j].colType == colTypeNow) { found =
true;
break; }
737 colType = colTypeNow;
741 DireSpaceEnd dipNow = DireSpaceEnd( sys, side, iRad, iRecNow, pTmax, colType,
742 chgType, weakType, MEtype, normalRecoil, weakPol,
743 DireSingleColChain(), iSpectator, mass);
745 dipNow.clearAllowedEmt();
747 if (updateAllowedEmissions(state, &dipNow)) {
748 dipEnds.push_back(dipNow);
757 vector<int> DireSpace::sharedColor(
const Particle& rad,
const Particle& rec) {
759 int radCol(rad.col()), radAcl(rad.acol()),
760 recCol(rec.col()), recAcl(rec.acol());
761 if ( rad.isFinal() && rec.isFinal() ) {
762 if (radCol != 0 && radCol == recAcl) ret.push_back(radCol);
763 if (radAcl != 0 && radAcl == recCol) ret.push_back(radAcl);
764 }
else if ( rad.isFinal() && !rec.isFinal() ) {
765 if (radCol != 0 && radCol == recCol) ret.push_back(radCol);
766 if (radAcl != 0 && radAcl == recAcl) ret.push_back(radAcl);
767 }
else if (!rad.isFinal() && rec.isFinal() ) {
768 if (radCol != 0 && radCol == recCol) ret.push_back(radCol);
769 if (radAcl != 0 && radAcl == recAcl) ret.push_back(radAcl);
770 }
else if (!rad.isFinal() && !rec.isFinal() ) {
771 if (radCol != 0 && radCol == recAcl) ret.push_back(radCol);
772 if (radAcl != 0 && radAcl == recCol) ret.push_back(radAcl);
781 void DireSpace::saveSiblings(
const Event& state,
int iSys) {
783 int sizeAllSys = partonSystemsPtr->sizeSys();
784 for (
int iSystem=0; iSystem < sizeAllSys; ++iSystem) {
786 if (iSys > -1 && iSystem != iSys)
continue;
788 vector<int> q, qb, g;
789 int sizeSystem(partonSystemsPtr->sizeAll(iSystem)), nFinal(0);
790 for (
int i = 0; i < sizeSystem; ++i) {
792 int iPos = partonSystemsPtr->getAll(iSystem, i);
793 if ( state[iPos].isFinal()) nFinal++;
795 if (!state[iPos].isFinal()
796 && state[iPos].mother1() != 1
797 && state[iPos].mother1() != 2)
continue;
799 if ( state[iPos].isFinal() && state[iPos].colType() == 1
800 && find(q.begin(),q.end(),iPos) == q.end())
803 if (!state[iPos].isFinal() && state[iPos].colType() ==-1
804 && find(q.begin(),q.end(),iPos) == q.end())
807 if ( state[iPos].isFinal() && state[iPos].colType() ==-1
808 && find(qb.begin(),qb.end(),iPos) == qb.end())
811 if (!state[iPos].isFinal() && state[iPos].colType() == 1
812 && find(qb.begin(),qb.end(),iPos) == qb.end())
815 if ( abs(state[iPos].colType()) == 2
816 && find(g.begin(),g.end(),iPos) == g.end())
821 DireColChains chains;
823 for (
int i = 0; i < int(q.size()); ++i) {
824 if (chains.chainOf(q[i]).size() != 0)
continue;
825 chains.addChain( DireSingleColChain(q[i],state, partonSystemsPtr));
828 for (
int i = 0; i < int(qb.size()); ++i) {
829 if (chains.chainOf(qb[i]).size() != 0)
continue;
830 chains.addChain( DireSingleColChain(qb[i],state, partonSystemsPtr));
833 for (
int i = 0; i < int(g.size()); ++i) {
834 if (chains.chainOf(g[i]).size() != 0)
continue;
835 chains.addChain( DireSingleColChain(g[i],state, partonSystemsPtr));
839 for (
int iDip = 0; iDip < int(dipEnd.size()); ++iDip) {
840 if (dipEnd[iDip].system != iSystem)
continue;
841 if (dipEnd[iDip].colType == 0) dipEnd[iDip].clearSiblings();
843 int col = dipEnd[iDip].colType > 0
844 ? state[dipEnd[iDip].iRadiator].col()
845 : state[dipEnd[iDip].iRadiator].acol();
846 dipEnd[iDip].setSiblings(chains.chainFromCol( dipEnd[iDip].iRadiator,
859 void DireSpace::updateDipoles(
const Event& state,
int iSys) {
863 for (
int iDip = 0; iDip < int(dipEnd.size()); ++iDip) {
864 if (!updateAllowedEmissions(state, &dipEnd[iDip])) iRemove.push_back(iDip);
865 dipEnd[iDip].init(state);
868 for (
int i = iRemove.size()-1; i >= 0; --i) {
869 dipEnd[iRemove[i]] = dipEnd.back();
874 for (
int iDip = 0; iDip < int(dipEnd.size()); ++iDip) {
875 DireSpaceEnd* dip = &dipEnd[iDip];
876 int iRad = dip->iRadiator;
877 int iRecNow = dip->iRecoiler;
879 if (dip->colType == 0 && state[iRad].colType() != 0) {
880 vector<int> shared = sharedColor(state[iRad], state[iRecNow]);
884 for (
int i=0; i < int(shared.size()); ++i) {
885 if ( state[iRad].colType() == 2 && state[iRad].col() == shared[i])
887 if ( state[iRad].colType() == 2 && state[iRad].acol() == shared[i])
889 if ( state[iRad].colType() == 1 && state[iRad].id() > 0
890 && state[iRad].col() == shared[i])
892 if ( state[iRad].colType() ==-1 && state[iRad].id() < 0
893 && state[iRad].acol() == shared[i])
896 dip->colType = colTypeNow;
900 saveSiblings(state, iSys);
906 bool DireSpace::updateAllowedEmissions(
const Event& state, DireSpaceEnd* dip) {
908 dip->clearAllowedEmt();
910 return appendAllowedEmissions(state, dip);
917 bool DireSpace::appendAllowedEmissions(
const Event& state, DireSpaceEnd* dip) {
921 bool isAllowed =
false;
922 int iRad(dip->iRadiator), iRecNow(dip->iRecoiler);
923 pair<int,int> iRadRec(make_pair(iRad, iRecNow));
924 pair<int,int> iRecRad(make_pair(iRecNow, iRad));
926 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
927 it != splits.end(); ++it ) {
930 bool allowed = it->second->useFastFunctions()
931 ? it->second->canRadiate(state,iRad,iRecNow)
932 : it->second->canRadiate(state,iRadRec,bool_settings);
933 if (!allowed)
continue;
936 vector<int> re = it->second->radAndEmt( state[iRad].
id(), dip->colType);
938 for (
int iEmtAft=1; iEmtAft < int(re.size()); ++iEmtAft) {
939 int idEmtAft = re[iEmtAft];
940 if (it->second->is_qcd) {
941 idEmtAft = abs(idEmtAft);
942 if (idEmtAft<10) idEmtAft = 1;
945 if (!it->second->isPartial()) {
946 dip->appendAllowedEmt(idEmtAft);
951 bool isPartialFractioned =
false;
952 for ( unordered_map<string,DireSplitting*>::iterator itRec =
954 itRec != splits.end(); ++itRec ) {
956 if ( isPartialFractioned )
break;
957 bool allowedRec = itRec->second->useFastFunctions()
958 ? itRec->second->canRadiate(state,iRecNow,iRad)
959 : itRec->second->canRadiate(state,iRecRad,bool_settings);
960 if (!allowedRec)
continue;
964 = state[iRecNow].isFinal() ? -dip->colType : dip->colType;
966 = itRec->second->radAndEmt( state[iRecNow].
id(), colTypeRec);
968 for (
int iEmtAftRec=1; iEmtAftRec<int(reRec.size()); ++iEmtAftRec) {
969 int idEmtAftRec = reRec[iEmtAftRec];
970 if (itRec->second->is_qcd) {
971 idEmtAftRec = abs(idEmtAftRec);
972 if (idEmtAftRec<10) idEmtAftRec = 1;
974 if (idEmtAftRec == idEmtAft) { isPartialFractioned =
true;
break;}
979 if (isPartialFractioned) {
980 dip->appendAllowedEmt(idEmtAft);
995 void DireSpace::update(
int iSys,
Event& event,
bool) {
998 int in1 = getInA(iSys);
999 int in2 = getInB(iSys);
1002 bool canRadiate1 = !(
event[in1].isRescatteredIncoming()) && doQCDshower;
1003 bool canRadiate2 = !(
event[in2].isRescatteredIncoming()) && doQCDshower;
1010 int colTag =
event[in1].col();
1011 if (canRadiate1 && colTag > 0) setupQCDdip( iSys, 1, colTag, 1, event,
1014 int acolTag =
event[in1].acol();
1015 if (canRadiate1 && acolTag > 0) setupQCDdip( iSys, 1, acolTag, -1, event,
1019 colTag =
event[in2].col();
1020 if (canRadiate2 && colTag > 0) setupQCDdip( iSys, 2, colTag, 1, event,
1023 acolTag =
event[in2].acol();
1024 if (canRadiate2 && acolTag > 0) setupQCDdip( iSys, 2, acolTag, -1, event,
1028 getGenDip( iSys, 1, event,
false, dipEnd);
1029 getGenDip( iSys, 2, event,
false, dipEnd);
1032 updateDipoles(event, iSys);
1040 double DireSpace::newPoint(
const Event& inevt) {
1042 int asOrderSave = alphaSorder;
1043 double pT2minMECsSave = pT2minMECs;
1046 double tFreeze = 1.;
1055 splittingNowName=
"";
1056 splittingSelName=
"";
1057 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
1058 it != splits.end(); ++it ) it->second->splitInfo.clear();
1059 splitInfoSel.clear();
1062 auxSel = overSel = auxNow = overNow = 0.;
1067 while (pT2sel==0.) {
1070 if (nTrials>=nTrialsMax)
break;
1073 for (
int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
1076 dipEndNow = &dipEnd[iDipEnd];
1077 iSysNow = dipEndNow->system;
1080 double pT2begDip = m2Max(iDipNow, event);
1082 if (pT2begDip > pT2sel) {
1083 double pT2endDip = 0.;
1086 sideA = ( abs(dipEndNow->side) == 1 );
1087 bool finalRec =
event[dipEndNow->iRecoiler].isFinal();
1088 BeamParticle& beamNow = (sideA) ? *beamAPtr : *beamBPtr;
1089 BeamParticle& beamRec = (sideA) ? *beamBPtr : *beamAPtr;
1090 iNow = beamNow[iSysNow].iPos();
1091 iRec = (finalRec) ? dipEndNow->iRecoiler
1092 : beamRec[iSysNow].iPos();
1093 idDaughter = beamNow[iSysNow].id();
1094 xDaughter = beamNow[iSysNow].x();
1095 x1Now = (sideA) ? xDaughter : beamRec[iSysNow].x();
1096 x2Now = (sideA) ? beamRec[iSysNow].x() : xDaughter;
1099 m2Rec = (dipEndNow->normalRecoil) ? 0. : event[iRec].m2();
1100 m2Dip = abs(2.*event[iNow].p()*event[iRec].p());
1103 dipEndNow->m2Dip = m2Dip;
1105 dipEndNow->pT2 = 0.0;
1106 dipEndNow->z = -1.0;
1107 dipEndNow->phi = -1.0;
1109 dipEndNow->sa1 = 0.0;
1110 dipEndNow->xa = -1.0;
1111 dipEndNow->phia1 = -1.0;
1114 if (pT2begDip > pT2endDip) {
1116 if ( dipEndNow->canEmit() ) pT2nextQCD( pT2begDip, 0.,
1117 *dipEndNow, event, 0., tFreeze,
true);
1120 if (dipEndNow->pT2 > pT2sel) {
1121 pT2sel = dipEndNow->pT2;
1124 dipEndSel = dipEndNow;
1125 splittingSelName = splittingNowName;
1126 splittingSel = splits[splittingSelName];
1127 splitInfoSel.store(splits[splittingSelName]->splitInfo);
1128 kernelSel = kernelNow;
1131 boostSel = boostNow;
1139 bool hasBranched =
false;
1140 if ( event[dipEndSel->iRecoiler].isFinal())
1141 hasBranched = branch_IF(event,
true, &splitInfoSel);
1142 else hasBranched = branch_II(event,
true, &splitInfoSel);
1147 splitInfoSel.clear();
1154 alphaSorder = asOrderSave;
1155 pT2minMECs = pT2minMECsSave;
1158 if (dipEndSel == 0)
return 0.;
1161 return sqrt(pT2sel);
1169 double DireSpace::pTnext(
Event& event,
double pTbegAll,
double pTendAll,
1170 int nRadIn,
bool doTrialIn) {
1172 direInfoPtr->message(1) <<
"Next ISR starting from " << pTbegAll << endl;
1175 sCM = m2( beamAPtr->p(), beamBPtr->p());
1177 pTbegRef = pTbegAll;
1181 double pT2sel = pow2(pTendAll);
1185 splittingNowName=
"";
1186 splittingSelName=
"";
1187 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
1188 it != splits.end(); ++it ) it->second->splitInfo.clear();
1189 splitInfoSel.clear();
1192 auxSel = overSel = auxNow = overNow = 0.;
1196 doTrialNow = doTrialIn;
1199 for (
int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
1201 dipEndNow = &dipEnd[iDipEnd];
1202 iSysNow = dipEndNow->system;
1203 double pTbegDip = min( pTbegAll, dipEndNow->pTmax );
1206 int nfmax = settingsPtr->mode(
"DireSpace:nFinalMax");
1208 for (
int i=0; i <
event.size(); ++i)
1209 if (event[i].isFinal()) nFinal++;
1210 if (nfmax > -10 && nFinal > nfmax)
continue;
1213 double pT2begDip = pow2(pTbegDip);
1215 if (pT2begDip > pT2sel) {
1216 double pT2endDip = 0.;
1219 pT2endDip = max( pT2cutMin(dipEndNow), pTendAll*pTendAll);
1220 pT2endDip = max(pT2endDip, pT2sel);
1223 sideA = ( abs(dipEndNow->side) == 1 );
1224 bool finalRec =
event[dipEndNow->iRecoiler].isFinal();
1225 BeamParticle& beamNow = (sideA) ? *beamAPtr : *beamBPtr;
1226 BeamParticle& beamRec = (sideA) ? *beamBPtr : *beamAPtr;
1227 iNow = beamNow[iSysNow].iPos();
1228 iRec = (finalRec) ? dipEndNow->iRecoiler
1229 : beamRec[iSysNow].iPos();
1230 idDaughter = beamNow[iSysNow].id();
1231 xDaughter = beamNow[iSysNow].x();
1232 x1Now = (sideA) ? xDaughter : beamRec[iSysNow].x();
1233 x2Now = (sideA) ? beamRec[iSysNow].x() : xDaughter;
1236 m2Rec = (dipEndNow->normalRecoil) ? 0. : event[iRec].m2();
1237 m2Dip = abs(2.*event[iNow].p()*event[iRec].p());
1240 dipEndNow->m2Dip = m2Dip;
1242 dipEndNow->pT2 = 0.0;
1243 dipEndNow->z = -1.0;
1244 dipEndNow->phi = -1.0;
1246 dipEndNow->sa1 = 0.0;
1247 dipEndNow->xa = -1.0;
1248 dipEndNow->phia1 = -1.0;
1251 if (pT2begDip > pT2endDip) {
1253 if ( dipEndNow->canEmit() ) pT2nextQCD( pT2begDip, pT2endDip,
1257 if (dipEndNow->pT2 > pT2sel) {
1258 pT2sel = dipEndNow->pT2;
1261 dipEndSel = dipEndNow;
1262 splittingSelName = splittingNowName;
1263 splittingSel = splits[splittingSelName];
1264 splitInfoSel.store(splits[splittingSelName]->splitInfo);
1265 kernelSel = kernelNow;
1268 boostSel = boostNow;
1276 for ( unordered_map<
string, multimap<double,double> >::iterator
1277 itR = rejectProbability.begin(); itR != rejectProbability.end(); ++itR)
1278 weights->insertWeights(acceptProbability[itR->first], itR->second,
1280 for ( unordered_map<
string, multimap<double,double> >::iterator
1281 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
1283 for ( unordered_map<
string, map<double,double> >::iterator
1284 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
1287 resetOverheadFactors();
1292 if (dipEndSel == 0) {
1293 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
1294 it != splits.end(); ++it ) {
1295 if (it->second->fsr && it->second->fsr->dipSel == 0) {
1296 it->second->fsr->finalize(event);
1303 return (dipEndSel == 0) ? 0. : sqrt(pT2sel);
1312 double DireSpace::noEmissionProbability(
double pTbegAll,
double pTendAll,
1313 double m2dip,
int idA,
int type,
double s,
double x) {
1320 pTbegRef = pTbegAll;
1323 double x2 = m2dip/s/x1;
1326 splittingNowName=
"";
1327 splittingSelName=
"";
1328 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
1329 it != splits.end(); ++it ) it->second->splitInfo.clear();
1333 state.init(
"(dummy event)", particleDataPtr);
1335 Vec4 pA(0., 0., 0.5*sqrt(m2dip), 0.5*sqrt(m2dip)), pB;
1336 if (type < 0) pB.p(0., 0.,-0.5*sqrt(m2dip), 0.5*sqrt(m2dip));
1337 if (type > 0) pB.p(0., 0.,0.5*sqrt(m2dip), 0.5*sqrt(m2dip));
1341 if (particleDataPtr->colType(idA) == 1) {colA = 1; acolA = 0;}
1342 if (particleDataPtr->colType(idA) ==-1) {colA = 0; acolA = 1;}
1345 state.append( 0, 0, 0, 0, 0, 0, 0, 0, pA+pB, 0.0, sqrt(m2dip) );
1346 state.append( idA, -21, 0, 0, 0, 0, colA, acolA, pA, 0.0, sqrt(m2dip) );
1349 int idB = (idA == 21) ? 21
1350 : ((type < 0) ? -idA : idA);
1351 vector<int> recids; recids.push_back(idB);
1354 for (
unsigned int i = 0; i < recids.size(); ++i) {
1355 int colB(2), acolB(1);
1357 && particleDataPtr->colType(idA) == 1
1358 && particleDataPtr->colType(recids[i]) ==-1) {colB = 0; acolB = colA;}
1360 && particleDataPtr->colType(idA) ==-1
1361 && particleDataPtr->colType(recids[i]) == 1) {colB = acolA; acolB = 0;}
1364 && particleDataPtr->colType(idA) == 2
1365 && particleDataPtr->colType(recids[i]) ==-1) {colB = 0; acolB = colA;}
1367 && particleDataPtr->colType(idA) == 2
1368 && particleDataPtr->colType(recids[i]) == 1) {colB = acolA; acolB = 0;}
1370 if (type < 0) state.append( recids[i], -21, 0, 0, 0, 0, colB, acolB,
1371 pB, 0.0, sqrt(m2dip) );
1372 if (type > 0) state.append( recids[i], 23, 0, 0, 0, 0, colB, acolB,
1373 pB, 0.0, sqrt(m2dip) );
1374 recpos.push_back(i+1);
1379 beamAPtr->append( 1, idA, x1);
1380 beamBPtr->append( 2, idB, x2);
1381 beamAPtr->xfISR( 0, idA, x1, pTbegAll*pTbegAll);
1382 int vsc1 = beamAPtr->pickValSeaComp();
1383 beamBPtr->xfISR( 0, idB, x2, pTbegAll*pTbegAll);
1384 int vsc2 = beamBPtr->pickValSeaComp();
1385 infoPtr->setValence( (vsc1 == -3), (vsc2 == -3));
1388 partonSystemsPtr->clear();
1389 partonSystemsPtr->addSys();
1390 partonSystemsPtr->setInA(0, 1);
1391 partonSystemsPtr->setInB(0, 2);
1392 partonSystemsPtr->setSHat( 0, m2dip);
1393 partonSystemsPtr->setPTHat( 0, pTbegAll);
1397 vector<DireSpaceEnd> dipEnds;
1398 int colTag = state[in1].col();
1399 if (colTag > 0) getQCDdip( in1, colTag, 1, state, dipEnds);
1400 int acolTag = state[in1].acol();
1401 if (acolTag > 0) getQCDdip( in1, acolTag, -1, state, dipEnds);
1404 double startingScale = pTbegAll;
1412 state.scale(startingScale);
1415 double minScale = pTendAll;
1417 mergingHooksPtr->setShowerStoppingScale(minScale);
1422 if (minScale >= startingScale)
break;
1425 double pTtrial = pTnext( dipEnds, state, startingScale, minScale,
1428 pair<double,double> wtShower
1429 = weights->getWeight( (pTtrial <= 0.) ? pow2(minScale) : pow2(pTtrial));
1431 double enhancement = 1.;
1432 if ( pTtrial > minScale) enhancement
1433 = weights->getTrialEnhancement( pow2(pTtrial));
1436 weights->clearTrialEnhancements();
1439 if ( pTtrial < minScale ) { wt *= wtShower.second;
break;}
1442 startingScale = pTtrial;
1444 if ( pTtrial > minScale) wt *= wtShower.first*wtShower.second
1445 * (1.-1./enhancement);
1446 if ( wt == 0.)
break;
1447 if ( pTtrial > minScale)
continue;
1456 partonSystemsPtr->clear();
1464 double DireSpace::pTnext( vector<DireSpaceEnd> dipEnds,
Event event,
1465 double pTbegAll,
double pTendAll,
double m2dip,
int,
double s,
1469 double x2 = m2dip/s/x1;
1473 double pT2sel = pow2(pTendAll);
1477 splittingNowName=
"";
1478 splittingSelName=
"";
1479 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
1480 it != splits.end(); ++it ) it->second->splitInfo.clear();
1481 splitInfoSel.clear();
1484 auxSel = overSel = auxNow = overNow = 0.;
1489 splits = splittingsPtr->getSplittings();
1491 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
1492 it != splits.end(); ++it ) overhead.insert(make_pair(it->first,1.));
1495 nProposedPT.clear();
1496 if ( nProposedPT.find(iSys) == nProposedPT.end() )
1497 nProposedPT.insert(make_pair(iSys,0));
1499 splittingSelName=
"";
1500 splittingNowName=
"";
1504 for ( unordered_map<
string, multimap<double,double> >::iterator
1505 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
1507 for ( unordered_map<
string, map<double,double> >::iterator
1508 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
1512 for (
int iDipEnd = 0; iDipEnd < int(dipEnds.size()); ++iDipEnd) {
1515 dipEndNow = &dipEnds[iDipEnd];
1516 double pTbegDip = min( pTbegAll, dipEndNow->pTmax );
1520 for (
int i=0; i <
event.size(); ++i)
1521 if (event[i].isFinal()) nFinal++;
1522 if (nFinalMax > -10 && nFinal > nFinalMax)
continue;
1525 double pT2begDip = pow2(pTbegDip);
1526 double pT2endDip = 0.;
1528 pT2endDip = max( pT2cutMin(dipEndNow), pTendAll*pTendAll);
1529 pT2endDip = max(pT2endDip, pT2sel);
1532 sideA = ( abs(dipEndNow->side) == 1 );
1533 iNow = dipEndNow->iRadiator;
1534 iRec = dipEndNow->iRecoiler;
1535 idDaughter =
event[dipEndNow->iRadiator].id();
1537 x1Now = (sideA) ? x1 : x2;
1538 x2Now = (sideA) ? x2 : x1;
1540 m2Rec = (dipEndNow->normalRecoil) ? 0. : event[iRec].m2();
1541 m2Dip = abs(2.*event[iNow].p()*event[iRec].p());
1544 dipEndNow->m2Dip = m2Dip;
1546 dipEndNow->pT2 = 0.0;
1547 dipEndNow->z = -1.0;
1548 dipEndNow->phi = -1.0;
1550 dipEndNow->sa1 = 0.0;
1551 dipEndNow->xa = -1.0;
1552 dipEndNow->phia1 = -1.0;
1555 if (pT2begDip > pT2endDip) {
1557 if ( dipEndNow->canEmit() ) pT2nextQCD( pT2begDip, pT2endDip,
1561 if (dipEndNow->pT2 > pT2sel) {
1562 pT2sel = dipEndNow->pT2;
1565 dipEndSel = dipEndNow;
1566 splittingSelName = splittingNowName;
1567 splittingSel = splits[splittingSelName];
1568 splitInfoSel.store(splits[splittingSelName]->splitInfo);
1569 kernelSel = kernelNow;
1572 boostSel = boostNow;
1580 for ( unordered_map<
string, multimap<double,double> >::iterator
1581 itR = rejectProbability.begin(); itR != rejectProbability.end(); ++itR){
1582 weights->insertWeights(acceptProbability[itR->first], itR->second,
1586 for ( unordered_map<
string, multimap<double,double> >::iterator
1587 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
1589 for ( unordered_map<
string, map<double,double> >::iterator
1590 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
1593 resetOverheadFactors();
1596 return (dipEndSel == 0) ? 0. : sqrt(pT2sel);
1602 double DireSpace::enhanceOverestimateFurther(
string name,
int,
1605 if (tOld < pT2minEnhance)
return 1.;
1606 double enhance = weights->enhanceOverestimate(name);
1613 double DireSpace::overheadFactors(
string name,
int idDau,
bool isValence,
1614 double m2dip,
double pT2Old ) {
1619 if (isValence && name.find(
"isr_qcd_1->1&21") != string::npos)
1620 factor *= log(max(2.71828,16/(pT2Old/m2dip)));
1623 if (name.find(
"isr_qcd_21->1&1") != string::npos)
1624 factor *= log(max(2.71828,log(max(2.71828,m2dip/pT2Old))
1625 + pow(m2dip/pT2Old,3./2.)));
1629 if (name.find(
"isr_qcd_1->1&21") != string::npos && !isValence)
1631 if (name.find(
"isr_qcd_1->21&1") != string::npos && !isValence)
1633 if (name.find(
"isr_qcd_21->1&1") != string::npos)
1635 if (name.find(
"isr_qcd_21->21&21a") != string::npos && pT2Old < 2.0)
1637 if (name.find(
"isr_qcd_21->21&21b") != string::npos && pT2Old < 2.0)
1641 if (pT2Old < pT2min*1.25) MARGIN = 1.0;
1646 if ( abs(idDau) == 4 && name.find(
"isr_qcd_21->1&1") != string::npos
1647 && pT2Old < 2.*m2cPhys) factor *= 1. / max(0.01, abs(pT2Old - m2cPhys));
1648 if ( abs(idDau) == 5 && name.find(
"isr_qcd_21->1&1") != string::npos
1649 && pT2Old < 2.*m2bPhys) factor *= 1. / max(0.01, abs(pT2Old - m2bPhys));
1652 if ( overhead.find(name) != overhead.end() ) factor *= overhead[name];
1662 void DireSpace::getNewOverestimates(
int idDau, DireSpaceEnd* dip,
1663 const Event& state,
double tOld,
double xDau,
double zMinAbs,
1664 double zMaxAbs, multimap<double,string>& newOverestimates ) {
1667 BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
1668 bool isValence = (usePDF) ? beam[iSysNow].isValence() :
false;
1669 pair<int,int> iRadRec(make_pair(dip->iRadiator, dip->iRecoiler));
1674 for ( unordered_map<string,DireSplitting*>::iterator it = splits.begin();
1675 it != splits.end(); ++it ) {
1677 string name = it->first;
1680 bool allowed = it->second->useFastFunctions()
1681 ? it->second->canRadiate(state,dip->iRadiator,dip->iRecoiler)
1682 : it->second->canRadiate(state,iRadRec,bool_settings);
1685 if (!allowed)
continue;
1688 vector<int> re = it->second->radAndEmt(state[dip->iRadiator].id(),
1690 if (
int(re.size()) < 2)
continue;
1692 for (
int iEmtAft=1; iEmtAft < int(re.size()); ++iEmtAft) {
1693 int idEmtAft = re[iEmtAft];
1694 if (it->second->is_qcd) {
1695 idEmtAft = abs(idEmtAft);
1696 if (idEmtAft<10) idEmtAft = 1;
1698 if (find(dip->allowedEmissions.begin(), dip->allowedEmissions.end(),
1699 idEmtAft) == dip->allowedEmissions.end() ) allowed =
false;
1701 if ( pT2cut(idEmtAft) > tOld) allowed =
false;
1704 if (!allowed)
continue;
1707 if ( tOld < 4.*m2bPhys && abs(idDau) == 5
1708 && it->second->nEmissions() == 2)
continue;
1709 else if ( tOld < 4.*m2cPhys && abs(idDau) == 4
1710 && it->second->nEmissions() == 2)
continue;
1713 int order = kernelOrder;
1715 bool hasInA = (getInA(dip->system) != 0);
1716 bool hasInB = (getInB(dip->system) != 0);
1717 if (dip->system != 0 && hasInA && hasInB) order = kernelOrderMPI;
1719 it->second->splitInfo.set_pT2Old ( tOld );
1720 it->second->splitInfo.storeRadBef(state[dip->iRadiator]);
1721 it->second->splitInfo.storeRecBef(state[dip->iRecoiler]);
1724 if (!it->second->aboveCutoff( tOld, state[dip->iRadiator],
1725 state[dip->iRecoiler], dip->system, partonSystemsPtr))
continue;
1728 double wt = it->second->overestimateInt(zMinAbs, zMaxAbs, tOld,
1733 double pdfRatio = getPDFOverestimates(idDau, tOld, xDau, it->first,
1734 false, -1., re[0], re[0]);
1741 overheadFactors(name, idDau, isValence, dip->m2Dip, tOld);
1745 double enhanceFurther = enhanceOverestimateFurther(name, idDau, tOld);
1746 wt *= enhanceFurther;
1748 if (!dryrun && it->second->hasMECBef(state, tOld)) wt *= KERNEL_HEADROOM;
1750 for (
int i=0; i < state.size(); ++i)
if (state[i].isFinal()) nFinal++;
1751 if (!dryrun) wt *= it->second->overhead
1752 (dip->m2Dip*xDau, state[dip->iRadiator].id(), nFinal);
1758 newOverestimates.insert(make_pair(sum,name));
1769 double DireSpace::getPDFOverestimates(
int idDau,
double tOld,
double xDau,
1770 string name,
bool pickMother,
double RN,
int& idMother,
int& idSister) {
1772 BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
1773 DireSplitting* splitNow = splits[name];
1776 double PDFscale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac*tOld;
1777 PDFscale2 = max(PDFscale2, pT2min);
1778 bool inD = (hasPDF(idDau)) ? beam.insideBounds(xDau, PDFscale2) :
1780 double xPDFdaughter = getXPDF( idDau, xDau, PDFscale2, iSysNow, &beam);
1782 if (abs(idDau) == 4 && m2cPhys > 0. && tOld < 4.*m2cPhys) {
1783 double xPDFthres = getXPDF( idDau, xDau, m2cPhys+0.1, iSysNow, &beam);
1784 xPDFdaughter = min(xPDFdaughter, xPDFthres);
1786 if (abs(idDau) == 5 && m2bPhys > 0. && tOld < 4.*m2bPhys) {
1787 double xPDFthres = getXPDF( idDau, xDau, m2bPhys, iSysNow, &beam);
1788 xPDFdaughter = min(xPDFdaughter, xPDFthres);
1793 double pdfRatio = 1.;
1794 if (splitNow->is(splittingsPtr->isrQCD_1_to_21_and_1)) {
1797 double xPDFmotherSum = 0.;
1798 double xPDFmother[21] = {0.};
1799 for (
int iQuark = -nQuarkIn; iQuark <= nQuarkIn; ++iQuark) {
1801 xPDFmother[10] = 0.;
1803 xPDFmother[iQuark+10] = getXPDF(iQuark,xDau,PDFscale2,iSysNow,&beam);
1804 xPDFmotherSum += xPDFmother[iQuark+10];
1811 double temp = xPDFmotherSum * RN;
1812 idMother = -nQuarkIn - 1;
1813 do { temp -= xPDFmother[(++idMother) + 10]; }
1814 while (temp > 0. && idMother < nQuarkIn);
1815 idSister = idMother;
1816 pdfRatio = xPDFmother[idMother + 10]/xPDFdaughter;
1818 pdfRatio = xPDFmotherSum/xPDFdaughter;
1822 }
else if (splitNow->is(splittingsPtr->isrQCD_21_to_1_and_1)){
1824 double xPDFmother = getXPDF(21, xDau, PDFscale2, iSysNow, &beam);
1825 if ( xPDFmother != 0. && abs(xPDFmother) < tinypdf(xDau) ) {
1826 int sign = (xPDFmother >= 0.) ? 1 : -1;
1827 xPDFmother = sign*tinypdf(xDau);
1829 pdfRatio = xPDFmother / xPDFdaughter;
1832 }
else if (splitNow->is(splittingsPtr->isrQCD_1_to_2_and_1_and_2)) {
1835 multimap<double, pair<int,double> > xPDFmother;
1836 double xPDFmotherSum = 0.;
1837 for (
int i =-nQuarkIn; i <= nQuarkIn; ++i)
1838 if (abs(i) != abs(idDau) && i != 0) {
1839 double temp = getXPDF( i, xDau, PDFscale2, iSysNow, &beam);
1841 if (particleDataPtr->isHadron(beam.id()) && (i == 1 || i == 2)) {
1842 if (abs(idDau) == 4 && m2cPhys > 0. && tOld < 4.*m2cPhys) {
1843 double xPDFval = getXPDF(i, 0.25, PDFscale2, iSysNow, &beam);
1844 temp = max(temp, xPDFval);
1846 if (abs(idDau) == 5 && m2bPhys > 0. && tOld < 4.*m2bPhys) {
1847 double xPDFval = getXPDF( i, 0.25, PDFscale2, iSysNow, &beam);
1848 temp = max(temp, xPDFval);
1851 xPDFmotherSum += temp;
1852 xPDFmother.insert(make_pair(xPDFmotherSum, make_pair(i, temp) ) );
1857 double R = xPDFmotherSum * RN;
1858 if (xPDFmother.lower_bound(R) == xPDFmother.end()) {
1859 idMother = xPDFmother.rbegin()->second.first;
1860 pdfRatio = xPDFmother.rbegin()->second.second / xPDFdaughter;
1862 idMother = xPDFmother.lower_bound(R)->second.first;
1863 pdfRatio = xPDFmother.lower_bound(R)->second.second / xPDFdaughter;
1865 idSister = idMother;
1867 pdfRatio = xPDFmotherSum / xPDFdaughter;
1871 }
else if (splitNow->is(splittingsPtr->isrQCD_1_to_1_and_1_and_1)) {
1873 double xPDFmother = getXPDF( -idDau, xDau, PDFscale2, iSysNow, &beam);
1874 if ( xPDFmother != 0. && abs(xPDFmother) < tinypdf(xDau) ) {
1875 int sign = (xPDFmother >= 0.) ? 1 : -1;
1876 xPDFmother = sign*tinypdf(xDau);
1878 pdfRatio = xPDFmother / xPDFdaughter;
1880 }
else if ( tOld < PT2_INCREASE_OVERESTIMATE
1881 && ( splitNow->is(splittingsPtr->isrQCD_21_to_21_and_21a)
1882 || splitNow->is(splittingsPtr->isrQCD_21_to_21_and_21b))) {
1883 double xPDFmother = xPDFdaughter;
1884 int NTSTEPS(3), NXSTEPS(3);
1885 for (
int i=1; i <= NTSTEPS; ++i) {
1886 double tNew = pT2min + double(i)/double(NTSTEPS)*(max(tOld, pT2min)
1888 for (
int j=1; j <= NXSTEPS; ++j) {
1889 double xNew = xDau + double(j)/double(NXSTEPS)*(0.999999-xDau);
1890 double xPDFnew = getXPDF( 21, xNew, tNew, iSysNow, &beam);
1891 xPDFmother = max(xPDFmother, xPDFnew);
1894 pdfRatio = xPDFmother/xPDFdaughter;
1898 double xPDFmother = getXPDF(idMother,xDau, PDFscale2, iSysNow, &beam);
1900 int NTSTEPS(3), NXSTEPS(3);
1901 for (
int i=0; i <= NTSTEPS; ++i) {
1902 double tNew = pT2min + double(i)/double(NTSTEPS)*(max(tOld, pT2min)
1904 for (
int j=1; j <= NXSTEPS; ++j) {
1905 double xNew = xDau + double(j)/double(NXSTEPS)*(0.999999-xDau);
1906 double xPDFnew = getXPDF(idMother, xNew, tNew, iSysNow, &beam);
1907 xPDFmother = max(xPDFmother, xPDFnew);
1910 pdfRatio = xPDFmother/xPDFdaughter;
1913 if (pdfRatio < 0. || abs(xPDFdaughter) < tinypdf(xDau)) pdfRatio = 0.;
1914 if (!inD) pdfRatio = 0.;
1924 void DireSpace::getNewSplitting(
const Event& state, DireSpaceEnd* dip,
1925 double tOld,
double xDau,
double t,
double zMinAbs,
double zMaxAbs,
1926 int idDau,
string name,
bool forceFixedAs,
int& idMother,
int& idSister,
1927 double& z,
double& wt, unordered_map<string,double>& full,
double& over ) {
1929 BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
1930 bool isValence = (usePDF) ? beam[iSysNow].isValence() :
false;
1932 DireSplitting* splitNow = splits[name];
1934 splitNow->splitInfo.storeRadBef(state[dip->iRadiator]);
1935 splitNow->splitInfo.storeRecBef(state[dip->iRecoiler]);
1936 bool canUseSplitInfo = splitNow->canUseForBranching();
1939 vector<int> re = splitNow->radAndEmt(idDau, dip->colType);
1944 if ( pT2cut(idSister) > t) { wt = over = 0.; full.clear();
return; }
1947 if(z< 0.) z = splitNow->zSplit(zMinAbs, zMaxAbs, dip->m2Dip);
1948 over = splitNow->overestimateDiff(z, dip->m2Dip);
1951 if (!splitNow->aboveCutoff( t, state[dip->iRadiator], state[dip->iRecoiler],
1952 dip->system, partonSystemsPtr)) { wt = over = 0.; full.clear();
return; }
1955 double RNflav = rndmPtr->flat();;
1956 double pdfRatio = getPDFOverestimates(idDau, tOld, xDau, name,
true, RNflav,
1957 idMother, idSister);
1963 int type = (state[dip->iRecoiler].isFinal()) ? 1 : -1;
1965 m2s = particleDataPtr->isResonance(state[dip->iRecoiler].id())
1966 ? getMass(state[dip->iRecoiler].id(),3,
1967 state[dip->iRecoiler].mCalc())
1968 : (state[dip->iRecoiler].idAbs() < 6)
1969 ? getMass(state[dip->iRecoiler].id(),2)
1970 : getMass(state[dip->iRecoiler].id(),1);
1977 if (type == 1 && (m2Bef > TINYMASS || m2r > TINYMASS || m2e > TINYMASS
1978 || m2s > TINYMASS)) type = 2;
1979 if (type ==-1 && (m2Bef > TINYMASS || m2r > TINYMASS || m2e > TINYMASS
1980 || m2s > TINYMASS)) type =-2;
1983 double m2dipCorr = dip->m2Dip - m2Bef + m2r + m2e;
1988 int kinType = splitNow->kinMap();
1994 dip->phi = 2.*M_PI*rndmPtr->flat();
1998 double m2i = getMass(idMother,2);
1999 double m2j = getMass(idSister,2);
2000 bool physical =
true;
2002 if ( splitNow->nEmissions() == 2 ) {
2003 dip->mass.push_back(m2r);
2004 dip->mass.push_back(m2i);
2005 dip->mass.push_back(m2j);
2006 dip->mass.push_back(m2s);
2008 zCollNextQCD( dip, dip->z, 1. );
2010 physical = virtNextQCD( dip, 0.0, dip->m2Dip);
2012 dip->phia1 = 2.*M_PI*rndmPtr->flat();
2016 vector <double> aux;
2017 if ( splitNow->nEmissions() == 2 ) {
2018 type = (state[dip->iRecoiler].isFinal()) ? 2 : -2;
2019 aux.push_back( dip->m2Dip );
2020 if (type > 0) aux.push_back( (state[dip->iRadiator].p()
2021 -state[dip->iRecoiler].p()).m2Calc() );
2022 else aux.push_back( (state[dip->iRadiator].p()
2023 +state[dip->iRecoiler].p()).m2Calc() );
2024 aux.push_back(dip->pT2);
2025 aux.push_back(dip->sa1);
2026 aux.push_back(dip->z);
2027 aux.push_back(dip->xa);
2028 aux.push_back(m2Bef);
2036 int nEmissions = splitNow->nEmissions();
2037 splitNow->splitInfo.storeInfo(name, type, dip->system, dip->system,
2038 dip->side, dip->iRadiator, dip->iRecoiler, state, idSister, idMother,
2039 nEmissions, m2dipCorr, dip->pT2, dip->pT2Old, dip->z, dip->phi, m2Bef,
2040 m2s, m2r, (nEmissions == 1 ? m2e : m2i), dip->sa1, dip->xa, dip->phia1,
2043 if (canUseSplitInfo) {
2044 splitNow->splitInfo.setRadAft(re[0]);
2045 splitNow->splitInfo.setEmtAft(re[1]);
2046 if (nEmissions==2) splitNow->splitInfo.setEmtAft2(re[2]);
2047 splitNow->splitInfo.canUseForBranching(
true);
2049 splitNow->splitInfo.setRadAft(idMother);
2050 splitNow->splitInfo.setEmtAft(idSister);
2051 if (nEmissions==2) splitNow->splitInfo.setEmtAft2(-idSister);
2055 double zcheck(z), tcheck(t);
2058 ? splitNow->zdire_ii(z, t, m2dipCorr)
2059 : splitNow->zdire_if(z, t, m2dipCorr);
2061 ? splitNow->tdire_ii(z, t, m2dipCorr)
2062 : splitNow->tdire_if(z, t, m2dipCorr);
2064 if ( !physical || !inAllowedPhasespace( kinType, zcheck, tcheck, m2dipCorr,
2065 xDau, type, m2Bef, m2r, m2s, m2e, aux ) )
2066 { wt = over = 0.; full.clear();
return; }
2069 int order = kernelOrder;
2071 bool hasInA = (getInA(dip->system) != 0);
2072 bool hasInB = (getInB(dip->system) != 0);
2073 if (dip->system != 0 && hasInA && hasInB) order = kernelOrderMPI;
2076 if (canUseSplitInfo) {
2077 vector< pair<int,int> > cols
2078 = splitNow->radAndEmtCols( dip->iRadiator, dip->colType, state);
2079 splitNow->splitInfo.setRadAft(re[0], cols[0].first, cols[0].second);
2080 splitNow->splitInfo.setEmtAft(re[1], cols[1].first, cols[1].second);
2081 if (nEmissions==2) splitNow->splitInfo.setEmtAft2(re[2], cols[2].first,
2086 over = splitNow->overestimateDiff(z, dip->m2Dip, order);
2089 if (splitNow->calc( state, order) ) full = splitNow->getKernelVals();
2091 if (!dryrun && splitNow->hasMECBef(state, tOld)) over *= KERNEL_HEADROOM;
2092 if (!dryrun && splitNow->hasMECBef(state, dip->pT2))
2093 for (unordered_map<string,double>::iterator it=full.begin();
2094 it != full.end(); ++it) it->second *= KERNEL_HEADROOM;
2096 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
2097 <<
" " << __LINE__ <<
" : New splitting "
2098 << setw(15) << name <<
" at pT="
2099 << setw(15) << sqrt(dip->pT2) <<
" z = "
2100 << setw(15) << dip->z <<
" prob = "
2101 << setw(15) << full[
"base"] << endl;
2104 double coupl = splitNow->coupling(dip->z, dip->pT2, m2dipCorr, -1.,
2105 make_pair(state[dip->iRadiator].id(), state[dip->iRadiator].isFinal()),
2106 make_pair(state[dip->iRecoiler].id(), state[dip->iRecoiler].isFinal()));
2108 double scale2 = splits[splittingNowName]->couplingScale2(dip->z, dip->pT2,
2110 make_pair (state[dip->iRadiator].id(), state[dip->iRadiator].isFinal()),
2111 make_pair (state[dip->iRecoiler].id(), state[dip->iRecoiler].isFinal()));
2112 if (scale2 < 0.) scale2 = dip->pT2;
2113 double talpha = max(scale2, pT2min);
2114 double renormMultFacNow = renormMultFac;
2115 if (forceFixedAs) renormMultFacNow = 1.0;
2118 full[
"base"] *= coupl / alphasNow(talpha, renormMultFacNow, dip->system);
2119 if (name.find(
"qcd") == string::npos) {
2120 for ( unordered_map<string,double>::iterator it = full.begin();
2121 it != full.end(); ++it ) {
2122 if (it->first ==
"base")
continue;
2123 it->second *= coupl / alphasNow(talpha, renormMultFacNow, dip->system);
2128 vector <int> in, out;
2129 for (
int i=0; i < state.size(); ++i) {
2130 if (i == dip->iRadiator)
continue;
2131 if (state[i].isFinal()) out.push_back(state[i].id());
2132 if (state[i].mother1() == 1 && state[i].mother2() == 0)
2133 in.push_back(state[i].id());
2134 if (state[i].mother1() == 2 && state[i].mother2() == 0)
2135 in.push_back(state[i].id());
2137 in.push_back(re[0]);
2138 for (
size_t i=1; i < re.size(); ++i) out.push_back(re[i]);
2139 bool hasME = dip->pT2 > pT2minMECs && doMEcorrections
2140 && weights->hasME(in,out);
2141 if (hasME)
for (unordered_map<string,double>::iterator it=full.begin();
2142 it != full.end(); ++it) it->second = abs(it->second);
2146 for (
int i=0; i < state.size(); ++i)
if (state[i].isFinal()) nFinal++;
2147 if (!dryrun) mecover = splitNow->
2148 overhead(dip->m2Dip*xDau, state[dip->iRadiator].id(), nFinal);
2149 for (unordered_map<string,double>::iterator it=full.begin();
2150 it != full.end(); ++it) it->second *= mecover;
2154 wt = full[
"base"]/over;
2157 if (pdfRatio != 0.) wt /= pdfRatio;
2162 double headRoom = overheadFactors(name, idDau, isValence, dip->m2Dip, tOld);
2173 pair<bool, pair<double, double> > DireSpace::getMEC (
const Event& state,
2174 DireSplitInfo* splitInfo) {
2176 double MECnum(1.0), MECden(1.0);
2179 = weights->hasME(makeHardEvent(max(0,splitInfo->system), state,
false));
2184 mergingHooksPtr->init();
2187 mergingHooksPtr->orderHistories(
false);
2190 if ( mergingHooksPtr->getProcessString().compare(
"pp>h") == 0)
2191 mergingHooksPtr->allowCutOnRecState(
true);
2196 Event newProcess( mergingHooksPtr->bareEvent(
2197 makeHardEvent(max(0,splitInfo->system), state,
false),
false) );
2199 mergingHooksPtr->storeHardProcessCandidates( newProcess );
2202 int nSteps = mergingHooksPtr->
2203 getNumberOfClusteringSteps( newProcess,
true);
2205 newProcess.scale(0.0);
2207 DireHistory myHistory( nSteps, 0.0, newProcess, DireClustering(),
2208 mergingHooksPtr, (*beamAPtr), (*beamBPtr), particleDataPtr, infoPtr,
2209 NULL, splits.begin()->second->fsr, splits.begin()->second->isr, weights,
2210 coupSMPtr,
true,
true, 1.0, 1.0, 1.0, 1.0, 0);
2212 myHistory.projectOntoDesiredHistories();
2214 MECnum = myHistory.MECnum;
2215 MECden = myHistory.MECden;
2218 mergingHooksPtr->init();
2222 if (abs(MECden) < 1e-15) direInfoPtr->message(1) << __FILE__ <<
" "
2224 <<
" " << __LINE__ <<
" : Small MEC denominator="
2225 << MECden <<
" for numerator=" << MECnum << endl;
2226 if (abs(MECnum/MECden) > 1e2) {direInfoPtr->message(1) << __FILE__ <<
" "
2228 <<
" " << __LINE__ <<
" : Large MEC. Denominator="
2229 << MECden <<
" Numerator=" << MECnum <<
" at pT="
2230 << sqrt(splitInfo->kinematics()->pT2) <<
" "
2234 return make_pair(hasME, make_pair(MECnum,MECden));
2240 bool DireSpace::applyMEC (
const Event& state, DireSplitInfo* splitInfo,
2241 vector<Event> auxState) {
2244 pair<bool, pair<double, double> > mec = getMEC ( state, splitInfo);
2245 bool hasME = mec.first;
2246 double MECnum = mec.second.first;
2247 double MECden = mec.second.second;
2248 double MECnumX = mec.second.first;
2249 double MECdenX = mec.second.second;
2251 if (!hasME)
return false;
2253 double kernel = kernelSel[
"base"];
2254 bool reject =
false;
2256 if (abs(MECnum/MECden) > 5e0 && auxState.size()>0) {
2257 pair<bool, pair<double, double> > mec1 = getMEC ( auxState[0], splitInfo);
2258 pair<bool, pair<double, double> > mec2 = getMEC ( auxState[1], splitInfo);
2259 double MECnum1 = mec1.second.first;
2260 double MECden1 = mec1.second.second;
2261 double MECnum2 = mec2.second.first;
2262 double MECden2 = mec2.second.second;
2263 if (MECnum/MECden > MECnum1/MECden1) {MECnum = MECnum1; MECden = MECden1;}
2264 if (MECnum/MECden > MECnum2/MECden2) {MECnum = MECnum2; MECden = MECden2;}
2265 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
2266 <<
" " << __LINE__ <<
" : Large MEC weight=" << MECnumX/MECdenX
2267 <<
" " << MECnum/MECden
2268 <<
"\t\t" << splitInfo->kinematics()->pT2/splitInfo->kinematics()->m2Dip
2269 <<
" " << splitInfo->kinematics()->z << endl;
2270 if (MECnum/MECden > (MECnum+MECnum1)/(MECden+MECden1))
2271 { MECnum += MECnum1; MECden += MECden1; }
2272 if (MECnum/MECden > (MECnum+MECnum2)/(MECden+MECden2))
2273 { MECnum += MECnum2; MECden += MECden2; }
2278 if (kernelSel.find(
"base_order_as2") != kernelSel.end() ) {
2279 oas2 = kernelSel[
"base_order_as2"];
2280 kernelSel.erase(kernelSel.find(
"base_order_as2"));
2282 double baseNew = ((kernel - oas2) * MECnum/MECden + oas2);
2285 double auxNew = kernel;
2286 double overNew = kernel;
2289 for (
int i=0; i < state.size(); ++i)
2290 if (state[i].isFinal()) nFinal++;
2292 if (dryrun) splittingSel->storeOverhead(
2293 splitInfo->kinematics()->m2Dip*splitInfo->kinematics()->xBef,
2294 splitInfo->kinematics()->xBef, state[splitInfo->iRadBef].id(), nFinal-1,
2295 max(baseNew/overNew,1.1));
2298 if (baseNew/auxNew < 0.) auxNew *= -1.;
2300 if (suppressLargeMECs)
while (baseNew/auxNew < 5e-2) auxNew /= 5.;
2303 if (baseNew/auxNew > 1.) {
2304 double rescale = baseNew/auxNew * 1.15;
2307 double wt = baseNew/auxNew;
2310 double wvNow = auxNew/overNew * (overNew - baseNew) / (auxNew - baseNew);
2313 double waNow = auxNew/overNew;
2315 if (abs(wvNow) > 1e0) {
2316 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
2317 <<
" " << __LINE__ <<
" : Large reject weight=" << wvNow
2318 <<
"\t for kernel=" << baseNew <<
" overestimate=" << overNew
2319 <<
"\t aux. overestimate=" << auxNew <<
" at pT2="
2320 << splitInfo->kinematics()->pT2
2321 <<
" for " << splittingSelName
2324 if (abs(waNow) > 1e0) {
2325 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
2326 <<
" " << __LINE__ <<
" : Large accept weight=" << waNow
2327 <<
"\t for kernel=" << baseNew <<
" overestimate=" << overNew
2328 <<
"\t aux. overestimate=" << auxNew <<
" at pT2="
2329 << splitInfo->kinematics()->pT2
2330 <<
" for " << splittingSelName
2334 if (wt < rndmPtr->flat()) {
2337 for (unordered_map<string,double>::iterator it= kernelSel.begin();
2338 it != kernelSel.end(); ++it) {
2340 double waOld = weights->getAcceptWeight( splitInfo->kinematics()->pT2,
2343 weights->eraseAcceptWeight(splitInfo->kinematics()->pT2, it->first);
2344 weights->resetRejectWeight(splitInfo->kinematics()->pT2, wvNow*waOld,
2350 for (unordered_map<string,double>::iterator it = kernelSel.begin();
2351 it != kernelSel.end(); ++it) {
2353 double waOld = weights->getAcceptWeight( splitInfo->kinematics()->pT2,
2356 weights->eraseRejectWeight(splitInfo->kinematics()->pT2, it->first);
2357 weights->resetAcceptWeight(splitInfo->kinematics()->pT2, waOld*waNow,
2375 bool DireSpace::inAllowedPhasespace(
int kinType,
double z,
double pT2,
2376 double m2dip,
double xOld,
int splitType,
double m2RadBef,
double m2r,
2377 double m2s,
double m2e, vector<double> aux) {
2379 double xIncoming = usePDF ? xOld : 0.;
2382 if (splitType == 1) {
2385 double kappa2 = pT2 / m2dip;
2387 double uCS = kappa2/(1-z);
2392 uCS = 0.5*xCS*( 1. - sqrt( 1. - 4.*xCS*kappa2 / pow2( 1- xCS)) );
2396 if ( xCS < xIncoming || xCS > 1. || uCS < 0. || uCS > 1. )
return false;
2399 }
else if (splitType == 2 && aux.size() == 0) {
2401 double kappa2 = pT2 / m2dip;
2403 double uCS = kappa2/(1-z);
2405 double pijpa = m2dip - m2r - m2e + m2RadBef;
2406 double mu2Rec = m2s / pijpa * xCS;
2407 double uCSmax = (1. - xCS) / (1. - xCS + mu2Rec );
2410 if (xCS < xIncoming || xCS > 1. || uCS < 0. || uCS > uCSmax)
return false;
2413 }
else if (splitType == 2 && aux.size() > 0) {
2416 if (
int(aux.size()) < 11)
return false;
2420 double sai = aux[3];
2423 double m2a = aux[7];
2424 double m2i = aux[8];
2425 double m2j = aux[9];
2426 double m2k = aux[10];
2427 double m2ai = -sai + m2a + m2i;
2430 double uCS = za*(m2ai-m2a-m2i) / q2;
2431 double xCS = uCS + xa - (t*za) / (q2*xa);
2432 double m2jk = t/xa + q2*( 1. - xa/za ) - m2ai;
2434 if ( m2jk < 0. )
return false;
2436 double mu2Rec = m2jk/(-q2+m2jk) * xCS;
2437 double uCSmax = (1. - xCS) / (1. - xCS + mu2Rec );
2440 if (xCS < xIncoming || xCS > 1. || uCS < 0. || uCS > uCSmax)
return false;
2443 double s_i_jk = (1. - 1./xCS)*(q2 - m2a) + (m2i + m2jk) / xCS;
2444 double zbar = (q2-s_i_jk-m2a) / bABC(q2,s_i_jk,m2a)
2445 *( uCS - m2a / gABC(q2,s_i_jk,m2a)
2446 * (s_i_jk + m2i - m2jk) / (q2 - s_i_jk - m2a));
2447 double kT2 = zbar*(1.-zbar)*s_i_jk - (1-zbar)*m2i - zbar*m2jk;
2448 if (kT2 < 0.)
return false;
2451 double zCS = t/xa / ( t/xa - q2*xa/za);
2452 double yCS = (m2jk - m2k - m2j)
2453 / (m2jk - m2k - m2j + t/xa - q2*xa/za);
2455 double q2_2 = m2ai + m2jk + t/xa - q2*xa/za;
2457 double sij = yCS * (q2_2 - m2ai) + (1.-yCS)*(m2j+m2k);
2458 zbar = (q2_2-sij-m2ai) / bABC(q2_2,sij,m2ai)
2459 * (zCS - m2ai/gABC(q2_2,sij,m2ai)
2460 *(sij + m2j - m2k)/(q2_2-sij-m2ai));
2461 kT2 = zbar*(1.-zbar)*sij - (1.-zbar)*m2j - zbar*m2k;
2463 if (kT2 < 0.)
return false;
2466 }
else if (splitType == -1) {
2469 double kappa2 = pT2 / m2dip;
2470 double xCS = (z*(1-z)- kappa2)/(1-z);
2471 double vCS = kappa2/(1-z);
2476 vCS = 0.5*xCS*( 1. - sqrt( 1. - 4.*xCS*kappa2 / pow2( 1- xCS)) );
2480 if (xCS < xIncoming || xCS > 1. || vCS < 0. || vCS > 1.)
return false;
2481 if (1.-xCS-vCS < 0.)
return false;
2484 }
else if (splitType == -2 && aux.size() == 0) {
2488 double q2 = m2dip + m2s + m2RadBef;
2489 double m2DipCorr = m2dip - m2RadBef + m2r + m2e;
2490 double kappa2 = pT2 / m2DipCorr;
2491 double xCS = (z*(1-z)- kappa2)/(1-z);
2492 double vCS = kappa2/(1-z);
2495 double sab = (q2 - m2e)/xCS + (m2r+m2s) * (1-1/xCS);
2496 double saj = -vCS*(sab - m2r-m2s) + m2r + m2e;
2497 double zbar = (sab - m2r - m2s) / bABC(sab,m2r,m2s)
2498 *( (xCS + vCS) - m2s / gABC(sab,m2r,m2s)
2499 * (saj + m2r - m2e) / (sab - m2r - m2s));
2500 double kT2 = zbar*(1.-zbar)*m2r - (1-zbar)*saj - zbar*m2e;
2503 if (kT2 < 0. || isnan(kT2))
return false;
2511 if (
int(aux.size()) < 11)
return false;
2514 double q2_1 = aux[1];
2516 double sai = aux[3];
2519 double m2a = aux[7];
2520 double m2i = aux[8];
2521 double m2j = aux[9];
2522 double m2k = aux[10];
2523 double m2ai = -sai + m2a + m2i;
2525 if (za < xIncoming || za > 1.)
return false;
2528 double p2ab = q2_1/za + m2a + m2k;
2529 double zbar = (p2ab - m2a - m2k) / bABC(p2ab,m2a,m2k)
2530 *( xa - m2k / gABC(p2ab,m2a,m2k)
2531 * (m2ai + m2a - m2i) / (p2ab - m2a - m2k));
2532 double kT2 = zbar*(1.-zbar)*m2a - (1-zbar)*m2ai - zbar*m2i;
2535 if (kT2 < 0. || isnan(kT2))
return false;
2538 double m2rec = m2ai;
2539 double m2emt = q2_1;
2541 double zCS = t/xa / (q2_1*xa/za + 2.*m2ai);
2542 double yCS = 1. / ( 1. + (q2_1*xa/za + 2.*m2ai)
2543 / (q2_1*(xa/za - 1.) + m2ai + m2k - m2j));
2544 double q2_2 = 4.*m2ai + 2.*q2_1*xa/za + m2k;
2547 if (yCS < 0. || yCS > 1. || zCS < 0. || zCS > 1.)
return false;
2550 double sij = yCS * (q2_2 - m2rec) + (1.-yCS)*(m2rad+m2emt);
2551 zbar = (q2_2-sij-m2rec) / bABC(q2_2,sij,m2rec)
2552 * (zCS - m2rec/gABC(q2_2,sij,m2rec)
2553 *(sij + m2rad - m2emt)/(q2_2-sij-m2rec));
2554 kT2 = zbar*(1.-zbar)*sij - (1.-zbar)*m2rad - zbar*m2emt;
2557 if (kT2 < 0. || isnan(kT2))
return false;
2569 void DireSpace::addNewOverestimates( multimap<double,string> newOverestimates,
2570 double& oldOverestimate ) {
2573 if (!newOverestimates.empty())
2574 oldOverestimate += newOverestimates.rbegin()->first;
2584 void DireSpace::alphasReweight(
double,
double talpha,
int iSys,
2585 bool forceFixedAs,
double& weight,
double& fullWeight,
double& overWeight,
2586 double renormMultFacNow) {
2589 overWeight *= alphaS2piOverestimate;
2590 weight *= alphasNow(pT2min, 1., iSys) / alphaS2piOverestimate;
2591 fullWeight *= alphasNow(pT2min, 1., iSys);
2594 talpha = max(talpha, pT2min);
2596 double scale = talpha*renormMultFacNow;
2597 scale = max(scale, pT2min);
2600 double asPT2piCorr = alphasNow(talpha, renormMultFacNow, iSys);
2604 if (usePDFalphas) asOver = alphaS2piOverestimate;
2605 else if (alphaSorder==0) asOver = alphaS2pi;
2606 else asOver = asPT2piCorr;
2609 if (alphaSorder == 0) asFull = alphaS2pi;
2610 else asFull = asPT2piCorr;
2612 fullWeight *= asFull;
2613 overWeight *= asOver;
2614 weight *= asFull/asOver;
2624 void DireSpace::pT2nextQCD(
double pT2begDip,
double pT2endDip,
2625 DireSpaceEnd& dip,
Event& event,
double pT2endForce,
double pT2freeze,
2626 bool forceBranching) {
2628 if (event[dip.iRecoiler].isFinal()) {
2629 pT2nextQCD_IF(pT2begDip, pT2endDip, dip, event, pT2endForce, pT2freeze,
2632 pT2nextQCD_II(pT2begDip, pT2endDip, dip, event, pT2endForce, pT2freeze,
2643 bool DireSpace::pT2nextQCD_II(
double pT2begDip,
double pT2sel,
2644 DireSpaceEnd& dip,
Event& event,
double pT2endForce,
double pT2freeze,
2645 bool forceBranching) {
2648 double pT2endDip = max( pT2sel, pT2cutMin(&dip));
2649 if (pT2endForce >= 0.) pT2endDip = pT2endForce;
2650 if (pT2begDip < pT2endDip)
return false;
2653 int iRadi = dip.iRadiator;
2654 int iReco = dip.iRecoiler;
2655 m2Dip = abs(2.*event[iRadi].p()*event[iReco].p());
2659 BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
2660 double tnow = pT2begDip;
2661 double xMaxAbs = beam.xMax(iSysNow);
2662 double zMinAbs = xDaughter;
2664 if (usePDF && xMaxAbs < 0.) {
2665 infoPtr->errorMsg(
"Warning in DireSpace::pT2nextQCD_II: "
2666 "xMaxAbs negative");
2672 double Lambda2 = Lambda3flav2;
2676 double zMaxAbs = 0.;
2677 double xPDFdaughter = 0.;
2678 double kernelPDF = 0.;
2679 double xMother = 0.;
2681 double mSister = 0.;
2682 double m2Sister = 0.;
2683 double pT2corr = 0.;
2684 double teval = pT2begDip;
2685 bool needNewPDF =
true;
2686 bool hasPDFdau = hasPDF(idDaughter);
2687 if (!hasPDFdau) zMinAbs = 0.;
2689 multimap<double,string> newOverestimates;
2690 unordered_map<string,double> fullWeightsNow;
2691 double fullWeightNow(0.), overWeightNow(0.), auxWeightNow(0.), daux(0.);
2694 int loopTinyPDFdau = 0;
2695 int nContinue(0), nContinueMax(10000);
2696 bool hasTinyPDFdau =
false;
2701 tnow = (!forceBranching) ? tnow : pT2begDip;
2707 if (forceBranching && nContinue >= nContinueMax) {
2708 wt = 0.0; dip.pT2 = tnow = 0.;
2713 if ( fullWeightNow != 0. && overWeightNow != 0. ) {
2714 double enhanceFurther
2715 = enhanceOverestimateFurther(splittingNowName, idDaughter, teval);
2716 if (doTrialNow) enhanceFurther = 1.;
2717 kernelNow = fullWeightsNow;
2718 auxNow = auxWeightNow;
2719 overNow = overWeightNow;
2720 boostNow = enhanceFurther;
2722 for ( unordered_map<string,double>::iterator it = fullWeightsNow.begin();
2723 it != fullWeightsNow.end(); ++it ) {
2726 if (it->first ==
"base_order_as2")
continue;
2728 double wv = auxWeightNow/overWeightNow
2729 * (overWeightNow- it->second/enhanceFurther)
2730 / (auxWeightNow - fullWeightNow);
2731 if (abs(wv) > 1e0) {
2732 direInfoPtr->message(1)
2733 << scientific << setprecision(15)
2734 << __FILE__ <<
" " << __func__
2735 <<
" " << __LINE__ <<
" : Large reject weight=" << wv
2736 <<
"\t for kernel=" << it->second
2737 <<
" " << fullWeightNow <<
" overestimate=" << overNow
2738 <<
"\t aux. overestimate=" << auxNow <<
" at pT2="
2740 <<
" for " << splittingNowName << endl;
2742 rejectProbability[it->first].insert( make_pair(tnow,wv));
2746 splittingNowName=
"";
2747 fullWeightsNow.clear();
2748 fullWeightNow = overWeightNow = auxWeightNow = 0.;
2751 if (abs(idDaughter)==4 && tnow <= m2cPhys) { dip.pT2 = 0.0;
return false;}
2752 if (abs(idDaughter)==5 && tnow <= m2bPhys) { dip.pT2 = 0.0;
return false;}
2755 double tnew = (useFixedFacScale) ? fixedFacScale2 : factorMultFac*tnow;
2756 tnew = max(tnew, pT2min);
2757 bool inNew = (hasPDFdau) ? beam.insideBounds(xDaughter, tnew) :
true;
2758 if (hasPDFdau && !inNew) { dip.pT2 = 0.0;
return false; }
2762 if (hasTinyPDFdau) ++loopTinyPDFdau;
2763 if (hasPDFdau && loopTinyPDFdau > MAXLOOPTINYPDF) {
2764 infoPtr->errorMsg(
"Warning in DireSpace::pT2nextQCD: "
2765 "small daughter PDF");
2774 || tnow < evalpdfstep(event[iRadi].
id(), tnow, m2cPhys, m2bPhys)*teval) {
2777 hasTinyPDFdau =
false;
2779 newOverestimates.clear();
2785 Lambda2 = Lambda5flav2;
2786 }
else if (tnow > m2c) {
2788 Lambda2 = Lambda4flav2;
2791 Lambda2 = Lambda3flav2;
2795 Lambda2 /= renormMultFac;
2796 zMinAbs = (hasPDFdau) ? xDaughter : 0.;
2800 pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac*tnow;
2801 pdfScale2 = max(pdfScale2, pT2min);
2802 xPDFdaughter = getXPDF(idDaughter, xDaughter, pdfScale2, iSysNow, &beam);
2803 if ( hasPDFdau && xPDFdaughter != 0.
2804 && abs(xPDFdaughter) < tinypdf(xDaughter)) {
2805 int sign = (xPDFdaughter > 0.) ? 1 : -1;
2806 xPDFdaughter = sign*tinypdf(xDaughter);
2807 hasTinyPDFdau =
true;
2811 getNewOverestimates( idDaughter, &dip, event, teval,
2812 xDaughter, zMinAbs, zMaxAbs, newOverestimates);
2813 addNewOverestimates(newOverestimates, kernelPDF);
2822 if (kernelPDF < TINYKERNELPDF) { dip.pT2 = 0.0;
return false; }
2823 if (newOverestimates.empty()) { dip.pT2 = 0.0;
return false; }
2826 bool forceFixedAs = (tnow < pT2min);
2827 tnow = tNextQCD( &dip, kernelPDF, tnow, pT2endDip, pT2freeze,
2828 (forceBranching ? -1 : 1));
2830 wt = dip.pT2 = tnow = 0.;
2831 double R0 = kernelPDF*rndmPtr->flat();
2832 if (!newOverestimates.empty()) {
2833 if (newOverestimates.lower_bound(R0) == newOverestimates.end())
2834 splittingNowName = newOverestimates.rbegin()->second;
2836 splittingNowName = newOverestimates.lower_bound(R0)->second;
2842 if ( tnow <= pT2endDip) { dip.pT2 = tnow = 0.;
break; }
2846 if (nFlavour == 5 && tnow <= m2bPhys) {
2849 }
else if (nFlavour == 4 && tnow <= m2cPhys) {
2854 if (abs(idDaughter) == 5 && tnow <= m2bPhys) {
2855 dip.pT2 = tnow = 0.;
break;
2856 }
else if (abs(idDaughter) == 4 && tnow <= m2cPhys) {
2857 dip.pT2 = tnow = 0.;
break;
2860 double xMin = xDaughter;
2861 if (!hasPDFdau) xMin = 0.;
2865 double R = kernelPDF*rndmPtr->flat();
2866 if (!newOverestimates.empty()) {
2867 if (newOverestimates.lower_bound(R) == newOverestimates.end())
2868 splittingNowName = newOverestimates.rbegin()->second;
2870 splittingNowName = newOverestimates.lower_bound(R)->second;
2871 getNewSplitting( event, &dip, teval, xMin, tnow, zMinAbs,
2872 zMaxAbs, idDaughter, splittingNowName, forceFixedAs, idMother,
2873 idSister, znow, wt, fullWeightsNow, overWeightNow);
2878 fullWeightsNow.clear();
2879 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
2880 nContinue++;
continue;
2883 fullWeightNow = fullWeightsNow[
"base"];
2887 if ( tnow <= 4.*m2bPhys
2888 && ( (abs(idDaughter) == 21 && abs(idSister) == 5)
2889 || (abs(idDaughter) == 5 && splits[splittingNowName]->nEmissions()==2)
2890 || (abs(idSister) == 5 && splits[splittingNowName]->nEmissions()==2))) {
2891 fullWeightsNow.clear();
2892 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
2893 nContinue++;
continue;
2894 }
else if ( tnow <= 4.*m2cPhys
2895 && ( (abs(idDaughter) == 21 && abs(idSister) == 4)
2896 || (abs(idDaughter) == 4 && splits[splittingNowName]->nEmissions()==2)
2897 || (abs(idSister) == 4 && splits[splittingNowName]->nEmissions()==2))) {
2898 fullWeightsNow.clear();
2899 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
2900 nContinue++;
continue;
2909 double q2 = (
event[iRadi].p()+
event[iReco].p()).m2Calc();
2915 if ( splits[splittingNowName]->nEmissions() == 2 )
2916 if ( (abs(idSister) == 4 && tnow < m2cPhys)
2917 || (abs(idSister) == 5 && tnow < m2bPhys)) {
2919 fullWeightsNow.clear();
2920 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
2921 nContinue++;
continue;
2925 double m2a(m2r), m2i(getMass(idMother,2)),
2926 m2j(getMass(-event[iRadi].
id(),2)), m2aij(m2Bef), m2k(m2s);
2929 double m2DipCorr = m2Dip - m2Bef + m2r + m2e;
2931 double kappa2 = tnow / m2DipCorr;
2932 double xCS = (znow*(1-znow)- kappa2)/(1-znow);
2935 double jacobian(1.);
2937 bool canUseSplitInfo = splits[splittingNowName]->canUseForBranching();
2938 if (canUseSplitInfo) {
2940 = splits[splittingNowName]->getJacobian(event,partonSystemsPtr);
2941 unordered_map<string,double> psvars
2942 = splits[splittingNowName]->
2943 getPhasespaceVars( event, partonSystemsPtr);
2944 xMother = psvars[
"xInAft"];
2946 if ( splits[splittingNowName]->nEmissions() == 2 ) {
2949 xCS = za * (q2 - m2a - m2i - m2j - m2k) / q2;
2952 double sab = q2/za + m2a + m2k;
2953 jacobian = (sab-m2a-m2k) / sqrt(lABC(sab, m2a, m2k) );
2954 double sai = dip.sa1;
2955 double m2ai = -sai + m2a + m2i;
2956 double sjq = q2*xa/za + m2ai + m2k;
2957 jacobian *= (sjq-m2ai-m2k) / sqrt(lABC(sjq, m2ai, m2k) );
2959 jacobian *= 1. / (1. - (m2ai + m2j - m2aij) / (dip.pT2/dip.xa)) ;
2961 xMother = xDaughter/xCS;
2965 pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * tnow;
2966 pdfScale2 = max(pdfScale2, tnow);
2967 double pdfScale2Old = pdfScale2;
2968 double pdfScale2New = pdfScale2;
2969 if (forceBranching) pdfScale2Old = pdfScale2New = infoPtr->Q2Fac();
2970 bool inD = (hasPDFdau) ? beam.insideBounds(xDaughter, pdfScale2Old) :
true;
2971 bool inM = (hasPDFdau) ? beam.insideBounds(xMother, pdfScale2New) :
true;
2972 double xPDFdaughterNew = getXPDF( idDaughter, xDaughter, pdfScale2Old,
2973 iSysNow, &beam,
false, znow, m2Dip);
2974 double xPDFmotherNew = getXPDF( idMother, xMother, pdfScale2New,
2975 iSysNow, &beam,
false, znow, m2Dip);
2976 if ( hasPDFdau && xPDFdaughterNew != 0.
2977 && abs(xPDFdaughterNew) < tinypdf(xDaughter) ) {
2978 hasTinyPDFdau =
true;
2980 fullWeightsNow.clear();
2981 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
2982 nContinue++;
continue;
2990 double xPDFdaughterLow = getXPDF( idDaughter, xDaughter,
2991 pdfScale2Old*pdfScale2Old/max(teval,pT2min), iSysNow, &beam);
2992 if ( hasPDFdau && idDaughter == 21
2993 && ( abs(xPDFdaughterNew/xPDFdaughter) < 1e-4
2994 || abs(xPDFdaughterLow/xPDFdaughterNew) < 1e-4) ) {
2995 hasTinyPDFdau =
true;
2997 fullWeightsNow.clear();
2998 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
2999 nContinue++;
continue;
3003 double pdfRatio = (inD && inM) ? xPDFmotherNew/xPDFdaughterNew : 0.;
3006 if (hasPDFdau && idDaughter == 21 && pdfScale2 == pT2min && pdfRatio>50.)
3009 fullWeightNow *= pdfRatio*jacobian;
3011 for ( unordered_map<string,double>::iterator it = fullWeightsNow.begin();
3012 it != fullWeightsNow.end(); ++it )
3013 it->second *= pdfRatio*jacobian;
3017 if ( splits[splittingNowName]->nEmissions() == 2 )
3018 dip.sa1 = splits[splittingNowName]->splitInfo.kinematics()->sai;
3020 if ( fullWeightNow == 0. ) {
3022 fullWeightsNow.clear();
3023 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3024 nContinue++;
continue;
3028 double scale2 = splits[splittingNowName]->couplingScale2 ( dip.z, tnow,
3030 make_pair (event[dip.iRadiator].id(),
event[dip.iRadiator].isFinal()),
3031 make_pair (event[dip.iRecoiler].id(),
event[dip.iRecoiler].isFinal()));
3032 if (scale2 < 0.) scale2 = tnow;
3033 double talpha = max(scale2, pT2min);
3037 alphasReweight(tnow, talpha, dip.system, forceFixedAs, wt, fullWeightNow,
3038 overWeightNow, renormMultFac);
3042 alphasReweight(tnow, talpha, dip.system, forceFixedAs, daux, asw, daux,
3044 fullWeightsNow[
"base"] *= asw;
3045 if (fullWeightsNow.find(
"base_order_as2") != fullWeightsNow.end())
3046 fullWeightsNow[
"base_order_as2"] *= asw;
3048 if ( splittingNowName.find(
"qcd") != string::npos
3049 && settingsPtr->parm(
"Variations:muRisrDown") != 1.) {
3051 alphasReweight(tnow, talpha, dip.system, forceFixedAs, daux, asw, daux,
3052 (tnow > pT2minVariations) ? settingsPtr->parm
3053 (
"Variations:muRisrDown")*renormMultFac
3055 fullWeightsNow[
"Variations:muRisrDown"] *= asw;
3056 }
else if ( splittingNowName.find(
"qcd") == string::npos )
3057 fullWeightsNow[
"Variations:muRisrDown"] *= asw;
3058 if ( splittingNowName.find(
"qcd") != string::npos
3059 && settingsPtr->parm(
"Variations:muRisrUp") != 1.) {
3061 alphasReweight(tnow, talpha, dip.system, forceFixedAs, daux, asw, daux,
3062 (tnow > pT2minVariations) ? settingsPtr->parm
3063 (
"Variations:muRisrUp")*renormMultFac
3065 fullWeightsNow[
"Variations:muRisrUp"] *= asw;
3066 }
else if ( splittingNowName.find(
"qcd") == string::npos )
3067 fullWeightsNow[
"Variations:muRisrUp"] *= asw;
3070 if (hasPDFdau && settingsPtr->flag(
"Variations:PDFup") ) {
3071 int valSea = (beam[iSysNow].isValence()) ? 1 : 0;
3072 if( beam[iSysNow].isUnmatched() ) valSea = 2;
3073 beam.calcPDFEnvelope( make_pair(idMother, idDaughter),
3074 make_pair(xMother, xDaughter), pdfScale2, valSea);
3075 PDF::PDFEnvelope ratioPDFEnv = beam.getPDFEnvelope();
3077 = min(ratioPDFEnv.errplusPDF / ratioPDFEnv.centralPDF, 10.);
3078 double deltaPDFminus
3079 = min(ratioPDFEnv.errminusPDF / ratioPDFEnv.centralPDF, 10.);
3080 fullWeightsNow[
"Variations:PDFup"] = fullWeightsNow[
"base"]
3081 * ((tnow > pT2minVariations) ? (1.0 + deltaPDFplus) : 1.0);
3082 fullWeightsNow[
"Variations:PDFdown"] = fullWeightsNow[
"base"]
3083 * ((tnow > pT2minVariations) ? (1.0 - deltaPDFminus) : 1.0);
3088 auxWeightNow = overWeightNow;
3091 if (fullWeightNow < 0.) {
3092 auxWeightNow *= -1.;
3096 if ( fullWeightNow/auxWeightNow > 1.) {
3097 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
3098 <<
" " << __LINE__ <<
" : Large acceptance weight="
3099 << fullWeightNow/auxWeightNow
3100 <<
" for splitting " << splittingNowName <<
" at pT2=" << tnow
3101 <<
" and z=" << znow <<
"\t(PDF ratio=" << pdfRatio <<
")" << endl;
3102 double rescale = fullWeightNow/auxWeightNow * 1.15;
3103 auxWeightNow *= rescale;
3104 infoPtr->errorMsg(
"Info in DireSpace::pT2nextQCD_II: Found large "
3105 "acceptance weight for " + splittingNowName);
3108 wt = fullWeightNow/auxWeightNow;
3111 }
while (wt < rndmPtr->flat()) ;
3114 if ( wt == 0.) { dip.pT2 = 0.0;
return false; }
3118 if ( fullWeightNow != 0. && overWeightNow != 0. ) {
3119 double enhanceFurther
3120 = enhanceOverestimateFurther(splittingNowName, idDaughter, teval);
3122 weights->addTrialEnhancement(tnow, enhanceFurther);
3123 enhanceFurther = 1.;
3125 kernelNow = fullWeightsNow;
3126 auxNow = auxWeightNow;
3127 overNow = overWeightNow;
3128 boostNow = enhanceFurther;
3129 for ( unordered_map<string,double>::iterator it = fullWeightsNow.begin();
3130 it != fullWeightsNow.end(); ++it ) {
3133 if (it->first ==
"base_order_as2")
continue;
3135 acceptProbability[it->first].insert(make_pair(tnow,
3136 auxWeightNow/overWeightNow * 1./enhanceFurther
3137 * it->second/fullWeightNow ) );
3138 if (auxWeightNow == fullWeightNow && overWeightNow == fullWeightNow)
3139 rejectProbability[it->first].insert( make_pair(tnow, 1.0));
3141 double wv = auxWeightNow/overWeightNow
3142 * (overWeightNow- it->second/enhanceFurther)
3143 / (auxWeightNow - fullWeightNow);
3144 if (abs(wv) > 1e0) {
3145 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
3146 <<
" " << __LINE__ <<
" : Large reject weight=" << wv
3147 <<
"\t for kernel=" << it->second <<
" overestimate=" << overNow
3148 <<
"\t aux. overestimate=" << auxNow <<
" at pT2="
3150 <<
" for " << splittingNowName << endl;
3152 rejectProbability[it->first].insert( make_pair(tnow, wv));
3157 double zStore = dip.z;
3158 double xaStore = dip.xa;
3159 double pT2Store = dip.pT2;
3160 double sa1Store = dip.sa1;
3161 double Q2store = 0.;
3164 dipEndNow->store( idDaughter,idMother, idSister,
3165 x1Now, x2Now, m2Dip, pT2Store, zStore, sa1Store, xaStore, xMother,
3166 Q2store, mSister, m2Sister, pT2corr, dip.phi, dip.phia1);
3177 bool DireSpace::pT2nextQCD_IF(
double pT2begDip,
double pT2sel,
3178 DireSpaceEnd& dip,
Event& event,
double pT2endForce,
double pT2freeze,
3179 bool forceBranching) {
3182 double pT2endDip = max( pT2sel, pT2cutMin(&dip));
3183 if (pT2endForce >= 0.) pT2endDip = pT2endForce;
3184 if (pT2begDip < pT2endDip)
return false;
3187 int iRadi = dip.iRadiator;
3188 int iReco = dip.iRecoiler;
3189 m2Dip = abs(2.*event[iRadi].p()*event[iReco].p());
3193 BeamParticle& beam = (sideA && particleDataPtr->isHadron(beamAPtr->id()))
3195 : (particleDataPtr->isHadron(beamBPtr->id()) ? *beamBPtr
3197 double tnow = pT2begDip;
3198 double xMaxAbs = beam.xMax(iSysNow);
3199 double zMinAbs = xDaughter;
3203 int iOther = sideA ? getInB(iSysNow) : getInA(iSysNow);
3204 Vec4 pOther(event[iOther].p());
3206 if (usePDF && xMaxAbs < 0.) {
3207 infoPtr->errorMsg(
"Warning in DireSpace::pT2nextQCD_IF: "
3208 "xMaxAbs negative");
3214 double Lambda2 = Lambda3flav2;
3218 double zMaxAbs = 0.;
3219 double xPDFdaughter = 0.;
3220 double kernelPDF = 0.;
3221 double xMother = 0.;
3223 double mSister = 0.;
3224 double m2Sister = 0.;
3225 double pT2corr = 0.;
3226 double teval = pT2begDip;
3227 bool needNewPDF =
true;
3228 bool hasPDFdau = hasPDF(idDaughter);
3229 if (!hasPDFdau) zMinAbs = 0.;
3231 multimap<double,string> newOverestimates;
3232 unordered_map<string,double> fullWeightsNow;
3233 double fullWeightNow(0.), overWeightNow(0.), auxWeightNow(0.), daux(0.);
3236 int loopTinyPDFdau = 0;
3237 int nContinue(0), nContinueMax(10000);
3238 bool hasTinyPDFdau =
false;
3245 tnow = (!forceBranching) ? tnow : pT2begDip;
3248 if (forceBranching && nContinue >= nContinueMax) {
3249 wt = 0.0; dip.pT2 = tnow = 0.;
3254 if ( fullWeightNow != 0. && overWeightNow != 0. ) {
3255 double enhanceFurther
3256 = enhanceOverestimateFurther(splittingNowName, idDaughter, teval);
3257 if (doTrialNow) enhanceFurther = 1.;
3258 kernelNow = fullWeightsNow;
3259 auxNow = auxWeightNow;
3260 overNow = overWeightNow;
3261 boostNow = enhanceFurther;
3262 for ( unordered_map<string,double>::iterator it = fullWeightsNow.begin();
3263 it != fullWeightsNow.end(); ++it ) {
3266 if (it->first ==
"base_order_as2")
continue;
3268 double wv = auxWeightNow/overWeightNow
3269 * (overWeightNow- it->second/enhanceFurther)
3270 / (auxWeightNow - fullWeightNow);
3271 if (abs(wv) > 1e0) {
3272 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
3273 <<
" " << __LINE__ <<
" : Large reject weight=" << wv
3274 <<
"\t for kernel=" << it->second <<
" overestimate=" << overNow
3275 <<
"\t aux. overestimate=" << auxNow <<
" at pT2="
3277 <<
" for " << splittingNowName << endl;
3279 rejectProbability[it->first].insert( make_pair(tnow,wv));
3283 splittingNowName=
"";
3284 fullWeightsNow.clear();
3285 fullWeightNow = overWeightNow = auxWeightNow = 0.;
3288 if (abs(idDaughter)==4 && tnow <= m2cPhys) { dip.pT2 = 0.0;
return false;}
3289 if (abs(idDaughter)==5 && tnow <= m2bPhys) { dip.pT2 = 0.0;
return false;}
3292 double tnew = (useFixedFacScale) ? fixedFacScale2 : factorMultFac*tnow;
3293 tnew = max(tnew, pT2min);
3294 bool inNew = (hasPDFdau) ? beam.insideBounds(xDaughter, tnew) :
true;
3295 if (!inNew && hasPDFdau) { dip.pT2 = 0.0;
return false; }
3299 if (hasTinyPDFdau) ++loopTinyPDFdau;
3300 if ( hasPDFdau && loopTinyPDFdau > MAXLOOPTINYPDF) {
3301 infoPtr->errorMsg(
"Warning in DireSpace::pT2nextQCD_IF: "
3302 "small daughter PDF");
3311 || tnow < evalpdfstep(event[iRadi].
id(), tnow, m2cPhys, m2bPhys)*teval) {
3313 hasTinyPDFdau =
false;
3315 newOverestimates.clear();
3321 Lambda2 = Lambda5flav2;
3322 }
else if (tnow > m2c) {
3324 Lambda2 = Lambda4flav2;
3327 Lambda2 = Lambda3flav2;
3331 Lambda2 /= renormMultFac;
3332 zMinAbs = (hasPDFdau) ? xDaughter : 0.;
3336 pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac*tnow;
3337 pdfScale2 = max(pdfScale2, pT2min);
3338 xPDFdaughter = getXPDF(idDaughter, xDaughter, pdfScale2, iSysNow, &beam);
3339 if ( hasPDFdau && xPDFdaughter != 0.
3340 && abs(xPDFdaughter) < tinypdf(xDaughter)) {
3341 int sign = (xPDFdaughter > 0.) ? 1 : -1;
3342 xPDFdaughter = sign*tinypdf(xDaughter);
3343 hasTinyPDFdau =
true;
3347 getNewOverestimates( idDaughter, &dip, event, teval,
3348 xDaughter, zMinAbs, zMaxAbs, newOverestimates );
3349 addNewOverestimates(newOverestimates, kernelPDF);
3359 if (kernelPDF < TINYKERNELPDF) { dip.pT2 = 0.0;
return false; }
3360 if (newOverestimates.empty())
return false;
3363 bool forceFixedAs = (tnow < pT2min);
3364 tnow = tNextQCD( &dip, kernelPDF, tnow, pT2endDip, pT2freeze,
3365 (forceBranching ? -1 : 1));
3367 wt = dip.pT2 = tnow = 0.;
3368 double R0 = kernelPDF*rndmPtr->flat();
3369 if (!newOverestimates.empty()) {
3370 if (newOverestimates.lower_bound(R0) == newOverestimates.end())
3371 splittingNowName = newOverestimates.rbegin()->second;
3373 splittingNowName = newOverestimates.lower_bound(R0)->second;
3379 if (tnow <= pT2endDip) { dip.pT2 = tnow = 0.;
break; }
3383 if (nFlavour == 5 && tnow <= m2bPhys) {
3386 }
else if (nFlavour == 4 && tnow <= m2cPhys) {
3391 if (abs(idDaughter) == 5 && tnow <= m2bPhys) {
3392 dip.pT2 = tnow = 0.;
break;
3393 }
else if (abs(idDaughter) == 4 && tnow <= m2cPhys) {
3394 dip.pT2 = tnow = 0.;
break;
3397 double xMin = xDaughter;
3398 if (!hasPDFdau) xMin = 0.;
3401 double R = kernelPDF*rndmPtr->flat();
3402 if (!newOverestimates.empty()) {
3403 if (newOverestimates.lower_bound(R) == newOverestimates.end())
3404 splittingNowName = newOverestimates.rbegin()->second;
3406 splittingNowName = newOverestimates.lower_bound(R)->second;
3407 getNewSplitting( event, &dip, teval, xMin, tnow, zMinAbs,
3408 zMaxAbs, idDaughter, splittingNowName, forceFixedAs, idMother,
3409 idSister, znow, wt, fullWeightsNow, overWeightNow);
3415 fullWeightsNow.clear();
3416 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3417 nContinue++;
continue;
3420 fullWeightNow = fullWeightsNow[
"base"];
3424 if ( tnow <= 4.*m2bPhys
3425 && ( (abs(idDaughter) == 21 && abs(idSister) == 5)
3426 || (abs(idDaughter) == 5 && splits[splittingNowName]->nEmissions()==2)
3427 || (abs(idSister) == 5 && splits[splittingNowName]->nEmissions()==2))) {
3428 fullWeightsNow.clear();
3429 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3430 nContinue++;
continue;
3431 }
else if ( tnow <= 4.*m2cPhys
3432 && ( (abs(idDaughter) == 21 && abs(idSister) == 4)
3433 || (abs(idDaughter) == 4 && splits[splittingNowName]->nEmissions()==2)
3434 || (abs(idSister) == 4 && splits[splittingNowName]->nEmissions()==2))) {
3435 fullWeightsNow.clear();
3436 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3437 nContinue++;
continue;
3444 double m2e = (abs(idSister)<6) ? getMass(idSister,2) : getMass(idSister,1);
3445 if (!forceMassiveMap) m2e = 0.0;
3446 double m2s = particleDataPtr->isResonance(event[dip.iRecoiler].id())
3447 ? getMass(event[dip.iRecoiler].id(),3,
3448 event[dip.iRecoiler].mCalc())
3449 : (event[dip.iRecoiler].idAbs() < 6)
3450 ? getMass(event[dip.iRecoiler].id(),2)
3451 : getMass(event[dip.iRecoiler].id(),1);
3457 if ( splits[splittingNowName]->nEmissions() == 2 )
3458 if ( (abs(idSister) == 4 && tnow < m2cPhys)
3459 || (abs(idSister) == 5 && tnow < m2bPhys)) {
3461 fullWeightsNow.clear();
3462 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3463 nContinue++;
continue;
3467 double jacobian(1.), m2aij(m2Bef), m2ai(0.), m2a(m2r),
3468 m2i(getMass(idMother,2)),
3469 m2j(getMass(-event[iRadi].
id(),2)), m2k(m2s);
3470 m2ai = -dip.sa1 + m2a + m2i;
3471 double q2 = (
event[iRadi].p()-
event[iReco].p()).m2Calc();
3473 bool canUseSplitInfo = splits[splittingNowName]->canUseForBranching();
3474 if (canUseSplitInfo) {
3476 = splits[splittingNowName]->getJacobian(event,partonSystemsPtr);
3477 unordered_map<string,double> psvars = splits[splittingNowName]->
3478 getPhasespaceVars( event, partonSystemsPtr);
3479 xMother = psvars[
"xInAft"];
3483 if ( splits[splittingNowName]->nEmissions() == 2 ) {
3484 double m2jk = dip.pT2/dip.xa + q2*( 1. - dip.xa/dip.z) - m2ai;
3488 double uCS = dip.z*(m2ai-m2a-m2i)/q2;
3489 double xCS = uCS + dip.xa - (dip.pT2*dip.z)/(q2*dip.xa);
3490 Vec4 q( event[iRadi].p() - event[iReco].p() );
3491 double sHatBef = (
event[iRadi].p() + pOther).m2Calc();
3492 double sijk = q2*(1.-1./dip.z) - m2a;
3498 if (!useGlobalMapIF) {
3501 Vec4 pTk_tilde(event[iReco].p().px(), event[iReco].p().py(), 0., 0.);
3502 Vec4 qpar( q + pTk_tilde );
3504 double q2par = qpar.m2Calc();
3505 double pT2k = -pTk_tilde.m2Calc();
3506 double s_i_jk = (1. - 1./xCS)*(q2 - m2a) + (m2i + m2jk) / xCS;
3508 Vec4 pa( ( event[iRadi].p() - 0.5*(q2-m2aij-m2k)/q2par * qpar )
3509 * sqrt( (lABC(q2,s_i_jk,m2a) - 4.*m2a*pT2k)
3510 / (lABC(q2,m2k,m2aij) - 4.*m2aij*pT2k))
3511 + qpar * 0.5 * (q2 + m2a - s_i_jk) / q2par);
3513 sHatAft = (pa + pOther).m2Calc();
3523 Vec4 pb_tilde( event[iReco].p() );
3524 Vec4 pa12_tilde( event[iRadi].p() );
3525 q.p(pb_tilde-pa12_tilde);
3528 double zbar = (q2-m2ai-m2jk) / bABC(q2,m2ai,m2jk)
3529 *( (xCS - 1)/(xCS-uCS) - m2jk / gABC(q2,m2ai,m2jk)
3530 * (m2ai + m2i - m2a) / (q2 - m2ai - m2jk));
3531 double kT2 = zbar*(1.-zbar)*m2ai - (1-zbar)*m2i - zbar*m2a;
3534 Vec4 pjk( (pb_tilde - q*pb_tilde/q2*q)
3535 *sqrt(lABC(q2,m2ai,m2jk)/lABC(q2,m2aij,m2k))
3536 + 0.5*(q2+m2jk-m2ai)/q2*q );
3542 pair<Vec4, Vec4> pTvecs = getTwoPerpendicular(pai, pjk);
3543 Vec4 kTmom( sqrt(kT2)*sin(dip.phi)*pTvecs.first
3544 + sqrt(kT2)*cos(dip.phi)*pTvecs.second);
3547 Vec4 pi( - zbar *(gABC(q2,m2ai,m2jk)*pai + m2ai*pjk)
3548 / bABC(q2,m2ai,m2jk)
3549 + ( (1.-zbar)*m2ai + m2i - m2a) / bABC(q2,m2ai,m2jk)
3550 * (pjk + m2jk/gABC(q2,m2ai,m2jk)*pai)
3557 sHatAft = (pa + pOther).m2Calc();
3561 double m2Other = pOther.m2Calc();
3562 double rho_aij = sqrt( lABC(sHatBef, m2a, m2Other)
3563 /lABC(sHatAft, m2a, m2Other));
3564 jacobian = rho_aij / dip.z * (sijk + m2a - q2)
3565 / sqrt(lABC(sijk, m2a, q2));
3567 jacobian *= -q2 * dip.xa / dip.z / sqrt(lABC(m2jk, m2ai, q2));
3569 jacobian *= 1. / (1. - (m2ai + m2j - m2aij) / (dip.pT2/dip.xa)) ;
3574 xMother = xDaughter/xCS;
3578 fullWeightNow *= jacobian;
3579 for ( unordered_map<string,double>::iterator it = fullWeightsNow.begin();
3580 it != fullWeightsNow.end(); ++it )
3581 it->second *= jacobian;
3584 double pdfRatio = 1.;
3585 pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * tnow;
3586 pdfScale2 = max(pdfScale2, pT2min);
3587 double pdfScale2Old = pdfScale2;
3588 double pdfScale2New = pdfScale2;
3589 if (forceBranching) pdfScale2Old = pdfScale2New = infoPtr->Q2Fac();
3590 bool inD = (hasPDFdau) ? beam.insideBounds(xDaughter, pdfScale2Old) :
true;
3591 bool inM = (hasPDFdau) ? beam.insideBounds(xMother, pdfScale2New) :
true;
3592 double xPDFdaughterNew = getXPDF( idDaughter, xDaughter, pdfScale2Old,
3593 iSysNow, &beam,
true, znow, m2Dip);
3594 double xPDFmotherNew = getXPDF( idMother, xMother, pdfScale2New,
3595 iSysNow, &beam,
true, znow, m2Dip);
3597 if ( hasPDFdau && xPDFdaughterNew != 0.
3598 && abs(xPDFdaughterNew) < tinypdf(xDaughter) ) {
3599 hasTinyPDFdau =
true;
3601 fullWeightsNow.clear();
3602 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3603 nContinue++;
continue;
3610 double xPDFdaughterLow = getXPDF( idDaughter, xDaughter,
3611 pdfScale2Old * pdfScale2Old/max(teval,pT2min), iSysNow, &beam);
3612 if ( hasPDFdau && idDaughter == 21
3613 && ( abs(xPDFdaughterNew/xPDFdaughter) < 1e-4
3614 || abs(xPDFdaughterLow/xPDFdaughterNew) < 1e-4) ) {
3615 hasTinyPDFdau =
true;
3617 fullWeightsNow.clear();
3618 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3619 nContinue++;
continue;
3623 pdfRatio = (hasPDFdau) ?
3624 ((inD && inM) ? xPDFmotherNew/xPDFdaughterNew : 0.) : 1.;
3627 if (idDaughter == 21 && pdfScale2 < 1.01 && pdfRatio > 50.) pdfRatio = 0.;
3629 fullWeightNow *= pdfRatio;
3630 for ( unordered_map<string,double>::iterator it = fullWeightsNow.begin();
3631 it != fullWeightsNow.end(); ++it )
3632 it->second *= pdfRatio;
3636 if ( splits[splittingNowName]->nEmissions() == 2 )
3637 dip.sa1 = splits[splittingNowName]->splitInfo.kinematics()->sai;
3639 if (fullWeightNow == 0.) {
3641 fullWeightsNow.clear();
3642 wt = fullWeightNow = overWeightNow = auxWeightNow = 0.;
3643 nContinue++;
continue;
3647 double m2DipCorr = dip.m2Dip - m2Bef + m2r + m2e;
3648 double scale2 = splits[splittingNowName]->couplingScale2 (
3649 dip.z, tnow, m2DipCorr,
3650 make_pair (event[dip.iRadiator].id(),
event[dip.iRadiator].isFinal()),
3651 make_pair (event[dip.iRecoiler].id(),
event[dip.iRecoiler].isFinal()));
3652 if (scale2 < 0.) scale2 = tnow;
3653 double talpha = max(scale2, pT2min);
3657 alphasReweight(tnow, talpha, dip.system, forceFixedAs, wt, fullWeightNow,
3658 overWeightNow, renormMultFac);
3662 alphasReweight(tnow, talpha, dip.system, forceFixedAs, daux, asw, daux,
3664 fullWeightsNow[
"base"] *= asw;
3665 if (fullWeightsNow.find(
"base_order_as2") != fullWeightsNow.end())
3666 fullWeightsNow[
"base_order_as2"] *= asw;
3668 if ( splittingNowName.find(
"qcd") != string::npos
3669 && settingsPtr->parm(
"Variations:muRisrDown") != 1.) {
3671 alphasReweight(tnow, talpha, dip.system, forceFixedAs, daux, asw, daux,
3672 (tnow > pT2minVariations) ? settingsPtr->parm
3673 (
"Variations:muRisrDown")*renormMultFac
3675 fullWeightsNow[
"Variations:muRisrDown"] *= asw;
3676 }
else if ( splittingNowName.find(
"qcd") == string::npos )
3677 fullWeightsNow[
"Variations:muRisrDown"] *= asw;
3678 if ( splittingNowName.find(
"qcd") != string::npos
3679 && settingsPtr->parm(
"Variations:muRisrUp") != 1.) {
3681 alphasReweight(tnow, talpha, dip.system, forceFixedAs, daux, asw, daux,
3682 (tnow > pT2minVariations) ? settingsPtr->parm
3683 (
"Variations:muRisrUp")*renormMultFac
3685 fullWeightsNow[
"Variations:muRisrUp"] *= asw;
3686 }
else if ( splittingNowName.find(
"qcd") == string::npos )
3687 fullWeightsNow[
"Variations:muRisrUp"] *= asw;
3690 if (hasPDFdau && settingsPtr->flag(
"Variations:PDFup") ) {
3691 int valSea = (beam[iSysNow].isValence()) ? 1 : 0;
3692 if( beam[iSysNow].isUnmatched() ) valSea = 2;
3693 beam.calcPDFEnvelope( make_pair(idMother, idDaughter),
3694 make_pair(xMother, xDaughter), pdfScale2, valSea);
3695 PDF::PDFEnvelope ratioPDFEnv = beam.getPDFEnvelope();
3697 = min(ratioPDFEnv.errplusPDF / ratioPDFEnv.centralPDF, 10.);
3698 double deltaPDFminus
3699 = min(ratioPDFEnv.errminusPDF / ratioPDFEnv.centralPDF, 10.);
3700 fullWeightsNow[
"Variations:PDFup"] = fullWeightsNow[
"base"]
3701 * ((tnow > pT2minVariations) ? (1.0 + deltaPDFplus) : 1.0);
3702 fullWeightsNow[
"Variations:PDFdown"] = fullWeightsNow[
"base"]
3703 * ((tnow > pT2minVariations) ? (1.0 - deltaPDFminus) : 1.0);
3707 auxWeightNow = overWeightNow;
3710 if (fullWeightNow < 0.) {
3711 auxWeightNow *= -1.;
3715 if ( fullWeightNow/auxWeightNow > 1.) {
3716 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
3717 <<
" " << __LINE__ <<
" : Large acceptance weight="
3718 << fullWeightNow/auxWeightNow
3719 <<
" for splitting " << splittingNowName <<
" at pT2=" << tnow
3720 <<
" and z=" << znow <<
"\t(PDF ratio=" << pdfRatio <<
")" <<
" "
3722 double rescale = fullWeightNow/auxWeightNow * 1.15;
3723 auxWeightNow *= rescale;
3725 infoPtr->errorMsg(
"Info in DireSpace::pT2nextQCD_IF: Found large "
3726 "acceptance weight for " + splittingNowName);
3729 wt = fullWeightNow/auxWeightNow;
3732 }
while (wt < rndmPtr->flat()) ;
3735 if ( wt == 0.) { dip.pT2 = 0.0;
return false; }
3739 if ( fullWeightNow != 0. && overWeightNow != 0. ) {
3740 double enhanceFurther
3741 = enhanceOverestimateFurther(splittingNowName, idDaughter, teval);
3743 weights->addTrialEnhancement(tnow, enhanceFurther);
3744 enhanceFurther = 1.;
3746 kernelNow = fullWeightsNow;
3747 auxNow = auxWeightNow;
3748 overNow = overWeightNow;
3749 boostNow = enhanceFurther;
3750 for ( unordered_map<string,double>::iterator it = fullWeightsNow.begin();
3751 it != fullWeightsNow.end(); ++it ) {
3754 if (it->first ==
"base_order_as2")
continue;
3756 acceptProbability[it->first].insert(make_pair(tnow,
3757 auxWeightNow/overWeightNow * 1./enhanceFurther
3758 * it->second/fullWeightNow ) );
3759 if (auxWeightNow == fullWeightNow && overWeightNow == fullWeightNow)
3760 rejectProbability[it->first].insert( make_pair(tnow, 1.0));
3762 double wv = auxWeightNow/overWeightNow
3763 * (overWeightNow- it->second/enhanceFurther)
3764 / (auxWeightNow - fullWeightNow);
3765 if (abs(wv) > 1e0) {
3766 direInfoPtr->message(1) << __FILE__ <<
" " << __func__
3767 <<
" " << __LINE__ <<
" : Large reject weight=" << wv
3768 <<
"\t for kernel=" << it->second <<
" overestimate=" << overNow
3769 <<
"\t aux. overestimate=" << auxNow <<
" at pT2="
3771 <<
" for " << splittingNowName << endl;
3773 rejectProbability[it->first].insert( make_pair(tnow, wv));
3778 double zStore = dip.z;
3779 double xaStore = dip.xa;
3780 double pT2Store = dip.pT2;
3781 double sa1Store = dip.sa1;
3782 double Q2store = 0.;
3785 dipEndNow->store( idDaughter, idMother, idSister,
3786 x1Now, x2Now, m2Dip, pT2Store, zStore, sa1Store, xaStore, xMother,
3787 Q2store, mSister, m2Sister, pT2corr, dip.phi, dip.phia1);
3796 double DireSpace::tNextQCD( DireSpaceEnd*,
double overestimateInt,
3797 double tOld,
double tMin,
double tFreeze,
int algoType) {
3799 bool forceFixedAs = (tOld < pT2min);
3800 double asOver = (usePDFalphas || forceFixedAs)
3801 ? alphaS2piOverestimate : alphaS2pi;
3803 double rnd = rndmPtr->flat();
3806 if (usePDFalphas || alphaSorder == 0) {
3807 double rndMin = pow( tMin/tOld, asOver * overestimateInt);
3808 if (rnd < rndMin)
return -1.*tMin;
3813 double Lambda2 = Lambda3flav2;
3816 Lambda2 = Lambda5flav2;
3817 }
else if (tOld > m2c) {
3819 Lambda2 = Lambda4flav2;
3822 Lambda2 = Lambda3flav2;
3825 Lambda2 /= renormMultFac;
3832 return pow(tMin+tFreeze,rnd) / pow(tnow+tFreeze,rnd-1) - tFreeze;
3834 if (usePDFalphas || forceFixedAs)
3835 tnow = (tnow+tFreeze) * pow( rnd,
3836 1. / (alphaS2piOverestimate * overestimateInt)) - tFreeze;
3838 else if (alphaSorder == 0)
3839 tnow = (tnow+tFreeze) * pow( rnd,
3840 1. / (alphaS2pi * overestimateInt) ) - tFreeze;
3842 else if (alphaSorder == 1)
3843 tnow = Lambda2 * pow( (tnow+tFreeze) / Lambda2,
3844 pow( rnd, b0 / overestimateInt) ) - tFreeze;
3848 tnow = Lambda2 * pow( (tnow+tFreeze) / Lambda2,
3849 pow(rndmPtr->flat(), b0 / overestimateInt) ) - tFreeze;
3850 Q2alphaS = renormMultFac * max( tnow+tFreeze,
3851 pow2(LAMBDA3MARGIN) * Lambda3flav2);
3852 }
while (alphaS.alphaS2OrdCorr(Q2alphaS) < rndmPtr->flat()
3865 bool DireSpace::zCollNextQCD( DireSpaceEnd* dip,
double zMin,
double zMax,
3869 dip->xa = zMax * pow( zMax/zMin, -rndmPtr->flat());
3880 bool DireSpace::virtNextQCD( DireSpaceEnd* dip,
double,
double,
3883 double v = (dip->z/dip->xa) * rndmPtr->flat();
3884 double m2j = dip->mass[2];
3885 dip->sa1 = v / (dip->z/dip->xa-v) * ( dip->pT2/dip->xa - m2j);
3886 if (abs(dip->z/dip->xa-v) < 1e-10)
return false;
3898 bool DireSpace::branch(
Event& event) {
3900 if (abs(dipEndSel->pT2 - pT2cutMin(dipEndSel)) < 1e-10)
return false;
3904 bool hasBranched =
false;
3905 if ( event[dipEndSel->iRecoiler].isFinal() )
3906 hasBranched = branch_IF(event,
false, &splitInfoSel);
3907 else hasBranched = branch_II(event,
false, &splitInfoSel);
3919 bool DireSpace::branch_II(
Event& event,
bool trial,
3920 DireSplitInfo* split ) {
3922 Event auxevent1 = event;
3923 Event auxevent2 = event;
3926 int side = (!trial) ? abs(dipEndSel->side) : split->side;
3929 int iDaughter = (!trial) ? dipEndSel->iRadiator : split->iRadBef;
3930 int iRecoiler = (!trial) ? dipEndSel->iRecoiler : split->iRecBef;
3931 int idDaughterNow = (!trial) ? dipEndSel->idDaughter : split->radBef()->id;
3932 int idMother = (!trial) ? dipEndSel->idMother : split->radAft()->id;
3933 int idSister = (!trial) ? dipEndSel->idSister : split->emtAft()->id;
3934 int colDaughter =
event[iDaughter].col();
3935 int acolDaughter =
event[iDaughter].acol();
3936 int colRecBef =
event[iRecoiler].col();
3937 int acolRecBef =
event[iRecoiler].acol();
3938 bool colMatch = (acolRecBef == colDaughter);
3939 bool acolMatch = (colRecBef == acolDaughter);
3942 string name = (!trial) ? splittingSelName : split->splittingSelName;
3943 int nEmissions = splits[name]->nEmissions();
3945 if ( nEmissions == 2 ) idSister = -
event[iDaughter].id();
3948 double pT2 = (!trial) ? dipEndSel->pT2 : split->kinematics()->pT2;
3949 double z = (!trial) ? dipEndSel->z : split->kinematics()->z;
3950 splits[name]->splitInfo.store(*split);
3951 unordered_map<string,double> psp(splits[name]->getPhasespaceVars
3952 (event, partonSystemsPtr));
3954 if (split->useForBranching) { pT2 = psp[
"pT2"]; z = psp[
"z"]; }
3957 double m2Bef = 0.0, m2r = 0.0, m2s = 0.0;
3959 double m2ex = (abs(idSister) < 6) ?
3960 getMass(idSister, 2) : getMass(idSister, 1);
3961 double m2e = (!trial) ? m2ex
3962 : ( (split->kinematics()->m2EmtAft > 0.) ? split->kinematics()->m2EmtAft
3965 if ( useMassiveBeams && (abs(idDaughter) == 11 || abs(idDaughter) == 13
3966 || abs(idDaughter) > 900000))
3967 m2Bef = getMass(idDaughter,1);
3968 if ( useMassiveBeams && (abs(idMother) == 11 || abs(idMother) == 13
3969 || abs(idMother) > 900000))
3970 m2r = getMass(idMother,1);
3972 if ( useMassiveBeams && (event[iRecoiler].idAbs() == 11
3973 || event[iRecoiler].idAbs() == 13
3974 || event[iRecoiler].idAbs() > 900000))
3975 m2s = getMass(event[iRecoiler].
id(),1);
3977 if ( useMassiveBeams && (abs(idSister) == 11 || abs(idSister) == 13
3978 || abs(idSister) > 900000))
3979 m2e = getMass(idSister,1);
3982 if (!forceMassiveMap) m2e = 0.0;
3983 double m2dip = (!trial) ? dipEndSel->m2Dip : split->kinematics()->m2Dip;
3984 double m2DipCorr = m2dip - m2Bef + m2r + m2e;
3985 double kappa2 = pT2 / m2DipCorr;
3986 double xCS = (z*(1-z)- kappa2)/(1-z);
3987 double vCS = kappa2/(1-z);
3988 double sai = (!trial) ? dipEndSel->sa1 : split->kinematics()->sai;
3989 double xa = (!trial) ? dipEndSel->xa : split->kinematics()->xa;
3991 if (split->useForBranching) { sai = psp[
"sai"]; xa = psp[
"xa"]; }
3996 Vec4 paj_tilde(event[iDaughter].p());
3997 Vec4 pb_tilde(event[iRecoiler].p());
3998 Vec4 q(pb_tilde+paj_tilde);
3999 double q2 = q.m2Calc();
4002 int eventSizeOld =
event.size();
4003 int iSysSelNow = (!trial) ? iSysSel : split->system;
4004 int systemSizeOld = partonSystemsPtr->sizeAll(iSysSelNow);
4007 int beamOff1 = 1 + beamOffset;
4008 int beamOff2 = 2 + beamOffset;
4009 int ev1Dau1V =
event[beamOff1].daughter1();
4010 int ev2Dau1V =
event[beamOff2].daughter1();
4011 vector<int> statusV, mother1V, mother2V, daughter1V, daughter2V;
4014 bool physical =
true;
4015 bool canMergeFirst = (mergingHooksPtr != 0)
4016 ? mergingHooksPtr->canVetoEmission() :
false;
4017 for (
int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
4018 int iOldCopy = partonSystemsPtr->getAll(iSysSelNow, iCopy);
4019 statusV.push_back( event[iOldCopy].status());
4020 mother1V.push_back( event[iOldCopy].mother1());
4021 mother2V.push_back( event[iOldCopy].mother2());
4022 daughter1V.push_back( event[iOldCopy].daughter1());
4023 daughter2V.push_back( event[iOldCopy].daughter2());
4028 int iMother(0), iNewRecoiler(0);
4029 for (
int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
4030 int iOldCopy = partonSystemsPtr->getAll(iSysSelNow, iCopy);
4031 int statusOld =
event[iOldCopy].status();
4032 int statusNew = (iOldCopy == iDaughter
4033 || iOldCopy == iRecoiler) ? statusOld : 44;
4034 int iNewCopy =
event.copy(iOldCopy, statusNew);
4035 auxevent1.copy(iOldCopy, statusNew);
4036 auxevent2.copy(iOldCopy, statusNew);
4037 if (iOldCopy == iDaughter) iMother = iNewCopy;
4038 if (iOldCopy == iRecoiler) iNewRecoiler = iNewCopy;
4039 if (statusOld < 0) {
4040 event[iNewCopy].statusNeg();
4041 auxevent1[iNewCopy].statusNeg();
4042 auxevent2[iNewCopy].statusNeg();
4047 int idMotherNow = idMother;
4048 if ( nEmissions == 2 ) idMotherNow = 21;
4052 int colMother = colDaughter;
4053 int acolMother = acolDaughter;
4057 double RN2 = rndmPtr->flat();
4058 if (idSister == 22 || idSister == 23 || idSister == 25) ;
4061 else if (idSister == 21 && ( (idMotherNow > 0 && idMotherNow < 9)
4062 || ( idMother == 21 && colMatch && acolMatch && RN2 < 0.5) ) ) {
4063 colMother =
event.nextColTag();
4064 colSister = colMother;
4065 acolSister = colDaughter;
4067 }
else if (idSister == 21 && ( (idMotherNow < 0 && idMotherNow > -9)
4068 || ( idMotherNow == 21 && colMatch && acolMatch) ) ) {
4069 acolMother =
event.nextColTag();
4070 acolSister = acolMother;
4071 colSister = acolDaughter;
4072 }
else if (idMotherNow == 21 && idSister == 21 && colMatch && !acolMatch) {
4073 colMother =
event.nextColTag();
4074 acolMother = acolDaughter;
4075 colSister = colMother;
4076 acolSister = colDaughter;
4077 }
else if (idMotherNow == 21 && idSister == 21 && !colMatch && acolMatch) {
4078 colMother = colDaughter;
4079 acolMother =
event.nextColTag();
4080 acolSister = acolMother;
4081 colSister = acolDaughter;
4083 }
else if (idDaughterNow == 21 && idMotherNow > 0) {
4084 colMother = colDaughter;
4086 colSister = acolDaughter;
4088 }
else if (idDaughterNow == 21) {
4089 acolMother = acolDaughter;
4091 acolSister = colDaughter;
4093 }
else if (idDaughterNow > 0 && idDaughterNow < 9) {
4094 acolMother =
event.nextColTag();
4095 acolSister = acolMother;
4097 }
else if (idDaughterNow < 0 && idDaughterNow > -9) {
4098 colMother =
event.nextColTag();
4099 colSister = colMother;
4103 int colSave = colSister, acolSave = acolSister;
4104 if ( splits[name]->is(splittingsPtr->isrQCD_21_to_21_and_21b)) {
4105 colSister = acolMother;
4106 acolSister = colMother;
4107 colMother = acolSave;
4108 acolMother = colSave;
4111 if (split->useForBranching) {
4112 idMotherNow=
event[iDaughter].id();
4113 colMother =
event[iDaughter].col();
4114 acolMother =
event[iDaughter].acol();
4118 if (split->radAft()->id != 0) idMotherNow= split->radAft()->id;
4119 if (split->emtAft()->id != 0) idSister = split->emtAft()->id;
4120 if (split->radAft()->col > -1) colMother = split->radAft()->col;
4121 if (split->radAft()->acol > -1) acolMother = split->radAft()->acol;
4122 if (split->emtAft()->col > -1) colSister = split->emtAft()->col;
4123 if (split->emtAft()->acol > -1) acolSister = split->emtAft()->acol;
4126 int colMother1, acolMother1;
4127 int colSister1, acolSister1;
4128 colMother1 = acolMother1 = colSister1 = acolSister1 = 0;
4129 if ( nEmissions == 2 ) {
4132 if (idMother*idDaughterNow > 0 && idMother > 0) {
4133 colMother1 = colDaughter;
4136 colSister1 = acolSister;
4140 if (idMother*idDaughterNow > 0 && idMother < 0) {
4142 acolMother1 = acolDaughter;
4143 acolSister1 = colSister;
4148 if (idMother*idDaughterNow < 0 && idMother > 0) {
4149 colMother1 = colSister;
4152 colSister1 = acolDaughter;
4154 acolMother = acolDaughter;
4155 colMother = colSister;
4159 if (idMother*idDaughterNow < 0 && idMother < 0) {
4161 acolMother1 = acolSister;
4162 acolSister1 = colDaughter;
4165 acolMother = acolSister;
4166 colMother = colDaughter;
4173 if ( nEmissions == 2 )
4174 iMother1 =
event.append( 0, 0, 0, 0, 0, 0, 0, 0, Vec4(0.,0.,0.,0.), 0.0,
4177 int iSister =
event.append( idSister, 43, iMother, 0, 0, 0,
4178 colSister, acolSister, Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
4179 auxevent1.append( idSister, 43, iMother, 0, 0, 0,
4180 colSister, acolSister, Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
4181 auxevent2.append( idSister, 43, iMother, 0, 0, 0,
4182 colSister, acolSister, Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
4186 if ( nEmissions == 2 )
4187 iSister1 =
event.append( 0, 0, 0, 0, 0, 0, 0, 0, Vec4(0.,0.,0.,0.), 0.0,
4191 Particle& daughter =
event[iDaughter];
4192 Particle& mother =
event[iMother];
4193 Particle& newRecoiler =
event[iNewRecoiler];
4194 Particle& sister =
event[iSister];
4195 Particle& mother1 =
event[iMother1];
4196 Particle& sister1 =
event[iSister1];
4199 mother.id( idMotherNow );
4200 mother.status( -41);
4201 mother.cols( colMother, acolMother);
4202 mother.p(0.,0., event[iDaughter].pz()/xCS, event[iDaughter].e()/xCS);
4203 if (mother.idAbs() == 21 || mother.idAbs() == 22) mother.pol(9);
4204 newRecoiler.status(-42);
4205 newRecoiler.p(0.,0., event[iRecoiler].pz(), event[iRecoiler].e());
4208 daughter.mothers( iMother, 0);
4209 mother.daughters( iSister, iDaughter);
4210 int iBeam1Dau1 =
event[beamOff1].daughter1();
4211 int iBeam2Dau1 =
event[beamOff2].daughter1();
4212 if (iSysSelNow == 0) {
4213 event[beamOff1].daughter1( (side == 1) ? iMother : iNewRecoiler );
4214 event[beamOff2].daughter1( (side == 2) ? iMother : iNewRecoiler );
4217 auxevent1[iMother].id( idMotherNow );
4218 auxevent1[iMother].status( -41);
4219 auxevent1[iMother].cols( colMother, acolMother);
4220 auxevent1[iMother].p(0.,0., event[iDaughter].pz()/xCS,
4221 event[iDaughter].e()/xCS);
4222 if (auxevent1[iMother].idAbs() == 21 || auxevent1[iMother].idAbs() == 22)
4223 auxevent1[iMother].pol(9);
4224 auxevent1[iNewRecoiler].status(-42);
4225 auxevent1[iNewRecoiler].p(0.,0., auxevent1[iRecoiler].pz(),
4226 auxevent1[iRecoiler].e());
4227 auxevent1[iDaughter].mothers( iMother, 0);
4228 auxevent1[iMother].daughters( iSister, iDaughter);
4229 if (iSysSelNow == 0) {
4230 auxevent1[beamOff1].daughter1( (side == 1) ? iMother : iNewRecoiler );
4231 auxevent1[beamOff2].daughter1( (side == 2) ? iMother : iNewRecoiler );
4234 auxevent2[iMother].id( idMotherNow );
4235 auxevent2[iMother].status( -41);
4236 auxevent2[iMother].cols( colMother, acolMother);
4237 auxevent2[iMother].p(0.,0., event[iDaughter].pz()/xCS,
4238 event[iDaughter].e()/xCS);
4239 if (auxevent2[iMother].idAbs() == 21 || auxevent2[iMother].idAbs() == 22)
4240 auxevent2[iMother].pol(9);
4241 auxevent2[iNewRecoiler].status(-42);
4242 auxevent2[iNewRecoiler].p(0.,0., auxevent2[iRecoiler].pz(),
4243 auxevent2[iRecoiler].e());
4244 auxevent2[iDaughter].mothers( iMother, 0);
4245 auxevent2[iMother].daughters( iSister, iDaughter);
4246 if (iSysSelNow == 0) {
4247 auxevent2[beamOff1].daughter1( (side == 1) ? iMother : iNewRecoiler );
4248 auxevent2[beamOff2].daughter1( (side == 2) ? iMother : iNewRecoiler );
4251 bool doVeto =
false;
4252 bool printWarnings = (!trial || forceMassiveMap);
4253 bool doMECreject =
false;
4256 if ( nEmissions != 2 ) {
4259 double sab = (q2 - m2Emt)/xCS + (m2Rad+m2s) * (1-1/xCS);
4260 double saj = -vCS*(sab - m2Rad-m2s) + m2Rad + m2Emt;
4261 double zbar = (sab - m2Rad - m2s) / bABC(sab,m2Rad,m2s)
4262 *( (xCS + vCS) - m2s / gABC(sab,m2Rad,m2s)
4263 * (saj + m2Rad - m2Emt) / (sab - m2Rad - m2s));
4264 double kT2 = zbar*(1.-zbar)*m2Rad - (1-zbar)*saj - zbar*m2Emt;
4269 infoPtr->errorMsg(
"Warning in DireSpace::branch_II: Reject state "
4270 "with kinematically forbidden kT^2.");
4275 Vec4 pRad = (paj_tilde - m2Bef/gABC(q2,m2Bef,m2s)*pb_tilde)
4276 *sqrt(lABC(sab,m2Rad,m2s)/lABC(q2,m2Bef,m2s))
4277 + m2Rad / gABC(sab,m2Rad,m2s)*pb_tilde;
4281 newRecoiler.p(pb_tilde);
4282 mother.m(sqrtpos(m2Rad));
4283 newRecoiler.m(sqrtpos(m2s));
4285 auxevent1[iMother].p(pRad);
4286 auxevent1[iNewRecoiler].p(pb_tilde);
4287 auxevent1[iMother].m(sqrtpos(m2Rad));
4288 auxevent1[iNewRecoiler].m(sqrtpos(m2s));
4290 auxevent2[iMother].p(pRad);
4291 auxevent2[iNewRecoiler].p(pb_tilde);
4292 auxevent2[iMother].m(sqrtpos(m2Rad));
4293 auxevent2[iNewRecoiler].m(sqrtpos(m2s));
4296 Vec4 kTilde(paj_tilde + pb_tilde);
4299 NewEvent.init(
"(hard process-modified)", particleDataPtr);
4301 for (
int i = 0; i <
event.size(); ++i)
4302 NewEvent.append( event[i] );
4305 Vec4 pSum(1e5,1e5,1e5,1e5);
4306 Vec4 pSumIn( beamAPtr->p() + beamBPtr->p() );
4312 while ( abs(pSum.px()-pSumIn.px()) > mTolErr
4313 || abs(pSum.py()-pSumIn.py()) > mTolErr ) {
4318 || (nTries > 1 && split->useForBranching)) {
4320 infoPtr->errorMsg(
"Warning in DireSpace::branch_II: Could not "
4321 "set up state after branching, thus reject.");
4322 physical =
false;
break;
4326 double phi_kt = (!trial) ? 2.*M_PI*rndmPtr->flat()
4327 : (split->kinematics()->phi < 0. ?
4328 2.*M_PI*rndmPtr->flat() : split->kinematics()->phi);
4331 if (split->useForBranching) { phi_kt = psp[
"phi"]; }
4333 double phi_kt1 = phi_kt+DPHI_II*M_PI;
4334 if (phi_kt1>2.*M_PI) phi_kt1 -= 2.*M_PI;
4335 double phi_kt2 = phi_kt-DPHI_II*M_PI;
4336 if (phi_kt2<0.) phi_kt2 += 2.*M_PI;
4337 if (phi_kt1<phi_kt2) swap(phi_kt1, phi_kt2);
4341 pair<Vec4, Vec4> pTvecs = getTwoPerpendicular(pRad, pb_tilde);
4342 Vec4 kTmom( sqrt(kT2)*sin(phi_kt)*pTvecs.first
4343 + sqrt(kT2)*cos(phi_kt)*pTvecs.second);
4346 Vec4 pEmt = (1-zbar) * (gABC(sab,m2Rad,m2s)*pRad - m2Rad*pb_tilde)
4347 / bABC(sab,m2Rad,m2s)
4348 + (m2Emt + kT2) / ((1-zbar)*bABC(sab,m2Rad,m2s))
4349 * (pb_tilde - m2s/gABC(sab,m2Rad,m2s)*pRad)
4354 sister.m(sqrtpos(m2Emt));
4356 kTmom.p( sqrt(kT2)*sin(phi_kt1)*pTvecs.first
4357 + sqrt(kT2)*cos(phi_kt1)*pTvecs.second);
4358 pEmt.p( (1-zbar) * (gABC(sab,m2Rad,m2s)*pRad - m2Rad*pb_tilde)
4359 / bABC(sab,m2Rad,m2s)
4360 + (m2Emt + kT2) / ((1-zbar)*bABC(sab,m2Rad,m2s))
4361 * (pb_tilde - m2s/gABC(sab,m2Rad,m2s)*pRad)
4363 auxevent1[iSister].p(pEmt);
4364 auxevent1[iSister].m(sqrtpos(m2Emt));
4366 kTmom.p( sqrt(kT2)*sin(phi_kt2)*pTvecs.first
4367 + sqrt(kT2)*cos(phi_kt2)*pTvecs.second);
4368 pEmt.p( (1-zbar) * (gABC(sab,m2Rad,m2s)*pRad - m2Rad*pb_tilde)
4369 / bABC(sab,m2Rad,m2s)
4370 + (m2Emt + kT2) / ((1-zbar)*bABC(sab,m2Rad,m2s))
4371 * (pb_tilde - m2s/gABC(sab,m2Rad,m2s)*pRad)
4373 auxevent2[iSister].p(pEmt);
4374 auxevent2[iSister].m(sqrtpos(m2Emt));
4377 vector<int> iPos(1,iSister);
4381 Vec4 k(mother.p() + newRecoiler.p() - sister.p());
4382 Vec4 kSum(kTilde + k);
4383 for (
int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i) {
4384 Vec4 pIn = NewEvent[i].p();
4385 double kSum2 = kSum.m2Calc();
4386 double k2 = k.m2Calc();
4387 double kTildeXp = kTilde*pIn;
4388 double kSumXp = kSum*pIn;
4389 Vec4 res = pIn - kSum * 2.0*( kSumXp / kSum2 )
4390 + k * 2.0 *( kTildeXp/k2);
4392 auxevent1[i].p(res);
4393 auxevent2[i].p(res);
4396 if (!validMomentum(event[i].p(), event[i].
id(), event[i].status()))
4397 pSum +=
event[i].p();
4399 if (event[i].status() > 0) pSum +=
event[i].p();
4403 for (
int i = 0; i <
event.size(); ++i)
4404 if ( event[i].isFinal()
4405 && partonSystemsPtr->getSystemOf(i,
true) == iSysSelNow
4406 && find(iPos.begin(), iPos.end(), i) == iPos.end() )
4407 pSum += event[i].p();
4411 if ( !validMomentum( mother.p(), idMother, -1)
4412 || !validMomentum( sister.p(), idSister, 1)
4413 || !validMomentum( newRecoiler.p(),
event[iNewRecoiler].id(), -1) )
4416 doVeto = (( canVetoEmission && userHooksPtr->doVetoISREmission(
4417 eventSizeOld, event, iSysSelNow))
4418 || ( canMergeFirst && mergingHooksPtr->doVetoEmission(event)) );
4420 double xm = 2.*mother.e() / (beamAPtr->e() + beamBPtr->e());
4423 double xAnew = (mother.mother1() == 1)
4424 ? 2.*mother.e() / (beamAPtr->e() + beamBPtr->e())
4425 : 2.*newRecoiler.e() / (beamAPtr->e() + beamBPtr->e());
4426 double iAold = (mother.mother1() == 1) ? iDaughter : iRecoiler;
4427 double iAnew = (mother.mother1() == 1) ? iMother : iNewRecoiler;
4428 double xBnew = (mother.mother1() == 1)
4429 ? 2.*newRecoiler.e() / (beamAPtr->e() + beamBPtr->e())
4430 : 2.*mother.e() / (beamAPtr->e() + beamBPtr->e());
4431 double iBold = (mother.mother1() == 1) ? iRecoiler : iDaughter;
4432 double iBnew = (mother.mother1() == 1) ? iNewRecoiler : iMother;
4433 if ( hasPDF(event[iAnew].
id()) && beamAPtr->size() > 0) {
4434 double xOld = (*beamAPtr)[iSysSelNow].x();
4435 (*beamAPtr)[iSysSelNow].iPos(iAnew);
4436 (*beamAPtr)[iSysSelNow].x(xAnew);
4437 if (beamAPtr->xMax(-1) < 0.0) {
4438 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_II: "
4439 "used up beam momentum; discard splitting.");
4443 (*beamAPtr)[iSysSelNow].iPos(iAold);
4444 (*beamAPtr)[iSysSelNow].x(xOld);
4446 if ( hasPDF(event[iBnew].
id()) && beamBPtr->size() > 0) {
4447 double xOld = (*beamBPtr)[iSysSelNow].x();
4448 (*beamBPtr)[iSysSelNow].iPos(iBnew);
4449 (*beamBPtr)[iSysSelNow].x(xBnew);
4450 if (beamBPtr->xMax(-1) < 0.0) {
4451 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_II: "
4452 "used up beam momentum; discard splitting.");
4456 (*beamBPtr)[iSysSelNow].iPos(iBold);
4457 (*beamBPtr)[iSysSelNow].x(xOld);
4461 bool isHardSystem = partonSystemsPtr->getSystemOf(iDaughter,
true) == 0
4462 && partonSystemsPtr->getSystemOf(iRecoiler,
true) == 0;
4463 if (isHardSystem && physical && doMEcorrections
4464 && pT2 > pT2minMECs && checkSIJ(event,Q2minMECs)) {
4465 int iA = getInA(iSysSelNow);
4466 int iB = getInB(iSysSelNow);
4467 vector<int> iOut(createvector<int>(0)(0));
4468 for (
int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
4469 iOut.push_back(partonSystemsPtr->getOut(iSysSelNow, iCopy - 2));
4470 bool motherHasPlusPz = (
event[iMother].pz() > 0.);
4471 if (motherHasPlusPz) {
4472 partonSystemsPtr->setInA(iSysSelNow, iMother);
4473 partonSystemsPtr->setInB(iSysSelNow, iNewRecoiler);
4475 partonSystemsPtr->setInA(iSysSelNow, iNewRecoiler);
4476 partonSystemsPtr->setInB(iSysSelNow, iMother);
4478 for (
int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
4479 partonSystemsPtr->setOut(iSysSelNow, iCopy - 2, eventSizeOld + iCopy);
4480 partonSystemsPtr->addOut(iSysSelNow, iSister);
4482 if ( nFinalMaxMECs < 0
4483 || nFinalMaxMECs > partonSystemsPtr->sizeOut(iSysSelNow))
4484 doMECreject = applyMEC (event, split,
4485 createvector<Event>(auxevent1)(auxevent2));
4487 partonSystemsPtr->setInA(iSysSelNow, iA);
4488 partonSystemsPtr->setInB(iSysSelNow, iB);
4489 for (
int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
4490 partonSystemsPtr->setOut(iSysSelNow, iCopy - 2, iOut[iCopy]);
4491 partonSystemsPtr->popBackOut(iSysSelNow);
4495 if (physical && !doVeto && !trial && !doMECreject) updateAfterII(
4497 iDipSel, eventSizeOld, systemSizeOld, event, iDaughter, iMother, iSister,
4498 iNewRecoiler, pT2, xm);
4505 double m2i = getMass(idMother,2);
4506 double m2j = getMass(idSister,2);
4507 double m2ai = -sai + m2a + m2i;
4510 q2 = (
event[iDaughter].p()
4511 +
event[iRecoiler].p()).m2Calc();
4517 double p2ab = q2/za + m2a + m2k;
4518 double zbar = (p2ab - m2a - m2k) / bABC(p2ab,m2a,m2k)
4519 *( xa - m2k / gABC(p2ab,m2a,m2k)
4520 * (m2ai + m2a - m2i) / (p2ab - m2a - m2k));
4521 double kT2 = zbar*(1.-zbar)*m2a - (1-zbar)*m2ai - zbar*m2i;
4524 if (kT2 < 0. || isnan(kT2)) physical =
false;
4527 Vec4 pa = (paj_tilde - m2aij/gABC(q2,m2aij,m2k)*pb_tilde)
4528 *sqrt(lABC(p2ab,m2a,m2k)/lABC(q2,m2aij,m2k))
4529 + m2a / gABC(p2ab,m2a,m2k)*pb_tilde;
4534 mother.m(sqrtpos(m2a));
4536 mother1.m(sqrtpos(m2a));
4537 newRecoiler.p(pb_tilde);
4540 double phi_kt = (!trial)
4541 ? ((dipEndSel->phi > 0.)
4542 ? dipEndSel->phi : 2.*M_PI*rndmPtr->flat())
4543 : ((split->kinematics()->phi > 0.)
4544 ? split->kinematics()->phi : 2.*M_PI*rndmPtr->flat());
4547 if (split->useForBranching) { phi_kt = psp[
"phi"]; }
4555 pair<Vec4, Vec4> pTvecs = getTwoPerpendicular(pa, pijb);
4556 Vec4 kTmom( sqrt(kT2)*sin(phi_kt)*pTvecs.first
4557 + sqrt(kT2)*cos(phi_kt)*pTvecs.second);
4560 Vec4 pi = (1-zbar) * (gABC(p2ab,m2a,m2k)*pa - m2a*pb_tilde)
4561 / bABC(p2ab,m2a,m2k)
4562 + (m2i + kT2) / ((1-zbar)*bABC(p2ab,m2a,m2k))
4563 * (pb_tilde - m2k/gABC(p2ab,m2a,m2k)*pa)
4568 sister1.m(sqrtpos(m2i));
4572 Vec4 pRadBef(pai+pb_tilde);
4575 double phiFF = (!trial)
4576 ? ((dipEndSel->phia1 > 0.)
4577 ? dipEndSel->phia1 : 2.*M_PI*rndmPtr->flat())
4578 : ((split->kinematics()->phi2 > 0.)
4579 ? split->kinematics()->phi2 : 2.*M_PI*rndmPtr->flat());
4582 if (split->useForBranching) { phiFF = psp[
"phi2"]; }
4585 double m2rec = m2ai;
4588 double zCS = pT2/xa / (q2*xa/za + 2.*m2ai);
4589 double yCS = 1. / ( 1. + (q2*xa/za + 2.*m2ai)
4590 / (q2*(xa/za - 1.) + m2ai + m2k - m2j));
4594 double q2tilde = qtilde.m2Calc();
4595 q.p(pRadBef + pRecBef);
4596 q2 = 4.*m2ai + 2.*q2tilde*xa/za + m2k;
4599 double sij = yCS * (q2 - m2rec) + (1.-yCS)*(m2rad+m2emt);
4600 zbar = (q2-sij-m2rec) / bABC(q2,sij,m2rec)
4601 * (zCS - m2rec/gABC(q2,sij,m2rec)
4602 *(sij + m2rad - m2emt)/(q2-sij-m2rec));
4603 kT2 = zbar*(1.-zbar)*sij - (1.-zbar)*m2rad - zbar*m2emt;
4608 infoPtr->errorMsg(
"Warning in DireSpace::branch_II: Reject state "
4609 "with kinematically forbidden kT^2.");
4615 if (physical && (kT2!=kT2 || abs(kT2-kT2) > 1e5) ) {
4617 infoPtr->errorMsg(
"Warning in DireSpace::branch_II: Reject state "
4618 "with not-a-number kT^2 for branching " + name);
4628 pTvecs = getTwoPerpendicular(pRec, pij);
4629 kTmom.p( sqrt(kT2)*sin(phiFF)*pTvecs.first
4630 + sqrt(kT2)*cos(phiFF)*pTvecs.second);
4633 Vec4 pj( zbar * (gABC(q2,sij,m2rec)*pij - sij*pRec) / bABC(q2,sij,m2rec)
4634 + (m2rad+kT2) / (zbar*bABC(q2,sij,m2rec))
4635 * (pRec - m2rec/gABC(q2,sij,m2rec)*pij)
4638 Vec4 qnew(q-pRec-pj);
4642 sister.m(sqrtpos(m2j));
4645 vector<int> iPos(1,iSister1);
4646 iPos.push_back(iSister);
4647 Vec4 pSum = mother1.p()+newRecoiler.p()-sister1.p()-sister.p();
4650 Vec4 kTilde(qtilde);
4652 Vec4 kSum(kTilde + k);
4653 for (
int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i) {
4654 Vec4 pIn =
event[i].p();
4655 double kSum2 = kSum.m2Calc();
4656 double k2 = k.m2Calc();
4657 double kTildeXp = kTilde*pIn;
4658 double kSumXp = kSum*pIn;
4659 Vec4 res = pIn - kSum * 2.0*( kSumXp / kSum2 )
4660 + k * 2.0 *( kTildeXp/k2);
4662 if (i != iSister && i != iSister1) {
4663 pSum -=
event[i].p();
4669 if ( !validMomentum( mother1.p(), idMother, -1)
4670 || !validMomentum( sister.p(), idSister, 1)
4671 || !validMomentum( sister1.p(), idMother, 1)
4672 || !validMomentum( newRecoiler.p(),
event[iNewRecoiler].id(), -1))
4677 Vec4 pk(event[iRecoiler].p());
4680 double saix(2.*pa*pi), sakx(2.*pa*pk), sajx(2.*pa*pj), sikx(2.*pi*pk),
4681 sjkx(2.*pj*pk), sijx(2.*pi*pj);
4682 double pptt = (sajx-sijx)*(sakx-sikx)/(sakx);
4683 double ssaaii = saix;
4684 double zzaa = 2.*
event[iDaughter].p()*
event[iRecoiler].p()/ ( sakx );
4685 double xxaa = (sakx-sikx) / ( sakx );
4687 (abs(pptt-pT2)/abs(pT2) > 1e-3 || abs(ssaaii-sai)/abs(sai) > 1e-3 ||
4688 abs(zzaa-za)/abs(za) > 1e-3 || abs(xxaa-xa)/abs(xa) > 1e-3 )) {
4689 cout << scientific << setprecision(8);
4690 cout <<
"Error in branch_II: Invariant masses after branching do not "
4691 <<
"match chosen values." << endl;
4693 <<
" Q2 " << (
event[iDaughter].p()+
event[iRecoiler].p()).m2Calc()
4697 <<
" xa " << xa << endl;
4698 cout <<
"Generated: "
4699 <<
" Q2 " << sakx-saix-sajx+sijx-sikx-sjkx+m2a+m2i+m2j+m2k
4701 <<
" sai " << ssaaii
4703 <<
" xa " << xxaa << endl;
4709 for (
int i = 0; i <
event.size(); ++i)
4710 if ( event[i].isFinal()
4711 && partonSystemsPtr->getSystemOf(i,
true) == iSysSelNow
4712 && find(iPos.begin(), iPos.end(), i) == iPos.end() )
4713 pSum -= event[i].p();
4714 if (abs(pSum.px()) > mTolErr || abs(pSum.py()) > mTolErr )
4718 double xAnew = (mother.mother1() == 1)
4719 ? 2.*mother1.e() / (beamAPtr->e() + beamBPtr->e())
4720 : 2.*newRecoiler.e() / (beamAPtr->e() + beamBPtr->e());
4721 double iAold = (mother.mother1() == 1) ? iDaughter : iRecoiler;
4722 double iAnew = (mother.mother1() == 1) ? iMother1 : iNewRecoiler;
4723 double xBnew = (mother.mother1() == 1)
4724 ? 2.*newRecoiler.e() / (beamAPtr->e() + beamBPtr->e())
4725 : 2.*mother1.e() / (beamAPtr->e() + beamBPtr->e());
4726 double iBold = (mother.mother1() == 1) ? iRecoiler : iDaughter;
4727 double iBnew = (mother.mother1() == 1) ? iNewRecoiler : iMother1;
4728 if ( hasPDF(event[iAnew].
id()) && beamAPtr->size() > 0) {
4729 double xOld = (*beamAPtr)[iSysSelNow].x();
4730 (*beamAPtr)[iSysSelNow].iPos(iAnew);
4731 (*beamAPtr)[iSysSelNow].x(xAnew);
4732 if (beamAPtr->xMax(-1) < 0.0) {
4733 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_II: "
4734 "used up beam momentum; discard splitting.");
4738 (*beamAPtr)[iSysSelNow].iPos(iAold);
4739 (*beamAPtr)[iSysSelNow].x(xOld);
4741 if ( hasPDF(event[iBnew].
id()) && beamBPtr->size() > 0) {
4742 double xOld = (*beamBPtr)[iSysSelNow].x();
4743 (*beamBPtr)[iSysSelNow].iPos(iBnew);
4744 (*beamBPtr)[iSysSelNow].x(xBnew);
4745 if (beamBPtr->xMax(-1) < 0.0) {
4746 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_II: "
4747 "used up beam momentum; discard splitting.");
4751 (*beamBPtr)[iSysSelNow].iPos(iBold);
4752 (*beamBPtr)[iSysSelNow].x(xOld);
4756 double xm = 2.*mother1.e() / (beamAPtr->e() + beamBPtr->e());
4757 if (physical && !trial) updateAfterII( iSysSelNow, side, iDipSel,
4758 eventSizeOld, systemSizeOld, event, iDaughter, iMother, iSister,
4759 iNewRecoiler, pT2, xm);
4762 mother1.id( idMother);
4763 mother1.status(-41);
4764 mother1.cols(colMother1, acolMother1);
4765 mother1.daughters( iSister1, iMother);
4766 mother1.mothers( mother.mother1(), mother.mother2());
4767 mother.mothers( iMother1, 0);
4768 sister1.id(idMother);
4770 sister1.mothers(iMother1,0);
4771 sister1.cols(colSister1, acolSister1);
4772 sister1.scale(sqrt(pT2));
4774 if (iSysSelNow == 0) {
4775 event[beamOff1].daughter1( (side == 1) ? iMother1 : iNewRecoiler );
4776 event[beamOff2].daughter1( (side == 2) ? iMother1 : iNewRecoiler );
4780 if (physical && !trial) updateAfterII( iSysSelNow, side, iDipSel, 0, 0,
4781 event, iMother, iMother1, iSister1, iNewRecoiler, pT2, xm);
4785 physical = physical && !doVeto;
4789 if (iSysSelNow > 0) {
4790 if (side == 1)
event[beamOff1].daughter1(0);
4791 if (side == 2)
event[beamOff2].daughter1(0);
4796 if ( physical && !trial && !doMECreject
4797 && !validMotherDaughter(event)) {
4799 infoPtr->errorMsg(
"Error in DireSpace::branch_II: Mother-daughter "
4800 "relations after branching not valid.");
4805 if (iSysSelNow > 0) {
4806 if (side == 1)
event[beamOff1].daughter1(iBeam1Dau1);
4807 if (side == 2)
event[beamOff2].daughter1(iBeam2Dau1);
4811 if ( !physical || doMECreject) {
4812 event.popBack( event.size() - eventSizeOld);
4813 if (iSysSelNow == 0) {
4814 event[beamOff1].daughter1( ev1Dau1V);
4815 event[beamOff2].daughter1( ev2Dau1V);
4817 for (
int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
4818 int iOldCopy = partonSystemsPtr->getAll(iSysSelNow, iCopy);
4819 event[iOldCopy].status( statusV[iCopy]);
4820 event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
4821 event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
4829 for ( unordered_map<
string, multimap<double,double> >::iterator
4830 it = rejectProbability.begin(); it != rejectProbability.end(); ++it){
4831 weights->eraseAcceptWeight(pT2, it->first);
4832 weights->eraseRejectWeight(pT2, it->first);
4836 if (!trial && doMECreject) {
4838 weights->calcWeight(pT2,
false,
true);
4841 for ( unordered_map<
string, multimap<double,double> >::iterator
4842 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
4844 for ( unordered_map<
string, map<double,double> >::iterator
4845 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
4853 if (trial) split->storePosAfter(
4854 (nEmissions < 2) ? iMother : iMother1, iNewRecoiler,
4855 iSister, (nEmissions < 2) ? 0 : iSister1);
4861 weights->calcWeight(pT2);
4864 direInfoPtr->updateSoftPos( iDaughter, iSister );
4865 direInfoPtr->updateSoftPosIfMatch( iRecoiler, iNewRecoiler );
4866 if (nEmissions > 1) direInfoPtr->addSoftPos( iSister1 );
4870 for ( unordered_map<
string, multimap<double,double> >::iterator
4871 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
4873 for ( unordered_map<
string, map<double,double> >::iterator
4874 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
4885 void DireSpace::updateAfterII(
int iSysSelNow,
int sideNow,
int iDipSelNow,
4886 int eventSizeOldNow,
int systemSizeOldNow,
Event& event,
int iDaughter,
4887 int iMother,
int iSister,
int iNewRecoiler,
double pT2,
double xNew) {
4890 if (nProposedPT.find(iSysSelNow) != nProposedPT.end())
4891 ++nProposedPT[iSysSelNow];
4893 int idMother =
event[iMother].id();
4894 int idDaughterNow =
event[iDaughter].id();
4895 bool motherHasPlusPz = (
event[iMother].pz() > 0.);
4898 if ( direInfoPtr->isRes(iDaughter) &&
4899 event[iMother].id() !=
event[iDaughter].id() )
4900 direInfoPtr->removeResPos(iDaughter);
4901 if ( particleDataPtr->isResonance(event[iMother].id()) ) {
4902 if ( direInfoPtr->isRes(iDaughter) )
4903 direInfoPtr->updateResPos(iDaughter,iMother);
4905 if ( particleDataPtr->isResonance(event[iNewRecoiler].id()) )
4906 direInfoPtr->addResPos(iNewRecoiler);
4907 if ( particleDataPtr->isResonance(event[iSister].id()) )
4908 direInfoPtr->addResPos(iSister);
4911 if (motherHasPlusPz) {
4912 partonSystemsPtr->setInA(iSysSelNow, iMother);
4913 partonSystemsPtr->setInB(iSysSelNow, iNewRecoiler);
4915 partonSystemsPtr->setInA(iSysSelNow, iNewRecoiler);
4916 partonSystemsPtr->setInB(iSysSelNow, iMother);
4918 for (
int iCopy = 2; iCopy < systemSizeOldNow; ++iCopy) {
4919 int iOut = partonSystemsPtr->getOut(iSysSelNow, iCopy - 2);
4921 direInfoPtr->updateResPosIfMatch ( iOut, eventSizeOldNow + iCopy);
4922 direInfoPtr->updateSoftPosIfMatch( iOut, eventSizeOldNow + iCopy);
4923 partonSystemsPtr->setOut(iSysSelNow, iCopy - 2, eventSizeOldNow + iCopy);
4925 partonSystemsPtr->addOut(iSysSelNow, iSister);
4928 int iA = getInA(iSysSelNow);
4929 int iB = getInB(iSysSelNow);
4930 double shat = (
event[iA].p() +
event[iB].p()).m2Calc();
4931 partonSystemsPtr->setSHat(iSysSelNow, shat);
4934 dipEndSel = &dipEnd[iDipSelNow];
4937 for (
int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
4938 if ( dipEnd[iDip].system == iSysSelNow) {
4939 if (abs(dipEnd[iDip].side) == sideNow) {
4940 dipEnd[iDip].iRadiator = iMother;
4941 dipEnd[iDip].iRecoiler = iNewRecoiler;
4942 if (dipEnd[iDip].colType != 0)
4943 dipEnd[iDip].colType =
event[iMother].colType();
4947 dipEnd[iDip].iRadiator = iNewRecoiler;
4948 dipEnd[iDip].iRecoiler = iMother;
4949 dipEnd[iDip].MEtype = 0;
4954 BeamParticle& beamNow = (sideNow == 1) ? *beamAPtr : *beamBPtr;
4955 beamNow[iSysSelNow].update( iMother, idMother, xNew);
4957 if (idMother != idDaughterNow) {
4958 pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
4959 pdfScale2 = max(pdfScale2, pT2min);
4960 beamNow.xfISR( iSysSelNow, idMother, xNew, pdfScale2);
4961 beamNow.pickValSeaComp();
4963 BeamParticle& beamRec = (sideNow == 1) ? *beamBPtr : *beamAPtr;
4964 beamRec[iSysSelNow].iPos( iNewRecoiler);
4967 update(iSysSelNow,event);
4980 bool DireSpace::branch_IF(
Event& event,
bool trial,
4981 DireSplitInfo* split ) {
4983 Event auxevent1 = event;
4984 Event auxevent2 = event;
4987 int side = (!trial) ? abs(dipEndSel->side) : split->side;
4990 int iDaughter = (!trial) ? dipEndSel->iRadiator : split->iRadBef;
4991 int iRecoiler = (!trial) ? dipEndSel->iRecoiler : split->iRecBef;
4992 int idDaughterNow = (!trial) ? dipEndSel->idDaughter : split->radBef()->id;
4993 int idMother = (!trial) ? dipEndSel->idMother : split->radAft()->id;
4994 int idSister = (!trial) ? dipEndSel->idSister : split->emtAft()->id;
4995 int colDaughter =
event[iDaughter].col();
4996 int acolDaughter =
event[iDaughter].acol();
4997 int colRecBef =
event[iRecoiler].col();
4998 int acolRecBef =
event[iRecoiler].acol();
4999 bool colMatch = (colRecBef == colDaughter);
5000 bool acolMatch = (acolRecBef == acolDaughter);
5001 int iSysSelNow = (!trial) ? iSysSel : 0;
5004 int iOldOther = (side==1) ? getInB(iSysSelNow) : getInA(iSysSelNow);
5005 string name = (!trial) ? splittingSelName : split->splittingSelName;
5006 int nEmissions = splits[name]->nEmissions();
5007 if ( nEmissions == 2) idSister = -
event[iDaughter].id();
5010 double pT2 = (!trial) ? dipEndSel->pT2 : split->kinematics()->pT2;
5011 double z = (!trial) ? dipEndSel->z : split->kinematics()->z;
5012 splits[name]->splitInfo.store(*split);
5013 unordered_map<string,double> psp(splits[name]->getPhasespaceVars
5014 (event, partonSystemsPtr));
5016 if (split->useForBranching) { pT2 = psp[
"pT2"]; z = psp[
"z"]; }
5022 int type = (
event[iRecoiler].isFinal()) ? 1 : -1;
5024 m2s = particleDataPtr->isResonance(event[iRecoiler].
id())
5025 ? getMass(event[iRecoiler].
id(),3,
5026 event[iRecoiler].mCalc())
5027 : (event[iRecoiler].idAbs() < 6)
5028 ? getMass(event[iRecoiler].id(),2)
5029 : getMass(event[iRecoiler].id(),1);
5033 double m2ex = (abs(idSister) < 6) ?
5034 getMass(idSister, 2) : getMass(idSister, 1);
5035 double m2e = (!trial) ? m2ex
5036 : ( (split->kinematics()->m2EmtAft > 0.) ? split->kinematics()->m2EmtAft
5040 if ( useMassiveBeams && (abs(idDaughter) == 11 || abs(idDaughter) == 13
5041 || abs(idDaughter) > 900000))
5042 m2Bef = getMass(idDaughter,1);
5043 if ( useMassiveBeams && (abs(idMother) == 11 || abs(idMother) == 13
5044 || abs(idMother) > 900000))
5045 m2r = getMass(idMother,1);
5047 if ( useMassiveBeams && (abs(idSister) == 11 || abs(idSister) == 13
5048 || abs(idSister) > 900000))
5049 m2e = getMass(idSister,1);
5052 if (!forceMassiveMap) m2e = 0.0;
5054 double m2dip = (!trial) ? dipEndSel->m2Dip : split->kinematics()->m2Dip;
5055 double m2dipCorr = m2dip - m2Bef + m2r + m2e;
5058 double uCS = (pT2/m2dipCorr)/(1-z);
5059 double sai = (!trial) ? dipEndSel->sa1 : split->kinematics()->sai;
5060 double xa = (!trial) ? dipEndSel->xa : split->kinematics()->xa;
5062 if (split->useForBranching) { sai = psp[
"sai"]; xa = psp[
"xa"]; }
5065 int eventSizeOld =
event.size();
5066 int systemSizeOld = partonSystemsPtr->sizeAll(iSysSelNow);
5069 int beamOff1 = 1 + beamOffset;
5070 int beamOff2 = 2 + beamOffset;
5071 int ev1Dau1V =
event[beamOff1].daughter1();
5072 int ev2Dau1V =
event[beamOff2].daughter1();
5073 vector<int> statusV, mother1V, mother2V, daughter1V, daughter2V;
5076 bool physical =
true;
5077 bool canMergeFirst = (mergingHooksPtr != 0)
5078 ? mergingHooksPtr->canVetoEmission() :
false;
5081 if (useGlobalMapIF) {
5082 for (
int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
5083 int iOldCopy = partonSystemsPtr->getAll(iSysSelNow, iCopy);
5084 statusV.push_back( event[iOldCopy].status());
5085 mother1V.push_back( event[iOldCopy].mother1());
5086 mother2V.push_back( event[iOldCopy].mother2());
5087 daughter1V.push_back( event[iOldCopy].daughter1());
5088 daughter2V.push_back( event[iOldCopy].daughter2());
5093 int iDauStatusV =
event[iDaughter].status();
5094 int iDauMot1V =
event[iDaughter].mother1();
5095 int iDauMot2V =
event[iDaughter].mother2();
5096 int iDauDau1V =
event[iDaughter].daughter1();
5097 int iDauDau2V =
event[iDaughter].daughter2();
5098 int iRecStatusV =
event[iRecoiler].status();
5099 int iRecMot1V =
event[iRecoiler].mother1();
5100 int iRecMot2V =
event[iRecoiler].mother2();
5101 int iRecDau1V =
event[iRecoiler].daughter1();
5102 int iRecDau2V =
event[iRecoiler].daughter2();
5106 int iMother(0), iNewRecoiler(0), iNewOther(0);
5107 if (useGlobalMapIF) {
5108 for (
int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
5109 int iOldCopy = partonSystemsPtr->getAll(iSysSelNow, iCopy);
5110 int statusOld =
event[iOldCopy].status();
5111 int statusNew = (iOldCopy == iDaughter
5112 || iOldCopy == iOldOther) ? statusOld :
5113 (iOldCopy == iRecoiler) ? 48 : 44;
5114 int iNewCopy =
event.copy(iOldCopy, statusNew);
5115 if (iOldCopy == iDaughter) iMother = iNewCopy;
5116 if (iOldCopy == iRecoiler) iNewRecoiler = iNewCopy;
5117 if (iOldCopy == iOldOther) iNewOther = iNewCopy;
5118 if (statusOld < 0)
event[iNewCopy].statusNeg();
5123 int idMotherNow = idMother;
5124 if ( nEmissions == 2) idMotherNow = 21;
5128 int colMother = colDaughter;
5129 int acolMother = acolDaughter;
5133 if (idSister == 22 || idSister == 23 || idSister == 25) ;
5136 else if (idSister == 21 && ( (idMotherNow > 0 && idMotherNow < 9)
5137 || ( idMotherNow == 21 && colMatch && acolMatch
5138 && rndmPtr->flat() < 0.5) ) ) {
5139 colMother =
event.nextColTag();
5140 colSister = colMother;
5141 acolSister = colDaughter;
5143 }
else if (idSister == 21 && ( (idMotherNow < 0 && idMotherNow > -9)
5144 || ( idMotherNow == 21 && colMatch && acolMatch) ) ) {
5145 acolMother =
event.nextColTag();
5146 acolSister = acolMother;
5147 colSister = acolDaughter;
5148 }
else if (idMotherNow == 21 && idSister == 21 && colMatch && !acolMatch) {
5149 colMother =
event.nextColTag();
5150 acolMother = acolDaughter;
5151 colSister = colMother;
5152 acolSister = colDaughter;
5153 }
else if (idMotherNow == 21 && idSister == 21 && !colMatch && acolMatch) {
5154 colMother = colDaughter;
5155 acolMother =
event.nextColTag();
5156 acolSister = acolMother;
5157 colSister = acolDaughter;
5159 }
else if (idDaughterNow == 21 && idMotherNow > 0) {
5160 colMother = colDaughter;
5162 colSister = acolDaughter;
5164 }
else if (idDaughterNow == 21) {
5165 acolMother = acolDaughter;
5167 acolSister = colDaughter;
5169 }
else if (idDaughterNow > 0 && idDaughterNow < 9) {
5170 acolMother =
event.nextColTag();
5171 acolSister = acolMother;
5173 }
else if (idDaughterNow < 0 && idDaughterNow > -9) {
5174 colMother =
event.nextColTag();
5175 colSister = colMother;
5179 int colSave = colSister, acolSave = acolSister;
5180 if ( splits[name]->is(splittingsPtr->isrQCD_21_to_21_and_21b)) {
5181 colSister = acolMother;
5182 acolSister = colMother;
5183 colMother = acolSave;
5184 acolMother = colSave;
5187 if (split->useForBranching) {
5188 idMotherNow=
event[iDaughter].id();
5189 colMother =
event[iDaughter].col();
5190 acolMother =
event[iDaughter].acol();
5194 if (split->radAft()->id != 0) idMotherNow= split->radAft()->id;
5195 if (split->emtAft()->id != 0) idSister = split->emtAft()->id;
5196 if (split->radAft()->col > -1) colMother = split->radAft()->col;
5197 if (split->radAft()->acol > -1) acolMother = split->radAft()->acol;
5198 if (split->emtAft()->col > -1) colSister = split->emtAft()->col;
5199 if (split->emtAft()->acol > -1) acolSister = split->emtAft()->acol;
5202 int colMother1, acolMother1;
5203 int colSister1, acolSister1;
5204 colMother1 = acolMother1 = colSister1 = acolSister1 = 0;
5205 if ( nEmissions == 2) {
5208 if (idMother*idDaughterNow > 0 && idMother > 0) {
5209 colMother1 = colDaughter;
5212 colSister1 = acolSister;
5216 if (idMother*idDaughterNow > 0 && idMother < 0) {
5218 acolMother1 = acolDaughter;
5219 acolSister1 = colSister;
5224 if (idMother*idDaughterNow < 0 && idMother > 0) {
5225 colMother1 = colSister;
5228 colSister1 = acolDaughter;
5230 acolMother = acolDaughter;
5231 colMother = colSister;
5235 if (idMother*idDaughterNow < 0 && idMother < 0) {
5237 acolMother1 = acolSister;
5238 acolSister1 = colDaughter;
5241 acolMother = acolSister;
5242 colMother = colDaughter;
5247 if (!useGlobalMapIF) {iMother =
event.append( idMotherNow, -41,
5248 event[iDaughter].mother1(), 0, 0, 0, colMother, acolMother,
5249 Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5250 auxevent1.append( idMotherNow, -41,
5251 event[iDaughter].mother1(), 0, 0, 0, colMother, acolMother,
5252 Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5253 auxevent2.append( idMotherNow, -41,
5254 event[iDaughter].mother1(), 0, 0, 0, colMother, acolMother,
5255 Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5258 if (!useGlobalMapIF) { iNewRecoiler =
5259 event.append( event[iRecoiler].
id(), 48,
5260 iRecoiler, iRecoiler, 0, 0, event[iRecoiler].col(),
5261 event[iRecoiler].acol(), Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5262 auxevent1.append( event[iRecoiler].
id(), 48,
5263 iRecoiler, iRecoiler, 0, 0, event[iRecoiler].col(),
5264 event[iRecoiler].acol(), Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5265 auxevent2.append( event[iRecoiler].
id(), 48,
5266 iRecoiler, iRecoiler, 0, 0, event[iRecoiler].col(),
5267 event[iRecoiler].acol(), Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5270 if (!useGlobalMapIF) iNewOther = iOldOther;
5273 if ( nEmissions == 2)
5274 iMother1 =
event.append( 0, 0, 0, 0, 0, 0, 0, 0, Vec4(0.,0.,0.,0.), 0.0,
5277 int iSister =
event.append( idSister, 43, iMother, 0, 0, 0,
5278 colSister, acolSister, Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5279 auxevent1.append( idSister, 43, iMother, 0, 0, 0,
5280 colSister, acolSister, Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5281 auxevent2.append( idSister, 43, iMother, 0, 0, 0,
5282 colSister, acolSister, Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5284 int iNewRecoiler1 = 0;
5285 if ( nEmissions == 2)
5286 iNewRecoiler1 =
event.append( event[iRecoiler].
id(), 48, iNewRecoiler,
5287 iNewRecoiler, 0, 0, event[iRecoiler].col(), event[iRecoiler].acol(),
5288 Vec4(0.,0.,0.,0.), 0.0, sqrt(pT2) );
5291 if ( nEmissions == 2)
5292 iSister1 =
event.append( 0, 0, 0, 0, 0, 0, 0, 0, Vec4(0.,0.,0.,0.), 0.0,
5296 Particle& daughter =
event[iDaughter];
5297 Particle& mother =
event[iMother];
5298 Particle& newRecoiler =
event[iNewRecoiler];
5299 Particle& sister =
event[iSister];
5300 Particle& mother1 =
event[iMother1];
5301 Particle& newRecoiler1 =
event[iNewRecoiler1];
5302 Particle& sister1 =
event[iSister1];
5305 event[iRecoiler].statusNeg();
5306 event[iRecoiler].daughters( iNewRecoiler, iNewRecoiler);
5307 auxevent1[iRecoiler].statusNeg();
5308 auxevent2[iRecoiler].statusNeg();
5309 auxevent1[iRecoiler].daughters( iNewRecoiler, iNewRecoiler);
5310 auxevent2[iRecoiler].daughters( iNewRecoiler, iNewRecoiler);
5311 if (mother.idAbs() == 21 || mother.idAbs() == 22) mother.pol(9);
5314 daughter.mothers( iMother, 0);
5315 mother.daughters( iSister, iDaughter);
5316 mother.cols( colMother, acolMother);
5317 mother.id( idMotherNow );
5318 int iBeam1Dau1 =
event[beamOff1].daughter1();
5319 int iBeam2Dau1 =
event[beamOff2].daughter1();
5320 if (iSysSelNow == 0) {
5321 event[beamOff1].daughter1( (side == 1) ? iMother : iNewOther );
5322 event[beamOff2].daughter1( (side == 2) ? iMother : iNewOther );
5325 auxevent1[iDaughter].mothers( iMother, 0);
5326 auxevent1[iMother].daughters( iSister, iDaughter);
5327 auxevent1[iMother].cols( colMother, acolMother);
5328 auxevent1[iMother].id( idMotherNow );
5329 if (iSysSelNow == 0) {
5330 auxevent1[beamOff1].daughter1( (side == 1) ? iMother : iNewOther );
5331 auxevent1[beamOff2].daughter1( (side == 2) ? iMother : iNewOther );
5334 auxevent2[iDaughter].mothers( iMother, 0);
5335 auxevent2[iMother].daughters( iSister, iDaughter);
5336 auxevent2[iMother].cols( colMother, acolMother);
5337 auxevent2[iMother].id( idMotherNow );
5338 if (iSysSelNow == 0) {
5339 auxevent2[beamOff1].daughter1( (side == 1) ? iMother : iNewOther );
5340 auxevent2[beamOff2].daughter1( (side == 2) ? iMother : iNewOther );
5343 bool doVeto =
false;
5344 bool printWarnings = (!trial || !forceMassiveMap);
5345 bool doMECreject =
false;
5348 if ( nEmissions != 2) {
5353 if (!useGlobalMapIF) {
5356 Vec4 paj_tilde(event[iDaughter].p());
5357 Vec4 pk_tilde(event[iRecoiler].p());
5358 Vec4 pTk_tilde(event[iRecoiler].px(),event[iRecoiler].py(),0.,0.);
5359 Vec4 q(paj_tilde-pk_tilde);
5360 Vec4 qpar(q+pTk_tilde);
5363 double q2 = q.m2Calc();
5364 double q2par = qpar.m2Calc();
5365 double pT2k = -pTk_tilde.m2Calc();
5366 double sjk = (1. - 1./xCS)*(q2 - m2r) + (m2e + m2s) / xCS;
5367 double zbar = (q2-sjk-m2r) / bABC(q2,sjk,m2r)
5368 *( uCS - m2r / gABC(q2,sjk,m2r)
5369 * (sjk + m2e - m2s) / (q2 - sjk - m2r));
5370 double kT2 = zbar*(1.-zbar)*sjk - (1-zbar)*m2e - zbar*m2s;
5374 infoPtr->errorMsg(
"Warning in DireSpace::branch_IF: Reject state "
5375 "with kinematically forbidden kT^2.");
5380 Vec4 pRad( ( paj_tilde - q*paj_tilde/q2par * qpar )
5381 * sqrt( (lABC(q2,sjk,m2r) - 4.*m2r*pT2k)
5382 / (lABC(q2,m2s,m2Bef) - 4.*m2Bef*pT2k))
5383 + qpar * 0.5 * (q2 + m2r - sjk) / q2par);
5386 Vec4 pSum(1e5,1e5,1e5,1e5);
5387 Vec4 pSumIn( beamAPtr->p() + beamBPtr->p() );
5393 Vec4 auxpRec1, auxpEmt1;
5394 Vec4 auxpRec2, auxpEmt2;
5396 while ( abs(pSum.px()-pSumIn.px()) > mTolErr
5397 || abs(pSum.py()-pSumIn.py()) > mTolErr ) {
5404 (
"Warning in DireSpace::branch_IF: Could not set up"
5405 " state after branching, thus reject.");
5406 physical =
false;
break;
5410 double phi_kt = (!trial) ? 2.*M_PI*rndmPtr->flat()
5411 : split->kinematics()->phi < 0. ?
5412 2.*M_PI*rndmPtr->flat() : split->kinematics()->phi;
5414 double phi_kt1 = phi_kt+DPHI_IF*M_PI;
5415 if (phi_kt1>2.*M_PI) phi_kt1 -= 2.*M_PI;
5416 double phi_kt2 = phi_kt-DPHI_IF*M_PI;
5417 if (phi_kt2<0.) phi_kt2 += 2.*M_PI;
5418 if (phi_kt1<phi_kt2) swap(phi_kt1, phi_kt2);
5421 if (split->useForBranching) { phi_kt = psp[
"phi"]; }
5428 pair<Vec4, Vec4> pTvecs = getTwoPerpendicular(pRad,pjk);
5429 Vec4 kTmom( sqrt(kT2)*sin(phi_kt)*pTvecs.first
5430 + sqrt(kT2)*cos(phi_kt)*pTvecs.second);
5433 pEmt.p( zbar * (gABC(q2,sjk,m2r)*pjk + sjk*pRad) / bABC(q2,sjk,m2r)
5434 - (m2e + kT2) / (zbar*bABC(q2,sjk,m2r))
5435 * (pRad + m2r/gABC(q2,sjk,m2r)*pjk)
5439 pRec.p(-q+pRad-pEmt);
5441 kTmom.p( sqrt(kT2)*sin(phi_kt1)*pTvecs.first
5442 + sqrt(kT2)*cos(phi_kt1)*pTvecs.second);
5443 auxpEmt1.p( zbar * (gABC(q2,sjk,m2r)*pjk + sjk*pRad) / bABC(q2,sjk,m2r)
5444 - (m2e + kT2) / (zbar*bABC(q2,sjk,m2r))
5445 * (pRad + m2r/gABC(q2,sjk,m2r)*pjk)
5447 auxpRec1.p(-q+pRad-auxpEmt1);
5449 kTmom.p( sqrt(kT2)*sin(phi_kt2)*pTvecs.first
5450 + sqrt(kT2)*cos(phi_kt2)*pTvecs.second);
5451 auxpEmt2.p( zbar * (gABC(q2,sjk,m2r)*pjk + sjk*pRad) / bABC(q2,sjk,m2r)
5452 - (m2e + kT2) / (zbar*bABC(q2,sjk,m2r))
5453 * (pRad + m2r/gABC(q2,sjk,m2r)*pjk)
5455 auxpRec2.p(-q+pRad-auxpEmt2);
5459 newRecoiler.p(pRec);
5461 auxevent1[iSister].p(auxpEmt1);
5462 auxevent1[iMother].p(pRad);
5463 auxevent1[iNewRecoiler].p(auxpRec1);
5465 auxevent2[iSister].p(auxpEmt2);
5466 auxevent2[iMother].p(pRad);
5467 auxevent2[iNewRecoiler].p(auxpRec2);
5470 vector<int> iPos(1,iSister);
5472 iPos.push_back(iNewRecoiler);
5473 pSum += newRecoiler.p();
5476 for (
int i = 0; i <
event.size(); ++i)
5477 if ( event[i].isFinal()
5478 && partonSystemsPtr->getSystemOf(i,
true) == iSysSelNow
5479 && find(iPos.begin(), iPos.end(), i) == iPos.end() )
5480 pSum += event[i].p();
5487 Vec4 paj_tilde(event[iDaughter].p());
5488 Vec4 pk_tilde(event[iRecoiler].p());
5489 Vec4 q(pk_tilde-paj_tilde);
5492 double q2 = q.m2Calc();
5494 double saj = uCS/xCS*(q2 - m2s) + (m2r+m2e) * (1-uCS/xCS);
5495 double zbar = (q2-saj-m2s) / bABC(q2,saj,m2s)
5496 *( (xCS - 1)/(xCS-uCS) - m2s / gABC(q2,saj,m2s)
5497 * (saj + m2e - m2r) / (q2 - saj - m2s));
5498 double kT2 = zbar*(1.-zbar)*saj - (1-zbar)*m2e - zbar*m2r;
5501 if (kT2 < 0. || isnan(kT2)) physical =
false;
5504 Vec4 pRec( (pk_tilde - q*pk_tilde/q2*q)
5505 *sqrt(lABC(q2,saj,m2s)/lABC(q2,m2Bef,m2s))
5506 + 0.5*(q2+m2s-saj)/q2*q );
5511 double phi_kt = (!trial) ? 2.*M_PI*rndmPtr->flat()
5512 : split->kinematics()->phi < 0. ?
5513 2.*M_PI*rndmPtr->flat() : split->kinematics()->phi;
5516 if (split->useForBranching) { phi_kt = psp[
"phi"]; }
5520 pair<Vec4, Vec4> pTvecs = getTwoPerpendicular(paj, pRec);
5521 Vec4 kTmom( sqrt(kT2)*sin(phi_kt)*pTvecs.first
5522 + sqrt(kT2)*cos(phi_kt)*pTvecs.second);
5525 Vec4 pEmt( - zbar * (gABC(q2,saj,m2s)*paj + saj*pRec) / bABC(q2,saj,m2s)
5526 + (m2e + kT2) / (zbar*bABC(q2,saj,m2s))
5527 * (pRec + m2s/gABC(q2,saj,m2s)*paj)
5531 Vec4 pRad(-q+pRec+pEmt);
5534 int iOther = getInB(iSysSelNow);
5535 if (side == 2) iOther = getInA(iSysSelNow);
5537 Vec4 pOther(event[iOther].p());
5540 RotBstMatrix toABCM;
5541 if (side == 1) toABCM.toCMframe( pRad, pOther);
5542 else toABCM.toCMframe( pOther, pRad);
5545 pRad.rotbst(toABCM);
5546 pOther.rotbst(toABCM);
5550 RotBstMatrix restoreB;
5551 restoreB.bst( pOther, event[iOther].p());
5554 pRad.rotbst(restoreB);
5555 pOther.rotbst(restoreB);
5560 newRecoiler.p(pRec);
5564 for (
int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i) {
5565 if ( event[i].isFinal()) {
5566 event[i].rotbst(toABCM);
5567 event[i].rotbst(restoreB);
5572 sister.rotbst(toABCM);
5573 sister.rotbst(restoreB);
5578 sister.m(sqrtpos(m2e));
5579 mother.m(sqrtpos(m2r));
5580 newRecoiler.m(sqrtpos(m2s));
5582 auxevent1[iSister].m(sqrtpos(m2e));
5583 auxevent1[iMother].m(sqrtpos(m2r));
5584 auxevent1[iNewRecoiler].m(sqrtpos(m2s));
5586 auxevent2[iSister].m(sqrtpos(m2e));
5587 auxevent2[iMother].m(sqrtpos(m2r));
5588 auxevent2[iNewRecoiler].m(sqrtpos(m2s));
5591 if ( !validMomentum( mother.p(), idMother, -1)
5592 || !validMomentum( sister.p(), idSister, 1)
5593 || !validMomentum( newRecoiler.p(),
event[iNewRecoiler].id(), 1))
5597 doVeto = (( canVetoEmission && userHooksPtr->doVetoISREmission(
5598 eventSizeOld, event, iSysSelNow))
5599 || ( canMergeFirst && mergingHooksPtr->doVetoEmission(event)) );
5601 double xm = 2.*mother.e() / (beamAPtr->e() + beamBPtr->e());
5604 int iOther = getInB(iSysSelNow);
5605 if (side == 2) iOther = getInA(iSysSelNow);
5606 double xAnew = (mother.mother1() == 1)
5607 ? 2.*mother.e() / (beamAPtr->e() + beamBPtr->e())
5608 : 2.*event[iOther].e() / (beamAPtr->e() + beamBPtr->e());
5609 double iAold = (mother.mother1() == 1) ? iDaughter : iOther;
5610 double iAnew = (mother.mother1() == 1) ? iMother : iOther;
5611 double xBnew = (mother.mother1() == 1)
5612 ? 2.*event[iOther].e() / (beamAPtr->e() + beamBPtr->e())
5613 : 2.*mother.e() / (beamAPtr->e() + beamBPtr->e());
5614 double iBold = (mother.mother1() == 1) ? iOther : iDaughter;
5615 double iBnew = (mother.mother1() == 1) ? iNewRecoiler : iOther;
5616 if ( hasPDF(event[iAnew].
id()) && beamAPtr->size() > 0) {
5617 double xOld = (*beamAPtr)[iSysSelNow].x();
5618 (*beamAPtr)[iSysSelNow].iPos(iAnew);
5619 (*beamAPtr)[iSysSelNow].x(xAnew);
5620 if (beamAPtr->xMax(-1) < 0.0) {
5621 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_IF: "
5622 "used up beam momentum; discard splitting.");
5626 (*beamAPtr)[iSysSelNow].iPos(iAold);
5627 (*beamAPtr)[iSysSelNow].x(xOld);
5629 if ( hasPDF(event[iBnew].
id()) && beamBPtr->size() > 0) {
5630 double xOld = (*beamBPtr)[iSysSelNow].x();
5631 (*beamBPtr)[iSysSelNow].iPos(iBnew);
5632 (*beamBPtr)[iSysSelNow].x(xBnew);
5633 if (beamBPtr->xMax(-1) < 0.0) {
5634 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_IF: "
5635 "used up beam momentum; discard splitting.");
5639 (*beamBPtr)[iSysSelNow].iPos(iBold);
5640 (*beamBPtr)[iSysSelNow].x(xOld);
5644 bool isHardSystem = partonSystemsPtr->getSystemOf(iDaughter,
true) == 0
5645 && partonSystemsPtr->getSystemOf(iRecoiler,
true) == 0;
5646 if (isHardSystem && physical && doMEcorrections
5647 && pT2 > pT2minMECs && checkSIJ(event,Q2minMECs)) {
5649 int iA = getInA(iSysSelNow);
5650 int iB = getInB(iSysSelNow);
5652 vector<int> iOut(createvector<int>(0)(0));
5653 if (useGlobalMapIF) {
5655 for (
int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
5656 iOut.push_back(partonSystemsPtr->getOut(iSysSel, iCopy - 2));
5659 bool motherHasPlusPz = (
event[iMother].pz() > 0.);
5660 if (motherHasPlusPz) {
5661 partonSystemsPtr->setInA(iSysSelNow, iMother);
5662 partonSystemsPtr->setInB(iSysSelNow, iNewOther);
5664 partonSystemsPtr->setInB(iSysSelNow, iMother);
5665 partonSystemsPtr->setInA(iSysSelNow, iNewOther);
5669 if (useGlobalMapIF) {
5671 for (
int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
5672 partonSystemsPtr->setOut(iSysSel, iCopy - 2, eventSizeOld + iCopy);
5673 partonSystemsPtr->addOut(iSysSelNow, iSister);
5674 partonSystemsPtr->replace(iSysSelNow, iRecoiler, iNewRecoiler);
5677 partonSystemsPtr->addOut(iSysSelNow, iSister);
5678 partonSystemsPtr->replace(iSysSelNow, iRecoiler, iNewRecoiler);
5681 if ( nFinalMaxMECs < 0
5682 || nFinalMaxMECs > partonSystemsPtr->sizeOut(iSysSelNow))
5683 doMECreject = applyMEC (event, split,
5684 createvector<Event>(auxevent1)(auxevent2));
5686 partonSystemsPtr->setInA(iSysSelNow, iA);
5687 partonSystemsPtr->setInB(iSysSelNow, iB);
5688 if (useGlobalMapIF) {
5689 for (
int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
5690 partonSystemsPtr->setOut(iSysSelNow, iCopy - 2, iOut[iCopy]);
5691 partonSystemsPtr->replace(iSysSelNow, iNewRecoiler, iRecoiler);
5692 partonSystemsPtr->popBackOut(iSysSelNow);
5694 partonSystemsPtr->replace(iSysSelNow, iNewRecoiler, iRecoiler);
5695 partonSystemsPtr->popBackOut(iSysSelNow);
5700 if (physical && !doVeto && !trial && !doMECreject) updateAfterIF(
5702 iDipSel, eventSizeOld, systemSizeOld, event, iDaughter, iRecoiler,
5703 iMother, iSister, iNewRecoiler, iNewOther, pT2, xm);
5712 Vec4 pa12_tilde(event[iDaughter].p());
5713 Vec4 pb_tilde(event[iRecoiler].p());
5720 if (!useGlobalMapIF) {
5728 double m2i = getMass(idMother,2);
5729 double m2j = getMass(idSister,2);
5732 double m2ai = -sai + m2a + m2i;
5735 Vec4 q( pa12_tilde - pb_tilde );
5736 double q2 = q.m2Calc();
5737 double m2jk = pT2/xa + q2*( 1. - xa/za) - m2ai;
5742 Vec4 pTk_tilde( pb_tilde.px(), pb_tilde.py(), 0., 0.);
5743 Vec4 qpar( q + pTk_tilde );
5745 uCS = za*(m2ai-m2a-m2i)/q2;
5746 xCS = uCS + xa - (pT2*za)/(q2*xa);
5749 double q2par = qpar.m2Calc();
5750 double pT2k = -pTk_tilde.m2Calc();
5751 double s_i_jk = (1. - 1./xCS)*(q2 - m2a) + (m2i + m2jk) / xCS;
5752 double zbar = (q2-s_i_jk-m2a) / bABC(q2,s_i_jk,m2a)
5753 *( uCS - m2a / gABC(q2,s_i_jk,m2a)
5754 * (s_i_jk + m2i - m2jk) / (q2 - s_i_jk - m2a));
5755 double kT2 = zbar*(1.-zbar)*s_i_jk - (1-zbar)*m2i - zbar*m2jk;
5758 if (kT2 < 0. || isnan(kT2)) physical =
false;
5761 Vec4 pa( ( pa12_tilde - 0.5*(q2-m2Bef-m2k)/q2par * qpar )
5762 * sqrt( (lABC(q2,s_i_jk,m2a) - 4.*m2a*pT2k)
5763 / (lABC(q2,m2k,m2Bef) - 4.*m2Bef*pT2k))
5764 + qpar * 0.5 * (q2 + m2a - s_i_jk) / q2par);
5769 double phi_kt = (!trial)
5770 ? ((dipEndSel->phi > 0.)
5771 ? dipEndSel->phi : 2.*M_PI*rndmPtr->flat())
5772 : ((split->kinematics()->phi > 0.)
5773 ? split->kinematics()->phi : 2.*M_PI*rndmPtr->flat());
5776 if (split->useForBranching) { phi_kt = psp[
"phi"]; }
5780 pair<Vec4, Vec4> pTvecs = getTwoPerpendicular(pa, pijk);
5781 Vec4 kTmom( sqrt(kT2)*sin(phi_kt)*pTvecs.first
5782 + sqrt(kT2)*cos(phi_kt)*pTvecs.second);
5785 Vec4 pi(zbar*(gABC(q2,s_i_jk,m2a)*pijk + s_i_jk*pa) / bABC(q2,s_i_jk,m2a)
5786 - ((1.-zbar)*s_i_jk - m2jk + m2i) / bABC(q2,s_i_jk,m2a)
5787 * (pa + m2a/gABC(q2,s_i_jk,m2a)*pijk)
5796 mother.m(sqrtpos(m2a));
5798 mother1.m(sqrtpos(m2a));
5802 sister1.m(sqrtpos(m2i));
5806 newRecoiler.m(sqrtpos(m2jk));
5813 double phiFF = (!trial)
5814 ? ((dipEndSel->phia1 > 0.)
5815 ? dipEndSel->phia1 : 2.*M_PI*rndmPtr->flat())
5816 : ((split->kinematics()->phi2 > 0.)
5817 ? split->kinematics()->phi2 : 2.*M_PI*rndmPtr->flat());
5820 if (split->useForBranching) { phiFF = psp[
"phi2"]; }
5829 double zCS = pT2/xa / ( pT2/xa - q2*xa/za);
5830 double yCS = (m2jk - m2Emt - m2Rad)
5831 / (m2jk - m2Emt - m2Rad + 2.*pai*pjk);
5836 double sij = yCS * (q2 - m2ai) + (1.-yCS)*(m2Rad+m2Emt);
5837 zbar = (q2-sij-m2ai) / bABC(q2,sij,m2ai)
5838 * (zCS - m2ai/gABC(q2,sij,m2ai)
5839 *(sij + m2Rad - m2Emt)/(q2-sij-m2ai));
5840 kT2 = zbar*(1.-zbar)*sij - (1.-zbar)*m2Rad - zbar*m2Emt;
5845 infoPtr->errorMsg(
"Warning in DireSpace::branch_IF: Reject state "
5846 "with kinematically forbidden kT^2.");
5852 if (physical && (kT2!=kT2 || abs(kT2-kT2) > 1e5) ) {
5854 infoPtr->errorMsg(
"Warning in DireSpace::branch_IF: Reject state "
5855 "with not-a-number kT^2 for branching " + name);
5864 pTvecs = getTwoPerpendicular(pai, pij);
5865 kTmom.p( sqrt(kT2)*sin(phiFF)*pTvecs.first
5866 + sqrt(kT2)*cos(phiFF)*pTvecs.second);
5869 Vec4 pj( zbar * (gABC(q2,sij,m2ai)*pij - sij*pai) / bABC(q2,sij,m2ai)
5870 + (m2Rad+kT2) / (zbar*bABC(q2,sij,m2ai))
5871 * (pai - m2ai/gABC(q2,sij,m2ai)*pij)
5879 sister.m(sqrtpos(m2j));
5881 newRecoiler1.m(sqrtpos(m2k));
5884 if ( !validMomentum( mother1.p(), idMother, -1)
5885 || !validMomentum( sister.p(), idSister, 1)
5886 || !validMomentum( sister1.p(), idMother, 1)
5887 || !validMomentum( newRecoiler1.p(),
event[iNewRecoiler1].id(), 1))
5892 double saix(2.*pa*pi), sakx(2.*pa*pk), sajx(2.*pa*pj), sikx(2.*pi*pk),
5893 sjkx(2.*pj*pk), sijx(2.*pi*pj);
5894 double pptt = (sajx-sijx)*(sakx-sikx)/(saix+sajx+sakx);
5895 double ssaaii = saix;
5896 double zzaa = -q2tot/ ( saix + sajx + sakx );
5897 double xxaa = (sakx-sikx) / ( saix + sajx + sakx );
5899 (abs(pptt-pT2) > 1e-5 || abs(ssaaii-sai) > 1e-5 ||
5900 abs(zzaa-za) > 1e-5 || abs(xxaa-xa) > 1e-5) ){
5901 cout <<
"Error in branch_IF: Invariant masses after branching do "
5902 <<
"not match chosen values." << endl;
5908 <<
" xa " << xa << endl;
5909 cout <<
"Generated: "
5910 <<
" Q2 " << saix+sajx+sakx-sijx-sikx-sjkx
5912 <<
" sai " << ssaaii
5914 <<
" xa " << xxaa << endl;
5920 int iOther = getInB(iSysSelNow);
5921 if (side == 2) iOther = getInA(iSysSelNow);
5922 double xAnew = (mother.mother1() == 1)
5923 ? 2.*mother1.e() / (beamAPtr->e() + beamBPtr->e())
5924 : 2.*event[iOther].e() / (beamAPtr->e() + beamBPtr->e());
5925 double iAold = (mother.mother1() == 1) ? iDaughter : iOther;
5926 double iAnew = (mother.mother1() == 1) ? iMother1 : iOther;
5927 double xBnew = (mother.mother1() == 1)
5928 ? 2.*event[iOther].e() / (beamAPtr->e() + beamBPtr->e())
5929 : 2.*mother1.e() / (beamAPtr->e() + beamBPtr->e());
5930 double iBold = (mother.mother1() == 1) ? iOther : iDaughter;
5931 double iBnew = (mother.mother1() == 1) ? iOther : iMother1;
5932 if ( hasPDF(event[iAnew].
id()) && beamAPtr->size() > 0) {
5933 double xOld = (*beamAPtr)[iSysSelNow].x();
5934 (*beamAPtr)[iSysSelNow].iPos(iAnew);
5935 (*beamAPtr)[iSysSelNow].x(xAnew);
5936 if (beamAPtr->xMax(-1) < 0.0) {
5937 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_II: "
5938 "used up beam momentum; discard splitting.");
5942 (*beamAPtr)[iSysSelNow].iPos(iAold);
5943 (*beamAPtr)[iSysSelNow].x(xOld);
5945 if ( hasPDF(event[iBnew].
id()) && beamBPtr->size() > 0) {
5946 double xOld = (*beamBPtr)[iSysSelNow].x();
5947 (*beamBPtr)[iSysSelNow].iPos(iBnew);
5948 (*beamBPtr)[iSysSelNow].x(xBnew);
5949 if (beamBPtr->xMax(-1) < 0.0) {
5950 if (!trial) infoPtr->errorMsg(
"Warning in DireSpace::branch_II: "
5951 "used up beam momentum; discard splitting.");
5955 (*beamBPtr)[iSysSelNow].iPos(iBold);
5956 (*beamBPtr)[iSysSelNow].x(xOld);
5959 double xm = 2.*mother1.e() / (beamAPtr->e() + beamBPtr->e());
5962 if (physical && !trial) updateAfterIF( iSysSelNow, side, iDipSel,
5963 eventSizeOld, systemSizeOld, event, iDaughter, iRecoiler, iMother,
5964 iSister, iNewRecoiler, iNewOther, pT2, xm);
5969 newRecoiler.status(-49);
5970 newRecoiler.statusNeg();
5971 newRecoiler.daughters( iNewRecoiler1, iNewRecoiler1);
5972 mother1.id( idMother);
5973 mother1.status(-41);
5974 mother1.cols(colMother1, acolMother1);
5975 mother1.daughters( iSister1, iMother);
5976 mother1.mothers( mother.mother1(), mother.mother2());
5977 mother.mothers( iMother1, 0);
5978 sister1.id(idMother);
5980 sister1.mothers(iMother1,0);
5981 sister1.cols(colSister1, acolSister1);
5982 sister1.scale(sqrt(pT2));
5984 if (iSysSelNow == 0) {
5985 event[beamOff1].daughter1( (side == 1) ? iMother1 : iNewOther );
5986 event[beamOff2].daughter1( (side == 2) ? iMother1 : iNewOther );
5990 if (physical && !trial) updateAfterIF( iSysSelNow, side, iDipSel, 0, 0,
5991 event, iMother, iNewRecoiler, iMother1, iSister1, iNewRecoiler1,
5992 iNewOther, pT2, xm);
5999 double m2i = getMass(idMother,2);
6000 double m2j = getMass(idSister,2);
6002 double m2ai = -sai + m2a + m2i;
6008 Vec4 q(pb_tilde-pa12_tilde);
6009 double q2 = q.m2Calc();
6011 double m2jk = pT2/xa + q2*( 1. - xa/za) - m2ai;
6012 uCS = za*(m2ai-m2a-m2i)/q2;
6013 xCS = uCS + xa - (pT2*za)/(q2*xa);
6016 double zbar = (q2-m2ai-m2jk) / bABC(q2,m2ai,m2jk)
6017 *( (xCS - 1)/(xCS-uCS) - m2jk / gABC(q2,m2ai,m2jk)
6018 * (m2ai + m2i - m2a) / (q2 - m2ai - m2jk));
6019 double kT2 = zbar*(1.-zbar)*m2ai - (1-zbar)*m2i - zbar*m2a;
6022 if (kT2 < 0. || isnan(kT2)) physical =
false;
6025 Vec4 pjk( (pb_tilde - q*pb_tilde/q2*q)
6026 *sqrt(lABC(q2,m2ai,m2jk)/lABC(q2,m2Bef,m2s))
6027 + 0.5*(q2+m2jk-m2ai)/q2*q );
6032 double phi_kt = (!trial)
6033 ? ((dipEndSel->phi > 0.)
6034 ? dipEndSel->phi : 2.*M_PI*rndmPtr->flat())
6035 : ((split->kinematics()->phi > 0.)
6036 ? split->kinematics()->phi : 2.*M_PI*rndmPtr->flat());
6039 if (split->useForBranching) { phi_kt = psp[
"phi"]; }
6043 pair<Vec4, Vec4> pTvecs = getTwoPerpendicular(pai, pjk);
6044 Vec4 kTmom( sqrt(kT2)*sin(phi_kt)*pTvecs.first
6045 + sqrt(kT2)*cos(phi_kt)*pTvecs.second);
6048 Vec4 pi( - zbar *(gABC(q2,m2ai,m2jk)*pai + m2ai*pjk) / bABC(q2,m2ai,m2jk)
6049 + ( (1.-zbar)*m2ai + m2i - m2a) / bABC(q2,m2ai,m2jk)
6050 * (pjk + m2jk/gABC(q2,m2ai,m2jk)*pai)
6057 int iOther = getInB(iSysSelNow);
6058 if (side == 2) iOther = getInA(iSysSelNow);
6059 Vec4 pOther(event[iOther].p());
6062 RotBstMatrix toABCM;
6063 if (side == 1) toABCM.toCMframe( pa, pOther);
6064 else toABCM.toCMframe( pOther, pa);
6068 pOther.rotbst(toABCM);
6072 RotBstMatrix restoreB;
6073 restoreB.bst( pOther, event[iOther].p());
6076 pa.rotbst(restoreB);
6077 pOther.rotbst(restoreB);
6081 pi.rotbst(restoreB);
6083 pjk.rotbst(restoreB);
6087 sister1.m(sqrtpos(m2i));
6090 mother.m(sqrtpos(m2a));
6092 mother1.m(sqrtpos(m2a));
6097 for (
int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i) {
6100 if ( i == iSister1 || i == iNewRecoiler)
continue;
6101 if ( event[i].isFinal()) {
6102 event[i].rotbst(toABCM);
6103 event[i].rotbst(restoreB);
6112 double phiFF = (!trial)
6113 ? ((dipEndSel->phia1 > 0.)
6114 ? dipEndSel->phia1 : 2.*M_PI*rndmPtr->flat())
6115 : ((split->kinematics()->phi2 > 0.)
6116 ? split->kinematics()->phi2 : 2.*M_PI*rndmPtr->flat());
6119 if (split->useForBranching) { phiFF = psp[
"phi2"]; }
6128 double zCS = pT2/xa / ( pT2/xa - q2*xa/za);
6129 double yCS = (m2jk - m2Emt - m2Rad)
6130 / (m2jk - m2Emt - m2Rad + 2.*pai*pjk);
6135 double sij = yCS * (q2 - m2ai) + (1.-yCS)*(m2Rad+m2Emt);
6136 zbar = (q2-sij-m2ai) / bABC(q2,sij,m2ai)
6137 * (zCS - m2ai/gABC(q2,sij,m2ai)
6138 *(sij + m2Rad - m2Emt)/(q2-sij-m2ai));
6139 kT2 = zbar*(1.-zbar)*sij - (1.-zbar)*m2Rad - zbar*m2Emt;
6144 infoPtr->errorMsg(
"Warning in DireSpace::branch_IF: Reject state "
6145 "with kinematically forbidden kT^2.");
6151 if (physical && (kT2!=kT2 || abs(kT2-kT2) > 1e5) ) {
6153 infoPtr->errorMsg(
"Warning in DireSpace::branch_IF: Reject state "
6154 "with not-a-number kT^2 for branching " + name);
6163 pTvecs = getTwoPerpendicular(pai, pij);
6164 kTmom.p( sqrt(kT2)*sin(phiFF)*pTvecs.first
6165 + sqrt(kT2)*cos(phiFF)*pTvecs.second);
6168 Vec4 pj( zbar * (gABC(q2,sij,m2ai)*pij - sij*pai) / bABC(q2,sij,m2ai)
6169 + (m2Rad+kT2) / (zbar*bABC(q2,sij,m2ai))
6170 * (pai - m2ai/gABC(q2,sij,m2ai)*pij)
6178 sister.m(sqrtpos(m2j));
6180 newRecoiler1.m(sqrtpos(m2k));
6183 if ( !validMomentum( mother1.p(), idMother, -1)
6184 || !validMomentum( sister.p(), idSister, 1)
6185 || !validMomentum( sister1.p(), idMother, 1)
6186 || !validMomentum( newRecoiler1.p(),
event[iNewRecoiler1].id(), 1))
6191 double saix(2.*pa*pi), sakx(2.*pa*pk), sajx(2.*pa*pj), sikx(2.*pi*pk),
6192 sjkx(2.*pj*pk), sijx(2.*pi*pj);
6193 double pptt = (sajx-sijx)*(sakx-sikx)/(saix+sajx+sakx);
6194 double ssaaii = saix;
6195 double zzaa = -q2tot/ ( saix + sajx + sakx );
6196 double xxaa = (sakx-sikx) / ( saix + sajx + sakx );
6198 (abs(pptt-pT2) > 1e-5 || abs(ssaaii-sai) > 1e-5 ||
6199 abs(zzaa-za) > 1e-5 || abs(xxaa-xa) > 1e-5) ){
6200 cout <<
"Error in branch_IF: Invariant masses after branching do "
6201 <<
"not match chosen values." << endl;
6207 <<
" xa " << xa << endl;
6208 cout <<
"Generated: "
6209 <<
" Q2 " << saix+sajx+sakx-sijx-sikx-sjkx
6211 <<
" sai " << ssaaii
6213 <<
" xa " << xxaa << endl;
6218 double xm = 2.*mother1.e() / (beamAPtr->e() + beamBPtr->e());
6221 if (physical && !trial) updateAfterIF( iSysSelNow, side, iDipSel,
6222 eventSizeOld, systemSizeOld, event, iDaughter, iRecoiler, iMother,
6223 iSister, iNewRecoiler, iNewOther, pT2, xm);
6228 newRecoiler.status(-49);
6229 newRecoiler.statusNeg();
6230 newRecoiler.daughters( iNewRecoiler1, iNewRecoiler1);
6231 mother1.id( idMother);
6232 mother1.status(-41);
6233 mother1.cols(colMother1, acolMother1);
6234 mother1.daughters( iSister1, iMother);
6235 mother1.mothers( mother.mother1(), mother.mother2());
6236 mother.mothers( iMother1, 0);
6237 sister1.id(idMother);
6239 sister1.mothers(iMother1,0);
6240 sister1.cols(colSister1, acolSister1);
6241 sister1.scale(sqrt(pT2));
6243 if (iSysSelNow == 0) {
6244 event[beamOff1].daughter1( (side == 1) ? iMother1 : iNewOther );
6245 event[beamOff2].daughter1( (side == 2) ? iMother1 : iNewOther );
6249 if (physical && !trial) updateAfterIF( iSysSelNow, side, iDipSel, 0, 0,
6250 event, iMother, iNewRecoiler, iMother1, iSister1, iNewRecoiler1,
6251 iNewOther, pT2, xm);
6259 physical = physical && !doVeto;
6263 if (iSysSelNow > 0) {
6264 if (side == 1)
event[beamOff1].daughter1(0);
6265 if (side == 2)
event[beamOff2].daughter1(0);
6270 if ( physical && !trial && !doMECreject && !validMotherDaughter(event)) {
6272 infoPtr->errorMsg(
"Error in DireSpace::branch_IF: Mother-daughter "
6273 "relations after branching not valid.");
6278 if (iSysSelNow > 0) {
6279 if (side == 1)
event[beamOff1].daughter1(iBeam1Dau1);
6280 if (side == 2)
event[beamOff2].daughter1(iBeam2Dau1);
6284 if ( !physical || doMECreject) {
6285 event.popBack( event.size() - eventSizeOld);
6287 if (iSysSelNow == 0) {
6288 event[beamOff1].daughter1( ev1Dau1V);
6289 event[beamOff2].daughter1( ev2Dau1V);
6292 if (useGlobalMapIF) {
6293 for (
int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
6294 int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
6295 event[iOldCopy].status( statusV[iCopy]);
6296 event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
6297 event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
6300 event[iDaughter].status( iDauStatusV);
6301 event[iDaughter].mothers(iDauMot1V, iDauMot2V);
6302 event[iDaughter].daughters(iDauDau1V, iDauDau2V);
6303 event[iRecoiler].status( iRecStatusV);
6304 event[iRecoiler].mothers( iRecMot1V, iRecMot2V);
6305 event[iRecoiler].daughters( iRecDau1V, iRecDau2V);
6313 for ( unordered_map<
string, multimap<double,double> >::iterator
6314 it = rejectProbability.begin(); it != rejectProbability.end(); ++it){
6315 weights->eraseAcceptWeight(pT2, it->first);
6316 weights->eraseRejectWeight(pT2, it->first);
6320 if (!trial && doMECreject) {
6321 weights->calcWeight(pT2,
false,
true);
6324 for ( unordered_map<
string, multimap<double,double> >::iterator
6325 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
6327 for ( unordered_map<
string, map<double,double> >::iterator
6328 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
6336 if (trial) split->storePosAfter(
6337 (nEmissions < 2) ? iMother : iMother1,
6338 (nEmissions < 2) ? iNewRecoiler : iNewRecoiler1,
6339 iSister, (nEmissions < 2) ? 0 : iSister1);
6344 weights->calcWeight(pT2);
6348 direInfoPtr->updateSoftPos( iDaughter, iSister );
6349 direInfoPtr->updateSoftPosIfMatch( iRecoiler, iNewRecoiler );
6350 if (nEmissions > 1) {
6351 direInfoPtr->updateSoftPosIfMatch( iNewRecoiler, iNewRecoiler1 );
6352 direInfoPtr->addSoftPos( iSister1 );
6357 for ( unordered_map<
string, multimap<double,double> >::iterator
6358 it = rejectProbability.begin(); it != rejectProbability.end(); ++it )
6360 for ( unordered_map<
string, map<double,double> >::iterator
6361 it = acceptProbability.begin(); it != acceptProbability.end(); ++it )
6372 void DireSpace::updateAfterIF(
int iSysSelNow,
int sideNow,
int iDipSelNow,
6373 int eventSizeOldNow,
int systemSizeOldNow,
Event& event,
int iDaughter,
6374 int iRecoiler,
int iMother,
int iSister,
int iNewRecoiler,
int iNewOther,
6375 double pT2,
double xNew) {
6378 if ( nProposedPT.find(iSysSelNow) != nProposedPT.end() )
6379 ++nProposedPT[iSysSelNow];
6381 int idMother =
event[iMother].id();
6382 int idDaughterNow =
event[iDaughter].id();
6383 bool motherHasPlusPz = (
event[iMother].pz() > 0.);
6386 if ( direInfoPtr->isRes(iDaughter) &&
6387 event[iMother].id() !=
event[iDaughter].id() )
6388 direInfoPtr->removeResPos(iDaughter);
6389 if ( particleDataPtr->isResonance(event[iMother].id()) ) {
6390 if ( direInfoPtr->isRes(iDaughter) )
6391 direInfoPtr->updateResPos(iDaughter,iMother);
6393 if ( direInfoPtr->isRes(iRecoiler) )
6394 direInfoPtr->updateResPos(iRecoiler,iNewRecoiler);
6395 if ( particleDataPtr->isResonance(event[iSister].id()) )
6396 direInfoPtr->addResPos(iSister);
6399 if (motherHasPlusPz) {
6400 partonSystemsPtr->setInA(iSysSelNow, iMother);
6401 partonSystemsPtr->setInB(iSysSelNow, iNewOther);
6403 partonSystemsPtr->setInB(iSysSelNow, iMother);
6404 partonSystemsPtr->setInA(iSysSelNow, iNewOther);
6408 if (useGlobalMapIF) {
6410 for (
int iCopy = 2; iCopy < systemSizeOldNow; ++iCopy) {
6411 int iOut = partonSystemsPtr->getOut(iSysSelNow, iCopy - 2);
6413 direInfoPtr->updateResPosIfMatch ( iOut, eventSizeOldNow + iCopy);
6414 direInfoPtr->updateSoftPosIfMatch( iOut, eventSizeOldNow + iCopy);
6415 partonSystemsPtr->setOut(iSysSelNow, iCopy - 2, eventSizeOldNow + iCopy);
6417 partonSystemsPtr->addOut(iSysSelNow, iSister);
6418 partonSystemsPtr->replace(iSysSelNow, iRecoiler, iNewRecoiler);
6421 partonSystemsPtr->addOut(iSysSelNow, iSister);
6422 partonSystemsPtr->replace(iSysSelNow, iRecoiler, iNewRecoiler);
6426 int iA = getInA(iSysSelNow);
6427 int iB = getInB(iSysSelNow);
6428 double shat = (
event[iA].p() +
event[iB].p()).m2Calc();
6429 partonSystemsPtr->setSHat(iSysSelNow, shat);
6432 dipEndSel = &dipEnd[iDipSelNow];
6435 for (
int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
6436 if ( dipEnd[iDip].system == iSysSelNow) {
6437 if (abs(dipEnd[iDip].side) == sideNow) {
6438 dipEnd[iDip].iRadiator = iMother;
6439 dipEnd[iDip].iRecoiler = iNewRecoiler;
6440 if (dipEnd[iDip].colType != 0)
6441 dipEnd[iDip].colType =
event[iMother].colType();
6445 dipEnd[iDip].iRadiator = iNewRecoiler;
6446 dipEnd[iDip].iRecoiler = iMother;
6447 dipEnd[iDip].MEtype = 0;
6452 BeamParticle& beamNow = (sideNow == 1) ? *beamAPtr : *beamBPtr;
6453 beamNow[iSysSelNow].update( iMother, idMother, xNew);
6455 if (idMother != idDaughterNow) {
6456 pdfScale2 = (useFixedFacScale) ? fixedFacScale2 : factorMultFac * pT2;
6457 pdfScale2 = max(pdfScale2, pT2min);
6458 beamNow.xfISR( iSysSelNow, idMother, xNew, pdfScale2);
6459 beamNow.pickValSeaComp();
6461 BeamParticle& beamRec = (sideNow == 1) ? *beamBPtr : *beamAPtr;
6462 beamRec[iSysSelNow].iPos( iNewOther);
6465 update(iSysSelNow,event);
6474 pair <Event, pair<int,int> > DireSpace::clustered_internal(
const Event& state,
6475 int iRad,
int iEmt,
int iRecAft,
string name ) {
6478 int radType = state[iRad].isFinal() ? 1 : -1;
6479 int recType = state[iRecAft].isFinal() ? 1 : -1;
6483 NewEvent.init(
"(hard process-modified)", particleDataPtr);
6486 for (
int i = 0; i < state.size(); ++i) {
6487 if ( i == iRad || i == iRecAft || i == iEmt )
continue;
6488 NewEvent.append( state[i] );
6492 for (
int i = 0; i < state.sizeJunction(); ++i)
6493 NewEvent.appendJunction( state.getJunction(i) );
6495 double mu = infoPtr->QFac();
6497 NewEvent.saveSize();
6498 NewEvent.saveJunctionSize();
6500 NewEvent.scaleSecond(mu);
6504 Particle RecBefore = Particle( state[iRecAft] );
6505 RecBefore.setEvtPtr(&NewEvent);
6506 RecBefore.daughters(0,0);
6508 int radID = splits[name]->radBefID(state[iRad].
id(), state[iEmt].
id());
6509 int recID = state[iRecAft].id();
6510 Particle RadBefore = Particle( state[iRad] );
6511 RadBefore.setEvtPtr(&NewEvent);
6512 RadBefore.id(radID);
6513 RadBefore.daughters(0,0);
6515 RadBefore.cols(RecBefore.acol(),RecBefore.col());
6518 if ( particleDataPtr->isResonance(radID) && radType == 1)
6519 RadBefore.status(state[iRad].status());
6522 double radMass = particleDataPtr->m0(radID);
6523 double recMass = particleDataPtr->m0(recID);
6524 if (radType == 1 ) RadBefore.m(radMass);
6525 else RadBefore.m(0.0);
6526 if (recType == 1 ) RecBefore.m(recMass);
6527 else RecBefore.m(0.0);
6530 bool isClustered =
false;
6531 if ( state[iRecAft].isFinal())
6532 isClustered = cluster_IF(state,iRad,iEmt,iRecAft,radID,RadBefore,
6533 RecBefore,NewEvent);
6535 isClustered = cluster_II(state,iRad,iEmt,iRecAft,radID,RadBefore,
6536 RecBefore,NewEvent);
6539 if (!isClustered) { NewEvent.clear();
return make_pair
6540 (NewEvent, make_pair(0,0));}
6543 RecBefore.scale(mu);
6544 RadBefore.scale(mu);
6547 NewEvent.append(RecBefore);
6550 pair<int,int> cols = splits[name]->radBefCols( state[iRad].col(),
6551 state[iRad].acol(), state[iEmt].col(), state[iEmt].acol());
6552 RadBefore.cols( cols.first, cols.second );
6556 outState.init(
"(hard process-modified)", particleDataPtr);
6560 for (
int i = 0; i < 3; ++i)
6561 outState.append( NewEvent[i] );
6563 for (
int i = 0; i < state.sizeJunction(); ++i)
6564 outState.appendJunction( state.getJunction(i) );
6566 outState.saveSize();
6567 outState.saveJunctionSize();
6569 outState.scaleSecond(mu);
6570 bool radAppended =
false;
6571 bool recAppended =
false;
6572 int size = int(outState.size());
6574 int radPos(0), recPos(0);
6577 if ( RecBefore.mother1() == 1) {
6578 recPos = outState.append( RecBefore );
6580 }
else if ( RadBefore.mother1() == 1 ) {
6581 radPos = outState.append( RadBefore );
6586 for(
int i=0; i < int(state.size()); ++i)
6587 if (state[i].mother1() == 1) in1 =i;
6588 outState.append( state[in1] );
6592 if ( RecBefore.mother1() == 2) {
6593 recPos = outState.append( RecBefore );
6595 }
else if ( RadBefore.mother1() == 2 ) {
6596 radPos = outState.append( RadBefore );
6601 for(
int i=0; i < int(state.size()); ++i)
6602 if (state[i].mother1() == 2) in2 =i;
6604 outState.append( state[in2] );
6609 if (!recAppended && !RecBefore.isFinal()) {
6611 recPos = outState.append( RecBefore);
6614 if (!radAppended && !RadBefore.isFinal()) {
6616 radPos = outState.append( RadBefore);
6620 outState[3].status(-21);
6621 outState[4].status(-21);
6625 for (
int i = 0; i < int(NewEvent.size()-1); ++i) {
6626 if (NewEvent[i].status() != -22)
continue;
6627 if ( NewEvent[i].daughter1() == NewEvent[i].daughter2()
6628 && NewEvent[i].daughter1() > 0)
continue;
6629 outState.append( NewEvent[i] );
6632 for (
int i = 0; i < int(NewEvent.size()-1); ++i)
6633 if (NewEvent[i].status() == 22) outState.append( NewEvent[i] );
6635 if (!radAppended && RadBefore.statusAbs() == 22)
6636 radPos = outState.append(RadBefore);
6638 recPos = outState.append(RecBefore);
6639 if (!radAppended && RadBefore.statusAbs() != 22)
6640 radPos = outState.append(RadBefore);
6642 for(
int i = 0; i < int(NewEvent.size()-1); ++i)
6643 if ( NewEvent[i].status() != 22
6644 && NewEvent[i].colType() != 0
6645 && NewEvent[i].isFinal()) {
6646 outState.append( NewEvent[i] );
6648 int status = particleDataPtr->isResonance(NewEvent[i].
id()) ? 22 : 23;
6649 outState.back().status(status);
6650 outState.back().mother1(3);
6651 outState.back().mother2(4);
6654 for(
int i = 0; i < int(NewEvent.size()-1); ++i)
6655 if ( NewEvent[i].status() != 22
6656 && NewEvent[i].colType() == 0
6657 && NewEvent[i].isFinal() ) {
6658 outState.append( NewEvent[i]);
6659 int status = particleDataPtr->isResonance(NewEvent[i].
id()) ? 22 : 23;
6660 outState.back().status(status);
6661 outState.back().mother1(3);
6662 outState.back().mother2(4);
6666 vector<int> PosIntermediate;
6667 vector<int> PosDaughter1;
6668 vector<int> PosDaughter2;
6669 for(
int i=0; i < int(outState.size()); ++i)
6670 if (outState[i].status() == -22) {
6671 PosIntermediate.push_back(i);
6672 int d1 = outState[i].daughter1();
6673 int d2 = outState[i].daughter2();
6675 int daughter1 = FindParticle( state[d1], outState);
6676 int daughter2 = FindParticle( state[d2], outState);
6681 PosDaughter1.push_back( daughter1);
6684 while(!outState[daughter1].isFinal() ) daughter1++;
6685 PosDaughter1.push_back( daughter1);
6688 PosDaughter2.push_back( daughter2);
6690 daughter2 = outState.size()-1;
6691 while(!outState[daughter2].isFinal() ) daughter2--;
6692 PosDaughter2.push_back( daughter2);
6696 for(
int i=0; i < int(PosIntermediate.size()); ++i) {
6697 outState[PosIntermediate[i]].daughters(PosDaughter1[i],PosDaughter2[i]);
6698 outState[PosDaughter1[i]].mother1(PosIntermediate[i]);
6699 outState[PosDaughter2[i]].mother1(PosIntermediate[i]);
6700 outState[PosDaughter1[i]].mother2(0);
6701 outState[PosDaughter2[i]].mother2(0);
6705 int minParFinal = int(outState.size());
6706 int maxParFinal = 0;
6707 for(
int i=0; i < int(outState.size()); ++i)
6708 if (outState[i].mother1() == 3 && outState[i].mother2() == 4) {
6709 minParFinal = min(i,minParFinal);
6710 maxParFinal = max(i,maxParFinal);
6713 if (minParFinal == maxParFinal) maxParFinal = 0;
6714 outState[3].daughters(minParFinal,maxParFinal);
6715 outState[4].daughters(minParFinal,maxParFinal);
6718 outState.saveSize();
6719 outState.saveJunctionSize();
6733 if ( radType == -1 && state[iRad].colType() == 1) {
6735 for(
int i=0; i < int(state.size()); ++i)
6736 if ( i != iRad && i != iEmt && i != iRecAft
6737 && state[i].status() == -22
6738 && state[i].col() == state[iRad].col() )
6740 }
else if ( radType == -1 && state[iRad].colType() == -1) {
6742 for(
int i=0; i < int(state.size()); ++i)
6743 if ( i != iRad && i != iEmt && i != iRecAft
6744 && state[i].status() == -22
6745 && state[i].acol() == state[iRad].acol() )
6747 }
else if ( radType == 1 && state[iRad].colType() == 1) {
6749 for(
int i=0; i < int(state.size()); ++i)
6750 if ( i != iRad && i != iEmt && i != iRecAft
6751 && state[i].status() == -22
6752 && state[i].col() == state[iRad].col() )
6754 }
else if ( radType == 1 && state[iRad].colType() == -1) {
6756 for(
int i=0; i < int(state.size()); ++i)
6757 if ( i != iRad && i != iEmt && i != iRecAft
6758 && state[i].status() == -22
6759 && state[i].acol() == state[iRad].acol() )
6765 int iColResNow = FindParticle( state[iColRes], outState);
6768 int radCol = outState[radPos].col();
6769 int radAcl = outState[radPos].acol();
6771 int resCol = outState[iColResNow].col();
6772 int resAcl = outState[iColResNow].acol();
6774 bool matchesRes = (radCol > 0
6775 && ( radCol == resCol || radCol == resAcl))
6777 && ( radAcl == resCol || radAcl == resAcl));
6781 if (!matchesRes && iColResNow > 0) {
6782 if ( radType == -1 && outState[radPos].colType() == 1)
6783 outState[iColResNow].col(radCol);
6784 else if ( radType ==-1 && outState[radPos].colType() ==-1)
6785 outState[iColResNow].acol(radAcl);
6786 else if ( radType == 1 && outState[radPos].colType() == 1)
6787 outState[iColResNow].col(radCol);
6788 else if ( radType == 1 && outState[radPos].colType() ==-1)
6789 outState[iColResNow].acol(radAcl);
6795 if (!matchesRes && iColResNow > 0 && iColRes != iColResNow)
6796 outState[radPos].mother1(iColResNow);
6801 int iDof1 = outState[3].mother1() == 1 ? 3 : 4;
6802 int iDof2 = outState[3].mother1() == 2 ? 3 : 4;
6804 if (outState[1].idAbs() == 11 && outState[iDof1].
id() == 22 ) {
6805 outState[1].id(outState[iDof1].
id());
6807 if (outState[2].idAbs() == 11 && outState[iDof2].
id() == 22 ) {
6808 outState[2].id(outState[iDof2].
id());
6810 if (outState[1].idAbs() == 11 && outState[iDof1].colType() != 0) {
6811 outState[1].id(2212);
6813 if (outState[2].idAbs() == 11 && outState[iDof2].colType() != 0) {
6814 outState[2].id(2212);
6818 for (
int i = 0; i < outState.size(); ++i) {
6819 if ( outState[i].status() == 23 && particleDataPtr->
6820 isResonance(outState[i].
id())) outState[i].status(22);
6824 if (!validEvent( outState,
true )) { outState.clear(); }
6827 return make_pair(outState, make_pair(radPos, recPos));
6833 bool DireSpace::cluster_II(
const Event& state,
6834 int iRad,
int iEmt,
int iRecAft,
int idRadBef, Particle& radBef,
6835 Particle& recBef,
Event& partial ) {
6837 if (
false) cout << idRadBef;
6840 double pT2 = pT2_II(state[iRad], state[iEmt], state[iRecAft]);
6841 double Q2 = 2.*state[iRad].p()*state[iRecAft].p()
6842 - 2.*state[iRad].p()*state[iEmt].p()
6843 - 2.*state[iEmt].p()*state[iRecAft].p();
6844 double z = z_II(state[iRad], state[iEmt], state[iRecAft]);
6846 double kappa2 = pT2 / Q2;
6847 double xCS = (z*(1-z)- kappa2)/(1-z);
6850 double m2Bef = 0.0, m2r = 0.0;
6851 double m2e = state[iEmt].p().m2Calc();
6852 double m2s = state[iRecAft].p().m2Calc();
6855 double xNew = 2.*state[iRad].e()/state[0].m();
6856 double xOld = xNew*xCS;
6858 double xMin = (particleDataPtr->colType(idRadBef) != 0) ? xOld : 0.;
6860 if ( !inAllowedPhasespace( 1, z, pT2, Q2, xMin, -2, m2Bef, m2r, m2s, m2e) ) {
6865 Vec4 q(state[iRad].p() - state[iEmt].p() + state[iRecAft].p());
6866 double q2 = q.m2Calc();
6867 double sab = (state[iRad].p() + state[iRecAft].p()).m2Calc();
6869 Vec4 pRad = ( state[iRad].p() - m2r/gABC(sab,m2r,m2s)*state[iRecAft].p())
6870 *sqrt(lABC(q2,m2Bef,m2s)/lABC(sab,m2r,m2s))
6871 + m2Bef / gABC(q2,m2Bef,m2s)*state[iRecAft].p();
6874 recBef.p( state[iRecAft].p() );
6880 Vec4 kTilde(radBef.p() + recBef.p());
6881 Vec4 k(state[iRad].p() + state[iRecAft].p() - state[iEmt].p());
6882 Vec4 kSum(kTilde + k);
6883 for (
int i = 0; i < partial.size(); ++i) {
6884 if ( !partial[i].isFinal() && partial[i].statusAbs() != 22 )
continue;
6885 Vec4 pIn = partial[i].p();
6886 double kSum2 = kSum.m2Calc();
6887 double k2 = k.m2Calc();
6889 double kSumXp = kSum*pIn;
6890 Vec4 res = pIn - kSum * 2.0*( kSumXp / kSum2 ) + kTilde * 2.0 *( kXp/k2);
6901 bool DireSpace::cluster_IF(
const Event& state,
6902 int iRad,
int iEmt,
int iRecAft,
int idRadBef, Particle& radBef,
6903 Particle& recBef,
Event& partial ) {
6905 if (
false) cout << idRadBef;
6908 double pT2 = pT2_IF(state[iRad], state[iEmt], state[iRecAft]);
6909 double z = z_IF(state[iRad], state[iEmt], state[iRecAft]);
6911 int side = (state[iRad].pz() > 0.) ? 1 : -1;
6914 double m2Bef = 0.0, m2r = 0.0;
6916 double m2e = state[iEmt].m2Calc();
6917 double m2s = state[iRecAft].p().m2Calc();
6920 double Q2 = 2.*state[iRad].p()*state[iEmt].p()
6921 + 2.*state[iRad].p()*state[iRecAft].p()
6922 - 2.*state[iRecAft].p()*state[iEmt].p();
6923 double xNew = 2.*state[iRad].e()/state[0].m();
6924 double xOld = xNew*xCS;
6927 double xMin = (particleDataPtr->colType(idRadBef) != 0) ? xOld : 0.;
6929 if ( !inAllowedPhasespace( 1, z, pT2, Q2, xMin, 2, m2Bef, m2r, m2s, m2e) )
6932 Vec4 pRadBef, pRecBef;
6935 if (!useGlobalMapIF) {
6938 Vec4 paj(state[iRad].p()-state[iEmt].p());
6939 Vec4 pk(state[iRecAft].p());
6940 Vec4 pTk(pk.px()+state[iEmt].px(),pk.py()+state[iEmt].py(),0.,0.);
6945 double q2 = q.m2Calc();
6946 double q2par = qpar.m2Calc();
6947 double pT2k = -pTk.m2Calc();
6949 double sjk = ( state[iRecAft].p() + state[iEmt].p()).m2Calc();
6952 pRadBef = ( state[iRad].p() - q*state[iRad].p()/q2par * qpar )
6953 * sqrt( (lABC(q2,m2s,m2Bef) - 4.*m2Bef*pT2k)
6954 / (lABC(q2,sjk,m2r) - 4.*m2r*pT2k))
6955 + qpar * 0.5 * (q2 + m2Bef - m2s) / q2par;
6959 pRecBef = -q+pRadBef;
6963 radBef.m(sqrtpos(m2r));
6964 recBef.m(sqrtpos(m2s));
6969 Vec4 paj(state[iRad].p()-state[iEmt].p());
6970 Vec4 pk(state[iRecAft].p());
6974 double q2 = q.m2Calc();
6975 double saj = 2.*state[iRad].p()*state[iEmt].p();
6978 pRecBef = (pk - q*pk/q2*q)
6979 *sqrt(lABC(q2,m2Bef,m2s)/lABC(q2,saj,m2s))
6980 + 0.5*(q2+m2s-m2Bef)/q2*q;
6983 pRadBef = -q+pRecBef;
6986 int iOther = getInB(iSysSel);
6987 if (side == -1) iOther = getInA(iSysSel);
6988 Vec4 pOther(state[iOther].p());
6991 RotBstMatrix toABCM;
6992 if (side == 1) toABCM.toCMframe( pRadBef, pOther);
6993 else toABCM.toCMframe( pOther, pRadBef);
6996 pRadBef.rotbst(toABCM);
6997 pOther.rotbst(toABCM);
7001 RotBstMatrix restoreB;
7002 restoreB.bst( pOther, state[iOther].p());
7005 pRadBef.rotbst(restoreB);
7006 pOther.rotbst(restoreB);
7011 radBef.m(sqrtpos(m2r));
7012 recBef.m(sqrtpos(m2s));
7016 for (
int i = 0; i < partial.size(); ++i) {
7017 if ( !partial[i].isFinal() && partial[i].statusAbs() != 22 )
continue;
7018 partial[i].rotbst(toABCM);
7019 partial[i].rotbst(restoreB);
7037 int DireSpace::FindParticle(
const Particle& particle,
const Event& event,
7038 bool checkStatus ) {
7042 for (
int i =
int(event.size()) - 1; i > 0; --i )
7043 if ( event[i].
id() == particle.id()
7044 &&
event[i].colType() == particle.colType()
7045 &&
event[i].chargeType() == particle.chargeType()
7046 &&
event[i].col() == particle.col()
7047 &&
event[i].acol() == particle.acol()
7048 &&
event[i].charge() == particle.charge() ) {
7053 if ( checkStatus && event[index].status() != particle.status() )
7061 double DireSpace::pT2_II (
const Particle& rad,
const Particle& emt,
7062 const Particle& rec) {
7063 double sai = -2.*rad.p()*emt.p();
7064 double sbi = -2.*rec.p()*emt.p();
7065 double sab = 2.*rad.p()*rec.p();
7067 return sai*sbi / sab * ( sai+sbi+sab ) / sab;
7073 double DireSpace::pT2_IF (
const Particle& rad,
const Particle& emt,
7074 const Particle& rec) {
7075 double sai = -2.*rad.p()*emt.p();
7076 double sik = 2.*rec.p()*emt.p();
7077 double sak = -2.*rad.p()*rec.p();
7078 return sai*sik / (sai+sak) * (sai+sik+sak) / (sai+sak);
7083 double DireSpace::z_II (
const Particle& rad,
const Particle& emt,
7084 const Particle& rec) {
7085 double sbi = -2.*rec.p()*emt.p();
7086 double sab = 2.*rad.p()*rec.p();
7087 return 1. + sbi/sab;
7092 double DireSpace::z_IF (
const Particle& rad,
const Particle& emt,
7093 const Particle& rec) {
7094 double sai = -2.*rad.p()*emt.p();
7095 double sik = 2.*rec.p()*emt.p();
7096 double sak = -2.*rad.p()*rec.p();
7097 return 1. + sik / (sai+sak);
7102 double DireSpace::m2dip_II (
const Particle& rad,
const Particle& emt,
7103 const Particle& rec) {
7104 double sai = -2.*rad.p()*emt.p();
7105 double sbi = -2.*rec.p()*emt.p();
7106 double sab = 2.*rad.p()*rec.p();
7107 return (sab + sai + sbi);
7112 double DireSpace::m2dip_IF (
const Particle& rad,
const Particle& emt,
7113 const Particle& rec) {
7114 double sai = -2.*rad.p()*emt.p();
7115 double sik = 2.*rec.p()*emt.p();
7116 double sak = -2.*rad.p()*rec.p();
7117 return -1.*(sai+sik+sak);
7126 map<string, double> DireSpace::getStateVariables (
const Event& state,
7127 int rad,
int emt,
int rec,
string name) {
7128 map<string,double> ret;
7131 if (rad > 0 && emt > 0 && rec > 0) {
7132 double pT2 = pT2Space ( state[rad], state[emt], state[rec]);
7133 double z = zSpace ( state[rad], state[emt], state[rec]);
7134 ret.insert(make_pair(
"t",pT2));
7135 ret.insert(make_pair(
"tRS",pT2));
7136 ret.insert(make_pair(
"scaleAS",pT2));
7137 ret.insert(make_pair(
"scaleEM",pT2));
7138 ret.insert(make_pair(
"scalePDF",pT2));
7139 ret.insert(make_pair(
"z",z));
7144 ? (*splittingsPtr)[name]->radBefID(state[rad].
id(), state[emt].
id())
7146 pair<int,int> radBefCols
7148 ? (*splittingsPtr)[name]->radBefCols(state[rad].col(),
7149 state[rad].acol(), state[emt].col(), state[emt].acol())
7151 ret.insert(make_pair(
"radBefID", radBefID));
7152 ret.insert(make_pair(
"radBefCol", radBefCols.first));
7153 ret.insert(make_pair(
"radBefAcol", radBefCols.second));
7157 ? (*splittingsPtr)[name]->couplingType(state[rad].
id(), state[emt].
id())
7159 double couplingValue
7161 ? (*splittingsPtr)[name]->coupling(z,pT2)
7163 ret.insert(make_pair
7164 (
"scaleForCoupling " + std::to_string(couplingType), pT2));
7165 ret.insert(make_pair(
"couplingType",couplingType));
7166 ret.insert(make_pair(
"couplingValue",couplingValue));
7168 double m2dip = m2dipSpace ( state[rad], state[emt], state[rec]);
7169 ret.insert(make_pair(
"m2dip", m2dip));
7175 ret.insert(make_pair(
"t",0.));
7176 ret.insert(make_pair(
"tRS",0.));
7177 ret.insert(make_pair(
"scaleAS",0.));
7178 ret.insert(make_pair(
"scaleEM",0.));
7179 ret.insert(make_pair(
"z",0.));
7180 ret.insert(make_pair(
"radBefID", 0));
7181 ret.insert(make_pair(
"radBefCol", 0));
7182 ret.insert(make_pair(
"radBefAcol", 0));
7183 ret.insert(make_pair(
"scaleForCoupling "+std::to_string(-1),0.));
7184 ret.insert(make_pair(
"couplingType",-1));
7185 ret.insert(make_pair(
"couplingValue",-1.));
7190 vector<DireSpaceEnd> dipEnds;
7192 int colTag = state[in1].col();
7193 if (colTag > 0) getQCDdip( in1, colTag, 1, state, dipEnds);
7195 int acolTag = state[in1].acol();
7196 if (acolTag > 0) getQCDdip( in1, acolTag, -1, state, dipEnds);
7198 colTag = state[in2].col();
7199 if (colTag > 0) getQCDdip( in2, colTag, 1, state, dipEnds);
7201 acolTag = state[in2].acol();
7202 if (acolTag > 0) getQCDdip( in2, acolTag, -1, state, dipEnds);
7205 getGenDip( -1, in1, state,
false, dipEnds);
7206 getGenDip( -1, in2, state,
false, dipEnds);
7209 double x1 = state[3].pPos() / state[0].m();
7210 double x2 = state[4].pNeg() / state[0].m();
7214 for (
int iDip = 0; iDip < int(dipEnds.size()); ++iDip) {
7215 double m2 = abs(2.*state[dipEnds[iDip].iRadiator].p()
7216 *state[dipEnds[iDip].iRecoiler].p());
7217 if ( dipEnds[iDip].iRadiator == in1) m2 /= x1;
7218 if ( dipEnds[iDip].iRecoiler == in1) m2 /= x1;
7219 if ( dipEnds[iDip].iRadiator == in2) m2 /= x2;
7220 if ( dipEnds[iDip].iRecoiler == in2) m2 /= x2;
7222 oss <<
"scalePDF-" << dipEnds[iDip].iRadiator
7223 <<
"-" << dipEnds[iDip].iRecoiler;
7224 ret.insert(make_pair(oss.str(),m2));
7235 double DireSpace::getSplittingProb(
const Event& state,
int iRad,
7236 int iEmt,
int iRecAft,
string name) {
7239 int order = atoi( (
char*)name.substr( name.find(
"-",0)+1,
7240 name.size() ).c_str() );
7241 name=name.substr( 0, name.size()-2);
7245 if ( splits[name]->splitInfo.extras.find(
"unitKernel")
7246 != splits[name]->splitInfo.extras.end() )
return 1.;
7248 double z = zSpace(state[iRad], state[iEmt], state[iRecAft]);
7249 double pT2 = pT2Space(state[iRad], state[iEmt], state[iRecAft]);
7250 double m2D = (state[iRecAft].isFinal())
7251 ? abs( 2.*state[iEmt].p()*state[iRad].p()
7252 -2.*state[iEmt].p()*state[iRecAft].p()
7253 +2.*state[iRad].p()*state[iRecAft].p())
7254 : abs(-2.*state[iEmt].p()*state[iRad].p()
7255 -2.*state[iEmt].p()*state[iRecAft].p()
7256 +2.*state[iRad].p()*state[iRecAft].p());
7257 int idRadBef = splits[name]->radBefID(state[iRad].
id(), state[iEmt].
id());
7258 double m2Bef = ( abs(idRadBef) < 6)
7259 ? getMass(idRadBef,2)
7260 : (idRadBef == state[iRad].id())
7261 ? getMass(idRadBef,3,state[iRad].mCalc())
7262 : getMass(idRadBef,2);
7263 double m2r = state[iRad].p().m2Calc();
7264 double m2e = state[iEmt].p().m2Calc();
7265 double m2s = state[iRecAft].p().m2Calc();
7266 int type = (state[iRecAft].isFinal()) ? 1 : -1;
7269 if ( pT2cut(state[iEmt].
id()) > pT2)
return 0.;
7270 if ( !splits[name]->aboveCutoff( pT2, state[iRad], state[iRecAft], 0,
7271 partonSystemsPtr))
return 0.;
7274 if (type == 1 && (m2Bef > TINYMASS || m2r > TINYMASS || m2e > TINYMASS
7275 || m2s > TINYMASS)) type = 2;
7276 if (type ==-1 && (m2Bef > TINYMASS || m2r > TINYMASS || m2e > TINYMASS
7277 || m2s > TINYMASS)) type =-2;
7280 double m2dipCorr = m2D;
7281 double kappa2 = pT2 / m2dipCorr;
7282 double xCS = (state[iRecAft].isFinal()
7283 ? z : (z*(1-z) - kappa2) / (1 -z));
7284 double xMot = 2.*state[iRad].e()/state[0].m();
7285 double xDau = xCS * xMot;
7288 double xMin = (particleDataPtr->colType(idRadBef) != 0) ? xDau : 0.;
7290 if ( !inAllowedPhasespace( 1, z, pT2, m2dipCorr, xMin, type, m2Bef, m2r, m2s,
7291 m2e) ) {
return 0.0;}
7297 int massSign = type > 0 ? 1 : -1;
7298 pair<Vec4, Vec4> pTdirection = getTwoPerpendicular(state[iRad].p(),
7299 massSign*state[iRecAft].p()+state[iEmt].p());
7300 double px= -pTdirection.first*state[iEmt].p();
7301 double py= -pTdirection.second*state[iEmt].p();
7302 double kT2 = pow2(px)+pow2(py);
7303 double phi1 = atan2(px/sqrt(kT2), py/sqrt(kT2));
7304 if (phi1 < 0.) phi1 = 2.*M_PI+phi1;
7307 pair <Event, pair<int,int> > born
7308 (clustered_internal( state, iRad, iEmt, iRecAft, name ));
7309 int nEmissions = splits[name]->nEmissions();
7310 double m2dipBef = abs(2.*born.first[born.second.first].p()
7311 *born.first[born.second.second].p());
7312 splits[name]->splitInfo.clear();
7313 splits[name]->splitInfo.storeInfo(name, type, 0, 0, 0,
7314 born.second.first, born.second.second,
7315 born.first, state[iEmt].id(), state[iRad].id(),
7316 nEmissions, m2dipBef, pT2, pT2, z, phi1, m2Bef, m2s,
7317 (nEmissions == 1 ? m2r : 0.0),(nEmissions == 1 ? m2e : 0.0),
7318 0.0, 0.0, 0.0, 0.0, xDau, xMot);
7321 unordered_map < string, double > kernels;
7323 if (splits[name]->calc(born.first, order) ) kernels =
7324 splits[name]->getKernelVals();
7325 if ( kernels.find(
"base") != kernels.end() ) p += kernels[
"base"];
7327 splits[name]->splitInfo.clear();
7335 bool hasME = pT2 > pT2minMECs && doMEcorrections && weights->hasME(state);
7336 if (hasME) p = abs(p);
7338 if (!dryrun && splits[name]->hasMECAft(state, pT2)) p *= KERNEL_HEADROOM;
7342 for (
int i=0; i < state.size(); ++i)
if (state[i].isFinal()) nFinal++;
7343 if (!dryrun) mecover = splits[name]->overhead(m2dipBef*xDau,
7344 state[iRad].
id(), nFinal-1);
7356 bool DireSpace::allowedSplitting(
const Event& state,
int iRad,
int iEmt) {
7358 bool isAP = state[iEmt].id() < 0;
7359 int idRad = state[iRad].id();
7360 int idEmt = state[iEmt].id();
7362 int colRad = state[iRad].col();
7363 int acolRad = state[iRad].acol();
7364 int colEmt = state[iEmt].col();
7365 int acolEmt = state[iEmt].acol();
7367 int colShared = (colRad > 0 && colRad == colEmt ) ? colRad
7368 : (acolRad > 0 && acolRad == acolEmt) ? acolRad : 0;
7371 if ( state[iRad].isFinal() )
return false;
7374 if (idEmt == 21 && colShared > 0)
7378 if ( abs(idRad) < 10 && idRad == idEmt && colShared == 0)
7382 if ( idRad == 21 && abs(idEmt) < 10
7383 && ( (isAP && acolEmt == acolRad) || (!isAP && colEmt == colRad) ) )
7387 if ( idEmt == 22 && abs(idRad) < 10)
7391 if ( idEmt == 22 && (abs(idRad) == 11 || abs(idRad) == 13
7392 || abs(idRad) == 15))
7396 if ( abs(idEmt) < 10 && idRad == idEmt && colShared > 0)
7400 if ( (abs(idEmt) == 11 || abs(idEmt) == 13 || abs(idEmt) == 15)
7405 if ( idRad == 22 && abs(idEmt) < 10 && idEmt == idRad && colShared > 0)
7409 if (idRad == 22 && (abs(idEmt) == 11 || abs(idEmt) == 13 || abs(idEmt) == 15)
7414 if ( idEmt == 23 && abs(idRad) < 10)
7418 if ( idEmt == 23 && (abs(idRad) == 11 || abs(idRad) == 13
7419 || abs(idRad) == 15))
7428 vector<int> DireSpace::getRecoilers(
const Event& state,
int iRad,
int iEmt,
7431 return splits[name]->recPositions(state, iRad, iEmt);
7436 Event DireSpace::makeHardEvent(
int iSys,
const Event& state,
bool isProcess) {
7438 bool hasSystems = !isProcess && partonSystemsPtr->sizeSys() > 0;
7439 int sizeSys = (hasSystems) ? partonSystemsPtr->sizeSys() : 1;
7443 event.init(
"(hard process-modified)", particleDataPtr );
7447 for (
int i = state.size()-1; i > 0; --i)
7448 if ( state[i].mother1() == 1 && state[i].mother2() == 0
7449 && (!hasSystems || partonSystemsPtr->getSystemOf(i,
true) == iSys))
7451 if (in1 == 0) in1 = getInA(iSys);
7453 for (
int i = state.size()-1; i > 0; --i)
7454 if ( state[i].mother1() == 2 && state[i].mother2() == 0
7455 && (!hasSystems || partonSystemsPtr->getSystemOf(i,
true) == iSys))
7457 if (in2 == 0) in2 = getInB(iSys);
7461 bool resonantIncoming =
false;
7462 if ( in1 == 0 && in2 == 0 ) {
7463 int iParentInOther = 0;
7464 int nSys = partonSystemsPtr->sizeAll(iSys);
7465 for (
int iInSys = 0; iInSys < nSys; ++iInSys){
7466 int iiNow = partonSystemsPtr->getAll(iSys,iInSys);
7467 bool hasOtherParent =
false;
7468 for (
int iOtherSys = 0; iOtherSys < sizeSys; ++iOtherSys){
7469 if (iOtherSys == iSys)
continue;
7470 int nOtherSys = partonSystemsPtr->sizeAll(iOtherSys);
7471 for (
int iInOtherSys = 0; iInOtherSys < nOtherSys; ++iInOtherSys){
7472 int iOtherNow = partonSystemsPtr->getAll(iOtherSys,iInOtherSys);
7473 if (state[iiNow].isAncestor(iOtherNow)) {
7474 iParentInOther = iOtherNow;
7475 hasOtherParent =
true;
7479 if (hasOtherParent)
break;
7481 if (hasOtherParent)
break;
7483 in1 = iParentInOther;
7484 if (iParentInOther) resonantIncoming =
true;
7487 event.append(state[0]);
7488 event.append(state[1]);
7489 event[1].daughters(3,0);
7490 event.append(state[2]);
7491 event[2].daughters(4,0);
7494 event.append(state[in1]);
7495 event[3].mothers(1,0);
7496 if (resonantIncoming)
event[3].status(-22);
7497 else event[3].status(-21);
7500 event.append(state[in2]);
7501 event[4].mothers(2,0);
7502 event[4].status(-21);
7504 for (
int i = 0; i < state.size(); ++i) {
7509 bool isFin = state[i].isFinal();
7510 bool isInSys = (partonSystemsPtr->getSystemOf(i) == iSys);
7512 bool isParentOfOther =
false;
7513 if (!isFin && isInSys) {
7514 for (
int iOtherSys = 0; iOtherSys < sizeSys; ++iOtherSys){
7515 if (iOtherSys == iSys)
continue;
7516 double nSys = partonSystemsPtr->sizeAll(iOtherSys);
7517 for (
int iInSys = 0; iInSys < nSys; ++iInSys){
7518 int iiNow = partonSystemsPtr->getAll(iOtherSys,iInSys);
7519 if (state[iiNow].isAncestor(i)) {isParentOfOther=
true;
break;}
7524 if ( (isFin || isParentOfOther) && (!hasSystems || isInSys) ) {
7525 int iN =
event.append(state[i]);
7526 event[iN].daughters(0,0);
7527 event[iN].mothers(3,4);
7528 int status = (state[i].statusAbs() == 22) ? state[i].statusAbs() : 23;
7529 if ( particleDataPtr->isResonance(state[i].id()) ) status = 22;
7530 event[iN].status(status);
7535 event[3].daughters(5,event.size()-1);
7536 event[4].daughters(5,event.size()-1);
7546 bool DireSpace::validMomentum(
const Vec4& p,
int id,
int status) {
7549 if (isnan(p) || isinf(p))
return false;
7552 double mNow = (status < 0) ? 0.
7553 : ((abs(
id) < 6) ? getMass(
id,2) : getMass(id,1));
7555 if (status < 0 && useMassiveBeams
7556 && (abs(
id) == 11 || abs(
id) == 13 || abs(
id) > 900000))
7557 mNow = getMass(
id,1);
7563 if ( abs(
id) == 6 || abs(
id) > 22) mNow = p.mCalc();
7564 double errMass = abs(p.mCalc() - mNow) / max( 1.0, p.e());
7565 if ( errMass > mTolErr )
return false;
7568 if ( p.e() < 0. )
return false;
7579 bool DireSpace::validEvent(
const Event& state,
bool isProcess ) {
7581 bool validColour =
true;
7582 bool validCharge =
true;
7583 bool validMomenta =
true;
7585 bool hasSystems = !isProcess && partonSystemsPtr->sizeSys() > 0;
7586 int sizeSys = (hasSystems) ? partonSystemsPtr->sizeSys() : 1;
7590 for (
int i = 0; i < state.size(); ++i)
7591 if (isnan(state[i].p()) || isinf(state[i].p()))
return false;
7593 for (
int iSys = 0; iSys < sizeSys; ++iSys) {
7596 if (!validColour || !validCharge )
break;
7599 e.init(
"(hard process-modified)", particleDataPtr );
7601 e = makeHardEvent(iSys, state, isProcess);
7604 for (
int i = 0; i < e.size(); ++i)
7606 if ( e[i].isFinal() && e[i].colType() == 1
7608 && ( FindCol(e[i].col(),vector<int>(1,i),e,1) == 0
7610 && FindCol(e[i].col(),vector<int>(1,i),e,2) == 0 )) {
7611 validColour =
false;
7614 }
else if ( e[i].isFinal() && e[i].colType() == -1
7616 && ( FindCol(e[i].acol(),vector<int>(1,i),e,2) == 0
7618 && FindCol(e[i].acol(),vector<int>(1,i),e,1) == 0 )) {
7619 validColour =
false;
7622 }
else if ( e[i].isFinal() && e[i].colType() == 2
7624 && ( FindCol(e[i].col(),vector<int>(1,i),e,1) == 0
7626 && FindCol(e[i].col(),vector<int>(1,i),e,2) == 0 )
7628 && ( FindCol(e[i].acol(),vector<int>(1,i),e,2) == 0
7630 && FindCol(e[i].acol(),vector<int>(1,i),e,1) == 0 )) {
7631 validColour =
false;
7636 double initCharge = e[3].charge() + e[4].charge();
7637 double finalCharge = 0.0;
7638 for(
int i = 0; i < e.size(); ++i)
7639 if (e[i].isFinal()) finalCharge += e[i].charge();
7640 if (abs(initCharge-finalCharge) > 1e-12) validCharge =
false;
7643 for (
int i = 0; i < e.size(); ++i) {
7644 if (i==3 || i==4 || e[i].isFinal()) {
7645 validMomenta = validMomenta
7646 && validMomentum(e[i].p(), e[i].
id(), (e[i].isFinal() ? 1 : -1));
7651 Vec4 pSum(0.,0.,0.,0.);
7652 for (
int i = 0; i < e.size(); ++i) {
7653 if ( e[i].status() == -21
7654 || e[i].status() == -22 ) pSum -= e[i].p();
7655 if ( e[i].isFinal() ) pSum += e[i].p();
7657 if ( abs(pSum.px()) > mTolErr || abs(pSum.py()) > mTolErr)
7658 validMomenta =
false;
7659 if ( e[3].status() == -21
7660 && (abs(e[3].px()) > mTolErr || abs(e[3].py()) > mTolErr))
7661 validMomenta =
false;
7662 if ( e[4].status() == -21
7663 && (abs(e[4].px()) > mTolErr || abs(e[4].py()) > mTolErr))
7664 validMomenta =
false;
7668 return (validColour && validCharge && validMomenta);
7674 bool DireSpace::validMotherDaughter(
const Event& event ) {
7678 vector< pair<int,int> > noMotDau;
7681 bool hasBeams =
false;
7682 for (
int i = 0; i <
event.size(); ++i) {
7683 int status =
event[i].status();
7684 if (abs(status) == 12) hasBeams =
true;
7687 vector<int> mList =
event[i].motherList();
7688 vector<int> dList =
event[i].daughterList();
7689 if (mList.size() == 0 && abs(status) != 11 && abs(status) != 12)
7691 if (dList.size() == 0 && status < 0 && status != -11)
7695 for (
int j = 0; j < int(mList.size()); ++j) {
7696 if ( event[mList[j]].daughter1() <= i
7697 &&
event[mList[j]].daughter2() >= i )
continue;
7698 vector<int> dmList =
event[mList[j]].daughterList();
7699 bool foundMatch =
false;
7700 for (
int k = 0; k < int(dmList.size()); ++k)
7701 if (dmList[k] == i) {
7705 if (!hasBeams && mList.size() == 1 && mList[0] == 0) foundMatch =
true;
7707 bool oldPair =
false;
7708 for (
int k = 0; k < int(noMotDau.size()); ++k)
7709 if (noMotDau[k].first == mList[j] && noMotDau[k].second == i) {
7713 if (!oldPair) noMotDau.push_back( make_pair( mList[j], i) );
7718 for (
int j = 0; j < int(dList.size()); ++j) {
7719 if ( event[dList[j]].statusAbs() > 80
7720 &&
event[dList[j]].statusAbs() < 90
7721 &&
event[dList[j]].mother1() <= i
7722 &&
event[dList[j]].mother2() >= i)
continue;
7723 vector<int> mdList =
event[dList[j]].motherList();
7724 bool foundMatch =
false;
7725 for (
int k = 0; k < int(mdList.size()); ++k)
7726 if (mdList[k] == i) {
7731 bool oldPair =
false;
7732 for (
int k = 0; k < int(noMotDau.size()); ++k)
7733 if (noMotDau[k].first == i && noMotDau[k].second == dList[j]) {
7737 if (!oldPair) noMotDau.push_back( make_pair( i, dList[j]) );
7744 if (noMot.size() > 0 || noDau.size() > 0 || noMotDau.size() > 0)
7756 int DireSpace::FindCol(
int col, vector<int> iExc,
const Event& event,
7757 int type,
int iSys) {
7761 int inA = 0, inB = 0;
7762 for (
int i=event.size()-1; i > 0; --i) {
7763 if ( event[i].mother1() == 1 &&
event[i].status() != -31
7764 &&
event[i].status() != -34) {
if (inA == 0) inA = i; }
7765 if ( event[i].mother1() == 2 &&
event[i].status() != -31
7766 &&
event[i].status() != -34) {
if (inB == 0) inB = i; }
7768 if (iSys >= 0) {inA = getInA(iSys); inB = getInB(iSys);}
7771 for(
int n = 0; n <
event.size(); ++n) {
7773 if ( find(iExc.begin(), iExc.end(), n) != iExc.end() )
continue;
7774 if ( event[n].colType() != 0 &&
event[n].status() > 0 ) {
7775 if ( event[n].acol() == col ) {
7779 if ( event[n].col() == col ) {
7786 for(
int n = event.size()-1; n > 0; --n) {
7788 if ( find(iExc.begin(), iExc.end(), n) != iExc.end() )
continue;
7789 if ( index == 0 && event[n].colType() != 0
7790 && ( n == inA || n == inB) ) {
7791 if ( event[n].acol() == col ) {
7795 if ( event[n].col() == col ) {
7803 if ( type == 1 && index < 0)
return abs(index);
7804 if ( type == 2 && index > 0)
return abs(index);
7813 void DireSpace::list()
const {
7816 cout <<
"\n -------- DIRE DireSpace Dipole Listing ------------------"
7817 <<
"--------------------------------------------------------------"
7819 <<
" i sys side rad rec pTmax col "
7820 <<
" m2Dip siblings allowedIDs\n"
7821 << fixed << setprecision(3);
7824 for (
int i = 0; i < int(dipEnd.size()); ++i) {
7825 cout << scientific << setprecision(4)
7826 << setw(4) << i <<
" | "
7827 << setw(4) << dipEnd[i].system <<
" | "
7828 << setw(4) << dipEnd[i].side <<
" | "
7829 << setw(4) << dipEnd[i].iRadiator <<
" | "
7830 << setw(4) << dipEnd[i].iRecoiler <<
" | "
7831 << setw(11) << dipEnd[i].pTmax <<
" | "
7832 << setw(3) << dipEnd[i].colType <<
" | "
7833 << setw(12) << dipEnd[i].m2Dip <<
" | ";
7835 os << dipEnd[i].iSiblings.listPos();
7836 cout << setw(15) << os.str() <<
" | ";
7838 for (
int j = 0; j < int(dipEnd[i].allowedEmissions.size()); ++j)
7839 os << setw(4) << dipEnd[i].allowedEmissions[j];
7840 cout << setw(15) << os.str() << endl;
7844 cout <<
"\n -------- End DIRE DireSpace Dipole Listing --------------"
7845 <<
"--------------------------------------------------------------"
7846 <<
"----------" << endl;
7849 for ( unordered_map<string,DireSplitting*>::const_iterator it
7850 = splits.begin(); it != splits.end(); ++it ) {
7851 multimap<double,OverheadInfo> bla = it->second->overhead_map;
7852 cout << it->first << endl;
7853 for ( multimap<double, OverheadInfo >::const_iterator itb = bla.begin();
7854 itb != bla.end(); ++itb )
7855 cout <<
" pT2=" << itb->first <<
" " << itb->second.list() << endl;
7866 double DireSpace::alphasNow(
double pT2,
double renormMultFacNow,
int iSys ) {
7869 BeamParticle* beam = (particleDataPtr->isHadron(beamAPtr->id()))
7871 : (particleDataPtr->isHadron(beamBPtr->id()) ? beamBPtr :
7873 if (usePDFalphas && beam == NULL) beam = beamAPtr;
7874 double scale = pT2*renormMultFacNow;
7875 scale = max(scale, pT2min);
7878 double asPT2pi = (usePDFalphas && beam != NULL)
7879 ? beam->alphaS(scale) / (2.*M_PI)
7880 : alphaS.alphaS(scale) / (2.*M_PI);
7883 int order = kernelOrder-1;
7885 bool hasInA = (getInA(iSys) != 0);
7886 bool hasInB = (getInB(iSys) != 0);
7887 if (iSys != 0 && hasInA && hasInB) order = kernelOrderMPI-1;
7891 double m2cNow(m2cPhys), m2bNow(m2bPhys);
7892 if ( !( (scale > m2cNow && pT2 < m2cNow)
7893 || (scale < m2cNow && pT2 > m2cNow) ) ) m2cNow = -1.;
7894 if ( !( (scale > m2bNow && pT2 < m2bNow)
7895 || (scale < m2bNow && pT2 > m2bNow) ) ) m2bNow = -1.;
7896 vector<double> scales;
7897 scales.push_back(scale);
7898 scales.push_back(pT2);
7899 if (m2cNow > 0.) scales.push_back(m2cNow);
7900 if (m2bNow > 0.) scales.push_back(m2bNow);
7901 sort( scales.begin(), scales.end());
7902 if (scale > pT2) reverse(scales.begin(), scales.end());
7904 double asPT2piCorr = asPT2pi;
7905 for (
int i = 1; i< int(scales.size()); ++i) {
7906 double NF = getNF( 0.5*(scales[i]+scales[i-1]) );
7907 double L = log( scales[i]/scales[i-1] );
7909 if (order > 0) subt += asPT2piCorr * beta0(NF) * L;
7910 if (order > 2) subt += pow2( asPT2piCorr ) * ( beta1(NF)*L
7911 - pow2(beta0(NF)*L) );
7912 if (order > 4) subt += pow( asPT2piCorr, 3) * ( beta2(NF)*L
7913 - 2.5 * beta0(NF)*beta1(NF)*L*L
7914 + pow( beta0(NF)*L, 3) );
7915 asPT2piCorr *= 1.0 - subt;
7927 double DireSpace::getNF(
double pT2) {
7929 BeamParticle* beam = (particleDataPtr->isHadron(beamAPtr->id()))
7931 : (particleDataPtr->isHadron(beamBPtr->id()) ? beamBPtr :
7934 if ( !usePDFalphas || beam == NULL ) {
7935 if ( pT2 > pow2( max(0., particleDataPtr->m0(5) ) )
7936 && pT2 < pow2( particleDataPtr->m0(6)) ) NF = 5.;
7937 else if ( pT2 > pow2( max( 0., particleDataPtr->m0(4)) ) ) NF = 4.;
7938 else if ( pT2 > pow2( max( 0., particleDataPtr->m0(3)) ) ) NF = 3.;
7940 if ( pT2 > pow2( max(0., beam->mQuarkPDF(5) ) )
7941 && pT2 < pow2( particleDataPtr->m0(6)) ) NF = 5.;
7942 else if ( pT2 > pow2( max( 0., beam->mQuarkPDF(4)) ) ) NF = 4.;
7943 else if ( pT2 > pow2( max( 0., beam->mQuarkPDF(3)) ) ) NF = 3.;