PYTHIA8/pythia8175/examples/main31.cc

   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.
   5
   6 #include "Pythia.h"
   7 using namespace Pythia8;
   8
   9 //==========================================================================
  10
  11 // Use userhooks to veto PYTHIA emissions above the POWHEG scale.
  12
  13 class PowhegHooks : public UserHooks {
  14
  15 public:
  16
  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) {};
  25   ~PowhegHooks() {}
  26
  27 //--------------------------------------------------------------------------
  28
  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.
  33
  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) {
  37
  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();
  43
  44     // Calculate virtuality of splitting
  45     double sign = (FSR) ? 1. : -1.;
  46     Vec4 Q(radVec + sign * emtVec);
  47     double Qsq = sign * Q.m2Calc();
  48
  49     // Mass term of radiator
  50     double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
  51                    pow2(particleDataPtr->m0(radID)) : 0.;
  52
  53     // z values for FSR and ISR
  54     double z, pTnow;
  55     if (FSR) {
  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;
  61       z     = x1 / (x1 + x3);
  62       pTnow = z * (1. - z);
  63
  64     } else {
  65       // Construct dipoles before/after splitting
  66       Vec4 qBR(radVec - emtVec + recVec);
  67       Vec4 qAR(radVec + recVec);
  68       z     = qBR.m2Calc() / qAR.m2Calc();
  69       pTnow = (1. - z);
  70     }
  71
  72     // Virtuality with correct sign
  73     pTnow *= (Qsq - sign * m2Rad);
  74
  75     // Can get negative pT for massive splittings
  76     if (pTnow < 0.) {
  77       cout << "Warning: pTpythia was negative" << endl;
  78       return -1.;
  79     }
  80
  81 #ifdef DBGOUTPUT
  82     cout << "pTpythia: rad = " << RadAfterBranch << ", emt = "
  83          << EmtAfterBranch << ", rec = " << RecAfterBranch
  84          << ", pTnow = " << sqrt(pTnow) << endl;
  85 #endif
  86
  87     // Return pT
  88     return sqrt(pTnow);
  89   }
  90
  91   // Compute the POWHEG pT separation between i and j
  92   double pTpowheg(const Event &e, int i, int j, bool FSR) {
  93
  94     // pT value for FSR and ISR
  95     double pTnow = 0.;
  96     if (FSR) {
  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
  99       // potential recoil).
 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()) );
 110
 111     } else {
 112       // POWHEG pT_ISR is just kinematic pT
 113       pTnow = e[j].pT();
 114     }
 115
 116     // Check result
 117     if (pTnow < 0.) {
 118       cout << "Warning: pTpowheg was negative" << endl;
 119       return -1.;
 120     }
 121
 122 #ifdef DBGOUTPUT
 123     cout << "pTpowheg: i = " << i << ", j = " << j
 124          << ", pTnow = " << pTnow << endl;
 125 #endif
 126
 127     return pTnow;
 128   }
 129
 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) {
 136
 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++) {
 142       // FSR flag
 143       bool FSR = (xSR == 0) ? false : true;
 144
 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);
 149
 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);
 153
 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);
 157
 158       // Otherwise need to try all possible combintations.
 159       } else {
 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(); }
 168
 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++) {
 173
 174           // Final-state and coloured jNow only
 175           if (!e[jNow].isFinal() || e[jNow].colType() == 0) continue;
 176
 177           // POWHEG
 178           if (pTdefMode == 0 || pTdefMode == 1) {
 179
 180             // ISR - only done once as just kinematical pT
 181             if (!FSR) {
 182               pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
 183               if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
 184
 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.
 188               } else {
 189
 190                 int outSize = partonSystemsPtr->sizeOut(0);
 191                 for (int iMem = 0; iMem < outSize; iMem++) {
 192                   int iNow = partonSystemsPtr->getOut(0, iMem);
 193
 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;
 198
 199                   pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0)
 200                     ? false : FSR);
 201                   if (pTnow > 0.) pTemt = (pTemt < 0)
 202                     ? pTnow : min(pTemt, pTnow);
 203                 } // for (iMem)
 204
 205               } // if (!FSR)
 206
 207           // Pythia
 208           } else if (pTdefMode == 2) {
 209
 210             // ISR - other incoming as recoiler
 211             if (!FSR) {
 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);
 216
 217             // FSR - try all final-state coloured partons as radiator
 218             //       after emission (k).
 219             } else {
 220               for (int kNow = 0; kNow < e.size(); kNow++) {
 221                 if (kNow == jNow || !e[kNow].isFinal() ||
 222                     e[kNow].colType() == 0) continue;
 223
 224                 // For this kNow, need to have a recoiler.
 225                 // Try two incoming.
 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);
 232
 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);
 240                 } // for (rNow)
 241
 242               } // for (kNow)
 243             } // if (!FSR)
 244           } // if (pTdefMode)
 245         } // for (j)
 246       }
 247     } // for (xSR)
 248
 249 #ifdef DBGOUTPUT
 250     cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k
 251          << ", r = " << r << ", xSR = " << xSRin
 252          << ", pTemt = " << pTemt << endl;
 253 #endif
 254
 255     return pTemt;
 256   }
 257
 258 //--------------------------------------------------------------------------
 259
 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.
 264
 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;
 270
 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
 274     int    count = 0;
 275     double pT1 = 0., pTsum = 0.;
 276     for (int i = e.size() - 1; i > 0; i--) {
 277       if (e[i].isFinal()) {
 278         count++;
 279         pT1    = e[i].pT();
 280         pTsum += e[i].pT();
 281       } else break;
 282     }
 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;
 286       exit(1);
 287     }
 288     // Flag if POWHEG radiation present and index
 289     bool isEmt = (count == nFinal) ? false : true;
 290     int  iEmt  = (isEmt) ? e.size() - 1 : -1;
 291
 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();
 295
 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);
 300
 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
 303     // taken
 304     } else if (pThardMode == 2) {
 305       pThard = pTcalc(e, -1, -1, -1, -1, -1);
 306
 307     }
 308
 309     // Find MPI veto pT if necessary
 310     if (MPIvetoMode == 1) {
 311       pTMPI = (isEmt) ? pTsum / 2. : pT1;
 312     }
 313
 314 #ifdef DBGOUTPUT
 315     cout << "doVetoMPIStep: Qfac = " << infoPtr->QFac()
 316          << ", pThard = " << pThard << endl << endl;
 317 #endif
 318
 319     // Initialise other variables
 320     accepted   = false;
 321     nAcceptSeq = nISRveto = nFSRveto = 0;
 322
 323     // Do not veto the event
 324     return false;
 325   }
 326
 327 //--------------------------------------------------------------------------
 328
 329   // ISR veto
 330
 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;
 335
 336     // If we already have accepted 'vetoCount' emissions in a row, do nothing
 337     if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
 338
 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;
 346     }
 347     if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
 348       e.list();
 349       cout << "Error: couldn't find Pythia ISR emission" << endl;
 350       exit(1);
 351     }
 352
 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;
 359     int k        = -1;
 360     int r        = (pTemtMode == 0) ? iRecAft : -1;
 361     double pTemt = pTcalc(e, i, j, k, r, xSR);
 362
 363 #ifdef DBGOUTPUT
 364     cout << "doVetoISREmission: pTemt = " << pTemt << endl << endl;
 365 #endif
 366
 367     // Veto if pTemt > pThard
 368     if (pTemt > pThard) {
 369       nAcceptSeq = 0;
 370       nISRveto++;
 371       return true;
 372     }
 373
 374     // Else mark that an emission has been accepted and continue
 375     nAcceptSeq++;
 376     accepted = true;
 377     return false;
 378   }
 379
 380 //--------------------------------------------------------------------------
 381
 382   // FSR veto
 383
 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;
 388
 389     // If we already have accepted 'vetoCount' emissions in a row, do nothing
 390     if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
 391
 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) {
 399       e.list();
 400       cout << "Error: couldn't find Pythia FSR emission" << endl;
 401       exit(1);
 402     }
 403
 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;
 412
 413     // When pTemtMode is 0 or 1, iEmt has been selected
 414     double pTemt = 0.;
 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
 421       int j = iRadAft;
 422       if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
 423
 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));
 427
 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;
 431       }
 432
 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);
 436
 437     }
 438
 439 #ifdef DBGOUTPUT
 440     cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl;
 441 #endif
 442
 443     // Veto if pTemt > pThard
 444     if (pTemt > pThard) {
 445       nAcceptSeq = 0;
 446       nFSRveto++;
 447       return true;
 448     }
 449
 450     // Else mark that an emission has been accepted and continue
 451     nAcceptSeq++;
 452     accepted = true;
 453     return false;
 454   }
 455
 456 //--------------------------------------------------------------------------
 457
 458   // MPI veto
 459
 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) {
 464 #ifdef DBGOUTPUT
 465         cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT()
 466              << ", pTMPI = " << pTMPI << endl << endl;
 467 #endif
 468         return true;
 469       }
 470     }
 471     return false;
 472   }
 473
 474 //--------------------------------------------------------------------------
 475
 476   // Functions to return information
 477
 478   int    getNISRveto() { return nISRveto; }
 479   int    getNFSRveto() { return nFSRveto; }
 480
 481 private:
 482   int    nFinal, vetoMode, vetoCount, pThardMode, pTemtMode,
 483          emittedMode, pTdefMode, MPIvetoMode;
 484   double pThard, pTMPI;
 485   bool   accepted;
 486   // The number of accepted emissions (in a row)
 487   int nAcceptSeq;
 488   // Statistics on vetos
 489   unsigned long int nISRveto, nFSRveto;
 490
 491 };
 492
 493 //==========================================================================
 494
 495 int main(int, char **) {
 496
 497   // Generator
 498   Pythia pythia;
 499
 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);
 509
 510   // Load configuration file
 511   pythia.readFile("main31.cmnd");
 512
 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);
 526
 527   // Add in user hooks for shower vetoing
 528   PowhegHooks *powhegHooks = NULL;
 529   if (loadHooks) {
 530
 531     // Set ISR and FSR to start at the kinematical limit
 532     if (vetoMode > 0) {
 533       pythia.readString("SpaceShower:pTmaxMatch = 2");
 534       pythia.readString("TimeShower:pTmaxMatch = 2");
 535     }
 536
 537     // Set MPI to start at the kinematical limit
 538     if (MPIvetoMode > 0) {
 539       pythia.readString("MultipartonInteractions:pTmaxMatch = 2");
 540     }
 541
 542     powhegHooks = new PowhegHooks(nFinal, vetoMode, vetoCount,
 543         pThardMode, pTemtMode, emittedMode, pTdefMode, MPIvetoMode);
 544     pythia.setUserHooksPtr((UserHooks *) powhegHooks);
 545   }
 546
 547   // Initialise and list settings
 548   pythia.init();
 549
 550   // Counters for number of ISR/FSR emissions vetoed
 551   unsigned long int nISRveto = 0, nFSRveto = 0;
 552
 553   // Begin event loop; generate until nEvent events are processed
 554   // or end of LHEF file
 555   int iEvent = 0, iError = 0;
 556   while (true) {
 557
 558     // Generate the next event
 559     if (!pythia.next()) {
 560
 561       // If failure because reached end of file then exit event loop
 562       if (pythia.info.atEndOfFile()) break;
 563
 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;
 568         break;
 569       }
 570
 571       continue;
 572     }
 573
 574     /*
 575      * Process dependent checks and analysis may be inserted here
 576      */
 577
 578     // Update ISR/FSR veto counters
 579     if (loadHooks) {
 580       nISRveto += powhegHooks->getNISRveto();
 581       nFSRveto += powhegHooks->getNFSRveto();
 582     }
 583
 584     // If nEvent is set, check and exit loop if necessary
 585     ++iEvent;
 586     if (nEvent != 0 && iEvent == nEvent) break;
 587
 588   } // End of event loop.
 589
 590   // Statistics, histograms and veto information
 591   pythia.stat();
 592   cout << "Number of ISR emissions vetoed: " << nISRveto << endl;
 593   cout << "Number of FSR emissions vetoed: " << nFSRveto << endl;
 594   cout << endl;
 595
 596   // Done.
 597   if (powhegHooks) delete powhegHooks;
 598   return 0;
 599 }
 600