gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
DarkBit_standalone_WIMP.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
22 
23 #include <iostream>
24 #include <fstream>
25 
30 
31 #include <boost/multi_array.hpp>
32 
33 using namespace DarkBit::Functown; // Functors wrapping the module's actual module functions
34 using namespace BackendIniBit::Functown; // Functors wrapping the backend initialisation functions
35 
36 QUICK_FUNCTION(DarkBit, TH_ProcessCatalog, OLD_CAPABILITY, TH_ProcessCatalog_WIMP, TH_ProcessCatalog, ())
38 QUICK_FUNCTION(DarkBit, DD_couplings, OLD_CAPABILITY, DD_couplings_WIMP, DM_nucleon_couplings, ())
39 
40 
41 void dump_array_to_file(const std::string & filename, const
42  boost::multi_array<double, 2> & a, const std::vector<double> & x, const
43  std::vector<double> & y)
44 {
45  std::fstream file;
46  file.open(filename, std::ios_base::out);
47  file << "0.0 ";
48  for (size_t i = 0; i < x.size(); i++)
49  file << x[i] << " ";
50  file << std::endl;
51  for (size_t j = 0; j < y.size(); j++)
52  {
53  file << y[j] << " ";
54  for (size_t i = 0; i < x.size(); i++)
55  {
56  file << a[i][j] << " ";
57  }
58  file << std::endl;
59  }
60  file.close();
61 }
62 
63 void dumpSpectrum(std::string filename, double mWIMP, double sv, std::vector<double> brList, double mPhi = -1)
64 {
65  DarkMatter_ID_WIMP.reset_and_calculate();
66  TH_ProcessCatalog_WIMP.setOption<std::vector<double>>("brList", brList);
67  TH_ProcessCatalog_WIMP.setOption<double>("mWIMP", mWIMP);
68  TH_ProcessCatalog_WIMP.setOption<double>("sv", sv);
69  if (mPhi != -1)
70  TH_ProcessCatalog_WIMP.setOption<double>("mPhi", mPhi);
71  TH_ProcessCatalog_WIMP.reset_and_calculate();
72  RD_fraction_one.reset_and_calculate();
73  SimYieldTable_DarkSUSY.reset_and_calculate();
74  SimYieldTable_MicrOmegas.reset_and_calculate();
75  GA_missingFinalStates.reset_and_calculate();
76  cascadeMC_FinalStates.reset_and_calculate();
77  cascadeMC_DecayTable.reset_and_calculate();
78  cascadeMC_LoopManager.reset_and_calculate();
79  cascadeMC_gammaSpectra.reset_and_calculate();
80  GA_AnnYield_General.reset_and_calculate();
81  dump_GammaSpectrum.setOption<std::string>("filename", filename);
82  dump_GammaSpectrum.reset_and_calculate();
83 }
84 
85 // ---- Set up basic internal structures for direct & indirect detection ----
86 
87 namespace Gambit
88 {
89  namespace DarkBit
90  {
91 
93  {
94  using namespace Pipes::TH_ProcessCatalog_WIMP;
95  using std::vector;
96  using std::string;
97 
98  // Initialize empty catalog and main annihilation process
99  TH_ProcessCatalog catalog;
100  TH_Process process_ann("WIMP", "WIMP");
101  TH_Process process_dec("phi");
102  TH_Process process_dec1("phi1");
103  TH_Process process_dec2("phi2");
104 
106  // Import particle masses and couplings
108 
109 #define addParticle(Name, Mass, spinX2) \
110  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
111  (Name , TH_ParticleProperty(Mass, spinX2)));
112 
114  double mWIMP = runOptions->getValue<double>("mWIMP");
116  double sv = runOptions->getValue<double>("sv");
117  double b = 0; // defined as sv(v) = sv(v=0) + b*(sv=0)*v**2
119  auto brList = runOptions->getValue<std::vector<double>>("brList");
121  double mPhi = runOptions->getValueOrDef<double>(59.0, "mPhi");
122 
123  addParticle("gamma", 0.0, 2)
124  addParticle("Z0", 91.2, 2)
125  addParticle("W+", 80.39, 2)
126  addParticle("W-", 80.39, 2)
127  addParticle("e+_3", 1.8, 1)
128  addParticle("e-_3", 1.8, 1)
129  addParticle("e+_1", 0.00051, 1)
130  addParticle("e-_1", 0.00051, 1)
131  addParticle("b", 4.9, 1)
132  addParticle("bbar", 4.9, 1)
133  addParticle("d_3", 4.9, 1)
134  addParticle("dbar_3", 4.9, 1)
135 
136  addParticle("WIMP", mWIMP, 0)
137  addParticle("phi", mPhi, 0)
138  addParticle("phi1", 100., 0)
139  addParticle("phi2", 100., 0)
140 #undef addParticle
141 
142  TH_Channel dec_channel(daFunk::vec<string>("gamma", "gamma"), daFunk::cnst(1.));
143  process_dec.channelList.push_back(dec_channel);
144 
145  TH_Channel dec_channel1(daFunk::vec<string>("e+_3", "e-_3"), daFunk::cnst(1.));
146  process_dec1.channelList.push_back(dec_channel1);
147 
148  TH_Channel dec_channel2(daFunk::vec<string>("d_3", "dbar_3"), daFunk::cnst(1.));
149  process_dec2.channelList.push_back(dec_channel2);
150 
151  process_ann.resonances_thresholds.threshold_energy.push_back(2*mWIMP);
152  auto p1 = daFunk::vec<string>("d_3", "gamma", "gamma", "e-_3", "W-", "e-_1", "phi");
153  auto p2 = daFunk::vec<string>("dbar_3", "Z0", "gamma", "e+_3", "W+", "e+_1", "phi2");
154  {
155  for ( unsigned int i = 0; i < brList.size()-1; i++ )
156  {
157  double mtot_final =
158  catalog.getParticleProperty(p1[i]).mass +
159  catalog.getParticleProperty(p2[i]).mass;
160  if ( mWIMP*2 > mtot_final && brList[i]!= 0.)
161  {
162  // std::cout << p1[i] << " " << p2[i] << " " << brList[i] << std::endl;
163  daFunk::Funk kinematicFunction = (daFunk::one("v")+pow(daFunk::var("v"), 2)*b)*sv*brList[i];
164  TH_Channel new_channel(
165  daFunk::vec<string>(p1[i], p2[i]), kinematicFunction
166  );
167  process_ann.channelList.push_back(new_channel);
168  }
169  else
170  {
172  push_back(mtot_final);
173  }
174  }
175  }
176 
177  if ( brList[7] > 0. )
178  {
179  auto E = daFunk::var("E");
180  // Note:: The below is an arbitrary form of the differential section for demonstration purposes
181  daFunk::Funk kinematicFunction = daFunk::one("v", "E1")/(pow(E-50, 4)+1)*sv*brList[7];
182  // Note: In the three body final states, the gamma yield from AnnYield currently is just the contribution
183  // from the first particle in the list (here the photon):
184  TH_Channel new_channel(daFunk::vec<string>("gamma", "e+_1", "e-_1"), kinematicFunction);
185  process_ann.channelList.push_back(new_channel);
186  }
187 
188  catalog.processList.push_back(process_ann);
189  catalog.processList.push_back(process_dec);
190  catalog.processList.push_back(process_dec1);
191  catalog.processList.push_back(process_dec2);
192 
193  catalog.validate();
194 
195  result = catalog;
196  } // function TH_ProcessCatalog_WIMP
197 
198  // Identifier for DM particle
199  void DarkMatter_ID_WIMP(std::string& result)
200  {
201  result = "WIMP";
202  }
203 
204  void DD_couplings_WIMP(DM_nucleon_couplings& result)
205  {
206  using namespace Pipes::DD_couplings_WIMP;
208  result.gps = runOptions->getValueOrDef<double>(0., "gps");
210  result.gns = runOptions->getValueOrDef<double>(0., "gns");
212  result.gpa = runOptions->getValueOrDef<double>(0., "gpa");
214  result.gna = runOptions->getValueOrDef<double>(0., "gna");
215  //std::cout << "DD_coupling says" << std::endl;
216  //std::cout << result.gps << std::endl;
217  }
218  }
219 }
220 
221 int main(int argc, char* argv[])
222 {
223  std::cout << std::endl;
224  std::cout << "Welcome to the DarkBit Generic WIMP standalone program!" << std::endl;
225  std::cout << std::endl;
226  std::cout << "**************************************************************************************" << std::endl;
227  std::cout << "This standalone example demonstrates how to calculate a range of observables and " << std::endl;
228  std::cout << "likelihoods for a generic WIMP model defined by the WIMP mass and an annihilation (or " << std::endl;
229  std::cout << "scattering) cross section. The model also contains three scalar particles which decay:" << std::endl;
230  std::cout << "phi -> gamma gamma phi_1 -> tau+ tau- phi_2 -> b bbar" << std::endl;
231  std::cout << std::endl;
232  std::cout << "Usage: DarkBit_standalone_WIMP mode" << std::endl;
233  std::cout << std::endl;
234  std::cout << "Mode Options: " << std::endl;
235  std::cout << " 0: Outputs spectrum of gamma rays from WIMP annihilation to b bbar (dPhi_dE0.dat)" << std::endl;
236  std::cout << " 1: Outputs spectrum of gamma rays from WIMP annihilation to gamma Z_0 (dPhi_dE1.dat)" << std::endl;
237  std::cout << " 2: Outputs spectrum of gamma rays from WIMP annihilation to gamma gamma (dPhi_dE2.dat)" << std::endl;
238  std::cout << " 3: Outputs spectrum of gamma rays from WIMP annihilation to tau+ tau- (dPhi_dE3.dat)" << std::endl;
239  std::cout << " 4: Outputs spectrum of gamma rays from WIMP annihilation to W+ W- (dPhi_dE4.dat)" << std::endl;
240  std::cout << " 5: Outputs spectrum of gamma rays from WIMP annihilation to gamma e+ e- " << std::endl;
241  std::cout << " (dPhi_dE5.dat)" << std::endl;
242  std::cout << " 6: Outputs tables of gamma-ray likelihoods and the relic density" << std::endl;
243  std::cout << " in <sigma v> / m_WIMP parameter space." << std::endl;
244  std::cout << " 7: Outputs tables of direct detection likelihoods in sigma / m_WIMP parameter" << std::endl;
245  std::cout << " space." << std::endl;
246  std::cout << " >=10: Outputs spectrum of gamma rays from WIMP annihilation to phi phi_2. The" << std::endl;
247  std::cout << " mode value is m_phi while m_phi_2=100 GeV (dPhi_dE_FCMC_(mode).dat)" << std::endl;
248  std::cout << " N.B. Here dPhi/dE = sigma v / m_chi^2 * dN/dE" << std::endl;
249  std::cout << "**************************************************************************************" << std::endl;
250  std::cout << std::endl;
251 
252  try
253  {
254  if (argc==1)
255  {
256  std::cout << "Please select test mode>=0" << std::endl;
257  exit(1);
258  }
259  int mode = std::stoi((std::string)argv[1]);
260  std::cout << "Starting with mode " << mode << std::endl;
261 
262 
263  // ---- Initialise logging and exceptions ----
264 
265  initialise_standalone_logs("runs/DarkBit_standalone_WIMP/logs/");
266  logger()<<"Running DarkBit standalone example"<<LogTags::info<<EOM;
267  model_warning().set_fatal(true);
268 
269 
270  // ---- Check that required backends are present ----
271 
272  if (not Backends::backendInfo().works["DarkSUSY_generic_wimp6.2.2"]) backend_error().raise(LOCAL_INFO, "DarkSUSY_generic_wimp_6.2.2 is missing!");
273  if (not Backends::backendInfo().works["gamLike1.0.1"]) backend_error().raise(LOCAL_INFO, "gamLike 1.0.1 is missing!");
274  if (not Backends::backendInfo().works["DDCalc2.2.0"]) backend_error().raise(LOCAL_INFO, "DDCalc 2.2.0 is missing!");
275  if (not Backends::backendInfo().works["MicrOmegas_MSSM3.6.9.2"]) backend_error().raise(LOCAL_INFO, "MicrOmegas 3.6.9.2 for MSSM is missing!");
276 
277  // ---- Initialize models ----
278 
279  // Initialize halo model
280  ModelParameters* Halo_primary_parameters = Models::Halo_Einasto::Functown::primary_parameters.getcontentsPtr();
281  Halo_primary_parameters->setValue("vrot", 235.); // Local properties
282  Halo_primary_parameters->setValue("v0", 235.);
283  Halo_primary_parameters->setValue("vesc", 550.);
284  Halo_primary_parameters->setValue("rho0", 0.4);
285  Halo_primary_parameters->setValue("r_sun", 8.5);
286 
287  Halo_primary_parameters->setValue("rs", 20.); // Global properties
288  Halo_primary_parameters->setValue("rhos", 0.08);
289  Halo_primary_parameters->setValue("alpha", 0.17);
290 
291 
292  // --- Resolve halo dependencies ---
293  ExtractLocalMaxwellianHalo.notifyOfModel("Halo_Einasto");
294  ExtractLocalMaxwellianHalo.resolveDependency(&Models::Halo_Einasto::Functown::primary_parameters);
295  ExtractLocalMaxwellianHalo.reset_and_calculate();
296 
297  GalacticHalo_Einasto.notifyOfModel("Halo_Einasto");
298  GalacticHalo_Einasto.resolveDependency(&Models::Halo_Einasto::Functown::primary_parameters);
299  GalacticHalo_Einasto.reset_and_calculate();
300 
301  // ---- Initialize backends ----
302 
303  // Assume for direct and indirect detection likelihoods that dark matter
304  // density is always the measured one (despite relic density results)
305  RD_fraction_one.reset_and_calculate();
306 
307  // Set up DDCalc backend initialization
308  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
309  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
310  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
311  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
312  // Assume for direct and indirect detection likelihoods that dark matter
313  // density is always the measured one (despite relic density results)
314  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
315  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
316  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_WIMP);
317 
318  // Initialize gamLike backend
319  gamLike_1_0_1_init.reset_and_calculate();
320 
321  // Initialize DarkSUSY backend
322  DarkSUSY_generic_wimp_6_2_2_init.reset_and_calculate();
323 
324  // Initialize MicrOmegas backend
325  // The below allows us to initialise MicrOmegas_MSSM without a particular MSSM model.
326  MicrOmegas_MSSM_3_6_9_2_init.notifyOfModel("Halo_Einasto");
327  MicrOmegas_MSSM_3_6_9_2_init.reset_and_calculate();
328 
329  // ---- Gamma-ray yields ----
330 
331  // Initialize tabulated gamma-ray yields
332  SimYieldTable_DarkSUSY.resolveBackendReq(&Backends::DarkSUSY_generic_wimp_6_2_2::Functown::dsanyield_sim);
333  SimYieldTable_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::dNdE);
334  SimYieldTable_DarkSUSY.setOption<bool>("allow_yield_extrapolation", true);
335  SimYieldTable_MicrOmegas.setOption<bool>("allow_yield_extrapolation", true);
336 
337  // Select SimYieldTable
338  //auto SimYieldTablePointer = &SimYieldTable_MicrOmegas;
339  auto SimYieldTablePointer = &SimYieldTable_DarkSUSY;
340 
341  // Collect missing final states for simulation in cascade MC
342  GA_missingFinalStates.resolveDependency(&TH_ProcessCatalog_WIMP);
343  GA_missingFinalStates.resolveDependency(SimYieldTablePointer);
344  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_WIMP);
345 
346  // Infer for which type of final states particles MC should be performed
347  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec((std::string)"gamma"));
348 
349  // Collect decay information for cascade MC
350  cascadeMC_DecayTable.resolveDependency(&TH_ProcessCatalog_WIMP);
351  cascadeMC_DecayTable.resolveDependency(SimYieldTablePointer);
352 
353  // Set up MC loop manager for cascade MC
354  cascadeMC_LoopManager.setOption<int>("cMC_maxEvents", 20000);
355  cascadeMC_Histograms.setOption<double>("cMC_endCheckFrequency", 25);
356  cascadeMC_Histograms.setOption<double>("cMC_gammaRelError", .05);
357  cascadeMC_Histograms.setOption<int>("cMC_numSpecSamples", 25);
358  cascadeMC_Histograms.setOption<int>("cMC_NhistBins", 300);
359  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
360  std::vector<functor*> nested_functions = initVector<functor*>(
362  cascadeMC_LoopManager.setNestedList(nested_functions);
363 
364  // Set up initial state for cascade MC step
365  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
366  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
367  //cascadeMC_InitialState.reset_and_calculate();
368 
369  // Perform MC step for cascade MC
370  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
371  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
372  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
373  //cascadeMC_GenerateChain.reset_and_calculate();
374 
375  // Generate histogram for cascade MC
376  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
377  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
378  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_WIMP);
379  cascadeMC_Histograms.resolveDependency(SimYieldTablePointer);
380  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
381  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
382  //cascadeMC_Histograms.reset_and_calculate();
383 
384  // Check convergence of cascade MC
385  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
386  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
387  //cascadeMC_EventCount.reset_and_calculate();
388 
389  // Start cascade MC loop
390 
391  // Infer gamma-ray spectra for recorded MC results
392  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
393  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
394  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
395  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
396 
397  // Calculate total gamma-ray yield (cascade MC + tabulated results)
398  GA_AnnYield_General.resolveDependency(&TH_ProcessCatalog_WIMP);
399  GA_AnnYield_General.resolveDependency(SimYieldTablePointer);
400  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_WIMP);
401  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
402 
403  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
404 
405  // Resolve Galactic halo requirements for gamLike
406  set_gamLike_GC_halo.resolveDependency(&GalacticHalo_Einasto);
407  set_gamLike_GC_halo.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::set_halo_profile);
408 
409  // Calculate Fermi LAT dwarf likelihood
411  // Assume for direct and indirect detection likelihoods that dark matter
412  // density is always the measured one (despite relic density results)
413  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
414  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
415 
416  lnL_HESSGC_gamLike.resolveDependency(&GA_AnnYield_General);
417  lnL_HESSGC_gamLike.resolveDependency(&RD_fraction_one);
418  lnL_HESSGC_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
419 
420  lnL_CTAGC_gamLike.resolveDependency(&GA_AnnYield_General);
421  lnL_CTAGC_gamLike.resolveDependency(&RD_fraction_one);
422  lnL_CTAGC_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
423 
424  lnL_FermiGC_gamLike.resolveDependency(&GA_AnnYield_General);
425  lnL_FermiGC_gamLike.resolveDependency(&RD_fraction_one);
426  lnL_FermiGC_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
427 
428 
429  // -- Calculate relic density --
430  // *any* of the models listed by "ALLOW_MODELS" in DarkBit_rollcall.hpp will work here
431  RD_eff_annrate_from_ProcessCatalog.notifyOfModel("ScalarSingletDM_Z2");
432  RD_eff_annrate_from_ProcessCatalog.resolveDependency(&TH_ProcessCatalog_WIMP);
433  RD_eff_annrate_from_ProcessCatalog.resolveDependency(&DarkMatter_ID_WIMP);
434 
437 
439 
440 
441  RD_oh2_DS_general.resolveDependency(&RD_spectrum_ordered_func);
442  RD_oh2_DS_general.resolveDependency(&RD_eff_annrate_from_ProcessCatalog);
443  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_generic_wimp_6_2_2::Functown::rdpars);
444  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_generic_wimp_6_2_2::Functown::rdtime);
445  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_generic_wimp_6_2_2::Functown::dsrdcom);
446  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_generic_wimp_6_2_2::Functown::dsrdstart);
447  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_generic_wimp_6_2_2::Functown::dsrdens);
448 
449 
450  // ---- Calculate direct detection constraints ----
451 
452  // Calculate direct detection rates for LZ, PandaX 2017, Xenon 1T and PICO-60
453  LZ_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
454  LZ_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
455  PandaX_2017_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
456  PandaX_2017_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
457  PICO_60_2017_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
458  PICO_60_2017_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
459  XENON1T_2017_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
460  XENON1T_2017_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
461 
462  // Calculate direct detection likelihood for LZ, PandaX 2017, Xenon 1T and PICO-60
463  LZ_GetLogLikelihood.resolveDependency(&LZ_Calc);
464  LZ_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
465  LZ_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
466  PandaX_2017_GetLogLikelihood.resolveDependency(&PandaX_2017_Calc);
467  PandaX_2017_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
468  PandaX_2017_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
469  XENON1T_2017_GetLogLikelihood.resolveDependency(&XENON1T_2017_Calc);
470  XENON1T_2017_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
471  XENON1T_2017_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
472  PICO_60_2017_GetLogLikelihood.resolveDependency(&PICO_60_2017_Calc);
473  PICO_60_2017_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
474  PICO_60_2017_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
475 
476  // Provide bin number in LZ
477  LZ_GetBinSignal.resolveDependency(&LZ_Calc);
478  LZ_GetBinSignal.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
479  LZ_GetBinSignal.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Bins);
480  LZ_GetBinSignal.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_BinSignal);
481 
482  // Set generic WIMP mass object
483  mwimp_generic.resolveDependency(&TH_ProcessCatalog_WIMP);
484  mwimp_generic.resolveDependency(&DarkMatter_ID_WIMP);
485  sigma_SI_p_simple.resolveDependency(&DD_couplings_WIMP);
486  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
487 
488  // Generate gamma-ray spectra for various final states
489  if ( (mode >= 0) and (mode < 6) )
490  {
491  std::cout << "Producing test spectra." << std::endl;
492  double mass = 100.;
493  double sv = 3e-26;
494  // The array that is being passed to dumpSpectrum give the branching fraction to various final states.
495  // They are (as defined in TH_ProcessCatalog_WIMP):
496  // 0: b bbar
497  // 1: gamma Z_0
498  // 2: gamma gamma
499  // 3: tau+ tau-
500  // 4: W+ W-
501  // 5: e+ e-
502  // 6: phi phi2
503  // 7: gamma e+ e-
504  if (mode==5) dumpSpectrum("dPhi_dE5.dat", mass, sv*0.1, daFunk::vec<double>(0., 0., 0., 0., 0., 0., 0., 1.));
505  if (mode==0) dumpSpectrum("dPhi_dE0.dat", mass, sv, daFunk::vec<double>(1., 0., 0., 0., 0., 0., 0., 0.));
506  if (mode==1) dumpSpectrum("dPhi_dE1.dat", mass, sv, daFunk::vec<double>(0., 1., 0., 0., 0., 0., 0., 0.));
507  if (mode==2) dumpSpectrum("dPhi_dE2.dat", mass, sv, daFunk::vec<double>(0., 0., 1., 0., 0., 0., 0., 0.));
508  if (mode==3) dumpSpectrum("dPhi_dE3.dat", mass, sv, daFunk::vec<double>(0., 0., 0., 1., 0., 0., 0., 0.));
509  if (mode==4) dumpSpectrum("dPhi_dE4.dat", mass, sv, daFunk::vec<double>(0., 0., 0., 0., 1., 0., 0., 0.));
510  }
511 
512  // Generate gamma-ray spectra for various masses
513  if (mode >= 10)
514  {
515  std::cout << "Producing test spectra." << std::endl;
516  double mass = 100.;
517  double sv = 3e-26;
518  std::string filename = "dPhi_dE_FCMC_" + std::to_string(mode) + ".dat";
519  dumpSpectrum(filename, mass, sv, daFunk::vec<double>(0., 0., 0., 0., 0., 0., 1., 0.), mode);
520  }
521 
522  // Systematic parameter maps annihilation
523  if (mode==6)
524  {
525  std::cout << "Producing gamma ray test maps." << std::endl;
526  int mBins = 60;
527  int svBins = 60;
528  double oh2, lnL;
529  std::vector<double> sv_list, m_list;
530 
531  GalacticHalo_Einasto.reset_and_calculate();
532  set_gamLike_GC_halo.reset_and_calculate();
533 
534  boost::multi_array<double, 2>
535  lnL_b_array{boost::extents[mBins][svBins]},
536  lnL_b_array2{boost::extents[mBins][svBins]},
537  lnL_b_array3{boost::extents[mBins][svBins]},
538  lnL_b_array4{boost::extents[mBins][svBins]},
539  lnL_tau_array{boost::extents[mBins][svBins]};
540  boost::multi_array<double, 2> oh2_array{boost::extents[mBins][svBins]};
541 
542  sv_list = daFunk::logspace(-28.0, -22.0, svBins);
543 
544  std::cout << "Calculating gamma-ray likelihood tables for annihilation to b bbar." << std::endl;
545  m_list = daFunk::logspace(log10(5.), 4., mBins);
546  for (size_t i = 0; i < m_list.size(); i++)
547  {
548  for (size_t j = 0; j < sv_list.size(); j++)
549  {
550  TH_ProcessCatalog_WIMP.setOption<double>("mWIMP", m_list[i]);
551  TH_ProcessCatalog_WIMP.setOption<double>("sv", sv_list[j]);
552  //std::cout << "Parameters: " << m_list[i] << " " << sv_list[j] << std::endl;
553 
554  TH_ProcessCatalog_WIMP.setOption<std::vector<double>>("brList", daFunk::vec<double>(1., 0., 0., 0., 0., 0., 0., 0.));
555  DarkMatter_ID_WIMP.reset_and_calculate();
556  TH_ProcessCatalog_WIMP.reset_and_calculate();
557  RD_fraction_one.reset_and_calculate();
558  SimYieldTable_DarkSUSY.reset_and_calculate();
559  SimYieldTable_MicrOmegas.reset_and_calculate();
560  GA_missingFinalStates.reset_and_calculate();
561  cascadeMC_FinalStates.reset_and_calculate();
562  cascadeMC_DecayTable.reset_and_calculate();
563  cascadeMC_LoopManager.reset_and_calculate();
564  cascadeMC_gammaSpectra.reset_and_calculate();
565  GA_AnnYield_General.reset_and_calculate();
566  lnL_FermiLATdwarfs_gamLike.setOption<std::string>("version", "pass8");
567  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
569  //std::cout << "Fermi dwarf likelihood: " << lnL << std::endl;
570  lnL_b_array[i][j] = lnL;
571  lnL_HESSGC_gamLike.setOption<std::string>("version", "integral_fixedJ");
572  lnL_HESSGC_gamLike.reset_and_calculate();
573  lnL = lnL_HESSGC_gamLike(0);
574  //std::cout << "HESS GC likelihood: " << lnL << std::endl;
575  lnL_b_array2[i][j] = lnL;
576  lnL_CTAGC_gamLike.reset_and_calculate();
577  lnL = lnL_CTAGC_gamLike(0);
578  //std::cout << "CTA GC likelihood: " << lnL << std::endl;
579  lnL_b_array3[i][j] = lnL;
580  lnL_FermiGC_gamLike.setOption<std::string>("version", "fixedJ");
581  lnL_FermiGC_gamLike.reset_and_calculate();
582  lnL = lnL_FermiGC_gamLike(0);
583  lnL_b_array4[i][j] = lnL;
584  //std::cout << "Fermi GC likelihood: " << lnL << std::endl;
585  }
586  }
587 
588  dump_array_to_file("FermiD_b_table.dat", lnL_b_array, m_list, sv_list);
589  dump_array_to_file("HESSGC_b_table.dat", lnL_b_array2, m_list, sv_list);
590  dump_array_to_file("CTAGC_b_table.dat", lnL_b_array3, m_list, sv_list);
591  dump_array_to_file("FermiGC_b_table.dat", lnL_b_array4, m_list, sv_list);
592 
593  std::cout << "Calculating Fermi-LAT dwarf spheroidal likehood table for annihilation to tau+ tau-." << std::endl;
594  m_list = daFunk::logspace(log10(1.9), 4., mBins);
595  for (size_t i = 0; i < m_list.size(); i++)
596  {
597  for (size_t j = 0; j < sv_list.size(); j++)
598  {
599  TH_ProcessCatalog_WIMP.setOption<double>("mWIMP", m_list[i]);
600  TH_ProcessCatalog_WIMP.setOption<double>("sv", sv_list[j]);
601  //std::cout << "Parameters: " << m_list[i] << " " << sv_list[j] << std::endl;
602 
603  TH_ProcessCatalog_WIMP.setOption<std::vector<double>>("brList", daFunk::vec<double>(0., 0., 0., 1., 0., 0., 0., 0.));
604  DarkMatter_ID_WIMP.reset_and_calculate();
605  TH_ProcessCatalog_WIMP.reset_and_calculate();
606  RD_fraction_one.reset_and_calculate();
607  SimYieldTable_DarkSUSY.reset_and_calculate();
608  SimYieldTable_MicrOmegas.reset_and_calculate();
609  GA_missingFinalStates.reset_and_calculate();
610  cascadeMC_FinalStates.reset_and_calculate();
611  cascadeMC_DecayTable.reset_and_calculate();
612  cascadeMC_LoopManager.reset_and_calculate();
613  cascadeMC_gammaSpectra.reset_and_calculate();
614  GA_AnnYield_General.reset_and_calculate();
615  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
617  //std::cout << "Fermi LAT likelihood: " << lnL << std::endl;
618  lnL_tau_array[i][j] = lnL;
619  }
620  }
621 
622  dump_array_to_file("FermiD_tau_table.dat", lnL_tau_array, m_list, sv_list);
623 
624  std::cout << "Calculating table of Omega h^2 values." << std::endl;
625  m_list = daFunk::logspace(-1.0, 4., mBins);
626  for (size_t i = 0; i < m_list.size(); i++)
627  {
628  for (size_t j = 0; j < sv_list.size(); j++)
629  {
630  TH_ProcessCatalog_WIMP.setOption<double>("mWIMP", m_list[i]);
631  TH_ProcessCatalog_WIMP.setOption<double>("sv", sv_list[j]);
632  //std::cout << "Parameters: " << m_list[i] << " " << sv_list[j] << std::endl;
633 
634  TH_ProcessCatalog_WIMP.setOption<std::vector<double>>("brList", daFunk::vec<double>(0., 0., 0., 0., 0., 1., 0., 0.));
635  DarkMatter_ID_WIMP.reset_and_calculate();
636  TH_ProcessCatalog_WIMP.reset_and_calculate();
637  RD_eff_annrate_from_ProcessCatalog.reset_and_calculate();
638  RD_spectrum_from_ProcessCatalog.reset_and_calculate();
639  RD_spectrum_ordered_func.reset_and_calculate();
640  RD_oh2_DS_general.reset_and_calculate();
641  oh2 = RD_oh2_DS_general(0);
642  //std::cout << "Omega h^2 = " << oh2 << std::endl;
643  oh2_array[i][j] = oh2;
644  }
645  }
646 
647  dump_array_to_file("oh2_table.dat", oh2_array, m_list, sv_list);
648  }
649 
650  // Systematic parameter maps scattering
651  if (mode==7)
652  {
653  std::cout << "Producing direct detection test maps." << std::endl;
654  double lnL1, lnL2, lnL3, lnL4;
655  int nbins;
656  double g, reduced_mass;
657  //int mBins = 300;
658  //int sBins = 200;
659  int mBins = 120;
660  int sBins = 80;
661  const double mN = (m_proton + m_neutron) / 2;
662  std::vector<double> m_list = daFunk::logspace(0.0, 4.0, mBins);
663  std::vector<double> s_list;
664  boost::multi_array<double, 2> lnL_array1{boost::extents[mBins][sBins]},
665  lnL_array2{boost::extents[mBins][sBins]}, lnL_array3{boost::extents[mBins][sBins]},
666  lnL_array4{boost::extents[mBins][sBins]};
667  TH_ProcessCatalog_WIMP.setOption<double>("sv", 0.);
668  TH_ProcessCatalog_WIMP.setOption<std::vector<double>>("brList", daFunk::vec<double>(1., 0., 0., 0., 0., 0., 0., 0.));
669 
670  s_list = daFunk::logspace(-47., -40., sBins);
671  // Calculate array of sigma_SI and lnL values for LZ, PandaX, XENON1T and PICO-60
672  // assuming gps=gns
673 
674  std::cout << "Calculating tables of SI likelihoods." << std::endl;
675  for (size_t i = 0; i < m_list.size(); i++)
676  {
677  for (size_t j = 0; j < s_list.size(); j++)
678  {
679  // Re-initialize DDCalc with LZ/Xenon/PandaX halo parameters
680  Halo_primary_parameters->setValue("rho0", 0.3);
681  Halo_primary_parameters->setValue("vrot", 232.7); // v_Earth = 245 km/s
682  Halo_primary_parameters->setValue("v0", 220.);
683  Halo_primary_parameters->setValue("vesc", 544.);
684  ExtractLocalMaxwellianHalo.reset_and_calculate();
685 
686  TH_ProcessCatalog_WIMP.setOption<double>("mWIMP", m_list[i]);
687  //std::cout << "Parameters: " << m_list[i] << " " << s_list[j] << std::endl;
688  reduced_mass = (m_list[i] * mN) / (mN + m_list[i]);
689  g = sqrt(s_list[j]*pi/gev2cm2) / (reduced_mass);
690  DarkMatter_ID_WIMP.reset_and_calculate();
691  TH_ProcessCatalog_WIMP.reset_and_calculate();
692  // Assume for direct and indirect detection likelihoods that dark matter
693  // density is always the measured one (despite relic density results)
694  RD_fraction_one.reset_and_calculate();
695  DD_couplings_WIMP.setOption<double>("gps", g);
696  DD_couplings_WIMP.setOption<double>("gns", g);
697  DD_couplings_WIMP.setOption<double>("gpa", 0.);
698  DD_couplings_WIMP.setOption<double>("gna", 0.);
699  DD_couplings_WIMP.reset_and_calculate();
700  mwimp_generic.reset_and_calculate();
701 
702  DDCalc_2_2_0_init.reset_and_calculate();
703  LZ_Calc.reset_and_calculate();
704  LZ_GetLogLikelihood.reset_and_calculate();
705 
706  XENON1T_2017_Calc.reset_and_calculate();
707  XENON1T_2017_GetLogLikelihood.reset_and_calculate();
708  PandaX_2017_Calc.reset_and_calculate();
709  PandaX_2017_GetLogLikelihood.reset_and_calculate();
710 
711  lnL1 = LZ_GetLogLikelihood(0);
712  lnL2 = PandaX_2017_GetLogLikelihood(0);
713  lnL3 = XENON1T_2017_GetLogLikelihood(0);
714 
715  // Set LocalHalo Model parameters to PICO-60 values
716  Halo_primary_parameters->setValue("rho0", 0.3);
717  Halo_primary_parameters->setValue("vrot", 220.);
718  Halo_primary_parameters->setValue("v0", 220.);
719  Halo_primary_parameters->setValue("vesc", 544.);
720  ExtractLocalMaxwellianHalo.reset_and_calculate();
721 
722  DDCalc_2_2_0_init.reset_and_calculate();
723  PICO_60_2017_Calc.reset_and_calculate();
724  PICO_60_2017_GetLogLikelihood.reset_and_calculate();
725  lnL4 = PICO_60_2017_GetLogLikelihood(0);
726 
727  //std::cout << "LZ SI lnL = " << lnL1 << std::endl;
728  //std::cout << "PandaX_2017 SI lnL = " << lnL2 << std::endl;
729  //std::cout << "XENON1T_2017 SI lnL = " << lnL3 << std::endl;
730  //std::cout << "PICO_60_2017 SI lnL = " << lnL4 << std::endl;
731 
732  DDCalc_2_2_0_init.reset_and_calculate();
733  LZ_Calc.reset_and_calculate();
734  std::vector<double> events;
735  LZ_GetBinSignal.reset_and_calculate();
736  events = LZ_GetBinSignal(0);
737  nbins = events.size();
738  std::cout << "Number of LZ bins: " << nbins << std::endl;
739  std::cout << "Predicted signal: ";
740  for (int ibin=0;ibin<=nbins-1;ibin++) {
741  std::cout << events[ibin] << " ";
742  }
743  std::cout << std::endl;
744 
745  lnL_array1[i][j] = lnL1;
746  lnL_array2[i][j] = lnL2;
747  lnL_array3[i][j] = lnL3;
748  lnL_array4[i][j] = lnL4;
749  }
750  }
751 
752  dump_array_to_file("LZ_SI_table.dat", lnL_array1, m_list, s_list);
753  dump_array_to_file("PandaX_2017_SI_table.dat", lnL_array2, m_list, s_list);
754  dump_array_to_file("XENON1T_2017_SI_table.dat", lnL_array3, m_list, s_list);
755  dump_array_to_file("PICO_60_2017_SI_table.dat", lnL_array4, m_list, s_list);
756 
757  s_list = daFunk::logspace(-42., -35., sBins);
758  // Calculate array of sigma_SI and lnL values for LZ, PandaX, XENON1T and PICO-60
759  // assuming gna=0 (proton-only)
760 
761  std::cout << "Calculating tables of SD likelihoods." << std::endl;
762  for (size_t i = 0; i < m_list.size(); i++)
763  {
764  for (size_t j = 0; j < s_list.size(); j++)
765  {
766  // Re-initialize DDCalc with LZ/Xenon/PandaX halo parameters
767  Halo_primary_parameters->setValue("rho0", 0.3);
768  Halo_primary_parameters->setValue("vrot", 232.7); // v_Earth = 245 km/s
769  Halo_primary_parameters->setValue("v0", 220.);
770  Halo_primary_parameters->setValue("vesc", 544.);
771  ExtractLocalMaxwellianHalo.reset_and_calculate();
772 
773  TH_ProcessCatalog_WIMP.setOption<double>("mWIMP", m_list[i]);
774  //std::cout << "Parameters: " << m_list[i] << " " << s_list[j] << std::endl;
775  reduced_mass = (m_list[i] * m_proton) / (m_proton + m_list[i]);
776  g = sqrt(s_list[j]*pi/(3*gev2cm2)) / (reduced_mass);
777  DarkMatter_ID_WIMP.reset_and_calculate();
778  TH_ProcessCatalog_WIMP.reset_and_calculate();
779  RD_fraction_one.reset_and_calculate();
780  DD_couplings_WIMP.setOption<double>("gps", 0.);
781  DD_couplings_WIMP.setOption<double>("gns", 0.);
782  DD_couplings_WIMP.setOption<double>("gpa", g);
783  DD_couplings_WIMP.setOption<double>("gna", 0.);
784  DD_couplings_WIMP.reset_and_calculate();
785  mwimp_generic.reset_and_calculate();
786 
787  DDCalc_2_2_0_init.reset_and_calculate();
788  LZ_Calc.reset_and_calculate();
789  LZ_GetLogLikelihood.reset_and_calculate();
790  XENON1T_2017_Calc.reset_and_calculate();
791  XENON1T_2017_GetLogLikelihood.reset_and_calculate();
792  PandaX_2017_Calc.reset_and_calculate();
793  PandaX_2017_GetLogLikelihood.reset_and_calculate();
794  lnL1 = LZ_GetLogLikelihood(0);
795  lnL2 = PandaX_2017_GetLogLikelihood(0);
796  lnL3 = XENON1T_2017_GetLogLikelihood(0);
797 
798  // Set LocalHalo Model parameters to PICO-60 values
799  Halo_primary_parameters->setValue("rho0", 0.3);
800  Halo_primary_parameters->setValue("vrot", 220.);
801  Halo_primary_parameters->setValue("v0", 220.);
802  Halo_primary_parameters->setValue("vesc", 544.);
803  ExtractLocalMaxwellianHalo.reset_and_calculate();
804 
805  DDCalc_2_2_0_init.reset_and_calculate();
806  PICO_60_2017_Calc.reset_and_calculate();
807  PICO_60_2017_GetLogLikelihood.reset_and_calculate();
808  lnL4 = PICO_60_2017_GetLogLikelihood(0);
809 
810  //std::cout << "LZ SD lnL = " << lnL1 << std::endl;
811  //std::cout << "PandaX_2017 SD lnL = " << lnL2 << std::endl;
812  //std::cout << "XENON1T_2017 SD lnL = " << lnL3 << std::endl;
813  //std::cout << "PICO_60_2017 SD lnL = " << lnL4 << std::endl;
814 
815  lnL_array1[i][j] = lnL1;
816  lnL_array2[i][j] = lnL2;
817  lnL_array3[i][j] = lnL3;
818  lnL_array4[i][j] = lnL4;
819  }
820  }
821 
822  dump_array_to_file("LZ_SD_table.dat", lnL_array1, m_list, s_list);
823  dump_array_to_file("PandaX_2017_SD_table.dat", lnL_array2, m_list, s_list);
824  dump_array_to_file("XENON1T_2017_SD_table.dat", lnL_array3, m_list, s_list);
825  dump_array_to_file("PICO_60_2017_SD_table.dat", lnL_array4, m_list, s_list);
826 
827  // Reset halo parameters to DarkBit defaults.
828  Halo_primary_parameters->setValue("rho0", 0.4);
829  Halo_primary_parameters->setValue("vrot", 235.);
830  Halo_primary_parameters->setValue("v0", 235.);
831  Halo_primary_parameters->setValue("vesc", 550.);
832  ExtractLocalMaxwellianHalo.reset_and_calculate();
833  GalacticHalo_Einasto.reset_and_calculate();
834 
835  }
836  }
837 
838  catch (std::exception& e)
839  {
840  std::cout << "DarkBit_standalone_WIMP has exited with fatal exception: " << e.what() << std::endl;
841  }
842 
843  return 0;
844 }
void lnL_FermiGC_gamLike(double &result)
Fermi LAT galactic center likelihoods, using gamLike backend.
std::vector< TH_Channel > channelList
List of channels.
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 dumpSpectrum(std::string filename, double mWIMP, double sv, std::vector< double > brList, double mPhi=-1)
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.
const double m_proton
void TH_ProcessCatalog_WIMP(TH_ProcessCatalog &result)
void sigma_SI_p_simple(double &result)
Simple calculator of the spin-independent WIMP-proton cross-section.
Includes everything needed to use a GAMBIT module as a standalone analysis code.
const double m_neutron
void RD_fraction_one(double &result)
#define LOCAL_INFO
Definition: local_info.hpp:34
void SimYieldTable_DarkSUSY(SimYieldTable &result)
SimYieldTable based on DarkSUSY6 tabulated results.
Definition: GamYields.cpp:558
STL namespace.
void cascadeMC_InitialState(std::string &pID)
Function selecting initial state for decay chain.
Definition: Cascades.cpp:140
void DD_couplings_WIMP(DM_nucleon_couplings &result)
void dump_GammaSpectrum(double &result)
Helper function to dump gamma-ray spectra.
void SimYieldTable_MicrOmegas(SimYieldTable &result)
SimYieldTable based on MicrOmegas tabulated results.
Definition: GamYields.cpp:691
Helpers for making simple spectra from SLHAea objects.
Funk var(std::string arg)
Definition: daFunk.hpp:921
START_MODEL b
Definition: demo.hpp:235
warning & model_warning()
Model warnings.
void RD_spectrum_from_ProcessCatalog(RD_spectrum_type &result)
Collects information about resonances and threshold energies directly from the ProcessCatalog [NB: th...
General small utility functions.
error & backend_error()
Backend errors.
double mass
Particle mass (GeV)
void cascadeMC_EventCount(std::map< std::string, int > &counts)
Event counter for cascade decays.
Definition: Cascades.cpp:175
void setValue(std::string const &inkey, double const &value)
Set single parameter value.
int main(int argc, char *argv[])
void cascadeMC_GenerateChain(DarkBit::DecayChain::ChainContainer &chain)
Function for generating decay chains.
Definition: Cascades.cpp:198
const double pi
void GA_AnnYield_General(daFunk::Funk &result)
General routine to derive annihilation yield.
Definition: GamYields.cpp:191
TH_resonances_thresholds resonances_thresholds
List of resonances and thresholds.
void initialise_standalone_logs(str)
Logger setup standalone utility function.
void lnL_FermiLATdwarfs_gamLike(double &result)
Fermi LAT dwarf likelihoods, using gamLike backend.
A container holding all annihilation and decay initial states relevant for DarkBit.
std::vector< double > logspace(double x0, double x1, unsigned int n)
Definition: daFunk.hpp:186
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
void set_gamLike_GC_halo(bool &result)
Fermi LAT dwarf likelihoods, based on arXiv:1108.2914.
QUICK_FUNCTION(DarkBit, TH_ProcessCatalog, OLD_CAPABILITY, TH_ProcessCatalog_WIMP, TH_ProcessCatalog,()) QUICK_FUNCTION(DarkBit
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
All data on a single annihilation or decay channel, e.g.
const double gev2cm2
void ExtractLocalMaxwellianHalo(LocalMaxwellianHalo &result)
Module function providing local density and velocity dispersion parameters.
Definition: DarkBit.cpp:152
Rollcall header for module DarkBit.
void lnL_HESSGC_gamLike(double &result)
Funk cnst(double x, Args... argss)
Definition: daFunk.hpp:606
double pow(const double &a)
Outputs a^i.
void mwimp_generic(double &result)
Retrieve the DM mass in GeV for generic models (GeV)
Definition: DarkBit.cpp:60
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
A container for a single process.
std::vector< T > vec(std::vector< T > vector)
Definition: daFunk.hpp:142
void RD_oh2_DS_general(double &result)
General routine for calculation of relic density, using DarkSUSY 6+ Boltzmann solver.
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
std::vector< TH_Process > processList
Vector of all processes in this catalog.
void GalacticHalo_Einasto(GalacticHaloProperties &result)
Module function to generate GalacticHaloProperties for Einasto profile.
Definition: DarkBit.cpp:140
void lnL_CTAGC_gamLike(double &result)
#define addParticle(Name, Mass, spinX2)
void RD_spectrum_ordered_func(RD_spectrum_type &result)
Order RD_spectrum object and derive coannihilation thresholds.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
double E(const double x, const double y)
Funk one(Args... argss)
Definition: daFunk.hpp:602