decay_table.hpp

Go to the documentation of this file.

//   GAMBIT: Global and Modular BSM Inference Tool
//   *********************************************

#ifndef __decay_table_hpp__
#define __decay_table_hpp__

#include <map>
#include <set>
#include <string>
#include <sstream>

#include "gambit/Elements/slhaea_spec_helpers.hpp"
#include "gambit/Elements/mssm_slhahelp.hpp"
#include "gambit/Utils/util_types.hpp"
#include "gambit/Utils/standalone_error_handlers.hpp"
#include "gambit/Models/partmap.hpp"

namespace Gambit
{

  class DecayTable
  {
    public:

      // Forward declaration of DecayTable entry class, for readability.  Holds the info on all decays of a given particle.
      class Entry;

      DecayTable() {}
      DecayTable(str slha, int context = 0, bool force_SM_fermion_gauge_eigenstates = false);
      DecayTable(str slha, const std::map<int, int>& PDG_map, int context = 0, bool force_SM_fermion_gauge_eigenstates = false);
      DecayTable(const SLHAstruct&, int context = 0, bool force_SM_fermion_gauge_eigenstates = false);
      DecayTable(const SLHAstruct&, const std::map<int, int>&, int context = 0, bool force_SM_fermion_gauge_eigenstates = false);

