Merge pull request #36546 from trackreco/mkFit-reincarnated

mkFit external integration

RecoTracker/MkFit/BuildFile.xml

@@ -3,7 +3,8 @@
 <use name="DataFormats/TrackerCommon"/>
 <use name="Geometry/TrackerGeometryBuilder"/>
 <use name="RecoTracker/TkDetLayers"/>
-<use name="mkfit"/>
+<use name="RecoTracker/MkFitCore"/>
+<use name="RecoTracker/MkFitCMS"/>
 <use name="rootcore"/>
 <export>
   <lib name="RecoTrackerMkFit"/>

RecoTracker/MkFit/plugins/BuildFile.xml

@@ -29,7 +29,6 @@
   <use name="TrackingTools/TrackFitters"/>
   <use name="TrackingTools/TrajectoryState"/>
   <use name="TrackingTools/TransientTrackingRecHit"/>
-  <use name="mkfit"/>
   <use name="rootmath"/>
   <flags EDM_PLUGIN="1"/>
 </library>

RecoTracker/MkFit/plugins/MkFitEventOfHitsProducer.cc

@@ -25,9 +25,9 @@
 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
 
 // mkFit includes
-#include "mkFit/HitStructures.h"
-#include "mkFit/MkStdSeqs.h"
-#include "LayerNumberConverter.h"
+#include "RecoTracker/MkFitCore/interface/HitStructures.h"
+#include "RecoTracker/MkFitCMS/interface/MkStdSeqs.h"
+#include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
 
 class MkFitEventOfHitsProducer : public edm::global::EDProducer<> {
 public:
@@ -95,7 +95,7 @@ void MkFitEventOfHitsProducer::produce(edm::StreamID iID, edm::Event& iEvent, co
   const auto& mkFitGeom = iSetup.getData(mkFitGeomToken_);
 
   auto eventOfHits = std::make_unique<mkfit::EventOfHits>(mkFitGeom.trackerInfo());
-  mkfit::StdSeq::Cmssw_LoadHits_Begin(*eventOfHits, {&pixelHits.hits(), &stripHits.hits()});
+  mkfit::StdSeq::cmssw_LoadHits_Begin(*eventOfHits, {&pixelHits.hits(), &stripHits.hits()});
 
   if (usePixelQualityDB_ || useStripStripQualityDB_) {
     std::vector<mkfit::DeadVec> deadvectors(mkFitGeom.layerNumberConverter().nLayers());
@@ -168,16 +168,16 @@ void MkFitEventOfHitsProducer::produce(edm::StreamID iID, edm::Event& iEvent, co
         }
       }
     }
-    mkfit::StdSeq::LoadDeads(*eventOfHits, deadvectors);
+    mkfit::StdSeq::loadDeads(*eventOfHits, deadvectors);
   }
 
   fill(iEvent.get(pixelClusterIndexToHitToken_).hits(), *eventOfHits, mkFitGeom);
   fill(iEvent.get(stripClusterIndexToHitToken_).hits(), *eventOfHits, mkFitGeom);
 
-  mkfit::StdSeq::Cmssw_LoadHits_End(*eventOfHits);
+  mkfit::StdSeq::cmssw_LoadHits_End(*eventOfHits);
 
   auto const bs = iEvent.get(beamSpotToken_);
-  eventOfHits->SetBeamSpot(
+  eventOfHits->setBeamSpot(
       mkfit::BeamSpot(bs.x0(), bs.y0(), bs.z0(), bs.sigmaZ(), bs.BeamWidthX(), bs.BeamWidthY(), bs.dxdz(), bs.dydz()));
 
   iEvent.emplace(putToken_, std::move(eventOfHits));
@@ -190,7 +190,7 @@ void MkFitEventOfHitsProducer::fill(const std::vector<const TrackingRecHit*>& hi
     const auto* hit = hits[i];
     if (hit != nullptr) {
       const auto ilay = mkFitGeom.mkFitLayerNumber(hit->geographicalId());
-      eventOfHits[ilay].RegisterHit(i);
+      eventOfHits[ilay].registerHit(i);
     }
   }
 }

RecoTracker/MkFit/plugins/MkFitGeometryESProducer.cc

@@ -13,9 +13,9 @@
 #include "createPhase1TrackerGeometry.h"
 
 // mkFit includes
