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
Log Message:
First comparison plotting script

File Contents

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