ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/scripts/MuonIPScut.C
Revision: 1.1
Committed: Tue May 18 04:36:56 2010 UTC (14 years, 11 months ago) by jengbou
Content type: text/plain
Branch: MAIN
Log Message:
Add file

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