ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/LimitCalculation.C
Revision: 1.1
Committed: Tue Jul 19 08:15:45 2011 UTC (13 years, 9 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Log Message:
Moved different limit calculation functions to separate module to unclutter Plotting_Functions; moved unneeded functions to the dusty module

File Contents

# Content
1 #include <iostream>
2 #include <vector>
3 #include <sys/stat.h>
4
5 #include <TCut.h>
6 #include <TROOT.h>
7 #include <TCanvas.h>
8 #include <TMath.h>
9 #include <TColor.h>
10 #include <TPaveText.h>
11 #include <TRandom.h>
12 #include <TH1.h>
13 #include <TH2.h>
14 #include <TF1.h>
15 #include <TSQLResult.h>
16 #include <TProfile.h>
17
18 //#include "TTbar_stuff.C"
19 using namespace std;
20
21 using namespace PlottingSetup;
22
23
24 void rediscover_the_top(string mcjzb, string datajzb) {
25 cout << "Hi! today we are going to (try to) rediscover the top!" << endl;
26 TCanvas *c3 = new TCanvas("c3","c3");
27 c3->SetLogy(1);
28 vector<float> binning;
29 //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
30 /*
31 binning.push_back(50);
32 binning.push_back(100);
33 binning.push_back(150);
34 binning.push_back(200);
35 binning.push_back(500);
36
37
38 TH1F *dataprediction = allsamples.Draw("dataprediction", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
39 TH1F *puresignal = allsamples.Draw("puresignal", datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
40 // TH1F *puresignal = allsamples.Draw("puresignal", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("TTJets"));
41 TH1F *observed = allsamples.Draw("observed", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
42 /*
43 ofstream myfile;
44 TH1F *ratio = (TH1F*)observed->Clone();
45 ratio->Divide(dataprediction);
46 ratio->GetYaxis()->SetTitle("Ratio obs/pred");
47 ratio->Draw();
48 c3->SaveAs("testratio.png");
49 myfile.open ("ShapeFit_log.txt");
50 establish_upper_limits(observed,dataprediction,puresignal,"LM4",myfile);
51 myfile.close();
52 */
53
54
55 int nbins=100;
56 float low=0;
57 float hi=500;
58 TCanvas *c4 = new TCanvas("c4","c4",900,900);
59 c4->Divide(2,2);
60 c4->cd(1);
61 c4->cd(1)->SetLogy(1);
62 TH1F *datapredictiont = allsamples.Draw("datapredictiont", "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
63 TH1F *datapredictiono = allsamples.Draw("datapredictiono", "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
64 datapredictiont->Add(datapredictiono,-1);
65 cout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl;
66 vector<TF1*> functions = do_cb_fit_to_plot(datapredictiont,10);
67 datapredictiont->SetMarkerColor(kRed);
68 datapredictiont->SetLineColor(kRed);
69 datapredictiont->Draw();
70 functions[1]->Draw("same");
71 TText *title1 = write_title("Top Background Prediction (JZB<0, with osof subtr)");
72 title1->Draw();
73
74 c4->cd(2);
75 c4->cd(2)->SetLogy(1);
76 TH1F *observedt = allsamples.Draw("observedt", datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
77 observedt->Draw();
78 datapredictiont->Draw("histo,same");
79 functions[1]->Draw("same");
80 TText *title2 = write_title("Observed and predicted background");
81 title2->Draw();
82
83 c4->cd(3);
84 c4->cd(3)->SetLogy(1);
85 // TH1F *ratio = (TH1F*)observedt->Clone();
86
87 TH1F *analytical_background_prediction= new TH1F("analytical_background_prediction","",nbins,low,hi);
88 for(int i=0;i<=nbins;i++) {
89 analytical_background_prediction->SetBinContent(i+1,functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5)));
90 analytical_background_prediction->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5))));
91 }
92 analytical_background_prediction->GetYaxis()->SetTitle("JZB [GeV]");
93 analytical_background_prediction->GetYaxis()->CenterTitle();
94 TH1F *analyticaldrawonly = (TH1F*)analytical_background_prediction->Clone();
95 analytical_background_prediction->SetFillColor(TColor::GetColor("#3399FF"));
96 analytical_background_prediction->SetMarkerSize(0);
97 analytical_background_prediction->Draw("e5");
98 analyticaldrawonly->Draw("histo,same");
99 functions[1]->Draw("same");
100 TText *title3 = write_title("Analytical bg pred histo");
101 title3->Draw();
102
103 c4->cd(4);
104 // c4->cd(4)->SetLogy(1);
105 vector<float> ratio_binning;
106 ratio_binning.push_back(0);
107 ratio_binning.push_back(5);
108 ratio_binning.push_back(10);
109 ratio_binning.push_back(20);
110 ratio_binning.push_back(50);
111 // ratio_binning.push_back(60);
112 /*
113 ratio_binning.push_back(51);
114 ratio_binning.push_back(52);
115 ratio_binning.push_back(53);
116 ratio_binning.push_back(54);
117 ratio_binning.push_back(55);
118 ratio_binning.push_back(56);
119 ratio_binning.push_back(57);
120 ratio_binning.push_back(58);
121 ratio_binning.push_back(59);
122 ratio_binning.push_back(60);
123 // ratio_binning.push_back(70);*/
124 // ratio_binning.push_back(80);
125 // ratio_binning.push_back(90);
126 ratio_binning.push_back(80);
127 // ratio_binning.push_back(110);
128 ratio_binning.push_back(500);
129
130 TH1F *observedtb = allsamples.Draw("observedtb", datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
131 TH1F *datapredictiontb = allsamples.Draw("datapredictiontb", "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
132 TH1F *datapredictiontbo = allsamples.Draw("datapredictiontbo", "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
133 datapredictiontb->Add(datapredictiontbo,-1);
134 TH1F *analytical_background_predictionb = allsamples.Draw("analytical_background_predictionb",datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"mll<2",data, luminosity);
135 for(int i=0;i<=ratio_binning.size();i++) {
136 analytical_background_predictionb->SetBinContent(i+1,functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i)));
137 analytical_background_predictionb->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i))));
138 }
139
140 TH1F *ratio = (TH1F*) observedtb->Clone();
141 ratio->Divide(datapredictiontb);
142
143 for (int i=0;i<=ratio_binning.size();i++) {
144 cout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
145 }
146
147 // ratio->Divide(datapredictiontb);
148 // ratio->Divide(analytical_background_predictionb);
149 // TGraphAsymmErrors *JZBratio= histRatio(observedtb,analytical_background_predictionb,data,ratio_binning);
150 // ratio->Divide(analytical_background_prediction);
151 // ratio->Divide(datapredictiont);
152 // ratio->GetYaxis()->SetTitle("obs/pred");
153 // JZBratio->Draw("AP");
154 ratio->GetYaxis()->SetRangeUser(0,10);
155 ratio->Draw();
156 //analytical_background_predictionb->Draw();
157 // JZBratio->SetTitle("");
158 TText *title4 = write_title("Ratio of observed to predicted");
159 title4->Draw();
160
161 // CompleteSave(c4,"test/ttbar_discovery_dataprediction___analytical_function");
162 CompleteSave(c4,"test/ttbar_discovery_dataprediction__analytical__new_binning_one_huge_bin_from_80");
163
164
165
166
167
168 }
169
170 void calculate_upper_limits(string mcjzb, string datajzb) {
171 write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
172 // write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
173 // rediscover_the_top(mcjzb,datajzb);
174 /*
175 TCanvas *c3 = new TCanvas("c3","c3");
176 c3->SetLogy(1);
177 vector<float> binning;
178 //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
179 binning.push_back(50);
180 binning.push_back(100);
181 binning.push_back(150);
182 binning.push_back(200);
183 binning.push_back(500);
184 TH1F *datapredictiona = allsamples.Draw("datapredictiona", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity);
185 TH1F *datapredictionb = allsamples.Draw("datapredictionb", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
186 TH1F *datapredictionc = allsamples.Draw("datapredictionc", datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
187 TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
188 dataprediction->Add(datapredictionb,-1);
189 dataprediction->Add(datapredictionc);
190 TH1F *puresignal = allsamples.Draw("puresignal", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
191 TH1F *signalpred = allsamples.Draw("signalpred", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
192 TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
193 TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
194 TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
195 signalpred->Add(signalpredlo,-1);
196 signalpred->Add(signalpredro);
197 puresignal->Add(signalpred,-1);//subtracting signal contamination
198 ofstream myfile;
199 myfile.open ("ShapeFit_log.txt");
200 establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
201 myfile.close();
202 */
203 }
204
205 void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts) {
206 cout << "Doing counting experiment ... " << endl;
207
208 for(int isample=0;isample<signalsamples.collection.size();isample++) {
209 cout << "Considering sample " << signalsamples.collection[isample].samplename << endl;
210 for(int ibin=0;ibin<jzbcuts.size();ibin++) {
211 cout << " Considering bin " << jzbcuts[ibin] << endl;
212 for(int isyst=0;isyst<uncertainties[isample*jzbcuts.size()+ibin].size();isyst++) {
213 if(isyst==0) {
214 if(uncertainties[isample*jzbcuts.size()+ibin][isyst]!=jzbcuts[ibin]) cout << "WATCH OUT THERE SEEMS TO BE A PROBLEM - THE TEST BIN DOESN'T CORRESPOND TO YOUR BIN!" << endl;
215 continue;
216 }
217 cout << " Considering syst " << isyst << " : " << uncertainties[isample*jzbcuts.size()+ibin][isyst] << endl;
218 }//end of syst loop
219 }//end of bin loop
220 }//end of sample loop
221 }
222
223 void susy_scan_axis_labeling(TH2F *histo) {
224 histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
225 histo->GetXaxis()->CenterTitle();
226 histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
227 histo->GetYaxis()->CenterTitle();
228 }
229
230 void scan_susy_space(string mcjzb, string datajzb) {
231 TCanvas *c3 = new TCanvas("c3","c3");
232 vector<float> binning;
233 binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
234 /*
235 binning.push_back(50);
236 binning.push_back(75);
237 binning.push_back(100);
238 binning.push_back(150);
239 binning.push_back(200);
240 binning.push_back(500);
241 */
242 float arrbinning[binning.size()];
243 for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
244 TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
245 puredata->SetMarkerSize(DataMarkerSize);
246 TH1F *allbgs = allsamples.Draw("allbgs", "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
247 TH1F *allbgsb = allsamples.Draw("allbgsb", "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
248 TH1F *allbgsc = allsamples.Draw("allbgsc", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
249 allbgs->Add(allbgsb,-1);
250 allbgs->Add(allbgsc);
251 int ndata=puredata->Integral();
252 ofstream myfile;
253 myfile.open ("susyscan_log.txt");
254 TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
255 TTree *suevents = (TTree*)susyscanfile->Get("events");
256 TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
257 TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
258 TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
259 TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
260
261 susy_scan_axis_labeling(exclusionmap);
262 susy_scan_axis_labeling(exclusionmap1s);
263 susy_scan_axis_labeling(exclusionmap2s);
264 susy_scan_axis_labeling(exclusionmap3s);
265
266 Int_t MyPalette[100];
267 Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
268 Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
269 Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
270 Double_t stop[] = {0., .25, .50, .75, 1.0};
271 Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
272 for (int i=0;i<100;i++) MyPalette[i] = FI+i;
273
274 gStyle->SetPalette(100, MyPalette);
275
276 for(int m23=50;m23<500;m23+=25) {
277 for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
278 {
279 c3->cd();
280 stringstream drawcondition;
281 drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
282 TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
283 TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
284 stringstream drawvar,drawvar2;
285 drawvar<<mcjzb<<">>puresignal";
286 drawvar2<<"-"<<mcjzb<<">>puresignall";
287 suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
288 suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
289 if(puresignal->Integral()<60) {
290 delete puresignal;
291 continue;
292 }
293 puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
294 puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
295 stringstream saveas;
296 saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
297 cout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
298 // TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
299 // TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
300 // TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
301 // signalpred->Add(signalpredlo,-1);
302 // signalpred->Add(signalpredro);
303 // puresignal->Add(signalpred,-1);//subtracting signal contamination
304 //---------------------
305 // cout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
306 // TH1F *puresignal = allsamples.Draw("puresignal",mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
307 vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
308 if(results.size()==0) {
309 delete puresignal;
310 continue;
311 }
312 exclusionmap->Fill(m23,m0,results[0]);
313 exclusionmap1s->Fill(m23,m0,results[1]);
314 exclusionmap2s->Fill(m23,m0,results[2]);
315 exclusionmap3s->Fill(m23,m0,results[3]);
316 delete puresignal;
317 cout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
318 }
319 }//end of model scan for loop
320
321 cout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
322 c3->cd();
323 exclusionmap->Draw("CONTZ");
324 CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
325 exclusionmap->Draw("COLZ");
326 CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
327
328 exclusionmap1s->Draw("CONTZ");
329 CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
330 exclusionmap1s->Draw("COLZ");
331 CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
332
333 exclusionmap2s->Draw("CONTZ");
334 CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
335 exclusionmap2s->Draw("COLZ");
336 CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
337
338 exclusionmap3s->Draw("CONTZ");
339 CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
340 exclusionmap3s->Draw("COLZ");
341 CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
342
343 TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
344 exclusionmap->Write();
345 exclusionmap1s->Write();
346 exclusionmap2s->Write();
347 exclusionmap3s->Write();
348 exclusion_limits->Close();
349 susyscanfile->Close();
350
351 myfile.close();
352 }
353
354
355
356