DatabasePDG.cc

Go to the documentation of this file.

/*
  Copyright   : The FASTMC and SPHMC Collaboration
  Author      : Ionut Cristian Arsene 
  Affiliation : Oslo University, Norway & Institute for Space Sciences, Bucharest, Romania
  e-mail      : i.c.arsene@fys.uio.no
  Date        : 2007/05/30
 
  This class is using the particle and decay lists provided by the 
  THERMINATOR (Computer Physics Communications 174 669 (2006)) and
  SHARE (Computer Physics Communications 167 229 (2005)) collaborations.
*/
 
#ifndef DATABASE_PDG
#include "GeneratorInterface/Hydjet2Interface/interface/DatabasePDG.h"
#endif
 
#include <cstring>
#include <iostream>
#include <string>
#include <fstream>
 
using namespace std;
using namespace edm;
 
const char *particlesDATAstr;
const char *tableDECAYstr;
 
DatabasePDG::DatabasePDG() {
  fNParticles = 0;
 
  std::string file1 = "GeneratorInterface/Hydjet2Interface/data/particles.data";
  edm::FileInPath f1(file1);
  particlesDATAstr = ((f1.fullPath()).c_str());
 
  std::string file2 = "GeneratorInterface/Hydjet2Interface/data/tabledecay.txt";
  edm::FileInPath f2(file2);
  tableDECAYstr = ((f2.fullPath()).c_str());
 
  strcpy(fParticleFilename, particlesDATAstr);
  strcpy(fDecayFilename, tableDECAYstr);
  for (int i = 0; i < kMaxParticles; i++) {
    fParticles[i] = new ParticlePDG();
    fStatus[i] = kFALSE;
  }
  fUseCharmParticles = kTRUE;
  fMinimumWidth = 0.;
  fMaximumWidth = 10.;
  fMinimumMass = 0.;
  fMaximumMass = 200.;
}
 
DatabasePDG::~DatabasePDG() {
  for (int i = 0; i < kMaxParticles; i++)
    delete fParticles[i];
}
 
void DatabasePDG::SetParticleFilename(char *filename) { strcpy(fParticleFilename, filename); }
 
void DatabasePDG::SetDecayFilename(char *filename) { strcpy(fDecayFilename, filename); }
 
bool DatabasePDG::LoadData() { return (LoadParticles() && LoadDecays()); }
 
bool DatabasePDG::LoadParticles() {
  ifstream particleFile;
  particleFile.open(fParticleFilename);
  if (!particleFile) {
    edm::LogError("DatabasePDG") << "The ASCII file containing the PDG particle list " << fParticleFilename
                                 << " was not found";
    return kFALSE;
  }
 
  char name[9];
  double mass, width, spin, isospin, isospinZ, q, s, aq, as, c, ac;
  int pdg;
  int goodStatusParticles = 0;
 
  edm::LogInfo("DatabasePDG") << "Start loading particles with the following criteria:" << endl
                              << "       Use particles containing charm quarks (1:yes;0:no) : " << fUseCharmParticles
                              << endl
                              << "       Mass range                                         : (" << fMinimumMass << "; "
                              << fMaximumMass << ")" << endl
                              << "       Width range                                        : (" << fMinimumWidth
                              << "; " << fMaximumWidth << ")";
 
  particleFile.exceptions(ios::failbit);
  while (!particleFile.eof()) {
    try {
      particleFile >> name >> mass >> width >> spin >> isospin >> isospinZ >> q >> s >> aq >> as >> c >> ac >> pdg;
    } catch (ios::failure const &problem) {
      LogDebug("DatabasePDG") << " ios:failure in particle file " << problem.what();
      break;
    }
 
    fParticles[fNParticles]->SetName(name);
    fParticles[fNParticles]->SetPDG(pdg);
    fParticles[fNParticles]->SetMass(mass);
    fParticles[fNParticles]->SetWidth(width);
    fParticles[fNParticles]->SetSpin(spin);
    fParticles[fNParticles]->SetIsospin(isospin);
    fParticles[fNParticles]->SetIsospinZ(isospinZ);
    fParticles[fNParticles]->SetLightQNumber(q);
    fParticles[fNParticles]->SetStrangeQNumber(s);
    fParticles[fNParticles]->SetLightAQNumber(aq);
    fParticles[fNParticles]->SetStrangeAQNumber(as);
    fParticles[fNParticles]->SetCharmQNumber(c);
    fParticles[fNParticles]->SetCharmAQNumber(ac);
    goodStatusParticles++;
    fStatus[fNParticles] = kTRUE;
    // check if we want charmed particles
    if (!fUseCharmParticles && (c > 0 || ac > 0)) {
      fStatus[fNParticles] = kFALSE;
      goodStatusParticles--;
    }
    // check that the particle mass is inside accepted limits
    if (!(fMinimumMass <= mass && mass <= fMaximumMass)) {
      fStatus[fNParticles] = kFALSE;
      goodStatusParticles--;
    }
    // check that the particle width is inside accepted limits
    if (!(fMinimumWidth <= width && width <= fMaximumWidth)) {
      fStatus[fNParticles] = kFALSE;
      goodStatusParticles--;
    }
 
    fNParticles++;
  }
  particleFile.close();
  if (fNParticles == 0) {
    LogWarning("DatabasePDG") << " No particles were found in the file specified!!";
    return kFALSE;
  }
  SortParticles();
  edm::LogInfo("DatabasePDG") << " Particle definitions found: " << fNParticles
                              << ". Good status particles: " << goodStatusParticles;
  return kTRUE;
}
 
