gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
RelicDensity.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
20 
21 #include <chrono>
22 
27 
28 
29 namespace Gambit
30 {
31  namespace DarkBit
32  {
33 
34 //#define DARKBIT_DEBUG
35 //#define DARKBIT_RD_DEBUG
36 
38  //
39  // DarkSUSY Relic density routines
40  //
42 
43 
48  {
49  using namespace Pipes::RD_spectrum_MSSM;
50 
51  std::string DMid = *Dep::DarkMatter_ID;
52  if ( DMid != "~chi0_1" )
53  {
55  "RD_spectrum_MSSM requires DMid to be ~chi0_1.");
56  }
57  // Neutralino DM is self-conjugate
58  result.isSelfConj = true;
59  // This function reports particle IDs in terms of PDG codes
60  result.particle_index_type = "PDG";
61 
62  // Import based on Spectrum objects
63  const Spectrum& matched_spectra = *Dep::MSSM_spectrum;
64  const SubSpectrum& spec = matched_spectra.get_HE();
65  const SubSpectrum& SMspec = matched_spectra.get_LE();
66  // Import based on decay table from DecayBit
67  const DecayTable* myDecays = &(*Dep::decay_rates);
68 
69  result.coannihilatingParticles.clear();
70  result.resonances.clear();
71  result.threshold_energy.clear();
72 
75  bool CoannCharginosNeutralinos = runOptions->getValueOrDef<bool>(true,
76  "CoannCharginosNeutralinos");
77 
80  bool CoannSfermions = runOptions->getValueOrDef<bool>(true,
81  "CoannSfermions");
82 
85  double CoannMaxMass = runOptions->getValueOrDef<double>(1.6,
86  "CoannMaxMass");
87 
88  // first add neutralino=WIMP=least massive 'coannihilating particle'
89  // Note: translation from PDG codes assumes context integer = 0 here and
90  // in the following (essentially "mass ordering"). For consistency the same
91  // has to be done when retrieving information from the RD_spectrum_type result
92  int ContInt = 0;
93  int PDGwimp = 1000022;
94  double mWIMP = std::abs(spec.get(Par::Pole_Mass,Models::ParticleDB().long_name(PDGwimp,ContInt)));
95  result.coannihilatingParticles.push_back(RD_coannihilating_particle(PDGwimp, 2, mWIMP));
96 
97  #ifdef DARKBIT_DEBUG
98  std::cout << "WIMP mass : "<< mWIMP << std::endl;
99  #endif
100 
101  // add all sparticles that are not too heavy
102  double msp;
103  auto addCoannParticle = [&](int pdg0, int dof)
104  {
105  msp = std::abs(spec.get(Par::Pole_Mass,Models::ParticleDB().long_name(pdg0,ContInt)));
106  if (msp/mWIMP < CoannMaxMass)
107  {
108  result.coannihilatingParticles.push_back(RD_coannihilating_particle(pdg0, dof, msp));
109  }
110  };
111 
112  if(CoannCharginosNeutralinos) // include neutralino & chargino coannihilation
113  {
114  addCoannParticle(1000023, 2); // "~chi0_2"
115  addCoannParticle(1000025, 2); // "~chi0_3"
116  addCoannParticle(1000035, 2); // "~chi0_4"
117  addCoannParticle(1000024, 4); // "~chi+_1"
118  addCoannParticle(1000037, 4); // "~chi+_2"
119  }
120  if(CoannSfermions) // include sfermion coannihilation
121  {
122  addCoannParticle(1000011, 2); // "~e-_1"
123  addCoannParticle(1000013, 2); // "~e-_2"
124  addCoannParticle(1000015, 2); // "~e-_3"
125  addCoannParticle(2000011, 2); // "~e-_4"
126  addCoannParticle(2000013, 2); // "~e-_5"
127  addCoannParticle(2000015, 2); // "~e-_6"
128  addCoannParticle(1000012, 1); // "~nu_1"
129  addCoannParticle(1000014, 1); // "~nu_2"
130  addCoannParticle(1000016, 1); // "~nu_3"
131  addCoannParticle(1000001, 6); // "~d_1"
132  addCoannParticle(1000003, 6); // "~d_2"
133  addCoannParticle(1000005, 6); // "~d_3"
134  addCoannParticle(2000001, 6); // "~d_4"
135  addCoannParticle(2000003, 6); // "~d_5"
136  addCoannParticle(2000005, 6); // "~d_6"
137  addCoannParticle(1000002, 6); // "~u_1"
138  addCoannParticle(1000004, 6); // "~u_2"
139  addCoannParticle(1000006, 6); // "~u_3"
140  addCoannParticle(2000002, 6); // "~u_4"
141  addCoannParticle(2000004, 6); // "~u_5"
142  addCoannParticle(2000006, 6); // "~u_6"
143  }
144 
145  // determine resonances for LSP annihilation
146  std::string SMreslist[] = {"Z0","W+"};
147  std::string reslist[] = {"h0_2","h0_1","A0","H+"};
148  int resSMmax=sizeof(SMreslist) / sizeof(SMreslist[0]);
149  int resmax=sizeof(reslist) / sizeof(reslist[0]);
150  // Charged resonances (last items in the lists) can only appear for coannihilations
151  if (result.coannihilatingParticles.size() == 1)
152  {
153  resSMmax -= 1;
154  resmax -= 1;
155  }
156  double mres,Gammares;
157  for (int i=0; i<resSMmax; i++)
158  {
159  mres = std::abs(SMspec.get(Par::Pole_Mass,SMreslist[i]));
160  Gammares = myDecays->at(SMreslist[i]).width_in_GeV;
161  result.resonances.push_back(TH_Resonance(mres,Gammares));
162  }
163  for (int i=0; i<resmax; i++)
164  {
165  mres = std::abs(spec.get(Par::Pole_Mass,reslist[i]));
166  Gammares = myDecays->at(reslist[i]).width_in_GeV;
167  result.resonances.push_back(TH_Resonance(mres,Gammares));
168  }
169 
170  // determine thresholds (coannihilation thresholds will be added later);
171  // lowest threshold = 2*WIMP rest mass (unlike DS convention!)
172  result.threshold_energy.push_back(2*mWIMP);
173  std::string thrlist[] = {"W+","Z0","t"};
174  int thrmax=sizeof(thrlist) / sizeof(thrlist[0]);
175  double mthr;
176  for (int i=0; i<thrmax; i++)
177  {
178  mthr = std::abs(SMspec.get(Par::Pole_Mass,thrlist[i]));
179  if (mthr > mWIMP)
180  {
181  result.threshold_energy.push_back(2*mthr);
182  }
183  }
184 
185  } // function RD_spectrum_MSSM
186 
187 
188 
193  {
194  using namespace Pipes::RD_spectrum_SUSY_DS5;
195 
196  std::vector<int> colist; //potential coannihilating particles (indices)
197  colist.clear();
198 
199  // Neutralino DM is self-conjugate
200  result.isSelfConj = true;
201  // This function reports particle IDs in terms of internal DarkSUSY codes
202  result.particle_index_type = "DarkSUSY";
203 
204  result.coannihilatingParticles.clear();
205  result.resonances.clear();
206  result.threshold_energy.clear();
207 
210  bool CoannCharginosNeutralinos = runOptions->getValueOrDef<bool>(true,
211  "CoannCharginosNeutralinos");
212 
215  bool CoannSfermions = runOptions->getValueOrDef<bool>(true,
216  "CoannSfermions");
217 
220  double CoannMaxMass = runOptions->getValueOrDef<double>(1.6,
221  "CoannMaxMass");
222 
223  // introduce pointers to DS mass spectrum and relevant particle info
224  DS5_PACODES *DSpart = BEreq::pacodes.pointer();
225  DS5_MSPCTM *mymspctm= BEreq::mspctm.pointer();
226  DS_INTDOF *myintdof= BEreq::intdof.pointer();
227 
228  // first add neutralino=WIMP=least massive 'coannihilating particle'
229  result.coannihilatingParticles.push_back(
230  RD_coannihilating_particle(DSpart->kn(1),
231  myintdof->kdof(DSpart->kn(1)),mymspctm->mass(DSpart->kn(1))));
232 
233  #ifdef DARKBIT_DEBUG
234  std::cout << "WIMP : "<< DSpart->kn(1) << " " <<
235  myintdof->kdof(DSpart->kn(1)) << " " << mymspctm->mass(DSpart->kn(1))
236  << std::endl;
237  #endif
238 
239  // include neutralino & chargino coannihilation
240  if(CoannCharginosNeutralinos)
241  {
242  for (int i=2; i<=4; i++)
243  colist.push_back(DSpart->kn(i));
244  colist.push_back(DSpart->kcha(1));
245  colist.push_back(DSpart->kcha(2));
246  }
247 
248  // include sfermion coannihilation
249  if(CoannSfermions)
250  {
251  for (int i=1; i<=6; i++)
252  colist.push_back(DSpart->ksl(i));
253  for (int i=1; i<=3; i++)
254  colist.push_back(DSpart->ksnu(i));
255  for (int i=1; i<=6; i++)
256  colist.push_back(DSpart->ksqu(i));
257  for (int i=1; i<=6; i++)
258  colist.push_back(DSpart->ksqd(i));
259  }
260 
261  // only keep sparticles that are not too heavy
262  for (size_t i=0; i<colist.size(); i++)
263  if (mymspctm->mass(colist[i])/mymspctm->mass(DSpart->kn(1))
264  < CoannMaxMass)
265  result.coannihilatingParticles.push_back(
266  RD_coannihilating_particle(colist[i], myintdof->kdof(colist[i]),
267  mymspctm->mass(colist[i])));
268 
269 
270  // determine resonances for LSP annihilation
271  int reslist[] = {BEreq::DS5particle_code("Z0"),
272  BEreq::DS5particle_code("h0_2"),
273  BEreq::DS5particle_code("h0_1"),
274  BEreq::DS5particle_code("A0"),
275  BEreq::DS5particle_code("W+"),
276  BEreq::DS5particle_code("H+")};
277  int resmax=sizeof(reslist) / sizeof(reslist[0]);
278  // the last 2 resonances in the list can only appear for coannihilations
279  if (result.coannihilatingParticles.size() == 1)
280  resmax -= 2;
281  // (Turns out resonances are never returned with DS5)
282 
283  // determine thresholds; lowest threshold = 2*WIMP rest mass (unlike DS
284  // convention!)
285  result.threshold_energy.push_back(
286  2*result.coannihilatingParticles[0].mass);
287  int thrlist[] = {BEreq::DS5particle_code("W+"),
288  BEreq::DS5particle_code("Z0"),
289  BEreq::DS5particle_code("u_3")};
290  int thrmax=sizeof(thrlist) / sizeof(thrlist[0]);
291  for (int i=0; i<thrmax; i++)
292  if (mymspctm->mass(thrlist[i])>result.coannihilatingParticles[0].mass)
293  result.threshold_energy.push_back(2*mymspctm->mass(thrlist[i]));
294 
295  } // function RD_spectrum_SUSY_DS5
296 
297 
303  {
305 
306  // retrieve annihilation processes and DM properties
307  std::string DMid= *Dep::DarkMatter_ID;
308  TH_Process annihilation =
309  (*Dep::TH_ProcessCatalog).getProcess(DMid, DMid);
310  TH_ParticleProperty DMproperty =
311  (*Dep::TH_ProcessCatalog).getParticleProperty(DMid);
312 
313  // Is DM self-conjugate ?
314  result.isSelfConj = annihilation.isSelfConj;
315  // This function reports particle IDs in terms of internal DarkSUSY codes
316  // NB: This should eventually be changed to PDG, which however requires updating
317  // all existing examples where the invariant rate (entering in the relic density)
318  // is calculated directly from the ProcessCatalogue!
319  result.particle_index_type = "DarkSUSY";
320 
321 
322  // get thresholds & resonances from process catalog
323  result.resonances = annihilation.resonances_thresholds.resonances;
325 
326  result.coannihilatingParticles.clear();
327  // add WIMP=least massive 'coannihilating particle'
328  // NB: particle code (1st entry) is irrelevant (unless Weff is obtained from DS)
329  result.coannihilatingParticles.push_back(
330  RD_coannihilating_particle(100,1+DMproperty.spin2,DMproperty.mass));
331 
332  #ifdef DARKBIT_DEBUG
333  std::cout << "DM dof = " << 1+ DMproperty.spin2 << std::endl;
334  // std::cout << "Test : " << BEreq::particle_code("d_3")
335  // << " " << BEreq::particle_code("u_3") << std::endl;
336  #endif
337 
338 
339  } // function RD_spectrum_from_ProcessCatalog
340 
341 
345  {
346  using namespace Pipes::RD_spectrum_ordered_func;
347 
348  result = *Dep::RD_spectrum;
349 
350  // NB: coannihilatingParticles does not necessarily have to be ordered,
351  // but it is assumed that coannihilatingParticles[0] is the DM particle
353  if (result.coannihilatingParticles.size() > 1)
354  for (std::size_t i=0; i<result.coannihilatingParticles.size()-1; i++)
355  {
356  for (std::size_t j=i+1; j<result.coannihilatingParticles.size(); j++)
357  {
358  if (result.coannihilatingParticles[j].mass<result.coannihilatingParticles[i].mass)
359  {
360  tmp_co=result.coannihilatingParticles[i];
362  result.coannihilatingParticles[j]=tmp_co;
363  }
364  }
365  }
366 
367 
368  // add coannihilation thresholds
369  if (result.coannihilatingParticles.size() > 1)
370  {
371  for (int i=0; i<(int)result.coannihilatingParticles.size(); i++)
372  {
373  for (int j=std::max(1,i); j<(int)result.coannihilatingParticles.size(); j++)
374  {
375  result.threshold_energy.push_back(
376  result.coannihilatingParticles[i].mass+result.coannihilatingParticles[j].mass);
377  }
378  }
379  }
380  //and order all thresholds
381  double tmp;
382  for (std::size_t i=0; i<result.threshold_energy.size()-1; i++)
383  {
384  for (std::size_t j=i+1; j<result.threshold_energy.size(); j++)
385  {
386  if (result.threshold_energy[j]<result.threshold_energy[i])
387  {
388  tmp=result.threshold_energy[i];
389  result.threshold_energy[i]=result.threshold_energy[j];
390  result.threshold_energy[j]=tmp;
391  }
392  }
393  }
394 
395  if (!result.resonances.empty()){
396  TH_Resonance tmp2;
397  for (std::size_t i=0; i<result.resonances.size()-1; i++)
398  {
399  for (std::size_t j=i+1; j<result.resonances.size(); j++)
400  {
401  if (result.resonances[j].energy < result.resonances[i].energy)
402  {
403  tmp2=result.resonances[i];
404  result.resonances[i]=result.resonances[j];
405  result.resonances[j]=tmp2;
406  }
407  }
408  }
409  }
410  } // function RD_spectrum_ordered_func
411 
412 
416  void RD_annrate_DS5prep_func(int &result)
417  {
418  using namespace Pipes::RD_annrate_DS5prep_func;
419 
420  // Read out coannihilating particles from RDspectrum.
421  RD_spectrum_type specres = *Dep::RD_spectrum;
422 
423  if ( specres.particle_index_type != "DarkSUSY" )
424  {
425  invalid_point().raise("RD_annrate_DS5prep_func is only optimized for use with "
426  "DarkSUSY5 and requires internal particle IDs. Try RD_annrate_DSprep_MSSM_func instead!");
427  }
428 
429  //write model-dependent info about coannihilating particles to DS common blocks
430  DS5_RDMGEV myrdmgev;
431  myrdmgev.nco = specres.coannihilatingParticles.size();
432  for (int i=1; i<=myrdmgev.nco; i++)
433  {
434  myrdmgev.mco(i)=fabs(specres.coannihilatingParticles[i-1].mass);
435  myrdmgev.mdof(i)=specres.coannihilatingParticles[i-1].degreesOfFreedom;
436  myrdmgev.kcoann(i)=specres.coannihilatingParticles[i-1].index;
437  }
438 
439  double tmp; int itmp;
440  for (int i=1; i<=myrdmgev.nco-1; i++)
441  {
442  for (int j=i+1; j<=myrdmgev.nco; j++)
443  {
444  if (myrdmgev.mco(j)<myrdmgev.mco(i))
445  {
446  tmp=myrdmgev.mco(i);
447  myrdmgev.mco(i)=myrdmgev.mco(j);
448  myrdmgev.mco(j)=tmp;
449  itmp=myrdmgev.mdof(i);
450  myrdmgev.mdof(i)=myrdmgev.mdof(j);
451  myrdmgev.mdof(j)=itmp;
452  itmp=myrdmgev.kcoann(i);
453  myrdmgev.kcoann(i)=myrdmgev.kcoann(j);
454  myrdmgev.kcoann(j)=itmp;
455  }
456  }
457  #ifdef DARKBIT_RD_DEBUG
458  std::cout << "DS5prep - co : "<< myrdmgev.kcoann(i) << " " <<
459  myrdmgev.mco(i) << " " << myrdmgev.mdof(i)
460  << std::endl;
461  #endif
462  }
463  #ifdef DARKBIT_RD_DEBUG
464  std::cout << "DS5prep - co : "<< myrdmgev.kcoann(myrdmgev.nco) << " " <<
465  myrdmgev.mco(myrdmgev.nco) << " " << myrdmgev.mdof(myrdmgev.nco)
466  << std::endl;
467  #endif
468 
469  *BEreq::rdmgev = myrdmgev;
470 
471  result=1; // everything OK
472 
473  } // function RD_annrate_DS5prep_func
474 
475 
479  void RD_annrate_DSprep_MSSM_func(int &result)
480  {
481  using namespace Pipes::RD_annrate_DSprep_MSSM_func;
482 
483  // Read out coannihilating particles from RDspectrum_ordered.
484  RD_spectrum_type specres = *Dep::RD_spectrum_ordered;
485 
486  if (specres.particle_index_type != "DarkSUSY" && specres.particle_index_type != "PDG")
487  {
488  invalid_point().raise("RD_annrate_DSprep_MSSM_func requires PDG or internal DS codes!");
489  }
490 
491  //write model-dependent info about coannihilating particles to DS common blocks
492  // Note: translation from PDG codes assumes for consistency context integer = 0 here
493  // (as done in RD_spectrum_MSSM)
494  int ContInt = 0;
495  DS_DSANCOANN mydsancoann;
496  mydsancoann.nco = specres.coannihilatingParticles.size();
497  int partID;
498  for (int i=1; i<=mydsancoann.nco; i++)
499  {
500  mydsancoann.mco(i)=fabs(specres.coannihilatingParticles[i-1].mass);
501  mydsancoann.mdof(i)=specres.coannihilatingParticles[i-1].degreesOfFreedom;
502  partID = specres.coannihilatingParticles[i-1].index;
503  mydsancoann.kco(i) = partID;
504  if (specres.particle_index_type == "PDG")
505  {
506  mydsancoann.kco(i) = BEreq::DSparticle_code(Models::ParticleDB().long_name(partID,ContInt));
507  };
508  #ifdef DARKBIT_RD_DEBUG
509  std::cout << "DS6prep_MSSM - co : "<< partID << " " << mydsancoann.kco(i) << " " <<
510  mydsancoann.mco(i) << " " << mydsancoann.mdof(i)
511  << std::endl;
512  #endif
513  }
514 
515  *BEreq::dsancoann = mydsancoann;
516 
517  result=1; // everything OK
518 
519  } // function RD_eff_annrate_DSprep_MSSM_func
520 
521 
526  void RD_eff_annrate_DS_MSSM(double(*&result)(double&))
527  {
528  using namespace Pipes::RD_eff_annrate_DS_MSSM;
529 
530  if (BEreq::dsanwx.origin() == "DarkSUSY_MSSM")
531  {
532  result=BEreq::dsanwx.pointer();
533  }
534  else DarkBit_error().raise(LOCAL_INFO, "Wrong DarkSUSY backend initialized?");
535  } // function RD_eff_annrate_DS_MSSM
536 
537  void RD_eff_annrate_DS5_MSSM(double(*&result)(double&))
538  {
539  using namespace Pipes::RD_eff_annrate_DS5_MSSM;
540 
541  if (BEreq::dsanwx.origin() == "DarkSUSY")
542  {
543  result=BEreq::dsanwx.pointer();
544  }
545  else DarkBit_error().raise(LOCAL_INFO, "Wrong DarkSUSY backend initialized?");
546  } // function RD_eff_annrate_DS5_MSSM
547 
548 
549 
552  // Carries pointer to Weff
553  DEF_FUNKTRAIT(RD_EFF_ANNRATE_FROM_PROCESSCATALOG_TRAIT)
554  void RD_eff_annrate_from_ProcessCatalog(double(*&result)(double&))
555  {
556  using namespace Pipes::RD_eff_annrate_from_ProcessCatalog;
557 
558  std::string DMid= *Dep::DarkMatter_ID;
559  TH_Process annProc = (*Dep::TH_ProcessCatalog).getProcess(DMid, DMid);
560  double mDM = (*Dep::TH_ProcessCatalog).getParticleProperty(DMid).mass;
561 
562  auto Weff = daFunk::zero("peff");
563  auto peff = daFunk::var("peff");
564  auto s = 4*(peff*peff + mDM*mDM);
565 
566  // Individual contributions to the invariant rate Weff. Note that no
567  // symmetry factor of 1/2 for non-identical initial state particles
568  // (non-self-conjugate DM) should appear here. This factor does explicitly
569  // enter, however, when calculating the relic density in RD_oh2_DS_general.
570  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
571  it != annProc.channelList.end(); ++it)
572  {
573  Weff = Weff +
574  it->genRate->set("v", 2*peff/sqrt(mDM*mDM+peff*peff))*s/gev2tocm3s1;
575  }
576  // Add genRateMisc to Weff
577  Weff = Weff + annProc.genRateMisc->set("v", 2*peff/sqrt(mDM*mDM+peff*peff))*s/gev2tocm3s1;
578  if ( Weff->getNArgs() != 1 )
579  DarkBit_error().raise(LOCAL_INFO,
580  "RD_eff_annrate_from_ProcessCatalog: Wrong number of arguments.\n"
581  "The probable cause are three-body final states, which are not supported for this function."
582  );
583  result = Weff->plain<RD_EFF_ANNRATE_FROM_PROCESSCATALOG_TRAIT>("peff");
584  } // function RD_eff_annrate_from_ProcessCatalog
585 
586 
594  void RD_oh2_DS_general(double &result)
595  {
596  using namespace Pipes::RD_oh2_DS_general;
597 
598  // Retrieve ordered list of resonances and thresholds from
599  // RD_thresholds_resonances.
600  RD_spectrum_type myRDspec = *Dep::RD_spectrum_ordered;
601  if (myRDspec.coannihilatingParticles.empty())
602  {
603  DarkBit_error().raise(LOCAL_INFO, "RD_oh2_DS_general: No DM particle!");
604  }
605  double mwimp=myRDspec.coannihilatingParticles[0].mass;
606 
607  // What follows below implements dsrdomega from DarkSUSY 6+
608  //We start by setting some general common block settings
609  BEreq::dsrdcom();
610  DS_RDPARS *myrdpars = BEreq::rdpars.pointer();
611 
614  int fast = runOptions->getValueOrDef<int>(1, "fast");
615 
618  BEreq::rdtime->rdt_max = runOptions->getValueOrDef<double>(600, "timeout");
619 
620  switch (fast)
621  {
622  case 0:
623  myrdpars->cosmin=0.996195;myrdpars->waccd=0.005;myrdpars->dpminr=1.0e-4;
624  myrdpars->dpthr=5.0e-4;myrdpars->wdiffr=0.05;myrdpars->wdifft=0.02;
625  break;
626  case 1:
627  myrdpars->cosmin=0.996195;myrdpars->waccd=0.05;myrdpars->dpminr=5.0e-4;
628  myrdpars->dpthr=2.5e-3;myrdpars->wdiffr=0.5;myrdpars->wdifft=0.1;
629  break;
630  default:
631  DarkBit_error().raise(LOCAL_INFO, "Invalid fast flag (should be 0 or 1). Fast > 1 not yet "
632  "supported in DarkBit::RD_oh2_DS_general. Please add relevant settings to this routine.");
633  }
634 
635  // now transfer information from myRDspec to DS common blocks
636  int tnco=myRDspec.coannihilatingParticles.size();
637  int tnrs=myRDspec.resonances.size();
638  int tnthr=myRDspec.threshold_energy.size();
639  double tmco[1000], tdof[1000], trm[1000], trw[1000], ttm[1000];
640  for (std::size_t i=0; i<((unsigned int)tnco); i++)
641  {
642  tmco[i] = myRDspec.coannihilatingParticles[i].mass;
643  tdof[i] = myRDspec.coannihilatingParticles[i].degreesOfFreedom;
644  #ifdef DARKBIT_RD_DEBUG
645  std::cout << "RD_oh2_DS_general - co : "<< tmco[i] << " " << tdof[i] << std::endl;
646  #endif
647  }
648  for (std::size_t i=0; i<((unsigned int)tnrs); i++)
649  {
650  trm[i] = myRDspec.resonances[i].energy;
651  trw[i] = myRDspec.resonances[i].width;
652  #ifdef DARKBIT_RD_DEBUG
653  std::cout << "RD_oh2_DS_general - res : "<< trm[i] << " " << trw[i] << std::endl;
654  #endif
655  }
656  //DS does not count 2* WIMP rest mass as thr, hence we start at i=1
657  tnthr -=tnthr;
658  for (std::size_t i=1; i<((unsigned int)tnthr+1); i++)
659  {
660  ttm[i] = myRDspec.threshold_energy[i];
661  #ifdef DARKBIT_RD_DEBUG
662  std::cout << "RD_oh2_DS_general - thr : "<< ttm[i] << std::endl;
663  #endif
664  }
665  #ifdef DARKBIT_RD_DEBUG
666  std::cout << "RD_oh2_DS_general - tnco,tnrs,tnthr : "<< tnco << " " << tnrs << " "
667  << tnthr << std::endl;
668  #endif
669  BEreq::dsrdstart(tnco,tmco,tdof,tnrs,trm,trw,tnthr,ttm);
670 
671  // always check that invariant rate is OK at least at one point
672  double peff = mwimp/100;
673  double weff = (*Dep::RD_eff_annrate)(peff);
674 
675  if (Utils::isnan(weff))
676  DarkBit_error().raise(LOCAL_INFO, "Weff is NaN in RD_oh2_DS_general. This means that the function\n"
677  "pointed to by RD_eff_annrate returned NaN for the invariant rate\n"
678  "entering the relic density calculation.");
679 
680  // Finally use DS Boltzmann solver with invariant rate
681  double oh2, xf;
682  int ierr=0; int iwar=0;
683  BEreq::dsrdens(byVal(*Dep::RD_eff_annrate),oh2,xf,fast,ierr,iwar);
684 
685  //Check for NAN result.
686  if ( Utils::isnan(oh2) ) DarkBit_error().raise(LOCAL_INFO, "DarkSUSY returned NaN for relic density!");
687 
688  // Check whether DarkSUSY threw an error
689  if (ierr == 1024)
690  {
691  invalid_point().raise("DarkSUSY invariant rate tabulation timed out.");
692  }
693  else if(ierr != 0)
694  {
695  DarkBit_error().raise(LOCAL_INFO, "DarkSUSY Boltzmann solver failed.");
696  }
697 
698  // If the DM particles are not their own antiparticles we need to add the relic
699  // density of anti-DM particles as well
700  result = (myRDspec.isSelfConj) ? oh2 : 2*oh2;
701 
702  logger() << LogTags::debug << "RD_oh2_DS_general: oh2 =" << result << EOM;
703 
704  #ifdef DARKBIT_DEBUG
705  std::cout << std::endl << "DM mass = " << mwimp<< std::endl;
706  std::cout << "Oh2 = " << result << std::endl << std::endl;
707  #endif
708 
709  } // function RD_oh2_DS_general
710 
711 
712 
720  void RD_oh2_DS5_general(double &result)
721  {
722  using namespace Pipes::RD_oh2_DS5_general;
723 
724  // Retrieve ordered list of resonances and thresholds from
725  // RD_thresholds_resonances.
726  RD_spectrum_type myRDspec = *Dep::RD_spectrum_ordered;
727  if (myRDspec.coannihilatingParticles.empty())
728  {
729  DarkBit_error().raise(LOCAL_INFO, "RD_oh2_DS5_general: No DM particle!");
730  }
731  double mwimp=myRDspec.coannihilatingParticles[0].mass;
732 
733  #ifdef DARKBIT_RD_DEBUG
734  bool tbtest=false;
735  #endif
736 
739  BEreq::rdtime->rdt_max = runOptions->getValueOrDef<double>(600, "timeout");
740 
741  // What follows below is the standard accurate calculation of oh2 in DS, in one of the
742  // following modes:
743  // fast = 0 - standard accurate calculation (accuracy better than 1%)
744  // 1 - faster calculation: sets parameters for when to add
745  // extra points less tough to avoid excessively
746  // adding extra points, expected accurarcy: 1% or better
747  // 2 - faster calculation: compared to fast=1, this option
748  // adds less points in Weff tabulation in general and
749  // is more elaborate in deciding when to include points
750  // close to thresholds and resonances
751  // expected accuracy: around 1%
752  // 3 - even more aggressive on trying minimize the number
753  // of tabulated points
754  // expected accuracy: 5-10%
755  // 9 - superfast. This method still makes sure to include
756  // resonances and threholds, but does not attempt to sample
757  // them very well. Should give an order of magnitude estimate
758  // expected accuracy: order of magnitude
759  // 10 - quick and dirty method, i.e. expand the annihilation
760  // cross section in x (not recommended)
761  // expected accuracy: can be orders of magnitude wrong
762  // for models with strong resonances or thresholds
763 
764  DS_RDPARS myrdpars;
767  int fast = runOptions->getValueOrDef<int>(1, "fast");
768  switch (fast)
769  {
770  case 0:
771  myrdpars.cosmin=0.996195;myrdpars.waccd=0.005;myrdpars.dpminr=1.0e-4;
772  myrdpars.dpthr=5.0e-4;myrdpars.wdiffr=0.05;myrdpars.wdifft=0.02;
773  break;
774  case 1:
775  myrdpars.cosmin=0.996195;myrdpars.waccd=0.05;myrdpars.dpminr=5.0e-4;
776  myrdpars.dpthr=2.5e-3;myrdpars.wdiffr=0.5;myrdpars.wdifft=0.1;
777  break;
778  default:
779  DarkBit_error().raise(LOCAL_INFO, "Invalid fast flag (should be 0 or 1). Fast > 1 not yet "
780  "supported in DarkBit::RD_oh2_DS5_general. Please add relevant settings to this routine.");
781  }
782 
783  myrdpars.hstep=0.01;myrdpars.hmin=1.0e-9;myrdpars.compeps=0.01;
784  myrdpars.xinit=2.0;myrdpars.xfinal=200.0;myrdpars.umax=10.0;
785  myrdpars.cfr=0.5;
786  *BEreq::rdpars = myrdpars;
787  DS_RDSWITCH myrdswitch;
788  myrdswitch.thavint=1;myrdswitch.rdprt=0;
789  *BEreq::rdswitch = myrdswitch;
790  DS_RDLUN myrdlun;
791  myrdlun.rdlulog=6;myrdlun.rdluerr=6;
792  *BEreq::rdlun = myrdlun;
793  DS_RDPADD myrdpadd;
794  myrdpadd.pdivr=2.0;myrdpadd.dpres=0.5;myrdpadd.nlow=20;
795  myrdpadd.nhigh=10;
796  myrdpadd.npres=4;myrdpadd.nthup=4;myrdpadd.cthtest=0;myrdpadd.spltest=1;
797  *BEreq::rdpadd = myrdpadd;
798 
799  DS_RDERRORS myrderrors;
800  myrderrors.rderr=0;myrderrors.rdwar=0;myrderrors.rdinit=1234;
801  *BEreq::rderrors = myrderrors;
802 
803  // write mass and dof of DM & coannihilating particle to DS common blocks
804  DS5_RDMGEV *myrdmgev = BEreq::rdmgev.pointer();
805 
806  myrdmgev->nco=myRDspec.coannihilatingParticles.size();
807  for (std::size_t i=1; i<=((unsigned int)myrdmgev->nco); i++)
808  {
809  myrdmgev->mco(i)=myRDspec.coannihilatingParticles[i-1].mass;
810  myrdmgev->mdof(i)=myRDspec.coannihilatingParticles[i-1].degreesOfFreedom;
811  myrdmgev->kcoann(i)=myRDspec.coannihilatingParticles[i-1].index;
812  #ifdef DARKBIT_RD_DEBUG
813  std::cout << "kcoann, mco, mdof: " << myrdmgev->kcoann(i) << " " << myrdmgev->mco(i)
814  << " " << myrdmgev->mdof(i) << std::endl;
815  #endif
816  }
817 
818  // write information about thresholds and resonances to DS common blocks
819  // [this is the model-independent part of dsrdstart]
820  myrdmgev->nres=0;
821  if (!myRDspec.resonances.empty())
822  {
823  myrdmgev->nres=myRDspec.resonances.size();
824  for (std::size_t i=1; i<=myRDspec.resonances.size(); i++)
825  {
826  myrdmgev->rgev(i)=myRDspec.resonances[i-1].energy;
827  myrdmgev->rwid(i)=myRDspec.resonances[i-1].width;
828  #ifdef DARKBIT_RD_DEBUG
829  std::cout << "rgev, rwid: " << myrdmgev->rgev(i) << " " << myrdmgev->rwid(i) << std::endl;
830  #endif
831  }
832  }
833  // convert to momenta and write to DS common blocks
834  DS_RDPTH myrdpth;
835  // NB: DS does not count 2* WIMP rest mass as thr
836  myrdpth.nth=myRDspec.threshold_energy.size()-1;
837  myrdpth.pth(0)=0; myrdpth.incth(0)=1;
838  for (std::size_t i=1; i<myRDspec.threshold_energy.size(); i++)
839  {
840  myrdpth.pth(i)=sqrt(pow(myRDspec.threshold_energy[i],2)/4-pow(mwimp,2));
841  myrdpth.incth(i)=1;
842  #ifdef DARKBIT_RD_DEBUG
843  std::cout << "pth, incth: " << myrdpth.pth(i) << " " << myrdpth.incth(i) << std::endl;
844  #endif
845  }
846  *BEreq::rdpth = myrdpth;
847 
848  // determine starting point for integration of Boltzmann eq and write
849  // to DS common blocks
850  DS_RDDOF *myrddof = BEreq::rddof.pointer();
851  double xstart=std::max(myrdpars.xinit,1.0001*mwimp/myrddof->tgev(1));
852  double tstart=mwimp/xstart;
853  int k; myrddof->khi=myrddof->nf; myrddof->klo=1;
854  while (myrddof->khi > myrddof->klo+1)
855  {
856  k=(myrddof->khi+myrddof->klo)/2;
857  if (myrddof->tgev(k) < tstart)
858  {
859  myrddof->khi=k;
860  }
861  else {
862  myrddof->klo=k;
863  }
864  }
865 
866  // follow wide res treatment for heavy Higgs adopted in DS
867  double widthheavyHiggs = BEreq::widths->width(BEreq::DS5particle_code("h0_2"));
868  if (widthheavyHiggs<0.1) BEreq::widths->width(BEreq::DS5particle_code("h0_2"))=0.1;
869 
870  // always check that invariant rate is OK at least at one point
871  double peff = mwimp/100;
872  double weff = (*Dep::RD_eff_annrate)(peff);
873  if (Utils::isnan(weff))
874  DarkBit_error().raise(LOCAL_INFO, "Weff is NaN in RD_oh2_DS5_general. This means that the function\n"
875  "pointed to by RD_eff_annrate returned NaN for the invariant rate\n"
876  "entering the relic density calculation.");
877 
878  #ifdef DARKBIT_RD_DEBUG
879  // Dump Weff info on screen
880  std::cout << "xstart = " << xstart << std::endl;
881  for ( peff = mwimp/1000; peff < mwimp; peff = peff*1.5 )
882  {
883  weff = (*Dep::RD_eff_annrate)(peff);
884  std::cout << "Weff(" << peff << ") = " << weff << std::endl;
885  // Check that the invariant rate is OK.
886  if (Utils::isnan(weff))
887  DarkBit_error().raise(LOCAL_INFO, "RD debug: Weff is NaN in RD_oh2_DS5_general.");
888  }
889  // Set up timing
890  std::chrono::time_point<std::chrono::system_clock> start, end;
891  start = std::chrono::system_clock::now();
892  logger() << "Tabulating RD_eff_annrate..." << EOM;
893  std::cout << "Starting dsrdtab..." << std::endl;
894  #endif
895 
896  // Tabulate the invariant rate
897  BEreq::dsrdtab(byVal(*Dep::RD_eff_annrate),xstart,fast);
898 
899  #ifdef DARKBIT_RD_DEBUG
900  logger() << LogTags::repeat_to_cout << "...done!" << EOM;
901 
902  // Get runtime
903  end = std::chrono::system_clock::now();
904  double runtime = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
905 
906  // Check if runtime too long
907  if ( runtime > 30. )
908  {
909  std::cout << "Duration [ms]: " << runtime << std::endl;
910  //SLHAstruct mySLHA = Dep::MSSM_spectrum->getSLHAea(2);
911  //std::ofstream ofs("RelicDensity_debug.slha");
912  //ofs << mySLHA;
913  //ofs.close();
914  tbtest=true;
915  }
916  #endif
917 
918  // Check whether DarkSUSY threw an error
919  if (BEreq::rderrors->rderr != 0)
920  {
921  if (BEreq::rderrors->rderr == 1024) invalid_point().raise("DarkSUSY invariant rate tabulation timed out.");
922  else DarkBit_error().raise(LOCAL_INFO, "DarkSUSY invariant rate tabulation failed.");
923  }
924 
925  // determine integration limit
926  BEreq::dsrdthlim();
927 
928  // now solve Boltzmann eqn using tabulated rate
929  double xend, yend, xf; int nfcn;
930  BEreq::dsrdeqn(byVal(BEreq::dsrdwintp.pointer()), xstart,xend,yend,xf,nfcn);
931  // using the untabulated rate gives the same result but is usually
932  // slower:
933  // BEreq::dsrdeqn(byVal(*Dep::RD_eff_annrate),xstart,xend,yend,xf,nfcn);
934 
935  // change heavy Higgs width in DS back to standard value
936  BEreq::widths->width(BEreq::DS5particle_code("h0_2")) = widthheavyHiggs;
937 
938  //Check for NAN result.
939  if ( Utils::isnan(yend) ) DarkBit_error().raise(LOCAL_INFO, "DarkSUSY returned NaN for relic density!");
940 
941  // Check whether DarkSUSY threw some other error
942  if (BEreq::rderrors->rderr != 0) DarkBit_error().raise(LOCAL_INFO, "DarkSUSY Boltzmann solver failed.");
943 
944  result = 0.70365e8*myrddof->fh(myrddof->nf)*mwimp*yend;
945 
946  // If the DM particles are not their own antiparticles we need to add the relic
947  // density of anti-DM particles as well
948  result = (myRDspec.isSelfConj) ? result : 2*result;
949 
950  logger() << LogTags::debug << "RD_oh2_DS5_general: oh2 =" << result << EOM;
951 
952  #ifdef DARKBIT_DEBUG
953  std::cout << std::endl << "DM mass = " << mwimp<< std::endl;
954  std::cout << "Oh2 = " << result << std::endl << std::endl;
955  #endif
956 
957  #ifdef DARKBIT_RD_DEBUG
958  if (tbtest) exit(1);
959  #endif
960 
961  } // function RD_oh2_DS5_general
962 
963 
964 
965 
967  //
968  // Simple relic density routines for cross-checks
969  // (MicrOmegas vs DarkSUSY)
970  //
972 
976  {
977  using namespace Pipes::RD_oh2_Xf_MicrOmegas;
978  // Input
979  int fast; // fast: 1, accurate: 0
980  double Beps; // Beps=1e-5 recommended, Beps=1 switches coannihilation off
981 
982  // Set options via ini-file (MicrOmegas-specific performance options)
983  fast = runOptions->getValueOrDef<int>(0, "fast");
984  Beps = runOptions->getValueOrDef<double>(1e-5, "Beps");
985 
986  logger() << LogTags::debug << "Using fast: " << fast << " Beps: " << Beps;
987 
988  // Output
989  double Xf;
990  double oh2 = BEreq::oh2(&Xf,byVal(fast), byVal(Beps));
991 
992  result.first = oh2;
993  result.second = Xf;
994 
995  logger() << LogTags::debug << "X_f = " << Xf << " Omega h^2 = " << oh2 << EOM;
996  }
997 
1000  void RD_oh2_DarkSUSY_DS5(double &result)
1001  {
1002  using namespace Pipes::RD_oh2_DarkSUSY_DS5;
1003  // Input
1004  int omtype; // 0: no coann; 1: all coann
1005  int fast; // 0: standard; 1: fast; 2: dirty
1006 
1007  // Set options via ini-file
1009  omtype = runOptions->getValueOrDef<int>(1, "omtype");
1011  fast = runOptions->getValueOrDef<int>(0, "fast");
1014  BEreq::rdtime->rdt_max = runOptions->getValueOrDef<double>(600, "timeout");
1015 
1016  // Output
1017  double xf; // freeze-out temperature
1018  int ierr; // error flag
1019  int iwar; // warming flag
1020  int nfc; // number of fnct calls to effective annihilation cross section
1021  logger() << LogTags::debug << "Starting DarkSUSY relic density calculation..." << EOM;
1022  double oh2 = BEreq::dsrdomega(omtype,fast,xf,ierr,iwar,nfc);
1023 
1024  // Check whether DarkSUSY threw an error
1025  if (BEreq::rderrors->rderr != 0)
1026  {
1027  if (BEreq::rderrors->rderr == 1024) invalid_point().raise("DarkSUSY invariant rate tabulation timed out.");
1028  else DarkBit_error().raise(LOCAL_INFO, "DarkSUSY relic density calculation failed.");
1029  }
1030 
1031  result = oh2;
1032  logger() << LogTags::debug << "RD_oh2_DarkSUSY_DS5: oh2 is " << oh2 << EOM;
1033  }
1034 
1035 
1036 
1037  void RD_oh2_MicrOmegas(double &result)
1038  {
1039  using namespace Pipes::RD_oh2_MicrOmegas;
1040 
1041  ddpair oh2_Xf = *Dep::RD_oh2_Xf;
1042  result = oh2_Xf.first;
1043  }
1044 
1045  void Xf_MicrOmegas(double &result)
1046  {
1047  using namespace Pipes::Xf_MicrOmegas;
1048 
1049  ddpair oh2_Xf = *Dep::RD_oh2_Xf;
1050  result = oh2_Xf.second;
1051  }
1052 
1053 
1055  {
1057 
1058  double Beps; // Beps=1e-5 recommended, Beps=1 switches coannihilation off
1059  Beps = runOptions->getValueOrDef<double>(1e-5, "Beps");
1060 
1061  double Xf = *Dep::Xf;
1062 
1063  double cut = runOptions->getValueOrDef<double>(1e-5, "cut");
1064 
1065  result = BEreq::momegas_print_channels(byVal(Xf),byVal(cut),byVal(Beps),byVal(1),byVal(stdout));
1066  }
1067 
1068 
1069  void get_semi_ann_MicrOmegas(double &result)
1070  {
1071  using namespace Pipes::get_semi_ann_MicrOmegas;
1072 
1073  double Beps; // Beps=1e-5 recommended, Beps=1 switches coannihilation off
1074  Beps = runOptions->getValueOrDef<double>(1e-5, "Beps");
1075 
1076  double Xf = *Dep::Xf;
1077 
1078  char*n1 = (char *)"~SS";
1079  char*n2 = (char *)"~SS";
1080  char*n3 = (char *)"h";
1081  char*n4 = (char *)"~ss";
1082 
1083  result = BEreq::get_oneChannel(byVal(Xf),byVal(Beps),byVal(n1),byVal(n2),byVal(n3),byVal(n4));
1084 
1085  }
1086 
1087 
1088 
1090  //
1091  // Infer fraction of Dark matter that is made up by scanned DM particles
1092  //
1094 
1095  void RD_fraction_one(double &result)
1096  {
1097  result = 1.0;
1098  logger() << LogTags::debug << "Fraction of dark matter that the scanned model accounts for: " << result << EOM;
1099  }
1100 
1101  void RD_fraction_leq_one(double &result)
1102  {
1103  using namespace Pipes::RD_fraction_leq_one;
1105  double oh2_obs = runOptions->getValueOrDef<double>(0.1188, "oh2_obs");
1106  double oh2_theory = *Dep::RD_oh2;
1107  result = std::min(1., oh2_theory/oh2_obs);
1108  logger() << LogTags::debug << "Fraction of dark matter that the scanned model accounts for: " << result << EOM;
1109  }
1110 
1111  void RD_fraction_rescaled(double &result)
1112  {
1113  using namespace Pipes::RD_fraction_rescaled;
1115  double oh2_obs = runOptions->getValueOrDef<double>(0.1188, "oh2_obs");
1116  double oh2_theory = *Dep::RD_oh2;
1117  result = oh2_theory/oh2_obs;
1118  logger() << LogTags::debug << "Fraction of dark matter that the scanned model accounts for: " << result << EOM;
1119  }
1120 
1121  }
1122 }
std::vector< TH_Channel > channelList
List of channels.
error & DarkBit_error()
void RD_oh2_DarkSUSY_DS5(double &result)
Relic density directly from a call of initialized DarkSUSY 5.
std::vector< TH_Resonance > resonances
Entry & at(std::pair< int, int >)
Get entry in decay table for a give particle, throwing an error if particle is absent.
void get_semi_ann_MicrOmegas(double &result)
void RD_annrate_DS5prep_func(int &result)
Some helper function to prepare evaluation of Weff from DarkSUSY 5.
void RD_fraction_one(double &result)
std::vector< RD_coannihilating_particle > coannihilatingParticles
#define LOCAL_INFO
Definition: local_info.hpp:34
str long_name(str, int) const
Retrieve the long name, from the short name and index.
Definition: partmap.cpp:116
void RD_spectrum_SUSY_DS5(RD_spectrum_type &result)
Collects spectrum information about coannihilating particles, resonances and threshold energies – di...
void RD_oh2_Xf_MicrOmegas(ddpair &result)
Relic density directly from a call of initialized MicrOmegas.
void print_channel_contributions_MicrOmegas(double &result)
Funk var(std::string arg)
Definition: daFunk.hpp:921
void RD_oh2_DS5_general(double &result)
General routine for calculation of relic density, using DarkSUSY 5 Boltzmann solver.
void Xf_MicrOmegas(double &result)
std::pair< double, double > ddpair
Shorthand for a pair of doubles.
Definition: util_types.hpp:66
void RD_spectrum_MSSM(RD_spectrum_type &result)
Collects spectrum information about coannihilating particles, resonances and threshold energies...
bool isSelfConj
Does the process contain self-conjugate DM? (accounting for correct factors of 1/2 in annihilation sp...
void RD_spectrum_from_ProcessCatalog(RD_spectrum_type &result)
Collects information about resonances and threshold energies directly from the ProcessCatalog [NB: th...
General small utility functions.
A single resonance of a given width at a given energy (both in GeV)
void RD_eff_annrate_DS5_MSSM(double(*&result)(double &))
double mass
Particle mass (GeV)
std::vector< double > threshold_energy
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
virtual double get(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
TH_resonances_thresholds resonances_thresholds
List of resonances and thresholds.
double width_in_GeV
Total particle width (in GeV)
partmap & ParticleDB()
Database accessor function.
Definition: partmap.cpp:32
DEF_FUNKTRAIT(RD_EFF_ANNRATE_FROM_PROCESSCATALOG_TRAIT) void RD_eff_annrate_from_ProcessCatalog(double(*&result)(double &))
Infer Weff from process catalog.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
Header file that includes all GAMBIT headers required for a module source file.
A container for the mass and spin of a particle.
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
void RD_fraction_leq_one(double &result)
unsigned int spin2
Twice the spin of the particle.
SubSpectrum & get_LE()
Standard getters Return references to internal data members.
Definition: spectrum.cpp:224
daFunk::Funk genRateMisc
Additional decay rate or sigmav (in addition to above channels)
DS5_MSPCTM DS_INTDOF int
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
Rollcall header for module DarkBit.
GAMBIT native decay table class.
Definition: decay_table.hpp:35
double pow(const double &a)
Outputs a^i.
void RD_oh2_MicrOmegas(double &result)
invalid_point_exception & invalid_point()
Invalid point exceptions.
Funk zero(Args... argss)
Definition: daFunk.hpp:604
A container for a single process.
T byVal(T t)
Redirection function to turn an lvalue into an rvalue, so that it is correctly passed by value when d...
const double gev2tocm3s1
void RD_fraction_rescaled(double &result)
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...
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
std::vector< TH_Resonance > resonances
SubSpectrum & get_HE()
Definition: spectrum.cpp:225
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110
Utility functions for DarkBit.