ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/various_assignments/QuickScan/QuickCompare.C
Revision: 1.1
Committed: Mon Mar 5 13:46:23 2012 UTC (13 years, 2 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Error occurred while calculating annotation data.
Log Message:
First comparison plotting script

File Contents

# Content
1 #include <iostream>
2 #include <sstream>
3
4 #include "TH2.h"
5 #include "TFile.h"
6 #include "TCanvas.h"
7 #include "TMath.h"
8 #include "TColor.h"
9 #include "../../Plotting/Modules/setTDRStyle.C"
10 #include "../../Plotting/Modules/GeneralToolBox.C"
11
12 using namespace std;
13
14 int titlecounter;
15 void write_this_title(string title,TVirtualPad *pad) {
16 pad->cd();
17 TText *text = write_title(title);
18 text->Draw();
19 }
20
21 void label_axis(TH2F *s) {
22 s->GetXaxis()->SetTitle("m_{0}");
23 s->GetXaxis()->CenterTitle();
24 s->GetYaxis()->SetTitleOffset(1.5);
25 s->GetYaxis()->SetTitle("m_{1/2}");
26 s->GetYaxis()->CenterTitle();
27 }
28
29 void Add(TH2F *h, TH2F *h2, int factor) {
30 for(int i=1;i<=h->GetNbinsX();i++) {
31 for(int j=1;j<=h->GetNbinsY();j++) {
32 h->SetBinContent(i,j,h->GetBinContent(i,j)+factor*h2->GetBinContent(i,j));
33 }
34 }
35 }
36
37 TH2F *Divide(TH2F *newh, TH2F *oldh) {
38 TH2F *divided = (TH2F*)newh->Clone("division");
39 for(int i=1;i<=newh->GetNbinsX();i++) {
40 for(int j=1;j<=newh->GetNbinsY();j++) {
41 if(oldh->GetBinContent(i,j)>0||oldh->GetBinContent(i,j)<0) divided->SetBinContent(i,j,(newh->GetBinContent(i,j))/oldh->GetBinContent(i,j));
42 else divided->SetBinContent(i,j,0);
43 }
44 }
45 return divided;
46 }
47
48 TH2F *Improvement(TH2F *newh, TH2F *oldh) {
49 TH2F *divided = (TH2F*)newh->Clone("division");
50 for(int i=1;i<=newh->GetNbinsX();i++) {
51 for(int j=1;j<=newh->GetNbinsY();j++) {
52 if(oldh->GetBinContent(i,j)>0) divided->SetBinContent(i,j,(newh->GetBinContent(i,j)-oldh->GetBinContent(i,j))/oldh->GetBinContent(i,j));
53 else divided->SetBinContent(i,j,0);
54 }
55 }
56 return divided;
57 }
58
59 TH2F *SQRT(TH2F *newh) {
60 TH2F *divided = (TH2F*)newh->Clone("division");
61 for(int i=1;i<=newh->GetNbinsX();i++) {
62 for(int j=1;j<=newh->GetNbinsY();j++) {
63 if(newh->GetBinContent(i,j)>0) divided->SetBinContent(i,j,TMath::Sqrt(newh->GetBinContent(i,j)));
64 else divided->SetBinContent(i,j,0);
65 }
66 }
67 return divided;
68 }
69
70 TH2F *get_ratio(string obs, string pred, TFile *f) {
71 TH2F *observed = (TH2F*)f->Get(obs.c_str());
72 TH2F *observedC = (TH2F*)observed->Clone("observedC");
73 TH2F *predicted= (TH2F*)f->Get(pred.c_str());
74 Add(observedC,predicted,-1);
75 TH2F *s = Divide(observedC,SQRT(observed));
76 s->SetStats(0);
77 s->GetZaxis()->SetRangeUser(-2,6);
78 label_axis(s);
79 return s;
80 }
81
82
83 void compare_onpeak_offpeak(float jzbcut, TFile *f, float mllcut) {
84 stringstream loadthis;
85 //efficiency_final_leptwo20to7000_mll_onpeak_jzb250
86 //onpeak
87 loadthis << "efficiency_final_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
88 TH2F *reference = (TH2F*) f->Get(loadthis.str().c_str());
89 reference->SetStats(0);
90
91 loadthis.str("");
92 loadthis << "efficiency_final_leptwo20to7000_mll_from" << mllcut << "_jzb" << jzbcut;
93 TH2F *hOP = (TH2F*)f->Get(loadthis.str().c_str());
94 hOP->SetStats(0);
95
96 TH2F *difference = (TH2F*)hOP->Clone("DIFF");
97 difference->SetTitle("Gain from going OffPeak");
98
99 Add(difference,reference,-1);
100
101 TH2F *improvement = Improvement(hOP,reference);
102 improvement->SetName("Relative gain from going offpeak");
103
104 improvement->GetZaxis()->SetRangeUser(0.0,1.5);
105
106 TCanvas *can = new TCanvas("can","can",1000,1000);
107 can->Divide(2,2);
108 can->cd(1);
109 reference->SetTitle("Situation on peak");
110 reference->Draw("COLZ");
111 can->cd(2);
112 stringstream hoptitle;
113 hoptitle << "Situation off peak (mll>" << mllcut << "&&pti>20)";
114 hOP->SetTitle(hoptitle.str().c_str());
115 hOP->Draw("COLZ");
116 can->cd(3);
117 difference->SetTitle("improvement off peak - on peak");
118 difference->Draw("COLZ");
119 can->cd(4);
120 improvement->SetTitle("relative improvement off peak - on peak");
121 improvement->Draw("COLZ");
122
123 stringstream saveas;
124 saveas << "GoingOffpeak___ComparingOnPeak_To_OffPeakPt2020Mll" << mllcut << "_atJZB_"<<jzbcut << ".png";
125 can->SaveAs(saveas.str().c_str());
126 delete can;
127 }
128
129 void compare_pt_cuts(float jzbcut, TFile*f) {
130 stringstream loadthis;
131 loadthis << "efficiency_final_leptwo20to7000_mll_from40_jzb" << jzbcut;
132 TH2F *reference = (TH2F*) f->Get(loadthis.str().c_str());
133 reference->SetStats(0);
134 TH2F *newone10 = (TH2F*)reference->Clone("newone10");
135 newone10->SetTitle("Going down to 10/10");
136 TH2F *newone15 = (TH2F*)reference->Clone("newone15");
137 newone15->SetTitle("Going down to 15/15");
138
139
140 stringstream loadthis2;
141 loadthis2 << "efficiency_final_lep15to20_mll_from40_jzb" << jzbcut;
142 TH2F *ontop15 = (TH2F*)f->Get(loadthis2.str().c_str());
143 ontop15->SetStats(0);
144 stringstream loadthis3;
145 loadthis3 << "efficiency_final_lep10to20_mll_from40_jzb" << jzbcut;
146 TH2F *ontop10 = (TH2F*)f->Get(loadthis3.str().c_str());
147 ontop10->SetStats(0);
148
149 Add(newone10,ontop10,1);
150 Add(newone15,ontop15,1);
151
152 TH2F *improvement10 = Improvement(newone10,reference);
153 TH2F *improvement15 = Improvement(newone15,reference);
154
155 improvement10->GetZaxis()->SetRangeUser(-0.5,0.5);
156 improvement15->GetZaxis()->SetRangeUser(-0.5,0.5);
157
158 TCanvas *can = new TCanvas("can","can",1000,1000);
159 can->Divide(2,2);
160 can->cd(1);
161 reference->Draw("COLZ");
162 can->cd(2);
163 ontop10->Draw("COLZ");
164 can->cd(3);
165 newone10->Draw("COLZ");
166 can->cd(4);
167 improvement10->Draw("COLZ");
168
169 stringstream saveas;
170 saveas << "Lowering_Cuts_To_1010_instead_of_2020_atJZB_"<<jzbcut << ".png";
171 can->SaveAs(saveas.str().c_str());
172 saveas.str("");
173 saveas << "Lowering_Cuts_To_1515_instead_of_2020_atJZB_"<<jzbcut << ".png";
174 can->cd(1);
175 reference->Draw("COLZ");
176 can->cd(2);
177 ontop15->Draw("COLZ");
178 can->cd(3);
179 newone15->Draw("COLZ");
180 can->cd(4);
181 improvement15->Draw("COLZ");
182 can->SaveAs(saveas.str().c_str());
183 saveas.str("");
184 saveas << "Comparison_of_improvement_1010_1515_atJZB_" << jzbcut << ".png";
185 can->cd(1);
186 improvement10->Draw("COLZ");
187 can->cd(2);
188 improvement15->Draw("COLZ");
189 can->cd(3)->Delete();
190 can->cd(4)->Delete();
191 can->SaveAs(saveas.str().c_str());
192
193 delete can;
194
195
196 cout << "hello" << endl;
197 }
198
199 void compare_stilde(float jzbcut, TFile *f) {
200 stringstream loadthis;
201 stringstream loadthis2;
202 loadthis << "efficiency_obs_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
203 loadthis2 << "efficiency_pred_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
204 TH2F *ratioONPEAK = get_ratio(loadthis.str(),loadthis2.str(),f);
205 loadthis.str("");
206 loadthis2.str("");
207 loadthis << "efficiency_obs_leptwo20to7000_mll_from50_jzb" << jzbcut;
208 loadthis2 << "efficiency_pred_leptwo20to7000_mll_from50_jzb" << jzbcut;
209 TH2F *ratioOFFPEAK = get_ratio(loadthis.str(),loadthis2.str(),f);
210
211 TH2F *difference = (TH2F*)ratioOFFPEAK->Clone("DIFF");
212 difference->SetTitle("Gain from going OffPeak");
213
214 Add(difference,ratioONPEAK,-1);
215
216 TH2F *improvement = Improvement(ratioOFFPEAK,ratioONPEAK);
217 improvement->SetName("Relative gain from going offpeak");
218
219 improvement->GetZaxis()->SetRangeUser(0.0,1.5);
220
221 TCanvas *can = new TCanvas("can","can",1000,1000);
222 can->Divide(2,2);
223 can->cd(1);
224 ratioONPEAK->SetTitle("Situation on peak");
225 ratioONPEAK->Draw("COLZ");
226 can->cd(2);
227 ratioOFFPEAK->SetTitle("Situation off peak (mll>40&&pti>20)");
228 ratioOFFPEAK->Draw("COLZ");
229 can->cd(3);
230 difference->SetTitle("improvement off peak - on peak");
231 difference->Draw("COLZ");
232 can->cd(4);
233 improvement->SetTitle("relative improvement off peak - on peak");
234 improvement->Draw("COLZ");
235
236 stringstream saveas;
237 saveas << "GoingOffpeak___Comparing_sTILDE_OnPeak_To_OffPeakPt2020Mll50_atJZB_"<<jzbcut << ".png";
238 can->SaveAs(saveas.str().c_str());
239 delete can;
240 }
241
242 void compare_obs_pred_final(float jzbcut, TFile *f, string mllexpression, string ptexpression) {
243 stringstream loadthis;
244 stringstream loadthis2;
245 loadthis << "efficiency_obs_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
246 loadthis2 << "efficiency_pred_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
247 TH2F *ratioONPEAK = get_ratio(loadthis.str(),loadthis2.str(),f);
248 loadthis.str("");
249 loadthis2.str("");
250 loadthis << "efficiency_obs_lep" << ptexpression << "_mll_" << mllexpression << "_jzb" << jzbcut;
251 loadthis2 << "efficiency_pred_lep" << ptexpression << "_mll_" << mllexpression << "_jzb" << jzbcut;
252 TH2F *ratioOFFPEAK = get_ratio(loadthis.str(),loadthis2.str(),f);
253
254 TH2F *difference = (TH2F*)ratioOFFPEAK->Clone("DIFF");
255 difference->SetTitle("Gain from going OffPeak");
256 Add(difference,ratioONPEAK,-1);
257 TH2F *improvement = Improvement(ratioOFFPEAK,ratioONPEAK);
258 improvement->SetName("Relative gain from going offpeak");
259 improvement->GetZaxis()->SetRangeUser(0.0,1.5);
260
261 TCanvas *can = new TCanvas("can","can",1600,1200);
262 can->Divide(4,3);
263 can->cd(4);
264 ratioONPEAK->SetTitle("s/(s+b) on peak");
265 ratioONPEAK->Draw("COLZ");
266 write_this_title("s/(s+b) on peak",can->cd(4));
267 can->cd(8);
268 ratioOFFPEAK->SetTitle("s/(s+b) off peak (mll>40&&pti>20)");
269 ratioOFFPEAK->Draw("COLZ");
270 write_this_title("s/(s+b) after",can->cd(8));
271 can->cd(12);
272 // improvement->SetTitle("Increase");
273 write_this_title("Increase",can->cd(12));
274 improvement->GetZaxis()->SetRangeUser(-0.3,0.3);
275 label_axis(improvement);
276 improvement->Draw("COLZ");
277 write_this_title("Increase",can->cd(12));
278
279 stringstream loadthis3;
280 loadthis.str("");
281 loadthis2.str("");
282 loadthis3.str("");
283 loadthis << "efficiency_obs_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
284 loadthis2 << "efficiency_pred_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
285 loadthis3 << "efficiency_final_leptwo20to7000_mll_onpeak_jzb" << jzbcut;
286 TH2F *obsb4 = (TH2F*)f->Get(loadthis.str().c_str());
287 TH2F *predb4 = (TH2F*)f->Get(loadthis2.str().c_str());
288 TH2F *finb4 = (TH2F*)f->Get(loadthis3.str().c_str());
289 obsb4->SetStats(0);
290 // obsb4->SetTitle("Observed, onpeak");
291 predb4->SetStats(0);
292 // predb4->SetTitle("Predicted, onpeak");
293 finb4->SetStats(0);
294 // finb4->SetTitle("obs-pred, onpeak");
295 loadthis.str("");
296 loadthis2.str("");
297 loadthis3.str("");
298 loadthis << "efficiency_obs_lep" << ptexpression << "_mll_" << mllexpression << "_jzb" << jzbcut;
299 loadthis2 << "efficiency_pred_lep" << ptexpression << "_mll_" << mllexpression << "_jzb" << jzbcut;
300 loadthis3 << "efficiency_final_lep" << ptexpression << "_mll_" << mllexpression << "_jzb" << jzbcut;
301 TH2F *obsaf = (TH2F*)f->Get(loadthis.str().c_str());
302 TH2F *predaf = (TH2F*)f->Get(loadthis2.str().c_str());
303 TH2F *finaf = (TH2F*)f->Get(loadthis3.str().c_str());
304 label_axis(obsaf);
305 label_axis(predaf);
306 label_axis(finaf);
307 label_axis(obsb4);
308 label_axis(predb4);
309 label_axis(finb4);
310
311 obsaf->SetStats(0);
312 // obsaf->SetTitle("Observed, offpeak");
313 predaf->SetStats(0);
314 // predaf->SetTitle("Predicted, offpeak");
315 finaf->SetStats(0);
316 // finaf->SetTitle("obs-pred, offpeak");
317
318 can->cd(1);
319 obsb4->Draw("COLZ");
320 write_this_title("observed (before)",can->cd(1));
321 can->cd(2);
322 predb4->Draw("COLZ");
323 write_this_title("predicted (before)",can->cd(2));
324 can->cd(3);
325 finb4->Draw("COLZ");
326 write_this_title("obs - pred (before)",can->cd(3));
327 can->cd(5);
328 obsaf->Draw("COLZ");
329 write_this_title("observed (after)",can->cd(5));
330 can->cd(6);
331 predaf->Draw("COLZ");
332 write_this_title("predicted (after)",can->cd(6));
333 can->cd(7);
334 finaf->Draw("COLZ");
335 write_this_title("obs - pred (after)",can->cd(7));
336 can->cd(9);
337 TH2F *a = Improvement(obsaf,obsb4);
338 // a->SetTitle("Increase");
339 a->GetZaxis()->SetRangeUser(-0.5,1.5);
340 label_axis(a);
341 a->Draw("COLZ");
342 write_this_title("Increase",can->cd(9));
343 can->cd(10);
344 TH2F *b = Improvement(predaf,predb4);
345 b->SetTitle("Increase");
346 b->GetZaxis()->SetRangeUser(-0.5,1.5);
347 label_axis(b);
348 b->Draw("COLZ");
349 write_this_title("Increase",can->cd(10));
350 can->cd(11);
351 TH2F *c = Improvement(finaf,finb4);
352 c->SetTitle("Increase");
353 c->GetZaxis()->SetRangeUser(-0.3,0.3);
354 label_axis(c);
355 c->Draw("COLZ");
356 write_this_title("Increase",can->cd(11));
357
358
359
360 stringstream saveas;
361 saveas << "FullOverview/GoingOffpeak___Comparing_Efficiencies_and_sTILDE_OnPeak_To_OffPeakPt" << ptexpression << "__" << mllexpression << "_atJZB_"<<jzbcut << ".png";
362 can->SaveAs(saveas.str().c_str());
363 delete can;
364 }
365
366 void compare_obs_pred_final(float jzbcut, TFile *f) {
367 vector<string> mllcuts;
368 mllcuts.push_back("from20");
369 mllcuts.push_back("from30");
370 mllcuts.push_back("from40");
371 mllcuts.push_back("from50");
372
373 vector<string> ptcuts;
374 ptcuts.push_back("two20to7000");
375 ptcuts.push_back("conf1010");
376 ptcuts.push_back("conf1510");
377 ptcuts.push_back("conf2010");
378 ptcuts.push_back("conf1515");
379 ptcuts.push_back("conf2015");
380
381 for(int i=0;i<mllcuts.size();i++) {
382 for(int j=0;j<ptcuts.size();j++) {
383 compare_obs_pred_final(jzbcut, f,mllcuts[i],ptcuts[j]);
384 }
385 }
386 }
387 void compare_for_jzb(float jzbcut, TFile *f) {
388 compare_pt_cuts(jzbcut,f);
389 compare_stilde(jzbcut,f);
390 compare_onpeak_offpeak(jzbcut,f,40.0);
391 compare_onpeak_offpeak(jzbcut,f,50.0);
392 compare_obs_pred_final(jzbcut,f);
393 }
394
395 int main() {
396 setTDRStyle(0);
397 //some refinement: nicer color gradient
398 Double_t stops[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
399 Double_t red[5] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
400 Double_t green[5] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
401 Double_t blue[5] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
402 int fi=TColor::CreateGradientColorTable(5, stops, red, green,blue, 255);
403 TFile *f = new TFile("NEWalloutputWithMll2complete.root");
404 compare_for_jzb(50,f);
405 compare_for_jzb(100,f);
406 compare_for_jzb(150,f);
407 compare_for_jzb(200,f);
408 compare_for_jzb(250,f);
409 return 0;
410 }
411