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