CkfPattern/plugins/GroupedCkfTrajectoryBuilder.cc

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