gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
SimpleLikelihoods.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
20 
24 
25 namespace Gambit {
26  namespace DarkBit {
27 
28 // //////////////////////////////////////////////////////////////////////////
29 // //
30 // // Simple likelihood functions
31 // //
32 // //////////////////////////////////////////////////////////////////////////
33 //
34 // /*! \brief Fermi LAT dwarf likelihoods, based on arXiv:1108.2914.
35 // */
36 // void lnL_FermiLATdwarfsSimple(double &result)
37 // {
38 // using namespace Pipes::lnL_FermiLATdwarfsSimple;
39 // // Koushiappas' limits [arXiv:1108.2914]
40 // //
41 // // This is the tabulated Phi-Likelihood function from Koushiappas et al.
42 // // Above L = 36, we use linear extrapolation up to L = 360000
43 // //
44 // // phi (defined as phi = sigmav/mDM**2*Ntot/8/pi * 1e26)
45 // double xgridArray [101] = { 0. , 6.74308086122e-05 , 0.000123192463137 ,
46 // 0.000171713798503 , 0.000215245918518 , 0.000255093268618 , 0.00029207805123 ,
47 // 0.000326751732695 , 0.000359503469472 , 0.000390620122006 , 0.000420321264006,
48 // 0.00044878042576 , 0.000476138421008 , 0.000502511975672 , 0.000527999496499,
49 // 0.000552685056887 , 0.000576641243501 , 0.000599931255273 ,
50 // 0.000622610497068 , 0.000644727821172 , 0.000666326515638 , 0.000687445105269,
51 // 0.000708118010141 , 0.000728376093388 , 0.000748247120993 , 0.00076775615078,
52 // 0.000786925863514 , 0.000805776846231 , 0.000824327835809 , 0.00084259592922,
53 // 0.000860596765645 , 0.000878344684789 , 0.000895852864914 ,
54 // 0.000913133443547 , 0.000930197623331 , 0.0009470557651 , 0.000963717469925 ,
55 // 0.00098019165163 , 0.000996486601006 , 0.00101261004288 , 0.00102856918685 ,
56 // 0.00104437077256 , 0.00106002111016 , 0.00107552611658 , 0.00109089134805 ,
57 // 0.00110612202935 , 0.00112122308019 , 0.00113619913897 , 0.00115105458439 ,
58 // 0.00116579355487 , 0.00118041996631 , 0.00119493752815 , 0.00120934975806 ,
59 // 0.00122365999528 , 0.00123787141289 , 0.00125198702892 , 0.00126600971667 ,
60 // 0.00127994221404 , 0.00129378713223 , 0.00130754696367 , 0.00132122408935 ,
61 // 0.00133482078559 , 0.00134833923028 , 0.00136178150869 , 0.0013751496188 ,
62 // 0.00138844547626 , 0.00140167091906 , 0.00141482771173 , 0.00142791754942 ,
63 // 0.00144094206154 , 0.0014539028153 , 0.00146680131887 , 0.00147963902447 ,
64 // 0.00149241733116 , 0.00150513758749 , 0.00151780109399 , 0.00153040910553 ,
65 // 0.00154296283341 , 0.00155546344754 , 0.00156791207827 , 0.00158030981824 ,
66 // 0.00159265772411 , 0.00160495681814 , 0.00161720808976 , 0.00162941249692 ,
67 // 0.00164157096757 , 0.00165368440081 , 0.00166575366823 , 0.00167777961494 ,
68 // 0.00168976306076 , 0.00170170480119 , 0.00171360560841 , 0.00172546623219 ,
69 // 0.00173728740083 , 0.00174906982191 , 0.00176081418314 , 0.00177252115315 ,
70 // 0.00178419138212 , 0.00179582550256 , 0.00180742412988 , 18.0 };
71 // //
72 // // Normalization w.r.t. p-value of phi=0
73 // //
74 // // chi^2
75 // double ygridArray [101] = { 0.0,
76 // 0.0513551, 0.177438, 0.35228, 0.561353, 0.795726, 1.04953, 1.3187, 1.60032,
77 // 1.89222, 2.19274, 2.50059, 2.81476, 3.13441, 3.45887, 3.78757, 4.12006,
78 // 4.45594, 4.79486, 5.13653, 5.48072, 5.82719, 6.17576, 6.52625, 6.87853,
79 // 7.23244, 7.58789, 7.94475, 8.30294, 8.66236, 9.02294, 9.38462, 9.74731,
80 // 10.111, 10.4755, 10.841, 11.2072, 11.5742, 11.9419, 12.3104, 12.6795, 13.0492,
81 // 13.4195, 13.7904, 14.1619, 14.5339, 14.9063, 15.2793, 15.6527, 16.0266,
82 // 16.4008, 16.7755, 17.1506, 17.5261, 17.9019, 18.2781, 18.6546, 19.0315,
83 // 19.4087, 19.7861, 20.1639, 20.542, 20.9203, 21.2989, 21.6778, 22.0569,
84 // 22.4362, 22.8158, 23.1957, 23.5757, 23.956, 24.3365, 24.7171, 25.098, 25.4791,
85 // 25.8604, 26.2418, 26.6235, 27.0053, 27.3872, 27.7694, 28.1517, 28.5342,
86 // 28.9168, 29.2996, 29.6825, 30.0655, 30.4487, 30.8321, 31.2155, 31.5992,
87 // 31.9829, 32.3667, 32.7507, 33.1348, 33.519, 33.9034, 34.2878, 34.6724,
88 // 35.0571, 350000.0 };
89 // // Convert arrays to vectors.
90 // std::vector<double> xgrid(xgridArray,
91 // xgridArray + sizeof xgridArray / sizeof xgridArray[0]);
92 // std::vector<double> ygrid(ygridArray,
93 // ygridArray + sizeof ygridArray / sizeof ygridArray[0]);
94 // // Construct interpolated function, using GAMBIT base functions.
95 // auto dwarf_likelihood = daFunk::interp("phi", xgrid, ygrid);
96 //
97 // double fraction = *Dep::RD_fraction;
98 //
99 // // Integate spectrum
100 // // (the zero velocity limit of the differential annihilation
101 // // cross-section as function of individual final state photons)
102 // double AnnYieldint = (*Dep::GA_AnnYield)->
103 // set("v", 0.)->gsl_integration("E", 1, 100)->set_epsabs(0)->set_epsrel(1e-3)->bind()->eval();
104 // logger() << "AnnYieldInt (1-100 GeV): " << AnnYieldint << EOM;
105 //
106 // // Calculate phi-value
107 // double phi = AnnYieldint / 8. / M_PI * 1e26 * fraction * fraction;
108 //
109 // // And return final likelihood
110 // result = 0.5*dwarf_likelihood->bind("phi")->eval(phi);
111 // logger() << "dwarf_likelihood: " << result << EOM;
112 // logger() << "phi: " << phi << EOM;
113 // }
114 
115 
116  // module function which sets the Milky Way halo profile for the gamlike backend
117  void set_gamLike_GC_halo(bool &result)
118  {
119  using namespace Pipes::set_gamLike_GC_halo;
120 
121  daFunk::Funk profile = (Dep::GalacticHalo)->DensityProfile;
122  auto r = daFunk::logspace(-3, 2, 100);
123  auto rho = daFunk::logspace(-3, 2, 100);
124  double dist = (Dep::GalacticHalo)->r_sun;
125  for ( size_t i = 0; i<r.size(); i++ )
126  {
127  rho[i] = profile->bind("r")->eval(r[i]);
128  }
129  BEreq::set_halo_profile(0, r, rho, byVal(dist));
130  result = true;
131  }
132 
133 
136  void lnL_FermiLATdwarfs_gamLike(double &result)
137  {
138  using namespace Pipes::lnL_FermiLATdwarfs_gamLike;
139 
140  double fraction = *Dep::RD_fraction;
141  int mode = 0;
142  result = 0;
143 
145  std::string version = runOptions->getValueOrDef<std::string>("pass8", "version");
146  if ( version == "pass8" ) mode = 1;
147  else if ( version == "pass7" ) mode = 0;
148  else DarkBit_error().raise(LOCAL_INFO, "Fermi LAT dwarf likelihood version unknown.");
149 
150  // from 0.5 to 500 GeV
151  std::vector<double> x = daFunk::logspace(-0.301, 2.699, 100);
152  x = daFunk::augmentSingl(x, (*Dep::GA_AnnYield)->set("v",0));
153 
154  std::vector<double> y = ((*Dep::GA_AnnYield)/8./M_PI*fraction*fraction)->
155  set("v", 0)->bind("E")->vect(x);
156 
157  result = BEreq::lnL(byVal(mode), x, y);
158 
159  logger() << LogTags::debug << "GamLike dSph likelihood is lnL = " << result << EOM;
160  }
161 
162  void lnL_HESSGC_gamLike(double &result)
163  {
164  using namespace Pipes::lnL_HESSGC_gamLike;
165 
166  double fraction = *Dep::RD_fraction;
167  int mode = 0;
168  result = 0;
169 
171  std::string version = runOptions->getValueOrDef<std::string>("spectral_externalJ", "version");
172  if ( version == "integral_fixedJ" ) mode = 6;
173  else if ( version == "spectral_fixedJ" ) mode = 7;
174  else if ( version == "integral_externalJ" ) mode = 9;
175  else if ( version == "spectral_externalJ" ) mode = 10;
176  else DarkBit_error().raise(LOCAL_INFO, "HESS GC likelihood version unknown.");
177 
178  // from 230(265) GeV to 30 TeV
179  std::vector<double> x = daFunk::logspace(2.36, 4.48, 100);
180  x = daFunk::augmentSingl(x, (*Dep::GA_AnnYield)->set("v",0));
181  std::vector<double> y = ((*Dep::GA_AnnYield)/8./M_PI*fraction*fraction)->
182  set("v", 0)->bind("E")->vect(x);
183 
184  result = BEreq::lnL(byVal(mode), x, y);
185 
186  logger() << LogTags::debug << "GamLike HESS GC likelihood is lnL = " << result << EOM;
187  }
188 
189  void lnL_CTAGC_gamLike(double &result)
190  {
191  using namespace Pipes::lnL_CTAGC_gamLike;
192 
193  double fraction = *Dep::RD_fraction;
194  result = 0;
195 
196  // from 25 GeV to 10 TeV
197  std::vector<double> x = daFunk::logspace(1.39, 4.00, 100);
198  x = daFunk::augmentSingl(x, (*Dep::GA_AnnYield)->set("v",0));
199  std::vector<double> y = ((*Dep::GA_AnnYield)/8./M_PI*fraction*fraction)->
200  set("v", 0)->bind("E")->vect(x);
201 
202  result = BEreq::lnL(5, x, y);
203 
204  logger() << LogTags::debug << "GamLike CTA GC likelihood is lnL = " << result << EOM;
205  }
206 
209  void lnL_FermiGC_gamLike(double &result)
210  {
211  using namespace Pipes::lnL_FermiGC_gamLike;
212 
213  double fraction = *Dep::RD_fraction;
214  int mode = 0;
215  result = 0;
216 
218  std::string version = runOptions->getValueOrDef<std::string>("externalJ", "version");
219  if ( version == "fixedJ" ) mode = 2;
220  else if ( version == "margJ" ) mode = 3;
221  else if ( version == "margJ_HEP" ) mode = 4;
222  else if ( version == "externalJ" ) mode = 8;
223  else DarkBit_error().raise(LOCAL_INFO, "Fermi LAT GC likelihood version unknown.");
224 
225  // from 0.3 to 500 GeV
226  std::vector<double> x = daFunk::logspace(-0.523, 2.699, 100);
227  x = daFunk::augmentSingl(x, (*Dep::GA_AnnYield)->set("v",0));
228  std::vector<double> y = ((*Dep::GA_AnnYield)/8./M_PI*fraction*fraction)->
229  set("v", 0)->bind("E")->vect(x);
230 
231  result = BEreq::lnL(byVal(mode), x, y);
232 
233  logger() << LogTags::debug << "GamLike Fermi GC likelihood is lnL = " << result << EOM;
234  }
235 
240  void lnL_oh2_Simple(double &result)
241  {
242  using namespace Pipes::lnL_oh2_Simple;
243  double oh2_theory = *Dep::RD_oh2;
246  double oh2_theoryerr = oh2_theory*runOptions->getValueOrDef<double>(0.05, "oh2_fractional_theory_err");
248  double oh2_obs = runOptions->getValueOrDef<double>(0.1188, "oh2_obs");
250  double oh2_obserr = runOptions->getValueOrDef<double>(0.001, "oh2_obserr");
252  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
253  result = Stats::gaussian_loglikelihood(oh2_theory, oh2_obs, oh2_theoryerr, oh2_obserr, profile);
254  logger() << LogTags::debug << "lnL_oh2_Simple yields " << result << EOM;
255  }
256 
261  void lnL_oh2_upperlimit(double &result)
262  {
263  using namespace Pipes::lnL_oh2_upperlimit;
264  double oh2_theory = *Dep::RD_oh2;
267  double oh2_theoryerr = oh2_theory*runOptions->getValueOrDef<double>(0.05, "oh2_fractional_theory_err");
269  double oh2_obs = runOptions->getValueOrDef<double>(0.1188, "oh2_obs");
271  double oh2_obserr = runOptions->getValueOrDef<double>(0.001, "oh2_obserr");
273  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
274  result = Stats::gaussian_upper_limit(oh2_theory, oh2_obs, oh2_theoryerr, oh2_obserr, profile);
275  logger() << LogTags::debug << "lnL_oh2_upperlimit yields " << result << EOM;
276  }
277 
278 
284 
285  void lnL_sigmas_sigmal(double &result)
286  {
287  using namespace Pipes::lnL_sigmas_sigmal;
288  double sigmas = *Param["sigmas"];
289  double sigmal = *Param["sigmal"];
291  double sigmas_obs = runOptions->getValueOrDef<double>(43., "sigmas_obs");
293  double sigmas_obserr = runOptions->getValueOrDef<double>(8., "sigmas_obserr");
295  double sigmal_obs = runOptions->getValueOrDef<double>(58., "sigmal_obs");
297  double sigmal_obserr = runOptions->getValueOrDef<double>(9., "sigmal_obserr");
299  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
300 
301  result = Stats::gaussian_loglikelihood(sigmas, sigmas_obs, 0, sigmas_obserr, profile)
302  + Stats::gaussian_loglikelihood(sigmal, sigmal_obs, 0, sigmal_obserr, profile);
303  logger() << LogTags::debug << "lnL for SI nuclear parameters is " << result << EOM;
304  }
305 
315  void lnL_deltaq(double &result)
316  {
317  using namespace Pipes::lnL_deltaq;
318  double deltad = *Param["deltad"];
319  double deltau = *Param["deltau"];
320  double deltas = *Param["deltas"];
321  double a3 = deltau - deltad;
322  double a8 = deltau + deltad - 2*deltas;
323 
325  double a3_obs = runOptions->getValueOrDef<double>(1.2723, "a3_obs");
327  double a3_obserr = runOptions->getValueOrDef<double>(0.0023, "a3_obserr");
329  double a8_obs = runOptions->getValueOrDef<double>(0.585, "a8_obs");
331  double a8_obserr = runOptions->getValueOrDef<double>(0.025, "a8_obserr");
333  double deltas_obs = runOptions->getValueOrDef<double>(-0.09, "deltas_obs");
335  double deltas_obserr = runOptions->getValueOrDef<double>(0.03, "deltas_obserr");
337  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
338 
339  result = Stats::gaussian_loglikelihood(a3, a3_obs, 0, a3_obserr, profile) +
340  Stats::gaussian_loglikelihood(a8, a8_obs, 0, a8_obserr, profile) +
341  Stats::gaussian_loglikelihood(deltas, deltas_obs, 0, deltas_obserr, profile);
342  }
343 
348 
349  void lnL_rho0_lognormal(double &result)
350  {
351  using namespace Pipes::lnL_rho0_lognormal;
352  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
353  double rho0 = LocalHaloParameters.rho0;
355  double rho0_obs = runOptions->getValueOrDef<double>(.4, "rho0_obs");
357  double rho0_obserr = runOptions->getValueOrDef<double>(.15, "rho0_obserr");
359  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
360 
361  result = Stats::lognormal_loglikelihood(rho0, rho0_obs, 0.,
362  rho0_obserr, profile);
363  logger() << LogTags::debug << "lnL_rho0 yields " << result << EOM;
364  }
365 
366  void lnL_vrot_gaussian(double &result)
367  {
368  using namespace Pipes::lnL_vrot_gaussian;
369  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
370  double vrot = LocalHaloParameters.vrot;
372  double vrot_obs = runOptions->getValueOrDef<double>(235, "vrot_obs");
374  double vrot_obserr = runOptions->getValueOrDef<double>(20, "vrot_obserr");
376  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
377  result = Stats::gaussian_loglikelihood(vrot, vrot_obs, 0., vrot_obserr, profile);
378  logger() << LogTags::debug << "lnL_vrot yields " << result << EOM;
379  }
380 
381  void lnL_v0_gaussian(double &result)
382  {
383  using namespace Pipes::lnL_v0_gaussian;
384  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
385  double v0 = LocalHaloParameters.v0;
387  double v0_obs = runOptions->getValueOrDef<double>(235, "v0_obs");
389  double v0_obserr = runOptions->getValueOrDef<double>(20, "v0_obserr");
391  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
392  result = Stats::gaussian_loglikelihood(v0, v0_obs, 0., v0_obserr, profile);
393  logger() << LogTags::debug << "lnL_v0 yields " << result << EOM;
394  }
395 
396  void lnL_vesc_gaussian(double &result)
397  {
398  using namespace Pipes::lnL_vesc_gaussian;
399  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
400  double vesc = LocalHaloParameters.vesc;
402  double vesc_obs = runOptions->getValueOrDef<double>(550, "vesc_obs");
404  double vesc_obserr = runOptions->getValueOrDef<double>(35, "vesc_obserr");
406  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
407  result = Stats::gaussian_loglikelihood(vesc, vesc_obs, 0., vesc_obserr, profile);
408  logger() << LogTags::debug << "lnL_vesc yields " << result << EOM;
409  }
410 
412  void dump_GammaSpectrum(double &result)
413  {
414  using namespace Pipes::dump_GammaSpectrum;
415  daFunk::Funk spectrum = (*Dep::GA_AnnYield)->set("v", 0.);
416  // Option filename<string>: Filename for gamma-ray spectrum dump
417  // (default: dNdE.dat)
418  std::string filename = runOptions->getValueOrDef<std::string>(
419  "dNdE.dat", "filename");
420  logger() << "FILENAME for gamma dump: " << filename << EOM;
421  std::ofstream myfile (filename);
422  if (myfile.is_open())
423  {
424  for (int i = 0; i<=1200; i++)
425  {
426  double energy = pow(10., i/200. - 4.);
427 
428  myfile << energy << " " << spectrum->bind("E")->eval(energy) << "\n";
429  }
430  myfile.close();
431  }
432  result = 0.;
433  }
434  }
435 }
void lnL_FermiGC_gamLike(double &result)
Fermi LAT galactic center likelihoods, using gamLike backend.
error & DarkBit_error()
void lnL_sigmas_sigmal(double &result)
Likelihoods for spin independent nuclear parameters.
void lnL_oh2_upperlimit(double &result)
Likelihood for cosmological relic density constraints, implemented as an upper limit only Default dat...
START_MODEL r_sun
#define LOCAL_INFO
Definition: local_info.hpp:34
double gaussian_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is Gaussian...
Definition: statistics.cpp:133
void dump_GammaSpectrum(double &result)
Helper function to dump gamma-ray spectra.
START_MODEL v0
void lnL_deltaq(double &result)
Likelihoods for spin dependent nuclear parameters.
void lnL_v0_gaussian(double &result)
START_MODEL vesc
void lnL_oh2_Simple(double &result)
Likelihood for cosmological relic density constraints.
Declarations of statistical utilities.
double gaussian_upper_limit(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a gaussian log-likelihood for an upper limit.
Definition: statistics.cpp:184
void lnL_FermiLATdwarfs_gamLike(double &result)
Fermi LAT dwarf likelihoods, using gamLike backend.
std::vector< double > logspace(double x0, double x1, unsigned int n)
Definition: daFunk.hpp:186
void lnL_rho0_lognormal(double &result)
Likelihoods for halo parameters.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
void set_gamLike_GC_halo(bool &result)
Fermi LAT dwarf likelihoods, based on arXiv:1108.2914.
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 deltau
START_MODEL dNur_CMB r
void lnL_vesc_gaussian(double &result)
START_MODEL vrot
void lnL_vrot_gaussian(double &result)
Rollcall header for module DarkBit.
void lnL_HESSGC_gamLike(double &result)
double pow(const double &a)
Outputs a^i.
double lognormal_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is log-norm...
Definition: statistics.cpp:154
std::vector< double > augmentSingl(const std::vector< double > &xgrid, Funk f, int N=100, double sigma=3)
Definition: daFunk.hpp:1714
T byVal(T t)
Redirection function to turn an lvalue into an rvalue, so that it is correctly passed by value when d...
def profile(file_name, frac_error=0.1, min_=0., max_=1., log_normal=True)
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
void lnL_CTAGC_gamLike(double &result)
TODO: see if we can use this one:
Definition: Analysis.hpp:33