]>
Commit | Line | Data |
---|---|---|
63ba5337 | 1 | // FragmentationFlavZpT.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 | // StringFlav, StringZ and StringPT classes. | |
8 | ||
9 | #include "FragmentationFlavZpT.h" | |
10 | ||
11 | namespace Pythia8 { | |
12 | ||
13 | //========================================================================== | |
14 | ||
15 | // The StringFlav class. | |
16 | ||
17 | //-------------------------------------------------------------------------- | |
18 | ||
19 | // Constants: could be changed here if desired, but normally should not. | |
20 | // These are of technical nature, as described for each. | |
21 | ||
22 | // Offset for different meson multiplet id values. | |
23 | const int StringFlav::mesonMultipletCode[6] | |
24 | = { 1, 3, 10003, 10001, 20003, 5}; | |
25 | ||
26 | // Clebsch-Gordan coefficients for baryon octet and decuplet are | |
27 | // fixed once and for all, so only weighted sum needs to be edited. | |
28 | // Order: ud0 + u, ud0 + s, uu1 + u, uu1 + d, ud1 + u, ud1 + s. | |
29 | const double StringFlav::baryonCGOct[6] | |
30 | = { 0.75, 0.5, 0., 0.1667, 0.0833, 0.1667}; | |
31 | const double StringFlav::baryonCGDec[6] | |
32 | = { 0., 0., 1., 0.3333, 0.6667, 0.3333}; | |
33 | ||
34 | //-------------------------------------------------------------------------- | |
35 | ||
36 | // Initialize data members of the flavour generation. | |
37 | ||
38 | void StringFlav::init(Settings& settings, Rndm* rndmPtrIn) { | |
39 | ||
40 | // Save pointer. | |
41 | rndmPtr = rndmPtrIn; | |
42 | ||
43 | // Basic parameters for generation of new flavour. | |
44 | probQQtoQ = settings.parm("StringFlav:probQQtoQ"); | |
45 | probStoUD = settings.parm("StringFlav:probStoUD"); | |
46 | probSQtoQQ = settings.parm("StringFlav:probSQtoQQ"); | |
47 | probQQ1toQQ0 = settings.parm("StringFlav:probQQ1toQQ0"); | |
48 | ||
49 | // Parameters derived from above. | |
50 | probQandQQ = 1. + probQQtoQ; | |
51 | probQandS = 2. + probStoUD; | |
52 | probQandSinQQ = 2. + probSQtoQQ * probStoUD; | |
53 | probQQ1corr = 3. * probQQ1toQQ0; | |
54 | probQQ1corrInv = 1. / probQQ1corr; | |
55 | probQQ1norm = probQQ1corr / (1. + probQQ1corr); | |
56 | ||
57 | // Parameters for normal meson production. | |
58 | for (int i = 0; i < 4; ++i) mesonRate[i][0] = 1.; | |
59 | mesonRate[0][1] = settings.parm("StringFlav:mesonUDvector"); | |
60 | mesonRate[1][1] = settings.parm("StringFlav:mesonSvector"); | |
61 | mesonRate[2][1] = settings.parm("StringFlav:mesonCvector"); | |
62 | mesonRate[3][1] = settings.parm("StringFlav:mesonBvector"); | |
63 | ||
64 | // Parameters for L=1 excited-meson production. | |
65 | mesonRate[0][2] = settings.parm("StringFlav:mesonUDL1S0J1"); | |
66 | mesonRate[1][2] = settings.parm("StringFlav:mesonSL1S0J1"); | |
67 | mesonRate[2][2] = settings.parm("StringFlav:mesonCL1S0J1"); | |
68 | mesonRate[3][2] = settings.parm("StringFlav:mesonBL1S0J1"); | |
69 | mesonRate[0][3] = settings.parm("StringFlav:mesonUDL1S1J0"); | |
70 | mesonRate[1][3] = settings.parm("StringFlav:mesonSL1S1J0"); | |
71 | mesonRate[2][3] = settings.parm("StringFlav:mesonCL1S1J0"); | |
72 | mesonRate[3][3] = settings.parm("StringFlav:mesonBL1S1J0"); | |
73 | mesonRate[0][4] = settings.parm("StringFlav:mesonUDL1S1J1"); | |
74 | mesonRate[1][4] = settings.parm("StringFlav:mesonSL1S1J1"); | |
75 | mesonRate[2][4] = settings.parm("StringFlav:mesonCL1S1J1"); | |
76 | mesonRate[3][4] = settings.parm("StringFlav:mesonBL1S1J1"); | |
77 | mesonRate[0][5] = settings.parm("StringFlav:mesonUDL1S1J2"); | |
78 | mesonRate[1][5] = settings.parm("StringFlav:mesonSL1S1J2"); | |
79 | mesonRate[2][5] = settings.parm("StringFlav:mesonCL1S1J2"); | |
80 | mesonRate[3][5] = settings.parm("StringFlav:mesonBL1S1J2"); | |
81 | ||
82 | // Store sum over multiplets for Monte Carlo generation. | |
83 | for (int i = 0; i < 4; ++i) mesonRateSum[i] | |
84 | = mesonRate[i][0] + mesonRate[i][1] + mesonRate[i][2] | |
85 | + mesonRate[i][3] + mesonRate[i][4] + mesonRate[i][5]; | |
86 | ||
87 | // Parameters for uubar - ddbar - ssbar meson mixing. | |
88 | for (int spin = 0; spin < 6; ++spin) { | |
89 | double theta; | |
90 | if (spin == 0) theta = settings.parm("StringFlav:thetaPS"); | |
91 | else if (spin == 1) theta = settings.parm("StringFlav:thetaV"); | |
92 | else if (spin == 2) theta = settings.parm("StringFlav:thetaL1S0J1"); | |
93 | else if (spin == 3) theta = settings.parm("StringFlav:thetaL1S1J0"); | |
94 | else if (spin == 4) theta = settings.parm("StringFlav:thetaL1S1J1"); | |
95 | else theta = settings.parm("StringFlav:thetaL1S1J2"); | |
96 | double alpha = (spin == 0) ? 90. - (theta + 54.7) : theta + 54.7; | |
97 | alpha *= M_PI / 180.; | |
98 | // Fill in (flavour, spin)-dependent probability of producing | |
99 | // the lightest or the lightest two mesons of the nonet. | |
100 | mesonMix1[0][spin] = 0.5; | |
101 | mesonMix2[0][spin] = 0.5 * (1. + pow2(sin(alpha))); | |
102 | mesonMix1[1][spin] = 0.; | |
103 | mesonMix2[1][spin] = pow2(cos(alpha)); | |
104 | } | |
105 | ||
106 | // Additional suppression of eta and etaPrime. | |
107 | etaSup = settings.parm("StringFlav:etaSup"); | |
108 | etaPrimeSup = settings.parm("StringFlav:etaPrimeSup"); | |
109 | ||
110 | // Sum of baryon octet and decuplet weights. | |
111 | decupletSup = settings.parm("StringFlav:decupletSup"); | |
112 | for (int i = 0; i < 6; ++i) baryonCGSum[i] | |
113 | = baryonCGOct[i] + decupletSup * baryonCGDec[i]; | |
114 | ||
115 | // Maximum SU(6) weight for ud0, ud1, uu1 types. | |
116 | baryonCGMax[0] = max( baryonCGSum[0], baryonCGSum[1]); | |
117 | baryonCGMax[1] = baryonCGMax[0]; | |
118 | baryonCGMax[2] = max( baryonCGSum[2], baryonCGSum[3]); | |
119 | baryonCGMax[3] = baryonCGMax[2]; | |
120 | baryonCGMax[4] = max( baryonCGSum[4], baryonCGSum[5]); | |
121 | baryonCGMax[5] = baryonCGMax[4]; | |
122 | ||
123 | // Popcorn baryon parameters. | |
124 | popcornRate = settings.parm("StringFlav:popcornRate"); | |
125 | popcornSpair = settings.parm("StringFlav:popcornSpair"); | |
126 | popcornSmeson = settings.parm("StringFlav:popcornSmeson"); | |
127 | ||
128 | // Suppression of leading (= first-rank) baryons. | |
129 | suppressLeadingB = settings.flag("StringFlav:suppressLeadingB"); | |
130 | lightLeadingBSup = settings.parm("StringFlav:lightLeadingBSup"); | |
131 | heavyLeadingBSup = settings.parm("StringFlav:heavyLeadingBSup"); | |
132 | ||
133 | // Begin calculation of derived parameters for baryon production. | |
134 | ||
135 | // Enumerate distinguishable diquark types (in diquark first is popcorn q). | |
136 | enum Diquark {ud0, ud1, uu1, us0, su0, us1, su1, ss1}; | |
137 | ||
138 | // Maximum SU(6) weight by diquark type. | |
139 | double barCGMax[8]; | |
140 | barCGMax[ud0] = baryonCGMax[0]; | |
141 | barCGMax[ud1] = baryonCGMax[4]; | |
142 | barCGMax[uu1] = baryonCGMax[2]; | |
143 | barCGMax[us0] = baryonCGMax[0]; | |
144 | barCGMax[su0] = baryonCGMax[0]; | |
145 | barCGMax[us1] = baryonCGMax[4]; | |
146 | barCGMax[su1] = baryonCGMax[4]; | |
147 | barCGMax[ss1] = baryonCGMax[2]; | |
148 | ||
149 | // Diquark SU(6) survival = Sum_quark (quark tunnel weight) * SU(6). | |
150 | double dMB[8]; | |
151 | dMB[ud0] = 2. * baryonCGSum[0] + probStoUD * baryonCGSum[1]; | |
152 | dMB[ud1] = 2. * baryonCGSum[4] + probStoUD * baryonCGSum[5]; | |
153 | dMB[uu1] = baryonCGSum[2] + (1. + probStoUD) * baryonCGSum[3]; | |
154 | dMB[us0] = (1. + probStoUD) * baryonCGSum[0] + baryonCGSum[1]; | |
155 | dMB[su0] = dMB[us0]; | |
156 | dMB[us1] = (1. + probStoUD) * baryonCGSum[4] + baryonCGSum[5]; | |
157 | dMB[su1] = dMB[us1]; | |
158 | dMB[ss1] = probStoUD * baryonCGSum[2] + 2. * baryonCGSum[3]; | |
159 | for (int i = 1; i < 8; ++i) dMB[i] = dMB[i] / dMB[0]; | |
160 | ||
161 | // Tunneling factors for diquark production; only half a pair = sqrt. | |
162 | double probStoUDroot = sqrt(probStoUD); | |
163 | double probSQtoQQroot = sqrt(probSQtoQQ); | |
164 | double probQQ1toQQ0root = sqrt(probQQ1toQQ0); | |
165 | double qBB[8]; | |
166 | qBB[ud1] = probQQ1toQQ0root; | |
167 | qBB[uu1] = probQQ1toQQ0root; | |
168 | qBB[us0] = probSQtoQQroot; | |
169 | qBB[su0] = probStoUDroot * probSQtoQQroot; | |
170 | qBB[us1] = probQQ1toQQ0root * qBB[us0]; | |
171 | qBB[su1] = probQQ1toQQ0root * qBB[su0]; | |
172 | qBB[ss1] = probStoUDroot * pow2(probSQtoQQroot) * probQQ1toQQ0root; | |
173 | ||
174 | // spin * (vertex factor) * (half-tunneling factor above). | |
175 | double qBM[8]; | |
176 | qBM[ud1] = 3. * qBB[ud1]; | |
177 | qBM[uu1] = 6. * qBB[uu1]; | |
178 | qBM[us0] = probStoUD * qBB[us0]; | |
179 | qBM[su0] = qBB[su0]; | |
180 | qBM[us1] = probStoUD * 3. * qBB[us1]; | |
181 | qBM[su1] = 3. * qBB[su1]; | |
182 | qBM[ss1] = probStoUD * 6. * qBB[ss1]; | |
183 | ||
184 | // Combine above two into total diquark weight for q -> B Bbar. | |
185 | for (int i = 1; i < 8; ++i) qBB[i] = qBB[i] * qBM[i]; | |
186 | ||
187 | // Suppression from having strange popcorn meson. | |
188 | qBM[us0] *= popcornSmeson; | |
189 | qBM[us1] *= popcornSmeson; | |
190 | qBM[ss1] *= popcornSmeson; | |
191 | ||
192 | // Suppression for a heavy quark of a diquark to fit into a baryon | |
193 | // on the other side of popcorn meson: (0) s/u for q -> B M; | |
194 | // (1) s/u for rank 0 diquark su -> M B; (2) ditto for s -> c/b. | |
195 | double uNorm = 1. + qBM[ud1] + qBM[uu1] + qBM[us0] + qBM[us1]; | |
196 | scbBM[0] = (2. * (qBM[su0] + qBM[su1]) + qBM[ss1]) / uNorm; | |
197 | scbBM[1] = scbBM[0] * popcornSpair * qBM[su0] / qBM[us0]; | |
198 | scbBM[2] = (1. + qBM[ud1]) * (2. + qBM[us0]) / uNorm; | |
199 | ||
200 | // Include maximum of Clebsch-Gordan coefficients. | |
201 | for (int i = 1; i < 8; ++i) dMB[i] *= qBM[i]; | |
202 | for (int i = 1; i < 8; ++i) qBM[i] *= barCGMax[i] / barCGMax[0]; | |
203 | for (int i = 1; i < 8; ++i) qBB[i] *= barCGMax[i] / barCGMax[0]; | |
204 | ||
205 | // Popcorn fraction for normal diquark production. | |
206 | double qNorm = uNorm * popcornRate / 3.; | |
207 | double sNorm = scbBM[0] * popcornSpair; | |
208 | popFrac = qNorm * (1. + qBM[ud1] + qBM[uu1] + qBM[us0] + qBM[us1] | |
209 | + sNorm * (qBM[su0] + qBM[su1] + 0.5 * qBM[ss1])) / (1. + qBB[ud1] | |
210 | + qBB[uu1] + 2. * (qBB[us0] + qBB[us1]) + 0.5 * qBB[ss1]); | |
211 | ||
212 | // Popcorn fraction for rank 0 diquarks, depending on number of s quarks. | |
213 | popS[0] = qNorm * qBM[ud1] / qBB[ud1]; | |
214 | popS[1] = qNorm * 0.5 * (qBM[us1] / qBB[us1] | |
215 | + sNorm * qBM[su1] / qBB[su1]); | |
216 | popS[2] = qNorm * sNorm * qBM[ss1] / qBB[ss1]; | |
217 | ||
218 | // Recombine diquark weights to flavour and spin ratios. Second index: | |
219 | // 0 = s/u popcorn quark ratio. | |
220 | // 1, 2 = s/u ratio for vertex quark if popcorn quark is u/d or s. | |
221 | // 3 = q/q' vertex quark ratio if popcorn quark is light and = q. | |
222 | // 4, 5, 6 = (spin 1)/(spin 0) ratio for su, us and ud. | |
223 | ||
224 | // Case 0: q -> B B. | |
225 | dWT[0][0] = (2. * (qBB[su0] + qBB[su1]) + qBB[ss1]) | |
226 | / (1. + qBB[ud1] + qBB[uu1] + qBB[us0] + qBB[us1]); | |
227 | dWT[0][1] = 2. * (qBB[us0] + qBB[us1]) / (1. + qBB[ud1] + qBB[uu1]); | |
228 | dWT[0][2] = qBB[ss1] / (qBB[su0] + qBB[su1]); | |
229 | dWT[0][3] = qBB[uu1] / (1. + qBB[ud1] + qBB[uu1]); | |
230 | dWT[0][4] = qBB[su1] / qBB[su0]; | |
231 | dWT[0][5] = qBB[us1] / qBB[us0]; | |
232 | dWT[0][6] = qBB[ud1]; | |
233 | ||
234 | // Case 1: q -> B M B. | |
235 | dWT[1][0] = (2. * (qBM[su0] + qBM[su1]) + qBM[ss1]) | |
236 | / (1. + qBM[ud1] + qBM[uu1] + qBM[us0] + qBM[us1]); | |
237 | dWT[1][1] = 2. * (qBM[us0] + qBM[us1]) / (1. + qBM[ud1] + qBM[uu1]); | |
238 | dWT[1][2] = qBM[ss1] / (qBM[su0] + qBM[su1]); | |
239 | dWT[1][3] = qBM[uu1] / (1. + qBM[ud1] + qBM[uu1]); | |
240 | dWT[1][4] = qBM[su1] / qBM[su0]; | |
241 | dWT[1][5] = qBM[us1] / qBM[us0]; | |
242 | dWT[1][6] = qBM[ud1]; | |
243 | ||
244 | // Case 2: qq -> M B; diquark inside chain. | |
245 | dWT[2][0] = (2. * (dMB[su0] + dMB[su1]) + dMB[ss1]) | |
246 | / (1. + dMB[ud1] + dMB[uu1] + dMB[us0] + dMB[us1]); | |
247 | dWT[2][1] = 2. * (dMB[us0] + dMB[us1]) / (1. + dMB[ud1] + dMB[uu1]); | |
248 | dWT[2][2] = dMB[ss1] / (dMB[su0] + dMB[su1]); | |
249 | dWT[2][3] = dMB[uu1] / (1. + dMB[ud1] + dMB[uu1]); | |
250 | dWT[2][4] = dMB[su1] / dMB[su0]; | |
251 | dWT[2][5] = dMB[us1] / dMB[us0]; | |
252 | dWT[2][6] = dMB[ud1]; | |
253 | ||
254 | } | |
255 | ||
256 | //-------------------------------------------------------------------------- | |
257 | ||
258 | // Pick a new flavour (including diquarks) given an incoming one. | |
259 | ||
260 | FlavContainer StringFlav::pick(FlavContainer& flavOld) { | |
261 | ||
262 | // Initial values for new flavour. | |
263 | FlavContainer flavNew; | |
264 | flavNew.rank = flavOld.rank + 1; | |
265 | ||
266 | // For original diquark assign popcorn quark and whether popcorn meson. | |
267 | int idOld = abs(flavOld.id); | |
268 | if (flavOld.rank == 0 && idOld > 1000) assignPopQ(flavOld); | |
269 | ||
270 | // Diquark exists, to be forced into baryon now. | |
271 | bool doOldBaryon = (idOld > 1000 && flavOld.nPop == 0); | |
272 | // Diquark exists, but do meson now. | |
273 | bool doPopcornMeson = flavOld.nPop > 0; | |
274 | // Newly created diquark gives baryon now, antibaryon later. | |
275 | bool doNewBaryon = false; | |
276 | ||
277 | // Choose whether to generate a new meson or a new baryon. | |
278 | if (!doOldBaryon && !doPopcornMeson && probQandQQ * rndmPtr->flat() > 1.) { | |
279 | doNewBaryon = true; | |
280 | if ((1. + popFrac) * rndmPtr->flat() > 1.) flavNew.nPop = 1; | |
281 | } | |
282 | ||
283 | // Optional suppression of first-rank baryon. | |
284 | if (flavOld.rank == 0 && doNewBaryon && suppressLeadingB) { | |
285 | double leadingBSup = (idOld < 4) ? lightLeadingBSup : heavyLeadingBSup; | |
286 | if (rndmPtr->flat() > leadingBSup) { | |
287 | doNewBaryon = false; | |
288 | flavNew.nPop = 0; | |
289 | } | |
290 | } | |
291 | ||
292 | // Single quark for new meson or for baryon where diquark already exists. | |
293 | if (!doPopcornMeson && !doNewBaryon) { | |
294 | flavNew.id = pickLightQ(); | |
295 | if ( (flavOld.id > 0 && flavOld.id < 9) || flavOld.id < -1000 ) | |
296 | flavNew.id = -flavNew.id; | |
297 | ||
298 | // Done for simple-quark case. | |
299 | return flavNew; | |
300 | } | |
301 | ||
302 | // Case: 0 = q -> B B, 1 = q -> B M B, 2 = qq -> M B. | |
303 | int iCase = flavNew.nPop; | |
304 | if (flavOld.nPop == 1) iCase = 2; | |
305 | ||
306 | // Flavour of popcorn quark (= q shared between B and Bbar). | |
307 | if (doNewBaryon) { | |
308 | double sPopWT = dWT[iCase][0]; | |
309 | if (iCase == 1) sPopWT *= scbBM[0] * popcornSpair; | |
310 | double rndmFlav = (2. + sPopWT) * rndmPtr->flat(); | |
311 | flavNew.idPop = 1; | |
312 | if (rndmFlav > 1.) flavNew.idPop = 2; | |
313 | if (rndmFlav > 2.) flavNew.idPop = 3; | |
314 | } else flavNew.idPop = flavOld.idPop; | |
315 | ||
316 | // Flavour of vertex quark. | |
317 | double sVtxWT = dWT[iCase][1]; | |
318 | if (flavNew.idPop >= 3) sVtxWT = dWT[iCase][2]; | |
319 | if (flavNew.idPop > 3) sVtxWT *= 0.5 * (1. + 1./dWT[iCase][4]); | |
320 | double rndmFlav = (2. + sVtxWT) * rndmPtr->flat(); | |
321 | flavNew.idVtx = 1; | |
322 | if (rndmFlav > 1.) flavNew.idVtx = 2; | |
323 | if (rndmFlav > 2.) flavNew.idVtx = 3; | |
324 | ||
325 | // Special case for light flavours, possibly identical. | |
326 | if (flavNew.idPop < 3 && flavNew.idVtx < 3) { | |
327 | flavNew.idVtx = flavNew.idPop; | |
328 | if (rndmPtr->flat() > dWT[iCase][3]) flavNew.idVtx = 3 - flavNew.idPop; | |
329 | } | |
330 | ||
331 | // Pick 2 * spin + 1. | |
332 | int spin = 3; | |
333 | if (flavNew.idVtx != flavNew.idPop) { | |
334 | double spinWT = dWT[iCase][6]; | |
335 | if (flavNew.idVtx == 3) spinWT = dWT[iCase][5]; | |
336 | if (flavNew.idPop >= 3) spinWT = dWT[iCase][4]; | |
337 | if ((1. + spinWT) * rndmPtr->flat() < 1.) spin = 1; | |
338 | } | |
339 | ||
340 | // Form outgoing diquark. Done. | |
341 | flavNew.id = 1000 * max(flavNew.idVtx, flavNew.idPop) | |
342 | + 100 * min(flavNew.idVtx, flavNew.idPop) + spin; | |
343 | if ( (flavOld.id < 0 && flavOld.id > -9) || flavOld.id > 1000 ) | |
344 | flavNew.id = -flavNew.id; | |
345 | return flavNew; | |
346 | ||
347 | } | |
348 | ||
349 | //-------------------------------------------------------------------------- | |
350 | ||
351 | // Combine two flavours (including diquarks) to produce a hadron. | |
352 | // The weighting of the combination may fail, giving output 0. | |
353 | ||
354 | int StringFlav::combine(FlavContainer& flav1, FlavContainer& flav2) { | |
355 | ||
356 | // Recognize largest and smallest flavour. | |
357 | int id1Abs = abs(flav1.id); | |
358 | int id2Abs = abs(flav2.id); | |
359 | int idMax = max(id1Abs, id2Abs); | |
360 | int idMin = min(id1Abs, id2Abs); | |
361 | ||
362 | // Construct a meson. | |
363 | if (idMax < 9 || idMin > 1000) { | |
364 | ||
365 | // Popcorn meson: use only vertex quarks. Fail if none. | |
366 | if (idMin > 1000) { | |
367 | id1Abs = flav1.idVtx; | |
368 | id2Abs = flav2.idVtx; | |
369 | idMax = max(id1Abs, id2Abs); | |
370 | idMin = min(id1Abs, id2Abs); | |
371 | if (idMin == 0) return 0; | |
372 | } | |
373 | ||
374 | // Pick spin state and preliminary code. | |
375 | int flav = (idMax < 3) ? 0 : idMax - 2; | |
376 | double rndmSpin = mesonRateSum[flav] * rndmPtr->flat(); | |
377 | int spin = -1; | |
378 | do rndmSpin -= mesonRate[flav][++spin]; | |
379 | while (rndmSpin > 0.); | |
380 | int idMeson = 100 * idMax + 10 * idMin + mesonMultipletCode[spin]; | |
381 | ||
382 | // For nondiagonal mesons distinguish particle/antiparticle. | |
383 | if (idMax != idMin) { | |
384 | int sign = (idMax%2 == 0) ? 1 : -1; | |
385 | if ( (idMax == id1Abs && flav1.id < 0) | |
386 | || (idMax == id2Abs && flav2.id < 0) ) sign = -sign; | |
387 | idMeson *= sign; | |
388 | ||
389 | // For light diagonal mesons include uubar - ddbar - ssbar mixing. | |
390 | } else if (flav < 2) { | |
391 | double rMix = rndmPtr->flat(); | |
392 | if (rMix < mesonMix1[flav][spin]) idMeson = 110; | |
393 | else if (rMix < mesonMix2[flav][spin]) idMeson = 220; | |
394 | else idMeson = 330; | |
395 | idMeson += mesonMultipletCode[spin]; | |
396 | ||
397 | // Additional suppression of eta and eta' may give failure. | |
398 | if (idMeson == 221 && etaSup < rndmPtr->flat()) return 0; | |
399 | if (idMeson == 331 && etaPrimeSup < rndmPtr->flat()) return 0; | |
400 | } | |
401 | ||
402 | // Finished for mesons. | |
403 | return idMeson; | |
404 | } | |
405 | ||
406 | // SU(6) factors for baryon production may give failure. | |
407 | int idQQ1 = idMax / 1000; | |
408 | int idQQ2 = (idMax / 100) % 10; | |
409 | int spinQQ = idMax % 10; | |
410 | int spinFlav = spinQQ - 1; | |
411 | if (spinFlav == 2 && idQQ1 != idQQ2) spinFlav = 4; | |
412 | if (idMin != idQQ1 && idMin != idQQ2) spinFlav++; | |
413 | if (baryonCGSum[spinFlav] < rndmPtr->flat() * baryonCGMax[spinFlav]) | |
414 | return 0; | |
415 | ||
416 | // Order quarks to form baryon. Pick spin. | |
417 | int idOrd1 = max( idMin, max( idQQ1, idQQ2) ); | |
418 | int idOrd3 = min( idMin, min( idQQ1, idQQ2) ); | |
419 | int idOrd2 = idMin + idQQ1 + idQQ2 - idOrd1 - idOrd3; | |
420 | int spinBar = (baryonCGSum[spinFlav] * rndmPtr->flat() | |
421 | < baryonCGOct[spinFlav]) ? 2 : 4; | |
422 | ||
423 | // Distinguish Lambda- and Sigma-like. | |
424 | bool LambdaLike = false; | |
425 | if (spinBar == 2 && idOrd1 > idOrd2 && idOrd2 > idOrd3) { | |
426 | LambdaLike = (spinQQ == 1); | |
427 | if (idOrd1 != idMin && spinQQ == 1) LambdaLike = (rndmPtr->flat() < 0.25); | |
428 | else if (idOrd1 != idMin) LambdaLike = (rndmPtr->flat() < 0.75); | |
429 | } | |
430 | ||
431 | // Form baryon code and return with sign. | |
432 | int idBaryon = (LambdaLike) | |
433 | ? 1000 * idOrd1 + 100 * idOrd3 + 10 * idOrd2 + spinBar | |
434 | : 1000 * idOrd1 + 100 * idOrd2 + 10 * idOrd3 + spinBar; | |
435 | return (flav1.id > 0) ? idBaryon : -idBaryon; | |
436 | ||
437 | } | |
438 | ||
439 | //-------------------------------------------------------------------------- | |
440 | ||
441 | // Assign popcorn quark inside an original (= rank 0) diquark. | |
442 | ||
443 | void StringFlav::assignPopQ(FlavContainer& flav) { | |
444 | ||
445 | // Safety check that intended to do something. | |
446 | int idAbs = abs(flav.id); | |
447 | if (flav.rank > 0 || idAbs < 1000) return; | |
448 | ||
449 | // Make choice of popcorn quark. | |
450 | int id1 = (idAbs/1000)%10; | |
451 | int id2 = (idAbs/100)%10; | |
452 | double pop2WT = 1.; | |
453 | if (id1 == 3) pop2WT = scbBM[1]; | |
454 | else if (id1 > 3) pop2WT = scbBM[2]; | |
455 | if (id2 == 3) pop2WT /= scbBM[1]; | |
456 | else if (id2 > 3) pop2WT /= scbBM[2]; | |
457 | // Agrees with Patrik code, but opposite to intention?? | |
458 | flav.idPop = ((1. + pop2WT) * rndmPtr->flat() > 1.) ? id2 : id1; | |
459 | flav.idVtx = id1 + id2 - flav.idPop; | |
460 | ||
461 | // Also determine if to produce popcorn meson. | |
462 | flav.nPop = 0; | |
463 | double popWT = popS[0]; | |
464 | if (id1 == 3) popWT = popS[1]; | |
465 | if (id2 == 3) popWT = popS[2]; | |
466 | if (idAbs%10 == 1) popWT *= sqrt(probQQ1toQQ0); | |
467 | if ((1. + popWT) * rndmPtr->flat() > 1.) flav.nPop = 1; | |
468 | ||
469 | } | |
470 | ||
471 | //-------------------------------------------------------------------------- | |
472 | ||
473 | // Combine two quarks to produce a diquark. | |
474 | // Normally according to production composition, but nonvanishing idHad | |
475 | // means diquark from known hadron content, so use SU(6) wave fucntion. | |
476 | ||
477 | int StringFlav::makeDiquark(int id1, int id2, int idHad) { | |
478 | ||
479 | // Initial values. | |
480 | int idMin = min( abs(id1), abs(id2)); | |
481 | int idMax = max( abs(id1), abs(id2)); | |
482 | int spin = 1; | |
483 | ||
484 | // Select spin of diquark formed from two valence quarks in proton. | |
485 | // (More hadron cases??) | |
486 | if (abs(idHad) == 2212) { | |
487 | if (idMin == 1 && idMax == 2 && rndmPtr->flat() < 0.75) spin = 0; | |
488 | ||
489 | // Else select spin of diquark according to production composition. | |
490 | } else { | |
491 | if (idMin != idMax && rndmPtr->flat() > probQQ1norm) spin = 0; | |
492 | } | |
493 | ||
494 | // Combined diquark code. | |
495 | int idNewAbs = 1000 * idMax + 100 * idMin + 2 * spin + 1; | |
496 | return (id1 > 0) ? idNewAbs : -idNewAbs; | |
497 | ||
498 | } | |
499 | ||
500 | //========================================================================== | |
501 | ||
502 | // The StringZ class. | |
503 | ||
504 | //-------------------------------------------------------------------------- | |
505 | ||
506 | // Constants: could be changed here if desired, but normally should not. | |
507 | // These are of technical nature, as described for each. | |
508 | ||
509 | // When a or c are close to special cases, default to these. | |
510 | const double StringZ::CFROMUNITY = 0.01; | |
511 | const double StringZ::AFROMZERO = 0.02; | |
512 | const double StringZ::AFROMC = 0.01; | |
513 | ||
514 | // Do not take exponent of too large or small number. | |
515 | const double StringZ::EXPMAX = 50.; | |
516 | ||
517 | //-------------------------------------------------------------------------- | |
518 | ||
519 | // Initialize data members of the string z selection. | |
520 | ||
521 | void StringZ::init(Settings& settings, ParticleData& particleData, | |
522 | Rndm* rndmPtrIn) { | |
523 | ||
524 | // Save pointer. | |
525 | rndmPtr = rndmPtrIn; | |
526 | ||
527 | // c and b quark masses. | |
528 | mc2 = pow2( particleData.m0(4)); | |
529 | mb2 = pow2( particleData.m0(5)); | |
530 | ||
531 | // Paramaters of Lund/Bowler symmetric fragmentation function. | |
532 | aLund = settings.parm("StringZ:aLund"); | |
533 | bLund = settings.parm("StringZ:bLund"); | |
534 | aExtraDiquark = settings.parm("StringZ:aExtraDiquark"); | |
535 | rFactC = settings.parm("StringZ:rFactC"); | |
536 | rFactB = settings.parm("StringZ:rFactB"); | |
537 | rFactH = settings.parm("StringZ:rFactH"); | |
538 | ||
539 | // Flags and parameters of Peterson/SLAC fragmentation function. | |
540 | usePetersonC = settings.flag("StringZ:usePetersonC"); | |
541 | usePetersonB = settings.flag("StringZ:usePetersonB"); | |
542 | usePetersonH = settings.flag("StringZ:usePetersonH"); | |
543 | epsilonC = settings.parm("StringZ:epsilonC"); | |
544 | epsilonB = settings.parm("StringZ:epsilonB"); | |
545 | epsilonH = settings.parm("StringZ:epsilonH"); | |
546 | ||
547 | // Parameters for joining procedure. | |
548 | stopM = settings.parm("StringFragmentation:stopMass"); | |
549 | stopNF = settings.parm("StringFragmentation:stopNewFlav"); | |
550 | stopS = settings.parm("StringFragmentation:stopSmear"); | |
551 | ||
552 | } | |
553 | ||
554 | //-------------------------------------------------------------------------- | |
555 | ||
556 | // Generate the fraction z that the next hadron will take, | |
557 | // using either Lund/Bowler or, for heavy, Peterson/SLAC functions. | |
558 | // Note: for a heavy new coloured particle we assume pT negligible. | |
559 | ||
560 | double StringZ::zFrag( int idOld, int idNew, double mT2) { | |
561 | ||
562 | // Find if old or new flavours correspond to diquarks. | |
563 | int idOldAbs = abs(idOld); | |
564 | int idNewAbs = abs(idNew); | |
565 | bool isOldDiquark = (idOldAbs > 1000 && idOldAbs < 10000); | |
566 | bool isNewDiquark = (idNewAbs > 1000 && idNewAbs < 10000); | |
567 | ||
568 | // Find heaviest quark in fragmenting parton/diquark. | |
569 | int idFrag = idOldAbs; | |
570 | if (isOldDiquark) idFrag = max( idOldAbs / 1000, (idOldAbs / 100) % 10); | |
571 | ||
572 | // Use Peterson where explicitly requested for heavy flavours. | |
573 | if (idFrag == 4 && usePetersonC) return zPeterson( epsilonC); | |
574 | if (idFrag == 5 && usePetersonB) return zPeterson( epsilonB); | |
575 | if (idFrag > 5 && usePetersonH) { | |
576 | double epsilon = epsilonH * mb2 / mT2; | |
577 | return zPeterson( epsilon); | |
578 | } | |
579 | ||
580 | // Shape parameters of Lund symmetric fragmentation function. | |
581 | double aShape = aLund; | |
582 | if (isOldDiquark) aShape += aExtraDiquark; | |
583 | double bShape = bLund * mT2; | |
584 | double cShape = 1.; | |
585 | if (isOldDiquark) cShape -= aExtraDiquark; | |
586 | if (isNewDiquark) cShape += aExtraDiquark; | |
587 | if (idFrag == 4) cShape += rFactC * bLund * mc2; | |
588 | if (idFrag == 5) cShape += rFactB * bLund * mb2; | |
589 | if (idFrag > 5) cShape += rFactH * bLund * mT2; | |
590 | return zLund( aShape, bShape, cShape); | |
591 | ||
592 | } | |
593 | ||
594 | //-------------------------------------------------------------------------- | |
595 | ||
596 | // Generate a random z according to the Lund/Bowler symmetric | |
597 | // fragmentation function f(z) = (1 -z)^a * exp(-b/z) / z^c. | |
598 | // Normalized so that f(z_max) = 1 it can also be written as | |
599 | // f(z) = exp( a * ln( (1 - z) / (1 - z_max) ) + b * (1/z_max - 1/z) | |
600 | // + c * ln(z_max/z) ). | |
601 | ||
602 | double StringZ::zLund( double a, double b, double c) { | |
603 | ||
604 | // Special cases for c = 1, a = 0 and a = c. | |
605 | bool cIsUnity = (abs( c - 1.) < CFROMUNITY); | |
606 | bool aIsZero = (a < AFROMZERO); | |
607 | bool aIsC = (abs(a - c) < AFROMC); | |
608 | ||
609 | // Determine position of maximum. | |
610 | double zMax; | |
611 | if (aIsZero) zMax = (c > b) ? b / c : 1.; | |
612 | else if (aIsC) zMax = b / (b + c); | |
613 | else { zMax = 0.5 * (b + c - sqrt( pow2(b - c) + 4. * a * b)) / (c - a); | |
614 | if (zMax > 0.9999 && b > 100.) zMax = min(zMax, 1. - a / b); } | |
615 | ||
616 | // Subdivide z range if distribution very peaked near either endpoint. | |
617 | bool peakedNearZero = (zMax < 0.1); | |
618 | bool peakedNearUnity = (zMax > 0.85 && b > 1.); | |
619 | ||
620 | // Find integral of trial function everywhere bigger than f. | |
621 | // (Dummy start values.) | |
622 | double fIntLow = 1.; | |
623 | double fIntHigh = 1.; | |
624 | double fInt = 2.; | |
625 | double zDiv = 0.5; | |
626 | double zDivC = 0.5; | |
627 | // When z_max is small use that f(z) | |
628 | // < 1 for z < z_div = 2.75 * z_max, | |
629 | // < (z_div/z)^c for z > z_div (=> logarithm for c = 1, else power). | |
630 | if (peakedNearZero) { | |
631 | zDiv = 2.75 * zMax; | |
632 | fIntLow = zDiv; | |
633 | if (cIsUnity) fIntHigh = -zDiv * log(zDiv); | |
634 | else { zDivC = pow( zDiv, 1. - c); | |
635 | fIntHigh = zDiv * (1. - 1./zDivC) / (c - 1.);} | |
636 | fInt = fIntLow + fIntHigh; | |
637 | // When z_max large use that f(z) | |
638 | // < exp( b * (z - z_div) ) for z < z_div with z_div messy expression, | |
639 | // < 1 for z > z_div. | |
640 | // To simplify expressions the integral is extended to z = -infinity. | |
641 | } else if (peakedNearUnity) { | |
642 | double rcb = sqrt(4. + pow2(c / b)); | |
643 | zDiv = rcb - 1./zMax - (c / b) * log( zMax * 0.5 * (rcb + c / b) ); | |
644 | if (!aIsZero) zDiv += (a/b) * log(1. - zMax); | |
645 | zDiv = min( zMax, max(0., zDiv)); | |
646 | fIntLow = 1. / b; | |
647 | fIntHigh = 1. - zDiv; | |
648 | fInt = fIntLow + fIntHigh; | |
649 | } | |
650 | ||
651 | // Choice of z, preweighted for peaks at low or high z. (Dummy start values.) | |
652 | double z = 0.5; | |
653 | double fPrel = 1.; | |
654 | double fVal = 1.; | |
655 | do { | |
656 | // Choice of z flat good enough for distribution peaked in the middle; | |
657 | // if not this z can be reused as a random number in general. | |
658 | z = rndmPtr->flat(); | |
659 | fPrel = 1.; | |
660 | // When z_max small use flat below z_div and 1/z^c above z_div. | |
661 | if (peakedNearZero) { | |
662 | if (fInt * rndmPtr->flat() < fIntLow) z = zDiv * z; | |
663 | else if (cIsUnity) {z = pow( zDiv, z); fPrel = zDiv / z;} | |
664 | else { z = pow( zDivC + (1. - zDivC) * z, 1. / (1. - c) ); | |
665 | fPrel = pow( zDiv / z, c); } | |
666 | // When z_max large use exp( b * (z -z_div) ) below z_div | |
667 | // and flat above it. | |
668 | } else if (peakedNearUnity) { | |
669 | if (fInt * rndmPtr->flat() < fIntLow) { | |
670 | z = zDiv + log(z) / b; | |
671 | fPrel = exp( b * (z - zDiv) ); | |
672 | } else z = zDiv + (1. - zDiv) * z; | |
673 | } | |
674 | ||
675 | // Evaluate actual f(z) (if in physical range) and correct. | |
676 | if (z > 0 && z < 1) { | |
677 | double fExp = b * (1. / zMax - 1. / z)+ c * log(zMax / z); | |
678 | if (!aIsZero) fExp += a * log( (1. - z) / (1. - zMax) ); | |
679 | fVal = exp( max( -EXPMAX, min( EXPMAX, fExp) ) ) ; | |
680 | } else fVal = 0.; | |
681 | } while (fVal < rndmPtr->flat() * fPrel); | |
682 | ||
683 | // Done. | |
684 | return z; | |
685 | ||
686 | } | |
687 | ||
688 | //-------------------------------------------------------------------------- | |
689 | ||
690 | // Generate a random z according to the Peterson/SLAC formula | |
691 | // f(z) = 1 / ( z * (1 - 1/z - epsilon/(1-z))^2 ) | |
692 | // = z * (1-z)^2 / ((1-z)^2 + epsilon * z)^2. | |
693 | ||
694 | double StringZ::zPeterson( double epsilon) { | |
695 | ||
696 | double z, fVal; | |
697 | ||
698 | // For large epsilon pick z flat and reject, | |
699 | // knowing that 4 * epsilon * f(z) < 1 everywhere. | |
700 | if (epsilon > 0.01) { | |
701 | do { | |
702 | z = rndmPtr->flat(); | |
703 | fVal = 4. * epsilon * z * pow2(1. - z) | |
704 | / pow2( pow2(1. - z) + epsilon * z); | |
705 | } while (fVal < rndmPtr->flat()); | |
706 | return z; | |
707 | } | |
708 | ||
709 | // Else split range, using that 4 * epsilon * f(z) | |
710 | // < 4 * epsilon / (1 - z)^2 for 0 < z < 1 - 2 * sqrt(epsilon) | |
711 | // < 1 for 1 - 2 * sqrt(epsilon) < z < 1 | |
712 | double epsRoot = sqrt(epsilon); | |
713 | double epsComb = 0.5 / epsRoot - 1.; | |
714 | double fIntLow = 4. * epsilon * epsComb; | |
715 | double fInt = fIntLow + 2. * epsRoot; | |
716 | do { | |
717 | if (rndmPtr->flat() * fInt < fIntLow) { | |
718 | z = 1. - 1. / (1. + rndmPtr->flat() * epsComb); | |
719 | fVal = z * pow2( pow2(1. - z) / (pow2(1. - z) + epsilon * z) ); | |
720 | } else { | |
721 | z = 1. - 2. * epsRoot * rndmPtr->flat(); | |
722 | fVal = 4. * epsilon * z * pow2(1. - z) | |
723 | / pow2( pow2(1. - z) + epsilon * z); | |
724 | } | |
725 | } while (fVal < rndmPtr->flat()); | |
726 | return z; | |
727 | ||
728 | } | |
729 | ||
730 | //========================================================================== | |
731 | ||
732 | // The StringPT class. | |
733 | ||
734 | //-------------------------------------------------------------------------- | |
735 | ||
736 | // Constants: could be changed here if desired, but normally should not. | |
737 | // These are of technical nature, as described for each. | |
738 | ||
739 | // To avoid division by zero one must have sigma > 0. | |
740 | const double StringPT::SIGMAMIN = 0.2; | |
741 | ||
742 | //-------------------------------------------------------------------------- | |
743 | ||
744 | // Initialize data members of the string pT selection. | |
745 | ||
746 | void StringPT::init(Settings& settings, ParticleData& , Rndm* rndmPtrIn) { | |
747 | ||
748 | // Save pointer. | |
749 | rndmPtr = rndmPtrIn; | |
750 | ||
751 | // Parameters of the pT width and enhancement. | |
752 | double sigma = settings.parm("StringPT:sigma"); | |
753 | sigmaQ = sigma / sqrt(2.); | |
754 | enhancedFraction = settings.parm("StringPT:enhancedFraction"); | |
755 | enhancedWidth = settings.parm("StringPT:enhancedWidth"); | |
756 | ||
757 | // Parameter for pT suppression in MiniStringFragmentation. | |
758 | sigma2Had = 2. * pow2( max( SIGMAMIN, sigma) ); | |
759 | ||
760 | } | |
761 | ||
762 | //-------------------------------------------------------------------------- | |
763 | ||
764 | // Generate Gaussian pT such that <p_x^2> = <p_x^2> = sigma^2 = width^2/2, | |
765 | // but with small fraction multiplied up to a broader spectrum. | |
766 | ||
767 | pair<double, double> StringPT::pxy() { | |
768 | ||
769 | double sigma = sigmaQ; | |
770 | if (rndmPtr->flat() < enhancedFraction) sigma *= enhancedWidth; | |
771 | pair<double, double> gauss2 = rndmPtr->gauss2(); | |
772 | return pair<double, double>(sigma * gauss2.first, sigma * gauss2.second); | |
773 | ||
774 | } | |
775 | ||
776 | //========================================================================== | |
777 | ||
778 | } // end namespace Pythia8 |