forked from h2gglobe/h2gglobe
-
Notifications
You must be signed in to change notification settings - Fork 0
/
EfficiencySmearer.cc
225 lines (188 loc) · 9.31 KB
/
EfficiencySmearer.cc
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
#include "EfficiencySmearer.h"
#include "PhotonReducedInfo.h"
#include "TRandom3.h"
#include <assert.h>
EfficiencySmearer::EfficiencySmearer(const efficiencySmearingParameters& par) : myParameters_(par), doPhoId_(false), doR9_(false)
{
rgen_ = new TRandom3(0);
name_="EfficiencySmearer_"+ par.categoryType + "_" + par.parameterSetName;
assert( ( myParameters_.categoryType == "2CatR9_EBEE" && myParameters_.n_categories == 4) ||
( myParameters_.categoryType == "3CatR9_EBEE" && myParameters_.n_categories == 6) ||
( myParameters_.categoryType == "2CatR9_EBEBm4EE" && myParameters_.n_categories == 6)
); // m4 is not treated separately as far as efficiencies are concerned
}
EfficiencySmearer::~EfficiencySmearer()
{
delete rgen_;
}
std::string EfficiencySmearer::photonCategory(PhotonReducedInfo & aPho) const
{
std::string myCategory="";
if (myParameters_.categoryType=="2CatR9_EBEE")
{
if (aPho.iDet()==1)
myCategory+="EB";
else
myCategory+="EE";
if (aPho.r9()>=0.94)
myCategory+="HighR9";
else
myCategory+="LowR9";
}
else if (myParameters_.categoryType=="3CatR9_EBEE")
{
if (aPho.iDet()==1)
myCategory+="EB";
else
myCategory+="EE";
if (aPho.r9()>=0.94)
myCategory+="HighR9";
else if (aPho.r9()>=0.90)
myCategory+="MidR9";
else
myCategory+="LowR9";
}
else if (myParameters_.categoryType=="EBEE")
{
if (aPho.iDet()==1)
myCategory+="EB";
else
myCategory+="EE";
}
else
{
std::cout << effName_ << " Unknown categorization. No category name is returned" << std::endl;
}
return myCategory;
}
bool EfficiencySmearer::smearPhoton(PhotonReducedInfo & aPho, float & weight, int run, float syst_shift) const
{
std::string category=photonCategory(aPho);
if (category == "")
{
std::cout << effName_ << " No category has been found associated with this photon. Giving Up" << std::endl;
return false;
}
/////////////////////// changing weigh of photon according to efficiencies ///////////////////////////////////////////
assert( !smearing_eff_graph_.empty() );
if( ! doPhoId_ || aPho.passId() ) {
weight = getWeight( ( aPho.energy() / cosh(aPho.caloPosition().PseudoRapidity()) ) ,category, syst_shift);
}
// since R9 is not a selection variable BUT a categorization variable
// introduce correlation between weights assigned to R9 and to !R9
if( doR9_ ) {
if ( aPho.r9()>0.94 ) weight = getWeight( ( aPho.energy() / cosh(aPho.caloPosition().PseudoRapidity()) ) ,category, syst_shift);
else weight = getWeight( ( aPho.energy() / cosh(aPho.caloPosition().PseudoRapidity()) ) ,category, -1*syst_shift);
}
return true;
}
bool EfficiencySmearer::init()
{
// if map is not empty, yuo're initilized and happy..
if( !smearing_eff_graph_.empty() ){
std:cout << "initialization of single photon efficiency smearer " << effName_ << " already done; proceed with usage. " << std::endl;
return true;
}
//otherwise, get smearing functions from file and set up map
std::cout << "\n>>>initializing one efficiency for single photon re-weighting: " << effName_ << std::endl;
// do basic sanity checks first
if( effName_.empty()){
std::cout << effName_ << "- you're initializing reweighting for efficiency but effName_ is empty" << std::endl; assert(false); }
if( myParameters_.efficiency_file.empty()){
std::cout << effName_ << "- you're initializing reweighting for efficiency: " << effName_ << " but input file with TGraphErrors is not specified; doing nothing. " << std::endl; assert(false); }
theEfficiencyFile_ = TFile::Open(myParameters_.efficiency_file.c_str());
// initialize formulas for the four categories;
std::string effTmpName; std::string photonCat; TGraphAsymmErrors *graphTmp, *graphClone;
photonCat = std::string("EBHighR9");
effTmpName = effName_+std::string("_")+photonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
//std::cout << "graphTmp: " << graphTmp << " " << effTmpName.c_str() <<std::endl;
assert(graphTmp!=0);
graphClone=(TGraphAsymmErrors*)graphTmp->Clone(); smearing_eff_graph_[photonCat]=graphClone;
photonCat = std::string("EBLowR9");
effTmpName = effName_+std::string("_")+photonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
//std::cout << "graphTmp: " << graphTmp << " " << effTmpName.c_str() <<std::endl;
assert(graphTmp!=0);
graphClone=(TGraphAsymmErrors*)graphTmp->Clone(); smearing_eff_graph_[photonCat]=graphClone;
photonCat = std::string("EBMidR9");
effTmpName = effName_+std::string("_")+photonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
//std::cout << "graphTmp: " << graphTmp << " " << effTmpName.c_str() <<std::endl;
if(graphTmp==0){
std::cout<<"when MidR9 is empty use LowR9"<<std::endl;
std::string fakePhotonCat = std::string("EBLowR9");
effTmpName = effName_+std::string("_")+fakePhotonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
}
graphClone=(TGraphAsymmErrors*)graphTmp->Clone(); smearing_eff_graph_[photonCat]=graphClone;
photonCat = std::string("EEHighR9");
effTmpName = effName_+std::string("_")+photonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
//std::cout << "graphTmp: " << graphTmp << " " << effTmpName.c_str() <<std::endl;
assert(graphTmp!=0);
graphClone=(TGraphAsymmErrors*)graphTmp->Clone(); smearing_eff_graph_[photonCat]=graphClone;
photonCat = std::string("EELowR9");
effTmpName = effName_+std::string("_")+photonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
//std::cout << "graphTmp: " << graphTmp << " " << effTmpName.c_str() <<std::endl;
assert(graphTmp!=0);
graphClone=(TGraphAsymmErrors*)graphTmp->Clone(); smearing_eff_graph_[photonCat]=graphClone;
photonCat = std::string("EEMidR9");
effTmpName = effName_+std::string("_")+photonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
//std::cout << "graphTmp: " << graphTmp << " " << effTmpName.c_str() <<std::endl;
if(graphTmp==0){
std::cout<<"when MidR9 is empty use LowR9"<<std::endl;
std::string fakePhotonCat = std::string("EELowR9");
effTmpName = effName_+std::string("_")+fakePhotonCat; graphTmp = (TGraphAsymmErrors*) theEfficiencyFile_->Get(effTmpName.c_str());
}
graphClone=(TGraphAsymmErrors*)graphTmp->Clone(); smearing_eff_graph_[photonCat]=graphClone;
theEfficiencyFile_->Close();
return true;
}
double EfficiencySmearer::getWeight(double pt, std::string theCategory, float syst_shift) const
{
std::map<std::string,TGraphAsymmErrors*>::const_iterator theIter = smearing_eff_graph_.find(theCategory);
if( theIter != smearing_eff_graph_.end() ) {
// determine the pair of bins between which you interpolate
int numPoints = ( theIter->second )->GetN();
double x, y;
int myBin = -1;
double xPrevious = -1e9;
for (int bin=0; bin<numPoints; bin++ ){
( theIter->second )->GetPoint(bin, x, y);
assert( xPrevious < x );
if(pt > x) {
myBin = bin; }
else break;
}
int binLow, binHigh; bool atBoundary(false);
if (myBin == -1) {binHigh = 0; binLow=0; atBoundary=true;}
else if (myBin == (numPoints-1)) {binHigh = numPoints-1; binLow=numPoints-1; atBoundary=true;}
else {binLow=myBin; binHigh=myBin+1;}
// if(myBin == -1) {binHigh = 0; binLow=0;}
// else if (myBin == (numPoints-1)) {binHigh = numPoints-1; binLow=numPoints-1;}
// else {binLow=myBin; binHigh=myBin+1;}
// get hold of efficiency ratio and error at either points
// low-high refer to the points ; up-down refers to the errors
double xLow, yLow; double xHigh, yHigh;
( theIter->second )->GetPoint(binLow, xLow, yLow);
( theIter->second )->GetPoint(binHigh, xHigh, yHigh);
double errLowYup = ( theIter->second )->GetErrorYhigh(binLow);
double errLowYdown = ( theIter->second )->GetErrorYlow(binLow);
double errHighYup = ( theIter->second )->GetErrorYhigh(binHigh);
double errHighYdown = ( theIter->second )->GetErrorYlow(binHigh);
double theErrorLow, theErrorHigh;
if(syst_shift>0) {theErrorLow = errLowYup; theErrorHigh = errHighYup;}
else {theErrorLow = errLowYdown; theErrorHigh = errHighYdown;}
double theWeight, theError;
if(!atBoundary) {
theWeight = yLow + (yHigh-yLow) / (xHigh-xLow) * (pt-xLow);
theError = theErrorLow + (theErrorHigh-theErrorLow) / (xHigh-xLow) * (pt-xLow);}
else // if instead you ARE at the boundaris of TGraphAsymmErrors, collapse on first or last of its points
{
if(myBin == (numPoints-1)) {
theWeight = yHigh; theError = theErrorHigh; }
else if (myBin == -1) {
theWeight = yLow; theError = theErrorLow; }
else { std::cout << effName_ << " ** you claim to be at boundaries of TGraphAsymmErrors but your not! This is a problem " << std::endl;}
}
return ( theWeight + (theError*syst_shift));
}
else { std::cout << effName_ << " - category asked: " << theCategory << " was not found - which is a problem. Returning weight 1. " << std::endl;
return 1.; }
}