StRoot  1
main19.cc
1 // main19.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 program runs four instances of Pythia simultaneously,
7 // one for signal events, one for pileup background ones, and two
8 // For beam-gas background ones. Note that Pythia does not do nuclear
9 // effects, so beam-gas is represented by "fixed-target" pp collisions.
10 // The = and += overloaded operators are used to join several
11 // event records into one, but should be used with caution.
12
13 // Note that each instance of Pythia is independent of any other,
14 // but with two important points to remember.
15 // 1) By default all generate the same random number sequence,
16 // which has to be corrected if they are to generate the same
17 // physics, like the two beam-gas ones below.
18 // 2) Interfaces to external Fortran programs are "by definition" static.
19 // Thus it is not a good idea to use LHAPDF to set different PDF's
20 // in different instances.
21
22 #include "Pythia.h"
23 using namespace Pythia8;
24
25 //==========================================================================
26
27 // Method to pick a number according to a Poissonian distribution.
28
29 int poisson(double nAvg, Rndm& rndm) {
30
31  // Set maximum to avoid overflow.
32  const int NMAX = 100;
33
34  // Random number.
35  double rPoisson = rndm.flat() * exp(nAvg);
36
37  // Initialize.
38  double rSum = 0.;
39  double rTerm = 1.;
40
41  // Add to sum and check whether done.
42  for (int i = 0; i < NMAX; ) {
43  rSum += rTerm;
44  if (rSum > rPoisson) return i;
45
46  // Evaluate next term.
47  ++i;
48  rTerm *= nAvg / i;
49  }
50
51  // Emergency return.
52  return NMAX;
53 }
54
55 //==========================================================================
56
57 int main() {
58
59  // Number of signal events to generate.
60  int nEvent = 100;
61
62  // Beam Energy.
63  double eBeam = 7000.;
64
65  // Average number of pileup events per signal event.
66  double nPileupAvg = 2.5;
67
68  // Average number of beam-gas events per signal ones, on two sides.
69  double nBeamAGasAvg = 0.5;
70  double nBeamBGasAvg = 0.5;
71
72  // Four generator instances.
73  Pythia pythiaSignal;
74  Pythia pythiaPileup;
75  Pythia pythiaBeamAGas;
76  Pythia pythiaBeamBGas;
77
78  // One object where all individual events are to be collected.
79  Event sumEvent;
80
81  // Switch off automatic event listing.
94
95  // Initialize generator for signal processes.
98  pythiaSignal.settings.parm("Beams:eCM", 2. * eBeam);
99  pythiaSignal.init();
100
101  // Initialize generator for pileup (background) processes.
105  pythiaPileup.settings.parm("Beams:eCM", 2. * eBeam);
106  pythiaPileup.init();
107
108  // Initialize generators for beam A - gas (background) processes.
113  pythiaBeamAGas.settings.parm("Beams:eA", eBeam);
114  pythiaBeamAGas.settings.parm("Beams:eB", 0.);
115  pythiaBeamAGas.init();
116
117  // Initialize generators for beam B - gas (background) processes.
122  pythiaBeamBGas.settings.parm("Beams:eA", 0.);
123  pythiaBeamBGas.settings.parm("Beams:eB", eBeam);
124  pythiaBeamBGas.init();
125
126  // Histograms: number of pileups, total charged multiplicity.
127  Hist nPileH("number of pileup events per signal event", 100, -0.5, 99.5);
128  Hist nAGH("number of beam A + gas events per signal event", 100, -0.5, 99.5);
129  Hist nBGH("number of beam B + gas events per signal event", 100, -0.5, 99.5);
130  Hist nChgH("number of charged multiplicity",100, -0.5, 1999.5);
131  Hist sumPZH("total pZ of system",100, -100000., 100000.);
132
133  // Loop over events.
134  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
135
136  // Generate a signal event. Copy this event into sumEvent.
137  if (!pythiaSignal.next()) continue;
138  sumEvent = pythiaSignal.event;
139
140  // Select the number of pileup events to generate.
141  int nPileup = poisson(nPileupAvg, pythiaPileup.rndm);
142  nPileH.fill( nPileup );
143
144  // Generate a number of pileup events. Add them to sumEvent.
145  for (int iPileup = 0; iPileup < nPileup; ++iPileup) {
146  pythiaPileup.next();
147  sumEvent += pythiaPileup.event;
148  }
149
150  // Select the number of beam A + gas events to generate.
151  int nBeamAGas = poisson(nBeamAGasAvg, pythiaBeamAGas.rndm);
152  nAGH.fill( nBeamAGas );
153
154  // Generate a number of beam A + gas events. Add them to sumEvent.
155  for (int iAG = 0; iAG < nBeamAGas; ++iAG) {
156  pythiaBeamAGas.next();
157  sumEvent += pythiaBeamAGas.event;
158  }
159
160  // Select the number of beam B + gas events to generate.
161  int nBeamBGas = poisson(nBeamBGasAvg, pythiaBeamBGas.rndm);
162  nBGH.fill( nBeamBGas );
163
164  // Generate a number of beam B + gas events. Add them to sumEvent.
165  for (int iBG = 0; iBG < nBeamBGas; ++iBG) {
166  pythiaBeamBGas.next();
167  sumEvent += pythiaBeamBGas.event;
168  }
169
170  // List first few events.
171  if (iEvent < 1) {
172  pythiaSignal.info.list();
173  pythiaSignal.process.list();
174  sumEvent.list();
175  }
176
177  // Find charged multiplicity.
178  int nChg = 0;
179  for (int i = 0; i < sumEvent.size(); ++i)
180  if (sumEvent[i].isFinal() && sumEvent[i].isCharged()) ++nChg;
181  nChgH.fill( nChg );
182
183  // Fill net pZ - nonvanishing owing to beam + gas.
184  sumPZH.fill( sumEvent[0].pz() );
185
186  // End of event loop
187  }
188
189  // Statistics. Histograms.
190  pythiaSignal.stat();
191  pythiaPileup.stat();
192  pythiaBeamAGas.stat();
193  pythiaBeamBGas.stat();
194  cout << nPileH << nAGH << nBGH << nChgH << sumPZH;
195
196  return 0;
197 }