StandardModel.h

Go to the documentation of this file.

/*
 * Copyright (C) 2012 HEPfit Collaboration
 *
 *
 * For the licensing terms see doc/COPYING.
 */
 
#ifndef STANDARDMODEL_H
#define STANDARDMODEL_H
 
#include "Flavour.h"
#include "QCD.h"
#include "CKM.h"
#include "PMNS.h"
#include "StandardModelMatching.h"
#include "Matching.h"
#include <gsl/gsl_integration.h>
 
class EWSMcache;
class EWSMOneLoopEW;
class EWSMTwoLoopQCD;
class EWSMTwoLoopEW;
class EWSMThreeLoopQCD;
class EWSMThreeLoopEW2QCD;
class EWSMThreeLoopEW;
class EWSMApproximateFormulae;
class LeptonFlavour;
 
 
class StandardModel : public QCD {
public:
 
 
    // Radiative Corrections for the LEP-II observables
    enum LEP2RCs {
        Weak = 0, 
        WeakBox, 
        ISR, 
        QEDFSR, 
        QCDFSR, 
        NUMofLEP2RCs
    };
    
    enum orders_EW {
        EW1 = 0, 
        EW1QCD1, 
        EW1QCD2, 
        EW2, 
        EW2QCD1, 
        EW3, 
        orders_EW_size 
    };
 
 
    
//    static const int NSMvars = 38; ///< The number of the model parameters in %StandardModel. !!! PMNS INCLUDED
    static const int NSMvars = 26; 
 
    static std::string SMvars[NSMvars];
 
    static const double GeVminus2_to_nb;
    
    static const double Mw_error;
 
    StandardModel();
 
    virtual ~StandardModel();
 
 
    // Initialization
 
    virtual bool InitializeModel();
 
 
    // Model parameters
 
    virtual bool Init(const std::map<std::string, double>& DPars);
 
    virtual bool PreUpdate();
 
    virtual bool Update(const std::map<std::string, double>& DPars);
 
    virtual bool PostUpdate();
    
    int getIterationNo() const
    {
        return iterationNo;
    }
 
    virtual bool CheckParameters(const std::map<std::string, double>& DPars);
 
 
    // Flags
 
    virtual bool setFlag(const std::string name, const bool value);
 
    virtual bool setFlagStr(const std::string name, const std::string value);
 
    virtual bool CheckFlags() const;
 
    bool IsFlagWithoutNonUniversalVC() const
    {
        return FlagWithoutNonUniversalVC;
    }
 
    bool IsFlagNoApproximateGammaZ() const
    {
        return FlagNoApproximateGammaZ;
    }
 
    void setFlagNoApproximateGammaZ(bool FlagNoApproximateGammaZ)
    {
        this->FlagNoApproximateGammaZ = FlagNoApproximateGammaZ;
    }
 
    std::string getFlagMw() const
    {
        return FlagMw;
    }
 
    std::string getFlagRhoZ() const
    {
        return FlagRhoZ;
    }
 
    std::string getFlagKappaZ() const
    {
        return FlagKappaZ;
    }
 
    void setFlagCacheInStandardModel(bool FlagCacheInStandardModel)
    {
        this->FlagCacheInStandardModel = FlagCacheInStandardModel;
    }
 
 
    // get and set methods for class members
 
    Particle getLeptons(const QCD::lepton p) const
    {
        return leptons[p];
    }
 
    double getMz() const
    {
        return Mz;
    }
 
    double getAlsMz() const
    {
        return AlsMz;
    }
 
    double getGF() const
    {
        return GF;
    }
 
    double getAle() const
    {
        return ale;
    }
 
    double getDAle5Mz() const
    {
        return dAle5Mz;
    }
 
    virtual double getMHl() const
    {
        return mHl;
    }
 
    double getDelMw() const
    {
        return delMw;
    }
 
    double getDelSin2th_l() const
    {
        return delSin2th_l;
    }
    
    double getDelSin2th_q() const
    {
        return delSin2th_q;
    }
    
