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 |
|
{ |
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; |
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(); |
87 |
> |
TVector3 bZ = msum.BoostVector(); |
88 |
|
|
89 |
< |
m1.Boost(-bZ); |
90 |
< |
m2.Boost(-bZ); |
89 |
> |
m1.Boost(-bZ); |
90 |
> |
m2.Boost(-bZ); |
91 |
|
|
92 |
< |
TVector3 b1; |
92 |
> |
TVector3 b1; |
93 |
|
|
94 |
|
|
95 |
< |
if((int) (pt) % 2 == 0) |
96 |
< |
b1 = m1.BoostVector(); |
97 |
< |
else |
98 |
< |
b1 = m2.BoostVector(); |
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); |
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; |
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(); |
111 |
> |
TVector3 bZ = msum.BoostVector(); |
112 |
|
|
113 |
< |
m1.Boost(-bZ); |
114 |
< |
m2.Boost(-bZ); |
113 |
> |
m1.Boost(-bZ); |
114 |
> |
m2.Boost(-bZ); |
115 |
|
|
116 |
< |
TVector3 b1; |
116 |
> |
TVector3 b1; |
117 |
|
|
118 |
< |
if((int) (pt) % 2 == 0) |
119 |
< |
b1 = m1.BoostVector(); |
120 |
< |
else |
121 |
< |
b1 = m2.BoostVector(); |
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 |
< |
} |
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 |
< |
{ |
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; |
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; |
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 |
< |
} |
144 |
> |
} |
145 |
|
|
146 |
< |
double metphiCorSysShift(double met, double metphi, int Nvtx, int EVENT_run) |
147 |
< |
{ |
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; |
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; |
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; |
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; |
169 |
> |
return phi1; |
170 |
|
else |
171 |
< |
return phi2; |
172 |
< |
} |
171 |
> |
return phi2; |
172 |
> |
} |
173 |
|
|
174 |
< |
double resolutionBias(double eta) |
175 |
< |
{ |
176 |
< |
// return 0;//Nominal! |
177 |
< |
if(eta< 1.1) return 0.05; |
178 |
< |
if(eta< 2.5) return 0.10; |
179 |
< |
if(eta< 5) return 0.30; |
180 |
< |
return 0; |
181 |
< |
} |
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 |
|
|
177 |
– |
double evalJERBias( double ptreco, double ptgen, double eta1){ |
178 |
– |
double eta = fabs(eta1); |
179 |
– |
double cor =1; |
180 |
– |
if ((fabs(ptreco - ptgen)/ ptreco)<0.5) { //Limit the effect to the core |
181 |
– |
cor = (ptreco +resolutionBias(eta) *(ptreco-ptgen))/ptreco; |
193 |
|
} |
194 |
< |
return ptreco*cor; |
195 |
< |
} |
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 |
< |
double evalEt( double pt, double eta, double phi, double e){ |
187 |
< |
TLorentzVector j; |
188 |
< |
j.SetPtEtaPhiE(pt,eta,phi, e ); |
189 |
< |
return j.Et(); |
197 |
> |
} |
198 |
|
|
191 |
– |
} |
199 |
|
|
200 |
< |
double evalMt( double pt, double eta, double phi, double e){ |
201 |
< |
TLorentzVector j; |
202 |
< |
j.SetPtEtaPhiE(pt,eta,phi, e ); |
196 |
< |
return j.Mt(); |
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 |
< |
} |
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 min(hJet_csvOld0,hJet_csvOld1); |
275 |
> |
else |
276 |
> |
return 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 max(hJet_csvOld0,hJet_csvOld1); |
283 |
> |
else |
284 |
> |
return 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 |
|
} |
300 |
|
|