ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Util/src/ElectronEnergySmearingScaling.cc
Revision: 1.4
Committed: Wed Oct 24 13:14:59 2012 UTC (12 years, 6 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.3: +1 -1 lines
Log Message:
fixed a bug

File Contents

# Content
1 #include "ElectronEnergySmearingScaling.h"
2 #include "CommonDefs.h"
3 #include <TRandom.h>
4 #include <TRandom1.h>
5
6
7 void correct_electron_energy(mithep::Electron * ele)
8 {
9 Float_t sigma;
10 Float_t scale;
11
12 //the scale is z_mass_data/z_mass_mc
13 //the sigma of the gaussian used for the smearing is
14 //sqrt(2*(sigma_data^2 - sigma_mc^2)/z_mass_mc^2)
15 //if the resolution is worse in mc, then we use 0 for the smearing factor
16 if(abs(ele->Eta()) < 1.479){
17 Float_t sigma = 0;
18 scale = (-0.0357085 + Z_MASS)/(0.455985 + Z_MASS);
19 }
20 else {
21 Float_t sigma = sqrt(2)*2.15011/(1.28248+91.1876);
22 scale = (1.09664+Z_MASS)/(1.2824+Z_MASS);
23 }
24 TRandom * rand = new TRandom(1024);
25 Float_t smear = rand->Gaus(0,sigma);
26 ele->SetPtEtaPhi((scale+smear)*ele->Pt(),ele->Eta(),ele->Phi());
27
28 }
29
30 //based on
31 //http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/EgammaAnalysis/ElectronTools/src/PatElectronEnergyCalibrator.cc?revision=1.6&view=markup
32
33 std::pair<float,float> apply_claudes_electron_energy_corrections(mithep::Electron * ele, const float energy, const float energy_error, int run, std::string dataset_, const bool debug_){
34
35 Bool_t isMC_ = (run == 1);
36
37 float corr=0., scale=1.;
38 float dsigMC=0., corrMC=0.;
39
40 if (!isMC_) {
41 if (dataset_=="Jan16ReReco") { // corrections for january 16 ReReco
42 // values from http://indico.cern.ch/getFile.py/access?contribId=2&resId=0&materialId=slides&confId=176520
43 if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()<0.94) {
44 if (run>=160431 && run<=167913) corr = -0.0014;
45 if (run>=170000 && run<=172619) corr = -0.0016;
46 if (run>=172620 && run<=173692) corr = -0.0017;
47 if (run>=175830 && run<=177139) corr = -0.0021;
48 if (run>=177140 && run<=178421) corr = -0.0025;
49 if (run>=178424 && run<=180252) corr = -0.0024;
50 } else if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()>=0.94) {
51 if (run>=160431 && run<=167913) corr = 0.0059;
52 if (run>=170000 && run<=172619) corr = 0.0046;
53 if (run>=172620 && run<=173692) corr = 0.0045;
54 if (run>=175830 && run<=177139) corr = 0.0042;
55 if (run>=177140 && run<=178421) corr = 0.0038;
56 if (run>=178424 && run<=180252) corr = 0.0039;
57 } else if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()<0.94) {
58 if (run>=160431 && run<=167913) corr = -0.0045;
59 if (run>=170000 && run<=172619) corr = -0.0066;
60 if (run>=172620 && run<=173692) corr = -0.0058;
61 if (run>=175830 && run<=177139) corr = -0.0073;
62 if (run>=177140 && run<=178421) corr = -0.0075;
63 if (run>=178424 && run<=180252) corr = -0.0071;
64 } else if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()>=0.94) {
65 if (run>=160431 && run<=167913) corr = 0.0084;
66 if (run>=170000 && run<=172619) corr = 0.0063;
67 if (run>=172620 && run<=173692) corr = 0.0071;
68 if (run>=175830 && run<=177139) corr = 0.0056;
69 if (run>=177140 && run<=178421) corr = 0.0054;
70 if (run>=178424 && run<=180252) corr = 0.0058;
71 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()<0.94) {
72 if (run>=160431 && run<=167913) corr = -0.0082;
73 if (run>=170000 && run<=172619) corr = -0.0025;
74 if (run>=172620 && run<=173692) corr = -0.0035;
75 if (run>=175830 && run<=177139) corr = -0.0017;
76 if (run>=177140 && run<=178421) corr = -0.0010;
77 if (run>=178424 && run<=180252) corr = 0.0030;
78 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()>=0.94) {
79 if (run>=160431 && run<=167913) corr = -0.0033;
80 if (run>=170000 && run<=172619) corr = 0.0024;
81 if (run>=172620 && run<=173692) corr = 0.0014;
82 if (run>=175830 && run<=177139) corr = 0.0032;
83 if (run>=177140 && run<=178421) corr = 0.0040;
84 if (run>=178424 && run<=180252) corr = 0.0079;
85 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()<0.94) {
86 if (run>=160431 && run<=167913) corr = -0.0064;
87 if (run>=170000 && run<=172619) corr = -0.0046;
88 if (run>=172620 && run<=173692) corr = -0.0029;
89 if (run>=175830 && run<=177139) corr = -0.0040;
90 if (run>=177140 && run<=178421) corr = -0.0050;
91 if (run>=178424 && run<=180252) corr = -0.0059;
92 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()>=0.94) {
93 if (run>=160431 && run<=167913) corr = 0.0042;
94 if (run>=170000 && run<=172619) corr = 0.0060;
95 if (run>=172620 && run<=173692) corr = 0.0077;
96 if (run>=175830 && run<=177139) corr = 0.0067;
97 if (run>=177140 && run<=178421) corr = 0.0056;
98 if (run>=178424 && run<=180252) corr = 0.0047;
99 }
100 }
101 else if (dataset_=="2012Jul13ReReco") {
102 // values from https://twiki.cern.ch/twiki/bin/view/CMS/ECALELF
103 if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()<0.94) {
104 if (run>=190645 && run<=190781) scale = 1.0020;
105 if (run>=190782 && run<=191042) scale = 1.0079;
106 if (run>=191043 && run<=193555) scale = 0.9989;
107 if (run>=193556 && run<=194150) scale = 0.9974;
108 if (run>=194151 && run<=194532) scale = 0.9980;
109 if (run>=194533 && run<=195113) scale = 0.9983;
110 if (run>=195114 && run<=195915) scale = 0.9984;
111 if (run>=195916 && run<=198115) scale = 0.9975;
112 if (run>=198116 && run<=199803) scale = 1.0010;
113 if (run>=199804 && run<=200048) scale = 1.0021;
114 if (run>=200049 && run<=200151) scale = 1.0035;
115 if (run>=200152 && run<=200490) scale = 1.0013;
116 if (run>=200491 && run<=200531) scale = 1.0035;
117 if (run>=200532 && run<=201656) scale = 1.0017;
118 if (run>=201657 && run<=202305) scale = 1.0026;
119 if (run>=202305 && run<=203002) scale = 1.0037;
120 } else if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()>=0.94) {
121 if (run>=190645 && run<=190781) scale = 0.9980;
122 if (run>=190782 && run<=191042) scale = 1.0039;
123 if (run>=191043 && run<=193555) scale = 0.9949;
124 if (run>=193556 && run<=194150) scale = 0.9934;
125 if (run>=194151 && run<=194532) scale = 0.9940;
126 if (run>=194533 && run<=195113) scale = 0.9943;
127 if (run>=195114 && run<=195915) scale = 0.9944;
128 if (run>=195916 && run<=198115) scale = 0.9936;
129 if (run>=198116 && run<=199803) scale = 0.9970;
130 if (run>=199804 && run<=200048) scale = 0.9982;
131 if (run>=200049 && run<=200151) scale = 0.9996;
132 if (run>=200152 && run<=200490) scale = 0.9973;
133 if (run>=200491 && run<=200531) scale = 0.9995;
134 if (run>=200532 && run<=201656) scale = 0.9978;
135 if (run>=201657 && run<=202305) scale = 0.9986;
136 if (run>=202305 && run<=203002) scale = 0.9998;
137 } else if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()<0.94) {
138 if (run>=190645 && run<=190781) scale = 1.0032;
139 if (run>=190782 && run<=191042) scale = 1.0063;
140 if (run>=191043 && run<=193555) scale = 0.9998;
141 if (run>=193556 && run<=194150) scale = 0.9954;
142 if (run>=194151 && run<=194532) scale = 0.9965;
143 if (run>=194533 && run<=195113) scale = 0.9984;
144 if (run>=195114 && run<=195915) scale = 0.9977;
145 if (run>=195916 && run<=198115) scale = 0.9965;
146 if (run>=198116 && run<=199803) scale = 0.9999;
147 if (run>=199804 && run<=200048) scale = 1.0008;
148 if (run>=200049 && run<=200151) scale = 1.0017;
149 if (run>=200152 && run<=200490) scale = 1.0003;
150 if (run>=200491 && run<=200531) scale = 1.0017;
151 if (run>=200532 && run<=201656) scale = 0.9999;
152 if (run>=201657 && run<=202305) scale = 1.0003;
153 if (run>=202305 && run<=203002) scale = 1.0010;
154 } else if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()>=0.94) {
155 if (run>=190645 && run<=190781) scale = 0.9919;
156 if (run>=190782 && run<=191042) scale = 0.9951;
157 if (run>=191043 && run<=193555) scale = 0.9885;
158 if (run>=193556 && run<=194150) scale = 0.9841;
159 if (run>=194151 && run<=194532) scale = 0.9852;
160 if (run>=194533 && run<=195113) scale = 0.9872;
161 if (run>=195114 && run<=195915) scale = 0.9864;
162 if (run>=195916 && run<=198115) scale = 0.9852;
163 if (run>=198116 && run<=199803) scale = 0.9886;
164 if (run>=199804 && run<=200048) scale = 0.9895;
165 if (run>=200049 && run<=200151) scale = 0.9905;
166 if (run>=200152 && run<=200490) scale = 0.9890;
167 if (run>=200491 && run<=200531) scale = 0.9905;
168 if (run>=200532 && run<=201656) scale = 0.9887;
169 if (run>=201657 && run<=202305) scale = 0.9891;
170 if (run>=202305 && run<=203002) scale = 0.9897;
171 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()<0.94) {
172 if (run>=190645 && run<=190781) scale = 0.9945;
173 if (run>=190782 && run<=191042) scale = 0.9996;
174 if (run>=191043 && run<=193555) scale = 0.9968;
175 if (run>=193556 && run<=194150) scale = 0.9969;
176 if (run>=194151 && run<=194532) scale = 0.9986;
177 if (run>=194533 && run<=195113) scale = 1.0006;
178 if (run>=195114 && run<=195915) scale = 1.0010;
179 if (run>=195916 && run<=198115) scale = 1.0020;
180 if (run>=198116 && run<=199803) scale = 0.9963;
181 if (run>=199804 && run<=200048) scale = 0.9965;
182 if (run>=200049 && run<=200151) scale = 0.9992;
183 if (run>=200152 && run<=200490) scale = 0.9991;
184 if (run>=200491 && run<=200531) scale = 0.9995;
185 if (run>=200532 && run<=201656) scale = 0.9978;
186 if (run>=201657 && run<=202305) scale = 0.9987;
187 if (run>=202305 && run<=203002) scale = 1.0003;
188 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()>=0.94) {
189 if (run>=190645 && run<=190781) scale = 0.9881;
190 if (run>=190782 && run<=191042) scale = 0.9932;
191 if (run>=191043 && run<=193555) scale = 0.9904;
192 if (run>=193556 && run<=194150) scale = 0.9905;
193 if (run>=194151 && run<=194532) scale = 0.9922;
194 if (run>=194533 && run<=195113) scale = 0.9943;
195 if (run>=195114 && run<=195915) scale = 0.9946;
196 if (run>=195916 && run<=198115) scale = 0.9956;
197 if (run>=198116 && run<=199803) scale = 0.9899;
198 if (run>=199804 && run<=200048) scale = 0.9901;
199 if (run>=200049 && run<=200151) scale = 0.9928;
200 if (run>=200152 && run<=200490) scale = 0.9927;
201 if (run>=200491 && run<=200531) scale = 0.9931;
202 if (run>=200532 && run<=201656) scale = 0.9914;
203 if (run>=201657 && run<=202305) scale = 0.9923;
204 if (run>=202305 && run<=203002) scale = 0.9940;
205 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()<0.94) {
206 if (run>=190645 && run<=190781) scale = 0.9965;
207 if (run>=190782 && run<=191042) scale = 1.0010;
208 if (run>=191043 && run<=193555) scale = 0.9987;
209 if (run>=193556 && run<=194150) scale = 0.9988;
210 if (run>=194151 && run<=194532) scale = 0.9994;
211 if (run>=194533 && run<=195113) scale = 0.9999;
212 if (run>=195114 && run<=195915) scale = 1.0004;
213 if (run>=195916 && run<=198115) scale = 0.9992;
214 if (run>=198116 && run<=199803) scale = 1.0044;
215 if (run>=199804 && run<=200048) scale = 1.0060;
216 if (run>=200049 && run<=200151) scale = 1.0101;
217 if (run>=200152 && run<=200490) scale = 1.0073;
218 if (run>=200491 && run<=200531) scale = 1.0106;
219 if (run>=200532 && run<=201656) scale = 1.0069;
220 if (run>=201657 && run<=202305) scale = 1.0121;
221 if (run>=202305 && run<=203002) scale = 1.0144;
222 } else if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()>=0.94) {
223 if (run>=190645 && run<=190781) scale = 0.9862;
224 if (run>=190782 && run<=191042) scale = 0.9907;
225 if (run>=191043 && run<=193555) scale = 0.9884;
226 if (run>=193556 && run<=194150) scale = 0.9885;
227 if (run>=194151 && run<=194532) scale = 0.9891;
228 if (run>=194533 && run<=195113) scale = 0.9896;
229 if (run>=195114 && run<=195915) scale = 0.9900;
230 if (run>=195916 && run<=198115) scale = 0.9889;
231 if (run>=198116 && run<=199803) scale = 0.9941;
232 if (run>=199804 && run<=200048) scale = 0.9957;
233 if (run>=200049 && run<=200151) scale = 0.9999;
234 if (run>=200152 && run<=200490) scale = 0.9970;
235 if (run>=200491 && run<=200531) scale = 1.0004;
236 if (run>=200532 && run<=201656) scale = 0.9967;
237 if (run>=201657 && run<=202305) scale = 1.0018;
238 if (run>=202305 && run<=203002) scale = 1.0042;
239 }
240 }
241 else
242 assert(0);
243 }
244 if(isMC_){
245 if (dataset_=="Summer12_DR53X_HCP2012") {
246 if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()<0.94) dsigMC = 0.0103;
247 if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0090;
248 if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()<0.94) dsigMC = 0.0190;
249 if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0156;
250 if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()<0.94) dsigMC = 0.0269;
251 if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0287;
252 if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()<0.94) dsigMC = 0.0364;
253 if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0321;
254 }
255 else if (dataset_=="Fall11") { // values from https://hypernews.cern.ch/HyperNews/CMS/get/higgs2g/634.html, consistant with Jan16ReReco corrections
256 if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()<0.94) dsigMC = 0.0096;
257 if (ele->IsEB() && fabs(ele->SCluster()->Eta())<1 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0074;
258 if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()<0.94) dsigMC = 0.0196;
259 if (ele->IsEB() && fabs(ele->SCluster()->Eta())>=1 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0141;
260 if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()<0.94) dsigMC = 0.0279;
261 if (ele->IsEE() && fabs(ele->SCluster()->Eta())<2 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0268;
262 if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()<0.94) dsigMC = 0.0301;
263 if (ele->IsEE() && fabs(ele->SCluster()->Eta())>=2 && ele->SCluster()->R9()>=0.94) dsigMC = 0.0293;
264 }
265
266 else
267 assert(0);
268 }
269
270 float corrected_energy = -999;
271
272 if (!isMC_){
273 if(dataset_ == "Jan16ReReco") corrected_energy = energy/(1+corr);
274 else if (!isMC_ && dataset_ == "2012Jul13ReReco") corrected_energy = energy*scale;
275 else assert(0);
276 }
277 if (isMC_) {
278 TRandom1 * rand = new TRandom1();
279 rand->SetSeed(0);
280 corrMC = rand->Gaus(1,dsigMC);
281 // corrMC = 1+dsigMC; //use this only for synchronization
282 if (debug_) std::cout << "[ElectronEnergyCalibrator] unsmeared energy " << energy << std::endl;
283 corrected_energy = energy*corrMC;
284 if (debug_) std::cout << "[ElectronEnergyCalibrator] smeared energy " << corrected_energy << std::endl;
285 }
286
287 float corrected_energy_error = sqrt(energy_error*energy_error + dsigMC*dsigMC*corrected_energy*corrected_energy) ;
288
289 return std::pair<float,float>(corrected_energy,corrected_energy_error);
290
291 }