[1] | 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 | |
---|
| 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. 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 | |
---|
| 443 | void 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 | |
---|
| 477 | int 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. |
---|
| 510 | const double StringZ::CFROMUNITY = 0.01; |
---|
| 511 | const double StringZ::AFROMZERO = 0.02; |
---|
| 512 | const double StringZ::AFROMC = 0.01; |
---|
| 513 | |
---|
| 514 | // Do not take exponent of too large or small number. |
---|
| 515 | const double StringZ::EXPMAX = 50.; |
---|
| 516 | |
---|
| 517 | //-------------------------------------------------------------------------- |
---|
| 518 | |
---|
| 519 | // Initialize data members of the string z selection. |
---|
| 520 | |
---|
| 521 | void 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 | |
---|
| 560 | double 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 | |
---|
| 602 | double 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 | |
---|
| 694 | double 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. |
---|
| 740 | const double StringPT::SIGMAMIN = 0.2; |
---|
| 741 | |
---|
| 742 | //-------------------------------------------------------------------------- |
---|
| 743 | |
---|
| 744 | // Initialize data members of the string pT selection. |
---|
| 745 | |
---|
| 746 | void 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 | |
---|
| 767 | pair<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 |
---|