GeneralTHDM.h

Go to the documentation of this file.

/* 
 * Copyright (C) 2016 HEPfit Collaboration
 *
 *
 * For the licensing terms see doc/COPYING.
 */
 
#ifndef GENERALTHDM_H
#define GENERALTHDM_H
 
#include "StandardModel.h"
#include "GeneralTHDMMatching.h"
#include "NPbase.h"
#include <gsl/gsl_errno.h>
 
 
class GeneralTHDMcache; //forward reference to GeneralTHDMcache class
 
class GeneralTHDM: public NPbase {
public:
 
 
    static const int NGeneralTHDMvars = 66;
    static std::string GeneralTHDMvars[NGeneralTHDMvars];
 
    GeneralTHDM();
 
    ~GeneralTHDM();
    
    virtual bool InitializeModel();
    
    
    
    virtual bool Update(const std::map<std::string, double>& DPars);
    
    
    
    virtual bool PostUpdate();
    
    
    //Function to debug the gsl errors
    static void custom_gsl_error_handler(const char *reason, const char *file, int line, int gsl_errno) {
    printf("GSL Error: %s\nFile: %s\nLine: %d\nError code: %d\n", reason, file, line, gsl_errno);
    // Optionally print a stack trace or additional debug info
    }
    
    
    virtual bool CheckParameters(const std::map<std::string, double>& DPars);
    
    double getMuw() const
    {
        return muw;
    }
 
 
    virtual GeneralTHDMMatching& getMatching() const
    {
        return GTHDMM.getObj();
    }
 
 
    // Flags
 
    virtual bool setFlagStr(const std::string name, const std::string value);
    
    virtual bool setFlag(const std::string, const bool);
 
    GeneralTHDMcache* getMyGTHDMCache() const
    {
        return myGTHDMcache;
    }
 
 
    //double getlogtb() const {
    //    return logtb;
    //}
 
    //double gettanb() const {
    //    return tanb;
    //}
 
    //double getsinb() const {
    //    return sinb;
    //}
 
    //double getcosb() const {
    //    return cosb;
    //}
 
    double getmHp2() const {
        if(flag_use_sq_masses) {
            return mHp2;
        }
        else
        {
            return mHp1*mHp1;
        }
    }
 
    double getmHp() const {
        if(flag_use_sq_masses) {
            if(mHp2 < 0.) {
                throw std::runtime_error("error in GeneralTHDM: mHp2 < 0!");
            }
            else
            {
                return sqrt(mHp2);
            }
        }
        else
        {
                return mHp1;
        }
    }
 
 
    double getmH1sq() const {
        if(flag_use_sq_masses) {
             if(mH1sq < 0.) {
                throw std::runtime_error("error in GeneralTHDM: mH1sq < 0!");
            }
                return mH1sq;
            }
         else
        {
            return mH1*mH1;
        }
    }
    
    
    
    double getmH2sq() const {
        if(flag_use_sq_masses) {
             if(mH2sq < 0.) {
                throw std::runtime_error("error in GeneralTHDM: mH2sq < 0!");
            }
                return mH2sq;
            }
         else
        {
            return mH21*mH21;
        }
    }
 
      double getmH3sq() const {
        if(flag_use_sq_masses) {
             if(mH3sq < 0.) {
                throw std::runtime_error("error in GeneralTHDM: mH3sq < 0!");
            }
                return mH3sq;
            }
         else
        {
            return mH31*mH31;
        }
    }
 
    double getalpha1() const {
        return alpha1;
    }
 
    double getcosalpha1() const {
        return cosalpha1;
    }
 
    double getsinalpha1() const {
        return sinalpha1;
    }
    
    double gettanalpha1() const {
        return tanalpha1;
    }
    
    
    double getalpha2() const {
        if(flag_CPconservation==true) {
            return 0.0;
        }
        else {
            return alpha2;
        }
    }
 
    double getcosalpha2() const {
        if(flag_CPconservation==true) {
            return 1.0;
        }
        else {
            return cosalpha2;
        }
    }
 
    double getsinalpha2() const {
        if(flag_CPconservation==true) {
            return 0.0;
        }
        else {
            return sinalpha2;
        }
    }
 
