]>
Commit | Line | Data |
---|---|---|
9419eeef | 1 | // main19.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2010 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 | // Initialize generator for signal processes. | |
82 | pythiaSignal.readString("HardQCD:all = on"); | |
83 | pythiaSignal.readString("PhaseSpace:pTHatMin = 50."); | |
84 | pythiaSignal.init( 2212, 2212, 2. * eBeam); | |
85 | ||
86 | // Initialize generator for pileup (background) processes. | |
87 | pythiaPileup.readString("Random:setSeed = on"); | |
88 | pythiaPileup.readString("Random:seed = 10000002"); | |
89 | pythiaPileup.readString("SoftQCD:all = on"); | |
90 | pythiaPileup.init( 2212, 2212, 2. * eBeam); | |
91 | ||
92 | // Initialize generators for beam A - gas (background) processes. | |
93 | pythiaBeamAGas.readString("Random:setSeed = on"); | |
94 | pythiaBeamAGas.readString("Random:seed = 10000003"); | |
95 | pythiaBeamAGas.readString("SoftQCD:all = on"); | |
96 | pythiaBeamAGas.init( 2212, 2212, eBeam, 0.); | |
97 | ||
98 | // Initialize generators for beam B - gas (background) processes. | |
99 | pythiaBeamBGas.readString("Random:setSeed = on"); | |
100 | pythiaBeamBGas.readString("Random:seed = 10000004"); | |
101 | pythiaBeamBGas.readString("SoftQCD:all = on"); | |
102 | pythiaBeamBGas.init( 2212, 2212, 0., eBeam); | |
103 | ||
104 | // Histograms: number of pileups, total charged multiplicity. | |
105 | Hist nPileH("number of pileup events per signal event", 100, -0.5, 99.5); | |
106 | Hist nAGH("number of beam A + gas events per signal event", 100, -0.5, 99.5); | |
107 | Hist nBGH("number of beam B + gas events per signal event", 100, -0.5, 99.5); | |
108 | Hist nChgH("number of charged multiplicity",100, -0.5, 1999.5); | |
109 | Hist sumPZH("total pZ of system",100, -100000., 100000.); | |
110 | ||
111 | // Loop over events. | |
112 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { | |
113 | ||
114 | // Generate a signal event. Copy this event into sumEvent. | |
115 | if (!pythiaSignal.next()) continue; | |
116 | sumEvent = pythiaSignal.event; | |
117 | ||
118 | // Select the number of pileup events to generate. | |
119 | int nPileup = poisson(nPileupAvg, pythiaPileup.rndm); | |
120 | nPileH.fill( nPileup ); | |
121 | ||
122 | // Generate a number of pileup events. Add them to sumEvent. | |
123 | for (int iPileup = 0; iPileup < nPileup; ++iPileup) { | |
124 | pythiaPileup.next(); | |
125 | sumEvent += pythiaPileup.event; | |
126 | } | |
127 | ||
128 | // Select the number of beam A + gas events to generate. | |
129 | int nBeamAGas = poisson(nBeamAGasAvg, pythiaBeamAGas.rndm); | |
130 | nAGH.fill( nBeamAGas ); | |
131 | ||
132 | // Generate a number of beam A + gas events. Add them to sumEvent. | |
133 | for (int iAG = 0; iAG < nBeamAGas; ++iAG) { | |
134 | pythiaBeamAGas.next(); | |
135 | sumEvent += pythiaBeamAGas.event; | |
136 | } | |
137 | ||
138 | // Select the number of beam B + gas events to generate. | |
139 | int nBeamBGas = poisson(nBeamBGasAvg, pythiaBeamBGas.rndm); | |
140 | nBGH.fill( nBeamBGas ); | |
141 | ||
142 | // Generate a number of beam B + gas events. Add them to sumEvent. | |
143 | for (int iBG = 0; iBG < nBeamBGas; ++iBG) { | |
144 | pythiaBeamBGas.next(); | |
145 | sumEvent += pythiaBeamBGas.event; | |
146 | } | |
147 | ||
148 | // List first few events. | |
149 | if (iEvent < 1) { | |
150 | pythiaSignal.info.list(); | |
151 | pythiaSignal.process.list(); | |
152 | sumEvent.list(); | |
153 | } | |
154 | ||
155 | // Find charged multiplicity. | |
156 | int nChg = 0; | |
157 | for (int i = 0; i < sumEvent.size(); ++i) | |
158 | if (sumEvent[i].isFinal() && sumEvent[i].isCharged()) ++nChg; | |
159 | nChgH.fill( nChg ); | |
160 | ||
161 | // Fill net pZ - nonvanishing owing to beam + gas. | |
162 | sumPZH.fill( sumEvent[0].pz() ); | |
163 | ||
164 | // End of event loop | |
165 | } | |
166 | ||
167 | // Statistics. Histograms. | |
168 | pythiaSignal.statistics(); | |
169 | pythiaPileup.statistics(); | |
170 | pythiaBeamAGas.statistics(); | |
171 | pythiaBeamBGas.statistics(); | |
172 | cout << nPileH << nAGH << nBGH << nChgH << sumPZH; | |
173 | ||
174 | return 0; | |
175 | } |