implement calculation of SM fermion pole masses

templates/standard_model_matching.cpp.in

@@ -246,28 +246,88 @@ Eigen::Matrix<double,3,3> calculate_MFe_DRbar_tree_level(const @ModelName@_mass_
    return mf;
 }
 
-Eigen::Matrix<double,3,3> calculate_MFu_pole_1loop(
-   const Standard_model& sm_0l,
-   const Standard_model& sm_1l)
+double calculate_MFu_pole_1loop(
+   unsigned i,
+   const Standard_model& sm_0l)
 {
-   // @todo implement me
-   return ZEROMATRIX(3,3);
+   const double p = sm_0l.get_MFu(i);
+   const auto self_energy_1  = Re(sm_0l.self_energy_Fu_1(p));
+   const auto self_energy_PL = Re(sm_0l.self_energy_Fu_PL(p));
+   const auto self_energy_PR = Re(sm_0l.self_energy_Fu_PR(p));
+   const auto M_tree = sm_0l.get_mass_matrix_Fu();
+   const auto M_loop = (M_tree - self_energy_PR * M_tree
+                        - M_tree * self_energy_PL - self_energy_1).eval();
+
+   Eigen::Array<double,3,1> MFu_pole;
+   fs_svd(M_loop, MFu_pole);
+
+   return MFu_pole(i);
 }
 
-Eigen::Matrix<double,3,3> calculate_MFd_pole_1loop(
-   const Standard_model& sm_0l,
-   const Standard_model& sm_1l)
+double calculate_MFd_pole_1loop(
+   unsigned i,
+   const Standard_model& sm_0l)
 {
-   // @todo implement me
-   return ZEROMATRIX(3,3);
+   const double p = sm_0l.get_MFd(i);
+   const auto self_energy_1  = Re(sm_0l.self_energy_Fd_1(p));
+   const auto self_energy_PL = Re(sm_0l.self_energy_Fd_PL(p));
+   const auto self_energy_PR = Re(sm_0l.self_energy_Fd_PR(p));
+   const auto M_tree = sm_0l.get_mass_matrix_Fd();
+   const auto M_loop = (M_tree - self_energy_PR * M_tree
+                        - M_tree * self_energy_PL - self_energy_1).eval();
+
+   Eigen::Array<double,3,1> MFd_pole;
+   fs_svd(M_loop, MFd_pole);
+
+   return MFd_pole(i);
 }
 
