198 double mu_p = muonP4.P();
199 double mu_pt = muonP4.Pt();
200 double mu_phi = muonP4.Phi();
201 double mu_eta = muonP4.Eta();
202 double mu_vz = muonVertex.z()/100.;
203 double mu_pz = muonP4.Pz();
205 double ecalPhi, ecalTheta;
206 double hcalPhi, hcalTheta;
207 double hoPhi, hoTheta;
214 ecalPhi = muonMETInfo.
ecalPos.Phi();
215 ecalTheta = muonMETInfo.
ecalPos.Theta();
216 hcalPhi = muonMETInfo.
hcalPos.Phi();
217 hcalTheta = muonMETInfo.
hcalPos.Theta();
218 hoPhi = muonMETInfo.
hoPos.Phi();
219 hoTheta = muonMETInfo.
hoPos.Theta();
230 double rEcal = 1.290;
233 if(
abs(mu_eta) > 0.3) rHo = 4.07;
235 double zFaceEcal = 3.209;
236 if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
238 double zFaceHcal = 3.88;
239 if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;
243 double bendr = mu_pt*1000/(300*bfield);
245 double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr));
246 double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr));
247 double tb_ho = TMath::ACos(1-rHo*rHo/(2*bendr*bendr));
248 double xEcal,yEcal,zEcal;
249 double xHcal,yHcal,zHcal;
253 if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) {
254 xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
255 yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
256 zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
259 double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
260 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
261 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
262 zEcal = fabs(zFaceEcal);
264 double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
265 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
266 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
267 zEcal = -fabs(zFaceEcal);
272 if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) {
273 xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
274 yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
275 zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
278 double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
279 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
280 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
281 zHcal = fabs(zFaceHcal);
283 double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
284 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
285 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
286 zHcal = -fabs(zFaceHcal);
291 xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
292 yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
293 zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;
295 ecalTheta = TMath::ACos(zEcal/
sqrt(
pow(xEcal,2) +
pow(yEcal,2)+
pow(zEcal,2)));
296 ecalPhi = atan2(yEcal,xEcal);
297 hcalTheta = TMath::ACos(zHcal/
sqrt(
pow(xHcal,2) +
pow(yHcal,2)+
pow(zHcal,2)));
298 hcalPhi = atan2(yHcal,xHcal);
299 hoTheta = TMath::ACos(zHo/
sqrt(
pow(xHo,2) +
pow(yHo,2)+
pow(zHo,2)));
300 hoPhi = atan2(yHo,xHo);
303 double r2dEcal =
sqrt(
pow(xEcal,2)+
pow(yEcal,2));
304 double r2dHcal =
sqrt(
pow(xHcal,2)+
pow(yHcal,2));
315 double dphi = mu_phi - ecalPhi;
318 ecalPhi = mu_phi - fabs(dphi);
321 ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
323 xEcal = r2dEcal*TMath::Cos(ecalPhi);
324 yEcal = r2dEcal*TMath::Sin(ecalPhi);
327 dphi = mu_phi - hcalPhi;
330 hcalPhi = mu_phi - fabs(dphi);
333 hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
335 xHcal = r2dHcal*TMath::Cos(hcalPhi);
336 yHcal = r2dHcal*TMath::Sin(hcalPhi);
340 dphi = mu_phi - hoPhi;
343 hoPhi = mu_phi - fabs(dphi);
346 hoPhi = -1*temp*hoPhi/fabs(hoPhi);
348 xHo = r2dHo*TMath::Cos(hoPhi);
349 yHo = r2dHo*TMath::Sin(hoPhi);
357 double mu_Ex = muonMETInfo.
ecalE*
sin(ecalTheta)*
cos(ecalPhi)
359 + muonMETInfo.
hoE*
sin(hoTheta)*
cos(hoPhi);
360 double mu_Ey = muonMETInfo.
ecalE*
sin(ecalTheta)*
sin(ecalPhi)
362 + muonMETInfo.
hoE*
sin(hoTheta)*
sin(hoPhi);
374 double dEdx_normalization = -(11.4 + 0.96*fabs(
log(50*2.8)) + 0.033*50*(1.0 -
pow(50, -0.33)) )*1
e-3;
375 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;
379 if(muonMETInfo.
useHO) {
381 if(fabs(mu_eta) < 0.2)
382 temp = 2.75*(1-0.00003*mu_p);
383 if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
384 temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
385 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
386 temp = 7.413-5.12*fabs(mu_eta);
387 if(fabs(mu_eta) > 1.3)
388 temp = 2.084-0.743*fabs(mu_eta);
390 if(fabs(mu_eta) < 1.0)
391 temp = 2.33*(1-0.0004*mu_p);
392 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
393 temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
394 if(fabs(mu_eta) > 1.3)
395 temp = 2.084-0.743*fabs(mu_eta);
398 double dep = temp*dEdx_normalization/dEdx_numerator;
399 if(dep < 0.5) dep = 0;
401 if(fabs(mu_eta) < 1.3) {
402 deltax += dep*
cos((ecalPhi+hcalPhi+hoPhi)/3);
403 deltay += dep*
sin((ecalPhi+hcalPhi+hoPhi)/3);
405 deltax += dep*
cos( (ecalPhi+hcalPhi)/2);
406 deltay += dep*
cos( (ecalPhi+hcalPhi)/2);
Sin< T >::type sin(const T &t)
Cos< T >::type cos(const T &t)
bool useTkAssociatorPositions
Log< T >::type log(const T &t)
Power< A, B >::type pow(const A &a, const B &b)