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