ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/scripts/comparePDF.C
Revision: 1.1
Committed: Sat May 23 08:01:08 2009 UTC (15 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 <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    
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].Draw("same");
68     }
69    
70     cv->cd(2);
71     for(int i=0; i < nmembers; i++){
72     TString suffix="PDF_";
73     suffix+=i;
74    
75     h2[i] = (TH1D*)f0->Get(TString("Mass/HadronicTop_mass_M3prime_")+suffix);
76     h2[i].SetLineColor(i+1);
77     h2[i].Scale(wttbar);
78    
79     if(i==0){
80     h2[i].GetXaxis()->SetTitle("M3' (#chi^{2}) [GeV/c^{2}]");
81     h2[i].GetXaxis()->SetLabelFont(42);
82     h2[i].GetXaxis()->SetLabelSize(0.04);
83     h2[i].GetXaxis()->SetTitleFont(42);
84     h2[i].GetXaxis()->SetTitleSize(0.04);
85     h2[i].GetXaxis()->SetTitleOffset(1.2);
86    
87     h2[i].GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
88     h2[i].SetMaximum(100.);
89     h2[i].GetYaxis()->SetLabelFont(42);
90     h2[i].GetYaxis()->SetLabelSize(0.04);
91     h2[i].GetYaxis()->SetTitleFont(42);
92     h2[i].GetYaxis()->SetTitleSize(0.04);
93     h2[i].GetYaxis()->SetTitleOffset(1.2);
94     h2[i].Draw();
95     }
96     else
97     h2[i].Draw("same");
98     }
99    
100     int nbin = h1[0].GetXaxis()->FindBin(500.)-1;
101     // cout<<"nbin = "<<nbin<<endl;
102     const int nbin1 = nbin+1;
103     double xMax=500.;
104     double xMin=100.;
105     double delta1[nbin1];
106     double delta2[nbin1];
107     double delta3[nbin1];
108     double delta4[nbin1];
109    
110     for (int i=1;i<nmembers;i++){
111     TString name1 = "hm3_";
112     TString name2 = "hm3prime_";
113     name1+=i;
114     name2+=i;
115     hm3[i] = new TH1D(name1,"M3 [GeV/c^{2}]",nbin,xMin,xMax);
116     hm3prime[i] = new TH1D(name2,"M3' (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
117    
118     for (int j=1;j<=nbin;j++){ // exclude over/underflow bins
119     ratio[0] = h1[i].GetBinContent(j)/h1[0].GetBinContent(j);
120     ratio[1] = h2[i].GetBinContent(j)/h2[0].GetBinContent(j);
121    
122     hm3[i].SetBinContent(j,ratio[0]);
123     hm3prime[i].SetBinContent(j,ratio[1]);
124     }
125    
126     cv->cd(3);
127     if (i==1){
128     hm3[i].GetXaxis()->SetTitle("M3 [GeV/c^{2}]");
129     hm3[i].GetXaxis()->SetTitleSize(0.04);
130     hm3[i].GetXaxis()->SetTitleOffset(1.2);
131     hm3[i].GetYaxis()->SetTitle("Ratio");
132     hm3[i].GetYaxis()->SetTitleSize(0.04);
133     hm3[i].GetYaxis()->SetTitleOffset(1.2);
134     hm3[i].SetMaximum(1.2);
135     hm3[i].SetMinimum(0.8);
136     hm3[i].Draw();
137     }
138     else{
139     hm3[i].SetLineColor(i+1);
140     hm3[i].Draw("same");
141     }
142    
143     cv->cd(4);
144     if (i==1){
145     hm3prime[i].GetXaxis()->SetTitle("M3' (#chi^{2}) [GeV/c^{2}]");
146     hm3prime[i].GetXaxis()->SetTitleSize(0.04);
147     hm3prime[i].GetXaxis()->SetTitleOffset(1.2);
148     hm3prime[i].GetYaxis()->SetTitle("Ratio");
149     hm3prime[i].GetYaxis()->SetTitleSize(0.04);
150     hm3prime[i].GetYaxis()->SetTitleOffset(1.2);
151     hm3prime[i].SetMaximum(1.2);
152     hm3prime[i].SetMinimum(0.8);
153     hm3prime[i].Draw();
154     }
155     else{
156     hm3prime[i].SetLineColor(i+1);
157     hm3prime[i].Draw("same");
158     }
159     }
160    
161     TH1D* hm3SystP = new TH1D("M3SystP","M3 [GeV/c^{2}]",nbin,xMin,xMax);
162     TH1D* hm3SystM = new TH1D("M3SystM","M3 [GeV/c^{2}]",nbin,xMin,xMax);
163     TH1D* hm3primeSystP = new TH1D("M3primeSystP","M3' (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
164     TH1D* hm3primeSystM = new TH1D("M3primeSystM","M3' (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
165    
166     // Calculate syst. errors
167     for (int i=1;i <= nbin;i++){
168     for (int j=1;j<nmembers;j++){
169     delta1[i-1] += (hm3[j].GetBinContent(i)>1)?pow(hm3[j].GetBinContent(i)-1,2):0;
170     delta2[i-1] += (hm3[j].GetBinContent(i)<1)?pow(hm3[j].GetBinContent(i)-1,2):0;
171     }
172     delta1[i-1] = sqrt(delta1[i-1]);
173     hm3SystP->SetBinContent(i,1+delta1[i-1]);
174     delta2[i-1] = sqrt(delta2[i-1]);
175     hm3SystM->SetBinContent(i,1-delta2[i-1]);
176     }
177     cv->cd(3);
178     hm3SystP->SetLineWidth(2);
179     hm3SystP->Draw("same");
180     hm3SystM->SetLineWidth(2);
181     hm3SystM->Draw("same");
182    
183     for (int i=1;i <= nbin;i++){
184     for (int j=1;j<nmembers;j++){
185     delta3[i-1] += (hm3prime[j].GetBinContent(i)>1)?pow(hm3prime[j].GetBinContent(i)-1,2):0;
186     delta4[i-1] += (hm3prime[j].GetBinContent(i)<1)?pow(hm3prime[j].GetBinContent(i)-1,2):0;
187     }
188     delta3[i-1] = sqrt(delta3[i-1]);
189     hm3primeSystP->SetBinContent(i,1+delta3[i-1]);
190     delta4[i-1] = sqrt(delta4[i-1]);
191     hm3primeSystM->SetBinContent(i,1-delta4[i-1]);
192     }
193     cv->cd(4);
194     hm3primeSystP->SetLineWidth(2);
195     hm3primeSystP->Draw("same");
196     hm3primeSystM->SetLineWidth(2);
197     hm3primeSystM->Draw("same");
198    
199     }