Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:22:05 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!