bool DatabasePDG::LoadDecays() {
  ifstream decayFile;
  decayFile.open(fDecayFilename);
  if (!decayFile) {
    edm::LogError("DatabasePDG") << "The ASCII file containing the decays list " << fDecayFilename << " was not found";
    return kFALSE;
  }
 
  int mother_pdg, daughter_pdg[3];
  double branching;
 
  decayFile.exceptions(ios::failbit);
  while (!decayFile.eof()) {
    mother_pdg = 0;
    for (int i = 0; i < 3; i++)
      daughter_pdg[i] = 0;
    branching = -1.0;
    try {
      decayFile >> mother_pdg;
      for (int i = 0; i < 3; i++)
        decayFile >> daughter_pdg[i];
      decayFile >> branching;
    } catch (ios::failure const &problem) {
      LogDebug("DatabasePDG") << " ios:failure in decay file " << problem.what();
      break;
    }
    if ((mother_pdg != 0) && (daughter_pdg[0] != 0) && (branching >= 0)) {
      int nDaughters = 0;
      for (int i = 0; i < 3; i++)
        if (daughter_pdg[i] != 0)
          nDaughters++;
      ParticlePDG *particle = GetPDGParticle(mother_pdg);
      if (!particle) {
        LogWarning("DatabasePDG") << " Mother particle PDG (" << mother_pdg
                                  << ") not found in the particle definition list:" << mother_pdg << " >>> ";
        for (int kk = 0; kk < nDaughters; kk++)
          LogWarning("DatabasePDG") << daughter_pdg[kk] << "  ";
        return kFALSE;
      }
      for (int kk = 0; kk < nDaughters; kk++) {
        if (!GetPDGParticle(daughter_pdg[kk])) {
          LogWarning("DatabasePDG") << "Daughter particle PDG (" << daughter_pdg[kk]
                                    << ") not found in the particle definition list: " << mother_pdg << ">>> ";
          for (int kkk = 0; kkk < nDaughters; kkk++)
            LogWarning("DatabasePDG") << daughter_pdg[kkk] << "  ";
        }
      }
      DecayChannel decay(mother_pdg, branching, nDaughters, daughter_pdg);
      particle->AddChannel(decay);
    }
  }
  decayFile.close();
  int nDecayChannels = 0;
  for (int i = 0; i < fNParticles; i++) {
    nDecayChannels += fParticles[i]->GetNDecayChannels();
  }
  edm::LogInfo("DatabasePDG") << "Number of decays found in the database is " << nDecayChannels;
  return kTRUE;
}
 
ParticlePDG *DatabasePDG::GetPDGParticleByIndex(int index) {
  if (index < 0 || index > fNParticles) {
    edm::LogWarning("DatabasePDG") << "Particle index is negative or too big !!" << endl
                                   << " It must be inside this range: (0, " << fNParticles - 1 << ")" << endl
                                   << " Returning null pointer!!";
    return nullptr;
  }
  return fParticles[index];
}
 
