NonMCBackground/interface/fake_defs.h

#ifndef FAKE_DEFS
#define FAKE_DEFS

//
// System headers
//
#include <vector>                   // STL vector class
#include <iostream>                 // standard I/O
#include <iomanip>                  // functions to format standard I/O
#include <fstream>                  // functions for file I/O
#include <string>                   // C++ string class
#include <iomanip>
#include <iostream>
#include <sstream>                  // class for parsing strings
#include <assert.h>
#include <stdlib.h>   
#include <getopt.h>
using namespace std;

//
// ROOT headers
//
#include <TROOT.h>                  // access to gROOT, entry point to ROOT system
#include <TSystem.h>                // interface to OS
#include <TFile.h>                  // file handle class
#include <TNtuple.h> 
#include <TTree.h>                  // class to access ntuples
#include <TChain.h>                 // 
#include <TBranch.h>                // class to access branches in TTree
#include <TClonesArray.h>           // ROOT array class
#include <TCanvas.h>                // class for drawing
#include <TH1D.h>                   // 1D histograms
#include <TBenchmark.h>             // class to track macro running statistics
#include <TLorentzVector.h>         // 4-vector class
#include <TVector3.h>               // 3D vector class
#include <TRandom.h>

//
// bambu format headers
//
#include "EventHeader.h"
#include "Electron.h"
#include "PFJet.h"
#include "Muon.h"
#include "Vertex.h"
#include "PFCandidate.h"
#include "PFCandidateCol.h"
#include "TriggerMask.h"
#include "TriggerTable.h"
#include "Names.h"
#include "BaseMod.h"
#include "TH2DAsymErr.h"

//
// utility headers
//
#include "ParseArgs.h"
#include "SampleWeight.h"
// #include "GlobalDataAndFuncs.h"
#include "MuonSelection.h"
#include "ElectronSelection.h"
#include "MathUtils.h"
#include "IsolationSelectionDefs.h"
#include "IsolationSelection.h"
#include "SimpleLepton.h"
#include "CEffUser2D.h"
#include "InfoStruct.h"
#include "KinematicsStruct.h"
#include "SelectionDefs.h"
#include "SelectionFuncs.h"
#include "FOArgs.h"
#include "NtuplerBase.h"
#include "CPlot.h"

//----------------------------------------------------------------------------------------
#define VAR_LIST_FOINFO "run/i:evt/i:lumi/i:type/I:npv/I:npu/I:pt/F:eta/F:phi/F:pass/F:mass/F:l1pt/F:l2pt/F:met/F";
typedef struct { 
  UInt_t run,evt,lumi;
  Int_t type,npv,npu;
  Float_t fopt,foeta,fophi,fopass,mass,l1pt,l2pt,met;
} foinfo;

// class to hold files and graphs with FR info
class FR_struct {
public:
  TFile *mufrfile,*elefrfile;
  TH2D *hFRmu,*hFRele;
  TH2D *hFRmuErrl,*hFReleErrl;
  TH2D *hFRmuErrh,*hFReleErrh;
  CEffUser2D mufr, elefr;
  bool absEtaMu,absEtaEle;
};