    double getDelSin2th_b() const
    {
        return delSin2th_b;
    }
 
    double getDelGammaZ() const
    {
        return delGammaZ;
    }
    
    double getDelSigma0H() const
    {
        return delsigma0H;
    }
    
    double getDelR0l() const
    {
        return delR0l;
    }
    
    double getDelR0c() const
    {
        return delR0c;
    }
    
    double getDelR0b() const
    {
        return delR0b;
    }
 
    gslpp::matrix<gslpp::complex> getVCKM() const // why don't we return a const reference?
    {
        return myCKM.getCKM();
    }
 
    CKM getCKM() const
    {
        return myCKM;
    }
 
    gslpp::matrix<gslpp::complex> getUPMNS() const
    {
        return myPMNS.getPMNS();
    }
 
    gslpp::matrix<gslpp::complex> getYu() const
    {
        return Yu;
    }
 
    gslpp::matrix<gslpp::complex> getYd() const
    {
        return Yd;
    }
 
    gslpp::matrix<gslpp::complex> getYn() const
    {
        return Yn;
    }
 
    gslpp::matrix<gslpp::complex> getYe() const
    {
        return Ye;
    }
 
    double getMuw() const
    {
        return muw;
    }
 
    virtual StandardModel getTrueSM() const
    {
        throw std::runtime_error("StandardModel::getTrueSM() must be overridden by the NP extension.");
    }
 
    virtual StandardModelMatching& getMatching() const
    {
        return SMM.getObj();
    }
 
    EWSMcache* getMyEWSMcache() const
    {
        return myEWSMcache;
    }
 
    EWSMOneLoopEW* getMyOneLoopEW() const
    {
        return myOneLoopEW;
    }
 
    EWSMApproximateFormulae* getMyApproximateFormulae() const
    {
        return myApproximateFormulae;
    }
 
    EWSMThreeLoopEW* getMyThreeLoopEW() const
    {
        return myThreeLoopEW;
    }
 
    EWSMThreeLoopEW2QCD* getMyThreeLoopEW2QCD() const
    {
        return myThreeLoopEW2QCD;
    }
 
    EWSMThreeLoopQCD* getMyThreeLoopQCD() const
    {
        return myThreeLoopQCD;
    }
 
    EWSMTwoLoopEW* getMyTwoLoopEW() const
    {
        return myTwoLoopEW;
    }
 
    EWSMTwoLoopQCD* getMyTwoLoopQCD() const
    {
        return myTwoLoopQCD;
    }
 
    const Flavour& getFlavour() const
    {
        return SMFlavour;
    }
 
    LeptonFlavour* getMyLeptonFlavour() const
    {
        return myLeptonFlavour;
    }
    // QED coupling
 
    double ale_OS(const double mu, orders order = FULLNLO) const;
 
    double Beta_s(int nm, unsigned int nf) const;
    
    double Beta_e(int nm, unsigned int nf) const;
 
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Woverloaded-virtual"
#endif
    double Als(double mu, orders order = FULLNLO, bool qed_flag = false, bool Nf_thr = true) const;
    double AlsByOrder(double mu, orders order = FULLNLO, bool qed_flag = false, bool Nf_thr = true) const;
#ifdef __clang__
#pragma clang diagnostic pop
#endif
 
    double Ale(double mu, orders order, bool Nf_thr = true) const;
 
    double DeltaAlphaLepton(const double s) const;
 
    double DeltaAlphaL5q() const;
 
    double DeltaAlphaTop(const double s) const;
 
    double DeltaAlpha() const;
 
    double alphaMz() const;
 
    double Alstilde5(const double mu) const;
    
    virtual double getCCC1() const { return 0.; };
    
    virtual double getCCC2() const { return 0.; };
    
    virtual double getCCC3() const { return 0.; };
    
    virtual double getCCC4() const { return 0.; };
    
    virtual double getCCC5() const { return 0.; };
 
 
 
    // Higgs VEV
 
    virtual double v() const;
 
 
    // The W-boson mass
 
    virtual double Mw_tree() const;
 
    double s02() const;
 
