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