bool DatabasePDG::GetPDGParticleStatusByIndex(int index) {
  if (index < 0 || index > fNParticles) {
    edm::LogWarning("DatabasePDG") << "Particle index is negative or too big !!" << endl
                                   << " It must be inside this range: (0, " << fNParticles - 1 << ")" << endl
                                   << " Returning null pointer!!";
    return kFALSE;
  }
  return fStatus[index];
}
 
ParticlePDG *DatabasePDG::GetPDGParticle(int pdg) {
  int nFindings = 0;
  int firstTimeIndex = 0;
  for (int i = 0; i < fNParticles; i++) {
    if (pdg == fParticles[i]->GetPDG()) {
      if (nFindings == 0)
        firstTimeIndex = i;
      nFindings++;
    }
  }
  if (nFindings == 1)
    return fParticles[firstTimeIndex];
  if (nFindings == 0) {
    edm::LogWarning("DatabasePDG") << "The particle required with PDG: " << pdg << " was not found in the database!!";
    return nullptr;
  }
  if (nFindings >= 2) {
    edm::LogWarning("DatabasePDG") << "The particle required with PDG: " << pdg << " was found with " << nFindings
                                   << " entries in the database. Check it out !!" << endl
                                   << "Returning the first instance found";
    return fParticles[firstTimeIndex];
  }
  return nullptr;
}
 
bool DatabasePDG::GetPDGParticleStatus(int pdg) {
  int nFindings = 0;
  int firstTimeIndex = 0;
  for (int i = 0; i < fNParticles; i++) {
    if (pdg == fParticles[i]->GetPDG()) {
      if (nFindings == 0)
        firstTimeIndex = i;
      nFindings++;
    }
  }
  if (nFindings == 1)
    return fStatus[firstTimeIndex];
  if (nFindings == 0) {
    edm::LogWarning("DatabasePDG") << "The particle required with PDG: " << pdg << " was not found in the database!!";
    return kFALSE;
  }
  if (nFindings >= 2) {
    edm::LogWarning("DatabasePDG") << "The particle status required for PDG: " << pdg << " was found with " << nFindings
                                   << " entries in the database. Check it out !!" << endl
                                   << "Returning the status of first instance found";
    return fStatus[firstTimeIndex];
  }
  return kFALSE;
}
 
ParticlePDG *DatabasePDG::GetPDGParticle(char *name) {
  int nFindings = 0;
  int firstTimeIndex = 0;
  for (int i = 0; i < fNParticles; i++) {
    if (!strcmp(name, fParticles[i]->GetName())) {
      if (nFindings == 0)
        firstTimeIndex = i;
      nFindings++;
    }
  }
  if (nFindings == 1)
    return fParticles[firstTimeIndex];
  if (nFindings == 0) {
    edm::LogWarning("DatabasePDG") << "The particle required with name (" << name
                                   << ") was not found in the database!!";
    return nullptr;
  }
  if (nFindings >= 2) {
    edm::LogWarning("DatabasePDG") << "The particle required with name (" << name << ") was found with " << nFindings
                                   << " entries in the database. Check it out !!" << endl
                                   << "Returning the first instance found";
    return fParticles[firstTimeIndex];
  }
  return nullptr;
}
 
bool DatabasePDG::GetPDGParticleStatus(char *name) {
  int nFindings = 0;
  int firstTimeIndex = 0;
  for (int i = 0; i < fNParticles; i++) {
    if (!strcmp(name, fParticles[i]->GetName())) {
      if (nFindings == 0)
        firstTimeIndex = i;
      nFindings++;
    }
  }
  if (nFindings == 1)
    return fStatus[firstTimeIndex];
  if (nFindings == 0) {
    edm::LogWarning("DatabasePDG") << "The particle required with name (" << name
                                   << ") was not found in the database!!";
    return kFALSE;
  }
  if (nFindings >= 2) {
    edm::LogWarning("DatabasePDG") << "The particle status required for name (" << name << ") was found with "
                                   << nFindings << " entries in the database. Check it out !!" << endl
                                   << "Returning the first instance found";
    return fStatus[firstTimeIndex];
  }
  return kFALSE;
}
 
