gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
SunNeutrinos.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
25 
29 
30 //#define DARKBIT_DEBUG
31 
32 namespace Gambit
33 {
34 
35  namespace DarkBit
36  {
37 
39  //
40  // Neutrino telescope likelihoods and observables
41  //
43 
47  void capture_rate_Sun_const_xsec_DS5(double &result)
48  {
50 
51  if (BEreq::cap_Sun_v0q0_isoscalar.origin()=="DarkSUSY")
52  if(!(*Dep::DarkSUSY5_PointInit_LocalHalo))
53  DarkBit_error().raise(LOCAL_INFO,"DarkSUSY halo model not initialized!");
54 
55  // When calculating the solar capture rate, DarkSUSY assumes that the
56  // proton and neutron scattering cross-sections are the same; we
57  // assume that whichever backend has been hooked up here does so too.
58 
59  result = BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp, *Dep::sigma_SI_p, *Dep::sigma_SD_p);
60 
61  }
62 
66  void capture_rate_Sun_const_xsec(double &result)
67  {
68  using namespace Pipes::capture_rate_Sun_const_xsec;
69 
70  if (BEreq::cap_Sun_v0q0_isoscalar.origin()=="DarkSUSY")
71  if(!(*Dep::DarkSUSY_PointInit_LocalHalo))
72  DarkBit_error().raise(LOCAL_INFO,"DarkSUSY halo model not initialized!");
73 
74  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
75  double rho0 = LocalHaloParameters.rho0;
76  double rho0_eff = (*Dep::RD_fraction)*rho0;
77 
78  // When calculating the solar capture rate, DarkSUSY assumes that the
79  // proton and neutron scattering cross-sections are the same; we
80  // assume that whichever backend has been hooked up here does so too.
81 
82  result = BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp, rho0_eff, *Dep::sigma_SI_p, *Dep::sigma_SD_p);
83 
84  }
85 
88  {
90  double resultSD;
91  double resultSI;
92  double maxcap;
93 
94  BEreq::cap_sun_saturation(*Dep::mwimp,maxcap);
95  BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp,*Dep::sigma_SD_p,*Dep::sigma_SI_p,resultSD,resultSI);
96  result = resultSI + resultSD;
97 
98  if (maxcap < result)
99  {
100  result = maxcap;
101  }
102 
103  }
104 
109  void capture_rate_Sun_vnqn(double &result)
110  {
111  using namespace Pipes::capture_rate_Sun_vnqn;
112 
113  double resultSD;
114  double resultSI;
115  double capped;
116  int qpow;
117  int vpow;
118  const int nelems = 29;
119  double maxcap;
120 
121  BEreq::cap_sun_saturation(*Dep::mwimp,maxcap);
122 
123  resultSI = 0e0;
124  resultSD = 0e0;
125  typedef map_intpair_dbl::const_iterator it_type;
126 
127  //Spin-dependent:
128  for(it_type iterator = Dep::sigma_SD_p->begin();
129  iterator != Dep::sigma_SD_p->end();
130  iterator++)
131  {
132  //don't capture anything if cross section is zero or all the DM is already capped
133  if((iterator->second > 1e-90) && (resultSD < maxcap))
134  {
135  qpow = (iterator->first).first/2 ;
136  vpow = (iterator->first).second/2;
137 
138  //Capture
139  BEreq::cap_Sun_vnqn_isoscalar(*Dep::mwimp,iterator->second,1,qpow,vpow,capped);
140  resultSD = resultSD+capped;
141  }
142  }
143 
144  //Spin independent:
145  for(it_type iterator = Dep::sigma_SI_p->begin();
146  iterator != Dep::sigma_SI_p->end();
147  iterator++)
148  {
149  if((iterator->second > 1e-90) && (resultSI+resultSD < maxcap))
150  {
151  qpow = (iterator->first).first/2 ;
152  vpow = (iterator->first).second/2;
153 
154  //Capture
155  BEreq::cap_Sun_vnqn_isoscalar(*Dep::mwimp,iterator->second,nelems,qpow,vpow,capped);
156  resultSI = resultSI+capped;
157  }
158  }
159  result = resultSI+resultSD;
160 
161  logger() << "Capgen captured: SI: " << resultSI << " SD: " << resultSD << " total: " << result << "max = " << maxcap << "\n" << EOM;
162 
163  // If capture is above saturation, return saturation value.
164  if (maxcap < result)
165  {
166  result = maxcap;
167  }
168 
169  //cout << "capture rate via capture_rate_Sun_vnqn = " << result << "\n";
170 
171  }
172 
176  void equilibration_time_Sun(double &result)
177  {
178  using namespace Pipes::equilibration_time_Sun;
179 
180  double sigmav = 0;
181  double T_Sun_core = 1.35e-6; // Sun's core temperature (GeV)
182 
183  std::string DMid = *Dep::DarkMatter_ID;
184  std::string DMbarid = *Dep::DarkMatterConj_ID;
185 
186  // Make sure that we're not trying to work with decaying DM.
187  const TH_Process* p = Dep::TH_ProcessCatalog->find(DMid, DMbarid);
188  if (p == NULL) DarkBit_error().raise(LOCAL_INFO, "Sorry, decaying DM is not supported yet by the DarkBit neutrino routines.");
189  TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMbarid);
190 
191  // Add all the regular channels
192  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
193  it != annProc.channelList.end(); ++it)
194  {
195  if ( it->nFinalStates == 2 )
196  {
197  // (sv)(v = sqrt(2T/mDM)) for two-body final state
198  sigmav += it->genRate->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
199  }
200  }
201 
202  // Add invisible contributions
203  sigmav += annProc.genRateMisc->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
204 
205  double ca = sigmav/6.6e28 * pow(*Dep::mwimp/20.0, 1.5);
206  // Scale the annihilation rate down by a factor of two if the DM is not self-conjugate
207  if (not (*Dep::TH_ProcessCatalog).getProcess(*Dep::DarkMatter_ID, *Dep::DarkMatterConj_ID).isSelfConj) ca *= 0.5;
208  result = pow(*Dep::capture_rate_Sun * ca, -0.5);
209 
210  // std::cout << "v = " << sqrt(2.0*T_Sun_core/(*Dep::mwimp)) << " and sigmav inside equilibration_time_Sun = " << sigmav << std::endl;
211  // std::cout << "capture_rate_Sun inside equilibration_time_Sun = " << *Dep::capture_rate_Sun << std::endl;
212  }
213 
215  void annihilation_rate_Sun(double &result)
216  {
217  using namespace Pipes::annihilation_rate_Sun;
218  double tt_sun = 1.5e17 / *Dep::equilibration_time_Sun;
219  result = *Dep::capture_rate_Sun * 0.5 * pow(tanh(tt_sun),2.0);
220  }
221 
224  {
225 
226  using namespace Pipes::nuyield_from_DS;
227  double annihilation_bf[29];
228  double Higgs_decay_BFs_neutral[29][3];
229  double Higgs_decay_BFs_charged[15];
230  double Higgs_masses_neutral[3];
231  double Higgs_mass_charged;
232 
233  // Set annihilation branching fractions
234  std::string DMid = *Dep::DarkMatter_ID;
235  std::string DMbarid = *Dep::DarkMatterConj_ID;
236  TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMbarid);
237  std::vector< std::vector<str> > neutral_channels = BEreq::get_DS_neutral_h_decay_channels();
238  // the missing channel
239  const std::vector<str> adhoc_chan = initVector<str>("W-", "H+");
240 
241  for (int i=0; i<29; i++)
242  {
243  const TH_Channel* channel = annProc.find(neutral_channels[i]);
244 
245  if (channel == NULL or i == 26) // Channel 26 has not been implemented in DarkSUSY.
246  {
247  annihilation_bf[i] = 0.;
248  }
249  else
250  {
251  annihilation_bf[i] = channel->genRate->bind("v")->eval(0.);
252  if (i == 10) // Add W- H+ for this channel
253  {
254  channel = annProc.find(adhoc_chan);
255  if (channel == NULL) DarkBit_error().raise(LOCAL_INFO,
256  "W+H- exists in process catalog but not W-H+."
257  " That's some suspiciously severe CP violation yo.");
258  annihilation_bf[i] += channel->genRate->bind("v")->eval(0.);
259  }
260  annihilation_bf[i] /= *Dep::sigmav;
261 
262  // Check that having this channel turned on makes sense at all.
263  #ifdef DARKBIT_DEBUG
264  double mtot = 0;
265  cout << "Particles and masses in DM annihilation final state: " << endl;
266  for (auto p = neutral_channels[i].begin(); p != neutral_channels[i].end(); ++p)
267  {
268  double m = Dep::TH_ProcessCatalog->getParticleProperty(*p).mass;
269  cout << " " << *p << " " << m << endl;
270  mtot += m;
271  }
272  cout << "Sqrt(s) vs total mass of final states: " << 2 * *Dep::mwimp << " vs. " << mtot << endl;
273  cout << "Branching fraction in v=0 limit: " << annihilation_bf[i] << endl << endl;
274  if (mtot > 2 * *Dep::mwimp and annihilation_bf[i] > 0.0)
275  DarkBit_error().raise(LOCAL_INFO, "Channel is open in process catalog but should not be kinematically allowed.");
276  #endif
277  }
278 
279  }
280 
281  // Set Higgs masses
282  if (Dep::TH_ProcessCatalog->hasParticleProperty("h0_1"))
283  Higgs_masses_neutral[1] = Dep::TH_ProcessCatalog->getParticleProperty("h0_1").mass;
284  else
285  DarkBit_error().raise(LOCAL_INFO, "No SM-like Higgs in ProcessCatalog!");
286  Higgs_masses_neutral[0] =
287  Dep::TH_ProcessCatalog->hasParticleProperty("h0_2") ?
288  Dep::TH_ProcessCatalog->getParticleProperty("h0_2").mass : 0.;
289  Higgs_masses_neutral[2] =
290  Dep::TH_ProcessCatalog->hasParticleProperty("A0") ?
291  Dep::TH_ProcessCatalog->getParticleProperty("A0").mass : 0.;
292  Higgs_mass_charged =
293  Dep::TH_ProcessCatalog->hasParticleProperty("H+") ?
294  Dep::TH_ProcessCatalog->getParticleProperty("H+").mass : 0.;
295 
296  // Find out which Higgs exist and have decay data in the process
297  // catalog.
298  const TH_Process* h0_decays[3];
299  h0_decays[0] = Dep::TH_ProcessCatalog->find("h0_2");
300  h0_decays[1] = Dep::TH_ProcessCatalog->find("h0_1");
301  h0_decays[2] = Dep::TH_ProcessCatalog->find("A0");
302  const TH_Process* Hplus_decays = Dep::TH_ProcessCatalog->find("H+");
303  const TH_Process* Hminus_decays = Dep::TH_ProcessCatalog->find("H-");
304  if (Hplus_decays != NULL and Hminus_decays == NULL) DarkBit_error().raise(
305  LOCAL_INFO, "H+ decays exist in process catalog but not H-.");
306  if (Hplus_decays == NULL and Hminus_decays != NULL) DarkBit_error().raise(
307  LOCAL_INFO, "H- decays exist in process catalog but not H+.");
308 
309  // Set the neutral Higgs decay branching fractions
310  for (int i=0; i<3; i++) // Loop over the three neutral Higgs
311  {
312 
313  // If this Higgs exists, set its decay properties.
314  if (h0_decays[i] != NULL)
315  {
316 
317  // Get the total decay width, for normalising partial widths to BFs.
318  // TODO: Replace when BFs become directly available.
319  double totalwidth = 0.0;
320  for (std::vector<TH_Channel>::const_iterator
321  it = h0_decays[i]->channelList.begin();
322  it != h0_decays[i]->channelList.end(); ++it)
323  {
324  // decay width in GeV for two-body final state
325  if ( it->nFinalStates == 2 ) totalwidth +=
326  it->genRate->bind()->eval();
327  }
328 
329  std::vector<str> neutral_channel;
330  // Loop over the decay channels for neutral scalars
331  for (int j=0; j<29; j++)
332  {
333  neutral_channel.clear();
334  neutral_channel.push_back(DarkBit_utils::str_flav_to_mass((neutral_channels[j])[0]));
335  neutral_channel.push_back(DarkBit_utils::str_flav_to_mass((neutral_channels[j])[1]));
336  const TH_Channel* channel = h0_decays[i]->find(neutral_channel);
337  // If this Higgs can decay into this channel, set the BF.
338  if (channel != NULL)
339  {
340  Higgs_decay_BFs_neutral[j][i] = channel->genRate->bind()->eval();
341  if (j == 10) // Add W- H+ for this channel
342  {
343  channel = h0_decays[i]->find(adhoc_chan);
344  if (channel == NULL) DarkBit_error().raise(LOCAL_INFO,
345  "W+H- exists in process catalog but not W-H+."
346  " That's some suspiciously severe CP violation yo.");
347  Higgs_decay_BFs_neutral[j][i]
348  += channel->genRate->bind()->eval();
349  }
350  // This channel has not been implemented in DarkSUSY.
351  if (j == 26) Higgs_decay_BFs_neutral[j][i] = 0.;
352  Higgs_decay_BFs_neutral[j][i] /= totalwidth;
353  }
354  else
355  {
356  Higgs_decay_BFs_neutral[j][i] = 0.;
357  }
358  }
359 
360  }
361 
362  else
363  {
364  // Loop over the decay channels for neutral scalars, setting all BFs
365  // for this Higgs to zero.
366  for (int j=0; j<29; j++) Higgs_decay_BFs_neutral[j][i] = 0.;
367  }
368 
369  }
370 
371  // If they exist, set the charged Higgs decay branching fractions
372  // (DarkSUSY assumes that H+/H- decays are CP-invariant)
373  if (Hplus_decays != NULL)
374  {
375 
376  // Define the charged Higgs decay channels
377  std::vector< std::vector<str> > charged_channels = BEreq::get_DS_charged_h_decay_channels();
378 
379  // Get the total decay width, for normalising partial widths to BFs.
380  // TODO: Replace when BFs become directly available.
381  double totalwidth = 0.0;
382  for (std::vector<TH_Channel>::const_iterator
383  it = Hplus_decays->channelList.begin();
384  it != Hplus_decays->channelList.end(); ++it)
385  {
386  // decay width in GeV for two-body final state
387  if (it->nFinalStates == 2) totalwidth += it->genRate->bind()->eval();
388  }
389 
390  std::vector<str> charged_channel;
391  // Loop over the decay channels for charged scalars
392  for (int j=0; j<15; j++)
393  {
394  charged_channel.clear();
395  charged_channel.push_back(DarkBit_utils::str_flav_to_mass((charged_channels[j])[0]));
396  charged_channel.push_back(DarkBit_utils::str_flav_to_mass((charged_channels[j])[1]));
397  const TH_Channel* channel = Hplus_decays->find(charged_channel);
398  // If this Higgs can decay into this channel, set the BF.
399  if (channel != NULL)
400  {
401  Higgs_decay_BFs_charged[j] = channel->genRate->bind()->eval();
402  Higgs_decay_BFs_charged[j] /= totalwidth;
403  }
404  else
405  {
406  Higgs_decay_BFs_charged[j] = 0.;
407  }
408  }
409 
410  }
411 
412  else
413  {
414  // Loop over the decay channels for charged scalars, setting all BFs
415  // for this Higgs to zero.
416  for (int j=0; j<15; j++) Higgs_decay_BFs_charged[j] = 0.;
417  }
418 
419  // Debug output
420  #ifdef DARKBIT_DEBUG
421  for (int j=0; j<29; j++)
422  {
423  cout<< "annihilation bfs: " << j << " " << annihilation_bf[j] << endl;
424  }
425  cout<< endl;
426  for (int i=0; i<3; i++)
427  {
428  for (int j=0; j<29; j++)
429  {
430  cout<< "higgs neutral bfs: " << i << " " << j << " "
431  << Higgs_decay_BFs_neutral[j][i] << endl;
432  }
433  }
434  cout<< endl;
435  for (int j=0; j<15; j++)
436  {
437  cout<< "higgs charged bfs: " << j << " "
438  << Higgs_decay_BFs_charged[j] << endl;
439  }
440  cout<< endl;
441  for (int j=0; j<3; j++)
442  {
443  cout<< "higgs masses neutral: " << j << " "
444  << Higgs_masses_neutral[j] << endl;
445  }
446  cout<< endl;
447  cout<< "higgs charged mass: " << Higgs_mass_charged << endl;
448  cout<< endl;
449  cout<< "*Dep::mwimp: " << *Dep::mwimp << endl;
450  cout<< endl;
451  cout<< "*Dep::sigmav: " << *Dep::sigmav << endl;
452  cout<< endl;
453  #endif
454 
455  // Set up DarkSUSY to do neutrino yields for this particular WIMP
456  BEreq::DS_nuyield_setup(annihilation_bf, Higgs_decay_BFs_neutral,
457  Higgs_decay_BFs_charged, Higgs_masses_neutral,
458  Higgs_mass_charged, *Dep::mwimp);
459 
460  // Hand back the pointer to the DarkSUSY neutrino yield function
461  result.pointer = BEreq::nuyield.pointer();
462 
463  //TODO: change below to >= when version numbers are available as ints
464  // Treat the yield function as threadsafe only if the loaded version of DarkSUSY supports it.
465  result.threadsafe = (BEreq::nuyield.version() == "5.1.3");
466 
467  // Avoid OpenMP with gfortran 6.x and later, as the implementation seems to be buggy (causes stack overflows).
468  #ifdef __GNUC__
469  #define GCC_VERSION (__GNUC__ * 10000 \
470  + __GNUC_MINOR__ * 100 \
471  + __GNUC_PATCHLEVEL__)
472  #if GCC_VERSION > 60000
473  result.threadsafe = false;
474  #endif
475  #undef GCC_VERSION
476  #endif
477 
478  }
479 
488 
491  void IC22_full(nudata &result)
492  {
493  using namespace Pipes::IC22_full;
494  static bool first = true;
495  const double bgloglike = -808.4581;
496  double sigpred, bgpred, lnLike, pval;
497  int totobs;
498  char experiment[300] = "IC-22";
499  void* context = NULL;
500  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
502  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
503  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
504  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
505  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
509  result.signal = sigpred;
510  result.bg = bgpred;
511  result.nobs = totobs;
512  result.loglike = lnLike;
513  result.pvalue = pval;
514  if (first)
515  {
516  result.bgloglike = bgloglike;
517  first = false;
518  }
519  }
520 
523  void IC79WH_full(nudata &result)
524  {
525  static bool first = true;
526  using namespace Pipes::IC79WH_full;
527  const double bgloglike = -11874.8689;
528  double sigpred, bgpred, lnLike, pval;
529  int totobs;
530  char experiment[300] = "IC-79 WH";
531  void* context = NULL;
532  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
534  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
535  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
536  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
537  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
541  result.signal = sigpred;
542  result.bg = bgpred;
543  result.nobs = totobs;
544  result.loglike = lnLike;
545  result.pvalue = pval;
546  if (first)
547  {
548  result.bgloglike = bgloglike;
549  first = false;
550  }
551  }
552 
555  void IC79WL_full(nudata &result)
556  {
557  static bool first = true;
558  using namespace Pipes::IC79WL_full;
559  const double bgloglike = -1813.4503;
560  double sigpred, bgpred, lnLike, pval;
561  int totobs;
562  char experiment[300] = "IC-79 WL";
563  void* context = NULL;
564  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
566  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
567  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
568  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
569  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
573  result.signal = sigpred;
574  result.bg = bgpred;
575  result.nobs = totobs;
576  result.loglike = lnLike;
577  result.pvalue = pval;
578  if (first)
579  {
580  result.bgloglike = bgloglike;
581  first = false;
582  }
583  }
584 
587  void IC79SL_full(nudata &result)
588  {
589  static bool first = true;
590  using namespace Pipes::IC79SL_full;
591  const double bgloglike = -5015.6474;
592  double sigpred, bgpred, lnLike, pval;
593  int totobs;
594  char experiment[300] = "IC-79 SL";
595  void* context = NULL;
596  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
598  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
599  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
600  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
601  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
605  result.signal = sigpred;
606  result.bg = bgpred;
607  result.nobs = totobs;
608  result.loglike = lnLike;
609  result.pvalue = pval;
610  if (first)
611  {
612  result.bgloglike = bgloglike;
613  first = false;
614  }
615  }
617 
620  void IC22_signal (double &result) {
621  result = Pipes::IC22_signal ::Dep::IC22_data->signal; }
622  void IC22_bg (double &result) {
623  result = Pipes::IC22_bg ::Dep::IC22_data->bg; }
624  void IC22_nobs (int &result) {
625  result = Pipes::IC22_nobs ::Dep::IC22_data->nobs; }
626  void IC22_loglike(double &result) {
627  result = Pipes::IC22_loglike::Dep::IC22_data->loglike; }
628  void IC22_bgloglike(double &result) {
629  result = Pipes::IC22_bgloglike::Dep::IC22_data->bgloglike; }
630  void IC22_pvalue (double &result) {
631  result = Pipes::IC22_pvalue ::Dep::IC22_data->pvalue; }
633 
636  void IC79WH_signal (double &result) {
637  result = Pipes::IC79WH_signal ::Dep::IC79WH_data->signal; }
638  void IC79WH_bg (double &result) {
639  result = Pipes::IC79WH_bg ::Dep::IC79WH_data->bg; }
640  void IC79WH_nobs (int &result) {
641  result = Pipes::IC79WH_nobs ::Dep::IC79WH_data->nobs; }
642  void IC79WH_loglike(double &result) {
643  result = Pipes::IC79WH_loglike::Dep::IC79WH_data->loglike; }
644  void IC79WH_bgloglike(double &result) {
645  result = Pipes::IC79WH_bgloglike::Dep::IC79WH_data->bgloglike; }
646  void IC79WH_pvalue (double &result) {
647  result = Pipes::IC79WH_pvalue ::Dep::IC79WH_data->pvalue; }
649 
652  void IC79WL_signal (double &result) {
653  result = Pipes::IC79WL_signal ::Dep::IC79WL_data->signal; }
654  void IC79WL_bg (double &result) {
655  result = Pipes::IC79WL_bg ::Dep::IC79WL_data->bg; }
656  void IC79WL_nobs (int &result) {
657  result = Pipes::IC79WL_nobs ::Dep::IC79WL_data->nobs; }
658  void IC79WL_loglike(double &result) {
659  result = Pipes::IC79WL_loglike::Dep::IC79WL_data->loglike; }
660  void IC79WL_bgloglike(double &result) {
661  result = Pipes::IC79WL_bgloglike::Dep::IC79WL_data->bgloglike; }
662  void IC79WL_pvalue (double &result) {
663  result = Pipes::IC79WL_pvalue ::Dep::IC79WL_data->pvalue; }
665 
668  void IC79SL_signal (double &result) {
669  result = Pipes::IC79SL_signal ::Dep::IC79SL_data->signal; }
670  void IC79SL_bg (double &result) {
671  result = Pipes::IC79SL_bg ::Dep::IC79SL_data->bg; }
672  void IC79SL_nobs (int &result) {
673  result = Pipes::IC79SL_nobs ::Dep::IC79SL_data->nobs; }
674  void IC79SL_loglike(double &result) {
675  result = Pipes::IC79SL_loglike::Dep::IC79SL_data->loglike; }
676  void IC79SL_bgloglike(double &result) {
677  result = Pipes::IC79SL_bgloglike::Dep::IC79SL_data->bgloglike; }
678  void IC79SL_pvalue (double &result) {
679  result = Pipes::IC79SL_pvalue ::Dep::IC79SL_data->pvalue; }
681 
683  void IC79_loglike(double &result)
684  {
685  using namespace Pipes::IC79_loglike;
689  if (result > 0.0) result = 0.0;
690  }
691 
693  void IC_loglike(double &result)
694  {
695  using namespace Pipes::IC_loglike;
700  if (result > 0.0) result = 0.0;
701 #ifdef DARKBIT_DEBUG
702  cout << "IC likelihood: " << result << endl;
703  cout << "IC79SL contribution: " << *Dep::IC79SL_loglike - *Dep::IC79SL_bgloglike << endl;
704  cout << "IC79WL contribution: " << *Dep::IC79WL_loglike - *Dep::IC79WL_bgloglike << endl;
705  cout << "IC79WH contribution: " << *Dep::IC79WH_loglike - *Dep::IC79WH_bgloglike << endl;
706  cout << "IC22 contribution: " << *Dep::IC22_loglike - *Dep::IC22_bgloglike << endl;
707 #endif
708  }
709 
712  {
714 
715  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
716 
717  double rho0 = LocalHaloParameters.rho0;
718  double rho0_eff = (*Dep::RD_fraction)*rho0;
719  double vrot = LocalHaloParameters.vrot;
720  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
721  double vesc = LocalHaloParameters.vesc;
723  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
724 
725  BEreq::dshmcom->rho0 = rho0;
726  BEreq::dshmcom->v_sun = vrot;
727  BEreq::dshmcom->v_earth = v_earth;
728  BEreq::dshmcom->rhox = rho0_eff;
729 
730  BEreq::dshmframevelcom->v_obs = vrot;
731 
732  BEreq::dshmisodf->vd_3d = vd_3d;
733  BEreq::dshmisodf->vgalesc = vesc;
734 
735  BEreq::dshmnoclue->vobs = vrot;
736 
738  << "Updating DarkSUSY halo parameters:" << std::endl
739  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
740  << " rho0_eff [GeV/cm^3] = " << rho0_eff << std::endl
741  << " v_sun [km/s] = " << vrot<< std::endl
742  << " v_earth [km/s] = " << v_earth << std::endl
743  << " v_obs [km/s] = " << vrot << std::endl
744  << " vd_3d [km/s] = " << vd_3d << std::endl
745  << " v_esc [km/s] = " << vesc << EOM;
746 
747  result = true;
748 
749  return;
750  }
751 
754  {
756 
757  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
758 
759  double rho0 = LocalHaloParameters.rho0;
760  double vrot = LocalHaloParameters.vrot;
761  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
762  double vesc = LocalHaloParameters.vesc;
764  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
765 
766  BEreq::dshmcom->rho0 = rho0;
767  BEreq::dshmcom->v_sun = vrot;
768  BEreq::dshmcom->v_earth = v_earth;
769 
770  BEreq::dshmframevelcom->v_obs = vrot;
771 
772  BEreq::dshmisodf->vd_3d = vd_3d;
773  BEreq::dshmisodf->vgalesc = vesc;
774 
775  BEreq::dshmnoclue->vobs = vrot;
776 
778  << "Updating DarkSUSY halo parameters:" << std::endl
779  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
780  << " v_sun [km/s] = " << vrot<< std::endl
781  << " v_earth [km/s] = " << v_earth << std::endl
782  << " v_obs [km/s] = " << vrot << std::endl
783  << " vd_3d [km/s] = " << vd_3d << std::endl
784  << " v_esc [km/s] = " << vesc << EOM;
785 
786  result = true;
787 
788  return;
789  }
790 
791  }
792 }
793 
std::vector< TH_Channel > channelList
List of channels.
error & DarkBit_error()
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
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)
void IC22_bgloglike(double &result)
void DarkSUSY_PointInit_LocalHalo_func(bool &result)
Function to set Local Halo Parameters in DarkSUSY (DS 6)
void IC79SL_bgloglike(double &result)
void IC79WH_nobs(int &result)
Neutrino telescope data container.
void IC79WL_bgloglike(double &result)
void IC79SL_signal(double &result)
79-string SL extractor module functions
void IC79WL_bg(double &result)
#define LOCAL_INFO
Definition: local_info.hpp:34
void IC79WL_pvalue(double &result)
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
void IC79WL_loglike(double &result)
void IC79SL_bg(double &result)
void nuyield_from_DS(nuyield_info &result)
Neutrino yield function pointer and setup.
void IC79SL_pvalue(double &result)
void annihilation_rate_Sun(double &result)
Annihilation rate of dark matter in the Sun (s^-1)
void IC79WH_loglike(double &result)
void IC22_signal(double &result)
22-string extractor module functions
void IC79SL_loglike(double &result)
void IC22_nobs(int &result)
Neutrino telescope yield info container.
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:544
START_MODEL vesc
void IC79SL_nobs(int &result)
void IC22_loglike(double &result)
void IC22_full(nudata &result)
Likelihood calculators for different IceCube event samples These functions all include the likelihood...
void capture_rate_Sun_vnqn(double &result)
Capture rate for v^n and q^n-dependent cross sections.
void IC_loglike(double &result)
Complete composite IceCube likelihood function.
daFunk::Funk genRate
Energy dependence of final state particles. Includes v_rel ("v") as last argument in case of annihila...
void IC79SL_full(nudata &result)
79-string IceCube SL sample: predicted signal and background counts, observed counts and likelihoods...
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 IC79WH_signal(double &result)
79-string WH extractor module functions
void IC22_pvalue(double &result)
const TH_Channel * find(std::vector< str >) const
Check for given channel. Return a pointer to it if found, NULL if not.
void IC79WH_bg(double &result)
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
Header file that includes all GAMBIT headers required for a module source file.
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
void IC79WL_nobs(int &result)
void IC79WL_signal(double &result)
79-string WL extractor module functions
void IC22_bg(double &result)
START_MODEL vrot
All data on a single annihilation or decay channel, e.g.
daFunk::Funk genRateMisc
Additional decay rate or sigmav (in addition to above channels)
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.
Rollcall header for module DarkBit.
void IC79WH_bgloglike(double &result)
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:564
void capture_rate_Sun_const_xsec_capgen(double &result)
Alternative to the darkSusy fct, using captn_specific from capgen instead.
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:481
warning & DarkBit_warning()
double pow(const double &a)
Outputs a^i.
std::string str_flav_to_mass(std::string flav)
nuyield_function_pointer pointer
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...
void IC79WH_pvalue(double &result)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Utility functions for DarkBit.