CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
MuonTrajectoryCleaner.cc
Go to the documentation of this file.
1 
6 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryCleaner.h"
7 #include "DataFormats/TrackReco/interface/Track.h"
8 #include "DataFormats/Math/interface/deltaR.h"
9 #include "FWCore/MessageLogger/interface/MessageLogger.h"
10 #include "DataFormats/Common/interface/RefToBase.h"
11 #include "DataFormats/Common/interface/AssociationMap.h"
12 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
13 #include "DataFormats/TrackingRecHit/interface/RecSegment.h"
14 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
15 
16 using namespace std;
17 
18 typedef edm::AssociationMap<edm::OneToMany<L2MuonTrajectorySeedCollection, L2MuonTrajectorySeedCollection> > L2SeedAssoc;
19 
20 void MuonTrajectoryCleaner::clean(TrajectoryContainer& trajC,
21  edm::Event& event,
22  const edm::Handle<edm::View<TrajectorySeed> >& seeds) {
23  const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
24 
25  LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
26 
27  TrajectoryContainer::iterator iter, jter;
28  Trajectory::DataContainer::const_iterator m1, m2;
29 
30  if (trajC.empty())
31  return;
32 
33  LogTrace(metname) << "Number of trajectories in the container: " << trajC.size() << endl;
34 
35  int i(0), j(0);
36  int match(0);
37 
38  // Map between chosen seed and ghost seeds (for trigger)
39  map<int, vector<int> > seedmap;
40 
41  // CAVEAT: vector<bool> is not a vector, its elements are not addressable!
42  // This is fine as long as only operator [] is used as in this case.
43  // cf. par 16.3.11
44  vector<bool> mask(trajC.size(), true);
45 
46  TrajectoryContainer result;
47 
48  for (iter = trajC.begin(); iter != trajC.end(); iter++) {
49  if (!mask[i]) {
50  i++;
51  continue;
52  }
53  if (!(*iter)->isValid() || (*iter)->empty()) {
54  mask[i] = false;
55  i++;
56  continue;
57  }
58  if (seedmap.count(i) == 0)
59  seedmap[i].push_back(i);
60  const Trajectory::DataContainer& meas1 = (*iter)->measurements();
61  j = i + 1;
62  bool skipnext = false;
63  for (jter = iter + 1; jter != trajC.end(); jter++) {
64  if (!mask[j]) {
65  j++;
66  continue;
67  }
68  if (seedmap.count(j) == 0)
69  seedmap[j].push_back(j);
70  const Trajectory::DataContainer& meas2 = (*jter)->measurements();
71  match = 0;
72  for (m1 = meas1.begin(); m1 != meas1.end(); m1++) {
73  for (m2 = meas2.begin(); m2 != meas2.end(); m2++) {
74  if (((*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag() < 10e-5)
75  match++;
76  } // end for( m2 ... )
77  } // end for( m1 ... )
78 
79  // FIXME Set Boff/on via cfg!
80  double chi2_dof_i = (*iter)->ndof() > 0 ? (*iter)->chiSquared() / (*iter)->ndof() : (*iter)->chiSquared() / 1e-10;
81  double chi2_dof_j = (*jter)->ndof() > 0 ? (*jter)->chiSquared() / (*jter)->ndof() : (*jter)->chiSquared() / 1e-10;
82 
83  LogTrace(metname) << " MuonTrajectoryCleaner:";
84  LogTrace(metname) << " * trajC " << i
85  << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
86  << " GeV) - chi2/nDOF = " << (*iter)->chiSquared() << "/" << (*iter)->ndof() << " = "
87  << chi2_dof_i;
88  LogTrace(metname) << " - valid RH = " << (*iter)->foundHits()
89  << " / total RH = " << (*iter)->measurements().size();
90  LogTrace(metname) << " * trajC " << j
91  << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
92  << " GeV) - chi2/nDOF = " << (*jter)->chiSquared() << "/" << (*jter)->ndof() << " = "
93  << chi2_dof_j;
94  LogTrace(metname) << " - valid RH = " << (*jter)->foundHits()
95  << " / total RH = " << (*jter)->measurements().size();
96  LogTrace(metname) << " *** Shared RecHits: " << match;
97 
98  int hit_diff = (*iter)->foundHits() - (*jter)->foundHits();
99  // If there are matches, reject the worst track
100  if (match > 0) {
101  // If the difference in # of rechits is null then compare the chi2/ndf of the two trajectories
102  if (abs(hit_diff) == 0) {
103  double minPt = 3.5;
104  double dPt =
105  7.; // i.e. considering 10% (conservative!) resolution at minPt it is ~ 10 sigma away from the central value
106 
107  double maxFraction = 0.95;
108 
109  double fraction = (2. * match) / ((*iter)->measurements().size() + (*jter)->measurements().size());
110  int belowLimit = 0;
111  int above = 0;
112 
113  if ((*jter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt)
114  ++belowLimit;
115  if ((*iter)->lastMeasurement().updatedState().globalMomentum().perp() <= minPt)
116  ++belowLimit;
117 
118  if ((*jter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt)
119  ++above;
120  if ((*iter)->lastMeasurement().updatedState().globalMomentum().perp() >= dPt)
121  ++above;
122 
123  if (fraction >= maxFraction && belowLimit == 1 && above == 1) {
124  if ((*iter)->lastMeasurement().updatedState().globalMomentum().perp() < minPt) {
125  mask[i] = false;
126  skipnext = true;
127  seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
128  seedmap.erase(i);
129  LogTrace(metname) << "Trajectory # " << i
130  << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
131  << " GeV) rejected because it has too low pt";
132  } else {
133  mask[j] = false;
134  seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
135  seedmap.erase(j);
136  LogTrace(metname) << "Trajectory # " << j
137  << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
138  << " GeV) rejected because it has too low pt";
139  }
140  } else {
141  if (chi2_dof_i > chi2_dof_j) {
142  mask[i] = false;
143  skipnext = true;
144  seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
145  seedmap.erase(i);
146  LogTrace(metname) << "Trajectory # " << i
147  << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
148  << " GeV) rejected";
149  } else {
150  mask[j] = false;
151  seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
152  seedmap.erase(j);
153  LogTrace(metname) << "Trajectory # " << j
154  << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
155  << " GeV) rejected";
156  }
157  }
158  } else { // different number of hits
159  if (hit_diff < 0) {
160  mask[i] = false;
161  skipnext = true;
162  seedmap[j].insert(seedmap[j].end(), seedmap[i].begin(), seedmap[i].end());
163  seedmap.erase(i);
164  LogTrace(metname) << "Trajectory # " << i
165  << " (pT=" << (*iter)->lastMeasurement().updatedState().globalMomentum().perp()
166  << " GeV) rejected";
167  } else {
168  mask[j] = false;
169  seedmap[i].insert(seedmap[i].end(), seedmap[j].begin(), seedmap[j].end());
170  seedmap.erase(j);
171  LogTrace(metname) << "Trajectory # " << j
172  << " (pT=" << (*jter)->lastMeasurement().updatedState().globalMomentum().perp()
173  << " GeV) rejected";
174  }
175  }
176  }
177  if (skipnext)
178  break;
179  j++;
180  }
181  i++;
182  if (skipnext)
183  continue;
184  }
185 
186  // Association map between the seed of the chosen trajectory and the seeds of the discarded trajectories
187  if (reportGhosts_) {
188  LogTrace(metname) << " Creating map between chosen seed and ghost seeds." << std::endl;
189 
190  auto seedToSeedsMap = std::make_unique<L2SeedAssoc>();
191  if (!seeds->empty()) {
192  edm::Ptr<TrajectorySeed> ptr = seeds->ptrAt(0);
193  edm::Handle<L2MuonTrajectorySeedCollection> seedsHandle;
194  event.get(ptr.id(), seedsHandle);
195  seedToSeedsMap = std::make_unique<L2SeedAssoc>(seedsHandle, seedsHandle);
196  }
197 
198  int seedcnt(0);
199 
200  for (map<int, vector<int> >::iterator itmap = seedmap.begin(); itmap != seedmap.end(); ++itmap, ++seedcnt) {
201  edm::RefToBase<TrajectorySeed> tmpSeedRef1 = trajC[(*itmap).first]->seedRef();
202  edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef1 =
203  tmpSeedRef1.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
204 
205  int tmpsize = (*itmap).second.size();
206  for (int cnt = 0; cnt != tmpsize; ++cnt) {
207  edm::RefToBase<TrajectorySeed> tmpSeedRef2 = trajC[(*itmap).second[cnt]]->seedRef();
208  edm::Ref<L2MuonTrajectorySeedCollection> tmpL2SeedRef2 =
209  tmpSeedRef2.castTo<edm::Ref<L2MuonTrajectorySeedCollection> >();
210  seedToSeedsMap->insert(tmpL2SeedRef1, tmpL2SeedRef2);
211  }
212  }
213 
214  event.put(std::move(seedToSeedsMap), "");
215 
216  } // end if(reportGhosts_)
217 
218  i = 0;
219  auto c = std::count(mask.begin(), mask.end(), true);
220  result.reserve(c);
221  for (iter = trajC.begin(); iter != trajC.end(); iter++) {
222  if (mask[i]) {
223  LogTrace(metname) << "Keep trajectory with pT = "
224  << (*iter)->lastMeasurement().updatedState().globalMomentum().perp() << " GeV";
225  result.push_back(std::move(*iter));
226  }
227  i++;
228  }
229 
230  trajC = std::move(result);
231 }
232 
233 //
234 // clean CandidateContainer
235 //
236 void MuonTrajectoryCleaner::clean(CandidateContainer& candC) {
237  const std::string metname = "Muon|RecoMuon|MuonTrajectoryCleaner";
238 
239  LogTrace(metname) << "Muon Trajectory Cleaner called" << endl;
240 
241  if (candC.size() < 2)
242  return;
243 
244  CandidateContainer::iterator iter, jter;
245  Trajectory::DataContainer::const_iterator m1, m2;
246 
247  const float deltaEta = 0.01;
248  const float deltaPhi = 0.01;
249  const float deltaPt = 1.0;
250 
251  LogTrace(metname) << "Number of muon candidates in the container: " << candC.size() << endl;
252 
253  int i(0), j(0);
254  int match(0);
255  //unused bool directionMatch = false;
256 
257  // CAVEAT: vector<bool> is not a vector, its elements are not addressable!
258  // This is fine as long as only operator [] is used as in this case.
259  // cf. par 16.3.11
260  vector<bool> mask(candC.size(), true);
261 
262  CandidateContainer result;
263 
264  for (iter = candC.begin(); iter != candC.end(); iter++) {
265  if (!mask[i]) {
266  i++;
267  continue;
268  }
269  const Trajectory::DataContainer& meas1 = (*iter)->trajectory()->measurements();
270  j = i + 1;
271  bool skipnext = false;
272 
273  TrajectoryStateOnSurface innerTSOS;
274 
275  if ((*iter)->trajectory()->direction() == alongMomentum) {
276  innerTSOS = (*iter)->trajectory()->firstMeasurement().updatedState();
277  } else if ((*iter)->trajectory()->direction() == oppositeToMomentum) {
278  innerTSOS = (*iter)->trajectory()->lastMeasurement().updatedState();
279  }
280  if (!(innerTSOS.isValid()))
281  continue;
282  float pt1 = innerTSOS.globalMomentum().perp();
283  float eta1 = innerTSOS.globalMomentum().eta();
284  float phi1 = innerTSOS.globalMomentum().phi();
285 
286  for (jter = iter + 1; jter != candC.end(); jter++) {
287  if (!mask[j]) {
288  j++;
289  continue;
290  }
291  //UNUSED: directionMatch = false;
292  const Trajectory::DataContainer& meas2 = (*jter)->trajectory()->measurements();
293  match = 0;
294  for (m1 = meas1.begin(); m1 != meas1.end(); m1++) {
295  for (m2 = meas2.begin(); m2 != meas2.end(); m2++) {
296  if ((*m1).recHit()->isValid() && (*m2).recHit()->isValid())
297  if (((*m1).recHit()->globalPosition() - (*m2).recHit()->globalPosition()).mag() < 10e-5)
298  match++;
299  }
300  }
301 
302  LogTrace(metname) << " MuonTrajectoryCleaner: candC " << i
303  << " chi2/nRH = " << (*iter)->trajectory()->chiSquared() << "/"
304  << (*iter)->trajectory()->foundHits() << " vs trajC " << j
305  << " chi2/nRH = " << (*jter)->trajectory()->chiSquared() << "/"
306  << (*jter)->trajectory()->foundHits() << " Shared RecHits: " << match;
307 
308  TrajectoryStateOnSurface innerTSOS2;
309  if ((*jter)->trajectory()->direction() == alongMomentum) {
310  innerTSOS2 = (*jter)->trajectory()->firstMeasurement().updatedState();
311  } else if ((*jter)->trajectory()->direction() == oppositeToMomentum) {
312  innerTSOS2 = (*jter)->trajectory()->lastMeasurement().updatedState();
313  }
314  if (!(innerTSOS2.isValid()))
315  continue;
316 
317  float pt2 = innerTSOS2.globalMomentum().perp();
318  float eta2 = innerTSOS2.globalMomentum().eta();
319  float phi2 = innerTSOS2.globalMomentum().phi();
320 
321  float deta(fabs(eta1 - eta2));
322  float dphi(fabs(Geom::Phi<float>(phi1) - Geom::Phi<float>(phi2)));
323  float dpt(fabs(pt1 - pt2));
324  if (dpt < deltaPt && deta < deltaEta && dphi < deltaPhi) {
325  //UNUSED: directionMatch = true;
326  LogTrace(metname) << " MuonTrajectoryCleaner: candC " << i << " and " << j
327  << " direction matched: " << innerTSOS.globalMomentum() << " and "
328  << innerTSOS2.globalMomentum();
329  }
330 
331  // If there are matches, reject the worst track
332  bool hitsMatch = ((match > 0) && ((match / ((*iter)->trajectory()->foundHits()) > 0.25) ||
333  (match / ((*jter)->trajectory()->foundHits()) > 0.25)))
334  ? true
335  : false;
336  bool tracksMatch =
337  (((*iter)->trackerTrack() == (*jter)->trackerTrack()) && (deltaR<double>((*iter)->muonTrack()->eta(),
338  (*iter)->muonTrack()->phi(),
339  (*jter)->muonTrack()->eta(),
340  (*jter)->muonTrack()->phi()) < 0.2));
341 
342  //if ( ( tracksMatch ) || (hitsMatch > 0) || directionMatch ) {
343  if ((tracksMatch) || (hitsMatch > 0)) {
344  if ((*iter)->trajectory()->foundHits() == (*jter)->trajectory()->foundHits()) {
345  if ((*iter)->trajectory()->chiSquared() > (*jter)->trajectory()->chiSquared()) {
346  mask[i] = false;
347  skipnext = true;
348  } else
349  mask[j] = false;
350  } else { // different number of hits
351  if ((*iter)->trajectory()->foundHits() < (*jter)->trajectory()->foundHits()) {
352  mask[i] = false;
353  skipnext = true;
354  } else
355  mask[j] = false;
356  }
357  }
358  if (skipnext)
359  break;
360  j++;
361  }
362  i++;
363  if (skipnext)
364  continue;
365  }
366 
367  auto c = std::count(mask.begin(), mask.end(), true);
368  i = 0;
369  result.reserve(c);
370  for (iter = candC.begin(); iter != candC.end(); iter++) {
371  if (mask[i]) {
372  result.push_back(std::move(*iter));
373  }
374  i++;
375  }
376 
377  candC = std::move(result);
378 }
c
const edm::EventSetup & c
Definition: SiStripLAProfileBooker.cc:66
MuonTrajectoryCleaner::clean
void clean(TrajectoryContainer &muonTrajectories, edm::Event &evt, const edm::Handle< edm::View< TrajectorySeed > > &seeds)
Clean the trajectories container, erasing the (worst) clone trajectory.
Definition: MuonTrajectoryCleaner.cc:20
mps_fire.i
i
Definition: mps_fire.py:428
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
metname
const std::string metname
Definition: MuonSeedOrcaPatternRecognition.cc:40
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
HLT_FULL_cff.eta2
tuple eta2
Definition: HLT_FULL_cff.py:9780
HLT_FULL_cff.fraction
tuple fraction
Definition: HLT_FULL_cff.py:53121
srCondWrite_cfg.deltaPhi
tuple deltaPhi
Definition: srCondWrite_cfg.py:54
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
spr::deltaEta
static const double deltaEta
Definition: CaloConstants.h:8
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
callgraph.m2
tuple m2
Definition: callgraph.py:57
L2SeedAssoc
edm::AssociationMap< edm::OneToMany< L2MuonTrajectorySeedCollection, L2MuonTrajectorySeedCollection > > L2SeedAssoc
Definition: MuonTrajectoryCleaner.cc:18
mps_fire.result
tuple result
Definition: mps_fire.py:311
Trajectory::DataContainer
std::vector< TrajectoryMeasurement > DataContainer
Definition: Trajectory.h:40
MuonTrajectoryCleaner::CandidateContainer
MuonCandidate::CandidateContainer CandidateContainer
Definition: MuonTrajectoryCleaner.h:19
eostools.move
def move
Definition: eostools.py:511
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
DetachedQuadStep_cff.seeds
string seeds
Definition: DetachedQuadStep_cff.py:212
HLT_FULL_cff.minPt
tuple minPt
Definition: HLT_FULL_cff.py:99
MuonTrajectoryCleaner::TrajectoryContainer
MuonCandidate::TrajectoryContainer TrajectoryContainer
Definition: MuonTrajectoryCleaner.h:18
HLT_FULL_cff.maxFraction
tuple maxFraction
Definition: HLT_FULL_cff.py:52278
edm::RefToBase::castTo
REF castTo() const
Definition: RefToBase.h:257
edm::Ptr::id
ProductID id() const
Accessor for product ID.
Definition: Ptr.h:158
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:73
TrajectoryStateOnSurface::globalMomentum
GlobalVector globalMomentum() const
Definition: TrajectoryStateOnSurface.h:66
dataset.end
string end
Definition: dataset.py:937
match
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
HLT_FULL_cff.eta1
tuple eta1
Definition: HLT_FULL_cff.py:9779
Generated for CMSSW Reference Manual by doxygen 1.8.5