gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
VectorSingletDM.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 VectorSingletDM
41  {
42  public:
44  VectorSingletDM(
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  ~VectorSingletDM() {}
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  // Note: Valid for mass > 45 GeV
77  double s = 4*mass*mass/(1-v*v/4);
78  double sqrt_s = sqrt(s);
79  if ( sqrt_s < 90 )
80  {
82  "VectorSingletDM sigmav called with sqrt_s < 90 GeV.");
83  return 0;
84  }
85 
86  if ( channel == "hh" )
87  {
88  if ( sqrt_s > mh*2 )
89  {
90  return sv_hh(lambda, mass, v);
91  }
92  else return 0;
93  }
94 
95  if ( channel == "bb" and sqrt_s < mb*2 ) return 0;
96  if ( channel == "cc" and sqrt_s < mc*2 ) return 0;
97  if ( channel == "tautau" and sqrt_s < mtau*2 ) return 0;
98  if ( channel == "tt" and sqrt_s < mt*2 ) return 0;
99  if ( channel == "ZZ" and sqrt_s < mZ0*2) return 0;
100  if ( channel == "WW" and sqrt_s < mW*2) return 0;
101 
102  if ( sqrt_s < 300 )
103  {
104  double br = virtual_SMHiggs_widths(channel,sqrt_s);
105  double Gamma_s = virtual_SMHiggs_widths("Gamma",sqrt_s);
106  double GeV2tocm3s1 = gev2cm2*s2cm;
107  double y = s/pow(mass, 2);
108 
109  // Explicitly close channel for off-shell top quarks
110  if ( channel == "tt" and sqrt_s < mt*2) return 0;
111 
112  double res = 2*lambda*lambda*v0*v0*(1-y/3 + pow(y,2)/12)/3/
113  sqrt_s*Dh2(s)*Gamma_s*GeV2tocm3s1*br;
114  return res;
115  }
116  else
117  {
118  if ( channel == "bb" ) return sv_ff(lambda, mass, v, mb, true);
119  if ( channel == "cc" ) return sv_ff(lambda, mass, v, mc, true);
120  if ( channel == "tautau" ) return sv_ff(lambda, mass, v, mtau, false);
121  if ( channel == "tt" ) return sv_ff(lambda, mass, v, mt, true);
122  if ( channel == "ZZ" ) return sv_ZZ(lambda, mass, v);
123  if ( channel == "WW" ) return sv_WW(lambda, mass, v);
124  }
125  return 0;
126  }
127 
128  // Annihilation into W bosons.
129  double sv_WW(double lambda, double mass, double v)
130  {
131  double s = 4*mass*mass/(1-v*v/4);
132  double x = pow(mW,2)/s;
133  double y = s/pow(mass, 2);
134  double GeV2tocm3s1 = gev2cm2*s2cm;
135  return pow(lambda,2)*s/24/M_PI*sqrt(1-4*x)*(1-y/3 + pow(y,2)/12)*
136  Dh2(s)*(1-4*x+12*pow(x,2))*GeV2tocm3s1;
137  }
138 
139  // Annihilation into Z bosons.
140  double sv_ZZ(double lambda, double mass, double v)
141  {
142  double s = 4*mass*mass/(1-v*v/4);
143  double x = pow(mZ0,2)/s;
144  double y = s/pow(mass, 2);
145  double GeV2tocm3s1 = gev2cm2*s2cm;
146  return pow(lambda,2)*s/48/M_PI*sqrt(1-4*x)*(1-y/3 + pow(y,2)/12)*
147  Dh2(s)*(1-4*x+12*pow(x,2))*GeV2tocm3s1;
148  }
149 
150  // Annihilation into fermions
151  double sv_ff(
152  double lambda, double mass, double v, double mf, bool is_quark)
153  {
154  double s = 4*mass*mass/(1-v*v/4);
155  double vf = sqrt(1-4*pow(mf,2)/s);
156  double y = s/pow(mass, 2);
157  double Xf = 1;
158  if ( is_quark ) Xf = 3 *
159  (1+(3/2*log(pow(mf,2)/s)+9/4)*4*alpha_s/3/M_PI);
160  double GeV2tocm3s1 = gev2cm2*s2cm;
161  return pow(lambda,2)*(1-y/3 + pow(y,2)/12)*
162  pow(mf,2)/12/M_PI*Xf*pow(vf,3) * Dh2(s) *GeV2tocm3s1;
163  }
164 
166  double sv_hh(double lambda, double mass, double v)
167  {
168  // Hardcode lower limit for velocity to avoid nan results.
169  v = std::max(v, 1e-6);
170  double GeV2tocm3s1 = gev2cm2*s2cm;
171  double s = 4*mass*mass/(1-v*v/4); // v is relative velocity
172  double xV = mass*mass/s;
173  double xH = mh*mh/s;
174  double xG = Gamma_mh*mh/s;
175  double beta = (s - 2*pow(mh,2))/sqrt((s - 4*pow(mh,2))*(s - 4*pow(mass,2)));
176  double vH = sqrt(1-4*xH);
177 
178  return GeV2tocm3s1*(vH*pow(lambda,2)*pow(s,-3)*pow(xV,-4)*pow(pow(xG,2) + pow(-1 + xH,2),-1)*
179  ((4*lambda*s*xV*(1 - 6*xV + xH*(2 + 4*xV))*pow(v0,2)*(1 + xH + pow(xG,2) - 2*pow(xH,2))*(xV - 4*xH*xV + pow(xH,2)) +
180  4*(1 + 4*xV*(-1 + 3*xV))*pow(s,2)*(xV - 4*xH*xV + pow(xH,2))*(pow(xG,2) + pow(1 + 2*xH,2))*pow(xV,2) +
181  pow(lambda,2)*pow(v0,4)*(pow(xG,2) + pow(-1 + xH,2))*(xV + (1 + 12*xV*(-1 + 4*xV))*pow(xH,2) + (4 - 32*xV)*pow(xH,3) + 6*pow(xH,4) - 64*xH*pow(xV,3) +
182  96*pow(xV,4)))*pow(xV - 4*xH*xV + pow(xH,2),-1) - 8*beta*lambda*(-1 + 4*xH)*atanh(1/beta)*pow(v0,2)*pow(1 - 2*xH,-1)*
183  (2*s*(-1 + 2*xH)*xV*((-1 + xH)*(1 + 2*xH) - pow(xG,2))*(2*xV*pow(-1 + xH,2) + pow(xH,2) - 4*(1 + 2*xH)*pow(xV,2) + 24*pow(xV,3)) -
184  lambda*pow(v0,2)*(pow(xG,2) + pow(-1 + xH,2))*(-8*xV*pow(xH,3) + 3*pow(xH,4) - pow(xH,2)*(1 + 8*pow(xV,2)) + 4*xH*(xV + 8*pow(xV,3)) -
185  2*xV*(1 - 2*xV + 24*pow(xV,3))))*pow(1 - 6*xH + 8*pow(xH,2),-1)))/(2304*M_PI);
186  }
187 
188  private:
189  double Gamma_mh, mh, v0, alpha_s, mb, mc, mtau, mt, mZ0, mW;
190  };
191 
192  void DarkMatter_ID_VectorSingletDM(std::string& result) { result = "V"; }
193 
195  void DD_couplings_VectorSingletDM_Z2(DM_nucleon_couplings &result)
196  {
199  const SubSpectrum& he = spec.get_HE();
200  double mass = spec.get(Par::Pole_Mass,"V");
201  double lambda = he.get(Par::dimensionless,"lambda_hV");
202  double mh = spec.get(Par::Pole_Mass,"h0_1");
203 
204  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
205  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
206  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
207 
208  result.gps = lambda*fp*m_proton/pow(mh,2)/mass/2;
209  result.gns = lambda*fn*m_neutron/pow(mh,2)/mass/2;
210  result.gpa = 0; // Only SI cross-section
211  result.gna = 0;
212 
213  logger() << LogTags::debug << "Vector DM DD couplings:" << std::endl;
214  logger() << " gps = " << result.gps << std::endl;
215  logger() << " gns = " << result.gns << std::endl;
216  logger() << " gpa = " << result.gpa << std::endl;
217  logger() << " gna = " << result.gna << EOM;
218 
219  } // function DD_couplings_VectorSingletDM_Z2
220 
223  {
225  using std::vector;
226  using std::string;
227 
228  // Initialize empty catalog
229  TH_ProcessCatalog catalog;
230  TH_Process process_ann("V", "V");
231 
232  // Explicitly state that Vector DM is self-conjugate
233  process_ann.isSelfConj = true;
234 
236  // Import particle masses and couplings
238 
239  // Convenience macros
240  #define getSMmass(Name, spinX2) \
241  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
242  (Name , TH_ParticleProperty(SM.get(Par::Pole_Mass,Name), spinX2)));
243  #define addParticle(Name, Mass, spinX2) \
244  catalog.particleProperties.insert(std::pair<string, TH_ParticleProperty> \
245  (Name , TH_ParticleProperty(Mass, spinX2)));
246 
247  // Import Spectrum objects
249  const SubSpectrum& he = spec.get_HE();
250  const SubSpectrum& SM = spec.get_LE();
251  const SMInputs& SMI = spec.get_SMInputs();
252 
253  // Import couplings
254  double lambda = he.get(Par::dimensionless,"lambda_hV");
255  double v = he.get(Par::mass1,"vev");
256  double alpha_s = SMI.alphaS; // alpha_s(mZ)^MSbar
257 
258  // Get SM pole masses
259  getSMmass("e-_1", 1)
260  getSMmass("e+_1", 1)
261  getSMmass("e-_2", 1)
262  getSMmass("e+_2", 1)
263  getSMmass("e-_3", 1)
264  getSMmass("e+_3", 1)
265  getSMmass("Z0", 2)
266  getSMmass("W+", 2)
267  getSMmass("W-", 2)
268  getSMmass("g", 2)
269  getSMmass("gamma", 2)
270  getSMmass("u_3", 1)
271  getSMmass("ubar_3", 1)
272  getSMmass("d_3", 1)
273  getSMmass("dbar_3", 1)
274 
275  // Pole masses not available for the light quarks.
276  addParticle("u_1" , SMI.mU, 1) // mu(2 GeV)^MS-bar, not pole mass
277  addParticle("ubar_1", SMI.mU, 1) // mu(2 GeV)^MS-bar, not pole mass
278  addParticle("d_1" , SMI.mD, 1) // md(2 GeV)^MS-bar, not pole mass
279  addParticle("dbar_1", SMI.mD, 1) // md(2 GeV)^MS-bar, not pole mass
280  addParticle("u_2" , SMI.mCmC,1) // mc(mc)^MS-bar, not pole mass
281  addParticle("ubar_2", SMI.mCmC,1) // mc(mc)^MS-bar, not pole mass
282  addParticle("d_2" , SMI.mS, 1) // ms(2 GeV)^MS-bar, not pole mass
283  addParticle("dbar_2", SMI.mS, 1) // ms(2 GeV)^MS-bar, not pole mass
284 
285  // Masses for neutrino flavour eigenstates. Set to zero.
286  // (presently not required)
287  addParticle("nu_e", 0.0, 1)
288  addParticle("nubar_e", 0.0, 1)
289  addParticle("nu_mu", 0.0, 1)
290  addParticle("nubar_mu", 0.0, 1)
291  addParticle("nu_tau", 0.0, 1)
292  addParticle("nubar_tau",0.0, 1)
293 
294  // Higgs-sector masses
295  double mV = spec.get(Par::Pole_Mass,"V");
296  double mH = spec.get(Par::Pole_Mass,"h0_1");
297  addParticle("V", mV, 2) // Vector DM
298  addParticle("h0_1", mH, 0) // SM-like Higgs
299  addParticle("pi0", meson_masses.pi0, 0)
300  addParticle("pi+", meson_masses.pi_plus, 0)
301  addParticle("pi-", meson_masses.pi_minus, 0)
302  addParticle("eta", meson_masses.eta, 0)
303  addParticle("rho0", meson_masses.rho0, 1)
304  addParticle("rho+", meson_masses.rho_plus, 1)
305  addParticle("rho-", meson_masses.rho_minus, 1)
306  addParticle("omega", meson_masses.omega, 1)
307 
308  // Get rid of convenience macros
309  #undef getSMmass
310  #undef addParticle
311 
313  // Import Decay information
315 
316  // Import decay table from DecayBit
317  const DecayTable* tbl = &(*Dep::decay_rates);
318 
319  // Save Higgs width for later
320  double gammaH = tbl->at("h0_1").width_in_GeV;
321 
322  // Set of imported decays
323  std::set<string> importedDecays;
324 
325  // Minimum branching ratio to include
326  double minBranching = 0;
327 
328  // Import relevant decays (only Higgs and subsequent decays)
330  // Notes: Virtual Higgs decays into offshell W+W- final states are not
331  // imported. All other channels are correspondingly rescaled. Decay
332  // into VV final states is accounted for, leading to zero photons.
333  ImportDecays("h0_1", catalog, importedDecays, tbl, minBranching,
334  daFunk::vec<std::string>("Z0", "W+", "W-", "e+_2", "e-_2", "e+_3", "e-_3"));
335 
336  // Instantiate new VectorSingletDM_Z2 object
337  auto vectorDM = boost::make_shared<VectorSingletDM>(&catalog, gammaH, v, alpha_s);
338 
339  // Populate annihilation channel list and add thresholds to threshold
340  // list.
341  // (remark: the lowest threshold is here = 2*mV, whereas in DS-internal
342  // conventions, this lowest threshold is not listed)
343  process_ann.resonances_thresholds.threshold_energy.push_back(2*mV);
344  auto channel =
345  daFunk::vec<string>("bb", "WW", "cc", "tautau", "ZZ", "tt", "hh");
346  auto p1 =
347  daFunk::vec<string>("d_3", "W+", "u_2", "e+_3", "Z0", "u_3", "h0_1");
348  auto p2 =
349  daFunk::vec<string>("dbar_3","W-", "ubar_2","e-_3", "Z0", "ubar_3","h0_1");
350  {
351  for ( unsigned int i = 0; i < channel.size(); i++ )
352  {
353  double mtot_final =
354  catalog.getParticleProperty(p1[i]).mass +
355  catalog.getParticleProperty(p2[i]).mass;
356  // Include final states that are open for T~m/20
357  if ( mV*2 > mtot_final*0.5 )
358  {
359  daFunk::Funk kinematicFunction = daFunk::funcM(vectorDM,
360  &VectorSingletDM::sv, channel[i], lambda, mV, daFunk::var("v"));
361  TH_Channel new_channel(
362  daFunk::vec<string>(p1[i], p2[i]), kinematicFunction
363  );
364  process_ann.channelList.push_back(new_channel);
365  }
366  if ( mV*2 > mtot_final )
367  {
369  push_back(mtot_final);
370  }
371  }
372  }
373 
374  // Populate resonance list
375  if ( mH >= mV*2 ) process_ann.resonances_thresholds.resonances.
376  push_back(TH_Resonance(mH, gammaH));
377 
378  catalog.processList.push_back(process_ann);
379 
380  // Validate
381  catalog.validate();
382 
383  result = catalog;
384 
385  } // function TH_ProcessCatalog_VectorSingletDM_Z2
386  }
387 }
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.
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
void TH_ProcessCatalog_VectorSingletDM_Z2(DarkBit::TH_ProcessCatalog &result)
Set up process catalog for the VectorSingletDM_Z2 model.
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) ...
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)
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)
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
#define addParticle(Name, Mass, spinX2)
void DD_couplings_VectorSingletDM_Z2(DM_nucleon_couplings &result)
Direct detection couplings for the VectorSingletDM_Z2 model.
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
Spectrum Spectrum Spectrum Spectrum VectorSingletDM_Z2_spectrum
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
void DarkMatter_ID_VectorSingletDM(std::string &result)
A container for a single process.
#define getSMmass(Name, spinX2)
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
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