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, 256);
  60 }
  61
  62 void DatabasePDG::SetDecayFilename(Char_t *filename) {
  63   strncpy(fDecayFilename, filename, 256);
  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 }