CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoTracker/CkfPattern/src/TrajectorySegmentBuilder.cc

Go to the documentation of this file.
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 // #define DBG_TSB
00023 
00024 using namespace std;
00025 
00026 TrajectorySegmentBuilder::TempTrajectoryContainer
00027 TrajectorySegmentBuilder::segments (const TSOS startingState)
00028 {
00029   //
00030   // create empty trajectory
00031   //
00032   theLockedHits.clear();
00033   TempTrajectory startingTrajectory(theFullPropagator.propagationDirection());
00034   //
00035   // get measurement groups
00036   //
00037   vector<TMG> measGroups = 
00038     theLayerMeasurements->groupedMeasurements(theLayer,startingState,theFullPropagator,theEstimator);
00039     //B.M. theLayer.groupedMeasurements(startingState,theFullPropagator,theEstimator);
00040 
00041 #ifdef DBG_TSB
00042   cout << "TSB: number of measurement groups = " << measGroups.size() << endl;
00043   //  theDbgFlg = measGroups.size()>1;
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 //   if ( measGroups.size()>4 ) {
00082     cout << typeid(theLayer).name() << endl;
00083     cout << startingState.localError().matrix() << endl;
00084 //     for (vector<TMG>::const_iterator ig=measGroups.begin();
00085 //       ig!=measGroups.end(); ig++) {
00086 //       cout << "Nr. of measurements = " << ig->measurements().size() << endl;
00087 //       const DetGroup& dg = ig->detGroup();
00088 //       for ( DetGroup::const_iterator id=dg.begin();
00089 //          id!=dg.end(); id++ ) {
00090 //      GlobalPoint p(id->det()->position());
00091 //      GlobalVector v(id->det()->toGlobal(LocalVector(0.,0.,1.)));
00092 //      cout << p.perp() << " " << p.phi() << " " << p.z() << " ; "
00093 //           << v.phi() << " " << v.z() << endl;
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   // clean
00104   //
00105   //
00106   // add invalid hit - try to get first detector hit by the extrapolation
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 //     TrajectoryMeasurement tm(traj.lastMeasurement());
00124 //     if ( tm.updatedState().isValid() ) {
00125 //       if ( !hit.det().surface().bounds().inside(tm.updatedState().localPosition(),
00126 //                                              tm.updatedState().localError().positionError(),3.) ) {
00127 //      cout << "Incompatibility after update for det at " << hit.det().position() << ":" << endl;
00128 //      cout << tm.predictedState().localPosition() << " " 
00129 //           << tm.predictedState().localError().positionError() << endl;
00130 //      cout << hit.localPosition() << " " << hit.localPositionError() << endl;
00131 //      cout << tm.updatedState().localPosition() << " "
00132 //           << tm.updatedState().localError().positionError() << endl;
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     //std::cout << "TrajectorySegmentBuilder::addGroup" << " traj.empty()=" << traj.empty() << "EMPTY" << std::endl;
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                    //        << " nr. of candidates = " << candidates.size() 
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); // a good upper bound
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         //B.M. to be ported later
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   //std::cout << "TrajectorySegmentBuilder::addGroup" << 
00209   //             " traj.empty()=" << traj.empty() << 
00210   //             " end-begin=" << (end-begin)  <<
00211   //             " #updated=" << updatedTrajectories.size() << 
00212   //             " #result=" << ret.size() << std::endl;
00213   return ret;
00214 }
00215 
00216 void
00217 TrajectorySegmentBuilder::updateCandidates (TempTrajectory& traj,
00218                                             const vector<TM>& measurements,
00219                                             TempTrajectoryContainer& candidates)
00220 {
00221   //
00222   // generate updated candidates with all valid hits
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   // keep old trajectory
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 //     else if (!bestIsValid && (im->recHit()->getType() != TrackingRecHit::missing)) {
00254 //         if (ibest==measurements.end()) ibest = im;
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   // keep old trajectorTempy
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   // loop over all dets
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     // ======== ask for measurements ==============       
00293     //B.M vector<TM> tmp = idet->det()->measurements(traj.lastMeasurement().updatedState(),
00294     //                                         theGeomPropagator,theEstimator);
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     //perhaps could be enough just:
00308     //vector<TM> tmp = mdet->fastMeasurements( idet->trajectoryState(),
00309     //                                       traj.lastMeasurement().updatedState(),
00310     //                                       theGeomPropagator, theEstimator);
00311     // ==================================================
00312 
00313     if ( tmp.empty() )  continue;
00314     //
00315     // only collect valid RecHits
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   // set layer in TM, because the Det cannot do it
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   // first try to find an inactive hit with dets crossed by the prediction,
00340   //   then take any inactive hit
00341   //
00342   // loop over groups
00343   for ( int iteration=0; iteration<2; iteration++ ) {
00344     for ( vector<TMG>::const_iterator ig=groups.begin(); ig!=groups.end(); ++ig) {
00345       // loop over measurements
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         // check, if the extrapolation traverses the Det or
00354         // if 2nd iteration
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             // add the hit
00362             /*TempTrajectory newTraj(traj);
00363             updateTrajectory(newTraj,*im);
00364             candidates.push_back(newTraj);  // FIXME: avoid useless copy */
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   // No suitable inactive hit: add a missing one
00377   //
00378   bool found(false);
00379   for ( int iteration=0; iteration<2; iteration++ ) {
00380     //
00381     // loop over groups
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         // only use invalid hits
00390         //
00391         ConstRecHitPointer hit = im->recHit();
00392         if ( hit->isValid() )  continue;
00393 
00394         //
00395         // check, if the extrapolation traverses the Det
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             // add invalid hit
00402             /*TempTrajectory newTraj(traj);
00403             updateTrajectory(newTraj,*im);
00404             candidates.push_back(newTraj);  // FIXME: avoid useless copy */
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             // add invalid hit
00415             /*TempTrajectory newTraj(traj);
00416             updateTrajectory(newTraj,*im);
00417             candidates.push_back(newTraj);  // FIXME: avoid useless copy */
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 //   if ( !theLockHits )  return measurements;
00441 
00442   //========== B.M. to be ported later ===============
00443   vector<TM> result;
00444   result.reserve(measurements.size());
00445 
00446   //RecHitEqualByChannels recHitEqual(false,true);
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   //return measurements; // temporary solution before RecHitEqualByChannels is ported
00467 }
00468 
00469 void
00470 TrajectorySegmentBuilder::lockMeasurement (const TM& measurement)
00471 {
00472   theLockedHits.push_back(measurement.recHit());
00473 }
00474 
00475 
00476 
00477 // ================= B.M. to be ported later ===============================
00478 void
00479 TrajectorySegmentBuilder::cleanCandidates (vector<TempTrajectory>& candidates) const
00480 {
00481   //
00482   // remove candidates which are subsets of others
00483   // assumptions: no invalid hits and no duplicates
00484   //
00485   if ( candidates.size()<=1 )  return;
00486   //RecHitEqualByChannels recHitEqual(false,true);
00487   //
00488   vector<TempTrajectory> sortedCandidates(candidates);
00489   sort(sortedCandidates.begin(),sortedCandidates.end(),TrajectoryLessByFoundHits());
00490 //   cout << "SortedCandidates.foundHits";
00491 //   for ( vector<Trajectory>::iterator i1=sortedCandidates.begin();
00492 //      i1!=sortedCandidates.end(); i1++ ) 
00493 //     cout << " " << i1->foundHits();
00494 //   cout << endl;
00495   //
00496   for ( vector<TempTrajectory>::iterator i1=sortedCandidates.begin();
00497         i1!=sortedCandidates.end()-1; ++i1 ) {
00498     // get measurements of candidate to be checked
00499     const TempTrajectory::DataContainer & measurements1 = i1->measurements();
00500     for ( vector<TempTrajectory>::iterator i2=i1+1;
00501           i2!=sortedCandidates.end(); ++i2 ) {
00502       // no duplicates: two candidates of same size are different
00503       if ( i2->foundHits()==i1->foundHits() )  continue;
00504       // get measurements of "reference"
00505       const TempTrajectory::DataContainer & measurements2 = i2->measurements();
00506       //
00507       // use the fact that TMs are ordered:
00508       // start search in trajectory#1 from last hit match found
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         // redundant protection - segments should not contain invalid RecHits
00515         if ( !im1->recHit()->isValid() )  continue;
00516         bool found(false);
00517         for ( TempTrajectory::DataContainer::const_iterator im2=from2;
00518               im2!=im2end; --im2 ) {
00519           // redundant protection - segments should not contain invalid RecHits
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   candidates.clear();
00537   for ( vector<Trajectory>::const_iterator i=sortedCandidates.begin();
00538         i!=sortedCandidates.end(); i++ ) {
00539     if ( i->isValid() )  candidates.push_back(*i);
00540   }
00541 */
00542   candidates.erase(std::remove_if( candidates.begin(),candidates.end(),
00543                                    std::not1(std::mem_fun_ref(&TempTrajectory::isValid))),
00544  //                                boost::bind(&TempTrajectory::isValid,_1)), 
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 //==================================================