]>
Commit | Line | Data |
---|---|---|
5ad4eb21 | 1 | // StandardModel.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2008 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 "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 | // Always evaluate running alpha_s above Lambda3 to avoid disaster. | |
25 | // Safety margin picked to freeze roughly for alpha_s = 10. | |
26 | const double AlphaStrong::SAFETYMARGIN1 = 1.07; | |
27 | const double AlphaStrong::SAFETYMARGIN2 = 1.33; | |
28 | ||
29 | //********* | |
30 | ||
31 | // Initialize alpha_strong calculation by finding Lambda values etc. | |
32 | ||
33 | void AlphaStrong::init( double valueIn, int orderIn) { | |
34 | ||
35 | // Order of alpha_s evaluation. Charm, bottom and Z0 masses. | |
36 | // Pick defaults if ParticleDataTable not properly initialized. | |
37 | valueRef = valueIn; | |
38 | order = max( 0, min( 2, orderIn ) ); | |
39 | mc = ParticleDataTable::m0(4); | |
40 | if (mc < 1. || mc > 2.) mc = 1.5; | |
41 | mb = ParticleDataTable::m0(5); | |
42 | if (mb < 3. || mb > 6.) mb = 4.8; | |
43 | mZ = ParticleDataTable::m0(23); | |
44 | if (mZ < 90. || mZ > 92.) mZ = 91.188; | |
45 | ||
46 | // Fix alpha_s. | |
47 | if (order == 0) { | |
48 | Lambda3Save = Lambda4Save = Lambda5Save = scale2Min = 0.; | |
49 | ||
50 | // First order alpha_s: match at flavour thresholds. | |
51 | } else if (order == 1) { | |
52 | Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueRef) ); | |
53 | Lambda4Save = Lambda5Save * pow(mb/Lambda5Save, 2./25.); | |
54 | Lambda3Save = Lambda4Save * pow(mc/Lambda4Save, 2./27.); | |
55 | scale2Min = pow2(SAFETYMARGIN1 * Lambda3Save); | |
56 | ||
57 | // Second order alpha_s: iterative match at flavour thresholds. | |
58 | } else { | |
59 | double b15 = 348. / 529.; | |
60 | double b14 = 462. / 625.; | |
61 | double b13 = 64. / 81.; | |
62 | double b25 = 224687. / 242208.; | |
63 | double b24 = 548575. / 426888.; | |
64 | double b23 = 938709. / 663552.; | |
65 | double logScale, loglogScale, correction, valueIter; | |
66 | ||
67 | // Find Lambda_5 at m_Z. | |
68 | Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueRef) ); | |
69 | for (int iter = 0; iter < NITER; ++iter) { | |
70 | logScale = 2. * log(mZ/Lambda5Save); | |
71 | loglogScale = log(logScale); | |
72 | correction = 1. - b15 * loglogScale / logScale | |
73 | + pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25); | |
74 | valueIter = valueRef / correction; | |
75 | Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueIter) ); | |
76 | } | |
77 | ||
78 | // Find Lambda_4 at m_b. | |
79 | double logScaleB = 2. * log(mb/Lambda5Save); | |
80 | double loglogScaleB = log(logScaleB); | |
81 | double valueB = 12. * M_PI / (23. * logScaleB) | |
82 | * (1. - b15 * loglogScaleB / logScaleB | |
83 | + pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) ); | |
84 | Lambda4Save = Lambda5Save; | |
85 | for (int iter = 0; iter < NITER; ++iter) { | |
86 | logScale = 2. * log(mb/Lambda4Save); | |
87 | loglogScale = log(logScale); | |
88 | correction = 1. - b14 * loglogScale / logScale | |
89 | + pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25); | |
90 | valueIter = valueB / correction; | |
91 | Lambda4Save = mb * exp( -6. * M_PI / (25. * valueIter) ); | |
92 | } | |
93 | ||
94 | // Find Lambda_3 at m_c. | |
95 | double logScaleC = 2. * log(mc/Lambda4Save); | |
96 | double loglogScaleC = log(logScaleC); | |
97 | double valueC = 12. * M_PI / (25. * logScaleC) | |
98 | * (1. - b14 * loglogScaleC / logScaleC | |
99 | + pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) ); | |
100 | Lambda3Save = Lambda4Save; | |
101 | for (int iter = 0; iter < NITER; ++iter) { | |
102 | logScale = 2. * log(mc/Lambda3Save); | |
103 | loglogScale = log(logScale); | |
104 | correction = 1. - b13 * loglogScale / logScale | |
105 | + pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25); | |
106 | valueIter = valueC / correction; | |
107 | Lambda3Save = mc * exp( -6. * M_PI / (27. * valueIter) ); | |
108 | } | |
109 | scale2Min = pow2(SAFETYMARGIN2 * Lambda3Save); | |
110 | } | |
111 | ||
112 | // Save squares of mass and Lambda values as well. | |
113 | mc2 = pow2(mc); | |
114 | mb2 = pow2(mb); | |
115 | Lambda3Save2 = pow2(Lambda3Save); | |
116 | Lambda4Save2 = pow2(Lambda4Save); | |
117 | Lambda5Save2 = pow2(Lambda5Save); | |
118 | valueNow = valueIn; | |
119 | scale2Now = mZ * mZ; | |
120 | isInit = true; | |
121 | ||
122 | } | |
123 | ||
124 | //********* | |
125 | ||
126 | // Calculate alpha_s value | |
127 | ||
128 | double AlphaStrong::alphaS( double scale2) { | |
129 | ||
130 | // Check for initialization and ensure minimal scale2 value. | |
131 | if (!isInit) return 0.; | |
132 | if (scale2 < scale2Min) scale2 = scale2Min; | |
133 | ||
134 | // If equal to old scale then same answer. | |
135 | if (scale2 == scale2Now && (order < 2 || lastCallToFull)) return valueNow; | |
136 | scale2Now = scale2; | |
137 | lastCallToFull = true; | |
138 | ||
139 | // Fix alpha_s. | |
140 | if (order == 0) { | |
141 | valueNow = valueRef; | |
142 | ||
143 | // First order alpha_s: differs by mass region. | |
144 | } else if (order == 1) { | |
145 | if (scale2 > mb2) | |
146 | valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2)); | |
147 | else if (scale2 > mc2) | |
148 | valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2)); | |
149 | else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2)); | |
150 | ||
151 | // Second order alpha_s: differs by mass region. | |
152 | } else { | |
153 | double Lambda2, b0, b1, b2; | |
154 | if (scale2 > mb2) { | |
155 | Lambda2 = Lambda5Save2; | |
156 | b0 = 23.; | |
157 | b1 = 348. / 529.; | |
158 | b2 = 224687. / 242208.; | |
159 | } else if (scale2 > mc2) { | |
160 | Lambda2 = Lambda4Save2; | |
161 | b0 = 25.; | |
162 | b1 = 462. / 625.; | |
163 | b2 = 548575. / 426888.; | |
164 | } else { | |
165 | Lambda2 = Lambda3Save2; | |
166 | b0 = 27.; | |
167 | b1 = 64. / 81.; | |
168 | b2 = 938709. / 663552.; | |
169 | } | |
170 | double logScale = log(scale2/Lambda2); | |
171 | double loglogScale = log(logScale); | |
172 | valueNow = 12. * M_PI / (b0 * logScale) | |
173 | * ( 1. - b1 * loglogScale / logScale | |
174 | + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) ); | |
175 | } | |
176 | ||
177 | // Done. | |
178 | return valueNow; | |
179 | ||
180 | } | |
181 | ||
182 | //********* | |
183 | ||
184 | // Calculate alpha_s value, but only use up to first-order piece. | |
185 | // (To be combined with alphaS2OrdCorr.) | |
186 | ||
187 | double AlphaStrong::alphaS1Ord( double scale2) { | |
188 | ||
189 | // Check for initialization and ensure minimal scale2 value. | |
190 | if (!isInit) return 0.; | |
191 | if (scale2 < scale2Min) scale2 = scale2Min; | |
192 | ||
193 | // If equal to old scale then same answer. | |
194 | if (scale2 == scale2Now && (order < 2 || !lastCallToFull)) return valueNow; | |
195 | scale2Now = scale2; | |
196 | lastCallToFull = false; | |
197 | ||
198 | // Fix alpha_S. | |
199 | if (order == 0) { | |
200 | valueNow = valueRef; | |
201 | ||
202 | // First/second order alpha_s: differs by mass region. | |
203 | } else { | |
204 | if (scale2 > mb2) | |
205 | valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2)); | |
206 | else if (scale2 > mc2) | |
207 | valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2)); | |
208 | else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2)); | |
209 | } | |
210 | ||
211 | // Done. | |
212 | return valueNow; | |
213 | } | |
214 | ||
215 | //********* | |
216 | ||
217 | // Calculates the second-order extra factor in alpha_s. | |
218 | // (To be combined with alphaS1Ord.) | |
219 | ||
220 | double AlphaStrong::alphaS2OrdCorr( double scale2) { | |
221 | ||
222 | // Check for initialization and ensure minimal scale2 value. | |
223 | if (!isInit) return 1.; | |
224 | if (scale2 < scale2Min) scale2 = scale2Min; | |
225 | ||
226 | // Only meaningful for second order calculations. | |
227 | if (order < 2) return 1.; | |
228 | ||
229 | // Second order correction term: differs by mass region. | |
230 | double Lambda2, b0, b1, b2; | |
231 | if (scale2 > mb2) { | |
232 | Lambda2 = Lambda5Save2; | |
233 | b0 = 23.; | |
234 | b1 = 348. / 529.; | |
235 | b2 = 224687. / 242208.; | |
236 | } else if (scale2 > mc2) { | |
237 | Lambda2 = Lambda4Save2; | |
238 | b0 = 25.; | |
239 | b1 = 462. / 625.; | |
240 | b2 = 548575. / 426888.; | |
241 | } else { | |
242 | Lambda2 = Lambda3Save2; | |
243 | b0 = 27.; | |
244 | b1 = 64. / 81.; | |
245 | b2 = 938709. / 663552.; | |
246 | } | |
247 | double logScale = log(scale2/Lambda2); | |
248 | double loglogScale = log(logScale); | |
249 | return ( 1. - b1 * loglogScale / logScale | |
250 | + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) ); | |
251 | ||
252 | } | |
253 | ||
254 | //************************************************************************** | |
255 | ||
256 | // The AlphaEM class. | |
257 | ||
258 | //********* | |
259 | ||
260 | // Definitions of static variables. | |
261 | double AlphaEM::alpEM0 = 0.00729735; | |
262 | double AlphaEM::alpEMmZ = 0.00781751; | |
263 | double AlphaEM::mZ2 = 8315.; | |
264 | // Effective thresholds for electron, muon, light quarks, tau+c, b. | |
265 | double AlphaEM::Q2step[5] = {0.26e-6, 0.011, 0.25, 3.5, 90.}; | |
266 | // Running coefficients are sum charge2 / 3 pi in pure QED, here slightly | |
267 | // enhanced for quarks to approximately account for QCD corrections. | |
268 | double AlphaEM::bRun[5] = {0.1061, 0.2122, 0.460, 0.700, 0.725}; | |
269 | double AlphaEM::alpEMstep[5] = {}; | |
270 | ||
271 | //********* | |
272 | ||
273 | // Initialize alpha_EM calculation. | |
274 | ||
275 | void AlphaEM::initStatic() { | |
276 | ||
277 | // Read in alpha_EM value at 0 and m_Z, and mass of Z. | |
278 | alpEM0 = Settings::parm("StandardModel:alphaEM0"); | |
279 | alpEMmZ = Settings::parm("StandardModel:alphaEMmZ"); | |
280 | double mZ = ParticleDataTable::m0(23); | |
281 | mZ2 = mZ * mZ; | |
282 | ||
283 | // AlphaEM values at matching scales and matching b value. | |
284 | ||
285 | // Step down from mZ to tau/charm threshold. | |
286 | alpEMstep[4] = alpEMmZ / ( 1. + alpEMmZ * bRun[4] | |
287 | * log(mZ2 / Q2step[4]) ); | |
288 | alpEMstep[3] = alpEMstep[4] / ( 1. - alpEMstep[4] * bRun[3] | |
289 | * log(Q2step[3] / Q2step[4]) ); | |
290 | ||
291 | // Step up from me to light-quark threshold. | |
292 | alpEMstep[0] = alpEM0; | |
293 | alpEMstep[1] = alpEMstep[0] / ( 1. - alpEMstep[0] * bRun[0] | |
294 | * log(Q2step[1] / Q2step[0]) ); | |
295 | alpEMstep[2] = alpEMstep[1] / ( 1. - alpEMstep[1] * bRun[1] | |
296 | * log(Q2step[2] / Q2step[1]) ); | |
297 | ||
298 | // Fit b in range between light-quark and tau/charm to join smoothly. | |
299 | bRun[2] = (1./alpEMstep[3] - 1./alpEMstep[2]) | |
300 | / log(Q2step[2] / Q2step[3]); | |
301 | ||
302 | } | |
303 | ||
304 | //********* | |
305 | ||
306 | // Calculate alpha_EM value | |
307 | ||
308 | double AlphaEM::alphaEM( double scale2) { | |
309 | ||
310 | // Fix alphaEM; for order = -1 fixed at m_Z. | |
311 | if (order == 0) return alpEM0; | |
312 | if (order < 0) return alpEMmZ; | |
313 | ||
314 | // Running alphaEM. | |
315 | for (int i = 4; i >= 0; --i) if (scale2 > Q2step[i]) | |
316 | return alpEMstep[i] / (1. - bRun[i] * alpEMstep[i] | |
317 | * log(scale2 / Q2step[i]) ); | |
318 | return alpEM0; | |
319 | ||
320 | } | |
321 | ||
322 | //************************************************************************** | |
323 | ||
324 | // The CoupEW class. | |
325 | ||
326 | //********* | |
327 | ||
328 | // Definitions of static variables. | |
329 | ||
330 | double CoupEW::s2tW = 0.232; | |
331 | double CoupEW::c2tW = 0.768; | |
332 | double CoupEW::s2tWbar = 0.232; | |
333 | double CoupEW::efSave[20] = { 0., -1./3., 2./3., -1./3., 2./3., -1./3., | |
334 | 2./3., -1./3., 2./3., 0., 0., -1., 0., -1., 0., -1., 0., -1., 0., 0.}; | |
335 | double CoupEW::vfSave[20] = { }; | |
336 | double CoupEW::afSave[20] = { 0., -1., 1., -1., 1., -1., 1., -1., 1., | |
337 | 0., 0., -1., 1., -1., 1., -1., 1., -1., 1., 0.}; | |
338 | double CoupEW::lfSave[20] = { }; | |
339 | double CoupEW::rfSave[20] = { }; | |
340 | double CoupEW::ef2Save[20] = { }; | |
341 | double CoupEW::vf2Save[20] = { }; | |
342 | double CoupEW::af2Save[20] = { }; | |
343 | double CoupEW::efvfSave[20] = { }; | |
344 | double CoupEW::vf2af2Save[20] = { }; | |
345 | ||
346 | //********* | |
347 | ||
348 | // Initialize electroweak mixing angle and couplings. | |
349 | ||
350 | void CoupEW::initStatic() { | |
351 | ||
352 | // Read in electroweak mixing angle. | |
353 | s2tW = Settings::parm("StandardModel:sin2thetaW"); | |
354 | c2tW = 1. - s2tW; | |
355 | s2tWbar = Settings::parm("StandardModel:sin2thetaWbar"); | |
356 | ||
357 | // Initialize electroweak couplings. | |
358 | for (int i = 0; i < 20; ++i) { | |
359 | vfSave[i] = afSave[i] - 4. * s2tWbar * efSave[i]; | |
360 | lfSave[i] = afSave[i] - 2. * s2tWbar * efSave[i]; | |
361 | rfSave[i] = - 2. * s2tWbar * efSave[i]; | |
362 | ef2Save[i] = pow2(efSave[i]); | |
363 | vf2Save[i] = pow2(vfSave[i]); | |
364 | af2Save[i] = pow2(afSave[i]); | |
365 | efvfSave[i] = efSave[i] * vfSave[i]; | |
366 | vf2af2Save[i] = vf2Save[i] + af2Save[i]; | |
367 | } | |
368 | ||
369 | } | |
370 | ||
371 | //************************************************************************** | |
372 | ||
373 | // The VCKM class. | |
374 | ||
375 | //********* | |
376 | ||
377 | // Definitions of static variables. Initialize to all elements zero. | |
378 | ||
379 | double VCKM::Vsave[5][5] = { }; | |
380 | double VCKM::V2save[5][5] = { }; | |
381 | double VCKM::V2out[20] = { }; | |
382 | ||
383 | //********* | |
384 | ||
385 | // Prepare to handle CKM matrix elements. | |
386 | ||
387 | void VCKM::initStatic() { | |
388 | ||
389 | // Read in matrix element values and store them. | |
390 | Vsave[1][1] = Settings::parm("StandardModel:Vud"); | |
391 | Vsave[1][2] = Settings::parm("StandardModel:Vus"); | |
392 | Vsave[1][3] = Settings::parm("StandardModel:Vub"); | |
393 | Vsave[2][1] = Settings::parm("StandardModel:Vcd"); | |
394 | Vsave[2][2] = Settings::parm("StandardModel:Vcs"); | |
395 | Vsave[2][3] = Settings::parm("StandardModel:Vcb"); | |
396 | Vsave[3][1] = Settings::parm("StandardModel:Vtd"); | |
397 | Vsave[3][2] = Settings::parm("StandardModel:Vts"); | |
398 | Vsave[3][3] = Settings::parm("StandardModel:Vtb"); | |
399 | ||
400 | // Also allow for the potential existence of a fourth generation. | |
401 | Vsave[1][4] = Settings::parm("FourthGeneration:VubPrime"); | |
402 | Vsave[2][4] = Settings::parm("FourthGeneration:VcbPrime"); | |
403 | Vsave[3][4] = Settings::parm("FourthGeneration:VtbPrime"); | |
404 | Vsave[4][1] = Settings::parm("FourthGeneration:VtPrimed"); | |
405 | Vsave[4][2] = Settings::parm("FourthGeneration:VtPrimes"); | |
406 | Vsave[4][3] = Settings::parm("FourthGeneration:VtPrimeb"); | |
407 | Vsave[4][4] = Settings::parm("FourthGeneration:VtPrimebPrime"); | |
408 | ||
409 | // Calculate squares of matrix elements. | |
410 | for(int i = 1; i < 5; ++i) for(int j = 1; j < 5; ++j) | |
411 | V2save[i][j] = pow2(Vsave[i][j]); | |
412 | ||
413 | // Sum VCKM^2_out sum for given incoming flavour, excluding top as partner. | |
414 | V2out[1] = V2save[1][1] + V2save[2][1]; | |
415 | V2out[2] = V2save[1][1] + V2save[1][2] + V2save[1][3]; | |
416 | V2out[3] = V2save[1][2] + V2save[2][2]; | |
417 | V2out[4] = V2save[2][1] + V2save[2][2] + V2save[2][3]; | |
418 | V2out[5] = V2save[1][3] + V2save[2][3]; | |
419 | V2out[6] = V2save[3][1] + V2save[3][2] + V2save[3][3]; | |
420 | V2out[7] = V2save[1][4] + V2save[2][4]; | |
421 | V2out[8] = V2save[4][1] + V2save[4][2] + V2save[4][3]; | |
422 | for (int i = 11; i <= 18; ++i) V2out[i] = 1.; | |
423 | ||
424 | } | |
425 | ||
426 | //********* | |
427 | ||
428 | // Return CKM value for incoming flavours (sign irrelevant). | |
429 | ||
430 | double VCKM::Vid(int id1, int id2) { | |
431 | ||
432 | // Use absolute sign (want to cover both f -> f' W and f fbar' -> W). | |
433 | int id1Abs = abs(id1); | |
434 | int id2Abs = abs(id2); | |
435 | if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.; | |
436 | ||
437 | // Ensure proper order before reading out from Vsave or lepton match. | |
438 | if (id1Abs%2 == 1) swap(id1Abs, id2Abs); | |
439 | if (id1Abs <= 8 && id2Abs <= 8) return Vsave[id1Abs/2][(id2Abs + 1)/2]; | |
440 | if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18) | |
441 | && id2Abs == id1Abs - 1 ) return 1.; | |
442 | ||
443 | // No more valid cases. | |
444 | return 0.; | |
445 | ||
446 | } | |
447 | ||
448 | //********* | |
449 | ||
450 | // Return squared CKM value for incoming flavours (sign irrelevant). | |
451 | ||
452 | double VCKM::V2id(int id1, int id2) { | |
453 | ||
454 | // Use absolute sign (want to cover both f -> f' W and f fbar' -> W). | |
455 | int id1Abs = abs(id1); | |
456 | int id2Abs = abs(id2); | |
457 | if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.; | |
458 | ||
459 | // Ensure proper order before reading out from V2save or lepton match. | |
460 | if (id1Abs%2 == 1) swap(id1Abs, id2Abs); | |
461 | if (id1Abs <= 8 && id2Abs <= 8) return V2save[id1Abs/2][(id2Abs + 1)/2]; | |
462 | if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18) | |
463 | && id2Abs == id1Abs - 1 ) return 1.; | |
464 | ||
465 | // No more valid cases. | |
466 | return 0.; | |
467 | ||
468 | } | |
469 | ||
470 | //********* | |
471 | ||
472 | // Pick an outgoing flavour for given incoming one, given CKM mixing. | |
473 | ||
474 | int VCKM::V2pick(int id) { | |
475 | ||
476 | // Initial values. | |
477 | int idIn = abs(id); | |
478 | int idOut = 0; | |
479 | ||
480 | // Quarks: need to make random choice. | |
481 | if (idIn >= 1 && idIn <= 8) { | |
482 | double V2rndm = Rndm::flat() * V2out[idIn]; | |
483 | if (idIn == 1) idOut = (V2rndm < V2save[1][1]) ? 2 : 4; | |
484 | else if (idIn == 2) idOut = (V2rndm < V2save[1][1]) ? 1 | |
485 | : ( (V2rndm < V2save[1][1] + V2save[1][2]) ? 3 : 5 ); | |
486 | else if (idIn == 3) idOut = (V2rndm < V2save[1][2]) ? 2 : 4; | |
487 | else if (idIn == 4) idOut = (V2rndm < V2save[2][1]) ? 1 | |
488 | : ( (V2rndm < V2save[2][1] + V2save[2][2]) ? 3 : 5 ); | |
489 | else if (idIn == 5) idOut = (V2rndm < V2save[1][3]) ? 2 : 4; | |
490 | else if (idIn == 6) idOut = (V2rndm < V2save[3][1]) ? 1 | |
491 | : ( (V2rndm < V2save[3][1] + V2save[3][2]) ? 3 : 5 ); | |
492 | else if (idIn == 7) idOut = (V2rndm < V2save[1][4]) ? 2 : 4; | |
493 | else if (idIn == 8) idOut = (V2rndm < V2save[4][1]) ? 1 | |
494 | : ( (V2rndm < V2save[4][1] + V2save[4][2]) ? 3 : 5 ); | |
495 | ||
496 | // Leptons: unambiguous. | |
497 | } else if (idIn >= 11 && idIn <= 18) { | |
498 | if (idIn%2 == 1) idOut = idIn + 1; | |
499 | else idOut = idIn - 1; | |
500 | } | |
501 | ||
502 | // Done. Return with sign. | |
503 | return ( (id > 0) ? idOut : -idOut ); | |
504 | ||
505 | } | |
506 | ||
507 | //************************************************************************** | |
508 | ||
509 | } // end namespace Pythia8 |