gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
MajoranaSingletDM.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
26 
31 #include "boost/make_shared.hpp"
33 
34 namespace Gambit
35 {
36 
37  namespace DarkBit
38  {
39 
40  class MajoranaSingletDM
41  {
42  public:
44  MajoranaSingletDM(
45  TH_ProcessCatalog* const catalog,
46  double gammaH,
47  double vev,
48  double alpha_strong)
49  : Gamma_mh(gammaH), v0 (vev),
50  alpha_s (alpha_strong)
51  {
52  mh = catalog->getParticleProperty("h0_1").mass;
53  mb = catalog->getParticleProperty("d_3").mass;
54  mc = catalog->getParticleProperty("u_2").mass;
55  mtau = catalog->getParticleProperty("e-_3").mass;
56  mt = catalog->getParticleProperty("u_3").mass;
57  mZ0 = catalog->getParticleProperty("Z0").mass;
58  mW = catalog->getParticleProperty("W+").mass;
59  };
60  ~MajoranaSingletDM() {}
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 cosXi, double v)
75  {
76  // Note: Valid for mass > 45 GeV
77 
78  // Hardcoded velocity avoids NaN results.
79  v = std::max(v, 1e-6);
80 
81  double s = 4*mass*mass/(1-v*v/4);
82  double sqrt_s = sqrt(s);
83  if ( sqrt_s < 90 )
84  {
86  "MajoranaSingletDM sigmav called with sqrt_s < 90 GeV.");
87  return 0;
88  }
89 
90  if ( channel == "hh" )
91  {
92  if ( sqrt_s > mh*2 )
93  {
94  return sv_hh(lambda, mass, v, cosXi);
95  }
96  else return 0;
97  }
98 
99  if ( channel == "bb" and sqrt_s < mb*2 ) return 0;
100  if ( channel == "cc" and sqrt_s < mc*2 ) return 0;
101  if ( channel == "tautau" and sqrt_s < mtau*2 ) return 0;
102  if ( channel == "tt" and sqrt_s < mt*2 ) return 0;
103  if ( channel == "ZZ" and sqrt_s < mZ0*2) return 0;
104  if ( channel == "WW" and sqrt_s < mW*2) return 0;
105 
106  if ( sqrt_s < 300 )
107  {
108  double br = virtual_SMHiggs_widths(channel,sqrt_s);
109  double Gamma_s = virtual_SMHiggs_widths("Gamma",sqrt_s);
110  double GeV2tocm3s1 = gev2cm2*s2cm;
111  double cos2Xi = cosXi*cosXi;
112  double sin2Xi = 1 - cos2Xi;
113  double numerator = (cos2Xi*v*v/4 + sin2Xi);
114 
115  // Explicitly close channel for off-shell top quarks
116  if ( channel == "tt" and sqrt_s < mt*2) return 0;
117 
118  double res = numerator*lambda*lambda*v0*v0*sqrt_s
119  *Dh2(s)*Gamma_s*GeV2tocm3s1*br;
120  return res;
121  }
122  else
123  {
124  if ( channel == "bb" ) return sv_ff(lambda, mass, v, mb, cosXi, true);
125  if ( channel == "cc" ) return sv_ff(lambda, mass, v, mc, cosXi, true);
126  if ( channel == "tautau" ) return sv_ff(lambda, mass, v, mtau, cosXi, false);
127  if ( channel == "tt" ) return sv_ff(lambda, mass, v, mt, cosXi, true);
128  if ( channel == "ZZ" ) return sv_ZZ(lambda, mass, v, cosXi);
129  if ( channel == "WW" ) return sv_WW(lambda, mass, v, cosXi);
130  }
131  return 0;
132  }
133 
134  // Annihilation into W bosons.
135  double sv_WW(double lambda, double mass, double v, double cosXi)
136  {
137  double s = 4*mass*mass/(1-v*v/4);
138  double cos2Xi = cosXi*cosXi;
139  double sin2Xi = 1 - cos2Xi;
140  double numerator = (cos2Xi*v*v/4 + sin2Xi);
141  double x = pow(mW,2)/s;
142  double GeV2tocm3s1 = gev2cm2*s2cm;
143  return pow(lambda,2)*pow(s,2)/16/M_PI*sqrt(1-4*x)*Dh2(s)*numerator*(1-4*x+12*pow(x,2))
144  *GeV2tocm3s1;
145  }
146 
147  // Annihilation into Z bosons.
148  double sv_ZZ(double lambda, double mass, double v, double cosXi)
149  {
150  double s = 4*mass*mass/(1-v*v/4);
151  double cos2Xi = cosXi*cosXi;
152  double sin2Xi = 1 - cos2Xi;
153  double numerator = (cos2Xi*v*v/4 + sin2Xi);
154  double x = pow(mZ0,2)/s;
155  double GeV2tocm3s1 = gev2cm2*s2cm;
156  return pow(lambda,2)*pow(s,2)/32/M_PI*sqrt(1-4*x)*Dh2(s)*numerator*(1-4*x+12*pow(x,2))
157  *GeV2tocm3s1;
158  }
159 
160  // Annihilation into fermions
161  double sv_ff(
162  double lambda, double mass, double v, double mf, double cosXi, bool is_quark)
163  {
164  double s = 4*mass*mass/(1-v*v/4);
165  double cos2Xi = cosXi*cosXi;
166  double sin2Xi = 1 - cos2Xi;
167  double numerator = (cos2Xi*v*v/4 + sin2Xi);
168  double vf = sqrt(1-4*pow(mf,2)/s);
169  double Xf = 1;
170  if ( is_quark ) Xf = 3 *
171  (1+(3/2*log(pow(mf,2)/s)+9/4)*4*alpha_s/3/M_PI);
172  double GeV2tocm3s1 = gev2cm2*s2cm;
173  return pow(lambda,2)*s*
174  pow(mf,2)/8/M_PI*Xf*pow(vf,3)*Dh2(s)*numerator*GeV2tocm3s1;
175  }
176 
178  double sv_hh(double lambda, double mass, double v, double cosXi)
179  {
180  // Hardcoded velocity avoids NAN results.
181  v = std::max(v, 1e-6);
182 
183  double s = 4*mass*mass/(1-v*v/4); // v is relative velocity
184  double GeV2tocm3s1 = gev2cm2*s2cm;
185  double xh = mh*mh/s;
186  double xpsi = mass*mass/s;
187  double xG = Gamma_mh*mh/s;
188 
189  double beta = (s - 2*pow(mh,2))/sqrt((s - 4*pow(mh,2))*(s - 4*pow(mass,2)));
190 
191  return (pow(lambda,2)*sqrt(1 - 4*xh)/(32.*M_PI*s)*(
192  s - 4*pow(cosXi,2)*s*xpsi - 8*cosXi*lambda*pow(v0,2)*mass +
193  (3*xh*(8*cosXi*lambda*pow(v0,2)*(-1 + xh)*sqrt(s*xpsi) - s*(2 + xh)*(-1 + 4*pow(cosXi,2)*xpsi)))/(pow(xG,2) + pow(-1 + xh,2))
194  - (2*pow(lambda,2)*pow(v0,4)*(3*pow(xh,2) - 8*(1 + pow(cosXi,2))*xh*xpsi + 2*xpsi*(1 + 8*pow(cosXi,4)*xpsi)))/(pow(xh,2) + xpsi - 4*xh*xpsi)
195  + (4*beta*lambda*pow(v0,2)*(2*cosXi*(-1 + 2*xh)*(-1 - pow(xG,2) + xh*(-1 + 2*xh))*sqrt(s*xpsi)*(-1 - 2*xh + 8*pow(cosXi,2)*xpsi) + lambda*pow(v0,2)*(pow(xG,2) + pow(-1 + xh,2))*
196  (1 - 4*xh + 6*pow(xh,2) - 16*pow(cosXi,2)*(-1 + xh)*xpsi - 32*pow(cosXi,4)*pow(xpsi,2)))*atanh(1/beta))/((pow(xG,2) + pow(-1 + xh,2))*pow(1 - 2*xh,2)))
197  )*GeV2tocm3s1;
198 
199  }
200 
201  private:
202  double Gamma_mh, mh, v0, alpha_s, mb, mc, mtau, mt, mZ0, mW;
203  };
204 
205  void DarkMatter_ID_MajoranaSingletDM(std::string & result) { result = "X"; }
206 
208  void DD_couplings_MajoranaSingletDM_Z2(DM_nucleon_couplings_fermionic_HP &result)
209  {
212  const SubSpectrum& he = spec.get_HE();
213  //double mass = spec.get(Par::Pole_Mass,"X");
214  double lambda = he.get(Par::dimensionless,"lX");
215  double cosXI = std::cos(he.get(Par::dimensionless,"xi"));
216  double sinXI = std::sin(he.get(Par::dimensionless,"xi"));
217  double mh = spec.get(Par::Pole_Mass,"h0_1");
218 
219  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
220  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
221  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
222 
223  // SI scalar and pseudoscalar couplings
224  result.gps = lambda*fp*m_proton*cosXI/pow(mh,2);
225  result.gns = lambda*fn*m_neutron*cosXI/pow(mh,2);
226  result.gp_q2 = lambda*fp*m_proton*sinXI/pow(mh,2);
227  result.gn_q2 = lambda*fn*m_neutron*sinXI/pow(mh,2);
228 
229  logger() << LogTags::debug << "Majorana DM DD couplings:" << std::endl;
230  logger() << " gps = " << result.gps << std::endl;
231  logger() << " gns = " << result.gns << std::endl;
232  logger() << " gp_q2 = " << result.gp_q2 << std::endl;
233  logger() << " gn_q2 = " << result.gn_q2 << EOM;
234 
235  } // function DD_couplings_MajoranaSingletDM_Z2
236 
239  {
241  using std::vector;
242  using std::string;
243 
244  // Initialize empty catalog
245  TH_ProcessCatalog catalog;
246  TH_Process process_ann("X", "X");
247 
248  // Explicitly state that Majorana DM is self-conjugate
249  process_ann.isSelfConj = true;
250 
252  // Import particle masses and couplings
254 
255  // Convenience macros
256  #define getSMmass(Name, spinX2) \
257  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
258  (Name , TH_ParticleProperty(SM.get(Par::Pole_Mass,Name), spinX2)));
259  #define addParticle(Name, Mass, spinX2) \
260  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
261  (Name , TH_ParticleProperty(Mass, spinX2)));
262 
263  // Import Spectrum objects
265  const SubSpectrum& he = spec.get_HE();
266  const SubSpectrum& SM = spec.get_LE();
267  const SMInputs& SMI = spec.get_SMInputs();
268 
269  // Import couplings
270  double lambda = he.get(Par::dimensionless,"lX");
271  double v = he.get(Par::mass1,"vev");
272  double alpha_s = SMI.alphaS; // alpha_s(mZ)^MSbar
273  double cosXi = std::cos(he.get(Par::dimensionless, "xi"));
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 
302  // Masses for neutrino flavour eigenstates. Set to zero.
303  // (presently not required)
304  addParticle("nu_e", 0.0, 1)
305  addParticle("nubar_e", 0.0, 1)
306  addParticle("nu_mu", 0.0, 1)
307  addParticle("nubar_mu", 0.0, 1)
308  addParticle("nu_tau", 0.0, 1)
309  addParticle("nubar_tau",0.0, 1)
310 
311  // Higgs-sector masses
312  double mX = spec.get(Par::Pole_Mass,"X");
313  double mH = spec.get(Par::Pole_Mass,"h0_1");
314  addParticle("X", mX, 1) // Majorana DM
315  addParticle("h0_1", mH, 0) // SM-like Higgs
316  addParticle("pi0", meson_masses.pi0, 0)
317  addParticle("pi+", meson_masses.pi_plus, 0)
318  addParticle("pi-", meson_masses.pi_minus, 0)
319  addParticle("eta", meson_masses.eta, 0)
320  addParticle("rho0", meson_masses.rho0, 1)
321  addParticle("rho+", meson_masses.rho_plus, 1)
322  addParticle("rho-", meson_masses.rho_minus, 1)
323  addParticle("omega", meson_masses.omega, 1)
324 
325  // Get rid of convenience macros
326  #undef getSMmass
327  #undef addParticle
328 
330  // Import Decay information
332 
333  // Import decay table from DecayBit
334  const DecayTable* tbl = &(*Dep::decay_rates);
335 
336  // Save Higgs width for later
337  double gammaH = tbl->at("h0_1").width_in_GeV;
338 
339  // Set of imported decays
340  std::set<string> importedDecays;
341 
342  // Minimum branching ratio to include
343  double minBranching = 0;
344 
345  // Import relevant decays (only Higgs and subsequent decays)
347  // Notes: Virtual Higgs decays into offshell W+W- final states are not
348  // imported. All other channels are correspondingly rescaled. Decay
349  // into XX final states is accounted for, leading to zero photons.
350  ImportDecays("h0_1", catalog, importedDecays, tbl, minBranching,
351  daFunk::vec<std::string>("Z0", "W+", "W-", "e+_2", "e-_2", "e+_3", "e-_3"));
352 
353  // Instantiate new MajoranaSingletDM_Z2 object
354  auto majoranaDM = boost::make_shared<MajoranaSingletDM>(&catalog, gammaH, v, alpha_s);
355 
356  // Populate annihilation channel list and add thresholds to threshold
357  // list.
358  // (remark: the lowest threshold is here = 2*mX, whereas in DS-internal
359  // conventions, this lowest threshold is not listed)
360  process_ann.resonances_thresholds.threshold_energy.push_back(2*mX);
361  auto channel =
362  daFunk::vec<string>("bb", "WW", "cc", "tautau", "ZZ", "tt", "hh");
363  auto p1 =
364  daFunk::vec<string>("d_3", "W+", "u_2", "e+_3", "Z0", "u_3", "h0_1");
365  auto p2 =
366  daFunk::vec<string>("dbar_3","W-", "ubar_2","e-_3", "Z0", "ubar_3","h0_1");
367  {
368  for ( unsigned int i = 0; i < channel.size(); i++ )
369  {
370  double mtot_final =
371  catalog.getParticleProperty(p1[i]).mass +
372  catalog.getParticleProperty(p2[i]).mass;
373  // Include final states that are open for T~m/20
374  if ( mX*2 > mtot_final*0.5 )
375  {
376  daFunk::Funk kinematicFunction = daFunk::funcM(majoranaDM,
377  &MajoranaSingletDM::sv, channel[i], lambda, mX, cosXi, daFunk::var("v"));
378  TH_Channel new_channel(
379  daFunk::vec<string>(p1[i], p2[i]), kinematicFunction
380  );
381  process_ann.channelList.push_back(new_channel);
382  }
383  if ( mX*2 > mtot_final )
384  {
386  push_back(mtot_final);
387  }
388  }
389  }
390 
391  // Populate resonance list
392  if ( mH >= mX*2 ) process_ann.resonances_thresholds.resonances.
393  push_back(TH_Resonance(mH, gammaH));
394 
395  // Add process to previous list
396  catalog.processList.push_back(process_ann);
397 
398  // Validate
399  catalog.validate();
400 
401  // Return the finished process catalog
402  result = catalog;
403 
404  } // function TH_ProcessCatalog_MajoranaSingletDM_Z2
405  }
406 }
std::vector< TH_Channel > channelList
List of channels.
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.
Spectrum Spectrum Spectrum Spectrum Spectrum MajoranaSingletDM_Z2_spectrum
#define addParticle(Name, Mass, spinX2)
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.
START_MODEL beta
Definition: Axions.hpp:35
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.
const double mt
Definition: topness.h:39
START_MODEL v0
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 DarkMatter_ID_MajoranaSingletDM(std::string &result)
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...
A single resonance of a given width at a given energy (both in GeV)
void TH_ProcessCatalog_MajoranaSingletDM_Z2(DarkBit::TH_ProcessCatalog &result)
Set up process catalog for the MajoranaSingletDM_Z2 model.
virtual double get(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
TH_resonances_thresholds resonances_thresholds
List of resonances and thresholds.
double width_in_GeV
Total particle width (in GeV)
#define getSMmass(Name, spinX2)
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
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.
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
TODO: see if we can use this one:
Definition: Analysis.hpp:33
void DD_couplings_MajoranaSingletDM_Z2(DM_nucleon_couplings_fermionic_HP &result)
Direct detection couplings for the MajoranaSingletDM_Z2 model.
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