1 | // main15.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 either |
---|
8 | // (a) B decays (sections marked "Repeated decays:), or |
---|
9 | // (b) all hadronization (sections marked "Repeated hadronization:") |
---|
10 | // could be repeated a number of times for each event, |
---|
11 | // to improve statistics when this could be a problem. |
---|
12 | // Option (a) is faster than (b), but less generic. |
---|
13 | |
---|
14 | #include "Pythia.h" |
---|
15 | using namespace Pythia8; |
---|
16 | |
---|
17 | int main() { |
---|
18 | |
---|
19 | // Main switches: redo B decays only or redo all hadronization, but not both. |
---|
20 | bool redoBDecays = false; |
---|
21 | bool redoHadrons = true; |
---|
22 | if (redoHadrons) redoBDecays = false; |
---|
23 | |
---|
24 | // Number of events. Number to list redone events. |
---|
25 | int nEvent = 100; |
---|
26 | int nListRedo = 1; |
---|
27 | |
---|
28 | // Number of times decays/hadronization should be redone for each event. |
---|
29 | int nRepeat = 10; |
---|
30 | if (!redoBDecays && !redoHadrons) nRepeat = 1; |
---|
31 | |
---|
32 | // Generator. Shorthand for event. |
---|
33 | Pythia pythia; |
---|
34 | Event& event = pythia.event; |
---|
35 | |
---|
36 | // Simulate b production above given pTmin scale. |
---|
37 | // Warning: these processes do not catch all possible production modes. |
---|
38 | // You would need to use HardQCD:all or even SoftQCD:minBias for that. |
---|
39 | pythia.readString("HardQCD:gg2bbbar = on"); |
---|
40 | pythia.readString("HardQCD:qqbar2bbbar = on"); |
---|
41 | pythia.readString("PhaseSpace:pTHatMin = 50."); |
---|
42 | |
---|
43 | // Repeated decays: list of weakly decaying B hadrons. |
---|
44 | // Note: this list is overkill; some will never be produced. |
---|
45 | int bCodes[28] = {511, 521, 531, 541, 5122, 5132, 5142, 5232, 5242, |
---|
46 | 5332, 5342, 5412, 5414, 5422, 5424, 5432, 5434, 5442, 5444, 5512, |
---|
47 | 5514, 5522, 5524, 5532, 5534, 5542, 5544, 5544 }; |
---|
48 | int nCodes = 28; |
---|
49 | |
---|
50 | // Repeated decays: location of B handrons. |
---|
51 | vector<int> iBHad; |
---|
52 | int nBHad = 0; |
---|
53 | |
---|
54 | // Repeated hadronization: spare copy of event. |
---|
55 | Event savedEvent; |
---|
56 | |
---|
57 | // Repeated hadronization: switch off normal HadronLevel call. |
---|
58 | if (redoHadrons) pythia.readString("HadronLevel:all = off"); |
---|
59 | |
---|
60 | // Initialize for LHC energies; default 14 TeV |
---|
61 | pythia.init(); |
---|
62 | |
---|
63 | // Histogram invariant mass of muon pairs. |
---|
64 | Hist nBperEvent("number of b quarks in an event", 10, -0.5, 9.5); |
---|
65 | Hist nSameEvent("number of times same event is used", 10, -0.5, 9.5); |
---|
66 | Hist oppSignMass("mass of opposite-sign muon pair", 100, 0.0, 100.0); |
---|
67 | Hist sameSignMass("mass of same-sign muon pair", 100, 0.0, 100.0); |
---|
68 | |
---|
69 | // Begin event loop. |
---|
70 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { |
---|
71 | |
---|
72 | // Repeated decays: switch off decays of weakly decaying B hadrons. |
---|
73 | // (More compact solution than repeated readString(..).) |
---|
74 | if (redoBDecays) for (int iC = 0; iC < nCodes; ++iC) |
---|
75 | pythia.particleData.mayDecay( bCodes[iC], false); |
---|
76 | |
---|
77 | // Generate event. Skip it if error. |
---|
78 | if (!pythia.next()) continue; |
---|
79 | |
---|
80 | // Find and histogram number of b quarks. |
---|
81 | int nBquark = 0; |
---|
82 | int stat; |
---|
83 | for (int i = 0; i < event.size(); ++i) { |
---|
84 | stat = event[i].statusAbs(); |
---|
85 | if (event[i].idAbs() == 5 && (stat == 62 || stat == 63)) ++nBquark; |
---|
86 | } |
---|
87 | nBperEvent.fill( nBquark ); |
---|
88 | |
---|
89 | // Repeated decays: find all locations where B hadrons are stored. |
---|
90 | if (redoBDecays) { |
---|
91 | iBHad.resize(0); |
---|
92 | for (int i = 0; i < event.size(); ++i) { |
---|
93 | int idAbs = event[i].idAbs(); |
---|
94 | for (int iC = 0; iC < 28; ++iC) |
---|
95 | if (idAbs == bCodes[iC]) { |
---|
96 | iBHad.push_back(i); |
---|
97 | break; |
---|
98 | } |
---|
99 | } |
---|
100 | |
---|
101 | // Repeated decays: check that #b = #B. |
---|
102 | nBHad = iBHad.size(); |
---|
103 | if (nBquark != nBHad) cout << " Warning: " << nBquark |
---|
104 | << " b quarks but " << nBHad << " B hadrons" << endl; |
---|
105 | |
---|
106 | // Repeated decays: store size of current event. |
---|
107 | event.saveSize(); |
---|
108 | |
---|
109 | // Repeated decays: switch back on weakly decaying B hadrons. |
---|
110 | for (int iC = 0; iC < nCodes; ++iC) |
---|
111 | pythia.particleData.mayDecay( bCodes[iC], true); |
---|
112 | |
---|
113 | // Repeated hadronization: copy event into spare position. |
---|
114 | } else if (redoHadrons) { |
---|
115 | savedEvent = event; |
---|
116 | } |
---|
117 | |
---|
118 | // Begin loop over rounds of decays / hadronization for same event. |
---|
119 | int nWithPair = 0; |
---|
120 | for (int iRepeat = 0; iRepeat < nRepeat; ++iRepeat) { |
---|
121 | |
---|
122 | // Repeated decays: remove B decay products from previous round. |
---|
123 | if (redoBDecays) { |
---|
124 | if (iRepeat > 0) { |
---|
125 | event.restoreSize(); |
---|
126 | |
---|
127 | // Repeated decays: mark decayed B hadrons as undecayed. |
---|
128 | for (int iB = 0; iB < nBHad; ++iB) event[ iBHad[iB] ].statusPos(); |
---|
129 | } |
---|
130 | |
---|
131 | // Repeated decays: do decays of B hadrons, sequentially for products. |
---|
132 | // Note: modeDecays does not work for bottomonium (or heavier) states, |
---|
133 | // since there decays like Upsilon -> g g g also need hadronization. |
---|
134 | // Also, there is no provision for Bose-Einstein effects. |
---|
135 | if (!pythia.moreDecays()) continue; |
---|
136 | |
---|
137 | |
---|
138 | // Repeated hadronization: restore saved event record. |
---|
139 | } else if (redoHadrons) { |
---|
140 | if (iRepeat > 0) event = savedEvent; |
---|
141 | |
---|
142 | // Repeated hadronization: do HadronLevel (repeatedly). |
---|
143 | // Note: argument false needed owing to bug in junction search?? |
---|
144 | if (!pythia.forceHadronLevel(false)) continue; |
---|
145 | } |
---|
146 | |
---|
147 | // List last repetition of first few events. |
---|
148 | if ( (redoBDecays || redoHadrons) && iEvent < nListRedo |
---|
149 | && iRepeat == nRepeat - 1) event.list(); |
---|
150 | |
---|
151 | // Look for muons among decay products (also from charm/tau/...). |
---|
152 | vector<int> iMuNeg, iMuPos; |
---|
153 | for (int i = 0; i < event.size(); ++i) { |
---|
154 | int id = event[i].id(); |
---|
155 | if (id == 13) iMuNeg.push_back(i); |
---|
156 | if (id == -13) iMuPos.push_back(i); |
---|
157 | } |
---|
158 | |
---|
159 | // Check whether pair(s) present. |
---|
160 | int nMuNeg = iMuNeg.size(); |
---|
161 | int nMuPos = iMuPos.size(); |
---|
162 | if (nMuNeg + nMuPos > 1) { |
---|
163 | ++nWithPair; |
---|
164 | |
---|
165 | // Fill masses of opposite-sign pairs. |
---|
166 | for (int iN = 0; iN < nMuNeg; ++iN) |
---|
167 | for (int iP = 0; iP < nMuPos; ++iP) |
---|
168 | oppSignMass.fill( |
---|
169 | (event[iMuNeg[iN]].p() + event[iMuPos[iP]].p()).mCalc() ); |
---|
170 | |
---|
171 | // Fill masses of same-sign pairs. |
---|
172 | for (int i1 = 0; i1 < nMuNeg - 1; ++i1) |
---|
173 | for (int i2 = i1 + 1; i2 < nMuNeg; ++i2) |
---|
174 | sameSignMass.fill( |
---|
175 | (event[iMuNeg[i1]].p() + event[iMuNeg[i2]].p()).mCalc() ); |
---|
176 | for (int i1 = 0; i1 < nMuPos - 1; ++i1) |
---|
177 | for (int i2 = i1 + 1; i2 < nMuPos; ++i2) |
---|
178 | sameSignMass.fill( |
---|
179 | (event[iMuPos[i1]].p() + event[iMuPos[i2]].p()).mCalc() ); |
---|
180 | |
---|
181 | // Finished analysis of current round. |
---|
182 | } |
---|
183 | |
---|
184 | // End of loop over many rounds. fill number of rounds with pairs. |
---|
185 | } |
---|
186 | nSameEvent.fill( nWithPair ); |
---|
187 | |
---|
188 | // End of event loop. |
---|
189 | } |
---|
190 | |
---|
191 | // Statistics. Histograms. |
---|
192 | pythia.stat(); |
---|
193 | cout << nBperEvent << nSameEvent << oppSignMass << sameSignMass << endl; |
---|
194 | |
---|
195 | // Done. |
---|
196 | return 0; |
---|
197 | } |
---|