StRoot  1
main21.cc
1 // main21.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 "Pythia8/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.mSel(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(2, 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.
220
221  // Optionally switch off resonance decays, or only showers in them.
224
225  // Optionally switch off ordinary decays.
227
228  // Switch off automatic event listing in favour of manual.
232
233  // Initialize.
234  pythia.init();
235
236  // Book histograms.
237  Hist epCons("deviation from energy-momentum conservation",100,0.,1e-4);
238  Hist chgCons("deviation from charge conservation",57,-9.5,9.5);
239  Hist nFinal("final particle multiplicity",100,-0.5,99.5);
240  Hist dnparticledp("dn/dp for particles",100,0.,ee);
241  Hist status85("multiplicity status code 85",50,-0.5,49.5);
242  Hist status86("multiplicity status code 86",50,-0.5,49.5);
243  Hist status83("multiplicity status code 83",50,-0.5,49.5);
244  Hist status84("multiplicity status code 84",50,-0.5,49.5);
245  Hist dndtheta("particle flow in event plane",100,-M_PI,M_PI);
246  Hist dedtheta("energy flow in event plane",100,-M_PI,M_PI);
247  Hist dpartondtheta("parton flow in event plane",100,-M_PI,M_PI);
248  Hist dndyAnti("dn/dy primaries antijunction",100, -10., 10.);
249  Hist dndyJun("dn/dy primaries junction",100, -10., 10.);
250  Hist dndySum("dn/dy all primaries",100, -10., 10.);
251
252  // Begin of event loop.
253  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
254
255  // Set up single particle, with random direction in solid angle.
256  if (type == 0 || type == 11) fillParticle( idGun, eeGun, -1., 0.,
257  event, pdt, pythia.rndm, atRest);
258
259  // Set up parton-level configuration.
260  else fillPartons( type, ee, event, pdt, pythia.rndm);
261
262  // Generate events. Quit if failure.
263  if (!pythia.next()) {
264  cout << " Event generation aborted prematurely, owing to error!\n";
265  break;
266  }
267
268  // List first few events.
269  if (iEvent < nList) {
270  event.list();
271  // Also list junctions.
272  event.listJunctions();
273  }
274
275  // Initialize statistics.
276  Vec4 pSum = - event[0].p();
277  double chargeSum = 0.;
278  if (type == 0) chargeSum = -event[1].charge();
279  if (type == 4 || type == 5) chargeSum = -1;
280  int nFin = 0;
281  int n85 = 0;
282  int n86 = 0;
283  int n83 = 0;
284  int n84 = 0;
285
286  // Loop over all particles.
287  for (int i = 0; i < event.size(); ++i) {
288  int status = event[i].statusAbs();
289
290  // Find any unrecognized particle codes.
291  int id = event[i].id();
292  if (id == 0 || !pdt.isParticle(id))
293  cout << " Error! Unknown code id = " << id << "\n";
294
295  // Find particles with E-p mismatch.
296  double eCalc = event[i].eCalc();
297  if (abs(eCalc/event[i].e() - 1.) > 1e-6) cout << " e mismatch, i = "
298  << i << " e_nominal = " << event[i].e() << " e-from-p = "
299  << eCalc << " m-from-e " << event[i].mCalc() << "\n";
300
301  // Parton flow in event plane.
302  if (status == 71 || status == 72) {
303  double thetaXZ = event[i].thetaXZ();
304  dpartondtheta.fill(thetaXZ);
305  }
306
307  // Origin of primary hadrons.
308  if (status == 85) ++n85;
309  if (status == 86) ++n86;
310  if (status == 83) ++n83;
311  if (status == 84) ++n84;
312
313  // Flow of primary hadrons in the event plane.
314  if (status > 80 && status < 90) {
315  double eAbs = event[i].e();
316  if (eAbs < 0.) {cout << " e < 0 line " << i; event.list();}
317  double thetaXZ = event[i].thetaXZ();
318  dndtheta.fill(thetaXZ);
319  dedtheta.fill(thetaXZ, eAbs);
320
321  // Rapidity distribution of primary hadrons.
322  double y = event[i].y();
323  dndySum.fill(y);
324  if (type >= 6) {
325  int motherId = event[event[i].mother1()].id();
326  if (motherId > 0 ) dndyJun.fill(event[i].y());
327  else dndyAnti.fill(event[i].y());
328  }
329  }
330
331  // Study final-state particles.
332  if (event[i].isFinal()) {
333  pSum += event[i].p();
334  chargeSum += event[i].charge();
335  nFin++;
336  double pAbs = event[i].pAbs();
337  dnparticledp.fill(pAbs);
338  }
339  }
340
341  // Fill histograms once for each event.
342  double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
343  + abs(pSum.pz());
344  epCons.fill(epDev);
345  chgCons.fill(chargeSum);
346  nFinal.fill(nFin);
347  status85.fill(n85);
348  status86.fill(n86);
349  status83.fill(n83);
350  status84.fill(n84);
351  if (epDev > 1e-3 || abs(chargeSum) > 0.1) event.list();
352
353  // End of event loop.
354  }
355
356  // Print statistics, histograms and done.
357  pythia.stat();
358  cout << epCons << chgCons << nFinal << dnparticledp
359  << dndtheta << dedtheta << dpartondtheta << dndySum;
360  if (type >= 4) cout << status85 << status86 << status83
361  << status84;
362  if (type >= 6) cout << dndyJun << dndyAnti;
363
364  // Done.
365  return 0;
366 }