    double c02() const;
 
    virtual double Mw() const;
 
    virtual double cW2(const double Mw_i) const;
    virtual double cW2() const;
 
    virtual double sW2(const double Mw_i) const;
    double sW2() const;
 
    virtual double DeltaR() const;
 
    void ComputeDeltaRho(const double Mw_i, double DeltaRho[orders_EW_size]) const;
 
    void ComputeDeltaR_rem(const double Mw_i, double DeltaR_rem[orders_EW_size]) const;
 
 
    // The W and Z masses in the complex-pole/fixed-width scheme
 
    double Mzbar() const;
 
    double MwbarFromMw(const double Mw) const;
 
    double MwFromMwbar(const double Mwbar) const;
 
    virtual double DeltaRbar() const;
 
 
    // The W-boson decay width
 
    virtual double rho_GammaW(const Particle fi, const Particle fj) const;
 
    virtual double GammaW(const Particle fi, const Particle fj) const;
 
    virtual double GammaW() const;
 
 
    // EWPO at Z-pole
 
    virtual double A_f(const Particle f) const;
 
    virtual double AFB(const Particle f) const;
 
    virtual double sin2thetaEff(const Particle f) const;
 
    virtual double GammaZ(const Particle f) const;
 
    virtual double Gamma_inv() const;
 
    virtual double Gamma_had() const;
 
    virtual double Gamma_Z() const;
 
    virtual double sigma0_had() const;
 
    virtual double R0_f(const Particle f) const;
    
    virtual double R_inv() const;
    
    virtual double N_nu() const;
 
 
    // Zff effective couplings
 
    virtual gslpp::complex gV_f(const Particle f) const;
 
    virtual gslpp::complex gA_f(const Particle f) const;
 
    virtual gslpp::complex rhoZ_f(const Particle f) const;
 
    virtual gslpp::complex kappaZ_f(const Particle f) const;
 
    virtual gslpp::complex deltaRhoZ_f(const Particle f) const;
 
    virtual gslpp::complex deltaKappaZ_f(const Particle f) const;
 
 
    // Epsilon parameters for EWPO
 
    virtual double epsilon1() const;
 
    virtual double epsilon2() const;
 
    virtual double epsilon3() const;
 
    virtual double epsilonb() const;
 
 
    // For EWPO caches
 
    static const int NumSMParamsForEWPO = 33;
 
    bool checkSMparamsForEWPO();
 
    // Several Higgs-related quantities used in Higgs coupling analysis
 
    double computeSigmaggH(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 16.83; // in pb for Mh=125.1 GeV
        } else if (sqrt_s == 8.0) {
            return 21.40; // in pb for Mh=125.1 GeV
        } else if (sqrt_s == 13.0) {
            return 48.61; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 14.0) {
            return 54.72; // in pb for Mh=125.09 GeV     
        } else if (sqrt_s == 27.0) {
            return 146.65; // in pb for Mh=125.09 GeV
        } else if (sqrt_s == 100.0) {
            return 740.3; // in pb for Mh=125. GeV            
        } else if (sqrt_s == 1.96) {
            return 0.9493; // in pb for Mh=125 GeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaggH()");
    }
 
