1 // main31.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2013 Richard Corke, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
7 using namespace Pythia8;
9 //==========================================================================
11 // Use userhooks to veto PYTHIA emissions above the POWHEG scale.
13 class PowhegHooks : public UserHooks {
17 // Constructor and destructor.
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) {};
27 //--------------------------------------------------------------------------
29 // Routines to calculate the pT (according to pTdefMode) in a splitting:
30 // ISR: i (radiator after) -> j (emitted after) k (radiator before)
31 // FSR: i (radiator before) -> j (emitted after) k (radiator after)
32 // For the Pythia pT definition, a recoiler (after) must be specified.
34 // Compute the Pythia pT separation. Based on pTLund function in History.cc
35 double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch,
36 int RecAfterBranch, bool FSR) {
38 // Convenient shorthands for later
39 Vec4 radVec = e[RadAfterBranch].p();
40 Vec4 emtVec = e[EmtAfterBranch].p();
41 Vec4 recVec = e[RecAfterBranch].p();
42 int radID = e[RadAfterBranch].id();
44 // Calculate virtuality of splitting
45 double sign = (FSR) ? 1. : -1.;
46 Vec4 Q(radVec + sign * emtVec);
47 double Qsq = sign * Q.m2Calc();
49 // Mass term of radiator
50 double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
51 pow2(particleDataPtr->m0(radID)) : 0.;
53 // z values for FSR and ISR
56 // Construct 2 -> 3 variables
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;
65 // Construct dipoles before/after splitting
66 Vec4 qBR(radVec - emtVec + recVec);
67 Vec4 qAR(radVec + recVec);
68 z = qBR.m2Calc() / qAR.m2Calc();
72 // Virtuality with correct sign
73 pTnow *= (Qsq - sign * m2Rad);
75 // Can get negative pT for massive splittings
77 cout << "Warning: pTpythia was negative" << endl;
82 cout << "pTpythia: rad = " << RadAfterBranch << ", emt = "
83 << EmtAfterBranch << ", rec = " << RecAfterBranch
84 << ", pTnow = " << sqrt(pTnow) << endl;
91 // Compute the POWHEG pT separation between i and j
92 double pTpowheg(const Event &e, int i, int j, bool FSR) {
94 // pT value for FSR and ISR
97 // POWHEG d_ij (in CM frame). Note that the incoming beams have not
98 // been updated in the parton systems pointer yet (i.e. prior to any
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()) );
112 // POWHEG pT_ISR is just kinematic pT
118 cout << "Warning: pTpowheg was negative" << endl;
123 cout << "pTpowheg: i = " << i << ", j = " << j
124 << ", pTnow = " << pTnow << endl;
130 // Calculate pT for a splitting based on pTdefMode.
131 // If j is -1, all final-state partons are tried.
132 // If i, k, r and xSR are -1, then all incoming and outgoing
133 // partons are tried.
134 // xSR set to 0 means ISR, while xSR set to 1 means FSR
135 double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) {
137 // Loop over ISR and FSR if necessary
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;
145 // If all necessary arguments have been given, then directly calculate.
146 // POWHEG ISR and FSR, need i and j.
147 if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
148 pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
150 // Pythia ISR, need i, j and r.
151 } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
152 pTemt = pTpythia(e, i, j, r, FSR);
154 // Pythia FSR, need k, j and r.
155 } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
156 pTemt = pTpythia(e, k, j, r, FSR);
158 // Otherwise need to try all possible combintations.
160 // Start by finding incoming legs to the hard system after
161 // branching (radiator after branching, i for ISR).
162 // Use partonSystemsPtr to find incoming just prior to the
163 // branching and track mothers.
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(); }
169 // If we do not have j, then try all final-state partons
170 int jNow = (j > 0) ? j : 0;
171 int jMax = (j > 0) ? j + 1 : e.size();
172 for (; jNow < jMax; jNow++) {
174 // Final-state and coloured jNow only
175 if (!e[jNow].isFinal() || e[jNow].colType() == 0) continue;
178 if (pTdefMode == 0 || pTdefMode == 1) {
180 // ISR - only done once as just kinematical pT
182 pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
183 if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
185 // FSR - try all outgoing partons from system before branching
186 // as i. Note that for the hard system, there is no
187 // "before branching" information.
190 int outSize = partonSystemsPtr->sizeOut(0);
191 for (int iMem = 0; iMem < outSize; iMem++) {
192 int iNow = partonSystemsPtr->getOut(0, iMem);
194 // Coloured only, i != jNow and no carbon copies
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) {
210 // ISR - other incoming as recoiler
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);
217 // FSR - try all final-state coloured partons as radiator
218 // after emission (k).
220 for (int kNow = 0; kNow < e.size(); kNow++) {
221 if (kNow == jNow || !e[kNow].isFinal() ||
222 e[kNow].colType() == 0) continue;
224 // For this kNow, need to have a recoiler.
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);
233 // Try all other outgoing.
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;
258 //--------------------------------------------------------------------------
260 // Extraction of pThard based on the incoming event.
261 // Assume that all the final-state particles are in a continuous block
262 // at the end of the event and the final entry is the POWHEG emission.
263 // If there is no POWHEG emission, then pThard is set to Qfac.
265 bool canVetoMPIStep() { return true; }
266 int numberVetoMPIStep() { return 1; }
267 bool doVetoMPIStep(int nMPI, const Event &e) {
268 // Extra check on nMPI
269 if (nMPI > 1) return false;
271 // Find if there is a POWHEG emission. Go backwards through the
272 // event record until there is a non-final particle. Also sum pT and
273 // find pT_1 for possible MPI vetoing
275 double pT1 = 0., pTsum = 0.;
276 for (int i = e.size() - 1; i > 0; i--) {
277 if (e[i].isFinal()) {
283 // Extra check that we have the correct final state
284 if (count != nFinal && count != nFinal + 1) {
285 cout << "Error: wrong number of final state particles in event" << endl;
288 // Flag if POWHEG radiation present and index
289 bool isEmt = (count == nFinal) ? false : true;
290 int iEmt = (isEmt) ? e.size() - 1 : -1;
292 // If there is no radiation or if pThardMode is 0 then set pThard to Qfac.
293 if (!isEmt || pThardMode == 0) {
294 pThard = infoPtr->QFac();
296 // If pThardMode is 1 then the pT of the POWHEG emission is checked against
297 // all other incoming and outgoing partons, with the minimal value taken
298 } else if (pThardMode == 1) {
299 pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
301 // If pThardMode is 2, then the pT of all final-state partons is checked
302 // against all other incoming and outgoing partons, with the minimal value
304 } else if (pThardMode == 2) {
305 pThard = pTcalc(e, -1, -1, -1, -1, -1);
309 // Find MPI veto pT if necessary
310 if (MPIvetoMode == 1) {
311 pTMPI = (isEmt) ? pTsum / 2. : pT1;
315 cout << "doVetoMPIStep: Qfac = " << infoPtr->QFac()
316 << ", pThard = " << pThard << endl << endl;
319 // Initialise other variables
321 nAcceptSeq = nISRveto = nFSRveto = 0;
323 // Do not veto the event
327 //--------------------------------------------------------------------------
331 bool canVetoISREmission() { return (vetoMode == 0) ? false : true; }
332 bool doVetoISREmission(int, const Event &e, int iSys) {
333 // Must be radiation from the hard system
334 if (iSys != 0) return false;
336 // If we already have accepted 'vetoCount' emissions in a row, do nothing
337 if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
339 // Pythia radiator after, emitted and recoiler after.
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;
353 // pTemtMode == 0: pT of emitted w.r.t. radiator
354 // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
355 // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
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;
367 // Veto if pTemt > pThard
368 if (pTemt > pThard) {
374 // Else mark that an emission has been accepted and continue
380 //--------------------------------------------------------------------------
384 bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; }
385 bool doVetoFSREmission(int, const Event &e, int iSys, bool) {
386 // Must be radiation from the hard system
387 if (iSys != 0) return false;
389 // If we already have accepted 'vetoCount' emissions in a row, do nothing
390 if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
392 // Pythia radiator (before and after), emitted and recoiler (after)
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;
404 // Behaviour based on pTemtMode:
405 // 0 - pT of emitted w.r.t. radiator before
406 // 1 - min(pT of emitted w.r.t. all incoming/outgoing)
407 // 2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
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;
413 // When pTemtMode is 0 or 1, iEmt has been selected
415 if (pTemtMode == 0 || pTemtMode == 1) {
416 // Which parton is emitted, based on emittedMode:
417 // 0 - Pythia definition of emitted
418 // 1 - Pythia definition of radiated after emission
419 // 2 - Random selection of emitted or radiated after emission
420 // 3 - Try both emitted and radiated after emission
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));
428 // For emittedMode == 3, have tried iRadAft, now try iEmt
429 if (emittedMode != 3) break;
430 if (k != -1) swap(j, k); else j = iEmt;
433 // If pTemtMode is 2, then try all final-state partons as emitted
434 } else if (pTemtMode == 2) {
435 pTemt = pTcalc(e, i, -1, k, r, xSR);
440 cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl;
443 // Veto if pTemt > pThard
444 if (pTemt > pThard) {
450 // Else mark that an emission has been accepted and continue
456 //--------------------------------------------------------------------------
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;
474 //--------------------------------------------------------------------------
476 // Functions to return information
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;
486 // The number of accepted emissions (in a row)
488 // Statistics on vetos
489 unsigned long int nISRveto, nFSRveto;
493 //==========================================================================
495 int main(int, char **) {
500 // Add further settings that can be set in the configuration file
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);
510 // Load configuration file
511 pythia.readFile("main31.cmnd");
513 // Read in main settings
514 int nEvent = pythia.settings.mode("Main:numberOfEvents");
515 int nError = pythia.settings.mode("Main:timesAllowErrors");
516 // Read in POWHEG settings
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);
527 // Add in user hooks for shower vetoing
528 PowhegHooks *powhegHooks = NULL;
531 // Set ISR and FSR to start at the kinematical limit
533 pythia.readString("SpaceShower:pTmaxMatch = 2");
534 pythia.readString("TimeShower:pTmaxMatch = 2");
537 // Set MPI to start at the kinematical limit
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);
547 // Initialise and list settings
550 // Counters for number of ISR/FSR emissions vetoed
551 unsigned long int nISRveto = 0, nFSRveto = 0;
553 // Begin event loop; generate until nEvent events are processed
554 // or end of LHEF file
555 int iEvent = 0, iError = 0;
558 // Generate the next event
559 if (!pythia.next()) {
561 // If failure because reached end of file then exit event loop
562 if (pythia.info.atEndOfFile()) break;
564 // Otherwise count event failure and continue/exit as necessary
565 cout << "Warning: event " << iEvent << " failed" << endl;
566 if (++iError == nError) {
567 cout << "Error: too many event failures.. exiting" << endl;
575 * Process dependent checks and analysis may be inserted here
578 // Update ISR/FSR veto counters
580 nISRveto += powhegHooks->getNISRveto();
581 nFSRveto += powhegHooks->getNFSRveto();
584 // If nEvent is set, check and exit loop if necessary
586 if (nEvent != 0 && iEvent == nEvent) break;
588 } // End of event loop.
590 // Statistics, histograms and veto information
592 cout << "Number of ISR emissions vetoed: " << nISRveto << endl;
593 cout << "Number of FSR emissions vetoed: " << nFSRveto << endl;
597 if (powhegHooks) delete powhegHooks;