File indexing completed on 2024-04-06 12:26:41
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include <cmath>
0012 #include <vector>
0013 #include "RecoMET/METAlgorithms/interface/MuonMETAlgo.h"
0014 #include "DataFormats/METReco/interface/CaloMET.h"
0015 #include "DataFormats/METReco/interface/SpecificCaloMETData.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0018
0019 #include "TMath.h"
0020
0021 using namespace std;
0022 using namespace reco;
0023
0024 typedef math::XYZTLorentzVector LorentzVector;
0025 typedef math::XYZPoint Point;
0026
0027
0028 CaloMET MuonMETAlgo::makeMET(const CaloMET& fMet,
0029 double fSumEt,
0030 const std::vector<CorrMETData>& fCorrections,
0031 const MET::LorentzVector& fP4) {
0032 return CaloMET(fMet.getSpecific(), fSumEt, fCorrections, fP4, fMet.vertex());
0033 }
0034
0035
0036 MET MuonMETAlgo::makeMET(const MET& fMet,
0037 double fSumEt,
0038 const std::vector<CorrMETData>& fCorrections,
0039 const MET::LorentzVector& fP4) {
0040 return MET(fSumEt, fCorrections, fP4, fMet.vertex());
0041 }
0042
0043
0044 template <class T>
0045 void MuonMETAlgo::MuonMETAlgo_run(const edm::View<reco::Muon>& inputMuons,
0046 const edm::ValueMap<reco::MuonMETCorrectionData>& vm_muCorrData,
0047 const edm::View<T>& v_uncorMET,
0048 std::vector<T>* v_corMET) {
0049 const T& uncorMETObj = v_uncorMET.front();
0050
0051 double corMETX = uncorMETObj.px();
0052 double corMETY = uncorMETObj.py();
0053
0054 double sumMuPx = 0.;
0055 double sumMuPy = 0.;
0056 double sumMuPt = 0.;
0057 double sumMuDepEx = 0.;
0058 double sumMuDepEy = 0.;
0059 double sumMuDepEt = 0.;
0060
0061 for (unsigned int i = 0; i < inputMuons.size(); ++i) {
0062 reco::MuonMETCorrectionData muCorrData = vm_muCorrData[inputMuons.refAt(i)];
0063
0064 if (muCorrData.type() == 0)
0065 continue;
0066
0067 float deltax = muCorrData.corrX();
0068 float deltay = muCorrData.corrY();
0069
0070 const reco::Muon* mu = &inputMuons[i];
0071 const LorentzVector& mup4 = mu->p4();
0072
0073 sumMuPx += mup4.px();
0074 sumMuPy += mup4.py();
0075 sumMuPt += mup4.pt();
0076 sumMuDepEx += deltax;
0077 sumMuDepEy += deltay;
0078 sumMuDepEt += sqrt(deltax * deltax + deltay * deltay);
0079 corMETX = corMETX - mup4.px() + deltax;
0080 corMETY = corMETY - mup4.py() + deltay;
0081 }
0082
0083 CorrMETData delta;
0084 delta.mex = sumMuDepEx - sumMuPx;
0085 delta.mey = sumMuDepEy - sumMuPy;
0086 delta.sumet = sumMuPt - sumMuDepEt;
0087 MET::LorentzVector correctedMET4vector(corMETX, corMETY, 0., sqrt(corMETX * corMETX + corMETY * corMETY));
0088 std::vector<CorrMETData> corrections = uncorMETObj.mEtCorr();
0089 corrections.push_back(delta);
0090
0091 T result = makeMET(uncorMETObj, uncorMETObj.sumEt() + delta.sumet, corrections, correctedMET4vector);
0092 v_corMET->push_back(result);
0093 }
0094
0095
0096 void MuonMETAlgo::GetMuDepDeltas(const reco::Muon* inputMuon,
0097 TrackDetMatchInfo& info,
0098 bool useTrackAssociatorPositions,
0099 bool useRecHits,
0100 bool useHO,
0101 double towerEtThreshold,
0102 double& deltax,
0103 double& deltay,
0104 double Bfield) {
0105 bool useAverage = false;
0106
0107
0108 double sumPt = inputMuon->isIsolationValid() ? inputMuon->isolationR03().sumPt : 0.0;
0109 double sumEtEcal = inputMuon->isIsolationValid() ? inputMuon->isolationR03().emEt : 0.0;
0110 double sumEtHcal = inputMuon->isIsolationValid() ? inputMuon->isolationR03().hadEt : 0.0;
0111
0112 if (sumPt > 3 || sumEtEcal + sumEtHcal > 5)
0113 useAverage = true;
0114
0115
0116
0117 MuonMETInfo muMETInfo;
0118 muMETInfo.useAverage = useAverage;
0119 muMETInfo.useTkAssociatorPositions = useTrackAssociatorPositions;
0120 muMETInfo.useHO = useHO;
0121
0122 TrackRef mu_track;
0123 if (inputMuon->isGlobalMuon()) {
0124 mu_track = inputMuon->globalTrack();
0125 } else if (inputMuon->isTrackerMuon() || inputMuon->isRPCMuon() || inputMuon->isGEMMuon() || inputMuon->isME0Muon()) {
0126 mu_track = inputMuon->innerTrack();
0127 } else
0128 mu_track = inputMuon->outerTrack();
0129
0130 if (useTrackAssociatorPositions) {
0131 muMETInfo.ecalPos = info.trkGlobPosAtEcal;
0132 muMETInfo.hcalPos = info.trkGlobPosAtHcal;
0133 muMETInfo.hoPos = info.trkGlobPosAtHO;
0134 }
0135
0136 if (!useAverage) {
0137 if (useRecHits) {
0138 muMETInfo.ecalE = inputMuon->calEnergy().emS9;
0139 muMETInfo.hcalE = inputMuon->calEnergy().hadS9;
0140 if (useHO)
0141 muMETInfo.hoE = inputMuon->calEnergy().hoS9;
0142 } else {
0143
0144
0145 std::vector<const CaloTower*> towers = info.crossedTowers;
0146 for (vector<const CaloTower*>::const_iterator it = towers.begin(); it != towers.end(); it++) {
0147 if ((*it)->et() < towerEtThreshold)
0148 continue;
0149 muMETInfo.ecalE += (*it)->emEnergy();
0150 muMETInfo.hcalE += (*it)->hadEnergy();
0151 if (useHO)
0152 muMETInfo.hoE += (*it)->outerEnergy();
0153 }
0154 }
0155 }
0156
0157
0158
0159 math::XYZTLorentzVector mup4;
0160 if (inputMuon->isGlobalMuon()) {
0161 if (inputMuon->globalTrack()->pt() < 200) {
0162 mup4 = LorentzVector(inputMuon->innerTrack()->px(),
0163 inputMuon->innerTrack()->py(),
0164 inputMuon->innerTrack()->pz(),
0165 inputMuon->innerTrack()->p());
0166 } else {
0167 mup4 = LorentzVector(inputMuon->globalTrack()->px(),
0168 inputMuon->globalTrack()->py(),
0169 inputMuon->globalTrack()->pz(),
0170 inputMuon->globalTrack()->p());
0171 }
0172 } else if (inputMuon->isTrackerMuon() || inputMuon->isRPCMuon() || inputMuon->isGEMMuon() || inputMuon->isME0Muon()) {
0173 mup4 = LorentzVector(inputMuon->innerTrack()->px(),
0174 inputMuon->innerTrack()->py(),
0175 inputMuon->innerTrack()->pz(),
0176 inputMuon->innerTrack()->p());
0177 } else
0178 mup4 = LorentzVector(inputMuon->outerTrack()->px(),
0179 inputMuon->outerTrack()->py(),
0180 inputMuon->outerTrack()->pz(),
0181 inputMuon->outerTrack()->p());
0182
0183
0184 correctMETforMuon(deltax, deltay, Bfield, inputMuon->charge(), mup4, inputMuon->vertex(), muMETInfo);
0185 }
0186
0187
0188 void MuonMETAlgo::correctMETforMuon(double& deltax,
0189 double& deltay,
0190 double bfield,
0191 int muonCharge,
0192 const math::XYZTLorentzVector& muonP4,
0193 const math::XYZPoint& muonVertex,
0194 MuonMETInfo& muonMETInfo) {
0195 double mu_p = muonP4.P();
0196 double mu_pt = muonP4.Pt();
0197 double mu_phi = muonP4.Phi();
0198 double mu_eta = muonP4.Eta();
0199 double mu_vz = muonVertex.z() / 100.;
0200 double mu_pz = muonP4.Pz();
0201
0202 double ecalPhi, ecalTheta;
0203 double hcalPhi, hcalTheta;
0204 double hoPhi, hoTheta;
0205
0206
0207
0208
0209 if (muonMETInfo.useTkAssociatorPositions) {
0210 ecalPhi = muonMETInfo.ecalPos.Phi();
0211 ecalTheta = muonMETInfo.ecalPos.Theta();
0212 hcalPhi = muonMETInfo.hcalPos.Phi();
0213 hcalTheta = muonMETInfo.hcalPos.Theta();
0214 hoPhi = muonMETInfo.hoPos.Phi();
0215 hoTheta = muonMETInfo.hoPos.Theta();
0216 } else {
0217
0218
0219
0220
0221
0222
0223
0224
0225 double rEcal = 1.290;
0226 double rHcal = 1.9;
0227 double rHo = 3.82;
0228 if (abs(mu_eta) > 0.3)
0229 rHo = 4.07;
0230
0231 double zFaceEcal = 3.209;
0232 if (mu_eta < 0)
0233 zFaceEcal = -1 * zFaceEcal;
0234
0235 double zFaceHcal = 3.88;
0236 if (mu_eta < 0)
0237 zFaceHcal = -1 * zFaceHcal;
0238
0239
0240
0241 double bendr = mu_pt * 1000 / (300 * bfield);
0242
0243 double tb_ecal = TMath::ACos(1 - rEcal * rEcal / (2 * bendr * bendr));
0244 double tb_hcal = TMath::ACos(1 - rHcal * rHcal / (2 * bendr * bendr));
0245 double tb_ho = TMath::ACos(1 - rHo * rHo / (2 * bendr * bendr));
0246 double xEcal, yEcal, zEcal;
0247 double xHcal, yHcal, zHcal;
0248 double xHo, yHo, zHo;
0249
0250
0251 if (fabs(mu_pz * bendr * tb_ecal / mu_pt + mu_vz) < fabs(zFaceEcal) && rEcal < 2 * bendr) {
0252 xEcal = bendr * (TMath::Sin(tb_ecal + mu_phi) - TMath::Sin(mu_phi));
0253 yEcal = bendr * (-TMath::Cos(tb_ecal + mu_phi) + TMath::Cos(mu_phi));
0254 zEcal = bendr * tb_ecal * mu_pz / mu_pt + mu_vz;
0255 } else {
0256 if (mu_pz > 0) {
0257 double te_ecal = (fabs(zFaceEcal) - mu_vz) * mu_pt / (bendr * mu_pz);
0258 xEcal = bendr * (TMath::Sin(te_ecal + mu_phi) - TMath::Sin(mu_phi));
0259 yEcal = bendr * (-TMath::Cos(te_ecal + mu_phi) + TMath::Cos(mu_phi));
0260 zEcal = fabs(zFaceEcal);
0261 } else {
0262 double te_ecal = -(fabs(zFaceEcal) + mu_vz) * mu_pt / (bendr * mu_pz);
0263 xEcal = bendr * (TMath::Sin(te_ecal + mu_phi) - TMath::Sin(mu_phi));
0264 yEcal = bendr * (-TMath::Cos(te_ecal + mu_phi) + TMath::Cos(mu_phi));
0265 zEcal = -fabs(zFaceEcal);
0266 }
0267 }
0268
0269
0270 if (fabs(mu_pz * bendr * tb_hcal / mu_pt + mu_vz) < fabs(zFaceHcal) && rEcal < 2 * bendr) {
0271 xHcal = bendr * (TMath::Sin(tb_hcal + mu_phi) - TMath::Sin(mu_phi));
0272 yHcal = bendr * (-TMath::Cos(tb_hcal + mu_phi) + TMath::Cos(mu_phi));
0273 zHcal = bendr * tb_hcal * mu_pz / mu_pt + mu_vz;
0274 } else {
0275 if (mu_pz > 0) {
0276 double te_hcal = (fabs(zFaceHcal) - mu_vz) * mu_pt / (bendr * mu_pz);
0277 xHcal = bendr * (TMath::Sin(te_hcal + mu_phi) - TMath::Sin(mu_phi));
0278 yHcal = bendr * (-TMath::Cos(te_hcal + mu_phi) + TMath::Cos(mu_phi));
0279 zHcal = fabs(zFaceHcal);
0280 } else {
0281 double te_hcal = -(fabs(zFaceHcal) + mu_vz) * mu_pt / (bendr * mu_pz);
0282 xHcal = bendr * (TMath::Sin(te_hcal + mu_phi) - TMath::Sin(mu_phi));
0283 yHcal = bendr * (-TMath::Cos(te_hcal + mu_phi) + TMath::Cos(mu_phi));
0284 zHcal = -fabs(zFaceHcal);
0285 }
0286 }
0287
0288
0289 xHo = bendr * (TMath::Sin(tb_ho + mu_phi) - TMath::Sin(mu_phi));
0290 yHo = bendr * (-TMath::Cos(tb_ho + mu_phi) + TMath::Cos(mu_phi));
0291 zHo = bendr * tb_ho * mu_pz / mu_pt + mu_vz;
0292
0293 ecalTheta = TMath::ACos(zEcal / sqrt(pow(xEcal, 2) + pow(yEcal, 2) + pow(zEcal, 2)));
0294 ecalPhi = atan2(yEcal, xEcal);
0295 hcalTheta = TMath::ACos(zHcal / sqrt(pow(xHcal, 2) + pow(yHcal, 2) + pow(zHcal, 2)));
0296 hcalPhi = atan2(yHcal, xHcal);
0297 hoTheta = TMath::ACos(zHo / sqrt(pow(xHo, 2) + pow(yHo, 2) + pow(zHo, 2)));
0298 hoPhi = atan2(yHo, xHo);
0299
0300
0301
0302
0303
0304
0305 if (muonCharge > 0) {
0306
0307 double dphi = mu_phi - ecalPhi;
0308 if (fabs(dphi) > TMath::Pi())
0309 dphi = 2 * TMath::Pi() - fabs(dphi);
0310 ecalPhi = mu_phi - fabs(dphi);
0311 if (fabs(ecalPhi) > TMath::Pi()) {
0312 double temp = 2 * TMath::Pi() - fabs(ecalPhi);
0313 ecalPhi = -1 * temp * ecalPhi / fabs(ecalPhi);
0314 }
0315
0316
0317 dphi = mu_phi - hcalPhi;
0318 if (fabs(dphi) > TMath::Pi())
0319 dphi = 2 * TMath::Pi() - fabs(dphi);
0320 hcalPhi = mu_phi - fabs(dphi);
0321 if (fabs(hcalPhi) > TMath::Pi()) {
0322 double temp = 2 * TMath::Pi() - fabs(hcalPhi);
0323 hcalPhi = -1 * temp * hcalPhi / fabs(hcalPhi);
0324 }
0325
0326
0327 dphi = mu_phi - hoPhi;
0328 if (fabs(dphi) > TMath::Pi())
0329 dphi = 2 * TMath::Pi() - fabs(dphi);
0330 hoPhi = mu_phi - fabs(dphi);
0331 if (fabs(hoPhi) > TMath::Pi()) {
0332 double temp = 2 * TMath::Pi() - fabs(hoPhi);
0333 hoPhi = -1 * temp * hoPhi / fabs(hoPhi);
0334 }
0335 }
0336 }
0337
0338
0339 if (!muonMETInfo.useAverage) {
0340 double mu_Ex = muonMETInfo.ecalE * sin(ecalTheta) * cos(ecalPhi) +
0341 muonMETInfo.hcalE * sin(hcalTheta) * cos(hcalPhi) + muonMETInfo.hoE * sin(hoTheta) * cos(hoPhi);
0342 double mu_Ey = muonMETInfo.ecalE * sin(ecalTheta) * sin(ecalPhi) +
0343 muonMETInfo.hcalE * sin(hcalTheta) * sin(hcalPhi) + muonMETInfo.hoE * sin(hoTheta) * sin(hoPhi);
0344
0345 deltax += mu_Ex;
0346 deltay += mu_Ey;
0347
0348 } else {
0349
0350
0351
0352
0353
0354
0355 double dEdx_normalization = -(11.4 + 0.96 * fabs(log(50 * 2.8)) + 0.033 * 50 * (1.0 - pow(50, -0.33))) * 1e-3;
0356 double dEdx_numerator = -(11.4 + 0.96 * fabs(log(mu_p * 2.8)) + 0.033 * mu_p * (1.0 - pow(mu_p, -0.33))) * 1e-3;
0357
0358 double temp = 0.0;
0359
0360 if (muonMETInfo.useHO) {
0361
0362 if (fabs(mu_eta) < 0.2)
0363 temp = 2.75 * (1 - 0.00003 * mu_p);
0364 if (fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
0365 temp = (2.38 + 0.0144 * fabs(mu_eta)) * (1 - 0.0003 * mu_p);
0366 if (fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0367 temp = 7.413 - 5.12 * fabs(mu_eta);
0368 if (fabs(mu_eta) > 1.3)
0369 temp = 2.084 - 0.743 * fabs(mu_eta);
0370 } else {
0371 if (fabs(mu_eta) < 1.0)
0372 temp = 2.33 * (1 - 0.0004 * mu_p);
0373 if (fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0374 temp = (7.413 - 5.12 * fabs(mu_eta)) * (1 - 0.0003 * mu_p);
0375 if (fabs(mu_eta) > 1.3)
0376 temp = 2.084 - 0.743 * fabs(mu_eta);
0377 }
0378
0379 double dep = temp * dEdx_normalization / dEdx_numerator;
0380 if (dep < 0.5)
0381 dep = 0;
0382
0383 if (fabs(mu_eta) < 1.3) {
0384 deltax += dep * cos((ecalPhi + hcalPhi + hoPhi) / 3);
0385 deltay += dep * sin((ecalPhi + hcalPhi + hoPhi) / 3);
0386 } else {
0387 deltax += dep * cos((ecalPhi + hcalPhi) / 2);
0388 deltay += dep * cos((ecalPhi + hcalPhi) / 2);
0389 }
0390 }
0391 }
0392
0393
0394 void MuonMETAlgo::run(const edm::View<reco::Muon>& inputMuons,
0395 const edm::ValueMap<reco::MuonMETCorrectionData>& vm_muCorrData,
0396 const edm::View<reco::MET>& uncorMET,
0397 METCollection* corMET) {
0398 MuonMETAlgo_run(inputMuons, vm_muCorrData, uncorMET, corMET);
0399 }
0400
0401
0402 void MuonMETAlgo::run(const edm::View<reco::Muon>& inputMuons,
0403 const edm::ValueMap<reco::MuonMETCorrectionData>& vm_muCorrData,
0404 const edm::View<reco::CaloMET>& uncorMET,
0405 CaloMETCollection* corMET) {
0406 MuonMETAlgo_run(inputMuons, vm_muCorrData, uncorMET, corMET);
0407 }
0408
0409