7 using namespace Pythia8;
18 PowhegHooks(
int nFinalIn,
int vetoModeIn,
int vetoCountIn,
19 int pThardModeIn,
int pTemtModeIn,
int emittedModeIn,
20 int pTdefModeIn,
int MPIvetoModeIn) : nFinal(nFinalIn),
21 vetoMode(vetoModeIn), vetoCount(vetoCountIn),
22 pThardMode(pThardModeIn), pTemtMode(pTemtModeIn),
23 emittedMode(emittedModeIn), pTdefMode(pTdefModeIn),
24 MPIvetoMode(MPIvetoModeIn) {};
35 double pTpythia(
const Event &e,
int RadAfterBranch,
int EmtAfterBranch,
36 int RecAfterBranch,
bool FSR) {
39 Vec4 radVec = e[RadAfterBranch].p();
40 Vec4 emtVec = e[EmtAfterBranch].p();
41 Vec4 recVec = e[RecAfterBranch].p();
42 int radID = e[RadAfterBranch].id();
45 double sign = (FSR) ? 1. : -1.;
46 Vec4 Q(radVec + sign * emtVec);
47 double Qsq = sign * Q.m2Calc();
50 double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
51 pow2(particleDataPtr->m0(radID)) : 0.;
57 Vec4 sum = radVec + recVec + emtVec;
58 double m2Dip = sum.m2Calc();
59 double x1 = 2. * (sum * radVec) / m2Dip;
60 double x3 = 2. * (sum * emtVec) / m2Dip;
66 Vec4 qBR(radVec - emtVec + recVec);
67 Vec4 qAR(radVec + recVec);
68 z = qBR.m2Calc() / qAR.m2Calc();
73 pTnow *= (Qsq - sign * m2Rad);
77 cout <<
"Warning: pTpythia was negative" << endl;
82 cout <<
"pTpythia: rad = " << RadAfterBranch <<
", emt = "
83 << EmtAfterBranch <<
", rec = " << RecAfterBranch
84 <<
", pTnow = " << sqrt(pTnow) << endl;
92 double pTpowheg(
const Event &e,
int i,
int j,
bool FSR) {
100 int iInA = partonSystemsPtr->getInA(0);
101 int iInB = partonSystemsPtr->getInB(0);
102 double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
103 ( e[iInA].e() + e[iInB].e() );
104 Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
105 iVecBst.bst(0., 0., betaZ);
106 jVecBst.bst(0., 0., betaZ);
107 pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
108 iVecBst.e() * jVecBst.e() /
109 pow2(iVecBst.e() + jVecBst.e()) );
118 cout <<
"Warning: pTpowheg was negative" << endl;
123 cout <<
"pTpowheg: i = " << i <<
", j = " << j
124 <<
", pTnow = " << pTnow << endl;
135 double pTcalc(
const Event &e,
int i,
int j,
int k,
int r,
int xSRin) {
138 double pTemt = -1., pTnow;
139 int xSR1 = (xSRin == -1) ? 0 : xSRin;
140 int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
141 for (
int xSR = xSR1; xSR < xSR2; xSR++) {
143 bool FSR = (xSR == 0) ?
false :
true;
147 if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
148 pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ?
false : FSR);
151 }
else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
152 pTemt = pTpythia(e, i, j, r, FSR);
155 }
else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
156 pTemt = pTpythia(e, k, j, r, FSR);
164 int iInA = partonSystemsPtr->getInA(0);
165 int iInB = partonSystemsPtr->getInB(0);
166 while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
167 while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
170 int jNow = (j > 0) ? j : 0;
171 int jMax = (j > 0) ? j + 1 : e.size();
172 for (; jNow < jMax; jNow++) {
175 if (!e[jNow].isFinal() || e[jNow].colType() == 0)
continue;
178 if (pTdefMode == 0 || pTdefMode == 1) {
182 pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ?
false : FSR);
183 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
190 int outSize = partonSystemsPtr->sizeOut(0);
191 for (
int iMem = 0; iMem < outSize; iMem++) {
192 int iNow = partonSystemsPtr->getOut(0, iMem);
195 if (iNow == jNow || e[iNow].colType() == 0)
continue;
196 if (jNow == e[iNow].daughter1()
197 && jNow == e[iNow].daughter2())
continue;
199 pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0)
201 if (pTnow > 0.) pTemt = (pTemt < 0)
202 ? pTnow : min(pTemt, pTnow);
208 }
else if (pTdefMode == 2) {
212 pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
213 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
214 pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
215 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
220 for (
int kNow = 0; kNow < e.size(); kNow++) {
221 if (kNow == jNow || !e[kNow].isFinal() ||
222 e[kNow].colType() == 0)
continue;
226 pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
227 if (pTnow > 0.) pTemt = (pTemt < 0)
228 ? pTnow : min(pTemt, pTnow);
229 pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
230 if (pTnow > 0.) pTemt = (pTemt < 0)
231 ? pTnow : min(pTemt, pTnow);
234 for (
int rNow = 0; rNow < e.size(); rNow++) {
235 if (rNow == kNow || rNow == jNow ||
236 !e[rNow].isFinal() || e[rNow].colType() == 0)
continue;
237 pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
238 if (pTnow > 0.) pTemt = (pTemt < 0)
239 ? pTnow : min(pTemt, pTnow);
250 cout <<
"pTcalc: i = " << i <<
", j = " << j <<
", k = " << k
251 <<
", r = " << r <<
", xSR = " << xSRin
252 <<
", pTemt = " << pTemt << endl;
265 bool canVetoMPIStep() {
return true; }
266 int numberVetoMPIStep() {
return 1; }
267 bool doVetoMPIStep(
int nMPI,
const Event &e) {
269 if (nMPI > 1)
return false;
275 double pT1 = 0., pTsum = 0.;
276 for (
int i = e.size() - 1; i > 0; i--) {
277 if (e[i].isFinal()) {
284 if (count != nFinal && count != nFinal + 1) {
285 cout <<
"Error: wrong number of final state particles in event" << endl;
289 bool isEmt = (count == nFinal) ?
false :
true;
290 int iEmt = (isEmt) ? e.size() - 1 : -1;
293 if (!isEmt || pThardMode == 0) {
294 pThard = infoPtr->QFac();
298 }
else if (pThardMode == 1) {
299 pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
304 }
else if (pThardMode == 2) {
305 pThard = pTcalc(e, -1, -1, -1, -1, -1);
310 if (MPIvetoMode == 1) {
311 pTMPI = (isEmt) ? pTsum / 2. : pT1;
315 cout <<
"doVetoMPIStep: Qfac = " << infoPtr->QFac()
316 <<
", pThard = " << pThard << endl << endl;
321 nAcceptSeq = nISRveto = nFSRveto = 0;
331 bool canVetoISREmission() {
return (vetoMode == 0) ?
false :
true; }
332 bool doVetoISREmission(
int,
const Event &e,
int iSys) {
334 if (iSys != 0)
return false;
337 if (vetoMode == 1 && nAcceptSeq >= vetoCount)
return false;
340 int iRadAft = -1, iEmt = -1, iRecAft = -1;
341 for (
int i = e.size() - 1; i > 0; i--) {
342 if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
343 else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
344 else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
345 if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
break;
347 if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
349 cout <<
"Error: couldn't find Pythia ISR emission" << endl;
356 int xSR = (pTemtMode == 0) ? 0 : -1;
357 int i = (pTemtMode == 0) ? iRadAft : -1;
358 int j = (pTemtMode != 2) ? iEmt : -1;
360 int r = (pTemtMode == 0) ? iRecAft : -1;
361 double pTemt = pTcalc(e, i, j, k, r, xSR);
364 cout <<
"doVetoISREmission: pTemt = " << pTemt << endl << endl;
368 if (pTemt > pThard) {
384 bool canVetoFSREmission() {
return (vetoMode == 0) ?
false :
true; }
385 bool doVetoFSREmission(
int,
const Event &e,
int iSys,
bool) {
387 if (iSys != 0)
return false;
390 if (vetoMode == 1 && nAcceptSeq >= vetoCount)
return false;
393 int iRecAft = e.size() - 1;
394 int iEmt = e.size() - 2;
395 int iRadAft = e.size() - 3;
396 int iRadBef = e[iEmt].mother1();
397 if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
398 e[iEmt].status() != 51 || e[iRadAft].status() != 51) {
400 cout <<
"Error: couldn't find Pythia FSR emission" << endl;
408 int xSR = (pTemtMode == 0) ? 1 : -1;
409 int i = (pTemtMode == 0) ? iRadBef : -1;
410 int k = (pTemtMode == 0) ? iRadAft : -1;
411 int r = (pTemtMode == 0) ? iRecAft : -1;
415 if (pTemtMode == 0 || pTemtMode == 1) {
422 if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
424 for (
int jLoop = 0; jLoop < 2; jLoop++) {
425 if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
426 else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
429 if (emittedMode != 3)
break;
430 if (k != -1) swap(j, k);
else j = iEmt;
434 }
else if (pTemtMode == 2) {
435 pTemt = pTcalc(e, i, -1, k, r, xSR);
440 cout <<
"doVetoFSREmission: pTemt = " << pTemt << endl << endl;
444 if (pTemt > pThard) {
460 bool canVetoMPIEmission() {
return (MPIvetoMode == 0) ?
false :
true; }
461 bool doVetoMPIEmission(
int,
const Event &e) {
462 if (MPIvetoMode == 1) {
463 if (e[e.size() - 1].pT() > pTMPI) {
465 cout <<
"doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT()
466 <<
", pTMPI = " << pTMPI << endl << endl;
478 int getNISRveto() {
return nISRveto; }
479 int getNFSRveto() {
return nFSRveto; }
482 int nFinal, vetoMode, vetoCount, pThardMode, pTemtMode,
483 emittedMode, pTdefMode, MPIvetoMode;
484 double pThard, pTMPI;
489 unsigned long int nISRveto, nFSRveto;
495 int main(
int,
char **) {
501 pythia.settings.addMode(
"POWHEG:nFinal", 2,
true,
false, 1, 0);
502 pythia.settings.addMode(
"POWHEG:veto", 0,
true,
true, 0, 2);
503 pythia.settings.addMode(
"POWHEG:vetoCount", 3,
true,
false, 0, 0);
504 pythia.settings.addMode(
"POWHEG:pThard", 0,
true,
true, 0, 2);
505 pythia.settings.addMode(
"POWHEG:pTemt", 0,
true,
true, 0, 2);
506 pythia.settings.addMode(
"POWHEG:emitted", 0,
true,
true, 0, 3);
507 pythia.settings.addMode(
"POWHEG:pTdef", 0,
true,
true, 0, 2);
508 pythia.settings.addMode(
"POWHEG:MPIveto", 0,
true,
true, 0, 1);
511 pythia.readFile(
"main31.cmnd");
514 int nEvent = pythia.settings.mode(
"Main:numberOfEvents");
515 int nError = pythia.settings.mode(
"Main:timesAllowErrors");
517 int nFinal = pythia.settings.mode(
"POWHEG:nFinal");
518 int vetoMode = pythia.settings.mode(
"POWHEG:veto");
519 int vetoCount = pythia.settings.mode(
"POWHEG:vetoCount");
520 int pThardMode = pythia.settings.mode(
"POWHEG:pThard");
521 int pTemtMode = pythia.settings.mode(
"POWHEG:pTemt");
522 int emittedMode = pythia.settings.mode(
"POWHEG:emitted");
523 int pTdefMode = pythia.settings.mode(
"POWHEG:pTdef");
524 int MPIvetoMode = pythia.settings.mode(
"POWHEG:MPIveto");
525 bool loadHooks = (vetoMode > 0 || MPIvetoMode > 0);
533 pythia.readString(
"SpaceShower:pTmaxMatch = 2");
534 pythia.readString(
"TimeShower:pTmaxMatch = 2");
538 if (MPIvetoMode > 0) {
539 pythia.readString(
"MultipartonInteractions:pTmaxMatch = 2");
542 powhegHooks =
new PowhegHooks(nFinal, vetoMode, vetoCount,
543 pThardMode, pTemtMode, emittedMode, pTdefMode, MPIvetoMode);
544 pythia.setUserHooksPtr((
UserHooks *) powhegHooks);
551 unsigned long int nISRveto = 0, nFSRveto = 0;
555 int iEvent = 0, iError = 0;
559 if (!pythia.next()) {
562 if (pythia.info.atEndOfFile())
break;
565 cout <<
"Warning: event " << iEvent <<
" failed" << endl;
566 if (++iError == nError) {
567 cout <<
"Error: too many event failures.. exiting" << endl;
580 nISRveto += powhegHooks->getNISRveto();
581 nFSRveto += powhegHooks->getNFSRveto();
586 if (nEvent != 0 && iEvent == nEvent)
break;
592 cout <<
"Number of ISR emissions vetoed: " << nISRveto << endl;
593 cout <<
"Number of FSR emissions vetoed: " << nFSRveto << endl;