flavour
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
mainbk.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 <cln/cln.h>
7 #include <cln/float.h>
8 int main(){
9  Digits=5;
10  cln::cl_inhibit_floating_point_underflow=1;
11 
12 
13 
14 
15  //TH1F * pdf1=new TH1F("pdf1","pdf1",npoints,10,500);
16 
17 
18  //TGraph * chi2=new TGraph(npoints);
19  double lmax=0;
20  double gaussmax=0;
21  double tanbmax=0,mcHmax=0, Btaunu=0,BDtaunu=0,BD2taunu=0;
22  uint gLm=0,gQm=0,lupm=0,qupm=0;
23 
24  for(uint gL=0;gL<3;gL++)
25  for(uint gQ=0;gQ<3;gQ++)
26  for(uint lup=0;lup<2;lup++)
27  for(uint qup=0;qup<2;qup++){
28  BGL* m=new BGL(gL,gQ,lup,qup);
29  char name[5]="0000";
30  name[0]+=gL;
31  name[1]+=gQ;
32  name[2]+=lup;
33  name[3]+=qup;
34 
35 
36 
37  /*
38  TF1 * f1 = new TF1("f1",m,&BGL::BranchingRatio,-3,3,0,"BGL","BranchingRatio");
39 
40  TCanvas * c1=new TCanvas("c1","",800,600);
41  f1->Draw();
42  c1->SaveAs("BR.png");
43  TF2 * f2 = new TF2("f2",m,&BGL::topBranchingRatio,-3,3,80,175,0,"BGL","topBranchingRatio");
44 
45  TCanvas * c3=new TCanvas("c3","",800,600);
46  f2->Draw("colz");
47  c3->SaveAs("topBR.png");
48  */
49  uint npoints=100;
50  //TH2F * pdf=new TH2F("pdf","pdf",npoints,-7/log(10.0),7/log(10.0),npoints,10,500);
51  double init=-3, final=3;
52  double init2=1, final2=4;
53 
54  uint steps=1e4;
55  Proposal prop1(m,m->iBtaunu);
56  prop1.findPeaks();
57  Proposal prop2(m,m->iBDtaunu);
58  prop2.findPeaks();
59  Proposal prop3(m,m->iBD2taunu);
60  prop3.findPeaks();
61  Proposal prop4(m,-1);
62  prop4.findPeaks();
63 
64  tanbmax=pow(10.0,prop4.floatPeak.pr[0].value);
65  mcHmax=pow(10.0,prop4.floatPeak.pr[1].value);
66 
67  /*
68 
69  TH2F * limits1=new TH2F("limits1","Likelihood",npoints,init,final,npoints,init2,final2);
70  TH2F * limits2=new TH2F("limits2","Likelihood",npoints,init,final,npoints,init2,final2);
71  TH2F * limits3=new TH2F("limits3","Likelihood",npoints,init,final,npoints,init2,final2);
72  TH2F * limits4=new TH2F("limits4","Likelihood",npoints,init,final,npoints,init2,final2);
73  TH2F * limits5=new TH2F("limits5","Likelihood",npoints,init,final,npoints,init2,final2);
74  TProfile2D * like1=new TProfile2D("like1","like",npoints,init,final,npoints,init2,final2);
75  TProfile2D * like2=new TProfile2D("like2","like",npoints,init,final,npoints,init2,final2);
76  TProfile2D * like3=new TProfile2D("like3","like",npoints,init,final,npoints,init2,final2);
77  TProfile2D * like4=new TProfile2D("like4","like",npoints,init,final,npoints,init2,final2);
78  THStack hs("hs","test stacked histograms");
79 
80 
81  for(uint i=steps;i;i--){
82  prop1.getNextPoint();
83  prop2.getNextPoint();
84  prop3.getNextPoint();
85  prop4.getNextPoint();
86 
87  double Btaunu=m->loglike(prop1.floatPeak.pr,m->iBtaunu);
88  double BDtaunu=m->loglike(prop2.floatPeak.pr,m->iBDtaunu);
89  double BD2taunu=m->loglike(prop3.floatPeak.pr,m->iBD2taunu);
90  double total=m->loglike(prop4.floatPeak.pr,-1);
91 
92  like1->Fill(prop1.floatPeak.pr[0].value, prop1.floatPeak.pr[1].value,Btaunu);
93  like2->Fill(prop2.floatPeak.pr[0].value, prop2.floatPeak.pr[1].value,BDtaunu);
94  like3->Fill(prop3.floatPeak.pr[0].value, prop3.floatPeak.pr[1].value,BD2taunu);
95  like4->Fill(prop4.floatPeak.pr[0].value, prop4.floatPeak.pr[1].value,total);
96 
97  if(i%(steps/10)==0) cout<<"Steps "<<i<<endl;
98  //pdf1->Fill(prop.floatPeak.pr[1].value);
99  }
100 
101  for(uint i=0;i<npoints;i++)
102  for(uint j=0;j<npoints;j++)
103  {
104  double bcmax=0;
105  int binmax=like1->GetBin(i+1,j+1);
106  //for(uint k=0;k<npoints;k++){
107  //int binmax=like1->GetBin(i+1,j+1,j+1);
108  //if(like1->GetBinEntries(bin)){
109  // double bc=like1->GetBinContent(bin);
110  // if(bc>bcmax || binmax==-1) {bcmax=bc; binmax=bin;}
111  // }
112  //}
113  double Btaunu=0;
114  double BDtaunu=0;
115  double BD2taunu=0;
116  double rest=0;
117  double sum=0;
118 
119  if(like1->GetBinEntries(binmax)) {
120  Btaunu=like1->GetBinContent(binmax);
121  //cout<<i<<" "<<j<<" "<<Btaunu<<" ";
122  sum+=Btaunu;
123  Btaunu=TMath::Prob(-2*Btaunu,1);
124  //cout<<Btaunu<<endl;
125  if(Btaunu<0.05) Btaunu=0;
126  }
127  if(like2->GetBinEntries(binmax)) {
128  BDtaunu=like2->GetBinContent(binmax);
129  sum+=BDtaunu;
130  BDtaunu=TMath::Prob(-2*BDtaunu,1);
131  if(BDtaunu<0.05) BDtaunu=0;
132  }
133  if(like3->GetBinEntries(binmax)) {
134  BD2taunu=like3->GetBinContent(binmax);
135  sum+=BD2taunu;
136  BD2taunu=TMath::Prob(-2*BD2taunu,1);
137  if(BD2taunu<0.05) BD2taunu=0;
138  }
139  if((like3->GetBinEntries(binmax) && like2->GetBinEntries(binmax) && like1->GetBinEntries(binmax) && sum>-20)){
140  sum=TMath::Prob(-2*sum,3);
141  if(sum>lmax){
142  lmax=sum;
143  tanbmax=pow(10.0,(i+0.5)*1.0/npoints*6-3);
144  mcHmax=pow(10.0,(j+0.5)*1.0/npoints*3+1);
145  }
146  }else sum=0;
147 
148  if(like4->GetBinEntries(binmax)){
149  rest=like4->GetBinContent(binmax);
150  rest=TMath::Prob(-2*rest,m->size());
151  }
152  int scale=100;
153  limits1->SetBinContent(i+1,j+1,Btaunu*scale);
154  limits2->SetBinContent(i+1,j+1,BDtaunu*scale);
155  limits3->SetBinContent(i+1,j+1,BD2taunu*scale);
156  limits4->SetBinContent(i+1,j+1,rest*scale);
157  limits5->SetBinContent(i+1,j+1,sum*scale);
158  }
159 
160  limits3->GetXaxis()->SetTitle("tan\\beta*MW/MH");
161  limits3->GetYaxis()->SetTitle("cotan\\beta*MW/MH");
162  limits1->SetStats(0);
163  limits2->SetStats(0);
164  limits3->SetStats(0);
165  limits4->SetStats(0);
166 
167  //limits4->SetMarkerStyle(7);
168  limits1->SetMarkerColor(1);
169  //limits1->SetMarkerStyle(7);
170  limits2->SetMarkerColor(2);
171  //limits2->SetMarkerStyle(7);
172  limits3->SetMarkerColor(3);
173  //limits3->SetMarkerStyle(7);
174  limits4->SetMarkerColor(9);
175  //hs.Add(limits4);
176  hs.Add(limits3);
177  hs.Add(limits2);
178  hs.Add(limits1);
179 
180  TCanvas * c2=new TCanvas("c2","",800,600);
181  limits4->Draw("colz");
182  //hs.Draw("nostack");
183  c2->SaveAs((string("pdf_")+string(name)+string(".png")).c_str());
184 
185  TCanvas * c3=new TCanvas("c3","",800,600);
186  limits1->Draw("colz");
187  //hs.Draw("nostack");
188  c3->SaveAs((string("pdf1_")+string(name)+string(".png")).c_str());
189 
190  TCanvas * c4=new TCanvas("c4","",800,600);
191  limits2->Draw("colz");
192  //hs.Draw("nostack");
193  c4->SaveAs((string("pdf2_")+string(name)+string(".png")).c_str());
194 
195 
196  TCanvas * c5=new TCanvas("c5","",800,600);
197  limits3->Draw("colz");
198  //hs.Draw("nostack");
199  c5->SaveAs((string("pdf3_")+string(name)+string(".png")).c_str());
200 
201  TCanvas * c6=new TCanvas("c6","",800,600);
202  limits5->Draw("colz");
203  //hs.Draw("nostack");
204  c6->SaveAs((string("pdf123_")+string(name)+string(".png")).c_str());
205 
206  TCanvas * c7=new TCanvas("c7","",800,600);
207  hs.Draw("nostack");
208  c7->SaveAs((string("pdfall_")+string(name)+string(".png")).c_str());
209  */
210  //cout<<"Lmax "<<lmax<<" "<<TMath::Prob(-2*log(lmax),m->size())<<" GaussMax "<<gaussmax<<endl;
211 
212  //cout<<"tanbmax "<<tanbmax<<" McHmax "<<mcHmax<<endl;
213  cout<<"gL "<<gL<<" gQ "<<gQ<<endl;
214  cout<<"lup "<<lup<<" qup "<<qup<<endl;
215 
216  cout<<"Btotaunu "<<m->BtotaunuR.subs(lst(m->tanb==tanbmax,m->McH==mcHmax))<<" exp "<<1.64/0.79<<" +/- "<<0.23*1.64/0.79<<endl;
217  cout<<"BtoDtaunu "<<m->BtoDtaunuR.subs(lst(m->tanb==tanbmax,m->McH==mcHmax))<<" exp "<<440.0/296<<" +/- "<<1.4*58.0/296<<endl;
218  cout<<"BtoD2taunu "<<m->BtoD2taunuR.subs(lst(m->tanb==tanbmax,m->McH==mcHmax))<<" exp "<<332.0/252<<" +/- "<<1.4*24.0/252<<endl;
219 
220 
221  delete m;
222 }
223 
224  return 0;
225 
226 }
227 
Proposal::findPeaks
void findPeaks(uint ns=1, int max=0)
Definition: MCMC.h:167
Proposal::floatPeak
Peak floatPeak
Definition: MCMC.h:223
uint
unsigned int uint
Definition: script.cpp:4
main
int main()
Definition: mainbk.cpp:8
Proposal
A class containing the parameters of a proposal for the next step in the Markov Chain.
Definition: MCMC.h:162
Peak::pr
parameters pr
Definition: MCMC.h:154