PYTHIA8/pythia8145/src/FragmentationFlavZpT.cxx

   1 // FragmentationFlavZpT.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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
  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   // 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
 260 FlavContainer 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
 354 int 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.
 366     if (idMin > 1000) {
 367       id1Abs = flav1.idVtx;
 368       id2Abs = flav2.idVtx;
 369       idMax = max(id1Abs, id2Abs);
 370       idMin = min(id1Abs, id2Abs);
 371     }
 372
 373     // Pick spin state and preliminary code.
 374     int flav = (idMax < 3) ? 0 : idMax - 2;
 375     double rndmSpin = mesonRateSum[flav] * rndmPtr->flat();
 376     int spin = -1;
 377     do rndmSpin -= mesonRate[flav][++spin];
 378     while (rndmSpin > 0.);
 379     int idMeson = 100 * idMax + 10 * idMin + mesonMultipletCode[spin];
 380
 381     // For nondiagonal mesons distinguish particle/antiparticle.
 382     if (idMax != idMin) {
 383       int sign = (idMax%2 == 0) ? 1 : -1;
 384       if ( (idMax == id1Abs && flav1.id < 0)
 385         || (idMax == id2Abs && flav2.id < 0) ) sign = -sign;
 386       idMeson *= sign;
 387
 388     // For light diagonal mesons include uubar - ddbar - ssbar mixing.
 389     } else if (flav < 2) {
 390       double rMix = rndmPtr->flat();
 391       if      (rMix < mesonMix1[flav][spin]) idMeson = 110;
 392       else if (rMix < mesonMix2[flav][spin]) idMeson = 220;
 393       else                                   idMeson = 330;
 394       idMeson += mesonMultipletCode[spin];
 395
 396       // Additional suppression of eta and eta' may give failure.
 397       if (idMeson == 221 && etaSup < rndmPtr->flat()) return 0;
 398       if (idMeson == 331 && etaPrimeSup < rndmPtr->flat()) return 0;
 399     }
 400
 401     // Finished for mesons.
 402     return idMeson;
 403   }
 404
 405   // SU(6) factors for baryon production may give failure.
 406   int idQQ1 = idMax / 1000;
 407   int idQQ2 = (idMax / 100) % 10;
 408   int spinQQ = idMax % 10;
 409   int spinFlav = spinQQ - 1;
 410   if (spinFlav == 2 && idQQ1 != idQQ2) spinFlav = 4;
 411   if (idMin != idQQ1 && idMin != idQQ2) spinFlav++;
 412   if (baryonCGSum[spinFlav] < rndmPtr->flat() * baryonCGMax[spinFlav])
 413     return 0;
 414
 415   // Order quarks to form baryon. Pick spin.
 416   int idOrd1 = max( idMin, max( idQQ1, idQQ2) );
 417   int idOrd3 = min( idMin, min( idQQ1, idQQ2) );
 418   int idOrd2 = idMin + idQQ1 + idQQ2 - idOrd1 - idOrd3;
 419   int spinBar = (baryonCGSum[spinFlav] * rndmPtr->flat()
 420     < baryonCGOct[spinFlav]) ? 2 : 4;
 421
 422   // Distinguish Lambda- and Sigma-like.
 423   bool LambdaLike = false;
 424   if (spinBar == 2 && idOrd1 > idOrd2 && idOrd2 > idOrd3) {
 425     LambdaLike = (spinQQ == 1);
 426     if (idOrd1 != idMin && spinQQ == 1) LambdaLike = (rndmPtr->flat() < 0.25);
 427     else if (idOrd1 != idMin)           LambdaLike = (rndmPtr->flat() < 0.75);
 428   }
 429
 430   // Form baryon code and return with sign.
 431   int idBaryon = (LambdaLike)
 432     ? 1000 * idOrd1 + 100 * idOrd3 + 10 * idOrd2 + spinBar
 433     : 1000 * idOrd1 + 100 * idOrd2 + 10 * idOrd3 + spinBar;
 434    return (flav1.id > 0) ? idBaryon : -idBaryon;
 435
 436 }
 437
 438 //--------------------------------------------------------------------------
 439
 440 // Assign popcorn quark inside an original (= rank 0) diquark.
 441
 442 void StringFlav::assignPopQ(FlavContainer& flav) {
 443
 444   // Safety check that intended to do something.
 445   int idAbs = abs(flav.id);
 446   if (flav.rank > 0 || idAbs < 1000) return;
 447
 448   // Make choice of popcorn quark.
 449   int id1 = (idAbs/1000)%10;
 450   int id2 = (idAbs/100)%10;
 451   double pop2WT = 1.;
 452        if (id1 == 3) pop2WT = scbBM[1];
 453   else if (id1 >  3) pop2WT = scbBM[2];
 454        if (id2 == 3) pop2WT /= scbBM[1];
 455   else if (id2 >  3) pop2WT /= scbBM[2];
 456   // Agrees with Patrik code, but opposite to intention??
 457   flav.idPop = ((1. + pop2WT) * rndmPtr->flat() > 1.) ? id2 : id1;
 458   flav.idVtx = id1 + id2 - flav.idPop;
 459
 460   // Also determine if to produce popcorn meson.
 461   flav.nPop = 0;
 462   double popWT = popS[0];
 463   if (id1 == 3) popWT = popS[1];
 464   if (id2 == 3) popWT = popS[2];
 465   if (idAbs%10 == 1) popWT *= sqrt(probQQ1toQQ0);
 466   if ((1. + popWT) * rndmPtr->flat() > 1.) flav.nPop = 1;
 467
 468 }
 469
 470 //--------------------------------------------------------------------------
 471
 472 // Combine two quarks to produce a diquark.
 473 // Normally according to production composition, but nonvanishing idHad
 474 // means diquark from known hadron content, so use SU(6) wave fucntion.
 475
 476 int StringFlav::makeDiquark(int id1, int id2, int idHad) {
 477
 478   // Initial values.
 479   int idMin = min( abs(id1), abs(id2));
 480   int idMax = max( abs(id1), abs(id2));
 481   int spin = 1;
 482
 483   // Select spin of diquark formed from two valence quarks in proton.
 484   // (More hadron cases??)
 485   if (abs(idHad) == 2212) {
 486     if (idMin == 1 && idMax == 2 && rndmPtr->flat() < 0.75) spin = 0;
 487
 488   // Else select spin of diquark according to production composition.
 489   } else {
 490     if (idMin != idMax && rndmPtr->flat() > probQQ1norm) spin = 0;
 491   }
 492
 493   // Combined diquark code.
 494   int idNewAbs = 1000 * idMax + 100 * idMin + 2 * spin + 1;
 495   return (id1 > 0) ? idNewAbs : -idNewAbs;
 496
 497 }
 498
 499 //==========================================================================
 500
 501 // The StringZ class.
 502
 503 //--------------------------------------------------------------------------
 504
 505 // Constants: could be changed here if desired, but normally should not.
 506 // These are of technical nature, as described for each.
 507
 508 // When a or c are close to special cases, default to these.
 509 const double StringZ::CFROMUNITY = 0.01;
 510 const double StringZ::AFROMZERO  = 0.02;
 511 const double StringZ::AFROMC     = 0.01;
 512
 513 // Do not take exponent of too large or small number.
 514 const double StringZ::EXPMAX     = 50.;
 515
 516 //--------------------------------------------------------------------------
 517
 518 // Initialize data members of the string z selection.
 519
 520 void StringZ::init(Settings& settings, ParticleData& particleData,
 521   Rndm* rndmPtrIn) {
 522
 523   // Save pointer.
 524   rndmPtr        = rndmPtrIn;
 525
 526   // c and b quark masses.
 527   mc2 = pow2( particleData.m0(4));
 528   mb2 = pow2( particleData.m0(5));
 529
 530   // Paramaters of Lund/Bowler symmetric fragmentation function.
 531   aLund = settings.parm("StringZ:aLund");
 532   bLund = settings.parm("StringZ:bLund");
 533   aExtraDiquark = settings.parm("StringZ:aExtraDiquark");
 534   rFactC = settings.parm("StringZ:rFactC");
 535   rFactB = settings.parm("StringZ:rFactB");
 536   rFactH = settings.parm("StringZ:rFactH");
 537
 538   // Flags and parameters of Peterson/SLAC fragmentation function.
 539   usePetersonC = settings.flag("StringZ:usePetersonC");
 540   usePetersonB = settings.flag("StringZ:usePetersonB");
 541   usePetersonH = settings.flag("StringZ:usePetersonH");
 542   epsilonC = settings.parm("StringZ:epsilonC");
 543   epsilonB = settings.parm("StringZ:epsilonB");
 544   epsilonH = settings.parm("StringZ:epsilonH");
 545
 546 }
 547
 548 //--------------------------------------------------------------------------
 549
 550 // Generate the fraction z that the next hadron will take,
 551 // using either Lund/Bowler or, for heavy, Peterson/SLAC functions.
 552 // Note: for a heavy new coloured particle we assume pT negligible.
 553
 554 double StringZ::zFrag( int idOld, int idNew, double mT2) {
 555
 556   // Find if old or new flavours correspond to diquarks.
 557   int idOldAbs = abs(idOld);
 558   int idNewAbs = abs(idNew);
 559   bool isOldDiquark = (idOldAbs > 1000 && idOldAbs < 10000);
 560   bool isNewDiquark = (idNewAbs > 1000 && idNewAbs < 10000);
 561
 562   // Find heaviest quark in fragmenting parton/diquark.
 563   int idFrag = idOldAbs;
 564   if (isOldDiquark) idFrag = max( idOldAbs / 1000, (idOldAbs / 100) % 10);
 565
 566   // Use Peterson where explicitly requested for heavy flavours.
 567   if (idFrag == 4 && usePetersonC) return zPeterson( epsilonC);
 568   if (idFrag == 5 && usePetersonB) return zPeterson( epsilonB);
 569   if (idFrag >  5 && usePetersonH) {
 570     double epsilon = epsilonH * mb2 / mT2;
 571     return zPeterson( epsilon);
 572   }
 573
 574   // Shape parameters of Lund symmetric fragmentation function.
 575   double aShape = aLund;
 576   if (isOldDiquark) aShape += aExtraDiquark;
 577   double bShape = bLund * mT2;
 578   double cShape = 1.;
 579   if (isOldDiquark) cShape -= aExtraDiquark;
 580   if (isNewDiquark) cShape += aExtraDiquark;
 581   if (idFrag == 4) cShape += rFactC * bLund * mc2;
 582   if (idFrag == 5) cShape += rFactB * bLund * mb2;
 583   if (idFrag >  5) cShape += rFactH * bLund * mT2;
 584   return zLund( aShape, bShape, cShape);
 585
 586 }
 587
 588 //--------------------------------------------------------------------------
 589
 590 // Generate a random z according to the Lund/Bowler symmetric
 591 // fragmentation function f(z) = (1 -z)^a * exp(-b/z) / z^c.
 592 // Normalized so that f(z_max) = 1  it can also be written as
 593 // f(z) = exp( a * ln( (1 - z) / (1 - z_max) ) + b * (1/z_max - 1/z)
 594 //           + c * ln(z_max/z) ).
 595
 596 double StringZ::zLund( double a, double b, double c) {
 597
 598   // Special cases for c = 1, a = 0 and a = c.
 599   bool cIsUnity = (abs( c - 1.) < CFROMUNITY);
 600   bool aIsZero = (a < AFROMZERO);
 601   bool aIsC = (abs(a - c) < AFROMC);
 602
 603   // Determine position of maximum.
 604   double zMax;
 605   if (aIsZero) zMax = (c > b) ? b / c : 1.;
 606   else if (aIsC) zMax = b / (b + c);
 607   else { zMax = 0.5 * (b + c - sqrt( pow2(b - c) + 4. * a * b)) / (c - a);
 608          if (zMax > 0.9999 && b > 100.) zMax = min(zMax, 1. - a / b); }
 609
 610   // Subdivide z range if distribution very peaked near either endpoint.
 611   bool peakedNearZero = (zMax < 0.1);
 612   bool peakedNearUnity = (zMax > 0.85 && b > 1.);
 613
 614   // Find integral of trial function everywhere bigger than f.
 615   // (Dummy start values.)
 616   double fIntLow = 1.;
 617   double fIntHigh = 1.;
 618   double fInt = 2.;
 619   double zDiv = 0.5;
 620   double zDivC = 0.5;
 621   // When z_max is small use that f(z)
 622   //   < 1     for z < z_div = 2.75 * z_max,
 623   //   < (z_div/z)^c for z > z_div (=> logarithm for c = 1, else power).
 624   if (peakedNearZero) {
 625     zDiv = 2.75 * zMax;
 626     fIntLow = zDiv;
 627     if (cIsUnity) fIntHigh = -zDiv * log(zDiv);
 628     else { zDivC = pow( zDiv, 1. - c);
 629            fIntHigh = zDiv * (1. - 1./zDivC) / (c - 1.);}
 630     fInt = fIntLow + fIntHigh;
 631   // When z_max large use that f(z)
 632   //   < exp( b * (z - z_div) ) for z < z_div with z_div messy expression,
 633   //   < 1   for z > z_div.
 634   // To simplify expressions the integral is extended to z =  -infinity.
 635   } else if (peakedNearUnity) {
 636     double rcb = sqrt(4. + pow2(c / b));
 637     zDiv = rcb - 1./zMax - (c / b) * log( zMax * 0.5 * (rcb + c / b) );
 638     if (!aIsZero) zDiv += (a/b) * log(1. - zMax);
 639     zDiv = min( zMax, max(0., zDiv));
 640     fIntLow = 1. / b;
 641     fIntHigh = 1. - zDiv;
 642     fInt = fIntLow + fIntHigh;
 643   }
 644
 645   // Choice of z, preweighted for peaks at low or high z. (Dummy start values.)
 646   double z = 0.5;
 647   double fPrel = 1.;
 648   double fVal = 1.;
 649   do {
 650     // Choice of z flat good enough for distribution peaked in the middle;
 651     // if not this z can be reused as a random number in general.
 652     z = rndmPtr->flat();
 653     fPrel = 1.;
 654     // When z_max small use flat below z_div and 1/z^c above z_div.
 655     if (peakedNearZero) {
 656       if (fInt * rndmPtr->flat() < fIntLow) z = zDiv * z;
 657       else if (cIsUnity) {z = pow( zDiv, z); fPrel = zDiv / z;}
 658       else { z = pow( zDivC + (1. - zDivC) * z, 1. / (1. - c) );
 659              fPrel = pow( zDiv / z, c); }
 660     // When z_max large use exp( b * (z -z_div) ) below z_div and flat above it.
 661     } else if (peakedNearUnity) {
 662       if (fInt * rndmPtr->flat() < fIntLow) {
 663         z = zDiv + log(z) / b;
 664         fPrel = exp( b * (z - zDiv) );
 665       } else z = zDiv + (1. - zDiv) * z;
 666     }
 667
 668     // Evaluate actual f(z) (if in physical range) and correct.
 669     if (z > 0 && z < 1) {
 670       double fExp = b * (1. / zMax - 1. / z)+ c * log(zMax / z);
 671       if (!aIsZero) fExp += a * log( (1. - z) / (1. - zMax) );
 672       fVal = exp( max( -EXPMAX, min( EXPMAX, fExp) ) ) ;
 673     } else fVal = 0.;
 674   } while (fVal < rndmPtr->flat() * fPrel);
 675
 676   // Done.
 677   return z;
 678
 679 }
 680
 681 //--------------------------------------------------------------------------
 682
 683 // Generate a random z according to the Peterson/SLAC formula
 684 // f(z) = 1 / ( z * (1 - 1/z - epsilon/(1-z))^2 )
 685 //      = z * (1-z)^2 / ((1-z)^2 + epsilon * z)^2.
 686
 687 double StringZ::zPeterson( double epsilon) {
 688
 689   double z, fVal;
 690
 691   // For large epsilon pick z flat and reject,
 692   // knowing that 4 * epsilon * f(z) < 1 everywhere.
 693   if (epsilon > 0.01) {
 694     do {
 695       z = rndmPtr->flat();
 696       fVal = 4. * epsilon * z * pow2(1. - z)
 697         / pow2( pow2(1. - z) + epsilon * z);
 698     } while (fVal < rndmPtr->flat());
 699     return z;
 700   }
 701
 702   // Else split range, using that 4 * epsilon * f(z)
 703   //   < 4 * epsilon / (1 - z)^2 for 0 < z < 1 - 2 * sqrt(epsilon)
 704   //   < 1                       for 1 - 2 * sqrt(epsilon) < z < 1
 705   double epsRoot = sqrt(epsilon);
 706   double epsComb = 0.5 / epsRoot - 1.;
 707   double fIntLow = 4. * epsilon * epsComb;
 708   double fInt = fIntLow + 2. * epsRoot;
 709   do {
 710     if (rndmPtr->flat() * fInt < fIntLow) {
 711       z = 1. - 1. / (1. + rndmPtr->flat() * epsComb);
 712       fVal = z * pow2( pow2(1. - z) / (pow2(1. - z) + epsilon * z) );
 713     } else {
 714       z = 1. - 2. * epsRoot * rndmPtr->flat();
 715       fVal = 4. * epsilon * z * pow2(1. - z)
 716         / pow2( pow2(1. - z) + epsilon * z);
 717     }
 718   } while (fVal < rndmPtr->flat());
 719   return z;
 720
 721 }
 722
 723 //==========================================================================
 724
 725 // The StringPT class.
 726
 727 //--------------------------------------------------------------------------
 728
 729 // Initialize data members of the string pT selection.
 730
 731 void StringPT::init(Settings& settings, Rndm* rndmPtrIn) {
 732
 733   // Save pointer.
 734   rndmPtr        = rndmPtrIn;
 735
 736   // Parameters of the pT width and enhancement.
 737   sigmaQ           = settings.parm("StringPT:sigma") / sqrt(2.);
 738   enhancedFraction = settings.parm("StringPT:enhancedFraction");
 739   enhancedWidth    = settings.parm("StringPT:enhancedWidth");
 740
 741 }
 742
 743 //--------------------------------------------------------------------------
 744
 745 // Generate Gaussian pT such that <p_x^2> = <p_x^2> = sigma^2 = width^2/2,
 746 // but with small fraction multiplied up to a broader spectrum.
 747
 748 pair<double, double> StringPT::pxy() {
 749
 750   double sigma = sigmaQ;
 751   if (rndmPtr->flat() < enhancedFraction) sigma *= enhancedWidth;
 752   pair<double, double> gauss2 = rndmPtr->gauss2();
 753   return pair<double, double>(sigma * gauss2.first, sigma * gauss2.second);
 754
 755 }
 756
 757 //==========================================================================
 758
 759 } // end namespace Pythia8