From cabdd5401ecd0b50098b6ea164c6282c9c1845ea Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 2 Sep 2016 19:51:01 +0200 Subject: [PATCH 1/5] use more non-universal NMSSM parameter to avoid Nan in the tadpoles --- test/test_NMSSM.hpp | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/test/test_NMSSM.hpp b/test/test_NMSSM.hpp index 7f220a0aa..852aedbc6 100644 --- a/test/test_NMSSM.hpp +++ b/test/test_NMSSM.hpp @@ -5,6 +5,7 @@ #include "wrappers.hpp" #include "ew_input.hpp" #include "nmssmsoftsusy.h" +#include "conversion.hpp" #include "NMSSM_two_scale_model.hpp" using namespace flexiblesusy; @@ -27,7 +28,7 @@ void setup_NMSSM_const(NMSSM& m, NmssmSoftsusy& s, const NMSSM_input_ const double cosBeta = cos(atan(tanBeta)); const double M12 = input.m12; const double m0 = input.m0; - const double a0 = input.Azero; + const double a0 = input.Azero + 100; const double root2 = sqrt(2.0); const double vev = 246.0; const double vu = vev * sinBeta; @@ -41,18 +42,23 @@ void setup_NMSSM_const(NMSSM& m, NmssmSoftsusy& s, const NMSSM_input_ Yu_SS(3,3) = 165.0 * root2 / (vev * sinBeta); Yd_SS(3,3) = 2.9 * root2 / (vev * cosBeta); Ye_SS(3,3) = 1.77699 * root2 / (vev * cosBeta); - DoubleMatrix ID(3, 3), mm0_SS(3, 3); - for (int i=1; i<=3; i++) ID(i, i) = 1.0; - mm0_SS = ID * sqr(m0); + DoubleMatrix mm0_SS(3, 3), mm0_SSb(3, 3); + mm0_SS(1,1) = sqr(m0 + 0); + mm0_SS(2,2) = sqr(m0 + 1); + mm0_SS(3,3) = sqr(m0 + 2); + mm0_SSb(1,1) = sqr(m0 + 3); + mm0_SSb(2,2) = sqr(m0 + 4); + mm0_SSb(3,3) = sqr(m0 + 5); Eigen::Matrix Yu(Eigen::Matrix::Zero()), Yd(Eigen::Matrix::Zero()), Ye(Eigen::Matrix::Zero()), - mm0(Eigen::Matrix::Zero()); + mm0, mm0b; Yu(2,2) = 165.0 * root2 / (vev * sinBeta); Yd(2,2) = 2.9 * root2 / (vev * cosBeta); Ye(2,2) = 1.77699 * root2 / (vev * cosBeta); - mm0 = Sqr(m0) * Eigen::Matrix::Identity(); + mm0 = ToEigenMatrix(mm0_SS); + mm0b = ToEigenMatrix(mm0_SSb); m.set_scale(scale); m.set_loops(1); @@ -69,9 +75,9 @@ void setup_NMSSM_const(NMSSM& m, NmssmSoftsusy& s, const NMSSM_input_ m.set_MassWB(M12); m.set_mq2(mm0); m.set_ml2(mm0); - m.set_md2(mm0); + m.set_md2(mm0b); m.set_mu2(mm0); - m.set_me2(mm0); + m.set_me2(mm0b); m.set_mHd2(Sqr(m0)); m.set_mHu2(Sqr(m0)); m.set_ms2(Sqr(m0)); @@ -101,9 +107,9 @@ void setup_NMSSM_const(NMSSM& m, NmssmSoftsusy& s, const NMSSM_input_ s.setGauginoMass(3, M12); s.setSoftMassMatrix(mQl, mm0_SS); s.setSoftMassMatrix(mUr, mm0_SS); - s.setSoftMassMatrix(mDr, mm0_SS); + s.setSoftMassMatrix(mDr, mm0_SSb); s.setSoftMassMatrix(mLl, mm0_SS); - s.setSoftMassMatrix(mEr, mm0_SS); + s.setSoftMassMatrix(mEr, mm0_SSb); s.setMh1Squared(sqr(m0)); s.setMh2Squared(sqr(m0)); s.setMsSquared(sqr(m0)); From 46131aa31aad098e0dd609b0d7efa0631c24518b Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 2 Sep 2016 19:51:47 +0200 Subject: [PATCH 2/5] relax test precision for new parameter point --- test/test_NMSSM_ewsb.cpp | 4 ++-- test/test_NMSSM_self_energies.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_NMSSM_ewsb.cpp b/test/test_NMSSM_ewsb.cpp index 192f9dc36..238941d57 100644 --- a/test/test_NMSSM_ewsb.cpp +++ b/test/test_NMSSM_ewsb.cpp @@ -165,7 +165,7 @@ BOOST_AUTO_TEST_CASE( test_NMSSM_one_loop_ewsb ) const double vS_fs = m.get_vS(); const double ms2_fs = m.get_ms2(); - BOOST_CHECK_CLOSE_FRACTION(kappa_ss, kappa_fs, 2.0e-9); + BOOST_CHECK_CLOSE_FRACTION(kappa_ss, kappa_fs, 2.1e-9); BOOST_CHECK_CLOSE_FRACTION(vS_ss , vS_fs , 3.0e-8); BOOST_CHECK_CLOSE_FRACTION(ms2_ss , ms2_fs , 5.0e-8); } @@ -282,7 +282,7 @@ BOOST_AUTO_TEST_CASE( test_NMSSM_two_loop_ewsb ) const double vS_fs = m.get_vS(); const double ms2_fs = m.get_ms2(); - BOOST_CHECK_CLOSE_FRACTION(kappa_ss, kappa_fs, 5.0e-9); + BOOST_CHECK_CLOSE_FRACTION(kappa_ss, kappa_fs, 5.5e-9); BOOST_CHECK_CLOSE_FRACTION(vS_ss , vS_fs , 4.0e-8); BOOST_CHECK_CLOSE_FRACTION(ms2_ss , ms2_fs , 8.0e-8); } diff --git a/test/test_NMSSM_self_energies.cpp b/test/test_NMSSM_self_energies.cpp index ff71d2e1d..57f18003e 100644 --- a/test/test_NMSSM_self_energies.cpp +++ b/test/test_NMSSM_self_energies.cpp @@ -74,7 +74,7 @@ BOOST_AUTO_TEST_CASE( test_NMSSM_self_energy_neutral_higgs ) const double vS_fs = m.get_vS(); const double ms2_fs = m.get_ms2(); - BOOST_CHECK_CLOSE_FRACTION(kappa_ss, kappa_fs, 5.0e-9); + BOOST_CHECK_CLOSE_FRACTION(kappa_ss, kappa_fs, 5.5e-9); BOOST_CHECK_CLOSE_FRACTION(vS_ss , vS_fs , 3.6e-8); BOOST_CHECK_CLOSE_FRACTION(ms2_ss , ms2_fs , 7.5e-8); From 81ff9f6bf68bf0fafcbcf99cc6a2e2943554bd78 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 2 Sep 2016 20:04:57 +0200 Subject: [PATCH 3/5] relax test precision after parameter point change Note: The tested point still leads to NaNs in the O(at*at) and O(atau*atau) corrections, which leads to different 2L tadpoles in FS and SS. --- test/test_NMSSM_ewsb.cpp | 11 ++++++++--- test/test_NMSSM_low_scale_constraint.cpp | 18 +++++++++--------- 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/test/test_NMSSM_ewsb.cpp b/test/test_NMSSM_ewsb.cpp index 238941d57..5ef5394d5 100644 --- a/test/test_NMSSM_ewsb.cpp +++ b/test/test_NMSSM_ewsb.cpp @@ -215,9 +215,14 @@ BOOST_AUTO_TEST_CASE( test_NMSSM_two_loop_tadpoles ) const double tu_1_and_2loop_ss = s.displayTadpole2Ms(); const double ts_1_and_2loop_ss = s.displayTadpoleSMs(); - BOOST_CHECK_CLOSE(two_loop_tadpole[0] / vd, td_1_and_2loop_ss - tadpole_ss_1, 1.0e-10); - BOOST_CHECK_CLOSE(two_loop_tadpole[1] / vu, tu_1_and_2loop_ss - tadpole_ss_2, 3.0e-11); - BOOST_CHECK_CLOSE(two_loop_tadpole[2] / vS, ts_1_and_2loop_ss - tadpole_ss_3, 1.0e-11); + // The following tests fail, because there are NaNs in the O(at*at) + // and O(atau*atau) 2L tadpoles. In this case Softsusy sets all + // tadpoles to 0. FlexibleSUSY sets only the ones to zero, which + // are NaN. + + // BOOST_CHECK_CLOSE(two_loop_tadpole[0] / vd, td_1_and_2loop_ss - tadpole_ss_1, 1.0e-10); + // BOOST_CHECK_CLOSE(two_loop_tadpole[1] / vu, tu_1_and_2loop_ss - tadpole_ss_2, 3.0e-11); + // BOOST_CHECK_CLOSE(two_loop_tadpole[2] / vS, ts_1_and_2loop_ss - tadpole_ss_3, 1.0e-11); } BOOST_AUTO_TEST_CASE( test_NMSSM_two_loop_ewsb ) diff --git a/test/test_NMSSM_low_scale_constraint.cpp b/test/test_NMSSM_low_scale_constraint.cpp index e67461a71..0fe5e7eba 100644 --- a/test/test_NMSSM_low_scale_constraint.cpp +++ b/test/test_NMSSM_low_scale_constraint.cpp @@ -129,17 +129,17 @@ BOOST_AUTO_TEST_CASE( test_low_energy_constraint ) // from the last run) to calculate the Yukawa couplings. BOOST_MESSAGE("testing diagonal yukawa elements"); - BOOST_CHECK_CLOSE_FRACTION(m.get_Yu()(0,0), s.displayYukawaMatrix(YU)(1,1), 0.002); - BOOST_CHECK_CLOSE_FRACTION(m.get_Yd()(0,0), s.displayYukawaMatrix(YD)(1,1), 0.002); - BOOST_CHECK_CLOSE_FRACTION(m.get_Ye()(0,0), s.displayYukawaMatrix(YE)(1,1), 0.002); + BOOST_CHECK_CLOSE_FRACTION(m.get_Yu()(0,0), s.displayYukawaMatrix(YU)(1,1), 0.005); + BOOST_CHECK_CLOSE_FRACTION(m.get_Yd()(0,0), s.displayYukawaMatrix(YD)(1,1), 0.005); + BOOST_CHECK_CLOSE_FRACTION(m.get_Ye()(0,0), s.displayYukawaMatrix(YE)(1,1), 0.005); - BOOST_CHECK_CLOSE_FRACTION(m.get_Yu()(1,1), s.displayYukawaMatrix(YU)(2,2), 0.002); - BOOST_CHECK_CLOSE_FRACTION(m.get_Yd()(1,1), s.displayYukawaMatrix(YD)(2,2), 0.002); - BOOST_CHECK_CLOSE_FRACTION(m.get_Ye()(1,1), s.displayYukawaMatrix(YE)(2,2), 0.002); + BOOST_CHECK_CLOSE_FRACTION(m.get_Yu()(1,1), s.displayYukawaMatrix(YU)(2,2), 0.005); + BOOST_CHECK_CLOSE_FRACTION(m.get_Yd()(1,1), s.displayYukawaMatrix(YD)(2,2), 0.005); + BOOST_CHECK_CLOSE_FRACTION(m.get_Ye()(1,1), s.displayYukawaMatrix(YE)(2,2), 0.005); - BOOST_CHECK_CLOSE_FRACTION(m.get_Yu()(2,2), s.displayYukawaMatrix(YU)(3,3), 0.002); - BOOST_CHECK_CLOSE_FRACTION(m.get_Yd()(2,2), s.displayYukawaMatrix(YD)(3,3), 0.002); - BOOST_CHECK_CLOSE_FRACTION(m.get_Ye()(2,2), s.displayYukawaMatrix(YE)(3,3), 0.002); + BOOST_CHECK_CLOSE_FRACTION(m.get_Yu()(2,2), s.displayYukawaMatrix(YU)(3,3), 0.005); + BOOST_CHECK_CLOSE_FRACTION(m.get_Yd()(2,2), s.displayYukawaMatrix(YD)(3,3), 0.005); + BOOST_CHECK_CLOSE_FRACTION(m.get_Ye()(2,2), s.displayYukawaMatrix(YE)(3,3), 0.005); BOOST_MESSAGE("testing running VEV"); const double running_vev = Sqrt(Sqr(m.get_vu()) + Sqr(m.get_vd())); From 6f6ae960d0b14af7375297dc8e0bbd273f8fa6df Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 2 Sep 2016 20:06:20 +0200 Subject: [PATCH 4/5] Bugfix: correcting Ben's catches for NaNs Fortran statements of the form if (condition) expression does not have an affect with my gfortran 4.9.2. --- src/mssm_twoloophiggs_impl.f | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/src/mssm_twoloophiggs_impl.f b/src/mssm_twoloophiggs_impl.f index fbad2167c..84694e0c2 100644 --- a/src/mssm_twoloophiggs_impl.f +++ b/src/mssm_twoloophiggs_impl.f @@ -5569,7 +5569,9 @@ function tauF1ab(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu) Nc = 1d0 mu2 = mu**2 - if(mu2.eq.0d0) mu2 = 1d-10 + if (mu2.eq.0d0) then + mu2 = 1d-10 + endif tauF1ab = $ (2.*BL*mu2/t*(mu2-BL+t)/delt(BL,mu2,t) @@ -5608,7 +5610,9 @@ function tauF1c(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu) Nc = 1d0 mu2 = mu**2 - if(mu2.eq.0d0) mu2 = 1d-10 + if (mu2.eq.0d0) then + mu2 = 1d-10 + endif Xt = s2t*(T1-T2)/2d0/Sqrt(t) Yt = s2t*(T1-T2)/2d0/Sqrt(t) - mu/sb/cb @@ -5721,7 +5725,9 @@ function tauF2c(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu) Nc = 1d0 mu2 = mu**2 - if(mu2.eq.0d0) mu2 = 1d-10 + if (mu2.eq.0d0) then + mu2 = 1d-10 + endif Xt = s2t*(T1-T2)/2d0/Sqrt(t) Yt = Xt - mu/cb/sb @@ -5837,7 +5843,9 @@ function tauF3c(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu) Nc = 1d0 mu2 = mu**2 - if(mu2.eq.0d0) mu2 = 1d-10 + if (mu2.eq.0d0) then + mu2 = 1d-10 + endif Xt = s2t*(T1-T2)/2d0/Sqrt(t) Yt = Xt - mu/cb/sb @@ -6316,12 +6324,16 @@ function tauFAc(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu) Nc = 1d0 mu2 = mu**2 - if(mu2.eq.0d0) mu2 = 1d-10 + if (mu2.eq.0d0) then + mu2 = 1d-10 + endif Xt = s2t*(T1-T2)/2d0/Sqrt(t) Yt = Xt - mu/cb/sb At = sb**2*Xt+cb**2*Yt - if (At.eq.0d0) At=1d-10 + if (At.eq.0d0) then + At=1d-10 + endif ct2 = (1d0+c2t)/2d0 st2 = (1d0-c2t)/2d0 @@ -6462,11 +6474,13 @@ subroutine tausqtad(t,A0,BL,T1,T2,st,ct,q,mu,tb,vv, c ADDED by BEN: guards against NANs when sin theta is zero! if (dabs(st).lt.1.0d-10) then - if (st.ge.0.0d0) st=1.0d-10 + if (st.ge.0.0d0) then + st=1.0d-10 else st=-1.0d-10 endif - ct = sqrt(1.d0-st*st) + endif + ct = sqrt(1.d0-st*st) c end of addition by BEN 13/6/12 s2t = 2d0*ct*st @@ -6523,12 +6537,15 @@ subroutine tautadfuncs(t,A0,BL,T1,T2,s2t,c2t,q,mu,sb,cb,F2l,G2l) cb2 = cb**2 mu2 = mu**2 - if(mu2.eq.0d0) mu2 = 1d-10 + if (mu2.eq.0d0) then + mu2 = 1d-10 + endif Xt = s2t*(T1-T2)/2d0/Sqrt(t) Yt = Xt - mu/cb/sb At = sb2*Xt+cb2*Yt + F2l = (delt(mu2,t,T1)+2*mu2*t)/T1*phi(mu2,t,T1) $ -(delt(mu2,t,T2)+2*mu2*t)/T2*phi(mu2,t,T2) $ +A0*cb2*(2*Sqrt(t)+s2t*Yt)/4/s2t/T1/(T1-T2)* From dd32401004c14c1880f1171fc913ee9dd8a81d8a Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Sat, 3 Sep 2016 00:37:20 +0200 Subject: [PATCH 5/5] ensure 4 characters wide unicode symbols --- meta/FSMathLink.m | 2 +- meta/Utils.m | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/meta/FSMathLink.m b/meta/FSMathLink.m index fdd00e09f..8a02964e1 100644 --- a/meta/FSMathLink.m +++ b/meta/FSMathLink.m @@ -134,7 +134,7 @@ ]; ToUTF8String[s_] := - StringJoin[("\\u0" <> IntegerString[#,16])& /@ ToCharacterCode[ToString[s]]]; + StringJoin[("\\u" <> Utils`FSStringPadLeft[IntegerString[#,16], 4, "0"])& /@ ToCharacterCode[ToString[s]]]; ToVaildWolframSymbolString[par_?CConversion`GreekQ] := ToUTF8String[par]; ToVaildWolframSymbolString[par_] := ToString[par]; diff --git a/meta/Utils.m b/meta/Utils.m index 8f21edfeb..0936894cd 100644 --- a/meta/Utils.m +++ b/meta/Utils.m @@ -84,6 +84,8 @@ occurrence of the given rule is replaced (if it exists) or added (if FSImportString::usage = "Returns the content of a file in form of a string. If the file does not exist, \"unknown\" is returned."; +FSStringPadLeft::usage = "StringPadLeft[] for Mathematica 9 and below."; + Begin["`Private`"]; ApplyAndConcatenate[Func_, l_List] := @@ -152,6 +154,9 @@ occurrence of the given rule is replaced (if it exists) or added (if ] ]; +FSStringPadLeft[str_String, width_, pad_String] := + StringJoin[PadLeft[Characters[str], width, pad]]; + ForceJoin[elem___] := Join[Sequence @@ Select[{elem}, (Head[#] === List)&]];