1 | // SigmaCompositeness.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 | // Function definitions (not found in the header) for the |
---|
7 | // compositeness simulation classes. |
---|
8 | |
---|
9 | #include "SigmaCompositeness.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // Sigma1qg2qStar class. |
---|
16 | // Cross section for q g -> q^* (excited quark state). |
---|
17 | |
---|
18 | //-------------------------------------------------------------------------- |
---|
19 | |
---|
20 | // Initialize process. |
---|
21 | |
---|
22 | void Sigma1qg2qStar::initProc() { |
---|
23 | |
---|
24 | // Set up process properties from the chosen quark flavour. |
---|
25 | idRes = 4000000 + idq; |
---|
26 | codeSave = 4000 + idq; |
---|
27 | if (idq == 1) nameSave = "d g -> d^*"; |
---|
28 | else if (idq == 2) nameSave = "u g -> u^*"; |
---|
29 | else if (idq == 3) nameSave = "s g -> s^*"; |
---|
30 | else if (idq == 4) nameSave = "c g -> c^*"; |
---|
31 | else nameSave = "b g -> b^*"; |
---|
32 | |
---|
33 | // Store q* mass and width for propagator. |
---|
34 | mRes = particleDataPtr->m0(idRes); |
---|
35 | GammaRes = particleDataPtr->mWidth(idRes); |
---|
36 | m2Res = mRes*mRes; |
---|
37 | GamMRat = GammaRes / mRes; |
---|
38 | |
---|
39 | // Locally stored properties and couplings. |
---|
40 | Lambda = settingsPtr->parm("ExcitedFermion:Lambda"); |
---|
41 | coupFcol = settingsPtr->parm("ExcitedFermion:coupFcol"); |
---|
42 | |
---|
43 | // Set pointer to particle properties and decay table. |
---|
44 | qStarPtr = particleDataPtr->particleDataEntryPtr(idRes); |
---|
45 | |
---|
46 | } |
---|
47 | |
---|
48 | //-------------------------------------------------------------------------- |
---|
49 | |
---|
50 | // Evaluate sigmaHat(sHat), part independent of incoming flavour. |
---|
51 | |
---|
52 | void Sigma1qg2qStar::sigmaKin() { |
---|
53 | |
---|
54 | // Incoming width for correct quark. |
---|
55 | widthIn = pow3(mH) * alpS * pow2(coupFcol) / (3. * pow2(Lambda)); |
---|
56 | |
---|
57 | // Set up Breit-Wigner. |
---|
58 | sigBW = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); |
---|
59 | |
---|
60 | } |
---|
61 | |
---|
62 | //-------------------------------------------------------------------------- |
---|
63 | |
---|
64 | // Evaluate sigmaHat(sHat) for specific incoming flavours. |
---|
65 | |
---|
66 | double Sigma1qg2qStar::sigmaHat() { |
---|
67 | |
---|
68 | // Identify whether correct incoming flavours. |
---|
69 | int idqNow = (id2 == 21) ? id1 : id2; |
---|
70 | if (abs(idqNow) != idq) return 0.; |
---|
71 | |
---|
72 | // Outgoing width and total sigma. Done. |
---|
73 | return widthIn * sigBW * qStarPtr->resWidthOpen(idqNow, mH); |
---|
74 | |
---|
75 | } |
---|
76 | |
---|
77 | //-------------------------------------------------------------------------- |
---|
78 | |
---|
79 | // Select identity, colour and anticolour. |
---|
80 | |
---|
81 | void Sigma1qg2qStar::setIdColAcol() { |
---|
82 | |
---|
83 | // Flavours. |
---|
84 | int idqNow = (id2 == 21) ? id1 : id2; |
---|
85 | int idqStar = (idqNow > 0) ? idRes : -idRes; |
---|
86 | setId( id1, id2, idqStar); |
---|
87 | |
---|
88 | // Colour flow topology. |
---|
89 | if (id1 == idqNow) setColAcol( 1, 0, 2, 1, 2, 0); |
---|
90 | else setColAcol( 2, 1, 1, 0, 2, 0); |
---|
91 | if (idqNow < 0) swapColAcol(); |
---|
92 | |
---|
93 | } |
---|
94 | |
---|
95 | //-------------------------------------------------------------------------- |
---|
96 | |
---|
97 | // Evaluate weight for q* decay angle. |
---|
98 | |
---|
99 | double Sigma1qg2qStar::weightDecay( Event& process, int iResBeg, |
---|
100 | int iResEnd) { |
---|
101 | |
---|
102 | // q* should sit in entry 5. Sequential Z/W decay assumed isotropic. |
---|
103 | if (iResBeg != 5 || iResEnd != 5) return 1.; |
---|
104 | |
---|
105 | // Sign of asymmetry. |
---|
106 | int sideIn = (process[3].idAbs() < 20) ? 1 : 2; |
---|
107 | int sideOut = (process[6].idAbs() < 20) ? 1 : 2; |
---|
108 | double eps = (sideIn == sideOut) ? 1. : -1.; |
---|
109 | |
---|
110 | // Phase space factors. |
---|
111 | double mr1 = pow2(process[6].m()) / sH; |
---|
112 | double mr2 = pow2(process[7].m()) / sH; |
---|
113 | double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); |
---|
114 | |
---|
115 | // Reconstruct decay angle. Default isotropic decay. |
---|
116 | double cosThe = (process[3].p() - process[4].p()) |
---|
117 | * (process[7].p() - process[6].p()) / (sH * betaf); |
---|
118 | double wt = 1.; |
---|
119 | double wtMax = 1.; |
---|
120 | |
---|
121 | // Decay q* -> q (g/gamma) or q (Z^0/W^+-). |
---|
122 | int idBoson = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs(); |
---|
123 | if (idBoson == 21 || idBoson == 22) { |
---|
124 | wt = 1. + eps * cosThe; |
---|
125 | wtMax = 2.; |
---|
126 | } else if (idBoson == 23 || idBoson == 24) { |
---|
127 | double mrB = (sideOut == 1) ? mr2 : mr1; |
---|
128 | double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB); |
---|
129 | wt = 1. + eps * cosThe * ratB; |
---|
130 | wtMax = 1. + ratB; |
---|
131 | } |
---|
132 | |
---|
133 | // Done. |
---|
134 | return (wt / wtMax); |
---|
135 | |
---|
136 | } |
---|
137 | |
---|
138 | //========================================================================== |
---|
139 | |
---|
140 | // Sigma1lgm2lStar class. |
---|
141 | // Cross section for l gamma -> l^* (excited lepton state). |
---|
142 | |
---|
143 | //-------------------------------------------------------------------------- |
---|
144 | |
---|
145 | // Initialize process. |
---|
146 | |
---|
147 | void Sigma1lgm2lStar::initProc() { |
---|
148 | |
---|
149 | // Set up process properties from the chosen lepton flavour. |
---|
150 | idRes = 4000000 + idl; |
---|
151 | codeSave = 4000 + idl; |
---|
152 | if (idl == 11) nameSave = "e gamma -> e^*"; |
---|
153 | else if (idl == 13) nameSave = "mu gamma -> mu^*"; |
---|
154 | else nameSave = "tau gamma -> tau^*"; |
---|
155 | |
---|
156 | // Store l* mass and width for propagator. |
---|
157 | mRes = particleDataPtr->m0(idRes); |
---|
158 | GammaRes = particleDataPtr->mWidth(idRes); |
---|
159 | m2Res = mRes*mRes; |
---|
160 | GamMRat = GammaRes / mRes; |
---|
161 | |
---|
162 | // Locally stored properties and couplings. |
---|
163 | Lambda = settingsPtr->parm("ExcitedFermion:Lambda"); |
---|
164 | double coupF = settingsPtr->parm("ExcitedFermion:coupF"); |
---|
165 | double coupFp = settingsPtr->parm("ExcitedFermion:coupFprime"); |
---|
166 | coupChg = -0.5 * coupF - 0.5 * coupFp; |
---|
167 | |
---|
168 | // Set pointer to particle properties and decay table. |
---|
169 | qStarPtr = particleDataPtr->particleDataEntryPtr(idRes); |
---|
170 | |
---|
171 | } |
---|
172 | |
---|
173 | //-------------------------------------------------------------------------- |
---|
174 | |
---|
175 | // Evaluate sigmaHat(sHat), part independent of incoming flavour. |
---|
176 | |
---|
177 | void Sigma1lgm2lStar::sigmaKin() { |
---|
178 | |
---|
179 | // Incoming width for correct lepton. |
---|
180 | widthIn = pow3(mH) * alpEM * pow2(coupChg) / pow2(Lambda); |
---|
181 | |
---|
182 | // Set up Breit-Wigner. |
---|
183 | sigBW = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); |
---|
184 | |
---|
185 | } |
---|
186 | |
---|
187 | //-------------------------------------------------------------------------- |
---|
188 | |
---|
189 | // Evaluate sigmaHat(sHat) for specific incoming flavours. |
---|
190 | |
---|
191 | double Sigma1lgm2lStar::sigmaHat() { |
---|
192 | |
---|
193 | // Identify whether correct incoming flavours. |
---|
194 | int idlNow = (id2 == 22) ? id1 : id2; |
---|
195 | if (abs(idlNow) != idl) return 0.; |
---|
196 | |
---|
197 | // Outgoing width and total sigma. Done. |
---|
198 | return widthIn * sigBW * qStarPtr->resWidthOpen(idlNow, mH); |
---|
199 | |
---|
200 | } |
---|
201 | |
---|
202 | //-------------------------------------------------------------------------- |
---|
203 | |
---|
204 | // Select identity, colour and anticolour. |
---|
205 | |
---|
206 | void Sigma1lgm2lStar::setIdColAcol() { |
---|
207 | |
---|
208 | // Flavours. |
---|
209 | int idlNow = (id2 == 22) ? id1 : id2; |
---|
210 | int idlStar = (idlNow > 0) ? idRes : -idRes; |
---|
211 | setId( id1, id2, idlStar); |
---|
212 | |
---|
213 | // No colour flow. |
---|
214 | setColAcol( 0, 0, 0, 0, 0, 0); |
---|
215 | |
---|
216 | } |
---|
217 | |
---|
218 | //-------------------------------------------------------------------------- |
---|
219 | |
---|
220 | // Evaluate weight for l* decay angle. |
---|
221 | |
---|
222 | double Sigma1lgm2lStar::weightDecay( Event& process, int iResBeg, |
---|
223 | int iResEnd) { |
---|
224 | |
---|
225 | // l* should sit in entry 5. Sequential Z/W decay assumed isotropic. |
---|
226 | if (iResBeg != 5 || iResEnd != 5) return 1.; |
---|
227 | |
---|
228 | // Sign of asymmetry. |
---|
229 | int sideIn = (process[3].idAbs() < 20) ? 1 : 2; |
---|
230 | int sideOut = (process[6].idAbs() < 20) ? 1 : 2; |
---|
231 | double eps = (sideIn == sideOut) ? 1. : -1.; |
---|
232 | |
---|
233 | // Phase space factors. |
---|
234 | double mr1 = pow2(process[6].m()) / sH; |
---|
235 | double mr2 = pow2(process[7].m()) / sH; |
---|
236 | double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); |
---|
237 | |
---|
238 | // Reconstruct decay angle. Default isotropic decay. |
---|
239 | double cosThe = (process[3].p() - process[4].p()) |
---|
240 | * (process[7].p() - process[6].p()) / (sH * betaf); |
---|
241 | double wt = 1.; |
---|
242 | double wtMax = 1.; |
---|
243 | |
---|
244 | // Decay l* -> l gamma or l (Z^0/W^+-). |
---|
245 | int idBoson = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs(); |
---|
246 | if (idBoson == 22) { |
---|
247 | wt = 1. + eps * cosThe; |
---|
248 | wtMax = 2.; |
---|
249 | } else if (idBoson == 23 || idBoson == 24) { |
---|
250 | double mrB = (sideOut == 1) ? mr2 : mr1; |
---|
251 | double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB); |
---|
252 | wt = 1. + eps * cosThe * ratB; |
---|
253 | wtMax = 1. + ratB; |
---|
254 | } |
---|
255 | |
---|
256 | // Done. |
---|
257 | return (wt / wtMax); |
---|
258 | |
---|
259 | } |
---|
260 | |
---|
261 | //========================================================================== |
---|
262 | |
---|
263 | // Sigma2qq2qStarq class. |
---|
264 | // Cross section for q q' -> q^* q' (excited quark state). |
---|
265 | |
---|
266 | //-------------------------------------------------------------------------- |
---|
267 | |
---|
268 | // Initialize process. |
---|
269 | |
---|
270 | void Sigma2qq2qStarq::initProc() { |
---|
271 | |
---|
272 | // Set up process properties from the chosen quark flavour. |
---|
273 | idRes = 4000000 + idq; |
---|
274 | codeSave = 4020 + idq; |
---|
275 | if (idq == 1) nameSave = "q q -> d^* q"; |
---|
276 | else if (idq == 2) nameSave = "q q -> u^* q"; |
---|
277 | else if (idq == 3) nameSave = "q q -> s^* q"; |
---|
278 | else if (idq == 4) nameSave = "q q -> c^* q"; |
---|
279 | else nameSave = "q q -> b^* q"; |
---|
280 | |
---|
281 | // Locally stored properties and couplings. |
---|
282 | Lambda = settingsPtr->parm("ExcitedFermion:Lambda"); |
---|
283 | preFac = M_PI / pow4(Lambda); |
---|
284 | |
---|
285 | // Secondary open width fractions. |
---|
286 | openFracPos = particleDataPtr->resOpenFrac( idRes); |
---|
287 | openFracNeg = particleDataPtr->resOpenFrac(-idRes); |
---|
288 | |
---|
289 | } |
---|
290 | |
---|
291 | //-------------------------------------------------------------------------- |
---|
292 | |
---|
293 | // Evaluate sigmaHat(sHat), part independent of incoming flavour. |
---|
294 | |
---|
295 | void Sigma2qq2qStarq::sigmaKin() { |
---|
296 | |
---|
297 | // Two possible expressions, for like or unlike sign. |
---|
298 | sigmaA = preFac * (1. - s3 / sH); |
---|
299 | sigmaB = preFac * (-uH) * (sH + tH) / sH2; |
---|
300 | |
---|
301 | } |
---|
302 | |
---|
303 | //-------------------------------------------------------------------------- |
---|
304 | |
---|
305 | // Evaluate sigmaHat(sHat) for specific incoming flavours. |
---|
306 | |
---|
307 | double Sigma2qq2qStarq::sigmaHat() { |
---|
308 | |
---|
309 | // Identify different allowed incoming flavour combinations. |
---|
310 | int id1Abs = abs(id1); |
---|
311 | int id2Abs = abs(id2); |
---|
312 | double open1 = (id1 > 0) ? openFracPos : openFracNeg; |
---|
313 | double open2 = (id2 > 0) ? openFracPos : openFracNeg; |
---|
314 | double sigma = 0.; |
---|
315 | if (id1 * id2 > 0) { |
---|
316 | if (id1Abs == idq) sigma += (4./3.) * sigmaA * open1; |
---|
317 | if (id2Abs == idq) sigma += (4./3.) * sigmaA * open2; |
---|
318 | } else if (id1Abs == idq && id2 == -id1) |
---|
319 | sigma = (8./3.) * sigmaB * (open1 + open2); |
---|
320 | else if (id2 == -id1) sigma = sigmaB * (open1 + open2); |
---|
321 | else if (id1Abs == idq) sigma = sigmaB * open1; |
---|
322 | else if (id2Abs == idq) sigma = sigmaB * open2; |
---|
323 | |
---|
324 | // Done. |
---|
325 | return sigma; |
---|
326 | |
---|
327 | } |
---|
328 | |
---|
329 | //-------------------------------------------------------------------------- |
---|
330 | |
---|
331 | // Select identity, colour and anticolour. |
---|
332 | |
---|
333 | void Sigma2qq2qStarq::setIdColAcol() { |
---|
334 | |
---|
335 | // Flavours: either side may have been excited. |
---|
336 | int idAbs1 = abs(id1); |
---|
337 | int idAbs2 = abs(id2); |
---|
338 | double open1 = 0.; |
---|
339 | double open2 = 0.; |
---|
340 | if (idAbs1 == idq) open1 = (id1 > 0) ? openFracPos : openFracNeg; |
---|
341 | if (idAbs2 == idq) open2 = (id2 > 0) ? openFracPos : openFracNeg; |
---|
342 | if (open1 == 0. && open2 == 0.) { |
---|
343 | open1 = (id1 > 0) ? openFracPos : openFracNeg; |
---|
344 | open2 = (id2 > 0) ? openFracPos : openFracNeg; |
---|
345 | } |
---|
346 | bool excite1 = (open1 > 0.); |
---|
347 | if (open1 > 0. && open2 > 0.) excite1 |
---|
348 | = (rndmPtr->flat() * (open1 + open2) < open1); |
---|
349 | |
---|
350 | // Always excited quark in slot 3 so colour flow flipped or not. |
---|
351 | if (excite1) { |
---|
352 | id3 = (id1 > 0) ? idRes : -idRes; |
---|
353 | id4 = id2; |
---|
354 | // Special case for s-channel like production. |
---|
355 | if ((idAbs1 == idAbs2) && (id1 * id2 < 0)) { |
---|
356 | id4 = (id3 > 0) ? -idq : idq; |
---|
357 | } |
---|
358 | if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0); |
---|
359 | else setColAcol( 1, 0, 0, 2, 1, 0, 0, 2); |
---|
360 | if (id1 < 0) swapColAcol(); |
---|
361 | } else { |
---|
362 | id3 = (id2 > 0) ? idRes : -idRes; |
---|
363 | id4 = id1; |
---|
364 | // Special case for s-channel like production. |
---|
365 | if ((idAbs1 == idAbs2) && (id1 * id2 < 0)) { |
---|
366 | id4 = (id3 > 0) ? -idq : idq; |
---|
367 | } |
---|
368 | swapTU = true; |
---|
369 | if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0); |
---|
370 | else setColAcol( 1, 0, 0, 2, 0, 2, 1, 0); |
---|
371 | if (id1 < 0) swapColAcol(); |
---|
372 | } |
---|
373 | setId( id1, id2, id3, id4); |
---|
374 | |
---|
375 | } |
---|
376 | |
---|
377 | //-------------------------------------------------------------------------- |
---|
378 | |
---|
379 | // Evaluate weight for q* decay angle. |
---|
380 | // SA: Angles dist. for decay q* -> q V, based on Eq. 1.7 |
---|
381 | // in CERN Yellow Reports 90-10 vol.2, p. 1014 to 1021. |
---|
382 | |
---|
383 | double Sigma2qq2qStarq::weightDecay( Event& process, int iResBeg, |
---|
384 | int iResEnd) { |
---|
385 | |
---|
386 | // q* should sit in entry 5. Sequential Z/W decay assumed isotropic. |
---|
387 | if (iResBeg != 5 && iResEnd != 5) return 1.; |
---|
388 | |
---|
389 | // Phase space factors. |
---|
390 | double mr1 = pow2(process[7].m() / process[5].m()); |
---|
391 | double mr2 = pow2(process[8].m() / process[5].m()); |
---|
392 | |
---|
393 | // Reconstruct decay angle in q* CoM frame. |
---|
394 | int idAbs3 = process[7].idAbs(); |
---|
395 | Vec4 pQStarCom = (idAbs3 < 20) ? process[7].p() : process[8].p(); |
---|
396 | pQStarCom.bstback(process[5].p()); |
---|
397 | double cosThe = costheta(pQStarCom, process[5].p()); |
---|
398 | double wt = 1.; |
---|
399 | |
---|
400 | // Decay q* -> q (g/gamma) or q (Z^0/W^+-). |
---|
401 | int idBoson = (idAbs3 < 20) ? process[8].idAbs() : process[7].idAbs(); |
---|
402 | if (idBoson == 21 || idBoson == 22) { |
---|
403 | wt = 0.5 * (1. + cosThe); |
---|
404 | } else if (idBoson == 23 || idBoson == 24) { |
---|
405 | double mrB = (idAbs3 < 20) ? mr2 : mr1; |
---|
406 | double kTrm = 0.5 * (mrB * (1. - cosThe)); |
---|
407 | wt = (1. + cosThe + kTrm) / (2 + mrB); |
---|
408 | } |
---|
409 | |
---|
410 | // Done. |
---|
411 | return wt; |
---|
412 | } |
---|
413 | |
---|
414 | //========================================================================== |
---|
415 | |
---|
416 | // Sigma2qqbar2lStarlbar class. |
---|
417 | // Cross section for q qbar -> l^* lbar (excited lepton state). |
---|
418 | |
---|
419 | //-------------------------------------------------------------------------- |
---|
420 | |
---|
421 | // Initialize process. |
---|
422 | |
---|
423 | void Sigma2qqbar2lStarlbar::initProc() { |
---|
424 | |
---|
425 | // Set up process properties from the chosen lepton flavour. |
---|
426 | idRes = 4000000 + idl; |
---|
427 | codeSave = 4020 + idl; |
---|
428 | if (idl == 11) nameSave = "q qbar -> e^*+- e^-+"; |
---|
429 | else if (idl == 12) nameSave = "q qbar -> nu_e^* nu_ebar"; |
---|
430 | else if (idl == 13) nameSave = "q qbar -> mu^*+- mu^-+"; |
---|
431 | else if (idl == 14) nameSave = "q qbar -> nu_mu^* nu_mubar"; |
---|
432 | else if (idl == 15) nameSave = "q qbar -> tau^*+- tau^-+"; |
---|
433 | else nameSave = "q qbar -> nu_tau^* nu_taubar"; |
---|
434 | |
---|
435 | // Secondary open width fractions. |
---|
436 | openFracPos = particleDataPtr->resOpenFrac( idRes); |
---|
437 | openFracNeg = particleDataPtr->resOpenFrac(-idRes); |
---|
438 | |
---|
439 | // Locally stored properties and couplings. |
---|
440 | Lambda = settingsPtr->parm("ExcitedFermion:Lambda"); |
---|
441 | preFac = (M_PI / pow4(Lambda)) * (openFracPos + openFracNeg) / 3.; |
---|
442 | |
---|
443 | } |
---|
444 | |
---|
445 | //-------------------------------------------------------------------------- |
---|
446 | |
---|
447 | // Evaluate sigmaHat(sHat), part independent of incoming flavour. |
---|
448 | |
---|
449 | void Sigma2qqbar2lStarlbar::sigmaKin() { |
---|
450 | |
---|
451 | // Only one possible expressions |
---|
452 | sigma = preFac * (-uH) * (sH + tH) / sH2; |
---|
453 | |
---|
454 | } |
---|
455 | |
---|
456 | //-------------------------------------------------------------------------- |
---|
457 | |
---|
458 | // Select identity, colour and anticolour. |
---|
459 | |
---|
460 | void Sigma2qqbar2lStarlbar::setIdColAcol() { |
---|
461 | |
---|
462 | // Flavours: either lepton or antilepton may be excited. |
---|
463 | if (rndmPtr->flat() * (openFracPos + openFracNeg) < openFracPos) { |
---|
464 | setId( id1, id2, idRes, -idl); |
---|
465 | if (id1 < 0) swapTU = true; |
---|
466 | } else { |
---|
467 | setId( id1, id2, -idRes, idl); |
---|
468 | if (id1 > 0) swapTU = true; |
---|
469 | } |
---|
470 | |
---|
471 | // Colour flow trivial. |
---|
472 | if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0); |
---|
473 | else setColAcol( 0, 1, 1, 0, 0, 0, 0, 0); |
---|
474 | |
---|
475 | } |
---|
476 | |
---|
477 | //-------------------------------------------------------------------------- |
---|
478 | |
---|
479 | // Evaluate weight for l* decay angle. |
---|
480 | // SA: Angles dist. for decay l* -> l V, based on Eq. 1.7 |
---|
481 | // in CERN Yellow Reports 90-10 vol.2, p. 1014 to 1021. |
---|
482 | |
---|
483 | double Sigma2qqbar2lStarlbar::weightDecay( Event& process, int iResBeg, |
---|
484 | int iResEnd) { |
---|
485 | |
---|
486 | // l* should sit in entry 5. Sequential Z/W decay assumed isotropic. |
---|
487 | if (iResBeg != 5 && iResEnd != 5) return 1.; |
---|
488 | |
---|
489 | // Phase space factors. |
---|
490 | double mr1 = pow2(process[7].m() / process[5].m()); |
---|
491 | double mr2 = pow2(process[8].m() / process[5].m()); |
---|
492 | |
---|
493 | // Reconstruct decay angle in l* CoM frame. |
---|
494 | int idAbs3 = process[7].idAbs(); |
---|
495 | Vec4 pLStarCom = (idAbs3 < 20) ? process[7].p() : process[8].p(); |
---|
496 | pLStarCom.bstback(process[5].p()); |
---|
497 | double cosThe = costheta(pLStarCom, process[5].p()); |
---|
498 | double wt = 1.; |
---|
499 | |
---|
500 | // Decay, l* -> l + gamma/Z^0/W^+-). |
---|
501 | int idBoson = (idAbs3 < 20) ? process[8].idAbs() : process[7].idAbs(); |
---|
502 | if (idBoson == 22) { |
---|
503 | wt = 0.5 * (1. + cosThe); |
---|
504 | } else if (idBoson == 23 || idBoson == 24) { |
---|
505 | double mrB = (idAbs3 < 20) ? mr2 : mr1; |
---|
506 | double kTrm = 0.5 * (mrB * (1. - cosThe)); |
---|
507 | wt = (1. + cosThe + kTrm) / (2 + mrB); |
---|
508 | } |
---|
509 | |
---|
510 | // Done. |
---|
511 | return wt; |
---|
512 | } |
---|
513 | |
---|
514 | //========================================================================== |
---|
515 | |
---|
516 | // Sigma2QCqq2qq class. |
---|
517 | // Cross section for q q -> q q (quark contact interactions). |
---|
518 | |
---|
519 | //-------------------------------------------------------------------------- |
---|
520 | |
---|
521 | // Initialize process. |
---|
522 | |
---|
523 | void Sigma2QCqq2qq::initProc() { |
---|
524 | |
---|
525 | qCLambda2 = settingsPtr->parm("ContactInteractions:Lambda"); |
---|
526 | qCetaLL = settingsPtr->mode("ContactInteractions:etaLL"); |
---|
527 | qCetaRR = settingsPtr->mode("ContactInteractions:etaRR"); |
---|
528 | qCetaLR = settingsPtr->mode("ContactInteractions:etaLR"); |
---|
529 | qCLambda2 *= qCLambda2; |
---|
530 | |
---|
531 | } |
---|
532 | |
---|
533 | //-------------------------------------------------------------------------- |
---|
534 | |
---|
535 | // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. |
---|
536 | |
---|
537 | void Sigma2QCqq2qq::sigmaKin() { |
---|
538 | |
---|
539 | // Calculate kinematics dependence for different terms. |
---|
540 | sigT = (4./9.) * (sH2 + uH2) / tH2; |
---|
541 | sigU = (4./9.) * (sH2 + tH2) / uH2; |
---|
542 | sigTU = - (8./27.) * sH2 / (tH * uH); |
---|
543 | sigST = - (8./27.) * uH2 / (sH * tH); |
---|
544 | |
---|
545 | sigQCSTU = sH2 * (1 / tH + 1 / uH); |
---|
546 | sigQCUTS = uH2 * (1 / tH + 1 / sH); |
---|
547 | |
---|
548 | } |
---|
549 | |
---|
550 | //-------------------------------------------------------------------------- |
---|
551 | |
---|
552 | // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. |
---|
553 | |
---|
554 | double Sigma2QCqq2qq::sigmaHat() { |
---|
555 | |
---|
556 | // Terms from QC contact interactions. |
---|
557 | double sigQCLL = 0; |
---|
558 | double sigQCRR = 0; |
---|
559 | double sigQCLR = 0; |
---|
560 | |
---|
561 | // Combine cross section terms; factor 1/2 when identical quarks. |
---|
562 | // q q -> q q |
---|
563 | if (id2 == id1) { |
---|
564 | |
---|
565 | // SM terms. |
---|
566 | sigSum = 0.5 * (sigT + sigU + sigTU); |
---|
567 | |
---|
568 | // Contact terms. |
---|
569 | sigQCLL = (8./9.) * alpS * (qCetaLL/qCLambda2) * sigQCSTU |
---|
570 | + (8./3.) * pow2(qCetaLL/qCLambda2) * sH2; |
---|
571 | sigQCRR = (8./9.) * alpS * (qCetaRR/qCLambda2) * sigQCSTU |
---|
572 | + (8./3.) * pow2(qCetaRR/qCLambda2) * sH2; |
---|
573 | sigQCLR = 2. * (uH2 + tH2) * pow2(qCetaLR/qCLambda2); |
---|
574 | |
---|
575 | sigQCLL /= 2; |
---|
576 | sigQCRR /= 2; |
---|
577 | sigQCLR /= 2; |
---|
578 | |
---|
579 | // q qbar -> q qbar, without pure s-channel term. |
---|
580 | } else if (id2 == -id1) { |
---|
581 | |
---|
582 | // SM terms. |
---|
583 | sigSum = sigT + sigST; |
---|
584 | |
---|
585 | // Contact terms, minus the terms included in qqbar2qqbar. |
---|
586 | sigQCLL = (8./9.) * alpS * (qCetaLL/qCLambda2) * sigQCUTS |
---|
587 | + (5./3.) * pow2(qCetaLL/qCLambda2) * uH2; |
---|
588 | sigQCRR = (8./9.) * alpS * (qCetaRR/qCLambda2) * sigQCUTS |
---|
589 | + (5./3.) * pow2(qCetaRR/qCLambda2) * uH2; |
---|
590 | sigQCLR = 2. * sH2 * pow2(qCetaLR/qCLambda2); |
---|
591 | |
---|
592 | // q q' -> q q' or q qbar' -> q qbar' |
---|
593 | } else { |
---|
594 | |
---|
595 | // SM terms. |
---|
596 | sigSum = sigT; |
---|
597 | |
---|
598 | // Contact terms. |
---|
599 | if (id1 * id2 > 0) { |
---|
600 | sigQCLL = pow2(qCetaLL/qCLambda2) * sH2; |
---|
601 | sigQCRR = pow2(qCetaRR/qCLambda2) * sH2; |
---|
602 | sigQCLR = 2 * pow2(qCetaLR/qCLambda2) * uH2; |
---|
603 | } else { |
---|
604 | sigQCLL = pow2(qCetaLL/qCLambda2) * uH2; |
---|
605 | sigQCRR = pow2(qCetaRR/qCLambda2) * uH2; |
---|
606 | sigQCLR = 2 * pow2(qCetaLR/qCLambda2) * sH2; |
---|
607 | } |
---|
608 | } |
---|
609 | |
---|
610 | // Answer. |
---|
611 | double sigma = (M_PI/sH2) * ( pow2(alpS) * sigSum |
---|
612 | + sigQCLL + sigQCRR + sigQCLR ); |
---|
613 | return sigma; |
---|
614 | |
---|
615 | } |
---|
616 | |
---|
617 | //-------------------------------------------------------------------------- |
---|
618 | |
---|
619 | // Select identity, colour and anticolour. |
---|
620 | |
---|
621 | void Sigma2QCqq2qq::setIdColAcol() { |
---|
622 | |
---|
623 | // Outgoing = incoming flavours. |
---|
624 | setId( id1, id2, id1, id2); |
---|
625 | |
---|
626 | // Colour flow topologies. Swap when antiquarks. |
---|
627 | if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0); |
---|
628 | else setColAcol( 1, 0, 0, 1, 2, 0, 0, 2); |
---|
629 | if (id2 == id1 && (sigT + sigU) * rndmPtr->flat() > sigT) |
---|
630 | setColAcol( 1, 0, 2, 0, 1, 0, 2, 0); |
---|
631 | if (id1 < 0) swapColAcol(); |
---|
632 | |
---|
633 | } |
---|
634 | |
---|
635 | //========================================================================== |
---|
636 | |
---|
637 | // Sigma2QCqqbar2qqbar class. |
---|
638 | // Cross section for q qbar -> q' qbar' (quark contact interactions). |
---|
639 | |
---|
640 | //-------------------------------------------------------------------------- |
---|
641 | |
---|
642 | // Initialize process. |
---|
643 | |
---|
644 | void Sigma2QCqqbar2qqbar::initProc() { |
---|
645 | |
---|
646 | qCnQuarkNew = settingsPtr->mode("ContactInteractions:nQuarkNew"); |
---|
647 | qCLambda2 = settingsPtr->parm("ContactInteractions:Lambda"); |
---|
648 | qCetaLL = settingsPtr->mode("ContactInteractions:etaLL"); |
---|
649 | qCetaRR = settingsPtr->mode("ContactInteractions:etaRR"); |
---|
650 | qCetaLR = settingsPtr->mode("ContactInteractions:etaLR"); |
---|
651 | qCLambda2 *= qCLambda2; |
---|
652 | |
---|
653 | } |
---|
654 | |
---|
655 | //-------------------------------------------------------------------------- |
---|
656 | |
---|
657 | // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence. |
---|
658 | |
---|
659 | void Sigma2QCqqbar2qqbar::sigmaKin() { |
---|
660 | |
---|
661 | // Pick new flavour. |
---|
662 | idNew = 1 + int( qCnQuarkNew * rndmPtr->flat() ); |
---|
663 | mNew = particleDataPtr->m0(idNew); |
---|
664 | m2New = mNew*mNew; |
---|
665 | |
---|
666 | // Calculate kinematics dependence. |
---|
667 | double sigQC = 0.; |
---|
668 | sigS = 0.; |
---|
669 | if (sH > 4. * m2New) { |
---|
670 | sigS = (4./9.) * (tH2 + uH2) / sH2; |
---|
671 | sigQC = pow2(qCetaLL/qCLambda2) * uH2 |
---|
672 | + pow2(qCetaRR/qCLambda2) * uH2 |
---|
673 | + 2 * pow2(qCetaLR/qCLambda2) * tH2; |
---|
674 | } |
---|
675 | |
---|
676 | // Answer is proportional to number of outgoing flavours. |
---|
677 | sigma = (M_PI / sH2) * qCnQuarkNew * ( pow2(alpS) * sigS + sigQC); |
---|
678 | |
---|
679 | } |
---|
680 | |
---|
681 | //-------------------------------------------------------------------------- |
---|
682 | |
---|
683 | // Select identity, colour and anticolour. |
---|
684 | |
---|
685 | void Sigma2QCqqbar2qqbar::setIdColAcol() { |
---|
686 | |
---|
687 | // Set outgoing flavours ones. |
---|
688 | id3 = (id1 > 0) ? idNew : -idNew; |
---|
689 | setId( id1, id2, id3, -id3); |
---|
690 | |
---|
691 | // Colour flow topologies. Swap when antiquarks. |
---|
692 | setColAcol( 1, 0, 0, 2, 1, 0, 0, 2); |
---|
693 | if (id1 < 0) swapColAcol(); |
---|
694 | |
---|
695 | } |
---|
696 | |
---|
697 | //========================================================================== |
---|
698 | |
---|
699 | // Sigma2QCffbar2llbar class. |
---|
700 | // Cross section for f fbar -> l lbar (contact interactions). |
---|
701 | |
---|
702 | //-------------------------------------------------------------------------- |
---|
703 | |
---|
704 | // Initialize process. |
---|
705 | |
---|
706 | void Sigma2QCffbar2llbar::initProc() { |
---|
707 | |
---|
708 | qCLambda2 = settingsPtr->parm("ContactInteractions:Lambda"); |
---|
709 | qCetaLL = settingsPtr->mode("ContactInteractions:etaLL"); |
---|
710 | qCetaRR = settingsPtr->mode("ContactInteractions:etaRR"); |
---|
711 | qCetaLR = settingsPtr->mode("ContactInteractions:etaLR"); |
---|
712 | qCLambda2 *= qCLambda2; |
---|
713 | |
---|
714 | // Process name. |
---|
715 | if (idNew == 11) nameNew = "f fbar -> (QC) -> e- e+"; |
---|
716 | if (idNew == 13) nameNew = "f fbar -> (QC) -> mu- mu+"; |
---|
717 | if (idNew == 15) nameNew = "f fbar -> (QC) -> tau- tau+"; |
---|
718 | |
---|
719 | // Kinematics. |
---|
720 | qCmNew = particleDataPtr->m0(idNew); |
---|
721 | qCmNew2 = qCmNew * qCmNew; |
---|
722 | qCmZ = particleDataPtr->m0(23); |
---|
723 | qCmZ2 = qCmZ * qCmZ; |
---|
724 | qCGZ = particleDataPtr->mWidth(23); |
---|
725 | qCGZ2 = qCGZ * qCGZ; |
---|
726 | |
---|
727 | } |
---|
728 | |
---|
729 | //-------------------------------------------------------------------------- |
---|
730 | |
---|
731 | // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence. |
---|
732 | |
---|
733 | void Sigma2QCffbar2llbar::sigmaKin() { |
---|
734 | |
---|
735 | qCPropGm = 1./sH; |
---|
736 | double denomPropZ = pow2(sH - qCmZ2) + qCmZ2 * qCGZ2; |
---|
737 | qCrePropZ = (sH - qCmZ2) / denomPropZ; |
---|
738 | qCimPropZ = -qCmZ * qCGZ / denomPropZ; |
---|
739 | |
---|
740 | sigma0 = 0.; |
---|
741 | if (sH > 4. * qCmNew2) sigma0 = 1./(16. * M_PI * sH2); |
---|
742 | |
---|
743 | } |
---|
744 | |
---|
745 | //-------------------------------------------------------------------------- |
---|
746 | |
---|
747 | // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence. |
---|
748 | |
---|
749 | double Sigma2QCffbar2llbar::sigmaHat() { |
---|
750 | |
---|
751 | // Incoming fermion flavor. |
---|
752 | int idAbs = abs(id1); |
---|
753 | |
---|
754 | // Couplings and constants. |
---|
755 | double tmPe2QfQl = 4. * M_PI * alpEM * couplingsPtr->ef(idAbs) |
---|
756 | * couplingsPtr->ef(idNew); |
---|
757 | double tmPgvf = 0.25 * couplingsPtr->vf(idAbs); |
---|
758 | double tmPgaf = 0.25 * couplingsPtr->af(idAbs); |
---|
759 | double tmPgLf = tmPgvf + tmPgaf; |
---|
760 | double tmPgRf = tmPgvf - tmPgaf; |
---|
761 | double tmPgvl = 0.25 * couplingsPtr->vf(idNew); |
---|
762 | double tmPgal = 0.25 * couplingsPtr->af(idNew); |
---|
763 | double tmPgLl = tmPgvl + tmPgal; |
---|
764 | double tmPgRl = tmPgvl - tmPgal; |
---|
765 | double tmPe2s2c2 = 4. * M_PI * alpEM |
---|
766 | / (couplingsPtr->sin2thetaW() * couplingsPtr->cos2thetaW()); |
---|
767 | |
---|
768 | // Complex amplitudes. |
---|
769 | complex I(0., 1.); |
---|
770 | complex meLL(0., 0.); |
---|
771 | complex meRR(0., 0.); |
---|
772 | complex meLR(0., 0.); |
---|
773 | complex meRL(0., 0.); |
---|
774 | |
---|
775 | // Amplitudes, M = gamma + Z + CI. |
---|
776 | meLL = tmPe2QfQl * qCPropGm |
---|
777 | + tmPe2s2c2 * tmPgLf * tmPgLl * (qCrePropZ + I * qCimPropZ) |
---|
778 | + 2. * M_PI * qCetaLL / qCLambda2; |
---|
779 | meRR = tmPe2QfQl * qCPropGm |
---|
780 | + tmPe2s2c2 * tmPgRf * tmPgRl * (qCrePropZ + I * qCimPropZ) |
---|
781 | + 2. * M_PI * qCetaRR / qCLambda2; |
---|
782 | meLR = tmPe2QfQl * qCPropGm |
---|
783 | + tmPe2s2c2 * tmPgLf * tmPgRl * (qCrePropZ + I * qCimPropZ) |
---|
784 | + 2. * M_PI * qCetaLR / qCLambda2; |
---|
785 | meRL = tmPe2QfQl * qCPropGm |
---|
786 | + tmPe2s2c2 * tmPgRf * tmPgLl * (qCrePropZ + I * qCimPropZ) |
---|
787 | + 2. * M_PI * qCetaLR / qCLambda2; |
---|
788 | |
---|
789 | double sigma = sigma0 * uH2 * real(meLL*conj(meLL)); |
---|
790 | sigma += sigma0 * uH2 * real(meRR*conj(meRR)); |
---|
791 | sigma += sigma0 * tH2 * real(meLR*conj(meLR)); |
---|
792 | sigma += sigma0 * tH2 * real(meRL*conj(meRL)); |
---|
793 | |
---|
794 | // If f fbar are quarks. |
---|
795 | if (idAbs < 9) sigma /= 3.; |
---|
796 | |
---|
797 | return sigma; |
---|
798 | } |
---|
799 | |
---|
800 | //-------------------------------------------------------------------------- |
---|
801 | |
---|
802 | // Select identity, colour and anticolour. |
---|
803 | |
---|
804 | void Sigma2QCffbar2llbar::setIdColAcol() { |
---|
805 | |
---|
806 | // Flavours trivial. |
---|
807 | setId(id1, id2, idNew, -idNew); |
---|
808 | |
---|
809 | // tH defined between f and f': must swap tHat <-> uHat if id1 is fbar. |
---|
810 | swapTU = (id2 > 0); |
---|
811 | |
---|
812 | // Colour flow topologies. Swap when antiquarks. |
---|
813 | if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0); |
---|
814 | else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0); |
---|
815 | if (id1 < 0) swapColAcol(); |
---|
816 | |
---|
817 | } |
---|
818 | |
---|
819 | //========================================================================== |
---|
820 | |
---|
821 | } // end namespace Pythia8 |
---|