Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main05.cc
1 // main05.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 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 "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 }