Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main16.cc
1 // main16.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 illustrates (a) how to collect the analysis code in a separate class
8 // and (b) how to provide the .cmnd filename on the command line
9
10 // Once you have linked the main program you can run it with a command line
11 // ./main16.exe main16.cmnd > out16
12
13 #include "Pythia.h"
14
15 using namespace Pythia8; 
16
17 //==========================================================================
18
19 // Put all your own analysis code in the myAnalysis class. 
20
21 class MyAnalysis {
22
23 public:
24
25   // Constructor can be empty.
26   MyAnalysis() {}
27
28   // Initialization actions.
29   void init();
30  
31   // Analysis of each new event.
32   void analyze(Event& event);
33
34   // Show final results.
35   void finish();
36
37 private:
38
39   // Declare variables and objects that span init - analyze - finish.
40   int  nEvt;
41   Hist yH, etaChg, mult; 
42
43 };
44
45 //--------------------------------------------------------------------------
46
47 // The initialization code. 
48
49 void MyAnalysis::init() {
50
51   // Initialize counter for number of events.
52   nEvt = 0;
53
54   // Book histograms.
55   yH.book("Higgs rapidity", 100, -10., 10.);
56   etaChg.book("charged pseudorapidity", 100, -10., 10.);
57   mult.book( "charged multiplicity", 100, -0.5, 799.5);
58
59
60
61 //--------------------------------------------------------------------------
62
63 // The event analysis code. 
64
65 void MyAnalysis::analyze(Event& event) {
66
67   // Increase counter.
68   ++nEvt;
69
70   // Find latest copy of Higgs and plot its rapidity.
71   int iH = 0;
72   for (int i = 0; i < event.size(); ++i) 
73     if (event[i].id() == 25) iH = i;
74   yH.fill( event[iH].y() );
75
76   // Plot pseudorapidity distribution. Sum up charged multiplicity.
77   int nChg = 0;
78   for (int i = 0; i < event.size(); ++i) 
79   if (event[i].isFinal() && event[i].isCharged()) {
80     etaChg.fill( event[i].eta() );
81     ++nChg;
82   }
83   mult.fill( nChg );
84
85
86
87 //--------------------------------------------------------------------------
88
89 // The finishing code. 
90
91 void MyAnalysis::finish() {
92
93   // Normalize histograms.
94   double binFactor = 5. / nEvt;
95   yH     *= binFactor;
96   etaChg *= binFactor;
97
98   // Print histograms.
99   cout << yH << etaChg << mult;
100
101
102
103 //==========================================================================
104
105 // You should not need to touch the main program: its actions are 
106 // determined by the .cmnd file and the rest belongs in MyAnalysis.
107
108 int main(int argc, char* argv[]) {
109
110   // Check that correct number of command-line arguments
111   if (argc != 2) {
112     cerr << " Unexpected number of command-line arguments. \n"
113          << " You are expected to provide a file name and nothing else. \n"
114          << " Program stopped! " << endl;
115     return 1;
116   }
117
118   // Check that the provided file name corresponds to an existing file.
119   ifstream is(argv[1]);  
120   if (!is) {
121     cerr << " Command-line file " << argv[1] << " was not found. \n"
122          << " Program stopped! " << endl;
123     return 1;
124   }
125
126   // Confirm that external file will be used for settings..
127   cout << " PYTHIA settings will be read from file " << argv[1] << endl;
128
129   // Declare generator. Read in commands from external file.
130   Pythia pythia;
131   pythia.readFile(argv[1]);
132  
133   // Initialization.
134   pythia.init();
135
136   // Declare user analysis class. Do initialization part of it.
137   MyAnalysis myAnalysis;
138   myAnalysis.init(); 
139
140   // Read in number of event and maximal number of aborts.
141   int nEvent = pythia.mode("Main:numberOfEvents");
142   int nAbort = pythia.mode("Main:timesAllowErrors");
143
144   // Begin event loop.
145   int iAbort = 0; 
146   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
147
148     // Generate events. Quit if too many failures.
149     if (!pythia.next()) {
150       if (++iAbort < nAbort) continue;
151       cout << " Event generation aborted prematurely, owing to error!\n"; 
152       break;
153     }
154
155     // User Analysis of current event.
156     myAnalysis.analyze( pythia.event);
157
158   // End of event loop.
159   }
160
161   // Final statistics.
162   pythia.stat();
163
164   // User finishing.
165   myAnalysis.finish();
166
167   // Done.
168   return 0;
169 }