Line data Source code
1 : // StandardModel.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 AlphaStrong class.
7 :
8 : #include "Pythia8/StandardModel.h"
9 :
10 : namespace Pythia8 {
11 :
12 : //==========================================================================
13 :
14 : // The AlphaStrong class.
15 :
16 : //--------------------------------------------------------------------------
17 :
18 : // Constants: could be changed here if desired, but normally should not.
19 : // These are of technical nature, as described for each.
20 :
21 : // Number of iterations to determine Lambda from given alpha_s.
22 : const int AlphaStrong::NITER = 10;
23 :
24 : // MZ normalization scale
25 : const double AlphaStrong::MZ = 91.188;
26 :
27 : // Always evaluate running alpha_s above Lambda3 to avoid disaster.
28 : // Safety margin picked to freeze roughly for alpha_s = 10.
29 : const double AlphaStrong::SAFETYMARGIN1 = 1.07;
30 : const double AlphaStrong::SAFETYMARGIN2 = 1.33;
31 :
32 : // CMW factor for 3, 4, 5, and 6 flavours.
33 : const double AlphaStrong::FACCMW3 = 1.661;
34 : const double AlphaStrong::FACCMW4 = 1.618;
35 : const double AlphaStrong::FACCMW5 = 1.569;
36 : const double AlphaStrong::FACCMW6 = 1.513;
37 :
38 : //--------------------------------------------------------------------------
39 :
40 : // Initialize alpha_strong calculation by finding Lambda values etc.
41 :
42 : void AlphaStrong::init( double valueIn, int orderIn, int nfmaxIn,
43 : bool useCMWIn) {
44 :
45 : // Set default mass thresholds if not already done
46 0 : if (mt <= 1.) setThresholds(1.5, 4.8, 171.0);
47 :
48 : // Order of alpha_s evaluation.Default values.
49 0 : valueRef = valueIn;
50 0 : order = max( 0, min( 2, orderIn ) );
51 0 : nfmax = max(5,min(6,nfmaxIn));
52 0 : useCMW = useCMWIn;
53 0 : lastCallToFull = false;
54 0 : Lambda3Save = Lambda4Save = Lambda5Save = Lambda6Save = scale2Min = 0.;
55 :
56 : // Fix alpha_s.
57 0 : if (order == 0) {
58 :
59 : // First order alpha_s: match at flavour thresholds.
60 0 : } else if (order == 1) {
61 0 : Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueRef) );
62 0 : Lambda6Save = Lambda5Save * pow(Lambda5Save/mt, 2./21.);
63 0 : Lambda4Save = Lambda5Save * pow(mb/Lambda5Save, 2./25.);
64 0 : Lambda3Save = Lambda4Save * pow(mc/Lambda4Save, 2./27.);
65 :
66 : // Second order alpha_s: iterative match at flavour thresholds.
67 0 : } else {
68 : // The one-loop coefficients: b1 / b0^2
69 : double b16 = 234. / 441.;
70 : double b15 = 348. / 529.;
71 : double b14 = 462. / 625.;
72 : double b13 = 576. / 729.;
73 : // The two-loop coefficients: b2 * b0 / b1^2
74 : double b26 = -36855. / 109512.;
75 : double b25 = 224687. / 242208.;
76 : double b24 = 548575. / 426888.;
77 : double b23 = 938709. / 663552.;
78 :
79 : double logScale, loglogScale, correction, valueIter;
80 : // Find Lambda_5 at m_Z, starting from one-loop value
81 0 : Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueRef) );
82 0 : for (int iter = 0; iter < NITER; ++iter) {
83 0 : logScale = 2. * log(MZ/Lambda5Save);
84 0 : loglogScale = log(logScale);
85 0 : correction = 1. - b15 * loglogScale / logScale
86 0 : + pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25);
87 0 : valueIter = valueRef / correction;
88 0 : Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueIter) );
89 : }
90 :
91 : // Find Lambda_6 at m_t, by requiring alphaS(nF=6,m_t) = alphaS(nF=5,m_t)
92 0 : double logScaleT = 2. * log(mt/Lambda5Save);
93 0 : double loglogScaleT = log(logScaleT);
94 0 : double valueT = 12. * M_PI / (23. * logScaleT)
95 0 : * (1. - b15 * loglogScaleT / logScaleT
96 0 : + pow2(b15 / logScaleT) * (pow2(loglogScaleT - 0.5) + b25 - 1.25) );
97 0 : Lambda6Save = Lambda5Save;
98 0 : for (int iter = 0; iter < NITER; ++iter) {
99 0 : logScale = 2. * log(mt/Lambda6Save);
100 0 : loglogScale = log(logScale);
101 0 : correction = 1. - b16 * loglogScale / logScale
102 0 : + pow2(b16 / logScale) * (pow2(loglogScale - 0.5) + b26 - 1.25);
103 0 : valueIter = valueT / correction;
104 0 : Lambda6Save = mt * exp( -6. * M_PI / (21. * valueIter) );
105 : }
106 :
107 : // Find Lambda_4 at m_b, by requiring alphaS(nF=4,m_b) = alphaS(nF=5,m_b)
108 0 : double logScaleB = 2. * log(mb/Lambda5Save);
109 0 : double loglogScaleB = log(logScaleB);
110 0 : double valueB = 12. * M_PI / (23. * logScaleB)
111 0 : * (1. - b15 * loglogScaleB / logScaleB
112 0 : + pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) );
113 0 : Lambda4Save = Lambda5Save;
114 0 : for (int iter = 0; iter < NITER; ++iter) {
115 0 : logScale = 2. * log(mb/Lambda4Save);
116 0 : loglogScale = log(logScale);
117 0 : correction = 1. - b14 * loglogScale / logScale
118 0 : + pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25);
119 0 : valueIter = valueB / correction;
120 0 : Lambda4Save = mb * exp( -6. * M_PI / (25. * valueIter) );
121 : }
122 : // Find Lambda_3 at m_c, by requiring alphaS(nF=3,m_c) = alphaS(nF=4,m_c)
123 0 : double logScaleC = 2. * log(mc/Lambda4Save);
124 0 : double loglogScaleC = log(logScaleC);
125 0 : double valueC = 12. * M_PI / (25. * logScaleC)
126 0 : * (1. - b14 * loglogScaleC / logScaleC
127 0 : + pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) );
128 0 : Lambda3Save = Lambda4Save;
129 0 : for (int iter = 0; iter < NITER; ++iter) {
130 0 : logScale = 2. * log(mc/Lambda3Save);
131 0 : loglogScale = log(logScale);
132 0 : correction = 1. - b13 * loglogScale / logScale
133 0 : + pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25);
134 0 : valueIter = valueC / correction;
135 0 : Lambda3Save = mc * exp( -6. * M_PI / (27. * valueIter) );
136 : }
137 : }
138 :
139 : // Optionally rescale Lambda values by CMW factor.
140 0 : if (useCMW) {
141 0 : Lambda3Save *= FACCMW3;
142 0 : Lambda4Save *= FACCMW4;
143 0 : Lambda5Save *= FACCMW5;
144 0 : Lambda6Save *= FACCMW6;
145 0 : }
146 :
147 : // Impose SAFETYMARGINs to prevent getting too close to LambdaQCD.
148 0 : if (order == 1) scale2Min = pow2(SAFETYMARGIN1 * Lambda3Save);
149 0 : else if (order == 2) scale2Min = pow2(SAFETYMARGIN2 * Lambda3Save);
150 :
151 : // Save squares of mass and Lambda values as well.
152 0 : Lambda3Save2 = pow2(Lambda3Save);
153 0 : Lambda4Save2 = pow2(Lambda4Save);
154 0 : Lambda5Save2 = pow2(Lambda5Save);
155 0 : Lambda6Save2 = pow2(Lambda6Save);
156 0 : mc2 = pow2(mc);
157 0 : mb2 = pow2(mb);
158 0 : mt2 = pow2(mt);
159 0 : valueNow = valueIn;
160 0 : scale2Now = MZ * MZ;
161 0 : isInit = true;
162 :
163 0 : }
164 :
165 : //--------------------------------------------------------------------------
166 :
167 : // Calculate alpha_s value.
168 :
169 : double AlphaStrong::alphaS( double scale2) {
170 :
171 : // Check for initialization and ensure minimal scale2 value.
172 0 : if (!isInit) return 0.;
173 0 : if (scale2 < scale2Min) scale2 = scale2Min;
174 :
175 : // If equal to old scale then same answer.
176 0 : if (scale2 == scale2Now && (order < 2 || lastCallToFull)) return valueNow;
177 0 : scale2Now = scale2;
178 0 : lastCallToFull = true;
179 :
180 : // Fix alpha_s.
181 0 : if (order == 0) {
182 0 : valueNow = valueRef;
183 :
184 : // First order alpha_s: differs by mass region.
185 0 : } else if (order == 1) {
186 0 : if (scale2 > mt2 && nfmax >= 6)
187 0 : valueNow = 12. * M_PI / (21. * log(scale2/Lambda6Save2));
188 0 : else if (scale2 > mb2)
189 0 : valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
190 0 : else if (scale2 > mc2)
191 0 : valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
192 0 : else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
193 :
194 : // Second order alpha_s: differs by mass region.
195 : } else {
196 : double Lambda2, b0, b1, b2;
197 0 : if (scale2 > mt2 && nfmax >= 6) {
198 0 : Lambda2 = Lambda6Save2;
199 : b0 = 21.;
200 : b1 = 234. / 441.;
201 : b2 = -36855. / 109512.;
202 0 : } else if (scale2 > mb2) {
203 0 : Lambda2 = Lambda5Save2;
204 : b0 = 23.;
205 : b1 = 348. / 529.;
206 : b2 = 224687. / 242208.;
207 0 : } else if (scale2 > mc2) {
208 0 : Lambda2 = Lambda4Save2;
209 : b0 = 25.;
210 : b1 = 462. / 625.;
211 : b2 = 548575. / 426888.;
212 0 : } else {
213 0 : Lambda2 = Lambda3Save2;
214 : b0 = 27.;
215 : b1 = 64. / 81.;
216 : b2 = 938709. / 663552.;
217 : }
218 0 : double logScale = log(scale2/Lambda2);
219 0 : double loglogScale = log(logScale);
220 0 : valueNow = 12. * M_PI / (b0 * logScale)
221 0 : * ( 1. - b1 * loglogScale / logScale
222 0 : + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
223 : }
224 :
225 : // Done.
226 0 : return valueNow;
227 :
228 0 : }
229 :
230 : //--------------------------------------------------------------------------
231 :
232 : // Calculate alpha_s value, but only use up to first-order piece.
233 : // (To be combined with alphaS2OrdCorr.)
234 :
235 : double AlphaStrong::alphaS1Ord( double scale2) {
236 :
237 : // Check for initialization and ensure minimal scale2 value.
238 0 : if (!isInit) return 0.;
239 0 : if (scale2 < scale2Min) scale2 = scale2Min;
240 :
241 : // If equal to old scale then same answer.
242 0 : if (scale2 == scale2Now && (order < 2 || !lastCallToFull)) return valueNow;
243 0 : scale2Now = scale2;
244 0 : lastCallToFull = false;
245 :
246 : // Fix alpha_S.
247 0 : if (order == 0) {
248 0 : valueNow = valueRef;
249 :
250 : // First/second order alpha_s: differs by mass region.
251 0 : } else {
252 0 : if (scale2 > mt2 && nfmax >= 6)
253 0 : valueNow = 12. * M_PI / (21. * log(scale2/Lambda6Save2));
254 0 : else if (scale2 > mb2)
255 0 : valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
256 0 : else if (scale2 > mc2)
257 0 : valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
258 0 : else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
259 : }
260 :
261 : // Done.
262 0 : return valueNow;
263 0 : }
264 :
265 : //--------------------------------------------------------------------------
266 :
267 : // Calculates the second-order extra factor in alpha_s.
268 : // (To be combined with alphaS1Ord.)
269 :
270 : double AlphaStrong::alphaS2OrdCorr( double scale2) {
271 :
272 : // Check for initialization and ensure minimal scale2 value.
273 0 : if (!isInit) return 1.;
274 0 : if (scale2 < scale2Min) scale2 = scale2Min;
275 :
276 : // Only meaningful for second-order calculations.
277 0 : if (order < 2) return 1.;
278 :
279 : // Second order correction term: differs by mass region.
280 : double Lambda2, b1, b2;
281 0 : if (scale2 > mt2 && nfmax >= 6) {
282 0 : Lambda2 = Lambda6Save2;
283 : b1 = 234. / 441.;
284 : b2 = -36855. / 109512.;
285 0 : } else if (scale2 > mb2) {
286 0 : Lambda2 = Lambda5Save2;
287 : b1 = 348. / 529.;
288 : b2 = 224687. / 242208.;
289 0 : } else if (scale2 > mc2) {
290 0 : Lambda2 = Lambda4Save2;
291 : b1 = 462. / 625.;
292 : b2 = 548575. / 426888.;
293 0 : } else {
294 0 : Lambda2 = Lambda3Save2;
295 : b1 = 64. / 81.;
296 : b2 = 938709. / 663552.;
297 : }
298 0 : double logScale = log(scale2/Lambda2);
299 0 : double loglogScale = log(logScale);
300 0 : return ( 1. - b1 * loglogScale / logScale
301 0 : + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
302 :
303 0 : }
304 :
305 : //--------------------------------------------------------------------------
306 :
307 : // muThres(2): tell what values of flavour thresholds are being used.
308 :
309 : double AlphaStrong::muThres( int idQ) {
310 0 : int idAbs = abs(idQ);
311 : // Return the scale of each flavour threshold included in running.
312 0 : if (idAbs == 4) return mc;
313 0 : else if (idAbs == 5) return mb;
314 0 : else if (idAbs == 6 && nfmax >= 6) return mt;
315 : // Else return -1 (indicates no such threshold is included in running).
316 0 : return -1.;
317 0 : }
318 :
319 : double AlphaStrong::muThres2( int idQ) {
320 0 : int idAbs = abs(idQ);
321 : // Return the scale of each flavour threshold included in running.
322 0 : if (idAbs == 4) return mc2;
323 0 : else if (idAbs == 5) return mb2;
324 0 : else if (idAbs == 6 && nfmax >=6) return mt2;
325 : // Else return -1 (indicates no such threshold is included in running).
326 0 : return -1.;
327 0 : }
328 :
329 : //--------------------------------------------------------------------------
330 :
331 : // facCMW: tells what values of the CMW factors are being used (if any).
332 :
333 : double AlphaStrong::facCMW( int NFIn) {
334 : // Return unity if we are not doing CMW rescaling..
335 0 : if (!isInit || !useCMW) return 1.0;
336 : // Else return the NF-dependent value of the CMW rescaling factor.
337 0 : else if (NFIn <= 3) return FACCMW3;
338 0 : else if (NFIn == 4) return FACCMW4;
339 0 : else if (NFIn == 5) return FACCMW5;
340 0 : else return FACCMW6;
341 0 : }
342 :
343 : //==========================================================================
344 :
345 : // The AlphaEM class.
346 :
347 : //--------------------------------------------------------------------------
348 :
349 : // Definitions of static variables.
350 :
351 : // Z0 mass. Used for normalization scale.
352 : const double AlphaEM::MZ = 91.188;
353 :
354 : // Effective thresholds for electron, muon, light quarks, tau+c, b.
355 : const double AlphaEM::Q2STEP[5] = {0.26e-6, 0.011, 0.25, 3.5, 90.};
356 :
357 : // Running coefficients are sum charge2 / 3 pi in pure QED, here slightly
358 : // enhanced for quarks to approximately account for QCD corrections.
359 : const double AlphaEM::BRUNDEF[5] = {0.1061, 0.2122, 0.460, 0.700, 0.725};
360 :
361 : //--------------------------------------------------------------------------
362 :
363 : // Initialize alpha_EM calculation.
364 :
365 : void AlphaEM::init(int orderIn, Settings* settingsPtr) {
366 :
367 : // Order. Read in alpha_EM value at 0 and m_Z, and mass of Z.
368 0 : order = orderIn;
369 0 : alpEM0 = settingsPtr->parm("StandardModel:alphaEM0");
370 0 : alpEMmZ = settingsPtr->parm("StandardModel:alphaEMmZ");
371 0 : mZ2 = MZ * MZ;
372 :
373 : // AlphaEM values at matching scales and matching b value.
374 0 : if (order <= 0) return;
375 0 : for (int i = 0; i < 5; ++i) bRun[i] = BRUNDEF[i];
376 :
377 : // Step down from mZ to tau/charm threshold.
378 0 : alpEMstep[4] = alpEMmZ / ( 1. + alpEMmZ * bRun[4]
379 0 : * log(mZ2 / Q2STEP[4]) );
380 0 : alpEMstep[3] = alpEMstep[4] / ( 1. - alpEMstep[4] * bRun[3]
381 0 : * log(Q2STEP[3] / Q2STEP[4]) );
382 :
383 : // Step up from me to light-quark threshold.
384 0 : alpEMstep[0] = alpEM0;
385 0 : alpEMstep[1] = alpEMstep[0] / ( 1. - alpEMstep[0] * bRun[0]
386 0 : * log(Q2STEP[1] / Q2STEP[0]) );
387 0 : alpEMstep[2] = alpEMstep[1] / ( 1. - alpEMstep[1] * bRun[1]
388 0 : * log(Q2STEP[2] / Q2STEP[1]) );
389 :
390 : // Fit b in range between light-quark and tau/charm to join smoothly.
391 0 : bRun[2] = (1./alpEMstep[3] - 1./alpEMstep[2])
392 0 : / log(Q2STEP[2] / Q2STEP[3]);
393 :
394 0 : }
395 :
396 : //--------------------------------------------------------------------------
397 :
398 : // Calculate alpha_EM value
399 :
400 : double AlphaEM::alphaEM( double scale2) {
401 :
402 : // Fix alphaEM; for order = -1 fixed at m_Z.
403 0 : if (order == 0) return alpEM0;
404 0 : if (order < 0) return alpEMmZ;
405 :
406 : // Running alphaEM.
407 0 : for (int i = 4; i >= 0; --i) if (scale2 > Q2STEP[i])
408 0 : return alpEMstep[i] / (1. - bRun[i] * alpEMstep[i]
409 0 : * log(scale2 / Q2STEP[i]) );
410 0 : return alpEM0;
411 :
412 0 : }
413 :
414 : //==========================================================================
415 :
416 : // The CoupSM class.
417 :
418 : //--------------------------------------------------------------------------
419 :
420 : // Definitions of static variables: charges and axial couplings.
421 : const double CoupSM::efSave[20] = { 0., -1./3., 2./3., -1./3., 2./3., -1./3.,
422 : 2./3., -1./3., 2./3., 0., 0., -1., 0., -1., 0., -1., 0., -1., 0., 0.};
423 : const double CoupSM::afSave[20] = { 0., -1., 1., -1., 1., -1., 1., -1., 1.,
424 : 0., 0., -1., 1., -1., 1., -1., 1., -1., 1., 0.};
425 :
426 : //--------------------------------------------------------------------------
427 :
428 : // Initialize electroweak mixing angle and couplings, and CKM matrix elements.
429 :
430 : void CoupSM::init(Settings& settings, Rndm* rndmPtrIn) {
431 :
432 : // Store input pointer;
433 0 : rndmPtr = rndmPtrIn;
434 :
435 : // Initialize the local AlphaStrong instance.
436 0 : double alphaSvalue = settings.parm("SigmaProcess:alphaSvalue");
437 0 : int alphaSorder = settings.mode("SigmaProcess:alphaSorder");
438 0 : int alphaSnfmax = settings.mode("StandardModel:alphaSnfmax");
439 0 : alphaSlocal.init( alphaSvalue, alphaSorder, alphaSnfmax, false);
440 :
441 : // Initialize the local AlphaEM instance.
442 0 : int order = settings.mode("SigmaProcess:alphaEMorder");
443 0 : alphaEMlocal.init( order, &settings);
444 :
445 : // Read in electroweak mixing angle and the Fermi constant.
446 0 : s2tW = settings.parm("StandardModel:sin2thetaW");
447 0 : c2tW = 1. - s2tW;
448 0 : s2tWbar = settings.parm("StandardModel:sin2thetaWbar");
449 0 : GFermi = settings.parm("StandardModel:GF");
450 :
451 : // Initialize electroweak couplings.
452 0 : for (int i = 0; i < 20; ++i) {
453 0 : vfSave[i] = afSave[i] - 4. * s2tWbar * efSave[i];
454 0 : lfSave[i] = afSave[i] - 2. * s2tWbar * efSave[i];
455 0 : rfSave[i] = - 2. * s2tWbar * efSave[i];
456 0 : ef2Save[i] = pow2(efSave[i]);
457 0 : vf2Save[i] = pow2(vfSave[i]);
458 0 : af2Save[i] = pow2(afSave[i]);
459 0 : efvfSave[i] = efSave[i] * vfSave[i];
460 0 : vf2af2Save[i] = vf2Save[i] + af2Save[i];
461 : }
462 :
463 : // Read in CKM matrix element values and store them.
464 0 : VCKMsave[1][1] = settings.parm("StandardModel:Vud");
465 0 : VCKMsave[1][2] = settings.parm("StandardModel:Vus");
466 0 : VCKMsave[1][3] = settings.parm("StandardModel:Vub");
467 0 : VCKMsave[2][1] = settings.parm("StandardModel:Vcd");
468 0 : VCKMsave[2][2] = settings.parm("StandardModel:Vcs");
469 0 : VCKMsave[2][3] = settings.parm("StandardModel:Vcb");
470 0 : VCKMsave[3][1] = settings.parm("StandardModel:Vtd");
471 0 : VCKMsave[3][2] = settings.parm("StandardModel:Vts");
472 0 : VCKMsave[3][3] = settings.parm("StandardModel:Vtb");
473 :
474 : // Also allow for the potential existence of a fourth generation.
475 0 : VCKMsave[1][4] = settings.parm("FourthGeneration:VubPrime");
476 0 : VCKMsave[2][4] = settings.parm("FourthGeneration:VcbPrime");
477 0 : VCKMsave[3][4] = settings.parm("FourthGeneration:VtbPrime");
478 0 : VCKMsave[4][1] = settings.parm("FourthGeneration:VtPrimed");
479 0 : VCKMsave[4][2] = settings.parm("FourthGeneration:VtPrimes");
480 0 : VCKMsave[4][3] = settings.parm("FourthGeneration:VtPrimeb");
481 0 : VCKMsave[4][4] = settings.parm("FourthGeneration:VtPrimebPrime");
482 :
483 : // Calculate squares of matrix elements.
484 0 : for(int i = 1; i < 5; ++i) for(int j = 1; j < 5; ++j)
485 0 : V2CKMsave[i][j] = pow2(VCKMsave[i][j]);
486 :
487 : // Sum VCKM^2_out sum for given incoming flavour, excluding top as partner.
488 0 : V2CKMout[1] = V2CKMsave[1][1] + V2CKMsave[2][1];
489 0 : V2CKMout[2] = V2CKMsave[1][1] + V2CKMsave[1][2] + V2CKMsave[1][3];
490 0 : V2CKMout[3] = V2CKMsave[1][2] + V2CKMsave[2][2];
491 0 : V2CKMout[4] = V2CKMsave[2][1] + V2CKMsave[2][2] + V2CKMsave[2][3];
492 0 : V2CKMout[5] = V2CKMsave[1][3] + V2CKMsave[2][3];
493 0 : V2CKMout[6] = V2CKMsave[3][1] + V2CKMsave[3][2] + V2CKMsave[3][3];
494 0 : V2CKMout[7] = V2CKMsave[1][4] + V2CKMsave[2][4];
495 0 : V2CKMout[8] = V2CKMsave[4][1] + V2CKMsave[4][2] + V2CKMsave[4][3];
496 0 : for (int i = 11; i <= 18; ++i) V2CKMout[i] = 1.;
497 :
498 0 : }
499 :
500 : //--------------------------------------------------------------------------
501 :
502 : // Return CKM value for incoming flavours (sign irrelevant).
503 :
504 : double CoupSM::VCKMid(int id1, int id2) {
505 :
506 : // Use absolute sign (want to cover both f -> f' W and f fbar' -> W).
507 0 : int id1Abs = abs(id1);
508 0 : int id2Abs = abs(id2);
509 0 : if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.;
510 :
511 : // Ensure proper order before reading out from VCKMsave or lepton match.
512 0 : if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
513 0 : if (id1Abs <= 8 && id2Abs <= 8) return VCKMsave[id1Abs/2][(id2Abs + 1)/2];
514 0 : if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
515 0 : && id2Abs == id1Abs - 1 ) return 1.;
516 :
517 : // No more valid cases.
518 0 : return 0.;
519 :
520 0 : }
521 :
522 : //--------------------------------------------------------------------------
523 :
524 : // Return squared CKM value for incoming flavours (sign irrelevant).
525 :
526 : double CoupSM::V2CKMid(int id1, int id2) {
527 :
528 : // Use absolute sign (want to cover both f -> f' W and f fbar' -> W).
529 0 : int id1Abs = abs(id1);
530 0 : int id2Abs = abs(id2);
531 0 : if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.;
532 :
533 : // Ensure proper order before reading out from V2CKMsave or lepton match.
534 0 : if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
535 0 : if (id1Abs <= 8 && id2Abs <= 8) return V2CKMsave[id1Abs/2][(id2Abs + 1)/2];
536 0 : if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
537 0 : && id2Abs == id1Abs - 1 ) return 1.;
538 :
539 : // No more valid cases.
540 0 : return 0.;
541 :
542 0 : }
543 :
544 : //--------------------------------------------------------------------------
545 :
546 : // Pick an outgoing flavour for given incoming one, given CKM mixing.
547 :
548 : int CoupSM::V2CKMpick(int id) {
549 :
550 : // Initial values.
551 0 : int idIn = abs(id);
552 : int idOut = 0;
553 :
554 : // Quarks: need to make random choice.
555 0 : if (idIn >= 1 && idIn <= 8) {
556 0 : double V2CKMrndm = rndmPtr->flat() * V2CKMout[idIn];
557 0 : if (idIn == 1) idOut = (V2CKMrndm < V2CKMsave[1][1]) ? 2 : 4;
558 0 : else if (idIn == 2) idOut = (V2CKMrndm < V2CKMsave[1][1]) ? 1
559 0 : : ( (V2CKMrndm < V2CKMsave[1][1] + V2CKMsave[1][2]) ? 3 : 5 );
560 0 : else if (idIn == 3) idOut = (V2CKMrndm < V2CKMsave[1][2]) ? 2 : 4;
561 0 : else if (idIn == 4) idOut = (V2CKMrndm < V2CKMsave[2][1]) ? 1
562 0 : : ( (V2CKMrndm < V2CKMsave[2][1] + V2CKMsave[2][2]) ? 3 : 5 );
563 0 : else if (idIn == 5) idOut = (V2CKMrndm < V2CKMsave[1][3]) ? 2 : 4;
564 0 : else if (idIn == 6) idOut = (V2CKMrndm < V2CKMsave[3][1]) ? 1
565 0 : : ( (V2CKMrndm < V2CKMsave[3][1] + V2CKMsave[3][2]) ? 3 : 5 );
566 0 : else if (idIn == 7) idOut = (V2CKMrndm < V2CKMsave[1][4]) ? 2 : 4;
567 0 : else if (idIn == 8) idOut = (V2CKMrndm < V2CKMsave[4][1]) ? 1
568 0 : : ( (V2CKMrndm < V2CKMsave[4][1] + V2CKMsave[4][2]) ? 3 : 5 );
569 :
570 : // Leptons: unambiguous.
571 0 : } else if (idIn >= 11 && idIn <= 18) {
572 0 : if (idIn%2 == 1) idOut = idIn + 1;
573 0 : else idOut = idIn - 1;
574 : }
575 :
576 : // Done. Return with sign.
577 0 : return ( (id > 0) ? idOut : -idOut );
578 :
579 : }
580 :
581 : //==========================================================================
582 :
583 : } // end namespace Pythia8
|
