gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
SpecBit_VS.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
20 
21 #include <string>
22 #include <sstream>
23 
24 #include <gsl/gsl_math.h>
25 #include <gsl/gsl_min.h>
26 
34 
35 #include "gambit/SpecBit/model_files_and_boxes.hpp" // #includes lots of flexiblesusy headers and defines interface classes
36 
37 // Flexible SUSY stuff (should not be needed by the rest of gambit)
38 #include "flexiblesusy/src/ew_input.hpp"
39 #include "flexiblesusy/src/lowe.h" // From softsusy; used by flexiblesusy
40 #include "flexiblesusy/src/numerics2.hpp"
41 
42 
43 // Switch for debug mode
44 //#define SPECBIT_DEBUG
45 
46 namespace Gambit
47 {
48 
49  namespace SpecBit
50  {
51 
52  using namespace LogTags;
53  using namespace flexiblesusy;
54 
56  {
57  // check that the electroweak scale stability conditions are satisfied
59 
61 
62  double lambda_h = fullspectrum.get(Par::dimensionless,"lambda_h");
63  double lambda_s = fullspectrum.get(Par::dimensionless,"lambda_S");
64  double lambda_hs = fullspectrum.get(Par::dimensionless,"lambda_hS");
65  double mu3 = fullspectrum.get(Par::mass1,"mu3");
66  double ms = fullspectrum.get(Par::Pole_Mass,"S");
67 
68  double check = 0;
69 
70  if ( (lambda_h*lambda_s > 0 ) )
71  {
72  check = 2 * pow( lambda_h * lambda_s , 0.5) + lambda_hs;
73  }
74  double check_2 = 2.*pow(abs(lambda_s),0.5)*ms - mu3;
75 
76  result = 0;
77 
78  // if any condition not satisfied set bad likelihood and invalidate point
79  if ( lambda_hs < 0 || lambda_s < 0 || check < 0 || check_2 < 0)
80  {
81  invalid_point().raise("Electroweak vacuum is unstable at low scale.");
82  result = -1e100;
83  }
84  }
85 
86  bool check_perturb_to_min_lambda(const Spectrum& spec,double scale,int pts
87  ,const std::vector<SpectrumParameter> required_parameters)
88  {
89 
90  std::unique_ptr<SubSpectrum> subspec = spec.clone_HE();
91  double step = log10(scale) / pts;
92  double runto;
93 
94  double ul=3.5449077018110318; // sqrt(4*Pi)
95  for (int i=0;i<pts;i++)
96  {
97  runto = pow(10,step*float(i+1.0)); // scale to run spectrum to
98  if (runto<100){runto=200.0;}// avoid running to low scales
99 
100  try
101  {
102  subspec -> RunToScale(runto);
103  }
104  catch (const Error& error)
105  {
106  return false;
107  }
108 
109  for(std::vector<SpectrumParameter>::const_iterator it = required_parameters.begin();
110  it != required_parameters.end(); ++it)
111  {
112  const Par::Tags tag = it->tag();
113  const std::string name = it->name();
114  const std::vector<int> shape = it->shape();
115  std::ostringstream label;
116  label << name <<" "<< Par::toString.at(tag);
117  if(shape.size()==1 and shape[0]==1)
118  {
119  if (abs(subspec->get(tag,name))>ul)
120  {
121  return false;
122  }
123  }
124 
125  else if(shape.size()==1 and shape[0]>1)
126  {
127  for(int k = 1; k<=shape[0]; ++k) {
128  if (abs(subspec->get(tag,name,k))>ul)
129  {
130  return false;
131  }
132 
133  }
134  }
135  else if(shape.size()==2)
136  {
137  for(int k = 1; k<=shape[0]; ++k) {
138  for(int j = 1; j<=shape[0]; ++j) {
139  if (abs(subspec->get(tag,name,k,j))>ul)
140  {
141  return false;
142  }
143  }
144  }
145  }
146  }
147  }
148 
149  return true;
150  }
151 
152  double run_lambda(double scale , void *params)
153  {
154 
155  std::unique_ptr<SubSpectrum>* spec=(std::unique_ptr<SubSpectrum>* )params;
156  SubSpectrum& speccloned = **spec;
157 
158  // clone the original spectrum incase the running takes the spectrum
159  // into a non-perturbative scale and thus the spectrum is no longer reliable
160  std::unique_ptr<SubSpectrum> speccloned2 = speccloned.clone();
161 
162  if (scale>1.0e21){scale=1.0e21;}// avoid running to high scales
163 
164  if (scale<100.0){scale=100.0;}// avoid running to very low scales
165 
166  double lambda;
167  try
168  {
169  speccloned2->RunToScale(scale);
170  lambda = speccloned2->get(Par::dimensionless,"lambda_h");
171  }
172  catch (const Error& error)
173  {
174  // ERROR(error.what());
175  //cout << "error encountered" << endl;
176  lambda = 0;
177  // return EXIT_FAILURE;
178  }
179 
180  return lambda;
181  }
182 
183 
184  void find_min_lambda_Helper(dbl_dbl_bool& vs_tuple, const Spectrum& fullspectrum,
185  double high_energy_limit, int check_perturb_pts,
186  const std::vector<SpectrumParameter> required_parameters)
187  {
188  std::unique_ptr<SubSpectrum> speccloned = fullspectrum.clone_HE();
189 
190  // three scales at which we choose to run the quartic coupling up to, and then use a Lagrange interpolating polynomial
191  // to get an estimate for the location of the minimum, this is an efficient way to narrow down over a huge energy range
192  double u_1 = 1, u_2 = 5, u_3 = 12;
193  double lambda_1,lambda_2,lambda_3;
194  double lambda_min = 0;
195 
196 
197  bool min_exists = 1;// check if gradient is positive at electroweak scale
198  if ( run_lambda(101.0, &speccloned ) > run_lambda(100.0, &speccloned ) )
199  {
200  // gradient is positive, the minimum is less than electroweak scale so
201  // lambda_h must be monotonally increasing
202  min_exists = 0;
203  lambda_min = run_lambda(100.0,&speccloned);
204  }
205 
206  double mu_min = 0;
207  if (min_exists)
208  {
209 
210  // fit parabola (in log space) to the 3 trial points and use this to estimate the minimum
211 
212  for (int i=1;i<2;i++)
213  {
214 
215  lambda_1 = run_lambda(pow(10,u_1),&speccloned);
216  lambda_2 = run_lambda(pow(10,u_2),&speccloned);
217  lambda_3 = run_lambda(pow(10,u_3),&speccloned);
218 
219  double min_u= (lambda_1*(pow(u_2,2)-pow(u_3,2)) - lambda_2*(pow(u_1,2)-pow(u_3,2)) + lambda_3*(pow(u_1,2)-pow(u_2,2)));
220  double denominator = ( lambda_1*(u_2-u_3)+ lambda_2*(u_3-u_1) +lambda_3*(u_1-u_2));
221 
222  min_u=0.5*(min_u/denominator);
223  u_1=min_u-2/(pow(float(i),0.01));
224  u_2=min_u;
225  u_3=min_u+2/(pow(float(i),0.01));
226 
227  }
228 
229  // run downhill minimization routine to find exact minimum
230 
231  double mu_lower = pow(10,u_1);
232  double mu_upper = pow(10,u_3);
233  mu_min = pow(10,u_2);
234 
235  gsl_function F;
236  F.function = &run_lambda;
237  F.params = &speccloned;
238 
239  int status;
240  int iteration = 0, max_iteration = 1000;
241 
242  const gsl_min_fminimizer_type *T;
243  gsl_min_fminimizer *s;
244 
245  T = gsl_min_fminimizer_brent;
246  s = gsl_min_fminimizer_alloc (T);
247  gsl_min_fminimizer_set (s, &F, mu_min, mu_lower, mu_upper);
248 
249  do
250  {
251  iteration++;
252  status = gsl_min_fminimizer_iterate (s);
253 
254  mu_min = gsl_min_fminimizer_x_minimum (s);
255  mu_lower = gsl_min_fminimizer_x_lower (s);
256  mu_upper = gsl_min_fminimizer_x_upper (s);
257 
258  status = gsl_min_test_interval (mu_lower, mu_upper, 0.0001, 0.0001);
259  // cout << "mu_lower = " << mu_lower << " mu_upper = " << mu_upper << endl;
260  }
261  while (status == GSL_CONTINUE && iteration < max_iteration);
262 
263  if (iteration == max_iteration)
264  {
265  SpecBit_error().raise(LOCAL_INFO,"The minimum of the quartic coupling could not be found");
266  }
267 
268  gsl_min_fminimizer_free (s);
269 
270  lambda_min = run_lambda(mu_min,&speccloned);
271 
272  }
273 
274  #ifdef SPECBIT_DEBUG
275  std::cout<< "minimum value of quartic coupling is "<< lambda_min << " at " << mu_min <<" GeV"<<std::endl;
276  #endif
277 
278  double lifetime,LB;
279  if (lambda_min<0) // second minimum exists, now determine stability and lifetime
280  {
281  LB=mu_min;
282  #ifdef SPECBIT_DEBUG
283  double p=exp(4*140-26/(abs(0.5*lambda_min)))*pow(LB/(1.2e19),4); // compute tunnelling rate
284  if (p>1)
285  {
286  cout<< "vacuum is unstable" << endl;
287  }
288  else
289  {
290  cout<< "vacuum is metastable" << endl;
291  }
292  #endif
293  double pi2_8_over3 = 8.* pow ( pi , 2 ) / 3.;
294  double conversion = (6.5821195e-25)/(31536000);
295  // top factor is hbar in units of GeV.s and bottom factor is number of seconds in a year
296  lifetime=conversion/(exp(3*140-pi2_8_over3/(abs(0.5*lambda_min)))*pow(1/(1.2e19),3)*pow(LB,4));
297 
298  }
299  else // quartic coupling always positive, set output to default values
300  {
301  LB=high_energy_limit;
302  lifetime=1e300;
303  #ifdef SPECBIT_DEBUG
304  cout << "vacuum is absolutely stable" << endl;
305  #endif
306  // vacuum is stable
307  }
308  // now do a check on the perturbativity of the couplings up to this scale
309  bool perturbative=check_perturb_to_min_lambda(fullspectrum,LB,check_perturb_pts,required_parameters);
310  double perturb=float(!perturbative);
311  #ifdef SPECBIT_DEBUG
312  cout << "perturbativity checked up to " << LB << " result = " << perturbative << endl;
313  cout << "Higgs pole mass = " << fullspectrum.get(Par::Pole_Mass, "h0_1") << endl;
314  #endif
315  vs_tuple = dbl_dbl_bool(lifetime,LB,perturb);
316  }
317 
318 
320  {
322  double high_energy_limit = myPipe::runOptions->getValueOrDef<double>(1.22e19,"set_high_scale");
323  int check_perturb_pts = myPipe::runOptions->getValueOrDef<double>(10,"check_perturb_pts");
324  static const SpectrumContents::ScalarSingletDM_Z2 contents;
325  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
326  find_min_lambda_Helper(vs_tuple, *myPipe::Dep::ScalarSingletDM_Z2_spectrum, high_energy_limit, check_perturb_pts, required_parameters);
327  }
328 
330  {
332  double high_energy_limit = myPipe::runOptions->getValueOrDef<double>(1.22e19,"set_high_scale");
333  int check_perturb_pts = myPipe::runOptions->getValueOrDef<double>(10,"check_perturb_pts");
334  static const SpectrumContents::ScalarSingletDM_Z2 contents;
335  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
336  find_min_lambda_Helper(vs_tuple, *myPipe::Dep::ScalarSingletDM_Z3_spectrum, high_energy_limit, check_perturb_pts, required_parameters);
337  }
338 
340  {
341  namespace myPipe = Pipes::find_min_lambda_MDM;
342  double high_energy_limit = myPipe::runOptions->getValueOrDef<double>(1.22e19,"set_high_scale");
343  int check_perturb_pts = myPipe::runOptions->getValueOrDef<double>(10,"check_perturb_pts");
344  static const SpectrumContents::MDM contents;
345  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters_with_tag(Par::dimensionless);
346  find_min_lambda_Helper(vs_tuple, *myPipe::Dep::MDM_spectrum, high_energy_limit, check_perturb_pts, required_parameters);
347  }
348 
349 
350  // the functions below are used to extract the desired outputs from find_min_lambda
351 
352  // gives expected lifetime in units of years, if stable give extremly large number (1e300)
353  void get_expected_vacuum_lifetime(double &lifetime)
354  {
355  namespace myPipe = Pipes::get_expected_vacuum_lifetime;
356  dbl_dbl_bool vs_tuple = *myPipe::Dep::high_scale_vacuum_info;
357 
358  if (vs_tuple.first<1e300)
359  {
360  lifetime=vs_tuple.first;
361  }
362  else
363  {
364  lifetime=1e300;
365  }
366  }
367 
368  // log of the likelihood
370  {
372  dbl_dbl_bool vs_tuple = *myPipe::Dep::high_scale_vacuum_info;
373 
374  const Options& runOptions=*myPipe::runOptions;
375  bool demand_stable = runOptions.getValueOrDef<bool>(false,"demand_stable");
376  double stability_scale = runOptions.getValueOrDef<double>(1.22e19,"set_stability_scale");
377 
378  if (demand_stable && (vs_tuple.second < stability_scale))
379  {
380  result = -1e100;
381  }
382  else
383  {
384  double conversion = (6.5821195e-25)/(31536000);
385  result=((- ( 1 / ( vs_tuple.first/conversion ) ) * exp(140) * (1/ (1.2e19) ) ) );
386  }
387 
388  }
389 
390  // Get the scale of the high-scale minimum
391  void get_lambdaB(double &result)
392  {
393  namespace myPipe = Pipes::get_lambdaB;
394  dbl_dbl_bool vs_tuple = *myPipe::Dep::high_scale_vacuum_info;
395  result=vs_tuple.second;
396  }
397 
398 
399  // Returns poor likelihood and invalidates point if couplings go non-perturbative at or below the scale of the high-scale minimum of the potential
400  void check_perturb_min_lambda(double &result)
401  {
402  namespace myPipe = Pipes::check_perturb_min_lambda;
403  dbl_dbl_bool vs_tuple = *myPipe::Dep::high_scale_vacuum_info;
404 
405  if (vs_tuple.flag)
406  {
407  invalid_point().raise("Couplings are non-perturbative before scale of vacuum instability");
408  result = -1e100;
409  }
410  else
411  {
412  result = 0;
413  }
414  }
415 
416 
417  } // end namespace SpecBit
418 } // end namespace Gambit
419 
double get(const Par::Tags partype, const std::string &mass) const
Definition: spectrum.cpp:249
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.
virtual std::unique_ptr< SubSpectrum > clone() const =0
Clone the SubSpectrum object.
Declarations of convenience (i.e.
bool check_perturb_to_min_lambda(const Spectrum &spec, double scale, int pts, const std::vector< SpectrumParameter > required_parameters)
Definition: SpecBit_VS.cpp:86
#define LOCAL_INFO
Definition: local_info.hpp:34
Flexiblesusy model header includes for SpecBit.
void find_min_lambda_ScalarSingletDM_Z3(dbl_dbl_bool &vs_tuple)
Definition: SpecBit_VS.cpp:329
void find_min_lambda_Helper(dbl_dbl_bool &vs_tuple, const Spectrum &fullspectrum, double high_energy_limit, int check_perturb_pts, const std::vector< SpectrumParameter > required_parameters)
Definition: SpecBit_VS.cpp:184
std::chrono::milliseconds ms
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum ScalarSingletDM_Z3_spectrum
void get_lambdaB(double &result)
Definition: SpecBit_VS.cpp:391
GAMBIT error class.
Definition: exceptions.hpp:136
virtual void raise(const std::string &)
Raise the exception, i.e. throw it.
Definition: exceptions.cpp:422
const double pi
void check_perturb_min_lambda(double &result)
Definition: SpecBit_VS.cpp:400
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.
std::unique_ptr< SubSpectrum > clone_HE() const
Definition: spectrum.cpp:235
TYPE getValueOrDef(TYPE def, const args &... keys) const
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
This class is used to wrap the QedQcd object used by SoftSUSY and FlexibleSUSY in a Gambit SubSpectru...
void find_min_lambda_ScalarSingletDM_Z2(dbl_dbl_bool &vs_tuple)
Definition: SpecBit_VS.cpp:319
void find_min_lambda_MDM(dbl_dbl_bool &vs_tuple)
Definition: SpecBit_VS.cpp:339
double pow(const double &a)
Outputs a^i.
invalid_point_exception & invalid_point()
Invalid point exceptions.
double run_lambda(double scale, void *params)
Definition: SpecBit_VS.cpp:152
void check_EW_stability_ScalarSingletDM_Z3(double &result)
Definition: SpecBit_VS.cpp:55
TODO: see if we can use this one:
Definition: Analysis.hpp:33
A small wrapper object for &#39;options&#39; nodes.
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110
void get_expected_vacuum_lifetime(double &lifetime)
Definition: SpecBit_VS.cpp:353
void lnL_highscale_vacuum_decay_single_field(double &result)
Definition: SpecBit_VS.cpp:369