gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
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  DarkMatterConj_ID_MSSM.resolveDependency(&createSpectrum);
250  DarkMatterConj_ID_MSSM.reset_and_calculate();
251 
252  // Assume for direct and indirect detection likelihoods that dark matter
253  // density is always the measured one (regardless of relic density results)
254  RD_fraction_one.reset_and_calculate();
255 
256  //
257  // ======= Initializations that can be done once =======
258  //
259 
260  // Initialize nulike backend
261  Backends::nulike_1_0_9::Functown::nulike_bounds.setStatus(2);
262  nulike_1_0_9_init.reset_and_calculate();
263 
264  // Initialize gamLike backend
265  gamLike_1_0_1_init.reset_and_calculate();
266 
267 
268 
269  //
270  // ======= Perform all calculations for backend DarkSUSY 5.1.3 =======
271  //
272 
273  current_backend = "DarkSUSY5.1.3";
274 
275  if (not Backends::backendInfo().works[current_backend])
276  {
277  backends_not_built.push_back(current_backend);
278  }
279  else
280  {
281  // Initialize DarkSUSY 5 backend
282  DarkSUSY_5_1_3_init.notifyOfModel("MSSM30atQ");
283  DarkSUSY_5_1_3_init.resolveDependency(&createSpectrum);
284  DarkSUSY_5_1_3_init.resolveDependency(&createDecays);
285  if (decays) DarkSUSY_5_1_3_init.setOption<bool>("use_dsSLHAread", false);
286  else DarkSUSY_5_1_3_init.setOption<bool>("use_dsSLHAread", true);
287  DarkSUSY_5_1_3_init.reset_and_calculate();
288 
289  // Initialize DarkSUSY 5 Local Halo Model
292  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmcom);
293  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmisodf);
294  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmframevelcom);
295  DarkSUSY5_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshmnoclue);
296  DarkSUSY5_PointInit_LocalHalo_func.reset_and_calculate();
297 
298 
299  // Relic density calculation with GAMBIT routines and DarkSUSY 5:
300  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::mspctm);
301  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::widths);
302  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::intdof);
303  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::pacodes);
304  RD_spectrum_SUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DSparticle_code);
305  // Below true if charginos and neutralinos are included in coannihilations:
306  RD_spectrum_SUSY_DS5.setOption<bool>("CoannCharginosNeutralinos", true);
307  // Below true if sfermions are included in coannihilations:
308  RD_spectrum_SUSY_DS5.setOption<bool>("CoannSfermions", true);
309  // Maximum sparticle mass to be icluded in coannihilations, in units of DM mass:
310  RD_spectrum_SUSY_DS5.setOption<double>("CoannMaxMass", 1.6);
311  RD_spectrum_SUSY_DS5.reset_and_calculate();
312 
314  RD_spectrum_ordered_func.reset_and_calculate();
315 
316  RD_annrate_DS5prep_func.resolveDependency(&RD_spectrum_SUSY_DS5);
317  RD_annrate_DS5prep_func.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdmgev);
318  RD_annrate_DS5prep_func.reset_and_calculate();
319 
320  RD_eff_annrate_DS5_MSSM.notifyOfModel("MSSM30atQ");
322  RD_eff_annrate_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsanwx);
323  RD_eff_annrate_DS5_MSSM.reset_and_calculate();
324 
325  RD_oh2_DS5_general.resolveDependency(&RD_spectrum_ordered_func);
326  RD_oh2_DS5_general.resolveDependency(&RD_eff_annrate_DS5_MSSM);
327  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdthlim);
328  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdtab);
329  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdeqn);
330  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsrdwintp);
331  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DSparticle_code);
332  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::widths);
333  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdmgev);
334  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdpth);
335  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdpars);
336  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdswitch);
337  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdlun);
338  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdpadd);
339  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rddof);
340  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rderrors);
341  RD_oh2_DS5_general.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::rdtime);
342  RD_oh2_DS5_general.setOption<int>("fast", 1); // 0: normal; 1: fast; 2: dirty
343  RD_oh2_DS5_general.reset_and_calculate();
344  results["oh2"][current_backend] = RD_oh2_DS5_general(0);
345 
346  lnL_oh2_Simple.resolveDependency(&RD_oh2_DS5_general);
347  lnL_oh2_Simple.reset_and_calculate();
348  // Save the result
349  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
350 
351 
352  // ---- Set up basic internal structures for direct & indirect detection ----
353 
354  // Set up process catalog based on DarkSUSY annihilation rates
355  TH_ProcessCatalog_DS5_MSSM.resolveDependency(&createSpectrum);
356  TH_ProcessCatalog_DS5_MSSM.resolveDependency(&createDecays);
358  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dssigmav);
359  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBffdxdy);
360  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBhhdxdy);
361  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBwhdxdy);
362  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsIBwwdxdy);
363  TH_ProcessCatalog_DS5_MSSM.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::IBintvars);
364  TH_ProcessCatalog_DS5_MSSM.reset_and_calculate();
365 
366  // Set generic WIMP mass object
367  mwimp_generic.resolveDependency(&TH_ProcessCatalog_DS5_MSSM);
368  mwimp_generic.resolveDependency(&DarkMatter_ID_MSSM);
369  mwimp_generic.reset_and_calculate();
370 
371  // Set generic annihilation rate in late universe (v->0 limit)
373  sigmav_late_universe.resolveDependency(&DarkMatter_ID_MSSM);
374  sigmav_late_universe.resolveDependency(&DarkMatterConj_ID_MSSM);
375  sigmav_late_universe.reset_and_calculate();
376  // Save the result
377  results["sigmav0"][current_backend] = sigmav_late_universe(0);
378 
379 
380  // ---- Direct detection -----
381 
382  // Calculate DD couplings with DarkSUSY
383  // DD_couplings_DarkSUSY_DS5.notifyOfModel("nuclear_params_fnq");
384  // DD_couplings_DarkSUSY_DS5.resolveDependency(&Models::nuclear_params_fnq::Functown::primary_parameters);
386  DD_couplings_DarkSUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::mspctm);
387  DD_couplings_DarkSUSY_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::ddcom);
388  // The below calculates the DD couplings using the full 1 loop calculation of
389  // Drees Nojiri Phys.Rev. D48 (1993) 3483
390  DD_couplings_DarkSUSY_DS5.setOption<bool>("loop", true);
391  // Setting the below to false approximates the squark propagator as 1/m_sq^2 to avoid poles.
392  // To reproduce numbers in Tables 11/12 of DarkBit paper (https://arxiv.org/abs/1705.07920), set "pole" to true.
393  DD_couplings_DarkSUSY_DS5.setOption<bool>("pole", false);
394  DD_couplings_DarkSUSY_DS5.reset_and_calculate();
395 
396  // Initialize DDCalc backend
397  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
398  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
399  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
400  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
401  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
402  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
403  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_DarkSUSY_DS5);
404  DDCalc_2_2_0_init.reset_and_calculate();
405 
406  // Calculate direct detection rates for LUX 2016
407  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
408  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
409  LUX_2016_Calc.reset_and_calculate();
410 
411  // Calculate direct detection likelihood for LUX 2016
412  LUX_2016_GetLogLikelihood.resolveDependency(&LUX_2016_Calc);
413  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
414  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
415  LUX_2016_GetLogLikelihood.reset_and_calculate();
416  // Save the result
417  results["LUX_2016_lnL"][current_backend] = LUX_2016_GetLogLikelihood(0);
418 
419  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
420  sigma_SI_p_simple.resolveDependency(&DD_couplings_DarkSUSY_DS5);
421  sigma_SI_p_simple.reset_and_calculate();
422  // Save the result
423  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
424 
425  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
426  sigma_SD_p_simple.resolveDependency(&DD_couplings_DarkSUSY_DS5);
427  sigma_SD_p_simple.reset_and_calculate();
428  // Save the result
429  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
430 
431 
432  // ---- Gamma-ray yields ----
433 
434  // Initialize tabulated gamma-ray yields
435  SimYieldTable_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dshayield);
436  SimYieldTable_DS5.reset_and_calculate();
437 
438  // Identify process as annihilation rather than decay
441  DM_process_from_ProcessCatalog.reset_and_calculate();
442 
443  // Collect missing final states for simulation in cascade MC
445  GA_missingFinalStates.resolveDependency(&SimYieldTable_DS5);
446  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_MSSM);
447  GA_missingFinalStates.resolveDependency(&DarkMatterConj_ID_MSSM);
449  GA_missingFinalStates.reset_and_calculate();
450 
451  // Infer for which type of final states particles MC should be performed
452  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec<std::string>("gamma"));
453  cascadeMC_FinalStates.reset_and_calculate();
454 
455  // Collect decay information for cascade MC
457  cascadeMC_DecayTable.resolveDependency(&SimYieldTable_DS5);
458  cascadeMC_DecayTable.reset_and_calculate();
459 
460  // cascadeMC_LoopManager.setOption<int>("cMC_maxEvents", 100000);
461  // cascadeMC_Histograms.setOption<double>("cMC_endCheckFrequency", 25);
462  // cascadeMC_Histograms.setOption<double>("cMC_gammaRelError", .05);
463  // cascadeMC_Histograms.setOption<int>("cMC_numSpecSamples", 25);
464  // cascadeMC_Histograms.setOption<int>("cMC_NhistBins", 300);
465 
466  // Set up MC loop manager for cascade MC
467  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
468  std::vector<functor*> nested_functions = initVector<functor*>(
470  cascadeMC_LoopManager.setNestedList(nested_functions);
471 
472  // Set up initial state for cascade MC step
473  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
474  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
475 
476  // Perform MC step for cascade MC
477  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
478  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
479  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
480 
481  // Generate histogram for cascade MC
482  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
483  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
484  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_DS5_MSSM);
485  cascadeMC_Histograms.resolveDependency(&SimYieldTable_DS5);
486  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
487  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
488 
489  // Check convergence of cascade MC
490  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
491  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
492 
493  // Start cascade MC loop
494  cascadeMC_LoopManager.reset_and_calculate();
495 
496  // Infer gamma-ray spectra for recorded MC results
497  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
498  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
499  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
500  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
501  cascadeMC_gammaSpectra.reset_and_calculate();
502 
503  // Calculate total gamma-ray yield (cascade MC + tabulated results)
505  GA_AnnYield_General.resolveDependency(&SimYieldTable_DS5);
506  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_MSSM);
507  GA_AnnYield_General.resolveDependency(&DarkMatterConj_ID_MSSM);
508  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
509  GA_AnnYield_General.reset_and_calculate();
510 
511  // Dump spectrum into file
512  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
513  dump_GammaSpectrum.setOption<std::string>("filename", current_backend + "_" + outname_dNdE_spectrum);
514  dump_GammaSpectrum.reset_and_calculate();
515 
516  // Calculate Fermi LAT dwarf likelihood
518  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
519  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
520  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
521  // Save the result
522  results["FermiLAT_dwarfsph_lnL"][current_backend] = lnL_FermiLATdwarfs_gamLike(0);
523 
524 
525  // ---- IceCube limits ----
526 
527  // Infer WIMP capture rate in Sun
531  capture_rate_Sun_const_xsec_DS5.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsntcapsuntab);
533  capture_rate_Sun_const_xsec_DS5.reset_and_calculate();
534 
535  // Infer WIMP equilibration time in Sun
537  equilibration_time_Sun.resolveDependency(&DarkMatter_ID_MSSM);
539  equilibration_time_Sun.resolveDependency(&mwimp_generic);
541  equilibration_time_Sun.reset_and_calculate();
542 
543  // Infer WIMP annihilation rate in Sun
544  annihilation_rate_Sun.resolveDependency(&equilibration_time_Sun);
546  annihilation_rate_Sun.reset_and_calculate();
547 
548  // Infer neutrino yield from Sun
549  nuyield_from_DS.resolveDependency(&TH_ProcessCatalog_DS5_MSSM);
550  nuyield_from_DS.resolveDependency(&mwimp_generic);
551  nuyield_from_DS.resolveDependency(&sigmav_late_universe);
552  nuyield_from_DS.resolveDependency(&DarkMatter_ID_MSSM);
553  nuyield_from_DS.resolveDependency(&DarkMatterConj_ID_MSSM);
554  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::dsgenericwimp_nusetup);
555  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::neutrino_yield);
556  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DS_neutral_h_decay_channels);
557  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_5_1_3::Functown::DS_charged_h_decay_channels);
558  nuyield_from_DS.reset_and_calculate();
559 
560  // Calculate number of events at IceCube
561  IC79WH_full.resolveDependency(&mwimp_generic);
562  IC79WH_full.resolveDependency(&annihilation_rate_Sun);
563  IC79WH_full.resolveDependency(&nuyield_from_DS);
564  IC79WH_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
565  IC79WH_full.reset_and_calculate();
566  IC79WL_full.resolveDependency(&mwimp_generic);
567  IC79WL_full.resolveDependency(&annihilation_rate_Sun);
568  IC79WL_full.resolveDependency(&nuyield_from_DS);
569  IC79WL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
570  IC79WL_full.reset_and_calculate();
571  IC79SL_full.resolveDependency(&mwimp_generic);
572  IC79SL_full.resolveDependency(&annihilation_rate_Sun);
573  IC79SL_full.resolveDependency(&nuyield_from_DS);
574  IC79SL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
575  IC79SL_full.reset_and_calculate();
576 
577  // Calculate IceCube likelihood
578  IC79WH_bgloglike.resolveDependency(&IC79WH_full);
579  IC79WH_bgloglike.reset_and_calculate();
580  IC79WH_loglike.resolveDependency(&IC79WH_full);
581  IC79WH_loglike.reset_and_calculate();
582  IC79WL_bgloglike.resolveDependency(&IC79WL_full);
583  IC79WL_bgloglike.reset_and_calculate();
584  IC79WL_loglike.resolveDependency(&IC79WL_full);
585  IC79WL_loglike.reset_and_calculate();
586  IC79SL_bgloglike.resolveDependency(&IC79SL_full);
587  IC79SL_bgloglike.reset_and_calculate();
588  IC79SL_loglike.resolveDependency(&IC79SL_full);
589  IC79SL_loglike.reset_and_calculate();
590  IC79_loglike.resolveDependency(&IC79WH_bgloglike);
591  IC79_loglike.resolveDependency(&IC79WH_loglike);
592  IC79_loglike.resolveDependency(&IC79WL_bgloglike);
593  IC79_loglike.resolveDependency(&IC79WL_loglike);
594  IC79_loglike.resolveDependency(&IC79SL_bgloglike);
595  IC79_loglike.resolveDependency(&IC79SL_loglike);
596  IC79_loglike.reset_and_calculate();
597  // Save the result
598  results["IceCube_79_lnL"][current_backend] = IC79_loglike(0);
599 
600  } // End of DarkSUSY 5.1.3 calculations
601 
602 
603 
604 
605  //
606  // ======= Perform all calculations for backend DarkSUSY_MSSM 6.1.1 =======
607  //
608 
609  current_backend = "DarkSUSY_MSSM6.1.1";
610 
611  if (not Backends::backendInfo().works[current_backend])
612  {
613  backends_not_built.push_back(current_backend);
614  }
615  else
616  {
617  // Initialize DarkSUSY 6 MSSM backend
618  DarkSUSY_MSSM_6_1_1_init.notifyOfModel("MSSM30atQ");
619  DarkSUSY_MSSM_6_1_1_init.resolveDependency(&createSpectrum);
620  DarkSUSY_MSSM_6_1_1_init.resolveDependency(&createDecays);
621  if (decays) DarkSUSY_MSSM_6_1_1_init.setOption<bool>("use_dsSLHAread", false);
622  else DarkSUSY_MSSM_6_1_1_init.setOption<bool>("use_dsSLHAread", true);
623  DarkSUSY_MSSM_6_1_1_init.reset_and_calculate();
624 
625  // Initialize DarkSUSY 6 Local Halo Model
628  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmcom);
629  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmisodf);
630  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmframevelcom);
631  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dshmnoclue);
632  DarkSUSY_PointInit_LocalHalo_func.reset_and_calculate();
633 
634  // Relic density calculation with GAMBIT routines and DarkSUSY 6:
635  RD_spectrum_MSSM.resolveDependency(&createDecays);
636  RD_spectrum_MSSM.resolveDependency(&createSpectrum);
637  RD_spectrum_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
638  // Below true if charginos and neutralinos are included in coannihilations:
639  RD_spectrum_MSSM.setOption<bool>("CoannCharginosNeutralinos", true);
640  // Below true if sfermions are included in coannihilations:
641  RD_spectrum_MSSM.setOption<bool>("CoannSfermions", true);
642  // Maximum sparticle mass to be icluded in coannihilations, in units of DM mass:
643  RD_spectrum_MSSM.setOption<double>("CoannMaxMass", 1.6);
644  RD_spectrum_MSSM.reset_and_calculate();
645 
646  RD_spectrum_ordered_func.resolveDependency(&RD_spectrum_MSSM);
647  RD_spectrum_ordered_func.reset_and_calculate();
648 
650  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsancoann);
651  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::DSparticle_code);
652  RD_annrate_DSprep_MSSM_func.reset_and_calculate();
653 
654  RD_eff_annrate_DS_MSSM.notifyOfModel("MSSM30atQ");
656  RD_eff_annrate_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsanwx);
657  RD_eff_annrate_DS_MSSM.reset_and_calculate();
658 
659  RD_oh2_DS_general.resolveDependency(&RD_spectrum_ordered_func);
660  RD_oh2_DS_general.resolveDependency(&RD_eff_annrate_DS_MSSM);
661  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::rdpars);
662  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::rdtime);
663  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsrdcom);
664  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsrdstart);
665  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsrdens);
666  RD_oh2_DS_general.setOption<int>("fast", 1); // 0: normal; 1: fast; 2: dirty
667  RD_oh2_DS_general.reset_and_calculate();
668  // Save the result
669  results["oh2"][current_backend] = RD_oh2_DS_general(0);
670 
671  lnL_oh2_Simple.resolveDependency(&RD_oh2_DS_general);
672  lnL_oh2_Simple.reset_and_calculate();
673  // Save the result
674  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
675 
676 
677  // Set up process catalog based on DarkSUSY annihilation rates
678  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createSpectrum);
679  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createDecays);
680  TH_ProcessCatalog_DS_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
681  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dssigmav0);
682  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dssigmav0tot);
683  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBffdxdy);
684  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBhhdxdy);
685  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBwhdxdy);
686  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsIBwwdxdy);
687  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::IBintvars);
688  TH_ProcessCatalog_DS_MSSM.reset_and_calculate();
689 
690  // Set generic WIMP mass object
691  mwimp_generic.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
692  mwimp_generic.resolveDependency(&DarkMatter_ID_MSSM);
693  mwimp_generic.reset_and_calculate();
694 
695  // Set generic annihilation rate in late universe (v->0 limit)
697  sigmav_late_universe.resolveDependency(&DarkMatter_ID_MSSM);
698  sigmav_late_universe.resolveDependency(&DarkMatterConj_ID_MSSM);
699  sigmav_late_universe.reset_and_calculate();
700  // Save the result
701  results["sigmav0"][current_backend] = sigmav_late_universe(0);
702 
703 
704  // ---- Gamma-ray yields ----
705 
706  // Initialize tabulated gamma-ray yields
707  SimYieldTable_DarkSUSY.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsanyield_sim);
708  SimYieldTable_DarkSUSY.reset_and_calculate();
709 
710  // Identify process as annihilation rather than decay
713  DM_process_from_ProcessCatalog.reset_and_calculate();
714 
715  // Collect missing final states for simulation in cascade MC
717  GA_missingFinalStates.resolveDependency(&SimYieldTable_DarkSUSY);
718  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_MSSM);
719  GA_missingFinalStates.resolveDependency(&DarkMatterConj_ID_MSSM);
721  GA_missingFinalStates.reset_and_calculate();
722 
723  // Infer for which type of final states particles MC should be performed
724  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec<std::string>("gamma"));
725  cascadeMC_FinalStates.reset_and_calculate();
726 
727  // Collect decay information for cascade MC
729  cascadeMC_DecayTable.resolveDependency(&SimYieldTable_DarkSUSY);
730  cascadeMC_DecayTable.reset_and_calculate();
731 
732  // Set up MC loop manager for cascade MC
733  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
734  std::vector<functor*> nested_functions = initVector<functor*>(
736  cascadeMC_LoopManager.setNestedList(nested_functions);
737 
738  // Set up initial state for cascade MC step
739  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
740  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
741 
742  // Perform MC step for cascade MC
743  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
744  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
745  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
746 
747  // Generate histogram for cascade MC
748  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
749  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
750  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
751  cascadeMC_Histograms.resolveDependency(&SimYieldTable_DarkSUSY);
752  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
753  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
754 
755  // Check convergence of cascade MC
756  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
757  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
758 
759  // Start cascade MC loop
760  cascadeMC_LoopManager.reset_and_calculate();
761 
762  // Infer gamma-ray spectra for recorded MC results
763  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
764  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
765  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
766  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
767  cascadeMC_gammaSpectra.reset_and_calculate();
768 
769  // Calculate total gamma-ray yield (cascade MC + tabulated results)
771  GA_AnnYield_General.resolveDependency(&SimYieldTable_DarkSUSY);
772  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_MSSM);
773  GA_AnnYield_General.resolveDependency(&DarkMatterConj_ID_MSSM);
774  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
775  GA_AnnYield_General.reset_and_calculate();
776 
777  // Dump spectrum into file
778  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
779  dump_GammaSpectrum.setOption<std::string>("filename", current_backend + "_" + outname_dNdE_spectrum);
780  dump_GammaSpectrum.reset_and_calculate();
781 
782  // Calculate Fermi LAT dwarf likelihood
784  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
785  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
786  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
787  // Save the result
788  results["FermiLAT_dwarfsph_lnL"][current_backend] = lnL_FermiLATdwarfs_gamLike(0);
789 
790 
791  // ---- Direct detection and IceCube limits ----
792 
793  // Calculate DD couplings with DarkSUSY
795  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::ddcomlegacy);
796  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::ddmssmcom);
797  // The below calculates the DD couplings using the full 1 loop calculation of
798  // Drees Nojiri Phys.Rev. D48 (1993) 3483
799  DD_couplings_DarkSUSY_MSSM.setOption<bool>("loop", true);
800  // Setting the below to false approximates the squark propagator as 1/m_sq^2 to avoid poles.
801  DD_couplings_DarkSUSY_MSSM.setOption<bool>("pole", false);
802  DD_couplings_DarkSUSY_MSSM.reset_and_calculate();
803 
804  // Initialize DDCalc backend
805  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
806  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
807  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
808  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
809  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
810  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
811  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
812  DDCalc_2_2_0_init.reset_and_calculate();
813 
814  // Calculate direct detection rates for LUX 2016
815  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
816  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
817  LUX_2016_Calc.reset_and_calculate();
818 
819  // Calculate direct detection likelihood for LUX 2016
820  LUX_2016_GetLogLikelihood.resolveDependency(&LUX_2016_Calc);
821  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
822  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
823  LUX_2016_GetLogLikelihood.reset_and_calculate();
824  // Save the result
825  results["LUX_2016_lnL"][current_backend] = LUX_2016_GetLogLikelihood(0);
826 
827 
828  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
829  sigma_SI_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
830  sigma_SI_p_simple.reset_and_calculate();
831  // Save the result
832  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
833 
834  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
835  sigma_SD_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
836  sigma_SD_p_simple.reset_and_calculate();
837  // Save the result
838  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
839 
840 
841  // Infer WIMP capture rate in Sun
842  capture_rate_Sun_const_xsec.resolveDependency(&mwimp_generic);
845  capture_rate_Sun_const_xsec.resolveDependency(&RD_fraction_one);
846  capture_rate_Sun_const_xsec.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dssenu_capsuntab);
849  capture_rate_Sun_const_xsec.reset_and_calculate();
850 
851  // Infer WIMP equilibration time in Sun
853  equilibration_time_Sun.resolveDependency(&DarkMatter_ID_MSSM);
855  equilibration_time_Sun.resolveDependency(&mwimp_generic);
857  equilibration_time_Sun.reset_and_calculate();
858 
859  // Infer WIMP annihilation rate in Sun
860  annihilation_rate_Sun.resolveDependency(&equilibration_time_Sun);
862  annihilation_rate_Sun.reset_and_calculate();
863 
864  // Infer neutrino yield from Sun
865  nuyield_from_DS.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
866  nuyield_from_DS.resolveDependency(&mwimp_generic);
867  nuyield_from_DS.resolveDependency(&sigmav_late_universe);
868  nuyield_from_DS.resolveDependency(&DarkMatter_ID_MSSM);
869  nuyield_from_DS.resolveDependency(&DarkMatterConj_ID_MSSM);
870  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::dsgenericwimp_nusetup);
871  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::neutrino_yield);
872  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::DS_neutral_h_decay_channels);
873  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_1_1::Functown::DS_charged_h_decay_channels);
874  nuyield_from_DS.reset_and_calculate();
875 
876 
877  // Calculate number of events at IceCube
878  IC79WH_full.resolveDependency(&mwimp_generic);
879  IC79WH_full.resolveDependency(&annihilation_rate_Sun);
880  IC79WH_full.resolveDependency(&nuyield_from_DS);
881  IC79WH_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
882  IC79WH_full.reset_and_calculate();
883  IC79WL_full.resolveDependency(&mwimp_generic);
884  IC79WL_full.resolveDependency(&annihilation_rate_Sun);
885  IC79WL_full.resolveDependency(&nuyield_from_DS);
886  IC79WL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
887  IC79WL_full.reset_and_calculate();
888  IC79SL_full.resolveDependency(&mwimp_generic);
889  IC79SL_full.resolveDependency(&annihilation_rate_Sun);
890  IC79SL_full.resolveDependency(&nuyield_from_DS);
891  IC79SL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
892  IC79SL_full.reset_and_calculate();
893 
894  // Calculate IceCube likelihood
895  IC79WH_bgloglike.resolveDependency(&IC79WH_full);
896  IC79WH_bgloglike.reset_and_calculate();
897  IC79WH_loglike.resolveDependency(&IC79WH_full);
898  IC79WH_loglike.reset_and_calculate();
899  IC79WL_bgloglike.resolveDependency(&IC79WL_full);
900  IC79WL_bgloglike.reset_and_calculate();
901  IC79WL_loglike.resolveDependency(&IC79WL_full);
902  IC79WL_loglike.reset_and_calculate();
903  IC79SL_bgloglike.resolveDependency(&IC79SL_full);
904  IC79SL_bgloglike.reset_and_calculate();
905  IC79SL_loglike.resolveDependency(&IC79SL_full);
906  IC79SL_loglike.reset_and_calculate();
907  IC79_loglike.resolveDependency(&IC79WH_bgloglike);
908  IC79_loglike.resolveDependency(&IC79WH_loglike);
909  IC79_loglike.resolveDependency(&IC79WL_bgloglike);
910  IC79_loglike.resolveDependency(&IC79WL_loglike);
911  IC79_loglike.resolveDependency(&IC79SL_bgloglike);
912  IC79_loglike.resolveDependency(&IC79SL_loglike);
913  IC79_loglike.reset_and_calculate();
914  // Save the result
915  results["IceCube_79_lnL"][current_backend] = IC79_loglike(0);
916 
917  } // End of DarkSUSY_MSSM 6.1.1 calculations
918 
919 
920 
921 
922  //
923  // ======= Perform all calculations for backend DarkSUSY_MSSM 6.2.2 =======
924  //
925 
926  current_backend = "DarkSUSY_MSSM6.2.2";
927 
928  if (not Backends::backendInfo().works[current_backend])
929  {
930  backends_not_built.push_back(current_backend);
931  }
932  else
933  {
934  // Initialize DarkSUSY 6 MSSM backend
935  DarkSUSY_MSSM_6_2_2_init.notifyOfModel("MSSM30atQ");
936  DarkSUSY_MSSM_6_2_2_init.resolveDependency(&createSpectrum);
937  DarkSUSY_MSSM_6_2_2_init.resolveDependency(&createDecays);
938  if (decays) DarkSUSY_MSSM_6_2_2_init.setOption<bool>("use_dsSLHAread", false);
939  else DarkSUSY_MSSM_6_2_2_init.setOption<bool>("use_dsSLHAread", true);
940  DarkSUSY_MSSM_6_2_2_init.reset_and_calculate();
941 
942  // Initialize DarkSUSY 6 Local Halo Model
945  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmcom);
946  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmisodf);
947  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmframevelcom);
948  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dshmnoclue);
949  DarkSUSY_PointInit_LocalHalo_func.reset_and_calculate();
950 
951  // Relic density calculation with GAMBIT routines and DarkSUSY 6:
952  RD_spectrum_MSSM.resolveDependency(&createDecays);
953  RD_spectrum_MSSM.resolveDependency(&createSpectrum);
954  RD_spectrum_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
955  // Below true if charginos and neutralinos are included in coannihilations:
956  RD_spectrum_MSSM.setOption<bool>("CoannCharginosNeutralinos", true);
957  // Below true if sfermions are included in coannihilations:
958  RD_spectrum_MSSM.setOption<bool>("CoannSfermions", true);
959  // Maximum sparticle mass to be icluded in coannihilations, in units of DM mass:
960  RD_spectrum_MSSM.setOption<double>("CoannMaxMass", 1.6);
961  RD_spectrum_MSSM.reset_and_calculate();
962 
963  RD_spectrum_ordered_func.resolveDependency(&RD_spectrum_MSSM);
964  RD_spectrum_ordered_func.reset_and_calculate();
965 
967  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsancoann);
968  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::DSparticle_code);
969  RD_annrate_DSprep_MSSM_func.reset_and_calculate();
970 
971  RD_eff_annrate_DS_MSSM.notifyOfModel("MSSM30atQ");
973  RD_eff_annrate_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsanwx);
974  RD_eff_annrate_DS_MSSM.reset_and_calculate();
975 
976  RD_oh2_DS_general.resolveDependency(&RD_spectrum_ordered_func);
977  RD_oh2_DS_general.resolveDependency(&RD_eff_annrate_DS_MSSM);
978  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::rdpars);
979  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::rdtime);
980  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsrdcom);
981  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsrdstart);
982  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsrdens);
983  RD_oh2_DS_general.setOption<int>("fast", 1); // 0: normal; 1: fast; 2: dirty
984  RD_oh2_DS_general.reset_and_calculate();
985  // Save the result
986  results["oh2"][current_backend] = RD_oh2_DS_general(0);
987 
988  lnL_oh2_Simple.resolveDependency(&RD_oh2_DS_general);
989  lnL_oh2_Simple.reset_and_calculate();
990  // Save the result
991  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
992 
993 
994  // Set up process catalog based on DarkSUSY annihilation rates
995  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createSpectrum);
996  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createDecays);
997  TH_ProcessCatalog_DS_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
998  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dssigmav0);
999  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dssigmav0tot);
1000  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBffdxdy);
1001  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBhhdxdy);
1002  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBwhdxdy);
1003  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsIBwwdxdy);
1004  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::IBintvars);
1005  TH_ProcessCatalog_DS_MSSM.reset_and_calculate();
1006 
1007  // Set generic WIMP mass object
1008  mwimp_generic.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1009  mwimp_generic.resolveDependency(&DarkMatter_ID_MSSM);
1010  mwimp_generic.reset_and_calculate();
1011 
1012  // Set generic annihilation rate in late universe (v->0 limit)
1013  sigmav_late_universe.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1014  sigmav_late_universe.resolveDependency(&DarkMatter_ID_MSSM);
1015  sigmav_late_universe.resolveDependency(&DarkMatterConj_ID_MSSM);
1016  sigmav_late_universe.reset_and_calculate();
1017  // Save the result
1018  results["sigmav0"][current_backend] = sigmav_late_universe(0);
1019 
1020 
1021  // ---- Gamma-ray yields ----
1022 
1023  // Initialize tabulated gamma-ray yields
1024  SimYieldTable_DarkSUSY.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsanyield_sim);
1025  SimYieldTable_DarkSUSY.reset_and_calculate();
1026 
1027  // Identify process as annihilation rather than decay
1030  DM_process_from_ProcessCatalog.reset_and_calculate();
1031 
1032  // Collect missing final states for simulation in cascade MC
1034  GA_missingFinalStates.resolveDependency(&SimYieldTable_DarkSUSY);
1035  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_MSSM);
1036  GA_missingFinalStates.resolveDependency(&DarkMatterConj_ID_MSSM);
1038  GA_missingFinalStates.reset_and_calculate();
1039 
1040  // Infer for which type of final states particles MC should be performed
1041  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec<std::string>("gamma"));
1042  cascadeMC_FinalStates.reset_and_calculate();
1043 
1044  // Collect decay information for cascade MC
1045  cascadeMC_DecayTable.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1046  cascadeMC_DecayTable.resolveDependency(&SimYieldTable_DarkSUSY);
1047  cascadeMC_DecayTable.reset_and_calculate();
1048 
1049  // Set up MC loop manager for cascade MC
1050  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
1051  std::vector<functor*> nested_functions = initVector<functor*>(
1053  cascadeMC_LoopManager.setNestedList(nested_functions);
1054 
1055  // Set up initial state for cascade MC step
1056  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
1057  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
1058 
1059  // Perform MC step for cascade MC
1060  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
1061  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
1062  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
1063 
1064  // Generate histogram for cascade MC
1065  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
1066  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
1067  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1068  cascadeMC_Histograms.resolveDependency(&SimYieldTable_DarkSUSY);
1069  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
1070  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
1071 
1072  // Check convergence of cascade MC
1073  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
1074  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
1075 
1076  // Start cascade MC loop
1077  cascadeMC_LoopManager.reset_and_calculate();
1078 
1079  // Infer gamma-ray spectra for recorded MC results
1080  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
1081  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
1082  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
1083  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
1084  cascadeMC_gammaSpectra.reset_and_calculate();
1085 
1086  // Calculate total gamma-ray yield (cascade MC + tabulated results)
1087  GA_AnnYield_General.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1088  GA_AnnYield_General.resolveDependency(&SimYieldTable_DarkSUSY);
1089  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_MSSM);
1090  GA_AnnYield_General.resolveDependency(&DarkMatterConj_ID_MSSM);
1091  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
1092  GA_AnnYield_General.reset_and_calculate();
1093 
1094  // Dump spectrum into file
1095  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
1096  dump_GammaSpectrum.setOption<std::string>("filename", current_backend + "_" + outname_dNdE_spectrum);
1097  dump_GammaSpectrum.reset_and_calculate();
1098 
1099  // Calculate Fermi LAT dwarf likelihood
1100  lnL_FermiLATdwarfs_gamLike.resolveDependency(&GA_AnnYield_General);
1101  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
1103  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
1104  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
1105  // Save the result
1106  results["FermiLAT_dwarfsph_lnL"][current_backend] = lnL_FermiLATdwarfs_gamLike(0);
1107 
1108 
1109  // ---- Direct detection and IceCube limits ----
1110 
1111  // Calculate DD couplings with DarkSUSY
1113  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::ddcomlegacy);
1114  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::ddmssmcom);
1115  // The below calculates the DD couplings using the full 1 loop calculation of
1116  // Drees Nojiri Phys.Rev. D48 (1993) 3483
1117  DD_couplings_DarkSUSY_MSSM.setOption<bool>("loop", true);
1118  // Setting the below to false approximates the squark propagator as 1/m_sq^2 to avoid poles.
1119  DD_couplings_DarkSUSY_MSSM.setOption<bool>("pole", false);
1120  DD_couplings_DarkSUSY_MSSM.reset_and_calculate();
1121 
1122  // Initialize DDCalc backend
1123  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
1124  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
1125  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
1126  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
1127  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
1128  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
1129  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1130  DDCalc_2_2_0_init.reset_and_calculate();
1131 
1132  // Calculate direct detection rates for LUX 2016
1133  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
1134  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
1135  LUX_2016_Calc.reset_and_calculate();
1136 
1137  // Calculate direct detection likelihood for LUX 2016
1138  LUX_2016_GetLogLikelihood.resolveDependency(&LUX_2016_Calc);
1139  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
1140  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
1141  LUX_2016_GetLogLikelihood.reset_and_calculate();
1142  // Save the result
1143  results["LUX_2016_lnL"][current_backend] = LUX_2016_GetLogLikelihood(0);
1144 
1145 
1146  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
1147  sigma_SI_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1148  sigma_SI_p_simple.reset_and_calculate();
1149  // Save the result
1150  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
1151 
1152  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
1153  sigma_SD_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1154  sigma_SD_p_simple.reset_and_calculate();
1155  // Save the result
1156  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
1157 
1158 
1159  // Infer WIMP capture rate in Sun
1160  capture_rate_Sun_const_xsec.resolveDependency(&mwimp_generic);
1161  capture_rate_Sun_const_xsec.resolveDependency(&sigma_SI_p_simple);
1162  capture_rate_Sun_const_xsec.resolveDependency(&sigma_SD_p_simple);
1163  capture_rate_Sun_const_xsec.resolveDependency(&RD_fraction_one);
1164  capture_rate_Sun_const_xsec.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dssenu_capsuntab);
1167  capture_rate_Sun_const_xsec.reset_and_calculate();
1168 
1169  // Infer WIMP equilibration time in Sun
1171  equilibration_time_Sun.resolveDependency(&DarkMatter_ID_MSSM);
1172  equilibration_time_Sun.resolveDependency(&DarkMatterConj_ID_MSSM);
1173  equilibration_time_Sun.resolveDependency(&mwimp_generic);
1175  equilibration_time_Sun.reset_and_calculate();
1176 
1177  // Infer WIMP annihilation rate in Sun
1178  annihilation_rate_Sun.resolveDependency(&equilibration_time_Sun);
1180  annihilation_rate_Sun.reset_and_calculate();
1181 
1182  // Infer neutrino yield from Sun
1183  nuyield_from_DS.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1184  nuyield_from_DS.resolveDependency(&mwimp_generic);
1185  nuyield_from_DS.resolveDependency(&sigmav_late_universe);
1186  nuyield_from_DS.resolveDependency(&DarkMatter_ID_MSSM);
1187  nuyield_from_DS.resolveDependency(&DarkMatterConj_ID_MSSM);
1188  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::dsgenericwimp_nusetup);
1189  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::neutrino_yield);
1190  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::DS_neutral_h_decay_channels);
1191  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_2::Functown::DS_charged_h_decay_channels);
1192  nuyield_from_DS.reset_and_calculate();
1193 
1194 
1195  // Calculate number of events at IceCube
1196  IC79WH_full.resolveDependency(&mwimp_generic);
1197  IC79WH_full.resolveDependency(&annihilation_rate_Sun);
1198  IC79WH_full.resolveDependency(&nuyield_from_DS);
1199  IC79WH_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1200  IC79WH_full.reset_and_calculate();
1201  IC79WL_full.resolveDependency(&mwimp_generic);
1202  IC79WL_full.resolveDependency(&annihilation_rate_Sun);
1203  IC79WL_full.resolveDependency(&nuyield_from_DS);
1204  IC79WL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1205  IC79WL_full.reset_and_calculate();
1206  IC79SL_full.resolveDependency(&mwimp_generic);
1207  IC79SL_full.resolveDependency(&annihilation_rate_Sun);
1208  IC79SL_full.resolveDependency(&nuyield_from_DS);
1209  IC79SL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1210  IC79SL_full.reset_and_calculate();
1211 
1212  // Calculate IceCube likelihood
1213  IC79WH_bgloglike.resolveDependency(&IC79WH_full);
1214  IC79WH_bgloglike.reset_and_calculate();
1215  IC79WH_loglike.resolveDependency(&IC79WH_full);
1216  IC79WH_loglike.reset_and_calculate();
1217  IC79WL_bgloglike.resolveDependency(&IC79WL_full);
1218  IC79WL_bgloglike.reset_and_calculate();
1219  IC79WL_loglike.resolveDependency(&IC79WL_full);
1220  IC79WL_loglike.reset_and_calculate();
1221  IC79SL_bgloglike.resolveDependency(&IC79SL_full);
1222  IC79SL_bgloglike.reset_and_calculate();
1223  IC79SL_loglike.resolveDependency(&IC79SL_full);
1224  IC79SL_loglike.reset_and_calculate();
1225  IC79_loglike.resolveDependency(&IC79WH_bgloglike);
1226  IC79_loglike.resolveDependency(&IC79WH_loglike);
1227  IC79_loglike.resolveDependency(&IC79WL_bgloglike);
1228  IC79_loglike.resolveDependency(&IC79WL_loglike);
1229  IC79_loglike.resolveDependency(&IC79SL_bgloglike);
1230  IC79_loglike.resolveDependency(&IC79SL_loglike);
1231  IC79_loglike.reset_and_calculate();
1232  // Save the result
1233  results["IceCube_79_lnL"][current_backend] = IC79_loglike(0);
1234 
1235  } // End of DarkSUSY_MSSM 6.2.2 calculations
1236 
1237 
1238  //
1239  // ======= Perform all calculations for backend DarkSUSY_MSSM 6.2.5 =======
1240  //
1241 
1242  current_backend = "DarkSUSY_MSSM6.2.5";
1243 
1244  if (not Backends::backendInfo().works[current_backend])
1245  {
1246  backends_not_built.push_back(current_backend);
1247  }
1248  else
1249  {
1250  // Initialize DarkSUSY 6 MSSM backend
1251  DarkSUSY_MSSM_6_2_5_init.notifyOfModel("MSSM30atQ");
1252  DarkSUSY_MSSM_6_2_5_init.resolveDependency(&createSpectrum);
1253  DarkSUSY_MSSM_6_2_5_init.resolveDependency(&createDecays);
1254  if (decays) DarkSUSY_MSSM_6_2_5_init.setOption<bool>("use_dsSLHAread", false);
1255  else DarkSUSY_MSSM_6_2_5_init.setOption<bool>("use_dsSLHAread", true);
1256  DarkSUSY_MSSM_6_2_5_init.reset_and_calculate();
1257 
1258  // Initialize DarkSUSY 6 Local Halo Model
1261  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dshmcom);
1262  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dshmisodf);
1263  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dshmframevelcom);
1264  DarkSUSY_PointInit_LocalHalo_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dshmnoclue);
1265  DarkSUSY_PointInit_LocalHalo_func.reset_and_calculate();
1266 
1267  // Relic density calculation with GAMBIT routines and DarkSUSY 6:
1268  RD_spectrum_MSSM.resolveDependency(&createDecays);
1269  RD_spectrum_MSSM.resolveDependency(&createSpectrum);
1270  RD_spectrum_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
1271  // Below true if charginos and neutralinos are included in coannihilations:
1272  RD_spectrum_MSSM.setOption<bool>("CoannCharginosNeutralinos", true);
1273  // Below true if sfermions are included in coannihilations:
1274  RD_spectrum_MSSM.setOption<bool>("CoannSfermions", true);
1275  // Maximum sparticle mass to be icluded in coannihilations, in units of DM mass:
1276  RD_spectrum_MSSM.setOption<double>("CoannMaxMass", 1.6);
1277  RD_spectrum_MSSM.reset_and_calculate();
1278 
1279  RD_spectrum_ordered_func.resolveDependency(&RD_spectrum_MSSM);
1280  RD_spectrum_ordered_func.reset_and_calculate();
1281 
1283  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsancoann);
1284  RD_annrate_DSprep_MSSM_func.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::DSparticle_code);
1285  RD_annrate_DSprep_MSSM_func.reset_and_calculate();
1286 
1287  RD_eff_annrate_DS_MSSM.notifyOfModel("MSSM30atQ");
1289  RD_eff_annrate_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsanwx);
1290  RD_eff_annrate_DS_MSSM.reset_and_calculate();
1291 
1292  RD_oh2_DS_general.resolveDependency(&RD_spectrum_ordered_func);
1293  RD_oh2_DS_general.resolveDependency(&RD_eff_annrate_DS_MSSM);
1294  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::rdpars);
1295  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::rdtime);
1296  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsrdcom);
1297  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsrdstart);
1298  RD_oh2_DS_general.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsrdens);
1299  RD_oh2_DS_general.setOption<int>("fast", 1); // 0: normal; 1: fast; 2: dirty
1300  RD_oh2_DS_general.reset_and_calculate();
1301  // Save the result
1302  results["oh2"][current_backend] = RD_oh2_DS_general(0);
1303 
1304  lnL_oh2_Simple.resolveDependency(&RD_oh2_DS_general);
1305  lnL_oh2_Simple.reset_and_calculate();
1306  // Save the result
1307  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
1308 
1309 
1310  // Set up process catalog based on DarkSUSY annihilation rates
1311  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createSpectrum);
1312  TH_ProcessCatalog_DS_MSSM.resolveDependency(&createDecays);
1313  TH_ProcessCatalog_DS_MSSM.resolveDependency(&DarkMatter_ID_MSSM);
1314  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dssigmav0);
1315  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dssigmav0tot);
1316  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsIBffdxdy);
1317  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsIBhhdxdy);
1318  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsIBwhdxdy);
1319  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsIBwwdxdy);
1320  TH_ProcessCatalog_DS_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::IBintvars);
1321  TH_ProcessCatalog_DS_MSSM.reset_and_calculate();
1322 
1323  // Set generic WIMP mass object
1324  mwimp_generic.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1325  mwimp_generic.resolveDependency(&DarkMatter_ID_MSSM);
1326  mwimp_generic.reset_and_calculate();
1327 
1328  // Set generic annihilation rate in late universe (v->0 limit)
1329  sigmav_late_universe.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1330  sigmav_late_universe.resolveDependency(&DarkMatter_ID_MSSM);
1331  sigmav_late_universe.reset_and_calculate();
1332  // Save the result
1333  results["sigmav0"][current_backend] = sigmav_late_universe(0);
1334 
1335 
1336  // ---- Gamma-ray yields ----
1337 
1338  // Initialize tabulated gamma-ray yields
1339  SimYieldTable_DarkSUSY.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsanyield_sim);
1340  SimYieldTable_DarkSUSY.reset_and_calculate();
1341 
1342  // Collect missing final states for simulation in cascade MC
1344  GA_missingFinalStates.resolveDependency(&SimYieldTable_DarkSUSY);
1345  GA_missingFinalStates.resolveDependency(&DarkMatter_ID_MSSM);
1346  GA_missingFinalStates.reset_and_calculate();
1347 
1348  // Infer for which type of final states particles MC should be performed
1349  cascadeMC_FinalStates.setOption<std::vector<std::string>>("cMC_finalStates", daFunk::vec<std::string>("gamma"));
1350  cascadeMC_FinalStates.reset_and_calculate();
1351 
1352  // Collect decay information for cascade MC
1353  cascadeMC_DecayTable.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1354  cascadeMC_DecayTable.resolveDependency(&SimYieldTable_DarkSUSY);
1355  cascadeMC_DecayTable.reset_and_calculate();
1356 
1357  // Set up MC loop manager for cascade MC
1358  cascadeMC_LoopManager.resolveDependency(&GA_missingFinalStates);
1359  std::vector<functor*> nested_functions = initVector<functor*>(
1361  cascadeMC_LoopManager.setNestedList(nested_functions);
1362 
1363  // Set up initial state for cascade MC step
1364  cascadeMC_InitialState.resolveDependency(&GA_missingFinalStates);
1365  cascadeMC_InitialState.resolveLoopManager(&cascadeMC_LoopManager);
1366 
1367  // Perform MC step for cascade MC
1368  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_InitialState);
1369  cascadeMC_GenerateChain.resolveDependency(&cascadeMC_DecayTable);
1370  cascadeMC_GenerateChain.resolveLoopManager(&cascadeMC_LoopManager);
1371 
1372  // Generate histogram for cascade MC
1373  cascadeMC_Histograms.resolveDependency(&cascadeMC_InitialState);
1374  cascadeMC_Histograms.resolveDependency(&cascadeMC_GenerateChain);
1375  cascadeMC_Histograms.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1376  cascadeMC_Histograms.resolveDependency(&SimYieldTable_DarkSUSY);
1377  cascadeMC_Histograms.resolveDependency(&cascadeMC_FinalStates);
1378  cascadeMC_Histograms.resolveLoopManager(&cascadeMC_LoopManager);
1379 
1380  // Check convergence of cascade MC
1381  cascadeMC_EventCount.resolveDependency(&cascadeMC_InitialState);
1382  cascadeMC_EventCount.resolveLoopManager(&cascadeMC_LoopManager);
1383 
1384  // Start cascade MC loop
1385  cascadeMC_LoopManager.reset_and_calculate();
1386 
1387  // Infer gamma-ray spectra for recorded MC results
1388  cascadeMC_gammaSpectra.resolveDependency(&GA_missingFinalStates);
1389  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_FinalStates);
1390  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_Histograms);
1391  cascadeMC_gammaSpectra.resolveDependency(&cascadeMC_EventCount);
1392  cascadeMC_gammaSpectra.reset_and_calculate();
1393 
1394  // Calculate total gamma-ray yield (cascade MC + tabulated results)
1395  GA_AnnYield_General.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1396  GA_AnnYield_General.resolveDependency(&SimYieldTable_DarkSUSY);
1397  GA_AnnYield_General.resolveDependency(&DarkMatter_ID_MSSM);
1398  GA_AnnYield_General.resolveDependency(&cascadeMC_gammaSpectra);
1399  GA_AnnYield_General.reset_and_calculate();
1400 
1401  // Dump spectrum into file
1402  dump_GammaSpectrum.resolveDependency(&GA_AnnYield_General);
1403  dump_GammaSpectrum.setOption<std::string>("filename", current_backend + "_" + outname_dNdE_spectrum);
1404  dump_GammaSpectrum.reset_and_calculate();
1405 
1406  // Calculate Fermi LAT dwarf likelihood
1407  lnL_FermiLATdwarfs_gamLike.resolveDependency(&GA_AnnYield_General);
1408  lnL_FermiLATdwarfs_gamLike.resolveDependency(&RD_fraction_one);
1409  lnL_FermiLATdwarfs_gamLike.resolveBackendReq(&Backends::gamLike_1_0_1::Functown::lnL);
1410  lnL_FermiLATdwarfs_gamLike.reset_and_calculate();
1411  // Save the result
1412  results["FermiLAT_dwarfsph_lnL"][current_backend] = lnL_FermiLATdwarfs_gamLike(0);
1413 
1414 
1415  // ---- Direct detection and IceCube limits ----
1416 
1417  // Calculate DD couplings with DarkSUSY
1419  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::ddcomlegacy);
1420  DD_couplings_DarkSUSY_MSSM.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::ddmssmcom);
1421  // The below calculates the DD couplings using the full 1 loop calculation of
1422  // Drees Nojiri Phys.Rev. D48 (1993) 3483
1423  DD_couplings_DarkSUSY_MSSM.setOption<bool>("loop", true);
1424  // Setting the below to false approximates the squark propagator as 1/m_sq^2 to avoid poles.
1425  DD_couplings_DarkSUSY_MSSM.setOption<bool>("pole", false);
1426  DD_couplings_DarkSUSY_MSSM.reset_and_calculate();
1427 
1428  // Initialize DDCalc backend
1429  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
1430  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
1431  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
1432  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
1433  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
1434  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
1435  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1436  DDCalc_2_2_0_init.reset_and_calculate();
1437 
1438  // Calculate direct detection rates for LUX 2016
1439  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
1440  LUX_2016_Calc.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple);
1441  LUX_2016_Calc.reset_and_calculate();
1442 
1443  // Calculate direct detection likelihood for LUX 2016
1444  LUX_2016_GetLogLikelihood.resolveDependency(&LUX_2016_Calc);
1445  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment);
1446  LUX_2016_GetLogLikelihood.resolveBackendReq(&Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood);
1447  LUX_2016_GetLogLikelihood.reset_and_calculate();
1448  // Save the result
1449  results["LUX_2016_lnL"][current_backend] = LUX_2016_GetLogLikelihood(0);
1450 
1451 
1452  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
1453  sigma_SI_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1454  sigma_SI_p_simple.reset_and_calculate();
1455  // Save the result
1456  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
1457 
1458  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
1459  sigma_SD_p_simple.resolveDependency(&DD_couplings_DarkSUSY_MSSM);
1460  sigma_SD_p_simple.reset_and_calculate();
1461  // Save the result
1462  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
1463 
1464 
1465  // Infer WIMP capture rate in Sun
1466  capture_rate_Sun_const_xsec.resolveDependency(&mwimp_generic);
1467  capture_rate_Sun_const_xsec.resolveDependency(&sigma_SI_p_simple);
1468  capture_rate_Sun_const_xsec.resolveDependency(&sigma_SD_p_simple);
1469  capture_rate_Sun_const_xsec.resolveDependency(&RD_fraction_one);
1470  capture_rate_Sun_const_xsec.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dssenu_capsuntab);
1473  capture_rate_Sun_const_xsec.reset_and_calculate();
1474 
1475  // Infer WIMP equilibration time in Sun
1477  equilibration_time_Sun.resolveDependency(&DarkMatter_ID_MSSM);
1478  equilibration_time_Sun.resolveDependency(&mwimp_generic);
1480  equilibration_time_Sun.reset_and_calculate();
1481 
1482  // Infer WIMP annihilation rate in Sun
1483  annihilation_rate_Sun.resolveDependency(&equilibration_time_Sun);
1485  annihilation_rate_Sun.reset_and_calculate();
1486 
1487  // Infer neutrino yield from Sun
1488  nuyield_from_DS.resolveDependency(&TH_ProcessCatalog_DS_MSSM);
1489  nuyield_from_DS.resolveDependency(&mwimp_generic);
1490  nuyield_from_DS.resolveDependency(&sigmav_late_universe);
1491  nuyield_from_DS.resolveDependency(&DarkMatter_ID_MSSM);
1492  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::dsgenericwimp_nusetup);
1493  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::neutrino_yield);
1494  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::DS_neutral_h_decay_channels);
1495  nuyield_from_DS.resolveBackendReq(&Backends::DarkSUSY_MSSM_6_2_5::Functown::DS_charged_h_decay_channels);
1496  nuyield_from_DS.reset_and_calculate();
1497 
1498 
1499  // Calculate number of events at IceCube
1500  IC79WH_full.resolveDependency(&mwimp_generic);
1501  IC79WH_full.resolveDependency(&annihilation_rate_Sun);
1502  IC79WH_full.resolveDependency(&nuyield_from_DS);
1503  IC79WH_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1504  IC79WH_full.reset_and_calculate();
1505  IC79WL_full.resolveDependency(&mwimp_generic);
1506  IC79WL_full.resolveDependency(&annihilation_rate_Sun);
1507  IC79WL_full.resolveDependency(&nuyield_from_DS);
1508  IC79WL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1509  IC79WL_full.reset_and_calculate();
1510  IC79SL_full.resolveDependency(&mwimp_generic);
1511  IC79SL_full.resolveDependency(&annihilation_rate_Sun);
1512  IC79SL_full.resolveDependency(&nuyield_from_DS);
1513  IC79SL_full.resolveBackendReq(&Backends::nulike_1_0_9::Functown::nulike_bounds);
1514  IC79SL_full.reset_and_calculate();
1515 
1516  // Calculate IceCube likelihood
1517  IC79WH_bgloglike.resolveDependency(&IC79WH_full);
1518  IC79WH_bgloglike.reset_and_calculate();
1519  IC79WH_loglike.resolveDependency(&IC79WH_full);
1520  IC79WH_loglike.reset_and_calculate();
1521  IC79WL_bgloglike.resolveDependency(&IC79WL_full);
1522  IC79WL_bgloglike.reset_and_calculate();
1523  IC79WL_loglike.resolveDependency(&IC79WL_full);
1524  IC79WL_loglike.reset_and_calculate();
1525  IC79SL_bgloglike.resolveDependency(&IC79SL_full);
1526  IC79SL_bgloglike.reset_and_calculate();
1527  IC79SL_loglike.resolveDependency(&IC79SL_full);
1528  IC79SL_loglike.reset_and_calculate();
1529  IC79_loglike.resolveDependency(&IC79WH_bgloglike);
1530  IC79_loglike.resolveDependency(&IC79WH_loglike);
1531  IC79_loglike.resolveDependency(&IC79WL_bgloglike);
1532  IC79_loglike.resolveDependency(&IC79WL_loglike);
1533  IC79_loglike.resolveDependency(&IC79SL_bgloglike);
1534  IC79_loglike.resolveDependency(&IC79SL_loglike);
1535  IC79_loglike.reset_and_calculate();
1536  // Save the result
1537  results["IceCube_79_lnL"][current_backend] = IC79_loglike(0);
1538 
1539  } // End of DarkSUSY_MSSM 6.2.5 calculations
1540 
1541 
1542 
1543 
1544  //
1545  // ======= Perform all calculations for backend MicrOmegas_MSSM 3.6.9.2 =======
1546  //
1547 
1548  current_backend = "MicrOmegas_MSSM3.6.9.2";
1549 
1550  if (not Backends::backendInfo().works[current_backend])
1551  {
1552  backends_not_built.push_back(current_backend);
1553  }
1554  else
1555  {
1556  // Initialize MicrOmegas backend
1557  MicrOmegas_MSSM_3_6_9_2_init.notifyOfModel("MSSM30atQ");
1558  MicrOmegas_MSSM_3_6_9_2_init.resolveDependency(&createSpectrum);
1559  MicrOmegas_MSSM_3_6_9_2_init.resolveDependency(&createDecays);
1560  MicrOmegas_MSSM_3_6_9_2_init.resolveDependency(&createSLHA1Names);
1561  // Use decay table if it is present:
1562  if (decays) MicrOmegas_MSSM_3_6_9_2_init.setOption<bool>("internal_decays", false);
1563  else MicrOmegas_MSSM_3_6_9_2_init.setOption<bool>("internal_decays", true);
1564  MicrOmegas_MSSM_3_6_9_2_init.reset_and_calculate();
1565  // For the below VXdecay = 0 - no 3 body final states via virtual X
1566  // 1 - annihilations to 3 body final states via virtual X
1567  // 2 - (co)annihilations to 3 body final states via virtual X
1568  MicrOmegas_MSSM_3_6_9_2_init.setOption<int>("VZdecay", 0);
1569  MicrOmegas_MSSM_3_6_9_2_init.setOption<int>("VWdecay", 0);
1570  MicrOmegas_MSSM_3_6_9_2_init.reset_and_calculate();
1571 
1572  // Relic density calculation with MicrOmegas
1573  RD_oh2_Xf_MicrOmegas.notifyOfModel("MSSM30atQ");
1574  RD_oh2_Xf_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::darkOmega);
1575  RD_oh2_Xf_MicrOmegas.setOption<int>("fast", 1); // 0: accurate; 1: fast
1576  RD_oh2_Xf_MicrOmegas.setOption<double>("Beps", 1e-5); // Beps=1e-5 recommended, Beps=1 switches coannihilation off
1577  RD_oh2_Xf_MicrOmegas.reset_and_calculate();
1578  RD_oh2_MicrOmegas.resolveDependency(&RD_oh2_Xf_MicrOmegas);
1579  RD_oh2_MicrOmegas.reset_and_calculate();
1580  // Save the result
1581  results["oh2"][current_backend] = RD_oh2_MicrOmegas(0);
1582 
1583  lnL_oh2_Simple.resolveDependency(&RD_oh2_MicrOmegas);
1584  lnL_oh2_Simple.reset_and_calculate();
1585  // Save the result
1586  results["oh2_lnL"][current_backend] = lnL_oh2_Simple(0);
1587 
1588  // <sigma v> (v->0 limit) self-annihilation calculation with MicrOmegas:
1589  sigmav_late_universe_MicrOmegas.notifyOfModel("MSSM30atQ");
1590  sigmav_late_universe_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::calcSpectrum);
1591  sigmav_late_universe_MicrOmegas.reset_and_calculate();
1592  results["sigmav0"][current_backend] = sigmav_late_universe_MicrOmegas(0);
1593 
1594  // Direct detection calculations with Micromegas
1595 
1596  // Calculate DD couplings with Micromegas
1597  // DD_couplings_MicrOmegas.notifyOfModel("MSSM30atQ");
1598  // DD_couplings_MicrOmegas.notifyOfModel("nuclear_params_fnq");
1599  // DD_couplings_MicrOmegas.resolveDependency(&Models::nuclear_params_fnq::Functown::primary_parameters);
1600  DD_couplings_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::nucleonAmplitudes);
1601  DD_couplings_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::FeScLoop);
1602  DD_couplings_MicrOmegas.resolveBackendReq(&Backends::MicrOmegas_MSSM_3_6_9_2::Functown::mocommon_);
1603  // The below includes neutralino-gluon scattering via a box diagram
1604  DD_couplings_MicrOmegas.setOption<bool>("box", true);
1605  DD_couplings_MicrOmegas.reset_and_calculate();
1606 
1607  // Initialize DDCalc backend
1608  Backends::DDCalc_2_2_0::Functown::DDCalc_CalcRates_simple.setStatus(2);
1609  Backends::DDCalc_2_2_0::Functown::DDCalc_Experiment.setStatus(2);
1610  Backends::DDCalc_2_2_0::Functown::DDCalc_LogLikelihood.setStatus(2);
1611  DDCalc_2_2_0_init.resolveDependency(&ExtractLocalMaxwellianHalo);
1612  DDCalc_2_2_0_init.resolveDependency(&RD_fraction_one);
1613  DDCalc_2_2_0_init.resolveDependency(&mwimp_generic);
1614  DDCalc_2_2_0_init.resolveDependency(&DD_couplings_MicrOmegas);
1615  DDCalc_2_2_0_init.reset_and_calculate();
1616 
1617 
1618  sigma_SI_p_simple.resolveDependency(&mwimp_generic);
1619  sigma_SI_p_simple.resolveDependency(&DD_couplings_MicrOmegas);
1620  sigma_SI_p_simple.reset_and_calculate();
1621  // Save the result
1622  results["sigma_SI_p"][current_backend] = sigma_SI_p_simple(0);
1623 
1624  sigma_SD_p_simple.resolveDependency(&mwimp_generic);
1625  sigma_SD_p_simple.resolveDependency(&DD_couplings_MicrOmegas);
1626  sigma_SD_p_simple.reset_and_calculate();
1627  // Save the result
1628  results["sigma_SD_p"][current_backend] = sigma_SD_p_simple(0);
1629 
1630  } // End of MicrOmegas_MSSM 3.6.9.2 calculations
1631 
1632 
1633 
1634  //
1635  // ======= Construct the output string =======
1636  //
1637 
1638  std::stringstream results_ss;
1639 
1640  for(const std::string& result_key : result_output_order)
1641  {
1642  const std::map<std::string,double>& backends_result_map = results.at(result_key);
1643  results_ss << result_key;
1644  if (results_units.at(result_key) != "") { results_ss << " [" << results_units.at(result_key) << "]"; }
1645  results_ss << " :" << endl;
1646 
1647  for(const auto& kv : backends_result_map)
1648  {
1649  const std::string& backendname = kv.first;
1650  const double& result = kv.second;
1651  results_ss << " " << backendname << ": " << result << " " << results_units.at(result_key) << endl;
1652  }
1653  results_ss << endl;
1654  }
1655 
1656 
1657  //
1658  // ======= Output the result string to screen =======
1659  //
1660 
1661  cout << endl;
1662  cout << "==== RESULTS ====" << endl;
1663  cout << endl;
1664  cout << results_ss.str();
1665  cout << endl;
1666 
1667  // Let the user know what they are missing...
1668  if (backends_not_built.size() > 0)
1669  {
1670  cout << endl;
1671  cout << "NOTE: The following backend(s) are not present:" << endl;
1672  for (const std::string& backend_name : backends_not_built)
1673  {
1674  cout << " - " << backend_name << endl;
1675  }
1676  cout << "If you want results from these backends you need to build them first." << endl;
1677  cout << endl;
1678  }
1679 
1680 
1681  //
1682  // ======= Output the result string to file =======
1683  //
1684 
1685  std::fstream file;
1686  file.open(outname_data, std::ios_base::out);
1687  file << results_ss.str();
1688  file.close();
1689 
1690  }
1691 
1692  catch (std::exception& e)
1693  {
1694  std::cout << "DarkBit_standalone_MSSM has exited with fatal exception: " << e.what() << std::endl;
1695  }
1696 
1697  return 0;
1698 
1699 }
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 DarkMatterConj_ID_MSSM(std::string &result)
Definition: MSSM.cpp:769
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:507
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:640
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.
void RD_eff_annrate_DS5_MSSM(double(*&result)(double &))
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:420
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:100
EXPORT_SYMBOLS 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:93
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:177
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:742
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 DM_process_from_ProcessCatalog(std::string &result)
Information about the nature of the DM process in question (i.e.
Definition: DarkBit.cpp:114
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:165
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