using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SigmaExtraDim.cxx
1 // SigmaExtraDim.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2008 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 // Function definitions (not found in the header) for the 
7 // extra-dimensional simulation classes. 
8
9 #include "SigmaExtraDim.h"
10
11 namespace Pythia8 {
12
13 //**************************************************************************
14
15 // Sigma1gg2GravitonStar class.
16 // Cross section for g g -> G* (excited graviton state). 
17
18 //*********
19
20 // Initialize process. 
21   
22 void Sigma1gg2GravitonStar::initProc() {
23
24   // Store G* mass and width for propagator. 
25   idGstar  = 5000039;
26   mRes     = ParticleDataTable::m0(idGstar);
27   GammaRes = ParticleDataTable::mWidth(idGstar);
28   m2Res    = mRes*mRes;
29   GamMRat  = GammaRes / mRes;
30
31   // Overall coupling strength kappa * m_G*.
32   kappaMG  = Settings::parm("ExtraDimensionsG*:kappaMG");
33
34   // Set pointer to particle properties and decay table.
35   gStarPtr = ParticleDataTable::particleDataPtr(idGstar);
36
37
38
39 //*********
40
41 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
42
43 void Sigma1gg2GravitonStar::sigmaKin() { 
44
45   // Incoming width for gluons.
46   double widthIn  = pow2(kappaMG) * mH / (160. * M_PI);
47
48   // Set up Breit-Wigner. Width out only includes open channels. 
49   double sigBW    = 5. * M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );    
50   double widthOut = gStarPtr->resWidthOpen(idGstar, mH);
51
52   // Modify cross section in wings of peak. Done.
53   sigma           = widthIn * sigBW * widthOut * pow2(sH / m2Res);    
54
55 }
56
57 //*********
58
59 // Select identity, colour and anticolour.
60
61 void Sigma1gg2GravitonStar::setIdColAcol() {
62
63   // Flavours trivial.
64   setId( 21, 21, idGstar);
65
66   // Colour flow topology.
67   setColAcol( 1, 2, 2, 1, 0, 0);
68
69 }
70
71 //*********
72
73 // Evaluate weight for G* decay angle.
74   
75 double Sigma1gg2GravitonStar::weightDecay( Event& process, int iResBeg, 
76   int iResEnd) {
77
78   // Identity of mother of decaying reseonance(s).
79   int idMother = process[process[iResBeg].mother1()].idAbs();
80
81   // For top decay hand over to standard routine.
82   if (idMother == 6) 
83     return weightTopDecay( process, iResBeg, iResEnd);
84
85   // G* should sit in entry 5.
86   if (iResBeg != 5 || iResEnd != 5) return 1.;
87
88   // Phase space factors. Reconstruct decay angle.
89   double mr1    = pow2(process[6].m()) / sH;
90   double mr2    = pow2(process[7].m()) / sH;
91   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
92   double cosThe = (process[3].p() - process[4].p()) 
93     * (process[7].p() - process[6].p()) / (sH * betaf);
94
95   // Default is isotropic decay.
96   double wt     = 1.;
97
98   // Angular weight for g + g -> G* -> f + fbar.
99   if (process[6].idAbs() < 19) wt = 1. - pow4(cosThe);
100
101   // Angular weight for g + g -> G* -> g + g or gamma + gamma.
102   else if (process[6].id() == 21 || process[6].id() == 22)
103     wt = (1. + 6. * pow2(cosThe) + pow4(cosThe)) / 8.;
104  
105   // Done.
106   return wt;
107
108 }
109
110 //**************************************************************************
111
112 // Sigma1ffbar2GravitonStar class.
113 // Cross section for f fbar -> G* (excited graviton state). 
114
115 //*********
116
117 // Initialize process. 
118   
119 void Sigma1ffbar2GravitonStar::initProc() {
120
121   // Store G* mass and width for propagator. 
122   idGstar  = 5000039;
123   mRes     = ParticleDataTable::m0(idGstar);
124   GammaRes = ParticleDataTable::mWidth(idGstar);
125   m2Res    = mRes*mRes;
126   GamMRat  = GammaRes / mRes;
127
128   // Overall coupling strength kappa * m_G*.
129   kappaMG  = Settings::parm("ExtraDimensionsG*:kappaMG");
130
131   // Set pointer to particle properties and decay table.
132   gStarPtr = ParticleDataTable::particleDataPtr(idGstar);
133
134
135
136 //*********
137
138 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
139
140 void Sigma1ffbar2GravitonStar::sigmaKin() { 
141
142   // Incoming width for fermions, disregarding colour factor.
143   double widthIn  = pow2(kappaMG) * mH / (80. * M_PI);
144
145   // Set up Breit-Wigner. Width out only includes open channels. 
146   double sigBW    = 5. * M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );    
147   double widthOut = gStarPtr->resWidthOpen(idGstar, mH);
148
149   // Modify cross section in wings of peak. Done.
150   sigma0          = widthIn * sigBW * widthOut * pow2(sH / m2Res);    
151
152 }
153
154 //*********
155
156 // Select identity, colour and anticolour.
157
158 void Sigma1ffbar2GravitonStar::setIdColAcol() {
159
160   // Flavours trivial.
161   setId( id1, id2, idGstar);
162
163   // Colour flow topologies. Swap when antiquarks.
164   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
165   else              setColAcol( 0, 0, 0, 0, 0, 0);
166   if (id1 < 0) swapColAcol();
167
168 }
169
170 //*********
171
172 // Evaluate weight for G* decay angle.
173   
174 double Sigma1ffbar2GravitonStar::weightDecay( Event& process, int iResBeg, 
175   int iResEnd) {
176
177   // Identity of mother of decaying reseonance(s).
178   int idMother = process[process[iResBeg].mother1()].idAbs();
179
180   // For top decay hand over to standard routine.
181   if (idMother == 6) 
182     return weightTopDecay( process, iResBeg, iResEnd);
183
184   // G* should sit in entry 5.
185   if (iResBeg != 5 || iResEnd != 5) return 1.;
186
187   // Phase space factors. Reconstruct decay angle.
188   double mr1    = pow2(process[6].m()) / sH;
189   double mr2    = pow2(process[7].m()) / sH;
190   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
191   double cosThe = (process[3].p() - process[4].p()) 
192     * (process[7].p() - process[6].p()) / (sH * betaf);
193
194   // Default is isotropic decay.
195   double wt     = 1.;
196
197   // Angular weight for f + fbar -> G* -> f + fbar.
198   if (process[6].idAbs() < 19)
199     wt = (1. - 3. * pow2(cosThe) + 4. * pow4(cosThe)) / 2.;
200
201   // Angular weight for f + fbar -> G* -> g + g or gamma + gamma.
202   else if (process[6].id() == 21 || process[6].id() == 22)
203     wt = 1. - pow4(cosThe);
204  
205   // Done.
206   return wt;
207
208 }
209
210 //**************************************************************************
211
212 // Sigma2gg2GravitonStarg class.
213 // Cross section for g g -> G* g (excited graviton state). 
214
215 //*********
216
217 // Initialize process. 
218   
219 void Sigma2gg2GravitonStarg::initProc() {
220
221   // Store G* mass and width for propagator. 
222   idGstar  = 5000039;
223   mRes     = ParticleDataTable::m0(idGstar);
224   GammaRes = ParticleDataTable::mWidth(idGstar);
225   m2Res    = mRes*mRes;
226   GamMRat  = GammaRes / mRes;
227
228   // Overall coupling strength kappa * m_G*.
229   kappaMG  = Settings::parm("ExtraDimensionsG*:kappaMG");
230
231    // Secondary open width fraction.
232   openFrac = ParticleDataTable::resOpenFrac(idGstar);
233
234
235
236 //*********
237
238 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
239
240 void Sigma2gg2GravitonStarg::sigmaKin() { 
241
242   //  Evaluate cross section. Secondary width for G*.
243   sigma = (3. * pow2(kappaMG) * alpS) / (32. * sH * s3)
244     * ( pow2(tH2 + tH * uH + uH2) / (sH2 * tH * uH) 
245     + 2. * (tH2 / uH + uH2 / tH) / sH + 3. * (tH / uH + uH / tH)
246     + 2. * (sH / uH + sH/tH) + sH2 / (tH * uH) );
247   sigma *= openFrac;
248
249 }
250
251 //*********
252
253 // Select identity, colour and anticolour.
254
255 void Sigma2gg2GravitonStarg::setIdColAcol() {
256
257   // Flavours trivial.
258   setId( 21, 21, idGstar, 21);
259
260   // Colour flow topologies: random choice between two mirrors.
261   if (Rndm::flat() < 0.5) setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
262   else                    setColAcol( 1, 2, 3, 1, 0, 0, 3, 2);
263
264 }
265
266 //*********
267
268 // Evaluate weight for decay angles: currently G* assumed isotropic.
269   
270 double Sigma2gg2GravitonStarg::weightDecay( Event& process, int iResBeg, 
271   int iResEnd) {
272
273   // Identity of mother of decaying reseonance(s).
274   int idMother = process[process[iResBeg].mother1()].idAbs();
275
276   // For top decay hand over to standard routine.
277   if (idMother == 6) 
278     return weightTopDecay( process, iResBeg, iResEnd);
279
280   // No equations for G* decay so assume isotropic.
281   return 1.;
282
283 }
284
285 //**************************************************************************
286
287 // Sigma2qg2GravitonStarq class.
288 // Cross section for q g -> G* q (excited graviton state). 
289
290 //*********
291
292 // Initialize process. 
293   
294 void Sigma2qg2GravitonStarq::initProc() {
295
296   // Store G* mass and width for propagator. 
297   idGstar  = 5000039;
298   mRes     = ParticleDataTable::m0(idGstar);
299   GammaRes = ParticleDataTable::mWidth(idGstar);
300   m2Res    = mRes*mRes;
301   GamMRat  = GammaRes / mRes;
302
303   // Overall coupling strength kappa * m_G*.
304   kappaMG  = Settings::parm("ExtraDimensionsG*:kappaMG");
305
306    // Secondary open width fraction.
307   openFrac = ParticleDataTable::resOpenFrac(idGstar);
308
309
310
311 //*********
312
313 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
314
315 void Sigma2qg2GravitonStarq::sigmaKin() { 
316
317   //  Evaluate cross section. Secondary width for G*.
318   sigma = -(pow2(kappaMG) * alpS) / (192. * sH * s3)
319     * ( 4. * (sH2 + uH2) / (tH * sH) + 9. * (sH + uH) / sH + sH / uH
320     + uH2 / sH2 + 3. * tH * (4. + sH / uH + uH / sH) / sH
321     + 4. * tH2 * (1. / uH + 1. / sH) / sH + 2. * tH2 * tH / (uH * sH2) );
322   sigma *= openFrac;
323
324 }
325
326 //*********
327
328 // Select identity, colour and anticolour.
329
330 void Sigma2qg2GravitonStarq::setIdColAcol() {
331
332   // Flavour set up for q g -> H q.
333   int idq = (id2 == 21) ? id1 : id2;
334   setId( id1, id2, idGstar, idq);
335
336   // tH defined between f and f': must swap tHat <-> uHat if q g in.
337   swapTU = (id2 == 21); 
338
339   // Colour flow topologies. Swap when antiquarks.
340   if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
341   else           setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
342   if (idq < 0) swapColAcol();
343
344 }
345
346 //*********
347
348 // Evaluate weight for decay angles: currently G* assumed isotropic.
349   
350 double Sigma2qg2GravitonStarq::weightDecay( Event& process, int iResBeg, 
351   int iResEnd) {
352
353   // Identity of mother of decaying reseonance(s).
354   int idMother = process[process[iResBeg].mother1()].idAbs();
355
356   // For top decay hand over to standard routine.
357   if (idMother == 6) 
358     return weightTopDecay( process, iResBeg, iResEnd);
359
360   // No equations for G* decay so assume isotropic.
361   return 1.;
362
363 }
364
365 //**************************************************************************
366
367 // Sigma2qqbar2GravitonStarg class.
368 // Cross section for q qbar -> G* g (excited graviton state). 
369
370 //*********
371
372 // Initialize process. 
373   
374 void Sigma2qqbar2GravitonStarg::initProc() {
375
376   // Store G* mass and width for propagator. 
377   idGstar  = 5000039;
378   mRes     = ParticleDataTable::m0(idGstar);
379   GammaRes = ParticleDataTable::mWidth(idGstar);
380   m2Res    = mRes*mRes;
381   GamMRat  = GammaRes / mRes;
382
383   // Overall coupling strength kappa * m_G*.
384   kappaMG  = Settings::parm("ExtraDimensionsG*:kappaMG");
385
386    // Secondary open width fraction.
387   openFrac = ParticleDataTable::resOpenFrac(idGstar);
388
389
390
391 //*********
392
393 // Evaluate sigmaHat(sHat), part independent of incoming flavour. 
394
395 void Sigma2qqbar2GravitonStarg::sigmaKin() { 
396
397   //  Evaluate cross section. Secondary width for G*.
398   sigma = (pow2(kappaMG) * alpS) / (72. * sH * s3)
399     * ( 4. * (tH2 + uH2) / sH2 + 9. * (tH + uH) / sH 
400     + (tH2 / uH + uH2 / tH) / sH + 3. * (4. + tH / uH + uH/ tH)
401     + 4. * (sH / uH + sH / tH) + 2. * sH2 / (tH * uH) );
402   sigma *= openFrac;
403
404 }
405
406 //*********
407
408 // Select identity, colour and anticolour.
409
410 void Sigma2qqbar2GravitonStarg::setIdColAcol() {
411
412   // Flavours trivial.
413   setId( id1, id2, idGstar, 21);
414
415   // Colour flow topologies. Swap when antiquarks.
416   setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
417   if (id1 < 0) swapColAcol();
418
419 }
420
421 //*********
422
423 // Evaluate weight for decay angles: currently G* assumed isotropic.
424   
425 double Sigma2qqbar2GravitonStarg::weightDecay( Event& process, int iResBeg, 
426   int iResEnd) {
427
428   // Identity of mother of decaying reseonance(s).
429   int idMother = process[process[iResBeg].mother1()].idAbs();
430
431   // For top decay hand over to standard routine.
432   if (idMother == 6) 
433     return weightTopDecay( process, iResBeg, iResEnd);
434
435   // No equations for G* decay so assume isotropic.
436   return 1.;
437
438 }
439
440
441 //**************************************************************************
442
443 } // end namespace Pythia8