StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main05.cc
1 // main05.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 // This is a simple test program.
7 // It studies jet production at the LHC, using SlowJet and CellJet.
8 // Note: the two finders are intended to construct approximately the same
9 // jet properties, but provides output in slightly different format,
10 // and have here not been optimized to show maximum possible agreement.
11 
12 #include "Pythia8/Pythia.h"
13 using namespace Pythia8;
14 
15 int main() {
16 
17  // Number of events, generated and listed ones.
18  int nEvent = 500;
19  int nListJets = 5;
20 
21  // Generator. LHC process and output selection. Initialization.
22  Pythia pythia;
23  pythia.readString("Beams:eCM = 14000.");
24  pythia.readString("HardQCD:all = on");
25  pythia.readString("PhaseSpace:pTHatMin = 200.");
26  pythia.readString("Next:numberShowInfo = 0");
27  pythia.readString("Next:numberShowProcess = 0");
28  pythia.readString("Next:numberShowEvent = 0");
29  pythia.init();
30 
31  // Common parameters for the two jet finders.
32  double etaMax = 4.;
33  double radius = 0.7;
34  double pTjetMin = 10.;
35  // Exclude neutrinos (and other invisible) from study.
36  int nSel = 2;
37  // Range and granularity of CellJet jet finder.
38  int nEta = 80;
39  int nPhi = 64;
40 
41  // Set up SlowJet jet finder, with anti-kT clustering
42  // and pion mass assumed for non-photons..
43  SlowJet slowJet( -1, radius, pTjetMin, etaMax, nSel, 1);
44 
45  // Set up CellJet jet finder.
46  CellJet cellJet( etaMax, nEta, nPhi, nSel);
47 
48  // Histograms. Note similarity in names, even when the two jet finders
49  // do not calculate identically the same property (pT vs. ET, y vs. eta).
50  Hist nJetsS("number of jets, SlowJet", 50, -0.5, 49.5);
51  Hist nJetsC("number of jets, CellJet", 50, -0.5, 49.5);
52  Hist nJetsD("number of jets, CellJet - SlowJet", 45, -22.5, 22.5);
53  Hist eTjetsS("pT for jets, SlowJet", 100, 0., 500.);
54  Hist eTjetsC("eT for jets, CellJet", 100, 0., 500.);
55  Hist etaJetsS("y for jets, SlowJet", 100, -5., 5.);
56  Hist etaJetsC("eta for jets, CellJet", 100, -5., 5.);
57  Hist phiJetsS("phi for jets, SlowJwt", 100, -M_PI, M_PI);
58  Hist phiJetsC("phi for jets, CellJet", 100, -M_PI, M_PI);
59  Hist distJetsS("R distance between jets, SlowJet", 100, 0., 10.);
60  Hist distJetsC("R distance between jets, CellJet", 100, 0., 10.);
61  Hist eTdiffS("pT difference, SlowJet", 100, -100., 400.);
62  Hist eTdiffC("eT difference, CellJet", 100, -100., 400.);
63 
64  // Begin event loop. Generate event. Skip if error.
65  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
66  if (!pythia.next()) continue;
67 
68  // Analyze Slowet jet properties. List first few.
69  slowJet. analyze( pythia.event );
70  if (iEvent < nListJets) slowJet.list();
71 
72  // Fill SlowJet inclusive jet distributions.
73  nJetsS.fill( slowJet.sizeJet() );
74  for (int i = 0; i < slowJet.sizeJet(); ++i) {
75  eTjetsS.fill( slowJet.pT(i) );
76  etaJetsS.fill( slowJet.y(i) );
77  phiJetsS.fill( slowJet.phi(i) );
78  }
79 
80  // Fill SlowJet distance between jets.
81  for (int i = 0; i < slowJet.sizeJet() - 1; ++i)
82  for (int j = i +1; j < slowJet.sizeJet(); ++j) {
83  double dEta = slowJet.y(i) - slowJet.y(j);
84  double dPhi = abs( slowJet.phi(i) - slowJet.phi(j) );
85  if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
86  double dR = sqrt( pow2(dEta) + pow2(dPhi) );
87  distJetsS.fill( dR );
88  }
89 
90  // Fill SlowJet pT-difference between jets (to check ordering of list).
91  for (int i = 1; i < slowJet.sizeJet(); ++i)
92  eTdiffS.fill( slowJet.pT(i-1) - slowJet.pT(i) );
93 
94  // Analyze CellJet jet properties. List first few.
95  cellJet. analyze( pythia.event, pTjetMin, radius );
96  if (iEvent < nListJets) cellJet.list();
97 
98  // Fill CellJet inclusive jet distributions.
99  nJetsC.fill( cellJet.size() );
100  for (int i = 0; i < cellJet.size(); ++i) {
101  eTjetsC.fill( cellJet.eT(i) );
102  etaJetsC.fill( cellJet.etaWeighted(i) );
103  phiJetsC.fill( cellJet.phiWeighted(i) );
104  }
105 
106  // Fill CellJet distance between jets.
107  for (int i = 0; i < cellJet.size() - 1; ++i)
108  for (int j = i +1; j < cellJet.size(); ++j) {
109  double dEta = cellJet.etaWeighted(i)
110  - cellJet.etaWeighted(j);
111  double dPhi = abs( cellJet.phiWeighted(i)
112  - cellJet.phiWeighted(j) );
113  if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
114  double dR = sqrt( pow2(dEta) + pow2(dPhi) );
115  distJetsC.fill( dR );
116  }
117 
118  // Fill CellJet ET-difference between jets (to check ordering of list).
119  for (int i = 1; i < cellJet.size(); ++i)
120  eTdiffC.fill( cellJet.eT(i-1) - cellJet.eT(i) );
121 
122  // Compare number of jets for the two finders.
123  nJetsD.fill( cellJet.size() - slowJet.sizeJet() );
124 
125  // End of event loop. Statistics. Histograms.
126  }
127  pythia.stat();
128  cout << nJetsS << nJetsC << nJetsD << eTjetsS << eTjetsC
129  << etaJetsS << etaJetsC << phiJetsS << phiJetsC
130  << distJetsS << distJetsC << eTdiffS << eTdiffC;
131 
132  // Done.
133  return 0;
134 }