      SLHAstruct getSLHAea(int SLHA_version, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;

      void writeSLHAfile(int SLHA_version, const str& filename, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;

      SLHAea::Block getSLHAea_block(int, std::pair<int,int>, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;
      SLHAea::Block getSLHAea_block(int, str, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;
      SLHAea::Block getSLHAea_block(int, str, int, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;

      Entry& operator()(std::pair<int,int>);
      Entry& operator()(str);
      Entry& operator()(str, int);
      const Entry& operator()(std::pair<int,int>) const;
      const Entry& operator()(str) const;
      const Entry& operator()(str, int) const;

      Entry& at(std::pair<int,int>);
      Entry& at(str);
      Entry& at(str, int);
      const Entry& at(std::pair<int,int>) const;
      const Entry& at(str) const;
      const Entry& at(str, int) const;

      std::map< std::pair<int,int>, Entry > particles;


      class Entry
      {
        private:

          void init(const SLHAea::Block&, int, bool force_SM_fermion_gauge_eigenstates = false);

          void check_particles_exist(std::multiset< std::pair<int,int> >&) const;

          void check_BF_validity(double, double, std::multiset< std::pair<int,int> >&) const;

          static void construct_key(std::multiset< std::pair<int,int> >&) {}
          template <typename... Args>
          static void construct_key(std::multiset< std::pair<int,int> >& key, str p1, Args... args)
          {
            construct_key(key, args...);
            key.insert(Models::ParticleDB().pdg_pair(p1));
          }
          template <typename... Args>
          static void construct_key(std::multiset< std::pair<int,int> >& key, str p1, int i1, Args... args)
          {
            construct_key(key, args...);
            key.insert(Models::ParticleDB().pdg_pair(p1, i1));
          }

        public:

          Entry() :  width_in_GeV(0.0), positive_error(0.0), negative_error(0.0), calculator(""), calculator_version(""), warnings(""), errors("") {}
          Entry(double width) :  width_in_GeV(width), calculator(""), calculator_version(""), warnings(""), errors("") {}
          Entry(const SLHAea::Block&, int context = 0, bool force_SM_fermion_gauge_eigenstates = false, str calculator = "", str calculator_version = "");
          Entry(const SLHAea::Block&, SLHAea::Block::const_iterator, int context = 0, bool force_SM_fermion_gauge_eigenstates = false,
           str calculator = "", str calculator_version = "");

          void set_BF(double, double, const std::vector<std::pair<int,int> >&);
          void set_BF(double, double, const std::vector<str>&);

          template <typename... Args>
          void set_BF(double BF, double error, std::pair<int,int> p1, Args... args)
          {
            std::pair<int,int> particles[] = {p1, args...};
            std::multiset< std::pair<int,int> > key(particles, particles+sizeof...(Args)+1);
            check_particles_exist(key);
            check_BF_validity(BF, error, key);
            channels[key] = std::pair<double, double>(BF, error);
          }

          template <typename... Args>
          void set_BF(double BF, double error, str p1, Args... args)
          {
            std::multiset< std::pair<int,int> > key;
            construct_key(key, p1, args...);
            check_BF_validity(BF, error, key);
            channels[key] = std::pair<double, double>(BF, error);
          }

          bool has_channel(const std::vector<std::pair<int,int> >&) const;
          bool has_channel(const std::vector<str>&) const;

          template <typename... Args>
          bool has_channel(std::pair<int,int> p1, Args... args) const
          {
            std::pair<int,int> particles[] = {p1, args...};
            std::multiset< std::pair<int,int> > key(particles, particles+sizeof...(Args)+1);
            check_particles_exist(key);
            return channels.find(key) != channels.end();
          }

          template <typename... Args>
          bool has_channel(str p1, Args... args) const
          {
            std::multiset< std::pair<int,int> > key;
            construct_key(key, p1, args...);
            return channels.find(key) != channels.end();
          }

          double BF(const std::vector<std::pair<int, int> >&) const;
          double BF(const std::vector<str>&) const;

          template <typename... Args>
          double BF(std::pair<int,int> p1, Args... args) const
          {
            std::pair<int,int> particles[] = {p1, args...};
            std::multiset< std::pair<int,int> > key(particles, particles+sizeof...(Args)+1);
            if (channels.find(key) == channels.end())
            {
              std::ostringstream err;
              err << "No branching fraction exists for the requested final states:";
              for (auto particle = key.begin(); particle != key.end(); ++particle)
              {
                err << " {" << particle->first << ", " << particle->second << "}";
              }
              model_error().raise(LOCAL_INFO,err.str());
            }
            return channels.at(key).first;
          }

          template <typename... Args>
          double BF(str p1, Args... args) const
          {
            std::multiset< std::pair<int,int> > key;
            construct_key(key, p1, args...);
            if (channels.find(key) == channels.end())
            {
              std::ostringstream err;
              err << "No branching fraction exists for the requested final states:";
              for (auto particle = key.begin(); particle != key.end(); ++particle)
              {
                err << " {" << particle->first << ", " << particle->second << "}";
              }
              model_error().raise(LOCAL_INFO,err.str());
            }
            return channels.at(key).first;
          }

          template <typename... Args>
          double BF_error(std::pair<int,int> p1, Args... args) const
          {
            std::pair<int,int> particles[] = {p1, args...};
            std::multiset< std::pair<int,int> > key(particles, particles+sizeof...(Args)+1);
            if (channels.find(key) == channels.end())
            {
              std::ostringstream err;
              err << "No branching fraction exists for the requested final states:";
              for (auto particle = key.begin(); particle != key.end(); ++particle)
              {
                err << " {" << particle->first << ", " << particle->second << "}";
              }
              model_error().raise(LOCAL_INFO,err.str());
            }
            return channels.at(key).second;
          }

          template <typename... Args>
          double BF_error(str p1, Args... args) const
          {
            std::multiset< std::pair<int,int> > key;
            construct_key(key, p1, args...);
            if (channels.find(key) == channels.end())
            {
              std::ostringstream err;
              err << "No branching fraction exists for the requested final states:";
              for (auto particle = key.begin(); particle != key.end(); ++particle)
              {
                err << " {" << particle->first << ", " << particle->second << "}";
              }
              model_error().raise(LOCAL_INFO,err.str());
            }
            return channels.at(key).second;
          }

          template <typename... Args>
          std::pair<double, double> BF_with_error(std::pair<int,int> p1, Args... args) const
          {
            std::pair<int,int> particles[] = {p1, args...};
            std::multiset< std::pair<int,int> > key(particles, particles+sizeof...(Args)+1);
            if (channels.find(key) == channels.end())
            {
              std::ostringstream err;
              err << "No branching fraction exists for the requested final states:";
              for (auto particle = key.begin(); particle != key.end(); ++particle)
              {
                err << " {" << particle->first << ", " << particle->second << "}";
              }
              model_error().raise(LOCAL_INFO,err.str());
            }
            return channels.at(key);
          }

          template <typename... Args>
          std::pair<double, double> BF_with_error(str p1, Args... args) const
          {
            std::multiset< std::pair<int,int> > key;
            construct_key(key, p1, args...);
            if (channels.find(key) == channels.end())
            {
              std::ostringstream err;
              err << "No branching fraction exists for the requested final states:";
              for (auto particle = key.begin(); particle != key.end(); ++particle)
              {
                err << " {" << particle->first << ", " << particle->second << "}";
              }
              model_error().raise(LOCAL_INFO,err.str());
            }
            return channels.at(key);
          }

          SLHAea::Block getSLHAea_block(int, str, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;
          SLHAea::Block getSLHAea_block(int, str, int, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;
          SLHAea::Block getSLHAea_block(int, std::pair<int,int>, bool include_zero_bfs=false, const mass_es_pseudonyms& psn=mass_es_pseudonyms()) const;

          double sum_BF() const;

          double width_in_GeV;

          double positive_error;

          double negative_error;

          str calculator;
          str calculator_version;
          str warnings;
          str errors;

          std::map< std::multiset< std::pair<int,int> >, std::pair<double, double> > channels;

      };

  };


}


#endif //#defined __decay_table_hpp__