[1] | 1 | // ParticleDecays.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 | // ParticleDecays class. |
---|
| 8 | |
---|
| 9 | #include "ParticleDecays.h" |
---|
| 10 | |
---|
| 11 | namespace Pythia8 { |
---|
| 12 | |
---|
| 13 | //========================================================================== |
---|
| 14 | |
---|
| 15 | // The ParticleDecays 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 | // Number of times one tries to let decay happen (for 2 nested loops). |
---|
| 23 | const int ParticleDecays::NTRYDECAY = 10; |
---|
| 24 | |
---|
| 25 | // Number of times one tries to pick valid hadronic content in decay. |
---|
| 26 | const int ParticleDecays::NTRYPICK = 100; |
---|
| 27 | |
---|
| 28 | // Number of times one tries to pick decay topology. |
---|
| 29 | const int ParticleDecays::NTRYMEWT = 1000; |
---|
| 30 | |
---|
| 31 | // Maximal loop count in Dalitz decay treatment. |
---|
| 32 | const int ParticleDecays::NTRYDALITZ = 1000; |
---|
| 33 | |
---|
| 34 | // Minimal Dalitz pair mass this factor above threshold. |
---|
| 35 | const double ParticleDecays::MSAFEDALITZ = 1.000001; |
---|
| 36 | |
---|
| 37 | // These numbers are hardwired empirical parameters, |
---|
| 38 | // intended to speed up the M-generator. |
---|
| 39 | const double ParticleDecays::WTCORRECTION[11] = { 1., 1., 1., |
---|
| 40 | 2., 5., 15., 60., 250., 1250., 7000., 50000. }; |
---|
| 41 | |
---|
| 42 | //-------------------------------------------------------------------------- |
---|
| 43 | |
---|
| 44 | // Initialize and save pointers. |
---|
| 45 | |
---|
| 46 | void ParticleDecays::init(Info* infoPtrIn, Settings& settings, |
---|
| 47 | ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, |
---|
| 48 | Couplings* couplingsPtrIn, TimeShower* timesDecPtrIn, |
---|
| 49 | StringFlav* flavSelPtrIn, DecayHandler* decayHandlePtrIn, |
---|
| 50 | vector<int> handledParticles) { |
---|
| 51 | |
---|
| 52 | // Save pointers to error messages handling and flavour generation. |
---|
| 53 | infoPtr = infoPtrIn; |
---|
| 54 | particleDataPtr = particleDataPtrIn; |
---|
| 55 | rndmPtr = rndmPtrIn; |
---|
| 56 | couplingsPtr = couplingsPtrIn; |
---|
| 57 | flavSelPtr = flavSelPtrIn; |
---|
| 58 | |
---|
| 59 | // Save pointer to timelike shower, as needed in some few decays. |
---|
| 60 | timesDecPtr = timesDecPtrIn; |
---|
| 61 | |
---|
| 62 | // Save pointer for external handling of some decays. |
---|
| 63 | decayHandlePtr = decayHandlePtrIn; |
---|
| 64 | |
---|
| 65 | // Set which particles should be handled externally. |
---|
| 66 | if (decayHandlePtr != 0) |
---|
| 67 | for (int i = 0; i < int(handledParticles.size()); ++i) |
---|
| 68 | particleDataPtr->doExternalDecay(handledParticles[i], true); |
---|
| 69 | |
---|
| 70 | // Safety margin in mass to avoid troubles. |
---|
| 71 | mSafety = settings.parm("ParticleDecays:mSafety"); |
---|
| 72 | |
---|
| 73 | // Lifetime and vertex rules for determining whether decay allowed. |
---|
| 74 | limitTau0 = settings.flag("ParticleDecays:limitTau0"); |
---|
| 75 | tau0Max = settings.parm("ParticleDecays:tau0Max"); |
---|
| 76 | limitTau = settings.flag("ParticleDecays:limitTau"); |
---|
| 77 | tauMax = settings.parm("ParticleDecays:tauMax"); |
---|
| 78 | limitRadius = settings.flag("ParticleDecays:limitRadius"); |
---|
| 79 | rMax = settings.parm("ParticleDecays:rMax"); |
---|
| 80 | limitCylinder = settings.flag("ParticleDecays:limitCylinder"); |
---|
| 81 | xyMax = settings.parm("ParticleDecays:xyMax"); |
---|
| 82 | zMax = settings.parm("ParticleDecays:zMax"); |
---|
| 83 | limitDecay = limitTau0 || limitTau || limitRadius || limitCylinder; |
---|
| 84 | |
---|
| 85 | // B-Bbar mixing parameters. |
---|
| 86 | mixB = settings.flag("ParticleDecays:mixB"); |
---|
| 87 | xBdMix = settings.parm("ParticleDecays:xBdMix"); |
---|
| 88 | xBsMix = settings.parm("ParticleDecays:xBsMix"); |
---|
| 89 | |
---|
| 90 | // Suppression of extra-hadron momenta in semileptonic decays. |
---|
| 91 | sigmaSoft = settings.parm("ParticleDecays:sigmaSoft"); |
---|
| 92 | |
---|
| 93 | // Selection of multiplicity and colours in "phase space" model. |
---|
| 94 | multIncrease = settings.parm("ParticleDecays:multIncrease"); |
---|
| 95 | multIncreaseWeak = settings.parm("ParticleDecays:multIncreaseWeak"); |
---|
| 96 | multRefMass = settings.parm("ParticleDecays:multRefMass"); |
---|
| 97 | multGoffset = settings.parm("ParticleDecays:multGoffset"); |
---|
| 98 | colRearrange = settings.parm("ParticleDecays:colRearrange"); |
---|
| 99 | |
---|
| 100 | // Minimum energy in system (+ m_q) from StringFragmentation. |
---|
| 101 | stopMass = settings.parm("StringFragmentation:stopMass"); |
---|
| 102 | |
---|
| 103 | // Parameters for Dalitz decay virtual gamma mass spectrum. |
---|
| 104 | sRhoDal = pow2(particleDataPtr->m0(113)); |
---|
| 105 | wRhoDal = pow2(particleDataPtr->mWidth(113)); |
---|
| 106 | |
---|
| 107 | // Allow showers in decays to qqbar/gg/ggg/gammagg. |
---|
| 108 | doFSRinDecays = settings.flag("ParticleDecays:FSRinDecays"); |
---|
| 109 | |
---|
| 110 | // Use standard decays or dedicated tau decay package |
---|
| 111 | sophisticatedTau = settings.mode("ParticleDecays:sophisticatedTau"); |
---|
| 112 | |
---|
| 113 | // Initialize the dedicated tau decay handler. |
---|
| 114 | if (sophisticatedTau) tauDecayer.init(infoPtr, &settings, |
---|
| 115 | particleDataPtr, rndmPtr, couplingsPtr); |
---|
| 116 | |
---|
| 117 | } |
---|
| 118 | |
---|
| 119 | //-------------------------------------------------------------------------- |
---|
| 120 | |
---|
| 121 | // Decay a particle; main method. |
---|
| 122 | |
---|
| 123 | bool ParticleDecays::decay( int iDec, Event& event) { |
---|
| 124 | |
---|
| 125 | // Check whether a decay is allowed, given the upcoming decay vertex. |
---|
| 126 | Particle& decayer = event[iDec]; |
---|
| 127 | hasPartons = false; |
---|
| 128 | keepPartons = false; |
---|
| 129 | if (limitDecay && !checkVertex(decayer)) return true; |
---|
| 130 | |
---|
| 131 | // Do not allow resonance decays (beyond handling capability). |
---|
| 132 | if (decayer.isResonance()) { |
---|
| 133 | infoPtr->errorMsg("Warning in ParticleDecays::decay: " |
---|
| 134 | "resonance left undecayed"); |
---|
| 135 | return true; |
---|
| 136 | } |
---|
| 137 | |
---|
| 138 | // Fill the decaying particle in slot 0 of arrays. |
---|
| 139 | idDec = decayer.id(); |
---|
| 140 | iProd.resize(0); |
---|
| 141 | idProd.resize(0); |
---|
| 142 | mProd.resize(0); |
---|
| 143 | iProd.push_back( iDec ); |
---|
| 144 | idProd.push_back( idDec ); |
---|
| 145 | mProd.push_back( decayer.m() ); |
---|
| 146 | |
---|
| 147 | // Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar. |
---|
| 148 | bool hasOscillated = (abs(idDec) == 511 || abs(idDec) == 531) |
---|
| 149 | ? oscillateB(decayer) : false; |
---|
| 150 | if (hasOscillated) {idDec = - idDec; idProd[0] = idDec;} |
---|
| 151 | |
---|
| 152 | // Particle data for decaying particle. |
---|
| 153 | decDataPtr = &decayer.particleDataEntry(); |
---|
| 154 | |
---|
| 155 | // Optionally send on to external decay program. |
---|
| 156 | bool doneExternally = false; |
---|
| 157 | if (decDataPtr->doExternalDecay()) { |
---|
| 158 | pProd.resize(0); |
---|
| 159 | pProd.push_back(decayer.p()); |
---|
| 160 | doneExternally = decayHandlePtr->decay(idProd, mProd, pProd, |
---|
| 161 | iDec, event); |
---|
| 162 | |
---|
| 163 | // If it worked, then store the decay products in the event record. |
---|
| 164 | if (doneExternally) { |
---|
| 165 | mult = idProd.size() - 1; |
---|
| 166 | int status = (hasOscillated) ? 94 : 93; |
---|
| 167 | for (int i = 1; i <= mult; ++i) { |
---|
| 168 | int iPos = event.append( idProd[i], status, iDec, 0, 0, 0, |
---|
| 169 | 0, 0, pProd[i], mProd[i]); |
---|
| 170 | iProd.push_back( iPos); |
---|
| 171 | } |
---|
| 172 | |
---|
| 173 | // Also mark mother decayed and store daughters. |
---|
| 174 | event[iDec].statusNeg(); |
---|
| 175 | event[iDec].daughters( iProd[1], iProd[mult]); |
---|
| 176 | } |
---|
| 177 | } |
---|
| 178 | |
---|
| 179 | // Check if the particle is tau and let the special tau decayer handle it. |
---|
| 180 | if (decayer.idAbs() == 15 && !doneExternally && sophisticatedTau) { |
---|
| 181 | doneExternally = tauDecayer.decay(iDec, event); |
---|
| 182 | if (doneExternally) return true; |
---|
| 183 | } |
---|
| 184 | |
---|
| 185 | // Now begin normal internal decay treatment. |
---|
| 186 | if (!doneExternally) { |
---|
| 187 | |
---|
| 188 | // Allow up to ten tries to pick a channel. |
---|
| 189 | if (!decDataPtr->preparePick(idDec, decayer.m())) return false; |
---|
| 190 | bool foundChannel = false; |
---|
| 191 | bool hasStored = false; |
---|
| 192 | for (int iTryChannel = 0; iTryChannel < NTRYDECAY; ++iTryChannel) { |
---|
| 193 | |
---|
| 194 | // Remove previous failed channel. |
---|
| 195 | if (hasStored) event.popBack(mult); |
---|
| 196 | hasStored = false; |
---|
| 197 | |
---|
| 198 | // Pick new channel. Read out basics. |
---|
| 199 | DecayChannel& channel = decDataPtr->pickChannel(); |
---|
| 200 | meMode = channel.meMode(); |
---|
| 201 | keepPartons = (meMode > 90 && meMode <= 100); |
---|
| 202 | mult = channel.multiplicity(); |
---|
| 203 | |
---|
| 204 | // Allow up to ten tries for each channel (e.g with different masses). |
---|
| 205 | bool foundMode = false; |
---|
| 206 | for (int iTryMode = 0; iTryMode < NTRYDECAY; ++iTryMode) { |
---|
| 207 | idProd.resize(1); |
---|
| 208 | mProd.resize(1); |
---|
| 209 | scale = 0.; |
---|
| 210 | |
---|
| 211 | // Extract and store the decay products in local arrays. |
---|
| 212 | hasPartons = false; |
---|
| 213 | for (int i = 0; i < mult; ++i) { |
---|
| 214 | int idNow = channel.product(i); |
---|
| 215 | int idAbs = abs(idNow); |
---|
| 216 | if ( idAbs < 10 || idAbs == 21 || idAbs == 81 || idAbs == 82 |
---|
| 217 | || idAbs == 83 || (idAbs > 1000 && idAbs < 10000 |
---|
| 218 | && (idAbs/10)%10 == 0) ) hasPartons = true; |
---|
| 219 | if (idDec < 0 && particleDataPtr->hasAnti(idNow)) idNow = -idNow; |
---|
| 220 | double mNow = particleDataPtr->mass(idNow); |
---|
| 221 | idProd.push_back( idNow); |
---|
| 222 | mProd.push_back( mNow); |
---|
| 223 | } |
---|
| 224 | |
---|
| 225 | // Decays into partons usually translate into hadrons. |
---|
| 226 | if (hasPartons && !keepPartons && !pickHadrons()) continue; |
---|
| 227 | |
---|
| 228 | // Need to set colour flow if explicit decay to partons. |
---|
| 229 | cols.resize(0); |
---|
| 230 | acols.resize(0); |
---|
| 231 | for (int i = 0; i <= mult; ++i) { |
---|
| 232 | cols.push_back(0); |
---|
| 233 | acols.push_back(0); |
---|
| 234 | } |
---|
| 235 | if (hasPartons && keepPartons && !setColours(event)) continue; |
---|
| 236 | |
---|
| 237 | // Check that enough phase space for decay. |
---|
| 238 | if (mult > 1) { |
---|
| 239 | double mDiff = mProd[0]; |
---|
| 240 | for (int i = 1; i <= mult; ++i) mDiff -= mProd[i]; |
---|
| 241 | if (mDiff < mSafety) continue; |
---|
| 242 | } |
---|
| 243 | |
---|
| 244 | // End of inner trial loops. Check if succeeded or not. |
---|
| 245 | foundMode = true; |
---|
| 246 | break; |
---|
| 247 | } |
---|
| 248 | if (!foundMode) continue; |
---|
| 249 | |
---|
| 250 | // Store decay products in the event record. |
---|
| 251 | int status = (hasOscillated) ? 92 : 91; |
---|
| 252 | for (int i = 1; i <= mult; ++i) { |
---|
| 253 | int iPos = event.append( idProd[i], status, iDec, 0, 0, 0, |
---|
| 254 | cols[i], acols[i], Vec4(0., 0., 0., 0.), mProd[i], scale); |
---|
| 255 | iProd.push_back( iPos); |
---|
| 256 | } |
---|
| 257 | hasStored = true; |
---|
| 258 | |
---|
| 259 | // Pick mass of Dalitz decay. Temporarily change multiplicity. |
---|
| 260 | if ( (meMode == 11 || meMode == 12 || meMode == 13) |
---|
| 261 | && !dalitzMass() ) continue; |
---|
| 262 | |
---|
| 263 | // Do a decay, split by multiplicity. |
---|
| 264 | bool decayed = false; |
---|
| 265 | if (mult == 1) decayed = oneBody(event); |
---|
| 266 | else if (mult == 2) decayed = twoBody(event); |
---|
| 267 | else if (mult == 3) decayed = threeBody(event); |
---|
| 268 | else decayed = mGenerator(event); |
---|
| 269 | if (!decayed) continue; |
---|
| 270 | |
---|
| 271 | // Kinematics of gamma* -> l- l+ in Dalitz decay. Restore multiplicity. |
---|
| 272 | if (meMode == 11 || meMode == 12 || meMode == 13) |
---|
| 273 | dalitzKinematics(event); |
---|
| 274 | |
---|
| 275 | // End of outer trial loops. |
---|
| 276 | foundChannel = true; |
---|
| 277 | break; |
---|
| 278 | } |
---|
| 279 | |
---|
| 280 | // If the decay worked, then mark mother decayed and store daughters. |
---|
| 281 | if (foundChannel) { |
---|
| 282 | event[iDec].statusNeg(); |
---|
| 283 | event[iDec].daughters( iProd[1], iProd[mult]); |
---|
| 284 | |
---|
| 285 | // Else remove unused daughters and return failure. |
---|
| 286 | } else { |
---|
| 287 | if (hasStored) event.popBack(mult); |
---|
| 288 | infoPtr->errorMsg("Error in ParticleDecays::decay: " |
---|
| 289 | "failed to find workable decay channel"); |
---|
| 290 | return false; |
---|
| 291 | } |
---|
| 292 | |
---|
| 293 | // Now finished normal internal decay treatment. |
---|
| 294 | } |
---|
| 295 | |
---|
| 296 | // Set decay vertex when this is displaced. |
---|
| 297 | if (event[iDec].hasVertex() || event[iDec].tau() > 0.) { |
---|
| 298 | Vec4 vDec = event[iDec].vDec(); |
---|
| 299 | for (int i = 1; i <= mult; ++i) event[iProd[i]].vProd( vDec ); |
---|
| 300 | } |
---|
| 301 | |
---|
| 302 | // Set lifetime of daughters. |
---|
| 303 | for (int i = 1; i <= mult; ++i) |
---|
| 304 | event[iProd[i]].tau( event[iProd[i]].tau0() * rndmPtr->exp() ); |
---|
| 305 | |
---|
| 306 | // In a decay explicitly to partons then optionally do a shower, |
---|
| 307 | // and always flag that partonic system should be fragmented. |
---|
| 308 | if (hasPartons && keepPartons && doFSRinDecays) |
---|
| 309 | timesDecPtr->shower( iProd[1], iProd.back(), event, mProd[0]); |
---|
| 310 | |
---|
| 311 | // For Hidden Valley particles also allow leptons to shower. |
---|
| 312 | else if (event[iDec].idAbs() > 4900000 && event[iDec].idAbs() < 5000000 |
---|
| 313 | && doFSRinDecays && mult == 2 && event[iProd[1]].isLepton()) { |
---|
| 314 | event[iProd[1]].scale(mProd[0]); |
---|
| 315 | event[iProd[2]].scale(mProd[0]); |
---|
| 316 | timesDecPtr->shower( iProd[1], iProd.back(), event, mProd[0]); |
---|
| 317 | } |
---|
| 318 | |
---|
| 319 | // Done. |
---|
| 320 | return true; |
---|
| 321 | |
---|
| 322 | } |
---|
| 323 | |
---|
| 324 | //-------------------------------------------------------------------------- |
---|
| 325 | |
---|
| 326 | // Check whether a decay is allowed, given the upcoming decay vertex. |
---|
| 327 | |
---|
| 328 | bool ParticleDecays::checkVertex(Particle& decayer) { |
---|
| 329 | |
---|
| 330 | // Check whether any of the conditions are not fulfilled. |
---|
| 331 | if (limitTau0 && decayer.tau0() > tau0Max) return false; |
---|
| 332 | if (limitTau && decayer.tau() > tauMax) return false; |
---|
| 333 | if (limitRadius && pow2(decayer.xDec()) + pow2(decayer.yDec()) |
---|
| 334 | + pow2(decayer.zDec()) > pow2(rMax)) return false; |
---|
| 335 | if (limitCylinder && (pow2(decayer.xDec()) + pow2(decayer.yDec()) |
---|
| 336 | > pow2(xyMax) || abs(decayer.zDec()) > zMax) ) return false; |
---|
| 337 | |
---|
| 338 | // Done. |
---|
| 339 | return true; |
---|
| 340 | |
---|
| 341 | } |
---|
| 342 | |
---|
| 343 | //-------------------------------------------------------------------------- |
---|
| 344 | |
---|
| 345 | // Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar. |
---|
| 346 | |
---|
| 347 | bool ParticleDecays::oscillateB(Particle& decayer) { |
---|
| 348 | |
---|
| 349 | // Extract relevant information and decide. |
---|
| 350 | if (!mixB) return false; |
---|
| 351 | double xBmix = (abs(decayer.id()) == 511) ? xBdMix : xBsMix; |
---|
| 352 | double tau = decayer.tau(); |
---|
| 353 | double tau0 = decayer.tau0(); |
---|
| 354 | double probosc = pow2(sin(0.5 * xBmix * tau / tau0)); |
---|
| 355 | return (probosc > rndmPtr->flat()); |
---|
| 356 | |
---|
| 357 | } |
---|
| 358 | |
---|
| 359 | //-------------------------------------------------------------------------- |
---|
| 360 | |
---|
| 361 | // Do a one-body decay. (Rare; e.g. for K0 -> K0_short.) |
---|
| 362 | |
---|
| 363 | bool ParticleDecays::oneBody(Event& event) { |
---|
| 364 | |
---|
| 365 | // References to the particles involved. |
---|
| 366 | Particle& decayer = event[iProd[0]]; |
---|
| 367 | Particle& prod = event[iProd[1]]; |
---|
| 368 | |
---|
| 369 | // Set momentum and expand mother information. |
---|
| 370 | prod.p( decayer.p() ); |
---|
| 371 | prod.m( decayer.m() ); |
---|
| 372 | prod.mother2( iProd[0] ); |
---|
| 373 | |
---|
| 374 | // Done. |
---|
| 375 | return true; |
---|
| 376 | |
---|
| 377 | } |
---|
| 378 | |
---|
| 379 | //-------------------------------------------------------------------------- |
---|
| 380 | |
---|
| 381 | // Do a two-body decay. |
---|
| 382 | |
---|
| 383 | bool ParticleDecays::twoBody(Event& event) { |
---|
| 384 | |
---|
| 385 | // References to the particles involved. |
---|
| 386 | Particle& decayer = event[iProd[0]]; |
---|
| 387 | Particle& prod1 = event[iProd[1]]; |
---|
| 388 | Particle& prod2 = event[iProd[2]]; |
---|
| 389 | |
---|
| 390 | // Masses. |
---|
| 391 | double m0 = mProd[0]; |
---|
| 392 | double m1 = mProd[1]; |
---|
| 393 | double m2 = mProd[2]; |
---|
| 394 | |
---|
| 395 | // Energies and absolute momentum in the rest frame. |
---|
| 396 | if (m1 + m2 + mSafety > m0) return false; |
---|
| 397 | double e1 = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0; |
---|
| 398 | double e2 = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0; |
---|
| 399 | double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2) |
---|
| 400 | * (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0; |
---|
| 401 | |
---|
| 402 | // When meMode = 2, for V -> PS2 + PS3 (V = vector, pseudoscalar), |
---|
| 403 | // need to check if production is PS0 -> PS1/gamma + V. |
---|
| 404 | int iMother = event[iProd[0]].mother1(); |
---|
| 405 | int idSister = 0; |
---|
| 406 | if (meMode == 2) { |
---|
| 407 | if (iMother <= 0 || iMother >= iProd[0]) meMode = 0; |
---|
| 408 | else { |
---|
| 409 | int iDaughter1 = event[iMother].daughter1(); |
---|
| 410 | int iDaughter2 = event[iMother].daughter2(); |
---|
| 411 | if (iDaughter2 != iDaughter1 + 1) meMode = 0; |
---|
| 412 | else { |
---|
| 413 | int idMother = abs( event[iMother].id() ); |
---|
| 414 | if (idMother <= 100 || idMother%10 !=1 |
---|
| 415 | || (idMother/1000)%10 != 0) meMode = 0; |
---|
| 416 | else { |
---|
| 417 | int iSister = (iProd[0] == iDaughter1) ? iDaughter2 : iDaughter1; |
---|
| 418 | idSister = abs( event[iSister].id() ); |
---|
| 419 | if ( (idSister <= 100 || idSister%10 !=1 |
---|
| 420 | || (idSister/1000)%10 != 0) && idSister != 22) meMode = 0; |
---|
| 421 | } |
---|
| 422 | } |
---|
| 423 | } |
---|
| 424 | } |
---|
| 425 | |
---|
| 426 | // Begin loop over matrix-element corrections. |
---|
| 427 | double wtME, wtMEmax; |
---|
| 428 | int loop = 0; |
---|
| 429 | do { |
---|
| 430 | wtME = 1.; |
---|
| 431 | wtMEmax = 1.; |
---|
| 432 | ++loop; |
---|
| 433 | |
---|
| 434 | // Isotropic angles give three-momentum. |
---|
| 435 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
| 436 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
| 437 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
| 438 | double pX = pAbs * sinTheta * cos(phi); |
---|
| 439 | double pY = pAbs * sinTheta * sin(phi); |
---|
| 440 | double pZ = pAbs * cosTheta; |
---|
| 441 | |
---|
| 442 | // Fill four-momenta and boost them away from mother rest frame. |
---|
| 443 | prod1.p( pX, pY, pZ, e1); |
---|
| 444 | prod2.p( -pX, -pY, -pZ, e2); |
---|
| 445 | prod1.bst( decayer.p(), decayer.m() ); |
---|
| 446 | prod2.bst( decayer.p(), decayer.m() ); |
---|
| 447 | |
---|
| 448 | // Matrix element for PS0 -> PS1 + V1 -> PS1 + PS2 + PS3 of form |
---|
| 449 | // cos**2(theta02) in V1 rest frame, and for PS0 -> gamma + V1 |
---|
| 450 | // -> gamma + PS2 + PS3 of form sin**2(theta02). |
---|
| 451 | if (meMode == 2) { |
---|
| 452 | double p10 = decayer.p() * event[iMother].p(); |
---|
| 453 | double p12 = decayer.p() * prod1.p(); |
---|
| 454 | double p02 = event[iMother].p() * prod1.p(); |
---|
| 455 | double s0 = pow2(event[iMother].m()); |
---|
| 456 | double s1 = pow2(decayer.m()); |
---|
| 457 | double s2 = pow2(prod1.m()); |
---|
| 458 | if (idSister != 22) wtME = pow2(p10 * p12 - s1 * p02); |
---|
| 459 | else wtME = s1 * (2. * p10 * p12 * p02 - s1 * p02*p02 |
---|
| 460 | - s0 * p12*p12 - s2 * p10*p10 + s1 * s0 * s2); |
---|
| 461 | wtME = max( wtME, 1e-6 * s1*s1 * s0 * s2); |
---|
| 462 | wtMEmax = (p10*p10 - s1 * s0) * (p12*p12 - s1 * s2); |
---|
| 463 | } |
---|
| 464 | |
---|
| 465 | // Break out of loop if no sensible ME weight. |
---|
| 466 | if(loop > NTRYMEWT) { |
---|
| 467 | infoPtr->errorMsg("ParticleDecays::twoBody: " |
---|
| 468 | "caught in infinite ME weight loop"); |
---|
| 469 | wtME = abs(wtMEmax); |
---|
| 470 | } |
---|
| 471 | |
---|
| 472 | // If rejected, try again with new invariant masses. |
---|
| 473 | } while ( wtME < rndmPtr->flat() * wtMEmax ); |
---|
| 474 | |
---|
| 475 | // Done. |
---|
| 476 | return true; |
---|
| 477 | |
---|
| 478 | } |
---|
| 479 | |
---|
| 480 | //-------------------------------------------------------------------------- |
---|
| 481 | |
---|
| 482 | // Do a three-body decay (except Dalitz decays). |
---|
| 483 | |
---|
| 484 | bool ParticleDecays::threeBody(Event& event) { |
---|
| 485 | |
---|
| 486 | // References to the particles involved. |
---|
| 487 | Particle& decayer = event[iProd[0]]; |
---|
| 488 | Particle& prod1 = event[iProd[1]]; |
---|
| 489 | Particle& prod2 = event[iProd[2]]; |
---|
| 490 | Particle& prod3 = event[iProd[3]]; |
---|
| 491 | |
---|
| 492 | // Mother and sum daughter masses. Fail if too close. |
---|
| 493 | double m0 = mProd[0]; |
---|
| 494 | double m1 = mProd[1]; |
---|
| 495 | double m2 = mProd[2]; |
---|
| 496 | double m3 = mProd[3]; |
---|
| 497 | double mSum = m1 + m2 + m3; |
---|
| 498 | double mDiff = m0 - mSum; |
---|
| 499 | if (mDiff < mSafety) return false; |
---|
| 500 | |
---|
| 501 | // Kinematical limits for 2+3 mass. Maximum phase-space weight. |
---|
| 502 | double m23Min = m2 + m3; |
---|
| 503 | double m23Max = m0 - m1; |
---|
| 504 | double p1Max = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min) |
---|
| 505 | * (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0; |
---|
| 506 | double p23Max = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3) |
---|
| 507 | * (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max; |
---|
| 508 | double wtPSmax = 0.5 * p1Max * p23Max; |
---|
| 509 | |
---|
| 510 | // Begin loop over matrix-element corrections. |
---|
| 511 | double wtME, wtMEmax, wtPS, m23, p1Abs, p23Abs; |
---|
| 512 | do { |
---|
| 513 | wtME = 1.; |
---|
| 514 | wtMEmax = 1.; |
---|
| 515 | |
---|
| 516 | // Pick an intermediate mass m23 flat in the allowed range. |
---|
| 517 | do { |
---|
| 518 | m23 = m23Min + rndmPtr->flat() * mDiff; |
---|
| 519 | |
---|
| 520 | // Translate into relative momenta and find phase-space weight. |
---|
| 521 | p1Abs = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23) |
---|
| 522 | * (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0; |
---|
| 523 | p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3) |
---|
| 524 | * (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23; |
---|
| 525 | wtPS = p1Abs * p23Abs; |
---|
| 526 | |
---|
| 527 | // If rejected, try again with new invariant masses. |
---|
| 528 | } while ( wtPS < rndmPtr->flat() * wtPSmax ); |
---|
| 529 | |
---|
| 530 | // Set up m23 -> m2 + m3 isotropic in its rest frame. |
---|
| 531 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
| 532 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
| 533 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
| 534 | double pX = p23Abs * sinTheta * cos(phi); |
---|
| 535 | double pY = p23Abs * sinTheta * sin(phi); |
---|
| 536 | double pZ = p23Abs * cosTheta; |
---|
| 537 | double e2 = sqrt( m2*m2 + p23Abs*p23Abs); |
---|
| 538 | double e3 = sqrt( m3*m3 + p23Abs*p23Abs); |
---|
| 539 | prod2.p( pX, pY, pZ, e2); |
---|
| 540 | prod3.p( -pX, -pY, -pZ, e3); |
---|
| 541 | |
---|
| 542 | // Set up m0 -> m1 + m23 isotropic in its rest frame. |
---|
| 543 | cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
| 544 | sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
| 545 | phi = 2. * M_PI * rndmPtr->flat(); |
---|
| 546 | pX = p1Abs * sinTheta * cos(phi); |
---|
| 547 | pY = p1Abs * sinTheta * sin(phi); |
---|
| 548 | pZ = p1Abs * cosTheta; |
---|
| 549 | double e1 = sqrt( m1*m1 + p1Abs*p1Abs); |
---|
| 550 | double e23 = sqrt( m23*m23 + p1Abs*p1Abs); |
---|
| 551 | prod1.p( pX, pY, pZ, e1); |
---|
| 552 | |
---|
| 553 | // Boost 2 + 3 to the 0 rest frame. |
---|
| 554 | Vec4 p23( -pX, -pY, -pZ, e23); |
---|
| 555 | prod2.bst( p23, m23 ); |
---|
| 556 | prod3.bst( p23, m23 ); |
---|
| 557 | |
---|
| 558 | // Matrix-element weight for omega/phi -> pi+ pi- pi0. |
---|
| 559 | if (meMode == 1) { |
---|
| 560 | double p1p2 = prod1.p() * prod2.p(); |
---|
| 561 | double p1p3 = prod1.p() * prod3.p(); |
---|
| 562 | double p2p3 = prod2.p() * prod3.p(); |
---|
| 563 | wtME = pow2(m1 * m2 * m3) - pow2(m1 * p2p3) - pow2(m2 * p1p3) |
---|
| 564 | - pow2(m3 * p1p2) + 2. * p1p2 * p1p3 * p2p3; |
---|
| 565 | wtMEmax = pow3(m0 * m0) / 150.; |
---|
| 566 | |
---|
| 567 | // Effective matrix element for nu spectrum in tau -> nu + hadrons. |
---|
| 568 | } else if (meMode == 21) { |
---|
| 569 | double x1 = 2. * prod1.e() / m0; |
---|
| 570 | wtME = x1 * (3. - 2. * x1); |
---|
| 571 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
| 572 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
| 573 | |
---|
| 574 | // Matrix element for weak decay (only semileptonic for c and b). |
---|
| 575 | } else if ((meMode == 22 || meMode == 23) && prod1.isLepton()) { |
---|
| 576 | wtME = m0 * prod1.e() * (prod2.p() * prod3.p()); |
---|
| 577 | wtMEmax = min( pow4(m0) / 16., m0 * (m0 - m1 - m2) * (m0 - m1 - m3) |
---|
| 578 | * (m0 - m2 - m3) ); |
---|
| 579 | |
---|
| 580 | // Effective matrix element for weak decay to hadrons (B -> D, D -> K). |
---|
| 581 | } else if (meMode == 22 || meMode == 23) { |
---|
| 582 | double x1 = 2. * prod1.pAbs() / m0; |
---|
| 583 | wtME = x1 * (3. - 2. * x1); |
---|
| 584 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
| 585 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
| 586 | |
---|
| 587 | // Effective matrix element for gamma spectrum in B -> gamma + hadrons. |
---|
| 588 | } else if (meMode == 31) { |
---|
| 589 | double x1 = 2. * prod1.e() / m0; |
---|
| 590 | wtME = pow3(x1); |
---|
| 591 | double x1Max = 1. - pow2(mSum / m0); |
---|
| 592 | wtMEmax = pow3(x1Max); |
---|
| 593 | |
---|
| 594 | // Matrix-element weight for "onium" -> g + g + g or gamma + g + g. |
---|
| 595 | } else if (meMode == 92) { |
---|
| 596 | double x1 = 2. * prod1.e() / m0; |
---|
| 597 | double x2 = 2. * prod2.e() / m0; |
---|
| 598 | double x3 = 2. * prod3.e() / m0; |
---|
| 599 | wtME = pow2( (1. - x1) / (x2 * x3) ) + pow2( (1. - x2) / (x1 * x3) ) |
---|
| 600 | + pow2( (1. - x3) / (x1 * x2) ); |
---|
| 601 | wtMEmax = 2.; |
---|
| 602 | // For gamma + g + g require minimum mass for g + g system. |
---|
| 603 | if (prod1.id() == 22 && sqrt(1. - x1) * m0 < 2. * stopMass) wtME = 0.; |
---|
| 604 | if (prod2.id() == 22 && sqrt(1. - x2) * m0 < 2. * stopMass) wtME = 0.; |
---|
| 605 | if (prod3.id() == 22 && sqrt(1. - x3) * m0 < 2. * stopMass) wtME = 0.; |
---|
| 606 | } |
---|
| 607 | |
---|
| 608 | // If rejected, try again with new invariant masses. |
---|
| 609 | } while ( wtME < rndmPtr->flat() * wtMEmax ); |
---|
| 610 | |
---|
| 611 | // Boost 1 + 2 + 3 to the current frame. |
---|
| 612 | prod1.bst( decayer.p(), decayer.m() ); |
---|
| 613 | prod2.bst( decayer.p(), decayer.m() ); |
---|
| 614 | prod3.bst( decayer.p(), decayer.m() ); |
---|
| 615 | |
---|
| 616 | // Done. |
---|
| 617 | return true; |
---|
| 618 | |
---|
| 619 | } |
---|
| 620 | |
---|
| 621 | //-------------------------------------------------------------------------- |
---|
| 622 | |
---|
| 623 | // Do a multibody decay using the M-generator algorithm. |
---|
| 624 | |
---|
| 625 | bool ParticleDecays::mGenerator(Event& event) { |
---|
| 626 | |
---|
| 627 | // Mother and sum daughter masses. Fail if too close or inconsistent. |
---|
| 628 | double m0 = mProd[0]; |
---|
| 629 | double mSum = mProd[1]; |
---|
| 630 | for (int i = 2; i <= mult; ++i) mSum += mProd[i]; |
---|
| 631 | double mDiff = m0 - mSum; |
---|
| 632 | if (mDiff < mSafety) return false; |
---|
| 633 | |
---|
| 634 | // Begin setup of intermediate invariant masses. |
---|
| 635 | mInv.resize(0); |
---|
| 636 | for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]); |
---|
| 637 | |
---|
| 638 | // Calculate the maximum weight in the decay. |
---|
| 639 | double wtPS, wtME, wtMEmax; |
---|
| 640 | double wtPSmax = 1. / WTCORRECTION[mult]; |
---|
| 641 | double mMax = mDiff + mProd[mult]; |
---|
| 642 | double mMin = 0.; |
---|
| 643 | for (int i = mult - 1; i > 0; --i) { |
---|
| 644 | mMax += mProd[i]; |
---|
| 645 | mMin += mProd[i+1]; |
---|
| 646 | double mNow = mProd[i]; |
---|
| 647 | wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow) |
---|
| 648 | * (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax; |
---|
| 649 | } |
---|
| 650 | |
---|
| 651 | // Begin loop over matrix-element corrections. |
---|
| 652 | do { |
---|
| 653 | wtME = 1.; |
---|
| 654 | wtMEmax = 1.; |
---|
| 655 | |
---|
| 656 | // Begin loop to find the set of intermediate invariant masses. |
---|
| 657 | do { |
---|
| 658 | wtPS = 1.; |
---|
| 659 | |
---|
| 660 | // Find and order random numbers in descending order. |
---|
| 661 | rndmOrd.resize(0); |
---|
| 662 | rndmOrd.push_back(1.); |
---|
| 663 | for (int i = 1; i < mult - 1; ++i) { |
---|
| 664 | double rndm = rndmPtr->flat(); |
---|
| 665 | rndmOrd.push_back(rndm); |
---|
| 666 | for (int j = i - 1; j > 0; --j) { |
---|
| 667 | if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] ); |
---|
| 668 | else break; |
---|
| 669 | } |
---|
| 670 | } |
---|
| 671 | rndmOrd.push_back(0.); |
---|
| 672 | |
---|
| 673 | // Translate into intermediate masses and find weight. |
---|
| 674 | for (int i = mult - 1; i > 0; --i) { |
---|
| 675 | mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff; |
---|
| 676 | wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i]) |
---|
| 677 | * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i]) |
---|
| 678 | * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i]; |
---|
| 679 | } |
---|
| 680 | |
---|
| 681 | // If rejected, try again with new invariant masses. |
---|
| 682 | } while ( wtPS < rndmPtr->flat() * wtPSmax ); |
---|
| 683 | |
---|
| 684 | // Perform two-particle decays in the respective rest frame. |
---|
| 685 | pInv.resize(mult + 1); |
---|
| 686 | for (int i = 1; i < mult; ++i) { |
---|
| 687 | double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i]) |
---|
| 688 | * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i]) |
---|
| 689 | * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i]; |
---|
| 690 | |
---|
| 691 | // Isotropic angles give three-momentum. |
---|
| 692 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
| 693 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
| 694 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
| 695 | double pX = pAbs * sinTheta * cos(phi); |
---|
| 696 | double pY = pAbs * sinTheta * sin(phi); |
---|
| 697 | double pZ = pAbs * cosTheta; |
---|
| 698 | |
---|
| 699 | // Calculate energies, fill four-momenta. |
---|
| 700 | double eHad = sqrt( mProd[i]*mProd[i] + pAbs*pAbs); |
---|
| 701 | double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs); |
---|
| 702 | event[iProd[i]].p( pX, pY, pZ, eHad); |
---|
| 703 | pInv[i+1].p( -pX, -pY, -pZ, eInv); |
---|
| 704 | } |
---|
| 705 | |
---|
| 706 | // Boost decay products to the mother rest frame. |
---|
| 707 | event[iProd[mult]].p( pInv[mult] ); |
---|
| 708 | for (int iFrame = mult - 1; iFrame > 1; --iFrame) |
---|
| 709 | for (int i = iFrame; i <= mult; ++i) |
---|
| 710 | event[iProd[i]].bst( pInv[iFrame], mInv[iFrame]); |
---|
| 711 | |
---|
| 712 | // Effective matrix element for nu spectrum in tau -> nu + hadrons. |
---|
| 713 | if (meMode == 21 && event[iProd[1]].isLepton()) { |
---|
| 714 | double x1 = 2. * event[iProd[1]].e() / m0; |
---|
| 715 | wtME = x1 * (3. - 2. * x1); |
---|
| 716 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
| 717 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
| 718 | |
---|
| 719 | // Effective matrix element for weak decay (only semileptonic for c and b). |
---|
| 720 | // Particles 4 onwards should be made softer explicitly? |
---|
| 721 | } else if ((meMode == 22 || meMode == 23) && event[iProd[1]].isLepton()) { |
---|
| 722 | Vec4 pRest = event[iProd[3]].p(); |
---|
| 723 | for (int i = 4; i <= mult; ++i) pRest += event[iProd[i]].p(); |
---|
| 724 | wtME = m0 * event[iProd[1]].e() * (event[iProd[2]].p() * pRest); |
---|
| 725 | for (int i = 4; i <= mult; ++i) wtME |
---|
| 726 | *= exp(- event[iProd[i]].pAbs2() / pow2(sigmaSoft) ); |
---|
| 727 | wtMEmax = pow4(m0) / 16.; |
---|
| 728 | |
---|
| 729 | // Effective matrix element for weak decay to hadrons (B -> D, D -> K). |
---|
| 730 | } else if (meMode == 22 || meMode == 23) { |
---|
| 731 | double x1 = 2. * event[iProd[1]].pAbs() / m0; |
---|
| 732 | wtME = x1 * (3. - 2. * x1); |
---|
| 733 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
| 734 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
| 735 | |
---|
| 736 | // Effective matrix element for gamma spectrum in B -> gamma + hadrons. |
---|
| 737 | } else if (meMode == 31) { |
---|
| 738 | double x1 = 2. * event[iProd[1]].e() / m0; |
---|
| 739 | wtME = pow3(x1); |
---|
| 740 | double x1Max = 1. - pow2(mSum / m0); |
---|
| 741 | wtMEmax = pow3(x1Max); |
---|
| 742 | } |
---|
| 743 | |
---|
| 744 | // If rejected, try again with new invariant masses. |
---|
| 745 | } while ( wtME < rndmPtr->flat() * wtMEmax ); |
---|
| 746 | |
---|
| 747 | // Boost decay products to the current frame. |
---|
| 748 | pInv[1].p( event[iProd[0]].p() ); |
---|
| 749 | for (int i = 1; i <= mult; ++i) event[iProd[i]].bst( pInv[1], mInv[1] ); |
---|
| 750 | |
---|
| 751 | // Done. |
---|
| 752 | return true; |
---|
| 753 | |
---|
| 754 | } |
---|
| 755 | |
---|
| 756 | //-------------------------------------------------------------------------- |
---|
| 757 | |
---|
| 758 | // Select mass of lepton pair in a Dalitz decay. |
---|
| 759 | |
---|
| 760 | bool ParticleDecays::dalitzMass() { |
---|
| 761 | |
---|
| 762 | // Mother and sum daughter masses. |
---|
| 763 | double mSum1 = 0; |
---|
| 764 | for (int i = 1; i <= mult - 2; ++i) mSum1 += mProd[i]; |
---|
| 765 | if (meMode == 13) mSum1 *= MSAFEDALITZ; |
---|
| 766 | double mSum2 = MSAFEDALITZ * (mProd[mult -1] + mProd[mult]); |
---|
| 767 | double mDiff = mProd[0] - mSum1 - mSum2; |
---|
| 768 | |
---|
| 769 | // Fail if too close or inconsistent. |
---|
| 770 | if (mDiff < mSafety) return false; |
---|
| 771 | if (idProd[mult - 1] + idProd[mult] != 0 |
---|
| 772 | || mProd[mult - 1] != mProd[mult]) { |
---|
| 773 | infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:" |
---|
| 774 | " inconsistent flavour/mass assignments"); |
---|
| 775 | return false; |
---|
| 776 | } |
---|
| 777 | if ( meMode == 13 && (idProd[1] + idProd[2] != 0 |
---|
| 778 | || mProd[1] != mProd[2]) ) { |
---|
| 779 | infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:" |
---|
| 780 | " inconsistent flavour/mass assignments"); |
---|
| 781 | return false; |
---|
| 782 | } |
---|
| 783 | |
---|
| 784 | // Case 1: one Dalitz pair. |
---|
| 785 | if (meMode == 11 || meMode == 12) { |
---|
| 786 | |
---|
| 787 | // Kinematical limits for gamma* squared mass. |
---|
| 788 | double sGamMin = pow2(mSum2); |
---|
| 789 | double sGamMax = pow2(mProd[0] - mSum1); |
---|
| 790 | // Select virtual gamma squared mass. Guessed form for meMode == 12. |
---|
| 791 | double sGam, wtGam; |
---|
| 792 | int loop = 0; |
---|
| 793 | do { |
---|
| 794 | if (++loop > NTRYDALITZ) return false; |
---|
| 795 | sGam = sGamMin * pow( sGamMax / sGamMin, rndmPtr->flat() ); |
---|
| 796 | wtGam = (1. + 0.5 * sGamMin / sGam) * sqrt(1. - sGamMin / sGam) |
---|
| 797 | * pow3(1. - sGam / sGamMax) * sRhoDal * (sRhoDal + wRhoDal) |
---|
| 798 | / ( pow2(sGam - sRhoDal) + sRhoDal * wRhoDal ); |
---|
| 799 | } while ( wtGam < rndmPtr->flat() ); |
---|
| 800 | |
---|
| 801 | // Store results in preparation for doing a one-less-body decay. |
---|
| 802 | --mult; |
---|
| 803 | mProd[mult] = sqrt(sGam); |
---|
| 804 | |
---|
| 805 | // Case 2: two Dalitz pairs. |
---|
| 806 | } else { |
---|
| 807 | |
---|
| 808 | // Kinematical limits for 1 + 2 and 3 + 4 gamma* masses. |
---|
| 809 | double s0 = pow2(mProd[0]); |
---|
| 810 | double s12Min = pow2(mSum1); |
---|
| 811 | double s12Max = pow2(mProd[0] - mSum2); |
---|
| 812 | double s34Min = pow2(mSum2); |
---|
| 813 | double s34Max = pow2(mProd[0] - mSum1); |
---|
| 814 | |
---|
| 815 | // Select virtual gamma squared masses. Guessed form for meMode == 13. |
---|
| 816 | double s12, s34, wt12, wt34, wtPAbs, wtAll; |
---|
| 817 | int loop = 0; |
---|
| 818 | do { |
---|
| 819 | if (++loop > NTRYDALITZ) return false; |
---|
| 820 | s12 = s12Min * pow( s12Max / s12Min, rndmPtr->flat() ); |
---|
| 821 | wt12 = (1. + 0.5 * s12Min / s12) * sqrt(1. - s12Min / s12) |
---|
| 822 | * sRhoDal * (sRhoDal + wRhoDal) |
---|
| 823 | / ( pow2(s12 - sRhoDal) + sRhoDal * wRhoDal ); |
---|
| 824 | s34 = s34Min * pow( s34Max / s34Min, rndmPtr->flat() ); |
---|
| 825 | wt34 = (1. + 0.5 * s34Min / s34) * sqrt(1. - s34Min / s34) |
---|
| 826 | * sRhoDal * (sRhoDal + wRhoDal) |
---|
| 827 | / ( pow2(s34 - sRhoDal) + sRhoDal * wRhoDal ); |
---|
| 828 | wtPAbs = sqrtpos( pow2(1. - (s12 + s34)/ s0) |
---|
| 829 | - 4. * s12 * s34 / (s0 * s0) ); |
---|
| 830 | wtAll = wt12 * wt34 * pow3(wtPAbs); |
---|
| 831 | if (wtAll > 1.) infoPtr->errorMsg( |
---|
| 832 | "Error in ParticleDecays::dalitzMass: weight > 1"); |
---|
| 833 | } while (wtAll < rndmPtr->flat()); |
---|
| 834 | |
---|
| 835 | // Store results in preparation for doing a two-body decay. |
---|
| 836 | mult = 2; |
---|
| 837 | mProd[1] = sqrt(s12); |
---|
| 838 | mProd[2] = sqrt(s34); |
---|
| 839 | } |
---|
| 840 | |
---|
| 841 | // Done. |
---|
| 842 | return true; |
---|
| 843 | |
---|
| 844 | } |
---|
| 845 | |
---|
| 846 | //-------------------------------------------------------------------------- |
---|
| 847 | |
---|
| 848 | // Do kinematics of gamma* -> l- l+ in Dalitz decay. |
---|
| 849 | |
---|
| 850 | bool ParticleDecays::dalitzKinematics(Event& event) { |
---|
| 851 | |
---|
| 852 | // Restore multiplicity. |
---|
| 853 | int nDal = (meMode < 13) ? 1 : 2; |
---|
| 854 | mult += nDal; |
---|
| 855 | |
---|
| 856 | // Loop over one or two lepton pairs. |
---|
| 857 | for (int iDal = 0; iDal < nDal; ++iDal) { |
---|
| 858 | |
---|
| 859 | // References to the particles involved. |
---|
| 860 | Particle& decayer = event[iProd[0]]; |
---|
| 861 | Particle& prodA = (iDal == 0) ? event[iProd[mult - 1]] |
---|
| 862 | : event[iProd[1]]; |
---|
| 863 | Particle& prodB = (iDal == 0) ? event[iProd[mult]] |
---|
| 864 | : event[iProd[2]]; |
---|
| 865 | |
---|
| 866 | // Reconstruct required rotations and boosts backwards. |
---|
| 867 | Vec4 pDec = decayer.p(); |
---|
| 868 | int iGam = (meMode < 13) ? mult - 1 : 2 - iDal; |
---|
| 869 | Vec4 pGam = event[iProd[iGam]].p(); |
---|
| 870 | pGam.bstback( pDec, decayer.m() ); |
---|
| 871 | double phiGam = pGam.phi(); |
---|
| 872 | pGam.rot( 0., -phiGam); |
---|
| 873 | double thetaGam = pGam.theta(); |
---|
| 874 | pGam.rot( -thetaGam, 0.); |
---|
| 875 | |
---|
| 876 | // Masses and phase space in gamma* rest frame. |
---|
| 877 | double mGam = (meMode < 13) ? mProd[mult - 1] : mProd[2 - iDal]; |
---|
| 878 | double mA = prodA.m(); |
---|
| 879 | double mB = prodB.m(); |
---|
| 880 | double mGamMin = MSAFEDALITZ * (mA + mB); |
---|
| 881 | double mGamRat = pow2(mGamMin / mGam); |
---|
| 882 | double pGamAbs = 0.5 * sqrtpos( (mGam - mA - mB) * (mGam + mA + mB) ); |
---|
| 883 | |
---|
| 884 | // Set up decay in gamma* rest frame, reference along +z axis. |
---|
| 885 | double cosTheta, cos2Theta; |
---|
| 886 | do { |
---|
| 887 | cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
| 888 | cos2Theta = cosTheta * cosTheta; |
---|
| 889 | } while ( 1. + cos2Theta + mGamRat * (1. - cos2Theta) |
---|
| 890 | < 2. * rndmPtr->flat() ); |
---|
| 891 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
| 892 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
| 893 | double pX = pGamAbs * sinTheta * cos(phi); |
---|
| 894 | double pY = pGamAbs * sinTheta * sin(phi); |
---|
| 895 | double pZ = pGamAbs * cosTheta; |
---|
| 896 | double eA = sqrt( mA*mA + pGamAbs*pGamAbs); |
---|
| 897 | double eB = sqrt( mB*mB + pGamAbs*pGamAbs); |
---|
| 898 | prodA.p( pX, pY, pZ, eA); |
---|
| 899 | prodB.p( -pX, -pY, -pZ, eB); |
---|
| 900 | |
---|
| 901 | // Boost to lab frame. |
---|
| 902 | prodA.bst( pGam, mGam); |
---|
| 903 | prodB.bst( pGam, mGam); |
---|
| 904 | prodA.rot( thetaGam, phiGam); |
---|
| 905 | prodB.rot( thetaGam, phiGam); |
---|
| 906 | prodA.bst( pDec, decayer.m() ); |
---|
| 907 | prodB.bst( pDec, decayer.m() ); |
---|
| 908 | } |
---|
| 909 | |
---|
| 910 | // Done. |
---|
| 911 | return true; |
---|
| 912 | |
---|
| 913 | } |
---|
| 914 | |
---|
| 915 | //-------------------------------------------------------------------------- |
---|
| 916 | |
---|
| 917 | // Translate a partonic content into a set of actual hadrons. |
---|
| 918 | |
---|
| 919 | bool ParticleDecays::pickHadrons() { |
---|
| 920 | |
---|
| 921 | // Find partonic decay products. Rest are known id's, mainly hadrons, |
---|
| 922 | // when necessary shuffled to beginning of idProd list. |
---|
| 923 | idPartons.resize(0); |
---|
| 924 | int nPartons = 0; |
---|
| 925 | int nKnown = 0; |
---|
| 926 | bool closedGLoop = false; |
---|
| 927 | for (int i = 1; i <= mult; ++i) { |
---|
| 928 | int idAbs = abs(idProd[i]); |
---|
| 929 | if ( idAbs < 9 || (idAbs > 1000 && idAbs < 10000 && (idAbs/10)%10 == 0) |
---|
| 930 | || idAbs == 81 || idAbs == 82 || idAbs == 83) { |
---|
| 931 | ++nPartons; |
---|
| 932 | idPartons.push_back(idProd[i]); |
---|
| 933 | if (idAbs == 83) closedGLoop = true; |
---|
| 934 | } else { |
---|
| 935 | ++nKnown; |
---|
| 936 | if (nPartons > 0) { |
---|
| 937 | idProd[nKnown] = idProd[i]; |
---|
| 938 | mProd[nKnown] = mProd[i]; |
---|
| 939 | } |
---|
| 940 | } |
---|
| 941 | } |
---|
| 942 | |
---|
| 943 | // Replace generic spectator flavour code by the actual one. |
---|
| 944 | for (int i = 0; i < nPartons; ++i) { |
---|
| 945 | int idPart = idPartons[i]; |
---|
| 946 | int idNew = idPart; |
---|
| 947 | if (idPart == 81) { |
---|
| 948 | int idAbs = abs(idDec); |
---|
| 949 | if ( (idAbs/1000)%10 == 0 ) { |
---|
| 950 | idNew = -(idAbs/10)%10; |
---|
| 951 | if ((idAbs/100)%2 == 1) idNew = -idNew; |
---|
| 952 | } else if ( (idAbs/100)%10 >= (idAbs/10)%10 ) |
---|
| 953 | idNew = 100 * ((idAbs/10)%100) + 3; |
---|
| 954 | else idNew = 1000 * ((idAbs/10)%10) + 100 * ((idAbs/100)%10) + 1; |
---|
| 955 | if (idDec < 0) idNew = -idNew; |
---|
| 956 | |
---|
| 957 | // Replace generic random flavour by a randomly selected one. |
---|
| 958 | } else if (idPart == 82 || idPart == 83) { |
---|
| 959 | double mFlav; |
---|
| 960 | do { |
---|
| 961 | int idDummy = -flavSelPtr->pickLightQ(); |
---|
| 962 | FlavContainer flavDummy(idDummy, idPart - 82); |
---|
| 963 | do idNew = flavSelPtr->pick(flavDummy).id; |
---|
| 964 | while (idNew == 0); |
---|
| 965 | mFlav = particleDataPtr->constituentMass(idNew); |
---|
| 966 | } while (2. * mFlav + stopMass > mProd[0]); |
---|
| 967 | } else if (idPart == -82 || idPart == -83) { |
---|
| 968 | idNew = -idPartons[i-1]; |
---|
| 969 | } |
---|
| 970 | idPartons[i] = idNew; |
---|
| 971 | } |
---|
| 972 | |
---|
| 973 | // Determine whether fixed multiplicity or to be selected at random. |
---|
| 974 | int nMin = max( 2, nKnown + nPartons / 2); |
---|
| 975 | if (meMode == 23) nMin = 3; |
---|
| 976 | if (meMode > 41 && meMode <= 50) nMin = meMode - 40; |
---|
| 977 | if (meMode > 51 && meMode <= 60) nMin = meMode - 50; |
---|
| 978 | int nFix = 0; |
---|
| 979 | if (meMode == 0) nFix = nMin; |
---|
| 980 | if (meMode == 11) nFix = 3; |
---|
| 981 | if (meMode == 12) nFix = 4; |
---|
| 982 | if (meMode > 61 && meMode <= 70) nFix = meMode - 60; |
---|
| 983 | if (meMode > 71 && meMode <= 80) nFix = meMode - 70; |
---|
| 984 | if (nFix > 0 && nKnown + nPartons/2 > nFix) return false; |
---|
| 985 | |
---|
| 986 | // Initial values for loop to set new hadronic content. |
---|
| 987 | int nFilled, nTotal, nNew, nSpec, nLeft; |
---|
| 988 | double mDiff; |
---|
| 989 | int nTry = 0; |
---|
| 990 | bool diquarkClash = false; |
---|
| 991 | bool usedChannel = false; |
---|
| 992 | |
---|
| 993 | // Begin loop; interrupt if multiple tries fail. |
---|
| 994 | do { |
---|
| 995 | ++nTry; |
---|
| 996 | if (nTry > NTRYPICK) return false; |
---|
| 997 | |
---|
| 998 | // Initialize variables inside new try. |
---|
| 999 | nFilled = nKnown + 1; |
---|
| 1000 | idProd.resize(nFilled); |
---|
| 1001 | mProd.resize(nFilled); |
---|
| 1002 | nTotal = nKnown; |
---|
| 1003 | nSpec = 0; |
---|
| 1004 | nLeft = nPartons; |
---|
| 1005 | mDiff = mProd[0]; |
---|
| 1006 | for (int i = 1; i < nFilled; ++i) mDiff -= mProd[i]; |
---|
| 1007 | diquarkClash = false; |
---|
| 1008 | usedChannel = false; |
---|
| 1009 | |
---|
| 1010 | // For weak decays collapse spectators to one particle. |
---|
| 1011 | if ( (meMode == 22 || meMode == 23) && nLeft > 1) { |
---|
| 1012 | FlavContainer flav1( idPartons[nPartons - 2] ); |
---|
| 1013 | FlavContainer flav2( idPartons[nPartons - 1] ); |
---|
| 1014 | int idHad; |
---|
| 1015 | do idHad = flavSelPtr->combine( flav1, flav2); |
---|
| 1016 | while (idHad == 0); |
---|
| 1017 | double mHad = particleDataPtr->mass(idHad); |
---|
| 1018 | mDiff -= mHad; |
---|
| 1019 | idProd.push_back( idHad); |
---|
| 1020 | mProd.push_back( mHad); |
---|
| 1021 | ++nFilled; |
---|
| 1022 | nSpec = 1; |
---|
| 1023 | nLeft -= 2; |
---|
| 1024 | } |
---|
| 1025 | |
---|
| 1026 | // If there are partons left, then determine how many hadrons to pick. |
---|
| 1027 | if (nLeft > 0) { |
---|
| 1028 | |
---|
| 1029 | // For B -> gamma + X use geometrical distribution. |
---|
| 1030 | if (meMode == 31) { |
---|
| 1031 | double geom = rndmPtr->flat(); |
---|
| 1032 | nTotal = 1; |
---|
| 1033 | do { |
---|
| 1034 | ++nTotal; |
---|
| 1035 | geom *= 2.; |
---|
| 1036 | } while (geom < 1. && nTotal < 10); |
---|
| 1037 | |
---|
| 1038 | // Calculate mass excess and from there average multiplicity. |
---|
| 1039 | } else if (nFix == 0) { |
---|
| 1040 | double multIncreaseNow = (meMode == 23) |
---|
| 1041 | ? multIncreaseWeak : multIncrease; |
---|
| 1042 | double mDiffPS = mDiff; |
---|
| 1043 | for (int i = 0; i < nLeft; ++i) |
---|
| 1044 | mDiffPS -= particleDataPtr->constituentMass( idPartons[i] ); |
---|
| 1045 | double average = 0.5 * (nKnown + nSpec) + 0.25 * nPartons |
---|
| 1046 | + multIncreaseNow * log( max( 1.1, mDiffPS / multRefMass ) ); |
---|
| 1047 | if (closedGLoop) average += multGoffset; |
---|
| 1048 | |
---|
| 1049 | // Pick multiplicity according to Poissonian. |
---|
| 1050 | double value = 1.; |
---|
| 1051 | double sum = 1.; |
---|
| 1052 | for (int nNow = nMin + 1; nNow <= 10; ++nNow) { |
---|
| 1053 | value *= average / nNow; |
---|
| 1054 | sum += value; |
---|
| 1055 | } |
---|
| 1056 | nTotal = nMin; |
---|
| 1057 | value = 1.; |
---|
| 1058 | sum *= rndmPtr->flat(); |
---|
| 1059 | sum -= value; |
---|
| 1060 | if (sum > 0.) do { |
---|
| 1061 | ++nTotal; |
---|
| 1062 | value *= average / nTotal; |
---|
| 1063 | sum -= value; |
---|
| 1064 | } while (sum > 0. && nTotal < 10); |
---|
| 1065 | |
---|
| 1066 | // Alternatively predetermined multiplicity. |
---|
| 1067 | } else { |
---|
| 1068 | nTotal = nFix; |
---|
| 1069 | } |
---|
| 1070 | nNew = nTotal - nKnown - nSpec; |
---|
| 1071 | |
---|
| 1072 | // Set up ends of fragmented system, as copy of idPartons. |
---|
| 1073 | flavEnds.resize(0); |
---|
| 1074 | for (int i = 0; i < nLeft; ++i) { |
---|
| 1075 | flavEnds.push_back( FlavContainer(idPartons[i]) ); |
---|
| 1076 | if (abs(idPartons[i]) > 100) flavSelPtr->assignPopQ( flavEnds[i] ); |
---|
| 1077 | } |
---|
| 1078 | |
---|
| 1079 | // Fragment off at random, but save nLeft/2 for final recombination. |
---|
| 1080 | if (nNew > nLeft/2) { |
---|
| 1081 | FlavContainer flavNew; |
---|
| 1082 | int idHad; |
---|
| 1083 | for (int i = 0; i < nNew - nLeft/2; ++i) { |
---|
| 1084 | // When four quarks consider last one to be spectator. |
---|
| 1085 | int iEnd = int( (nLeft - 1.) * rndmPtr->flat() ); |
---|
| 1086 | // Pick new flavour and form a new hadron. |
---|
| 1087 | do { |
---|
| 1088 | flavNew = flavSelPtr->pick( flavEnds[iEnd] ); |
---|
| 1089 | idHad = flavSelPtr->combine( flavEnds[iEnd], flavNew); |
---|
| 1090 | } while (idHad == 0); |
---|
| 1091 | // Store new hadron and endpoint flavour. |
---|
| 1092 | idProd.push_back( idHad); |
---|
| 1093 | flavEnds[iEnd].anti(flavNew); |
---|
| 1094 | } |
---|
| 1095 | } |
---|
| 1096 | |
---|
| 1097 | // When only two quarks left, combine to form final hadron. |
---|
| 1098 | if (nLeft == 2) { |
---|
| 1099 | int idHad; |
---|
| 1100 | if ( abs(flavEnds[0].id) > 8 && abs(flavEnds[1].id) > 8) |
---|
| 1101 | diquarkClash = true; |
---|
| 1102 | else { |
---|
| 1103 | do idHad = flavSelPtr->combine( flavEnds[0], flavEnds[1]); |
---|
| 1104 | while (idHad == 0); |
---|
| 1105 | idProd.push_back( idHad); |
---|
| 1106 | } |
---|
| 1107 | |
---|
| 1108 | // If four quarks, decide how to pair them up. |
---|
| 1109 | } else { |
---|
| 1110 | int iEnd1 = 0; |
---|
| 1111 | int iEnd2 = 1; |
---|
| 1112 | int iEnd3 = 2; |
---|
| 1113 | int iEnd4 = 3; |
---|
| 1114 | if ( rndmPtr->flat() < colRearrange) iEnd2 = 3; |
---|
| 1115 | int relColSign = |
---|
| 1116 | ( (flavEnds[iEnd1].id > 0 && flavEnds[iEnd1].id < 9) |
---|
| 1117 | || flavEnds[iEnd1].id < -10 ) ? 1 : -1; |
---|
| 1118 | if ( (flavEnds[iEnd2].id < 0 && flavEnds[iEnd2].id > -9) |
---|
| 1119 | || flavEnds[iEnd2].id > 10 ) relColSign *= -1; |
---|
| 1120 | if (relColSign == 1) iEnd2 = 2; |
---|
| 1121 | if (iEnd2 == 2) iEnd3 = 1; |
---|
| 1122 | if (iEnd2 == 3) iEnd4 = 1; |
---|
| 1123 | |
---|
| 1124 | // Then combine to get final two hadrons. |
---|
| 1125 | int idHad; |
---|
| 1126 | if ( abs(flavEnds[iEnd1].id) > 8 && abs(flavEnds[iEnd2].id) > 8) |
---|
| 1127 | diquarkClash = true; |
---|
| 1128 | else { |
---|
| 1129 | do idHad = flavSelPtr->combine( flavEnds[iEnd1], flavEnds[iEnd2]); |
---|
| 1130 | while (idHad == 0); |
---|
| 1131 | idProd.push_back( idHad); |
---|
| 1132 | } |
---|
| 1133 | if ( abs(flavEnds[iEnd3].id) > 8 && abs(flavEnds[iEnd4].id) > 8) |
---|
| 1134 | diquarkClash = true; |
---|
| 1135 | else { |
---|
| 1136 | do idHad = flavSelPtr->combine( flavEnds[iEnd3], flavEnds[iEnd4]); |
---|
| 1137 | while (idHad == 0); |
---|
| 1138 | idProd.push_back( idHad); |
---|
| 1139 | } |
---|
| 1140 | } |
---|
| 1141 | |
---|
| 1142 | // Find masses of the new hadrons. |
---|
| 1143 | for (int i = nFilled; i < int(idProd.size()) ; ++i) { |
---|
| 1144 | double mHad = particleDataPtr->mass(idProd[i]); |
---|
| 1145 | mProd.push_back( mHad); |
---|
| 1146 | mDiff -= mHad; |
---|
| 1147 | } |
---|
| 1148 | } |
---|
| 1149 | |
---|
| 1150 | // Optional: check that this decay mode is not explicitly defined. |
---|
| 1151 | if ( (meMode > 61 && meMode <= 80) && mDiff > mSafety && !diquarkClash ) { |
---|
| 1152 | int idMatch[10], idPNow; |
---|
| 1153 | usedChannel = false; |
---|
| 1154 | bool matched = false; |
---|
| 1155 | // Loop through all channels. Done if not same multiplicity. |
---|
| 1156 | for (int i = 0; i < decDataPtr->sizeChannels(); ++i) { |
---|
| 1157 | DecayChannel& channel = decDataPtr->channel(i); |
---|
| 1158 | if (channel.multiplicity() != nTotal) continue; |
---|
| 1159 | for (int k = 0; k < nTotal; ++k) idMatch[k] = channel.product(k); |
---|
| 1160 | // Match particles one by one until fail. |
---|
| 1161 | // Do not distinguish K0/K0bar/K0short/K0long at this stage. |
---|
| 1162 | for (int j = 0; j < nTotal; ++j) { |
---|
| 1163 | matched = false; |
---|
| 1164 | idPNow = idProd[j + 1]; |
---|
| 1165 | if (idPNow == -311 || idPNow == 130 || idPNow == 310) idPNow = 311; |
---|
| 1166 | for (int k = 0; k < nTotal - j; ++k) |
---|
| 1167 | if (idMatch[k] == idPNow || (idMatch[k] == -311 && idPNow == 311)) { |
---|
| 1168 | // Compress list of unmatched when matching worked. |
---|
| 1169 | idMatch[k] = idMatch[nTotal - 1 - j]; |
---|
| 1170 | matched = true; |
---|
| 1171 | break; |
---|
| 1172 | } |
---|
| 1173 | if (!matched) break; |
---|
| 1174 | } |
---|
| 1175 | // If matching worked, then chosen channel to be rejected. |
---|
| 1176 | if (matched) {usedChannel = true; break;} |
---|
| 1177 | } |
---|
| 1178 | } |
---|
| 1179 | |
---|
| 1180 | // Keep on trying until enough phase space and no clash. |
---|
| 1181 | } while (mDiff < mSafety || diquarkClash || usedChannel); |
---|
| 1182 | |
---|
| 1183 | // Update particle multiplicity. |
---|
| 1184 | mult = idProd.size() - 1; |
---|
| 1185 | |
---|
| 1186 | // For Dalitz decays shuffle Dalitz pair to the end of the list. |
---|
| 1187 | if (meMode == 11 || meMode == 12) { |
---|
| 1188 | int iL1 = 0; |
---|
| 1189 | int iL2 = 0; |
---|
| 1190 | for (int i = 1; i <= mult; ++i) { |
---|
| 1191 | if (idProd[i] == 11 || idProd[i] == 13 || idProd[i] == 15) iL1 = i; |
---|
| 1192 | if (idProd[i] == -11 || idProd[i] == -13 || idProd[i] == -15) iL2 = i; |
---|
| 1193 | } |
---|
| 1194 | if (iL1 > 0 && iL2 > 0) { |
---|
| 1195 | int idL1 = idProd[iL1]; |
---|
| 1196 | int idL2 = idProd[iL2]; |
---|
| 1197 | double mL1 = mProd[iL1]; |
---|
| 1198 | double mL2 = mProd[iL2]; |
---|
| 1199 | int iMove = 0; |
---|
| 1200 | for (int i = 1; i <= mult; ++i) if (i != iL1 && i != iL2) { |
---|
| 1201 | ++iMove; |
---|
| 1202 | idProd[iMove] = idProd[i]; |
---|
| 1203 | mProd[iMove] = mProd[i]; |
---|
| 1204 | } |
---|
| 1205 | idProd[mult - 1] = idL1; |
---|
| 1206 | idProd[mult] = idL2; |
---|
| 1207 | mProd[mult - 1] = mL1; |
---|
| 1208 | mProd[mult] = mL2; |
---|
| 1209 | } |
---|
| 1210 | } |
---|
| 1211 | |
---|
| 1212 | // Done. |
---|
| 1213 | return true; |
---|
| 1214 | |
---|
| 1215 | } |
---|
| 1216 | |
---|
| 1217 | //-------------------------------------------------------------------------- |
---|
| 1218 | |
---|
| 1219 | // Set colour flow and scale in a decay explicitly to partons. |
---|
| 1220 | |
---|
| 1221 | bool ParticleDecays::setColours(Event& event) { |
---|
| 1222 | |
---|
| 1223 | // Decay to q qbar (or qbar q). |
---|
| 1224 | if (meMode == 91 && idProd[1] > 0 && idProd[1] < 9) { |
---|
| 1225 | int newCol = event.nextColTag(); |
---|
| 1226 | cols[1] = newCol; |
---|
| 1227 | acols[2] = newCol; |
---|
| 1228 | } else if (meMode == 91 && idProd[1] < 0 && idProd[1] > -9) { |
---|
| 1229 | int newCol = event.nextColTag(); |
---|
| 1230 | cols[2] = newCol; |
---|
| 1231 | acols[1] = newCol; |
---|
| 1232 | |
---|
| 1233 | // Decay to g g. |
---|
| 1234 | } else if (meMode == 91 && idProd[1] == 21) { |
---|
| 1235 | int newCol1 = event.nextColTag(); |
---|
| 1236 | int newCol2 = event.nextColTag(); |
---|
| 1237 | cols[1] = newCol1; |
---|
| 1238 | acols[1] = newCol2; |
---|
| 1239 | cols[2] = newCol2; |
---|
| 1240 | acols[2] = newCol1; |
---|
| 1241 | |
---|
| 1242 | // Decay to g g g. |
---|
| 1243 | } else if (meMode == 92 && idProd[1] == 21 && idProd[2] == 21 |
---|
| 1244 | && idProd[3] == 21) { |
---|
| 1245 | int newCol1 = event.nextColTag(); |
---|
| 1246 | int newCol2 = event.nextColTag(); |
---|
| 1247 | int newCol3 = event.nextColTag(); |
---|
| 1248 | cols[1] = newCol1; |
---|
| 1249 | acols[1] = newCol2; |
---|
| 1250 | cols[2] = newCol2; |
---|
| 1251 | acols[2] = newCol3; |
---|
| 1252 | cols[3] = newCol3; |
---|
| 1253 | acols[3] = newCol1; |
---|
| 1254 | |
---|
| 1255 | // Decay to g g gamma: locate which is gamma. |
---|
| 1256 | } else if (meMode == 92) { |
---|
| 1257 | int iGlu1 = (idProd[1] == 21) ? 1 : 3; |
---|
| 1258 | int iGlu2 = (idProd[2] == 21) ? 2 : 3; |
---|
| 1259 | int newCol1 = event.nextColTag(); |
---|
| 1260 | int newCol2 = event.nextColTag(); |
---|
| 1261 | cols[iGlu1] = newCol1; |
---|
| 1262 | acols[iGlu1] = newCol2; |
---|
| 1263 | cols[iGlu2] = newCol2; |
---|
| 1264 | acols[iGlu2] = newCol1; |
---|
| 1265 | |
---|
| 1266 | // Unknown decay mode means failure. |
---|
| 1267 | } else return false; |
---|
| 1268 | |
---|
| 1269 | // Set maximum scale to be mass of decaying particle. |
---|
| 1270 | scale = mProd[0]; |
---|
| 1271 | |
---|
| 1272 | // Done. |
---|
| 1273 | return true; |
---|
| 1274 | |
---|
| 1275 | } |
---|
| 1276 | |
---|
| 1277 | //========================================================================== |
---|
| 1278 | |
---|
| 1279 | } // end namespace Pythia8 |
---|
| 1280 | |
---|