    double getalpha3() const {
        if(flag_CPconservation==true) {
            return 0.0;
        }
        else {
            return alpha3;
        }
    }
 
    double getcosalpha3() const {
        if(flag_CPconservation==true) {
            return 1.0;
        }
        else {
            return cosalpha3;
        }
    }
 
    double getsinalpha3() const {
        if(flag_CPconservation==true) {
            return 0.0;
        }
        else {
            return sinalpha3;
        }
    }
    
    
    double getlambda1() const {
        return lambda1;
    }
    
    
    double getlambda2() const {
        return lambda2;
    }
    
    
    double getlambda3() const {
        return lambda3;
    }
    
    double getlambda4() const {
        return lambda4;
    }
    
    double getRelambda5() const {
        return Relambda5;
    }
    
    double getImlambda5() const {
        if(flag_CPconservation==true) {
            return 0.0;
        }
        else {
            return Imlambda5;
        }
    }
    
    double getRelambda6() const {
        return Relambda6;
    }
 
    
    double getImlambda6() const {
        if(flag_CPconservation==true) {
            return 0.0;
        }
        else {
            return Imlambda6;
        }
    }
    
    double getRelambda7() const {
        return Relambda7;
    }
    
    
    double getImlambda7() const {
        if(flag_CPconservation==true) {
            return 0.0;
        }
        else {
            return Imlambda7;
        }
    }
 
    gslpp::complex getNu_11() const {
        gslpp::complex i = gslpp::complex::i();
         if(flag_sigma) {
        return (Nu_11r + i*Nu_11i);
        }
         else
         {
             if(alpha1==0)
               return 0 +i*0;
             else
                return (yu1R_GTHDM - cosalpha1)/(sinalpha1);
         }
    }
       
