StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
FragmentationFlavZpT.cc
1 // FragmentationFlavZpT.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 // Functions for unnormalised and average Lund FF.
16 
17 //--------------------------------------------------------------------------
18 
19 // The unnormalised Lund FF
20 
21 double LundFFRaw(double z, double a, double b, double c, double mT2) {
22  if (z <= 0. || z >= 1.) return 0.;
23  return pow(1. - z, a) / pow(z, c) * exp(-b * mT2 / z);
24 }
25 
26 //--------------------------------------------------------------------------
27 
28 // Average, <z>, of Lund FF.
29 
30 double LundFFAvg(double a, double b, double c,
31  double mT2, double tol = 1.e-6) {
32 
33  // Checks whether the integration succeeded
34  bool check;
35 
36  // Define lundFF as a function of only z, fixing a, b, c, mT2 as parameters
37  // Note that c must be captured by reference, since it is modified later
38  auto lundFF = [=, &c](double z) { return LundFFRaw(z, a, b, c, mT2); };
39 
40  // Get denominator
41  double denominator = 1.;
42  check = integrateGauss(denominator, lundFF, 0, 1, tol);
43  if (!check || denominator <= 0.) return -1.;
44 
45  // Get numerator
46  c -= 1;
47  double numerator = 0.;
48  check = integrateGauss(numerator, lundFF, 0., 1., tol);
49  if (!check || numerator <= 0.) return -1.;
50 
51  // Done.
52  return numerator / denominator;
53 }
54 
55 //==========================================================================
56 
57 // The StringFlav class.
58 
59 //--------------------------------------------------------------------------
60 
61 // Constants: could be changed here if desired, but normally should not.
62 // These are of technical nature, as described for each.
63 
64 // Offset for different meson multiplet id values.
65 const int StringFlav::mesonMultipletCode[6]
66  = { 1, 3, 10003, 10001, 20003, 5};
67 
68 // Clebsch-Gordan coefficients for baryon octet and decuplet are
69 // fixed once and for all, so only weighted sum needs to be edited.
70 // Order: ud0 + u, ud0 + s, uu1 + u, uu1 + d, ud1 + u, ud1 + s.
71 const double StringFlav::baryonCGOct[6]
72  = { 0.75, 0.5, 0., 0.1667, 0.0833, 0.1667};
73 const double StringFlav::baryonCGDec[6]
74  = { 0., 0., 1., 0.3333, 0.6667, 0.3333};
75 
76 //--------------------------------------------------------------------------
77 
78 // Initialize data members of the flavour generation.
79 
80 void StringFlav::init() {
81 
82  // Save pointers.
83  // Basic parameters for generation of new flavour.
84  probQQtoQ = parm("StringFlav:probQQtoQ");
85  probStoUD = parm("StringFlav:probStoUD");
86  probSQtoQQ = parm("StringFlav:probSQtoQQ");
87  probQQ1toQQ0 = parm("StringFlav:probQQ1toQQ0");
88 
89  // Parameters derived from above.
90  probQandQQ = 1. + probQQtoQ;
91  probQandS = 2. + probStoUD;
92  probQandSinQQ = 2. + probSQtoQQ * probStoUD;
93  probQQ1corr = 3. * probQQ1toQQ0;
94  probQQ1corrInv = 1. / probQQ1corr;
95  probQQ1norm = probQQ1corr / (1. + probQQ1corr);
96 
97  // Spin parameters for combining two quarks to a diquark.
98  vector<double> pQQ1tmp = settingsPtr->pvec("StringFlav:probQQ1toQQ0join");
99  for (int i = 0; i < 4; ++i)
100  probQQ1join[i] = 3. * pQQ1tmp[i] / (1. + 3. * pQQ1tmp[i]);
101 
102  // Parameters for normal meson production.
103  for (int i = 0; i < 4; ++i) mesonRate[i][0] = 1.;
104  mesonRate[0][1] = parm("StringFlav:mesonUDvector");
105  mesonRate[1][1] = parm("StringFlav:mesonSvector");
106  mesonRate[2][1] = parm("StringFlav:mesonCvector");
107  mesonRate[3][1] = parm("StringFlav:mesonBvector");
108 
109  // Parameters for L=1 excited-meson production.
110  mesonRate[0][2] = parm("StringFlav:mesonUDL1S0J1");
111  mesonRate[1][2] = parm("StringFlav:mesonSL1S0J1");
112  mesonRate[2][2] = parm("StringFlav:mesonCL1S0J1");
113  mesonRate[3][2] = parm("StringFlav:mesonBL1S0J1");
114  mesonRate[0][3] = parm("StringFlav:mesonUDL1S1J0");
115  mesonRate[1][3] = parm("StringFlav:mesonSL1S1J0");
116  mesonRate[2][3] = parm("StringFlav:mesonCL1S1J0");
117  mesonRate[3][3] = parm("StringFlav:mesonBL1S1J0");
118  mesonRate[0][4] = parm("StringFlav:mesonUDL1S1J1");
119  mesonRate[1][4] = parm("StringFlav:mesonSL1S1J1");
120  mesonRate[2][4] = parm("StringFlav:mesonCL1S1J1");
121  mesonRate[3][4] = parm("StringFlav:mesonBL1S1J1");
122  mesonRate[0][5] = parm("StringFlav:mesonUDL1S1J2");
123  mesonRate[1][5] = parm("StringFlav:mesonSL1S1J2");
124  mesonRate[2][5] = parm("StringFlav:mesonCL1S1J2");
125  mesonRate[3][5] = parm("StringFlav:mesonBL1S1J2");
126 
127  // Store sum over multiplets for Monte Carlo generation.
128  for (int i = 0; i < 4; ++i) mesonRateSum[i]
129  = mesonRate[i][0] + mesonRate[i][1] + mesonRate[i][2]
130  + mesonRate[i][3] + mesonRate[i][4] + mesonRate[i][5];
131 
132  // Parameters for uubar - ddbar - ssbar meson mixing.
133  for (int spin = 0; spin < 6; ++spin) {
134  double theta;
135  if (spin == 0) theta = parm("StringFlav:thetaPS");
136  else if (spin == 1) theta = parm("StringFlav:thetaV");
137  else if (spin == 2) theta = parm("StringFlav:thetaL1S0J1");
138  else if (spin == 3) theta = parm("StringFlav:thetaL1S1J0");
139  else if (spin == 4) theta = parm("StringFlav:thetaL1S1J1");
140  else theta = parm("StringFlav:thetaL1S1J2");
141  double alpha = (spin == 0) ? 90. - (theta + 54.7) : theta + 54.7;
142  alpha *= M_PI / 180.;
143  // Fill in (flavour, spin)-dependent probability of producing
144  // the lightest or the lightest two mesons of the nonet.
145  mesonMix1[0][spin] = 0.5;
146  mesonMix2[0][spin] = 0.5 * (1. + pow2(sin(alpha)));
147  mesonMix1[1][spin] = 0.;
148  mesonMix2[1][spin] = pow2(cos(alpha));
149  // Fill in rates for multiplication.
150  mesMixRate1[0][spin] = mesonMix1[0][spin];
151  mesMixRate2[0][spin] = mesonMix2[0][spin] - mesonMix1[0][spin];
152  mesMixRate3[0][spin] = 1.0 - mesMixRate1[0][spin] - mesMixRate2[0][spin];
153  mesMixRate1[1][spin] = mesonMix1[1][spin];
154  mesMixRate2[1][spin] = mesonMix2[1][spin] - mesonMix1[1][spin];
155  mesMixRate3[1][spin] = 1.0 - mesMixRate1[1][spin] - mesMixRate2[1][spin];
156  }
157 
158  // Additional suppression of eta and etaPrime.
159  etaSup = parm("StringFlav:etaSup");
160  etaPrimeSup = parm("StringFlav:etaPrimeSup");
161 
162  // Sum of baryon octet and decuplet weights.
163  decupletSup = parm("StringFlav:decupletSup");
164  for (int i = 0; i < 6; ++i) baryonCGSum[i]
165  = baryonCGOct[i] + decupletSup * baryonCGDec[i];
166 
167  // Maximum SU(6) weight for ud0, ud1, uu1 types.
168  baryonCGMax[0] = max( baryonCGSum[0], baryonCGSum[1]);
169  baryonCGMax[1] = baryonCGMax[0];
170  baryonCGMax[2] = max( baryonCGSum[2], baryonCGSum[3]);
171  baryonCGMax[3] = baryonCGMax[2];
172  baryonCGMax[4] = max( baryonCGSum[4], baryonCGSum[5]);
173  baryonCGMax[5] = baryonCGMax[4];
174 
175  // Popcorn baryon parameters.
176  popcornRate = parm("StringFlav:popcornRate");
177  popcornSpair = parm("StringFlav:popcornSpair");
178  popcornSmeson = parm("StringFlav:popcornSmeson");
179 
180  // Suppression of leading (= first-rank) baryons.
181  suppressLeadingB = flag("StringFlav:suppressLeadingB");
182  lightLeadingBSup = parm("StringFlav:lightLeadingBSup");
183  heavyLeadingBSup = parm("StringFlav:heavyLeadingBSup");
184 
185  // Use Gaussian model but with mT2 suppression?
186  mT2suppression = flag("StringPT:mT2suppression");
187  sigmaHad = (sqrt(2.0)*parm("StringPT:sigma"));
188  widthPreStrange = parm("StringPT:widthPreStrange");
189  widthPreDiquark = parm("StringPT:widthPreDiquark");
190  useWidthPre = (widthPreStrange > 1.0) || (widthPreDiquark > 1.0);
191 
192  // Enhanded-rate prefactor for MPIs and/or nearby string pieces.
193  closePacking = flag("StringPT:closePacking");
194  exponentMPI = parm("StringPT:expMPI");
195  exponentNSP = parm("StringPT:expNSP");
196 
197  // Begin calculation of derived parameters for baryon production.
198 
199  // Enumerate distinguishable diquark types (in diquark first is popcorn q).
200  enum Diquark {ud0, ud1, uu1, us0, su0, us1, su1, ss1};
201 
202  // Maximum SU(6) weight by diquark type.
203  double barCGMax[8];
204  barCGMax[ud0] = baryonCGMax[0];
205  barCGMax[ud1] = baryonCGMax[4];
206  barCGMax[uu1] = baryonCGMax[2];
207  barCGMax[us0] = baryonCGMax[0];
208  barCGMax[su0] = baryonCGMax[0];
209  barCGMax[us1] = baryonCGMax[4];
210  barCGMax[su1] = baryonCGMax[4];
211  barCGMax[ss1] = baryonCGMax[2];
212 
213  // Diquark SU(6) survival = Sum_quark (quark tunnel weight) * SU(6).
214  double dMB[8];
215  dMB[ud0] = 2. * baryonCGSum[0] + probStoUD * baryonCGSum[1];
216  dMB[ud1] = 2. * baryonCGSum[4] + probStoUD * baryonCGSum[5];
217  dMB[uu1] = baryonCGSum[2] + (1. + probStoUD) * baryonCGSum[3];
218  dMB[us0] = (1. + probStoUD) * baryonCGSum[0] + baryonCGSum[1];
219  dMB[su0] = dMB[us0];
220  dMB[us1] = (1. + probStoUD) * baryonCGSum[4] + baryonCGSum[5];
221  dMB[su1] = dMB[us1];
222  dMB[ss1] = probStoUD * baryonCGSum[2] + 2. * baryonCGSum[3];
223  for (int i = 1; i < 8; ++i) dMB[i] = dMB[i] / dMB[0];
224 
225  // Tunneling factors for diquark production; only half a pair = sqrt.
226  double probStoUDroot = sqrt(probStoUD);
227  double probSQtoQQroot = sqrt(probSQtoQQ);
228  double probQQ1toQQ0root = sqrt(probQQ1toQQ0);
229  double qBB[8];
230  qBB[ud1] = probQQ1toQQ0root;
231  qBB[uu1] = probQQ1toQQ0root;
232  qBB[us0] = probSQtoQQroot;
233  qBB[su0] = probStoUDroot * probSQtoQQroot;
234  qBB[us1] = probQQ1toQQ0root * qBB[us0];
235  qBB[su1] = probQQ1toQQ0root * qBB[su0];
236  qBB[ss1] = probStoUDroot * pow2(probSQtoQQroot) * probQQ1toQQ0root;
237 
238  // spin * (vertex factor) * (half-tunneling factor above).
239  double qBM[8];
240  qBM[ud1] = 3. * qBB[ud1];
241  qBM[uu1] = 6. * qBB[uu1];
242  qBM[us0] = probStoUD * qBB[us0];
243  qBM[su0] = qBB[su0];
244  qBM[us1] = probStoUD * 3. * qBB[us1];
245  qBM[su1] = 3. * qBB[su1];
246  qBM[ss1] = probStoUD * 6. * qBB[ss1];
247 
248  // Combine above two into total diquark weight for q -> B Bbar.
249  for (int i = 1; i < 8; ++i) qBB[i] = qBB[i] * qBM[i];
250 
251  // Suppression from having strange popcorn meson.
252  qBM[us0] *= popcornSmeson;
253  qBM[us1] *= popcornSmeson;
254  qBM[ss1] *= popcornSmeson;
255 
256  // Suppression for a heavy quark of a diquark to fit into a baryon
257  // on the other side of popcorn meson: (0) s/u for q -> B M;
258  // (1) s/u for rank 0 diquark su -> M B; (2) ditto for s -> c/b.
259  double uNorm = 1. + qBM[ud1] + qBM[uu1] + qBM[us0] + qBM[us1];
260  scbBM[0] = (2. * (qBM[su0] + qBM[su1]) + qBM[ss1]) / uNorm;
261  scbBM[1] = scbBM[0] * popcornSpair * qBM[su0] / qBM[us0];
262  scbBM[2] = (1. + qBM[ud1]) * (2. + qBM[us0]) / uNorm;
263 
264  // Include maximum of Clebsch-Gordan coefficients.
265  for (int i = 1; i < 8; ++i) dMB[i] *= qBM[i];
266  for (int i = 1; i < 8; ++i) qBM[i] *= barCGMax[i] / barCGMax[0];
267  for (int i = 1; i < 8; ++i) qBB[i] *= barCGMax[i] / barCGMax[0];
268 
269  // Popcorn fraction for normal diquark production.
270  double qNorm = uNorm * popcornRate / 3.;
271  double sNorm = scbBM[0] * popcornSpair;
272  popFrac = qNorm * (1. + qBM[ud1] + qBM[uu1] + qBM[us0] + qBM[us1]
273  + sNorm * (qBM[su0] + qBM[su1] + 0.5 * qBM[ss1])) / (1. + qBB[ud1]
274  + qBB[uu1] + 2. * (qBB[us0] + qBB[us1]) + 0.5 * qBB[ss1]);
275 
276  // Popcorn fraction for rank 0 diquarks, depending on number of s quarks.
277  popS[0] = qNorm * qBM[ud1] / qBB[ud1];
278  popS[1] = qNorm * 0.5 * (qBM[us1] / qBB[us1]
279  + sNorm * qBM[su1] / qBB[su1]);
280  popS[2] = qNorm * sNorm * qBM[ss1] / qBB[ss1];
281 
282  // Recombine diquark weights to flavour and spin ratios. Second index:
283  // 0 = s/u popcorn quark ratio.
284  // 1, 2 = s/u ratio for vertex quark if popcorn quark is u/d or s.
285  // 3 = q/q' vertex quark ratio if popcorn quark is light and = q.
286  // 4, 5, 6 = (spin 1)/(spin 0) ratio for su, us and ud.
287 
288  // Case 0: q -> B B.
289  dWT[0][0] = (2. * (qBB[su0] + qBB[su1]) + qBB[ss1])
290  / (1. + qBB[ud1] + qBB[uu1] + qBB[us0] + qBB[us1]);
291  dWT[0][1] = 2. * (qBB[us0] + qBB[us1]) / (1. + qBB[ud1] + qBB[uu1]);
292  dWT[0][2] = qBB[ss1] / (qBB[su0] + qBB[su1]);
293  dWT[0][3] = qBB[uu1] / (1. + qBB[ud1] + qBB[uu1]);
294  dWT[0][4] = qBB[su1] / qBB[su0];
295  dWT[0][5] = qBB[us1] / qBB[us0];
296  dWT[0][6] = qBB[ud1];
297 
298  // Case 1: q -> B M B.
299  dWT[1][0] = (2. * (qBM[su0] + qBM[su1]) + qBM[ss1])
300  / (1. + qBM[ud1] + qBM[uu1] + qBM[us0] + qBM[us1]);
301  dWT[1][1] = 2. * (qBM[us0] + qBM[us1]) / (1. + qBM[ud1] + qBM[uu1]);
302  dWT[1][2] = qBM[ss1] / (qBM[su0] + qBM[su1]);
303  dWT[1][3] = qBM[uu1] / (1. + qBM[ud1] + qBM[uu1]);
304  dWT[1][4] = qBM[su1] / qBM[su0];
305  dWT[1][5] = qBM[us1] / qBM[us0];
306  dWT[1][6] = qBM[ud1];
307 
308  // Case 2: qq -> M B; diquark inside chain.
309  dWT[2][0] = (2. * (dMB[su0] + dMB[su1]) + dMB[ss1])
310  / (1. + dMB[ud1] + dMB[uu1] + dMB[us0] + dMB[us1]);
311  dWT[2][1] = 2. * (dMB[us0] + dMB[us1]) / (1. + dMB[ud1] + dMB[uu1]);
312  dWT[2][2] = dMB[ss1] / (dMB[su0] + dMB[su1]);
313  dWT[2][3] = dMB[uu1] / (1. + dMB[ud1] + dMB[uu1]);
314  dWT[2][4] = dMB[su1] / dMB[su0];
315  dWT[2][5] = dMB[us1] / dMB[us0];
316  dWT[2][6] = dMB[ud1];
317 
318  // Use thermal model?
319  thermalModel = flag("StringPT:thermalModel");
320  if (thermalModel || mT2suppression) {
321 
322  // Temperature parameters for thermal model.
323  temperature = parm("StringPT:temperature");
324  tempPreFactor = parm("StringPT:tempPreFactor");
325 
326  // Hadron multiplets in thermal model.
327  mesonNonetL1 = flag("StringFlav:mesonNonetL1");
328  nNewQuark = mode("StringFlav:nQuark");
329 
330  // Fill list of possible hadrons that are allowed to be produced.
331  // Also include a list of "emergency" hadrons that are needed to get
332  // rid of all possible endpoint (di)quarks.
333  vector<int> hadIDsProd, hadIDsHvyC, hadIDsHvyB;
334 
335  // Baryon octet and decuplet.
336  int baryonLight[18] = { 2112, 2212, 3112, 3122, 3212, 3222, 3312, 3322,
337  1114, 2114, 2214, 2224, 3114, 3214, 3224, 3314,
338  3324, 3334 };
339  int baryonHvyC[22] = { 4112, 4122, 4132, 4212, 4222, 4232, 4312, 4322,
340  4332, 4412, 4422, 4432,
341  4114, 4214, 4224, 4314, 4324, 4334, 4414, 4424,
342  4434, 4444 };
343  int baryonHvyB[35] = { 5112, 5122, 5132, 5142, 5212, 5222, 5232, 5242,
344  5312, 5322, 5332, 5342, 5412, 5422, 5432, 5442,
345  5512, 5522, 5532, 5542,
346  5114, 5214, 5224, 5314, 5324, 5334, 5414, 5424,
347  5434, 5444, 5514, 5524, 5534, 5544, 5554 };
348  for (int i = 0; i < 18; i++) hadIDsProd.push_back( baryonLight[i] );
349  // Check how many heavy baryons to include.
350  if (nNewQuark > 4) {
351  for (int i = 0; i < 35; i++) hadIDsProd.push_back( baryonHvyB[i] );
352  } else {
353  // Only include lightest combinations.
354  int bBar[9] = { 5112, 5122, 5132, 5212, 5222, 5232, 5312, 5322, 5332 };
355  for (int i = 0; i < 9; i++) {
356  hadIDsHvyB.push_back( bBar[i] );
357  hadIDsHvyB.push_back( -bBar[i] );
358  }
359  }
360  if (nNewQuark > 3) {
361  for (int i = 0; i < 22; i++) hadIDsProd.push_back( baryonHvyC[i] );
362  } else {
363  // Only include lightest combinations.
364  int cBar[9] = { 4112, 4122, 4132, 4212, 4222, 4232, 4312, 4322, 4332 };
365  for (int i = 0; i < 9; i++) {
366  hadIDsHvyC.push_back( cBar[i] );
367  hadIDsHvyC.push_back( -cBar[i] );
368  }
369  }
370  // Antibaryons.
371  int sizeNow = int(hadIDsProd.size());
372  for (int i = 0; i < sizeNow; i++) hadIDsProd.push_back( -hadIDsProd[i] );
373 
374  // Mesons nonets. Take pseudoscalar PDG codes as basis.
375  int mesonPSLight[9] = { 311, 321, 211, -311, -321, -211, 111, 221, 331 };
376  int mesonPSHvyC[7] = { 411, 421, 431, -411, -421, -431, 441 };
377  int mesonPSHvyB[9] = { 511, 521, 531, 541, -511, -521, -531, -541, 551 };
378  vector<int> mesonPS;
379  for (int i = 0; i < 9; i++) mesonPS.push_back( mesonPSLight[i] );
380  // Check how many heavy mesons to include. If not included in ordinary
381  // production, fill minimal list with "emergency" hadrons
382  if (nNewQuark > 4) {
383  for (int i = 0; i < 9; i++) mesonPS.push_back( mesonPSHvyB[i] );
384  } else {
385  // Include all possible combinations, only pseudoscalar as they
386  // are the lightest ones.
387  int bMes[10] = { 511, 521, 531, 541, -511, -521, -531, -541, 551 };
388  for (int i = 0; i < 10; i++) hadIDsHvyB.push_back( bMes[i] );
389  }
390  if (nNewQuark > 3) {
391  for (int i = 0; i < 7; i++) mesonPS.push_back( mesonPSHvyC[i] );
392  } else {
393  // Include all possible combinations, only pseudoscalar as they
394  // are the lightest ones.
395  int cMes[8] = { 411, 421, 431, -411, -421, -431, 441 };
396  for (int i = 0; i < 8; i++) hadIDsHvyC.push_back( cMes[i] );
397  }
398  int nMeson = int(mesonPS.size());
399  // Pseudoscalar nonet J=0, S=0, L=0.
400  for (int i = 0; i < nMeson; i++)
401  hadIDsProd.push_back( mesonPS[i] );
402  // Vector nonet J=1, S=1, L=0.
403  for (int i = 0; i < nMeson; i++)
404  hadIDsProd.push_back( mesonPS[i] + (mesonPS[i] > 0 ? 2 : -2) );
405  // Include L=1 nonets?
406  if (mesonNonetL1) {
407  // Pseudovector nonet J=1, S=0, L=1.
408  for (int i = 0; i < nMeson; i++)
409  hadIDsProd.push_back( mesonPS[i] + (mesonPS[i] > 0 ? 10002 : -10002) );
410  // Scalar nonet J=0, S=1, L=1.
411  for (int i = 0; i < nMeson; i++)
412  hadIDsProd.push_back( mesonPS[i] + (mesonPS[i] > 0 ? 10000 : -10000) );
413  // Pseudovector nonet J=1, S=1, L=1.
414  for (int i = 0; i < nMeson; i++)
415  hadIDsProd.push_back( mesonPS[i] + (mesonPS[i] > 0 ? 20002 : -20002) );
416  // Tensor nonet J=2, S=1, L=1.
417  for (int i = 0; i < nMeson; i++)
418  hadIDsProd.push_back( mesonPS[i] + (mesonPS[i] > 0 ? 4 : -4) );
419  }
420 
421  // Fill list of all hadrons ids (ordinary and "emergency").
422  vector<int> hadIDsAll;
423  for (int i = 0; i < int(hadIDsProd.size()); i++)
424  hadIDsAll.push_back( hadIDsProd[i] );
425  for (int i = 0; i < int(hadIDsHvyC.size()); i++)
426  hadIDsAll.push_back( hadIDsHvyC[i] );
427  for (int i = 0; i < int(hadIDsHvyB.size()); i++)
428  hadIDsAll.push_back( hadIDsHvyB[i] );
429 
430  // Fill map with IDs of hadron constituents for all hadrons.
431  for (int i = 0; i < int(hadIDsAll.size()); i++) {
432  int id = hadIDsAll[i];
433  int idAbs = abs(id);
434  vector< pair<int,int> > quarkCombis;
435  // Baryon can be split into q + qq in several different ways.
436  if (particleDataPtr->isBaryon(id)) {
437  bool isOctet = ( (idAbs % 10) == 2 );
438  int q3 = (idAbs/10) % 10;
439  int q2 = (idAbs/100) % 10;
440  int q1 = (idAbs/1000) % 10;
441  bool threeFlav = q1 != q2 && q1 != q3 && q2 != q3;
442  // Baryon octet J=1/2.
443  if (isOctet) {
444  if (threeFlav) {
445  // Add (q2+q3)_0/1 + q1.
446  // if (q2 < q3) (q2+q3)_0 and if (q2 > q3) (q2+q3)_1.
447  int j = (q2 < q3) ? 1 : 3;
448  int qn[2] = { min( q3, q2), max( q3, q2) };
449  addQuarkDiquark(quarkCombis, q1,
450  1000 * qn[1] + 100 * qn[0] + j, id);
451  // Add other combinations. Can be both, J=0 or J=1.
452  for (j = 1; j < 4; j += 2) {
453  // (q1+q3)j + q2
454  addQuarkDiquark(quarkCombis, q2, 1000 * q1 + 100 * q3 + j, id);
455  // (q1+q2)j + q3
456  addQuarkDiquark(quarkCombis, q3, 1000 * q1 + 100 * q2 + j, id);
457  }
458  } else {
459  // Quarks with the same flavour form J=1,
460  // all other combinations can be both, J=0 or J=1.
461  for (int j = 1; j < 4; j += 2) {
462  // (q1+q2)1 + q3
463  if ( j == 3 || q1 != q2 )
464  addQuarkDiquark(quarkCombis, q3, 1000 * q1 + 100 * q2 + j, id);
465  // (q1+q3)1 + q2
466  if ( j == 3 || q1 != q3 )
467  addQuarkDiquark(quarkCombis, q2, 1000 * q1 + 100 * q3 + j, id);
468  // (q2+q3)1 + q1
469  if ( j == 3 || q2 != q3 )
470  addQuarkDiquark(quarkCombis, q1, 1000 * q2 + 100 * q3 + j, id);
471  }
472  }
473  }
474  // Baryon decuplet J=3/2.
475  else {
476  // All quark pairs form diquarks with J=1.
477  // (q1+q2)1 + q3
478  addQuarkDiquark(quarkCombis, q3, 1000 * q1 + 100 * q2 + 3, id);
479  // (q1+q3)1 + q2
480  addQuarkDiquark(quarkCombis, q2, 1000 * q1 + 100 * q3 + 3, id);
481  // (q2+q3)1 + q1
482  addQuarkDiquark(quarkCombis, q1, 1000 * q2 + 100 * q3 + 3, id);
483  }
484  // Mesons usually have a trivial subdivision into quark + antiquark.
485  // Mixing of diagonal mesons is taken into account later.
486  } else {
487  int q1 = (idAbs/100) % 10;
488  bool uptype1 = (q1 % 2 == 0);
489  int q2 = (idAbs/10) % 10;
490  bool uptype2 = (q2 % 2 == 0);
491  int quark = q1;
492  int antiQuark = q2;
493  // id > 0: downtype+uptype: up = quark, down = antiquark (default)
494  // id > 0: same type -> larger id decides
495  if ( uptype2 && !uptype1 ) swap( quark, antiQuark);
496  if ( (q1 > q2 && !uptype1 && !uptype2)
497  || (q2 > q1 && uptype2 && uptype1) ) swap( quark, antiQuark);
498  if (id < 0) swap( quark, antiQuark);
499  quarkCombis.push_back( make_pair( quark, -antiQuark) );
500  }
501  hadronConstIDs[id] = quarkCombis;
502  }
503 
504  // Copy into smaller versions (one for ordinary production, two for
505  // "emergency")
506  map< int, vector< pair<int,int> > > hadConstIDsC, hadConstIDsB,
507  hadConstIDsProd;
508  for (int i=0; i<int(hadIDsAll.size()); i++) {
509  int id = hadIDsAll[i];
510  if (find(hadIDsProd.begin(), hadIDsProd.end(), id) != hadIDsProd.end())
511  hadConstIDsProd[id] = hadronConstIDs[id];
512  if (find(hadIDsHvyC.begin(), hadIDsHvyC.end(), id) != hadIDsHvyC.end())
513  hadConstIDsC[id] = hadronConstIDs[id];
514  if (find(hadIDsHvyB.begin(), hadIDsHvyB.end(), id) != hadIDsHvyB.end())
515  hadConstIDsB[id] = hadronConstIDs[id];
516  }
517  map< int, map< int, vector< pair<int,int> > > > hadConstIDsHvy;
518  hadConstIDsHvy[4] = hadConstIDsC;
519  hadConstIDsHvy[5] = hadConstIDsB;
520 
521  // List with all possible initial (di)quarks we could get.
522  int inIDs[26] = { 1, 2, 3, 4, 5, 1103, 2203, 3303, 2101, 2103, 3101,
523  3103, 3201, 3203, 4101, 4103, 4201, 4203, 4301,
524  4303, 5101, 5103, 5201, 5203, 5301, 5303 };
525  int inIDsHvyC[2] = { 4403, -4403 };
526  int inIDsHvyB[6] = { 5503, -5503, 5401, -5401, 5403, -5403 };
527  vector<int> incomingIDs;
528  for (int i = 0; i < 26; i++) {
529  incomingIDs.push_back( inIDs[i]);
530  incomingIDs.push_back(-inIDs[i]);
531  }
532  // If we include heavy quark hadrons we include the following diquarks in
533  // addition.
534  if (nNewQuark > 3) {
535  for (int i = 0; i < 2; i++) incomingIDs.push_back(inIDsHvyC[i]);
536  if (nNewQuark > 4) {
537  for (int i = 0; i < 6; i++) incomingIDs.push_back( inIDsHvyB[i]);
538  }
539  }
540  int nIncome = int(incomingIDs.size());
541 
542  // Loop over list with all possible initial (di)quarks.
543  // Fill map possibleHadrons with
544  // key = initial (di)quark id, value = list of possible hadron ids
545  // + nr in hadronConstIDs.
546  for (int iIDin = 0; iIDin < nIncome; iIDin++) {
547  int idIn = incomingIDs[iIDin];
548  int idInAbs = abs(idIn);
549  map< int, vector< pair<int,int> > > hadConstIDsNow = hadConstIDsProd;
550  // For heavy quarks add "emergency" list, if needed.
551  for (int iHvy = nNewQuark+1; iHvy <= 5; iHvy++) {
552  if (particleDataPtr->nQuarksInCode(idInAbs, iHvy) > 0)
553  for (map< int, vector< pair<int,int> > >::iterator
554  it = hadConstIDsHvy[iHvy].begin();
555  it != hadConstIDsHvy[iHvy].end(); ++it)
556  hadConstIDsNow[it->first] = it->second;
557  }
558  // Fill list: first parameter of pair is hadron ID, second is nr of
559  // hadron constituents in the list.
560  vector< pair<int,int> > possibleHadronIDs;
561  // Loop through list with hadrons and their (di)quark content,
562  // check if possible to produce given the choice of initial (di)quark.
563  for (map< int, vector< pair<int,int> > >::iterator
564  it = hadConstIDsNow.begin(); it != hadConstIDsNow.end(); ++it) {
565  vector< pair<int,int> > constituentIDs = it->second;
566  int nConst = int(constituentIDs.size());
567  int hadronID = it->first;
568  // Loop over constituent IDs.
569  for (int iConst = 0; iConst < nConst; iConst++) {
570  int ID1 = constituentIDs[iConst].first;
571  int ID2 = constituentIDs[iConst].second;
572  if ( (ID1 == idIn) || (ID2 == idIn) ) {
573  possibleHadronIDs.push_back( make_pair(hadronID,iConst) );
574  // To include uubar-ddbar-ssbar mixing include all diagonal mesons.
575  if ( (idInAbs < 4) && (ID1 == -ID2) ) {
576  if (idInAbs == 1) {
577  possibleHadronIDs.push_back( make_pair(hadronID+110,iConst) );
578  possibleHadronIDs.push_back( make_pair(hadronID+220,iConst) );
579  } else if (idInAbs == 2) {
580  possibleHadronIDs.push_back( make_pair(hadronID-110,iConst) );
581  possibleHadronIDs.push_back( make_pair(hadronID+110,iConst) );
582  } else if (idInAbs == 3) {
583  possibleHadronIDs.push_back( make_pair(hadronID-220,iConst) );
584  possibleHadronIDs.push_back( make_pair(hadronID-110,iConst) );
585  }
586  }
587  }
588  }
589  }
590  if (int(possibleHadronIDs.size()) < 1)
591  infoPtr->errorMsg("Error in StringFlav::init: no possible "
592  "hadrons found");
593  possibleHadrons[idIn] = possibleHadronIDs;
594  }
595 
596  // Calculate baryon octet and decuplet weighting factors
597  // based on Clebsch-Gordan coefficients and spin counting.
598  // Parameters: qDi1 qDi2 q3 spin.
599  // Zero for flavour=0 and same flavour diquarks with J=0.
600  for (int q1 = 0; q1 < 6; q1++) {
601  for (int q2 = 0; q2 < 6; q2++) {
602  baryonOctWeight[q1][q1][q2][0] = 0.0; // qq0 + r
603  baryonDecWeight[q1][q1][q2][0] = 0.0; // qq0 + r
604  for (int spin = 0; spin < 1; spin++) {
605  baryonOctWeight[ 0][q1][q2][spin] = 0.0;
606  baryonOctWeight[q1][ 0][q2][spin] = 0.0;
607  baryonOctWeight[q1][q2][ 0][spin] = 0.0;
608  baryonDecWeight[ 0][q1][q2][spin] = 0.0;
609  baryonDecWeight[q1][ 0][q2][spin] = 0.0;
610  baryonDecWeight[q1][q2][ 0][spin] = 0.0;
611  }
612  }
613  }
614  // Clebsch-Gordon for the rest.
615  for (int q1 = 1; q1 < 6; q1++) {
616  baryonOctWeight[q1][q1][q1][1] = 0.0; // qq1 + q
617  baryonDecWeight[q1][q1][q1][1] = 1.0;
618  for (int q2 = 1; q2 < 6; q2++) if (q1!=q2) {
619  baryonOctWeight[q1][q1][q2][1] = 0.1667; // qq1 + r
620  baryonDecWeight[q1][q1][q2][1] = 0.3333;
621  baryonOctWeight[q1][q2][q1][0] = 0.75; // qr0 + q
622  baryonDecWeight[q1][q2][q1][0] = 0.0;
623  baryonOctWeight[q2][q1][q1][0] = 0.75; // rq0 + q
624  baryonDecWeight[q2][q1][q1][0] = 0.0;
625  baryonOctWeight[q1][q2][q1][1] = 0.0833; // qr1 + q
626  baryonDecWeight[q1][q2][q1][1] = 0.6667;
627  baryonOctWeight[q2][q1][q1][1] = 0.0833; // rq1 + q
628  baryonDecWeight[q2][q1][q1][1] = 0.6667;
629  for (int q3 = 0; q3 < 6; q3++) if ((q1 != q3) && (q2 != q3)) {
630  baryonOctWeight[q1][q2][q3][0] = 0.5; // qr0 + s
631  baryonDecWeight[q1][q2][q3][0] = 0.0;
632  baryonOctWeight[q1][q2][q3][1] = 0.1667; // qr1 + s
633  baryonDecWeight[q1][q2][q3][1] = 0.3333;
634  }
635  }
636  }
637  // Spin 1 diquarks get extra factor of 3. And all factors
638  // get relative baryon-to-meson ratio.
639  double BtoMratio = parm("StringFlav:BtoMratio");
640  for (int q1 = 0; q1 < 6; q1++) {
641  for (int q2 = 0; q2 < 6; q2++) {
642  for (int q3 = 0; q3 < 6; q3++) {
643  for (int spin = 0; spin < 2; spin++) {
644  baryonOctWeight[q1][q2][q3][spin] *= BtoMratio;
645  baryonDecWeight[q1][q2][q3][spin] *= BtoMratio;
646  if (spin == 1) {
647  baryonOctWeight[q1][q2][q3][1] *= 3.0;
648  baryonDecWeight[q1][q2][q3][1] *= 3.0;
649  }
650  }
651  }
652  }
653  }
654 
655  // Go through the list of possible hadrons and calculate the prefactor
656  // that will multiply the rate.
657  double strSup = parm("StringFlav:StrangeSuppression");
658  for (int iIDin = 0; iIDin < nIncome; iIDin++) {
659  int idIn = incomingIDs[iIDin];
660  int idInAbs = abs(idIn);
661  vector< pair<int,int> > possibleHadronsNow = possibleHadrons[idIn];
662  vector<double> prefactors;
663  for (int iHad = 0; iHad < int(possibleHadronsNow.size()); iHad++) {
664  double prefacNow = 1.0;
665  // Get hadron and constituents.
666  int hadronID = possibleHadronsNow[iHad].first;
667  int hadronIDabs = abs(hadronID);
668  int iConst = possibleHadronsNow[iHad].second;
669  int ID1 = hadronConstIDs[hadronID][iConst].first;
670  int ID2 = hadronConstIDs[hadronID][iConst].second;
671  // Extra suppression factor for s/c/b quarks.
672  double nHeavy = 0.0;
673  for (int i = 3; i <= 5; i++) {
674  nHeavy += particleDataPtr->nQuarksInCode( ID1, i);
675  nHeavy += particleDataPtr->nQuarksInCode( ID2, i);
676  }
677  prefacNow *= pow(strSup, nHeavy);
678  if (particleDataPtr->isMeson(hadronID)) {
679  // Extra factor according to last digit for spin counting.
680  prefacNow *= (abs(hadronID) % 10);
681  // Include correct uubar-ddbar-ssbar mixing factor;
682  if ( (idInAbs < 4) && (ID1 == -ID2) ) {
683  int flav = ( (idInAbs < 3) ? 0 : 1 );
684  // Get spin used as counter for the different multiplets
685  int spin = getMesonSpinCounter(hadronID);
686  double mesonMix[3] = { mesMixRate1[flav][spin],
687  mesMixRate2[flav][spin],
688  mesMixRate3[flav][spin] };
689  prefacNow *= mesonMix[abs(ID1)-1];
690  }
691  } else {
692  // Check if baryon is octet or decuplet.
693  bool isOct = ((hadronIDabs % 10) == 2);
694  // Make sure ID2 is diquark.
695  if (abs(ID2) < abs(ID1)) swap(ID1,ID2);
696  // Extract quark flavours and spin from diquark.
697  int Q1 = ( (abs(ID2)/1000) % 10 );
698  int Q2 = ( (abs(ID2)/100) % 10 );
699  if (Q1 > 5 || Q2 > 5) {
700  infoPtr->errorMsg("Error in StringFlav::init: invalid quark "
701  "content flavours for diquark");
702  continue;
703  }
704  int diqSpin = ( ((abs(ID2) % 10) == 1) ? 0 : 1 );
705  // Single quark.
706  int Q3 = abs(ID1);
707  // Find Clebsch-Gordan: q1 in DQ | q2 in DQ | q3 | S of DQ
708  if (isOct) prefacNow *= baryonOctWeight[Q1][Q2][Q3][diqSpin];
709  else prefacNow *= baryonDecWeight[Q1][Q2][Q3][diqSpin];
710  // Special cases for Lamda (312) and Sigma (321) or the like.
711  if ( isOct && (Q1!=Q2) && (Q1!=Q3) && (Q2!=Q3) ) {
712  // Extract the two lightest quarks from hadron.
713  int Qhad1 = ( (hadronIDabs/10) % 10 );
714  int Qhad2 = ( (hadronIDabs/100) % 10 );
715  int QhadMin = min(Qhad1,Qhad2);
716  int QhadMax = max(Qhad1,Qhad2);
717  // Extract the two quarks from the diquark.
718  int QdiqMin = min(Q1,Q2);
719  int QdiqMax = max(Q1,Q2);
720  // Don't do anything if (12) or (21) is diquark.
721  if ( !((QdiqMin == QhadMin) && (QdiqMax == QhadMax)) ) {
722  // Sigma (321)
723  if (Qhad2 > Qhad1) prefacNow *= ( (diqSpin == 0) ? 0.75 : 0.25 );
724  // Lamda (312)
725  else prefacNow *= ( (diqSpin == 0) ? 0.25 : 0.27 );
726  }
727  }
728  }
729  // Save prefactor.
730  prefactors.push_back(prefacNow);
731  }
732  possibleRatePrefacs[idIn] = prefactors;
733  }
734 
735  // Now the same again for joining the last two (di)quarks into hadron.
736  for (int iIDin1 = 0; iIDin1 < nIncome; iIDin1++) {
737  int idIn1 = incomingIDs[iIDin1];
738  int idIn1Abs = abs(idIn1);
739  // Loop over possible partners, start with next quark.
740  for (int iIDin2 = iIDin1+1; iIDin2 < nIncome; iIDin2++) {
741  int idIn2 = incomingIDs[iIDin2];
742  int idIn2Abs = abs(idIn2);
743  int idInNow[2] = { min(idIn1,idIn2), max(idIn1,idIn2) };
744  pair<int,int> inPair = pair<int,int>(idInNow[0], idInNow[1]);
745  // Skip all combinations with two diquarks.
746  if ( (idIn1Abs > 1000) && (idIn2Abs > 1000) ) continue;
747  // Skip all combinations with two quarks or two antiquarks.
748  if ( ( ((idIn1 > 0) && (idIn2 > 0)) || ((idIn1 < 0) && (idIn2 < 0)) )
749  && (idIn1Abs < 10) && (idIn2Abs < 10) ) continue;
750  // Skip all combinations with quark-antidiquark and
751  // antiquark-diquark. (1 = diquark, 2 = quark not possible).
752  if ( ((idIn2 > 1000) && (idIn1Abs < 10) && (idIn1 < 0)) ||
753  ((idIn2 < -1000) && (idIn1Abs < 10) && (idIn1 > 0)) ) continue;
754  // If we are not including heavy quarks skip combinations
755  // of heavy quark - diquark with heavy quark.
756  if ((idIn1Abs < 10) && (idIn2Abs > 1000)) {
757  vector< pair<int,int> > hvyCombs;
758  if (nNewQuark < 5) {
759  hvyCombs.push_back(make_pair(4,4));
760  if (nNewQuark < 4) {
761  hvyCombs.push_back(make_pair(5,4));
762  hvyCombs.push_back(make_pair(4,5));
763  hvyCombs.push_back(make_pair(5,5));
764  }
765  }
766  bool skip = false;
767  for (int iComb = 0; iComb < int(hvyCombs.size()); iComb++) {
768  int idNow[2] = { hvyCombs[iComb].first, hvyCombs[iComb].second };
769  if ( (particleDataPtr->nQuarksInCode(idIn2Abs,idNow[0]) > 0) &&
770  (idIn1Abs == idNow[1]) ) skip = true;
771  }
772  if (skip) continue;
773  }
774  // Now decide which list of possible hadrons to use.
775  // As we might have to use the special list for heavy quarks we
776  // use the maximum of the absolute ids in case of two quarks and
777  // check the maximum flavour in case of quark - diquark pair.
778  int idUse;
779  if ( (idIn1Abs < 10) && (idIn2Abs < 10) ) { // quark - quark
780  idUse = ( (idIn1Abs > idIn2Abs) ? idIn1 : idIn2 );
781  } else { // quark - diquark
782  // Check if diquark contains a heavier flavour then the quark.
783  bool useDiquark = false;
784  for (int plus = 1; plus < 5; plus++)
785  if (particleDataPtr->nQuarksInCode(idIn2Abs, idIn1Abs + plus) > 0)
786  useDiquark = true;
787  idUse = ( useDiquark ? idIn2 : idIn1 );
788  }
789  vector<double> possibleRatePrefacsNow = possibleRatePrefacs[idUse];
790  vector< pair<int,int> > possibleHadronsNow = possibleHadrons[idUse];
791  // New list to fill.
792  vector< pair<int,int> > possibleHadronsNew;
793  vector<double> possibleRatePrefacsNew;
794  // Now loop over possible hadrons and check if other (di)quark
795  // in constituents matches idIn2.
796  for (int iHad = 0; iHad < int(possibleHadronsNow.size()); iHad++) {
797  // Get constituents.
798  int hadronID = possibleHadronsNow[iHad].first;
799  int iConst = possibleHadronsNow[iHad].second;
800  int ID1 = hadronConstIDs[hadronID][iConst].first;
801  int ID2 = hadronConstIDs[hadronID][iConst].second;
802  if ( ((ID1 == idIn1) && (ID2 == idIn2)) ||
803  ((ID1 == idIn2) && (ID2 == idIn1)) ) {
804  // Can take this combination.
805  possibleHadronsNew.push_back(possibleHadronsNow[iHad]);
806  possibleRatePrefacsNew.push_back(possibleRatePrefacsNow[iHad]);
807  }
808  }
809  if (int(possibleHadronsNew.size()) < 1)
810  infoPtr->errorMsg("Error in StringFlav::init: no possible "
811  "hadrons found for last two");
812  // Save.
813  possibleRatePrefacsLast[inPair] = possibleRatePrefacsNew;
814  possibleHadronsLast[inPair] = possibleHadronsNew;
815  }
816  }
817  }
818 
819  // Initialize winning parameters.
820  hadronIDwin = 0;
821  idNewWin = 0;
822  hadronMassWin = -1.0;
823 
824 }
825 
826 //--------------------------------------------------------------------------
827 
828 // Pick a new flavour (including diquarks) given an incoming one for
829 // Gaussian pTq^2 distribution.
830 
831 FlavContainer StringFlav::pickGauss(FlavContainer& flavOld, bool allowPop) {
832 
833  // Initial values for new flavour.
834  FlavContainer flavNew;
835  flavNew.rank = flavOld.rank + 1;
836 
837  // For original diquark assign popcorn quark and whether popcorn meson.
838  int idOld = abs(flavOld.id);
839  if (flavOld.rank == 0 && idOld > 1000 && allowPop) assignPopQ(flavOld);
840 
841  // Diquark exists, to be forced into baryon now.
842  bool doOldBaryon = (idOld > 1000 && flavOld.nPop == 0);
843  // Diquark exists, but do meson now.
844  bool doPopcornMeson = flavOld.nPop > 0;
845  // Newly created diquark gives baryon now, antibaryon later.
846  bool doNewBaryon = false;
847 
848  // Choose whether to generate a new meson or a new baryon.
849  if (!doOldBaryon && !doPopcornMeson && probQandQQ * rndmPtr->flat() > 1.) {
850  doNewBaryon = true;
851  if ((1. + popFrac) * rndmPtr->flat() > 1.) flavNew.nPop = 1;
852  }
853 
854  // Optional suppression of first-rank baryon.
855  if (flavOld.rank == 0 && doNewBaryon && suppressLeadingB) {
856  double leadingBSup = (idOld < 4) ? lightLeadingBSup : heavyLeadingBSup;
857  if (rndmPtr->flat() > leadingBSup) {
858  doNewBaryon = false;
859  flavNew.nPop = 0;
860  }
861  }
862 
863  // Single quark for new meson or for baryon where diquark already exists.
864  if (!doPopcornMeson && !doNewBaryon) {
865  flavNew.id = pickLightQ();
866  if ( (flavOld.id > 0 && flavOld.id < 9) || flavOld.id < -1000 )
867  flavNew.id = -flavNew.id;
868 
869  // Done for simple-quark case.
870  return flavNew;
871  }
872 
873  // Case: 0 = q -> B B, 1 = q -> B M B, 2 = qq -> M B.
874  int iCase = flavNew.nPop;
875  if (flavOld.nPop == 1) iCase = 2;
876 
877  // Flavour of popcorn quark (= q shared between B and Bbar).
878  if (doNewBaryon) {
879  double sPopWT = dWT[iCase][0];
880  if (iCase == 1) sPopWT *= scbBM[0] * popcornSpair;
881  double rndmFlav = (2. + sPopWT) * rndmPtr->flat();
882  flavNew.idPop = 1;
883  if (rndmFlav > 1.) flavNew.idPop = 2;
884  if (rndmFlav > 2.) flavNew.idPop = 3;
885  } else flavNew.idPop = flavOld.idPop;
886 
887  // Flavour of vertex quark.
888  double sVtxWT = dWT[iCase][1];
889  if (flavNew.idPop >= 3) sVtxWT = dWT[iCase][2];
890  if (flavNew.idPop > 3) sVtxWT *= 0.5 * (1. + 1./dWT[iCase][4]);
891  double rndmFlav = (2. + sVtxWT) * rndmPtr->flat();
892  flavNew.idVtx = 1;
893  if (rndmFlav > 1.) flavNew.idVtx = 2;
894  if (rndmFlav > 2.) flavNew.idVtx = 3;
895 
896  // Special case for light flavours, possibly identical.
897  if (flavNew.idPop < 3 && flavNew.idVtx < 3) {
898  flavNew.idVtx = flavNew.idPop;
899  if (rndmPtr->flat() > dWT[iCase][3]) flavNew.idVtx = 3 - flavNew.idPop;
900  }
901 
902  // Pick 2 * spin + 1.
903  int spin = 3;
904  if (flavNew.idVtx != flavNew.idPop) {
905  double spinWT = dWT[iCase][6];
906  if (flavNew.idVtx == 3) spinWT = dWT[iCase][5];
907  if (flavNew.idPop >= 3) spinWT = dWT[iCase][4];
908  if ((1. + spinWT) * rndmPtr->flat() < 1.) spin = 1;
909  }
910 
911  // Form outgoing diquark. Done.
912  flavNew.id = 1000 * max(flavNew.idVtx, flavNew.idPop)
913  + 100 * min(flavNew.idVtx, flavNew.idPop) + spin;
914  if ( (flavOld.id < 0 && flavOld.id > -9) || flavOld.id > 1000 )
915  flavNew.id = -flavNew.id;
916  return flavNew;
917 
918 }
919 
920 //--------------------------------------------------------------------------
921 
922 // Pick a hadron, based on generated pT value and initial (di)quark.
923 // Check all possible hadrons and calculate their relative suppression
924 // based on exp(-mThadron/T), possibly multiplied by spin counting, meson
925 // mixing or baryon weighting factors.
926 // First return value is hadron ID, second new (di)quark ID.
927 
928 FlavContainer StringFlav::pickThermal(FlavContainer& flavOld,
929  double pT, double nNSP) {
930 
931  // Initial values for new flavour.
932  FlavContainer flavNew;
933  flavNew.rank = flavOld.rank + 1;
934 
935  int idIn = flavOld.id;
936  int idInAbs = abs(idIn);
937  double temprNow = temperature;
938  // Temperature increase to work against asymmetry. Apply for
939  // s/c/b and diquarks.
940  if (idInAbs > 2) temprNow *= tempPreFactor;
941  // Enhanded-rate prefactor for MPIs and/or nearby string pieces.
942  if (closePacking) {
943  temprNow *= pow(max(1.0,double(infoPtr->nMPI())), exponentMPI);
944  temprNow *= pow(max(1.0,nNSP), exponentNSP);
945  }
946  // Get Gaussian width in case of mT2 suppression.
947  double sigmaNow = sigmaHad;
948 
949  // Prefactor for strange quarks and diquarks.
950  if (useWidthPre) {
951  if (abs(idIn) > 10) sigmaNow *= widthPreDiquark;
952  sigmaNow *= pow(widthPreStrange,
953  particleDataPtr->nQuarksInCode(idIn,3) );
954  }
955 
956  // Enhanded-rate prefactor for MPIs and/or nearby string pieces.
957  if (closePacking) {
958  sigmaNow *= pow(max(1.0,double(infoPtr->nMPI())), exponentMPI);
959  sigmaNow *= pow(max(1.0,nNSP), exponentNSP);
960  }
961 
962  // Get the list of allowed hadrons and constituents for that
963  // initial (di)quark. First parameter of pair is hadron ID, second
964  // is nr of hadron constituents in the list.
965  vector<double> possibleRatePrefacsNow = possibleRatePrefacs[idIn];
966  vector< pair<int,int> > possibleHadronsNow = possibleHadrons[idIn];
967  int nPossHads = int(possibleHadronsNow.size());
968  if (nPossHads < 1) {
969  infoPtr->errorMsg("Error in StringFlav::pickThermal: no possible "
970  "hadrons found");
971  return 0;
972  }
973 
974  // Vector with hadron masses. Is -1.0 if m0 is use for calculating
975  // the suppression rate and mSel if mSel is used.
976  vector<double> possibleHadronMasses;
977 
978  // Calculate rates/suppression factors for given pT.
979  vector<double> rates;
980  double rateSum = 0.0;
981  for (int iHad = 0; iHad < nPossHads; iHad++) {
982  int hadronID = possibleHadronsNow[iHad].first;
983  // Pick mass and calculate suppression factor.
984  double mass = particleDataPtr->mSel(hadronID);
985  possibleHadronMasses.push_back(mass);
986  double rate = exp( -sqrt(pow2(pT)+pow2(mass))/temprNow );
987  // mT2 suppression with Gaussian pT?
988  if (mT2suppression) rate = exp( -(pow2(pT)+pow2(mass))/pow2(sigmaNow) );
989  // Multiply rate with prefactor.
990  rate *= possibleRatePrefacsNow[iHad];
991  // Save rate and add to sum
992  rates.push_back(rate);
993  rateSum += rate;
994  }
995  // Normalize rates
996  for (int iHad = 0; iHad < nPossHads; iHad++) rates[iHad] /= rateSum;
997 
998  // Get accumulated rates
999  vector<double> accumRates;
1000  for (int iHad = 0; iHad < nPossHads; iHad++) accumRates.push_back(0);
1001  for (int iHad1 = 0; iHad1 < nPossHads; iHad1++)
1002  for (int iHad2 = 0; iHad2 <= iHad1; iHad2++)
1003  accumRates[iHad1] += rates[iHad2];
1004 
1005  // Random number to decide which hadron to pick
1006  double rand = rndmPtr->flat();
1007  int hadronID = 0;
1008  int iConst = 0;
1009  double hadronMass = -1.0;
1010  for (int iHad = 0; iHad < nPossHads; iHad++) {
1011  if (rand <= accumRates[iHad]) {
1012  hadronID = possibleHadronsNow[iHad].first;
1013  iConst = possibleHadronsNow[iHad].second;
1014  hadronMass = possibleHadronMasses[iHad];
1015  break;
1016  }
1017  }
1018 
1019  // Get flavour of (di)quark to use next time.
1020  int idNext = 0;
1021  vector< pair<int,int> > constituentIDs = hadronConstIDs[hadronID];
1022  // Mesons
1023  if (particleDataPtr->isMeson(hadronID)) {
1024  int ID1 = constituentIDs[0].first;
1025  int ID2 = constituentIDs[0].second;
1026  // Special case for diagonal meson, flavour remains
1027  if (ID1 == -ID2) idNext = idIn;
1028  else idNext = (idIn == ID1 ? -ID2 : -ID1);
1029  }
1030  // Baryons
1031  else {
1032  int ID1 = constituentIDs[iConst].first;
1033  int ID2 = constituentIDs[iConst].second;
1034  if (ID1 == idIn) idNext = -ID2;
1035  if (ID2 == idIn) idNext = -ID1;
1036  }
1037 
1038  // Save new flavour and hadron.
1039  flavNew.id = -idNext; // id used to build hadron
1040  hadronIDwin = hadronID;
1041  idNewWin = idNext; // id used in next step
1042  hadronMassWin = hadronMass;
1043 
1044  // Done.
1045  return flavNew;
1046 
1047 }
1048 
1049 //--------------------------------------------------------------------------
1050 
1051 // Combine two flavours (including diquarks) to produce a hadron.
1052 // The weighting of the combination may fail, giving output 0.
1053 
1054 int StringFlav::combine(FlavContainer& flav1, FlavContainer& flav2) {
1055 
1056  // Recognize largest and smallest flavour.
1057  int id1Abs = abs(flav1.id);
1058  int id2Abs = abs(flav2.id);
1059  int idMax = max(id1Abs, id2Abs);
1060  int idMin = min(id1Abs, id2Abs);
1061 
1062  // Construct a meson.
1063  if (idMax < 9 || idMin > 1000) {
1064 
1065  // Popcorn meson: use only vertex quarks. Fail if none.
1066  if (idMin > 1000) {
1067  id1Abs = flav1.idVtx;
1068  id2Abs = flav2.idVtx;
1069  idMax = max(id1Abs, id2Abs);
1070  idMin = min(id1Abs, id2Abs);
1071  if (idMin == 0) return 0;
1072  }
1073 
1074  // Pick spin state and preliminary code.
1075  int flav = (idMax < 3) ? 0 : idMax - 2;
1076  double rndmSpin = mesonRateSum[flav] * rndmPtr->flat();
1077  int spin = -1;
1078  do rndmSpin -= mesonRate[flav][++spin];
1079  while (rndmSpin > 0.);
1080  int idMeson = 100 * idMax + 10 * idMin + mesonMultipletCode[spin];
1081 
1082  // For nondiagonal mesons distinguish particle/antiparticle.
1083  if (idMax != idMin) {
1084  int sign = (idMax%2 == 0) ? 1 : -1;
1085  if ( (idMax == id1Abs && flav1.id < 0)
1086  || (idMax == id2Abs && flav2.id < 0) ) sign = -sign;
1087  idMeson *= sign;
1088 
1089  // For light diagonal mesons include uubar - ddbar - ssbar mixing.
1090  } else if (flav < 2) {
1091  double rMix = rndmPtr->flat();
1092  if (rMix < mesonMix1[flav][spin]) idMeson = 110;
1093  else if (rMix < mesonMix2[flav][spin]) idMeson = 220;
1094  else idMeson = 330;
1095  idMeson += mesonMultipletCode[spin];
1096 
1097  // Additional suppression of eta and eta' may give failure.
1098  if (idMeson == 221 && etaSup < rndmPtr->flat()) return 0;
1099  if (idMeson == 331 && etaPrimeSup < rndmPtr->flat()) return 0;
1100  }
1101 
1102  // Finished for mesons.
1103  return idMeson;
1104  }
1105 
1106  // SU(6) factors for baryon production may give failure.
1107  int idQQ1 = idMax / 1000;
1108  int idQQ2 = (idMax / 100) % 10;
1109  int spinQQ = idMax % 10;
1110  int spinFlav = spinQQ - 1;
1111  if (spinFlav == 2 && idQQ1 != idQQ2) spinFlav = 4;
1112  if (idMin != idQQ1 && idMin != idQQ2) spinFlav++;
1113  if (spinFlav < 0 || spinFlav > 5) return 0;
1114  if (baryonCGSum[spinFlav] < rndmPtr->flat() * baryonCGMax[spinFlav])
1115  return 0;
1116 
1117  // Order quarks to form baryon. Pick spin.
1118  int idOrd1 = max( idMin, max( idQQ1, idQQ2) );
1119  int idOrd3 = min( idMin, min( idQQ1, idQQ2) );
1120  int idOrd2 = idMin + idQQ1 + idQQ2 - idOrd1 - idOrd3;
1121  int spinBar = (baryonCGSum[spinFlav] * rndmPtr->flat()
1122  < baryonCGOct[spinFlav]) ? 2 : 4;
1123 
1124  // Distinguish Lambda- and Sigma-like.
1125  bool LambdaLike = false;
1126  if (spinBar == 2 && idOrd1 > idOrd2 && idOrd2 > idOrd3) {
1127  LambdaLike = (spinQQ == 1);
1128  if (idOrd1 != idMin && spinQQ == 1) LambdaLike = (rndmPtr->flat() < 0.25);
1129  else if (idOrd1 != idMin) LambdaLike = (rndmPtr->flat() < 0.75);
1130  }
1131 
1132  // Form baryon code and return with sign.
1133  int idBaryon = (LambdaLike)
1134  ? 1000 * idOrd1 + 100 * idOrd3 + 10 * idOrd2 + spinBar
1135  : 1000 * idOrd1 + 100 * idOrd2 + 10 * idOrd3 + spinBar;
1136  return (flav1.id > 0) ? idBaryon : -idBaryon;
1137 
1138 }
1139 
1140 //--------------------------------------------------------------------------
1141 
1142 // Combine two flavours (including diquarks) to produce the lightest hadron
1143 // allowed for that flavour content. No popcorn flavours.
1144 
1145 int StringFlav::combineToLightest( int id1, int id2) {
1146 
1147  // Recognize largest and smallest flavour.
1148  int id1Abs = abs(id1);
1149  int id2Abs = abs(id2);
1150  int idMax = max(id1Abs, id2Abs);
1151  int idMin = min(id1Abs, id2Abs);
1152 
1153  // Construct a meson. Preliminary code.
1154  if (idMax < 9) {
1155  int idMeson = 100 * idMax + 10 * idMin + 1;
1156 
1157  // For nondiagonal mesons distinguish particle/antiparticle.
1158  if (idMax != idMin) {
1159  int sign = (idMax%2 == 0) ? 1 : -1;
1160  if ( (idMax == id1Abs && id1 < 0)
1161  || (idMax == id2Abs && id2 < 0) ) sign = -sign;
1162  idMeson *= sign;
1163  }
1164 
1165  // For light diagonal mesons pick pi0 or eta.
1166  else if (idMax < 3) idMeson = 111;
1167  else if (idMax == 3) idMeson = 221;
1168 
1169  // Finished for mesons.
1170  return idMeson;
1171  }
1172 
1173  // Split up diquark and order quarks.
1174  int idQQ1 = idMax / 1000;
1175  int idQQ2 = (idMax / 100) % 10;
1176  int idOrd1 = max( idMin, max( idQQ1, idQQ2) );
1177  int idOrd3 = min( idMin, min( idQQ1, idQQ2) );
1178  int idOrd2 = idMin + idQQ1 + idQQ2 - idOrd1 - idOrd3;
1179 
1180  // Create baryon. Special cases with spin 3/2 and lambdalike.
1181  int idBaryon = 1000 * idOrd1 + 100 * idOrd2 + 10 * idOrd3 + 2;
1182  if (idOrd3 == idOrd1) idBaryon += 2;
1183  else if (idOrd2 != idOrd1 && idOrd3 != idOrd2)
1184  idBaryon = 1000 * idOrd1 + 100 * idOrd3 + 10 * idOrd2 + 2;
1185 
1186  // Finished for baryons.
1187  return (id1 > 0) ? idBaryon : -idBaryon;
1188 
1189 }
1190 
1191 //--------------------------------------------------------------------------
1192 
1193 // Combine two flavours (including diquarks) to produce a hadron. Function
1194 // called in case of combining the two remaining flavours into last hadron.
1195 
1196 int StringFlav::combineLastThermal(FlavContainer& flav1, FlavContainer& flav2,
1197  double pT, double nNSP) {
1198 
1199  // Decide randomly on whether to treat flav1 or flav2 as incoming.
1200  int idIn[2] = { flav1.id, flav2.id };
1201  if (rndmPtr->flat() < 0.5) swap(idIn[0], idIn[1]);
1202  int idInNow[2] = { min(idIn[0],idIn[1]), max(idIn[0],idIn[1]) };
1203 
1204  int idInAbs = abs(idIn[0]);
1205  double temprNow = temperature;
1206  // Temperature increase to work against asymmetry. Apply for
1207  // s/c/b and diquarks.
1208  if (idInAbs > 2) temprNow *= tempPreFactor;
1209 
1210  // Enhanded-rate prefactor for MPIs and/or nearby string pieces.
1211  if (closePacking) {
1212  temprNow *= pow(max(1.0,double(infoPtr->nMPI())), exponentMPI);
1213  temprNow *= pow(max(1.0,nNSP), exponentNSP);
1214  }
1215 
1216  // Get Gaussian width in case of mT2 suppression.
1217  double sigmaNow = sigmaHad;
1218  // Prefactor for strange quarks and diquarks.
1219  if (useWidthPre) {
1220  if (abs(idInAbs) > 10) sigmaNow *= widthPreDiquark;
1221  sigmaNow *= pow(widthPreStrange,
1222  particleDataPtr->nQuarksInCode(idInAbs,3) );
1223  }
1224 
1225  // Enhanded-rate prefactor for MPIs and/or nearby string pieces.
1226  if (closePacking) {
1227  sigmaNow *= pow(max(1.0,double(infoPtr->nMPI())), exponentMPI);
1228  sigmaNow *= pow(max(1.0,nNSP), exponentNSP);
1229  }
1230 
1231  // Get the list of allowed hadrons and constituents for that combination
1232  // of (di)quarks. First parameter of pair is hadron ID, second
1233  // is nr of hadron constituents in the list.
1234  pair<int,int> inPr = pair<int,int>(idInNow[0], idInNow[1]);
1235  vector<double> possibleRatePrefacsNow = possibleRatePrefacsLast[inPr];
1236  vector< pair<int,int> > possibleHadronsNow = possibleHadronsLast[inPr];
1237  int nPossHads = int(possibleHadronsNow.size());
1238  if (nPossHads < 1) {
1239  infoPtr->errorMsg("Error in StringFlav::combineLastThermal: no "
1240  "possible hadrons found for last two");
1241  return 0;
1242  }
1243 
1244  // Vector with hadron masses. Is -1.0 if m0 is use for calculating
1245  // the suppression rate and mSel if mSel is used.
1246  vector<double> possibleHadronMasses;
1247 
1248  // Calculate rates/suppression factors for given pT.
1249  vector<double> rates;
1250  double rateSum = 0.0;
1251  for (int iHad = 0; iHad < nPossHads; iHad++) {
1252  int hadronID = possibleHadronsNow[iHad].first;
1253  // Pick mass and calculate suppression factor.
1254  double mass = particleDataPtr->mSel(hadronID);
1255  possibleHadronMasses.push_back(mass);
1256  double rate = exp( -sqrt(pow2(pT)+pow2(mass))/temprNow );
1257  // mT2 suppression with Gaussian pT?
1258  if (mT2suppression) rate = exp( -(pow2(pT)+pow2(mass))/pow2(sigmaNow) );
1259  // Multiply rate with prefactor.
1260  rate *= possibleRatePrefacsNow[iHad];
1261  // Save rate and add to sum
1262  rates.push_back(rate);
1263  rateSum += rate;
1264  }
1265  // Normalize rates
1266  for (int iHad = 0; iHad < nPossHads; iHad++) rates[iHad] /= rateSum;
1267 
1268  // Get accumulated rates
1269  vector<double> accumRates;
1270  for (int iHad = 0; iHad < nPossHads; iHad++) accumRates.push_back(0);
1271  for (int iHad1 = 0; iHad1 < nPossHads; iHad1++)
1272  for (int iHad2 = 0; iHad2 <= iHad1; iHad2++)
1273  accumRates[iHad1] += rates[iHad2];
1274 
1275  // Random number to decide which hadron to pick
1276  double rand = rndmPtr->flat();
1277  int hadronID = 0;
1278  double hadronMass = -1.0;
1279  for (int iHad = 0; iHad < nPossHads; iHad++) {
1280  if (rand <= accumRates[iHad]) {
1281  hadronID = possibleHadronsNow[iHad].first;
1282  hadronMass = possibleHadronMasses[iHad];
1283  break;
1284  }
1285  }
1286 
1287  // Save hadron.
1288  hadronIDwin = hadronID;
1289  hadronMassWin = hadronMass;
1290 
1291  // Done.
1292  return hadronIDwin;
1293 }
1294 
1295 //--------------------------------------------------------------------------
1296 
1297 // Assign popcorn quark inside an original (= rank 0) diquark.
1298 
1299 void StringFlav::assignPopQ(FlavContainer& flav) {
1300 
1301  // Safety check that intended to do something.
1302  int idAbs = abs(flav.id);
1303  if (flav.rank > 0 || idAbs < 1000) return;
1304 
1305  // Make choice of popcorn quark.
1306  int id1 = (idAbs/1000)%10;
1307  int id2 = (idAbs/100)%10;
1308  double pop2WT = 1.;
1309  if (id1 == 3) pop2WT = scbBM[1];
1310  else if (id1 > 3) pop2WT = scbBM[2];
1311  if (id2 == 3) pop2WT /= scbBM[1];
1312  else if (id2 > 3) pop2WT /= scbBM[2];
1313  // Agrees with Patrik code, but opposite to intention??
1314  flav.idPop = ((1. + pop2WT) * rndmPtr->flat() > 1.) ? id2 : id1;
1315  flav.idVtx = id1 + id2 - flav.idPop;
1316 
1317  // Also determine if to produce popcorn meson.
1318  flav.nPop = 0;
1319  double popWT = popS[0];
1320  if (id1 == 3) popWT = popS[1];
1321  if (id2 == 3) popWT = popS[2];
1322  if (idAbs%10 == 1) popWT *= sqrt(probQQ1toQQ0);
1323  if ((1. + popWT) * rndmPtr->flat() > 1.) flav.nPop = 1;
1324 
1325 }
1326 
1327 //--------------------------------------------------------------------------
1328 
1329 // Combine two quarks to produce a diquark.
1330 // Normally according to production composition, but nonvanishing idHad
1331 // means diquark from known hadron content, so use SU(6) wave function.
1332 
1333 int StringFlav::makeDiquark(int id1, int id2, int idHad) {
1334 
1335  // Initial values.
1336  int idMin = min( abs(id1), abs(id2));
1337  int idMax = max( abs(id1), abs(id2));
1338  int spin = 1;
1339 
1340  // Select spin of diquark formed from two valence quarks in proton.
1341  // (More hadron cases??)
1342  if (abs(idHad) == 2212 || abs(idHad) == 2112) {
1343  if (idMin == 1 && idMax == 2 && rndmPtr->flat() < 0.75) spin = 0;
1344 
1345  // Else select spin of diquark according to assumed spin-1 suppression.
1346  } else if (idMin != idMax) {
1347  if (rndmPtr->flat() > probQQ1join[min(idMax,5) - 2]) spin = 0;
1348  }
1349 
1350  // Combined diquark code.
1351  int idNewAbs = 1000 * idMax + 100 * idMin + 2 * spin + 1;
1352  return (id1 > 0) ? idNewAbs : -idNewAbs;
1353 
1354 }
1355 
1356 //==========================================================================
1357 
1358 // The StringZ class.
1359 
1360 //--------------------------------------------------------------------------
1361 
1362 // Constants: could be changed here if desired, but normally should not.
1363 // These are of technical nature, as described for each.
1364 
1365 // When a or c are close to special cases, default to these.
1366 const double StringZ::CFROMUNITY = 0.01;
1367 const double StringZ::AFROMZERO = 0.02;
1368 const double StringZ::AFROMC = 0.01;
1369 
1370 // Do not take exponent of too large or small number.
1371 const double StringZ::EXPMAX = 50.;
1372 
1373 //--------------------------------------------------------------------------
1374 
1375 // Initialize data members of the string z selection.
1376 
1377 void StringZ::init() {
1378 
1379  // c and b quark masses.
1380  mc2 = pow2( particleDataPtr->m0(4));
1381  mb2 = pow2( particleDataPtr->m0(5));
1382 
1383  // Paramaters of Lund/Bowler symmetric fragmentation function.
1384  aLund = parm("StringZ:aLund");
1385  bLund = parm("StringZ:bLund");
1386  aExtraSQuark = parm("StringZ:aExtraSQuark");
1387  aExtraDiquark = parm("StringZ:aExtraDiquark");
1388  rFactC = parm("StringZ:rFactC");
1389  rFactB = parm("StringZ:rFactB");
1390  rFactH = parm("StringZ:rFactH");
1391 
1392  // Alternative parameterisation of Lund FF using average z(rho) instead of b.
1393  if (flag("StringZ:deriveBLund")) {
1394  if (!deriveBLund()) {
1395  infoPtr->errorMsg("Error in StringZ::init: Derivation of b parameter "
1396  " failed. Reverting to default.");
1397  settingsPtr->resetParm("StringZ:bLund");
1398  }
1399  }
1400 
1401  // Flags and parameters of nonstandard Lund fragmentation functions.
1402  useNonStandC = flag("StringZ:useNonstandardC");
1403  useNonStandB = flag("StringZ:useNonstandardB");
1404  useNonStandH = flag("StringZ:useNonstandardH");
1405  aNonC = parm("StringZ:aNonstandardC");
1406  aNonB = parm("StringZ:aNonstandardB");
1407  aNonH = parm("StringZ:aNonstandardH");
1408  bNonC = parm("StringZ:bNonstandardC");
1409  bNonB = parm("StringZ:bNonstandardB");
1410  bNonH = parm("StringZ:bNonstandardH");
1411 
1412  // Flags and parameters of Peterson/SLAC fragmentation function.
1413  usePetersonC = flag("StringZ:usePetersonC");
1414  usePetersonB = flag("StringZ:usePetersonB");
1415  usePetersonH = flag("StringZ:usePetersonH");
1416  epsilonC = parm("StringZ:epsilonC");
1417  epsilonB = parm("StringZ:epsilonB");
1418  epsilonH = parm("StringZ:epsilonH");
1419 
1420  // Parameters for joining procedure.
1421  stopM = parm("StringFragmentation:stopMass");
1422  stopNF = parm("StringFragmentation:stopNewFlav");
1423  stopS = parm("StringFragmentation:stopSmear");
1424 
1425 }
1426 
1427 //--------------------------------------------------------------------------
1428 
1429 // Alternative parameterisation of the Lund function. Derive the bLund
1430 // parameter given the average z for fixed a and mT2.
1431 
1432 bool StringZ::deriveBLund() {
1433 
1434  // Set up using reference mT2 = mRho^2 + 2*sigmaPT^2
1435  double mRef = particleDataPtr->m0(113);
1436  double mT2ref = pow2(mRef) + 2.*pow2(parm("stringPT:sigma"));
1437  double avgZ = parm("StringZ:avgZLund");
1438  double a = parm("StringZ:aLund");
1439 
1440  // Define lundFF as a function of only b, fixing a, c and mT2 as parameters
1441  auto lundFF = [=](double b) { return LundFFAvg(a, b, 1., mT2ref, 1.e-6); };
1442 
1443  // Solve for b
1444  double bNow;
1445  bool check = brent(bNow, lundFF, avgZ, 0.01, 20.0, 1.e-6);
1446 
1447  // Check if derived b fell inside the nominal range for bLund
1448  if (check) {
1449  settingsPtr->parm("StringZ:bLund", bNow, false);
1450 
1451  // Print out derived value for b (and mT2ref), noting if outside range.
1452  stringstream msg;
1453  msg << fixed << setprecision(2) << "\n <z(rho)> = " << setw(5)
1454  << avgZ << " for aLund = "<< a <<" & mT2ref = " << setw(5) << mT2ref
1455  << " GeV^2 gave bLund = " << setw(5) << bNow << " GeV^-2:";
1456  if ( bNow == parm("StringZ:bLund") ) {
1457  if (!settingsPtr->parm("Print:quiet"))
1458  cout << msg.str() << " accepted" << endl;
1459  } else {
1460  // If outside range, tell user but force anyway so fits can see
1461  // behaviour.
1462  msg << " accepted (forced)";
1463  infoPtr->errorMsg("Warning in StringZ::deriveBLund", msg.str());
1464  settingsPtr->parm("StringZ:bLund", bNow, true);
1465  }
1466 
1467  // No further calls needed since b parameter updated in settings database.
1468  settingsPtr->flag("StringZ:deriveBLund", false);
1469  }
1470  return check;
1471 }
1472 
1473 //--------------------------------------------------------------------------
1474 
1475 // Generate the fraction z that the next hadron will take,
1476 // using either Lund/Bowler or, for heavy, Peterson/SLAC functions.
1477 // Note: for a heavy new coloured particle we assume pT negligible.
1478 
1479 double StringZ::zFrag( int idOld, int idNew, double mT2) {
1480 
1481  // Find if old or new flavours correspond to diquarks.
1482  int idOldAbs = abs(idOld);
1483  int idNewAbs = abs(idNew);
1484  bool isOldSQuark = (idOldAbs == 3);
1485  bool isNewSQuark = (idNewAbs == 3);
1486  bool isOldDiquark = (idOldAbs > 1000 && idOldAbs < 10000);
1487  bool isNewDiquark = (idNewAbs > 1000 && idNewAbs < 10000);
1488 
1489  // Find heaviest quark in fragmenting parton/diquark.
1490  int idFrag = idOldAbs;
1491  if (isOldDiquark) idFrag = max( idOldAbs / 1000, (idOldAbs / 100) % 10);
1492 
1493  // Use Peterson where explicitly requested for heavy flavours.
1494  if (idFrag == 4 && usePetersonC) return zPeterson( epsilonC);
1495  if (idFrag == 5 && usePetersonB) return zPeterson( epsilonB);
1496  if (idFrag > 5 && usePetersonH) {
1497  double epsilon = epsilonH * mb2 / mT2;
1498  return zPeterson( epsilon);
1499  }
1500 
1501  // Nonstandard a and b values implemented for heavy flavours.
1502  double aNow = aLund;
1503  double bNow = bLund;
1504  if (idFrag == 4 && useNonStandC) {
1505  aNow = aNonC;
1506  bNow = bNonC;
1507  } else if (idFrag == 5 && useNonStandB) {
1508  aNow = aNonB;
1509  bNow = bNonB;
1510  } else if (idFrag > 5 && useNonStandH) {
1511  aNow = aNonH;
1512  bNow = bNonH;
1513  }
1514 
1515  // Shape parameters of Lund symmetric fragmentation function.
1516  double aShape = aNow;
1517  if (isOldSQuark) aShape += aExtraSQuark;
1518  if (isOldDiquark) aShape += aExtraDiquark;
1519  double bShape = bNow * mT2;
1520  double cShape = 1.;
1521  if (isOldSQuark) cShape -= aExtraSQuark;
1522  if (isNewSQuark) cShape += aExtraSQuark;
1523  if (isOldDiquark) cShape -= aExtraDiquark;
1524  if (isNewDiquark) cShape += aExtraDiquark;
1525  if (idFrag == 4) cShape += rFactC * bNow * mc2;
1526  if (idFrag == 5) cShape += rFactB * bNow * mb2;
1527  if (idFrag > 5) cShape += rFactH * bNow * mT2;
1528  return zLund( aShape, bShape, cShape);
1529 
1530 }
1531 
1532 //--------------------------------------------------------------------------
1533 
1534 // Generate a random z according to the Lund/Bowler symmetric
1535 // fragmentation function f(z) = (1 -z)^a * exp(-b/z) / z^c.
1536 // Normalized so that f(z_max) = 1 it can also be written as
1537 // f(z) = exp( a * ln( (1 - z) / (1 - z_max) ) + b * (1/z_max - 1/z)
1538 // + c * ln(z_max/z) ).
1539 
1540 double StringZ::zLund( double a, double b, double c) {
1541 
1542  // Special cases for c = 1, a = 0 and a = c.
1543  bool cIsUnity = (abs( c - 1.) < CFROMUNITY);
1544  bool aIsZero = (a < AFROMZERO);
1545  bool aIsC = (abs(a - c) < AFROMC);
1546 
1547  // Determine position of maximum.
1548  double zMax;
1549  if (aIsZero) zMax = (c > b) ? b / c : 1.;
1550  else if (aIsC) zMax = b / (b + c);
1551  else { zMax = 0.5 * (b + c - sqrt( pow2(b - c) + 4. * a * b)) / (c - a);
1552  if (zMax > 0.9999 && b > 100.) zMax = min(zMax, 1. - a / b); }
1553 
1554  // Subdivide z range if distribution very peaked near either endpoint.
1555  bool peakedNearZero = (zMax < 0.1);
1556  bool peakedNearUnity = (zMax > 0.85 && b > 1.);
1557 
1558  // Find integral of trial function everywhere bigger than f.
1559  // (Dummy start values.)
1560  double fIntLow = 1.;
1561  double fIntHigh = 1.;
1562  double fInt = 2.;
1563  double zDiv = 0.5;
1564  double zDivC = 0.5;
1565  // When z_max is small use that f(z)
1566  // < 1 for z < z_div = 2.75 * z_max,
1567  // < (z_div/z)^c for z > z_div (=> logarithm for c = 1, else power).
1568  if (peakedNearZero) {
1569  zDiv = 2.75 * zMax;
1570  fIntLow = zDiv;
1571  if (cIsUnity) fIntHigh = -zDiv * log(zDiv);
1572  else { zDivC = pow( zDiv, 1. - c);
1573  fIntHigh = zDiv * (1. - 1./zDivC) / (c - 1.);}
1574  fInt = fIntLow + fIntHigh;
1575  // When z_max large use that f(z)
1576  // < exp( b * (z - z_div) ) for z < z_div with z_div messy expression,
1577  // < 1 for z > z_div.
1578  // To simplify expressions the integral is extended to z = -infinity.
1579  } else if (peakedNearUnity) {
1580  double rcb = sqrt(4. + pow2(c / b));
1581  zDiv = rcb - 1./zMax - (c / b) * log( zMax * 0.5 * (rcb + c / b) );
1582  if (!aIsZero) zDiv += (a/b) * log(1. - zMax);
1583  zDiv = min( zMax, max(0., zDiv));
1584  fIntLow = 1. / b;
1585  fIntHigh = 1. - zDiv;
1586  fInt = fIntLow + fIntHigh;
1587  }
1588 
1589  // Choice of z, preweighted for peaks at low or high z. (Dummy start values.)
1590  double z = 0.5;
1591  double fPrel = 1.;
1592  double fVal = 1.;
1593  do {
1594  // Choice of z flat good enough for distribution peaked in the middle;
1595  // if not this z can be reused as a random number in general.
1596  z = rndmPtr->flat();
1597  fPrel = 1.;
1598  // When z_max small use flat below z_div and 1/z^c above z_div.
1599  if (peakedNearZero) {
1600  if (fInt * rndmPtr->flat() < fIntLow) z = zDiv * z;
1601  else if (cIsUnity) {z = pow( zDiv, z); fPrel = zDiv / z;}
1602  else { z = pow( zDivC + (1. - zDivC) * z, 1. / (1. - c) );
1603  fPrel = pow( zDiv / z, c); }
1604  // When z_max large use exp( b * (z -z_div) ) below z_div
1605  // and flat above it.
1606  } else if (peakedNearUnity) {
1607  if (fInt * rndmPtr->flat() < fIntLow) {
1608  z = zDiv + log(z) / b;
1609  fPrel = exp( b * (z - zDiv) );
1610  } else z = zDiv + (1. - zDiv) * z;
1611  }
1612 
1613  // Evaluate actual f(z) (if in physical range) and correct.
1614  if (z > 0 && z < 1) {
1615  double fExp = b * (1. / zMax - 1. / z)+ c * log(zMax / z);
1616  if (!aIsZero) fExp += a * log( (1. - z) / (1. - zMax) );
1617  fVal = exp( max( -EXPMAX, min( EXPMAX, fExp) ) ) ;
1618  } else fVal = 0.;
1619  } while (fVal < rndmPtr->flat() * fPrel);
1620 
1621  // Done.
1622  return z;
1623 
1624 }
1625 
1626 //--------------------------------------------------------------------------
1627 
1628 // Generate a random z according to the Peterson/SLAC formula
1629 // f(z) = 1 / ( z * (1 - 1/z - epsilon/(1-z))^2 )
1630 // = z * (1-z)^2 / ((1-z)^2 + epsilon * z)^2.
1631 
1632 double StringZ::zPeterson( double epsilon) {
1633 
1634  double z, fVal;
1635 
1636  // For large epsilon pick z flat and reject,
1637  // knowing that 4 * epsilon * f(z) < 1 everywhere.
1638  if (epsilon > 0.01) {
1639  do {
1640  z = rndmPtr->flat();
1641  fVal = 4. * epsilon * z * pow2(1. - z)
1642  / pow2( pow2(1. - z) + epsilon * z);
1643  } while (fVal < rndmPtr->flat());
1644  return z;
1645  }
1646 
1647  // Else split range, using that 4 * epsilon * f(z)
1648  // < 4 * epsilon / (1 - z)^2 for 0 < z < 1 - 2 * sqrt(epsilon)
1649  // < 1 for 1 - 2 * sqrt(epsilon) < z < 1
1650  double epsRoot = sqrt(epsilon);
1651  double epsComb = 0.5 / epsRoot - 1.;
1652  double fIntLow = 4. * epsilon * epsComb;
1653  double fInt = fIntLow + 2. * epsRoot;
1654  do {
1655  if (rndmPtr->flat() * fInt < fIntLow) {
1656  z = 1. - 1. / (1. + rndmPtr->flat() * epsComb);
1657  fVal = z * pow2( pow2(1. - z) / (pow2(1. - z) + epsilon * z) );
1658  } else {
1659  z = 1. - 2. * epsRoot * rndmPtr->flat();
1660  fVal = 4. * epsilon * z * pow2(1. - z)
1661  / pow2( pow2(1. - z) + epsilon * z);
1662  }
1663  } while (fVal < rndmPtr->flat());
1664  return z;
1665 
1666 }
1667 
1668 //==========================================================================
1669 
1670 // The StringPT class.
1671 
1672 //--------------------------------------------------------------------------
1673 
1674 // Constants: could be changed here if desired, but normally should not.
1675 // These are of technical nature, as described for each.
1676 
1677 // To avoid division by zero one must have sigma > 0.
1678 const double StringPT::SIGMAMIN = 0.2;
1679 
1680 //--------------------------------------------------------------------------
1681 
1682 // Initialize data members of the string pT selection.
1683 
1684 void StringPT::init() {
1685 
1686  // Parameters of the pT width and enhancement.
1687  double sigma = parm("StringPT:sigma");
1688  sigmaQ = sigma / sqrt(2.);
1689  enhancedFraction = parm("StringPT:enhancedFraction");
1690  enhancedWidth = parm("StringPT:enhancedWidth");
1691  widthPreStrange = parm("StringPT:widthPreStrange");
1692  widthPreDiquark = parm("StringPT:widthPreDiquark");
1693  useWidthPre = (widthPreStrange > 1.0) || (widthPreDiquark > 1.0);
1694 
1695  // Temperature for thermal model.
1696  thermalModel = flag("StringPT:thermalModel");
1697  temperature = parm("StringPT:temperature");
1698  tempPreFactor = parm("StringPT:tempPreFactor");
1699 
1700  // Upper estimate of thermal spectrum: fraction at x = pT_quark/T < 1.
1701  fracSmallX = 0.6 / (0.6 + (1.2/0.9) * exp(-0.9));
1702 
1703  // Enhanded-width prefactor for MPIs and/or nearby string pieces.
1704  closePacking = flag("StringPT:closePacking");
1705  exponentMPI = parm("StringPT:expMPI");
1706  exponentNSP = parm("StringPT:expNSP");
1707 
1708  // Parameter for pT suppression in MiniStringFragmentation.
1709  sigma2Had = 2. * pow2( max( SIGMAMIN, sigma) );
1710 
1711 }
1712 
1713 //--------------------------------------------------------------------------
1714 
1715 // Generate quark pT according to fitting functions, such that
1716 // hadron pT is generated according to exp(-pT/T) d^2pT.
1717 
1718 pair<double, double> StringPT::pxyThermal(int idIn, double nNSP) {
1719 
1720  double temprNow = temperature;
1721  // Temperature increase to work against asymmetry. Apply for
1722  // s/c/b and diquarks.
1723  if (abs(idIn) > 2) temprNow *= tempPreFactor;
1724 
1725  // Enhanded-width prefactor for MPIs and/or nearby string pieces.
1726  if (closePacking) {
1727  temprNow *= pow(max(1.0,double(infoPtr->nMPI())), exponentMPI);
1728  temprNow *= pow(max(1.0,nNSP), exponentNSP);
1729  }
1730 
1731  // Pick x = pT_quark/T according to K_{1/4}(x)/x^{1/4} * x dx.
1732  double xrand, approx, wanted;
1733  do {
1734  xrand = (rndmPtr->flat() < fracSmallX) ? rndmPtr->flat()
1735  : 1. - log(rndmPtr->flat()) / 0.9;
1736  approx = (xrand < 1.) ? 0.6 : 1.2 * exp(-0.9 * xrand);
1737  wanted = BesselK14(xrand) * pow( xrand, 0.75);
1738  } while (rndmPtr->flat() * approx > wanted);
1739 
1740  // Find pT_quark. Random number to decide on angle.
1741  double pTquark = xrand * temprNow;
1742  double phi = 2.0 * M_PI * rndmPtr->flat();
1743 
1744  // Done.
1745  return pair<double, double>( pTquark * cos(phi), pTquark * sin(phi) );
1746 
1747 }
1748 
1749 //--------------------------------------------------------------------------
1750 
1751 // Generate Gaussian pT such that <p_x^2> = <p_x^2> = sigma^2 = width^2/2,
1752 // but with small fraction multiplied up to a broader spectrum.
1753 
1754 pair<double, double> StringPT::pxyGauss(int idIn, double nNSP) {
1755 
1756  // Normal (classical) width selection.
1757  double sigma = sigmaQ;
1758  if (rndmPtr->flat() < enhancedFraction) sigma *= enhancedWidth;
1759 
1760  // Prefactor for strange quarks and diquarks.
1761  if (useWidthPre) {
1762  if (abs(idIn) > 10) sigma *= widthPreDiquark;
1763  sigma *= pow(widthPreStrange, particleDataPtr->nQuarksInCode(idIn, 3) );
1764  }
1765 
1766  // Enhanded-width prefactor for MPIs and/or nearby string pieces.
1767  if (closePacking) {
1768  sigma *= pow(max(1.0,double(infoPtr->nMPI())), exponentMPI);
1769  sigma *= pow(max(1.0,nNSP), exponentNSP);
1770  }
1771 
1772  // Generate (p_x, p_y) pair.
1773  pair<double, double> gauss2 = rndmPtr->gauss2();
1774  return pair<double, double>(sigma * gauss2.first, sigma * gauss2.second);
1775 
1776 }
1777 
1778 //--------------------------------------------------------------------------
1779 
1780 // Evaluate Bessel function K_{1/4}(x).
1781 // Use power series for x < 2.5 and asymptotic expansion for x > 2.5.
1782 // Number of terms picked to have accuracy better than 1 per mille.
1783 // Based on M. Abramowitz and I.A. Stegun, eqs. 9.6.2, 9.6.10, 9.7.2.
1784 
1785 double StringPT::BesselK14(double x) {
1786 
1787  // Power series expansion of K_{1/4} : k = 0 term.
1788  if (x < 2.5) {
1789  double xRat = 0.25 * x * x;
1790  double prodP = pow( 0.5 * x, -0.25) / 1.2254167024;
1791  double prodN = pow( 0.5 * x, 0.25) / 0.9064024771;
1792  double sum = prodP - prodN;
1793 
1794  // Power series expansion of K_{1/4} : m > 0 terms.
1795  for (int k = 1; k < 6; ++k) {
1796  prodP *= xRat / (k * (k - 0.25));
1797  prodN *= xRat / (k * (k + 0.25));
1798  sum += prodP - prodN;
1799  }
1800  sum *= M_PI * sqrt(0.5);
1801  return sum;
1802 
1803  // Asymptotic expansion of K_{1/4}.
1804  } else {
1805  double asym = sqrt(M_PI * 0.5 / x) * exp(-x);
1806  double term1 = - 0.75 / ( 8. * x);
1807  double term2 = -term1 * 8.75 / (16. * x);
1808  double term3 = -term2 * 24.75 / (24. * x);
1809  double term4 = -term3 * 48.75 / (32. * x);
1810  asym *= 1. + term1 + term2 + term3 + term4;
1811  return asym;
1812  }
1813 }
1814 
1815 //==========================================================================
1816 
1817 } // end namespace Pythia8
Generated by   doxygen 1.8.5