gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
ProcessCatalog.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
17 
21 
22 //#define DARKBIT_DEBUG
23 
24 namespace Gambit {
25  namespace DarkBit {
26 
27  // TH_ParticleProperty definitions
28 
29  TH_ParticleProperty::TH_ParticleProperty(double mass, unsigned int spin2)
30  : mass(mass), spin2(spin2)
31  {}
32 
33 
34  // TH_Channel definitions
35 
37  TH_Channel::TH_Channel(std::vector<str> finalStateIDs, daFunk::Funk genRate)
38  : finalStateIDs(finalStateIDs), nFinalStates(finalStateIDs.size()),
39  genRate(genRate)
40  {
41  if ( nFinalStates < 2 )
42  {
43  DarkBit_error().raise(LOCAL_INFO, "Need at least two final state particles.");
44  }
45  }
46 
49  {
50  logger() << "Channel: ";
51  for ( auto it = finalStateIDs.begin(); it != finalStateIDs.end(); it++ )
52  {
53  logger() << *it << " ";
54  }
55  logger() << EOM;
56  }
57 
62  {
63  return std::find(finalStateIDs.begin(),
64  finalStateIDs.end(), p) != finalStateIDs.end();
65  }
66 
70  bool TH_Channel::isChannel(str p0, str p1, str p2, str p3) const
71  {
72  str vals[] = {p0, p1};
73  std::vector<str> inIDs(std::begin(vals), std::end(vals));
74  if (p2 != "") inIDs.push_back(p2);
75  if (p3 != "") inIDs.push_back(p3);
76  return isChannel(inIDs);
77  }
78 
82  bool TH_Channel::isChannel(std::vector<str> particles) const
83  {
84  if (particles.size() != finalStateIDs.size()) return false;
85  return std::is_permutation(
86  finalStateIDs.begin(), finalStateIDs.end(), particles.begin());
87  }
88 
89 
90  // TH_Process definitions
91 
93  TH_Process::TH_Process(const str & particle1ID)
94  : isAnnihilation (false),
95  particle1ID (particle1ID),
96  particle2ID (""),
97  genRateMisc (daFunk::zero())
98  {}
99 
102  : isAnnihilation (true),
103  particle1ID(particle1ID),
104  particle2ID(particle2ID),
105  genRateMisc(daFunk::zero("v"))
106  {
107  if (particle1ID.compare(particle2ID) > 0)
108  {
109  DarkBit_error().raise(LOCAL_INFO,
110  "Particle identifiers should be in alphabetical order.");
111  }
112  }
113 
115  bool TH_Process::isProcess(const str & p1, const str & p2) const
116  {
117  sspair candidate_process(p1, p2);
118  sspair this_process(this->particle1ID, this->particle2ID);
119  return candidate_process == this_process;
120  }
121 
123  const TH_Channel* TH_Process::find(std::vector<str> final_states) const
124  {
125  for (auto it = channelList.begin(); it != channelList.end(); ++it)
126  {
127  if (it->isChannel(final_states)) return &(*it);
128  }
129  return NULL;
130  }
131 
132 
133  // TH_ProcessCatalog definitions
134 
137  {
138  const TH_Process* temp = find(id1, id2);
139  if (temp == NULL)
140  {
141  DarkBit_error().raise(LOCAL_INFO,
142  "Process with initial states " + id1 + " " + id2 + " missing.");
143  }
144  return *temp;
145  }
146 
148  const TH_Process* TH_ProcessCatalog::find(str id1, str id2) const
149  {
150  for (std::vector<TH_Process>::const_iterator it = processList.begin();
151  it != processList.end(); ++it)
152  {
153  if ( it -> isProcess(id1, id2) ) return &(*it);
154  }
155  return NULL;
156  }
157 
162  {
163  auto it = particleProperties.find(id);
164  if ( it == particleProperties.end() )
165  {
166  DarkBit_error().raise(LOCAL_INFO,
167  "Cannot find entry for " + id +
168  " in particleProperties of TH_ProcessCatalog.");
169  }
170  return it->second;
171  }
172 
174  {
175  auto it = particleProperties.find(id);
176  return (it != particleProperties.end());
177  }
178 
179 
181  {
182 #ifdef DARKBIG_DEBUG
183  std::cout << std::endl;
184  std::cout << "********************************" << std::endl;
185  std::cout << "*** Validate Process catalog ***" << std::endl;
186  std::cout << "********************************" << std::endl;
187  std::cout << std::endl;
188 #endif
189 
190  for (auto it = processList.begin(); it != processList.end(); it++)
191  {
192  if (it->isAnnihilation) // Annihilation
193  {
194  std::string p1ID = it->particle1ID;
195  std::string p2ID = it->particle2ID;
196  double m1 = getParticleProperty(p1ID).mass;
197  double m2 = getParticleProperty(p2ID).mass;
198  double E_in = m1 + m2;
199 #ifdef DARKBIT_DEBUG
200  std::cout << "Annihilation process for: " << p1ID << " " << p2ID << std::endl;
201  std::cout << " Masses (GeV): " << m1 << " " << m2 << std::endl;
202 #endif
203 
204  if ((it->genRateMisc->getNArgs() != 1) or not(it->genRateMisc->hasArg("v")))
205  DarkBit_error().raise(LOCAL_INFO,
206  "Invalid TH_ProcessCatalog annihilation entry for " + it->particle1ID + " " + it->particle2ID + "\n"
207  " genRateMisc must have relative velocity v as only argument.");
208 
209 #ifdef DARKBIT_DEBUG
210  std::cout << " genRateMisc(v=0) = " << it->genRateMisc->bind("v")->eval(0) << std::endl;
211 #endif
212 
213  for (auto it2 = it->channelList.begin(); it2 != it->channelList.end(); it2++)
214  {
215  double E_out = 0;
216  std::string outstring = "";
217  for (size_t i = 0; i < it2->finalStateIDs.size(); i++)
218  {
219  E_out += getParticleProperty(it2->finalStateIDs[i]).mass;
220  outstring += it2->finalStateIDs[i] + " ";
221  }
222 #ifdef DARKBIT_DEBUG
223  std::cout << " Channel: " << outstring << std::endl;
224  std::cout << " Total mass: " << E_out << std::endl;
225 #endif
226  if (it2->nFinalStates != it2->finalStateIDs.size())
227  DarkBit_error().raise(LOCAL_INFO,
228  "Invalid TH_ProcessCatalog annihilation entry for " + it->particle1ID + " " + it->particle2ID + "\n"
229  " inconsistent nFinalStates entry for final states " + outstring);
230  if ((it2->finalStateIDs.size() < 2) or (it2->finalStateIDs.size() > 3))
231  DarkBit_error().raise(LOCAL_INFO,
232  "Invalid TH_ProcessCatalog annihilation entry for " + it->particle1ID + " " + it->particle2ID + "\n"
233  " number of final states not supported for annihilation into " + outstring);
234  if (it2->finalStateIDs.size() == 2)
235  if ((it2->genRate->getNArgs() != 1) or not it2->genRate->hasArg("v"))
236  DarkBit_error().raise(LOCAL_INFO,
237  "Invalid TH_ProcessCatalog annihilation entry for " + it->particle1ID + " " + it->particle2ID + "\n"
238  " genRate must have relative velocity v as only argument for annihilation into " + outstring);
239  if (it2->finalStateIDs.size() == 3)
240  if (it2->genRate->getNArgs() != 3 or not it2->genRate->hasArg("E") or not it2->genRate->hasArg("E1") or not it2->genRate->hasArg("v"))
241  DarkBit_error().raise(LOCAL_INFO,
242  "Invalid TH_ProcessCatalog annihilation entry for " + it->particle1ID + " " + it->particle2ID + "\n"
243  " genRate must have three arguments (v, E and E1) for annihilation into " + outstring);
244  if ((it2->finalStateIDs.size() == 3) and E_out > E_in)
245  {
246  // TODO: This could be defined more generally (allowing three-body final states that are closed for v=0)
247  DarkBit_error().raise(LOCAL_INFO,
248  "Invalid TH_ProcessCatalog annihilation entry for " + it->particle1ID + " " + it->particle2ID + "\n"
249  " three-body final states kinematically not allowed for v=0 for annihilation into " + outstring);
250  }
251  // Test whether channels kinematically closed for v=0 are indeed closed
252  if ((it2->finalStateIDs.size() == 2) and (E_out > E_in))
253  {
254  double v_min = sqrt(1-1/pow(E_out/E_in,2));
255  if ( it2->genRate->bind("v")->eval(v_min*0.999) != 0 )
256  DarkBit_error().raise(LOCAL_INFO,
257  "Invalid TH_ProcessCatalog annihilation entry for " + it->particle1ID + " " + it->particle2ID + "\n"
258  " genRate must be zero for values of v that are below kinematic threshold for annihilation into " + outstring);
259  }
260 #ifdef DARKBIT_DEBUG
261  if (it2->finalStateIDs.size() == 2)
262  std::cout << " genRate(v=0) = " << it2->genRate->bind("v")->eval() << std::endl;
263  else
264  std::cout << " genRate(v=0) = f(E, E1)" << std::endl;
265 #endif
266  }
267  }
268  else // Decay
269  {
270  double E_in = getParticleProperty(it->particle1ID).mass;
271 
272 #ifdef DARKBIT_DEBUG
273  std::cout << "Decay process for: " << it->particle1ID << std::endl;
274  std::cout << " Mass (GeV): " << E_in << std::endl;
275  std::cout << " genRateMisc = " << it->genRateMisc->bind()->eval() << std::endl;
276 #endif
277 
278  if (it->genRateMisc->getNArgs() != 0)
279  DarkBit_error().raise(LOCAL_INFO,
280  "Invalid TH_ProcessCatalog decay entry for " + it->particle1ID + "\n"
281  " genRateMisc must have zero arguments.");
282  for (auto it2 = it->channelList.begin(); it2 != it->channelList.end(); it2++)
283  {
284  double E_out = 0;
285  std::string outstring = "";
286  for (size_t i = 0; i < it2->finalStateIDs.size(); i++)
287  {
288  E_out += getParticleProperty(it2->finalStateIDs[i]).mass;
289  outstring += it2->finalStateIDs[i] + " ";
290  }
291 #ifdef DARKBIT_DEBUG
292  std::cout << " Channel: " << outstring << std::endl;
293  std::cout << " Total mass: " << E_out << std::endl;
294 #endif
295  if (it2->nFinalStates != it2->finalStateIDs.size())
296  DarkBit_error().raise(LOCAL_INFO,
297  "Invalid TH_ProcessCatalog decay entry for " + it->particle1ID + "\n"
298  " inconsistent nFinalStates entry for final states " + outstring);
299  if ((it2->finalStateIDs.size() < 2) or (it2->finalStateIDs.size() > 3))
300  DarkBit_error().raise(LOCAL_INFO,
301  "Invalid TH_ProcessCatalog decay entry for " + it->particle1ID + "\n"
302  " number of final states not supported for decay into " + outstring);
303  if (it2->finalStateIDs.size() == 2)
304  if (it2->genRate->getNArgs() != 0)
305  DarkBit_error().raise(LOCAL_INFO,
306  "Invalid TH_ProcessCatalog decay entry for " + it->particle1ID + "\n"
307  " genRate must have zero arguments for decay into " + outstring);
308  if (it2->finalStateIDs.size() == 3)
309  if (it2->genRate->getNArgs() != 2 or not it2->genRate->hasArg("E") or not it2->genRate->hasArg("E1"))
310  DarkBit_error().raise(LOCAL_INFO,
311  "Invalid TH_ProcessCatalog decay entry for " + it->particle1ID + "\n"
312  " genRate must have two arguments (E and E1) for decay into " + outstring);
313  if (E_out > E_in)
314  DarkBit_error().raise(LOCAL_INFO,
315  "Invalid TH_ProcessCatalog decay entry for " + it->particle1ID + "\n"
316  " kinematically forbidden decay into " + outstring);
317 #ifdef DARKBIT_DEBUG
318  if (it2->finalStateIDs.size() == 2)
319  std::cout << " genRate = " << it2->genRate->bind()->eval() << std::endl;
320  else
321  std::cout << " genRate = f(E, E1)" << std::endl;
322 #endif
323  }
324  }
325  }
326 #ifdef DARKBIT_DEBUG
327  std::cout << std::endl;
328  std::cout << "*****************" << std::endl;
329  std::cout << "*** Validated ***" << std::endl;
330  std::cout << "*****************" << std::endl;
331  std::cout << std::endl;
332 #endif
333  }
334  } // namespace DarkBit
335 } // namespace Gambit
336 
337 #undef DARKBIT_DEBUG
std::vector< TH_Channel > channelList
List of channels.
error & DarkBit_error()
TH_ParticleProperty getParticleProperty(str) const
Retrieve properties of a given particle involved in one or more processes in this catalog...
void validate()
Validate kinematics and entries.
bool isChannel(str, str, str="", str="") const
Indicate whether or not this channel is the one defined by some specific final states. Particle name version.
#define LOCAL_INFO
Definition: local_info.hpp:34
TH_Channel(std::vector< str > finalStateIDs, daFunk::Funk genRate)
Constructor.
TH_Process(const str &particle1ID)
Constructor for decay process.
bool hasParticleProperty(str) const
Check whether particle is in particle properties catalog.
TH_Process getProcess(str, str="") const
Retrieve a specific process from the catalog.
std::pair< str, str > sspair
Shorthand for a pair of standard strings.
Definition: util_types.hpp:64
unsigned int nFinalStates
Number of final state particles in this channel.
const TH_Channel * find(std::vector< str >) const
Check for given channel. Return a pointer to it if found, NULL if not.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
bool channelContains(str p) const
Indicate whether or not the final states of this channel contain a specific particle.
Header file that includes all GAMBIT headers required for a module source file.
A container for the mass and spin of a particle.
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
Type definition header for DarkBit Process Catalog constituents types.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
TH_ParticleProperty(double mass, unsigned int spin2)
Constructor.
All data on a single annihilation or decay channel, e.g.
daFunk::Funk genRateMisc
Additional decay rate or sigmav (in addition to above channels)
bool isAnnihilation
Annihilation or decay?
std::vector< std::string > finalStateIDs
Final state identifiers.
Rollcall header for module DarkBit.
void printChannel() const
Print information about this channel.
double pow(const double &a)
Outputs a^i.
Funk zero(Args... argss)
Definition: daFunk.hpp:604
A container for a single process.
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
const TH_Process * find(str, str="") const
Check for a specific process in the catalog.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
str particle1ID
Decaying particle or particle pair.
bool isProcess(const str &, const str &=std::string()) const
Compare initial states.