Add some initial unit tests for the VCMSSM

test/module.mk

@@ -292,6 +292,13 @@ TEST_SRC += \
 		$(DIR)/test_NUTNMSSM_spectrum.cpp
 endif
 
+ifeq ($(WITH_VCMSSM) $(WITH_CMSSM),yes yes)
+TEST_SRC += \
+		$(DIR)/test_VCMSSM_ewsb.cpp
+TEST_SH += \
+		$(DIR)/test_VCMSSM_spectrum.sh
+endif
+
 ifeq ($(WITH_HSSUSY),yes)
 TEST_SH += \
 		$(DIR)/test_HSSUSY_SUSYHD.sh
@@ -734,6 +741,8 @@ $(DIR)/test_SMHighPrecision_two_loop_spectrum.x: $(LIBSMHighPrecision) $(LIBFLEX
 
 $(DIR)/test_NSM_low_scale_constraint.x: $(LIBNSM) $(LIBFLEXI) $(LIBLEGACY) $(filter-out -%,$(LOOPFUNCLIBS))
 
+$(DIR)/test_VCMSSM_ewsb.x: $(LIBVCMSSM) $(LIBCMSSM) $(LIBFLEXI) $(LIBLEGACY) $(filter-out -%,$(LOOPFUNCLIBS))
+
 # general test rule which links all libraries needed for a generated model
 $(DIR)/test_%.x: $(DIR)/test_%.o
 		$(CXX) -o $@ $(call abspathx,$^) $(filter -%,$(LOOPFUNCLIBS)) $(BOOSTTESTLIBS) $(BOOSTTHREADLIBS) $(THREADLIBS) $(GSLLIBS) $(LAPACKLIBS) $(BLASLIBS) $(FLIBS)

test/test_VCMSSM.hpp

@@ -0,0 +1,120 @@
+#ifndef TEST_VCMSSM_H
+#define TEST_VCMSSM_H
+
+#include "VCMSSM_two_scale_model.hpp"
+#include "CMSSM_two_scale_model.hpp"
+
+#include "ew_input.hpp"
+#include "wrappers.hpp"
+
+void setup_VCMSSM_const(flexiblesusy::VCMSSM<flexiblesusy::Two_scale>& m,
+                         const flexiblesusy::VCMSSM_input_parameters& input)
+{
+   using namespace flexiblesusy;
+
+   const double ALPHASMZ = 0.1176;
+   const double ALPHAMZ = 1.0 / 127.918;
+   const double sinthWsq = 0.23122;
+   const double alpha1 = 5.0 * ALPHAMZ / (3.0 * (1.0 - sinthWsq));
+   const double alpha2 = ALPHAMZ / sinthWsq;
+   const double g1 = Sqrt(4.0 * Pi * alpha1);
+   const double g2 = Sqrt(4.0 * Pi * alpha2);
+   const double g3 = Sqrt(4.0 * Pi * ALPHASMZ);
+   const double M12 = input.m12;
+   const double m0 = input.m0;
+   const double a0 = input.Azero + 100.;
+   const double root2 = Sqrt(2.0);
+   const double vev = 246.0;
+   const double vu = vev * Sin(ArcTan(input.TBGuess));
+   const double vd = vev * Cos(ArcTan(input.TBGuess));
+   const double scale = Electroweak_constants::MZ;
+   const double mu_guess = Electroweak_constants::MZ;
+
+   Eigen::Matrix<double,3,3> mm0(Eigen::Matrix<double,3,3>::Zero()),
+      mm0b(Eigen::Matrix<double,3,3>::Zero()),
+      Yu(Eigen::Matrix<double,3,3>::Zero()),
+      Yd(Eigen::Matrix<double,3,3>::Zero()),
+      Ye(Eigen::Matrix<double,3,3>::Zero());
+
+   mm0(0,0) = Sqr(m0);
+   mm0(1,1) = Sqr(m0 + 1);
+   mm0(2,2) = Sqr(m0 + 2);
+   mm0b(0,0) = Sqr(m0 + 3);
+   mm0b(1,1) = Sqr(m0 + 4);
+   mm0b(2,2) = Sqr(m0 + 5);
+
+   Yu(2,2) = 165.0   * root2 / vu;
+   Yd(2,2) = 2.9     * root2 / vd;
+   Ye(2,2) = 1.77699 * root2 / vd;
+
+   m.set_scale(scale);
+   m.set_loops(1);
+   m.set_g1(g1);
+   m.set_g2(g2);
+   m.set_g3(g3);
+   m.set_Yu(Yu);
+   m.set_Yd(Yd);
+   m.set_Ye(Ye);
+   m.set_Mu(mu_guess);
+   m.set_MassB(M12);
+   m.set_MassWB(M12);
+   m.set_MassG(M12);
+   m.set_mq2(mm0);
+   m.set_ml2(mm0);
+   m.set_md2(mm0b);
+   m.set_mu2(mm0);
+   m.set_me2(mm0b);
+   m.set_mHd2(Sqr(m0));
+   m.set_mHu2(Sqr(m0));
+   m.set_TYu(a0 * Yu);
+   m.set_TYd(a0 * Yd);
+   m.set_TYe(a0 * Ye);
+   m.set_BMu((a0 + m0) * m.get_Mu());
+   m.set_vu(vu);
+   m.set_vd(vd);
+   m.set_vev(vev);
+}
+
+void setup_VCMSSM(flexiblesusy::VCMSSM<flexiblesusy::Two_scale>& m,
+                  flexiblesusy::VCMSSM_input_parameters& input)
+{
+   input.m0 = 125.;
+   input.m12 = 500.;
+   input.Azero = 0.;
+   input.SignMu = 1.;
+   input.TBGuess = 10.0;
+
+   setup_VCMSSM_const(m, input);
+}
+
+void match_CMSSM_to_VCMSSM(flexiblesusy::CMSSM<flexiblesusy::Two_scale>& cmssm,
+                           const flexiblesusy::VCMSSM<flexiblesusy::Two_scale>& vcmssm)
+{
+   cmssm.set_scale(vcmssm.get_scale());
+   cmssm.set_loops(vcmssm.get_loops());
+   cmssm.set_g1(vcmssm.get_g1());
+   cmssm.set_g2(vcmssm.get_g2());
+   cmssm.set_g3(vcmssm.get_g3());
+   cmssm.set_Yu(vcmssm.get_Yu());
+   cmssm.set_Yd(vcmssm.get_Yd());
+   cmssm.set_Ye(vcmssm.get_Ye());
+   cmssm.set_Mu(vcmssm.get_Mu());
+   cmssm.set_MassB(vcmssm.get_MassB());
+   cmssm.set_MassWB(vcmssm.get_MassWB());
+   cmssm.set_MassG(vcmssm.get_MassG());
+   cmssm.set_mq2(vcmssm.get_mq2());
+   cmssm.set_ml2(vcmssm.get_ml2());
+   cmssm.set_mu2(vcmssm.get_mu2());
+   cmssm.set_md2(vcmssm.get_md2());
+   cmssm.set_me2(vcmssm.get_me2());
+   cmssm.set_mHd2(vcmssm.get_mHd2());
+   cmssm.set_mHu2(vcmssm.get_mHu2());
+   cmssm.set_TYu(vcmssm.get_TYu());
+   cmssm.set_TYd(vcmssm.get_TYd());
+   cmssm.set_TYe(vcmssm.get_TYe());
+   cmssm.set_BMu(vcmssm.get_BMu());
+   cmssm.set_vu(vcmssm.get_vu());
+   cmssm.set_vd(vcmssm.get_vd());
+}
+
+#endif

test/test_VCMSSM_ewsb.cpp

@@ -0,0 +1,124 @@
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE test_VCMSSM_ewsb
+
+#include <boost/test/unit_test.hpp>
+
+#include "test_VCMSSM.hpp"
+
+using namespace flexiblesusy;
+
+BOOST_AUTO_TEST_CASE( test_VCMSSM_ewsb_tree_level )
+{
+   VCMSSM_input_parameters input;
+   VCMSSM<Two_scale> vcmssm(input);
+   const double precision = 1.0e-5;
+   setup_VCMSSM(vcmssm, input);
+
+   // initial guess
+   const double vev = Sqrt(Sqr(vcmssm.get_vd()) + Sqr(vcmssm.get_vu()));
+   vcmssm.set_Mu(100.);
+   vcmssm.set_vu(vev * Sin(ArcTan(input.TBGuess)));
+   vcmssm.set_vd(vev * Cos(ArcTan(input.TBGuess)));
+
+   vcmssm.set_ewsb_iteration_precision(precision);
+   const int vcmssm_error = vcmssm.solve_ewsb_tree_level();
+
+   BOOST_CHECK_EQUAL(vcmssm_error, 0);
+
+   BOOST_CHECK_SMALL(vcmssm.get_ewsb_eq_hh_1(), precision);
+   BOOST_CHECK_SMALL(vcmssm.get_ewsb_eq_hh_2(), precision);
+
+   const double vcmssm_Mu_soln = vcmssm.get_Mu();
+   const double vcmssm_BMu = vcmssm.get_BMu();
+
+   CMSSM<Two_scale> cmssm;
+   match_CMSSM_to_VCMSSM(cmssm, vcmssm);
+
+   // initial guess: VCMSSM solution with small perturbation
+   const double shift = 5.;
+   cmssm.set_Mu(vcmssm_Mu_soln + shift);
+   cmssm.set_BMu(vcmssm_BMu + shift);
+
+   cmssm.set_ewsb_iteration_precision(precision);
+   const int cmssm_error = cmssm.solve_ewsb_tree_level();
+
+   BOOST_CHECK_EQUAL(cmssm_error, 0);
+
+   BOOST_CHECK_SMALL(cmssm.get_ewsb_eq_hh_1(), precision);
+   BOOST_CHECK_SMALL(cmssm.get_ewsb_eq_hh_2(), precision);
+
+   const double cmssm_Mu_soln = cmssm.get_Mu();
+   const double cmssm_BMu_soln = cmssm.get_BMu();
+
+   BOOST_CHECK_CLOSE_FRACTION(vcmssm_Mu_soln, cmssm_Mu_soln, precision);
+   BOOST_CHECK_CLOSE_FRACTION(vcmssm_BMu, cmssm_BMu_soln, precision);
+}
+
+BOOST_AUTO_TEST_CASE( test_VCMSSM_ewsb_one_loop )
+{
+   VCMSSM_input_parameters vcmssm_input;
+   VCMSSM<Two_scale> vcmssm(vcmssm_input);
+   const double precision = 1.0e-5;
+   setup_VCMSSM(vcmssm, vcmssm_input);
+
+   // initial guess
+   const double vev = Sqrt(Sqr(vcmssm.get_vd()) + Sqr(vcmssm.get_vu()));
+   vcmssm.set_Mu(100.);
+   vcmssm.set_vu(vev * Sin(ArcTan(vcmssm_input.TBGuess)));
+   vcmssm.set_vd(vev * Cos(ArcTan(vcmssm_input.TBGuess)));
+
+   vcmssm.calculate_DRbar_masses();
+
+   vcmssm.set_ewsb_iteration_precision(precision);
+   const int vcmssm_error = vcmssm.solve_ewsb_one_loop();
+
+   BOOST_CHECK_EQUAL(vcmssm_error, 0);
+
+   const std::complex<double> vcmssm_tadpole_hh_1(vcmssm.tadpole_hh(0));
+   const std::complex<double> vcmssm_tadpole_hh_2(vcmssm.tadpole_hh(1));
+
+   BOOST_CHECK_SMALL(Im(vcmssm_tadpole_hh_1), 1.0e-12);
+   BOOST_CHECK_SMALL(Im(vcmssm_tadpole_hh_2), 1.0e-12);
+
+   BOOST_CHECK_SMALL(vcmssm.get_ewsb_eq_hh_1() - Re(vcmssm_tadpole_hh_1), 1.0);
+   BOOST_CHECK_SMALL(vcmssm.get_ewsb_eq_hh_2() - Re(vcmssm_tadpole_hh_2), 1.0);
+
+   const double vcmssm_Mu_soln = vcmssm.get_Mu();
+   const double vcmssm_BMu = vcmssm.get_BMu();
+
+   CMSSM_input_parameters cmssm_input;
+   cmssm_input.m12 = vcmssm_input.m12;
+   cmssm_input.m0 = vcmssm_input.m0;
+   cmssm_input.Azero = vcmssm_input.Azero;
+   cmssm_input.TanBeta = vcmssm.get_vu() / vcmssm.get_vd();
+   cmssm_input.SignMu = vcmssm_input.SignMu;
+   CMSSM<Two_scale> cmssm(cmssm_input);
+   match_CMSSM_to_VCMSSM(cmssm, vcmssm);
+
+   cmssm.calculate_DRbar_masses();
+
+   // initial guess: VCMSSM solution with small perturbation
+   const double shift = 5.;
+   cmssm.set_Mu(vcmssm_Mu_soln + shift);
+   cmssm.set_BMu(vcmssm_BMu + shift);
+
+   cmssm.set_ewsb_iteration_precision(precision);
+   const int cmssm_error = cmssm.solve_ewsb_one_loop();
+
+   BOOST_CHECK_EQUAL(cmssm_error, 0);
+
+   const std::complex<double> cmssm_tadpole_hh_1(cmssm.tadpole_hh(0));
+   const std::complex<double> cmssm_tadpole_hh_2(cmssm.tadpole_hh(1));
+
+   BOOST_CHECK_SMALL(Im(cmssm_tadpole_hh_1), 1.0e-12);
+   BOOST_CHECK_SMALL(Im(cmssm_tadpole_hh_2), 1.0e-12);
+
+   BOOST_CHECK_SMALL(cmssm.get_ewsb_eq_hh_1() - Re(cmssm_tadpole_hh_1), 1.0);
+   BOOST_CHECK_SMALL(cmssm.get_ewsb_eq_hh_2() - Re(cmssm_tadpole_hh_2), 1.0);
+
+   const double cmssm_Mu_soln = cmssm.get_Mu();
+   const double cmssm_BMu_soln = cmssm.get_BMu();
+
+   BOOST_CHECK_CLOSE_FRACTION(vcmssm_Mu_soln, cmssm_Mu_soln, precision);
+   BOOST_CHECK_CLOSE_FRACTION(vcmssm_BMu, cmssm_BMu_soln, precision);
+}