]>
Commit | Line | Data |
---|---|---|
c6b60c38 | 1 | // main22.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2013 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 | // Note: Pythia will automatically delete this pointer, | |
242 | // along with all other resonances. | |
243 | ResonanceWidths* resonanceTheta = new ResonanceTheta(663); | |
244 | pythia.setResonancePtr(resonanceTheta); | |
245 | ||
246 | // Create instance of a class to generate the q qbar -> Theta process | |
247 | // from an external matrix element. Hand in pointer to Pythia. | |
248 | SigmaProcess* sigma1Theta = new Sigma1qqbar2Theta(); | |
249 | pythia.setSigmaPtr(sigma1Theta); | |
250 | ||
251 | // Optionally only compare cross sections. | |
252 | //pythia.readString("PartonLevel:all = off"); | |
253 | pythia.readString("Check:nErrList = 2"); | |
254 | ||
255 | // Initialization for LHC. | |
256 | pythia.init(); | |
257 | ||
258 | // Book histogram. | |
259 | Hist mTheta("Theta mass", 100, 300., 400.); | |
260 | ||
261 | // Begin event loop. | |
262 | int iAbort = 0; | |
263 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { | |
264 | ||
265 | // Generate events. Quit if many failures. | |
266 | if (!pythia.next()) { | |
267 | if (++iAbort < nAbort) continue; | |
268 | cout << " Event generation aborted prematurely, owing to error!\n"; | |
269 | break; | |
270 | } | |
271 | ||
272 | // Fill Theta mass. End of event loop. | |
273 | mTheta.fill( pythia.process[5].m() ); | |
274 | } | |
275 | ||
276 | // Final statistics. Print histogram. | |
277 | pythia.stat(); | |
278 | cout << mTheta; | |
279 | ||
280 | // Done. | |
281 | delete sigma1Theta; | |
282 | return 0; | |
283 | } |