CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DynamicTruncation.cc
Go to the documentation of this file.
1 
33 
34 
35 #define MAX_THR 1e7
36 
37 using namespace edm;
38 using namespace std;
39 using namespace reco;
40 
42  DTThr(0), CSCThr(0)
43 {
44  theEvent = &event;
45  theSetup = &theService.eventSetup();
46  propagator = theService.propagator("SmartPropagator");
47  propagatorCompatibleDet = theService.propagator("SteppingHelixPropagatorAny");
48  theG = theService.trackingGeometry();
49  theService.eventSetup().get<TransientRecHitRecord>().get("MuonRecHitBuilder",theMuonRecHitBuilder);
50  theService.eventSetup().get<TrackingComponentsRecord>().get("KFUpdator",updatorHandle);
51  theService.eventSetup().get<MuonGeometryRecord>().get(cscGeom);
52  theService.eventSetup().get<MuonRecoGeometryRecord>().get(navMuon);
53  theService.eventSetup().get<IdealMagneticFieldRecord>().get(magfield);
54  navigation = new DirectMuonNavigation(theService.detLayerGeometry());
55 }
56 
57 
58 
60  if (navigation) delete navigation;
61 }
62 
63 
64 
66  result.clear();
67  // Put the tracker hits in the final vector and get the last tracker valid measure
68  std::vector<TrajectoryMeasurement> muonMeasurements = traj.measurements();
69  TrajectoryMeasurement lastTKm = muonMeasurements.front();
70  for (std::vector<TrajectoryMeasurement>::const_iterator imT = muonMeasurements.begin(); imT != muonMeasurements.end(); imT++ ) {
71  if ( !(*imT).recHit()->isValid() ) continue;
72  const TransientTrackingRecHit* hit = &(*(*imT).recHit());
73  if (hit->geographicalId().det() == DetId::Tracker) {
74  result.push_back((*imT).recHit());
75  if (!(*imT).forwardPredictedState().isValid()) continue;
76  if ((*imT).forwardPredictedState().globalPosition().mag()
77  > lastTKm.forwardPredictedState().globalPosition().mag()) lastTKm = *imT;
78  }
79  }
80 
81  // get the last (forward) predicted state for the tracker
83 
84  // update the state with the last tracker measure
85  update(currentState, lastTKm.recHit());
86 
87  // use the navigation to get the list of compatible dets
88  map<int, std::vector<DetId> > detMap;
90 
91  // selects the muon hits for the final refit
92  filteringAlgo(detMap);
93 
94  return result;
95 }
96 
97 
98 
99 void DynamicTruncation::setThr(int bThr, int eThr) {
100  if (bThr <= MAX_THR && bThr >= 0) DTThr = bThr; // DT thr
101  else DTThr = MAX_THR;
102  if (eThr <= MAX_THR && eThr >= 0) CSCThr = eThr; // CSC thr
103  else CSCThr = MAX_THR;
104 }
105 
106 
107 
108 double DynamicTruncation::getBest(std::vector<CSCSegment>& segs, TrajectoryStateOnSurface& tsos, CSCSegment& bestCSCSeg) {
109  unsigned int i = 0;
110  double val = MAX_THR;
111  std::vector<CSCSegment>::size_type sz = segs.size();
112  for (i=0; i<sz; i++) {
113  StateSegmentMatcher estim(&tsos, &segs[i]);
114  double tmp = estim.value();
115  if (tmp < val) {
116  bestCSCSeg = segs[i];
117  val = tmp;
118  }
119  }
120  return val;
121 }
122 
123 
124 
125 double DynamicTruncation::getBest(std::vector<DTRecSegment4D>& segs, TrajectoryStateOnSurface& tsos, DTRecSegment4D& bestDTSeg) {
126  unsigned int i = 0;
127  double val = MAX_THR;
129  for (i=0; i<sz; i++) {
130  StateSegmentMatcher estim(&tsos, &segs[i]);
131  double tmp = estim.value();
132  if (tmp < val) {
133  bestDTSeg = segs[i];
134  val = tmp;
135  }
136  }
137  return val;
138 }
139 
140 
141 
142 void DynamicTruncation::compatibleDets(TrajectoryStateOnSurface& tsos, map<int, std::vector<DetId> >& detMap) {
143  // SteppingHelixPropagator prop(magfield.product(), anyDirection);
144  // MuonPatternRecoDumper dumper;
145  MeasurementEstimator *theEstimator = new Chi2MeasurementEstimator(100, 3);
146  std::vector<const DetLayer *> navLayers;
148  unsigned int nlayer = 0;
149  for ( unsigned int ilayer=0; ilayer<navLayers.size(); ilayer++ ) {
150  // Skip RPC layers
151  if (navLayers[ilayer]->subDetector() == GeomDetEnumerators::RPCEndcap
152  || navLayers[ilayer]->subDetector() == GeomDetEnumerators::RPCBarrel) continue;
153  std::vector<DetLayer::DetWithState> comps = navLayers[ilayer]->compatibleDets(currentState,
154  *propagatorCompatibleDet, *theEstimator);
155  // cout << comps.size() << " compatible Dets with " << navLayers[ilayer]->subDetector() << " Layer " << ilayer << " "
156  // << dumper.dumpLayer(navLayers[ilayer]) << " " << endl;
157  if (comps.size() > 0) {
158  DetId id(comps.front().first->geographicalId().rawId());
159  detMap[nlayer].push_back(id);
160  }
161  }
162  if (theEstimator) delete theEstimator;
163 }
164 
165 
166 
167 void DynamicTruncation::filteringAlgo(map<int, std::vector<DetId> >& detMap) {
169  for (unsigned int iDet = 0; iDet < detMap.size(); ++iDet) {
170  double bestLayerValue = MAX_THR;
171  ConstRecHitContainer layerRH;
172  std::vector<DetId> chamber = detMap[iDet];
173  for (unsigned int j = 0; j < chamber.size(); ++j) {
174  DetId id = chamber[j];
175 
176  if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::DT) {
177  DTChamberId DTid(id);
178 
179  std::vector<DTRecSegment4D> allDTsegs;
180  std::map<int, std::vector<DTRecSegment4D> >::const_iterator dtIter = getSegs.getDTlist().find(DTid.station());
181  if (dtIter != getSegs.getDTlist().end()){
182  allDTsegs = dtIter->second;
183  }
184 
185  std::vector<DTRecSegment4D>::size_type sz = allDTsegs.size();
186  for (unsigned int iSeg=0; iSeg<sz; ++iSeg) {
187 
188  // Propagate the state to the current chamber
189  TrajectoryStateOnSurface tsosdt = propagator->propagate(currentState,
190  theG->idToDet(allDTsegs[iSeg].chamberId())->surface());
191  if (!tsosdt.isValid()) continue;
192 
193  std::vector<DTRecSegment4D> DTsegs;
194  DTsegs.push_back(allDTsegs[iSeg]);
195 
196  DTRecSegment4D bestDTSeg;
197  double bestChamberValue = getBest(DTsegs, tsosdt, bestDTSeg);
198  if (bestChamberValue < bestLayerValue) bestLayerValue = bestChamberValue;
199 
200  // Check if the best estimator value is below the THR and then get the RH componing the segment
201  if (bestChamberValue >= DTThr || bestChamberValue > bestLayerValue) continue;
202  layerRH.clear();
203  std::vector<DTRecHit1D> DTrh = getSegs.getDTRHmap(bestDTSeg);
204  for (std::vector<DTRecHit1D>::iterator it = DTrh.begin(); it != DTrh.end(); it++) {
205  layerRH.push_back(theMuonRecHitBuilder->build(&*it));
206  }
207  }
208  }
209  if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
210  CSCDetId CSCid(id);
211 
212  std::vector<CSCSegment> allCSCsegs;
213  std::map<int, std::vector<CSCSegment> >::const_iterator cscIter = getSegs.getCSClist().find(CSCid.station());
214  if (cscIter != getSegs.getCSClist().end()){
215  allCSCsegs = cscIter->second;
216  }
217 
218  std::vector<CSCSegment>::size_type sz = allCSCsegs.size();
219  for (unsigned int iSeg=0; iSeg<sz; ++iSeg) {
220 
221  // Propagate the state to the current chamber
222  TrajectoryStateOnSurface tsoscsc = propagator->propagate(currentState,
223  theG->idToDet(allCSCsegs[iSeg].cscDetId())->surface());
224  if (!tsoscsc.isValid()) continue;
225 
226  std::vector<CSCSegment> CSCsegs;
227  CSCsegs.push_back(allCSCsegs[iSeg]);
228 
229  CSCSegment bestCSCSeg;
230  double bestChamberValue = getBest(CSCsegs, tsoscsc, bestCSCSeg);
231  if (bestChamberValue < bestLayerValue) bestLayerValue = bestChamberValue;
232 
233  // Check if the best estimator value is below the THR and then get the RH componing the segment
234  if (bestChamberValue >= CSCThr || bestChamberValue > bestLayerValue) continue;
235  layerRH.clear();
236 
237  std::vector<CSCRecHit2D> CSCrh = getSegs.getCSCRHmap(bestCSCSeg);
238  for (std::vector<CSCRecHit2D>::iterator it = CSCrh.begin(); it != CSCrh.end(); ++it) {
239  layerRH.push_back(theMuonRecHitBuilder->build(&*it));
240  }
241  }
242  }
243  }
244 
245  if (layerRH.size() > 0) {
246  for (ConstRecHitContainer::iterator it = layerRH.begin(); it != layerRH.end(); ++it) {
247  result.push_back((*it));
248  }
249 
250  // sort the vector
251  layerRH = sort(layerRH);
252 
253  // update the currentState using all rh
254  DetId id = layerRH.front()->geographicalId();
255  if (id.subdetId() == MuonSubdetId::DT) updateWithDThits(layerRH);
256  else updateWithCSChits(layerRH);
257  }
258  layerRH.clear();
259  }
260 }
261 
262 
263 
265  TrajectoryStateOnSurface temp = updatorHandle->update(tsos, *rechit);
266  if (temp.isValid()) currentState = updatorHandle->update(tsos, *rechit);
267 }
268 
269 
270 
272  for (ConstRecHitContainer::const_iterator it = recHits.begin(); it != recHits.end(); ++it) {
273  DTLayerId layid((*it)->det()->geographicalId());
274  TrajectoryStateOnSurface temp = propagator->propagate(currentState, theG->idToDet(layid)->surface());
275  if (temp.isValid()) currentState = updatorHandle->update(temp, **it);
276  }
277 }
278 
279 
280 
282  for (ConstRecHitContainer::const_iterator it = recHits.begin(); it != recHits.end(); ++it) {
283  const CSCLayer* cscChamber = cscGeom->layer((*it)->det()->geographicalId());
284  CSCDetId layid = cscChamber->geographicalId();
285  TrajectoryStateOnSurface temp = propagator->propagate(currentState, theG->idToDet(layid)->surface());
286  if (temp.isValid()) currentState = updatorHandle->update(temp, **it);
287  }
288 }
289 
290 
291 
293  unsigned int i=0;
294  unsigned int j=0;
295  ConstRecHitContainer::size_type n = recHits.size();
296  for(i=1; i<n; ++i)
297  for(j=n-1; j>=i; --j)
298  {
299  if(recHits[j-1]->globalPosition().mag() > recHits[j]->globalPosition().mag())
300  {
301  swap (recHits[j-1],recHits[j]);
302  }
303  }
304  return recHits;
305 }
306 
void swap(ora::Record &rh, ora::Record &lh)
Definition: Record.h:74
DirectMuonNavigation * navigation
int i
Definition: DBlmapReader.cc:9
void compatibleDets(TrajectoryStateOnSurface &, std::map< int, std::vector< DetId > > &)
TrajectoryStateOnSurface currentState
edm::ESHandle< Propagator > propagatorCompatibleDet
TrajectoryStateOnSurface forwardPredictedState() const
Access to forward predicted state (from fitter or builder)
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
GlobalPoint globalPosition() const
ConstRecHitPointer recHit() const
std::vector< CSCRecHit2D > getCSCRHmap(CSCSegment)
uint16_t size_type
std::vector< const DetLayer * > compatibleLayers(const FreeTrajectoryState &fts, PropagationDirection timeDirection) const
DataContainer const & measurements() const
Definition: Trajectory.h:169
void filteringAlgo(std::map< int, std::vector< DetId > > &)
static const int CSC
Definition: MuonSubdetId.h:15
T mag() const
Definition: PV3DBase.h:61
std::vector< DTRecHit1D > getDTRHmap(DTRecSegment4D)
edm::ESHandle< MuonDetLayerGeometry > navMuon
const std::map< int, std::vector< CSCSegment > > & getCSClist() const
ConstRecHitContainer result
FreeTrajectoryState * freeState(bool withErrors=true) const
const edm::EventSetup * theSetup
TransientTrackingRecHit::ConstRecHitContainer ConstRecHitContainer
double getBest(std::vector< CSCSegment > &, TrajectoryStateOnSurface &, CSCSegment &)
int j
Definition: DBlmapReader.cc:9
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:74
edm::ESHandle< TrajectoryStateUpdator > updatorHandle
TransientTrackingRecHit::ConstRecHitContainer filter(const Trajectory &)
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
const edm::Event * theEvent
edm::ESHandle< TransientTrackingRecHitBuilder > theMuonRecHitBuilder
#define MAX_THR
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: DetId.h:20
DynamicTruncation(const edm::Event &, const MuonServiceProxy &)
void update(TrajectoryStateOnSurface &, ConstRecHitPointer)
ConstRecHitContainer sort(ConstRecHitContainer &)
edm::ESHandle< Propagator > propagator
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
void setThr(int, int)
const std::map< int, std::vector< DTRecSegment4D > > & getDTlist() const
void updateWithDThits(ConstRecHitContainer &)
void updateWithCSChits(ConstRecHitContainer &)
int station() const
Definition: CSCDetId.h:88
static const int DT
Definition: MuonSubdetId.h:14
DetId geographicalId() const
int station() const
Return the station number.
Definition: DTChamberId.h:53
Detector det() const
get the detector field from this detid
Definition: DetId.h:37
edm::ESHandle< MagneticField > magfield
edm::ESHandle< GlobalTrackingGeometry > theG
edm::ESHandle< CSCGeometry > cscGeom