Add producers for mkFit

Also include customization functions to easily include it in matrix
workflows. In this context amend the list of initialStep tracking
modules for MTV timing plots.

RecoTracker/MkFit/BuildFile.xml

@@ -0,0 +1,8 @@
+<use name="DataFormats/Provenance"/>
+<use name="FWCore/Framework"/>
+<use name="FWCore/MessageLogger"/>
+<use name="mkfit"/>
+<use name="rootcore"/>
+<export>
+  <lib name="RecoTrackerMkFit"/>
+</export>

RecoTracker/MkFit/README.md

@@ -0,0 +1,27 @@
+# mkFit
+
+This package holds the glue modules for running
+[mkFit](http://trackreco.github.io/) within CMSSW.
+
+Note that at the moment there may be only one `MkFitProducer` in a
+single job. This restriction will be removed in the future.
+
+## Customize functions for runTheMatrix workflows (offline reconstruction)
+
+* `RecoTracker/MkFit/customizeInitialStepToMkFit.customizeInitialStepToMkFit`
+  * Replaces initialStep track building module with `mkFit`.
+* `RecoTracker/MkFit/customizeInitialStepOnly.customizeInitialStepOnly`
+  * Run only the initialStep tracking. In practice this configuration
+    runs the initialStepPreSplitting iteration, but named as
+    initialStep. MultiTrackValidator is included, and configured to
+    monitor initialStep. Intended to provide the minimal configuration
+    for CMSSW tests.
+* `RecoTracker/MkFit/customizeInitialStepOnly.customizeInitialStepOnlyNoMTV`
+  * Otherwise same as `customizeInitialStepOnly` except drops
+    MultiTrackValidator. Intended for profiling.
+
+
+These can be used with e.g.
+```bash
+$ runTheMatrix.py -l <workflow(s)> --apply 2 --command "--customise RecoTracker/MkFit/customizeInitialStepToMkFit.customizeInitialStepToMkFit"
+```

RecoTracker/MkFit/interface/MkFitIndexLayer.h

@@ -0,0 +1,45 @@
+#ifndef RecoTracker_MkFit_MkFitIndexLayer_h
+#define RecoTracker_MkFit_MkFitIndexLayer_h
+
+#include "DataFormats/Provenance/interface/ProductID.h"
+
+#include <vector>
+
+class TrackingRecHit;
+
+class MkFitIndexLayer {
+public:
+  struct HitInfo {
+    HitInfo() : index(-1), layer(-1) {}
+    HitInfo(int i, int l) : index(i), layer(l) {}
+    int index;
+    int layer;
+  };
+
+  struct Coll {
+    explicit Coll(edm::ProductID id) : productID(id) {}
+    edm::ProductID productID;
+    std::vector<HitInfo> infos;  // indexed by cluster index
+  };
+
+  MkFitIndexLayer() = default;
+
+  void insert(edm::ProductID id, size_t clusterIndex, int hit, int layer, const TrackingRecHit *hitPtr);
+
+  const HitInfo &get(edm::ProductID id, size_t clusterIndex) const;
+
+  const TrackingRecHit *getHitPtr(int layer, int hit) const { return hits_.at(layer).at(hit).ptr; }
+
+  size_t getClusterIndex(int layer, int hit) const { return hits_.at(layer).at(hit).clusterIndex; }
+
+private:
+  struct CMSSWHit {
+    const TrackingRecHit *ptr = nullptr;
+    size_t clusterIndex = 0;
+  };
+
+  std::vector<Coll> colls_;                   // mapping from CMSSW(ProductID, index) -> mkfit(index, layer)
+  std::vector<std::vector<CMSSWHit> > hits_;  // reverse mapping, mkfit(layer, index) -> CMSSW hit
+};
+
+#endif

RecoTracker/MkFit/interface/MkFitInputWrapper.h

@@ -0,0 +1,44 @@
+#ifndef RecoTracker_MkFit_MkFitInputWrapper_h
+#define RecoTracker_MkFit_MkFitInputWrapper_h
+
+#include "RecoTracker/MkFit/interface/MkFitIndexLayer.h"
+
+#include <memory>
+#include <vector>
+
+namespace mkfit {
+  class Hit;
+  class Track;
+  class LayerNumberConverter;
+  using HitVec = std::vector<Hit>;
+  using TrackVec = std::vector<Track>;
+}  // namespace mkfit
+
+class MkFitInputWrapper {
+public:
+  MkFitInputWrapper();
+  MkFitInputWrapper(MkFitIndexLayer&& indexLayers,
+                    std::vector<mkfit::HitVec>&& hits,
+                    mkfit::TrackVec&& seeds,
+                    mkfit::LayerNumberConverter&& lnc);
+  ~MkFitInputWrapper();
+
+  MkFitInputWrapper(MkFitInputWrapper const&) = delete;
+  MkFitInputWrapper& operator=(MkFitInputWrapper const&) = delete;
+  MkFitInputWrapper(MkFitInputWrapper&&) = default;
+  MkFitInputWrapper& operator=(MkFitInputWrapper&&) = default;
+
+  MkFitIndexLayer const& indexLayers() const { return indexLayers_; }
+  mkfit::TrackVec const& seeds() const { return *seeds_; }
+  std::vector<mkfit::HitVec> const& hits() const { return hits_; }
+  mkfit::LayerNumberConverter const& layerNumberConverter() const { return *lnc_; }
+  unsigned int nlayers() const;
+
+private:
+  MkFitIndexLayer indexLayers_;
+  std::vector<mkfit::HitVec> hits_;
+  std::unique_ptr<mkfit::TrackVec> seeds_;            // for pimpl pattern
+  std::unique_ptr<mkfit::LayerNumberConverter> lnc_;  // for pimpl pattern
+};
+
+#endif

RecoTracker/MkFit/interface/MkFitOutputWrapper.h

@@ -0,0 +1,30 @@
+#ifndef RecoTracker_MkFit_MkFitOutputWrapper_h
+#define RecoTracker_MkFit_MkFitOutputWrapper_h
+
+#include <vector>
+
+namespace mkfit {
+  class Track;
+  using TrackVec = std::vector<Track>;
+}  // namespace mkfit
+
+class MkFitOutputWrapper {
+public:
+  MkFitOutputWrapper();
+  MkFitOutputWrapper(mkfit::TrackVec&& candidateTracks, mkfit::TrackVec&& fitTracks);
+  ~MkFitOutputWrapper();
+
+  MkFitOutputWrapper(MkFitOutputWrapper const&) = delete;
+  MkFitOutputWrapper& operator=(MkFitOutputWrapper const&) = delete;
+  MkFitOutputWrapper(MkFitOutputWrapper&&) = default;
+  MkFitOutputWrapper& operator=(MkFitOutputWrapper&&) = default;
+
+  mkfit::TrackVec const& candidateTracks() const { return candidateTracks_; }
+  mkfit::TrackVec const& fitTracks() const { return fitTracks_; }
+
+private:
+  mkfit::TrackVec candidateTracks_;
+  mkfit::TrackVec fitTracks_;
+};
+
+#endif

RecoTracker/MkFit/plugins/BuildFile.xml

@@ -0,0 +1,27 @@
+<library file="*.cc" name="RecoTrackerMkFitPlugins">
+  <use name="FWCore/Framework"/>
+  <use name="FWCore/MessageLogger"/>
+  <use name="FWCore/ParameterSet"/> 
+  <use name="FWCore/PluginManager"/>
+  <use name="DataFormats/TrajectorySeed"/>
+  <use name="DataFormats/TrackerRecHit2D"/>
+  <use name="DataFormats/TrackCandidate"/>
+  <use name="DataFormats/TrackReco"/>
+  <use name="DataFormats/TrackerCommon"/>
+  <use name="Geometry/Records"/>
+  <use name="Geometry/TrackerGeometryBuilder"/>
+  <use name="RecoTracker/TransientTrackingRecHit"/>
+  <use name="TrackingTools/GeomPropagators"/>
+  <use name="TrackingTools/KalmanUpdators"/>
+  <use name="TrackingTools/MaterialEffects"/>
+  <use name="RecoTracker/MeasurementDet"/>
+  <use name="TrackingTools/Records"/>
+  <use name="RecoTracker/TkDetLayers"/>
+  <use name="TrackingTools/TrackFitters"/>
+  <use name="TrackingTools/TrajectoryState"/>
+  <use name="TrackingTools/TransientTrackingRecHit"/>
+  <use name="RecoTracker/MkFit"/>
+  <use name="rootmath"/>
+  <use name="mkfit"/>
+  <flags EDM_PLUGIN="1"/>
+</library>

RecoTracker/MkFit/plugins/MkFitInputConverter.cc

@@ -0,0 +1,211 @@
+#include "FWCore/Framework/interface/global/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
+#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
+#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
+#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
+
+#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "Geometry/Records/interface/TrackerTopologyRcd.h"
+
+#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
+#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
+#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
+
+#include "MagneticField/Engine/interface/MagneticField.h"
+#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
+
+#include "RecoTracker/MkFit/interface/MkFitInputWrapper.h"
+
+// ROOT
+#include "Math/SVector.h"
+#include "Math/SMatrix.h"
+
+// MkFit includes
+#include "Hit.h"
+#include "Track.h"
+#include "LayerNumberConverter.h"
+
+class MkFitInputConverter : public edm::global::EDProducer<> {
+public:
+  explicit MkFitInputConverter(edm::ParameterSet const& iConfig);
+  ~MkFitInputConverter() override = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
+
+  template <typename HitCollection>
+  void convertHits(const HitCollection& hits,
+                   std::vector<mkfit::HitVec>& mkfitHits,
+                   MkFitIndexLayer& indexLayers,
+                   int& totalHits,
+                   const TrackerTopology& ttopo,
+                   const TransientTrackingRecHitBuilder& ttrhBuilder,
+                   const mkfit::LayerNumberConverter& lnc) const;
+
+  bool passCCC(const SiStripRecHit2D& hit, const DetId hitId) const;
+  bool passCCC(const SiPixelRecHit& hit, const DetId hitId) const;
+
+  mkfit::TrackVec convertSeeds(const edm::View<TrajectorySeed>& seeds,
+                               const MkFitIndexLayer& indexLayers,
+                               const TransientTrackingRecHitBuilder& ttrhBuilder,
+                               const MagneticField& mf) const;
+
+  using SVector3 = ROOT::Math::SVector<float, 3>;
+  using SMatrixSym33 = ROOT::Math::SMatrix<float, 3, 3, ROOT::Math::MatRepSym<float, 3>>;
+  using SMatrixSym66 = ROOT::Math::SMatrix<float, 6, 6, ROOT::Math::MatRepSym<float, 6>>;
+
+  edm::EDGetTokenT<SiPixelRecHitCollection> pixelRecHitToken_;
+  edm::EDGetTokenT<SiStripRecHit2DCollection> stripRphiRecHitToken_;
+  edm::EDGetTokenT<SiStripRecHit2DCollection> stripStereoRecHitToken_;
+  edm::EDGetTokenT<edm::View<TrajectorySeed>> seedToken_;
+  edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> ttrhBuilderToken_;
+  edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttopoToken_;
+  edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mfToken_;
+  edm::EDPutTokenT<MkFitInputWrapper> putToken_;
+};
+
+MkFitInputConverter::MkFitInputConverter(edm::ParameterSet const& iConfig)
+    : pixelRecHitToken_{consumes<SiPixelRecHitCollection>(iConfig.getParameter<edm::InputTag>("pixelRecHits"))},
+      stripRphiRecHitToken_{
+          consumes<SiStripRecHit2DCollection>(iConfig.getParameter<edm::InputTag>("stripRphiRecHits"))},
+      stripStereoRecHitToken_{
+          consumes<SiStripRecHit2DCollection>(iConfig.getParameter<edm::InputTag>("stripStereoRecHits"))},
+      seedToken_{consumes<edm::View<TrajectorySeed>>(iConfig.getParameter<edm::InputTag>("seeds"))},
+      ttrhBuilderToken_{esConsumes<TransientTrackingRecHitBuilder, TransientRecHitRecord>(
+          iConfig.getParameter<edm::ESInputTag>("ttrhBuilder"))},
+      ttopoToken_{esConsumes<TrackerTopology, TrackerTopologyRcd>()},
+      mfToken_{esConsumes<MagneticField, IdealMagneticFieldRecord>()},
+      putToken_{produces<MkFitInputWrapper>()} {}
+
+void MkFitInputConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+
+  desc.add("pixelRecHits", edm::InputTag{"siPixelRecHits"});
+  desc.add("stripRphiRecHits", edm::InputTag{"siStripMatchedRecHits", "rphiRecHit"});
+  desc.add("stripStereoRecHits", edm::InputTag{"siStripMatchedRecHits", "stereoRecHit"});
+  desc.add("seeds", edm::InputTag{"initialStepSeeds"});
+  desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"});
+
+  descriptions.addWithDefaultLabel(desc);
+}
+
+void MkFitInputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
+  mkfit::LayerNumberConverter lnc{mkfit::TkLayout::phase1};
+
+  // Then import hits
+  const auto& ttrhBuilder = iSetup.getData(ttrhBuilderToken_);
+  const auto& ttopo = iSetup.getData(ttopoToken_);
+
+  std::vector<mkfit::HitVec> mkfitHits(lnc.nLayers());
+  MkFitIndexLayer indexLayers;
+  int totalHits = 0;  // I need to have a global hit index in order to have the hit remapping working?
+  convertHits(iEvent.get(pixelRecHitToken_), mkfitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc);
+  convertHits(iEvent.get(stripRphiRecHitToken_), mkfitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc);
+  convertHits(iEvent.get(stripStereoRecHitToken_), mkfitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc);
+
+  // Then import seeds
+  auto mkfitSeeds = convertSeeds(iEvent.get(seedToken_), indexLayers, ttrhBuilder, iSetup.getData(mfToken_));
+
+  iEvent.emplace(putToken_, std::move(indexLayers), std::move(mkfitHits), std::move(mkfitSeeds), std::move(lnc));
+}
+
+bool MkFitInputConverter::passCCC(const SiStripRecHit2D& hit, const DetId hitId) const {
+  return (siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster()) >= 1620);
+}
+
+bool MkFitInputConverter::passCCC(const SiPixelRecHit& hit, const DetId hitId) const { return true; }
+
+template <typename HitCollection>
+void MkFitInputConverter::convertHits(const HitCollection& hits,
+                                      std::vector<mkfit::HitVec>& mkfitHits,
+                                      MkFitIndexLayer& indexLayers,
+                                      int& totalHits,
+                                      const TrackerTopology& ttopo,
+                                      const TransientTrackingRecHitBuilder& ttrhBuilder,
+                                      const mkfit::LayerNumberConverter& lnc) const {
+  for (const auto& detset : hits) {
+    const DetId detid = detset.detId();
+    const auto subdet = detid.subdetId();
+    const auto layer = ttopo.layer(detid);
+    const auto isStereo = ttopo.isStereo(detid);
+
+    for (const auto& hit : detset) {
+      if (!passCCC(hit, detid))
+        continue;
+
+      const auto& gpos = hit.globalPosition();
+      SVector3 pos(gpos.x(), gpos.y(), gpos.z());
+      const auto& gerr = hit.globalPositionError();
+      SMatrixSym33 err;
+      err.At(0, 0) = gerr.cxx();
+      err.At(1, 1) = gerr.cyy();
+      err.At(2, 2) = gerr.czz();
+      err.At(0, 1) = gerr.cyx();
+      err.At(0, 2) = gerr.czx();
+      err.At(1, 2) = gerr.czy();
+
+      const auto ilay = lnc.convertLayerNumber(subdet, layer, false, isStereo, gpos.z() > 0);
+      LogTrace("MkFitInputConverter") << "Adding hit detid " << detid.rawId() << " subdet " << subdet << " layer "
+                                      << layer << " isStereo " << isStereo << " zplus " << (gpos.z() > 0) << " ilay "
+                                      << ilay;
+
+      indexLayers.insert(hit.firstClusterRef().id(), hit.firstClusterRef().index(), mkfitHits[ilay].size(), ilay, &hit);
+      mkfitHits[ilay].emplace_back(pos, err, totalHits);
+      ++totalHits;
+    }
+  }
+}
+
+mkfit::TrackVec MkFitInputConverter::convertSeeds(const edm::View<TrajectorySeed>& seeds,
+                                                  const MkFitIndexLayer& indexLayers,
+                                                  const TransientTrackingRecHitBuilder& ttrhBuilder,
+                                                  const MagneticField& mf) const {
+  mkfit::TrackVec ret;
+  ret.reserve(seeds.size());
+  int index = 0;
+  for (const auto& seed : seeds) {
+    const auto hitRange = seed.recHits();
+    const auto lastRecHit = ttrhBuilder.build(&*(hitRange.second - 1));
+    const auto tsos = trajectoryStateTransform::transientState(seed.startingState(), lastRecHit->surface(), &mf);
+    const auto& stateGlobal = tsos.globalParameters();
+    const auto& gpos = stateGlobal.position();
+    const auto& gmom = stateGlobal.momentum();
+    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) {
+        err.At(i, j) = cov[i][j];
+      }
+    }
+
+    mkfit::TrackState state(tsos.charge(), pos, mom, err);
+    state.convertFromCartesianToCCS();
+    ret.emplace_back(state, 0, index, 0, nullptr);
+
+    // Add hits
+    for (auto iHit = hitRange.first; iHit != hitRange.second; ++iHit) {
+      const auto* hit = dynamic_cast<const BaseTrackerRecHit*>(&*iHit);
+      if (hit == nullptr) {
+        throw cms::Exception("Assert") << "Encountered a seed with a hit which is not BaseTrackerRecHit";
+      }
+
+      const auto& info = indexLayers.get(hit->firstClusterRef().id(), hit->firstClusterRef().index());
+      ret.back().addHitIdx(info.index, info.layer, 0);  // per-hit chi2 is not known
+    }
+    ++index;
+  }
+  return ret;
+}
+
+DEFINE_FWK_MODULE(MkFitInputConverter);