gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
SpecBit_MSSM.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
32 
33 #include <string>
34 #include <sstream>
35 #include <cmath>
36 #include <complex>
37 
42 #include "gambit/Utils/stream_overloads.hpp" // Just for more convenient output to logger
48 #include "gambit/SpecBit/model_files_and_boxes.hpp" // #includes lots of flexiblesusy headers and defines interface classes
49 
50 // Flexible SUSY stuff (should not be needed by the rest of gambit)
51 #include "flexiblesusy/src/ew_input.hpp"
52 #include "flexiblesusy/src/lowe.h" // From softsusy; used by flexiblesusy
53 #include "flexiblesusy/src/numerics2.hpp"
54 //#include "flexiblesusy/src/mssm_twoloophiggs.hpp"
55 #include "flexiblesusy/src/spectrum_generator_settings.hpp"
56 
57 // Switch for debug mode
58 //#define SPECBIT_DEBUG
59 
60 namespace Gambit
61 {
62 
63  namespace SpecBit
64  {
65  using namespace LogTags;
66  using namespace flexiblesusy;
67 
68  // To check if a model is currently being scanned:
69  // bool Pipes::<fname>::ModelInUse(str model_name)
70 
72  // =======================================================================
73  // These are not known to Gambit, but they do basically all the real work.
74  // The Gambit module functions merely wrap the functions here and hook
75  // them up to their dependencies, and input parameters.
76 
78  // In GAMBIT there are THREE flexiblesusy MSSM spectrum generators currently in
79  // use, for each of three possible boundary condition types:
80  // - GUT scale input
81  // - Electroweak symmetry breaking scale input
82  // - Intermediate scale Q input
83  // These each require slightly different setup, but once that is done the rest
84  // of the code required to run them is the same; this is what is contained in
85  // the below template function.
86  // MI for Model Interface, as defined in model_files_and_boxes.hpp
87  template <class MI>
89  ( const typename MI::InputParameters& input
90  , const SMInputs& sminputs
91  , const Options& runOptions
92  , const std::map<str, safe_ptr<const double> >& input_Param
93  )
94  {
95  // SoftSUSY object used to set quark and lepton masses and gauge
96  // couplings in QEDxQCD effective theory
97  // Will be initialised by default using values in lowe.h, which we will
98  // override next.
99  softsusy::QedQcd oneset;
100 
101  // Fill QedQcd object with SMInputs values
102  setup_QedQcd(oneset,sminputs);
103 
104  // Run everything to Mz
105  //oneset.toMz();
106 
107  // Create spectrum generator object
108  typename MI::SpectrumGenerator spectrum_generator;
109 
110  // Spectrum generator settings
111  // Default options copied from flexiblesusy/src/spectrum_generator_settings.hpp
112  //
113  // | enum | possible values | default value |
114  // |----------------------------------|------------------------------|-----------------|
115  // | precision | any positive double | 1.0e-4 |
116  // | max_iterations | any positive double | 0 (= automatic) |
117  // | algorithm | 0 (two-scale) or 1 (lattice) | 0 (= two-scale) |
118  // | calculate_sm_masses | 0 (no) or 1 (yes) | 0 (= no) |
119  // | pole_mass_loop_order | 0, 1, 2 | 2 (= 2-loop) |
120  // | ewsb_loop_order | 0, 1, 2 | 2 (= 2-loop) |
121  // | beta_loop_order | 0, 1, 2 | 2 (= 2-loop) |
122  // | threshold_corrections_loop_order | 0, 1 | 1 (= 1-loop) |
123  // | higgs_2loop_correction_at_as | 0, 1 | 1 (= enabled) |
124  // | higgs_2loop_correction_ab_as | 0, 1 | 1 (= enabled) |
125  // | higgs_2loop_correction_at_at | 0, 1 | 1 (= enabled) |
126  // | higgs_2loop_correction_atau_atau | 0, 1 | 1 (= enabled) |
127 
128  Spectrum_generator_settings settings;
129  settings.set(Spectrum_generator_settings::precision, runOptions.getValueOrDef<double>(1.0e-4,"precision_goal"));
130  settings.set(Spectrum_generator_settings::max_iterations, runOptions.getValueOrDef<double>(0,"max_iterations"));
131  settings.set(Spectrum_generator_settings::calculate_sm_masses, runOptions.getValueOrDef<bool> (false, "calculate_sm_masses"));
132  settings.set(Spectrum_generator_settings::pole_mass_loop_order, runOptions.getValueOrDef<int>(2,"pole_mass_loop_order"));
133  settings.set(Spectrum_generator_settings::pole_mass_loop_order, runOptions.getValueOrDef<int>(2,"ewsb_loop_order"));
134  settings.set(Spectrum_generator_settings::beta_loop_order, runOptions.getValueOrDef<int>(2,"beta_loop_order"));
135  settings.set(Spectrum_generator_settings::threshold_corrections_loop_order, runOptions.getValueOrDef<int>(2,"threshold_corrections_loop_order"));
136  settings.set(Spectrum_generator_settings::higgs_2loop_correction_at_as, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_at_as"));
137  settings.set(Spectrum_generator_settings::higgs_2loop_correction_ab_as, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_ab_as"));
138  settings.set(Spectrum_generator_settings::higgs_2loop_correction_at_at, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_at_at"));
139  settings.set(Spectrum_generator_settings::higgs_2loop_correction_atau_atau, runOptions.getValueOrDef<int>(1,"higgs_2loop_correction_atau_atau"));
140  settings.set(Spectrum_generator_settings::top_pole_qcd_corrections, runOptions.getValueOrDef<int>(1,"top_pole_qcd_corrections"));
141  settings.set(Spectrum_generator_settings::beta_zero_threshold, runOptions.getValueOrDef<double>(1.000000000e-14,"beta_zero_threshold"));
142  settings.set(Spectrum_generator_settings::eft_matching_loop_order_up, runOptions.getValueOrDef<int>(1,"eft_matching_loop_order_up"));
143  settings.set(Spectrum_generator_settings::eft_matching_loop_order_down, runOptions.getValueOrDef<int>(1,"eft_matching_loop_order_down"));
144  settings.set(Spectrum_generator_settings::threshold_corrections, runOptions.getValueOrDef<int>(123111321,"threshold_corrections"));
145 
146 
147  spectrum_generator.set_settings(settings);
148 
149  // Generate spectrum
150  spectrum_generator.run(oneset, input);
151 
152  // Extract report on problems...
153  const typename MI::Problems& problems = spectrum_generator.get_problems();
154 
155  // Create Model_interface to carry the input and results, and know
156  // how to access the flexiblesusy routines.
157  // Note: Output of spectrum_generator.get_model() returns type, e.g. CMSSM.
158  // Need to convert it to type CMSSM_slha (which alters some conventions of
159  // parameters into SLHA format)
160  MI model_interface(spectrum_generator,oneset,input);
161 
162  // Create SubSpectrum object to wrap flexiblesusy data
163  // THIS IS STATIC so that it lives on once we leave this module function. We
164  // therefore cannot run the same spectrum generator twice in the same loop and
165  // maintain the spectrum resulting from both. But we should never want to do
166  // this.
167  // A pointer to this object is what gets turned into a SubSpectrum pointer and
168  // passed around Gambit.
169  //
170  // This object will COPY the interface data members into itself, so it is now the
171  // one-stop-shop for all spectrum information, including the model interface object.
172  MSSMSpec<MI> mssmspec(model_interface, "FlexibleSUSY", "2.0.beta");
173 
174  // Add extra information about the scales used to the wrapper object
175  // (last parameter turns on the 'allow_new' option for the override setter, which allows
176  // us to set parameters that don't previously exist)
177  mssmspec.set_override(Par::mass1,spectrum_generator.get_high_scale(),"high_scale",true);
178  mssmspec.set_override(Par::mass1,spectrum_generator.get_susy_scale(),"susy_scale",true);
179  mssmspec.set_override(Par::mass1,spectrum_generator.get_low_scale(), "low_scale", true);
180 
181 
182  // Has the user chosen to override any pole mass values?
183  // This will typically break consistency, but may be useful in some special cases
184  if (runOptions.hasKey("override_FS_pole_masses"))
185  {
186  std::vector<str> particle_names = runOptions.getNames("override_FS_pole_masses");
187  for (auto& name : particle_names)
188  {
189  double mass = runOptions.getValue<double>("override_FS_pole_masses", name);
190  mssmspec.set_override(Par::Pole_Mass, mass, name);
191  }
192  }
193 
194  // Add theory errors
195  static const MSSM_strs ms;
196 
197  static const std::vector<int> i12 = initVector(1,2);
198  static const std::vector<int> i123 = initVector(1,2,3);
199  static const std::vector<int> i1234 = initVector(1,2,3,4);
200  static const std::vector<int> i123456 = initVector(1,2,3,4,5,6);
201 
202  // 3% theory "error"
205  mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, ms.pole_mass_strs_1_6, i123456, true);
206  mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, ms.pole_mass_strs_1_6, i123456, true);
207  mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, "~chi0", i1234, true);
208  mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, "~chi0", i1234, true);
209  mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, ms.pole_mass_strs_1_3, i123, true);
210  mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, ms.pole_mass_strs_1_3, i123, true);
211  mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, ms.pole_mass_strs_1_2, i12, true);
212  mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, ms.pole_mass_strs_1_2, i12, true);
213 
214  // Do the lightest Higgs mass separately. The default in most codes is 3 GeV. That seems like
215  // an underestimate if the stop masses are heavy enough, but an overestimate for most points.
216  double rd_mh1 = 2.0 / mssmspec.get(Par::Pole_Mass, ms.h0, 1);
217  mssmspec.set_override(Par::Pole_Mass_1srd_high, rd_mh1, ms.h0, 1, true);
218  mssmspec.set_override(Par::Pole_Mass_1srd_low, rd_mh1, ms.h0, 1, true);
219 
220  // Do the W mass separately. Here we use 10 MeV based on the size of corrections from two-loop papers and advice from Dominik Stockinger.
221  double rd_mW = 0.01 / mssmspec.get(Par::Pole_Mass, "W+");
222  mssmspec.set_override(Par::Pole_Mass_1srd_high, rd_mW, "W+", true);
223  mssmspec.set_override(Par::Pole_Mass_1srd_low, rd_mW, "W+", true);
224 
225  // Save the input value of TanBeta
226  // Probably need to make it a full requirement of the MSSM SpectrumContents
227  if(input_Param.find("TanBeta") != input_Param.end())
228  {
229  mssmspec.set_override(Par::dimensionless, *input_Param.at("TanBeta"), "tanbeta(mZ)", true);
230  }
231 
232  // Create a second SubSpectrum object to wrap the qedqcd object used to initialise the spectrum generator
233  // Attach the sminputs object as well, so that SM pole masses can be passed on (these aren't easily
234  // extracted from the QedQcd object, so use the values that we put into it.)
235  QedQcdWrapper qedqcdspec(oneset,sminputs);
236 
237  // Deal with points where spectrum generator encountered a problem
238  #ifdef SPECBIT_DEBUG
239  std::cout<<"Problem? "<<problems.have_problem()<<std::endl;
240  #endif
241  if( problems.have_problem() )
242  {
243  if( runOptions.getValueOrDef<bool>(false,"invalid_point_fatal") )
244  {
247  std::ostringstream errmsg;
248  errmsg << "A serious problem was encountered during spectrum generation!; ";
249  errmsg << "Message from FlexibleSUSY below:" << std::endl;
250  problems.print_problems(errmsg);
251  problems.print_warnings(errmsg);
252  SpecBit_error().raise(LOCAL_INFO,errmsg.str());
253  }
254  else
255  {
258  std::ostringstream msg;
259  //msg << "";
260  //if( have_bad_mass() ) msg << "bad mass " << std::endl; // TODO: check which one
261  //if( have_tachyon() ) msg << "tachyon" << std::endl;
262  //if( have_thrown() ) msg << "error" << std::endl;
263  //if( have_non_perturbative_parameter() ) msg << "non-perturb. param" << std::endl; // TODO: check which
264  //if( have_failed_pole_mass_convergence() ) msg << "fail pole mass converg." << std::endl; // TODO: check which
265  //if( no_ewsb() ) msg << "no ewsb" << std::endl;
266  //if( no_convergence() ) msg << "no converg." << std::endl;
267  //if( no_perturbative() ) msg << "no pertub." << std::endl;
268  //if( no_rho_convergence() ) msg << "no rho converg." << std::endl;
269  //if( msg.str()=="" ) msg << " Unrecognised problem! ";
270 
272  problems.print_problems(msg);
273  invalid_point().raise(msg.str()); //TODO: This message isn't ending up in the logs.
274  }
275  }
276 
277  if( problems.have_warning() )
278  {
279  std::ostringstream msg;
280  problems.print_warnings(msg);
281  SpecBit_warning().raise(LOCAL_INFO,msg.str()); //TODO: Is a warning the correct thing to do here?
282  }
283 
284  // Write SLHA file (for debugging purposes...)
285  #ifdef SPECBIT_DEBUG
286  typename MI::SlhaIo slha_io;
287  slha_io.set_spinfo(problems);
288  slha_io.set_sminputs(oneset);
289  slha_io.set_minpar(input);
290  slha_io.set_extpar(input);
291  slha_io.set_spectrum(mssmspec.model_interface.model);
292  slha_io.write_to_file("SpecBit/initial_CMSSM_spectrum->slha");
293  #endif
294 
295  // Retrieve any mass cuts
296  static const Spectrum::mc_info mass_cut = runOptions.getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
297  static const Spectrum::mr_info mass_ratio_cut = runOptions.getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
298 
299  // Package QedQcd SubSpectrum object, MSSM SubSpectrum object, and SMInputs struct into a 'full' Spectrum object
300  return Spectrum(qedqcdspec,mssmspec,sminputs,&input_Param,mass_cut,mass_ratio_cut);
301  }
302 
303  //Version for 1.5.1 commented out because we should make it possible to support FS versions in parallel.
304 
305  // template <class MI>
306  // Spectrum run_FS1_5_1_spectrum_generator
307  // ( const typename MI::InputParameters& input
308  // , const SMInputs& sminputs
309  // , const Options& runOptions
310  // , const std::map<str, safe_ptr<const double> >& input_Param
311  // )
312  // {
313  // // SoftSUSY object used to set quark and lepton masses and gauge
314  // // couplings in QEDxQCD effective theory
315  // // Will be initialised by default using values in lowe.h, which we will
316  // // override next.
317  // softsusy::QedQcd oneset;
318 
319  // // Fill QedQcd object with SMInputs values
320  // setup_QedQcd(oneset,sminputs);
321 
322  // // Run everything to Mz
323  // oneset.toMz();
324 
325  // // Create spectrum generator object
326  // typename MI::SpectrumGenerator spectrum_generator;
327 
328  // // Spectrum generator settings
329  // // Default options copied from flexiblesusy/src/spectrum_generator_settings.hpp
330  // //
331  // // | enum | possible values | default value |
332  // // |----------------------------------|------------------------------|-----------------|
333  // // | precision | any positive double | 1.0e-4 |
334  // // | max_iterations | any positive double | 0 (= automatic) |
335  // // | algorithm | 0 (two-scale) or 1 (lattice) | 0 (= two-scale) |
336  // // | calculate_sm_masses | 0 (no) or 1 (yes) | 0 (= no) |
337  // // | pole_mass_loop_order | 0, 1, 2 | 2 (= 2-loop) |
338  // // | ewsb_loop_order | 0, 1, 2 | 2 (= 2-loop) |
339  // // | beta_loop_order | 0, 1, 2 | 2 (= 2-loop) |
340  // // | threshold_corrections_loop_order | 0, 1 | 1 (= 1-loop) |
341  // // | higgs_2loop_correction_at_as | 0, 1 | 1 (= enabled) |
342  // // | higgs_2loop_correction_ab_as | 0, 1 | 1 (= enabled) |
343  // // | higgs_2loop_correction_at_at | 0, 1 | 1 (= enabled) |
344  // // | higgs_2loop_correction_atau_atau | 0, 1 | 1 (= enabled) |
345 
346  // spectrum_generator.set_precision_goal (runOptions.getValueOrDef<double>(1.0e-4,"precision_goal"));
347  // spectrum_generator.set_max_iterations (runOptions.getValueOrDef<double>(0, "max_iterations"));
348  // spectrum_generator.set_calculate_sm_masses (runOptions.getValueOrDef<bool> (false, "calculate_sm_masses"));
349  // spectrum_generator.set_pole_mass_loop_order (runOptions.getValueOrDef<int> (2, "pole_mass_loop_order"));
350  // spectrum_generator.set_ewsb_loop_order (runOptions.getValueOrDef<int> (2, "ewsb_loop_order"));
351  // spectrum_generator.set_beta_loop_order (runOptions.getValueOrDef<int> (2, "beta_loop_order"));
352  // spectrum_generator.set_threshold_corrections_loop_order(runOptions.getValueOrDef<int> (2, "threshold_corrections_loop_order"));
353 
354  // // Higgs loop corrections are a little different... sort them out now
355  // Two_loop_corrections two_loop_settings;
356 
357  // // alpha_t alpha_s
358  // // alpha_b alpha_s
359  // // alpha_t^2 + alpha_t alpha_b + alpha_b^2
360  // // alpha_tau^2
361  // two_loop_settings.higgs_at_as = runOptions.getValueOrDef<bool>(true,"use_higgs_2loop_at_as");
362  // two_loop_settings.higgs_ab_as = runOptions.getValueOrDef<bool>(true,"use_higgs_2loop_ab_as");
363  // two_loop_settings.higgs_at_at = runOptions.getValueOrDef<bool>(true,"use_higgs_2loop_at_at");
364  // two_loop_settings.higgs_atau_atau = runOptions.getValueOrDef<bool>(true,"use_higgs_2loop_atau_atau");
365 
366  // spectrum_generator.set_two_loop_corrections(two_loop_settings);
367 
368  // // Generate spectrum
369  // spectrum_generator.run(oneset, input);
370 
371  // // Extract report on problems...
372  // const typename MI::Problems& problems = spectrum_generator.get_problems();
373 
374  // // Create Model_interface to carry the input and results, and know
375  // // how to access the flexiblesusy routines.
376  // // Note: Output of spectrum_generator.get_model() returns type, e.g. CMSSM.
377  // // Need to convert it to type CMSSM_slha (which alters some conventions of
378  // // parameters into SLHA format)
379  // MI model_interface(spectrum_generator,oneset,input);
380 
381  // // Create SubSpectrum object to wrap flexiblesusy data
382  // // THIS IS STATIC so that it lives on once we leave this module function. We
383  // // therefore cannot run the same spectrum generator twice in the same loop and
384  // // maintain the spectrum resulting from both. But we should never want to do
385  // // this.
386  // // A pointer to this object is what gets turned into a SubSpectrum pointer and
387  // // passed around Gambit.
388  // //
389  // // This object will COPY the interface data members into itself, so it is now the
390  // // one-stop-shop for all spectrum information, including the model interface object.
391  // MSSMSpec<MI> mssmspec(model_interface, "FlexibleSUSY", "1.5.1");
392 
393  // // Add extra information about the scales used to the wrapper object
394  // // (last parameter turns on the 'allow_new' option for the override setter, which allows
395  // // us to set parameters that don't previously exist)
396  // mssmspec.set_override(Par::mass1,spectrum_generator.get_high_scale(),"high_scale",true);
397  // mssmspec.set_override(Par::mass1,spectrum_generator.get_susy_scale(),"susy_scale",true);
398  // mssmspec.set_override(Par::mass1,spectrum_generator.get_low_scale(), "low_scale", true);
399 
400  // // Add theory errors
401  // static const MSSM_strs ms;
402 
403  // static const std::vector<int> i12 = initVector(1,2);
404  // static const std::vector<int> i123 = initVector(1,2,3);
405  // static const std::vector<int> i1234 = initVector(1,2,3,4);
406  // static const std::vector<int> i123456 = initVector(1,2,3,4,5,6);
407 
408  // // 3% theory "error"
409  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, ms.pole_mass_pred, true);
410  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, ms.pole_mass_pred, true);
411  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, ms.pole_mass_strs_1_6, i123456, true);
412  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, ms.pole_mass_strs_1_6, i123456, true);
413  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, "~chi0", i1234, true);
414  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, "~chi0", i1234, true);
415  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, ms.pole_mass_strs_1_3, i123, true);
416  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, ms.pole_mass_strs_1_3, i123, true);
417  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_high, 0.03, ms.pole_mass_strs_1_2, i12, true);
418  // mssmspec.set_override_vector(Par::Pole_Mass_1srd_low, 0.03, ms.pole_mass_strs_1_2, i12, true);
419 
420  // // Do the lightest Higgs mass separately. The default in most codes is 3 GeV. That seems like
421  // // an underestimate if the stop masses are heavy enough, but an overestimate for most points.
422  // double rd_mh1 = 2.0 / mssmspec.get(Par::Pole_Mass, ms.h0, 1);
423  // mssmspec.set_override(Par::Pole_Mass_1srd_high, rd_mh1, ms.h0, 1, true);
424  // mssmspec.set_override(Par::Pole_Mass_1srd_low, rd_mh1, ms.h0, 1, true);
425 
426  // // Do the W mass separately. Here we use 10 MeV based on the size of corrections from two-loop papers and advice from Dominik Stockinger.
427  // double rd_mW = 0.01 / mssmspec.get(Par::Pole_Mass, "W+");
428  // mssmspec.set_override(Par::Pole_Mass_1srd_high, rd_mW, "W+", true);
429  // mssmspec.set_override(Par::Pole_Mass_1srd_low, rd_mW, "W+", true);
430 
431  // // Save the input value of TanBeta
432  // // Probably need to make it a full requirement of the MSSM SpectrumContents
433  // if(input_Param.find("TanBeta") != input_Param.end())
434  // {
435  // mssmspec.set_override(Par::dimensionless, *input_Param.at("TanBeta"), "tanbeta(mZ)", true);
436  // }
437 
438  // // Create a second SubSpectrum object to wrap the qedqcd object used to initialise the spectrum generator
439  // // Attach the sminputs object as well, so that SM pole masses can be passed on (these aren't easily
440  // // extracted from the QedQcd object, so use the values that we put into it.)
441  // QedQcdWrapper qedqcdspec(oneset,sminputs);
442 
443  // // Deal with points where spectrum generator encountered a problem
444  // #ifdef SPECBIT_DEBUG
445  // std::cout<<"Problem? "<<problems.have_problem()<<std::endl;
446  // #endif
447  // if( problems.have_problem() )
448  // {
449  // if( runOptions.getValueOrDef<bool>(false,"invalid_point_fatal") )
450  // {
451  // ///TODO: Need to tell gambit that the spectrum is not viable somehow. For now
452  // /// just die.
453  // std::ostringstream errmsg;
454  // errmsg << "A serious problem was encountered during spectrum generation!; ";
455  // errmsg << "Message from FlexibleSUSY below:" << std::endl;
456  // problems.print_problems(errmsg);
457  // problems.print_warnings(errmsg);
458  // SpecBit_error().raise(LOCAL_INFO,errmsg.str());
459  // }
460  // else
461  // {
462  // /// Check what the problem was
463  // /// see: contrib/MassSpectra/flexiblesusy/src/problems.hpp
464  // std::ostringstream msg;
465  // //msg << "";
466  // //if( have_bad_mass() ) msg << "bad mass " << std::endl; // TODO: check which one
467  // //if( have_tachyon() ) msg << "tachyon" << std::endl;
468  // //if( have_thrown() ) msg << "error" << std::endl;
469  // //if( have_non_perturbative_parameter() ) msg << "non-perturb. param" << std::endl; // TODO: check which
470  // //if( have_failed_pole_mass_convergence() ) msg << "fail pole mass converg." << std::endl; // TODO: check which
471  // //if( no_ewsb() ) msg << "no ewsb" << std::endl;
472  // //if( no_convergence() ) msg << "no converg." << std::endl;
473  // //if( no_perturbative() ) msg << "no pertub." << std::endl;
474  // //if( no_rho_convergence() ) msg << "no rho converg." << std::endl;
475  // //if( msg.str()=="" ) msg << " Unrecognised problem! ";
476 
477  // /// Fast way for now:
478  // problems.print_problems(msg);
479  // invalid_point().raise(msg.str()); //TODO: This message isn't ending up in the logs.
480  // }
481  // }
482 
483  // if( problems.have_warning() )
484  // {
485  // std::ostringstream msg;
486  // problems.print_warnings(msg);
487  // SpecBit_warning().raise(LOCAL_INFO,msg.str()); //TODO: Is a warning the correct thing to do here?
488  // }
489 
490  // // Write SLHA file (for debugging purposes...)
491  // #ifdef SPECBIT_DEBUG
492  // typename MI::SlhaIo slha_io;
493  // slha_io.set_spinfo(problems);
494  // slha_io.set_sminputs(oneset);
495  // slha_io.set_minpar(input);
496  // slha_io.set_extpar(input);
497  // slha_io.set_spectrum(mssmspec.model_interface.model);
498  // slha_io.write_to_file("SpecBit/initial_CMSSM_spectrum->slha");
499  // #endif
500 
501  // // Retrieve any mass cuts
502  // static const Spectrum::mc_info mass_cut = runOptions.getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
503  // static const Spectrum::mr_info mass_ratio_cut = runOptions.getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
504 
505  // // Package QedQcd SubSpectrum object, MSSM SubSpectrum object, and SMInputs struct into a 'full' Spectrum object
506  // return Spectrum(qedqcdspec,mssmspec,sminputs,&input_Param,mass_cut,mass_ratio_cut);
507  // }
508 
510  // Names must conform to convention "<parname>_ij"
511  Eigen::Matrix<double,3,3> fill_3x3_parameter_matrix(const std::string& rootname, const std::map<str, safe_ptr<const double> >& Param)
512  {
513  Eigen::Matrix<double,3,3> output;
514  for(int i=0; i<3; ++i) for(int j=0; j<3; ++j)
515  {
516  output(i,j) = *Param.at(rootname + "_" + std::to_string(i+1) + std::to_string(j+1));
517  }
518  return output;
519  }
520 
522  Eigen::Matrix<double,3,3> fill_3x3_symmetric_parameter_matrix(const std::string& rootname, const std::map<str, safe_ptr<const double> >& Param)
523  {
524  Eigen::Matrix<double,3,3> output;
525  for(int i=0; i<3; ++i) for(int j=0; j<3; ++j)
526  {
527  str parname = rootname + "_" + ( i < j ? std::to_string(i+1) + std::to_string(j+1) : std::to_string(j+1) + std::to_string(i+1));
528  output(i,j) = *Param.at(parname);
529  }
530  return output;
531  }
532 
535  template <class T>
536  void fill_MSSM63_input(T& input, const std::map<str, safe_ptr<const double> >& Param )
537  {
538  //double valued parameters
539  input.TanBeta = *Param.at("TanBeta");
540  input.MassBInput = *Param.at("M1");
541  input.MassWBInput = *Param.at("M2");
542  input.MassGInput = *Param.at("M3");
543 
544  // Sanity checks
545  if(input.TanBeta<0)
546  {
547  std::ostringstream msg;
548  msg << "Tried to set TanBeta parameter to a negative value ("<<input.TanBeta<<")! This parameter must be positive. Please check your inifile and try again.";
549  SpecBit_error().raise(LOCAL_INFO,msg.str());
550  }
551 
552  //3x3 matrices; filled with the help of a convenience function
553  input.mq2Input = fill_3x3_symmetric_parameter_matrix("mq2", Param);
554  input.ml2Input = fill_3x3_symmetric_parameter_matrix("ml2", Param);
555  input.md2Input = fill_3x3_symmetric_parameter_matrix("md2", Param);
556  input.mu2Input = fill_3x3_symmetric_parameter_matrix("mu2", Param);
557  input.me2Input = fill_3x3_symmetric_parameter_matrix("me2", Param);
558  input.Aeij = fill_3x3_parameter_matrix("Ae", Param);
559  input.Adij = fill_3x3_parameter_matrix("Ad", Param);
560  input.Auij = fill_3x3_parameter_matrix("Au", Param);
561 
562  #ifdef SPECBIT_DEBUG
563  #define INPUT(p) input.p
564  #define ostr std::cout
565  #define oend std::endl
566  ostr << "TanBeta = " << INPUT(TanBeta) << ", " << oend ;
567  ostr << "mq2Input = " << oend << INPUT(mq2Input) << ", " << oend;
568  ostr << "ml2Input = " << oend << INPUT(ml2Input) << ", " << oend;
569  ostr << "md2Input = " << oend << INPUT(md2Input) << ", " << oend;
570  ostr << "mu2Input = " << oend << INPUT(mu2Input) << ", " << oend;
571  ostr << "me2Input = " << oend << INPUT(me2Input) << ", " << oend;
572  ostr << "MassBInput = " << INPUT(MassBInput) << ", " << oend;
573  ostr << "MassWBInput = " << INPUT(MassWBInput) << ", " << oend;
574  ostr << "MassGInput = " << INPUT(MassGInput) << ", " << oend;
575  ostr << "Aeij = " << oend << INPUT(Aeij) << ", " << oend;
576  ostr << "Adij = " << oend << INPUT(Adij) << ", " << oend;
577  ostr << "Auij = " << oend << INPUT(Auij) << ", " << oend;
578  #undef INPUT
579  #undef ostr
580  #undef oend
581  #endif
582  }
583 
584 
585  // Similar to above, except this is for MSSMEFTHiggs spectrum
586  // generator and a few others where different names for inputs for many
587  // parameters are used. This should be standardised.
588 
589  template <class T>
590  void fill_MSSM63_input_altnames(T& input, const std::map<str, safe_ptr<const double> >& Param )
591  {
592  //double valued parameters
593  input.TanBeta = *Param.at("TanBeta");
594  input.M1Input = *Param.at("M1");
595  input.M2Input = *Param.at("M2");
596  input.M3Input = *Param.at("M3");
597 
598  // Sanity checks
599  if(input.TanBeta<0)
600  {
601  std::ostringstream msg;
602  msg << "Tried to set TanBeta parameter to a negative value ("<<input.TanBeta<<")! This parameter must be positive. Please check your inifile and try again.";
603  SpecBit_error().raise(LOCAL_INFO,msg.str());
604  }
605 
606  //3x3 matrices; filled with the help of a convenience function
607  input.mq2Input = fill_3x3_symmetric_parameter_matrix("mq2", Param);
608  input.ml2Input = fill_3x3_symmetric_parameter_matrix("ml2", Param);
609  input.md2Input = fill_3x3_symmetric_parameter_matrix("md2", Param);
610  input.mu2Input = fill_3x3_symmetric_parameter_matrix("mu2", Param);
611  input.me2Input = fill_3x3_symmetric_parameter_matrix("me2", Param);
612  input.AeInput = fill_3x3_parameter_matrix("Ae", Param);
613  input.AdInput= fill_3x3_parameter_matrix("Ad", Param);
614  input.AuInput = fill_3x3_parameter_matrix("Au", Param);
615 
616  #ifdef SPECBIT_DEBUG
617  #define INPUT(p) input.p
618  #define ostr std::cout
619  #define oend std::endl
620  ostr << "TanBeta = " << INPUT(TanBeta) << ", " << oend ;
621  ostr << "mq2Input = " << oend << INPUT(mq2Input) << ", " << oend;
622  ostr << "ml2Input = " << oend << INPUT(ml2Input) << ", " << oend;
623  ostr << "md2Input = " << oend << INPUT(md2Input) << ", " << oend;
624  ostr << "mu2Input = " << oend << INPUT(mu2Input) << ", " << oend;
625  ostr << "me2Input = " << oend << INPUT(me2Input) << ", " << oend;
626  ostr << "M1Input = " << INPUT(M1Input) << ", " << oend;
627  ostr << "M2Input = " << INPUT(M2Input) << ", " << oend;
628  ostr << "M3Input = " << INPUT(M3Input) << ", " << oend;
629  ostr << "AeInput = " << oend << INPUT(AeInput) << ", " << oend;
630  ostr << "AdInput = " << oend << INPUT(AdInput) << ", " << oend;
631  ostr << "AuInput = " << oend << INPUT(AuInput) << ", " << oend;
632  #undef INPUT
633  #undef ostr
634  #undef oend
635  #endif
636 
637  }
638 
639 
641  bool has_neutralino_LSP(const Spectrum& result)
642  {
643  double msqd = result.get(Par::Pole_Mass, 1000001, 0);
644  double msqu = result.get(Par::Pole_Mass, 1000002, 0);
645  double msl = result.get(Par::Pole_Mass, 1000011, 0);
646  double msneu = result.get(Par::Pole_Mass, 1000012, 0);
647  double mglui = result.get(Par::Pole_Mass, 1000021, 0);
648  double mchi0 = std::abs(result.get(Par::Pole_Mass, 1000022, 0));
649  double mchip = std::abs(result.get(Par::Pole_Mass, 1000024, 0));
650 
651  return mchi0 < mchip &&
652  mchi0 < mglui &&
653  mchi0 < msl &&
654  mchi0 < msneu &&
655  mchi0 < msqu &&
656  mchi0 < msqd;
657  }
659  bool has_neutralino_LSP(const Spectrum* result)
660  {
661  return has_neutralino_LSP(*result);
662  }
663 
664 
666 
667 
669  // =======================================================================
670  // These are wrapped up in Gambit functor objects according to the
671  // instructions in the rollcall header
672 
673  // Functions to change the capability associated with a Spectrum object to "SM_spectrum"
676  void convert_E6MSSM_to_SM (Spectrum &result) {result = *Pipes::convert_E6MSSM_to_SM::Dep::E6MSSM_spectrum;}
677 
679  {
680  namespace myPipe = Pipes::get_MSSM_spectrum_SPheno;
681  const SMInputs &sminputs = *myPipe::Dep::SMINPUTS;
682 
683  // Set up the input structure
684  Finputs inputs;
685  inputs.sminputs = sminputs;
686  inputs.param = myPipe::Param;
687  inputs.options = myPipe::runOptions;
688 
689  // Retrieve any mass cuts
690  static const Spectrum::mc_info mass_cut = myPipe::runOptions->getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
691  static const Spectrum::mr_info mass_ratio_cut = myPipe::runOptions->getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
692 
693  // Get the spectrum from the Backend
694  myPipe::BEreq::SPheno_MSSMspectrum(spectrum, inputs);
695 
696  // Only allow neutralino LSPs.
697  if (not has_neutralino_LSP(spectrum)) invalid_point().raise("Neutralino is not LSP.");
698 
699  // Drop SLHA files if requested
700  spectrum.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
701 
702  }
703 
704  // Runs FlexibleSUSY MSSMEFTHiggs model spectrum generator with SUSY
705  // scale boundary conditions, ie accepts MSSM parameters at MSUSY,
706  // and has DRbar mA and mu as an input and mHu2 and mHd2 as EWSB
707  // outputs, so it is for the MSSMatMSUSY_mA model.
708  #if(FS_MODEL_MSSMatMSUSYEFTHiggs_mAmu_IS_BUILT)
709  void get_MSSMatMSUSY_mA_spectrum_FlexibleEFTHiggs (Spectrum& result)
710  {
711  // Access the pipes for this function to get model and parameter information
712  namespace myPipe = Pipes::get_MSSMatMSUSY_mA_spectrum_FlexibleEFTHiggs;
713 
714  // Get SLHA2 SMINPUTS values
715  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
716 
717  // Get input parameters (from flexiblesusy namespace)
718  MSSMatMSUSYEFTHiggs_mAmu_input_parameters input;
719  input.MuInput = *myPipe::Param.at("mu");
720  // This FS spectrum generator has mA as the parameter
721  input.mAInput = *myPipe::Param.at("mA");
722  fill_MSSM63_input_altnames(input,myPipe::Param); // Fill the rest
723  result = run_FS_spectrum_generator<MSSMatMSUSYEFTHiggs_mAmu_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
724 
725  // Only allow neutralino LSPs.
726  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
727 
728  // Drop SLHA files if requested
729  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
730 
731  }
732  #endif
733 
734  // Runs FlexibleSUSY MSSMEFTHiggs model spectrum generator
735  // and has m3^2 and mu as EWSB outputs, so it is for the
736  // MSSMatQ_model.
737  #if(FS_MODEL_MSSMEFTHiggs_IS_BUILT)
738  void get_MSSMatQ_spectrum_FlexibleEFTHiggs (Spectrum& result)
739  {
740  // Access the pipes for this function to get model and parameter information
741  namespace myPipe = Pipes::get_MSSMatQ_spectrum_FlexibleEFTHiggs;
742 
743  // Get SLHA2 SMINPUTS values
744  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
745 
746  // Get input parameters (from flexiblesusy namespace)
747  MSSMEFTHiggs_input_parameters input;
748  // MSSMatQ also requires input scale to be supplied with name MSUSY
749  input.MSUSY = *myPipe::Param.at("Qin");
750  input.mHu2IN = *myPipe::Param.at("mHu2");
751  input.mHd2IN = *myPipe::Param.at("mHd2");
752  input.SignMu = *myPipe::Param.at("SignMu");
753  fill_MSSM63_input_altnames(input,myPipe::Param); // Fill the rest
754  result = run_FS_spectrum_generator<MSSMEFTHiggs_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
755 
756  // Only allow neutralino LSPs.
757  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
758 
759  // Drop SLHA files if requested
760  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
761 
762  }
763  #endif
764 
765  // Runs FlexibleSUSY MSSMEFTHiggs_mAmu spectrum generator with
766  // boundary conditions at a user specified scale, ie accepts MSSM
767  // parameters at Q, and has DRbar mA and mu as an input and mHu2
768  // and mHd2 as EWSB outputs, so it is for the MSSMatMSUSY_mA model.
769  #if(FS_MODEL_MSSMEFTHiggs_mAmu_IS_BUILT)
770  void get_MSSMatQ_mA_spectrum_FlexibleEFTHiggs (Spectrum& result)
771  {
772  // Access the pipes for this function to get model and parameter information
773  namespace myPipe = Pipes::get_MSSMatQ_mA_spectrum_FlexibleEFTHiggs;
774 
775  // Get SLHA2 SMINPUTS values
776  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
777 
778  // Get input parameters (from flexiblesusy namespace)
779  MSSMEFTHiggs_mAmu_input_parameters input;
780  input.MuInput = *myPipe::Param.at("mu");
781  // This FS spectrum generator has mA as the parameter
782  input.mAInput = *myPipe::Param.at("mA");
783  // Note: Qin has been named MSUSY inside the spectrum generator
784  // but it is a user-input scale in this case.
785  input.MSUSY = *myPipe::Param.at("Qin");
786  // Fill the rest.
787  // Note: This particular spectrum generator has been created with
788  // different names for parameter inputs. We should standardise this
789  fill_MSSM63_input_altnames(input,myPipe::Param);
790  result = run_FS_spectrum_generator<MSSMEFTHiggs_mAmu_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
791 
792  // Only allow neutralino LSPs.
793  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
794 
795  // Drop SLHA files if requested
796  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
797 
798  }
799  #endif
800 
801 
802 
803  // Runs FlexibleSUSY MSSM spectrum generator with CMSSM (GUT scale) boundary conditions
804  // In principle an identical spectrum can be obtained from the function
805  // get_MSSMatGUT_spectrum_FS
806  // by setting the input parameters to match the CMSSM assumptions
807  #if(FS_MODEL_CMSSM_IS_BUILT)
808  void get_CMSSM_spectrum_FS (Spectrum& result)
809  {
810 
811  // Access the pipes for this function to get model and parameter information
812  namespace myPipe = Pipes::get_CMSSM_spectrum_FS;
813 
814  // Get SLHA2 SMINPUTS values
815  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
816 
817  // Get input parameters (from flexiblesusy namespace)
818  CMSSM_input_parameters input;
819 
820  input.m0 = *myPipe::Param["M0"];
821  input.m12 = *myPipe::Param["M12"];
822  input.TanBeta = *myPipe::Param["TanBeta"];
823  input.SignMu = *myPipe::Param["SignMu"];
824  input.Azero = *myPipe::Param["A0"];
825 
826  // Sanity checks
827  if(input.TanBeta<0)
828  {
829  std::ostringstream msg;
830  msg << "Tried to set TanBeta parameter to a negative value ("<<input.TanBeta<<")! This parameter must be positive. Please check your inifile and try again.";
831  SpecBit_error().raise(LOCAL_INFO,msg.str());
832  }
833  if(input.SignMu!=-1 and input.SignMu!=1)
834  {
835  std::ostringstream msg;
836  msg << "Tried to set SignMu parameter to a value that is not a sign! ("<<input.SignMu<<")! This parameter must be set to either 1 or -1. Please check your inifile and try again.";
837  SpecBit_error().raise(LOCAL_INFO,msg.str());
838  }
839 
840  // Run spectrum generator
841  result = run_FS_spectrum_generator<CMSSM_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
842 
843  // Only allow neutralino LSPs.
844  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
845 
846  // Drop SLHA files if requested
847  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
848 
849  }
850  #endif
851 
852  // Runs FlexibleSUSY MSSM spectrum generator with EWSB scale input (boundary conditions)
853  #if(FS_MODEL_MSSM_IS_BUILT)
854  void get_MSSMatQ_spectrum_FS (Spectrum& result)
855  {
856  using namespace softsusy;
857  namespace myPipe = Pipes::get_MSSMatQ_spectrum_FS;
858  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
859  MSSM_input_parameters input;
860  input.Qin = *myPipe::Param.at("Qin"); // MSSMatQ also requires input scale to be supplied
861  input.mHu2IN = *myPipe::Param.at("mHu2");
862  input.mHd2IN = *myPipe::Param.at("mHd2");
863  input.SignMu = *myPipe::Param.at("SignMu");
864  if(input.SignMu!=-1 and input.SignMu!=1)
865  {
866  std::ostringstream msg;
867  msg << "Tried to set SignMu parameter to a value that is not a sign! ("<<input.SignMu<<")! This parameter must be set to either 1 or -1. Please check your inifile and try again.";
868  SpecBit_error().raise(LOCAL_INFO,msg.str());
869  }
870  fill_MSSM63_input(input,myPipe::Param);
871  result = run_FS_spectrum_generator<MSSM_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
872  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
873  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
874  }
875  #endif
876 
877  // Runs FlexibleSUSY MSSM spectrum generator with EWSB scale input (boundary conditions)
878  // but with mA and mu as parameters instead of mHu2 and mHd2
879  #if(FS_MODEL_MSSM_mAmu_IS_BUILT)
880  void get_MSSMatQ_mA_spectrum_FS (Spectrum& result)
881  {
882  using namespace softsusy;
883  namespace myPipe = Pipes::get_MSSMatQ_mA_spectrum_FS;
884  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
885  MSSM_mAmu_input_parameters input;
886  input.Qin = *myPipe::Param.at("Qin"); // MSSMatQ also requires input scale to be supplied
887  input.MuInput = *myPipe::Param.at("mu");
888  // Note this spectrum generator mA2 is the parameter.
889  // However this freedom is not used in GAMBIT
890  // and not needed as mA is a DRbar mass eigenstate for a scalar
891  double mA = *myPipe::Param.at("mA");
892  input.mA2Input = mA*mA;
893  fill_MSSM63_input(input,myPipe::Param); // Fill the rest
894  result = run_FS_spectrum_generator<MSSM_mAmu_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
895  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
896  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
897  }
898  #endif
899 
900  // Runs FlexibleSUSY MSSM spectrum generator with GUT scale input (boundary conditions)
901  #if(FS_MODEL_MSSMatMGUT_IS_BUILT)
902  void get_MSSMatMGUT_spectrum_FS (Spectrum& result)
903  {
904  using namespace softsusy;
905  namespace myPipe = Pipes::get_MSSMatMGUT_spectrum_FS;
906  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
907  MSSMatMGUT_input_parameters input;
908  input.mHu2IN = *myPipe::Param.at("mHu2");
909  input.mHd2IN = *myPipe::Param.at("mHd2");
910  input.SignMu = *myPipe::Param.at("SignMu");
911  if(input.SignMu!=-1 and input.SignMu!=1)
912  {
913  std::ostringstream msg;
914  msg << "Tried to set SignMu parameter to a value that is not a sign! ("<<input.SignMu<<")! This parameter must be set to either 1 or -1. Please check your inifile and try again.";
915  SpecBit_error().raise(LOCAL_INFO,msg.str());
916  }
917  fill_MSSM63_input(input,myPipe::Param);
918  result = run_FS_spectrum_generator<MSSMatMGUT_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
919  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
920  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
921  }
922  #endif
923 
924  // Runs FlexibleSUSY MSSMatMGUTEFTHiggs model spectrum generator
925  // and has m3^2 and mu as EWSB outputs, so it is for the
926  // MSSMatMGUT_model.
927  #if(FS_MODEL_MSSMatMGUTEFTHiggs_IS_BUILT)
928  void get_MSSMatMGUT_spectrum_FlexibleEFTHiggs (Spectrum& result)
929  {
930  // Access the pipes for this function to get model and parameter information
931  namespace myPipe = Pipes::get_MSSMatMGUT_spectrum_FlexibleEFTHiggs;
932 
933  // Get SLHA2 SMINPUTS values
934  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
935 
936  // Get input parameters (from flexiblesusy namespace)
937  MSSMatMGUTEFTHiggs_input_parameters input;
938  input.mHu2IN = *myPipe::Param.at("mHu2");
939  input.mHd2IN = *myPipe::Param.at("mHd2");
940  input.SignMu = *myPipe::Param.at("SignMu");
941  if(input.SignMu!=-1 and input.SignMu!=1)
942  {
943  std::ostringstream msg;
944  msg << "Tried to set SignMu parameter to a value that is not a sign! ("<<input.SignMu<<")! This parameter must be set to either 1 or -1. Please check your inifile and try again.";
945  SpecBit_error().raise(LOCAL_INFO,msg.str());
946  }
947 
948  fill_MSSM63_input(input,myPipe::Param); // Fill the rest
949  result = run_FS_spectrum_generator<MSSMatMGUTEFTHiggs_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
950 
951  // Only allow neutralino LSPs.
952  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
953 
954  // Drop SLHA files if requested
955  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
956 
957  }
958  #endif
959 
960 
961  // Runs FlexibleSUSY MSSM spectrum generator with GUT scale input (boundary conditions)
962  // but with mA and mu as parameters instead of mHu2 and mHd2
963  #if(FS_MODEL_MSSMatMGUT_mAmu_IS_BUILT)
964  void get_MSSMatMGUT_mA_spectrum_FS (Spectrum& result)
965  {
966  using namespace softsusy;
967  namespace myPipe = Pipes::get_MSSMatMGUT_mA_spectrum_FS;
968  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
969  MSSMatMGUT_mAmu_input_parameters input;
970  input.MuInput = *myPipe::Param.at("mu");
971  // Note this spectrum generator mA2 is the parameter.
972  // However this freedom is not used in GAMBIT
973  // and not needed as mA is a DRbar mass eigenstate for a scalar
974  double mA = *myPipe::Param.at("mA");
975  input.mA2Input = mA*mA;
976  fill_MSSM63_input(input,myPipe::Param); // Fill the rest
977  result = run_FS_spectrum_generator<MSSMatMGUT_mAmu_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
978  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
979  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
980  }
981  #endif
982 
983  // Runs FlexibleSUSY MSSMatMGUTEFTHiggs model spectrum generator
984  // and has m3^2 and mu as EWSB outputs, so it is for the
985  // MSSMatMGUT_model.
986  #if(FS_MODEL_MSSMatMGUTEFTHiggs_mAmu_IS_BUILT)
987  void get_MSSMatMGUT_mA_spectrum_FlexibleEFTHiggs (Spectrum& result)
988  {
989  // Access the pipes for this function to get model and parameter information
990  namespace myPipe = Pipes::get_MSSMatMGUT_mA_spectrum_FlexibleEFTHiggs;
991 
992  // Get SLHA2 SMINPUTS values
993  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
994 
995  // Get input parameters (from flexiblesusy namespace)
996  MSSMatMGUTEFTHiggs_mAmu_input_parameters input;
997  input.MuInput = *myPipe::Param.at("mu");
998  // Note this spectrum generator mA2 is the parameter.
999  // However this freedom is not used in GAMBIT
1000  // and not needed as mA is a DRbar mass eigenstate for a scalar
1001  double mA = *myPipe::Param.at("mA");
1002  input.mA2Input = mA*mA;
1003 
1004  fill_MSSM63_input(input,myPipe::Param); // Fill the rest
1005  result = run_FS_spectrum_generator<MSSMatMGUTEFTHiggs_mAmu_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
1006 
1007  // Only allow neutralino LSPs.
1008  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
1009 
1010  // Drop SLHA files if requested
1011  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
1012 
1013  }
1014  #endif
1015 
1016 
1017  // Runs FlexibleSUSY MSSM spectrum generator with SUSY scale input (boundary conditions)
1018  // but with mA and mu as parameters instead of mHu2 and mHd2
1019  #if(FS_MODEL_MSSMatMSUSY_mAmu_IS_BUILT)
1020  void get_MSSMatMSUSY_mA_spectrum_FS (Spectrum& result)
1021  {
1022  using namespace softsusy;
1023  namespace myPipe = Pipes::get_MSSMatMSUSY_mA_spectrum_FS;
1024  const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
1025  MSSMatMSUSY_mAmu_input_parameters input;
1026  input.MuInput = *myPipe::Param.at("mu");
1027  // Note this spectrum generator mA2 is the parameter.
1028  // However this freedom is not used in GAMBIT
1029  // and not needed as mA is a DRbar mass eigenstate for a scalar
1030  double mA = *myPipe::Param.at("mA");
1031  input.mA2Input = mA*mA; // FS has mA^2 as the parameter
1032  fill_MSSM63_input(input,myPipe::Param); // Fill the rest
1033  result = run_FS_spectrum_generator<MSSMatMSUSY_mAmu_interface<ALGORITHM1>>(input,sminputs,*myPipe::runOptions,myPipe::Param);
1034  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
1035  result.drop_SLHAs_if_requested(myPipe::runOptions, "GAMBIT_unimproved_spectrum");
1036  }
1037  #endif
1038 
1039  void get_GUTMSSMB_spectrum (Spectrum &/*result*/)
1040  {
1041  // Placeholder
1042  }
1043 
1046 
1051  {
1053  const Spectrum& matched_spectra(*myPipe::Dep::unimproved_MSSM_spectrum);
1054  result = &matched_spectra.get_LE();
1055  }
1056 
1059  {
1061  }
1062 
1065  {
1067  }
1068 
1074  {
1075  // Static counter running in a loop over all filenames
1076  static unsigned int counter = 0;
1077  static long int ncycle = 1;
1078  SLHAstruct input_slha;
1079 
1080  namespace myPipe = Pipes::get_MSSM_spectrum_from_SLHAfile;
1081 
1082  // Read filename from yaml file
1083  std::vector<std::string> filenames =
1084  myPipe::runOptions->getValue<std::vector<std::string>>("filenames");
1085 
1086  // Check how many loop over the input files we are doing.
1087  long int cycles = myPipe::runOptions->getValueOrDef<int>(-1,"cycles");
1088 
1089  // Check if we have completed the requested number of cycles
1090  if(cycles>0 and ncycle>cycles)
1091  {
1092  std::ostringstream msg;
1093  msg << "Preset number of loops through input files reached! Stopping. (tried to start cycle "<<ncycle<<" of "<<cycles<<")";
1094  SpecBit_error().raise(LOCAL_INFO,msg.str());
1095  }
1096 
1097  std::string filename = filenames[counter];
1098 
1099  logger() << "Reading SLHA file: " << filename << EOM;
1100  std::ifstream ifs(filename.c_str());
1101  if(!ifs.good()){ SpecBit_error().raise(LOCAL_INFO,"ERROR: SLHA file not found."); }
1102  ifs >> input_slha;
1103  ifs.close();
1104  counter++;
1105  if( counter >= filenames.size() )
1106  {
1107  logger() << "Returning to start of input SLHA file list (finished "<<ncycle<<" cycles)" << EOM;
1108  counter = 0;
1109  ncycle++;
1110  }
1111 
1112  // Retrieve any mass cuts
1113  static const Spectrum::mc_info mass_cut = myPipe::runOptions->getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
1114  static const Spectrum::mr_info mass_ratio_cut = myPipe::runOptions->getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
1115 
1116  // Create Spectrum object from the slhaea object
1117  result = spectrum_from_SLHAea<MSSMSimpleSpec, SLHAstruct>(input_slha, input_slha, mass_cut, mass_ratio_cut);
1118 
1119  // Add getter for susy scale if option set for this
1120  bool add_susy_scale = myPipe::runOptions->getValueOrDef<bool>(false,"assume_Q_is_MSUSY");
1121  if(add_susy_scale)
1122  {
1123  result.get_HE().set_override(Par::mass1,result.get_HE().GetScale(),"susy_scale",true);
1124  }
1125 
1126  // No sneaking in charged LSPs via SLHA, jävlar.
1127  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
1128  }
1129 
1134  {
1135  namespace myPipe = Pipes::get_MSSM_spectrum_from_SLHAstruct;
1136  const SLHAstruct& input_slha_tmp = *myPipe::Dep::unimproved_MSSM_spectrum; // Retrieve dependency on SLHAstruct
1137 
1139  SLHAstruct input_slha(input_slha_tmp); // Copy struct (for demo adding of GAMBIT block only)
1140  // For example; add this to your input SLHAstruct:
1141  input_slha["GAMBIT"][""] << "BLOCK" << "GAMBIT";
1142  input_slha["GAMBIT"][""] << 1 << 1e99 << "# Input scale";
1143 
1144  // Retrieve any mass cuts
1145  static const Spectrum::mc_info mass_cut = myPipe::runOptions->getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
1146  static const Spectrum::mr_info mass_ratio_cut = myPipe::runOptions->getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");
1147 
1148  // Create Spectrum object from the slhaea object
1149  result = spectrum_from_SLHAea<MSSMSimpleSpec, SLHAstruct>(input_slha, input_slha, mass_cut, mass_ratio_cut);
1150 
1151  // No sneaking in charged LSPs via SLHA, jävlar.
1152  if (not has_neutralino_LSP(result)) invalid_point().raise("Neutralino is not LSP.");
1153 
1154  // In order to translate from e.g. MSSM63atMGUT to MSSM63atQ, we need
1155  // to know that input scale Q. This is generally not stored in SLHA format,
1156  // but we need it, so if you want to produce a Spectrum object this way you
1157  // will need to add this information to your SLHAstruct:
1158  // BLOCK GAMBIT
1159  // 1 <high_scale> # Input scale of (upper) boundary conditions, e.g. GUT scale
1160 
1161  // Need to check if this information exists:
1162  SLHAstruct::const_iterator block = input_slha.find("GAMBIT");
1163  std::vector<std::string> key(1, "1");
1164  if(block == input_slha.end() or block->find(key) == block->end())
1165  {
1166  // Big problem
1167  std::ostringstream errmsg;
1168  errmsg << "Error constructing Spectrum object from a pre-existing SLHAstruct! " << endl
1169  << "The supplied SLHAstruct did not have the special GAMBIT block added. " << endl
1170  << "This block carries extra information from the spectrum generator " << endl
1171  << "that is usually thrown away, but which is needed for properly creating" << endl
1172  << "a Spectrum object. In whatever module function created the SLHAstruct " << endl
1173  << "that you want to use, please add code that adds the following " << endl
1174  << "information to the SLHAstruct (SLHAea::Coll): " << endl
1175  << " BLOCK GAMBIT " << endl
1176  << " 1 <high_scale> # Input scale of (upper) boundary conditions, e.g. GUT scale\n";
1177  SpecBit_error().raise(LOCAL_INFO,errmsg.str());
1178  }
1179 
1180  // OK the GAMBIT block exists, add the data to the MSSM SubSpectrum object.
1181  result.get_HE().set_override(Par::mass1,SLHAea::to<double>(input_slha.at("GAMBIT").at(1).at(1)), "high_scale", false);
1182  }
1183 
1185  void FH_MSSMMasses(fh_MSSMMassObs &result)
1186  {
1187  using namespace Pipes::FH_MSSMMasses;
1188 
1189  #ifdef SPECBIT_DEBUG
1190  cout << "****** calling FH_MSSMMasses ******" << endl;
1191  #endif
1192 
1193  // zero if minimal, non-zero if non-minimal flavour violation
1194  int nmfv;
1195 
1196  // MSf(s,t,g) MFV squark masses with indices
1197  // s = 1..2 sfermion index
1198  // t = 1..5 sfermion type nu,e,u,d,?
1199  // g = 1..3 generation index
1201 
1202  // USf(s1,s2,t,g) MFV squark mixing matrices with indices
1203  // s1 = 1..2 sfermion index (mass eigenstates)
1204  // s2 = 1..2 sfermion index (gauge eigenstates, L/R)
1205  // t = 1..5 sfermion type nu,e,u,d,?
1206  // g = 1..3 generation index
1208 
1209  // NMFV squark masses, with indices
1210  // a = 1..6 extended sfermion index
1211  // t = 1..5 sfermion type
1213 
1214  // NMFV squark mixing matrices, with indices
1215  // a1 = 1..6 extended sfermion index (mass eigenstates)
1216  // a2 = 1..6 extended sfermion index (gauge eigenstates)
1217  // t = 1..5 sftermion type nu,e,u,d,?
1219 
1220  // chargino masses
1221  Farray<fh_real, 1,2> MCha;
1222 
1223  // chargino mixing matrices (mass,gauge) eigenstates (2 x 2)
1226 
1227  // neutralino masses
1228  Farray<fh_real, 1,4> MNeu;
1229 
1230  // neutralino mixing matrices (mass,gauge) eigenstates (4 x 4)
1232 
1233  // correction to bottom Yukawa coupling
1234  fh_complex DeltaMB;
1235 
1236  // gluino mass
1237  fh_real MGl;
1238 
1239  // tree-level Higgs masses (Mh, MH, MA, MHpm)
1240  Farray<fh_real, 1,4> MHtree;
1241 
1242  // tree-level Higgs mixing parameters sin alpha
1243  fh_real SAtree;
1244 
1245  #ifdef SPECBIT_DEBUG
1246  cout << "****** calling FHGetPara ******" << endl;
1247  #endif
1248 
1249  int error = 1;
1250  BEreq::FHGetPara(error, nmfv, MSf, USf, MASf, UASf,
1251  MCha, UCha, VCha, MNeu, ZNeu,
1252  DeltaMB, MGl, MHtree, SAtree);
1253  if (error != 0)
1254  {
1255  std::ostringstream err;
1256  err << "BEreq::FHGetPara raised error flag: " << error << ".";
1257  invalid_point().raise(err.str());
1258  }
1259 
1260  fh_MSSMMassObs MassObs;
1261  for(int i = 0; i < 2; i++)
1262  for(int j = 0; j < 5; j++)
1263  for(int k = 0; k < 3; k++)
1264  MassObs.MSf[i][j][k] = MSf(i+1,j+1,k+1);
1265  for(int i = 0; i < 2; i++)
1266  for(int j = 0; j < 2; j++)
1267  for(int k = 0; k < 5; k++)
1268  for(int l = 0; l < 3; l++)
1269  MassObs.USf[i][j][k][l] = USf(i+1,j+1,k+1,l+1);
1270  for(int i = 0; i < 6; i++)
1271  for(int j = 0; j < 5; j++)
1272  MassObs.MASf[i][j] = MASf(i+1,j+1);
1273  for(int i = 0; i < 36; i++)
1274  for(int j = 0; j < 5; j++)
1275  MassObs.UASf[i][j] = UASf(i+1,j+1);
1276  for(int i = 0; i < 2; i++)
1277  MassObs.MCha[i] = MCha(i+1);
1278  for(int i = 0; i < 4; i++)
1279  {
1280  MassObs.UCha[i] = UCha(i+1);
1281  MassObs.VCha[i] = VCha(i+1);
1282  }
1283  for(int i = 0; i < 4; i++)
1284  MassObs.MNeu[i] = MNeu(i+1);
1285  for(int i = 0; i < 16; i++)
1286  MassObs.ZNeu[i] = ZNeu(i+1);
1287  MassObs.deltaMB = DeltaMB;
1288  MassObs.MGl = MGl;
1289  for(int i = 0; i < 4; i++)
1290  MassObs.MHtree[i] = MHtree(i+1);
1291  MassObs.SinAlphatree = SAtree;
1292 
1293  result = MassObs;
1294  }
1295 
1296 
1298  void FH_AllHiggsMasses(fh_HiggsMassObs &result)
1299  {
1300  using namespace Pipes::FH_AllHiggsMasses;
1301 
1302  #ifdef SPECBIT_DEBUG
1303  cout << "****** calling FH_HiggsMasses ******" << endl;
1304  #endif
1305 
1306  // Higgs mass with
1307  // 0 - m1 (Mh in rMSSM)
1308  // 1 - m2 (MH in rMSSM)
1309  // 2 - m3 (MA in rMSSM)
1310  // 3 - MHpm
1311  Farray<fh_real, 1,4> MHiggs;
1312  Farray<fh_real, 1,4> DeltaMHiggs;
1313 
1314  // sine of effective Higgs mixing angle, alpha_eff
1315  fh_complex SAeff;
1316  fh_complex DeltaSAeff;
1317 
1318  // matrix needed to rotate Higgs
1319  // mass matrix to diagonal form
1321  Farray<fh_complex, 1,3, 1,3> DeltaUHiggs;
1322 
1323  // matrix of Z-factors needed to combine
1324  // amplitudes involving on-shell Higgs
1326  Farray<fh_complex, 1,3, 1,3> DeltaZHiggs;
1327 
1328  #ifdef SPECBIT_DEBUG
1329  cout << "****** calling FHHiggsCorr ******" << endl;
1330  #endif
1331 
1332  int error = 1;
1333  BEreq::FHHiggsCorr(error, MHiggs, SAeff, UHiggs, ZHiggs);
1334  if (error != 0)
1335  {
1336  std::ostringstream err;
1337  err << "BEreq::FHHiggsCorr raised error flag: " << error << ".";
1338  invalid_point().raise(err.str());
1339  }
1340 
1341  #ifdef SPECBIT_DEBUG
1342  cout << "****** calling FHUncertainties ******" << endl;
1343  #endif
1344 
1345  error = 1;
1346  BEreq::FHUncertainties(error, DeltaMHiggs, DeltaSAeff, DeltaUHiggs, DeltaZHiggs);
1347  if (error != 0)
1348  {
1349  std::ostringstream err;
1350  err << "BEreq::FHUncertainties raised error flag: " << error << ".";
1351  invalid_point().raise(err.str());
1352  }
1353 
1354  fh_HiggsMassObs HiggsMassObs;
1355  for(int i = 0; i < 4; i++)
1356  {
1357  HiggsMassObs.MH[i] = MHiggs(i+1);
1358  HiggsMassObs.deltaMH[i] = DeltaMHiggs(i+1);
1359  }
1360  HiggsMassObs.SinAlphaEff = SAeff;
1361  HiggsMassObs.deltaSinAlphaEff = DeltaSAeff;
1362  for(int i = 0; i < 3; i++)
1363  for(int j = 0; j < 3; j++)
1364  {
1365  HiggsMassObs.UH[i][j] = UHiggs(i+1,j+1);
1366  HiggsMassObs.deltaUH[i][j] = DeltaUHiggs(i+1,j+1);
1367  HiggsMassObs.ZH[i][j] = ZHiggs(i+1,j+1);
1368  HiggsMassObs.deltaZH[i][j] = DeltaZHiggs(i+1,j+1);
1369  }
1370 
1371  result = HiggsMassObs;
1372  }
1373 
1377 
1379  void FH_Couplings(fh_Couplings &result)
1380  {
1381  using namespace Pipes::FH_Couplings;
1382 
1383  #ifdef SPECBIT_DEBUG
1384  cout << "****** calling FH_Couplings ******" << endl;
1385  #endif
1386 
1387  // what to use for internal Higgs mixing
1388  // (ex. in couplings)
1389  // (default = 1)
1390  // 0 - no mixing
1391  // 1 - UHiggs
1392  // 2 - ZHiggs
1393  int uzint = 2;
1394  // what to use for external Higgs mixing
1395  // (ex. in decays)
1396  // (default = 2)
1397  // 0 - no mixing
1398  // 1 - UHiggs
1399  // 2 - ZHiggs
1400  int uzext = 2;
1401  // which effective bottom mass to use
1402  int mfeff = 1;
1403 
1404  #ifdef SPECBIT_DEBUG
1405  cout << "****** calling FHSelectUZ ******" << endl;
1406  #endif
1407 
1408  int error = 1;
1409  BEreq::FHSelectUZ(error, uzint, uzext, mfeff);
1410  if (error != 0)
1411  {
1412  std::ostringstream err;
1413  err << "BEreq::FHSelectUZ raised error flag: " << error << ".";
1414  invalid_point().raise(err.str());
1415  }
1416 
1417  Farray<fh_complex, 1,ncouplings> couplings; // MSSM Higgs couplings
1418  Farray<fh_complex, 1,ncouplingsms> couplings_sm; // SM Higgs couplings
1419  Farray<fh_real, 1,ngammas> gammas; // Higgs decay widths and BR's (MSSM)
1420  Farray<fh_real, 1,ngammasms> gammas_sm; // Higgs decay widths and BR's (SM)
1421  int fast = 1; // include off-diagonal fermion decays? (1 = no)
1422 
1423  #ifdef SPECBIT_DEBUG
1424  cout << "****** calling FHCouplings ******" << endl;
1425  #endif
1426 
1427  error = 1;
1428  BEreq::FHCouplings(error, couplings, couplings_sm,
1429  gammas, gammas_sm, fast);
1430  if (error != 0)
1431  {
1432  std::ostringstream err;
1433  err << "BEreq::FHCouplings raised error flag: " << error << ".";
1434  invalid_point().raise(err.str());
1435  }
1436 
1437  fh_Couplings Couplings;
1438  for(int i = 0; i < ncouplings; i++) Couplings.couplings[i] = couplings(i+1);
1439  for(int i = 0; i < ncouplingsms; i++) Couplings.couplings_sm[i] = couplings_sm(i+1);
1440  for(int i = 0; i < ngammas; i++) Couplings.gammas[i] = gammas(i+1);
1441  for(int i = 0; i < ngammasms; i++) Couplings.gammas_sm[i] = gammas_sm(i+1);
1442  Couplings.calculator = BEreq::FHCouplings.origin();
1443  Couplings.calculator_version = BEreq::FHCouplings.version();
1444 
1445  result = Couplings;
1446  }
1447 
1448 
1450  std::vector<str> get_invisibles(const SubSpectrum& spec)
1451  {
1452  // Get the lighter of the lightest neutralino and the lightest sneutrino
1453  std::pair<str,double> neutralino("~chi0_1", spec.get(Par::Pole_Mass,"~chi0",1));
1454  std::pair<str,double> sneutrino("~nu_1", spec.get(Par::Pole_Mass,"~nu",1));
1455  std::pair<str,double> lnp = (neutralino.second < sneutrino.second ? neutralino : sneutrino);
1456 
1457  // Work out if this is indeed the LSP, and if decays of at least one neutral higgs to it are kinematically possible.
1458  bool inv_lsp = spec.get(Par::Pole_Mass,"~chi+",1) > lnp.second and
1459  spec.get(Par::Pole_Mass,"~g") > lnp.second and
1460  spec.get(Par::Pole_Mass,"~d",1) > lnp.second and
1461  spec.get(Par::Pole_Mass,"~u",1) > lnp.second and
1462  spec.get(Par::Pole_Mass,"~e-",1) > lnp.second and
1463  (spec.get(Par::Pole_Mass,"h0",2) > 2.*lnp.second or
1464  spec.get(Par::Pole_Mass,"A0") > 2.*lnp.second);
1465 
1466  // Create a vector containing all invisible products of higgs decays.
1467  if (inv_lsp) return initVector<str>(lnp.first);
1468  return std::vector<str>();
1469  }
1470 
1473  {
1474  using namespace Pipes::MSSM_higgs_couplings_pwid;
1475 
1476  // Retrieve spectrum contents
1477  const SubSpectrum& spec = Dep::MSSM_spectrum->get_HE();
1478 
1479  // Set up neutral Higgses
1480  static const std::vector<str> sHneut = initVector<str>("h0_1", "h0_2", "A0");
1481 
1482  // Set the CP of the Higgs states. Note that this would need to be more sophisticated to deal with the complex MSSM!
1483  result.CP[0] = 1; //h0_1
1484  result.CP[1] = 1; //h0_2
1485  result.CP[2] = -1; //A0
1486 
1487  // Work out which SM values correspond to which SUSY Higgs
1488  int higgs = (SMlike_higgs_PDG_code(spec) == 25 ? 0 : 1);
1489  int other_higgs = (higgs == 0 ? 1 : 0);
1490 
1491  // Set the decays
1492  result.set_neutral_decays_SM(higgs, sHneut[higgs], *Dep::Reference_SM_Higgs_decay_rates);
1493  result.set_neutral_decays_SM(other_higgs, sHneut[other_higgs], *Dep::Reference_SM_other_Higgs_decay_rates);
1494  result.set_neutral_decays_SM(2, sHneut[2], *Dep::Reference_SM_A0_decay_rates);
1495  result.set_neutral_decays(0, sHneut[0], *Dep::Higgs_decay_rates);
1496  result.set_neutral_decays(1, sHneut[1], *Dep::h0_2_decay_rates);
1497  result.set_neutral_decays(2, sHneut[2], *Dep::A0_decay_rates);
1498  result.set_charged_decays(0, "H+", *Dep::H_plus_decay_rates);
1499  result.set_t_decays(*Dep::t_decay_rates);
1500 
1501  // Use them to compute effective couplings for all neutral higgses, except for hhZ.
1502  for (int i = 0; i < 3; i++)
1503  {
1504  result.C_WW2[i] = result.compute_effective_coupling(i, std::pair<int,int>(24, 0), std::pair<int,int>(-24, 0));
1505  result.C_ZZ2[i] = result.compute_effective_coupling(i, std::pair<int,int>(23, 0), std::pair<int,int>(23, 0));
1506  result.C_tt2[i] = result.compute_effective_coupling(i, std::pair<int,int>(6, 1), std::pair<int,int>(-6, 1));
1507  result.C_bb2[i] = result.compute_effective_coupling(i, std::pair<int,int>(5, 1), std::pair<int,int>(-5, 1));
1508  result.C_cc2[i] = result.compute_effective_coupling(i, std::pair<int,int>(4, 1), std::pair<int,int>(-4, 1));
1509  result.C_tautau2[i] = result.compute_effective_coupling(i, std::pair<int,int>(15, 1), std::pair<int,int>(-15, 1));
1510  result.C_gaga2[i] = result.compute_effective_coupling(i, std::pair<int,int>(22, 0), std::pair<int,int>(22, 0));
1511  result.C_gg2[i] = result.compute_effective_coupling(i, std::pair<int,int>(21, 0), std::pair<int,int>(21, 0));
1512  result.C_mumu2[i] = result.compute_effective_coupling(i, std::pair<int,int>(13, 1), std::pair<int,int>(-13, 1));
1513  result.C_Zga2[i] = result.compute_effective_coupling(i, std::pair<int,int>(23, 0), std::pair<int,int>(21, 0));
1514  result.C_ss2[i] = result.compute_effective_coupling(i, std::pair<int,int>(3, 1), std::pair<int,int>(-3, 1));
1515  }
1516 
1517  // Calculate hhZ effective couplings. Here we scale out the kinematic prefactor
1518  // of the decay width, assuming we are well above threshold if the channel is open.
1519  // If not, we simply assume SM couplings.
1520  const double mZ = Dep::MSSM_spectrum->get(Par::Pole_Mass,23,0);
1521  const double scaling = 8.*sqrt(2.)*pi/Dep::MSSM_spectrum->get_SMInputs().GF;
1522  for(int i = 0; i < 3; i++)
1523  for(int j = 0; j < 3; j++)
1524  {
1525  double mhi = spec.get(Par::Pole_Mass, sHneut[i]);
1526  double mhj = spec.get(Par::Pole_Mass, sHneut[j]);
1527  if (mhi > mhj + mZ and result.get_neutral_decays(i).has_channel(sHneut[j], "Z0"))
1528  {
1529  double gamma = result.get_neutral_decays(i).width_in_GeV*result.get_neutral_decays(i).BF(sHneut[j], "Z0");
1530  double k[2] = {(mhj + mZ)/mhi, (mhj - mZ)/mhi};
1531  for (int l = 0; l < 2; l++) k[l] = (1.0 - k[l]) * (1.0 + k[l]);
1532  double K = mhi*sqrt(k[0]*k[1]);
1533  result.C_hiZ2[i][j] = scaling / (K*K*K) * gamma;
1534  }
1535  else // If the channel is missing from the decays or kinematically disallowed, just return the SM result.
1536  {
1537  result.C_hiZ2[i][j] = 1.;
1538  }
1539  }
1540 
1541  // Work out which invisible decays are possible
1542  result.invisibles = get_invisibles(spec);
1543  }
1544 
1545 
1548  {
1549  using namespace Pipes::MSSM_higgs_couplings_FH;
1550 
1551  // Retrieve spectrum contents
1552  const SubSpectrum& spec = Dep::MSSM_spectrum->get_HE();
1553  const SMInputs& sminputs = Dep::MSSM_spectrum->get_SMInputs();
1554 
1555  // Set up neutral Higgses
1556  static const std::vector<str> sHneut = initVector<str>("h0_1", "h0_2", "A0");
1557 
1558  // Work out which SM values correspond to which SUSY Higgs
1559  int higgs = (SMlike_higgs_PDG_code(spec) == 25 ? 0 : 1);
1560  int other_higgs = (higgs == 0 ? 1 : 0);
1561 
1562  // Set the decays
1563  result.set_neutral_decays_SM(higgs, sHneut[higgs], *Dep::Reference_SM_Higgs_decay_rates);
1564  result.set_neutral_decays_SM(other_higgs, sHneut[other_higgs], *Dep::Reference_SM_other_Higgs_decay_rates);
1565  result.set_neutral_decays_SM(2, sHneut[2], *Dep::Reference_SM_A0_decay_rates);
1566  result.set_neutral_decays(0, sHneut[0], *Dep::Higgs_decay_rates);
1567  result.set_neutral_decays(1, sHneut[1], *Dep::h0_2_decay_rates);
1568  result.set_neutral_decays(2, sHneut[2], *Dep::A0_decay_rates);
1569  result.set_charged_decays(0, "H+", *Dep::H_plus_decay_rates);
1570  result.set_t_decays(*Dep::t_decay_rates);
1571 
1572  // Use the branching fractions to compute gluon, gamma/Z and second generation fermionic effective couplings
1573  for (int i = 0; i < 3; i++)
1574  {
1575  result.C_gg2[i] = result.compute_effective_coupling(i, std::pair<int,int>(21, 0), std::pair<int,int>(21, 0));
1576  result.C_gaga2[i] = result.compute_effective_coupling(i, std::pair<int,int>(22, 0), std::pair<int,int>(22, 0));
1577  result.C_Zga2[i] = result.compute_effective_coupling(i, std::pair<int,int>(23, 0), std::pair<int,int>(22, 0));
1578  result.C_mumu2[i] = result.compute_effective_coupling(i, std::pair<int,int>(13, 1), std::pair<int,int>(-13, 1));
1579  result.C_ss2[i] = result.compute_effective_coupling(i, std::pair<int,int>(3, 1), std::pair<int,int>(-3, 1));
1580  result.C_cc2[i] = result.compute_effective_coupling(i, std::pair<int,int>(4, 1), std::pair<int,int>(-4, 1));
1581  }
1582 
1583  // Use couplings to get effective third-generation couplings
1584  for(int i = 0; i < 3; i++)
1585  {
1586  // Compute effective couplings
1587  double g2_s[3], g2_p[3];
1588  for (int j = 0; j < 3; j++) // j=0,1,2 => tau, t, b
1589  {
1590  fh_complex fh_L = Dep::FH_Couplings_output->couplings[H0FF(i+1,j+2,3,3)-1];
1591  fh_complex fh_R = Dep::FH_Couplings_output->couplings[H0FF(i+1,j+2,3,3)+Roffset-1];
1592  fh_complex fh_SM_L = Dep::FH_Couplings_output->couplings_sm[H0FF(i+1,j+2,3,3)-1];
1593  fh_complex fh_SM_R = Dep::FH_Couplings_output->couplings_sm[H0FF(i+1,j+2,3,3)+RSMoffset-1];
1594  std::complex<double> L(fh_L.re,fh_L.im);
1595  std::complex<double> R(fh_R.re,fh_R.im);
1596  std::complex<double> SM_L(fh_SM_L.re,fh_SM_L.im);
1597  std::complex<double> SM_R(fh_SM_R.re,fh_SM_R.im);
1598  g2_s[j] = 0.25*pow(std::abs(R/SM_R + L/SM_L), 2.);
1599  g2_p[j] = 0.25*pow(std::abs(R/SM_R - L/SM_L), 2.);
1600  }
1601  result.C_tautau2[i] = g2_s[0] + g2_p[0];
1602  result.C_tt2[i] = g2_s[1] + g2_p[1];
1603  result.C_bb2[i] = g2_s[2] + g2_p[2];
1604 
1605  // Calculate CP of state
1606  if(g2_p[2] < 1e-10)
1607  result.CP[i] = 1;
1608  else if(g2_s[2] < 1e-10)
1609  result.CP[i] = -1;
1610  else
1611  result.CP[i] = 0.;
1612  }
1613 
1614  // Use couplings to get di-boson effective couplings
1615  for(int i = 0; i < 3; i++)
1616  {
1617  fh_complex c_gWW = Dep::FH_Couplings_output->couplings[H0VV(i+1,4)-1];
1618  fh_complex c_gWW_SM = Dep::FH_Couplings_output->couplings_sm[H0VV(i+1,4)-1];
1619  fh_complex c_gZZ = Dep::FH_Couplings_output->couplings[H0VV(i+1,3)-1];
1620  fh_complex c_gZZ_SM = Dep::FH_Couplings_output->couplings_sm[H0VV(i+1,3)-1];
1621  std::complex<double> WW(c_gWW.re,c_gWW.im);
1622  std::complex<double> WW_SM(c_gWW_SM.re,c_gWW_SM.im);
1623  std::complex<double> ZZ(c_gZZ.re,c_gZZ.im);
1624  std::complex<double> ZZ_SM(c_gZZ_SM.re,c_gZZ_SM.im);
1625  result.C_WW2[i] = pow(std::abs(WW/WW_SM), 2.);
1626  result.C_ZZ2[i] = pow(std::abs(ZZ/ZZ_SM), 2.);
1627  }
1628 
1629  // Use couplings to get hhZ effective couplings
1630  double norm = sminputs.GF*sqrt(2.)*sminputs.mZ*sminputs.mZ;
1631  for(int i = 0; i < 3; i++)
1632  for(int j = 0; j < 3; j++)
1633  {
1634  fh_complex c_gHV = Dep::FH_Couplings_output->couplings[H0HV(i+1,j+1)-1];
1635  double g2HV = c_gHV.re*c_gHV.re+c_gHV.im*c_gHV.im;
1636  result.C_hiZ2[i][j] = g2HV/norm;
1637  }
1638 
1639  // Work out which invisible decays are possible
1640  result.invisibles = get_invisibles(spec);
1641  }
1642 
1643 
1647 
1649  template<class Contents>
1650  void fill_map_from_subspectrum(std::map<std::string,double>&, const SubSpectrum&);
1651 
1653  void add_extra_MSSM_parameter_combinations(std::map<std::string,double>& specmap, const SubSpectrum& mssm)
1654  {
1655  double At = 0;
1656  double Ab = 0;
1657  const double Yt = mssm.get(Par::dimensionless, "Yu", 3, 3);
1658  const double Yb = mssm.get(Par::dimensionless, "Yd", 3, 3);
1659  if(std::abs(Yt) > 1e-12)
1660  {
1661  At = mssm.get(Par::mass1, "TYu", 3, 3) / Yt;
1662  }
1663  if(std::abs(Yb) > 1e-12)
1664  {
1665  Ab = mssm.get(Par::mass1, "TYd", 3, 3) / Yb;
1666  }
1667 
1668  const double MuSUSY = mssm.get(Par::mass1, "Mu");
1669  const double tb = mssm.get(Par::dimensionless, "tanbeta");
1670 
1671  specmap["Xt"] = At - MuSUSY / tb;
1672  specmap["Xb"] = Ab - MuSUSY * tb;
1674  str msf1, msf2;
1678  const static double tol = 0.5;
1679  const static bool pt_error = true;
1680  slhahelp::family_state_mix_matrix("~u", 3, msf1, msf2, mssm, tol,
1681  LOCAL_INFO, pt_error);
1682  specmap["mstop1"] = mssm.get(Par::Pole_Mass, msf1);
1683  specmap["mstop2"] = mssm.get(Par::Pole_Mass, msf2);
1684  slhahelp::family_state_mix_matrix("~d", 3, msf1, msf2, mssm, tol,
1685  LOCAL_INFO, pt_error);
1686  specmap["msbottom1"] = mssm.get(Par::Pole_Mass, msf1);
1687  specmap["msbottom2"] = mssm.get(Par::Pole_Mass, msf2);
1688  slhahelp::family_state_mix_matrix("~e-", 3, msf1, msf2, mssm, tol,
1689  LOCAL_INFO, pt_error);
1690  specmap["mstau1"] = mssm.get(Par::Pole_Mass, msf1);
1691  specmap["mstau2"] = mssm.get(Par::Pole_Mass, msf2);
1694  const str gs_suL = slhahelp::mass_es_from_gauge_es("~u_L", mssm, tol,
1695  LOCAL_INFO, pt_error);
1696  specmap["msupL"] = mssm.get(Par::Pole_Mass,gs_suL);
1697  const str gs_scL = slhahelp::mass_es_from_gauge_es("~c_L", mssm, tol,
1698  LOCAL_INFO, pt_error);
1699  specmap["mscharmL"] = mssm.get(Par::Pole_Mass,gs_scL);
1700  const str gs_sdL = slhahelp::mass_es_from_gauge_es("~d_L", mssm, tol,
1701  LOCAL_INFO, pt_error);
1702  specmap["msdownL"] = mssm.get(Par::Pole_Mass,gs_sdL);
1703  const str gs_ssL = slhahelp::mass_es_from_gauge_es("~s_L", mssm, tol,
1704  LOCAL_INFO, pt_error);
1705  specmap["msstrangeL"] = mssm.get(Par::Pole_Mass,gs_ssL);
1706  const str gs_suR = slhahelp::mass_es_from_gauge_es("~u_R", mssm, tol,
1707  LOCAL_INFO, pt_error);
1708  specmap["msupR"] = mssm.get(Par::Pole_Mass,gs_suR);
1709  const str gs_scR = slhahelp::mass_es_from_gauge_es("~c_R", mssm, tol,
1710  LOCAL_INFO, pt_error);
1711  specmap["mscharmR"] = mssm.get(Par::Pole_Mass,gs_scR);
1712  const str gs_sdR = slhahelp::mass_es_from_gauge_es("~d_R", mssm, tol,
1713  LOCAL_INFO, pt_error);
1714  specmap["msdownR"] = mssm.get(Par::Pole_Mass,gs_sdR);
1715  const str gs_ssR = slhahelp::mass_es_from_gauge_es("~s_R", mssm, tol,
1716  LOCAL_INFO, pt_error);
1717  specmap["msstrangeR"] = mssm.get(Par::Pole_Mass,gs_ssR);
1718  const str gs_seL = slhahelp::mass_es_from_gauge_es("~e_L", mssm, tol,
1719  LOCAL_INFO, pt_error);
1720  specmap["mselectronL"] = mssm.get(Par::Pole_Mass,gs_seL);
1721  const str gs_sMuL = slhahelp::mass_es_from_gauge_es("~mu_L", mssm, tol,
1722  LOCAL_INFO, pt_error);
1723  specmap["msmuonL"] = mssm.get(Par::Pole_Mass,gs_sMuL);
1724  const str gs_seR = slhahelp::mass_es_from_gauge_es("~e_R", mssm, tol,
1725  LOCAL_INFO, pt_error);
1726  specmap["mselectronR"] = mssm.get(Par::Pole_Mass,gs_seR);
1727  const str gs_sMuR = slhahelp::mass_es_from_gauge_es("~mu_R", mssm, tol,
1728  LOCAL_INFO, pt_error);
1729  specmap["msmuonR"] = mssm.get(Par::Pole_Mass,gs_sMuR);
1730  const str gs_snu1 = slhahelp::mass_es_from_gauge_es("~nu_e_L", mssm, tol,
1731  LOCAL_INFO, pt_error);
1732  specmap["msnue"] = mssm.get(Par::Pole_Mass,gs_snu1);
1733  const str gs_snu2 = slhahelp::mass_es_from_gauge_es("~nu_mu_L", mssm, tol,
1734  LOCAL_INFO, pt_error);
1735  specmap["msnumu"] = mssm.get(Par::Pole_Mass,gs_snu2);
1736  const str gs_snu3 = slhahelp::mass_es_from_gauge_es("~nu_tau_L", mssm, tol,
1737  LOCAL_INFO, pt_error);
1738  specmap["msnutau"] = mssm.get(Par::Pole_Mass,gs_snu3);
1739 
1740  }
1741 
1742  void get_MSSM_spectrum_as_map (std::map<std::string,double>& specmap)
1743  {
1744  namespace myPipe = Pipes::get_MSSM_spectrum_as_map;
1745  const Spectrum& mssmspec(*myPipe::Dep::MSSM_spectrum);
1746  fill_map_from_subspectrum<SpectrumContents::SM> (specmap, mssmspec.get_LE());
1747  fill_map_from_subspectrum<SpectrumContents::MSSM>(specmap, mssmspec.get_HE());
1748  add_extra_MSSM_parameter_combinations(specmap, mssmspec.get_HE());
1749  }
1750  void get_unimproved_MSSM_spectrum_as_map (std::map<std::string,double>& specmap)
1751  {
1752  namespace myPipe = Pipes::get_unimproved_MSSM_spectrum_as_map;
1754  fill_map_from_subspectrum<SpectrumContents::SM> (specmap, mssmspec.get_LE());
1755  fill_map_from_subspectrum<SpectrumContents::MSSM>(specmap, mssmspec.get_HE());
1756  add_extra_MSSM_parameter_combinations(specmap, mssmspec.get_HE());
1757  }
1759 
1761  template<class Contents>
1762  void fill_map_from_subspectrum(std::map<std::string,double>& specmap, const SubSpectrum& subspec)
1763  {
1765  static const Contents contents;
1766  static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters();
1767 
1768  for(std::vector<SpectrumParameter>::const_iterator it = required_parameters.begin();
1769  it != required_parameters.end(); ++it)
1770  {
1771  const Par::Tags tag = it->tag();
1772  const std::string name = it->name();
1773  const std::vector<int> shape = it->shape();
1774 
1776 
1777  // Check scalar case
1778  if(shape.size()==1 and shape[0]==1)
1779  {
1780  std::ostringstream label;
1781  label << name <<" "<< Par::toString.at(tag);
1782  specmap[label.str()] = subspec.get(tag,name);
1783  // Check again ignoring overrides (if the value has an override defined)
1784  if(subspec.has(tag,name,overrides_only) and
1785  subspec.has(tag,name,ignore_overrides))
1786  {
1787  label << " (unimproved)";
1788  specmap[label.str()] = subspec.get(tag,name,ignore_overrides);
1789  }
1790  }
1791  // Check vector case
1792  else if(shape.size()==1 and shape[0]>1)
1793  {
1794  for(int i = 1; i<=shape[0]; ++i) {
1795  std::ostringstream label;
1796  label << name <<"_"<<i<<" "<< Par::toString.at(tag);
1797  specmap[label.str()] = subspec.get(tag,name,i);
1798  // Check again ignoring overrides
1799  if(subspec.has(tag,name,i,overrides_only) and
1800  subspec.has(tag,name,i,ignore_overrides))
1801  {
1802  label << " (unimproved)";
1803  specmap[label.str()] = subspec.get(tag,name,i,ignore_overrides);
1804  }
1805  }
1806  }
1807  // Check matrix case
1808  else if(shape.size()==2)
1809  {
1810  for(int i = 1; i<=shape[0]; ++i) {
1811  for(int j = 1; j<=shape[0]; ++j) {
1812  std::ostringstream label;
1813  label << name <<"_("<<i<<","<<j<<") "<<Par::toString.at(tag);
1814  specmap[label.str()] = subspec.get(tag,name,i,j);
1815  // Check again ignoring overrides
1816  if(subspec.has(tag,name,i,j,overrides_only) and
1817  subspec.has(tag,name,i,j,ignore_overrides))
1818  {
1819  label << " (unimproved)";
1820  specmap[label.str()] = subspec.get(tag,name,i,j,ignore_overrides);
1821  }
1822  }
1823  }
1824  }
1825  // Deal with all other cases
1826  else
1827  {
1828  // ERROR
1829  std::ostringstream errmsg;
1830  errmsg << "Error, invalid parameter received while converting SubSpectrum with contents \""<<contents.getName()<<"\" 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.";
1831  errmsg << "Problematic parameter was: "<< tag <<", " << name << ", shape="<< shape;
1832  SpecBit_error().forced_throw(LOCAL_INFO,errmsg.str());
1833  }
1834  }
1835  // add the scale!
1836  specmap["scale(Q)"] = subspec.GetScale();
1837  }
1838 
1840  {
1841  using namespace Pipes::FH_HiggsMass;
1842  //FH indices: 0=h0_1, 1=h0_2
1843  int i = 0;
1844  const SubSpectrum& spec = Dep::unimproved_MSSM_spectrum->get_HE();
1845  int higgs = SMlike_higgs_PDG_code(spec);
1846  if (higgs == 25) i = 0;
1847  else if (higgs == 35) i = 1;
1848  else SpecBit_error().raise(LOCAL_INFO, "Urecognised SM-like Higgs PDG code!");
1849  result.central = Dep::FH_HiggsMasses->MH[i];
1850  result.upper = Dep::FH_HiggsMasses->deltaMH[i];
1851  result.lower = Dep::FH_HiggsMasses->deltaMH[i];
1852  }
1853 
1855  {
1856  using namespace Pipes::FH_HeavyHiggsMasses;
1857  const int neutrals[2] = {25, 35};
1858  int i = -1;
1859  const SubSpectrum& spec = Dep::unimproved_MSSM_spectrum->get_HE();
1860  int higgs = SMlike_higgs_PDG_code(spec);
1861  if (higgs == neutrals[0]) i = 1;
1862  else if (higgs == neutrals[1]) i = 0;
1863  else SpecBit_error().raise(LOCAL_INFO, "Urecognised SM-like Higgs PDG code!");
1864  result.clear();
1865  result[neutrals[i]].central = Dep::FH_HiggsMasses->MH[i];
1866  result[neutrals[i]].upper = Dep::FH_HiggsMasses->deltaMH[i];
1867  result[neutrals[i]].lower = Dep::FH_HiggsMasses->deltaMH[i];
1868  result[36].central = Dep::FH_HiggsMasses->MH[2];
1869  result[36].upper = Dep::FH_HiggsMasses->deltaMH[2];
1870  result[36].lower = Dep::FH_HiggsMasses->deltaMH[2];
1871  result[37].central = Dep::FH_HiggsMasses->MH[3];
1872  result[37].upper = Dep::FH_HiggsMasses->deltaMH[3];
1873  result[37].lower = Dep::FH_HiggsMasses->deltaMH[3];
1874  }
1875 
1877  {
1878  using namespace Pipes::SHD_HiggsMass;
1879 
1880  const Spectrum& fullspectrum = *Dep::unimproved_MSSM_spectrum;
1881  SLHAea::Coll slhaea = fullspectrum.getSLHAea(1);
1882 
1883  #ifdef SPECBIT_DEBUG
1884  cout << "****** calling SHD_HiggsMass ******" << endl;
1885  #endif
1886 
1887  MList<MReal> parameterList = {
1888  SLHAea::to<double>(slhaea.at("HMIX").at(2).at(1)), // tanbeta
1889  SLHAea::to<double>(slhaea.at("MSOFT").at(1).at(1)), // M1
1890  SLHAea::to<double>(slhaea.at("MSOFT").at(2).at(1)), // M2
1891  SLHAea::to<double>(slhaea.at("MSOFT").at(3).at(1)), // M3
1892  SLHAea::to<double>(slhaea.at("HMIX").at(1).at(1)), // mu
1893  SLHAea::to<double>(slhaea.at("AU").at(3).at(2)), // At
1894  SLHAea::to<double>(slhaea.at("MSOFT").at(43).at(1)), // mQ3
1895  SLHAea::to<double>(slhaea.at("MSOFT").at(46).at(1)), // mU3
1896  SLHAea::to<double>(slhaea.at("MSOFT").at(49).at(1)), // mD3
1897  SLHAea::to<double>(slhaea.at("MSOFT").at(42).at(1)), // mQ2
1898  SLHAea::to<double>(slhaea.at("MSOFT").at(45).at(1)), // mU2
1899  SLHAea::to<double>(slhaea.at("MSOFT").at(48).at(1)), // mD2
1900  SLHAea::to<double>(slhaea.at("MSOFT").at(41).at(1)), // mQ1
1901  SLHAea::to<double>(slhaea.at("MSOFT").at(44).at(1)), // mU1
1902  SLHAea::to<double>(slhaea.at("MSOFT").at(47).at(1)), // mD1
1903  SLHAea::to<double>(slhaea.at("MSOFT").at(33).at(1)), // mL3
1904  SLHAea::to<double>(slhaea.at("MSOFT").at(36).at(1)), // mE3
1905  SLHAea::to<double>(slhaea.at("MSOFT").at(32).at(1)), // mL2
1906  SLHAea::to<double>(slhaea.at("MSOFT").at(35).at(1)), // mE2
1907  SLHAea::to<double>(slhaea.at("MSOFT").at(31).at(1)), // mL1
1908  SLHAea::to<double>(slhaea.at("MSOFT").at(34).at(1)), // mE1
1909  sqrt(SLHAea::to<double>(slhaea.at("HMIX").at(4).at(1))) // mA
1910  };
1911 
1912  MReal MHiggs = BEreq::SUSYHD_MHiggs(parameterList);
1913 
1914  #ifdef SPECBIT_DEBUG
1915  cout << "****** calling SHD_DeltaHiggsMass ******" << endl;
1916  #endif
1917 
1918  MReal DeltaMHiggs = BEreq::SUSYHD_DeltaMHiggs(parameterList);
1919 
1920  result.central = MHiggs;
1921  result.upper = DeltaMHiggs;
1922  result.lower = DeltaMHiggs;
1923 
1924  }
1925 
1926 
1928 
1929  } // end namespace SpecBit
1930 } // end namespace Gambit
1931 
void drop_SLHAs_if_requested(const safe_ptr< Options > &, const str &)
Helper function to drop SLHA files.
Definition: spectrum.cpp:429
double get(const Par::Tags partype, const std::string &mass) const
Definition: spectrum.cpp:249
Define overloadings of the stream operator for various containers.
Array class that matches the memory structure and functionality of arrays in Fortran codes Syntax: Fa...
Definition: util_types.hpp:306
void MSSM_higgs_couplings_pwid(HiggsCouplingsTable &result)
Put together the Higgs couplings for the MSSM, from partial widths only.
std::vector< double > family_state_mix_matrix(str type, int generation, str &mass_es1, str &mass_es2, const SubSpectrum &mssm, double tol, str context, bool pterror)
Get the family mixing matrix and corresponding mass eigenstates, then check for interfamily mixing...
void MSSM_higgs_couplings_FH(HiggsCouplingsTable &result)
Put together the Higgs couplings for the MSSM, using FeynHiggs.
Helper function to determine which Higgs is most SM-like.
START_CAPABILITY NMSSM_spectrum
double C_ZZ2[max_neutral_higgses]
Rollcall header for module SpecBit.
double BF(const std::vector< std::pair< int, int > > &) const
Retrieve branching fraction for decay to a given final state.
General small utility macros.
std::vector< str > invisibles
Particles that higgses can decay invisibly to.
double C_tt2[max_neutral_higgses]
Declarations of convenience (i.e.
MSSM derived version of SubSpectrum class.
void set_t_decays(const DecayTable::Entry &)
Assign decay entries to top.
void get_unimproved_MSSM_spectrum_as_map(std::map< std::string, double > &specmap)
void get_MSSM_spectrum_as_SLHAea_SLHA1(SLHAstruct &result)
Extract an SLHAea version of the spectrum contained in a Spectrum object, in SLHA1 format...
double C_cc2[max_neutral_higgses]
double CP[max_neutral_higgses]
CP of neutral higgses.
double C_WW2[max_neutral_higgses]
Effective couplings for neutral higgses.
void FH_HeavyHiggsMasses(map_int_triplet_dbl &result)
This struct contains all the strings we will use for the MSSM in the maps.
const DecayTable::Entry & get_neutral_decays(int) const
Retrieve decays of a specific neutral Higgs, by index.
const std::vector< str > getNames(const args &... keys) const
Retrieve values from key-value pairs in options node.
#define LOCAL_INFO
Definition: local_info.hpp:34
SLHAstruct getSLHAea(int) const
SLHAea object getter First constructs an SLHAea object from the SMINPUTS object, then adds the info f...
Definition: spectrum.cpp:409
GAMBIT native higgs coupling table class.
virtual double GetScale() const
Returns the renormalisation scale of parameters.
void fill_MSSM63_input_altnames(T &input, const std::map< str, safe_ptr< const double > > &Param)
Helpers for making simple spectra from SLHAea objects.
Flexiblesusy model header includes for SpecBit.
void set_charged_decays(int, const str &, const DecayTable::Entry &)
Assign decay entries to charged higgses.
void get_MSSM_spectrum_from_SLHAstruct(Spectrum &result)
Get an MSSMSpectrum object from an SLHAstruct Wraps it up in MSSMSimpleSpec; i.e. ...
void FH_AllHiggsMasses(fh_HiggsMassObs &result)
Higgs masses and mixings with theoretical uncertainties.
double C_ss2[max_neutral_higgses]
std::chrono::milliseconds ms
bool hasKey(const args &... keys) const
Getters for key/value pairs (which is all the options node should contain)
void set_override_vector(const Par::Tags, const double, const std::vector< str > &, const bool allow_new=false, const bool decouple=false)
Vector override functions.
TYPE getValue(const args &... keys) const
static const std::vector< str > pole_mass_strs_1_3
pole mass strings with 1 index and three entries
void FH_MSSMMasses(fh_MSSMMassObs &result)
FeynHiggs SUSY masses and mixings.
SLHAea::Coll SLHAstruct
Less confusing name for SLHAea container class.
void get_GUTMSSMB_spectrum(Spectrum &)
GAMBIT error class.
Definition: exceptions.hpp:136
double C_mumu2[max_neutral_higgses]
void setup_QedQcd(softsusy::QedQcd &oneset, const SMInputs &sminputs)
Non-Gambit helper functions.
Definition: SpecBit.cpp:47
int SMlike_higgs_PDG_code(const SubSpectrum &)
Determine which MSSM higgs is most SM-like.
static const std::vector< str > pole_mass_strs_1_2
pole mass strings with 1 index and two entries
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
double C_bb2[max_neutral_higgses]
void convert_NMSSM_to_SM(Spectrum &result)
const double pi
Eigen::Matrix< double, 3, 3 > fill_3x3_parameter_matrix(const std::string &rootname, const std::map< str, safe_ptr< const double > > &Param)
Helper function for setting 3x3 matrix-valued parameters.
double C_tautau2[max_neutral_higgses]
double compute_effective_coupling(int index, const T &p1, const T &p2)
Compute a neutral higgs effective coupling from the current two-body neutral higgs decays...
double width_in_GeV
Total particle width (in GeV)
void add_extra_MSSM_parameter_combinations(std::map< std::string, double > &specmap, const SubSpectrum &mssm)
Adds additional information from interesting combinations of MSSM parameters.
bool has_channel(const std::vector< std::pair< int, int > > &) const
Check if a given final state exists in this DecayTable::Entry.
static const std::vector< str > pole_mass_pred
void SHD_HiggsMass(triplet< double > &result)
bool has_neutralino_LSP(const Spectrum &result)
Check that the spectrum has a neutralino LSP.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
void fill_MSSM63_input(T &input, const std::map< str, safe_ptr< const double > > &Param)
Helper function for filling MSSM63-compatible input parameter objects Leaves out mHu2, mHd2, SignMu, (mA, mu) because we use two different parameterisations of these.
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
Spectrum Spectrum
std::vector< T > MList
Definition: util_types.hpp:651
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
void convert_E6MSSM_to_SM(Spectrum &result)
static const std::vector< str > pole_mass_strs_1_6
pole mass strings with 1 index and six entries
std::vector< T > initVector(std::vector< T > vector)
SubSpectrum & get_LE()
Standard getters Return references to internal data members.
Definition: spectrum.cpp:224
TYPE getValueOrDef(TYPE def, const args &... keys) const
void fill_map_from_subspectrum(std::map< std::string, double > &, const SubSpectrum &)
Convert MSSM type Spectrum object into a map, so it can be printed.
void get_MSSM_spectrum_as_SLHAea_SLHA2(SLHAstruct &result)
Extract an SLHAea version of the spectrum contained in a Spectrum object, in SLHA2 format...
Eigen::Matrix< double, 3, 3 > fill_3x3_symmetric_parameter_matrix(const std::string &rootname, const std::map< str, safe_ptr< const double > > &Param)
As above, but for symmetric input (i.e. 6 entries, assumed to be the upper triangle) ...
void convert_MSSM_to_SM(Spectrum &result)
Gambit module functions.
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
double C_Zga2[max_neutral_higgses]
void get_MSSM_spectrum_from_SLHAfile(Spectrum &result)
Get an MSSMSpectrum object from an SLHA file Wraps it up in MSSMSimpleSpec; i.e.
This class is used to wrap the QedQcd object used by SoftSUSY and FlexibleSUSY in a Gambit SubSpectru...
void FH_Couplings(fh_Couplings &result)
Call FH_Couplings from FeynHiggs and collect the output.
void set_neutral_decays(int, const str &, const DecayTable::Entry &)
Assign decay entries to neutral higgses.
double pow(const double &a)
Outputs a^i.
void set_override(const Par::Tags, const double, const str &, const bool allow_new=false, const bool decouple=false)
Parameter override functions.
double get(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool=SafeBool(true)) const
Definition: spec.hpp:72
std::vector< YAML::sdd > mc_info
Typedefs for making it easier to manipulate mass cut and mass ratio cut info.
Definition: spectrum.hpp:119
double C_gaga2[max_neutral_higgses]
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.
void get_MSSM_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
std::map< int, triplet< double > > map_int_triplet_dbl
Shorthand for an int-to-double triplet map.
Definition: util_types.hpp:127
str mass_es_from_gauge_es(str gauge_es, double &max_mixing, std::vector< double > &gauge_composition, const SubSpectrum &mssm)
indentifies the state with largest gauge_es content also fills largest max_mixing and full gauge_comp...
double C_gg2[max_neutral_higgses]
std::vector< str > get_invisibles(const SubSpectrum &spec)
Helper function to work out if the LSP is invisible, and if so, which particle it is...
TODO: see if we can use this one:
Definition: Analysis.hpp:33
A small wrapper object for &#39;options&#39; nodes.
virtual bool has(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
Getters/Setters etc.
double MReal
Definition: util_types.hpp:647
SubSpectrum & get_HE()
Definition: spectrum.cpp:225
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110
double C_hiZ2[max_neutral_higgses][max_neutral_higgses]
void get_SM_SubSpectrum_from_MSSM_Spectrum(const SubSpectrum *&result)
Retrieve SubSpectrum* to SM LE model from Spectrum object DEPENDENCY(MSSM_spectrum, Spectrum)
Container class for Standard Model input information (defined as in SLHA2)
Definition: sminputs.hpp:29
void FH_HiggsMass(triplet< double > &result)
void get_MSSM_spectrum_SPheno(Spectrum &spectrum)