PYTHIA8/pythia8145/examples/main61.cc

   1 // main61.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 Richard Corke.
   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 /*
   7  * Simple example of fastjet analysis. Roughly follows analysis of:
   8  * T. Aaltonen et al. [CDF Collaboration],
   9  * Measurement of the cross section for W-boson production in association
  10  * with jets in ppbar collisions at sqrt(s)=1.96$ TeV
  11  * Phys. Rev. D 77 (2008) 011108
  12  * arXiv:0711.4044 [hep-ex]
  13  *
  14  * Cuts:
  15  *   ET(elec)     > 20GeV
  16  *   |eta(elec)|  < 1.1
  17  *   ET(missing)  > 30GeV
  18  *   ET(jet)      > 20GeV
  19  *   |eta(jet)|   < 2.0
  20  *   deltaR(elec, jet) > 0.52
  21  * Not used:
  22  *   mT(W)        > 20GeV
  23  */
  24
  25 #include "Pythia.h"
  26
  27 #include "fastjet/PseudoJet.hh"
  28 #include "fastjet/ClusterSequence.hh"
  29
  30 using namespace Pythia8;
  31
  32 // Experimental cross section
  33 // sigma(W -> ev + >= n-jet; ET(n'th-jet) > 25GeV), n = 0, 1, 2, 3, 4
  34 const double expCrossSec[] = { 798.0, 53.5, 6.8, 0.84, 0.074 };
  35
  36 int main() {
  37   // Settings
  38   int  nEvent = 10000;
  39   bool doMI   = true;
  40
  41   // Generator
  42   Pythia pythia;
  43
  44   // Single W production
  45   pythia.readString("WeakSingleBoson:ffbar2W = on");
  46   // Force decay W->ev
  47   pythia.readString("24:onMode = off");
  48   pythia.readString("24:onIfAny = 11 12");
  49   // Multiple Interactions
  50   if (doMI == false)
  51     pythia.readString("PartonLevel:MI = off");
  52
  53   // Initialisation, p pbar @ 1.96 TeV
  54   pythia.init( 2212, -2212, 1960.);
  55
  56   // Histograms
  57   Hist dSigma1("1-jet cross-section (E_jet1 > 20 GeV)", 70, 0.0, 350.0);
  58   Hist dSigma2("2-jet cross-section (E_jet2 > 20 GeV)", 38, 0.0, 190.0);
  59   Hist dSigma3("3-jet cross-section (E_jet3 > 20 GeV)", 16, 0.0, 80.0);
  60   Hist dSigma4("4-jet cross-section (E_jet4 > 20 GeV)",  7, 0.0, 35.0);
  61   Hist *dSigmaHist[5] = { NULL, &dSigma1, &dSigma2, &dSigma3, &dSigma4 };
  62   double dSigmaBin[5] = { 0.0, 350.0 / 70.0, 190.0 / 38.0,
  63                           80.0 / 16.0, 35.0 / 7.0 };
  64
  65   // Fastjet analysis - select algorithm and parameters
  66   double Rparam = 0.4;
  67   fastjet::Strategy               strategy = fastjet::Best;
  68   fastjet::RecombinationScheme    recombScheme = fastjet::E_scheme;
  69   fastjet::JetDefinition         *jetDef = NULL;
  70   jetDef = new fastjet::JetDefinition(fastjet::kt_algorithm, Rparam,
  71                                       recombScheme, strategy);
  72
  73   // Fastjet input
  74   std::vector <fastjet::PseudoJet> fjInputs;
  75
  76   // Statistics for later
  77   int nEventAccept25[5] = { 0, 0, 0, 0, 0 };
  78   int vetoCount[4] = { 0, 0, 0, 0 };
  79   const char *vetoStr[] = { "ET(elec)", "|eta(elec)|",
  80                             "ET(missing)", "deltaR(elec, jet)" };
  81   bool firstEvent = true;
  82
  83   // Begin event loop. Generate event. Skip if error. List first one.
  84   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
  85     if (iEvent % 10000 == 0) printf("Beginning event %d\n", iEvent);
  86     if (!pythia.next()) continue;
  87     if (iEvent < 1) {
  88       pythia.info.list();
  89       pythia.event.list();
  90     }
  91
  92     // Need to find the electron from the W decay - cheat a bit here
  93     // and find it from the W in the event record
  94     int idxW = -1;
  95     for (int i = pythia.event.size() - 1; i > 0; i--) {
  96       if (pythia.event[i].idAbs() == 24) {
  97         idxW = i;
  98         break;
  99       }
 100     }
 101     if (idxW == -1) {
 102       cout << "Error: Could not find W" << endl;
 103       continue;
 104     }
 105
 106     // Find the electron from the W decay
 107     int idxElec = idxW;
 108     while(true) {
 109       int daughter = pythia.event[idxElec].daughter1();
 110       if   (daughter == 0) break;
 111       else                 idxElec = daughter;
 112     }
 113     if (pythia.event[idxElec].idAbs() != 11 ||
 114        !pythia.event[idxElec].isFinal()) {
 115       cout << "Error: Found incorrect decay product of the W" << endl;
 116       continue;
 117     }
 118
 119     // Electron cuts
 120     if (pythia.event[idxElec].pT() < 20.0) {
 121       vetoCount[0]++;
 122       continue;
 123     }
 124     if (abs(pythia.event[idxElec].eta()) > 1.1) {
 125       vetoCount[1]++;
 126       continue;
 127     }
 128
 129     // Reset Fastjet input
 130     fjInputs.resize(0);
 131
 132     // Keep track of missing ET
 133     Vec4 missingETvec;
 134
 135     // Loop over event record to decide what to pass to FastJet
 136     for (int i = 0; i < pythia.event.size(); ++i) {
 137       // Final state only
 138       if (!pythia.event[i].isFinal())        continue;
 139
 140       // No neutrinos
 141       if (pythia.event[i].idAbs() == 12 || pythia.event[i].idAbs() == 14 ||
 142           pythia.event[i].idAbs() == 16)     continue;
 143
 144       // Only |eta| < 3.6
 145       if (fabs(pythia.event[i].eta()) > 3.6) continue;
 146
 147       // Missing ET
 148       missingETvec += pythia.event[i].p();
 149
 150       // Do not include the electron from the W decay
 151       if (i == idxElec)                      continue;
 152
 153       // Store as input to Fastjet
 154       fjInputs.push_back( fastjet::PseudoJet (pythia.event[i].px(),
 155                                                 pythia.event[i].py(), pythia.event[i].pz(),
 156                                                 pythia.event[i].e() ) );
 157     }
 158
 159     if (fjInputs.size() == 0) {
 160       cout << "Error: event with no final state particles" << endl;
 161       continue;
 162     }
 163
 164     // Run Fastjet algorithm
 165     vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 166     fastjet::ClusterSequence clustSeq(fjInputs, *jetDef);
 167
 168     // For the first event, print the FastJet details
 169     if (firstEvent) {
 170       cout << "Ran " << jetDef->description() << endl;
 171       cout << "Strategy adopted by FastJet was "
 172            << clustSeq.strategy_string() << endl << endl;
 173       firstEvent = false;
 174     }
 175
 176     // Extract inclusive jets sorted by pT (note minimum pT of 20.0 GeV)
 177     inclusiveJets = clustSeq.inclusive_jets(20.0);
 178     sortedJets    = sorted_by_pt(inclusiveJets);
 179
 180     // Missing ET cut
 181     double missingET = missingETvec.pT();
 182     if (missingET < 30.0) {
 183       vetoCount[2]++;
 184       continue;
 185     }
 186
 187     // Keep track of jets with pT > 20/25 GeV
 188     int  jetCount20 = 0, jetCount25 = 0;
 189     // For the deltaR calculation below
 190     bool vetoEvent = false;
 191     fastjet::PseudoJet fjElec(pythia.event[idxElec].px(),
 192                               pythia.event[idxElec].py(),
 193                               pythia.event[idxElec].pz(),
 194                               pythia.event[idxElec].e());
 195
 196     for (unsigned int i = 0; i < sortedJets.size(); i++) {
 197       // Only count jets that have |eta| < 2.0
 198       if (fabs(sortedJets[i].rap()) > 2.0) continue;
 199       // Check distance between W decay electron and jets
 200       if (fjElec.squared_distance(sortedJets[i]) < 0.52 * 0.52)
 201         { vetoEvent = true; break; }
 202
 203       // Fill dSigma histograms and count jets with ET > 25.0
 204       if (sortedJets[i].perp() > 25.0)
 205         jetCount25++;
 206
 207       if (jetCount20 <= 3)
 208         dSigmaHist[++jetCount20]->fill(sortedJets[i].perp());
 209     }
 210     if (vetoEvent) { vetoCount[3]++; continue; }
 211
 212     if (jetCount25 > 4) jetCount25 = 4;
 213     for (int i = jetCount25; i >= 0; i--)
 214       nEventAccept25[i]++;
 215
 216   // End of event loop.
 217   }
 218
 219   // Statistics
 220   pythia.statistics();
 221
 222   // Output histograms
 223   double sigmapb = pythia.info.sigmaGen() * 1.0E9;
 224
 225   for (int i = 1; i <= 4; i++)
 226     (*dSigmaHist[i]) = ((*dSigmaHist[i]) * sigmapb) / nEvent / dSigmaBin[i];
 227   cout << dSigma1 << dSigma2 << dSigma3 << dSigma4 << endl;
 228
 229   // Output cross-sections
 230   cout << "Jet algorithm is kT" << endl;
 231   cout << "Multiple interactions are switched "
 232        << ( (doMI) ? "on" : "off" ) << endl;
 233   cout << endl << nEvent << " events generated. " << nEventAccept25[0]
 234        << " events passed cuts." << endl;
 235   cout << "Vetos:" << endl;
 236   for (int i = 0; i < 4; i++)
 237     cout << "  " << vetoStr[i] << " = " << vetoCount[i] << endl;
 238
 239   cout << endl << "Inclusive cross-sections (pb):" << endl;
 240   for (int i = 0; i < 5; i++) {
 241     cout << scientific << setprecision(3)
 242          << "  " << i << "-jet - Pythia = "
 243          << ((double) nEventAccept25[i] / (double) nEvent) * sigmapb;
 244     cout << ", Experimental = " << expCrossSec[i];
 245     if (i != 0) {
 246       cout << scientific << setprecision(3)
 247            << ", Pythia ratio to " << i - 1 << "-jet = "
 248            << ((double) nEventAccept25[i] / (double) nEventAccept25[i - 1]);
 249       cout << scientific << setprecision(3)
 250            << ", Experimental ratio to " << i - 1 << "-jet = "
 251            << expCrossSec[i] / expCrossSec[i - 1];
 252     }
 253     cout << endl;
 254   }
 255
 256   return 0;
 257 }
 258