void DatabasePDG::DumpData(bool dumpAll) {
  cout << "***********************************************************************************************************"
       << endl;
  cout << "Dumping all the information contained in the database..." << endl;
  int nDecays = 0;
  int nGoodStatusDecays = 0;
  int nGoodStatusParticles = 0;
  for (int currPart = 0; currPart < fNParticles; currPart++) {
    nGoodStatusParticles += (fStatus[currPart] ? 1 : 0);
    nGoodStatusDecays += (fStatus[currPart] ? fParticles[currPart]->GetNDecayChannels() : 0);
    nDecays += fParticles[currPart]->GetNDecayChannels();
    if (!(dumpAll || (!dumpAll && fStatus[currPart])))
      continue;
    cout << "###### Particle: " << fParticles[currPart]->GetName() << " with PDG code "
         << fParticles[currPart]->GetPDG() << endl;
    cout << "   status          = " << fStatus[currPart] << endl;
    cout << "   mass            = " << fParticles[currPart]->GetMass() << " GeV" << endl;
    cout << "   width           = " << fParticles[currPart]->GetWidth() << " GeV" << endl;
    cout << "   2*spin          = " << int(2. * fParticles[currPart]->GetSpin()) << endl;
    cout << "   2*isospin       = " << int(2. * fParticles[currPart]->GetIsospin()) << endl;
    cout << "   2*isospin3      = " << int(2. * fParticles[currPart]->GetIsospinZ()) << endl;
    cout << "   u,d quarks      = " << int(fParticles[currPart]->GetLightQNumber()) << endl;
    cout << "   s quarks        = " << int(fParticles[currPart]->GetStrangeQNumber()) << endl;
    cout << "   c quarks        = " << int(fParticles[currPart]->GetCharmQNumber()) << endl;
    cout << "   anti u,d quarks = " << int(fParticles[currPart]->GetLightAQNumber()) << endl;
    cout << "   anti s quarks   = " << int(fParticles[currPart]->GetStrangeAQNumber()) << endl;
    cout << "   anti c quarks   = " << int(fParticles[currPart]->GetCharmAQNumber()) << endl;
    cout << "   baryon number   = " << int(fParticles[currPart]->GetBaryonNumber()) << endl;
    cout << "   strangeness     = " << int(fParticles[currPart]->GetStrangeness()) << endl;
    cout << "   charmness       = " << int(fParticles[currPart]->GetCharmness()) << endl;
    cout << "   electric charge = " << int(fParticles[currPart]->GetElectricCharge()) << endl;
    cout << "   full branching  = " << fParticles[currPart]->GetFullBranching() << endl;
    cout << "   decay modes     = " << fParticles[currPart]->GetNDecayChannels() << endl;
    for (int currChannel = 0; currChannel < fParticles[currPart]->GetNDecayChannels(); currChannel++) {
      cout << "   channel " << currChannel + 1 << " with branching "
           << fParticles[currPart]->GetDecayChannel(currChannel)->GetBranching() << endl;
      cout << "   daughters PDG codes: ";
      double daughtersMass = 0.0;
      for (int currDaughter = 0; currDaughter < fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters();
           currDaughter++) {
        cout << fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter) << "\t";
        ParticlePDG *daughter =
            GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
        daughtersMass += daughter->GetMass();
      }
      cout << endl;
      cout << "   daughters sum mass = " << daughtersMass << endl;
    }
  }
  if (dumpAll) {
    cout << "Finished dumping information for " << fNParticles << " particles with " << nDecays
         << " decay channels in total." << endl;
    cout << "**********************************************************************************************************"
            "***"
         << endl;
  } else {
    cout << "Finished dumping information for " << nGoodStatusParticles << "(" << fNParticles << ")"
         << " particles with " << nGoodStatusDecays << "(" << nDecays << ")"
         << " decay channels in total." << endl;
    cout << "**********************************************************************************************************"
            "***"
         << endl;
  }
}
 
