SchieferD/jetcalib/mcfakerate_x.cpp

////////////////////////////////////////////////////////////////////////////////
//
// mcfakerate_x
// ------------
//
//            06/28/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
////////////////////////////////////////////////////////////////////////////////


#include "utils/cmdline.h"

#include <TROOT.h>
#include <TStyle.h>
#include <TRint.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TLegend.h>

#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include <vector>


using namespace std;


////////////////////////////////////////////////////////////////////////////////
// declare global functions
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
void setBinomialErrors(TH1F* hEff,const TH1F* hFake, const TH1F* hAll);

////////////////////////////////////////////////////////////////////////////////
// main
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
int main(int argc,char**argv)
{
  cmdline cl;
  cl.parse(argc,argv);

  vector<string> input    = cl.get_vector<string>("input");
  string         treename = cl.get_value<string> ("treename",   "t");
  double         drmax    = cl.get_value<double> ("drmax",      0.3);
  int            nptbins  = cl.get_value<int>    ("nptbins",     50);
  int            netabins = cl.get_value<int>    ("netabins",    25);
  int            nphibins = cl.get_value<int>    ("nphibins",    25);
  double         ptmin    = cl.get_value<double> ("ptmin",      0.0);
  double         ptmax    = cl.get_value<double> ("ptmax",    200.0);
  double         etamin   = cl.get_value<double> ("etamin",    -5.0);
  double         etamax   = cl.get_value<double> ("etamax",     5.0);
  double         phimin   = cl.get_value<double> ("phimin",    -3.2);
  double         phimax   = cl.get_value<double> ("phimax",     3.2);
  
  if (!cl.check()) return 0;
  cl.print();

  
  // jetalgs
  vector<string> jetalgs;
  for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
    string jetalg = *it;
    string::size_type pos;
    while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
    jetalg = jetalg.substr(0,jetalg.find(".root"));
    jetalgs.push_back(jetalg);
  }
  
  // loop over all input samples
  vector<TH1F*> ptHistos;
  vector<TH1F*> etaHistos;
  vector<TH1F*> phiHistos;
  for (unsigned int i=0;i<input.size();++i) {
    TFile* f = new TFile(input[i].c_str(),"READ");
    TTree* t = (TTree*)f->Get(treename.c_str());
    cout<<jetalgs[i]<<" sample has "<<t->GetEntries()<<" entries."<<endl;


    string namePtAll = "PtAll_"+jetalgs[i];
    string namePtFake = "PtFake_"+jetalgs[i];
    string namePtEff = "PtEff_"+jetalgs[i];
    
    TH1F* hPtAll  = new TH1F(namePtAll.c_str(),namePtAll.c_str(),
                             nptbins,ptmin,ptmax);
    TH1F* hPtFake = new TH1F(namePtFake.c_str(),namePtFake.c_str(),
                             nptbins,ptmin,ptmax);
    TH1F* hPtEff  = new TH1F(namePtEff.c_str(),namePtEff.c_str(),
                             nptbins,ptmin,ptmax);
    
    hPtAll->Sumw2();
    hPtFake->Sumw2();
    hPtEff->Sumw2();
    
    string nameEtaAll = "EtaAll_"+jetalgs[i];
    string nameEtaFake = "EtaFake_"+jetalgs[i];
    string nameEtaEff = "EtaEff_"+jetalgs[i];
    
    TH1F* hEtaAll  = new TH1F(nameEtaAll.c_str(),nameEtaAll.c_str(),
                              netabins,etamin,etamax);
    TH1F* hEtaFake = new TH1F(nameEtaFake.c_str(),nameEtaFake.c_str(),
                              netabins,etamin,etamax);
    TH1F* hEtaEff  = new TH1F(nameEtaEff.c_str(),nameEtaEff.c_str(),
                              netabins,etamin,etamax);
    
    hEtaAll->Sumw2();
    hEtaFake->Sumw2();
    hEtaEff->Sumw2();
    
    string namePhiAll = "PhiAll_"+jetalgs[i];
    string namePhiFake = "PhiFake_"+jetalgs[i];
    string namePhiEff = "PhiEff_"+jetalgs[i];
    
    TH1F* hPhiAll = new TH1F(namePhiAll.c_str(),namePhiAll.c_str(),
                             nphibins,phimin,phimax);
    TH1F* hPhiFake = new TH1F(namePhiFake.c_str(),namePhiFake.c_str(),
                              nphibins,phimin,phimax);
    TH1F* hPhiEff = new TH1F(namePhiEff.c_str(),namePhiEff.c_str(),
                             nphibins,phimin,phimax);
    
    hPhiAll->Sumw2();
    hPhiFake->Sumw2();
    hPhiEff->Sumw2();
    
    float weight;       t->SetBranchAddress("weight",&weight);
    char  njt;          t->SetBranchAddress("njt",   &njt);
    float jtpt[100];    t->SetBranchAddress("jtpt",   jtpt);
    float jteta[100];   t->SetBranchAddress("jteta",  jteta);
    float jtphi[100];   t->SetBranchAddress("jtphi",  jtphi);
    float jtgendr[100]; t->SetBranchAddress("jtgendr",jtgendr);
    
    int nevts=t->GetEntries();
    
    for (int ievt=0;ievt<nevts;ievt++) {
      t->GetEntry(ievt);
      for (int ijt=0;ijt<njt;ijt++) {
        hPtAll ->Fill(jtpt[ijt], weight);
        hEtaAll->Fill(jteta[ijt],weight);
        hPhiAll->Fill(jtphi[ijt],weight);
        
        if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) {
          hPtFake ->Fill(jtpt[ijt], weight);
          hEtaFake->Fill(jteta[ijt],weight);
          hPhiFake->Fill(jtphi[ijt],weight);
        }
      }
    }
    
    hPtEff->Divide(hPtFake,hPtAll,1.0,1.0,"B");
    setBinomialErrors(hPtEff,hPtFake,hPtAll);
    ptHistos.push_back(hPtEff);
    
    hEtaEff->Divide(hEtaFake,hEtaAll,1.0,1.0,"B");
    setBinomialErrors(hEtaEff,hEtaFake,hEtaAll);
    etaHistos.push_back(hEtaEff);
 
    hPhiEff->Divide(hPhiFake,hPhiAll,1.0,1.0,"B");
    setBinomialErrors(hPhiEff,hPhiFake,hPhiAll);
    phiHistos.push_back(hPhiEff);
  }
  
  
  // make plots
  argc=1;
  TRint* app = new TRint(argv[0],&argc,argv);
  
  gStyle->SetOptStat(0);
  
  Color_t color;
  
  // phi canvas
  TCanvas* cPhiEff  = new TCanvas("FakeVsPhi", "FakeVsPhi",830,10,400,400);
  TLegend* lPhiEff  = new TLegend(0.65,0.85-phiHistos.size()*0.05,0.86,0.85);
  cPhiEff->cd();
  lPhiEff->SetFillColor(10);
  
  color = kBlack;
  for (unsigned int i=0;i<phiHistos.size();i++) {
    TH1F* h = phiHistos[i];
    lPhiEff->AddEntry(h,jetalgs[i].c_str(),"l");
    h->SetMarkerColor(color);
    h->SetLineColor(color++);
    h->SetMarkerStyle(kFullCircle);
    h->SetMarkerSize(0.5);
    if (i==0) {
      h->SetTitle("#epsilon_{fake} vs #phi");
      h->SetXTitle("jet #phi");
      h->SetYTitle("#epsilon_{fake}");
      h->Draw("E");
      h->SetMinimum(0.0);
      h->SetMaximum(1.05);
    }
    else h->Draw("ESAME");
  }
  lPhiEff->Draw();

  // eta canvas
  TCanvas* cEtaEff  = new TCanvas("EffVsEta", "EffVsEta",420,10,400,400);
  TLegend* lEtaEff  = new TLegend(0.65,0.85-etaHistos.size()*0.05,0.86,0.85);
  cEtaEff->cd();
  lEtaEff->SetFillColor(10);
  
  color = kBlack;
  for (unsigned int i=0;i<etaHistos.size();i++) {
    TH1F* h = etaHistos[i];
    lEtaEff->AddEntry(h,jetalgs[i].c_str(),"l");
    h->SetMarkerColor(color);
    h->SetLineColor(color++);
    h->SetMarkerStyle(kFullCircle);
    h->SetMarkerSize(0.5);
    if (i==0) {
      h->SetTitle("#epsilon_{fake} vs #eta");
      h->SetXTitle("jet #eta");
      h->SetYTitle("#epsilon_{fake}");
      h->Draw("E");
      h->SetMinimum(0.0);
      h->SetMaximum(1.05);
    }
    else h->Draw("ESAME");
    
  }
  lEtaEff->Draw();

  // pt canvas
  TCanvas* cPtEff  = new TCanvas("EffVsPt", "EffVsPt",10,10,400,400);
  TLegend* lPtEff  = new TLegend(0.65,0.85-ptHistos.size()*0.05,0.86,0.85);
  cPtEff->cd();
  lPtEff->SetFillColor(10);
  
  color = kBlack;
  for (unsigned int i=0;i<ptHistos.size();i++) {
    TH1F* h = ptHistos[i];
    lPtEff->AddEntry(h,jetalgs[i].c_str(),"l");
    h->SetMarkerColor(color);
    h->SetLineColor(color++);
    h->SetMarkerStyle(kFullCircle);
    h->SetMarkerSize(0.5);
    if (i==0) {
      h->SetTitle("#epsilon_{fake} vs p_{T}");
      h->SetXTitle("jet p_{T} [GeV]");
      h->SetYTitle("#epsilon_{fake}");
      h->Draw("E");
      h->SetMinimum(0.0);
      h->SetMaximum(1.05);
    }
    else h->Draw("ESAME");
  }
  lPtEff->Draw();
  
  app->Run();
  
  return 0;
}


