Jeng/scripts/comparePDF.C

#include <memory>
#include <TROOT.h>
#include <TString.h>
#include <TStyle.h>
#include <TH1.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TFile.h>
#include <TTree.h>
#include <Riostream.h>
// #include "TLorentzVector.h"
// #include <map>
// #include <vector>
// #include <iostream>
using namespace std;

void comparePDF(){
  gROOT->SetStyle("CMS");
  //   gStyle->SetOptStat("n");
  //   gStyle->SetOptTitle(1);
  //   gStyle->SetOptFit(11111);
  gStyle->SetPadTopMargin(0.11);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadLeftMargin(0.11);
  gStyle->SetPadBottomMargin(0.11);
  
  const double wttbar = 0.0081;
  const int nmembers = 45;
  TH1D *h1 = new TH1D[nmembers];
  TH1D *h2 = new TH1D[nmembers];
  TH1D *hm3 = new TH1D[nmembers-1];
  TH1D *hm3prime = new TH1D[nmembers-1];
  double ratio[2] = {0.,0.};
  double systMax1 = 0.,systMax2 = 0.;
  
  TFile *f0=TFile::Open("TopAnalysis_TTJets-madgraph_Fall08_all_all_pdf.root");
  
  TCanvas *cv = new TCanvas("cv","M3 & M3' systematics",1000,800);
  cv->Divide(2,2,0.001,0.001);
  cv->cd(1);
  for(int i=0; i < nmembers; i++){
    TString suffix="PDF_";
    suffix+=i;
    
    h1[i] = (TH1D*)f0->Get(TString("Mass/HadronicTop_mass_M3_")+suffix);
    h1[i].SetLineColor(i+1);
    h1[i].Scale(wttbar);
    
    if(i==0){
      h1[i].GetXaxis()->SetTitle("M3 [GeV/c^{2}]");
      h1[i].GetXaxis()->SetLabelFont(42);
      h1[i].GetXaxis()->SetLabelSize(0.04);
      h1[i].GetXaxis()->SetTitleFont(42);
      h1[i].GetXaxis()->SetTitleSize(0.04);
      h1[i].GetXaxis()->SetTitleOffset(1.2);
      h1[i].GetXaxis()->SetRangeUser(100.,450.);
      h1[i].GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
      h1[i].SetMaximum(50.);
      h1[i].GetYaxis()->SetLabelFont(42);
      h1[i].GetYaxis()->SetLabelSize(0.04);
      h1[i].GetYaxis()->SetTitleFont(42);
      h1[i].GetYaxis()->SetTitleSize(0.04);
      h1[i].GetYaxis()->SetTitleOffset(1.2);
      h1[i].Draw();
    }
    else
      h1[i].SetLineWidth(0.5);
      h1[i].Draw("same");
  }
  
  cv->cd(2);
  for(int i=0; i < nmembers; i++){
    TString suffix="PDF_";
    suffix+=i;
    
    h2[i] = (TH1D*)f0->Get(TString("Mass/HadronicTop_mass_M3prime_")+suffix);
    h2[i].SetLineColor(i+1);
    h2[i].Scale(wttbar);
    
    if(i==0){
      h2[i].GetXaxis()->SetTitle("M3' (#chi^{2}) [GeV/c^{2}]");
      h2[i].GetXaxis()->SetLabelFont(42);
      h2[i].GetXaxis()->SetLabelSize(0.04);
      h2[i].GetXaxis()->SetTitleFont(42);
      h2[i].GetXaxis()->SetTitleSize(0.04);
      h2[i].GetXaxis()->SetTitleOffset(1.2);
      h2[i].GetXaxis()->SetRangeUser(100.,450.);
      h2[i].GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
      h2[i].SetMaximum(100.);
      h2[i].GetYaxis()->SetLabelFont(42);
      h2[i].GetYaxis()->SetLabelSize(0.04);
      h2[i].GetYaxis()->SetTitleFont(42);
      h2[i].GetYaxis()->SetTitleSize(0.04);
      h2[i].GetYaxis()->SetTitleOffset(1.2);
      h2[i].Draw();
    }
    else
      h2[i].SetLineWidth(0.5);
      h2[i].Draw("same");
  }
  
  int nbin = h1[0].GetXaxis()->FindBin(500.)-1;
  //   cout<<"nbin = "<<nbin<<endl;
  const int nbin1 = nbin+1;
  double xMax=500.;
  double xMin=100.;
  double delta1[nbin1];
  double delta2[nbin1];
  double delta3[nbin1];
  double delta4[nbin1];
  
  for (int i=1;i<nmembers;i++){
    TString name1 = "hm3_";
    TString name2 = "hm3prime_";
    name1+=i;
    name2+=i;
    hm3[i] = new TH1D(name1,"M3 [GeV/c^{2}]",nbin,xMin,xMax);
    hm3prime[i] = new TH1D(name2,"M3'  (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
    
    for (int j=1;j<=nbin;j++){ // exclude over/underflow bins
      ratio[0] = h1[i].GetBinContent(j)/h1[0].GetBinContent(j);
      ratio[1] = h2[i].GetBinContent(j)/h2[0].GetBinContent(j);
      
      hm3[i].SetBinContent(j,ratio[0]);
      hm3prime[i].SetBinContent(j,ratio[1]);
    }
    
    cv->cd(3);
    if (i==1){
      hm3[i].GetXaxis()->SetTitle("M3 [GeV/c^{2}]");
      hm3[i].GetXaxis()->SetTitleSize(0.04);
      hm3[i].GetXaxis()->SetTitleOffset(1.2);
      hm3[i].GetXaxis()->SetRangeUser(100.,450.);
      hm3[i].GetYaxis()->SetTitle("Ratio");
      hm3[i].GetYaxis()->SetTitleSize(0.04);
      hm3[i].GetYaxis()->SetTitleOffset(1.2);
      hm3[i].SetMaximum(1.2);
      hm3[i].SetMinimum(0.8);
      hm3[i].Draw();
    }
    else{
      hm3[i].SetLineWidth(0.5);
      hm3[i].SetLineColor(i+1);
      hm3[i].Draw("same");
    }
    
    cv->cd(4);
    if (i==1){
      hm3prime[i].GetXaxis()->SetTitle("M3' (#chi^{2}) [GeV/c^{2}]");
      hm3prime[i].GetXaxis()->SetTitleSize(0.04);
      hm3prime[i].GetXaxis()->SetTitleOffset(1.2);
      hm3prime[i].GetXaxis()->SetRangeUser(100.,450.);
      hm3prime[i].GetYaxis()->SetTitle("Ratio");
      hm3prime[i].GetYaxis()->SetTitleSize(0.04);
      hm3prime[i].GetYaxis()->SetTitleOffset(1.2);
      hm3prime[i].SetMaximum(1.2);
      hm3prime[i].SetMinimum(0.8);
      hm3prime[i].Draw();
    }
    else{
      hm3prime[i].SetLineWidth(0.5);
      hm3prime[i].SetLineColor(i+1);
      hm3prime[i].Draw("same");
    }
  }
  
  TH1D* hm3SystP = new TH1D("M3SystP","M3 [GeV/c^{2}]",nbin,xMin,xMax);
  TH1D* hm3SystM = new TH1D("M3SystM","M3 [GeV/c^{2}]",nbin,xMin,xMax);
  TH1D* hm3primeSystP = new TH1D("M3primeSystP","M3'  (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);
  TH1D* hm3primeSystM = new TH1D("M3primeSystM","M3'  (#chi^{2}) [GeV/c^{2}]",nbin,xMin,xMax);

  // Calculate syst. errors
  for (int i=1;i <= nbin;i++){
    for (int j=1;j<nmembers;j++){
      delta1[i-1] += (hm3[j].GetBinContent(i)>1)?pow(hm3[j].GetBinContent(i)-1,2):0;
      delta2[i-1] += (hm3[j].GetBinContent(i)<1)?pow(hm3[j].GetBinContent(i)-1,2):0;
    }
    delta1[i-1] = sqrt(delta1[i-1]);
    hm3SystP->SetBinContent(i,1+delta1[i-1]);
    delta2[i-1] = sqrt(delta2[i-1]);
    hm3SystM->SetBinContent(i,1-delta2[i-1]);
  }
  cv->cd(3);
  hm3SystP->SetLineWidth(2);
  hm3SystP->Draw("same");
  hm3SystM->SetLineWidth(2);
  hm3SystM->Draw("same");
  
  for (int i=1;i <= nbin;i++){
    for (int j=1;j<nmembers;j++){
      delta3[i-1] += (hm3prime[j].GetBinContent(i)>1)?pow(hm3prime[j].GetBinContent(i)-1,2):0;
      delta4[i-1] += (hm3prime[j].GetBinContent(i)<1)?pow(hm3prime[j].GetBinContent(i)-1,2):0;
    }
    delta3[i-1] = sqrt(delta3[i-1]);
    hm3primeSystP->SetBinContent(i,1+delta3[i-1]);
    delta4[i-1] = sqrt(delta4[i-1]);
    hm3primeSystM->SetBinContent(i,1-delta4[i-1]);
  }
  cv->cd(4);
  hm3primeSystP->SetLineWidth(2);
  hm3primeSystP->Draw("same");
  hm3primeSystM->SetLineWidth(2);
  hm3primeSystM->Draw("same");
  
}
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
#	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	jengbou	1.2	h1[i].GetXaxis()->SetRangeUser(100.,450.);
57	jengbou	1.1	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	jengbou	1.2	h1[i].SetLineWidth(0.5);
68	jengbou	1.1	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	jengbou	1.2	h2[i].GetXaxis()->SetRangeUser(100.,450.);
88	jengbou	1.1	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	jengbou	1.2	h2[i].SetLineWidth(0.5);
99	jengbou	1.1	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	jengbou	1.2	hm3[i].GetXaxis()->SetRangeUser(100.,450.);
134	jengbou	1.1	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	jengbou	1.2	hm3[i].SetLineWidth(0.5);
143	jengbou	1.1	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	jengbou	1.2	hm3prime[i].GetXaxis()->SetRangeUser(100.,450.);
153	jengbou	1.1	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	jengbou	1.2	hm3prime[i].SetLineWidth(0.5);
162	jengbou	1.1	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			}