int DatabasePDG::CheckImpossibleDecays(bool dump) {
  // Check the database for impossible decays
  int nImpossibleDecays = 0;
  for (int currPart = 0; currPart < fNParticles; currPart++) {
    if (!fStatus[currPart])
      continue;
    int allChannels = fParticles[currPart]->GetNDecayChannels();
    int allowedChannels = GetNAllowedChannels(fParticles[currPart], fParticles[currPart]->GetMass());
    if (dump) {
      cout << "Particle " << fParticles[currPart]->GetPDG() << " has " << allChannels
           << " decay channels specified in the database" << endl;
      cout << " Allowed channels assuming table mass = " << allowedChannels << endl;
    }
    if (dump && allChannels > 0 && allowedChannels == 0) {
      cout << "**********************************************************************" << endl;
      cout << "       All channels for this particles are not allowed" << endl;
      cout << "**********************************************************************" << endl;
    }
    if (dump && fParticles[currPart]->GetWidth() > 0. && allChannels == 0) {
      cout << "**********************************************************************" << endl;
      cout << "    Particle has finite width but no decay channels specified" << endl;
      cout << "**********************************************************************" << endl;
    }
    for (int currChannel = 0; currChannel < fParticles[currPart]->GetNDecayChannels(); currChannel++) {
      double motherMass = fParticles[currPart]->GetMass();
      double daughtersSumMass = 0.;
      for (int currDaughter = 0; currDaughter < fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters();
           currDaughter++) {
        ParticlePDG *daughter =
            GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
        daughtersSumMass += daughter->GetMass();
      }
      if (daughtersSumMass >= motherMass) {
        nImpossibleDecays++;
        if (dump) {
          cout << "Imposible decay for particle " << fParticles[currPart]->GetPDG() << endl;
          cout << "  Channel: " << fParticles[currPart]->GetPDG() << " --> ";
          for (int currDaughter = 0; currDaughter < fParticles[currPart]->GetDecayChannel(currChannel)->GetNDaughters();
               currDaughter++) {
            ParticlePDG *daughter =
                GetPDGParticle(fParticles[currPart]->GetDecayChannel(currChannel)->GetDaughterPDG(currDaughter));
            cout << daughter->GetPDG() << " ";
          }
          cout << endl;
          cout << "  Mother particle mass = " << motherMass << endl;
          cout << "  Daughters sum mass   = " << daughtersSumMass << endl;
        }
      }
    }
  }
  return nImpossibleDecays;
}
 
void DatabasePDG::SetUseCharmParticles(bool flag) {
  if (fNParticles > 0) {
    fUseCharmParticles = flag;
    for (int i = 0; i < fNParticles; i++) {
      if (fParticles[i]->GetCharmQNumber() > 0 || fParticles[i]->GetCharmAQNumber())
        fStatus[i] = flag;
    }
    SortParticles();
    return;
  } else
    fUseCharmParticles = flag;
  return;
}
 
void DatabasePDG::SetMinimumWidth(double value) {
  if (fNParticles > 0) {
    fMinimumWidth = value;
    for (int i = 0; i < fNParticles; i++) {
      if (fParticles[i]->GetWidth() < fMinimumWidth)
        fStatus[i] = kFALSE;
    }
    SortParticles();
    return;
  } else
    fMinimumWidth = value;
  return;
}
 
void DatabasePDG::SetMaximumWidth(double value) {
  if (fNParticles > 0) {
    fMaximumWidth = value;
    for (int i = 0; i < fNParticles; i++) {
      if (fParticles[i]->GetWidth() > fMaximumWidth)
        fStatus[i] = kFALSE;
    }
    SortParticles();
    return;
  } else
    fMaximumWidth = value;
  return;
}
 
void DatabasePDG::SetWidthRange(double min, double max) {
  if (fNParticles > 0) {
    fMinimumWidth = min;
    fMaximumWidth = max;
    for (int i = 0; i < fNParticles; i++) {
      if ((fParticles[i]->GetWidth() < fMinimumWidth) || (fParticles[i]->GetWidth() > fMaximumWidth))
        fStatus[i] = kFALSE;
    }
    SortParticles();
 
    return;
  } else {
    fMinimumWidth = min;
    fMaximumWidth = max;
  }
 
  return;
}
 
void DatabasePDG::SetMinimumMass(double value) {
  if (fNParticles > 0) {
    fMinimumMass = value;
    for (int i = 0; i < fNParticles; i++) {
      if (fParticles[i]->GetMass() < fMinimumMass)
        fStatus[i] = kFALSE;
    }
    SortParticles();
    return;
  } else
    fMinimumMass = value;
  return;
}
 
void DatabasePDG::SetMaximumMass(double value) {
  if (fNParticles > 0) {
    fMaximumMass = value;
    for (int i = 0; i < fNParticles; i++) {
      if (fParticles[i]->GetMass() > fMaximumMass)
        fStatus[i] = kFALSE;
    }
    SortParticles();
    return;
  } else
    fMaximumMass = value;
  return;
}
 
