CkfPattern/plugins/GroupedCkfTrajectoryBuilder.cc

0001 #include <algorithm>
0002 #include <array>
0003
0004 #include "GroupedCkfTrajectoryBuilder.h"
0005 #include "TrajectorySegmentBuilder.h"
0006
0007 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0008 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0009
0010 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0011 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0012 #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
0013 #include "GroupedTrajCandLess.h"
0014 #include "TrackingTools/TrajectoryFiltering/interface/RegionalTrajectoryFilter.h"
0015 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilterFactory.h"
0016 #include "TrackingTools/PatternTools/interface/TempTrajectory.h"
0017 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0018 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0019 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0020 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0021 #include "TrackingTools/DetLayers/interface/DetGroup.h"
0022 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0023 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0024 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0025 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0026
0027 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
0028 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0029 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0030 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0031 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 #include "FWCore/ParameterSet/interface/PluginDescription.h"
0034 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0035 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0036
0037 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0038
0039 // only included for RecHit comparison operator:
0040 #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h"
0041
0042 #include "TrackingTools/TrajectoryCleaning/interface/FastTrajectoryCleaner.h"
0043
0044 // for looper reconstruction
0045 #include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
0046 #include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
0047 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0048
0049 #include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
0050
0051 // #define STAT_TSB
0052
0053 namespace {
0054 #ifdef STAT_TSB
0055   struct StatCount {
0056     long long totSeed;
0057     long long totTraj;
0058     long long totRebuilt;
0059     long long totInvCand;
0060     void zero() { totSeed = totTraj = totRebuilt = totInvCand = 0; }
0061     void traj(long long t) { totTraj += t; }
0062     void seed() { ++totSeed; }
0063     void rebuilt(long long t) { totRebuilt += t; }
0064     void invalid() { ++totInvCand; }
0065     void print() const {
0066       std::cout << "GroupedCkfTrajectoryBuilder stat\nSeed/Traj/Rebuilt " << totSeed << '/' << totTraj << '/'
0067                 << totRebuilt << std::endl;
0068     }
0069     StatCount() { zero(); }
0070     ~StatCount() { print(); }
0071   };
0072   StatCount statCount;
0073
0074 #else
0075   struct StatCount {
0076     void traj(long long) {}
0077     void seed() {}
0078     void rebuilt(long long) {}
0079     void invalid() {}
0080   };
0081   CMS_THREAD_SAFE StatCount statCount;
0082 #endif
0083
0084 }  // namespace
0085
0086 using namespace std;
0087
0088 //#define DBG2_GCTB
0089
0090 //#define STANDARD_INTERMEDIARYCLEAN
0091
0092 #ifdef STANDARD_INTERMEDIARYCLEAN
0093 #include "RecoTracker/CkfPattern/interface/IntermediateTrajectoryCleaner.h"
0094 #endif
0095
0096 /* ====== B.M. to be ported layer ===========
0097 #ifdef DBG_GCTB
0098 #include "RecoTracker/CkfPattern/src/ShowCand.h"
0099 #endif
0100 // #define DBG2_GCTB
0101 #ifdef DBG2_GCTB
0102 #include "RecoTracker/CkfPattern/src/SimIdPrinter.h"
0103 #include "Tracker/TkDebugTools/interface/LayerFinderByDet.h"
0104 #include "Tracker/TkLayout/interface/TkLayerName.h"
0105 #endif
0106 =================================== */
0107
0108 namespace {
0109
0110   thread_local GroupedCkfTrajectoryBuilder::TempTrajectoryContainer work_;
0111
0112 }
0113
0114 GroupedCkfTrajectoryBuilder::GroupedCkfTrajectoryBuilder(const edm::ParameterSet& conf, edm::ConsumesCollector& iC)
0115     : BaseCkfTrajectoryBuilder(conf,
0116                                iC,
0117                                BaseCkfTrajectoryBuilder::createTrajectoryFilter(
0118                                    conf.getParameter<edm::ParameterSet>("trajectoryFilter"), iC),
0119                                conf.getParameter<bool>("useSameTrajFilter")
0120                                    ? BaseCkfTrajectoryBuilder::createTrajectoryFilter(
0121                                          conf.getParameter<edm::ParameterSet>("trajectoryFilter"), iC)
0122                                    : BaseCkfTrajectoryBuilder::createTrajectoryFilter(
0123                                          conf.getParameter<edm::ParameterSet>("inOutTrajectoryFilter"), iC)) {
0124   // fill data members from parameters (eventually data members could be dropped)
0125   //
0126   theMaxCand = conf.getParameter<int>("maxCand");
0127   theLostHitPenalty = conf.getParameter<double>("lostHitPenalty");
0128   theFoundHitBonus = conf.getParameter<double>("foundHitBonus");
0129   theIntermediateCleaning = conf.getParameter<bool>("intermediateCleaning");
0130   theAlwaysUseInvalid = conf.getParameter<bool>("alwaysUseInvalidHits");
0131   theLockHits = conf.getParameter<bool>("lockHits");
0132   theBestHitOnly = conf.getParameter<bool>("bestHitOnly");
0133   theMinNrOf2dHitsForRebuild = 2;
0134   theRequireSeedHitsInRebuild = conf.getParameter<bool>("requireSeedHitsInRebuild");
0135   theKeepOriginalIfRebuildFails = conf.getParameter<bool>("keepOriginalIfRebuildFails");
0136   theMinNrOfHitsForRebuild = max(0, conf.getParameter<int>("minNrOfHitsForRebuild"));
0137   maxPt2ForLooperReconstruction = conf.getParameter<double>("maxPtForLooperReconstruction");
0138   maxPt2ForLooperReconstruction *= maxPt2ForLooperReconstruction;
0139   maxDPhiForLooperReconstruction = conf.getParameter<double>("maxDPhiForLooperReconstruction");
0140
0141   /* ======= B.M. to be ported layer ===========
0142   bool setOK = thePropagator->setMaxDirectionChange(1.6);
0143   if (!setOK)
0144     cout  << "GroupedCkfTrajectoryBuilder WARNING: "
0145       << "propagator does not support setMaxDirectionChange"
0146       << endl;
0147   //   addStopCondition(theMinPtStopCondition);
0148
0149   theConfigurableCondition = createAlgo<TrajectoryFilter>(componentConfig("StopCondition"));
0150   ===================================== */
0151 }
0152
0153 void GroupedCkfTrajectoryBuilder::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0154   BaseCkfTrajectoryBuilder::fillPSetDescription(iDesc);
0155
0156   iDesc.add<bool>("useSameTrajFilter", true);
0157   iDesc.add<int>("maxCand", 5);
0158   iDesc.add<double>("lostHitPenalty", 30.);
0159   iDesc.add<double>("foundHitBonus", 10.);
0160   iDesc.add<bool>("intermediateCleaning", true);
0161   iDesc.add<bool>("alwaysUseInvalidHits", true);
0162   iDesc.add<bool>("lockHits", true);
0163   iDesc.add<bool>("bestHitOnly", true);
0164   iDesc.add<bool>("requireSeedHitsInRebuild", true);
0165   iDesc.add<bool>("keepOriginalIfRebuildFails", false);
0166   iDesc.add<int>("minNrOfHitsForRebuild", 5);
0167   iDesc.add<double>("maxPtForLooperReconstruction", 0.);
0168   iDesc.add<double>("maxDPhiForLooperReconstruction", 2.0);
0169
0170   edm::ParameterSetDescription psdTJ1;
0171   psdTJ1.addNode(edm::PluginDescription<TrajectoryFilterFactory>("ComponentType", true));
0172   iDesc.add<edm::ParameterSetDescription>("trajectoryFilter", psdTJ1);
0173
0174   edm::ParameterSetDescription psdTJ2;
0175   psdTJ2.addNode(edm::PluginDescription<TrajectoryFilterFactory>("ComponentType", true));
0176   iDesc.add<edm::ParameterSetDescription>("inOutTrajectoryFilter", psdTJ2);
0177 }
0178
0179 /*
0180   void GroupedCkfTrajectoryBuilder::setEvent(const edm::Event& event) const
0181   {
0182   theMeasurementTracker->update(event);
0183 }
0184 */
0185
0186 void GroupedCkfTrajectoryBuilder::setEvent_(const edm::Event& event, const edm::EventSetup& iSetup) {}
0187
0188 GroupedCkfTrajectoryBuilder::TrajectoryContainer GroupedCkfTrajectoryBuilder::trajectories(
0189     const TrajectorySeed& seed) const {
0190   TrajectoryContainer ret;
0191   ret.reserve(10);
0192   unsigned int tmp;
0193   buildTrajectories(seed, ret, tmp, nullptr);
0194   return ret;
0195 }
0196
0197 GroupedCkfTrajectoryBuilder::TrajectoryContainer GroupedCkfTrajectoryBuilder::trajectories(
0198     const TrajectorySeed& seed, const TrackingRegion& region) const {
0199   TrajectoryContainer ret;
0200   ret.reserve(10);
0201   unsigned int tmp;
0202   RegionalTrajectoryFilter regionalCondition(region);
0203   buildTrajectories(seed, ret, tmp, &regionalCondition);
0204   return ret;
0205 }
0206
0207 void GroupedCkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed,
0208                                                GroupedCkfTrajectoryBuilder::TrajectoryContainer& ret) const {
0209   unsigned int tmp;
0210   buildTrajectories(seed, ret, tmp, nullptr);
0211 }
0212
0213 void GroupedCkfTrajectoryBuilder::trajectories(const TrajectorySeed& seed,
0214                                                GroupedCkfTrajectoryBuilder::TrajectoryContainer& ret,
0215                                                const TrackingRegion& region) const {
0216   RegionalTrajectoryFilter regionalCondition(region);
0217   unsigned int tmp;
0218   buildTrajectories(seed, ret, tmp, &regionalCondition);
0219 }
0220
0221 void GroupedCkfTrajectoryBuilder::rebuildSeedingRegion(const TrajectorySeed& seed, TrajectoryContainer& result) const {
0222   rebuildTrajectories(seed, result);
0223 }
0224
0225 void GroupedCkfTrajectoryBuilder::rebuildTrajectories(const TrajectorySeed& seed, TrajectoryContainer& result) const {
0226   TempTrajectoryContainer work;
0227
0228   TrajectoryContainer final;
0229
0230   // better the seed to be always the same...
0231   std::shared_ptr<const TrajectorySeed> sharedSeed;
0232   if (result.empty())
0233     sharedSeed.reset(new TrajectorySeed(seed));
0234   else
0235     sharedSeed = result.front().sharedSeed();
0236
0237   work.reserve(result.size());
0238   for (auto&& traj : result)
0239     if (traj.isValid())
0240       work.emplace_back(std::move(traj));
0241
0242   rebuildSeedingRegion(seed, work);
0243
0244   // we clean here now
0245   FastTrajectoryCleaner cleaner(theFoundHitBonus, theLostHitPenalty, false);
0246   cleaner.clean(work);
0247
0248   for (auto const& it : work)
0249     if (it.isValid()) {
0250       final.push_back(it.toTrajectory());
0251       final.back().setSharedSeed(sharedSeed);
0252     }
0253
0254   result.swap(final);
0255
0256   statCount.rebuilt(result.size());
0257 }
0258
0259 void GroupedCkfTrajectoryBuilder::buildTrajectories(const TrajectorySeed& seed,
0260                                                     GroupedCkfTrajectoryBuilder::TrajectoryContainer& result,
0261                                                     unsigned int& nCandPerSeed,
0262                                                     const TrajectoryFilter* regionalCondition) const {
0263   if (theMeasurementTracker == nullptr) {
0264     throw cms::Exception("LogicError")
0265         << "Asking to create trajectories to an un-initialized GroupedCkfTrajectoryBuilder.\nYou have to call "
0266            "clone(const MeasurementTrackerEvent *data) and then call trajectories on it instead.\n";
0267   }
0268
0269   statCount.seed();
0270   //
0271   // Build trajectory outwards from seed
0272   //
0273
0274   analyseSeed(seed);
0275
0276   TempTrajectory const& startingTraj = createStartingTrajectory(seed);
0277
0278   work_.clear();
0279   const bool inOut = true;
0280   nCandPerSeed = groupedLimitedCandidates(seed, startingTraj, regionalCondition, forwardPropagator(seed), inOut, work_);
0281   if (work_.empty())
0282     return;
0283
0284   // cleaning now done here...
0285   FastTrajectoryCleaner cleaner(theFoundHitBonus, theLostHitPenalty);
0286   cleaner.clean(work_);
0287
0288   std::shared_ptr<const TrajectorySeed> pseed(new TrajectorySeed(seed));
0289   for (auto const& it : work_)
0290     if (it.isValid()) {
0291       result.push_back(it.toTrajectory());
0292       result.back().setSharedSeed(pseed);
0293     }
0294   work_.clear();
0295   if (work_.capacity() > work_MaxSize_) {
0296     TempTrajectoryContainer().swap(work_);
0297     work_.reserve(work_MaxSize_ / 2);
0298   }
0299
0300   analyseResult(result);
0301
0302   LogDebug("CkfPattern") << "GroupedCkfTrajectoryBuilder: returning result of size " << result.size();
0303   statCount.traj(result.size());
0304
0305 #ifdef VI_DEBUG
0306   int kt = 0;
0307   for (auto const& traj : result) {
0308     int chit[7] = {};
0309     for (auto const& tm : traj.measurements()) {
0310       auto const& hit = tm.recHitR();
0311       if (!hit.isValid())
0312         ++chit[0];
0313       if (hit.det() == nullptr)
0314         ++chit[1];
0315       if (trackerHitRTTI::isUndef(hit))
0316         continue;
0317       if (hit.dimension() != 2) {
0318         ++chit[2];
0319       } else if (trackerHitRTTI::isFromDet(hit)) {
0320         auto const& thit = static_cast<BaseTrackerRecHit const&>(hit);
0321         auto const& clus = thit.firstClusterRef();
0322         if (clus.isPixel())
0323           ++chit[3];
0324         else if (thit.isMatched()) {
0325           ++chit[4];
0326         } else if (thit.isProjected()) {
0327           ++chit[5];
0328         } else {
0329           ++chit[6];
0330         }
0331       }
0332     }
0333
0334     std::cout << "ckf " << kt++ << ": ";
0335     for (auto c : chit)
0336       std::cout << c << '/';
0337     std::cout << std::endl;
0338   }
0339 #endif
0340 }
0341
0342 unsigned int GroupedCkfTrajectoryBuilder::groupedLimitedCandidates(const TrajectorySeed& seed,
0343                                                                    TempTrajectory const& startingTraj,
0344                                                                    const TrajectoryFilter* regionalCondition,
0345                                                                    const Propagator* propagator,
0346                                                                    bool inOut,
0347                                                                    TempTrajectoryContainer& result) const {
0348   unsigned int nIter = 1;
0349   unsigned int nCands = 0;  // ignore startingTraj
0350   unsigned int prevNewCandSize = 0;
0351   TempTrajectoryContainer candidates;
0352   TempTrajectoryContainer newCand;
0353   newCand.reserve(theMaxCand);
0354   candidates.push_back(startingTraj);
0355
0356   while (!candidates.empty()) {
0357     newCand.clear();
0358     for (TempTrajectoryContainer::iterator traj = candidates.begin(); traj != candidates.end(); traj++) {
0359       if (!advanceOneLayer(seed, *traj, regionalCondition, propagator, inOut, newCand, result)) {
0360         LogDebug("CkfPattern") << "GCTB: terminating after advanceOneLayer==false";
0361         continue;
0362       }
0363
0364       LogDebug("CkfPattern") << "newCand(1): after advanced one layer:\n" << PrintoutHelper::dumpCandidates(newCand);
0365       // account only new candidates, i.e.
0366       // - 1 candidate -> 1 candidate, don't increase count
0367       // - 1 candidate -> 2 candidates, increase count by 1
0368       nCands += newCand.size() - prevNewCandSize;
0369       prevNewCandSize = newCand.size();
0370
0371       assert((int)newCand.size() <= theMaxCand);
0372     }
0373
0374     LogDebug("CkfPattern") << "newCand.size() at end = " << newCand.size();
0375
0376     if (theIntermediateCleaning) {
0377 #ifdef STANDARD_INTERMEDIARYCLEAN
0378       IntermediateTrajectoryCleaner::clean(newCand);
0379 #else
0380       groupedIntermediaryClean(newCand);
0381 #endif
0382     }
0383     candidates.swap(newCand);
0384
0385     LogDebug("CkfPattern") << "candidates(3): " << result.size() << " candidates after " << nIter++
0386                            << " groupedCKF iteration: \n"
0387                            << PrintoutHelper::dumpCandidates(result) << "\n " << candidates.size()
0388                            << " running candidates are: \n"
0389                            << PrintoutHelper::dumpCandidates(candidates);
0390   }
0391
0392   return nCands;
0393 }
0394
0395 #ifdef EDM_ML_DEBUG
0396 std::string whatIsTheNextStep(TempTrajectory const& traj,
0397                               std::pair<TrajectoryStateOnSurface, std::vector<const DetLayer*> >& stateAndLayers) {
0398   std::stringstream buffer;
0399   vector<const DetLayer*>& nl = stateAndLayers.second;
0400   // #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0401   // #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0402   //B.M. TkLayerName layerName;
0403   //B.M. buffer << "Started from " << layerName(traj.lastLayer())
0404   const BarrelDetLayer* sbdl = dynamic_cast<const BarrelDetLayer*>(traj.lastLayer());
0405   const ForwardDetLayer* sfdl = dynamic_cast<const ForwardDetLayer*>(traj.lastLayer());
0406   if (sbdl) {
0407     buffer << "Started from " << traj.lastLayer() << " r=" << sbdl->specificSurface().radius()
0408            << " phi=" << sbdl->specificSurface().phi() << endl;
0409   } else if (sfdl) {
0410     buffer << "Started from " << traj.lastLayer() << " z " << sfdl->specificSurface().position().z() << " phi "
0411            << sfdl->specificSurface().phi() << endl;
0412   }
0413   buffer << "Trying to go to";
0414   for (vector<const DetLayer*>::iterator il = nl.begin(); il != nl.end(); il++) {
0415     //B.M. buffer << " " << layerName(*il)  << " " << *il << endl;
0416     const BarrelDetLayer* bdl = dynamic_cast<const BarrelDetLayer*>(*il);
0417     const ForwardDetLayer* fdl = dynamic_cast<const ForwardDetLayer*>(*il);
0418
0419     if (bdl)
0420       buffer << " r " << bdl->specificSurface().radius() << endl;
0421     if (fdl)
0422       buffer << " z " << fdl->specificSurface().position().z() << endl;
0423     //buffer << " " << *il << endl;
0424   }
0425   return buffer.str();
0426 }
0427
0428 std::string whatIsTheStateToUse(TrajectoryStateOnSurface& initial,
0429                                 TrajectoryStateOnSurface& stateToUse,
0430                                 const DetLayer* l) {
0431   std::stringstream buffer;
0432   buffer << "GCTB: starting from "
0433          << " r / phi / z = " << stateToUse.globalPosition().perp() << " / " << stateToUse.globalPosition().phi()
0434          << " / " << stateToUse.globalPosition().z()
0435          << " , pt / phi / pz /charge = " << stateToUse.globalMomentum().perp() << " / "
0436          << stateToUse.globalMomentum().phi() << " / " << stateToUse.globalMomentum().z() << " / "
0437          << stateToUse.charge() << " for layer at " << l << endl;
0438   buffer << "     errors:";
0439   for (int i = 0; i < 5; i++)
0440     buffer << " " << sqrt(stateToUse.curvilinearError().matrix()(i, i));
0441   buffer << endl;
0442
0443   //buffer << "GCTB: starting from r / phi / z = " << initial.globalPosition().perp()
0444   //<< " / " << initial.globalPosition().phi()
0445   //<< " / " << initial.globalPosition().z() << " , pt / pz = "
0446   //<< initial.globalMomentum().perp() << " / "
0447   //<< initial.globalMomentum().z() << " for layer at "
0448   //<< l << endl;
0449   //buffer << "     errors:";
0450   //for ( int i=0; i<5; i++ )  buffer << " " << sqrt(initial.curvilinearError().matrix()(i,i));
0451   //buffer << endl;
0452   return buffer.str();
0453 }
0454 #endif
0455
0456 bool GroupedCkfTrajectoryBuilder::advanceOneLayer(const TrajectorySeed& seed,
0457                                                   TempTrajectory& traj,
0458                                                   const TrajectoryFilter* regionalCondition,
0459                                                   const Propagator* propagator,
0460                                                   bool inOut,
0461                                                   TempTrajectoryContainer& newCand,
0462                                                   TempTrajectoryContainer& result) const {
0463   std::pair<TSOS, std::vector<const DetLayer*> >&& stateAndLayers = findStateAndLayers(seed, traj);
0464
0465   bool full = (int)newCand.size() == theMaxCand;
0466   auto lessTraj = GroupedTrajCandLess(theLostHitPenalty, theFoundHitBonus);
0467
0468   if (maxPt2ForLooperReconstruction > 0) {
0469     if (
0470         //stateAndLayers.second.size()==0 &&
0471         traj.lastLayer()->location() == 0) {
0472       float pt2 = stateAndLayers.first.globalMomentum().perp2();
0473       if (pt2 < maxPt2ForLooperReconstruction && pt2 > (0.3f * 0.3f))
0474         stateAndLayers.second.push_back(traj.lastLayer());
0475     }
0476   }
0477
0478   auto layerBegin = stateAndLayers.second.begin();
0479   auto layerEnd = stateAndLayers.second.end();
0480
0481   //   if (nl.empty()) {
0482   //     addToResult(traj,result,inOut);
0483   //     return false;
0484   //   }
0485
0486 #ifdef EDM_ML_DEBUG
0487   LogDebug("CkfPattern") << whatIsTheNextStep(traj, stateAndLayers);
0488 #endif
0489
0490   bool foundSegments(false);
0491   bool foundNewCandidates(false);
0492   for (auto il = layerBegin; il != layerEnd; il++) {
0493     TSOS stateToUse = stateAndLayers.first;
0494
0495     double dPhiCacheForLoopersReconstruction(0);
0496     if UNLIKELY (!traj.empty() && (*il) == traj.lastLayer()) {
0497       if (maxPt2ForLooperReconstruction > 0) {
0498         // ------ For loopers reconstruction
0499         //cout<<" self propagating in advanceOneLayer (for loopers) \n";
0500         const BarrelDetLayer* sbdl = dynamic_cast<const BarrelDetLayer*>(traj.lastLayer());
0501         if (sbdl) {
0502           HelixBarrelCylinderCrossing cylinderCrossing(stateToUse.globalPosition(),
0503                                                        stateToUse.globalMomentum(),
0504                                                        stateToUse.transverseCurvature(),
0505                                                        propagator->propagationDirection(),
0506                                                        sbdl->specificSurface());
0507           if (!cylinderCrossing.hasSolution())
0508             continue;
0509           GlobalPoint starting = stateToUse.globalPosition();
0510           GlobalPoint target1 = cylinderCrossing.position1();
0511           GlobalPoint target2 = cylinderCrossing.position2();
0512
0513           GlobalPoint farther =
0514               fabs(starting.phi() - target1.phi()) > fabs(starting.phi() - target2.phi()) ? target1 : target2;
0515
0516           const Bounds& bounds(sbdl->specificSurface().bounds());
0517           float length = 0.5f * bounds.length();
0518
0519           /*
0520           cout << "starting: " << starting << endl;
0521           cout << "target1: " << target1 << endl;
0522           cout << "target2: " << target2 << endl;
0523           cout << "dphi: " << (target1.phi()-target2.phi()) << endl;
0524         cout << "length: " << length << endl;
0525         */
0526
0527           /*
0528           float deltaZ = bounds.thickness()/2.f/fabs(tan(stateToUse.globalDirection().theta()) ) ;
0529           if(stateToUse.hasError())
0530           deltaZ += 3*sqrt(stateToUse.cartesianError().position().czz());
0531           if( fabs(farther.z()) > length + deltaZ ) continue;
0532         */
0533           if (fabs(farther.z()) * 0.95f > length)
0534             continue;
0535
0536           Geom::Phi<float> tmpDphi = target1.phi() - target2.phi();
0537           if (std::abs(tmpDphi) > maxDPhiForLooperReconstruction)
0538             continue;
0539           GlobalPoint target(0.5f * (target1.basicVector() + target2.basicVector()));
0540           //cout << "target: " << target << endl;
0541
0542           TransverseImpactPointExtrapolator extrapolator;
0543           stateToUse = extrapolator.extrapolate(stateToUse, target, *propagator);
0544           if (!stateToUse.isValid())
0545             continue;  //SK: consider trying the original? probably not
0546
0547           //dPhiCacheForLoopersReconstruction = fabs(target1.phi()-target2.phi())*2.;
0548           dPhiCacheForLoopersReconstruction = std::abs(tmpDphi);
0549           traj.incrementLoops();
0550         } else {  // not barrel
0551           continue;
0552         }
0553       } else {  // loopers not requested (why else???)
0554                 // ------ For cosmics reconstruction
0555         LogDebug("CkfPattern") << " self propagating in advanceOneLayer.\n from: \n" << stateToUse;
0556         //self navigation case
0557         // go to a middle point first
0558         TransverseImpactPointExtrapolator middle;
0559         GlobalPoint center(0, 0, 0);
0560         stateToUse = middle.extrapolate(stateToUse, center, *(forwardPropagator(seed)));
0561
0562         if (!stateToUse.isValid())
0563           continue;
0564         LogDebug("CkfPattern") << "to: " << stateToUse;
0565       }
0566     }  // last layer...
0567
0568     //unsigned int maxCandidates = theMaxCand > 21 ? theMaxCand*2 : 42; //limit the number of returned segments
0569     LayerMeasurements layerMeasurements(theMeasurementTracker->measurementTracker(), *theMeasurementTracker);
0570     TrajectorySegmentBuilder layerBuilder(
0571         &layerMeasurements, **il, *propagator, *theUpdator, *theEstimator, theLockHits, theBestHitOnly, theMaxCand);
0572
0573 #ifdef EDM_ML_DEBUG
0574     LogDebug("CkfPattern") << whatIsTheStateToUse(stateAndLayers.first, stateToUse, *il);
0575 #endif
0576
0577     auto&& segments = layerBuilder.segments(stateToUse);
0578
0579     LogDebug("CkfPattern") << "GCTB: number of segments = " << segments.size();
0580
0581     if (!segments.empty())
0582       foundSegments = true;
0583
0584     for (auto is = segments.begin(); is != segments.end(); is++) {
0585       //
0586       // assume "invalid hit only" segment is last in list
0587       //
0588       auto const& measurements = is->measurements();
0589       if (!theAlwaysUseInvalid && is != segments.begin() && measurements.size() == 1 &&
0590           (measurements.front().recHit()->getType() == TrackingRecHit::missing))
0591         break;
0592
0593       //----  avoid to add the same hits more than once in the trajectory ----
0594       bool toBeRejected(false);
0595       for (auto revIt = measurements.rbegin(); revIt != measurements.rend(); --revIt) {
0596         // int tmpCounter(0);
0597         for (auto newTrajMeasIt = traj.measurements().rbegin(); newTrajMeasIt != traj.measurements().rend();
0598              --newTrajMeasIt) {
0599           //if(tmpCounter==2) break;
0600           if (revIt->recHitR().geographicalId() == newTrajMeasIt->recHitR().geographicalId() &&
0601               (revIt->recHitR().geographicalId() != DetId(0))) {
0602             toBeRejected = true;
0603             goto rejected;  //break;  // see http://stackoverflow.com/questions/1257744/can-i-use-break-to-exit-multiple-nested-for-loops
0604           }
0605           // tmpCounter++;
0606         }
0607       }
0608
0609     rejected:;  // http://xkcd.com/292/
0610       if (toBeRejected) {
0611 #ifdef VI_DEBUG
0612         cout << "WARNING: neglect candidate because it contains the same hit twice \n";
0613         cout << "-- discarded track's pt,eta,#found/lost: "
0614              << traj.lastMeasurement().updatedState().globalMomentum().perp() << " , "
0615              << traj.lastMeasurement().updatedState().globalMomentum().eta() << " , " << traj.foundHits() << '/'
0616              << traj.lostHits() << "\n";
0617 #endif
0618         traj.setDPhiCacheForLoopersReconstruction(dPhiCacheForLoopersReconstruction);
0619         continue;  //Are we sure about this????
0620       }
0621       // ------------------------
0622
0623       //
0624       // create new candidate
0625       //
0626       TempTrajectory newTraj(traj);
0627       traj.setDPhiCacheForLoopersReconstruction(dPhiCacheForLoopersReconstruction);
0628       newTraj.join(*is);
0629
0630       //std::cout << "DEBUG: newTraj after push found,lost: "
0631       //      << newTraj.foundHits() << " , "
0632       //      << newTraj.lostHits() << " , "
0633       //      << newTraj.measurements().size() << std::endl;
0634
0635       //GIO// for ( vector<TM>::const_iterator im=measurements.begin();
0636       //GIO//        im!=measurements.end(); im++ )  newTraj.push(*im);
0637       //if ( toBeContinued(newTraj,regionalCondition) ) { TOBE FIXED
0638       if (toBeContinued(newTraj, inOut)) {
0639         // Have added one more hit to track candidate
0640
0641         LogDebug("CkfPattern") << "GCTB: adding updated trajectory to candidates: inOut=" << inOut
0642                                << " hits=" << newTraj.foundHits();
0643
0644         newTraj.setStopReason(StopReason::NOT_STOPPED);
0645         if (full) {
0646           bool better = lessTraj(newTraj, newCand.front());
0647           if (better) {
0648             // replace worst
0649             foundNewCandidates = true;
0650             std::pop_heap(newCand.begin(), newCand.end(), lessTraj);
0651             newCand.back().swap(newTraj);
0652             std::push_heap(newCand.begin(), newCand.end(), lessTraj);
0653           }  // else? no need to add it just to remove it later!
0654         } else {
0655           newCand.push_back(std::move(newTraj));
0656           foundNewCandidates = true;
0657           full = (int)newCand.size() == theMaxCand;
0658           if (full)
0659             std::make_heap(newCand.begin(), newCand.end(), lessTraj);
0660         }
0661       } else {
0662         // Have finished building this track. Check if it passes cuts.
0663
0664         LogDebug("CkfPattern") << "GCTB: adding completed trajectory to results if passes cuts: inOut=" << inOut
0665                                << " hits=" << newTraj.foundHits();
0666         moveToResult(std::move(newTraj), result, inOut);
0667       }
0668     }  // loop over segs
0669   }    // loop over layers
0670
0671   if (!foundSegments) {
0672     LogDebug("CkfPattern") << "GCTB: adding input trajectory to result";
0673     if (!stateAndLayers.second.empty())
0674       traj.setStopReason(StopReason::NO_SEGMENTS_FOR_VALID_LAYERS);
0675     addToResult(traj, result, inOut);
0676   }
0677   return foundNewCandidates;
0678 }
0679
0680 namespace {
0681   /// fills in a list of layers from a container of TrajectoryMeasurements
0682   ///
0683   struct LayersInTraj {
0684     static constexpr int N = 3;
0685     TempTrajectory* traj;
0686     std::array<DetLayer const*, N> layers;
0687     int tot;
0688     void fill(TempTrajectory& t) {
0689       traj = &t;
0690       tot = 0;
0691       const TempTrajectory::DataContainer& measurements = traj->measurements();
0692
0693       auto currl = layers[tot] = measurements.back().layer();
0694       TempTrajectory::DataContainer::const_iterator ifirst = measurements.rbegin();
0695       --ifirst;
0696       for (TempTrajectory::DataContainer::const_iterator im = ifirst; im != measurements.rend(); --im) {
0697         if (im->layer() != currl) {
0698           ++tot;
0699           currl = im->layer();
0700           if (tot < N)
0701             layers[tot] = currl;
0702         }
0703       }
0704       ++tot;
0705     }
0706
0707     //void verify() {
0708     //  for (vector<const DetLayer*>::const_iterator iter = result.begin(); iter != result.end(); iter++)
0709     //  if (!*iter) edm::LogWarning("CkfPattern")<< "Warning: null det layer!! ";
0710     // }
0711   };
0712 }  // namespace
0713
0714 //TempTrajectoryContainer
0715 void GroupedCkfTrajectoryBuilder::groupedIntermediaryClean(TempTrajectoryContainer& theTrajectories) const {
0716   //if (theTrajectories.empty()) return TrajectoryContainer();
0717   //TrajectoryContainer result;
0718   if (theTrajectories.empty())
0719     return;
0720   //RecHitEqualByChannels recHitEqualByChannels(false, false);
0721   LayersInTraj layers[theTrajectories.size()];
0722   int ntraj = 0;
0723   for (auto& t : theTrajectories) {
0724     if (t.isValid() && t.lastMeasurement().recHitR().isValid())
0725       layers[ntraj++].fill(t);
0726   }
0727
0728   if (ntraj < 2)
0729     return;
0730
0731   for (int ifirst = 0; ifirst != ntraj - 1; ++ifirst) {
0732     auto firstTraj = layers[ifirst].traj;
0733     if (!firstTraj->isValid())
0734       continue;
0735     const TempTrajectory::DataContainer& firstMeasurements = firstTraj->measurements();
0736
0737     int firstLayerSize = layers[ifirst].tot;
0738     if (firstLayerSize < 4)
0739       continue;
0740     auto const& firstLayers = layers[ifirst].layers;
0741
0742     for (int isecond = ifirst + 1; isecond != ntraj; ++isecond) {
0743       auto secondTraj = layers[isecond].traj;
0744       if (!secondTraj->isValid())
0745         continue;
0746
0747       const TempTrajectory::DataContainer& secondMeasurements = secondTraj->measurements();
0748
0749       int secondLayerSize = layers[isecond].tot;
0750       //
0751       // only candidates using the same last 3 layers are compared
0752       //
0753       if (firstLayerSize != secondLayerSize)
0754         continue;  // V.I.  why equal???
0755       auto const& secondLayers = layers[isecond].layers;
0756       if (firstLayers[0] != secondLayers[0] || firstLayers[1] != secondLayers[1] || firstLayers[2] != secondLayers[2])
0757         continue;
0758
0759       TempTrajectory::DataContainer::const_iterator im1 = firstMeasurements.rbegin();
0760       TempTrajectory::DataContainer::const_iterator im2 = secondMeasurements.rbegin();
0761       //
0762       // check for identical hits in the last layer
0763       //
0764       bool unequal(false);
0765       const DetLayer* layerPtr = firstLayers[0];
0766       while (im1 != firstMeasurements.rend() && im2 != secondMeasurements.rend()) {
0767         if (im1->layer() != layerPtr || im2->layer() != layerPtr)
0768           break;
0769         if (!(im1->recHit()->isValid()) || !(im2->recHit()->isValid()) ||
0770             !im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::some)) {
0771           //!im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::all) ) {
0772           unequal = true;
0773           break;
0774         }
0775         --im1;
0776         --im2;
0777       }
0778       if (im1 == firstMeasurements.rend() || im2 == secondMeasurements.rend() || im1->layer() == layerPtr ||
0779           im2->layer() == layerPtr || unequal)
0780         continue;
0781       //
0782       // check for invalid hits in the layer -2
0783       // compare only candidates with invalid / valid combination
0784       //
0785       layerPtr = firstLayers[1];
0786       bool firstValid(true);
0787       while (im1 != firstMeasurements.rend() && im1->layer() == layerPtr) {
0788         if (!im1->recHit()->isValid())
0789           firstValid = false;
0790         --im1;
0791       }
0792       bool secondValid(true);
0793       while (im2 != secondMeasurements.rend() && im2->layer() == layerPtr) {
0794         if (!im2->recHit()->isValid())
0795           secondValid = false;
0796         --im2;
0797       }
0798       if (!tkxor(firstValid, secondValid))
0799         continue;
0800       //
0801       // ask for identical hits in layer -3
0802       //
0803       unequal = false;
0804       layerPtr = firstLayers[2];
0805       while (im1 != firstMeasurements.rend() && im2 != secondMeasurements.rend()) {
0806         if (im1->layer() != layerPtr || im2->layer() != layerPtr)
0807           break;
0808         if (!(im1->recHit()->isValid()) || !(im2->recHit()->isValid()) ||
0809             !im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::some)) {
0810           //!im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::all) ) {
0811           unequal = true;
0812           break;
0813         }
0814         --im1;
0815         --im2;
0816       }
0817       if (im1 == firstMeasurements.rend() || im2 == secondMeasurements.rend() || im1->layer() == layerPtr ||
0818           im2->layer() == layerPtr || unequal)
0819         continue;
0820
0821       if (!firstValid) {
0822         firstTraj->invalidate();
0823         break;
0824       } else {
0825         secondTraj->invalidate();  // V.I. why break?
0826         break;
0827       }
0828     }  // second
0829   }    // first
0830        /*
0831   for (TempTrajectoryContainer::const_iterator it = theTrajectories.begin();
0832        it != theTrajectories.end(); it++) {
0833     if(it->isValid()) result.push_back( *it);
0834   }
0835
0836   return result;
0837 */
0838   theTrajectories.erase(
0839       std::remove_if(theTrajectories.begin(), theTrajectories.end(), std::not_fn(&TempTrajectory::isValid)),
0840       theTrajectories.end());
0841 }
0842
0843 void GroupedCkfTrajectoryBuilder::rebuildSeedingRegion(const TrajectorySeed& seed,
0844                                                        TempTrajectoryContainer& result) const {
0845   //
0846   // Rebuilding of trajectories. Candidates are taken from result,
0847   // which will be replaced with the solutions after rebuild
0848   // (assume vector::swap is more efficient than building new container)
0849   //
0850   LogDebug("CkfPattern") << "Starting to rebuild " << result.size() << " tracks";
0851   //
0852   // Fitter (need to create it here since the propagation direction
0853   // might change between different starting trajectories)
0854   //
0855   auto hitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(hitBuilder())->cloner();
0856   KFTrajectoryFitter fitter(backwardPropagator(seed), &updator(), &estimator(), 3, nullptr, &hitCloner);
0857   //
0858   std::vector<const TrackingRecHit*> seedHits;
0859
0860   unsigned int nSeed = seed.nHits();
0861   //seedHits.reserve(nSeed);
0862   TempTrajectoryContainer rebuiltTrajectories;
0863
0864   for (TempTrajectoryContainer::iterator it = result.begin(); it != result.end(); it++) {
0865     // Refit - keep existing trajectory in case fit is not possible
0866     // or fails
0867     //
0868     assert(it->isValid());
0869     auto&& reFitted = backwardFit(*it, nSeed, fitter, seedHits);
0870     if UNLIKELY (!reFitted.isValid()) {
0871       rebuiltTrajectories.push_back(std::move(*it));
0872       LogDebug("CkfPattern") << "RebuildSeedingRegion skipped as backward fit failed";
0873       //                << "after reFitted.size() " << reFitted.size();
0874       continue;
0875     }
0876     //LogDebug("CkfPattern")<<"after reFitted.size() " << reFitted.size();
0877     //
0878     // Rebuild seeding part. In case it fails: keep initial trajectory
0879     // (better to drop it??)
0880     //
0881     int nRebuilt = rebuildSeedingRegion(seed, seedHits, reFitted, rebuiltTrajectories);
0882     // Loop over the last nRebuilt trajectories and propagate back the
0883     // real cause that stopped the original in-out trajectory, since
0884     // that's the one we want to monitor
0885     for (size_t i = rebuiltTrajectories.size() - 1; i < rebuiltTrajectories.size() - nRebuilt - 1; --i) {
0886       rebuiltTrajectories[i].setStopReason(it->stopReason());
0887     }
0888
0889     if (nRebuilt == 0 && !theKeepOriginalIfRebuildFails)
0890       it->invalidate();  // won't use original in-out track
0891
0892     if (nRebuilt < 0)
0893       rebuiltTrajectories.push_back(std::move(*it));
0894   }
0895   //
0896   // Replace input trajectories with new ones
0897   //
0898   result.swap(rebuiltTrajectories);
0899   result.erase(std::remove_if(result.begin(), result.end(), std::not_fn(&TempTrajectory::isValid)), result.end());
0900 }
0901
0902 int GroupedCkfTrajectoryBuilder::rebuildSeedingRegion(const TrajectorySeed& seed,
0903                                                       const std::vector<const TrackingRecHit*>& seedHits,
0904                                                       TempTrajectory& candidate,
0905                                                       TempTrajectoryContainer& result) const {
0906   //
0907   // Starting from track found by in-out tracking phase, extrapolate it inwards through
0908   // the seeding region if possible in towards smaller Tracker radii, searching for additional
0909   // hits.
0910   // The resulting trajectories are returned in result,
0911   // the count is the return value.
0912   //
0913   TempTrajectoryContainer rebuiltTrajectories;
0914 #ifdef DBG2_GCTB
0915   /*  const LayerFinderByDet layerFinder;
0916   if ( !seedHits.empty() && seedHits.front().isValid() ) {
0917     DetLayer* seedLayer = layerFinder(seedHits.front().det());
0918     cout << "Seed hit at " << seedHits.front().globalPosition()
0919      << " " << seedLayer << endl;
0920     cout << "Started from "
0921      << candidate.lastMeasurement().updatedState().globalPosition().perp() << " "
0922      << candidate.lastMeasurement().updatedState().globalPosition().z() << endl;
0923     pair<bool,TrajectoryStateOnSurface> layerComp(false,TrajectoryStateOnSurface());
0924     if ( seedLayer ) layerComp =
0925       seedLayer->compatible(candidate.lastMeasurement().updatedState(),
0926                   propagator(),estimator());
0927     pair<bool,TrajectoryStateOnSurface> detComp =
0928       seedHits.front().det().compatible(candidate.lastMeasurement().updatedState(),
0929                     propagator(),estimator());
0930     cout << "  layer compatibility = " << layerComp.first;
0931     cout << "  det compatibility = " << detComp.first;
0932     if ( detComp.first ) {
0933       cout << "  estimate = "
0934        << estimator().estimate(detComp.second,seedHits.front()).second ;
0935     }
0936     cout << endl;
0937   }*/
0938   cout << "Before backward building: #measurements = " << candidate.measurements().size();  //<< endl;;
0939 #endif
0940   //
0941   // Use standard building with standard cuts. Maybe better to use different
0942   // cuts from "forward" building (e.g. no check on nr. of invalid hits)?
0943   //
0944   const bool inOut = false;
0945   groupedLimitedCandidates(seed, candidate, nullptr, backwardPropagator(seed), inOut, rebuiltTrajectories);
0946
0947   LogDebug("CkfPattern") << " After backward building: " << PrintoutHelper::dumpCandidates(rebuiltTrajectories);
0948
0949   //
0950   // Check & count resulting candidates
0951   //
0952   int nrOfTrajectories(0);
0953   bool orig_ok = false;
0954   //const RecHitEqualByChannels recHitEqual(false,false);
0955   //vector<TM> oldMeasurements(candidate.measurements());
0956   for (TempTrajectoryContainer::iterator it = rebuiltTrajectories.begin(); it != rebuiltTrajectories.end(); it++) {
0957     TempTrajectory::DataContainer newMeasurements(it->measurements());
0958     //
0959     // Verify presence of seeding hits?
0960     //
0961     if (theRequireSeedHitsInRebuild) {
0962       orig_ok = true;
0963       // no hits found (and possibly some invalid hits discarded): drop track
0964       if (newMeasurements.size() <= candidate.measurements().size()) {
0965         LogDebug("CkfPattern") << "newMeasurements.size()<=candidate.measurements().size()";
0966         continue;
0967       }
0968       // verify presence of hits
0969       //GIO//if ( !verifyHits(newMeasurements.begin()+candidate.measurements().size(),
0970       //GIO//              newMeasurements.end(),seedHits) ){
0971       if (!verifyHits(newMeasurements.rbegin(), newMeasurements.size() - candidate.measurements().size(), seedHits)) {
0972         LogDebug("CkfPattern") << "seed hits not found in rebuild";
0973         continue;
0974       }
0975     }
0976     //
0977     // construct final trajectory in the right order
0978     //
0979     // save & count result
0980     nrOfTrajectories++;
0981     result.emplace_back(seed.direction(), seed.nHits());
0982     TempTrajectory& reversedTrajectory = result.back();
0983     reversedTrajectory.setNLoops(it->nLoops());
0984     for (TempTrajectory::DataContainer::const_iterator im = newMeasurements.rbegin(), ed = newMeasurements.rend();
0985          im != ed;
0986          --im) {
0987       reversedTrajectory.push(*im);
0988     }
0989
0990     LogDebug("CkgPattern") << "New traj direction = " << reversedTrajectory.direction() << "\n"
0991                            << PrintoutHelper::dumpMeasurements(reversedTrajectory.measurements());
0992   }  // rebuiltTrajectories
0993
0994   // If nrOfTrajectories = 0 and orig_ok = true, this means that a track was actually found on the
0995   // out-in step (meeting those requirements) but did not have the seed hits in it.
0996   // In this case when we return we will go ahead and use the original in-out track.
0997
0998   // If nrOfTrajectories = 0 and orig_ok = false, this means that the out-in step failed to
0999   // find any track.  Two cases are a technical failure in fitting the original seed hits or
1000   // because the track did not meet the out-in criteria (which may be stronger than the out-in
1001   // criteria).  In this case we will NOT allow the original in-out track to be used.
1002
1003   if ((nrOfTrajectories == 0) && orig_ok) {
1004     nrOfTrajectories = -1;
1005   }
1006   return nrOfTrajectories;
1007 }
1008
1009 TempTrajectory GroupedCkfTrajectoryBuilder::backwardFit(TempTrajectory& candidate,
1010                                                         unsigned int nSeed,
1011                                                         const TrajectoryFitter& fitter,
1012                                                         std::vector<const TrackingRecHit*>& remainingHits) const {
1013   //
1014   // clear array of non-fitted hits
1015   //
1016   remainingHits.clear();
1017
1018   LogDebug("CkfPattern") << "nSeed " << nSeed << endl
1019                          << "Old traj direction = " << candidate.direction() << endl
1020                          << PrintoutHelper::dumpMeasurements(candidate.measurements());
1021
1022   //
1023   // backward fit trajectory.
1024   // (Will try to fit only hits outside the seeding region. However,
1025   // if there are not enough of these, it will also use the seeding hits).
1026   //
1027   //   const unsigned int nHitAllMin(5);
1028   //   const unsigned int nHit2dMin(2);
1029   unsigned int nHit(0);  // number of valid hits after seeding region
1030   //unsigned int nHit2d(0);  // number of valid hits after seeding region with 2D info
1031   // use all hits except the first n (from seed), but require minimum
1032   // specified in configuration.
1033   //  Swapped over next two lines.
1034   unsigned int nHitMin = std::max(candidate.foundHits() - nSeed, theMinNrOfHitsForRebuild);
1035   //  unsigned int nHitMin = oldMeasurements.size()-nSeed;
1036   // we want to rebuild only if the number of VALID measurements excluding the seed measurements is higher than the cut
1037   if UNLIKELY (nHitMin < theMinNrOfHitsForRebuild)
1038     return TempTrajectory();
1039
1040   LogDebug("CkfPattern") /* << "nHitMin " << nHitMin*/ << "Sizes: " << candidate.measurements().size() << " / ";
1041   //
1042   // create input trajectory for backward fit
1043   //
1044   Trajectory fwdTraj(oppositeDirection(candidate.direction()));
1045   fwdTraj.setNLoops(candidate.nLoops());
1046   //const TrajectorySeed seed = TrajectorySeed(PTrajectoryStateOnDet(), TrajectorySeed::recHitContainer(), oppositeDirection(candidate.direction()));
1047   //Trajectory fwdTraj(seed, oppositeDirection(candidate.direction()));
1048
1049   const DetLayer* bwdDetLayer[candidate.measurements().size()];
1050   int nl = 0;
1051   for (auto const& tm : candidate.measurements()) {
1052     const TrackingRecHit* hit = tm.recHitR().hit();
1053     //
1054     // add hits until required number is reached
1055     //
1056     if (nHit < nHitMin) {  //|| nHit2d<theMinNrOf2dHitsForRebuild ) {
1057       fwdTraj.push(tm);
1058       bwdDetLayer[nl++] = tm.layer();
1059       //
1060       // count valid / 2D hits
1061       //
1062       if LIKELY (hit->isValid()) {
1063         nHit++;
1064         //if ( hit.isMatched() ||
1065         //     hit.det().detUnits().front()->type().module()==pixel )
1066         //nHit2d++;
1067       }
1068     }
1069     //if (nHit==nHitMin) lastBwdDetLayer=im->layer();
1070     //
1071     // keep remaining (valid) hits for verification
1072     //
1073     else if (hit->isValid()) {
1074       //std::cout << "Adding a remaining hit" << std::endl;
1075       remainingHits.push_back(hit);
1076     }
1077   }
1078   //
1079   // Fit only if required number of valid hits can be used
1080   //
1081   if UNLIKELY (nHit < nHitMin)
1082     return TempTrajectory();
1083
1084   //
1085   // Do the backward fit (important: start from scaled, not random cov. matrix!)
1086   //
1087   TrajectoryStateOnSurface firstTsos(fwdTraj.firstMeasurement().updatedState());
1088   //cout << "firstTsos "<< firstTsos << endl;
1089   firstTsos.rescaleError(10.);
1090   //TrajectoryContainer bwdFitted(fitter.fit(fwdTraj.seed(),fwdTraj.recHits(),firstTsos));
1091   Trajectory&& bwdFitted =
1092       fitter.fitOne(TrajectorySeed({}, {}, oppositeDirection(candidate.direction())), fwdTraj.recHits(), firstTsos);
1093   if UNLIKELY (!bwdFitted.isValid())
1094     return TempTrajectory();
1095
1096   LogDebug("CkfPattern") << "Obtained bwdFitted trajectory with measurement size " << bwdFitted.measurements().size();
1097   TempTrajectory fitted(fwdTraj.direction(), nSeed);
1098   fitted.setNLoops(fwdTraj.nLoops());
1099   vector<TM> const& tmsbf = bwdFitted.measurements();
1100   int iDetLayer = 0;
1101   //this is ugly but the TM in the fitted track do not contain the DetLayer.
1102   //So we have to cache the detLayer pointers and replug them in.
1103   //For the backward building it would be enaugh to cache the last DetLayer,
1104   //but for the intermediary cleaning we need all
1105   for (vector<TM>::const_iterator im = tmsbf.begin(); im != tmsbf.end(); im++) {
1106     fitted.emplace((*im).forwardPredictedState(),
1107                    (*im).backwardPredictedState(),
1108                    (*im).updatedState(),
1109                    (*im).recHit(),
1110                    (*im).estimate(),
1111                    bwdDetLayer[iDetLayer]);
1112
1113     LogDebug("CkfPattern") << PrintoutHelper::dumpMeasurement(*im);
1114     iDetLayer++;
1115   }
1116   /*
1117     TM lastMeas = bwdFitted.front().lastMeasurement();
1118     fitted.pop();
1119     fitted.push(TM(lastMeas.forwardPredictedState(),
1120     lastMeas.backwardPredictedState(),
1121     lastMeas.updatedState(),
1122     lastMeas.recHit(),
1123     lastMeas.estimate(),
1124     lastBwdDetLayer));
1125   */
1126
1127   return fitted;
1128 }
1129
1130 bool GroupedCkfTrajectoryBuilder::verifyHits(TempTrajectory::DataContainer::const_iterator rbegin,
1131                                              size_t maxDepth,
1132                                              const std::vector<const TrackingRecHit*>& hits) const {
1133   //
1134   // verify presence of the seeding hits
1135   //
1136   LogDebug("CkfPattern") << "Checking for " << hits.size() << " hits in " << maxDepth << " measurements" << endl;
1137
1138   auto rend = rbegin;
1139   while (maxDepth > 0) {
1140     --maxDepth;
1141     --rend;
1142   }
1143   for (auto ir = hits.begin(); ir != hits.end(); ir++) {
1144     // assume that all seeding hits are valid!
1145     bool foundHit(false);
1146     for (auto im = rbegin; im != rend; --im) {
1147       if (im->recHit()->isValid() && (*ir)->sharesInput(im->recHit()->hit(), TrackingRecHit::some)) {
1148         foundHit = true;
1149         break;
1150       }
1151     }
1152     if (!foundHit)
1153       return false;
1154   }
1155   return true;
1156 }