source: HiSusy/trunk/Pythia8/pythia8170/src/FragmentationFlavZpT.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

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