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