1 | // main21.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 how to feed in a single particle (including a resonance) |
---|
8 | // or a toy parton-level configurations. |
---|
9 | |
---|
10 | #include "Pythia.h" |
---|
11 | using namespace Pythia8; |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // Single-particle gun. The particle must be a colour singlet. |
---|
16 | // Input: flavour, energy, direction (theta, phi). |
---|
17 | // If theta < 0 then random choice over solid angle. |
---|
18 | // Optional final argument to put particle at rest => E = m. |
---|
19 | |
---|
20 | void fillParticle(int id, double ee, double thetaIn, double phiIn, |
---|
21 | Event& event, ParticleData& pdt, Rndm& rndm, bool atRest = false) { |
---|
22 | |
---|
23 | // Reset event record to allow for new event. |
---|
24 | event.reset(); |
---|
25 | |
---|
26 | // Select particle mass; where relevant according to Breit-Wigner. |
---|
27 | double mm = pdt.mass(id); |
---|
28 | double pp = sqrtpos(ee*ee - mm*mm); |
---|
29 | |
---|
30 | // Special case when particle is supposed to be at rest. |
---|
31 | if (atRest) { |
---|
32 | ee = mm; |
---|
33 | pp = 0.; |
---|
34 | } |
---|
35 | |
---|
36 | // Angles as input or uniform in solid angle. |
---|
37 | double cThe, sThe, phi; |
---|
38 | if (thetaIn >= 0.) { |
---|
39 | cThe = cos(thetaIn); |
---|
40 | sThe = sin(thetaIn); |
---|
41 | phi = phiIn; |
---|
42 | } else { |
---|
43 | cThe = 2. * rndm.flat() - 1.; |
---|
44 | sThe = sqrtpos(1. - cThe * cThe); |
---|
45 | phi = 2. * M_PI * rndm.flat(); |
---|
46 | } |
---|
47 | |
---|
48 | // Store the particle in the event record. |
---|
49 | event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), |
---|
50 | pp * cThe, ee, mm); |
---|
51 | |
---|
52 | } |
---|
53 | |
---|
54 | //========================================================================== |
---|
55 | |
---|
56 | // Simple method to do the filling of partons into the event record. |
---|
57 | |
---|
58 | void fillPartons(int type, double ee, Event& event, ParticleData& pdt, |
---|
59 | Rndm& rndm) { |
---|
60 | |
---|
61 | // Reset event record to allow for new event. |
---|
62 | event.reset(); |
---|
63 | |
---|
64 | // Information on a q qbar system, to be hadronized. |
---|
65 | if (type == 1) { |
---|
66 | int id = 2; |
---|
67 | double mm = pdt.m0(id); |
---|
68 | double pp = sqrtpos(ee*ee - mm*mm); |
---|
69 | event.append( id, 23, 101, 0, 0., 0., pp, ee, mm); |
---|
70 | event.append( -id, 23, 0, 101, 0., 0., -pp, ee, mm); |
---|
71 | |
---|
72 | // Information on a g g system, to be hadronized. |
---|
73 | } else if (type == 2) { |
---|
74 | event.append( 21, 23, 101, 102, 0., 0., ee, ee); |
---|
75 | event.append( 21, 23, 102, 101, 0., 0., -ee, ee); |
---|
76 | |
---|
77 | // Information on a g g g system, to be hadronized. |
---|
78 | } else if (type == 3) { |
---|
79 | event.append( 21, 23, 101, 102, 0., 0., ee, ee); |
---|
80 | event.append( 21, 23, 102, 103, 0.8 * ee, 0., -0.6 * ee, ee); |
---|
81 | event.append( 21, 23, 103, 101, -0.8 * ee, 0., -0.6 * ee, ee); |
---|
82 | |
---|
83 | // Information on a q q q junction system, to be hadronized. |
---|
84 | } else if (type == 4 || type == 5) { |
---|
85 | |
---|
86 | // Need a colour singlet mother parton to define junction origin. |
---|
87 | event.append( 1000022, -21, 0, 0, 2, 4, 0, 0, |
---|
88 | 0., 0., 1.01 * ee, 1.01 * ee); |
---|
89 | |
---|
90 | // The three endpoint q q q; the minimal system. |
---|
91 | double rt75 = sqrt(0.75); |
---|
92 | event.append( 2, 23, 1, 0, 0, 0, 101, 0, |
---|
93 | 0., 0., 1.01 * ee, 1.01 * ee); |
---|
94 | event.append( 2, 23, 1, 0, 0, 0, 102, 0, |
---|
95 | rt75 * ee, 0., -0.5 * ee, ee ); |
---|
96 | event.append( 1, 23, 1, 0, 0, 0, 103, 0, |
---|
97 | -rt75 * ee, 0., -0.5 * ee, ee ); |
---|
98 | |
---|
99 | // Define the qqq configuration as starting point for adding gluons. |
---|
100 | if (type == 5) { |
---|
101 | int colNow[4] = {0, 101, 102, 103}; |
---|
102 | Vec4 pQ[4]; |
---|
103 | pQ[1] = Vec4(0., 0., 1., 0.); |
---|
104 | pQ[2] = Vec4( rt75, 0., -0.5, 0.); |
---|
105 | pQ[3] = Vec4(-rt75, 0., -0.5, 0.); |
---|
106 | |
---|
107 | // Minimal cos(q-g opening angle), allows more or less nasty events. |
---|
108 | double cosThetaMin =0.; |
---|
109 | |
---|
110 | // Add a few gluons (almost) at random. |
---|
111 | for (int nglu = 0; nglu < 5; ++nglu) { |
---|
112 | int iq = 1 + int( 2.99999 * rndm.flat() ); |
---|
113 | double px, py, pz, e, prod; |
---|
114 | do { |
---|
115 | e = ee * rndm.flat(); |
---|
116 | double cThe = 2. * rndm.flat() - 1.; |
---|
117 | double phi = 2. * M_PI * rndm.flat(); |
---|
118 | px = e * sqrt(1. - cThe*cThe) * cos(phi); |
---|
119 | py = e * sqrt(1. - cThe*cThe) * sin(phi); |
---|
120 | pz = e * cThe; |
---|
121 | prod = ( px * pQ[iq].px() + py * pQ[iq].py() + pz * pQ[iq].pz() ) |
---|
122 | / e; |
---|
123 | } while (prod < cosThetaMin); |
---|
124 | int colNew = 104 + nglu; |
---|
125 | event.append( 21, 23, 1, 0, 0, 0, colNew, colNow[iq], |
---|
126 | px, py, pz, e, 0.); |
---|
127 | colNow[iq] = colNew; |
---|
128 | } |
---|
129 | // Update daughter range of mother. |
---|
130 | event[1].daughters(1, event.size() - 1); |
---|
131 | |
---|
132 | } |
---|
133 | |
---|
134 | // Information on a q q qbar qbar dijunction system, to be hadronized. |
---|
135 | } else if (type >= 6) { |
---|
136 | |
---|
137 | // The two fictitious beam remnant particles; needed for junctions. |
---|
138 | event.append( 2212, -12, 0, 0, 3, 5, 0, 0, 0., 0., ee, ee, 0.); |
---|
139 | event.append(-2212, -12, 0, 0, 6, 8, 0, 0, 0., 0., ee, ee, 0.); |
---|
140 | |
---|
141 | // Opening angle between "diquark" legs. |
---|
142 | double theta = 0.2; |
---|
143 | double cThe = cos(theta); |
---|
144 | double sThe = sin(theta); |
---|
145 | |
---|
146 | // Set one colour depending on whether more gluons or not. |
---|
147 | int acol = (type == 6) ? 103 : 106; |
---|
148 | |
---|
149 | // The four endpoint q q qbar qbar; the minimal system. |
---|
150 | // Two additional fictitious partons to make up original beams. |
---|
151 | event.append( 2, 23, 1, 0, 0, 0, 101, 0, |
---|
152 | ee * sThe, 0., ee * cThe, ee, 0.); |
---|
153 | event.append( 1, 23, 1, 0, 0, 0, 102, 0, |
---|
154 | -ee * sThe, 0., ee * cThe, ee, 0.); |
---|
155 | event.append( 2, -21, 1, 0, 0, 0, 103, 0, |
---|
156 | 0., 0., ee , ee, 0.); |
---|
157 | event.append( -2, 23, 2, 0, 0, 0, 0, 104, |
---|
158 | ee * sThe, 0., -ee * cThe, ee, 0.); |
---|
159 | event.append( -1, 23, 2, 0, 0, 0, 0, 105, |
---|
160 | -ee * sThe, 0., -ee * cThe, ee, 0.); |
---|
161 | event.append( -2, -21, 2, 0, 0, 0, 0, acol, |
---|
162 | 0., 0., -ee , ee, 0.); |
---|
163 | |
---|
164 | // Add extra gluons on string between junctions. |
---|
165 | if (type == 7) { |
---|
166 | event.append( 21, 23, 5, 8, 0, 0, 103, 106, 0., ee, 0., ee, 0.); |
---|
167 | } else if (type == 8) { |
---|
168 | event.append( 21, 23, 5, 8, 0, 0, 103, 108, 0., ee, 0., ee, 0.); |
---|
169 | event.append( 21, 23, 5, 8, 0, 0, 108, 106, 0.,-ee, 0., ee, 0.); |
---|
170 | } else if (type == 9) { |
---|
171 | event.append( 21, 23, 5, 8, 0, 0, 103, 107, 0., ee, 0., ee, 0.); |
---|
172 | event.append( 21, 23, 5, 8, 0, 0, 107, 108, ee, 0., 0., ee, 0.); |
---|
173 | event.append( 21, 23, 5, 8, 0, 0, 108, 106, 0.,-ee, 0., ee, 0.); |
---|
174 | } else if (type == 10) { |
---|
175 | event.append( 21, 23, 5, 8, 0, 0, 103, 107, 0., ee, 0., ee, 0.); |
---|
176 | event.append( 21, 23, 5, 8, 0, 0, 107, 108, ee, 0., 0., ee, 0.); |
---|
177 | event.append( 21, 23, 5, 8, 0, 0, 108, 109, 0.,-ee, 0., ee, 0.); |
---|
178 | event.append( 21, 23, 5, 8, 0, 0, 109, 106,-ee, 0., 0., ee, 0.); |
---|
179 | } |
---|
180 | |
---|
181 | // No more cases: done. |
---|
182 | } |
---|
183 | } |
---|
184 | |
---|
185 | //========================================================================== |
---|
186 | |
---|
187 | int main() { |
---|
188 | |
---|
189 | // Pick kind of events to generate: |
---|
190 | // 0 = single-particle gun. |
---|
191 | // 1 = q qbar. |
---|
192 | // 2 = g g. |
---|
193 | // 3 = g g g. |
---|
194 | // 4 = minimal q q q junction topology. |
---|
195 | // 5 = q q q junction topology with gluons on the strings. |
---|
196 | // 6 = q q qbar qbar dijunction topology, no gluons. |
---|
197 | // 7 - 10 : ditto, but with 1 - 4 gluons on string between junctions. |
---|
198 | // 11 : single-resonance gun. |
---|
199 | int type = 11; |
---|
200 | |
---|
201 | // Set particle species and energy for single-particle gun. |
---|
202 | int idGun = (type == 0) ? 15 : 25; |
---|
203 | double eeGun = (type == 0) ? 20. : 125.; |
---|
204 | bool atRest = (type == 0) ? false : true; |
---|
205 | |
---|
206 | // Set typical energy per parton. |
---|
207 | double ee = 20.0; |
---|
208 | |
---|
209 | // Set number of events to generate and to list. |
---|
210 | int nEvent = 10000; |
---|
211 | int nList = 3; |
---|
212 | |
---|
213 | // Generator; shorthand for event and particleData. |
---|
214 | Pythia pythia; |
---|
215 | Event& event = pythia.event; |
---|
216 | ParticleData& pdt = pythia.particleData; |
---|
217 | |
---|
218 | // Key requirement: switch off ProcessLevel, and thereby also PartonLevel. |
---|
219 | pythia.readString("ProcessLevel:all = off"); |
---|
220 | |
---|
221 | // Also allow resonance decays, with showers in them |
---|
222 | pythia.readString("Standalone:allowResDec = on"); |
---|
223 | |
---|
224 | // Optionally switch off decays. |
---|
225 | //pythia.readString("HadronLevel:Decay = off"); |
---|
226 | |
---|
227 | // Switch off automatic event listing in favour of manual. |
---|
228 | pythia.readString("Next:numberShowInfo = 0"); |
---|
229 | pythia.readString("Next:numberShowProcess = 0"); |
---|
230 | pythia.readString("Next:numberShowEvent = 0"); |
---|
231 | |
---|
232 | // Initialize. |
---|
233 | pythia.init(); |
---|
234 | |
---|
235 | // Book histograms. |
---|
236 | Hist epCons("deviation from energy-momentum conservation",100,0.,1e-4); |
---|
237 | Hist chgCons("deviation from charge conservation",57,-9.5,9.5); |
---|
238 | Hist nFinal("final particle multiplicity",100,-0.5,99.5); |
---|
239 | Hist dnparticledp("dn/dp for particles",100,0.,ee); |
---|
240 | Hist status85("multiplicity status code 85",50,-0.5,49.5); |
---|
241 | Hist status86("multiplicity status code 86",50,-0.5,49.5); |
---|
242 | Hist status83("multiplicity status code 83",50,-0.5,49.5); |
---|
243 | Hist status84("multiplicity status code 84",50,-0.5,49.5); |
---|
244 | Hist dndtheta("particle flow in event plane",100,-M_PI,M_PI); |
---|
245 | Hist dedtheta("energy flow in event plane",100,-M_PI,M_PI); |
---|
246 | Hist dpartondtheta("parton flow in event plane",100,-M_PI,M_PI); |
---|
247 | Hist dndyAnti("dn/dy primaries antijunction",100, -10., 10.); |
---|
248 | Hist dndyJun("dn/dy primaries junction",100, -10., 10.); |
---|
249 | Hist dndySum("dn/dy all primaries",100, -10., 10.); |
---|
250 | |
---|
251 | // Begin of event loop. |
---|
252 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { |
---|
253 | |
---|
254 | // Set up single particle, with random direction in solid angle. |
---|
255 | if (type == 0 || type == 11) fillParticle( idGun, eeGun, -1., 0., |
---|
256 | event, pdt, pythia.rndm, atRest); |
---|
257 | |
---|
258 | // Set up parton-level configuration. |
---|
259 | else fillPartons( type, ee, event, pdt, pythia.rndm); |
---|
260 | |
---|
261 | // Generate events. Quit if failure. |
---|
262 | if (!pythia.next()) { |
---|
263 | cout << " Event generation aborted prematurely, owing to error!\n"; |
---|
264 | break; |
---|
265 | } |
---|
266 | |
---|
267 | // List first few events. |
---|
268 | if (iEvent < nList) { |
---|
269 | event.list(); |
---|
270 | // Also list junctions. |
---|
271 | event.listJunctions(); |
---|
272 | } |
---|
273 | |
---|
274 | // Initialize statistics. |
---|
275 | Vec4 pSum = - event[0].p(); |
---|
276 | double chargeSum = 0.; |
---|
277 | if (type == 0) chargeSum = -event[1].charge(); |
---|
278 | if (type == 4 || type == 5) chargeSum = -1; |
---|
279 | int nFin = 0; |
---|
280 | int n85 = 0; |
---|
281 | int n86 = 0; |
---|
282 | int n83 = 0; |
---|
283 | int n84 = 0; |
---|
284 | |
---|
285 | // Loop over all particles. |
---|
286 | for (int i = 0; i < event.size(); ++i) { |
---|
287 | int status = event[i].statusAbs(); |
---|
288 | |
---|
289 | // Find any unrecognized particle codes. |
---|
290 | int id = event[i].id(); |
---|
291 | if (id == 0 || !pdt.isParticle(id)) |
---|
292 | cout << " Error! Unknown code id = " << id << "\n"; |
---|
293 | |
---|
294 | // Find particles with E-p mismatch. |
---|
295 | double eCalc = event[i].eCalc(); |
---|
296 | if (abs(eCalc/event[i].e() - 1.) > 1e-6) cout << " e mismatch, i = " |
---|
297 | << i << " e_nominal = " << event[i].e() << " e-from-p = " |
---|
298 | << eCalc << " m-from-e " << event[i].mCalc() << "\n"; |
---|
299 | |
---|
300 | // Parton flow in event plane. |
---|
301 | if (status == 71 || status == 72) { |
---|
302 | double thetaXZ = event[i].thetaXZ(); |
---|
303 | dpartondtheta.fill(thetaXZ); |
---|
304 | } |
---|
305 | |
---|
306 | // Origin of primary hadrons. |
---|
307 | if (status == 85) ++n85; |
---|
308 | if (status == 86) ++n86; |
---|
309 | if (status == 83) ++n83; |
---|
310 | if (status == 84) ++n84; |
---|
311 | |
---|
312 | // Flow of primary hadrons in the event plane. |
---|
313 | if (status > 80 && status < 90) { |
---|
314 | double eAbs = event[i].e(); |
---|
315 | if (eAbs < 0.) {cout << " e < 0 line " << i; event.list();} |
---|
316 | double thetaXZ = event[i].thetaXZ(); |
---|
317 | dndtheta.fill(thetaXZ); |
---|
318 | dedtheta.fill(thetaXZ, eAbs); |
---|
319 | |
---|
320 | // Rapidity distribution of primary hadrons. |
---|
321 | double y = event[i].y(); |
---|
322 | dndySum.fill(y); |
---|
323 | if (type >= 6) { |
---|
324 | int motherId = event[event[i].mother1()].id(); |
---|
325 | if (motherId > 0 ) dndyJun.fill(event[i].y()); |
---|
326 | else dndyAnti.fill(event[i].y()); |
---|
327 | } |
---|
328 | } |
---|
329 | |
---|
330 | // Study final-state particles. |
---|
331 | if (event[i].isFinal()) { |
---|
332 | pSum += event[i].p(); |
---|
333 | chargeSum += event[i].charge(); |
---|
334 | nFin++; |
---|
335 | double pAbs = event[i].pAbs(); |
---|
336 | dnparticledp.fill(pAbs); |
---|
337 | } |
---|
338 | } |
---|
339 | |
---|
340 | // Fill histograms once for each event. |
---|
341 | double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py()) |
---|
342 | + abs(pSum.pz()); |
---|
343 | epCons.fill(epDev); |
---|
344 | chgCons.fill(chargeSum); |
---|
345 | nFinal.fill(nFin); |
---|
346 | status85.fill(n85); |
---|
347 | status86.fill(n86); |
---|
348 | status83.fill(n83); |
---|
349 | status84.fill(n84); |
---|
350 | if (epDev > 1e-3 || abs(chargeSum) > 0.1) event.list(); |
---|
351 | |
---|
352 | // End of event loop. |
---|
353 | } |
---|
354 | |
---|
355 | // Print statistics, histograms and done. |
---|
356 | pythia.stat(); |
---|
357 | cout << epCons << chgCons << nFinal << dnparticledp |
---|
358 | << dndtheta << dedtheta << dpartondtheta << dndySum; |
---|
359 | if (type >= 4) cout << status85 << status86 << status83 |
---|
360 | << status84; |
---|
361 | if (type >= 6) cout << dndyJun << dndyAnti; |
---|
362 | |
---|
363 | // Done. |
---|
364 | return 0; |
---|
365 | } |
---|