gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
Cascades.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
20 
23 
24 //#define DARKBIT_DEBUG
25 
26 namespace Gambit
27 {
28  namespace DarkBit
29  {
30 
32  //
33  // Cascade Decays
34  //
36 
39 
41  void cascadeMC_FinalStates(std::vector<std::string> &list)
42  {
43  list.clear();
44  list.push_back("gamma");
45  using namespace Pipes::cascadeMC_FinalStates;
47  list = runOptions->getValueOrDef<std::vector<std::string>>(list, "cMC_finalStates");
48  #ifdef DARKBIT_DEBUG
49  std::cout << "Final states to generate: " << list.size() << std::endl;
50  for(size_t i=0; i < list.size(); i++)
51  {
52  std::cout << " " << list[i] << std::endl;
53  }
54  #endif
55  }
56 
59  {
60  using namespace DecayChain;
61  using namespace Pipes::cascadeMC_DecayTable;
62  std::set<std::string> disabled;
63  // Note: One could add to "disabled" particles decays in that are in the
64  // process catalog but should for some reason not propagate to the FCMC
65  // DecayTable.
66  try
67  {
68  table = DecayTable(*Dep::TH_ProcessCatalog, *Dep::SimYieldTable, disabled);
69  }
71  {
72  DarkBit_error().raise(err.first,err.second);
73  }
74 #ifdef DARKBIT_DEBUG
75  table.printTable();
76 #endif
77  }
78 
81  {
82  using namespace Pipes::cascadeMC_LoopManager;
83  std::vector<std::string> chainList = *Dep::GA_missingFinalStates;
84  int cMC_minEvents = 2; // runOptions->getValueOrDef<int>(2, "cMC_minEvents");
85  // Get YAML options
87  int cMC_maxEvents = runOptions->getValueOrDef<int>(20000, "cMC_maxEvents");
88 
89 
90  // Initialization run
91  Loop::executeIteration(MC_INIT);
92 
93  // Check whether there is anything to do
94  if ( chainList.size() == 0 )
95  {
96  return;
97  }
98 
99  // Iterate over initial state particles
100  for(std::vector<std::string>::const_iterator
101  cit =chainList.begin(); cit != chainList.end(); cit++)
102  {
103  int it;
104  int counter = 0;
105  bool finished = false;
106  // Set next initial state
107  Loop::executeIteration(MC_NEXT_STATE);
108  // Event generation loop
109  #pragma omp parallel private(it) shared(counter, finished)
110  {
111  while (!finished)
112  {
113  #pragma omp critical (cascadeMC_Counter)
114  {
115  counter++;
116  it = counter;
117  }
118  Loop::executeIteration(it);
119  #pragma omp critical (cascadeMC_Counter)
120  {
121  if((*Loop::done and ((counter >= cMC_minEvents) or piped_errors.inquire()))
122  or (counter >= cMC_maxEvents))
123  finished=true;
124  if (counter >= cMC_maxEvents)
125  DarkBit_warning().raise(LOCAL_INFO,
126  "WARNING FCMC: cMC_maxEvents reached without convergence.");
127  }
128  }
129  }
130  // Raise any exceptions
134  Loop::reset();
135  }
136  Loop::executeIteration(MC_FINALIZE);
137  }
138 
140  void cascadeMC_InitialState(std::string &pID)
141  {
142  using namespace DecayChain;
143  using namespace Pipes::cascadeMC_InitialState;
144  std::vector<std::string> chainList = *Dep::GA_missingFinalStates;
145  static int iteration;
146  switch(*Loop::iteration)
147  {
148  case MC_INIT:
149  iteration = -1;
150  return;
151  case MC_NEXT_STATE:
152  iteration++;
153  break;
154  case MC_FINALIZE:
155  return;
156  }
157  if(int(chainList.size()) <= iteration)
158  {
159  Loop::wrapup();
161  "Desync between cascadeMC_LoopManager and cascadeMC_InitialState");
162  }
163  else
164  pID = chainList[iteration];
165 #ifdef DARKBIT_DEBUG
166  std::cout << "cascadeMC_InitialState" << std::endl;
167  std::cout << " Iteration: " << *Loop::iteration << std::endl;
168  std::cout << " Number of states to simulate: "
169  << chainList.size() << std::endl;
170  std::cout << " Current state: " << pID << std::endl;
171 #endif
172  }
173 
175  void cascadeMC_EventCount(std::map<std::string, int> &counts)
176  {
177  using namespace Pipes::cascadeMC_EventCount;
178  static std::map<std::string, int> counters;
179  switch(*Loop::iteration)
180  {
181  case MC_INIT:
182  counters.clear();
183  break;
184  case MC_NEXT_STATE:
185  counters[*Dep::cascadeMC_InitialState] = 0;
186  break;
187  case MC_FINALIZE:
188  // For performance, only return the actual result once finished
189  counts=counters;
190  break;
191  default:
192 #pragma omp atomic
193  counters[*Dep::cascadeMC_InitialState]++;
194  }
195  }
196 
200  {
201  using namespace DecayChain;
202  using namespace Pipes::cascadeMC_GenerateChain;
203  static int cMC_maxChainLength;
204  static double cMC_Emin;
205  switch(*Loop::iteration)
206  {
207  case MC_INIT:
209  cMC_maxChainLength = runOptions->getValueOrDef<int> (-1, "cMC_maxChainLength");
211  cMC_Emin = runOptions->getValueOrDef<double> (-1, "cMC_Emin");
212  return;
213  case MC_NEXT_STATE:
214  case MC_FINALIZE:
215  return;
216  }
217  shared_ptr<ChainParticle> chn;
218  try
219  {
220  chn.reset(new ChainParticle( vec3(0), &(*Dep::cascadeMC_DecayTable),
222  chn->generateDecayChainMC(cMC_maxChainLength,cMC_Emin);
223  }
225  {
226  Loop::wrapup();
227  piped_errors.request(err);
228  }
229  chain=ChainContainer(chn);
230  }
231 
236  const DarkBit::DecayChain::ChainParticle* endpoint,
237  std::string finalState,
238  const TH_ProcessCatalog &catalog,
239  std::map<std::string, std::map<std::string, SimpleHist> > &histList,
240  std::string initialState,
241  double weight, int cMC_numSpecSamples
242  )
243  {
244 #ifdef DARKBIT_DEBUG
245  std::cout << "SampleSimYield" << std::endl;
246 #endif
247  std::string p1,p2;
248  double gamma,beta;
249  double M;
250  switch(endpoint->getnChildren())
251  {
252  case 0:
253  {
254  p1 = endpoint->getpID();
255  p2 = "";
256  const DarkBit::DecayChain::ChainParticle* parent = endpoint->getParent();
257  if(parent == NULL)
258  {
259  endpoint->getBoost(gamma,beta);
260  M = endpoint->m;
261  }
262  else
263  {
264  parent->getBoost(gamma,beta);
265  M = endpoint->E_parentFrame();
266  }
267  break;
268  }
269  case 2:
270  {
271  p1=(*endpoint)[0]->getpID();
272  p2=(*endpoint)[1]->getpID();
273  endpoint->getBoost(gamma,beta);
274  M = endpoint->m;
275  break;
276  }
277  default:
279  "cascadeMC_sampleSimYield called with invalid endpoint state.");
280  return;
281  }
282  const SimYieldChannel &chn = table.getChannel(p1 , p2, finalState);
283  // Get Lorentz boost information
284 
285  const double gammaBeta = gamma*beta;
286  // Mass of final state squared
287  const double m = catalog.getParticleProperty(finalState).mass;
288  const double msq = m*m;
289  // Get histogram edges
290  double histEmin, histEmax;
291  histList[initialState][finalState].getEdges(histEmin, histEmax);
292 
293  // Calculate energies to sample between. A particle decaying
294  // isotropically in its rest frame will give a box spectrum. This is
295  // assumed and used here to add box contributions rather than points to
296  // the histograms. Limits are chosen such that we only sample energies
297  // that can contribute to histogram bins.
298  const double Ecmin = std::max( gamma*histEmin
299  - gammaBeta*sqrt(histEmin*histEmin-msq) , 0*chn.Ecm_min ); // CW: chn.Ecm_min refers to initial not final energies
300  const double Ecmax = std::min(std::min(
301  // Highest energy that can contribute to the histogram
302  gamma*histEmax + gammaBeta*sqrt(histEmax*histEmax-msq),
303  // Highest energy in SimYieldChannel
304  chn.Ecm_max ),
305  // Estimate for highest kinematically allowed CoM energy
306  0.5*(M*M + msq)/M );
307  if(Ecmin>=Ecmax) return;
308  const double logmin = log(Ecmin);
309  const double logmax = log(Ecmax);
310  const double dlogE=logmax-logmin;
311 
312 #ifdef DARKBIT_DEBUG
313  std::cout << "M = " << M << std::endl;
314  std::cout << "E_lab = " << endpoint->E_Lab() << std::endl;
315  std::cout << "p_lab = " << endpoint->p_Lab() << std::endl;
316  std::cout << "Lorentz factors gamma, beta: " << gamma << ", "
317  << beta << std::endl;
318  std::cout << "Initial state: " << initialState << std::endl;
319  std::cout << "Channel: " << p1 << " " << p2 << std::endl;
320  std::cout << "Final particles: " << finalState << std::endl;
321  std::cout << "Event weight: " << weight << std::endl;
322  std::cout << "histEmin/histEmax: " << histEmin << " " << histEmax
323  << std::endl;
324  std::cout << "chn.Ecm_min/max: " << chn.Ecm_min << " " << chn.Ecm_max
325  << std::endl;
326  std::cout << "Ecmin/max: " << Ecmin << " " << Ecmax << std::endl;
327  std::cout << "Final state mass^2: " << msq << std::endl;
328 #endif
329 
330  double specSum=0;
331  int Nsampl=0;
332  SimpleHist spectrum(histList[initialState][finalState].binLower);
333  while(Nsampl<cMC_numSpecSamples)
334  {
335  // Draw an energy in the CoM frame of the endpoint. Logarithmic
336  // sampling.
337  double E_CoM= exp(logmin+(logmax-logmin)*Random::draw());
338  double dN_dE = chn.dNdE_bound->eval(E_CoM, M);
339 
340  double weight2 = E_CoM*dlogE*dN_dE;
341  specSum += weight2;
342  weight2 *= weight;
343  double tmp1 = gamma*E_CoM;
344  double tmp2 = gammaBeta*sqrt(E_CoM*E_CoM-msq);
345  // Add box spectrum to histogram
346  spectrum.addBox(tmp1-tmp2,tmp1+tmp2,weight2);
347  Nsampl++;
348  }
349 #ifdef DARKBIT_DEBUG
350  std::cout << "Number of samples = " << Nsampl << std::endl;
351 #endif
352  if(Nsampl>0)
353  {
354  spectrum.multiply(1.0/Nsampl);
355  // Add bin contents of spectrum histogram to main histogram as weighted
356  // events
357  #pragma omp critical (cascadeMC_histList)
358  histList[initialState][finalState].addHistAsWeights_sameBin(spectrum);
359  }
360  }
361 
363  void cascadeMC_Histograms(std::map<std::string, std::map<std::string,
364  SimpleHist> > &result)
365  {
366  using namespace DecayChain;
367  using namespace Pipes::cascadeMC_Histograms;
368 
369  // YAML options
370  static int cMC_numSpecSamples;
371  static int cMC_endCheckFrequency;
372  static double cMC_gammaBGPower;
373  static double cMC_gammaRelError;
374  static int cMC_NhistBins;
375  static double cMC_binLow;
376  static double cMC_binHigh;
377  // Histogram list shared between all threads
378  static std::map<std::string, std::map<std::string, SimpleHist> > histList;
379 
380  switch(*Loop::iteration)
381  {
382  case MC_INIT:
383  // Initialization
386  cMC_numSpecSamples = runOptions->getValueOrDef<int> (25, "cMC_numSpecSamples");
389  cMC_endCheckFrequency =
390  runOptions->getValueOrDef<int> (25, "cMC_endCheckFrequency");
393  cMC_gammaBGPower =
394  runOptions->getValueOrDef<double>(-2.5, "cMC_gammaBGPower");
397  cMC_gammaRelError =
398  runOptions->getValueOrDef<double>(0.20, "cMC_gammaRelError");
399 
400  // Note: use same binning for all particle species
402  cMC_NhistBins = runOptions->getValueOrDef<int> (140, "cMC_NhistBins");
404  cMC_binLow = runOptions->getValueOrDef<double>(0.001, "cMC_binLow");
406  cMC_binHigh = runOptions->getValueOrDef<double>(10000.0,"cMC_binHigh");
407  histList.clear();
408  return;
409  case MC_NEXT_STATE:
410  // Initialize histograms
411  for(std::vector<std::string>::const_iterator it =
413  it!=Dep::cascadeMC_FinalStates->end(); ++it)
414  {
415 #ifdef DARKBIT_DEBUG
416  std::cout << "Defining new histList entry!!!" << std::endl;
417  std::cout << "for: " << *Dep::cascadeMC_InitialState
418  << " " << *it << std::endl;
419 #endif
420  histList[*Dep::cascadeMC_InitialState][*it]=
421  SimpleHist(cMC_NhistBins,cMC_binLow,cMC_binHigh,true);
422  }
423  return;
424  case MC_FINALIZE:
425  // For performance, only return the actual result once finished
426  result = histList;
427  return;
428  }
429 
430  // Get list of endpoint states for this chain
431  vector<const ChainParticle*> endpoints;
432  (*Dep::cascadeMC_ChainEvent).chain->
433  collectEndpointStates(endpoints, false);
434  // Iterate over final states of interest
435  for(std::vector<std::string>::const_iterator pit =
437  pit!=Dep::cascadeMC_FinalStates->end(); ++pit)
438  {
439  // Iterate over all endpoint states of the decay chain. These can
440  // either be final state particles themselves or parents of final state
441  // particles. The reason for not using only final state particles is
442  // that certain endpoints (e.g. quark-antiquark pairs) cannot be
443  // handled as separate particles.
444  for(vector<const ChainParticle*>::const_iterator it =endpoints.begin();
445  it != endpoints.end(); it++)
446  {
447 #ifdef DARKBIT_DEBUG
448  std::cout << " working on endpoint: " << (*it)->getpID() << std::endl;
449  (*it)->printChain();
450 #endif
451 
452  // Weighting factor (correction for mismatch between decay width
453  // of available decay channels and total decay width)
454  double weight;
455  // Analyze single particle endpoints
456  bool ignored = true;
457  if((*it)->getnChildren() ==0)
458  {
459  weight = (*it)->getWeight();
460  // Check if the final state itself is the particle we are looking
461  // for.
462  if((*it)->getpID()==*pit)
463  {
464  double E = (*it)->E_Lab();
465 #pragma omp critical (cascadeMC_histList)
466  histList[*Dep::cascadeMC_InitialState][*pit].addEvent(E,weight);
467  ignored = false;
468  }
469  // Check if tabulated spectra exist for this final state
470  else if((*Dep::SimYieldTable).hasChannel( (*it)->getpID(), *pit ))
471  {
473  *Dep::SimYieldTable, *it, *pit, *Dep::TH_ProcessCatalog,
474  histList, *Dep::cascadeMC_InitialState, weight,
475  cMC_numSpecSamples
476  );
477  // Check if an error was raised
478  ignored = false;
479  if(piped_errors.inquire())
480  {
481  Loop::wrapup();
482  return;
483  }
484  }
485  }
486  // Analyze multiparticle endpoints (the endpoint particle is here the
487  // parent of final state particles).
488  else
489  {
490  weight = (*(*it))[0]->getWeight();
491  bool hasTabulated = false;
492  if((*it)->getnChildren() == 2)
493  {
494 #ifdef DARKBIT_DEBUG
495  std::cout << " check whether two-body final state is tabulated: "
496  << (*(*it))[0]->getpID() << " " << (*(*it))[1]->getpID() <<
497  std::endl;
498 #endif
499  // Check if tabulated spectra exist for this final state
500  if((*Dep::SimYieldTable).hasChannel(
501  (*(*it))[0]->getpID() , (*(*it))[1]->getpID(), *pit ))
502  {
503  hasTabulated = true;
504  cascadeMC_sampleSimYield(*Dep::SimYieldTable, *it, *pit,
505  *Dep::TH_ProcessCatalog, histList,
507  cMC_numSpecSamples
508  );
509  // Check if an error was raised
510  ignored = false;
511  if(piped_errors.inquire())
512  {
513  Loop::wrapup();
514  return;
515  }
516  }
517  }
518  if(!hasTabulated)
519  {
520  for(int i=0; i<((*it)->getnChildren()); i++)
521  {
522  const ChainParticle* child = (*(*it))[i];
523  // Check if the child particle is the particle we are looking
524  // for.
525  if(child->getpID()==*pit)
526  {
527  double E = child->E_Lab();
528 #pragma omp critical (cascadeMC_histList)
529  histList[
530  *Dep::cascadeMC_InitialState][*pit].addEvent(E,weight);
531  ignored = false;
532  }
533  // Check if tabulated spectra exist for this final state
534  else if((*Dep::SimYieldTable).hasChannel( child->getpID(),
535  *pit))
536  {
537  cascadeMC_sampleSimYield(*Dep::SimYieldTable, child, *pit,
538  *Dep::TH_ProcessCatalog, histList,
540  cMC_numSpecSamples
541  );
542  // Check if an error was raised
543  ignored = false;
544  if(piped_errors.inquire())
545  {
546  Loop::wrapup();
547  return;
548  }
549  }
550  }
551  }
552  }
553  if (ignored)
554  {
555  DarkBit_warning().raise(LOCAL_INFO,
556  "WARNING FCMC: Missing complete decay information for "
557  + (*it)->getpID() + ". This state is ignored.");
558  }
559  }
560  }
561  // Check if finished every cMC_endCheckFrequency events
562  if((*Loop::iteration % cMC_endCheckFrequency) == 0)
563  {
564  enum status{untouched,unfinished,finished};
565  status cond = untouched;
566  for(std::vector<std::string>::const_iterator it =
568  it != Dep::cascadeMC_FinalStates->end(); ++it)
569  {
570  // End conditions currently only implemented for gamma final state
571  if(*it=="gamma")
572  {
573  SimpleHist hist;
574 #pragma omp critical (cascadeMC_histList)
575  hist = histList[*Dep::cascadeMC_InitialState][*it];
576 #ifdef DARKBIT_DEBUG
577  std::cout << "Checking whether convergence is reached" << std::endl;
578  for ( int i = 0; i < hist.nBins; i++ )
579  std::cout << "Estimated error at " << hist.binCenter(i) << " GeV : " << hist.getRelError(i) << std::endl;
580 #endif
581  double sbRatioMax=-1.0;
582  int maxBin=0;
583  for(int i=0; i<hist.nBins; i++)
584  {
585  double E = hist.binCenter(i);
586  double background = pow(E,cMC_gammaBGPower);
587  double sbRatio = hist.binVals[i]/background;
588  if(sbRatio>sbRatioMax)
589  {
590  sbRatioMax = sbRatio;
591  maxBin=i;
592  }
593  }
594 #ifdef DARKBIT_DEBUG
595  std::cout << "Estimated maxBin: " << maxBin << std::endl;
596  std::cout << "Energy at maxBin: " << hist.binCenter(maxBin) << std::endl;
597  std::cout << "Estimated error at maxBin: " << hist.getRelError(maxBin) << std::endl;
598  std::cout << "Value at maxBin: " << hist.getBinValues()[maxBin];
599 #endif
600  // Check if end condition is fulfilled. If not, set cond to
601  // unfinished.
602  if(hist.getRelError(maxBin) > cMC_gammaRelError) cond = unfinished;
603 
604  // If end condition is fulfilled, set cond to finished, unless
605  // already set to unfinished by another condition.
606  else if(cond != unfinished) cond = finished;
607  }
608  }
609  // Break Monte Carlo loop if all end conditions are fulfilled.
610  if(cond==finished)
611  {
612 #ifdef DARKBIT_DEBUG
613  std::cout << "!! wrapping up !!" << std::endl;
614  std::cout << "Performed iterations: " << *Loop::iteration << std::endl;
615 #endif
616  Loop::wrapup();
617  }
618  }
619  }
620 
625  void cascadeMC_fetchSpectra(std::map<std::string, daFunk::Funk> &spectra,
626  std::string finalState,
627  const std::vector<std::string> &ini,
628  const std::vector<std::string> &fin,
629  const std::map<std::string, std::map<std::string,SimpleHist> > &h,
630  const std::map<std::string,int> &eventCounts)
631  {
632  spectra.clear();
633  // Check if final state has been calculated
634  bool calculated = (
635  std::find(fin.begin(), fin.end(), finalState) != fin.end());
636  for(std::vector<std::string>::const_iterator it = ini.begin();
637  it != ini.end(); ++it )
638  {
639 #ifdef DARKBIT_DEBUG
640  std::cout << "Trying to get cascade spectra for initial state: " << *it << std::endl;
641 #endif
642  if(calculated)
643  {
644 #ifdef DARKBIT_DEBUG
645  std::cout << finalState << "...was calculated!" << std::endl;
646  std::cout << eventCounts.at(*it) << " events generated" << std::endl;
647 #endif
648  SimpleHist hist = h.at(*it).at(finalState);
649  hist.divideByBinSize();
650  std::vector<double> E = hist.getBinCenters();
651  std::vector<double> dN_dE = hist.getBinValues();
652  // Normalize to per-event spectrum
653  int i = 0;
654  for (std::vector<double>::iterator it2=dN_dE.begin();
655  it2!=dN_dE.end();++it2)
656  {
657  *it2 /= eventCounts.at(*it);
658 #ifdef DARKBIT_DEBUG
659  std::cout << E[i] << " " << *it2 << std::endl;
660 #endif
661  i++;
662  }
663  // Default values provide 1-2% accuracy for singular integrals
664  // Make this optional.
665  spectra[*it] = daFunk::Funk(new daFunk::FunkInterp("E", E, dN_dE, "lin"));
666 
667  for (size_t i = 1; i<E.size()-1; i++)
668  {
669  if (dN_dE[i]/(dN_dE[i-1]+dN_dE[i+1]+dN_dE[i]*1e-4) > 1e2)
670  {
671 #ifdef DARKBIT_DEBUG
672  std::cout << "Set singularity at " << E[i] << " with width "
673  << E[i+1]-E[i] << endl;
674 #endif
675  spectra[*it]->set_singularity("E", E[i], (E[i+1]-E[i]));
676  }
677  }
678  }
679  else
680  {
681  spectra[*it] = daFunk::zero("E");
682  }
683  }
684  }
685 
687  void cascadeMC_gammaSpectra(std::map<std::string, daFunk::Funk> &spectra)
688  {
689  using namespace Pipes::cascadeMC_gammaSpectra;
693 #ifdef DARKBIT_DEBUG
694  std::cout << "Retrieving cascade spectra for gamma final states" << std::endl;
695  std::cout << "Number of simulated final states: " << spectra.size() << std::endl;
696  for ( auto it = spectra.begin(); it != spectra.end(); it ++ )
697  {
698  std::cout << "Particle: " << it->first << std::endl;
699  auto f= it->second;
700  for ( double E = 0.1; E < 1000; E*=1.5 )
701  {
702  std::cout << " " << E << " " << f->bind("E")->eval(E) << std::endl;
703  }
704  std::cout << " Integrated spectrum: " << f->gsl_integration("E", 0, 1000)->bind()->eval() << std::endl;
705  }
706 #endif
707  }
708 
709  /*
710  void cascadeMC_PrintResult(bool &dummy)
711  {
712  dummy=true;
713  using namespace Pipes::cascadeMC_PrintResult;
714  logger() << "************************" << std::endl;
715  logger() << "Cascade decay results:" << std::endl;
716  logger() << "------------------------" << EOM;
717  std::map<std::string, std::map<std::string,SimpleHist> >
718  cascadeMC_HistList = *Dep::cascadeMC_Histograms;
719 
720  for(std::map<std::string, std::map<std::string,SimpleHist> >::iterator
721  it = cascadeMC_HistList.begin();
722  it != cascadeMC_HistList.end(); ++it )
723  {
724  logger() << "Initial state: " << (it->first) << ":" << EOM;
725  int nEvents = (*Dep::cascadeMC_EventCount).at(it->first);
726  logger() << "Number of events: " << nEvents << EOM;
727  for(std::map<std::string,SimpleHist>::iterator
728  it2 = (it->second).begin(); it2 != (it->second).end(); ++it2 )
729  {
730  logger() << (it2->first) << ": ";
731  //(it2->second).divideByBinSize();
732  (it2->second).multiply(1.0/nEvents);
733  for(int i=0;i<50;i++)
734  {
735  logger() << (it2->second).binVals[i] << " ";
736  }
737  logger() << std::endl;
738  }
739  logger() << "------------------------" << std::endl;
740  }
741  logger() << "************************" << EOM;
742  }
743  */
744  }
745 }
746 
747 #undef DARKBIT_DEBUG
error & DarkBit_error()
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
void cascadeMC_FinalStates(std::vector< std::string > &list)
Function for retrieving list of final states for cascade decays.
Definition: Cascades.cpp:41
void GA_missingFinalStates(std::vector< std::string > &result)
Identification of final states that are not yet tabulated.
Definition: GamYields.cpp:61
void divideByBinSize()
Divide all histogram bins by the respective bin size.
Definition: SimpleHist.cpp:173
TH_ParticleProperty getParticleProperty(str) const
Retrieve properties of a given particle involved in one or more processes in this catalog...
double getRelError(int bin) const
Get relative error for a specified bin.
Definition: SimpleHist.cpp:168
std::vector< double > binVals
Definition: SimpleHist.hpp:93
std::pair< std::string, std::string > description
Definition: exceptions.hpp:278
START_MODEL beta
Definition: Axions.hpp:35
void cascadeMC_fetchSpectra(std::map< std::string, daFunk::Funk > &spectra, std::string finalState, const std::vector< std::string > &ini, const std::vector< std::string > &fin, const std::map< std::string, std::map< std::string, SimpleHist > > &h, const std::map< std::string, int > &eventCounts)
Convenience function for getting a daFunk::Funk object of a given spectrum.
Definition: Cascades.cpp:625
void getBoost(double &gamma, double &beta) const
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:527
Histogram class for cascade decays.
Definition: SimpleHist.hpp:37
#define LOCAL_INFO
Definition: local_info.hpp:34
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
Annihilation/decay channel.
void cascadeMC_InitialState(std::string &pID)
Function selecting initial state for decay chain.
Definition: Cascades.cpp:140
double binCenter(int bin) const
Get central value of bin.
Definition: SimpleHist.cpp:204
bool inquire()
Check whether any exceptions were requested without handling them.
Definition: exceptions.cpp:587
void cascadeMC_gammaSpectra(std::map< std::string, daFunk::Funk > &spectra)
Function requesting and returning gamma ray spectra from cascade decays.
Definition: Cascades.cpp:687
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
const std::vector< double > & getBinValues() const
Definition: SimpleHist.cpp:219
double mass
Particle mass (GeV)
void cascadeMC_EventCount(std::map< std::string, int > &counts)
Event counter for cascade decays.
Definition: Cascades.cpp:175
static double draw()
Draw a single uniform random deviate from the interval (0,1) using the chosen RNG engine...
const double logmin
Minimum finite result returnable from log(double x);.
Definition: statistics.cpp:37
void cascadeMC_GenerateChain(DarkBit::DecayChain::ChainContainer &chain)
Function for generating decay chains.
Definition: Cascades.cpp:198
dictionary fin
void addBox(double Emin, double Emax, double weight=1.0)
Add a box spectrum to the histogram.
Definition: SimpleHist.cpp:94
void cascadeMC_sampleSimYield(const SimYieldTable &table, const DarkBit::DecayChain::ChainParticle *endpoint, std::string finalState, const TH_ProcessCatalog &catalog, std::map< std::string, std::map< std::string, SimpleHist > > &histList, std::string initialState, double weight, int cMC_numSpecSamples)
Function for sampling SimYieldTables (tabulated spectra).
Definition: Cascades.cpp:235
A container holding all annihilation and decay initial states relevant for DarkBit.
Header file that includes all GAMBIT headers required for a module source file.
void multiply(double x)
Multiply all bin contents by x.
Definition: SimpleHist.cpp:182
std::vector< double > getBinCenters() const
Double get central values of all bins.
Definition: SimpleHist.cpp:209
START_MODEL M
cascadeMC_SpecialEvents
Special events for event loop.
Definition: Cascades.cpp:38
Rollcall header for module DarkBit.
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:544
GAMBIT native decay table class.
Definition: decay_table.hpp:35
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:473
warning & DarkBit_warning()
double pow(const double &a)
Outputs a^i.
Channel container Object containing tabularized yields for particle decay and two-body final states...
void cascadeMC_LoopManager()
Loop manager for cascade decays.
Definition: Cascades.cpp:80
void cascadeMC_Histograms(std::map< std::string, std::map< std::string, SimpleHist > > &result)
Function responsible for histogramming, and evaluating end conditions for event loop.
Definition: Cascades.cpp:363
void cascadeMC_DecayTable(DarkBit::DecayChain::DecayTable &table)
Function setting up the decay table used in decay chains.
Definition: Cascades.cpp:58
Funk zero(Args... argss)
Definition: daFunk.hpp:604
const SimYieldChannel & getChannel(std::string p1, std::string p2, std::string finalState) const
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
TODO: see if we can use this one:
Definition: Analysis.hpp:33
double E(const double x, const double y)
const ChainParticle * getParent() const