Line data Source code
1 : // ParticleData.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 :
6 : // Function definitions (not found in the header) for the
7 : // DecayChannel, ParticleDataEntry and ParticleData classes.
8 :
9 : #include "Pythia8/ParticleData.h"
10 : #include "Pythia8/ResonanceWidths.h"
11 : #include "Pythia8/StandardModel.h"
12 : #include "Pythia8/SusyResonanceWidths.h"
13 :
14 : // Allow string and character manipulation.
15 : #include <cctype>
16 :
17 : namespace Pythia8 {
18 :
19 : //==========================================================================
20 :
21 : // DecayChannel class.
22 : // This class holds info on a single decay channel.
23 :
24 : //--------------------------------------------------------------------------
25 :
26 : // Check whether id1 occurs anywhere in product list.
27 :
28 : bool DecayChannel::contains(int id1) const {
29 :
30 : bool found1 = false;
31 0 : for (int i = 0; i < nProd; ++i) if (prod[i] == id1) found1 = true;
32 0 : return found1;
33 :
34 : }
35 :
36 : //--------------------------------------------------------------------------
37 :
38 : // Check whether id1 and id2 occur anywhere in product list.
39 : // iF ID1 == ID2 then two copies of this particle must be present.
40 :
41 : bool DecayChannel::contains(int id1, int id2) const {
42 :
43 : bool found1 = false;
44 : bool found2 = false;
45 0 : for (int i = 0; i < nProd; ++i) {
46 0 : if (!found1 && prod[i] == id1) {found1 = true; continue;}
47 0 : if (!found2 && prod[i] == id2) {found2 = true; continue;}
48 : }
49 0 : return found1 && found2;
50 :
51 : }
52 :
53 : //--------------------------------------------------------------------------
54 :
55 : // Check whether id1, id2 and id3 occur anywhere in product list.
56 : // iF ID1 == ID2 then two copies of this particle must be present, etc.
57 :
58 : bool DecayChannel::contains(int id1, int id2, int id3) const {
59 :
60 : bool found1 = false;
61 : bool found2 = false;
62 : bool found3 = false;
63 0 : for (int i = 0; i < nProd; ++i) {
64 0 : if (!found1 && prod[i] == id1) {found1 = true; continue;}
65 0 : if (!found2 && prod[i] == id2) {found2 = true; continue;}
66 0 : if (!found3 && prod[i] == id3) {found3 = true; continue;}
67 : }
68 0 : return found1 && found2 && found3;
69 :
70 : }
71 :
72 : //==========================================================================
73 :
74 : // ParticleDataEntry class.
75 : // This class holds info on a single particle species.
76 :
77 : //--------------------------------------------------------------------------
78 :
79 : // Constants: could be changed here if desired, but normally should not.
80 : // These are of technical nature, as described for each.
81 :
82 : // A particle is invisible if it has neither strong nor electric charge,
83 : // and is not made up of constituents that have it. Only relevant for
84 : // long-lived particles. This list may need to be extended.
85 : const int ParticleDataEntry::INVISIBLENUMBER = 52;
86 : const int ParticleDataEntry::INVISIBLETABLE[52] = {
87 : 12, 14, 16, 18, 23, 25, 32, 33,
88 : 35, 36, 39, 41, 45, 46, 1000012, 1000014,
89 : 1000016, 1000018, 1000022, 1000023, 1000025, 1000035, 1000045, 1000039,
90 : 2000012, 2000014, 2000016, 2000018, 4900012, 4900014, 4900016, 4900021,
91 : 4900022, 4900101, 4900102, 4900103, 4900104, 4900105, 4900106, 4900107,
92 : 4900108, 4900111, 4900113, 4900211, 4900213, 4900991, 5000039, 5100039,
93 : 9900012, 9900014, 9900016, 9900023 };
94 :
95 : // For some particles we know it is necessary to switch off width,
96 : // although they do have one, so do not warn.
97 : const int ParticleDataEntry::KNOWNNOWIDTH[3] = {10313, 10323, 10333};
98 :
99 : // Particles with a read-in tau0 (in mm/c) below this mayDecay by default.
100 : const double ParticleDataEntry::MAXTAU0FORDECAY = 1000.;
101 :
102 : // Particles with a read-in m0 above this isResonance by default.
103 : const double ParticleDataEntry::MINMASSRESONANCE = 20.;
104 :
105 : // Narrow states are assigned nominal mass.
106 : const double ParticleDataEntry::NARROWMASS = 1e-6;
107 :
108 : // Constituent quark masses (d, u, s, c, b, -, -, -, g).
109 : const double ParticleDataEntry::CONSTITUENTMASSTABLE[10]
110 : = {0., 0.325, 0.325, 0.50, 1.60, 5.00, 0., 0., 0., 0.7};
111 :
112 : //--------------------------------------------------------------------------
113 :
114 : // Destructor: delete any ResonanceWidths object.
115 :
116 0 : ParticleDataEntry::~ParticleDataEntry() {
117 0 : if (resonancePtr != 0) delete resonancePtr;
118 0 : }
119 :
120 : //--------------------------------------------------------------------------
121 :
122 : // Set initial default values for some quantities.
123 :
124 : void ParticleDataEntry::setDefaults() {
125 :
126 : // A particle is a resonance if it is heavy enough.
127 0 : isResonanceSave = (m0Save > MINMASSRESONANCE);
128 :
129 : // A particle may decay if it is shortlived enough.
130 0 : mayDecaySave = (tau0Save < MAXTAU0FORDECAY);
131 :
132 : // A particle by default has no external decays.
133 0 : doExternalDecaySave = false;
134 :
135 : // A particle is invisible if in current table of such.
136 0 : isVisibleSave = true;
137 0 : for (int i = 0; i < INVISIBLENUMBER; ++i)
138 0 : if (idSave == INVISIBLETABLE[i]) isVisibleSave = false;
139 :
140 : // Normally a resonance should not have width forced to fixed value.
141 0 : doForceWidthSave = false;
142 :
143 : // Set up constituent masses.
144 0 : setConstituentMass();
145 :
146 : // No Breit-Wigner mass selection before initialized.
147 0 : modeBWnow = 0;
148 :
149 0 : }
150 :
151 : //--------------------------------------------------------------------------
152 :
153 : // Find out if a particle is a hadron.
154 : // Only covers normal hadrons, not e.g. R-hadrons.
155 :
156 : bool ParticleDataEntry::isHadron() const {
157 :
158 0 : if (idSave <= 100 || (idSave >= 1000000 && idSave <= 9000000)
159 0 : || idSave >= 9900000) return false;
160 0 : if (idSave == 130 || idSave == 310) return true;
161 0 : if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0)
162 0 : return false;
163 0 : return true;
164 :
165 0 : }
166 :
167 : //--------------------------------------------------------------------------
168 :
169 : // Find out if a particle is a meson.
170 : // Only covers normal hadrons, not e.g. R-hadrons.
171 :
172 : bool ParticleDataEntry::isMeson() const {
173 :
174 0 : if (idSave <= 100 || (idSave >= 1000000 && idSave <= 9000000)
175 0 : || idSave >= 9900000) return false;
176 0 : if (idSave == 130 || idSave == 310) return true;
177 0 : if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0
178 0 : || (idSave/1000)%10 != 0) return false;
179 0 : return true;
180 :
181 0 : }
182 :
183 : //--------------------------------------------------------------------------
184 :
185 : // Find out if a particle is a baryon.
186 : // Only covers normal hadrons, not e.g. R-hadrons.
187 :
188 : bool ParticleDataEntry::isBaryon() const {
189 :
190 0 : if (idSave <= 1000 || (idSave >= 1000000 && idSave <= 9000000)
191 0 : || idSave >= 9900000) return false;
192 0 : if (idSave%10 == 0 || (idSave/10)%10 == 0 || (idSave/100)%10 == 0
193 0 : || (idSave/1000)%10 == 0) return false;
194 0 : return true;
195 :
196 :
197 0 : }
198 :
199 : //--------------------------------------------------------------------------
200 :
201 : // Extract the heaviest (= largest id) quark in a hadron.
202 :
203 : int ParticleDataEntry::heaviestQuark(int idIn) const {
204 :
205 0 : if (!isHadron()) return 0;
206 : int hQ = 0;
207 :
208 : // Meson.
209 0 : if ( (idSave/1000)%10 == 0 ) {
210 0 : hQ = (idSave/100)%10;
211 0 : if (idSave == 130) hQ = 3;
212 0 : if (hQ%2 == 1) hQ = -hQ;
213 :
214 : // Baryon.
215 : } else hQ = (idSave/1000)%10;
216 :
217 : // Done.
218 0 : return (idIn > 0) ? hQ : -hQ;
219 :
220 0 : }
221 :
222 : //--------------------------------------------------------------------------
223 :
224 : // Calculate three times baryon number, i.e. net quark - antiquark number.
225 :
226 : int ParticleDataEntry::baryonNumberType(int idIn) const {
227 :
228 : // Quarks.
229 0 : if (isQuark()) return (idIn > 0) ? 1 : -1;
230 :
231 : // Diquarks
232 0 : if (isDiquark()) return (idIn > 0) ? 2 : -2;
233 :
234 : // Baryons.
235 0 : if (isBaryon()) return (idIn > 0) ? 3 : -3;
236 :
237 : // Done.
238 0 : return 0;
239 :
240 0 : }
241 :
242 : //--------------------------------------------------------------------------
243 :
244 : // Prepare the Breit-Wigner mass selection by precalculating
245 : // frequently-used expressions.
246 :
247 : void ParticleDataEntry::initBWmass() {
248 :
249 : // Find Breit-Wigner mode for current particle.
250 0 : modeBWnow = particleDataPtr->modeBreitWigner;
251 0 : if ( m0Save < NARROWMASS ) mWidthSave = 0.;
252 0 : if ( mWidthSave < NARROWMASS || (mMaxSave > mMinSave
253 0 : && mMaxSave - mMinSave < NARROWMASS) ) modeBWnow = 0;
254 0 : if (modeBWnow == 0) return;
255 :
256 : // Find atan expressions to be used in random mass selection.
257 0 : if (modeBWnow < 3) {
258 0 : atanLow = atan( 2. * (mMinSave - m0Save) / mWidthSave );
259 0 : double atanHigh = (mMaxSave > mMinSave)
260 0 : ? atan( 2. * (mMaxSave - m0Save) / mWidthSave ) : 0.5 * M_PI;
261 0 : atanDif = atanHigh - atanLow;
262 0 : } else {
263 0 : atanLow = atan( (pow2(mMinSave) - pow2(m0Save))
264 0 : / (m0Save * mWidthSave) );
265 0 : double atanHigh = (mMaxSave > mMinSave)
266 0 : ? atan( (pow2(mMaxSave) - pow2(m0Save)) / (m0Save * mWidthSave) )
267 : : 0.5 * M_PI;
268 0 : atanDif = atanHigh - atanLow;
269 : }
270 :
271 : // Done if no threshold factor.
272 0 : if (modeBWnow%2 == 1) return;
273 :
274 : // Find average mass threshold for threshold-factor correction.
275 : double bRatSum = 0.;
276 : double mThrSum = 0;
277 0 : for (int i = 0; i < int(channels.size()); ++i)
278 0 : if (channels[i].onMode() > 0) {
279 0 : bRatSum += channels[i].bRatio();
280 : double mChannelSum = 0.;
281 0 : for (int j = 0; j < channels[i].multiplicity(); ++j)
282 0 : mChannelSum += particleDataPtr->m0( channels[i].product(j) );
283 0 : mThrSum += channels[i].bRatio() * mChannelSum;
284 0 : }
285 0 : mThr = (bRatSum == 0.) ? 0. : mThrSum / bRatSum;
286 :
287 : // Switch off Breit-Wigner if very close to threshold.
288 0 : if (mThr + NARROWMASS > m0Save) {
289 0 : modeBWnow = 0;
290 : bool knownProblem = false;
291 0 : for (int i = 0; i < 3; ++i) if (idSave == KNOWNNOWIDTH[i])
292 0 : knownProblem = true;
293 0 : if (!knownProblem) {
294 0 : ostringstream osWarn;
295 0 : osWarn << "for id = " << idSave;
296 0 : particleDataPtr->infoPtr->errorMsg("Warning in ParticleDataEntry::"
297 0 : "initBWmass: switching off width", osWarn.str(), true);
298 0 : }
299 0 : }
300 :
301 0 : }
302 :
303 : //--------------------------------------------------------------------------
304 :
305 : // Function to give mass of a particle, either at the nominal value
306 : // or picked according to a (linear or quadratic) Breit-Wigner.
307 :
308 : double ParticleDataEntry::mSel() {
309 :
310 : // Nominal value. (Width check should not be needed, but just in case.)
311 0 : if (modeBWnow == 0 || mWidthSave < NARROWMASS) return m0Save;
312 0 : double mNow, m2Now;
313 :
314 : // Mass according to a Breit-Wigner linear in m.
315 0 : if (modeBWnow == 1) {
316 0 : mNow = m0Save + 0.5 * mWidthSave
317 0 : * tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
318 :
319 : // Ditto, but make Gamma proportional to sqrt(m^2 - m_threshold^2).
320 0 : } else if (modeBWnow == 2) {
321 : double mWidthNow, fixBW, runBW;
322 0 : double m0ThrS = m0Save*m0Save - mThr*mThr;
323 0 : do {
324 0 : mNow = m0Save + 0.5 * mWidthSave
325 0 : * tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
326 0 : mWidthNow = mWidthSave * sqrtpos( (mNow*mNow - mThr*mThr) / m0ThrS );
327 0 : fixBW = mWidthSave / (pow2(mNow - m0Save) + pow2(0.5 * mWidthSave));
328 0 : runBW = mWidthNow / (pow2(mNow - m0Save) + pow2(0.5 * mWidthNow));
329 0 : } while (runBW < particleDataPtr->rndmPtr->flat()
330 0 : * particleDataPtr->maxEnhanceBW * fixBW);
331 :
332 : // Mass according to a Breit-Wigner quadratic in m.
333 0 : } else if (modeBWnow == 3) {
334 0 : m2Now = m0Save*m0Save + m0Save * mWidthSave
335 0 : * tan( atanLow + atanDif * particleDataPtr->rndmPtr->flat() );
336 0 : mNow = sqrtpos( m2Now);
337 :
338 : // Ditto, but m_0 Gamma_0 -> m Gamma(m) with threshold factor as above.
339 0 : } else {
340 0 : double mwNow, fixBW, runBW;
341 0 : double m2Ref = m0Save * m0Save;
342 0 : double mwRef = m0Save * mWidthSave;
343 0 : double m2Thr = mThr * mThr;
344 0 : do {
345 0 : m2Now = m2Ref + mwRef * tan( atanLow + atanDif
346 0 : * particleDataPtr->rndmPtr->flat() );
347 0 : mNow = sqrtpos( m2Now);
348 0 : mwNow = mNow * mWidthSave
349 0 : * sqrtpos( (m2Now - m2Thr) / (m2Ref - m2Thr) );
350 0 : fixBW = mwRef / (pow2(m2Now - m2Ref) + pow2(mwRef));
351 0 : runBW = mwNow / (pow2(m2Now - m2Ref) + pow2(mwNow));
352 0 : } while (runBW < particleDataPtr->rndmPtr->flat()
353 0 : * particleDataPtr->maxEnhanceBW * fixBW);
354 0 : }
355 :
356 : // Done.
357 : return mNow;
358 0 : }
359 :
360 : //--------------------------------------------------------------------------
361 :
362 : // Function to calculate running mass at given mass scale.
363 :
364 : double ParticleDataEntry::mRun(double mHat) {
365 :
366 : // Except for six quarks return nominal mass.
367 0 : if (idSave > 6) return m0Save;
368 0 : double mQRun = particleDataPtr->mQRun[idSave];
369 0 : double Lam5 = particleDataPtr->Lambda5Run;
370 :
371 : // For d, u, s quarks start running at 2 GeV (RPP 2006 p. 505).
372 0 : if (idSave < 4) return mQRun * pow ( log(2. / Lam5)
373 0 : / log(max(2., mHat) / Lam5), 12./23.);
374 :
375 : // For c, b and t quarks start running at respective mass.
376 0 : return mQRun * pow ( log(mQRun / Lam5)
377 0 : / log(max(mQRun, mHat) / Lam5), 12./23.);
378 :
379 0 : }
380 :
381 : //--------------------------------------------------------------------------
382 :
383 : // Rescale all branching ratios to assure normalization to unity.
384 :
385 : void ParticleDataEntry::rescaleBR(double newSumBR) {
386 :
387 : // Sum up branching ratios. Find rescaling factor. Rescale.
388 : double oldSumBR = 0.;
389 0 : for ( int i = 0; i < int(channels.size()); ++ i)
390 0 : oldSumBR += channels[i].bRatio();
391 0 : double rescaleFactor = newSumBR / oldSumBR;
392 0 : for ( int i = 0; i < int(channels.size()); ++ i)
393 0 : channels[i].rescaleBR(rescaleFactor);
394 :
395 0 : }
396 :
397 : //--------------------------------------------------------------------------
398 :
399 : // Prepare to pick a decay channel.
400 :
401 : bool ParticleDataEntry::preparePick(int idSgn, double mHat, int idInFlav) {
402 :
403 : // Reset sum of allowed widths/branching ratios.
404 0 : currentBRSum = 0.;
405 :
406 : // For resonances the widths are calculated dynamically.
407 0 : if (isResonanceSave && resonancePtr != 0) {
408 0 : resonancePtr->widthStore(idSgn, mHat, idInFlav);
409 0 : for (int i = 0; i < int(channels.size()); ++i)
410 0 : currentBRSum += channels[i].currentBR();
411 :
412 : // Else use normal fixed branching ratios.
413 0 : } else {
414 : int onMode;
415 : double currentBRNow;
416 0 : for (int i = 0; i < int(channels.size()); ++i) {
417 0 : onMode = channels[i].onMode();
418 : currentBRNow = 0.;
419 0 : if ( idSgn > 0 && (onMode == 1 || onMode == 2) )
420 0 : currentBRNow = channels[i].bRatio();
421 0 : else if ( idSgn < 0 && (onMode == 1 || onMode == 3) )
422 0 : currentBRNow = channels[i].bRatio();
423 0 : channels[i].currentBR(currentBRNow);
424 0 : currentBRSum += currentBRNow;
425 : }
426 : }
427 :
428 : // Failure if no channels found with positive branching ratios.
429 0 : return (currentBRSum > 0.);
430 :
431 : }
432 :
433 : //--------------------------------------------------------------------------
434 :
435 : // Pick a decay channel according to branching ratios from preparePick.
436 :
437 : DecayChannel& ParticleDataEntry::pickChannel() {
438 :
439 : // Find channel in table.
440 0 : int size = channels.size();
441 0 : double rndmBR = currentBRSum * particleDataPtr->rndmPtr->flat();
442 : int i = -1;
443 0 : do rndmBR -= channels[++i].currentBR();
444 0 : while (rndmBR > 0. && i < size);
445 :
446 : // Emergency if no channel found. Done.
447 0 : if (i == size) i = 0;
448 0 : return channels[i];
449 :
450 : }
451 :
452 : //--------------------------------------------------------------------------
453 :
454 : // Access methods stored in ResonanceWidths. Could have been
455 : // inline in .h, except for problems with forward declarations.
456 :
457 : void ParticleDataEntry::setResonancePtr(
458 : ResonanceWidths* resonancePtrIn) {
459 0 : if (resonancePtr == resonancePtrIn) return;
460 0 : if (resonancePtr != 0) delete resonancePtr;
461 0 : resonancePtr = resonancePtrIn;
462 0 : }
463 :
464 : void ParticleDataEntry::resInit(Info* infoPtrIn, Settings* settingsPtrIn,
465 : ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
466 0 : if (resonancePtr != 0) resonancePtr->init(infoPtrIn, settingsPtrIn,
467 : particleDataPtrIn, couplingsPtrIn);
468 0 : }
469 :
470 : double ParticleDataEntry::resWidth(int idSgn, double mHat, int idIn,
471 : bool openOnly, bool setBR) {
472 0 : return (resonancePtr != 0) ? resonancePtr->width( idSgn, mHat,
473 0 : idIn, openOnly, setBR) : 0.;
474 : }
475 :
476 : double ParticleDataEntry::resWidthOpen(int idSgn, double mHat, int idIn) {
477 0 : return (resonancePtr != 0) ? resonancePtr->widthOpen( idSgn, mHat, idIn)
478 : : 0.;
479 : }
480 :
481 : double ParticleDataEntry::resWidthStore(int idSgn, double mHat, int idIn) {
482 0 : return (resonancePtr != 0) ? resonancePtr->widthStore( idSgn, mHat, idIn)
483 : : 0.;
484 : }
485 :
486 : double ParticleDataEntry::resOpenFrac(int idSgn) {
487 0 : return (resonancePtr != 0) ? resonancePtr->openFrac(idSgn) : 1.;
488 : }
489 :
490 : double ParticleDataEntry::resWidthRescaleFactor() {
491 0 : return (resonancePtr != 0) ? resonancePtr->widthRescaleFactor() : 1.;
492 : }
493 :
494 : double ParticleDataEntry::resWidthChan(double mHat, int idAbs1,
495 : int idAbs2) {
496 0 : return (resonancePtr != 0) ? resonancePtr->widthChan( mHat, idAbs1,
497 : idAbs2) : 0.;
498 : }
499 :
500 : //--------------------------------------------------------------------------
501 :
502 : // Constituent masses for (d, u, s, c, b) quarks and diquarks.
503 : // Hardcoded in CONSTITUENTMASSTABLE so that they are not overwritten
504 : // by mistake, and separated from the "normal" masses.
505 : // Called both by setDefaults and setM0 so kept as separate method.
506 :
507 : void ParticleDataEntry::setConstituentMass() {
508 :
509 : // Equate with the normal masses as default guess.
510 0 : constituentMassSave = m0Save;
511 :
512 : // Quark masses trivial. Also gluon mass.
513 0 : if (idSave < 6) constituentMassSave = CONSTITUENTMASSTABLE[idSave];
514 0 : if (idSave == 21) constituentMassSave = CONSTITUENTMASSTABLE[9];
515 :
516 : // Diquarks as simple sum of constituent quarks.
517 0 : if (idSave > 1000 && idSave < 10000 && (idSave/10)%10 == 0) {
518 0 : int id1 = idSave/1000;
519 0 : int id2 = (idSave/100)%10;
520 0 : if (id1 <6 && id2 < 6) constituentMassSave
521 0 : = CONSTITUENTMASSTABLE[id1] + CONSTITUENTMASSTABLE[id2];
522 0 : }
523 :
524 0 : }
525 :
526 : //==========================================================================
527 :
528 : // ParticleData class.
529 : // This class holds a map of all ParticleDataEntries,
530 : // each entry containing info on a particle species.
531 :
532 : //--------------------------------------------------------------------------
533 :
534 : // Get data to be distributed among particles during setup.
535 : // Note: this routine is called twice. Firstly from init(...), but
536 : // the data should not be used at that point, so is likely overkill.
537 : // Secondly, from initWidths, after user has had time to change.
538 :
539 : void ParticleData::initCommon() {
540 :
541 : // Mass generation: fixed mass or linear/quadratic Breit-Wigner.
542 0 : modeBreitWigner = settingsPtr->mode("ParticleData:modeBreitWigner");
543 :
544 : // Maximum tail enhancement when adding threshold factor to Breit-Wigner.
545 0 : maxEnhanceBW = settingsPtr->parm("ParticleData:maxEnhanceBW");
546 :
547 : // Find initial MSbar masses for five light flavours.
548 0 : mQRun[1] = settingsPtr->parm("ParticleData:mdRun");
549 0 : mQRun[2] = settingsPtr->parm("ParticleData:muRun");
550 0 : mQRun[3] = settingsPtr->parm("ParticleData:msRun");
551 0 : mQRun[4] = settingsPtr->parm("ParticleData:mcRun");
552 0 : mQRun[5] = settingsPtr->parm("ParticleData:mbRun");
553 0 : mQRun[6] = settingsPtr->parm("ParticleData:mtRun");
554 :
555 : // Find Lambda5 value to use in running of MSbar masses.
556 0 : double alphaSvalue = settingsPtr->parm("ParticleData:alphaSvalueMRun");
557 0 : AlphaStrong alphaS;
558 0 : alphaS.init( alphaSvalue, 1, 5, false);
559 0 : Lambda5Run = alphaS.Lambda5();
560 :
561 0 : }
562 :
563 : //--------------------------------------------------------------------------
564 :
565 : // Initialize pointer for particles to the full database, the Breit-Wigners
566 : // of normal hadrons and the ResonanceWidths of resonances. For the latter
567 : // the order of initialization is essential to get secondary widths right.
568 :
569 : void ParticleData::initWidths( vector<ResonanceWidths*> resonancePtrs) {
570 :
571 : // Initialize some common data.
572 0 : initCommon();
573 :
574 : // Pointer to database and Breit-Wigner mass initialization for each
575 : // particle.
576 : ResonanceWidths* resonancePtr = 0;
577 0 : for (map<int, ParticleDataEntry>::iterator pdtEntry = pdt.begin();
578 0 : pdtEntry != pdt.end(); ++pdtEntry) {
579 0 : ParticleDataEntry& pdtNow = pdtEntry->second;
580 0 : pdtNow.initBWmass();
581 :
582 : // Remove any existing resonances.
583 0 : resonancePtr = pdtNow.getResonancePtr();
584 0 : if (resonancePtr != 0) pdtNow.setResonancePtr(0);
585 : }
586 :
587 : // Begin set up new resonance objects.
588 : // Z0, W+- and top are almost always needed.
589 0 : resonancePtr = new ResonanceGmZ(23);
590 0 : setResonancePtr( 23, resonancePtr);
591 0 : resonancePtr = new ResonanceW(24);
592 0 : setResonancePtr( 24, resonancePtr);
593 0 : resonancePtr = new ResonanceTop(6);
594 0 : setResonancePtr( 6, resonancePtr);
595 :
596 : // Higgs in SM.
597 0 : if (!settingsPtr->flag("Higgs:useBSM")) {
598 0 : resonancePtr = new ResonanceH(0, 25);
599 0 : setResonancePtr( 25, resonancePtr);
600 :
601 : // Higgses in BSM.
602 0 : } else {
603 0 : resonancePtr = new ResonanceH(1, 25);
604 0 : setResonancePtr( 25, resonancePtr);
605 0 : resonancePtr = new ResonanceH(2, 35);
606 0 : setResonancePtr( 35, resonancePtr);
607 0 : resonancePtr = new ResonanceH(3, 36);
608 0 : setResonancePtr( 36, resonancePtr);
609 0 : resonancePtr = new ResonanceHchg(37);
610 0 : setResonancePtr( 37, resonancePtr);
611 0 : resonancePtr = new ResonanceH(4, 45);
612 0 : setResonancePtr( 45, resonancePtr);
613 0 : resonancePtr = new ResonanceH(5, 46);
614 0 : setResonancePtr( 46, resonancePtr);
615 : }
616 :
617 : // A fourth generation: b', t', tau', nu'_tau.
618 0 : resonancePtr = new ResonanceFour(7);
619 0 : setResonancePtr( 7, resonancePtr);
620 0 : resonancePtr = new ResonanceFour(8);
621 0 : setResonancePtr( 8, resonancePtr);
622 0 : resonancePtr = new ResonanceFour(17);
623 0 : setResonancePtr( 17, resonancePtr);
624 0 : resonancePtr = new ResonanceFour(18);
625 0 : setResonancePtr( 18, resonancePtr);
626 :
627 : // New gauge bosons: Z', W', R.
628 0 : resonancePtr = new ResonanceZprime(32);
629 0 : setResonancePtr( 32, resonancePtr);
630 0 : resonancePtr = new ResonanceWprime(34);
631 0 : setResonancePtr( 34, resonancePtr);
632 0 : resonancePtr = new ResonanceRhorizontal(41);
633 0 : setResonancePtr( 41, resonancePtr);
634 :
635 : // A leptoquark.
636 0 : resonancePtr = new ResonanceLeptoquark(42);
637 0 : setResonancePtr( 42, resonancePtr);
638 :
639 : // 93 = Z0copy and 94 = W+-copy used to pick decay channels
640 : // for W/Z production in parton showers.
641 0 : resonancePtr = new ResonanceGmZ(93);
642 0 : setResonancePtr( 93, resonancePtr);
643 0 : resonancePtr = new ResonanceW(94);
644 0 : setResonancePtr( 94, resonancePtr);
645 :
646 : // Supersymmetry
647 : // - Squarks
648 0 : for(int i = 1; i < 7; i++){
649 0 : resonancePtr = new ResonanceSquark(1000000 + i);
650 0 : setResonancePtr( 1000000 + i, resonancePtr);
651 0 : resonancePtr = new ResonanceSquark(2000000 + i);
652 0 : setResonancePtr( 2000000 + i, resonancePtr);
653 : }
654 :
655 : // - Sleptons and sneutrinos
656 0 : for(int i = 1; i < 7; i++){
657 0 : resonancePtr = new ResonanceSlepton(1000010 + i);
658 0 : setResonancePtr( 1000010 + i, resonancePtr);
659 0 : resonancePtr = new ResonanceSlepton(2000010 + i);
660 0 : setResonancePtr( 2000010 + i, resonancePtr);
661 : }
662 :
663 : // - Gluino
664 0 : resonancePtr = new ResonanceGluino(1000021);
665 0 : setResonancePtr( 1000021, resonancePtr);
666 :
667 : // - Charginos
668 0 : resonancePtr = new ResonanceChar(1000024);
669 0 : setResonancePtr( 1000024, resonancePtr);
670 0 : resonancePtr = new ResonanceChar(1000037);
671 0 : setResonancePtr( 1000037, resonancePtr);
672 :
673 : // - Neutralinos
674 0 : resonancePtr = new ResonanceNeut(1000022);
675 0 : setResonancePtr( 1000022, resonancePtr);
676 0 : resonancePtr = new ResonanceNeut(1000023);
677 0 : setResonancePtr( 1000023, resonancePtr);
678 0 : resonancePtr = new ResonanceNeut(1000025);
679 0 : setResonancePtr( 1000025, resonancePtr);
680 0 : resonancePtr = new ResonanceNeut(1000035);
681 0 : setResonancePtr( 1000035, resonancePtr);
682 0 : resonancePtr = new ResonanceNeut(1000045);
683 0 : setResonancePtr( 1000045, resonancePtr);
684 :
685 : // Excited quarks and leptons.
686 0 : for (int i = 1; i < 7; ++i) {
687 0 : resonancePtr = new ResonanceExcited(4000000 + i);
688 0 : setResonancePtr( 4000000 + i, resonancePtr);
689 : }
690 0 : for (int i = 11; i < 17; ++i) {
691 0 : resonancePtr = new ResonanceExcited(4000000 + i);
692 0 : setResonancePtr( 4000000 + i, resonancePtr);
693 : }
694 :
695 : // An excited graviton/gluon in extra-dimensional scenarios.
696 0 : resonancePtr = new ResonanceGraviton(5100039);
697 0 : setResonancePtr( 5100039, resonancePtr);
698 0 : resonancePtr = new ResonanceKKgluon(5100021);
699 0 : setResonancePtr( 5100021, resonancePtr);
700 :
701 : // A left-right-symmetric scenario with new righthanded neutrinos,
702 : // righthanded gauge bosons and doubly charged Higgses.
703 0 : resonancePtr = new ResonanceNuRight(9900012);
704 0 : setResonancePtr( 9900012, resonancePtr);
705 0 : resonancePtr = new ResonanceNuRight(9900014);
706 0 : setResonancePtr( 9900014, resonancePtr);
707 0 : resonancePtr = new ResonanceNuRight(9900016);
708 0 : setResonancePtr( 9900016, resonancePtr);
709 0 : resonancePtr = new ResonanceZRight(9900023);
710 0 : setResonancePtr( 9900023, resonancePtr);
711 0 : resonancePtr = new ResonanceWRight(9900024);
712 0 : setResonancePtr( 9900024, resonancePtr);
713 0 : resonancePtr = new ResonanceHchgchgLeft(9900041);
714 0 : setResonancePtr( 9900041, resonancePtr);
715 0 : resonancePtr = new ResonanceHchgchgRight(9900042);
716 0 : setResonancePtr( 9900042, resonancePtr);
717 :
718 : // Attach user-defined external resonances and do basic initialization.
719 0 : for (int i = 0; i < int(resonancePtrs.size()); ++i) {
720 0 : int idNow = resonancePtrs[i]->id();
721 0 : resonancePtrs[i]->initBasic(idNow);
722 0 : setResonancePtr( idNow, resonancePtrs[i]);
723 : }
724 :
725 : // Set up lists to order resonances in ascending mass.
726 0 : vector<int> idOrdered;
727 0 : vector<double> m0Ordered;
728 :
729 : // Put Z0 and W+- first, since known to be SM and often off-shell.
730 0 : idOrdered.push_back(23);
731 0 : m0Ordered.push_back(m0(23));
732 0 : idOrdered.push_back(24);
733 0 : m0Ordered.push_back(m0(24));
734 :
735 : // Loop through particle table to find resonances.
736 0 : for (map<int, ParticleDataEntry>::iterator pdtEntry = pdt.begin();
737 0 : pdtEntry != pdt.end(); ++pdtEntry) {
738 0 : ParticleDataEntry& pdtNow = pdtEntry->second;
739 0 : int idNow = pdtNow.id();
740 :
741 : // Set up a simple default object for uninitialized resonances.
742 0 : if (pdtNow.isResonance() && pdtNow.getResonancePtr() == 0) {
743 0 : resonancePtr = new ResonanceGeneric(idNow);
744 0 : setResonancePtr( idNow, resonancePtr);
745 : }
746 :
747 : // Insert resonances in ascending mass, to respect decay hierarchies.
748 0 : if (pdtNow.getResonancePtr() != 0 && idNow != 23 && idNow != 24) {
749 0 : double m0Now = pdtNow.m0();
750 0 : idOrdered.push_back(idNow);
751 0 : m0Ordered.push_back(m0Now);
752 0 : for (int i = int(idOrdered.size()) - 2; i > 1; --i) {
753 0 : if (m0Ordered[i] < m0Now) break;
754 0 : swap( idOrdered[i], idOrdered[i + 1]);
755 0 : swap( m0Ordered[i], m0Ordered[i + 1]);
756 : }
757 0 : }
758 0 : }
759 :
760 : // Initialize the resonances in ascending mass order. Reset mass generation.
761 0 : for (int i = 0; i < int(idOrdered.size()); ++i) {
762 0 : resInit( idOrdered[i]);
763 0 : ParticleDataEntry* pdtPtrNow = particleDataEntryPtr( idOrdered[i]);
764 0 : pdtPtrNow->initBWmass();
765 : }
766 :
767 0 : }
768 :
769 : //--------------------------------------------------------------------------
770 :
771 : // Read in database from specific XML file (which may refer to others).
772 :
773 : bool ParticleData::readXML(string inFile, bool reset) {
774 :
775 : // Normally reset whole database before beginning.
776 0 : if (reset) {pdt.clear(); isInit = false;}
777 :
778 : // List of files to be checked.
779 0 : vector<string> files;
780 0 : files.push_back(inFile);
781 :
782 : // Loop over files. Open them for read.
783 0 : for (int i = 0; i < int(files.size()); ++i) {
784 0 : const char* cstring = files[i].c_str();
785 0 : ifstream is(cstring);
786 :
787 : // Check that instream is OK.
788 0 : if (!is.good()) {
789 0 : infoPtr->errorMsg("Error in ParticleData::readXML:"
790 0 : " did not find file", files[i]);
791 0 : return false;
792 : }
793 :
794 : // Read in one line at a time.
795 0 : particlePtr = 0;
796 0 : string line;
797 0 : while ( getline(is, line) ) {
798 :
799 : // Get first word of a line.
800 0 : istringstream getfirst(line);
801 0 : string word1;
802 0 : getfirst >> word1;
803 :
804 : // Check for occurence of a particle. Add any continuation lines.
805 0 : if (word1 == "<particle") {
806 0 : while (line.find(">") == string::npos) {
807 0 : string addLine;
808 0 : getline(is, addLine);
809 0 : line += addLine;
810 0 : }
811 :
812 : // Read in particle properties.
813 0 : int idTmp = intAttributeValue( line, "id");
814 0 : string nameTmp = attributeValue( line, "name");
815 0 : string antiNameTmp = attributeValue( line, "antiName");
816 0 : if (antiNameTmp == "") antiNameTmp = "void";
817 0 : int spinTypeTmp = intAttributeValue( line, "spinType");
818 0 : int chargeTypeTmp = intAttributeValue( line, "chargeType");
819 0 : int colTypeTmp = intAttributeValue( line, "colType");
820 0 : double m0Tmp = doubleAttributeValue( line, "m0");
821 0 : double mWidthTmp = doubleAttributeValue( line, "mWidth");
822 0 : double mMinTmp = doubleAttributeValue( line, "mMin");
823 0 : double mMaxTmp = doubleAttributeValue( line, "mMax");
824 0 : double tau0Tmp = doubleAttributeValue( line, "tau0");
825 :
826 : // Erase if particle already exists.
827 0 : if (isParticle(idTmp)) pdt.erase(idTmp);
828 :
829 : // Store new particle. Save pointer, to be used for decay channels.
830 0 : addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
831 : colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
832 0 : particlePtr = particleDataEntryPtr(idTmp);
833 :
834 : // Check for occurence of a decay channel. Add any continuation lines.
835 0 : } else if (word1 == "<channel") {
836 0 : while (line.find(">") == string::npos) {
837 0 : string addLine;
838 0 : getline(is, addLine);
839 0 : line += addLine;
840 0 : }
841 :
842 : // Read in channel properties - products so far only as a string.
843 0 : int onMode = intAttributeValue( line, "onMode");
844 0 : double bRatio = doubleAttributeValue( line, "bRatio");
845 0 : int meMode = intAttributeValue( line, "meMode");
846 0 : string products = attributeValue( line, "products");
847 :
848 : // Read in decay products from stream. Must have at least one.
849 0 : istringstream prodStream(products);
850 0 : int prod0 = 0; int prod1 = 0; int prod2 = 0; int prod3 = 0;
851 0 : int prod4 = 0; int prod5 = 0; int prod6 = 0; int prod7 = 0;
852 0 : prodStream >> prod0 >> prod1 >> prod2 >> prod3 >> prod4 >> prod5
853 0 : >> prod6 >> prod7;
854 0 : if (prod0 == 0) {
855 0 : infoPtr->errorMsg("Error in ParticleData::readXML:"
856 0 : " incomplete decay channel", line);
857 0 : return false;
858 : }
859 :
860 : // Store new channel (if particle already known).
861 0 : if (particlePtr == 0) {
862 0 : infoPtr->errorMsg("Error in ParticleData::readXML:"
863 0 : " orphan decay channel", line);
864 0 : return false;
865 : }
866 0 : particlePtr->addChannel(onMode, bRatio, meMode, prod0, prod1,
867 0 : prod2, prod3, prod4, prod5, prod6, prod7);
868 :
869 : // Check for occurence of a file also to be read.
870 0 : } else if (word1 == "<file") {
871 0 : string file = attributeValue(line, "name");
872 0 : if (file == "") {
873 0 : infoPtr->errorMsg("Error in ParticleData::readXML:"
874 0 : " skip unrecognized file name", line);
875 0 : } else files.push_back(file);
876 0 : }
877 :
878 : // End of loop over lines in input file and loop over files.
879 0 : };
880 0 : };
881 :
882 : // All particle data at this stage defines baseline original.
883 0 : if (reset) for (map<int, ParticleDataEntry>::iterator pdtEntry
884 0 : = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
885 0 : particlePtr = &pdtEntry->second;
886 0 : particlePtr->setHasChanged(false);
887 0 : }
888 :
889 : // Done.
890 0 : isInit = true;
891 0 : return true;
892 :
893 0 : }
894 :
895 : //--------------------------------------------------------------------------
896 :
897 : // Print out complete database in numerical order as an XML file.
898 :
899 : void ParticleData::listXML(string outFile) {
900 :
901 : // Convert file name to ofstream.
902 0 : const char* cstring = outFile.c_str();
903 0 : ofstream os(cstring);
904 :
905 : // Iterate through the particle data table.
906 0 : for (map<int, ParticleDataEntry>::iterator pdtEntry
907 0 : = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
908 0 : particlePtr = &pdtEntry->second;
909 :
910 : // Print particle properties.
911 0 : os << "<particle id=\"" << particlePtr->id() << "\""
912 0 : << " name=\"" << particlePtr->name() << "\"";
913 0 : if (particlePtr->hasAnti())
914 0 : os << " antiName=\"" << particlePtr->name(-1) << "\"";
915 0 : os << " spinType=\"" << particlePtr->spinType() << "\""
916 0 : << " chargeType=\"" << particlePtr->chargeType() << "\""
917 0 : << " colType=\"" << particlePtr->colType() << "\"\n";
918 : // Pick format for mass and width based on mass value.
919 0 : double m0Now = particlePtr->m0();
920 0 : if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
921 0 : os << fixed << setprecision(5);
922 0 : else os << scientific << setprecision(3);
923 0 : os << " m0=\"" << m0Now << "\"";
924 0 : if (particlePtr->mWidth() > 0.)
925 0 : os << " mWidth=\"" << particlePtr->mWidth() << "\""
926 0 : << " mMin=\"" << particlePtr->mMin() << "\""
927 0 : << " mMax=\"" << particlePtr->mMax() << "\"";
928 0 : if (particlePtr->tau0() > 0.) os << scientific << setprecision(5)
929 0 : << " tau0=\"" << particlePtr->tau0() << "\"";
930 0 : os << ">\n";
931 :
932 : // Loop through the decay channel table for each particle.
933 0 : if (particlePtr->sizeChannels() > 0) {
934 0 : for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
935 0 : const DecayChannel& channel = particlePtr->channel(i);
936 0 : int mult = channel.multiplicity();
937 :
938 : // Print decay channel properties.
939 0 : os << " <channel onMode=\"" << channel.onMode() << "\""
940 0 : << fixed << setprecision(7)
941 0 : << " bRatio=\"" << channel.bRatio() << "\"";
942 0 : if (channel.meMode() > 0)
943 0 : os << " meMode=\"" << channel.meMode() << "\"";
944 0 : os << " products=\"";
945 0 : for (int j = 0; j < mult; ++j) {
946 0 : os << channel.product(j);
947 0 : if (j < mult - 1) os << " ";
948 : }
949 :
950 : // Finish off line and loop over allowed decay channels.
951 0 : os << "\"/>\n";
952 : }
953 0 : }
954 :
955 : // Finish off existing particle.
956 0 : os << "</particle>\n\n";
957 :
958 : }
959 :
960 0 : }
961 :
962 : //--------------------------------------------------------------------------
963 :
964 : // Read in database from specific free format file.
965 :
966 : bool ParticleData::readFF(string inFile, bool reset) {
967 :
968 : // Normally reset whole database before beginning.
969 0 : if (reset) {pdt.clear(); isInit = false;}
970 :
971 : // Open file for read and check that instream is OK.
972 0 : const char* cstring = inFile.c_str();
973 0 : ifstream is(cstring);
974 0 : if (!is.good()) {
975 0 : infoPtr->errorMsg("Error in ParticleData::readFF:"
976 0 : " did not find file", inFile);
977 0 : return false;
978 : }
979 :
980 : // Read in one line at a time.
981 0 : particlePtr = 0;
982 0 : string line;
983 : bool readParticle = false;
984 0 : while ( getline(is, line) ) {
985 :
986 : // Empty lines begins new particle.
987 0 : if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) {
988 : readParticle = true;
989 0 : continue;
990 : }
991 :
992 : // Prepare to use standard read from line.
993 0 : istringstream readLine(line);
994 :
995 : // Read in a line with particle information.
996 0 : if (readParticle) {
997 :
998 : // Properties to be read.
999 0 : int idTmp;
1000 0 : string nameTmp, antiNameTmp;
1001 0 : int spinTypeTmp, chargeTypeTmp, colTypeTmp;
1002 0 : double m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp;
1003 0 : string mayTmp;
1004 :
1005 : // Do the reading.
1006 0 : readLine >> idTmp >> nameTmp >> antiNameTmp >> spinTypeTmp
1007 0 : >> chargeTypeTmp >> colTypeTmp >> m0Tmp >> mWidthTmp
1008 0 : >> mMinTmp >> mMaxTmp >> tau0Tmp;
1009 :
1010 : // Error printout if something went wrong.
1011 0 : if (!readLine) {
1012 0 : infoPtr->errorMsg("Error in ParticleData::readFF:"
1013 0 : " incomplete particle", line);
1014 0 : return false;
1015 : }
1016 :
1017 : // Erase if particle already exists.
1018 0 : if (isParticle(idTmp)) pdt.erase(idTmp);
1019 :
1020 : // Store new particle. Save pointer, to be used for decay channels.
1021 0 : addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1022 0 : colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1023 0 : particlePtr = particleDataEntryPtr(idTmp);
1024 : readParticle = false;
1025 :
1026 : // Read in a line with decay channel information.
1027 0 : } else {
1028 :
1029 : // Properties to be read.
1030 0 : int onMode = 0;
1031 0 : double bRatio = 0.;
1032 0 : int meMode = 0;
1033 0 : int prod0 = 0; int prod1 = 0; int prod2 = 0; int prod3 = 0;
1034 0 : int prod4 = 0; int prod5 = 0; int prod6 = 0; int prod7 = 0;
1035 :
1036 : // Read in data from stream. Need at least one decay product.
1037 0 : readLine >> onMode >> bRatio >> meMode >> prod0;
1038 0 : if (!readLine) {
1039 0 : infoPtr->errorMsg("Error in ParticleData::readFF:"
1040 0 : " incomplete decay channel", line);
1041 0 : return false;
1042 : }
1043 0 : readLine >> prod1 >> prod2 >> prod3 >> prod4 >> prod5
1044 0 : >> prod6 >> prod7;
1045 :
1046 : // Store new channel.
1047 0 : if (particlePtr == 0) {
1048 0 : infoPtr->errorMsg("Error in ParticleData::readFF:"
1049 0 : " orphan decay channel", line);
1050 0 : return false;
1051 : }
1052 0 : particlePtr->addChannel(onMode, bRatio, meMode, prod0, prod1,
1053 0 : prod2, prod3, prod4, prod5, prod6, prod7);
1054 :
1055 0 : }
1056 :
1057 : // End of while loop over lines in the file.
1058 0 : }
1059 :
1060 :
1061 : // Done.
1062 0 : isInit = true;
1063 0 : return true;
1064 :
1065 0 : }
1066 :
1067 : //--------------------------------------------------------------------------
1068 :
1069 : // Print out complete database in numerical order as a free format file.
1070 :
1071 : void ParticleData::listFF(string outFile) {
1072 :
1073 : // Convert file name to ofstream.
1074 0 : const char* cstring = outFile.c_str();
1075 0 : ofstream os(cstring);
1076 :
1077 : // Iterate through the particle data table.
1078 0 : for (map<int, ParticleDataEntry>::iterator pdtEntry
1079 0 : = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1080 0 : particlePtr = &pdtEntry->second;
1081 :
1082 : // Pick format for mass and width based on mass value.
1083 0 : double m0Now = particlePtr->m0();
1084 0 : if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1085 0 : os << fixed << setprecision(5);
1086 0 : else os << scientific << setprecision(3);
1087 :
1088 : // Print particle properties.
1089 0 : os << "\n" << setw(8) << particlePtr->id() << " "
1090 0 : << left << setw(16) << particlePtr->name() << " "
1091 0 : << setw(16) << particlePtr->name(-1) << " "
1092 0 : << right << setw(2) << particlePtr->spinType() << " "
1093 0 : << setw(2) << particlePtr->chargeType() << " "
1094 0 : << setw(2) << particlePtr->colType() << " "
1095 0 : << setw(10) << particlePtr->m0() << " "
1096 0 : << setw(10) << particlePtr->mWidth() << " "
1097 0 : << setw(10) << particlePtr->mMin() << " "
1098 0 : << setw(10) << particlePtr->mMax() << " "
1099 0 : << scientific << setprecision(5)
1100 0 : << setw(12) << particlePtr->tau0() << "\n";
1101 :
1102 : // Loop through the decay channel table for each particle.
1103 0 : if (particlePtr->sizeChannels() > 0) {
1104 0 : for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1105 0 : const DecayChannel& channel = particlePtr->channel(i);
1106 0 : os << " " << setw(6) << channel.onMode()
1107 0 : << " " << fixed << setprecision(7) << setw(10)
1108 0 : << channel.bRatio() << " "
1109 0 : << setw(3) << channel.meMode() << " ";
1110 0 : for (int j = 0; j < channel.multiplicity(); ++j)
1111 0 : os << setw(8) << channel.product(j) << " ";
1112 0 : os << "\n";
1113 : }
1114 0 : }
1115 :
1116 : }
1117 :
1118 0 : }
1119 :
1120 : //--------------------------------------------------------------------------
1121 :
1122 : // Read in updates from a character string, like a line of a file.
1123 : // Is used by readString (and readFile) in Pythia.
1124 :
1125 : bool ParticleData::readString(string lineIn, bool warn, ostream& os) {
1126 :
1127 : // If empty line then done.
1128 0 : if (lineIn.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return true;
1129 :
1130 : // Take copy that will be modified.
1131 0 : string line = lineIn;
1132 :
1133 : // If first character is not a digit then taken to be a comment.
1134 0 : int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
1135 0 : if (!isdigit(line[firstChar])) return true;
1136 :
1137 : // Replace colons and equal signs by blanks to make parsing simpler.
1138 0 : for ( int j = 0; j < int(line.size()); ++ j)
1139 0 : if (line[j] == ':' || line[j] == '=') line[j] = ' ';
1140 :
1141 : // Get particle id and property name.
1142 0 : int idTmp;
1143 0 : string property;
1144 0 : istringstream getWord(line);
1145 0 : getWord >> idTmp >> property;
1146 0 : property = toLower(property);
1147 :
1148 : // Check that valid particle.
1149 0 : if ( (!isParticle(idTmp) && property != "all" && property != "new")
1150 0 : || idTmp <= 0) {
1151 0 : if (warn) os << "\n PYTHIA Error: input particle not found in Particle"
1152 0 : << " Data Table:\n " << lineIn << "\n";
1153 0 : readingFailedSave = true;
1154 0 : return false;
1155 : }
1156 :
1157 : // Identify particle property and read + set its value, case by case.
1158 0 : if (property == "name") {
1159 0 : string nameTmp;
1160 0 : getWord >> nameTmp;
1161 0 : pdt[idTmp].setName(nameTmp);
1162 : return true;
1163 0 : }
1164 0 : if (property == "antiname") {
1165 0 : string antiNameTmp;
1166 0 : getWord >> antiNameTmp;
1167 0 : pdt[idTmp].setAntiName(antiNameTmp);
1168 : return true;
1169 0 : }
1170 0 : if (property == "names") {
1171 0 : string nameTmp, antiNameTmp;
1172 0 : getWord >> nameTmp >> antiNameTmp;
1173 0 : pdt[idTmp].setNames(nameTmp, antiNameTmp);
1174 : return true;
1175 0 : }
1176 0 : if (property == "spintype") {
1177 0 : int spinTypeTmp;
1178 0 : getWord >> spinTypeTmp;
1179 0 : pdt[idTmp].setSpinType(spinTypeTmp);
1180 : return true;
1181 0 : }
1182 0 : if (property == "chargetype") {
1183 0 : int chargeTypeTmp;
1184 0 : getWord >> chargeTypeTmp;
1185 0 : pdt[idTmp].setChargeType(chargeTypeTmp);
1186 : return true;
1187 0 : }
1188 0 : if (property == "coltype") {
1189 0 : int colTypeTmp;
1190 0 : getWord >> colTypeTmp;
1191 0 : pdt[idTmp].setColType(colTypeTmp);
1192 : return true;
1193 0 : }
1194 0 : if (property == "m0") {
1195 0 : double m0Tmp;
1196 0 : getWord >> m0Tmp;
1197 0 : pdt[idTmp].setM0(m0Tmp);
1198 : return true;
1199 0 : }
1200 0 : if (property == "mwidth") {
1201 0 : double mWidthTmp;
1202 0 : getWord >> mWidthTmp;
1203 0 : pdt[idTmp].setMWidth(mWidthTmp);
1204 : return true;
1205 0 : }
1206 0 : if (property == "mmin") {
1207 0 : double mMinTmp;
1208 0 : getWord >> mMinTmp;
1209 0 : pdt[idTmp].setMMin(mMinTmp);
1210 : return true;
1211 0 : }
1212 0 : if (property == "mmax") {
1213 0 : double mMaxTmp;
1214 0 : getWord >> mMaxTmp;
1215 0 : pdt[idTmp].setMMax(mMaxTmp);
1216 : return true;
1217 0 : }
1218 0 : if (property == "tau0") {
1219 0 : double tau0Tmp;
1220 0 : getWord >> tau0Tmp;
1221 0 : pdt[idTmp].setTau0(tau0Tmp);
1222 : return true;
1223 0 : }
1224 0 : if (property == "isresonance") {
1225 0 : string isresTmp;
1226 0 : getWord >> isresTmp;
1227 0 : bool isResonanceTmp = boolString(isresTmp);
1228 0 : pdt[idTmp].setIsResonance(isResonanceTmp);
1229 : return true;
1230 0 : }
1231 0 : if (property == "maydecay") {
1232 0 : string mayTmp;
1233 0 : getWord >> mayTmp;
1234 0 : bool mayDecayTmp = boolString(mayTmp);
1235 0 : pdt[idTmp].setMayDecay(mayDecayTmp);
1236 : return true;
1237 0 : }
1238 0 : if (property == "doexternaldecay") {
1239 0 : string extdecTmp;
1240 0 : getWord >> extdecTmp;
1241 0 : bool doExternalDecayTmp = boolString(extdecTmp);
1242 0 : pdt[idTmp].setDoExternalDecay(doExternalDecayTmp);
1243 : return true;
1244 0 : }
1245 0 : if (property == "isvisible") {
1246 0 : string isvisTmp;
1247 0 : getWord >> isvisTmp;
1248 0 : bool isVisibleTmp = boolString(isvisTmp);
1249 0 : pdt[idTmp].setIsVisible(isVisibleTmp);
1250 : return true;
1251 0 : }
1252 0 : if (property == "doforcewidth") {
1253 0 : string doforceTmp;
1254 0 : getWord >> doforceTmp;
1255 0 : bool doForceWidthTmp = boolString(doforceTmp);
1256 0 : pdt[idTmp].setDoForceWidth(doForceWidthTmp);
1257 : return true;
1258 0 : }
1259 :
1260 : // Addition or complete replacement of a particle.
1261 0 : if (property == "all" || property == "new") {
1262 :
1263 : // Default values for properties to be read.
1264 0 : string nameTmp = "void";
1265 0 : string antiNameTmp = "void";
1266 0 : int spinTypeTmp = 0;
1267 0 : int chargeTypeTmp = 0;
1268 0 : int colTypeTmp = 0;
1269 0 : double m0Tmp = 0.;
1270 0 : double mWidthTmp = 0.;
1271 0 : double mMinTmp = 0.;
1272 0 : double mMaxTmp = 0.;
1273 0 : double tau0Tmp = 0.;
1274 :
1275 : // Read in data from stream.
1276 0 : getWord >> nameTmp >> antiNameTmp >> spinTypeTmp >> chargeTypeTmp
1277 0 : >> colTypeTmp >> m0Tmp >> mWidthTmp >> mMinTmp >> mMaxTmp
1278 0 : >> tau0Tmp;
1279 :
1280 : // To keep existing decay channels, only overwrite particle data.
1281 0 : if (property == "all" && isParticle(idTmp)) {
1282 0 : setAll( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1283 0 : colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1284 :
1285 : // Else start over completely from scratch.
1286 0 : } else {
1287 0 : if (isParticle(idTmp)) pdt.erase(idTmp);
1288 0 : addParticle( idTmp, nameTmp, antiNameTmp, spinTypeTmp, chargeTypeTmp,
1289 0 : colTypeTmp, m0Tmp, mWidthTmp, mMinTmp, mMaxTmp, tau0Tmp);
1290 : }
1291 : return true;
1292 0 : }
1293 :
1294 : // Set onMode of all decay channels in one go.
1295 0 : if (property == "onmode") {
1296 0 : int onMode = 0;
1297 0 : string onModeIn;
1298 0 : getWord >> onModeIn;
1299 : // For onMode allow the optional possibility of Bool input.
1300 0 : if (isdigit(onModeIn[0])) {
1301 0 : istringstream getOnMode(onModeIn);
1302 0 : getOnMode >> onMode;
1303 0 : } else onMode = (boolString(onModeIn)) ? 1 : 0;
1304 0 : for (int i = 0; i < pdt[idTmp].sizeChannels(); ++i)
1305 0 : pdt[idTmp].channel(i).onMode(onMode);
1306 : return true;
1307 0 : }
1308 :
1309 : // Selective search for matching decay channels.
1310 : int matchKind = 0;
1311 0 : if (property == "offifany" || property == "onifany" ||
1312 0 : property == "onposifany" || property == "onnegifany") matchKind = 1;
1313 0 : if (property == "offifall" || property == "onifall" ||
1314 0 : property == "onposifall" || property == "onnegifall") matchKind = 2;
1315 0 : if (property == "offifmatch" || property == "onifmatch" ||
1316 0 : property == "onposifmatch" || property == "onnegifmatch") matchKind = 3;
1317 0 : if (matchKind > 0) {
1318 : int onMode = 0;
1319 0 : if (property == "onifany" || property == "onifall"
1320 0 : || property == "onifmatch") onMode = 1;
1321 0 : if (property == "onposifany" || property == "onposifall"
1322 0 : || property == "onposifmatch") onMode = 2;
1323 0 : if (property == "onnegifany" || property == "onnegifall"
1324 0 : || property == "onnegifmatch") onMode = 3;
1325 :
1326 : // Read in particles to match.
1327 0 : vector<int> idToMatch;
1328 0 : int idRead;
1329 0 : for ( ; ; ) {
1330 0 : getWord >> idRead;
1331 0 : if (!getWord) break;
1332 0 : idToMatch.push_back(abs(idRead));
1333 : }
1334 0 : int nToMatch = idToMatch.size();
1335 :
1336 : // Loop over all decay channels.
1337 0 : for (int i = 0; i < pdt[idTmp].sizeChannels(); ++i) {
1338 0 : int multi = pdt[idTmp].channel(i).multiplicity();
1339 :
1340 : // Look for any match at all.
1341 0 : if (matchKind == 1) {
1342 : bool foundMatch = false;
1343 0 : for (int j = 0; j < multi; ++j) {
1344 0 : int idNow = abs(pdt[idTmp].channel(i).product(j));
1345 0 : for (int k = 0; k < nToMatch; ++k)
1346 0 : if (idNow == idToMatch[k]) {foundMatch = true; break;}
1347 0 : if (foundMatch) break;
1348 0 : }
1349 0 : if (foundMatch) pdt[idTmp].channel(i).onMode(onMode);
1350 :
1351 : // Look for match of all products provided.
1352 0 : } else {
1353 : int nUnmatched = nToMatch;
1354 0 : if (multi < nToMatch);
1355 0 : else if (multi > nToMatch && matchKind == 3);
1356 : else {
1357 0 : vector<int> idUnmatched;
1358 0 : for (int k = 0; k < nToMatch; ++k)
1359 0 : idUnmatched.push_back(idToMatch[k]);
1360 0 : for (int j = 0; j < multi; ++j) {
1361 0 : int idNow = abs(pdt[idTmp].channel(i).product(j));
1362 0 : for (int k = 0; k < nUnmatched; ++k)
1363 0 : if (idNow == idUnmatched[k]) {
1364 0 : idUnmatched[k] = idUnmatched[--nUnmatched];
1365 0 : break;
1366 : }
1367 0 : if (nUnmatched == 0) break;
1368 0 : }
1369 0 : }
1370 0 : if (nUnmatched == 0) pdt[idTmp].channel(i).onMode(onMode);
1371 : }
1372 : }
1373 : return true;
1374 0 : }
1375 :
1376 : // Rescale all branching ratios by common factor.
1377 0 : if (property == "rescalebr") {
1378 0 : double factor;
1379 0 : getWord >> factor;
1380 0 : pdt[idTmp].rescaleBR(factor);
1381 : return true;
1382 0 : }
1383 :
1384 : // Reset decay table in preparation for new input.
1385 0 : if (property == "onechannel") pdt[idTmp].clearChannels();
1386 :
1387 : // Add or change a decay channel: get channel number and new property.
1388 0 : if (property == "addchannel" || property == "onechannel"
1389 0 : || isdigit(property[0])) {
1390 0 : int channel;
1391 0 : if (property == "addchannel" || property == "onechannel") {
1392 0 : pdt[idTmp].addChannel();
1393 0 : channel = pdt[idTmp].sizeChannels() - 1;
1394 0 : property = "all";
1395 : } else{
1396 0 : istringstream getChannel(property);
1397 0 : getChannel >> channel;
1398 0 : getWord >> property;
1399 0 : property = toLower(property);
1400 0 : }
1401 :
1402 : // Check that channel exists.
1403 0 : if (channel < 0 || channel >= pdt[idTmp].sizeChannels()) return false;
1404 :
1405 : // Find decay channel property and value, case by case.
1406 : // At same time also do case where all should be replaced.
1407 0 : if (property == "onmode" || property == "all") {
1408 0 : int onMode = 0;
1409 0 : string onModeIn;
1410 0 : getWord >> onModeIn;
1411 : // For onMode allow the optional possibility of Bool input.
1412 0 : if (isdigit(onModeIn[0])) {
1413 0 : istringstream getOnMode(onModeIn);
1414 0 : getOnMode >> onMode;
1415 0 : } else onMode = (boolString(onModeIn)) ? 1 : 0;
1416 0 : pdt[idTmp].channel(channel).onMode(onMode);
1417 0 : if (property == "onmode") return true;
1418 0 : }
1419 0 : if (property == "bratio" || property == "all") {
1420 0 : double bRatio;
1421 0 : getWord >> bRatio;
1422 0 : pdt[idTmp].channel(channel).bRatio(bRatio);
1423 0 : if (property == "bratio") return true;
1424 0 : }
1425 0 : if (property == "memode" || property == "all") {
1426 0 : int meMode;
1427 0 : getWord >> meMode;
1428 0 : pdt[idTmp].channel(channel).meMode(meMode);
1429 0 : if (property == "memode") return true;
1430 0 : }
1431 :
1432 : // Scan for products until end of line.
1433 0 : if (property == "products" || property == "all") {
1434 : int nProd = 0;
1435 0 : for (int iProd = 0; iProd < 8; ++iProd) {
1436 0 : int idProd;
1437 0 : getWord >> idProd;
1438 0 : if (!getWord) break;
1439 0 : pdt[idTmp].channel(channel).product(iProd, idProd);
1440 0 : ++nProd;
1441 0 : }
1442 0 : for (int iProd = nProd; iProd < 8; ++iProd)
1443 0 : pdt[idTmp].channel(channel).product(iProd, 0);
1444 0 : pdt[idTmp].channel(channel).multiplicity(nProd);
1445 : return true;
1446 : }
1447 :
1448 : // Rescale an existing branching ratio.
1449 0 : if (property == "rescalebr") {
1450 0 : double factor;
1451 0 : getWord >> factor;
1452 0 : pdt[idTmp].channel(channel).rescaleBR(factor);
1453 : return true;
1454 0 : }
1455 0 : }
1456 :
1457 : // Return false if failed to recognize property.
1458 0 : if (warn) os << "\n PYTHIA Error: input property not found in Particle"
1459 0 : << " Data Table:\n " << lineIn << "\n";
1460 0 : readingFailedSave = true;
1461 0 : return false;
1462 :
1463 0 : }
1464 :
1465 : //--------------------------------------------------------------------------
1466 :
1467 : // Print out complete or changed table of database in numerical order.
1468 :
1469 : void ParticleData::list(bool changedOnly, bool changedRes, ostream& os) {
1470 :
1471 : // Table header; output for bool as off/on.
1472 0 : if (!changedOnly) {
1473 0 : os << "\n -------- PYTHIA Particle Data Table (complete) --------"
1474 0 : << "------------------------------------------------------------"
1475 0 : << "--------------\n \n";
1476 :
1477 0 : } else {
1478 0 : os << "\n -------- PYTHIA Particle Data Table (changed only) ----"
1479 0 : << "------------------------------------------------------------"
1480 0 : << "--------------\n \n";
1481 : }
1482 0 : os << " id name antiName spn chg col m0"
1483 0 : << " mWidth mMin mMax tau0 res dec ext "
1484 0 : << "vis wid\n no onMode bRatio meMode products \n";
1485 :
1486 : // Iterate through the particle data table. Option to skip unchanged.
1487 : int nList = 0;
1488 0 : for (map<int, ParticleDataEntry>::iterator pdtEntry
1489 0 : = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1490 0 : particlePtr = &pdtEntry->second;
1491 0 : if ( !changedOnly || particlePtr->hasChanged() ||
1492 0 : ( changedRes && particlePtr->getResonancePtr() != 0 ) ) {
1493 :
1494 : // Pick format for mass and width based on mass value.
1495 0 : double m0Now = particlePtr->m0();
1496 0 : if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1497 0 : os << fixed << setprecision(5);
1498 0 : else os << scientific << setprecision(3);
1499 :
1500 : // Print particle properties.
1501 0 : ++nList;
1502 0 : os << "\n" << setw(8) << particlePtr->id() << " " << left;
1503 0 : if (particlePtr->name(-1) == "void")
1504 0 : os << setw(33) << particlePtr->name() << " ";
1505 0 : else os << setw(16) << particlePtr->name() << " "
1506 0 : << setw(16) << particlePtr->name(-1) << " ";
1507 0 : os << right << setw(2) << particlePtr->spinType() << " "
1508 0 : << setw(2) << particlePtr->chargeType() << " "
1509 0 : << setw(2) << particlePtr->colType() << " "
1510 0 : << setw(10) << particlePtr->m0() << " "
1511 0 : << setw(10) << particlePtr->mWidth() << " "
1512 0 : << setw(10) << particlePtr->mMin() << " "
1513 0 : << setw(10) << particlePtr->mMax() << " "
1514 0 : << scientific << setprecision(5)
1515 0 : << setw(12) << particlePtr->tau0() << " " << setw(2)
1516 0 : << particlePtr->isResonance() << " " << setw(2)
1517 0 : << (particlePtr->mayDecay() && particlePtr->canDecay())
1518 0 : << " " << setw(2) << particlePtr->doExternalDecay() << " "
1519 0 : << setw(2) << particlePtr->isVisible()<< " "
1520 0 : << setw(2) << particlePtr->doForceWidth() << "\n";
1521 :
1522 : // Loop through the decay channel table for each particle.
1523 0 : if (particlePtr->sizeChannels() > 0) {
1524 0 : for (int i = 0; i < int(particlePtr->sizeChannels()); ++i) {
1525 0 : const DecayChannel& channel = particlePtr->channel(i);
1526 0 : os << " " << setprecision(7)
1527 0 : << setw(5) << i
1528 0 : << setw(6) << channel.onMode()
1529 0 : << fixed<< setw(12) << channel.bRatio()
1530 0 : << setw(5) << channel.meMode() << " ";
1531 0 : for (int j = 0; j < channel.multiplicity(); ++j)
1532 0 : os << setw(8) << channel.product(j) << " ";
1533 0 : os << "\n";
1534 : }
1535 0 : }
1536 0 : }
1537 :
1538 : }
1539 :
1540 : // End of loop over database contents.
1541 0 : if (changedOnly && nList == 0) os << "\n no particle data has been "
1542 0 : << "changed from its default value \n";
1543 0 : os << "\n -------- End PYTHIA Particle Data Table -----------------"
1544 0 : << "--------------------------------------------------------------"
1545 0 : << "----------\n" << endl;
1546 :
1547 0 : }
1548 :
1549 : //--------------------------------------------------------------------------
1550 :
1551 : // Print out partial table of database in input order.
1552 :
1553 : void ParticleData::list(vector<int> idList, ostream& os) {
1554 :
1555 : // Table header; output for bool as off/on.
1556 0 : os << "\n -------- PYTHIA Particle Data Table (partial) ---------"
1557 0 : << "------------------------------------------------------------"
1558 0 : << "--------------\n \n";
1559 0 : os << " id name antiName spn chg col m0"
1560 0 : << " mWidth mMin mMax tau0 res dec ext "
1561 0 : << "vis wid\n no onMode bRatio meMode products \n";
1562 :
1563 : // Iterate through the given list of input particles.
1564 0 : for (int i = 0; i < int(idList.size()); ++i) {
1565 0 : particlePtr = particleDataEntryPtr(idList[i]);
1566 :
1567 : // Pick format for mass and width based on mass value.
1568 0 : double m0Now = particlePtr->m0();
1569 0 : if (m0Now == 0 || (m0Now > 0.1 && m0Now < 1000.))
1570 0 : os << fixed << setprecision(5);
1571 0 : else os << scientific << setprecision(3);
1572 :
1573 : // Print particle properties.
1574 0 : os << "\n" << setw(8) << particlePtr->id() << " " << left;
1575 0 : if (particlePtr->name(-1) == "void")
1576 0 : os << setw(33) << particlePtr->name() << " ";
1577 0 : else os << setw(16) << particlePtr->name() << " "
1578 0 : << setw(16) << particlePtr->name(-1) << " ";
1579 0 : os << right << setw(2) << particlePtr->spinType() << " "
1580 0 : << setw(2) << particlePtr->chargeType() << " "
1581 0 : << setw(2) << particlePtr->colType() << " "
1582 0 : << setw(10) << particlePtr->m0() << " "
1583 0 : << setw(10) << particlePtr->mWidth() << " "
1584 0 : << setw(10) << particlePtr->mMin() << " "
1585 0 : << setw(10) << particlePtr->mMax() << " "
1586 0 : << scientific << setprecision(5)
1587 0 : << setw(12) << particlePtr->tau0() << " " << setw(2)
1588 0 : << particlePtr->isResonance() << " " << setw(2)
1589 0 : << (particlePtr->mayDecay() && particlePtr->canDecay())
1590 0 : << " " << setw(2) << particlePtr->doExternalDecay() << " "
1591 0 : << setw(2) << particlePtr->isVisible() << " "
1592 0 : << setw(2) << particlePtr->doForceWidth() << "\n";
1593 :
1594 : // Loop through the decay channel table for each particle.
1595 0 : if (particlePtr->sizeChannels() > 0) {
1596 0 : for (int j = 0; j < int(particlePtr->sizeChannels()); ++j) {
1597 0 : const DecayChannel& channel = particlePtr->channel(j);
1598 0 : os << " " << setprecision(7)
1599 0 : << setw(5) << j
1600 0 : << setw(6) << channel.onMode()
1601 0 : << fixed<< setw(12) << channel.bRatio()
1602 0 : << setw(5) << channel.meMode() << " ";
1603 0 : for (int k = 0; k < channel.multiplicity(); ++k)
1604 0 : os << setw(8) << channel.product(k) << " ";
1605 0 : os << "\n";
1606 : }
1607 0 : }
1608 :
1609 : }
1610 :
1611 : // End of loop over database contents.
1612 0 : os << "\n -------- End PYTHIA Particle Data Table -----------------"
1613 0 : << "--------------------------------------------------------------"
1614 0 : << "----------\n" << endl;
1615 :
1616 0 : }
1617 :
1618 : //--------------------------------------------------------------------------
1619 :
1620 : // Check that table makes sense: e.g. consistent names and mass ranges,
1621 : // that branching ratios sum to unity, that charge is conserved and
1622 : // that phase space is open in each channel.
1623 : // verbosity = 0: mimimal amount of checks, e.g. that no channels open.
1624 : // = 1: further warning if individual channels closed
1625 : // (except for resonances).
1626 : // = 2: also print branching-ratio-averaged threshold mass.
1627 : // = 11, 12: as 1, 2, but include resonances in detailed checks.
1628 :
1629 : void ParticleData::checkTable(int verbosity, ostream& os) {
1630 :
1631 : // Header.
1632 0 : os << "\n -------- PYTHIA Check of Particle Data Table ------------"
1633 0 : <<"------\n\n";
1634 : int nErr = 0;
1635 :
1636 : // Loop through all particles.
1637 0 : for (map<int, ParticleDataEntry>::iterator pdtEntry
1638 0 : = pdt.begin(); pdtEntry != pdt.end(); ++pdtEntry) {
1639 0 : particlePtr = &pdtEntry->second;
1640 :
1641 : // Extract some particle properties. Set some flags.
1642 0 : int idNow = particlePtr->id();
1643 0 : bool hasAntiNow = particlePtr->hasAnti();
1644 0 : int spinTypeNow = particlePtr->spinType();
1645 0 : int chargeTypeNow = particlePtr->chargeType();
1646 0 : int baryonTypeNow = particlePtr->baryonNumberType();
1647 0 : double m0Now = particlePtr->m0();
1648 0 : double mMinNow = particlePtr->mMin();
1649 0 : double mMaxNow = particlePtr->mMax();
1650 0 : double mWidthNow = particlePtr->mWidth();
1651 0 : double tau0Now = particlePtr->tau0();
1652 0 : bool isResonanceNow = particlePtr->isResonance();
1653 : bool hasPrinted = false;
1654 0 : bool studyCloser = verbosity > 10 || !isResonanceNow;
1655 :
1656 : // Check that particle name consistent with charge information.
1657 0 : string particleName = particlePtr->name(1);
1658 0 : if (particleName.size() > 16) {
1659 0 : os << " Warning: particle " << idNow << " has name " << particleName
1660 0 : << " of length " << particleName.size() << "\n";
1661 : hasPrinted = true;
1662 0 : ++nErr;
1663 0 : }
1664 : int nPos = 0;
1665 : int nNeg = 0;
1666 0 : for (int i = 0; i < int(particleName.size()); ++i) {
1667 0 : if (particleName[i] == '+') ++nPos;
1668 0 : if (particleName[i] == '-') ++nNeg;
1669 : }
1670 0 : if ( (nPos > 0 && nNeg > 0) || ( nPos + nNeg > 0
1671 0 : && 3 * (nPos - nNeg) != chargeTypeNow )) {
1672 0 : os << " Warning: particle " << idNow << " has name " << particleName
1673 0 : << " inconsistent with charge type " << chargeTypeNow << "\n";
1674 : hasPrinted = true;
1675 0 : ++nErr;
1676 0 : }
1677 :
1678 : // Check that antiparticle name consistent with charge information.
1679 0 : if (hasAntiNow) {
1680 0 : particleName = particlePtr->name(-1);
1681 0 : if (particleName.size() > 16) {
1682 0 : os << " Warning: particle " << idNow << " has name " << particleName
1683 0 : << " of length " << particleName.size() << "\n";
1684 : hasPrinted = true;
1685 0 : ++nErr;
1686 0 : }
1687 : nPos = 0;
1688 : nNeg = 0;
1689 0 : for (int i = 0; i < int(particleName.size()); ++i) {
1690 0 : if (particleName[i] == '+') ++nPos;
1691 0 : if (particleName[i] == '-') ++nNeg;
1692 : }
1693 0 : if ( (nPos > 0 && nNeg > 0) || ( nPos + nNeg > 0
1694 0 : && 3 * (nPos - nNeg) != -chargeTypeNow )) {
1695 0 : os << " Warning: particle " << -idNow << " has name "
1696 0 : << particleName << " inconsistent with charge type "
1697 0 : << -chargeTypeNow << "\n";
1698 : hasPrinted = true;
1699 0 : ++nErr;
1700 0 : }
1701 : }
1702 :
1703 : // Check that mass, mass range and width are consistent.
1704 0 : if (particlePtr->useBreitWigner()) {
1705 0 : if (mMinNow > m0Now) {
1706 0 : os << " Error: particle " << idNow << " has mMin "
1707 0 : << fixed << setprecision(5) << mMinNow
1708 0 : << " larger than m0 " << m0Now << "\n";
1709 : hasPrinted = true;
1710 0 : ++nErr;
1711 0 : }
1712 0 : if (mMaxNow > mMinNow && mMaxNow < m0Now) {
1713 0 : os << " Error: particle " << idNow << " has mMax "
1714 0 : << fixed << setprecision(5) << mMaxNow
1715 0 : << " smaller than m0 " << m0Now << "\n";
1716 : hasPrinted = true;
1717 0 : ++nErr;
1718 0 : }
1719 0 : if (mMaxNow > mMinNow && mMaxNow - mMinNow < mWidthNow) {
1720 0 : os << " Warning: particle " << idNow << " has mMax - mMin "
1721 0 : << fixed << setprecision(5) << mMaxNow - mMinNow
1722 0 : << " smaller than mWidth " << mWidthNow << "\n";
1723 : hasPrinted = true;
1724 0 : ++nErr;
1725 0 : }
1726 : }
1727 :
1728 : // Check that particle does not both have width and lifetime.
1729 0 : if (mWidthNow > 0. && tau0Now > 0.) {
1730 0 : os << " Warning: particle " << idNow << " has both nonvanishing width "
1731 0 : << scientific << setprecision(5) << mWidthNow << " and lifetime "
1732 0 : << tau0Now << "\n";
1733 : hasPrinted = true;
1734 0 : ++nErr;
1735 0 : }
1736 :
1737 : // Begin study decay channels.
1738 0 : if (particlePtr->sizeChannels() > 0) {
1739 :
1740 : // Loop through all decay channels.
1741 : double bRatioSum = 0.;
1742 : double bRatioPos = 0.;
1743 : double bRatioNeg = 0.;
1744 : bool hasPosBR = false;
1745 : bool hasNegBR = false;
1746 : double threshMass = 0.;
1747 : bool openChannel = false;
1748 0 : for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1749 :
1750 : // Extract channel properties.
1751 0 : int onMode = particlePtr->channel(i).onMode();
1752 0 : double bRatio = particlePtr->channel(i).bRatio();
1753 0 : int meMode = particlePtr->channel(i).meMode();
1754 0 : int mult = particlePtr->channel(i).multiplicity();
1755 0 : int prod[8];
1756 0 : for (int j = 0; j < 8; ++j)
1757 0 : prod[j] = particlePtr->channel(i).product(j);
1758 :
1759 : // Sum branching ratios. Include off-channels.
1760 0 : if (onMode == 0 || onMode == 1) bRatioSum += bRatio;
1761 0 : else if (onMode == 2) {bRatioPos += bRatio; hasPosBR = true;}
1762 0 : else if (onMode == 3) {bRatioNeg += bRatio; hasNegBR = true;}
1763 :
1764 : // Error printout when unknown decay product code.
1765 0 : for (int j = 0; j < 8; ++j) {
1766 0 : if ( prod[j] != 0 && !isParticle(prod[j]) ) {
1767 0 : os << " Error: unknown decay product for " << idNow
1768 0 : << " -> " << prod[j] << "\n";
1769 : hasPrinted = true;
1770 0 : ++nErr;
1771 0 : continue;
1772 : }
1773 : }
1774 :
1775 : // Error printout when no decay products or 0 interspersed.
1776 : int nLast = 0;
1777 0 : for (int j = 0; j < 8; ++j)
1778 0 : if (prod[j] != 0) nLast = j + 1;
1779 0 : if (mult == 0 || mult != nLast) {
1780 0 : os << " Error: corrupted decay product list for "
1781 0 : << particlePtr->id() << " -> ";
1782 0 : for (int j = 0; j < 8; ++j) os << prod[j] << " ";
1783 0 : os << "\n";
1784 : hasPrinted = true;
1785 0 : ++nErr;
1786 0 : continue;
1787 : }
1788 :
1789 : // Check charge conservation and open phase space in decay channel.
1790 0 : int chargeTypeSum = -chargeTypeNow;
1791 0 : int baryonTypeSum = -baryonTypeNow;
1792 : double avgFinalMass = 0.;
1793 : double minFinalMass = 0.;
1794 : bool canHandle = true;
1795 0 : for (int j = 0; j < mult; ++j) {
1796 0 : chargeTypeSum += chargeType( prod[j] );
1797 0 : baryonTypeSum += baryonNumberType( prod[j] );
1798 0 : avgFinalMass += m0( prod[j] );
1799 0 : minFinalMass += m0Min( prod[j] );
1800 0 : if (prod[j] == 81 || prod[j] == 82 || prod[j] == 83)
1801 0 : canHandle = false;
1802 : }
1803 0 : threshMass += bRatio * avgFinalMass;
1804 :
1805 : // Error printout when charge or baryon number not conserved.
1806 0 : if (chargeTypeSum != 0 && canHandle) {
1807 0 : os << " Error: 3*charge changed by " << chargeTypeSum
1808 0 : << " in " << idNow << " -> ";
1809 0 : for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1810 0 : os << "\n";
1811 : hasPrinted = true;
1812 0 : ++nErr;
1813 0 : continue;
1814 : }
1815 0 : if ( baryonTypeSum != 0 && canHandle && particlePtr->isHadron() ) {
1816 0 : os << " Error: 3*baryon number changed by " << baryonTypeSum
1817 0 : << " in " << idNow << " -> ";
1818 0 : for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1819 0 : os << "\n";
1820 : hasPrinted = true;
1821 0 : ++nErr;
1822 0 : continue;
1823 : }
1824 :
1825 : // Begin check that some matrix elements are used correctly.
1826 : bool correctME = true;
1827 :
1828 : // Check matrix element mode 0: recommended not into partons.
1829 0 : if (meMode == 0 && !isResonanceNow) {
1830 : bool hasPartons = false;
1831 0 : for (int j = 0; j < mult; ++j) {
1832 0 : int idAbs = abs(prod[j]);
1833 0 : if ( idAbs < 10 || idAbs == 21 || idAbs == 81 || idAbs == 82
1834 0 : || idAbs == 83 || (idAbs > 1000 && idAbs < 10000
1835 0 : && (idAbs/10)%10 == 0) ) hasPartons = true;
1836 : }
1837 0 : if (hasPartons) correctME = false;
1838 0 : }
1839 :
1840 : // Check matrix element mode 1: rho/omega -> pi+ pi- pi0.
1841 0 : bool useME1 = ( mult == 3 && spinTypeNow == 3 && idNow > 100
1842 0 : && idNow < 1000
1843 0 : && particlePtr->channel(i).contains(211, -211, 111) );
1844 0 : if ( meMode == 1 && !useME1 ) correctME = false;
1845 0 : if ( meMode != 1 && useME1 ) correctME = false;
1846 :
1847 : // Check matrix element mode 2: polarization in V -> PS + PS.
1848 0 : bool useME2 = ( mult == 2 && spinTypeNow == 3 && idNow > 100
1849 0 : && idNow < 1000 && spinType(prod[0]) == 1
1850 0 : && spinType(prod[1]) == 1 );
1851 0 : if ( meMode == 2 && !useME2 ) correctME = false;
1852 0 : if ( meMode != 2 && useME2 ) correctME = false;
1853 :
1854 : // Check matrix element mode 11, 12 and 13: Dalitz decay with
1855 : // one or more particles in addition to the lepton pair,
1856 : // or double Dalitz decay.
1857 0 : bool useME11 = ( mult == 3 && !isResonanceNow
1858 0 : && (prod[1] == 11 || prod[1] == 13 || prod[1] == 15)
1859 0 : && prod[2] == -prod[1] );
1860 0 : bool useME12 = ( mult > 3 && !isResonanceNow
1861 0 : && (prod[mult - 2] == 11 || prod[mult - 2] == 13
1862 0 : || prod[mult - 2] == 15) && prod[mult - 1] == -prod[mult - 2] );
1863 0 : bool useME13 = ( mult == 4 && !isResonanceNow
1864 0 : && (prod[0] == 11 || prod[0] == 13) && prod[1] == -prod[0]
1865 0 : && (prod[2] == 11 || prod[2] == 13) && prod[3] == -prod[2] );
1866 0 : if (useME13) useME12 = false;
1867 0 : if ( meMode == 11 && !useME11 ) correctME = false;
1868 0 : if ( meMode != 11 && useME11 ) correctME = false;
1869 0 : if ( meMode == 12 && !useME12 ) correctME = false;
1870 0 : if ( meMode != 12 && useME12 ) correctME = false;
1871 0 : if ( meMode == 13 && !useME13 ) correctME = false;
1872 0 : if ( meMode != 13 && useME13 ) correctME = false;
1873 :
1874 : // Check matrix element mode 21: tau -> nu_tau hadrons.
1875 0 : bool useME21 = (idNow == 15 && mult > 2 && prod[0] == 16
1876 0 : && abs(prod[1]) > 100);
1877 0 : if ( meMode == 21 && !useME21 ) correctME = false;
1878 0 : if ( meMode != 21 && useME21 ) correctME = false;
1879 :
1880 : // Check matrix element mode 22, but only for semileptonic decay.
1881 : // For a -> b c d require types u = 2, ubar = -2, d = 1, dbar = -1.
1882 0 : if ( isLepton(prod[0]) && isLepton(prod[1]) ) {
1883 : bool useME22 = false;
1884 : int typeA = 0;
1885 : int typeB = 0;
1886 : int typeC = 0;
1887 0 : if (particlePtr->isLepton()) {
1888 0 : typeA = (idNow > 0) ? 1 + (idNow-1)%2 : -1 - (1-idNow)%2;
1889 0 : } else if (particlePtr->isHadron()) {
1890 0 : int hQ = particlePtr->heaviestQuark();
1891 : // Special case: for B_c either bbar or c decays.
1892 0 : if (idNow == 541 && heaviestQuark(prod[2]) == -5) hQ = 4;
1893 0 : typeA = (hQ > 0) ? 1 + (hQ-1)%2 : -1 - (1-hQ)%2;
1894 0 : }
1895 0 : typeB = (prod[0] > 0) ? 1 + (prod[0]-1)%2 : -1 - (1-prod[0])%2;
1896 0 : typeC = (prod[1] > 0) ? 1 + (prod[1]-1)%2 : -1 - (1-prod[1])%2;
1897 : // Special cases.
1898 0 : if ( (idNow == 130 || idNow == 310) && typeC * typeA < 0)
1899 0 : typeA = -typeA;
1900 0 : if (mult == 3 && idNow == 2112 && prod[2] == 2212)
1901 0 : typeA = 3 - typeA;
1902 0 : if (mult == 3 && idNow == 3222 && prod[2] == 3122)
1903 0 : typeA = 3 - typeA;
1904 0 : if (mult > 2 && typeC == typeA && typeB * typeC < 0
1905 0 : && (typeB + typeC)%2 != 0) useME22 = true;
1906 0 : if ( meMode == 22 && !useME22 ) correctME = false;
1907 0 : if ( meMode != 22 && useME22 ) correctME = false;
1908 0 : }
1909 :
1910 : // Check for matrix element mode 31.
1911 0 : if (meMode == 31) {
1912 : int nGamma = 0;
1913 0 : for (int j = 0; j < mult; ++j) if (prod[j] == 22) ++nGamma;
1914 0 : if (nGamma != 1) correctME = false;
1915 0 : }
1916 :
1917 : // Check for unknown mode, or resonance-only mode.
1918 0 : if ( !isResonanceNow && (meMode < 0 || (meMode > 2 && meMode <= 10)
1919 0 : || (meMode > 13 && meMode <= 20) || (meMode > 23 && meMode <= 30)
1920 0 : || (meMode > 31 && meMode <= 41) || meMode == 51 || meMode == 61
1921 0 : || meMode == 71 || (meMode > 80 && meMode <= 90)
1922 0 : || (!particlePtr->isOctetHadron() && meMode > 92) ) )
1923 0 : correctME = false;
1924 :
1925 : // Print if incorrect matrix element mode.
1926 0 : if ( !correctME ) {
1927 0 : os << " Warning: meMode " << meMode << " used for "
1928 0 : << idNow << " -> ";
1929 0 : for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1930 0 : os << "\n";
1931 : hasPrinted = true;
1932 0 : ++nErr;
1933 0 : }
1934 :
1935 : // Warning printout when no phase space for decay.
1936 0 : if ( studyCloser && verbosity > 0 && canHandle && onMode > 0
1937 0 : && particlePtr->m0Min() - minFinalMass < 0. ) {
1938 0 : if (particlePtr->m0Max() - minFinalMass < 0.)
1939 0 : os << " Error: decay never possible for ";
1940 0 : else os << " Warning: decay sometimes not possible for ";
1941 0 : os << idNow << " -> ";
1942 0 : for (int j = 0; j < mult; ++j) os << prod[j] << " ";
1943 0 : os << "\n";
1944 : hasPrinted = true;
1945 0 : ++nErr;
1946 0 : }
1947 :
1948 : // End loop over decay channels.
1949 0 : if (onMode > 0 && bRatio > 0.) openChannel = true;
1950 0 : }
1951 :
1952 : // Optional printout of threshold.
1953 0 : if (verbosity%10 > 1 && particlePtr->useBreitWigner()) {
1954 0 : threshMass /= bRatioSum;
1955 0 : os << " Info: particle " << idNow << fixed << setprecision(5)
1956 0 : << " has average mass threshold " << threshMass
1957 0 : << " while mMin is " << mMinNow << "\n";
1958 : hasPrinted = true;
1959 0 : }
1960 :
1961 : // Error printout when no acceptable decay channels found.
1962 0 : if (studyCloser && !openChannel) {
1963 0 : os << " Error: no acceptable decay channel found for particle "
1964 0 : << idNow << "\n";
1965 : hasPrinted = true;
1966 0 : ++nErr;
1967 0 : }
1968 :
1969 : // Warning printout when branching ratios do not sum to unity.
1970 0 : if (studyCloser && (!hasAntiNow || (!hasPosBR && !hasNegBR))
1971 0 : && abs(bRatioSum + bRatioPos - 1.) > 1e-8) {
1972 0 : os << " Warning: particle " << idNow << fixed << setprecision(8)
1973 0 : << " has branching ratio sum " << bRatioSum << "\n";
1974 : hasPrinted = true;
1975 0 : ++nErr;
1976 0 : } else if (studyCloser && hasAntiNow
1977 0 : && (abs(bRatioSum + bRatioPos - 1.) > 1e-8
1978 0 : || abs(bRatioSum + bRatioNeg - 1.) > 1e-8)) {
1979 0 : os << " Warning: particle " << idNow << fixed << setprecision(8)
1980 0 : << " has branching ratio sum " << bRatioSum + bRatioPos
1981 0 : << " while its antiparticle has " << bRatioSum + bRatioNeg
1982 0 : << "\n";
1983 : hasPrinted = true;
1984 0 : ++nErr;
1985 0 : }
1986 :
1987 : // End study of decay channels and loop over particles.
1988 0 : }
1989 0 : if (hasPrinted) os << "\n";
1990 0 : }
1991 :
1992 : // Final output. Done.
1993 0 : os << " Total number of errors and warnings is " << nErr << "\n";
1994 0 : os << "\n -------- End PYTHIA Check of Particle Data Table --------"
1995 0 : << "------\n" << endl;
1996 :
1997 0 : }
1998 :
1999 : //--------------------------------------------------------------------------
2000 :
2001 : // Return the id of the sequentially next particle stored in table.
2002 :
2003 : int ParticleData::nextId(int idIn) {
2004 :
2005 : // Return 0 for negative or unknown codes. Return first for 0.
2006 0 : if (idIn < 0 || (idIn > 0 && !isParticle(idIn))) return 0;
2007 0 : if (idIn == 0) return pdt.begin()->first;
2008 :
2009 : // Find pointer to current particle and step up. Return 0 if impossible.
2010 0 : map<int, ParticleDataEntry>::const_iterator pdtIn = pdt.find(idIn);
2011 0 : if (pdtIn == pdt.end()) return 0;
2012 0 : ++pdtIn;
2013 0 : if (pdtIn == pdt.end()) return 0;
2014 0 : return pdtIn->first;
2015 :
2016 0 : }
2017 :
2018 : //--------------------------------------------------------------------------
2019 :
2020 : // Fractional width associated with open channels of one or two resonances.
2021 :
2022 : double ParticleData::resOpenFrac(int id1In, int id2In, int id3In) {
2023 :
2024 : // Default value.
2025 : double answer = 1.;
2026 :
2027 : // First resonance.
2028 0 : if (isParticle(id1In)) answer = pdt[abs(id1In)].resOpenFrac(id1In);
2029 :
2030 : // Possibly second resonance.
2031 0 : if (isParticle(id2In)) answer *= pdt[abs(id2In)].resOpenFrac(id2In);
2032 :
2033 : // Possibly third resonance.
2034 0 : if (isParticle(id3In)) answer *= pdt[abs(id2In)].resOpenFrac(id3In);
2035 :
2036 : // Done.
2037 0 : return answer;
2038 :
2039 : }
2040 :
2041 : //--------------------------------------------------------------------------
2042 :
2043 : // Extract XML value string following XML attribute.
2044 :
2045 : string ParticleData::attributeValue(string line, string attribute) {
2046 :
2047 0 : if (line.find(attribute) == string::npos) return "";
2048 0 : int iBegAttri = line.find(attribute);
2049 0 : int iBegQuote = line.find("\"", iBegAttri + 1);
2050 0 : int iEndQuote = line.find("\"", iBegQuote + 1);
2051 0 : return line.substr(iBegQuote + 1, iEndQuote - iBegQuote - 1);
2052 :
2053 0 : }
2054 :
2055 : //--------------------------------------------------------------------------
2056 :
2057 : // Extract XML bool value following XML attribute.
2058 :
2059 : bool ParticleData::boolAttributeValue(string line, string attribute) {
2060 :
2061 0 : string valString = attributeValue(line, attribute);
2062 0 : if (valString == "") return false;
2063 0 : return boolString(valString);
2064 :
2065 0 : }
2066 :
2067 : //--------------------------------------------------------------------------
2068 :
2069 : // Extract XML int value following XML attribute.
2070 :
2071 : int ParticleData::intAttributeValue(string line, string attribute) {
2072 0 : string valString = attributeValue(line, attribute);
2073 0 : if (valString == "") return 0;
2074 0 : istringstream valStream(valString);
2075 0 : int intVal;
2076 0 : valStream >> intVal;
2077 0 : return intVal;
2078 :
2079 0 : }
2080 :
2081 : //--------------------------------------------------------------------------
2082 :
2083 : // Extract XML double value following XML attribute.
2084 :
2085 : double ParticleData::doubleAttributeValue(string line, string attribute) {
2086 0 : string valString = attributeValue(line, attribute);
2087 0 : if (valString == "") return 0.;
2088 0 : istringstream valStream(valString);
2089 0 : double doubleVal;
2090 0 : valStream >> doubleVal;
2091 0 : return doubleVal;
2092 :
2093 0 : }
2094 :
2095 : //==========================================================================
2096 :
2097 : } // end namespace Pythia8
|
