ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/scripts/comparePDF.C
Revision: 1.2
Committed: Sat May 23 19:35:08 2009 UTC (15 years, 11 months ago) by jengbou
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +8 -2 lines
Log Message:
update

File Contents

# Content
1 #include <memory>
2 #include <TROOT.h>
3 #include <TString.h>
4 #include <TStyle.h>
5 #include <TH1.h>
6 #include <TCanvas.h>
7 #include <TLegend.h>
8 #include <TFile.h>
9 #include <TTree.h>
10 #include <Riostream.h>
11 // #include "TLorentzVector.h"
12 // #include <map>
13 // #include <vector>
14 // #include <iostream>
15 using namespace std;
16
17 void comparePDF(){
18 gROOT->SetStyle("CMS");
19 // gStyle->SetOptStat("n");
20 // gStyle->SetOptTitle(1);
21 // gStyle->SetOptFit(11111);
22 gStyle->SetPadTopMargin(0.11);
23 gStyle->SetPadRightMargin(0.05);
24 gStyle->SetPadLeftMargin(0.11);
25 gStyle->SetPadBottomMargin(0.11);
26
27 const double wttbar = 0.0081;
28 const int nmembers = 45;
29 TH1D *h1 = new TH1D[nmembers];
30 TH1D *h2 = new TH1D[nmembers];
31 TH1D *hm3 = new TH1D[nmembers-1];
32 TH1D *hm3prime = new TH1D[nmembers-1];
33 double ratio[2] = {0.,0.};
34 double systMax1 = 0.,systMax2 = 0.;
35
36 TFile *f0=TFile::Open("TopAnalysis_TTJets-madgraph_Fall08_all_all_pdf.root");
37
38 TCanvas *cv = new TCanvas("cv","M3 & M3' systematics",1000,800);
39 cv->Divide(2,2,0.001,0.001);
40 cv->cd(1);
41 for(int i=0; i < nmembers; i++){
42 TString suffix="PDF_";
43 suffix+=i;
44
45 h1[i] = (TH1D*)f0->Get(TString("Mass/HadronicTop_mass_M3_")+suffix);
46 h1[i].SetLineColor(i+1);
47 h1[i].Scale(wttbar);
48
49 if(i==0){
50 h1[i].GetXaxis()->SetTitle("M3 [GeV/c^{2}]");
51 h1[i].GetXaxis()->SetLabelFont(42);
52 h1[i].GetXaxis()->SetLabelSize(0.04);
53 h1[i].GetXaxis()->SetTitleFont(42);
54 h1[i].GetXaxis()->SetTitleSize(0.04);
55 h1[i].GetXaxis()->SetTitleOffset(1.2);
56 h1[i].GetXaxis()->SetRangeUser(100.,450.);
57 h1[i].GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
58 h1[i].SetMaximum(50.);
59 h1[i].GetYaxis()->SetLabelFont(42);
60 h1[i].GetYaxis()->SetLabelSize(0.04);
61 h1[i].GetYaxis()->SetTitleFont(42);
62 h1[i].GetYaxis()->SetTitleSize(0.04);
63 h1[i].GetYaxis()->SetTitleOffset(1.2);
64 h1[i].Draw();
65 }
66 else
67 h1[i].SetLineWidth(0.5);
68 h1[i].Draw("same");
69 }
70
71 cv->cd(2);
72 for(int i=0; i < nmembers; i++){
73 TString suffix="PDF_";
74 suffix+=i;
75
76 h2[i] = (TH1D*)f0->Get(TString("Mass/HadronicTop_mass_M3prime_")+suffix);
77 h2[i].SetLineColor(i+1);
78 h2[i].Scale(wttbar);
79
80 if(i==0){
81 h2[i].GetXaxis()->SetTitle("M3' (#chi^{2}) [GeV/c^{2}]");
82 h2[i].GetXaxis()->SetLabelFont(42);
83 h2[i].GetXaxis()->SetLabelSize(0.04);
84 h2[i].GetXaxis()->SetTitleFont(42);
85 h2[i].GetXaxis()->SetTitleSize(0.04);
86 h2[i].GetXaxis()->SetTitleOffset(1.2);
87 h2[i].GetXaxis()->SetRangeUser(100.,450.);
88 h2[i].GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
89 h2[i].SetMaximum(100.);
90 h2[i].GetYaxis()->SetLabelFont(42);
91 h2[i].GetYaxis()->SetLabelSize(0.04);
92 h2[i].GetYaxis()->SetTitleFont(42);
93 h2[i].GetYaxis()->SetTitleSize(0.04);
94 h2[i].GetYaxis()->SetTitleOffset(1.2);
95 h2[i].Draw();
96 }
97 else
98 h2[i].SetLineWidth(0.5);
99 h2[i].Draw("same");
100 }
101
102 int nbin = h1[0].GetXaxis()->FindBin(500.)-1;
103 // cout<<"nbin = "<<nbin<<endl;
104 const int nbin1 = nbin+1;
105 double xMax=500.;
106 double xMin=100.;
107 double delta1[nbin1];
108 double delta2[nbin1];
109 double delta3[nbin1];
110 double delta4[nbin1];
111
112 for (int i=1;i<nmembers;i++){
113 TString name1 = "hm3_";
114 TString name2 = "hm3prime_";
115 name1+=i;
116 name2+=i;
117 hm3[i] = new TH1D(name1,"M3 [GeV/c^{2}]",nbin,xMin,xMax);
118 hm3prime[i] = new TH1D(name2,"M3' (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
119
120 for (int j=1;j<=nbin;j++){ // exclude over/underflow bins
121 ratio[0] = h1[i].GetBinContent(j)/h1[0].GetBinContent(j);
122 ratio[1] = h2[i].GetBinContent(j)/h2[0].GetBinContent(j);
123
124 hm3[i].SetBinContent(j,ratio[0]);
125 hm3prime[i].SetBinContent(j,ratio[1]);
126 }
127
128 cv->cd(3);
129 if (i==1){
130 hm3[i].GetXaxis()->SetTitle("M3 [GeV/c^{2}]");
131 hm3[i].GetXaxis()->SetTitleSize(0.04);
132 hm3[i].GetXaxis()->SetTitleOffset(1.2);
133 hm3[i].GetXaxis()->SetRangeUser(100.,450.);
134 hm3[i].GetYaxis()->SetTitle("Ratio");
135 hm3[i].GetYaxis()->SetTitleSize(0.04);
136 hm3[i].GetYaxis()->SetTitleOffset(1.2);
137 hm3[i].SetMaximum(1.2);
138 hm3[i].SetMinimum(0.8);
139 hm3[i].Draw();
140 }
141 else{
142 hm3[i].SetLineWidth(0.5);
143 hm3[i].SetLineColor(i+1);
144 hm3[i].Draw("same");
145 }
146
147 cv->cd(4);
148 if (i==1){
149 hm3prime[i].GetXaxis()->SetTitle("M3' (#chi^{2}) [GeV/c^{2}]");
150 hm3prime[i].GetXaxis()->SetTitleSize(0.04);
151 hm3prime[i].GetXaxis()->SetTitleOffset(1.2);
152 hm3prime[i].GetXaxis()->SetRangeUser(100.,450.);
153 hm3prime[i].GetYaxis()->SetTitle("Ratio");
154 hm3prime[i].GetYaxis()->SetTitleSize(0.04);
155 hm3prime[i].GetYaxis()->SetTitleOffset(1.2);
156 hm3prime[i].SetMaximum(1.2);
157 hm3prime[i].SetMinimum(0.8);
158 hm3prime[i].Draw();
159 }
160 else{
161 hm3prime[i].SetLineWidth(0.5);
162 hm3prime[i].SetLineColor(i+1);
163 hm3prime[i].Draw("same");
164 }
165 }
166
167 TH1D* hm3SystP = new TH1D("M3SystP","M3 [GeV/c^{2}]",nbin,xMin,xMax);
168 TH1D* hm3SystM = new TH1D("M3SystM","M3 [GeV/c^{2}]",nbin,xMin,xMax);
169 TH1D* hm3primeSystP = new TH1D("M3primeSystP","M3' (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
170 TH1D* hm3primeSystM = new TH1D("M3primeSystM","M3' (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
171
172 // Calculate syst. errors
173 for (int i=1;i <= nbin;i++){
174 for (int j=1;j<nmembers;j++){
175 delta1[i-1] += (hm3[j].GetBinContent(i)>1)?pow(hm3[j].GetBinContent(i)-1,2):0;
176 delta2[i-1] += (hm3[j].GetBinContent(i)<1)?pow(hm3[j].GetBinContent(i)-1,2):0;
177 }
178 delta1[i-1] = sqrt(delta1[i-1]);
179 hm3SystP->SetBinContent(i,1+delta1[i-1]);
180 delta2[i-1] = sqrt(delta2[i-1]);
181 hm3SystM->SetBinContent(i,1-delta2[i-1]);
182 }
183 cv->cd(3);
184 hm3SystP->SetLineWidth(2);
185 hm3SystP->Draw("same");
186 hm3SystM->SetLineWidth(2);
187 hm3SystM->Draw("same");
188
189 for (int i=1;i <= nbin;i++){
190 for (int j=1;j<nmembers;j++){
191 delta3[i-1] += (hm3prime[j].GetBinContent(i)>1)?pow(hm3prime[j].GetBinContent(i)-1,2):0;
192 delta4[i-1] += (hm3prime[j].GetBinContent(i)<1)?pow(hm3prime[j].GetBinContent(i)-1,2):0;
193 }
194 delta3[i-1] = sqrt(delta3[i-1]);
195 hm3primeSystP->SetBinContent(i,1+delta3[i-1]);
196 delta4[i-1] = sqrt(delta4[i-1]);
197 hm3primeSystM->SetBinContent(i,1-delta4[i-1]);
198 }
199 cv->cd(4);
200 hm3primeSystP->SetLineWidth(2);
201 hm3primeSystP->Draw("same");
202 hm3primeSystM->SetLineWidth(2);
203 hm3primeSystM->Draw("same");
204
205 }