File indexing completed on 2024-12-20 03:14:09
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
0041 #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h"
0042
0043 #include "TrackingTools/TrajectoryCleaning/interface/FastTrajectoryCleaner.h"
0044
0045
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
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 }
0086
0087 using namespace std;
0088
0089
0090
0091
0092
0093 #ifdef STANDARD_INTERMEDIARYCLEAN
0094 #include "RecoTracker/CkfPattern/interface/IntermediateTrajectoryCleaner.h"
0095 #endif
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
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
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
0143
0144
0145
0146
0147
0148
0149
0150
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
0182
0183
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, ®ionalCondition);
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, ®ionalCondition);
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
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
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
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
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;
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
0367
0368
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
0402
0403
0404
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
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
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
0445
0446
0447
0448
0449
0450
0451
0452
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
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
0483
0484
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
0500
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
0522
0523
0524
0525
0526
0527
0528
0529
0530
0531
0532
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
0542
0543 TransverseImpactPointExtrapolator extrapolator;
0544 stateToUse = extrapolator.extrapolate(stateToUse, target, *propagator);
0545 if (!stateToUse.isValid())
0546 continue;
0547
0548
0549 dPhiCacheForLoopersReconstruction = std::abs(tmpDphi);
0550 traj.incrementLoops();
0551 } else {
0552 continue;
0553 }
0554 } else {
0555
0556 LogDebug("CkfPattern") << " self propagating in advanceOneLayer.\n from: \n" << stateToUse;
0557
0558
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 }
0568
0569
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
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
0595 bool toBeRejected(false);
0596 for (auto revIt = measurements.rbegin(); revIt != measurements.rend(); --revIt) {
0597
0598 for (auto newTrajMeasIt = traj.measurements().rbegin(); newTrajMeasIt != traj.measurements().rend();
0599 --newTrajMeasIt) {
0600
0601 if (revIt->recHitR().geographicalId() == newTrajMeasIt->recHitR().geographicalId() &&
0602 (revIt->recHitR().geographicalId() != DetId(0))) {
0603 toBeRejected = true;
0604 goto rejected;
0605 }
0606
0607 }
0608 }
0609
0610 rejected:;
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;
0621 }
0622
0623
0624
0625
0626
0627 TempTrajectory newTraj(traj);
0628 traj.setDPhiCacheForLoopersReconstruction(dPhiCacheForLoopersReconstruction);
0629 newTraj.join(*is);
0630
0631
0632
0633
0634
0635
0636
0637
0638
0639 if (toBeContinued(newTraj, inOut)) {
0640
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
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 }
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
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 }
0670 }
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
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
0709
0710
0711
0712 };
0713 }
0714
0715
0716 void GroupedCkfTrajectoryBuilder::groupedIntermediaryClean(TempTrajectoryContainer& theTrajectories) const {
0717
0718
0719 if (theTrajectories.empty())
0720 return;
0721
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
0753
0754 if (firstLayerSize != secondLayerSize)
0755 continue;
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
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
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
0784
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
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
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();
0827 break;
0828 }
0829 }
0830 }
0831
0832
0833
0834
0835
0836
0837
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
0848
0849
0850
0851 LogDebug("CkfPattern") << "Starting to rebuild " << result.size() << " tracks";
0852
0853
0854
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
0863 TempTrajectoryContainer rebuiltTrajectories;
0864
0865 for (TempTrajectoryContainer::iterator it = result.begin(); it != result.end(); it++) {
0866
0867
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
0875 continue;
0876 }
0877
0878
0879
0880
0881
0882 int nRebuilt = rebuildSeedingRegion(seed, seedHits, reFitted, rebuiltTrajectories);
0883
0884
0885
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();
0892
0893 if (nRebuilt < 0)
0894 rebuiltTrajectories.push_back(std::move(*it));
0895 }
0896
0897
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
0909
0910
0911
0912
0913
0914 TempTrajectoryContainer rebuiltTrajectories;
0915 #ifdef DBG2_GCTB
0916
0917
0918
0919
0920
0921
0922
0923
0924
0925
0926
0927
0928
0929
0930
0931
0932
0933
0934
0935
0936
0937
0938
0939 cout << "Before backward building: #measurements = " << candidate.measurements().size();
0940 #endif
0941
0942
0943
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
0952
0953 int nrOfTrajectories(0);
0954 bool orig_ok = false;
0955
0956
0957 for (TempTrajectoryContainer::iterator it = rebuiltTrajectories.begin(); it != rebuiltTrajectories.end(); it++) {
0958 TempTrajectory::DataContainer newMeasurements(it->measurements());
0959
0960
0961
0962 if (theRequireSeedHitsInRebuild) {
0963 orig_ok = true;
0964
0965 if (newMeasurements.size() <= candidate.measurements().size()) {
0966 LogDebug("CkfPattern") << "newMeasurements.size()<=candidate.measurements().size()";
0967 continue;
0968 }
0969
0970
0971
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
0979
0980
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 }
0994
0995
0996
0997
0998
0999
1000
1001
1002
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
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
1025
1026
1027
1028
1029
1030 unsigned int nHit(0);
1031
1032
1033
1034
1035 unsigned int nHitMin = std::max(candidate.foundHits() - nSeed, theMinNrOfHitsForRebuild);
1036
1037
1038 if UNLIKELY (nHitMin < theMinNrOfHitsForRebuild)
1039 return TempTrajectory();
1040
1041 LogDebug("CkfPattern") << "Sizes: " << candidate.measurements().size() << " / ";
1042
1043
1044
1045 Trajectory fwdTraj(oppositeDirection(candidate.direction()));
1046 fwdTraj.setNLoops(candidate.nLoops());
1047
1048
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
1056
1057 if (nHit < nHitMin) {
1058 fwdTraj.push(tm);
1059 bwdDetLayer[nl++] = tm.layer();
1060
1061
1062
1063 if LIKELY (hit->isValid()) {
1064 nHit++;
1065
1066
1067
1068 }
1069 }
1070
1071
1072
1073
1074 else if (hit->isValid()) {
1075
1076 remainingHits.push_back(hit);
1077 }
1078 }
1079
1080
1081
1082 if UNLIKELY (nHit < nHitMin)
1083 return TempTrajectory();
1084
1085
1086
1087
1088 TrajectoryStateOnSurface firstTsos(fwdTraj.firstMeasurement().updatedState());
1089
1090 firstTsos.rescaleError(10.);
1091
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
1103
1104
1105
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
1119
1120
1121
1122
1123
1124
1125
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
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
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 }