gambit is hosted by Hepforge, IPPP Durham
 GAMBIT  v1.5.0-2191-ga4742ac a Global And Modular Bsm Inference Tool
DarkBit_utils.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
19
23
24 //#define DARKBIT_DEBUG
25
26 namespace Gambit
27 {
28  namespace DarkBit
29  {
30  namespace DarkBit_utils
31  {
40  template <int i>
41  double gamma3bdy_limits(double Eg, double M_DM, double m1, double m2)
42  {
43  double x = Eg/M_DM;
44  double x0, x1;
45  // Check if kinematic constraints are satisfied
46  if((1.-pow(m1+m2,2)/(4*M_DM*M_DM))<=x)
47  {
48  x0 = x1 = 0;
49  }
50  else
51  {
52  double eta = pow(m1/M_DM,2);
53  double diffeta=pow(m2/M_DM,2);
54  diffeta = 0.25*(eta-diffeta);
55  double f1 = 0.25*eta + diffeta*x/(2*(1-x));
56  double f2 = sqrt(pow(1+diffeta/(1-x),2)-eta/(1-x));
57  double aint = f1 + 0.5*(1-f2)*x;
58  double bint = f1 + 0.5*(1+f2)*x;
59  // Now convert these limits to limits on E1
60  double f3 = pow(0.5*m2/M_DM,2);
61  x0 = M_DM*(1-x+aint-f3);
62  x1 = M_DM*(1-x+bint-f3);
63  }
64  if ( i == 0 ) return x0;
65  else return x1;
66  }
67  template double gamma3bdy_limits<0>(double, double, double, double);
68  template double gamma3bdy_limits<1>(double, double, double, double);
69
70
71  // Convert between mass and flavour eigenstate identifiers
72  std::string str_flav_to_mass(std::string flav)
73  {
74  if (flav=="e+") return "e+_1";
75  if (flav=="mu+") return "e+_2";
76  if (flav=="tau+") return "e+_3";
77  if (flav=="u") return "u_1";
78  if (flav=="d") return "d_1";
79  if (flav=="c") return "u_2";
80  if (flav=="s") return "d_2";
81  if (flav=="t") return "u_3";
82  if (flav=="b") return "d_3";
83  if (flav=="e-") return "e-_1";
84  if (flav=="mu-") return "e-_2";
85  if (flav=="tau-") return "e-_3";
86  if (flav=="ubar") return "ubar_1";
87  if (flav=="dbar") return "dbar_1";
88  if (flav=="cbar") return "ubar_2";
89  if (flav=="sbar") return "dbar_2";
90  if (flav=="tbar") return "ubar_3";
91  if (flav=="bbar") return "dbar_3";
92  return flav;
93  }
94
95  std::string str_mass_to_flav(std::string mass)
96  {
97  if (mass=="e+_1" ) return "e+";
98  if (mass=="e+_2" ) return "mu+";
99  if (mass=="e+_3" ) return "tau+";
100  if (mass=="u_1" ) return "u";
101  if (mass=="d_1" ) return "d";
102  if (mass=="u_2" ) return "c";
103  if (mass=="d_2" ) return "s";
104  if (mass=="u_3" ) return "t";
105  if (mass=="d_3" ) return "b";
106  if (mass=="e-_1" ) return "e-";
107  if (mass=="e-_2" ) return "mu-";
108  if (mass=="e-_3" ) return "tau-";
109  if (mass=="ubar_1" ) return "ubar";
110  if (mass=="dbar_1" ) return "dbar";
111  if (mass=="ubar_2" ) return "cbar";
112  if (mass=="dbar_2" ) return "sbar";
113  if (mass=="ubar_3" ) return "tbar";
114  if (mass=="dbar_3" ) return "bbar";
115  return mass;
116  }
117
118
119  void ImportDecays(std::string pID, TH_ProcessCatalog &catalog,
120  std::set<std::string> &importedDecays,
121  const DecayTable* tbl, double minBranching,
122  std::vector<std::string> excludeDecays)
123  {
124  if(importedDecays.count(pID) == 1) return;
125 #ifdef DARKBIT_DEBUG
126  std::cout << "Importing decay information for: " << pID << std::endl;
127 #endif
128  importedDecays.insert(pID);
129  const double m_init = catalog.getParticleProperty(pID).mass;
130  const DecayTable::Entry* entry;
131  try{entry = &(tbl->at(pID));}
132  catch(const std::out_of_range& oor)
133  {
134 #ifdef DARKBIT_DEBUG
135  std::cout << " ! No entry found in DecayTable object (nothing imported) !" << std::endl;
136 #endif
137  return;
138  }
139  double totalWidth = entry->width_in_GeV;
140 #ifdef DARKBIT_DEBUG
141  std::cout << " totalWidth = " << totalWidth << std::endl;
142 #endif
143  if(totalWidth>0)
144  {
145  TH_Process process(pID);
146  process.genRateMisc = daFunk::cnst(0.);
147  for(auto fState_it = entry->channels.begin();
148  fState_it!= entry->channels.end(); ++fState_it)
149  {
150  double bFraction = (fState_it->second).first;
151  if(bFraction>minBranching)
152  {
153  std::vector<std::string> pIDs;
154  double m_final = 0;
155 #ifdef DARKBIT_DEBUG
156  std::cout << " Final state: ";
157 #endif
158  for(auto pit = fState_it->first.begin();
159  pit != fState_it->first.end(); ++pit)
160  {
161  std::string name = Models::ParticleDB().long_name(*pit);
162  name = DarkBit_utils::str_flav_to_mass(name);
163  m_final += catalog.getParticleProperty(name).mass;
164  pIDs.push_back(name);
165 #ifdef DARKBIT_DEBUG
166  std::cout << name << " ";
167 #endif
168  }
169  double partialWidth = totalWidth * bFraction;
170 #ifdef DARKBIT_DEBUG
171  std::cout << std::endl;
172  std::cout << " m_init = " << m_init << std::endl;
173  std::cout << " m_final = " << m_final << std::endl;
174  std::cout << " bFraction = " << bFraction << std::endl;
175 #endif
176  if(m_final<=m_init)
177  {
178  process.channelList.push_back(
179  TH_Channel(pIDs, daFunk::cnst(partialWidth)));
180  // Recursively import decays of final states (for cascades)
181  for(auto f_it=pIDs.begin(); f_it!=pIDs.end();++f_it)
182  {
183  if (std::find(excludeDecays.begin(), excludeDecays.end(), *f_it) == excludeDecays.end())
184  ImportDecays(*f_it, catalog, importedDecays, tbl, minBranching, excludeDecays);
185  }
186  }
187 #ifdef DARKBIT_DEBUG
188  else
189  std::cout << " ! Process ignored since it is offshell !" << std::endl;
190 #endif
191  }
192  }
193  catalog.processList.push_back(process);
194  }
195  }
196
197  } // namespace DarkBit_utils
198  } // namespace DarkBit
199 } // namespace Gambit
200
201 #undef DARKBIT_DEBUG
std::vector< TH_Channel > channelList
List of channels.
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
DecayTable entry class. Holds the info on all decays of a given particle.
Definition: decay_table.hpp:96
TH_ParticleProperty getParticleProperty(str) const
Retrieve properties of a given particle involved in one or more processes in this catalog...
template double gamma3bdy_limits< 0 >(double, double, double, double)
Entry & at(std::pair< int, int >)
Get entry in decay table for a give particle, throwing an error if particle is absent.
void ImportDecays(std::string pID, TH_ProcessCatalog &catalog, std::set< std::string > &importedDecays, const DecayTable *tbl, double minBranching, std::vector< std::string > excludeDecays=std::vector< std::string >())
str long_name(str, int) const
Retrieve the long name, from the short name and index.
Definition: partmap.cpp:124
double mass
Particle mass (GeV)
double width_in_GeV
Total particle width (in GeV)
partmap & ParticleDB()
Database accessor function.
Definition: partmap.cpp:36
A container holding all annihilation and decay initial states relevant for DarkBit.
template double gamma3bdy_limits< 1 >(double, double, double, double)
Header file that includes all GAMBIT headers required for a module source file.
All data on a single annihilation or decay channel, e.g.
daFunk::Funk genRateMisc
double gamma3bdy_limits(double Eg, double M_DM, double m1, double m2)
Calculate kinematical limits for three-body final states.
GAMBIT native decay table class.
Definition: decay_table.hpp:35
Funk cnst(double x, Args... argss)
Definition: daFunk.hpp:606
double pow(const double &a)
Outputs a^i.
std::string str_flav_to_mass(std::string flav)
A container for a single process.
std::map< std::multiset< std::pair< int, int > >, std::pair< double, double > > channels
The actual underlying map of channels to their BFs.
std::vector< TH_Process > processList
Vector of all processes in this catalog.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Utility functions for DarkBit.
std::string str_mass_to_flav(std::string mass)