    gslpp::complex getNu_12() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_12r + i*Nu_12i);
    }
    
    gslpp::complex getNu_13() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_13r + i*Nu_13i);
    }
    
    gslpp::complex getNu_21() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_21r + i*Nu_21i);
    }
    
    gslpp::complex getNu_22() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_22r + i*Nu_22i);
    }
      
    gslpp::complex getNu_23() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_23r + i*Nu_23i);
    }
    
    gslpp::complex getNu_31() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_31r + i*Nu_31i);
    }
    
    gslpp::complex getNu_32() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_32r + i*Nu_32i);
    }
    
    gslpp::complex getNu_33() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nu_33r + i*Nu_33i);
    }
    
    gslpp::complex getNd_11() const {
        gslpp::complex i = gslpp::complex::i();
        if(flag_sigma) {
        return (Nd_11r + i*Nd_11i);
        }
     else
         {
            if(alpha1==0)
                return 0 +i*0;
            else
             return (yd1R_GTHDM - cosalpha1)/(sinalpha1);
         }
    }
    
 
    gslpp::complex getNd_12() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_12r + i*Nd_12i);
    }
    
    gslpp::complex getNd_13() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_13r + i*Nd_13i);
    }
    
    gslpp::complex getNd_21() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_21r + i*Nd_21i);
    }
    
    gslpp::complex getNd_22() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_22r + i*Nd_22i);
    }
    
    gslpp::complex getNd_23() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_23r + i*Nd_23i);
    }
    
    gslpp::complex getNd_31() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_31r + i*Nd_31i);
    }
    
    gslpp::complex getNd_32() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_32r + i*Nd_32i);
    }
    
    gslpp::complex getNd_33() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nd_33r + i*Nd_33i);
    }
    
    gslpp::complex getNl_11() const {
        gslpp::complex i = gslpp::complex::i();
        if(flag_sigma) {
        return (Nl_11r + i*Nl_11i);
        }
        else
         {
            if(alpha1==0)
                return 0 +i*0;
            else
             return (yl1R_GTHDM  -cosalpha1)/(sinalpha1);
         }
    }
    
    gslpp::complex getNl_12() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_12r + i*Nl_12i);
    }
    
    gslpp::complex getNl_13() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_13r + i*Nl_13i);
    }
    
    gslpp::complex getNl_21() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_21r + i*Nl_21i);
    }
    
    gslpp::complex getNl_22() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_22r + i*Nl_22i);
    }
    
    gslpp::complex getNl_23() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_23r + i*Nl_23i);
    }
    
    gslpp::complex getNl_31() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_31r + i*Nl_31i);
    }
    
    gslpp::complex getNl_32() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_32r + i*Nl_32i);
    }
    
    gslpp::complex getNl_33() const {
        gslpp::complex i = gslpp::complex::i();
        return (Nl_33r + i*Nl_33i);
    }
    
    
    
    gslpp::complex getyu1() const {
        if(!flag_sigma) {
            return (yu1R_GTHDM);
        }
        else {
            gslpp::complex i = gslpp::complex::i();
            double R11 = cosalpha1*cosalpha2;
            double R12 = sinalpha1*cosalpha2;
            double R13 = -sinalpha2;
            gslpp::complex yu1 = (R11 + (R12 - i*R13)*(Nu_11r + i*Nu_11i).conjugate());
            // if(yu1.abs2() > 4.) return std::numeric_limits<double>::quiet_NaN();
            return yu1;
        }
    }
        
     
    double getyu1R() const {
        if(!flag_sigma) {
            return (yu1R_GTHDM);
        }
        else {
            double R11 = cosalpha1*cosalpha2;
            double R12 = sinalpha1*cosalpha2;
            double yu1R = (R11 + R12*Nu_11r);
            // if(yu1R*yu1R > 4.) return std::numeric_limits<double>::quiet_NaN();
            return yu1R;
        }
    }    
        
        
    gslpp::complex getyd1() const {
        if(!flag_sigma) {
            return (yd1R_GTHDM);
        }
        else {
            gslpp::complex i = gslpp::complex::i();
            double R11 = cosalpha1*cosalpha2;
            double R12 = sinalpha1*cosalpha2;
            double R13 = -sinalpha2;
            gslpp::complex yd1 = (R11 + (R12 - i*R13)*(Nd_11r + i*Nd_11i));
            // if(yd1.abs2() > 4.) return std::numeric_limits<double>::quiet_NaN();
            return yd1;
        }
    }
    
    
    double getyd1R() const {
        if(!flag_sigma) {
            return (yd1R_GTHDM);
        }
        else {
            double R11 = cosalpha1*cosalpha2;
            double R12 = sinalpha1*cosalpha2;
            double yd1R = (R11 + R12*Nd_11r);
            // if(yd1R*yd1R > 4.) return std::numeric_limits<double>::quiet_NaN();
            return yd1R;
        }
    }
    
    gslpp::complex getyl1() const {
        if(!flag_sigma) {
            return (yl1R_GTHDM);
        }
        else {
            gslpp::complex i = gslpp::complex::i();
            double R11 = cosalpha1*cosalpha2;
            double R12 = sinalpha1*cosalpha2;
            double R13 = -sinalpha2;
            gslpp::complex yl1 = (R11 + (R12 - i*R13)*(Nl_11r + i*Nl_11i)); 
            // if(yl1.abs2() > 4.) return std::numeric_limits<double>::quiet_NaN();
            return yl1;
        }
    }
    
    
    double getyl1R() const {
        if(!flag_sigma) {
            return (yl1R_GTHDM);
        }
         else {
            double R11 = cosalpha1*cosalpha2;
            double R12 = sinalpha1*cosalpha2;            
            double yl1R = (R11 + R12*Nl_11r);
            // if(yl1R*yl1R > 4.) return std::numeric_limits<double>::quiet_NaN();
            return yl1R;
        }
    }
    
    double getQ_GTHDM() const {
        return Q_GTHDM;
    }
 
    double getRpepsGTHDM() const {
        return RpepsGTHDM;
    }
 
    double getNLOuniscaleGTHDM() const {
        return NLOuniscaleGTHDM;
    }
 
    bool getATHDMflag() const {
        return flag_ATHDM;
    }
 
    bool getCPconservationflag() const {
        return flag_CPconservation;
    }
 
    std::string getRGEorderflag() const {
        return flag_RGEorder;
    }
    
    bool getsqmassesflag() const {
        return flag_use_sq_masses;
    }
 
    virtual const double muggH(const double sqrt_s) const;
    virtual const double muVBF(const double sqrt_s) const;
    virtual const double mueeWBF(const double sqrt_s, const double Pol_em, const double Pol_ep) const;
    virtual const double muWH(const double sqrt_s) const;
    virtual const double muZH(const double sqrt_s) const;
    virtual const double mueeZH(const double sqrt_s, const double Pol_em, const double Pol_ep) const;
    virtual const double muVH(const double sqrt_s) const;
    virtual const double muVBFpVH(const double sqrt_s) const;
    virtual const double muttH(const double sqrt_s) const;
    virtual const double mutH(const double sqrt_s) const;
    virtual const double mubbH(const double sqrt_s) const;
    virtual const double muttHptH(const double sqrt_s) const;
    virtual const double muggHpVBFpbbH(const double sqrt_s) const;
    virtual const double muggHpttHptHpbbH(const double sqrt_s) const;
    virtual const double muggHpbbH(const double sqrt_s) const;
    
    
    virtual const double GammaTotal() const;
    virtual const double BrHggRatio() const;
    virtual const double BrHWWRatio() const;
    virtual const double BrHZZRatio() const;
    virtual const double BrHZgaRatio() const;
    virtual const double BrHgagaRatio() const;
    virtual const double BrHmumuRatio() const;
    virtual const double BrHtautauRatio() const;
    virtual const double BrHccRatio() const;
    virtual const double BrHbbRatio() const;
    virtual const double muggHgaga(const double sqrt_s) const;
    virtual const double muggHpbbH_Hgaga(const double sqrt_s) const;
    virtual const double muVBFHgaga(const double sqrt_s) const;
    virtual const double muVHgaga(const double sqrt_s) const;
    virtual const double muWHgaga(const double sqrt_s) const;
    virtual const double muZHgaga(const double sqrt_s) const;
    virtual const double muttHgaga(const double sqrt_s) const;
    virtual const double mutHgaga(const double sqrt_s) const;
    virtual const double muttHptH_Hgaga(const double sqrt_s) const;
    
    
    
    virtual const double muggHZZ(const double sqrt_s) const;
    virtual const double muggHpbbH_HZZ(const double sqrt_s) const;
    virtual const double muVBFHZZ(const double sqrt_s) const;
    virtual const double muVHZZ(const double sqrt_s) const;
    virtual const double muZHZZ(const double sqrt_s) const;
    virtual const double muWHZZ(const double sqrt_s) const;
    virtual const double muttHZZ(const double sqrt_s) const;
    virtual const double muttHptH_HZZ(const double sqrt_s) const;
    virtual const double muggHWW(const double sqrt_s) const;
    virtual const double muggHpbbH_HWW(const double sqrt_s) const;
    virtual const double muVBFHWW(const double sqrt_s) const;
    virtual const double muWHWW(const double sqrt_s) const;
    virtual const double muZHWW(const double sqrt_s) const;
    virtual const double muVHWW(const double sqrt_s) const;
    virtual const double muttHWW(const double sqrt_s) const;
    virtual const double muttHptH_HWW(const double sqrt_s) const;
    virtual const double muggHtautau(const double sqrt_s) const;
    virtual const double muggHpbbH_Htautau(const double sqrt_s) const;
    virtual const double muVBFHtautau(const double sqrt_s) const;
    virtual const double muWHtautau(const double sqrt_s) const;
    virtual const double muZHtautau(const double sqrt_s) const;
    virtual const double muVHtautau(const double sqrt_s) const;
    virtual const double muttHtautau(const double sqrt_s) const;
    virtual const double muttHptH_Htautau(const double sqrt_s) const;
    virtual const double muggHbb(const double sqrt_s) const;
    virtual const double muggHpVBFpbbH_Hbb(const double sqrt_s) const;
    virtual const double muVBFHbb(const double sqrt_s) const;
    virtual const double muWHbb(const double sqrt_s) const;
    virtual const double muZHbb(const double sqrt_s) const;
    virtual const double muVHbb(const double sqrt_s) const;
    virtual const double muttHbb(const double sqrt_s) const;
    virtual const double muttHptH_Hbb(const double sqrt_s) const;
    virtual const double muVHcc(const double sqrt_s) const;
    virtual const double muppHmumu(const double sqrt_s) const;
    virtual const double muggHmumu(const double sqrt_s) const;
    virtual const double muggHpttHptHpbbH_Hmumu(const double sqrt_s) const;
    virtual const double muVBFHmumu(const double sqrt_s) const;
    virtual const double muVHmumu(const double sqrt_s) const;
    virtual const double muttHptH_Hmumu(const double sqrt_s) const;
    virtual const double muVBFpVH_Hmumu(const double sqrt_s) const;
    virtual const double muggHZga(const double sqrt_s) const;
    virtual const double muVBFHZga(const double sqrt_s) const;
    virtual const double muppHZga(const double sqrt_s) const;
    virtual const double computeGammaTotalRatio() const;
 
    
    const double F(const double m02, const double m12) const;
    
    
    
    const double GTHDMDeltaS() const;
    
    const double GTHDMDeltaT() const;
    
    
    const double GTHDMDeltaU() const;
    
    
    //  Why was this done??? It doesn't seem to make any sense
    /* 
     * @brief Method overriding the NPbase Mw with the tree-level W mass
     */
    //virtual const double Mw() const;
    
    //We substitute the tree-level value with the real SM value but the oblique corrections are missing
    /* 
     * @brief Method overriding the NPbase Mw with the SM W mass. We
     * did not include the oblique corrections for the moment but it 
     * should be done for consistency
     */
    virtual const double Mw() const;
 
