flavour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
draw.cpp
Go to the documentation of this file.
1 #include "MCMC.h"
2 #include "BGL.h"
3 #include "TF2.h"
4 #include "TProfile3D.h"
5 #include "THStack.h"
6 #include "TColor.h"
7 #include "TROOT.h"
8 #include "TStyle.h"
9 #include "TGraph.h"
10 #include "TLatex.h"
11 #include "TFile.h"
12 #include "TH2F.h"
13 #include "TVector.h"
14 #include "TCanvas.h"
15 #include "TMath.h"
16 #include <iostream>
17 #include <fstream>
18 
19 
20 #include <cln/cln.h>
21 #include <cln/float.h>
22 
23 using namespace BGLmodels;
24 
28 class BGL2: public Model{
29 public:
30 
31 BGL2(int genL=2,int genQ=2, int lup=0, int qup=0, int mssm=0):
32  planck(6.58211928e-25),
33  GF("G_F"),
34  MZ("M_Z"),
35  MW("M_W"),
36  Mpip("Mpip",0.1396,"M_{\\pi^+}",domain::real),
37  Mpi0("Mpi0",0.1349766,"M_{\\pi^0}",domain::real),
38  MBp("MBp",5.279,"M_{B^+}",domain::real),
39  MB0("MB0",5.2795,"M_{B^0}",domain::real),
40  MBs0("MBs0",5.3663,"M_{B_s^0}",domain::real),
41  MKp("MKp",0.493677,"MKp",domain::real),
42  MK0("MK0",0.497614,"MK0",domain::real),
43  MDp("MDp",1.86957,"MDp",domain::real),
44  MD0("MD0",1.86480,"MD0",domain::real),
45  MDsp("MDsp",1.96845,"MDsp",domain::real),
46  MDs0("MDs0",0),
47  Fpi("Fpi",0.132,"Fpi",domain::real),
48  FB("FB",0.189,"FB",domain::real),
49  FBs("FBs",0.225,"FBs",domain::real),
50  FK("FK",0.159,"FK",domain::real),
51  FD("FD",0.208,"FD",domain::real),
52  FDs("FDs",0.248,"FDs",domain::real),
53  //alpha(7.297352e-3*4*M_PI),
54  cos2(pow(MW/MZ,2)),
55  g(sqrt(GF*8/sqrt(ex(2)))*MW),
56  //g(sqrt(4*Pi*alpha/(1-cos2))),
57  tanb("tg\\beta"),
58  cp("cp"),
59  McH("M_{H^+}"),
60  MR("M_{R}"),
61  MI("M_{I}"),
62  Tparam("T_param"),
63  Sparam("S_param"),
64  QCD1("QCD_1"),
65  QCD2("QCD_2"),
66  mixes(tanb,cp, genL,genQ, lup, qup, mssm),
67  mu("\\mu"),
68  BGLtype(4,0),
69  mmmax(1000),
70  stepsize(1e-2)
71  {
72  alpha=pow(g,2)*(1-cos2)/(4*Pi);
73  replacements.append(GF==1.166371e-5);
74  replacements.append(MZ==M_MZ);
75  replacements.append(MW==M_MW);
76 
77  mixes.appendtolst(replacements);
78 
79  replacements.append(Pi==M_PI);
80  replacements.append(sqrt(ex(2))==sqrt(2));
81  replacements.append(Pi==M_PI);
82  replacements.append(sqrt(ex(2))==sqrt(2));
83 
84  Boson boson;
85 
86  realsymbol q3("q3");
87  ex vq3=dirac_slash(q3,4);
88  varidx jmu(mu,4,1);
89 
90  for(uint i=0;i<2;i++)
91  for(uint j=0;j<3;j++)
92  for(uint k=0;k<3;k++){
93  conjtoabs.append(conjugate(mixes.V[i][j][k])==pow(abs(mixes.V[i][j][k]),2)/mixes.V[i][j][k]);
94  }
95 
96  //W+ boson
97  boson.mass=MW;
98  boson.s=sVector;
99 
100  for(uint t=tLepton;t<=tQuark;t++) boson.C[t][iUp][iDown][hLeft]=mixes.V[t]*Matrix(g/sqrt(ex(2)));
101  Boson wboson=boson;
102  bosons.push_back(boson);
103  boson.reset();
104 
105  //H+ boson
106  boson.mass=McH;
107  boson.s=sScalar;
108 
109  for(uint t=tLepton;t<=tQuark;t++)
110  for(uint i=iUp;i<=iDown;i++) boson.C[t][iUp][iDown][i]=mixes.VN[t][i]*Matrix(g/MW/sqrt(ex(2)));
111  Boson chiggs=boson;
112  bosons.push_back(boson);
113  boson.reset();
114 
115  for(int b=bosons.size()-1;b>=0;b--){
116  boson.mass=bosons[b].mass;
117  boson.s=bosons[b].s;
118  if(boson.s==sVector)
119  for(uint t=tLepton;t<=tQuark;t++)
120  for(uint i=iUp;i<=iDown;i++)
121  for(uint j=iUp;j<=iDown;j++)
122  for(uint h=hLeft;h<=hRight;h++){
123  boson.C[t][i][j][h]=bosons[b].C[t][j][i][h].conjugate();
124  }
125  else for(uint t=tLepton;t<=tQuark;t++)
126  for(uint i=iUp;i<=iDown;i++)
127  for(uint j=iUp;j<=iDown;j++)
128  for(uint h=hLeft;h<=hRight;h++){
129  boson.C[t][i][j][hLeft]=bosons[b].C[t][j][i][hRight].conjugate();
130  boson.C[t][i][j][hRight]=bosons[b].C[t][j][i][hLeft].conjugate();
131  }
132  bosons.push_back(boson);
133  boson.reset();
134  }
135 
136  //(R+iI)/sqrt(2) boson
137  boson.mass=MR;
138  boson.s=sScalar;
139 
140  for(uint t=tLepton;t<=tQuark;t++){
141  boson.C[t][iDown][iDown][hRight]=mixes.N[t][iDown]*Matrix(g/MW/ex(2));
142  boson.C[t][iUp][iUp][hLeft]=mixes.N[t][iUp].conjugate()*Matrix(g/MW/ex(2));
143  boson.C[t][iDown][iDown][hLeft]=mixes.N[t][iDown].conjugate()*Matrix(g/MW/ex(2));
144  boson.C[t][iUp][iUp][hRight]=mixes.N[t][iUp]*Matrix(g/MW/ex(2));
145  }
146  bosons.push_back(boson);
147  boson.reset();
148 
149  //(R+iI)/sqrt(2) boson
150  boson.mass=MI;
151  boson.s=sScalar;
152 
153  for(uint t=tLepton;t<=tQuark;t++){
154  boson.C[t][iDown][iDown][hRight]=mixes.N[t][iDown]*Matrix(I*g/MW/ex(2));
155  boson.C[t][iUp][iUp][hLeft]=mixes.N[t][iUp].conjugate()*Matrix(I*g/MW/ex(2));
156  boson.C[t][iDown][iDown][hLeft]=mixes.N[t][iDown].conjugate()*Matrix(-I*g/MW/ex(2));
157  boson.C[t][iUp][iUp][hRight]=mixes.N[t][iUp]*Matrix(-I*g/MW/ex(2));
158  }
159  bosons.push_back(boson);
160  boson.reset();
161 
162  Fermion electron(tLepton,iDown,fElectron);
163  Fermion electronR(tLepton,iDown,fElectron,cParticle,hRight);
164 
165  Fermion muon(tLepton,iDown,fMuon);
166  Fermion muonR(tLepton,iDown,fMuon,cParticle,hRight);
167 
168  Fermion tau(tLepton,iDown,fTau);
169  Fermion tauR(tLepton,iDown,fTau,cParticle,hRight);
170  Fermion neutrino(tLepton,iUp);
171  Fermion neutrinotau(tLepton,iUp,fTau);
172  Fermion neutrinomuon(tLepton,iUp,fMuon);
173  Fermion neutrinoe(tLepton,iUp,fElectron);
174 
175  Fermion up(tQuark,iUp,fElectron);
176  Fermion down(tQuark,iDown,fElectron);
177  Fermion bottom(tQuark,iDown,fTau);
178  Fermion strange(tQuark,iDown,fMuon);
179  Fermion charm(tQuark,iUp,fMuon);
180  Fermion top(tQuark,iUp,fTau);
181 
182  Meson Pi0d(down,down,Mpi0,Fpi);
183  Meson Pi0u(down,down,Mpi0,Fpi);
184  Meson Pip(up,down,Mpip,Fpi);
185  Meson Pim(down,up,Mpip,Fpi);
186 
187  Meson K0(down,strange,MK0,FK);
188  Meson Kp(up,strange,MKp,FK);
189 
190  Meson D0(charm,up,MD0,FD);
191  Meson Dp(charm,down,MDp,FD);
192  Meson Dsp(charm,strange,MDsp,FDs);
193 
194  Meson B0(down,bottom,MB0,FB);
195  Meson Bp(up,bottom,MBp,FB);
196  Meson Bs0(strange,bottom,MBs0,FBs);
197 
198  lst sb;
199  //sb.append(mixes.M[tQuark][iUp][0][0]==0);
200  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));
201  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));
202 
203  //cout<<pow(sqrt(2)/8*pow(g/MW,2),2)<<endl;
204  //cout<<pow(1.166,2)<<endl;
205  double fK=0.156;
206  ex KKbar=ex(std::pow(fK,2))*mesonmixing(MK0,strange,down);
207  ex eK=0.94*imag_part(KKbar)/3.5e-15/sqrt(2);
208  cout<<"KKbar "<<KKbar<<endl;
209  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))));
210  KKbar=expand(KKbar.evalf());
211  ex aKKbar=sqrt(KKbar.real_part()*KKbar.real_part()+KKbar.imag_part()*KKbar.imag_part());
212 
213  eK=0.94*imag_part(KKbar)/3.5e-15/sqrt(2);
214  eK=eK.subs(replacements).real_part();
215  eK=collect_common_factors(expand(eK.evalf()));
216  cout<<"eK"<<eK<<endl;
217  //add("a_eK",abs(eK),new limitedobs(2.2e-3));
218  epsilonK=new calcuex(new limitedobs(2*0.011e-3),abs(eK));
219 
220 
221 }
222 
223 ~BGL2(){epsilonK->~calcuex();}
224 
225 parameters generateparameters(int max=0) const{
226  parameters p;
227  //x=log_10(tanb)
228  p.push_back(freeparameter(-3,3,r,stepsize));
229  //y=log_10(McH)
230  if(max==1) p.push_back(freeparameter(10,10000,r,stepsize));
231  else p.push_back(freeparameter(10,mmmax,r,stepsize));
232  //log_10(massR)
233  p.push_back(freeparameter(-200,200,r,stepsize));
234  //log_10(massI)
235  p.push_back(freeparameter(-50,50,r,stepsize));
236 
237  return p;
238 }
239 
240 
241 parameters getlist(const parameters & p) const{
242  //cout<<aux<<endl;
243  //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;
244 
245  double x=pow(10.0,p[0].value);
246  //double y=pow(10.0,p[1].value);
247  //double z=pow(10.0,p[2].value);
248  //double w=pow(10.0,p[3].value);
249 
250  double y=p[1].value;
251  double z=y+p[2].value;
252  double w=z+p[3].value;
253 
254  parameters pp(p);
255  pp[0].value=x;
256  pp[2].value+=pp[1].value;
257  pp[3].value+=pp[2].value;
258  pp.values=vector<double>();
259  for(uint i=0; i<4; i++) pp.values.push_back(pp[i].value);
260  lst &l=pp.p;
261  l=lst(tanb==x,McH==y,MR==z,MI==w);
262  l.append(QCD1==inter1.Eval(y));
263  l.append(QCD2==inter2.Eval(y));
264 
265  for(uint i=0;i<3;i++){
266  l.append(Mu[i]==Mu_[i].Eval(log(y)));
267  l.append(Md[i]==Md_[i].Eval(log(y)));
268  }
269  return pp;
270 }
271 
272 ex mesonmixing(ex mesonmass, const Fermion& f1, const Fermion& f2) const{
273 
274  ex ret=0;
275 
276  ex v1=0, v2=0;
277  ex mq1=mixes.mass(f1),mq2=mixes.mass(f2);
278  ex m2q1=mq1*mq1, m2q2=mq2*mq2;
279 
280  for(uint i=0;i<bosons.size();i++)
281  if(bosons[i].s==0){
282  ex a=(bosons[i].couplingdaggerR(f2,f1)-bosons[i].couplingdaggerL(f2,f1));
283  v1=v1+pow(a/bosons[i].mass,2);
284 
285  ex b=(bosons[i].couplingdaggerR(f2,f1)+bosons[i].couplingdaggerL(f2,f1));
286  v2=v2+pow(b/bosons[i].mass,2);
287  }
288 
289  ex fc=mesonmass/(mixes.massnum(f1)+mixes.massnum(f2));
290  fc=pow(fc,2);
291 
292  ret=2*(-v1*(1+11*fc)+v2*(1+fc))*mesonmass/96;
293 
294  return collect_common_factors(ret.subs(conjtoabs));
295  //return expand(ret.subs(lst(exp(-I*wild())==1/exp(I*wild()),sin(wild())==sqrt(1-pow(cos(wild()),2)))));
296 }
297 
298 double bsgammawidth(double tanb_,double McH_,double MR_,double MI_, int option=0){
299  parameters p=generateparameters();
300  p[0].value=pow(10.0,tanb_);
301  p[1].value=McH_;
302  p[2].value=MR_;
303  p[3].value=MI_;
304  calcubtosgamma2 cal(mixes);
305 
306  return cal.width(p,option);
307 }
308 
309 double epsK(double tanb_,double McH_,double MR_,double MI_, int option=0){
310  parameters p=generateparameters();
311  p[0].value=pow(10.0,tanb_);
312  p[1].value=McH_;
313  p[2].value=MR_;
314  p[3].value=MI_;
315  p.p=lst(tanb==p[0].value,McH==p[1].value,MR==p[2].value,MI==p[3].value);
316 
317 
318  return epsilonK->error(p);
319 }
320 
321 const double planck;
322 const possymbol GF, MZ, MW, Mh;
323 const constant Mpip, Mpi0, MBp,MB0,MBs0, MKp,MK0,MDp,MD0,MDsp,MDs0;
324 const constant Fpi, FB,FBs, FK,FD,FDs;
325 ex cos2, g, alpha;
326 const possymbol tanb, cp, McH, MR, MI, rho;
327 const realsymbol Tparam, Sparam, QCD1, QCD2;
328 possymbol Mu[3],Md[3];
329 vector< Boson > bosons;
330 
331 lst replacements;
332 ex Btaunu;
333 ex BR_Htotaunu;
334 ex BR_toptoHq;
335 ex BtotaunuR;
336 ex BtoDtaunuR;
337 ex BtoD2taunuR;
338 
339 const Mixes mixes;
340 lst conjtoabs;
341 realsymbol mu;
342 
343 int iBtaunu, iBDtaunu, iBD2taunu;
344 vector<int> BGLtype;
345 ROOT::Math::Interpolator inter1, inter2;
346 ROOT::Math::Interpolator Mu_[3],Md_[3];
347 double mmmax, stepsize;
348 
349 calcuex * epsilonK;
350 
351 };
352 
353 
357 int main(int argc, char* argv[]){
358  // Check the number of parameters
359 
360  if(argc<6){
361  std::cerr<<"Usage: "<<argv[0]<<" inputfile gL gQ lup qup"<<std::endl;
362  return 1;}
363  CD cmu=conj(Vud[1][0])*Vud[1][1];
364  CD umu=conj(Vud[0][0])*Vud[0][1];
365 
366 
367  cout<<"RATIO "<<cmu*cmu<<endl;
368  cout<<umu*umu<<endl;
369 
370  int gL=atoi(argv[2]);
371  int gQ=atoi(argv[3]);
372  int lup=atoi(argv[4]);
373  int qup=atoi(argv[5]);
374  char name[5]="0000";
375 
376  name[0]+=gL;
377  name[1]+=gQ;
378  name[2]+=lup;
379  name[3]+=qup;
380  string ll[2][3]={{"#nu_{1}","#nu_{2}","#nu_{3}"},{"e","#mu","#tau"}};
381  string qq[2][3]={{"u","c","t"},{"d","s","b"}};
382  //Int_t MyPalette[100];
383  Double_t r[] = {1, 0.3};
384  Double_t g[] = {1, 0.3};
385  Double_t b[] = {1, 0.3};
386  Double_t stop[] = {0., 1.0};
387  TColor::CreateGradientColorTable(2, stop, r, g,b, 100);
388  //TH1F * pdf1=new TH1F("pdf1","pdf1",npoints,10,500);
389  //TGraph * chi2=new TGraph(npoints);
390 
391  uint npoints=200;
392  double init1=-3, final1=3;
393  double init2=10, final2=1000;
394  double initBmumu=0, finalBmumu=3;
395  double initBsmumu=0, finalBsmumu=6;
396 
397  double llmax=-1000,McHmax=1000,MRmax=1000,MImax=1000,tbmax=1;
398 
399  TFile *f=new TFile(argv[1],"read");
400  if(!f->IsOpen()) cout<<"NOFILE"<<endl;
401  //f->ShowStreamerInfo();
402 
403  TH2F *limits4,*Bmumu_Bsmumu,*limits_tb_MR,*limits_tb_MI;
404  TH2F *limits_MR_MI,*limits_MR_McH,*limits_MI_McH;
405 
406  f->GetObject("limits4;1",limits4);
407  f->GetObject("Bmumu_Bsmumu;1",Bmumu_Bsmumu);
408  f->GetObject("limits_tb_MR;1",limits_tb_MR);
409  f->GetObject("limits_tb_MI;1",limits_tb_MI);
410  f->GetObject("limits_MR_MI;1",limits_MR_MI);
411  f->GetObject("limits_MR_McH;1",limits_MR_McH);
412  f->GetObject("limits_MI_McH;1",limits_MI_McH);
413 
414  TVectorD* vllmax=NULL;
415 
416  f->GetObject("vllmax;1",vllmax);
417  if(!vllmax) cout<<"ERROR"<<endl;
418  llmax=(*vllmax)[0];
419  //tbmax=(*vllmax)[1];
420  //McHmax=(*vllmax)[2];
421  //MRmax=McHmax+(*vllmax)[3];
422  //MImax=MRmax+(*vllmax)[4];
423  cout<<llmax<<" "<<tbmax<<" "<<McHmax<<" "<<MRmax<<" "<<MImax<<" "<<endl;
424 
425  BGL2* m=new BGL2(gL,gQ,lup,qup);
426  double sm_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,5));
427  double charged_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,1));
428  double neutral_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,2));
429  double neutralR_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,3));
430  double neutralI_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,4));
431  double eK_(m->epsK(tbmax,McHmax,MRmax,MImax));
432 
433  double all_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,0));
434 
435  //for(int gL=2;gL>=0;gL--)
436  //for(int gQ=2;gQ>=0;gQ--)
437  //for(uint lup=0;lup<2;lup++)
438  //for(uint qup=0;qup<2;qup++)
439  uint min1=npoints, min2=npoints, min3=npoints;
440  uint min11=npoints, min21=npoints, min31=npoints;
441  uint min12=npoints, min22=npoints, min32=npoints;
442 
443  for(uint i=0;i<npoints;i++)
444  for(uint j=0;j<npoints;j++){
445  int binmax=limits4->GetBin(i+1,j+1);
446  double rest=limits4->GetBinContent(binmax);
447  if(rest>=llmax) rest=1;
448  else rest=TMath::Prob(-2*(rest-llmax),2);
449  if(rest>=0.05 && j<min1){min1=j;}
450  if(rest>=0.05 && j<min11 && j>uint(180*npoints/990)){min11=j;}
451  if(rest>=0.05 && j<min12 && j>uint(400*npoints/990)){min12=j;}
452  limits4->SetBinContent(i+1,j+1,rest);
453 
454  rest=Bmumu_Bsmumu->GetBinContent(binmax);
455  if(rest>=llmax) rest=1;
456  else rest=TMath::Prob(-2*(rest-llmax),2);
457  //int nn=4;
458  //int ii=(i/nn)*nn, jj=(j/nn)*nn;
459  //for(int iii=ii;iii<ii+n;++iii)
460  //for(int iii=ii;iii<ii+n;++iii)
461  Bmumu_Bsmumu->SetBinContent(i+1,j+1,rest);
462 
463  rest=limits_MR_MI->GetBinContent(binmax);
464  if(rest>=llmax) rest=1;
465  else rest=TMath::Prob(-2*(rest-llmax),2);
466  limits_MR_MI->SetBinContent(i+1,j+1,rest);
467 
468  rest=limits_MR_McH->GetBinContent(binmax);
469  if(rest>=llmax) rest=1;
470  else rest=TMath::Prob(-2*(rest-llmax),2);
471  if(rest>=0.05 && i<min2){min2=i;}
472  if(rest>=0.05 && i<min21 && j>uint(180*npoints/990)){min21=i;}
473  if(rest>=0.05 && i<min22 && j>uint(400*npoints/990)){min22=i;}
474  limits_MR_McH->SetBinContent(i+1,j+1,rest);
475 
476  rest=limits_MI_McH->GetBinContent(binmax);
477  if(rest>=llmax) rest=1;
478  else rest=TMath::Prob(-2*(rest-llmax),2);
479  if(rest>=0.05 && i<min3){min3=i;}
480  if(rest>=0.05 && i<min31 && j>uint(180*npoints/990)){min31=i;}
481  if(rest>=0.05 && i<min32 && j>uint(400*npoints/990)){min32=i;}
482  limits_MI_McH->SetBinContent(i+1,j+1,rest);
483 
484 
485  rest=limits_tb_MR->GetBinContent(binmax);
486  if(rest>=llmax) rest=1;
487  else rest=TMath::Prob(-2*(rest-llmax),2);
488  limits_tb_MR->SetBinContent(i+1,j+1,rest);
489  }
490  double mmin1=init2+((min1+0.5)*(final2-init2))/npoints;
491  double mmin2=init2+((min2+0.5)*(final2-init2))/npoints;
492  double mmin3=init2+((min3+0.5)*(final2-init2))/npoints;
493 
494  double mmin11=init2+((min11+0.5)*(final2-init2))/npoints;
495  double mmin21=init2+((min21+0.5)*(final2-init2))/npoints;
496  double mmin31=init2+((min31+0.5)*(final2-init2))/npoints;
497 
498  double mmin12=init2+((min12+0.5)*(final2-init2))/npoints;
499  double mmin22=init2+((min22+0.5)*(final2-init2))/npoints;
500  double mmin32=init2+((min32+0.5)*(final2-init2))/npoints;
501 
502  //mass<<name<<" "<<mmin1<<" "<<mmin2<<" "<<mmin3<<endl;
503 
504  ofstream maxs((string("maxs_")+string(name)+string(".out")).c_str());
505  maxs<<mmin1<<" "<<mmin2<<" "<<mmin3<<endl;
506  maxs<<mmin11<<" "<<mmin21<<" "<<mmin31<<endl;
507  maxs<<mmin12<<" "<<mmin22<<" "<<mmin32<<endl;
508  maxs<<llmax<<" "<<tbmax<<" "<<McHmax<<" "<<MRmax<<" "<<MImax<<" "<<endl;
509  maxs<<eK_<<endl;
510  maxs<<sm_<<" "<<charged_<<" "<<neutral_<<" "<<neutralR_<<" "<<neutralI_<<" "<<all_<<endl;
511 
512  //for(uint j=0;j<npoints;j++)
513  //for(uint i=0;i<npoints;i++){
514  // int binmax=limits4->GetBin(i+1,j+1);
515  // maxs<<"("<<i<<","<<j<<"):"<<limits4->GetBinContent(binmax)<<endl;
516  // }
517 
518  maxs.close();
519 
520  double ma=0,me=.2, x0=1,y0=120;
521  gStyle->SetOptTitle(0);
522  gStyle->SetPaperSize(10.,10.);
523  TCanvas * c21=new TCanvas("c21","",int(800*(1+ma+me)),int(600*(1+ma+me)));
524  c21->SetMargin(me,ma,me,ma);
525  c21->SetGrid();
526 
527  limits4->SetStats(0);
528  limits4->GetXaxis()->SetTitle("log_{10}(tan\\beta)");
529  limits4->GetYaxis()->SetTitle("M_{H+} (GeV)");
530 
531  Bmumu_Bsmumu->SetStats(0);
532  Bmumu_Bsmumu->GetXaxis()->SetTitle("Br(B\\to\\mu\\mu)/10^{-10}");
533  Bmumu_Bsmumu->GetYaxis()->SetTitle("Br(B_{s}\\to\\mu\\mu)/10^{-9}");
534 
535  limits_MR_MI->SetStats(0);
536  limits_MR_MI->GetYaxis()->SetTitle("M_{I} (GeV)");
537  limits_MR_MI->GetXaxis()->SetTitle("M_{R} (GeV)");
538 
539 
540  limits_MR_McH->SetStats(0);
541  limits_MR_McH->GetYaxis()->SetTitle("M_{H+} (GeV)");
542  limits_MR_McH->GetXaxis()->SetTitle("M_{R} (GeV)");
543 
544  limits_MI_McH->SetStats(0);
545  limits_MI_McH->GetYaxis()->SetTitle("M_{H+} (GeV)");
546  limits_MI_McH->GetXaxis()->SetTitle("M_{I} (GeV)");
547 
548  limits_tb_MR->SetStats(0);
549  limits_tb_MR->GetXaxis()->SetTitle("log_{10}(tan\\beta)");
550  limits_tb_MR->GetYaxis()->SetTitle("M_{R} (GeV)");
551 
552 
553 
554  Double_t contours[3];
555  contours[0] = 0.003;
556  contours[1] = 0.05;
557  contours[2] = 0.32;
558 
559 
560 
561 
562  limits4->SetContour(3, contours);
563  //limits4->GetYaxis()->SetLabelOffset(0.02);
564  limits4->GetYaxis()->SetLabelSize(0.08);
565  limits4->GetYaxis()->SetTitleSize(0.08);
566  limits4->GetYaxis()->SetTitleOffset(1.2);
567  limits4->GetYaxis()->SetLimits(1,999);
568 
569 
570 
571  //limits4->GetXaxis()->SetLabelOffset(0.02);
572  limits4->GetXaxis()->SetLabelSize(0.08);
573  limits4->GetXaxis()->SetTitleSize(0.08);
574  limits4->GetXaxis()->SetTitleOffset(1.2);
575  limits4->GetXaxis()->SetLimits(-2.99,2.99);
576 
577 
578 
579  TLatex l;
580  l.SetTextSize(0.08);
581  string ss=qq[qup][gQ]+","+ll[lup][gL];
582 
583 
584  limits4->Draw("CONT Z LIST");
585  //limits4->Draw("CONT LIST");
586  //limits4->Draw("colz");
587 
588  l.DrawLatex(x0,y0,ss.c_str());
589 
590  c21->SaveAs((string("pdf_")+string(name)+string(".png")).c_str());
591 
592  delete c21;
593  //Bmumu_Bsmumu->SetBit(TH1::kCanRebin);
594  Bmumu_Bsmumu->Rebin2D(2,2);
595 
596  //Bmumu_Bsmumu->SetContour(3, contours);
597  //limits4->GetYaxis()->SetLabelOffset(0.02);
598  Bmumu_Bsmumu->GetYaxis()->SetLabelSize(0.08);
599  Bmumu_Bsmumu->GetYaxis()->SetTitleSize(0.08);
600  Bmumu_Bsmumu->GetYaxis()->SetTitleOffset(0.8);
601  Bmumu_Bsmumu->GetYaxis()->SetLimits(0.01,4.99);
602  //Bmumu_Bsmumu->GetYaxis()->SetRangeUser(0.01, 3.49);
603  // Bmumu_Bsmumu->GetYaxis()->SetLimits(0.01,3.49);
604  //limits4->GetXaxis()->SetLabelOffset(0.02);
605  Bmumu_Bsmumu->GetYaxis()->SetNdivisions(5, kTRUE);
606  Bmumu_Bsmumu->GetXaxis()->SetNdivisions(5, kTRUE);
607 
608  Bmumu_Bsmumu->GetXaxis()->SetLabelSize(0.08);
609  Bmumu_Bsmumu->GetXaxis()->SetTitleSize(0.08);
610  Bmumu_Bsmumu->GetXaxis()->SetTitleOffset(1.2);
611  Bmumu_Bsmumu->GetXaxis()->SetLimits(0.01,1.99);
612  //Bmumu_Bsmumu->GetXaxis()->SetRangeUser(0., 2);
613 
614  TCanvas * cB=new TCanvas("cB","",int(800*(1+ma+me)),int(600*(1+ma+me)));
615  cB->SetMargin(.14,ma,me,ma);
616  cB->SetGrid();
617  //limits4->Draw("CONT Z LIST");
618  Bmumu_Bsmumu->Draw("COLZ");
619 
620  l.DrawLatex(1.5,1,ss.c_str());
621 
622  cB->SaveAs((string("Bmumu_Bsmumu_")+string(name)+string(".png")).c_str());
623 
624  delete cB;
625 
626  limits_MR_MI->SetContour(3, contours);
627  //limits4->GetYaxis()->SetLabelOffset(0.02);
628  limits_MR_MI->GetYaxis()->SetLabelSize(0.06);
629  limits_MR_MI->GetYaxis()->SetTitleSize(0.06);
630  limits_MR_MI->GetYaxis()->SetTitleOffset(1.1);
631  limits_MR_MI->GetYaxis()->SetLimits(1,999);
632  //limits4->GetXaxis()->SetLabelOffset(0.02);
633  limits_MR_MI->GetXaxis()->SetLabelSize(0.06);
634  limits_MR_MI->GetXaxis()->SetTitleSize(0.06);
635  limits_MR_MI->GetXaxis()->SetTitleOffset(1.1);
636  limits_MR_MI->GetXaxis()->SetLimits(1,999);
637 
638  TCanvas * c3=new TCanvas("c3","",800,600);
639  c3->SetMargin(me,ma,me,ma);
640  c3->SetGrid();
641 
642  limits_MR_MI->Draw("CONT LIST");
643 
644  c3->SaveAs((string("pdf_")+string(name)+string("_MRMI.png")).c_str());
645 
646  limits_MR_McH->SetContour(3, contours);
647  //limits4->GetYaxis()->SetLabelOffset(0.02);
648  limits_MR_McH->GetYaxis()->SetLabelSize(0.06);
649  limits_MR_McH->GetYaxis()->SetTitleSize(0.06);
650  limits_MR_McH->GetYaxis()->SetTitleOffset(1.1);
651  limits_MR_McH->GetYaxis()->SetLimits(1,999);
652  //limits4->GetXaxis()->SetLabelOffset(0.02);
653  limits_MR_McH->GetXaxis()->SetLabelSize(0.06);
654  limits_MR_McH->GetXaxis()->SetTitleSize(0.06);
655  limits_MR_McH->GetXaxis()->SetTitleOffset(1.1);
656  limits_MR_McH->GetXaxis()->SetLimits(1,999);
657 
658  TCanvas * c4=new TCanvas("c4","",800,600);
659  limits_MR_McH->Draw("CONT LIST");
660  c4->SetMargin(me,ma,me,ma);
661  c4->SetGrid();
662  c4->SaveAs((string("pdf_")+string(name)+string("_MRMcH.png")).c_str());
663 
664  limits_MI_McH->SetContour(3, contours);
665  //limits4->GetYaxis()->SetLabelOffset(0.02);
666  limits_MI_McH->GetYaxis()->SetLabelSize(0.06);
667  limits_MI_McH->GetYaxis()->SetTitleSize(0.06);
668  limits_MI_McH->GetYaxis()->SetTitleOffset(1.1);
669  limits_MI_McH->GetYaxis()->SetLimits(1,999);
670  //limits4->GetXaxis()->SetLabelOffset(0.02);
671  limits_MI_McH->GetXaxis()->SetLabelSize(0.06);
672  limits_MI_McH->GetXaxis()->SetTitleSize(0.06);
673  limits_MI_McH->GetXaxis()->SetTitleOffset(1.1);
674  limits_MI_McH->GetXaxis()->SetLimits(1,999);
675 
676  TCanvas * c6=new TCanvas("c6","",800,600);
677  limits_MI_McH->Draw("CONT LIST");
678  c6->SetMargin(me,ma,me,ma);
679  c6->SetGrid();
680  c6->SaveAs((string("pdf_")+string(name)+string("_MIMcH.png")).c_str());
681 
682  TCanvas * c5=new TCanvas("c5","",800,600);
683  limits_tb_MR->Draw("colz");
684 
685  c5->SaveAs((string("pdf_")+string(name)+string("_tbMR.png")).c_str());
686 
687  delete m;
688  //mass.close();
689  f->Close();
690  delete f;
691  return 0;
692 
693 }
694 
BGL2::mu
realsymbol mu
Definition: draw.cpp:341
BGLmodels::Boson
Gauge boson.
Definition: BGL.h:21
BGL2::iBtaunu
int iBtaunu
Definition: draw.cpp:343
BGLmodels::calcubtosgamma2
calculus of the constraints coming from the b->s gamma decay
Definition: Formulas.h:279
BGLmodels::Meson
a meson properties
Definition: Formulas.h:47
BGL2::conjtoabs
lst conjtoabs
Definition: draw.cpp:340
BGL2::BtotaunuR
ex BtotaunuR
Definition: draw.cpp:335
BGLmodels::Boson::reset
void reset()
Definition: BGL.h:63
BGL2::planck
const double planck
Definition: draw.cpp:321
BGL2::inter2
ROOT::Math::Interpolator inter2
Definition: draw.cpp:345
BGL2::generateparameters
parameters generateparameters(int max=0) const
Definition: draw.cpp:225
BGL2::MZ
const possymbol MZ
Definition: draw.cpp:322
BGLmodels::M_MZ
constexpr double M_MZ
Definition: Formulas.h:58
BGL2::BtoD2taunuR
ex BtoD2taunuR
Definition: draw.cpp:337
BGL2::mixes
const Mixes mixes
Definition: draw.cpp:339
BGL2::replacements
lst replacements
Definition: draw.cpp:331
BGL2::Fpi
const constant Fpi
Definition: draw.cpp:324
BGL2::Btaunu
ex Btaunu
Definition: draw.cpp:332
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
BGL2
A second implementation of the BGL model, for testing purposes.
Definition: draw.cpp:28
parameters
vector of parameters
Definition: model.h:177
BGL2::BR_toptoHq
ex BR_toptoHq
Definition: draw.cpp:334
BGL2::BtoDtaunuR
ex BtoDtaunuR
Definition: draw.cpp:336
BGLmodels::calcubtosgamma2::width
double width(const parameters &p, int option=0) const
Definition: Formulas.h:517
BGL2::mesonmixing
ex mesonmixing(ex mesonmass, const Fermion &f1, const Fermion &f2) const
Definition: draw.cpp:272
BGLmodels::Boson::C
multivector< Matrix, 4 > C
Definition: BGL.h:68
calcuex
class to do the calculus of a constraint based on a GiNaC symbolic expression
Definition: model.h:282
BGL2::Tparam
const realsymbol Tparam
Definition: draw.cpp:327
parameters::p
lst p
Definition: model.h:231
BGL2::BGLtype
vector< int > BGLtype
Definition: draw.cpp:344
Model
Abstract class for a model.
Definition: model.h:361
BGL2::g
ex g
Definition: draw.cpp:325
main
int main(int argc, char *argv[])
the main function takes the arguments inputfile gL gQ lup qup which specify the file containing the s...
Definition: draw.cpp:357
BGL2::BGL2
BGL2(int genL=2, int genQ=2, int lup=0, int qup=0, int mssm=0)
Definition: draw.cpp:31
BGLmodels
Definition: BGL.h:16
BGLmodels::Mixes
definition of the couplings for the different BGL models
Definition: Formulas.h:105
BGL2::~BGL2
~BGL2()
Definition: draw.cpp:223
BGL2::bsgammawidth
double bsgammawidth(double tanb_, double McH_, double MR_, double MI_, int option=0)
Definition: draw.cpp:298
BGL2::epsK
double epsK(double tanb_, double McH_, double MR_, double MI_, int option=0)
Definition: draw.cpp:309
freeparameter
A parameter which will be fitted in the simulation.
Definition: model.h:124
BGL2::getlist
parameters getlist(const parameters &p) const
Definition: draw.cpp:241
BGL2::bosons
vector< Boson > bosons
Definition: draw.cpp:329
BGLmodels::CD
std::complex< double > CD
Definition: Formulas.h:65
BGLmodels::M_MW
constexpr double M_MW
Definition: Formulas.h:59
BGL2::epsilonK
calcuex * epsilonK
Definition: draw.cpp:349
BGL2::Mpip
const constant Mpip
Definition: draw.cpp:323
BGLmodels::Fermion
a fermion properties
Definition: Formulas.h:32
BGLmodels::Vud
const Matrixx Vud(13.04 *M_PI/180, 0.201 *M_PI/180, 2.38 *M_PI/180, 1.2)
BGL2::BR_Htotaunu
ex BR_Htotaunu
Definition: draw.cpp:333
BGL2::tanb
const possymbol tanb
Definition: draw.cpp:326
parameters::values
vector< double > values
Definition: model.h:232
BGLmodels::Boson::s
BSpin s
Definition: BGL.h:66
BGL2::stepsize
double stepsize
Definition: draw.cpp:347