197 double mu_p = muonP4.P();
198 double mu_pt = muonP4.Pt();
199 double mu_phi = muonP4.Phi();
200 double mu_eta = muonP4.Eta();
201 double mu_vz = muonVertex.z()/100.;
202 double mu_pz = muonP4.Pz();
204 double ecalPhi, ecalTheta;
205 double hcalPhi, hcalTheta;
206 double hoPhi, hoTheta;
213 ecalPhi = muonMETInfo.
ecalPos.Phi();
214 ecalTheta = muonMETInfo.
ecalPos.Theta();
215 hcalPhi = muonMETInfo.
hcalPos.Phi();
216 hcalTheta = muonMETInfo.
hcalPos.Theta();
217 hoPhi = muonMETInfo.
hoPos.Phi();
218 hoTheta = muonMETInfo.
hoPos.Theta();
229 double rEcal = 1.290;
232 if(
abs(mu_eta) > 0.3) rHo = 4.07;
234 double zFaceEcal = 3.209;
235 if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
237 double zFaceHcal = 3.88;
238 if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;
242 double bendr = mu_pt*1000/(300*bfield);
244 double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr));
245 double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr));
246 double tb_ho = TMath::ACos(1-rHo*rHo/(2*bendr*bendr));
247 double xEcal,yEcal,zEcal;
248 double xHcal,yHcal,zHcal;
252 if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) {
253 xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
254 yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
255 zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
258 double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
259 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
260 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
261 zEcal = fabs(zFaceEcal);
263 double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
264 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
265 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
266 zEcal = -fabs(zFaceEcal);
271 if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) {
272 xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
273 yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
274 zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
277 double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
278 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
279 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
280 zHcal = fabs(zFaceHcal);
282 double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
283 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
284 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
285 zHcal = -fabs(zFaceHcal);
290 xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
291 yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
292 zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;
294 ecalTheta = TMath::ACos(zEcal/
sqrt(
pow(xEcal,2) +
pow(yEcal,2)+
pow(zEcal,2)));
295 ecalPhi = atan2(yEcal,xEcal);
296 hcalTheta = TMath::ACos(zHcal/
sqrt(
pow(xHcal,2) +
pow(yHcal,2)+
pow(zHcal,2)));
297 hcalPhi = atan2(yHcal,xHcal);
298 hoTheta = TMath::ACos(zHo/
sqrt(
pow(xHo,2) +
pow(yHo,2)+
pow(zHo,2)));
299 hoPhi = atan2(yHo,xHo);
302 double r2dEcal =
sqrt(
pow(xEcal,2)+
pow(yEcal,2));
303 double r2dHcal =
sqrt(
pow(xHcal,2)+
pow(yHcal,2));
314 double dphi = mu_phi - ecalPhi;
317 ecalPhi = mu_phi - fabs(dphi);
320 ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
322 xEcal = r2dEcal*TMath::Cos(ecalPhi);
323 yEcal = r2dEcal*TMath::Sin(ecalPhi);
326 dphi = mu_phi - hcalPhi;
329 hcalPhi = mu_phi - fabs(dphi);
332 hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
334 xHcal = r2dHcal*TMath::Cos(hcalPhi);
335 yHcal = r2dHcal*TMath::Sin(hcalPhi);
339 dphi = mu_phi - hoPhi;
342 hoPhi = mu_phi - fabs(dphi);
345 hoPhi = -1*temp*hoPhi/fabs(hoPhi);
347 xHo = r2dHo*TMath::Cos(hoPhi);
348 yHo = r2dHo*TMath::Sin(hoPhi);
356 double mu_Ex = muonMETInfo.
ecalE*
sin(ecalTheta)*
cos(ecalPhi)
358 + muonMETInfo.
hoE*
sin(hoTheta)*
cos(hoPhi);
359 double mu_Ey = muonMETInfo.
ecalE*
sin(ecalTheta)*
sin(ecalPhi)
361 + muonMETInfo.
hoE*
sin(hoTheta)*
sin(hoPhi);
373 double dEdx_normalization = -(11.4 + 0.96*fabs(
log(50*2.8)) + 0.033*50*(1.0 -
pow(50, -0.33)) )*1
e-3;
374 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;
378 if(muonMETInfo.
useHO) {
380 if(fabs(mu_eta) < 0.2)
381 temp = 2.75*(1-0.00003*mu_p);
382 if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
383 temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
384 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
385 temp = 7.413-5.12*fabs(mu_eta);
386 if(fabs(mu_eta) > 1.3)
387 temp = 2.084-0.743*fabs(mu_eta);
389 if(fabs(mu_eta) < 1.0)
390 temp = 2.33*(1-0.0004*mu_p);
391 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
392 temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
393 if(fabs(mu_eta) > 1.3)
394 temp = 2.084-0.743*fabs(mu_eta);
397 double dep = temp*dEdx_normalization/dEdx_numerator;
398 if(dep < 0.5) dep = 0;
400 if(fabs(mu_eta) < 1.3) {
401 deltax += dep*
cos((ecalPhi+hcalPhi+hoPhi)/3);
402 deltay += dep*
sin((ecalPhi+hcalPhi+hoPhi)/3);
404 deltax += dep*
cos( (ecalPhi+hcalPhi)/2);
405 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)