7270477f84ded227c869f646e351bd49fd6b264b
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / DatabasePDG.cxx
1 // DatabasePDG stores and handles PDG information
2 // The PDG particle definitions and decay channels are read
3 // in the begining from ASCII files
4 // PDG definitions loaded can be selected according to their
5 // mass and decay width
6
7 /*
8   Copyright   : The FASTMC and SPHMC Collaboration
9   Author      : Ionut Cristian Arsene 
10   Affiliation : Oslo University, Norway & Institute for Space Sciences, Bucharest, Romania
11   e-mail      : i.c.arsene@fys.uio.no
12   Date        : 2007/05/30
13
14   This class is using the particle and decay lists provided by the 
15   THERMINATOR (Computer Physics Communications 174 669 (2006)) and
16   SHARE (Computer Physics Communications 167 229 (2005)) collaborations.
17 */
18
19
20 #ifndef DATABASEPDG_H
21 #include "DatabasePDG.h"
22 #endif
23
24 #include <iostream>
25 #include <cstring>
26 #include <fstream>
27 using std::cout;
28 using std::endl;
29 using std::strcmp;
30 using std::strcpy;
31 using std::ifstream;
32 using std::ios;
33
34 DatabasePDG::DatabasePDG():
35   fNParticles(0),
36   fUseCharmParticles(kTRUE),
37   fMinimumWidth(0.0),
38   fMaximumWidth(10.),
39   fMinimumMass(0.0001),
40   fMaximumMass(10.)
41 {
42   // Default constructor, initialize members, set input files  
43   
44   strncpy(fParticleFilename, "particles.data", 256);
45   strncpy(fDecayFilename, "tabledecay.txt", 256);
46   for(Int_t i=0; i<kMaxParticles; i++) {
47     fParticles[i] = new ParticlePDG();
48     fStatus[i] = kFALSE;
49   }
50 }
51
52 DatabasePDG::~DatabasePDG() {
53   for(Int_t i=0; i<kMaxParticles; i++)
54     if(fParticles[i]) 
55       delete fParticles[i];
56 }
57
58 void DatabasePDG::SetParticleFilename(Char_t *filename) {
59   if(strlen(filename)<255)
60     strncpy(fParticleFilename, filename, 256);
61 }
62
63 void DatabasePDG::SetDecayFilename(Char_t *filename) {
64   if(strlen(filename)<255)
65     strncpy(fDecayFilename, filename, 256);
66 }
67
68 Bool_t DatabasePDG::LoadData() {
69   return (LoadParticles() && LoadDecays());
70 }
71
72 Bool_t DatabasePDG::LoadParticles() {
73   // Read particle definitions from the ascii file  
74
75   ifstream particleFile;
76   particleFile.open(fParticleFilename);
77   if(!particleFile) {
78     cout << "ERROR in DatabasePDG::LoadParticles() : The ASCII file containing the PDG particle list (\""
79          << fParticleFilename << "\") was not found !! Exiting..." << endl;
80     return kFALSE;
81   }
82   
83   Char_t name[9];
84   Double_t mass, width, spin, isospin, isospinZ, q, s, aq, as, c, ac;
85   Int_t pdg;
86   Int_t goodStatusParticles = 0;
87
88   cout << "Info in DatabasePDG::LoadParticles() : Start loading particles with the following criteria:" << endl
89        << "       Use particles containing charm quarks (1-yes;0-no) : " << fUseCharmParticles << endl
90        << "       Mass range                                         : (" << fMinimumMass << "; " << fMaximumMass << ")" << endl
91        << "       Width range                                        : (" << fMinimumWidth << "; " << fMaximumWidth << ")" << endl;
92   
93   particleFile.exceptions(ios::failbit);
94   while(!particleFile.eof()) {
95     try {
96       particleFile >> name >> mass >> width >> spin >> isospin >> isospinZ >> q >> s >> aq >> as >> c >> ac >> pdg;
97     }
98     catch (ios::failure const &problem) {
99       cout << problem.what() << endl;
100       break;
101     }
102         
103     fParticles[fNParticles]->SetName(name);
104     fParticles[fNParticles]->SetPDG(pdg);
105     fParticles[fNParticles]->SetMass(mass);
106     fParticles[fNParticles]->SetWidth(width);
107     fParticles[fNParticles]->SetSpin(spin);
108     fParticles[fNParticles]->SetIsospin(isospin);
109     fParticles[fNParticles]->SetIsospinZ(isospinZ);
110     fParticles[fNParticles]->SetLightQNumber(q);
111     fParticles[fNParticles]->SetStrangeQNumber(s);
112     fParticles[fNParticles]->SetLightAQNumber(aq);
113     fParticles[fNParticles]->SetStrangeAQNumber(as);
114     fParticles[fNParticles]->SetCharmQNumber(c);
115     fParticles[fNParticles]->SetCharmAQNumber(ac);
116     
117     fStatus[fNParticles] = kTRUE;
118     // check if we want charmed particles
119     if(!fUseCharmParticles && (c>0 || ac>0)) {
120       fStatus[fNParticles] = kFALSE;
121     }
122     // check that the particle mass is inside accepted limits
123     if(!(fMinimumMass<=mass && mass<=fMaximumMass)) {
124       fStatus[fNParticles] = kFALSE;
125     }
126     // check that the particle width is inside accepted limits
127     if(!(fMinimumWidth<=width && width<=fMaximumWidth)) {
128       fStatus[fNParticles] = kFALSE;
129     }
130     if(fStatus[fNParticles]) goodStatusParticles++;
131     fNParticles++;
132   }
133   particleFile.close();
134   if(fNParticles==0) {
135     cout << "Warning in DatabasePDG::LoadParticles(): No particles were found in the file specified!!" << endl;
136     return kFALSE;
137   }
138   SortParticles();
139   cout << "Info in DatabasePDG::LoadParticles(): Particle definitions found = " << fNParticles << endl
140        << "                                      Good status particles      = " << goodStatusParticles << endl;
141   return kTRUE;
142 }
143
144 Bool_t DatabasePDG::LoadDecays() {
145   // Read the decay channel definitions from the ascii file
146
147   ifstream decayFile;
148   decayFile.open(fDecayFilename);
149   if(!decayFile) {
150     cout << "ERROR in DatabasePDG::LoadDecays() : The ASCII file containing the decays list (\""
151          << fDecayFilename << "\") was not found !! Exiting..." << endl;
152     return kFALSE;
153   }
154   
155   Int_t motherPdg, daughterPdg[3];
156   Double_t branching;
157   
158   decayFile.exceptions(ios::failbit);
159   while(!decayFile.eof()) {
160     motherPdg = 0;
161     for(Int_t i=0; i<3; i++) daughterPdg[i] = 0;
162     branching = -1.0;
163     try {
164       decayFile >> motherPdg;
165       for(Int_t i=0; i<3; i++) 
166         decayFile >> daughterPdg[i];
167       decayFile >> branching;
168     }
169     catch (ios::failure const &problem) {
170       cout << problem.what() << endl;
171       break;
172     }
173     if((motherPdg!=0) && (daughterPdg[0]!=0) && (branching>=0)) {
174       Int_t nDaughters = 0;
175       for(Int_t i=0; i<3; i++)
176         if(daughterPdg[i]!=0)
177           nDaughters++;
178       ParticlePDG* particle = GetPDGParticle(motherPdg);
179       DecayChannel decay(motherPdg, branching, nDaughters, daughterPdg);
180       particle->AddChannel(&decay);
181     }
182   }
183   decayFile.close();
184   Int_t nDecayChannels = 0;
185   for(Int_t i=0; i<fNParticles; i++) {
186     nDecayChannels += fParticles[i]->GetNDecayChannels();
187   }
188   cout << "Info in DatabasePDG::LoadDecays(): Number of decays found in the database is " << nDecayChannels << endl;
189   return kTRUE;
190 }
191
192 ParticlePDG* DatabasePDG::GetPDGParticleByIndex(Int_t index) const {
193   // Return a PDG particle definition based on its index in the particle list 
194
195   if(index<0 || index>fNParticles) {
196     cout << "Warning in DatabasePDG::GetPDGParticleByIndex(Int_t): Particle index is negative or too big !!" << endl
197          << " It must be inside this range: [0, " << fNParticles-1 << "]" << endl
198          << " Returning null pointer!!" << endl;
199     return 0x0;
200   }
201   return fParticles[index];
202 }
203
204 Bool_t DatabasePDG::GetPDGParticleStatusByIndex(Int_t index) const {
205   // Return the status of a PDG particle definition based on its index in the particle list
206   // The status is kTRUE when a particle passed the mass, width and charm criteria
207   // and kFALSE otherwise
208
209   if(index<0 || index>fNParticles) {
210     cout << "Warning in DatabasePDG::GetPDGParticleStatusByIndex(Int_t): Particle index is negative or too big !!" << endl
211          << " It must be inside this range: [0, " << fNParticles-1 << "]" << endl
212          << " Returning null pointer!!" << endl;
213     return kFALSE;
214   }
215   return fStatus[index];
216 }
217
218 ParticlePDG* DatabasePDG::GetPDGParticle(Int_t pdg) const {
219   // Return a PDG particle definition based on the PDG code (PYTHIA convention) 
220   // If more than 1 definition with the asked PDG code is found then a warning is issued 
221
222   Int_t nFindings = 0;
223   Int_t firstTimeIndex = 0;
224   for(Int_t i=0; i<fNParticles; i++) {
225     if(pdg == fParticles[i]->GetPDG()) {
226       if(nFindings == 0) firstTimeIndex = i;
227       nFindings++;
228     }
229   }
230   if(nFindings == 1) return fParticles[firstTimeIndex];
231   if(nFindings == 0) {
232     //cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
233     //     << " was not found in the database!!" << endl;
234     return 0x0;
235   }
236   if(nFindings >= 2) {
237     cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
238          << " was found with " << nFindings << " entries in the database. Check it out !!" << endl
239          << "Returning the first instance found" << endl;
240     return fParticles[firstTimeIndex];
241   }
242   return 0x0;
243 }
244
245 Bool_t DatabasePDG::GetPDGParticleStatus(Int_t pdg) const {
246   // Return the status of a PDG particle definition based on its PDG code (PYTHIA convention)
247   // The status is kTRUE when a particle passed the mass, width and charm criteria
248   // and kFALSE otherwise
249
250   Int_t nFindings = 0;
251   Int_t firstTimeIndex = 0;
252   for(Int_t i=0; i<fNParticles; i++) {
253     if(pdg == fParticles[i]->GetPDG()) {
254       if(nFindings == 0) firstTimeIndex = i;
255       nFindings++;
256     }
257   }
258   if(nFindings == 1) return fStatus[firstTimeIndex];
259   if(nFindings == 0) {
260     //cout << "Warning in DatabasePDG::GetPDGParticle(Int_t): The particle required with PDG = " << pdg
261     //     << " was not found in the database!!" << endl;
262     return kFALSE;
263   }
264   if(nFindings >= 2) {
265     cout << "Warning in DatabasePDG::GetPDGParticleStatus(Int_t): The particle status required for PDG = " << pdg
266          << " was found with " << nFindings << " entries in the database. Check it out !!" << endl
267          << "Returning the status of first instance found" << endl;
268     return fStatus[firstTimeIndex];
269   }
270   return kFALSE;
271 }
272
273 ParticlePDG* DatabasePDG::GetPDGParticle(Char_t* name) const {
274   // Return a PDG particle definition based on its name
275   // If more than 1 definition with the asked PDG code is found then a warning is issued
276
277   Int_t nFindings = 0;
278   Int_t firstTimeIndex = 0;
279   for(Int_t i=0; i<fNParticles; i++) {
280     if(!strcmp(name, fParticles[i]->GetName())) {
281       if(nFindings == 0) firstTimeIndex = i;
282       nFindings++;
283     }
284   }
285   if(nFindings == 1) return fParticles[firstTimeIndex];
286   if(nFindings == 0) {
287     //cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
288     //     << "\" was not found in the database!!" << endl;
289     return 0x0;
290   }
291   if(nFindings >= 2) {
292     cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
293          << "\" was found with " << nFindings << " entries in the database. Check it out !!" << endl
294          << "Returning the first instance found" << endl;
295     return fParticles[firstTimeIndex];
296   }
297   return 0x0;
298 }
299
300 Bool_t DatabasePDG::GetPDGParticleStatus(Char_t* name) const {
301   // Return the status of a PDG particle definition based on its name
302   // The status is kTRUE when a particle passed the mass, width and charm criteria
303   // and kFALSE otherwise
304
305   Int_t nFindings = 0;
306   Int_t firstTimeIndex = 0;
307   for(Int_t i=0; i<fNParticles; i++) {
308     if(!strcmp(name, fParticles[i]->GetName())) {
309       if(nFindings == 0) firstTimeIndex = i;
310       nFindings++;
311     }
312   }
313   if(nFindings == 1) return fStatus[firstTimeIndex];
314   if(nFindings == 0) {
315     //cout << "Warning in DatabasePDG::GetPDGParticle(Char_t*): The particle required with name \"" << name
316     //     << "\" was not found in the database!!" << endl;
317     return kFALSE;
318   }
319   if(nFindings >= 2) {
320     cout << "Warning in DatabasePDG::GetPDGParticleStatus(Char_t*): The particle status required for name \"" << name
321          << "\" was found with " << nFindings << " entries in the database. Check it out !!" << endl
322          << "Returning the first instance found" << endl;
323     return fStatus[firstTimeIndex];
324   }
325   return kFALSE;
326 }
327
328 void DatabasePDG::DumpData(Bool_t dumpAll) const {
329   // Printout all the information in the PDG database
330
331   cout << "***********************************************************************************************************" << endl;
332   cout << "Dumping all the information contained in the database..." << endl;
333   Int_t nDecays = 0;
334   Int_t nGoodStatusDecays = 0;
335   Int_t nGoodStatusParticles = 0;
336   for(Int_t currPart=0; currPart<fNParticles; currPart++) {
337     nGoodStatusParticles += (fStatus[currPart] ? 1:0);
338     nGoodStatusDecays += (fStatus[currPart] ? fParticles[currPart]->GetNDecayChannels() : 0);
339     nDecays += fParticles[currPart]->GetNDecayChannels();
340     if(!(dumpAll || (!dumpAll && fStatus[currPart]))) continue;
341     cout << "###### Particle: " << fParticles[currPart]->GetName() << " with PDG code " << fParticles[currPart]->GetPDG() << endl;
342     cout << "   status          = " << fStatus[currPart] << endl;
343     cout << "   mass            = " << fParticles[currPart]->GetMass() << " GeV" << endl;
344     cout << "   width           = " << fParticles[currPart]->GetWidth() << " GeV" << endl;
345     cout << "   2*spin          = " << Int_t(2.*fParticles[currPart]->GetSpin()) << endl;
346     cout << "   2*isospin       = " << Int_t(2.*fParticles[currPart]->GetIsospin()) << endl;
347     cout << "   2*isospin3      = " << Int_t(2.*fParticles[currPart]->GetIsospinZ()) << endl;
348     cout << "   u,d quarks      = " << Int_t(fParticles[currPart]->GetLightQNumber()) << endl;
349     cout << "   s quarks        = " << Int_t(fParticles[currPart]->GetStrangeQNumber()) << endl;
350     cout << "   c quarks        = " << Int_t(fParticles[currPart]->GetCharmQNumber()) << endl;
351     cout << "   anti u,d quarks = " << Int_t(fParticles[currPart]->GetLightAQNumber()) << endl;
352     cout << "   anti s quarks   = " << Int_t(fParticles[currPart]->GetStrangeAQNumber()) << endl;
353     cout << "   anti c quarks   = " << Int_t(fParticles[currPart]->GetCharmAQNumber()) << endl;
354     cout << "   baryon number   = " << Int_t(fParticles[currPart]->GetBaryonNumber()) << endl;
355     cout << "   strangeness     = " << Int_t(fParticles[currPart]->GetStrangeness()) << endl;
356     cout << "   charmness       = " << Int_t(fParticles[currPart]->GetCharmness()) << endl;
357     cout << "   electric charge = " << Int_t(fParticles[currPart]->GetElectricCharge()) << endl;
358     cout << "   full branching  = " << fParticles[currPart]->GetFullBranching() << endl;
359     cout << "   decay modes     = " << fParticles[currPart]->GetNDecayChannels() << endl;
360     for(Int_t currChannel=0; currChannel<fParticles[currPart]->GetNDecayChannels(); currChannel++) {
361       cout << "   channel " << currChannel+1 << " with branching " << fParticles[currPart]->GetDecayChannel(currChannel)->GetBranching() << endl;
362       cout << "   daughters PDG codes: ";
363       Double_t daughtersMass = 0.0;
364       for(Int_t currDaughter=0; currDaughter<fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters(); currDaughter++) {
365         cout << fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter) << "\t";
366         ParticlePDG *daughter = GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
367         daughtersMass += daughter->GetMass();
368       }
369       cout << endl;
370       cout << "   daughters sum mass = " << daughtersMass << endl;
371     }
372   }
373   if(dumpAll) {
374     cout << "Finished dumping information for " << fNParticles << " particles with " << nDecays << " decay channels in total." << endl;
375     cout << "*************************************************************************************************************" << endl;
376   }
377   else {
378     cout << "Finished dumping information for " << nGoodStatusParticles << "(" << fNParticles << ")" 
379          << " particles with " << nGoodStatusDecays << "(" << nDecays << ")" << " decay channels in total." << endl;
380     cout << "*************************************************************************************************************" << endl;
381   }
382 }
383
384 Int_t DatabasePDG::CheckImpossibleDecays(Bool_t dump) const {
385   // Check the database for impossible decays
386
387   Int_t nImpossibleDecays = 0;
388   for(Int_t currPart=0; currPart<fNParticles; currPart++) {
389     if(!fStatus[currPart]) continue;
390     Int_t allChannels = fParticles[currPart]->GetNDecayChannels();
391     Int_t allowedChannels = GetNAllowedChannels(fParticles[currPart], fParticles[currPart]->GetMass());
392     if(dump) {
393       cout << "Particle " << fParticles[currPart]->GetPDG() << " has " << allChannels << " decay channels specified in the database" << endl;
394       cout << " Allowed channels assuming table mass = " << allowedChannels << endl;
395     }
396     if(dump && allChannels>0 && allowedChannels == 0) {
397       cout << "**********************************************************************" << endl;
398       cout << "       All channels for this particles are not allowed" << endl;
399       cout << "**********************************************************************" << endl;
400     }
401     if(dump && fParticles[currPart]->GetWidth() > 0. && allChannels == 0) {
402       cout << "**********************************************************************" << endl;
403       cout << "    Particle has finite width but no decay channels specified" << endl;
404       cout << "**********************************************************************" << endl;
405     }
406     for(Int_t currChannel=0; currChannel<fParticles[currPart]->GetNDecayChannels(); currChannel++) {
407       Double_t motherMass = fParticles[currPart]->GetMass();
408       Double_t daughtersSumMass = 0.;
409       for(Int_t currDaughter=0; currDaughter<fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters(); currDaughter++) {
410         ParticlePDG *daughter = GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
411         daughtersSumMass += daughter->GetMass();
412       }
413       if(daughtersSumMass >= motherMass) {
414         nImpossibleDecays++;
415         if(dump) {
416           cout << "Imposible decay for particle " << fParticles[currPart]->GetPDG() << endl;
417           cout << "  Channel: " << fParticles[currPart]->GetPDG() << " --> ";
418           for(Int_t currDaughter=0; currDaughter<fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters(); currDaughter++) {
419             ParticlePDG *daughter = GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
420             cout << daughter->GetPDG() << " ";
421           }
422           cout << endl;
423           cout << "  Mother particle mass = " << motherMass << endl;
424           cout << "  Daughters sum mass   = " << daughtersSumMass << endl;
425         }
426       }
427     }
428   }
429   return nImpossibleDecays;
430 }
431
432 void DatabasePDG::SetUseCharmParticles(Bool_t flag) {
433   // Switch on/off the use of charmed particles
434
435   if(fNParticles>0) {
436     fUseCharmParticles = flag;
437     for(Int_t i=0; i<fNParticles; i++) {
438       if(fParticles[i]->GetCharmQNumber()>0 || fParticles[i]->GetCharmAQNumber())                 
439         fStatus[i] = flag;
440     }
441     SortParticles();
442     return;
443   }
444   else
445     fUseCharmParticles = flag;
446   return;
447 }
448
449 void DatabasePDG::SetMinimumWidth(Double_t value) {
450   // Set the minimum decay width for the particle definitions to be generated in the soft fireball
451   if(fNParticles>0) {
452     fMinimumWidth = value;
453     for(Int_t i=0; i<fNParticles; i++) {
454       if(fParticles[i]->GetWidth() < fMinimumWidth)               
455         fStatus[i] = kFALSE;
456     }
457     SortParticles();
458     return;
459   }
460   else
461     fMinimumWidth = value;
462   return;
463 }
464
465 void DatabasePDG::SetMaximumWidth(Double_t value) {
466   // Set the maximum decay width for the particle definitions to be generated in the soft fireball
467
468   if(fNParticles>0) {
469     fMaximumWidth = value;
470     for(Int_t i=0; i<fNParticles; i++) {
471       if(fParticles[i]->GetWidth() > fMaximumWidth)               
472         fStatus[i] = kFALSE;
473     }
474     SortParticles();
475     return;
476   }
477   else
478     fMaximumWidth = value;
479   return;
480 }
481
482 void DatabasePDG::SetWidthRange(Double_t min, Double_t max) {
483   // Set the decay width range for the particle definitions to be generated in the soft fireball
484
485   if(fNParticles>0) {
486     fMinimumWidth = min;
487     fMaximumWidth = max;
488     for(Int_t i=0; i<fNParticles; i++) {
489       if((fParticles[i]->GetWidth()<fMinimumWidth) || (fParticles[i]->GetWidth()>fMaximumWidth))  
490         fStatus[i] = kFALSE;
491     }
492     SortParticles();
493     return;
494   }
495   else {
496     fMinimumWidth = min;
497     fMaximumWidth = max;
498   }
499   return;
500 }
501
502 void DatabasePDG::SetMinimumMass(Double_t value) {
503   // Set the minimum mass for the particle definitions to be generated in the soft fireball
504
505   if(fNParticles>0) {
506     fMinimumMass = value;
507     for(Int_t i=0; i<fNParticles; i++) {
508       if(fParticles[i]->GetMass() < fMinimumMass)                 
509         fStatus[i] = kFALSE;
510     }
511     SortParticles();
512     return;
513   }
514   else
515     fMinimumMass = value;
516   return;
517 }
518
519 void DatabasePDG::SetMaximumMass(Double_t value) {
520   // Set the maximum mass for the particle definitions to be generated in the soft fireball
521
522   if(fNParticles>0) {
523     fMaximumMass = value;
524     for(Int_t i=0; i<fNParticles; i++) {
525       if(fParticles[i]->GetMass() > fMaximumMass)                 
526         fStatus[i] = kFALSE;
527     }
528     SortParticles();
529     return;
530   }
531   else
532     fMaximumMass = value;
533   return;
534 }
535
536 void DatabasePDG::SetMassRange(Double_t min, Double_t max) {
537   // Set the mass range for the particle definitions to be generated in the soft fireball
538
539   if(fNParticles>0) {
540     fMinimumMass = min;
541     fMaximumMass = max;
542     for(Int_t i=0; i<fNParticles; i++) {
543       if((fParticles[i]->GetMass()<fMinimumMass) || (fParticles[i]->GetMass()>fMaximumMass))  
544         fStatus[i] = kFALSE;
545     }
546     SortParticles();
547     return;
548   }
549   else {
550     fMinimumMass = min;
551     fMaximumMass = max;
552   }
553   return;
554 }
555
556 void DatabasePDG::SortParticles() {
557   // Sort the particle list so that those with kTRUE status will be always on top of the list
558
559   if(fNParticles<2) {
560     cout << "Warning in DatabasePDG::SortParticles() : No particles to sort. Load data first!!" << endl;
561     return;
562   }
563
564   Int_t nGoodStatus = 0;
565   for(Int_t i=0; i<fNParticles; i++)
566     if(fStatus[i]) nGoodStatus++;
567   if(nGoodStatus==fNParticles)    // if all particles have good status then there is nothing to do
568     return;
569   if(nGoodStatus==0)              // no good status particles, again nothing to do
570     return;
571
572   Int_t shifts = 1;
573   while(shifts) {
574     shifts = 0;
575     for(Int_t i=0; i<fNParticles-1; i++) {
576       if(!fStatus[i] && fStatus[i+1]) {   // switch if false status is imediately before a true status particle
577         ParticlePDG *temporaryPointer = fParticles[i];
578         fParticles[i] = fParticles[i+1];
579         fParticles[i+1] = temporaryPointer;
580         Bool_t temporaryStatus = fStatus[i];
581         fStatus[i] = fStatus[i+1];
582         fStatus[i+1] = temporaryStatus;
583         shifts++;
584       }
585     }
586   }
587   return;
588 }
589
590 Int_t DatabasePDG::GetNParticles(Bool_t all) const {
591   // Return the number of particle definitions in the database
592   // If all is kTRUE then return number of all particle
593   // If all is kFALSE then return the number of good status (kTRUE) particles
594
595   if(all)
596     return fNParticles;
597
598   Int_t nGoodStatus = 0;
599   for(Int_t i=0; i<fNParticles; i++)
600     if(fStatus[i]) nGoodStatus++;
601   return nGoodStatus;
602 }
603
604 void DatabasePDG::UseThisListOfParticles(Char_t *filename, Bool_t exclusive) {
605   // Read a list of PDG codes from the file "filename" and mark them with good status (kTRUE)
606   // while all the other will be marked kFALSE (only if exclusive = kTRUE)     
607
608   if(fNParticles<1) {
609     cout << "Error in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : You must load the data before calling this function!!" << endl;
610     return;
611   }
612
613   ifstream listFile;
614   listFile.open(filename);
615   if(!listFile) {
616     cout << "ERROR in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The ASCII file containing the PDG codes list (\""
617          << filename << "\") was not found !! Exiting..." << endl;
618     return;
619   }
620
621   Bool_t flaggedIndexes[kMaxParticles];
622   for(Int_t i=0; i<kMaxParticles; i++)
623     flaggedIndexes[i] = kFALSE;
624   Int_t pdg = 0;
625   listFile.exceptions(ios::failbit);
626   while(!listFile.eof()) {
627     try {
628       listFile >> pdg;
629     }
630     catch (ios::failure const &problem) {
631       cout << problem.what() << endl;
632       break;
633     }
634     Int_t found = 0;
635     for(Int_t i=0; i<fNParticles; i++) {
636       if(fParticles[i]->GetPDG()==pdg) {
637         found++;
638         flaggedIndexes[i] = kTRUE;
639       }
640     }
641     if(!found) {
642       cout << "Warning in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The particle with PDG code "
643            << pdg << " was asked but not found in the database!!" << endl;
644     }
645     if(found>1) {
646       cout << "Warning in DatabasePDG::UseThisListOfParticles(Char_t*, Bool_t) : The particle with PDG code "
647            << pdg << " was found more than once in the database!!" << endl;
648     }
649   }
650
651   if(exclusive) {
652     for(Int_t i=0; i<kMaxParticles; i++)
653       fStatus[i] = flaggedIndexes[i];
654   }
655   else {
656     for(Int_t i=0; i<kMaxParticles; i++)
657       fStatus[i] = (fStatus[i] && flaggedIndexes[i]);
658   }
659   SortParticles();
660
661   return;
662 }
663
664 Bool_t DatabasePDG::IsChannelAllowed(DecayChannel *channel, Double_t motherMass) const {
665   // Check if the decay channel "channel" is allowed by using the mother particle mass "motherMass"
666
667   Double_t daughtersSumMass = 0.0;
668   for(Int_t i=0; i<channel->GetNDaughters(); i++)
669     daughtersSumMass += GetPDGParticle(channel->GetDaughterPDG(i))->GetMass();
670   if(daughtersSumMass<=motherMass)
671     return kTRUE;
672   return kFALSE;
673 }
674
675 Int_t DatabasePDG::GetNAllowedChannels(ParticlePDG *particle, Double_t motherMass) const {
676   // Check how many decay channels are allowed for a given particle definition at a given mass
677
678   Int_t nAllowedChannels = 0;
679   for(Int_t i=0; i<particle->GetNDecayChannels(); i++) 
680     nAllowedChannels += (IsChannelAllowed(particle->GetDecayChannel(i), motherMass) ? 1:0);
681   
682   return nAllowedChannels;
683 }