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 |
---|