gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
SpecBit_MDM.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
18 
19 #include <string>
20 #include <sstream>
21 
23 
27 
33 
35 
36 // Flexible SUSY stuff (should not be needed by the rest of gambit)
37 #include "flexiblesusy/src/ew_input.hpp"
38 #include "flexiblesusy/src/lowe.h"
39 #include "flexiblesusy/src/numerics2.hpp"
40 #include "flexiblesusy/src/spectrum_generator_settings.hpp"
41 // Switch for debug mode
42 //#define SPECBIT_DEBUG
43 
44 namespace Gambit
45 {
46 
47  namespace SpecBit
48  {
49  using namespace LogTags;
50  using namespace flexiblesusy;
51 
52 
53  template<class MI,class SI>
55  ( const typename MI::InputParameters& input
56  , const SMInputs& sminputs
57  , const Options& runOptions
58  , const std::map<str, safe_ptr<const double> >& input_Param
59  )
60  {
61  softsusy::QedQcd oneset;
62 
63  // Fill QedQcd object with SMInputs values
64  setup_QedQcd(oneset,sminputs);
65 
66  // Run everything to Mz
67  oneset.toMz();
68 
69  // Create spectrum generator object
70  typename MI::SpectrumGenerator spectrum_generator;
71 
72  Spectrum_generator_settings settings;
73  settings.set(Spectrum_generator_settings::precision, runOptions.getValueOrDef<double>(1.0e-4,"precision_goal"));
74  settings.set(Spectrum_generator_settings::max_iterations, runOptions.getValueOrDef<double>(0,"max_iterations"));
75  settings.set(Spectrum_generator_settings::calculate_sm_masses, runOptions.getValueOrDef<bool> (true, "calculate_sm_masses"));
76  settings.set(Spectrum_generator_settings::pole_mass_loop_order, runOptions.getValueOrDef<int>(2,"pole_mass_loop_order"));
77  settings.set(Spectrum_generator_settings::pole_mass_loop_order, runOptions.getValueOrDef<int>(2,"ewsb_loop_order"));
78  settings.set(Spectrum_generator_settings::beta_loop_order, runOptions.getValueOrDef<int>(2,"beta_loop_order"));
79  settings.set(Spectrum_generator_settings::threshold_corrections_loop_order, runOptions.getValueOrDef<int>(2,"threshold_corrections_loop_order"));
80  settings.set(Spectrum_generator_settings::higgs_2loop_correction_at_as, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_at_as"));
81  settings.set(Spectrum_generator_settings::higgs_2loop_correction_ab_as, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_ab_as"));
82  settings.set(Spectrum_generator_settings::higgs_2loop_correction_at_at, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_at_at"));
83  settings.set(Spectrum_generator_settings::higgs_2loop_correction_atau_atau, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_atau_atau"));
84 
85  spectrum_generator.set_settings(settings);
86 
87  // Generate spectrum
88  spectrum_generator.run(oneset, input);
89  const typename MI::Problems& problems = spectrum_generator.get_problems();
90  MI model_interface(spectrum_generator,oneset,input);
91  SI mdmspec(model_interface, "FlexibleSUSY", "1.5.1");
92 
93 
94  mdmspec.set_override(Par::mass1,spectrum_generator.get_high_scale(),"high_scale",true);
95  mdmspec.set_override(Par::mass1,spectrum_generator.get_susy_scale(),"susy_scale",true);
96  mdmspec.set_override(Par::mass1,spectrum_generator.get_low_scale(), "low_scale", true);
97 
98  mdmspec.set_override(Par::Pole_Mass_1srd_high,0.0, "h0_1", true);
99  mdmspec.set_override(Par::Pole_Mass_1srd_low,0.0,"h0_1", true);
100 
101  QedQcdWrapper qedqcdspec(oneset,sminputs);
102 
103  // Deal with points where spectrum generator encountered a problem
104  #ifdef SPECBIT_DEBUG
105  std::cout<<"Problem? "<<problems.have_problem()<<std::endl;
106  std::cout<<"Warning? "<<problems.have_warning()<<std::endl;
107  #endif
108  if( problems.have_problem() )
109  {
110  if( runOptions.getValueOrDef<bool>(false,"invalid_point_fatal") )
111  {
112  std::ostringstream errmsg;
113  errmsg << "A serious problem was encountered during spectrum generation!";
114  errmsg << "Message from FlexibleSUSY:" << std::endl;
115  problems.print_problems(errmsg);
116  problems.print_warnings(errmsg);
117  SpecBit_error().raise(LOCAL_INFO,errmsg.str());
118  }
119  else
120  {
121  // Check what the problem was
122  // see: contrib/MassSpectra/flexiblesusy/src/problems.hpp
123  std::ostringstream msg;
124  problems.print_problems(msg);
125  invalid_point().raise(msg.str()); //TODO: This message isn't ending up in the logs.
126  }
127  }
128  double QEWSB = *input_Param.at("QEWSB");
129  mdmspec.RunToScaleOverride(QEWSB);
130 
131  // Retrieve any mass cuts
132  static const Spectrum::mc_info mass_cut = runOptions.getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
133  static const Spectrum::mr_info mass_ratio_cut = runOptions.getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
134 
135  return Spectrum(qedqcdspec,mdmspec,sminputs,&input_Param, mass_cut, mass_ratio_cut);
136 
137  }
138 
139 
140  template <class T>
141  void fill_MDM_input(T& input, const std::map<str, safe_ptr<const double> >& Param,SMInputs sminputs)
142  {
143  double mH = *Param.at("mH");
144  double mChi = *Param.at("mChi");
145 
146  //double QEFT = mChi;
147  /*
148  if (QEFT < sminputs.mT)
149  {
150  QEFT = sminputs.mT;
151  }
152  */
153  double QEFT = sminputs.mT;
154 
155  input.HiggsIN = 0.5*pow(mH,2);
156 
157  input.YcIN = 0.5*mChi;
158 
159  input.QEWSB = QEFT; // scale where EWSB conditions are applied
160  input.Qin = QEFT; //scale; // highest scale at which model is run to
161 
162  }
163 
164  bool check_perturb_MDM(const Spectrum& spec,double scale,int pts)
165  {
166  using namespace flexiblesusy;
167  using namespace Gambit;
168  using namespace SpecBit;
169  std::unique_ptr<SubSpectrum> MDM = spec.clone_HE();
170  double step = log10(scale) / pts;
171  double runto;
172 
173  //const double ul = std::sqrt(4.0 * pi); // Maximum value for perturbative couplings, same perturbativity bound that FlexibleSUSY uses
174  double ul = 4.0 * pi;
175  for (int i=0;i<pts;i++)
176  {
177  runto = pow(10,step*float(i+1.0)); // scale to run spectrum to
178  if (runto<100){runto=100.0;}// avoid running to low scales
179 
180  try
181  {
182  MDM -> RunToScale(runto);
183  }
184  catch (const Error& error)
185  {
186  return false;
187  };
188 
189 
190 
191 
192  static const SpectrumContents::MDM contents;
193  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
194 
195  for(std::vector<SpectrumParameter>::const_iterator it = required_parameters.begin();
196  it != required_parameters.end(); ++it)
197  {
198  const Par::Tags tag = it->tag();
199  const std::string name = it->name();
200  const std::vector<int> shape = it->shape();
201  std::ostringstream label;
202  label << name <<" "<< Par::toString.at(tag);
203 
204  if (name == "lambda_h"){ul = 2*pi;}
205  else {ul = 4.0 * pi;}
206 
207  if(shape.size()==1 and shape[0]==1)
208  {
209  if (abs(MDM->get(tag,name))>ul)
210  {
211  return false;
212  }
213  }
214  else if(shape.size()==1 and shape[0]>1)
215  {
216  for(int k = 1; k<=shape[0]; ++k)
217  {
218  if (abs(MDM->get(tag,name,k))>ul) return false;
219  }
220  }
221  else if(shape.size()==2)
222  {
223  for(int k = 1; k<=shape[0]; ++k)
224  {
225  for(int j = 1; j<=shape[0]; ++j)
226  {
227  if (abs(MDM->get(tag,name,k,j))>ul) return false;
228  }
229  }
230  }
231 
232  }
233  }
234 
235  return true;
236 
237  }
238 
239  #if(FS_MODEL_MDM_IS_BUILT)
240  void get_MDM_spectrum(Spectrum& result)
241  {
242  using namespace softsusy;
243  namespace myPipe = Pipes::get_MDM_spectrum;
244  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
245  const Options& runOptions=*myPipe::runOptions;
246  //double scale = runOptions.getValueOrDef<double>(173.34,"FS_high_scale");
247  MDM_input_parameters input;
248  fill_MDM_input(input,myPipe::Param,sminputs);
249  result = run_FS_spectrum_generator<MDM_interface<ALGORITHM1>,MDMSpec<MDM_interface<ALGORITHM1>>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
250 
251  int check_perturb_pts = runOptions.getValueOrDef<double>(10,"check_perturb_pts");
252  double do_check_perturb = runOptions.getValueOrDef<bool>(false,"check_perturb");
253  double check_perturb_scale = runOptions.getValueOrDef<double>(1.22e19,"check_high_scale");
254  double input_scale_tolerance = runOptions.getValueOrDef<double>(1e-3,"input_scale_tolerance");
255 
256  // Check that the Higgs MSbar parameter has been provided at the scale of the MDM MSbar multiplet mass
257  double Qin = *myPipe::Param.at("Qin");
258  if (abs((Qin - *myPipe::Param.at("mChi"))/Qin) > input_scale_tolerance) SpecBit_error().raise(LOCAL_INFO, "SM_Higgs_running::Qin must equal MDM::mChi");
259 
260  if (do_check_perturb)
261  {
262  if (!check_perturb_MDM(result,check_perturb_scale,check_perturb_pts))
263  {
264  // invalidate point as spectrum not perturbative up to scale
265  std::ostringstream msg;
266  msg << "Spectrum not perturbative up to scale = " << check_perturb_scale << std::endl;
267  #ifdef SPECBIT_DEBUG
268  cout << "Spectrum not perturbative up to scale = " << check_perturb_scale << endl;
269  #endif
270  invalid_point().raise(msg.str());
271  }
272  }
273 
274 
275  }
276  #endif
277 
278 
279  void find_non_perturb_scale_MDM(double &result)
280  {
281  using namespace flexiblesusy;
282  using namespace softsusy;
283  namespace myPipe = Pipes::find_non_perturb_scale_MDM;
284  using namespace Gambit;
285  using namespace SpecBit;
286 
287  const Spectrum& fullspectrum = *myPipe::Dep::MDM_spectrum;
288 
289  // bound x by (a,b)
290  // do all this is log space please
291 
292  double ms = *myPipe::Param.at("mS");
293 
294  double a = log10(ms);
295 
296  if (a > 20.0)
297  {
298  std::ostringstream msg;
299  msg << "Scalar mass larger than 10^20 GeV " << std::endl;
300  invalid_point().raise(msg.str());
301  }
302 
303  double b = 20.0;
304  double x = 0.5 * ( b + ms );
305 
306  while (abs(a-b)>1e-10)
307  {
308  //cout<< "\r" << "(a,b) = " << a << " " << b << endl;
309  //std::cout << std::flush;
310 
311  x=0.5*(b-a)+a;
312 
313 
314  if (!check_perturb_MDM(fullspectrum,pow(10,x),3))
315  {
316  b=x;
317  }
318  else
319  {
320  a=x;
321  }
322  }
323  result = pow(10,0.5*(a+b));
324  }
325 
326 
328  void fill_map_from_MDMspectrum(std::map<std::string,double>&, const Spectrum&);
329 
330  void get_MDM_spectrum_as_map (std::map<std::string,double>& specmap)
331  {
332  namespace myPipe = Pipes::get_MDM_spectrum_as_map;
333  const Spectrum& mdmspec(*myPipe::Dep::MDM_spectrum);
334  fill_map_from_MDMspectrum(specmap, mdmspec);
335  }
336 
337  void fill_map_from_MDMspectrum(std::map<std::string,double>& specmap, const Spectrum& mdmspec)
338  {
340  static const SpectrumContents::MDM contents;
341  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters();
342 
343  for(std::vector<SpectrumParameter>::const_iterator it = required_parameters.begin();
344  it != required_parameters.end(); ++it)
345  {
346  const Par::Tags tag = it->tag();
347  const std::string name = it->name();
348  const std::vector<int> shape = it->shape();
349 
351 
352  // Check scalar case
353  if(shape.size()==1 and shape[0]==1)
354  {
355  std::ostringstream label;
356  label << name <<" "<< Par::toString.at(tag);
357  specmap[label.str()] = mdmspec.get_HE().get(tag,name);
358  }
359  // Check vector case
360  else if(shape.size()==1 and shape[0]>1)
361  {
362  for(int i = 1; i<=shape[0]; ++i) {
363  std::ostringstream label;
364  label << name <<"_"<<i<<" "<< Par::toString.at(tag);
365  specmap[label.str()] = mdmspec.get_HE().get(tag,name,i);
366  }
367  }
368  // Check matrix case
369  else if(shape.size()==2)
370  {
371  for(int i = 1; i<=shape[0]; ++i) {
372  for(int j = 1; j<=shape[0]; ++j) {
373  std::ostringstream label;
374  label << name <<"_("<<i<<","<<j<<") "<<Par::toString.at(tag);
375  specmap[label.str()] = mdmspec.get_HE().get(tag,name,i,j);
376  }
377  }
378  }
379  // Deal with all other cases
380  else
381  {
382  // ERROR
383  std::ostringstream errmsg;
384  errmsg << "Error, invalid parameter received while converting MDMspectrum 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.";
385  errmsg << "Problematic parameter was: "<< tag <<", " << name << ", shape="<< shape;
386  utils_error().forced_throw(LOCAL_INFO,errmsg.str());
387  }
388  }
389 
390  }
391 
392  } // end namespace SpecBit
393 } // end namespace Gambit
394 
MDM derived version of SubSpectrum class.
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.
Declarations of convenience (i.e.
EXPORT_SYMBOLS error & utils_error()
Utility errors.
void find_non_perturb_scale_MDM(double &result)
bool check_perturb_MDM(const Spectrum &spec, double scale, int pts)
A simple SubSpectrum wrapper for Standard Model Higgs information.
#define LOCAL_INFO
Definition: local_info.hpp:34
Flexiblesusy model header includes for SpecBit.
START_MODEL b
Definition: demo.hpp:235
std::chrono::milliseconds ms
GAMBIT error class.
Definition: exceptions.hpp:136
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
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
void fill_MDM_input(T &input, const std::map< str, safe_ptr< const double > > &Param, SMInputs sminputs)
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
std::unique_ptr< SubSpectrum > clone_HE() const
Definition: spectrum.cpp:235
void fill_map_from_MDMspectrum(std::map< std::string, double > &, const Spectrum &)
Print MDM spectrum out. Stripped down copy from MSSM version with variable names changed.
TYPE getValueOrDef(TYPE def, const args &... keys) const
This class is used to wrap the QedQcd object used by SoftSUSY and FlexibleSUSY in a Gambit SubSpectru...
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
invalid_point_exception & invalid_point()
Invalid point exceptions.
void get_MDM_spectrum_as_map(std::map< std::string, double > &specmap)
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.