StandardModel.h

Go to the documentation of this file.

/*
 * Copyright (C) 2012 HEPfit Collaboration
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 */

#ifndef STANDARDMODEL_H
#define STANDARDMODEL_H

#include <gslpp.h>
#include "QCD.h"
#include "CKM.h"
#include "WilsonCoefficient.h"
#include "StandardModelMatching.h"

class EWSMcache;
class EWSMOneLoopEW;
class EWSMTwoLoopQCD;
class EWSMTwoLoopEW;
class EWSMThreeLoopQCD;
class EWSMThreeLoopEW2QCD;
class EWSMThreeLoopEW;
class EWSMApproximateFormulae;
class Flavour;
class LeptonFlavour;
class EWSMTwoFermionsLEP2;
class StandardModel : public QCD {
public:

    enum lepton {
        NEUTRINO_1, 
        ELECTRON, 
        NEUTRINO_2, 
        MU, 
        NEUTRINO_3, 
        TAU 
    };

    enum orders_EW {
        EW1 = 0, 
        EW1QCD1, 
        EW1QCD2, 
        EW2, 
        EW2QCD1, 
        EW3, 
        orders_EW_size 
    };

    static const int NSMvars = 26; 

    static const 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 StandardModel::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;
    }

    double getMHl() const
    {
        return mHl;
    }

    double getDelMw() const
    {
        return delMw;
    }

    double getDelSin2th_l() const
    {
        return delSin2th_l;
    }

    double getDelGammaZ() const
    {
        return delGammaZ;
    }

    gslpp::matrix<gslpp::complex> getVCKM() const
    {
        return VCKM;
    }

    CKM getCKM() const
    {
        return myCKM;
    }

    double getLambda() const
    {
        return lambda;
    }

    double getA() const
    {
        return A;
    }

    double getRhob() const
    {
        return rhob;
    }

    double getEtab() const
    {
        return etab;
    }

    gslpp::matrix<gslpp::complex> getUPMNS() const
    {
        return UPMNS;
    }

    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;
    }

    double getEpsK() const
    {
        return EpsK;
    }

    double getphiEpsK() const
    {
        return phiEpsK;
    }

    double getKbarEpsK() const
    {
        return KbarEpsK;
    }

    double getDeltaMK() const
    {
        return DeltaMK;
    }

    double getDmk() const
    {
        return Dmk;
    }

    double getSM_M12D() const
    {
        return SM_M12D;
    }

    virtual StandardModel getTrueSM() const
    {
        throw std::runtime_error("StandardModel::getTrueSM() must be overridden by the NP extension.");
    }

    virtual StandardModelMatching* getMyMatching() const
    {
        return myStandardModelMatching;
    }

    EWSMcache* getMyEWSMcache() const
    {
        return myEWSMcache;
    }

    EWSMOneLoopEW* getMyOneLoopEW() const
    {
        return myOneLoopEW;
    }

    EWSMApproximateFormulae* getMyApproximateFormulae() const
    {
        return myApproximateFormulae;
    }

    EWSMTwoFermionsLEP2* getMyTwoFermionsLEP2() const
    {
        return myTwoFermionsLEP2;
    }
    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;
    }

    Flavour* getMyFlavour() const
    {
        return myFlavour;
    }

    LeptonFlavour* getMyLeptonFlavour() const
    {
        return myLeptonFlavour;
    }


    // CKM parameters

    double computeBeta() const;

    double computeGamma() const;

    double computeAlpha() const;

    double computeBetas() const;

    gslpp::complex computelamt() const;

    gslpp::complex computelamc() const;

    gslpp::complex computelamu() const;

    gslpp::complex computelamt_d() const;

    gslpp::complex computelamc_d() const;

    gslpp::complex computelamu_d() const;

    gslpp::complex computelamt_s() const;

    gslpp::complex computelamc_s() const;

    gslpp::complex computelamu_s() const;

    double computeRt() const;

    double computeRts() const;

    double computeRb() const;


    // QED coupling

    double ale_OS(const double mu, orders order = FULLNLO) const;

    double DeltaAlphaLepton(const double s) const;

    double DeltaAlphaL5q() const;

    double DeltaAlphaTop(const double s) const;

    double DeltaAlpha() const;

    double alphaMz() const;


    // 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;


    // 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 = 27;

    bool checkSMparamsForEWPO();

    // Several Higgs-related quantities used in Higgs coupling analysis

    double computeSigmaggH(const double sqrt_s) const
    {
        if (sqrt_s == 7.0) {
            return 15.11; // in pb for Mh=125.1 GeV
            //return 14.99; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 8.0) {
            return 19.24; // in pb for Mh=125.1 GeV
            //return 19.09; // in pb for Mh=125.6 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 14.355; // in pb for Mh=125 GeV with a scale choice of Mh (corrected from the value in 1307.1347v2)
            return 15.89; // in pb for Mh=125 GeV with a scale choice of Mh/2
            //return 15.76; // in pb for Mh=125.5 GeV with a scale choice of Mh/2
        } else if (sqrt_s == 8.0) {
            //return 18.31; // in pb for Mh=125 GeV with a scale choice of Mh
            return 20.18; // in pb for Mh=125 GeV with a scale choice of Mh/2
            //return 20.02; // in pb for Mh=125.5 GeV with a scale choice of Mh/2
        } 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.09528; // in pb for Mh=125 GeV with a scale choice of Mh
            return 0.1079; // in pb for Mh=125 GeV with a scale choice of Mh/2
            //return 0.1057; // in pb for Mh=125.5 GeV with a scale choice of Mh/2
        } else if (sqrt_s == 8.0) {
            //return 0.1206; // in pb for Mh=125 GeV with a scale choice of Mh
            return 0.1357; // in pb for Mh=125 GeV with a scale choice of Mh/2
            //return 0.1330; // in pb for Mh=125.5 GeV with a scale choice of Mh/2
        } 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.8970; // in pb for Mh=125 GeV with a scale choice of Mh
            return -0.9726; // in pb for Mh=125 GeV with a scale choice of Mh/2
            //return -0.9609; // in pb for Mh=125.5 GeV with a scale choice of Mh/2
        } else if (sqrt_s == 8.0) {
            //return -1.125; // in pb for Mh=125 GeV with a scale choice of Mh
            return -1.206; // in pb for Mh=125 GeV with a scale choice of Mh/2
            //return -1.192; // in pb for Mh=125.5 GeV with a scale choice of Mh/2
        } 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.222; // in pb for Mh=125.1 GeV
            //return 1.214; // in pb for Mh=125.6 GeV
        } else if (sqrt_s == 8.0) {
            return 1.579; // in pb for Mh=125.1 GeV
            //return 1.572; // in pb for Mh=125.6 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.938; // in pb for Mh=125 GeV
        } else if (sqrt_s == 8.0) {
            return 1.210; // in pb for Mh=125 GeV
        } 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.321; // in pb for Mh=125 GeV
        } else if (sqrt_s == 8.0) {
            return 0.417; // in pb for Mh=125 GeV
        } 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 == 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 == 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 == 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.56e-2; // Mh=125.1 GeV
    }

    double computeBrHtoWW() const
    {
        //return 2.23e-1; // Mh=125.5 GeV
        return 2.16e-1; // Mh=125.1 GeV
    }

    double computeBrHtoZZ() const
    {
        return 2.67e-2; // Mh=125.1 GeV
        //return 2.79e-2; // Mh=125.6 GeV
    }

    double computeBrHtoZga() const
    {
        return 1.55e-3; // Mh=125.1 GeV
        //return 1.59e-3; // Mh=125.6 GeV
    }

    double computeBrHtogaga() const
    {
        return 2.28e-3; // Mh=125.1 GeV
    }

    double computeBrHtotautau() const
    {
        return 6.30e-2; // Mh=125.1 GeV
        //return 6.22e-2; // Mh=125.6 GeV
    }

    double computeBrHtocc() const
    {
        return 2.90e-2; // Mh=125.1 GeV
        //return 2.86e-2; // Mh=125.6 GeV
    }

    double computeBrHtobb() const
    {
        return 5.75e-1; // Mh=125.1 GeV
        //return 5.67e-1; // Mh=125.6 GeV
    }

    double computeGammaHTotal() const
    {
        return 4.08e-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
    }

