Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:18 

0001 /**
0002  *  A selector for muon tracks
0003  *
0004  *  \author R.Bellan - INFN Torino
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   // Map between chosen seed and ghost seeds (for trigger)
0039   map<int, vector<int> > seedmap;
0040 
0041   // CAVEAT: vector<bool> is not a vector, its elements are not addressable!
0042   // This is fine as long as only operator [] is used as in this case.
0043   // cf. par 16.3.11
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         }  // end for( m2 ... )
0077       }    // end for( m1 ... )
0078 
0079       // FIXME Set Boff/on via cfg!
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       // If there are matches, reject the worst track
0100       if (match > 0) {
0101         // If the difference in # of rechits is null then compare the chi2/ndf of the two trajectories
0102         if (abs(hit_diff) == 0) {
0103           double minPt = 3.5;
0104           double dPt =
0105               7.;  // i.e. considering 10% (conservative!) resolution at minPt it is ~ 10 sigma away from the central value
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 {  // different number of hits
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   // Association map between the seed of the chosen trajectory and the seeds of the discarded trajectories
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   }  // end if(reportGhosts_)
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 // clean CandidateContainer
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   //unused  bool directionMatch = false;
0256 
0257   // CAVEAT: vector<bool> is not a vector, its elements are not addressable!
0258   // This is fine as long as only operator [] is used as in this case.
0259   // cf. par 16.3.11
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       //UNUSED:      directionMatch = false;
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         //UNUSED:        directionMatch = true;
0326         LogTrace(metname) << " MuonTrajectoryCleaner: candC " << i << " and " << j
0327                           << " direction matched: " << innerTSOS.globalMomentum() << " and "
0328                           << innerTSOS2.globalMomentum();
0329       }
0330 
0331       // If there are matches, reject the worst track
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       //if ( ( tracksMatch ) || (hitsMatch > 0) || directionMatch ) {
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 {  // different number of hits
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!