44 list.push_back(
"gamma");
47 list = runOptions->getValueOrDef<std::vector<std::string>>(list,
"cMC_finalStates");
49 std::cout <<
"Final states to generate: " << list.size() << std::endl;
50 for(
size_t i=0; i < list.size(); i++)
52 std::cout <<
" " << list[i] << std::endl;
60 using namespace DecayChain;
62 std::set<std::string> disabled;
68 table =
DecayTable(*Dep::TH_ProcessCatalog, *Dep::SimYieldTable, disabled);
84 int cMC_minEvents = 2;
87 int cMC_maxEvents = runOptions->getValueOrDef<
int>(20000,
"cMC_maxEvents");
91 Loop::executeIteration(
MC_INIT);
94 if ( chainList.size() == 0 )
100 for(std::vector<std::string>::const_iterator
101 cit =chainList.begin(); cit != chainList.end(); cit++)
105 bool finished =
false;
109 #pragma omp parallel private(it) shared(counter, finished) 113 #pragma omp critical (cascadeMC_Counter) 118 Loop::executeIteration(it);
119 #pragma omp critical (cascadeMC_Counter) 122 or (counter >= cMC_maxEvents))
124 if (counter >= cMC_maxEvents)
126 "WARNING FCMC: cMC_maxEvents reached without convergence.");
142 using namespace DecayChain;
145 static int iteration;
146 switch(*Loop::iteration)
157 if(
int(chainList.size()) <= iteration)
161 "Desync between cascadeMC_LoopManager and cascadeMC_InitialState");
164 pID = chainList[iteration];
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;
178 static std::map<std::string, int> counters;
179 switch(*Loop::iteration)
201 using namespace DecayChain;
203 static int cMC_maxChainLength;
204 static double cMC_Emin;
205 switch(*Loop::iteration)
209 cMC_maxChainLength = runOptions->getValueOrDef<
int> (-1,
"cMC_maxChainLength");
211 cMC_Emin = runOptions->getValueOrDef<
double> (-1,
"cMC_Emin");
217 shared_ptr<ChainParticle> chn;
222 chn->generateDecayChainMC(cMC_maxChainLength,cMC_Emin);
229 chain=ChainContainer(chn);
237 std::string finalState,
239 std::map<std::string, std::map<std::string, SimpleHist> > &histList,
240 std::string initialState,
241 double weight,
int cMC_numSpecSamples
245 std::cout <<
"SampleSimYield" << std::endl;
271 p1=(*endpoint)[0]->getpID();
272 p2=(*endpoint)[1]->getpID();
279 "cascadeMC_sampleSimYield called with invalid endpoint state.");
285 const double gammaBeta = gamma*
beta;
288 const double msq = m*m;
290 double histEmin, histEmax;
291 histList[initialState][finalState].getEdges(histEmin, histEmax);
298 const double Ecmin = std::max( gamma*histEmin
299 - gammaBeta*sqrt(histEmin*histEmin-msq) , 0*chn.
Ecm_min );
300 const double Ecmax = std::min(std::min(
302 gamma*histEmax + gammaBeta*sqrt(histEmax*histEmax-msq),
307 if(Ecmin>=Ecmax)
return;
308 const double logmin = log(Ecmin);
309 const double logmax = log(Ecmax);
310 const double dlogE=logmax-
logmin;
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
324 std::cout <<
"chn.Ecm_min/max: " << chn.
Ecm_min <<
" " << chn.
Ecm_max 326 std::cout <<
"Ecmin/max: " << Ecmin <<
" " << Ecmax << std::endl;
327 std::cout <<
"Final state mass^2: " << msq << std::endl;
332 SimpleHist spectrum(histList[initialState][finalState].binLower);
333 while(Nsampl<cMC_numSpecSamples)
337 double E_CoM= exp(logmin+(logmax-logmin)*
Random::draw());
338 double dN_dE = chn.
dNdE_bound->eval(E_CoM, M);
340 double weight2 = E_CoM*dlogE*dN_dE;
343 double tmp1 = gamma*E_CoM;
344 double tmp2 = gammaBeta*sqrt(E_CoM*E_CoM-msq);
346 spectrum.
addBox(tmp1-tmp2,tmp1+tmp2,weight2);
350 std::cout <<
"Number of samples = " << Nsampl << std::endl;
357 #pragma omp critical (cascadeMC_histList) 358 histList[initialState][finalState].addHistAsWeights_sameBin(spectrum);
366 using namespace DecayChain;
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;
378 static std::map<std::string, std::map<std::string, SimpleHist> > histList;
380 switch(*Loop::iteration)
386 cMC_numSpecSamples = runOptions->getValueOrDef<
int> (25,
"cMC_numSpecSamples");
389 cMC_endCheckFrequency =
390 runOptions->getValueOrDef<
int> (25,
"cMC_endCheckFrequency");
394 runOptions->getValueOrDef<
double>(-2.5,
"cMC_gammaBGPower");
398 runOptions->getValueOrDef<
double>(0.20,
"cMC_gammaRelError");
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");
411 for(std::vector<std::string>::const_iterator it =
416 std::cout <<
"Defining new histList entry!!!" << std::endl;
418 <<
" " << *it << std::endl;
421 SimpleHist(cMC_NhistBins,cMC_binLow,cMC_binHigh,
true);
431 vector<const ChainParticle*> endpoints;
432 (*Dep::cascadeMC_ChainEvent).chain->
433 collectEndpointStates(endpoints,
false);
435 for(std::vector<std::string>::const_iterator pit =
444 for(vector<const ChainParticle*>::const_iterator it =endpoints.begin();
445 it != endpoints.end(); it++)
448 std::cout <<
" working on endpoint: " << (*it)->getpID() << std::endl;
457 if((*it)->getnChildren() ==0)
459 weight = (*it)->getWeight();
462 if((*it)->getpID()==*pit)
464 double E = (*it)->E_Lab();
465 #pragma omp critical (cascadeMC_histList) 470 else if((*Dep::SimYieldTable).hasChannel( (*it)->getpID(), *pit ))
473 *Dep::SimYieldTable, *it, *pit, *Dep::TH_ProcessCatalog,
490 weight = (*(*it))[0]->getWeight();
491 bool hasTabulated =
false;
492 if((*it)->getnChildren() == 2)
495 std::cout <<
" check whether two-body final state is tabulated: " 496 << (*(*it))[0]->getpID() <<
" " << (*(*it))[1]->getpID() <<
500 if((*Dep::SimYieldTable).hasChannel(
501 (*(*it))[0]->getpID() , (*(*it))[1]->getpID(), *pit ))
505 *Dep::TH_ProcessCatalog, histList,
520 for(
int i=0; i<((*it)->getnChildren()); i++)
522 const ChainParticle* child = (*(*it))[i];
525 if(child->getpID()==*pit)
527 double E = child->E_Lab();
528 #pragma omp critical (cascadeMC_histList) 534 else if((*Dep::SimYieldTable).hasChannel( child->getpID(),
538 *Dep::TH_ProcessCatalog, histList,
556 "WARNING FCMC: Missing complete decay information for " 557 + (*it)->getpID() +
". This state is ignored.");
562 if((*Loop::iteration % cMC_endCheckFrequency) == 0)
564 enum status{untouched,unfinished,finished};
565 status cond = untouched;
566 for(std::vector<std::string>::const_iterator it =
574 #pragma omp critical (cascadeMC_histList) 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;
581 double sbRatioMax=-1.0;
583 for(
int i=0; i<hist.
nBins; i++)
586 double background =
pow(E,cMC_gammaBGPower);
587 double sbRatio = hist.
binVals[i]/background;
588 if(sbRatio>sbRatioMax)
590 sbRatioMax = sbRatio;
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];
602 if(hist.
getRelError(maxBin) > cMC_gammaRelError) cond = unfinished;
606 else if(cond != unfinished) cond = finished;
613 std::cout <<
"!! wrapping up !!" << std::endl;
614 std::cout <<
"Performed iterations: " << *Loop::iteration << std::endl;
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)
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 )
640 std::cout <<
"Trying to get cascade spectra for initial state: " << *it << std::endl;
645 std::cout << finalState <<
"...was calculated!" << std::endl;
646 std::cout << eventCounts.at(*it) <<
" events generated" << std::endl;
654 for (std::vector<double>::iterator it2=dN_dE.begin();
655 it2!=dN_dE.end();++it2)
657 *it2 /= eventCounts.at(*it);
659 std::cout << E[i] <<
" " << *it2 << std::endl;
665 spectra[*it] =
daFunk::Funk(
new daFunk::FunkInterp(
"E", E, dN_dE,
"lin"));
667 for (
size_t i = 1; i<E.size()-1; i++)
669 if (dN_dE[i]/(dN_dE[i-1]+dN_dE[i+1]+dN_dE[i]*1e-4) > 1e2)
672 std::cout <<
"Set singularity at " << E[i] <<
" with width " 673 << E[i+1]-E[i] << endl;
675 spectra[*it]->set_singularity(
"E", E[i], (E[i+1]-E[i]));
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 ++ )
698 std::cout <<
"Particle: " << it->first << std::endl;
700 for (
double E = 0.1;
E < 1000;
E*=1.5 )
702 std::cout <<
" " <<
E <<
" " <<
f->bind(
"E")->eval(
E) << std::endl;
704 std::cout <<
" Integrated spectrum: " <<
f->gsl_integration(
"E", 0, 1000)->bind()->eval() << std::endl;
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.
void GA_missingFinalStates(std::vector< std::string > &result)
Identification of final states that are not yet tabulated.
void divideByBinSize()
Divide all histogram bins by the respective bin size.
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.
std::vector< double > binVals
std::pair< std::string, std::string > description
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.
void getBoost(double &gamma, double &beta) const
daFunk::BoundFunk dNdE_bound
void request(std::string origin, std::string message)
Request an exception.
Histogram class for cascade decays.
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.
double binCenter(int bin) const
Get central value of bin.
bool inquire()
Check whether any exceptions were requested without handling them.
void cascadeMC_gammaSpectra(std::map< std::string, daFunk::Funk > &spectra)
Function requesting and returning gamma ray spectra from cascade decays.
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
const std::vector< double > & getBinValues() const
double mass
Particle mass (GeV)
void cascadeMC_EventCount(std::map< std::string, int > &counts)
Event counter for cascade decays.
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);.
void cascadeMC_GenerateChain(DarkBit::DecayChain::ChainContainer &chain)
Function for generating decay chains.
void addBox(double Emin, double Emax, double weight=1.0)
Add a box spectrum to the histogram.
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).
A container holding all annihilation and decay initial states relevant for DarkBit.
void multiply(double x)
Multiply all bin contents by x.
std::vector< double > getBinCenters() const
Double get central values of all bins.
cascadeMC_SpecialEvents
Special events for event loop.
Rollcall header for module DarkBit.
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
GAMBIT native decay table class.
void check()
Check whether an exception was requested, and throw it if necessary.
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.
double E_parentFrame() const
void cascadeMC_Histograms(std::map< std::string, std::map< std::string, SimpleHist > > &result)
Function responsible for histogramming, and evaluating end conditions for event loop.
void cascadeMC_DecayTable(DarkBit::DecayChain::DecayTable &table)
Function setting up the decay table used in decay chains.
const SimYieldChannel & getChannel(std::string p1, std::string p2, std::string finalState) const
shared_ptr< FunkBase > Funk
TODO: see if we can use this one:
double E(const double x, const double y)
const ChainParticle * getParent() const