protected:

    virtual void setParameter(const std::string name, const double& value);

    virtual void computeCKM();

    virtual void computeYukawas();

    Particle leptons[6]; 
    CKM myCKM; 
    gslpp::matrix<gslpp::complex> VCKM; 
    gslpp::matrix<gslpp::complex> UPMNS; 
    gslpp::matrix<gslpp::complex> Yu; 
    gslpp::matrix<gslpp::complex> Yd; 
    gslpp::matrix<gslpp::complex> Yn; 
    gslpp::matrix<gslpp::complex> Ye; 

    // model parameters
    double AlsMz; 
    double Mz; 
    double GF; 
    double ale; 
    double dAle5Mz; 
    double mHl; 
    double delMw; 
    double delSin2th_l; 
    double delGammaZ; 
    double lambda; 
    double A; 
    double rhob; 
    double etab; 
    double muw; 

    double EpsK;
    double phiEpsK;
    double DeltaMK;
    double KbarEpsK;
    double Dmk;
    double SM_M12D;


    // 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 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; 


private:

    EWSMcache* myEWSMcache; 
    EWSMOneLoopEW* myOneLoopEW; 
    EWSMTwoLoopQCD* myTwoLoopQCD; 
    EWSMThreeLoopQCD* myThreeLoopQCD; 
    EWSMTwoLoopEW* myTwoLoopEW; 
    EWSMThreeLoopEW2QCD* myThreeLoopEW2QCD; 
    EWSMThreeLoopEW* myThreeLoopEW; 
    EWSMApproximateFormulae* myApproximateFormulae; 
    Flavour* myFlavour; 
    LeptonFlavour* myLeptonFlavour; 

    EWSMTwoFermionsLEP2* myTwoFermionsLEP2; 

    StandardModelMatching* myStandardModelMatching; 

    bool FlagWithoutNonUniversalVC; 
    bool FlagNoApproximateGammaZ; 
    std::string FlagMw; 
    std::string FlagRhoZ; 
    std::string FlagKappaZ; 



    // 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;

};

#endif  /* STANDARDMODEL_H */