file examples/SpecBit_examples.cpp

[No description available] More…

Namespaces

Name
Gambit
TODO: see if we can use this one:
Gambit::SpecBit

Defines

Name
ECHO(COMMAND)
PDB
GET_MIX_MATRIX_EL(r, PRODUCT)
GET_MIX_MATRIX(NAME, BLOCK, __IND1, __IND2)
GET_MATRIX_EL(r, PRODUCT)
GET_MATRIX(NAME, BLOCK, __IND1, __IND2)
GET_M1_MATRIX_EL(r, PRODUCT)
GET_M1_MATRIX(NAME, BLOCK, __IND1, __IND2)
GET_M2_MATRIX_EL(r, PRODUCT)
GET_M2_MATRIX(NAME, BLOCK, __IND1, __IND2)

Detailed Description

Author: Ben Farmer (ben.farmer@gmail.com)

Date: 2014 Dec

Example retrieval of data from Spectrum object.


Authors (add name and date if you modify):


Macros Documentation

define ECHO

#define ECHO(
    COMMAND
)
         {                                                    \
             try {                                            \
                cout << "  " << STRINGIFY(COMMAND) << " = " << COMMAND << endl;\
             }                                                \
             catch (const std::exception& e)                  \
             { add_error(report,e,STRINGIFY(COMMAND)); }      \
         }

define PDB

#define PDB Models::ParticleDB()

define GET_MIX_MATRIX_EL

