|
|
|
#include <MuonTrajectoryCleaner.h>
Public Types | |
| typedef MuonCandidate::CandidateContainer | CandidateContainer |
| typedef MuonCandidate::TrajectoryContainer | TrajectoryContainer |
Public Member Functions | |
| void | clean (TrajectoryContainer &muonTrajectories, edm::Event &evt) |
| Clean the trajectories container, erasing the (worst) clone trajectory. | |
| void | clean (CandidateContainer &muonTrajectories) |
| Clean the candidates container, erasing the (worst) clone trajectory. | |
| MuonTrajectoryCleaner (bool reportGhosts) | |
| Constructor for L2 muons (enable reportGhosts) | |
| MuonTrajectoryCleaner () | |
| Constructor. | |
| virtual | ~MuonTrajectoryCleaner () |
| Destructor. | |
Private Attributes | |
| bool | reportGhosts_ |
No description available.
Definition at line 18 of file MuonTrajectoryCleaner.h.
Definition at line 21 of file MuonTrajectoryCleaner.h.
Definition at line 20 of file MuonTrajectoryCleaner.h.
| MuonTrajectoryCleaner::MuonTrajectoryCleaner | ( | ) | [inline] |
| MuonTrajectoryCleaner::MuonTrajectoryCleaner | ( | bool | reportGhosts | ) | [inline] |
Constructor for L2 muons (enable reportGhosts)
Definition at line 28 of file MuonTrajectoryCleaner.h.
: reportGhosts_(reportGhosts) {}
| virtual MuonTrajectoryCleaner::~MuonTrajectoryCleaner | ( | ) | [inline, virtual] |
| void MuonTrajectoryCleaner::clean | ( | TrajectoryContainer & | muonTrajectories, |
| edm::Event & | evt | ||
| ) |
Clean the trajectories container, erasing the (worst) clone trajectory.
Definition at line 22 of file MuonTrajectoryCleaner.cc.
References abs, begin, MuonTransientTrackingRecHitBreaker::breakInSubRecHits(), edm::RefToBase< T >::castTo(), end, i, j, LogTrace, mag(), match(), metname, gsfElectronCkfTrackCandidateMaker_cff::minPt, and query::result.
Referenced by MuonTrackFinder::reconstruct().
{
const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
TrajectoryContainer::iterator iter, jter;
Trajectory::DataContainer::const_iterator m1, m2;
if ( trajC.size() < 1 ) return;
LogTrace(metname) << "Number of trajectories in the container: " <<trajC.size()<< endl;
int i(0), j(0);
int match(0);
int nTot1DHits_i(0), nTot1DHits_j(0); // tot number of 1D/2D hits (including invalid)
// Map between chosen seed and ghost seeds (for trigger)
map<int, vector<int> > seedmap;
// CAVEAT: vector<bool> is not a vector, its elements are not addressable!
// This is fine as long as only operator [] is used as in this case.
// cf. par 16.3.11
vector<bool> mask(trajC.size(),true);
TrajectoryContainer result;
for ( iter = trajC.begin(); iter != trajC.end(); iter++ ) {
if ( !mask[i] ) { i++; continue; }
if ( !(*iter)->isValid() || (*iter)->empty() ) {mask[i] = false; i++; continue; }
if(seedmap.count(i)==0) seedmap[i].push_back(i);
const Trajectory::DataContainer& meas1 = (*iter)->measurements();
j = i+1;
bool skipnext=false;
for ( jter = iter+1; jter != trajC.end(); jter++ ) {
if ( !mask[j] ) { j++; continue; }
if(seedmap.count(j)==0) seedmap[j].push_back(j);
const Trajectory::DataContainer& meas2 = (*jter)->measurements();
match = 0;
nTot1DHits_i = 0;
for( m1 = meas1.begin(); m1 != meas1.end(); m1++ ) {
TransientTrackingRecHit::ConstRecHitContainer trhC1;
if((*m1).recHit()->dimension()==4) trhC1 = MuonTransientTrackingRecHitBreaker::breakInSubRecHits((*m1).recHit(), 2);
else trhC1.push_back((*m1).recHit());
for( TransientTrackingRecHit::ConstRecHitContainer::const_iterator trh1=trhC1.begin(); trh1!=trhC1.end(); ++trh1, ++nTot1DHits_i) {
nTot1DHits_j = 0;
for( m2 = meas2.begin(); m2 != meas2.end(); m2++ ) {
TransientTrackingRecHit::ConstRecHitContainer trhC2;
if((*m2).recHit()->dimension()==4) trhC2 = MuonTransientTrackingRecHitBreaker::breakInSubRecHits((*m2).recHit(), 2);
else trhC2.push_back((*m2).recHit());
for( TransientTrackingRecHit::ConstRecHitContainer::const_iterator trh2=trhC2.begin(); trh2!=trhC2.end(); ++trh2, ++nTot1DHits_j) {
if( ( (*trh1)->globalPosition() - (*trh2)->globalPosition()).mag() < 10e-5 ) match++;
// if ( ( (*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag()< 10e-5 ) match++;
} // end for( trh2 ... )
} // end for( m2 ... )
} // end for( trh1 ... )
} // end for( m1 ... )
int nTotHits_i = (*iter)->measurements().size();
int nTotHits_j = (*jter)->measurements().size();
// FIXME Set Boff/on via cfg!
double chi2_dof_i = (*iter)->ndof() > 0 ? (*iter)->chiSquared()/(*iter)->ndof() : (*iter)->chiSquared()/1e-10;
double chi2_dof_j = (*jter)->ndof() > 0 ? (*jter)->chiSquared()/(*jter)->ndof() : (*jter)->chiSquared()/1e-10;
LogTrace(metname) << " MuonTrajectoryCleaner:";
LogTrace(metname) << " * trajC " << i
<< " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp()
<< " GeV) - chi2/nDOF = " << (*iter)->chiSquared() << "/" << (*iter)->ndof() << " = " << chi2_dof_i;
LogTrace(metname) << " - valid RH = " << (*iter)->foundHits() << " / total RH = " << nTotHits_i << " / total 1D RH = " << nTot1DHits_i;
LogTrace(metname) << " * trajC " << j
<< " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp()
<< " GeV) - chi2/nDOF = " << (*jter)->chiSquared() << "/" << (*jter)->ndof() << " = " << chi2_dof_j;
LogTrace(metname) << " - valid RH = " << (*jter)->foundHits() << " / total RH = " << nTotHits_j << " / total 1D RH = " << nTot1DHits_j;
LogTrace(metname) << " *** Shared 1D RecHits: " << match;
int hit_diff = (*iter)->foundHits() - (*jter)->foundHits() ;
// If there are matches, reject the worst track
if ( match > 0 ) {
// If the difference of # of rechits is less than 4, compare the chi2/ndf
if ( abs(hit_diff) <= 4 ) {
double minPt = 3.5;
double dPt = 7.; // i.e. considering 10% (conservative!) resolution at minPt it is ~ 10 sigma away from the central value
double maxFraction = 0.95;
double fraction = (2.*match)/((*iter)->measurements().size()+(*jter)->measurements().size());
int belowLimit = 0;
int above = 0;
if((*jter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt) ++belowLimit;
if((*iter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt) ++belowLimit;
if((*jter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt) ++above;
if((*iter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt) ++above;
if(fraction >= maxFraction && belowLimit == 1 && above == 1){
if((*iter)->lastMeasurement().updatedState().globalMomentum().perp() < minPt){
mask[i] = false;
skipnext=true;
seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
seedmap.erase(i);
LogTrace(metname) << "Trajectory # " << i << " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp()
<< " GeV) rejected because it has too low pt";
}
else {
mask[j] = false;
seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
seedmap.erase(j);
LogTrace(metname) << "Trajectory # " << j << " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp()
<< " GeV) rejected because it has too low pt";
}
}
else{
if (chi2_dof_i > chi2_dof_j) {
mask[i] = false;
skipnext=true;
seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
seedmap.erase(i);
LogTrace(metname) << "Trajectory # " << i << " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
}
else{
mask[j] = false;
seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
seedmap.erase(j);
LogTrace(metname) << "Trajectory # " << j << " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
}
}
}
else { // different number of hits
if ( hit_diff < 0 ) {
mask[i] = false;
skipnext=true;
seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
seedmap.erase(i);
LogTrace(metname) << "Trajectory # " << i << " (pT="<<(*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
}
else {
mask[j] = false;
seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
seedmap.erase(j);
LogTrace(metname) << "Trajectory # " << j << " (pT="<<(*jter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV) rejected";
}
}
}
if(skipnext) break;
j++;
}
i++;
if(skipnext) continue;
}
// Association map between the seed of the chosen trajectory and the seeds of the discarded trajectories
if(reportGhosts_) {
LogTrace(metname) << " Creating map between chosen seed and ghost seeds." << std::endl;
auto_ptr<L2SeedAssoc> seedToSeedsMap(new L2SeedAssoc);
int seedcnt(0);
for(map<int, vector<int> >::iterator itmap=seedmap.begin(); itmap!=seedmap.end(); ++itmap, ++seedcnt) {
edm::RefToBase<TrajectorySeed> tmpSeedRef1 = trajC[(*itmap).first]->seedRef();
edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef1 = tmpSeedRef1.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
int tmpsize=(*itmap).second.size();
for(int cnt=0; cnt!=tmpsize; ++cnt) {
edm::RefToBase<TrajectorySeed> tmpSeedRef2 = trajC[(*itmap).second[cnt]]->seedRef();
edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef2 = tmpSeedRef2.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
seedToSeedsMap->insert(tmpL2SeedRef1, tmpL2SeedRef2);
}
}
event.put(seedToSeedsMap, "");
} // end if(reportGhosts_)
i = 0;
for ( iter = trajC.begin(); iter != trajC.end(); iter++ ) {
if ( mask[i] ){
result.push_back(*iter);
LogTrace(metname) << "Keep trajectory with pT = " << (*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV";
}
else delete *iter;
i++;
}
trajC.clear();
trajC = result;
}
| void MuonTrajectoryCleaner::clean | ( | CandidateContainer & | muonTrajectories | ) |
Clean the candidates container, erasing the (worst) clone trajectory.
Definition at line 217 of file MuonTrajectoryCleaner.cc.
References alongMomentum, Geom::deltaPhi(), PV3DBase< T, PVType, FrameType >::eta(), TrajectoryStateOnSurface::globalMomentum(), i, TrajectoryStateOnSurface::isValid(), j, LogTrace, mag(), match(), metname, oppositeToMomentum, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), and query::result.
{
const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
if ( candC.size() < 2 ) return;
CandidateContainer::iterator iter, jter;
Trajectory::DataContainer::const_iterator m1, m2;
const float deltaEta = 0.01;
const float deltaPhi = 0.01;
const float deltaPt = 1.0;
LogTrace(metname) << "Number of muon candidates in the container: " <<candC.size()<< endl;
int i(0), j(0);
int match(0);
bool directionMatch = false;
// CAVEAT: vector<bool> is not a vector, its elements are not addressable!
// This is fine as long as only operator [] is used as in this case.
// cf. par 16.3.11
vector<bool> mask(candC.size(),true);
CandidateContainer result;
for ( iter = candC.begin(); iter != candC.end(); iter++ ) {
if ( !mask[i] ) { i++; continue; }
const Trajectory::DataContainer& meas1 = (*iter)->trajectory()->measurements();
j = i+1;
bool skipnext=false;
TrajectoryStateOnSurface innerTSOS;
if ((*iter)->trajectory()->direction() == alongMomentum) {
innerTSOS = (*iter)->trajectory()->firstMeasurement().updatedState();
}
else if ((*iter)->trajectory()->direction() == oppositeToMomentum) {
innerTSOS = (*iter)->trajectory()->lastMeasurement().updatedState();
}
if ( !(innerTSOS.isValid()) ) continue;
float pt1 = innerTSOS.globalMomentum().perp();
float eta1 = innerTSOS.globalMomentum().eta();
float phi1 = innerTSOS.globalMomentum().phi();
for ( jter = iter+1; jter != candC.end(); jter++ ) {
if ( !mask[j] ) { j++; continue; }
directionMatch = false;
const Trajectory::DataContainer& meas2 = (*jter)->trajectory()->measurements();
match = 0;
for ( m1 = meas1.begin(); m1 != meas1.end(); m1++ ) {
for ( m2 = meas2.begin(); m2 != meas2.end(); m2++ ) {
if ( (*m1).recHit()->isValid() && (*m2).recHit()->isValid() )
if ( ( (*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag()< 10e-5 ) match++;
}
}
LogTrace(metname)
<< " MuonTrajectoryCleaner: candC " << i << " chi2/nRH = "
<< (*iter)->trajectory()->chiSquared() << "/" << (*iter)->trajectory()->foundHits() <<
" vs trajC " << j << " chi2/nRH = " << (*jter)->trajectory()->chiSquared() <<
"/" << (*jter)->trajectory()->foundHits() << " Shared RecHits: " << match;
TrajectoryStateOnSurface innerTSOS2;
if ((*jter)->trajectory()->direction() == alongMomentum) {
innerTSOS2 = (*jter)->trajectory()->firstMeasurement().updatedState();
}
else if ((*jter)->trajectory()->direction() == oppositeToMomentum) {
innerTSOS2 = (*jter)->trajectory()->lastMeasurement().updatedState();
}
if ( !(innerTSOS2.isValid()) ) continue;
float pt2 = innerTSOS2.globalMomentum().perp();
float eta2 = innerTSOS2.globalMomentum().eta();
float phi2 = innerTSOS2.globalMomentum().phi();
float deta(fabs(eta1-eta2));
float dphi(fabs(Geom::Phi<float>(phi1)-Geom::Phi<float>(phi2)));
float dpt(fabs(pt1-pt2));
if ( dpt < deltaPt && deta < deltaEta && dphi < deltaPhi ) {
directionMatch = true;
LogTrace(metname)
<< " MuonTrajectoryCleaner: candC " << i<<" and "<<j<< " direction matched: "
<<innerTSOS.globalMomentum()<<" and " <<innerTSOS2.globalMomentum();
}
// If there are matches, reject the worst track
bool hitsMatch = ((match > 0) && ((match/((*iter)->trajectory()->foundHits()) > 0.25) || (match/((*jter)->trajectory()->foundHits()) > 0.25))) ? true : false;
bool tracksMatch =
( ( (*iter)->trackerTrack() == (*jter)->trackerTrack() ) &&
( deltaR<double>((*iter)->muonTrack()->eta(),(*iter)->muonTrack()->phi(), (*jter)->muonTrack()->eta(),(*jter)->muonTrack()->phi()) < 0.2) );
//if ( ( tracksMatch ) || (hitsMatch > 0) || directionMatch ) {
if ( ( tracksMatch ) || (hitsMatch > 0) ) {
if ( (*iter)->trajectory()->foundHits() == (*jter)->trajectory()->foundHits() ) {
if ( (*iter)->trajectory()->chiSquared() > (*jter)->trajectory()->chiSquared() ) {
mask[i] = false;
skipnext=true;
}
else mask[j] = false;
}
else { // different number of hits
if ( (*iter)->trajectory()->foundHits() < (*jter)->trajectory()->foundHits() ) {
mask[i] = false;
skipnext=true;
}
else mask[j] = false;
}
}
if(skipnext) break;
j++;
}
i++;
if(skipnext) continue;
}
i = 0;
for ( iter = candC.begin(); iter != candC.end(); iter++ ) {
if ( mask[i] ) {
result.push_back(*iter);
} else {
delete (*iter)->trajectory();
delete (*iter)->trackerTrajectory();
delete *iter;
}
i++;
}
candC.clear();
candC = result;
}
bool MuonTrajectoryCleaner::reportGhosts_ [private] |
Definition at line 44 of file MuonTrajectoryCleaner.h.
1.7.3