flavour
BGL.h
Go to the documentation of this file.
1 #ifndef BGL_H
2 #define BGL_H
3 
4 #include "widthcalc.h"
5 
6 #include "TH2F.h"
7 #include "TProfile2D.h"
8 #include "TCanvas.h"
9 #include <iostream>
10 
11 #include "Math/Polynomial.h"
12 #include "Math/Interpolator.h"
13 #include "Formulas.h"
14 
15 
16 namespace BGLmodels{
17 
18 /**
19 * @brief Gauge boson
20 */
21 class Boson {
22 public:
23 
24  Boson(): C(Matrix(),2,2,2,2){}
25 
26  ex couplingL(const Fermion& f2,const Fermion& f1) const {
27  bool quiralfilter=0;
28  if(f1.type!=f2.type) return 0;
29  if(s==sScalar){
30  if(f1.helicity!=hRight && f2.helicity!=hLeft) quiralfilter=1;
31  }
32  else if(s==sVector){
33  if(f1.helicity!=hRight && f2.helicity!=hRight) quiralfilter=1;
34  }
35 
36  if(quiralfilter) return C[f2.type][f2.isospin][f1.isospin][hLeft][f2.flavour][f1.flavour];
37  return 0;
38  }
39  ex couplingR(const Fermion& f2,const Fermion& f1) const {
40  bool quiralfilter=0;
41  if(f1.type!=f2.type) return 0;
42  if(s==sScalar){
43  if(f2.helicity!=hRight && f1.helicity!=hLeft) quiralfilter=1;
44  }
45  else if(s==sVector){
46  if(f1.helicity!=hLeft && f2.helicity!=hLeft) quiralfilter=1;
47  }
48  if(quiralfilter) return C[f2.type][f2.isospin][f1.isospin][hRight][f2.flavour][f1.flavour];
49  return 0;
50  }
51  ex couplingdaggerL(const Fermion& f2,const Fermion& f1) const {
52  if(s==sScalar) return couplingR(f1,f2).conjugate();
53  return couplingL(f1,f2).conjugate();
54  }
55  ex couplingdaggerR(const Fermion& f2,const Fermion& f1) const {
56  if(s==sScalar) return couplingL(f1,f2).conjugate();
57  return couplingR(f1,f2).conjugate();
58  }
59  ex coupling(const Fermion& f2,const Fermion& f1, ex mu){
60  if(s==sScalar) return couplingL(f2,f1)*dirac_gammaL()+couplingR(f2,f1)*dirac_gammaR();
61  else return couplingL(f2,f1)*dirac_gammaL()+couplingR(f2,f1)*dirac_gammaR();
62  }
63  void reset(){
64  C=multivector<Matrix,4>(Matrix(),2,2,2,2);
65  }
66  BSpin s;
67  ex mass;
68  multivector<Matrix,4> C;
69 };
70 
71 
72 
73 //Boson::Type Boson::dagger[2][2]={Boson::scalarright,Boson::scalarleft,Boson::vectorleft,Boson::vectorright};
74 
75 /**
76 * @brief Implementation of the BGL model
77 */
78 class BGL: public Model{
79 public:
80 
81 BGL(int genL=2,int genQ=2, int lup=0, int qup=0, int mssm=0):
82  planck(6.58211928e-25),
83  GF("G_F"),
84  MZ("M_Z"),
85  MW("M_W"),
86  Mpip("Mpip",0.1396,"M_{\\pi^+}",domain::real),
87  Mpi0("Mpi0",0.1349766,"M_{\\pi^0}",domain::real),
88  MBp("MBp",5.279,"M_{B^+}",domain::real),
89  MB0("MB0",5.2795,"M_{B^0}",domain::real),
90  MBs0("MBs0",5.3663,"M_{B_s^0}",domain::real),
91  MKp("MKp",0.493677,"MKp",domain::real),
92  MK0("MK0",0.497614,"MK0",domain::real),
93  MDp("MDp",1.86957,"MDp",domain::real),
94  MD0("MD0",1.86480,"MD0",domain::real),
95  MDsp("MDsp",1.96845,"MDsp",domain::real),
96  MDs0("MDs0",0),
97  Fpi("Fpi",0.132,"Fpi",domain::real),
98  FB("FB",0.189,"FB",domain::real),
99  FBs("FBs",0.225,"FBs",domain::real),
100  FK("FK",0.159,"FK",domain::real),
101  FD("FD",0.208,"FD",domain::real),
102  FDs("FDs",0.248,"FDs",domain::real),
103  //alpha(7.297352e-3*4*M_PI),
104  cos2(pow(MW/MZ,2)),
105  g(sqrt(GF*8/sqrt(ex(2)))*MW),
106  //g(sqrt(4*Pi*alpha/(1-cos2))),
107  tanb("tg\\beta"),
108  cp("cp"),
109  McH("M_{H^+}"),
110  MR("M_{R}"),
111  MI("M_{I}"),
112  mixes(tanb,cp, genL,genQ, lup, qup, mssm),
113  mu("\\mu"),
114  BGLtype(4,0),
115  mmmax(1000),
116  stepsize(1e-2)
117  //muwidth(planck/2.197034e-6)
118  {
119  alpha=pow(g,2)*(1-cos2)/(4*Pi);
120  replacements.append(GF==1.166371e-5);
121  replacements.append(MZ==M_MZ);
122  replacements.append(MW==M_MW);
123 
124  mixes.appendtolst(replacements);
125 
126  replacements.append(Pi==M_PI);
127  replacements.append(sqrt(ex(2))==sqrt(2));
128 
129  //cout<<pow(sqrt(2)/8*pow(g/MW,2),2)<<endl;
130  //cout<<pow(1.166,2)<<endl;
131 
132  Boson boson;
133 
134  realsymbol q3("q3");
135  ex vq3=dirac_slash(q3,4);
136  varidx jmu(mu,4,1);
137 
138  for(uint i=0;i<2;i++)
139  for(uint j=0;j<3;j++)
140  for(uint k=0;k<3;k++){
141  conjtoabs.append(conjugate(mixes.V[i][j][k])==pow(abs(mixes.V[i][j][k]),2)/mixes.V[i][j][k]);
142  }
143  /*
144  //Gamma boson
145  boson.mass=0;
146  boson.s=Boson::vector;
147 
148  boson.coupsL[0][0]=Matrix(g*sqrt(1-cos2)*0);
149  boson.coupsL[1][1]=Matrix(g*sqrt(1-cos2)*(-1));
150  boson.coupsL[2][2]=Matrix(g*sqrt(1-cos2)*ex(2)/3);
151  boson.coupsL[3][3]=Matrix(g*sqrt(1-cos2)*ex(-1)/3);
152 
153  boson.coupsR[0][0]=Matrix(g*sqrt(1-cos2)*0);
154  boson.coupsR[1][1]=Matrix(g*sqrt(1-cos2)*(-1));
155  boson.coupsR[2][2]=Matrix(g*sqrt(1-cos2)*ex(2)/3);
156  boson.coupsR[3][3]=Matrix(g*sqrt(1-cos2)*ex(-1)/3);
157 
158  bosons.push_back(boson);
159  boson.reset();
160  */
161  //W+ boson
162  boson.mass=MW;
163  boson.s=sVector;
164 
165  for(uint t=tLepton;t<=tQuark;t++) boson.C[t][iUp][iDown][hLeft]=mixes.V[t]*Matrix(g/sqrt(ex(2)));
166  Boson wboson=boson;
167  bosons.push_back(boson);
168  boson.reset();
169 
170  //H+ boson
171  boson.mass=McH;
172  boson.s=sScalar;
173 
174  for(uint t=tLepton;t<=tQuark;t++)
175  for(uint i=iUp;i<=iDown;i++) boson.C[t][iUp][iDown][i]=mixes.VN[t][i]*Matrix(g/MW/sqrt(ex(2)));
176  Boson chiggs=boson;
177  bosons.push_back(boson);
178  boson.reset();
179 
180  for(int b=bosons.size()-1;b>=0;b--){
181  boson.mass=bosons[b].mass;
182  boson.s=bosons[b].s;
183  if(boson.s==sVector)
184  for(uint t=tLepton;t<=tQuark;t++)
185  for(uint i=iUp;i<=iDown;i++)
186  for(uint j=iUp;j<=iDown;j++)
187  for(uint h=hLeft;h<=hRight;h++){
188  boson.C[t][i][j][h]=bosons[b].C[t][j][i][h].conjugate();
189  }
190  else for(uint t=tLepton;t<=tQuark;t++)
191  for(uint i=iUp;i<=iDown;i++)
192  for(uint j=iUp;j<=iDown;j++)
193  for(uint h=hLeft;h<=hRight;h++){
194  boson.C[t][i][j][hLeft]=bosons[b].C[t][j][i][hRight].conjugate();
195  boson.C[t][i][j][hRight]=bosons[b].C[t][j][i][hLeft].conjugate();
196  }
197  bosons.push_back(boson);
198  boson.reset();
199  }
200 
201  //(R+iI)/sqrt(2) boson
202  boson.mass=MR;
203  boson.s=sScalar;
204 
205  for(uint t=tLepton;t<=tQuark;t++){
206  boson.C[t][iDown][iDown][hRight]=mixes.N[t][iDown]*Matrix(g/MW/ex(2));
207  boson.C[t][iUp][iUp][hLeft]=mixes.N[t][iUp].conjugate()*Matrix(g/MW/ex(2));
208  boson.C[t][iDown][iDown][hLeft]=mixes.N[t][iDown].conjugate()*Matrix(g/MW/ex(2));
209  boson.C[t][iUp][iUp][hRight]=mixes.N[t][iUp]*Matrix(g/MW/ex(2));
210  }
211  bosons.push_back(boson);
212  boson.reset();
213 
214  //(R+iI)/sqrt(2) boson
215  boson.mass=MI;
216  boson.s=sScalar;
217 
218  for(uint t=tLepton;t<=tQuark;t++){
219  boson.C[t][iDown][iDown][hRight]=mixes.N[t][iDown]*Matrix(I*g/MW/ex(2));
220  boson.C[t][iUp][iUp][hLeft]=mixes.N[t][iUp].conjugate()*Matrix(I*g/MW/ex(2));
221  boson.C[t][iDown][iDown][hLeft]=mixes.N[t][iDown].conjugate()*Matrix(-I*g/MW/ex(2));
222  boson.C[t][iUp][iUp][hRight]=mixes.N[t][iUp]*Matrix(-I*g/MW/ex(2));
223  }
224  bosons.push_back(boson);
225  boson.reset();
226 
227  Fermion electron(tLepton,iDown,fElectron);
228  Fermion electronR(tLepton,iDown,fElectron,cParticle,hRight);
229 
230  Fermion muon(tLepton,iDown,fMuon);
231  Fermion muonR(tLepton,iDown,fMuon,cParticle,hRight);
232 
233  Fermion tau(tLepton,iDown,fTau);
234  Fermion tauR(tLepton,iDown,fTau,cParticle,hRight);
235  Fermion neutrino(tLepton,iUp);
236  Fermion neutrinotau(tLepton,iUp,fTau);
237  Fermion neutrinomuon(tLepton,iUp,fMuon);
238  Fermion neutrinoe(tLepton,iUp,fElectron);
239 
240  Fermion up(tQuark,iUp,fElectron);
241  Fermion down(tQuark,iDown,fElectron);
242  Fermion bottom(tQuark,iDown,fTau);
243  Fermion strange(tQuark,iDown,fMuon);
244  Fermion charm(tQuark,iUp,fMuon);
245  Fermion top(tQuark,iUp,fTau);
246 
247  Meson Pi0d(down,down,Mpi0,Fpi);
248  Meson Pi0u(down,down,Mpi0,Fpi);
249  Meson Pip(up,down,Mpip,Fpi);
250  Meson Pim(down,up,Mpip,Fpi);
251 
252  Meson K0(down,strange,MK0,FK);
253  Meson Kp(up,strange,MKp,FK);
254 
255  Meson D0(charm,up,MD0,FD);
256  Meson Dp(charm,down,MDp,FD);
257  Meson Dsp(charm,strange,MDsp,FDs);
258 
259  Meson B0(down,bottom,MB0,FB);
260  Meson Bp(up,bottom,MBp,FB);
261  Meson Bs0(strange,bottom,MBs0,FBs);
262 
263  lst sb;
264  //sb.append(mixes.M[tQuark][iUp][0][0]==0);
265  sb.append(pow(abs(mixes.V[0][2][2]),2)==1-pow(abs(mixes.V[0][1][2]),2)-pow(abs(mixes.V[0][0][2]),2));
266  sb.append(pow(abs(mixes.V[0][2][1]),2)==1-pow(abs(mixes.V[0][1][1]),2)-pow(abs(mixes.V[0][0][1]),2));
267 
268  //cout<<"Btaunu "<<collect_common_factors(expand(Btaunu.subs(sb).subs(conjtoabs)))<<endl;
269 
270  cout<<latex;
271 
272  ex mutoenunu=decaywidth(muon,neutrino,electron,neutrino);
273 
274  //cout<<"mutoenunu "<<mutoenunu<<endl;
275  //add("mutoenunu",decaywidth(muon,neutrino,electron,neutrino),new gaussobs(planck/2.197034e-6,0.03));
276 
277  add("muRtoeRnunu",gRR2(muon,electron),new limitedobs(std::pow(0.035,2),0.95));
278 
279  //add("tautoenunu",decaywidth(tau,neutrino,electron,neutrino),new gaussobs(planck/290.6e-15*0.1782,0.03));
280  add("tauRtoeRnunu",gRR2(tau,electron),new limitedobs(std::pow(0.7,2),0.95));
281 
282  //add("tautomununu",decaywidth(tau,neutrino,muon,neutrino),new gaussobs(planck/290.6e-15*0.1739,0.03));
283  add("tauRtomuRnunu",gRR2(tau,muon),new limitedobs(std::pow(0.72,2),0.95));
284 
285  add("tautomu_tautoe",tautomu_tautoe(),new gaussobs(1.0018,0.0014/1.0018)); //PROBLEM!!!
286  cout<<"tautomu_tautoe: "<<1/1.0018<< "ERROR: "<<0.0014/1.0018<<endl;
287  cout<<"ratio1 "<<tautomu_tautoe().subs(replacements)<<endl;
288  cout<<"ratio2 "<<(decaywidth(tau,neutrino,muon,neutrino,sVector)/decaywidth(tau,neutrino,electron,neutrino,sVector)).subs(replacements)<<endl;
289  //muto3e
290  ex mu3e=decaywidth(muon,electron,electron,electron);
291  add("muto3e", mu3e,new limitedobs(planck/2.197034e-6*1e-12));
292  cout<<"mu3e "<<decaywidthtest2(muon)<<endl;
293 
294  //tauto3e
295  add("tauto3e", decaywidth(tau,electron,electron,electron),new limitedobs(planck/290.6e-15*2.7e-8));
296  //tauto2e1mu+
297  add("tauto2e1mup", decaywidth(tau,electron,electron,muon), new limitedobs(planck/290.6e-15*1.5e-8));
298  //tauto2e1mu-
299  add("tauto2e1mu", decaywidth(tau,electron,muon,electron), new limitedobs(planck/290.6e-15*1.8e-8));
300  //tauto2mu1e+
301 
302  add("tauto2mu1ep", decaywidth(tau,muon,muon,electron), new limitedobs(planck/290.6e-15*1.7e-8));
303  cout<<"tauto2mu1ep "<<decaywidthtest2(tau)<<endl;
304  //tauto2mu1e-
305  add("tauto2mu1ep", decaywidth(tau,muon,electron,muon), new limitedobs(planck/290.6e-15*2.7e-8));
306  cout<<"tauto2mu1e "<<decaywidthtest2(tau)<<endl;
307 
308  //tauto3mu
309  add("tauto3mu", decaywidth(tau,muon,muon,muon),new limitedobs(planck/290.6e-15*2.1e-8));
310  //cout<<"tauto3mu "<<collect_common_factors(expand(decaywidth(tau,muon,muon,muon)))<<endl;
311 
312 
313  ex piratio=1.2352e-4/(mesondw(Pip,neutrino,electron,sVector)/mesondw(Pip,neutrino,muon,sVector));
314  ex picorrection=piratio.subs(replacements);
315  ex pierror=picorrection*0.0001/1.2352;
316  cout<<"PiRatio "<<picorrection-1<<" +/- "<<pierror<<endl;
317  piratio*=mesondw(Pip,neutrino,electron)/mesondw(Pip,neutrino,muon);
318  add("piontoenu_munu",piratio,new gaussobs(1.230e-4,0.003)); //PROBLEM!!!
319  cout<<"piontoenu_munu: "<<1.2352e-4/1.230e-4<<" ERROR: "<<0.003<<endl;
320 
321 
322  add("tautopinu_pitomunu",(1+0.16e-2)*fermiontomeson(tau,neutrino,Pip)/mesondw(Pip,neutrino,muon),new gaussobs((10.83e-2/290.6e-15/(0.9998770/2.6033e-8)),0.06/10.83));
323  cout<<"tautopinu/pitomunu: "<<(1+0.16e-2)*(fermiontomeson(tau,neutrino,Pip,sVector)/mesondw(Pip,neutrino,muon,sVector)).subs(replacements)/((10.83e-2/290.6e-15/(0.9998770/2.6033e-8)))<<" ERROR: "<<0.06/10.83<<endl;
324  cout<<"tautopinu: "<<fermiontomesontest(tau,neutrino,Pip)<<endl;
325  cout<<"tautopinu_pitomunu: "<<10.83e-2/290.6e-15/(0.9998770/2.6033e-8)<<" +/- "<<0.06e-2/290.6e-15/(0.9998770/2.6033e-8)<<endl;
326 
327  add("tautoKnu_Ktomunu",(1+0.9e-2)*fermiontomeson(tau,neutrino,Kp)/mesondw(Kp,neutrino,muon),new gaussobs((7e-3/290.6e-15)/(0.6355/1.238e-8),0.1/7));
328  cout<<"tautoKnu/Ktomunu: "<<(1+0.9e-2)*(fermiontomeson(tau,neutrino,Kp,sVector)/mesondw(Kp,neutrino,muon,sVector)).subs(replacements)/((7e-3/290.6e-15)/(0.6355/1.238e-8))<<" ERROR: "<<0.1/7<<endl;
329  cout<<"tautoKnu/Ktomunu: "<<(7e-3/290.6e-15)/(0.6355/1.238e-8)<<" +/- "<<(0.1e-3/290.6e-15)/(0.6355/1.238e-8)<<endl;
330 
331  //ex pi0toemu=(mesondecaywidth(Mpi0,down,down,electron,muon)+mesondecaywidth(Mpi0,up,up,electron,muon)+mesondecaywidth(Mpi0,down,down,muon,electron)+mesondecaywidth(Mpi0,up,up,muon,electron))/2;
332  ex pi0toemu=(mesondw(Pi0d,electron,muon)+mesondw(Pi0d,muon,electron)+mesondw(Pi0u,electron,muon)+mesondw(Pi0u,muon,electron))/2;
333  add("pi0toemu",pi0toemu,new limitedobs(3.6e-10*planck/8.52e-17));
334 
335  ex Kratio=2.477e-5/(mesondw(Kp,neutrino,electron,sVector)/mesondw(Kp,neutrino,muon,sVector));
336  ex Kcorrection=Kratio.subs(replacements);
337  ex Kerror=Kcorrection*0.001/2.477;
338  cout<<"KRatio "<<Kcorrection-1<<" +/- "<<Kerror<<endl;
339  Kratio*=mesondw(Kp,neutrino,electron)/mesondw(Kp,neutrino,muon);
340  add("Ktoenu_munu",Kratio,new gaussobs(2.488e-5,0.005));
341  cout<<"Ktoenu_munu: "<<2.477e-5/2.488e-5<<" ERROR: "<<0.005<<endl;
342 
343  ex k0Ltoemu=mesondw(K0,electron,muon)+mesondw(K0,muon,electron);
344  add("K0Ltoemu",k0Ltoemu,new limitedobs((4.7e-12*planck/5.116e-8)));
345  add("K0Ltoee",mesondw(K0,electron,electron),new limitedobs((9e-12*planck/5.116e-8)));
346 
347  //add("K0Ltomumu",mesondw(K0,muon,muon),new limitedobs((6.84e-9*planck/5.116e-8)));
348 
349  add("Dtoenu",mesondw(Dp,neutrino,electron),new limitedobs(8.8e-6*planck/1040e-15));
350  add("Dtomunu",mesondw(Dp,neutrino,muon),new gaussobs(3.82e-4*planck/1040e-15,0.1)); //PROBLEM!!!
351  cout<<"Dtomunu: "<<mesondw(Dp,neutrino,muon,sVector).subs(replacements)/(3.82e-4*planck/1040e-15)<<" ERROR: "<<0.1<<endl;
352 
353  add("Dtotaunu",mesondw(Dp,neutrino,tau),new limitedobs(1.2e-3*planck/1040e-15)); //PROBLEM!!!
354  cout<<"Dtotaunu: "<<mesondw(Dp,neutrino,tau,sVector).subs(replacements)/(1.2e-3*planck/1040e-15)<<" LIMIT"<<endl;
355 
356  //D0 2.6e-7/410.1e-15
357 
358  ex D0toemu=mesondw(D0,electron,muon)+mesondw(D0,muon,electron);
359  add("D0toemu",D0toemu,new limitedobs((2.6e-7*planck/410.1e-15)));
360  ex D0toee=mesondw(D0,electron,electron);
361  add("D0toee",D0toee,new limitedobs((7.9e-8*planck/410.1e-15)));
362  ex D0tomumu=mesondw(D0,muon,muon);
363  add("D0tomumu",D0tomumu,new limitedobs((1.4e-7*planck/410.1e-15)));
364 
365  //ex Dstomunu=mesondecaywidth(MDsp,strange,charm,muon,neutrino); //500e-15
366  add("Dstomunu",mesondw(Dsp,neutrino,muon),new gaussobs(5.9e-3*planck/500e-15,0.33/5.9)); //PROBLEM!!!
367  cout<<"Dstomunu: "<<mesondw(Dsp,neutrino,muon,sVector).subs(replacements)/(5.9e-3*planck/500e-15)<<" ERROR: "<<0.33/5.9<<endl;
368 
369  add("Dstoenu",mesondw(Dsp,neutrino,electron),new limitedobs(1.2e-4*planck/500e-15));
370  add("Dstotaunu",mesondw(Dsp,neutrino,tau),new gaussobs(5.43e-2*planck/500e-15,0.31/5.43)); //PROBLEM!!!
371  cout<<"Dstotaunu: "<<mesondw(Dsp,neutrino,tau,sVector).subs(replacements)/(5.43e-2*planck/500e-15)<<" ERROR: "<<0.31/5.43<<endl;
372 
373  add("Btomunu",mesondw(Bp,neutrino,muon),new limitedobs(9.8e-7*planck/1.641e-12));
374  add("Btoenu",mesondw(Bp,neutrino,electron),new limitedobs(1e-6*planck/1.641e-12));
375 
376  add("Btotaunu",mesondw(Bp,neutrino,tau),new gaussobs(1.15e-4*planck/1.641e-12,0.23/1.15));
377  //add("Btotaunu",mesondw(Bp,neutrino,tau),new gaussobs(0.79e-4*planck/1.641e-12,0.23/1.15));
378 
379  //calcuBmumu calcutest(mixes,Bs0,muon,muon,2.9e-9*planck/1.516e-12,0.7e-9*planck/1.516e-12,"Bs_to_mumu");
380  //calcuBmumu calcutest2(mixes,B0,muon,muon,3.6e-10*planck/1.519e-12,1.6e-10*planck/1.516e-19,"B_to_mumu");
381  //cout<<"TESTE "<<endl;
382  //double ps[4]={1,1e16,1e16,1e16};
383  //double resteste=0,resteste2=0;
384  //int nt=4,mt=1;
385  //calcutest.fp(&nt,ps,&mt,&resteste);
386  //calcutest2.fp(&nt,ps,&mt,&resteste2);
387  //cout<<"TESTE "<<resteste/(2.9e-9*planck/1.516e-12)<<" "<<resteste2/(3.6e-10*planck/1.519e-12)<<endl;
388  //ex B0tomumu=mesondw(B0,muon,muon);
389  //cout<<"B0tomumu "<<collect_common_factors(B0tomumu)<<endl;
390  //1.65e-4
391  //add("B0tomumu",B0tomumu,new limitedobs((8e-10*planck/1.519e-12)));
392 
393  push_back(prediction(new calcuBmumu(mixes,B0,muon,muon,new limitedobs(6.3e-10*planck/1.519e-12),"B_to_mumu")));
394  //push_back(prediction(new calcuBmumu(mixes,B0,muon,muon,new gauss2obs(3.6e-10*planck/1.519e-12,1.6e-10*planck/1.519e-12),"B_to_mumu")));
395  push_back(prediction(new calcuBmumu(mixes,Bs0,muon,muon,new gauss2obs(2.9e-9*planck/1.516e-12,0.7e-9*planck/1.516e-12),"Bs_to_mumu")));
396  push_back(prediction(new calcuBmumu(mixes,K0,muon,muon,new limitedobs(2.3e-9*planck/5.116e-8),"K0L_to_mumu")));
397 
398  cBmumu=new calcuBmumu(mixes,B0,muon,muon,new limitedobs(6.3e-10*planck/1.519e-12),"B_to_mumu");
399  cBsmumu=new calcuBmumu(mixes,Bs0,muon,muon,new gauss2obs(2.9e-9*planck/1.516e-12,0.7e-9*planck/1.516e-12),"Bs_to_mumu");
400 
401  ex B0toetau=mesondw(B0,electron,tau)+mesondw(B0,tau,electron);
402  add("B0toetau",B0toetau,new limitedobs((2.8e-5*planck/1.519e-12)));
403  ex B0tomutau=mesondw(B0,muon,tau)+mesondw(B0,tau,muon);
404  add("B0tomutau",B0tomutau,new limitedobs((2.2e-5*planck/1.519e-12)));
405  ex B0toee=mesondw(B0,electron,electron);
406  add("B0toee",B0toee,new limitedobs((8.3e-8*planck/1.519e-12)));
407  ex B0totautau=mesondw(B0,tau,tau);
408  add("B0totautau",B0totautau,new limitedobs((4.1e-3*planck/1.519e-12)));
409 
410  //B0s m=5.3663, life=1.472e-12 emu=2e-7, ee=2.8e-7 mumu=4.2e-8
411  ex Bs0toemu=mesondw(Bs0,electron,muon)+mesondw(Bs0,muon,electron);
412  add("Bs0toemu",Bs0toemu,new limitedobs((2e-7*planck/1.516e-12)));
413  ex Bs0toee=mesondw(Bs0,electron,electron);
414  add("Bs0toee",Bs0toee,new limitedobs((2.8e-7*planck/1.516e-12)));
415 // ex Bs0tomumu=mesondw(Bs0,muon,muon);
416 // add("Bs0tomumu",Bs0tomumu,new limitedobs((3.2e-9*planck/1.516e-12)));
417 
418  // add("chargedHiggs",pow(McH,-2),new limitedobs(std::pow(80.0,-2),0.9));
419 
420  cout<<"Bs0tomumu: "<<mesondwtest(Bs0,muon,muon)<<endl;
421  //add("chargedHiggs",1/McH,new limitedobs(1/80.0,0));
422 
423  /*
424  Matrix llgamma2loop=Matrix(sqrt(ex(2))*mixes.N[tQuark][iUp][2][2]*mixes.M[tQuark][iUp][fTau][fTau]*pow(1/McH*log(mixes.M[tQuark][iUp][fTau][fTau]/McH),2))*mixes.N[tLepton][iDown];
425  for(uint i=0;i<3;i++)
426  for(uint j=0;j<3;j++)
427  if(j<i) llgamma2loop[i][j]=ex(3)*pow(g*g*(1-cos2)/4/Pi/Pi,3)*llgamma2loop[j][i]*llgamma2loop[j][i].conjugate()/pow(mixes.M[tLepton][iDown][i][i],2);
428  else llgamma2loop[i][j]=0;
429  add("mutoegamma",llgamma2loop[1][0],new limitedobs(1.2e-11));
430  add("tautoegamma",llgamma2loop[2][0],new limitedobs(3.3e-8));
431  add("tautomugamma",llgamma2loop[2][1],new limitedobs(4.4e-8));
432  //add("mutoegamma",llgamma2loop[1][0],new limitedobs(planck/2.197034e-6*1.2e-11));
433  //add("tautoegamma",llgamma2loop[2][0],new limitedobs(planck/290.6e-15*3.3e-8));
434  //add("tautomugamma",llgamma2loop[2][1],new limitedobs(planck/290.6e-15*4.4e-8));
435  */
436 
437  Matrix llgammaCH;
438  //Matrix llgammaCH=Matrix(g*g/(16*Pi*4*Pi*MW*MW*McH*McH*12))*mixes.M[tLepton][iDown]*mixes.VN[tLepton][1].conjugate()*mixes.VN[tLepton][1];
439  Matrix llgammaH0M,llgammaH0E;
440  /*for(uint i=0;i<3;i++)
441  for(uint j=0;j<3;j++)
442  for(uint k=0;k<3;k++){
443  ex z=pow(fmasses[1][i][i]/McH,2); mixes.M[tQuark][iUp][i][i]/MR,2);
444  llgammaH0M[j][k]=llgammaH0M[j][k]+(mixes.VN[0][1][j][i].conjugate()*mixes.VN[0][1][k][i]+mixes.VN[0][1][i][j]*mixes.VN[0][1][i][k].conjugate())/pow(mixes.M[tLepton][iDown][i][i],2)*(2*z+6*z*z*log(z))/6;
445  llgammaH0M[j][k]=llgammaH0M[j][k]+(mixes.VN[0][1][i][j]*mixes.VN[0][1][k][i])/mixes.M[tLepton][iDown][i][i]/mixes.M[tLepton][iDown][j][j]*(3*z+2*z*log(z));
446 
447  llgammaH0E[j][k]=llgammaH0E[j][k]+(mixes.VN[0][1][j][i].conjugate()*mixes.VN[0][1][k][i]-mixes.VN[0][1][i][j]*mixes.VN[0][1][i][k].conjugate())/pow(mixes.M[tLepton][iDown][i][i],2)*(2*z+6*z*z*log(z))/6;
448  llgammaH0E[j][k]=llgammaH0E[j][k]+(mixes.VN[0][1][i][j]*mixes.VN[0][1][k][i])/mixes.M[tLepton][iDown][i][i]/mixes.M[tLepton][iDown][j][j]*(3*z+2*z*log(z));
449  }
450  */
451  llgammaH0M=Matrix(g*g/(16*Pi*4*Pi*MW*MW*McH*McH))*mixes.M[tLepton][iDown]*llgammaH0M;
452  llgammaH0E=Matrix(g*g/(16*Pi*4*Pi*MW*MW*McH*McH))*mixes.M[tLepton][iDown]*llgammaH0E;
453 
454  Matrix llgamma, llgamma2;
455 
456  for(uint i=0;i<3;i++)
457  for(uint j=0;j<3;j++){
458  //if(j<i) llgamma[i][j]=(llgammaCH[i][j]*llgammaCH[i][j].conjugate()+llgammaH0E[i][j]*llgammaH0E[i][j].conjugate()+llgammaH0M[i][j]*llgammaH0M[i][j].conjugate())*g*g*(1-cos2)*pow((pow(mixes.M[tLepton][iDown][i][i],2)-pow(mixes.M[tLepton][iDown][j][j],2))/mixes.M[tLepton][iDown][i][i],3)/(4*Pi);
459  ex mmuon=mixes.M[tLepton][iDown][i][i];
460  ex A,B;
461 
462  if(j<i){ for(uint k=0;k<3;k++){
463  ex mtau=mixes.M[tLepton][iDown][k][k];
464  B+=-mixes.VN[tLepton][1][k][j].conjugate()*mixes.VN[tLepton][1][k][i]/(12*pow(McH,2));
465  B+=mixes.N[tLepton][1][k][j].conjugate()*mixes.N[tLepton][1][k][i]/12*(pow(MR,-2)+pow(MI,-2));
466 
467  A+=mixes.N[tLepton][1][j][k]*mixes.N[tLepton][1][i][k].conjugate()*(pow(MR,-2)+pow(MI,-2))/12;
468  A+=mixes.N[tLepton][1][j][k]*mixes.N[tLepton][1][k][i]/mtau/mmuon*(Fh2(pow(mtau/MR,2))-Fh2(pow(mtau/MI,2)))/4;
469  }
470  llgamma[i][j]=(A*A.conjugate()+B*B.conjugate())*alpha*pow(mmuon,5)*GF*GF/(128*pow(Pi,4));
471  }
472  else if(j==i){
473  for(uint k=0;k<3;k++){
474  ex mtau=mixes.M[tLepton][iDown][k][k];
475  B+=-mixes.VN[tLepton][1][k][j].conjugate()*mixes.VN[tLepton][1][k][i]/(12*pow(McH,2));
476  B+=mixes.N[tLepton][1][k][j].conjugate()*mixes.N[tLepton][1][k][i]/12*(pow(MR,-2)+pow(MI,-2));
477  B+=mixes.N[tLepton][1][j][k].conjugate()*mixes.N[tLepton][1][i][k]/12*(pow(MR,-2)+pow(MI,-2));
478  }
479  llgamma[i][j]=-B*GF*sqrt(1/2)/(8*pow(Pi,2))*2*mmuon; //e (GeV)^-1=1/(51e6)(e cm) where e=sqrt(alpha*4*Pi)
480  }
481  }
482  add("mutoegamma",llgamma[1][0],new limitedobs(planck/2.197034e-6*2.4e-12));
483  add("tautoegamma",llgamma[2][0],new limitedobs(planck/290.6e-15*3.3e-8));
484  add("tautomugamma",llgamma[2][1],new limitedobs(planck/290.6e-15*4.4e-8));
485 
486  /*add("d_e",abs(llgamma[0][0].imag_part()),new limitedobs(10.5e-28*51e6));
487  add("d_mu",abs(llgamma[1][1].imag_part()),new limitedobs(1.9e-19*51e6));
488  add("d_tau",llgamma[2][2].imag_part(),new gaussobs(-0.85e-17*51e6,0.825/0.85));
489  cout<<"EDM: "<<llgamma[0][0].subs(conjtoabs).subs(replacements).imag_part()<<endl;
490  add("a_mu",-llgamma[1][1].real_part()*2*mixes.M[tLepton][iDown][1][1],new gaussobs(3e-9,1.0/3.0));
491  */
492  /*llgammaCH=Matrix(g*g/(16*Pi*4*Pi*MW*MW*McH*McH*12))*mixes.M[tQuark][iDown]*mixes.VN[1][1].conjugate()*mixes.VN[1][1]; //4+1
493  //Matrix llgammaH0M,llgammaH0E;
494  for(uint i=0;i<3;i++)
495  for(uint j=0;j<3;j++)
496  for(uint k=0;k<3;k++){
497  ex z=pow(mixes.M[tQuark][iUp][i][i]/MR,2);
498  llgammaH0M[j][k]=llgammaH0M[j][k]+(mixes.VN[1][1][j][i].conjugate()*mixes.VN[1][1][k][i]+mixes.VN[1][1][i][j]*mixes.VN[1][1][i][k].conjugate())/pow(mixes.M[tQuark][iDown][i][i],2)*(2*z+6*z*z*log(z))/6;
499  llgammaH0M[j][k]=llgammaH0M[j][k]+(mixes.VN[1][1][i][j]*mixes.VN[1][1][k][i])/mixes.M[tQuark][iDown][i][i]/mixes.M[tQuark][iDown][j][j]*(3*z+2*z*log(z));
500 
501  llgammaH0E[j][k]=llgammaH0E[j][k]+(mixes.VN[1][1][j][i].conjugate()*mixes.VN[1][1][k][i]-mixes.VN[1][1][i][j]*mixes.VN[1][1][i][k].conjugate())/pow(mixes.M[tQuark][iDown][i][i],2)*(2*z+6*z*z*log(z))/6;
502  llgammaH0E[j][k]=llgammaH0E[j][k]+(mixes.VN[1][1][i][j]*mixes.VN[1][1][k][i])/mixes.M[tQuark][iDown][i][i]/mixes.M[tQuark][iDown][j][j]*(3*z+2*z*log(z));
503  }
504  llgammaH0M=Matrix(g*g/(16*Pi*4*Pi*MW*MW*McH*McH))*mixes.M[tQuark][iDown]*llgammaH0M;
505  llgammaH0E=Matrix(g*g/(16*Pi*4*Pi*MW*MW*McH*McH))*mixes.M[tQuark][iDown]*llgammaH0E;
506 
507  //Matrix llgamma;
508  for(uint i=0;i<3;i++)
509  for(uint j=0;j<3;j++)
510  if(j<i) {llgamma[i][j]=(llgammaCH[i][j]*llgammaCH[i][j].conjugate()+llgammaH0E[i][j]*llgammaH0E[i][j].conjugate()+llgammaH0M[i][j]*llgammaH0M[i][j].conjugate())*g*g*(1-cos2)*pow((pow(mixes.M[tQuark][iDown][i][i],2)-pow(mixes.M[tQuark][iDown][j][j],2))/mixes.M[tQuark][iDown][i][i],3)/(4*Pi);
511  //llgamma[i][j]=llgamma[i][j].subs(lst(abs(wild()*pow(MH0,-2))==abs(wild())*pow(MH0,-2)));
512  }
513  else llgamma[i][j]=0;
514  */
515 
516 
517  push_back(prediction(new calcubtosgamma2(mixes)));
518 
519  //add("btosgamma",llgamma[2][1],new gaussobs(3.55e-4,sqrt(2)*0.25/3.55),1);
520  //cout<<csrc<<llgamma[2][1]<<endl;
521  //cout<<latex;
522 
523 
524  BR_Htotaunu=(CHdecaycoupling(chiggs,tau,neutrino)+3*CHdecaycoupling(chiggs,strange,charm))/factor(CHdecaycoupling(chiggs,Fermion(tLepton,iDown),neutrino)+3*CHdecaycoupling(chiggs,Fermion(tQuark,iDown),charm)+3*CHdecaycoupling(chiggs,Fermion(tQuark,iDown),up));
525  BR_Htotaunu=BR_Htotaunu.subs(replacements);
526 
527  //BR_toptoHq=decaywidth(top,bottom,chiggs);
528  //ex toptoWb=decaywidth(top,bottom,wboson);
529  //BR_toptoHq=BR_toptoHq/(BR_toptoHq+toptoWb);
530  //BR_toptoHq=BR_toptoHq.subs(replacements);
531 
532  //cout<<"toptoWb "<<toptoWb.subs(replacements).evalf()<<endl;
533 
534  //b to c tau- nu/b to c e- nu
535  //ex btocR=decaywidth(bottom,charm,tau,neutrino,sVector)/(decaywidth(bottom,charm,electron,neutrino,sVector)+decaywidth(bottom,charm,muon,neutrino,sVector));
536  //cout<<btocR.subs(replacements)<<endl;
537 
538  ex BtoDtaunu,BtoD2taunu, BtoDtaunuSM, KtoPi;
539  for(uint i=0; i<3; i++){
540  ex Wcoup=wboson.couplingL(charm,bottom)*wboson.couplingdaggerL(tau,Fermion(tLepton,iUp,FFlavour(i)));
541  if(Wcoup.subs(replacements)==ex(0)) continue;
542  ex chcoup_Wcoup=-pow(MW/McH,2)*(chiggs.couplingR(charm,bottom)+chiggs.couplingL(charm,bottom))*chiggs.couplingdaggerL(tau,Fermion(tLepton,iUp,FFlavour(i)))/Wcoup;
543  ex chcoup2_Wcoup=-pow(MW/McH,2)*(chiggs.couplingR(charm,bottom)-chiggs.couplingL(charm,bottom))*chiggs.couplingdaggerL(tau,Fermion(tLepton,iUp,FFlavour(i)))/Wcoup;
544 
545  BtoDtaunuSM+=Wcoup*Wcoup.conjugate();
546  BtoDtaunu+=Wcoup*Wcoup.conjugate()*(1+1.5*chcoup_Wcoup.real_part()+chcoup_Wcoup.conjugate()*chcoup_Wcoup);
547  BtoD2taunu+=Wcoup*Wcoup.conjugate()*(1+0.12*chcoup2_Wcoup.real_part()+0.05*chcoup2_Wcoup.conjugate()*chcoup2_Wcoup);
548  }
549  lst r2(pow(mixes.V[1][1][2].imag_part(),2)==pow(abs(mixes.V[1][1][2]),2)-pow(mixes.V[1][1][2].real_part(),2));
550  r2.append(pow(mixes.V[0][2][2].imag_part(),2)==pow(abs(mixes.V[0][2][2]),2)-pow(mixes.V[0][2][2].real_part(),2));
551 
552  r2.append(mixes.M[1][0][1][1]==0);
553  r2.append(pow(abs(mixes.V[0][2][2]),2)==1-pow(abs(mixes.V[0][1][2]),2)-pow(abs(mixes.V[0][0][2]),2));
554  r2.append(pow(abs(mixes.V[0][2][1]),2)==1-pow(abs(mixes.V[0][1][1]),2)-pow(abs(mixes.V[0][0][1]),2));
555  r2.append(abs(sqrt(ex(2))* GF)==sqrt(ex(2))* GF);
556 
557  BtoDtaunuSM=collect_common_factors(BtoDtaunuSM.subs(conjtoabs).subs(r2));
558  BtoDtaunu=collect_common_factors(BtoDtaunu.subs(conjtoabs).subs(r2));
559 
560  BtoDtaunuR=(BtoDtaunu/BtoDtaunuSM).subs(replacements).real_part();
561 
562  BtoD2taunu=BtoD2taunu.subs(conjtoabs).subs(r2);
563  BtoD2taunuR=(BtoD2taunu/BtoDtaunuSM).subs(replacements).real_part();
564 
565 
566  //cout<<"BtoDtaunu/BtoDtaunuSM "<<expand(BtoDtaunu/BtoDtaunuSM)<<endl;
567  iBDtaunu=size();
568  add("BtoDtaunu_BtoDtaunuSM",BtoDtaunu/BtoDtaunuSM,new gaussobs(440.0/296, 1.4*58.0/440));
569 
570  iBD2taunu=size();
571  //cout<<"BtoD2taunu/BtoD2taunuSM "<<1+collect_common_factors(expand(BtoD2taunu/BtoDtaunuSM-1))<<endl;
572  add("BtoD2taunu_BtoD2taunuSM",BtoD2taunu/BtoDtaunuSM,new gaussobs(332.0/252, 1.4*24.0/332.0));
573 
574 
575 
576  for(uint j=0; j<2; j++){
577  ex KtoPimunu, KtoPimunuSM;
578  for(uint i=0; i<3; i++){
579  ex Wcoup=wboson.couplingL(up,strange)*wboson.couplingdaggerL(Fermion(tLepton,iDown,FFlavour(j)),Fermion(tLepton,iUp,FFlavour(i)));
580  if(Wcoup.subs(replacements)==ex(0)) continue;
581  ex chcoup_Wcoup=-pow(MW/McH,2)*(chiggs.couplingR(up,strange)+chiggs.couplingL(up,strange))\
582  *chiggs.couplingdaggerL(muon,Fermion(tLepton,iUp,FFlavour(i)))/Wcoup*(pow(MKp,2)-pow(Mpip,2))\
583  /(mixes.mass(Fermion(tLepton,iDown,FFlavour(j)))*(mixes.mass(strange)-mixes.mass(up)));
584  chcoup_Wcoup=collect_common_factors(expand(chcoup_Wcoup));
585  KtoPimunuSM+=collect_common_factors(expand(Wcoup*Wcoup.conjugate()));
586  KtoPimunu+=collect_common_factors(expand(Wcoup*Wcoup.conjugate()*pow(1+chcoup_Wcoup,2)));
587  }
588  KtoPimunuSM=collect_common_factors(expand(KtoPimunuSM.subs(conjtoabs).subs(r2)));
589  KtoPimunu=collect_common_factors(expand(KtoPimunu.subs(conjtoabs).subs(r2)));
590  KtoPimunu=expand(KtoPimunu.subs(replacements).real_part().subs(lst(abs(wild()*pow(MR,-2))==abs(wild())*pow(MR,-2))).subs(lst(log(wild()*pow(MR,-2))==log(wild())-2*log(MR))));
591  KtoPimunu=expand(KtoPimunu.evalf());
592  KtoPimunuSM=expand(KtoPimunuSM.subs(replacements).real_part().subs(lst(abs(wild()*pow(MR,-2))==abs(wild())*pow(MR,-2))).subs(lst(log(wild()*pow(MR,-2))==log(wild())-2*log(MR))));
593  KtoPimunuSM=expand(KtoPimunuSM.evalf());
594  KtoPi+=0.5*log(KtoPimunu/KtoPimunuSM);
595  }
596 
597  add("KtoPi",KtoPi/(pow(MKp,2)-pow(Mpip,2)),new gaussobs(0.08, 0.11/0.08));
598 
599 
600  //add("b to c tau- nu/b to c e- nu", decaywidth(bottom,charm,electron,neutrino), new limitedobs(planck/290.6e-15*2.7e-8));
601 
602  double fD=0.207;
603  ex DDbar=ex(std::pow(fD,2))*mesonmixing(MD0,charm,up);
604  DDbar=expand(DDbar.subs(replacements).subs(lst(abs(wild()*pow(MR,-2))==abs(wild())*pow(MR,-2))).subs(lst(log(wild()*pow(MR,-2))==log(wild())-2*log(MR))));
605  DDbar=expand(DDbar.evalf());
606  ex aDDbar=sqrt(DDbar.real_part()*DDbar.real_part()+DDbar.imag_part()*DDbar.imag_part());
607  add("DDbar",aDDbar,new limitedobs(9.47e-15));
608  cout<<DDbar<<endl;
609 //2|M12|<6.6e-15GeV
610 
611  double fK=0.156;
612  ex KKbar=ex(std::pow(fK,2))*mesonmixing(MK0,strange,down);
613  KKbar=expand(KKbar.subs(replacements).subs(lst(abs(wild()*pow(MR,-2))==abs(wild())*pow(MR,-2))).subs(lst(log(wild()*pow(MR,-2))==log(wild())-2*log(MR))));
614  KKbar=expand(KKbar.evalf());
615  ex aKKbar=sqrt(KKbar.real_part()*KKbar.real_part()+KKbar.imag_part()*KKbar.imag_part());
616  add("KKbar",aKKbar,new limitedobs(3.5e-15));
617  ex eK=0.94*imag_part(KKbar)/3.5e-15/sqrt(2);
618  //add("a_eK",abs(eK),new limitedobs(2.2e-3));
619  add("a_eK",abs(eK),new limitedobs(20*0.0114e-3));
620  cout<<abs(KKbar)<<endl;
621 
622  double fB=0.189;
623  ex Vtb=mixes.V[tQuark][2][2]/mixes.V[tQuark][2][2].conjugate();
624  ex Vtd=mixes.V[tQuark][2][0]/mixes.V[tQuark][2][0].conjugate();
625  ex Vts=mixes.V[tQuark][2][1]/mixes.V[tQuark][2][1].conjugate();
626 
627  ex BBbar=1+ex(std::pow(fB,2))*mesonmixing(MB0,bottom,down)/(3.337e-13*Vtb*Vtd.conjugate());
628  add("BBbarimag",imag_part(BBbar),new gauss2obs(-0.199,0.062));
629  add("BBbarreal",real_part(BBbar),new gauss2obs(0.823,0.143));
630  cout<<BBbar<<endl;
631  BBbar=3.337e-13*Vtb*Vtd.conjugate();
632  cout<<"Bbar "<<(abs(imag_part(BBbar))/abs(BBbar)).subs(replacements)<<endl;
633  double fBs=0.225;
634  ex BsBsbar=1+ex(std::pow(fBs,2))*mesonmixing(MBs0,bottom,strange)/(1.186e-11*Vtb*Vts.conjugate());
635  add("BsBsbarimag",imag_part(BsBsbar),new gauss2obs(0,0.1));
636  add("BsBsbarreal",real_part(BsBsbar),new gauss2obs(0.965,0.133));
637  cout<<BsBsbar<<endl;
638  BsBsbar=1.186e-11*Vtb*Vts.conjugate();
639  cout<<"Bbar "<<(abs(imag_part(BsBsbar))/abs(BsBsbar)).subs(replacements)<<endl;
640 
641  ex McH2=McH*McH;
642  ex MR2=MR*MR;
643  ex MI2=MI*MI;
644 
645  ex cu=collect_common_factors(expand(chiggs.couplingL(top,bottom)))/mixes.mass(top)/(g/MW/sqrt(ex(2)))/mixes.V[1][2][2];
646  cout<<"cu "<<cu<<endl;
647  ex Zbb=(cu-0.72)/McH;
648  add("Zbb",Zbb,new limitedobs(0.0024));
649  cout<<"Zbb "<<Zbb<<endl;
650  cout<<"SIZE "<<size()<<endl;
651 
652  push_back(prediction(new calcuOblique()));
653 }
654 
655 ~BGL(){
656  delete cBmumu;
657  delete cBsmumu;
658  }
659 
660 ex Y(ex x) const{
661  return 1.0113*x/8/(1-x)*(4-x+3*x*log(x)/(1-x));
662  }
663 
664 ex GW(ex x) const{
665  return (94-x*(179+x*(-55+12*x)))/(36*pow(1-x,3))+x*(16+x*(-32+9*x))*log(x)/(6*pow(1-x,4));
666  }
667 
668 ex GH1(ex x) const{
669  return x*(x*((39-14*x)*x-6)+6*x*(3*x-8)*log(x)-19)/(36*pow(x-1,4));
670  //return -x/12;
671  }
672 
673 ex GH2(ex x) const{
674  return x*((x-1)*(11*x-21)+(16-6*x)*log(x))/(6*pow(x-1,3));
675  //return x/2;
676  }
677 
678 ex FW(ex x) const{
679  return (94-x*(179+x*(-55+12*x)))/(36*pow(1-x,3))+x*(16+x*(-32+9*x))*log(x)/(6*pow(1-x,4));
680  }
681 
682 ex FH1(ex x) const{
683  return -x/12;
684  }
685 
686 ex FH2(ex x) const{
687  return x/2;
688  }
689 
690 ex Fh1(ex x) const{
691  //return (2*x+3*pow(x,2)-6*pow(x,3)+pow(x,4)+6*pow(x,2)*log(x))/(6*pow(1-x,4));
692  return x/3;
693  }
694 
695 ex Fh2(ex x) const{
696  //return (-3*x+4*pow(x,2)-pow(x,3)-2*x*log(x))/pow(1-x,3);
697  return -2*(3/2+log(x))*x;
698  }
699 
700 ex A0(ex x) const{
701  return x*(2+3*x-6*x*x+ x*x*x+6*x*log(x))/(24*pow(1-x,4));
702  }
703 
704 ex A1(ex x) const{
705  return x*(-3+4*x-x*x-2*log(x))/(4*pow(1-x,3));
706  }
707 
708 ex A2(ex x) const{
709  return x/(6*pow(1-x,3))*((-7+5*x-8*x*x)/6.0+x*log(x)/(1-x)*(-2+3*x));
710  }
711 
712 ex A3(ex x) const{
713  return (-3+8*x-5*x*x+(6*x-4)*log(x))*x/(6*pow(1-x,3));
714  }
715 
716 void add(const char * s, ex pred, observable * ob, bool sb=0){
717  //cout<<s<<endl;
718  //cout<<"prediction symb"<<pred<<endl;
719  //,pow(sin(wild()), 2) == 1-pow(cos(wild()), 2)
720  //ex p=expand(pred.subs(replacements).real_part().subs(lst(abs(wild()*pow(MR,-2))==abs(wild())*pow(MR,-2))).subs(lst(log(wild()*pow(MR,-2))==log(wild())-2*log(MR))));
721 
722  ex p=pred.subs(replacements).real_part();
723  p=collect_common_factors(expand(p.evalf()));
724  FUNCP_CUBA fp;
725 
726  lst l(tanb,McH,MR,MI);
727 
728  for(uint i=0;i<3;i++){
729  l.append(Mu[i]);
730  l.append(Md[i]);
731  }
732  if(sb) push_back(prediction(ob,p));
733  else {
734  compile_ex(lst(p), l, fp);
735  //cout<<"prediction numeric"<<p<<endl;
736  //cout<<"exp "<<ob->expected()<<endl<<endl;
737  push_back(prediction(ob,fp));
738  }
739  }
740 
741 int veto(const parameters & p, int max=0) const{
742  if(!p.isvalid()) return 1;
743  if(max==1){
744  double mr=p[1].value+p[2].value;
745  if(mr<10 || mr>10000) return 1;
746  mr+=p[3].value;
747  if(mr<10 || mr>10000) return 1;
748  return 0;
749  }
750  else{
751  double mr=p[1].value+p[2].value;
752  if(mr<10 || mr>mmmax) return 1;
753  mr+=p[3].value;
754  if(mr<10 || mr>mmmax) return 1;
755  return 0;
756  }
757  }
758 
759 parameters generateparameters(int max=0) const{
760  parameters p;
761  //x=log_10(tanb)
762  p.push_back(freeparameter(-3,3,r,stepsize));
763  //y=log_10(McH)
764  if(max==1) p.push_back(freeparameter(10,10000,r,stepsize));
765  else p.push_back(freeparameter(10,mmmax,r,stepsize));
766  //log_10(massR)
767  p.push_back(freeparameter(-200,200,r,stepsize));
768  //log_10(massI)
769  p.push_back(freeparameter(-50,50,r,stepsize));
770 
771  return p;
772 }
773 
774 parameters getlist(const parameters & p) const{
775  //cout<<aux<<endl;
776  //double c2=(1+sqrt(1-4*sqrt(ex_to<numeric>(mudecay.subs(lst(tanb==exp(p[0].value),McH==p[1].value))).to_double())))/2;
777 
778  double x=pow(10.0,p[0].value);
779  //double y=pow(10.0,p[1].value);
780  //double z=pow(10.0,p[2].value);
781  //double w=pow(10.0,p[3].value);
782 
783  double y=p[1].value;
784  double z=y+p[2].value;
785  double w=z+p[3].value;
786 
787  parameters pp(p);
788  pp[0].value=x;
789  pp[2].value+=pp[1].value;
790  pp[3].value+=pp[2].value;
791  pp.values=vector<double>();
792  for(uint i=0; i<4; i++) pp.values.push_back(pp[i].value);
793  lst &l=pp.p;
794  l=lst(tanb==x,McH==y,MR==z,MI==w);
795 
796  return pp;
797 }
798 
799 double bsgammawidth(double tanb,double McH,double MR,double MI, int option=0){
800  parameters p=generateparameters();
801  p[0].value=pow(10.0,tanb);
802  p[1].value=McH;
803  p[2].value=MR;
804  p[3].value=MI;
805 
806  calcubtosgamma2 cal(mixes);
807 
808  return cal.width(p,option);
809 }
810 /*
811 ex decaywidth2(const Fermion& f1, const Fermion& ff2, const Fermion& ff3, const Fermion& ff4, BSpin s=sAny) const{
812 
813  Fermion f2=ff2,f3=ff3, f4=ff4;
814 
815  ex ret=0;
816 
817  realsymbol q1("q1"), q2("q2"), q3("q3"), q4("q2"), s3("s3");
818  ex s2=pow(mixes.mass(f1),2);
819 
820  for(uint j=fElectron;j<=fTau;j++)
821  if(ff2.flavour==fAny || ff2.flavour==j){
822  f2.flavour=(FFlavour)j;
823  for(uint k=fElectron;k<=fTau;k++)
824  if(ff3.flavour==fAny || ff3.flavour==k){
825  f3.flavour=(FFlavour)k;
826  for(uint l=fElectron;l<=fTau;l++)
827  if(ff4.flavour==fAny || ff4.flavour==l){
828  f4.flavour=(FFlavour)l;
829  ex v1=0, v2=0;
830  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2),mq3=mixes.mass(f3),mq4=mixes.mass(f4);
831  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
832  ex s4=m2q1+m2q2+m2q3+m2q4-s2-s3;
833 
834  scalar_products sp;
835  sp.add(q2, q3, (s4-m2q2-m2q3)/2);
836  sp.add(q4, q3, (s2-m2q4-m2q3)/2);
837  sp.add(q2, q4, (s3-m2q2-m2q4)/2);
838 
839  sp.add(q2, q2, m2q2);
840  sp.add(q3, q3, m2q3);
841  sp.add(q4, q4, m2q4);
842 
843  ex q10=(s2+mq1*mq1-mq2*mq2)/(2*sqrt(s2)), lq1l=sqrt(q10*q10-mq1*mq1);
844  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
845 
846  for(uint i=0;i<bosons.size();i++)if(bosons[i].s==s || s==sAny){
847  if(bosons[i].s==0){
848  ex a=-(bosons[i].couplingR(f2,f1)-bosons[i].couplingL(f2,f1))*bosons[i].couplingdaggerL(f3,f4)*s2/(mq1+mq2)/pow(bosons[i].mass,2);
849  v1=v1+a*dirac_gammaL();
850  v2=v2+a.conjugate()*dirac_gammaR();
851  a=-(bosons[i].couplingR(f2,f1)-bosons[i].couplingL(f2,f1))*bosons[i].couplingdaggerR(f3,f4)*s2/(mq1+mq2)/pow(bosons[i].mass,2);
852  v1=v1+a*dirac_gammaR();
853  v2=v2+a.conjugate()*dirac_gammaL();
854  }
855  else{
856  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
857  ex a=(bosons[i].couplingR(f2,f1)-bosons[i].couplingL(f2,f1))*bosons[i].couplingdaggerL(f3,f4)/pow(bosons[i].mass,2);
858  v1=v1+a*sl*dirac_gammaL();
859  v2=v2+a.conjugate()*sl*dirac_gammaL();
860  a=(bosons[i].couplingR(f2,f1)-bosons[i].couplingL(f2,f1))*bosons[i].couplingdaggerR(f3,f4)/pow(bosons[i].mass,2);
861  v1=v1+a*sl*dirac_gammaR();
862  v2=v2+a.conjugate()*sl*dirac_gammaR();
863  }
864  }
865  ex vq3=dirac_slash(q3,4)-mq3*dirac_ONE(), vq4=dirac_slash(q4,4)+mq4*dirac_ONE();
866  ex dt=dirac_trace(vq3*v1*vq4*v2).simplify_indexed(sp);
867  cout<<"dt: "<<dt<<endl;
868  ex result=expand(dt*4*lq3l/s2/Pi/128);
869 
870  ret+=result;
871  }
872  }
873 
874  return collect_common_factors(ret.subs(conjtoabs));
875  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
876 }
877 */
878 
879 ex decaywidth(const Fermion& ff1, const Fermion& ff2, const Fermion& ff3, const Fermion& ff4, BSpin s=sAny) const{
880  multivector<ex,4> a(0,bosons.size(),2,2,2);
881  vector<ex> mass(bosons.size(),0);
882  vector<int> op(bosons.size(),0);
883  ex ret=0;
884  Fermion f1=ff1, f2=ff2,f3=ff3, f4=ff4;
885 
886 
887  for(uint i=fElectron;i<=fTau;i=i+1)
888  if(ff1.flavour==fAny || ff1.flavour==i){
889  f1.flavour=(FFlavour)i;
890  for(uint j=fElectron;j<=fTau;j++)
891  if(ff2.flavour==fAny || ff2.flavour==j){
892  f2.flavour=(FFlavour)j;
893  for(uint k=fElectron;k<=fTau;k++)
894  if(ff3.flavour==fAny || ff3.flavour==k){
895  f3.flavour=(FFlavour)k;
896  for(uint l=fElectron;l<=fTau;l++)
897  if(ff4.flavour==fAny || ff4.flavour==l){
898  f4.flavour=(FFlavour)l;
899  for(uint i=0;i<bosons.size();i++) if(bosons[i].s==s || s==sAny){
900  op[i]=bosons[i].s;
901  mass[i]=bosons[i].mass;
902  a[i][0][0][0]=bosons[i].couplingdaggerL(f2,f1)*bosons[i].couplingL(f3,f4);
903  a[i][0][0][1]=bosons[i].couplingdaggerL(f2,f1)*bosons[i].couplingR(f3,f4);
904  a[i][0][1][0]=bosons[i].couplingdaggerR(f2,f1)*bosons[i].couplingL(f3,f4);
905  a[i][0][1][1]=bosons[i].couplingdaggerR(f2,f1)*bosons[i].couplingR(f3,f4);
906 
907  a[i][1][0][0]=bosons[i].couplingdaggerL(f3,f1)*bosons[i].couplingL(f2,f4);
908  a[i][1][0][1]=bosons[i].couplingdaggerL(f3,f1)*bosons[i].couplingR(f2,f4);
909  a[i][1][1][0]=bosons[i].couplingdaggerR(f3,f1)*bosons[i].couplingL(f2,f4);
910  a[i][1][1][1]=bosons[i].couplingdaggerR(f3,f1)*bosons[i].couplingR(f2,f4);
911  }
912 
913  ret+=wc.get_integral_symb(a,mass,op,mixes.mass(f1));
914  //ret+=wc.get_integral(a,mass,op,mixes.massnum(f1),mixes.massnum(f2),mixes.massnum(f3),mixes.massnum(f4))/pow(mixes.massnum(f1),5)*pow(mixes.mass(f1),5);
915  }}}}
916  if(ff2.flavour==ff4.flavour) ret=ret/2;
917  return collect_common_factors(ret.subs(conjtoabs));
918  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
919 }
920 
921 ex get_integral_symb(const multivector<ex,3>& a, ex m1) const{
922  realsymbol s2("s2"), s3("s3");
923  realsymbol q1("q1"), q2("q2"), q3("q3"), q4("q4");
924 
925  ex m2q1=m1*m1;
926 
927  ex vq1=dirac_slash(q2,4)+dirac_slash(q3,4)+dirac_slash(q4,4)+m1*dirac_ONE(), vq2=dirac_slash(q2,4);
928  ex vq3=dirac_slash(q3,4), vq4=dirac_slash(q4,4);
929 
930  ex s4=m2q1-s2-s3;
931  scalar_products sp;
932  sp.add(q2, q3, (s4)/2);
933  sp.add(q4, q3, (s2)/2);
934  sp.add(q2, q4, (s3)/2);
935 
936  sp.add(q2, q2, 0);
937  sp.add(q3, q3, 0);
938  sp.add(q4, q4, 0);
939 
940  multivector<ex,2> v(0,2,2);
941  v[0][0]=dirac_gammaL(); v[0][1]=dirac_gammaR();
942  v[1][0]=dirac_gammaR(); v[1][1]=dirac_gammaL();
943 
944  multivector<ex,5> traces(0,2,2,2,2,2);
945  for(uint k=0;k<2;k++)
946  for(uint l=0;l<2;l++)
947  for(uint m=0;m<2;m++)
948  for(uint n=0;n<2;n++){
949  ex vk=v[k][0];
950  ex vm=v[m][0];
951  ex vl=v[l][1];
952  ex vn=v[n][1];
953 
954  traces[k][l][m][n][0]=dirac_trace(vq2*vk*vq1*vl)*dirac_trace(vq3*vm*vq4*vn);
955  traces[k][l][m][n][1]=-dirac_trace(vq2*vk*vq1*vl*vq3*vm*vq4*vn);
956  }
957 
958  for(uint k=0;k<2;k++)
959  for(uint l=0;l<2;l++)
960  for(uint m=0;m<2;m++)
961  for(uint n=0;n<2;n++)
962  for(uint o=0;o<2;o++)
963  {
964  traces[k][l][m][n][o]=(traces[k][l][m][n][o]).simplify_indexed(sp);
965  }
966 
967  ex q10=(s2+m1*m1)/(2*sqrt(s2)), lq1l=(m1*m1-s2)/(2*sqrt(s2));
968  ex q30=sqrt(s2)/2, lq3l=q30;
969  ex q20=(m1*m1-s2)/(2*m1), lq2l=q20;
970 
971  ex total=0;
972  for(uint k=0;k<2;k++)
973  for(uint l=0;l<2;l++)
974  for(uint m=0;m<2;m++)
975  for(uint n=0;n<2;n++)
976  for(uint r=0;r<2;r++)
977  for(uint s=0;s<2;s++){
978  ex coup=a[r][k][m]*a[s][l][n].conjugate();
979  ex integrand=traces[k][l][m][n][(r+s)%2];
980  integrand=expand(integral(s3, 0, m1*m1-s2, integrand).eval_integ()/lq1l/sqrt(s2)*lq2l/m1/m1);
981  //double mm2=0, mm3=0, m4=0;
982  ex result=integral(s2,0,m1*m1,integrand).eval_integ()/pow(Pi,3)/512;
983  ex partial=result*coup;
984  total=total+partial;
985  }
986  return total;
987 }
988 
989 ex decaywidthtest2(const Fermion& ff1) const{
990  multivector<ex,3> a(0,2,2,2);
991  symbol gLL("g_{LL}"),gLR("g_{LR}"),gRL("g_{RL}"),gRR("g_{RR}"),cLL("c_{LL}"),cLR("c_{LR}"),cRL("c_{RL}"),cRR("c_{RR}");
992 
993  a[0][0][0]=gLL;
994  a[0][0][1]=gLR;
995  a[0][1][0]=gRL;
996  a[0][1][1]=gRR;
997 
998  a[1][0][0]=cLL;
999  a[1][0][1]=cLR;
1000  a[1][1][0]=cRL;
1001  a[1][1][1]=cRR;
1002 
1003  ex ret=get_integral_symb(a,mixes.mass(ff1));
1004  //ret+=wc.get_integral(a,mass,op,mixes.massnum(f1),mixes.massnum(f2),mixes.massnum(f3),mixes.massnum(f4))/pow(mixes.massnum(f1),5)*pow(mixes.mass(f1),5);
1005 
1006  return collect_common_factors(ret.subs(conjtoabs));
1007  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
1008 }
1009 
1010 /*
1011 ex decaywidthtest(const Fermion& f1, const Fermion& f2, const Fermion& ff3, const Fermion& ff4, BSpin s=sAny) const{
1012 
1013  Fermion f1=ff1, f2=ff2,f3=ff3, f4=ff4;
1014 
1015  ex ret=0;
1016 
1017  realsymbol q1("q1"), q2("q2"), q3("q3"), q4("q4");
1018  symbol gL("gL"), gR("gR");
1019  ex s2("s2"), s3("s3");
1020 
1021  for(uint i=fElectron;i<=fTau;i=i+1)
1022  if(ff1.flavour==fAny || ff1.flavour==i){
1023  f1.flavour=(FFlavour)i;
1024  for(uint j=fElectron;j<=fTau;j++)
1025  if(ff2.flavour==fAny || ff2.flavour==j){
1026  f2.flavour=(FFlavour)j;
1027  for(uint k=fElectron;k<=fTau;k++)
1028  if(ff3.flavour==fAny || ff3.flavour==k){
1029  f3.flavour=(FFlavour)k;
1030  for(uint l=fElectron;l<=fTau;l++)
1031  if(ff4.flavour==fAny || ff4.flavour==l){
1032  f4.flavour=(FFlavour)l;
1033  ex v1=0, v2=0;
1034  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2),mq3=mixes.mass(f3),mq4=mixes.mass(f4);
1035  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
1036 
1037  ex s4=m2q1+m2q2+m2q3+m2q4-s2-s3;
1038 
1039  scalar_products sp;
1040 
1041  sp.add(q4, q3, (s2-m2q4-m2q3)/2);
1042  sp.add(q3, q3, m2q3);
1043  sp.add(q4, q4, m2q4);
1044  ex q10=(s2+mq1*mq1-mq2*mq2)/(2*sqrt(s2)), lq1l=sqrt(q10*q10-mq1*mq1);
1045  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
1046 
1047  ex vq1=dirac_slash(q2,4)+dirac_slash(q3,4)+dirac_slash(q4,4)+mq1*dirac_ONE(), vq2=dirac_slash(q2,4)+mq2*dirac_ONE();
1048  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)-mq4*dirac_ONE();
1049 
1050  ex a;
1051 
1052  a=gL;
1053  v1=v1+a*dirac_gammaL();
1054  v2=v2+a.conjugate()*dirac_gammaR();
1055  a=gR;
1056  v1=v1+a*dirac_gammaR();
1057  v2=v2+a.conjugate()*dirac_gammaL();
1058 
1059  / *}
1060  else{
1061  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
1062  ex a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingL(f3,f4)/pow(bosons[i].mass,2);
1063  v1=v1+a*sl*dirac_gammaL();
1064  v2=v2+a.conjugate()*sl*dirac_gammaL();
1065  a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingR(f3,f4)/pow(bosons[i].mass,2);
1066  v1=v1+a*sl*dirac_gammaR();
1067  v2=v2+a.conjugate()*sl*dirac_gammaR();
1068  }* /
1069 
1070  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)-mq4*dirac_ONE();
1071  ex dt=dirac_trace(vq3*v1*vq4*v2).simplify_indexed(sp);
1072  ex result=expand(dt*4*lq3l/s2/Pi/128);
1073 
1074  ret+=result;
1075  }
1076  }
1077  lst ltest;
1078  ltest.append(conjugate(gL)==pow(abs(gL),2)/gL);
1079  ltest.append(conjugate(gR)==pow(abs(gR),2)/gR);
1080  return pow(meson.decay_factor,2)*collect_common_factors(ret.subs(conjtoabs).subs(ltest));
1081  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
1082 }
1083 */
1084 
1085 ex gRR2(const Fermion& f1, const Fermion& f3) const{
1086 
1087  ex ret1=0,ret2=0;
1088  Fermion f2(tLepton,iUp);
1089  Fermion f4(tLepton,iUp);
1090 
1091  for(uint k=fElectron;k<=fTau;k++){
1092  f2.flavour=(FFlavour)k;
1093  for(uint l=fElectron;l<=fTau;l++){
1094  f4.flavour=(FFlavour)l;
1095  for(uint i=0;i<bosons.size();i++)
1096  if(bosons[i].s==sScalar) {
1097  ex x=bosons[i].couplingdaggerR(f2,f1)*bosons[i].couplingL(f3,f4)/pow(bosons[i].mass,2);
1098  ret1+=x*x.conjugate();
1099  }
1100  else if(bosons[i].s==sVector) {
1101  ex x=bosons[i].couplingdaggerL(f2,f1)*bosons[i].couplingL(f3,f4)/pow(bosons[i].mass,2);
1102  ret2+=x*x.conjugate();
1103  }
1104  }}
1105  //r2.append();
1106  ret2=ret2.subs(conjtoabs);
1107  ret1=ret1.subs(conjtoabs);
1108  for(uint i=0;i<3;i++){
1109  ret1=collect_common_factors(expand(ret1.subs(pow(abs(mixes.V[0][2][i]),2)==1-pow(abs(mixes.V[0][1][i]),2)-pow(abs(mixes.V[0][0][i]),2))));
1110  ret2=collect_common_factors(expand(ret2.subs(pow(abs(mixes.V[0][2][i]),2)==1-pow(abs(mixes.V[0][1][i]),2)-pow(abs(mixes.V[0][0][i]),2))));
1111  }
1112 
1113  cout<<ret2<<endl;
1114  return collect_common_factors(ret1/ret2);
1115 }
1116 
1117 ex tautomu_tautoe() const{
1118 
1119  ex ret1=0,ret2=0, rety1=0, rety2=0;
1120 
1121  Fermion f1(tLepton,iDown,fTau);
1122  Fermion f31(tLepton,iDown,fMuon);
1123  Fermion f32(tLepton,iDown,fElectron);
1124 
1125  Fermion f2(tLepton,iUp);
1126  Fermion f4(tLepton,iUp);
1127 
1128 
1129  for(uint k=fElectron;k<=fTau;k++){
1130  f2.flavour=(FFlavour)k;
1131  for(uint l=fElectron;l<=fTau;l++){
1132  f4.flavour=(FFlavour)l;
1133  ex x1=0, x2=0, y1=0, y2=0;
1134  for(uint i=0;i<bosons.size();i++){
1135  if(bosons[i].s==sScalar) {
1136  x1+=bosons[i].couplingdaggerR(f2,f1)*bosons[i].couplingL(f31,f4)/pow(bosons[i].mass,2);
1137  x2+=bosons[i].couplingdaggerR(f2,f1)*bosons[i].couplingL(f32,f4)/pow(bosons[i].mass,2);
1138  }
1139  else if(bosons[i].s==sVector) {
1140  y1+=bosons[i].couplingdaggerL(f2,f1)*bosons[i].couplingL(f31,f4)/pow(bosons[i].mass,2);
1141  y2+=bosons[i].couplingdaggerL(f2,f1)*bosons[i].couplingL(f32,f4)/pow(bosons[i].mass,2);
1142  }
1143  }
1144  ret1+=(x1*y1.conjugate()).real_part();
1145  ret2+=(x2*y2.conjugate()).real_part();
1146  rety1+=y1*y1.conjugate();
1147  rety2+=y2*y2.conjugate();
1148  }}
1149  ret2=(ret2/rety2*mixes.mass(f32)/mixes.mass(f1)).subs(conjtoabs);
1150  ret1=(ret1/rety1*mixes.mass(f31)/mixes.mass(f1)).subs(conjtoabs);
1151  for(uint i=0;i<3;i++){
1152  ret1=collect_common_factors(expand(ret1.subs(pow(abs(mixes.V[0][2][i]),2)==1-pow(abs(mixes.V[0][1][i]),2)-pow(abs(mixes.V[0][0][i]),2))));
1153  ret2=collect_common_factors(expand(ret2.subs(pow(abs(mixes.V[0][2][i]),2)==1-pow(abs(mixes.V[0][1][i]),2)-pow(abs(mixes.V[0][0][i]),2))));
1154  }
1155 
1156  ex x=pow(mixes.mass(f31)/mixes.mass(f1),2);
1157  ex F1=1-8*x+8*pow(x,3)-pow(x,4)-12*pow(x,2)*log(x);
1158  ex g1=1+9*x-9*pow(x,2)-pow(x,3)+6*x*(1+x)*log(x);
1159  ex N1=1+gRR2(f1,f31)/4;
1160 
1161  x=pow(mixes.mass(f32)/mixes.mass(f1),2);
1162 
1163  ex F2=1-8*x+8*pow(x,3)-pow(x,4)-12*pow(x,2)*log(x);
1164  ex g2=1+9*x-9*pow(x,2)-pow(x,3)+6*x*(1+x)*log(x);
1165  ex N2=1+gRR2(f1,f32)/4;
1166 
1167  return collect_common_factors(N1*(F1+2/N1*ret1*g1)/N2/(F2+2/N2*ret2*g2)*F2/F1);
1168 }
1169 
1170 ex mesondw(const Meson & meson, const Fermion& ff3, const Fermion& ff4, BSpin s=sAny) const{
1171 
1172  const Fermion& f1(meson.q1), f2(meson.q2);
1173  ex mesonmass=meson.mass;
1174 
1175  Fermion f3=ff3, f4=ff4;
1176 
1177  ex ret=0;
1178 
1179  realsymbol q3("q3"), q4("q4");
1180  ex s2=pow(mesonmass,2);
1181 
1182  for(uint k=fElectron;k<=fTau;k++)
1183  if(ff3.flavour==fAny || ff3.flavour==k){
1184  f3.flavour=(FFlavour)k;
1185  for(uint l=fElectron;l<=fTau;l++)
1186  if(ff4.flavour==fAny || ff4.flavour==l){
1187  f4.flavour=(FFlavour)l;
1188  ex v1=0, v2=0;
1189  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2),mq3=mixes.mass(f3),mq4=mixes.mass(f4);
1190  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
1191  scalar_products sp;
1192  sp.add(q4, q3, (s2-m2q4-m2q3)/2);
1193  sp.add(q3, q3, m2q3);
1194  sp.add(q4, q4, m2q4);
1195  ex q10=(s2+mq1*mq1-mq2*mq2)/(2*sqrt(s2)), lq1l=sqrt(q10*q10-mq1*mq1);
1196  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
1197 
1198  for(uint i=0;i<bosons.size();i++) if(bosons[i].s==s || s==sAny){
1199  if(bosons[i].s==0){
1200  ex a=-(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingL(f3,f4)*s2/(mq1+mq2)/pow(bosons[i].mass,2);
1201  v1=v1+a*dirac_gammaL();
1202  v2=v2+a.conjugate()*dirac_gammaR();
1203  a=-(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingR(f3,f4)*s2/(mq1+mq2)/pow(bosons[i].mass,2);
1204  v1=v1+a*dirac_gammaR();
1205  v2=v2+a.conjugate()*dirac_gammaL();
1206  }
1207  else{
1208  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
1209  ex a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingL(f3,f4)/pow(bosons[i].mass,2);
1210  v1=v1+a*sl*dirac_gammaL();
1211  v2=v2+a.conjugate()*sl*dirac_gammaL();
1212  a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingR(f3,f4)/pow(bosons[i].mass,2);
1213  v1=v1+a*sl*dirac_gammaR();
1214  v2=v2+a.conjugate()*sl*dirac_gammaR();
1215  }
1216  }
1217  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)-mq4*dirac_ONE();
1218  ex dt=dirac_trace(vq3*v1*vq4*v2).simplify_indexed(sp);
1219  ex result=expand(dt*4*lq3l/s2/Pi/128);
1220 
1221  ret+=result;
1222  }
1223  }
1224 
1225  return pow(meson.decay_factor,2)*collect_common_factors(ret.subs(conjtoabs));
1226  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
1227 }
1228 
1229 
1230 ex mesondwtest(const Meson & meson, const Fermion& ff3, const Fermion& ff4, BSpin s=sAny) const{
1231 
1232  const Fermion& f1(meson.q1), f2(meson.q2);
1233  ex mesonmass=meson.mass;
1234 
1235  Fermion f3=ff3, f4=ff4;
1236 
1237  ex ret=0;
1238 
1239  realsymbol q3("q3"), q4("q4");
1240  symbol gL("gL"), gR("gR"),gVL("gVL"), gVR("gVR");
1241  symbol gS("gS"), gP("gP"), gA("gA");
1242 
1243  ex s2=pow(mesonmass,2);
1244 
1245  for(uint k=fElectron;k<=fTau;k++)
1246  if(ff3.flavour==fAny || ff3.flavour==k){
1247  f3.flavour=(FFlavour)k;
1248  for(uint l=fElectron;l<=fTau;l++)
1249  if(ff4.flavour==fAny || ff4.flavour==l){
1250  f4.flavour=(FFlavour)l;
1251  ex v1=0, v2=0;
1252  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2),mq3=mixes.mass(f3),mq4=mixes.mass(f4);
1253  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
1254  scalar_products sp;
1255  sp.add(q4, q3, (s2-m2q4-m2q3)/2);
1256  sp.add(q3, q3, m2q3);
1257  sp.add(q4, q4, m2q4);
1258  ex q10=(s2+mq1*mq1-mq2*mq2)/(2*sqrt(s2)), lq1l=sqrt(q10*q10-mq1*mq1);
1259  ex q30=(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
1260 
1261  ex a;
1262  /* a=-gL*s2/(mq1+mq2);
1263  v1=v1+a*dirac_gammaL();
1264  v2=v2+a.conjugate()*dirac_gammaR();
1265  a=-gR*s2/(mq1+mq2);
1266  v1=v1+a*dirac_gammaR();
1267  v2=v2+a.conjugate()*dirac_gammaL();
1268 
1269  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
1270  a=gA;
1271  v1=v1+a*sl*dirac_gamma5();
1272  v2=v2+a.conjugate()*sl*dirac_gamma5();
1273 */
1274  a=-gS*s2/(mq1+mq2);
1275  v1=v1+a*dirac_ONE();
1276  v2=v2+a.conjugate()*dirac_ONE();
1277  a=-gP*s2/(mq1+mq2);
1278  v1=v1+a*dirac_gamma5();
1279  v2=v2-a.conjugate()*dirac_gamma5();
1280  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
1281  // a=gA;
1282  // v1=v1+a*sl*dirac_gamma5();
1283  // v2=v2+a.conjugate()*sl*dirac_gamma5();
1284 
1285  /*}
1286  else{
1287  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
1288  ex a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingL(f3,f4)/pow(bosons[i].mass,2);
1289  v1=v1+a*sl*dirac_gammaL();
1290  v2=v2+a.conjugate()*sl*dirac_gammaL();
1291  a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingR(f3,f4)/pow(bosons[i].mass,2);
1292  v1=v1+a*sl*dirac_gammaR();
1293  v2=v2+a.conjugate()*sl*dirac_gammaR();
1294  }*/
1295 
1296  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)-mq4*dirac_ONE();
1297  ex dt=dirac_trace(vq3*v1*vq4*v2).simplify_indexed(sp);
1298  ex result=expand(dt*4*lq3l/s2/Pi/128);
1299 
1300  ret+=result;
1301  }
1302  }
1303  lst ltest;
1304  ltest.append(conjugate(gL)==pow(abs(gL),2)/gL);
1305  ltest.append(conjugate(gR)==pow(abs(gR),2)/gR);
1306  ltest.append(conjugate(gS)==pow(abs(gS),2)/gS);
1307  ltest.append(conjugate(gP)==pow(abs(gP),2)/gP);
1308  ltest.append(conjugate(gA)==pow(abs(gA),2)/gA);
1309 
1310  return pow(meson.decay_factor,2)*collect_common_factors(expand(ret.subs(conjtoabs).subs(ltest)));
1311  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
1312 }
1313 
1314 ex fermiontomeson(const Fermion& ff4, const Fermion& ff3, const Meson & meson, BSpin s=sAny) const{
1315 
1316  const Fermion& f1(meson.q1), f2(meson.q2);
1317  ex mesonmass=meson.mass;
1318 
1319  Fermion f3=ff3, f4=ff4;
1320 
1321  ex ret=0;
1322 
1323  realsymbol q3("q3"), q4("q4");
1324  ex s2=pow(mesonmass,2);
1325 
1326  for(uint k=fElectron;k<=fTau;k++)
1327  if(ff3.flavour==fAny || ff3.flavour==k){
1328  f3.flavour=(FFlavour)k;
1329  for(uint l=fElectron;l<=fTau;l++)
1330  if(ff4.flavour==fAny || ff4.flavour==l){
1331  f4.flavour=(FFlavour)l;
1332  ex v1=0, v2=0;
1333  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2),mq3=mixes.mass(f3),mq4=mixes.mass(f4);
1334  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
1335  scalar_products sp;
1336  sp.add(q4, q3, -(s2-m2q4-m2q3)/2);
1337  sp.add(q3, q3, m2q3);
1338  sp.add(q4, q4, m2q4);
1339  //ex qm0=(s2-mq3*mq3+mq4*mq4)/(2*mq4), lqml=sqrt(qm0*qm0-s2);
1340  ex q30=(-s2+mq3*mq3+mq4*mq4)/(2*mq4), lq3l=sqrt(q30*q30-mq3*mq3);
1341  //ex q30=-(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
1342 
1343  for(uint i=0;i<bosons.size();i++)if(bosons[i].s==s || s==sAny){
1344  if(bosons[i].s==0){
1345  ex a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingL(f3,f4)*s2/(mq1+mq2)/pow(bosons[i].mass,2);
1346  v1=v1+a*dirac_gammaL();
1347  v2=v2+a.conjugate()*dirac_gammaR();
1348  a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingR(f3,f4)*s2/(mq1+mq2)/pow(bosons[i].mass,2);
1349  v1=v1+a*dirac_gammaR();
1350  v2=v2+a.conjugate()*dirac_gammaL();
1351  }
1352  else{
1353  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
1354  ex a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingL(f3,f4)/pow(bosons[i].mass,2);
1355  v1=v1+a*sl*dirac_gammaL();
1356  v2=v2+a.conjugate()*sl*dirac_gammaL();
1357  a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1))*bosons[i].couplingR(f3,f4)/pow(bosons[i].mass,2);
1358  v1=v1+a*sl*dirac_gammaR();
1359  v2=v2+a.conjugate()*sl*dirac_gammaR();
1360  }
1361  }
1362  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)+mq4*dirac_ONE();
1363  ex dt=dirac_trace(vq3*v1*vq4*v2).simplify_indexed(sp);
1364  ex result=expand(dt*2*lq3l/mq4/mq4/Pi/128);
1365 
1366  ret+=result;
1367  }
1368  }
1369 
1370 
1371 
1372  return pow(meson.decay_factor,2)*collect_common_factors(ret.subs(conjtoabs));
1373  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
1374 }
1375 
1376 ex fermiontomesontest(const Fermion& ff4, const Fermion& ff3, const Meson & meson, BSpin s=sAny) const{
1377 
1378  const Fermion& f1(meson.q1), f2(meson.q2);
1379  ex mesonmass=meson.mass;
1380 
1381  Fermion f3=ff3, f4=ff4;
1382 
1383  ex ret=0;
1384 
1385  realsymbol q3("q3"), q4("q4");
1386 
1387  symbol sL("sL"), sR("sR"), vL("vL"), vR("vR");
1388  ex s2=pow(mesonmass,2);
1389 
1390  for(uint k=fElectron;k<=fTau;k++)
1391  if(ff3.flavour==fAny || ff3.flavour==k){
1392  f3.flavour=(FFlavour)k;
1393  for(uint l=fElectron;l<=fTau;l++)
1394  if(ff4.flavour==fAny || ff4.flavour==l){
1395  f4.flavour=(FFlavour)l;
1396  ex v1=0, v2=0;
1397  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2),mq3=mixes.mass(f3),mq4=mixes.mass(f4);
1398  ex m2q1=mq1*mq1, m2q2=mq2*mq2, m2q3=mq3*mq3, m2q4=mq4*mq4;
1399  scalar_products sp;
1400  sp.add(q4, q3, -(s2-m2q4-m2q3)/2);
1401  sp.add(q3, q3, m2q3);
1402  sp.add(q4, q4, m2q4);
1403  //ex qm0=(s2-mq3*mq3+mq4*mq4)/(2*mq4), lqml=sqrt(qm0*qm0-s2);
1404  ex q30=(-s2+mq3*mq3+mq4*mq4)/(2*mq4), lq3l=sqrt(q30*q30-mq3*mq3);
1405  //ex q30=-(s2+mq3*mq3-mq4*mq4)/(2*sqrt(s2)), lq3l=sqrt(q30*q30-mq3*mq3);
1406 
1407 
1408  ex a=sL;
1409  v1=v1+a*dirac_gammaL();
1410  v2=v2+a.conjugate()*dirac_gammaR();
1411  a=sR;
1412  v1=v1+a*dirac_gammaR();
1413  v2=v2+a.conjugate()*dirac_gammaL();
1414 
1415  ex sl=(dirac_slash(q3,4)+dirac_slash(q4,4));
1416  a=vL;
1417  v1=v1+a*sl*dirac_gammaL();
1418  v2=v2+a.conjugate()*sl*dirac_gammaL();
1419  a=vR;
1420  v1=v1+a*sl*dirac_gammaR();
1421  v2=v2+a.conjugate()*sl*dirac_gammaR();
1422 
1423  ex vq3=dirac_slash(q3,4)+mq3*dirac_ONE(), vq4=dirac_slash(q4,4)+mq4*dirac_ONE();
1424  ex dt=dirac_trace(vq3*v1*vq4*v2).simplify_indexed(sp);
1425  ex result=expand(dt*2*lq3l/mq4/mq4/Pi/128);
1426 
1427  ret+=result;
1428  }
1429  }
1430 
1431  return pow(meson.decay_factor,2)*collect_common_factors(ret.subs(conjtoabs));
1432  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
1433 }
1434 
1435 ex mesonmixing(ex mesonmass, const Fermion& f1, const Fermion& f2) const{
1436 
1437  ex ret=0;
1438 
1439  ex v1=0, v2=0;
1440  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2);
1441  ex m2q1=mq1*mq1, m2q2=mq2*mq2;
1442 
1443  for(uint i=0;i<bosons.size();i++)
1444  if(bosons[i].s==0){
1445  ex a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1));
1446  v1=v1+pow(a/bosons[i].mass,2);
1447 
1448  ex b=(bosons[i].couplingdaggerR(f2,f1)+bosons[i].couplingdaggerL(f2,f1));
1449  v2=v2+pow(b/bosons[i].mass,2);
1450  }
1451 
1452  ex fc=mesonmass/(mixes.massnum(f1)+mixes.massnum(f2));
1453  fc=pow(fc,2);
1454 
1455  ret=2*(-v1*(1+11*fc)+v2*(1+fc))*mesonmass/96;
1456 
1457  return collect_common_factors(ret.subs(conjtoabs));
1458  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
1459 }
1460 
1461 ex CHdecaycoupling(Boson higgs, const Fermion& ff3, const Fermion& ff4) const{
1462 
1463  Fermion f3=ff3, f4=ff4;
1464  ex ret=0;
1465  for(uint k=fElectron;k<=fTau;k++)
1466  if(ff3.flavour==fAny || ff3.flavour==k){
1467  f3.flavour=(FFlavour)k;
1468  for(uint l=fElectron;l<=fTau;l++)
1469  if(ff4.flavour==fAny || ff4.flavour==l){
1470  f4.flavour=(FFlavour)l;
1471  ret+=higgs.couplingdaggerL(f3,f4)*higgs.couplingdaggerL(f3,f4).conjugate()+higgs.couplingdaggerR(f3,f4)*higgs.couplingdaggerR(f3,f4).conjugate();
1472  }}
1473  return collect_common_factors(ret.subs(conjtoabs));
1474 }
1475 
1476 
1477 double BranchingRatio(double * xx,double * p){
1478  return ex_to<numeric>(BR_Htotaunu.subs(tanb==pow(10.0,xx[0])).evalf()).to_double();
1479 }
1480 
1481 
1482 double topBranchingRatio(double * xx,double * p){
1483  return ex_to<numeric>(BR_toptoHq.subs(lst(tanb==pow(10.0,xx[0]),McH==xx[1])).evalf()).to_double();
1484 }
1485 
1486 widthcalc wc;
1487 
1488 const double planck;
1489 const possymbol GF, MZ, MW, Mh;
1490 const constant Mpip, Mpi0, MBp,MB0,MBs0, MKp,MK0,MDp,MD0,MDsp,MDs0;
1491 const constant Fpi, FB,FBs, FK,FD,FDs;
1492 ex cos2, g, alpha;
1493 const possymbol tanb, cp, McH, MR, MI, rho;
1494 possymbol Mu[3],Md[3];
1495 vector< Boson > bosons;
1496 
1497 lst replacements;
1498 ex Btaunu;
1499 ex BR_Htotaunu;
1500 ex BR_toptoHq;
1501 ex BtotaunuR;
1502 ex BtoDtaunuR;
1503 ex BtoD2taunuR;
1504 
1505 const Mixes mixes;
1506 lst conjtoabs;
1507 realsymbol mu;
1508 
1509 int iBtaunu, iBDtaunu, iBD2taunu;
1510 vector<int> BGLtype;
1511 
1512 double mmmax, stepsize;
1513 
1514 calcuBmumu * cBmumu;
1515 calcuBmumu * cBsmumu;
1516 
1517 
1518 };
1519 
1520 
1521 }
1522 #endif
BGLmodels::Boson
Gauge boson.
Definition: BGL.h:21
widthcalc
this class calculates decay widths of one lepton to 3 leptons
Definition: widthcalc.h:21
std::Matrix::conjugate
Matrix conjugate() const
computes the hermitian conjugate of the matrix
Definition: multivector.h:125
BGLmodels::calcubtosgamma2
calculus of the constraints coming from the b->s gamma decay
Definition: Formulas.h:279
BGLmodels::BGL::BranchingRatio
double BranchingRatio(double *xx, double *p)
Definition: BGL.h:1477
BGLmodels::BGL::BtoD2taunuR
ex BtoD2taunuR
Definition: BGL.h:1503
observable
A base class representing an experimental measure.
Definition: model.h:35
BGLmodels::Boson::couplingdaggerL
ex couplingdaggerL(const Fermion &f2, const Fermion &f1) const
Definition: BGL.h:51
BGLmodels::Meson
a meson properties
Definition: Formulas.h:47
BGLmodels::BGL::decaywidthtest2
ex decaywidthtest2(const Fermion &ff1) const
Definition: BGL.h:989
BGLmodels::BGL::mu
realsymbol mu
Definition: BGL.h:1507
BGLmodels::BGL::GW
ex GW(ex x) const
Definition: BGL.h:664
BGLmodels::BGL::get_integral_symb
ex get_integral_symb(const multivector< ex, 3 > &a, ex m1) const
Definition: BGL.h:921
BGLmodels::Fermion::isospin
FIsospin isospin
Definition: Formulas.h:38
BGLmodels::BGL::FW
ex FW(ex x) const
Definition: BGL.h:678
BGLmodels::Fermion::helicity
FHelicity helicity
Definition: Formulas.h:41
BGLmodels::Boson::reset
void reset()
Definition: BGL.h:63
BGLmodels::BGL::GH2
ex GH2(ex x) const
Definition: BGL.h:673
BGLmodels::BGL::Fpi
const constant Fpi
Definition: BGL.h:1491
BGLmodels::BGL::tautomu_tautoe
ex tautomu_tautoe() const
Definition: BGL.h:1117
BGLmodels::BGL::A1
ex A1(ex x) const
Definition: BGL.h:704
BGLmodels::BGL::BGL
BGL(int genL=2, int genQ=2, int lup=0, int qup=0, int mssm=0)
Definition: BGL.h:81
BGLmodels::BGL::BGLtype
vector< int > BGLtype
Definition: BGL.h:1510
BGLmodels::BGL::cBmumu
calcuBmumu * cBmumu
Definition: BGL.h:1514
BGLmodels::BGL::Mpip
const constant Mpip
Definition: BGL.h:1490
BGLmodels::BGL::Fh2
ex Fh2(ex x) const
Definition: BGL.h:695
BGLmodels::BGL::mesondw
ex mesondw(const Meson &meson, const Fermion &ff3, const Fermion &ff4, BSpin s=sAny) const
Definition: BGL.h:1170
BGLmodels::BGL::BR_toptoHq
ex BR_toptoHq
Definition: BGL.h:1500
BGLmodels::M_MZ
constexpr double M_MZ
Definition: Formulas.h:58
BGLmodels::BGL::replacements
lst replacements
Definition: BGL.h:1497
BGLmodels::BGL::GH1
ex GH1(ex x) const
Definition: BGL.h:668
BGLmodels::BGL::getlist
parameters getlist(const parameters &p) const
Definition: BGL.h:774
BGLmodels::BGL::topBranchingRatio
double topBranchingRatio(double *xx, double *p)
Definition: BGL.h:1482
BGLmodels::BGL::iBtaunu
int iBtaunu
Definition: BGL.h:1509
BGLmodels::BGL::fermiontomesontest
ex fermiontomesontest(const Fermion &ff4, const Fermion &ff3, const Meson &meson, BSpin s=sAny) const
Definition: BGL.h:1376
BGLmodels::calcuBmumu
calculus of the constraints coming from the B->mu mu decay
Definition: Formulas.h:672
BGLmodels::BGL::decaywidth
ex decaywidth(const Fermion &ff1, const Fermion &ff2, const Fermion &ff3, const Fermion &ff4, BSpin s=sAny) const
Definition: BGL.h:879
BGLmodels::BGL::BtoDtaunuR
ex BtoDtaunuR
Definition: BGL.h:1502
BGLmodels::Boson::couplingL
ex couplingL(const Fermion &f2, const Fermion &f1) const
Definition: BGL.h:26
BGLmodels::BGL::A0
ex A0(ex x) const
Definition: BGL.h:700
BGLmodels::BGL::A3
ex A3(ex x) const
Definition: BGL.h:712
uint
unsigned int uint
Definition: script.cpp:4
limitedobs
An experimental measure which is an upper limit on a parameter with a given Confidence Level...
Definition: model.h:52
BGLmodels::BGL::cBsmumu
calcuBmumu * cBsmumu
Definition: BGL.h:1515
BGLmodels::BGL::fermiontomeson
ex fermiontomeson(const Fermion &ff4, const Fermion &ff3, const Meson &meson, BSpin s=sAny) const
Definition: BGL.h:1314
parameters
vector of parameters
Definition: model.h:177
BGLmodels::BGL::CHdecaycoupling
ex CHdecaycoupling(Boson higgs, const Fermion &ff3, const Fermion &ff4) const
Definition: BGL.h:1461
BGLmodels::calcubtosgamma2::width
double width(const parameters &p, int option=0) const
Definition: Formulas.h:517
gaussobs
An experimental measure of a parameter which is a mean value and a standard deviation.
Definition: model.h:78
BGLmodels::Boson::C
multivector< Matrix, 4 > C
Definition: BGL.h:68
BGLmodels::Fermion::flavour
FFlavour flavour
Definition: Formulas.h:39
BGLmodels::BGL::FH1
ex FH1(ex x) const
Definition: BGL.h:682
BGLmodels::BGL::tanb
const possymbol tanb
Definition: BGL.h:1493
BGLmodels::BGL::mesondwtest
ex mesondwtest(const Meson &meson, const Fermion &ff3, const Fermion &ff4, BSpin s=sAny) const
Definition: BGL.h:1230
parameters::p
lst p
Definition: model.h:231
BGLmodels::BGL::MZ
const possymbol MZ
Definition: BGL.h:1489
BGLmodels::BGL::mixes
const Mixes mixes
Definition: BGL.h:1505
Model
Abstract class for a model.
Definition: model.h:361
BGLmodels::Boson::couplingdaggerR
ex couplingdaggerR(const Fermion &f2, const Fermion &f1) const
Definition: BGL.h:55
BGLmodels::BGL::wc
widthcalc wc
Definition: BGL.h:1486
BGLmodels::BGL::generateparameters
parameters generateparameters(int max=0) const
Definition: BGL.h:759
BGLmodels::Boson::coupling
ex coupling(const Fermion &f2, const Fermion &f1, ex mu)
Definition: BGL.h:59
BGLmodels::Meson::q2
Fermion q2
Definition: Formulas.h:52
BGLmodels::BGL::FH2
ex FH2(ex x) const
Definition: BGL.h:686
BGLmodels::BGL::Y
ex Y(ex x) const
Definition: BGL.h:660
BGLmodels::calcuOblique
calculus of the constraints coming from the oblique parameters
Definition: Formulas.h:221
BGLmodels::Meson::q1
Fermion q1
Definition: Formulas.h:52
BGLmodels
Definition: BGL.h:16
BGLmodels::Mixes
definition of the couplings for the different BGL models
Definition: Formulas.h:105
BGLmodels::BGL::A2
ex A2(ex x) const
Definition: BGL.h:708
BGLmodels::BGL::veto
int veto(const parameters &p, int max=0) const
Definition: BGL.h:741
BGLmodels::BGL::bosons
vector< Boson > bosons
Definition: BGL.h:1495
BGLmodels::BGL::Fh1
ex Fh1(ex x) const
Definition: BGL.h:690
BGLmodels::BGL::BR_Htotaunu
ex BR_Htotaunu
Definition: BGL.h:1499
BGLmodels::BGL::BtotaunuR
ex BtotaunuR
Definition: BGL.h:1501
freeparameter
A parameter which will be fitted in the simulation.
Definition: model.h:124
BGLmodels::BGL::gRR2
ex gRR2(const Fermion &f1, const Fermion &f3) const
Definition: BGL.h:1085
gauss2obs
the same as gaussobs but with a different initializer, such that the uncertainty sigma is absolute ...
Definition: model.h:100
BGLmodels::BGL::mesonmixing
ex mesonmixing(ex mesonmass, const Fermion &f1, const Fermion &f2) const
Definition: BGL.h:1435
BGLmodels::BGL::add
void add(const char *s, ex pred, observable *ob, bool sb=0)
Definition: BGL.h:716
BGLmodels::M_MW
constexpr double M_MW
Definition: Formulas.h:59
prediction
theoretical expression for an experimental measure
Definition: model.h:344
BGLmodels::BGL::stepsize
double stepsize
Definition: BGL.h:1512
BGLmodels::Fermion
a fermion properties
Definition: Formulas.h:32
parameters::isvalid
bool isvalid() const
checks if all the values are between their minimums and maximums
Definition: model.h:193
BGLmodels::BGL::conjtoabs
lst conjtoabs
Definition: BGL.h:1506
parameters::values
vector< double > values
Definition: model.h:232
BGLmodels::BGL::planck
const double planck
Definition: BGL.h:1488
BGLmodels::Boson::s
BSpin s
Definition: BGL.h:66
BGLmodels::Boson::couplingR
ex couplingR(const Fermion &f2, const Fermion &f1) const
Definition: BGL.h:39
BGLmodels::BGL
Implementation of the BGL model.
Definition: BGL.h:78
BGLmodels::BGL::bsgammawidth
double bsgammawidth(double tanb, double McH, double MR, double MI, int option=0)
Definition: BGL.h:799