1 | // main22.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 | // Simple illustration how to provide (a) your own resonance-width class, |
---|
7 | // and (b) your own cross-section class, with instances handed in to Pythia. |
---|
8 | // The hypothetical scenario is that top would have been so long-lived |
---|
9 | // that a toponium resonance Theta could form. Then production could |
---|
10 | // proceed via q qbar -> gamma*/Z* -> Theta, with decay either to |
---|
11 | // a fermion pair or (dominantly) to three gluons. |
---|
12 | // The implementation is not physically correct in any number of ways, |
---|
13 | // but should exemplify the strategy needed for realistic cases. |
---|
14 | |
---|
15 | #include "Pythia.h" |
---|
16 | |
---|
17 | using namespace Pythia8; |
---|
18 | |
---|
19 | //========================================================================== |
---|
20 | |
---|
21 | // The ResonanceTheta class handles a toponium resonance. |
---|
22 | |
---|
23 | class ResonanceTheta : public ResonanceWidths { |
---|
24 | |
---|
25 | public: |
---|
26 | |
---|
27 | // Constructor. |
---|
28 | ResonanceTheta(int idResIn) {initBasic(idResIn);} |
---|
29 | |
---|
30 | private: |
---|
31 | |
---|
32 | // Locally stored properties and couplings. |
---|
33 | double normTheta2qqbar, normTheta2llbar, normTheta2ggg; |
---|
34 | |
---|
35 | // Initialize constants. |
---|
36 | virtual void initConstants(); |
---|
37 | |
---|
38 | // Calculate various common prefactors for the current mass. |
---|
39 | // Superfluous here, so skipped. |
---|
40 | //virtual void calcPreFac(bool = false); |
---|
41 | |
---|
42 | // Calculate width for currently considered channel. |
---|
43 | virtual void calcWidth(bool = false); |
---|
44 | |
---|
45 | }; |
---|
46 | |
---|
47 | //-------------------------------------------------------------------------- |
---|
48 | |
---|
49 | // Initialize constants. |
---|
50 | |
---|
51 | void ResonanceTheta::initConstants() { |
---|
52 | |
---|
53 | // Dummy normalization of couplings to the allowed decay channels. |
---|
54 | normTheta2qqbar = 0.0001; |
---|
55 | normTheta2llbar = 0.0001; |
---|
56 | normTheta2ggg = 0.001; |
---|
57 | } |
---|
58 | |
---|
59 | //-------------------------------------------------------------------------- |
---|
60 | |
---|
61 | // Calculate width for currently considered channel. |
---|
62 | |
---|
63 | void ResonanceTheta::calcWidth(bool) { |
---|
64 | |
---|
65 | // Expression for Theta -> q qbar (q up to b). Colour factor. |
---|
66 | if (id1Abs < 6) widNow = 3. * normTheta2qqbar * mHat; |
---|
67 | |
---|
68 | // Expression for Theta -> l lbar (l = e, mu, tau). |
---|
69 | else if (id1Abs == 11 || id1Abs == 13 || id1Abs == 15) |
---|
70 | widNow = normTheta2llbar * mHat; |
---|
71 | |
---|
72 | // Expression for Theta -> g g g. Colour factor. |
---|
73 | else if (id1Abs == 21) widNow = 8. * normTheta2ggg * mHat; |
---|
74 | |
---|
75 | } |
---|
76 | |
---|
77 | //========================================================================== |
---|
78 | |
---|
79 | // A derived class for q qbar -> Theta (toponium bound state). |
---|
80 | |
---|
81 | class Sigma1qqbar2Theta : public Sigma1Process { |
---|
82 | |
---|
83 | public: |
---|
84 | |
---|
85 | // Constructor. |
---|
86 | Sigma1qqbar2Theta() {} |
---|
87 | |
---|
88 | // Initialize process. |
---|
89 | virtual void initProc(); |
---|
90 | |
---|
91 | // Calculate flavour-independent parts of cross section. |
---|
92 | virtual void sigmaKin(); |
---|
93 | |
---|
94 | // Evaluate sigmaHat(sHat). Assumed flavour-independent so simple. |
---|
95 | virtual double sigmaHat() {return sigma;} |
---|
96 | |
---|
97 | // Select flavour, colour and anticolour. |
---|
98 | virtual void setIdColAcol(); |
---|
99 | |
---|
100 | // Evaluate weight for decay angles. |
---|
101 | virtual double weightDecay( Event& process, int iResBeg, int iResEnd); |
---|
102 | |
---|
103 | // Info on the subprocess. |
---|
104 | virtual string name() const {return "q qbar -> Theta";} |
---|
105 | virtual int code() const {return 621;} |
---|
106 | virtual string inFlux() const {return "qqbarSame";} |
---|
107 | virtual int resonanceA() const {return 663;} |
---|
108 | |
---|
109 | private: |
---|
110 | |
---|
111 | // Store flavour-specific process information and standard prefactor. |
---|
112 | int idTheta; |
---|
113 | double mRes, GammaRes, m2Res, GamMRat, normTheta2qqbar, sigma; |
---|
114 | |
---|
115 | // Pointer to properties of Theta, to access decay width. |
---|
116 | ParticleDataEntry* particlePtr; |
---|
117 | |
---|
118 | }; |
---|
119 | |
---|
120 | //-------------------------------------------------------------------------- |
---|
121 | |
---|
122 | // Initialize process. |
---|
123 | |
---|
124 | void Sigma1qqbar2Theta::initProc() { |
---|
125 | |
---|
126 | // Store Theta mass and width for propagator. |
---|
127 | idTheta = 663; |
---|
128 | mRes = particleDataPtr->m0(idTheta); |
---|
129 | GammaRes = particleDataPtr->mWidth(idTheta); |
---|
130 | m2Res = mRes*mRes; |
---|
131 | GamMRat = GammaRes / mRes; |
---|
132 | |
---|
133 | // Same normlization as in ResonanceTheta for coupling strength. |
---|
134 | normTheta2qqbar = 0.0001; |
---|
135 | |
---|
136 | // Set pointer to particle properties and decay table. |
---|
137 | particlePtr = particleDataPtr->particleDataEntryPtr(idTheta); |
---|
138 | |
---|
139 | } |
---|
140 | |
---|
141 | //-------------------------------------------------------------------------- |
---|
142 | |
---|
143 | // Evaluate sigmaHat(sHat); first step when inflavours unknown. |
---|
144 | |
---|
145 | void Sigma1qqbar2Theta::sigmaKin() { |
---|
146 | |
---|
147 | // Incoming width with colour factor. |
---|
148 | double widthIn = normTheta2qqbar * mH / 3.; |
---|
149 | |
---|
150 | // Breit-Wigner, including some (guessed) spin factors. |
---|
151 | double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); |
---|
152 | |
---|
153 | // Outgoing width: only includes channels left open. |
---|
154 | double widthOut = particlePtr->resWidthOpen(663, mH); |
---|
155 | |
---|
156 | // Total answer. |
---|
157 | sigma = widthIn * sigBW * widthOut; |
---|
158 | |
---|
159 | } |
---|
160 | |
---|
161 | //-------------------------------------------------------------------------- |
---|
162 | |
---|
163 | // Select identity, colour and anticolour. |
---|
164 | |
---|
165 | void Sigma1qqbar2Theta::setIdColAcol() { |
---|
166 | |
---|
167 | // Flavours trivial. |
---|
168 | setId( id1, id2, idTheta); |
---|
169 | |
---|
170 | // Colour flow topologies. Swap when antiquarks. |
---|
171 | setColAcol( 1, 0, 0, 1, 0, 0); |
---|
172 | if (id1 < 0) swapColAcol(); |
---|
173 | |
---|
174 | } |
---|
175 | |
---|
176 | //-------------------------------------------------------------------------- |
---|
177 | |
---|
178 | // Evaluate weight for Theta -> g g g. |
---|
179 | |
---|
180 | double Sigma1qqbar2Theta::weightDecay( Event& process, int iResBeg, |
---|
181 | int iResEnd) { |
---|
182 | |
---|
183 | // Should be Theta decay. (This is only option here, so overkill.) |
---|
184 | if (iResEnd != iResBeg || process[iResBeg].idAbs() != idTheta) |
---|
185 | return 1.; |
---|
186 | |
---|
187 | // Should be decay to three gluons. |
---|
188 | int i1 = process[iResBeg].daughter1(); |
---|
189 | int i2 = i1 + 1; |
---|
190 | int i3 = i2 + 1; |
---|
191 | if (i3 != process[iResBeg].daughter2() || process[i1].id() != 21) |
---|
192 | return 1.; |
---|
193 | |
---|
194 | // Energy fractions x_i = 2 E_i/m_Theta of gluons in Theta rest frame. |
---|
195 | double x1 = 2. * process[i1].p() * process[iResBeg].p() |
---|
196 | / process[iResBeg].m2(); |
---|
197 | double x2 = 2. * process[i2].p() * process[iResBeg].p() |
---|
198 | / process[iResBeg].m2(); |
---|
199 | double x3 = 2. * process[i3].p() * process[iResBeg].p() |
---|
200 | / process[iResBeg].m2(); |
---|
201 | |
---|
202 | // Matrix-element expression for Theta -> g g g. |
---|
203 | double wtME = pow2( (1. - x1) / (x2 * x3) ) |
---|
204 | + pow2( (1. - x2) / (x1 * x3) ) + pow2( (1. - x3) / (x1 * x2) ); |
---|
205 | double wtMEmax = 2.; |
---|
206 | return wtME / wtMEmax; |
---|
207 | |
---|
208 | } |
---|
209 | |
---|
210 | //========================================================================== |
---|
211 | |
---|
212 | int main() { |
---|
213 | |
---|
214 | // Number of events to generate. Max number of errors. |
---|
215 | // Warning: generation of complete events is much slower than if you use |
---|
216 | // PartonLevel:all = off to only get cross sections, so adjust nEvent. |
---|
217 | int nEvent = 1000; |
---|
218 | int nAbort = 5; |
---|
219 | |
---|
220 | // Pythia generator. |
---|
221 | Pythia pythia; |
---|
222 | |
---|
223 | // Create the toponium resonance and a few production/decay channels. |
---|
224 | // Warning: many more exist, e.g. weak ones of one top quark. |
---|
225 | // Note: to obtain the correct width for the Breit-Wigner you must |
---|
226 | // include all channels, but you only need leave those on that you |
---|
227 | // want to study. |
---|
228 | pythia.readString("663:new = Theta void 3 0 0 342.0 0.2 300. 400. 0."); |
---|
229 | pythia.readString("663:addChannel = 1 0. 0 1 -1"); |
---|
230 | pythia.readString("663:addChannel = 1 0. 0 2 -2"); |
---|
231 | pythia.readString("663:addChannel = 1 0. 0 3 -3"); |
---|
232 | pythia.readString("663:addChannel = 1 0. 0 4 -4"); |
---|
233 | pythia.readString("663:addChannel = 1 0. 0 5 -5"); |
---|
234 | pythia.readString("663:addChannel = 1 0. 0 11 -11"); |
---|
235 | pythia.readString("663:addChannel = 1 0. 0 13 -13"); |
---|
236 | pythia.readString("663:addChannel = 1 0. 0 15 -15"); |
---|
237 | pythia.readString("663:addChannel = 1 0. 0 21 21 21"); |
---|
238 | |
---|
239 | // Create instance of a class to calculate the width of Theta to the |
---|
240 | // above channels. Hand in pointer to Pythia. |
---|
241 | ResonanceWidths* resonanceTheta = new ResonanceTheta(663); |
---|
242 | pythia.setResonancePtr(resonanceTheta); |
---|
243 | |
---|
244 | // Create instance of a class to generate the q qbar -> Theta process |
---|
245 | // from an external matrix element. Hand in pointer to Pythia. |
---|
246 | SigmaProcess* sigma1Theta = new Sigma1qqbar2Theta(); |
---|
247 | pythia.setSigmaPtr(sigma1Theta); |
---|
248 | |
---|
249 | // Optionally only compare cross sections. |
---|
250 | //pythia.readString("PartonLevel:all = off"); |
---|
251 | pythia.readString("Check:nErrList = 2"); |
---|
252 | |
---|
253 | // Initialization for LHC. |
---|
254 | pythia.init(); |
---|
255 | |
---|
256 | // Book histogram. |
---|
257 | Hist mTheta("Theta mass", 100, 300., 400.); |
---|
258 | |
---|
259 | // Begin event loop. |
---|
260 | int iAbort = 0; |
---|
261 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { |
---|
262 | |
---|
263 | // Generate events. Quit if many failures. |
---|
264 | if (!pythia.next()) { |
---|
265 | if (++iAbort < nAbort) continue; |
---|
266 | cout << " Event generation aborted prematurely, owing to error!\n"; |
---|
267 | break; |
---|
268 | } |
---|
269 | |
---|
270 | // Fill Theta mass. End of event loop. |
---|
271 | mTheta.fill( pythia.process[5].m() ); |
---|
272 | } |
---|
273 | |
---|
274 | // Final statistics. Print histogram. |
---|
275 | pythia.stat(); |
---|
276 | cout << mTheta; |
---|
277 | |
---|
278 | // Done. |
---|
279 | return 0; |
---|
280 | } |
---|