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 |
}
|