void DatabasePDG::SetMassRange(double min, double max) {
  if (fNParticles > 0) {
    fMinimumMass = min;
    fMaximumMass = max;
    for (int i = 0; i < fNParticles; i++) {
      if ((fParticles[i]->GetMass() < fMinimumMass) || (fParticles[i]->GetMass() > fMaximumMass))
        fStatus[i] = kFALSE;
    }
    SortParticles();
 
    return;
  } else {
    fMinimumMass = min;
    fMaximumMass = max;
  }
 
  return;
}
 
void DatabasePDG::SortParticles() {
  if (fNParticles < 2) {
    edm::LogWarning("DatabasePDG") << "No particles to sort. Load data first!!";
    return;
  }
 
  int nGoodStatus = 0;
  for (int i = 0; i < fNParticles; i++)
    if (fStatus[i])
      nGoodStatus++;
 
  if (nGoodStatus == fNParticles)  // if all particles have good status then there is nothing to do
    return;
 
  if (nGoodStatus == 0)  // no good status particles, again nothing to do
    return;
 
  int shifts = 1;
  while (shifts) {
    shifts = 0;
    for (int i = 0; i < fNParticles - 1; i++) {
      if (!fStatus[i] && fStatus[i + 1]) {  // switch if false status is imediately before a true status particle
        ParticlePDG *temporaryPointer = fParticles[i];
        fParticles[i] = fParticles[i + 1];
        fParticles[i + 1] = temporaryPointer;
        bool temporaryStatus = fStatus[i];
        fStatus[i] = fStatus[i + 1];
        fStatus[i + 1] = temporaryStatus;
        shifts++;
      }
    }
  }
 
  return;
}
 
int DatabasePDG::GetNParticles(bool all) {
  if (all)
    return fNParticles;
 
  int nGoodStatus = 0;
  for (int i = 0; i < fNParticles; i++)
    if (fStatus[i])
      nGoodStatus++;
  return nGoodStatus;
}
 
void DatabasePDG::UseThisListOfParticles(char *filename, bool exclusive) {
  if (fNParticles < 1) {
    edm::LogError("DatabasePDG") << "You must load the data before calling this function!!";
    return;
  }
 
  ifstream listFile;
  listFile.open(filename);
  if (!listFile) {
    edm::LogError("DatabasePDG") << "The ASCII file containing the PDG codes list (" << filename
                                 << ") was not found !!";
    return;
  }
 
  bool flaggedIndexes[kMaxParticles];
  for (int i = 0; i < kMaxParticles; i++)
    flaggedIndexes[i] = kFALSE;
  int pdg = 0;
  listFile.exceptions(ios::failbit);
  while (!listFile.eof()) {
    try {
      listFile >> pdg;
    } catch (ios::failure const &problem) {
      LogDebug("DatabasePDG") << "ios:failure in list file" << problem.what();
      break;
    }
    int found = 0;
    for (int i = 0; i < fNParticles; i++) {
      if (fParticles[i]->GetPDG() == pdg) {
        found++;
        flaggedIndexes[i] = kTRUE;
      }
    }
    if (!found) {
      edm::LogWarning("DatabasePDG") << "The particle with PDG code " << pdg
                                     << " was asked but not found in the database!!";
    }
    if (found > 1) {
      edm::LogWarning("DatabasePDG") << "The particle with PDG code " << pdg
                                     << " was found more than once in the database!!";
    }
  }
 
  if (exclusive) {
    for (int i = 0; i < kMaxParticles; i++)
      fStatus[i] = flaggedIndexes[i];
  } else {
    for (int i = 0; i < kMaxParticles; i++)
      fStatus[i] = (fStatus[i] && flaggedIndexes[i]);
  }
  SortParticles();
 
  return;
}
 
bool DatabasePDG::IsChannelAllowed(DecayChannel *channel, double motherMass) {
  double daughtersSumMass = 0.0;
  for (int i = 0; i < channel->GetNDaughters(); i++)
    daughtersSumMass += GetPDGParticle(channel->GetDaughterPDG(i))->GetMass();
  if (daughtersSumMass <= motherMass)
    return kTRUE;
  return kFALSE;
}
 
int DatabasePDG::GetNAllowedChannels(ParticlePDG *particle, double motherMass) {
  int nAllowedChannels = 0;
  for (int i = 0; i < particle->GetNDecayChannels(); i++)
    nAllowedChannels += (IsChannelAllowed(particle->GetDecayChannel(i), motherMass) ? 1 : 0);
  return nAllowedChannels;
}