14 #ifndef Pythia8_PowhegHooks_H 15 #define Pythia8_PowhegHooks_H 18 #include "Pythia8/Pythia.h" 19 #include "Pythia8/Plugins.h" 49 showerModel =
settingsPtr->mode(
"PartonShowers:model");
61 inline double pT(
const Event& e,
int RadAfterBranch,
62 int EmtAfterBranch,
int RecAfterBranch,
bool FSR) {
65 return pTvincia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch);
67 return pTpythia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch, FSR);
75 int EmtAfterBranch,
int RecAfterBranch,
bool FSR) {
78 Vec4 radVec = e[RadAfterBranch].p();
79 Vec4 emtVec = e[EmtAfterBranch].p();
80 Vec4 recVec = e[RecAfterBranch].p();
81 int radID = e[RadAfterBranch].id();
84 double sign = (FSR) ? 1. : -1.;
85 Vec4 Q(radVec + sign * emtVec);
86 double Qsq = sign * Q.m2Calc();
89 double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
96 Vec4 sum = radVec + recVec + emtVec;
97 double m2Dip = sum.m2Calc();
98 double x1 = 2. * (sum * radVec) / m2Dip;
99 double x3 = 2. * (sum * emtVec) / m2Dip;
101 pTnow = z * (1. - z);
105 Vec4 qBR(radVec - emtVec + recVec);
106 Vec4 qAR(radVec + recVec);
107 z = qBR.m2Calc() / qAR.m2Calc();
112 pTnow *= (Qsq - sign * m2Rad);
131 Vec4 p1 =
event[i1].p();
132 Vec4 p3 =
event[i3].p();
133 Vec4 p2 =
event[i2].p();
136 int iMoth1 =
event[i1].mother1();
137 int iMoth2 =
event[i2].mother1();
138 if (iMoth1 == 0 || iMoth2 == 0) {
139 loggerPtr->ABORT_MSG(
"could not find mothers of particles");
144 double mMoth1Sq =
event[iMoth1].m2();
145 double mMoth2Sq =
event[iMoth2].m2();
146 double sgn1 =
event[i1].isFinal() ? 1. : -1.;
147 double sgn2 =
event[i2].isFinal() ? 1. : -1.;
148 double qSq13 = sgn1*(
m2(sgn1*p1+p3) - mMoth1Sq);
149 double qSq23 = sgn2*(
m2(sgn2*p2+p3) - mMoth2Sq);
153 if (event[i1].isFinal() &&
event[i2].isFinal()) {
155 sMax =
m2(p1+p2+p3) - mMoth1Sq - mMoth2Sq;
156 }
else if ((event[i1].isResonance() && event[i2].isFinal())
157 || (!event[i1].isFinal() && event[i2].isFinal())) {
159 sMax = 2.*p1*p3 + 2.*p1*p2;
160 }
else if ((event[i1].isFinal() && event[i2].isResonance())
161 || (event[i1].isFinal() && !event[i2].isFinal())) {
163 sMax = 2.*p2*p3 + 2.*p1*p2;
164 }
else if (!event[i1].isFinal() || !event[i2].isFinal()) {
168 loggerPtr->ABORT_MSG(
"could not determine branching type");
173 double pT2now = qSq13*qSq23/sMax;
198 double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
199 ( e[iInA].e() + e[iInB].e() );
200 Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
201 iVecBst.bst(0., 0., betaZ);
202 jVecBst.
bst(0., 0., betaZ);
203 pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
204 iVecBst.e() * jVecBst.e() /
205 pow2(iVecBst.e() + jVecBst.e()) );
228 inline double pTcalc(
const Event &e,
int i,
int j,
int k,
int r,
int xSRin) {
231 double pTemt = -1., pTnow;
232 int xSR1 = (xSRin == -1) ? 0 : xSRin;
233 int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
234 for (
int xSR = xSR1; xSR < xSR2; xSR++) {
236 bool FSR = (xSR == 0) ?
false :
true;
240 if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
241 pTemt =
pTpowheg(e, i, j, (pTdefMode == 0) ?
false : FSR);
244 }
else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
245 pTemt =
pT(e, i, j, r, FSR);
248 }
else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
249 pTemt =
pT(e, k, j, r, FSR);
259 while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
260 while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
263 int jNow = (j > 0) ? j : 0;
264 int jMax = (j > 0) ? j + 1 : e.
size();
265 for (; jNow < jMax; jNow++) {
268 if (!e[jNow].isFinal())
continue;
270 if (QEDvetoMode==0 && e[jNow].colType() == 0)
continue;
273 if (pTdefMode == 0 || pTdefMode == 1) {
277 pTnow =
pTpowheg(e, iInA, jNow, (pTdefMode == 0) ?
false : FSR);
278 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
286 for (
int iMem = 0; iMem < outSize; iMem++) {
290 if (iNow == jNow )
continue;
292 if (QEDvetoMode==0 && e[iNow].colType() == 0)
continue;
293 if (jNow == e[iNow].daughter1()
294 && jNow == e[iNow].daughter2())
continue;
296 pTnow =
pTpowheg(e, iNow, jNow, (pTdefMode == 0)
298 if (pTnow > 0.) pTemt = (pTemt < 0)
299 ? pTnow : min(pTemt, pTnow);
305 }
else if (pTdefMode == 2) {
309 pTnow =
pT(e, iInA, jNow, iInB, FSR);
310 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
311 pTnow =
pT(e, iInB, jNow, iInA, FSR);
312 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
317 for (
int kNow = 0; kNow < e.
size(); kNow++) {
318 if (kNow == jNow || !e[kNow].isFinal())
continue;
319 if (QEDvetoMode==0 && e[kNow].colType() == 0)
continue;
323 pTnow =
pT(e, kNow, jNow, iInA, FSR);
324 if (pTnow > 0.) pTemt = (pTemt < 0)
325 ? pTnow : min(pTemt, pTnow);
326 pTnow =
pT(e, kNow, jNow, iInB, FSR);
327 if (pTnow > 0.) pTemt = (pTemt < 0)
328 ? pTnow : min(pTemt, pTnow);
331 for (
int rNow = 0; rNow < e.
size(); rNow++) {
332 if (rNow == kNow || rNow == jNow ||
333 !e[rNow].isFinal())
continue;
334 if(QEDvetoMode==0 && e[rNow].colType() == 0)
continue;
335 pTnow =
pT(e, kNow, jNow, rNow, FSR);
336 if (pTnow > 0.) pTemt = (pTemt < 0)
337 ? pTnow : min(pTemt, pTnow);
366 if (nMPI > 1)
return false;
368 double pT1(0.), pTsum(0.);
378 for (
int i = e.
size() - 1; i > 0; i--) {
379 if (e[i].isFinal()) {
386 if (count != nFinal && count != nFinal + 1) {
387 loggerPtr->ABORT_MSG(
"wrong number of final state particles in event");
391 isEmt = (count == nFinal) ?
false :
true;
392 iEmt = (isEmt) ? e.
size() - 1 : -1;
401 for (
int i = e.
size() - 1; i > 0; i--) {
402 if (e[i].isFinal()) {
403 if ( e[i].isParton() && iEmt < 0
404 && e[e[i].mother1()].isParton() ) {
415 if (!isEmt || pThardMode == 0) {
420 }
else if (pThardMode == 1) {
423 "pThard == 1 not available for nFinal == -1, reset pThard = 0.");
427 pThard =
pTcalc(e, -1, iEmt, -1, -1, -1);
433 }
else if (pThardMode == 2) {
436 "pThard == 2 not available for nFinal == -1, reset pThard = 0.");
440 pThard =
pTcalc(e, -1, -1, -1, -1, -1);
445 if (MPIvetoMode == 1) {
446 pTMPI = (isEmt) ? pTsum / 2. : pT1;
451 nAcceptSeq = nISRveto = nFSRveto = 0;
464 if (iSys != 0)
return false;
467 if (vetoMode == 1 && nAcceptSeq >= vetoCount)
return false;
470 int iRadAft = -1, iEmt = -1, iRecAft = -1;
471 for (
int i = e.
size() - 1; i > 0; i--) {
472 if (showerModel == 1) {
474 if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
475 else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
476 else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
477 }
else if (showerModel == 2) {
479 if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
480 else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
481 else if (iRecAft == -1
482 && (e[i].status() == -41 || e[i].status() == 44)) iRecAft = i;
484 if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
break;
486 if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
487 loggerPtr->ABORT_MSG(
"could not find ISR emission");
494 int xSR = (pTemtMode == 0) ? 0 : -1;
495 int i = (pTemtMode == 0) ? iRadAft : -1;
496 int j = (pTemtMode != 2) ? iEmt : -1;
498 int r = (pTemtMode == 0) ? iRecAft : -1;
499 double pTemt =
pTcalc(e, i, j, k, r, xSR);
503 bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
507 if (pTemt > pThard) {
510 nAcceptSeq = vetoCount-1;
531 if (iSys != 0)
return false;
534 if (vetoMode == 1 && nAcceptSeq >= vetoCount)
return false;
537 int iRecAft = e.
size() - 1;
538 int iEmt = e.
size() - 2;
539 int iRadAft = e.
size() - 3;
540 int iRadBef = e[iEmt].mother1();
542 if (showerModel == 2) {
544 if ( (e[iRecAft].status() != 51 && e[iRecAft].status() != 52) ||
545 e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop =
true;
548 if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
549 e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop =
true;
552 loggerPtr->ABORT_MSG(
"could not find FSR emission");
560 int xSR = (pTemtMode == 0) ? 1 : -1;
561 int i = (pTemtMode == 0) ? iRadBef : -1;
562 int k = (pTemtMode == 0) ? iRadAft : -1;
563 int r = (pTemtMode == 0) ? iRecAft : -1;
567 if (pTemtMode == 0 || pTemtMode == 1) {
574 if (emittedMode == 0 || (emittedMode == 2 &&
rndmPtr->
flat() < 0.5)) j++;
576 for (
int jLoop = 0; jLoop < 2; jLoop++) {
577 if (jLoop == 0) pTemt =
pTcalc(e, i, j, k, r, xSR);
578 else if (jLoop == 1) pTemt = min(pTemt,
pTcalc(e, i, j, k, r, xSR));
581 if (emittedMode != 3)
break;
582 if (k != -1) swap(j, k);
else j = iEmt;
586 }
else if (pTemtMode == 2) {
587 pTemt =
pTcalc(e, i, -1, k, r, xSR);
593 bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
597 if (pTemt > pThard) {
600 nAcceptSeq = vetoCount-1;
620 if (MPIvetoMode == 1) {
621 if (e[e.
size() - 1].pT() > pTMPI)
return true;
631 inline int getNFSRveto() {
return nFSRveto; }
636 double pThard, pTMPI;
639 int showerModel, nFinal, vetoMode, MPIvetoMode, QEDvetoMode, vetoCount;
640 int pThardMode, pTemtMode, emittedMode, pTdefMode;
641 bool accepted, isEmt;
646 unsigned long int nISRveto, nFSRveto;
655 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
double pTpowheg(const Event &e, int i, int j, bool FSR)
Compute the POWHEG pT separation between i and j.
Definition: PowhegHooks.h:188
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:173
bool initAfterBeams()
Initialize settings, detailing merging strategy to use.
Definition: PowhegHooks.h:39
Settings * settingsPtr
Pointer to the settings database.
Definition: PhysicsBase.h:85
The Event class holds all info on the generated event.
Definition: Event.h:408
bool canVetoMPIStep()
Definition: PowhegHooks.h:362
PartonSystems * partonSystemsPtr
Pointer to information on subcollision parton locations.
Definition: PhysicsBase.h:116
Rndm * rndmPtr
Pointer to the random number generator.
Definition: PhysicsBase.h:97
Use userhooks to veto PYTHIA emissions above the POWHEG scale.
Definition: PowhegHooks.h:27
ParticleData * particleDataPtr
Pointer to the particle data table.
Definition: PhysicsBase.h:88
bool canVetoMPIEmission()
MPI veto.
Definition: PowhegHooks.h:618
double pTvincia(const Event &event, int i1, int i3, int i2)
Definition: PowhegHooks.h:128
bool doVetoMPIStep(int nMPI, const Event &e)
Definition: PowhegHooks.h:364
int getNISRveto()
Functions to return information.
Definition: PowhegHooks.h:630
void bst(double betaX, double betaY, double betaZ)
Boost (simple).
Definition: Basics.cc:433
bool doVetoFSREmission(int, const Event &e, int iSys, bool)
Definition: PowhegHooks.h:529
bool doVetoISREmission(int, const Event &e, int iSys)
Definition: PowhegHooks.h:462
int size() const
Event record size.
Definition: Event.h:459
bool canVetoISREmission()
ISR veto.
Definition: PowhegHooks.h:461
double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR)
Definition: PowhegHooks.h:74
double m2(const Vec4 &v1)
The squared invariant mass of one or more four-vectors.
Definition: Basics.cc:597
Logger * loggerPtr
Pointer to logger.
Definition: PhysicsBase.h:91
bool doVetoMPIEmission(int, const Event &e)
Definition: PowhegHooks.h:619
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The Pythia class contains the top-level routines to generate an event.
Definition: Pythia.h:72
Info * infoPtr
Definition: PhysicsBase.h:82
UserHooks is base class for user access to program execution.
Definition: UserHooks.h:32
int numberVetoMPIStep()
Up to how many MPI steps should be checked.
Definition: PowhegHooks.h:363
double flat()
Generate next random number uniformly between 0 and 1.
Definition: Basics.cc:189
double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin)
Definition: PowhegHooks.h:228
bool canVetoFSREmission()
FSR veto.
Definition: PowhegHooks.h:528
double pT(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR)
Top-level wrapper for shower pTs.
Definition: PowhegHooks.h:61
PowhegHooks()
Constructors and destructor.
Definition: PowhegHooks.h:32
Definition: Settings.h:196