1
// ParticleData.cc is a part of the PYTHIA event generator.
2
// Copyright (C) 2012 Torbjorn Sjostrand.
3
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4
// Please respect the MCnet Guidelines, see GUIDELINES for details.
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 = 50;
85
const int ParticleDataEntry::INVISIBLETABLE[50] = { 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, 4900102,
89
4900103, 4900104, 4900105, 4900106, 4900107, 4900108, 4900111, 4900113,
90
4900211, 4900213, 4900991, 5000039, 5100039, 9900012, 9900014, 9900016,
93
// Particles with a read-in tau0 (in mm/c) below this mayDecay by default.
94
const double ParticleDataEntry::MAXTAU0FORDECAY = 1000.;
96
// Particles with a read-in m0 above this isResonance by default.
97
const double ParticleDataEntry::MINMASSRESONANCE = 20.;
99
// Narrow states are assigned nominal mass.
100
const double ParticleDataEntry::NARROWMASS = 1e-6;
102
// Constituent quark masses (d, u, s, c, b, -, -, -, g).
103
const double ParticleDataEntry::CONSTITUENTMASSTABLE[10]
104
= {0., 0.325, 0.325, 0.50, 1.60, 5.00, 0., 0., 0., 0.7};
106
//--------------------------------------------------------------------------
108
// Destructor: delete any ResonanceWidths object.
110
ParticleDataEntry::~ParticleDataEntry() {
111
if (resonancePtr != 0) delete resonancePtr;
114
//--------------------------------------------------------------------------
116
// Set initial default values for some quantities.
118
void ParticleDataEntry::setDefaults() {
120
// A particle is a resonance if it is heavy enough.
121
isResonanceSave = (m0Save > MINMASSRESONANCE);
123
// A particle may decay if it is shortlived enough.
124
mayDecaySave = (tau0Save < MAXTAU0FORDECAY);
126
// A particle by default has no external decays.
127
doExternalDecaySave = false;
129
// A particle is invisible if in current table of such.
130
isVisibleSave = true;
131
for (int i = 0; i < INVISIBLENUMBER; ++i)
132
if (idSave == INVISIBLETABLE[i]) isVisibleSave = false;
134
// Normally a resonance should not have width forced to fixed value.
135
doForceWidthSave = false;
137
// Set up constituent masses.
138
setConstituentMass();
140
// No Breit-Wigner mass selection before initialized.
145
//--------------------------------------------------------------------------
147
// Find out if a particle is a hadron.
148
// Only covers normal hadrons, not e.g. R-hadrons.
150
bool ParticleDataEntry::isHadron() const {
152
if (idSave <= 100 || (idSave >= 1000000 && idSave <= 9000000)
153
|| idSave >= 9900000) return false;
154
if (idSave == 130 || idSave == 310) return true;
155
if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0)
161
//--------------------------------------------------------------------------
163
// Find out if a particle is a meson.
164
// Only covers normal hadrons, not e.g. R-hadrons.
166
bool ParticleDataEntry::isMeson() const {
168
if (idSave <= 100 || (idSave >= 1000000 && idSave <= 9000000)
169
|| idSave >= 9900000) return false;
170
if (idSave == 130 || idSave == 310) return true;
171
if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0
172
|| (idSave/1000)%10 != 0) return false;
177
//--------------------------------------------------------------------------
179
// Find out if a particle is a baryon.
180
// Only covers normal hadrons, not e.g. R-hadrons.
182
bool ParticleDataEntry::isBaryon() const {
184
if (idSave <= 1000 || (idSave >= 1000000 && idSave <= 9000000)
185
|| idSave >= 9900000) return false;
186
if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0
187
|| (idSave/1000)%10 == 0) return false;
193
//--------------------------------------------------------------------------
195
// Extract the heaviest (= largest id) quark in a hadron.
197
int ParticleDataEntry::heaviestQuark(int idIn) const {
199
if (!isHadron()) return 0;
203
if ( (idSave/1000)%10 == 0 ) {
204
hQ = (idSave/100)%10;
205
if (idSave == 130) hQ = 3;
206
if (hQ%2 == 1) hQ = -hQ;
209
} else hQ = (idSave/1000)%10;
212
return (idIn > 0) ? hQ : -hQ;
216
//--------------------------------------------------------------------------
218
// Calculate three times baryon number, i.e. net quark - antiquark number.
220
int ParticleDataEntry::baryonNumberType(int idIn) const {
223
if (isQuark()) return (idIn > 0) ? 1 : -1;
226
if (isDiquark()) return (idIn > 0) ? 2 : -2;
229
if (isBaryon()) return (idIn > 0) ? 3 : -3;
236
//--------------------------------------------------------------------------
238
// Prepare the Breit-Wigner mass selection by precalculating
239
// frequently-used expressions.
241
void ParticleDataEntry::initBWmass() {
243
// Find Breit-Wigner mode for current particle.
244
modeBWnow = particleDataPtr->modeBreitWigner;
245
if ( mWidthSave < NARROWMASS || (mMaxSave > mMinSave
246
&& mMaxSave - mMinSave < NARROWMASS) ) modeBWnow = 0;
247
if (modeBWnow == 0) return;
249
// Find atan expressions to be used in random mass selection.
251
atanLow = atan( 2. * (mMinSave - m0Save) / mWidthSave );
252
double atanHigh = (mMaxSave > mMinSave)
253
? atan( 2. * (mMaxSave - m0Save) / mWidthSave ) : 0.5 * M_PI;
254
atanDif = atanHigh - atanLow;
256
atanLow = atan( (pow2(mMinSave) - pow2(m0Save))
257
/ (m0Save * mWidthSave) );
258
double atanHigh = (mMaxSave > mMinSave)
259
? atan( (pow2(mMaxSave) - pow2(m0Save)) / (m0Save * mWidthSave) )
261
atanDif = atanHigh - atanLow;
264
// Done if no threshold factor.
265
if (modeBWnow%2 == 1) return;
267
// Find average mass threshold for threshold-factor correction.
270
for (int i = 0; i < int(channels.size()); ++i)
271
if (channels[i].onMode() >= 0) {
272
bRatSum += channels[i].bRatio();
273
double mChannelSum = 0.;
274
for (int j = 0; j < channels[i].multiplicity(); ++j)
275
mChannelSum += particleDataPtr->m0( channels[i].product(j) );
276
mThrSum += channels[i].bRatio() * mChannelSum;
278
mThr = (bRatSum == 0.) ? 0. : mThrSum / bRatSum;
280
// Switch off Breit-Wigner if very close to threshold.
281
if (mThr + NARROWMASS > m0Save) {
282
ostringstream osWarn;
283
osWarn << "for id = " << idSave;
284
particleDataPtr->infoPtr->errorMsg("Warning in ParticleDataEntry::"
285
"initBWmass: switching off width", osWarn.str(), true);
291
//--------------------------------------------------------------------------
293
// Function to give mass of a particle, either at the nominal value
294
// or picked according to a (linear or quadratic) Breit-Wigner.
296
double ParticleDataEntry::mass() {
299
if (modeBWnow == 0) return m0Save;
302
// Mass according to a Breit-Wigner linear in m.
303
if (modeBWnow == 1) {
304
mNow = m0Save + 0.5 * mWidthSave
305
* tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
307
// Ditto, but make Gamma proportional to sqrt(m^2 - m_threshold^2).
308
} else if (modeBWnow == 2) {
309
double mWidthNow, fixBW, runBW;
310
double m0ThrS = m0Save*m0Save - mThr*mThr;
312
mNow = m0Save + 0.5 * mWidthSave
313
* tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
314
mWidthNow = mWidthSave * sqrtpos( (mNow*mNow - mThr*mThr) / m0ThrS );
315
fixBW = mWidthSave / (pow2(mNow - m0Save) + pow2(0.5 * mWidthSave));
316
runBW = mWidthNow / (pow2(mNow - m0Save) + pow2(0.5 * mWidthNow));
317
} while (runBW < particleDataPtr->rndmPtr->flat()
318
* particleDataPtr->maxEnhanceBW * fixBW);
320
// Mass according to a Breit-Wigner quadratic in m.
321
} else if (modeBWnow == 3) {
322
m2Now = m0Save*m0Save + m0Save * mWidthSave
323
* tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
324
mNow = sqrtpos( m2Now);
326
// Ditto, but m_0 Gamma_0 -> m Gamma(m) with threshold factor as above.
328
double mwNow, fixBW, runBW;
329
double m2Ref = m0Save * m0Save;
330
double mwRef = m0Save * mWidthSave;
331
double m2Thr = mThr * mThr;
333
m2Now = m2Ref + mwRef * tan( atanLow + atanDif
334
* particleDataPtr->rndmPtr->flat() );
335
mNow = sqrtpos( m2Now);
336
mwNow = mNow * mWidthSave
337
* sqrtpos( (m2Now - m2Thr) / (m2Ref - m2Thr) );
338
fixBW = mwRef / (pow2(m2Now - m2Ref) + pow2(mwRef));
339
runBW = mwNow / (pow2(m2Now - m2Ref) + pow2(mwNow));
340
} while (runBW < particleDataPtr->rndmPtr->flat()
341
* particleDataPtr->maxEnhanceBW * fixBW);
348
//--------------------------------------------------------------------------
350
// Function to calculate running mass at given mass scale.
352
double ParticleDataEntry::mRun(double mHat) {
354
// Except for six quarks return nominal mass.
355
if (idSave > 6) return m0Save;
356
double mQRun = particleDataPtr->mQRun[idSave];
357
double Lam5 = particleDataPtr->Lambda5Run;
359
// For d, u, s quarks start running at 2 GeV (RPP 2006 p. 505).
360
if (idSave < 4) return mQRun * pow ( log(2. / Lam5)
361
/ log(max(2., mHat) / Lam5), 12./23.);
363
// For c, b and t quarks start running at respective mass.
364
return mQRun * pow ( log(mQRun / Lam5)
365
/ log(max(mQRun, mHat) / Lam5), 12./23.);
369
//--------------------------------------------------------------------------
371
// Rescale all branching ratios to assure normalization to unity.
373
void ParticleDataEntry::rescaleBR(double newSumBR) {
375
// Sum up branching ratios. Find rescaling factor. Rescale.
376
double oldSumBR = 0.;
377
for ( int i = 0; i < int(channels.size()); ++ i)
378
oldSumBR += channels[i].bRatio();
379
double rescaleFactor = newSumBR / oldSumBR;
380
for ( int i = 0; i < int(channels.size()); ++ i)
381
channels[i].rescaleBR(rescaleFactor);
385
//--------------------------------------------------------------------------
387
// Prepare to pick a decay channel.
389
bool ParticleDataEntry::preparePick(int idSgn, double mHat, int idInFlav) {
391
// Reset sum of allowed widths/branching ratios.
394
// For resonances the widths are calculated dynamically.
395
if (resonancePtr != 0) {
396
resonancePtr->widthStore(idSgn, mHat, idInFlav);
397
for (int i = 0; i < int(channels.size()); ++i)
398
currentBRSum += channels[i].currentBR();
400
// Else use normal fixed branching ratios.
404
for (int i = 0; i < int(channels.size()); ++i) {
405
onMode = channels[i].onMode();
407
if ( idSgn > 0 && (onMode == 1 || onMode == 2) )
408
currentBRNow = channels[i].bRatio();
409
else if ( idSgn < 0 && (onMode == 1 || onMode == 3) )
410
currentBRNow = channels[i].bRatio();
411
channels[i].currentBR(currentBRNow);
412
currentBRSum += currentBRNow;
416
// Failure if no channels found with positive branching ratios.
417
return (currentBRSum > 0.);
421
//--------------------------------------------------------------------------
423
// Pick a decay channel according to branching ratios from preparePick.
425
DecayChannel& ParticleDataEntry::pickChannel() {
427
// Find channel in table.
428
int size = channels.size();
429
double rndmBR = currentBRSum * particleDataPtr->rndmPtr->flat();
431
do rndmBR -= channels[++i].currentBR();
432
while (rndmBR > 0. && i < size);
434
// Emergency if no channel found. Done.
435
if (i == size) i = 0;
440
//--------------------------------------------------------------------------
442
// Access methods stored in ResonanceWidths. Could have been
443
// inline in .h, except for problems with forward declarations.
445
void ParticleDataEntry::setResonancePtr(
446
ResonanceWidths* resonancePtrIn) {
447
if (resonancePtr == resonancePtrIn) return;
448
if (resonancePtr != 0) delete resonancePtr;
449
resonancePtr = resonancePtrIn;
452
void ParticleDataEntry::resInit(Info* infoPtrIn, Settings* settingsPtrIn,
453
ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
454
if (resonancePtr != 0) resonancePtr->init(infoPtrIn, settingsPtrIn,
455
particleDataPtrIn, couplingsPtrIn);
458
double ParticleDataEntry::resWidth(int idSgn, double mHat, int idIn,
459
bool openOnly, bool setBR) {
460
return (resonancePtr != 0) ? resonancePtr->width( idSgn, mHat,
461
idIn, openOnly, setBR) : 0.;
464
double ParticleDataEntry::resWidthOpen(int idSgn, double mHat, int idIn) {
465
return (resonancePtr != 0) ? resonancePtr->widthOpen( idSgn, mHat, idIn)
469
double ParticleDataEntry::resWidthStore(int idSgn, double mHat, int idIn) {
470
return (resonancePtr != 0) ? resonancePtr->widthStore( idSgn, mHat, idIn)
474
double ParticleDataEntry::resOpenFrac(int idSgn) {
475
return (resonancePtr != 0) ? resonancePtr->openFrac(idSgn) : 1.;
478
double ParticleDataEntry::resWidthRescaleFactor() {
479
return (resonancePtr != 0) ? resonancePtr->widthRescaleFactor() : 1.;
482
double ParticleDataEntry::resWidthChan(double mHat, int idAbs1,
484
return (resonancePtr != 0) ? resonancePtr->widthChan( mHat, idAbs1,
488
//--------------------------------------------------------------------------
490
// Constituent masses for (d, u, s, c, b) quarks and diquarks.
491
// Hardcoded in CONSTITUENTMASSTABLE so that they are not overwritten
492
// by mistake, and separated from the "normal" masses.
493
// Called both by setDefaults and setM0 so kept as separate method.
495
void ParticleDataEntry::setConstituentMass() {
497
// Equate with the normal masses as default guess.
498
constituentMassSave = m0Save;
500
// Quark masses trivial. Also gluon mass.
501
if (idSave < 6) constituentMassSave = CONSTITUENTMASSTABLE[idSave];
502
if (idSave == 21) constituentMassSave = CONSTITUENTMASSTABLE[9];
504
// Diquarks as simple sum of constituent quarks.
505
if (idSave > 1000 && idSave < 10000 && (idSave/10)%10 == 0) {
506
int id1 = idSave/1000;
507
int id2 = (idSave/100)%10;
508
if (id1 <6 && id2 < 6) constituentMassSave
509
= CONSTITUENTMASSTABLE[id1] + CONSTITUENTMASSTABLE[id2];
514
//--------------------------------------------------------------------------
516
// Convert string to lowercase for case-insensitive comparisons.
518
string ParticleDataEntry::toLower(const string& nameConv) {
520
string temp(nameConv);
521
for (int i = 0; i < int(temp.length()); ++i)
522
temp[i] = std::tolower(temp[i]);
527
//==========================================================================
529
// ParticleData class.
530
// This class holds a map of all ParticleDataEntries,
531
// each entry containing info on a particle species.
533
//--------------------------------------------------------------------------
535
// Get data to be distributed among particles during setup.
536
// Note: this routine is called twice. Firstly from init(...), but
537
// the data should not be used at that point, so is likely overkill.
538
// Secondly, from initWidths, after user has had time to change.
540
void ParticleData::initCommon() {
542
// Mass generation: fixed mass or linear/quadratic Breit-Wigner.
543
modeBreitWigner = settingsPtr->mode("ParticleData:modeBreitWigner");
545
// Maximum tail enhancement when adding threshold factor to Breit-Wigner.
546
maxEnhanceBW = settingsPtr->parm("ParticleData:maxEnhanceBW");
548
// Find initial MSbar masses for five light flavours.
549
mQRun[1] = settingsPtr->parm("ParticleData:mdRun");
550
mQRun[2] = settingsPtr->parm("ParticleData:muRun");
551
mQRun[3] = settingsPtr->parm("ParticleData:msRun");
552
mQRun[4] = settingsPtr->parm("ParticleData:mcRun");
553
mQRun[5] = settingsPtr->parm("ParticleData:mbRun");
554
mQRun[6] = settingsPtr->parm("ParticleData:mtRun");
556
// Find Lambda5 value to use in running of MSbar masses.
557
double alphaSvalue = settingsPtr->parm("ParticleData:alphaSvalueMRun");
559
alphaS.init( alphaSvalue, 1);
560
Lambda5Run = alphaS.Lambda5();
564
//--------------------------------------------------------------------------
566
// Initialize pointer for particles to the full database, the Breit-Wigners
567
// of normal hadrons and the ResonanceWidths of resonances. For the latter
568
// the order of initialization is essential to get secondary widths right.
570
void ParticleData::initWidths( vector<ResonanceWidths*> resonancePtrs) {
572
// Initialize some common data.
575
// Pointer to database and Breit-Wigner mass initialization for each
577
ResonanceWidths* resonancePtr = 0;
578
for (map<int, ParticleDataEntry>::iterator pdtEntry = pdt.begin();
579
pdtEntry != pdt.end(); ++pdtEntry) {
580
ParticleDataEntry& pdtNow = pdtEntry->second;
581
pdtNow.initPtr( this);
584
// Remove any existing resonances.
585
resonancePtr = pdtNow.getResonancePtr();
586
if (resonancePtr != 0) pdtNow.setResonancePtr(0);
589
// Begin set up new resonance objects.
590
// Z0, W+- and top are almost always needed.
591
resonancePtr = new ResonanceGmZ(23);
592
setResonancePtr( 23, resonancePtr);
593
resonancePtr = new ResonanceW(24);
594
setResonancePtr( 24, resonancePtr);
595
resonancePtr = new ResonanceTop(6);
596
setResonancePtr( 6, resonancePtr);
599
if (!settingsPtr->flag("Higgs:useBSM")) {
600
resonancePtr = new ResonanceH(0, 25);
601
setResonancePtr( 25, resonancePtr);
605
resonancePtr = new ResonanceH(1, 25);
606
setResonancePtr( 25, resonancePtr);
607
resonancePtr = new ResonanceH(2, 35);
608
setResonancePtr( 35, resonancePtr);
609
resonancePtr = new ResonanceH(3, 36);
610
setResonancePtr( 36, resonancePtr);
611
resonancePtr = new ResonanceHchg(37);
612
setResonancePtr( 37, resonancePtr);
615
// A fourth generation: b', t', tau', nu'_tau.
616
resonancePtr = new ResonanceFour(7);
617
setResonancePtr( 7, resonancePtr);
618
resonancePtr = new ResonanceFour(8);
619
setResonancePtr( 8, resonancePtr);
620
resonancePtr = new ResonanceFour(17);
621
setResonancePtr( 17, resonancePtr);
622
resonancePtr = new ResonanceFour(18);
623
setResonancePtr( 18, resonancePtr);
625
// New gauge bosons: Z', W', R.
626
resonancePtr = new ResonanceZprime(32);
627
setResonancePtr( 32, resonancePtr);
628
resonancePtr = new ResonanceWprime(34);
629
setResonancePtr( 34, resonancePtr);
630
resonancePtr = new ResonanceRhorizontal(41);
631
setResonancePtr( 41, resonancePtr);
634
resonancePtr = new ResonanceLeptoquark(42);
635
setResonancePtr( 42, resonancePtr);
639
for(int i = 1; i < 7; i++){
640
resonancePtr = new ResonanceSquark(1000000 + i);
641
setResonancePtr( 1000000 + i, resonancePtr);
642
resonancePtr = new ResonanceSquark(2000000 + i);
643
setResonancePtr( 2000000 + i, resonancePtr);
646
// - Sleptons and sneutrinos
647
for(int i = 1; i < 7; i++){
648
resonancePtr = new ResonanceSlepton(1000010 + i);
649
setResonancePtr( 1000010 + i, resonancePtr);
650
resonancePtr = new ResonanceSlepton(2000010 + i);
651
setResonancePtr( 2000010 + i, resonancePtr);
655
resonancePtr = new ResonanceGluino(1000021);
656
setResonancePtr( 1000021, resonancePtr);
659
resonancePtr = new ResonanceChar(1000024);
660
setResonancePtr( 1000024, resonancePtr);
661
resonancePtr = new ResonanceChar(1000037);
662
setResonancePtr( 1000037, resonancePtr);
665
resonancePtr = new ResonanceNeut(1000022);
666
setResonancePtr( 1000022, resonancePtr);
667
resonancePtr = new ResonanceNeut(1000023);
668
setResonancePtr( 1000023, resonancePtr);
669
resonancePtr = new ResonanceNeut(1000025);
670
setResonancePtr( 1000025, resonancePtr);
671
resonancePtr = new ResonanceNeut(1000035);
672
setResonancePtr( 1000035, resonancePtr);
673
resonancePtr = new ResonanceNeut(1000045);
674
setResonancePtr( 1000045, resonancePtr);
676
// Excited quarks and leptons.
677
for (int i = 1; i < 7; ++i) {
678
resonancePtr = new ResonanceExcited(4000000 + i);
679
setResonancePtr( 4000000 + i, resonancePtr);
681
for (int i = 11; i < 17; ++i) {
682
resonancePtr = new ResonanceExcited(4000000 + i);
683
setResonancePtr( 4000000 + i, resonancePtr);
686
// An excited graviton/gluon in extra-dimensional scenarios.
687
resonancePtr = new ResonanceGraviton(5100039);
688
setResonancePtr( 5100039, resonancePtr);
689
resonancePtr = new ResonanceKKgluon(5100021);
690
setResonancePtr( 5100021, resonancePtr);
692
// A left-right-symmetric scenario with new righthanded neutrinos,
693
// righthanded gauge bosons and doubly charged Higgses.
694
resonancePtr = new ResonanceNuRight(9900012);
695
setResonancePtr( 9900012, resonancePtr);
696
resonancePtr = new ResonanceNuRight(9900014);
697
setResonancePtr( 9900014, resonancePtr);
698
resonancePtr = new ResonanceNuRight(9900016);
699
setResonancePtr( 9900016, resonancePtr);
700
resonancePtr = new ResonanceZRight(9900023);
701
setResonancePtr( 9900023, resonancePtr);
702
resonancePtr = new ResonanceWRight(9900024);
703
setResonancePtr( 9900024, resonancePtr);
704
resonancePtr = new ResonanceHchgchgLeft(9900041);
705
setResonancePtr( 9900041, resonancePtr);
706
resonancePtr = new ResonanceHchgchgRight(9900042);
707
setResonancePtr( 9900042, resonancePtr);
709
// Attach user-defined external resonances and do basic initialization.
710
for (int i = 0; i < int(resonancePtrs.size()); ++i) {
711
int idNow = resonancePtrs[i]->id();
712
resonancePtrs[i]->initBasic(idNow);
713
setResonancePtr( idNow, resonancePtrs[i]);
716
// Set up lists to order resonances in ascending mass.
717
vector<int> idOrdered;
718
vector<double> m0Ordered;
720
// Put Z0 and W+- first, since known to be SM and often off-shell.
721
idOrdered.push_back(23);
722
m0Ordered.push_back(m0(23));
723
idOrdered.push_back(24);
724
m0Ordered.push_back(m0(24));
726
// Loop through particle table to find resonances.
727
for (map<int, ParticleDataEntry>::iterator pdtEntry = pdt.begin();
728
pdtEntry != pdt.end(); ++pdtEntry) {
729
ParticleDataEntry& pdtNow = pdtEntry->second;
730
int idNow = pdtNow.id();
732
// Set up a simple default object for uninitialized resonances.
733
if (pdtNow.isResonance() && pdtNow.getResonancePtr() == 0) {
734
resonancePtr = new ResonanceGeneric(idNow);
735
setResonancePtr( idNow, resonancePtr);
738
// Insert resonances in ascending mass, to respect decay hierarchies.
739
if (pdtNow.getResonancePtr() != 0 && idNow != 23 && idNow != 24) {
740
double m0Now = pdtNow.m0();
741
idOrdered.push_back(idNow);
742
m0Ordered.push_back(m0Now);
743
for (int i = int(idOrdered.size()) - 2; i > 1; --i) {
744
if (m0Ordered[i] < m0Now) break;
745
swap( idOrdered[i], idOrdered[i + 1]);
746
swap( m0Ordered[i], m0Ordered[i + 1]);
751
// Initialize the resonances in ascending mass order.
752
for (int i = 0; i < int(idOrdered.size()); ++i) resInit( idOrdered[i]);
756
//--------------------------------------------------------------------------
758
// Read in database from specific XML file (which may refer to others).
760
bool ParticleData::readXML(string inFile, bool reset) {
762
// Normally reset whole database before beginning.
763
if (reset) {pdt.clear(); isInit = false;}
765
// List of files to be checked.
766
vector<string> files;
767
files.push_back(inFile);
769
// Loop over files. Open them for read.
770
for (int i = 0; i < int(files.size()); ++i) {
771
const char* cstring = files[i].c_str();
772
ifstream is(cstring);
774
// Check that instream is OK.
776
infoPtr->errorMsg("Error in ParticleData::readXML:"
777
" did not find file", files[i]);
781
// Read in one line at a time.
784
while ( getline(is, line) ) {
786
// Get first word of a line.
787
istringstream getfirst(line);
791
// Check for occurence of a particle. Add any continuation lines.
792
if (word1 == "<particle") {
793
while (line.find(">") == string::npos) {
795
getline(is, addLine);
799
// Read in particle properties.
800
int idTmp = intAttributeValue( line, "id");
801
string nameTmp = attributeValue( line, "name");
802
string antiNameTmp = attributeValue( line, "antiName");
803
if (antiNameTmp == "") antiNameTmp = "void";
804
int spinTypeTmp = intAttributeValue( line, "spinType");
805
int chargeTypeTmp = intAttributeValue( line, "chargeType");
806
int colTypeTmp = intAttributeValue( line, "colType");
807
double m0Tmp = doubleAttributeValue( line, "m0");
808
double mWidthTmp = doubleAttributeValue( line, "mWidth");
809
double mMinTmp = doubleAttributeValue( line, "mMin");
810
double mMaxTmp = doubleAttributeValue( line, "mMax");
811
double tau0Tmp = doubleAttributeValue( line, "tau0");
813
// Erase if particle already exists.
814
if (isParticle(idTmp)) pdt.erase(idTmp);
816
// Store new particle. Save pointer, to be used for decay channels.
817
addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
818
colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
819
particlePtr = particleDataEntryPtr(idTmp);
821
// Check for occurence of a decay channel. Add any continuation lines.
822
} else if (word1 == "<channel") {
823
while (line.find(">") == string::npos) {
825
getline(is, addLine);
829
// Read in channel properties - products so far only as a string.
830
int onMode = intAttributeValue( line, "onMode");
831
double bRatio = doubleAttributeValue( line, "bRatio");
832
int meMode = intAttributeValue( line, "meMode");
833
string products = attributeValue( line, "products");
835
// Read in decay products from stream. Must have at least one.
836
istringstream prodStream(products);
837
int prod0 = 0; int prod1 = 0; int prod2 = 0; int prod3 = 0;
838
int prod4 = 0; int prod5 = 0; int prod6 = 0; int prod7 = 0;
839
prodStream >> prod0 >> prod1 >> prod2 >> prod3 >> prod4 >> prod5
842
infoPtr->errorMsg("Error in ParticleData::readXML:"
843
" incomplete decay channel", line);
847
// Store new channel (if particle already known).
848
if (particlePtr == 0) {
849
infoPtr->errorMsg("Error in ParticleData::readXML:"
850
" orphan decay channel", line);
853
particlePtr->addChannel(onMode, bRatio, meMode, prod0, prod1,
854
prod2, prod3, prod4, prod5, prod6, prod7);
856
// Check for occurence of a file also to be read.
857
} else if (word1 == "<file") {
858
string file = attributeValue(line, "name");
860
infoPtr->errorMsg("Warning in ParticleData::readXML:"
861
" skip unrecognized file name", line);
862
} else files.push_back(file);
865
// End of loop over lines in input file and loop over files.
869
// All particle data at this stage defines baseline original.
870
if (reset) for (map<int, ParticleDataEntry>::iterator pdtEntry
871
= pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
872
particlePtr = &pdtEntry->second;
873
particlePtr->setHasChanged(false);
882
//--------------------------------------------------------------------------
884
// Print out complete database in numerical order as an XML file.
886
void ParticleData::listXML(string outFile) {
888
// Convert file name to ofstream.
889
const char* cstring = outFile.c_str();
890
ofstream os(cstring);
892
// Iterate through the particle data table.
893
for (map<int, ParticleDataEntry>::iterator pdtEntry
894
= pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
895
particlePtr = &pdtEntry->second;
897
// Print particle properties.
898
os << "<particle id=\"" << particlePtr->id() << "\""
899
<< " name=\"" << particlePtr->name() << "\"";
900
if (particlePtr->hasAnti())
901
os << " antiName=\"" << particlePtr->name(-1) << "\"";
902
os << " spinType=\"" << particlePtr->spinType() << "\""
903
<< " chargeType=\"" << particlePtr->chargeType() << "\""
904
<< " colType=\"" << particlePtr->colType() << "\"\n";
905
// Pick format for mass and width based on mass value.
906
double m0Now = particlePtr->m0();
907
if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
908
os << fixed << setprecision(5);
909
else os << scientific << setprecision(3);
910
os << " m0=\"" << m0Now << "\"";
911
if (particlePtr->mWidth() > 0.)
912
os << " mWidth=\"" << particlePtr->mWidth() << "\""
913
<< " mMin=\"" << particlePtr->mMin() << "\""
914
<< " mMax=\"" << particlePtr->mMax() << "\"";
915
if (particlePtr->tau0() > 0.) os << scientific << setprecision(5)
916
<< " tau0=\"" << particlePtr->tau0() << "\"";
919
// Loop through the decay channel table for each particle.
920
if (particlePtr->sizeChannels() > 0) {
921
for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
922
const DecayChannel& channel = particlePtr->channel(i);
923
int mult = channel.multiplicity();
925
// Print decay channel properties.
926
os << " <channel onMode=\"" << channel.onMode() << "\""
927
<< fixed << setprecision(7)
928
<< " bRatio=\"" << channel.bRatio() << "\"";
929
if (channel.meMode() > 0)
930
os << " meMode=\"" << channel.meMode() << "\"";
931
os << " products=\"";
932
for (int j = 0; j < mult; ++j) {
933
os << channel.product(j);
934
if (j < mult - 1) os << " ";
937
// Finish off line and loop over allowed decay channels.
942
// Finish off existing particle.
943
os << "</particle>\n\n";
949
//--------------------------------------------------------------------------
951
// Read in database from specific free format file.
953
bool ParticleData::readFF(string inFile, bool reset) {
955
// Normally reset whole database before beginning.
956
if (reset) {pdt.clear(); isInit = false;}
958
// Open file for read and check that instream is OK.
959
const char* cstring = inFile.c_str();
960
ifstream is(cstring);
962
infoPtr->errorMsg("Error in ParticleData::readFF:"
963
" did not find file", inFile);
967
// Read in one line at a time.
970
bool readParticle = false;
971
while ( getline(is, line) ) {
973
// Empty lines begins new particle.
974
if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) {
979
// Prepare to use standard read from line.
980
istringstream readLine(line);
982
// Read in a line with particle information.
985
// Properties to be read.
987
string nameTmp, antiNameTmp;
988
int spinTypeTmp, chargeTypeTmp, colTypeTmp;
989
double m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp;
993
readLine >> idTmp >> nameTmp >> antiNameTmp >> spinTypeTmp
994
>> chargeTypeTmp >> colTypeTmp >> m0Tmp >> mWidthTmp
995
>> mMinTmp >> mMaxTmp >> tau0Tmp;
997
// Error printout if something went wrong.
999
infoPtr->errorMsg("Error in ParticleData::readFF:"
1000
" incomplete particle", line);
1004
// Erase if particle already exists.
1005
if (isParticle(idTmp)) pdt.erase(idTmp);
1007
// Store new particle. Save pointer, to be used for decay channels.
1008
addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1009
colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1010
particlePtr = particleDataEntryPtr(idTmp);
1011
readParticle = false;
1013
// Read in a line with decay channel information.
1016
// Properties to be read.
1020
int prod0 = 0; int prod1 = 0; int prod2 = 0; int prod3 = 0;
1021
int prod4 = 0; int prod5 = 0; int prod6 = 0; int prod7 = 0;
1023
// Read in data from stream. Need at least one decay product.
1024
readLine >> onMode >> bRatio >> meMode >> prod0;
1026
infoPtr->errorMsg("Error in ParticleData::readFF:"
1027
" incomplete decay channel", line);
1030
readLine >> prod1 >> prod2 >> prod3 >> prod4 >> prod5
1033
// Store new channel.
1034
if (particlePtr == 0) {
1035
infoPtr->errorMsg("Error in ParticleData::readFF:"
1036
" orphan decay channel", line);
1039
particlePtr->addChannel(onMode, bRatio, meMode, prod0, prod1,
1040
prod2, prod3, prod4, prod5, prod6, prod7);
1044
// End of while loop over lines in the file.
1054
//--------------------------------------------------------------------------
1056
// Print out complete database in numerical order as a free format file.
1058
void ParticleData::listFF(string outFile) {
1060
// Convert file name to ofstream.
1061
const char* cstring = outFile.c_str();
1062
ofstream os(cstring);
1064
// Iterate through the particle data table.
1065
for (map<int, ParticleDataEntry>::iterator pdtEntry
1066
= pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1067
particlePtr = &pdtEntry->second;
1069
// Pick format for mass and width based on mass value.
1070
double m0Now = particlePtr->m0();
1071
if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1072
os << fixed << setprecision(5);
1073
else os << scientific << setprecision(3);
1075
// Print particle properties.
1076
os << "\n" << setw(8) << particlePtr->id() << " "
1077
<< left << setw(16) << particlePtr->name() << " "
1078
<< setw(16) << particlePtr->name(-1) << " "
1079
<< right << setw(2) << particlePtr->spinType() << " "
1080
<< setw(2) << particlePtr->chargeType() << " "
1081
<< setw(2) << particlePtr->colType() << " "
1082
<< setw(10) << particlePtr->m0() << " "
1083
<< setw(10) << particlePtr->mWidth() << " "
1084
<< setw(10) << particlePtr->mMin() << " "
1085
<< setw(10) << particlePtr->mMax() << " "
1086
<< scientific << setprecision(5)
1087
<< setw(12) << particlePtr->tau0() << "\n";
1089
// Loop through the decay channel table for each particle.
1090
if (particlePtr->sizeChannels() > 0) {
1091
for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1092
const DecayChannel& channel = particlePtr->channel(i);
1093
os << " " << setw(6) << channel.onMode()
1094
<< " " << fixed << setprecision(7) << setw(10)
1095
<< channel.bRatio() << " "
1096
<< setw(3) << channel.meMode() << " ";
1097
for (int j = 0; j < channel.multiplicity(); ++j)
1098
os << setw(8) << channel.product(j) << " ";
1107
//--------------------------------------------------------------------------
1109
// Read in updates from a character string, like a line of a file.
1110
// Is used by readString (and readFile) in Pythia.
1112
bool ParticleData::readString(string lineIn, bool warn, ostream& os) {
1114
// If empty line then done.
1115
if (lineIn.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return true;
1117
// Take copy that will be modified.
1118
string line = lineIn;
1120
// If first character is not a digit then taken to be a comment.
1121
int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
1122
if (!isdigit(line[firstChar])) return true;
1124
// Replace colons and equal signs by blanks to make parsing simpler.
1125
for ( int j = 0; j < int(line.size()); ++ j)
1126
if (line[j] == ':' || line[j] == '=') line[j] = ' ';
1128
// Get particle id and property name.
1131
istringstream getWord(line);
1132
getWord >> idTmp >> property;
1133
property = toLower(property);
1135
// Check that valid particle.
1136
if ( (!isParticle(idTmp) && property != "all" && property != "new")
1138
if (warn) os << "\n Warning: input particle not found in Particle"
1139
<< " Data Table; skip:\n " << lineIn << "\n";
1143
// Identify particle property and read + set its value, case by case.
1144
if (property == "name") {
1147
pdt[idTmp].setName(nameTmp);
1150
if (property == "antiname") {
1152
getWord >> antiNameTmp;
1153
pdt[idTmp].setAntiName(antiNameTmp);
1156
if (property == "names") {
1157
string nameTmp, antiNameTmp;
1158
getWord >> nameTmp >> antiNameTmp;
1159
pdt[idTmp].setNames(nameTmp, antiNameTmp);
1162
if (property == "spintype") {
1164
getWord >> spinTypeTmp;
1165
pdt[idTmp].setSpinType(spinTypeTmp);
1168
if (property == "chargetype") {
1170
getWord >> chargeTypeTmp;
1171
pdt[idTmp].setChargeType(chargeTypeTmp);
1174
if (property == "coltype") {
1176
getWord >> colTypeTmp;
1177
pdt[idTmp].setColType(colTypeTmp);
1180
if (property == "m0") {
1183
pdt[idTmp].setM0(m0Tmp);
1186
if (property == "mwidth") {
1188
getWord >> mWidthTmp;
1189
pdt[idTmp].setMWidth(mWidthTmp);
1192
if (property == "mmin") {
1195
pdt[idTmp].setMMin(mMinTmp);
1198
if (property == "mmax") {
1201
pdt[idTmp].setMMax(mMaxTmp);
1204
if (property == "tau0") {
1207
pdt[idTmp].setTau0(tau0Tmp);
1210
if (property == "isresonance") {
1212
getWord >> isresTmp;
1213
bool isResonanceTmp = boolString(isresTmp);
1214
pdt[idTmp].setIsResonance(isResonanceTmp);
1217
if (property == "maydecay") {
1220
bool mayDecayTmp = boolString(mayTmp);
1221
pdt[idTmp].setMayDecay(mayDecayTmp);
1224
if (property == "doexternaldecay") {
1226
getWord >> extdecTmp;
1227
bool doExternalDecayTmp = boolString(extdecTmp);
1228
pdt[idTmp].setDoExternalDecay(doExternalDecayTmp);
1231
if (property == "isvisible") {
1233
getWord >> isvisTmp;
1234
bool isVisibleTmp = boolString(isvisTmp);
1235
pdt[idTmp].setIsVisible(isVisibleTmp);
1238
if (property == "doforcewidth") {
1240
getWord >> doforceTmp;
1241
bool doForceWidthTmp = boolString(doforceTmp);
1242
pdt[idTmp].setDoForceWidth(doForceWidthTmp);
1246
// Addition or complete replacement of a particle.
1247
if (property == "all" || property == "new") {
1249
// Default values for properties to be read.
1250
string nameTmp = "void";
1251
string antiNameTmp = "void";
1252
int spinTypeTmp = 0;
1253
int chargeTypeTmp = 0;
1256
double mWidthTmp = 0.;
1257
double mMinTmp = 0.;
1258
double mMaxTmp = 0.;
1259
double tau0Tmp = 0.;
1261
// Read in data from stream.
1262
getWord >> nameTmp >> antiNameTmp >> spinTypeTmp >> chargeTypeTmp
1263
>> colTypeTmp >> m0Tmp >> mWidthTmp >> mMinTmp >> mMaxTmp
1266
// To keep existing decay channels, only overwrite particle data.
1267
if (property == "all" && isParticle(idTmp)) {
1268
setAll( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1269
colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1271
// Else start over completely from scratch.
1273
if (isParticle(idTmp)) pdt.erase(idTmp);
1274
addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1275
colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1280
// Set onMode of all decay channels in one go.
1281
if (property == "onmode") {
1284
getWord >> onModeIn;
1285
// For onMode allow the optional possibility of Bool input.
1286
if (isdigit(onModeIn[0])) {
1287
istringstream getOnMode(onModeIn);
1288
getOnMode >> onMode;
1289
} else onMode = (boolString(onModeIn)) ? 1 : 0;
1290
for (int i = 0; i < pdt[idTmp].sizeChannels(); ++i)
1291
pdt[idTmp].channel(i).onMode(onMode);
1295
// Selective search for matching decay channels.
1297
if (property == "offifany" || property == "onifany" ||
1298
property == "onposifany" || property == "onnegifany") matchKind = 1;
1299
if (property == "offifall" || property == "onifall" ||
1300
property == "onposifall" || property == "onnegifall") matchKind = 2;
1301
if (property == "offifmatch" || property == "onifmatch" ||
1302
property == "onposifmatch" || property == "onnegifmatch") matchKind = 3;
1303
if (matchKind > 0) {
1305
if (property == "onifany" || property == "onifall"
1306
|| property == "onifmatch") onMode = 1;
1307
if (property == "onposifany" || property == "onposifall"
1308
|| property == "onposifmatch") onMode = 2;
1309
if (property == "onnegifany" || property == "onnegifall"
1310
|| property == "onnegifmatch") onMode = 3;
1312
// Read in particles to match.
1313
vector<int> idToMatch;
1317
if (!getWord) break;
1318
idToMatch.push_back(abs(idRead));
1320
int nToMatch = idToMatch.size();
1322
// Loop over all decay channels.
1323
for (int i = 0; i < pdt[idTmp].sizeChannels(); ++i) {
1324
int multi = pdt[idTmp].channel(i).multiplicity();
1326
// Look for any match at all.
1327
if (matchKind == 1) {
1328
bool foundMatch = false;
1329
for (int j = 0; j < multi; ++j) {
1330
int idNow = abs(pdt[idTmp].channel(i).product(j));
1331
for (int k = 0; k < nToMatch; ++k)
1332
if (idNow == idToMatch[k]) {foundMatch = true; break;}
1333
if (foundMatch) break;
1335
if (foundMatch) pdt[idTmp].channel(i).onMode(onMode);
1337
// Look for match of all products provided.
1339
int nUnmatched = nToMatch;
1340
if (multi < nToMatch);
1341
else if (multi > nToMatch && matchKind == 3);
1343
vector<int> idUnmatched;
1344
for (int k = 0; k < nToMatch; ++k)
1345
idUnmatched.push_back(idToMatch[k]);
1346
for (int j = 0; j < multi; ++j) {
1347
int idNow = abs(pdt[idTmp].channel(i).product(j));
1348
for (int k = 0; k < nUnmatched; ++k)
1349
if (idNow == idUnmatched[k]) {
1350
idUnmatched[k] = idUnmatched[--nUnmatched];
1353
if (nUnmatched == 0) break;
1356
if (nUnmatched == 0) pdt[idTmp].channel(i).onMode(onMode);
1362
// Rescale all branching ratios by common factor.
1363
if (property == "rescalebr") {
1366
pdt[idTmp].rescaleBR(factor);
1370
// Reset decay table in preparation for new input.
1371
if (property == "onechannel") pdt[idTmp].clearChannels();
1373
// Add or change a decay channel: get channel number and new property.
1374
if (property == "addchannel" || property == "onechannel"
1375
|| isdigit(property[0])) {
1377
if (property == "addchannel" || property == "onechannel") {
1378
pdt[idTmp].addChannel();
1379
channel = pdt[idTmp].sizeChannels() - 1;
1382
istringstream getChannel(property);
1383
getChannel >> channel;
1384
getWord >> property;
1385
property = toLower(property);
1388
// Check that channel exists.
1389
if (channel < 0 || channel >= pdt[idTmp].sizeChannels()) return false;
1391
// Find decay channel property and value, case by case.
1392
// At same time also do case where all should be replaced.
1393
if (property == "onmode" || property == "all") {
1396
getWord >> onModeIn;
1397
// For onMode allow the optional possibility of Bool input.
1398
if (isdigit(onModeIn[0])) {
1399
istringstream getOnMode(onModeIn);
1400
getOnMode >> onMode;
1401
} else onMode = (boolString(onModeIn)) ? 1 : 0;
1402
pdt[idTmp].channel(channel).onMode(onMode);
1403
if (property == "onmode") return true;
1405
if (property == "bratio" || property == "all") {
1408
pdt[idTmp].channel(channel).bRatio(bRatio);
1409
if (property == "bratio") return true;
1411
if (property == "memode" || property == "all") {
1414
pdt[idTmp].channel(channel).meMode(meMode);
1415
if (property == "memode") return true;
1418
// Scan for products until end of line.
1419
if (property == "products" || property == "all") {
1421
for (int iProd = 0; iProd < 8; ++iProd) {
1424
if (!getWord) break;
1425
pdt[idTmp].channel(channel).product(iProd, idProd);
1428
for (int iProd = nProd; iProd < 8; ++iProd)
1429
pdt[idTmp].channel(channel).product(iProd, 0);
1430
pdt[idTmp].channel(channel).multiplicity(nProd);
1434
// Rescale an existing branching ratio.
1435
if (property == "rescalebr") {
1438
pdt[idTmp].channel(channel).rescaleBR(factor);
1443
// Return false if failed to recognize property.
1444
if (warn) os << "\n Warning: input property not found in Particle"
1445
<< " Data Table; skip:\n " << lineIn << "\n";
1450
//--------------------------------------------------------------------------
1452
// Print out complete or changed table of database in numerical order.
1454
void ParticleData::list(bool changedOnly, bool changedRes, ostream& os) {
1456
// Table header; output for bool as off/on.
1458
os << "\n -------- PYTHIA Particle Data Table (complete) --------"
1459
<< "------------------------------------------------------------"
1460
<< "--------------\n \n";
1463
os << "\n -------- PYTHIA Particle Data Table (changed only) ----"
1464
<< "------------------------------------------------------------"
1465
<< "--------------\n \n";
1467
os << " id name antiName spn chg col m0"
1468
<< " mWidth mMin mMax tau0 res dec ext "
1469
<< "vis wid\n no onMode bRatio meMode products \n";
1471
// Iterate through the particle data table. Option to skip unchanged.
1473
for (map<int, ParticleDataEntry>::iterator pdtEntry
1474
= pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1475
particlePtr = &pdtEntry->second;
1476
if ( !changedOnly || particlePtr->hasChanged() ||
1477
( changedRes && particlePtr->getResonancePtr() != 0 ) ) {
1479
// Pick format for mass and width based on mass value.
1480
double m0Now = particlePtr->m0();
1481
if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1482
os << fixed << setprecision(5);
1483
else os << scientific << setprecision(3);
1485
// Print particle properties.
1487
string antiNameTmp = particlePtr->name(-1);
1488
if (antiNameTmp == "void") antiNameTmp = " ";
1489
os << "\n" << setw(8) << particlePtr->id() << " "
1490
<< left << setw(16) << particlePtr->name() << " "
1491
<< setw(16) << antiNameTmp << " "
1492
<< right << setw(2) << particlePtr->spinType() << " "
1493
<< setw(2) << particlePtr->chargeType() << " "
1494
<< setw(2) << particlePtr->colType() << " "
1495
<< setw(10) << particlePtr->m0() << " "
1496
<< setw(10) << particlePtr->mWidth() << " "
1497
<< setw(10) << particlePtr->mMin() << " "
1498
<< setw(10) << particlePtr->mMax() << " "
1499
<< scientific << setprecision(5)
1500
<< setw(12) << particlePtr->tau0() << " " << setw(2)
1501
<< particlePtr->isResonance() << " " << setw(2)
1502
<< (particlePtr->mayDecay() && particlePtr->canDecay())
1503
<< " " << setw(2) << particlePtr->doExternalDecay() << " "
1504
<< setw(2) << particlePtr->isVisible()<< " "
1505
<< setw(2) << particlePtr->doForceWidth() << "\n";
1507
// Loop through the decay channel table for each particle.
1508
if (particlePtr->sizeChannels() > 0) {
1509
for (int i = 0; i < int(particlePtr->sizeChannels()); ++i) {
1510
const DecayChannel& channel = particlePtr->channel(i);
1511
os << " " << setprecision(7)
1513
<< setw(6) << channel.onMode()
1514
<< fixed<< setw(12) << channel.bRatio()
1515
<< setw(5) << channel.meMode() << " ";
1516
for (int j = 0; j < channel.multiplicity(); ++j)
1517
os << setw(8) << channel.product(j) << " ";
1525
// End of loop over database contents.
1526
if (changedOnly && nList == 0) os << "\n no particle data has been "
1527
<< "changed from its default value \n";
1528
os << "\n -------- End PYTHIA Particle Data Table -----------------"
1529
<< "--------------------------------------------------------------"
1530
<< "----------\n" << endl;
1534
//--------------------------------------------------------------------------
1536
// Print out partial table of database in input order.
1538
void ParticleData::list(vector<int> idList, ostream& os) {
1540
// Table header; output for bool as off/on.
1541
os << "\n -------- PYTHIA Particle Data Table (partial) ---------"
1542
<< "------------------------------------------------------------"
1543
<< "--------------\n \n";
1544
os << " id name antiName spn chg col m0"
1545
<< " mWidth mMin mMax tau0 res dec ext "
1546
<< "vis wid\n no onMode bRatio meMode products \n";
1548
// Iterate through the given list of input particles.
1549
for (int i = 0; i < int(idList.size()); ++i) {
1550
particlePtr = particleDataEntryPtr(idList[i]);
1552
// Pick format for mass and width based on mass value.
1553
double m0Now = particlePtr->m0();
1554
if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1555
os << fixed << setprecision(5);
1556
else os << scientific << setprecision(3);
1558
// Print particle properties.
1559
string antiNameTmp = particlePtr->name(-1);
1560
if (antiNameTmp == "void") antiNameTmp = " ";
1561
os << "\n" << setw(8) << particlePtr->id() << " "
1562
<< left << setw(16) << particlePtr->name() << " "
1563
<< setw(16) << antiNameTmp << " "
1564
<< right << setw(2) << particlePtr->spinType() << " "
1565
<< setw(2) << particlePtr->chargeType() << " "
1566
<< setw(2) << particlePtr->colType() << " "
1567
<< setw(10) << particlePtr->m0() << " "
1568
<< setw(10) << particlePtr->mWidth() << " "
1569
<< setw(10) << particlePtr->mMin() << " "
1570
<< setw(10) << particlePtr->mMax() << " "
1571
<< scientific << setprecision(5)
1572
<< setw(12) << particlePtr->tau0() << " " << setw(2)
1573
<< particlePtr->isResonance() << " " << setw(2)
1574
<< (particlePtr->mayDecay() && particlePtr->canDecay())
1575
<< " " << setw(2) << particlePtr->doExternalDecay() << " "
1576
<< setw(2) << particlePtr->isVisible() << " "
1577
<< setw(2) << particlePtr->doForceWidth() << "\n";
1579
// Loop through the decay channel table for each particle.
1580
if (particlePtr->sizeChannels() > 0) {
1581
for (int j = 0; j < int(particlePtr->sizeChannels()); ++j) {
1582
const DecayChannel& channel = particlePtr->channel(j);
1583
os << " " << setprecision(7)
1585
<< setw(6) << channel.onMode()
1586
<< fixed<< setw(12) << channel.bRatio()
1587
<< setw(5) << channel.meMode() << " ";
1588
for (int k = 0; k < channel.multiplicity(); ++k)
1589
os << setw(8) << channel.product(k) << " ";
1596
// End of loop over database contents.
1597
os << "\n -------- End PYTHIA Particle Data Table -----------------"
1598
<< "--------------------------------------------------------------"
1599
<< "----------\n" << endl;
1603
//--------------------------------------------------------------------------
1605
// Check that table makes sense: e.g. consistent names and mass ranges,
1606
// that branching ratios sum to unity, that charge is conserved and
1607
// that phase space is open in each channel.
1608
// verbosity = 0: mimimal amount of checks, e.g. that no channels open.
1609
// = 1: further warning if individual channels closed
1610
// (except for resonances).
1611
// = 2: also print branching-ratio-averaged threshold mass.
1612
// = 11, 12: as 1, 2, but include resonances in detailed checks.
1614
void ParticleData::checkTable(int verbosity, ostream& os) {
1617
os << "\n -------- PYTHIA Check of Particle Data Table ------------"
1621
// Loop through all particles.
1622
for (map<int, ParticleDataEntry>::iterator pdtEntry
1623
= pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1624
particlePtr = &pdtEntry->second;
1626
// Extract some particle properties. Set some flags.
1627
int idNow = particlePtr->id();
1628
bool hasAntiNow = particlePtr->hasAnti();
1629
int spinTypeNow = particlePtr->spinType();
1630
int chargeTypeNow = particlePtr->chargeType();
1631
int baryonTypeNow = particlePtr->baryonNumberType();
1632
double m0Now = particlePtr->m0();
1633
double mMinNow = particlePtr->mMin();
1634
double mMaxNow = particlePtr->mMax();
1635
double mWidthNow = particlePtr->mWidth();
1636
double tau0Now = particlePtr->tau0();
1637
bool isResonanceNow = particlePtr->isResonance();
1638
bool hasPrinted = false;
1639
bool studyCloser = verbosity > 10 || !isResonanceNow;
1641
// Check that particle name consistent with charge information.
1642
string particleName = particlePtr->name(1);
1643
if (particleName.size() > 16) {
1644
os << " Warning: particle " << idNow << " has name " << particleName
1645
<< " of length " << particleName.size() << "\n";
1651
for (int i = 0; i < int(particleName.size()); ++i) {
1652
if (particleName[i] == '+') ++nPos;
1653
if (particleName[i] == '-') ++nNeg;
1655
if ( (nPos > 0 && nNeg > 0) || ( nPos + nNeg > 0
1656
&& 3 * (nPos - nNeg) != chargeTypeNow )) {
1657
os << " Warning: particle " << idNow << " has name " << particleName
1658
<< " inconsistent with charge type " << chargeTypeNow << "\n";
1663
// Check that antiparticle name consistent with charge information.
1665
particleName = particlePtr->name(-1);
1666
if (particleName.size() > 16) {
1667
os << " Warning: particle " << idNow << " has name " << particleName
1668
<< " of length " << particleName.size() << "\n";
1674
for (int i = 0; i < int(particleName.size()); ++i) {
1675
if (particleName[i] == '+') ++nPos;
1676
if (particleName[i] == '-') ++nNeg;
1678
if ( (nPos > 0 && nNeg > 0) || ( nPos + nNeg > 0
1679
&& 3 * (nPos - nNeg) != -chargeTypeNow )) {
1680
os << " Warning: particle " << -idNow << " has name "
1681
<< particleName << " inconsistent with charge type "
1682
<< -chargeTypeNow << "\n";
1688
// Check that mass, mass range and width are consistent.
1689
if (particlePtr->useBreitWigner()) {
1690
if (mMinNow > m0Now) {
1691
os << " Error: particle " << idNow << " has mMin "
1692
<< fixed << setprecision(5) << mMinNow
1693
<< " larger than m0 " << m0Now << "\n";
1697
if (mMaxNow > mMinNow && mMaxNow < m0Now) {
1698
os << " Error: particle " << idNow << " has mMax "
1699
<< fixed << setprecision(5) << mMaxNow
1700
<< " smaller than m0 " << m0Now << "\n";
1704
if (mMaxNow > mMinNow && mMaxNow - mMinNow < mWidthNow) {
1705
os << " Warning: particle " << idNow << " has mMax - mMin "
1706
<< fixed << setprecision(5) << mMaxNow - mMinNow
1707
<< " smaller than mWidth " << mWidthNow << "\n";
1713
// Check that particle does not both have width and lifetime.
1714
if (mWidthNow > 0. && tau0Now > 0.) {
1715
os << " Warning: particle " << idNow << " has both nonvanishing width "
1716
<< scientific << setprecision(5) << mWidthNow << " and lifetime "
1722
// Begin study decay channels.
1723
if (particlePtr->sizeChannels() > 0) {
1725
// Loop through all decay channels.
1726
double bRatioSum = 0.;
1727
double bRatioPos = 0.;
1728
double bRatioNeg = 0.;
1729
bool hasPosBR = false;
1730
bool hasNegBR = false;
1731
double threshMass = 0.;
1732
bool openChannel = false;
1733
for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1735
// Extract channel properties.
1736
int onMode = particlePtr->channel(i).onMode();
1737
double bRatio = particlePtr->channel(i).bRatio();
1738
int meMode = particlePtr->channel(i).meMode();
1739
int mult = particlePtr->channel(i).multiplicity();
1741
for (int j = 0; j < 8; ++j)
1742
prod[j] = particlePtr->channel(i).product(j);
1744
// Sum branching ratios. Include off-channels.
1745
if (onMode == 0 || onMode == 1) bRatioSum += bRatio;
1746
else if (onMode == 2) {bRatioPos += bRatio; hasPosBR = true;}
1747
else if (onMode == 3) {bRatioNeg += bRatio; hasNegBR = true;}
1749
// Error printout when unknown decay product code.
1750
for (int j = 0; j < 8; ++j) {
1751
if ( prod[j] != 0 && !isParticle(prod[j]) ) {
1752
os << " Error: unknown decay product for " << idNow
1753
<< " -> " << prod[j] << "\n";
1760
// Error printout when no decay products or 0 interspersed.
1762
for (int j = 0; j < 8; ++j)
1763
if (prod[j] != 0) nLast = j + 1;
1764
if (mult == 0 || mult != nLast) {
1765
os << " Error: corrupted decay product list for "
1766
<< particlePtr->id() << " -> ";
1767
for (int j = 0; j < 8; ++j) os << prod[j] << " ";
1774
// Check charge conservation and open phase space in decay channel.
1775
int chargeTypeSum = -chargeTypeNow;
1776
int baryonTypeSum = -baryonTypeNow;
1777
double avgFinalMass = 0.;
1778
double minFinalMass = 0.;
1779
bool canHandle = true;
1780
for (int j = 0; j < mult; ++j) {
1781
chargeTypeSum += chargeType( prod[j] );
1782
baryonTypeSum += baryonNumberType( prod[j] );
1783
avgFinalMass += m0( prod[j] );
1784
minFinalMass += m0Min( prod[j] );
1785
if (prod[j] == 81 || prod[j] == 82 || prod[j] == 83)
1788
threshMass += bRatio * avgFinalMass;
1790
// Error printout when charge or baryon number not conserved.
1791
if (chargeTypeSum != 0 && canHandle) {
1792
os << " Error: 3*charge changed by " << chargeTypeSum
1793
<< " in " << idNow << " -> ";
1794
for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1800
if ( baryonTypeSum != 0 && canHandle && particlePtr->isHadron() ) {
1801
os << " Error: 3*baryon number changed by " << baryonTypeSum
1802
<< " in " << idNow << " -> ";
1803
for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1810
// Begin check that some matrix elements are used correctly.
1811
bool correctME = true;
1813
// Check matrix element mode 0: recommended not into partons.
1814
if (meMode == 0 && !isResonanceNow) {
1815
bool hasPartons = false;
1816
for (int j = 0; j < mult; ++j) {
1817
int idAbs = abs(prod[j]);
1818
if ( idAbs < 10 || idAbs == 21 || idAbs == 81 || idAbs == 82
1819
|| idAbs == 83 || (idAbs > 1000 && idAbs < 10000
1820
&& (idAbs/10)%10 == 0) ) hasPartons = true;
1822
if (hasPartons) correctME = false;
1825
// Check matrix element mode 1: rho/omega -> pi+ pi- pi0.
1826
bool useME1 = ( mult == 3 && spinTypeNow == 3 && idNow > 100
1828
&& particlePtr->channel(i).contains(211, -211, 111) );
1829
if ( meMode == 1 && !useME1 ) correctME = false;
1830
if ( meMode != 1 && useME1 ) correctME = false;
1832
// Check matrix element mode 2: polarization in V -> PS + PS.
1833
bool useME2 = ( mult == 2 && spinTypeNow == 3 && idNow > 100
1834
&& idNow < 1000 && spinType(prod[0]) == 1
1835
&& spinType(prod[1]) == 1 );
1836
if ( meMode == 2 && !useME2 ) correctME = false;
1837
if ( meMode != 2 && useME2 ) correctME = false;
1839
// Check matrix element mode 11, 12 and 13: Dalitz decay with
1840
// one or more particles in addition to the lepton pair,
1841
// or double Dalitz decay.
1842
bool useME11 = ( mult == 3 && !isResonanceNow
1843
&& (prod[1] == 11 || prod[1] == 13 || prod[1] == 15)
1844
&& prod[2] == -prod[1] );
1845
bool useME12 = ( mult > 3 && !isResonanceNow
1846
&& (prod[mult - 2] == 11 || prod[mult - 2] == 13
1847
|| prod[mult - 2] == 15) && prod[mult - 1] == -prod[mult - 2] );
1848
bool useME13 = ( mult == 4 && !isResonanceNow
1849
&& (prod[0] == 11 || prod[0] == 13) && prod[1] == -prod[0]
1850
&& (prod[2] == 11 || prod[2] == 13) && prod[3] == -prod[2] );
1851
if (useME13) useME12 = false;
1852
if ( meMode == 11 && !useME11 ) correctME = false;
1853
if ( meMode != 11 && useME11 ) correctME = false;
1854
if ( meMode == 12 && !useME12 ) correctME = false;
1855
if ( meMode != 12 && useME12 ) correctME = false;
1856
if ( meMode == 13 && !useME13 ) correctME = false;
1857
if ( meMode != 13 && useME13 ) correctME = false;
1859
// Check matrix element mode 21: tau -> nu_tau hadrons.
1860
bool useME21 = (idNow == 15 && mult > 2 && prod[0] == 16
1861
&& abs(prod[1]) > 100);
1862
if ( meMode == 21 && !useME21 ) correctME = false;
1863
if ( meMode != 21 && useME21 ) correctME = false;
1865
// Check matrix element mode 22, but only for semileptonic decay.
1866
// For a -> b c d require types u = 2, ubar = -2, d = 1, dbar = -1.
1867
if ( isLepton(prod[0]) && isLepton(prod[1]) ) {
1868
bool useME22 = false;
1872
if (particlePtr->isLepton()) {
1873
typeA = (idNow > 0) ? 1 + (idNow-1)%2 : -1 - (1-idNow)%2;
1874
} else if (particlePtr->isHadron()) {
1875
int hQ = particlePtr->heaviestQuark();
1876
// Special case: for B_c either bbar or c decays.
1877
if (idNow == 541 && heaviestQuark(prod[2]) == -5) hQ = 4;
1878
typeA = (hQ > 0) ? 1 + (hQ-1)%2 : -1 - (1-hQ)%2;
1880
typeB = (prod[0] > 0) ? 1 + (prod[0]-1)%2 : -1 - (1-prod[0])%2;
1881
typeC = (prod[1] > 0) ? 1 + (prod[1]-1)%2 : -1 - (1-prod[1])%2;
1883
if ( (idNow == 130 || idNow == 310) && typeC * typeA < 0)
1885
if (mult == 3 && idNow == 2112 && prod[2] == 2212)
1887
if (mult == 3 && idNow == 3222 && prod[2] == 3122)
1889
if (mult > 2 && typeC == typeA && typeB * typeC < 0
1890
&& (typeB + typeC)%2 != 0) useME22 = true;
1891
if ( meMode == 22 && !useME22 ) correctME = false;
1892
if ( meMode != 22 && useME22 ) correctME = false;
1895
// Check for matrix element mode 31.
1898
for (int j = 0; j < mult; ++j) if (prod[j] == 22) ++nGamma;
1899
if (nGamma != 1) correctME = false;
1902
// Check for unknown mode, or resonance-only mode.
1903
if ( !isResonanceNow && (meMode < 0 || (meMode > 2 && meMode <= 10)
1904
|| (meMode > 13 && meMode <= 20) || (meMode > 23 && meMode <= 30)
1905
|| (meMode > 31 && meMode <= 41) || meMode == 51 || meMode == 61
1906
|| meMode == 71 || (meMode > 80 && meMode <= 90)
1907
|| (!particlePtr->isOctetHadron() && meMode > 92) ) )
1910
// Print if incorrect matrix element mode.
1912
os << " Warning: meMode " << meMode << " used for "
1914
for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1920
// Warning printout when no phase space for decay.
1921
if ( studyCloser && verbosity > 0 && canHandle && onMode > 0
1922
&& particlePtr->m0Min() - minFinalMass < 0. ) {
1923
if (particlePtr->m0Max() - minFinalMass < 0.)
1924
os << " Error: decay never possible for ";
1925
else os << " Warning: decay sometimes not possible for ";
1926
os << idNow << " -> ";
1927
for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1933
// End loop over decay channels.
1934
if (onMode > 0 && bRatio > 0.) openChannel = true;
1937
// Optional printout of threshold.
1938
if (verbosity%10 > 1 && particlePtr->useBreitWigner()) {
1939
threshMass /= bRatioSum;
1940
os << " Info: particle " << idNow << fixed << setprecision(5)
1941
<< " has average mass threshold " << threshMass
1942
<< " while mMin is " << mMinNow << "\n";
1946
// Error printout when no acceptable decay channels found.
1947
if (studyCloser && !openChannel) {
1948
os << " Error: no acceptable decay channel found for particle "
1954
// Warning printout when branching ratios do not sum to unity.
1955
if (studyCloser && (!hasAntiNow || (!hasPosBR && !hasNegBR))
1956
&& abs(bRatioSum + bRatioPos - 1.) > 1e-8) {
1957
os << " Warning: particle " << idNow << fixed << setprecision(8)
1958
<< " has branching ratio sum " << bRatioSum << "\n";
1961
} else if (studyCloser && hasAntiNow
1962
&& (abs(bRatioSum + bRatioPos - 1.) > 1e-8
1963
|| abs(bRatioSum + bRatioNeg - 1.) > 1e-8)) {
1964
os << " Warning: particle " << idNow << fixed << setprecision(8)
1965
<< " has branching ratio sum " << bRatioSum + bRatioPos
1966
<< " while its antiparticle has " << bRatioSum + bRatioNeg
1972
// End study of decay channels and loop over particles.
1974
if (hasPrinted) os << "\n";
1977
// Final output. Done.
1978
os << " Total number of errors and warnings is " << nErr << "\n";
1979
os << "\n -------- End PYTHIA Check of Particle Data Table --------"
1980
<< "------\n" << endl;
1984
//--------------------------------------------------------------------------
1986
// Return the id of the sequentially next particle stored in table.
1988
int ParticleData::nextId(int idIn) {
1990
// Return 0 for negative or unknown codes. Return first for 0.
1991
if (idIn < 0 || (idIn > 0 && !isParticle(idIn))) return 0;
1992
if (idIn == 0) return pdt.begin()->first;
1994
// Find pointer to current particle and step up. Return 0 if impossible.
1995
map<int, ParticleDataEntry>::const_iterator pdtIn = pdt.find(idIn);
1996
if (pdtIn == pdt.end()) return 0;
1997
return (++pdtIn)->first;
2001
//--------------------------------------------------------------------------
2003
// Fractional width associated with open channels of one or two resonances.
2005
double ParticleData::resOpenFrac(int id1In, int id2In, int id3In) {
2011
if (isParticle(id1In)) answer = pdt[abs(id1In)].resOpenFrac(id1In);
2013
// Possibly second resonance.
2014
if (isParticle(id2In)) answer *= pdt[abs(id2In)].resOpenFrac(id2In);
2016
// Possibly third resonance.
2017
if (isParticle(id3In)) answer *= pdt[abs(id2In)].resOpenFrac(id3In);
2024
//--------------------------------------------------------------------------
2026
// Convert string to lowercase for case-insensitive comparisons.
2028
string ParticleData::toLower(const string& nameConv) {
2030
string temp(nameConv);
2031
for (int i = 0; i < int(temp.length()); ++i)
2032
temp[i] = std::tolower(temp[i]);
2037
//--------------------------------------------------------------------------
2039
// Allow several alternative inputs for true/false.
2041
bool ParticleData::boolString(string tag) {
2043
string tagLow = toLower(tag);
2044
return ( tagLow == "true" || tagLow == "1" || tagLow == "on"
2045
|| tagLow == "yes" || tagLow == "ok" );
2049
//--------------------------------------------------------------------------
2051
// Extract XML value string following XML attribute.
2053
string ParticleData::attributeValue(string line, string attribute) {
2055
if (line.find(attribute) == string::npos) return "";
2056
int iBegAttri = line.find(attribute);
2057
int iBegQuote = line.find("\"", iBegAttri + 1);
2058
int iEndQuote = line.find("\"", iBegQuote + 1);
2059
return line.substr(iBegQuote + 1, iEndQuote - iBegQuote - 1);
2063
//--------------------------------------------------------------------------
2065
// Extract XML bool value following XML attribute.
2067
bool ParticleData::boolAttributeValue(string line, string attribute) {
2069
string valString = attributeValue(line, attribute);
2070
if (valString == "") return false;
2071
return boolString(valString);
2075
//--------------------------------------------------------------------------
2077
// Extract XML int value following XML attribute.
2079
int ParticleData::intAttributeValue(string line, string attribute) {
2080
string valString = attributeValue(line, attribute);
2081
if (valString == "") return 0;
2082
istringstream valStream(valString);
2084
valStream >> intVal;
2089
//--------------------------------------------------------------------------
2091
// Extract XML double value following XML attribute.
2093
double ParticleData::doubleAttributeValue(string line, string attribute) {
2094
string valString = attributeValue(line, attribute);
2095
if (valString == "") return 0.;
2096
istringstream valStream(valString);
2098
valStream >> doubleVal;
2103
//==========================================================================
2105
} // end namespace Pythia8
added
removed
src/ParticleData.cc
