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