From 50c1f8852c2815770acfd475554802810784590b Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Mon, 28 Oct 2024 20:19:36 +0100 Subject: [PATCH 01/15] Trilinear sign fix for C2HDM, WIP --- src/models/ClassPotentialC2HDM.cpp | 325 ++++++++++++++++++++++++++++- 1 file changed, 324 insertions(+), 1 deletion(-) diff --git a/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index e1e9aac39..42a6c216f 100644 --- a/src/models/ClassPotentialC2HDM.cpp +++ b/src/models/ClassPotentialC2HDM.cpp @@ -2396,6 +2396,329 @@ std::vector Class_Potential_C2HDM::calc_CT() const return parCT; } + +/** + * Ensures the correct rotation matrix convention + */ +void Class_Potential_C2HDM::AdjustRotationMatrix() +{ + // Here you implement the rotation matrix convention of your model + // and define HiggsRotationMatrixEnsuredConvention, use then + // HiggsRotationMatrixEnsuredConvention in TripleHiggsCouplings + + if (!SetCurvatureDone) SetCurvatureArrays(); + 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++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + } + } + + // 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; + + + const double ZeroThreshold = 1e-5; + + // 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 + + int row_G0 = -1; + for (std::size_t i = 0; i < 3; i++) + { + double psi1psi2 = std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)); + double non_psi1psi2 = std::abs(HiggsRot(i, pos_rho1)) + + std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_rho2)) + + std::abs(HiggsRot(i, pos_eta2)) + + std::abs(HiggsRot(i, pos_zeta1)) + + std::abs(HiggsRot(i, pos_zeta2)); + + if (psi1psi2 > ZeroThreshold and non_psi1psi2 < ZeroThreshold) + { + row_G0 = i; + break; + } + } + + if (row_G0 == -1) + { + throw std::runtime_error("Error. Something went wrong with finding + the row for the neutral Goldstone."); + } + + // Swap rows to ensure that G0 is always in the first row + // for the following manipulations + MatrixXd MoveNeutralGoldstoneFirst(NHiggs, NHiggs); + MoveNeutralGoldstoneFirst.setIdentity(NHiggs, NHiggs); + if (row_G0 != 0) + { + MoveNeutralGoldstoneFirst(0, 0) = 0.; + MoveNeutralGoldstoneFirst(row_G0, row_G0) = 0.; + MoveNeutralGoldstoneFirst(0, row_G0) = 1.; + MoveNeutralGoldstoneFirst(row_G0, 0) = 1.; + } + + + 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; + + // 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)) + + MatrixXd RotGoldstoneMassBasis(NHiggs, NHiggs); + RotGoldstoneMassBasis = MoveNeutralGoldstoneFirst*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., col1 = 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 + if (std::abs(std::abs(RotGoldstoneMassBasis(0, 0)) - 1.) > 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 + // CB: 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; + + // basis = {G^0, rho1, eta1, rho2, eta2, zeta1, zeta2, zeta3} + // 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 (pos_si_G1 == -1) + { + pos_si_G1 = i; + } + else + { + pos_si_H1 = i; + } + } + if (std::abs(RotGoldstoneMassBasis(i, pos_si_eta1)) + + std::abs(RotGoldstoneMassBasis(i, pos_si_eta2)) > ZeroThreshold) + // use that 0 = mGpm < mHpm + { + if (pos_si_G2 == -1) + { + pos_si_G2 = i; + } + else + { + pos_si_H2 = i; + } + } + + // Neutral submatrix (mixed CP-even and CP-odd states); + // neutral Goldstone already rotated out + 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 + { + if (pos_si_h1 == -1) { + pos_si_h1 = i; + } else if (pos_si_h2 == -1) { + pos_si_h2 = i; + } else { + pos_si_h3 = i; + } + } + } + + // 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 + bool 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++) + { + 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; + } + // Check against ZeroThreshold + if (zero_element and std::abs(RotGoldstoneMassBasis(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) = RotGoldstoneMassBasis.row(i); + } + + // Neutral Goldstone; flip sign if its 1 element 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 = MoveNeutralGoldstoneFirst*HiggsRotFixedGoldstone*RotRemoveGoldstone; + + + 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_C2HDM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); @@ -2427,7 +2750,7 @@ void Class_Potential_C2HDM::TripleHiggsCouplings() { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } MatrixXd HiggsRotSort(NHiggs, NHiggs); From 904c0e1d1d1a12864de0d4a7687954183922dcba Mon Sep 17 00:00:00 2001 From: Christoph Borschensky Date: Mon, 28 Oct 2024 23:26:24 +0100 Subject: [PATCH 02/15] Sign fix for rotation matrices for all currently implemented models --- include/BSMPT/models/ClassPotentialC2HDM.h | 1 + .../BSMPT/models/ClassPotentialCPintheDark.h | 1 + include/BSMPT/models/ClassPotentialCxSM.h | 1 + include/BSMPT/models/ClassPotentialN2HDM.h | 1 + include/BSMPT/models/ClassPotentialOrigin.h | 20 ++ include/BSMPT/models/ClassPotentialR2HDM.h | 1 + include/BSMPT/models/ClassPotentialSM.h | 1 + include/BSMPT/models/ClassTemplate.h | 1 + src/models/ClassPotentialC2HDM.cpp | 120 ++++----- src/models/ClassPotentialCPintheDark.cpp | 212 +++++++++++++++- src/models/ClassPotentialCxSM.cpp | 165 +++++++++++- src/models/ClassPotentialN2HDM.cpp | 239 +++++++++++++++++- src/models/ClassPotentialOrigin.cpp | 60 +++++ src/models/ClassPotentialR2HDM.cpp | 190 +++++++++++++- src/models/ClassPotentialSM.cpp | 7 + src/models/ClassTemplate.cpp | 13 + 16 files changed, 969 insertions(+), 64 deletions(-) diff --git a/include/BSMPT/models/ClassPotentialC2HDM.h b/include/BSMPT/models/ClassPotentialC2HDM.h index d2b50763d..a496f6ac5 100644 --- a/include/BSMPT/models/ClassPotentialC2HDM.h +++ b/include/BSMPT/models/ClassPotentialC2HDM.h @@ -121,6 +121,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..faa12b285 100644 --- a/include/BSMPT/models/ClassPotentialCPintheDark.h +++ b/include/BSMPT/models/ClassPotentialCPintheDark.h @@ -106,6 +106,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..38b05f19d 100644 --- a/include/BSMPT/models/ClassPotentialCxSM.h +++ b/include/BSMPT/models/ClassPotentialCxSM.h @@ -109,6 +109,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..cc7c1351f 100644 --- a/include/BSMPT/models/ClassPotentialN2HDM.h +++ b/include/BSMPT/models/ClassPotentialN2HDM.h @@ -125,6 +125,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..ee09dbf13 100644 --- a/include/BSMPT/models/ClassPotentialOrigin.h +++ b/include/BSMPT/models/ClassPotentialOrigin.h @@ -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,21 @@ class Class_Potential_Origin */ Eigen::MatrixXcd LeptonMassMatrix(const std::vector &v) const; + /** + * Ensures the correct rotation matrix convention + */ + virtual void AdjustRotationMatrix() = 0; + /** + * Returns true if two values are the same given some relative precision + */ + bool almost_the_same(double a, double b, double rel_precision = 0.01); + bool almost_the_same(std::complex a, + std::complex b, + double rel_precision = 0.01); + /** + * 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..38f6850cd 100644 --- a/include/BSMPT/models/ClassPotentialR2HDM.h +++ b/include/BSMPT/models/ClassPotentialR2HDM.h @@ -112,6 +112,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/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index 42a6c216f..a19cc2c2e 100644 --- a/src/models/ClassPotentialC2HDM.cpp +++ b/src/models/ClassPotentialC2HDM.cpp @@ -2402,12 +2402,10 @@ std::vector Class_Potential_C2HDM::calc_CT() const */ void Class_Potential_C2HDM::AdjustRotationMatrix() { - // Here you implement the rotation matrix convention of your model - // and define HiggsRotationMatrixEnsuredConvention, use then - // HiggsRotationMatrixEnsuredConvention in TripleHiggsCouplings + const double ZeroThreshold = 1e-5; if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsDone) CalculatePhysicalCouplings(); + if (!CalcCouplingsdone) CalculatePhysicalCouplings(); if (!CheckRotationMatrix()) // Check whether generically generated rotation // matrix is proper rotation matrix @@ -2424,13 +2422,28 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() } } - // interaction basis + // 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; - - const double ZeroThreshold = 1e-5; + // 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 @@ -2461,25 +2474,11 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() if (row_G0 == -1) { - throw std::runtime_error("Error. Something went wrong with finding - the row for the neutral Goldstone."); + throw std::runtime_error("Error. Something went wrong with finding " + "the row for the neutral Goldstone."); } - // Swap rows to ensure that G0 is always in the first row - // for the following manipulations - MatrixXd MoveNeutralGoldstoneFirst(NHiggs, NHiggs); - MoveNeutralGoldstoneFirst.setIdentity(NHiggs, NHiggs); - if (row_G0 != 0) - { - MoveNeutralGoldstoneFirst(0, 0) = 0.; - MoveNeutralGoldstoneFirst(row_G0, row_G0) = 0.; - MoveNeutralGoldstoneFirst(0, row_G0) = 1.; - MoveNeutralGoldstoneFirst(row_G0, 0) = 1.; - } - - 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.; @@ -2489,34 +2488,28 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() RotGoldstone.row(6) << 0., 0., 0., 0., 0., 0., 1., 0.; RotGoldstone.row(7) << 0., 0., 0., 0., 0., -C_SinBeta, 0., C_CosBeta; - // 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)) + // 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 (row_G0 != 0) + { + MoveGoldstoneFirst(0, 0) = 0.; + MoveGoldstoneFirst(row_G0, row_G0) = 0.; + MoveGoldstoneFirst(0, row_G0) = 1.; + MoveGoldstoneFirst(row_G0, 0) = 1.; + } MatrixXd RotGoldstoneMassBasis(NHiggs, NHiggs); - RotGoldstoneMassBasis = MoveNeutralGoldstoneFirst*HiggsRot*RotGoldstone.transpose(); + RotGoldstoneMassBasis = MoveGoldstoneFirst*HiggsRot*RotGoldstone.transpose(); // Semi-interaction basis (neutral Goldstone rotated out) - // {G^0, rho1, eta1, rho2, eta2, zeta1, zeta2, zeta3} + // 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., col1 = 0.; + 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++) @@ -2526,24 +2519,25 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() } // Consistency check that the Goldstone was rotated out properly - if (std::abs(std::abs(RotGoldstoneMassBasis(0, 0)) - 1.) > ZeroThreshold - or std::abs(row1) > ZeroThreshold or std::abs(col1) > ZeroThreshold) + 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."); + 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 - // CB: position of neutral Goldstone is fixed to 0, see above + // 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; - // basis = {G^0, rho1, eta1, rho2, eta2, zeta1, zeta2, zeta3} // 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 + 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 + // 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 @@ -2615,7 +2609,7 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() { zero_element = true; } - // Check against ZeroThreshold + if (zero_element and std::abs(RotGoldstoneMassBasis(i, j)) > ZeroThreshold) { throw std::runtime_error("Error. Invalid rotation matrix detected."); @@ -2630,8 +2624,8 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() HiggsRotFixed.row(i) = RotGoldstoneMassBasis.row(i); } - // Neutral Goldstone; flip sign if its 1 element is negative - if (HiggsRotFixed(pos_si_G0, pos_si_G0) < 0) // G0 G0 (= +1) + // 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; } @@ -2654,7 +2648,6 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() 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): @@ -2698,14 +2691,21 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() 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); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixedGoldstone.row(i) = HiggsRotFixed.row(i); + } - HiggsRotFixed = MoveNeutralGoldstoneFirst*HiggsRotFixedGoldstone*RotRemoveGoldstone; + HiggsRotFixed = MoveGoldstoneFirst*HiggsRotFixedGoldstone*RotGoldstone; + // Extract the fixed mixing angles + double sina2 = HiggsRotFixedGoldstone(pos_si_h1, pos_zeta3); // +sin(a2) + double cosa2 = std::sqrt(1.0 - sina2*sina2); + alpha1 = std::asin(HiggsRotFixedGoldstone(pos_si_h1, pos_zeta2)/cosa2); // +sin(a1) cos(a2) + alpha2 = std::asin(sina2); + alpha3 = std::asin(HiggsRotFixedGoldstone(pos_si_h2, pos_zeta3)/cosa2); // +cos(a2) sin(a3) for (std::size_t i = 0; i < NHiggs; i++) { diff --git a/src/models/ClassPotentialCPintheDark.cpp b/src/models/ClassPotentialCPintheDark.cpp index d13d99b4c..db8dbd6a2 100644 --- a/src/models/ClassPotentialCPintheDark.cpp +++ b/src/models/ClassPotentialCPintheDark.cpp @@ -1314,12 +1314,222 @@ std::vector Class_Potential_CPintheDark::calc_CT() const return parCT; } +/** + * Ensures the correct rotation matrix convention + */ +void Class_Potential_CPintheDark::AdjustRotationMatrix() +{ + const double ZeroThreshold = 1e-5; + + if (!SetCurvatureDone) SetCurvatureArrays(); + 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++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + } + } + + // 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 semi-interaction to mass basis + int pos_Gp = -1, pos_Gm = -1, pos_Hp = -1, pos_Hm = -1, pos_HSM = -1; + int pos_G0 = -1, pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; + + + // 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, pos_zeta2)) + std::abs(HiggsRot(i, pos_psi2)) + + std::abs(HiggsRot(i, pos_rhoS))) > ZeroThreshold) + { + // use that scalars are sorted by mass + 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_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; +} + // mass basis triple couplings void Class_Potential_CPintheDark::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + // position indices store the position of the physical fields std::size_t posGp = 0; std::size_t posGm = 0; @@ -1336,7 +1546,7 @@ void Class_Potential_CPintheDark::TripleHiggsCouplings() { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } diff --git a/src/models/ClassPotentialCxSM.cpp b/src/models/ClassPotentialCxSM.cpp index c5f5baa40..e2422ceb6 100644 --- a/src/models/ClassPotentialCxSM.cpp +++ b/src/models/ClassPotentialCxSM.cpp @@ -765,11 +765,174 @@ std::vector Class_CxSM::calc_CT() const return parCT; } +/** + * Ensures the correct rotation matrix convention + */ +void Class_CxSM::AdjustRotationMatrix() +{ + const double ZeroThreshold = 1e-5; + + if (!SetCurvatureDone) SetCurvatureArrays(); + 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++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + } + } + + // CxSM interaction basis + // 0 1 2 3 4 5 + // Gp, Gm, G0, zeta1, zeta2, zeta3 + int pos_Gp = 0, pos_Gm = 1, pos_G0 = 2, + pos_zeta1 = 3, pos_zeta2 = 4, pos_zeta3 = 5; + + // Indices of mass eigenstates for rotation from interaction to mass basis + // Only neutral/physical part necessary, as Goldstones are already diagonal + int pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; + + // 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; + } + } + + 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; + 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 @@ -780,7 +943,7 @@ void Class_CxSM::TripleHiggsCouplings() { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } diff --git a/src/models/ClassPotentialN2HDM.cpp b/src/models/ClassPotentialN2HDM.cpp index 2caed38ef..bc6dd3d2f 100644 --- a/src/models/ClassPotentialN2HDM.cpp +++ b/src/models/ClassPotentialN2HDM.cpp @@ -772,11 +772,248 @@ std::vector Class_Potential_N2HDM::calc_CT() const return parCT; } +/** + * Ensures the correct rotation matrix convention + */ +void Class_Potential_N2HDM::AdjustRotationMatrix() +{ + const double ZeroThreshold = 1e-5; + + if (!SetCurvatureDone) SetCurvatureArrays(); + 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++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + } + } + + // 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; + + // Unlike in the C2HDM, there is no CP-mixing, so we do not have to rotate out the Goldstone + + // Indices of mass eigenstates for rotation from interaction to mass basis + int pos_G0 = -1, pos_G1 = -1, pos_G2 = -1, pos_H1 = -1, pos_H2 = -1; + int 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 + { + // Charged submatrix + if (std::abs(HiggsRot(i, pos_rho1)) + + std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) + { + pos_G1 = i; + } + if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_G2 = i; + } + + // Neutral CP-odd submatrix + if (std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) + { + pos_G0 = i; + } + } + + // 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; + } + if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_H2 = i; + } + + // Neutral CP-odd submatrix + if (std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) + { + pos_A = i; + } + + // Neutral CP-even submatrix + 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++) + { + 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_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; + } + } + + MatrixXd HiggsRotFixed(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixed.row(i) = HiggsRot.row(i); + } + + // Neutral CP-odd submatrix + if (HiggsRotFixed(pos_G0, pos_psi1) < 0) // G0 psi1 (+ cos(beta)) + { + HiggsRotFixed.row(pos_G0) *= -1; + } + if (HiggsRotFixed(pos_A, pos_psi2) < 0) // A psi2 (+ cos(beta)) + { + HiggsRotFixed.row(pos_A) *= -1; + } + + // 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-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 + { + // if negative, flip sign of h2 + HiggsRotFixed.row(pos_h2) *= -1; + } + + // Extract the fixed mixing angles + double sina2 = HiggsRotFixed(pos_si_h1, pos_zeta3); // +sin(a2) + double cosa2 = std::sqrt(1.0 - sina2*sina2); + alpha1 = std::asin(HiggsRotFixed(pos_si_h1, pos_zeta2)/cosa2); // +sin(a1) cos(a2) + alpha2 = std::asin(sina2); + alpha3 = std::asin(HiggsRotFixed(pos_si_h2, pos_zeta3)/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; + std::vector TripleDeriv; TripleDeriv = WeinbergThirdDerivative(); std::vector>> GaugeBasis( @@ -799,7 +1036,7 @@ void Class_Potential_N2HDM::TripleHiggsCouplings() { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } diff --git a/src/models/ClassPotentialOrigin.cpp b/src/models/ClassPotentialOrigin.cpp index 271e7d059..e3e8594fc 100644 --- a/src/models/ClassPotentialOrigin.cpp +++ b/src/models/ClassPotentialOrigin.cpp @@ -959,6 +959,64 @@ Class_Potential_Origin::SecondDerivativeOfEigenvaluesNonRepeated( return res; } +bool Class_Potential_Origin::almost_the_same(double a, + double b, + double rel_precision) +{ + if (std::abs(a) < 1e-10 and std::abs(b) < 1e-10) + { + return true; + } + return std::abs(a - b) < std::abs(a + b) / 2 * rel_precision; +} + +bool Class_Potential_Origin::almost_the_same(std::complex a, + std::complex b, + double rel_precision) +{ + bool real_part = almost_the_same(a.real(), b.real(), rel_precision); + bool imag_part = almost_the_same(a.imag(), b.imag(), rel_precision); + return (real_part and imag_part); +} + +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 DetIsOne = almost_the_same(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 (DetIsOne and InvEqTrans) + { + return true; + } + return false; +} + void Class_Potential_Origin::CalculatePhysicalCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); @@ -3778,6 +3836,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..93eddeb8b 100644 --- a/src/models/ClassPotentialR2HDM.cpp +++ b/src/models/ClassPotentialR2HDM.cpp @@ -742,6 +742,191 @@ std::vector Class_Potential_R2HDM::calc_CT() const return parCT; } +/** + * Ensures the correct rotation matrix convention + */ +void Class_Potential_R2HDM::AdjustRotationMatrix() +{ + const double ZeroThreshold = 1e-5; + + if (!SetCurvatureDone) SetCurvatureArrays(); + 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++) + { + for (std::size_t j = 0; j < NHiggs; j++) + { + HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + } + } + + // initialize position indices (new initialization for each point in multiline + // files) + int posG1 = -1, posG2 = -1, posH1 = -1, posH2 = -1, posG0 = -1, posA = -1, + posH = -1, posh = -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 (posG1 == -1) + { + posG1 = i; + } + else + { + posH1 = i; + } + } + if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > + ZeroThreshold) // use that mGpm < mHpm + { + if (posG2 == -1) + { + posG2 = i; + } + else + { + posH2 = i; + } + } + if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) > + ZeroThreshold) // use that mh < mH + { + if (posh == -1) + { + posh = i; + } + else + { + posH = i; + } + } + if (std::abs(HiggsRot(i, pos_psi1)) + std::abs(HiggsRot(i, pos_psi2)) > + ZeroThreshold) // use that 0 = mG0 < mA + { + if (posG0 == -1) + { + posG0 = i; + } + else + { + posA = i; + } + } + } + + // check if all position indices are set + if (posG1 == -1 or posG2 == -1 or posH1 == -1 or posH2 == -1 or posG0 == -1 or + posA == -1 or posH == -1 or posh == -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 == posG1 or ii == posH1)) or + (jj == pos_eta1 and (ii == posG2 or ii == posH2)) or + (jj == pos_zeta1 and (ii == posh or ii == posH)) or + (jj == pos_psi1 and (ii == posG0 or ii == posA)) or + (jj == pos_rho2 and (ii == posG1 or ii == posH1)) or + (jj == pos_eta2 and (ii == posG2 or ii == posH2)) or + (jj == pos_zeta2 and (ii == posh or ii == posH)) or + (jj == pos_psi2 and (ii == posG0 or ii == posA)))) + { + 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 HiggsRotFixed(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) + { + HiggsRotFixed.row(i) = HiggsRot.row(i); + } + + // charged submatrix + if (HiggsRotFixed(posG1, pos_rho1) < 0) // G1 rho1 (+ cos(beta)) + { + HiggsRotFixed.row(posG1) *= -1; + } + if (HiggsRotFixed(posG2, pos_eta1) < 0) // G2 eta1 (+ cos(beta)) + { + HiggsRotFixed.row(posG2) *= -1; + } + if (HiggsRotFixed(posH1, pos_rho2) < 0) // H1 rho2 (+ cos(beta)) + { + HiggsRotFixed.row(posH1) *= -1; + } + if (HiggsRotFixed(posH2, pos_eta2) < 0) // H2 eta2 (+ cos(beta)) + { + HiggsRotFixed.row(posH2) *= -1; + } + + // check neutral, CP-odd submatrix + if (HiggsRotFixed(posG0, pos_psi1) < 0) // G0 psi1 (+ cos(beta)) + { + HiggsRotFixed.row(posG0) *= -1; // G0 + } + if (HiggsRotFixed(posA, pos_psi2) < 0) // A psi2 (+ cos(beta)) + { + HiggsRotFixed.row(posA) *= -1; // A + } + + // // check neutral, CP-even submatrix + if (HiggsRotFixed(posH, pos_zeta1) < 0) // H zeta1 (+ cos(alpha)) + { + // if negative, rotate H + HiggsRotFixed.row(posH) *= -1; // H + } + if (HiggsRotFixed(posh, pos_zeta2) < 0) // h zeta2 (+ cos(alpha)) + { + // if negative, rotate h + HiggsRotFixed.row(posh) *= -1; // h + } + + // Extract the fixed mixing angle + alpha = std::asin(HiggsRotFixed(posH, pos_zeta2)); // H zeta2 (+ sin(alpha)) + + 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; +} + /** * Calculates the corrections to the Triple higgs couplings in the mass basis. * @@ -753,6 +938,9 @@ void Class_Potential_R2HDM::TripleHiggsCouplings() if (!SetCurvatureDone) SetCurvatureArrays(); if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + std::vector TripleDeriv; TripleDeriv = WeinbergThirdDerivative(); std::vector>> GaugeBasis( @@ -775,7 +963,7 @@ void Class_Potential_R2HDM::TripleHiggsCouplings() { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } diff --git a/src/models/ClassPotentialSM.cpp b/src/models/ClassPotentialSM.cpp index 37e95c747..8244a8b70 100644 --- a/src/models/ClassPotentialSM.cpp +++ b/src/models/ClassPotentialSM.cpp @@ -268,11 +268,18 @@ std::vector Class_SM::calc_CT() const return parCT; } +void Class_SM::AdjustRotationMatrix() +{ +} + void Class_SM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + std::vector HiggsOrder(NHiggs); for (std::size_t i = 0; i < NHiggs; i++) diff --git a/src/models/ClassTemplate.cpp b/src/models/ClassTemplate.cpp index 938cafcd7..d934d1d8c 100644 --- a/src/models/ClassTemplate.cpp +++ b/src/models/ClassTemplate.cpp @@ -279,11 +279,24 @@ 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 (CalculatedTripleCopulings) return; + CalculatedTripleCopulings = true; + 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 From 447a15b499b0bab2565ae317b3a285f709cd93b6 Mon Sep 17 00:00:00 2001 From: Christoph Borschensky Date: Mon, 28 Oct 2024 23:28:18 +0100 Subject: [PATCH 03/15] Renamed variable CalcCouplingsdone to CalcCouplingsDone for consistent naming --- include/BSMPT/models/ClassPotentialOrigin.h | 4 ++-- src/models/ClassPotentialC2HDM.cpp | 8 ++++---- src/models/ClassPotentialCPintheDark.cpp | 6 +++--- src/models/ClassPotentialCxSM.cpp | 8 ++++---- src/models/ClassPotentialN2HDM.cpp | 6 +++--- src/models/ClassPotentialOrigin.cpp | 12 ++++++------ src/models/ClassPotentialR2HDM.cpp | 6 +++--- src/models/ClassPotentialSM.cpp | 4 ++-- src/models/ClassTemplate.cpp | 4 ++-- 9 files changed, 29 insertions(+), 29 deletions(-) diff --git a/include/BSMPT/models/ClassPotentialOrigin.h b/include/BSMPT/models/ClassPotentialOrigin.h index ee09dbf13..6193bdc0a 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 diff --git a/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index a19cc2c2e..a8df7d1fd 100644 --- a/src/models/ClassPotentialC2HDM.cpp +++ b/src/models/ClassPotentialC2HDM.cpp @@ -2128,7 +2128,7 @@ 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++) @@ -2294,7 +2294,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"; @@ -2405,7 +2405,7 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() const double ZeroThreshold = 1e-5; if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (!CheckRotationMatrix()) // Check whether generically generated rotation // matrix is proper rotation matrix @@ -2722,7 +2722,7 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() void Class_Potential_C2HDM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; diff --git a/src/models/ClassPotentialCPintheDark.cpp b/src/models/ClassPotentialCPintheDark.cpp index db8dbd6a2..8c7b49499 100644 --- a/src/models/ClassPotentialCPintheDark.cpp +++ b/src/models/ClassPotentialCPintheDark.cpp @@ -1258,7 +1258,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"; @@ -1322,7 +1322,7 @@ void Class_Potential_CPintheDark::AdjustRotationMatrix() const double ZeroThreshold = 1e-5; if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (!CheckRotationMatrix()) // Check whether generically generated rotation // matrix is proper rotation matrix @@ -1525,7 +1525,7 @@ void Class_Potential_CPintheDark::AdjustRotationMatrix() void Class_Potential_CPintheDark::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; diff --git a/src/models/ClassPotentialCxSM.cpp b/src/models/ClassPotentialCxSM.cpp index e2422ceb6..f554f0a29 100644 --- a/src/models/ClassPotentialCxSM.cpp +++ b/src/models/ClassPotentialCxSM.cpp @@ -589,7 +589,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"; @@ -773,7 +773,7 @@ void Class_CxSM::AdjustRotationMatrix() const double ZeroThreshold = 1e-5; if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (!CheckRotationMatrix()) // Check whether generically generated rotation // matrix is proper rotation matrix @@ -928,7 +928,7 @@ void Class_CxSM::AdjustRotationMatrix() void Class_CxSM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; @@ -1468,7 +1468,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 bc6dd3d2f..2cf024633 100644 --- a/src/models/ClassPotentialN2HDM.cpp +++ b/src/models/ClassPotentialN2HDM.cpp @@ -623,7 +623,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"; @@ -780,7 +780,7 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() const double ZeroThreshold = 1e-5; if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (!CheckRotationMatrix()) // Check whether generically generated rotation // matrix is proper rotation matrix @@ -1009,7 +1009,7 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() void Class_Potential_N2HDM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; diff --git a/src/models/ClassPotentialOrigin.cpp b/src/models/ClassPotentialOrigin.cpp index e3e8594fc..d6b85d49e 100644 --- a/src/models/ClassPotentialOrigin.cpp +++ b/src/models/ClassPotentialOrigin.cpp @@ -1453,7 +1453,7 @@ void Class_Potential_Origin::CalculatePhysicalCouplings() } } - CalcCouplingsdone = true; + CalcCouplingsDone = true; return; } @@ -1461,7 +1461,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__; @@ -1555,7 +1555,7 @@ std::vector Class_Potential_Origin::WeinbergFirstDerivative() const Eigen::MatrixXd Class_Potential_Origin::WeinbergSecondDerivativeAsMatrixXd() const { - if (!CalcCouplingsdone) + if (!CalcCouplingsDone) { // CalculatePhysicalCouplings(); std::string retmes = __func__; @@ -1714,7 +1714,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."; @@ -1960,7 +1960,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."; @@ -3530,7 +3530,7 @@ void Class_Potential_Origin::sym4Dim( void Class_Potential_Origin::resetbools() { SetCurvatureDone = false; - CalcCouplingsdone = false; + CalcCouplingsDone = false; CalculatedTripleCopulings = false; parStored.clear(); parCTStored.clear(); diff --git a/src/models/ClassPotentialR2HDM.cpp b/src/models/ClassPotentialR2HDM.cpp index 93eddeb8b..b110ee00a 100644 --- a/src/models/ClassPotentialR2HDM.cpp +++ b/src/models/ClassPotentialR2HDM.cpp @@ -642,7 +642,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"; @@ -750,7 +750,7 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() const double ZeroThreshold = 1e-5; if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (!CheckRotationMatrix()) // Check whether generically generated rotation // matrix is proper rotation matrix @@ -936,7 +936,7 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() void Class_Potential_R2HDM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; diff --git a/src/models/ClassPotentialSM.cpp b/src/models/ClassPotentialSM.cpp index 8244a8b70..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"; @@ -275,7 +275,7 @@ void Class_SM::AdjustRotationMatrix() void Class_SM::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; diff --git a/src/models/ClassTemplate.cpp b/src/models/ClassTemplate.cpp index d934d1d8c..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"; @@ -292,7 +292,7 @@ void Class_Template::AdjustRotationMatrix() void Class_Template::TripleHiggsCouplings() { if (!SetCurvatureDone) SetCurvatureArrays(); - if (!CalcCouplingsdone) CalculatePhysicalCouplings(); + if (!CalcCouplingsDone) CalculatePhysicalCouplings(); if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; From a912bc5a818176fedf9741e77d0bd90c9293c2ae Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Tue, 29 Oct 2024 09:07:56 +0100 Subject: [PATCH 04/15] Fixed some comments in AdjustMixingMatrix() of the C2HDM --- src/models/ClassPotentialC2HDM.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index a8df7d1fd..3a1b68f09 100644 --- a/src/models/ClassPotentialC2HDM.cpp +++ b/src/models/ClassPotentialC2HDM.cpp @@ -2478,6 +2478,7 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() "the row for the neutral Goldstone."); } + // 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.; @@ -2500,6 +2501,8 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() MoveGoldstoneFirst(row_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(); @@ -2518,7 +2521,8 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() col1 += std::abs(RotGoldstoneMassBasis(i, 0)); } - // Consistency check that the Goldstone was rotated out properly + // 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) { @@ -2531,7 +2535,7 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() 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 + // 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 { From 5004581230c434d50f22224c7b40a22a7a00d412 Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Tue, 29 Oct 2024 09:24:49 +0100 Subject: [PATCH 05/15] Fixed variable names for mass/gauge state indices in AdjustRotationMatrix() --- src/models/ClassPotentialC2HDM.cpp | 6 +++--- src/models/ClassPotentialN2HDM.cpp | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index 3a1b68f09..e2b105e11 100644 --- a/src/models/ClassPotentialC2HDM.cpp +++ b/src/models/ClassPotentialC2HDM.cpp @@ -2705,11 +2705,11 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() HiggsRotFixed = MoveGoldstoneFirst*HiggsRotFixedGoldstone*RotGoldstone; // Extract the fixed mixing angles - double sina2 = HiggsRotFixedGoldstone(pos_si_h1, pos_zeta3); // +sin(a2) + 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_zeta2)/cosa2); // +sin(a1) cos(a2) + 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_zeta3)/cosa2); // +cos(a2) sin(a3) + alpha3 = std::asin(HiggsRotFixedGoldstone(pos_si_h2, pos_si_zeta3)/cosa2); // +cos(a2) sin(a3) for (std::size_t i = 0; i < NHiggs; i++) { diff --git a/src/models/ClassPotentialN2HDM.cpp b/src/models/ClassPotentialN2HDM.cpp index 2cf024633..d242d7bf7 100644 --- a/src/models/ClassPotentialN2HDM.cpp +++ b/src/models/ClassPotentialN2HDM.cpp @@ -989,11 +989,11 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() } // Extract the fixed mixing angles - double sina2 = HiggsRotFixed(pos_si_h1, pos_zeta3); // +sin(a2) + double sina2 = HiggsRotFixed(pos_h1, pos_rhoS); // +sin(a2) double cosa2 = std::sqrt(1.0 - sina2*sina2); - alpha1 = std::asin(HiggsRotFixed(pos_si_h1, pos_zeta2)/cosa2); // +sin(a1) cos(a2) + alpha1 = std::asin(HiggsRotFixed(pos_h1, pos_zeta2)/cosa2); // +sin(a1) cos(a2) alpha2 = std::asin(sina2); - alpha3 = std::asin(HiggsRotFixed(pos_si_h2, pos_zeta3)/cosa2); // +cos(a2) sin(a3) + alpha3 = std::asin(HiggsRotFixed(pos_h2, pos_rhoS)/cosa2); // +cos(a2) sin(a3) for (std::size_t i = 0; i < NHiggs; i++) { From 322549c43822f7072b8e5552d8c0ef55f080c357 Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Tue, 29 Oct 2024 09:25:25 +0100 Subject: [PATCH 06/15] Added empty AdjustRotationMatrix() to standalone/GenericModel.cpp --- standalone/GenericModel.cpp | 1 + 1 file changed, 1 insertion(+) 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, From 4d24c9380bc0dabc91880dbc4d67c16f63d0361c Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Tue, 29 Oct 2024 10:00:20 +0100 Subject: [PATCH 07/15] Modified CheckRotationMatrix() to allow for det=+1 or -1; initialised HiggsRotationMatrixEnsuredConvention vector --- src/models/ClassPotentialOrigin.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/models/ClassPotentialOrigin.cpp b/src/models/ClassPotentialOrigin.cpp index d6b85d49e..7fde4749a 100644 --- a/src/models/ClassPotentialOrigin.cpp +++ b/src/models/ClassPotentialOrigin.cpp @@ -979,6 +979,9 @@ bool Class_Potential_Origin::almost_the_same(std::complex a, return (real_part and imag_part); } +// 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); @@ -992,7 +995,8 @@ bool Class_Potential_Origin::CheckRotationMatrix() double precision = 1e-10; - bool DetIsOne = almost_the_same(mat.determinant(), 1., precision); + bool AbsDetIsOne = almost_the_same(std::abs(mat.determinant()), 1., + precision); bool InvEqTrans = true; auto inv = mat.inverse(); @@ -1010,7 +1014,7 @@ bool Class_Potential_Origin::CheckRotationMatrix() } } - if (DetIsOne and InvEqTrans) + if (AbsDetIsOne and InvEqTrans) { return true; } @@ -3445,6 +3449,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( From cc1db38a6dc5dee2d9bc1c3809b7cc8dc7466546 Mon Sep 17 00:00:00 2001 From: Christoph Borschensky Date: Thu, 31 Oct 2024 01:00:46 +0100 Subject: [PATCH 08/15] Calculate mass indices only once, WIP --- include/BSMPT/models/ClassPotentialN2HDM.h | 3 + include/BSMPT/models/ClassPotentialR2HDM.h | 2 + src/models/ClassPotentialN2HDM.cpp | 106 ++++------- src/models/ClassPotentialR2HDM.cpp | 193 +++++++-------------- 4 files changed, 100 insertions(+), 204 deletions(-) diff --git a/include/BSMPT/models/ClassPotentialN2HDM.h b/include/BSMPT/models/ClassPotentialN2HDM.h index cc7c1351f..96917ed72 100644 --- a/include/BSMPT/models/ClassPotentialN2HDM.h +++ b/include/BSMPT/models/ClassPotentialN2HDM.h @@ -99,6 +99,9 @@ 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; + void ReadAndSet(const std::string &linestr, std::vector &par) override; std::vector addLegendCT() const override; diff --git a/include/BSMPT/models/ClassPotentialR2HDM.h b/include/BSMPT/models/ClassPotentialR2HDM.h index 38f6850cd..77b747c07 100644 --- a/include/BSMPT/models/ClassPotentialR2HDM.h +++ b/include/BSMPT/models/ClassPotentialR2HDM.h @@ -90,6 +90,8 @@ 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; + void ReadAndSet(const std::string &linestr, std::vector &par) override; std::vector addLegendCT() const override; diff --git a/src/models/ClassPotentialN2HDM.cpp b/src/models/ClassPotentialN2HDM.cpp index d242d7bf7..4d8ad31c6 100644 --- a/src/models/ClassPotentialN2HDM.cpp +++ b/src/models/ClassPotentialN2HDM.cpp @@ -807,8 +807,8 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() // Unlike in the C2HDM, there is no CP-mixing, so we do not have to rotate out the Goldstone // Indices of mass eigenstates for rotation from interaction to mass basis - int pos_G0 = -1, pos_G1 = -1, pos_G2 = -1, pos_H1 = -1, pos_H2 = -1; - int pos_h1 = -1, pos_h2 = -1, pos_h3 = -1, pos_A = -1; + 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 @@ -1040,93 +1040,49 @@ 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; - - 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) 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; - } - 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, 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) - { - posN.at(countposN) = i; - countposN++; - } - testsum = std::abs(HiggsRot(i, 4)) + std::abs(HiggsRot(i, 5)); - if (testsum > ZeroThreshold) posA = i; - } + int 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++) + double MSM = 125.09; + double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h1]) - MSM); + double diff2 = std::abs(std::sqrt(HiggsMasses[pos_h2]) - MSM); + double diff3 = std::abs(std::sqrt(HiggsMasses[pos_h3]) - MSM); + 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) { - 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]); + pos_h_SM = pos_h2; + pos_h_l = pos_h1; + pos_h_H = pos_h3; } else { - MhUp = std::sqrt(NeutralHiggs[1]); - MhDown = std::sqrt(NeutralHiggs[0]); + pos_h_SM = pos_h3; + pos_h_l = pos_h1; + pos_h_H = pos_h2; } - 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; - } + 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; - 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++) + MatrixXd HiggsRotSort(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) { - HiggsRotSort.row(i) = HiggsRot.row(posN[i - 6]); + HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); } TripleHiggsCorrectionsCWPhysical.resize(NHiggs); diff --git a/src/models/ClassPotentialR2HDM.cpp b/src/models/ClassPotentialR2HDM.cpp index b110ee00a..e964eddf8 100644 --- a/src/models/ClassPotentialR2HDM.cpp +++ b/src/models/ClassPotentialR2HDM.cpp @@ -567,64 +567,34 @@ void Class_Potential_R2HDM::write() const } } - 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++; - } - 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"; + << "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()); } @@ -769,8 +739,8 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() // initialize position indices (new initialization for each point in multiline // files) - int posG1 = -1, posG2 = -1, posH1 = -1, posH2 = -1, posG0 = -1, posA = -1, - posH = -1, posh = -1; + 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 @@ -786,56 +756,56 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() if (std::abs(HiggsRot(i, pos_rho1)) + std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) // use that mGpm < mHpm { - if (posG1 == -1) + if (pos_G1 == -1) { - posG1 = i; + pos_G1 = i; } else { - posH1 = i; + pos_H1 = i; } } if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) // use that mGpm < mHpm { - if (posG2 == -1) + if (pos_G2 == -1) { - posG2 = i; + pos_G2 = i; } else { - posH2 = i; + pos_H2 = i; } } if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) > ZeroThreshold) // use that mh < mH { - if (posh == -1) + if (pos_h == -1) { - posh = i; + pos_h = i; } else { - posH = i; + pos_H = i; } } if (std::abs(HiggsRot(i, pos_psi1)) + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) // use that 0 = mG0 < mA { - if (posG0 == -1) + if (pos_G0 == -1) { - posG0 = i; + pos_G0 = i; } else { - posA = i; + pos_A = i; } } } // check if all position indices are set - if (posG1 == -1 or posG2 == -1 or posH1 == -1 or posH2 == -1 or posG0 == -1 or - posA == -1 or posH == -1 or posh == -1) + 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."); } @@ -848,14 +818,14 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() { int ii = int(i); int jj = int(j); - if (not((jj == pos_rho1 and (ii == posG1 or ii == posH1)) or - (jj == pos_eta1 and (ii == posG2 or ii == posH2)) or - (jj == pos_zeta1 and (ii == posh or ii == posH)) or - (jj == pos_psi1 and (ii == posG0 or ii == posA)) or - (jj == pos_rho2 and (ii == posG1 or ii == posH1)) or - (jj == pos_eta2 and (ii == posG2 or ii == posH2)) or - (jj == pos_zeta2 and (ii == posh or ii == posH)) or - (jj == pos_psi2 and (ii == posG0 or ii == posA)))) + 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; } @@ -874,47 +844,47 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() } // charged submatrix - if (HiggsRotFixed(posG1, pos_rho1) < 0) // G1 rho1 (+ cos(beta)) + if (HiggsRotFixed(pos_G1, pos_rho1) < 0) // G1 rho1 (+ cos(beta)) { - HiggsRotFixed.row(posG1) *= -1; + HiggsRotFixed.row(pos_G1) *= -1; } - if (HiggsRotFixed(posG2, pos_eta1) < 0) // G2 eta1 (+ cos(beta)) + if (HiggsRotFixed(pos_G2, pos_eta1) < 0) // G2 eta1 (+ cos(beta)) { - HiggsRotFixed.row(posG2) *= -1; + HiggsRotFixed.row(pos_G2) *= -1; } - if (HiggsRotFixed(posH1, pos_rho2) < 0) // H1 rho2 (+ cos(beta)) + if (HiggsRotFixed(pos_H1, pos_rho2) < 0) // H1 rho2 (+ cos(beta)) { - HiggsRotFixed.row(posH1) *= -1; + HiggsRotFixed.row(pos_H1) *= -1; } - if (HiggsRotFixed(posH2, pos_eta2) < 0) // H2 eta2 (+ cos(beta)) + if (HiggsRotFixed(pos_H2, pos_eta2) < 0) // H2 eta2 (+ cos(beta)) { - HiggsRotFixed.row(posH2) *= -1; + HiggsRotFixed.row(pos_H2) *= -1; } // check neutral, CP-odd submatrix - if (HiggsRotFixed(posG0, pos_psi1) < 0) // G0 psi1 (+ cos(beta)) + if (HiggsRotFixed(pos_G0, pos_psi1) < 0) // G0 psi1 (+ cos(beta)) { - HiggsRotFixed.row(posG0) *= -1; // G0 + HiggsRotFixed.row(pos_G0) *= -1; // G0 } - if (HiggsRotFixed(posA, pos_psi2) < 0) // A psi2 (+ cos(beta)) + if (HiggsRotFixed(pos_A, pos_psi2) < 0) // A psi2 (+ cos(beta)) { - HiggsRotFixed.row(posA) *= -1; // A + HiggsRotFixed.row(pos_A) *= -1; // A } // // check neutral, CP-even submatrix - if (HiggsRotFixed(posH, pos_zeta1) < 0) // H zeta1 (+ cos(alpha)) + if (HiggsRotFixed(pos_H, pos_zeta1) < 0) // H zeta1 (+ cos(alpha)) { // if negative, rotate H - HiggsRotFixed.row(posH) *= -1; // H + HiggsRotFixed.row(pos_H) *= -1; // H } - if (HiggsRotFixed(posh, pos_zeta2) < 0) // h zeta2 (+ cos(alpha)) + if (HiggsRotFixed(pos_h, pos_zeta2) < 0) // h zeta2 (+ cos(alpha)) { // if negative, rotate h - HiggsRotFixed.row(posh) *= -1; // h + HiggsRotFixed.row(pos_h) *= -1; // h } // Extract the fixed mixing angle - alpha = std::asin(HiggsRotFixed(posH, pos_zeta2)); // H zeta2 (+ sin(alpha)) + alpha = std::asin(HiggsRotFixed(pos_H, pos_zeta2)); // H zeta2 (+ sin(alpha)) for (std::size_t i = 0; i < NHiggs; i++) { @@ -967,52 +937,17 @@ void Class_Potential_R2HDM::TripleHiggsCouplings() } } - 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++) - { - 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, 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) - { - posN[countposN] = i; - countposN++; - } - } - - posh = posN[0]; - posH = posN[1]; - 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]); From a81ffafaaf11b3d44f074d7cae287ac4104f1ff4 Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Thu, 31 Oct 2024 17:51:56 +0100 Subject: [PATCH 09/15] Mass index determination for N2HDM in AdjustRotationMatrix() --- include/BSMPT/models/ClassPotentialN2HDM.h | 1 + src/models/ClassPotentialN2HDM.cpp | 59 ++++++++++++---------- src/models/ClassPotentialR2HDM.cpp | 2 +- 3 files changed, 33 insertions(+), 29 deletions(-) diff --git a/include/BSMPT/models/ClassPotentialN2HDM.h b/include/BSMPT/models/ClassPotentialN2HDM.h index 96917ed72..da629186a 100644 --- a/include/BSMPT/models/ClassPotentialN2HDM.h +++ b/include/BSMPT/models/ClassPotentialN2HDM.h @@ -101,6 +101,7 @@ class Class_Potential_N2HDM : public Class_Potential_Origin 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; diff --git a/src/models/ClassPotentialN2HDM.cpp b/src/models/ClassPotentialN2HDM.cpp index 4d8ad31c6..6173609ba 100644 --- a/src/models/ClassPotentialN2HDM.cpp +++ b/src/models/ClassPotentialN2HDM.cpp @@ -911,6 +911,37 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() } } + // Determine the additional indices relative to the SM Higgs + 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 MSM = 125.09; + double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h1]) - MSM); + double diff2 = std::abs(std::sqrt(HiggsMasses[pos_h2]) - MSM); + double diff3 = std::abs(std::sqrt(HiggsMasses[pos_h3]) - MSM); + 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_SM = pos_h2; + pos_h_l = pos_h1; + pos_h_H = pos_h3; + } + else + { + pos_h_SM = pos_h3; + pos_h_l = pos_h1; + pos_h_H = pos_h2; + } + MatrixXd HiggsRotFixed(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { @@ -1040,34 +1071,6 @@ void Class_Potential_N2HDM::TripleHiggsCouplings() } } - int pos_h_SM = -1, pos_h_l = -1, pos_h_H = -1; - - std::vector HiggsMasses; - HiggsMasses = HiggsMassesSquared(vevTree, 0); - - double MSM = 125.09; - double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h1]) - MSM); - double diff2 = std::abs(std::sqrt(HiggsMasses[pos_h2]) - MSM); - double diff3 = std::abs(std::sqrt(HiggsMasses[pos_h3]) - MSM); - 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_SM = pos_h2; - pos_h_l = pos_h1; - pos_h_H = pos_h3; - } - else - { - pos_h_SM = pos_h3; - pos_h_l = pos_h1; - pos_h_H = pos_h2; - } - std::vector HiggsOrder(NHiggs); HiggsOrder[0] = pos_G1; HiggsOrder[1] = pos_G2; diff --git a/src/models/ClassPotentialR2HDM.cpp b/src/models/ClassPotentialR2HDM.cpp index e964eddf8..127f3826c 100644 --- a/src/models/ClassPotentialR2HDM.cpp +++ b/src/models/ClassPotentialR2HDM.cpp @@ -570,7 +570,7 @@ void Class_Potential_R2HDM::write() const std::vector HiggsMasses; HiggsMasses = HiggsMassesSquared(vevTree, 0); - ss << "The mass spectrum is given by :\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" From 720373f79c14a4d4a5b7353cca1825cfd3ea2c3b Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Tue, 5 Nov 2024 17:04:10 +0100 Subject: [PATCH 10/15] Do the index calculation just once for all implemented models --- include/BSMPT/models/ClassPotentialC2HDM.h | 4 + .../BSMPT/models/ClassPotentialCPintheDark.h | 3 + include/BSMPT/models/ClassPotentialCxSM.h | 4 + include/BSMPT/models/ClassPotentialR2HDM.h | 1 + src/models/ClassPotentialC2HDM.cpp | 391 ++++++++---------- src/models/ClassPotentialCPintheDark.cpp | 150 ++----- src/models/ClassPotentialCxSM.cpp | 178 +++----- src/models/ClassPotentialN2HDM.cpp | 132 ++---- src/models/ClassPotentialR2HDM.cpp | 60 ++- 9 files changed, 352 insertions(+), 571 deletions(-) diff --git a/include/BSMPT/models/ClassPotentialC2HDM.h b/include/BSMPT/models/ClassPotentialC2HDM.h index a496f6ac5..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; diff --git a/include/BSMPT/models/ClassPotentialCPintheDark.h b/include/BSMPT/models/ClassPotentialCPintheDark.h index faa12b285..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; diff --git a/include/BSMPT/models/ClassPotentialCxSM.h b/include/BSMPT/models/ClassPotentialCxSM.h index 38b05f19d..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; diff --git a/include/BSMPT/models/ClassPotentialR2HDM.h b/include/BSMPT/models/ClassPotentialR2HDM.h index 77b747c07..c80167930 100644 --- a/include/BSMPT/models/ClassPotentialR2HDM.h +++ b/include/BSMPT/models/ClassPotentialR2HDM.h @@ -91,6 +91,7 @@ class Class_Potential_R2HDM : public Class_Potential_Origin 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; diff --git a/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index e2b105e11..6ddb9b2a9 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"; @@ -2135,90 +2133,20 @@ void Class_Potential_C2HDM::write() const { 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 "; @@ -2427,6 +2355,140 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() 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 + { + // Charged submatrix + if (std::abs(HiggsRot(i, pos_rho1)) + + std::abs(HiggsRot(i, pos_rho2)) > ZeroThreshold) + { + pos_G1 = i; + } + if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_G2 = i; + } + + // Neutral submatrix + if (std::abs(HiggsRot(i, pos_psi1)) + + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) + { + pos_G0 = i; + } + } + + // 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; + } + if (std::abs(HiggsRot(i, pos_eta1)) + + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) + { + pos_H2 = i; + } + + // Neutral CP-mixed submatrix + 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; + } + } + } + + // 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) + { + 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_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; + } + } + + // 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): @@ -2453,31 +2515,6 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() // => Look for row which has psi1 and psi2 mixing components =/= 0, // and the rest = 0 - int row_G0 = -1; - for (std::size_t i = 0; i < 3; i++) - { - double psi1psi2 = std::abs(HiggsRot(i, pos_psi1)) - + std::abs(HiggsRot(i, pos_psi2)); - double non_psi1psi2 = std::abs(HiggsRot(i, pos_rho1)) - + std::abs(HiggsRot(i, pos_eta1)) - + std::abs(HiggsRot(i, pos_rho2)) - + std::abs(HiggsRot(i, pos_eta2)) - + std::abs(HiggsRot(i, pos_zeta1)) - + std::abs(HiggsRot(i, pos_zeta2)); - - if (psi1psi2 > ZeroThreshold and non_psi1psi2 < ZeroThreshold) - { - row_G0 = i; - break; - } - } - - if (row_G0 == -1) - { - throw std::runtime_error("Error. Something went wrong with finding " - "the row for the neutral Goldstone."); - } - // 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; @@ -2493,12 +2530,12 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() // for the following manipulations MatrixXd MoveGoldstoneFirst(NHiggs, NHiggs); MoveGoldstoneFirst.setIdentity(NHiggs, NHiggs); - if (row_G0 != 0) + if (pos_G0 != 0) { MoveGoldstoneFirst(0, 0) = 0.; - MoveGoldstoneFirst(row_G0, row_G0) = 0.; - MoveGoldstoneFirst(0, row_G0) = 1.; - MoveGoldstoneFirst(row_G0, 0) = 1.; + 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 @@ -2595,7 +2632,7 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() } // Check if all other elements of rotation matrix are zero - bool zero_element = false; + zero_element = false; // Start with i, j = 1, skip neutral Goldstone for (std::size_t i = 1; i < NHiggs; i++) { @@ -2731,24 +2768,6 @@ void Class_Potential_C2HDM::TripleHiggsCouplings() if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; - 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); - } - } - } - MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { @@ -2757,100 +2776,50 @@ void Class_Potential_C2HDM::TripleHiggsCouplings() HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[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++) + + std::vector HiggsOrder(NHiggs); + HiggsOrder[0] = pos_G1; + HiggsOrder[1] = pos_G2; + HiggsOrder[2] = pos_H1; + HiggsOrder[3] = pos_H2; + HiggsOrder[4] = pos_G0; + if (UseHsmNotationInTripleHiggs) { - 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; + HiggsOrder[5] = pos_h_SM; + HiggsOrder[6] = pos_h_l; + HiggsOrder[7] = pos_h_H; } - for (std::size_t i = 3; i < NHiggs; i++) + else { - 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.at(countposN) = i; - countposN++; - } + HiggsOrder[5] = pos_h1; + HiggsOrder[6] = pos_h2; + HiggsOrder[7] = pos_h3; } - if (UseHsmNotationInTripleHiggs) + MatrixXd HiggsRotSort(NHiggs, NHiggs); + for (std::size_t i = 0; i < NHiggs; i++) { - std::vector HiggsMasses; - double MhUp = 0, MhDown = 0, MSM = 0; - HiggsMasses = HiggsMassesSquared(vevTree, 0); - - double NeutralHiggs[3]; - for (int i = 0; i < 3; i++) - { - NeutralHiggs[i] = HiggsMasses[posN.at(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; - } + HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); } - 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]; - + 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++) { - HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[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); + } + } } - PosSM = 5; - TripleHiggsCorrectionsCWPhysical.resize(NHiggs); TripleHiggsCorrectionsTreePhysical.resize(NHiggs); TripleHiggsCorrectionsCTPhysical.resize(NHiggs); diff --git a/src/models/ClassPotentialCPintheDark.cpp b/src/models/ClassPotentialCPintheDark.cpp index 8c7b49499..744cdb69e 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()); } @@ -1346,9 +1294,9 @@ void Class_Potential_CPintheDark::AdjustRotationMatrix() pos_zeta1 = 4, pos_psi1 = 5, pos_zeta2 = 6, pos_psi2 = 7, pos_rhoS = 8; - // Indices of mass eigenstates for rotation from semi-interaction to mass basis - int pos_Gp = -1, pos_Gm = -1, pos_Hp = -1, pos_Hm = -1, pos_HSM = -1; - int pos_G0 = -1, pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; + // 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; // basis = {rho1, eta1, rho2, eta2, zeta1, psi1, zeta2, psi2, rhoS} @@ -1530,17 +1478,6 @@ void Class_Potential_CPintheDark::TripleHiggsCouplings() if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; - // 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; - MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { @@ -1550,60 +1487,23 @@ void Class_Potential_CPintheDark::TripleHiggsCouplings() } } - const double ZeroThreshold = 1e-5; + 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; - for (size_t i = 0; i < NHiggs; i++) + MatrixXd HiggsRotSort(NHiggs, NHiggs); + for (std::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; - - // the neutral dark scalars mix - if ((std::abs(HiggsRot(i, 6)) + std::abs(HiggsRot(i, 7)) + - std::abs(HiggsRot(i, 8))) > ZeroThreshold) - { - // use that scalars are sorted by mass - if (posh1 == 0) - { - posh1 = i; - } - else - { - if (posh2 == 0) - { - posh2 = i; - } - else - { - posh3 = i; - } - } - } + HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); } - // new rotation matrix with - MatrixXd HiggsRotSort(NHiggs, NHiggs); - - 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); - std::vector TripleDeriv; TripleDeriv = WeinbergThirdDerivative(); std::vector>> GaugeBasis( diff --git a/src/models/ClassPotentialCxSM.cpp b/src/models/ClassPotentialCxSM.cpp index f554f0a29..272e76fe2 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]; - } - } - - 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; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } - 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 "; @@ -793,12 +736,13 @@ void Class_CxSM::AdjustRotationMatrix() // CxSM interaction basis // 0 1 2 3 4 5 // Gp, Gm, G0, zeta1, zeta2, zeta3 - int pos_Gp = 0, pos_Gm = 1, pos_G0 = 2, - pos_zeta1 = 3, pos_zeta2 = 4, pos_zeta3 = 5; + pos_Gp = 0, pos_Gm = 1, pos_G0 = 2; + + int pos_zeta1 = 3, pos_zeta2 = 4, pos_zeta3 = 5; // Indices of mass eigenstates for rotation from interaction to mass basis // Only neutral/physical part necessary, as Goldstones are already diagonal - int pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; + pos_h1 = -1, pos_h2 = -1, pos_h3 = -1; // Starting from 3 to NHiggs, i.e. fixing the physical Higgs bosons for (std::size_t i = 3; i < NHiggs; i++) @@ -851,6 +795,40 @@ void Class_CxSM::AdjustRotationMatrix() } } + // 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; + } + MatrixXd HiggsRotFixed(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { @@ -933,11 +911,6 @@ void Class_CxSM::TripleHiggsCouplings() if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; - 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) - MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { @@ -947,51 +920,21 @@ 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; + 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++) { - // 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) - { - // use that scalars are sorted by mass - if (posH1 == 0) - { - posH1 = i; - } - else - { - if (posH2 == 0) - { - posH2 = i; - } - else - { - posH3 = i; - } - } - } + HiggsRotSort.row(i) = HiggsRot.row(HiggsOrder[i]); } - // 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 TripleDeriv; TripleDeriv = WeinbergThirdDerivative(); std::vector>> GaugeBasis( @@ -1009,13 +952,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); diff --git a/src/models/ClassPotentialN2HDM.cpp b/src/models/ClassPotentialN2HDM.cpp index 6173609ba..6b45b922f 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; @@ -911,7 +848,8 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() } } - // Determine the additional indices relative to the SM Higgs + // 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; @@ -919,10 +857,12 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() // Due to the masses being ordered, we will always have // HiggsMasses[pos_h1] <= HiggsMasses[pos_h2] <= HiggsMasses[pos_h3] - double MSM = 125.09; - double diff1 = std::abs(std::sqrt(HiggsMasses[pos_h1]) - MSM); - double diff2 = std::abs(std::sqrt(HiggsMasses[pos_h2]) - MSM); - double diff3 = std::abs(std::sqrt(HiggsMasses[pos_h3]) - MSM); + 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; @@ -931,15 +871,15 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() } else if (diff2 < diff1 and diff2 < diff3) { - pos_h_SM = pos_h2; pos_h_l = pos_h1; + pos_h_SM = pos_h2; pos_h_H = pos_h3; } else { - pos_h_SM = pos_h3; pos_h_l = pos_h1; pos_h_H = pos_h2; + pos_h_SM = pos_h3; } MatrixXd HiggsRotFixed(NHiggs, NHiggs); @@ -1045,23 +985,6 @@ void Class_Potential_N2HDM::TripleHiggsCouplings() if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; - 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); - } - } - } - MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { @@ -1088,6 +1011,23 @@ void Class_Potential_N2HDM::TripleHiggsCouplings() 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/ClassPotentialR2HDM.cpp b/src/models/ClassPotentialR2HDM.cpp index 127f3826c..9dc7e710c 100644 --- a/src/models/ClassPotentialR2HDM.cpp +++ b/src/models/ClassPotentialR2HDM.cpp @@ -563,7 +563,7 @@ void Class_Potential_R2HDM::write() const { for (std::size_t j = 0; j < NHiggs; j++) { - HiggsRot(i, j) = HiggsRotationMatrix[i][j]; + HiggsRot(i, j) = HiggsRotationMatrixEnsuredConvention[i][j]; } } @@ -837,6 +837,30 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() } } + // 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++) { @@ -911,23 +935,6 @@ void Class_Potential_R2HDM::TripleHiggsCouplings() if (CalculatedTripleCopulings) return; CalculatedTripleCopulings = true; - 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); - } - } - } - MatrixXd HiggsRot(NHiggs, NHiggs); for (std::size_t i = 0; i < NHiggs; i++) { @@ -953,6 +960,23 @@ void Class_Potential_R2HDM::TripleHiggsCouplings() 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); From 1ef271ea17cce9f63f241727ac7531ed05eecd91 Mon Sep 17 00:00:00 2001 From: Christoph Borschensky Date: Wed, 6 Nov 2024 10:49:57 +0100 Subject: [PATCH 11/15] Added a check in case a non-Goldstone appears in a row of the mixing matrices where only massless Goldstones should be (which can happen for a faulty parameter point with a negative mass) --- src/models/ClassPotentialC2HDM.cpp | 26 ++++++++++-------- src/models/ClassPotentialCPintheDark.cpp | 3 +- src/models/ClassPotentialCxSM.cpp | 35 ++++++++++++++++++++++-- src/models/ClassPotentialN2HDM.cpp | 19 +++++++------ src/models/ClassPotentialR2HDM.cpp | 6 ++-- 5 files changed, 63 insertions(+), 26 deletions(-) diff --git a/src/models/ClassPotentialC2HDM.cpp b/src/models/ClassPotentialC2HDM.cpp index 6ddb9b2a9..3c1b29b70 100644 --- a/src/models/ClassPotentialC2HDM.cpp +++ b/src/models/ClassPotentialC2HDM.cpp @@ -2371,18 +2371,23 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() { pos_G1 = i; } - if (std::abs(HiggsRot(i, pos_eta1)) + else if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) { pos_G2 = i; } - // Neutral submatrix - if (std::abs(HiggsRot(i, pos_psi1)) + 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 @@ -2395,16 +2400,16 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() { pos_H1 = i; } - if (std::abs(HiggsRot(i, pos_eta1)) + else if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) { pos_H2 = i; } - // Neutral CP-mixed submatrix - 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) + 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) { @@ -2592,7 +2597,7 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() pos_si_H1 = i; } } - if (std::abs(RotGoldstoneMassBasis(i, pos_si_eta1)) + else if (std::abs(RotGoldstoneMassBasis(i, pos_si_eta1)) + std::abs(RotGoldstoneMassBasis(i, pos_si_eta2)) > ZeroThreshold) // use that 0 = mGpm < mHpm { @@ -2605,10 +2610,9 @@ void Class_Potential_C2HDM::AdjustRotationMatrix() pos_si_H2 = i; } } - // Neutral submatrix (mixed CP-even and CP-odd states); // neutral Goldstone already rotated out - if (std::abs(RotGoldstoneMassBasis(i, pos_si_zeta1)) + 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 diff --git a/src/models/ClassPotentialCPintheDark.cpp b/src/models/ClassPotentialCPintheDark.cpp index 744cdb69e..f3eb24d07 100644 --- a/src/models/ClassPotentialCPintheDark.cpp +++ b/src/models/ClassPotentialCPintheDark.cpp @@ -1327,9 +1327,8 @@ void Class_Potential_CPintheDark::AdjustRotationMatrix() { pos_G0 = i; } - // the neutral dark scalars mix - if ((std::abs(HiggsRot(i, pos_zeta2)) + std::abs(HiggsRot(i, pos_psi2)) + + 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 diff --git a/src/models/ClassPotentialCxSM.cpp b/src/models/ClassPotentialCxSM.cpp index 272e76fe2..fa907cb64 100644 --- a/src/models/ClassPotentialCxSM.cpp +++ b/src/models/ClassPotentialCxSM.cpp @@ -736,14 +736,42 @@ void Class_CxSM::AdjustRotationMatrix() // CxSM interaction basis // 0 1 2 3 4 5 // Gp, Gm, G0, zeta1, zeta2, zeta3 - pos_Gp = 0, pos_Gm = 1, pos_G0 = 2; - + int pos_i_G1 = 0, pos_i_G2 = 1, pos_i_G0 = 2; int pos_zeta1 = 3, pos_zeta2 = 4, pos_zeta3 = 5; // 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++) + { + if (std::abs(HiggsRot(i, pos_i_G1)) + + std::abs(HiggsRot(i, pos_i_G2)) > ZeroThreshold) + { + if (pos_Gp == -1) + { + pos_Gp = i; + } + else + { + 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 @@ -802,6 +830,9 @@ void Class_CxSM::AdjustRotationMatrix() 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]) diff --git a/src/models/ClassPotentialN2HDM.cpp b/src/models/ClassPotentialN2HDM.cpp index 6b45b922f..835dea12d 100644 --- a/src/models/ClassPotentialN2HDM.cpp +++ b/src/models/ClassPotentialN2HDM.cpp @@ -759,18 +759,23 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() { pos_G1 = i; } - if (std::abs(HiggsRot(i, pos_eta1)) + else if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) { pos_G2 = i; } - // Neutral CP-odd submatrix - if (std::abs(HiggsRot(i, pos_psi1)) + 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 @@ -783,21 +788,19 @@ void Class_Potential_N2HDM::AdjustRotationMatrix() { pos_H1 = i; } - if (std::abs(HiggsRot(i, pos_eta1)) + else if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) { pos_H2 = i; } - // Neutral CP-odd submatrix - if (std::abs(HiggsRot(i, pos_psi1)) + else if (std::abs(HiggsRot(i, pos_psi1)) + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) { pos_A = i; } - // Neutral CP-even submatrix - if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) + 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 { diff --git a/src/models/ClassPotentialR2HDM.cpp b/src/models/ClassPotentialR2HDM.cpp index 9dc7e710c..285d5b6e8 100644 --- a/src/models/ClassPotentialR2HDM.cpp +++ b/src/models/ClassPotentialR2HDM.cpp @@ -765,7 +765,7 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() pos_H1 = i; } } - if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > + else if (std::abs(HiggsRot(i, pos_eta1)) + std::abs(HiggsRot(i, pos_eta2)) > ZeroThreshold) // use that mGpm < mHpm { if (pos_G2 == -1) @@ -777,7 +777,7 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() pos_H2 = i; } } - if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) > + else if (std::abs(HiggsRot(i, pos_zeta1)) + std::abs(HiggsRot(i, pos_zeta2)) > ZeroThreshold) // use that mh < mH { if (pos_h == -1) @@ -789,7 +789,7 @@ void Class_Potential_R2HDM::AdjustRotationMatrix() pos_H = i; } } - if (std::abs(HiggsRot(i, pos_psi1)) + std::abs(HiggsRot(i, pos_psi2)) > + else if (std::abs(HiggsRot(i, pos_psi1)) + std::abs(HiggsRot(i, pos_psi2)) > ZeroThreshold) // use that 0 = mG0 < mA { if (pos_G0 == -1) From 3a7f73c96de26336fee40c3d07c952168f0f6667 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 14 Nov 2024 08:41:45 +0100 Subject: [PATCH 12/15] Bump JamesIves/github-pages-deploy-action from 4.6.4 to 4.6.9 (#176) Bumps [JamesIves/github-pages-deploy-action](https://github.com/jamesives/github-pages-deploy-action) from 4.6.4 to 4.6.9. - [Release notes](https://github.com/jamesives/github-pages-deploy-action/releases) - [Commits](https://github.com/jamesives/github-pages-deploy-action/compare/v4.6.4...v4.6.9) --- updated-dependencies: - dependency-name: JamesIves/github-pages-deploy-action dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/doc.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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. From 7e1583adb4e7a4451e982b2699c5daeafa68a13c Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 19 Nov 2024 00:35:10 +0100 Subject: [PATCH 13/15] Bump codecov/codecov-action from 4 to 5 (#177) Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 4 to 5. - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/v4...v5) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/codecov.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 8e6943164e699a6f867db5bdcfe4ceb1650db13d Mon Sep 17 00:00:00 2001 From: Philipp Basler <28863303+phbasler@users.noreply.github.com> Date: Mon, 9 Dec 2024 11:04:28 +0100 Subject: [PATCH 14/15] Combine cmake and c++ linter in one job (#185) * Combine linters * Add whitelines to trigger formatting * Automatically applied linter * add linter PAT * Automatically applied linter * Remove new lines * rename linter job --------- Co-authored-by: GitHub Actions Bot <> --- .github/workflows/cmake-format-linter.yml | 29 ------------------- .../workflows/{cpp-linter.yml => linter.yml} | 24 ++++++++++----- 2 files changed, 16 insertions(+), 37 deletions(-) delete mode 100644 .github/workflows/cmake-format-linter.yml rename .github/workflows/{cpp-linter.yml => linter.yml} (80%) 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/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: | From 25acde2620384dca623ea93296de3ac8420c4623 Mon Sep 17 00:00:00 2001 From: Borschensky Christoph Date: Wed, 11 Dec 2024 00:57:09 +0100 Subject: [PATCH 15/15] Moved all almost_the_same functions to utility.cpp/.h --- include/BSMPT/minimum_tracer/minimum_tracer.h | 20 +------ include/BSMPT/models/ClassPotentialOrigin.h | 7 --- include/BSMPT/utility/utility.h | 23 ++++++++ src/minimum_tracer/minimum_tracer.cpp | 46 --------------- src/models/ClassPotentialOrigin.cpp | 20 ------- src/utility/utility.cpp | 59 +++++++++++++++++++ 6 files changed, 83 insertions(+), 92 deletions(-) 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/ClassPotentialOrigin.h b/include/BSMPT/models/ClassPotentialOrigin.h index 6193bdc0a..9bff3b634 100644 --- a/include/BSMPT/models/ClassPotentialOrigin.h +++ b/include/BSMPT/models/ClassPotentialOrigin.h @@ -978,13 +978,6 @@ class Class_Potential_Origin * Ensures the correct rotation matrix convention */ virtual void AdjustRotationMatrix() = 0; - /** - * Returns true if two values are the same given some relative precision - */ - bool almost_the_same(double a, double b, double rel_precision = 0.01); - bool almost_the_same(std::complex a, - std::complex b, - double rel_precision = 0.01); /** * Checks whether rotation matrix is properly set after implying conventions */ 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/ClassPotentialOrigin.cpp b/src/models/ClassPotentialOrigin.cpp index 7fde4749a..6d1ae7947 100644 --- a/src/models/ClassPotentialOrigin.cpp +++ b/src/models/ClassPotentialOrigin.cpp @@ -959,26 +959,6 @@ Class_Potential_Origin::SecondDerivativeOfEigenvaluesNonRepeated( return res; } -bool Class_Potential_Origin::almost_the_same(double a, - double b, - double rel_precision) -{ - if (std::abs(a) < 1e-10 and std::abs(b) < 1e-10) - { - return true; - } - return std::abs(a - b) < std::abs(a + b) / 2 * rel_precision; -} - -bool Class_Potential_Origin::almost_the_same(std::complex a, - std::complex b, - double rel_precision) -{ - bool real_part = almost_the_same(a.real(), b.real(), rel_precision); - bool imag_part = almost_the_same(a.imag(), b.imag(), rel_precision); - return (real_part and imag_part); -} - // 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 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