#define GET_MIX_MATRIX_EL(
    r,
    PRODUCT
)
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
           try{                                                                  \
            std::ostringstream echo1;                                            \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_Pole_Mixing("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_Pole_Mixing(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           }                                                                     \
           catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
         }

define GET_MIX_MATRIX

#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)))

define GET_MATRIX_EL

#define GET_MATRIX_EL(
    r,
    PRODUCT
)
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
           try{                                                                  \
            std::ostringstream echo1;                                            \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_dimensionless_parameter("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_dimensionless_parameter(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           } catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
         }

define GET_MATRIX

#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)))

define GET_M1_MATRIX_EL

#define GET_M1_MATRIX_EL(
    r,
    PRODUCT
)
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
           try{                                                                  \
            std::ostringstream echo1;                                            \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_mass_parameter("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_mass_parameter(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           } catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
         }

define GET_M1_MATRIX

#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)))

define GET_M2_MATRIX_EL

#define GET_M2_MATRIX_EL(
    r,
    PRODUCT
)
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
            std::ostringstream echo1;                                            \
           try {                                                                 \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_mass2_parameter("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_mass2_parameter(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           } catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
        }

define GET_M2_MATRIX

#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)))

Source code

//   GAMBIT: Global and Modular BSM Inference Tool
//   *********************************************
///  \file
///
///  Example retrieval of data from Spectrum
///  object.
///
///  *********************************************
///
///  Authors (add name and date if you modify):
///
///  \author Ben Farmer
///          (ben.farmer@gmail.com)
///    \date 2014 Dec
///
///  *********************************************

#include <string>
#include <sstream>
#include <cmath>

#include "gambit/Utils/gambit_module_headers.hpp"
#include "gambit/Elements/spectrum.hpp"
#include "gambit/SpecBit/SpecBit_rollcall.hpp"

#include "SLHAea/slhaea.h"
#include <boost/preprocessor/tuple/to_seq.hpp>
#include <boost/preprocessor/seq/for_each.hpp>
#include <boost/preprocessor/seq/elem.hpp>
#include <boost/preprocessor/seq/for_each_product.hpp>

// For the demo output only
#include "gambit/Utils/stream_overloads.hpp"

namespace Gambit
{

  namespace SpecBit
  {

    using namespace LogTags;

    // Helper function to add error information to report
    void add_error(std::ostringstream& out, const std::exception& e, const std::string& msg)
    {
       out << "------------------------------" << std::endl;
       out << "TEST FAIL: " << msg << std::endl;
       out << "Exception thrown was: "<<e.what()<<std::endl;
       return;
    }

    // Example showing reading of SLHA information from Spectrum object, using
    // both SLHAea objects and direct Spectrum getter functions.
    void exampleRead (bool &result)
    {

      // Retreive pointer to Spectrum object, delivered by dependency resolver
      // Module function asks for Spectrum* with capability MSSM_spectrum.
      // i.e. has DEPENDENCY(MSSM_spectrum, Spectrum*)
      namespace myPipe = Pipes::exampleRead;
      /*TAG*/ Spectrum fullspectrum = *myPipe::Dep::MSSM_spectrum;
      const SubSpectrum* spec = fullspectrum->get_HE(); // MSSMSpec SubSpectrum object
      const SubSpectrum* SM   = fullspectrum->get_LE(); // QedQcdWrapper SubSpectrum object

      std::ostringstream report; // Information about any problems encountered

      // Extract SLHAea object
      // This copies the data out. Could possible change it to pass out a
      // reference instead, or have another function to do that.
      SLHAea::Coll slhaea = fullspectrum->getSLHAea(2);
      // for testing, write this to file
      std::ofstream out1;
      out1.open("SpecBit/exampleRead_test.slha");
      out1 << slhaea;
      out1.close();

      // SLHAea::Coll slhaea = spec->getSLHAea(2); // The above is just a wrapper for this.

      // If this is a valid model point, return true and dump information, else false

      // SLHAea objects behave mostly like maps, but with special kinds of keys. For
      // "at" and "operator[]", it does automatic conversion, but for "find" it does
      // not, so we have to manually do it.
      SLHAea::Block spinfo = slhaea.at("SPINFO");
      //std::vector<std::string> k3(1, "3");
      std::vector<std::string> k4(1, "4");

      // See if error code entries exist
      //if(spinfo.find(k3) == spinfo.end() and spinfo.find(k4) == spinfo.end())
      if(spinfo.find(k4) == spinfo.end())
      {
         std::cout << "Good spectrum found! Inspecting contents..." << std::endl;
         std::cout << std::endl << slhaea << std::endl;

         // Write to file so we can check it
         std::string filename(myPipe::runOptions->getValue<std::string>("output_slha_filename"));
         spec->writeSLHAfile(2,filename);

         // ---------------------------------------------------------
         // BEGIN DEMO OF SPECTRUM OBJECT AND PARTICLE DATABASE
         // ---------------------------------------------------------

         //#define ECHO(COMMAND) cout << "  " << STRINGIFY(COMMAND) << " = " << COMMAND << endl;
         //Replacing with a version that deals with SLHAea access errors
         #define ECHO(COMMAND)                                \
         {                                                    \
             try {                                            \
                cout << "  " << STRINGIFY(COMMAND) << " = " << COMMAND << endl;\
             }                                                \
             catch (const std::exception& e)                  \
             { add_error(report,e,STRINGIFY(COMMAND)); }      \
         }

         /* ----------Test particle database access ---------------- */
         #define PDB Models::ParticleDB()

         // First check out what is actually in the database
         PDB.check_contents();

         // Demo a couple of particle name retrievals
         cout<<endl;
         cout<<"Demo retrieval of lightest neutralino info from particle database"<<endl;
         cout<<"-----------------------------------------------------------------"<<endl;
         ECHO( PDB.pdg_pair("~chi0_1")                   ) // Input long name, retrieve PDG code + context integer
         ECHO( PDB.pdg_pair("~chi0",1)                   ) // Input short name + index, retrieve PDG code + context integer
         ECHO( PDB.long_name("~chi0",1)                  ) // Input short name + index, retrieve long name
         ECHO( PDB.long_name(std::make_pair(1000022,0))  ) // Input PDG code + context integer, retrieve long name
         ECHO( PDB.long_name(1000022,0)                  ) // Input PDG code + context integer, retrieve long name
         ECHO( PDB.short_name_pair("~chi0_1")            ) // Input long name, retrieve short name + index
         ECHO( PDB.short_name_pair(std::make_pair(1000022,0)) ) // Input PDG code plus context integer, retrieve short name + index
         ECHO( PDB.short_name_pair(1000022,0)            ) // Input PDG code plus context integer, retrieve short name + index

         cout<<endl;
         cout<<"Demo retrieval of antiparticle names/codes using particle names/codes"<<endl;
         cout<<"-----------------------------------------------------------------"<<endl;
         // Check existence in various ways
         ECHO( PDB.has_antiparticle("~chi0_1")                 )
         ECHO( PDB.has_antiparticle("~chi0",1)                 )
         ECHO( PDB.has_antiparticle(std::make_pair("~chi0",1)) )
         ECHO( PDB.has_antiparticle(1000022,0)                 )
         ECHO( PDB.has_antiparticle(std::make_pair(1000022,0)) )
         ECHO( PDB.get_antiparticle("~chi0_1")                 ) // Input long name, retrieve antiparticle long name
         ECHO( PDB.get_antiparticle("~chi0",1)                 ) // Input short name + index, retrieve antiparticle short name + index
         ECHO( PDB.get_antiparticle(std::make_pair("~chi0",1)) ) // Input short name + index, retrieve antiparticle short name + index
         ECHO( PDB.get_antiparticle(1000022,0)                 ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
         ECHO( PDB.get_antiparticle(std::make_pair(1000022,0)) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
         ECHO( PDB.has_antiparticle("~chi+_1")                 )
         ECHO( PDB.has_antiparticle("~chi+",1)                 )
         ECHO( PDB.has_antiparticle(std::make_pair("~chi+",1)) )
         ECHO( PDB.has_antiparticle(1000024,0)                 )
         ECHO( PDB.has_antiparticle(std::make_pair(1000024,0)) )
         ECHO( PDB.get_antiparticle("~chi+_1")                 ) // Input long name, retrieve antiparticle long name
         ECHO( PDB.get_antiparticle("~chi+",1)                 ) // Input short name + index, retrieve antiparticle short name + index
         ECHO( PDB.get_antiparticle(std::make_pair("~chi+",1)) ) // Input short name + index, retrieve antiparticle short name + index
         ECHO( PDB.get_antiparticle(1000024,0)                 ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
         ECHO( PDB.get_antiparticle(std::make_pair(1000024,0)) ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
         ECHO( PDB.has_antiparticle("u_1")                 )
         ECHO( PDB.has_antiparticle("u",1)                 )
         ECHO( PDB.has_antiparticle(std::make_pair("u",1)) )
         ECHO( PDB.has_antiparticle(2,0)                   )
         ECHO( PDB.has_antiparticle(std::make_pair(2,0))   )
         ECHO( PDB.get_antiparticle("u_1")                 ) // Input long name, retrieve antiparticle long name
         ECHO( PDB.get_antiparticle("u",1)                 ) // Input short name + index, retrieve antiparticle short name + index
         ECHO( PDB.get_antiparticle(std::make_pair("u",1)) ) // Input short name + index, retrieve antiparticle short name + index
         ECHO( PDB.get_antiparticle(2,0)                   ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer
         ECHO( PDB.get_antiparticle(std::make_pair(2,0))   ) // Input PDG code + context integet, retrieve antiparticle PDG code + context integer


         cout<<endl;
         cout<<"Demo retrieval when no short name exists"<<endl;
         cout<<"-----------------------------------------------------------------"<<endl;
         ECHO( PDB.pdg_pair("H+")                        )
         //ECHO( PDB.pdg_pair("H+",1)                      ) // Error!
         //ECHO( PDB.long_name("H+",1)                     ) // Error!
         ECHO( PDB.long_name(std::make_pair(37,0))       )
         ECHO( PDB.long_name(37,0)                       )
         //ECHO( PDB.short_name_pair("H+")                 ) // Error!
         //ECHO( PDB.short_name_pair(std::make_pair(37,0)) ) // Error!
         //ECHO( PDB.short_name_pair(37,0)                 ) // Error!
         //ECHO( PDB.short_name_pair(37)                   ) // Error!
         cout<<endl;
         /* ----------------- Pole masses --------------------------- */

         cout<<"Begin demo retrievals from Spectrum and SubSpectrum objects"<<endl;
         cout<<"-----------------------------------------------------------------"<<endl;
         cout<<endl;
         cout<<"First, general methods for accessing different sorts of information."<<endl;
         cout<<endl;
         // At the moment it is only pole masses which have getters overloaded to use
         // the particle database information. It is only the MASS block in the
         // spectrum generator output SLHA files which use PDG numbers anyway, so I
         // think this makes sense.
         cout<<"Lightest neutral Higgs boson pole mass:"<<endl;
         ECHO(  fullspectrum->get_Pole_Mass( PDB.short_name_pair(25,0) )   )
         ECHO(  fullspectrum->get_Pole_Mass( PDB.long_name(25,0) )         )
         ECHO(  fullspectrum->get_Pole_Mass(25,0)                          )
         ECHO(  fullspectrum->get_Pole_Mass( PDB.pdg_pair("h0",1) )        )
         ECHO(  fullspectrum->get_Pole_Mass("h0",1)                        )
         ECHO(  fullspectrum->get_Pole_Mass("h0_1")                        )

         ECHO(  spec->get_Pole_Mass( PDB.short_name_pair(25,0) )   )
         ECHO(  spec->get_Pole_Mass( PDB.long_name(25,0) )         )
         ECHO(  spec->get_Pole_Mass(25,0)                          )
         ECHO(  spec->get_Pole_Mass( PDB.pdg_pair("h0",1) )        )
         ECHO(  spec->get_Pole_Mass("h0",1)                        )
         ECHO(  spec->get_Pole_Mass("h0_1")                        )

         cout<<endl;
         cout<<"Retrieval of Spectrum object contents, with"<<endl;
         cout<<"correspondence to SLHAea object entries"<<endl;
         cout<<"-----------------------------------------------------------------"<<endl;

         // MZ was a bad first example; it is empty unless you switch on the SM pole mass
         // calculator for flexiblesusy. We do not yet pass any input value of MZ though
         // Similar issues with other gauge boson masses. So don't use these yet or
         // you'll get zero for all these masses.
         cout<<endl;
         cout<<"Gauge boson pole masses:"<<endl;
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("Z0") )
         ECHO(  SM->get_Pole_Mass("Z0")           )
         ECHO(  slhaea.at("SMINPUTS").at(4).at(1)      )
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("gamma")    )
         ECHO(  SM->get_Pole_Mass("gamma")    )
         cout<<"  ***Not in slha***"<<endl;
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("W+")       )
         ECHO(  SM->get_Pole_Mass("W+")         )
         ECHO(  slhaea.at("MASS").at(24).at(1)    )
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("g")       )
         ECHO(  SM->get_Pole_Mass("g")          )
         cout<<"  ***Not in slha***"<<endl;
         cout<<endl;
         cout<<"Quark pole masses (actually the slha entries are MSbar except the top mass):"<<endl;
         cout<<endl;
         // I'm a little unclear on what the pole masses for the lighter quarks mean, since I thought
         // that non-perturbative effects made definining them difficult... well anyway will have
         // to ask Peter what is being computed here.

         //ECHO(  spec->get_Pole_Mass("u",1)      )  // i.e. up (mass eigenstate)
         cout<<"  ***u Pole mass not well defined***"<<endl;
         ECHO(  slhaea.at("SMINPUTS").at(22).at(1)   )  // mu(2 GeV)^MS-bar, not pole mass
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("u",2)      )  // i.e. charm
         cout<<"  ***c Pole mass not well defined***"<<endl;
         ECHO(  slhaea.at("SMINPUTS").at(24).at(1)   )  // mc(mc)^MS-bar, not pole mass
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("u",3)      )  // i.e. top
         ECHO(  fullspectrum->get_Pole_Mass("t")       )
         ECHO(  SM->get_Pole_Mass("u",3)      )  // i.e. top
         ECHO(  slhaea.at("SMINPUTS").at(6).at(1)    )
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("d",1)      )  // i.e. down
         cout<<"  ***d Pole mass not well defined***"<<endl;
         ECHO(  slhaea.at("SMINPUTS").at(21).at(1)   )  // md(2 GeV)^MS-bar, not pole mass
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("d",2)      )  // i.e. strange
         cout<<"  ***s Pole mass not well defined***"<<endl;
         ECHO(  slhaea.at("SMINPUTS").at(23).at(1)   )  // ms(2 GeV)^MS-bar, not pole mass
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("d",3)      )  // i.e. bottom
         ECHO(  fullspectrum->get_Pole_Mass("b")       )
         ECHO(  SM->get_Pole_Mass("d",3)      )  // i.e. bottom
         ECHO(  slhaea.at("SMINPUTS").at(5).at(1)    )  //  mb(mb)^MS-bar, not pole mass.
         cout<<endl;
         cout<<"Charged fermions pole masses:"<<endl;
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("e-",1)     )  // i.e. electron
         ECHO(  fullspectrum->get_Pole_Mass("e-")       )
         ECHO(  SM->get_Pole_Mass("e-",1)     )  // i.e. electron
         ECHO(  slhaea.at("SMINPUTS").at(11).at(1)   )
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("e-",2)     )  // i.e. muon
         ECHO(  fullspectrum->get_Pole_Mass("mu-")       )
         ECHO(  SM->get_Pole_Mass("e-",2)     )  // i.e. muon
         ECHO(  slhaea.at("SMINPUTS").at(13).at(1)   )
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("e-",3)     )  // i.e. tau
         ECHO(  fullspectrum->get_Pole_Mass("tau-")       )
         ECHO(  SM->get_Pole_Mass("e-",3)     )  // i.e. tau
         ECHO(  slhaea.at("SMINPUTS").at(7).at(1)    )
         cout<<endl;
         cout<<"Neutrinos pole masses:"<<endl;
         cout<<endl;
         // These will produce errors because currently no neutrino mass getters are hooked up
         //ECHO(  spec->get_Pole_Mass("nu",1)     )  // Just mass ordered (if there is mixing)
         ECHO(  fullspectrum->get_Pole_Mass("nu",1)       )
         ECHO(  SM->get_Pole_Mass("nu",1)     )  // Just mass ordered (if there is mixing)
         ECHO(  slhaea.at("SMINPUTS").at(12).at(1)   )
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("nu",2)     )
         ECHO(  fullspectrum->get_Pole_Mass("nu",2)       )
         ECHO(  SM->get_Pole_Mass("nu",2)     )
         ECHO(  slhaea.at("SMINPUTS").at(14).at(1)   )
         cout<<endl;
         //ECHO(  spec->get_Pole_Mass("nu",3)     )
         ECHO(  fullspectrum->get_Pole_Mass("nu",3)       )
         ECHO(  SM->get_Pole_Mass("nu",3)     )
         ECHO(  slhaea.at("SMINPUTS").at(8).at(1)    )
         cout<<endl;
         // Now for SUSY particles
         cout<<endl;
         cout<<"MSSM Higgs sector pole masses:"<<endl;
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("h0",1)       )
         ECHO(  spec->get_Pole_Mass("h0",1)     )  // Lightest neutral Higgs boson
         ECHO(  slhaea.at("MASS").at(25).at(1)       )
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("h0",2)       )
         ECHO(  spec->get_Pole_Mass("h0",2)     )  // Heavy neutral Higgs boson
         ECHO(  slhaea.at("MASS").at(35).at(1)       )
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("H+")       )
         ECHO(  spec->get_Pole_Mass("H+")       )  // Charged Higgs
         ECHO(  slhaea.at("MASS").at(37).at(1)       )
         cout<<endl;
         ECHO(  fullspectrum->get_Pole_Mass("A0")       )
         ECHO(  spec->get_Pole_Mass("A0")       )  // Pseudoscalar neutral Higgs
         ECHO(  slhaea.at("MASS").at(36).at(1)       )
         cout<<endl;

         // I'm going to use these nested functors to save lots of typing for the rest. It is just the
         // same as the examples above, except that the PDG codes are retrieved from the particle database.
         // The PDG code - string name correspondences are defined in 'Models/src/particle_database.cpp'

         struct get_polemass_functor
         {
           // Single mass
           void operator()(const std::string& longname)
           {
             std::ostringstream echo1, echo2, echo3;
             echo1 <<     "  fullspectrum->get_Pole_Mass("<<longname<<") = ";
             double value1 = fullspectrum->get_Pole_Mass(longname);
             echo2 <<     "  spec->get_Pole_Mass("<<longname<<") = ";
             double value2 = spec->get_Pole_Mass(longname);
             echo3 <<  "  slhaea.at(\"MASS\").at("<<PDB.pdg_pair(longname).first<<").at(1) = ";
             str value3 = slhaea.at("MASS").at( PDB.pdg_pair(longname).first ).at(1);
             cout << echo1.str() << value1 << endl;
             cout << echo2.str() << value2 << endl;
             cout << echo3.str() << value3 << endl;
             cout<<endl;
           }
           // Range of indexes masses
           void operator()(const std::string& longname, int from, int to)
           {
             for(int i=from; i<=to; ++i)
             {
               std::ostringstream echo1;
               std::ostringstream echo2;
               echo1 <<     "  spec->get_Pole_Mass("<<longname<<","<<i<<") = ";
               double value1 = spec->get_Pole_Mass(longname,i);
               echo2 <<  "  slhaea.at(\"MASS\").at("<<PDB.pdg_pair(longname,i).first<<").at(1) = ";
               str value2 = slhaea.at("MASS").at( PDB.pdg_pair(longname,i).first ).at(1);
               cout << echo1.str() << value1 << endl;
               cout << echo2.str() << value2 << endl;
               cout<<endl;
             }
           }

           get_polemass_functor(std::ostringstream& report, /*TAG*/ Spectrum fullin, const SubSpectrum* specin, SLHAea::Coll& slhaeain)
             : report(report)
             , fullspectrum(fullin)
             , spec(specin)
             , slhaea(slhaeain)
           {}

           private:
             std::ostringstream& report;
             /*TAG*/ Spectrum fullspectrum;
             const SubSpectrum* spec;
             SLHAea::Coll slhaea;
         };

         get_polemass_functor get_polemass(report,fullspectrum,spec,slhaea);

         cout<<endl<<"Gaugino pole masses:"<<endl<<endl;
         get_polemass("~g");
         get_polemass("~chi+",1,2);
         get_polemass("~chi0",1,4);
         cout<<endl<<"Squark pole masses:"<<endl<<endl;
         get_polemass("~d",1,6);
         get_polemass("~u",1,6);
         cout<<endl<<"Slepton pole masses:"<<endl<<endl;
         get_polemass("~e-",1,6);
         get_polemass("~nu",1,3);

         cout << endl << "Mixing matrices:" << endl << endl;

         // Note, currently we are not using a matrix object or any such thing, so you have to
         // extract the elements of each matrix one at a time. It would probably be handy to
         // add such a return type though.

         #define GET_MIX_MATRIX_EL(r, PRODUCT)                                   \
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
           try{                                                                  \
            std::ostringstream echo1;                                            \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_Pole_Mixing("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_Pole_Mixing(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           }                                                                     \
           catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
         }

         #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)))

         // The names here could perhaps be improved. They are not so immediately obvious to me.

         GET_MIX_MATRIX("~chi-","UMIX",(1,2),(1,2)) cout<<endl;
         GET_MIX_MATRIX("~chi+","VMIX",(1,2),(1,2)) cout<<endl;
         GET_MIX_MATRIX("A0","PSEUDOSCALARMIX",(1,2),(1,2)) cout<<endl;
         GET_MIX_MATRIX("~d","DSQMIX",(1,2,3,4,5,6),(1,2,3,4,5,6)) cout<<endl;
         GET_MIX_MATRIX("~e-","SELMIX",(1,2,3,4,5,6),(1,2,3,4,5,6)) cout<<endl;
         GET_MIX_MATRIX("h0","SCALARMIX",(1,2),(1,2)) cout<<endl;
         GET_MIX_MATRIX("~chi0","NMIX",(1,2,3,4),(1,2,3,4)) cout<<endl;
         GET_MIX_MATRIX("H+","CHARGEMIX",(1,2),(1,2)) cout<<endl;
         GET_MIX_MATRIX("~u","USQMIX",(1,2,3,4,5,6),(1,2,3,4,5,6)) cout<<endl;
         GET_MIX_MATRIX("~nu","SNUMIX",(1,2,3),(1,2,3)) cout<<endl;

         cout<<endl;
         cout << "Next up: running parameters" << endl;
         cout << "These are all given in the DRbar scheme, at least when running FlexibleSUSY or SoftSUSY. ";
         cout << "There may be some switching or converting once other spectrum generator are added." << endl;
         cout<<endl;
         cout << "Spectrum object running parameters are currently defined at scale Q="
              << spec->GetScale() << " [GeV]" << endl << endl;
         cout<<endl;
         cout << "-- Dimensionless parameters --" <<endl;
         cout << endl << "Gauge couplings:" << endl << endl;
         ECHO(  spec->get_dimensionless_parameter("g1")  )  // U_Y(1) gauge coupling in SU(5) normalisation
         ECHO(  slhaea.at("GAUGE").at(1).at(1)  ) // This is in the Standard Model normalisation as per SLHA conventions
         cout << "Note: " << spec->get_dimensionless_parameter("g1") << " * sqrt(3/5) = "
                          << spec->get_dimensionless_parameter("g1")*sqrt(3./5.) << endl;
         cout<<endl;
         ECHO(  spec->get_dimensionless_parameter("g2")  )  // SU(2) gauge coupling
         ECHO(  slhaea.at("GAUGE").at(2).at(1)  )
         cout<<endl;
         ECHO(  spec->get_dimensionless_parameter("g3")  )  // SU(3) gauge coupling
         ECHO(  slhaea.at("GAUGE").at(3).at(1)  )
         cout<<endl;

         cout << endl << "Yukawa matrices:" << endl << endl;

         // Note, currently we are not using a matrix object or any such thing, so you have to
         // extract the elements of the matrix one at a time. It would probably be handy to
         // add such a return type though.

         #define GET_MATRIX_EL(r, PRODUCT)                                       \
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
           try{                                                                  \
            std::ostringstream echo1;                                            \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_dimensionless_parameter("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_dimensionless_parameter(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           } catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
         }

         #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)))

         GET_MATRIX("Yu","YU",(1,2,3),(1,2,3)) cout << endl;
         GET_MATRIX("Yd","YD",(1,2,3),(1,2,3)) cout << endl;
         GET_MATRIX("Ye","YE",(1,2,3),(1,2,3)) cout << endl;

         // Mass dimension 1 parameters

         cout<<endl;
         cout<<"MSSM mass dimension 1 running parameters"<<endl;
         cout<<endl;
         ECHO(  spec->get_mass_parameter("M1")    )  // Gaugino mass parameter "MassB"
         ECHO(  slhaea.at("MSOFT").at(1).at(1)                )
         cout<<endl;
         ECHO(  spec->get_mass_parameter("M2")    )  // Gaugino mass parameter "MassWB"
         ECHO(  slhaea.at("MSOFT").at(2).at(1)                )
         cout<<endl;
         ECHO(  spec->get_mass_parameter("M3")    )  // Gaugino mass parameter "MassG"
         ECHO(  slhaea.at("MSOFT").at(3).at(1)                )
         cout<<endl;
         ECHO(  spec->get_mass_parameter("Mu")    )  // Superpotential mu parameter
         ECHO(  slhaea.at("HMIX").at(1).at(1)                 )
         cout<<endl;
         ECHO(  spec->get_mass_parameter("vd")    )  // Down-type Higgs vev
         ECHO(  slhaea.at("HMIX").at(102).at(1)               )
         cout<<endl;
         ECHO(  spec->get_mass_parameter("vu")    )  // Up-type Higgs vev
         ECHO(  slhaea.at("HMIX").at(103).at(1)               )
         cout<<endl;

         // Matrices

         #define GET_M1_MATRIX_EL(r, PRODUCT)                                    \
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
           try{                                                                  \
            std::ostringstream echo1;                                            \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_mass_parameter("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_mass_parameter(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           } catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
         }

         #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)))

         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;

         // Seem to be the trilinears, and TYu and au etc. seem to be equal. Ask Peter...

         GET_M1_MATRIX("TYu","TU",(1,2,3),(1,2,3)) cout << endl;
         GET_M1_MATRIX("TYd","TD",(1,2,3),(1,2,3)) cout << endl;
         GET_M1_MATRIX("TYe","TE",(1,2,3),(1,2,3)) cout << endl;
         cout << endl;
         GET_M1_MATRIX("au","TU",(1,2,3),(1,2,3)) cout << endl;
         GET_M1_MATRIX("ad","TD",(1,2,3),(1,2,3)) cout << endl;
         GET_M1_MATRIX("ae","TE",(1,2,3),(1,2,3)) cout << endl;

         // Mass dimension 2 parameters

         cout<<endl;
         cout<<"MSSM mass dimension 2 running parameters"<<endl;
         cout<<endl;
         ECHO(  spec->get_mass2_parameter("mHd2")  )  // Down-type Higgs soft mass
         ECHO(  slhaea.at("MSOFT").at(21).at(1)               )
         cout<<endl;
         ECHO(  spec->get_mass2_parameter("mHu2")  )  // Up-type Higgs soft mass
         ECHO(  slhaea.at("MSOFT").at(22).at(1)               )
         cout<<endl;
         ECHO(  spec->get_mass2_parameter("BMu")   )  // Higgs bilinear soft parameter
         ECHO(  slhaea.at("HMIX").at(101).at(1)               )
         cout<<endl;

         // Matrices

         #define GET_M2_MATRIX_EL(r, PRODUCT)                                    \
         {                                                                       \
            str label = BOOST_PP_SEQ_ELEM(0, PRODUCT);                           \
            str block = BOOST_PP_SEQ_ELEM(1, PRODUCT);                           \
            int i = BOOST_PP_SEQ_ELEM(2, PRODUCT);                               \
            int j = BOOST_PP_SEQ_ELEM(3, PRODUCT);                               \
            std::ostringstream echo1;                                            \
           try {                                                                 \
            std::ostringstream echo2;                                            \
            echo1 <<     "  spec->get_mass2_parameter("<<label<<","<<i<<","<<j<<") = "; \
            double value1 = spec->get_mass2_parameter(label,i,j); \
            echo2 <<     "  SLHAea::to<double>( slhaea.at("<<block<<").at("<<i<<","<<j<<").at(2) ) = "; \
            double value2 = SLHAea::to<double>( slhaea.at(block).at(i,j).at(2) ); \
            cout << echo1.str() << value1 <<endl;                                \
            cout << echo2.str() << value2 <<endl;                                \
            cout << endl;                                                        \
           } catch (const std::exception& e)                                       \
           { add_error(report,e,label+": "+block); }                             \
        }

         #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)))

         cout << endl << "Mass matrices:" << endl << endl;

         GET_M2_MATRIX("mq2","MSQ2",(1,2,3),(1,2,3)) cout << endl;
         GET_M2_MATRIX("mu2","MSU2",(1,2,3),(1,2,3)) cout << endl;
         GET_M2_MATRIX("md2","MSD2",(1,2,3),(1,2,3)) cout << endl;
         GET_M2_MATRIX("me2","MSE2",(1,2,3),(1,2,3)) cout << endl;
         GET_M2_MATRIX("ml2","MSL2",(1,2,3),(1,2,3)) cout << endl;

         cout << endl;

         /// NEW! Tests of override setters

         cout << "Testing set_override functions" << endl;

         cout << "Original M1:" << spec->get_mass_parameter("M1") << endl;
         spec->set_override_mass_parameter(-666,"M1");
         cout << "Override M1:" << spec->get_mass_parameter("M1") << endl;

         cout << "Original ~e-(1):" << spec->get_Pole_Mass("~e-",1) << endl;
         spec->set_override_Pole_Mass(-667,"~e-",1);
         cout << "Override ~e-(1):" << spec->get_Pole_Mass("~e-",1) << endl;

         cout << "Original ml2(1,1):" << spec->get_mass2_parameter("ml2",1,1) << endl;
         spec->set_override_mass2_parameter(-668,"ml2",1,1);
         cout << "Override ml2(1,1):" << spec->get_mass2_parameter("ml2",1,1) << endl;

         /// Now add some entry that didn't exist before
         cout << "has 'new_entry'? " << spec->has_mass_parameter("new_entry") << endl;
         cout << "..." << endl;
         /// Note: if we try to do it like this, it should fail:
         //spec->set_override_mass2_parameter(-1234,"new_entry"); // incorrect: safety still on
         spec->set_override_mass_parameter(-1234,"new_entry",false); // correct: safety check turned off
         cout << "has 'new_entry'? " << spec->has_mass_parameter("new_entry") << endl;
         cout << "new_entry = " << spec->get_mass_parameter("new_entry") << endl;
         cout << endl;


         // Things left to add to demo:
         // - Tree-level masses
         // - Standard Model mixings (need to be added to MSSMSpec as well)
         // - some extra stuff to do regarding organising Standard Model input
         //   parameters (and computing/transferring the pole masses)

         cout << "Test report:" << std::endl << report.str();

         SpecBit_warning().raise(LOCAL_INFO,"\n *** Stopped on purpose to examine spectrum contents ***");
         result = 0;
      }
      //else if(spinfo.find(k3) != spinfo.end())
      //{
      //   std::cout << "Bad spectrum: " << spinfo.at(3) << std::endl;
      //   result = 1;
      //}
      else if(spinfo.find(k4) != spinfo.end())
      {
         std::cout << "Bad spectrum: " << spinfo.at(4) << std::endl;
         result = 1;
      }
      else
      {
         SpecBit_error().raise(LOCAL_INFO,"Bug in exampleRead 'if' logic");
      }
   }

  }
}

Updated on 2024-07-18 at 13:53:32 +0000