-#include "ConfigWrapper.h"
-#include "TrackerInfo.h"
-#include "mkFit/IterationConfig.h"
+#include "RecoTracker/MkFitCore/interface/ConfigWrapper.h"
+#include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
+#include "RecoTracker/MkFitCore/interface/IterationConfig.h"
 
 #include <atomic>
 

RecoTracker/MkFit/plugins/MkFitIterationConfigESProducer.cc

@@ -6,10 +6,10 @@
 #include "RecoTracker/MkFit/interface/MkFitGeometry.h"
 
 // mkFit includes
-#include "Track.h"
-#include "TrackerInfo.h"
-#include "mkFit/HitStructures.h"
-#include "mkFit/IterationConfig.h"
+#include "RecoTracker/MkFitCore/interface/Track.h"
+#include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
+#include "RecoTracker/MkFitCore/interface/HitStructures.h"
+#include "RecoTracker/MkFitCore/interface/IterationConfig.h"
 
 namespace {
   using namespace mkfit;
@@ -26,7 +26,7 @@ namespace {
       const bool z_dir_pos = S.pz() > 0;
 
       const auto &hot = S.getLastHitOnTrack();
-      const float eta = eoh[hot.layer].GetHit(hot.index).eta();
+      const float eta = eoh[hot.layer].refHit(hot.index).eta();
 
       // Region to be defined by propagation / intersection tests
       TrackerInfo::EtaRegion reg;
@@ -42,39 +42,39 @@ namespace {
       constexpr int tecp1_id = 27;
       constexpr int tecn1_id = 54;
 
-      const LayerInfo &tib1 = trk_info.m_layers[tib1_id];
-      const LayerInfo &tob1 = trk_info.m_layers[tob1_id];
+      const LayerInfo &tib1 = trk_info.layer(tib1_id);
+      const LayerInfo &tob1 = trk_info.layer(tob1_id);
 
-      const LayerInfo &tecp1 = trk_info.m_layers[tecp1_id];
-      const LayerInfo &tecn1 = trk_info.m_layers[tecn1_id];
+      const LayerInfo &tecp1 = trk_info.layer(tecp1_id);
+      const LayerInfo &tecn1 = trk_info.layer(tecn1_id);
 
       const LayerInfo &tec_first = z_dir_pos ? tecp1 : tecn1;
 
       const float maxR = S.maxReachRadius();
       float z_at_maxr;
 
-      bool can_reach_outer_brl = S.canReachRadius(outer_brl.m_rout);
+      bool can_reach_outer_brl = S.canReachRadius(outer_brl.rout());
       float z_at_outer_brl;
       bool misses_first_tec;
       if (can_reach_outer_brl) {
-        z_at_outer_brl = S.zAtR(outer_brl.m_rout);
+        z_at_outer_brl = S.zAtR(outer_brl.rout());
         if (z_dir_pos)
-          misses_first_tec = z_at_outer_brl < tec_first.m_zmin;
+          misses_first_tec = z_at_outer_brl < tec_first.zmin();
         else
-          misses_first_tec = z_at_outer_brl > tec_first.m_zmax;
+          misses_first_tec = z_at_outer_brl > tec_first.zmax();
       } else {
         z_at_maxr = S.zAtR(maxR);
         if (z_dir_pos)
-          misses_first_tec = z_at_maxr < tec_first.m_zmin;
+          misses_first_tec = z_at_maxr < tec_first.zmin();
         else
-          misses_first_tec = z_at_maxr > tec_first.m_zmax;
+          misses_first_tec = z_at_maxr > tec_first.zmax();
       }
 
       if (misses_first_tec) {
         reg = TrackerInfo::Reg_Barrel;
       } else {
-        if ((S.canReachRadius(tib1.m_rin) && tib1.is_within_z_limits(S.zAtR(tib1.m_rin))) ||
-            (S.canReachRadius(tob1.m_rin) && tob1.is_within_z_limits(S.zAtR(tob1.m_rin)))) {
+        if ((S.canReachRadius(tib1.rin()) && tib1.is_within_z_limits(S.zAtR(tib1.rin()))) ||
+            (S.canReachRadius(tob1.rin()) && tob1.is_within_z_limits(S.zAtR(tob1.rin())))) {
           reg = z_dir_pos ? TrackerInfo::Reg_Transition_Pos : TrackerInfo::Reg_Transition_Neg;
         } else {
           reg = z_dir_pos ? TrackerInfo::Reg_Endcap_Pos : TrackerInfo::Reg_Endcap_Neg;
@@ -101,26 +101,26 @@ namespace {
     constexpr int tecp1_id = 27;
     constexpr int tecn1_id = 54;
 
-    const LayerInfo &tib1 = trk_info.m_layers[tib1_id];
-    const LayerInfo &tob1 = trk_info.m_layers[tob1_id];
+    const LayerInfo &tib1 = trk_info.layer(tib1_id);
+    const LayerInfo &tob1 = trk_info.layer(tob1_id);
 
-    const LayerInfo &tidp1 = trk_info.m_layers[tidp1_id];
-    const LayerInfo &tidn1 = trk_info.m_layers[tidn1_id];
+    const LayerInfo &tidp1 = trk_info.layer(tidp1_id);
+    const LayerInfo &tidn1 = trk_info.layer(tidn1_id);
 
-    const LayerInfo &tecp1 = trk_info.m_layers[tecp1_id];
-    const LayerInfo &tecn1 = trk_info.m_layers[tecn1_id];
+    const LayerInfo &tecp1 = trk_info.layer(tecp1_id);
+    const LayerInfo &tecn1 = trk_info.layer(tecn1_id);
 
     // Merge first two layers to account for mono/stereo coverage.
     // TrackerInfo could hold joint limits for sub-detectors.
-    const auto &L = trk_info.m_layers;
-    const float tidp_rin = std::min(L[tidp1_id].m_rin, L[tidp1_id + 1].m_rin);
-    const float tidp_rout = std::max(L[tidp1_id].m_rout, L[tidp1_id + 1].m_rout);
-    const float tecp_rin = std::min(L[tecp1_id].m_rin, L[tecp1_id + 1].m_rin);
-    const float tecp_rout = std::max(L[tecp1_id].m_rout, L[tecp1_id + 1].m_rout);
-    const float tidn_rin = std::min(L[tidn1_id].m_rin, L[tidn1_id + 1].m_rin);
-    const float tidn_rout = std::max(L[tidn1_id].m_rout, L[tidn1_id + 1].m_rout);
-    const float tecn_rin = std::min(L[tecn1_id].m_rin, L[tecn1_id + 1].m_rin);
-    const float tecn_rout = std::max(L[tecn1_id].m_rout, L[tecn1_id + 1].m_rout);
+    const auto &L = trk_info;
+    const float tidp_rin = std::min(L[tidp1_id].rin(), L[tidp1_id + 1].rin());
+    const float tidp_rout = std::max(L[tidp1_id].rout(), L[tidp1_id + 1].rout());
+    const float tecp_rin = std::min(L[tecp1_id].rin(), L[tecp1_id + 1].rin());
+    const float tecp_rout = std::max(L[tecp1_id].rout(), L[tecp1_id + 1].rout());
+    const float tidn_rin = std::min(L[tidn1_id].rin(), L[tidn1_id + 1].rin());
+    const float tidn_rout = std::max(L[tidn1_id].rout(), L[tidn1_id + 1].rout());
+    const float tecn_rin = std::min(L[tecn1_id].rin(), L[tecn1_id + 1].rin());
+    const float tecn_rout = std::max(L[tecn1_id].rout(), L[tecn1_id + 1].rout());
 
     // Bias towards more aggressive transition-region assignemnts.
     // With current tunning it seems to make things a bit worse.
@@ -153,7 +153,7 @@ namespace {
       const Track &S = in_seeds[i];
 
       const auto &hot = S.getLastHitOnTrack();
-      const float eta = eoh[hot.layer].GetHit(hot.index).eta();
+      const float eta = eoh[hot.layer].refHit(hot.index).eta();
 
       // Region to be defined by propagation / intersection tests
       TrackerInfo::EtaRegion reg;
@@ -165,14 +165,14 @@ namespace {
       const float maxR = S.maxReachRadius();
 
       if (z_dir_pos) {
-        const bool in_tib = barrel_pos_check(S, maxR, tib1.m_rin, tib1.m_zmax);
-        const bool in_tob = barrel_pos_check(S, maxR, tob1.m_rin, tob1.m_zmax);
+        const bool in_tib = barrel_pos_check(S, maxR, tib1.rin(), tib1.zmax());
+        const bool in_tob = barrel_pos_check(S, maxR, tob1.rin(), tob1.zmax());
 
         if (!in_tib && !in_tob) {
           reg = TrackerInfo::Reg_Endcap_Pos;
         } else {
-          const bool in_tid = endcap_pos_check(S, maxR, tidp_rout, tidp_rin, tidp1.m_zmin - tid_z_extra);
-          const bool in_tec = endcap_pos_check(S, maxR, tecp_rout, tecp_rin, tecp1.m_zmin - tec_z_extra);
+          const bool in_tid = endcap_pos_check(S, maxR, tidp_rout, tidp_rin, tidp1.zmin() - tid_z_extra);
+          const bool in_tec = endcap_pos_check(S, maxR, tecp_rout, tecp_rin, tecp1.zmin() - tec_z_extra);
 
           if (!in_tid && !in_tec) {
             reg = TrackerInfo::Reg_Barrel;
@@ -181,14 +181,14 @@ namespace {
           }
         }
       } else {
-        const bool in_tib = barrel_neg_check(S, maxR, tib1.m_rin, tib1.m_zmin);
-        const bool in_tob = barrel_neg_check(S, maxR, tob1.m_rin, tob1.m_zmin);
+        const bool in_tib = barrel_neg_check(S, maxR, tib1.rin(), tib1.zmin());
+        const bool in_tob = barrel_neg_check(S, maxR, tob1.rin(), tob1.zmin());
 
         if (!in_tib && !in_tob) {
           reg = TrackerInfo::Reg_Endcap_Neg;
         } else {
-          const bool in_tid = endcap_neg_check(S, maxR, tidn_rout, tidn_rin, tidn1.m_zmax + tid_z_extra);
-          const bool in_tec = endcap_neg_check(S, maxR, tecn_rout, tecn_rin, tecn1.m_zmax + tec_z_extra);
+          const bool in_tid = endcap_neg_check(S, maxR, tidn_rout, tidn_rin, tidn1.zmax() + tid_z_extra);
+          const bool in_tec = endcap_neg_check(S, maxR, tecn_rout, tecn_rin, tecn1.zmax() + tec_z_extra);
 
           if (!in_tid && !in_tec) {
             reg = TrackerInfo::Reg_Barrel;
@@ -233,7 +233,8 @@ void MkFitIterationConfigESProducer::fillDescriptions(edm::ConfigurationDescript
 
 std::unique_ptr<mkfit::IterationConfig> MkFitIterationConfigESProducer::produce(
     const TrackerRecoGeometryRecord &iRecord) {
-  auto it_conf = mkfit::ConfigJson_Load_File(configFile_);
+  mkfit::ConfigJson cj;
+  auto it_conf = cj.load_File(configFile_);
   it_conf->m_partition_seeds = partitionSeeds1;
   return it_conf;
 }

RecoTracker/MkFit/plugins/MkFitOutputConverter.cc

@@ -37,9 +37,9 @@
 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
 
 // mkFit indludes
-#include "LayerNumberConverter.h"
-#include "Track.h"
-#include "mkFit/HitStructures.h"
+#include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
+#include "RecoTracker/MkFitCore/interface/Track.h"
+#include "RecoTracker/MkFitCore/interface/HitStructures.h"
 
 namespace {
   template <typename T>

RecoTracker/MkFit/plugins/MkFitProducer.cc

@@ -18,11 +18,11 @@
 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
 
 // mkFit includes
-#include "ConfigWrapper.h"
-#include "LayerNumberConverter.h"
-#include "mkFit/buildtestMPlex.h"
-#include "mkFit/IterationConfig.h"
-#include "mkFit/MkBuilderWrapper.h"
+#include "RecoTracker/MkFitCore/interface/ConfigWrapper.h"
+#include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
+#include "RecoTracker/MkFitCMS/interface/runFunctions.h"
+#include "RecoTracker/MkFitCore/interface/IterationConfig.h"
+#include "RecoTracker/MkFitCore/interface/MkBuilderWrapper.h"
 
 // TBB includes
 #include "oneapi/tbb/task_arena.h"
@@ -54,7 +54,6 @@ class MkFitProducer : public edm::global::EDProducer<edm::StreamCache<mkfit::MkB
   const edm::ESGetToken<MkFitGeometry, TrackerRecoGeometryRecord> mkFitGeomToken_;
   const edm::ESGetToken<mkfit::IterationConfig, TrackerRecoGeometryRecord> mkFitIterConfigToken_;
   const edm::EDPutTokenT<MkFitOutputWrapper> putToken_;
-  std::function<double(mkfit::Event&, mkfit::MkBuilder&)> buildFunction_;
   const float minGoodStripCharge_;
   const bool seedCleaning_;
   const bool backwardFitInCMSSW_;
@@ -101,7 +100,7 @@ MkFitProducer::MkFitProducer(edm::ParameterSet const& iConfig)
 
   // TODO: what to do when we have multiple instances of MkFitProducer in a job?
   mkfit::MkBuilderWrapper::populate();
-  mkfit::ConfigWrapper::initializeForCMSSW(mkFitSilent_);
+  mkfit::ConfigWrapper::initializeForCMSSW();
 }
 
 void MkFitProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
@@ -134,12 +133,9 @@ void MkFitProducer::fillDescriptions(edm::ConfigurationDescriptions& description
 }
 
 std::unique_ptr<mkfit::MkBuilderWrapper> MkFitProducer::beginStream(edm::StreamID iID) const {
-  return std::make_unique<mkfit::MkBuilderWrapper>();
+  return std::make_unique<mkfit::MkBuilderWrapper>(mkFitSilent_);
 }
 
-namespace {
-  std::once_flag geometryFlag;
-}
 void MkFitProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
   const auto& pixelHits = iEvent.get(pixelHitsToken_);
   const auto& stripHits = iEvent.get(stripHitsToken_);
@@ -180,13 +176,6 @@ void MkFitProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::Ev
     stripClusterChargeCut(iEvent.get(stripClusterChargeToken_), stripMask);
   }
 
-  // Initialize the number of layers, has to be done exactly once in
-  // the whole program.
-  // TODO: the mechanism needs to be improved...
-  std::call_once(geometryFlag, [nlayers = mkFitGeom.layerNumberConverter().nLayers()]() {
-    mkfit::ConfigWrapper::setNTotalLayers(nlayers);
-  });
-
   // seeds need to be mutable because of the possible cleaning
   auto seeds_mutable = seeds.seeds();
   mkfit::TrackVec tracks;

RecoTracker/MkFit/plugins/MkFitSeedConverter.cc

@@ -31,8 +31,8 @@
 #include "Math/SMatrix.h"
 
 // mkFit includes
-#include "LayerNumberConverter.h"
-#include "Track.h"
+#include "RecoTracker/MkFitCMS/interface/LayerNumberConverter.h"
+#include "RecoTracker/MkFitCore/interface/Track.h"
 
 class MkFitSeedConverter : public edm::global::EDProducer<> {
 public:

RecoTracker/MkFit/plugins/convertHits.h

@@ -17,8 +17,8 @@
 #include "Math/SMatrix.h"
 
 // mkFit includes
-#include "Hit.h"
-#include "mkFit/HitStructures.h"
+#include "RecoTracker/MkFitCore/interface/Hit.h"
+#include "RecoTracker/MkFitCore/interface/HitStructures.h"
 
 namespace mkfit {
   template <typename Traits, typename HitCollection>

RecoTracker/MkFit/plugins/createPhase1TrackerGeometry.cc

@@ -2,9 +2,8 @@
 // Phase1 tracker geometry
 //-------------------
 
-#include "Config.h"
-#include "Debug.h"
-#include "TrackerInfo.h"
+#include "RecoTracker/MkFitCore/interface/Config.h"
+#include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
 
 #include <functional>
 
@@ -16,20 +15,6 @@ namespace {
 
 namespace mkfit {
   void createPhase1TrackerGeometry(TrackerInfo &ti, bool verbose) {
-    // TODO: these writes to global variables need to be removed
-    Config::nTotalLayers = 18 + 2 * 27;
-
-    Config::useCMSGeom = true;
-
-    Config::finding_requires_propagation_to_hit_pos = true;
-    Config::finding_inter_layer_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
-    Config::finding_intra_layer_pflags = PropagationFlags(PF_none);
-    Config::backward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
-    Config::forward_fit_pflags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
-    Config::seed_fit_pflags = PropagationFlags(PF_none);
-    Config::pca_prop_pflags = PropagationFlags(PF_none);
-
-    ti.set_eta_regions(0.9, 1.7, 2.45, false);
     ti.create_layers(18, 27, 27);
     createPhase1TrackerGeometryAutoGen(ti);
 
@@ -38,8 +23,8 @@ namespace mkfit {
       printf("==========================================================================================\n");
       printf("Phase1 tracker -- Create_TrackerInfo finished\n");
       printf("==========================================================================================\n");
-      for (auto &i : ti.m_layers)
-        i.print_layer();
+      for (int ii = 0; ii < ti.n_layers(); ++ii)
+        ti.layer(ii).print_layer();
       printf("==========================================================================================\n");
     }
   }