flavour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
main.cpp
Go to the documentation of this file.
1 
14 #include "MCMC.h"
15 #include "BGL.h"
16 #include "TF2.h"
17 #include "TProfile3D.h"
18 #include "THStack.h"
19 #include "TColor.h"
20 #include "TROOT.h"
21 #include "TStyle.h"
22 #include "TGraph.h"
23 #include "TLatex.h"
24 #include "TFile.h"
25 #include "TVector.h"
26 
27 #include <cln/cln.h>
28 #include <cln/float.h>
29 
30 using namespace BGLmodels;
31 
57 int main(int argc, char* argv[]){
58  // Check the number of parameters
59  if(argc < 5){
60  std::cerr<<"Usage: "<<argv[0]<<" gL gQ lup qup [mssm]"<<std::endl;
61  return 1;}
62 
63 
64  int gL=atoi(argv[1]);
65  int gQ=atoi(argv[2]);
66  int lup=atoi(argv[3]);
67  int qup=atoi(argv[4]);
68  int mssm=0;
69  if(argc>5) mssm=atoi(argv[5]);
70 
71 
72  Digits=5;
73  cln::cl_inhibit_floating_point_underflow=1;
74 
75  string ll[2][3]={{"#nu_{1}","#nu_{2}","#nu_{3}"},{"e","#mu","#tau"}};
76  string qq[2][3]={{"u","c","t"},{"d","s","b"}};
77  //Int_t MyPalette[100];
78  Double_t r[] = {1, 0.3};
79  Double_t g[] = {1, 0.3};
80  Double_t b[] = {1, 0.3};
81  Double_t stop[] = {0., 1.0};
82  TColor::CreateGradientColorTable(2, stop, r, g,b, 100);
83  //TH1F * pdf1=new TH1F("pdf1","pdf1",npoints,10,500);
84  //TGraph * chi2=new TGraph(npoints);
85 
86  uint npoints=200;
87  double init1=-3, final1=3;
88  double init2=10, final2=1000;
89  double initBmumu=0, finalBmumu=2;
90  double initBsmumu=0, finalBsmumu=5;
91 
92  multivector<BGL *,4> ms(0,3,3,2,2);
93  multivector<parameters,2> plots(parameters(),npoints,npoints);
94  multivector<double,2> likely(-1000,npoints,npoints);
95  multivector<double,2> likelyB(-1000,npoints,npoints);
96  multivector<double,2> likely2(-1000,npoints,npoints);
97  multivector<double,2> likely3(-1000,npoints,npoints);
98  multivector<double,2> likely4(-1000,npoints,npoints);
99  multivector<double,2> likely5(-1000,npoints,npoints);
100  multivector<double,2> likely6(-1000,npoints,npoints);
101 
102  //ofstream mass("mass.out");
103 
104  //for(int gL=2;gL>=0;gL--)
105  //for(int gQ=2;gQ>=0;gQ--)
106  //for(int lup=0;lup<2;lup++)
107  //for(int qup=0;qup<2;qup++){
108 
109 // for(uint qup=0;qup<2;qup++)
110 //for(uint gL=0;gL<3;gL++)
111 //for(uint lup=0;lup<2;lup++)
112 //for(uint gQ=0;gQ<3;gQ++){
113  //if(gL==0 && gQ==2 && lup==0 && qup==1) {t=1; continue;}
114  //if(t==0) continue;
115  double llmax=-1000,gmax=0, McHmax=1000,MRmax=1000,MImax=1000,tbmax=1;
116  BGL* m=new BGL(gL,gQ,lup,qup);
117  m->mmmax=300;
118  ms[gL][gQ][lup][qup]=m;
119 
120  char name[5]="0000";
121  name[0]+=gL;
122  name[1]+=gQ;
123  name[2]+=lup;
124  name[3]+=qup;
125 
126  /*
127  TF1 * f1 = new TF1("f1",m,&BGL::BranchingRatio,-3,3,0,"BGL","BranchingRatio");
128 
129  TCanvas * c1=new TCanvas("c1","",800,600);
130  f1->Draw();
131  c1->SaveAs("BR.png");
132  TF2 * f2 = new TF2("f2",m,&BGL::topBranchingRatio,-3,3,80,175,0,"BGL","topBranchingRatio");
133 
134  TCanvas * c3=new TCanvas("c3","",800,600);
135  f2->Draw("colz");
136  c3->SaveAs("topBR.png");
137  */
138 
139  Proposal prop4(m);
140  //m->stepsize=1e-4;
141  prop4.findPeaks(100,1);
142  llmax=log(prop4.floatPeak.lmax);
143  tbmax=prop4.floatPeak.pr[0].value;
144  McHmax=prop4.floatPeak.pr[1].value;
145  MRmax=prop4.floatPeak.pr[2].value;
146  MImax=prop4.floatPeak.pr[3].value;
147 
148  m->stepsize=1e-4;
149  prop4.findPeaks(100);
150  double llmax_=log(prop4.floatPeak.lmax);
151  if(llmax_>llmax) {llmax=llmax_;
152  tbmax=prop4.floatPeak.pr[0].value;
153  McHmax=prop4.floatPeak.pr[1].value;
154  MRmax=prop4.floatPeak.pr[2].value;
155  MImax=prop4.floatPeak.pr[3].value;
156  }
157 
158  TH2F * limits4=new TH2F("limits4","Likelihood",npoints,init1,final1,npoints,init2,final2);
159  TH2F * Bmumu_Bsmumu=new TH2F("Bmumu_Bsmumu","Likelihood",npoints,initBmumu,finalBmumu,npoints,initBsmumu,finalBsmumu);
160  TH2F * limits_tb_MR=new TH2F("limits_tb_MR","Likelihood",npoints,init1,final1,npoints,init2,final2);
161  TH2F * limits_tb_MI=new TH2F("limits_tb_MI","Likelihood",npoints,init1,final1,npoints,init2,final2);
162  TH2F * limits_MR_MI=new TH2F("limits_MR_MI","Likelihood",npoints,init2,final2,npoints,init2,final2);
163  TH2F * limits_MR_McH=new TH2F("limits_MR_McH","Likelihood",npoints,init2,final2,npoints,init2,final2);
164  TH2F * limits_MI_McH=new TH2F("limits_MI_McH","Likelihood",npoints,init2,final2,npoints,init2,final2);
165 
166 
167  for(uint i=0;i<npoints;i++)
168  for(uint j=0;j<npoints;j++) {
169  limits4->SetBinContent(i+1,j+1,-1000);
170  Bmumu_Bsmumu->SetBinContent(i+1,j+1,-1000);
171  limits_MR_McH->SetBinContent(i+1,j+1,-1000);
172  limits_MI_McH->SetBinContent(i+1,j+1,-1000);
173  limits_MR_MI->SetBinContent(i+1,j+1,-1000);
174  limits_tb_MR->SetBinContent(i+1,j+1,-1000);
175  plots[i][j]=m->generateparameters();
176  }
177  uint steps=40e6;
178  //uint steps=npoints*npoints;
179  double brtaunu=1;
180  for(uint i=steps;i;i--){
181  //prop1.getNextPoint();
182  //prop2.getNextPoint();
183  //prop3.getNextPoint();
184  double total=0;
185  double gp=0;
186  //if(i==steps) cout<<" total "<<m->loglike(prop4.floatPeak.pr)<<endl;
187  //else{
188  if(i>npoints*npoints){
189  if(i==steps/2) {
190  m->mmmax=1000;
191  m->stepsize=1e-2;
192  prop4.findPeaks(100);
193  llmax_=log(prop4.floatPeak.lmax);
194  if(llmax_>llmax) {llmax=llmax_;
195  tbmax=prop4.floatPeak.pr[0].value;
196  McHmax=prop4.floatPeak.pr[1].value;
197  MRmax=prop4.floatPeak.pr[2].value;
198  MImax=prop4.floatPeak.pr[3].value;
199  }
200  }
201  if(i<steps) prop4.getNextPoint();
202  total=log(prop4.floatPeak.lmax);
203  //total=m->loglike(prop4.floatPeak.pr);
204  //gp=m->gaussprob(prop4.floatPeak.pr);
205  }
206  else{
207  uint x=(i-1)/npoints, y=(i-1)%npoints;
208  prop4.floatPeak.pr[0].value=init1+((x+0.5)*(final1-init1))/npoints;
209  prop4.floatPeak.pr[1].value=init2+((y+0.5)*(final2-init2))/npoints;
210  prop4.floatPeak.pr[2].value=0;
211  prop4.floatPeak.pr[3].value=0;
212  total=m->loglike(prop4.floatPeak.pr);
213  }
214  //}
215 
216 
217  //gp=m->gaussprob_noT(prop4.floatPeak.pr);
218  //double Btaunu=m->loglike(prop1.floatPeak.pr,m->iBtaunu);
219  //double BDtaunu=m->loglike(prop2.floatPeak.pr,m->iBDtaunu);
220  //double BD2taunu=m->loglike(prop3.floatPeak.pr,m->iBD2taunu);
221  double brtaunu0=ex_to<numeric>(m->BR_Htotaunu.subs(m->getlist(prop4.floatPeak.pr).p).evalf()).to_double();
222  if(brtaunu0<brtaunu) brtaunu=brtaunu0;
223  if(total>llmax) {llmax=total; gmax=gp;
224  tbmax=prop4.floatPeak.pr[0].value;
225  McHmax=prop4.floatPeak.pr[1].value;
226  MRmax=prop4.floatPeak.pr[2].value;
227  MImax=prop4.floatPeak.pr[3].value;
228  }
229 
230 
231 
232  uint p1=uint((prop4.floatPeak.pr[0].value-init1)/(final1-init1)*npoints);
233  double mr=prop4.floatPeak.pr[1].value;
234  uint p2=uint((mr-init2)/(final2-init2)*npoints);
235  mr+=prop4.floatPeak.pr[2].value;
236  uint p3=uint((prop4.floatPeak.pr[2].value+prop4.floatPeak.pr[1].value-init2)/(final2-init2)*npoints);
237  mr+=prop4.floatPeak.pr[3].value;
238  uint p4=uint((mr-init2)/(final2-init2)*npoints);
239  if(p1<npoints && p2<npoints && p3<npoints && p4<npoints){
240  if(total>likely[p1][p2]){
241  likely[p1][p2]=total;
242  plots[p1][p2].setvalues(prop4.floatPeak.pr);
243  limits4->SetBinContent(p1+1,p2+1,total);
244  }
245  if(total>likely2[p3][p4]){
246  likely2[p3][p4]=total;
247  //plots[gL][gQ][lup][qup][p1][p2].setvalues(prop4.floatPeak.pr);
248  limits_MR_MI->SetBinContent(p3+1,p4+1,total);
249  }
250 
251  if(total>likely3[p3][p2]){
252  likely3[p3][p2]=total;
253  //plots[gL][gQ][lup][qup][p1][p2].setvalues(prop4.floatPeak.pr);
254  limits_MR_McH->SetBinContent(p3+1,p2+1,total);
255  }
256 
257  if(total>likely4[p1][p3]){
258  likely4[p1][p3]=total;
259  //plots[gL][gQ][lup][qup][p1][p2].setvalues(prop4.floatPeak.pr);
260  limits_tb_MR->SetBinContent(p1+1,p3+1,total);
261  }
262  if(total>likely5[p4][p2]){
263  likely5[p4][p2]=total;
264  //plots[gL][gQ][lup][qup][p1][p2].setvalues(prop4.floatPeak.pr);
265  limits_MI_McH->SetBinContent(p4+1,p2+1,total);
266  }
267  parameters ppr=m->getlist(prop4.floatPeak.pr);
268 
269  double Bmumu=m->cBmumu->obsvalue(ppr)/(1e-10*m->planck/1.519e-12);
270  double Bsmumu=m->cBsmumu->obsvalue(ppr)/(1e-9*m->planck/1.516e-12);
271  //cout<<Bmumu<<" "<<Bsmumu<<endl;
272  uint pB=uint((Bmumu-initBmumu)/(finalBmumu-initBmumu)*npoints);
273  uint pBs=uint((Bsmumu-initBsmumu)/(finalBsmumu-initBsmumu)*npoints);
274  if(pB<npoints && pBs<npoints)
275  if(total>likelyB[pB][pBs]){
276  likely[pB][pBs]=total;
277  Bmumu_Bsmumu->SetBinContent(pB+1,pBs+1,total);
278  }
279 
280  }else cout<<"OUT!!!"<<endl;
281 
282  if(i%100000==0){
283  cout<<"Steps "<<i<<" logtb "<<prop4.floatPeak.pr[0].value<<" logMcH "<<prop4.floatPeak.pr[1].value;
284  cout<<" total "<<total<<" gmax "<<gp<<" "<<brtaunu<<endl;
285  }
286 
287  }
288  //}
289 
290  TFile f((string("h")+string(name)+string(".root")).c_str(),"recreate");
291  TVectorD v(5);
292  v[0] = llmax;
293  v[1]=tbmax;
294  v[2]=McHmax;
295  v[3]=MRmax;
296  v[4]=MImax;
297  v.Write("vllmax");
298  limits4->Write();
299  Bmumu_Bsmumu->Write();
300  limits_MR_MI->Write();
301  limits_MR_McH->Write();
302  limits_MI_McH->Write();
303  limits_tb_MR->Write();
304  f.Close();
305 
306  //for(int gL=2;gL>=0;gL--)
307  //for(int gQ=2;gQ>=0;gQ--)
308  //for(uint lup=0;lup<2;lup++)
309  //for(uint qup=0;qup<2;qup++)
310  uint min1=npoints, min2=npoints, min3=npoints;
311  uint min11=npoints, min21=npoints, min31=npoints;
312  uint min12=npoints, min22=npoints, min32=npoints;
313 
314  for(uint i=0;i<npoints;i++)
315  for(uint j=0;j<npoints;j++){
316  int binmax=limits4->GetBin(i+1,j+1);
317  double rest=limits4->GetBinContent(binmax);
318  if(rest>=llmax) rest=1;
319  else rest=TMath::Prob(-2*(rest-llmax),2);
320  if(rest>=0.05 && j<min1){min1=j;}
321  if(rest>=0.05 && j<min11 && j>uint(180*npoints/990)){min11=j;}
322  if(rest>=0.05 && j<min12 && j>uint(400*npoints/990)){min12=j;}
323  limits4->SetBinContent(i+1,j+1,rest);
324 
325  rest=Bmumu_Bsmumu->GetBinContent(binmax);
326  if(rest>=llmax) rest=1;
327  else rest=TMath::Prob(-2*(rest-llmax),2);
328  Bmumu_Bsmumu->SetBinContent(i+1,j+1,rest);
329 
330  rest=limits_MR_MI->GetBinContent(binmax);
331  if(rest>=llmax) rest=1;
332  else rest=TMath::Prob(-2*(rest-llmax),2);
333  limits_MR_MI->SetBinContent(i+1,j+1,rest);
334 
335  rest=limits_MR_McH->GetBinContent(binmax);
336  if(rest>=llmax) rest=1;
337  else rest=TMath::Prob(-2*(rest-llmax),2);
338  if(rest>=0.05 && i<min2){min2=i;}
339  if(rest>=0.05 && i<min21 && j>uint(180*npoints/990)){min21=i;}
340  if(rest>=0.05 && i<min22 && j>uint(400*npoints/990)){min22=i;}
341  limits_MR_McH->SetBinContent(i+1,j+1,rest);
342 
343  rest=limits_MI_McH->GetBinContent(binmax);
344  if(rest>=llmax) rest=1;
345  else rest=TMath::Prob(-2*(rest-llmax),2);
346  if(rest>=0.05 && i<min3){min3=i;}
347  if(rest>=0.05 && i<min31 && j>uint(180*npoints/990)){min31=i;}
348  if(rest>=0.05 && i<min32 && j>uint(400*npoints/990)){min32=i;}
349  limits_MI_McH->SetBinContent(i+1,j+1,rest);
350 
351 
352  rest=limits_tb_MR->GetBinContent(binmax);
353  if(rest>=llmax) rest=1;
354  else rest=TMath::Prob(-2*(rest-llmax),2);
355  limits_tb_MR->SetBinContent(i+1,j+1,rest);
356  }
357  double mmin1=init2+((min1+0.5)*(final2-init2))/npoints;
358  double mmin2=init2+((min2+0.5)*(final2-init2))/npoints;
359  double mmin3=init2+((min3+0.5)*(final2-init2))/npoints;
360 
361  double mmin11=init2+((min11+0.5)*(final2-init2))/npoints;
362  double mmin21=init2+((min21+0.5)*(final2-init2))/npoints;
363  double mmin31=init2+((min31+0.5)*(final2-init2))/npoints;
364 
365  double mmin12=init2+((min12+0.5)*(final2-init2))/npoints;
366  double mmin22=init2+((min22+0.5)*(final2-init2))/npoints;
367  double mmin32=init2+((min32+0.5)*(final2-init2))/npoints;
368 
369  //mass<<name<<" "<<mmin1<<" "<<mmin2<<" "<<mmin3<<endl;
370 
371  ofstream maxs((string("maxs_")+string(name)+string(".out")).c_str());
372  maxs<<mmin1<<" "<<mmin2<<" "<<mmin3<<endl;
373  maxs<<mmin11<<" "<<mmin21<<" "<<mmin31<<endl;
374  maxs<<mmin12<<" "<<mmin22<<" "<<mmin32<<endl;
375  maxs<<llmax<<endl;
376  //for(uint j=0;j<npoints;j++)
377  //for(uint i=0;i<npoints;i++){
378  // int binmax=limits4->GetBin(i+1,j+1);
379  // maxs<<"("<<i<<","<<j<<"):"<<limits4->GetBinContent(binmax)<<endl;
380  // }
381 
382  maxs.close();
383 
384  limits4->SetStats(0);
385  limits4->GetXaxis()->SetTitle("log_{10}(tan\\beta)");
386  limits4->GetYaxis()->SetTitle("M_{H+} (GeV)");
387 
388  Bmumu_Bsmumu->SetStats(0);
389  Bmumu_Bsmumu->GetXaxis()->SetTitle("Br(B\\to\\mu\\mu)/10^{-10}");
390  Bmumu_Bsmumu->GetYaxis()->SetTitle("Br(B_{s}\\to\\mu\\mu)/10^{-9}");
391 
392 
393  limits_MR_MI->SetStats(0);
394  limits_MR_MI->GetYaxis()->SetTitle("M_{I} (GeV)");
395  limits_MR_MI->GetXaxis()->SetTitle("M_{R} (GeV)");
396 
397 
398 
399  limits_MR_McH->SetStats(0);
400  limits_MR_McH->GetYaxis()->SetTitle("M_{H+} (GeV)");
401  limits_MR_McH->GetXaxis()->SetTitle("M_{R} (GeV)");
402 
403  limits_MI_McH->SetStats(0);
404  limits_MI_McH->GetYaxis()->SetTitle("M_{H+} (GeV)");
405  limits_MI_McH->GetXaxis()->SetTitle("M_{I} (GeV)");
406 
407  limits_tb_MR->SetStats(0);
408  limits_tb_MR->GetXaxis()->SetTitle("log_{10}(tan\\beta)");
409  limits_tb_MR->GetYaxis()->SetTitle("M_{R} (GeV)");
410 
411 
412  gStyle->SetOptTitle(0);
413 
414  Double_t contours[3];
415  contours[0] = 0.003;
416  contours[1] = 0.05;
417  contours[2] = 0.32;
418 
419 
420  double ma=0,me=.2, x0=1,y0=120;
421 
422  limits4->SetContour(3, contours);
423  //limits4->GetYaxis()->SetLabelOffset(0.02);
424  limits4->GetYaxis()->SetLabelSize(0.08);
425  limits4->GetYaxis()->SetTitleSize(0.08);
426  limits4->GetYaxis()->SetTitleOffset(1.2);
427  limits4->GetYaxis()->SetLimits(1,999);
428  //limits4->GetXaxis()->SetLabelOffset(0.02);
429  limits4->GetXaxis()->SetLabelSize(0.08);
430  limits4->GetXaxis()->SetTitleSize(0.08);
431  limits4->GetXaxis()->SetTitleOffset(1.2);
432  limits4->GetXaxis()->SetLimits(-2.99,2.99);
433 
434 
435 
436 
437  gStyle->SetPaperSize(10.,10.);
438 
439  TLatex l;
440  l.SetTextSize(0.08);
441  string ss=qq[qup][gQ]+","+ll[lup][gL];
442 
443  TCanvas * c21=new TCanvas("c21","",int(800*(1+ma+me)),int(600*(1+ma+me)));
444  c21->SetMargin(me,ma,me,ma);
445  c21->SetGrid();
446  //limits4->Draw("CONT Z LIST");
447  limits4->Draw("CONT Z LIST");
448  if(!mssm) l.DrawLatex(x0,y0,ss.c_str());
449 
450  c21->SaveAs((string("pdf_")+string(name)+string(".png")).c_str());
451 
452  delete c21;
453 
454  // Bmumu_Bsmumu->SetContour(3, contours);
455  //limits4->GetYaxis()->SetLabelOffset(0.02);
456  Bmumu_Bsmumu->GetYaxis()->SetLabelSize(0.08);
457  Bmumu_Bsmumu->GetYaxis()->SetTitleSize(0.08);
458  Bmumu_Bsmumu->GetYaxis()->SetTitleOffset(1.2);
459  Bmumu_Bsmumu->GetYaxis()->SetLimits(0.01,4.99);
460  //limits4->GetXaxis()->SetLabelOffset(0.02);
461  Bmumu_Bsmumu->GetXaxis()->SetLabelSize(0.08);
462  Bmumu_Bsmumu->GetXaxis()->SetTitleSize(0.08);
463  Bmumu_Bsmumu->GetXaxis()->SetTitleOffset(1.2);
464  Bmumu_Bsmumu->GetXaxis()->SetLimits(0.01,1.99);
465 
466  TCanvas * cB=new TCanvas("cB","",int(800*(1+ma+me)),int(600*(1+ma+me)));
467  cB->SetMargin(me,ma,me,ma);
468  cB->SetGrid();
469  //limits4->Draw("CONT Z LIST");
470  Bmumu_Bsmumu->Draw("COLZ");
471  l.DrawLatex(x0,y0,ss.c_str());
472 
473  cB->SaveAs((string("Bmumu_Bsmumu_")+string(name)+string(".png")).c_str());
474 
475  delete cB;
476 
477 
478  limits_MR_MI->SetContour(3, contours);
479  //limits4->GetYaxis()->SetLabelOffset(0.02);
480  limits_MR_MI->GetYaxis()->SetLabelSize(0.06);
481  limits_MR_MI->GetYaxis()->SetTitleSize(0.06);
482  limits_MR_MI->GetYaxis()->SetTitleOffset(1.1);
483  limits_MR_MI->GetYaxis()->SetLimits(1,999);
484  //limits4->GetXaxis()->SetLabelOffset(0.02);
485  limits_MR_MI->GetXaxis()->SetLabelSize(0.06);
486  limits_MR_MI->GetXaxis()->SetTitleSize(0.06);
487  limits_MR_MI->GetXaxis()->SetTitleOffset(1.1);
488  limits_MR_MI->GetXaxis()->SetLimits(1,999);
489 
490  TCanvas * c3=new TCanvas("c3","",800,600);
491  c3->SetMargin(me,ma,me,ma);
492  c3->SetGrid();
493 
494  limits_MR_MI->Draw("CONT LIST");
495 
496  c3->SaveAs((string("pdf_")+string(name)+string("_MRMI.png")).c_str());
497 
498  limits_MR_McH->SetContour(3, contours);
499  //limits4->GetYaxis()->SetLabelOffset(0.02);
500  limits_MR_McH->GetYaxis()->SetLabelSize(0.06);
501  limits_MR_McH->GetYaxis()->SetTitleSize(0.06);
502  limits_MR_McH->GetYaxis()->SetTitleOffset(1.1);
503  limits_MR_McH->GetYaxis()->SetLimits(1,999);
504  //limits4->GetXaxis()->SetLabelOffset(0.02);
505  limits_MR_McH->GetXaxis()->SetLabelSize(0.06);
506  limits_MR_McH->GetXaxis()->SetTitleSize(0.06);
507  limits_MR_McH->GetXaxis()->SetTitleOffset(1.1);
508  limits_MR_McH->GetXaxis()->SetLimits(1,999);
509 
510  TCanvas * c4=new TCanvas("c4","",800,600);
511  limits_MR_McH->Draw("CONT LIST");
512  c4->SetMargin(me,ma,me,ma);
513  c4->SetGrid();
514  c4->SaveAs((string("pdf_")+string(name)+string("_MRMcH.png")).c_str());
515 
516  limits_MI_McH->SetContour(3, contours);
517  //limits4->GetYaxis()->SetLabelOffset(0.02);
518  limits_MI_McH->GetYaxis()->SetLabelSize(0.06);
519  limits_MI_McH->GetYaxis()->SetTitleSize(0.06);
520  limits_MI_McH->GetYaxis()->SetTitleOffset(1.1);
521  limits_MI_McH->GetYaxis()->SetLimits(1,999);
522  //limits4->GetXaxis()->SetLabelOffset(0.02);
523  limits_MI_McH->GetXaxis()->SetLabelSize(0.06);
524  limits_MI_McH->GetXaxis()->SetTitleSize(0.06);
525  limits_MI_McH->GetXaxis()->SetTitleOffset(1.1);
526  limits_MI_McH->GetXaxis()->SetLimits(1,999);
527 
528  TCanvas * c6=new TCanvas("c6","",800,600);
529  limits_MI_McH->Draw("CONT LIST");
530  c6->SetMargin(me,ma,me,ma);
531  c6->SetGrid();
532  c6->SaveAs((string("pdf_")+string(name)+string("_MIMcH.png")).c_str());
533 
534  TCanvas * c5=new TCanvas("c5","",800,600);
535  limits_tb_MR->Draw("colz");
536 
537  c5->SaveAs((string("pdf_")+string(name)+string("_tbMR.png")).c_str());
538 
539 
540  cout<<"gL "<<gL<<" gQ "<<gQ<<endl;
541  cout<<"lup "<<lup<<" qup "<<qup<<endl;
542  cout<<"llmax "<<llmax<<" gmax "<<gmax<<endl;
543 
544 
545  delete m;
546  delete limits4;
547  delete Bmumu_Bsmumu;
548  delete limits_tb_MR;
549  delete limits_tb_MI;
550  delete limits_MR_MI;
551  delete limits_MI_McH;
552  delete limits_MR_McH;
553 
554  //mass.close();
555  return 0;
556 
557 }
558 
Proposal::findPeaks
void findPeaks(uint ns=1, int max=0)
Definition: MCMC.h:167
BGLmodels::calcuBmumu::obsvalue
double obsvalue(const parameters &p) const
Definition: Formulas.h:716
std::multivector
A vector of vectors of vectors of... (N times) of class T objects.
Definition: multivector.h:8
BGLmodels::BGL::cBmumu
calcuBmumu * cBmumu
Definition: BGL.h:1514
Proposal::floatPeak
Peak floatPeak
Definition: MCMC.h:223
BGLmodels::BGL::mmmax
double mmmax
Definition: BGL.h:1512
BGLmodels::BGL::getlist
parameters getlist(const parameters &p) const
Definition: BGL.h:774
uint
unsigned int uint
Definition: script.cpp:4
BGLmodels::BGL::cBsmumu
calcuBmumu * cBsmumu
Definition: BGL.h:1515
parameters
vector of parameters
Definition: model.h:177
parameters::p
lst p
Definition: model.h:231
Proposal::getNextPoint
void getNextPoint()
Definition: MCMC.h:210
Model::loglike
double loglike(const parameters &p, bool check=1, int max=0) const
Definition: model.h:381
BGLmodels::BGL::generateparameters
parameters generateparameters(int max=0) const
Definition: BGL.h:759
BGLmodels
Definition: BGL.h:16
BGLmodels::BGL::BR_Htotaunu
ex BR_Htotaunu
Definition: BGL.h:1499
Peak::lmax
double lmax
Definition: MCMC.h:155
main
int main(int argc, char *argv[])
the main function takes the arguments gL gQ lup qup which specify a BGL model and runs the simulation...
Definition: main.cpp:57
BGLmodels::BGL::stepsize
double stepsize
Definition: BGL.h:1512
Proposal
A class containing the parameters of a proposal for the next step in the Markov Chain.
Definition: MCMC.h:162
BGLmodels::BGL::planck
const double planck
Definition: BGL.h:1488
Peak::pr
parameters pr
Definition: MCMC.h:154
BGLmodels::BGL
Implementation of the BGL model.
Definition: BGL.h:78