200 double mu_p = muonP4.P();
201 double mu_pt = muonP4.Pt();
202 double mu_phi = muonP4.Phi();
203 double mu_eta = muonP4.Eta();
204 double mu_vz = muonVertex.z()/100.;
205 double mu_pz = muonP4.Pz();
207 double ecalPhi, ecalTheta;
208 double hcalPhi, hcalTheta;
209 double hoPhi, hoTheta;
216 ecalPhi = muonMETInfo.
ecalPos.Phi();
217 ecalTheta = muonMETInfo.
ecalPos.Theta();
218 hcalPhi = muonMETInfo.
hcalPos.Phi();
219 hcalTheta = muonMETInfo.
hcalPos.Theta();
220 hoPhi = muonMETInfo.
hoPos.Phi();
221 hoTheta = muonMETInfo.
hoPos.Theta();
232 double rEcal = 1.290;
235 if(
abs(mu_eta) > 0.3) rHo = 4.07;
237 double zFaceEcal = 3.209;
238 if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
240 double zFaceHcal = 3.88;
241 if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;
245 double bendr = mu_pt*1000/(300*bfield);
247 double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr));
248 double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr));
249 double tb_ho = TMath::ACos(1-rHo*rHo/(2*bendr*bendr));
250 double xEcal,yEcal,zEcal;
251 double xHcal,yHcal,zHcal;
255 if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) {
256 xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
257 yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
258 zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
261 double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
262 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
263 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
264 zEcal = fabs(zFaceEcal);
266 double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
267 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
268 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
269 zEcal = -fabs(zFaceEcal);
274 if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) {
275 xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
276 yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
277 zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
280 double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
281 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
282 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
283 zHcal = fabs(zFaceHcal);
285 double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
286 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
287 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
288 zHcal = -fabs(zFaceHcal);
293 xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
294 yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
295 zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;
297 ecalTheta = TMath::ACos(zEcal/
sqrt(
pow(xEcal,2) +
pow(yEcal,2)+
pow(zEcal,2)));
298 ecalPhi = atan2(yEcal,xEcal);
299 hcalTheta = TMath::ACos(zHcal/
sqrt(
pow(xHcal,2) +
pow(yHcal,2)+
pow(zHcal,2)));
300 hcalPhi = atan2(yHcal,xHcal);
301 hoTheta = TMath::ACos(zHo/
sqrt(
pow(xHo,2) +
pow(yHo,2)+
pow(zHo,2)));
302 hoPhi = atan2(yHo,xHo);
305 double r2dEcal =
sqrt(
pow(xEcal,2)+
pow(yEcal,2));
306 double r2dHcal =
sqrt(
pow(xHcal,2)+
pow(yHcal,2));
317 double dphi = mu_phi - ecalPhi;
320 ecalPhi = mu_phi - fabs(dphi);
323 ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
325 xEcal = r2dEcal*TMath::Cos(ecalPhi);
326 yEcal = r2dEcal*TMath::Sin(ecalPhi);
329 dphi = mu_phi - hcalPhi;
332 hcalPhi = mu_phi - fabs(dphi);
335 hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
337 xHcal = r2dHcal*TMath::Cos(hcalPhi);
338 yHcal = r2dHcal*TMath::Sin(hcalPhi);
342 dphi = mu_phi - hoPhi;
345 hoPhi = mu_phi - fabs(dphi);
348 hoPhi = -1*temp*hoPhi/fabs(hoPhi);
350 xHo = r2dHo*TMath::Cos(hoPhi);
351 yHo = r2dHo*TMath::Sin(hoPhi);
359 double mu_Ex = muonMETInfo.
ecalE*
sin(ecalTheta)*
cos(ecalPhi)
361 + muonMETInfo.
hoE*
sin(hoTheta)*
cos(hoPhi);
362 double mu_Ey = muonMETInfo.
ecalE*
sin(ecalTheta)*
sin(ecalPhi)
364 + muonMETInfo.
hoE*
sin(hoTheta)*
sin(hoPhi);
376 double dEdx_normalization = -(11.4 + 0.96*fabs(
log(50*2.8)) + 0.033*50*(1.0 -
pow(50, -0.33)) )*1
e-3;
377 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;
381 if(muonMETInfo.
useHO) {
383 if(fabs(mu_eta) < 0.2)
384 temp = 2.75*(1-0.00003*mu_p);
385 if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
386 temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
387 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
388 temp = 7.413-5.12*fabs(mu_eta);
389 if(fabs(mu_eta) > 1.3)
390 temp = 2.084-0.743*fabs(mu_eta);
392 if(fabs(mu_eta) < 1.0)
393 temp = 2.33*(1-0.0004*mu_p);
394 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
395 temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
396 if(fabs(mu_eta) > 1.3)
397 temp = 2.084-0.743*fabs(mu_eta);
400 double dep = temp*dEdx_normalization/dEdx_numerator;
401 if(dep < 0.5) dep = 0;
403 if(fabs(mu_eta) < 1.3) {
404 deltax += dep*
cos((ecalPhi+hcalPhi+hoPhi)/3);
405 deltay += dep*
sin((ecalPhi+hcalPhi+hoPhi)/3);
407 deltax += dep*
cos( (ecalPhi+hcalPhi)/2);
408 deltay += dep*
cos( (ecalPhi+hcalPhi)/2);
Sin< T >::type sin(const T &t)
Cos< T >::type cos(const T &t)
Abs< T >::type abs(const T &t)
bool useTkAssociatorPositions
Power< A, B >::type pow(const A &a, const B &b)