gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
DarkBit_standalone_MSSM.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
34 
41 
42 using namespace DarkBit::Functown; // Functors wrapping the module's actual module functions
43 using namespace BackendIniBit::Functown; // Functors wrapping the backend initialisation functions
44 
47 QUICK_FUNCTION(DarkBit, SLHA_pseudonyms, NEW_CAPABILITY, createSLHA1Names, mass_es_pseudonyms, (), (MSSM_spectrum, Spectrum))
48 QUICK_FUNCTION(DarkBit, cascadeMC_gammaSpectra, OLD_CAPABILITY, CMC_dummy, DarkBit::stringFunkMap, ())
49 
50 
51 namespace Gambit
52 {
53  namespace DarkBit
54  {
55  // Dummy functor for returning empty cascadeMC result spectra
57  {
59  result = sfm;
60  }
61 
62  // Create spectrum object from SLHA file
63  void createSpectrum(Spectrum& outSpec)
64  {
65  using namespace Pipes::createSpectrum;
67  std::string inputFileName = runOptions->getValue<std::string>("filename");
68  std::cout << "Loading spectrum from: " << inputFileName << std::endl;
69  outSpec = spectrum_from_SLHA<MSSMSimpleSpec>(inputFileName, Spectrum::mc_info(), Spectrum::mr_info());
70  }
71 
72  // Create decay object from SLHA file
73  void createDecays(DecayTable& outDecays)
74  {
75  using namespace Pipes::createDecays;
77  std::string inputFileName = runOptions->getValue<std::string>("filename");
78  std::cout << "Loading decays from: " << inputFileName << std::endl;
79  outDecays = DecayTable(inputFileName, 0, true);
80  }
81 
82  // Create SLHA1 pseudonyms from Spectrum object
83  void createSLHA1Names(mass_es_pseudonyms& names)
84  {
85  const double gauge_mixing_tol = 0.5;
86  const bool tol_invalidates_pt = true;
87  const bool debug = false;
88  names.refill(Pipes::createSLHA1Names::Dep::MSSM_spectrum->get_HE(), gauge_mixing_tol, tol_invalidates_pt, debug);
89  }
90 
91  }
92 }
93 
94 int main(int argc, char* argv[])
95 {
96  std::cout << std::endl;
97  std::cout << "Welcome to the DarkBit MSSM standalone program!" << std::endl;
98  std::cout << std::endl;
99  std::cout << "********************************************************************************" << std::endl;
100  std::cout << "Usage: DarkBit_standalone_MSSM SLHA_file (spectrum) (output)" << std::endl;
101  std::cout << std::endl;
102  std::cout << "SLHA_file: SLHA file used to intialise the program (required)" << std::endl;
103  std::cout << "(spectrum): name of output file for gamma-ray spectrum (default: BACKENDNAME_dNdE.dat)" << std::endl;
104  std::cout << "(output): name of output file for observables and likelihoods (default: DarkBit_standalone_MSSM.out)" << std::endl;
105  std::cout << std::endl;
106  std::cout << "The SLHA files for the MSSM-7 benchmarks in the DarkBit paper are located in " << std::endl;
107  std::cout << "DarkBit/data/benchmarks/ " << std::endl;
108  std::cout << "********************************************************************************" << std::endl;
109  std::cout << std::endl;
110 
111  try
112  {
113  if (argc == 1)
114  {
115  std::cout << "Please provide name of SLHA file at command line." << std::endl;
116  exit(1);
117  }
118  std::string filename = argv[1];
119  std::string outname_dNdE_spectrum = "dNdE.dat";
120  if (argc >= 3) outname_dNdE_spectrum = argv[2];
121  std::string outname_data = "DarkBit_standalone_MSSM.out";
122  if (argc >= 4) outname_data = argv[3];
123 
124 
125  // ---- Initialise logging and exceptions ----
126 
127  initialise_standalone_logs("runs/DarkBit_standalone_MSSM/logs/");
128  logger()<<"Running DarkBit standalone example"<<LogTags::info<<EOM;
129  model_warning().set_fatal(true);
130 
131  // ---- Check which backends are present ----
132  if (not Backends::backendInfo().works["gamLike1.0.1"]) backend_error().raise(LOCAL_INFO, "gamLike 1.0.1 is missing!");
133  if (not Backends::backendInfo().works["DDCalc2.2.0"]) backend_error().raise(LOCAL_INFO, "DDCalc 2.2.0 is missing!");
134  if (not Backends::backendInfo().works["nulike1.0.9"]) backend_error().raise(LOCAL_INFO, "nulike 1.0.9 is missing!");
135 
136  // ---- Useful variables ----
137  //
138  // Prepare a str-str-double map of maps to hold the results.
139  // We will add results for the individual backends as
140  // results["oh2"]["MicrOmegas_MSSM3.6.9.2"] = ...
141  // results["oh2"]["DarkSUSY_MSSM6.1.1"] = ...
142  //
143  std::map<std::string, std::map<std::string,double> > results;
144  results["oh2"] = std::map<std::string,double>();
145  results["oh2_lnL"] = std::map<std::string,double>();
146  results["sigma_SI_p"] = std::map<std::string,double>();
147  results["sigma_SD_p"] = std::map<std::string,double>();
148  results["LUX_2016_lnL"] = std::map<std::string,double>();
149  results["IceCube_79_lnL"] = std::map<std::string,double>();
150  results["sigmav0"] = std::map<std::string,double>();
151  results["FermiLAT_dwarfsph_lnL"] = std::map<std::string,double>();
152 
153  std::map<std::string, std::string> results_units;
154  results_units["oh2"] = "";
155  results_units["oh2_lnL"] = "";
156  results_units["sigma_SI_p"] = "cm^2";
157  results_units["sigma_SD_p"] = "cm^2";
158  results_units["LUX_2016_lnL"] = "";
159  results_units["IceCube_79_lnL"] = "";
160  results_units["sigmav0"] = "cm^3/s";
161  results_units["FermiLAT_dwarfsph_lnL"] = "";
162 
163  std::vector<std::string> result_output_order = {
164  "oh2",
165  "oh2_lnL",
166  "sigma_SI_p",
167  "sigma_SD_p",
168  "LUX_2016_lnL",
169  "IceCube_79_lnL",
170  "sigmav0",
171  "FermiLAT_dwarfsph_lnL",
172  };
173 
174 
175  // A string to refer to the current backend being used for calculations
176  std::string current_backend = "";
177 
178  // Keep track of which backends are not installed
179  std::vector<std::string> backends_not_built;
180 
181  // ---- Initialize models ----
182 
183  // Initialize halo model
184  ModelParameters* Halo_primary_parameters = Models::Halo_Einasto::Functown::primary_parameters.getcontentsPtr();
185  Halo_primary_parameters->setValue("rho0", 0.4);
186  Halo_primary_parameters->setValue("rhos", 0.08);
187  Halo_primary_parameters->setValue("vrot", 235.);
188  Halo_primary_parameters->setValue("v0", 235.);
189  Halo_primary_parameters->setValue("vesc", 550.);
190  Halo_primary_parameters->setValue("rs", 20.);
191  Halo_primary_parameters->setValue("r_sun", 8.5);
192  Halo_primary_parameters->setValue("alpha", 0.17);
193 
194  // --- Resolve halo dependencies ---
195  ExtractLocalMaxwellianHalo.notifyOfModel("Halo_Einasto");
196  ExtractLocalMaxwellianHalo.resolveDependency(&Models::Halo_Einasto::Functown::primary_parameters);
197  ExtractLocalMaxwellianHalo.reset_and_calculate();
198 
199  GalacticHalo_Einasto.notifyOfModel("Halo_Einasto");
200  GalacticHalo_Einasto.resolveDependency(&Models::Halo_Einasto::Functown::primary_parameters);
201  GalacticHalo_Einasto.reset_and_calculate();
202 
203  // Initialize nuclear_params_fnq model
204  ModelParameters* nuclear_params_fnq = Models::nuclear_params_fnq::Functown::primary_parameters.getcontentsPtr();
205  nuclear_params_fnq->setValue("fpd", 0.034);
206  nuclear_params_fnq->setValue("fpu", 0.023);
207  nuclear_params_fnq->setValue("fps", 0.14);
208  nuclear_params_fnq->setValue("fnd", 0.041);
209  nuclear_params_fnq->setValue("fnu", 0.019);
210  nuclear_params_fnq->setValue("fns", 0.14);
211  nuclear_params_fnq->setValue("deltad", -0.40);
212  nuclear_params_fnq->setValue("deltau", 0.74);
213  nuclear_params_fnq->setValue("deltas", -0.12);
214 
215  // Resolve other dependencies related directly to the GAMBIT Models
216  DD_couplings_DarkSUSY_DS5.notifyOfModel("nuclear_params_fnq");
217  DD_couplings_DarkSUSY_DS5.resolveDependency(&Models::nuclear_params_fnq::Functown::primary_parameters);
218 
219  DD_couplings_DarkSUSY_MSSM.notifyOfModel("nuclear_params_fnq");
220  DD_couplings_DarkSUSY_MSSM.resolveDependency(&Models::nuclear_params_fnq::Functown::primary_parameters);
221 
222  DD_couplings_MicrOmegas.notifyOfModel("MSSM30atQ");
223  DD_couplings_MicrOmegas.notifyOfModel("nuclear_params_fnq");
224  DD_couplings_MicrOmegas.resolveDependency(&Models::nuclear_params_fnq::Functown::primary_parameters);
225 
226 
227  // ---- Initialize spectrum and decays from SLHA file ----
228  createSpectrum.setOption<std::string>("filename", filename);
229  createSpectrum.reset_and_calculate();
230  createSLHA1Names.resolveDependency(&createSpectrum);
231  createSLHA1Names.reset_and_calculate();
232  createDecays.setOption<std::string>("filename", filename);
233  createDecays.reset_and_calculate();
234 
235  // Check that the decay table contains ~chi0_2 (if it doesn't,
236  // we do not use information from the SLHA decay block)
237  bool decays = true;
238  try { createDecays(0).at("~chi0_2"); }
239  catch(std::exception& e)
240  {
241  decays = false;
242  cout << "It appears that the decay block is missing from the SLHA file. Decay widths\n"
243  "will be determined by the backends." << endl;
244  }
245 
246  // Set identifier for DM particle
247  DarkMatter_ID_MSSM.resolveDependency(&createSpectrum);
248  DarkMatter_ID_MSSM.reset_and_calculate();
249 
250  // Assume for direct and indirect detection likelihoods that dark matter
251  // density is always the measured one (regardless of relic density results)
252  RD_fraction_one.reset_and_calculate();
253 
254  //
255  // ======= Initializations that can be done once =======
256  //
257 
258  // Initialize nulike backend
259  Backends::nulike_1_0_9::Functown::nulike_bounds.setStatus(2);
260  nulike_1_0_9_init.reset_and_calculate();
261 
262  // Initialize gamLike backend
263  gamLike_1_0_1_init.reset_and_calculate();
264 
265 
266 
267  //
268  // ======= Perform all calculations for backend DarkSUSY 5.1.3 =======
269  //
270 
271  current_backend = "DarkSUSY5.1.3";
272 
273  if (not Backends::backendInfo().works[current_backend])
274  {
275  backends_not_built.push_back(current_backend);
276  }
277  else
278  {
279  // Initialize DarkSUSY 5 backend
280  DarkSUSY_5_1_3_init.notifyOfModel("MSSM30atQ");
281  DarkSUSY_5_1_3_init.resolveDependency(&createSpectrum);
282  DarkSUSY_5_1_3_init.resolveDependency(&createDecays);
283  if (decays) DarkSUSY_5_1_3_init.setOption<bool>("use_dsSLHAread", false);
284  else DarkSUSY_5_1_3_init.setOption<bool>("use_dsSLHAread", true);
285  DarkSUSY_5_1_3_init.reset_and_calculate();
286 
287  // Initialize DarkSUSY 5 Local Halo Model
290  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmcom);
291  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmisodf);
292  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmframevelcom);
293  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmnoclue);
294  DarkSUSY5_PointInit_LocalHalo_func.reset_and_calculate();
295 
296 
297  // Relic density calculation with GAMBIT routines and DarkSUSY 5:
298  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::mspctm);
299  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::widths);
300  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::intdof);
301  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::pacodes);
302  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DSparticle_code);
303  // Below true if charginos and neutralinos are included in coannihilations:
304  RD_spectrum_SUSY_DS5.setOption<bool>("CoannCharginosNeutralinos", true);
305  // Below true if sfermions are included in coannihilations:
306  RD_spectrum_SUSY_DS5.setOption<bool>("CoannSfermions", true);
307  // Maximum sparticle mass to be icluded in coannihilations, in units of DM mass:
308  RD_spectrum_SUSY_DS5.setOption<double>("CoannMaxMass", 1.6);
309  RD_spectrum_SUSY_DS5.reset_and_calculate();
310 
312  RD_spectrum_ordered_func.reset_and_calculate();
313 
314  RD_annrate_DS5prep_func.resolveDependency(&RD_spectrum_SUSY_DS5);
315  RD_annrate_DS5prep_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdmgev);
316  RD_annrate_DS5prep_func.reset_and_calculate();
317 
318  RD_eff_annrate_DS_MSSM.notifyOfModel("MSSM30atQ");
320  RD_eff_annrate_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsanwx);
321  RD_eff_annrate_DS_MSSM.reset_and_calculate();
322 
323  RD_oh2_DS5_general.resolveDependency(&RD_spectrum_ordered_func);
324  RD_oh2_DS5_general.resolveDependency(&RD_eff_annrate_DS_MSSM);
325  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdthlim);
326  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdtab);
327  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdeqn);
328  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdwintp);
329  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DSparticle_code);
330  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::widths);
331  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdmgev);
332  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdpth);
333  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdpars);
334  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdswitch);
335  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdlun);
336  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdpadd);
337  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rddof);
338  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rderrors);
339  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdtime);
340  RD_oh2_DS5_general.setOption<int>("fast", 1); // 0: normal; 1: fast; 2: dirty
341  RD_oh2_DS5_general.reset_and_calculate();
342  results["oh2"][current_backend] = RD_oh2_DS5_general(0);
343 
344  lnL_oh2_Simple.resolveDependency(&RD_oh2_DS5_general);
345  lnL_oh2_Simple.reset_and_calculate();
346  // Save the result
347  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
348 
349 
350  // ---- Set up basic internal structures for direct & indirect detection ----
351 
352  // Set up process catalog based on DarkSUSY annihilation rates
353  TH_ProcessCatalog_DS5_MSSM.resolveDependency(&createSpectrum);
354  TH_ProcessCatalog_DS5_MSSM.resolveDependency(&createDecays);
356  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dssigmav);
357  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBffdxdy);
358  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBhhdxdy);
359  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBwhdxdy);
360  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBwwdxdy);
361  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::IBintvars);
362  TH_ProcessCatalog_DS5_MSSM.reset_and_calculate();
363 
364  // Set generic WIMP mass object
365  mwimp_generic.resolveDependency(&TH_ProcessCatalog_DS5_MSSM);
366  mwimp_generic.resolveDependency(&DarkMatter_ID_MSSM);
367  mwimp_generic.reset_and_calculate();
368 
369  // Set generic annihilation rate in late universe (v->0 limit)
371  sigmav_late_universe.resolveDependency(&DarkMatter_ID_MSSM);
372  sigmav_late_universe.reset_and_calculate();
373  // Save the result
374  results["sigmav0"][current_backend] = sigmav_late_universe(0);
375 
376 
377  // ---- Direct detection -----
378 
379  // Calculate DD couplings with DarkSUSY
380  // DD_couplings_DarkSUSY_DS5.notifyOfModel("nuclear_params_fnq");
381  // DD_couplings_DarkSUSY_DS5.resolveDependency(&Models::nuclear_params_fnq::Functown::primary_parameters);
382  DD_couplings_DarkSUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DD_couplings);
383  DD_couplings_DarkSUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::mspctm);
384  DD_couplings_DarkSUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::ddcom);
385  // The below calculates the DD couplings using the full 1 loop calculation of
386  // Drees Nojiri Phys.Rev. D48 (1993) 3483
387  DD_couplings_DarkSUSY_DS5.setOption<bool>("loop", true);
388  // Setting the below to false approximates the squark propagator as 1/m_sq^2 to avoid poles.
389  // To reproduce numbers in Tables 11/12 of DarkBit paper (https://arxiv.org/abs/1705.07920), set "pole" to true.
390  DD_couplings_DarkSUSY_DS5.setOption<bool>("pole", false);
391  DD_couplings_DarkSUSY_DS5.reset_and_calculate();
392 
393  // Initialize DDCalc backend
394  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
395  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
396  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
397  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
398  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
399  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
400  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_DarkSUSY_DS5);
401  DDCalc_2_2_0_init.reset_and_calculate();
402 
403  // Calculate direct detection rates for LUX 2016
404  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
405  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
406  LUX_2016_Calc.reset_and_calculate();
407 
408  // Calculate direct detection likelihood for LUX 2016
409  LUX_2016_GetLogLikelihood.resolveDependency(&LUX_2016_Calc);
410  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
411  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
412  LUX_2016_GetLogLikelihood.reset_and_calculate();
413  // Save the result
414  results["LUX_2016_lnL"][current_backend] = LUX_2016_GetLogLikelihood(0);
415 
416  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
417  sigma_SI_p_simple.resolveDependency(&DD_couplings_DarkSUSY_DS5);
418  sigma_SI_p_simple.reset_and_calculate();
419  // Save the result
420  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
421 
422  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
423  sigma_SD_p_simple.resolveDependency(&DD_couplings_DarkSUSY_DS5);
424  sigma_SD_p_simple.reset_and_calculate();
425  // Save the result
426  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
427 
428 
429  // ---- Gamma-ray yields ----
430 
431  // Initialize tabulated gamma-ray yields
432  SimYieldTable_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshayield);
433  SimYieldTable_DS5.reset_and_calculate();
434 
435  // Collect missing final states for simulation in cascade MC
437  GA_missingFinalStates.resolveDependency(&SimYieldTable_DS5);
438  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_MSSM);
439  GA_missingFinalStates.reset_and_calculate();
440 
441  // Infer for which type of final states particles MC should be performed
442  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec<std::string>("gamma"));
443  cascadeMC_FinalStates.reset_and_calculate();
444 
445  // Collect decay information for cascade MC
447  cascadeMC_DecayTable.resolveDependency(&SimYieldTable_DS5);
448  cascadeMC_DecayTable.reset_and_calculate();
449 
450  // cascadeMC_LoopManager.setOption<int>("cMC_maxEvents", 100000);
451  // cascadeMC_Histograms.setOption<double>("cMC_endCheckFrequency", 25);
452  // cascadeMC_Histograms.setOption<double>("cMC_gammaRelError", .05);
453  // cascadeMC_Histograms.setOption<int>("cMC_numSpecSamples", 25);
454  // cascadeMC_Histograms.setOption<int>("cMC_NhistBins", 300);
455 
456  // Set up MC loop manager for cascade MC
457  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
458  std::vector<functor*> nested_functions = initVector<functor*>(
460  cascadeMC_LoopManager.setNestedList(nested_functions);
461 
462  // Set up initial state for cascade MC step
463  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
464  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
465 
466  // Perform MC step for cascade MC
467  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
468  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
469  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
470 
471  // Generate histogram for cascade MC
472  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
473  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
474  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_DS5_MSSM);
475  cascadeMC_Histograms.resolveDependency(&SimYieldTable_DS5);
476  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
477  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
478 
479  // Check convergence of cascade MC
480  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
481  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
482 
483  // Start cascade MC loop
484  cascadeMC_LoopManager.reset_and_calculate();
485 
486  // Infer gamma-ray spectra for recorded MC results
487  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
488  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
489  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
490  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
491  cascadeMC_gammaSpectra.reset_and_calculate();
492 
493  // Calculate total gamma-ray yield (cascade MC + tabulated results)
495  GA_AnnYield_General.resolveDependency(&SimYieldTable_DS5);
496  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_MSSM);
497  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
498  GA_AnnYield_General.reset_and_calculate();
499 
500  // Dump spectrum into file
501  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
502  dump_GammaSpectrum.setOption<std::string>("filename", current_backend + "_" + outname_dNdE_spectrum);
503  dump_GammaSpectrum.reset_and_calculate();
504 
505  // Calculate Fermi LAT dwarf likelihood
507  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
508  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
509  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
510  // Save the result
511  results["FermiLAT_dwarfsph_lnL"][current_backend] = lnL_FermiLATdwarfs_gamLike(0);
512 
513 
514  // ---- IceCube limits ----
515 
516  // Infer WIMP capture rate in Sun
520  capture_rate_Sun_const_xsec_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsntcapsuntab);
522  capture_rate_Sun_const_xsec_DS5.reset_and_calculate();
523 
524  // Infer WIMP equilibration time in Sun
526  equilibration_time_Sun.resolveDependency(&DarkMatter_ID_MSSM);
527  equilibration_time_Sun.resolveDependency(&mwimp_generic);
529  equilibration_time_Sun.reset_and_calculate();
530 
531  // Infer WIMP annihilation rate in Sun
532  annihilation_rate_Sun.resolveDependency(&equilibration_time_Sun);
534  annihilation_rate_Sun.reset_and_calculate();
535 
536  // Infer neutrino yield from Sun
537  nuyield_from_DS.resolveDependency(&TH_ProcessCatalog_DS5_MSSM);
538  nuyield_from_DS.resolveDependency(&mwimp_generic);
539  nuyield_from_DS.resolveDependency(&sigmav_late_universe);
540  nuyield_from_DS.resolveDependency(&DarkMatter_ID_MSSM);
541  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsgenericwimp_nusetup);
542  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::neutrino_yield);
543  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DS_neutral_h_decay_channels);
544  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DS_charged_h_decay_channels);
545  nuyield_from_DS.reset_and_calculate();
546 
547  // Calculate number of events at IceCube
548  IC79WH_full.resolveDependency(&mwimp_generic);
549  IC79WH_full.resolveDependency(&annihilation_rate_Sun);
550  IC79WH_full.resolveDependency(&nuyield_from_DS);
551  IC79WH_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
552  IC79WH_full.reset_and_calculate();
553  IC79WL_full.resolveDependency(&mwimp_generic);
554  IC79WL_full.resolveDependency(&annihilation_rate_Sun);
555  IC79WL_full.resolveDependency(&nuyield_from_DS);
556  IC79WL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
557  IC79WL_full.reset_and_calculate();
558  IC79SL_full.resolveDependency(&mwimp_generic);
559  IC79SL_full.resolveDependency(&annihilation_rate_Sun);
560  IC79SL_full.resolveDependency(&nuyield_from_DS);
561  IC79SL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
562  IC79SL_full.reset_and_calculate();
563 
564  // Calculate IceCube likelihood
565  IC79WH_bgloglike.resolveDependency(&IC79WH_full);
566  IC79WH_bgloglike.reset_and_calculate();
567  IC79WH_loglike.resolveDependency(&IC79WH_full);
568  IC79WH_loglike.reset_and_calculate();
569  IC79WL_bgloglike.resolveDependency(&IC79WL_full);
570  IC79WL_bgloglike.reset_and_calculate();
571  IC79WL_loglike.resolveDependency(&IC79WL_full);
572  IC79WL_loglike.reset_and_calculate();
573  IC79SL_bgloglike.resolveDependency(&IC79SL_full);
574  IC79SL_bgloglike.reset_and_calculate();
575  IC79SL_loglike.resolveDependency(&IC79SL_full);
576  IC79SL_loglike.reset_and_calculate();
577  IC79_loglike.resolveDependency(&IC79WH_bgloglike);
578  IC79_loglike.resolveDependency(&IC79WH_loglike);
579  IC79_loglike.resolveDependency(&IC79WL_bgloglike);
580  IC79_loglike.resolveDependency(&IC79WL_loglike);
581  IC79_loglike.resolveDependency(&IC79SL_bgloglike);
582  IC79_loglike.resolveDependency(&IC79SL_loglike);
583  IC79_loglike.reset_and_calculate();
584  // Save the result
585  results["IceCube_79_lnL"][current_backend] = IC79_loglike(0);
586 
587  } // End of DarkSUSY 5.1.3 calculations
588 
589 
590 
591 
592  //
593  // ======= Perform all calculations for backend DarkSUSY_MSSM 6.1.1 =======
594  //
595 
596  current_backend = "DarkSUSY_MSSM6.1.1";
597 
598  if (not Backends::backendInfo().works[current_backend])
599  {
600  backends_not_built.push_back(current_backend);
601  }
602  else
603  {
604  // Initialize DarkSUSY 6 MSSM backend
605  DarkSUSY_MSSM_6_1_1_init.notifyOfModel("MSSM30atQ");
606  DarkSUSY_MSSM_6_1_1_init.resolveDependency(&createSpectrum);
607  DarkSUSY_MSSM_6_1_1_init.resolveDependency(&createDecays);
608  if (decays) DarkSUSY_MSSM_6_1_1_init.setOption<bool>("use_dsSLHAread", false);
609  else DarkSUSY_MSSM_6_1_1_init.setOption<bool>("use_dsSLHAread", true);
610  DarkSUSY_MSSM_6_1_1_init.reset_and_calculate();
611 
612  // Initialize DarkSUSY 6 Local Halo Model
615  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmcom);
616  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmisodf);
617  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmframevelcom);
618  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmnoclue);
619  DarkSUSY_PointInit_LocalHalo_func.reset_and_calculate();
620 
621  // Relic density calculation with GAMBIT routines and DarkSUSY 6:
622  RD_spectrum_MSSM.resolveDependency(&createDecays);
623  RD_spectrum_MSSM.resolveDependency(&createSpectrum);
624  RD_spectrum_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
625  // Below true if charginos and neutralinos are included in coannihilations:
626  RD_spectrum_MSSM.setOption<bool>("CoannCharginosNeutralinos", true);
627  // Below true if sfermions are included in coannihilations:
628  RD_spectrum_MSSM.setOption<bool>("CoannSfermions", true);
629  // Maximum sparticle mass to be icluded in coannihilations, in units of DM mass:
630  RD_spectrum_MSSM.setOption<double>("CoannMaxMass", 1.6);
631  RD_spectrum_MSSM.reset_and_calculate();
632 
633  RD_spectrum_ordered_func.resolveDependency(&RD_spectrum_MSSM);
634  RD_spectrum_ordered_func.reset_and_calculate();
635 
637  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsancoann);
638  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::DSparticle_code);
639  RD_annrate_DSprep_MSSM_func.reset_and_calculate();
640 
641  RD_eff_annrate_DS_MSSM.notifyOfModel("MSSM30atQ");
643  RD_eff_annrate_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsanwx);
644  RD_eff_annrate_DS_MSSM.reset_and_calculate();
645 
646  RD_oh2_DS_general.resolveDependency(&RD_spectrum_ordered_func);
647  RD_oh2_DS_general.resolveDependency(&RD_eff_annrate_DS_MSSM);
648  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::rdpars);
649  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::rdtime);
650  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsrdcom);
651  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsrdstart);
652  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsrdens);
653  RD_oh2_DS_general.setOption<int>("fast", 1); // 0: normal; 1: fast; 2: dirty
654  RD_oh2_DS_general.reset_and_calculate();
655  // Save the result
656  results["oh2"][current_backend] = RD_oh2_DS_general(0);
657 
658  lnL_oh2_Simple.resolveDependency(&RD_oh2_DS_general);
659  lnL_oh2_Simple.reset_and_calculate();
660  // Save the result
661  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
662 
663 
664  // Set up process catalog based on DarkSUSY annihilation rates
665  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createSpectrum);
666  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createDecays);
667  TH_ProcessCatalog_DS_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
668  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dssigmav0);
669  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dssigmav0tot);
670  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBffdxdy);
671  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBhhdxdy);
672  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBwhdxdy);
673  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBwwdxdy);
674  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::IBintvars);
675  TH_ProcessCatalog_DS_MSSM.reset_and_calculate();
676 
677  // Set generic WIMP mass object
678  mwimp_generic.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
679  mwimp_generic.resolveDependency(&DarkMatter_ID_MSSM);
680  mwimp_generic.reset_and_calculate();
681 
682  // Set generic annihilation rate in late universe (v->0 limit)
684  sigmav_late_universe.resolveDependency(&DarkMatter_ID_MSSM);
685  sigmav_late_universe.reset_and_calculate();
686  // Save the result
687  results["sigmav0"][current_backend] = sigmav_late_universe(0);
688 
689 
690  // ---- Gamma-ray yields ----
691 
692  // Initialize tabulated gamma-ray yields
693  SimYieldTable_DarkSUSY.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsanyield_sim);
694  SimYieldTable_DarkSUSY.reset_and_calculate();
695 
696  // Collect missing final states for simulation in cascade MC
698  GA_missingFinalStates.resolveDependency(&SimYieldTable_DarkSUSY);
699  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_MSSM);
700  GA_missingFinalStates.reset_and_calculate();
701 
702  // Infer for which type of final states particles MC should be performed
703  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec<std::string>("gamma"));
704  cascadeMC_FinalStates.reset_and_calculate();
705 
706  // Collect decay information for cascade MC
708  cascadeMC_DecayTable.resolveDependency(&SimYieldTable_DarkSUSY);
709  cascadeMC_DecayTable.reset_and_calculate();
710 
711  // Set up MC loop manager for cascade MC
712  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
713  std::vector<functor*> nested_functions = initVector<functor*>(
715  cascadeMC_LoopManager.setNestedList(nested_functions);
716 
717  // Set up initial state for cascade MC step
718  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
719  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
720 
721  // Perform MC step for cascade MC
722  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
723  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
724  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
725 
726  // Generate histogram for cascade MC
727  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
728  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
729  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
730  cascadeMC_Histograms.resolveDependency(&SimYieldTable_DarkSUSY);
731  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
732  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
733 
734  // Check convergence of cascade MC
735  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
736  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
737 
738  // Start cascade MC loop
739  cascadeMC_LoopManager.reset_and_calculate();
740 
741  // Infer gamma-ray spectra for recorded MC results
742  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
743  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
744  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
745  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
746  cascadeMC_gammaSpectra.reset_and_calculate();
747 
748  // Calculate total gamma-ray yield (cascade MC + tabulated results)
750  GA_AnnYield_General.resolveDependency(&SimYieldTable_DarkSUSY);
751  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_MSSM);
752  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
753  GA_AnnYield_General.reset_and_calculate();
754 
755  // Dump spectrum into file
756  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
757  dump_GammaSpectrum.setOption<std::string>("filename", current_backend + "_" + outname_dNdE_spectrum);
758  dump_GammaSpectrum.reset_and_calculate();
759 
760  // Calculate Fermi LAT dwarf likelihood
762  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
763  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
764  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
765  // Save the result
766  results["FermiLAT_dwarfsph_lnL"][current_backend] = lnL_FermiLATdwarfs_gamLike(0);
767 
768 
769  // ---- Direct detection and IceCube limits ----
770 
771  // Calculate DD couplings with DarkSUSY
772  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::DD_couplings);
773  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::ddcomlegacy);
774  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::ddmssmcom);
775  // The below calculates the DD couplings using the full 1 loop calculation of
776  // Drees Nojiri Phys.Rev. D48 (1993) 3483
777  DD_couplings_DarkSUSY_MSSM.setOption<bool>("loop", true);
778  // Setting the below to false approximates the squark propagator as 1/m_sq^2 to avoid poles.
779  DD_couplings_DarkSUSY_MSSM.setOption<bool>("pole", false);
780  DD_couplings_DarkSUSY_MSSM.reset_and_calculate();
781 
782  // Initialize DDCalc backend
783  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
784  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
785  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
786  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
787  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
788  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
789  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
790  DDCalc_2_2_0_init.reset_and_calculate();
791 
792  // Calculate direct detection rates for LUX 2016
793  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
794  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
795  LUX_2016_Calc.reset_and_calculate();
796 
797  // Calculate direct detection likelihood for LUX 2016
798  LUX_2016_GetLogLikelihood.resolveDependency(&LUX_2016_Calc);
799  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
800  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
801  LUX_2016_GetLogLikelihood.reset_and_calculate();
802  // Save the result
803  results["LUX_2016_lnL"][current_backend] = LUX_2016_GetLogLikelihood(0);
804 
805 
806  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
807  sigma_SI_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
808  sigma_SI_p_simple.reset_and_calculate();
809  // Save the result
810  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
811 
812  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
813  sigma_SD_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
814  sigma_SD_p_simple.reset_and_calculate();
815  // Save the result
816  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
817 
818 
819  // Infer WIMP capture rate in Sun
820  capture_rate_Sun_const_xsec.resolveDependency(&mwimp_generic);
823  capture_rate_Sun_const_xsec.resolveDependency(&RD_fraction_one);
824  capture_rate_Sun_const_xsec.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dssenu_capsuntab);
827  capture_rate_Sun_const_xsec.reset_and_calculate();
828 
829  // Infer WIMP equilibration time in Sun
831  equilibration_time_Sun.resolveDependency(&DarkMatter_ID_MSSM);
832  equilibration_time_Sun.resolveDependency(&mwimp_generic);
834  equilibration_time_Sun.reset_and_calculate();
835 
836  // Infer WIMP annihilation rate in Sun
837  annihilation_rate_Sun.resolveDependency(&equilibration_time_Sun);
839  annihilation_rate_Sun.reset_and_calculate();
840 
841  // Infer neutrino yield from Sun
842  nuyield_from_DS.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
843  nuyield_from_DS.resolveDependency(&mwimp_generic);
844  nuyield_from_DS.resolveDependency(&sigmav_late_universe);
845  nuyield_from_DS.resolveDependency(&DarkMatter_ID_MSSM);
846  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsgenericwimp_nusetup);
847  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::neutrino_yield);
848  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::DS_neutral_h_decay_channels);
849  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::DS_charged_h_decay_channels);
850  nuyield_from_DS.reset_and_calculate();
851 
852 
853  // Calculate number of events at IceCube
854  IC79WH_full.resolveDependency(&mwimp_generic);
855  IC79WH_full.resolveDependency(&annihilation_rate_Sun);
856  IC79WH_full.resolveDependency(&nuyield_from_DS);
857  IC79WH_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
858  IC79WH_full.reset_and_calculate();
859  IC79WL_full.resolveDependency(&mwimp_generic);
860  IC79WL_full.resolveDependency(&annihilation_rate_Sun);
861  IC79WL_full.resolveDependency(&nuyield_from_DS);
862  IC79WL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
863  IC79WL_full.reset_and_calculate();
864  IC79SL_full.resolveDependency(&mwimp_generic);
865  IC79SL_full.resolveDependency(&annihilation_rate_Sun);
866  IC79SL_full.resolveDependency(&nuyield_from_DS);
867  IC79SL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
868  IC79SL_full.reset_and_calculate();
869 
870  // Calculate IceCube likelihood
871  IC79WH_bgloglike.resolveDependency(&IC79WH_full);
872  IC79WH_bgloglike.reset_and_calculate();
873  IC79WH_loglike.resolveDependency(&IC79WH_full);
874  IC79WH_loglike.reset_and_calculate();
875  IC79WL_bgloglike.resolveDependency(&IC79WL_full);
876  IC79WL_bgloglike.reset_and_calculate();
877  IC79WL_loglike.resolveDependency(&IC79WL_full);
878  IC79WL_loglike.reset_and_calculate();
879  IC79SL_bgloglike.resolveDependency(&IC79SL_full);
880  IC79SL_bgloglike.reset_and_calculate();
881  IC79SL_loglike.resolveDependency(&IC79SL_full);
882  IC79SL_loglike.reset_and_calculate();
883  IC79_loglike.resolveDependency(&IC79WH_bgloglike);
884  IC79_loglike.resolveDependency(&IC79WH_loglike);
885  IC79_loglike.resolveDependency(&IC79WL_bgloglike);
886  IC79_loglike.resolveDependency(&IC79WL_loglike);
887  IC79_loglike.resolveDependency(&IC79SL_bgloglike);
888  IC79_loglike.resolveDependency(&IC79SL_loglike);
889  IC79_loglike.reset_and_calculate();
890  // Save the result
891  results["IceCube_79_lnL"][current_backend] = IC79_loglike(0);
892 
893  } // End of DarkSUSY_MSSM 6.1.1 calculations
894 
895 
896 
897 
898  //
899  // ======= Perform all calculations for backend DarkSUSY_MSSM 6.2.2 =======
900  //
901 
902  current_backend = "DarkSUSY_MSSM6.2.2";
903 
904  if (not Backends::backendInfo().works[current_backend])
905  {
906  backends_not_built.push_back(current_backend);
907  }
908  else
909  {
910  // Initialize DarkSUSY 6 MSSM backend
911  DarkSUSY_MSSM_6_2_2_init.notifyOfModel("MSSM30atQ");
912  DarkSUSY_MSSM_6_2_2_init.resolveDependency(&createSpectrum);
913  DarkSUSY_MSSM_6_2_2_init.resolveDependency(&createDecays);
914  if (decays) DarkSUSY_MSSM_6_2_2_init.setOption<bool>("use_dsSLHAread", false);
915  else DarkSUSY_MSSM_6_2_2_init.setOption<bool>("use_dsSLHAread", true);
916  DarkSUSY_MSSM_6_2_2_init.reset_and_calculate();
917 
918  // Initialize DarkSUSY 6 Local Halo Model
921  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmcom);
922  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmisodf);
923  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmframevelcom);
924  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmnoclue);
925  DarkSUSY_PointInit_LocalHalo_func.reset_and_calculate();
926 
927  // Relic density calculation with GAMBIT routines and DarkSUSY 6:
928  RD_spectrum_MSSM.resolveDependency(&createDecays);
929  RD_spectrum_MSSM.resolveDependency(&createSpectrum);
930  RD_spectrum_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
931  // Below true if charginos and neutralinos are included in coannihilations:
932  RD_spectrum_MSSM.setOption<bool>("CoannCharginosNeutralinos", true);
933  // Below true if sfermions are included in coannihilations:
934  RD_spectrum_MSSM.setOption<bool>("CoannSfermions", true);
935  // Maximum sparticle mass to be icluded in coannihilations, in units of DM mass:
936  RD_spectrum_MSSM.setOption<double>("CoannMaxMass", 1.6);
937  RD_spectrum_MSSM.reset_and_calculate();
938 
939  RD_spectrum_ordered_func.resolveDependency(&RD_spectrum_MSSM);
940  RD_spectrum_ordered_func.reset_and_calculate();
941 
943  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsancoann);
944  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::DSparticle_code);
945  RD_annrate_DSprep_MSSM_func.reset_and_calculate();
946 
947  RD_eff_annrate_DS_MSSM.notifyOfModel("MSSM30atQ");
949  RD_eff_annrate_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsanwx);
950  RD_eff_annrate_DS_MSSM.reset_and_calculate();
951 
952  RD_oh2_DS_general.resolveDependency(&RD_spectrum_ordered_func);
953  RD_oh2_DS_general.resolveDependency(&RD_eff_annrate_DS_MSSM);
954  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::rdpars);
955  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::rdtime);
956  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsrdcom);
957  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsrdstart);
958  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsrdens);
959  RD_oh2_DS_general.setOption<int>("fast", 1); // 0: normal; 1: fast; 2: dirty
960  RD_oh2_DS_general.reset_and_calculate();
961  // Save the result
962  results["oh2"][current_backend] = RD_oh2_DS_general(0);
963 
964  lnL_oh2_Simple.resolveDependency(&RD_oh2_DS_general);
965  lnL_oh2_Simple.reset_and_calculate();
966  // Save the result
967  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
968 
969 
970  // Set up process catalog based on DarkSUSY annihilation rates
971  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createSpectrum);
972  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createDecays);
973  TH_ProcessCatalog_DS_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
974  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dssigmav0);
975  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dssigmav0tot);
976  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBffdxdy);
977  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBhhdxdy);
978  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBwhdxdy);
979  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBwwdxdy);
980  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::IBintvars);
981  TH_ProcessCatalog_DS_MSSM.reset_and_calculate();
982 
983  // Set generic WIMP mass object
984  mwimp_generic.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
985  mwimp_generic.resolveDependency(&DarkMatter_ID_MSSM);
986  mwimp_generic.reset_and_calculate();
987 
988  // Set generic annihilation rate in late universe (v->0 limit)
990  sigmav_late_universe.resolveDependency(&DarkMatter_ID_MSSM);
991  sigmav_late_universe.reset_and_calculate();
992  // Save the result
993  results["sigmav0"][current_backend] = sigmav_late_universe(0);
994 
995 
996  // ---- Gamma-ray yields ----
997 
998  // Initialize tabulated gamma-ray yields
999  SimYieldTable_DarkSUSY.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsanyield_sim);
1000  SimYieldTable_DarkSUSY.reset_and_calculate();
1001 
1002  // Collect missing final states for simulation in cascade MC
1004  GA_missingFinalStates.resolveDependency(&SimYieldTable_DarkSUSY);
1005  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_MSSM);
1006  GA_missingFinalStates.reset_and_calculate();
1007 
1008  // Infer for which type of final states particles MC should be performed
1009  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec<std::string>("gamma"));
1010  cascadeMC_FinalStates.reset_and_calculate();
1011 
1012  // Collect decay information for cascade MC
1013  cascadeMC_DecayTable.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1014  cascadeMC_DecayTable.resolveDependency(&SimYieldTable_DarkSUSY);
1015  cascadeMC_DecayTable.reset_and_calculate();
1016 
1017  // Set up MC loop manager for cascade MC
1018  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
1019  std::vector<functor*> nested_functions = initVector<functor*>(
1021  cascadeMC_LoopManager.setNestedList(nested_functions);
1022 
1023  // Set up initial state for cascade MC step
1024  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
1025  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
1026 
1027  // Perform MC step for cascade MC
1028  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
1029  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
1030  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
1031 
1032  // Generate histogram for cascade MC
1033  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
1034  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
1035  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1036  cascadeMC_Histograms.resolveDependency(&SimYieldTable_DarkSUSY);
1037  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
1038  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
1039 
1040  // Check convergence of cascade MC
1041  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
1042  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
1043 
1044  // Start cascade MC loop
1045  cascadeMC_LoopManager.reset_and_calculate();
1046 
1047  // Infer gamma-ray spectra for recorded MC results
1048  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
1049  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
1050  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
1051  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
1052  cascadeMC_gammaSpectra.reset_and_calculate();
1053 
1054  // Calculate total gamma-ray yield (cascade MC + tabulated results)
1055  GA_AnnYield_General.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1056  GA_AnnYield_General.resolveDependency(&SimYieldTable_DarkSUSY);
1057  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_MSSM);
1058  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
1059  GA_AnnYield_General.reset_and_calculate();
1060 
1061  // Dump spectrum into file
1062  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
1063  dump_GammaSpectrum.setOption<std::string>("filename", current_backend + "_" + outname_dNdE_spectrum);
1064  dump_GammaSpectrum.reset_and_calculate();
1065 
1066  // Calculate Fermi LAT dwarf likelihood
1067  lnL_FermiLATdwarfs_gamLike.resolveDependency(&GA_AnnYield_General);
1068  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
1069  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
1070  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
1071  // Save the result
1072  results["FermiLAT_dwarfsph_lnL"][current_backend] = lnL_FermiLATdwarfs_gamLike(0);
1073 
1074 
1075  // ---- Direct detection and IceCube limits ----
1076 
1077  // Calculate DD couplings with DarkSUSY
1078  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::DD_couplings);
1079  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::ddcomlegacy);
1080  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::ddmssmcom);
1081  // The below calculates the DD couplings using the full 1 loop calculation of
1082  // Drees Nojiri Phys.Rev. D48 (1993) 3483
1083  DD_couplings_DarkSUSY_MSSM.setOption<bool>("loop", true);
1084  // Setting the below to false approximates the squark propagator as 1/m_sq^2 to avoid poles.
1085  DD_couplings_DarkSUSY_MSSM.setOption<bool>("pole", false);
1086  DD_couplings_DarkSUSY_MSSM.reset_and_calculate();
1087 
1088  // Initialize DDCalc backend
1089  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
1090  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
1091  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
1092  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
1093  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
1094  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
1095  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1096  DDCalc_2_2_0_init.reset_and_calculate();
1097 
1098  // Calculate direct detection rates for LUX 2016
1099  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
1100  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
1101  LUX_2016_Calc.reset_and_calculate();
1102 
1103  // Calculate direct detection likelihood for LUX 2016
1104  LUX_2016_GetLogLikelihood.resolveDependency(&LUX_2016_Calc);
1105  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
1106  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
1107  LUX_2016_GetLogLikelihood.reset_and_calculate();
1108  // Save the result
1109  results["LUX_2016_lnL"][current_backend] = LUX_2016_GetLogLikelihood(0);
1110 
1111 
1112  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
1113  sigma_SI_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1114  sigma_SI_p_simple.reset_and_calculate();
1115  // Save the result
1116  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
1117 
1118  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
1119  sigma_SD_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1120  sigma_SD_p_simple.reset_and_calculate();
1121  // Save the result
1122  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
1123 
1124 
1125  // Infer WIMP capture rate in Sun
1126  capture_rate_Sun_const_xsec.resolveDependency(&mwimp_generic);
1127  capture_rate_Sun_const_xsec.resolveDependency(&sigma_SI_p_simple);
1128  capture_rate_Sun_const_xsec.resolveDependency(&sigma_SD_p_simple);
1129  capture_rate_Sun_const_xsec.resolveDependency(&RD_fraction_one);
1130  capture_rate_Sun_const_xsec.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dssenu_capsuntab);
1133  capture_rate_Sun_const_xsec.reset_and_calculate();
1134 
1135  // Infer WIMP equilibration time in Sun
1137  equilibration_time_Sun.resolveDependency(&DarkMatter_ID_MSSM);
1138  equilibration_time_Sun.resolveDependency(&mwimp_generic);
1140  equilibration_time_Sun.reset_and_calculate();
1141 
1142  // Infer WIMP annihilation rate in Sun
1143  annihilation_rate_Sun.resolveDependency(&equilibration_time_Sun);
1145  annihilation_rate_Sun.reset_and_calculate();
1146 
1147  // Infer neutrino yield from Sun
1148  nuyield_from_DS.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1149  nuyield_from_DS.resolveDependency(&mwimp_generic);
1150  nuyield_from_DS.resolveDependency(&sigmav_late_universe);
1151  nuyield_from_DS.resolveDependency(&DarkMatter_ID_MSSM);
1152  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsgenericwimp_nusetup);
1153  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::neutrino_yield);
1154  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::DS_neutral_h_decay_channels);
1155  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::DS_charged_h_decay_channels);
1156  nuyield_from_DS.reset_and_calculate();
1157 
1158 
1159  // Calculate number of events at IceCube
1160  IC79WH_full.resolveDependency(&mwimp_generic);
1161  IC79WH_full.resolveDependency(&annihilation_rate_Sun);
1162  IC79WH_full.resolveDependency(&nuyield_from_DS);
1163  IC79WH_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1164  IC79WH_full.reset_and_calculate();
1165  IC79WL_full.resolveDependency(&mwimp_generic);
1166  IC79WL_full.resolveDependency(&annihilation_rate_Sun);
1167  IC79WL_full.resolveDependency(&nuyield_from_DS);
1168  IC79WL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1169  IC79WL_full.reset_and_calculate();
1170  IC79SL_full.resolveDependency(&mwimp_generic);
1171  IC79SL_full.resolveDependency(&annihilation_rate_Sun);
1172  IC79SL_full.resolveDependency(&nuyield_from_DS);
1173  IC79SL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1174  IC79SL_full.reset_and_calculate();
1175 
1176  // Calculate IceCube likelihood
1177  IC79WH_bgloglike.resolveDependency(&IC79WH_full);
1178  IC79WH_bgloglike.reset_and_calculate();
1179  IC79WH_loglike.resolveDependency(&IC79WH_full);
1180  IC79WH_loglike.reset_and_calculate();
1181  IC79WL_bgloglike.resolveDependency(&IC79WL_full);
1182  IC79WL_bgloglike.reset_and_calculate();
1183  IC79WL_loglike.resolveDependency(&IC79WL_full);
1184  IC79WL_loglike.reset_and_calculate();
1185  IC79SL_bgloglike.resolveDependency(&IC79SL_full);
1186  IC79SL_bgloglike.reset_and_calculate();
1187  IC79SL_loglike.resolveDependency(&IC79SL_full);
1188  IC79SL_loglike.reset_and_calculate();
1189  IC79_loglike.resolveDependency(&IC79WH_bgloglike);
1190  IC79_loglike.resolveDependency(&IC79WH_loglike);
1191  IC79_loglike.resolveDependency(&IC79WL_bgloglike);
1192  IC79_loglike.resolveDependency(&IC79WL_loglike);
1193  IC79_loglike.resolveDependency(&IC79SL_bgloglike);
1194  IC79_loglike.resolveDependency(&IC79SL_loglike);
1195  IC79_loglike.reset_and_calculate();
1196  // Save the result
1197  results["IceCube_79_lnL"][current_backend] = IC79_loglike(0);
1198 
1199  } // End of DarkSUSY_MSSM 6.2.2 calculations
1200 
1201 
1202 
1203 
1204  //
1205  // ======= Perform all calculations for backend MicrOmegas_MSSM 3.6.9.2 =======
1206  //
1207 
1208  current_backend = "MicrOmegas_MSSM3.6.9.2";
1209 
1210  if (not Backends::backendInfo().works[current_backend])
1211  {
1212  backends_not_built.push_back(current_backend);
1213  }
1214  else
1215  {
1216  // Initialize MicrOmegas backend
1217  MicrOmegas_MSSM_3_6_9_2_init.notifyOfModel("MSSM30atQ");
1218  MicrOmegas_MSSM_3_6_9_2_init.resolveDependency(&createSpectrum);
1219  MicrOmegas_MSSM_3_6_9_2_init.resolveDependency(&createDecays);
1220  MicrOmegas_MSSM_3_6_9_2_init.resolveDependency(&createSLHA1Names);
1221  // Use decay table if it is present:
1222  if (decays) MicrOmegas_MSSM_3_6_9_2_init.setOption<bool>("internal_decays", false);
1223  else MicrOmegas_MSSM_3_6_9_2_init.setOption<bool>("internal_decays", true);
1224  MicrOmegas_MSSM_3_6_9_2_init.reset_and_calculate();
1225  // For the below VXdecay = 0 - no 3 body final states via virtual X
1226  // 1 - annihilations to 3 body final states via virtual X
1227  // 2 - (co)annihilations to 3 body final states via virtual X
1228  MicrOmegas_MSSM_3_6_9_2_init.setOption<int>("VZdecay", 0);
1229  MicrOmegas_MSSM_3_6_9_2_init.setOption<int>("VWdecay", 0);
1230  MicrOmegas_MSSM_3_6_9_2_init.reset_and_calculate();
1231 
1232  // Relic density calculation with MicrOmegas
1233  RD_oh2_Xf_MicrOmegas.notifyOfModel("MSSM30atQ");
1234  RD_oh2_Xf_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::darkOmega);
1235  RD_oh2_Xf_MicrOmegas.setOption<int>("fast", 1); // 0: accurate; 1: fast
1236  RD_oh2_Xf_MicrOmegas.setOption<double>("Beps", 1e-5); // Beps=1e-5 recommended, Beps=1 switches coannihilation off
1237  RD_oh2_Xf_MicrOmegas.reset_and_calculate();
1238  RD_oh2_MicrOmegas.resolveDependency(&RD_oh2_Xf_MicrOmegas);
1239  RD_oh2_MicrOmegas.reset_and_calculate();
1240  // Save the result
1241  results["oh2"][current_backend] = RD_oh2_MicrOmegas(0);
1242 
1243  lnL_oh2_Simple.resolveDependency(&RD_oh2_MicrOmegas);
1244  lnL_oh2_Simple.reset_and_calculate();
1245  // Save the result
1246  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
1247 
1248  // <sigma v> (v->0 limit) self-annihilation calculation with MicrOmegas:
1249  sigmav_late_universe_MicrOmegas.notifyOfModel("MSSM30atQ");
1250  sigmav_late_universe_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::calcSpectrum);
1251  sigmav_late_universe_MicrOmegas.reset_and_calculate();
1252  results["sigmav0"][current_backend] = sigmav_late_universe_MicrOmegas(0);
1253 
1254  // Direct detection calculations with Micromegas
1255 
1256  // Calculate DD couplings with Micromegas
1257  // DD_couplings_MicrOmegas.notifyOfModel("MSSM30atQ");
1258  // DD_couplings_MicrOmegas.notifyOfModel("nuclear_params_fnq");
1259  // DD_couplings_MicrOmegas.resolveDependency(&Models::nuclear_params_fnq::Functown::primary_parameters);
1260  DD_couplings_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::nucleonAmplitudes);
1261  DD_couplings_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::FeScLoop);
1262  DD_couplings_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::mocommon_);
1263  // The below includes neutralino-gluon scattering via a box diagram
1264  DD_couplings_MicrOmegas.setOption<bool>("box", true);
1265  DD_couplings_MicrOmegas.reset_and_calculate();
1266 
1267  // Initialize DDCalc backend
1268  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
1269  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
1270  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
1271  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
1272  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
1273  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
1274  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_MicrOmegas);
1275  DDCalc_2_2_0_init.reset_and_calculate();
1276 
1277 
1278  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
1279  sigma_SI_p_simple.resolveDependency(&DD_couplings_MicrOmegas);
1280  sigma_SI_p_simple.reset_and_calculate();
1281  // Save the result
1282  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
1283 
1284  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
1285  sigma_SD_p_simple.resolveDependency(&DD_couplings_MicrOmegas);
1286  sigma_SD_p_simple.reset_and_calculate();
1287  // Save the result
1288  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
1289 
1290  } // End of MicrOmegas_MSSM 3.6.9.2 calculations
1291 
1292 
1293 
1294  //
1295  // ======= Construct the output string =======
1296  //
1297 
1298  std::stringstream results_ss;
1299 
1300  for(const std::string& result_key : result_output_order)
1301  {
1302  const std::map<std::string,double>& backends_result_map = results.at(result_key);
1303  results_ss << result_key;
1304  if (results_units.at(result_key) != "") { results_ss << " [" << results_units.at(result_key) << "]"; }
1305  results_ss << " :" << endl;
1306 
1307  for(const auto& kv : backends_result_map)
1308  {
1309  const std::string& backendname = kv.first;
1310  const double& result = kv.second;
1311  results_ss << " " << backendname << ": " << result << " " << results_units.at(result_key) << endl;
1312  }
1313  results_ss << endl;
1314  }
1315 
1316 
1317  //
1318  // ======= Output the result string to screen =======
1319  //
1320 
1321  cout << endl;
1322  cout << "==== RESULTS ====" << endl;
1323  cout << endl;
1324  cout << results_ss.str();
1325  cout << endl;
1326 
1327  // Let the user know what they are missing...
1328  if (backends_not_built.size() > 0)
1329  {
1330  cout << endl;
1331  cout << "NOTE: The following backend(s) are not present:" << endl;
1332  for (const std::string& backend_name : backends_not_built)
1333  {
1334  cout << " - " << backend_name << endl;
1335  }
1336  cout << "If you want results from these backends you need to build them first." << endl;
1337  cout << endl;
1338  }
1339 
1340 
1341  //
1342  // ======= Output the result string to file =======
1343  //
1344 
1345  std::fstream file;
1346  file.open(outname_data, std::ios_base::out);
1347  file << results_ss.str();
1348  file.close();
1349 
1350  }
1351 
1352  catch (std::exception& e)
1353  {
1354  std::cout << "DarkBit_standalone_MSSM has exited with fatal exception: " << e.what() << std::endl;
1355  }
1356 
1357  return 0;
1358 
1359 }
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 IC79WL_full(nudata &result)
79-string IceCube WL sample: predicted signal and background counts, observed counts and likelihoods...
void DarkSUSY5_PointInit_LocalHalo_func(bool &result)
Function to set Local Halo Parameters in DarkSUSY (DS5 only)
void equilibration_time_Sun(double &result)
Equilibration time for capture and annihilation of dark matter in the Sun (s)
int main(int argc, char *argv[])
void sigmav_late_universe(double &result)
Retrieve the total thermally-averaged annihilation cross-section for indirect detection (cm^3 / s)...
Definition: DarkBit.cpp:70
void DarkSUSY_PointInit_LocalHalo_func(bool &result)
Function to set Local Halo Parameters in DarkSUSY (DS 6)
void IC79SL_bgloglike(double &result)
#define NEW_CAPABILITY
void sigma_SI_p_simple(double &result)
Simple calculator of the spin-independent WIMP-proton cross-section.
void SimYieldTable_DS5(SimYieldTable &result)
SimYieldTable based on DarkSUSY5 tabulated results. (DS6 below)
Definition: GamYields.cpp:425
void RD_annrate_DS5prep_func(int &result)
Some helper function to prepare evaluation of Weff from DarkSUSY 5.
Includes everything needed to use a GAMBIT module as a standalone analysis code.
void IC79WL_bgloglike(double &result)
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
void cascadeMC_InitialState(std::string &pID)
Function selecting initial state for decay chain.
Definition: Cascades.cpp:140
void RD_spectrum_SUSY_DS5(RD_spectrum_type &result)
Collects spectrum information about coannihilating particles, resonances and threshold energies – di...
void IC79WL_loglike(double &result)
void dump_GammaSpectrum(double &result)
Helper function to dump gamma-ray spectra.
void DD_couplings_DarkSUSY_DS5(DM_nucleon_couplings &result)
Get direct detection couplings from initialized DarkSUSY 5.
void RD_oh2_Xf_MicrOmegas(ddpair &result)
Relic density directly from a call of initialized MicrOmegas.
void createDecays(DecayTable &outDecays)
Helpers for making simple spectra from SLHAea objects.
void nuyield_from_DS(nuyield_info &result)
Neutrino yield function pointer and setup.
void annihilation_rate_Sun(double &result)
Annihilation rate of dark matter in the Sun (s^-1)
void IC79WH_loglike(double &result)
void RD_oh2_DS5_general(double &result)
General routine for calculation of relic density, using DarkSUSY 5 Boltzmann solver.
void IC79SL_loglike(double &result)
void RD_spectrum_MSSM(RD_spectrum_type &result)
Collects spectrum information about coannihilating particles, resonances and threshold energies...
warning & model_warning()
Model warnings.
Routines to help translate between SLHA2 sfermions and SLHA1 (or similar) sfermions.
std::map< str, daFunk::Funk > stringFunkMap
Definition: SimpleHist.hpp:101
General small utility functions.
void lnL_oh2_Simple(double &result)
Likelihood for cosmological relic density constraints.
void DD_couplings_MicrOmegas(DM_nucleon_couplings &result)
Get direct detection couplings from initialized MicrOmegas.
error & backend_error()
Backend errors.
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.
void IC79SL_full(nudata &result)
79-string IceCube SL sample: predicted signal and background counts, observed counts and likelihoods...
void cascadeMC_GenerateChain(DarkBit::DecayChain::ChainContainer &chain)
Function for generating decay chains.
Definition: Cascades.cpp:198
void capture_rate_Sun_const_xsec_DS5(double &result)
Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1)...
void capture_rate_Sun_const_xsec(double &result)
Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1)...
void GA_AnnYield_General(daFunk::Funk &result)
General routine to derive annihilation yield.
Definition: GamYields.cpp:191
void initialise_standalone_logs(str)
Logger setup standalone utility function.
void lnL_FermiLATdwarfs_gamLike(double &result)
Fermi LAT dwarf likelihoods, using gamLike backend.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
void DD_couplings_DarkSUSY_MSSM(DM_nucleon_couplings &result)
Get direct detection couplings from DarkSUSY 6 initialized with MSSM module.
void sigmav_late_universe_MicrOmegas(double &result)
Definition: DarkBit.cpp:90
void IC79WH_full(nudata &result)
79-string IceCube WH sample: predicted signal and background counts, observed counts and likelihoods...
void IC79_loglike(double &result)
Composite IceCube 79-string likelihood function.
void ExtractLocalMaxwellianHalo(LocalMaxwellianHalo &result)
Module function providing local density and velocity dispersion parameters.
Definition: DarkBit.cpp:152
Rollcall header for module DarkBit.
void IC79WH_bgloglike(double &result)
void TH_ProcessCatalog_DS_MSSM(DarkBit::TH_ProcessCatalog &result)
Initialization of Process Catalog based on DarkSUSY 6 calculations.
Definition: MSSM.cpp:417
GAMBIT native decay table class.
Definition: decay_table.hpp:35
void createSLHA1Names(mass_es_pseudonyms &names)
std::vector< YAML::sdd > mc_info
Typedefs for making it easier to manipulate mass cut and mass ratio cut info.
Definition: spectrum.hpp:119
void mwimp_generic(double &result)
Retrieve the DM mass in GeV for generic models (GeV)
Definition: DarkBit.cpp:60
void sigma_SD_p_simple(double &result)
Simple calculator of the spin-dependent WIMP-proton cross-section.
void cascadeMC_LoopManager()
Loop manager for cascade decays.
Definition: Cascades.cpp:80
void DarkMatter_ID_MSSM(std::string &result)
Definition: MSSM.cpp:743
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 RD_oh2_MicrOmegas(double &result)
void cascadeMC_DecayTable(DarkBit::DecayChain::DecayTable &table)
Function setting up the decay table used in decay chains.
Definition: Cascades.cpp:58
void RD_oh2_DS_general(double &result)
General routine for calculation of relic density, using DarkSUSY 6+ Boltzmann solver.
void RD_eff_annrate_DS_MSSM(double(*&result)(double &))
Get Weff directly from initialized DarkSUSY. Note that these functions do not (and should not) correc...
std::vector< YAML::ssdd > mr_info
Definition: spectrum.hpp:120
QUICK_FUNCTION(DarkBit, decay_rates, NEW_CAPABILITY, createDecays, DecayTable,()) QUICK_FUNCTION(DarkBit
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry decay_rates
void TH_ProcessCatalog_DS5_MSSM(DarkBit::TH_ProcessCatalog &result)
Initialization of Process Catalog based on DarkSUSY 5 calculations.
Definition: MSSM.cpp:100
void GalacticHalo_Einasto(GalacticHaloProperties &result)
Module function to generate GalacticHaloProperties for Einasto profile.
Definition: DarkBit.cpp:140
void RD_annrate_DSprep_MSSM_func(int &result)
Some helper function to prepare evaluation of Weff from DarkSUSY 6.
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
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110