gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
DarkBit.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
43 
47 
48 namespace Gambit
49 {
50  namespace DarkBit
51  {
52 
54  //
55  // Simple DM property extractors
56  //
58 
60  void mwimp_generic(double &result)
61  {
62  using namespace Pipes::mwimp_generic;
63  result = Dep::TH_ProcessCatalog->getParticleProperty(*Dep::DarkMatter_ID).mass;
64  if (result < 0.0) DarkBit_error().raise(LOCAL_INFO, "Negative WIMP mass detected.");
65  }
66 
70  void sigmav_late_universe(double &result)
71  {
72  using namespace Pipes::sigmav_late_universe;
73  std::string DMid = *Dep::DarkMatter_ID;
74  TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMid);
75  result = 0.0;
76  // Add all the regular channels
77  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
78  it != annProc.channelList.end(); ++it)
79  {
80  if ( it->nFinalStates == 2 )
81  {
82  // (sv)(v=0) for two-body final state
83  result += it->genRate->bind("v")->eval(0.);
84  }
85  }
86  // Add invisible contributions
87  result += annProc.genRateMisc->bind("v")->eval(0.);
88  }
89 
90  void sigmav_late_universe_MicrOmegas(double &result)
91  {
93 
94  // set requested spectra to NULL since we don't need them
95  double * SpNe=NULL,*SpNm=NULL,*SpNl=NULL;
96  double * SpA=NULL,*SpE=NULL,*SpP=NULL;
97  int err;
98 
99  // TODO: Add option to choose options for calculation from MicrOmegas (right now all effects are turned on)
100  result = BEreq::calcSpectrum(byVal(3),byVal(SpA),byVal(SpE),byVal(SpP),byVal(SpNe),byVal(SpNm),byVal(SpNl) ,byVal(&err));
101 
102  if (err != 0 )
103  {
104  DarkBit_error().raise(LOCAL_INFO, "MicrOmegas spectrum calculation returned error code = " + std::to_string(err));
105  }
106 
107  }
108 
109 
111  //
112  // Extraction of local and global dark matter halo properties
113  //
115 
116 
118  double profile_gNFW(double rhos, double rs, double alpha, double beta, double gamma, double r)
119  { return pow(2, (beta-gamma)/alpha)*rhos/pow(r/rs, gamma)/pow(1+pow(r/rs, alpha), (beta-gamma)/alpha); }
120 
122  double profile_Einasto(double rhos, double rs, double alpha, double r)
123  { return rhos*exp((-2.0/alpha)*(pow(r/rs, alpha)-1)); }
124 
127  {
128  using namespace Pipes::GalacticHalo_gNFW;
129  double rhos = *Param["rhos"];
130  double rs = *Param["rs"];
131  double r_sun = *Param["r_sun"];
132  double alpha = *Param["alpha"];
133  double beta = *Param["beta"];
134  double gamma = *Param["gamma"];
135  result.DensityProfile = daFunk::func(profile_gNFW, rhos, rs, alpha, beta, gamma, daFunk::var("r"));
136  result.r_sun = r_sun;
137  }
138 
141  {
142  using namespace Pipes::GalacticHalo_Einasto;
143  double rhos = *Param["rhos"];
144  double rs = *Param["rs"];
145  double r_sun = *Param["r_sun"];
146  double alpha = *Param["alpha"];
147  result.DensityProfile = daFunk::func(profile_Einasto, rhos, rs, alpha, daFunk::var("r"));
148  result.r_sun = r_sun;
149  }
150 
153  {
154  using namespace Pipes::ExtractLocalMaxwellianHalo;
155  double rho0 = *Param["rho0"];
156  double v0 = *Param["v0"];
157  double vesc = *Param["vesc"];
158  double vrot = *Param["vrot"];
159  result.rho0 = rho0;
160  result.v0 = v0;
161  result.vesc = vesc;
162  result.vrot = vrot;
163  }
164 
165 
166 
168  //
169  // DarkBit Unit Test
170  //
172 
177  void UnitTest_DarkBit(int &result)
178  {
179  using namespace Pipes::UnitTest_DarkBit;
181  /* This function depends on all relevant DM observables (indirect and
182  * direct) and dumps them into convenient files in YAML format, which
183  * afterwards can be checked against the expectations.
184  */
185 
186  double M_DM =
187  Dep::TH_ProcessCatalog->getParticleProperty(*Dep::DarkMatter_ID).mass;
188  double Gps = (*Dep::DD_couplings).gps;
189  double Gpa = (*Dep::DD_couplings).gpa;
190  double Gns = (*Dep::DD_couplings).gns;
191  double Gna = (*Dep::DD_couplings).gna;
192  double oh2 = *Dep::RD_oh2;
193 
194  std::string DMid = *Dep::DarkMatter_ID;
195  TH_Process annProc = (*Dep::TH_ProcessCatalog).getProcess(DMid, DMid);
196  daFunk::Funk spectrum = (*Dep::GA_AnnYield)->set("v", 0.);
197 
198  std::ostringstream filename;
199  /*
200  static unsigned int counter = 0;
201  filename << runOptions->getValueOrDef<std::string>("UnitTest_DarkBit",
202  "fileroot");
203  filename << "_" << counter << ".yml";
204  counter++;
205  */
207  filename << runOptions->getValueOrDef<std::string>("UnitTest.yaml", "filename");
208 
209  std::ofstream os;
210  os.open(filename.str());
211  if(os)
212  {
213  // Standard output.
214  os << "# Direct detection couplings\n";
215  os << "DDcouplings:\n";
216  os << " gps: " << Gps << "\n";
217  os << " gpa: " << Gpa << "\n";
218  os << " gns: " << Gns << "\n";
219  os << " gna: " << Gna << "\n";
220  os << "\n";
221  os << "# Particle masses [GeV] \n";
222  os << "ParticleMasses:\n";
223  os << " Mchi: " << M_DM << "\n";
224  os << "\n";
225  os << "# Relic density Omega h^2\n";
226  os << "RelicDensity:\n";
227  os << " oh2: " << oh2 << "\n";
228  os << "\n";
229 
230  // Output gamma-ray spectrum (grid be set in YAML file).
231  double x_min =
233  runOptions->getValueOrDef<double>(0.1, "GA_AnnYield", "Emin");
234  double x_max =
236  runOptions->getValueOrDef<double>(10000, "GA_AnnYield", "Emax");
238  int n = runOptions->getValueOrDef<double>(26, "GA_AnnYield", "nbins");
239  // from 0.1 to 500 GeV
240  std::vector<double> x = daFunk::logspace(log10(x_min), log10(x_max), n);
241  std::vector<double> y = spectrum->bind("E")->vect(x);
242  os << "# Annihilation spectrum dNdE [1/GeV]\n";
243  os << "GammaRaySpectrum:\n";
244  os << " E: [";
245  for (std::vector<double>::iterator it = x.begin(); it != x.end(); it++)
246  os << *it << ", ";
247  os << "]\n";
248  os << " dNdE: [";
249  for (std::vector<double>::iterator it = y.begin(); it != y.end(); it++)
250  os << *it << ", ";
251  os << "]\n";
252  os << std::endl;
253 
254  os << "# Annihilation rates\n";
255  os << "AnnihilationRates:\n";
256  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
257  it != annProc.channelList.end(); ++it)
258  {
259  os << " ";
260  for (std::vector<std::string>::iterator
261  jt = it->finalStateIDs.begin(); jt!=it->finalStateIDs.end(); jt++)
262  {
263  os << *jt << "";
264  }
265  if (it->finalStateIDs.size() == 2)
266  os << ": " << it->genRate->bind("v")->eval(0);
267  if (it->finalStateIDs.size() == 3)
268  {
269  double m1 = (*Dep::TH_ProcessCatalog).getParticleProperty(
270  it->finalStateIDs[1]).mass;
271  double m2 = (*Dep::TH_ProcessCatalog).getParticleProperty(
272  it->finalStateIDs[2]).mass;
273  double mass = M_DM;
274  daFunk::Funk E1_low = daFunk::func(gamma3bdy_limits<0>, daFunk::var("E"), mass, m1, m2);
275  daFunk::Funk E1_high = daFunk::func(gamma3bdy_limits<1>, daFunk::var("E"), mass, m1, m2);
276  daFunk::Funk dsigmavde = it->genRate->gsl_integration("E1", E1_low, E1_high)->set_epsrel(1e-3)->set("v", 0);
277  auto xgrid = daFunk::logspace(-2, 3, 1000);
278  auto ygrid = daFunk::logspace(-2, 3, 1000);
279  for ( size_t i = 0; i<xgrid.size(); i++ )
280  {
281  ygrid[i] = dsigmavde->bind("E")->eval(xgrid[i]);
282  }
283  auto interp = daFunk::interp("E", xgrid, ygrid);
284  double svTOT = interp->gsl_integration("E", 10., 20.)->set_epsabs(1e-3)->bind()->eval();
285  os << ": " << svTOT;
286  }
287  os << "\n";
288  }
289  os << std::endl;
290  }
291  else
292  {
293  logger() << "Warning: outputfile not open for writing." << EOM;
294  }
295  os.close();
296  result = 0;
297  }
298  }
299 }
std::vector< TH_Channel > channelList
List of channels.
error & DarkBit_error()
template double gamma3bdy_limits< 0 >(double, double, double, double)
double profile_Einasto(double rhos, double rs, double alpha, double r)
Einasto dark matter halo profile function.
Definition: DarkBit.cpp:122
START_MODEL beta
Definition: Axions.hpp:35
void sigmav_late_universe(double &result)
Retrieve the total thermally-averaged annihilation cross-section for indirect detection (cm^3 / s)...
Definition: DarkBit.cpp:70
START_MODEL r_sun
START_MODEL rs
#define LOCAL_INFO
Definition: local_info.hpp:34
START_MODEL v0
Funk var(std::string arg)
Definition: daFunk.hpp:921
void GalacticHalo_gNFW(GalacticHaloProperties &result)
Module function to generate GalacticHaloProperties for gNFW profile.
Definition: DarkBit.cpp:126
START_MODEL vesc
START_MODEL alpha
MATH_OPERATION(Dif,-) MATH_OPERATION(pow) MATH_OPERATION(fmin) MATH_OPERATION(fmax) class FunkInterp shared_ptr< FunkInterp > interp(T f, std::vector< double > x, std::vector< double > y)
Definition: daFunk.hpp:1349
std::vector< double > logspace(double x0, double x1, unsigned int n)
Definition: daFunk.hpp:186
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
Funk func(double(*f)(funcargs...), Args... args)
Definition: daFunk.hpp:768
template double gamma3bdy_limits< 1 >(double, double, double, double)
Header file that includes all GAMBIT headers required for a module source file.
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
START_MODEL dNur_CMB r
double profile_gNFW(double rhos, double rs, double alpha, double beta, double gamma, double r)
Generalized NFW dark matter halo profile function.
Definition: DarkBit.cpp:118
void sigmav_late_universe_MicrOmegas(double &result)
Definition: DarkBit.cpp:90
START_MODEL vrot
daFunk::Funk genRateMisc
Additional decay rate or sigmav (in addition to above channels)
double gamma3bdy_limits(double Eg, double M_DM, double m1, double m2)
Calculate kinematical limits for three-body final states.
void UnitTest_DarkBit(int &result)
Central unit test routine.
Definition: DarkBit.cpp:177
void ExtractLocalMaxwellianHalo(LocalMaxwellianHalo &result)
Module function providing local density and velocity dispersion parameters.
Definition: DarkBit.cpp:152
Rollcall header for module DarkBit.
double pow(const double &a)
Outputs a^i.
void mwimp_generic(double &result)
Retrieve the DM mass in GeV for generic models (GeV)
Definition: DarkBit.cpp:60
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...
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
void GalacticHalo_Einasto(GalacticHaloProperties &result)
Module function to generate GalacticHaloProperties for Einasto profile.
Definition: DarkBit.cpp:140
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Utility functions for DarkBit.