ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbb/interface/VHbbNameSpace.h
Revision: 1.14
Committed: Tue Apr 16 16:20:30 2013 UTC (12 years ago) by nmohr
Content type: text/plain
Branch: MAIN
CVS Tags: lhcp_UnblindFix, hcp_Unblind, HEAD
Changes since 1.13: +61 -8 lines
Log Message:
New electron corrections

File Contents

# Content
1 #include "TLorentzVector.h"
2 #include "TVector3.h"
3 #include "TVector2.h"
4 #include "TMath.h"
5 /*#if !defined(__CINT__) && !defined(__MAKECINT__)
6 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
7 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
8 #endif*/
9
10 namespace VHbb {
11
12
13 double deltaPhi(double phi1,double phi2)
14 {
15 double result = phi1 - phi2;
16 while (result > TMath::Pi()) result -= 2*TMath::Pi();
17 while (result <= -TMath::Pi()) result += 2*TMath::Pi();
18 return result;
19 }
20
21 inline double deltaR(double eta1,double phi1,double eta2,double phi2)
22 {
23 double deta = eta1 - eta2;
24 double dphi = deltaPhi(phi1, phi2);
25 return TMath::Sqrt(deta*deta + dphi*dphi);
26 }
27
28
29 double Hmass( double V_eta,double V_phi,double V_pt,
30 double hJet1_eta,double hJet1_phi,double hJet1_pt,
31 double hJet2_eta,double hJet2_phi,double hJet2_pt ){
32
33 TVector3 V(1,1,1);
34 V.SetPtEtaPhi(V_pt,V_eta,V_phi);
35
36 TVector3 H1(1,1,1);
37 H1.SetPtEtaPhi(hJet1_pt,hJet1_eta,hJet1_phi);
38 H1.SetMag(1/sin(H1.Theta()));
39
40 TVector3 H2(1,1,1);
41 H2.SetPtEtaPhi(hJet2_pt,hJet2_eta,hJet2_phi);
42 H2.SetMag(1/sin(H2.Theta()));
43
44 TVector3 n1(H1);
45 TVector3 n2(H2);
46
47 double det= n1.Px() * n2.Py() - n2.Px() * n1.Py();
48
49 H1.SetMag( ( - n2.Py() * V.Px() + n2.Px() * V.Py() ) / (sin(n1.Theta()) *det ) );
50 H2.SetMag( ( + n1.Py() * V.Px() - n1.Px() * V.Py() ) / (sin(n2.Theta()) *det ) );
51
52 double mass=TMath::Sqrt( TMath::Power( (H1.Mag()+H2.Mag()),2 ) - TMath::Power(( ( H1+H2 ).Mag()),2) );
53
54 return mass;
55
56 }
57
58 double Hmass_comb(double hJet1_eta,double hJet1_phi,double hJet1_pt, double hJet1_mass,
59 double hJet2_eta,double hJet2_phi,double hJet2_pt, double hJet2_mass){
60
61 TLorentzVector H1, H2;
62 H1.SetPtEtaPhiM(hJet1_pt,hJet1_eta,hJet1_phi, hJet1_mass);;
63 H2.SetPtEtaPhiM(hJet2_pt,hJet2_eta,hJet2_phi, hJet2_mass);
64
65 return (H1 + H2).M();
66
67 }
68
69 double Hmass_3j(double h_eta,double h_phi,double h_pt, double h_mass,
70 double aJet_eta,double aJet_phi,double aJet_pt, double aJet_mass){
71
72 TLorentzVector H, H3;
73 H.SetPtEtaPhiM( h_pt,h_eta,h_phi, h_mass);;
74 H3.SetPtEtaPhiM(aJet_pt,aJet_eta,aJet_phi, aJet_mass);
75
76 return (H + H3).M();
77
78
79 }
80
81 double ANGLELZ(double pt, double eta, double phi, double mass, double pt2, double eta2, double phi2, double mass2){
82 TLorentzVector m1, m2, msum;
83 m1.SetPtEtaPhiM(pt, eta, phi, mass);
84 m2.SetPtEtaPhiM(pt2, eta2, phi2, mass2);
85 msum = m1 + m2;
86
87 TVector3 bZ = msum.BoostVector();
88
89 m1.Boost(-bZ);
90 m2.Boost(-bZ);
91
92 TVector3 b1;
93
94
95 if((int) (pt) % 2 == 0)
96 b1 = m1.BoostVector();
97 else
98 b1 = m2.BoostVector();
99
100 double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
101 return(cosTheta);
102 }
103
104
105 double ANGLEHB(double pt, double eta, double phi, double e, double pt2, double eta2, double phi2, double e2){
106 TLorentzVector m1, m2, msum;
107 m1.SetPtEtaPhiE(pt, eta, phi, e);
108 m2.SetPtEtaPhiE(pt2, eta2, phi2, e2);
109 msum = m1 + m2;
110
111 TVector3 bZ = msum.BoostVector();
112
113 m1.Boost(-bZ);
114 m2.Boost(-bZ);
115
116 TVector3 b1;
117
118 if((int) (pt) % 2 == 0)
119 b1 = m1.BoostVector();
120 else
121 b1 = m2.BoostVector();
122
123 double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
124 return(cosTheta);
125 }
126
127 double metCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
128 {
129 double metx = met * cos(metphi);
130 double mety = met * sin(metphi);
131 double px = 0.0, py = 0.0;
132 if (EVENT_run!=1) {
133 //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
134 px = +1.68804e-01 + 3.37139e-01*Nvtx;
135 py = -1.72555e-01 - 1.79594e-01*Nvtx;
136 } else {
137 //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
138 px = +2.22335e-02 - 6.59183e-02*Nvtx;
139 py = +1.52720e-01 - 1.28052e-01*Nvtx;
140 }
141 metx -= px;
142 mety -= py;
143 return std::sqrt(metx*metx + mety*mety);
144 }
145
146 double metphiCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
147 {
148 double metx = met * cos(metphi);
149 double mety = met * sin(metphi);
150 double px = 0.0, py = 0.0;
151 if (EVENT_run!=1) {
152
153 //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
154 px = +1.68804e-01 + 3.37139e-01*Nvtx;
155 py = -1.72555e-01 - 1.79594e-01*Nvtx;
156 } else {
157 //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
158 px = +2.22335e-02 - 6.59183e-02*Nvtx;
159 py = +1.52720e-01 - 1.28052e-01*Nvtx;
160 }
161 metx -= px;
162 mety -= py;
163 if (metx == 0.0 && mety == 0.0)
164 return 0.0;
165
166 double phi1 = std::atan2(mety,metx);
167 double phi2 = std::atan2(mety,metx)-2.0*M_PI;
168 if (std::abs(phi1-metphi) < std::abs(phi2-metphi)+0.5*M_PI)
169 return phi1;
170 else
171 return phi2;
172 }
173
174 TVector2 metType1Reg(double met, double metphi, double corr1, double corr2, double pt1, double eta1, double phi1, double e1, double pt2, double eta2, double phi2, double e2)
175 {
176 double metx = met * cos(metphi);
177 double mety = met * sin(metphi);
178 TLorentzVector j1;
179 TLorentzVector j2;
180 j1.SetPtEtaPhiE(pt1,eta1,phi1, e1 );
181 j2.SetPtEtaPhiE(pt2,eta2,phi2, e2 );
182 metx += j1.Px()*(1-corr1);
183 metx += j2.Px()*(1-corr2);
184 mety += j1.Py()*(1-corr1);
185 mety += j2.Py()*(1-corr2);
186 TVector2 corrMET(metx, mety);
187 return corrMET;
188 }
189
190 double metType1Phi(double met, double metphi, double corr1, double corr2, double pt1, double eta1, double phi1, double e1, double pt2, double eta2, double phi2, double e2){
191 return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Phi();
192
193 }
194 double metType1Et(double met, double metphi, double corr1, double corr2, double pt1, double eta1, double phi1, double e1, double pt2, double eta2, double phi2, double e2){
195 return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Mod();
196
197 }
198
199
200 double met_MPF(double met, double metphi, double pt, double phi)
201 {
202 return 1.+met*pt*std::cos( deltaPhi(metphi,phi) ) / (pt*pt);
203
204 }
205
206 double resolutionBias(double eta)
207 {
208 // return 0;//Nominal!
209 if(eta< 0.5) return 0.052;
210 if(eta< 1.1) return 0.057;
211 if(eta< 1.7) return 0.096;
212 if(eta< 2.3) return 0.134;
213 if(eta< 5) return 0.28;
214 return 0;
215 }
216
217 double evalJERBias( double ptreco, double ptgen, double eta1){
218 double eta = fabs(eta1);
219 double cor =1;
220 if ((fabs(ptreco - ptgen)/ ptreco)<0.5) { //Limit the effect to the core
221 cor = (ptreco +resolutionBias(eta) *(ptreco-ptgen))/ptreco;
222 }
223 if (ptgen > 0.) return ptreco*cor;
224 else return ptreco;
225 }
226
227 double evalEt( double pt, double eta, double phi, double e){
228 TLorentzVector j;
229 j.SetPtEtaPhiE(pt,eta,phi, e );
230 return j.Et();
231
232 }
233
234 double evalMt( double pt, double eta, double phi, double e){
235 TLorentzVector j;
236 j.SetPtEtaPhiE(pt,eta,phi, e );
237 return j.Mt();
238
239 }
240 /*double evalJECUnc( double pt, double eta){
241 // Total uncertainty for reference
242 JetCorrectionUncertainty *total = new JetCorrectionUncertainty("/shome/nmohr/CMSSW_5_2_6_patch1/src/UserCode/VHbb/data/START53_V15MC_Uncertainty_AK5PFchs.txt");
243
244 total->setJetPt(pt);
245 total->setJetEta(eta);
246 double uncert = total->getUncertainty(true);
247 delete total;
248 return uncert;
249 }*/
250
251 double ptWeightDY( double lheV_pt, double sign = 1.)
252 {
253 double SF = 1.;
254 if (50. < lheV_pt && lheV_pt < 100.){
255 SF = 0.873885+0.00175853*sign*lheV_pt;
256 }
257 else if (lheV_pt > 100){
258 SF = 1.10651-0.000705265*sign*lheV_pt;
259 }
260 return SF;
261 }
262
263 // weights correction for EWK NLO correction
264 double ptWeightZllH( double genHPt){
265 double SF = 1.;
266 if ( genHPt > 1.) SF = 0.999757-0.000174527*genHPt;
267 return SF>0?SF:0;
268 }
269
270
271 double minCSVold(int EVENT_run, double hJet_csv0, double hJet_csv1, double hJet_csvOld0, double hJet_csvOld1){
272
273 if( EVENT_run < 2 )
274 return std::min(hJet_csvOld0,hJet_csvOld1);
275 else
276 return std::min(hJet_csv0,hJet_csv1);
277 }
278
279 double maxCSVold(int EVENT_run, double hJet_csv0, double hJet_csv1, double hJet_csvOld0, double hJet_csvOld1){
280
281 if( EVENT_run < 2 )
282 return std::max(hJet_csvOld0,hJet_csvOld1);
283 else
284 return std::max(hJet_csv0,hJet_csv1);
285 }
286
287
288 // weights correction for EWK NLO correction from Atlas
289 double ewkAtlas8TeVZllH( double genHPt, double genZPt){
290 if (genHPt < 1.) return 1.;
291 double hll8_contents[95] = {0.000664024, -0.00357095, -0.00767076, -0.00967366, -0.0134844, -0.0157148, -0.0181885, -0.0209647, -0.0232788, -0.0252373, -0.0265634, -0.0275069, -0.0285776, -0.0281683, -0.0294206, -0.0299975, -0.0308047, -0.0311716, -0.030913, -0.0324821, -0.0323192, -0.0324639, -0.0319356, -0.0322621, -0.0331146, -0.0338905, -0.0345189, -0.0358591, -0.0358407, -0.040018, -0.0396389, -0.0407177, -0.0445103, -0.0441406, -0.0471215, -0.0463301, -0.0513777, -0.0536773, -0.0546446, -0.0568508, -0.0590333, -0.0612157, -0.0633981, -0.0655805, -0.067763, -0.0699454, -0.0721278, -0.0743103, -0.0764927, -0.0786751, -0.0808575, -0.08304, -0.0852224, -0.0874048, -0.0895872, -0.0917697, -0.0939521, -0.0961345, -0.098317, -0.100499, -0.102682, -0.104864, -0.107047, -0.109229, -0.111412, -0.113594, -0.115776, -0.117959, -0.120141, -0.122324, -0.124506, -0.126689, -0.128871, -0.131053, -0.133236, -0.135418, -0.137601, -0.139783, -0.141965, -0.144148, -0.14633, -0.148513, -0.150695, -0.152878, -0.15506, -0.157242, -0.159425, -0.161607, -0.16379, -0.165972, -0.168155, -0.170337, -0.172519, -0.174702, -0.176884};
292 int corrBin = (int) ( (genZPt - 25.) / 5. );
293 if (corrBin < 0) corrBin = 0;
294 if (corrBin > 94) corrBin = 94;
295 double SF = 1.+hll8_contents[corrBin];
296 return SF;
297 }
298
299 double eleEff(double pt, double eta){
300 if (pt >= 20.0 && pt < 25.0 && eta < -1.57000005245 && eta > -2.5) return 0.955429255962 ;
301 if (pt >= 20.0 && pt < 25.0 && eta < -1.44000005722 && eta > -1.57000005245) return 0.927380621433 ;
302 if (pt >= 20.0 && pt < 25.0 && eta < -0.800000011921 && eta > -1.44000005722) return 0.752811849117 ;
303 if (pt >= 20.0 && pt < 25.0 && eta < 0.0 && eta > -0.800000011921) return 0.758898377419 ;
304 if (pt >= 20.0 && pt < 25.0 && eta < 0.800000011921 && eta > 0.0) return 0.755395233631 ;
305 if (pt >= 20.0 && pt < 25.0 && eta < 1.44000005722 && eta > 0.800000011921) return 0.765102207661 ;
306 if (pt >= 20.0 && pt < 25.0 && eta < 1.57000005245 && eta > 1.44000005722) return 0.889689266682 ;
307 if (pt >= 20.0 && pt < 25.0 && eta < 2.5 && eta > 1.57000005245) return 0.960859179497 ;
308 if (pt >= 25.0 && pt < 30.0 && eta < -1.57000005245 && eta > -2.5) return 0.973749041557 ;
309 if (pt >= 25.0 && pt < 30.0 && eta < -1.44000005722 && eta > -1.57000005245) return 0.95562505722 ;
310 if (pt >= 25.0 && pt < 30.0 && eta < -0.800000011921 && eta > -1.44000005722) return 0.854979991913 ;
311 if (pt >= 25.0 && pt < 30.0 && eta < 0.0 && eta > -0.800000011921) return 0.859907507896 ;
312 if (pt >= 25.0 && pt < 30.0 && eta < 0.800000011921 && eta > 0.0) return 0.856850981712 ;
313 if (pt >= 25.0 && pt < 30.0 && eta < 1.44000005722 && eta > 0.800000011921) return 0.854210674763 ;
314 if (pt >= 25.0 && pt < 30.0 && eta < 1.57000005245 && eta > 1.44000005722) return 0.961790204048 ;
315 if (pt >= 25.0 && pt < 30.0 && eta < 2.5 && eta > 1.57000005245) return 0.984032571316 ;
316 if (pt >= 30.0 && pt < 35.0 && eta < -1.57000005245 && eta > -2.5) return 0.975343048573 ;
317 if (pt >= 30.0 && pt < 35.0 && eta < -1.44000005722 && eta > -1.57000005245) return 0.971944391727 ;
318 if (pt >= 30.0 && pt < 35.0 && eta < -0.800000011921 && eta > -1.44000005722) return 0.916224777699 ;
319 if (pt >= 30.0 && pt < 35.0 && eta < 0.0 && eta > -0.800000011921) return 0.918354153633 ;
320 if (pt >= 30.0 && pt < 35.0 && eta < 0.800000011921 && eta > 0.0) return 0.914323210716 ;
321 if (pt >= 30.0 && pt < 35.0 && eta < 1.44000005722 && eta > 0.800000011921) return 0.918937265873 ;
322 if (pt >= 30.0 && pt < 35.0 && eta < 1.57000005245 && eta > 1.44000005722) return 0.949005782604 ;
323 if (pt >= 30.0 && pt < 35.0 && eta < 2.5 && eta > 1.57000005245) return 0.978435099125 ;
324 if (pt >= 35.0 && pt < 40.0 && eta < -1.57000005245 && eta > -2.5) return 0.977416038513 ;
325 if (pt >= 35.0 && pt < 40.0 && eta < -1.44000005722 && eta > -1.57000005245) return 0.96676659584 ;
326 if (pt >= 35.0 && pt < 40.0 && eta < -0.800000011921 && eta > -1.44000005722) return 0.949674129486 ;
327 if (pt >= 35.0 && pt < 40.0 && eta < 0.0 && eta > -0.800000011921) return 0.955099225044 ;
328 if (pt >= 35.0 && pt < 40.0 && eta < 0.800000011921 && eta > 0.0) return 0.95435655117 ;
329 if (pt >= 35.0 && pt < 40.0 && eta < 1.44000005722 && eta > 0.800000011921) return 0.949129760265 ;
330 if (pt >= 35.0 && pt < 40.0 && eta < 1.57000005245 && eta > 1.44000005722) return 0.966977357864 ;
331 if (pt >= 35.0 && pt < 40.0 && eta < 2.5 && eta > 1.57000005245) return 0.974912643433 ;
332 if (pt >= 40.0 && pt < 45.0 && eta < -1.57000005245 && eta > -2.5) return 0.976379692554 ;
333 if (pt >= 40.0 && pt < 45.0 && eta < -1.44000005722 && eta > -1.57000005245) return 0.96400809288 ;
334 if (pt >= 40.0 && pt < 45.0 && eta < -0.800000011921 && eta > -1.44000005722) return 0.967367231846 ;
335 if (pt >= 40.0 && pt < 45.0 && eta < 0.0 && eta > -0.800000011921) return 0.976113498211 ;
336 if (pt >= 40.0 && pt < 45.0 && eta < 0.800000011921 && eta > 0.0) return 0.970276534557 ;
337 if (pt >= 40.0 && pt < 45.0 && eta < 1.44000005722 && eta > 0.800000011921) return 0.965076506138 ;
338 if (pt >= 40.0 && pt < 45.0 && eta < 1.57000005245 && eta > 1.44000005722) return 0.962332487106 ;
339 if (pt >= 40.0 && pt < 45.0 && eta < 2.5 && eta > 1.57000005245) return 0.975438177586 ;
340 if (pt >= 45.0 && pt < 50.0 && eta < -1.57000005245 && eta > -2.5) return 0.974642693996 ;
341 if (pt >= 45.0 && pt < 50.0 && eta < -1.44000005722 && eta > -1.57000005245) return 0.976564764977 ;
342 if (pt >= 45.0 && pt < 50.0 && eta < -0.800000011921 && eta > -1.44000005722) return 0.975869596004 ;
343 if (pt >= 45.0 && pt < 50.0 && eta < 0.0 && eta > -0.800000011921) return 0.984652400017 ;
344 if (pt >= 45.0 && pt < 50.0 && eta < 0.800000011921 && eta > 0.0) return 0.976347208023 ;
345 if (pt >= 45.0 && pt < 50.0 && eta < 1.44000005722 && eta > 0.800000011921) return 0.974825322628 ;
346 if (pt >= 45.0 && pt < 50.0 && eta < 1.57000005245 && eta > 1.44000005722) return 0.968582391739 ;
347 if (pt >= 45.0 && pt < 50.0 && eta < 2.5 && eta > 1.57000005245) return 0.973198652267 ;
348 if (pt >= 50.0 && pt < 200.0 && eta < -1.57000005245 && eta > -2.5) return 0.976521909237 ;
349 if (pt >= 50.0 && pt < 200.0 && eta < -1.44000005722 && eta > -1.57000005245) return 0.968289792538 ;
350 if (pt >= 50.0 && pt < 200.0 && eta < -0.800000011921 && eta > -1.44000005722) return 0.980251729488 ;
351 if (pt >= 50.0 && pt < 200.0 && eta < 0.0 && eta > -0.800000011921) return 0.988307952881 ;
352 if (pt >= 50.0 && pt < 200.0 && eta < 0.800000011921 && eta > 0.0) return 0.983105957508 ;
353 if (pt >= 50.0 && pt < 200.0 && eta < 1.44000005722 && eta > 0.800000011921) return 0.976873219013 ;
354 if (pt >= 50.0 && pt < 200.0 && eta < 1.57000005245 && eta > 1.44000005722) return 0.970499277115 ;
355 if (pt >= 50.0 && pt < 200.0 && eta < 2.5 && eta > 1.57000005245) return 0.972752988338 ;
356 return 1.;
357 }
358
359 double mueEff(int Vtype, double eta0, double eta1, double pt0, double pt1){
360 if (Vtype == 0) return 1.;
361 if (Vtype == 1) {
362 return eleEff(pt0,eta0)*eleEff(pt1,eta1);
363 }
364 return 1.;
365 }
366
367 }
368