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
0040 #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h"
0041
0042 #include "TrackingTools/TrajectoryCleaning/interface/FastTrajectoryCleaner.h"
0043
0044
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
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 }
0085
0086 using namespace std;
0087
0088
0089
0090
0091
0092 #ifdef STANDARD_INTERMEDIARYCLEAN
0093 #include "RecoTracker/CkfPattern/interface/IntermediateTrajectoryCleaner.h"
0094 #endif
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
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
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
0142
0143
0144
0145
0146
0147
0148
0149
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
0181
0182
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, ®ionalCondition);
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, ®ionalCondition);
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
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
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
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
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;
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
0370
0371
0372 nCands += newCand.size() - prevNewCandSize;
0373 prevNewCandSize = newCand.size();
0374
0375 if ((int)newCand.size() > theMaxCand) {
0376
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
0391
0392
0393
0394
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
0422
0423
0424
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
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
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
0465
0466
0467
0468
0469
0470
0471
0472
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
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
0500
0501
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
0517
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
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
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
0559
0560 TransverseImpactPointExtrapolator extrapolator;
0561 stateToUse = extrapolator.extrapolate(stateToUse, target, *propagator);
0562 if (!stateToUse.isValid())
0563 continue;
0564
0565
0566 dPhiCacheForLoopersReconstruction = std::abs(tmpDphi);
0567 traj.incrementLoops();
0568 } else {
0569 continue;
0570 }
0571 } else {
0572
0573 LogDebug("CkfPattern") << " self propagating in advanceOneLayer.\n from: \n" << stateToUse;
0574
0575
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 }
0585
0586
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
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
0612 bool toBeRejected(false);
0613 for (auto revIt = measurements.rbegin(); revIt != measurements.rend(); --revIt) {
0614
0615 for (auto newTrajMeasIt = traj.measurements().rbegin(); newTrajMeasIt != traj.measurements().rend();
0616 --newTrajMeasIt) {
0617
0618 if (revIt->recHitR().geographicalId() == newTrajMeasIt->recHitR().geographicalId() &&
0619 (revIt->recHitR().geographicalId() != DetId(0))) {
0620 toBeRejected = true;
0621 goto rejected;
0622 }
0623
0624 }
0625 }
0626
0627 rejected:;
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;
0638 }
0639
0640
0641
0642
0643
0644 TempTrajectory newTraj(traj);
0645 traj.setDPhiCacheForLoopersReconstruction(dPhiCacheForLoopersReconstruction);
0646 newTraj.join(*is);
0647
0648
0649
0650
0651
0652
0653
0654
0655
0656 if (toBeContinued(newTraj, inOut)) {
0657
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
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 }
0673 }
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
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
0712
0713
0714
0715 };
0716 }
0717
0718
0719 void GroupedCkfTrajectoryBuilder::groupedIntermediaryClean(TempTrajectoryContainer& theTrajectories) const {
0720
0721
0722 if (theTrajectories.empty())
0723 return;
0724
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
0756
0757 if (firstLayerSize != secondLayerSize)
0758 continue;
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
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
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
0787
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
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
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();
0830 break;
0831 }
0832 }
0833 }
0834
0835
0836
0837
0838
0839
0840
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
0852
0853
0854
0855 LogDebug("CkfPattern") << "Starting to rebuild " << result.size() << " tracks";
0856
0857
0858
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
0867 TempTrajectoryContainer rebuiltTrajectories;
0868
0869 for (TempTrajectoryContainer::iterator it = result.begin(); it != result.end(); it++) {
0870
0871
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
0879 continue;
0880 }
0881
0882
0883
0884
0885
0886 int nRebuilt = rebuildSeedingRegion(seed, seedHits, reFitted, rebuiltTrajectories);
0887
0888
0889
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();
0896
0897 if (nRebuilt < 0)
0898 rebuiltTrajectories.push_back(std::move(*it));
0899 }
0900
0901
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
0913
0914
0915
0916
0917
0918 TempTrajectoryContainer rebuiltTrajectories;
0919 #ifdef DBG2_GCTB
0920
0921
0922
0923
0924
0925
0926
0927
0928
0929
0930
0931
0932
0933
0934
0935
0936
0937
0938
0939
0940
0941
0942
0943 cout << "Before backward building: #measurements = " << candidate.measurements().size();
0944 #endif
0945
0946
0947
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
0956
0957 int nrOfTrajectories(0);
0958 bool orig_ok = false;
0959
0960
0961 for (TempTrajectoryContainer::iterator it = rebuiltTrajectories.begin(); it != rebuiltTrajectories.end(); it++) {
0962 TempTrajectory::DataContainer newMeasurements(it->measurements());
0963
0964
0965
0966 if (theRequireSeedHitsInRebuild) {
0967 orig_ok = true;
0968
0969 if (newMeasurements.size() <= candidate.measurements().size()) {
0970 LogDebug("CkfPattern") << "newMeasurements.size()<=candidate.measurements().size()";
0971 continue;
0972 }
0973
0974
0975
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
0983
0984
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 }
0998
0999
1000
1001
1002
1003
1004
1005
1006
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
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
1029
1030
1031
1032
1033
1034 unsigned int nHit(0);
1035
1036
1037
1038
1039 unsigned int nHitMin = std::max(candidate.foundHits() - nSeed, theMinNrOfHitsForRebuild);
1040
1041
1042 if UNLIKELY (nHitMin < theMinNrOfHitsForRebuild)
1043 return TempTrajectory();
1044
1045 LogDebug("CkfPattern") << "Sizes: " << candidate.measurements().size() << " / ";
1046
1047
1048
1049 Trajectory fwdTraj(oppositeDirection(candidate.direction()));
1050 fwdTraj.setNLoops(candidate.nLoops());
1051
1052
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
1060
1061 if (nHit < nHitMin) {
1062 fwdTraj.push(tm);
1063 bwdDetLayer[nl++] = tm.layer();
1064
1065
1066
1067 if LIKELY (hit->isValid()) {
1068 nHit++;
1069
1070
1071
1072 }
1073 }
1074
1075
1076
1077
1078 else if (hit->isValid()) {
1079
1080 remainingHits.push_back(hit);
1081 }
1082 }
1083
1084
1085
1086 if UNLIKELY (nHit < nHitMin)
1087 return TempTrajectory();
1088
1089
1090
1091
1092 TrajectoryStateOnSurface firstTsos(fwdTraj.firstMeasurement().updatedState());
1093
1094 firstTsos.rescaleError(10.);
1095
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
1107
1108
1109
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
1123
1124
1125
1126
1127
1128
1129
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
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
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 }