1 // ParticleData.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2010 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 // Function definitions (not found in the header) for the
7 // DecayChannel, ParticleDataEntry and ParticleData classes.
9 #include "ParticleData.h"
10 #include "StandardModel.h"
11 #include "SusyResonanceWidths.h"
13 // Allow string and character manipulation.
18 //==========================================================================
20 // DecayChannel class.
21 // This class holds info on a single decay channel.
23 //--------------------------------------------------------------------------
25 // Check whether id1 occurs anywhere in product list.
27 bool DecayChannel::contains(int id1) const {
30 for (int i = 0; i < nProd; ++i) if (prod[i] == id1) found1 = true;
35 //--------------------------------------------------------------------------
37 // Check whether id1 and id2 occur anywhere in product list.
38 // iF ID1 == ID2 then two copies of this particle must be present.
40 bool DecayChannel::contains(int id1, int id2) const {
44 for (int i = 0; i < nProd; ++i) {
45 if (!found1 && prod[i] == id1) {found1 = true; continue;}
46 if (!found2 && prod[i] == id2) {found2 = true; continue;}
48 return found1 && found2;
52 //--------------------------------------------------------------------------
54 // Check whether id1, id2 and id3 occur anywhere in product list.
55 // iF ID1 == ID2 then two copies of this particle must be present, etc.
57 bool DecayChannel::contains(int id1, int id2, int id3) const {
62 for (int i = 0; i < nProd; ++i) {
63 if (!found1 && prod[i] == id1) {found1 = true; continue;}
64 if (!found2 && prod[i] == id2) {found2 = true; continue;}
65 if (!found3 && prod[i] == id3) {found3 = true; continue;}
67 return found1 && found2 && found3;
71 //==========================================================================
73 // ParticleDataEntry class.
74 // This class holds info on a single particle species.
76 //--------------------------------------------------------------------------
78 // Constants: could be changed here if desired, but normally should not.
79 // These are of technical nature, as described for each.
81 // A particle is invisible if it has neither strong nor electric charge,
82 // and is not made up of constituents that have it. Only relevant for
83 // long-lived particles. This list may need to be extended.
84 const int ParticleDataEntry::INVISIBLENUMBER = 38;
85 const int ParticleDataEntry::INVISIBLETABLE[38] = { 12, 14, 16, 18, 23, 25,
86 32, 33, 35, 36, 39, 41, 1000012, 1000014, 1000016, 1000018, 1000022,
87 1000023, 1000025, 1000035, 1000045, 1000039, 2000012, 2000014, 2000016,
88 2000018, 4900012, 4900014, 4900016, 4900021, 4900022, 4900101, 5000039,
89 5100039, 9900012, 9900014, 9900016, 9900023};
91 // Particles with a read-in tau0 (in mm/c) below this mayDecay by default.
92 const double ParticleDataEntry::MAXTAU0FORDECAY = 1000.;
94 // Particles with a read-in m0 above this isResonance by default.
95 const double ParticleDataEntry::MINMASSRESONANCE = 20.;
97 // Narrow states are assigned nominal mass.
98 const double ParticleDataEntry::NARROWMASS = 1e-6;
100 // Constituent quark masses (d, u, s, c, b).
101 const double ParticleDataEntry::CONSTITUENTMASSTABLE[6]
102 = {0., 0.325, 0.325, 0.50, 1.60, 5.00};
104 //--------------------------------------------------------------------------
106 // Destructor: delete any ResonanceWidths object.
108 ParticleDataEntry::~ParticleDataEntry() {
109 if (resonancePtr != 0) delete resonancePtr;
112 //--------------------------------------------------------------------------
114 // Set initial default values for some quantities.
116 void ParticleDataEntry::setDefaults() {
118 // A particle is a resonance if it is heavy enough.
119 isResonanceSave = (m0Save > MINMASSRESONANCE);
121 // A particle may decay if it is shortlived enough.
122 mayDecaySave = (tau0Save < MAXTAU0FORDECAY);
124 // A particle by default has no external decays.
125 doExternalDecaySave = false;
127 // A particle is invisible if in current table of such.
128 isVisibleSave = true;
129 for (int i = 0; i < INVISIBLENUMBER; ++i)
130 if (idSave == INVISIBLETABLE[i]) isVisibleSave = false;
132 // Normally a resonance should not have width forced to fixed value.
133 doForceWidthSave = false;
135 // Set up constituent masses.
136 setConstituentMass();
138 // No Breit-Wigner mass selection before initialized.
143 //--------------------------------------------------------------------------
145 // Find out if a particle is a hadron.
146 // Only covers normal hadrons, not e.g. R-hadrons.
148 bool ParticleDataEntry::isHadron() const {
150 if (idSave <= 100 || (idSave >= 1000000 && idSave <= 9000000)
151 || idSave >= 9900000) return false;
152 if (idSave == 130 || idSave == 310) return true;
153 if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0)
159 //--------------------------------------------------------------------------
161 // Find out if a particle is a meson.
162 // Only covers normal hadrons, not e.g. R-hadrons.
164 bool ParticleDataEntry::isMeson() const {
166 if (idSave <= 100 || (idSave >= 1000000 && idSave <= 9000000)
167 || idSave >= 9900000) return false;
168 if (idSave == 130 || idSave == 310) return true;
169 if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0
170 || (idSave/1000)%10 != 0) return false;
175 //--------------------------------------------------------------------------
177 // Find out if a particle is a baryon.
178 // Only covers normal hadrons, not e.g. R-hadrons.
180 bool ParticleDataEntry::isBaryon() const {
182 if (idSave <= 1000 || (idSave >= 1000000 && idSave <= 9000000)
183 || idSave >= 9900000) return false;
184 if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0
185 || (idSave/1000)%10 == 0) return false;
191 //--------------------------------------------------------------------------
193 // Extract the heaviest (= largest id) quark in a hadron.
195 int ParticleDataEntry::heaviestQuark(int idIn) const {
197 if (!isHadron()) return 0;
201 if ( (idSave/1000)%10 == 0 ) {
202 hQ = (idSave/100)%10;
203 if (idSave == 130) hQ = 3;
204 if (hQ%2 == 1) hQ = -hQ;
207 } else hQ = (idSave/1000)%10;
210 return (idIn > 0) ? hQ : -hQ;
214 //--------------------------------------------------------------------------
216 // Calculate three times baryon number, i.e. net quark - antiquark number.
218 int ParticleDataEntry::baryonNumberType(int idIn) const {
221 if (isQuark()) return (idIn > 0) ? 1 : -1;
224 if (isDiquark()) return (idIn > 0) ? 2 : -2;
227 if (isBaryon()) return (idIn > 0) ? 3 : -3;
234 //--------------------------------------------------------------------------
236 // Prepare the Breit-Wigner mass selection by precalculating
237 // frequently-used expressions.
239 void ParticleDataEntry::initBWmass() {
241 // Find Breit-Wigner mode for current particle.
242 modeBWnow = particleDataPtr->modeBreitWigner;
243 if ( mWidthSave < NARROWMASS || (mMaxSave > mMinSave
244 && mMaxSave - mMinSave < NARROWMASS) ) modeBWnow = 0;
245 if (modeBWnow == 0) return;
247 // Find atan expressions to be used in random mass selection.
249 atanLow = atan( 2. * (mMinSave - m0Save) / mWidthSave );
250 double atanHigh = (mMaxSave > mMinSave)
251 ? atan( 2. * (mMaxSave - m0Save) / mWidthSave ) : 0.5 * M_PI;
252 atanDif = atanHigh - atanLow;
254 atanLow = atan( (pow2(mMinSave) - pow2(m0Save))
255 / (m0Save * mWidthSave) );
256 double atanHigh = (mMaxSave > mMinSave)
257 ? atan( (pow2(mMaxSave) - pow2(m0Save)) / (m0Save * mWidthSave) )
259 atanDif = atanHigh - atanLow;
262 // Done if no threshold factor.
263 if (modeBWnow%2 == 1) return;
265 // Find average mass threshold for threshold-factor correction.
268 for (int i = 0; i < int(channels.size()); ++i)
269 if (channels[i].onMode() >= 0) {
270 bRatSum += channels[i].bRatio();
271 double mChannelSum = 0.;
272 for (int j = 0; j < channels[i].multiplicity(); ++j)
273 mChannelSum += particleDataPtr->m0( channels[i].product(j) );
274 mThrSum += channels[i].bRatio() * mChannelSum;
276 mThr = (bRatSum == 0.) ? 0. : mThrSum / bRatSum;
278 // Switch off Breit-Wigner if very close to threshold.
279 if (mThr + NARROWMASS > m0Save) {
280 ostringstream osWarn;
281 osWarn << "for id = " << idSave;
282 particleDataPtr->infoPtr->errorMsg("Warning in ParticleDataEntry::"
283 "initBWmass: switching off width", osWarn.str(), true);
289 //--------------------------------------------------------------------------
291 // Function to give mass of a particle, either at the nominal value
292 // or picked according to a (linear or quadratic) Breit-Wigner.
294 double ParticleDataEntry::mass() {
297 if (modeBWnow == 0) return m0Save;
300 // Mass according to a Breit-Wigner linear in m.
301 if (modeBWnow == 1) {
302 mNow = m0Save + 0.5 * mWidthSave
303 * tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
305 // Ditto, but make Gamma proportional to sqrt(m^2 - m_threshold^2).
306 } else if (modeBWnow == 2) {
307 double mWidthNow, fixBW, runBW;
308 double m0ThrS = m0Save*m0Save - mThr*mThr;
310 mNow = m0Save + 0.5 * mWidthSave
311 * tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
312 mWidthNow = mWidthSave * sqrtpos( (mNow*mNow - mThr*mThr) / m0ThrS );
313 fixBW = mWidthSave / (pow2(mNow - m0Save) + pow2(0.5 * mWidthSave));
314 runBW = mWidthNow / (pow2(mNow - m0Save) + pow2(0.5 * mWidthNow));
315 } while (runBW < particleDataPtr->rndmPtr->flat()
316 * particleDataPtr->maxEnhanceBW * fixBW);
318 // Mass according to a Breit-Wigner quadratic in m.
319 } else if (modeBWnow == 3) {
320 m2Now = m0Save*m0Save + m0Save * mWidthSave
321 * tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
322 mNow = sqrtpos( m2Now);
324 // Ditto, but m_0 Gamma_0 -> m Gamma(m) with threshold factor as above.
326 double mwNow, fixBW, runBW;
327 double m2Ref = m0Save * m0Save;
328 double mwRef = m0Save * mWidthSave;
329 double m2Thr = mThr * mThr;
331 m2Now = m2Ref + mwRef * tan( atanLow + atanDif
332 * particleDataPtr->rndmPtr->flat() );
333 mNow = sqrtpos( m2Now);
334 mwNow = mNow * mWidthSave
335 * sqrtpos( (m2Now - m2Thr) / (m2Ref - m2Thr) );
336 fixBW = mwRef / (pow2(m2Now - m2Ref) + pow2(mwRef));
337 runBW = mwNow / (pow2(m2Now - m2Ref) + pow2(mwNow));
338 } while (runBW < particleDataPtr->rndmPtr->flat()
339 * particleDataPtr->maxEnhanceBW * fixBW);
346 //--------------------------------------------------------------------------
348 // Function to calculate running mass at given mass scale.
350 double ParticleDataEntry::mRun(double mHat) {
352 // Except for six quarks return nominal mass.
353 if (idSave > 6) return m0Save;
354 double mQRun = particleDataPtr->mQRun[idSave];
355 double Lam5 = particleDataPtr->Lambda5Run;
357 // For d, u, s quarks start running at 2 GeV (RPP 2006 p. 505).
358 if (idSave < 4) return mQRun * pow ( log(2. / Lam5)
359 / log(max(2., mHat) / Lam5), 12./23.);
361 // For c, b and t quarks start running at respective mass.
362 return mQRun * pow ( log(mQRun / Lam5)
363 / log(max(mQRun, mHat) / Lam5), 12./23.);
367 //--------------------------------------------------------------------------
369 // Rescale all branching ratios to assure normalization to unity.
371 void ParticleDataEntry::rescaleBR(double newSumBR) {
373 // Sum up branching ratios. Find rescaling factor. Rescale.
374 double oldSumBR = 0.;
375 for ( int i = 0; i < int(channels.size()); ++ i)
376 oldSumBR += channels[i].bRatio();
377 double rescaleFactor = newSumBR / oldSumBR;
378 for ( int i = 0; i < int(channels.size()); ++ i)
379 channels[i].rescaleBR(rescaleFactor);
383 //--------------------------------------------------------------------------
385 // Prepare to pick a decay channel.
387 bool ParticleDataEntry::preparePick(int idSgn, double mHat, int idInFlav) {
389 // Reset sum of allowed widths/branching ratios.
392 // For resonances the widths are calculated dynamically.
393 if (resonancePtr != 0) {
394 resonancePtr->widthStore(idSgn, mHat, idInFlav);
395 for (int i = 0; i < int(channels.size()); ++i)
396 currentBRSum += channels[i].currentBR();
398 // Else use normal fixed branching ratios.
402 for (int i = 0; i < int(channels.size()); ++i) {
403 onMode = channels[i].onMode();
405 if ( idSgn > 0 && (onMode == 1 || onMode == 2) )
406 currentBRNow = channels[i].bRatio();
407 else if ( idSgn < 0 && (onMode == 1 || onMode == 3) )
408 currentBRNow = channels[i].bRatio();
409 channels[i].currentBR(currentBRNow);
410 currentBRSum += currentBRNow;
414 // Failure if no channels found with positive branching ratios.
415 return (currentBRSum > 0.);
419 //--------------------------------------------------------------------------
421 // Pick a decay channel according to branching ratios from preparePick.
423 DecayChannel& ParticleDataEntry::pickChannel() {
425 // Find channel in table.
426 int size = channels.size();
427 double rndmBR = currentBRSum * particleDataPtr->rndmPtr->flat();
429 do rndmBR -= channels[++i].currentBR();
430 while (rndmBR > 0. && i < size);
432 // Emergency if no channel found. Done.
433 if (i == size) i = 0;
438 //--------------------------------------------------------------------------
440 // Access methods stored in ResonanceWidths. Could have been
441 // inline in .h, except for problems with forward declarations.
443 void ParticleDataEntry::setResonancePtr(
444 ResonanceWidths* resonancePtrIn) {
445 if (resonancePtr == resonancePtrIn) return;
446 if (resonancePtr != 0) delete resonancePtr;
447 resonancePtr = resonancePtrIn;
450 void ParticleDataEntry::resInit(Info* infoPtrIn, Settings* settingsPtrIn,
451 ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
452 if (resonancePtr != 0) resonancePtr->init(infoPtrIn, settingsPtrIn,
453 particleDataPtrIn, couplingsPtrIn);
456 double ParticleDataEntry::resWidth(int idSgn, double mHat, int idIn,
457 bool openOnly, bool setBR) {
458 return (resonancePtr != 0) ? resonancePtr->width( idSgn, mHat,
459 idIn, openOnly, setBR) : 0.;
462 double ParticleDataEntry::resWidthOpen(int idSgn, double mHat, int idIn) {
463 return (resonancePtr != 0) ? resonancePtr->widthOpen( idSgn, mHat, idIn)
467 double ParticleDataEntry::resWidthStore(int idSgn, double mHat, int idIn) {
468 return (resonancePtr != 0) ? resonancePtr->widthStore( idSgn, mHat, idIn)
472 double ParticleDataEntry::resOpenFrac(int idSgn) {
473 return (resonancePtr != 0) ? resonancePtr->openFrac(idSgn) : 1.;
476 double ParticleDataEntry::resWidthRescaleFactor() {
477 return (resonancePtr != 0) ? resonancePtr->widthRescaleFactor() : 1.;
480 double ParticleDataEntry::resWidthChan(double mHat, int idAbs1,
482 return (resonancePtr != 0) ? resonancePtr->widthChan( mHat, idAbs1,
486 //--------------------------------------------------------------------------
488 // Constituent masses for (d, u, s, c, b) quarks and diquarks.
489 // Hardcoded in CONSTITUENTMASSTABLE so that they are not overwritten
490 // by mistake, and separated from the "normal" masses.
491 // Called both by setDefaults and setM0 so kept as separate method.
493 void ParticleDataEntry::setConstituentMass() {
495 // Equate with the normal masses as default guess.
496 constituentMassSave = m0Save;
498 // Quark masses trivial.
499 if (idSave < 6) constituentMassSave = CONSTITUENTMASSTABLE[idSave];
501 // Diquarks as simple sum of constituent quarks.
502 if (idSave > 1000 && idSave < 10000 && (idSave/10)%10 == 0) {
503 int id1 = idSave/1000;
504 int id2 = (idSave/100)%10;
505 if (id1 <6 && id2 < 6) constituentMassSave
506 = CONSTITUENTMASSTABLE[id1] + CONSTITUENTMASSTABLE[id2];
511 //--------------------------------------------------------------------------
513 // Convert string to lowercase for case-insensitive comparisons.
515 string ParticleDataEntry::toLower(const string& nameConv) {
517 string temp(nameConv);
518 for (int i = 0; i < int(temp.length()); ++i)
519 temp[i] = std::tolower(temp[i]);
524 //==========================================================================
526 // ParticleData class.
527 // This class holds a map of all ParticleDataEntries,
528 // each entry containing info on a particle species.
530 //--------------------------------------------------------------------------
532 // Get data to be distributed among particles during setp.
534 void ParticleData::initCommon() {
536 // Mass generation: fixed mass or linear/quadratic Breit-Wigner.
537 modeBreitWigner = settingsPtr->mode("ParticleData:modeBreitWigner");
539 // Maximum tail enhancement when adding threshold factor to Breit-Wigner.
540 maxEnhanceBW = settingsPtr->parm("ParticleData:maxEnhanceBW");
542 // Find initial MSbar masses for five light flavours.
543 mQRun[1] = settingsPtr->parm("ParticleData:mdRun");
544 mQRun[2] = settingsPtr->parm("ParticleData:muRun");
545 mQRun[3] = settingsPtr->parm("ParticleData:msRun");
546 mQRun[4] = settingsPtr->parm("ParticleData:mcRun");
547 mQRun[5] = settingsPtr->parm("ParticleData:mbRun");
548 mQRun[6] = settingsPtr->parm("ParticleData:mtRun");
550 // Find Lambda5 value to use in running of MSbar masses.
551 double alphaSvalue = settingsPtr->parm("ParticleData:alphaSvalueMRun");
553 alphaS.init( alphaSvalue, 1);
554 Lambda5Run = alphaS.Lambda5();
558 //--------------------------------------------------------------------------
560 // Initialize pointer for particles to the full database, the Breit-Wigners
561 // of normal hadrons and the ResonanceWidths of resonances. For the latter
562 // the order of initialization is essential to get secondary widths right.
564 void ParticleData::initWidths( vector<ResonanceWidths*> resonancePtrs) {
566 // Pointer to database and Breit-Wigner mass initialization for each
568 ResonanceWidths* resonancePtr = 0;
569 for (map<int, ParticleDataEntry>::iterator pdtEntry = pdt.begin();
570 pdtEntry != pdt.end(); ++pdtEntry) {
571 ParticleDataEntry& pdtNow = pdtEntry->second;
572 pdtNow.initPtr( this);
575 // Remove any existing resonances.
576 resonancePtr = pdtNow.getResonancePtr();
577 if (resonancePtr != 0) pdtNow.setResonancePtr(0);
580 // Begin set up new resonance objects.
581 // Z0, W+- and top are almost always needed.
582 resonancePtr = new ResonanceGmZ(23);
583 setResonancePtr( 23, resonancePtr);
584 resonancePtr = new ResonanceW(24);
585 setResonancePtr( 24, resonancePtr);
586 resonancePtr = new ResonanceTop(6);
587 setResonancePtr( 6, resonancePtr);
590 if (!settingsPtr->flag("Higgs:useBSM")) {
591 resonancePtr = new ResonanceH(0, 25);
592 setResonancePtr( 25, resonancePtr);
596 resonancePtr = new ResonanceH(1, 25);
597 setResonancePtr( 25, resonancePtr);
598 resonancePtr = new ResonanceH(2, 35);
599 setResonancePtr( 35, resonancePtr);
600 resonancePtr = new ResonanceH(3, 36);
601 setResonancePtr( 36, resonancePtr);
602 resonancePtr = new ResonanceHchg(37);
603 setResonancePtr( 37, resonancePtr);
606 // A fourth generation: b', t', tau', nu'_tau.
607 resonancePtr = new ResonanceFour(7);
608 setResonancePtr( 7, resonancePtr);
609 resonancePtr = new ResonanceFour(8);
610 setResonancePtr( 8, resonancePtr);
611 resonancePtr = new ResonanceFour(17);
612 setResonancePtr( 17, resonancePtr);
613 resonancePtr = new ResonanceFour(18);
614 setResonancePtr( 18, resonancePtr);
616 // New gauge bosons: Z', W', R.
617 resonancePtr = new ResonanceZprime(32);
618 setResonancePtr( 32, resonancePtr);
619 resonancePtr = new ResonanceWprime(34);
620 setResonancePtr( 34, resonancePtr);
621 resonancePtr = new ResonanceRhorizontal(41);
622 setResonancePtr( 41, resonancePtr);
625 resonancePtr = new ResonanceLeptoquark(42);
626 setResonancePtr( 42, resonancePtr);
630 for(int i = 1; i < 7; i++){
631 resonancePtr = new ResonanceSquark(1000000 + i);
632 setResonancePtr( 1000000 + i, resonancePtr);
633 resonancePtr = new ResonanceSquark(2000000 + i);
634 setResonancePtr( 2000000 + i, resonancePtr);
638 resonancePtr = new ResonanceGluino(1000021);
639 setResonancePtr( 1000021, resonancePtr);
642 resonancePtr = new ResonanceChar(1000024);
643 setResonancePtr( 1000024, resonancePtr);
644 resonancePtr = new ResonanceChar(1000037);
645 setResonancePtr( 1000037, resonancePtr);
648 resonancePtr = new ResonanceNeut(1000022);
649 setResonancePtr( 1000022, resonancePtr);
650 resonancePtr = new ResonanceNeut(1000023);
651 setResonancePtr( 1000023, resonancePtr);
652 resonancePtr = new ResonanceNeut(1000025);
653 setResonancePtr( 1000025, resonancePtr);
654 resonancePtr = new ResonanceNeut(1000035);
655 setResonancePtr( 1000035, resonancePtr);
656 resonancePtr = new ResonanceNeut(1000045);
657 setResonancePtr( 1000045, resonancePtr);
659 // Excited quarks and leptons.
660 for (int i = 1; i < 7; ++i) {
661 resonancePtr = new ResonanceExcited(4000000 + i);
662 setResonancePtr( 4000000 + i, resonancePtr);
664 for (int i = 11; i < 17; ++i) {
665 resonancePtr = new ResonanceExcited(4000000 + i);
666 setResonancePtr( 4000000 + i, resonancePtr);
669 // An excited graviton/gluon in extra-dimensional scenarios.
670 resonancePtr = new ResonanceGraviton(5100039);
671 setResonancePtr( 5100039, resonancePtr);
672 resonancePtr = new ResonanceKKgluon(5100021);
673 setResonancePtr( 5100021, resonancePtr);
675 // A left-right-symmetric scenario with new righthanded neutrinos,
676 // righthanded gauge bosons and doubly charged Higgses.
677 resonancePtr = new ResonanceNuRight(9900012);
678 setResonancePtr( 9900012, resonancePtr);
679 resonancePtr = new ResonanceNuRight(9900014);
680 setResonancePtr( 9900014, resonancePtr);
681 resonancePtr = new ResonanceNuRight(9900016);
682 setResonancePtr( 9900016, resonancePtr);
683 resonancePtr = new ResonanceZRight(9900023);
684 setResonancePtr( 9900023, resonancePtr);
685 resonancePtr = new ResonanceWRight(9900024);
686 setResonancePtr( 9900024, resonancePtr);
687 resonancePtr = new ResonanceHchgchgLeft(9900041);
688 setResonancePtr( 9900041, resonancePtr);
689 resonancePtr = new ResonanceHchgchgRight(9900042);
690 setResonancePtr( 9900042, resonancePtr);
692 // Attach user-defined external resonances and do basic initialization.
693 for (int i = 0; i < int(resonancePtrs.size()); ++i) {
694 int idNow = resonancePtrs[i]->id();
695 resonancePtrs[i]->initBasic(idNow);
696 setResonancePtr( idNow, resonancePtrs[i]);
699 // Set up lists to order resonances in ascending mass.
700 vector<int> idOrdered;
701 vector<double> m0Ordered;
703 // Put Z0 and W+- first, since known to be SM and often off-shell.
704 idOrdered.push_back(23);
705 m0Ordered.push_back(m0(23));
706 idOrdered.push_back(24);
707 m0Ordered.push_back(m0(24));
709 // Loop through particle table to find resonances.
710 for (map<int, ParticleDataEntry>::iterator pdtEntry = pdt.begin();
711 pdtEntry != pdt.end(); ++pdtEntry) {
712 ParticleDataEntry& pdtNow = pdtEntry->second;
713 int idNow = pdtNow.id();
715 // Set up a simple default object for uninitialized resonances.
716 if (pdtNow.isResonance() && pdtNow.getResonancePtr() == 0) {
717 resonancePtr = new ResonanceGeneric(idNow);
718 setResonancePtr( idNow, resonancePtr);
721 // Insert resonances in ascending mass, to respect decay hierarchies.
722 if (pdtNow.getResonancePtr() != 0 && idNow != 23 && idNow != 24) {
723 double m0Now = pdtNow.m0();
724 idOrdered.push_back(idNow);
725 m0Ordered.push_back(m0Now);
726 for (int i = int(idOrdered.size()) - 2; i > 1; --i) {
727 if (m0Ordered[i] < m0Now) break;
728 swap( idOrdered[i], idOrdered[i + 1]);
729 swap( m0Ordered[i], m0Ordered[i + 1]);
734 // Initialize the resonances in ascending mass order.
735 for (int i = 0; i < int(idOrdered.size()); ++i) resInit( idOrdered[i]);
739 //--------------------------------------------------------------------------
741 // Read in database from specific XML file (which may refer to others).
743 bool ParticleData::readXML(string inFile, bool reset) {
745 // Normally reset whole database before beginning.
746 if (reset) {pdt.clear(); isInit = false;}
748 // List of files to be checked.
749 vector<string> files;
750 files.push_back(inFile);
752 // Loop over files. Open them for read.
753 for (int i = 0; i < int(files.size()); ++i) {
754 const char* cstring = files[i].c_str();
755 ifstream is(cstring);
757 // Check that instream is OK.
759 infoPtr->errorMsg("Error in ParticleData::readXML:"
760 " did not find file", files[i]);
764 // Read in one line at a time.
767 while ( getline(is, line) ) {
769 // Get first word of a line.
770 istringstream getfirst(line);
774 // Check for occurence of a particle. Add any continuation lines.
775 if (word1 == "<particle") {
776 while (line.find(">") == string::npos) {
778 getline(is, addLine);
782 // Read in particle properties.
783 int idTmp = intAttributeValue( line, "id");
784 string nameTmp = attributeValue( line, "name");
785 string antiNameTmp = attributeValue( line, "antiName");
786 if (antiNameTmp == "") antiNameTmp = "void";
787 int spinTypeTmp = intAttributeValue( line, "spinType");
788 int chargeTypeTmp = intAttributeValue( line, "chargeType");
789 int colTypeTmp = intAttributeValue( line, "colType");
790 double m0Tmp = doubleAttributeValue( line, "m0");
791 double mWidthTmp = doubleAttributeValue( line, "mWidth");
792 double mMinTmp = doubleAttributeValue( line, "mMin");
793 double mMaxTmp = doubleAttributeValue( line, "mMax");
794 double tau0Tmp = doubleAttributeValue( line, "tau0");
796 // Erase if particle already exists.
797 if (isParticle(idTmp)) pdt.erase(idTmp);
799 // Store new particle. Save pointer, to be used for decay channels.
800 addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
801 colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
802 particlePtr = particleDataEntryPtr(idTmp);
804 // Check for occurence of a decay channel. Add any continuation lines.
805 } else if (word1 == "<channel") {
806 while (line.find(">") == string::npos) {
808 getline(is, addLine);
812 // Read in channel properties - products so far only as a string.
813 int onMode = intAttributeValue( line, "onMode");
814 double bRatio = doubleAttributeValue( line, "bRatio");
815 int meMode = intAttributeValue( line, "meMode");
816 string products = attributeValue( line, "products");
818 // Read in decay products from stream. Must have at least one.
819 istringstream prodStream(products);
820 int prod0 = 0; int prod1 = 0; int prod2 = 0; int prod3 = 0;
821 int prod4 = 0; int prod5 = 0; int prod6 = 0; int prod7 = 0;
822 prodStream >> prod0 >> prod1 >> prod2 >> prod3 >> prod4 >> prod5
825 infoPtr->errorMsg("Error in ParticleData::readXML:"
826 " incomplete decay channel", line);
830 // Store new channel (if particle already known).
831 if (particlePtr == 0) {
832 infoPtr->errorMsg("Error in ParticleData::readXML:"
833 " orphan decay channel", line);
836 particlePtr->addChannel(onMode, bRatio, meMode, prod0, prod1,
837 prod2, prod3, prod4, prod5, prod6, prod7);
839 // Check for occurence of a file also to be read.
840 } else if (word1 == "<file") {
841 string file = attributeValue(line, "name");
843 infoPtr->errorMsg("Warning in ParticleData::readXML:"
844 " skip unrecognized file name", line);
845 } else files.push_back(file);
848 // End of loop over lines in input file and loop over files.
852 // All particle data at this stage defines baseline original.
853 if (reset) for (map<int, ParticleDataEntry>::iterator pdtEntry
854 = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
855 particlePtr = &pdtEntry->second;
856 particlePtr->setHasChanged(false);
865 //--------------------------------------------------------------------------
867 // Print out complete database in numerical order as an XML file.
869 void ParticleData::listXML(string outFile) {
871 // Convert file name to ofstream.
872 const char* cstring = outFile.c_str();
873 ofstream os(cstring);
875 // Iterate through the particle data table.
876 for (map<int, ParticleDataEntry>::iterator pdtEntry
877 = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
878 particlePtr = &pdtEntry->second;
880 // Print particle properties.
881 os << "<particle id=\"" << particlePtr->id() << "\""
882 << " name=\"" << particlePtr->name() << "\"";
883 if (particlePtr->hasAnti())
884 os << " antiName=\"" << particlePtr->name(-1) << "\"";
885 os << " spinType=\"" << particlePtr->spinType() << "\""
886 << " chargeType=\"" << particlePtr->chargeType() << "\""
887 << " colType=\"" << particlePtr->colType() << "\"\n";
888 // Pick format for mass and width based on mass value.
889 double m0Now = particlePtr->m0();
890 if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
891 os << fixed << setprecision(5);
892 else os << scientific << setprecision(3);
893 os << " m0=\"" << m0Now << "\"";
894 if (particlePtr->mWidth() > 0.)
895 os << " mWidth=\"" << particlePtr->mWidth() << "\""
896 << " mMin=\"" << particlePtr->mMin() << "\""
897 << " mMax=\"" << particlePtr->mMax() << "\"";
898 if (particlePtr->tau0() > 0.) os << scientific << setprecision(5)
899 << " tau0=\"" << particlePtr->tau0() << "\"";
902 // Loop through the decay channel table for each particle.
903 if (particlePtr->sizeChannels() > 0) {
904 for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
905 const DecayChannel& channel = particlePtr->channel(i);
906 int mult = channel.multiplicity();
908 // Print decay channel properties.
909 os << " <channel onMode=\"" << channel.onMode() << "\""
910 << fixed << setprecision(7)
911 << " bRatio=\"" << channel.bRatio() << "\"";
912 if (channel.meMode() > 0)
913 os << " meMode=\"" << channel.meMode() << "\"";
914 os << " products=\"";
915 for (int j = 0; j < mult; ++j) {
916 os << channel.product(j);
917 if (j < mult - 1) os << " ";
920 // Finish off line and loop over allowed decay channels.
925 // Finish off existing particle.
926 os << "</particle>\n\n";
932 //--------------------------------------------------------------------------
934 // Read in database from specific free format file.
936 bool ParticleData::readFF(string inFile, bool reset) {
938 // Normally reset whole database before beginning.
939 if (reset) {pdt.clear(); isInit = false;}
941 // Open file for read and check that instream is OK.
942 const char* cstring = inFile.c_str();
943 ifstream is(cstring);
945 infoPtr->errorMsg("Error in ParticleData::readFF:"
946 " did not find file", inFile);
950 // Read in one line at a time.
953 bool readParticle = false;
954 while ( getline(is, line) ) {
956 // Empty lines begins new particle.
957 if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) {
962 // Prepare to use standard read from line.
963 istringstream readLine(line);
965 // Read in a line with particle information.
968 // Properties to be read.
970 string nameTmp, antiNameTmp;
971 int spinTypeTmp, chargeTypeTmp, colTypeTmp;
972 double m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp;
976 readLine >> idTmp >> nameTmp >> antiNameTmp >> spinTypeTmp
977 >> chargeTypeTmp >> colTypeTmp >> m0Tmp >> mWidthTmp
978 >> mMinTmp >> mMaxTmp >> tau0Tmp;
980 // Error printout if something went wrong.
982 infoPtr->errorMsg("Error in ParticleData::readFF:"
983 " incomplete particle", line);
987 // Erase if particle already exists.
988 if (isParticle(idTmp)) pdt.erase(idTmp);
990 // Store new particle. Save pointer, to be used for decay channels.
991 addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
992 colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
993 particlePtr = particleDataEntryPtr(idTmp);
994 readParticle = false;
996 // Read in a line with decay channel information.
999 // Properties to be read.
1003 int prod0 = 0; int prod1 = 0; int prod2 = 0; int prod3 = 0;
1004 int prod4 = 0; int prod5 = 0; int prod6 = 0; int prod7 = 0;
1006 // Read in data from stream. Need at least one decay product.
1007 readLine >> onMode >> bRatio >> meMode >> prod0;
1009 infoPtr->errorMsg("Error in ParticleData::readFF:"
1010 " incomplete decay channel", line);
1013 readLine >> prod1 >> prod2 >> prod3 >> prod4 >> prod5
1016 // Store new channel.
1017 if (particlePtr == 0) {
1018 infoPtr->errorMsg("Error in ParticleData::readFF:"
1019 " orphan decay channel", line);
1022 particlePtr->addChannel(onMode, bRatio, meMode, prod0, prod1,
1023 prod2, prod3, prod4, prod5, prod6, prod7);
1027 // End of while loop over lines in the file.
1037 //--------------------------------------------------------------------------
1039 // Print out complete database in numerical order as a free format file.
1041 void ParticleData::listFF(string outFile) {
1043 // Convert file name to ofstream.
1044 const char* cstring = outFile.c_str();
1045 ofstream os(cstring);
1047 // Iterate through the particle data table.
1048 for (map<int, ParticleDataEntry>::iterator pdtEntry
1049 = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1050 particlePtr = &pdtEntry->second;
1052 // Pick format for mass and width based on mass value.
1053 double m0Now = particlePtr->m0();
1054 if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1055 os << fixed << setprecision(5);
1056 else os << scientific << setprecision(3);
1058 // Print particle properties.
1059 os << "\n" << setw(8) << particlePtr->id() << " "
1060 << left << setw(16) << particlePtr->name() << " "
1061 << setw(16) << particlePtr->name(-1) << " "
1062 << right << setw(2) << particlePtr->spinType() << " "
1063 << setw(2) << particlePtr->chargeType() << " "
1064 << setw(2) << particlePtr->colType() << " "
1065 << setw(10) << particlePtr->m0() << " "
1066 << setw(10) << particlePtr->mWidth() << " "
1067 << setw(10) << particlePtr->mMin() << " "
1068 << setw(10) << particlePtr->mMax() << " "
1069 << scientific << setprecision(5)
1070 << setw(12) << particlePtr->tau0() << "\n";
1072 // Loop through the decay channel table for each particle.
1073 if (particlePtr->sizeChannels() > 0) {
1074 for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1075 const DecayChannel& channel = particlePtr->channel(i);
1076 os << " " << setw(6) << channel.onMode()
1077 << " " << fixed << setprecision(7) << setw(10)
1078 << channel.bRatio() << " "
1079 << setw(3) << channel.meMode() << " ";
1080 for (int j = 0; j < channel.multiplicity(); ++j)
1081 os << setw(8) << channel.product(j) << " ";
1090 //--------------------------------------------------------------------------
1092 // Read in updates from a character string, like a line of a file.
1093 // Is used by readString (and readFile) in Pythia.
1095 bool ParticleData::readString(string lineIn, bool warn, ostream& os) {
1097 // If empty line then done.
1098 if (lineIn.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return true;
1100 // Take copy that will be modified.
1101 string line = lineIn;
1103 // If first character is not a digit then taken to be a comment.
1104 int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
1105 if (!isdigit(line[firstChar])) return true;
1107 // Replace colons and equal signs by blanks to make parsing simpler.
1108 for ( int j = 0; j < int(line.size()); ++ j)
1109 if (line[j] == ':' || line[j] == '=') line[j] = ' ';
1111 // Get particle id and property name.
1114 istringstream getWord(line);
1115 getWord >> idTmp >> property;
1116 property = toLower(property);
1118 // Check that valid particle.
1119 if ( (!isParticle(idTmp) && property != "all" && property != "new")
1121 if (warn) os << "\n Warning: input particle not found in Particle"
1122 << " Data Table; skip:\n " << lineIn << "\n";
1126 // Identify particle property and read + set its value, case by case.
1127 if (property == "name") {
1130 pdt[idTmp].setName(nameTmp);
1133 if (property == "antiname") {
1135 getWord >> antiNameTmp;
1136 pdt[idTmp].setAntiName(antiNameTmp);
1139 if (property == "names") {
1140 string nameTmp, antiNameTmp;
1141 getWord >> nameTmp >> antiNameTmp;
1142 pdt[idTmp].setNames(nameTmp, antiNameTmp);
1145 if (property == "spintype") {
1147 getWord >> spinTypeTmp;
1148 pdt[idTmp].setSpinType(spinTypeTmp);
1151 if (property == "chargetype") {
1153 getWord >> chargeTypeTmp;
1154 pdt[idTmp].setChargeType(chargeTypeTmp);
1157 if (property == "coltype") {
1159 getWord >> colTypeTmp;
1160 pdt[idTmp].setColType(colTypeTmp);
1163 if (property == "m0") {
1166 pdt[idTmp].setM0(m0Tmp);
1169 if (property == "mwidth") {
1171 getWord >> mWidthTmp;
1172 pdt[idTmp].setMWidth(mWidthTmp);
1175 if (property == "mmin") {
1178 pdt[idTmp].setMMin(mMinTmp);
1181 if (property == "mmax") {
1184 pdt[idTmp].setMMax(mMaxTmp);
1187 if (property == "tau0") {
1190 pdt[idTmp].setTau0(tau0Tmp);
1193 if (property == "isresonance") {
1195 getWord >> isresTmp;
1196 bool isResonanceTmp = boolString(isresTmp);
1197 pdt[idTmp].setIsResonance(isResonanceTmp);
1200 if (property == "maydecay") {
1203 bool mayDecayTmp = boolString(mayTmp);
1204 pdt[idTmp].setMayDecay(mayDecayTmp);
1207 if (property == "doexternaldecay") {
1209 getWord >> extdecTmp;
1210 bool doExternalDecayTmp = boolString(extdecTmp);
1211 pdt[idTmp].setDoExternalDecay(doExternalDecayTmp);
1214 if (property == "isvisible") {
1216 getWord >> isvisTmp;
1217 bool isVisibleTmp = boolString(isvisTmp);
1218 pdt[idTmp].setIsVisible(isVisibleTmp);
1221 if (property == "doforcewidth") {
1223 getWord >> doforceTmp;
1224 bool doForceWidthTmp = boolString(doforceTmp);
1225 pdt[idTmp].setDoForceWidth(doForceWidthTmp);
1229 // Addition or complete replacement of a particle.
1230 if (property == "all" || property == "new") {
1232 // Default values for properties to be read.
1233 string nameTmp = "void";
1234 string antiNameTmp = "void";
1235 int spinTypeTmp = 0;
1236 int chargeTypeTmp = 0;
1239 double mWidthTmp = 0.;
1240 double mMinTmp = 0.;
1241 double mMaxTmp = 0.;
1242 double tau0Tmp = 0.;
1244 // Read in data from stream.
1245 getWord >> nameTmp >> antiNameTmp >> spinTypeTmp >> chargeTypeTmp
1246 >> colTypeTmp >> m0Tmp >> mWidthTmp >> mMinTmp >> mMaxTmp
1249 // To keep existing decay channels, only overwrite particle data.
1250 if (property == "all" && isParticle(idTmp)) {
1251 setAll( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1252 colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1254 // Else start over completely from scratch.
1256 if (isParticle(idTmp)) pdt.erase(idTmp);
1257 addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1258 colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1263 // Set onMode of all decay channels in one go.
1264 if (property == "onmode") {
1267 getWord >> onModeIn;
1268 // For onMode allow the optional possibility of Bool input.
1269 if (isdigit(onModeIn[0])) {
1270 istringstream getOnMode(onModeIn);
1271 getOnMode >> onMode;
1272 } else onMode = (boolString(onModeIn)) ? 1 : 0;
1273 for (int i = 0; i < pdt[idTmp].sizeChannels(); ++i)
1274 pdt[idTmp].channel(i).onMode(onMode);
1278 // Selective search for matching decay channels.
1280 if (property == "offifany" || property == "onifany" ||
1281 property == "onposifany" || property == "onnegifany") matchKind = 1;
1282 if (property == "offifall" || property == "onifall" ||
1283 property == "onposifall" || property == "onnegifall") matchKind = 2;
1284 if (property == "offifmatch" || property == "onifmatch" ||
1285 property == "onposifmatch" || property == "onnegifmatch") matchKind = 3;
1286 if (matchKind > 0) {
1288 if (property == "onifany" || property == "onifall"
1289 || property == "onifmatch") onMode = 1;
1290 if (property == "onposifany" || property == "onposifall"
1291 || property == "onposifmatch") onMode = 2;
1292 if (property == "onnegifany" || property == "onnegifall"
1293 || property == "onnegifmatch") onMode = 3;
1295 // Read in particles to match.
1296 vector<int> idToMatch;
1300 if (!getWord) break;
1301 idToMatch.push_back(abs(idRead));
1303 int nToMatch = idToMatch.size();
1305 // Loop over all decay channels.
1306 for (int i = 0; i < pdt[idTmp].sizeChannels(); ++i) {
1307 int multi = pdt[idTmp].channel(i).multiplicity();
1309 // Look for any match at all.
1310 if (matchKind == 1) {
1311 bool foundMatch = false;
1312 for (int j = 0; j < multi; ++j) {
1313 int idNow = abs(pdt[idTmp].channel(i).product(j));
1314 for (int k = 0; k < nToMatch; ++k)
1315 if (idNow == idToMatch[k]) {foundMatch = true; break;}
1316 if (foundMatch) break;
1318 if (foundMatch) pdt[idTmp].channel(i).onMode(onMode);
1320 // Look for match of all products provided.
1322 int nUnmatched = nToMatch;
1323 if (multi < nToMatch);
1324 else if (multi > nToMatch && matchKind == 3);
1326 vector<int> idUnmatched;
1327 for (int k = 0; k < nToMatch; ++k)
1328 idUnmatched.push_back(idToMatch[k]);
1329 for (int j = 0; j < multi; ++j) {
1330 int idNow = abs(pdt[idTmp].channel(i).product(j));
1331 for (int k = 0; k < nUnmatched; ++k)
1332 if (idNow == idUnmatched[k]) {
1333 idUnmatched[k] = idUnmatched[--nUnmatched];
1336 if (nUnmatched == 0) break;
1339 if (nUnmatched == 0) pdt[idTmp].channel(i).onMode(onMode);
1345 // Rescale all branching ratios by common factor.
1346 if (property == "rescalebr") {
1349 pdt[idTmp].rescaleBR(factor);
1353 // Reset decay table in preparation for new input.
1354 if (property == "onechannel") pdt[idTmp].clearChannels();
1356 // Add or change a decay channel: get channel number and new property.
1357 if (property == "addchannel" || property == "onechannel"
1358 || isdigit(property[0])) {
1360 if (property == "addchannel" || property == "onechannel") {
1361 pdt[idTmp].addChannel();
1362 channel = pdt[idTmp].sizeChannels() - 1;
1365 istringstream getChannel(property);
1366 getChannel >> channel;
1367 getWord >> property;
1368 property = toLower(property);
1371 // Check that channel exists.
1372 if (channel < 0 || channel >= pdt[idTmp].sizeChannels()) return false;
1374 // Find decay channel property and value, case by case.
1375 // At same time also do case where all should be replaced.
1376 if (property == "onmode" || property == "all") {
1379 getWord >> onModeIn;
1380 // For onMode allow the optional possibility of Bool input.
1381 if (isdigit(onModeIn[0])) {
1382 istringstream getOnMode(onModeIn);
1383 getOnMode >> onMode;
1384 } else onMode = (boolString(onModeIn)) ? 1 : 0;
1385 pdt[idTmp].channel(channel).onMode(onMode);
1386 if (property == "onmode") return true;
1388 if (property == "bratio" || property == "all") {
1391 pdt[idTmp].channel(channel).bRatio(bRatio);
1392 if (property == "bratio") return true;
1394 if (property == "memode" || property == "all") {
1397 pdt[idTmp].channel(channel).meMode(meMode);
1398 if (property == "memode") return true;
1401 // Scan for products until end of line.
1402 if (property == "products" || property == "all") {
1404 for (int iProd = 0; iProd < 8; ++iProd) {
1407 if (!getWord) break;
1408 pdt[idTmp].channel(channel).product(iProd, idProd);
1411 for (int iProd = nProd; iProd < 8; ++iProd)
1412 pdt[idTmp].channel(channel).product(iProd, 0);
1413 pdt[idTmp].channel(channel).multiplicity(nProd);
1417 // Rescale an existing branching ratio.
1418 if (property == "rescalebr") {
1421 pdt[idTmp].channel(channel).rescaleBR(factor);
1426 // Return false if failed to recognize property.
1427 if (warn) os << "\n Warning: input property not found in Particle"
1428 << " Data Table; skip:\n " << lineIn << "\n";
1433 //--------------------------------------------------------------------------
1435 // Print out complete or changed table of database in numerical order.
1437 void ParticleData::list(bool changedOnly, bool changedRes, ostream& os) {
1439 // Table header; output for bool as off/on.
1441 os << "\n -------- PYTHIA Particle Data Table (complete) --------"
1442 << "------------------------------------------------------------"
1443 << "--------------\n \n";
1446 os << "\n -------- PYTHIA Particle Data Table (changed only) ----"
1447 << "------------------------------------------------------------"
1448 << "--------------\n \n";
1450 os << " id name antiName spn chg col m0"
1451 << " mWidth mMin mMax tau0 res dec ext "
1452 << "vis wid\n no onMode bRatio meMode products \n";
1454 // Iterate through the particle data table. Option to skip unchanged.
1456 for (map<int, ParticleDataEntry>::iterator pdtEntry
1457 = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1458 particlePtr = &pdtEntry->second;
1459 if ( !changedOnly || particlePtr->hasChanged() ||
1460 ( changedRes && particlePtr->getResonancePtr() != 0 ) ) {
1462 // Pick format for mass and width based on mass value.
1463 double m0Now = particlePtr->m0();
1464 if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1465 os << fixed << setprecision(5);
1466 else os << scientific << setprecision(3);
1468 // Print particle properties.
1470 string antiNameTmp = particlePtr->name(-1);
1471 if (antiNameTmp == "void") antiNameTmp = " ";
1472 os << "\n" << setw(8) << particlePtr->id() << " "
1473 << left << setw(16) << particlePtr->name() << " "
1474 << setw(16) << antiNameTmp << " "
1475 << right << setw(2) << particlePtr->spinType() << " "
1476 << setw(2) << particlePtr->chargeType() << " "
1477 << setw(2) << particlePtr->colType() << " "
1478 << setw(10) << particlePtr->m0() << " "
1479 << setw(10) << particlePtr->mWidth() << " "
1480 << setw(10) << particlePtr->mMin() << " "
1481 << setw(10) << particlePtr->mMax() << " "
1482 << scientific << setprecision(5)
1483 << setw(12) << particlePtr->tau0() << " " << setw(2)
1484 << particlePtr->isResonance() << " " << setw(2)
1485 << (particlePtr->mayDecay() && particlePtr->canDecay())
1486 << " " << setw(2) << particlePtr->doExternalDecay() << " "
1487 << setw(2) << particlePtr->isVisible()<< " "
1488 << setw(2) << particlePtr->doForceWidth() << "\n";
1490 // Loop through the decay channel table for each particle.
1491 if (particlePtr->sizeChannels() > 0) {
1492 for (int i = 0; i < int(particlePtr->sizeChannels()); ++i) {
1493 const DecayChannel& channel = particlePtr->channel(i);
1494 os << " " << setprecision(7)
1496 << setw(6) << channel.onMode()
1497 << fixed<< setw(12) << channel.bRatio()
1498 << setw(5) << channel.meMode() << " ";
1499 for (int j = 0; j < channel.multiplicity(); ++j)
1500 os << setw(8) << channel.product(j) << " ";
1508 // End of loop over database contents.
1509 if (changedOnly && nList == 0) os << "\n no particle data has been "
1510 << "changed from its default value \n";
1511 os << "\n -------- End PYTHIA Particle Data Table -----------------"
1512 << "--------------------------------------------------------------"
1513 << "----------\n" << endl;
1517 //--------------------------------------------------------------------------
1519 // Print out partial table of database in input order.
1521 void ParticleData::list(vector<int> idList, ostream& os) {
1523 // Table header; output for bool as off/on.
1524 os << "\n -------- PYTHIA Particle Data Table (partial) ---------"
1525 << "------------------------------------------------------------"
1526 << "--------------\n \n";
1527 os << " id name antiName spn chg col m0"
1528 << " mWidth mMin mMax tau0 res dec ext "
1529 << "vis wid\n no onMode bRatio meMode products \n";
1531 // Iterate through the given list of input particles.
1532 for (int i = 0; i < int(idList.size()); ++i) {
1533 particlePtr = particleDataEntryPtr(idList[i]);
1535 // Pick format for mass and width based on mass value.
1536 double m0Now = particlePtr->m0();
1537 if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1538 os << fixed << setprecision(5);
1539 else os << scientific << setprecision(3);
1541 // Print particle properties.
1542 string antiNameTmp = particlePtr->name(-1);
1543 if (antiNameTmp == "void") antiNameTmp = " ";
1544 os << "\n" << setw(8) << particlePtr->id() << " "
1545 << left << setw(16) << particlePtr->name() << " "
1546 << setw(16) << antiNameTmp << " "
1547 << right << setw(2) << particlePtr->spinType() << " "
1548 << setw(2) << particlePtr->chargeType() << " "
1549 << setw(2) << particlePtr->colType() << " "
1550 << setw(10) << particlePtr->m0() << " "
1551 << setw(10) << particlePtr->mWidth() << " "
1552 << setw(10) << particlePtr->mMin() << " "
1553 << setw(10) << particlePtr->mMax() << " "
1554 << scientific << setprecision(5)
1555 << setw(12) << particlePtr->tau0() << " " << setw(2)
1556 << particlePtr->isResonance() << " " << setw(2)
1557 << (particlePtr->mayDecay() && particlePtr->canDecay())
1558 << " " << setw(2) << particlePtr->doExternalDecay() << " "
1559 << setw(2) << particlePtr->isVisible() << " "
1560 << setw(2) << particlePtr->doForceWidth() << "\n";
1562 // Loop through the decay channel table for each particle.
1563 if (particlePtr->sizeChannels() > 0) {
1564 for (int j = 0; j < int(particlePtr->sizeChannels()); ++j) {
1565 const DecayChannel& channel = particlePtr->channel(j);
1566 os << " " << setprecision(7)
1568 << setw(6) << channel.onMode()
1569 << fixed<< setw(12) << channel.bRatio()
1570 << setw(5) << channel.meMode() << " ";
1571 for (int k = 0; k < channel.multiplicity(); ++k)
1572 os << setw(8) << channel.product(k) << " ";
1579 // End of loop over database contents.
1580 os << "\n -------- End PYTHIA Particle Data Table -----------------"
1581 << "--------------------------------------------------------------"
1582 << "----------\n" << endl;
1586 //--------------------------------------------------------------------------
1588 // Check that table makes sense: e.g. consistent names and mass ranges,
1589 // that branching ratios sum to unity, that charge is conserved and
1590 // that phase space is open in each channel.
1591 // verbosity = 0: mimimal amount of checks, e.g. that no channels open.
1592 // = 1: further warning if individual channels closed
1593 // (except for resonances).
1594 // = 2: also print branching-ratio-averaged threshold mass.
1595 // = 11, 12: as 1, 2, but include resonances in detailed checks.
1597 void ParticleData::checkTable(int verbosity, ostream& os) {
1600 os << "\n -------- PYTHIA Check of Particle Data Table ------------"
1604 // Loop through all particles.
1605 for (map<int, ParticleDataEntry>::iterator pdtEntry
1606 = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1607 particlePtr = &pdtEntry->second;
1609 // Extract some particle properties. Set some flags.
1610 int idNow = particlePtr->id();
1611 bool hasAntiNow = particlePtr->hasAnti();
1612 int spinTypeNow = particlePtr->spinType();
1613 int chargeTypeNow = particlePtr->chargeType();
1614 int baryonTypeNow = particlePtr->baryonNumberType();
1615 double m0Now = particlePtr->m0();
1616 double mMinNow = particlePtr->mMin();
1617 double mMaxNow = particlePtr->mMax();
1618 double mWidthNow = particlePtr->mWidth();
1619 double tau0Now = particlePtr->tau0();
1620 bool isResonanceNow = particlePtr->isResonance();
1621 bool hasPrinted = false;
1622 bool studyCloser = verbosity > 10 || !isResonanceNow;
1624 // Check that particle name consistent with charge information.
1625 string particleName = particlePtr->name(1);
1626 if (particleName.size() > 16) {
1627 os << " Warning: particle " << idNow << " has name " << particleName
1628 << " of length " << particleName.size() << "\n";
1634 for (int i = 0; i < int(particleName.size()); ++i) {
1635 if (particleName[i] == '+') ++nPos;
1636 if (particleName[i] == '-') ++nNeg;
1638 if ( (nPos > 0 && nNeg > 0) || ( nPos + nNeg > 0
1639 && 3 * (nPos - nNeg) != chargeTypeNow )) {
1640 os << " Warning: particle " << idNow << " has name " << particleName
1641 << " inconsistent with charge type " << chargeTypeNow << "\n";
1646 // Check that antiparticle name consistent with charge information.
1648 particleName = particlePtr->name(-1);
1649 if (particleName.size() > 16) {
1650 os << " Warning: particle " << idNow << " has name " << particleName
1651 << " of length " << particleName.size() << "\n";
1657 for (int i = 0; i < int(particleName.size()); ++i) {
1658 if (particleName[i] == '+') ++nPos;
1659 if (particleName[i] == '-') ++nNeg;
1661 if ( (nPos > 0 && nNeg > 0) || ( nPos + nNeg > 0
1662 && 3 * (nPos - nNeg) != -chargeTypeNow )) {
1663 os << " Warning: particle " << -idNow << " has name "
1664 << particleName << " inconsistent with charge type "
1665 << -chargeTypeNow << "\n";
1671 // Check that mass, mass range and width are consistent.
1672 if (particlePtr->useBreitWigner()) {
1673 if (mMinNow > m0Now) {
1674 os << " Error: particle " << idNow << " has mMin "
1675 << fixed << setprecision(5) << mMinNow
1676 << " larger than m0 " << m0Now << "\n";
1680 if (mMaxNow > mMinNow && mMaxNow < m0Now) {
1681 os << " Error: particle " << idNow << " has mMax "
1682 << fixed << setprecision(5) << mMaxNow
1683 << " smaller than m0 " << m0Now << "\n";
1687 if (mMaxNow > mMinNow && mMaxNow - mMinNow < mWidthNow) {
1688 os << " Warning: particle " << idNow << " has mMax - mMin "
1689 << fixed << setprecision(5) << mMaxNow - mMinNow
1690 << " smaller than mWidth " << mWidthNow << "\n";
1696 // Check that particle does not both have width and lifetime.
1697 if (mWidthNow > 0. && tau0Now > 0.) {
1698 os << " Warning: particle " << idNow << " has both nonvanishing width "
1699 << scientific << setprecision(5) << mWidthNow << " and lifetime "
1705 // Begin study decay channels.
1706 if (particlePtr->sizeChannels() > 0) {
1708 // Loop through all decay channels.
1709 double bRatioSum = 0.;
1710 double bRatioPos = 0.;
1711 double bRatioNeg = 0.;
1712 bool hasPosBR = false;
1713 bool hasNegBR = false;
1714 double threshMass = 0.;
1715 bool openChannel = false;
1716 for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1718 // Extract channel properties.
1719 int onMode = particlePtr->channel(i).onMode();
1720 double bRatio = particlePtr->channel(i).bRatio();
1721 int meMode = particlePtr->channel(i).meMode();
1722 int mult = particlePtr->channel(i).multiplicity();
1724 for (int j = 0; j < 8; ++j)
1725 prod[j] = particlePtr->channel(i).product(j);
1727 // Sum branching ratios. Include off-channels.
1728 if (onMode == 0 || onMode == 1) bRatioSum += bRatio;
1729 else if (onMode == 2) {bRatioPos += bRatio; hasPosBR = true;}
1730 else if (onMode == 3) {bRatioNeg += bRatio; hasNegBR = true;}
1732 // Error printout when unknown decay product code.
1733 for (int j = 0; j < 8; ++j) {
1734 if ( prod[j] != 0 && !isParticle(prod[j]) ) {
1735 os << " Error: unknown decay product for " << idNow
1736 << " -> " << prod[j] << "\n";
1743 // Error printout when no decay products or 0 interspersed.
1745 for (int j = 0; j < 8; ++j)
1746 if (prod[j] != 0) nLast = j + 1;
1747 if (mult == 0 || mult != nLast) {
1748 os << " Error: corrupted decay product list for "
1749 << particlePtr->id() << " -> ";
1750 for (int j = 0; j < 8; ++j) os << prod[j] << " ";
1757 // Check charge conservation and open phase space in decay channel.
1758 int chargeTypeSum = -chargeTypeNow;
1759 int baryonTypeSum = -baryonTypeNow;
1760 double avgFinalMass = 0.;
1761 double minFinalMass = 0.;
1762 bool canHandle = true;
1763 for (int j = 0; j < mult; ++j) {
1764 chargeTypeSum += chargeType( prod[j] );
1765 baryonTypeSum += baryonNumberType( prod[j] );
1766 avgFinalMass += m0( prod[j] );
1767 minFinalMass += m0Min( prod[j] );
1768 if (prod[j] == 81 || prod[j] == 82 || prod[j] == 83)
1771 threshMass += bRatio * avgFinalMass;
1773 // Error printout when charge or baryon number not conserved.
1774 if (chargeTypeSum != 0 && canHandle) {
1775 os << " Error: 3*charge changed by " << chargeTypeSum
1776 << " in " << idNow << " -> ";
1777 for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1783 if ( baryonTypeSum != 0 && canHandle && particlePtr->isHadron() ) {
1784 os << " Error: 3*baryon number changed by " << baryonTypeSum
1785 << " in " << idNow << " -> ";
1786 for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1793 // Begin check that some matrix elements are used correctly.
1794 bool correctME = true;
1796 // Check matrix element mode 0: recommended not into partons.
1797 if (meMode == 0 && !isResonanceNow) {
1798 bool hasPartons = false;
1799 for (int j = 0; j < mult; ++j) {
1800 int idAbs = abs(prod[j]);
1801 if ( idAbs < 10 || idAbs == 21 || idAbs == 81 || idAbs == 82
1802 || idAbs == 83 || (idAbs > 1000 && idAbs < 10000
1803 && (idAbs/10)%10 == 0) ) hasPartons = true;
1805 if (hasPartons) correctME = false;
1808 // Check matrix element mode 1: rho/omega -> pi+ pi- pi0.
1809 bool useME1 = ( mult == 3 && spinTypeNow == 3 && idNow > 100
1810 && idNow < 1000 && particlePtr->channel(i).contains(211, -211, 111) );
1811 if ( meMode == 1 && !useME1 ) correctME = false;
1812 if ( meMode != 1 && useME1 ) correctME = false;
1814 // Check matrix element mode 2: polarization in V -> PS + PS.
1815 bool useME2 = ( mult == 2 && spinTypeNow == 3 && idNow > 100
1816 && idNow < 1000 && spinType(prod[0]) == 1
1817 && spinType(prod[1]) == 1 );
1818 if ( meMode == 2 && !useME2 ) correctME = false;
1819 if ( meMode != 2 && useME2 ) correctME = false;
1821 // Check matrix element mode 11, 12 and 13: Dalitz decay with
1822 // one or more particles in addition to the lepton pair,
1823 // or double Dalitz decay.
1824 bool useME11 = ( mult == 3 && !isResonanceNow
1825 && (prod[1] == 11 || prod[1] == 13 || prod[1] == 15)
1826 && prod[2] == -prod[1] );
1827 bool useME12 = ( mult > 3 && !isResonanceNow
1828 && (prod[mult - 2] == 11 || prod[mult - 2] == 13
1829 || prod[mult - 2] == 15) && prod[mult - 1] == -prod[mult - 2] );
1830 bool useME13 = ( mult == 4 && !isResonanceNow
1831 && (prod[0] == 11 || prod[0] == 13) && prod[1] == -prod[0]
1832 && (prod[2] == 11 || prod[2] == 13) && prod[3] == -prod[2] );
1833 if (useME13) useME12 = false;
1834 if ( meMode == 11 && !useME11 ) correctME = false;
1835 if ( meMode != 11 && useME11 ) correctME = false;
1836 if ( meMode == 12 && !useME12 ) correctME = false;
1837 if ( meMode != 12 && useME12 ) correctME = false;
1838 if ( meMode == 13 && !useME13 ) correctME = false;
1839 if ( meMode != 13 && useME13 ) correctME = false;
1841 // Check matrix element mode 21: tau -> nu_tau hadrons.
1842 bool useME21 = (idNow == 15 && mult > 2 && prod[0] == 16
1843 && abs(prod[1]) > 100);
1844 if ( meMode == 21 && !useME21 ) correctME = false;
1845 if ( meMode != 21 && useME21 ) correctME = false;
1847 // Check matrix element mode 22, but only for semileptonic decay.
1848 // For a -> b c d require types u = 2, ubar = -2, d = 1, dbar = -1.
1849 if ( isLepton(prod[0]) && isLepton(prod[1]) ) {
1850 bool useME22 = false;
1854 if (particlePtr->isLepton()) {
1855 typeA = (idNow > 0) ? 1 + (idNow-1)%2 : -1 - (1-idNow)%2;
1856 } else if (particlePtr->isHadron()) {
1857 int hQ = particlePtr->heaviestQuark();
1858 // Special case: for B_c either bbar or c decays.
1859 if (idNow == 541 && heaviestQuark(prod[2]) == -5) hQ = 4;
1860 typeA = (hQ > 0) ? 1 + (hQ-1)%2 : -1 - (1-hQ)%2;
1862 typeB = (prod[0] > 0) ? 1 + (prod[0]-1)%2 : -1 - (1-prod[0])%2;
1863 typeC = (prod[1] > 0) ? 1 + (prod[1]-1)%2 : -1 - (1-prod[1])%2;
1865 if ( (idNow == 130 || idNow == 310) && typeC * typeA < 0)
1867 if (mult == 3 && idNow == 2112 && prod[2] == 2212)
1869 if (mult == 3 && idNow == 3222 && prod[2] == 3122)
1871 if (mult > 2 && typeC == typeA && typeB * typeC < 0
1872 && (typeB + typeC)%2 != 0) useME22 = true;
1873 if ( meMode == 22 && !useME22 ) correctME = false;
1874 if ( meMode != 22 && useME22 ) correctME = false;
1877 // Check for matrix element mode 31.
1880 for (int j = 0; j < mult; ++j) if (prod[j] == 22) ++nGamma;
1881 if (nGamma != 1) correctME = false;
1884 // Check for unknown mode, or resonance-only mode.
1885 if ( !isResonanceNow && (meMode < 0 || (meMode > 2 && meMode <= 10)
1886 || (meMode > 13 && meMode <= 20) || (meMode > 23 && meMode <= 30)
1887 || (meMode > 31 && meMode <= 41) || meMode == 51 || meMode == 61
1888 || meMode == 71 || (meMode > 80 && meMode <= 90)
1889 || (!particlePtr->isOctetHadron() && meMode > 92) ) )
1892 // Print if incorrect matrix element mode.
1894 os << " Warning: meMode " << meMode << " used for "
1896 for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1902 // Warning printout when no phase space for decay.
1903 if ( studyCloser && verbosity > 0 && canHandle && onMode > 0
1904 && particlePtr->m0Min() - minFinalMass < 0. ) {
1905 if (particlePtr->m0Max() - minFinalMass < 0.)
1906 os << " Error: decay never possible for ";
1907 else os << " Warning: decay sometimes not possible for ";
1908 os << idNow << " -> ";
1909 for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1915 // End loop over decay channels.
1916 if (onMode > 0 && bRatio > 0.) openChannel = true;
1919 // Optional printout of threshold.
1920 if (verbosity%10 > 1 && particlePtr->useBreitWigner()) {
1921 threshMass /= bRatioSum;
1922 os << " Info: particle " << idNow << fixed << setprecision(5)
1923 << " has average mass threshold " << threshMass
1924 << " while mMin is " << mMinNow << "\n";
1928 // Error printout when no acceptable decay channels found.
1929 if (studyCloser && !openChannel) {
1930 os << " Error: no acceptable decay channel found for particle "
1936 // Warning printout when branching ratios do not sum to unity.
1937 if (studyCloser && (!hasAntiNow || (!hasPosBR && !hasNegBR))
1938 && abs(bRatioSum + bRatioPos - 1.) > 1e-8) {
1939 os << " Warning: particle " << idNow << fixed << setprecision(8)
1940 << " has branching ratio sum " << bRatioSum << "\n";
1943 } else if (studyCloser && hasAntiNow
1944 && (abs(bRatioSum + bRatioPos - 1.) > 1e-8
1945 || abs(bRatioSum + bRatioNeg - 1.) > 1e-8)) {
1946 os << " Warning: particle " << idNow << fixed << setprecision(8)
1947 << " has branching ratio sum " << bRatioSum + bRatioPos
1948 << " while its antiparticle has " << bRatioSum + bRatioNeg
1954 // End study of decay channels and loop over particles.
1956 if (hasPrinted) os << "\n";
1959 // Final output. Done.
1960 os << " Total number of errors and warnings is " << nErr << "\n";
1961 os << "\n -------- End PYTHIA Check of Particle Data Table --------"
1962 << "------\n" << endl;
1966 //--------------------------------------------------------------------------
1968 // Return the id of the sequentially next particle stored in table.
1970 int ParticleData::nextId(int idIn) {
1972 // Return 0 for negative or unknown codes. Return first for 0.
1973 if (idIn < 0 || (idIn > 0 && !isParticle(idIn))) return 0;
1974 if (idIn == 0) return pdt.begin()->first;
1976 // Find pointer to current particle and step up. Return 0 if impossible.
1977 map<int, ParticleDataEntry>::const_iterator pdtIn = pdt.find(idIn);
1978 if (pdtIn == pdt.end()) return 0;
1979 return (++pdtIn)->first;
1983 //--------------------------------------------------------------------------
1985 // Fractional width associated with open channels of one or two resonances.
1987 double ParticleData::resOpenFrac(int id1In, int id2In, int id3In) {
1993 if (isParticle(id1In)) answer = pdt[abs(id1In)].resOpenFrac(id1In);
1995 // Possibly second resonance.
1996 if (isParticle(id2In)) answer *= pdt[abs(id2In)].resOpenFrac(id2In);
1998 // Possibly third resonance.
1999 if (isParticle(id3In)) answer *= pdt[abs(id2In)].resOpenFrac(id3In);
2006 //--------------------------------------------------------------------------
2008 // Convert string to lowercase for case-insensitive comparisons.
2010 string ParticleData::toLower(const string& nameConv) {
2012 string temp(nameConv);
2013 for (int i = 0; i < int(temp.length()); ++i)
2014 temp[i] = std::tolower(temp[i]);
2019 //--------------------------------------------------------------------------
2021 // Allow several alternative inputs for true/false.
2023 bool ParticleData::boolString(string tag) {
2025 string tagLow = toLower(tag);
2026 return ( tagLow == "true" || tagLow == "1" || tagLow == "on"
2027 || tagLow == "yes" || tagLow == "ok" );
2031 //--------------------------------------------------------------------------
2033 // Extract XML value string following XML attribute.
2035 string ParticleData::attributeValue(string line, string attribute) {
2037 if (line.find(attribute) == string::npos) return "";
2038 int iBegAttri = line.find(attribute);
2039 int iBegQuote = line.find("\"", iBegAttri + 1);
2040 int iEndQuote = line.find("\"", iBegQuote + 1);
2041 return line.substr(iBegQuote + 1, iEndQuote - iBegQuote - 1);
2045 //--------------------------------------------------------------------------
2047 // Extract XML bool value following XML attribute.
2049 bool ParticleData::boolAttributeValue(string line, string attribute) {
2051 string valString = attributeValue(line, attribute);
2052 if (valString == "") return false;
2053 return boolString(valString);
2057 //--------------------------------------------------------------------------
2059 // Extract XML int value following XML attribute.
2061 int ParticleData::intAttributeValue(string line, string attribute) {
2062 string valString = attributeValue(line, attribute);
2063 if (valString == "") return 0;
2064 istringstream valStream(valString);
2066 valStream >> intVal;
2071 //--------------------------------------------------------------------------
2073 // Extract XML double value following XML attribute.
2075 double ParticleData::doubleAttributeValue(string line, string attribute) {
2076 string valString = attributeValue(line, attribute);
2077 if (valString == "") return 0.;
2078 istringstream valStream(valString);
2080 valStream >> doubleVal;
2085 //==========================================================================
2087 } // end namespace Pythia8