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.
6
// Function definitions (not found in the header) for the
7
// StringFlav, StringZ and StringPT classes.
9
#include "FragmentationFlavZpT.h"
13
//==========================================================================
15
// The StringFlav class.
17
//--------------------------------------------------------------------------
19
// Constants: could be changed here if desired, but normally should not.
20
// These are of technical nature, as described for each.
22
// Offset for different meson multiplet id values.
23
const int StringFlav::mesonMultipletCode[6]
24
= { 1, 3, 10003, 10001, 20003, 5};
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};
34
//--------------------------------------------------------------------------
36
// Initialize data members of the flavour generation.
38
void StringFlav::init(Settings& settings, Rndm* rndmPtrIn) {
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");
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);
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");
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");
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];
87
// Parameters for uubar - ddbar - ssbar meson mixing.
88
for (int spin = 0; spin < 6; ++spin) {
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;
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));
106
// Additional suppression of eta and etaPrime.
107
etaSup = settings.parm("StringFlav:etaSup");
108
etaPrimeSup = settings.parm("StringFlav:etaPrimeSup");
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];
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];
123
// Popcorn baryon parameters.
124
popcornRate = settings.parm("StringFlav:popcornRate");
125
popcornSpair = settings.parm("StringFlav:popcornSpair");
126
popcornSmeson = settings.parm("StringFlav:popcornSmeson");
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");
133
// Begin calculation of derived parameters for baryon production.
135
// Enumerate distinguishable diquark types (in diquark first is popcorn q).
136
enum Diquark {ud0, ud1, uu1, us0, su0, us1, su1, ss1};
138
// Maximum SU(6) weight by diquark type.
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];
149
// Diquark SU(6) survival = Sum_quark (quark tunnel weight) * SU(6).
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];
156
dMB[us1] = (1. + probStoUD) * baryonCGSum[4] + baryonCGSum[5];
158
dMB[ss1] = probStoUD * baryonCGSum[2] + 2. * baryonCGSum[3];
159
for (int i = 1; i < 8; ++i) dMB[i] = dMB[i] / dMB[0];
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);
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;
174
// spin * (vertex factor) * (half-tunneling factor above).
176
qBM[ud1] = 3. * qBB[ud1];
177
qBM[uu1] = 6. * qBB[uu1];
178
qBM[us0] = probStoUD * qBB[us0];
180
qBM[us1] = probStoUD * 3. * qBB[us1];
181
qBM[su1] = 3. * qBB[su1];
182
qBM[ss1] = probStoUD * 6. * qBB[ss1];
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];
187
// Suppression from having strange popcorn meson.
188
qBM[us0] *= popcornSmeson;
189
qBM[us1] *= popcornSmeson;
190
qBM[ss1] *= popcornSmeson;
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;
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];
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]);
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];
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.
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];
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];
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];
256
//--------------------------------------------------------------------------
258
// Pick a new flavour (including diquarks) given an incoming one.
260
FlavContainer StringFlav::pick(FlavContainer& flavOld) {
262
// Initial values for new flavour.
263
FlavContainer flavNew;
264
flavNew.rank = flavOld.rank + 1;
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);
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;
277
// Choose whether to generate a new meson or a new baryon.
278
if (!doOldBaryon && !doPopcornMeson && probQandQQ * rndmPtr->flat() > 1.) {
280
if ((1. + popFrac) * rndmPtr->flat() > 1.) flavNew.nPop = 1;
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) {
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;
298
// Done for simple-quark case.
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;
306
// Flavour of popcorn quark (= q shared between B and Bbar).
308
double sPopWT = dWT[iCase][0];
309
if (iCase == 1) sPopWT *= scbBM[0] * popcornSpair;
310
double rndmFlav = (2. + sPopWT) * rndmPtr->flat();
312
if (rndmFlav > 1.) flavNew.idPop = 2;
313
if (rndmFlav > 2.) flavNew.idPop = 3;
314
} else flavNew.idPop = flavOld.idPop;
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();
322
if (rndmFlav > 1.) flavNew.idVtx = 2;
323
if (rndmFlav > 2.) flavNew.idVtx = 3;
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;
331
// Pick 2 * spin + 1.
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;
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;
349
//--------------------------------------------------------------------------
351
// Combine two flavours (including diquarks) to produce a hadron.
352
// The weighting of the combination may fail, giving output 0.
354
int StringFlav::combine(FlavContainer& flav1, FlavContainer& flav2) {
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);
362
// Construct a meson.
363
if (idMax < 9 || idMin > 1000) {
365
// Popcorn meson: use only vertex quarks. Fail if none.
367
id1Abs = flav1.idVtx;
368
id2Abs = flav2.idVtx;
369
idMax = max(id1Abs, id2Abs);
370
idMin = min(id1Abs, id2Abs);
371
if (idMin == 0) return 0;
374
// Pick spin state and preliminary code.
375
int flav = (idMax < 3) ? 0 : idMax - 2;
376
double rndmSpin = mesonRateSum[flav] * rndmPtr->flat();
378
do rndmSpin -= mesonRate[flav][++spin];
379
while (rndmSpin > 0.);
380
int idMeson = 100 * idMax + 10 * idMin + mesonMultipletCode[spin];
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;
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;
395
idMeson += mesonMultipletCode[spin];
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;
402
// Finished for mesons.
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])
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;
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);
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;
439
//--------------------------------------------------------------------------
441
// Assign popcorn quark inside an original (= rank 0) diquark.
443
void StringFlav::assignPopQ(FlavContainer& flav) {
445
// Safety check that intended to do something.
446
int idAbs = abs(flav.id);
447
if (flav.rank > 0 || idAbs < 1000) return;
449
// Make choice of popcorn quark.
450
int id1 = (idAbs/1000)%10;
451
int id2 = (idAbs/100)%10;
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;
461
// Also determine if to produce popcorn meson.
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;
471
//--------------------------------------------------------------------------
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.
477
int StringFlav::makeDiquark(int id1, int id2, int idHad) {
480
int idMin = min( abs(id1), abs(id2));
481
int idMax = max( abs(id1), abs(id2));
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;
489
// Else select spin of diquark according to production composition.
491
if (idMin != idMax && rndmPtr->flat() > probQQ1norm) spin = 0;
494
// Combined diquark code.
495
int idNewAbs = 1000 * idMax + 100 * idMin + 2 * spin + 1;
496
return (id1 > 0) ? idNewAbs : -idNewAbs;
500
//==========================================================================
502
// The StringZ class.
504
//--------------------------------------------------------------------------
506
// Constants: could be changed here if desired, but normally should not.
507
// These are of technical nature, as described for each.
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;
514
// Do not take exponent of too large or small number.
515
const double StringZ::EXPMAX = 50.;
517
//--------------------------------------------------------------------------
519
// Initialize data members of the string z selection.
521
void StringZ::init(Settings& settings, ParticleData& particleData,
527
// c and b quark masses.
528
mc2 = pow2( particleData.m0(4));
529
mb2 = pow2( particleData.m0(5));
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");
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");
547
// Parameters for joining procedure.
548
stopM = settings.parm("StringFragmentation:stopMass");
549
stopNF = settings.parm("StringFragmentation:stopNewFlav");
550
stopS = settings.parm("StringFragmentation:stopSmear");
554
//--------------------------------------------------------------------------
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.
560
double StringZ::zFrag( int idOld, int idNew, double mT2) {
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);
568
// Find heaviest quark in fragmenting parton/diquark.
569
int idFrag = idOldAbs;
570
if (isOldDiquark) idFrag = max( idOldAbs / 1000, (idOldAbs / 100) % 10);
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);
580
// Shape parameters of Lund symmetric fragmentation function.
581
double aShape = aLund;
582
if (isOldDiquark) aShape += aExtraDiquark;
583
double bShape = bLund * mT2;
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);
594
//--------------------------------------------------------------------------
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) ).
602
double StringZ::zLund( double a, double b, double c) {
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);
609
// Determine position of maximum.
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); }
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.);
620
// Find integral of trial function everywhere bigger than f.
621
// (Dummy start values.)
623
double fIntHigh = 1.;
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) {
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));
647
fIntHigh = 1. - zDiv;
648
fInt = fIntLow + fIntHigh;
651
// Choice of z, preweighted for peaks at low or high z. (Dummy start values.)
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.
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;
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) ) ) ;
681
} while (fVal < rndmPtr->flat() * fPrel);
688
//--------------------------------------------------------------------------
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.
694
double StringZ::zPeterson( double epsilon) {
698
// For large epsilon pick z flat and reject,
699
// knowing that 4 * epsilon * f(z) < 1 everywhere.
700
if (epsilon > 0.01) {
703
fVal = 4. * epsilon * z * pow2(1. - z)
704
/ pow2( pow2(1. - z) + epsilon * z);
705
} while (fVal < rndmPtr->flat());
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;
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) );
721
z = 1. - 2. * epsRoot * rndmPtr->flat();
722
fVal = 4. * epsilon * z * pow2(1. - z)
723
/ pow2( pow2(1. - z) + epsilon * z);
725
} while (fVal < rndmPtr->flat());
730
//==========================================================================
732
// The StringPT class.
734
//--------------------------------------------------------------------------
736
// Constants: could be changed here if desired, but normally should not.
737
// These are of technical nature, as described for each.
739
// To avoid division by zero one must have sigma > 0.
740
const double StringPT::SIGMAMIN = 0.2;
742
//--------------------------------------------------------------------------
744
// Initialize data members of the string pT selection.
746
void StringPT::init(Settings& settings, ParticleData& , Rndm* rndmPtrIn) {
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");
757
// Parameter for pT suppression in MiniStringFragmentation.
758
sigma2Had = 2. * pow2( max( SIGMAMIN, sigma) );
762
//--------------------------------------------------------------------------
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.
767
pair<double, double> StringPT::pxy() {
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);
776
//==========================================================================
778
} // end namespace Pythia8
added
removed
src/FragmentationFlavZpT.cc
