Update mkFit to support additional tracking iterations and other technical developments

Configuration/ProcessModifiers/python/trackingMkFitPixelPairStep_cff.py

@@ -0,0 +1,4 @@
+import FWCore.ParameterSet.Config as cms
+
+# This modifier sets replaces the default pattern recognition with mkFit for pixelPairStep
+trackingMkFitPixelPairStep = cms.Modifier()

Configuration/ProcessModifiers/python/trackingMkFit_cff.py

@@ -8,6 +8,7 @@
 from Configuration.ProcessModifiers.trackingMkFitLowPtTripletStep_cff import *
 from Configuration.ProcessModifiers.trackingMkFitDetachedQuadStep_cff import *
 from Configuration.ProcessModifiers.trackingMkFitDetachedTripletStep_cff import *
+from Configuration.ProcessModifiers.trackingMkFitPixelPairStep_cff import *
 from Configuration.ProcessModifiers.trackingMkFitMixedTripletStep_cff import *
 from Configuration.ProcessModifiers.trackingMkFitPixelLessStep_cff import *
 from Configuration.ProcessModifiers.trackingMkFitTobTecStep_cff import *
@@ -22,6 +23,7 @@
     trackingMkFitLowPtTripletStep,
     trackingMkFitDetachedQuadStep,
 #    trackingMkFitDetachedTripletStep, # to be enabled later
+#    trackingMkFitPixelPairStep,       # to be enabled later
 #    trackingMkFitMixedTripletStep,    # to be enabled later
     trackingMkFitPixelLessStep,
     trackingMkFitTobTecStep,

RecoTracker/IterativeTracking/python/InitialStepPreSplitting_cff.py

@@ -233,8 +233,12 @@
 _InitialStepPreSplittingTask_trackingPhase1.replace(initialStepHitTripletsPreSplitting, cms.Task(initialStepHitTripletsPreSplitting,initialStepHitQuadrupletsPreSplitting))
 trackingPhase1.toReplaceWith(InitialStepPreSplittingTask, _InitialStepPreSplittingTask_trackingPhase1.copyAndExclude([initialStepHitTripletsPreSplitting]))
 
+from Configuration.ProcessModifiers.trackingMkFitCommon_cff import trackingMkFitCommon
+_InitialStepPreSplittingTask_trackingMkFitCommon = InitialStepPreSplittingTask.copy()
+_InitialStepPreSplittingTask_trackingMkFitCommon.add(mkFitSiStripHits)
+trackingMkFitCommon.toReplaceWith(InitialStepPreSplittingTask, _InitialStepPreSplittingTask_trackingMkFitCommon)
 _InitialStepPreSplittingTask_trackingMkFit = InitialStepPreSplittingTask.copy()
-_InitialStepPreSplittingTask_trackingMkFit.add(mkFitSiPixelHitsPreSplitting, mkFitSiStripHits, mkFitEventOfHitsPreSplitting, initialStepTrackCandidatesMkFitSeedsPreSplitting, initialStepTrackCandidatesMkFitPreSplitting, initialStepTrackCandidatesMkFitConfigPreSplitting)
+_InitialStepPreSplittingTask_trackingMkFit.add(mkFitSiPixelHitsPreSplitting, mkFitEventOfHitsPreSplitting, initialStepTrackCandidatesMkFitSeedsPreSplitting, initialStepTrackCandidatesMkFitPreSplitting, initialStepTrackCandidatesMkFitConfigPreSplitting)
 trackingMkFitInitialStepPreSplitting.toReplaceWith(InitialStepPreSplittingTask, _InitialStepPreSplittingTask_trackingMkFit)
 
 

RecoTracker/IterativeTracking/python/PixelPairStep_cff.py

@@ -300,6 +300,30 @@
     onlyPixelHitsForSeedCleaner = cms.bool(True),
 
 )
