TUHKMgen/UHKM/DatabasePDG.cxx

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