put calculation of tree-level Higgs mass matrix into separate function

templates/standard_model_matching.cpp.in

@@ -168,6 +168,21 @@ double calculate_Mh2_pole(const Standard_model& sm, const Standard_model& sm_0l)
    return (tadpole_lo - tadpole_lom1).eval();
 }
 
+/**
+ * Calculates tree-level Higgs mass matrix in the @ModelName@ without
+ * including tadpoles implicitly.  This is achieved by solving the
+ * EWSB equations at tree-level in order to avoid the inclusion of
+ * loop tadpoles.
+ *
+ * @param model @ModelName@ parameters
+ * @return tree-level Higgs mass matrix in the @ModelName@ w/o tadpoles
+ */
+auto calculate_mh2_tree_level(@ModelName@_mass_eigenstates model) -> decltype(model.get_mass_matrix_@HiggsBoson@())
+{
+   model.solve_ewsb_tree_level();
+   return model.get_mass_matrix_@HiggsBoson@();
+}
+
 /**
  * Calculates squared Higgs pole mass in the @ModelName@,
  * \f$(M_h^{\text{@ModelName@}})^2\f$.
@@ -179,14 +194,12 @@ double calculate_Mh2_pole(const Standard_model& sm, const Standard_model& sm_0l)
  * @return squared Higgs pole mass in the @ModelName@
  */
 double calculate_Mh2_pole(
-   @ModelName@_mass_eigenstates model,
+   const @ModelName@_mass_eigenstates& model,
    const @ModelName@_mass_eigenstates& model_0l,
    unsigned idx)
 {
-   // solve at tree-level to obtain plain tree-level mass matrix
-   model.solve_ewsb_tree_level();
-   model.calculate_DRbar_masses();
-   const auto mh2_tree = model.get_mass_matrix_@HiggsBoson@();
+   // calculate tree-level mass matrix
+   const auto mh2_tree = calculate_mh2_tree_level(model);
 
    // calculate 1L self-energy using tree-level parameters
    const double p = model_0l.get_M@HiggsBoson(idx)@;