CMS 3D CMS Logo

/data/git/CMSSW_5_3_11_patch5/src/RecoMuon/TrackingTools/src/MuonTrajectoryCleaner.cc

Go to the documentation of this file.
00001 
00008 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryCleaner.h"
00009 #include "DataFormats/TrackReco/interface/Track.h"
00010 #include "DataFormats/Math/interface/deltaR.h"
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012 #include "DataFormats/Common/interface/RefToBase.h"
00013 #include "DataFormats/Common/interface/AssociationMap.h"
00014 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
00015 #include "DataFormats/TrackingRecHit/interface/RecSegment.h"                                // Daniele
00016 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"  // Daniele
00017 
00018 using namespace std; 
00019 
00020 typedef edm::AssociationMap<edm::OneToMany<L2MuonTrajectorySeedCollection, L2MuonTrajectorySeedCollection> > L2SeedAssoc;
00021 
00022 void MuonTrajectoryCleaner::clean(TrajectoryContainer& trajC, edm::Event& event) {
00023   const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
00024 
00025   LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
00026 
00027   TrajectoryContainer::iterator iter, jter;
00028   Trajectory::DataContainer::const_iterator m1, m2;
00029 
00030   if ( trajC.size() < 1 ) return;
00031   
00032   LogTrace(metname) << "Number of trajectories in the container: " <<trajC.size()<< endl;
00033 
00034   int i(0), j(0);
00035   int match(0);
00036   int nTot1DHits_i(0), nTot1DHits_j(0);  // tot number of 1D/2D hits (including invalid)
00037 
00038   // Map between chosen seed and ghost seeds (for trigger)
00039   map<int, vector<int> > seedmap;
00040 
00041   // CAVEAT: vector<bool> is not a vector, its elements are not addressable!
00042   // This is fine as long as only operator [] is used as in this case.
00043   // cf. par 16.3.11
00044   vector<bool> mask(trajC.size(),true);
00045   
00046   TrajectoryContainer result;
00047   
00048   for ( iter = trajC.begin(); iter != trajC.end(); iter++ ) {
00049     if ( !mask[i] ) { i++; continue; }
00050     if ( !(*iter)->isValid() || (*iter)->empty() ) {mask[i] = false; i++; continue; }
00051     if(seedmap.count(i)==0) seedmap[i].push_back(i);
00052     const Trajectory::DataContainer& meas1 = (*iter)->measurements();
00053     j = i+1;
00054     bool skipnext=false;
00055     for ( jter = iter+1; jter != trajC.end(); jter++ ) {
00056       if ( !mask[j] ) { j++; continue; }
00057       if(seedmap.count(j)==0) seedmap[j].push_back(j);
00058       const Trajectory::DataContainer& meas2 = (*jter)->measurements();
00059       match = 0;
00060       nTot1DHits_i = 0;
00061       for( m1 = meas1.begin(); m1 != meas1.end(); m1++ ) {
00062         TransientTrackingRecHit::ConstRecHitContainer trhC1;
00063         if((*m1).recHit()->dimension()==4) trhC1 = MuonTransientTrackingRecHitBreaker::breakInSubRecHits((*m1).recHit(), 2);
00064         else trhC1.push_back((*m1).recHit());
00065         for( TransientTrackingRecHit::ConstRecHitContainer::const_iterator trh1=trhC1.begin(); trh1!=trhC1.end(); ++trh1, ++nTot1DHits_i) {
00066           nTot1DHits_j = 0;
00067           for( m2 = meas2.begin(); m2 != meas2.end(); m2++ ) {
00068             TransientTrackingRecHit::ConstRecHitContainer trhC2;
00069             if((*m2).recHit()->dimension()==4) trhC2 = MuonTransientTrackingRecHitBreaker::breakInSubRecHits((*m2).recHit(), 2);
00070             else trhC2.push_back((*m2).recHit());
00071             for( TransientTrackingRecHit::ConstRecHitContainer::const_iterator trh2=trhC2.begin(); trh2!=trhC2.end(); ++trh2, ++nTot1DHits_j) {
00072               if( ( (*trh1)->globalPosition() - (*trh2)->globalPosition()).mag() < 10e-5 ) match++;
00073               //          if ( ( (*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag()< 10e-5 ) match++;
00074             } // end for( trh2 ... )
00075           } // end for( m2 ... )
00076         } // end for( trh1 ... )
00077       } // end for( m1 ... )
00078       
00079 
00080       int nTotHits_i = (*iter)->measurements().size();
00081       int nTotHits_j = (*jter)->measurements().size();
00082 
00083       // FIXME Set Boff/on via cfg!
00084       double chi2_dof_i = (*iter)->ndof() > 0 ? (*iter)->chiSquared()/(*iter)->ndof() : (*iter)->chiSquared()/1e-10;
00085       double chi2_dof_j = (*jter)->ndof() > 0 ? (*jter)->chiSquared()/(*jter)->ndof() : (*jter)->chiSquared()/1e-10;
00086 
00087       LogTrace(metname) << " MuonTrajectoryCleaner:";
00088       LogTrace(metname) << " * trajC " << i 
00089                         << " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp() 
00090                         << " GeV) - chi2/nDOF = " << (*iter)->chiSquared() << "/" << (*iter)->ndof() << " = " << chi2_dof_i;
00091       LogTrace(metname) << "     - valid RH = " << (*iter)->foundHits() << " / total RH = " <<  nTotHits_i << " / total 1D RH = " << nTot1DHits_i;
00092       LogTrace(metname) << " * trajC " << j 
00093                         << " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp() 
00094                         << " GeV) - chi2/nDOF = " << (*jter)->chiSquared() << "/" << (*jter)->ndof() << " = " << chi2_dof_j;
00095       LogTrace(metname) << "     - valid RH = " << (*jter)->foundHits() << " / total RH = " <<  nTotHits_j << " / total 1D RH = " << nTot1DHits_j;
00096       LogTrace(metname) << " *** Shared 1D RecHits: " << match;
00097 
00098       int hit_diff =  (*iter)->foundHits() - (*jter)->foundHits() ;       
00099       // If there are matches, reject the worst track
00100       if ( match > 0 ) {
00101         // If the difference of # of rechits is less than 4, compare the chi2/ndf
00102         if ( abs(hit_diff) <= 4  ) {
00103 
00104           double minPt = 3.5;
00105           double dPt = 7.;  // i.e. considering 10% (conservative!) resolution at minPt it is ~ 10 sigma away from the central value
00106 
00107           double maxFraction = 0.95;
00108 
00109           double fraction = (2.*match)/((*iter)->measurements().size()+(*jter)->measurements().size());
00110           int belowLimit = 0;
00111           int above = 0;
00112 
00113           if((*jter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt) ++belowLimit; 
00114           if((*iter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt) ++belowLimit; 
00115          
00116           if((*jter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt) ++above; 
00117           if((*iter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt) ++above; 
00118           
00119 
00120           if(fraction >= maxFraction && belowLimit == 1 && above == 1){
00121             if((*iter)->lastMeasurement().updatedState().globalMomentum().perp() < minPt){
00122               mask[i] = false;
00123               skipnext=true;
00124               seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
00125               seedmap.erase(i);
00126               LogTrace(metname) << "Trajectory # " << i << " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp() 
00127                                 << " GeV) rejected because it has too low pt";
00128             } 
00129             else {
00130               mask[j] = false;
00131               seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
00132               seedmap.erase(j);
00133               LogTrace(metname) << "Trajectory # " << j << " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp() 
00134                                 << " GeV) rejected because it has too low pt";
00135             }
00136           }
00137           else{   
00138             if (chi2_dof_i  > chi2_dof_j) {
00139               mask[i] = false;
00140               skipnext=true;
00141               seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
00142               seedmap.erase(i);
00143               LogTrace(metname) << "Trajectory # " << i << " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
00144             }
00145             else{
00146               mask[j] = false;
00147               seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
00148               seedmap.erase(j);
00149               LogTrace(metname) << "Trajectory # " << j << " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
00150             }
00151           }
00152         }
00153         else { // different number of hits
00154           if ( hit_diff < 0 ) {
00155             mask[i] = false;
00156             skipnext=true;
00157             seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
00158             seedmap.erase(i);
00159             LogTrace(metname) << "Trajectory # " << i << " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
00160           }
00161           else { 
00162             mask[j] = false;
00163             seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
00164             seedmap.erase(j);
00165             LogTrace(metname) << "Trajectory # " << j << " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
00166           }
00167         } 
00168       }
00169       if(skipnext) break;
00170       j++;
00171     }
00172     i++;
00173     if(skipnext) continue;
00174   }
00175   
00176 
00177   // Association map between the seed of the chosen trajectory and the seeds of the discarded trajectories
00178   if(reportGhosts_) {
00179     LogTrace(metname) << " Creating map between chosen seed and ghost seeds." << std::endl;
00180 
00181     auto_ptr<L2SeedAssoc> seedToSeedsMap(new L2SeedAssoc);
00182     int seedcnt(0);
00183 
00184     for(map<int, vector<int> >::iterator itmap=seedmap.begin(); itmap!=seedmap.end(); ++itmap, ++seedcnt) {
00185       edm::RefToBase<TrajectorySeed> tmpSeedRef1 = trajC[(*itmap).first]->seedRef();
00186       edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef1 = tmpSeedRef1.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
00187 
00188       int tmpsize=(*itmap).second.size();
00189       for(int cnt=0; cnt!=tmpsize; ++cnt) {
00190         edm::RefToBase<TrajectorySeed> tmpSeedRef2 = trajC[(*itmap).second[cnt]]->seedRef();
00191         edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef2 = tmpSeedRef2.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
00192         seedToSeedsMap->insert(tmpL2SeedRef1, tmpL2SeedRef2);
00193       }
00194     }
00195 
00196     event.put(seedToSeedsMap, "");
00197 
00198   } // end if(reportGhosts_)
00199 
00200   i = 0;
00201   for ( iter = trajC.begin(); iter != trajC.end(); iter++ ) {
00202     if ( mask[i] ){
00203       result.push_back(*iter);
00204       LogTrace(metname) << "Keep trajectory with pT = " << (*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV";
00205     }
00206     else delete *iter;
00207     i++;
00208   }
00209 
00210   trajC.clear();
00211   trajC = result;
00212 }
00213 
00214 //
00215 // clean CandidateContainer
00216 //
00217 void MuonTrajectoryCleaner::clean(CandidateContainer& candC){ 
00218   const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
00219 
00220   LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
00221 
00222   if ( candC.size() < 2 ) return;
00223 
00224   CandidateContainer::iterator iter, jter;
00225   Trajectory::DataContainer::const_iterator m1, m2;
00226 
00227   const float deltaEta = 0.01;
00228   const float deltaPhi = 0.01;
00229   const float deltaPt  = 1.0;
00230   
00231   LogTrace(metname) << "Number of muon candidates in the container: " <<candC.size()<< endl;
00232 
00233   int i(0), j(0);
00234   int match(0);
00235   //unused  bool directionMatch = false;
00236 
00237   // CAVEAT: vector<bool> is not a vector, its elements are not addressable!
00238   // This is fine as long as only operator [] is used as in this case.
00239   // cf. par 16.3.11
00240   vector<bool> mask(candC.size(),true);
00241   
00242   CandidateContainer result;
00243   
00244   for ( iter = candC.begin(); iter != candC.end(); iter++ ) {
00245     if ( !mask[i] ) { i++; continue; }
00246     const Trajectory::DataContainer& meas1 = (*iter)->trajectory()->measurements();
00247     j = i+1;
00248     bool skipnext=false;
00249 
00250     TrajectoryStateOnSurface innerTSOS;
00251 
00252     if ((*iter)->trajectory()->direction() == alongMomentum) {
00253       innerTSOS = (*iter)->trajectory()->firstMeasurement().updatedState();
00254     } 
00255     else if ((*iter)->trajectory()->direction() == oppositeToMomentum) { 
00256       innerTSOS = (*iter)->trajectory()->lastMeasurement().updatedState();
00257          }
00258     if ( !(innerTSOS.isValid()) ) continue;
00259     float pt1 = innerTSOS.globalMomentum().perp();
00260     float eta1 = innerTSOS.globalMomentum().eta();
00261     float phi1 = innerTSOS.globalMomentum().phi();
00262 
00263     for ( jter = iter+1; jter != candC.end(); jter++ ) {
00264       if ( !mask[j] ) { j++; continue; }
00265       //UNUSED:      directionMatch = false;
00266       const Trajectory::DataContainer& meas2 = (*jter)->trajectory()->measurements();
00267       match = 0;
00268       for ( m1 = meas1.begin(); m1 != meas1.end(); m1++ ) {
00269         for ( m2 = meas2.begin(); m2 != meas2.end(); m2++ ) {
00270           if ( (*m1).recHit()->isValid() && (*m2).recHit()->isValid() ) 
00271             if ( ( (*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag()< 10e-5 ) match++;
00272         }
00273       }
00274       
00275       LogTrace(metname) 
00276         << " MuonTrajectoryCleaner: candC " << i << " chi2/nRH = " 
00277         << (*iter)->trajectory()->chiSquared() << "/" << (*iter)->trajectory()->foundHits() <<
00278         " vs trajC " << j << " chi2/nRH = " << (*jter)->trajectory()->chiSquared() <<
00279         "/" << (*jter)->trajectory()->foundHits() << " Shared RecHits: " << match;
00280 
00281       TrajectoryStateOnSurface innerTSOS2;       
00282       if ((*jter)->trajectory()->direction() == alongMomentum) {
00283         innerTSOS2 = (*jter)->trajectory()->firstMeasurement().updatedState();
00284       }
00285       else if ((*jter)->trajectory()->direction() == oppositeToMomentum) {
00286         innerTSOS2 = (*jter)->trajectory()->lastMeasurement().updatedState();
00287       }
00288       if ( !(innerTSOS2.isValid()) ) continue;
00289 
00290       float pt2 = innerTSOS2.globalMomentum().perp();
00291       float eta2 = innerTSOS2.globalMomentum().eta();
00292       float phi2 = innerTSOS2.globalMomentum().phi();
00293 
00294       float deta(fabs(eta1-eta2));
00295       float dphi(fabs(Geom::Phi<float>(phi1)-Geom::Phi<float>(phi2)));
00296       float dpt(fabs(pt1-pt2));
00297       if ( dpt < deltaPt && deta < deltaEta && dphi < deltaPhi ) {
00298         //UNUSED:        directionMatch = true;
00299         LogTrace(metname)
00300         << " MuonTrajectoryCleaner: candC " << i<<" and "<<j<< " direction matched: "
00301         <<innerTSOS.globalMomentum()<<" and " <<innerTSOS2.globalMomentum();
00302       }
00303 
00304       // If there are matches, reject the worst track
00305       bool hitsMatch = ((match > 0) && ((match/((*iter)->trajectory()->foundHits()) > 0.25) || (match/((*jter)->trajectory()->foundHits()) > 0.25))) ? true : false;
00306       bool tracksMatch = 
00307       ( ( (*iter)->trackerTrack() == (*jter)->trackerTrack() ) && 
00308         ( deltaR<double>((*iter)->muonTrack()->eta(),(*iter)->muonTrack()->phi(), (*jter)->muonTrack()->eta(),(*jter)->muonTrack()->phi()) < 0.2) );
00309 
00310       //if ( ( tracksMatch ) || (hitsMatch > 0) || directionMatch ) {
00311                         if ( ( tracksMatch ) || (hitsMatch > 0) ) { 
00312         if (  (*iter)->trajectory()->foundHits() == (*jter)->trajectory()->foundHits() ) {
00313           if ( (*iter)->trajectory()->chiSquared() > (*jter)->trajectory()->chiSquared() ) {
00314             mask[i] = false;
00315             skipnext=true;
00316           }
00317           else mask[j] = false;
00318         }
00319         else { // different number of hits
00320           if ( (*iter)->trajectory()->foundHits() < (*jter)->trajectory()->foundHits() ) {
00321             mask[i] = false;
00322             skipnext=true;
00323           }
00324           else mask[j] = false;
00325         }
00326       }
00327       if(skipnext) break;
00328       j++;
00329     }
00330     i++;
00331     if(skipnext) continue;
00332   }
00333   
00334   i = 0;
00335   for ( iter = candC.begin(); iter != candC.end(); iter++ ) {
00336     if ( mask[i] ) {
00337        result.push_back(*iter);
00338     } else {
00339        delete (*iter)->trajectory();
00340        delete (*iter)->trackerTrajectory();
00341        delete *iter;
00342     } 
00343     i++;
00344   }
00345   
00346   candC.clear();
00347   candC = result;
00348 }
00349