195 double mu_p = muonP4.P();
196 double mu_pt = muonP4.Pt();
197 double mu_phi = muonP4.Phi();
198 double mu_eta = muonP4.Eta();
199 double mu_vz = muonVertex.z() / 100.;
200 double mu_pz = muonP4.Pz();
203 double hcalPhi, hcalTheta;
204 double hoPhi, hoTheta;
211 ecalTheta = muonMETInfo.
ecalPos.Theta();
212 hcalPhi = muonMETInfo.
hcalPos.Phi();
213 hcalTheta = muonMETInfo.
hcalPos.Theta();
214 hoPhi = muonMETInfo.
hoPos.Phi();
215 hoTheta = muonMETInfo.
hoPos.Theta();
225 double rEcal = 1.290;
228 if (
abs(mu_eta) > 0.3)
231 double zFaceEcal = 3.209;
233 zFaceEcal = -1 * zFaceEcal;
235 double zFaceHcal = 3.88;
237 zFaceHcal = -1 * zFaceHcal;
241 double bendr = mu_pt * 1000 / (300 * bfield);
243 double tb_ecal = TMath::ACos(1 - rEcal * rEcal / (2 * bendr * bendr));
244 double tb_hcal = TMath::ACos(1 - rHcal * rHcal / (2 * bendr * bendr));
245 double tb_ho = TMath::ACos(1 - rHo * rHo / (2 * bendr * bendr));
246 double xEcal, yEcal, zEcal;
247 double xHcal, yHcal, zHcal;
248 double xHo, yHo, zHo;
251 if (fabs(mu_pz * bendr * tb_ecal / mu_pt + mu_vz) < fabs(zFaceEcal) && rEcal < 2 * bendr) {
252 xEcal = bendr * (TMath::Sin(tb_ecal + mu_phi) - TMath::Sin(mu_phi));
253 yEcal = bendr * (-TMath::Cos(tb_ecal + mu_phi) + TMath::Cos(mu_phi));
254 zEcal = bendr * tb_ecal * mu_pz / mu_pt + mu_vz;
257 double te_ecal = (fabs(zFaceEcal) - mu_vz) * mu_pt / (bendr * mu_pz);
258 xEcal = bendr * (TMath::Sin(te_ecal + mu_phi) - TMath::Sin(mu_phi));
259 yEcal = bendr * (-TMath::Cos(te_ecal + mu_phi) + TMath::Cos(mu_phi));
260 zEcal = fabs(zFaceEcal);
262 double te_ecal = -(fabs(zFaceEcal) + mu_vz) * mu_pt / (bendr * mu_pz);
263 xEcal = bendr * (TMath::Sin(te_ecal + mu_phi) - TMath::Sin(mu_phi));
264 yEcal = bendr * (-TMath::Cos(te_ecal + mu_phi) + TMath::Cos(mu_phi));
265 zEcal = -fabs(zFaceEcal);
270 if (fabs(mu_pz * bendr * tb_hcal / mu_pt + mu_vz) < fabs(zFaceHcal) && rEcal < 2 * bendr) {
271 xHcal = bendr * (TMath::Sin(tb_hcal + mu_phi) - TMath::Sin(mu_phi));
272 yHcal = bendr * (-TMath::Cos(tb_hcal + mu_phi) + TMath::Cos(mu_phi));
273 zHcal = bendr * tb_hcal * mu_pz / mu_pt + mu_vz;
276 double te_hcal = (fabs(zFaceHcal) - mu_vz) * mu_pt / (bendr * mu_pz);
277 xHcal = bendr * (TMath::Sin(te_hcal + mu_phi) - TMath::Sin(mu_phi));
278 yHcal = bendr * (-TMath::Cos(te_hcal + mu_phi) + TMath::Cos(mu_phi));
279 zHcal = fabs(zFaceHcal);
281 double te_hcal = -(fabs(zFaceHcal) + mu_vz) * mu_pt / (bendr * mu_pz);
282 xHcal = bendr * (TMath::Sin(te_hcal + mu_phi) - TMath::Sin(mu_phi));
283 yHcal = bendr * (-TMath::Cos(te_hcal + mu_phi) + TMath::Cos(mu_phi));
284 zHcal = -fabs(zFaceHcal);
289 xHo = bendr * (TMath::Sin(tb_ho + mu_phi) - TMath::Sin(mu_phi));
290 yHo = bendr * (-TMath::Cos(tb_ho + mu_phi) + TMath::Cos(mu_phi));
291 zHo = bendr * tb_ho * mu_pz / mu_pt + mu_vz;
293 ecalTheta = TMath::ACos(zEcal /
sqrt(
pow(xEcal, 2) +
pow(yEcal, 2) +
pow(zEcal, 2)));
295 hcalTheta = TMath::ACos(zHcal /
sqrt(
pow(xHcal, 2) +
pow(yHcal, 2) +
pow(zHcal, 2)));
296 hcalPhi = atan2(yHcal, xHcal);
297 hoTheta = TMath::ACos(zHo /
sqrt(
pow(xHo, 2) +
pow(yHo, 2) +
pow(zHo, 2)));
298 hoPhi = atan2(yHo, xHo);
305 if (muonCharge > 0) {
307 double dphi = mu_phi -
ecalPhi;
317 dphi = mu_phi - hcalPhi;
320 hcalPhi = mu_phi - fabs(dphi);
323 hcalPhi = -1 *
temp * hcalPhi / fabs(hcalPhi);
327 dphi = mu_phi - hoPhi;
330 hoPhi = mu_phi - fabs(dphi);
333 hoPhi = -1 *
temp * hoPhi / fabs(hoPhi);
341 muonMETInfo.
hcalE *
sin(hcalTheta) *
cos(hcalPhi) + muonMETInfo.
hoE *
sin(hoTheta) *
cos(hoPhi);
343 muonMETInfo.
hcalE *
sin(hcalTheta) *
sin(hcalPhi) + muonMETInfo.
hoE *
sin(hoTheta) *
sin(hoPhi);
355 double dEdx_normalization = -(11.4 + 0.96 * fabs(
log(50 * 2.8)) + 0.033 * 50 * (1.0 -
pow(50, -0.33))) * 1
e-3;
356 double dEdx_numerator = -(11.4 + 0.96 * fabs(
log(mu_p * 2.8)) + 0.033 * mu_p * (1.0 -
pow(mu_p, -0.33))) * 1
e-3;
360 if (muonMETInfo.
useHO) {
362 if (fabs(mu_eta) < 0.2)
363 temp = 2.75 * (1 - 0.00003 * mu_p);
364 if (fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
365 temp = (2.38 + 0.0144 * fabs(mu_eta)) * (1 - 0.0003 * mu_p);
366 if (fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
367 temp = 7.413 - 5.12 * fabs(mu_eta);
368 if (fabs(mu_eta) > 1.3)
369 temp = 2.084 - 0.743 * fabs(mu_eta);
371 if (fabs(mu_eta) < 1.0)
372 temp = 2.33 * (1 - 0.0004 * mu_p);
373 if (fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
374 temp = (7.413 - 5.12 * fabs(mu_eta)) * (1 - 0.0003 * mu_p);
375 if (fabs(mu_eta) > 1.3)
376 temp = 2.084 - 0.743 * fabs(mu_eta);
379 double dep =
temp * dEdx_normalization / dEdx_numerator;
383 if (fabs(mu_eta) < 1.3) {
384 deltax += dep *
cos((
ecalPhi + hcalPhi + hoPhi) / 3);
385 deltay += dep *
sin((
ecalPhi + hcalPhi + hoPhi) / 3);