-Eigen::Matrix<double,3,3> calculate_MFe_pole_1loop(
-   const Standard_model& sm_0l,
-   const Standard_model& sm_1l)
+double calculate_MFe_pole_1loop(
+   unsigned i,
+   const Standard_model& sm_0l)
 {
-   // @todo implement me
-   return ZEROMATRIX(3,3);
+   const double p = sm_0l.get_MFe(i);
+   const auto self_energy_1  = Re(sm_0l.self_energy_Fe_1(p));
+   const auto self_energy_PL = Re(sm_0l.self_energy_Fe_PL(p));
+   const auto self_energy_PR = Re(sm_0l.self_energy_Fe_PR(p));
+   const auto M_tree = sm_0l.get_mass_matrix_Fe();
+   const auto M_loop = (M_tree - self_energy_PR * M_tree
+                        - M_tree * self_energy_PL - self_energy_1).eval();
+
+   Eigen::Array<double,3,1> MFe_pole;
+   fs_svd(M_loop, MFe_pole);
+
+   return MFe_pole(i);
+}
+
+Eigen::Matrix<double,3,3> calculate_MFu_pole_1loop(const Standard_model& sm_0l)
+{
+   Eigen::Matrix<double,3,1> MFu_pole;
+   MFu_pole << calculate_MFu_pole_1loop(0, sm_0l),
+               calculate_MFu_pole_1loop(1, sm_0l),
+               calculate_MFu_pole_1loop(2, sm_0l);
+
+   return MFu_pole.asDiagonal();
+}
+
+Eigen::Matrix<double,3,3> calculate_MFd_pole_1loop(const Standard_model& sm_0l)
+{
+   Eigen::Matrix<double,3,1> MFd_pole;
+   MFd_pole << calculate_MFd_pole_1loop(0, sm_0l),
+               calculate_MFd_pole_1loop(1, sm_0l),
+               calculate_MFd_pole_1loop(2, sm_0l);
+
+   return MFd_pole.asDiagonal();
+}
+
+Eigen::Matrix<double,3,3> calculate_MFe_pole_1loop(const Standard_model& sm_0l)
+{
+   Eigen::Matrix<double,3,1> MFe_pole;
+   MFe_pole << calculate_MFe_pole_1loop(0, sm_0l),
+               calculate_MFe_pole_1loop(1, sm_0l),
+               calculate_MFe_pole_1loop(2, sm_0l);
+
+   return MFe_pole.asDiagonal();
 }
 
 Eigen::Matrix<double,3,3> calculate_MFu_pole_1loop(
@@ -296,13 +356,12 @@ Eigen::Matrix<double,3,3> calculate_MFe_pole_1loop(
 
 Eigen::Matrix<double,3,3> calculate_MFu_DRbar_1loop(
    const Standard_model& sm_0l,
-   const Standard_model& sm_1l,
    const @ModelName@_mass_eigenstates& model_0l,
    const @ModelName@_mass_eigenstates& model_1l)
 {
    Eigen::Matrix<double,3,3> mf_sm = ZEROMATRIX(3,3);
 
-   const auto Mf_sm  = calculate_MFu_pole_1loop(sm_0l, sm_1l);
+   const auto Mf_sm  = calculate_MFu_pole_1loop(sm_0l);
    const auto Mf_bsm = calculate_MFu_pole_1loop(model_0l, model_1l);
    const auto mf_bsm = calculate_MFu_DRbar_tree_level(model_1l);
 
@@ -313,13 +372,12 @@ Eigen::Matrix<double,3,3> calculate_MFu_DRbar_1loop(
 
 Eigen::Matrix<double,3,3> calculate_MFd_DRbar_1loop(
    const Standard_model& sm_0l,
-   const Standard_model& sm_1l,
    const @ModelName@_mass_eigenstates& model_0l,
    const @ModelName@_mass_eigenstates& model_1l)
 {
    Eigen::Matrix<double,3,3> mf_sm = ZEROMATRIX(3,3);
 
-   const auto Mf_sm  = calculate_MFd_pole_1loop(sm_0l, sm_1l);
+   const auto Mf_sm  = calculate_MFd_pole_1loop(sm_0l);
    const auto Mf_bsm = calculate_MFd_pole_1loop(model_0l, model_1l);
    const auto mf_bsm = calculate_MFd_DRbar_tree_level(model_1l);
 
@@ -330,13 +388,12 @@ Eigen::Matrix<double,3,3> calculate_MFd_DRbar_1loop(
 
 Eigen::Matrix<double,3,3> calculate_MFe_DRbar_1loop(
    const Standard_model& sm_0l,
-   const Standard_model& sm_1l,
    const @ModelName@_mass_eigenstates& model_0l,
    const @ModelName@_mass_eigenstates& model_1l)
 {
    Eigen::Matrix<double,3,3> mf_sm = ZEROMATRIX(3,3);
 
-   const auto Mf_sm  = calculate_MFe_pole_1loop(sm_0l, sm_1l);
+   const auto Mf_sm  = calculate_MFe_pole_1loop(sm_0l);
    const auto Mf_bsm = calculate_MFe_pole_1loop(model_0l, model_1l);
    const auto mf_bsm = calculate_MFe_DRbar_tree_level(model_1l);
 
@@ -448,9 +505,9 @@ double calculate_MZ_DRbar_1loop(
    const double mZ2_1L = Sqr(calculate_MZ_DRbar_1loop(sm_0l, sm, model_0l, model));
 
    // running @ModelName@ quark and lepton masses (via 1L matching)
-   const Eigen::Matrix<double, 3, 3> upQuarksDRbar    = calculate_MFu_DRbar_1loop(sm_0l, sm, model_0l, model);
-   const Eigen::Matrix<double, 3, 3> downQuarksDRbar  = calculate_MFd_DRbar_1loop(sm_0l, sm, model_0l, model);
-   const Eigen::Matrix<double, 3, 3> downLeptonsDRbar = calculate_MFe_DRbar_1loop(sm_0l, sm, model_0l, model);
+   const Eigen::Matrix<double, 3, 3> upQuarksDRbar    = calculate_MFu_DRbar_1loop(sm_0l, model_0l, model);
+   const Eigen::Matrix<double, 3, 3> downQuarksDRbar  = calculate_MFd_DRbar_1loop(sm_0l, model_0l, model);
+   const Eigen::Matrix<double, 3, 3> downLeptonsDRbar = calculate_MFe_DRbar_1loop(sm_0l, model_0l, model);
 
    // running @ModelName@ gauge couplings (via 1L matching)
    const double g1_1L = AbsSqrt(4. * Pi * alpha_em * (1. + delta_alpha_em) * mZ2_1L / mW2_1L) / @ModelName@_info::normalization_@hyperchargeCoupling@;