+
+from Configuration.ProcessModifiers.trackingMkFitPixelPairStep_cff import trackingMkFitPixelPairStep
+import RecoTracker.MkFit.mkFitSeedConverter_cfi as _mkFitSeedConverter_cfi
+import RecoTracker.MkFit.mkFitIterationConfigESProducer_cfi as _mkFitIterationConfigESProducer_cfi
+import RecoTracker.MkFit.mkFitProducer_cfi as _mkFitProducer_cfi
+import RecoTracker.MkFit.mkFitOutputConverter_cfi as _mkFitOutputConverter_cfi
+pixelPairStepTrackCandidatesMkFitSeeds = _mkFitSeedConverter_cfi.mkFitSeedConverter.clone(
+    seeds = 'pixelPairStepSeeds',
+)
+pixelPairStepTrackCandidatesMkFitConfig = _mkFitIterationConfigESProducer_cfi.mkFitIterationConfigESProducer.clone(
+    ComponentName = 'pixelPairStepTrackCandidatesMkFitConfig',
+    config = 'RecoTracker/MkFit/data/mkfit-phase1-pixelPairStep.json',
+)
+pixelPairStepTrackCandidatesMkFit = _mkFitProducer_cfi.mkFitProducer.clone(
+    seeds = 'pixelPairStepTrackCandidatesMkFitSeeds',
+    config = ('', 'pixelPairStepTrackCandidatesMkFitConfig'),
+    clustersToSkip = 'pixelPairStepClusters',
+)
+trackingMkFitPixelPairStep.toReplaceWith(pixelPairStepTrackCandidates, _mkFitOutputConverter_cfi.mkFitOutputConverter.clone(
+    seeds = 'pixelPairStepSeeds',
+    mkFitSeeds = 'pixelPairStepTrackCandidatesMkFitSeeds',
+    tracks = 'pixelPairStepTrackCandidatesMkFit',
+))
+
 trackingPhase2PU140.toModify(pixelPairStepTrackCandidates,
     clustersToSkip       = None,
     phase2clustersToSkip = cms.InputTag('pixelPairStepClusters'),
@@ -435,6 +459,10 @@
                          pixelPairStep)
 PixelPairStep = cms.Sequence(PixelPairStepTask)
 
+_PixelPairStepTask_trackingMkFit = PixelPairStepTask.copy()
+_PixelPairStepTask_trackingMkFit.add(pixelPairStepTrackCandidatesMkFitSeeds, pixelPairStepTrackCandidatesMkFit, pixelPairStepTrackCandidatesMkFit)
+trackingMkFitPixelPairStep.toReplaceWith(PixelPairStepTask, _PixelPairStepTask_trackingMkFit)
+
 _PixelPairStepTask_LowPU_Phase2PU140 = PixelPairStepTask.copy()
 _PixelPairStepTask_LowPU_Phase2PU140.replace(pixelPairStep, pixelPairStepSelector)
 trackingLowPU.toReplaceWith(PixelPairStepTask, _PixelPairStepTask_LowPU_Phase2PU140)

RecoTracker/MkFit/README.md

@@ -39,6 +39,7 @@ $ runTheMatrix.py -l <workflow(s)> --apply 2 --command "--procModifiers tracking
 * *m_track_algorithm:* CMSSW track algorithm (used internally for reporting and consistency checks)
 * *m_requires_seed_hit_sorting:* do hits on seed tracks need to be sorted (required for seeds that include strip layers)
 * *m_require_quality_filter:* is additional post-processing required for result tracks
+* *m_require_dupclean_tight:* is tight duplicate removal post-processing required for result tracks
 * *m_params:* IterationParams structure for this iteration
 * *m_layer_configs:* std::vector of per-layer parameters
 
@@ -68,6 +69,9 @@ $ runTheMatrix.py -l <workflow(s)> --apply 2 --command "--procModifiers tracking
 
 * *minHitsQF:* min number of hits on track candidate to apply duplicate cleaning based on fraction of shared hits
 * *fracSharedHits:* min fraction of shared hits to determine duplicate track candidate
+* *drth_central:* dR cut used to identify duplicate candidates if std::abs(cotan(theta))<1.99 (abs(eta)<1.44)
+* *drth_obarrel:* dR cut used to identify duplicate candidates if 1.99<std::abs(cotan(theta))<6.05 (1.44<abs(eta)<2.5)
+* *drth_forward:* dR cut used to identify duplicate candidates if std::abs(cotan(theta))>6.05 (abs(eta)>2.5)
 
 ### Per-layer parameters [class IterationLayerConfig]
 

RecoTracker/MkFit/interface/MkFitGeometry.h

