smash::DecayModes Class Reference

#include <decaymodes.h>

The DecayModes class is used to store and update information about decay branches (i.e.

the possible children a mother particle can decay into), including their relative weights.

If you want to find a DecayModes object for a specific particle type use ParticleType::decay_modes().

Definition at line 29 of file decaymodes.h.

Collaboration diagram for smash::DecayModes:

Classes
struct	InvalidDecay
struct	LoadFailure
struct	MissingDecays
struct	ReferencedParticleNotFound

Public Member Functions
void	add_mode (ParticleTypePtr mother, double ratio, int L, ParticleTypePtrList particle_types)
	Add a decay mode using all necessary information. More...
void	add_mode (DecayBranchPtr branch)
	Add a decay mode from an already existing decay branch. More...
bool	renormalize (const std::string &name)
	Renormalize the branching ratios to add up to 1. More...
bool	is_empty () const
const DecayBranchList &	decay_mode_list () const

Static Public Member Functions
static void	load_decaymodes (const std::string &input)
	Loads the DecayModes map as described in the input string. More...
static DecayType *	get_decay_type (ParticleTypePtr mother, ParticleTypePtrList particle_types, int L)
	Retrieve a decay type. More...

Private Attributes
DecayBranchList	decay_modes_
	Vector of decay modes. More...

Static Private Attributes
static std::vector< DecayModes > *	all_decay_modes = nullptr
	A list of all DecayModes objects using the same indexing as all_particle_types. More...

Friends
const DecayModes &	ParticleType::decay_modes () const
	allow ParticleType::decay_modes to access all_decay_modes More...

Member Function Documentation

add_mode() [1/2]

void smash::DecayModes::add_mode	(	ParticleTypePtr	mother,
		double	ratio,
		int	L,
		ParticleTypePtrList	particle_types
	)

Add a decay mode using all necessary information.

Parameters

[in]	mother	the particle which decays
[in]	ratio	the weight to add to the current mode
[in]	L	angular momentum
[in]	particle_types	a list of the products of the decay

Definition at line 29 of file decaymodes.cc.

                                                              {
  DecayType *type = get_decay_type(mother, particle_types, L);
  // Check if mode already exists: if yes, add weight.
  for (auto &mode : decay_modes_) {
    if (type == &mode->type()) {
      mode->set_weight(mode->weight() + ratio);
      return;
    }
  }
  // Add new mode.
  decay_modes_.push_back(make_unique<DecayBranch>(*type, ratio));
}

Here is the call graph for this function:

Here is the caller graph for this function:

add_mode() [2/2]

void smash::DecayModes::add_mode ( DecayBranchPtr  branch )

inline

Add a decay mode from an already existing decay branch.

Parameters

[in]

branch

the decay branch to add

Definition at line 47 of file decaymodes.h.

                                       {
    decay_modes_.push_back(std::move(branch));
  }

renormalize()

bool smash::DecayModes::renormalize

(

const std::string &

name

)

Renormalize the branching ratios to add up to 1.

Parameters

[in]

name

the name of the decaying particle

Returns

if the branching ratios were renormalized by more than 1%

Definition at line 98 of file decaymodes.cc.

                                                  {
  const auto &log = logger<LogArea::DecayModes>();
  double sum = 0.;
  bool is_large_renormalization = false;
  for (auto &mode : decay_modes_) {
    sum += mode->weight();
  }
  if (std::abs(sum - 1.) < really_small) {
    log.debug("Particle ", name,
              ": Extremely small renormalization constant: ", sum,
              "\n=> Skipping the renormalization.");
  } else {
    is_large_renormalization = (std::abs(sum - 1.) > 0.01);
    log.debug("Particle ", name, ": Renormalizing decay modes with ", sum);
    double new_sum = 0.0;
    for (auto &mode : decay_modes_) {
      mode->set_weight(mode->weight() / sum);
      new_sum += mode->weight();
    }
    log.debug("After renormalization sum of ratios is ", new_sum);
  }
  return is_large_renormalization;
}

is_empty()

bool smash::DecayModes::is_empty ( ) const

inline

Returns

true if empty (i.e. no decay modes)

Definition at line 60 of file decaymodes.h.

{ return decay_modes_.empty(); }

Here is the caller graph for this function:

decay_mode_list()

const DecayBranchList& smash::DecayModes::decay_mode_list ( ) const

inline

Returns

pass out the decay modes list

Definition at line 63 of file decaymodes.h.