void printSiEff(TH2DAsymErr *h, TH2DAsymErr::EErrType errType)
{
  const Int_t nx = h->GetNbinsX();
  const Int_t ny = h->GetNbinsY();

  for(Int_t iy=0; iy<=ny; iy++) {
    for(Int_t ix=0; ix<=nx; ix++) {
      if(ix==0 && iy==0) {
        cout << setw(11) << "";
      } else if(ix==0) {
        cout << "[" << setw(4) << h->GetYaxis()->GetBinLowEdge(iy) << "," << setw(4) << h->GetYaxis()->GetBinLowEdge(iy+1) << "]";
      } else if(iy==0) { 
        cout << "[" << setw(4) << h->GetXaxis()->GetBinLowEdge(ix) << "," << setw(4) << h->GetXaxis()->GetBinLowEdge(ix+1) << "]";
      } else {
        ios_base::fmtflags flags = cout.flags();
        cout.precision(7);
        cout << " " << setw(9) << fixed << h->GetBinError(ix,iy,errType) << " ";
        cout.flags(flags);
      }
    }
    cout << endl;
  }
}
//----------------------------------------------------------------------------------------
void setErrors(TH2DAsymErr *h, CEffUser2D *fr)
{
  assert(fr->getHerrl() && fr->getHerrl());
  for(int xbin=0; xbin<h->GetXaxis()->GetNbins()+2; xbin++) {
    for(int ybin=0; ybin<h->GetYaxis()->GetNbins()+2; ybin++) {
      int ncell = h->GetBin(xbin,ybin);
      int ncellChk = fr->getHerrl()->GetBin(xbin,ybin);
      assert(ncell==ncellChk);

      double errl = h->GetBinError(xbin,ybin,TH2DAsymErr::kStatErrLow);
      fr->getHerrl()->SetBinContent(ncell,errl);
      double errh = h->GetBinError(xbin,ybin,TH2DAsymErr::kStatErrHigh);
      fr->getHerrh()->SetBinContent(ncell,errh);
    }
  }
}
//----------------------------------------------------------------------------------------
FR_struct initFRs(TString mufile, TString elefile)
{
  FR_struct fr;
  cout << "loading fake file: " << mufile << endl;
  // kevin's format
  fr.mufrfile = TFile::Open(mufile); assert(fr.mufrfile->IsOpen());
  fr.hFRmu     = (TH2D*)fr.mufrfile->Get("frEtaPt");
  fr.hFRmuErrl = (TH2D*)fr.mufrfile->Get("errlEtaPt");
  fr.hFRmuErrh = (TH2D*)fr.mufrfile->Get("errhEtaPt");
  cout << "loading fake file: " << elefile << endl;
  fr.elefrfile = TFile::Open(elefile); assert(fr.elefrfile->IsOpen());
  fr.hFRele     = (TH2D*)fr.elefrfile->Get("frEtaPt");
  fr.hFReleErrl = (TH2D*)fr.elefrfile->Get("errlEtaPt");
  fr.hFReleErrh = (TH2D*)fr.elefrfile->Get("errhEtaPt");
  // si's format
  if(!fr.hFRmu) {
    TH2DAsymErr *hAsymMu = (TH2DAsymErr*)fr.mufrfile->Get("MuonFakeRateDenominatorV1_EtaPt");
    fr.hFRmu = (TH2D*)hAsymMu;
    fr.hFRmuErrl = new TH2D(*fr.hFRmu);     fr.hFRmuErrl->Reset();
    fr.hFRmuErrh = new TH2D(*fr.hFRmu);     fr.hFRmuErrh->Reset();
    fr.mufr.loadEff(fr.hFRmu,fr.hFRmuErrl,fr.hFRmuErrh);
    setErrors(hAsymMu,&fr.mufr);
  }
  if(!fr.hFRele) {
    TH2DAsymErr *hAsymEle = (TH2DAsymErr*)fr.elefrfile->Get("ElectronFakeRateDenominatorV1_EtaPt");
    fr.hFRele = (TH2D*)hAsymEle;
    fr.hFReleErrl = new TH2D(*fr.hFRele);     fr.hFReleErrl->Reset();
    fr.hFReleErrh = new TH2D(*fr.hFRele);     fr.hFReleErrh->Reset();
    fr.elefr.loadEff(fr.hFRele,fr.hFReleErrl,fr.hFReleErrh);
    setErrors(hAsymEle,&fr.elefr);
  }
  fr.mufr.loadEff(fr.hFRmu,fr.hFRmuErrl,fr.hFRmuErrh);
  fr.elefr.loadEff(fr.hFRele,fr.hFReleErrl,fr.hFReleErrh);

  // figure out whether to use abs eta
  if(fr.hFRmu->GetXaxis()->GetBinLowEdge(1) < 0) fr.absEtaMu = false;
  else fr.absEtaMu = true;
  if(fr.absEtaMu) cout << "\nUsing abs eta in muon fake rates: " <<  setprecision(3) << fr.hFRmu->GetXaxis()->GetBinLowEdge(1) << endl;
  else cout << "\nNot using abs eta in muon fake rates: " <<  setprecision(3) << fr.hFRmu->GetXaxis()->GetBinLowEdge(1) << endl;
  if(fr.hFRele->GetXaxis()->GetBinLowEdge(1) < 0) fr.absEtaEle = false;
  else fr.absEtaEle = true;
  if(fr.absEtaEle) cout << "\nUsing abs eta in electron fake rates: " <<  setprecision(3) << fr.hFRele->GetXaxis()->GetBinLowEdge(1) << endl;
  else cout << "\nNot using abs eta in electron fake rates: " <<  setprecision(3) << fr.hFRele->GetXaxis()->GetBinLowEdge(1) << endl;

  // fr.mufr.printEff(cout);
  // fr.elefr.printEff(cout);
  // fr.mufr.printErrLow(cout);
  // fr.elefr.printErrLow(cout);
  // fr.mufr.printErrHigh(cout);
  // fr.elefr.printErrHigh(cout);

  return fr;
}
//----------------------------------------------------------------------------------------
double get_fake_weight(TString uncert, SimpleLepton fake_lep, FR_struct fr)
{
  double rate=0;
  assert(uncert=="" || uncert=="lo" || uncert=="hi");
  double x = fake_lep.vec.Eta();
  double y = fake_lep.vec.Pt();
  if(abs(fake_lep.type) == 13) { // fake muons:
    if(fr.absEtaMu) x = fabs(x);
    rate = fr.mufr.getEff(x,y,true);
    if(uncert == "lo") rate -= fr.mufr.getErrLow(x,y,true);
    if(uncert == "hi") rate += fr.mufr.getErrHigh(x,y,true);
  } else if(abs(fake_lep.type) == 11) {
    if(fr.absEtaEle) x = fabs(x);
    rate = fr.elefr.getEff(x,y,true);
    if(uncert == "lo") rate -= fr.elefr.getErrLow(x,y,true);
    if(uncert == "hi") rate += fr.elefr.getErrHigh(x,y,true);
  } else { cout << "Bad fake type: " << fake_lep.type << endl; assert(0); }
  
  return rate/(1-rate);
}
#endif
Revision:	1.8
Committed:	Fri Jul 6 15:26:50 2012 UTC (12 years, 10 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.7:	+39 -15 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	dkralph	1.1	#ifndef FAKE_DEFS
2			#define FAKE_DEFS
3
4			//
5			// System headers
6			//
7			#include <vector> // STL vector class
8			#include <iostream> // standard I/O
9			#include <iomanip> // functions to format standard I/O
10			#include <fstream> // functions for file I/O
11			#include <string> // C++ string class
12	dkralph	1.6	#include <iomanip>
13			#include <iostream>
14	dkralph	1.1	#include <sstream> // class for parsing strings
15			#include <assert.h>
16			#include <stdlib.h>
17			#include <getopt.h>
18			using namespace std;
19
20			//
21			// ROOT headers
22			//
23			#include <TROOT.h> // access to gROOT, entry point to ROOT system
24			#include <TSystem.h> // interface to OS
25			#include <TFile.h> // file handle class
26			#include <TNtuple.h>
27			#include <TTree.h> // class to access ntuples
28			#include <TChain.h> //
29			#include <TBranch.h> // class to access branches in TTree
30			#include <TClonesArray.h> // ROOT array class
31			#include <TCanvas.h> // class for drawing
32			#include <TH1D.h> // 1D histograms
33			#include <TBenchmark.h> // class to track macro running statistics
34			#include <TLorentzVector.h> // 4-vector class
35			#include <TVector3.h> // 3D vector class
36			#include <TRandom.h>
37
38			//
39			// bambu format headers
40			//
41			#include "EventHeader.h"
42			#include "Electron.h"
43			#include "PFJet.h"
44			#include "Muon.h"
45			#include "Vertex.h"
46			#include "PFCandidate.h"
47			#include "PFCandidateCol.h"
48			#include "TriggerMask.h"
49			#include "TriggerTable.h"
50			#include "Names.h"
51			#include "BaseMod.h"
52	dkralph	1.6	#include "TH2DAsymErr.h"
53	dkralph	1.1
54			//
55			// utility headers
56			//
57			#include "ParseArgs.h"
58			#include "SampleWeight.h"
59			// #include "GlobalDataAndFuncs.h"
60			#include "MuonSelection.h"
61			#include "ElectronSelection.h"
62			#include "MathUtils.h"
63			#include "IsolationSelectionDefs.h"
64			#include "IsolationSelection.h"
65			#include "SimpleLepton.h"
66			#include "CEffUser2D.h"
67			#include "InfoStruct.h"
68	dkralph	1.5	#include "KinematicsStruct.h"
69	dkralph	1.1	#include "SelectionDefs.h"
70	dkralph	1.6	#include "SelectionFuncs.h"
71			#include "FOArgs.h"
72	dkralph	1.8	#include "NtuplerBase.h"
73	dkralph	1.6	#include "CPlot.h"
74	dkralph	1.1
75	dkralph	1.8	//----------------------------------------------------------------------------------------
76			#define VAR_LIST_FOINFO "run/i:evt/i:lumi/i:type/I:npv/I:npu/I:pt/F:eta/F:phi/F:pass/F:mass/F:l1pt/F:l2pt/F:met/F";
77			typedef struct {
78			UInt_t run,evt,lumi;
79			Int_t type,npv,npu;
80			Float_t fopt,foeta,fophi,fopass,mass,l1pt,l2pt,met;
81			} foinfo;
82
83	dkralph	1.1	// class to hold files and graphs with FR info
84			class FR_struct {
85			public:
86			TFile mufrfile,elefrfile;
87			TH2D hFRmu,hFRele;
88			TH2D hFRmuErrl,hFReleErrl;
89			TH2D hFRmuErrh,hFReleErrh;
90			CEffUser2D mufr, elefr;
91	dkralph	1.6	bool absEtaMu,absEtaEle;
92	dkralph	1.1	};
93	dkralph	1.6
94			void printSiEff(TH2DAsymErr *h, TH2DAsymErr::EErrType errType)
95			{
96			const Int_t nx = h->GetNbinsX();
97			const Int_t ny = h->GetNbinsY();
98
99			for(Int_t iy=0; iy<=ny; iy++) {
100			for(Int_t ix=0; ix<=nx; ix++) {
101			if(ix==0 && iy==0) {
102			cout << setw(11) << "";
103			} else if(ix==0) {
104			cout << "[" << setw(4) << h->GetYaxis()->GetBinLowEdge(iy) << "," << setw(4) << h->GetYaxis()->GetBinLowEdge(iy+1) << "]";
105			} else if(iy==0) {
106			cout << "[" << setw(4) << h->GetXaxis()->GetBinLowEdge(ix) << "," << setw(4) << h->GetXaxis()->GetBinLowEdge(ix+1) << "]";
107			} else {
108			ios_base::fmtflags flags = cout.flags();
109			cout.precision(7);
110			cout << " " << setw(9) << fixed << h->GetBinError(ix,iy,errType) << " ";
111			cout.flags(flags);
112			}
113			}
114			cout << endl;
115			}
116			}
117			//----------------------------------------------------------------------------------------
118			void setErrors(TH2DAsymErr h, CEffUser2D fr)
119			{
120			assert(fr->getHerrl() && fr->getHerrl());
121	dkralph	1.7	for(int xbin=0; xbin<h->GetXaxis()->GetNbins()+2; xbin++) {
122			for(int ybin=0; ybin<h->GetYaxis()->GetNbins()+2; ybin++) {
123	dkralph	1.6	int ncell = h->GetBin(xbin,ybin);
124			int ncellChk = fr->getHerrl()->GetBin(xbin,ybin);
125			assert(ncell==ncellChk);
126
127			double errl = h->GetBinError(xbin,ybin,TH2DAsymErr::kStatErrLow);
128			fr->getHerrl()->SetBinContent(ncell,errl);
129			double errh = h->GetBinError(xbin,ybin,TH2DAsymErr::kStatErrHigh);
130			fr->getHerrh()->SetBinContent(ncell,errh);
131			}
132			}
133			}
134	dkralph	1.1	//----------------------------------------------------------------------------------------
135			FR_struct initFRs(TString mufile, TString elefile)
136			{
137			FR_struct fr;
138			cout << "loading fake file: " << mufile << endl;
139	dkralph	1.6	// kevin's format
140	dkralph	1.1	fr.mufrfile = TFile::Open(mufile); assert(fr.mufrfile->IsOpen());
141			fr.hFRmu = (TH2D*)fr.mufrfile->Get("frEtaPt");
142			fr.hFRmuErrl = (TH2D*)fr.mufrfile->Get("errlEtaPt");
143			fr.hFRmuErrh = (TH2D*)fr.mufrfile->Get("errhEtaPt");
144			cout << "loading fake file: " << elefile << endl;
145			fr.elefrfile = TFile::Open(elefile); assert(fr.elefrfile->IsOpen());
146			fr.hFRele = (TH2D*)fr.elefrfile->Get("frEtaPt");
147			fr.hFReleErrl = (TH2D*)fr.elefrfile->Get("errlEtaPt");
148			fr.hFReleErrh = (TH2D*)fr.elefrfile->Get("errhEtaPt");
149	dkralph	1.6	// si's format
150			if(!fr.hFRmu) {
151			TH2DAsymErr hAsymMu = (TH2DAsymErr)fr.mufrfile->Get("MuonFakeRateDenominatorV1_EtaPt");
152			fr.hFRmu = (TH2D*)hAsymMu;
153			fr.hFRmuErrl = new TH2D(*fr.hFRmu); fr.hFRmuErrl->Reset();
154			fr.hFRmuErrh = new TH2D(*fr.hFRmu); fr.hFRmuErrh->Reset();
155			fr.mufr.loadEff(fr.hFRmu,fr.hFRmuErrl,fr.hFRmuErrh);
156			setErrors(hAsymMu,&fr.mufr);
157			}
158	dkralph	1.1	if(!fr.hFRele) {
159	dkralph	1.6	TH2DAsymErr hAsymEle = (TH2DAsymErr)fr.elefrfile->Get("ElectronFakeRateDenominatorV1_EtaPt");
160			fr.hFRele = (TH2D*)hAsymEle;
161			fr.hFReleErrl = new TH2D(*fr.hFRele); fr.hFReleErrl->Reset();
162			fr.hFReleErrh = new TH2D(*fr.hFRele); fr.hFReleErrh->Reset();
163			fr.elefr.loadEff(fr.hFRele,fr.hFReleErrl,fr.hFReleErrh);
164			setErrors(hAsymEle,&fr.elefr);
165	dkralph	1.1	}
166	dkralph	1.6	fr.mufr.loadEff(fr.hFRmu,fr.hFRmuErrl,fr.hFRmuErrh);
167	dkralph	1.1	fr.elefr.loadEff(fr.hFRele,fr.hFReleErrl,fr.hFReleErrh);
168
169	dkralph	1.6	// figure out whether to use abs eta
170			if(fr.hFRmu->GetXaxis()->GetBinLowEdge(1) < 0) fr.absEtaMu = false;
171			else fr.absEtaMu = true;
172	dkralph	1.8	if(fr.absEtaMu) cout << "\nUsing abs eta in muon fake rates: " << setprecision(3) << fr.hFRmu->GetXaxis()->GetBinLowEdge(1) << endl;
173			else cout << "\nNot using abs eta in muon fake rates: " << setprecision(3) << fr.hFRmu->GetXaxis()->GetBinLowEdge(1) << endl;
174	dkralph	1.6	if(fr.hFRele->GetXaxis()->GetBinLowEdge(1) < 0) fr.absEtaEle = false;
175			else fr.absEtaEle = true;
176	dkralph	1.8	if(fr.absEtaEle) cout << "\nUsing abs eta in electron fake rates: " << setprecision(3) << fr.hFRele->GetXaxis()->GetBinLowEdge(1) << endl;
177			else cout << "\nNot using abs eta in electron fake rates: " << setprecision(3) << fr.hFRele->GetXaxis()->GetBinLowEdge(1) << endl;
178	dkralph	1.6
179	dkralph	1.8	// fr.mufr.printEff(cout);
180			// fr.elefr.printEff(cout);
181			// fr.mufr.printErrLow(cout);
182			// fr.elefr.printErrLow(cout);
183			// fr.mufr.printErrHigh(cout);
184			// fr.elefr.printErrHigh(cout);
185	dkralph	1.3
186	dkralph	1.1	return fr;
187			}
188			//----------------------------------------------------------------------------------------
189	dkralph	1.8	double get_fake_weight(TString uncert, SimpleLepton fake_lep, FR_struct fr)
190			{
191			double rate=0;
192			assert(uncert=="" \|\| uncert=="lo" \|\| uncert=="hi");
193			double x = fake_lep.vec.Eta();
194			double y = fake_lep.vec.Pt();
195			if(abs(fake_lep.type) == 13) { // fake muons:
196			if(fr.absEtaMu) x = fabs(x);
197			rate = fr.mufr.getEff(x,y,true);
198			if(uncert == "lo") rate -= fr.mufr.getErrLow(x,y,true);
199			if(uncert == "hi") rate += fr.mufr.getErrHigh(x,y,true);
200			} else if(abs(fake_lep.type) == 11) {
201			if(fr.absEtaEle) x = fabs(x);
202			rate = fr.elefr.getEff(x,y,true);
203			if(uncert == "lo") rate -= fr.elefr.getErrLow(x,y,true);
204			if(uncert == "hi") rate += fr.elefr.getErrHigh(x,y,true);
205			} else { cout << "Bad fake type: " << fake_lep.type << endl; assert(0); }
206
207			return rate/(1-rate);
208			}
209	dkralph	1.1	#endif