a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
SUSYSpectrum.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 #include <iostream>
9 #include <stdexcept>
10 #include <limits>
11 #include "SUSYSpectrum.h"
12 #include "SUSY.h"
13 
14 
15 SUSYSpectrum::SUSYSpectrum(const SUSY & SUSY_in)
16 : mySUSY(SUSY_in), Mchargino(2,2,0.), Mneutralino(4,4,0.),
17  U(2,2,0.), V(2,2,0.), N(4,4,0.),
18  Msup2(6,0.), Msdown2(6,0.), Msneutrino2(6,0.), Mselectron2(6,0.),
19  mch(2,0.), mneu(4,0.), m_su2(6,0.), m_sd2(6,0.), m_sn2(6,0.), m_se2(6,0.),
20  Ru(6,6,0.), Rd(6,6,0.), Rn(6,6,0.), Rl(6,6,0.)
21 {
22 }
23 
24 bool SUSYSpectrum::CalcHiggs(double mh[4], gslpp::complex& saeff_i)
25 {
26  double Mw = mySUSY.Mw_tree();
27  double Mz = mySUSY.getMz();
28  double cos2b = 2.0 * mySUSY.getCosb() * mySUSY.getCosb() - 1.0;
29 
30  /* charged Higgs mass */
31  mh[3] = mySUSY.mHptree;
32 
33  /* pseudo-scalar Higgs mass */
34  mh[2] = sqrt(mh[3] * mh[3] - Mw * Mw);
35 
36  double temp = mh[2] * mh[2] + Mz * Mz;
37  double temp1 = 2.0 * mh[2] * Mz * cos2b;
38  double temp2 = sqrt(temp * temp - temp1 * temp1);
39 
40  /* light Higgs mass */
41  mh[0] = sqrt((temp - temp2)/2.0);
42 
43  /* heavy Higgs mass */
44  mh[1] = sqrt((temp + temp2)/2.0);
45 
46  /* CP-even Higgs mixing angle*/
47  saeff_i = sin(atan((mh[2]*mh[2]+Mz*Mz)*sqrt(1-cos2b*cos2b)/((mh[2]*mh[2]-Mz*Mz)*cos2b))/2.0);
48 
49  return true;
50 }
51 
52 bool SUSYSpectrum::CalcChargino(gslpp::matrix<gslpp::complex>& U_i, gslpp::matrix<gslpp::complex>& V_i, gslpp::vector<double>& mch_i)
53 {
54  double Mw = mySUSY.Mw_tree();
55 
56  Mchargino.assign(0, 0, mySUSY.getM2());
57  Mchargino.assign(0, 1, sqrt(2) * Mw * mySUSY.getSinb());
58  //
59  Mchargino.assign(1, 0, sqrt(2) * Mw * mySUSY.getCosb());
60  Mchargino.assign(1, 1, mySUSY.getMuH());
61 
62  /*
63  * M.singularvalue(&U, &V, &S) decompose M into
64  * M = U diag(S) V^+.
65  */
66  gslpp::matrix<gslpp::complex> Utmp(2, 2, 0.), Vtmp(2, 2, 0.);
67  Mchargino.singularvalue(Utmp, Vtmp, mch_i);
68 
72  U_i = Utmp.transpose();
73  V_i = Vtmp.hconjugate();
74 
75  return true;
76 }
77 
78 bool SUSYSpectrum::CalcNeutralino(gslpp::matrix<gslpp::complex>& N_i, gslpp::vector<double>& mneu_i)
79 {
80  double Mw = mySUSY.Mw_tree();
81  double Mz = mySUSY.getMz();
82  double cW2 = Mw*Mw/Mz/Mz;
83  double cW = sqrt(cW2);
84  double sW = sqrt(1.0 - cW2);
85  double sb = mySUSY.getSinb();
86  double cb = mySUSY.getCosb();
87 
88  Mneutralino.assign(0, 0, mySUSY.getM1());
89  Mneutralino.assign(0, 2, - Mz * sW * cb);
90  Mneutralino.assign(0, 3, Mz * sW * sb);
91  //
92  Mneutralino.assign(1, 1, mySUSY.getM2());
93  Mneutralino.assign(1, 2, Mz * cW * cb);
94  Mneutralino.assign(1, 3, - Mz * cW * sb);
95  //
96  Mneutralino.assign(2, 0, Mneutralino(0,2));
97  Mneutralino.assign(2, 1, Mneutralino(1,2));
98  Mneutralino.assign(2, 3, - mySUSY.getMuH());
99  //
100  Mneutralino.assign(3, 0, Mneutralino(0,3));
101  Mneutralino.assign(3, 1, Mneutralino(1,3));
102  Mneutralino.assign(3, 2, Mneutralino(2,3));
103 
104  /*
105  * M.singularvalue(&U, &V, &S) decompose M into
106  * M = U diag(S) V^+.
107  */
108  gslpp::matrix<gslpp::complex> Ntemp1(4, 4, 0.), Ntemp2(4, 4, 0.);
109  Mneutralino.singularvalue(Ntemp1, Ntemp2, mneu_i);
110 
111  //gslpp::matrix<gslpp::complex> Nleft = Ntemp1.transpose();
112  gslpp::matrix<gslpp::complex> Nright = Ntemp2.hconjugate();
113  //std::cout << "Nleft = " << Nleft << std::endl;
114  //std::cout << "Nright = " << Nright << std::endl;
115 
116  /* adopt N=Nright as N^* M N^+ = diag, not as N M N^+.
117  * As a result, a phase rotation is required for N. */
118  gslpp::matrix<gslpp::complex> Mdiag_tmp(4, 4, 0.);
119  Mdiag_tmp = Nright.hconjugate().transpose() * Mneutralino * Nright.hconjugate();
120  //std::cout << "Mdiag_tmp = " << Mdiag_tmp << std::endl;
121  gslpp::vector<gslpp::complex> v1(4, 0.);
122  for(int i = 0; i < 4; i++)
123  v1.assign(i, gslpp::complex(1., Mdiag_tmp(i,i).arg()/2.0, true));
124  Nright = gslpp::matrix<gslpp::complex>(v1) * Nright;
125 
126  N_i = Nright;
127 
128  return true;
129 }
130 
131 bool SUSYSpectrum::CalcSup(gslpp::matrix<gslpp::complex>& Ru_i, gslpp::vector<double>& m_su2_i)
132 {
133  double Mw = mySUSY.Mw_tree();
134  double Mz2 = mySUSY.getMz()*mySUSY.getMz();
135  double sW2 = 1.0 - Mw*Mw/Mz2;
136  double cos2b = 2.0 * mySUSY.getCosb() * mySUSY.getCosb() - 1.0;
137  gslpp::matrix<gslpp::complex> CKM(mySUSY.getVCKM());
138  gslpp::matrix<gslpp::complex> Id3 = gslpp::matrix<gslpp::complex>::Id(3);
139 
140  /* DRbar up-type quark masses at scale Q */
141  gslpp::matrix<double> Mu(3, 3, 0.);
142  Mu(0, 0) = mySUSY.Mq_Q(mySUSY.UP);
143  Mu(1, 1) = mySUSY.Mq_Q(mySUSY.CHARM);
144  Mu(2, 2) = mySUSY.Mq_Q(mySUSY.TOP);
145 
146  gslpp::matrix<gslpp::complex> uLL( CKM * mySUSY.msQhat2 * CKM.hconjugate() + Mu * Mu
147  + cos2b * Mz2 * (1.0/2.0 - 2.0/3.0 * sW2) * Id3 );
148  gslpp::matrix<gslpp::complex> uLR( mySUSY.v2()/sqrt(2.0) * mySUSY.getTUhat().hconjugate()
149  - mySUSY.getMuH() * Mu / mySUSY.getTanb() );
150  gslpp::matrix<gslpp::complex> uRR( mySUSY.msUhat2 + Mu * Mu + cos2b * Mz2 * 2.0/3.0 * sW2 * Id3 );
151  for(int i = 0; i < 3; i++)
152  for(int j = 0; j < 3; j++) {
153  Msup2.assign(i, j, uLL(i,j));
154  Msup2.assign(i, j+3, uLR(i,j));
155  Msup2.assign(i+3, j, uLR(j,i).conjugate());
156  Msup2.assign(i+3, j+3, uRR(i,j));
157  }
158 
159  /*
160  * M.singularvalue(&U, &V, &S) decompose M into
161  * M = U diag(S) V^+.
162  */
163  gslpp::matrix<gslpp::complex> RuTmp(6, 6, 0.);
164  Msup2.eigensystem(RuTmp, m_su2_i);
165 
166  Ru_i = RuTmp.hconjugate();
167 
168  return true;
169 
170 }
171 
172 bool SUSYSpectrum::CalcSdown(gslpp::matrix<gslpp::complex>& Rd_i, gslpp::vector<double>& m_sd2_i)
173 {
174  double Mw = mySUSY.Mw_tree();
175  double Mz2 = mySUSY.getMz()*mySUSY.getMz();
176  double sW2 = 1.0 - Mw*Mw/Mz2;
177  double cos2b = 2.0 * mySUSY.getCosb() * mySUSY.getCosb() - 1.0;
178  gslpp::matrix<gslpp::complex> Id3 = gslpp::matrix<gslpp::complex>::Id(3);
179 
180  /* DRbar down-type quark masses at scale Q */
181  gslpp::matrix<double> Md(3, 3, 0.);
182  Md(0, 0) = mySUSY.Mq_Q(mySUSY.DOWN);
183  Md(1, 1) = mySUSY.Mq_Q(mySUSY.STRANGE);
184  Md(2, 2) = mySUSY.Mq_Q(mySUSY.BOTTOM);
185 
186  gslpp::matrix<gslpp::complex> dLL( mySUSY.msQhat2 + Md * Md
187  + cos2b * Mz2 * (- 1.0/2.0 + 1.0/3.0 * sW2) * Id3 );
188  gslpp::matrix<gslpp::complex> dLR( mySUSY.v1()/sqrt(2.0) * mySUSY.getTDhat().hconjugate()
189  - mySUSY.getMuH() * Md * mySUSY.getTanb() );
190  gslpp::matrix<gslpp::complex> dRR( mySUSY.msDhat2 + Md * Md - cos2b * Mz2 /3.0 * sW2 * Id3 );
191  for(int i = 0; i < 3; i++)
192  for(int j = 0; j < 3; j++) {
193  Msdown2.assign(i, j, dLL(i,j));
194  Msdown2.assign(i, j+3, dLR(i,j));
195  Msdown2.assign(i+3, j, dLR(j,i).conjugate());
196  Msdown2.assign(i+3, j+3, dRR(i,j));
197  }
198 
199  /*
200  * M.singularvalue(&U, &V, &S) decompose M into
201  * M = U diag(S) V^+.
202  */
203  gslpp::matrix<gslpp::complex> RdTmp(6, 6, 0.);
204  Msdown2.eigensystem(RdTmp, m_sd2_i);
205 
206  Rd_i = RdTmp.hconjugate();
207 
208  return true;
209 
210 }
211 
212 bool SUSYSpectrum::CalcSneutrino(gslpp::matrix<gslpp::complex>& Rn_i, gslpp::vector<double>& m_sn2_i)
213 {
214  double Mz2 = mySUSY.getMz()*mySUSY.getMz();
215  double cos2b = 2.0 * mySUSY.getCosb() * mySUSY.getCosb() - 1.0;
216  gslpp::matrix<gslpp::complex> Id3 = gslpp::matrix<gslpp::complex>::Id(3);
217 
218  gslpp::matrix<gslpp::complex> nLL( mySUSY.msLhat2 + cos2b * Mz2 /2.0 * Id3 );
219  gslpp::matrix<gslpp::complex> nLR( mySUSY.v2()/sqrt(2.0) * mySUSY.getTNhat().hconjugate() );
220  gslpp::matrix<gslpp::complex> nRR( mySUSY.msNhat2 );
221  for(int i = 0; i < 3; i++) {
222  for(int j = 0; j < 3; j++) {
223  Msneutrino2.assign(i, j, nLL(i,j));
224  Msneutrino2.assign(i, j+3, nLR(i,j));
225  Msneutrino2.assign(i+3, j, nLR(j,i).conjugate());
226  Msneutrino2.assign(i+3, j+3, nRR(i,j));
227  }
228  }
229  /*
230  * M.singularvalue(&U, &V, &S) decompose M into
231  * M = U diag(S) V^+.
232  */
233  gslpp::matrix<gslpp::complex> RnTmp(6, 6, 0.);
234  Msneutrino2.eigensystem(RnTmp, m_sn2_i);
235 
236  m_sn2_i(0)=std::numeric_limits<double>::max();
237  m_sn2_i(1)=std::numeric_limits<double>::max();
238  m_sn2_i(2)=std::numeric_limits<double>::max();
239  Rn_i = RnTmp.hconjugate();
240 
241  return true;
242 
243 }
244 
245 bool SUSYSpectrum::CalcSelectron(gslpp::matrix<gslpp::complex>& Rl_i, gslpp::vector<double>& m_se2_i)
246 {
247  double Mw = mySUSY.Mw_tree();
248  double Mz2 = mySUSY.getMz()*mySUSY.getMz();
249  double sW2 = 1.0 - Mw*Mw/Mz2;
250  double cos2b = 2.0 * mySUSY.getCosb() * mySUSY.getCosb() - 1.0;
251  gslpp::matrix<gslpp::complex> Id3 = gslpp::matrix<gslpp::complex>::Id(3);
252 
253  gslpp::matrix<double> Me(3, 3, 0.);
254  Me(0, 0) = mySUSY.Ml_Q(mySUSY.ELECTRON);
255  Me(1, 1) = mySUSY.Ml_Q(mySUSY.MU);
256  Me(2, 2) = mySUSY.Ml_Q(mySUSY.TAU);
257 
258  gslpp::matrix<gslpp::complex> eLL( mySUSY.msLhat2 + Me * Me
259  + cos2b * Mz2 * (- 1.0/2.0 + sW2) * Id3 );
260  gslpp::matrix<gslpp::complex> eLR( mySUSY.v1()/sqrt(2.0) * mySUSY.getTEhat().hconjugate()
261  - mySUSY.getMuH() * Me * mySUSY.getTanb() );
262  gslpp::matrix<gslpp::complex> eRR( mySUSY.msEhat2 + Me * Me - cos2b * Mz2 * sW2 * Id3 );
263 
264  for(int i = 0; i < 3; i++)
265  {
266  for(int j = 0; j < 3; j++)
267  {
268  Mselectron2.assign(i, j, eLL(i,j));
269  Mselectron2.assign(i, j+3, eLR(i,j));
270  Mselectron2.assign(i+3, j, eLR(j,i).conjugate());
271  Mselectron2.assign(i+3, j+3, eRR(i,j));
272  }
273  }
274 
275  /*
276  * M.singularvalue(&U, &V, &S) decompose M into
277  * M = U diag(S) V^+.
278  */
279  gslpp::matrix<gslpp::complex> RlTmp(6, 6, 0.);
280  Mselectron2.eigensystem(RlTmp, m_se2_i);
281 
282  Rl_i = RlTmp.hconjugate();
283 
284  return true;
285 
286 }
287 
288 void SUSYSpectrum::SortSfermionMasses(gslpp::vector<double>& m_sf2, gslpp::matrix<gslpp::complex>& Rf) const
289 {
290  int newIndex[6];
291  for (int i = 0; i < 6; i++)
292  newIndex[i] = i;
293 
294  /* sort sfermion masses in increasing order */
295  for (int i = 0; i < 5; i++)
296  for (int k = i + 1; k < 6; k++)
297  if (m_sf2(i) > m_sf2(k)) {
298  std::swap(m_sf2(i), m_sf2(k));
299  std::swap(newIndex[i], newIndex[k]);
300  }
301 
302  /* sort the corresponding rotation matrix, where the first(second) index
303  * denotes mass(gauge) eigenstates. */
304  gslpp::matrix<gslpp::complex> myRf(6, 6, 0.);
305  for (int i = 0; i < 6; i++)
306  for (int k = 0; k < 6; k++)
307  myRf.assign(k, i, Rf(newIndex[k], i));
308  Rf = myRf;
309 }
310 //
311 //bool SUSYSpectrum::CalcSpectrum()
312 //{
313 // CalcHiggs();
314 // CalcChargino();
315 // CalcNeutralino();
316 // CalcSup();
317 // CalcSdown();
319 // CalcSelectron();
320 // return true;
321 //}
QCD::TAU
Definition: QCD.h:316
SUSYSpectrum::CalcSdown
bool CalcSdown(gslpp::matrix< gslpp::complex > &Rd_i, gslpp::vector< double > &m_sd2_i)
Computes the down-type squark spectrum at tree level.
Definition: SUSYSpectrum.cpp:172
SUSYSpectrum::CalcSneutrino
bool CalcSneutrino(gslpp::matrix< gslpp::complex > &Rn_i, gslpp::vector< double > &m_sn2_i)
Computes the sneutrino spectrum at tree level.
Definition: SUSYSpectrum.cpp:212
QCD::BOTTOM
Definition: QCD.h:329
SUSY::msNhat2
gslpp::matrix< gslpp::complex > msNhat2
Definition: SUSY.h:549
SUSYSpectrum::Msdown2
gslpp::matrix< gslpp::complex > Msdown2
Definition: SUSYSpectrum.h:325
SUSYSpectrum::Msup2
gslpp::matrix< gslpp::complex > Msup2
Stores the tree-level Up-squark, Down-squark, Sneutrino, and Slepton mass matrix.
Definition: SUSYSpectrum.h:325
SUSY::getTDhat
gslpp::matrix< gslpp::complex > getTDhat() const
Gets the trilinear-coupling matrix for down-type squarks.
Definition: SUSY.h:340
gslpp::matrix< double >
A class for constructing and defining operations on real matrices.
Definition: gslpp_matrix_double.h:48
SUSY::Ml_Q
double Ml_Q(const lepton l) const
Definition: SUSY.h:493
SUSY::msLhat2
gslpp::matrix< gslpp::complex > msLhat2
Definition: SUSY.h:549
SUSY::msQhat2
gslpp::matrix< gslpp::complex > msQhat2
Definition: SUSY.h:549
SUSYSpectrum::CalcHiggs
bool CalcHiggs(double mh[4], gslpp::complex &saeff_i)
Computes the Higgs spectrum at tree level.
Definition: SUSYSpectrum.cpp:24
SUSY::getSinb
double getSinb() const
Gets .
Definition: SUSY.h:186
gslpp::sin
complex sin(const complex &z)
Definition: gslpp_complex.cpp:420
SUSY::getM2
gslpp::complex getM2() const
Gets the wino mass.
Definition: SUSY.h:118
QCD::UP
Definition: QCD.h:324
SUSY::msDhat2
gslpp::matrix< gslpp::complex > msDhat2
Definition: SUSY.h:549
QCD::CHARM
Definition: QCD.h:326
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
gslpp::matrix< gslpp::complex >
QCD::ELECTRON
Definition: QCD.h:312
SUSY::v1
double v1() const
Definition: SUSY.cpp:356
SUSYSpectrum::SortSfermionMasses
void SortSfermionMasses(gslpp::vector< double > &m_sf2, gslpp::matrix< gslpp::complex > &Rf) const
Definition: SUSYSpectrum.cpp:288
SUSYSpectrum::Msneutrino2
gslpp::matrix< gslpp::complex > Msneutrino2
Definition: SUSYSpectrum.h:325
SUSYSpectrum::mh
double mh[4]
Stores the tree-level Higgs spectrum.
Definition: SUSYSpectrum.h:305
SUSY::mHptree
double mHptree
Definition: SUSY.h:542
SUSY
A base class for SUSY models.
Definition: SUSY.h:26
SUSYSpectrum::CalcSup
bool CalcSup(gslpp::matrix< gslpp::complex > &Ru_i, gslpp::vector< double > &m_su2_i)
Computes the up-type squark spectrum at tree level.
Definition: SUSYSpectrum.cpp:131
QCD::TOP
Definition: QCD.h:328
SUSYSpectrum::mySUSY
const SUSY & mySUSY
Definition: SUSYSpectrum.h:300
SUSYSpectrum::CalcNeutralino
bool CalcNeutralino(gslpp::matrix< gslpp::complex > &N_i, gslpp::vector< double > &mneu_i)
Computes the neutralino spectrum at tree level.
Definition: SUSYSpectrum.cpp:78
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
SUSYSpectrum.h
SUSY::v2
double v2() const
Definition: SUSY.cpp:361
SUSY::getTEhat
gslpp::matrix< gslpp::complex > getTEhat() const
Gets the trilinear-coupling matrix for charged sleptons.
Definition: SUSY.h:432
SUSY::msEhat2
gslpp::matrix< gslpp::complex > msEhat2
Definition: SUSY.h:549
SUSY::getM1
gslpp::complex getM1() const
Gets the bino mass.
Definition: SUSY.h:109
Mchargino
A class for the chargino masses.
Definition: Mchargino.h:22
StandardModel::getVCKM
gslpp::matrix< gslpp::complex > getVCKM() const
A get method to retrieve the CKM matrix.
Definition: StandardModel.h:870
SUSYSpectrum::SUSYSpectrum
SUSYSpectrum(const SUSY &SUSY_in)
A SUSYSpectrum constructor.
Definition: SUSYSpectrum.cpp:15
gslpp::vector< double >
A class for constructing and defining operations on real vectors.
Definition: gslpp_vector_double.h:33
SUSY::getTanb
double getTanb() const
Gets .
Definition: SUSY.h:154
SUSY::getMuH
gslpp::complex getMuH() const
Gets the parameter in the superpotential.
Definition: SUSY.h:145
SUSY::Mq_Q
double Mq_Q(const quark q) const
Definition: SUSY.h:477
SUSYSpectrum::Mselectron2
gslpp::matrix< gslpp::complex > Mselectron2
Definition: SUSYSpectrum.h:325
SUSYSpectrum::Mneutralino
gslpp::matrix< gslpp::complex > Mneutralino
Definition: SUSYSpectrum.h:315
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:721
StandardModel::Mw_tree
virtual double Mw_tree() const
The tree-level mass of the boson, .
Definition: StandardModel.cpp:951
Mw
An observable class for the -boson mass.
Definition: Mw.h:22
QCD::STRANGE
Definition: QCD.h:327
CKM
A class for the CKM matrix elements.
Definition: CKM.h:23
SUSY::msUhat2
gslpp::matrix< gslpp::complex > msUhat2
Definition: SUSY.h:549
SUSY.h
SUSY::getCosb
double getCosb() const
Gets .
Definition: SUSY.h:195
SUSY::getTUhat
gslpp::matrix< gslpp::complex > getTUhat() const
Gets the trilinear-coupling matrix for up-type squarks.
Definition: SUSY.h:331
QCD::DOWN
Definition: QCD.h:325
Mneutralino
A class for the neutralino masses.
Definition: Mneutralino.h:24
SUSYSpectrum::CalcChargino
bool CalcChargino(gslpp::matrix< gslpp::complex > &U_i, gslpp::matrix< gslpp::complex > &V_i, gslpp::vector< double > &mch_i)
Computes the chargino spectrum at tree level.
Definition: SUSYSpectrum.cpp:52
gslpp::vector< gslpp::complex >
SUSY::getTNhat
gslpp::matrix< gslpp::complex > getTNhat() const
Gets the trilinear-coupling matrix for sneutrinos.
Definition: SUSY.h:423
QCD::MU
Definition: QCD.h:314
SUSYSpectrum::CalcSelectron
bool CalcSelectron(gslpp::matrix< gslpp::complex > &Rl_i, gslpp::vector< double > &m_se2_i)
Computes the charged-slepton spectrum at tree level.
Definition: SUSYSpectrum.cpp:245