StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main06.cc
1 // main06.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 studies event properties of LEP1 events.
8 
9 #include "Pythia.h"
10 
11 using namespace Pythia8;
12 
13 int main() {
14 
15  // Generator.
16  Pythia pythia;
17 
18  // Allow no substructure in e+- beams: normal for corrected LEP data.
19  pythia.readString("PDF:lepton = off");
20  // Process selection.
21  pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
22  // Switch off all Z0 decays and then switch back on those to quarks.
23  pythia.readString("23:onMode = off");
24  pythia.readString("23:onIfAny = 1 2 3 4 5");
25 
26  // LEP1 initialization at Z0 mass.
27  pythia.readString("Beams:idA = 11");
28  pythia.readString("Beams:idB = -11");
29  double mZ = pythia.particleData.m0(23);
30  pythia.settings.parm("Beams:eCM", mZ);
31  pythia.init();
32 
33  // Histograms.
34  Hist nCharge("charged multiplicity", 100, -0.5, 99.5);
35  Hist spheri("Sphericity", 100, 0., 1.);
36  Hist linea("Linearity", 100, 0., 1.);
37  Hist thrust("thrust", 100, 0.5, 1.);
38  Hist oblateness("oblateness", 100, 0., 1.);
39  Hist sAxis("cos(theta_Sphericity)", 100, -1., 1.);
40  Hist lAxis("cos(theta_Linearity)", 100, -1., 1.);
41  Hist tAxis("cos(theta_Thrust)", 100, -1., 1.);
42  Hist nLund("Lund jet multiplicity", 40, -0.5, 39.5);
43  Hist nJade("Jade jet multiplicity", 40, -0.5, 39.5);
44  Hist nDurham("Durham jet multiplicity", 40, -0.5, 39.5);
45  Hist eDifLund("Lund e_i - e_{i+1}", 100, -5.,45.);
46  Hist eDifJade("Jade e_i - e_{i+1}", 100, -5.,45.);
47  Hist eDifDurham("Durham e_i - e_{i+1}", 100, -5.,45.);
48 
49  // Set up Sphericity, "Linearity", Thrust and cluster jet analyses.
50  Sphericity sph;
51  Sphericity lin(1.);
52  Thrust thr;
53  ClusterJet lund("Lund");
54  ClusterJet jade("Jade");
55  ClusterJet durham("Durham");
56 
57  // Begin event loop. Generate event. Skip if error. List first few.
58  for (int iEvent = 0; iEvent < 10000; ++iEvent) {
59  if (!pythia.next()) continue;
60 
61  // Find and histogram charged multiplicity.
62  int nCh = 0;
63  for (int i = 0; i < pythia.event.size(); ++i)
64  if (pythia.event[i].isFinal() && pythia.event[i].isCharged()) ++nCh;
65  nCharge.fill( nCh );
66 
67  // Find and histogram sphericity.
68  if (sph.analyze( pythia.event )) {
69  if (iEvent < 3) sph.list();
70  spheri.fill( sph.sphericity() );
71  sAxis.fill( sph.eventAxis(1).pz() );
72  double e1 = sph.eigenValue(1);
73  double e2 = sph.eigenValue(2);
74  double e3 = sph.eigenValue(3);
75  if (e2 > e1 || e3 > e2) cout << "eigenvalues out of order: "
76  << e1 << " " << e2 << " " << e3 << endl;
77  }
78 
79  // Find and histogram linearized sphericity.
80  if (lin.analyze( pythia.event )) {
81  if (iEvent < 3) lin.list();
82  linea.fill( lin.sphericity() );
83  lAxis.fill( lin.eventAxis(1).pz() );
84  double e1 = lin.eigenValue(1);
85  double e2 = lin.eigenValue(2);
86  double e3 = lin.eigenValue(3);
87  if (e2 > e1 || e3 > e2) cout << "eigenvalues out of order: "
88  << e1 << " " << e2 << " " << e3 << endl;
89  }
90 
91  // Find and histogram thrust.
92  if (thr.analyze( pythia.event )) {
93  if (iEvent < 3) thr.list();
94  thrust.fill( thr.thrust() );
95  oblateness.fill( thr.oblateness() );
96  tAxis.fill( thr.eventAxis(1).pz() );
97  if ( abs(thr.eventAxis(1).pAbs() - 1.) > 1e-8
98  || abs(thr.eventAxis(2).pAbs() - 1.) > 1e-8
99  || abs(thr.eventAxis(3).pAbs() - 1.) > 1e-8
100  || abs(thr.eventAxis(1) * thr.eventAxis(2)) > 1e-8
101  || abs(thr.eventAxis(1) * thr.eventAxis(3)) > 1e-8
102  || abs(thr.eventAxis(2) * thr.eventAxis(3)) > 1e-8 ) {
103  cout << " suspicious Thrust eigenvectors " << endl;
104  thr.list();
105  }
106  }
107 
108  // Find and histogram cluster jets: Lund, Jade and Durham distance.
109  if (lund.analyze( pythia.event, 0.01, 0.)) {
110  if (iEvent < 3) lund.list();
111  nLund.fill( lund.size() );
112  for (int j = 0; j < lund.size() - 1; ++j)
113  eDifLund.fill( lund.p(j).e() - lund.p(j+1).e() );
114  }
115  if (jade.analyze( pythia.event, 0.01, 0.)) {
116  if (iEvent < 3) jade.list();
117  nJade.fill( jade.size() );
118  for (int j = 0; j < jade.size() - 1; ++j)
119  eDifJade.fill( jade.p(j).e() - jade.p(j+1).e() );
120  }
121  if (durham.analyze( pythia.event, 0.01, 0.)) {
122  if (iEvent < 3) durham.list();
123  nDurham.fill( durham.size() );
124  for (int j = 0; j < durham.size() - 1; ++j)
125  eDifDurham.fill( durham.p(j).e() - durham.p(j+1).e() );
126  }
127 
128  // End of event loop. Statistics. Output histograms.
129  }
130  pythia.stat();
131  cout << nCharge << spheri << linea << thrust << oblateness
132  << sAxis << lAxis << tAxis
133  << nLund << nJade << nDurham
134  << eDifLund << eDifJade << eDifDurham;
135 
136  // Done.
137  return 0;
138 }