CkfPattern/src/CkfTrajectoryBuilder.cc

0001 #include "RecoTracker/CkfPattern/interface/CkfTrajectoryBuilder.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ParameterSet/interface/PluginDescription.h"
0005
0006 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0007
0008 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0009 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
0010 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0011
0012 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0013 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0014 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
0015 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0016 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0017 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0018
0019 #include "RecoTracker/CkfPattern/interface/TrajCandLess.h"
0020
0021 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0022
0023 #include "RecoTracker/CkfPattern/interface/IntermediateTrajectoryCleaner.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0025 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0026 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilter.h"
0027 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilterFactory.h"
0028
0029 using namespace std;
0030
0031 CkfTrajectoryBuilder::CkfTrajectoryBuilder(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0032     : CkfTrajectoryBuilder(conf,
0033                            iC,
0034                            BaseCkfTrajectoryBuilder::createTrajectoryFilter(
0035                                conf.getParameter<edm::ParameterSet>("trajectoryFilter"), iC)) {}
0036
0037 CkfTrajectoryBuilder::CkfTrajectoryBuilder(const edm::ParameterSet& conf,
0038                                            edm::ConsumesCollector iC,
0039                                            std::unique_ptr<TrajectoryFilter> filter)
0040     : BaseCkfTrajectoryBuilder(conf, iC, std::move(filter)) {
0041   theMaxCand = conf.getParameter<int>("maxCand");
0042   theLostHitPenalty = conf.getParameter<double>("lostHitPenalty");
0043   theIntermediateCleaning = conf.getParameter<bool>("intermediateCleaning");
0044   theAlwaysUseInvalidHits = conf.getParameter<bool>("alwaysUseInvalidHits");
0045   /*
0046     theSharedSeedCheck = conf.getParameter<bool>("SharedSeedCheck");
0047     std::stringstream ss;
0048     ss<<"CkfTrajectoryBuilder_"<<conf.getParameter<std::string>("ComponentName")<<"_"<<this;
0049     theUniqueName = ss.str();
0050     LogDebug("CkfPattern")<<"my unique name is: "<<theUniqueName;
0051   */
0052 }
0053
0054 void CkfTrajectoryBuilder::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0055   BaseCkfTrajectoryBuilder::fillPSetDescription(iDesc);
0056   iDesc.add<int>("maxCand", 5);
0057   iDesc.add<double>("lostHitPenalty", 30.);
0058   iDesc.add<bool>("intermediateCleaning", true);
0059   iDesc.add<bool>("alwaysUseInvalidHits", true);
0060
0061   edm::ParameterSetDescription psdTF;
0062   psdTF.addNode(edm::PluginDescription<TrajectoryFilterFactory>("ComponentType", true));
0063   iDesc.add<edm::ParameterSetDescription>("trajectoryFilter", psdTF);
0064 }
0065
0066 /*
0067   void CkfTrajectoryBuilder::setEvent(const edm::Event& event) const
0068   {
0069   theMeasurementTracker->update(event);
0070   }
0071 */
0072
0073 void CkfTrajectoryBuilder::setEvent_(const edm::Event& event, const edm::EventSetup& iSetup) {}
0074
0075 CkfTrajectoryBuilder::TrajectoryContainer CkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed) const {
0076   TrajectoryContainer result;
0077   result.reserve(5);
0078   trajectories(seed, result);
0079   return result;
0080 }
0081
0082 /*
0083   void CkfTrajectoryBuilder::rememberSeedAndTrajectories(const TrajectorySeed& seed,
0084   CkfTrajectoryBuilder::TrajectoryContainer &result) const
0085   {
0086
0087   //result ----> theCachedTrajectories
0088   //every first iteration on event. forget about everything that happened before
0089   if (edm::Service<UpdaterService>()->checkOnce(theUniqueName))
0090   theCachedTrajectories.clear();
0091
0092   if (seed.nHits()==0) return;
0093
0094   //then remember those trajectories
0095   for (TrajectoryContainer::iterator traj=result.begin();
0096   traj!=result.end(); ++traj) {
0097   theCachedTrajectories.insert(std::make_pair(seed.recHits().first->geographicalId(),*traj));
0098   }
0099   }
0100
0101   bool CkfTrajectoryBuilder::sharedSeed(const TrajectorySeed& s1,const TrajectorySeed& s2) const{
0102   //quit right away on nH=0
0103   if (s1.nHits()==0 || s2.nHits()==0) return false;
0104   //quit right away if not the same number of hits
0105   if (s1.nHits()!=s2.nHits()) return false;
0106   TrajectorySeed::range r1=s1.recHits();
0107   TrajectorySeed::range r2=s2.recHits();
0108   TrajectorySeed::const_iterator i1,i2;
0109   TrajectorySeed::const_iterator & i1_e=r1.second,&i2_e=r2.second;
0110   TrajectorySeed::const_iterator & i1_b=r1.first,&i2_b=r2.first;
0111   //quit right away if first detId does not match. front exist because of ==0 ->quit test
0112   if(i1_b->geographicalId() != i2_b->geographicalId()) return false;
0113   //then check hit by hit if they are the same
0114   for (i1=i1_b,i2=i2_b;i1!=i1_e,i2!=i2_e;++i1,++i2){
0115   if (!i1->sharesInput(&(*i2),TrackingRecHit::all)) return false;
0116   }
0117   return true;
0118   }
0119   bool CkfTrajectoryBuilder::seedAlreadyUsed(const TrajectorySeed& seed,
0120   CkfTrajectoryBuilder::TempTrajectoryContainer &candidates) const
0121   {
0122   //theCachedTrajectories ---> candidates
0123   if (seed.nHits()==0) return false;
0124   bool answer=false;
0125   pair<SharedTrajectory::const_iterator, SharedTrajectory::const_iterator> range =
0126   theCachedTrajectories.equal_range(seed.recHits().first->geographicalId());
0127   SharedTrajectory::const_iterator trajP;
0128   for (trajP = range.first; trajP!=range.second;++trajP){
0129   //check whether seeds are identical
0130   if (sharedSeed(trajP->second.seed(),seed)){
0131   candidates.push_back(trajP->second);
0132   answer=true;
0133   }//already existing trajectory shares the seed.
0134   }//loop already made trajectories
0135
0136   return answer;
0137   }
0138 */
0139
0140 void CkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed,
0141                                         CkfTrajectoryBuilder::TrajectoryContainer& result) const {
0142   // analyseSeed( seed);
0143   /*
0144     if (theSharedSeedCheck){
0145     TempTrajectoryContainer candidates;
0146     if (seedAlreadyUsed(seed,candidates))
0147     {
0148     //start with those candidates already made before
0149     limitedCandidates(candidates,result);
0150     //and quit
0151     return;
0152     }
0153     }
0154   */
0155
0156   unsigned int tmp;
0157   buildTrajectories(seed, result, tmp, nullptr);
0158 }
0159
0160 TempTrajectory CkfTrajectoryBuilder::buildTrajectories(const TrajectorySeed& seed,
0161                                                        TrajectoryContainer& result,
0162                                                        unsigned int& nCandPerSeed,
0163                                                        const TrajectoryFilter*) const {
0164   if (theMeasurementTracker == nullptr) {
0165     throw cms::Exception("LogicError")
0166         << "Asking to create trajectories to an un-initialized CkfTrajectoryBuilder.\nYou have to call clone(const "
0167            "MeasurementTrackerEvent *data) and then call trajectories on it instead.\n";
0168   }
0169
0170   TempTrajectory startingTraj = createStartingTrajectory(seed);
0171
0172   /// limitedCandidates( startingTraj, regionalCondition, result);
0173   /// FIXME: restore regionalCondition
0174   nCandPerSeed = limitedCandidates(seed, startingTraj, result);
0175
0176   return startingTraj;
0177
0178   /*
0179   //and remember what you just did
0180   if (theSharedSeedCheck)  rememberSeedAndTrajectories(seed,result);
0181   */
0182
0183   // analyseResult(result);
0184 }
0185
0186 unsigned int CkfTrajectoryBuilder::limitedCandidates(const TrajectorySeed& seed,
0187                                                      TempTrajectory& startingTraj,
0188                                                      TrajectoryContainer& result) const {
0189   TempTrajectoryContainer candidates;
0190   candidates.push_back(startingTraj);
0191   std::shared_ptr<const TrajectorySeed> sharedSeed(new TrajectorySeed(seed));
0192   return limitedCandidates(sharedSeed, candidates, result);
0193 }
0194
0195 unsigned int CkfTrajectoryBuilder::limitedCandidates(const std::shared_ptr<const TrajectorySeed>& sharedSeed,
0196                                                      TempTrajectoryContainer& candidates,
0197                                                      TrajectoryContainer& result) const {
0198   unsigned int nIter = 1;
0199   unsigned int nCands = 0;  // ignore startingTraj
0200   unsigned int prevNewCandSize = 0;
0201   TempTrajectoryContainer newCand;  // = TrajectoryContainer();
0202   newCand.reserve(2 * theMaxCand);
0203
0204   auto trajCandLess = [&](TempTrajectory const& a, TempTrajectory const& b) {
0205     return (a.chiSquared() + a.lostHits() * theLostHitPenalty) < (b.chiSquared() + b.lostHits() * theLostHitPenalty);
0206   };
0207
0208   while (!candidates.empty()) {
0209     newCand.clear();
0210     for (auto traj = candidates.begin(); traj != candidates.end(); traj++) {
0211       std::vector<TM> meas;
0212       findCompatibleMeasurements(*sharedSeed, *traj, meas);
0213
0214       // --- method for debugging
0215       if (!analyzeMeasurementsDebugger(
0216               *traj, meas, theMeasurementTracker, forwardPropagator(*sharedSeed), theEstimator, theTTRHBuilder))
0217         return nCands;
0218       // ---
0219
0220       if (meas.empty()) {
0221         addToResult(sharedSeed, *traj, result);
0222       } else {
0223         std::vector<TM>::const_iterator last;
0224         if (theAlwaysUseInvalidHits)
0225           last = meas.end();
0226         else {
0227           if (meas.front().recHit()->isValid()) {
0228             last = find_if(meas.begin(), meas.end(), [](auto const& meas) { return !meas.recHit()->isValid(); });
0229           } else
0230             last = meas.end();
0231         }
0232
0233         for (auto itm = meas.begin(); itm != last; itm++) {
0234           TempTrajectory newTraj = *traj;
0235           updateTrajectory(newTraj, std::move(*itm));
0236
0237           if (toBeContinued(newTraj)) {
0238             newCand.push_back(std::move(newTraj));
0239             std::push_heap(newCand.begin(), newCand.end(), trajCandLess);
0240           } else {
0241             addToResult(sharedSeed, newTraj, result);
0242             //// don't know yet
0243           }
0244         }
0245       }
0246
0247       // account only new candidates, i.e.
0248       // - 1 candidate -> 1 candidate, don't increase count
0249       // - 1 candidate -> 2 candidates, increase count by 1
0250       nCands += newCand.size() - prevNewCandSize;
0251       prevNewCandSize = newCand.size();
0252
0253       /*
0254       auto trajVal = [&](TempTrajectory const & a) {
0255         return  a.chiSquared() + a.lostHits()*theLostHitPenalty;
0256       };
0257
0258       // safe (stable?) logig: always sort, kill exceeding only if worse than last to keep
0259       // if ((int)newCand.size() > theMaxCand) std::cout << "TrajVal " << theMaxCand  << ' ' << newCand.size() << ' ' <<  trajVal(newCand.front());
0260       int toCut = int(newCand.size()) - int(theMaxCand);
0261       if (toCut>0) {
0262         // move largest "toCut" to the end
0263         for (int i=0; i<toCut; ++i)
0264           std::pop_heap(newCand.begin(),newCand.end()-i,trajCandLess);
0265         auto fval = trajVal(newCand.front());
0266         // remove till equal to highest to keep
0267         for (int i=0; i<toCut; ++i) {
0268            if (fval==trajVal(newCand.back())) break;
0269            newCand.pop_back();
0270         }
0271     //assert((int)newCand.size() >= theMaxCand);
0272     //std::cout << "; " << newCand.size() << ' ' << trajVal(newCand.front())  << " " << trajVal(newCand.back());
0273
0274     // std::make_heap(newCand.begin(),newCand.end(),trajCandLess);
0275         // push_heap again the one left
0276         for (auto iter = newCand.begin()+theMaxCand+1; iter<=newCand.end(); ++iter  )
0277       std::push_heap(newCand.begin(),iter,trajCandLess);
0278
0279     // std::cout << "; " << newCand.size() << ' ' << trajVal(newCand.front())  << " " << trajVal(newCand.back()) << std::endl;
0280       }
0281
0282       */
0283
0284       // intermedeate login: always sort,  kill all exceeding
0285       while ((int)newCand.size() > theMaxCand) {
0286         std::pop_heap(newCand.begin(), newCand.end(), trajCandLess);
0287         // if ((int)newCand.size() == theMaxCand+1) std::cout << " " << trajVal(newCand.front())  << " " << trajVal(newCand.back()) << std::endl;
0288         newCand.pop_back();
0289       }
0290
0291       /*
0292       //   original logic: sort only if > theMaxCand, kill all exceeding
0293       if ((int)newCand.size() > theMaxCand) {
0294     std::sort( newCand.begin(), newCand.end(), TrajCandLess<TempTrajectory>(theLostHitPenalty));
0295     // std::partial_sort( newCand.begin(), newCand.begin()+theMaxCand, newCand.end(), TrajCandLess<TempTrajectory>(theLostHitPenalty));
0296     std::cout << "TrajVal " << theMaxCand  << ' ' << newCand.size() << ' '
0297     << trajVal(newCand.back()) << ' ' << trajVal(newCand[theMaxCand-1]) << ' ' << trajVal(newCand[theMaxCand])  << std::endl;
0298     newCand.resize(theMaxCand);
0299       }
0300       */
0301
0302     }  // end loop on candidates
0303
0304     std::sort_heap(newCand.begin(), newCand.end(), trajCandLess);
0305     if (theIntermediateCleaning)
0306       IntermediateTrajectoryCleaner::clean(newCand);
0307
0308     candidates.swap(newCand);
0309
0310     LogDebug("CkfPattern") << result.size() << " candidates after " << nIter++ << " CKF iteration: \n"
0311                            << PrintoutHelper::dumpCandidates(result) << "\n " << candidates.size()
0312                            << " running candidates are: \n"
0313                            << PrintoutHelper::dumpCandidates(candidates);
0314   }
0315   return nCands;
0316 }
0317
0318 void CkfTrajectoryBuilder::updateTrajectory(TempTrajectory& traj, TM&& tm) const {
0319   auto&& predictedState = tm.predictedState();
0320   auto&& hit = tm.recHit();
0321   if (hit->isValid()) {
0322     auto&& upState = theUpdator->update(predictedState, *hit);
0323     traj.emplace(predictedState, std::move(upState), hit, tm.estimate(), tm.layer());
0324   } else {
0325     traj.emplace(predictedState, hit, 0, tm.layer());
0326   }
0327 }
0328
0329 void CkfTrajectoryBuilder::findCompatibleMeasurements(const TrajectorySeed& seed,
0330                                                       const TempTrajectory& traj,
0331                                                       std::vector<TrajectoryMeasurement>& result) const {
0332   int invalidHits = 0;
0333   //Use findStateAndLayers which handles the hitless seed use case
0334   std::pair<TSOS, std::vector<const DetLayer*> >&& stateAndLayers = findStateAndLayers(seed, traj);
0335   if (stateAndLayers.second.empty())
0336     return;
0337
0338   auto layerBegin = stateAndLayers.second.begin();
0339   auto layerEnd = stateAndLayers.second.end();
0340   LogDebug("CkfPattern") << "looping on " << stateAndLayers.second.size() << " layers.";
0341   const Propagator* fwdPropagator = forwardPropagator(seed);
0342   for (auto il = layerBegin; il != layerEnd; il++) {
0343     LogDebug("CkfPattern") << "looping on a layer in findCompatibleMeasurements.\n last layer: " << traj.lastLayer()
0344                            << " current layer: " << (*il);
0345
0346     TSOS stateToUse = stateAndLayers.first;
0347     //Added protection before asking for the lastLayer on the trajectory
0348     if UNLIKELY (!traj.empty() && (*il) == traj.lastLayer()) {
0349       LogDebug("CkfPattern") << " self propagating in findCompatibleMeasurements.\n from: \n" << stateToUse;
0350       //self navigation case
0351       // go to a middle point first
0352       TransverseImpactPointExtrapolator middle;
0353       GlobalPoint center(0, 0, 0);
0354       stateToUse = middle.extrapolate(stateToUse, center, *fwdPropagator);
0355
0356       if (!stateToUse.isValid())
0357         continue;
0358       LogDebug("CkfPattern") << "to: " << stateToUse;
0359     }
0360
0361     LayerMeasurements layerMeasurements(theMeasurementTracker->measurementTracker(), *theMeasurementTracker);
0362     std::vector<TrajectoryMeasurement>&& tmp =
0363         layerMeasurements.measurements((**il), stateToUse, *fwdPropagator, *theEstimator);
0364
0365     if (!tmp.empty()) {
0366       if (result.empty())
0367         result.swap(tmp);
0368       else {
0369         // keep one dummy TM at the end, skip the others
0370         result.insert(
0371             result.end() - invalidHits, std::make_move_iterator(tmp.begin()), std::make_move_iterator(tmp.end()));
0372       }
0373       invalidHits++;
0374     }
0375   }
0376
0377   // sort the final result, keep dummy measurements at the end
0378   if (result.size() > 1) {
0379     std::sort(result.begin(), result.end() - invalidHits, TrajMeasLessEstim());
0380   }
0381
0382   LogDebug("CkfPattern") << "starting from:\n"
0383                          << "x: " << stateAndLayers.first.globalPosition() << "\n"
0384                          << "p: " << stateAndLayers.first.globalMomentum() << "\n"
0385                          << PrintoutHelper::dumpMeasurements(result);
0386
0387 #ifdef DEBUG_INVALID
0388   bool afterInvalid = false;
0389   for (vector<TM>::const_iterator i = result.begin(); i != result.end(); i++) {
0390     if (!i->recHit().isValid())
0391       afterInvalid = true;
0392     if (afterInvalid && i->recHit().isValid()) {
0393       edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: valid hit after invalid!";
0394     }
0395   }
0396 #endif
0397
0398   //analyseMeasurements( result, traj);
0399 }