protected:
 
    virtual void setParameter(const std::string, const double&);
 
    bool checkRGEorder(const std::string RGEorder) const
    {
        if (getModelName() == "GeneralTHDMZ2") {
            if (RGEorder.compare("LO") == 0
                || RGEorder.compare("approxNLO") == 0
                || RGEorder.compare("NLO") == 0)
                return true;
            else
                return false;
        }
        else {
            if (RGEorder.compare("LO") == 0)
//                || RGEorder.compare("approxNLO") == 0)
                return true;
            else
                return false;
        }
    }
 
    mutable Matching<GeneralTHDMMatching,GeneralTHDM> GTHDMM; 
 
private:
 
 
    
    GeneralTHDMcache* myGTHDMcache;
 
    //Let's define here all the parameters of the model, including those who are linearly dependent
    //We'll also set all the parameters in the setParameter functions in such a way that we can use the formulae written in terms of any parameter
    //These parameters are private but used in the inherited classes, think a bit about if this is the best way of including them
    double /*logtb, tanb, sinb, cosb,*/mH1sq, mHp2, mH2sq, mH3sq, mH1, mHp1, mH21, mH31, alpha1, cosalpha1, sinalpha1, tanalpha1, alpha2, cosalpha2, sinalpha2,
            alpha3, cosalpha3, sinalpha3, mu2, lambda1, lambda2, lambda3, lambda4, Relambda5, Imlambda5, Relambda6, Imlambda6, Relambda7, Imlambda7,//we can trade one of the Im by the Im of lambda7, we include it any way since we've not decided yet which phase to keep
            yu1R_GTHDM, yd1R_GTHDM, yl1R_GTHDM, 
            Nu_11r, Nu_11i, Nu_12r, Nu_12i, Nu_13r, Nu_13i, 
            Nu_21r, Nu_21i, Nu_22r, Nu_22i, Nu_23r, Nu_23i, 
            Nu_31r, Nu_31i, Nu_32r, Nu_32i, Nu_33r, Nu_33i, 
            Nd_11r, Nd_11i, Nd_12r, Nd_12i, Nd_13r, Nd_13i, 
            Nd_21r, Nd_21i, Nd_22r, Nd_22i, Nd_23r, Nd_23i, 
            Nd_31r, Nd_31i, Nd_32r, Nd_32i, Nd_33r, Nd_33i, 
            Nl_11r, Nl_11i, Nl_12r, Nl_12i, Nl_13r, Nl_13i, 
            Nl_21r, Nl_21i, Nl_22r, Nl_22i, Nl_23r, Nl_23i, 
            Nl_31r, Nl_31i, Nl_32r, Nl_32i, Nl_33r, Nl_33i, 
            Q_GTHDM, RpepsGTHDM, NLOuniscaleGTHDM;
    bool flag_ATHDM, flag_CPconservation, flag_use_sq_masses, flag_sigma, flag_SM_Higgs;
    std::string flag_RGEorder;
};
 
#endif  /* GENERALTHDM_H */