101 #include "Pythia8/Pythia.h"
104 #if PYTHIA_VERSION_INTEGER >= 8310
108 # include "Pythia8/Plugins.h"
109 # undef PYTHIA8_PLUGIN_CLASS
110 # undef PYTHIA8_PLUGIN_VERSIONS
111 # define PYTHIA8_PLUGIN_CLASS(BASE, CLASS, PYTHIA, SETTINGS, LOGGER)
112 # define PYTHIA8_PLUGIN_VERSIONS(...)
127 namespace Pythia8 {
class PowhegBB4Ldlsl; }
139 : m_debug(
"Powheg:bb4l:DEBUG", false),
140 m_vetoFSREmission(
"Powheg:bb4l:FSREmission:veto", true),
141 m_vetoQED(
"Powheg:bb4l:vetoQED", false),
142 m_topresscale(-1.), m_atopresscale(-1.),
143 m_wpresscale(-1.),m_wmresscale(-1.),
144 m_nInResonanceFSRveto(0),
145 m_pTmin(
"Powheg:bb4l:pTminVeto", 0.5),
146 m_vetoProduction(
"Powheg:veto", 1),
147 m_pTpythiaVeto(
"Powheg:bb4l:pTpythiaVeto", 0),
148 m_vetoDipoleFrame(
"Powheg:bb4l:FSREmission:vetoDipoleFrame", 0),
149 m_scaleResonanceVeto(
"Powheg:bb4l:ScaleResonance:veto", 0.),
150 m_pThardMode(
"POWHEG:pThard",0)
152 std::cout <<
"**********************************************************" << std::endl;
153 std::cout <<
"* *" << std::endl;
154 std::cout <<
"* Applying Powheg BB4L UserHook (PowhegBB4Ldlsl)! *" << std::endl;
155 std::cout <<
"* Must run on dedicated Powheg LHE input file *" << std::endl;
156 std::cout <<
"* (on your own responsibility) *" << std::endl;
157 std::cout <<
"* *" << std::endl;
158 std::cout <<
"**********************************************************" << std::endl;
180 ss << infoPtr->getEventComments();
183 assert (
temp ==
"#rwgt");
187 if (
radtype_.radtype==2)
return false;
189 int i_top = -1, i_atop = -1, i_wp = -1, i_wm = -1;
190 for (
int i = 0;
i <
e.size();
i++) {
191 if (
e[
i].
id() == 6) i_top =
i;
192 if (
e[
i].
id() == -6) i_atop =
i;
193 if (
e[
i].
id() == 24) i_wp =
i;
194 if (
e[
i].
id() == -24) i_wm =
i;
197 m_topresscale = findresscale(i_top,
e);
199 m_atopresscale = findresscale(i_atop,
e);
201 m_wpresscale = findresscale(i_wp,
e);
202 m_wmresscale = findresscale(i_wm,
e);
216 inline virtual bool canVetoFSREmission()
override {
return m_vetoFSREmission(settingsPtr) || m_vetoProduction(settingsPtr); }
220 if (inResonance && m_vetoFSREmission(settingsPtr)) {
223 int iRecAft =
e.size() - 1;
224 int iEmt =
e.size() - 2;
225 int iRadAft =
e.size() - 3;
226 int iRadBef =
e[iEmt].mother1();
229 int iRes =
e[iRadBef].mother1();
230 while ( iRes > 0 && (std::abs(
e[iRes].
id()) !=6 && std::abs(
e[iRes].
id()) != 24) ) {
231 iRes =
e[iRes].mother1();
234 #if PYTHIA_VERSION_INTEGER >= 8310
235 loggerPtr->ERROR_MSG(
"Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from the top quark or from the W boson, not vetoing");
237 infoPtr->errorMsg(
"Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from the top quark or from the W boson, not vetoing");
239 return doVetoFSR(
false);
241 int iResId =
e[iRes].id();
246 if(m_pTpythiaVeto(settingsPtr))
247 scale = pTpythia(
e, iRadAft, iEmt, iRecAft);
250 Vec4
pr(
e[iRadAft].
p()),
pe(
e[iEmt].
p()), pres(
e[iRes].
p()),
prec(
e[iRecAft].
p()), psystem;
254 if(m_vetoDipoleFrame(settingsPtr))
260 if (
e[iRadBef].
id() == 21)
261 scale = gSplittingScale(psystem,
pr,
pe);
263 else if (std::abs(
e[iRadBef].
id()) <= 5 && ((
e[iEmt].
id() == 21) && ! m_vetoQED(settingsPtr)) )
264 scale = qSplittingScale(psystem,
pr,
pe);
273 if ( m_debug(settingsPtr) &&
scale > m_topresscale)
274 cout << iResId <<
": " <<
e[iRadBef].id() <<
" > " <<
e[iRadAft].id() <<
" + " <<
e[iEmt].id() <<
"; " <<
scale << endl;
275 return doVetoFSR(
scale > m_topresscale);
277 else if (iResId == -6){
278 if ( m_debug(settingsPtr) &&
scale > m_atopresscale)
279 cout << iResId <<
": " <<
e[iRadBef].id() <<
" > " <<
e[iRadAft].id() <<
" + " <<
e[iEmt].id() <<
"; " <<
scale << endl;
280 return doVetoFSR(
scale > m_atopresscale);
282 else if (iResId == 24) {
283 if ( m_debug(settingsPtr) &&
scale > m_wpresscale)
284 cout << iResId <<
": " <<
e[iRadBef].id() <<
" > " <<
e[iRadAft].id() <<
" + " <<
e[iEmt].id() <<
"; " <<
scale << endl;
285 return doVetoFSR(
scale > m_wpresscale);
287 else if (iResId == -24){
288 if ( m_debug(settingsPtr) &&
scale > m_wmresscale)
289 cout << iResId <<
": " <<
e[iRadBef].id() <<
" > " <<
e[iRadAft].id() <<
" + " <<
e[iEmt].id() <<
"; " <<
scale << endl;
290 return doVetoFSR(
scale > m_wmresscale);
293 #if PYTHIA_VERSION_INTEGER >= 8310
294 loggerPtr->ERROR_MSG(
"Error in PowhegHooksBB4L::doVetoFSREmission: unimplemented case");
296 infoPtr->errorMsg(
"Error in PowhegHooksBB4L::doVetoFSREmission: unimplemented case");
302 else if(!inResonance && m_vetoProduction(settingsPtr)){
312 if (
radtype_.radtype==2)
return false;
314 m_nInResonanceFSRveto++;
331 return sqrt(
e[iRes].m2Calc());
333 if (
e[iRes].
id() == 6)
334 return m_topresscale;
335 else if (
e[iRes].
id() == -6)
336 return m_atopresscale;
337 else if (
e[iRes].
id() == 24)
339 else if (
e[iRes].
id() == -24)
352 if (iRes < 0)
return 1e30;
355 int nDau =
event[iRes].daughterList().size();
365 return m_pTmin(settingsPtr);
368 else if (std::abs(
event[iRes].
id()) == 6) {
370 int idw = -1,
idb = -1, idg = -1;
371 for (
int i = 0;
i < nDau;
i++) {
372 int iDau =
event[iRes].daughterList()[
i];
373 if (std::abs(
event[iDau].
id()) == 24) idw = iDau;
374 if (std::abs(
event[iDau].
id()) == 5)
idb = iDau;
375 if (std::abs(
event[iDau].
id()) == 21) idg = iDau;
380 pw.bstback(
event[iRes].
p());
390 else if (std::abs(
event[iRes].
id()) == 24) {
392 int idq = -1, ida = -1, idg = -1;
393 for (
int i = 0;
i < nDau;
i++) {
394 int iDau =
event[iRes].daughterList()[
i];
395 if (
event[iDau].
id() == 21) idg = iDau;
396 else if (
event[iDau].
id() > 0) idq = iDau;
397 else if (
event[iDau].
id() < 0) ida = iDau;
399 if (idq<0 or ida <0){
400 throw std::out_of_range(
"idq or ida out of range in PowhegHooksBB4Ldlsl.cxx");
404 pq.bstback(
event[iRes].
p());
406 pa.bstback(
event[iRes].
p());
411 Vec4 pw = pq + pa +
pg;
413 double csi = 2*
pg.e()/sqrt(q2);
414 double yq = 1 -
pg*pq/(
pg.e()*pq.e());
415 double ya = 1 -
pg*pa/(
pg.e()*pa.e());
417 return sqrt(
min(1-yq,1-ya)*pow2(csi)*q2/2);
428 inline bool match_decay(
int iparticle,
const Event &
e,
const vector<int> &
ids, vector<int> &positions, vector<Vec4> &momenta,
bool exitOnExtraLegs =
true){
430 if (
e[iparticle].daughterList().
size() !=
ids.size()) {
431 if (exitOnExtraLegs &&
e[iparticle].daughterList().
size() >
ids.size()) {
432 cout <<
"extra leg" << endl;
438 for (
long unsigned int i = 0;
i <
e[iparticle].daughterList().
size();
i++) {
439 int di =
e[iparticle].daughterList()[
i];
447 for (
long unsigned int i = 0;
i <
e[iparticle].daughterList().
size();
i++) {
448 int di =
e[iparticle].daughterList()[
i];
449 positions.push_back(di);
450 momenta.push_back(
e[di].
p());
458 return sqrt( 2*
p1*
p2*
p2.e()/
p1.e() );
476 Vec4 radVec =
e[RadAfterBranch].p();
477 Vec4 emtVec =
e[EmtAfterBranch].p();
478 Vec4 recVec =
e[RecAfterBranch].p();
479 int radID =
e[RadAfterBranch].id();
482 Vec4 Q(radVec + emtVec);
483 double Qsq = Q.m2Calc();
486 double m2Rad = (std::abs(radID) >= 4 && std::abs(radID) < 7) ?
487 pow2(particleDataPtr->m0(radID)) : 0.;
492 Vec4
sum = radVec + recVec + emtVec;
493 double m2Dip =
sum.m2Calc();
495 double x1 = 2. * (
sum * radVec) / m2Dip;
496 double x3 = 2. * (
sum * emtVec) / m2Dip;
498 pTnow =
z * (1. -
z);
502 pTnow *= (Qsq - m2Rad);
505 cout <<
"Warning: pTpythia was negative" << endl;
531 if (m_pThardMode(settingsPtr) == 1 || m_pThardMode(settingsPtr) == 2){
533 if (nMPI > 1)
return false;
536 bool top_exists =
false;
537 bool antitop_exists =
false;
538 for (
int i =
e.size() - 1;
i > 0;
i--) {
539 if (
e[
i].isFinal()) {
542 if (
e[
i].
id() == 6) top_exists =
true;
543 if (
e[
i].
id() == -6) antitop_exists =
true;
548 int nCurrentFinal = -1;
550 if (top_exists && antitop_exists) nCurrentFinal = 2;
552 else if (top_exists || antitop_exists) nCurrentFinal = 3;
554 std::cout <<
"Error: neither top nor anti-top found - something went wrong" << std::endl;
558 if (
count != nCurrentFinal &&
count != nCurrentFinal + 1) {
559 cout <<
"Error: wrong number of final state particles in event: count " <<
count <<
" nCurrentFinal: "<< nCurrentFinal << endl;
560 for (
int i =
e.size() - 1;
i > 0;
i--) {
561 if (
e[
i].isFinal()) { std:: cout <<
"FS particle " <<
i <<
" PID: " <<
e[
i].id() << std::endl;
567 bool isEmt = (
count == nCurrentFinal) ?
false :
true;
568 int iEmt = (isEmt) ?
e.size() - 1 : -1;
571 if (m_pThardMode(settingsPtr) == 1) {
576 }
else if (m_pThardMode(settingsPtr) == 2) {
579 std::cout <<
"Error: m_pThardMode neither 1 or 2 - something went wrong" << std::endl;