gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
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  TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMid);
185 
186  // Add all the regular channels
187  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
188  it != annProc.channelList.end(); ++it)
189  {
190  if ( it->nFinalStates == 2 )
191  {
192  // (sv)(v = sqrt(2T/mDM)) for two-body final state
193  sigmav += it->genRate->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
194  }
195  }
196 
197  // Add invisible contributions
198  sigmav += annProc.genRateMisc->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
199 
200  double ca = sigmav/6.6e28 * pow(*Dep::mwimp/20.0, 1.5);
201  // Scale the annihilation rate down by a factor of two if the DM is not self-conjugate
202  if (not (*Dep::TH_ProcessCatalog).getProcess(*Dep::DarkMatter_ID, *Dep::DarkMatter_ID).isSelfConj) ca *= 0.5;
203  result = pow(*Dep::capture_rate_Sun * ca, -0.5);
204 
205  // std::cout << "v = " << sqrt(2.0*T_Sun_core/(*Dep::mwimp)) << " and sigmav inside equilibration_time_Sun = " << sigmav << std::endl;
206  // std::cout << "capture_rate_Sun inside equilibration_time_Sun = " << *Dep::capture_rate_Sun << std::endl;
207  }
208 
210  void annihilation_rate_Sun(double &result)
211  {
212  using namespace Pipes::annihilation_rate_Sun;
213  double tt_sun = 1.5e17 / *Dep::equilibration_time_Sun;
214  result = *Dep::capture_rate_Sun * 0.5 * pow(tanh(tt_sun),2.0);
215  }
216 
219  {
220 
221  using namespace Pipes::nuyield_from_DS;
222  double annihilation_bf[29];
223  double Higgs_decay_BFs_neutral[29][3];
224  double Higgs_decay_BFs_charged[15];
225  double Higgs_masses_neutral[3];
226  double Higgs_mass_charged;
227 
228  // Set annihilation branching fractions
229  std::string DMid = *Dep::DarkMatter_ID;
230  TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMid);
231  std::vector< std::vector<str> > neutral_channels = BEreq::get_DS_neutral_h_decay_channels();
232  // the missing channel
233  const std::vector<str> adhoc_chan = initVector<str>("W-", "H+");
234 
235  for (int i=0; i<29; i++)
236  {
237  const TH_Channel* channel = annProc.find(neutral_channels[i]);
238 
239  if (channel == NULL or i == 26) // Channel 26 has not been implemented in DarkSUSY.
240  {
241  annihilation_bf[i] = 0.;
242  }
243  else
244  {
245  annihilation_bf[i] = channel->genRate->bind("v")->eval(0.);
246  if (i == 10) // Add W- H+ for this channel
247  {
248  channel = annProc.find(adhoc_chan);
249  if (channel == NULL) DarkBit_error().raise(LOCAL_INFO,
250  "W+H- exists in process catalog but not W-H+."
251  " That's some suspiciously severe CP violation yo.");
252  annihilation_bf[i] += channel->genRate->bind("v")->eval(0.);
253  }
254  annihilation_bf[i] /= *Dep::sigmav;
255 
256  // Check that having this channel turned on makes sense at all.
257  #ifdef DARKBIT_DEBUG
258  double mtot = 0;
259  cout << "Particles and masses in DM annihilation final state: " << endl;
260  for (auto p = neutral_channels[i].begin(); p != neutral_channels[i].end(); ++p)
261  {
262  double m = Dep::TH_ProcessCatalog->getParticleProperty(*p).mass;
263  cout << " " << *p << " " << m << endl;
264  mtot += m;
265  }
266  cout << "Sqrt(s) vs total mass of final states: " << 2 * *Dep::mwimp << " vs. " << mtot << endl;
267  cout << "Branching fraction in v=0 limit: " << annihilation_bf[i] << endl << endl;
268  if (mtot > 2 * *Dep::mwimp and annihilation_bf[i] > 0.0)
269  DarkBit_error().raise(LOCAL_INFO, "Channel is open in process catalog but should not be kinematically allowed.");
270  #endif
271  }
272 
273  }
274 
275  // Set Higgs masses
276  if (Dep::TH_ProcessCatalog->hasParticleProperty("h0_1"))
277  Higgs_masses_neutral[1] = Dep::TH_ProcessCatalog->getParticleProperty("h0_1").mass;
278  else
279  DarkBit_error().raise(LOCAL_INFO, "No SM-like Higgs in ProcessCatalog!");
280  Higgs_masses_neutral[0] =
281  Dep::TH_ProcessCatalog->hasParticleProperty("h0_2") ?
282  Dep::TH_ProcessCatalog->getParticleProperty("h0_2").mass : 0.;
283  Higgs_masses_neutral[2] =
284  Dep::TH_ProcessCatalog->hasParticleProperty("A0") ?
285  Dep::TH_ProcessCatalog->getParticleProperty("A0").mass : 0.;
286  Higgs_mass_charged =
287  Dep::TH_ProcessCatalog->hasParticleProperty("H+") ?
288  Dep::TH_ProcessCatalog->getParticleProperty("H+").mass : 0.;
289 
290  // Find out which Higgs exist and have decay data in the process
291  // catalog.
292  const TH_Process* h0_decays[3];
293  h0_decays[0] = Dep::TH_ProcessCatalog->find("h0_2");
294  h0_decays[1] = Dep::TH_ProcessCatalog->find("h0_1");
295  h0_decays[2] = Dep::TH_ProcessCatalog->find("A0");
296  const TH_Process* Hplus_decays = Dep::TH_ProcessCatalog->find("H+");
297  const TH_Process* Hminus_decays = Dep::TH_ProcessCatalog->find("H-");
298  if (Hplus_decays != NULL and Hminus_decays == NULL) DarkBit_error().raise(
299  LOCAL_INFO, "H+ decays exist in process catalog but not H-.");
300  if (Hplus_decays == NULL and Hminus_decays != NULL) DarkBit_error().raise(
301  LOCAL_INFO, "H- decays exist in process catalog but not H+.");
302 
303  // Set the neutral Higgs decay branching fractions
304  for (int i=0; i<3; i++) // Loop over the three neutral Higgs
305  {
306 
307  // If this Higgs exists, set its decay properties.
308  if (h0_decays[i] != NULL)
309  {
310 
311  // Get the total decay width, for normalising partial widths to BFs.
312  // TODO: Replace when BFs become directly available.
313  double totalwidth = 0.0;
314  for (std::vector<TH_Channel>::const_iterator
315  it = h0_decays[i]->channelList.begin();
316  it != h0_decays[i]->channelList.end(); ++it)
317  {
318  // decay width in GeV for two-body final state
319  if ( it->nFinalStates == 2 ) totalwidth +=
320  it->genRate->bind()->eval();
321  }
322 
323  std::vector<str> neutral_channel;
324  // Loop over the decay channels for neutral scalars
325  for (int j=0; j<29; j++)
326  {
327  neutral_channel.clear();
328  neutral_channel.push_back(DarkBit_utils::str_flav_to_mass((neutral_channels[j])[0]));
329  neutral_channel.push_back(DarkBit_utils::str_flav_to_mass((neutral_channels[j])[1]));
330  const TH_Channel* channel = h0_decays[i]->find(neutral_channel);
331  // If this Higgs can decay into this channel, set the BF.
332  if (channel != NULL)
333  {
334  Higgs_decay_BFs_neutral[j][i] = channel->genRate->bind()->eval();
335  if (j == 10) // Add W- H+ for this channel
336  {
337  channel = h0_decays[i]->find(adhoc_chan);
338  if (channel == NULL) DarkBit_error().raise(LOCAL_INFO,
339  "W+H- exists in process catalog but not W-H+."
340  " That's some suspiciously severe CP violation yo.");
341  Higgs_decay_BFs_neutral[j][i]
342  += channel->genRate->bind()->eval();
343  }
344  // This channel has not been implemented in DarkSUSY.
345  if (j == 26) Higgs_decay_BFs_neutral[j][i] = 0.;
346  Higgs_decay_BFs_neutral[j][i] /= totalwidth;
347  }
348  else
349  {
350  Higgs_decay_BFs_neutral[j][i] = 0.;
351  }
352  }
353 
354  }
355 
356  else
357  {
358  // Loop over the decay channels for neutral scalars, setting all BFs
359  // for this Higgs to zero.
360  for (int j=0; j<29; j++) Higgs_decay_BFs_neutral[j][i] = 0.;
361  }
362 
363  }
364 
365  // If they exist, set the charged Higgs decay branching fractions
366  // (DarkSUSY assumes that H+/H- decays are CP-invariant)
367  if (Hplus_decays != NULL)
368  {
369 
370  // Define the charged Higgs decay channels
371  std::vector< std::vector<str> > charged_channels = BEreq::get_DS_charged_h_decay_channels();
372 
373  // Get the total decay width, for normalising partial widths to BFs.
374  // TODO: Replace when BFs become directly available.
375  double totalwidth = 0.0;
376  for (std::vector<TH_Channel>::const_iterator
377  it = Hplus_decays->channelList.begin();
378  it != Hplus_decays->channelList.end(); ++it)
379  {
380  // decay width in GeV for two-body final state
381  if (it->nFinalStates == 2) totalwidth += it->genRate->bind()->eval();
382  }
383 
384  std::vector<str> charged_channel;
385  // Loop over the decay channels for charged scalars
386  for (int j=0; j<15; j++)
387  {
388  charged_channel.clear();
389  charged_channel.push_back(DarkBit_utils::str_flav_to_mass((charged_channels[j])[0]));
390  charged_channel.push_back(DarkBit_utils::str_flav_to_mass((charged_channels[j])[1]));
391  const TH_Channel* channel = Hplus_decays->find(charged_channel);
392  // If this Higgs can decay into this channel, set the BF.
393  if (channel != NULL)
394  {
395  Higgs_decay_BFs_charged[j] = channel->genRate->bind()->eval();
396  Higgs_decay_BFs_charged[j] /= totalwidth;
397  }
398  else
399  {
400  Higgs_decay_BFs_charged[j] = 0.;
401  }
402  }
403 
404  }
405 
406  else
407  {
408  // Loop over the decay channels for charged scalars, setting all BFs
409  // for this Higgs to zero.
410  for (int j=0; j<15; j++) Higgs_decay_BFs_charged[j] = 0.;
411  }
412 
413  // Debug output
414  #ifdef DARKBIT_DEBUG
415  for (int j=0; j<29; j++)
416  {
417  cout<< "annihilation bfs: " << j << " " << annihilation_bf[j] << endl;
418  }
419  cout<< endl;
420  for (int i=0; i<3; i++)
421  {
422  for (int j=0; j<29; j++)
423  {
424  cout<< "higgs neutral bfs: " << i << " " << j << " "
425  << Higgs_decay_BFs_neutral[j][i] << endl;
426  }
427  }
428  cout<< endl;
429  for (int j=0; j<15; j++)
430  {
431  cout<< "higgs charged bfs: " << j << " "
432  << Higgs_decay_BFs_charged[j] << endl;
433  }
434  cout<< endl;
435  for (int j=0; j<3; j++)
436  {
437  cout<< "higgs masses neutral: " << j << " "
438  << Higgs_masses_neutral[j] << endl;
439  }
440  cout<< endl;
441  cout<< "higgs charged mass: " << Higgs_mass_charged << endl;
442  cout<< endl;
443  cout<< "*Dep::mwimp: " << *Dep::mwimp << endl;
444  cout<< endl;
445  cout<< "*Dep::sigmav: " << *Dep::sigmav << endl;
446  cout<< endl;
447  #endif
448 
449  // Set up DarkSUSY to do neutrino yields for this particular WIMP
450  BEreq::DS_nuyield_setup(annihilation_bf, Higgs_decay_BFs_neutral,
451  Higgs_decay_BFs_charged, Higgs_masses_neutral,
452  Higgs_mass_charged, *Dep::mwimp);
453 
454  // Hand back the pointer to the DarkSUSY neutrino yield function
455  result.pointer = BEreq::nuyield.pointer();
456 
457  //TODO: change below to >= when version numbers are available as ints
458  // Treat the yield function as threadsafe only if the loaded version of DarkSUSY supports it.
459  result.threadsafe = (BEreq::nuyield.version() == "5.1.3");
460 
461  // Avoid OpenMP with gfortran 6.x and later, as the implementation seems to be buggy (causes stack overflows).
462  #ifdef __GNUC__
463  #define GCC_VERSION (__GNUC__ * 10000 \
464  + __GNUC_MINOR__ * 100 \
465  + __GNUC_PATCHLEVEL__)
466  #if GCC_VERSION > 60000
467  result.threadsafe = false;
468  #endif
469  #undef GCC_VERSION
470  #endif
471 
472  }
473 
482 
485  void IC22_full(nudata &result)
486  {
487  using namespace Pipes::IC22_full;
488  static bool first = true;
489  const double bgloglike = -808.4581;
490  double sigpred, bgpred, lnLike, pval;
491  int totobs;
492  char experiment[300] = "IC-22";
493  void* context = NULL;
494  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
496  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
497  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
498  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
499  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
503  result.signal = sigpred;
504  result.bg = bgpred;
505  result.nobs = totobs;
506  result.loglike = lnLike;
507  result.pvalue = pval;
508  if (first)
509  {
510  result.bgloglike = bgloglike;
511  first = false;
512  }
513  }
514 
517  void IC79WH_full(nudata &result)
518  {
519  static bool first = true;
520  using namespace Pipes::IC79WH_full;
521  const double bgloglike = -11874.8689;
522  double sigpred, bgpred, lnLike, pval;
523  int totobs;
524  char experiment[300] = "IC-79 WH";
525  void* context = NULL;
526  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
528  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
529  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
530  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
531  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
535  result.signal = sigpred;
536  result.bg = bgpred;
537  result.nobs = totobs;
538  result.loglike = lnLike;
539  result.pvalue = pval;
540  if (first)
541  {
542  result.bgloglike = bgloglike;
543  first = false;
544  }
545  }
546 
549  void IC79WL_full(nudata &result)
550  {
551  static bool first = true;
552  using namespace Pipes::IC79WL_full;
553  const double bgloglike = -1813.4503;
554  double sigpred, bgpred, lnLike, pval;
555  int totobs;
556  char experiment[300] = "IC-79 WL";
557  void* context = NULL;
558  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
560  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
561  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
562  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
563  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
567  result.signal = sigpred;
568  result.bg = bgpred;
569  result.nobs = totobs;
570  result.loglike = lnLike;
571  result.pvalue = pval;
572  if (first)
573  {
574  result.bgloglike = bgloglike;
575  first = false;
576  }
577  }
578 
581  void IC79SL_full(nudata &result)
582  {
583  static bool first = true;
584  using namespace Pipes::IC79SL_full;
585  const double bgloglike = -5015.6474;
586  double sigpred, bgpred, lnLike, pval;
587  int totobs;
588  char experiment[300] = "IC-79 SL";
589  void* context = NULL;
590  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
592  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
593  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
594  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
595  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
599  result.signal = sigpred;
600  result.bg = bgpred;
601  result.nobs = totobs;
602  result.loglike = lnLike;
603  result.pvalue = pval;
604  if (first)
605  {
606  result.bgloglike = bgloglike;
607  first = false;
608  }
609  }
611 
614  void IC22_signal (double &result) {
615  result = Pipes::IC22_signal ::Dep::IC22_data->signal; }
616  void IC22_bg (double &result) {
617  result = Pipes::IC22_bg ::Dep::IC22_data->bg; }
618  void IC22_nobs (int &result) {
619  result = Pipes::IC22_nobs ::Dep::IC22_data->nobs; }
620  void IC22_loglike(double &result) {
621  result = Pipes::IC22_loglike::Dep::IC22_data->loglike; }
622  void IC22_bgloglike(double &result) {
623  result = Pipes::IC22_bgloglike::Dep::IC22_data->bgloglike; }
624  void IC22_pvalue (double &result) {
625  result = Pipes::IC22_pvalue ::Dep::IC22_data->pvalue; }
627 
630  void IC79WH_signal (double &result) {
631  result = Pipes::IC79WH_signal ::Dep::IC79WH_data->signal; }
632  void IC79WH_bg (double &result) {
633  result = Pipes::IC79WH_bg ::Dep::IC79WH_data->bg; }
634  void IC79WH_nobs (int &result) {
635  result = Pipes::IC79WH_nobs ::Dep::IC79WH_data->nobs; }
636  void IC79WH_loglike(double &result) {
637  result = Pipes::IC79WH_loglike::Dep::IC79WH_data->loglike; }
638  void IC79WH_bgloglike(double &result) {
639  result = Pipes::IC79WH_bgloglike::Dep::IC79WH_data->bgloglike; }
640  void IC79WH_pvalue (double &result) {
641  result = Pipes::IC79WH_pvalue ::Dep::IC79WH_data->pvalue; }
643 
646  void IC79WL_signal (double &result) {
647  result = Pipes::IC79WL_signal ::Dep::IC79WL_data->signal; }
648  void IC79WL_bg (double &result) {
649  result = Pipes::IC79WL_bg ::Dep::IC79WL_data->bg; }
650  void IC79WL_nobs (int &result) {
651  result = Pipes::IC79WL_nobs ::Dep::IC79WL_data->nobs; }
652  void IC79WL_loglike(double &result) {
653  result = Pipes::IC79WL_loglike::Dep::IC79WL_data->loglike; }
654  void IC79WL_bgloglike(double &result) {
655  result = Pipes::IC79WL_bgloglike::Dep::IC79WL_data->bgloglike; }
656  void IC79WL_pvalue (double &result) {
657  result = Pipes::IC79WL_pvalue ::Dep::IC79WL_data->pvalue; }
659 
662  void IC79SL_signal (double &result) {
663  result = Pipes::IC79SL_signal ::Dep::IC79SL_data->signal; }
664  void IC79SL_bg (double &result) {
665  result = Pipes::IC79SL_bg ::Dep::IC79SL_data->bg; }
666  void IC79SL_nobs (int &result) {
667  result = Pipes::IC79SL_nobs ::Dep::IC79SL_data->nobs; }
668  void IC79SL_loglike(double &result) {
669  result = Pipes::IC79SL_loglike::Dep::IC79SL_data->loglike; }
670  void IC79SL_bgloglike(double &result) {
671  result = Pipes::IC79SL_bgloglike::Dep::IC79SL_data->bgloglike; }
672  void IC79SL_pvalue (double &result) {
673  result = Pipes::IC79SL_pvalue ::Dep::IC79SL_data->pvalue; }
675 
677  void IC79_loglike(double &result)
678  {
679  using namespace Pipes::IC79_loglike;
683  if (result > 0.0) result = 0.0;
684  }
685 
687  void IC_loglike(double &result)
688  {
689  using namespace Pipes::IC_loglike;
694  if (result > 0.0) result = 0.0;
695 #ifdef DARKBIT_DEBUG
696  cout << "IC likelihood: " << result << endl;
697  cout << "IC79SL contribution: " << *Dep::IC79SL_loglike - *Dep::IC79SL_bgloglike << endl;
698  cout << "IC79WL contribution: " << *Dep::IC79WL_loglike - *Dep::IC79WL_bgloglike << endl;
699  cout << "IC79WH contribution: " << *Dep::IC79WH_loglike - *Dep::IC79WH_bgloglike << endl;
700  cout << "IC22 contribution: " << *Dep::IC22_loglike - *Dep::IC22_bgloglike << endl;
701 #endif
702  }
703 
706  {
708 
709  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
710 
711  double rho0 = LocalHaloParameters.rho0;
712  double rho0_eff = (*Dep::RD_fraction)*rho0;
713  double vrot = LocalHaloParameters.vrot;
714  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
715  double vesc = LocalHaloParameters.vesc;
717  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
718 
719  BEreq::dshmcom->rho0 = rho0;
720  BEreq::dshmcom->v_sun = vrot;
721  BEreq::dshmcom->v_earth = v_earth;
722  BEreq::dshmcom->rhox = rho0_eff;
723 
724  BEreq::dshmframevelcom->v_obs = vrot;
725 
726  BEreq::dshmisodf->vd_3d = vd_3d;
727  BEreq::dshmisodf->vgalesc = vesc;
728 
729  BEreq::dshmnoclue->vobs = vrot;
730 
732  << "Updating DarkSUSY halo parameters:" << std::endl
733  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
734  << " rho0_eff [GeV/cm^3] = " << rho0_eff << std::endl
735  << " v_sun [km/s] = " << vrot<< std::endl
736  << " v_earth [km/s] = " << v_earth << std::endl
737  << " v_obs [km/s] = " << vrot << std::endl
738  << " vd_3d [km/s] = " << vd_3d << std::endl
739  << " v_esc [km/s] = " << vesc << EOM;
740 
741  result = true;
742 
743  return;
744  }
745 
748  {
750 
751  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
752 
753  double rho0 = LocalHaloParameters.rho0;
754  double vrot = LocalHaloParameters.vrot;
755  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
756  double vesc = LocalHaloParameters.vesc;
758  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
759 
760  BEreq::dshmcom->rho0 = rho0;
761  BEreq::dshmcom->v_sun = vrot;
762  BEreq::dshmcom->v_earth = v_earth;
763 
764  BEreq::dshmframevelcom->v_obs = vrot;
765 
766  BEreq::dshmisodf->vd_3d = vd_3d;
767  BEreq::dshmisodf->vgalesc = vesc;
768 
769  BEreq::dshmnoclue->vobs = vrot;
770 
772  << "Updating DarkSUSY halo parameters:" << std::endl
773  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
774  << " v_sun [km/s] = " << vrot<< std::endl
775  << " v_earth [km/s] = " << v_earth << std::endl
776  << " v_obs [km/s] = " << vrot << std::endl
777  << " vd_3d [km/s] = " << vd_3d << std::endl
778  << " v_esc [km/s] = " << vesc << EOM;
779 
780  result = true;
781 
782  return;
783  }
784 
785  }
786 }
787 
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:524
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:99
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
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:544
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:473
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.