Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main22.cc
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 }