CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CosmicMuonTrajectoryBuilder.cc
Go to the documentation of this file.
13 
14 /* Collaborating Class Header */
38 
39 #include <algorithm>
40 
41 using namespace edm;
42 using namespace std;
43 
45 
46  thePropagatorName = par.getParameter<string>("Propagator");
47 
48  bool enableDTMeasurement = par.getParameter<bool>("EnableDTMeasurement");
49  bool enableCSCMeasurement = par.getParameter<bool>("EnableCSCMeasurement");
50  bool enableRPCMeasurement = par.getParameter<bool>("EnableRPCMeasurement");
51 
52 
53 // if(enableDTMeasurement)
54  InputTag DTRecSegmentLabel = par.getParameter<InputTag>("DTRecSegmentLabel");
55 
56 // if(enableCSCMeasurement)
57  InputTag CSCRecSegmentLabel = par.getParameter<InputTag>("CSCRecSegmentLabel");
58 
59 // if(enableRPCMeasurement)
60  InputTag RPCRecSegmentLabel = par.getParameter<InputTag>("RPCRecSegmentLabel");
61 
62 // if(enableGEMMeasurement)
63 // InputTag GEMRecSegmentLabel = par.getParameter<InputTag>("GEMRecSegmentLabel");
64 
65  theLayerMeasurements= new MuonDetLayerMeasurements(DTRecSegmentLabel,
66  CSCRecSegmentLabel,
67  RPCRecSegmentLabel,
68  edm::InputTag(),
69  edm::InputTag(),
70  iC,
71  enableDTMeasurement,
72  enableCSCMeasurement,
73  enableRPCMeasurement,
74  false,
75  false);
76 
77  ParameterSet muonUpdatorPSet = par.getParameter<ParameterSet>("MuonTrajectoryUpdatorParameters");
78 
79  theNavigation = 0; // new DirectMuonNavigation(theService->detLayerGeometry());
80  theUpdator = new MuonTrajectoryUpdator(muonUpdatorPSet, insideOut);
81 
83 
84  ParameterSet muonBackwardUpdatorPSet = par.getParameter<ParameterSet>("BackwardMuonTrajectoryUpdatorParameters");
85 
86  theBKUpdator = new MuonTrajectoryUpdator(muonBackwardUpdatorPSet, outsideIn);
87 
88  theTraversingMuonFlag = par.getParameter<bool>("BuildTraversingMuon");
89 
90  theStrict1LegFlag = par.getParameter<bool>("Strict1Leg");
91 
92  ParameterSet smootherPSet = par.getParameter<ParameterSet>("MuonSmootherParameters");
93 
94  theNavigationPSet = par.getParameter<ParameterSet>("MuonNavigationParameters");
95 
96  theSmoother = new CosmicMuonSmoother(smootherPSet,theService);
97 
98  theNTraversing = 0;
99  theNSuccess = 0;
100  theCacheId_DG = 0;
101  category_ = "Muon|RecoMuon|CosmicMuon|CosmicMuonTrajectoryBuilder";
102 
103 }
104 
105 
107 
108  LogTrace(category_)<< "CosmicMuonTrajectoryBuilder dtor called";
109  if (theUpdator) delete theUpdator;
110  if (theBKUpdator) delete theBKUpdator;
112  if (theSmoother) delete theSmoother;
113  if (theNavigation) delete theNavigation;
115 
116  LogTrace(category_)<< "CosmicMuonTrajectoryBuilder Traversing: "<<theNSuccess<<"/"<<theNTraversing;
117 
118 }
119 
120 
122 
124 
125  // DetLayer Geometry
126  unsigned long long newCacheId_DG = theService->eventSetup().get<MuonRecoGeometryRecord>().cacheIdentifier();
127  if ( newCacheId_DG != theCacheId_DG ) {
128  LogTrace(category_) << "Muon Reco Geometry changed!";
129  theCacheId_DG = newCacheId_DG;
130  if (theNavigation) delete theNavigation;
132  }
133 
134 
135 }
136 
137 
140 
141  vector<Trajectory*> trajL = vector<Trajectory*>();
142 
144 
145  PTrajectoryStateOnDet ptsd1(seed.startingState());
146  DetId did(ptsd1.detId());
147  const BoundPlane& bp = theService->trackingGeometry()->idToDet(did)->surface();
148  TrajectoryStateOnSurface lastTsos = trajectoryStateTransform::transientState(ptsd1,&bp,&*theService->magneticField());
149  LogTrace(category_) << "Seed: mom "<<lastTsos.globalMomentum()
150  <<"pos: " <<lastTsos.globalPosition();
151  LogTrace(category_) << "Seed: mom eta "<<lastTsos.globalMomentum().eta()
152  <<"pos eta: " <<lastTsos.globalPosition().eta();
153 
154  bool beamhaloFlag = ( (did.subdetId() == MuonSubdetId::CSC) && fabs(lastTsos.globalMomentum().eta()) > 4.0);
155 
156  vector<const DetLayer*> navLayers;
157 
158  if (did.subdetId() == MuonSubdetId::DT) {
159  //DT
160  navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
161  }
162  else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
163  //CSC
164  navLayers = navigation()->compatibleEndcapLayers(*(lastTsos.freeState()), alongMomentum);
165  } else {
166  navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
167  }
168 
169  LogTrace(category_) <<"found "<<navLayers.size()<<" compatible DetLayers for the Seed";
170 
171  if (navLayers.empty()) return trajL;
172 
173  vector<DetWithState> detsWithStates;
174  LogTrace(category_) << "Compatible layers: ";
175  for ( vector<const DetLayer*>::const_iterator layer = navLayers.begin();
176  layer != navLayers.end(); layer++) {
177  LogTrace(category_) << debug.dumpMuonId((*layer)->basicComponents().front()->geographicalId())
178  << debug.dumpLayer(*layer);
179  }
180 
181  detsWithStates = navLayers.front()->compatibleDets(lastTsos, *propagator(), *(updator()->estimator()));
182  LogTrace(category_) << "Number of compatible dets: " << detsWithStates.size() << endl;
183 
184  if ( !detsWithStates.empty() ) {
185  // get the updated TSOS
186  if ( detsWithStates.front().second.isValid() ) {
187  LogTrace(category_) << "New starting TSOS is on det: " << endl;
188  LogTrace(category_) << debug.dumpMuonId(detsWithStates.front().first->geographicalId())
189  << debug.dumpLayer(navLayers.front());
190  lastTsos = detsWithStates.front().second;
191  LogTrace(category_) << "Seed after extrapolation: mom " << lastTsos.globalMomentum()
192  << "pos: " << lastTsos.globalPosition();
193  }
194  }
195  detsWithStates.clear();
196  if ( !lastTsos.isValid() ) return trajL;
197 
198  TrajectoryStateOnSurface secondLast = lastTsos;
199 
200  lastTsos.rescaleError(10.0);
201 
202  Trajectory* theTraj = new Trajectory(seed,alongMomentum);
203 
204  navLayers.clear();
205 
206  if (fabs(lastTsos.globalMomentum().eta()) < 1.0) {
207  //DT
208  navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
209  } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
210  //CSC
211  navLayers = navigation()->compatibleEndcapLayers(*(lastTsos.freeState()), alongMomentum);
212  } else {
213  navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
214  }
215 
216 
217  int DTChamberUsedBack = 0;
218  int CSCChamberUsedBack = 0;
219  int RPCChamberUsedBack = 0;
220  int TotalChamberUsedBack = 0;
222  vector<TrajectoryMeasurement> measL;
223 
224  LogTrace(category_) << "Begin forward fit " << navLayers.size();
225 
226  for ( vector<const DetLayer*>::const_iterator rnxtlayer = navLayers.begin(); rnxtlayer!= navLayers.end(); ++rnxtlayer) {
227  LogTrace(category_) << "new layer ";
228  measL.clear();
229  LogTrace(category_) << debug.dumpMuonId((*rnxtlayer)->basicComponents().front()->geographicalId())
230  << debug.dumpLayer(*rnxtlayer);
231  LogTrace(category_) << "from lastTsos " << lastTsos.globalMomentum()<< " at " <<lastTsos.globalPosition();
232 
233  measL = findBestMeasurements(*rnxtlayer, lastTsos, propagator(), (updator()->estimator()));
234 
235  if ( measL.empty() && (fabs(theService->magneticField()->inTesla(GlobalPoint(0,0,0)).z()) < 0.01) && (theService->propagator("StraightLinePropagator").isValid() ) ) {
236  LogTrace(category_) << "try straight line propagator ";
237  measL = findBestMeasurements(*rnxtlayer, lastTsos, &*theService->propagator("StraightLinePropagator"), (updator()->estimator()));
238  }
239  if ( measL.empty() ) continue;
240 
241  for (vector<TrajectoryMeasurement>::const_iterator theMeas = measL.begin(); theMeas != measL.end(); ++theMeas) {
242  pair<bool,TrajectoryStateOnSurface> result = updator()->update((&*theMeas), *theTraj, propagator());
243 
244  if (result.first ) {
245  LogTrace(category_) << "update ok ";
246  incrementChamberCounters((*rnxtlayer), DTChamberUsedBack, CSCChamberUsedBack, RPCChamberUsedBack, TotalChamberUsedBack);
247  secondLast = lastTsos;
248  if ( (!theTraj->empty()) && result.second.isValid() ) {
249  lastTsos = result.second;
250  LogTrace(category_) << "get new lastTsos here " << lastTsos.globalMomentum() << " at " << lastTsos.globalPosition();
251  } else if ((theMeas)->predictedState().isValid()) lastTsos = (theMeas)->predictedState();
252  }
253  }
254  }
255  measL.clear();
256  while (!theTraj->empty()) {
257  theTraj->pop();
258  }
259 
260  if (!theTraj->isValid() || TotalChamberUsedBack < 2 || (DTChamberUsedBack+CSCChamberUsedBack) == 0 || !lastTsos.isValid()) {
261  delete theTraj;
262  return trajL;
263  }
264  delete theTraj;
265 
266 
267  // if got good trajectory, then do backward refitting
268  DTChamberUsedBack = 0;
269  CSCChamberUsedBack = 0;
270  RPCChamberUsedBack = 0;
271  TotalChamberUsedBack = 0;
272 
273  Trajectory myTraj(seed, oppositeToMomentum);
274 
275  // set starting navigation direction for MuonTrajectoryUpdator
276 
277  GlobalPoint lastPos = lastTsos.globalPosition();
278  GlobalPoint secondLastPos = secondLast.globalPosition();
279  GlobalVector momDir = secondLastPos - lastPos;
280 
281  if ( lastPos.basicVector().dot(momDir.basicVector()) > 0 ) {
282 // LogTrace("CosmicMuonTrajectoryBuilder")<<"Fit direction changed to insideOut";
285 
286  if (fabs(lastTsos.globalMomentum().eta()) < 1.0) {
287  //DT
288  navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), oppositeToMomentum);
289  } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
290  //CSC
291  std::reverse(navLayers.begin(), navLayers.end());
292  } else {
293  navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), oppositeToMomentum);
294  }
295 
296  LogTrace(category_) << "Begin backward refitting, with " << navLayers.size() << " layers" << endl;
297 
298  for (vector<const DetLayer*>::const_iterator rnxtlayer = navLayers.begin();
299  rnxtlayer!= navLayers.end(); ++rnxtlayer) {
300 
301  measL.clear();
302 
303  measL = findBestMeasurements(*rnxtlayer, lastTsos, propagator(), (backwardUpdator()->estimator()));
304 
305  if ( measL.empty() ) {
307  for (MuonRecHitContainer::const_iterator ihit = tmpHits.begin();
308  ihit != tmpHits.end(); ++ihit ) {
309  allUnusedHits.push_back(*ihit);
310  }
311  continue;
312  }
313 
314  for (vector<TrajectoryMeasurement>::const_iterator theMeas = measL.begin(); theMeas != measL.end(); ++theMeas) {
315 
316  // if the part change, we need to reconsider the fit direction
317  if (rnxtlayer != navLayers.begin()) {
318 
319  vector<const DetLayer*>::const_iterator lastlayer = rnxtlayer;
320  lastlayer--;
321 
322  if ( (*rnxtlayer)->location() != (*lastlayer)->location() ) {
323 
324  lastPos = lastTsos.globalPosition();
325  GlobalPoint thisPos = (theMeas)->predictedState().globalPosition();
326  GlobalVector momDir = thisPos - lastPos;
327 // LogTrace("CosmicMuonTrajectoryBuilder")<<"momDir "<<momDir;
328 
329  if ( momDir.mag() > 0.01 ) { //if lastTsos is on the surface, no need
330  if ( thisPos.basicVector().dot(momDir.basicVector()) > 0 ) {
333  }
334  }
335  if ( ((*lastlayer)->location() == GeomDetEnumerators::endcap) &&
336  ((*rnxtlayer)->location() == GeomDetEnumerators::endcap) &&
337  (lastTsos.globalPosition().z() * (theMeas)->predictedState().globalPosition().z() < 0) ) {
339  }
340  }
341 
342 // if (theBKUpdator->fitDirection() == insideOut)
343 // LogTrace("CosmicMuonTrajectoryBuilder")<<"Fit direction insideOut";
344 // else LogTrace("CosmicMuonTrajectoryBuilder")<<"Fit direction outsideIn";
345  pair<bool,TrajectoryStateOnSurface> bkresult
346  = backwardUpdator()->update((&*theMeas), myTraj, propagator());
347 
348  if (bkresult.first ) {
349 
350  incrementChamberCounters((*rnxtlayer), DTChamberUsedBack, CSCChamberUsedBack, RPCChamberUsedBack, TotalChamberUsedBack);
351 
352  if ( theTraversingMuonFlag ) {
353 
354  MuonRecHitContainer tmpUnusedHits = unusedHits(*rnxtlayer,*theMeas);
355  allUnusedHits.insert(allUnusedHits.end(),tmpUnusedHits.begin(),tmpUnusedHits.end());
356  }
357  if ( (!myTraj.empty()) && bkresult.second.isValid() )
358  lastTsos = bkresult.second;
359  else if ((theMeas)->predictedState().isValid())
360  lastTsos = (theMeas)->predictedState();
361  }
362  }
363  }
364 
365  for ( vector<Trajectory*>::iterator t = trajL.begin(); t != trajL.end(); ++t ) delete *t;
366 
367  trajL.clear();
368 
369  if (( !myTraj.isValid() ) || ( myTraj.empty() ) || ( (selfDuplicate(myTraj)) )|| TotalChamberUsedBack < 2 || (DTChamberUsedBack+CSCChamberUsedBack) < 1) {
370  return trajL;
371  }
372 
373  if ( theTraversingMuonFlag && ( allUnusedHits.size() >= 2 )) {
374 // LogTrace(category_)<<utilities()->print(allUnusedHits);
375  LogTrace(category_)<<"Building trajectory in second hemisphere...";
376  buildSecondHalf(myTraj);
377  // check if traversing trajectory has hits in both hemispheres
378 
379  if ( theStrict1LegFlag && !utilities()->isTraversing(myTraj) ) {trajL.clear(); return trajL;}
380  } else if (theStrict1LegFlag && theTraversingMuonFlag) {trajL.clear(); return trajL;}
381 
382  LogTrace(category_) <<" traj ok ";
383 
384 // getDirectionByTime(myTraj);
385  if (beamhaloFlag) estimateDirection(myTraj);
386  if ( myTraj.empty() ) return trajL;
387 
388  // try to smooth it
389  vector<Trajectory> smoothed = theSmoother->trajectories(myTraj);
390 
391  if ( !smoothed.empty() && smoothed.front().foundHits()> 3 ) {
392  LogTrace(category_) <<" Smoothed successfully.";
393  myTraj = smoothed.front();
394  }
395  else {
396  LogTrace(category_) <<" Smooth failed.";
397  }
398 
399  LogTrace(category_) <<"first "<< myTraj.firstMeasurement().updatedState()
400  <<"\n last "<<myTraj.lastMeasurement().updatedState();
401  if ( myTraj.direction() == alongMomentum ) LogTrace(category_)<<"alongMomentum";
402  else if (myTraj.direction() == oppositeToMomentum ) LogTrace(category_)<<"oppositeMomentum";
403  else LogTrace(category_)<<"anyDirection";
404 
405  if (!beamhaloFlag) {
406  if ( myTraj.lastMeasurement().updatedState().globalMomentum().y() > 0 ) {
407  LogTrace(category_)<<"flip trajectory ";
408  flipTrajectory(myTraj);
409  }
410 
411  if ( ( myTraj.direction() == alongMomentum &&
413  < myTraj.lastMeasurement().updatedState().globalPosition().y()))
414  || (myTraj.direction() == oppositeToMomentum &&
416  > myTraj.lastMeasurement().updatedState().globalPosition().y())) ) {
417  LogTrace(category_)<<"reverse propagation direction";
419  }
420  }
421 // getDirectionByTime(myTraj);
422  if ( !myTraj.isValid() ) return trajL;
423 
424  // check direction agree with position!
426  GlobalVector dirFromPos = myTraj.measurements().back().recHit()->globalPosition() - myTraj.measurements().front().recHit()->globalPosition();
427 
428  if ( theStrict1LegFlag && !utilities()->isTraversing(myTraj) ) {trajL.clear(); return trajL;}
429 
430  LogTrace(category_)<< "last hit " <<myTraj.measurements().back().recHit()->globalPosition()<<endl;
431  LogTrace(category_)<< "first hit " <<myTraj.measurements().front().recHit()->globalPosition()<<endl;
432 
433  LogTrace(category_)<< "last tsos " <<myTraj.measurements().back().updatedState().globalPosition()<<" mom "<<myTraj.measurements().back().updatedState().globalMomentum()<<endl;
434  LogTrace(category_)<< "first tsos " <<myTraj.measurements().front().updatedState().globalPosition()<<" mom "<<myTraj.measurements().front().updatedState().globalMomentum()<<endl;
435 
436  PropagationDirection propDir =
438  LogTrace(category_)<<" dir "<<dir <<" propDir "<<propDir<<endl;
439 
440  LogTrace(category_)<<"chi2 "<<myTraj.chiSquared() <<endl;
441 
442  if (dir != propDir ) {
443  LogTrace(category_)<< "reverse propagation direction ";
445  }
446  if ( myTraj.empty() ) return trajL;
447 
448  trajL.push_back(new Trajectory(myTraj));
449  navLayers.clear();
450  return trajL;
451 }
452 
453 
454 //
455 //
456 //
459 
462  for (MuonRecHitContainer::const_iterator ihit = tmpHits.begin();
463  ihit != tmpHits.end(); ++ihit ) {
464  if ((*ihit)->geographicalId() != meas.recHit()->geographicalId() ){
465  result.push_back(*ihit);
466  LogTrace(category_) << "Unused hit: " << (*ihit)->globalPosition() << endl;
467  }
468  }
469 
470  return result;
471 
472 }
473 
474 
475 //
476 // continue to build a trajectory starting from given trajectory state
477 //
479  const NavigationDirection& startingDir,
480  Trajectory& traj) {
481 
482  if ( !ts.isValid() ) return;
483 
484  const FreeTrajectoryState* fts = ts.freeState();
485  if ( !fts ) return;
486 
487  vector<const DetLayer*> navLayers;
488 
489  if (fabs(fts->momentum().basicVector().eta()) < 1.0) {
490  //DT
491  if (fts->position().basicVector().dot(fts->momentum().basicVector())>0){
492  navLayers = navigation()->compatibleLayers((*fts), alongMomentum);
493  } else {
494  navLayers = navigation()->compatibleLayers((*fts), oppositeToMomentum);
495  }
496 
498  //CSC
499  if (fts->position().basicVector().dot(fts->momentum().basicVector())>0){
500  navLayers = navigation()->compatibleEndcapLayers((*fts), alongMomentum);
501  } else {
502  navLayers = navigation()->compatibleEndcapLayers((*fts), oppositeToMomentum);
503  }
504  } else {
505 
506  if (fts->position().basicVector().dot(fts->momentum().basicVector())>0){
507  navLayers = navigation()->compatibleLayers((*fts), alongMomentum);
508  } else {
509  navLayers = navigation()->compatibleLayers((*fts), oppositeToMomentum);
510  }
511 
512 }
513 
514  if (navLayers.empty()) return;
515 
516  theBKUpdator->setFitDirection(startingDir);
517 
518  int DTChamberUsedBack = 0;
519  int CSCChamberUsedBack = 0;
520  int RPCChamberUsedBack = 0;
521  int TotalChamberUsedBack = 0;
522 
523  TrajectoryStateOnSurface lastTsos =
526  propagatorAlong()->propagate((*fts),navLayers.front()->surface()) : propagatorOpposite()->propagate((*fts),navLayers.front()->surface());
527 
528  if ( !lastTsos.isValid() ) {
529  LogTrace(category_)<<"propagation failed from fts to inner cylinder";
530  return;
531  }
532  LogTrace(category_)<<"tsos "<<lastTsos.globalPosition();
533  lastTsos.rescaleError(10.);
534  vector<TrajectoryMeasurement> measL;
535  for (vector<const DetLayer*>::const_iterator rnxtlayer = navLayers.begin();
536  rnxtlayer!= navLayers.end(); ++rnxtlayer) {
537 
538  measL.clear();
539  measL = findBestMeasurements(*rnxtlayer, lastTsos, propagator(), (backwardUpdator()->estimator()));
540 
541  if ( measL.empty() ) continue;
542 
543  for (vector<TrajectoryMeasurement>::const_iterator theMeas = measL.begin(); theMeas != measL.end(); ++theMeas) {
544 
545  pair<bool,TrajectoryStateOnSurface> bkresult
546  = backwardUpdator()->update((&*theMeas), traj, propagator());
547  if (bkresult.first ) {
548  LogTrace(category_)<<"update ok : "<<(theMeas)->recHit()->globalPosition() ;
549 
550  incrementChamberCounters((*rnxtlayer), DTChamberUsedBack, CSCChamberUsedBack, RPCChamberUsedBack, TotalChamberUsedBack);
551 
552  if ( (!traj.empty()) && bkresult.second.isValid() )
553  lastTsos = bkresult.second;
554  else if ((theMeas)->predictedState().isValid())
555  lastTsos = (theMeas)->predictedState();
556  }
557  }
558  }
559  navLayers.clear();
560  updator()->makeFirstTime();
562 
563  measL.clear();
564 
565  return;
566 
567 }
568 
569 
570 //
571 // build trajectory in second hemisphere with pattern
572 // recognition starting from an intermediate state
573 //
575 
576  if ( (traj.firstMeasurement().recHit()->globalPosition().perp()
577  < traj.lastMeasurement().recHit()->globalPosition().perp()) ) {
578  LogTrace(category_)<<"inside-out: reverseTrajectory";
579  reverseTrajectory(traj);
580  }
581  if (traj.empty()) return;
583  if ( !tsos.isValid() ) tsos = traj.lastMeasurement().predictedState();
584  LogTrace(category_)<<"last tsos on traj: pos: "<< tsos.globalPosition()<<" mom: "<< tsos.globalMomentum();
585  if ( !tsos.isValid() ) {
586  LogTrace(category_)<<"last tsos on traj invalid";
587  return;
588  }
589 
590  build(intermediateState(tsos),insideOut,traj);
591 
592  return;
593 
594 }
595 
596 
597 //
598 //
599 //
601 
602  PerpendicularBoundPlaneBuilder planeBuilder;
603  GlobalPoint pos(0.0, 0.0, 0.0);
604  BoundPlane* SteppingPlane = planeBuilder(pos,tsos.globalDirection());
605 
606  TrajectoryStateOnSurface predTsos = propagator()->propagate(tsos, *SteppingPlane);
607  if ( predTsos.isValid() )
608  LogTrace(category_)<<"intermediateState: a intermediate state: pos: "<<predTsos.globalPosition() << "mom: " << predTsos.globalMomentum();
609 
610  return predTsos;
611 
612 }
613 
614 
615 //
616 //
617 //
619 
620  if ( hits.size() < 2 ) return;
621 
623  vector<bool> keep(hits.size(),true);
624  int i(0);
625  int j(0);
626 
627  for (MuonRecHitContainer::const_iterator ihit = hits.begin();
628  ihit != hits.end(); ++ihit ) {
629  if ( !keep[i] ) { i++; continue; };
630  j = i + 1;
631  for (MuonRecHitContainer::const_iterator ihit2 = ihit + 1;
632  ihit2 != hits.end(); ++ihit2 ) {
633  if ( !keep[j] ) { j++; continue; }
634  if ((*ihit)->geographicalId() == (*ihit2)->geographicalId() ) {
635  if ( (*ihit)->dimension() > (*ihit2)->dimension() ) {
636  keep[j] = false;
637  } else if ( (*ihit)->dimension() < (*ihit2)->dimension() ) {
638  keep[i] = false;
639  } else {
640  if ( (*ihit)->transientHits().size()>(*ihit2)->transientHits().size() ) {
641  keep[j] = false;
642  } else if ( (*ihit)->transientHits().size()<(*ihit2)->transientHits().size() ) {
643  keep[i] = false;
644  }
645  else if ( (*ihit)->degreesOfFreedom() != 0 && (*ihit2)->degreesOfFreedom() != 0) {
646  if (((*ihit)->chi2()/(*ihit)->degreesOfFreedom()) > ((*ihit2)->chi2()/(*ihit)->degreesOfFreedom())) keep[i] = false;
647  else keep[j] = false;
648  }
649  }
650  } // if same geomid
651  j++;
652  }
653  i++;
654  }
655 
656  i = 0;
657  for (MuonRecHitContainer::const_iterator ihit = hits.begin();
658  ihit != hits.end(); ++ihit ) {
659  if (keep[i] ) tmp.push_back(*ihit);
660  i++;
661  }
662 
663  hits.clear();
664  hits.swap(tmp);
665 
666  return;
667 
668 }
669 
670 
671 //
672 //
673 //
675 
677 
678  if (traj.empty()) return true;
679 
680  bool result = false;
681  for (ConstRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++ ) {
682  if ( !(*ir)->isValid() ) continue;
683  for (ConstRecHitContainer::const_iterator ir2 = ir+1; ir2 != hits.end(); ir2++ ) {
684  if ( !(*ir2)->isValid() ) continue;
685  if ( (*ir) == (*ir2) ) result = true;
686  }
687  }
688 
689  return result;
690 
691 }
692 
693 
694 //
695 // reverse a trajectory without refitting
696 // this can only be used for cosmic muons that come from above
697 //
699 
700  PropagationDirection newDir = (traj.firstMeasurement().recHit()->globalPosition().y()
701  < traj.lastMeasurement().recHit()->globalPosition().y())
703  Trajectory newTraj(traj.seed(), newDir);
704 
705  /* does not work in gcc4.8?)
706  std::vector<TrajectoryMeasurement> & meas = traj.measurements();
707  for (auto itm = meas.rbegin(); itm != meas.rend(); ++itm ) {
708  newTraj.push(std::move(*itm));
709  }
710  traj = std::move(newTraj);
711  */
712 
713  std::vector<TrajectoryMeasurement> const & meas = traj.measurements();
714  for (auto itm = meas.rbegin(); itm != meas.rend(); ++itm ) {
715  newTraj.push(*itm);
716  }
717  traj = newTraj;
718 
719 
720 }
721 
722 
723 //
724 // reverse a trajectory momentum direction and then refit
725 //
727 
729  if ( !lastTSOS.isValid() ) {
730  LogTrace(category_) << "Error: last TrajectoryState invalid.";
731  }
733  std::reverse(hits.begin(), hits.end());
734 
735  LogTrace(category_) << "last tsos before flipping "<<lastTSOS;
736  utilities()->reverseDirection(lastTSOS,&*theService->magneticField());
737  LogTrace(category_) << "last tsos after flipping "<<lastTSOS;
738 
739  vector<Trajectory> refittedback = theSmoother->fit(traj.seed(),hits,lastTSOS);
740  if ( refittedback.empty() ) {
741  LogTrace(category_) <<"flipTrajectory fail. "<<endl;
742  return;
743  }
744  LogTrace(category_) <<"flipTrajectory: first "<< refittedback.front().firstMeasurement().updatedState()
745  <<"\nflipTrajectory: last "<<refittedback.front().lastMeasurement().updatedState();
746 
747  traj = refittedback.front();
748 
749  return;
750 
751 }
752 
753 
754 //
755 //
756 //
758 
759  if ( traj.direction() == anyDirection ) return;
761  Trajectory newTraj(traj.seed(), newDir);
762  const std::vector<TrajectoryMeasurement>& meas = traj.measurements();
763 
764  for (std::vector<TrajectoryMeasurement>::const_iterator itm = meas.begin(); itm != meas.end(); ++itm) {
765  newTraj.push(*itm);
766  }
767 
768  while (!traj.empty()) {
769  traj.pop();
770  }
771 
772  traj = newTraj;
773 
774 }
775 
776 
777 //
778 // guess the direction by normalized chi2
779 //
781 
783 
785 
787 
788  if ( !firstTSOS.isValid() || !lastTSOS.isValid() ) return;
789 
790  LogTrace(category_) <<"Two ends of the traj "<<firstTSOS.globalPosition()
791  <<", "<<lastTSOS.globalPosition();
792 
793  LogTrace(category_) <<"Their mom: "<<firstTSOS.globalMomentum()
794  <<", "<<lastTSOS.globalMomentum();
795 
796  LogTrace(category_) <<"Their mom eta: "<<firstTSOS.globalMomentum().eta()
797  <<", "<<lastTSOS.globalMomentum().eta();
798 
799  // momentum eta can be used to estimate direction
800  // the beam-halo muon seems enter with a larger |eta|
801 
802  if ( fabs(firstTSOS.globalMomentum().eta()) > fabs(lastTSOS.globalMomentum().eta()) ) {
803 
804  vector<Trajectory> refitted = theSmoother->trajectories(traj.seed(),hits,firstTSOS);
805  if ( !refitted.empty() ) traj = refitted.front();
806 
807  } else {
808  std::reverse(hits.begin(), hits.end());
809  utilities()->reverseDirection(lastTSOS,&*theService->magneticField());
810  vector<Trajectory> refittedback = theSmoother->trajectories(traj.seed(),hits,lastTSOS);
811  if ( !refittedback.empty() ) traj = refittedback.front();
812 
813  }
814 
815  return;
816 
817 }
818 
819 
820 //
821 // get direction from timing information of rechits and segments
822 //
824 
826  LogTrace(category_) << "getDirectionByTime"<<endl;
827  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++ ) {
828  if ( !(*ir)->isValid() ) {
829  LogTrace(category_) << "invalid RecHit"<<endl;
830  continue;
831  }
832 
833  const GlobalPoint& pos = (*ir)->globalPosition();
835  << "pos" << pos
836  << "radius " << pos.perp()
837  << " dim " << (*ir)->dimension()
838  << " det " << (*ir)->det()->geographicalId().det()
839  << " sub det " << (*ir)->det()->subDetector()<<endl;
840 
841  if ((*ir)->det()->geographicalId().det() == 2 && (*ir)->det()->subDetector() == 6) {
842 // const CSCRecHit2D* iCSC = dynamic_cast<const CSCRecHit2D*>(&**ir);
843 // if (iCSC) LogTrace(category_)<<"csc from cast tpeak "<<iCSC->tpeak();
844  CSCRecHit2DCollection::range thisrange = cschits_->get(CSCDetId((*ir)->geographicalId()));
845  for (CSCRecHit2DCollection::const_iterator rechit = thisrange.first; rechit!=thisrange.second;++rechit) {
846  if ((*rechit).isValid()) LogTrace(category_)<<"csc from collection tpeak "<<(*rechit).tpeak();
847  }
848  }
849  if ((*ir)->det()->geographicalId().det() == 2 && (*ir)->det()->subDetector() == 7) {
850 // const DTRecHit1D* iDT = dynamic_cast<const DTRecHit1D*>(&**ir);
851 // if (iDT) LogTrace(category_)<<"dt digitime "<<iDT->digiTime();
852  DTRecHitCollection::range thisrange = dthits_->get(DTLayerId((*ir)->geographicalId()));
853  for (DTRecHitCollection::const_iterator rechit = thisrange.first; rechit!=thisrange.second;++rechit) {
854  if ((*rechit).isValid()) LogTrace(category_)<<"dt from collection digitime "<<(*rechit).digiTime();
855  }
856  }
857  }
858 
859  return;
860 
861 }
862 
863 
864 //
865 //
866 //
867 std::vector<TrajectoryMeasurement>
869  const TrajectoryStateOnSurface& tsos,
870  const Propagator* propagator,
872 
873  std::vector<TrajectoryMeasurement> result;
874  std::vector<TrajectoryMeasurement> measurements;
875 
876  if ( layer->hasGroups() ) {
877  std::vector<TrajectoryMeasurementGroup> measurementGroups =
878  theLayerMeasurements->groupedMeasurements(layer, tsos, *propagator, *estimator);
879 
880  for (std::vector<TrajectoryMeasurementGroup>::const_iterator tmGroupItr = measurementGroups.begin();
881  tmGroupItr != measurementGroups.end(); ++tmGroupItr) {
882 
883  measurements = tmGroupItr->measurements();
884  const TrajectoryMeasurement* bestMeasurement
885  = theBestMeasurementFinder->findBestMeasurement(measurements, propagator);
886 
887  if (bestMeasurement) result.push_back(*bestMeasurement);
888  }
889  }
890  else {
891  measurements = theLayerMeasurements->measurements(layer, tsos, *propagator, *estimator);
892  const TrajectoryMeasurement* bestMeasurement
893  = theBestMeasurementFinder->findBestMeasurement(measurements, propagator);
894 
895  if (bestMeasurement) result.push_back(*bestMeasurement);
896  }
897  measurements.clear();
898 
899  return result;
900 
901 }
902 
903 
904 //
905 //
906 //
908  int& dtChambers,
909  int& cscChambers,
910  int& rpcChambers,
911  int& totalChambers) {
912 
913  if (layer->subDetector()==GeomDetEnumerators::DT) dtChambers++;
914  else if (layer->subDetector()==GeomDetEnumerators::CSC) cscChambers++;
915  else if (layer->subDetector()==GeomDetEnumerators::RPCBarrel || layer->subDetector()==GeomDetEnumerators::RPCEndcap) rpcChambers++;
916  totalChambers++;
917 
918 }
919 
920 
921 //
922 //
923 //
925 
926  if ( (dtseg == 0) || (!dtseg->hasPhi()) ) return 0;
927  // timing information
928  double result = 0;
929  if ( dtseg->phiSegment() == 0 ) return 0;
930  int phiHits = dtseg->phiSegment()->specificRecHits().size();
931  LogTrace(category_) << "phiHits " << phiHits;
932  if ( phiHits > 5 ) {
933  if(dtseg->phiSegment()->ist0Valid()) result = dtseg->phiSegment()->t0();
934  if (dtseg->phiSegment()->ist0Valid()){
935  LogTrace(category_) << " Phi t0: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
936  } else {
937  LogTrace(category_) << " Phi t0 is invalid: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
938  }
939  }
940 
941  return result;
942 
943 }
944 
945 
946 //
947 //
948 //
950  const DTRecSegment4D* dtseg2) const {
951 
952  LogTrace(category_) << "comparing dtseg: " << dtseg1 << " " << dtseg2 << endl;
953  if (dtseg1 == dtseg2 || t0(dtseg1) == t0(dtseg2)) return anyDirection;
954 
956  ( t0(dtseg1) < t0(dtseg2) ) ? alongMomentum : oppositeToMomentum;
957 
958  return result;
959 
960 }
bool empty() const
True if trajectory has no measurements.
Definition: Trajectory.h:299
void reverseTrajectoryPropagationDirection(Trajectory &) const
reverse the propagation direction of a trajectory
T getParameter(std::string const &) const
const Propagator * propagator() const
int i
Definition: DBlmapReader.cc:9
virtual TrajectoryContainer trajectories(const Trajectory &traj) const override
double t0(const DTRecSegment4D *deseg) const
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
TrajectoryStateOnSurface const & predictedState() const
const MeasurementEstimator * estimator() const
accasso at the propagator
T perp() const
Definition: PV3DBase.h:72
ConstRecHitPointer const & recHit() const
TrajectorySeed const & seed() const
Access to the seed used to reconstruct the Trajectory.
Definition: Trajectory.h:330
std::vector< const DetLayer * > compatibleEndcapLayers(const FreeTrajectoryState &fts, PropagationDirection timeDirection) const
void estimateDirection(Trajectory &) const
check the direction of trajectory by checking eta spread
void makeFirstTime()
reset the theFirstTSOSFlag
TrajectoryStateOnSurface intermediateState(const TrajectoryStateOnSurface &) const
edm::Handle< CSCRecHit2DCollection > cschits_
const DTChamberRecSegment2D * phiSegment() const
The superPhi segment: 0 if no phi projection available.
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
virtual SubDetector subDetector() const =0
The type of detector (PixelBarrel, PixelEndcap, TIB, TOB, TID, TEC, CSC, DT, RPCBarrel, RPCEndcap)
MuonTrajectoryUpdator * backwardUpdator() const
GlobalPoint globalPosition() const
ConstRecHitContainer recHits() const
Definition: Trajectory.h:258
std::vector< Trajectory * > trajectories(const TrajectorySeed &)
build trajectories from seed
MuonCandidate::TrajectoryContainer TrajectoryContainer
edm::Handle< DTRecHitCollection > dthits_
const Propagator * propagatorOpposite() const
PropagationDirection
MuonDetLayerMeasurements * theLayerMeasurements
const CosmicMuonUtilities * utilities() const
std::vector< TrajectoryMeasurementGroup > groupedMeasurements(const DetLayer *layer, const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est, const edm::Event &iEvent)
std::vector< const DetLayer * > compatibleLayers(const FreeTrajectoryState &fts, PropagationDirection timeDirection) const
tuple result
Definition: mps_fire.py:84
PropagationDirection const & direction() const
Definition: Trajectory.cc:125
DataContainer const & measurements() const
Definition: Trajectory.h:250
const int keep
DirectMuonNavigation * navigation() const
static const int CSC
Definition: MuonSubdetId.h:13
T mag() const
Definition: PV3DBase.h:67
void flipTrajectory(Trajectory &) const
flip a trajectory with refit (the momentum direction is opposite)
MuonTrajectoryUpdator * updator() const
MuonRecHitContainer recHits(const DetLayer *layer, const edm::Event &iEvent)
returns the rechits which are on the layer
void reverseTrajectory(Trajectory &) const
reverse a trajectory without refit (out the measurements order changed)
MuonBestMeasurementFinder * theBestMeasurementFinder
void setFitDirection(NavigationDirection fitDirection)
set fit direction
TrajectoryMeasurement const & lastMeasurement() const
Definition: Trajectory.h:228
FreeTrajectoryState const * freeState(bool withErrors=true) const
T z() const
Definition: PV3DBase.h:64
int j
Definition: DBlmapReader.cc:9
void incrementChamberCounters(const DetLayer *layer, int &dtChambers, int &cscChambers, int &rpcChambers, int &totalChambers)
bool hasPhi() const
Does it have the Phi projection?
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
MeasurementContainer measurements(const DetLayer *layer, const GeomDet *det, const TrajectoryStateOnSurface &stateOnDet, const MeasurementEstimator &est, const edm::Event &iEvent)
GlobalVector momentum() const
#define LogTrace(id)
std::vector< DTRecHit1D > specificRecHits() const
Access to specific components.
bool ist0Valid() const
std::vector< ConstRecHitPointer > ConstRecHitContainer
void getDirectionByTime(Trajectory &) const
check the direction of trajectory by checking the timing
Definition: DetId.h:18
GlobalPoint position() const
PropagationDirection checkDirectionByT0(const DTRecSegment4D *, const DTRecSegment4D *) const
bool isValid() const
Definition: Trajectory.h:324
TrajectoryMeasurement const & firstMeasurement() const
Definition: Trajectory.h:241
PTrajectoryStateOnDet const & startingState() const
#define debug
Definition: HDRShower.cc:19
TrajectoryStateOnSurface transientState(const PTrajectoryStateOnDet &ts, const Surface *surface, const MagneticField *field)
void pop()
Definition: Trajectory.cc:16
std::vector< Trajectory > fit(const Trajectory &) const
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
float chiSquared() const
Definition: Trajectory.h:307
virtual ~CosmicMuonTrajectoryBuilder()
Destructor.
T eta() const
Definition: PV3DBase.h:76
void selectHits(MuonTransientTrackingRecHit::MuonRecHitContainer &) const
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
void setEvent(const edm::Event &)
set event
GlobalVector globalMomentum() const
if(dp >Float(M_PI)) dp-
CosmicMuonTrajectoryBuilder(const edm::ParameterSet &, const MuonServiceProxy *service, edm::ConsumesCollector &iC)
Constructor.
const Propagator * propagatorAlong() const
TrajectoryMeasurement * findBestMeasurement(std::vector< TrajectoryMeasurement > &measC, const Propagator *propagator)
return the Tm with the best chi2: no cut applied.
MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer
std::vector< TrajectoryMeasurement > findBestMeasurements(const DetLayer *, const TrajectoryStateOnSurface &, const Propagator *, const MeasurementEstimator *)
virtual std::pair< bool, TrajectoryStateOnSurface > update(const TrajectoryMeasurement *measurement, Trajectory &trajectory, const Propagator *propagator)
update the Trajectory with the TrajectoryMeasurement
bool selfDuplicate(const Trajectory &) const
check if the trajectory iterates the same hit more than once
static const int DT
Definition: MuonSubdetId.h:12
TrajectoryStateOnSurface const & updatedState() const
dbl *** dir
Definition: mlp_gen.cc:35
MuonTransientTrackingRecHit::MuonRecHitContainer unusedHits(const DetLayer *, const TrajectoryMeasurement &) const
double t0() const
Get the segment t0 (if recomputed, 0 is returned otherwise)
void reverseDirection(TrajectoryStateOnSurface &, const MagneticField *) const
void build(const TrajectoryStateOnSurface &, const NavigationDirection &, Trajectory &)
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:56
std::vector< MuonRecHitPointer > MuonRecHitContainer
GlobalVector globalDirection() const
virtual void setEvent(const edm::Event &)
pass the Event to the algo at each event
T dot(const Basic3DVector &rh) const
Scalar product, or &quot;dot&quot; product, with a vector of same type.