ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/scripts/MuonIPScut.C
Revision: 1.2
Committed: Thu May 20 22:06:34 2010 UTC (14 years, 11 months ago) by jengbou
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +19 -17 lines
Log Message:
update

File Contents

# User Rev Content
1 jengbou 1.1 #include "TROOT.h"
2     #include "TH1.h"
3     #include "TF1.h"
4     #include "TMath.h"
5     #include "TAxis.h"
6     #include "TString.h"
7     #include "TFile.h"
8     #include "TCanvas.h"
9     #include "TLegend.h"
10     #include "map.h"
11    
12     std::map<TString,TString> histName;
13     std::map<TString,TH1*> hs;
14 jengbou 1.2 std::map<TString,TCanvas*> ts;
15 jengbou 1.1
16     void MuonIPScut()
17     {
18     gROOT->SetStyle("CMS");
19     //gStyle->SetOptFit(1111);
20 jengbou 1.2 TFile *fsig = TFile::Open("TTbarJets_madgraph_Spring10_5_all.root");
21     TFile *fbkg = TFile::Open("InclusiveMu15_Spring10_5_all.root");
22 jengbou 1.1
23     // LOAD HISTOGRAMS FROM FILES
24     const int nDim = 2;
25     const int nRef = 2;
26     const int nSel = 3;
27     const int nVar = 3;
28     double weight[2] = {0.001305,0.1255};
29     string dim[nDim] = {"2D","3D"};
30     string ref[nRef] = {"bs","pv"};
31     string sel[nSel] = {"cut0","cut1","cut2"};
32     string sel1[nSel] = {"excl0","excl1","excl2"};
33     string varName[3] = {"IPS","IP","IPerr"};
34    
35     for (int i = 0; i < nDim; i++) {
36     for (int j = 0; j < nRef; j++) {
37     for (int k = 0; k < nSel; k++) {
38    
39     for (int vi = 0; vi < nVar; vi++) {
40     TString histComm = TString(varName[vi] + "_" + dim[i] + "_" + ref[j] + "_" + sel[k]);
41     TString histIncl = histComm + "_incl";
42     TString histTagIncl = TString(dim[i] + "_" + ref[j] + "_" + sel[k]);
43     histName[histIncl] = TString("Muons/muon_" + varName[vi] + "_" + histTagIncl);
44     TString histExcl = histComm + "_excl";
45     TString histTagExcl = TString(dim[i] + "_" + ref[j] + "_" + sel1[k]);
46     histName[histExcl] = TString("Muons/muon_" + varName[vi] + "_" + histTagExcl);
47    
48     cout << "histComm = " << histComm << endl;
49     cout << "histName[histIncl] = " << histName[histIncl] << endl;
50     cout << "histName[histExcl] = " << histName[histExcl] << endl;
51    
52     TString hTT = histIncl + "_TT";
53     hs[hTT] = (TH1F*) fsig->Get(histName[histIncl]);
54     TString hQCD = histIncl + "_QCD";
55     hs[hQCD] = (TH1F*) fbkg->Get(histName[histIncl]);
56     TString hQCDbc = histExcl + "_QCDbc";
57     hs[hQCDbc] = (TH1F*) fbkg->Get(histName[histExcl]);
58    
59     // For weighted distributions
60     hs[histComm + "_sig"] = (TH1F*) hs[hTT]->Clone();
61     hs[histComm + "_bkg"] = (TH1F*) hs[hQCD]->Clone(); // QCD other decays
62     hs[histComm + "_bkgbc"] = (TH1F*) hs[hQCDbc]->Clone(); //QCD b/c decays
63    
64     //ttbar
65     hs[histComm + "_sig"]->Scale(weight[0]);
66     // QCD
67     hs[histComm + "_bkg"]->Scale(weight[1]);
68     hs[histComm + "_bkgbc"]->Scale(weight[1]);
69    
70     // Normalized to area 1.
71     hs[hTT]->Scale(1./hs[hTT]->Integral());
72     hs[hQCD]->Scale(1./hs[hQCD]->Integral());
73    
74     // Calculating S/B
75     hs["SoverB_"+histComm] = new TH1F(TString("SoverB_"+histComm),"S/B",
76     hs[histComm+"_sig"]->GetNbinsX(),
77     hs[histComm+"_sig"]->GetXaxis()->GetXmin(),
78     hs[histComm+"_sig"]->GetXaxis()->GetXmax());
79     for (int bi = 1; bi <= hs[histComm+"_sig"]->GetNbinsX(); bi++) {
80     double ss = hs[histComm+"_sig"]->Integral(hs[histComm+"_sig"]->GetXaxis()->GetXmin(),bi);
81     double bb = hs[histComm+"_bkg"]->Integral(hs[histComm+"_sig"]->GetXaxis()->GetXmin(),bi)
82     + hs[histComm+"_bkgbc"]->Integral(hs[histComm+"_sig"]->GetXaxis()->GetXmin(),bi);
83     double val = ((bb > 0) ? (ss/bb) : 0.);
84     hs["SoverB_"+histComm]->SetBinContent(bi,val);
85     }
86    
87     // Total weighted bkgshis
88     hs[histComm+"_QCDtot"] = new TH1F(TString(histComm + "_QCDtot"),"B1+B2",
89     hs[histComm+"_bkg"]->GetNbinsX(),hs[histComm+"_bkg"]->GetXaxis()->GetXmin(),
90     hs[histComm+"_bkg"]->GetXaxis()->GetXmax());
91     hs[histComm+"_QCDtot"]->Add(hs[histComm+"_bkg"],hs[histComm+"_bkgbc"],1.,1.);
92     // Normalized to area 1
93     hs[hQCDbc]->Scale((hs[hQCDbc]->GetEntries()/hs[histComm+"_QCDtot"]->GetEntries()*1.)/hs[hQCDbc]->Integral());
94     hs[histComm+"_QCDtot"]->Scale(1./hs[histComm+"_QCDtot"]->Integral());
95    
96    
97     // weighted sig plus bkg
98     hs["SplusB_"+histComm] = new TH1F(TString("SplusB_"+histComm),"S+B",
99     hs[histComm+"_sig"]->GetNbinsX(),
100     hs[histComm+"_sig"]->GetXaxis()->GetXmin(),
101     hs[histComm+"_sig"]->GetXaxis()->GetXmax());
102     hs["SplusB_"+histComm]->Add(hs[histComm+"_sig"],hs[histComm+"_bkg"],1.,1.);
103     hs["SplusB_"+histComm]->Add(hs[histComm+"_bkgbc"]);
104    
105    
106     // Weighted
107     hs["SoverSqrtSplusB_"+histComm] = new TH1F(TString("SoverSqrtSplusB_"+histComm),
108     "#frac{S}{#sqrt{S+B}}",
109     hs[histComm+"_sig"]->GetNbinsX(),
110     hs[histComm+"_sig"]->GetXaxis()->GetXmin(),
111     hs[histComm+"_sig"]->GetXaxis()->GetXmax());
112     double maxCut = 0.;
113     double maxSig = 0.;
114    
115     for (int bi = 1; bi <= hs[histComm+"_sig"]->GetNbinsX(); bi++) {
116     double ss = hs[histComm+"_sig"]->Integral(hs[histComm+"_sig"]->GetXaxis()->GetXmin(),bi);
117     double bb = hs[histComm+"_bkg"]->Integral(hs[histComm+"_sig"]->GetXaxis()->GetXmin(),bi)
118     + hs[histComm+"_bkgbc"]->Integral(hs[histComm+"_sig"]->GetXaxis()->GetXmin(),bi);
119     double val = ((ss+bb > 0) ? ss/sqrt(ss+bb) : 0.);
120     if (val > maxSig) {
121     maxSig = val;
122     maxCut = hs[histComm+"_sig"]->GetXaxis()->GetXmin() +
123     1.*bi*(hs[histComm+"_sig"]->GetXaxis()->GetXmax() -
124     hs[histComm+"_sig"]->GetXaxis()->GetXmin())/hs[histComm+"_sig"]->GetNbinsX();
125     }
126     hs["SoverSqrtSplusB_"+histComm]->SetBinContent(bi,val);
127     }
128     std::cout << "Max " << varName[vi] << " = " << maxSig << " ; cut at: " << maxCut << std::endl;
129     }
130     // }
131     // }
132     // }
133    
134     // // Plot
135     // for (int i = 0; i < nDim; i++) {
136     // for (int j = 0; j < nRef; j++) {
137     // for (int k = 0; k < nSel; k++) {
138    
139     TString CommonTag = TString("IPS_" + dim[i] + "_" + ref[j] + "_" + sel[k]);
140     TString CommonTag1 = TString("IP_" + dim[i] + "_" + ref[j] + "_" + sel[k]);
141     //cout << "CommonTag = " << CommonTag << endl;
142 jengbou 1.2 TString fileName = TString("Summery_"+dim[i]+"_"+ref[j]+"_"+sel[k]);
143     ts[fileName] = new TCanvas(fileName,fileName,800,800);
144     ts[fileName]->Divide(2,3);
145 jengbou 1.1
146 jengbou 1.2 ts[fileName]->cd(1);
147 jengbou 1.1 hs["SplusB_"+CommonTag]->SetLineColor(3);
148 jengbou 1.2 hs["SplusB_"+CommonTag]->GetXaxis()->SetTitle(TString(dim[i]+" Impact Parameter Significance w.r.t. "+ref[j]+"-"+sel[k]));
149 jengbou 1.1 hs["SplusB_"+CommonTag]->GetYaxis()->SetTitle("Normalized to 5/pb");
150     // hs["SplusB_"+CommonTag]->GetYaxis()->SetTitle("Num of Events @ L_{int} = 5 pb^{-1}");
151     if (ref[j] == "pv") hs["SplusB_"+CommonTag]->GetXaxis()->SetRangeUser(0,6);
152     hs["SplusB_"+CommonTag]->SetMinimum(0);
153    
154     hs["SplusB_"+CommonTag]->Draw("hist");
155     hs[CommonTag+"_bkg"]->SetLineColor(9);
156     hs[CommonTag+"_bkg"]->SetLineWidth(2);
157     hs[CommonTag+"_bkg"]->Draw("histsame");
158     hs[CommonTag+"_bkgbc"]->SetLineColor(6);
159     hs[CommonTag+"_bkgbc"]->Draw("histsame");
160     hs[CommonTag+"_bkgbc"]->SetLineWidth(2);
161     hs[CommonTag+"_sig"]->SetLineColor(2);
162     hs[CommonTag+"_sig"]->SetLineWidth(2);
163     hs[CommonTag+"_sig"]->Draw("histsame");
164    
165     TLegend *l0 = new TLegend(0.45,0.4,0.8,0.7);
166     l0->SetFillColor(0);
167     l0->AddEntry(hs["SplusB_"+CommonTag],"S+B_{1}+B_{2}","l");
168     l0->AddEntry(hs[CommonTag+"_sig"],"t#bar{t} (S)","l");
169     l0->AddEntry(hs[CommonTag+"_bkgbc"],"QCD - b/c decays (B_{1})","l");
170     l0->AddEntry(hs[CommonTag+"_bkg"],"QCD - others (B_{2})","l");
171     l0->Draw();
172    
173 jengbou 1.2 ts[fileName]->cd(2);
174 jengbou 1.1 hs[CommonTag+"_incl_TT"]->SetLineColor(2);
175     hs[CommonTag+"_incl_TT"]->SetLineWidth(2);
176 jengbou 1.2 hs[CommonTag+"_incl_TT"]->GetXaxis()->SetTitle(TString(dim[i]+" Impact Parameter Significance w.r.t. "+ref[j]+"-"+sel[k]));
177 jengbou 1.1 hs[CommonTag+"_incl_TT"]->GetYaxis()->SetTitle("Arbitrary Units");
178     if (ref[j] == "pv") hs[CommonTag+"_incl_TT"]->GetXaxis()->SetRangeUser(0,6);
179    
180     hs[CommonTag+"_incl_TT"]->SetLabelFont(42,"xyz");
181     hs[CommonTag+"_incl_TT"]->SetTitleFont(42,"xyz");
182     hs[CommonTag+"_incl_TT"]->SetLabelSize(0.032,"x");
183     hs[CommonTag+"_incl_TT"]->SetLabelSize(0.035,"xyz");//0.035
184     hs[CommonTag+"_incl_TT"]->SetTitleSize(0.05,"xyz");
185     hs[CommonTag+"_incl_TT"]->SetTitleOffset(1.,"xy");
186    
187     hs[CommonTag+"_incl_TT"]->Draw("hist");
188    
189     // hs[CommonTag+"_incl_TT"]->Fit("gaus","","QRM",0.,10.);
190     // TF1 *f0 = (TF1*)hs[CommonTag+"_incl_TT"]->GetFunction("gaus");
191     // f0->SetLineColor(7);
192     // f0->Draw("same");
193     // hs[CommonTag+"_QCD"]->SetLineColor(4);
194     // hs[CommonTag+"_QCD"]->Draw("histsame");
195     hs[CommonTag+"_QCDtot"]->SetLineColor(4);
196     hs[CommonTag+"_QCDtot"]->Draw("histsame");
197     hs[CommonTag+"_excl_QCDbc"]->SetLineColor(6);
198     hs[CommonTag+"_excl_QCDbc"]->SetLineWidth(2);
199     hs[CommonTag+"_excl_QCDbc"]->Draw("histsame");
200     TLegend *l1 = new TLegend(0.5,0.5,0.85,0.7);
201     l1->SetFillColor(0);
202     l1->AddEntry(hs[CommonTag+"_incl_TT"],"t#bar{t}","l");
203     // l1->AddEntry(hs[CommonTag+"_QCD"],"QCD","l");
204     l1->AddEntry(hs[CommonTag+"_QCDtot"],"QCD all","l");
205     l1->AddEntry(hs[CommonTag+"_excl_QCDbc"],"QCD - b/c decays","l");
206     l1->Draw();
207    
208 jengbou 1.2 ts[fileName]->cd(3);
209 jengbou 1.1 hs["SplusB_"+CommonTag1]->SetLineColor(3);
210 jengbou 1.2 hs["SplusB_"+CommonTag1]->GetXaxis()->SetTitle(TString(dim[i]+" Impact Parameter w.r.t. "+ref[j]+"-"+sel[k]));
211 jengbou 1.1 hs["SplusB_"+CommonTag1]->GetYaxis()->SetTitle("Normalized to 5/pb");
212     // hs["SplusB"+CommonTag1]->GetYaxis()->SetTitle("Num of Events @ L_{int} = 5 pb^{-1}");
213     if (ref[j] == "pv") hs["SplusB_"+CommonTag1]->GetXaxis()->SetRangeUser(0,0.04);
214     hs["SplusB_"+CommonTag1]->Draw("hist");
215     hs[CommonTag1+"_bkg"]->SetLineColor(9);
216     hs[CommonTag1+"_bkg"]->SetLineWidth(2);
217     hs[CommonTag1+"_bkg"]->Draw("histsame");
218     hs[CommonTag1+"_bkgbc"]->SetLineColor(6);
219     hs[CommonTag1+"_bkgbc"]->Draw("histsame");
220     hs[CommonTag1+"_bkgbc"]->SetLineWidth(2);
221     hs[CommonTag1+"_sig"]->SetLineColor(2);
222     hs[CommonTag1+"_sig"]->SetLineWidth(2);
223     hs[CommonTag1+"_sig"]->Draw("histsame");
224    
225     TLegend *l01 = new TLegend(0.45,0.4,0.8,0.7);
226     l01->SetFillColor(0);
227     l01->AddEntry(hs["SplusB_"+CommonTag1],"S+B_{1}+B_{2}","l");
228     l01->AddEntry(hs[CommonTag1+"_sig"],"t#bar{t} (S)","l");
229     l01->AddEntry(hs[CommonTag1+"_bkgbc"],"QCD - b/c decays (B_{1})","l");
230     l01->AddEntry(hs[CommonTag1+"_bkg"],"QCD - others (B_{2})","l");
231     l01->Draw();
232    
233 jengbou 1.2 ts[fileName]->cd(4);
234 jengbou 1.1 hs[CommonTag1+"_incl_TT"]->SetLineColor(2);
235     hs[CommonTag1+"_incl_TT"]->SetLineWidth(2);
236 jengbou 1.2 hs[CommonTag1+"_incl_TT"]->GetXaxis()->SetTitle(TString(dim[i]+" Impact Parameter w.r.t. "+ref[j]+"-"+sel[k]));
237 jengbou 1.1 hs[CommonTag1+"_incl_TT"]->GetYaxis()->SetTitle("Arbitrary Units");
238     if (ref[j] == "pv") hs[CommonTag1+"_incl_TT"]->GetXaxis()->SetRangeUser(0,0.04);
239    
240     hs[CommonTag1+"_incl_TT"]->SetLabelFont(42,"xyz");
241     hs[CommonTag1+"_incl_TT"]->SetTitleFont(42,"xyz");
242     hs[CommonTag1+"_incl_TT"]->SetLabelSize(0.032,"x");
243     hs[CommonTag1+"_incl_TT"]->SetLabelSize(0.035,"xyz");//0.035
244     hs[CommonTag1+"_incl_TT"]->SetTitleSize(0.05,"xyz");
245     hs[CommonTag1+"_incl_TT"]->SetTitleOffset(1.,"xy");
246    
247     hs[CommonTag1+"_incl_TT"]->Draw("hist");
248    
249     // hs[CommonTag1+"_incl_TT"]->Fit("gaus","","QRM",0.,10.);
250     // TF1 *f0 = (TF1*)hs[CommonTag1+"_incl_TT"]->GetFunction("gaus");
251     // f0->SetLineColor(7);
252     // f0->Draw("same");
253     // hs[CommonTag1+"_incl_QCD"]->SetLineColor(4);
254     // hs[CommonTag1+"_incl_QCD"]->Draw("histsame");
255     hs[CommonTag1+"_QCDtot"]->SetLineColor(4);
256     hs[CommonTag1+"_QCDtot"]->Draw("histsame");
257     hs[CommonTag1+"_excl_QCDbc"]->SetLineColor(6);
258     hs[CommonTag1+"_excl_QCDbc"]->SetLineWidth(2);
259     hs[CommonTag1+"_excl_QCDbc"]->Draw("histsame");
260     TLegend *l11 = new TLegend(0.5,0.5,0.85,0.7);
261     l11->SetFillColor(0);
262     l11->AddEntry(hs[CommonTag1+"_incl_TT"],"t#bar{t}","l");
263     // l11->AddEntry(hs[CommonTag1+"_incl_QCD"],"QCD","l");
264     l11->AddEntry(hs[CommonTag1+"_QCDtot"],"QCD all","l");
265     l11->AddEntry(hs[CommonTag1+"_excl_QCDbc"],"QCD - b/c decays","l");
266     l11->Draw();
267    
268 jengbou 1.2 ts[fileName]->cd(5);
269 jengbou 1.1
270     hs["SoverB_"+CommonTag]->SetLineColor(4);
271 jengbou 1.2 hs["SoverB_"+CommonTag]->GetXaxis()->SetTitle(TString(dim[i]+" Impact Parameter Significance w.r.t. "+ref[j]+"-"+sel[k]));
272 jengbou 1.1 hs["SoverB_"+CommonTag]->GetYaxis()->SetTitle("Arbitrary Units");
273     if (ref[j] == "pv") hs["SoverB_"+CommonTag]->GetXaxis()->SetRangeUser(0,6);
274     hs["SoverB_"+CommonTag]->Draw("hist");
275    
276     TLegend *l2 = new TLegend(0.6,0.6,0.68,0.68);
277     l2->SetFillColor(0);
278     l2->AddEntry(hs["SoverB_IPS_"+CommonTag],"#frac{S}{B}","");
279     l2->Draw();
280    
281 jengbou 1.2 ts[fileName]->cd(6);
282 jengbou 1.1 hs["SoverSqrtSplusB_"+CommonTag]->SetLineColor(4);
283 jengbou 1.2 hs["SoverSqrtSplusB_"+CommonTag]->GetXaxis()->SetTitle(TString(dim[i]+" Impact Parameter Significance w.r.t. "+ref[j]+"-"+sel[k]));
284 jengbou 1.1 hs["SoverSqrtSplusB_"+CommonTag]->GetYaxis()->SetTitle("Arbitrary Units");
285     if (ref[j] == "pv") hs["SoverSqrtSplusB_"+CommonTag]->GetXaxis()->SetRangeUser(0,6);
286     hs["SoverSqrtSplusB_"+CommonTag]->Draw("hist");
287    
288     TLegend *l3 = new TLegend(0.5,0.4,0.65,0.5);
289     l3->SetFillColor(0);
290     l3->AddEntry(hs["SoverSqrtSplusB_"+CommonTag],"#frac{S}{#sqrt{S+B}}","");
291     l3->Draw();
292    
293 jengbou 1.2 ts[fileName]->Print(TString("Summary_" + dim[i] + "_" + ref[j] + "_" + sel[k])+".png");
294 jengbou 1.1 //hs.clear();
295     }
296     }
297     }
298     }