-
Notifications
You must be signed in to change notification settings - Fork 1
/
ScaleFactor.h
241 lines (189 loc) · 12.3 KB
/
ScaleFactor.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
#ifndef ScaleFactor_h
#define ScaleFactor_h
#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TGraphAsymmErrors.h"
#include <iostream>
#include <map>
#include <cmath>
#include <string>
class ScaleFactor {
private:
std::map<std::string, TGraphAsymmErrors *> eff_data;
std::map<std::string, TGraphAsymmErrors *> eff_mc;
TH1D * etaBinsH;
void SetAxisBins(TGraphAsymmErrors* graph){
int NPOINTS = graph->GetN();
double AXISBINS [NPOINTS+1] = {};
for (int i=0; i<NPOINTS; i++) { AXISBINS[i] = (graph->GetX()[i] - graph->GetErrorXlow(i)); }
AXISBINS[NPOINTS] = (graph->GetX()[NPOINTS-1] + graph->GetErrorXhigh(NPOINTS-1));
graph->GetXaxis()->Set(NPOINTS, AXISBINS);
return;
}
bool check_SameBinning(TGraphAsymmErrors* graph1, TGraphAsymmErrors* graph2){
bool haveSameBins = false;
int n1 = graph1->GetXaxis()->GetNbins();
int n2 = graph2->GetXaxis()->GetNbins();
if (n1 != n2 ) {return false;}
else {
haveSameBins = true;
const int nbins = n1;
double x1, x2;
for (int i=0; i<nbins; i++){
x1 = (graph1->GetXaxis()->GetXbins())->GetArray()[i];
x2 = (graph2->GetXaxis()->GetXbins())->GetArray()[i];
haveSameBins = haveSameBins and (x1== x2) ;
}
}
return haveSameBins;
}
//std::string FindEtaLabel(double);
//int FindPtBin( std::map<std::string, TGraphAsymmErrors *>, std::string, double);
public:
ScaleFactor(){};
void init_ScaleFactor(TString inputRootFile, std::string HistoBaseName){
TFile * fileIn = new TFile(inputRootFile, "read");
if (fileIn->IsZombie() ) { std::cout << "ERROR in ScaleFactor::init_ScaleFactor(TString inputRootFile) from LepEffInterface/src/ScaleFactor.cc : File " <<inputRootFile << " does not exist. Please check. " <<std::endl; exit(1); };
etaBinsH = (TH1D*)fileIn->Get("etaBinsH");
std::string etaLabel, GraphName;
int nEtaBins = etaBinsH->GetNbinsX();
for (int iBin=0; iBin<nEtaBins; iBin++){
etaLabel = etaBinsH->GetXaxis()->GetBinLabel(iBin+1);
GraphName = HistoBaseName+etaLabel+"_Data";
eff_data[etaLabel] = (TGraphAsymmErrors*)fileIn->Get(TString(GraphName));
SetAxisBins(eff_data[etaLabel]);
GraphName = HistoBaseName+etaLabel+"_MC";
eff_mc[etaLabel] = (TGraphAsymmErrors*)fileIn->Get(TString(GraphName));
SetAxisBins(eff_mc[etaLabel]);
bool sameBinning = check_SameBinning(eff_data[etaLabel], eff_mc[etaLabel]);
if (!sameBinning) {std::cout<< "ERROR in ScaleFactor::init_ScaleFactor(TString inputRootFile) from LepEffInterface/src/ScaleFactor.cc . Can not proceed because ScaleFactor::check_SameBinning returned different pT binning for data and MC for eta label " << etaLabel << std::endl; exit(1); };
}
return;
}
void init_ScaleFactor(TString inputRootFile){
TFile * fileIn = new TFile(inputRootFile, "read");
if (fileIn->IsZombie() ) { std::cout << "ERROR in ScaleFactor::init_ScaleFactor(TString inputRootFile) from NTupleMaker/src/ScaleFactor.cc : #File " <<inputRootFile << " does not exist. Please check. " <<std::endl; exit(1); };
std::string HistoBaseName = "ZMass";
etaBinsH = (TH1D*)fileIn->Get("etaBinsH");
std::string etaLabel, GraphName;
int nEtaBins = etaBinsH->GetNbinsX();
for (int iBin=0; iBin<nEtaBins; iBin++){
etaLabel = etaBinsH->GetXaxis()->GetBinLabel(iBin+1);
GraphName = HistoBaseName+etaLabel+"_Data";
if (fileIn->GetListOfKeys()->Contains(TString(GraphName))){
eff_data[etaLabel] = (TGraphAsymmErrors*)fileIn->Get(TString(GraphName));
SetAxisBins(eff_data[etaLabel]);
}
else eff_data[etaLabel] = 0;
GraphName = HistoBaseName+etaLabel+"_MC";
if (fileIn->GetListOfKeys()->Contains(TString(GraphName))){
eff_mc[etaLabel] = (TGraphAsymmErrors*)fileIn->Get(TString(GraphName));
SetAxisBins(eff_mc[etaLabel]);
}
else eff_mc[etaLabel] =0;
if (eff_data[etaLabel] != 0 && eff_mc[etaLabel] != 0 ) {
bool sameBinning = check_SameBinning(eff_data[etaLabel], eff_mc[etaLabel]);
if (!sameBinning) {std::cout<< "ERROR in ScaleFactor::init_ScaleFactor(TString inputRootFile) from LepEffInterface/src/ScaleFactor.cc . Can not proceed because ScaleFactor::check_SameBinning returned different pT binning for data and MC for eta label " << etaLabel << std::endl; exit(1); };
}
}
return;
}
~ ScaleFactor(){};
std::string FindEtaLabel(double Eta, std::string Which){
Eta = fabs(Eta);
int binNumber = etaBinsH->GetXaxis()->FindFixBin(Eta);
std::string EtaLabel = etaBinsH->GetXaxis()->GetBinLabel(binNumber);
std::map<std::string, TGraphAsymmErrors*>::iterator it;
if (Which == "data"){
it = eff_data.find(EtaLabel);
if ( it == eff_data.end()) {
std::cout << "ERROR in ScaleFactor::get_EfficiencyData(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc : no object corresponding to eta label "<< EtaLabel << " for data " << std::endl; exit(1);
}
}
else if (Which == "mc"){
it = eff_mc.find(EtaLabel);
if (it == eff_mc.end()) {
std::cout << "ERROR in ScaleFactor::get_EfficiencyData(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc : no object corresponding to eta label "<< EtaLabel << " for MC " << std::endl; exit(1);
}
}
return EtaLabel;
}
int FindPtBin( std::map<std::string, TGraphAsymmErrors *> eff_map, std::string EtaLabel, double Pt){
int Npoints = eff_map[EtaLabel]->GetN();
double ptMAX = (eff_map[EtaLabel]->GetX()[Npoints-1])+(eff_map[EtaLabel]->GetErrorXhigh(Npoints-1));
double ptMIN = (eff_map[EtaLabel]->GetX()[0])-(eff_map[EtaLabel]->GetErrorXlow(0));
if (Pt >= ptMAX ) return Npoints;
else if (Pt < ptMIN){
std::cout<< "WARNING in ScaleFactor::get_EfficiencyData(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc: pT too low (pt = " << Pt << "), min value is " << ptMIN << ". Returned efficiency =1. Weight will be 1. " << std::endl;
return -99;}
else {return eff_map[EtaLabel]->GetXaxis()->FindFixBin(Pt);}
}
double get_EfficiencyData(double pt, double eta){
double eff;
std::string label = FindEtaLabel(eta,"data");
int ptbin = FindPtBin(eff_data, label, pt);
if (ptbin == -99){eff =1;} // if pt is underflow
else eff = eff_data[label]->GetY()[ptbin-1];
if (eff > 1.) {std::cout<< "WARNING in ScaleFactor::get_EfficiencyData(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc: Returned efficiency in data > 1. " << std::endl;}
if (eff < 0 ) {std::cout<<"WARNING in ScaleFactor::get_EfficiencyData(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc: Returned negative efficiency in data" <<std::endl;}
return eff;
}
double get_EfficiencyMC(double pt, double eta){
double eff;
std::string label = FindEtaLabel(eta,"mc");
int ptbin = FindPtBin(eff_mc, label, pt);
if (ptbin == -99){eff =1;} // if pt is underflow
else eff= eff_mc[label]->GetY()[ptbin-1];
if (eff > 1. ) {std::cout << "WARNING in ScaleFactor::get_EfficiencyMC(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc : Returned efficiency in MC > 1. " << std::endl;}
if (eff < 0 ) {std::cout<<"WARNING in ScaleFactor::get_EfficiencyMC(double pt, double eta) from LepEffIntrface/src/ScaleFactor.cc : Returned negative efficiency in MC. " <<std::endl;}
return eff;
}
double get_ScaleFactor(double pt, double eta){
double efficiency_data = get_EfficiencyData(pt, eta);
double efficiency_mc = get_EfficiencyMC(pt, eta);
double SF;
if ( efficiency_mc != 0) {SF = efficiency_data/efficiency_mc;}
else {
SF=0.; //std::cout << "WARNING in ScaleFactor::get_ScaleFactor(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc : MC efficiency = 0. Scale Factor set to 0. ";
}
return SF;
}
double get_EfficiencyDataError(double pt, double eta){
double eff_error;
std::string label = FindEtaLabel(eta,"data");
int ptbin = FindPtBin(eff_data, label, pt);
if (ptbin == -99){eff_error =0.;} // if pt is underflow
else eff_error= eff_data[label]->GetErrorYhigh(ptbin-1);
double effData = get_EfficiencyData(pt,eta);
if (eff_error > effData) eff_error = 0.5*effData;
return eff_error;
}
double get_EfficiencyMCError(double pt, double eta){
double eff_error;
std::string label = FindEtaLabel(eta,"mc");
int ptbin = FindPtBin(eff_mc, label, pt);
if (ptbin == -99){eff_error =0.;}
else eff_error= eff_mc[label]->GetErrorYhigh(ptbin-1);
double effMC = get_EfficiencyMC(pt,eta);
if (eff_error > effMC ) eff_error = 0.5*effMC;
return eff_error;
}
double get_ScaleFactorError(double pt, double eta){
double SF_error = 0.;
double effData = get_EfficiencyData(pt, eta);
double effMC = get_EfficiencyMC(pt, eta);
double errData = get_EfficiencyDataError(pt, eta);
double errMC = get_EfficiencyMCError(pt, eta);
if (errData == 0) {std::cout<<"WARNING in ScaleFactor::get_ScaleFactorError(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc: uncertainty on data point = 0, can not calculate uncertainty on scale factor. Uncertainty set to 0." << std::endl;}
if (errMC ==0) {std::cout<<"WARNING in ScaleFactor::get_ScaleFactorError(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc: uncertainty on MC = 0, can not calculate uncerttainty on scale factor. Uncertainty set to 0." << std::endl;}
if (effData ==0) {std::cout<<"WARNING in ScaleFactor::get_ScaleFactorError(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc: efficiency in data = 0, can not calculate uncertainty on scale factor. Uncertainty set to 0." << std::endl;}
if (effMC ==0) {std::cout<<"WARNING in ScaleFactor::get_ScaleFactorError(double pt, double eta) from LepEffInterface/src/ScaleFactor.cc: efficiency in MC = 0, can not calculate uncertainty on scale factor. Uncertainty set to 0." << std::endl;}
else {
SF_error = pow((errData/effData),2) + pow((errMC/effMC),2);
SF_error = pow(SF_error, 0.5)*(effData/effMC);
}
return SF_error;
}
};
#endif