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