File indexing completed on 2024-09-07 04:37:48
0001
0002
0003
0004
0005
0006 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryCleaner.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "DataFormats/Common/interface/RefToBase.h"
0011 #include "DataFormats/Common/interface/AssociationMap.h"
0012 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
0013 #include "DataFormats/TrackingRecHit/interface/RecSegment.h"
0014 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
0015
0016 using namespace std;
0017
0018 typedef edm::AssociationMap<edm::OneToMany<L2MuonTrajectorySeedCollection, L2MuonTrajectorySeedCollection> > L2SeedAssoc;
0019
0020 void MuonTrajectoryCleaner::clean(TrajectoryContainer& trajC,
0021 edm::Event& event,
0022 const edm::Handle<edm::View<TrajectorySeed> >& seeds) {
0023 const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
0024
0025 LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
0026
0027 TrajectoryContainer::iterator iter, jter;
0028 Trajectory::DataContainer::const_iterator m1, m2;
0029
0030 if (trajC.empty())
0031 return;
0032
0033 LogTrace(metname) << "Number of trajectories in the container: " << trajC.size() << endl;
0034
0035 int i(0), j(0);
0036 int match(0);
0037
0038
0039 map<int, vector<int> > seedmap;
0040
0041
0042
0043
0044 vector<bool> mask(trajC.size(), true);
0045
0046 TrajectoryContainer result;
0047
0048 for (iter = trajC.begin(); iter != trajC.end(); iter++) {
0049 if (!mask[i]) {
0050 i++;
0051 continue;
0052 }
0053 if (!(*iter)->isValid() || (*iter)->empty()) {
0054 mask[i] = false;
0055 i++;
0056 continue;
0057 }
0058 if (seedmap.count(i) == 0)
0059 seedmap[i].push_back(i);
0060 const Trajectory::DataContainer& meas1 = (*iter)->measurements();
0061 j = i + 1;
0062 bool skipnext = false;
0063 for (jter = iter + 1; jter != trajC.end(); jter++) {
0064 if (!mask[j]) {
0065 j++;
0066 continue;
0067 }
0068 if (seedmap.count(j) == 0)
0069 seedmap[j].push_back(j);
0070 const Trajectory::DataContainer& meas2 = (*jter)->measurements();
0071 match = 0;
0072 for (m1 = meas1.begin(); m1 != meas1.end(); m1++) {
0073 for (m2 = meas2.begin(); m2 != meas2.end(); m2++) {
0074 if (((*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag() < 10e-5)
0075 match++;
0076 }
0077 }
0078
0079
0080 double chi2_dof_i = (*iter)->ndof() > 0 ? (*iter)->chiSquared() / (*iter)->ndof() : (*iter)->chiSquared() / 1e-10;
0081 double chi2_dof_j = (*jter)->ndof() > 0 ? (*jter)->chiSquared() / (*jter)->ndof() : (*jter)->chiSquared() / 1e-10;
0082
0083 LogTrace(metname) << " MuonTrajectoryCleaner:";
0084 LogTrace(metname) << " * trajC " << i
0085 << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
0086 << " GeV) - chi2/nDOF = " << (*iter)->chiSquared() << "/" << (*iter)->ndof() << " = "
0087 << chi2_dof_i;
0088 LogTrace(metname) << " - valid RH = " << (*iter)->foundHits()
0089 << " / total RH = " << (*iter)->measurements().size();
0090 LogTrace(metname) << " * trajC " << j
0091 << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
0092 << " GeV) - chi2/nDOF = " << (*jter)->chiSquared() << "/" << (*jter)->ndof() << " = "
0093 << chi2_dof_j;
0094 LogTrace(metname) << " - valid RH = " << (*jter)->foundHits()
0095 << " / total RH = " << (*jter)->measurements().size();
0096 LogTrace(metname) << " *** Shared RecHits: " << match;
0097
0098 int hit_diff = (*iter)->foundHits() - (*jter)->foundHits();
0099
0100 if (match > 0) {
0101
0102 if (abs(hit_diff) == 0) {
0103 double minPt = 3.5;
0104 double dPt =
0105 7.;
0106
0107 double maxFraction = 0.95;
0108
0109 double fraction = (2. * match) / ((*iter)->measurements().size() + (*jter)->measurements().size());
0110 int belowLimit = 0;
0111 int above = 0;
0112
0113 if ((*jter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt)
0114 ++belowLimit;
0115 if ((*iter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt)
0116 ++belowLimit;
0117
0118 if ((*jter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt)
0119 ++above;
0120 if ((*iter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt)
0121 ++above;
0122
0123 if (fraction >= maxFraction && belowLimit == 1 && above == 1) {
0124 if ((*iter)->lastMeasurement().updatedState().globalMomentum().perp() < minPt) {
0125 mask[i] = false;
0126 skipnext = true;
0127 seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
0128 seedmap.erase(i);
0129 LogTrace(metname) << "Trajectory # " << i
0130 << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
0131 << " GeV) rejected because it has too low pt";
0132 } else {
0133 mask[j] = false;
0134 seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
0135 seedmap.erase(j);
0136 LogTrace(metname) << "Trajectory # " << j
0137 << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
0138 << " GeV) rejected because it has too low pt";
0139 }
0140 } else {
0141 if (chi2_dof_i > chi2_dof_j) {
0142 mask[i] = false;
0143 skipnext = true;
0144 seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
0145 seedmap.erase(i);
0146 LogTrace(metname) << "Trajectory # " << i
0147 << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
0148 << " GeV) rejected";
0149 } else {
0150 mask[j] = false;
0151 seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
0152 seedmap.erase(j);
0153 LogTrace(metname) << "Trajectory # " << j
0154 << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
0155 << " GeV) rejected";
0156 }
0157 }
0158 } else {
0159 if (hit_diff < 0) {
0160 mask[i] = false;
0161 skipnext = true;
0162 seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
0163 seedmap.erase(i);
0164 LogTrace(metname) << "Trajectory # " << i
0165 << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
0166 << " GeV) rejected";
0167 } else {
0168 mask[j] = false;
0169 seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
0170 seedmap.erase(j);
0171 LogTrace(metname) << "Trajectory # " << j
0172 << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
0173 << " GeV) rejected";
0174 }
0175 }
0176 }
0177 if (skipnext)
0178 break;
0179 j++;
0180 }
0181 i++;
0182 if (skipnext)
0183 continue;
0184 }
0185
0186
0187 if (reportGhosts_) {
0188 LogTrace(metname) << " Creating map between chosen seed and ghost seeds." << std::endl;
0189
0190 auto seedToSeedsMap = std::make_unique<L2SeedAssoc>();
0191 if (!seeds->empty()) {
0192 edm::Ptr<TrajectorySeed> ptr = seeds->ptrAt(0);
0193 edm::Handle<L2MuonTrajectorySeedCollection> seedsHandle;
0194 event.get(ptr.id(), seedsHandle);
0195 seedToSeedsMap = std::make_unique<L2SeedAssoc>(seedsHandle, seedsHandle);
0196 }
0197
0198 int seedcnt(0);
0199
0200 for (map<int, vector<int> >::iterator itmap = seedmap.begin(); itmap != seedmap.end(); ++itmap, ++seedcnt) {
0201 edm::RefToBase<TrajectorySeed> tmpSeedRef1 = trajC[(*itmap).first]->seedRef();
0202 edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef1 =
0203 tmpSeedRef1.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
0204
0205 int tmpsize = (*itmap).second.size();
0206 for (int cnt = 0; cnt != tmpsize; ++cnt) {
0207 edm::RefToBase<TrajectorySeed> tmpSeedRef2 = trajC[(*itmap).second[cnt]]->seedRef();
0208 edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef2 =
0209 tmpSeedRef2.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
0210 seedToSeedsMap->insert(tmpL2SeedRef1, tmpL2SeedRef2);
0211 }
0212 }
0213
0214 event.put(std::move(seedToSeedsMap), "");
0215
0216 }
0217
0218 i = 0;
0219 auto c = std::count(mask.begin(), mask.end(), true);
0220 result.reserve(c);
0221 for (iter = trajC.begin(); iter != trajC.end(); iter++) {
0222 if (mask[i]) {
0223 LogTrace(metname) << "Keep trajectory with pT = "
0224 << (*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV";
0225 result.push_back(std::move(*iter));
0226 }
0227 i++;
0228 }
0229
0230 trajC = std::move(result);
0231 }
0232
0233
0234
0235
0236 void MuonTrajectoryCleaner::clean(CandidateContainer& candC) {
0237 const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
0238
0239 LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
0240
0241 if (candC.size() < 2)
0242 return;
0243
0244 CandidateContainer::iterator iter, jter;
0245 Trajectory::DataContainer::const_iterator m1, m2;
0246
0247 const float deltaEta = 0.01;
0248 const float deltaPhi = 0.01;
0249 const float deltaPt = 1.0;
0250
0251 LogTrace(metname) << "Number of muon candidates in the container: " << candC.size() << endl;
0252
0253 int i(0), j(0);
0254 int match(0);
0255
0256
0257
0258
0259
0260 vector<bool> mask(candC.size(), true);
0261
0262 CandidateContainer result;
0263
0264 for (iter = candC.begin(); iter != candC.end(); iter++) {
0265 if (!mask[i]) {
0266 i++;
0267 continue;
0268 }
0269 const Trajectory::DataContainer& meas1 = (*iter)->trajectory()->measurements();
0270 j = i + 1;
0271 bool skipnext = false;
0272
0273 TrajectoryStateOnSurface innerTSOS;
0274
0275 if ((*iter)->trajectory()->direction() == alongMomentum) {
0276 innerTSOS = (*iter)->trajectory()->firstMeasurement().updatedState();
0277 } else if ((*iter)->trajectory()->direction() == oppositeToMomentum) {
0278 innerTSOS = (*iter)->trajectory()->lastMeasurement().updatedState();
0279 }
0280 if (!(innerTSOS.isValid()))
0281 continue;
0282 float pt1 = innerTSOS.globalMomentum().perp();
0283 float eta1 = innerTSOS.globalMomentum().eta();
0284 float phi1 = innerTSOS.globalMomentum().phi();
0285
0286 for (jter = iter + 1; jter != candC.end(); jter++) {
0287 if (!mask[j]) {
0288 j++;
0289 continue;
0290 }
0291
0292 const Trajectory::DataContainer& meas2 = (*jter)->trajectory()->measurements();
0293 match = 0;
0294 for (m1 = meas1.begin(); m1 != meas1.end(); m1++) {
0295 for (m2 = meas2.begin(); m2 != meas2.end(); m2++) {
0296 if ((*m1).recHit()->isValid() && (*m2).recHit()->isValid())
0297 if (((*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag() < 10e-5)
0298 match++;
0299 }
0300 }
0301
0302 LogTrace(metname) << " MuonTrajectoryCleaner: candC " << i
0303 << " chi2/nRH = " << (*iter)->trajectory()->chiSquared() << "/"
0304 << (*iter)->trajectory()->foundHits() << " vs trajC " << j
0305 << " chi2/nRH = " << (*jter)->trajectory()->chiSquared() << "/"
0306 << (*jter)->trajectory()->foundHits() << " Shared RecHits: " << match;
0307
0308 TrajectoryStateOnSurface innerTSOS2;
0309 if ((*jter)->trajectory()->direction() == alongMomentum) {
0310 innerTSOS2 = (*jter)->trajectory()->firstMeasurement().updatedState();
0311 } else if ((*jter)->trajectory()->direction() == oppositeToMomentum) {
0312 innerTSOS2 = (*jter)->trajectory()->lastMeasurement().updatedState();
0313 }
0314 if (!(innerTSOS2.isValid()))
0315 continue;
0316
0317 float pt2 = innerTSOS2.globalMomentum().perp();
0318 float eta2 = innerTSOS2.globalMomentum().eta();
0319 float phi2 = innerTSOS2.globalMomentum().phi();
0320
0321 float deta(fabs(eta1 - eta2));
0322 float dphi(fabs(Geom::Phi<float>(phi1) - Geom::Phi<float>(phi2)));
0323 float dpt(fabs(pt1 - pt2));
0324 if (dpt < deltaPt && deta < deltaEta && dphi < deltaPhi) {
0325
0326 LogTrace(metname) << " MuonTrajectoryCleaner: candC " << i << " and " << j
0327 << " direction matched: " << innerTSOS.globalMomentum() << " and "
0328 << innerTSOS2.globalMomentum();
0329 }
0330
0331
0332 bool hitsMatch = ((match > 0) && ((match / ((*iter)->trajectory()->foundHits()) > 0.25) ||
0333 (match / ((*jter)->trajectory()->foundHits()) > 0.25)))
0334 ? true
0335 : false;
0336 bool tracksMatch =
0337 (((*iter)->trackerTrack() == (*jter)->trackerTrack()) && (deltaR<double>((*iter)->muonTrack()->eta(),
0338 (*iter)->muonTrack()->phi(),
0339 (*jter)->muonTrack()->eta(),
0340 (*jter)->muonTrack()->phi()) < 0.2));
0341
0342
0343 if ((tracksMatch) || (hitsMatch > 0)) {
0344 if ((*iter)->trajectory()->foundHits() == (*jter)->trajectory()->foundHits()) {
0345 if ((*iter)->trajectory()->chiSquared() > (*jter)->trajectory()->chiSquared()) {
0346 mask[i] = false;
0347 skipnext = true;
0348 } else
0349 mask[j] = false;
0350 } else {
0351 if ((*iter)->trajectory()->foundHits() < (*jter)->trajectory()->foundHits()) {
0352 mask[i] = false;
0353 skipnext = true;
0354 } else
0355 mask[j] = false;
0356 }
0357 }
0358 if (skipnext)
0359 break;
0360 j++;
0361 }
0362 i++;
0363 if (skipnext)
0364 continue;
0365 }
0366
0367 auto c = std::count(mask.begin(), mask.end(), true);
0368 i = 0;
0369 result.reserve(c);
0370 for (iter = candC.begin(); iter != candC.end(); iter++) {
0371 if (mask[i]) {
0372 result.push_back(std::move(*iter));
0373 }
0374 i++;
0375 }
0376
0377 candC = std::move(result);
0378 }