00001
00002 #include "RecoTracker/CkfPattern/interface/TrajectorySegmentBuilder.h"
00003
00004 #include "RecoTracker/CkfPattern/src/RecHitIsInvalid.h"
00005 #include "TrackingTools/PatternTools/interface/TempTrajectory.h"
00006
00007 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
00008 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00009 #include "TrackingTools/PatternTools/interface/MeasurementEstimator.h"
00010 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00011 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00012 #include "TrackingTools/MeasurementDet/interface/TrajectoryMeasurementGroup.h"
00013 #include "TrackingTools/DetLayers/interface/DetGroup.h"
00014 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00015 #include "RecoTracker/CkfPattern/src/TrajectoryLessByFoundHits.h"
00016 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
00017 #include "TrackingTools/DetLayers/interface/GeomDetCompatibilityChecker.h"
00018 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
00019
00020 #include <algorithm>
00021
00022
00023
00024 using namespace std;
00025
00026 TrajectorySegmentBuilder::TempTrajectoryContainer
00027 TrajectorySegmentBuilder::segments (const TSOS startingState)
00028 {
00029
00030
00031
00032 theLockedHits.clear();
00033 TempTrajectory startingTrajectory(theFullPropagator.propagationDirection());
00034
00035
00036
00037 vector<TMG> measGroups =
00038 theLayerMeasurements->groupedMeasurements(theLayer,startingState,theFullPropagator,theEstimator);
00039
00040
00041 #ifdef DBG_TSB
00042 cout << "TSB: number of measurement groups = " << measGroups.size() << endl;
00043
00044 theDbgFlg = true;
00045 #else
00046 theDbgFlg = false;
00047 #endif
00048
00049 if ( theDbgFlg ) {
00050 int ntot(1);
00051 for (vector<TMG>::const_iterator ig=measGroups.begin();
00052 ig!=measGroups.end(); ++ig) {
00053 int ngrp(0);
00054 const vector<TM>& measurements = ig->measurements();
00055 for ( vector<TM>::const_iterator im=measurements.begin();
00056 im!=measurements.end(); ++im ) {
00057 if ( im->recHit()->isValid() ) ngrp++;
00058 }
00059 cout << " " << ngrp;
00060 if ( ngrp>0 ) ntot *= ngrp;
00061 }
00062 cout << endl;
00063 cout << "TrajectorySegmentBuilder::partialTrajectories:: det ids & hit types / group" << endl;
00064 for (vector<TMG>::const_iterator ig=measGroups.begin();
00065 ig!=measGroups.end(); ++ig) {
00066 const vector<TM>& measurements = ig->measurements();
00067 for ( vector<TM>::const_iterator im=measurements.begin();
00068 im!=measurements.end(); ++im ) {
00069 if ( im!=measurements.begin() ) cout << " / ";
00070 if ( im->recHit()->det() )
00071 cout << im->recHit()->det()->geographicalId().rawId() << " "
00072 << im->recHit()->getType();
00073 else
00074 cout << "no det";
00075 }
00076 cout << endl;
00077 }
00078 }
00079
00080 #ifdef DBG_TSB
00081
00082 cout << typeid(theLayer).name() << endl;
00083 cout << startingState.localError().matrix() << endl;
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 #endif
00098
00099 TempTrajectoryContainer candidates =
00100 addGroup(startingTrajectory,measGroups.begin(),measGroups.end());
00101
00102 if (theDbgFlg) cout << "TSB: back with " << candidates.size() << " candidates" << endl;
00103
00104
00105
00106
00107
00108
00109 updateWithInvalidHit(startingTrajectory,measGroups,candidates);
00110 if (theDbgFlg) cout << "TSB: " << candidates.size() << " candidates after invalid hit" << endl;
00111 return candidates;
00112 }
00113
00114 void TrajectorySegmentBuilder::updateTrajectory (TempTrajectory& traj,
00115 const TM& tm) const
00116 {
00117 TSOS predictedState = tm.predictedState();
00118 ConstRecHitPointer hit = tm.recHit();
00119
00120 if ( hit->isValid()) {
00121 traj.push( TM( predictedState, theUpdator.update( predictedState, *hit),
00122 hit, tm.estimate(), tm.layer()));
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135 }
00136 else {
00137 traj.push( TM( predictedState, hit,0, tm.layer()));
00138 }
00139 }
00140
00141
00142 TrajectorySegmentBuilder::TempTrajectoryContainer
00143 TrajectorySegmentBuilder::addGroup (TempTrajectory& traj,
00144 vector<TMG>::const_iterator begin,
00145 vector<TMG>::const_iterator end)
00146 {
00147 vector<TempTrajectory> ret;
00148 if ( begin==end ) {
00149
00150 if (theDbgFlg) cout << "TSB::addGroup : no groups left" << endl;
00151 if ( !traj.empty() )
00152 ret.push_back(traj);
00153 return ret;
00154 }
00155
00156 if (theDbgFlg) cout << "TSB::addGroup : traj.size() = " << traj.measurements().size()
00157 << " first group at " << &(*begin)
00158
00159 << endl;
00160
00161
00162 TempTrajectoryContainer updatedTrajectories; updatedTrajectories.reserve(2);
00163 if ( traj.measurements().empty() ) {
00164 vector<TM> firstMeasurements = unlockedMeasurements(begin->measurements());
00165 if ( theBestHitOnly )
00166 updateCandidatesWithBestHit(traj,firstMeasurements,updatedTrajectories);
00167 else
00168 updateCandidates(traj,begin->measurements(),updatedTrajectories);
00169 if (theDbgFlg) cout << "TSB::addGroup : updating with first group - "
00170 << updatedTrajectories.size() << " trajectories" << endl;
00171 }
00172 else {
00173 if ( theBestHitOnly )
00174 updateCandidatesWithBestHit(traj,redoMeasurements(traj,begin->detGroup()),
00175 updatedTrajectories);
00176 else
00177 updateCandidates(traj,redoMeasurements(traj,begin->detGroup()),
00178 updatedTrajectories);
00179 if (theDbgFlg) cout << "TSB::addGroup : updating"
00180 << updatedTrajectories.size() << " trajectories" << endl;
00181 }
00182
00183 if (begin+1 != end) {
00184 ret.reserve(4);
00185 for ( TempTrajectoryContainer::iterator it=updatedTrajectories.begin();
00186 it!=updatedTrajectories.end(); ++it ) {
00187 if (theDbgFlg) cout << "TSB::addGroup : trying to extend candidate at "
00188 << &(*it) << " size " << it->measurements().size() << endl;
00189 vector<TempTrajectory> finalTrajectories = addGroup(*it,begin+1,end);
00190 if (theDbgFlg) cout << "TSB::addGroup : " << finalTrajectories.size()
00191 << " finalised candidates before cleaning" << endl;
00192
00193 cleanCandidates(finalTrajectories);
00194
00195 if (theDbgFlg) cout << "TSB::addGroup : got " << finalTrajectories.size()
00196 << " finalised candidates" << endl;
00197 ret.insert(ret.end(),finalTrajectories.begin(),
00198 finalTrajectories.end());
00199 }
00200 } else {
00201 ret.reserve(updatedTrajectories.size());
00202 for (TempTrajectoryContainer::iterator it=updatedTrajectories.begin();
00203 it!=updatedTrajectories.end(); ++it ) {
00204 if (!it->empty()) ret.push_back(*it);
00205 }
00206 }
00207
00208
00209
00210
00211
00212
00213 return ret;
00214 }
00215
00216 void
00217 TrajectorySegmentBuilder::updateCandidates (TempTrajectory& traj,
00218 const vector<TM>& measurements,
00219 TempTrajectoryContainer& candidates)
00220 {
00221
00222
00223
00224 for ( vector<TM>::const_iterator im=measurements.begin();
00225 im!=measurements.end(); ++im ) {
00226 if ( im->recHit()->isValid() ) {
00227 candidates.push_back(traj);
00228 updateTrajectory(candidates.back(),*im);
00229 if ( theLockHits ) lockMeasurement(*im);
00230 }
00231 }
00232
00233
00234
00235 candidates.push_back(traj);
00236 }
00237
00238 void
00239 TrajectorySegmentBuilder::updateCandidatesWithBestHit (TempTrajectory& traj,
00240 const vector<TM>& measurements,
00241 TempTrajectoryContainer& candidates)
00242 {
00243 vector<TM>::const_iterator ibest = measurements.end();
00244 bool bestIsValid = false;
00245 for ( vector<TM>::const_iterator im=measurements.begin();
00246 im!=measurements.end(); ++im ) {
00247 if ( im->recHit()->isValid() ) {
00248 if (!bestIsValid || (ibest==measurements.end()) || (im->estimate()<ibest->estimate()) ) {
00249 ibest = im;
00250 bestIsValid = true;
00251 }
00252 }
00253
00254
00255
00256 }
00257 if ( ibest!=measurements.end() ) {
00258 if ( theLockHits ) lockMeasurement(*ibest);
00259 candidates.push_back(traj);
00260 updateTrajectory(candidates.back(),*ibest);
00261 if ( theDbgFlg ) {
00262 if (bestIsValid) {
00263 cout << "TSB: found best measurement at "
00264 << ibest->recHit()->globalPosition().perp() << " "
00265 << ibest->recHit()->globalPosition().phi() << " "
00266 << ibest->recHit()->globalPosition().z() << endl;
00267 } else {
00268 cout << "TSB: found best measurement at invalid hit on det " << ibest->recHit()->geographicalId().rawId() << endl;
00269 }
00270 }
00271 }
00272
00273
00274
00275 candidates.push_back(traj);
00276 }
00277
00278 vector<TrajectoryMeasurement>
00279 TrajectorySegmentBuilder::redoMeasurements (const TempTrajectory& traj,
00280 const DetGroup& detGroup) const
00281 {
00282 vector<TM> result;
00283 vector<TM> tmpResult;
00284
00285
00286
00287 if (theDbgFlg) cout << "TSB::redoMeasurements : nr. of measurements / group =";
00288
00289 for (DetGroup::const_iterator idet=detGroup.begin();
00290 idet!=detGroup.end(); ++idet) {
00291
00292
00293
00294
00295
00296 pair<bool, TrajectoryStateOnSurface> compat =
00297 GeomDetCompatibilityChecker().isCompatible(idet->det(),
00298 traj.lastMeasurement().updatedState(),
00299 theGeomPropagator,theEstimator);
00300
00301 vector<TM> tmp;
00302 if(compat.first){
00303 const MeasurementDet* mdet = theMeasurementTracker->idToDet(idet->det()->geographicalId());
00304 tmp = mdet->fastMeasurements( compat.second, idet->trajectoryState(), theGeomPropagator, theEstimator);
00305 }
00306
00307
00308
00309
00310
00311
00312
00313 if ( tmp.empty() ) continue;
00314
00315
00316
00317 vector<TM>::iterator end = (tmp.back().recHit()->getType() != TrackingRecHit::missing ? tmp.end() : tmp.end()-1);
00318 if (theDbgFlg) cout << " " << tmp.size();
00319 tmpResult.insert(tmpResult.end(),tmp.begin(),end);
00320 }
00321 if (theDbgFlg) cout << endl;
00322
00323
00324
00325 for(vector<TM>::const_iterator tmpIt=tmpResult.begin();tmpIt!=tmpResult.end();++tmpIt){
00326 if ( tmpIt->recHit()->isValid() )
00327 result.push_back( TrajectoryMeasurement(tmpIt->predictedState(),tmpIt->recHit(),tmpIt->estimate(),&theLayer) );
00328 }
00329
00330 return result;
00331 }
00332
00333 void
00334 TrajectorySegmentBuilder::updateWithInvalidHit (TempTrajectory& traj,
00335 const vector<TMG>& groups,
00336 TempTrajectoryContainer& candidates) const
00337 {
00338
00339
00340
00341
00342
00343 for ( int iteration=0; iteration<2; iteration++ ) {
00344 for ( vector<TMG>::const_iterator ig=groups.begin(); ig!=groups.end(); ++ig) {
00345
00346 const vector<TM>& measurements = ig->measurements();
00347 for ( vector<TM>::const_reverse_iterator im=measurements.rbegin();
00348 im!=measurements.rend(); ++im ) {
00349 ConstRecHitPointer hit = im->recHit();
00350 if ( hit->getType()==TrackingRecHit::valid ||
00351 hit->getType()==TrackingRecHit::missing ) continue;
00352
00353
00354
00355
00356 if ( hit->det() ) {
00357 TSOS predState(im->predictedState());
00358 if ( iteration>0 ||
00359 (predState.isValid() &&
00360 hit->det()->surface().bounds().inside(predState.localPosition())) ) {
00361
00362
00363
00364
00365 candidates.push_back(traj);
00366 updateTrajectory(candidates.back(), *im);
00367 if ( theDbgFlg ) cout << "TrajectorySegmentBuilder::updateWithInvalidHit "
00368 << "added inactive hit" << endl;
00369 return;
00370 }
00371 }
00372 }
00373 }
00374 }
00375
00376
00377
00378 bool found(false);
00379 for ( int iteration=0; iteration<2; iteration++ ) {
00380
00381
00382
00383 for ( vector<TMG>::const_iterator ig=groups.begin();
00384 ig!=groups.end(); ++ig) {
00385 const vector<TM>& measurements = ig->measurements();
00386 for ( vector<TM>::const_reverse_iterator im=measurements.rbegin();
00387 im!=measurements.rend(); ++im ) {
00388
00389
00390
00391 ConstRecHitPointer hit = im->recHit();
00392 if ( hit->isValid() ) continue;
00393
00394
00395
00396
00397 TSOS predState(im->predictedState());
00398 if(hit->det()){
00399 if ( iteration>0 || (predState.isValid() &&
00400 hit->det()->surface().bounds().inside(predState.localPosition())) ) {
00401
00402
00403
00404
00405 candidates.push_back(traj);
00406 updateTrajectory(candidates.back(), *im);
00407 found = true;
00408 break;
00409 }
00410
00411 }else{
00412 if ( iteration>0 || (predState.isValid() &&
00413 im->layer()->surface().bounds().inside(predState.localPosition())) ){
00414
00415
00416
00417
00418 candidates.push_back(traj);
00419 updateTrajectory(candidates.back(), *im);
00420 found = true;
00421 break;
00422 }
00423 }
00424 }
00425 if ( found ) break;
00426 }
00427 if ( theDbgFlg && !found ) cout << "TrajectorySegmentBuilder::updateWithInvalidHit: "
00428 << " did not find invalid hit on 1st iteration" << endl;
00429 if ( found ) break;
00430 }
00431 if ( !found ) {
00432 if (theDbgFlg) cout << "TrajectorySegmentBuilder::updateWithInvalidHit: "
00433 << " did not find invalid hit" << endl;
00434 }
00435 }
00436
00437 vector<TrajectoryMeasurement>
00438 TrajectorySegmentBuilder::unlockedMeasurements (const vector<TM>& measurements) const
00439 {
00440
00441
00442
00443 vector<TM> result;
00444 result.reserve(measurements.size());
00445
00446
00447
00448 for ( vector<TM>::const_iterator im=measurements.begin();
00449 im!=measurements.end(); ++im ) {
00450 ConstRecHitPointer testHit = im->recHit();
00451 if ( !testHit->isValid() ) continue;
00452 bool found(false);
00453 if ( theLockHits ) {
00454 for ( ConstRecHitContainer::const_iterator ih=theLockedHits.begin();
00455 ih!=theLockedHits.end(); ++ih ) {
00456 if ( (*ih)->hit()->sharesInput(testHit->hit(), TrackingRecHit::all) ) {
00457 found = true;
00458 break;
00459 }
00460 }
00461 }
00462 if ( !found ) result.push_back(*im);
00463 }
00464 return result;
00465
00466
00467 }
00468
00469 void
00470 TrajectorySegmentBuilder::lockMeasurement (const TM& measurement)
00471 {
00472 theLockedHits.push_back(measurement.recHit());
00473 }
00474
00475
00476
00477
00478 void
00479 TrajectorySegmentBuilder::cleanCandidates (vector<TempTrajectory>& candidates) const
00480 {
00481
00482
00483
00484
00485 if ( candidates.size()<=1 ) return;
00486
00487
00488 vector<TempTrajectory> sortedCandidates(candidates);
00489 sort(sortedCandidates.begin(),sortedCandidates.end(),TrajectoryLessByFoundHits());
00490
00491
00492
00493
00494
00495
00496 for ( vector<TempTrajectory>::iterator i1=sortedCandidates.begin();
00497 i1!=sortedCandidates.end()-1; ++i1 ) {
00498
00499 const TempTrajectory::DataContainer & measurements1 = i1->measurements();
00500 for ( vector<TempTrajectory>::iterator i2=i1+1;
00501 i2!=sortedCandidates.end(); ++i2 ) {
00502
00503 if ( i2->foundHits()==i1->foundHits() ) continue;
00504
00505 const TempTrajectory::DataContainer & measurements2 = i2->measurements();
00506
00507
00508
00509
00510 bool allFound(true);
00511 TempTrajectory::DataContainer::const_iterator from2 = measurements2.rbegin(), im2end = measurements2.rend();
00512 for ( TempTrajectory::DataContainer::const_iterator im1=measurements1.rbegin(),im1end = measurements1.rend();
00513 im1!=im1end; --im1 ) {
00514
00515 if ( !im1->recHit()->isValid() ) continue;
00516 bool found(false);
00517 for ( TempTrajectory::DataContainer::const_iterator im2=from2;
00518 im2!=im2end; --im2 ) {
00519
00520 if ( !im2->recHit()->isValid() ) continue;
00521 if ( im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::all) ) {
00522 found = true;
00523 from2 = im2; --from2;
00524 break;
00525 }
00526 }
00527 if ( !found ) {
00528 allFound = false;
00529 break;
00530 }
00531 }
00532 if ( allFound ) i1->invalidate();
00533 }
00534 }
00535
00536
00537
00538
00539
00540
00541
00542 candidates.erase(std::remove_if( candidates.begin(),candidates.end(),
00543 std::not1(std::mem_fun_ref(&TempTrajectory::isValid))),
00544
00545 candidates.end());
00546 #ifdef DBG_TSB
00547 cout << "TSB: cleanCandidates: reduced from " << sortedCandidates.size()
00548 << " to " << candidates.size() << " candidates" << endl;
00549 #endif
00550 return;
00551 }
00552
00553