gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
mssm_slhahelp.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
23 
27 
28 namespace Gambit
29 {
30 
31  namespace slhahelp
32  {
33 
36  const std::map<str, p_int_string> gauge_label_to_index_type = init_gauge_label_to_index_type();
37  const std::map<str, p_int_string> mass_label_to_index_type = init_mass_label_to_index_type();
38  const std::map<str, pair_string_ints> familystate_label = init_familystate_label();
39  const std::map<p_int_string, std::vector<str> > type_family_to_gauge_states = init_type_family_to_gauge_states();
40  const std::map<str,std::vector<str> > family_state_to_gauge_state = init_family_state_to_gauge_state();
41  const std::map<str,std::vector<str> > gauge_es_to_family_states = init_gauge_es_to_family_states() ;
42  const std::map<str,std::vector<str> > type_to_vec_of_mass_es = init_type_to_vec_of_mass_es();
43  const std::map<str,std::vector<str> > type_to_vec_of_gauge_es = init_type_to_vec_of_gauge_es();
45 
46  // FIXME: these two should be made members of the spectrum object itself
47  std::vector<double> get_Pole_Mixing_col(str type, int gauge_index, const SubSpectrum& mssm)
48  {
49  //extract info about indices for type using map
50  std::vector<str> mass_es_strs = type_to_vec_of_mass_es.at(type);
51  double col_length = mass_es_strs.size();
52  std::vector<double> mass_state_content(col_length);
53  //iterate over column in some way, e..g
54  for(std::vector<int>::size_type i = 1; i <= col_length; i++)
55  {
56  //Mix_{row, col}. Iterate through row index with column index fixed
57  mass_state_content[i - 1] = mssm.get(Par::Pole_Mixing,type, i, gauge_index);
58  }
59  return mass_state_content;
60  }
61  std::vector<double> get_Pole_Mixing_row(str type, int mass_index, const SubSpectrum& mssm)
62  {
63  std::vector<str> gauge_es_strs = type_to_vec_of_gauge_es.at(type);
64  double row_length = gauge_es_strs.size();
65  std::vector<double> gauge_state_content(row_length);
66  for(std::vector<int>::size_type i = 1; i <= row_length; i++)
67  {
69  gauge_state_content.at(i - 1) = mssm.get(Par::Pole_Mixing,type, mass_index, i);
70  }
71  return gauge_state_content;
72  }
73 
75  void add_MODSEL_disclaimer(SLHAstruct& slha, const str& object)
76  {
77  slha.push_front("# depend on which calculator you intend this object or file to be used with.");
78  slha.push_front("# Note that block MODSEL is not automatically emitted, as its contents");
79  slha.push_front("# This SLHA(ea) object was created from a GAMBIT "+object+" object.");
80  }
81 
83  void attempt_to_add_SLHA1_mixing(const str& block, SLHAstruct& slha, const str& type,
84  const SubSpectrum& spec, double tol, str& s1, str& s2, bool pterror)
85  {
86  if (slha.find(block) == slha.end())
87  {
88  std::vector<double> matmix = slhahelp::family_state_mix_matrix(type, 3, s1, s2, spec, tol, LOCAL_INFO, pterror);
89  SLHAea_add_matrix(slha, block, matmix, 2, 2);
90  }
91  else
92  {
93  std::map<str,str> family_to_3gen; // TODO: make const or something
94  family_to_3gen["~u"] = "~t";
95  family_to_3gen["~d"] = "~b";
96  family_to_3gen["~e-"] = "~tau";
97  s1 = slhahelp::mass_es_closest_to_family(family_to_3gen.at(type)+"_1", spec, tol, LOCAL_INFO, pterror);
98  s2 = slhahelp::mass_es_closest_to_family(family_to_3gen.at(type)+"_2", spec, tol, LOCAL_INFO, pterror);
99  }
100  }
101 
105  std::vector<double> get_mass_comp_for_gauge(str gauge_es,
106  const SubSpectrum& mssm)
107  {
109  p_int_string index_type = gauge_label_to_index_type.at(gauge_es);
110  str type = index_type.second;
111  int gauge_index = index_type.first;
112 
113  std::vector<double> mass_state_content =
114  get_Pole_Mixing_col(type, gauge_index, mssm);
115 
116  return mass_state_content;
117  }
118 
121  double get_mixing_element(str gauge_es, str mass_es, const SubSpectrum& mssm)
122  {
124  p_int_string mass_es_index_type = mass_label_to_index_type.at(mass_es);
125  p_int_string gauge_es_index_type = gauge_label_to_index_type.at(gauge_es);
126  int gauge_index = gauge_es_index_type.first;
127  int mass_index = mass_es_index_type.first;
129  str type = mass_es_index_type.second;
130  str type_gauge = gauge_es_index_type.second;
131  if(type!=type_gauge)
132  {
134  utils_error().raise(LOCAL_INFO, "function get_mixing_element "
135  "called with types for the gauge eigenstate and mass eigenstate that don't match.");
136  }
138  double admix = mssm.get(Par::Pole_Mixing,type, mass_index,
139  gauge_index);
140  return admix;
141  }
142 
147  std::vector<double> get_gauge_comp_for_mass(str mass_es, const SubSpectrum& mssm)
148  {
150  p_int_string index_type = mass_label_to_index_type.at(mass_es);
151  int mass_index = index_type.first;
152  str type = index_type.second;
153  //fill vector with mixings
154  std::vector<double> mass_state_content =
155  get_Pole_Mixing_row(type, mass_index, mssm);
156 
157  return mass_state_content;
158  }
159 
162  str mass_es_from_gauge_es(str gauge_es, double & max_mixing,
163  std::vector<double> & gauge_composition,
164  const SubSpectrum& mssm)
165  {
167  double temp_admix = 0.0;
169  max_mixing = 0;
171  str type = (gauge_label_to_index_type.at(gauge_es)).second;
172  str mass_es, temp_mass_es;
174  std::vector<str> mass_es_set = type_to_vec_of_mass_es.at(type);
175  typedef std::vector<str>::iterator iter;
176  for(iter it = mass_es_set.begin(); it != mass_es_set.end(); ++it){
177  temp_mass_es = *it;
178  temp_admix = get_mixing_element(gauge_es, temp_mass_es,
179  mssm);
180  gauge_composition.push_back(temp_admix);
181  //select largest
182  if(fabs(temp_admix) > fabs(max_mixing))
183  {
184  max_mixing = temp_admix;
185  mass_es = temp_mass_es;
186  }
187  } //end iteration over temp_mass_es
188 
189  return mass_es;
190  }
191 
196  str mass_es_from_gauge_es(str gauge_es, double & max_mixing,
197  const SubSpectrum& mssm)
198  {
199  std::vector<double> gauge_composition;
200  str mass_es = mass_es_from_gauge_es(gauge_es, max_mixing,
201  gauge_composition, mssm);
202  return mass_es;
203 
204  }
205 
211  std::vector<double> & gauge_composition,
212  const SubSpectrum& mssm)
213  {
214  double max_mixing = 0;
215  str mass_es = mass_es_from_gauge_es(gauge_es, max_mixing,
216  gauge_composition, mssm);
217 
218  return mass_es;
219  }
220 
226  const SubSpectrum& mssm)
227  {
228  double max_mixing = 0;
229  std::vector<double> gauge_composition;
230  str mass_es = mass_es_from_gauge_es(gauge_es, max_mixing,
231  gauge_composition, mssm);
232  return mass_es;
233  }
234 
236  str mass_es_from_gauge_es(str gauge_es, const SubSpectrum& mssm,
237  double tol, str context, bool pterror)
238  {
239  double max_mixing = 0;
240  std::vector<double> gauge_composition;
241  str mass_es = mass_es_from_gauge_es(gauge_es, max_mixing,
242  gauge_composition, mssm);
243  if((max_mixing*max_mixing) <= 1-tol)
244  {
245  const str errmsg = "Mass_es_from_gauge_es requested when mixing "
246  "away from closest gauge_es is greater than tol.";
247  str full_context = LOCAL_INFO + " called from " + context;
248  if (pterror)
249  {
250  invalid_point().raise(errmsg+" Raised at: "+full_context);
251  }
252  else
253  {
254  utils_error().raise(full_context, errmsg);
255  }
256  }
257  return mass_es;
258  }
259 
262  str gauge_es_from_mass_es(str mass_es, double & max_mixing,
263  std::vector<double> & mass_composition,
264  const SubSpectrum& mssm)
265  {
267  double temp_admix = 0.0;
269  max_mixing = 0;
271  str type = (mass_label_to_index_type.at(mass_es)).second;
272  str gauge_es, temp_gauge_es;
274  std::vector<str> gauge_es_vec = type_to_vec_of_gauge_es.at(type);
275  typedef std::vector<str>::iterator iter;
276  for(iter it = gauge_es_vec.begin(); it != gauge_es_vec.end(); ++it)
277  {
278  temp_gauge_es = *it;
279  temp_admix = get_mixing_element(temp_gauge_es, mass_es, mssm);
280  mass_composition.push_back(temp_admix);
281  //select largest
282  if(fabs(temp_admix) > fabs(max_mixing))
283  {
284  max_mixing = temp_admix;
285  gauge_es = temp_gauge_es;
286  }
287  } //end iteration over temp_mass_es
288 
289  //return string for closest gauge_es
290  return gauge_es;
291  }
292 
297  str gauge_es_from_mass_es(str mass_es, double & max_mixing,
298  const SubSpectrum& mssm)
299  {
300  std::vector<double> mass_composition;
301  str gauge_es = gauge_es_from_mass_es(mass_es, max_mixing,
302  mass_composition, mssm);
303  return gauge_es;
304 
305  }
306 
312  std::vector<double> & mass_composition,
313  const SubSpectrum& mssm)
314  {
315  double max_mixing;
316  str gauge_es = gauge_es_from_mass_es(mass_es, max_mixing,
317  mass_composition, mssm);
318 
319  return gauge_es;
320  }
321 
327  const SubSpectrum& mssm)
328  {
329  double max_mixing;
330  std::vector<double> mass_composition;
331  str gauge_es = gauge_es_from_mass_es(mass_es, max_mixing,
332  mass_composition, mssm);
333 
334  return gauge_es;
335  }
336 
338  str gauge_es_from_mass_es(str mass_es, const SubSpectrum& mssm,
339  double tol, str context, bool pterror)
340  {
341  double max_mixing;
342  std::vector<double> mass_composition;
343  str gauge_es = gauge_es_from_mass_es(mass_es, max_mixing,
344  mass_composition, mssm);
345  if((max_mixing*max_mixing) <= 1-tol)
346  {
347  const str errmsg = "Gauge_es from mass_es requested when mxing away "
348  "from closest mass_es is greater than tol";
349  str full_context = LOCAL_INFO + " called from " + context;
350  if (pterror)
351  {
352  invalid_point().raise(errmsg+" Raised at: "+full_context);
353  }
354  else
355  {
356  utils_error().raise(full_context, errmsg);
357  }
358  }
359  return gauge_es;
360  }
361 
372  int family,
373  const SubSpectrum& mssm)
374  {
377  p_int_string gen_type(family,type);
378  std::vector<str> gauge_states;
379  try { gauge_states = type_family_to_gauge_states.at(gen_type); }
380  catch (std::out_of_range&) { utils_error().raise(LOCAL_INFO, "Sfermion type or generation index not recognised; use type=~u,~d,~e-, gen=1,2,3."); }
381  str gauge_esL=gauge_states[0];
382  str gauge_esR=gauge_states[1];
383 
386  str mass_esL = mass_es_from_gauge_es(gauge_esL, mssm);
387  str mass_esR = mass_es_from_gauge_es(gauge_esR, mssm);
388 
389  sspair answer;
390  int mass_index_L = (mass_label_to_index_type.at(mass_esL)).first;
391  int mass_index_R = (mass_label_to_index_type.at(mass_esR)).first;
392  // order pair by mass
393  if(mass_index_L < mass_index_R)
394  answer = std::make_pair(mass_esL,mass_esR);
395  else answer = std::make_pair(mass_esR,mass_esL);
396 
397  return answer;
398  }
399 
404  const SubSpectrum& mssm)
405  {
406  std::vector<str> family_gauge_states;
407  try { family_gauge_states = family_state_to_gauge_state.at(familystate); }
408  catch (std::out_of_range&) { utils_error().raise(LOCAL_INFO, "Unrecognised family state. ('"+familystate+"' was requested)"); }
409  str gauge_esL = family_gauge_states[0];
410  str gauge_esR = family_gauge_states[1];
411 
412  // finds the mass_es with the largets mixing to
413  // passed gauge_es
414  str mass_esL = mass_es_from_gauge_es(gauge_esL, mssm);
415  str mass_esR = mass_es_from_gauge_es(gauge_esR, mssm);
416 
417  // extract mass order (1 or 2) from string via map
418  pair_string_ints type_family_massorder =
419  familystate_label.at(familystate);
420  pair_ints family_massorder = type_family_massorder.second;
421  int mass_order = family_massorder.second;
422  // if massorder is 1 choose select from masstateL and mass_esR the one
423  // with the lowest index else take highest
424  int massorderL = (mass_label_to_index_type.at(mass_esL)).first;
425  int massorderR = (mass_label_to_index_type.at(mass_esR)).first;
426  str answer;
427  if( (mass_order == 1 && massorderL < massorderR) ||
428  (mass_order == 2 && massorderL > massorderR) ) answer = mass_esL;
429  else answer = mass_esR;
430 
431  return answer;
432  }
433 
437  std::vector<double> get_gauge_comp_for_family_state(str familystate,
438  str & mass_es,
439  const SubSpectrum& mssm)
440  {
441  //get mass_es using one of our routines
442  mass_es = mass_es_closest_to_family(familystate, mssm);
444  int mass_index = (mass_label_to_index_type.at(mass_es)).first;
445  pair_string_ints state_info = familystate_label.at(familystate);
446  str type = state_info.first;
447  std::vector<double> gauge_es_content =
448  get_Pole_Mixing_row(type, mass_index,mssm);
449 
450  return gauge_es_content;
451 
452  }
453 
459  std::vector<double> & gauge_composition,
460  std::vector<double> & off_family_mixing,
461  const SubSpectrum& mssm)
462  {
463  //get mass_es using one of our routines
464  str mass_es = mass_es_closest_to_family(familystate, mssm);
466  std::vector<str> gauge_states;
467  try { gauge_states = family_state_to_gauge_state.at(familystate); }
468  catch (std::out_of_range&) { utils_error().raise(LOCAL_INFO, "Unrecognised family state. ('"+familystate+"' was requested)"); }
469  str gauge_state_L = gauge_states[0];
470  str gauge_state_R = gauge_states[1];
471 
472  p_int_string gauge_Lindex_type =
473  gauge_label_to_index_type.at(gauge_state_L);
474  unsigned int gauge_L_index = gauge_Lindex_type.first;
475  str type = gauge_Lindex_type.second;
476  unsigned int gauge_R_index
477  = (gauge_label_to_index_type.at(gauge_state_R)).first;
478  int mass_index = (mass_label_to_index_type.at(mass_es)).first;
479  std::vector<str> gauge_es_strs = type_to_vec_of_gauge_es.at(type);
480  double row_length = gauge_es_strs.size();
481  for(std::vector<int>::size_type i = 1; i <= row_length; i++)
482  {
483  double temp = mssm.get(Par::Pole_Mixing,type, mass_index, i);
484  if(i == gauge_L_index || i == gauge_R_index)
485  gauge_composition.push_back(temp);
486  else off_family_mixing.push_back(temp);
487  }
488 
489  return mass_es;
490 
491  }
492 
497  std::vector<double> & gauge_composition,
498  const SubSpectrum& mssm)
499  {
500  std::vector<double> off_family_mixing;
501  str mass_es = mass_es_closest_to_family(familystate, gauge_composition,
502  off_family_mixing, mssm);
503  return mass_es;
504 
505  }
506 
511  double & sqr_sum_mix,
512  const SubSpectrum& mssm)
513  {
514  std::vector<double> off_family_mixing;
515  std::vector<double> gauge_composition;
516  str mass_es = mass_es_closest_to_family(familystate, gauge_composition,
517  off_family_mixing, mssm);
518  sqr_sum_mix = gauge_composition[0] * gauge_composition[0];
519  sqr_sum_mix += gauge_composition[1] * gauge_composition[1];
520  return mass_es;
521 
522  }
523 
526  str mass_es_closest_to_family(str familystate, const SubSpectrum& mssm,
527  double tol, str context, bool pterror)
528  {
529  std::vector<double> off_family_mixing;
530  std::vector<double> gauge_composition;
531  str mass_es = mass_es_closest_to_family(familystate, gauge_composition,
532  off_family_mixing, mssm);
533  double sqr_sum_mix = gauge_composition[0] * gauge_composition[0];
534  sqr_sum_mix += gauge_composition[1] * gauge_composition[1];
535  if(sqr_sum_mix <= 1-tol)
536  {
537  const str errmsg = "Mass_es_closest_to_family requested when family "
538  "mixing away from closest mass_es is greater than tol";
539  str full_context = LOCAL_INFO + " called from " + context;
540  if (pterror)
541  {
542  invalid_point().raise(errmsg+" Raised at: "+full_context);
543  }
544  else
545  {
546  utils_error().raise(full_context, errmsg);
547  }
548  }
549 
550  return mass_es;
551 
552  }
553 
555  std::vector<double> family_state_mix_matrix(str type /*"~u", "~d" or "~e-"*/, int generation,
556  str & mass_es1, str & mass_es2, const SubSpectrum& mssm,
557  double tol, str context, bool pterror)
558  {
559  std::vector<double> m = family_state_mix_matrix(type, generation, mass_es1, mass_es2, mssm);
560  if (m[0]*m[0] + m[1]*m[1] < 1-tol || m[2]*m[2] + m[3]*m[3] < 1-tol)
561  {
562  const str errmsg = "Too much interfamily mixing to safely determine "
563  "intrafamily mixing matrix.";
564  str full_context = LOCAL_INFO + " called from " + context;
565  if (pterror)
566  {
567  invalid_point().raise(errmsg+" Raised at: "+full_context);
568  }
569  else
570  {
571  utils_error().raise(full_context, errmsg);
572  }
573  }
574  return m;
575  }
576 
585  std::vector<double> family_state_mix_matrix(str type,
586  int family,
587  str & mass_es1,
588  str & mass_es2,
589  const SubSpectrum& mssm)
590  {
592  sspair mass_ess = identify_mass_ess_for_family(type, family, mssm);
593  mass_es1 = mass_ess.first;
594  mass_es2 = mass_ess.second;
595 
598  p_int_string gen_type(family,type);
599  std::vector<str> gauge_states;
600  try { gauge_states = type_family_to_gauge_states.at(gen_type); }
601  catch (std::out_of_range&) { utils_error().raise(LOCAL_INFO, "Sfermion type or generation index not recognised; use type=~u,~d,~e-, gen=1,2,3."); }
602  str gauge_es_L=gauge_states[0];
603  str gauge_es_R=gauge_states[1];
606  p_int_string gauge_Lindex_type =
607  gauge_label_to_index_type.at(gauge_es_L);
608  unsigned int gauge_L_index = gauge_Lindex_type.first;
609  unsigned int gauge_R_index
610  = (gauge_label_to_index_type.at(gauge_es_R)).first;
611 
612  str type_L = gauge_Lindex_type.second;
613  int mass_index1 = (mass_label_to_index_type.at(mass_es1)).first;
614  int mass_index2 = (mass_label_to_index_type.at(mass_es2)).first;
615  std::vector<double> mix_row_1;
616  std::vector<double> mix_row_2;
617  std::vector<str> gauge_es_strs = type_to_vec_of_gauge_es.at(type);
618  double row_length = gauge_es_strs.size();
619  for(std::vector<int>::size_type i = 1; i <= row_length; i++)
620  {
621  double temp1 = mssm.get(Par::Pole_Mixing,type, mass_index1, i);
622  double temp2 = mssm.get(Par::Pole_Mixing,type, mass_index2, i);
623  if(i == gauge_L_index || i == gauge_R_index)
624  {
625  mix_row_1.push_back(temp1);
626  mix_row_2.push_back(temp2);
627  }
628  }
629 
631  mix_row_1.insert(mix_row_1.end(), mix_row_2.begin(), mix_row_2.end());
632 
633  return mix_row_1;
634  }
635 
640  str gauge_es,
641  str & mass_es,
642  const SubSpectrum& mssm)
643  {
644  pair_string_ints type_family_massorder;
645  try { type_family_massorder = familystate_label.at(familystate); }
646  catch (std::out_of_range&) { utils_error().raise(LOCAL_INFO, "Unrecognised family state."); }
647  str family_type = type_family_massorder.first;
648  p_int_string gauge_es_index_type = gauge_label_to_index_type.at(gauge_es);
649  int gauge_index = gauge_es_index_type.first;
651  str type_gauge = gauge_es_index_type.second;
652  if(family_type!=type_gauge)
653  {
654  utils_error().raise(LOCAL_INFO, "function get_gauge_admix_for_family_state "
655  "called with types for the family state and mass eigenstate that don't match.");
656  }
658  mass_es = mass_es_closest_to_family(familystate, mssm);
660  int mass_index = (mass_label_to_index_type.at(mass_es)).first;
661  double admix = mssm.get(Par::Pole_Mixing,type_gauge, mass_index,
662  gauge_index);
663  return admix;
664  }
665 
670  str family_state_closest_to_mass_es(str mass_es, double & sum_sq_mix,
671  std::vector<double> & mass_comp,
672  const SubSpectrum& mssm)
673  {
675  str gauge_es = gauge_es_from_mass_es(mass_es, mass_comp, mssm);
677  std::vector<str> family_states = gauge_es_to_family_states.at(gauge_es);
678  str family_state1 = family_states[0];
679  str family_state2 = family_states[1];
680  std::vector<str> gauge_states = family_state_to_gauge_state.at(family_state1);
681  str gauge_es_L = gauge_states[0];
682  str gauge_es_R = gauge_states[1];
683  str mass_es_other;
684  if(gauge_es == gauge_es_L)
685  mass_es_other = mass_es_from_gauge_es(gauge_es_R, mssm);
686  else mass_es_other = mass_es_from_gauge_es(gauge_es_L, mssm);
688  int mass_index = (mass_label_to_index_type.at(mass_es)).first;
689  int mass_index_other = (mass_label_to_index_type.at(mass_es_other)).first;
690  str fam_state;
693  if(mass_index < mass_index_other) fam_state = family_state1;
694  else fam_state = family_state2;
695 
696  //get gauge_indices to sum correct mixing elements
697  int gauge_index_L = (gauge_label_to_index_type.at(gauge_es_L)).first;
698  int gauge_index_R = (gauge_label_to_index_type.at(gauge_es_R)).first;
700  sum_sq_mix = mass_comp.at(gauge_index_L-1) * mass_comp.at(gauge_index_L-1);
701  sum_sq_mix += mass_comp.at(gauge_index_R-1) * mass_comp.at(gauge_index_R-1);
702 
703  return fam_state;
704  }
705 
710  str family_state_closest_to_mass_es(str mass_es, double & sum_sq_mix,
711  const SubSpectrum& mssm)
712  {
713  std::vector<double> mass_comp;
714  str fs = family_state_closest_to_mass_es(mass_es, sum_sq_mix,
715  mass_comp, mssm);
716  return fs;
717  }
718 
723  std::vector<double> & mass_comp,
724  const SubSpectrum& mssm)
725  {
726  double sum_sq_mix;
727  str fs = family_state_closest_to_mass_es(mass_es, sum_sq_mix, mass_comp,
728  mssm);
729  return fs;
730  }
731 
736  double tol, str context, bool pterror)
737  {
738  double sum_sq_mix;
739  std::vector<double> mass_comp;
740  str fs = family_state_closest_to_mass_es(mass_es, sum_sq_mix,
741  mass_comp, mssm);
742  if(sum_sq_mix <= 1-tol)
743  {
744  const str errmsg = "Family_state_closest_to_mass_es called when family "
745  "mixing away from closest mass_es is greater than tol.";
746  str full_context = LOCAL_INFO + " called from " + context;
747  if (pterror)
748  {
749  invalid_point().raise(errmsg+" Raised at: "+full_context);
750  }
751  else
752  {
753  utils_error().raise(full_context, errmsg);
754  }
755  }
756  return fs;
757  }
758 
760  // Here we assume that all SM input info comes from the SMINPUT object,
761  // and all low-E stuff (quark pole masses and the like) come from the LE subspectrum.
762  // In other words all those things should be added to the SLHAea object via
763  // different functions to this one. Here we add only MSSM information
764  // NOTE: check the below statement:
765  // Note that the SMINPUT object's dump-to-SLHAea function does not know how to discriminate
766  // between SLHA1 and SLHA2, but that doesn't matter, as the SM parameters defined in SLHA2
767  // just constitute additional blocks/block entries, not replacements for SLHA1 blocks. In the
768  // MSSM sector, this is not true, and we take care to write version-specific blocks here.
769  //
770  // slha_version - should be 1 or 2. Specifies whether to output closest-matching SLHA1 format
771  // entries, or to maintain SLHA2 as is used internally.
772  void add_MSSM_spectrum_to_SLHAea(const SubSpectrum& mssmspec, SLHAstruct& slha, int slha_version)
773  {
774  std::ostringstream comment;
775 
776  // Make sure to overwrite all entries if they exist already (from say a "hurriedly" copied SM subspectrum + unknown extra MSSM junk)
777 
778  //SPINFO block should be added separately.
779  // MINPAR block; some programs need tanbeta(mZ), so we should output it here if possible
780  SLHAea_check_block(slha, "MINPAR");
781  if(mssmspec.has(Par::dimensionless,"tanbeta(mZ)"))
782  {
783  SLHAea_add_from_subspec(slha,LOCAL_INFO,mssmspec,Par::dimensionless,"tanbeta(mZ)","MINPAR",3,"# tanbeta(mZ)^DRbar");
784  }
785  int sgnmu = sgn(mssmspec.get(Par::mass1,"Mu")); // Mu isn't at the input scale anymore, but sign(mu) doesn't change with running.
786  SLHAea_add(slha,"MINPAR",4,sgnmu,"# sign(mu)", true);
787 
788  // HMIX block
789  SLHAea_delete_block(slha, "HMIX");
790  SLHAea_add_block (slha, "HMIX", mssmspec.GetScale());
791  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::mass1,"Mu","HMIX",1,"mu DRbar");
792  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::dimensionless,"tanbeta","HMIX",2,"tan(beta) = vu/vd DRbar");
793  if (not mssmspec.has(Par::mass1,"vu")) utils_error().raise(LOCAL_INFO, "MSSM subspectrum does not contain vu!");
794  if (not mssmspec.has(Par::mass1,"vd")) utils_error().raise(LOCAL_INFO, "MSSM subspectrum does not contain vd!");
795  double vu = mssmspec.get(Par::mass1,"vu");
796  double vd = mssmspec.get(Par::mass1,"vd");
797  SLHAea_add(slha,"HMIX",3,sqrt(vu*vu + vd*vd),"v = sqrt(vd^2 + vu^2) DRbar", true);
798  SLHAea_add_from_subspec(slha,LOCAL_INFO,mssmspec,Par::mass2,"mA2","HMIX",4,"m^2_A (tree)");
799  SLHAea_add_from_subspec(slha,LOCAL_INFO,mssmspec,Par::mass2,"BMu","HMIX",101,"Bmu DRbar");
800  SLHAea_add(slha,"HMIX",102,vd,"vd DRbar", true);
801  SLHAea_add(slha,"HMIX",103,vu,"vu DRbar", true);
802  // GAUGE block
803  SLHAea_delete_block(slha, "GAUGE");
804  SLHAea_add_block (slha, "GAUGE", mssmspec.GetScale());
805  // Scale gY is in SU(5)/GUT normalisation internally; convert it to SM normalisation for SLHA output by multiplying by sqrt(3/5).
806  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::dimensionless,"g1","GAUGE",1,"g' = g1 = gY DRbar", true, sqrt(3./5.));
807  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::dimensionless,"g2","GAUGE",2,"g = g2 DRbar");
808  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::dimensionless,"g3","GAUGE",3,"g_s = g3 DRbar");
809 
810  const int pdg_codes[33] = {24,25,35,37,36,1000021,1000024,1000037,1000022,1000023,1000025,1000035,1000006,2000006,1000005,2000005,1000015,2000015,1000012,1000014,1000016,1000001,1000003,
811  2000001,2000003,1000011,1000013,2000011,2000013,1000002,1000004,2000002,2000004,};
812 
813  // Here we add the SLHA1-specific stuff, for backwards compatibility with backwards backends.
814  if (slha_version == 1)
815  {
816  const str slha1_sfermions[21] = {"~t_1", "~t_2", "~b_1", "~b_2", "~tau_1", "~tau_2",
817  "~nu_e_L", "~nu_mu_L", "~nu_tau_L",
818  "~d_L", "~s_L", "~d_R", "~s_R", "~e_L", "~mu_L",
819  "~e_R", "~mu_R", "~u_L", "~c_L", "~u_R", "~c_R"};
820  str slha2_sfermions[21];
821 
822  // STOPMIX, SBOTMIX and STAUMIX blocks
823  slhahelp::attempt_to_add_SLHA1_mixing("STOPMIX", slha, "~u", mssmspec, 1.0, slha2_sfermions[0], slha2_sfermions[1], true);
824  slhahelp::attempt_to_add_SLHA1_mixing("SBOTMIX", slha, "~d", mssmspec, 1.0, slha2_sfermions[2], slha2_sfermions[3], true);
825  slhahelp::attempt_to_add_SLHA1_mixing("STAUMIX", slha, "~e-", mssmspec, 1.0, slha2_sfermions[4], slha2_sfermions[5], true);
826 
827  // MASS block. Do everything except sfermions the same way as SLHA2.
828  for(int i=0;i<12;i++)
829  {
830  str comment(Models::ParticleDB().long_name(pdg_codes[i], 0));
831  SLHAea_add_from_subspec(slha, LOCAL_INFO, mssmspec, Par::Pole_Mass, std::pair<int, int>(pdg_codes[i],0), "MASS", comment);
832  }
833  for(int i=0;i<21;i++)
834  {
835  if (i > 5)
836  {
837  double max_mixing; // Don't actually care about this; we're going to SLHA1 whether it is a good approximation or not.
838  slha2_sfermions[i] = slhahelp::mass_es_from_gauge_es(slha1_sfermions[i], max_mixing, mssmspec);
839  }
840  SLHAea_add(slha, "MASS", pdg_codes[i+12], mssmspec.get(Par::Pole_Mass, slha2_sfermions[i]), slha1_sfermions[i], true);
841  }
842  }
843  else if (slha_version == 2)
844  {
845  // MASS block
846  for(int i=0;i<33;i++)
847  {
848  str comment(Models::ParticleDB().long_name(pdg_codes[i], 0));
849  SLHAea_add_from_subspec(slha, LOCAL_INFO, mssmspec, Par::Pole_Mass, std::pair<int, int>(pdg_codes[i],0), "MASS", comment);
850  }
851 
852  // USQMIX, DSQMIX, SELMIX
853  sspair S[3] = {sspair("USQMIX","~u"), sspair("DSQMIX","~d"), sspair("SELMIX","~e-")};
854  for (int k=0;k<3;k++)
855  {
856  SLHAea_delete_block(slha, S[k].first);
857  SLHAea_add_block (slha, S[k].first, mssmspec.GetScale());
858  for(int i=1;i<7;i++) for(int j=1;j<7;j++)
859  {
860  comment.str(""); comment << S[k].second << "-type sfermion mixing (" << i << "," << j << ")";
861  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec, Par::Pole_Mixing, S[k].second, i, j, S[k].first, i, j, comment.str());
862  }
863  }
864 
865  // SNUMIX block
866  sspair V("SNUMIX","~nu");
867  SLHAea_delete_block(slha, V.first);
868  SLHAea_add_block (slha, V.first, mssmspec.GetScale());
869  for(int i=1;i<4;i++) for(int j=1;j<4;j++)
870  {
871  comment.str(""); comment << V.second << " mixing matrix (" << i << "," << j << ")";
872  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec, Par::Pole_Mixing, V.second, i, j, V.first, i, j, comment.str());
873  }
874 
875  }
876  else
877  {
878  utils_error().raise(LOCAL_INFO, "Unrecognised version of SLHA standard; only SLHA1 and SLHA2 are permitted.");
879  }
880 
881  // MSOFT block (SLHA1 and SLHA2) plus MSL2, MSE2, MSQ2, MSU2 and MSD2 blocks (SLHA2 only)
882  if(not SLHAea_block_exists(slha,"MSOFT"))
883  {
884  SLHAea_add_block(slha, "MSOFT", mssmspec.GetScale());
885  }
886  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::mass1,"M1","MSOFT",1,"bino mass parameter M1");
887  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::mass1,"M2","MSOFT",2,"wino mass parameter M2");
888  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::mass1,"M3","MSOFT",3,"gluino mass parameter M3");
889  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::mass2,"mHd2","MSOFT",21,"d-type Higgs mass parameter mHd2");
890  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec,Par::mass2,"mHu2","MSOFT",22,"u-type Higgs mass parameter mHu2");
891  sspair M[5] = {sspair("MSL2","ml2"), sspair("MSE2","me2"), sspair("MSQ2","mq2"), sspair("MSU2","mu2"), sspair("MSD2","md2")};
892  for (int k=0;k<5;k++)
893  {
894  SLHAea_delete_block(slha, M[k].first);
895  if (slha_version == 2) SLHAea_add_block(slha, M[k].first, mssmspec.GetScale());
896  for(int i=1;i<4;i++) for(int j=1;j<4;j++)
897  {
898  comment.str(""); comment << M[k].second << "(" << i << "," << j << ")";
899  if (slha_version == 2) SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec, Par::mass2, M[k].second, i, j, M[k].first, i, j, comment.str());
900  if (i== j)
901  {
902  double entry = mssmspec.get(Par::mass2, M[k].second, i, j);
903  SLHAea_add(slha,"MSOFT",30+3*k+i+(k>1?4:0),sgn(entry)*sqrt(std::abs(entry)),"sqrt("+comment.str()+")", true);
904  }
905  }
906  }
907 
908  // Yukawa and trilinear blocks. YU, YD and YE, plus [YU, YD and YE; SLHA1 only], or [TU, TD and TE; SLHA2 only].
909  sspair A[3] = {sspair("AU","Au"), sspair("AD","Ad"), sspair("AE","Ae")};
910  sspair Y[3] = {sspair("YU","Yu"), sspair("YD","Yd"), sspair("YE","Ye")};
911  sspair T[3] = {sspair("TU","TYu"), sspair("TD","TYd"), sspair("TE","TYe")};
912  for (int k=0;k<3;k++)
913  {
914  // Erase these blocks if they already exist; we need them in the right format
915  SLHAea_delete_block(slha, Y[k].first);
916  SLHAea_delete_block(slha, A[k].first);
917  SLHAea_delete_block(slha, T[k].first);
918  // Now add them back
919  SLHAea_add_block(slha, Y[k].first, mssmspec.GetScale());
920  if (slha_version == 1) SLHAea_add_block(slha, A[k].first, mssmspec.GetScale());
921  if (slha_version == 2) SLHAea_add_block(slha, T[k].first, mssmspec.GetScale());
922  for(int i=1;i<4;i++)
923  {
924  // SLHA2 says to output only diagonal of Yukawa matrices, since we should be in a basis in which they are diagonal.
925  // SLHA1 says to give only the 3,3 element, but we'll give the whole diagonal anyway, codes will ignore them if not
926  // needed.
927  comment.str(""); comment << Y[k].second << "(" << i << "," << i << ")";
928  SLHAea_add_from_subspec(slha,LOCAL_INFO,mssmspec,Par::dimensionless,Y[k].second, i, i, Y[k].first, i, i, comment.str());
929 
930  if (slha_version == 1)
931  {
932  comment.str(""); comment << A[k].second << "(" << i << "," << i << ")";
933  double invTii = 1.0/mssmspec.get(Par::dimensionless,Y[k].second,i,i);
934  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec, Par::mass1, T[k].second, i, i, A[k].first, i, i, comment.str(), true, invTii); // last argument is a rescaling factor
935  }
936  else if (slha_version == 2)
937  {
938  for(int j=1;j<4;j++)
939  {
940  comment.str(""); comment << T[k].second << "(" << i << "," << j << ")";
941  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec, Par::mass1, T[k].second, i, j, T[k].first, i, j, comment.str());
942  }
943  }
944  else
945  {
946  std::ostringstream msg;
947  msg << "Tried to ouput SLHA data, but SLHA version requested was invalid (should be 1 or 2, but received "<<slha_version<<")";
948  utils_error().raise(LOCAL_INFO,msg.str());
949  }
950  }
951  }
952 
953  // ALPHA block
954  // if this exists already, delete it entirely
955  auto it = slha.find("ALPHA");
956  if(it!=slha.end()) slha.erase(it); // TODO: format of this call is confusing, perhaps write a wrapper for it.
957  // ...and now add it back
958  SLHAea_add_block(slha, "ALPHA", mssmspec.GetScale());
959  slha["ALPHA"][""] << asin(mssmspec.get(Par::Pole_Mixing, "h0", 2, 2)) << "# sin^-1(SCALARMIX(2,2))";
960 
961  // UMIX and VMIX blocks, plus some FlexibleSUSY-only extensions: PSEUDOSCALARMIX, SCALARMIX and CHARGEMIX.
962  sspair U[5] = {sspair("UMIX","~chi-"), sspair("VMIX","~chi+"), sspair("PSEUDOSCALARMIX","A0"), sspair("SCALARMIX","h0"), sspair("CHARGEMIX","H+")};
963  for (int k=0;k<5;k++)
964  {
965  SLHAea_delete_block(slha, U[k].first);
966  SLHAea_add_block(slha, U[k].first, mssmspec.GetScale());
967  for(int i=1;i<3;i++) for(int j=1;j<3;j++)
968  {
969  comment.str(""); comment << U[k].second << " mixing matrix (" << i << "," << j << ")";
970  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec, Par::Pole_Mixing, U[k].second, i, j, U[k].first, i, j, comment.str());
971  }
972  }
973 
974  // NMIX block
975  sspair N("NMIX","~chi0");
976  SLHAea_delete_block(slha, N.first);
977  SLHAea_add_block(slha, N.first, mssmspec.GetScale());
978  for(int i=1;i<5;i++) for(int j=1;j<5;j++)
979  {
980  comment.str(""); comment << N.second << " mixing matrix (" << i << "," << j << ")";
981  SLHAea_add_from_subspec(slha, LOCAL_INFO,mssmspec, Par::Pole_Mixing, N.second, i, j, N.first, i, j, comment.str());
982  }
983  }
984 
985  } // namespace slhahelp
986 
987 
988  // Definitions of methods for struct mass_es_pseudonyms
989 
991  void mass_es_pseudonyms::refill(const SubSpectrum& mssm, double tol, bool pt_error, bool debug)
992  {
993  filled = false;
994  fill(mssm, tol, pt_error, debug);
995  }
996 
998  void mass_es_pseudonyms::fill(const SubSpectrum& mssm, double tol, bool pt_error, bool debug)
999  {
1000  if(filled == true) return; // Don't refill unnecessarily
1001 
1002  fill_mass_es_psn_gauge (isdl, isdlbar, "~d_L", mssm, tol, pt_error, debug);
1003  fill_mass_es_psn_gauge (isul, isulbar, "~u_L", mssm, tol, pt_error, debug);
1004  fill_mass_es_psn_gauge (issl, isslbar, "~s_L", mssm, tol, pt_error, debug);
1005  fill_mass_es_psn_gauge (iscl, isclbar, "~c_L", mssm, tol, pt_error, debug);
1006  fill_mass_es_psn_family(isb1, isb1bar, "~b_1", mssm, tol, pt_error, debug);
1007  fill_mass_es_psn_family(ist1, ist1bar, "~t_1", mssm, tol, pt_error, debug);
1008 
1009  fill_mass_es_psn_gauge (isell, isellbar, "~e_L", mssm, tol, pt_error, debug);
1010  fill_mass_es_psn_gauge (ismul, ismulbar, "~mu_L", mssm, tol, pt_error, debug);
1011  fill_mass_es_psn_family(istau1, istau1bar, "~tau_1", mssm, tol, pt_error, debug);
1012 
1013  fill_mass_es_psn_gauge(isnel, isnelbar, "~nu_e_L", mssm, tol, pt_error, debug);
1014  fill_mass_es_psn_gauge(isnmul, isnmulbar, "~nu_mu_L", mssm, tol, pt_error, debug);
1015  fill_mass_es_psn_gauge(isntaul, isntaulbar, "~nu_tau_L", mssm, tol, pt_error, debug);
1016 
1017  fill_mass_es_psn_gauge (isdr, isdrbar, "~d_R", mssm, tol, pt_error, debug);
1018  fill_mass_es_psn_gauge (isur, isurbar, "~u_R", mssm, tol, pt_error, debug);
1019  fill_mass_es_psn_gauge (issr, issrbar, "~s_R", mssm, tol, pt_error, debug);
1020  fill_mass_es_psn_gauge (iscr, iscrbar, "~c_R", mssm, tol, pt_error, debug);
1021  fill_mass_es_psn_family(isb2, isb2bar, "~b_2", mssm, tol, pt_error, debug);
1022  fill_mass_es_psn_family(ist2, ist2bar, "~t_2", mssm, tol, pt_error, debug);
1023 
1024  fill_mass_es_psn_gauge (iselr, iselrbar, "~e_R", mssm, tol, pt_error, debug);
1025  fill_mass_es_psn_gauge (ismur, ismurbar, "~mu_R", mssm, tol, pt_error, debug);
1026  fill_mass_es_psn_family(istau2, istau2bar, "~tau_2", mssm, tol, pt_error, debug);
1027 
1028  filled = true;
1029 
1030  if (debug) debug_print(mssm);
1031 
1032  }
1033 
1035  void mass_es_pseudonyms::fill_mass_es_psn_gauge(str& is, str& isbar, str gauge_es, const SubSpectrum& mssm,
1036  double tol, bool pt_error_on_mixing_failure, bool debug)
1037  {
1038  double max_mix = 0;
1039  str mass_es = slhahelp::mass_es_from_gauge_es(gauge_es, max_mix, mssm);
1040 
1041  if((max_mix*max_mix) >= 1-tol)
1042  {
1043  is = mass_es;
1044  isbar = Models::ParticleDB().get_antiparticle(mass_es);
1045  gauge_family_eigenstates[is] = gauge_es;
1046  gauge_family_eigenstates[isbar] = gauge_es.substr(0,gauge_es.length()-2)+"bar"+gauge_es.substr(gauge_es.length()-2);
1047  mass_eigenstates[gauge_family_eigenstates[isbar]] = is;
1048  mass_eigenstates[gauge_family_eigenstates[isbar]] = isbar;
1049  }
1050  else
1051  {
1052  std::stringstream ss;
1053  ss << "MSSM mass(-squared) eigenstate " << mass_es << " is only " << max_mix*max_mix*100 << "% gauge eigenstate " << gauge_es << "." << endl
1054  << "The requested tolerance is " << tol*100 << " => too much sfermion mixing to assume that this is a pure gauge eigenstate.";
1055  if (pt_error_on_mixing_failure)
1056  {
1057  invalid_point().raise(ss.str());
1058  }
1059  else
1060  {
1061  utils_error().raise(LOCAL_INFO, ss.str());
1062  }
1063  }
1064 
1065  if (debug) debug_print_gauge(mssm, gauge_es, mass_es, max_mix);
1066  }
1067 
1069  void mass_es_pseudonyms::fill_mass_es_psn_family(str& is, str& isbar, str family_state, const SubSpectrum& mssm,
1070  double tol, bool pt_error_on_mixing_failure, bool debug)
1071  {
1074  std::vector<double> right_fam_gauge_comp;
1075  str mass_es = slhahelp::mass_es_closest_to_family(family_state, right_fam_gauge_comp, mssm);
1076 
1078  double mix_mag_sq = 0.0;
1079  for(auto i = right_fam_gauge_comp.begin(); i != right_fam_gauge_comp.end(); i++)
1080  {
1081  double mix = *i;
1082  mix_mag_sq += mix*mix;
1083  }
1084 
1085  if(mix_mag_sq > 1-tol)
1086  {
1087  is = mass_es;
1088  isbar = Models::ParticleDB().get_antiparticle(mass_es);
1089  gauge_family_eigenstates[is] = family_state;
1090  gauge_family_eigenstates[isbar] = family_state.substr(0,family_state.length()-2)+"bar"+family_state.substr(family_state.length()-2);
1091  mass_eigenstates[gauge_family_eigenstates[isbar]] = is;
1092  mass_eigenstates[gauge_family_eigenstates[isbar]] = isbar;
1093  }
1094  else
1095  {
1096  std::stringstream ss;
1097  ss << "MSSM mass(-squared) eigenstate " << mass_es << " is only " << mix_mag_sq*100 << "% the same family as "
1098  << "family state " << family_state << "." << endl << "The requested tolerance is " << tol*100
1099  << "% => too much inter-family sfermion mixing to assume that this is a pure family state.";
1100  if (pt_error_on_mixing_failure)
1101  {
1102  invalid_point().raise(ss.str());
1103  }
1104  else
1105  {
1106  utils_error().raise(LOCAL_INFO, ss.str());
1107  }
1108  }
1109 
1110  if (debug) debug_print_family(mssm, family_state, mass_es, mix_mag_sq, tol);
1111  }
1112 
1114  void mass_es_pseudonyms::debug_print(const SubSpectrum& mssm)
1115  {
1116  std::cout.precision(8);
1117  std::cout << "Dmix :" << std::endl;;
1118  for(int i = 1; i <=6; i++)
1119  {
1120  for(int j = 1; j <=6; j++)
1121  {
1122  std::cout << " " << i << j << " = "
1123  << std::scientific << std::setw(10)
1124  << mssm.get(Par::Pole_Mixing,"~d", i, j);
1125  }
1126  std::cout << std::endl;
1127  }
1128 
1129  std::cout << "Umix :" << std::endl;;
1130  for(int i = 1; i <=6; i++)
1131  {
1132  for(int j = 1; j <=6; j++)
1133  {
1134  std::cout << " " << i << j << " = "
1135  << mssm.get(Par::Pole_Mixing,"~u", i, j);
1136  }
1137  std::cout << std::endl;
1138  }
1139 
1140  std::cout << "Emix :" << std::endl;;
1141  for(int i = 1; i <=6; i++)
1142  {
1143  for(int j = 1; j <=6; j++)
1144  {
1145  std::cout << " " << i << j << " = "
1146  << mssm.get(Par::Pole_Mixing,"~e-", i, j);
1147  }
1148  std::cout << std::endl;
1149  }
1150 
1151  std::cout << "NUmix :" << std::endl;;
1152  for(int i = 1; i <=3; i++)
1153  {
1154  for(int j = 1; j <=3; j++)
1155  {
1156  std::cout << " " << i << j << " = "
1157  << mssm.get(Par::Pole_Mixing,"~nu", i, j);
1158  }
1159  std::cout << std::endl;
1160  }
1161 
1162  std::cout << "isdl = " << isdl << std::endl;
1163  std::cout << "isdlbar = " << isdlbar << std::endl;
1164  std::cout << "isdr = " << isdr << std::endl;
1165  std::cout << "isdrbar = " << isdrbar << std::endl;
1166 
1167  std::cout << "isul = " << isul << std::endl;
1168  std::cout << "isulbar = " << isulbar << std::endl;
1169  std::cout << "isur = " << isur << std::endl;
1170  std::cout << "isurbar = " << isurbar << std::endl;
1171 
1172  std::cout << "issl = " << issl << std::endl;
1173  std::cout << "isslbar = " << isslbar << std::endl;
1174  std::cout << "issr = " << issr << std::endl;
1175  std::cout << "issrbar = " << issrbar << std::endl;
1176 
1177  std::cout << "iscl = " << iscl << std::endl;
1178  std::cout << "isclbar = " << isclbar << std::endl;
1179  std::cout << "iscr = " << iscr << std::endl;
1180  std::cout << "iscrbar = " << iscrbar << std::endl;
1181 
1182  std::cout << "isb1 = " << isb1 << std::endl;
1183  std::cout << "isb1bar = " << isb1bar << std::endl;
1184  std::cout << "isb2 = " << isb2 << std::endl;
1185  std::cout << "isb2bar = " << isb2bar << std::endl;
1186 
1187  std::cout << "ist1 = " << ist1 << std::endl;
1188  std::cout << "ist1bar = " << ist1bar << std::endl;
1189  std::cout << "ist2 = " << ist2 << std::endl;
1190  std::cout << "ist2bar = " << ist2bar << std::endl;
1191 
1192  std::cout << "isell = " << isell << std::endl;
1193  std::cout << "isellbar = " << isellbar << std::endl;
1194  std::cout << "iselr = " << iselr << std::endl;
1195  std::cout << "iselrbar = " << iselrbar << std::endl;
1196 
1197  std::cout << "isnel = " << isnel << std::endl;
1198  std::cout << "isnelbar = " << isnelbar << std::endl;
1199 
1200  std::cout << "ismul = " << ismul << std::endl;
1201  std::cout << "ismulbar = " << ismulbar << std::endl;
1202  std::cout << "ismur = " << ismur << std::endl;
1203  std::cout << "ismurbar = " << ismurbar << std::endl;
1204 
1205  std::cout << "isnmull = " << isnmul << std::endl;
1206  std::cout << "isnmullbar = " << isnmulbar << std::endl;
1207 
1208  std::cout << "istau1 = " << istau1 << std::endl;
1209  std::cout << "istau1bar = " << istau1bar << std::endl;
1210  std::cout << "istau2 = " << istau2 << std::endl;
1211  std::cout << "istau2bar = " << istau2bar << std::endl;
1212 
1213  std::cout << "isntaul = " << isntaul << std::endl;
1214  std::cout << "isntaulbar = " << isntaulbar << std::endl;
1215 
1216  }
1217 
1219  void mass_es_pseudonyms::debug_print_gauge(const SubSpectrum& mssm, str& gauge_es, str& mass_es, double& max_mix)
1220  {
1221  std::cout << "******** Extra tests ********* " << std::endl;
1222  std::cout << "gauge_es = " << gauge_es << std::endl;
1223  std::cout << "mass_es = " << mass_es << std::endl;
1224 
1225  double max_mix_r = 0.0;
1226  str g_es = slhahelp::gauge_es_from_mass_es(mass_es, max_mix_r, mssm);
1227  std::cout << "g_es = " << g_es << std::endl;
1228  std::cout << "max_mix = " << max_mix<< std::endl;
1229  std::cout << "max_mix_r = " << max_mix_r << std::endl;
1230  if(g_es != gauge_es) std::cout << "g_s error! " << std::endl;
1231  if(max_mix_r != max_mix) std::cout << "g_s max_mix_r error! " << std::endl;
1232 
1233  str ges = slhahelp::gauge_es_from_mass_es(mass_es, mssm, 1e-3, LOCAL_INFO, false);
1234  std::cout << "ges = " << ges << std::endl;
1235  if(ges != gauge_es) std::cout << "ges error! " << std::endl;
1236  str mes = slhahelp::mass_es_from_gauge_es(gauge_es, mssm, 1e-3, LOCAL_INFO, false);
1237  std::cout << "mes = " << ges << std::endl;
1238  if(mes != mass_es) std::cout << "mes error! " << std::endl;
1239  }
1240 
1242  void mass_es_pseudonyms::debug_print_family(const SubSpectrum& mssm, str& family_state, str& mass_es, double& mix_mag_sq, double& tol)
1243  {
1244  std::cout << "******** Extra tests ********* " << std::endl;
1245  std::cout << "family_state = " << family_state <<std::endl;
1246  std::cout << "mass_es obtained from family_state = " << mass_es
1247  << std::endl;
1248  double sum_sq_mix;
1249  str fs = slhahelp::family_state_closest_to_mass_es(mass_es, sum_sq_mix,
1250  mssm);
1251  std::cout << "fs obtained from mass_es = " << fs << std::endl;
1252  std::cout << "sum_sq_mix = " << sum_sq_mix << std::endl;
1253  std::cout << "mix_mag_sq = " << mix_mag_sq << std::endl;
1254  if(fs != family_state) std::cout << "fs error! = " << std::endl;
1255  str f_s = slhahelp::family_state_closest_to_mass_es(mass_es, mssm, 1e-3, LOCAL_INFO, false);
1256  std::cout << "f_s obtained from mass_es = " << f_s << std::endl;
1257  if(f_s != family_state) std::cout << "f_s error! = " << std::endl;
1258 
1259  str m_es = slhahelp::mass_es_closest_to_family(family_state, mssm, tol, LOCAL_INFO, false);
1260  std::cout << "m_es = " << m_es << std::endl;
1261  if(m_es != mass_es) std::cout << "m_es error! = " << std::endl;
1262 
1263  std::cout << "******** Special family_state_mix_matrix tests ********"
1264  << std::endl;
1265  str mass_es1, mass_es2, type;
1266  str types[] = {"~u","~d", "~e-"};
1267  std::set<str> set_type = {types, types+3};
1268  std::set<str>::iterator it;
1269  for (it = set_type.begin(); it != set_type.end(); ++it)
1270  {
1271  type = *it;
1272  for(int gen = 1; gen <=3; gen++)
1273  {
1274  std::cout << "entering type = " << type << " and gen "
1275  << gen << std::endl;
1276  std::vector<double> f_mix_matrix = slhahelp::family_state_mix_matrix(type,gen, mass_es1, mass_es2, mssm);
1277  std::cout << "mass_es1 = " << mass_es1 << std::endl;
1278  std::cout << "mass_es2 = " << mass_es2 << std::endl;
1279  for(int i = 0; i<=3; i++)
1280  {
1281  std::cout << "f_mix_matrix[" << i << "] = "
1282  << f_mix_matrix[i] << std::endl;
1283  }
1284  double row1sq = f_mix_matrix[0] * f_mix_matrix[0];
1285  row1sq += f_mix_matrix[1] * f_mix_matrix[1];
1286  double row2sq = f_mix_matrix[2] * f_mix_matrix[2];
1287  row2sq += f_mix_matrix[3] * f_mix_matrix[3];
1288  std::cout << "row1sq = " << row1sq << " row2sq = "
1289  << row2sq << std::endl;
1290  }
1291  }
1292  }
1293 
1294 
1295 } //namespace gambit
std::pair< int, str > p_int_string
Typedefs for pairs that we will use in maps.
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...
std::map< str, std::vector< str > > init_family_state_to_gauge_state()
maps directly from family string to left_right gauge_pairs helps us reuse other routines that take st...
std::vector< double > get_gauge_comp_for_mass(str mass_es, const SubSpectrum &mssm)
returns vector representing composition of requested mass eigenstate in terms of the slha2 gauge eige...
EXPORT_SYMBOLS error & utils_error()
Utility errors.
void add_MODSEL_disclaimer(SLHAstruct &slha, const str &object)
Add a disclaimer about the absence of a MODSEL block in a generated SLHAea object.
void attempt_to_add_SLHA1_mixing(const str &block, SLHAstruct &slha, const str &type, const SubSpectrum &spec, double tol, str &s1, str &s2, bool pterror)
Simple helper function for for adding missing SLHA1 2x2 family mixing matrices to an SLHAea object...
const std::map< str, pair_string_ints > familystate_label
void SLHAea_add_block(SLHAstruct &, const str &name, const double scale=-1)
Add a new block to an SLHAea object, with our without a scale.
str family_state_closest_to_mass_es(str mass_es, double &sum_sq_mix, std::vector< double > &mass_comp, const SubSpectrum &mssm)
returns family state that best matches the given mass_es fills a double with the sum of the square mi...
std::pair< str, pair_ints > pair_string_ints
#define LOCAL_INFO
Definition: local_info.hpp:34
virtual double GetScale() const
Returns the renormalisation scale of parameters.
std::vector< double > get_mass_comp_for_gauge(str gauge_es, const SubSpectrum &mssm)
returns vector representing composition of requested gauge state in terms of the slha2 mass eigenstat...
void SLHAea_add(SLHAstruct &slha, const str &block, const int index, const double value, const str &comment="", const bool overwrite=false)
Add an entry to an SLHAea object (if overwrite=false, only if it doesn&#39;t already exist) ...
std::vector< double > get_Pole_Mixing_row(str type, int mass_index, const SubSpectrum &mssm)
Routines to help translate between SLHA2 sfermions and SLHA1 (or similar) sfermions.
std::pair< str, str > sspair
Shorthand for a pair of standard strings.
Definition: util_types.hpp:64
std::vector< double > get_Pole_Mixing_col(str type, int gauge_index, const SubSpectrum &mssm)
double get_mixing_element(str gauge_es, str mass_es, const SubSpectrum &mssm)
routine to return mass state admixure for given gauge state in the end this is a trival routine but m...
General small utility functions.
sspair identify_mass_ess_for_family(str type, int family, const SubSpectrum &mssm)
identify the two mass eigenstate corresponding to the approximate family states, e.g.
std::map< p_int_string, std::vector< str > > init_type_family_to_gauge_states()
map to obtain left_right gauge_pairs from state info helps us reuse other routiones with string argum...
Functions for triggering initialisation code.
const std::map< p_int_string, std::vector< str > > type_family_to_gauge_states
std::map< str, p_int_string > init_mass_label_to_index_type()
map from mass eigenstate strings to string, index pairs
const std::map< str, p_int_string > gauge_label_to_index_type
Known maps filled at initialisation.
str get_antiparticle(str) const
Get the matching anti-particle long name for a particle in the database, using the long name...
Definition: partmap.cpp:231
void add_MSSM_spectrum_to_SLHAea(const SubSpectrum &mssmspec, SLHAstruct &slha, int slha_version)
Add an entire MSSM spectrum to an SLHAea object.
SLHAea::Coll SLHAstruct
Less confusing name for SLHAea container class.
double get_gauge_admix_for_family_state(str familystate, str gauge_es, str &mass_es, const SubSpectrum &mssm)
returns admix of gauge eigenstate in the mass eigenstate closest to the given family state and stores...
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
std::pair< int, int > pair_ints
bool SLHAea_block_exists(SLHAstruct &slha, const str &block)
Check if a block exists in an SLHAea object.
const std::map< str, p_int_string > mass_label_to_index_type
partmap & ParticleDB()
Database accessor function.
Definition: partmap.cpp:32
const std::map< str, std::vector< str > > type_to_vec_of_mass_es
std::map< str, std::vector< str > > init_type_to_vec_of_mass_es()
map from string representing type (ie up-squarks, down-squarks or charged sleptons) to appropriate se...
const std::map< str, std::vector< str > > gauge_es_to_family_states
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
const std::map< str, std::vector< str > > type_to_vec_of_gauge_es
std::map< str, std::vector< str > > init_type_to_vec_of_gauge_es()
map from string representing type (ie up-squarks, down-squarks or charged sleptons) to appropriate se...
const std::map< str, std::vector< str > > family_state_to_gauge_state
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
START_MODEL M
void SLHAea_delete_block(SLHAstruct &slha, const std::string &block)
Delete an entire block from an SLHAea object, if it exists (actually just the first block matching th...
str mass_es_closest_to_family(str familystate, const SubSpectrum &mssm)
identify the mass eigenstate corresponding to family state takes string and returns only requested st...
std::map< str, pair_string_ints > init_familystate_label()
map to extract info from family state
invalid_point_exception & invalid_point()
Invalid point exceptions.
void SLHAea_add_matrix(SLHAstruct &slha, const str &block, const std::vector< T > &matrix, const int rows, const int cols, const str &comment="", const bool overwrite=false)
Add a whole matrix to an SLHAea object if it doesn&#39;t already exist.
str gauge_es_from_mass_es(str mass_es, double &max_mixing, std::vector< double > &mass_composition, const SubSpectrum &mssm)
identifies gauge_es with largest mass_es content also fills largest max_mixing and full mass_composit...
std::map< str, p_int_string > init_gauge_label_to_index_type()
map from gauge eigenstate strings to string, index pairs
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...
int sgn(T val)
void SLHAea_add_from_subspec(SLHAstruct &slha, const str local_info, const SubSpectrum &subspec, const Par::Tags partype, const std::pair< int, int > &pdg_pair, const str &block, const str &comment, const bool error_if_missing=true, const double rescale=1.0)
Add an entry from a subspectrum getter to an SLHAea object; SLHA index given by pdg code...
TODO: see if we can use this one:
Definition: Analysis.hpp:33
std::vector< double > get_gauge_comp_for_family_state(str familystate, str &mass_es, const SubSpectrum &mssm)
returns vector with composition of closest the mass eigenstate to give family state in terms of gauge...
virtual bool has(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
Getters/Setters etc.
bool SLHAea_check_block(SLHAstruct &slha, const str &block)
Check if a block exists in an SLHAea object, add it if not.
std::map< str, std::vector< str > > init_gauge_es_to_family_states()
maps directly from gauge_es string to familystates helps us reuse other routines that take string arg...