StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaGeneric.cc
1 // SigmaGeneric.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 Johan Bijnens, 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 various generic
7 // production processes, to be used as building blocks for some BSM processes.
8 // Currently represented by QCD pair production of colour triplet objects,
9 // with spin either 0, 1/2 or 1.
10 
11 // Cross sections are only provided for fixed m3 = m4, so do some gymnastics:
12 // i) s34Avg picked so that beta34 same when s3, s4 -> s34Avg.
13 // ii) tHQ = tH - mQ^2 = -0.5 sH (1 - beta34 cos(thetaH)) for m3 = m4 = mQ,
14 // but tH - uH = sH beta34 cos(thetaH) also for m3 != m4, so use
15 // tH, uH selected for m3 != m4 to derive tHQ, uHQ valid for m3 = m4.
16 
17 #include "SigmaGeneric.h"
18 
19 namespace Pythia8 {
20 
21 //==========================================================================
22 
23 // Sigma2gg2qGqGbar class.
24 // Cross section for g g -> qG qGbar (generic quark of spin 0, 1/2 or 1).
25 
26 //--------------------------------------------------------------------------
27 
28 // Initialize process.
29 
30 void Sigma2gg2qGqGbar::initProc() {
31 
32  // Number of colours. Anomalous coupling kappa - 1 used for vector state.
33  nCHV = settingsPtr->mode("HiddenValley:Ngauge");
34  kappam1 = settingsPtr->parm("HiddenValley:kappa") - 1.;
35  hasKappa = (abs(kappam1) > 1e-8);
36 
37  // Secondary open width fraction.
38  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
39 
40 }
41 
42 //--------------------------------------------------------------------------
43 
44 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
45 
46 void Sigma2gg2qGqGbar::sigmaKin() {
47 
48  // Modified Mandelstam variables for massive kinematics with m3 = m4.
49  double delta = 0.25 * pow2(s3 - s4) / sH;
50  double s34Avg = 0.5 * (s3 + s4) - delta;
51  double tHavg = tH - delta;
52  double uHavg = uH - delta;
53  double tHQ = -0.5 * (sH - tH + uH);
54  double uHQ = -0.5 * (sH + tH - uH);
55  double tHQ2 = tHQ * tHQ;
56  double uHQ2 = uHQ * uHQ;
57 
58  // Evaluate cross section for spin 0 colour triplet.
59  if (spinSave == 0) {
60  sigSum = 0.5 * ( 7. / 48. + 3. * pow2(uHavg - tHavg) / (16. * sH2) )
61  * ( 1. + 2. * s34Avg * tHavg / pow2(tHavg - s34Avg)
62  + 2. * s34Avg * uHavg / pow2(uHavg - s34Avg)
63  + 4. * pow2(s34Avg) / ((tHavg - s34Avg) * (uHavg - s34Avg)) );
64 
65  // Equal probability for two possible colour flows.
66  sigTS = 0.5 * sigSum;
67  sigUS = sigTS;
68  }
69 
70  // Evaluate cross section for spin 1/2 colour triplet.
71  else if (spinSave == 1) {
72  double tumHQ = tHQ * uHQ - s34Avg * sH;
73  sigTS = ( uHQ / tHQ - 2.25 * uHQ2 / sH2 + 4.5 * s34Avg * tumHQ
74  / ( sH * tHQ2) + 0.5 * s34Avg * (tHQ + s34Avg) / tHQ2
75  - s34Avg*s34Avg / (sH * tHQ) ) / 6.;
76  sigUS = ( tHQ / uHQ - 2.25 * tHQ2 / sH2 + 4.5 * s34Avg * tumHQ
77  / ( sH * uHQ2) + 0.5 * s34Avg * (uHQ + s34Avg) / uHQ2
78  - s34Avg*s34Avg / (sH * uHQ) ) / 6.;
79  sigSum = sigTS + sigUS;
80  }
81 
82  // Evaluate cross section for spin 1 colour triplet.
83  else {
84  double tmu = tHavg - uHavg;
85  double s34Pos = s34Avg / sH;
86  double s34Pos2 = s34Pos * s34Pos;
87  double s34Neg = sH / s34Avg;
88  double s34Neg2 = s34Neg * s34Neg;
89  sigSum = pow2(tmu) * sH2 * (241./1536. - 1./32. * s34Pos
90  + 9./16. * s34Pos2)
91  + pow4(tmu) * (37./512. + 9./64. * s34Pos)
92  + pow6(tmu) * (9./512. / sH2)
93  + sH2 * sH2 * (133./1536. - 7./64. * s34Pos + 7./16. * s34Pos2);
94 
95  // Anomalous coupling.
96  if (hasKappa)
97  sigSum += pow2(tmu) * sH2 * (kappam1 * (143./384. - 7./3072 * s34Neg)
98  + pow2(kappam1) * (- 1./768. * s34Neg + 185./768.)
99  + pow3(kappam1) * (- 7./3072. * s34Neg2
100  - 25./3072. * s34Neg + 67./1536.)
101  + pow4(kappam1) * (- 37./49152. * s34Neg2
102  - 25./6144. * s34Neg + 5./1536.) )
103  + pow4(tmu) * (kappam1 * 3./32.
104  + pow2(kappam1) * (7./6144. * s34Neg2 - 7./768. * s34Neg + 3./128.)
105  + pow3(kappam1) * (7./6144. * s34Neg2 - 7./1536. * s34Neg)
106  + pow4(kappam1) * (- 1./49152. * s34Neg2 + 5./6144. * s34Neg) )
107  + pow6(tmu) * pow4(kappam1) * 13./49152. / pow2(s34Avg)
108  + sH2 * sH2 * ( kappam1 * 77./384.
109  + pow2(kappam1) * (7./6144. * s34Neg2 + 1./96.* s34Neg + 39./256.)
110  + pow3(kappam1) * (7./6144. * s34Neg2 + 13./1024. * s34Neg + 61./1536.)
111  + pow4(kappam1) * (25./49152. * s34Neg2 + 5./1536. * s34Neg + 1./512.)
112  );
113 
114  // Equal probability for two possible colour flows.
115  sigSum /= pow2( (uHavg-s34Avg) * (tHavg-s34Avg) );
116  sigTS = 0.5 * sigSum;
117  sigUS = sigTS;
118  }
119 
120  // Final answer, with common factors.
121  sigma = (M_PI / sH2) * pow2(alpS) * sigSum * nCHV * openFracPair;
122 
123 }
124 
125 //--------------------------------------------------------------------------
126 
127 // Select identity, colour and anticolour.
128 
129 void Sigma2gg2qGqGbar::setIdColAcol() {
130 
131  // Flavours trivial.
132  setId( 21, 21, idNew, -idNew);
133 
134  // Two colour flow topologies.
135  double sigRand = sigSum * rndmPtr->flat();
136  if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
137  else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
138 
139 }
140 
141 //==========================================================================
142 
143 // Sigma2qqbar2qGqGbar class.
144 // Cross section for q qbar -> qG qGbar (generic quark of spin 0, 1/2 or 1).
145 
146 //--------------------------------------------------------------------------
147 
148 // Initialize process.
149 
150 void Sigma2qqbar2qGqGbar::initProc() {
151 
152  // Number of colours. Coupling kappa used for vector state.
153  nCHV = settingsPtr->mode("HiddenValley:Ngauge");
154  kappa = settingsPtr->parm("HiddenValley:kappa");
155 
156  // Secondary open width fraction.
157  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
158 
159 }
160 
161 //--------------------------------------------------------------------------
162 
163 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
164 
165 void Sigma2qqbar2qGqGbar::sigmaKin() {
166 
167  // Modified Mandelstam variables for massive kinematics with m3 = m4.
168  double delta = 0.25 * pow2(s3 - s4) / sH;
169  double s34Avg = 0.5 * (s3 + s4) - delta;
170  double tHavg = tH - delta;
171  double uHavg = uH - delta;
172  double tHQ = -0.5 * (sH - tH + uH);
173  double uHQ = -0.5 * (sH + tH - uH);
174  double tHQ2 = tHQ * tHQ;
175  double uHQ2 = uHQ * uHQ;
176 
177  // Evaluate cross section for spin 0 colour triplet.
178  if (spinSave == 0) {
179  sigSum = (1./9.) * (sH * (sH - 4. * s34Avg)
180  - pow2(uHavg - tHavg)) / sH2;
181  }
182 
183  // Evaluate cross section for spin 1/2 colour triplet.
184  else if (spinSave == 1) {
185  sigSum = (4./9.) * ((tHQ2 + uHQ2) / sH2 + 2. * s34Avg / sH);
186  }
187 
188  // Evaluate cross section for spin 1 colour triplet.
189  else {
190  double tuH34 = (tHavg + uHavg) / s34Avg;
191  sigSum = (1./9.) * (
192  pow2(1. + kappa) * sH * s34Avg * (pow2(tuH34) - 4.)
193  + (tHavg * uHavg - pow2(s34Avg)) * (8. + 2. * (1. - pow2(kappa)) * tuH34
194  + pow2(kappa) * pow2(tuH34)) ) / sH2;
195  }
196 
197  // Final answer, with common factors.
198  sigma = (M_PI / sH2) * pow2(alpS) * sigSum * nCHV * openFracPair;
199 
200 }
201 
202 //--------------------------------------------------------------------------
203 
204 // Select identity, colour and anticolour.
205 
206 void Sigma2qqbar2qGqGbar::setIdColAcol() {
207 
208  // Flavours trivial.
209  setId( id1, id2, idNew, -idNew);
210 
211  // tH defined between f and qG: must swap tHat <-> uHat if qbar q in.
212  swapTU = (id1 < 0);
213 
214  // Colour flow topologies.
215  if (id1 > 0) setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
216  else setColAcol( 0, 2, 1, 0, 1, 0, 0, 2);
217 
218 }
219 
220 //==========================================================================
221 
222 // Sigma2ffbar2fGfGbar class.
223 // Cross section for f fbar -> qG qGbar (generic quark of spin 0, 1/2 or 1)
224 // via gamma^*/Z^* s-channel exchange. Still under development!! ??
225 
226 //--------------------------------------------------------------------------
227 
228 // Initialize process.
229 
230 void Sigma2ffbar2fGfGbar::initProc() {
231 
232  // Charge and number of colours. Coupling kappa used for vector state.
233  if (settingsPtr->flag("HiddenValley:doKinMix"))
234  eQHV2 = pow2(settingsPtr->parm("HiddenValley:kinMix"));
235  else
236  eQHV2 = pow2( particleDataPtr->charge(idNew) );
237  nCHV = settingsPtr->mode("HiddenValley:Ngauge");
238  kappa = settingsPtr->parm("HiddenValley:kappa");
239 
240  // Coloured or uncoloured particle.
241  hasColour = (particleDataPtr->colType(idNew) != 0);
242  colFac = (hasColour) ? 3. : 1.;
243 
244  // Secondary open width fraction.
245  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
246 
247 }
248 
249 //--------------------------------------------------------------------------
250 
251 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
252 
253 void Sigma2ffbar2fGfGbar::sigmaKin() {
254 
255  // Modified Mandelstam variables for massive kinematics with m3 = m4.
256  double delta = 0.25 * pow2(s3 - s4) / sH;
257  double s34Avg = 0.5 * (s3 + s4) - delta;
258  double tHavg = tH - delta;
259  double uHavg = uH - delta;
260  double tHQ = -0.5 * (sH - tH + uH);
261  double uHQ = -0.5 * (sH + tH - uH);
262  double tHQ2 = tHQ * tHQ;
263  double uHQ2 = uHQ * uHQ;
264 
265  // Evaluate cross section for spin 0 colour triplet.
266  if (spinSave == 0) {
267  sigSum = 0.5 * (sH * (sH - 4. * s34Avg) - pow2(uHavg - tHavg)) / sH2;
268  }
269 
270  // Evaluate cross section for spin 1/2 colour triplet.
271  else if (spinSave == 1) {
272  sigSum = 2. * ((tHQ2 + uHQ2) / sH2 + 2. * s34Avg / sH);
273  }
274 
275  // Evaluate cross section for spin 1 colour triplet.
276  else {
277  double tuH34 = (tHavg + uHavg) / s34Avg;
278  sigSum = 0.5 * ( pow2(1. + kappa) * sH * s34Avg * (pow2(tuH34) - 4.)
279  + (tHavg * uHavg - pow2(s34Avg)) * (8. + 2. * (1. - pow2(kappa)) * tuH34
280  + pow2(kappa) * pow2(tuH34)) ) / sH2;
281  }
282 
283  // Final-state charge factors.
284  sigSum *= colFac * eQHV2 * (1. + alpS / M_PI);
285 
286  // Final answer, except for initial-state weight
287  sigma0 = (M_PI / sH2) * pow2(alpEM) * sigSum * nCHV * openFracPair;
288 
289 }
290 
291 //--------------------------------------------------------------------------
292 
293 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
294 
295 double Sigma2ffbar2fGfGbar::sigmaHat() {
296 
297  // Charge and colour factors.
298  double eNow = couplingsPtr->ef( abs(id1) );
299  double sigma = sigma0 * pow2(eNow);
300  if (abs(id1) < 9) sigma /= 3.;
301 
302  // Answer.
303  return sigma;
304 
305 }
306 
307 //--------------------------------------------------------------------------
308 
309 // Select identity, colour and anticolour.
310 
311 void Sigma2ffbar2fGfGbar::setIdColAcol() {
312 
313  // Flavours trivial.
314  setId( id1, id2, idNew, -idNew);
315 
316  // tH defined between f and qG: must swap tHat <-> uHat if fbar f in.
317  swapTU = (id1 < 0);
318 
319  // Colour flow topologies.
320  if (hasColour) {
321  if (id1 > 0 && id1 < 7) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
322  else if (id1 > -7 && id1 < 0) setColAcol( 0, 1, 1, 0, 2, 0, 0, 2);
323  else setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
324  } else {
325  if (id1 > 0 && id1 < 7) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
326  else if (id1 > -7 && id1 < 0) setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
327  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
328  }
329 
330 }
331 
332 //==========================================================================
333 
334 // Sigma1ffbar2Zv class.
335 // Cross section for f fbar -> Zv, where Zv couples both to the SM and
336 // to a hidden sector. Primitive coupling structure.
337 
338 //--------------------------------------------------------------------------
339 
340 // Initialize process.
341 
342 void Sigma1ffbar2Zv::initProc() {
343 
344  // Store Zv mass and width for propagator.
345  idZv = 4900023;
346  mRes = particleDataPtr->m0(idZv);
347  GammaRes = particleDataPtr->mWidth(idZv);
348  m2Res = mRes*mRes;
349  GamMRat = GammaRes / mRes;
350 
351  // Set pointer to particle properties and decay table.
352  particlePtr = particleDataPtr->particleDataEntryPtr(idZv);
353 
354 }
355 
356 //--------------------------------------------------------------------------
357 
358 // Evaluate sigmaHat(sHat); first step when inflavours unknown.
359 
360 void Sigma1ffbar2Zv::sigmaKin() {
361 
362  // Breit-Wigner, including some (guessed) spin factors.
363  double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
364 
365  // Outgoing width: only includes channels left open.
366  double widthOut = particlePtr->resWidthOpen(663, mH);
367 
368  // Temporary answer.
369  sigOut = sigBW * widthOut;
370 
371 }
372 
373 //--------------------------------------------------------------------------
374 
375 // Evaluate sigmaHat(sHat); second step when inflavours known.
376 
377 double Sigma1ffbar2Zv::sigmaHat() {
378 
379  // Incoming quark or lepton; for former need two 1/3 colour factors.
380  int id1Abs = abs(id1);
381  double widthIn = particlePtr->resWidthChan( mH, id1Abs, -id1Abs);
382  if (id1Abs < 6) widthIn /= 9.;
383  return widthIn * sigOut;
384 
385 }
386 
387 //--------------------------------------------------------------------------
388 
389 // Select identity, colour and anticolour.
390 
391 void Sigma1ffbar2Zv::setIdColAcol() {
392 
393  // Flavours trivial.
394  setId( id1, id2, idZv);
395 
396  // Colour flow topologies. Swap when antiquarks.
397  if (abs(id1) < 6) setColAcol( 1, 0, 0, 1, 0, 0);
398  else setColAcol( 0, 0, 0, 0, 0, 0);
399  if (id1 < 0) swapColAcol();
400 
401 }
402 
403 //--------------------------------------------------------------------------
404 
405 // Evaluate weight for decay angles.
406 
407 double Sigma1ffbar2Zv::weightDecay( Event& process, int iResBeg,
408  int iResEnd) {
409 
410  // Identity of mother of decaying resonance(s).
411  int idMother = process[process[iResBeg].mother1()].idAbs();
412 
413  // For Z' itself angular distribution as if gamma*.
414  if (iResBeg == 5 && iResEnd == 5) {
415  double mr = 4. * pow2(process[6].m()) / sH;
416  double cosThe = (process[3].p() - process[4].p())
417  * (process[7].p() - process[6].p()) / (sH * sqrtpos(1. - mr));
418  double wt = 1. + pow2(cosThe) + mr * (1. - pow2(cosThe));
419  return 0.5 * wt;
420  }
421 
422  // For top decay hand over to standard routine.
423  if (idMother == 6)
424  return weightTopDecay( process, iResBeg, iResEnd);
425 
426  // Else done.
427  return 1.;
428 
429 }
430 
431 //==========================================================================
432 
433 } // end namespace Pythia8
Definition: AgUStep.h:26