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

# User Rev Content
1 buchmann 1.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