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