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