1 #include "Pythia8/Pythia.h" 21 const Pythia8::Event &
e,
int RadAfterBranch,
int EmtAfterBranch,
int RecAfterBranch,
bool FSR) {
26 int radID =
e[RadAfterBranch].id();
29 double sign = (FSR) ? 1. : -1.;
31 double Qsq =
sign * Q.m2Calc();
34 double m2Rad = (
abs(radID) >= 4 &&
abs(radID) < 7) ?
Pythia8::pow2(particleDataPtr->m0(radID)) : 0.;
41 double m2Dip = sum.m2Calc();
42 double x1 = 2. * (sum * radVec) / m2Dip;
43 double x3 = 2. * (sum * emtVec) / m2Dip;
51 z = qBR.m2Calc() / qAR.m2Calc();
56 pTnow *= (Qsq -
sign * m2Rad);
60 std::cout <<
"Warning: pTpythia was negative" << std::endl;
65 std::cout <<
"pTpythia: rad = " << RadAfterBranch <<
", emt = " << EmtAfterBranch <<
", rec = " << RecAfterBranch
66 <<
", pTnow = " <<
sqrt(pTnow) << std::endl;
81 int iInA = partonSystemsPtr->getInA(0);
82 int iInB = partonSystemsPtr->getInB(0);
83 double betaZ = -(
e[iInA].pz() +
e[iInB].pz()) / (
e[iInA].
e() +
e[iInB].e());
85 iVecBst.bst(0., 0., betaZ);
86 jVecBst.bst(0., 0., betaZ);
87 pTnow =
sqrt((iVecBst + jVecBst).m2Calc() * iVecBst.e() * jVecBst.e() /
Pythia8::pow2(iVecBst.e() + jVecBst.e()));
96 std::cout <<
"Warning: pTpowheg was negative" << std::endl;
101 std::cout <<
"pTpowheg: i = " <<
i <<
", j = " <<
j <<
", pTnow = " << pTnow << std::endl;
114 double pTemt = -1., pTnow;
115 int xSR1 = (xSRin == -1) ? 0 : xSRin;
116 int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
117 for (
int xSR = xSR1; xSR < xSR2; xSR++) {
119 bool FSR = (xSR == 0) ?
false :
true;
123 if ((pTdefMode == 0 || pTdefMode == 1) &&
i > 0 &&
j > 0) {
124 pTemt = pTpowheg(
e,
i,
j, (pTdefMode == 0) ?
false : FSR);
127 }
else if (!FSR && pTdefMode == 2 &&
i > 0 &&
j > 0 && r > 0) {
128 pTemt = pTpythia(
e,
i,
j, r, FSR);
131 }
else if (FSR && pTdefMode == 2 &&
j > 0 &&
k > 0 && r > 0) {
132 pTemt = pTpythia(
e,
k,
j, r, FSR);
140 int iInA = partonSystemsPtr->getInA(0);
141 int iInB = partonSystemsPtr->getInB(0);
142 while (
e[iInA].mother1() != 1) {
143 iInA =
e[iInA].mother1();
145 while (
e[iInB].mother1() != 2) {
146 iInB =
e[iInB].mother1();
150 int jNow = (
j > 0) ?
j : 0;
151 int jMax = (
j > 0) ?
j + 1 :
e.size();
152 for (; jNow < jMax; jNow++) {
154 if (!
e[jNow].isFinal())
157 if (QEDvetoMode == 0 &&
e[jNow].colType() == 0)
161 if (pTdefMode == 0 || pTdefMode == 1) {
164 pTnow = pTpowheg(
e, iInA, jNow, (pTdefMode == 0) ?
false : FSR);
166 pTemt = (pTemt < 0) ? pTnow :
std::min(pTemt, pTnow);
172 int outSize = partonSystemsPtr->sizeOut(0);
173 for (
int iMem = 0; iMem < outSize; iMem++) {
174 int iNow = partonSystemsPtr->getOut(0, iMem);
180 if (QEDvetoMode == 0 &&
e[iNow].colType() == 0)
182 if (jNow ==
e[iNow].daughter1() && jNow ==
e[iNow].daughter2())
185 pTnow = pTpowheg(
e, iNow, jNow, (pTdefMode == 0) ?
false : FSR);
187 pTemt = (pTemt < 0) ? pTnow :
std::min(pTemt, pTnow);
193 }
else if (pTdefMode == 2) {
196 pTnow = pTpythia(
e, iInA, jNow, iInB, FSR);
198 pTemt = (pTemt < 0) ? pTnow :
std::min(pTemt, pTnow);
199 pTnow = pTpythia(
e, iInB, jNow, iInA, FSR);
201 pTemt = (pTemt < 0) ? pTnow :
std::min(pTemt, pTnow);
206 for (
int kNow = 0; kNow <
e.size(); kNow++) {
207 if (kNow == jNow || !
e[kNow].isFinal())
209 if (QEDvetoMode == 0 &&
e[kNow].colType() == 0)
214 pTnow = pTpythia(
e, kNow, jNow, iInA, FSR);
216 pTemt = (pTemt < 0) ? pTnow :
std::min(pTemt, pTnow);
217 pTnow = pTpythia(
e, kNow, jNow, iInB, FSR);
219 pTemt = (pTemt < 0) ? pTnow :
std::min(pTemt, pTnow);
222 for (
int rNow = 0; rNow <
e.size(); rNow++) {
223 if (rNow == kNow || rNow == jNow || !
e[rNow].isFinal())
225 if (QEDvetoMode == 0 &&
e[rNow].colType() == 0)
227 pTnow = pTpythia(
e, kNow, jNow, rNow, FSR);
229 pTemt = (pTemt < 0) ? pTnow :
std::min(pTemt, pTnow);
240 std::cout <<
"pTcalc: i = " <<
i <<
", j = " <<
j <<
", k = " <<
k <<
", r = " << r <<
", xSR = " << xSRin
241 <<
", pTemt = " << pTemt << std::endl;
254 if (nFinalMode == 3 && pThardMode != 0)
256 "When nFinalMode is set to 3, ptHardMode should be set to 0, since the emission variables in doVetoMPIStep are " 257 "not set correctly case when there are three possible particle Born particle counts."));
267 int count = 0, inonfinal = 0;
268 double pT1 = 0., pTsum = 0.;
272 for (
int i =
e.size() - 1;
i > 0;
i--) {
274 if (
e[
i].isFinal()) {
283 if (nFinalMode == 2) {
284 if ((
abs(
e[
i].
id()) < 6 ||
e[
i].
id() == 21) &&
285 (
abs(
e[
e[
i].mother1()].
id()) < 6 ||
e[
e[
i].mother1()].
id() == 21)) {
295 if (nFinal < 0 || nFinalMode == 1) {
296 int first = -1, myid;
298 for (
int ip = 2; ip <
e.size(); ip++) {
300 if (
abs(myid) < 6 ||
abs(myid) == 21)
306 fatalEmissionVeto(
std::string(
"signal particles not found"));
307 for (
int ip =
first; ip <
e.size(); ip++) {
309 if (
abs(myid) < 6 ||
abs(myid) == 21)
313 nFinal =
last - inonfinal;
317 if (nFinalMode != 2) {
321 if (nFinalMode == 3) {
323 fatalEmissionVeto(
std::string(
"Wrong number of final state particles in event"));
326 fatalEmissionVeto(
std::string(
"Wrong number of final state particles in event"));
329 if (
count == nFinal + 1)
336 if (!isEmt || pThardMode == 0) {
337 pThard = infoPtr->QRen();
341 }
else if (pThardMode == 1) {
342 pThard = pTcalc(
e, -1, iEmt, -1, -1, -1);
347 }
else if (pThardMode == 2) {
348 pThard = pTcalc(
e, -1, -1, -1, -1, -1);
353 pTMPI = (isEmt) ? pTsum / 2. :
pT1;
357 std::cout <<
"doVetoMPIStep: QFac = " << infoPtr->QFac() <<
", QRen = " << infoPtr->QRen()
358 <<
", pThard = " << pThard << std::endl
379 if (vetoOn && nAcceptSeq >= vetoCount)
383 int iRadAft = -1, iEmt = -1, iRecAft = -1;
384 for (
int i =
e.size() - 1;
i > 0;
i--) {
385 if (iRadAft == -1 &&
e[
i].
status() == -41)
387 else if (iEmt == -1 &&
e[
i].
status() == 43)
389 else if (iRecAft == -1 &&
e[
i].
status() == -42)
391 if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
394 if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
396 fatalEmissionVeto(
std::string(
"Couldn't find Pythia ISR emission"));
402 int xSR = (pTemtMode == 0) ? 0 : -1;
403 int i = (pTemtMode == 0) ? iRadAft : -1;
404 int j = (pTemtMode != 2) ? iEmt : -1;
406 int r = (pTemtMode == 0) ? iRecAft : -1;
407 double pTemt = pTcalc(
e,
i,
j,
k, r, xSR);
410 std::cout <<
"doVetoISREmission: pTemt = " << pTemt << std::endl << std::endl;
415 bool vetoParton = (!isEmt &&
e[iEmt].colType() == 0 && QEDvetoMode == 2) ?
false :
true;
418 if (pTemt > pThard) {
421 nAcceptSeq = vetoCount - 1;
445 if (inResonance && settingsPtr->flag(
"POWHEG:bb4l:FSREmission:veto") == 1)
449 if (vetoOn && nAcceptSeq >= vetoCount)
453 int iRecAft =
e.size() - 1;
454 int iEmt =
e.size() - 2;
455 int iRadAft =
e.size() - 3;
456 int iRadBef =
e[iEmt].mother1();
460 fatalEmissionVeto(
std::string(
"Couldn't find Pythia FSR emission"));
467 int xSR = (pTemtMode == 0) ? 1 : -1;
468 int i = (pTemtMode == 0) ? iRadBef : -1;
469 int k = (pTemtMode == 0) ? iRadAft : -1;
470 int r = (pTemtMode == 0) ? iRecAft : -1;
474 if (pTemtMode == 0 || pTemtMode == 1) {
481 if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5))
484 for (
int jLoop = 0; jLoop < 2; jLoop++) {
486 pTemt = pTcalc(
e,
i,
j,
k, r, xSR);
491 if (emittedMode != 3)
500 }
else if (pTemtMode == 2) {
501 pTemt = pTcalc(
e,
i, -1,
k, r, xSR);
505 std::cout <<
"doVetoFSREmission: pTemt = " << pTemt << std::endl << std::endl;
510 bool vetoParton = (!isEmt &&
e[iEmt].colType() == 0 && QEDvetoMode == 2) ?
false :
true;
513 if (pTemt > pThard) {
516 nAcceptSeq = vetoCount - 1;
536 if (
e[
e.size() - 1].pT() > pTMPI) {
538 std::cout <<
"doVetoMPIEmission: pTnow = " <<
e[
e.size() - 1].pT() <<
", pTMPI = " << pTMPI << std::endl
bool doVetoFSREmission(int, const Pythia8::Event &e, int iSys, bool) override
void swap(Association< C > &lhs, Association< C > &rhs)
bool doVetoISREmission(int, const Pythia8::Event &e, int iSys) override
Abs< T >::type abs(const T &t)
double pTpowheg(const Pythia8::Event &e, int i, int j, bool FSR)
bool doVetoMPIStep(int nMPI, const Pythia8::Event &e) override
double pTcalc(const Pythia8::Event &e, int i, int j, int k, int r, int xSRin)
double pTpythia(const Pythia8::Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR)
constexpr int pow2(int x)
void fatalEmissionVeto(std::string message)
bool doVetoMPIEmission(int, const Pythia8::Event &e) override