diff --git a/.github/workflows/cmake-format-linter.yml b/.github/workflows/cmake-format-linter.yml deleted file mode 100644 index ee9073aec..000000000 --- a/.github/workflows/cmake-format-linter.yml +++ /dev/null @@ -1,29 +0,0 @@ -name: Run cmake-format linter - -on: - pull_request: - branches: [ master, develop ] - -concurrency: - group: ${{ github.workflow }}-${{ github.ref }} - cancel-in-progress: true - -jobs: - build: - runs-on: ubuntu-latest - - steps: - - name: Checkout repository - uses: actions/checkout@v4 - - - name: Format CMake files - id: cmake-format - uses: PuneetMatharu/cmake-format-lint-action@v1.0.4 - with: - args: --config-files .cmake-format.py --in-place - - - name: Commit changes - uses: stefanzweifel/git-auto-commit-action@v5 - with: - commit_user_name: cmake-format-bot - commit_message: 'Automated commit of cmake-format changes.' \ No newline at end of file diff --git a/.github/workflows/codecov.yml b/.github/workflows/codecov.yml index 675020cf1..67a3aa4d0 100644 --- a/.github/workflows/codecov.yml +++ b/.github/workflows/codecov.yml @@ -61,7 +61,7 @@ jobs: - name: Generate Coverage run: cmake --build --preset $GeneratedCMakeProfile -j${{ steps.cpu-cores.outputs.count }} -t coverage - name: Upload coverage - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 with: token: ${{ secrets.BSMPT_CODECOV_UPLOAD_TOKEN }} flags: unittests # optional diff --git a/.github/workflows/doc.yml b/.github/workflows/doc.yml index 7f4d27b09..697bf3faa 100644 --- a/.github/workflows/doc.yml +++ b/.github/workflows/doc.yml @@ -54,7 +54,7 @@ jobs: name: documentation path: documentation - name: Deploy - uses: JamesIves/github-pages-deploy-action@v4.6.4 + uses: JamesIves/github-pages-deploy-action@v4.6.9 with: branch: gh-pages # The branch the action should deploy to. folder: documentation # The folder the action should deploy. diff --git a/.github/workflows/cpp-linter.yml b/.github/workflows/linter.yml similarity index 80% rename from .github/workflows/cpp-linter.yml rename to .github/workflows/linter.yml index 50b092f27..905fd7eca 100644 --- a/.github/workflows/cpp-linter.yml +++ b/.github/workflows/linter.yml @@ -1,21 +1,23 @@ -name: cpp-linter +name: Run all linters -on: +on: pull_request: - branches: [master] - + branches: [ master, develop ] + concurrency: - group: ${{ github.workflow }}-${{ github.ref }} - cancel-in-progress: true - + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: - cpp-linter: + linter: runs-on: ubuntu-latest + steps: - uses: actions/checkout@v4 with: fetch-depth: 0 ref: ${{ github.head_ref }} + token: ${{ secrets.LINTER_PAT }} - name: Install clang-format if: inputs.apply_clang_format @@ -28,6 +30,12 @@ jobs: run: echo "branchName=$GITHUB_BASE_REF" >> $GITHUB_OUTPUT + - name: Format CMake files + id: cmake-format + uses: PuneetMatharu/cmake-format-lint-action@v1.0.4 + with: + args: --config-files .cmake-format.py --in-place + - name: Define base git diff args id: git-diff-args run: | diff --git a/include/BSMPT/minimum_tracer/minimum_tracer.h b/include/BSMPT/minimum_tracer/minimum_tracer.h index eafb84ce0..d8fc28f79 100644 --- a/include/BSMPT/minimum_tracer/minimum_tracer.h +++ b/include/BSMPT/minimum_tracer/minimum_tracer.h @@ -580,24 +580,6 @@ Create1DimGrid(const std::vector &min_start, const std::vector &min_end, const int npoints = 100); -/** - * Returns true if two values are the same given some relative precision - */ -bool almost_the_same(const double &a, - const double &b, - const double &rel_precision = 0.01, - const double &num_zero = 1e-10); - -/** - * Returns true if two vectors are the element-wise the same given some relative - * precision - */ -bool almost_the_same(const std::vector &a, - const std::vector &b, - const bool &allow_for_sign_flip = false, - const double &rel_precision = 0.01, - const double &num_zero = 1e-10); - /** * @brief Phase object * @@ -950,4 +932,4 @@ struct Vacuum void PrintPhasesDiagram(int size = 100); }; -} // namespace BSMPT \ No newline at end of file +} // namespace BSMPT diff --git a/include/BSMPT/models/ClassPotentialC2HDM.h b/include/BSMPT/models/ClassPotentialC2HDM.h index d2b50763d..dfe1e1997 100644 --- a/include/BSMPT/models/ClassPotentialC2HDM.h +++ b/include/BSMPT/models/ClassPotentialC2HDM.h @@ -98,6 +98,10 @@ class Class_Potential_C2HDM : public Class_Potential_Origin double CTempC1 = 0, CTempC2 = 0, CTempCS = 0; double R_Hh_1 = 0, R_Hh_2 = 0, R_Hh_3 = 0, R_Hl_1 = 0, R_Hl_2 = 0, R_Hl_3 = 0, R_Hsm_1 = 0, R_Hsm_2 = 0, R_Hsm_3 = 0; + + int pos_G0, pos_G1, pos_G2, pos_H1, pos_H2, pos_h1, pos_h2, pos_h3; + int pos_h_SM, pos_h_l, pos_h_H; + void ReadAndSet(const std::string &linestr, std::vector &par) override; std::vector addLegendCT() const override; @@ -121,6 +125,7 @@ class Class_Potential_C2HDM : public Class_Potential_Origin void set_CT_Pot_Par(const std::vector &par) override; void write() const override; + void AdjustRotationMatrix() override; void TripleHiggsCouplings() override; std::vector calc_CT() const override; diff --git a/include/BSMPT/models/ClassPotentialCPintheDark.h b/include/BSMPT/models/ClassPotentialCPintheDark.h index dc9aff646..8673b833d 100644 --- a/include/BSMPT/models/ClassPotentialCPintheDark.h +++ b/include/BSMPT/models/ClassPotentialCPintheDark.h @@ -95,6 +95,9 @@ class Class_Potential_CPintheDark : public Class_Potential_Origin // vev double v1; + int pos_Gp, pos_Gm, pos_Hp, pos_Hm, pos_HSM; + int pos_G0, pos_h1, pos_h2, pos_h3; + void ReadAndSet(const std::string &linestr, std::vector &par) override; std::vector addLegendCT() const override; @@ -106,6 +109,7 @@ class Class_Potential_CPintheDark : public Class_Potential_Origin void set_CT_Pot_Par(const std::vector &par) override; void write() const override; + void AdjustRotationMatrix() override; void TripleHiggsCouplings() override; std::vector calc_CT() const override; diff --git a/include/BSMPT/models/ClassPotentialCxSM.h b/include/BSMPT/models/ClassPotentialCxSM.h index 80fa16b36..e51f54929 100644 --- a/include/BSMPT/models/ClassPotentialCxSM.h +++ b/include/BSMPT/models/ClassPotentialCxSM.h @@ -83,6 +83,10 @@ class Class_CxSM : public Class_Potential_Origin double vh, vs, va; + int pos_Gp, pos_Gm, pos_G0; + int pos_h1, pos_h2, pos_h3; + int pos_h_SM, pos_h_l, pos_h_H; + void ReadAndSet(const std::string &linestr, std::vector &par) override; std::vector addLegendCT() const override; @@ -109,6 +113,7 @@ class Class_CxSM : public Class_Potential_Origin void set_CT_Pot_Par(const std::vector &par) override; void write() const override; + void AdjustRotationMatrix() override; void TripleHiggsCouplings() override; std::vector calc_CT() const override; diff --git a/include/BSMPT/models/ClassPotentialN2HDM.h b/include/BSMPT/models/ClassPotentialN2HDM.h index 80da9df5f..da629186a 100644 --- a/include/BSMPT/models/ClassPotentialN2HDM.h +++ b/include/BSMPT/models/ClassPotentialN2HDM.h @@ -99,6 +99,10 @@ class Class_Potential_N2HDM : public Class_Potential_Origin double NDus = 0, NDL6 = 0, NDL7 = 0, NDL8 = 0, NDvs = 0, NDTS = 0; double DTCharged = 0; + int pos_G0, pos_G1, pos_G2, pos_H1, pos_H2; + int pos_h1, pos_h2, pos_h3, pos_A; + int pos_h_SM, pos_h_l, pos_h_H; + void ReadAndSet(const std::string &linestr, std::vector &par) override; std::vector addLegendCT() const override; @@ -125,6 +129,7 @@ class Class_Potential_N2HDM : public Class_Potential_Origin void set_CT_Pot_Par(const std::vector &par) override; void write() const override; + void AdjustRotationMatrix() override; void TripleHiggsCouplings() override; std::vector calc_CT() const override; diff --git a/include/BSMPT/models/ClassPotentialOrigin.h b/include/BSMPT/models/ClassPotentialOrigin.h index 4c791f3ff..9bff3b634 100644 --- a/include/BSMPT/models/ClassPotentialOrigin.h +++ b/include/BSMPT/models/ClassPotentialOrigin.h @@ -182,10 +182,10 @@ class Class_Potential_Origin */ bool SetCurvatureDone = false; /** - * @brief CalcCouplingsdone Used to check if CalculatePhysicalCouplings has + * @brief CalcCouplingsDone Used to check if CalculatePhysicalCouplings has * already been called */ - bool CalcCouplingsdone = false; + bool CalcCouplingsDone = false; /** * @brief CalculatedTripleCopulings Used to check if TripleHiggsCouplings has * already been called @@ -287,6 +287,11 @@ class Class_Potential_Origin * tree-level Vacuum */ std::vector> HiggsRotationMatrix; + /** + * Storage of the model-specific Higgs rotation matrix for the Higgs mass + * matrix at the tree-level Vacuum + */ + std::vector> HiggsRotationMatrixEnsuredConvention; /** * @brief Couplings_Higgs_Quartic Stores the quartic Higgs couplings in the * mass base @@ -969,6 +974,14 @@ class Class_Potential_Origin */ Eigen::MatrixXcd LeptonMassMatrix(const std::vector &v) const; + /** + * Ensures the correct rotation matrix convention + */ + virtual void AdjustRotationMatrix() = 0; + /** + * Checks whether rotation matrix is properly set after implying conventions + */ + bool CheckRotationMatrix(); /** * Calculates the triple Higgs couplings at NLO in the mass basis. * diff --git a/include/BSMPT/models/ClassPotentialR2HDM.h b/include/BSMPT/models/ClassPotentialR2HDM.h index 9f8bf3caa..c80167930 100644 --- a/include/BSMPT/models/ClassPotentialR2HDM.h +++ b/include/BSMPT/models/ClassPotentialR2HDM.h @@ -90,6 +90,9 @@ class Class_Potential_R2HDM : public Class_Potential_Origin int Type = 0; double CTempC1 = 0, CTempC2 = 0, CTempCS = 0; + int pos_G1, pos_G2, pos_H1, pos_H2, pos_G0, pos_A, pos_H, pos_h; + int pos_h_SM, pos_h_H; + void ReadAndSet(const std::string &linestr, std::vector &par) override; std::vector addLegendCT() const override; @@ -112,6 +115,7 @@ class Class_Potential_R2HDM : public Class_Potential_Origin void set_CT_Pot_Par(const std::vector &par) override; void write() const override; + void AdjustRotationMatrix() override; void TripleHiggsCouplings() override; std::vector calc_CT() const override; diff --git a/include/BSMPT/models/ClassPotentialSM.h b/include/BSMPT/models/ClassPotentialSM.h index c657e8c24..9574fe866 100644 --- a/include/BSMPT/models/ClassPotentialSM.h +++ b/include/BSMPT/models/ClassPotentialSM.h @@ -83,6 +83,7 @@ class Class_SM : public Class_Potential_Origin void set_CT_Pot_Par(const std::vector &par) override; void write() const override; + void AdjustRotationMatrix() override; void TripleHiggsCouplings() override; std::vector calc_CT() const override; diff --git a/include/BSMPT/models/ClassTemplate.h b/include/BSMPT/models/ClassTemplate.h index d596eca37..7f5c6ad0c 100644 --- a/include/BSMPT/models/ClassTemplate.h +++ b/include/BSMPT/models/ClassTemplate.h @@ -48,6 +48,7 @@ class Class_Template : public Class_Potential_Origin void set_CT_Pot_Par(const std::vector &par) override; void write() const override; + void AdjustRotationMatrix() override; void TripleHiggsCouplings() override; std::vector calc_CT() const override; diff --git a/include/BSMPT/utility/utility.h b/include/BSMPT/utility/utility.h index 8d06ce5b2..374a0fd55 100644 --- a/include/BSMPT/utility/utility.h +++ b/include/BSMPT/utility/utility.h @@ -15,6 +15,8 @@ #include #include #include +#include + #ifdef Boost_FOUND #include @@ -284,6 +286,27 @@ double Li2(const double &x); */ double EllipIntSecond(const double &x); +/** + * @brief Checks if two double numbers are (almost) the same with a given + * relative precision; if both numbers are smaller than num_zero, then they + * are considered to be zero and the function always returns true; with + * additional versions of the function for std::complex and + * std::vector input. + */ +bool almost_the_same(const double &a, + const double &b, + const double &rel_precision = 0.01, + const double &num_zero = 1e-10); +bool almost_the_same(const std::complex &a, + const std::complex &b, + const double &rel_precision = 0.01, + const double &num_zero = 1e-10); +bool almost_the_same(const std::vector &a, + const std::vector &b, + const bool &allow_for_sign_flip = false, + const double &rel_precision = 0.01, + const double &num_zero = 1e-10); + /** * @brief operator << overload for the model parameter */ diff --git a/src/minimum_tracer/minimum_tracer.cpp b/src/minimum_tracer/minimum_tracer.cpp index e806ad2f7..a705f5e54 100644 --- a/src/minimum_tracer/minimum_tracer.cpp +++ b/src/minimum_tracer/minimum_tracer.cpp @@ -1425,52 +1425,6 @@ Create1DimGrid(const std::vector &min_start, return res_vec; } -bool almost_the_same(const double &a, - const double &b, - const double &rel_precision, - const double &num_zero) -{ - if (std::abs(a) < num_zero and std::abs(b) < num_zero) - { - return true; - } - return std::abs(a - b) < std::abs(a + b) / 2 * rel_precision; -} - -bool almost_the_same(const std::vector &a, - const std::vector &b, - const bool &allow_for_sign_flip, - const double &rel_precision, - const double &num_zero) -{ - if (a.size() != b.size()) - { - throw std::runtime_error("Error. Vectors must have the same size."); - } - int count_true = 0; - for (std::size_t i = 0; i < a.size(); i++) - { - if (allow_for_sign_flip) - { - count_true += - int(almost_the_same(a.at(i), b.at(i), rel_precision, num_zero)); - } - else - { - count_true += int(almost_the_same( - std::abs(a.at(i)), std::abs(b.at(i)), rel_precision, num_zero)); - } - } - if (std::size_t(count_true) == a.size()) - { - return true; - } - else - { - return false; - } -} - Phase::Phase() { } diff --git a/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index e1e9aac39..3c1b29b70 100644 --- a/src/models/ClassPotentialC2HDM.cpp +++ b/src/models/ClassPotentialC2HDM.cpp @@ -2088,8 +2088,6 @@ void Class_Potential_C2HDM::write() const std::stringstream ss; ss.precision(std::numeric_limits::max_digits10); - double MSM = 0, MhUp = 0, MhDown = 0; - ss << "scale = " << scale << std::endl; ss << "The parameters are : \n"; @@ -2128,97 +2126,27 @@ void Class_Potential_C2HDM::write() const ss << "DT3:= " << DT3 << ";\n"; ss << "DIL6:= " << DIL6CT << ";\n"; - if (CalcCouplingsdone) + if (CalcCouplingsDone) { MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; - } - } - - int posMHCS1 = 0; - int posN[3]; - int countposN = 0; - int posG0 = 0; - double testsum = 0; - const double ZeroThreshold = 1e-5; - for (int i = 0; i < 3; i++) - { - // testsum = std::abs(HiggsRot(i,0)) + std::abs(HiggsRot(i,2)); - // if(testsum > ZeroThreshold) posG1 = i; - // testsum = std::abs(HiggsRot(i,1)) + std::abs(HiggsRot(i,3)); - // if(testsum > ZeroThreshold) posG2 = i; - testsum = std::abs(HiggsRot(i, 5)) + std::abs(HiggsRot(i, 7)); - if (testsum > ZeroThreshold) posG0 = i; - } - for (std::size_t i = 3; i < NHiggs; i++) - { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posMHCS1 = i; - // testsum = std::abs(HiggsRot(i,1)) + std::abs(HiggsRot(i,3)); - // if(testsum > ZeroThreshold) posMHCS2 = i; - testsum = 0; - for (int k = 4; k < 8; k++) - testsum += std::abs(HiggsRot(i, k)); - if (testsum > ZeroThreshold) - { - posN[countposN] = i; - countposN++; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } std::vector HiggsMasses; HiggsMasses = HiggsMassesSquared(vevTree, 0); - double NeutralHiggs[3]; - for (int i = 0; i < 3; i++) - { - NeutralHiggs[i] = HiggsMasses[posN[i]]; - } - for (int i = 0; i < 3; i++) - { - if (std::sqrt(NeutralHiggs[i]) < 126 and std::sqrt(NeutralHiggs[i]) > 124) - MSM = std::sqrt(NeutralHiggs[i]); - } - if (std::sqrt(NeutralHiggs[0]) == MSM) - { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[1]); - } - else if (std::sqrt(NeutralHiggs[1]) == MSM) - { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[0]); - } - else - { - MhUp = std::sqrt(NeutralHiggs[1]); - MhDown = std::sqrt(NeutralHiggs[0]); - } - - if (MSM > MhUp) - { - double tmp = posN[1]; - posN[1] = posN[2]; - posN[2] = tmp; - } - if (MSM > MhDown) - { - double tmp = posN[0]; - posN[0] = posN[1]; - posN[1] = tmp; - } - MatrixXd NeutralMatrix(4, 4); for (int j = 0; j < 4; j++) - NeutralMatrix(0, j) = HiggsRot(posG0, j + 4); - for (int i = 1; i < 4; i++) { - for (int j = 0; j < 4; j++) - NeutralMatrix(i, j) = HiggsRot(posN[i - 1], j + 4); + NeutralMatrix(0, j) = HiggsRot(pos_G0, j + 4); + NeutralMatrix(1, j) = HiggsRot(pos_h_SM, j + 4); + NeutralMatrix(2, j) = HiggsRot(pos_h_l, j + 4); + NeutralMatrix(3, j) = HiggsRot(pos_h_H, j + 4); } MatrixXd MassMixing(3, 3); @@ -2231,10 +2159,10 @@ void Class_Potential_C2HDM::write() const } ss << "The mass spectrum is given by :\n"; - ss << "m_{H^+} = " << std::sqrt(HiggsMasses[posMHCS1]) << " GeV \n" - << "m_{H_SM} = " << MSM << " GeV \n" - << "m_{H_l} = " << MhDown << " GeV \n" - << "m_{H_h} = " << MhUp << " GeV \n"; + ss << "m_{H^+} = " << std::sqrt(HiggsMasses[pos_H1]) << " GeV \n" + << "m_{H_SM} = " << std::sqrt(HiggsMasses[pos_h_SM]) << " GeV \n" + << "m_{H_l} = " << std::sqrt(HiggsMasses[pos_h_l]) << " GeV \n" + << "m_{H_h} = " << std::sqrt(HiggsMasses[pos_h_H]) << " GeV \n"; ss << "The neutral mixing Matrix is given by :\n"; bool IsNegative = MassMixing(0, 1) < 0; ss << "H_{SM} = " << MassMixing(0, 0) << " zeta_1 "; @@ -2294,7 +2222,7 @@ std::vector Class_Potential_C2HDM::calc_CT() const retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; @@ -2396,30 +2324,21 @@ std::vector Class_Potential_C2HDM::calc_CT() const return parCT; } -void Class_Potential_C2HDM::TripleHiggsCouplings() + +/** + * Ensures the correct rotation matrix convention + */ +void Class_Potential_C2HDM::AdjustRotationMatrix() { - if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + const double ZeroThreshold = 1e-5; - if (CalculatedTripleCopulings) return; - CalculatedTripleCopulings = true; + if (!SetCurvatureDone) SetCurvatureArrays(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); - std::vector TripleDeriv; - TripleDeriv = WeinbergThirdDerivative(); - std::vector>> GaugeBasis( - NHiggs, - std::vector>(NHiggs, std::vector(NHiggs))); - // double GaugeBasis[NHiggs][NHiggs][NHiggs]; - for (std::size_t i = 0; i < NHiggs; i++) + if (!CheckRotationMatrix()) // Check whether generically generated rotation + // matrix is proper rotation matrix { - for (std::size_t j = 0; j < NHiggs; j++) - { - for (std::size_t k = 0; k < NHiggs; k++) - { - GaugeBasis[i][j][k] = - TripleDeriv.at(i + j * NHiggs + k * NHiggs * NHiggs); - } - } + throw std::runtime_error("Error in rotation matrix."); } MatrixXd HiggsRot(NHiggs, NHiggs); @@ -2430,99 +2349,480 @@ void Class_Potential_C2HDM::TripleHiggsCouplings() HiggsRot(i, j) = HiggsRotationMatrix[i][j]; } } - MatrixXd HiggsRotSort(NHiggs, NHiggs); - int posMHCS1 = 0, posMHCS2 = 0; - std::vector posN(3); - std::size_t countposN = 0; - std::size_t posG1 = 0, posG2 = 0, posG0 = 0; - double testsum = 0; - const double ZeroThreshold = 1e-5; - for (int i = 0; i < 3; i++) + + // C2HDM interaction basis + // rho1, eta1, rho2, eta2, zeta1, psi1, zeta2, psi2 + int pos_rho1 = 0, pos_eta1 = 1, pos_rho2 = 2, pos_eta2 = 3, pos_zeta1 = 4, + pos_psi1 = 5, pos_zeta2 = 6, pos_psi2 = 7; + + // Indices of mass eigenstates for rotation from interaction to mass basis + pos_G0 = -1, pos_G1 = -1, pos_G2 = -1, pos_H1 = -1, pos_H2 = -1, + pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; + + // Going from 0 to 2, i.e. fixing the Goldstone indices first + // (using that the Goldstone masses appear before the physical Higgs bosons + // since they are the smallest mass eigenvalues (= 0 in the Landau gauge)) + for (std::size_t i = 0; i < 3; i++) + // mass base index i corresponds to mass vector sorted in ascending mass { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posG1 = i; - testsum = std::abs(HiggsRot(i, 1)) + std::abs(HiggsRot(i, 3)); - if (testsum > ZeroThreshold) posG2 = i; - testsum = std::abs(HiggsRot(i, 5)) + std::abs(HiggsRot(i, 7)); - if (testsum > ZeroThreshold) posG0 = i; + // Charged submatrix + if (std::abs(HiggsRot(i, pos_rho1)) + + std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) + { + pos_G1 = i; + } + else if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_G2 = i; + } + // Neutral submatrix + else if (std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) + { + pos_G0 = i; + } + else + { + throw std::runtime_error("Error. Non-Goldstone in Goldstone submatrix." + " Particle with negative mass? Check your " + "parameter point."); + } } + + // Now going from 3 to NHiggs, i.e. fixing the physical Higgs bosons for (std::size_t i = 3; i < NHiggs; i++) + // mass base index i corresponds to mass vector sorted in ascending mass { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posMHCS1 = i; - testsum = std::abs(HiggsRot(i, 1)) + std::abs(HiggsRot(i, 3)); - if (testsum > ZeroThreshold) posMHCS2 = i; - testsum = 0; - for (int k = 4; k < 8; k++) - testsum += std::abs(HiggsRot(i, k)); - if (testsum > ZeroThreshold) + // Charged submatrix + if (std::abs(HiggsRot(i, pos_rho1)) + + std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) + { + pos_H1 = i; + } + else if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) { - posN.at(countposN) = i; - countposN++; + pos_H2 = i; + } + // Neutral CP-mixed submatrix + else if (std::abs(HiggsRot(i, pos_zeta1)) + + std::abs(HiggsRot(i, pos_zeta2)) + + std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) + // use that mh1 < mh2 < mh3 + { + if (pos_h1 == -1) { + pos_h1 = i; + } else if (pos_h2 == -1) { + pos_h2 = i; + } else { + pos_h3 = i; + } } } - if (UseHsmNotationInTripleHiggs) + // check if all position indices are set + if (pos_G0 == -1 or pos_G1 == -1 or pos_G2 == -1 or + pos_H1 == -1 or pos_H2 == -1 or + pos_h1 == -1 or pos_h2 == -1 or pos_h3 == -1) { - std::vector HiggsMasses; - double MhUp = 0, MhDown = 0, MSM = 0; - HiggsMasses = HiggsMassesSquared(vevTree, 0); + throw std::runtime_error("Error. Not all position indices are set."); + } - double NeutralHiggs[3]; - for (int i = 0; i < 3; i++) + // check if all other elements of rotation matrix are zero + bool zero_element = false; + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) { - NeutralHiggs[i] = HiggsMasses[posN.at(i)]; + int ii = int(i); + int jj = int(j); + if (not((jj == pos_rho1 and (ii == pos_G1 or ii == pos_H1)) or + (jj == pos_rho2 and (ii == pos_G1 or ii == pos_H1)) or + (jj == pos_eta1 and (ii == pos_G2 or ii == pos_H2)) or + (jj == pos_eta2 and (ii == pos_G2 or ii == pos_H2)) or + (jj == pos_psi1 and (ii == pos_G0 or + ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_psi2 and (ii == pos_G0 or + ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_zeta1 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_zeta2 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)))) + { + zero_element = true; + } + + if (zero_element and std::abs(HiggsRot(i, j)) > ZeroThreshold) + { + throw std::runtime_error("Error. Invalid rotation matrix detected."); + } + zero_element = false; } - for (int i = 0; i < 3; i++) + } + + // Determine the additional indices for the SM-like + // and lighter/heavier Higgses + pos_h_SM = -1, pos_h_l = -1, pos_h_H = -1; + + std::vector HiggsMasses; + HiggsMasses = HiggsMassesSquared(vevTree, 0); + + // Due to the masses being ordered, we will always have + // HiggsMasses[pos_h1] <= HiggsMasses[pos_h2] <= HiggsMasses[pos_h3] + double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h1]) + - SMConstants.C_MassSMHiggs); + double diff2 = std::abs(std::sqrt(HiggsMasses[pos_h2]) + - SMConstants.C_MassSMHiggs); + double diff3 = std::abs(std::sqrt(HiggsMasses[pos_h3]) + - SMConstants.C_MassSMHiggs); + if (diff1 < diff2 and diff1 < diff3) + { + pos_h_SM = pos_h1; + pos_h_l = pos_h2; + pos_h_H = pos_h3; + } + else if (diff2 < diff1 and diff2 < diff3) + { + pos_h_l = pos_h1; + pos_h_SM = pos_h2; + pos_h_H = pos_h3; + } + else + { + pos_h_l = pos_h1; + pos_h_H = pos_h2; + pos_h_SM = pos_h3; + } + + // Steps: + // (1) Rotate mass matrix from interaction to semi-interaction basis + // (i.e. interaction basis with neutral Goldstone rotated out): + // + // From interaction basis + // 0 1 2 3 4 5 6 7 + // rho1, eta1, rho2, eta2, zeta1, psi1, zeta2, psi2 + // to semi-interaction basis (Goldstone rotated out) + // 0 1 2 3 4 5 6 7 + // G^0, rho1, eta1, rho2, eta2, zeta1, zeta2, zeta3 + // + // (2) Diagonalise mass matrix in semi-interaction basis + // -> obtain rotation matrix from semi-interaction to mass basis + // + // (3) Neutral part of rotation matrix from step 2 must have the 3x3 form + // as in arXiv:1803.02846 Eq. (3.91), so this is the one to check for + // R11 > 0, R33 > 0, det > 0 (see arXiv:2007.02985 Eq. (6)) + + // Find position of neutral Goldstone: + // * Mass eigenvalues are ordered from smallest to largest + // => First three rows correspond to Goldstones + // (= massless in Landau gauge) + // * Charged and neutral Goldstones do not mix + // => Look for row which has psi1 and psi2 mixing components =/= 0, + // and the rest = 0 + + // Matrix to "rotate out" the neutral Goldstone boson, see arXiv:1803.02846 Eq. (3.89) + MatrixXd RotGoldstone(NHiggs, NHiggs); + RotGoldstone.row(0) << 0., 0., 0., 0., 0., C_CosBeta, 0., C_SinBeta; + RotGoldstone.row(1) << 1., 0., 0., 0., 0., 0., 0., 0.; + RotGoldstone.row(2) << 0., 1., 0., 0., 0., 0., 0., 0.; + RotGoldstone.row(3) << 0., 0., 1., 0., 0., 0., 0., 0.; + RotGoldstone.row(4) << 0., 0., 0., 1., 0., 0., 0., 0.; + RotGoldstone.row(5) << 0., 0., 0., 0., 1., 0., 0., 0.; + RotGoldstone.row(6) << 0., 0., 0., 0., 0., 0., 1., 0.; + RotGoldstone.row(7) << 0., 0., 0., 0., 0., -C_SinBeta, 0., C_CosBeta; + + // Swap rows to ensure that G0 is always in the first row + // for the following manipulations + MatrixXd MoveGoldstoneFirst(NHiggs, NHiggs); + MoveGoldstoneFirst.setIdentity(NHiggs, NHiggs); + if (pos_G0 != 0) + { + MoveGoldstoneFirst(0, 0) = 0.; + MoveGoldstoneFirst(pos_G0, pos_G0) = 0.; + MoveGoldstoneFirst(0, pos_G0) = 1.; + MoveGoldstoneFirst(pos_G0, 0) = 1.; + } + + // Compute rotation matrix from the "semi-interaction" (with G0 rotated out) to + // the mass basis, to get the same rotation as in arXiv:1803.02846 Eqs. (3.90)-(3.91) + MatrixXd RotGoldstoneMassBasis(NHiggs, NHiggs); + RotGoldstoneMassBasis = MoveGoldstoneFirst*HiggsRot*RotGoldstone.transpose(); + + // Semi-interaction basis (neutral Goldstone rotated out) + // G^0, rho1, eta1, rho2, eta2, zeta1, zeta2, zeta3 + int pos_si_G0 = 0, pos_si_rho1 = 1, pos_si_eta1 = 2, + pos_si_rho2 = 3, pos_si_eta2 = 4, + pos_si_zeta1 = 5, pos_si_zeta2 = 6, pos_si_zeta3 = 7; + + double row1 = 0.0, col1 = 0.0; + // Sum only over index starting from 1 (i.e. don't include the (0,0) element, + // i.e. the upper left element, of the matrix); the sum should be very small + for (std::size_t i = 1; i < NHiggs; i++) + { + row1 += std::abs(RotGoldstoneMassBasis(0, i)); + col1 += std::abs(RotGoldstoneMassBasis(i, 0)); + } + + // Consistency check that the Goldstone was rotated out properly: + // first row/column should contain only zeroes except for the upper left element + if (std::abs(std::abs(RotGoldstoneMassBasis(0, 0)) - 1.0) > ZeroThreshold + or std::abs(row1) > ZeroThreshold or std::abs(col1) > ZeroThreshold) + { + throw std::runtime_error("Error. Something went wrong after rotating " + "out the neutral Goldstone."); + } + + // Indices of mass eigenstates for rotation from semi-interaction to mass basis; + // position of neutral Goldstone is fixed to 0, see above + int pos_si_G1 = -1, pos_si_G2 = -1, pos_si_H1 = -1, pos_si_H2 = -1; + int pos_si_h1 = -1, pos_si_h2 = -1, pos_si_h3 = -1; + + // Start with i = 1, i.e. skip over the neutral Goldstone + for (std::size_t i = 1; i < NHiggs; i++) + // mass base index i corresponds to mass vector sorted in ascending mass + { + // Charged submatrices + // Check if the field with index i has a rho1 or rho2 component; + // since Goldstone mass is zero (Landau gauge), it appears before the charged Higgs + if (std::abs(RotGoldstoneMassBasis(i, pos_si_rho1)) + + std::abs(RotGoldstoneMassBasis(i, pos_si_rho2)) > ZeroThreshold) + // use that 0 = mGpm < mHpm { - if (std::sqrt(NeutralHiggs[i]) < 126 and std::sqrt(NeutralHiggs[i]) > 124) - MSM = std::sqrt(NeutralHiggs[i]); + if (pos_si_G1 == -1) + { + pos_si_G1 = i; + } + else + { + pos_si_H1 = i; + } } - if (std::sqrt(NeutralHiggs[0]) == MSM) + else if (std::abs(RotGoldstoneMassBasis(i, pos_si_eta1)) + + std::abs(RotGoldstoneMassBasis(i, pos_si_eta2)) > ZeroThreshold) + // use that 0 = mGpm < mHpm { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[1]); + if (pos_si_G2 == -1) + { + pos_si_G2 = i; + } + else + { + pos_si_H2 = i; + } } - else if (std::sqrt(NeutralHiggs[1]) == MSM) + // Neutral submatrix (mixed CP-even and CP-odd states); + // neutral Goldstone already rotated out + else if (std::abs(RotGoldstoneMassBasis(i, pos_si_zeta1)) + + std::abs(RotGoldstoneMassBasis(i, pos_si_zeta2)) + + std::abs(RotGoldstoneMassBasis(i, pos_si_zeta3)) > ZeroThreshold) + // use that mh1 < mh2 < mh3 { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[0]); + if (pos_si_h1 == -1) { + pos_si_h1 = i; + } else if (pos_si_h2 == -1) { + pos_si_h2 = i; + } else { + pos_si_h3 = i; + } } - else + } + + // Check if all position indices are set + if (pos_si_G1 == -1 or pos_si_G2 == -1 or + pos_si_H1 == -1 or pos_si_H2 == -1 or + pos_si_h1 == -1 or pos_si_h2 == -1 or pos_si_h3 == -1) + { + throw std::runtime_error("Error. Not all position indices are set."); + } + + // Check if all other elements of rotation matrix are zero + zero_element = false; + // Start with i, j = 1, skip neutral Goldstone + for (std::size_t i = 1; i < NHiggs; i++) + { + for (std::size_t j = 1; j < NHiggs; j++) { - MhUp = std::sqrt(NeutralHiggs[1]); - MhDown = std::sqrt(NeutralHiggs[0]); + int ii = int(i); + int jj = int(j); + if (not((jj == pos_si_rho1 and (ii == pos_si_G1 or ii == pos_si_H1)) or + (jj == pos_si_eta1 and (ii == pos_si_G2 or ii == pos_si_H2)) or + (jj == pos_si_rho2 and (ii == pos_si_G1 or ii == pos_si_H1)) or + (jj == pos_si_eta2 and (ii == pos_si_G2 or ii == pos_si_H2)) or + (jj == pos_si_zeta1 and (ii == pos_si_h1 or ii == pos_si_h2 or ii == pos_si_h3)) or + (jj == pos_si_zeta2 and (ii == pos_si_h1 or ii == pos_si_h2 or ii == pos_si_h3)) or + (jj == pos_si_zeta3 and (ii == pos_si_h1 or ii == pos_si_h2 or ii == pos_si_h3)))) + { + zero_element = true; + } + + if (zero_element and std::abs(RotGoldstoneMassBasis(i, j)) > ZeroThreshold) + { + throw std::runtime_error("Error. Invalid rotation matrix detected."); + } + zero_element = false; } + } - if (MSM > MhUp) + MatrixXd HiggsRotFixed(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixed.row(i) = RotGoldstoneMassBasis.row(i); + } + + // Neutral Goldstone; flip sign if its element, which should be "1", is negative + if (HiggsRotFixed(pos_si_G0, pos_si_G0) < 0) // G0 G0 (+1) + { + HiggsRotFixed.row(pos_si_G0) *= -1; + } + + // charged submatrix + if (HiggsRotFixed(pos_si_G1, pos_si_rho1) < 0) // G1 rho1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_si_G1) *= -1; + } + if (HiggsRotFixed(pos_si_G2, pos_si_eta1) < 0) // G2 eta1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_si_G2) *= -1; + } + if (HiggsRotFixed(pos_si_H1, pos_si_rho2) < 0) // H1 rho2 (+ cos(beta)) + { + HiggsRotFixed.row(pos_si_H1) *= -1; + } + if (HiggsRotFixed(pos_si_H2, pos_si_eta2) < 0) // H2 eta2 (+ cos(beta)) + { + HiggsRotFixed.row(pos_si_H2) *= -1; + } + + // Check neutral submatrix + // Use the "ScannerS" criteria from arXiv:2007.02985 Eq. (6) + // (since they use the same parametrisation of the angles as BSMPT): + // * (1) if R[1][1] < 0: h1 -> -h1 (i.e. multiply the h1 row with -1) + // * (2) if R[3][3] < 0: h3 -> -h3 (i.e. multiply the h3 row with -1) + // * (3) if det R < 0: h2 -> -h2 (i.e. multiply the h2 row with -1) + + // check neutral, CP-even submatrix + if (HiggsRotFixed(pos_si_h1, pos_si_zeta1) < 0) + // h1 zeta1 (condition (1) above, R11 < 0) + { + // if negative, flip sign of h1 + HiggsRotFixed.row(pos_si_h1) *= -1; + } + + if (HiggsRotFixed(pos_si_h3, pos_si_zeta3) < 0) + // h3 zeta3 (condition (2) above, R33 < 0) + { + // if negative, flip sign of h3 + HiggsRotFixed.row(pos_si_h3) *= -1; + } + + // Calculate the determinant AFTER flipping the signs for rows 1 and 3 above + MatrixXd HiggsRotFixedNeutral(3, 3); + HiggsRotFixedNeutral(0, 0) = HiggsRotFixed(pos_si_h1, pos_si_zeta1); + HiggsRotFixedNeutral(0, 1) = HiggsRotFixed(pos_si_h1, pos_si_zeta2); + HiggsRotFixedNeutral(0, 2) = HiggsRotFixed(pos_si_h1, pos_si_zeta3); + + HiggsRotFixedNeutral(1, 0) = HiggsRotFixed(pos_si_h2, pos_si_zeta1); + HiggsRotFixedNeutral(1, 1) = HiggsRotFixed(pos_si_h2, pos_si_zeta2); + HiggsRotFixedNeutral(1, 2) = HiggsRotFixed(pos_si_h2, pos_si_zeta3); + + HiggsRotFixedNeutral(2, 0) = HiggsRotFixed(pos_si_h3, pos_si_zeta1); + HiggsRotFixedNeutral(2, 1) = HiggsRotFixed(pos_si_h3, pos_si_zeta2); + HiggsRotFixedNeutral(2, 2) = HiggsRotFixed(pos_si_h3, pos_si_zeta3); + + if (HiggsRotFixedNeutral.determinant() < 0) + // condition (3) above, det(R) < 0 + { + // if negative, flip sign of h2 + HiggsRotFixed.row(pos_si_h2) *= -1; + } + + // Undo Goldstone rotation and flip that made G0 the first element + MatrixXd HiggsRotFixedGoldstone(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixedGoldstone.row(i) = HiggsRotFixed.row(i); + } + + HiggsRotFixed = MoveGoldstoneFirst*HiggsRotFixedGoldstone*RotGoldstone; + + // Extract the fixed mixing angles + double sina2 = HiggsRotFixedGoldstone(pos_si_h1, pos_si_zeta3); // +sin(a2) + double cosa2 = std::sqrt(1.0 - sina2*sina2); + alpha1 = std::asin(HiggsRotFixedGoldstone(pos_si_h1, pos_si_zeta2)/cosa2); // +sin(a1) cos(a2) + alpha2 = std::asin(sina2); + alpha3 = std::asin(HiggsRotFixedGoldstone(pos_si_h2, pos_si_zeta3)/cosa2); // +cos(a2) sin(a3) + + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) { - double tmp = posN[1]; - posN[1] = posN[2]; - posN[2] = tmp; + HiggsRotationMatrixEnsuredConvention[i][j] = HiggsRotFixed(i, j); } - if (MSM > MhDown) + } + + return; +} + + +void Class_Potential_C2HDM::TripleHiggsCouplings() +{ + if (!SetCurvatureDone) SetCurvatureArrays(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + + MatrixXd HiggsRot(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) { - double tmp = posN[0]; - posN[0] = posN[1]; - posN[1] = tmp; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } std::vector HiggsOrder(NHiggs); - HiggsOrder[0] = posG1; - HiggsOrder[1] = posG2; - HiggsOrder[2] = posMHCS1; - HiggsOrder[3] = posMHCS2; - HiggsOrder[4] = posG0; - for (int i = 5; i < 8; i++) - HiggsOrder[i] = posN[i - 5]; + HiggsOrder[0] = pos_G1; + HiggsOrder[1] = pos_G2; + HiggsOrder[2] = pos_H1; + HiggsOrder[3] = pos_H2; + HiggsOrder[4] = pos_G0; + if (UseHsmNotationInTripleHiggs) + { + HiggsOrder[5] = pos_h_SM; + HiggsOrder[6] = pos_h_l; + HiggsOrder[7] = pos_h_H; + } + else + { + HiggsOrder[5] = pos_h1; + HiggsOrder[6] = pos_h2; + HiggsOrder[7] = pos_h3; + } + MatrixXd HiggsRotSort(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); } - PosSM = 5; + std::vector TripleDeriv; + TripleDeriv = WeinbergThirdDerivative(); + std::vector>> GaugeBasis( + NHiggs, + std::vector>(NHiggs, std::vector(NHiggs))); + // double GaugeBasis[NHiggs][NHiggs][NHiggs]; + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + for (std::size_t k = 0; k < NHiggs; k++) + { + GaugeBasis[i][j][k] = + TripleDeriv.at(i + j * NHiggs + k * NHiggs * NHiggs); + } + } + } TripleHiggsCorrectionsCWPhysical.resize(NHiggs); TripleHiggsCorrectionsTreePhysical.resize(NHiggs); diff --git a/src/models/ClassPotentialCPintheDark.cpp b/src/models/ClassPotentialCPintheDark.cpp index d13d99b4c..f3eb24d07 100644 --- a/src/models/ClassPotentialCPintheDark.cpp +++ b/src/models/ClassPotentialCPintheDark.cpp @@ -1175,71 +1175,19 @@ void Class_Potential_CPintheDark::write() const ss << "The scale is given by mu = " << scale << " GeV \n"; - const double ZeroThreshold = 1e-5; - - std::size_t posGp = 0; - std::size_t posGm = 0; - std::size_t posHp = 0; - std::size_t posHm = 0; - std::size_t posHSM = 0; - std::size_t posG0 = 0; - std::size_t posh1 = 0; - std::size_t posh2 = 0; - std::size_t posh3 = 0; - - for (size_t i = 0; i < NHiggs; i++) - { - // the rotation matrix is diagonal besides for the neutral dark scalars - if (std::abs(HiggsRotationMatrix[i][0]) > ZeroThreshold) - posGp = i; - else if (std::abs(HiggsRotationMatrix[i][1]) > ZeroThreshold) - posGm = i; - else if (std::abs(HiggsRotationMatrix[i][2]) > ZeroThreshold) - posHp = i; - else if (std::abs(HiggsRotationMatrix[i][3]) > ZeroThreshold) - posHm = i; - else if (std::abs(HiggsRotationMatrix[i][4]) > ZeroThreshold) - posHSM = i; - else if (std::abs(HiggsRotationMatrix[i][5]) > ZeroThreshold) - posG0 = i; - - // the neutral dark scalars mix - if ((std::abs(HiggsRotationMatrix[i][6]) + - std::abs(HiggsRotationMatrix[i][7]) + - std::abs(HiggsRotationMatrix[i][8])) > ZeroThreshold) - { - // use that scalars are sorted by mass - if (posh1 == 0) - { - posh1 = i; - } - else - { - if (posh2 == 0) - { - posh2 = i; - } - else - { - posh3 = i; - } - } - } - } - std::vector HiggsMasses; HiggsMasses = HiggsMassesSquared(vevTree, 0); ss << "The mass spectrum is given by :\n"; - ss << "m_{G^+} = " << std::sqrt(HiggsMasses[posGp]) << " GeV \n" - << "m_{G^-} = " << std::sqrt(HiggsMasses[posGm]) << " GeV \n" - << "m_{H^+} = " << std::sqrt(HiggsMasses[posHp]) << " GeV \n" - << "m_{H^-} = " << std::sqrt(HiggsMasses[posHm]) << " GeV \n" - << "m_{hSM} = " << std::sqrt(HiggsMasses[posHSM]) << " GeV \n" - << "m_{G^0} = " << std::sqrt(HiggsMasses[posG0]) << " GeV \n" - << "m_{h_1} = " << std::sqrt(HiggsMasses[posh1]) << " GeV \n" - << "m_{h_2} = " << std::sqrt(HiggsMasses[posh2]) << " GeV \n" - << "m_{h_3} = " << std::sqrt(HiggsMasses[posh3]) << " GeV \n"; + ss << "m_{G^+} = " << std::sqrt(HiggsMasses[pos_Gp]) << " GeV \n" + << "m_{G^-} = " << std::sqrt(HiggsMasses[pos_Gm]) << " GeV \n" + << "m_{H^+} = " << std::sqrt(HiggsMasses[pos_Hp]) << " GeV \n" + << "m_{H^-} = " << std::sqrt(HiggsMasses[pos_Hm]) << " GeV \n" + << "m_{hSM} = " << std::sqrt(HiggsMasses[pos_HSM]) << " GeV \n" + << "m_{G^0} = " << std::sqrt(HiggsMasses[pos_G0]) << " GeV \n" + << "m_{h_1} = " << std::sqrt(HiggsMasses[pos_h1]) << " GeV \n" + << "m_{h_2} = " << std::sqrt(HiggsMasses[pos_h2]) << " GeV \n" + << "m_{h_3} = " << std::sqrt(HiggsMasses[pos_h3]) << " GeV \n"; Logger::Write(LoggingLevel::Default, ss.str()); } @@ -1258,7 +1206,7 @@ std::vector Class_Potential_CPintheDark::calc_CT() const retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; @@ -1314,22 +1262,21 @@ std::vector Class_Potential_CPintheDark::calc_CT() const return parCT; } -// mass basis triple couplings -void Class_Potential_CPintheDark::TripleHiggsCouplings() +/** + * Ensures the correct rotation matrix convention + */ +void Class_Potential_CPintheDark::AdjustRotationMatrix() { + const double ZeroThreshold = 1e-5; + if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); - - // position indices store the position of the physical fields - std::size_t posGp = 0; - std::size_t posGm = 0; - std::size_t posHp = 0; - std::size_t posHm = 0; - std::size_t posHSM = 0; - std::size_t posG0 = 0; - std::size_t posh1 = 0; - std::size_t posh2 = 0; - std::size_t posh3 = 0; + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (!CheckRotationMatrix()) // Check whether generically generated rotation + // matrix is proper rotation matrix + { + throw std::runtime_error("Error in rotation matrix."); + } MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) @@ -1340,59 +1287,221 @@ void Class_Potential_CPintheDark::TripleHiggsCouplings() } } - const double ZeroThreshold = 1e-5; + // CP in the Dark interaction basis + // 0 1 2 3 4 5 6 7 8 + // rho1, eta1, rho2, eta2, zeta1, psi1, zeta2, psi2, rhoS + int pos_rho1 = 0, pos_eta1 = 1, pos_rho2 = 2, pos_eta2 = 3, + pos_zeta1 = 4, pos_psi1 = 5, pos_zeta2 = 6, pos_psi2 = 7, + pos_rhoS = 8; + + // Indices of mass eigenstates for rotation from interaction to mass basis + pos_Gp = -1, pos_Gm = -1, pos_Hp = -1, pos_Hm = -1, pos_HSM = -1, + pos_G0 = -1, pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; - for (size_t i = 0; i < NHiggs; i++) - { - // the rotation matrix is diagonal besides for the neutral dark scalars - if (std::abs(HiggsRot(i, 0)) > ZeroThreshold) - posGp = i; - else if (std::abs(HiggsRot(i, 1)) > ZeroThreshold) - posGm = i; - else if (std::abs(HiggsRot(i, 2)) > ZeroThreshold) - posHp = i; - else if (std::abs(HiggsRot(i, 3)) > ZeroThreshold) - posHm = i; - else if (std::abs(HiggsRot(i, 4)) > ZeroThreshold) - posHSM = i; - else if (std::abs(HiggsRot(i, 5)) > ZeroThreshold) - posG0 = i; + // basis = {rho1, eta1, rho2, eta2, zeta1, psi1, zeta2, psi2, rhoS} + // the rotation matrix is diagonal besides for the neutral dark scalars + for (std::size_t i = 0; i < NHiggs; i++) // mass base index i corresponds to mass vector sorted in ascending mass + { + if (std::abs(HiggsRot(i, pos_rho1)) > ZeroThreshold) + { + pos_Gp = i; + } + else if (std::abs(HiggsRot(i, pos_eta1)) > ZeroThreshold) + { + pos_Gm = i; + } + else if (std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) + { + pos_Hp = i; + } + else if (std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_Hm = i; + } + else if (std::abs(HiggsRot(i, pos_zeta1)) > ZeroThreshold) + { + pos_HSM = i; + } + else if (std::abs(HiggsRot(i, pos_psi1)) > ZeroThreshold) + { + pos_G0 = i; + } // the neutral dark scalars mix - if ((std::abs(HiggsRot(i, 6)) + std::abs(HiggsRot(i, 7)) + - std::abs(HiggsRot(i, 8))) > ZeroThreshold) + else if ((std::abs(HiggsRot(i, pos_zeta2)) + std::abs(HiggsRot(i, pos_psi2)) + + std::abs(HiggsRot(i, pos_rhoS))) > ZeroThreshold) { // use that scalars are sorted by mass - if (posh1 == 0) + if (pos_h1 == -1) + { + pos_h1 = i; + } + else if (pos_h2 == -1) { - posh1 = i; + pos_h2 = i; } else { - if (posh2 == 0) - { - posh2 = i; - } - else - { - posh3 = i; - } + pos_h3 = i; } } } - // new rotation matrix with - MatrixXd HiggsRotSort(NHiggs, NHiggs); + // check if all position indices are set + if (pos_Gp == -1 or pos_Gm == -1 or pos_Hp == -1 or pos_Hm == -1 or + pos_HSM == -1 or pos_G0 == -1 or + pos_h1 == -1 or pos_h2 == -1 or pos_h3 == -1) + { + throw std::runtime_error("Error. Not all position indices are set."); + } + + // check if all other elements of rotation matrix are zero + bool zero_element = false; + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + int ii = int(i); + int jj = int(j); + if (not((jj == pos_rho1 and ii == pos_Gp) or + (jj == pos_eta1 and ii == pos_Gm) or + (jj == pos_rho2 and ii == pos_Hp) or + (jj == pos_eta2 and ii == pos_Hm) or + (jj == pos_zeta1 and ii == pos_HSM) or + (jj == pos_psi1 and ii == pos_G0) or + (jj == pos_zeta2 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_psi2 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_rhoS and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)))) + { + zero_element = true; + } + + if (zero_element and std::abs(HiggsRot(i, j)) > ZeroThreshold) + { + throw std::runtime_error("Error. Invalid rotation matrix detected."); + } + zero_element = false; + } + } + + MatrixXd HiggsRotFixed(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixed.row(i) = HiggsRot.row(i); + } + + if (HiggsRotFixed(pos_Gp, pos_rho1) < 0) // Gp rho1 (= +1) + { + HiggsRotFixed.row(pos_Gp) *= -1; + } + if (HiggsRotFixed(pos_Gm, pos_eta1) < 0) // Gm eta1 (= +1) + { + HiggsRotFixed.row(pos_Gm) *= -1; + } + if (HiggsRotFixed(pos_Hp, pos_rho2) < 0) // Hp rho2 (= +1) + { + HiggsRotFixed.row(pos_Hp) *= -1; + } + if (HiggsRotFixed(pos_Hm, pos_eta2) < 0) // Hm eta2 (= +1) + { + HiggsRotFixed.row(pos_Hm) *= -1; + } + if (HiggsRotFixed(pos_HSM, pos_zeta1) < 0) // HSM zeta1 (= +1) + { + HiggsRotFixed.row(pos_HSM) *= -1; + } + if (HiggsRotFixed(pos_G0, pos_psi1) < 0) // G0 psi1 (= +1) + { + HiggsRotFixed.row(pos_G0) *= -1; + } + + // Check dark neutral submatrix + // Use the "ScannerS" criteria from arXiv:2007.02985 Eq. (6) + // (since ScannerS uses the same parametrisation of the angles as BSMPT): + // * (1) if R[1][1] < 0: h1 -> -h1 (i.e. multiply the h1 row with -1) + // * (2) if R[3][3] < 0: h3 -> -h3 (i.e. multiply the h3 row with -1) + // * (3) if det R < 0: h2 -> -h2 (i.e. multiply the h2 row with -1) + + if (HiggsRotFixed(pos_h1, pos_zeta2) < 0) + // h1 zeta2 (condition (1) above, R11 < 0) + { + // if negative, flip sign of h1 + HiggsRotFixed.row(pos_h1) *= -1; + } + + if (HiggsRotFixed(pos_h3, pos_rhoS) < 0) + // h3 rhoS (condition (2) above, R33 < 0) + { + // if negative, flip sign of h3 + HiggsRotFixed.row(pos_h3) *= -1; + } + + // Calculate the determinant AFTER flipping the signs for rows 1 and 3 above + MatrixXd HiggsRotFixedNeutral(3, 3); + HiggsRotFixedNeutral(0, 0) = HiggsRotFixed(pos_h1, pos_zeta2); + HiggsRotFixedNeutral(0, 1) = HiggsRotFixed(pos_h1, pos_psi2); + HiggsRotFixedNeutral(0, 2) = HiggsRotFixed(pos_h1, pos_rhoS); + + HiggsRotFixedNeutral(1, 0) = HiggsRotFixed(pos_h2, pos_zeta2); + HiggsRotFixedNeutral(1, 1) = HiggsRotFixed(pos_h2, pos_psi2); + HiggsRotFixedNeutral(1, 2) = HiggsRotFixed(pos_h2, pos_rhoS); + + HiggsRotFixedNeutral(2, 0) = HiggsRotFixed(pos_h3, pos_zeta2); + HiggsRotFixedNeutral(2, 1) = HiggsRotFixed(pos_h3, pos_psi2); + HiggsRotFixedNeutral(2, 2) = HiggsRotFixed(pos_h3, pos_rhoS); + + if (HiggsRotFixedNeutral.determinant() < 0) + // condition (3) above, det(R) < 0 + { + // if negative, flip sign of h2 + HiggsRotFixed.row(pos_h2) *= -1; + } + + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRotationMatrixEnsuredConvention[i][j] = HiggsRotFixed(i, j); + } + } + + return; +} - HiggsRotSort.row(0) = HiggsRot.row(posGp); - HiggsRotSort.row(1) = HiggsRot.row(posGm); - HiggsRotSort.row(2) = HiggsRot.row(posHp); - HiggsRotSort.row(3) = HiggsRot.row(posHm); - HiggsRotSort.row(4) = HiggsRot.row(posHSM); - HiggsRotSort.row(5) = HiggsRot.row(posG0); - HiggsRotSort.row(6) = HiggsRot.row(posh1); - HiggsRotSort.row(7) = HiggsRot.row(posh2); - HiggsRotSort.row(8) = HiggsRot.row(posh3); +// mass basis triple couplings +void Class_Potential_CPintheDark::TripleHiggsCouplings() +{ + if (!SetCurvatureDone) SetCurvatureArrays(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + + MatrixXd HiggsRot(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; + } + } + + std::vector HiggsOrder(NHiggs); + HiggsOrder[0] = pos_Gp; + HiggsOrder[1] = pos_Gm; + HiggsOrder[2] = pos_Hp; + HiggsOrder[3] = pos_Hm; + HiggsOrder[4] = pos_HSM; + HiggsOrder[5] = pos_G0; + HiggsOrder[6] = pos_h1; + HiggsOrder[7] = pos_h2; + HiggsOrder[8] = pos_h3; + + MatrixXd HiggsRotSort(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); + } std::vector TripleDeriv; TripleDeriv = WeinbergThirdDerivative(); diff --git a/src/models/ClassPotentialCxSM.cpp b/src/models/ClassPotentialCxSM.cpp index c5f5baa40..fa907cb64 100644 --- a/src/models/ClassPotentialCxSM.cpp +++ b/src/models/ClassPotentialCxSM.cpp @@ -453,84 +453,27 @@ void Class_CxSM::write() const { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } - int posN[3]; - posN[0] = 3; - posN[1] = 4; - posN[2] = 5; - int posGCharged = 0, posG0 = 0; - double testsum = 0; - const double ZeroThreshold = 1e-5; - for (int i = 0; i < 3; i++) - { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 1)); - if (testsum > ZeroThreshold and posGCharged == 0) - { - posGCharged = i; - } - testsum = std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posG0 = i; - } - std::vector HiggsMasses; HiggsMasses = HiggsMassesSquared(vevTree, 0); - double MhUp = 0, MhDown = 0, MSM = 0; - double NeutralHiggs[3]; - for (int i = 0; i < 3; i++) - { - NeutralHiggs[i] = HiggsMasses[posN[i]]; - } - for (int i = 0; i < 3; i++) - { - if (std::sqrt(NeutralHiggs[i]) < 126 and std::sqrt(NeutralHiggs[i]) > 124) - MSM = std::sqrt(NeutralHiggs[i]); - } - if (std::sqrt(NeutralHiggs[0]) == MSM) - { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[1]); - } - else if (std::sqrt(NeutralHiggs[1]) == MSM) - { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[0]); - } - else - { - MhUp = std::sqrt(NeutralHiggs[1]); - MhDown = std::sqrt(NeutralHiggs[0]); - } - - if (MSM > MhUp) - { - double tmp = posN[1]; - posN[1] = posN[2]; - posN[2] = tmp; - } - if (MSM > MhDown) - { - double tmp = posN[0]; - posN[0] = posN[1]; - posN[1] = tmp; - } - MatrixXd NeutralMatrix(3, 3); for (int j = 0; j < 3; j++) { - for (int i = 0; i < 3; i++) - NeutralMatrix(i, j) = HiggsRot(posN[i], j + 3); + NeutralMatrix(0, j) = HiggsRot(pos_h_SM, j + 3); + NeutralMatrix(1, j) = HiggsRot(pos_h_l, j + 3); + NeutralMatrix(2, j) = HiggsRot(pos_h_H, j + 3); } ss << "The mass spectrum is given by :\n"; - ss << "m_{G^+}^2 = " << HiggsMasses[posGCharged] << " GeV^2 \n" - << "m_{G^0}^2 = " << HiggsMasses[posG0] << " GeV^2 \n" - << "m_{H_SM} = " << MSM << " GeV \n" - << "m_{H_l} = " << MhDown << " GeV \n" - << "m_{H_h} = " << MhUp << " GeV \n"; + ss << "m_{G^+}^2 = " << HiggsMasses[pos_Gp] << " GeV^2 \n" + << "m_{G^0}^2 = " << HiggsMasses[pos_G0] << " GeV^2 \n" + << "m_{H_SM} = " << std::sqrt(HiggsMasses[pos_h_SM]) << " GeV \n" + << "m_{H_l} = " << std::sqrt(HiggsMasses[pos_h_l]) << " GeV \n" + << "m_{H_h} = " << std::sqrt(HiggsMasses[pos_h_H]) << " GeV \n"; ss << "The neutral mixing Matrix is given by :\n"; bool IsNegative = NeutralMatrix(0, 1) < 0; ss << "H_{SM} = " << NeutralMatrix(0, 0) << " h "; @@ -589,7 +532,7 @@ std::vector Class_CxSM::calc_CT() const retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; @@ -765,15 +708,21 @@ std::vector Class_CxSM::calc_CT() const return parCT; } -void Class_CxSM::TripleHiggsCouplings() +/** + * Ensures the correct rotation matrix convention + */ +void Class_CxSM::AdjustRotationMatrix() { + const double ZeroThreshold = 1e-5; + if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); - std::vector HiggsOrder(NHiggs); - // Here you have to set the vector HiggsOrder. By telling e.g. HiggsOrder[0] = - // 5 you always want your 6th lightest particle to be the first particle in - // the vector (which has the index 5 because they are sorted by mass) + if (!CheckRotationMatrix()) // Check whether generically generated rotation + // matrix is proper rotation matrix + { + throw std::runtime_error("Error in rotation matrix."); + } MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) @@ -784,50 +733,238 @@ void Class_CxSM::TripleHiggsCouplings() } } - std::size_t posGp = 0, posGm = 0, posG0 = 0; - std::size_t posH1 = 0, posH2 = 0, posH3 = 0; - const double ZeroThreshold = 1e-5; + // CxSM interaction basis + // 0 1 2 3 4 5 + // Gp, Gm, G0, zeta1, zeta2, zeta3 + int pos_i_G1 = 0, pos_i_G2 = 1, pos_i_G0 = 2; + int pos_zeta1 = 3, pos_zeta2 = 4, pos_zeta3 = 5; - for (std::size_t i = 0; i < NHiggs; i++) + // Indices of mass eigenstates for rotation from interaction to mass basis + // Only neutral/physical part necessary, as Goldstones are already diagonal + pos_Gp = -1, pos_Gm = -1, pos_G0 = -1, + pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; + + // First finding the correct rows for the Goldstone bosons + for (std::size_t i = 0; i < 3; i++) { - // the rotation matrix is diagonal besides for the neutral scalars - if (std::abs(HiggsRot(i, 0)) > ZeroThreshold) - posGp = i; - else if (std::abs(HiggsRot(i, 1)) > ZeroThreshold) - posGm = i; - else if (std::abs(HiggsRot(i, 2)) > ZeroThreshold) - posG0 = i; - - // the neutral scalars mix - if ((std::abs(HiggsRot(i, 3)) + std::abs(HiggsRot(i, 4)) + - std::abs(HiggsRot(i, 5))) > ZeroThreshold) + if (std::abs(HiggsRot(i, pos_i_G1)) + + std::abs(HiggsRot(i, pos_i_G2)) > ZeroThreshold) { - // use that scalars are sorted by mass - if (posH1 == 0) + if (pos_Gp == -1) { - posH1 = i; + pos_Gp = i; } else { - if (posH2 == 0) - { - posH2 = i; - } - else - { - posH3 = i; - } + pos_Gm = i; + } + } + else if (std::abs(HiggsRot(i, pos_i_G0)) > ZeroThreshold) + { + pos_G0 = i; + } + else + { + throw std::runtime_error("Error. Non-Goldstone in Goldstone submatrix." + " Particle with negative mass? Check your " + "parameter point."); + } + } + + + // Starting from 3 to NHiggs, i.e. fixing the physical Higgs bosons + for (std::size_t i = 3; i < NHiggs; i++) + // mass base index i corresponds to mass vector sorted in ascending mass + { + // Neutral submatrix + if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) + + std::abs(HiggsRot(i, pos_zeta3)) > ZeroThreshold) + // use that mh1 < mh2 < mh3 + { + if (pos_h1 == -1) { + pos_h1 = i; + } else if (pos_h2 == -1) { + pos_h2 = i; + } else { + pos_h3 = i; + } + } + } + + // check if all position indices are set + if (pos_h1 == -1 or pos_h2 == -1 or pos_h3 == -1) + { + throw std::runtime_error("Error. Not all position indices are set."); + } + + // check if all other elements of rotation matrix are zero + bool zero_element = false; + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + int ii = int(i); + int jj = int(j); + if (not((jj == pos_Gp and ii == pos_Gp) or + (jj == pos_Gm and ii == pos_Gm) or + (jj == pos_G0 and ii == pos_G0) or + (jj == pos_zeta1 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_zeta2 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_zeta3 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)))) + { + zero_element = true; + } + + if (zero_element and std::abs(HiggsRot(i, j)) > ZeroThreshold) + { + throw std::runtime_error("Error. Invalid rotation matrix detected."); } + zero_element = false; + } + } + + // Determine the additional indices for the SM-like + // and lighter/heavier Higgses + pos_h_SM = -1, pos_h_l = -1, pos_h_H = -1; + + std::vector HiggsMasses; + HiggsMasses = HiggsMassesSquared(vevTree, 0); + + std::cout << "HiggsMasses=" << std::endl; + std::cout << HiggsMasses << std::endl; + + // Due to the masses being ordered, we will always have + // HiggsMasses[pos_h1] <= HiggsMasses[pos_h2] <= HiggsMasses[pos_h3] + double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h1]) + - SMConstants.C_MassSMHiggs); + double diff2 = std::abs(std::sqrt(HiggsMasses[pos_h2]) + - SMConstants.C_MassSMHiggs); + double diff3 = std::abs(std::sqrt(HiggsMasses[pos_h3]) + - SMConstants.C_MassSMHiggs); + if (diff1 < diff2 and diff1 < diff3) + { + pos_h_SM = pos_h1; + pos_h_l = pos_h2; + pos_h_H = pos_h3; + } + else if (diff2 < diff1 and diff2 < diff3) + { + pos_h_l = pos_h1; + pos_h_SM = pos_h2; + pos_h_H = pos_h3; + } + else + { + pos_h_l = pos_h1; + pos_h_H = pos_h2; + pos_h_SM = pos_h3; + } + + MatrixXd HiggsRotFixed(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixed.row(i) = HiggsRot.row(i); + } + + // Goldstone elements (all diagonal) + if (HiggsRotFixed(pos_Gp, pos_Gp) < 0) // Gp Gp (+1) + { + HiggsRotFixed.row(pos_Gp) *= -1; + } + if (HiggsRotFixed(pos_Gm, pos_Gm) < 0) // Gm Gm (+1) + { + HiggsRotFixed.row(pos_Gm) *= -1; + } + if (HiggsRotFixed(pos_G0, pos_G0) < 0) // G0 G0 (+1) + { + HiggsRotFixed.row(pos_G0) *= -1; + } + + // Check neutral submatrix + // Use the "ScannerS" criteria from arXiv:2007.02985 Eq. (6) + // (since ScannerS uses the same parametrisation of the angles as BSMPT): + // * (1) if R[1][1] < 0: h1 -> -h1 (i.e. multiply the h1 row with -1) + // * (2) if R[3][3] < 0: h3 -> -h3 (i.e. multiply the h3 row with -1) + // * (3) if det R < 0: h2 -> -h2 (i.e. multiply the h2 row with -1) + + // check neutral, CP-even submatrix + if (HiggsRotFixed(pos_h1, pos_zeta1) < 0) + // h1 zeta1 (condition (1) above, R11 < 0) + { + // if negative, flip sign of h1 + HiggsRotFixed.row(pos_h1) *= -1; + } + + if (HiggsRotFixed(pos_h3, pos_zeta3) < 0) + // h3 zeta3 (condition (2) above, R33 < 0) + { + // if negative, flip sign of h3 + HiggsRotFixed.row(pos_h3) *= -1; + } + + // Calculate the determinant AFTER flipping the signs for rows 1 and 3 above + MatrixXd HiggsRotFixedNeutral(3, 3); + HiggsRotFixedNeutral(0, 0) = HiggsRotFixed(pos_h1, pos_zeta1); + HiggsRotFixedNeutral(0, 1) = HiggsRotFixed(pos_h1, pos_zeta2); + HiggsRotFixedNeutral(0, 2) = HiggsRotFixed(pos_h1, pos_zeta3); + + HiggsRotFixedNeutral(1, 0) = HiggsRotFixed(pos_h2, pos_zeta1); + HiggsRotFixedNeutral(1, 1) = HiggsRotFixed(pos_h2, pos_zeta2); + HiggsRotFixedNeutral(1, 2) = HiggsRotFixed(pos_h2, pos_zeta3); + + HiggsRotFixedNeutral(2, 0) = HiggsRotFixed(pos_h3, pos_zeta1); + HiggsRotFixedNeutral(2, 1) = HiggsRotFixed(pos_h3, pos_zeta2); + HiggsRotFixedNeutral(2, 2) = HiggsRotFixed(pos_h3, pos_zeta3); + + if (HiggsRotFixedNeutral.determinant() < 0) + // condition (3) above, det(R) < 0 + { + // if negative, flip sign of h2 + HiggsRotFixed.row(pos_h2) *= -1; + } + + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRotationMatrixEnsuredConvention[i][j] = HiggsRotFixed(i, j); + } + } + + return; +} + +void Class_CxSM::TripleHiggsCouplings() +{ + if (!SetCurvatureDone) SetCurvatureArrays(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + + MatrixXd HiggsRot(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } - // mass order: Gp, Gm, G0, H1, H2, H3 - HiggsOrder[0] = posGp; - HiggsOrder[1] = posGm; - HiggsOrder[2] = posG0; - HiggsOrder[3] = posH1; - HiggsOrder[4] = posH2; - HiggsOrder[5] = posH3; + std::vector HiggsOrder(NHiggs); + // mass order: Gp, Gm, G0, h1, h2, h3 + HiggsOrder[0] = pos_Gp; + HiggsOrder[1] = pos_Gm; + HiggsOrder[2] = pos_G0; + HiggsOrder[3] = pos_h1; + HiggsOrder[4] = pos_h2; + HiggsOrder[5] = pos_h3; + + MatrixXd HiggsRotSort(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); + } std::vector TripleDeriv; TripleDeriv = WeinbergThirdDerivative(); @@ -846,13 +983,6 @@ void Class_CxSM::TripleHiggsCouplings() } } - MatrixXd HiggsRotSort(NHiggs, NHiggs); - - for (std::size_t i = 0; i < NHiggs; i++) - { - HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); - } - TripleHiggsCorrectionsCWPhysical.resize(NHiggs); TripleHiggsCorrectionsTreePhysical.resize(NHiggs); TripleHiggsCorrectionsCTPhysical.resize(NHiggs); @@ -1305,7 +1435,7 @@ void Class_CxSM::Debugging(const std::vector &input, retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; diff --git a/src/models/ClassPotentialN2HDM.cpp b/src/models/ClassPotentialN2HDM.cpp index 2caed38ef..835dea12d 100644 --- a/src/models/ClassPotentialN2HDM.cpp +++ b/src/models/ClassPotentialN2HDM.cpp @@ -485,90 +485,27 @@ void Class_Potential_N2HDM::write() const { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; - } - } - - MatrixXd HiggsRotSort(NHiggs, NHiggs); - int posMHCS1 = 0; - int posA = 0; - int posN[3]; - int countposN = 0; - double testsum = 0; - const double ZeroThreshold = 1e-5; - - for (std::size_t i = 3; i < NHiggs; i++) - { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 1)); - if (testsum > ZeroThreshold) posMHCS1 = i; - testsum = std::abs(HiggsRot(i, 6)) + std::abs(HiggsRot(i, 7)) + - std::abs(HiggsRot(i, 8)); - if (testsum > ZeroThreshold) - { - posN[countposN] = i; - countposN++; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } - testsum = std::abs(HiggsRot(i, 4)) + std::abs(HiggsRot(i, 5)); - if (testsum > ZeroThreshold) posA = i; } std::vector HiggsMasses; HiggsMasses = HiggsMassesSquared(vevTree, 0); - double NeutralHiggs[3]; - double MSMlocal = 0, MhUplocal = 0, MhDownlocal = 0; - for (int i = 0; i < 3; i++) - { - NeutralHiggs[i] = HiggsMasses[posN[i]]; - } - for (int i = 0; i < 3; i++) - { - if (std::sqrt(NeutralHiggs[i]) < 126 and std::sqrt(NeutralHiggs[i]) > 124) - MSMlocal = std::sqrt(NeutralHiggs[i]); - } - if (std::sqrt(NeutralHiggs[0]) == MSMlocal) - { - MhUplocal = std::sqrt(NeutralHiggs[2]); - MhDownlocal = std::sqrt(NeutralHiggs[1]); - } - else if (std::sqrt(NeutralHiggs[1]) == MSMlocal) - { - MhUplocal = std::sqrt(NeutralHiggs[2]); - MhDownlocal = std::sqrt(NeutralHiggs[0]); - } - else - { - MhUplocal = std::sqrt(NeutralHiggs[1]); - MhDownlocal = std::sqrt(NeutralHiggs[0]); - } - - if (MSMlocal > MhUplocal) - { - double tmp = posN[1]; - posN[1] = posN[2]; - posN[2] = tmp; - } - if (MSMlocal > MhDownlocal) - { - double tmp = posN[0]; - posN[0] = posN[1]; - posN[1] = tmp; - } - MatrixXd NeutralMixing(3, 3); for (int i = 0; i < 3; i++) { - NeutralMixing(0, i) = HiggsRot(posN[0], i + 6); - NeutralMixing(1, i) = HiggsRot(posN[1], i + 6); - NeutralMixing(2, i) = HiggsRot(posN[2], i + 6); + NeutralMixing(0, i) = HiggsRot(pos_h_SM, i + 6); + NeutralMixing(1, i) = HiggsRot(pos_h_l, i + 6); + NeutralMixing(2, i) = HiggsRot(pos_h_H, i + 6); } ss << "The mass spectrum is given by :\n"; - ss << "m_{H^+} = " << std::sqrt(HiggsMasses[posMHCS1]) << " GeV \n" - << "m_A = " << std::sqrt(HiggsMasses[posA]) << " GeV \n" - << "m_{H_SM} = " << MSMlocal << " GeV \n" - << "m_{H_l} = " << MhDownlocal << " GeV \n" - << "m_{H_h} = " << MhUplocal << " GeV \n"; + ss << "m_{H^+} = " << std::sqrt(HiggsMasses[pos_H1]) << " GeV \n" + << "m_A = " << std::sqrt(HiggsMasses[pos_A]) << " GeV \n" + << "m_{H_SM} = " << std::sqrt(HiggsMasses[pos_h_SM]) << " GeV \n" + << "m_{H_l} = " << std::sqrt(HiggsMasses[pos_h_l]) << " GeV \n" + << "m_{H_h} = " << std::sqrt(HiggsMasses[pos_h_H]) << " GeV \n"; ss << "The neutral mixing Matrix is given by :\n"; bool IsNegative = NeutralMixing(0, 1) < 0; @@ -623,7 +560,7 @@ std::vector Class_Potential_N2HDM::calc_CT() const retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; @@ -772,26 +709,20 @@ std::vector Class_Potential_N2HDM::calc_CT() const return parCT; } -void Class_Potential_N2HDM::TripleHiggsCouplings() +/** + * Ensures the correct rotation matrix convention + */ +void Class_Potential_N2HDM::AdjustRotationMatrix() { + const double ZeroThreshold = 1e-5; + if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); - std::vector TripleDeriv; - TripleDeriv = WeinbergThirdDerivative(); - std::vector>> GaugeBasis( - NHiggs, - std::vector>(NHiggs, std::vector(NHiggs))); - for (std::size_t i = 0; i < NHiggs; i++) + if (!CheckRotationMatrix()) // Check whether generically generated rotation + // matrix is proper rotation matrix { - for (std::size_t j = 0; j < NHiggs; j++) - { - for (std::size_t k = 0; k < NHiggs; k++) - { - GaugeBasis[i][j][k] = - TripleDeriv.at(i + j * NHiggs + k * NHiggs * NHiggs); - } - } + throw std::runtime_error("Error in rotation matrix."); } MatrixXd HiggsRot(NHiggs, NHiggs); @@ -803,93 +734,301 @@ void Class_Potential_N2HDM::TripleHiggsCouplings() } } - MatrixXd HiggsRotSort(NHiggs, NHiggs); - int posMHCS1 = 0, posMHCS2 = 0; - int posA = 0; - std::vector posN(3); - std::size_t countposN = 0; - std::size_t posG1 = 0, posG2 = 0, posG0 = 0; + // N2HDM interaction basis + // 0 1 2 3 4 5 6 7 8 + // rho1, rho2, eta1, eta2, psi1, psi2, zeta1, zeta2, rhoS + int pos_rho1 = 0, pos_rho2 = 1, pos_eta1 = 2, pos_eta2 = 3, + pos_psi1 = 4, pos_psi2 = 5, pos_zeta1 = 6, pos_zeta2 = 7, + pos_rhoS = 8; - double testsum = 0; - const double ZeroThreshold = 1e-5; + // Unlike in the C2HDM, there is no CP-mixing, so we do not have to rotate out the Goldstone - for (int i = 0; i < 3; i++) + // Indices of mass eigenstates for rotation from interaction to mass basis + pos_G0 = -1, pos_G1 = -1, pos_G2 = -1, pos_H1 = -1, pos_H2 = -1, + pos_h1 = -1, pos_h2 = -1, pos_h3 = -1, pos_A = -1; + + // Going from 0 to 2, i.e. fixing the Goldstone indices first + // (using that the Goldstone masses appear before the physical Higgs bosons + // since they are the smallest mass eigenvalues (= 0 in the Landau gauge)) + for (std::size_t i = 0; i < 3; i++) + // mass base index i corresponds to mass vector sorted in ascending mass { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 1)); - if (testsum > ZeroThreshold) posG1 = i; - testsum = std::abs(HiggsRot(i, 2)) + std::abs(HiggsRot(i, 3)); - if (testsum > ZeroThreshold) posG2 = i; - testsum = std::abs(HiggsRot(i, 4)) + std::abs(HiggsRot(i, 5)); - if (testsum > ZeroThreshold) posG0 = i; + // Charged submatrix + if (std::abs(HiggsRot(i, pos_rho1)) + + std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) + { + pos_G1 = i; + } + else if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_G2 = i; + } + // Neutral CP-odd submatrix + else if (std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) + { + pos_G0 = i; + } + else + { + throw std::runtime_error("Error. Non-Goldstone in Goldstone submatrix." + " Particle with negative mass? Check your " + "parameter point."); + } } + + // Now going from 3 to NHiggs, i.e. fixing the physical Higgs bosons for (std::size_t i = 3; i < NHiggs; i++) + // mass base index i corresponds to mass vector sorted in ascending mass + { + // Charged submatrix + if (std::abs(HiggsRot(i, pos_rho1)) + + std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) + { + pos_H1 = i; + } + else if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_H2 = i; + } + // Neutral CP-odd submatrix + else if (std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) + { + pos_A = i; + } + // Neutral CP-even submatrix + else if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) + + std::abs(HiggsRot(i, pos_rhoS)) > ZeroThreshold) + // use that mh1 < mh2 < mh3 + { + if (pos_h1 == -1) { + pos_h1 = i; + } else if (pos_h2 == -1) { + pos_h2 = i; + } else { + pos_h3 = i; + } + } + } + + // check if all position indices are set + if (pos_G0 == -1 or pos_G1 == -1 or pos_G2 == -1 or + pos_H1 == -1 or pos_H2 == -1 or + pos_h1 == -1 or pos_h2 == -1 or pos_h3 == -1 or pos_A == -1) + { + throw std::runtime_error("Error. Not all position indices are set."); + } + + // check if all other elements of rotation matrix are zero + bool zero_element = false; + for (std::size_t i = 0; i < NHiggs; i++) { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 1)); - if (testsum > ZeroThreshold) posMHCS1 = i; - testsum = std::abs(HiggsRot(i, 2)) + std::abs(HiggsRot(i, 3)); - if (testsum > ZeroThreshold) posMHCS2 = i; - testsum = std::abs(HiggsRot(i, 6)) + std::abs(HiggsRot(i, 7)) + - std::abs(HiggsRot(i, 8)); - if (testsum > ZeroThreshold) + for (std::size_t j = 0; j < NHiggs; j++) { - posN.at(countposN) = i; - countposN++; + int ii = int(i); + int jj = int(j); + if (not((jj == pos_rho1 and (ii == pos_G1 or ii == pos_H1)) or + (jj == pos_rho2 and (ii == pos_G1 or ii == pos_H1)) or + (jj == pos_eta1 and (ii == pos_G2 or ii == pos_H2)) or + (jj == pos_eta2 and (ii == pos_G2 or ii == pos_H2)) or + (jj == pos_psi1 and (ii == pos_G0 or ii == pos_A)) or + (jj == pos_psi2 and (ii == pos_G0 or ii == pos_A)) or + (jj == pos_zeta1 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_zeta2 and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)) or + (jj == pos_rhoS and (ii == pos_h1 or ii == pos_h2 or ii == pos_h3)))) + { + zero_element = true; + } + + if (zero_element and std::abs(HiggsRot(i, j)) > ZeroThreshold) + { + throw std::runtime_error("Error. Invalid rotation matrix detected."); + } + zero_element = false; } - testsum = std::abs(HiggsRot(i, 4)) + std::abs(HiggsRot(i, 5)); - if (testsum > ZeroThreshold) posA = i; } + // Determine the additional indices for the SM-like + // and lighter/heavier Higgses + pos_h_SM = -1, pos_h_l = -1, pos_h_H = -1; + std::vector HiggsMasses; HiggsMasses = HiggsMassesSquared(vevTree, 0); - double NeutralHiggs[3]; - for (int i = 0; i < 3; i++) + // Due to the masses being ordered, we will always have + // HiggsMasses[pos_h1] <= HiggsMasses[pos_h2] <= HiggsMasses[pos_h3] + double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h1]) + - SMConstants.C_MassSMHiggs); + double diff2 = std::abs(std::sqrt(HiggsMasses[pos_h2]) + - SMConstants.C_MassSMHiggs); + double diff3 = std::abs(std::sqrt(HiggsMasses[pos_h3]) + - SMConstants.C_MassSMHiggs); + if (diff1 < diff2 and diff1 < diff3) { - NeutralHiggs[i] = HiggsMasses[posN.at(i)]; + pos_h_SM = pos_h1; + pos_h_l = pos_h2; + pos_h_H = pos_h3; } - for (int i = 0; i < 3; i++) + else if (diff2 < diff1 and diff2 < diff3) + { + pos_h_l = pos_h1; + pos_h_SM = pos_h2; + pos_h_H = pos_h3; + } + else { - if (std::sqrt(NeutralHiggs[i]) < 126 and std::sqrt(NeutralHiggs[i]) > 124) - MSM = std::sqrt(NeutralHiggs[i]); + pos_h_l = pos_h1; + pos_h_H = pos_h2; + pos_h_SM = pos_h3; } - if (std::sqrt(NeutralHiggs[0]) == MSM) + + MatrixXd HiggsRotFixed(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[1]); + HiggsRotFixed.row(i) = HiggsRot.row(i); } - else if (std::sqrt(NeutralHiggs[1]) == MSM) + + // Neutral CP-odd submatrix + if (HiggsRotFixed(pos_G0, pos_psi1) < 0) // G0 psi1 (+ cos(beta)) { - MhUp = std::sqrt(NeutralHiggs[2]); - MhDown = std::sqrt(NeutralHiggs[0]); + HiggsRotFixed.row(pos_G0) *= -1; } - else + if (HiggsRotFixed(pos_A, pos_psi2) < 0) // A psi2 (+ cos(beta)) { - MhUp = std::sqrt(NeutralHiggs[1]); - MhDown = std::sqrt(NeutralHiggs[0]); + HiggsRotFixed.row(pos_A) *= -1; } - if (MSM > MhUp) + // charged submatrix + if (HiggsRotFixed(pos_G1, pos_rho1) < 0) // G1 rho1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_G1) *= -1; + } + if (HiggsRotFixed(pos_G2, pos_eta1) < 0) // G2 eta1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_G2) *= -1; + } + if (HiggsRotFixed(pos_H1, pos_rho2) < 0) // H1 rho2 (+ cos(beta)) { - double tmp = posN[1]; - posN[1] = posN[2]; - posN[2] = tmp; + HiggsRotFixed.row(pos_H1) *= -1; } - if (MSM > MhDown) + if (HiggsRotFixed(pos_H2, pos_eta2) < 0) // H2 eta2 (+ cos(beta)) { - double tmp = posN[0]; - posN[0] = posN[1]; - posN[1] = tmp; + HiggsRotFixed.row(pos_H2) *= -1; } - HiggsRotSort.row(0) = HiggsRot.row(posG1); - HiggsRotSort.row(1) = HiggsRot.row(posG2); - HiggsRotSort.row(2) = HiggsRot.row(posMHCS1); - HiggsRotSort.row(3) = HiggsRot.row(posMHCS2); - HiggsRotSort.row(4) = HiggsRot.row(posG0); - HiggsRotSort.row(5) = HiggsRot.row(posA); - for (int i = 6; i < 9; i++) + // Check neutral CP-even submatrix + // Use the "ScannerS" criteria from arXiv:2007.02985 Eq. (6) + // (since ScannerS uses the same parametrisation of the angles as BSMPT): + // * (1) if R[1][1] < 0: h1 -> -h1 (i.e. multiply the h1 row with -1) + // * (2) if R[3][3] < 0: h3 -> -h3 (i.e. multiply the h3 row with -1) + // * (3) if det R < 0: h2 -> -h2 (i.e. multiply the h2 row with -1) + + // check neutral, CP-even submatrix + if (HiggsRotFixed(pos_h1, pos_zeta1) < 0) + // h1 zeta1 (condition (1) above, R11 < 0) + { + // if negative, flip sign of h1 + HiggsRotFixed.row(pos_h1) *= -1; + } + + if (HiggsRotFixed(pos_h3, pos_rhoS) < 0) + // h3 rhoS (condition (2) above, R33 < 0) + { + // if negative, flip sign of h3 + HiggsRotFixed.row(pos_h3) *= -1; + } + + // Calculate the determinant AFTER flipping the signs for rows 1 and 3 above + MatrixXd HiggsRotFixedNeutral(3, 3); + HiggsRotFixedNeutral(0, 0) = HiggsRotFixed(pos_h1, pos_zeta1); + HiggsRotFixedNeutral(0, 1) = HiggsRotFixed(pos_h1, pos_zeta2); + HiggsRotFixedNeutral(0, 2) = HiggsRotFixed(pos_h1, pos_rhoS); + + HiggsRotFixedNeutral(1, 0) = HiggsRotFixed(pos_h2, pos_zeta1); + HiggsRotFixedNeutral(1, 1) = HiggsRotFixed(pos_h2, pos_zeta2); + HiggsRotFixedNeutral(1, 2) = HiggsRotFixed(pos_h2, pos_rhoS); + + HiggsRotFixedNeutral(2, 0) = HiggsRotFixed(pos_h3, pos_zeta1); + HiggsRotFixedNeutral(2, 1) = HiggsRotFixed(pos_h3, pos_zeta2); + HiggsRotFixedNeutral(2, 2) = HiggsRotFixed(pos_h3, pos_rhoS); + + if (HiggsRotFixedNeutral.determinant() < 0) + // condition (3) above, det(R) < 0 { - HiggsRotSort.row(i) = HiggsRot.row(posN[i - 6]); + // if negative, flip sign of h2 + HiggsRotFixed.row(pos_h2) *= -1; + } + + // Extract the fixed mixing angles + double sina2 = HiggsRotFixed(pos_h1, pos_rhoS); // +sin(a2) + double cosa2 = std::sqrt(1.0 - sina2*sina2); + alpha1 = std::asin(HiggsRotFixed(pos_h1, pos_zeta2)/cosa2); // +sin(a1) cos(a2) + alpha2 = std::asin(sina2); + alpha3 = std::asin(HiggsRotFixed(pos_h2, pos_rhoS)/cosa2); // +cos(a2) sin(a3) + + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRotationMatrixEnsuredConvention[i][j] = HiggsRotFixed(i, j); + } + } + + return; +} + +void Class_Potential_N2HDM::TripleHiggsCouplings() +{ + if (!SetCurvatureDone) SetCurvatureArrays(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + + MatrixXd HiggsRot(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; + } + } + + std::vector HiggsOrder(NHiggs); + HiggsOrder[0] = pos_G1; + HiggsOrder[1] = pos_G2; + HiggsOrder[2] = pos_H1; + HiggsOrder[3] = pos_H2; + HiggsOrder[4] = pos_G0; + HiggsOrder[5] = pos_A; + HiggsOrder[6] = pos_h_SM; + HiggsOrder[7] = pos_h_l; + HiggsOrder[8] = pos_h_H; + + MatrixXd HiggsRotSort(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); + } + + std::vector TripleDeriv; + TripleDeriv = WeinbergThirdDerivative(); + std::vector>> GaugeBasis( + NHiggs, + std::vector>(NHiggs, std::vector(NHiggs))); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + for (std::size_t k = 0; k < NHiggs; k++) + { + GaugeBasis[i][j][k] = + TripleDeriv.at(i + j * NHiggs + k * NHiggs * NHiggs); + } + } } TripleHiggsCorrectionsCWPhysical.resize(NHiggs); diff --git a/src/models/ClassPotentialOrigin.cpp b/src/models/ClassPotentialOrigin.cpp index 271e7d059..6d1ae7947 100644 --- a/src/models/ClassPotentialOrigin.cpp +++ b/src/models/ClassPotentialOrigin.cpp @@ -959,6 +959,48 @@ Class_Potential_Origin::SecondDerivativeOfEigenvaluesNonRepeated( return res; } +// Sanity check to make sure HiggsRotationMatrix is a proper rotation +// matrix, i.e. its inverse should correspond to its transpose, and its +// determinant should be +1 or -1 +bool Class_Potential_Origin::CheckRotationMatrix() +{ + MatrixXd mat(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + mat(i, j) = HiggsRotationMatrix[i][j]; + } + } + + double precision = 1e-10; + + bool AbsDetIsOne = almost_the_same(std::abs(mat.determinant()), 1., + precision); + bool InvEqTrans = true; + + auto inv = mat.inverse(); + auto transp = mat.transpose(); + + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + if (!almost_the_same(inv(i, j), transp(i, j), precision)) + { + InvEqTrans = false; + break; + } + } + } + + if (AbsDetIsOne and InvEqTrans) + { + return true; + } + return false; +} + void Class_Potential_Origin::CalculatePhysicalCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); @@ -1395,7 +1437,7 @@ void Class_Potential_Origin::CalculatePhysicalCouplings() } } - CalcCouplingsdone = true; + CalcCouplingsDone = true; return; } @@ -1403,7 +1445,7 @@ void Class_Potential_Origin::CalculatePhysicalCouplings() std::vector Class_Potential_Origin::WeinbergFirstDerivative() const { std::vector res; - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { // CalculatePhysicalCouplings(); std::string retmes = __func__; @@ -1497,7 +1539,7 @@ std::vector Class_Potential_Origin::WeinbergFirstDerivative() const Eigen::MatrixXd Class_Potential_Origin::WeinbergSecondDerivativeAsMatrixXd() const { - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { // CalculatePhysicalCouplings(); std::string retmes = __func__; @@ -1656,7 +1698,7 @@ std::vector Class_Potential_Origin::WeinbergSecondDerivative() const std::vector Class_Potential_Origin::WeinbergThirdDerivative() const { - if (not CalcCouplingsdone) + if (not CalcCouplingsDone) { std::string retmes = __func__; retmes += " tries to use Physical couplings but they are not initialised."; @@ -1902,7 +1944,7 @@ std::vector Class_Potential_Origin::WeinbergThirdDerivative() const std::vector Class_Potential_Origin::WeinbergForthDerivative() const { - if (not CalcCouplingsdone) + if (not CalcCouplingsDone) { std::string retmes = __func__; retmes += " tries to use Physical couplings but they are not initialised."; @@ -3387,6 +3429,9 @@ void Class_Potential_Origin::initVectors() vec3Complex{NQuarks, vec2Complex{NHiggs, vec1Complex(NHiggs, 0)}}}; HiggsVev = std::vector(NHiggs, 0); + + HiggsRotationMatrixEnsuredConvention = + std::vector>{NHiggs, std::vector(NHiggs, 0)}; } void Class_Potential_Origin::sym2Dim( @@ -3472,7 +3517,7 @@ void Class_Potential_Origin::sym4Dim( void Class_Potential_Origin::resetbools() { SetCurvatureDone = false; - CalcCouplingsdone = false; + CalcCouplingsDone = false; CalculatedTripleCopulings = false; parStored.clear(); parCTStored.clear(); @@ -3778,6 +3823,8 @@ Class_Potential_Origin::initModel(const std::vector &par) CalculateDebye(); CalculateDebyeGauge(); + AdjustRotationMatrix(); + parStored = par; parCTStored = parCT; diff --git a/src/models/ClassPotentialR2HDM.cpp b/src/models/ClassPotentialR2HDM.cpp index 12c40c4cf..285d5b6e8 100644 --- a/src/models/ClassPotentialR2HDM.cpp +++ b/src/models/ClassPotentialR2HDM.cpp @@ -563,68 +563,38 @@ void Class_Potential_R2HDM::write() const { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; - } - } - - int posMHCS1 = 0; - int posN[2]; - int countposN = 0; - int posA = 0; - int posG1 = 0, posG0 = 0; - double testsum = 0; - const double ZeroThreshold = 1e-5; - for (std::size_t i = 0; i < 3; i++) - { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posG1 = i; - // testsum = std::abs(HiggsRot(i,1)) + std::abs(HiggsRot(i,3)); - // if(testsum > ZeroThreshold) posG2 = i; - testsum = std::abs(HiggsRot(i, 5)) + std::abs(HiggsRot(i, 7)); - if (testsum > ZeroThreshold) posG0 = i; - } - for (std::size_t i = 3; i < NHiggs; i++) - { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posMHCS1 = i; - // testsum = std::abs(HiggsRot(i,1)) + std::abs(HiggsRot(i,3)); - // if(testsum > ZeroThreshold) posMHCS2 = i; - testsum = std::abs(HiggsRot(i, 4)) + std::abs(HiggsRot(i, 6)); - if (testsum > ZeroThreshold) - { - posN[countposN] = i; - countposN++; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } - testsum = std::abs(HiggsRot(i, 5)) + std::abs(HiggsRot(i, 7)); - if (testsum > ZeroThreshold) posA = i; } std::vector HiggsMasses; HiggsMasses = HiggsMassesSquared(vevTree, 0); - ss << "The mass spectrum is given by :\n"; - ss << "m_{G^+}^2 = " << HiggsMasses[posG1] << " GeV^2 \n"; - ss << "m_{G^0}^2 = " << HiggsMasses[posG0] << " GeV^2 \n"; - ss << "m_{H^+} = " << std::sqrt(HiggsMasses[posMHCS1]) << " GeV \n" - << "m_h = " << std::sqrt(HiggsMasses[posN[0]]) << " GeV \n" - << "m_H = " << std::sqrt(HiggsMasses[posN[1]]) << " GeV \n" - << "m_A = " << std::sqrt(HiggsMasses[posA]) << " GeV \n"; + ss << "The mass spectrum is given by :\n" + << "m_{G^+}^2 = " << HiggsMasses[pos_G1] << " GeV^2 \n" + << "m_{G^-}^2 = " << HiggsMasses[pos_G2] << " GeV^2 \n" + << "m_{G^0}^2 = " << HiggsMasses[pos_G0] << " GeV^2 \n" + << "m_{H^+} = " << std::sqrt(HiggsMasses[pos_H1]) << " GeV \n" + << "m_{H^-} = " << std::sqrt(HiggsMasses[pos_H2]) << " GeV \n" + << "m_h = " << std::sqrt(HiggsMasses[pos_h]) << " GeV \n" + << "m_H = " << std::sqrt(HiggsMasses[pos_H]) << " GeV \n" + << "m_A = " << std::sqrt(HiggsMasses[pos_A]) << " GeV \n"; ss << "The neutral mixing Matrix is given by :\n"; - ss << "h = " << HiggsRot(posN[0], 4) << " zeta_1 "; - bool IsNegative = HiggsRot(posN[0], 6) < 0; + ss << "h = " << HiggsRot(pos_h, 4) << " zeta_1 "; + bool IsNegative = HiggsRot(pos_h, 6) < 0; if (IsNegative) ss << "-"; else ss << "+"; - ss << std::abs(HiggsRot(posN[0], 6)) << " zeta_2\n" - << "H = " << HiggsRot(posN[1], 4) << " zeta_1 "; - IsNegative = HiggsRot(posN[1], 6) < 0; + ss << std::abs(HiggsRot(pos_h, 6)) << " zeta_2\n" + << "H = " << HiggsRot(pos_H, 4) << " zeta_1 "; + IsNegative = HiggsRot(pos_H, 6) < 0; if (IsNegative) ss << "-"; else ss << "+"; - ss << std::abs(HiggsRot(posN[1], 6)) << " zeta_2" << std::endl; + ss << std::abs(HiggsRot(pos_H, 6)) << " zeta_2" << std::endl; Logger::Write(LoggingLevel::Default, ss.str()); } @@ -642,7 +612,7 @@ std::vector Class_Potential_R2HDM::calc_CT() const retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; @@ -743,93 +713,270 @@ std::vector Class_Potential_R2HDM::calc_CT() const } /** - * Calculates the corrections to the Triple higgs couplings in the mass basis. - * - * Use the vector TripleHiggsCorrectionsCWPhysical to save your couplings and - * set the nTripleCouplings to the number of couplings you want as output. + * Ensures the correct rotation matrix convention */ -void Class_Potential_R2HDM::TripleHiggsCouplings() +void Class_Potential_R2HDM::AdjustRotationMatrix() { + const double ZeroThreshold = 1e-5; + if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); - std::vector TripleDeriv; - TripleDeriv = WeinbergThirdDerivative(); - std::vector>> GaugeBasis( - NHiggs, - std::vector>(NHiggs, std::vector(NHiggs))); + if (!CheckRotationMatrix()) // Check whether generically generated rotation + // matrix is proper rotation matrix + { + throw std::runtime_error("Error in rotation matrix."); + } + + MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { for (std::size_t j = 0; j < NHiggs; j++) { - for (std::size_t k = 0; k < NHiggs; k++) + HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + } + } + + // initialize position indices (new initialization for each point in multiline + // files) + pos_G1 = -1, pos_G2 = -1, pos_H1 = -1, pos_H2 = -1, pos_G0 = -1, pos_A = -1, + pos_H = -1, pos_h = -1; + + // interaction basis + // rho1, eta1, rho2, eta2, zeta1, psi1, zeta2, psi2 + int pos_rho1 = 0, pos_eta1 = 1, pos_rho2 = 2, pos_eta2 = 3, pos_zeta1 = 4, + pos_psi1 = 5, pos_zeta2 = 6, pos_psi2 = 7; + + // higgsbasis = {rho1, eta1, rho2, eta2, zeta1, psi1, zeta2, psi2} + for (std::size_t i = 0; i < NHiggs; + i++) // mass base index i corresponds to mass vector sorted in ascending + // mass + { + // charged submatrices + if (std::abs(HiggsRot(i, pos_rho1)) + std::abs(HiggsRot(i, pos_rho2)) > + ZeroThreshold) // use that mGpm < mHpm + { + if (pos_G1 == -1) { - GaugeBasis[i][j][k] = - TripleDeriv.at(i + j * NHiggs + k * NHiggs * NHiggs); + pos_G1 = i; + } + else + { + pos_H1 = i; + } + } + else if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > + ZeroThreshold) // use that mGpm < mHpm + { + if (pos_G2 == -1) + { + pos_G2 = i; + } + else + { + pos_H2 = i; + } + } + else if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) > + ZeroThreshold) // use that mh < mH + { + if (pos_h == -1) + { + pos_h = i; + } + else + { + pos_H = i; + } + } + else if (std::abs(HiggsRot(i, pos_psi1)) + std::abs(HiggsRot(i, pos_psi2)) > + ZeroThreshold) // use that 0 = mG0 < mA + { + if (pos_G0 == -1) + { + pos_G0 = i; + } + else + { + pos_A = i; } } } - MatrixXd HiggsRot(NHiggs, NHiggs); + // check if all position indices are set + if (pos_G1 == -1 or pos_G2 == -1 or pos_H1 == -1 or pos_H2 == -1 or pos_G0 == -1 or + pos_A == -1 or pos_H == -1 or pos_h == -1) + { + throw std::runtime_error("Error. Not all position indices are set."); + } + + // check if all other elements of rotation matrix are zero + bool zero_element = false; for (std::size_t i = 0; i < NHiggs; i++) { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + int ii = int(i); + int jj = int(j); + if (not((jj == pos_rho1 and (ii == pos_G1 or ii == pos_H1)) or + (jj == pos_eta1 and (ii == pos_G2 or ii == pos_H2)) or + (jj == pos_zeta1 and (ii == pos_h or ii == pos_H)) or + (jj == pos_psi1 and (ii == pos_G0 or ii == pos_A)) or + (jj == pos_rho2 and (ii == pos_G1 or ii == pos_H1)) or + (jj == pos_eta2 and (ii == pos_G2 or ii == pos_H2)) or + (jj == pos_zeta2 and (ii == pos_h or ii == pos_H)) or + (jj == pos_psi2 and (ii == pos_G0 or ii == pos_A)))) + { + zero_element = true; + } + if (zero_element and std::abs(HiggsRot(i, j)) > 0) + { + throw std::runtime_error("Error. Invalid rotation matrix detected."); + } + zero_element = false; } } - MatrixXd HiggsRotSort(NHiggs, NHiggs); - int posMHCS1 = 0, posMHCS2 = 0; - int posN[2] = {-1, -1}; - int countposN = 0; - int posG1 = 0, posG2 = 0, posG0 = 0; - int posA = 0, posh = 0, posH = 0; - double testsum = 0; - const double ZeroThreshold = 1e-5; - for (std::size_t i = 0; i < 3; i++) + // Determine the additional indices for the SM-like + // and lighter/heavier Higgses + pos_h_SM = -1, pos_h_H = -1; + + std::vector HiggsMasses; + HiggsMasses = HiggsMassesSquared(vevTree, 0); + + // Due to the masses being ordered, we will always have + // HiggsMasses[pos_h] <= HiggsMasses[pos_H] + double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h]) + - SMConstants.C_MassSMHiggs); + double diff2 = std::abs(std::sqrt(HiggsMasses[pos_H]) + - SMConstants.C_MassSMHiggs); + if (diff1 < diff2) + { + pos_h_SM = pos_h; + pos_h_H = pos_H; + } + else + { + pos_h_H = pos_h; + pos_h_SM = pos_H; + } + + MatrixXd HiggsRotFixed(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixed.row(i) = HiggsRot.row(i); + } + + // charged submatrix + if (HiggsRotFixed(pos_G1, pos_rho1) < 0) // G1 rho1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_G1) *= -1; + } + if (HiggsRotFixed(pos_G2, pos_eta1) < 0) // G2 eta1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_G2) *= -1; + } + if (HiggsRotFixed(pos_H1, pos_rho2) < 0) // H1 rho2 (+ cos(beta)) + { + HiggsRotFixed.row(pos_H1) *= -1; + } + if (HiggsRotFixed(pos_H2, pos_eta2) < 0) // H2 eta2 (+ cos(beta)) + { + HiggsRotFixed.row(pos_H2) *= -1; + } + + // check neutral, CP-odd submatrix + if (HiggsRotFixed(pos_G0, pos_psi1) < 0) // G0 psi1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_G0) *= -1; // G0 + } + if (HiggsRotFixed(pos_A, pos_psi2) < 0) // A psi2 (+ cos(beta)) + { + HiggsRotFixed.row(pos_A) *= -1; // A + } + + // // check neutral, CP-even submatrix + if (HiggsRotFixed(pos_H, pos_zeta1) < 0) // H zeta1 (+ cos(alpha)) { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posG1 = i; - testsum = std::abs(HiggsRot(i, 1)) + std::abs(HiggsRot(i, 3)); - if (testsum > ZeroThreshold) posG2 = i; - testsum = std::abs(HiggsRot(i, 5)) + std::abs(HiggsRot(i, 7)); - if (testsum > ZeroThreshold) posG0 = i; + // if negative, rotate H + HiggsRotFixed.row(pos_H) *= -1; // H } - for (std::size_t i = 3; i < NHiggs; i++) + if (HiggsRotFixed(pos_h, pos_zeta2) < 0) // h zeta2 (+ cos(alpha)) { - testsum = std::abs(HiggsRot(i, 0)) + std::abs(HiggsRot(i, 2)); - if (testsum > ZeroThreshold) posMHCS1 = i; - testsum = std::abs(HiggsRot(i, 1)) + std::abs(HiggsRot(i, 3)); - if (testsum > ZeroThreshold) posMHCS2 = i; - testsum = std::abs(HiggsRot(i, 5)) + std::abs(HiggsRot(i, 7)); - if (testsum > ZeroThreshold) posA = i; - testsum = std::abs(HiggsRot(i, 4)) + std::abs(HiggsRot(i, 6)); - if (testsum > ZeroThreshold) + // if negative, rotate h + HiggsRotFixed.row(pos_h) *= -1; // h + } + + // Extract the fixed mixing angle + alpha = std::asin(HiggsRotFixed(pos_H, pos_zeta2)); // H zeta2 (+ sin(alpha)) + + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) { - posN[countposN] = i; - countposN++; + HiggsRotationMatrixEnsuredConvention[i][j] = HiggsRotFixed(i, j); } } - posh = posN[0]; - posH = posN[1]; + return; +} + +/** + * Calculates the corrections to the Triple higgs couplings in the mass basis. + * + * Use the vector TripleHiggsCorrectionsCWPhysical to save your couplings and + * set the nTripleCouplings to the number of couplings you want as output. + */ +void Class_Potential_R2HDM::TripleHiggsCouplings() +{ + if (!SetCurvatureDone) SetCurvatureArrays(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + + MatrixXd HiggsRot(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; + } + } std::vector HiggsOrder(NHiggs); - HiggsOrder[0] = posG1; - HiggsOrder[1] = posG2; - HiggsOrder[2] = posMHCS1; - HiggsOrder[3] = posMHCS2; - HiggsOrder[4] = posG0; - HiggsOrder[5] = posA; - HiggsOrder[6] = posh; - HiggsOrder[7] = posH; + HiggsOrder[0] = pos_G1; + HiggsOrder[1] = pos_G2; + HiggsOrder[2] = pos_H1; + HiggsOrder[3] = pos_H2; + HiggsOrder[4] = pos_G0; + HiggsOrder[5] = pos_A; + HiggsOrder[6] = pos_h; + HiggsOrder[7] = pos_H; + MatrixXd HiggsRotSort(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); } + std::vector TripleDeriv; + TripleDeriv = WeinbergThirdDerivative(); + std::vector>> GaugeBasis( + NHiggs, + std::vector>(NHiggs, std::vector(NHiggs))); + for (std::size_t i = 0; i < NHiggs; i++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + for (std::size_t k = 0; k < NHiggs; k++) + { + GaugeBasis[i][j][k] = + TripleDeriv.at(i + j * NHiggs + k * NHiggs * NHiggs); + } + } + } + TripleHiggsCorrectionsCWPhysical.resize(NHiggs); TripleHiggsCorrectionsTreePhysical.resize(NHiggs); TripleHiggsCorrectionsCTPhysical.resize(NHiggs); diff --git a/src/models/ClassPotentialSM.cpp b/src/models/ClassPotentialSM.cpp index 37e95c747..0ec2cc42c 100644 --- a/src/models/ClassPotentialSM.cpp +++ b/src/models/ClassPotentialSM.cpp @@ -236,7 +236,7 @@ std::vector Class_SM::calc_CT() const retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; @@ -268,10 +268,17 @@ std::vector Class_SM::calc_CT() const return parCT; } +void Class_SM::AdjustRotationMatrix() +{ +} + void Class_SM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; std::vector HiggsOrder(NHiggs); diff --git a/src/models/ClassTemplate.cpp b/src/models/ClassTemplate.cpp index 938cafcd7..06a8f8609 100644 --- a/src/models/ClassTemplate.cpp +++ b/src/models/ClassTemplate.cpp @@ -245,7 +245,7 @@ std::vector Class_Template::calc_CT() const retmes += " was called before SetCurvatureArrays()!\n"; throw std::runtime_error(retmes); } - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { std::string retmes = __func__; retmes += " was called before CalculatePhysicalCouplings()!\n"; @@ -279,10 +279,23 @@ std::vector Class_Template::calc_CT() const return parCT; } +/** + * Ensures the correct rotation matrix convention + */ +void Class_Template::AdjustRotationMatrix() +{ + // Here you implement the rotation matrix convention of your model + // and define HiggsRotationMatrixEnsuredConvention, use then HiggsRotationMatrixEnsuredConvention in + // TripleHiggsCouplings +} + void Class_Template::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; std::vector HiggsOrder(NHiggs); // Here you have to set the vector HiggsOrder. By telling e.g. HiggsOrder[0] = diff --git a/src/utility/utility.cpp b/src/utility/utility.cpp index 300066a89..354123bb2 100644 --- a/src/utility/utility.cpp +++ b/src/utility/utility.cpp @@ -9,6 +9,7 @@ #include #include #include +#include /** * @file @@ -129,4 +130,62 @@ double EllipIntSecond(const double &x) return result; } +bool almost_the_same(const double &a, + const double &b, + const double &rel_precision, + const double &num_zero) +{ + if (std::abs(a) < num_zero and std::abs(b) < num_zero) + { + return true; + } + return std::abs(a - b) < std::abs(a + b) / 2 * rel_precision; +} + +bool almost_the_same(const std::complex &a, + const std::complex &b, + const double &rel_precision, + const double &num_zero) +{ + bool real_part = almost_the_same(a.real(), b.real(), rel_precision, + num_zero); + bool imag_part = almost_the_same(a.imag(), b.imag(), rel_precision, + num_zero); + return (real_part and imag_part); +} + +bool almost_the_same(const std::vector &a, + const std::vector &b, + const bool &allow_for_sign_flip, + const double &rel_precision, + const double &num_zero) +{ + if (a.size() != b.size()) + { + throw std::runtime_error("Error. Vectors must have the same size."); + } + int count_true = 0; + for (std::size_t i = 0; i < a.size(); i++) + { + if (allow_for_sign_flip) + { + count_true += + int(almost_the_same(a.at(i), b.at(i), rel_precision, num_zero)); + } + else + { + count_true += int(almost_the_same( + std::abs(a.at(i)), std::abs(b.at(i)), rel_precision, num_zero)); + } + } + if (std::size_t(count_true) == a.size()) + { + return true; + } + else + { + return false; + } +} + } // namespace BSMPT diff --git a/standalone/GenericModel.cpp b/standalone/GenericModel.cpp index 767f5cf3b..1f2bd8c0b 100644 --- a/standalone/GenericModel.cpp +++ b/standalone/GenericModel.cpp @@ -58,6 +58,7 @@ class Class_Potential_OriginDerived : public Class_Potential_Origin (void)v; return 0; }; + void AdjustRotationMatrix() override { return; }; void TripleHiggsCouplings() override { return; }; std::vector calc_CT() const override { return {0}; }; void Debugging(const std::vector &input,