StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StandardModel.cc
1 // StandardModel.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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  if (mt <= 1.) setThresholds(1.5, 4.8, 171.0);
47 
48  // Order of alpha_s evaluation.Default values.
49  valueRef = valueIn;
50  order = max( 0, min( 2, orderIn ) );
51  nfmax = max(5,min(6,nfmaxIn));
52  useCMW = useCMWIn;
53  lastCallToFull = false;
54  Lambda3Save = Lambda4Save = Lambda5Save = Lambda6Save = scale2Min = 0.;
55 
56  // Fix alpha_s.
57  if (order == 0) {
58 
59  // First order alpha_s: match at flavour thresholds.
60  } else if (order == 1) {
61  Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueRef) );
62  Lambda6Save = Lambda5Save * pow(Lambda5Save/mt, 2./21.);
63  Lambda4Save = Lambda5Save * pow(mb/Lambda5Save, 2./25.);
64  Lambda3Save = Lambda4Save * pow(mc/Lambda4Save, 2./27.);
65 
66  // Second order alpha_s: iterative match at flavour thresholds.
67  } 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  Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueRef) );
82  for (int iter = 0; iter < NITER; ++iter) {
83  logScale = 2. * log(MZ/Lambda5Save);
84  loglogScale = log(logScale);
85  correction = 1. - b15 * loglogScale / logScale
86  + pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25);
87  valueIter = valueRef / correction;
88  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  double logScaleT = 2. * log(mt/Lambda5Save);
93  double loglogScaleT = log(logScaleT);
94  double valueT = 12. * M_PI / (23. * logScaleT)
95  * (1. - b15 * loglogScaleT / logScaleT
96  + pow2(b15 / logScaleT) * (pow2(loglogScaleT - 0.5) + b25 - 1.25) );
97  Lambda6Save = Lambda5Save;
98  for (int iter = 0; iter < NITER; ++iter) {
99  logScale = 2. * log(mt/Lambda6Save);
100  loglogScale = log(logScale);
101  correction = 1. - b16 * loglogScale / logScale
102  + pow2(b16 / logScale) * (pow2(loglogScale - 0.5) + b26 - 1.25);
103  valueIter = valueT / correction;
104  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  double logScaleB = 2. * log(mb/Lambda5Save);
109  double loglogScaleB = log(logScaleB);
110  double valueB = 12. * M_PI / (23. * logScaleB)
111  * (1. - b15 * loglogScaleB / logScaleB
112  + pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) );
113  Lambda4Save = Lambda5Save;
114  for (int iter = 0; iter < NITER; ++iter) {
115  logScale = 2. * log(mb/Lambda4Save);
116  loglogScale = log(logScale);
117  correction = 1. - b14 * loglogScale / logScale
118  + pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25);
119  valueIter = valueB / correction;
120  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  double logScaleC = 2. * log(mc/Lambda4Save);
124  double loglogScaleC = log(logScaleC);
125  double valueC = 12. * M_PI / (25. * logScaleC)
126  * (1. - b14 * loglogScaleC / logScaleC
127  + pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) );
128  Lambda3Save = Lambda4Save;
129  for (int iter = 0; iter < NITER; ++iter) {
130  logScale = 2. * log(mc/Lambda3Save);
131  loglogScale = log(logScale);
132  correction = 1. - b13 * loglogScale / logScale
133  + pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25);
134  valueIter = valueC / correction;
135  Lambda3Save = mc * exp( -6. * M_PI / (27. * valueIter) );
136  }
137  }
138 
139  // Optionally rescale Lambda values by CMW factor.
140  if (useCMW) {
141  Lambda3Save *= FACCMW3;
142  Lambda4Save *= FACCMW4;
143  Lambda5Save *= FACCMW5;
144  Lambda6Save *= FACCMW6;
145  }
146 
147  // Impose SAFETYMARGINs to prevent getting too close to LambdaQCD.
148  if (order == 1) scale2Min = pow2(SAFETYMARGIN1 * Lambda3Save);
149  else if (order == 2) scale2Min = pow2(SAFETYMARGIN2 * Lambda3Save);
150 
151  // Save squares of mass and Lambda values as well.
152  Lambda3Save2 = pow2(Lambda3Save);
153  Lambda4Save2 = pow2(Lambda4Save);
154  Lambda5Save2 = pow2(Lambda5Save);
155  Lambda6Save2 = pow2(Lambda6Save);
156  mc2 = pow2(mc);
157  mb2 = pow2(mb);
158  mt2 = pow2(mt);
159  valueNow = valueIn;
160  scale2Now = MZ * MZ;
161  isInit = true;
162 
163 }
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  if (!isInit) return 0.;
173  if (scale2 < scale2Min) scale2 = scale2Min;
174 
175  // If equal to old scale then same answer.
176  if (scale2 == scale2Now && (order < 2 || lastCallToFull)) return valueNow;
177  scale2Now = scale2;
178  lastCallToFull = true;
179 
180  // Fix alpha_s.
181  if (order == 0) {
182  valueNow = valueRef;
183 
184  // First order alpha_s: differs by mass region.
185  } else if (order == 1) {
186  if (scale2 > mt2 && nfmax >= 6)
187  valueNow = 12. * M_PI / (21. * log(scale2/Lambda6Save2));
188  else if (scale2 > mb2)
189  valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
190  else if (scale2 > mc2)
191  valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
192  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  if (scale2 > mt2 && nfmax >= 6) {
198  Lambda2 = Lambda6Save2;
199  b0 = 21.;
200  b1 = 234. / 441.;
201  b2 = -36855. / 109512.;
202  } else if (scale2 > mb2) {
203  Lambda2 = Lambda5Save2;
204  b0 = 23.;
205  b1 = 348. / 529.;
206  b2 = 224687. / 242208.;
207  } else if (scale2 > mc2) {
208  Lambda2 = Lambda4Save2;
209  b0 = 25.;
210  b1 = 462. / 625.;
211  b2 = 548575. / 426888.;
212  } else {
213  Lambda2 = Lambda3Save2;
214  b0 = 27.;
215  b1 = 64. / 81.;
216  b2 = 938709. / 663552.;
217  }
218  double logScale = log(scale2/Lambda2);
219  double loglogScale = log(logScale);
220  valueNow = 12. * M_PI / (b0 * logScale)
221  * ( 1. - b1 * loglogScale / logScale
222  + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
223  }
224 
225  // Done.
226  return valueNow;
227 
228 }
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  if (!isInit) return 0.;
239  if (scale2 < scale2Min) scale2 = scale2Min;
240 
241  // If equal to old scale then same answer.
242  if (scale2 == scale2Now && (order < 2 || !lastCallToFull)) return valueNow;
243  scale2Now = scale2;
244  lastCallToFull = false;
245 
246  // Fix alpha_S.
247  if (order == 0) {
248  valueNow = valueRef;
249 
250  // First/second order alpha_s: differs by mass region.
251  } else {
252  if (scale2 > mt2 && nfmax >= 6)
253  valueNow = 12. * M_PI / (21. * log(scale2/Lambda6Save2));
254  else if (scale2 > mb2)
255  valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
256  else if (scale2 > mc2)
257  valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
258  else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
259  }
260 
261  // Done.
262  return valueNow;
263 }
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  if (!isInit) return 1.;
274  if (scale2 < scale2Min) scale2 = scale2Min;
275 
276  // Only meaningful for second-order calculations.
277  if (order < 2) return 1.;
278 
279  // Second order correction term: differs by mass region.
280  double Lambda2, b1, b2;
281  if (scale2 > mt2 && nfmax >= 6) {
282  Lambda2 = Lambda6Save2;
283  b1 = 234. / 441.;
284  b2 = -36855. / 109512.;
285  } else if (scale2 > mb2) {
286  Lambda2 = Lambda5Save2;
287  b1 = 348. / 529.;
288  b2 = 224687. / 242208.;
289  } else if (scale2 > mc2) {
290  Lambda2 = Lambda4Save2;
291  b1 = 462. / 625.;
292  b2 = 548575. / 426888.;
293  } else {
294  Lambda2 = Lambda3Save2;
295  b1 = 64. / 81.;
296  b2 = 938709. / 663552.;
297  }
298  double logScale = log(scale2/Lambda2);
299  double loglogScale = log(logScale);
300  return ( 1. - b1 * loglogScale / logScale
301  + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
302 
303 }
304 
305 //--------------------------------------------------------------------------
306 
307 // muThres(2): tell what values of flavour thresholds are being used.
308 
309 double AlphaStrong::muThres( int idQ) {
310  int idAbs = abs(idQ);
311  // Return the scale of each flavour threshold included in running.
312  if (idAbs == 4) return mc;
313  else if (idAbs == 5) return mb;
314  else if (idAbs == 6 && nfmax >= 6) return mt;
315  // Else return -1 (indicates no such threshold is included in running).
316  return -1.;
317 }
318 
319 double AlphaStrong::muThres2( int idQ) {
320  int idAbs = abs(idQ);
321  // Return the scale of each flavour threshold included in running.
322  if (idAbs == 4) return mc2;
323  else if (idAbs == 5) return mb2;
324  else if (idAbs == 6 && nfmax >=6) return mt2;
325  // Else return -1 (indicates no such threshold is included in running).
326  return -1.;
327 }
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  if (!isInit || !useCMW) return 1.0;
336  // Else return the NF-dependent value of the CMW rescaling factor.
337  else if (NFIn <= 3) return FACCMW3;
338  else if (NFIn == 4) return FACCMW4;
339  else if (NFIn == 5) return FACCMW5;
340  else return FACCMW6;
341 }
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  order = orderIn;
369  alpEM0 = settingsPtr->parm("StandardModel:alphaEM0");
370  alpEMmZ = settingsPtr->parm("StandardModel:alphaEMmZ");
371  mZ2 = MZ * MZ;
372 
373  // AlphaEM values at matching scales and matching b value.
374  if (order <= 0) return;
375  for (int i = 0; i < 5; ++i) bRun[i] = BRUNDEF[i];
376 
377  // Step down from mZ to tau/charm threshold.
378  alpEMstep[4] = alpEMmZ / ( 1. + alpEMmZ * bRun[4]
379  * log(mZ2 / Q2STEP[4]) );
380  alpEMstep[3] = alpEMstep[4] / ( 1. - alpEMstep[4] * bRun[3]
381  * log(Q2STEP[3] / Q2STEP[4]) );
382 
383  // Step up from me to light-quark threshold.
384  alpEMstep[0] = alpEM0;
385  alpEMstep[1] = alpEMstep[0] / ( 1. - alpEMstep[0] * bRun[0]
386  * log(Q2STEP[1] / Q2STEP[0]) );
387  alpEMstep[2] = alpEMstep[1] / ( 1. - alpEMstep[1] * bRun[1]
388  * log(Q2STEP[2] / Q2STEP[1]) );
389 
390  // Fit b in range between light-quark and tau/charm to join smoothly.
391  bRun[2] = (1./alpEMstep[3] - 1./alpEMstep[2])
392  / log(Q2STEP[2] / Q2STEP[3]);
393 
394 }
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  if (order == 0) return alpEM0;
404  if (order < 0) return alpEMmZ;
405 
406  // Running alphaEM.
407  for (int i = 4; i >= 0; --i) if (scale2 > Q2STEP[i])
408  return alpEMstep[i] / (1. - bRun[i] * alpEMstep[i]
409  * log(scale2 / Q2STEP[i]) );
410  return alpEM0;
411 
412 }
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  rndmPtr = rndmPtrIn;
434 
435  // Initialize the local AlphaStrong instance.
436  double alphaSvalue = settings.parm("SigmaProcess:alphaSvalue");
437  int alphaSorder = settings.mode("SigmaProcess:alphaSorder");
438  int alphaSnfmax = settings.mode("StandardModel:alphaSnfmax");
439  alphaSlocal.init( alphaSvalue, alphaSorder, alphaSnfmax, false);
440 
441  // Initialize the local AlphaEM instance.
442  int order = settings.mode("SigmaProcess:alphaEMorder");
443  alphaEMlocal.init( order, &settings);
444 
445  // Read in electroweak mixing angle and the Fermi constant.
446  s2tW = settings.parm("StandardModel:sin2thetaW");
447  c2tW = 1. - s2tW;
448  s2tWbar = settings.parm("StandardModel:sin2thetaWbar");
449  GFermi = settings.parm("StandardModel:GF");
450 
451  // Initialize electroweak couplings.
452  for (int i = 0; i < 20; ++i) {
453  vfSave[i] = afSave[i] - 4. * s2tWbar * efSave[i];
454  lfSave[i] = afSave[i] - 2. * s2tWbar * efSave[i];
455  rfSave[i] = - 2. * s2tWbar * efSave[i];
456  ef2Save[i] = pow2(efSave[i]);
457  vf2Save[i] = pow2(vfSave[i]);
458  af2Save[i] = pow2(afSave[i]);
459  efvfSave[i] = efSave[i] * vfSave[i];
460  vf2af2Save[i] = vf2Save[i] + af2Save[i];
461  }
462 
463  // Read in CKM matrix element values and store them.
464  VCKMsave[1][1] = settings.parm("StandardModel:Vud");
465  VCKMsave[1][2] = settings.parm("StandardModel:Vus");
466  VCKMsave[1][3] = settings.parm("StandardModel:Vub");
467  VCKMsave[2][1] = settings.parm("StandardModel:Vcd");
468  VCKMsave[2][2] = settings.parm("StandardModel:Vcs");
469  VCKMsave[2][3] = settings.parm("StandardModel:Vcb");
470  VCKMsave[3][1] = settings.parm("StandardModel:Vtd");
471  VCKMsave[3][2] = settings.parm("StandardModel:Vts");
472  VCKMsave[3][3] = settings.parm("StandardModel:Vtb");
473 
474  // Also allow for the potential existence of a fourth generation.
475  VCKMsave[1][4] = settings.parm("FourthGeneration:VubPrime");
476  VCKMsave[2][4] = settings.parm("FourthGeneration:VcbPrime");
477  VCKMsave[3][4] = settings.parm("FourthGeneration:VtbPrime");
478  VCKMsave[4][1] = settings.parm("FourthGeneration:VtPrimed");
479  VCKMsave[4][2] = settings.parm("FourthGeneration:VtPrimes");
480  VCKMsave[4][3] = settings.parm("FourthGeneration:VtPrimeb");
481  VCKMsave[4][4] = settings.parm("FourthGeneration:VtPrimebPrime");
482 
483  // Calculate squares of matrix elements.
484  for(int i = 1; i < 5; ++i) for(int j = 1; j < 5; ++j)
485  V2CKMsave[i][j] = pow2(VCKMsave[i][j]);
486 
487  // Sum VCKM^2_out sum for given incoming flavour, excluding top as partner.
488  V2CKMout[1] = V2CKMsave[1][1] + V2CKMsave[2][1];
489  V2CKMout[2] = V2CKMsave[1][1] + V2CKMsave[1][2] + V2CKMsave[1][3];
490  V2CKMout[3] = V2CKMsave[1][2] + V2CKMsave[2][2];
491  V2CKMout[4] = V2CKMsave[2][1] + V2CKMsave[2][2] + V2CKMsave[2][3];
492  V2CKMout[5] = V2CKMsave[1][3] + V2CKMsave[2][3];
493  V2CKMout[6] = V2CKMsave[3][1] + V2CKMsave[3][2] + V2CKMsave[3][3];
494  V2CKMout[7] = V2CKMsave[1][4] + V2CKMsave[2][4];
495  V2CKMout[8] = V2CKMsave[4][1] + V2CKMsave[4][2] + V2CKMsave[4][3];
496  for (int i = 11; i <= 18; ++i) V2CKMout[i] = 1.;
497 
498 }
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  int id1Abs = abs(id1);
508  int id2Abs = abs(id2);
509  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  if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
513  if (id1Abs <= 8 && id2Abs <= 8) return VCKMsave[id1Abs/2][(id2Abs + 1)/2];
514  if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
515  && id2Abs == id1Abs - 1 ) return 1.;
516 
517  // No more valid cases.
518  return 0.;
519 
520 }
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  int id1Abs = abs(id1);
530  int id2Abs = abs(id2);
531  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  if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
535  if (id1Abs <= 8 && id2Abs <= 8) return V2CKMsave[id1Abs/2][(id2Abs + 1)/2];
536  if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
537  && id2Abs == id1Abs - 1 ) return 1.;
538 
539  // No more valid cases.
540  return 0.;
541 
542 }
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  int idIn = abs(id);
552  int idOut = 0;
553 
554  // Quarks: need to make random choice.
555  if (idIn >= 1 && idIn <= 8) {
556  double V2CKMrndm = rndmPtr->flat() * V2CKMout[idIn];
557  if (idIn == 1) idOut = (V2CKMrndm < V2CKMsave[1][1]) ? 2 : 4;
558  else if (idIn == 2) idOut = (V2CKMrndm < V2CKMsave[1][1]) ? 1
559  : ( (V2CKMrndm < V2CKMsave[1][1] + V2CKMsave[1][2]) ? 3 : 5 );
560  else if (idIn == 3) idOut = (V2CKMrndm < V2CKMsave[1][2]) ? 2 : 4;
561  else if (idIn == 4) idOut = (V2CKMrndm < V2CKMsave[2][1]) ? 1
562  : ( (V2CKMrndm < V2CKMsave[2][1] + V2CKMsave[2][2]) ? 3 : 5 );
563  else if (idIn == 5) idOut = (V2CKMrndm < V2CKMsave[1][3]) ? 2 : 4;
564  else if (idIn == 6) idOut = (V2CKMrndm < V2CKMsave[3][1]) ? 1
565  : ( (V2CKMrndm < V2CKMsave[3][1] + V2CKMsave[3][2]) ? 3 : 5 );
566  else if (idIn == 7) idOut = (V2CKMrndm < V2CKMsave[1][4]) ? 2 : 4;
567  else if (idIn == 8) idOut = (V2CKMrndm < V2CKMsave[4][1]) ? 1
568  : ( (V2CKMrndm < V2CKMsave[4][1] + V2CKMsave[4][2]) ? 3 : 5 );
569 
570  // Leptons: unambiguous.
571  } else if (idIn >= 11 && idIn <= 18) {
572  if (idIn%2 == 1) idOut = idIn + 1;
573  else idOut = idIn - 1;
574  }
575 
576  // Done. Return with sign.
577  return ( (id > 0) ? idOut : -idOut );
578 
579 }
580 
581 //==========================================================================
582 
583 } // end namespace Pythia8