gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
SpecBit_ScalarSingletDM.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
19 
20 #include <string>
21 #include <sstream>
22 
24 
28 
35 
38 
39 // Flexible SUSY stuff (should not be needed by the rest of gambit)
40 #include "flexiblesusy/src/ew_input.hpp"
41 #include "flexiblesusy/src/lowe.h"
42 #include "flexiblesusy/src/numerics2.hpp"
43 #include "flexiblesusy/src/spectrum_generator_settings.hpp"
44 // Switch for debug mode
45 //#define SPECBIT_DEBUG
46 
47 namespace Gambit
48 {
49 
50  namespace SpecBit
51  {
52  using namespace LogTags;
53  using namespace flexiblesusy;
54 
57  {
58  namespace myPipe = Pipes::get_ScalarSingletDM_Z2_spectrum;
59  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
60 
61  // Initialise an object to carry the Singlet plus Higgs sector information
63 
64  // quantities needed to fill container spectrum, intermediate calculations
65  double alpha_em = 1.0 / sminputs.alphainv;
66  double C = alpha_em * pi / (sminputs.GF * pow(2,0.5));
67  double sinW2 = 0.5 - pow( 0.25 - C/pow(sminputs.mZ,2) , 0.5);
68  double cosW2 = 0.5 + pow( 0.25 - C/pow(sminputs.mZ,2) , 0.5);
69  double e = pow( 4*pi*( alpha_em ),0.5) ;
70 
71  // Higgs sector
72  double mh = *myPipe::Param.at("mH");
73  singletmodel.HiggsPoleMass = mh;
74 
75  double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
76  singletmodel.HiggsVEV = vev;
77 
78  // Scalar singlet sector
79  singletmodel.SingletPoleMass = *myPipe::Param.at("mS");
80  singletmodel.SingletLambda = *myPipe::Param.at("lambda_hS");
81  singletmodel.SingletLambdaS = 0;
82 
83  // Standard model
84  singletmodel.sinW2 = sinW2;
85 
86  // gauge couplings
87  singletmodel.g1 = sqrt(5/3) * e / sqrt(cosW2);
88  singletmodel.g2 = e / sqrt(sinW2);
89  singletmodel.g3 = pow( 4*pi*( sminputs.alphaS ),0.5) ;
90 
91  // Yukawas
92  double sqrt2v = pow(2.0,0.5)/vev;
93  singletmodel.Yu[0] = sqrt2v * sminputs.mU;
94  singletmodel.Yu[1] = sqrt2v * sminputs.mCmC;
95  singletmodel.Yu[2] = sqrt2v * sminputs.mT;
96  singletmodel.Ye[0] = sqrt2v * sminputs.mE;
97  singletmodel.Ye[1] = sqrt2v * sminputs.mMu;
98  singletmodel.Ye[2] = sqrt2v * sminputs.mTau;
99  singletmodel.Yd[0] = sqrt2v * sminputs.mD;
100  singletmodel.Yd[1] = sqrt2v * sminputs.mS;
101  singletmodel.Yd[2] = sqrt2v * sminputs.mBmB;
102 
103  // Create a SubSpectrum object to wrap the EW sector information
104  Models::ScalarSingletDM_Z2SimpleSpec singletspec(singletmodel);
105 
106  // Retrieve any mass cuts
107  static const Spectrum::mc_info mass_cut = myPipe::runOptions->getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
108  static const Spectrum::mr_info mass_ratio_cut = myPipe::runOptions->getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
109 
110  // We don't supply a LE subspectrum here; an SMSimpleSpec will therefore be automatically created from 'sminputs'
111  result = Spectrum(singletspec,sminputs,&myPipe::Param,mass_cut,mass_ratio_cut);
112 
113  }
114 
117  {
118  namespace myPipe = Pipes::get_ScalarSingletDM_Z3_spectrum;
119  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
120 
121  // Initialise an object to carry the Singlet plus Higgs sector information
122  Models::ScalarSingletDM_Z3Model singletmodel;
123 
124  // quantities needed to fill container spectrum, intermediate calculations
125  double alpha_em = 1.0 / sminputs.alphainv;
126  double C = alpha_em * pi / (sminputs.GF * pow(2,0.5));
127  double sinW2 = 0.5 - pow( 0.25 - C/pow(sminputs.mZ,2) , 0.5);
128  double cosW2 = 0.5 + pow( 0.25 - C/pow(sminputs.mZ,2) , 0.5);
129  double e = pow( 4*pi*( alpha_em ),0.5) ;
130 
131  // Higgs sector
132  double mh = *myPipe::Param.at("mH");
133  singletmodel.HiggsPoleMass = mh;
134 
135  double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
136  singletmodel.HiggsVEV = vev;
137 
138  // Scalar singlet sector
139  singletmodel.SingletPoleMass = *myPipe::Param.at("mS");
140  singletmodel.SingletLambda = *myPipe::Param.at("lambda_hS");
141  singletmodel.SingletLambdaS = *myPipe::Param.at("lambda_S");
142 
143  singletmodel.mu3 = *myPipe::Param.at("mu3");
144 
145  // Standard model
146  singletmodel.sinW2 = sinW2;
147 
148  // gauge couplings
149  singletmodel.g1 = sqrt(5/3) * e / sqrt(cosW2);
150  singletmodel.g2 = e / sqrt(sinW2);
151  singletmodel.g3 = pow( 4*pi*( sminputs.alphaS ),0.5) ;
152 
153  // Yukawas
154  double sqrt2v = pow(2.0,0.5)/vev;
155  singletmodel.Yu[0] = sqrt2v * sminputs.mU;
156  singletmodel.Yu[1] = sqrt2v * sminputs.mCmC;
157  singletmodel.Yu[2] = sqrt2v * sminputs.mT;
158  singletmodel.Ye[0] = sqrt2v * sminputs.mE;
159  singletmodel.Ye[1] = sqrt2v * sminputs.mMu;
160  singletmodel.Ye[2] = sqrt2v * sminputs.mTau;
161  singletmodel.Yd[0] = sqrt2v * sminputs.mD;
162  singletmodel.Yd[1] = sqrt2v * sminputs.mS;
163  singletmodel.Yd[2] = sqrt2v * sminputs.mBmB;
164 
165  // Create a SubSpectrum object to wrap the EW sector information
166  Models::ScalarSingletDM_Z3SimpleSpec singletspec(singletmodel);
167 
168  // Retrieve any mass cuts
169  static const Spectrum::mc_info mass_cut = myPipe::runOptions->getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
170  static const Spectrum::mr_info mass_ratio_cut = myPipe::runOptions->getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
171 
172  // We don't supply a LE subspectrum here; an SMSimpleSpec will therefore be automatically created from 'sminputs'
173  result = Spectrum(singletspec,sminputs,&myPipe::Param,mass_cut,mass_ratio_cut);
174 
175  }
176 
177 
178  template<class MI,class SI>
180  ( const typename MI::InputParameters& input
181  , const SMInputs& sminputs
182  , const Options& runOptions
183  , const std::map<str, safe_ptr<const double> >& input_Param
184  )
185  {
186  // SoftSUSY object used to set quark and lepton masses and gauge
187  // couplings in QEDxQCD effective theory
188  // Will be initialised by default using values in lowe.h, which we will
189  // override next.
190  softsusy::QedQcd oneset;
191 
192  // Fill QedQcd object with SMInputs values
193  setup_QedQcd(oneset,sminputs);
194 
195  // Run everything to Mz
196  oneset.toMz();
197 
198  // Create spectrum generator object
199  typename MI::SpectrumGenerator spectrum_generator;
200 
201  // Spectrum generator settings
202  // Default options copied from flexiblesusy/src/spectrum_generator_settings.hpp
203  //
204  // | enum | possible values | default value |
205  // |----------------------------------|------------------------------|-----------------|
206  // | precision | any positive double | 1.0e-4 |
207  // | max_iterations | any positive double | 0 (= automatic) |
208  // | algorithm | 0 (two-scale) or 1 (lattice) | 0 (= two-scale) |
209  // | calculate_sm_masses | 0 (no) or 1 (yes) | 0 (= no) |
210  // | pole_mass_loop_order | 0, 1, 2 | 2 (= 2-loop) |
211  // | ewsb_loop_order | 0, 1, 2 | 2 (= 2-loop) |
212  // | beta_loop_order | 0, 1, 2 | 2 (= 2-loop) |
213  // | threshold_corrections_loop_order | 0, 1 | 1 (= 1-loop) |
214  // | higgs_2loop_correction_at_as | 0, 1 | 1 (= enabled) |
215  // | higgs_2loop_correction_ab_as | 0, 1 | 1 (= enabled) |
216  // | higgs_2loop_correction_at_at | 0, 1 | 1 (= enabled) |
217  // | higgs_2loop_correction_atau_atau | 0, 1 | 1 (= enabled) |
218 
219  Spectrum_generator_settings settings;
220  settings.set(Spectrum_generator_settings::precision, runOptions.getValueOrDef<double>(1.0e-4,"precision_goal"));
221  settings.set(Spectrum_generator_settings::max_iterations, runOptions.getValueOrDef<double>(0,"max_iterations"));
222  settings.set(Spectrum_generator_settings::calculate_sm_masses, runOptions.getValueOrDef<bool> (true, "calculate_sm_masses"));
223  settings.set(Spectrum_generator_settings::pole_mass_loop_order, runOptions.getValueOrDef<int>(2,"pole_mass_loop_order"));
224  settings.set(Spectrum_generator_settings::pole_mass_loop_order, runOptions.getValueOrDef<int>(2,"ewsb_loop_order"));
225  settings.set(Spectrum_generator_settings::beta_loop_order, runOptions.getValueOrDef<int>(2,"beta_loop_order"));
226  settings.set(Spectrum_generator_settings::threshold_corrections_loop_order, runOptions.getValueOrDef<int>(2,"threshold_corrections_loop_order"));
227  settings.set(Spectrum_generator_settings::higgs_2loop_correction_at_as, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_at_as"));
228  settings.set(Spectrum_generator_settings::higgs_2loop_correction_ab_as, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_ab_as"));
229  settings.set(Spectrum_generator_settings::higgs_2loop_correction_at_at, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_at_at"));
230  settings.set(Spectrum_generator_settings::higgs_2loop_correction_atau_atau, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_atau_atau"));
231 
232  spectrum_generator.set_settings(settings);
233 
234  // Generate spectrum
235  spectrum_generator.run(oneset, input);
236  const typename MI::Problems& problems = spectrum_generator.get_problems();
237  MI model_interface(spectrum_generator,oneset,input);
238  SI singletdmspec(model_interface, "FlexibleSUSY", "1.5.1");
239 
240  singletdmspec.set_override(Par::mass1,spectrum_generator.get_high_scale(),"high_scale",true);
241  singletdmspec.set_override(Par::mass1,spectrum_generator.get_susy_scale(),"susy_scale",true);
242  singletdmspec.set_override(Par::mass1,spectrum_generator.get_low_scale(), "low_scale", true);
243 
244  singletdmspec.set_override(Par::Pole_Mass_1srd_high,0.0, "h0_1", true);
245  singletdmspec.set_override(Par::Pole_Mass_1srd_low,0.0,"h0_1", true);
246 
247  // Create a second SubSpectrum object to wrap the qedqcd object used to initialise the spectrum generator
248  // Attach the sminputs object as well, so that SM pole masses can be passed on (these aren't easily
249  // extracted from the QedQcd object, so use the values that we put into it.)
250  QedQcdWrapper qedqcdspec(oneset,sminputs);
251 
252  // Deal with points where spectrum generator encountered a problem
253  #ifdef SPECBIT_DEBUG
254  std::cout<<"Problem? "<<problems.have_problem()<<std::endl;
255  std::cout<<"Warning? "<<problems.have_warning()<<std::endl;
256  #endif
257  if( problems.have_problem() )
258  {
259  if( runOptions.getValueOrDef<bool>(false,"invalid_point_fatal") )
260  {
261  std::ostringstream errmsg;
262  errmsg << "A serious problem was encountered during spectrum generation!";
263  errmsg << "Message from FlexibleSUSY:" << std::endl;
264  problems.print_problems(errmsg);
265  problems.print_warnings(errmsg);
266  SpecBit_error().raise(LOCAL_INFO,errmsg.str());
267  }
268  else
269  {
270  // Check what the problem was
271  // see: contrib/MassSpectra/flexiblesusy/src/problems.hpp
272  std::ostringstream msg;
273  problems.print_problems(msg);
274  invalid_point().raise(msg.str()); //TODO: This message isn't ending up in the logs.
275  }
276  }
277  double QEWSB = *input_Param.at("QEWSB");
278  singletdmspec.RunToScaleOverride(QEWSB);
279 
280  // Retrieve any mass cuts
281  static const Spectrum::mc_info mass_cut = runOptions.getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
282  static const Spectrum::mr_info mass_ratio_cut = runOptions.getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
283 
284  return Spectrum(qedqcdspec,singletdmspec,sminputs,&input_Param, mass_cut, mass_ratio_cut);
285 
286  }
287 
288 
289  template <class T>
290  void fill_ScalarSingletDM_input(T& input, const std::map<str, safe_ptr<const double> >& Param,SMInputs sminputs)//,double scale)
291  {
292  double mH = *Param.at("mH");
293  double mS = *Param.at("mS");
294  double lambda_hs = *Param.at("lambda_hS");
295  double lambda_s = *Param.at("lambda_S");
296 
297  double QEFT = mS;
298 
299  if (QEFT < sminputs.mT)
300  {
301  //QEFT = *Param.at("mZ");
302  QEFT = sminputs.mT;
303  }
304 
305  input.HiggsIN = -0.5*pow(mH,2);
306  double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
307 
308  input.muSInput = pow(mS,2)-0.5*lambda_hs*pow(vev,2);
309  input.LamSHInput = lambda_hs;
310  input.LamSInput = lambda_s;
311  input.QEWSB = QEFT; // scale where EWSB conditions are applied
312  input.Qin = QEFT; //scale; // highest scale at which model is run to
313 
314  }
315 
316  template <class T>
317  void fill_extra_input(T& input, const std::map<str, safe_ptr<const double> >& Param )
318  {
319  input.mu3Input=*Param.at("mu3");
320  }
321 
322  bool check_perturb(const Spectrum& spec, const std::vector<SpectrumParameter>& required_parameters, double scale, int pts)
323  {
324  std::unique_ptr<SubSpectrum> ScalarSingletDM = spec.clone_HE();
325  double step = log10(scale) / pts;
326  double runto;
327  double ul = 4.0 * pi;
328 
329  for (int i=0;i<pts;i++)
330  {
331  runto = pow(10,step*float(i+1.0)); // scale to run spectrum to
332  if (runto<100){runto=100.0;}// avoid running to low scales
333 
334  ScalarSingletDM -> RunToScale(runto);
335 
336  for(std::vector<SpectrumParameter>::const_iterator it = required_parameters.begin();
337  it != required_parameters.end(); ++it)
338  {
339  const Par::Tags tag = it->tag();
340  const std::string name = it->name();
341  const std::vector<int> shape = it->shape();
342  std::ostringstream label;
343  label << name <<" "<< Par::toString.at(tag);
344 
345  if (name == "lambda_S"){ul = pi;}
346  else if (name == "lambda_h"){ul = 2*pi;}
347  else if (name == "lambda_hS"){ul = 4*pi;}
348  else {ul = 100;}
349 
350  if(shape.size()==1 and shape[0]==1)
351  {
352  if (abs(ScalarSingletDM->get(tag,name))>ul)
353  {
354  return false;
355  }
356  }
357  else if(shape.size()==1 and shape[0]>1)
358  {
359  for(int k = 1; k<=shape[0]; ++k)
360  {
361  if (abs(ScalarSingletDM->get(tag,name,k))>ul) return false;
362  }
363  }
364  else if(shape.size()==2)
365  {
366  for(int k = 1; k<=shape[0]; ++k)
367  {
368  for(int j = 1; j<=shape[0]; ++j)
369  {
370  if (abs(ScalarSingletDM->get(tag,name,k,j))>ul) return false;
371  }
372  }
373  }
374 
375  // check stability condition
376 
377  double lamhs = ScalarSingletDM->get(Par::dimensionless,"lambda_hS");
378  double lams = ScalarSingletDM->get(Par::dimensionless,"lambda_S");
379  double lamh = ScalarSingletDM->get(Par::dimensionless,"lambda_h");
380 
381  double stability_condition = 2.0 * pow(0.5* 0.25 * lamh*lams,0.5) + 0.5*lamhs;
382 
383  if (!(stability_condition > 0) && (lams>0) && (lamh>0))
384  {
385  //cout << "EW stability condition violated at Q = " << scale <<" , lambda_hs = "<< lamhs << "lamh = " << lamh << " lams = " << lams << endl;
386  return false;
387  }
388 
389  }
390  }
391 
392  return true;
393 
394  }
395 
396  #if(FS_MODEL_ScalarSingletDM_Z2_IS_BUILT)
397  void get_ScalarSingletDM_Z2_spectrum_pole(Spectrum& result)
398  {
399  using namespace softsusy;
400  namespace myPipe = Pipes::get_ScalarSingletDM_Z2_spectrum_pole;
401  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
402  const Options& runOptions=*myPipe::runOptions;
403  ScalarSingletDM_Z2_input_parameters input;
404  fill_ScalarSingletDM_input(input,myPipe::Param,sminputs);
405  result = run_FS_spectrum_generator<ScalarSingletDM_Z2_interface<ALGORITHM1>,ScalarSingletDM_Z2Spec<ScalarSingletDM_Z2_interface<ALGORITHM1>>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
406 
407  int check_perturb_pts = runOptions.getValueOrDef<double>(10,"check_perturb_pts");
408  double do_check_perturb = runOptions.getValueOrDef<bool>(false,"check_perturb");
409  double check_perturb_scale = runOptions.getValueOrDef<double>(1.22e19,"check_high_scale");
410  double input_scale_tolerance = runOptions.getValueOrDef<double>(1e-3,"input_scale_tolerance");
411 
412  // Check that the Higgs MSbar parameter has been provided at the scale of the singlet MSbar mass
413  double Qin = *myPipe::Param.at("Qin");
414  if (abs((Qin - *myPipe::Param.at("mS"))/Qin) > input_scale_tolerance) SpecBit_error().raise(LOCAL_INFO, "SM_Higgs_running::Qin must equal ScalarSingletDM_Z2_running::mS");
415 
416  if (do_check_perturb)
417  {
418  static const SpectrumContents::ScalarSingletDM_Z2 contents;
419  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
420  if (!check_perturb(result,required_parameters,check_perturb_scale,check_perturb_pts))
421  {
422  // invalidate point as spectrum not perturbative up to scale
423  std::ostringstream msg;
424  msg << "Spectrum not perturbative up to scale = " << check_perturb_scale << std::endl;
425  #ifdef SPECBIT_DEBUG
426  cout << "Spectrum not perturbative up to scale = " << check_perturb_scale << endl;
427  #endif
428  invalid_point().raise(msg.str());
429  }
430  }
431 
432  }
433  #endif
434 
435  #if(FS_MODEL_ScalarSingletDM_Z3_IS_BUILT)
436  void get_ScalarSingletDM_Z3_spectrum_pole(Spectrum& result)
437  {
438  using namespace softsusy;
439  namespace myPipe = Pipes::get_ScalarSingletDM_Z3_spectrum_pole;
440  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
441  const Options& runOptions=*myPipe::runOptions;
442  //double scale = runOptions.getValueOrDef<double>(173.34,"FS_high_scale");
443  ScalarSingletDM_Z3_input_parameters input;
444  fill_ScalarSingletDM_input(input,myPipe::Param,sminputs);
445  fill_extra_input(input,myPipe::Param);
446  result = run_FS_spectrum_generator<ScalarSingletDM_Z3_interface<ALGORITHM1>,ScalarSingletDM_Z3Spec<ScalarSingletDM_Z3_interface<ALGORITHM1>>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
447 
448  int check_perturb_pts = runOptions.getValueOrDef<double>(10,"check_perturb_pts");
449  double do_check_perturb = runOptions.getValueOrDef<bool>(false,"check_perturb");
450  double check_perturb_scale = runOptions.getValueOrDef<double>(1.22e19,"check_high_scale");
451  double input_scale_tolerance = runOptions.getValueOrDef<double>(1e-3,"input_scale_tolerance");
452 
453  // Check that the Higgs MSbar parameter has been provided at the scale of the singlet MSbar mass
454  double Qin = *myPipe::Param.at("Qin");
455  if (abs((Qin - *myPipe::Param.at("mS"))/Qin) > input_scale_tolerance) SpecBit_error().raise(LOCAL_INFO, "SM_Higgs_running::Qin must equal ScalarSingletDM_Z3_running::mS");
456 
457  if (do_check_perturb)
458  {
459  static const SpectrumContents::ScalarSingletDM_Z3 contents;
460  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
461  if (!check_perturb(result,required_parameters,check_perturb_scale,check_perturb_pts))
462  {
463  // invalidate point as spectrum not perturbative up to scale
464  std::ostringstream msg;
465 
466  msg << "Spectrum not perturbative up to scale = " << check_perturb_scale << std::endl;
467  #ifdef SPECBIT_DEBUG
468  cout << "Spectrum not perturbative up to scale = " << check_perturb_scale << endl;
469  #endif
470  invalid_point().raise(msg.str());
471  }
472  }
473 
474  }
475  #endif
476 
477  #if(FS_MODEL_ScalarSingletDM_Z2_IS_BUILT)
478  void find_non_perturb_scale_ScalarSingletDM_Z2(double &result)
479  {
480  using namespace flexiblesusy;
481  using namespace softsusy;
482  namespace myPipe = Pipes::find_non_perturb_scale_ScalarSingletDM_Z2;
483  using namespace Gambit;
484  using namespace SpecBit;
485 
486  const Spectrum& fullspectrum = *myPipe::Dep::ScalarSingletDM_Z2_spectrum;
487 
488  // bound x by (a,b)
489 
490  double ms = *myPipe::Param.at("mS");
491 
492  double a = log10(ms);
493 
494  if (a > 20.0)
495  {
496  std::ostringstream msg;
497  msg << "Scalar mass larger than 10^20 GeV " << std::endl;
498  invalid_point().raise(msg.str());
499  }
500 
501  double b = 20.0;
502  double x = 0.5 * ( b + ms );
503 
504  while (abs(a-b)>1e-10)
505  {
506  x=0.5*(b-a)+a;
507  static const SpectrumContents::ScalarSingletDM_Z2 contents;
508  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
509  if (!check_perturb(fullspectrum,required_parameters,pow(10,x),3))
510  {
511  b=x;
512  }
513  else
514  {
515  a=x;
516  }
517  }
518  result = pow(10,0.5*(a+b));
519  }
520  #endif
521 
522  #if(FS_MODEL_ScalarSingletDM_Z3_IS_BUILT)
523  void find_non_perturb_scale_ScalarSingletDM_Z3(double &result)
524  {
525  using namespace flexiblesusy;
526  using namespace softsusy;
527  namespace myPipe = Pipes::find_non_perturb_scale_ScalarSingletDM_Z3;
528  using namespace Gambit;
529  using namespace SpecBit;
530 
531  const Spectrum& fullspectrum = *myPipe::Dep::ScalarSingletDM_Z3_spectrum;
532 
533  // bound x by (a,b)
534 
535  double ms = *myPipe::Param.at("mS");
536 
537  double a = log10(ms);
538 
539  if (a > 20.0)
540  {
541  std::ostringstream msg;
542  msg << "Scalar mass larger than 10^20 GeV " << std::endl;
543  invalid_point().raise(msg.str());
544  }
545 
546  double b = 20.0;
547  double x = 0.5 * ( b + ms );
548 
549  while (abs(a-b)>1e-10)
550  {
551  x=0.5*(b-a)+a;
552  static const SpectrumContents::ScalarSingletDM_Z3 contents;
553  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
554  if (!check_perturb(fullspectrum,required_parameters,pow(10,x),3))
555  {
556  b=x;
557  }
558  else
559  {
560  a=x;
561  }
562  }
563  result = pow(10,0.5*(a+b));
564  }
565  #endif
566 
569  {
571  dep_bucket<Spectrum>* spectrum_dependency = nullptr;
572  if (ModelInUse("ScalarSingletDM_Z2") or ModelInUse("ScalarSingletDM_Z2_running"))
573  {
574  spectrum_dependency = &Dep::ScalarSingletDM_Z2_spectrum;
575  }
576  else if (ModelInUse("ScalarSingletDM_Z3") or ModelInUse("ScalarSingletDM_Z3_running"))
577  {
578  spectrum_dependency = &Dep::ScalarSingletDM_Z3_spectrum;
579  }
580  else SpecBit_error().raise(LOCAL_INFO, "No valid model for ScalarSingletDM_higgs_couplings_pwid.");
581  const SubSpectrum& spec = (*spectrum_dependency)->get_HE();
582 
583  // Set the CP of the Higgs.
584  result.CP[0] = 1;
585  // Set the decays
586  result.set_neutral_decays_SM(0, "h0_1", *Dep::Reference_SM_Higgs_decay_rates);
587  result.set_neutral_decays(0, "h0_1", *Dep::Higgs_decay_rates);
588 
589  // Identify the singlet as the only possible invisible particle
590  if (spec.get(Par::Pole_Mass, "S") * 2.0 < spec.get(Par::Pole_Mass, "h0_1"))
591  result.invisibles = initVector<str>("S");
592  else
593  result.invisibles.clear();
594  // Leave all the effective couplings for all neutral higgses set to unity (done at construction).
595  }
596 
598  void fill_map_from_ScalarSingletDM_spectrum(std::map<std::string,double>& specmap, const Spectrum& singletdmspec,
599  const std::vector<SpectrumParameter>& required_parameters)
600  {
601  for(std::vector<SpectrumParameter>::const_iterator it = required_parameters.begin();
602  it != required_parameters.end(); ++it)
603  {
604  const Par::Tags tag = it->tag();
605  const std::string name = it->name();
606  const std::vector<int> shape = it->shape();
607 
609 
610  // Check scalar case
611  if(shape.size()==1 and shape[0]==1)
612  {
613  std::ostringstream label;
614  label << name <<" "<< Par::toString.at(tag);
615  specmap[label.str()] = singletdmspec.get_HE().get(tag,name);
616  }
617  // Check vector case
618  else if(shape.size()==1 and shape[0]>1)
619  {
620  for(int i = 1; i<=shape[0]; ++i) {
621  std::ostringstream label;
622  label << name <<"_"<<i<<" "<< Par::toString.at(tag);
623  specmap[label.str()] = singletdmspec.get_HE().get(tag,name,i);
624  }
625  }
626  // Check matrix case
627  else if(shape.size()==2)
628  {
629  for(int i = 1; i<=shape[0]; ++i) {
630  for(int j = 1; j<=shape[0]; ++j) {
631  std::ostringstream label;
632  label << name <<"_("<<i<<","<<j<<") "<<Par::toString.at(tag);
633  specmap[label.str()] = singletdmspec.get_HE().get(tag,name,i,j);
634  }
635  }
636  }
637  // Deal with all other cases
638  else
639  {
640  // ERROR
641  std::ostringstream errmsg;
642  errmsg << "Error, invalid parameter received while converting SingletDMspectrum to map of strings! This should no be possible if the spectrum content verification routines were working correctly; they must be buggy, please report this.";
643  errmsg << "Problematic parameter was: "<< tag <<", " << name << ", shape="<< shape;
644  utils_error().forced_throw(LOCAL_INFO,errmsg.str());
645  }
646  }
647 
648  }
649 
650  void get_ScalarSingletDM_Z2_spectrum_as_map (std::map<std::string,double>& specmap)
651  {
654  static const SpectrumContents::ScalarSingletDM_Z2 contents;
655  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters();
656  fill_map_from_ScalarSingletDM_spectrum(specmap, spec, required_parameters);
657  }
658 
659  void get_ScalarSingletDM_Z3_spectrum_as_map (std::map<std::string,double>& specmap)
660  {
663  static const SpectrumContents::ScalarSingletDM_Z3 contents;
664  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters();
665  fill_map_from_ScalarSingletDM_spectrum(specmap, spec, required_parameters);
666  }
667 
668  } // end namespace SpecBit
669 } // end namespace Gambit
670 
Define overloadings of the stream operator for various containers.
This class is used to deliver both information defined in the Standard Model (or potentially just QED...
Rollcall header for module SpecBit.
General small utility macros.
std::vector< str > invisibles
Particles that higgses can decay invisibly to.
Declarations of convenience (i.e.
EXPORT_SYMBOLS error & utils_error()
Utility errors.
void ScalarSingletDM_higgs_couplings_pwid(HiggsCouplingsTable &result)
Put together the Higgs couplings for the ScalarSingletDM models, from partial widths only...
double CP[max_neutral_higgses]
CP of neutral higgses.
A simple SubSpectrum wrapper for Standard Model Higgs information.
void get_ScalarSingletDM_Z3_spectrum_as_map(std::map< std::string, double > &specmap)
A simple SubSpectrum wrapper for the scalar singlet dark matter model.
#define LOCAL_INFO
Definition: local_info.hpp:34
GAMBIT native higgs coupling table class.
void fill_map_from_ScalarSingletDM_spectrum(std::map< std::string, double > &specmap, const Spectrum &singletdmspec, const std::vector< SpectrumParameter > &required_parameters)
Print ScalarSingletDM spectra out. Stripped down copy of MSSM version with variable names changed...
Flexiblesusy model header includes for SpecBit.
void get_ScalarSingletDM_Z3_spectrum(Spectrum &result)
Get a (simple) Spectrum object wrapper for the ScalarSingletDM_Z3 model.
ScalarSingletDM_Z3 derived version of SubSpectrum class.
START_MODEL b
Definition: demo.hpp:235
void get_ScalarSingletDM_Z2_spectrum_as_map(std::map< std::string, double > &specmap)
void fill_ScalarSingletDM_input(T &input, const std::map< str, safe_ptr< const double > > &Param, SMInputs sminputs)
ScalarSingletDM_Z2 derived version of SubSpectrum class.
std::chrono::milliseconds ms
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum ScalarSingletDM_Z3_spectrum
bool check_perturb(const Spectrum &spec, const std::vector< SpectrumParameter > &required_parameters, double scale, int pts)
void setup_QedQcd(softsusy::QedQcd &oneset, const SMInputs &sminputs)
Non-Gambit helper functions.
Definition: SpecBit.cpp:47
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
virtual double get(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
const double pi
An interface class for module dependencies.
void fill_extra_input(T &input, const std::map< str, safe_ptr< const double > > &Param)
std::vector< SpectrumParameter > all_parameters_with_tag(Par::Tags tag) const
Function to retreive all parameters matching a certain tag.
Header file that includes all GAMBIT headers required for a module source file.
Spectrum Spectrum
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
std::unique_ptr< SubSpectrum > clone_HE() const
Definition: spectrum.cpp:235
TYPE getValueOrDef(TYPE def, const args &... keys) const
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
This class is used to wrap the QedQcd object used by SoftSUSY and FlexibleSUSY in a Gambit SubSpectru...
void set_neutral_decays(int, const str &, const DecayTable::Entry &)
Assign decay entries to neutral higgses.
double pow(const double &a)
Outputs a^i.
std::vector< YAML::sdd > mc_info
Typedefs for making it easier to manipulate mass cut and mass ratio cut info.
Definition: spectrum.hpp:119
void get_ScalarSingletDM_Z2_spectrum(Spectrum &result)
Get a (simple) Spectrum object wrapper for the ScalarSingletDM_Z2 model.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum mh
Simple extension of the SMHiggsSimpleSpec "model object" to include scalar singlet DM parameters We c...
void set_neutral_decays_SM(int, const str &, const DecayTable::Entry &)
Assign decay entries to the various table components.
invalid_point_exception & invalid_point()
Invalid point exceptions.
Spectrum run_FS_spectrum_generator(const typename MI::InputParameters &input, const SMInputs &sminputs, const Options &runOptions, const std::map< str, safe_ptr< const double > > &input_Param)
Definition: SpecBit_MDM.cpp:55
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum SMINPUTS
std::vector< YAML::ssdd > mr_info
Definition: spectrum.hpp:120
TODO: see if we can use this one:
Definition: Analysis.hpp:33
A small wrapper object for &#39;options&#39; nodes.
SubSpectrum & get_HE()
Definition: spectrum.cpp:225
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110
Container class for Standard Model input information (defined as in SLHA2)
Definition: sminputs.hpp:29
std::vector< SpectrumParameter > all_parameters() const
Function to retreive all parameters.