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. |
---|
82 | pythiaSignal.readString("Next:numberShowInfo = 0"); |
---|
83 | pythiaSignal.readString("Next:numberShowProcess = 0"); |
---|
84 | pythiaSignal.readString("Next:numberShowEvent = 0"); |
---|
85 | pythiaPileup.readString("Next:numberShowInfo = 0"); |
---|
86 | pythiaPileup.readString("Next:numberShowProcess = 0"); |
---|
87 | pythiaPileup.readString("Next:numberShowEvent = 0"); |
---|
88 | pythiaBeamAGas.readString("Next:numberShowInfo = 0"); |
---|
89 | pythiaBeamAGas.readString("Next:numberShowProcess = 0"); |
---|
90 | pythiaBeamAGas.readString("Next:numberShowEvent = 0"); |
---|
91 | pythiaBeamBGas.readString("Next:numberShowInfo = 0"); |
---|
92 | pythiaBeamBGas.readString("Next:numberShowProcess = 0"); |
---|
93 | pythiaBeamBGas.readString("Next:numberShowEvent = 0"); |
---|
94 | |
---|
95 | // Initialize generator for signal processes. |
---|
96 | pythiaSignal.readString("HardQCD:all = on"); |
---|
97 | pythiaSignal.readString("PhaseSpace:pTHatMin = 50."); |
---|
98 | pythiaSignal.settings.parm("Beams:eCM", 2. * eBeam); |
---|
99 | pythiaSignal.init(); |
---|
100 | |
---|
101 | // Initialize generator for pileup (background) processes. |
---|
102 | pythiaPileup.readString("Random:setSeed = on"); |
---|
103 | pythiaPileup.readString("Random:seed = 10000002"); |
---|
104 | pythiaPileup.readString("SoftQCD:all = on"); |
---|
105 | pythiaPileup.settings.parm("Beams:eCM", 2. * eBeam); |
---|
106 | pythiaPileup.init(); |
---|
107 | |
---|
108 | // Initialize generators for beam A - gas (background) processes. |
---|
109 | pythiaBeamAGas.readString("Random:setSeed = on"); |
---|
110 | pythiaBeamAGas.readString("Random:seed = 10000003"); |
---|
111 | pythiaBeamAGas.readString("SoftQCD:all = on"); |
---|
112 | pythiaBeamAGas.readString("Beams:frameType = 2"); |
---|
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. |
---|
118 | pythiaBeamBGas.readString("Random:setSeed = on"); |
---|
119 | pythiaBeamBGas.readString("Random:seed = 10000004"); |
---|
120 | pythiaBeamBGas.readString("SoftQCD:all = on"); |
---|
121 | pythiaBeamBGas.readString("Beams:frameType = 2"); |
---|
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 | } |
---|