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