@@ -36,6 +36,7 @@ class MkFitGeometry {
   mkfit::IterationsInfo const& iterationsInfo() const { return *iterationsInfo_; }
   const std::vector<const DetLayer*>& detLayers() const { return dets_; }
   unsigned int uniqueIdInLayer(int layer, unsigned int detId) const { return detIdToShortId_.at(layer).at(detId); }
+  const TrackerTopology* topology() const { return ttopo_; }
 
 private:
   const TrackerTopology* ttopo_;

RecoTracker/MkFit/plugins/MkFitEventOfHitsProducer.cc

@@ -12,9 +12,15 @@
 #include "RecoTracker/MkFit/interface/MkFitHitWrapper.h"
 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
 
+#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
+#include "CalibTracker/Records/interface/SiStripQualityRcd.h"
+#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
+#include "DataFormats/TrackerCommon/interface/TrackerDetSide.h"
+
 // mkFit includes
 #include "mkFit/HitStructures.h"
 #include "mkFit/MkStdSeqs.h"
+#include "LayerNumberConverter.h"
 
 class MkFitEventOfHitsProducer : public edm::global::EDProducer<> {
 public:
@@ -35,7 +41,10 @@ class MkFitEventOfHitsProducer : public edm::global::EDProducer<> {
   const edm::EDGetTokenT<MkFitClusterIndexToHit> pixelClusterIndexToHitToken_;
   const edm::EDGetTokenT<MkFitClusterIndexToHit> stripClusterIndexToHitToken_;
   const edm::ESGetToken<MkFitGeometry, TrackerRecoGeometryRecord> mkFitGeomToken_;
+  edm::ESGetToken<SiStripQuality, SiStripQualityRcd> qualityToken_;
+  edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
   const edm::EDPutTokenT<MkFitEventOfHits> putToken_;
+  const bool useStripStripQualityDB_;
 };
 
 MkFitEventOfHitsProducer::MkFitEventOfHitsProducer(edm::ParameterSet const& iConfig)
@@ -44,13 +53,20 @@ MkFitEventOfHitsProducer::MkFitEventOfHitsProducer(edm::ParameterSet const& iCon
       pixelClusterIndexToHitToken_{consumes(iConfig.getParameter<edm::InputTag>("pixelHits"))},
       stripClusterIndexToHitToken_{consumes(iConfig.getParameter<edm::InputTag>("stripHits"))},
       mkFitGeomToken_{esConsumes()},
-      putToken_{produces<MkFitEventOfHits>()} {}
+      putToken_{produces<MkFitEventOfHits>()},
+      useStripStripQualityDB_{iConfig.getParameter<bool>("useStripStripQualityDB")} {
+  if (useStripStripQualityDB_) {
+    qualityToken_ = esConsumes();
+    geomToken_ = esConsumes();
+  }
+}
 
 void MkFitEventOfHitsProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   edm::ParameterSetDescription desc;
 
   desc.add("pixelHits", edm::InputTag{"mkFitSiPixelHits"});
   desc.add("stripHits", edm::InputTag{"mkFitSiStripHits"});
+  desc.add("useStripStripQualityDB", false)->setComment("Use SiStrip quality DB information");
 
   descriptions.addWithDefaultLabel(desc);
 }
@@ -63,6 +79,24 @@ void MkFitEventOfHitsProducer::produce(edm::StreamID iID, edm::Event& iEvent, co
   auto eventOfHits = std::make_unique<mkfit::EventOfHits>(mkFitGeom.trackerInfo());
   mkfit::StdSeq::Cmssw_LoadHits_Begin(*eventOfHits, {&pixelHits.hits(), &stripHits.hits()});
 
+  if (useStripStripQualityDB_) {
+    std::vector<mkfit::DeadVec> deadvectors(mkFitGeom.layerNumberConverter().nLayers());
+    const auto& siStripQuality = iSetup.getData(qualityToken_);
+    const auto& trackerGeom = iSetup.getData(geomToken_);
+    const auto& badStrips = siStripQuality.getBadComponentList();
+    for (const auto& bs : badStrips) {
+      const auto& surf = trackerGeom.idToDet(DetId(bs.detid))->surface();
+      const DetId detid(bs.detid);
+      bool isBarrel = (mkFitGeom.topology()->side(detid) == static_cast<unsigned>(TrackerDetSide::Barrel));
+      const auto ilay = mkFitGeom.mkFitLayerNumber(detid);
+      deadvectors[ilay].push_back({surf.phiSpan().first,
+                                   surf.phiSpan().second,
+                                   (isBarrel ? surf.zSpan().first : surf.rSpan().first),
+                                   (isBarrel ? surf.zSpan().second : surf.rSpan().second)});
+    }
+    mkfit::StdSeq::LoadDeads(*eventOfHits, deadvectors);
+  }
+
   fill(iEvent.get(pixelClusterIndexToHitToken_).hits(), *eventOfHits, mkFitGeom);
   fill(iEvent.get(stripClusterIndexToHitToken_).hits(), *eventOfHits, mkFitGeom);
 

RecoTracker/MkFit/plugins/MkFitOutputConverter.cc

@@ -272,46 +272,26 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp
     // seed
     const auto seedIndex = cand.label();
     LogTrace("MkFitOutputConverter") << " from seed " << seedIndex << " seed hits";
-    const auto& mkseed = mkFitSeeds.at(cand.label());
-    for (int i = 0; i < mkseed.nTotalHits(); ++i) {
-      const auto& hitOnTrack = mkseed.getHitOnTrack(i);
-      LogTrace("MkFitOutputConverter") << "  hit on layer " << hitOnTrack.layer << " index " << hitOnTrack.index;
-      // sanity check for now
-      const auto& candHitOnTrack = cand.getHitOnTrack(i);
-      if (hitOnTrack.layer != candHitOnTrack.layer) {
-        throw cms::Exception("LogicError")
-            << "Candidate " << candIndex << " from seed " << seedIndex << " hit " << i
-            << " has different layer in candidate (" << candHitOnTrack.layer << ") and seed (" << hitOnTrack.layer
-            << ")."
-            << " Hit indices are " << candHitOnTrack.index << " and " << hitOnTrack.index << ", respectively";
-      }
-      if (hitOnTrack.index != candHitOnTrack.index) {
-        throw cms::Exception("LogicError") << "Candidate " << candIndex << " from seed " << seedIndex << " hit " << i
-                                           << " has different hit index in candidate (" << candHitOnTrack.index
-                                           << ") and seed (" << hitOnTrack.index << ") on layer " << hitOnTrack.layer;
-      }
-    }
 
     // state
     auto state = cand.state();  // copy because have to modify
-    state.convertFromCCSToCartesian();
+    state.convertFromCCSToGlbCurvilinear();
     const auto& param = state.parameters;
     const auto& err = state.errors;
-    AlgebraicSymMatrix66 cov;
-    for (int i = 0; i < 6; ++i) {
-      for (int j = i; j < 6; ++j) {
+    AlgebraicSymMatrix55 cov;
+    for (int i = 0; i < 5; ++i) {
+      for (int j = i; j < 5; ++j) {
         cov[i][j] = err.At(i, j);
       }
     }
 
     auto fts = FreeTrajectoryState(
         GlobalTrajectoryParameters(
             GlobalPoint(param[0], param[1], param[2]), GlobalVector(param[3], param[4], param[5]), state.charge, &mf),
-        CartesianTrajectoryError(cov));
+        CurvilinearTrajectoryError(cov));
     if (!fts.curvilinearError().posDef()) {
       edm::LogWarning("MkFitOutputConverter") << "Curvilinear error not pos-def\n"
-                                              << fts.curvilinearError().matrix() << "\noriginal 6x6 covariance matrix\n"
-                                              << cov << "\ncandidate ignored";
+                                              << fts.curvilinearError().matrix() << "\ncandidate ignored";
       continue;
     }
 

RecoTracker/MkFit/plugins/MkFitSeedConverter.cc

@@ -111,17 +111,16 @@ mkfit::TrackVec MkFitSeedConverter::convertSeeds(const edm::View<TrajectorySeed>
     SVector3 pos(gpos.x(), gpos.y(), gpos.z());
     SVector3 mom(gmom.x(), gmom.y(), gmom.z());
 
-    const auto cartError = tsos.cartesianError();  // returns a temporary, so can't chain with the following line
-    const auto& cov = cartError.matrix();
-    SMatrixSym66 err;
-    for (int i = 0; i < 6; ++i) {
-      for (int j = i; j < 6; ++j) {
+    const auto& cov = tsos.curvilinearError().matrix();
+    SMatrixSym66 err;  //fill a sub-matrix, mkfit::TrackState will convert internally
+    for (int i = 0; i < 5; ++i) {
+      for (int j = i; j < 5; ++j) {
         err.At(i, j) = cov[i][j];
       }
     }
 
     mkfit::TrackState state(tsos.charge(), pos, mom, err);
-    state.convertFromCartesianToCCS();
+    state.convertFromGlbCurvilinearToCCS();
     ret.emplace_back(state, 0, seed_index, 0, nullptr);
     LogTrace("MkFitSeedConverter") << "Inserted seed with index " << seed_index;
 

RecoTracker/MkFit/plugins/createPhase1TrackerGeometry.cc

@@ -193,7 +193,7 @@ namespace mkfit {
     ti.set_eta_regions(0.9, 1.7, 2.45, false);
     ti.create_layers(18, 27, 27);
 
-    ii.resize(9);
+    ii.resize(10);
     ii[0].set_iteration_index_and_track_algorithm(0, (int)TrackBase::TrackAlgorithm::initialStep);
     ii[0].set_num_regions_layers(5, 72);
 
@@ -226,7 +226,10 @@ namespace mkfit {
     ii[8].Clone(ii[0]);
     ii[8].set_iteration_index_and_track_algorithm(8, (int)TrackBase::TrackAlgorithm::tobTecStep);
 
-    //for the latter 2 iter investing in maxCand & stop condition (for time) + QF and Dupl. cleaning (for quality)
+    ii[9].Clone(ii[0]);
+    ii[9].set_iteration_index_and_track_algorithm(9, (int)TrackBase::TrackAlgorithm::pixelPairStep);
+
+    // for iters [7] and [8]: investing in maxCand & stop condition (for time) + QF and Dupl. cleaning (for quality)
 
     // TODO: replace with MessageLogger
     if (verbose) {