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