    double computeSigmaggH_tt(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 16.69; // in pb for Mh=125.09 GeV 
        } else if (sqrt_s == 8.0) {
            return 21.20; // in pb for Mh=125.09 GeV
        } else if (sqrt_s == 13.0) {
            return 47.94; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 14.0) {
            return 53.93; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 27.0) {
            return computeSigmaggH(sqrt_s) / computeSigmaggH(14.) * computeSigmaggH_tt(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else if (sqrt_s == 100.0) {
            return computeSigmaggH(sqrt_s) / computeSigmaggH(14.) * computeSigmaggH_tt(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaggH_tt()");
    }
 
    double computeSigmaggH_bb(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 0.04; // in pb for Mh=125.09 GeV 
        } else if (sqrt_s == 8.0) {
            return 0.05; // in pb for Mh=125.09 GeV 
        } else if (sqrt_s == 13.0) {
            return 0.10; // in pb for Mh=125.09 GeV             
        } else if (sqrt_s == 14.0) {
            return 0.11; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 27.0) {
            return computeSigmaggH(sqrt_s) / computeSigmaggH(14.) * computeSigmaggH_bb(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else if (sqrt_s == 100.0) {
            return computeSigmaggH(sqrt_s) / computeSigmaggH(14.) * computeSigmaggH_bb(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaggH_bb()");
    }
 
    double computeSigmaggH_tb(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return -0.66; // in pb for Mh=125.09 GeV 
        } else if (sqrt_s == 8.0) {
            return -0.82; // in pb for Mh=125.09 GeV 
        } else if (sqrt_s == 13.0) {
            return -1.73; // in pb for Mh=125.09 GeV 
        } else if (sqrt_s == 14.0) {
            return -1.92; // in pb for Mh=125.09 GeV 
        } else if (sqrt_s == 27.0) {
            return computeSigmaggH(sqrt_s) / computeSigmaggH(14.) * computeSigmaggH_tb(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else if (sqrt_s == 100.0) {
            return computeSigmaggH(sqrt_s) / computeSigmaggH(14.) * computeSigmaggH_tb(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaggH_tb()");
    }
 
    double computeSigmaVBF(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 1.241; // in pb for Mh=125.09 GeV
        } else if (sqrt_s == 8.0) {
            return 1.601; // in pb for Mh=125.09 GeV
        } else if (sqrt_s == 13.0) {
            return 3.766; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 14.0) {
            return 4.260; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 27.0) {
            return 11.838; // in pb for Mh=125.09 GeV
        } else if (sqrt_s == 100.0) {
            return 82.0; // in pb for Mh=125. GeV
        } else if (sqrt_s == 1.96) {
            return 0.0653; // in pb for Mh=125 GeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaVBF()");
    }
 
    double computeSigmaWF(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 0.946; // in pb for Mh=125 GeV
        } else if (sqrt_s == 8.0) {
            return 1.220; // in pb for Mh=125 GeV
        } else if (sqrt_s == 13.0) {
            return 2.882; // in pb for Mh=125 GeV            
        } else if (sqrt_s == 14.0) {
            return 3.260; // in pb for Mh=125 GeV            
        } else if (sqrt_s == 27.0) {
            return computeSigmaVBF(sqrt_s) / computeSigmaVBF(14.) * computeSigmaWF(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else if (sqrt_s == 100.0) {
            return computeSigmaVBF(sqrt_s) / computeSigmaVBF(14.) * computeSigmaWF(14.); // in the absence of this value we rescale the LHC result at 14 TeV 
        } else if (sqrt_s == 1.96) {
            return computeSigmaVBF(sqrt_s) / computeSigmaVBF(7.) * computeSigmaWF(7.); // in the absence of individual cross sections for TeVatron we rescale the LHC ones
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaWF()");
    }
 
    double computeSigmaZF(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 0.333; // in pb for Mh=125 GeV
        } else if (sqrt_s == 8.0) {
            return 0.432; // in pb for Mh=125 GeV
        } else if (sqrt_s == 13.0) {
            return 1.049; // in pb for Mh=125 GeV
        } else if (sqrt_s == 14.0) {
            return 1.191; // in pb for Mh=125 GeV
        } else if (sqrt_s == 27.0) {
            return computeSigmaVBF(sqrt_s) / computeSigmaVBF(14.) * computeSigmaZF(14.); // in the absence of this value we rescale the LHC result at 14 TeV
        } else if (sqrt_s == 100.0) {
            return computeSigmaVBF(sqrt_s) / computeSigmaVBF(14.) * computeSigmaZF(14.); // in the absence of this value we rescale the LHC result at 14 TeV 
        } else if (sqrt_s == 1.96) {
            return computeSigmaVBF(sqrt_s) / computeSigmaVBF(7.) * computeSigmaZF(7.); // in the absence of individual cross sections for TeVatron we rescale the LHC ones
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaZF()");
    }
 
    double computeSigmaZWF(const double sqrt_s) const
    {
        return 0.;
    }
 
    double computeSigmaWH(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 0.577; // in pb for Mh=125.1 GeV
            //return 0.5688; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 8.0) {
            return 0.7027; // in pb for Mh=125.1 GeV
            //return 0.6931; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 13.0) {
            return 1.358; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 14.0) {
            return 1.498; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 27.0) {
            return 3.397; // in pb for Mh=125.09 GeV
        } else if (sqrt_s == 100.0) {
            return 15.9; // in pb for Mh=125. GeV
        } else if (sqrt_s == 1.96) {
            return 0.1295; // in pb for Mh=125 GeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaWH()");
    }
 
    double computeSigmaZH(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 0.3341; // in pb for Mh=125.1 GeV
            //return 0.3299; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 8.0) {
            return 0.4142; // in pb for Mh=125.1 GeV
            //return 0.4091; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 13.0) {
            return 0.880; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 14.0) {
            return 0.981; // in pb for Mh=125.09 GeV   
        } else if (sqrt_s == 27.0) {
            return 2.463; // in pb for Mh=125.09 GeV        
        } else if (sqrt_s == 100.0) {
            return 11.26; // in pb for Mh=125. GeV
        } else if (sqrt_s == 1.96) {
            return 0.0785; // in pb for Mh=125 GeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmaZH()");
    }
 
    double computeSigmattH(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 0.0861; // in pb for Mh=125.1 GeV
            //return 0.0851; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 8.0) {
            return 0.129; // in pb for Mh=125.1 GeV
            //return 0.1274; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 13.0) {
            return 0.5060; // in pb for Mh=125.1 GeV            
        } else if (sqrt_s == 14.0) {
            return 0.6128; // in pb for Mh=125.09 GeV            
        } else if (sqrt_s == 27.0) {
            return 2.86; // in pb for Mh=125.09 GeV
        } else if (sqrt_s == 100.0) {
            return 37.9; // in pb for Mh=125. GeV
        } else if (sqrt_s == 1.96) {
            return 0.0043; // in pb for Mh=125 GeV
        } else
            throw std::runtime_error("Bad argument in StandardModel::computeSigmattH()");
    }
 
    double computeBrHtogg() const
    {
        return 8.179e-2; // Mh=125.1 GeV
    }
 
    double computeBrHtoWW() const
    {
        //return 2.23e-1; // Mh=125.5 GeV
        return 2.154e-1; // Mh=125.1 GeV
    }
 
    double computeBrHtoZZ() const
    {
        return 2.643e-2; // Mh=125.1 GeV
        //return 2.79e-2; // Mh=125.6 GeV
    }
    
    double computeBrHtoZZinv() const
    {
        return 1.06e-3;
    }
 
    double computeBrHtoZga() const
    {
        return 1.541e-3; // Mh=125.1 GeV
        //return 1.59e-3; // Mh=125.6 GeV
    }
 
    double computeBrHtogaga() const
    {
        return 2.27e-3; // Mh=125.1 GeV
    }
    
    double computeBrHtomumu() const
    {
        return 2.17e-4; // Mh=125.1 GeV
    }
 
    double computeBrHtotautau() const
    {
        return 6.256e-2; // Mh=125.1 GeV
        //return 6.22e-2; // Mh=125.6 GeV
    }
 
    double computeBrHtocc() const
    {
        return 2.883e-2; // Mh=125.1 GeV
        //return 2.86e-2; // Mh=125.6 GeV
    }
    
    double computeBrHtoss() const
    {
        return 4.0e-4;
    }
 
    double computeBrHtobb() const
    {
        return 5.807e-1; // Mh=125.1 GeV
        //return 5.67e-1; // Mh=125.6 GeV
    }
 
    double computeGammaHTotal() const
    {
        return 4.101e-3; // Mh=125.1 GeV
        //return 4.15e-3; // Mh=125.6 GeV
    }
 
    double computeGammaHgg_tt() const
    {
        return 380.8; // in keV for Mh=125 GeV
        //return 389.6; // in keV for Mh=126 GeV
    }
 
    double computeGammaHgg_bb() const
    {
        return 3.96; // in keV for Mh=125 GeV
        //return 3.95; // in keV for Mh=126 GeV
    }
 
    double computeGammaHgg_tb() const
    {
        return -42.1; // in keV for Mh=125 GeV
        //return -42.7; // in keV for Mh=126 GeV
    }
 
    double computeGammaHZga_tt() const
    {
        return 21.74; // in eV for Mh=125 GeV
        //return 23.51; // in eV for Mh=126 GeV
    }
 
    double computeGammaHZga_WW() const
    {
        return 7005.6; // in eV for Mh=125 GeV
        //return 7648.4; // in eV for Mh=126 GeV
    }
 
    double computeGammaHZga_tW() const
    {
        return -780.4; // in eV for Mh=125 GeV
        //return -848.1; // in eV for Mh=126 GeV
    }
 
    double computeGammaHgaga_tt() const
    {
        return 662.84; // in eV for Mh=125 GeV
        //return 680.39; // in eV for Mh=126 GeV
    }
 
    double computeGammaHgaga_WW() const
    {
        return 14731.86; // in eV for Mh=125 GeV
        //return 15221.98; // in eV for Mh=126 GeV
    }
 
    double computeGammaHgaga_tW() const
    {
        return -6249.93; // in eV for Mh=125 GeV
        //return -6436.35; // in eV for Mh=126 GeV
    }
 
    virtual double getCBd() const
    {
        return 1.;
    }
 
    virtual double getCBs() const
    {
        return 1.;
    }
 
    virtual double getCDMK() const
    {
        return 1.;
    }
 
    virtual double getCepsK() const
    {
        return 1.;
    }
 
    virtual double getPhiBs() const
    {
        return 0.;
    }
 
    virtual double getPhiBd() const
    {
        return 0.;
    }
 
    virtual double getmq(const QCD::quark q, const double mu) const
    {
        return m_q(q, mu, FULLNLO);
    }
protected:
 
    virtual void setParameter(const std::string name, const double& value);
 
    virtual void computeCKM();
 
    virtual void computeYukawas();
 
    Particle leptons[6]; 
    CKM myCKM; 
    PMNS myPMNS;
//    gslpp::matrix<gslpp::complex> VCKM; ///< The %CKM matrix.
//    gslpp::matrix<gslpp::complex> UPMNS; ///<  The %PMNS matrix.
    gslpp::matrix<gslpp::complex> Yu; 
    gslpp::matrix<gslpp::complex> Yd; 
    gslpp::matrix<gslpp::complex> Yn; 
    gslpp::matrix<gslpp::complex> Ye; 
 
    mutable Matching<StandardModelMatching,StandardModel> SMM; 
    
    // model parameters
    double AlsMz; 
    double Mz; 
    double GF; 
    double ale; 
    double dAle5Mz; 
    double mHl; 
    double delMw; 
    double delSin2th_l; 
    double delSin2th_q; 
    double delSin2th_b; 
    double delGammaZ; 
    double delsigma0H; 
    double delR0l; 
    double delR0c; 
    double delR0b; 
    double lambda; 
    double A; 
    double rhob; 
    double etab; 
    double Vus; 
    double Vcb; 
    double Vub; 
    double gamma; 
    double muw; 
    double s12, s13, s23, delta, alpha21, alpha31;
 
 
    // For EWPO
 
    bool flag_order[orders_EW_size];
 
    double SchemeToDouble(const std::string scheme) const
    {
        if (scheme.compare("NORESUM") == 0)
            return 0.0;
        else if (scheme.compare("OMSI") == 0)
            return 1.0;
        else if (scheme.compare("INTERMEDIATE") == 0)
            return 2.0;
        else if (scheme.compare("OMSII") == 0)
            return 3.0;
        else if (scheme.compare("APPROXIMATEFORMULA") == 0)
            return 4.0;
        else
            throw std::runtime_error("EWSM::SchemeToDouble: bad scheme");
    }
 
    bool checkEWPOscheme(const std::string scheme) const
    {
        if (scheme.compare("NORESUM") == 0
                || scheme.compare("OMSI") == 0
                || scheme.compare("INTERMEDIATE") == 0
                || scheme.compare("OMSII") == 0
                || scheme.compare("APPROXIMATEFORMULA") == 0)
            return true;
        else
            return false;
    }
    
    
    double m_q(const QCD::quark q, const double mu, const orders order=FULLNLO) const 
    {
        switch(q) {
            case QCD::UP:
            case QCD::DOWN:
            case QCD::STRANGE:
                return Mrun(mu, getQuarks(q).getMass_scale(), 
                                         getQuarks(q).getMass(), order);
            case QCD::CHARM:
            case QCD::BOTTOM:
                return Mrun(mu, getQuarks(q).getMass(), order);
            case QCD::TOP:
                return getMtpole(); // the pole mass
            default:
                throw std::runtime_error("Error in StandardModel::m_q()"); 
        }
    }
 
    double resumMw(const double Mw_i, const double DeltaRho[orders_EW_size],
            const double DeltaR_rem[orders_EW_size]) const;
 
    double resumRhoZ(const double DeltaRho[orders_EW_size],
            const double deltaRho_rem[orders_EW_size],
            const double DeltaRbar_rem, const bool bool_Zbb) const;
 
    double resumKappaZ(const double DeltaRho[orders_EW_size],
            const double deltaKappa_rem[orders_EW_size],
            const double DeltaRbar_rem, const bool bool_Zbb) const;
 
    double taub() const;
 
    double Delta_EWQCD(const QCD::quark q) const;
 
    double RVq(const QCD::quark q) const;
 
    double RAq(const QCD::quark q) const;
 
    double RVh() const;
 
    bool requireCKM; 
    bool requireYe; 
    bool requireYn; 
 
    Flavour SMFlavour; 
        
private:
    EWSMcache* myEWSMcache; 
    EWSMOneLoopEW* myOneLoopEW; 
    EWSMTwoLoopQCD* myTwoLoopQCD; 
    EWSMThreeLoopQCD* myThreeLoopQCD; 
    EWSMTwoLoopEW* myTwoLoopEW; 
    EWSMThreeLoopEW2QCD* myThreeLoopEW2QCD; 
    EWSMThreeLoopEW* myThreeLoopEW; 
    EWSMApproximateFormulae* myApproximateFormulae; 
    LeptonFlavour* myLeptonFlavour; 
 
    bool FlagWithoutNonUniversalVC; 
    bool FlagNoApproximateGammaZ; 
    std::string FlagMw; 
    std::string FlagRhoZ; 
    std::string FlagKappaZ; 
    bool FlagWolfenstein; 
    
    bool FlagSMAux; 
 
    // Caches for EWPO
 
    bool FlagCacheInStandardModel; 
    mutable double SMparamsForEWPO_cache[NumSMParamsForEWPO];
    mutable double DeltaAlphaLepton_cache; 
    mutable double DeltaAlpha_cache; 
    mutable double Mw_cache; 
    mutable double GammaW_cache; 
    mutable gslpp::complex rhoZ_f_cache[12]; 
    mutable gslpp::complex kappaZ_f_cache[12]; 
    mutable bool useDeltaAlphaLepton_cache;
    mutable bool useDeltaAlpha_cache;
    mutable bool useMw_cache;
    mutable bool useGammaW_cache;
    mutable bool useRhoZ_f_cache[12];
    mutable bool useKappaZ_f_cache[12];
    
    int iterationNo;
    
    double AlsWithInit(double mu, double alsi, double mu_i, orders order, bool qed_flag) const;
    double AleWithInit(double mu, double alsi, double mu_i, orders order) const;
    static const int CacheSize = 5; 
    mutable double als_cache[11][CacheSize]; 
    mutable double ale_cache[10][CacheSize]; 
    mutable orders realorder;
 
 
};
 
#endif  /* STANDARDMODEL_H */