StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main22.cc
1 // main22.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 "Pythia8/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.
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 }