17 #ifndef Pythia8_PowhegHooksSetterMethod_H
18 #define Pythia8_PowhegHooksSetterMethod_H
21 #include "Pythia8/Pythia.h"
22 #include "Pythia8/Plugins.h"
81 m_nFinal = settingsPtr->mode(
"POWHEG:nFinal");
83 m_vetoCount = settingsPtr->mode(
"POWHEG:vetoCount");
102 inline double pT(
const Event&
e,
int RadAfterBranch,
103 int EmtAfterBranch,
int RecAfterBranch,
bool FSR) {
106 return pTvincia(
e, RadAfterBranch, EmtAfterBranch, RecAfterBranch);
109 return pTdire(
e, RadAfterBranch, EmtAfterBranch, RecAfterBranch);
110 return pTpythia(
e, RadAfterBranch, EmtAfterBranch, RecAfterBranch,
FSR);
118 int EmtAfterBranch,
int RecAfterBranch,
bool FSR) {
121 Vec4 radVec =
e[RadAfterBranch].p();
122 Vec4 emtVec =
e[EmtAfterBranch].p();
123 Vec4 recVec =
e[RecAfterBranch].p();
124 int radID =
e[RadAfterBranch].id();
127 double sign = (
FSR) ? 1. : -1.;
128 Vec4 Q(radVec +
sign * emtVec);
129 double Qsq =
sign * Q.m2Calc();
132 double m2Rad = (std::abs(radID) >= 4 && std::abs(radID) < 7) ?
133 pow2(particleDataPtr->m0(radID)) : 0.;
139 Vec4
sum = radVec + recVec + emtVec;
140 double m2Dip =
sum.m2Calc();
141 double x1 = 2. * (
sum * radVec) / m2Dip;
142 double x3 = 2. * (
sum * emtVec) / m2Dip;
144 pTnow =
z * (1. -
z);
148 Vec4 qBR(radVec - emtVec + recVec);
149 Vec4 qAR(radVec + recVec);
150 z = qBR.m2Calc() / qAR.m2Calc();
155 pTnow *= (Qsq -
sign * m2Rad);
159 loggerPtr->WARNING_MSG(
"negative pT");
174 Vec4
p1 =
event[i1].p();
175 Vec4
p3 =
event[i3].p();
176 Vec4
p2 =
event[i2].p();
179 int iMoth1 =
event[i1].mother1();
180 int iMoth2 =
event[i2].mother1();
181 if (iMoth1 == 0 || iMoth2 == 0) {
182 loggerPtr->ABORT_MSG(
"could not find mothers of particles");
187 double mMoth1Sq =
event[iMoth1].m2();
188 double mMoth2Sq =
event[iMoth2].m2();
189 double sgn1 =
event[i1].isFinal() ? 1. : -1.;
190 double sgn2 =
event[i2].isFinal() ? 1. : -1.;
191 double qSq13 = sgn1*(
m2(sgn1*
p1+
p3) - mMoth1Sq);
192 double qSq23 = sgn2*(
m2(sgn2*
p2+
p3) - mMoth2Sq);
196 if (
event[i1].isFinal() &&
event[i2].isFinal()) {
198 sMax =
m2(
p1+
p2+
p3) - mMoth1Sq - mMoth2Sq;
200 || (!
event[i1].isFinal() &&
event[i2].isFinal())) {
204 || (
event[i1].isFinal() && !
event[i2].isFinal())) {
207 }
else if (!
event[i1].isFinal() || !
event[i2].isFinal()) {
211 loggerPtr->ABORT_MSG(
"could not determine branching type");
216 double pT2now = qSq13*qSq23/sMax;
220 loggerPtr->WARNING_MSG(
"negative pT");
242 if (
rad.isFinal() && rec.isFinal()) {
244 const double sij = 2.*
rad.p()*emt.p();
245 const double sik = 2.*
rad.p()*rec.p();
246 const double sjk = 2.*rec.p()*emt.p();
247 pT2 = sij*sjk/(sij+sik+sjk);
248 }
else if (
rad.isFinal() && !rec.isFinal()) {
250 const double sij = 2.*
rad.p()*emt.p();
251 const double sai = -2.*rec.p()*
rad.p();
252 const double saj = -2.*rec.p()*emt.p();
253 pT2 = sij*saj/(sai+saj)*(sij+saj+sai)/(sai+saj);
254 if (sij+saj+sai < 1
e-5 && std::abs(sij+saj+sai) < 1
e-5) pT2 = sij;
255 }
else if (!
rad.isFinal() && rec.isFinal()) {
257 const double sai = -2.*
rad.p()*emt.p();
258 const double sik = 2.*rec.p()*emt.p();
259 const double sak = -2.*
rad.p()*rec.p();
260 pT2 = sai*sik/(sai+sak)*(sai+sik+sak)/(sai+sak);
261 }
else if (!
rad.isFinal() || !rec.isFinal()) {
263 const double sai = -2.*
rad.p()*emt.p();
264 const double sbi = -2.*rec.p()*emt.p();
265 const double sab = 2.*
rad.p()*rec.p();
266 pT2 = sai*sbi/sab*(sai+sbi+sab)/sab;
268 loggerPtr->ABORT_MSG(
"could not determine branching type");
274 loggerPtr->WARNING_MSG(
"negative pT");
293 int iInA = partonSystemsPtr->getInA(0);
294 int iInB = partonSystemsPtr->getInB(0);
295 double betaZ = - (
e[iInA].pz() +
e[iInB].pz() ) /
296 (
e[iInA].
e() +
e[iInB].e() );
297 Vec4 iVecBst(
e[
i].
p()), jVecBst(
e[j].
p());
298 iVecBst.bst(0., 0., betaZ);
299 jVecBst.bst(0., 0., betaZ);
300 pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
301 iVecBst.e() * jVecBst.e() /
302 pow2(iVecBst.e() + jVecBst.e()) );
311 loggerPtr->WARNING_MSG(
"negative pT");
328 double pTemt = -1., pTnow;
329 int xSR1 = (xSRin == -1) ? 0 : xSRin;
330 int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
331 for (
int xSR = xSR1; xSR < xSR2; xSR++) {
333 bool FSR = (xSR == 0) ?
false :
true;
354 int iInA = partonSystemsPtr->getInA(0);
355 int iInB = partonSystemsPtr->getInB(0);
356 while (
e[iInA].mother1() != 1) { iInA =
e[iInA].mother1(); }
357 while (
e[iInB].mother1() != 2) { iInB =
e[iInB].mother1(); }
360 int jNow = (j > 0) ? j : 0;
361 int jMax = (j > 0) ? j + 1 :
e.size();
362 for (; jNow < jMax; jNow++) {
365 if (!
e[jNow].isFinal())
continue;
375 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow :
min(pTemt, pTnow);
382 int outSize = partonSystemsPtr->sizeOut(0);
383 for (
int iMem = 0; iMem < outSize; iMem++) {
384 int iNow = partonSystemsPtr->getOut(0, iMem);
387 if (iNow == jNow )
continue;
390 if (jNow ==
e[iNow].daughter1()
391 && jNow ==
e[iNow].daughter2())
continue;
395 if (pTnow > 0.) pTemt = (pTemt < 0)
396 ? pTnow :
min(pTemt, pTnow);
406 pTnow =
pT(
e, iInA, jNow, iInB,
FSR);
407 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow :
min(pTemt, pTnow);
408 pTnow =
pT(
e, iInB, jNow, iInA,
FSR);
409 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow :
min(pTemt, pTnow);
414 for (
int kNow = 0; kNow <
e.size(); kNow++) {
415 if (kNow == jNow || !
e[kNow].isFinal())
continue;
420 pTnow =
pT(
e, kNow, jNow, iInA,
FSR);
421 if (pTnow > 0.) pTemt = (pTemt < 0)
422 ? pTnow :
min(pTemt, pTnow);
423 pTnow =
pT(
e, kNow, jNow, iInB,
FSR);
424 if (pTnow > 0.) pTemt = (pTemt < 0)
425 ? pTnow :
min(pTemt, pTnow);
428 for (
int rNow = 0; rNow <
e.size(); rNow++) {
429 if (rNow == kNow || rNow == jNow ||
430 !
e[rNow].isFinal())
continue;
432 pTnow =
pT(
e, kNow, jNow, rNow,
FSR);
433 if (pTnow > 0.) pTemt = (pTemt < 0)
434 ? pTnow :
min(pTemt, pTnow);
463 if (nMPI > 1)
return false;
465 double pT1(0.), pTsum(0.);
475 for (
int i =
e.size() - 1;
i > 0;
i--) {
476 if (
e[
i].isFinal()) {
484 loggerPtr->ABORT_MSG(
"wrong number of final state particles in event");
489 iEmt = (
m_isEmt) ?
e.size() - 1 : -1;
498 for (
int i =
e.size() - 1;
i > 0;
i--) {
499 if (
e[
i].isFinal()) {
501 &&
e[
e[
i].mother1()].isParton() ) {
519 loggerPtr->WARNING_MSG(
520 "pThard == 1 not available for m_nFinal == -1, reset pThard = 0. (But might be reset later by child class)");
532 loggerPtr->WARNING_MSG(
533 "pThard == 2 not available for m_nFinal == -1, reset pThard = 0. (But might be reset later by child class)");
561 if (iSys != 0)
return false;
567 int iRadAft = -1, iEmt = -1, iRecAft = -1;
568 for (
int i =
e.size() - 1;
i > 0;
i--) {
571 if (iRadAft == -1 &&
e[
i].
status() == -41) iRadAft =
i;
572 else if (iEmt == -1 &&
e[
i].
status() == 43) iEmt =
i;
573 else if (iRecAft == -1 &&
e[
i].
status() == -42) iRecAft =
i;
576 if (iRadAft == -1 &&
e[
i].
status() == -41) iRadAft =
i;
577 else if (iEmt == -1 &&
e[
i].
status() == 43) iEmt =
i;
578 else if (iRecAft == -1
582 if (iRadAft == -1 &&
e[
i].
status() == -41) iRadAft =
i;
583 else if (iEmt == -1 &&
e[
i].
status() == 43) iEmt =
i;
584 else if (iRecAft == -1
588 if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
break;
590 if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
591 loggerPtr->ABORT_MSG(
"could not find ISR emission");
635 if (iSys != 0)
return false;
641 int iRecAft =
e.size() - 1;
642 int iEmt =
e.size() - 2;
643 int iRadAft =
e.size() - 3;
644 int iRadBef =
e[iEmt].mother1();
656 loggerPtr->ABORT_MSG(
"could not find FSR emission");
680 for (
int jLoop = 0; jLoop < 2; jLoop++) {
681 if (jLoop == 0) pTemt =
pTcalc(
e,
i, j,
k,
r, xSR);
682 else if (jLoop == 1) pTemt =
min(pTemt,
pTcalc(
e,
i, j,
k,
r, xSR));
686 if (
k != -1)
swap(j,
k);
else j = iEmt;
725 if (
e[
e.size() - 1].pT() >
m_pTMPI)
return true;
762 #endif // end Pythia8_PowhegHooksSetterMethod_H