{ return decay_modes_; }

Here is the caller graph for this function:

load_decaymodes()

void smash::DecayModes::load_decaymodes ( const std::string &  input )

static

Loads the DecayModes map as described in the input string.

It does sanity checking - that the particles it talks about are in the ParticleType map.

Parameters

[in]

input

the full decaymodes input file, as a string

Exceptions

MissingDecays	if there are no decays specified for an unstable particle
LoadFailure	if there are duplicate entries detailing decaymodes for the same particle, or if the angular momentum is smaller than 0 or larger than 4
InvalidDecay	if product of the decay does not exist in the particle list or as an isospin multiplet, or if decay is forbidden by isospin or charge conservation, or if sum of the minimum mass of products is less than the pole mass of the mother particle
runtime_error	if there are less than 2 or more than 3 products to a given decay branch, or if a branch cannot be found to be part of an isospin multiplet

Definition at line 163 of file decaymodes.cc.

                                                       {
  const auto &log = logger<LogArea::DecayModes>();
  // create the DecayType vector first, then it outlives the DecayModes vector,
  // which references the DecayType objects.
  static std::vector<DecayTypePtr> decaytypes;
  decaytypes.clear();  // in case an exception was thrown and should try again
  // ten decay types per decay mode should be a good guess.
  decaytypes.reserve(10 * ParticleType::list_all().size());
  all_decay_types = &decaytypes;

  static std::vector<DecayModes> decaymodes;
  decaymodes.clear();  // in case an exception was thrown and should try again
  decaymodes.resize(ParticleType::list_all().size());
  all_decay_modes = &decaymodes;

  const IsoParticleType *isotype_mother = nullptr;
  ParticleTypePtrList mother_states;
  std::vector<DecayModes> decay_modes_to_add;  // one for each mother state
  int total_large_renormalized = 0;

  const auto end_of_decaymodes = [&]() {
    if (isotype_mother == nullptr) {  // at the start of the file
      return;
    }
    // Loop over all states in the mother multiplet and add modes
    for (size_t m = 0; m < mother_states.size(); m++) {
      if (decay_modes_to_add[m].is_empty() && !mother_states[m]->is_stable()) {
        throw MissingDecays("No decay modes found for particle " +
                            mother_states[m]->name());
      }
      bool is_large_renorm =
          decay_modes_to_add[m].renormalize(mother_states[m]->name());
      total_large_renormalized += is_large_renorm;
      PdgCode pdgcode = mother_states[m]->pdgcode();
      /* Add the list of decay modes for this particle type */
      decaymodes[find_offset(pdgcode)] = std::move(decay_modes_to_add[m]);
    }
    if (isotype_mother->has_anti_multiplet()) {
      /* Construct the decay modes for the anti-multiplet.  */
      log.debug("generating decay modes for anti-multiplet: " +
                isotype_mother->name());
      for (const auto &state : mother_states) {
        PdgCode pdg = state->pdgcode();
        PdgCode pdg_anti = pdg.get_antiparticle();
        const ParticleType &type_anti = ParticleType::find(pdg_anti);
        DecayModes &decay_modes_orig = decaymodes[find_offset(pdg)];
        DecayModes &decay_modes_anti = decaymodes[find_offset(pdg_anti)];
        for (const auto &mode : decay_modes_orig.decay_mode_list()) {
          ParticleTypePtrList list = mode->particle_types();
          for (auto &type : list) {
            if (type->has_antiparticle()) {
              type = type->get_antiparticle();
            }
          }
          decay_modes_anti.add_mode(&type_anti, mode->weight(),
                                    mode->angular_momentum(), list);
        }
      }
    }
  };

  // Track the line number for better error messages.
  // FIXME: At the moment this does not include comments and empty lines.
  uint64_t linenumber = 1;
  for (const Line &line : line_parser(input)) {
    const auto trimmed = trim(line.text);
    assert(!trimmed.empty());  // trim(line.text) is never empty,
                               // else line_parser is broken
    if (trimmed.find_first_of(" \t") == std::string::npos) {
      // a single record on one line signifies a new decay mode section
      end_of_decaymodes();
      std::string name = trim(line.text);
      isotype_mother = &IsoParticleType::find(name);
      mother_states = isotype_mother->get_states();
      decay_modes_to_add.clear();
      decay_modes_to_add.resize(mother_states.size());
      log.debug("reading decay modes for " + name);
      // check if any of the states have decay modes already
      for (size_t m = 0; m < mother_states.size(); m++) {
        PdgCode pdgcode = mother_states[m]->pdgcode();
        if (!decaymodes[find_offset(pdgcode)].is_empty()) {
          throw LoadFailure("Duplicate entry for " + name +
                            " in decaymodes.txt:" + std::to_string(linenumber));
        }
      }
    } else {
      std::istringstream lineinput(line.text);
      std::vector<std::string> decay_particles;
      decay_particles.reserve(3);
      double ratio;
      lineinput >> ratio;

      int L;
      lineinput >> L;
      if (L < 0) {
        throw LoadFailure("Invalid angular momentum '" + std::to_string(L) +
                          "' in decaymodes.txt:" + std::to_string(line.number) +
                          ": '" + line.text + "'");
      }

      std::string name;
      lineinput >> name;
      bool multi = true;
      while (lineinput) {
        decay_particles.emplace_back(name);
        const auto isotype = IsoParticleType::try_find(name);
        const bool is_multiplet = isotype;
        const bool is_state = ParticleType::exists(name);
        if (!is_multiplet && !is_state) {
          throw InvalidDecay(
              "Daughter " + name +
              " is neither an isospin multiplet nor a particle." + " (line " +
              std::to_string(linenumber) + ": \"" + trimmed + "\")");
        }
        const bool is_hadronic_multiplet =
            is_multiplet && isotype->get_states()[0]->is_hadron();
        multi &= is_hadronic_multiplet;
        lineinput >> name;
      }
      Parity parity;
      bool is_strong_decay;
      const int s0 = isotype_mother->spin();
      int min_L = 0;
      int max_L = 0;
      if (multi) {
        /* References to isospin multiplets: Automatically determine all valid
         * combinations and calculate Clebsch-Gordan factors */
        switch (decay_particles.size()) {
          case 2: {
            const IsoParticleType &isotype_daughter_1 =
                IsoParticleType::find(decay_particles[0]);
            const IsoParticleType &isotype_daughter_2 =
                IsoParticleType::find(decay_particles[1]);
            parity = isotype_daughter_1.parity() * isotype_daughter_2.parity();
            is_strong_decay = isotype_daughter_1.is_hadron() &&
                              isotype_daughter_2.is_hadron();
            const int s1 = isotype_daughter_1.spin();
            const int s2 = isotype_daughter_2.spin();
            min_L = min_angular_momentum(s0, s1, s2);
            max_L = (s0 + s1 + s2) / 2;
            // loop through multiplets
            bool forbidden_by_isospin = true;
            for (size_t m = 0; m < mother_states.size(); m++) {
              for (const auto &daughter1 : isotype_daughter_1.get_states()) {
                for (const auto &daughter2 : isotype_daughter_2.get_states()) {
                  // calculate Clebsch-Gordan factor
                  const double cg_sqr = isospin_clebsch_gordan_sqr_2to1(
                      *daughter1, *daughter2, *mother_states[m]);
                  if (cg_sqr > 0.) {
                    // add mode
                    log.debug(
                        "decay mode generated: " + mother_states[m]->name() +
                        " -> " + daughter1->name() + " " + daughter2->name() +
                        " (" + std::to_string(ratio * cg_sqr) + ")");
                    decay_modes_to_add[m].add_mode(mother_states[m],
                                                   ratio * cg_sqr, L,
                                                   {daughter1, daughter2});
                    forbidden_by_isospin = false;
                  }
                }
              }
            }
            if (forbidden_by_isospin) {
              std::stringstream s;
              s << ",\nwhere isospin mother: " << isotype_mother->isospin()
                << ", daughters: " << isotype_daughter_1.isospin() << " "
                << isotype_daughter_2.isospin();
              throw InvalidDecay(isotype_mother->name() +
                                 " decay mode is forbidden by isospin: \"" +
                                 line.text + "\"" + s.str());
            }
            break;
          }
          case 3: {
            const IsoParticleType &isotype_daughter_1 =
                IsoParticleType::find(decay_particles[0]);
            const IsoParticleType &isotype_daughter_2 =
                IsoParticleType::find(decay_particles[1]);
            const IsoParticleType &isotype_daughter_3 =
                IsoParticleType::find(decay_particles[2]);
            parity = isotype_daughter_1.parity() * isotype_daughter_2.parity() *
                     isotype_daughter_3.parity();
            is_strong_decay = isotype_daughter_1.is_hadron() &&
                              isotype_daughter_2.is_hadron() &&
                              isotype_daughter_3.is_hadron();
            const int s1 = isotype_daughter_1.spin();
            const int s2 = isotype_daughter_2.spin();
            const int s3 = isotype_daughter_2.spin();
            min_L = min_angular_momentum(s0, s1, s2, s3);
            max_L = (s0 + s1 + s2 + s3) / 2;
            // loop through multiplets
            for (size_t m = 0; m < mother_states.size(); m++) {
              for (const auto &daughter1 : isotype_daughter_1.get_states()) {
                for (const auto &daughter2 : isotype_daughter_2.get_states()) {
                  for (const auto &daughter3 :
                       isotype_daughter_3.get_states()) {
                    const double cg_sqr = isospin_clebsch_gordan_sqr_3to1(
                        *daughter1, *daughter2, *daughter3, *mother_states[m]);
                    if (cg_sqr > 0.) {
                      // add mode
                      log.debug(
                          "decay mode generated: " + mother_states[m]->name() +
                          " -> " + daughter1->name() + " " + daughter2->name() +
                          " " + daughter3->name() + " (" +
                          std::to_string(ratio * cg_sqr) + ")");
                      decay_modes_to_add[m].add_mode(
                          mother_states[m], ratio * cg_sqr, L,
                          {daughter1, daughter2, daughter3});
                    }
                  }
                }
              }
            }
            break;
          }
          default:
            throw std::runtime_error(
                "References to isospin multiplets only "
                "allowed in two-body or three-body decays: " +
                line.text + " (line " + std::to_string(linenumber) + ": \"" +
                trimmed + "\")");
        }
      } else {
        /* References to specific states, not multiplets:
         * Loop over all mother states and check charge conservation. */
        ParticleTypePtrList types;
        int charge = 0;
        parity = Parity::Pos;
        is_strong_decay = true;
        for (auto part : decay_particles) {
          try {
            types.push_back(IsoParticleType::find_state(part));
          } catch (std::runtime_error &e) {
            throw std::runtime_error(std::string() + e.what() + " (line " +
                                     std::to_string(linenumber) + ": \"" +
                                     trimmed + "\")");
          }
          charge += types.back()->charge();
          parity *= types.back()->parity();
          is_strong_decay &= types.back()->is_hadron();
        }
        if (types.size() == 2) {
          const int s1 = types[0]->spin();
          const int s2 = types[1]->spin();
          min_L = min_angular_momentum(s0, s1, s2);
          max_L = (s0 + s1 + s2) / 2;
        } else if (types.size() == 3) {
          const int s1 = types[0]->spin();
          const int s2 = types[1]->spin();
          const int s3 = types[2]->spin();
          min_L = min_angular_momentum(s0, s1, s2, s3);
          max_L = (s0 + s1 + s2 + s3) / 2;
        } else {
          throw InvalidDecay(isotype_mother->name() +
                             " decay mode has an invalid number of particles"
                             " in the final state " +
                             "(line " + std::to_string(linenumber) + ": \"" +
                             trimmed + "\")");
        }
        bool no_decays = true;
        for (size_t m = 0; m < mother_states.size(); m++) {
          if (mother_states[m]->charge() == charge) {
            log.debug("decay mode found: " + mother_states[m]->name() + " -> " +
                      std::to_string(decay_particles.size()));
            decay_modes_to_add[m].add_mode(mother_states[m], ratio, L, types);
            no_decays = false;
          }
        }
        if (no_decays) {
          throw InvalidDecay(isotype_mother->name() +
                             " decay mode violates charge conservation " +
                             "(line " + std::to_string(linenumber) + ": \"" +
                             trimmed + "\")");
        }
      }
      // Take angular momentum into account.
      if (L % 2 == 1) {
        parity = -parity;
      }
      // Make sure the decay has the correct parity.
      if (is_strong_decay && parity != mother_states[0]->parity()) {
        throw InvalidDecay(mother_states[0]->name() +
                           " decay mode violates parity conservation " +
                           "(line " + std::to_string(linenumber) + ": \"" +
                           trimmed + "\")");
      }
      // Make sure the decay has a correct angular momentum.
      if (L < min_L || L > max_L) {
        throw InvalidDecay(
            mother_states[0]->name() +
            " decay mode violates angular momentum conservation: " +
            std::to_string(L) + " not in [" + std::to_string(min_L) + ", " +
            std::to_string(max_L) + "] (line " + std::to_string(linenumber) +
            ": \"" + trimmed + "\")");
      }
    }
    linenumber++;
  }
  end_of_decaymodes();

  // Check whether the mother's pole mass is strictly larger than the minimal
  // masses of the daughters. This is required by the Manley-Saleski ansatz.
  const auto &particles = ParticleType::list_all();
  for (const auto &mother : particles) {
    if (mother.is_stable()) {
      continue;
    }
    const auto &decays = mother.decay_modes().decay_mode_list();
    for (const auto &decay : decays) {
      if (mother.mass() <= decay->threshold()) {
        std::stringstream s;
        s << mother.name() << " →  ";
        for (const auto p : decay->particle_types()) {
          s << p->name();
        }
        s << " with " << mother.mass() << " ≤ " << decay->threshold();
        throw InvalidDecay(
            "For all decays, the minimum mass of daughters"
            "must be smaller\nthan the mother's pole mass "
            "(Manley-Saleski Ansatz)\n"
            "Violated by the following decay: " +
            s.str());
      }
    }
  }
  if (total_large_renormalized > 0) {
    log.warn("Branching ratios of ", total_large_renormalized,
             " hadrons were renormalized by more than 1% to have sum 1.");
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

get_decay_type()

DecayType * smash::DecayModes::get_decay_type ( ParticleTypePtr  mother, ParticleTypePtrList  particle_types, int  L  )

static

Retrieve a decay type.

Parameters

[in]	mother	the decaying particle
[in]	particle_types	the products of the decay
[in]	L	the angular momentum

Returns

the corresponding DecayType object

Exceptions

InvalidDecay

if there are less than 2 or more than 3 products

Definition at line 43 of file decaymodes.cc.

                                             {
  assert(all_decay_types != nullptr);

  // check if the decay type already exisits
  for (const auto &type : *all_decay_types) {
    if (type->has_mother(mother) && type->has_particles(particle_types) &&
        type->angular_momentum() == L) {
      return type.get();
    }
  }

  // if the type does not exist yet, create a new one
  switch (particle_types.size()) {
    case 2:
      if (is_dilepton(particle_types[0]->pdgcode(),
                      particle_types[1]->pdgcode())) {
        all_decay_types->emplace_back(
            make_unique<TwoBodyDecayDilepton>(particle_types, L));
      } else if (particle_types[0]->is_stable() &&
                 particle_types[1]->is_stable()) {
        all_decay_types->emplace_back(
            make_unique<TwoBodyDecayStable>(particle_types, L));
      } else if (particle_types[0]->is_stable() ||
                 particle_types[1]->is_stable()) {
        all_decay_types->emplace_back(
            make_unique<TwoBodyDecaySemistable>(particle_types, L));
      } else {
        all_decay_types->emplace_back(
            make_unique<TwoBodyDecayUnstable>(particle_types, L));
      }
      break;
    case 3:
      if (has_lepton_pair(particle_types[0]->pdgcode(),
                          particle_types[1]->pdgcode(),
                          particle_types[2]->pdgcode())) {
        all_decay_types->emplace_back(
            make_unique<ThreeBodyDecayDilepton>(mother, particle_types, L));
      } else {
        all_decay_types->emplace_back(
            make_unique<ThreeBodyDecay>(particle_types, L));
      }
      break;
    default:
      throw InvalidDecay(
          "DecayModes::get_decay_type was instructed to add a decay mode "
          "with " +
          std::to_string(particle_types.size()) +
          " particles. This is an invalid input.");
  }

  return all_decay_types->back().get();
}

Here is the call graph for this function:

Here is the caller graph for this function:

Friends And Related Function Documentation

ParticleType::decay_modes

const DecayModes& ParticleType::decay_modes ( ) const

friend

allow ParticleType::decay_modes to access all_decay_modes

Member Data Documentation

decay_modes_

DecayBranchList smash::DecayModes::decay_modes_

private

Vector of decay modes.

Each mode consists of a vector of the pdg codes of decay products and a ratio of this decay mode compared to all possible modes

Definition at line 123 of file decaymodes.h.

all_decay_modes

std::vector< DecayModes > * smash::DecayModes::all_decay_modes = nullptr

staticprivate

A list of all DecayModes objects using the same indexing as all_particle_types.

Global pointer to the decay modes list.

Definition at line 132 of file decaymodes.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/decaymodes.h
• src/decaymodes.cc