File indexing completed on 2026-03-31 08:26:53
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #ifndef Pythia8_PowhegHooks_H
0015 #define Pythia8_PowhegHooks_H
0016
0017
0018 #include "Pythia8/Pythia.h"
0019 #include "Pythia8/Plugins.h"
0020
0021 namespace Pythia8 {
0022
0023
0024
0025
0026
0027 class PowhegHooks : public UserHooks {
0028
0029 public:
0030
0031
0032 PowhegHooks() {}
0033 PowhegHooks(Pythia*, Settings*, Logger*) {}
0034 ~PowhegHooks() {}
0035
0036
0037
0038
0039 bool initAfterBeams() {
0040 nFinal = settingsPtr->mode("POWHEG:nFinal");
0041 vetoMode = settingsPtr->mode("POWHEG:veto");
0042 vetoCount = settingsPtr->mode("POWHEG:vetoCount");
0043 pThardMode = settingsPtr->mode("POWHEG:pThard");
0044 pTemtMode = settingsPtr->mode("POWHEG:pTemt");
0045 emittedMode = settingsPtr->mode("POWHEG:emitted");
0046 pTdefMode = settingsPtr->mode("POWHEG:pTdef");
0047 MPIvetoMode = settingsPtr->mode("POWHEG:MPIveto");
0048 QEDvetoMode = settingsPtr->mode("POWHEG:QEDveto");
0049 showerModel = settingsPtr->mode("PartonShowers:model");
0050 return true;
0051 }
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061 inline double pT(const Event& e, int RadAfterBranch,
0062 int EmtAfterBranch, int RecAfterBranch, bool FSR) {
0063
0064 if (showerModel == 2)
0065 return pTvincia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch);
0066
0067 return pTpythia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch, FSR);
0068 }
0069
0070
0071
0072
0073
0074 inline double pTpythia(const Event& e, int RadAfterBranch,
0075 int EmtAfterBranch, int RecAfterBranch, bool FSR) {
0076
0077
0078 Vec4 radVec = e[RadAfterBranch].p();
0079 Vec4 emtVec = e[EmtAfterBranch].p();
0080 Vec4 recVec = e[RecAfterBranch].p();
0081 int radID = e[RadAfterBranch].id();
0082
0083
0084 double sign = (FSR) ? 1. : -1.;
0085 Vec4 Q(radVec + sign * emtVec);
0086 double Qsq = sign * Q.m2Calc();
0087
0088
0089 double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
0090 pow2(particleDataPtr->m0(radID)) : 0.;
0091
0092
0093 double z, pTnow;
0094 if (FSR) {
0095
0096 Vec4 sum = radVec + recVec + emtVec;
0097 double m2Dip = sum.m2Calc();
0098 double x1 = 2. * (sum * radVec) / m2Dip;
0099 double x3 = 2. * (sum * emtVec) / m2Dip;
0100 z = x1 / (x1 + x3);
0101 pTnow = z * (1. - z);
0102
0103 } else {
0104
0105 Vec4 qBR(radVec - emtVec + recVec);
0106 Vec4 qAR(radVec + recVec);
0107 z = qBR.m2Calc() / qAR.m2Calc();
0108 pTnow = (1. - z);
0109 }
0110
0111
0112 pTnow *= (Qsq - sign * m2Rad);
0113
0114
0115 if (pTnow < 0.) {
0116 loggerPtr->WARNING_MSG("negative pT");
0117 return -1.;
0118 }
0119
0120
0121 return sqrt(pTnow);
0122 }
0123
0124
0125
0126
0127
0128 inline double pTvincia(const Event& event, int i1, int i3, int i2) {
0129
0130
0131 Vec4 p1 = event[i1].p();
0132 Vec4 p3 = event[i3].p();
0133 Vec4 p2 = event[i2].p();
0134
0135
0136 int iMoth1 = event[i1].mother1();
0137 int iMoth2 = event[i2].mother1();
0138 if (iMoth1 == 0 || iMoth2 == 0) {
0139 loggerPtr->ABORT_MSG("could not find mothers of particles");
0140 exit(1);
0141 }
0142
0143
0144 double mMoth1Sq = event[iMoth1].m2();
0145 double mMoth2Sq = event[iMoth2].m2();
0146 double sgn1 = event[i1].isFinal() ? 1. : -1.;
0147 double sgn2 = event[i2].isFinal() ? 1. : -1.;
0148 double qSq13 = sgn1*(m2(sgn1*p1+p3) - mMoth1Sq);
0149 double qSq23 = sgn2*(m2(sgn2*p2+p3) - mMoth2Sq);
0150
0151
0152 double sMax = -1.;
0153 if (event[i1].isFinal() && event[i2].isFinal()) {
0154
0155 sMax = m2(p1+p2+p3) - mMoth1Sq - mMoth2Sq;
0156 } else if ((event[i1].isResonance() && event[i2].isFinal())
0157 || (!event[i1].isFinal() && event[i2].isFinal())) {
0158
0159 sMax = 2.*p1*p3 + 2.*p1*p2;
0160 } else if ((event[i1].isFinal() && event[i2].isResonance())
0161 || (event[i1].isFinal() && !event[i2].isFinal())) {
0162
0163 sMax = 2.*p2*p3 + 2.*p1*p2;
0164 } else if (!event[i1].isFinal() || !event[i2].isFinal()) {
0165
0166 sMax = 2.*p1*p2;
0167 } else {
0168 loggerPtr->ABORT_MSG("could not determine branching type");
0169 exit(1);
0170 }
0171
0172
0173 double pT2now = qSq13*qSq23/sMax;
0174
0175
0176 if (pT2now < 0.) {
0177 loggerPtr->WARNING_MSG("negative pT");
0178 return -1.;
0179 }
0180
0181
0182 return sqrt(pT2now);
0183 }
0184
0185
0186
0187
0188 inline double pTpowheg(const Event &e, int i, int j, bool FSR) {
0189
0190
0191 double pTnow = 0.;
0192 if (FSR) {
0193
0194
0195
0196 int iInA = partonSystemsPtr->getInA(0);
0197 int iInB = partonSystemsPtr->getInB(0);
0198 double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
0199 ( e[iInA].e() + e[iInB].e() );
0200 Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
0201 iVecBst.bst(0., 0., betaZ);
0202 jVecBst.bst(0., 0., betaZ);
0203 pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
0204 iVecBst.e() * jVecBst.e() /
0205 pow2(iVecBst.e() + jVecBst.e()) );
0206
0207 } else {
0208
0209 pTnow = e[j].pT();
0210 }
0211
0212
0213 if (pTnow < 0.) {
0214 loggerPtr->WARNING_MSG("negative pT");
0215 return -1.;
0216 }
0217
0218 return pTnow;
0219 }
0220
0221
0222
0223
0224
0225
0226
0227
0228 inline double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) {
0229
0230
0231 double pTemt = -1., pTnow;
0232 int xSR1 = (xSRin == -1) ? 0 : xSRin;
0233 int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
0234 for (int xSR = xSR1; xSR < xSR2; xSR++) {
0235
0236 bool FSR = (xSR == 0) ? false : true;
0237
0238
0239
0240 if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
0241 pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
0242
0243
0244 } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
0245 pTemt = pT(e, i, j, r, FSR);
0246
0247
0248 } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
0249 pTemt = pT(e, k, j, r, FSR);
0250
0251
0252 } else {
0253
0254
0255
0256
0257 int iInA = partonSystemsPtr->getInA(0);
0258 int iInB = partonSystemsPtr->getInB(0);
0259 while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
0260 while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
0261
0262
0263 int jNow = (j > 0) ? j : 0;
0264 int jMax = (j > 0) ? j + 1 : e.size();
0265 for (; jNow < jMax; jNow++) {
0266
0267
0268 if (!e[jNow].isFinal()) continue;
0269
0270 if (QEDvetoMode==0 && e[jNow].colType() == 0) continue;
0271
0272
0273 if (pTdefMode == 0 || pTdefMode == 1) {
0274
0275
0276 if (!FSR) {
0277 pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
0278 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
0279
0280
0281
0282
0283 } else {
0284
0285 int outSize = partonSystemsPtr->sizeOut(0);
0286 for (int iMem = 0; iMem < outSize; iMem++) {
0287 int iNow = partonSystemsPtr->getOut(0, iMem);
0288
0289
0290 if (iNow == jNow ) continue;
0291
0292 if (QEDvetoMode==0 && e[iNow].colType() == 0) continue;
0293 if (jNow == e[iNow].daughter1()
0294 && jNow == e[iNow].daughter2()) continue;
0295
0296 pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0)
0297 ? false : FSR);
0298 if (pTnow > 0.) pTemt = (pTemt < 0)
0299 ? pTnow : min(pTemt, pTnow);
0300 }
0301
0302 }
0303
0304
0305 } else if (pTdefMode == 2) {
0306
0307
0308 if (!FSR) {
0309 pTnow = pT(e, iInA, jNow, iInB, FSR);
0310 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
0311 pTnow = pT(e, iInB, jNow, iInA, FSR);
0312 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
0313
0314
0315
0316 } else {
0317 for (int kNow = 0; kNow < e.size(); kNow++) {
0318 if (kNow == jNow || !e[kNow].isFinal()) continue;
0319 if (QEDvetoMode==0 && e[kNow].colType() == 0) continue;
0320
0321
0322
0323 pTnow = pT(e, kNow, jNow, iInA, FSR);
0324 if (pTnow > 0.) pTemt = (pTemt < 0)
0325 ? pTnow : min(pTemt, pTnow);
0326 pTnow = pT(e, kNow, jNow, iInB, FSR);
0327 if (pTnow > 0.) pTemt = (pTemt < 0)
0328 ? pTnow : min(pTemt, pTnow);
0329
0330
0331 for (int rNow = 0; rNow < e.size(); rNow++) {
0332 if (rNow == kNow || rNow == jNow ||
0333 !e[rNow].isFinal()) continue;
0334 if(QEDvetoMode==0 && e[rNow].colType() == 0) continue;
0335 pTnow = pT(e, kNow, jNow, rNow, FSR);
0336 if (pTnow > 0.) pTemt = (pTemt < 0)
0337 ? pTnow : min(pTemt, pTnow);
0338 }
0339
0340 }
0341
0342 }
0343
0344 }
0345
0346 }
0347
0348 }
0349 }
0350
0351
0352 return pTemt;
0353 }
0354
0355
0356
0357
0358
0359
0360
0361
0362 inline bool canVetoMPIStep() { return true; }
0363 inline int numberVetoMPIStep() { return 1; }
0364 inline bool doVetoMPIStep(int nMPI, const Event &e) {
0365
0366 if (nMPI > 1) return false;
0367 int iEmt = -1;
0368 double pT1(0.), pTsum(0.);
0369
0370
0371
0372
0373 if (nFinal > 0) {
0374
0375
0376
0377 int count = 0;
0378 for (int i = e.size() - 1; i > 0; i--) {
0379 if (e[i].isFinal()) {
0380 count++;
0381 pT1 = e[i].pT();
0382 pTsum += e[i].pT();
0383 } else break;
0384 }
0385
0386 if (count != nFinal && count != nFinal + 1) {
0387 loggerPtr->ABORT_MSG("wrong number of final state particles in event");
0388 exit(1);
0389 }
0390
0391 isEmt = (count == nFinal) ? false : true;
0392 iEmt = (isEmt) ? e.size() - 1 : -1;
0393
0394
0395
0396
0397 } else {
0398
0399
0400 isEmt = false;
0401 for (int i = e.size() - 1; i > 0; i--) {
0402 if (e[i].isFinal()) {
0403 if ( e[i].isParton() && iEmt < 0
0404 && e[e[i].mother1()].isParton() ) {
0405 isEmt = true;
0406 iEmt = i;
0407 }
0408 pT1 = e[i].pT();
0409 pTsum += e[i].pT();
0410 } else break;
0411 }
0412 }
0413
0414
0415 if (!isEmt || pThardMode == 0) {
0416 pThard = infoPtr->scalup();
0417
0418
0419
0420 } else if (pThardMode == 1) {
0421 if (nFinal < 0) {
0422 loggerPtr->WARNING_MSG(
0423 "pThard == 1 not available for nFinal == -1, reset pThard = 0.");
0424 pThardMode = 0;
0425 pThard = infoPtr->scalup();
0426 } else {
0427 pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
0428 }
0429
0430
0431
0432
0433 } else if (pThardMode == 2) {
0434 if (nFinal < 0) {
0435 loggerPtr->WARNING_MSG(
0436 "pThard == 2 not available for nFinal == -1, reset pThard = 0.");
0437 pThardMode = 0;
0438 pThard = infoPtr->scalup();
0439 } else {
0440 pThard = pTcalc(e, -1, -1, -1, -1, -1);
0441 }
0442 }
0443
0444
0445 if (MPIvetoMode == 1) {
0446 pTMPI = (isEmt) ? pTsum / 2. : pT1;
0447 }
0448
0449
0450 accepted = false;
0451 nAcceptSeq = nISRveto = nFSRveto = 0;
0452
0453
0454 return false;
0455 }
0456
0457
0458
0459
0460
0461 inline bool canVetoISREmission() { return (vetoMode == 0) ? false : true; }
0462 inline bool doVetoISREmission(int, const Event &e, int iSys) {
0463
0464 if (iSys != 0) return false;
0465
0466
0467 if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
0468
0469
0470 int iRadAft = -1, iEmt = -1, iRecAft = -1;
0471 for (int i = e.size() - 1; i > 0; i--) {
0472 if (showerModel == 1) {
0473
0474 if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
0475 else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
0476 else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
0477 } else if (showerModel == 2) {
0478
0479 if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
0480 else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
0481 else if (iRecAft == -1
0482 && (e[i].status() == -41 || e[i].status() == 44)) iRecAft = i;
0483 }
0484 if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break;
0485 }
0486 if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
0487 loggerPtr->ABORT_MSG("could not find ISR emission");
0488 exit(1);
0489 }
0490
0491
0492
0493
0494 int xSR = (pTemtMode == 0) ? 0 : -1;
0495 int i = (pTemtMode == 0) ? iRadAft : -1;
0496 int j = (pTemtMode != 2) ? iEmt : -1;
0497 int k = -1;
0498 int r = (pTemtMode == 0) ? iRecAft : -1;
0499 double pTemt = pTcalc(e, i, j, k, r, xSR);
0500
0501
0502
0503 bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
0504 ? false: true;
0505
0506
0507 if (pTemt > pThard) {
0508 if(!vetoParton) {
0509
0510 nAcceptSeq = vetoCount-1;
0511 } else {
0512 nAcceptSeq = 0;
0513 nISRveto++;
0514 return true;
0515 }
0516 }
0517
0518
0519 nAcceptSeq++;
0520 accepted = true;
0521 return false;
0522 }
0523
0524
0525
0526
0527
0528 inline bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; }
0529 inline bool doVetoFSREmission(int, const Event &e, int iSys, bool) {
0530
0531 if (iSys != 0) return false;
0532
0533
0534 if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
0535
0536
0537 int iRecAft = e.size() - 1;
0538 int iEmt = e.size() - 2;
0539 int iRadAft = e.size() - 3;
0540 int iRadBef = e[iEmt].mother1();
0541 bool stop = false;
0542 if (showerModel == 2) {
0543
0544 if ( (e[iRecAft].status() != 51 && e[iRecAft].status() != 52) ||
0545 e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop = true;
0546 } else {
0547
0548 if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
0549 e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop = true;
0550 }
0551 if (stop) {
0552 loggerPtr->ABORT_MSG("could not find FSR emission");
0553 exit(1);
0554 }
0555
0556
0557
0558
0559
0560 int xSR = (pTemtMode == 0) ? 1 : -1;
0561 int i = (pTemtMode == 0) ? iRadBef : -1;
0562 int k = (pTemtMode == 0) ? iRadAft : -1;
0563 int r = (pTemtMode == 0) ? iRecAft : -1;
0564
0565
0566 double pTemt = 0.;
0567 if (pTemtMode == 0 || pTemtMode == 1) {
0568
0569
0570
0571
0572
0573 int j = iRadAft;
0574 if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
0575
0576 for (int jLoop = 0; jLoop < 2; jLoop++) {
0577 if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
0578 else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
0579
0580
0581 if (emittedMode != 3) break;
0582 if (k != -1) swap(j, k); else j = iEmt;
0583 }
0584
0585
0586 } else if (pTemtMode == 2) {
0587 pTemt = pTcalc(e, i, -1, k, r, xSR);
0588
0589 }
0590
0591
0592
0593 bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
0594 ? false: true;
0595
0596
0597 if (pTemt > pThard) {
0598 if(!vetoParton) {
0599
0600 nAcceptSeq = vetoCount-1;
0601 } else {
0602 nAcceptSeq = 0;
0603 nFSRveto++;
0604 return true;
0605 }
0606 }
0607
0608
0609 nAcceptSeq++;
0610 accepted = true;
0611 return false;
0612 }
0613
0614
0615
0616
0617
0618 inline bool canVetoMPIEmission() {return (MPIvetoMode == 0) ? false : true;}
0619 inline bool doVetoMPIEmission(int, const Event &e) {
0620 if (MPIvetoMode == 1) {
0621 if (e[e.size() - 1].pT() > pTMPI) return true;
0622 }
0623 return false;
0624 }
0625
0626
0627
0628
0629
0630 inline int getNISRveto() { return nISRveto; }
0631 inline int getNFSRveto() { return nFSRveto; }
0632
0633
0634
0635 protected:
0636 double pThard, pTMPI;
0637
0638 private:
0639 int showerModel, nFinal, vetoMode, MPIvetoMode, QEDvetoMode, vetoCount;
0640 int pThardMode, pTemtMode, emittedMode, pTdefMode;
0641 bool accepted, isEmt;
0642
0643
0644 int nAcceptSeq;
0645
0646 unsigned long int nISRveto, nFSRveto;
0647
0648 };
0649
0650
0651
0652
0653
0654 PYTHIA8_PLUGIN_CLASS(UserHooks, PowhegHooks, false, false, false)
0655 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
0656
0657
0658
0659 }
0660
0661 #endif