////////////////////////////////////////////////////////////////////////////////
// implementation of global function
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
void setBinomialErrors(TH1F* hEff,const TH1F* hFake, const TH1F* hAll)
{
  for (int i=1;i<hEff->GetNbinsX();i++) {
    double nfake = hFake->GetBinContent(i);
    double nall  = hAll ->GetBinContent(i);
    double eeff  = (nall>0.0) ? std::sqrt(nfake/(nall*nall)*(1-nfake/nall)) : 0.0;
    hEff->SetBinError(i,eeff);
  }
}
Revision:	1.1
Committed:	Wed Aug 15 17:20:06 2007 UTC (17 years, 8 months ago) by schiefer
Branch:	MAIN
Log Message:	first import
#	User	Rev	Content
1	schiefer	1.1	////////////////////////////////////////////////////////////////////////////////
2			//
3			// mcfakerate_x
4			// ------------
5			//
6			// 06/28/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
7			////////////////////////////////////////////////////////////////////////////////
8
9
10			#include "utils/cmdline.h"
11
12			#include <TROOT.h>
13			#include <TStyle.h>
14			#include <TRint.h>
15			#include <TFile.h>
16			#include <TTree.h>
17			#include <TH1F.h>
18			#include <TCanvas.h>
19			#include <TLegend.h>
20
21			#include <iostream>
22			#include <iomanip>
23			#include <cmath>
24			#include <string>
25			#include <vector>
26
27
28			using namespace std;
29
30
31			////////////////////////////////////////////////////////////////////////////////
32			// declare global functions
33			////////////////////////////////////////////////////////////////////////////////
34
35			//______________________________________________________________________________
36			void setBinomialErrors(TH1F* hEff,const TH1F* hFake, const TH1F* hAll);
37
38			////////////////////////////////////////////////////////////////////////////////
39			// main
40			////////////////////////////////////////////////////////////////////////////////
41
42			//______________________________________________________________________________
43			int main(int argc,char**argv)
44			{
45			cmdline cl;
46			cl.parse(argc,argv);
47
48			vector<string> input = cl.get_vector<string>("input");
49			string treename = cl.get_value<string> ("treename", "t");
50			double drmax = cl.get_value<double> ("drmax", 0.3);
51			int nptbins = cl.get_value<int> ("nptbins", 50);
52			int netabins = cl.get_value<int> ("netabins", 25);
53			int nphibins = cl.get_value<int> ("nphibins", 25);
54			double ptmin = cl.get_value<double> ("ptmin", 0.0);
55			double ptmax = cl.get_value<double> ("ptmax", 200.0);
56			double etamin = cl.get_value<double> ("etamin", -5.0);
57			double etamax = cl.get_value<double> ("etamax", 5.0);
58			double phimin = cl.get_value<double> ("phimin", -3.2);
59			double phimax = cl.get_value<double> ("phimax", 3.2);
60
61			if (!cl.check()) return 0;
62			cl.print();
63
64
65			// jetalgs
66			vector<string> jetalgs;
67			for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
68			string jetalg = *it;
69			string::size_type pos;
70			while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
71			jetalg = jetalg.substr(0,jetalg.find(".root"));
72			jetalgs.push_back(jetalg);
73			}
74
75			// loop over all input samples
76			vector<TH1F*> ptHistos;
77			vector<TH1F*> etaHistos;
78			vector<TH1F*> phiHistos;
79			for (unsigned int i=0;i<input.size();++i) {
80			TFile* f = new TFile(input[i].c_str(),"READ");
81			TTree* t = (TTree*)f->Get(treename.c_str());
82			cout<<jetalgs[i]<<" sample has "<<t->GetEntries()<<" entries."<<endl;
83
84
85			string namePtAll = "PtAll_"+jetalgs[i];
86			string namePtFake = "PtFake_"+jetalgs[i];
87			string namePtEff = "PtEff_"+jetalgs[i];
88
89			TH1F* hPtAll = new TH1F(namePtAll.c_str(),namePtAll.c_str(),
90			nptbins,ptmin,ptmax);
91			TH1F* hPtFake = new TH1F(namePtFake.c_str(),namePtFake.c_str(),
92			nptbins,ptmin,ptmax);
93			TH1F* hPtEff = new TH1F(namePtEff.c_str(),namePtEff.c_str(),
94			nptbins,ptmin,ptmax);
95
96			hPtAll->Sumw2();
97			hPtFake->Sumw2();
98			hPtEff->Sumw2();
99
100			string nameEtaAll = "EtaAll_"+jetalgs[i];
101			string nameEtaFake = "EtaFake_"+jetalgs[i];
102			string nameEtaEff = "EtaEff_"+jetalgs[i];
103
104			TH1F* hEtaAll = new TH1F(nameEtaAll.c_str(),nameEtaAll.c_str(),
105			netabins,etamin,etamax);
106			TH1F* hEtaFake = new TH1F(nameEtaFake.c_str(),nameEtaFake.c_str(),
107			netabins,etamin,etamax);
108			TH1F* hEtaEff = new TH1F(nameEtaEff.c_str(),nameEtaEff.c_str(),
109			netabins,etamin,etamax);
110
111			hEtaAll->Sumw2();
112			hEtaFake->Sumw2();
113			hEtaEff->Sumw2();
114
115			string namePhiAll = "PhiAll_"+jetalgs[i];
116			string namePhiFake = "PhiFake_"+jetalgs[i];
117			string namePhiEff = "PhiEff_"+jetalgs[i];
118
119			TH1F* hPhiAll = new TH1F(namePhiAll.c_str(),namePhiAll.c_str(),
120			nphibins,phimin,phimax);
121			TH1F* hPhiFake = new TH1F(namePhiFake.c_str(),namePhiFake.c_str(),
122			nphibins,phimin,phimax);
123			TH1F* hPhiEff = new TH1F(namePhiEff.c_str(),namePhiEff.c_str(),
124			nphibins,phimin,phimax);
125
126			hPhiAll->Sumw2();
127			hPhiFake->Sumw2();
128			hPhiEff->Sumw2();
129
130			float weight; t->SetBranchAddress("weight",&weight);
131			char njt; t->SetBranchAddress("njt", &njt);
132			float jtpt[100]; t->SetBranchAddress("jtpt", jtpt);
133			float jteta[100]; t->SetBranchAddress("jteta", jteta);
134			float jtphi[100]; t->SetBranchAddress("jtphi", jtphi);
135			float jtgendr[100]; t->SetBranchAddress("jtgendr",jtgendr);
136
137			int nevts=t->GetEntries();
138
139			for (int ievt=0;ievt<nevts;ievt++) {
140			t->GetEntry(ievt);
141			for (int ijt=0;ijt<njt;ijt++) {
142			hPtAll ->Fill(jtpt[ijt], weight);
143			hEtaAll->Fill(jteta[ijt],weight);
144			hPhiAll->Fill(jtphi[ijt],weight);
145
146			if (jtgendr[ijt]<0.0\|\|jtgendr[ijt]>drmax) {
147			hPtFake ->Fill(jtpt[ijt], weight);
148			hEtaFake->Fill(jteta[ijt],weight);
149			hPhiFake->Fill(jtphi[ijt],weight);
150			}
151			}
152			}
153
154			hPtEff->Divide(hPtFake,hPtAll,1.0,1.0,"B");
155			setBinomialErrors(hPtEff,hPtFake,hPtAll);
156			ptHistos.push_back(hPtEff);
157
158			hEtaEff->Divide(hEtaFake,hEtaAll,1.0,1.0,"B");
159			setBinomialErrors(hEtaEff,hEtaFake,hEtaAll);
160			etaHistos.push_back(hEtaEff);
161
162			hPhiEff->Divide(hPhiFake,hPhiAll,1.0,1.0,"B");
163			setBinomialErrors(hPhiEff,hPhiFake,hPhiAll);
164			phiHistos.push_back(hPhiEff);
165			}
166
167
168			// make plots
169			argc=1;
170			TRint* app = new TRint(argv[0],&argc,argv);
171
172			gStyle->SetOptStat(0);
173
174			Color_t color;
175
176			// phi canvas
177			TCanvas* cPhiEff = new TCanvas("FakeVsPhi", "FakeVsPhi",830,10,400,400);
178			TLegend* lPhiEff = new TLegend(0.65,0.85-phiHistos.size()*0.05,0.86,0.85);
179			cPhiEff->cd();
180			lPhiEff->SetFillColor(10);
181
182			color = kBlack;
183			for (unsigned int i=0;i<phiHistos.size();i++) {
184			TH1F* h = phiHistos[i];
185			lPhiEff->AddEntry(h,jetalgs[i].c_str(),"l");
186			h->SetMarkerColor(color);
187			h->SetLineColor(color++);
188			h->SetMarkerStyle(kFullCircle);
189			h->SetMarkerSize(0.5);
190			if (i==0) {
191			h->SetTitle("#epsilon_{fake} vs #phi");
192			h->SetXTitle("jet #phi");
193			h->SetYTitle("#epsilon_{fake}");
194			h->Draw("E");
195			h->SetMinimum(0.0);
196			h->SetMaximum(1.05);
197			}
198			else h->Draw("ESAME");
199			}
200			lPhiEff->Draw();
201
202			// eta canvas
203			TCanvas* cEtaEff = new TCanvas("EffVsEta", "EffVsEta",420,10,400,400);
204			TLegend* lEtaEff = new TLegend(0.65,0.85-etaHistos.size()*0.05,0.86,0.85);
205			cEtaEff->cd();
206			lEtaEff->SetFillColor(10);
207
208			color = kBlack;
209			for (unsigned int i=0;i<etaHistos.size();i++) {
210			TH1F* h = etaHistos[i];
211			lEtaEff->AddEntry(h,jetalgs[i].c_str(),"l");
212			h->SetMarkerColor(color);
213			h->SetLineColor(color++);
214			h->SetMarkerStyle(kFullCircle);
215			h->SetMarkerSize(0.5);
216			if (i==0) {
217			h->SetTitle("#epsilon_{fake} vs #eta");
218			h->SetXTitle("jet #eta");
219			h->SetYTitle("#epsilon_{fake}");
220			h->Draw("E");
221			h->SetMinimum(0.0);
222			h->SetMaximum(1.05);
223			}
224			else h->Draw("ESAME");
225
226			}
227			lEtaEff->Draw();
228
229			// pt canvas
230			TCanvas* cPtEff = new TCanvas("EffVsPt", "EffVsPt",10,10,400,400);
231			TLegend* lPtEff = new TLegend(0.65,0.85-ptHistos.size()*0.05,0.86,0.85);
232			cPtEff->cd();
233			lPtEff->SetFillColor(10);
234
235			color = kBlack;
236			for (unsigned int i=0;i<ptHistos.size();i++) {
237			TH1F* h = ptHistos[i];
238			lPtEff->AddEntry(h,jetalgs[i].c_str(),"l");
239			h->SetMarkerColor(color);
240			h->SetLineColor(color++);
241			h->SetMarkerStyle(kFullCircle);
242			h->SetMarkerSize(0.5);
243			if (i==0) {
244			h->SetTitle("#epsilon_{fake} vs p_{T}");
245			h->SetXTitle("jet p_{T} [GeV]");
246			h->SetYTitle("#epsilon_{fake}");
247			h->Draw("E");
248			h->SetMinimum(0.0);
249			h->SetMaximum(1.05);
250			}
251			else h->Draw("ESAME");
252			}
253			lPtEff->Draw();
254
255			app->Run();
256
257			return 0;
258			}
259
260
261			////////////////////////////////////////////////////////////////////////////////
262			// implementation of global function
263			////////////////////////////////////////////////////////////////////////////////
264
265			//______________________________________________________________________________
266			void setBinomialErrors(TH1F* hEff,const TH1F* hFake, const TH1F* hAll)
267			{
268			for (int i=1;i<hEff->GetNbinsX();i++) {
269			double nfake = hFake->GetBinContent(i);
270			double nall = hAll ->GetBinContent(i);
271			double eeff = (nall>0.0) ? std::sqrt(nfake/(nallnall)(1-nfake/nall)) : 0.0;
272			hEff->SetBinError(i,eeff);
273			}
274			}