ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Util/src/ElectronEnergySmearingScaling.cc
Revision: 1.3
Committed: Fri Oct 19 18:17:43 2012 UTC (12 years, 6 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.2: +1 -1 lines
Log Message: 
updates for sync on 2012 data

File Contents

# User Rev Content
1 anlevin 1.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 anlevin 1.3 else if (dataset_=="2012Jul13ReReco") {
102 anlevin 1.1 // 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 anlevin 1.2 corrMC = 1+dsigMC;
282 anlevin 1.1 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     }