PYTHIA8/pythia8145/examples/main71.cc

   1 // main71.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 // Test program that illustrates how a Les Houches Event File
   7 // written by POWHEG-hvq (top pair production) can processed with
   8 // PYTHIA. For both initial- and final-state radiation, the parton
   9 // showers are started at the kinematical limit, and emissions
  10 // above the POWHEG scale are vetoed.
  11 // For a detailed physics discussion see R. Corke and T. Sjostrand,
  12 // LU-TP-10-07, MCNET-10-04, arXiv:1003.2384 [hep-ph].
  13
  14 #include "Pythia.h"
  15 using namespace Pythia8;
  16
  17 //==========================================================================
  18
  19 // pT of emissions and vetoing if necessary.
  20
  21 class MyUserHooks : public UserHooks {
  22
  23 public:
  24
  25    // Constructor and destructor.
  26    MyUserHooks() : firstNoRad(true) { }
  27   ~MyUserHooks() { }
  28
  29   // Use VetoMIStep to analyse the incoming LHEF event and
  30   // extract the veto scale
  31   bool canVetoMIStep()    { return true; }
  32   int  numberVetoMIStep() { return 1; }
  33   bool doVetoMIStep(int, const Event &e) {
  34     // Check that partons 5 and 6 are the t/tbar pair
  35     if (e[5].id() != 6 || e[6].id() != -6) {
  36       cout << "Error: could not find t/tbar pair" << endl;
  37       e.list();
  38       exit(1);
  39     }
  40
  41     // Some events may not have radiation from POWHEG
  42     switch (e.size()) {
  43     case 7:
  44       // No radiation - veto scale is given by the factorisation scale
  45       pTpowheg = -1.;
  46       pTveto   = infoPtr->QFac();
  47       noRad    = true;
  48
  49       // If this is the first no radiation event, then print scale
  50       if (firstNoRad) {
  51         cout << "Info: no POWHEG radiation, Q = " << pTveto
  52              << " GeV" << endl;
  53         firstNoRad = false;
  54       }
  55       break;
  56
  57     case 8:
  58       // Radiation is parton 7 - first check that it is as expected
  59       if (e[7].id() != 21 && e[7].idAbs() > 5) {
  60         cout << "Error: jet is not quark/gluon" << endl;
  61         e.list();
  62         exit(1);
  63       }
  64       // Veto scale is given by jet pT
  65       pTveto = pTpowheg = e[7].pT();
  66       noRad  = false;
  67       break;
  68     }
  69
  70     // Initialise other variables
  71     nISRveto = nFSRveto = 0;
  72     pTshower = -1.;
  73
  74     // Do not veto the event
  75     return false;
  76   }
  77
  78   // For subsequent ISR/FSR emissions, find the pT of the shower
  79   // emission and veto as necessary
  80   bool canVetoISREmission() { return true; }
  81   bool doVetoISREmission(int, const Event &e) {
  82     // ISR - next shower emission is given status 43
  83     int i;
  84     for (i = e.size() - 1; i > 6; i--)
  85       if (e[i].isFinal() && e[i].status() == 43) break;
  86     if (i == 6) {
  87       cout << "Error: couldn't find ISR emission" << endl;
  88       e.list();
  89       exit(1);
  90     }
  91
  92     // Veto if above the POWHEG scale
  93     if (e[i].pT() > pTveto) {
  94       nISRveto++;
  95       return true;
  96     }
  97     // Store the first shower emission pT
  98     if (pTshower < 0.) pTshower = e[i].pT();
  99
 100     return false;
 101   }
 102
 103   bool canVetoFSREmission() { return true; }
 104   bool doVetoFSREmission(int, const Event &e) {
 105     // FSR - shower emission will have status 51 and not be t/tbar
 106     int i;
 107     for (i = e.size() - 1; i > 6; i--)
 108       if (e[i].isFinal() && e[i].status() == 51 &&
 109           e[i].idAbs() != 6) break;
 110     if (i == 6) {
 111       cout << "Error: couldn't find FSR emission" << endl;
 112       e.list();
 113       exit(1);
 114     }
 115
 116     // Veto if above the POWHEG scale
 117     if (e[i].pT() > pTveto) {
 118       nFSRveto++;
 119       return true;
 120     }
 121     // Store the first shower emission pT
 122     if (pTshower < 0.) pTshower = e[i].pT();
 123
 124     return false;
 125   }
 126
 127   // Functions to return information
 128   double getPTpowheg() { return pTpowheg; }
 129   double getPTshower() { return pTshower; }
 130   int    getNISRveto() { return nISRveto; }
 131   int    getNFSRveto() { return nFSRveto; }
 132   bool   getNoRad()    { return noRad;    }
 133
 134 private:
 135
 136   double pTveto, pTpowheg, pTshower;
 137   int    nISRveto, nFSRveto;
 138   bool   noRad, firstNoRad;
 139 };
 140
 141 //==========================================================================
 142
 143 int main(int, char **) {
 144
 145   // Generator
 146   Pythia pythia;
 147
 148   // Set the starting shower scales to the kinematic limit
 149   pythia.readString("SpaceShower:pTmaxMatch = 2");
 150   pythia.readString("TimeShower:pTmaxMatch  = 2");
 151
 152   // Turn off MPI, hadronisation and top decays
 153   // To be commented-out for complete run!
 154   pythia.readString("PartonLevel:MI = off");
 155   pythia.readString("HadronLevel:All = off");
 156   pythia.particleData.readString("6:mayDecay = off");
 157
 158   // Add in user hooks for shower vetoing
 159   MyUserHooks *myUserHooks = new MyUserHooks();
 160   UserHooks* pythiaUserHooks = myUserHooks;
 161   pythia.setUserHooksPtr(pythiaUserHooks);
 162
 163   // Initialize Les Houches Event File run.
 164   pythia.init("powheg-hvq.lhe");
 165   // List initialization information.
 166   pythia.settings.listChanged();
 167
 168   // Histograms
 169   Hist pTpowhegH("pT of POWHEG emission", 100, 0., 100.);
 170   Hist pTshowerH("pT of first shower emission", 100, 0., 100.);
 171
 172   // Counters for events with no POWHEG radiation and number of
 173   // ISR/FSR emissions vetoed
 174   int nNoRad = 0, nISRveto = 0, nFSRveto = 0;
 175
 176   // Begin event loop; generate until none left in input file
 177   int iEvent = 0;
 178   while (true) {
 179     if (iEvent % 10 == 0)
 180       cout << "Begin event " << iEvent << endl;
 181
 182     if (!pythia.next()) {
 183       // If failure because reached end of file then exit event loop
 184       if (pythia.info.atEndOfFile()) {
 185         cout << "Info: end of Les Houches file reached" << endl;
 186         break;
 187       }
 188       cout << "Warning: event " << iEvent << "failed" << endl;
 189       continue;
 190     }
 191
 192     // List first event: Les Houches, hard process and complete
 193     if (iEvent < 1) {
 194       pythia.LHAeventList();
 195       pythia.info.list();
 196       pythia.process.list();
 197       pythia.event.list();
 198     }
 199
 200     // Fill in histograms with information from the UserHooks class
 201     double pTpowheg = myUserHooks->getPTpowheg();
 202     if (pTpowheg > 0.)
 203       pTpowhegH.fill(pTpowheg);
 204
 205     double pTshower = myUserHooks->getPTshower();
 206     if (pTshower > 0.)
 207       pTshowerH.fill(pTshower);
 208
 209     // Other information from the UserHooks class
 210     if (myUserHooks->getNoRad() == true) nNoRad++;
 211     nISRveto += myUserHooks->getNISRveto();
 212     nFSRveto += myUserHooks->getNFSRveto();
 213
 214     ++iEvent;
 215   } // End of event loop.
 216
 217   // Statistics, histograms and veto information
 218   pTpowhegH /= double(iEvent);
 219   pTshowerH /= double(iEvent);
 220   pythia.statistics();
 221   cout << pTpowhegH << pTshowerH << endl;
 222   cout << "Number of events with no POWHEG radiation: "
 223        << nNoRad << endl;
 224   cout << "Number of ISR emissions vetoed: " << nISRveto << endl;
 225   cout << "Number of FSR emissions vetoed: " << nFSRveto << endl;
 226   cout << endl;
 227
 228   // Done.
 229   return 0;
 230 }