gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
SpecBit_examples.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
17 
18 #include <string>
19 #include <sstream>
20 #include <cmath>
21 
22 #include "gambit/Utils/gambit_module_headers.hpp"
25 
26 #include "SLHAea/slhaea.h"
27 #include <boost/preprocessor/tuple/to_seq.hpp>
28 #include <boost/preprocessor/seq/for_each.hpp>
29 #include <boost/preprocessor/seq/elem.hpp>
30 #include <boost/preprocessor/seq/for_each_product.hpp>
31 
32 // For the demo output only
34 
35 namespace Gambit
36 {
37 
38  namespace SpecBit
39  {
40 
41  using namespace LogTags;
42 
43  // Helper function to add error information to report
44  void add_error(std::ostringstream& out, const std::exception& e, const std::string& msg)
45  {
46  out << "------------------------------" << std::endl;
47  out << "TEST FAIL: " << msg << std::endl;
48  out << "Exception thrown was: "<<e.what()<<std::endl;
49  return;
50  }
51 
52  // Example showing reading of SLHA information from Spectrum object, using
53  // both SLHAea objects and direct Spectrum getter functions.
54  void exampleRead (bool &result)
55  {
56 
57  // Retreive pointer to Spectrum object, delivered by dependency resolver
58  // Module function asks for Spectrum* with capability MSSM_spectrum.
59  // i.e. has DEPENDENCY(MSSM_spectrum, Spectrum*)
60  namespace myPipe = Pipes::exampleRead;
61  /*TAG*/ Spectrum fullspectrum = *myPipe::Dep::MSSM_spectrum;
62  const SubSpectrum* spec = fullspectrum->get_HE(); // MSSMSpec SubSpectrum object
63  const SubSpectrum* SM = fullspectrum->get_LE(); // QedQcdWrapper SubSpectrum object
64 
65  std::ostringstream report; // Information about any problems encountered
66 
67  // Extract SLHAea object
68  // This copies the data out. Could possible change it to pass out a
69  // reference instead, or have another function to do that.
70  SLHAea::Coll slhaea = fullspectrum->getSLHAea(2);
71  // for testing, write this to file
72  std::ofstream out1;
73  out1.open("SpecBit/exampleRead_test.slha");
74  out1 << slhaea;
75  out1.close();
76 
77  // SLHAea::Coll slhaea = spec->getSLHAea(2); // The above is just a wrapper for this.
78 
79  // If this is a valid model point, return true and dump information, else false
80 
81  // SLHAea objects behave mostly like maps, but with special kinds of keys. For
82  // "at" and "operator[]", it does automatic conversion, but for "find" it does
83  // not, so we have to manually do it.
84  SLHAea::Block spinfo = slhaea.at("SPINFO");
85  //std::vector<std::string> k3(1, "3");
86  std::vector<std::string> k4(1, "4");
87 
88  // See if error code entries exist
89  //if(spinfo.find(k3) == spinfo.end() and spinfo.find(k4) == spinfo.end())
90  if(spinfo.find(k4) == spinfo.end())
91  {
92  std::cout << "Good spectrum found! Inspecting contents..." << std::endl;
93  std::cout << std::endl << slhaea << std::endl;
94 
95  // Write to file so we can check it
96  std::string filename(myPipe::runOptions->getValue<std::string>("output_slha_filename"));
97  spec->writeSLHAfile(2,filename);
98 
99  // ---------------------------------------------------------
100  // BEGIN DEMO OF SPECTRUM OBJECT AND PARTICLE DATABASE
101  // ---------------------------------------------------------
102 
103  //#define ECHO(COMMAND) cout << " " << STRINGIFY(COMMAND) << " = " << COMMAND << endl;
104  //Replacing with a version that deals with SLHAea access errors
105  #define ECHO(COMMAND) \
106  { \
107  try { \
108  cout << " " << STRINGIFY(COMMAND) << " = " << COMMAND << endl;\
109  } \
110  catch (const std::exception& e) \
111  { add_error(report,e,STRINGIFY(COMMAND)); } \
112  }
113 
114  /* ----------Test particle database access ---------------- */
115  #define PDB Models::ParticleDB()
116 
117  // First check out what is actually in the database
118  PDB.check_contents();
119 
120  // Demo a couple of particle name retrievals
121  cout<<endl;
122  cout<<"Demo retrieval of lightest neutralino info from particle database"<<endl;
123  cout<<"-----------------------------------------------------------------"<<endl;
124  ECHO( PDB.pdg_pair("~chi0_1") ) // Input long name, retrieve PDG code + context integer
125  ECHO( PDB.pdg_pair("~chi0",1) ) // Input short name + index, retrieve PDG code + context integer
126  ECHO( PDB.long_name("~chi0",1) ) // Input short name + index, retrieve long name
127  ECHO( PDB.long_name(std::make_pair(1000022,0)) ) // Input PDG code + context integer, retrieve long name
128  ECHO( PDB.long_name(1000022,0) ) // Input PDG code + context integer, retrieve long name
129  ECHO( PDB.short_name_pair("~chi0_1") ) // Input long name, retrieve short name + index
130  ECHO( PDB.short_name_pair(std::make_pair(1000022,0)) ) // Input PDG code plus context integer, retrieve short name + index
131  ECHO( PDB.short_name_pair(1000022,0) ) // Input PDG code plus context integer, retrieve short name + index
132 
133  cout<<endl;
134  cout<<"Demo retrieval of antiparticle names/codes using particle names/codes"<<endl;
135  cout<<"-----------------------------------------------------------------"<<endl;
136  // Check existence in various ways
137  ECHO( PDB.has_antiparticle("~chi0_1") )
138  ECHO( PDB.has_antiparticle("~chi0",1) )
139  ECHO( PDB.has_antiparticle(std::make_pair("~chi0",1)) )
140  ECHO( PDB.has_antiparticle(1000022,0) )
141  ECHO( PDB.has_antiparticle(std::make_pair(1000022,0)) )
142  ECHO( PDB.get_antiparticle("~chi0_1") ) // Input long name, retrieve antiparticle long name
143  ECHO( PDB.get_antiparticle("~chi0",1) ) // Input short name + index, retrieve antiparticle short name + index
144  ECHO( PDB.get_antiparticle(std::make_pair("~chi0",1)) ) // Input short name + index, retrieve antiparticle short name + index
145  ECHO( PDB.get_antiparticle(1000022,0) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
146  ECHO( PDB.get_antiparticle(std::make_pair(1000022,0)) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
147  ECHO( PDB.has_antiparticle("~chi+_1") )
148  ECHO( PDB.has_antiparticle("~chi+",1) )
149  ECHO( PDB.has_antiparticle(std::make_pair("~chi+",1)) )
150  ECHO( PDB.has_antiparticle(1000024,0) )
151  ECHO( PDB.has_antiparticle(std::make_pair(1000024,0)) )
152  ECHO( PDB.get_antiparticle("~chi+_1") ) // Input long name, retrieve antiparticle long name
153  ECHO( PDB.get_antiparticle("~chi+",1) ) // Input short name + index, retrieve antiparticle short name + index
154  ECHO( PDB.get_antiparticle(std::make_pair("~chi+",1)) ) // Input short name + index, retrieve antiparticle short name + index
155  ECHO( PDB.get_antiparticle(1000024,0) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
156  ECHO( PDB.get_antiparticle(std::make_pair(1000024,0)) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
157  ECHO( PDB.has_antiparticle("u_1") )
158  ECHO( PDB.has_antiparticle("u",1) )
159  ECHO( PDB.has_antiparticle(std::make_pair("u",1)) )
160  ECHO( PDB.has_antiparticle(2,0) )
161  ECHO( PDB.has_antiparticle(std::make_pair(2,0)) )
162  ECHO( PDB.get_antiparticle("u_1") ) // Input long name, retrieve antiparticle long name
163  ECHO( PDB.get_antiparticle("u",1) ) // Input short name + index, retrieve antiparticle short name + index
164  ECHO( PDB.get_antiparticle(std::make_pair("u",1)) ) // Input short name + index, retrieve antiparticle short name + index
165  ECHO( PDB.get_antiparticle(2,0) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
166  ECHO( PDB.get_antiparticle(std::make_pair(2,0)) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
167 
168 
169  cout<<endl;
170  cout<<"Demo retrieval when no short name exists"<<endl;
171  cout<<"-----------------------------------------------------------------"<<endl;
172  ECHO( PDB.pdg_pair("H+") )
173  //ECHO( PDB.pdg_pair("H+",1) ) // Error!
174  //ECHO( PDB.long_name("H+",1) ) // Error!
175  ECHO( PDB.long_name(std::make_pair(37,0)) )
176  ECHO( PDB.long_name(37,0) )
177  //ECHO( PDB.short_name_pair("H+") ) // Error!
178  //ECHO( PDB.short_name_pair(std::make_pair(37,0)) ) // Error!
179  //ECHO( PDB.short_name_pair(37,0) ) // Error!
180  //ECHO( PDB.short_name_pair(37) ) // Error!
181  cout<<endl;
182  /* ----------------- Pole masses --------------------------- */
183 
184  cout<<"Begin demo retrievals from Spectrum and SubSpectrum objects"<<endl;
185  cout<<"-----------------------------------------------------------------"<<endl;
186  cout<<endl;
187  cout<<"First, general methods for accessing different sorts of information."<<endl;
188  cout<<endl;
189  // At the moment it is only pole masses which have getters overloaded to use
190  // the particle database information. It is only the MASS block in the
191  // spectrum generator output SLHA files which use PDG numbers anyway, so I
192  // think this makes sense.
193  cout<<"Lightest neutral Higgs boson pole mass:"<<endl;
194  ECHO( fullspectrum->get_Pole_Mass( PDB.short_name_pair(25,0) ) )
195  ECHO( fullspectrum->get_Pole_Mass( PDB.long_name(25,0) ) )
196  ECHO( fullspectrum->get_Pole_Mass(25,0) )
197  ECHO( fullspectrum->get_Pole_Mass( PDB.pdg_pair("h0",1) ) )
198  ECHO( fullspectrum->get_Pole_Mass("h0",1) )
199  ECHO( fullspectrum->get_Pole_Mass("h0_1") )
200 
201  ECHO( spec->get_Pole_Mass( PDB.short_name_pair(25,0) ) )
202  ECHO( spec->get_Pole_Mass( PDB.long_name(25,0) ) )
203  ECHO( spec->get_Pole_Mass(25,0) )
204  ECHO( spec->get_Pole_Mass( PDB.pdg_pair("h0",1) ) )
205  ECHO( spec->get_Pole_Mass("h0",1) )
206  ECHO( spec->get_Pole_Mass("h0_1") )
207 
208  cout<<endl;
209  cout<<"Retrieval of Spectrum object contents, with"<<endl;
210  cout<<"correspondence to SLHAea object entries"<<endl;
211  cout<<"-----------------------------------------------------------------"<<endl;
212 
213  // MZ was a bad first example; it is empty unless you switch on the SM pole mass
214  // calculator for flexiblesusy. We do not yet pass any input value of MZ though
215  // Similar issues with other gauge boson masses. So don't use these yet or
216  // you'll get zero for all these masses.
217  cout<<endl;
218  cout<<"Gauge boson pole masses:"<<endl;
219  cout<<endl;
220  ECHO( fullspectrum->get_Pole_Mass("Z0") )
221  ECHO( SM->get_Pole_Mass("Z0") )
222  ECHO( slhaea.at("SMINPUTS").at(4).at(1) )
223  cout<<endl;
224  ECHO( fullspectrum->get_Pole_Mass("gamma") )
225  ECHO( SM->get_Pole_Mass("gamma") )
226  cout<<" ***Not in slha***"<<endl;
227  cout<<endl;
228  ECHO( fullspectrum->get_Pole_Mass("W+") )
229  ECHO( SM->get_Pole_Mass("W+") )
230  ECHO( slhaea.at("MASS").at(24).at(1) )
231  cout<<endl;
232  ECHO( fullspectrum->get_Pole_Mass("g") )
233  ECHO( SM->get_Pole_Mass("g") )
234  cout<<" ***Not in slha***"<<endl;
235  cout<<endl;
236  cout<<"Quark pole masses (actually the slha entries are MSbar except the top mass):"<<endl;
237  cout<<endl;
238  // I'm a little unclear on what the pole masses for the lighter quarks mean, since I thought
239  // that non-perturbative effects made definining them difficult... well anyway will have
240  // to ask Peter what is being computed here.
241 
242  //ECHO( spec->get_Pole_Mass("u",1) ) // i.e. up (mass eigenstate)
243  cout<<" ***u Pole mass not well defined***"<<endl;
244  ECHO( slhaea.at("SMINPUTS").at(22).at(1) ) // mu(2 GeV)^MS-bar, not pole mass
245  cout<<endl;
246  //ECHO( spec->get_Pole_Mass("u",2) ) // i.e. charm
247  cout<<" ***c Pole mass not well defined***"<<endl;
248  ECHO( slhaea.at("SMINPUTS").at(24).at(1) ) // mc(mc)^MS-bar, not pole mass
249  cout<<endl;
250  //ECHO( spec->get_Pole_Mass("u",3) ) // i.e. top
251  ECHO( fullspectrum->get_Pole_Mass("t") )
252  ECHO( SM->get_Pole_Mass("u",3) ) // i.e. top
253  ECHO( slhaea.at("SMINPUTS").at(6).at(1) )
254  cout<<endl;
255  //ECHO( spec->get_Pole_Mass("d",1) ) // i.e. down
256  cout<<" ***d Pole mass not well defined***"<<endl;
257  ECHO( slhaea.at("SMINPUTS").at(21).at(1) ) // md(2 GeV)^MS-bar, not pole mass
258  cout<<endl;
259  //ECHO( spec->get_Pole_Mass("d",2) ) // i.e. strange
260  cout<<" ***s Pole mass not well defined***"<<endl;
261  ECHO( slhaea.at("SMINPUTS").at(23).at(1) ) // ms(2 GeV)^MS-bar, not pole mass
262  cout<<endl;
263  //ECHO( spec->get_Pole_Mass("d",3) ) // i.e. bottom
264  ECHO( fullspectrum->get_Pole_Mass("b") )
265  ECHO( SM->get_Pole_Mass("d",3) ) // i.e. bottom
266  ECHO( slhaea.at("SMINPUTS").at(5).at(1) ) // mb(mb)^MS-bar, not pole mass.
267  cout<<endl;
268  cout<<"Charged fermions pole masses:"<<endl;
269  cout<<endl;
270  //ECHO( spec->get_Pole_Mass("e-",1) ) // i.e. electron
271  ECHO( fullspectrum->get_Pole_Mass("e-") )
272  ECHO( SM->get_Pole_Mass("e-",1) ) // i.e. electron
273  ECHO( slhaea.at("SMINPUTS").at(11).at(1) )
274  cout<<endl;
275  //ECHO( spec->get_Pole_Mass("e-",2) ) // i.e. muon
276  ECHO( fullspectrum->get_Pole_Mass("mu-") )
277  ECHO( SM->get_Pole_Mass("e-",2) ) // i.e. muon
278  ECHO( slhaea.at("SMINPUTS").at(13).at(1) )
279  cout<<endl;
280  //ECHO( spec->get_Pole_Mass("e-",3) ) // i.e. tau
281  ECHO( fullspectrum->get_Pole_Mass("tau-") )
282  ECHO( SM->get_Pole_Mass("e-",3) ) // i.e. tau
283  ECHO( slhaea.at("SMINPUTS").at(7).at(1) )
284  cout<<endl;
285  cout<<"Neutrinos pole masses:"<<endl;
286  cout<<endl;
287  // These will produce errors because currently no neutrino mass getters are hooked up
288  //ECHO( spec->get_Pole_Mass("nu",1) ) // Just mass ordered (if there is mixing)
289  ECHO( fullspectrum->get_Pole_Mass("nu",1) )
290  ECHO( SM->get_Pole_Mass("nu",1) ) // Just mass ordered (if there is mixing)
291  ECHO( slhaea.at("SMINPUTS").at(12).at(1) )
292  cout<<endl;
293  //ECHO( spec->get_Pole_Mass("nu",2) )
294  ECHO( fullspectrum->get_Pole_Mass("nu",2) )
295  ECHO( SM->get_Pole_Mass("nu",2) )
296  ECHO( slhaea.at("SMINPUTS").at(14).at(1) )
297  cout<<endl;
298  //ECHO( spec->get_Pole_Mass("nu",3) )
299  ECHO( fullspectrum->get_Pole_Mass("nu",3) )
300  ECHO( SM->get_Pole_Mass("nu",3) )
301  ECHO( slhaea.at("SMINPUTS").at(8).at(1) )
302  cout<<endl;
303  // Now for SUSY particles
304  cout<<endl;
305  cout<<"MSSM Higgs sector pole masses:"<<endl;
306  cout<<endl;
307  ECHO( fullspectrum->get_Pole_Mass("h0",1) )
308  ECHO( spec->get_Pole_Mass("h0",1) ) // Lightest neutral Higgs boson
309  ECHO( slhaea.at("MASS").at(25).at(1) )
310  cout<<endl;
311  ECHO( fullspectrum->get_Pole_Mass("h0",2) )
312  ECHO( spec->get_Pole_Mass("h0",2) ) // Heavy neutral Higgs boson
313  ECHO( slhaea.at("MASS").at(35).at(1) )
314  cout<<endl;
315  ECHO( fullspectrum->get_Pole_Mass("H+") )
316  ECHO( spec->get_Pole_Mass("H+") ) // Charged Higgs
317  ECHO( slhaea.at("MASS").at(37).at(1) )
318  cout<<endl;
319  ECHO( fullspectrum->get_Pole_Mass("A0") )
320  ECHO( spec->get_Pole_Mass("A0") ) // Pseudoscalar neutral Higgs
321  ECHO( slhaea.at("MASS").at(36).at(1) )
322  cout<<endl;
323 
324  // I'm going to use these nested functors to save lots of typing for the rest. It is just the
325  // same as the examples above, except that the PDG codes are retrieved from the particle database.
326  // The PDG code - string name correspondences are defined in 'Models/src/particle_database.cpp'
327 
328  struct get_polemass_functor
329  {
330  // Single mass
331  void operator()(const std::string& longname)
332  {
333  std::ostringstream echo1, echo2, echo3;
334  echo1 << " fullspectrum->get_Pole_Mass("<<longname<<") = ";
335  double value1 = fullspectrum->get_Pole_Mass(longname);
336  echo2 << " spec->get_Pole_Mass("<<longname<<") = ";
337  double value2 = spec->get_Pole_Mass(longname);
338  echo3 << " slhaea.at(\"MASS\").at("<<PDB.pdg_pair(longname).first<<").at(1) = ";
339  str value3 = slhaea.at("MASS").at( PDB.pdg_pair(longname).first ).at(1);
340  cout << echo1.str() << value1 << endl;
341  cout << echo2.str() << value2 << endl;
342  cout << echo3.str() << value3 << endl;
343  cout<<endl;
344  }
345  // Range of indexes masses
346  void operator()(const std::string& longname, int from, int to)
347  {
348  for(int i=from; i<=to; ++i)
349  {
350  std::ostringstream echo1;
351  std::ostringstream echo2;
352  echo1 << " spec->get_Pole_Mass("<<longname<<","<<i<<") = ";
353  double value1 = spec->get_Pole_Mass(longname,i);
354  echo2 << " slhaea.at(\"MASS\").at("<<PDB.pdg_pair(longname,i).first<<").at(1) = ";
355  str value2 = slhaea.at("MASS").at( PDB.pdg_pair(longname,i).first ).at(1);
356  cout << echo1.str() << value1 << endl;
357  cout << echo2.str() << value2 << endl;
358  cout<<endl;
359  }
360  }
361 
362  get_polemass_functor(std::ostringstream& report, /*TAG*/ Spectrum fullin, const SubSpectrum* specin, SLHAea::Coll& slhaeain)
363  : report(report)
364  , fullspectrum(fullin)
365  , spec(specin)
366  , slhaea(slhaeain)
367  {}
368 
369  private:
370  std::ostringstream& report;
371  /*TAG*/ Spectrum fullspectrum;
372  const SubSpectrum* spec;
373  SLHAea::Coll slhaea;
374  };
375 
376  get_polemass_functor get_polemass(report,fullspectrum,spec,slhaea);
377 
378  cout<<endl<<"Gaugino pole masses:"<<endl<<endl;
379  get_polemass("~g");
380  get_polemass("~chi+",1,2);
381  get_polemass("~chi0",1,4);
382  cout<<endl<<"Squark pole masses:"<<endl<<endl;
383  get_polemass("~d",1,6);
384  get_polemass("~u",1,6);
385  cout<<endl<<"Slepton pole masses:"<<endl<<endl;
386  get_polemass("~e-",1,6);
387  get_polemass("~nu",1,3);
388 
389  cout << endl << "Mixing matrices:" << endl << endl;
390 
391  // Note, currently we are not using a matrix object or any such thing, so you have to
392  // extract the elements of each matrix one at a time. It would probably be handy to
393  // add such a return type though.
394 
395  #define GET_MIX_MATRIX_EL(r, PRODUCT) \
396  { \
397  str label = BOOST_PP_SEQ_ELEM(0, PRODUCT); \
398  str block = BOOST_PP_SEQ_ELEM(1, PRODUCT); \
399  int i = BOOST_PP_SEQ_ELEM(2, PRODUCT); \
400  int j = BOOST_PP_SEQ_ELEM(3, PRODUCT); \
401  try{ \
402  std::ostringstream echo1; \
403  std::ostringstream echo2; \
404  echo1 << " spec->get_Pole_Mixing("<<label<<","<<i<<","<<j<<") = "; \
405  double value1 = spec->get_Pole_Mixing(label,i,j); \
406  echo2 << " SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
407  double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
408  cout << echo1.str() << value1 <<endl; \
409  cout << echo2.str() << value2 <<endl; \
410  cout << endl; \
411  } \
412  catch (const std::exception& e) \
413  { add_error(report,e,label+": "+block); } \
414  }
415 
416  #define GET_MIX_MATRIX(NAME,BLOCK,__IND1,__IND2) BOOST_PP_SEQ_FOR_EACH_PRODUCT(GET_MIX_MATRIX_EL, ((NAME))((BLOCK))(BOOST_PP_TUPLE_TO_SEQ(__IND1))(BOOST_PP_TUPLE_TO_SEQ(__IND2)))
417 
418  // The names here could perhaps be improved. They are not so immediately obvious to me.
419 
420  GET_MIX_MATRIX("~chi-","UMIX",(1,2),(1,2)) cout<<endl;
421  GET_MIX_MATRIX("~chi+","VMIX",(1,2),(1,2)) cout<<endl;
422  GET_MIX_MATRIX("A0","PSEUDOSCALARMIX",(1,2),(1,2)) cout<<endl;
423  GET_MIX_MATRIX("~d","DSQMIX",(1,2,3,4,5,6),(1,2,3,4,5,6)) cout<<endl;
424  GET_MIX_MATRIX("~e-","SELMIX",(1,2,3,4,5,6),(1,2,3,4,5,6)) cout<<endl;
425  GET_MIX_MATRIX("h0","SCALARMIX",(1,2),(1,2)) cout<<endl;
426  GET_MIX_MATRIX("~chi0","NMIX",(1,2,3,4),(1,2,3,4)) cout<<endl;
427  GET_MIX_MATRIX("H+","CHARGEMIX",(1,2),(1,2)) cout<<endl;
428  GET_MIX_MATRIX("~u","USQMIX",(1,2,3,4,5,6),(1,2,3,4,5,6)) cout<<endl;
429  GET_MIX_MATRIX("~nu","SNUMIX",(1,2,3),(1,2,3)) cout<<endl;
430 
431  cout<<endl;
432  cout << "Next up: running parameters" << endl;
433  cout << "These are all given in the DRbar scheme, at least when running FlexibleSUSY or SoftSUSY. ";
434  cout << "There may be some switching or converting once other spectrum generator are added." << endl;
435  cout<<endl;
436  cout << "Spectrum object running parameters are currently defined at scale Q="
437  << spec->GetScale() << " [GeV]" << endl << endl;
438  cout<<endl;
439  cout << "-- Dimensionless parameters --" <<endl;
440  cout << endl << "Gauge couplings:" << endl << endl;
441  ECHO( spec->get_dimensionless_parameter("g1") ) // U_Y(1) gauge coupling in SU(5) normalisation
442  ECHO( slhaea.at("GAUGE").at(1).at(1) ) // This is in the Standard Model normalisation as per SLHA conventions
443  cout << "Note: " << spec->get_dimensionless_parameter("g1") << " * sqrt(3/5) = "
444  << spec->get_dimensionless_parameter("g1")*sqrt(3./5.) << endl;
445  cout<<endl;
446  ECHO( spec->get_dimensionless_parameter("g2") ) // SU(2) gauge coupling
447  ECHO( slhaea.at("GAUGE").at(2).at(1) )
448  cout<<endl;
449  ECHO( spec->get_dimensionless_parameter("g3") ) // SU(3) gauge coupling
450  ECHO( slhaea.at("GAUGE").at(3).at(1) )
451  cout<<endl;
452 
453  cout << endl << "Yukawa matrices:" << endl << endl;
454 
455  // Note, currently we are not using a matrix object or any such thing, so you have to
456  // extract the elements of the matrix one at a time. It would probably be handy to
457  // add such a return type though.
458 
459  #define GET_MATRIX_EL(r, PRODUCT) \
460  { \
461  str label = BOOST_PP_SEQ_ELEM(0, PRODUCT); \
462  str block = BOOST_PP_SEQ_ELEM(1, PRODUCT); \
463  int i = BOOST_PP_SEQ_ELEM(2, PRODUCT); \
464  int j = BOOST_PP_SEQ_ELEM(3, PRODUCT); \
465  try{ \
466  std::ostringstream echo1; \
467  std::ostringstream echo2; \
468  echo1 << " spec->get_dimensionless_parameter("<<label<<","<<i<<","<<j<<") = "; \
469  double value1 = spec->get_dimensionless_parameter(label,i,j); \
470  echo2 << " SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
471  double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
472  cout << echo1.str() << value1 <<endl; \
473  cout << echo2.str() << value2 <<endl; \
474  cout << endl; \
475  } catch (const std::exception& e) \
476  { add_error(report,e,label+": "+block); } \
477  }
478 
479  #define GET_MATRIX(NAME,BLOCK,__IND1,__IND2) BOOST_PP_SEQ_FOR_EACH_PRODUCT(GET_MATRIX_EL, ((NAME))((BLOCK))(BOOST_PP_TUPLE_TO_SEQ(__IND1))(BOOST_PP_TUPLE_TO_SEQ(__IND2)))
480 
481  GET_MATRIX("Yu","YU",(1,2,3),(1,2,3)) cout << endl;
482  GET_MATRIX("Yd","YD",(1,2,3),(1,2,3)) cout << endl;
483  GET_MATRIX("Ye","YE",(1,2,3),(1,2,3)) cout << endl;
484 
485  // Mass dimension 1 parameters
486 
487  cout<<endl;
488  cout<<"MSSM mass dimension 1 running parameters"<<endl;
489  cout<<endl;
490  ECHO( spec->get_mass_parameter("M1") ) // Gaugino mass parameter "MassB"
491  ECHO( slhaea.at("MSOFT").at(1).at(1) )
492  cout<<endl;
493  ECHO( spec->get_mass_parameter("M2") ) // Gaugino mass parameter "MassWB"
494  ECHO( slhaea.at("MSOFT").at(2).at(1) )
495  cout<<endl;
496  ECHO( spec->get_mass_parameter("M3") ) // Gaugino mass parameter "MassG"
497  ECHO( slhaea.at("MSOFT").at(3).at(1) )
498  cout<<endl;
499  ECHO( spec->get_mass_parameter("Mu") ) // Superpotential mu parameter
500  ECHO( slhaea.at("HMIX").at(1).at(1) )
501  cout<<endl;
502  ECHO( spec->get_mass_parameter("vd") ) // Down-type Higgs vev
503  ECHO( slhaea.at("HMIX").at(102).at(1) )
504  cout<<endl;
505  ECHO( spec->get_mass_parameter("vu") ) // Up-type Higgs vev
506  ECHO( slhaea.at("HMIX").at(103).at(1) )
507  cout<<endl;
508 
509  // Matrices
510 
511  #define GET_M1_MATRIX_EL(r, PRODUCT) \
512  { \
513  str label = BOOST_PP_SEQ_ELEM(0, PRODUCT); \
514  str block = BOOST_PP_SEQ_ELEM(1, PRODUCT); \
515  int i = BOOST_PP_SEQ_ELEM(2, PRODUCT); \
516  int j = BOOST_PP_SEQ_ELEM(3, PRODUCT); \
517  try{ \
518  std::ostringstream echo1; \
519  std::ostringstream echo2; \
520  echo1 << " spec->get_mass_parameter("<<label<<","<<i<<","<<j<<") = "; \
521  double value1 = spec->get_mass_parameter(label,i,j); \
522  echo2 << " SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
523  double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
524  cout << echo1.str() << value1 <<endl; \
525  cout << echo2.str() << value2 <<endl; \
526  cout << endl; \
527  } catch (const std::exception& e) \
528  { add_error(report,e,label+": "+block); } \
529  }
530 
531  #define GET_M1_MATRIX(NAME,BLOCK,__IND1,__IND2) BOOST_PP_SEQ_FOR_EACH_PRODUCT(GET_M1_MATRIX_EL, ((NAME))((BLOCK))(BOOST_PP_TUPLE_TO_SEQ(__IND1))(BOOST_PP_TUPLE_TO_SEQ(__IND2)))
532 
533  cout << endl << "Triliner coupling matrices? SLHA says these blocks should be called AU,AD,AE, not TU,TD,TE though, so I'm not sure. Need to check with Peter." << endl << endl;
534 
535  // Seem to be the trilinears, and TYu and au etc. seem to be equal. Ask Peter...
536 
537  GET_M1_MATRIX("TYu","TU",(1,2,3),(1,2,3)) cout << endl;
538  GET_M1_MATRIX("TYd","TD",(1,2,3),(1,2,3)) cout << endl;
539  GET_M1_MATRIX("TYe","TE",(1,2,3),(1,2,3)) cout << endl;
540  cout << endl;
541  GET_M1_MATRIX("au","TU",(1,2,3),(1,2,3)) cout << endl;
542  GET_M1_MATRIX("ad","TD",(1,2,3),(1,2,3)) cout << endl;
543  GET_M1_MATRIX("ae","TE",(1,2,3),(1,2,3)) cout << endl;
544 
545  // Mass dimension 2 parameters
546 
547  cout<<endl;
548  cout<<"MSSM mass dimension 2 running parameters"<<endl;
549  cout<<endl;
550  ECHO( spec->get_mass2_parameter("mHd2") ) // Down-type Higgs soft mass
551  ECHO( slhaea.at("MSOFT").at(21).at(1) )
552  cout<<endl;
553  ECHO( spec->get_mass2_parameter("mHu2") ) // Up-type Higgs soft mass
554  ECHO( slhaea.at("MSOFT").at(22).at(1) )
555  cout<<endl;
556  ECHO( spec->get_mass2_parameter("BMu") ) // Higgs bilinear soft parameter
557  ECHO( slhaea.at("HMIX").at(101).at(1) )
558  cout<<endl;
559 
560  // Matrices
561 
562  #define GET_M2_MATRIX_EL(r, PRODUCT) \
563  { \
564  str label = BOOST_PP_SEQ_ELEM(0, PRODUCT); \
565  str block = BOOST_PP_SEQ_ELEM(1, PRODUCT); \
566  int i = BOOST_PP_SEQ_ELEM(2, PRODUCT); \
567  int j = BOOST_PP_SEQ_ELEM(3, PRODUCT); \
568  std::ostringstream echo1; \
569  try { \
570  std::ostringstream echo2; \
571  echo1 << " spec->get_mass2_parameter("<<label<<","<<i<<","<<j<<") = "; \
572  double value1 = spec->get_mass2_parameter(label,i,j); \
573  echo2 << " SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
574  double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
575  cout << echo1.str() << value1 <<endl; \
576  cout << echo2.str() << value2 <<endl; \
577  cout << endl; \
578  } catch (const std::exception& e) \
579  { add_error(report,e,label+": "+block); } \
580  }
581 
582  #define GET_M2_MATRIX(NAME,BLOCK,__IND1,__IND2) BOOST_PP_SEQ_FOR_EACH_PRODUCT(GET_M2_MATRIX_EL, ((NAME))((BLOCK))(BOOST_PP_TUPLE_TO_SEQ(__IND1))(BOOST_PP_TUPLE_TO_SEQ(__IND2)))
583 
584  cout << endl << "Mass matrices:" << endl << endl;
585 
586  GET_M2_MATRIX("mq2","MSQ2",(1,2,3),(1,2,3)) cout << endl;
587  GET_M2_MATRIX("mu2","MSU2",(1,2,3),(1,2,3)) cout << endl;
588  GET_M2_MATRIX("md2","MSD2",(1,2,3),(1,2,3)) cout << endl;
589  GET_M2_MATRIX("me2","MSE2",(1,2,3),(1,2,3)) cout << endl;
590  GET_M2_MATRIX("ml2","MSL2",(1,2,3),(1,2,3)) cout << endl;
591 
592  cout << endl;
593 
595 
596  cout << "Testing set_override functions" << endl;
597 
598  cout << "Original M1:" << spec->get_mass_parameter("M1") << endl;
599  spec->set_override_mass_parameter(-666,"M1");
600  cout << "Override M1:" << spec->get_mass_parameter("M1") << endl;
601 
602  cout << "Original ~e-(1):" << spec->get_Pole_Mass("~e-",1) << endl;
603  spec->set_override_Pole_Mass(-667,"~e-",1);
604  cout << "Override ~e-(1):" << spec->get_Pole_Mass("~e-",1) << endl;
605 
606  cout << "Original ml2(1,1):" << spec->get_mass2_parameter("ml2",1,1) << endl;
607  spec->set_override_mass2_parameter(-668,"ml2",1,1);
608  cout << "Override ml2(1,1):" << spec->get_mass2_parameter("ml2",1,1) << endl;
609 
611  cout << "has 'new_entry'? " << spec->has_mass_parameter("new_entry") << endl;
612  cout << "..." << endl;
614  //spec->set_override_mass2_parameter(-1234,"new_entry"); // incorrect: safety still on
615  spec->set_override_mass_parameter(-1234,"new_entry",false); // correct: safety check turned off
616  cout << "has 'new_entry'? " << spec->has_mass_parameter("new_entry") << endl;
617  cout << "new_entry = " << spec->get_mass_parameter("new_entry") << endl;
618  cout << endl;
619 
620 
621  // Things left to add to demo:
622  // - Tree-level masses
623  // - Standard Model mixings (need to be added to MSSMSpec as well)
624  // - some extra stuff to do regarding organising Standard Model input
625  // parameters (and computing/transferring the pole masses)
626 
627  cout << "Test report:" << std::endl << report.str();
628 
629  SpecBit_warning().raise(LOCAL_INFO,"\n *** Stopped on purpose to examine spectrum contents ***");
630  result = 0;
631  }
632  //else if(spinfo.find(k3) != spinfo.end())
633  //{
634  // std::cout << "Bad spectrum: " << spinfo.at(3) << std::endl;
635  // result = 1;
636  //}
637  else if(spinfo.find(k4) != spinfo.end())
638  {
639  std::cout << "Bad spectrum: " << spinfo.at(4) << std::endl;
640  result = 1;
641  }
642  else
643  {
644  SpecBit_error().raise(LOCAL_INFO,"Bug in exampleRead 'if' logic");
645  }
646  }
647 
648  }
649 }
650 
651 
652 
653 
654 
655 
656 
657 
658 
659 
660 
661 
662 
Define overloadings of the stream operator for various containers.
This class is used to deliver both information defined in the Standard Model (or potentially just QED...
Rollcall header for module SpecBit.
#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:432
virtual double GetScale() const
Returns the renormalisation scale of parameters.
virtual void writeSLHAfile(int, const str &) const
Dump out spectrum information to an SLHA file (if possible)
Definition: subspectrum.cpp:40
#define PDB
void exampleRead(bool &result)
#define GET_MATRIX(NAME, BLOCK, __IND1, __IND2)
void add_error(std::ostringstream &out, const std::exception &e, const std::string &msg)
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
#define GET_MIX_MATRIX(NAME, BLOCK, __IND1, __IND2)
SubSpectrum & get_LE()
Standard getters Return references to internal data members.
Definition: spectrum.cpp:224
#define GET_M2_MATRIX(NAME, BLOCK, __IND1, __IND2)
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
#define GET_M1_MATRIX(NAME, BLOCK, __IND1, __IND2)
#define ECHO(COMMAND)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
SubSpectrum & get_HE()
Definition: spectrum.cpp:225
"Standard Model" (low-energy) plus high-energy model container class
Definition: spectrum.hpp:110