CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrajectorySegmentBuilder.cc
Go to the documentation of this file.
1 
3 
6 
19 
20 #include <algorithm>
21 
22 // #define DBG_TSB
23 
24 using namespace std;
25 
28 {
29  //
30  // create empty trajectory
31  //
32  theLockedHits.clear();
33  TempTrajectory startingTrajectory(theFullPropagator.propagationDirection());
34  //
35  // get measurement groups
36  //
37  vector<TMG> measGroups =
38  theLayerMeasurements->groupedMeasurements(theLayer,startingState,theFullPropagator,theEstimator);
39  //B.M. theLayer.groupedMeasurements(startingState,theFullPropagator,theEstimator);
40 
41 #ifdef DBG_TSB
42  cout << "TSB: number of measurement groups = " << measGroups.size() << endl;
43  // theDbgFlg = measGroups.size()>1;
44  theDbgFlg = true;
45 #else
46  theDbgFlg = false;
47 #endif
48 
49  if ( theDbgFlg ) {
50  int ntot(1);
51  for (vector<TMG>::const_iterator ig=measGroups.begin();
52  ig!=measGroups.end(); ++ig) {
53  int ngrp(0);
54  const vector<TM>& measurements = ig->measurements();
55  for ( vector<TM>::const_iterator im=measurements.begin();
56  im!=measurements.end(); ++im ) {
57  if ( im->recHit()->isValid() ) ngrp++;
58  }
59  cout << " " << ngrp;
60  if ( ngrp>0 ) ntot *= ngrp;
61  }
62  cout << endl;
63  cout << "TrajectorySegmentBuilder::partialTrajectories:: det ids & hit types / group" << endl;
64  for (vector<TMG>::const_iterator ig=measGroups.begin();
65  ig!=measGroups.end(); ++ig) {
66  const vector<TM>& measurements = ig->measurements();
67  for ( vector<TM>::const_iterator im=measurements.begin();
68  im!=measurements.end(); ++im ) {
69  if ( im!=measurements.begin() ) cout << " / ";
70  if ( im->recHit()->det() )
71  cout << im->recHit()->det()->geographicalId().rawId() << " "
72  << im->recHit()->getType();
73  else
74  cout << "no det";
75  }
76  cout << endl;
77  }
78  }
79 
80 #ifdef DBG_TSB
81 // if ( measGroups.size()>4 ) {
82  cout << typeid(theLayer).name() << endl;
83  cout << startingState.localError().matrix() << endl;
84 // for (vector<TMG>::const_iterator ig=measGroups.begin();
85 // ig!=measGroups.end(); ig++) {
86 // cout << "Nr. of measurements = " << ig->measurements().size() << endl;
87 // const DetGroup& dg = ig->detGroup();
88 // for ( DetGroup::const_iterator id=dg.begin();
89 // id!=dg.end(); id++ ) {
90 // GlobalPoint p(id->det()->position());
91 // GlobalVector v(id->det()->toGlobal(LocalVector(0.,0.,1.)));
92 // cout << p.perp() << " " << p.phi() << " " << p.z() << " ; "
93 // << v.phi() << " " << v.z() << endl;
94 // }
95 // }
96 // }
97 #endif
98 
99  TempTrajectoryContainer candidates =
100  addGroup(startingTrajectory,measGroups.begin(),measGroups.end());
101 
102  if (theDbgFlg) cout << "TSB: back with " << candidates.size() << " candidates" << endl;
103  // clean
104  //
105  //
106  // add invalid hit - try to get first detector hit by the extrapolation
107  //
108 
109  updateWithInvalidHit(startingTrajectory,measGroups,candidates);
110  if (theDbgFlg) cout << "TSB: " << candidates.size() << " candidates after invalid hit" << endl;
111  return candidates;
112 }
113 
115  const TM& tm) const
116 {
117  TSOS predictedState = tm.predictedState();
119 
120  if ( hit->isValid()) {
121  traj.push( TM( predictedState, theUpdator.update( predictedState, *hit),
122  hit, tm.estimate(), tm.layer()));
123 // TrajectoryMeasurement tm(traj.lastMeasurement());
124 // if ( tm.updatedState().isValid() ) {
125 // if ( !hit.det().surface().bounds().inside(tm.updatedState().localPosition(),
126 // tm.updatedState().localError().positionError(),3.) ) {
127 // cout << "Incompatibility after update for det at " << hit.det().position() << ":" << endl;
128 // cout << tm.predictedState().localPosition() << " "
129 // << tm.predictedState().localError().positionError() << endl;
130 // cout << hit.localPosition() << " " << hit.localPositionError() << endl;
131 // cout << tm.updatedState().localPosition() << " "
132 // << tm.updatedState().localError().positionError() << endl;
133 // }
134 // }
135  }
136  else {
137  traj.push( TM( predictedState, hit,0, tm.layer()));
138  }
139 }
140 
141 
144  vector<TMG>::const_iterator begin,
145  vector<TMG>::const_iterator end)
146 {
147  vector<TempTrajectory> ret;
148  if ( begin==end ) {
149  //std::cout << "TrajectorySegmentBuilder::addGroup" << " traj.empty()=" << traj.empty() << "EMPTY" << std::endl;
150  if (theDbgFlg) cout << "TSB::addGroup : no groups left" << endl;
151  if ( !traj.empty() )
152  ret.push_back(traj);
153  return ret;
154  }
155 
156  if (theDbgFlg) cout << "TSB::addGroup : traj.size() = " << traj.measurements().size()
157  << " first group at " << &(*begin)
158  // << " nr. of candidates = " << candidates.size()
159  << endl;
160 
161 
162  TempTrajectoryContainer updatedTrajectories; updatedTrajectories.reserve(2);
163  if ( traj.measurements().empty() ) {
164  vector<TM> firstMeasurements = unlockedMeasurements(begin->measurements());
165  if ( theBestHitOnly )
166  updateCandidatesWithBestHit(traj,firstMeasurements,updatedTrajectories);
167  else
168  updateCandidates(traj,begin->measurements(),updatedTrajectories);
169  if (theDbgFlg) cout << "TSB::addGroup : updating with first group - "
170  << updatedTrajectories.size() << " trajectories" << endl;
171  }
172  else {
173  if ( theBestHitOnly )
174  updateCandidatesWithBestHit(traj,redoMeasurements(traj,begin->detGroup()),
175  updatedTrajectories);
176  else
177  updateCandidates(traj,redoMeasurements(traj,begin->detGroup()),
178  updatedTrajectories);
179  if (theDbgFlg) cout << "TSB::addGroup : updating"
180  << updatedTrajectories.size() << " trajectories" << endl;
181  }
182 
183  if (begin+1 != end) {
184  ret.reserve(4); // a good upper bound
185  for ( TempTrajectoryContainer::iterator it=updatedTrajectories.begin();
186  it!=updatedTrajectories.end(); ++it ) {
187  if (theDbgFlg) cout << "TSB::addGroup : trying to extend candidate at "
188  << &(*it) << " size " << it->measurements().size() << endl;
189  vector<TempTrajectory> finalTrajectories = addGroup(*it,begin+1,end);
190  if (theDbgFlg) cout << "TSB::addGroup : " << finalTrajectories.size()
191  << " finalised candidates before cleaning" << endl;
192  //B.M. to be ported later
193  cleanCandidates(finalTrajectories);
194 
195  if (theDbgFlg) cout << "TSB::addGroup : got " << finalTrajectories.size()
196  << " finalised candidates" << endl;
197  ret.insert(ret.end(),finalTrajectories.begin(),
198  finalTrajectories.end());
199  }
200  } else {
201  ret.reserve(updatedTrajectories.size());
202  for (TempTrajectoryContainer::iterator it=updatedTrajectories.begin();
203  it!=updatedTrajectories.end(); ++it ) {
204  if (!it->empty()) ret.push_back(*it);
205  }
206  }
207 
208  //std::cout << "TrajectorySegmentBuilder::addGroup" <<
209  // " traj.empty()=" << traj.empty() <<
210  // " end-begin=" << (end-begin) <<
211  // " #updated=" << updatedTrajectories.size() <<
212  // " #result=" << ret.size() << std::endl;
213  return ret;
214 }
215 
216 void
218  const vector<TM>& measurements,
219  TempTrajectoryContainer& candidates)
220 {
221  //
222  // generate updated candidates with all valid hits
223  //
224  for ( vector<TM>::const_iterator im=measurements.begin();
225  im!=measurements.end(); ++im ) {
226  if ( im->recHit()->isValid() ) {
227  candidates.push_back(traj);
228  updateTrajectory(candidates.back(),*im);
229  if ( theLockHits ) lockMeasurement(*im);
230  }
231  }
232  //
233  // keep old trajectory
234  //
235  candidates.push_back(traj);
236 }
237 
238 void
240  const vector<TM>& measurements,
241  TempTrajectoryContainer& candidates)
242 {
243  vector<TM>::const_iterator ibest = measurements.end();
244  bool bestIsValid = false;
245  for ( vector<TM>::const_iterator im=measurements.begin();
246  im!=measurements.end(); ++im ) {
247  if ( im->recHit()->isValid() ) {
248  if (!bestIsValid || (ibest==measurements.end()) || (im->estimate()<ibest->estimate()) ) {
249  ibest = im;
250  bestIsValid = true;
251  }
252  }
253 // else if (!bestIsValid && (im->recHit()->getType() != TrackingRecHit::missing)) {
254 // if (ibest==measurements.end()) ibest = im;
255 // }
256  }
257  if ( ibest!=measurements.end() ) {
258  if ( theLockHits ) lockMeasurement(*ibest);
259  candidates.push_back(traj);
260  updateTrajectory(candidates.back(),*ibest);
261  if ( theDbgFlg ) {
262  if (bestIsValid) {
263  cout << "TSB: found best measurement at "
264  << ibest->recHit()->globalPosition().perp() << " "
265  << ibest->recHit()->globalPosition().phi() << " "
266  << ibest->recHit()->globalPosition().z() << endl;
267  } else {
268  cout << "TSB: found best measurement at invalid hit on det " << ibest->recHit()->geographicalId().rawId() << endl;
269  }
270  }
271  }
272  //
273  // keep old trajectorTempy
274  //
275  candidates.push_back(traj);
276 }
277 
278 vector<TrajectoryMeasurement>
280  const DetGroup& detGroup) const
281 {
282  vector<TM> result;
283  vector<TM> tmpResult;
284  //
285  // loop over all dets
286  //
287  if (theDbgFlg) cout << "TSB::redoMeasurements : nr. of measurements / group =";
288 
289  for (DetGroup::const_iterator idet=detGroup.begin();
290  idet!=detGroup.end(); ++idet) {
291  //
292  // ======== ask for measurements ==============
293  //B.M vector<TM> tmp = idet->det()->measurements(traj.lastMeasurement().updatedState(),
294  // theGeomPropagator,theEstimator);
295 
296  pair<bool, TrajectoryStateOnSurface> compat =
298  traj.lastMeasurement().updatedState(),
299  theGeomPropagator,theEstimator);
300 
301  vector<TM> tmp;
302  if(compat.first){
303  const MeasurementDet* mdet = theMeasurementTracker->idToDet(idet->det()->geographicalId());
304  tmp = mdet->fastMeasurements( compat.second, idet->trajectoryState(), theGeomPropagator, theEstimator);
305  }
306 
307  //perhaps could be enough just:
308  //vector<TM> tmp = mdet->fastMeasurements( idet->trajectoryState(),
309  // traj.lastMeasurement().updatedState(),
310  // theGeomPropagator, theEstimator);
311  // ==================================================
312 
313  if ( tmp.empty() ) continue;
314  //
315  // only collect valid RecHits
316  //
317  vector<TM>::iterator end = (tmp.back().recHit()->getType() != TrackingRecHit::missing ? tmp.end() : tmp.end()-1);
318  if (theDbgFlg) cout << " " << tmp.size();
319  tmpResult.insert(tmpResult.end(),tmp.begin(),end);
320  }
321  if (theDbgFlg) cout << endl;
322  //
323  // set layer in TM, because the Det cannot do it
324  //
325  for(vector<TM>::const_iterator tmpIt=tmpResult.begin();tmpIt!=tmpResult.end();++tmpIt){
326  if ( tmpIt->recHit()->isValid() )
327  result.push_back( TrajectoryMeasurement(tmpIt->predictedState(),tmpIt->recHit(),tmpIt->estimate(),&theLayer) );
328  }
329 
330  return result;
331 }
332 
333 void
335  const vector<TMG>& groups,
336  TempTrajectoryContainer& candidates) const
337 {
338  //
339  // first try to find an inactive hit with dets crossed by the prediction,
340  // then take any inactive hit
341  //
342  // loop over groups
343  for ( int iteration=0; iteration<2; iteration++ ) {
344  for ( vector<TMG>::const_iterator ig=groups.begin(); ig!=groups.end(); ++ig) {
345  // loop over measurements
346  const vector<TM>& measurements = ig->measurements();
347  for ( vector<TM>::const_reverse_iterator im=measurements.rbegin();
348  im!=measurements.rend(); ++im ) {
349  ConstRecHitPointer hit = im->recHit();
350  if ( hit->getType()==TrackingRecHit::valid ||
351  hit->getType()==TrackingRecHit::missing ) continue;
352  //
353  // check, if the extrapolation traverses the Det or
354  // if 2nd iteration
355  //
356  if ( hit->det() ) {
357  TSOS predState(im->predictedState());
358  if ( iteration>0 ||
359  (predState.isValid() &&
360  hit->det()->surface().bounds().inside(predState.localPosition())) ) {
361  // add the hit
362  /*TempTrajectory newTraj(traj);
363  updateTrajectory(newTraj,*im);
364  candidates.push_back(newTraj); // FIXME: avoid useless copy */
365  candidates.push_back(traj);
366  updateTrajectory(candidates.back(), *im);
367  if ( theDbgFlg ) cout << "TrajectorySegmentBuilder::updateWithInvalidHit "
368  << "added inactive hit" << endl;
369  return;
370  }
371  }
372  }
373  }
374  }
375  //
376  // No suitable inactive hit: add a missing one
377  //
378  bool found(false);
379  for ( int iteration=0; iteration<2; iteration++ ) {
380  //
381  // loop over groups
382  //
383  for ( vector<TMG>::const_iterator ig=groups.begin();
384  ig!=groups.end(); ++ig) {
385  const vector<TM>& measurements = ig->measurements();
386  for ( vector<TM>::const_reverse_iterator im=measurements.rbegin();
387  im!=measurements.rend(); ++im ) {
388  //
389  // only use invalid hits
390  //
391  ConstRecHitPointer hit = im->recHit();
392  if ( hit->isValid() ) continue;
393 
394  //
395  // check, if the extrapolation traverses the Det
396  //
397  TSOS predState(im->predictedState());
398  if(hit->det()){
399  if ( iteration>0 || (predState.isValid() &&
400  hit->det()->surface().bounds().inside(predState.localPosition())) ) {
401  // add invalid hit
402  /*TempTrajectory newTraj(traj);
403  updateTrajectory(newTraj,*im);
404  candidates.push_back(newTraj); // FIXME: avoid useless copy */
405  candidates.push_back(traj);
406  updateTrajectory(candidates.back(), *im);
407  found = true;
408  break;
409  }
410 
411  }else{
412  if ( iteration>0 || (predState.isValid() &&
413  im->layer()->surface().bounds().inside(predState.localPosition())) ){
414  // add invalid hit
415  /*TempTrajectory newTraj(traj);
416  updateTrajectory(newTraj,*im);
417  candidates.push_back(newTraj); // FIXME: avoid useless copy */
418  candidates.push_back(traj);
419  updateTrajectory(candidates.back(), *im);
420  found = true;
421  break;
422  }
423  }
424  }
425  if ( found ) break;
426  }
427  if ( theDbgFlg && !found ) cout << "TrajectorySegmentBuilder::updateWithInvalidHit: "
428  << " did not find invalid hit on 1st iteration" << endl;
429  if ( found ) break;
430  }
431  if ( !found ) {
432  if (theDbgFlg) cout << "TrajectorySegmentBuilder::updateWithInvalidHit: "
433  << " did not find invalid hit" << endl;
434  }
435 }
436 
437 vector<TrajectoryMeasurement>
438 TrajectorySegmentBuilder::unlockedMeasurements (const vector<TM>& measurements) const
439 {
440 // if ( !theLockHits ) return measurements;
441 
442  //========== B.M. to be ported later ===============
443  vector<TM> result;
444  result.reserve(measurements.size());
445 
446  //RecHitEqualByChannels recHitEqual(false,true);
447 
448  for ( vector<TM>::const_iterator im=measurements.begin();
449  im!=measurements.end(); ++im ) {
450  ConstRecHitPointer testHit = im->recHit();
451  if ( !testHit->isValid() ) continue;
452  bool found(false);
453  if ( theLockHits ) {
454  for ( ConstRecHitContainer::const_iterator ih=theLockedHits.begin();
455  ih!=theLockedHits.end(); ++ih ) {
456  if ( (*ih)->hit()->sharesInput(testHit->hit(), TrackingRecHit::all) ) {
457  found = true;
458  break;
459  }
460  }
461  }
462  if ( !found ) result.push_back(*im);
463  }
464  return result;
465  //=================================
466  //return measurements; // temporary solution before RecHitEqualByChannels is ported
467 }
468 
469 void
471 {
472  theLockedHits.push_back(measurement.recHit());
473 }
474 
475 
476 
477 // ================= B.M. to be ported later ===============================
478 void
479 TrajectorySegmentBuilder::cleanCandidates (vector<TempTrajectory>& candidates) const
480 {
481  //
482  // remove candidates which are subsets of others
483  // assumptions: no invalid hits and no duplicates
484  //
485  if ( candidates.size()<=1 ) return;
486  //RecHitEqualByChannels recHitEqual(false,true);
487  //
488  vector<TempTrajectory> sortedCandidates(candidates);
489  sort(sortedCandidates.begin(),sortedCandidates.end(),TrajectoryLessByFoundHits());
490 // cout << "SortedCandidates.foundHits";
491 // for ( vector<Trajectory>::iterator i1=sortedCandidates.begin();
492 // i1!=sortedCandidates.end(); i1++ )
493 // cout << " " << i1->foundHits();
494 // cout << endl;
495  //
496  for ( vector<TempTrajectory>::iterator i1=sortedCandidates.begin();
497  i1!=sortedCandidates.end()-1; ++i1 ) {
498  // get measurements of candidate to be checked
499  const TempTrajectory::DataContainer & measurements1 = i1->measurements();
500  for ( vector<TempTrajectory>::iterator i2=i1+1;
501  i2!=sortedCandidates.end(); ++i2 ) {
502  // no duplicates: two candidates of same size are different
503  if ( i2->foundHits()==i1->foundHits() ) continue;
504  // get measurements of "reference"
505  const TempTrajectory::DataContainer & measurements2 = i2->measurements();
506  //
507  // use the fact that TMs are ordered:
508  // start search in trajectory#1 from last hit match found
509  //
510  bool allFound(true);
511  TempTrajectory::DataContainer::const_iterator from2 = measurements2.rbegin(), im2end = measurements2.rend();
512  for ( TempTrajectory::DataContainer::const_iterator im1=measurements1.rbegin(),im1end = measurements1.rend();
513  im1!=im1end; --im1 ) {
514  // redundant protection - segments should not contain invalid RecHits
515  if ( !im1->recHit()->isValid() ) continue;
516  bool found(false);
518  im2!=im2end; --im2 ) {
519  // redundant protection - segments should not contain invalid RecHits
520  if ( !im2->recHit()->isValid() ) continue;
521  if ( im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::all) ) {
522  found = true;
523  from2 = im2; --from2;
524  break;
525  }
526  }
527  if ( !found ) {
528  allFound = false;
529  break;
530  }
531  }
532  if ( allFound ) i1->invalidate();
533  }
534  }
535 /*
536  candidates.clear();
537  for ( vector<Trajectory>::const_iterator i=sortedCandidates.begin();
538  i!=sortedCandidates.end(); i++ ) {
539  if ( i->isValid() ) candidates.push_back(*i);
540  }
541 */
542  candidates.erase(std::remove_if( candidates.begin(),candidates.end(),
543  std::not1(std::mem_fun_ref(&TempTrajectory::isValid))),
544  // boost::bind(&TempTrajectory::isValid,_1)),
545  candidates.end());
546 #ifdef DBG_TSB
547  cout << "TSB: cleanCandidates: reduced from " << sortedCandidates.size()
548  << " to " << candidates.size() << " candidates" << endl;
549 #endif
550  return;
551 }
552 
553 //==================================================
const_iterator rend() const
Definition: bqueue.h:142
void updateCandidatesWithBestHit(TempTrajectory &traj, const std::vector< TM > &measurements, TempTrajectoryContainer &candidates)
creation of a new candidate from a segment and the best hit out of a collection
bool empty() const
True if trajectory has no measurements.
void update(const LocalTrajectoryParameters &p, const Surface &aSurface, const MagneticField *field, const SurfaceSide side=SurfaceSideDefinition::atCenterOfSurface)
bool empty() const
Definition: bqueue.h:144
std::vector< TempTrajectory > TempTrajectoryContainer
void cleanCandidates(std::vector< TempTrajectory > &candidates) const
clean a set of candidates
void updateCandidates(TempTrajectory &traj, const std::vector< TM > &measurements, TempTrajectoryContainer &candidates)
creation of new candidates from a segment and a collection of hits
const DataContainer & measurements() const
std::vector< TrajectoryMeasurement > unlockedMeasurements(const std::vector< TM > &measurements) const
get list of unused hits
ConstRecHitPointer recHit() const
static std::pair< bool, TrajectoryStateOnSurface > isCompatible(const GeomDet *theDet, const TrajectoryStateOnSurface &ts, const Propagator &prop, const MeasurementEstimator &est)
TempTrajectoryContainer segments(const TSOS startingState)
new segments within layer
const TrajectoryMeasurement & lastMeasurement() const
tuple iteration
Definition: align_cfg.py:5
void lockMeasurement(const TM &measurement)
mark a hit as used
tuple result
Definition: query.py:137
std::vector< TrajectoryMeasurement > redoMeasurements(const TempTrajectory &traj, const DetGroup &detGroup) const
retrieve compatible hits from a DetGroup
virtual std::vector< TrajectoryMeasurement > fastMeasurements(const TrajectoryStateOnSurface &stateOnThisDet, const TrajectoryStateOnSurface &startingState, const Propagator &, const MeasurementEstimator &) const =0
const DetLayer * layer() const
TrajectoryStateOnSurface updatedState() const
const AlgebraicSymMatrix55 & matrix() const
#define end
Definition: vmac.h:38
TrajectoryStateOnSurface predictedState() const
const LocalTrajectoryError & localError() const
std::vector< TempTrajectory > addGroup(TempTrajectory &traj, std::vector< TrajectoryMeasurementGroup >::const_iterator begin, std::vector< TrajectoryMeasurementGroup >::const_iterator end)
bool isValid() const
iterator rbegin()
Definition: bqueue.h:140
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
#define begin
Definition: vmac.h:31
tuple cout
Definition: gather_cfg.py:41
size_type size() const
Definition: bqueue.h:143
void updateTrajectory(TempTrajectory &traj, const TM &tm) const
update of a trajectory with a hit
void push(const TrajectoryMeasurement &tm)
void updateWithInvalidHit(TempTrajectory &traj, const std::vector< TMG > &groups, TempTrajectoryContainer &candidates) const