Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main72.cc
1 // main72.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 compares SlowJet and FastJet, showing that they find the same jets.
8
9 #include "Pythia.h"
10
11 // The FastJet3.h header enables automatic initialisation of 
12 // fastjet::PseudoJet objects from Pythia8 Particle and Vec4 objects, 
13 // as well as advanced features such as access to (a copy of) 
14 // the original Pythia 8 Particle directly from the PseudoJet, 
15 // and fastjet selectors that make use of the Particle properties.
16 // See the extensive comments in the header file for further details 
17 // and examples.
18 #include "FastJet3.h"
19
20 using namespace Pythia8;
21  
22 int main() {
23
24   // Number of events, generated and listed ones (for jets).
25   int nEvent    = 1000;
26   int nListJets = 3;
27
28   // Select common parameters for SlowJet and FastJet analyses. 
29   int    power   = -1;     // -1 = anti-kT; 0 = C/A; 1 = kT.
30   double R       = 0.7;    // Jet size.
31   double pTMin   = 5.0;    // Min jet pT.
32   double etaMax  = 5.0;    // Pseudorapidity range of detector.
33   int    select  = 2;      // Which particles are included?
34   int    massSet = 2;      // Which mass are they assumed to have?
35
36   // Generator. Shorthand for event.
37   Pythia pythia;
38   Event& event = pythia.event;
39
40   // Process selection.
41   pythia.readString("HardQCD:all = on");    
42   pythia.readString("PhaseSpace:pTHatMin = 200.");  
43
44   // No event record printout.
45   pythia.readString("Next:numberShowInfo = 0"); 
46   pythia.readString("Next:numberShowProcess = 0"); 
47   pythia.readString("Next:numberShowEvent = 0"); 
48
49   // LHC initialization.
50   pythia.readString("Beams:eCM = 14000.");  
51   pythia.init();
52
53   // Set up SlowJet jet finder.
54   SlowJet slowJet( power, R, pTMin, etaMax, select, massSet);
55
56   // Set up FastJet jet finder.
57   //   one can use either explicitly use antikt, cambridge, etc., or
58   //   just use genkt_algorithm with specification of power
59   //fastjet::JetAlgorithm algorithm;
60   //if (power == -1)      algorithm = fastjet::antikt_algorithm;
61   //if (power ==  0)      algorithm = fastjet::cambridge_algorithm;
62   //if (power ==  1)      algorithm = fastjet::kt_algorithm;
63   //fastjet::JetDefinition jetDef(algorithm, R);
64   // there's no need for a pointer to the jetDef (it's a fairly small object)
65   fastjet::JetDefinition jetDef(fastjet::genkt_algorithm, R, power);
66   std::vector <fastjet::PseudoJet> fjInputs;
67
68   // Histograms.
69   Hist nJetsS("number of jets, SlowJet", 100, -0.5, 99.5);
70   Hist nJetsF("number of jets, FastJet", 100, -0.5, 99.5);
71   Hist nJetsD("number of jets, SlowJet-FastJet", 99, -49.5, 49.5);
72   Hist pTjetsS("pT for jets, SlowJet", 100, 0., 500.);
73   Hist pTjetsF("pT for jets, FastJet", 100, 0., 500.);
74   Hist pTjetsD("pT for jets, SlowJet - FastJet", 100, 0., 500.);
75   Hist RdistD("R distance SlowJet to FastJet", 100, 0., 1.);
76   Hist etaJets("eta for jets", 100, -5., 5.);
77   Hist phiJets("phi for jets", 100, -M_PI, M_PI);  
78   Hist distJets("R distance between jets", 100, 0., 10.);
79   Hist pTdiff("pT difference between consecutive jets", 100, -100., 400.);
80   Hist nAna("multiplicity of analyzed event", 100, -0.5, 999.5);
81   Hist tGen("generation time as fn of multiplicity", 100, -0.5, 999.5);
82   Hist tSlow("SlowJet time as fn of multiplicity", 100, -0.5, 999.5);
83   Hist tFast("FastJet time as fn of multiplicity", 100, -0.5, 999.5);
84   Hist tSlowGen("SlowJet/generation time as fn of multiplicity", 
85     100, -0.5, 999.5);
86   Hist tFastGen("FastJet/generation time as fn of multiplicity", 
87     100, -0.5, 999.5);
88
89   // Begin event loop. Generate event. Skip if error.
90   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
91     clock_t befGen = clock();
92     if (!pythia.next()) continue;
93     clock_t aftGen = clock();
94
95     // Begin SlowJet analysis of jet properties. List first few. 
96     clock_t befSlow = clock();
97     slowJet.analyze( pythia.event );
98     clock_t aftSlow = clock();
99     if (iEvent < nListJets) slowJet.list();
100
101     // Fill inclusive SlowJet jet distributions.
102     int nSlow = slowJet.sizeJet();
103     nJetsS.fill( nSlow );
104     for (int i = 0; i < nSlow; ++i) {
105       pTjetsS.fill( slowJet.pT(i) );
106       etaJets.fill( slowJet.y(i) );
107       phiJets.fill( slowJet.phi(i) );
108     }
109
110     // Fill distance between SlowJet jets.
111     for (int i = 0; i < nSlow - 1; ++i)
112     for (int j = i + 1; j < nSlow; ++j) {
113       double dY = slowJet.y(i)  - slowJet.y(j);
114       double dPhi = abs( slowJet.phi(i) - slowJet.phi(j) );
115       if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
116       double dR = sqrt( pow2(dY) + pow2(dPhi) );
117       distJets.fill( dR );
118     }
119
120     // Fill pT-difference between SlowJet jets (to check ordering of list).
121     for (int i = 1; i < nSlow; ++i) 
122       pTdiff.fill( slowJet.pT(i-1)- slowJet.pT(i) );
123
124     // Begin FastJet analysis: extract particles from event record. 
125     clock_t befFast = clock();
126     fjInputs.resize(0);
127     Vec4   pTemp;
128     double mTemp;
129     int nAnalyze = 0;
130     for (int i = 0; i < event.size(); ++i) if (event[i].isFinal()) {
131
132       // Require visible/charged particles inside detector.
133       if      (select > 2 &&  event[i].isNeutral() ) continue;
134       else if (select == 2 && !event[i].isVisible() ) continue;
135       if (etaMax < 20. && abs(event[i].eta()) > etaMax) continue;
136      
137       // Create a PseudoJet from the complete Pythia particle.
138       fastjet::PseudoJet particleTemp = event[i];
139      
140       // Optionally modify mass and energy.
141       pTemp = event[i].p();
142       mTemp = event[i].m();
143       if (massSet < 2) {
144         mTemp = (massSet == 0 || event[i].id() == 22) ? 0. : 0.13957; 
145         pTemp.e( sqrt(pTemp.pAbs2() + mTemp*mTemp) );
146         particleTemp.reset_momentum( pTemp.px(), pTemp.py(),
147            pTemp.pz(), pTemp.e() ); 
148       }
149
150       // Store acceptable particles as input to Fastjet. 
151       // Conversion to PseudoJet is performed automatically 
152       // with the help of the code in FastJet3.h.
153       fjInputs.push_back( particleTemp);
154       ++nAnalyze;
155     }
156
157     // Run Fastjet algorithm and sort jets in pT order.
158     vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
159     fastjet::ClusterSequence clustSeq(fjInputs, jetDef);
160     inclusiveJets = clustSeq.inclusive_jets(pTMin);
161     sortedJets    = sorted_by_pt(inclusiveJets);
162     clock_t aftFast = clock();
163
164     // List first few FastJet jets and some info about them.
165     // Note: the final few columns are illustrative of what information
166     // can be extracted, but does not exhaust the possibilities.
167     if (iEvent < nListJets) {
168       cout << "\n --------  FastJet jets, p = " << setw(2) << power 
169            << "  --------------------------------------------------\n\n "
170            << "  i         pT        y      phi  mult chgmult photons"
171            << "      hardest  pT in neutral " << endl
172            << "                                                       "
173            << "  constituent        hadrons " << endl;
174       for (int i = 0; i < int(sortedJets.size()); ++i) {
175         vector<fastjet::PseudoJet> constituents 
176           = sortedJets[i].constituents();
177         fastjet::PseudoJet hardest 
178           = fastjet::SelectorNHardest(1)(constituents)[0];
179         vector<fastjet::PseudoJet> neutral_hadrons 
180           = ( fastjet::SelectorIsHadron()
181            && fastjet::SelectorIsNeutral())(constituents);
182         double neutral_hadrons_pt = join(neutral_hadrons).perp();
183         cout << setw(4) << i << fixed << setprecision(3) << setw(11) 
184              << sortedJets[i].perp() << setw(9)  << sortedJets[i].rap()
185              << setw(9) << sortedJets[i].phi_std() 
186              << setw(6) << constituents.size() 
187              << setw(8) << fastjet::SelectorIsCharged().count(constituents)
188              << setw(8) << fastjet::SelectorId(22).count(constituents)
189              << setw(13) << hardest.user_info<Particle>().name() 
190              << "     " << setw(10) << neutral_hadrons_pt << endl;
191       }
192       cout << "\n --------  End FastJet Listing  ------------------"
193            << "---------------------------------" << endl;
194     }
195  
196     // Fill inclusive FastJet jet distributions.
197     int nFast         = sortedJets.size();  
198     nJetsF.fill( nFast );
199     for (int i = 0; i < int(sortedJets.size()); ++i) 
200       pTjetsF.fill( sortedJets[i].perp() );
201
202     // Comparison of SlowJet and FastJet.
203     nJetsD.fill( nSlow - nFast);
204     if (nFast == nSlow) {
205       for (int i = 0; i < nSlow; ++i) {
206         double dist2 = pow2( slowJet.y(i) - sortedJets[i].rap())
207           + pow2( slowJet.phi(i) - sortedJets[i].phi_std()); 
208         RdistD.fill( sqrt(dist2) );    
209       }
210     }
211     
212     // Comparison of time consumption by analyzed multiplicity.
213     nAna.fill( nAnalyze);
214     tGen.fill( nAnalyze, aftGen - befGen); 
215     tSlow.fill( nAnalyze, aftSlow - befSlow); 
216     tFast.fill( nAnalyze, aftFast - befFast); 
217
218   // End of event loop.
219   }
220
221   // Statistics. Histograms. 
222   pythia.stat();
223   pTjetsD = pTjetsS - pTjetsF;
224   tSlowGen = tSlow / tGen; 
225   tFastGen = tFast / tGen; 
226   tGen  /= nAna;
227   tSlow /= nAna;
228   tFast /= nAna;
229
230   cout << nJetsS << nJetsF << nJetsD << pTjetsS << pTjetsF << pTjetsD 
231        << RdistD << etaJets << phiJets << distJets << pTdiff
232        << nAna << tGen << tSlow << tFast << tSlowGen << tFastGen;
233
234   // Done. 
235   return 0;
236 }