gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
ScalarSingletDM.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
23 
28 #include "boost/make_shared.hpp"
30 
31 namespace Gambit
32 {
33 
34  namespace DarkBit
35  {
36  //#define DARKBIT_DEBUG
37 
38  class ScalarSingletDM
39  {
40  public:
42  ScalarSingletDM(
43  TH_ProcessCatalog* const catalog,
44  double gammaH,
45  double vev,
46  double alpha_strong,
47  double vSigma_s)
48  : Gamma_mh(gammaH), v0 (vev),
49  alpha_s (alpha_strong), vSigma_semi (vSigma_s)
50  {
51  mh = catalog->getParticleProperty("h0_1").mass;
52  mb = catalog->getParticleProperty("d_3").mass;
53  mc = catalog->getParticleProperty("u_2").mass;
54  mtau = catalog->getParticleProperty("e-_3").mass;
55  mt = catalog->getParticleProperty("u_3").mass;
56  mZ0 = catalog->getParticleProperty("Z0").mass;
57  mW = catalog->getParticleProperty("W+").mass;
58  mS = catalog->getParticleProperty("S").mass;
59  };
60  ~ScalarSingletDM() {}
61 
63  double Dh2 (double s)
64  {
65  return 1/((s-mh*mh)*(s-mh*mh)+mh*mh*Gamma_mh*Gamma_mh);
66  }
67 
74  double sv(std::string channel, double lambda, double mass, double v)
75  {
76 
77  // Note: Valid for mass > 45 GeV
78  double s = 4*mass*mass/(1-v*v/4);
79  double sqrt_s = sqrt(s);
80  if ( sqrt_s < 90 )
81  {
83  "SingletDM sigmav called with sqrt_s < 90 GeV.");
84  return 0;
85  }
86 
87  if ( channel == "Sh" )
88  {
89  if (sqrt_s > (mh+mS)) return vSigma_semi;
90  else return 0;
91  }
92 
93  if ( channel == "hh" )
94  {
95  if ( sqrt_s > mh*2 )
96  {
97  double GeV2tocm3s1 = gev2cm2*s2cm;
98  return sv_hh(lambda, mass, v)*GeV2tocm3s1;
99  }
100  else return 0;
101  }
102 
103  if ( channel == "bb" and sqrt_s < mb*2 ) return 0;
104  if ( channel == "cc" and sqrt_s < mc*2 ) return 0;
105  if ( channel == "tautau" and sqrt_s < mtau*2 ) return 0;
106  if ( channel == "tt" and sqrt_s < mt*2 ) return 0;
107  if ( channel == "ZZ" and sqrt_s < mZ0*2) return 0;
108  if ( channel == "WW" and sqrt_s < mW*2) return 0;
109 
110  if ( sqrt_s < 300 )
111  {
112  double br = virtual_SMHiggs_widths(channel,sqrt_s);
113  double Gamma_s = virtual_SMHiggs_widths("Gamma",sqrt_s);
114  double GeV2tocm3s1 = gev2cm2*s2cm;
115 
116  // Explicitly close channel for off-shell top quarks
117  if ( channel == "tt" and sqrt_s < mt*2) return 0;
118 
119  double res = 2*lambda*lambda*v0*v0/
120  sqrt_s*Dh2(s)*Gamma_s*GeV2tocm3s1*br;
121  return res;
122  }
123  else
124  {
125  if ( channel == "bb" ) return sv_ff(lambda, mass, v, mb, true);
126  if ( channel == "cc" ) return sv_ff(lambda, mass, v, mc, true);
127  if ( channel == "tautau" ) return sv_ff(lambda, mass, v, mtau, false);
128  if ( channel == "tt" ) return sv_ff(lambda, mass, v, mt, true);
129  if ( channel == "ZZ" ) return sv_ZZ(lambda, mass, v);
130  if ( channel == "WW" ) return sv_WW(lambda, mass, v);
131  }
132  return 0;
133  }
134 
135  // Annihilation into W bosons.
136  double sv_WW(double lambda, double mass, double v)
137  {
138  double s = 4*mass*mass/(1-v*v/4);
139  double x = pow(mW,2)/s;
140  double GeV2tocm3s1 = gev2cm2*s2cm;
141  return pow(lambda,2)*s/8/M_PI*sqrt(1-4*x)*Dh2(s)*(1-4*x+12*pow(x,2))
142  *GeV2tocm3s1;
143  }
144 
145  // Annihilation into Z bosons.
146  double sv_ZZ(double lambda, double mass, double v)
147  {
148  double s = 4*mass*mass/(1-v*v/4);
149  double x = pow(mZ0,2)/s;
150  double GeV2tocm3s1 = gev2cm2*s2cm;
151  return pow(lambda,2)*s/16/M_PI*sqrt(1-4*x)*Dh2(s)*(1-4*x+12*pow(x,2))
152  * GeV2tocm3s1;
153  }
154 
155  // Annihilation into fermions
156  double sv_ff(
157  double lambda, double mass, double v, double mf, bool is_quark)
158  {
159  double s = 4*mass*mass/(1-v*v/4);
160  double vf = sqrt(1-4*pow(mf,2)/s);
161  double Xf = 1;
162  if ( is_quark ) Xf = 3 *
163  (1+(3/2*log(pow(mf,2)/s)+9/4)*4*alpha_s/3/M_PI);
164  double GeV2tocm3s1 = gev2cm2*s2cm;
165  return pow(lambda,2)*
166  pow(mf,2)/4/M_PI*Xf*pow(vf,3) * Dh2(s) * GeV2tocm3s1;
167  }
168 
170  double sv_hh(double lambda, double mass, double v)
171  {
172 
173  double s = 4*mass*mass/(1-v*v/4); // v is relative velocity
174  double vh = sqrt(1-4*mh*mh/s); // vh and vs are lab velocities
175  // Hardcoded lower velocity avoids nan results
176  double vs = std::max(v/2, 1e-6);
177  double tp = pow(mass,2)+pow(mh,2)-0.5*s*(1-vs*vh);
178  double tm = pow(mass,2)+pow(mh,2)-0.5*s*(1+vs*vh);
179 
180  double aR = 1+3*mh*mh*(s-mh*mh)*Dh2(s);
181  double aI = 3*mh*mh*sqrt(s)*Gamma_mh*Dh2(s);
182 
183  return pow(lambda,2)/16/M_PI/pow(s,2)/vs *
184  (
185  (pow(aR,2)+pow(aI,2))*s*vh*vs
186  +4*lambda*pow(v0,2)*(aR-lambda*pow(v0,2)/(s-2*pow(mh,2)))
187  *log(std::abs(pow(mass,2)-tp)/std::abs(pow(mass,2)-tm))
188  +(2*pow(lambda,2)*pow(v0,4)*s*vh*vs)
189  /(pow(mass,2)-tm)/(pow(mass,2)-tp));
190  }
191 
192  private:
193  double Gamma_mh, mh, v0, alpha_s, mb, mc, mtau, mt, mZ0, mW, mS, vSigma_semi;
194  };
195 
196  void DarkMatter_ID_ScalarSingletDM(std::string & result) { result = "S"; }
197 
199  void get_ScalarSingletDM_DD_couplings(const Spectrum &spec, DM_nucleon_couplings &result, Models::safe_param_map<safe_ptr<const double> > &Param)
200  {
201  const SubSpectrum& he = spec.get_HE();
202  double mass = spec.get(Par::Pole_Mass,"S");
203  double lambda = he.get(Par::dimensionless,"lambda_hS");
204  double mh = spec.get(Par::Pole_Mass,"h0_1");
205 
206  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
207  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
208  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
209 
210  result.gps = lambda*fp*m_proton/pow(mh,2)/mass/2;
211  result.gns = lambda*fn*m_neutron/pow(mh,2)/mass/2;
212  result.gpa = 0; // Only SI cross-section
213  result.gna = 0;
214 
215  logger() << LogTags::debug << "Singlet DM DD couplings:" << std::endl;
216  logger() << " gps = " << result.gps << std::endl;
217  logger() << " gns = " << result.gns << std::endl;
218  logger() << " gpa = " << result.gpa << std::endl;
219  logger() << " gna = " << result.gna << EOM;
220  }
221 
223  void DD_couplings_ScalarSingletDM_Z2(DM_nucleon_couplings &result)
224  {
227  get_ScalarSingletDM_DD_couplings(spec, result, Param);
228  }
229 
231  void DD_couplings_ScalarSingletDM_Z3(DM_nucleon_couplings &result)
232  {
235  get_ScalarSingletDM_DD_couplings(spec, result, Param);
236  }
237 
238 
241  {
243  using std::vector;
244  using std::string;
245 
246  // Initialize empty catalog and main annihilation process
247  TH_ProcessCatalog catalog;
248  TH_Process process_ann("S", "S");
249 
250  // Explicitly state that Z2 Scalar DM is self-conjugate
251  process_ann.isSelfConj = true;
252 
254  // Import particle masses and couplings
256 
257  // Convenience macros
258  #define getSMmass(Name, spinX2) \
259  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
260  (Name , TH_ParticleProperty(SM.get(Par::Pole_Mass,Name), spinX2)));
261  #define addParticle(Name, Mass, spinX2) \
262  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
263  (Name , TH_ParticleProperty(Mass, spinX2)));
264 
265  // Import Spectrum objects
267  const SubSpectrum& he = spec.get_HE();
268  const SubSpectrum& SM = spec.get_LE();
269  const SMInputs& SMI = spec.get_SMInputs();
270 
271  // Import couplings
272  double lambda = he.get(Par::dimensionless,"lambda_hS");
273  double v = he.get(Par::mass1,"vev");
274 
275  // Get SM pole masses
276  getSMmass("e-_1", 1)
277  getSMmass("e+_1", 1)
278  getSMmass("e-_2", 1)
279  getSMmass("e+_2", 1)
280  getSMmass("e-_3", 1)
281  getSMmass("e+_3", 1)
282  getSMmass("Z0", 2)
283  getSMmass("W+", 2)
284  getSMmass("W-", 2)
285  getSMmass("g", 2)
286  getSMmass("gamma", 2)
287  getSMmass("u_3", 1)
288  getSMmass("ubar_3", 1)
289  getSMmass("d_3", 1)
290  getSMmass("dbar_3", 1)
291 
292  // Pole masses not available for the light quarks.
293  addParticle("u_1" , SMI.mU, 1) // mu(2 GeV)^MS-bar, not pole mass
294  addParticle("ubar_1", SMI.mU, 1) // mu(2 GeV)^MS-bar, not pole mass
295  addParticle("d_1" , SMI.mD, 1) // md(2 GeV)^MS-bar, not pole mass
296  addParticle("dbar_1", SMI.mD, 1) // md(2 GeV)^MS-bar, not pole mass
297  addParticle("u_2" , SMI.mCmC,1) // mc(mc)^MS-bar, not pole mass
298  addParticle("ubar_2", SMI.mCmC,1) // mc(mc)^MS-bar, not pole mass
299  addParticle("d_2" , SMI.mS, 1) // ms(2 GeV)^MS-bar, not pole mass
300  addParticle("dbar_2", SMI.mS, 1) // ms(2 GeV)^MS-bar, not pole mass
301  double alpha_s = SMI.alphaS; // alpha_s(mZ)^MSbar
302 
303  // Masses for neutrino flavour eigenstates. Set to zero.
304  // (presently not required)
305  addParticle("nu_e", 0.0, 1)
306  addParticle("nubar_e", 0.0, 1)
307  addParticle("nu_mu", 0.0, 1)
308  addParticle("nubar_mu", 0.0, 1)
309  addParticle("nu_tau", 0.0, 1)
310  addParticle("nubar_tau",0.0, 1)
311 
312  // Higgs-sector masses
313  double mS = spec.get(Par::Pole_Mass,"S");
314  double mH = spec.get(Par::Pole_Mass,"h0_1");
315  addParticle("S", mS, 0) // Scalar Singlet DM
316  addParticle("h0_1", mH, 0) // SM-like Higgs
317  addParticle("pi0", meson_masses.pi0, 0)
318  addParticle("pi+", meson_masses.pi_plus, 0)
319  addParticle("pi-", meson_masses.pi_minus, 0)
320  addParticle("eta", meson_masses.eta, 0)
321  addParticle("rho0", meson_masses.rho0, 1)
322  addParticle("rho+", meson_masses.rho_plus, 1)
323  addParticle("rho-", meson_masses.rho_minus, 1)
324  addParticle("omega", meson_masses.omega, 1)
325 
326  // Get rid of convenience macros
327  #undef getSMmass
328  #undef addParticle
329 
330 
332  // Import Decay information
334 
335  // Import decay table from DecayBit
336  const DecayTable* tbl = &(*Dep::decay_rates);
337 
338  // Save Higgs width for later
339  double gammaH = tbl->at("h0_1").width_in_GeV;
340 
341  // Set of imported decays
342  std::set<string> importedDecays;
343 
344  // Minimum branching ratio to include
345  double minBranching = 0;
346 
347  // Import relevant decays (only Higgs and subsequent decays)
349  // Notes: Virtual Higgs decays into offshell W+W- final states are not
350  // imported. All other channels are correspondingly rescaled. Decay
351  // into SS final states is accounted for, leading to zero photons.
352  ImportDecays("h0_1", catalog, importedDecays, tbl, minBranching,
353  daFunk::vec<std::string>("Z0", "W+", "W-", "e+_2", "e-_2", "e+_3", "e-_3"));
354 
355  // Instantiate new ScalarSingletDM object
356  auto singletDM = boost::make_shared<ScalarSingletDM>(&catalog, gammaH, v, alpha_s, 0.0);
357 
358  // Populate annihilation channel list and add thresholds to threshold
359  // list.
360  // (remark: the lowest threshold is here = 2*mS, whereas in DS-internal
361  // conventions, this lowest threshold is not listed)
362  process_ann.resonances_thresholds.threshold_energy.push_back(2*mS);
363  auto channel =
364  daFunk::vec<string>("bb", "WW", "cc", "tautau", "ZZ", "tt", "hh");
365  auto p1 =
366  daFunk::vec<string>("d_3", "W+", "u_2", "e+_3", "Z0", "u_3", "h0_1");
367  auto p2 =
368  daFunk::vec<string>("dbar_3","W-", "ubar_2","e-_3", "Z0", "ubar_3","h0_1");
369  {
370  for ( unsigned int i = 0; i < channel.size(); i++ )
371  {
372  double mtot_final =
373  catalog.getParticleProperty(p1[i]).mass +
374  catalog.getParticleProperty(p2[i]).mass;
375  // Include final states that are open for T~m/20
376  if ( mS*2 > mtot_final*0.5 )
377  {
378  daFunk::Funk kinematicFunction = daFunk::funcM(singletDM,
379  &ScalarSingletDM::sv, channel[i], lambda, mS, daFunk::var("v"));
380  TH_Channel new_channel(
381  daFunk::vec<string>(p1[i], p2[i]), kinematicFunction
382  );
383  process_ann.channelList.push_back(new_channel);
384  }
385  if ( mS*2 > mtot_final )
386  {
388  push_back(mtot_final);
389  }
390  }
391  }
392 
393  // Populate resonance list
394  if ( mH >= mS*2 ) process_ann.resonances_thresholds.resonances.
395  push_back(TH_Resonance(mH, gammaH));
396 
397  catalog.processList.push_back(process_ann);
398 
399  // Validate
400  catalog.validate();
401 
402 
403  #ifdef DARKBIT_DEBUG
404  // get DarkBit computed vSigma total
405  // so we can compare with that from MO
406  ScalarSingletDM test = *singletDM;
407  int nc = 7;
408  double total = 0;
409  for (int ii=0; ii < nc ; ii++)
410  {
411  total = total + test.sv(channel[ii], lambda, mS, 0.0);
412  }
413  cout << " --- Testing process catalouge --- " << endl;
414  cout << "Total sigma V from process catalouge = " << total << endl;
415  cout << " ---------- " << endl;
416  #endif
417 
418 
419 
420  result = catalog;
421  } // function TH_ProcessCatalog_ScalarSingletDM_Z2
422 
423 
424 
427  {
429  using std::vector;
430  using std::string;
431 
432  // Initialize empty catalog and main annihilation process
433  TH_ProcessCatalog catalog;
434  TH_Process process_ann("S", "S");
435 
436  // Explicitly state that Z3 Scalar DM is not self-conjugate
437  process_ann.isSelfConj = false;
438 
439  // get semi-annihilation cross-section from MicrOmegas
440 
441  // set requested spectra to NULL since we don't need them
442  double * SpNe=NULL,*SpNm=NULL,*SpNl=NULL;
443  double * SpA=NULL,*SpE=NULL,*SpP=NULL;
444  int err;
445 
446  double vSigma_total = BEreq::calcSpectrum(byVal(3),byVal(SpA),byVal(SpE),byVal(SpP),byVal(SpNe),byVal(SpNm),byVal(SpNl) ,byVal(&err));
447 
448  if (err != 0 )
449  {
450  DarkBit_error().raise(LOCAL_INFO, "MicrOmegas spectrum calculation returned error code = " + std::to_string(err));
451  }
452 
453  // get BR for each channel as filled by calcSpectrum
454  MicrOmegas::aChannel* vSigmaCh;
455  vSigmaCh = *BEreq::vSigmaCh;
456 
457  int n_channels = sizeof(vSigmaCh);
458  double BR_semi = 0; // semi-annihilation BR
459  #ifdef DARKBIT_DEBUG
460  cout << "--- Semi-annihilation processes and BRs --- " << endl;
461  #endif
462  // find semi-annihilation channels
463  int i = 0;
464  while ((vSigmaCh[i].weight!=0) && (i < n_channels))
465  {
466  // get final states
467  const char* p1 = vSigmaCh[i].prtcl[2];
468  const char* p2 = vSigmaCh[i].prtcl[3];
469  // semi-annihilation final states are h + ~ss or h + ~SS
470 
471  if ( (strcmp(p1,"h")==0) && ( (strcmp(p2,"~ss")==0) || (strcmp(p2,"~SS")==0) ) )
472  {
473  BR_semi = BR_semi + vSigmaCh[i].weight;
474  #ifdef DARKBIT_DEBUG
475  cout << "process: " << vSigmaCh[i].prtcl[0] << " + ";
476  cout << vSigmaCh[i].prtcl[1] << " -> " << p1 << " + " << p2;
477  cout << ", BR = " << vSigmaCh[i].weight << endl;
478  #endif
479  }
480  i++;
481  }
482 
483  // Get the semi-annihilation cross-section, accounting for the fact that MO already scales
484  // vSigma_total down by a factor of 2 to account for the fact that Z3 scalar singlet DM is not self-conjugate.
485  double vSigma_semi = BR_semi * 2.0 * vSigma_total;
486 
487  #ifdef DARKBIT_DEBUG
488  cout << "--- --- " << endl;
489  cout << "Total sigma v from MicrOmegas = " << vSigma_total << " cm^3/s" << endl;
490  cout << "semi-annihilation BR = " << BR_semi << endl;
491  cout << "semi-annihilation sigma v from MicrOmegas = " << 0.5*vSigma_semi << " cm^3/s" << endl;
492  cout << "Total sigma v from MicrOmegas excluding semi-annihilations = " << vSigma_total - 0.5*vSigma_semi << endl;
493  cout << "--------- " << endl;
494  #endif
495 
496 
498  // Import particle masses and couplings
500 
501  // Convenience macros
502  #define getSMmass(Name, spinX2) \
503  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
504  (Name , TH_ParticleProperty(SM.get(Par::Pole_Mass,Name), spinX2)));
505  #define addParticle(Name, Mass, spinX2) \
506  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
507  (Name , TH_ParticleProperty(Mass, spinX2)));
508 
509  // Import Spectrum objects
511  const SubSpectrum& he = spec.get_HE();
512  const SubSpectrum& SM = spec.get_LE();
513  const SMInputs& SMI = spec.get_SMInputs();
514 
515  // Import couplings
516  double lambda = he.get(Par::dimensionless,"lambda_hS");
517  double v = he.get(Par::mass1,"vev");
518 
519  // Get SM pole masses
520  getSMmass("e-_1", 1)
521  getSMmass("e+_1", 1)
522  getSMmass("e-_2", 1)
523  getSMmass("e+_2", 1)
524  getSMmass("e-_3", 1)
525  getSMmass("e+_3", 1)
526  getSMmass("Z0", 2)
527  getSMmass("W+", 2)
528  getSMmass("W-", 2)
529  getSMmass("g", 2)
530  getSMmass("gamma", 2)
531  getSMmass("u_3", 1)
532  getSMmass("ubar_3", 1)
533  getSMmass("d_3", 1)
534  getSMmass("dbar_3", 1)
535 
536  // Pole masses not available for the light quarks.
537  addParticle("u_1" , SMI.mU, 1) // mu(2 GeV)^MS-bar, not pole mass
538  addParticle("ubar_1", SMI.mU, 1) // mu(2 GeV)^MS-bar, not pole mass
539  addParticle("d_1" , SMI.mD, 1) // md(2 GeV)^MS-bar, not pole mass
540  addParticle("dbar_1", SMI.mD, 1) // md(2 GeV)^MS-bar, not pole mass
541  addParticle("u_2" , SMI.mCmC,1) // mc(mc)^MS-bar, not pole mass
542  addParticle("ubar_2", SMI.mCmC,1) // mc(mc)^MS-bar, not pole mass
543  addParticle("d_2" , SMI.mS, 1) // ms(2 GeV)^MS-bar, not pole mass
544  addParticle("dbar_2", SMI.mS, 1) // ms(2 GeV)^MS-bar, not pole mass
545  double alpha_s = SMI.alphaS; // alpha_s(mZ)^MSbar
546 
547  // Masses for neutrino flavour eigenstates. Set to zero.
548  // (presently not required)
549  addParticle("nu_e", 0.0, 1)
550  addParticle("nubar_e", 0.0, 1)
551  addParticle("nu_mu", 0.0, 1)
552  addParticle("nubar_mu", 0.0, 1)
553  addParticle("nu_tau", 0.0, 1)
554  addParticle("nubar_tau",0.0, 1)
555 
556  // Higgs-sector masses
557  double mS = spec.get(Par::Pole_Mass,"S");
558  double mH = spec.get(Par::Pole_Mass,"h0_1");
559  addParticle("S", mS, 0) // Scalar Singlet DM
560  addParticle("h0_1", mH, 0) // SM-like Higgs
561  addParticle("pi0", meson_masses.pi0, 0)
562  addParticle("pi+", meson_masses.pi_plus, 0)
563  addParticle("pi-", meson_masses.pi_minus, 0)
564  addParticle("eta", meson_masses.eta, 0)
565  addParticle("rho0", meson_masses.rho0, 1)
566  addParticle("rho+", meson_masses.rho_plus, 1)
567  addParticle("rho-", meson_masses.rho_minus, 1)
568  addParticle("omega", meson_masses.omega, 1)
569 
570  // Get rid of convenience macros
571  #undef getSMmass
572  #undef addParticle
573 
574 
576  // Import Decay information
578 
579  // Import decay table from DecayBit
580  const DecayTable* tbl = &(*Dep::decay_rates);
581 
582  // Save Higgs width for later
583  double gammaH = tbl->at("h0_1").width_in_GeV;
584 
585  // Set of imported decays
586  std::set<string> importedDecays;
587 
588  // Minimum branching ratio to include
589  double minBranching = 0;
590 
591  // Import relevant decays (only Higgs and subsequent decays)
593  // Notes: Virtual Higgs decays into offshell W+W- final states are not
594  // imported. All other channels are correspondingly rescaled. Decay
595  // into SS final states is accounted for, leading to zero photons.
596  ImportDecays("h0_1", catalog, importedDecays, tbl, minBranching,
597  daFunk::vec<std::string>("Z0", "W+", "W-", "e+_2", "e-_2", "e+_3", "e-_3"));
598 
599  // Instantiate new ScalarSingletDM object
600  auto singletDM = boost::make_shared<ScalarSingletDM>(&catalog, gammaH, v, alpha_s, vSigma_semi);
601 
602  // Populate annihilation channel list and add thresholds to threshold
603  // list.
604  // (remark: the lowest threshold is here = 2*mS, whereas in DS-internal
605  // conventions, this lowest threshold is not listed)
606  process_ann.resonances_thresholds.threshold_energy.push_back(2*mS);
607  auto channel =
608  daFunk::vec<string>("bb", "WW", "cc", "tautau", "ZZ", "tt", "hh", "Sh");
609  auto p1 =
610  daFunk::vec<string>("d_3", "W+", "u_2", "e+_3", "Z0", "u_3", "h0_1", "S");
611  auto p2 =
612  daFunk::vec<string>("dbar_3","W-", "ubar_2","e-_3", "Z0", "ubar_3", "h0_1", "h0_1");
613  {
614  for ( unsigned int i = 0; i < channel.size(); i++ )
615  {
616  double mtot_final =
617  catalog.getParticleProperty(p1[i]).mass +
618  catalog.getParticleProperty(p2[i]).mass;
619  // Include final states that are open for T~m/20
620  if ( mS*2 > mtot_final*0.5 )
621  {
622  daFunk::Funk kinematicFunction = daFunk::funcM(singletDM,
623  &ScalarSingletDM::sv, channel[i], lambda, mS, daFunk::var("v"));
624  TH_Channel new_channel(
625  daFunk::vec<string>(p1[i], p2[i]), kinematicFunction
626  );
627  process_ann.channelList.push_back(new_channel);
628  }
629  if ( mS*2 > mtot_final )
630  {
632  push_back(mtot_final);
633  }
634  }
635  }
636 
637  // Populate resonance list
638  if ( mH >= mS*2 ) process_ann.resonances_thresholds.resonances.
639  push_back(TH_Resonance(mH, gammaH));
640 
641  catalog.processList.push_back(process_ann);
642 
643  // Validate
644  catalog.validate();
645 
646  #ifdef DARKBIT_DEBUG
647  // get DarkBit computed vSigma total
648  // so we can compare with that from MO
649  ScalarSingletDM test = *singletDM;
650  int nc = 7;
651  double total = 0;
652  for (int ii=0; ii < nc ; ii++)
653  {
654  total = total + test.sv(channel[ii], lambda, mS, 0.0);
655  }
656  cout << " --- Testing process catalouge --- " << endl;
657  cout << "Total sigma V from process catalouge excluding semi-annihilations = " << total << endl;
658  total = total + test.sv("Sh", lambda, mS, 0.0);
659  cout << "Total including semi-annihilations = " << total << endl;
660  cout << " ---------- " << endl;
661  #endif
662 
663  result = catalog;
664  } // function TH_ProcessCatalog_ScalarSingletDM_Z3
665 
666 
667 
668  }
669 }
std::vector< TH_Channel > channelList
List of channels.
error & DarkBit_error()
double get(const Par::Tags partype, const std::string &mass) const
Definition: spectrum.cpp:249
Simple function for returning SM Higgs partial and total widths as a function of virtuality, as read from data files from CERN Yellow Paper arXiv:1101.0593.
std::chrono::steady_clock::time_point tp
void DD_couplings_ScalarSingletDM_Z2(DM_nucleon_couplings &result)
Direct detection couplings for Z2 scalar singlet DM.
const double m_proton
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 ImportDecays(std::string pID, TH_ProcessCatalog &catalog, std::set< std::string > &importedDecays, const DecayTable *tbl, double minBranching, std::vector< std::string > excludeDecays=std::vector< std::string >())
Funk funcM(O *obj, double(O::*f)(funcargs...), Args... args)
Definition: daFunk.hpp:861
const double s2cm
const double m_neutron
Simple reader for ASCII tables.
#define addParticle(Name, Mass, spinX2)
#define LOCAL_INFO
Definition: local_info.hpp:34
const double mt
Definition: topness.h:39
START_MODEL v0
void TH_ProcessCatalog_ScalarSingletDM_Z2(DarkBit::TH_ProcessCatalog &result)
Set up process catalog for Z2 scalar singlet DM.
Funk var(std::string arg)
Definition: daFunk.hpp:921
const double mW
Definition: topness.h:40
double virtual_SMHiggs_widths(str, double)
Higgs branching ratios and total width Gamma [GeV], as function of mass [GeV] (90 - 300 GeV) ...
void DD_couplings_ScalarSingletDM_Z3(DM_nucleon_couplings &result)
Direct detection couplings for Z3 scalar singlet DM.
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
bool isSelfConj
Does the process contain self-conjugate DM? (accounting for correct factors of 1/2 in annihilation sp...
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum ScalarSingletDM_Z3_spectrum
A single resonance of a given width at a given energy (both in GeV)
virtual double get(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
void DarkMatter_ID_ScalarSingletDM(std::string &result)
TH_resonances_thresholds resonances_thresholds
List of resonances and thresholds.
double width_in_GeV
Total particle width (in GeV)
void TH_ProcessCatalog_ScalarSingletDM_Z3(DarkBit::TH_ProcessCatalog &result)
Set up process catalog for Z3 scalar singlet DM.
A container holding all annihilation and decay initial states relevant for DarkBit.
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
SubSpectrum & get_LE()
Standard getters Return references to internal data members.
Definition: spectrum.cpp:224
All data on a single annihilation or decay channel, e.g.
const double gev2cm2
double lambda(double x, double y, double z)
Definition: MSSM_H.hpp:38
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
#define getSMmass(Name, spinX2)
Rollcall header for module DarkBit.
GAMBIT native decay table class.
Definition: decay_table.hpp:35
double pow(const double &a)
Outputs a^i.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum mh
void request(std::string message)
Request an exception.
Definition: exceptions.cpp:463
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 get_ScalarSingletDM_DD_couplings(const Spectrum &spec, DM_nucleon_couplings &result, Models::safe_param_map< safe_ptr< const double > > &Param)
Common code for different scalar singlet direct detection coupling routines.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
SubSpectrum & get_HE()
Definition: spectrum.cpp:225
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110
SMInputs & get_SMInputs()
Definition: spectrum.cpp:226
Utility functions for DarkBit.
Container class for Standard Model input information (defined as in SLHA2)
Definition: sminputs.hpp:29