From 0bd86e7a46f4232bcfa356a892e8ea390cc89a71 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 21 Mar 2018 16:27:24 +0100 Subject: [PATCH 01/20] 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 | 8 + RecoTracker/MkFit/README.md | 27 + RecoTracker/MkFit/interface/MkFitIndexLayer.h | 45 ++ .../MkFit/interface/MkFitInputWrapper.h | 44 ++ .../MkFit/interface/MkFitOutputWrapper.h | 30 + RecoTracker/MkFit/plugins/BuildFile.xml | 27 + .../MkFit/plugins/MkFitInputConverter.cc | 211 +++++++ .../MkFit/plugins/MkFitOutputConverter.cc | 580 ++++++++++++++++++ RecoTracker/MkFit/plugins/MkFitProducer.cc | 148 +++++ .../MkFit/python/customizeInitialStepOnly.py | 146 +++++ .../python/customizeInitialStepToMkFit.py | 19 + RecoTracker/MkFit/src/MkFitIndexLayer.cc | 45 ++ RecoTracker/MkFit/src/MkFitInputWrapper.cc | 20 + RecoTracker/MkFit/src/MkFitOutputWrapper.cc | 11 + RecoTracker/MkFit/src/classes.h | 8 + RecoTracker/MkFit/src/classes_def.xml | 7 + .../python/plotting/trackingPlots.py | 3 + 17 files changed, 1379 insertions(+) create mode 100644 RecoTracker/MkFit/BuildFile.xml create mode 100644 RecoTracker/MkFit/README.md create mode 100644 RecoTracker/MkFit/interface/MkFitIndexLayer.h create mode 100644 RecoTracker/MkFit/interface/MkFitInputWrapper.h create mode 100644 RecoTracker/MkFit/interface/MkFitOutputWrapper.h create mode 100644 RecoTracker/MkFit/plugins/BuildFile.xml create mode 100644 RecoTracker/MkFit/plugins/MkFitInputConverter.cc create mode 100644 RecoTracker/MkFit/plugins/MkFitOutputConverter.cc create mode 100644 RecoTracker/MkFit/plugins/MkFitProducer.cc create mode 100644 RecoTracker/MkFit/python/customizeInitialStepOnly.py create mode 100644 RecoTracker/MkFit/python/customizeInitialStepToMkFit.py create mode 100644 RecoTracker/MkFit/src/MkFitIndexLayer.cc create mode 100644 RecoTracker/MkFit/src/MkFitInputWrapper.cc create mode 100644 RecoTracker/MkFit/src/MkFitOutputWrapper.cc create mode 100644 RecoTracker/MkFit/src/classes.h create mode 100644 RecoTracker/MkFit/src/classes_def.xml diff --git a/RecoTracker/MkFit/BuildFile.xml b/RecoTracker/MkFit/BuildFile.xml new file mode 100644 index 0000000000000..46570ba3ef629 --- /dev/null +++ b/RecoTracker/MkFit/BuildFile.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/RecoTracker/MkFit/README.md b/RecoTracker/MkFit/README.md new file mode 100644 index 0000000000000..7bbfde5274be7 --- /dev/null +++ b/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 --apply 2 --command "--customise RecoTracker/MkFit/customizeInitialStepToMkFit.customizeInitialStepToMkFit" +``` \ No newline at end of file diff --git a/RecoTracker/MkFit/interface/MkFitIndexLayer.h b/RecoTracker/MkFit/interface/MkFitIndexLayer.h new file mode 100644 index 0000000000000..44035ecb17918 --- /dev/null +++ b/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 + +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 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 colls_; // mapping from CMSSW(ProductID, index) -> mkfit(index, layer) + std::vector > hits_; // reverse mapping, mkfit(layer, index) -> CMSSW hit +}; + +#endif diff --git a/RecoTracker/MkFit/interface/MkFitInputWrapper.h b/RecoTracker/MkFit/interface/MkFitInputWrapper.h new file mode 100644 index 0000000000000..a281b66c3e415 --- /dev/null +++ b/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 +#include + +namespace mkfit { + class Hit; + class Track; + class LayerNumberConverter; + using HitVec = std::vector; + using TrackVec = std::vector; +} // namespace mkfit + +class MkFitInputWrapper { +public: + MkFitInputWrapper(); + MkFitInputWrapper(MkFitIndexLayer&& indexLayers, + std::vector&& 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 const& hits() const { return hits_; } + mkfit::LayerNumberConverter const& layerNumberConverter() const { return *lnc_; } + unsigned int nlayers() const; + +private: + MkFitIndexLayer indexLayers_; + std::vector hits_; + std::unique_ptr seeds_; // for pimpl pattern + std::unique_ptr lnc_; // for pimpl pattern +}; + +#endif diff --git a/RecoTracker/MkFit/interface/MkFitOutputWrapper.h b/RecoTracker/MkFit/interface/MkFitOutputWrapper.h new file mode 100644 index 0000000000000..d03a4af250556 --- /dev/null +++ b/RecoTracker/MkFit/interface/MkFitOutputWrapper.h @@ -0,0 +1,30 @@ +#ifndef RecoTracker_MkFit_MkFitOutputWrapper_h +#define RecoTracker_MkFit_MkFitOutputWrapper_h + +#include + +namespace mkfit { + class Track; + using TrackVec = std::vector; +} // 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 diff --git a/RecoTracker/MkFit/plugins/BuildFile.xml b/RecoTracker/MkFit/plugins/BuildFile.xml new file mode 100644 index 0000000000000..facea9b6356da --- /dev/null +++ b/RecoTracker/MkFit/plugins/BuildFile.xml @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc new file mode 100644 index 0000000000000..d6e70527bc6ee --- /dev/null +++ b/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 + void convertHits(const HitCollection& hits, + std::vector& 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& seeds, + const MkFitIndexLayer& indexLayers, + const TransientTrackingRecHitBuilder& ttrhBuilder, + const MagneticField& mf) const; + + using SVector3 = ROOT::Math::SVector; + using SMatrixSym33 = ROOT::Math::SMatrix>; + using SMatrixSym66 = ROOT::Math::SMatrix>; + + edm::EDGetTokenT pixelRecHitToken_; + edm::EDGetTokenT stripRphiRecHitToken_; + edm::EDGetTokenT stripStereoRecHitToken_; + edm::EDGetTokenT> seedToken_; + edm::ESGetToken ttrhBuilderToken_; + edm::ESGetToken ttopoToken_; + edm::ESGetToken mfToken_; + edm::EDPutTokenT putToken_; +}; + +MkFitInputConverter::MkFitInputConverter(edm::ParameterSet const& iConfig) + : pixelRecHitToken_{consumes(iConfig.getParameter("pixelRecHits"))}, + stripRphiRecHitToken_{ + consumes(iConfig.getParameter("stripRphiRecHits"))}, + stripStereoRecHitToken_{ + consumes(iConfig.getParameter("stripStereoRecHits"))}, + seedToken_{consumes>(iConfig.getParameter("seeds"))}, + ttrhBuilderToken_{esConsumes( + iConfig.getParameter("ttrhBuilder"))}, + ttopoToken_{esConsumes()}, + mfToken_{esConsumes()}, + putToken_{produces()} {} + +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 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 +void MkFitInputConverter::convertHits(const HitCollection& hits, + std::vector& 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& 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(&*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); diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc new file mode 100644 index 0000000000000..4df7d011cd0d4 --- /dev/null +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -0,0 +1,580 @@ +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Utilities/interface/do_nothing_deleter.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" +#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h" +#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h" +#include "DataFormats/TrackReco/interface/SeedStopInfo.h" +#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h" + +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.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 "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" +#include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h" +#include "TrackingTools/MaterialEffects/src/PropagatorWithMaterial.cc" + +#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" +#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h" + +#include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" +#include "RecoTracker/MkFit/interface/MkFitOutputWrapper.h" + +// MkFit indludes +#include "LayerNumberConverter.h" +#include "Track.h" + +namespace { + template + bool isBarrel(T subdet) { + return subdet == PixelSubdetector::PixelBarrel || subdet == StripSubdetector::TIB || + subdet == StripSubdetector::TOB; + } + + template + bool isEndcap(T subdet) { + return subdet == PixelSubdetector::PixelEndcap || subdet == StripSubdetector::TID || + subdet == StripSubdetector::TEC; + } +} // namespace + +class MkFitOutputConverter : public edm::global::EDProducer<> { +public: + explicit MkFitOutputConverter(edm::ParameterSet const& iConfig); + ~MkFitOutputConverter() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; + + std::vector createDetLayers(const mkfit::LayerNumberConverter& lnc, + const GeometricSearchTracker& tracker, + const TrackerTopology& ttopo) const; + + TrackCandidateCollection convertCandidates(const MkFitOutputWrapper& mkfitOutput, + const MkFitIndexLayer& indexLayers, + const edm::View& seeds, + const TrackerGeometry& geom, + const MagneticField& mf, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite, + const TkClonerImpl& hitCloner, + const std::vector& detLayers, + const mkfit::TrackVec& mkfitSeeds) const; + + std::pair backwardFit(const FreeTrajectoryState& fts, + const edm::OwnVector& hits, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite, + const TkClonerImpl& hitCloner, + bool lastHitWasInvalid, + bool lastHitWasChanged) const; + + std::pair backwardFitImpl( + const FreeTrajectoryState& fts, + const TransientTrackingRecHit::ConstRecHitContainer& firstHits, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite, + const TkClonerImpl& hitCloner, + bool lastHitWasInvalid, + bool lastHitWasChanged) const; + + std::pair convertInnermostState(const FreeTrajectoryState& fts, + const edm::OwnVector& hits, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite) const; + + edm::EDGetTokenT hitsSeedsToken_; + edm::EDGetTokenT tracksToken_; + edm::EDGetTokenT> seedToken_; + edm::EDGetTokenT mteToken_; + edm::ESGetToken geomToken_; + edm::ESGetToken propagatorAlongToken_; + edm::ESGetToken propagatorOppositeToken_; + edm::ESGetToken ttopoToken_; + edm::ESGetToken mfToken_; + edm::ESGetToken ttrhBuilderToken_; + edm::EDPutTokenT putTrackCandidateToken_; + edm::EDPutTokenT> putSeedStopInfoToken_; + std::string ttrhBuilderName_; + std::string propagatorAlongName_; + std::string propagatorOppositeName_; + bool backwardFitInCMSSW_; +}; + +MkFitOutputConverter::MkFitOutputConverter(edm::ParameterSet const& iConfig) + : hitsSeedsToken_{consumes(iConfig.getParameter("hitsSeeds"))}, + tracksToken_{consumes(iConfig.getParameter("tracks"))}, + seedToken_{consumes>(iConfig.getParameter("seeds"))}, + mteToken_{consumes(iConfig.getParameter("measurementTrackerEvent"))}, + geomToken_{esConsumes()}, + propagatorAlongToken_{ + esConsumes(iConfig.getParameter("propagatorAlong"))}, + propagatorOppositeToken_{esConsumes( + iConfig.getParameter("propagatorOpposite"))}, + ttopoToken_{esConsumes()}, + mfToken_{esConsumes()}, + ttrhBuilderToken_{esConsumes( + iConfig.getParameter("ttrhBuilder"))}, + putTrackCandidateToken_{produces()}, + putSeedStopInfoToken_{produces>()}, + backwardFitInCMSSW_{iConfig.getParameter("backwardFitInCMSSW")} {} + +void MkFitOutputConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + + desc.add("hitsSeeds", edm::InputTag{"mkFitInputConverter"}); + desc.add("tracks", edm::InputTag{"mkfitProducer"}); + desc.add("seeds", edm::InputTag{"initialStepSeeds"}); + desc.add("measurementTrackerEvent", edm::InputTag{"MeasurementTrackerEvent"}); + desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"}); + desc.add("propagatorAlong", edm::ESInputTag{"", "PropagatorWithMaterial"}); + desc.add("propagatorOpposite", edm::ESInputTag{"", "PropagatorWithMaterialOpposite"}); + desc.add("backwardFitInCMSSW", false) + ->setComment("Do backward fit (to innermost hit) in CMSSW (true) or mkFit (false)"); + + descriptions.addWithDefaultLabel(desc); +} + +void MkFitOutputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + const auto& seeds = iEvent.get(seedToken_); + const auto& hitsSeeds = iEvent.get(hitsSeedsToken_); + const auto& mte = iEvent.get(mteToken_); + + const auto& ttrhBuilder = iSetup.getData(ttrhBuilderToken_); + const auto* tkBuilder = dynamic_cast(&ttrhBuilder); + if (!tkBuilder) { + throw cms::Exception("LogicError") << "TTRHBuilder must be of type TkTransientTrackingRecHitBuilder"; + } + + // Convert mkfit presentation back to CMSSW + const auto detlayers = + createDetLayers(hitsSeeds.layerNumberConverter(), *(mte.geometricSearchTracker()), iSetup.getData(ttopoToken_)); + iEvent.emplace(putTrackCandidateToken_, + convertCandidates(iEvent.get(tracksToken_), + hitsSeeds.indexLayers(), + seeds, + iSetup.getData(geomToken_), + iSetup.getData(mfToken_), + iSetup.getData(propagatorAlongToken_), + iSetup.getData(propagatorOppositeToken_), + tkBuilder->cloner(), + detlayers, + hitsSeeds.seeds())); + + // TODO: SeedStopInfo is currently unfilled + iEvent.emplace(putSeedStopInfoToken_, seeds.size()); +} + +std::vector MkFitOutputConverter::createDetLayers(const mkfit::LayerNumberConverter& lnc, + const GeometricSearchTracker& tracker, + const TrackerTopology& ttopo) const { + std::vector dets(lnc.nLayers(), nullptr); + + auto set = [&](unsigned int index, DetId id) { + auto layer = tracker.idToLayer(id); + if (layer == nullptr) { + throw cms::Exception("LogicError") << "No layer for DetId " << id.rawId(); + } + LogTrace("MkFitOutputConverter") << "Setting DetLayer for index " << index << " subdet " << id.subdetId() + << " layer " << ttopo.layer(id) << " ptr " << layer; + + dets[index] = layer; + }; + + // TODO: currently hardcoded... + // Logic copied from mkfit::LayerNumberConverter + unsigned int off = 1; + // BPix + set(0, ttopo.pxbDetId(1, 0, 0)); + set(1, ttopo.pxbDetId(2, 0, 0)); + set(2, ttopo.pxbDetId(3, 0, 0)); + set(3, ttopo.pxbDetId(4, 0, 0)); + // TIB + set(off + 3, ttopo.tibDetId(1, 0, 0, 0, 0, 0)); + set(off + 4, ttopo.tibDetId(1, 0, 0, 0, 0, 1)); + set(off + 5, ttopo.tibDetId(2, 0, 0, 0, 0, 0)); + set(off + 6, ttopo.tibDetId(2, 0, 0, 0, 0, 1)); + set(off + 7, ttopo.tibDetId(3, 0, 0, 0, 0, 0)); + set(off + 8, ttopo.tibDetId(4, 0, 0, 0, 0, 0)); + // TOB + set(off + 9, ttopo.tobDetId(1, 0, 0, 0, 0)); + set(off + 10, ttopo.tobDetId(1, 0, 0, 0, 1)); + set(off + 11, ttopo.tobDetId(2, 0, 0, 0, 0)); + set(off + 12, ttopo.tobDetId(2, 0, 0, 0, 1)); + set(off + 13, ttopo.tobDetId(3, 0, 0, 0, 0)); + set(off + 14, ttopo.tobDetId(4, 0, 0, 0, 0)); + set(off + 15, ttopo.tobDetId(5, 0, 0, 0, 0)); + set(off + 16, ttopo.tobDetId(6, 0, 0, 0, 0)); + + auto setForward = [&](unsigned int side) { + // FPix + set(off + 0, ttopo.pxfDetId(side, 1, 0, 0, 0)); + set(off + 1, ttopo.pxfDetId(side, 2, 0, 0, 0)); + set(off + 2, ttopo.pxfDetId(side, 3, 0, 0, 0)); + // TID+ + off += 1; + set(off + 2, ttopo.tidDetId(side, 1, 0, 0, 0, 0)); + set(off + 3, ttopo.tidDetId(side, 1, 0, 0, 0, 1)); + set(off + 4, ttopo.tidDetId(side, 2, 0, 0, 0, 0)); + set(off + 5, ttopo.tidDetId(side, 2, 0, 0, 0, 1)); + set(off + 6, ttopo.tidDetId(side, 3, 0, 0, 0, 0)); + set(off + 7, ttopo.tidDetId(side, 3, 0, 0, 0, 1)); + // TEC + set(off + 8, ttopo.tecDetId(side, 1, 0, 0, 0, 0, 0)); + set(off + 9, ttopo.tecDetId(side, 1, 0, 0, 0, 0, 1)); + set(off + 10, ttopo.tecDetId(side, 2, 0, 0, 0, 0, 0)); + set(off + 11, ttopo.tecDetId(side, 2, 0, 0, 0, 0, 1)); + set(off + 12, ttopo.tecDetId(side, 3, 0, 0, 0, 0, 0)); + set(off + 13, ttopo.tecDetId(side, 3, 0, 0, 0, 0, 1)); + set(off + 14, ttopo.tecDetId(side, 4, 0, 0, 0, 0, 0)); + set(off + 15, ttopo.tecDetId(side, 4, 0, 0, 0, 0, 1)); + set(off + 16, ttopo.tecDetId(side, 5, 0, 0, 0, 0, 0)); + set(off + 17, ttopo.tecDetId(side, 5, 0, 0, 0, 0, 1)); + set(off + 18, ttopo.tecDetId(side, 6, 0, 0, 0, 0, 0)); + set(off + 19, ttopo.tecDetId(side, 6, 0, 0, 0, 0, 1)); + set(off + 20, ttopo.tecDetId(side, 7, 0, 0, 0, 0, 0)); + set(off + 21, ttopo.tecDetId(side, 7, 0, 0, 0, 0, 1)); + set(off + 22, ttopo.tecDetId(side, 8, 0, 0, 0, 0, 0)); + set(off + 23, ttopo.tecDetId(side, 8, 0, 0, 0, 0, 1)); + set(off + 24, ttopo.tecDetId(side, 9, 0, 0, 0, 0, 0)); + set(off + 25, ttopo.tecDetId(side, 9, 0, 0, 0, 0, 1)); + }; + + // plus + off = 17 + 1; + setForward(2); + + // minus + off = 17 + 1 + 25 + 2; + setForward(1); + + return dets; +} + +TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutputWrapper& mkfitOutput, + const MkFitIndexLayer& indexLayers, + const edm::View& seeds, + const TrackerGeometry& geom, + const MagneticField& mf, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite, + const TkClonerImpl& hitCloner, + const std::vector& detLayers, + const mkfit::TrackVec& mkfitSeeds) const { + TrackCandidateCollection output; + const auto& candidates = backwardFitInCMSSW_ ? mkfitOutput.candidateTracks() : mkfitOutput.fitTracks(); + output.reserve(candidates.size()); + + LogTrace("MkFitOutputConverter") << "Number of candidates " << mkfitOutput.candidateTracks().size(); + + int candIndex = -1; + for (const auto& cand : candidates) { + ++candIndex; + LogTrace("MkFitOutputConverter") << "Candidate " << candIndex << " pT " << cand.pT() << " eta " << cand.momEta() + << " phi " << cand.momPhi() << " chi2 " << cand.chi2(); + + // hits + edm::OwnVector recHits; + const int nhits = cand.nTotalHits(); // what exactly is the difference between nTotalHits() and nFoundHits()? + bool lastHitInvalid = false; + for (int i = 0; i < nhits; ++i) { + const auto& hitOnTrack = cand.getHitOnTrack(i); + LogTrace("MkFitOutputConverter") << " hit on layer " << hitOnTrack.layer << " index " << hitOnTrack.index; + if (hitOnTrack.index < 0) { + // What is the exact meaning of -1, -2, -3? + // In order to use the regular InvalidTrackingRecHit I'd need + // a GeomDet (and "unfortunately" that is needed in + // TrackProducer). + // + // I guess we could take the track state and propagate it to + // each layer to find the actual module the track crosses, and + // check whether it is active or not to be able to mark + // inactive hits + const auto* detLayer = detLayers.at(hitOnTrack.layer); + if (detLayer == nullptr) { + throw cms::Exception("LogicError") << "DetLayer for layer index " << hitOnTrack.layer << " is null!"; + } + // Actually it is necessary to leave dealing with invalid hits to the TrackProducer? + //recHits.push_back(new InvalidTrackingRecHitNoDet(detLayer->surface(), TrackingRecHit::missing)); // let's put them all as missing for now + lastHitInvalid = true; + } else { + recHits.push_back(indexLayers.getHitPtr(hitOnTrack.layer, hitOnTrack.index)->clone()); + LogTrace("MkFitOutputConverter") << " pos " << recHits.back().globalPosition().x() << " " + << recHits.back().globalPosition().y() << " " + << recHits.back().globalPosition().z() << " mag2 " + << recHits.back().globalPosition().mag2() << " detid " + << recHits.back().geographicalId().rawId() << " cluster " + << indexLayers.getClusterIndex(hitOnTrack.layer, hitOnTrack.index); + lastHitInvalid = false; + } + } + + const auto lastHitId = recHits.back().geographicalId(); + + // MkFit hits are *not* in the order of propagation, sort by 3D radius for now (as we don't have loopers) + // TODO: Improve the sorting (extract keys? maybe even bubble sort would work well as the hits are almost in the correct order) + recHits.sort([](const auto& a, const auto& b) { + const auto aid = a.geographicalId(); + const auto bid = b.geographicalId(); + + const auto asub = aid.subdetId(); + const auto bsub = bid.subdetId(); + if (asub != bsub) { + // Subdetector order (BPix, FPix, TIB, TID, TOB, TEC) corresponds also the navigation + return asub < bsub; + } + + const auto& apos = a.globalPosition(); + const auto& bpos = b.globalPosition(); + + if (isBarrel(asub)) { + return apos.perp2() < bpos.perp2(); + } + return std::abs(apos.z()) < std::abs(bpos.z()); + }); + + const bool lastHitChanged = (recHits.back().geographicalId() != lastHitId); // TODO: make use of the bools + + // 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") + edm::LogError("MkFitOutputConverter") + << "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") + edm::LogError("MkFitOutputConverter") + << "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(); + 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) { + 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)); + if (!fts.curvilinearError().posDef()) { + edm::LogWarning("MkFitOutputConverter") << "Curvilinear error not pos-def\n" + << fts.curvilinearError().matrix() << "\noriginal 6x6 covariance matrix\n" + << cov << "\ncandidate ignored"; + continue; + } + + auto tsosDet = + backwardFitInCMSSW_ + ? backwardFit(fts, recHits, propagatorAlong, propagatorOpposite, hitCloner, lastHitInvalid, lastHitChanged) + : convertInnermostState(fts, recHits, propagatorAlong, propagatorOpposite); + if (!tsosDet.first.isValid()) { + edm::LogWarning("MkFitOutputConverter") + << "Backward fit of candidate " << candIndex << " failed, ignoring the candidate"; + continue; + } + + // convert to persistent, from CkfTrackCandidateMakerBase + auto pstate = trajectoryStateTransform::persistentState(tsosDet.first, tsosDet.second->geographicalId().rawId()); + + output.emplace_back( + recHits, + seeds.at(seedIndex), + pstate, + seeds.refAt(seedIndex), + 0, // TODO: nloops, let's ignore for now + static_cast(StopReason::UNINITIALIZED) // TODO: ignore details of stopping reason as well for now + ); + } + return output; +} + +std::pair MkFitOutputConverter::backwardFit( + const FreeTrajectoryState& fts, + const edm::OwnVector& hits, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite, + const TkClonerImpl& hitCloner, + bool lastHitWasInvalid, + bool lastHitWasChanged) const { + // First filter valid hits as in TransientInitialStateEstimator + TransientTrackingRecHit::ConstRecHitContainer firstHits; + + for (int i = hits.size() - 1; i >= 0; --i) { + if (hits[i].det()) { + // TransientTrackingRecHit::ConstRecHitContainer has shared_ptr, + // and it is passed to backFitter below so it is really needed + // to keep the interface. Since we keep the ownership in hits, + // let's disable the deleter. + firstHits.emplace_back(&(hits[i]), edm::do_nothing_deleter{}); + } + } + + auto ret = backwardFitImpl( + fts, firstHits, propagatorAlong, propagatorOpposite, hitCloner, lastHitWasInvalid, lastHitWasChanged); + return ret; +} + +std::pair MkFitOutputConverter::backwardFitImpl( + const FreeTrajectoryState& fts, + const TransientTrackingRecHit::ConstRecHitContainer& firstHits, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite, + const TkClonerImpl& hitCloner, + bool lastHitWasInvalid, + bool lastHitWasChanged) const { + // Then propagate along to the surface of the last hit to get a TSOS + const auto& lastHitSurface = firstHits.front()->det()->surface(); + + const Propagator* tryFirst = &propagatorAlong; + const Propagator* trySecond = &propagatorOpposite; + if (lastHitWasInvalid || lastHitWasChanged) { + LogTrace("MkFitOutputConverter") << "Propagating first opposite, then along, because lastHitWasInvalid? " + << lastHitWasInvalid << " or lastHitWasChanged? " << lastHitWasChanged; + std::swap(tryFirst, trySecond); + } else { + const auto lastHitSubdet = firstHits.front()->geographicalId().subdetId(); + const auto& surfacePos = lastHitSurface.position(); + const auto& lastHitPos = firstHits.front()->globalPosition(); + bool doSwitch = false; + if (isBarrel(lastHitSubdet)) { + doSwitch = (surfacePos.perp2() < lastHitPos.perp2()); + } else { + doSwitch = (surfacePos.z() < lastHitPos.z()); + } + if (doSwitch) { + LogTrace("MkFitOutputConverter") + << "Propagating first opposite, then along, because surface is inner than the hit; surface perp2 " + << surfacePos.perp() << " hit " << lastHitPos.perp2() << " surface z " << surfacePos.z() << " hit " + << lastHitPos.z(); + + std::swap(tryFirst, trySecond); + } + } + + auto tsosDouble = tryFirst->propagateWithPath(fts, lastHitSurface); + if (!tsosDouble.first.isValid()) { + LogDebug("MkFitOutputConverter") << "Propagating to startingState failed, trying in another direction next"; + tsosDouble = trySecond->propagateWithPath(fts, lastHitSurface); + } + auto& startingState = tsosDouble.first; + + if (!startingState.isValid()) { + edm::LogWarning("MkFitOutputConverter") + << "startingState is not valid, FTS was\n" + << fts << " last hit surface surface:" + << "\n position " << lastHitSurface.position() << "\n phiSpan " << lastHitSurface.phiSpan().first << "," + << lastHitSurface.phiSpan().first << "\n rSpan " << lastHitSurface.rSpan().first << "," + << lastHitSurface.rSpan().first << "\n zSpan " << lastHitSurface.zSpan().first << "," + << lastHitSurface.zSpan().first; + return std::pair(); + } + + // Then return back to the logic from TransientInitialStateEstimator + startingState.rescaleError(100.); + + // avoid cloning + KFUpdator const aKFUpdator; + Chi2MeasurementEstimator const aChi2MeasurementEstimator(100., 3); + KFTrajectoryFitter backFitter( + &propagatorAlong, &aKFUpdator, &aChi2MeasurementEstimator, firstHits.size(), nullptr, &hitCloner); + + PropagationDirection backFitDirection = + oppositeToMomentum; // assume for now that the propagation in mkfit always alongMomentum + + // only direction matters in this contest + TrajectorySeed fakeSeed(PTrajectoryStateOnDet(), edm::OwnVector(), backFitDirection); + + Trajectory fitres = + backFitter.fitOne(fakeSeed, firstHits, startingState, TrajectoryFitter::standard); // ignore loopers for now + + LogDebug("MkFitOutputConverter") << "using a backward fit of :" << firstHits.size() << " hits, starting from:\n" + << startingState << " to get the estimate of the initial state of the track."; + + if (!fitres.isValid()) { + edm::LogWarning("MkFitOutputConverter") << "FitTester: first hits fit failed"; + return std::pair(); + } + + TrajectoryMeasurement const& firstMeas = fitres.lastMeasurement(); + + // magnetic field can be different! + TrajectoryStateOnSurface firstState(firstMeas.updatedState().localParameters(), + firstMeas.updatedState().localError(), + firstMeas.updatedState().surface(), + propagatorAlong.magneticField()); + + firstState.rescaleError(100.); + + LogDebug("MkFitOutputConverter") << "the initial state is found to be:\n:" << firstState + << "\n it's field pointer is: " << firstState.magneticField() + << "\n the pointer from the state of the back fit was: " + << firstMeas.updatedState().magneticField(); + + return std::make_pair(firstState, firstMeas.recHit()->det()); +} + +std::pair MkFitOutputConverter::convertInnermostState( + const FreeTrajectoryState& fts, + const edm::OwnVector& hits, + const Propagator& propagatorAlong, + const Propagator& propagatorOpposite) const { + auto det = hits[0].det(); + if (det == nullptr) { + throw cms::Exception("LogicError") << "Got nullptr from the first hit det()"; + } + + const auto& firstHitSurface = det->surface(); + + auto tsosDouble = propagatorAlong.propagateWithPath(fts, firstHitSurface); + if (!tsosDouble.first.isValid()) { + LogDebug("MkFitOutputConverter") << "Propagating to startingState along momentum failed, trying opposite next"; + tsosDouble = propagatorOpposite.propagateWithPath(fts, firstHitSurface); + } + + //return std::make_pair(TrajectoryStateOnSurface(fts, det->surface()), det); + return std::make_pair(tsosDouble.first, det); +} + +DEFINE_FWK_MODULE(MkFitOutputConverter); diff --git a/RecoTracker/MkFit/plugins/MkFitProducer.cc b/RecoTracker/MkFit/plugins/MkFitProducer.cc new file mode 100644 index 0000000000000..06c40228da85e --- /dev/null +++ b/RecoTracker/MkFit/plugins/MkFitProducer.cc @@ -0,0 +1,148 @@ +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" + +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" + +#include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" +#include "RecoTracker/MkFit/interface/MkFitOutputWrapper.h" + +// MkFit includes +#include "ConfigWrapper.h" +#include "Event.h" +#include "mkFit/buildtestMPlex.h" +#include "mkFit/MkBuilderWrapper.h" + +// TBB includes +#include "tbb/task_arena.h" + +// std includes +#include + +class MkFitProducer : public edm::global::EDProducer > { +public: + explicit MkFitProducer(edm::ParameterSet const& iConfig); + ~MkFitProducer() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + std::unique_ptr beginStream(edm::StreamID) const override; + +private: + void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; + + edm::EDGetTokenT hitsSeedsToken_; + edm::ESGetToken geomToken_; + edm::EDPutTokenT putToken_; + std::function buildFunction_; + bool backwardFitInCMSSW_; + bool mkfitSilent_; +}; + +MkFitProducer::MkFitProducer(edm::ParameterSet const& iConfig) + : hitsSeedsToken_{consumes(iConfig.getParameter("hitsSeeds"))}, + geomToken_{esConsumes()}, + putToken_{produces()}, + backwardFitInCMSSW_{iConfig.getParameter("backwardFitInCMSSW")}, + mkfitSilent_{iConfig.getUntrackedParameter("mkfitSilent")} { + const auto build = iConfig.getParameter("buildingRoutine"); + bool isFV = false; + if (build == "bestHit") { + buildFunction_ = mkfit::runBuildingTestPlexBestHit; + } else if (build == "standard") { + buildFunction_ = mkfit::runBuildingTestPlexStandard; + } else if (build == "cloneEngine") { + buildFunction_ = mkfit::runBuildingTestPlexCloneEngine; + } else if (build == "fullVector") { + isFV = true; + buildFunction_ = mkfit::runBuildingTestPlexFV; + } else { + throw cms::Exception("Configuration") << "Invalid value for parameter 'buildingRoutine' " << build + << ", allowed are bestHit, standard, cloneEngine, fullVector"; + } + + const auto seedClean = iConfig.getParameter("seedCleaning"); + auto seedCleanOpt = mkfit::ConfigWrapper::SeedCleaningOpts::noCleaning; + if (seedClean == "none") { + seedCleanOpt = mkfit::ConfigWrapper::SeedCleaningOpts::noCleaning; + } else if (seedClean == "N2") { + seedCleanOpt = mkfit::ConfigWrapper::SeedCleaningOpts::cleanSeedsN2; + } else { + throw cms::Exception("Configuration") + << "Invalida value for parameter 'seedCleaning' " << seedClean << ", allowed are none, N2"; + } + + auto backwardFitOpt = + backwardFitInCMSSW_ ? mkfit::ConfigWrapper::BackwardFit::noFit : mkfit::ConfigWrapper::BackwardFit::toFirstLayer; + + // TODO: what to do when we have multiple instances of MkFitProducer in a job? + mkfit::MkBuilderWrapper::populate(isFV); + mkfit::ConfigWrapper::initializeForCMSSW(seedCleanOpt, backwardFitOpt, mkfitSilent_); +} + +void MkFitProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + + desc.add("hitsSeeds", edm::InputTag("mkFitInputConverter")); + desc.add("buildingRoutine", "cloneEngine") + ->setComment("Valid values are: 'bestHit', 'standard', 'cloneEngine', 'fullVector'"); + desc.add("seedCleaning", "N2")->setComment("Valid values are: 'none', 'N2'"); + desc.add("backwardFitInCMSSW", false) + ->setComment("Do backward fit (to innermost hit) in CMSSW (true) or mkFit (false)"); + desc.addUntracked("mkfitSilent", true)->setComment("Allows to enables printouts from mkfit with 'False'"); + + descriptions.add("mkFitProducer", desc); +} + +std::unique_ptr MkFitProducer::beginStream(edm::StreamID iID) const { + return std::make_unique(); +} + +namespace { + std::once_flag geometryFlag; +} +void MkFitProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + const auto& hitsSeeds = iEvent.get(hitsSeedsToken_); + const auto& geom = iSetup.getData(geomToken_); + + // First do some initialization + // TODO: the mechanism needs to be improved... + std::call_once(geometryFlag, [&geom, nlayers = hitsSeeds.nlayers()]() { + // TODO: eventually automatize fully + // For now it is easier to use purely the infrastructure from mkfit + /* + const auto barrelLayers = geom.numberOfLayers(PixelSubdetector::PixelBarrel) + geom.numberOfLayers(StripSubdetector::TIB) + geom.numberOfLayers(StripSubdetector::TOB); + const auto endcapLayers = geom.numberOfLayers(PixelSubdetector::PixelEndcap) + geom.numberOfLayers(StripSubdetector::TID) + geom.numberOfLayers(StripSubdetector::TEC); + // TODO: Number of stereo layers is hardcoded for now, so this won't work for phase2 tracker + const auto barrelStereo = 2 + 2; + const auto endcapStereo = geom.numberOfLayers(StripSubdetector::TID) + geom.numberOfLayers(StripSubdetector::TEC); + const auto nlayers = barrelLayers + barrelStereo + 2*(endcapLayers + endcapStereo); + LogDebug("MkFitProducer") << "Total number of tracker layers (stereo counted separately) " << nlayers; + */ + if (geom.numberOfLayers(PixelSubdetector::PixelBarrel) != 4 || + geom.numberOfLayers(PixelSubdetector::PixelEndcap) != 3) { + throw cms::Exception("Assert") << "For now this code works only with phase1 tracker, you have something else"; + } + mkfit::ConfigWrapper::setNTotalLayers(nlayers); + }); + + // CMSSW event ID (64-bit unsigned) does not fit in int + // In addition, unique ID requires also lumi and run + // But does the event ID really matter within mkfit? + mkfit::Event ev(iEvent.id().event()); + + ev.setInputFromCMSSW(hitsSeeds.hits(), hitsSeeds.seeds()); + + tbb::this_task_arena::isolate([&]() { buildFunction_(ev, streamCache(iID)->get()); }); + + iEvent.emplace(putToken_, std::move(ev.candidateTracks_), std::move(ev.fitTracks_)); +} + +DEFINE_FWK_MODULE(MkFitProducer); diff --git a/RecoTracker/MkFit/python/customizeInitialStepOnly.py b/RecoTracker/MkFit/python/customizeInitialStepOnly.py new file mode 100644 index 0000000000000..5a7fdff65e571 --- /dev/null +++ b/RecoTracker/MkFit/python/customizeInitialStepOnly.py @@ -0,0 +1,146 @@ +import FWCore.ParameterSet.Config as cms + +def removePath(process, pname): + if hasattr(process, pname): + process.schedule.remove(getattr(process, pname)) + delattr(process, pname) + +def customizeInitialStepOnly(process): + # Customize reconstruction + process.trackerClusterCheck.PixelClusterCollectionLabel = 'siPixelClustersPreSplitting' + process.initialStepSeedLayers.FPix.HitProducer = 'siPixelRecHitsPreSplitting' + process.initialStepSeedLayers.BPix.HitProducer = 'siPixelRecHitsPreSplitting' + process.initialStepHitQuadruplets.SeedComparitorPSet.clusterShapeCacheSrc = "siPixelClusterShapeCachePreSplitting" + process.initialStepSeeds.SeedComparitorPSet.ClusterShapeCacheSrc = "siPixelClusterShapeCachePreSplitting" + if hasattr(process.initialStepTrackCandidates, "measurementTrackerEvent"): + # MkFit case + process.initialStepTrackCandidates.measurementTrackerEvent = 'MeasurementTrackerEventPreSplitting' + process.initialStepTrackCandidatesMkFitInput.pixelRecHits = "siPixelRecHitsPreSplitting" + else: + process.initialStepTrackCandidates.MeasurementTrackerEvent = 'MeasurementTrackerEventPreSplitting' + process.initialStepTracks.MeasurementTrackerEvent = 'MeasurementTrackerEventPreSplitting' + process.iterTrackingTask = cms.Task(process.trackerClusterCheck, + process.InitialStepTask) + + # Customize MTV + for selector in [ + process.cutsRecoTracksInitialStep, + process.cutsRecoTracksPt09InitialStep, + process.cutsRecoTracksFromPVInitialStep, + process.cutsRecoTracksFromPVPt09InitialStep, + ]: + selector.algorithm = [] + selector.quality = [] + selector.src = "initialStepTracks" + selector.vertexTag = "firstStepPrimaryVertices" + + process.trackingParticleRecoTrackAsssociationPreSplitting = process.trackingParticleRecoTrackAsssociation.clone( + label_tr = "initialStepTracks", + associator = "quickTrackAssociatorByHitsPreSplitting", + ) + process.VertexAssociatorByPositionAndTracksPreSplitting = process.VertexAssociatorByPositionAndTracks.clone( + trackAssociation = "trackingParticleRecoTrackAsssociationPreSplitting" + ) + + + def setInput(mtvs, labels): + for mtv in mtvs: + mod = getattr(process, mtv) + mod.label = labels + mod.label_vertex = "firstStepPrimaryVertices" + if mod.UseAssociators.value(): + mod.associators = ["quickTrackAssociatorByHitsPreSplitting"] + else: + mod.associators = ["trackingParticleRecoTrackAsssociationPreSplitting"] + mod.vertexAssociator = "VertexAssociatorByPositionAndTracksPreSplitting" + mod.trackCollectionForDrCalculation = "initialStepTracks" + mod.dodEdxPlots = False + mod.doResolutionPlotsForLabels = [] + + setInput(["trackValidatorTrackingOnly", "trackValidatorAllTPEfficStandalone", "trackValidatorTPPtLess09Standalone", "trackValidatorBHadronTrackingOnly"], + ["cutsRecoTracksInitialStep", "cutsRecoTracksPt09InitialStep"]) + setInput(["trackValidatorFromPVStandalone", "trackValidatorFromPVAllTPStandalone"], ["cutsRecoTracksFromPVInitialStep", "cutsRecoTracksFromPVPt09InitialStep"]) + setInput(["trackValidatorSeedingTrackingOnly"], ["seedTracksinitialStepSeeds"]) + setInput(["trackValidatorBuilding"], ["initialStepTracks"]) + process.trackValidatorBuilding.mvaLabels = cms.untracked.PSet(initialStepTracks = cms.untracked.vstring('initialStep')) + + process.tracksPreValidationTrackingOnly = cms.Task( + process.cutsRecoTracksInitialStep, + process.cutsRecoTracksPt09InitialStep, + process.cutsRecoTracksFromPVInitialStep, + process.cutsRecoTracksFromPVPt09InitialStep, + process.tracksValidationTruth, + process.trackingParticlesSignal, + process.trackingParticlesBHadron, + process.trackingParticleRecoTrackAsssociationPreSplitting, + process.VertexAssociatorByPositionAndTracksPreSplitting, + ) + process.trackValidatorsTrackingOnly.remove(process.trackValidatorConversionTrackingOnly) + process.trackValidatorsTrackingOnly.remove(process.trackValidatorSeedingPreSplittingTrackingOnly) + process.trackValidatorsTrackingOnly.remove(process.trackValidatorBuildingPreSplitting) + + # Remove vertex validation + process.globalPrevalidationTrackingOnly.remove(process.vertexValidationTrackingOnly) + + # Remove DQM + removePath(process, "dqmoffline_step") + removePath(process, "dqmofflineOnPAT_step") + + # Remove RECO output if it exists + removePath(process, "RECOSIMoutput_step") + + return process + + +def customizeInitialStepOnlyNoMTV(process): + process = customizeInitialStepOnly(process) + + process.options.wantSummary = cms.untracked.bool(True) + + # Remove validation + removePath(process, "prevalidation_step") + removePath(process, "validation_step") + removePath(process, "DQMoutput_step") + + # Minimize the rest + process.RawToDigiTask = cms.Task( + process.siPixelDigis, + process.siStripDigis + ) + process.reconstruction_trackingOnly = cms.Sequence( + process.trackerlocalreco + + process.offlineBeamSpot + + process.siPixelClusterShapeCachePreSplitting + + process.MeasurementTrackerEventPreSplitting, + process.iterTrackingTask + ) + process.trackerlocalrecoTask.remove(process.clusterSummaryProducer) + process.iterTrackingTask.remove(process.ak4CaloJetsForTrk) + process.iterTrackingTask.remove(process.caloTowerForTrk) + process.iterTrackingTask.remove(process.firstStepPrimaryVertices) + process.iterTrackingTask.remove(process.firstStepPrimaryVerticesUnsorted) + process.iterTrackingTask.remove(process.initialStepTrackRefsForJets) + process.iterTrackingTask.remove(process.initialStepClassifier1) + process.iterTrackingTask.remove(process.initialStep) + + # Add a dummy output module to trigger the (minimal) prefetching + process.out = cms.OutputModule("AsciiOutputModule", + outputCommands = cms.untracked.vstring( + "keep *_initialStepTracks_*_*", + ), + verbosity = cms.untracked.uint32(0) + ) + process.outPath = cms.EndPath(process.out) + process.schedule = cms.Schedule(process.raw2digi_step, process.reconstruction_step, process.outPath) + + # Minimize printouts + process.MessageLogger.cerr.FwkReport.reportEvery = 100 + process.MessageLogger.cerr.default.limit = 1 + process.MessageLogger.suppressWarning.extend([ + "initialStepTrackCandidatesMkFitInput", + "initialStepTrackCandidatesMkFit", + "initialStepTrackCandidates", + "initialStepTracks", + ]) + + return process diff --git a/RecoTracker/MkFit/python/customizeInitialStepToMkFit.py b/RecoTracker/MkFit/python/customizeInitialStepToMkFit.py new file mode 100644 index 0000000000000..84eb01fef3700 --- /dev/null +++ b/RecoTracker/MkFit/python/customizeInitialStepToMkFit.py @@ -0,0 +1,19 @@ +import RecoTracker.MkFit.mkFitInputConverter_cfi as mkFitInputConverter_cfi +import RecoTracker.MkFit.mkFitProducer_cfi as mkFitProducer_cfi +import RecoTracker.MkFit.mkFitOutputConverter_cfi as mkFitOutputConverter_cfi + +def customizeInitialStepToMkFit(process): + process.initialStepTrackCandidatesMkFitInput = mkFitInputConverter_cfi.mkFitInputConverter.clone( + seeds = "initialStepSeeds", + ) + process.initialStepTrackCandidatesMkFit = mkFitProducer_cfi.mkFitProducer.clone( + hitsSeeds = "initialStepTrackCandidatesMkFitInput", + ) + process.initialStepTrackCandidates = mkFitOutputConverter_cfi.mkFitOutputConverter.clone( + seeds = "initialStepSeeds", + hitsSeeds = "initialStepTrackCandidatesMkFitInput", + tracks = "initialStepTrackCandidatesMkFit", + ) + process.InitialStepTask.add(process.initialStepTrackCandidatesMkFitInput, + process.initialStepTrackCandidatesMkFit) + return process diff --git a/RecoTracker/MkFit/src/MkFitIndexLayer.cc b/RecoTracker/MkFit/src/MkFitIndexLayer.cc new file mode 100644 index 0000000000000..a9b8e091297f4 --- /dev/null +++ b/RecoTracker/MkFit/src/MkFitIndexLayer.cc @@ -0,0 +1,45 @@ +#include "RecoTracker/MkFit/interface/MkFitIndexLayer.h" + +#include "FWCore/Utilities/interface/Exception.h" + +#include + +void MkFitIndexLayer::insert(edm::ProductID id, size_t clusterIndex, int hit, int layer, const TrackingRecHit* hitPtr) { + // mapping CMSSW->mkfit + auto found = std::find_if(colls_.begin(), colls_.end(), [&](const auto& item) { return item.productID == id; }); + if (found == colls_.end()) { + found = colls_.emplace(colls_.end(), id); + } + if (found->infos.size() <= clusterIndex) { + found->infos.resize(clusterIndex + 1); + } + found->infos[clusterIndex] = HitInfo(hit, layer); + + // mapping mkfit->CMSSW + if (layer >= static_cast(hits_.size())) { + hits_.resize(layer + 1); + } + if (hit >= static_cast(hits_[layer].size())) { + hits_[layer].resize(hit + 1); + } + hits_[layer][hit].ptr = hitPtr; + hits_[layer][hit].clusterIndex = clusterIndex; +} + +const MkFitIndexLayer::HitInfo& MkFitIndexLayer::get(edm::ProductID id, size_t clusterIndex) const { + auto found = std::find_if(colls_.begin(), colls_.end(), [&](const auto& item) { return item.productID == id; }); + if (found == colls_.end()) { + auto exp = cms::Exception("Assert"); + exp << "Encountered a seed with a hit having productID " << id + << " which is not any of the input hit collections: "; + for (const auto& elem : colls_) { + exp << elem.productID << " "; + } + throw exp; + } + const HitInfo& ret = found->infos[clusterIndex]; + if (ret.index < 0) { + throw cms::Exception("Assert") << "No hit index for cluster " << clusterIndex << " of collection " << id; + } + return ret; +} diff --git a/RecoTracker/MkFit/src/MkFitInputWrapper.cc b/RecoTracker/MkFit/src/MkFitInputWrapper.cc new file mode 100644 index 0000000000000..a4f6ccc211865 --- /dev/null +++ b/RecoTracker/MkFit/src/MkFitInputWrapper.cc @@ -0,0 +1,20 @@ +#include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" + +// MkFit includes +#include "Hit.h" +#include "LayerNumberConverter.h" +#include "Track.h" + +MkFitInputWrapper::MkFitInputWrapper() = default; +MkFitInputWrapper::MkFitInputWrapper(MkFitIndexLayer&& indexLayers, + std::vector&& hits, + mkfit::TrackVec&& seeds, + mkfit::LayerNumberConverter&& lnc) + : indexLayers_{std::move(indexLayers)}, + hits_{std::move(hits)}, + seeds_{std::make_unique(std::move(seeds))}, + lnc_{std::make_unique(std::move(lnc))} {} + +MkFitInputWrapper::~MkFitInputWrapper() = default; + +unsigned int MkFitInputWrapper::nlayers() const { return lnc_->nLayers(); } diff --git a/RecoTracker/MkFit/src/MkFitOutputWrapper.cc b/RecoTracker/MkFit/src/MkFitOutputWrapper.cc new file mode 100644 index 0000000000000..6dadcd5c61337 --- /dev/null +++ b/RecoTracker/MkFit/src/MkFitOutputWrapper.cc @@ -0,0 +1,11 @@ +#include "RecoTracker/MkFit/interface/MkFitOutputWrapper.h" + +// MkFit includes +#include "Track.h" + +MkFitOutputWrapper::MkFitOutputWrapper() = default; + +MkFitOutputWrapper::MkFitOutputWrapper(mkfit::TrackVec&& candidateTracks, mkfit::TrackVec&& fitTracks) + : candidateTracks_{std::move(candidateTracks)}, fitTracks_{std::move(fitTracks)} {} + +MkFitOutputWrapper::~MkFitOutputWrapper() = default; diff --git a/RecoTracker/MkFit/src/classes.h b/RecoTracker/MkFit/src/classes.h new file mode 100644 index 0000000000000..c092179fb499f --- /dev/null +++ b/RecoTracker/MkFit/src/classes.h @@ -0,0 +1,8 @@ +#ifndef RecoTracker_MkFit_classes_h +#define RecoTracker_MkFit_classes_h + +#include "DataFormats/Common/interface/Wrapper.h" +#include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" +#include "RecoTracker/MkFit/interface/MkFitOutputWrapper.h" + +#endif diff --git a/RecoTracker/MkFit/src/classes_def.xml b/RecoTracker/MkFit/src/classes_def.xml new file mode 100644 index 0000000000000..e6363f8945f61 --- /dev/null +++ b/RecoTracker/MkFit/src/classes_def.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/Validation/RecoTrack/python/plotting/trackingPlots.py b/Validation/RecoTrack/python/plotting/trackingPlots.py index af9e8331e9bad..e0b184dedd214 100644 --- a/Validation/RecoTrack/python/plotting/trackingPlots.py +++ b/Validation/RecoTrack/python/plotting/trackingPlots.py @@ -1465,6 +1465,9 @@ def modules(self): "initialStepClassifier3", "initialStep", "initialStepSelector"], + building=["initialStepTrackCandidatesMkFitInput", + "initialStepTrackCandidatesMkFit", + "initialStepTrackCandidates"], other=["firstStepPrimaryVerticesUnsorted", "initialStepTrackRefsForJets", "caloTowerForTrk", From b557c23f450a46b37124407050f17d727ebf2006 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 14 Aug 2019 19:06:30 +0200 Subject: [PATCH 02/20] Add a matrix workflow testing MkFit in initialStep --- Configuration/PyReleaseValidation/python/relval_2017.py | 2 ++ Configuration/PyReleaseValidation/python/relval_steps.py | 7 +++++++ .../PyReleaseValidation/python/relval_upgrade.py | 4 +++- .../python/upgradeWorkflowComponents.py | 9 +++++++++ 4 files changed, 21 insertions(+), 1 deletion(-) diff --git a/Configuration/PyReleaseValidation/python/relval_2017.py b/Configuration/PyReleaseValidation/python/relval_2017.py index 55f33022e2e40..68bb0323078a1 100644 --- a/Configuration/PyReleaseValidation/python/relval_2017.py +++ b/Configuration/PyReleaseValidation/python/relval_2017.py @@ -23,6 +23,7 @@ # he collapse: TTbar, TTbar PU, TTbar design # ParkingBPH: TTbar # 2021 (ZMM, TTbar, ZEE, MinBias, TTbar PU, TTbar PU premix, ZEE PU, TTbar design) +# (TTbar trackingMkFit) # 2023 (TTbar, TTbar PU, TTbar PU premix) # 2024 (TTbar, TTbar PU, TTbar PU premix) numWFIB = [10001.0,10002.0,10003.0,10004.0,10005.0,10006.0,10007.0,10008.0,10009.0,10059.0,10071.0, @@ -34,6 +35,7 @@ 10824.6,11024.6,11224.6, 10824.8, 11650.0,11634.0,11646.0,11640.0,11834.0,11834.99,11846.0,12024.0, + 11634.7, 12434.0,12634.0,12634.99, 12834.0,13034.0,13034.99] for numWF in numWFIB: diff --git a/Configuration/PyReleaseValidation/python/relval_steps.py b/Configuration/PyReleaseValidation/python/relval_steps.py index f619829abf873..fe78cc98c8726 100644 --- a/Configuration/PyReleaseValidation/python/relval_steps.py +++ b/Configuration/PyReleaseValidation/python/relval_steps.py @@ -2090,6 +2090,9 @@ def gen2018HiMix(fragment,howMuch): step3_trackingLowPU = { '--era': 'Run2_2016_trackingLowPU' } +step3_trackingMkFit = { + '--customise': 'RecoTracker/MkFit/customizeInitialStepToMkFit.customizeInitialStepToMkFit' +} step3_HIPM = { '--era': 'Run2_2016_HIPM' } @@ -3240,6 +3243,10 @@ def gen2018HiMix(fragment,howMuch): if 'Reco' in step and upgradeStepDict[step][k]['--era']=='Run2_2017': upgradeStepDict[stepName][k] = merge([{'--era': 'Run2_2017_trackingLowPU'}, upgradeStepDict[step][k]]) + for step in upgradeSteps['trackingMkFit']['steps']: + stepName = step + upgradeSteps['trackingMkFit']['suffix'] + if 'Reco' in step: upgradeStepDict[stepName][k] = merge([step3_trackingMkFit, upgradeStepDict[step][k]]) + for step in upgradeSteps['Neutron']['steps']: if 'GenSim' in step: custNew = "SimG4Core/Application/NeutronBGforMuonsXS_cff.customise" diff --git a/Configuration/PyReleaseValidation/python/relval_upgrade.py b/Configuration/PyReleaseValidation/python/relval_upgrade.py index c48247d4e08e8..3c22488676614 100644 --- a/Configuration/PyReleaseValidation/python/relval_upgrade.py +++ b/Configuration/PyReleaseValidation/python/relval_upgrade.py @@ -73,7 +73,7 @@ def makeStepName(key,frag,step,suffix): # special workflows for tracker if (upgradeDatasetFromFragment[frag]=="TTbar_13" or upgradeDatasetFromFragment[frag]=="TTbar_14TeV") and not 'PU' in key and hasHarvest: # skip ALCA and Nano - trackingVariations = ['trackingOnly','trackingRun2','trackingOnlyRun2','trackingLowPU','pixelTrackingOnly'] + trackingVariations = ['trackingOnly','trackingRun2','trackingOnlyRun2','trackingLowPU','pixelTrackingOnly','trackingMkFit'] for tv in trackingVariations: stepList[tv] = [s for s in stepList[tv] if (("ALCA" not in s) and ("Nano" not in s))] workflows[numWF+upgradeSteps['trackingOnly']['offset']] = [ upgradeDatasetFromFragment[frag], stepList['trackingOnly']] @@ -82,6 +82,8 @@ def makeStepName(key,frag,step,suffix): workflows[numWF+upgradeSteps[tv]['offset']] = [ upgradeDatasetFromFragment[frag], stepList[tv]] elif '2018' in key: workflows[numWF+upgradeSteps['pixelTrackingOnly']['offset']] = [ upgradeDatasetFromFragment[frag], stepList['pixelTrackingOnly']] + elif '2021' in key: + workflows[numWF+upgradeSteps['trackingMkFit']['offset']] = [ upgradeDatasetFromFragment[frag], stepList['trackingMkFit']] # special workflows for HGCAL/TICL if (upgradeDatasetFromFragment[frag]=="CloseByParticleGun") and ('2026' in key): diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index 30179cd68e939..93896bb950664 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -188,6 +188,15 @@ 'suffix' : '_heCollapse', 'offset' : 0.6, } +upgradeSteps['trackingMkFit'] = { + 'steps' : [ + 'RecoFull', + 'RecoFullGlobal', + ], + 'PU' : [], + 'suffix' : '_trackingMkFit', + 'offset' : 0.7, +} upgradeSteps['ParkingBPH'] = { 'steps' : [ 'RecoFull', From a10f90ecaf81fed2cdd62a259c06883e9f3c31cd Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Fri, 30 Aug 2019 20:54:54 +0200 Subject: [PATCH 03/20] Address review comments - Replace dynamic_cast on TrackerRecHit with trackerHitRTTI::isFromDet() + static_cast - Make the offset business more clear in the layer structure building - Reduce temporaries - Restore exceptions The seed-candidate hit index mismatch seems not to be occurring anymore (at least on a test of 100k 10mu and 10k ttbar nopu+PU50+PU70 events). --- RecoTracker/MkFit/README.md | 3 ++ .../MkFit/plugins/MkFitInputConverter.cc | 10 ++-- .../MkFit/plugins/MkFitOutputConverter.cc | 48 ++++++++++--------- 3 files changed, 33 insertions(+), 28 deletions(-) diff --git a/RecoTracker/MkFit/README.md b/RecoTracker/MkFit/README.md index 7bbfde5274be7..11bc66edf6800 100644 --- a/RecoTracker/MkFit/README.md +++ b/RecoTracker/MkFit/README.md @@ -6,6 +6,9 @@ This package holds the glue modules for running Note that at the moment there may be only one `MkFitProducer` in a single job. This restriction will be removed in the future. +Also note that at the moment the mkFit works only with the CMS phase1 +tracker detector. Support for the phase2 tracker will be added later. + ## Customize functions for runTheMatrix workflows (offline reconstruction) * `RecoTracker/MkFit/customizeInitialStepToMkFit.customizeInitialStepToMkFit` diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc index d6e70527bc6ee..f740574988d6b 100644 --- a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -8,6 +8,7 @@ #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h" #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h" +#include "DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" @@ -195,12 +196,11 @@ mkfit::TrackVec MkFitInputConverter::convertSeeds(const edm::View(&*iHit); - if (hit == nullptr) { - throw cms::Exception("Assert") << "Encountered a seed with a hit which is not BaseTrackerRecHit"; + if (not trackerHitRTTI::isFromDet(*iHit)) { + throw cms::Exception("Assert") << "Encountered a seed with a hit which is not trackerHitRTTI::isFromDet()"; } - - const auto& info = indexLayers.get(hit->firstClusterRef().id(), hit->firstClusterRef().index()); + const auto& clusterRef = static_cast(*iHit).firstClusterRef(); + const auto& info = indexLayers.get(clusterRef.id(), clusterRef.index()); ret.back().addHitIdx(info.index, info.layer, 0); // per-hit chi2 is not known } ++index; diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 4df7d011cd0d4..b3d8bbefa0893 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -204,12 +204,17 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: // TODO: currently hardcoded... // Logic copied from mkfit::LayerNumberConverter - unsigned int off = 1; + unsigned int off = 0; // BPix - set(0, ttopo.pxbDetId(1, 0, 0)); - set(1, ttopo.pxbDetId(2, 0, 0)); - set(2, ttopo.pxbDetId(3, 0, 0)); - set(3, ttopo.pxbDetId(4, 0, 0)); + set(off + 0, ttopo.pxbDetId(1, 0, 0)); + set(off + 1, ttopo.pxbDetId(2, 0, 0)); + set(off + 2, ttopo.pxbDetId(3, 0, 0)); + set(off + 3, ttopo.pxbDetId(4, 0, 0)); + // offset needs to be increased by 1 here to accommodate the 4th + // pixel barrel layer, while keeping the off+N numbering consistent + // with mkfit::LayerNumberConverter (that, to limited degree, + // supports the layer numbering for phase0 pixel as well) + off += 1; // TIB set(off + 3, ttopo.tibDetId(1, 0, 0, 0, 0, 0)); set(off + 4, ttopo.tibDetId(1, 0, 0, 0, 0, 1)); @@ -227,13 +232,13 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: set(off + 15, ttopo.tobDetId(5, 0, 0, 0, 0)); set(off + 16, ttopo.tobDetId(6, 0, 0, 0, 0)); - auto setForward = [&](unsigned int side) { + auto setForward = [&set, &ttopo](unsigned int off, unsigned int side) { // FPix set(off + 0, ttopo.pxfDetId(side, 1, 0, 0, 0)); set(off + 1, ttopo.pxfDetId(side, 2, 0, 0, 0)); set(off + 2, ttopo.pxfDetId(side, 3, 0, 0, 0)); // TID+ - off += 1; + off += 1; // see comment above for barrel set(off + 2, ttopo.tidDetId(side, 1, 0, 0, 0, 0)); set(off + 3, ttopo.tidDetId(side, 1, 0, 0, 0, 1)); set(off + 4, ttopo.tidDetId(side, 2, 0, 0, 0, 0)); @@ -261,13 +266,16 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: set(off + 25, ttopo.tecDetId(side, 9, 0, 0, 0, 0, 1)); }; + constexpr unsigned int nlay_barrel = 16 + 1; // 16 in phase0, 1 more in phase1 + constexpr unsigned int nlay_endcap = 25 + 1; // 25 in phase0, 1 more in phase1 + // plus - off = 17 + 1; - setForward(2); + off = nlay_barrel + 1; // +1 to move to next slot + setForward(off, 2); // minus - off = 17 + 1 + 25 + 2; - setForward(1); + off += nlay_endcap + 1; // +1 to move to next slot + setForward(off, 1); return dets; } @@ -335,11 +343,8 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp // MkFit hits are *not* in the order of propagation, sort by 3D radius for now (as we don't have loopers) // TODO: Improve the sorting (extract keys? maybe even bubble sort would work well as the hits are almost in the correct order) recHits.sort([](const auto& a, const auto& b) { - const auto aid = a.geographicalId(); - const auto bid = b.geographicalId(); - - const auto asub = aid.subdetId(); - const auto bsub = bid.subdetId(); + const auto asub = a.geographicalId().subdetId(); + const auto bsub = b.geographicalId().subdetId(); if (asub != bsub) { // Subdetector order (BPix, FPix, TIB, TID, TOB, TEC) corresponds also the navigation return asub < bsub; @@ -366,19 +371,16 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp // sanity check for now const auto& candHitOnTrack = cand.getHitOnTrack(i); if (hitOnTrack.layer != candHitOnTrack.layer) { - //throw cms::Exception("LogicError") - edm::LogError("MkFitOutputConverter") + 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") - edm::LogError("MkFitOutputConverter") - << "Candidate " << candIndex << " from seed " << seedIndex << " hit " << i - << " has different hit index in candidate (" << candHitOnTrack.index << ") and seed (" << hitOnTrack.index - << ") on layer " << hitOnTrack.layer; + 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; } } From 6f032c408925126b6ba4276f948d9eaacb3b1cdf Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 17:02:17 +0200 Subject: [PATCH 04/20] Move move constructors and assignment operators out-of-line The other constructors already were out-of-line, and this way the move operations can be called everywhere without requiring the #including the mkFit headers. --- RecoTracker/MkFit/interface/MkFitInputWrapper.h | 4 ++-- RecoTracker/MkFit/interface/MkFitOutputWrapper.h | 4 ++-- RecoTracker/MkFit/src/MkFitInputWrapper.cc | 4 ++++ RecoTracker/MkFit/src/MkFitOutputWrapper.cc | 3 +++ 4 files changed, 11 insertions(+), 4 deletions(-) diff --git a/RecoTracker/MkFit/interface/MkFitInputWrapper.h b/RecoTracker/MkFit/interface/MkFitInputWrapper.h index a281b66c3e415..2645716b9efda 100644 --- a/RecoTracker/MkFit/interface/MkFitInputWrapper.h +++ b/RecoTracker/MkFit/interface/MkFitInputWrapper.h @@ -25,8 +25,8 @@ class MkFitInputWrapper { MkFitInputWrapper(MkFitInputWrapper const&) = delete; MkFitInputWrapper& operator=(MkFitInputWrapper const&) = delete; - MkFitInputWrapper(MkFitInputWrapper&&) = default; - MkFitInputWrapper& operator=(MkFitInputWrapper&&) = default; + MkFitInputWrapper(MkFitInputWrapper&&); + MkFitInputWrapper& operator=(MkFitInputWrapper&&); MkFitIndexLayer const& indexLayers() const { return indexLayers_; } mkfit::TrackVec const& seeds() const { return *seeds_; } diff --git a/RecoTracker/MkFit/interface/MkFitOutputWrapper.h b/RecoTracker/MkFit/interface/MkFitOutputWrapper.h index d03a4af250556..8297b7c511f77 100644 --- a/RecoTracker/MkFit/interface/MkFitOutputWrapper.h +++ b/RecoTracker/MkFit/interface/MkFitOutputWrapper.h @@ -16,8 +16,8 @@ class MkFitOutputWrapper { MkFitOutputWrapper(MkFitOutputWrapper const&) = delete; MkFitOutputWrapper& operator=(MkFitOutputWrapper const&) = delete; - MkFitOutputWrapper(MkFitOutputWrapper&&) = default; - MkFitOutputWrapper& operator=(MkFitOutputWrapper&&) = default; + MkFitOutputWrapper(MkFitOutputWrapper&&); + MkFitOutputWrapper& operator=(MkFitOutputWrapper&&); mkfit::TrackVec const& candidateTracks() const { return candidateTracks_; } mkfit::TrackVec const& fitTracks() const { return fitTracks_; } diff --git a/RecoTracker/MkFit/src/MkFitInputWrapper.cc b/RecoTracker/MkFit/src/MkFitInputWrapper.cc index a4f6ccc211865..f5274009fa179 100644 --- a/RecoTracker/MkFit/src/MkFitInputWrapper.cc +++ b/RecoTracker/MkFit/src/MkFitInputWrapper.cc @@ -6,6 +6,7 @@ #include "Track.h" MkFitInputWrapper::MkFitInputWrapper() = default; + MkFitInputWrapper::MkFitInputWrapper(MkFitIndexLayer&& indexLayers, std::vector&& hits, mkfit::TrackVec&& seeds, @@ -17,4 +18,7 @@ MkFitInputWrapper::MkFitInputWrapper(MkFitIndexLayer&& indexLayers, MkFitInputWrapper::~MkFitInputWrapper() = default; +MkFitInputWrapper::MkFitInputWrapper(MkFitInputWrapper&&) = default; +MkFitInputWrapper& MkFitInputWrapper::operator=(MkFitInputWrapper&&) = default; + unsigned int MkFitInputWrapper::nlayers() const { return lnc_->nLayers(); } diff --git a/RecoTracker/MkFit/src/MkFitOutputWrapper.cc b/RecoTracker/MkFit/src/MkFitOutputWrapper.cc index 6dadcd5c61337..9088ad23644e0 100644 --- a/RecoTracker/MkFit/src/MkFitOutputWrapper.cc +++ b/RecoTracker/MkFit/src/MkFitOutputWrapper.cc @@ -9,3 +9,6 @@ MkFitOutputWrapper::MkFitOutputWrapper(mkfit::TrackVec&& candidateTracks, mkfit: : candidateTracks_{std::move(candidateTracks)}, fitTracks_{std::move(fitTracks)} {} MkFitOutputWrapper::~MkFitOutputWrapper() = default; + +MkFitOutputWrapper::MkFitOutputWrapper(MkFitOutputWrapper&&) = default; +MkFitOutputWrapper& MkFitOutputWrapper::operator=(MkFitOutputWrapper&&) = default; From 9b81d92b285b92677a7fadeb9fa3e6406a11fc89 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 17:07:55 +0200 Subject: [PATCH 05/20] Sort BuildFile --- RecoTracker/MkFit/plugins/BuildFile.xml | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/RecoTracker/MkFit/plugins/BuildFile.xml b/RecoTracker/MkFit/plugins/BuildFile.xml index facea9b6356da..7574ba951ba0c 100644 --- a/RecoTracker/MkFit/plugins/BuildFile.xml +++ b/RecoTracker/MkFit/plugins/BuildFile.xml @@ -1,27 +1,27 @@ - - - - - - + + + + + + + + + - - - - + From 169fdfad3767e56793b664b1cdfbfe433ef4aa97 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 17:15:34 +0200 Subject: [PATCH 06/20] More consistent capitalization of mkFit --- RecoTracker/MkFit/interface/MkFitIndexLayer.h | 4 ++-- .../MkFit/plugins/MkFitInputConverter.cc | 22 +++++++++---------- .../MkFit/plugins/MkFitOutputConverter.cc | 18 +++++++-------- RecoTracker/MkFit/plugins/MkFitProducer.cc | 14 ++++++------ .../MkFit/python/customizeInitialStepOnly.py | 2 +- RecoTracker/MkFit/src/MkFitInputWrapper.cc | 2 +- RecoTracker/MkFit/src/MkFitOutputWrapper.cc | 2 +- 7 files changed, 32 insertions(+), 32 deletions(-) diff --git a/RecoTracker/MkFit/interface/MkFitIndexLayer.h b/RecoTracker/MkFit/interface/MkFitIndexLayer.h index 44035ecb17918..6482809b96b4c 100644 --- a/RecoTracker/MkFit/interface/MkFitIndexLayer.h +++ b/RecoTracker/MkFit/interface/MkFitIndexLayer.h @@ -38,8 +38,8 @@ class MkFitIndexLayer { size_t clusterIndex = 0; }; - std::vector colls_; // mapping from CMSSW(ProductID, index) -> mkfit(index, layer) - std::vector > hits_; // reverse mapping, mkfit(layer, index) -> CMSSW hit + std::vector colls_; // mapping from CMSSW(ProductID, index) -> mkFit(index, layer) + std::vector > hits_; // reverse mapping, mkFit(layer, index) -> CMSSW hit }; #endif diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc index f740574988d6b..740c0240c9272 100644 --- a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -26,7 +26,7 @@ #include "Math/SVector.h" #include "Math/SMatrix.h" -// MkFit includes +// mkFit includes #include "Hit.h" #include "Track.h" #include "LayerNumberConverter.h" @@ -43,7 +43,7 @@ class MkFitInputConverter : public edm::global::EDProducer<> { template void convertHits(const HitCollection& hits, - std::vector& mkfitHits, + std::vector& mkFitHits, MkFitIndexLayer& indexLayers, int& totalHits, const TrackerTopology& ttopo, @@ -104,17 +104,17 @@ void MkFitInputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const e const auto& ttrhBuilder = iSetup.getData(ttrhBuilderToken_); const auto& ttopo = iSetup.getData(ttopoToken_); - std::vector mkfitHits(lnc.nLayers()); + std::vector 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); + 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_)); + 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)); + 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 { @@ -125,7 +125,7 @@ bool MkFitInputConverter::passCCC(const SiPixelRecHit& hit, const DetId hitId) c template void MkFitInputConverter::convertHits(const HitCollection& hits, - std::vector& mkfitHits, + std::vector& mkFitHits, MkFitIndexLayer& indexLayers, int& totalHits, const TrackerTopology& ttopo, @@ -157,8 +157,8 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, << 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); + indexLayers.insert(hit.firstClusterRef().id(), hit.firstClusterRef().index(), mkFitHits[ilay].size(), ilay, &hit); + mkFitHits[ilay].emplace_back(pos, err, totalHits); ++totalHits; } } diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index b3d8bbefa0893..277b08c74abab 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -39,7 +39,7 @@ #include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" #include "RecoTracker/MkFit/interface/MkFitOutputWrapper.h" -// MkFit indludes +// mkFit indludes #include "LayerNumberConverter.h" #include "Track.h" @@ -71,7 +71,7 @@ class MkFitOutputConverter : public edm::global::EDProducer<> { const GeometricSearchTracker& tracker, const TrackerTopology& ttopo) const; - TrackCandidateCollection convertCandidates(const MkFitOutputWrapper& mkfitOutput, + TrackCandidateCollection convertCandidates(const MkFitOutputWrapper& mkFitOutput, const MkFitIndexLayer& indexLayers, const edm::View& seeds, const TrackerGeometry& geom, @@ -80,7 +80,7 @@ class MkFitOutputConverter : public edm::global::EDProducer<> { const Propagator& propagatorOpposite, const TkClonerImpl& hitCloner, const std::vector& detLayers, - const mkfit::TrackVec& mkfitSeeds) const; + const mkfit::TrackVec& mkFitSeeds) const; std::pair backwardFit(const FreeTrajectoryState& fts, const edm::OwnVector& hits, @@ -144,7 +144,7 @@ void MkFitOutputConverter::fillDescriptions(edm::ConfigurationDescriptions& desc edm::ParameterSetDescription desc; desc.add("hitsSeeds", edm::InputTag{"mkFitInputConverter"}); - desc.add("tracks", edm::InputTag{"mkfitProducer"}); + desc.add("tracks", edm::InputTag{"mkFitProducer"}); desc.add("seeds", edm::InputTag{"initialStepSeeds"}); desc.add("measurementTrackerEvent", edm::InputTag{"MeasurementTrackerEvent"}); desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"}); @@ -280,7 +280,7 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: return dets; } -TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutputWrapper& mkfitOutput, +TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutputWrapper& mkFitOutput, const MkFitIndexLayer& indexLayers, const edm::View& seeds, const TrackerGeometry& geom, @@ -289,12 +289,12 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp const Propagator& propagatorOpposite, const TkClonerImpl& hitCloner, const std::vector& detLayers, - const mkfit::TrackVec& mkfitSeeds) const { + const mkfit::TrackVec& mkFitSeeds) const { TrackCandidateCollection output; - const auto& candidates = backwardFitInCMSSW_ ? mkfitOutput.candidateTracks() : mkfitOutput.fitTracks(); + const auto& candidates = backwardFitInCMSSW_ ? mkFitOutput.candidateTracks() : mkFitOutput.fitTracks(); output.reserve(candidates.size()); - LogTrace("MkFitOutputConverter") << "Number of candidates " << mkfitOutput.candidateTracks().size(); + LogTrace("MkFitOutputConverter") << "Number of candidates " << mkFitOutput.candidateTracks().size(); int candIndex = -1; for (const auto& cand : candidates) { @@ -364,7 +364,7 @@ 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()); + 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; diff --git a/RecoTracker/MkFit/plugins/MkFitProducer.cc b/RecoTracker/MkFit/plugins/MkFitProducer.cc index 06c40228da85e..af51bde5843e6 100644 --- a/RecoTracker/MkFit/plugins/MkFitProducer.cc +++ b/RecoTracker/MkFit/plugins/MkFitProducer.cc @@ -14,7 +14,7 @@ #include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" #include "RecoTracker/MkFit/interface/MkFitOutputWrapper.h" -// MkFit includes +// mkFit includes #include "ConfigWrapper.h" #include "Event.h" #include "mkFit/buildtestMPlex.h" @@ -43,7 +43,7 @@ class MkFitProducer : public edm::global::EDProducer putToken_; std::function buildFunction_; bool backwardFitInCMSSW_; - bool mkfitSilent_; + bool mkFitSilent_; }; MkFitProducer::MkFitProducer(edm::ParameterSet const& iConfig) @@ -51,7 +51,7 @@ MkFitProducer::MkFitProducer(edm::ParameterSet const& iConfig) geomToken_{esConsumes()}, putToken_{produces()}, backwardFitInCMSSW_{iConfig.getParameter("backwardFitInCMSSW")}, - mkfitSilent_{iConfig.getUntrackedParameter("mkfitSilent")} { + mkFitSilent_{iConfig.getUntrackedParameter("mkFitSilent")} { const auto build = iConfig.getParameter("buildingRoutine"); bool isFV = false; if (build == "bestHit") { @@ -84,7 +84,7 @@ MkFitProducer::MkFitProducer(edm::ParameterSet const& iConfig) // TODO: what to do when we have multiple instances of MkFitProducer in a job? mkfit::MkBuilderWrapper::populate(isFV); - mkfit::ConfigWrapper::initializeForCMSSW(seedCleanOpt, backwardFitOpt, mkfitSilent_); + mkfit::ConfigWrapper::initializeForCMSSW(seedCleanOpt, backwardFitOpt, mkFitSilent_); } void MkFitProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { @@ -96,7 +96,7 @@ void MkFitProducer::fillDescriptions(edm::ConfigurationDescriptions& description desc.add("seedCleaning", "N2")->setComment("Valid values are: 'none', 'N2'"); desc.add("backwardFitInCMSSW", false) ->setComment("Do backward fit (to innermost hit) in CMSSW (true) or mkFit (false)"); - desc.addUntracked("mkfitSilent", true)->setComment("Allows to enables printouts from mkfit with 'False'"); + desc.addUntracked("mkFitSilent", true)->setComment("Allows to enables printouts from mkFit with 'False'"); descriptions.add("mkFitProducer", desc); } @@ -116,7 +116,7 @@ void MkFitProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::Ev // TODO: the mechanism needs to be improved... std::call_once(geometryFlag, [&geom, nlayers = hitsSeeds.nlayers()]() { // TODO: eventually automatize fully - // For now it is easier to use purely the infrastructure from mkfit + // For now it is easier to use purely the infrastructure from mkFit /* const auto barrelLayers = geom.numberOfLayers(PixelSubdetector::PixelBarrel) + geom.numberOfLayers(StripSubdetector::TIB) + geom.numberOfLayers(StripSubdetector::TOB); const auto endcapLayers = geom.numberOfLayers(PixelSubdetector::PixelEndcap) + geom.numberOfLayers(StripSubdetector::TID) + geom.numberOfLayers(StripSubdetector::TEC); @@ -135,7 +135,7 @@ void MkFitProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::Ev // CMSSW event ID (64-bit unsigned) does not fit in int // In addition, unique ID requires also lumi and run - // But does the event ID really matter within mkfit? + // But does the event ID really matter within mkFit? mkfit::Event ev(iEvent.id().event()); ev.setInputFromCMSSW(hitsSeeds.hits(), hitsSeeds.seeds()); diff --git a/RecoTracker/MkFit/python/customizeInitialStepOnly.py b/RecoTracker/MkFit/python/customizeInitialStepOnly.py index 5a7fdff65e571..1b42f710c3e35 100644 --- a/RecoTracker/MkFit/python/customizeInitialStepOnly.py +++ b/RecoTracker/MkFit/python/customizeInitialStepOnly.py @@ -13,7 +13,7 @@ def customizeInitialStepOnly(process): process.initialStepHitQuadruplets.SeedComparitorPSet.clusterShapeCacheSrc = "siPixelClusterShapeCachePreSplitting" process.initialStepSeeds.SeedComparitorPSet.ClusterShapeCacheSrc = "siPixelClusterShapeCachePreSplitting" if hasattr(process.initialStepTrackCandidates, "measurementTrackerEvent"): - # MkFit case + # mkFit case process.initialStepTrackCandidates.measurementTrackerEvent = 'MeasurementTrackerEventPreSplitting' process.initialStepTrackCandidatesMkFitInput.pixelRecHits = "siPixelRecHitsPreSplitting" else: diff --git a/RecoTracker/MkFit/src/MkFitInputWrapper.cc b/RecoTracker/MkFit/src/MkFitInputWrapper.cc index f5274009fa179..54ee931bbc5ac 100644 --- a/RecoTracker/MkFit/src/MkFitInputWrapper.cc +++ b/RecoTracker/MkFit/src/MkFitInputWrapper.cc @@ -1,6 +1,6 @@ #include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" -// MkFit includes +// mkFit includes #include "Hit.h" #include "LayerNumberConverter.h" #include "Track.h" diff --git a/RecoTracker/MkFit/src/MkFitOutputWrapper.cc b/RecoTracker/MkFit/src/MkFitOutputWrapper.cc index 9088ad23644e0..0c670e22a7660 100644 --- a/RecoTracker/MkFit/src/MkFitOutputWrapper.cc +++ b/RecoTracker/MkFit/src/MkFitOutputWrapper.cc @@ -1,6 +1,6 @@ #include "RecoTracker/MkFit/interface/MkFitOutputWrapper.h" -// MkFit includes +// mkFit includes #include "Track.h" MkFitOutputWrapper::MkFitOutputWrapper() = default; From 5cca5e4cfb722c7cbb7b0c9dbbd922658f58bdbd Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 17:17:00 +0200 Subject: [PATCH 07/20] Improve comments Improve content, and also move comments for better formatting --- .../MkFit/plugins/MkFitOutputConverter.cc | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 277b08c74abab..d25aac51c9db3 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -304,13 +304,17 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp // hits edm::OwnVector recHits; - const int nhits = cand.nTotalHits(); // what exactly is the difference between nTotalHits() and nFoundHits()? + // nTotalHits() gives sum of valid hits (nFoundHits()) and + // invalid/missing hits (up to a maximum of 32 inside mkFit, + // restriction to be lifted in the future) + const int nhits = cand.nTotalHits(); bool lastHitInvalid = false; for (int i = 0; i < nhits; ++i) { const auto& hitOnTrack = cand.getHitOnTrack(i); LogTrace("MkFitOutputConverter") << " hit on layer " << hitOnTrack.layer << " index " << hitOnTrack.index; if (hitOnTrack.index < 0) { - // What is the exact meaning of -1, -2, -3? + // See index-desc.txt file in mkFit for description of negative values + // // In order to use the regular InvalidTrackingRecHit I'd need // a GeomDet (and "unfortunately" that is needed in // TrackProducer). @@ -425,7 +429,7 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp seeds.at(seedIndex), pstate, seeds.refAt(seedIndex), - 0, // TODO: nloops, let's ignore for now + 0, // mkFit does not produce loopers, so set nLoops=0 static_cast(StopReason::UNINITIALIZED) // TODO: ignore details of stopping reason as well for now ); } @@ -522,14 +526,14 @@ std::pair MkFitOutputConverter::backwa KFTrajectoryFitter backFitter( &propagatorAlong, &aKFUpdator, &aChi2MeasurementEstimator, firstHits.size(), nullptr, &hitCloner); - PropagationDirection backFitDirection = - oppositeToMomentum; // assume for now that the propagation in mkfit always alongMomentum + // assume for now that the propagation in mkfit always alongMomentum + PropagationDirection backFitDirection = oppositeToMomentum; - // only direction matters in this contest + // only direction matters in this context TrajectorySeed fakeSeed(PTrajectoryStateOnDet(), edm::OwnVector(), backFitDirection); - Trajectory fitres = - backFitter.fitOne(fakeSeed, firstHits, startingState, TrajectoryFitter::standard); // ignore loopers for now + // ignore loopers for now + Trajectory fitres = backFitter.fitOne(fakeSeed, firstHits, startingState, TrajectoryFitter::standard); LogDebug("MkFitOutputConverter") << "using a backward fit of :" << firstHits.size() << " hits, starting from:\n" << startingState << " to get the estimate of the initial state of the track."; From 83aedfe3b83d3ce292d8ec253a37353b11115014 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 18:24:38 +0200 Subject: [PATCH 08/20] Move number of layers check out of the call_once, remove commented code --- RecoTracker/MkFit/plugins/MkFitProducer.cc | 27 +++++++--------------- 1 file changed, 8 insertions(+), 19 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitProducer.cc b/RecoTracker/MkFit/plugins/MkFitProducer.cc index af51bde5843e6..33a2ac0f33853 100644 --- a/RecoTracker/MkFit/plugins/MkFitProducer.cc +++ b/RecoTracker/MkFit/plugins/MkFitProducer.cc @@ -112,26 +112,15 @@ void MkFitProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::Ev const auto& hitsSeeds = iEvent.get(hitsSeedsToken_); const auto& geom = iSetup.getData(geomToken_); - // First do some initialization + if (geom.numberOfLayers(PixelSubdetector::PixelBarrel) != 4 || + geom.numberOfLayers(PixelSubdetector::PixelEndcap) != 3) { + throw cms::Exception("Assert") << "For now this code works only with phase1 tracker, you have something else"; + } + + // 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, [&geom, nlayers = hitsSeeds.nlayers()]() { - // TODO: eventually automatize fully - // For now it is easier to use purely the infrastructure from mkFit - /* - const auto barrelLayers = geom.numberOfLayers(PixelSubdetector::PixelBarrel) + geom.numberOfLayers(StripSubdetector::TIB) + geom.numberOfLayers(StripSubdetector::TOB); - const auto endcapLayers = geom.numberOfLayers(PixelSubdetector::PixelEndcap) + geom.numberOfLayers(StripSubdetector::TID) + geom.numberOfLayers(StripSubdetector::TEC); - // TODO: Number of stereo layers is hardcoded for now, so this won't work for phase2 tracker - const auto barrelStereo = 2 + 2; - const auto endcapStereo = geom.numberOfLayers(StripSubdetector::TID) + geom.numberOfLayers(StripSubdetector::TEC); - const auto nlayers = barrelLayers + barrelStereo + 2*(endcapLayers + endcapStereo); - LogDebug("MkFitProducer") << "Total number of tracker layers (stereo counted separately) " << nlayers; - */ - if (geom.numberOfLayers(PixelSubdetector::PixelBarrel) != 4 || - geom.numberOfLayers(PixelSubdetector::PixelEndcap) != 3) { - throw cms::Exception("Assert") << "For now this code works only with phase1 tracker, you have something else"; - } - mkfit::ConfigWrapper::setNTotalLayers(nlayers); - }); + std::call_once(geometryFlag, [nlayers = hitsSeeds.nlayers()]() { mkfit::ConfigWrapper::setNTotalLayers(nlayers); }); // CMSSW event ID (64-bit unsigned) does not fit in int // In addition, unique ID requires also lumi and run From 148ce4b2c78d21b88210ff556f7c1a6daae9dfac Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 18:38:05 +0200 Subject: [PATCH 09/20] Merge MkFitOutputConverter::backwardFitImpl() with backwardFit() --- .../MkFit/plugins/MkFitOutputConverter.cc | 22 ------------------- 1 file changed, 22 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index d25aac51c9db3..84978929e8e24 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -90,15 +90,6 @@ class MkFitOutputConverter : public edm::global::EDProducer<> { bool lastHitWasInvalid, bool lastHitWasChanged) const; - std::pair backwardFitImpl( - const FreeTrajectoryState& fts, - const TransientTrackingRecHit::ConstRecHitContainer& firstHits, - const Propagator& propagatorAlong, - const Propagator& propagatorOpposite, - const TkClonerImpl& hitCloner, - bool lastHitWasInvalid, - bool lastHitWasChanged) const; - std::pair convertInnermostState(const FreeTrajectoryState& fts, const edm::OwnVector& hits, const Propagator& propagatorAlong, @@ -457,19 +448,6 @@ std::pair MkFitOutputConverter::backwa } } - auto ret = backwardFitImpl( - fts, firstHits, propagatorAlong, propagatorOpposite, hitCloner, lastHitWasInvalid, lastHitWasChanged); - return ret; -} - -std::pair MkFitOutputConverter::backwardFitImpl( - const FreeTrajectoryState& fts, - const TransientTrackingRecHit::ConstRecHitContainer& firstHits, - const Propagator& propagatorAlong, - const Propagator& propagatorOpposite, - const TkClonerImpl& hitCloner, - bool lastHitWasInvalid, - bool lastHitWasChanged) const { // Then propagate along to the surface of the last hit to get a TSOS const auto& lastHitSurface = firstHits.front()->det()->surface(); From ed7bfead3d0045fd92df3a7a4d377671aeb4bf33 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 18:39:20 +0200 Subject: [PATCH 10/20] Remove unnecessary commented code --- RecoTracker/MkFit/plugins/MkFitOutputConverter.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 84978929e8e24..1d38ed91cd8fc 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -318,8 +318,9 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp if (detLayer == nullptr) { throw cms::Exception("LogicError") << "DetLayer for layer index " << hitOnTrack.layer << " is null!"; } - // Actually it is necessary to leave dealing with invalid hits to the TrackProducer? - //recHits.push_back(new InvalidTrackingRecHitNoDet(detLayer->surface(), TrackingRecHit::missing)); // let's put them all as missing for now + // In principle an InvalidTrackingRecHitNoDet could be + // inserted here, but it seems that it is best to deal with + // them in the TrackProducer. lastHitInvalid = true; } else { recHits.push_back(indexLayers.getHitPtr(hitOnTrack.layer, hitOnTrack.index)->clone()); @@ -557,7 +558,6 @@ std::pair MkFitOutputConverter::conver tsosDouble = propagatorOpposite.propagateWithPath(fts, firstHitSurface); } - //return std::make_pair(TrajectoryStateOnSurface(fts, det->surface()), det); return std::make_pair(tsosDouble.first, det); } From 37cb82e43f85a22145e37f465d78c7d655543400 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 20:30:06 +0200 Subject: [PATCH 11/20] Replace hardcoding with calls to mkfit::LayerNumberConverter --- .../MkFit/plugins/MkFitOutputConverter.cc | 122 ++++++++---------- 1 file changed, 51 insertions(+), 71 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 1d38ed91cd8fc..5ac030e3a4539 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -15,6 +15,7 @@ #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/TrackerCommon/interface/TrackerDetSide.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" #include "TrackingTools/Records/interface/TransientRecHitRecord.h" @@ -182,7 +183,7 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: const TrackerTopology& ttopo) const { std::vector dets(lnc.nLayers(), nullptr); - auto set = [&](unsigned int index, DetId id) { + auto setDet = [&](unsigned int index, DetId id) { auto layer = tracker.idToLayer(id); if (layer == nullptr) { throw cms::Exception("LogicError") << "No layer for DetId " << id.rawId(); @@ -193,80 +194,59 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: dets[index] = layer; }; - // TODO: currently hardcoded... - // Logic copied from mkfit::LayerNumberConverter - unsigned int off = 0; - // BPix - set(off + 0, ttopo.pxbDetId(1, 0, 0)); - set(off + 1, ttopo.pxbDetId(2, 0, 0)); - set(off + 2, ttopo.pxbDetId(3, 0, 0)); - set(off + 3, ttopo.pxbDetId(4, 0, 0)); - // offset needs to be increased by 1 here to accommodate the 4th - // pixel barrel layer, while keeping the off+N numbering consistent - // with mkfit::LayerNumberConverter (that, to limited degree, - // supports the layer numbering for phase0 pixel as well) - off += 1; - // TIB - set(off + 3, ttopo.tibDetId(1, 0, 0, 0, 0, 0)); - set(off + 4, ttopo.tibDetId(1, 0, 0, 0, 0, 1)); - set(off + 5, ttopo.tibDetId(2, 0, 0, 0, 0, 0)); - set(off + 6, ttopo.tibDetId(2, 0, 0, 0, 0, 1)); - set(off + 7, ttopo.tibDetId(3, 0, 0, 0, 0, 0)); - set(off + 8, ttopo.tibDetId(4, 0, 0, 0, 0, 0)); - // TOB - set(off + 9, ttopo.tobDetId(1, 0, 0, 0, 0)); - set(off + 10, ttopo.tobDetId(1, 0, 0, 0, 1)); - set(off + 11, ttopo.tobDetId(2, 0, 0, 0, 0)); - set(off + 12, ttopo.tobDetId(2, 0, 0, 0, 1)); - set(off + 13, ttopo.tobDetId(3, 0, 0, 0, 0)); - set(off + 14, ttopo.tobDetId(4, 0, 0, 0, 0)); - set(off + 15, ttopo.tobDetId(5, 0, 0, 0, 0)); - set(off + 16, ttopo.tobDetId(6, 0, 0, 0, 0)); - - auto setForward = [&set, &ttopo](unsigned int off, unsigned int side) { - // FPix - set(off + 0, ttopo.pxfDetId(side, 1, 0, 0, 0)); - set(off + 1, ttopo.pxfDetId(side, 2, 0, 0, 0)); - set(off + 2, ttopo.pxfDetId(side, 3, 0, 0, 0)); - // TID+ - off += 1; // see comment above for barrel - set(off + 2, ttopo.tidDetId(side, 1, 0, 0, 0, 0)); - set(off + 3, ttopo.tidDetId(side, 1, 0, 0, 0, 1)); - set(off + 4, ttopo.tidDetId(side, 2, 0, 0, 0, 0)); - set(off + 5, ttopo.tidDetId(side, 2, 0, 0, 0, 1)); - set(off + 6, ttopo.tidDetId(side, 3, 0, 0, 0, 0)); - set(off + 7, ttopo.tidDetId(side, 3, 0, 0, 0, 1)); - // TEC - set(off + 8, ttopo.tecDetId(side, 1, 0, 0, 0, 0, 0)); - set(off + 9, ttopo.tecDetId(side, 1, 0, 0, 0, 0, 1)); - set(off + 10, ttopo.tecDetId(side, 2, 0, 0, 0, 0, 0)); - set(off + 11, ttopo.tecDetId(side, 2, 0, 0, 0, 0, 1)); - set(off + 12, ttopo.tecDetId(side, 3, 0, 0, 0, 0, 0)); - set(off + 13, ttopo.tecDetId(side, 3, 0, 0, 0, 0, 1)); - set(off + 14, ttopo.tecDetId(side, 4, 0, 0, 0, 0, 0)); - set(off + 15, ttopo.tecDetId(side, 4, 0, 0, 0, 0, 1)); - set(off + 16, ttopo.tecDetId(side, 5, 0, 0, 0, 0, 0)); - set(off + 17, ttopo.tecDetId(side, 5, 0, 0, 0, 0, 1)); - set(off + 18, ttopo.tecDetId(side, 6, 0, 0, 0, 0, 0)); - set(off + 19, ttopo.tecDetId(side, 6, 0, 0, 0, 0, 1)); - set(off + 20, ttopo.tecDetId(side, 7, 0, 0, 0, 0, 0)); - set(off + 21, ttopo.tecDetId(side, 7, 0, 0, 0, 0, 1)); - set(off + 22, ttopo.tecDetId(side, 8, 0, 0, 0, 0, 0)); - set(off + 23, ttopo.tecDetId(side, 8, 0, 0, 0, 0, 1)); - set(off + 24, ttopo.tecDetId(side, 9, 0, 0, 0, 0, 0)); - set(off + 25, ttopo.tecDetId(side, 9, 0, 0, 0, 0, 1)); + auto setBarrel = [&](int det, auto detIdFunc, int layer, int stereo) { + setDet(lnc.convertLayerNumber(det, layer, false, stereo, true), detIdFunc(layer, stereo)); + }; + auto setForward = [&](int det, auto detIdFunc, int disk, int stereo) { + // minus side + setDet(lnc.convertLayerNumber(det, disk, false, stereo, false), + detIdFunc(static_cast(TrackerDetSide::NegEndcap), disk, stereo)); + // plus side + setDet(lnc.convertLayerNumber(det, disk, false, stereo, true), + detIdFunc(static_cast(TrackerDetSide::PosEndcap), disk, stereo)); }; - constexpr unsigned int nlay_barrel = 16 + 1; // 16 in phase0, 1 more in phase1 - constexpr unsigned int nlay_endcap = 25 + 1; // 25 in phase0, 1 more in phase1 + auto pxbDetId = [&ttopo](int layer, int stereo) { return ttopo.pxbDetId(layer, 0, 0); }; + auto tibDetId = [&ttopo](int layer, int stereo) { return ttopo.tibDetId(layer, 0, 0, 0, 0, stereo); }; + auto tobDetId = [&ttopo](int layer, int stereo) { return ttopo.tobDetId(layer, 0, 0, 0, stereo); }; - // plus - off = nlay_barrel + 1; // +1 to move to next slot - setForward(off, 2); + auto pxfDetId = [&ttopo](int side, int disk, int stereo) { return ttopo.pxfDetId(side, disk, 0, 0, 0); }; + auto tidDetId = [&ttopo](int side, int disk, int stereo) { return ttopo.tidDetId(side, disk, 0, 0, 0, stereo); }; + auto tecDetId = [&ttopo](int side, int disk, int stereo) { return ttopo.tecDetId(side, disk, 0, 0, 0, 0, stereo); }; - // minus - off += nlay_endcap + 1; // +1 to move to next slot - setForward(off, 1); + // TODO: detector structure currently hardcoded... + // BPix + for (int layer = 1; layer <= 4; ++layer) { + setBarrel(PixelSubdetector::PixelBarrel, pxbDetId, layer, 0); + } + // TIB + for (int layer = 1; layer <= 4; ++layer) { + setBarrel(StripSubdetector::TIB, tibDetId, layer, 0); + if (layer == 1 or layer == 2) { + setBarrel(StripSubdetector::TIB, tibDetId, layer, 1); + } + } + // TOB + for (int layer = 1; layer <= 6; ++layer) { + setBarrel(StripSubdetector::TOB, tobDetId, layer, 0); + if (layer == 1 or layer == 2) { + setBarrel(StripSubdetector::TOB, tobDetId, layer, 1); + } + } + // FPix + for (int disk = 1; disk <= 3; ++disk) { + setForward(PixelSubdetector::PixelEndcap, pxfDetId, disk, 0); + } + // TID + for (int disk = 1; disk <= 3; ++disk) { + setForward(StripSubdetector::TID, tidDetId, disk, 0); + setForward(StripSubdetector::TID, tidDetId, disk, 1); + } + // TEC + for (int disk = 1; disk <= 9; ++disk) { + setForward(StripSubdetector::TEC, tecDetId, disk, 0); + setForward(StripSubdetector::TEC, tecDetId, disk, 1); + } return dets; } From 0d27626fc7418ec3ac6824c9f2264b425f059511 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 20:34:59 +0200 Subject: [PATCH 12/20] Shorten lines --- RecoTracker/MkFit/python/customizeInitialStepOnly.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/RecoTracker/MkFit/python/customizeInitialStepOnly.py b/RecoTracker/MkFit/python/customizeInitialStepOnly.py index 1b42f710c3e35..9c4e0ef7f4e3d 100644 --- a/RecoTracker/MkFit/python/customizeInitialStepOnly.py +++ b/RecoTracker/MkFit/python/customizeInitialStepOnly.py @@ -57,9 +57,11 @@ def setInput(mtvs, labels): mod.dodEdxPlots = False mod.doResolutionPlotsForLabels = [] - setInput(["trackValidatorTrackingOnly", "trackValidatorAllTPEfficStandalone", "trackValidatorTPPtLess09Standalone", "trackValidatorBHadronTrackingOnly"], + setInput(["trackValidatorTrackingOnly", "trackValidatorAllTPEfficStandalone", + "trackValidatorTPPtLess09Standalone", "trackValidatorBHadronTrackingOnly"], ["cutsRecoTracksInitialStep", "cutsRecoTracksPt09InitialStep"]) - setInput(["trackValidatorFromPVStandalone", "trackValidatorFromPVAllTPStandalone"], ["cutsRecoTracksFromPVInitialStep", "cutsRecoTracksFromPVPt09InitialStep"]) + setInput(["trackValidatorFromPVStandalone", "trackValidatorFromPVAllTPStandalone"], + ["cutsRecoTracksFromPVInitialStep", "cutsRecoTracksFromPVPt09InitialStep"]) setInput(["trackValidatorSeedingTrackingOnly"], ["seedTracksinitialStepSeeds"]) setInput(["trackValidatorBuilding"], ["initialStepTracks"]) process.trackValidatorBuilding.mvaLabels = cms.untracked.PSet(initialStepTracks = cms.untracked.vstring('initialStep')) From b5b6944767e5027f2cada3b88e0ffa9b6c9cd147 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 20:59:17 +0200 Subject: [PATCH 13/20] Deduce layer side from detid, and avoid unnecessarily repeated call to convertLayerNumber() --- RecoTracker/MkFit/plugins/MkFitInputConverter.cc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc index 740c0240c9272..a5f49a5b9d2c4 100644 --- a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -11,6 +11,7 @@ #include "DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/TrackerCommon/interface/TrackerDetSide.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" #include "TrackingTools/Records/interface/TransientRecHitRecord.h" @@ -131,11 +132,15 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, const TrackerTopology& ttopo, const TransientTrackingRecHitBuilder& ttrhBuilder, const mkfit::LayerNumberConverter& lnc) const { + auto isPlusSide = [&ttopo](const DetId& detid) { + return ttopo.side(detid) == static_cast(TrackerDetSide::PosEndcap); + }; 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); + const auto ilay = lnc.convertLayerNumber(subdet, layer, false, isStereo, isPlusSide(detid)); for (const auto& hit : detset) { if (!passCCC(hit, detid)) @@ -152,9 +157,8 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, 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 " + << layer << " isStereo " << isStereo << " zplus " << isPlusSide(detid) << " ilay " << ilay; indexLayers.insert(hit.firstClusterRef().id(), hit.firstClusterRef().index(), mkFitHits[ilay].size(), ilay, &hit); From d685b81b246aaab2bd5b296227dd11fd560afff3 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 21:15:46 +0200 Subject: [PATCH 14/20] Make strip cluster charge cut configurable, default being 'SiStripClusterChargeCutLoose' --- RecoTracker/MkFit/plugins/MkFitInputConverter.cc | 13 ++++++++++--- RecoTracker/MkFit/python/mkFitInputConverter_cfi.py | 10 ++++++++++ 2 files changed, 20 insertions(+), 3 deletions(-) create mode 100644 RecoTracker/MkFit/python/mkFitInputConverter_cfi.py diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc index a5f49a5b9d2c4..c85de4ab943e0 100644 --- a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -71,6 +71,7 @@ class MkFitInputConverter : public edm::global::EDProducer<> { edm::ESGetToken ttopoToken_; edm::ESGetToken mfToken_; edm::EDPutTokenT putToken_; + const float minGoodStripCharge_; }; MkFitInputConverter::MkFitInputConverter(edm::ParameterSet const& iConfig) @@ -84,7 +85,9 @@ MkFitInputConverter::MkFitInputConverter(edm::ParameterSet const& iConfig) iConfig.getParameter("ttrhBuilder"))}, ttopoToken_{esConsumes()}, mfToken_{esConsumes()}, - putToken_{produces()} {} + putToken_{produces()}, + minGoodStripCharge_{static_cast( + iConfig.getParameter("minGoodStripCharge").getParameter("value"))} {} void MkFitInputConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; @@ -95,7 +98,11 @@ void MkFitInputConverter::fillDescriptions(edm::ConfigurationDescriptions& descr desc.add("seeds", edm::InputTag{"initialStepSeeds"}); desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"}); - descriptions.addWithDefaultLabel(desc); + edm::ParameterSetDescription descCCC; + descCCC.add("value"); + desc.add("minGoodStripCharge", descCCC); + + descriptions.add("mkFitInputConverterDefault", desc); } void MkFitInputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { @@ -119,7 +126,7 @@ void MkFitInputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const e } bool MkFitInputConverter::passCCC(const SiStripRecHit2D& hit, const DetId hitId) const { - return (siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster()) >= 1620); + return (siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster()) > minGoodStripCharge_); } bool MkFitInputConverter::passCCC(const SiPixelRecHit& hit, const DetId hitId) const { return true; } diff --git a/RecoTracker/MkFit/python/mkFitInputConverter_cfi.py b/RecoTracker/MkFit/python/mkFitInputConverter_cfi.py new file mode 100644 index 0000000000000..bdd39c281f910 --- /dev/null +++ b/RecoTracker/MkFit/python/mkFitInputConverter_cfi.py @@ -0,0 +1,10 @@ +import FWCore.ParameterSet.Config as cms + +from RecoTracker.MkFit.mkFitInputConverterDefault_cfi import mkFitInputConverterDefault as _mkFitInputConverterDefault +from RecoLocalTracker.SiStripClusterizer.SiStripClusterChargeCut_cfi import * + +mkFitInputConverter = _mkFitInputConverterDefault.clone( + minGoodStripCharge = cms.PSet( + refToPSet_ = cms.string('SiStripClusterChargeCutLoose') + ) +) From 7352bc2eb7263385772e79ffae30a5ca6c8c86a6 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 23:07:58 +0200 Subject: [PATCH 15/20] Preallocate memory for hit index mapping structure --- RecoTracker/MkFit/interface/MkFitIndexLayer.h | 2 ++ .../MkFit/plugins/MkFitInputConverter.cc | 17 ++++++++- RecoTracker/MkFit/src/MkFitIndexLayer.cc | 35 +++++++++++++++---- 3 files changed, 46 insertions(+), 8 deletions(-) diff --git a/RecoTracker/MkFit/interface/MkFitIndexLayer.h b/RecoTracker/MkFit/interface/MkFitIndexLayer.h index 6482809b96b4c..44f0fd6a33cca 100644 --- a/RecoTracker/MkFit/interface/MkFitIndexLayer.h +++ b/RecoTracker/MkFit/interface/MkFitIndexLayer.h @@ -24,6 +24,8 @@ class MkFitIndexLayer { MkFitIndexLayer() = default; + void resizeByClusterIndex(edm::ProductID id, size_t clusterIndex); + void increaseLayerSize(int layer, size_t additionalSize); 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; diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc index c85de4ab943e0..c49d733d2b5f8 100644 --- a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -115,9 +115,10 @@ void MkFitInputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const e std::vector 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); + // Process strips first for better memory allocation pattern convertHits(iEvent.get(stripRphiRecHitToken_), mkFitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc); convertHits(iEvent.get(stripStereoRecHitToken_), mkFitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc); + convertHits(iEvent.get(pixelRecHitToken_), mkFitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc); // Then import seeds auto mkFitSeeds = convertSeeds(iEvent.get(seedToken_), indexLayers, ttrhBuilder, iSetup.getData(mfToken_)); @@ -139,15 +140,29 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, const TrackerTopology& ttopo, const TransientTrackingRecHitBuilder& ttrhBuilder, const mkfit::LayerNumberConverter& lnc) const { + if (hits.empty()) + return; auto isPlusSide = [&ttopo](const DetId& detid) { return ttopo.side(detid) == static_cast(TrackerDetSide::PosEndcap); }; + + { + const DetId detid{hits.ids().back()}; + const auto ilay = + lnc.convertLayerNumber(detid.subdetId(), ttopo.layer(detid), false, ttopo.isStereo(detid), isPlusSide(detid)); + // Do initial reserves to minimize further memory allocations + const auto& lastClusterRef = hits.data().back().firstClusterRef(); + indexLayers.resizeByClusterIndex(lastClusterRef.id(), lastClusterRef.index()); + indexLayers.increaseLayerSize(ilay, hits.detsetSize(hits.ids().size() - 1)); + } + 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); const auto ilay = lnc.convertLayerNumber(subdet, layer, false, isStereo, isPlusSide(detid)); + indexLayers.increaseLayerSize(ilay, detset.size()); // to minimize memory allocations for (const auto& hit : detset) { if (!passCCC(hit, detid)) diff --git a/RecoTracker/MkFit/src/MkFitIndexLayer.cc b/RecoTracker/MkFit/src/MkFitIndexLayer.cc index a9b8e091297f4..44778b27bc4e5 100644 --- a/RecoTracker/MkFit/src/MkFitIndexLayer.cc +++ b/RecoTracker/MkFit/src/MkFitIndexLayer.cc @@ -4,18 +4,39 @@ #include -void MkFitIndexLayer::insert(edm::ProductID id, size_t clusterIndex, int hit, int layer, const TrackingRecHit* hitPtr) { - // mapping CMSSW->mkfit - auto found = std::find_if(colls_.begin(), colls_.end(), [&](const auto& item) { return item.productID == id; }); - if (found == colls_.end()) { - found = colls_.emplace(colls_.end(), id); +namespace { + template + auto resizeByClusterIndexImpl(T& cmsswToMkFit, edm::ProductID id, size_t clusterIndex) -> typename T::iterator { + auto found = std::find_if(cmsswToMkFit.begin(), cmsswToMkFit.end(), [&](const auto& item) { return item.productID == id; }); + if (found == cmsswToMkFit.end()) { + found = cmsswToMkFit.emplace(cmsswToMkFit.end(), id); + } + if (found->infos.size() <= clusterIndex) { + found->infos.resize(clusterIndex + 1); + } + return found; } - if (found->infos.size() <= clusterIndex) { - found->infos.resize(clusterIndex + 1); +} + +void MkFitIndexLayer::resizeByClusterIndex(edm::ProductID id, size_t clusterIndex) { + resizeByClusterIndexImpl(colls_, id, clusterIndex); +} + +void MkFitIndexLayer::increaseLayerSize(int layer, size_t additionalSize) { + if (layer >= static_cast(hits_.size())) { + hits_.resize(layer + 1); } + hits_[layer].resize(hits_[layer].size() + additionalSize); +} + +void MkFitIndexLayer::insert(edm::ProductID id, size_t clusterIndex, int hit, int layer, const TrackingRecHit* hitPtr) { + // mapping CMSSW->mkfit + auto found = resizeByClusterIndexImpl(colls_, id, clusterIndex); found->infos[clusterIndex] = HitInfo(hit, layer); // mapping mkfit->CMSSW + // when client calls increaseLayerSize() the two checks below are + // redundant, but better to keep them if (layer >= static_cast(hits_.size())) { hits_.resize(layer + 1); } From 48e563fac25455ed89364303a127428a491d3261 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Wed, 4 Sep 2019 23:38:17 +0200 Subject: [PATCH 16/20] Rename MkFitIndexLayer to MkFitHitIndexMap --- .../{MkFitIndexLayer.h => MkFitHitIndexMap.h} | 8 ++--- .../MkFit/interface/MkFitInputWrapper.h | 8 ++--- .../MkFit/plugins/MkFitInputConverter.cc | 30 +++++++++---------- .../MkFit/plugins/MkFitOutputConverter.cc | 10 +++---- ...MkFitIndexLayer.cc => MkFitHitIndexMap.cc} | 10 +++---- RecoTracker/MkFit/src/MkFitInputWrapper.cc | 4 +-- 6 files changed, 35 insertions(+), 35 deletions(-) rename RecoTracker/MkFit/interface/{MkFitIndexLayer.h => MkFitHitIndexMap.h} (89%) rename RecoTracker/MkFit/src/{MkFitIndexLayer.cc => MkFitHitIndexMap.cc} (81%) diff --git a/RecoTracker/MkFit/interface/MkFitIndexLayer.h b/RecoTracker/MkFit/interface/MkFitHitIndexMap.h similarity index 89% rename from RecoTracker/MkFit/interface/MkFitIndexLayer.h rename to RecoTracker/MkFit/interface/MkFitHitIndexMap.h index 44f0fd6a33cca..f59a242f6e6dd 100644 --- a/RecoTracker/MkFit/interface/MkFitIndexLayer.h +++ b/RecoTracker/MkFit/interface/MkFitHitIndexMap.h @@ -1,5 +1,5 @@ -#ifndef RecoTracker_MkFit_MkFitIndexLayer_h -#define RecoTracker_MkFit_MkFitIndexLayer_h +#ifndef RecoTracker_MkFit_MkFitHitIndexMap_h +#define RecoTracker_MkFit_MkFitHitIndexMap_h #include "DataFormats/Provenance/interface/ProductID.h" @@ -7,7 +7,7 @@ class TrackingRecHit; -class MkFitIndexLayer { +class MkFitHitIndexMap { public: struct HitInfo { HitInfo() : index(-1), layer(-1) {} @@ -22,7 +22,7 @@ class MkFitIndexLayer { std::vector infos; // indexed by cluster index }; - MkFitIndexLayer() = default; + MkFitHitIndexMap() = default; void resizeByClusterIndex(edm::ProductID id, size_t clusterIndex); void increaseLayerSize(int layer, size_t additionalSize); diff --git a/RecoTracker/MkFit/interface/MkFitInputWrapper.h b/RecoTracker/MkFit/interface/MkFitInputWrapper.h index 2645716b9efda..ae17cf6e32584 100644 --- a/RecoTracker/MkFit/interface/MkFitInputWrapper.h +++ b/RecoTracker/MkFit/interface/MkFitInputWrapper.h @@ -1,7 +1,7 @@ #ifndef RecoTracker_MkFit_MkFitInputWrapper_h #define RecoTracker_MkFit_MkFitInputWrapper_h -#include "RecoTracker/MkFit/interface/MkFitIndexLayer.h" +#include "RecoTracker/MkFit/interface/MkFitHitIndexMap.h" #include #include @@ -17,7 +17,7 @@ namespace mkfit { class MkFitInputWrapper { public: MkFitInputWrapper(); - MkFitInputWrapper(MkFitIndexLayer&& indexLayers, + MkFitInputWrapper(MkFitHitIndexMap&& hitIndexMap, std::vector&& hits, mkfit::TrackVec&& seeds, mkfit::LayerNumberConverter&& lnc); @@ -28,14 +28,14 @@ class MkFitInputWrapper { MkFitInputWrapper(MkFitInputWrapper&&); MkFitInputWrapper& operator=(MkFitInputWrapper&&); - MkFitIndexLayer const& indexLayers() const { return indexLayers_; } + MkFitHitIndexMap const& hitIndexMap() const { return hitIndexMap_; } mkfit::TrackVec const& seeds() const { return *seeds_; } std::vector const& hits() const { return hits_; } mkfit::LayerNumberConverter const& layerNumberConverter() const { return *lnc_; } unsigned int nlayers() const; private: - MkFitIndexLayer indexLayers_; + MkFitHitIndexMap hitIndexMap_; std::vector hits_; std::unique_ptr seeds_; // for pimpl pattern std::unique_ptr lnc_; // for pimpl pattern diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc index c49d733d2b5f8..7b7245623a581 100644 --- a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -45,7 +45,7 @@ class MkFitInputConverter : public edm::global::EDProducer<> { template void convertHits(const HitCollection& hits, std::vector& mkFitHits, - MkFitIndexLayer& indexLayers, + MkFitHitIndexMap& hitIndexMap, int& totalHits, const TrackerTopology& ttopo, const TransientTrackingRecHitBuilder& ttrhBuilder, @@ -55,7 +55,7 @@ class MkFitInputConverter : public edm::global::EDProducer<> { bool passCCC(const SiPixelRecHit& hit, const DetId hitId) const; mkfit::TrackVec convertSeeds(const edm::View& seeds, - const MkFitIndexLayer& indexLayers, + const MkFitHitIndexMap& hitIndexMap, const TransientTrackingRecHitBuilder& ttrhBuilder, const MagneticField& mf) const; @@ -113,17 +113,17 @@ void MkFitInputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const e const auto& ttopo = iSetup.getData(ttopoToken_); std::vector mkFitHits(lnc.nLayers()); - MkFitIndexLayer indexLayers; + MkFitHitIndexMap hitIndexMap; int totalHits = 0; // I need to have a global hit index in order to have the hit remapping working? // Process strips first for better memory allocation pattern - convertHits(iEvent.get(stripRphiRecHitToken_), mkFitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc); - convertHits(iEvent.get(stripStereoRecHitToken_), mkFitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc); - convertHits(iEvent.get(pixelRecHitToken_), mkFitHits, indexLayers, totalHits, ttopo, ttrhBuilder, lnc); + convertHits(iEvent.get(stripRphiRecHitToken_), mkFitHits, hitIndexMap, totalHits, ttopo, ttrhBuilder, lnc); + convertHits(iEvent.get(stripStereoRecHitToken_), mkFitHits, hitIndexMap, totalHits, ttopo, ttrhBuilder, lnc); + convertHits(iEvent.get(pixelRecHitToken_), mkFitHits, hitIndexMap, totalHits, ttopo, ttrhBuilder, lnc); // Then import seeds - auto mkFitSeeds = convertSeeds(iEvent.get(seedToken_), indexLayers, ttrhBuilder, iSetup.getData(mfToken_)); + auto mkFitSeeds = convertSeeds(iEvent.get(seedToken_), hitIndexMap, ttrhBuilder, iSetup.getData(mfToken_)); - iEvent.emplace(putToken_, std::move(indexLayers), std::move(mkFitHits), std::move(mkFitSeeds), std::move(lnc)); + iEvent.emplace(putToken_, std::move(hitIndexMap), std::move(mkFitHits), std::move(mkFitSeeds), std::move(lnc)); } bool MkFitInputConverter::passCCC(const SiStripRecHit2D& hit, const DetId hitId) const { @@ -135,7 +135,7 @@ bool MkFitInputConverter::passCCC(const SiPixelRecHit& hit, const DetId hitId) c template void MkFitInputConverter::convertHits(const HitCollection& hits, std::vector& mkFitHits, - MkFitIndexLayer& indexLayers, + MkFitHitIndexMap& hitIndexMap, int& totalHits, const TrackerTopology& ttopo, const TransientTrackingRecHitBuilder& ttrhBuilder, @@ -152,8 +152,8 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, lnc.convertLayerNumber(detid.subdetId(), ttopo.layer(detid), false, ttopo.isStereo(detid), isPlusSide(detid)); // Do initial reserves to minimize further memory allocations const auto& lastClusterRef = hits.data().back().firstClusterRef(); - indexLayers.resizeByClusterIndex(lastClusterRef.id(), lastClusterRef.index()); - indexLayers.increaseLayerSize(ilay, hits.detsetSize(hits.ids().size() - 1)); + hitIndexMap.resizeByClusterIndex(lastClusterRef.id(), lastClusterRef.index()); + hitIndexMap.increaseLayerSize(ilay, hits.detsetSize(hits.ids().size() - 1)); } for (const auto& detset : hits) { @@ -162,7 +162,7 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, const auto layer = ttopo.layer(detid); const auto isStereo = ttopo.isStereo(detid); const auto ilay = lnc.convertLayerNumber(subdet, layer, false, isStereo, isPlusSide(detid)); - indexLayers.increaseLayerSize(ilay, detset.size()); // to minimize memory allocations + hitIndexMap.increaseLayerSize(ilay, detset.size()); // to minimize memory allocations for (const auto& hit : detset) { if (!passCCC(hit, detid)) @@ -183,7 +183,7 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, << layer << " isStereo " << isStereo << " zplus " << isPlusSide(detid) << " ilay " << ilay; - indexLayers.insert(hit.firstClusterRef().id(), hit.firstClusterRef().index(), mkFitHits[ilay].size(), ilay, &hit); + hitIndexMap.insert(hit.firstClusterRef().id(), hit.firstClusterRef().index(), mkFitHits[ilay].size(), ilay, &hit); mkFitHits[ilay].emplace_back(pos, err, totalHits); ++totalHits; } @@ -191,7 +191,7 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, } mkfit::TrackVec MkFitInputConverter::convertSeeds(const edm::View& seeds, - const MkFitIndexLayer& indexLayers, + const MkFitHitIndexMap& hitIndexMap, const TransientTrackingRecHitBuilder& ttrhBuilder, const MagneticField& mf) const { mkfit::TrackVec ret; @@ -226,7 +226,7 @@ mkfit::TrackVec MkFitInputConverter::convertSeeds(const edm::View(*iHit).firstClusterRef(); - const auto& info = indexLayers.get(clusterRef.id(), clusterRef.index()); + const auto& info = hitIndexMap.get(clusterRef.id(), clusterRef.index()); ret.back().addHitIdx(info.index, info.layer, 0); // per-hit chi2 is not known } ++index; diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 5ac030e3a4539..d967ddfc3127c 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -73,7 +73,7 @@ class MkFitOutputConverter : public edm::global::EDProducer<> { const TrackerTopology& ttopo) const; TrackCandidateCollection convertCandidates(const MkFitOutputWrapper& mkFitOutput, - const MkFitIndexLayer& indexLayers, + const MkFitHitIndexMap& hitIndexMap, const edm::View& seeds, const TrackerGeometry& geom, const MagneticField& mf, @@ -164,7 +164,7 @@ void MkFitOutputConverter::produce(edm::StreamID iID, edm::Event& iEvent, const createDetLayers(hitsSeeds.layerNumberConverter(), *(mte.geometricSearchTracker()), iSetup.getData(ttopoToken_)); iEvent.emplace(putTrackCandidateToken_, convertCandidates(iEvent.get(tracksToken_), - hitsSeeds.indexLayers(), + hitsSeeds.hitIndexMap(), seeds, iSetup.getData(geomToken_), iSetup.getData(mfToken_), @@ -252,7 +252,7 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: } TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutputWrapper& mkFitOutput, - const MkFitIndexLayer& indexLayers, + const MkFitHitIndexMap& hitIndexMap, const edm::View& seeds, const TrackerGeometry& geom, const MagneticField& mf, @@ -303,13 +303,13 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp // them in the TrackProducer. lastHitInvalid = true; } else { - recHits.push_back(indexLayers.getHitPtr(hitOnTrack.layer, hitOnTrack.index)->clone()); + recHits.push_back(hitIndexMap.getHitPtr(hitOnTrack.layer, hitOnTrack.index)->clone()); LogTrace("MkFitOutputConverter") << " pos " << recHits.back().globalPosition().x() << " " << recHits.back().globalPosition().y() << " " << recHits.back().globalPosition().z() << " mag2 " << recHits.back().globalPosition().mag2() << " detid " << recHits.back().geographicalId().rawId() << " cluster " - << indexLayers.getClusterIndex(hitOnTrack.layer, hitOnTrack.index); + << hitIndexMap.getClusterIndex(hitOnTrack.layer, hitOnTrack.index); lastHitInvalid = false; } } diff --git a/RecoTracker/MkFit/src/MkFitIndexLayer.cc b/RecoTracker/MkFit/src/MkFitHitIndexMap.cc similarity index 81% rename from RecoTracker/MkFit/src/MkFitIndexLayer.cc rename to RecoTracker/MkFit/src/MkFitHitIndexMap.cc index 44778b27bc4e5..f68dee6a83c85 100644 --- a/RecoTracker/MkFit/src/MkFitIndexLayer.cc +++ b/RecoTracker/MkFit/src/MkFitHitIndexMap.cc @@ -1,4 +1,4 @@ -#include "RecoTracker/MkFit/interface/MkFitIndexLayer.h" +#include "RecoTracker/MkFit/interface/MkFitHitIndexMap.h" #include "FWCore/Utilities/interface/Exception.h" @@ -18,18 +18,18 @@ namespace { } } -void MkFitIndexLayer::resizeByClusterIndex(edm::ProductID id, size_t clusterIndex) { +void MkFitHitIndexMap::resizeByClusterIndex(edm::ProductID id, size_t clusterIndex) { resizeByClusterIndexImpl(colls_, id, clusterIndex); } -void MkFitIndexLayer::increaseLayerSize(int layer, size_t additionalSize) { +void MkFitHitIndexMap::increaseLayerSize(int layer, size_t additionalSize) { if (layer >= static_cast(hits_.size())) { hits_.resize(layer + 1); } hits_[layer].resize(hits_[layer].size() + additionalSize); } -void MkFitIndexLayer::insert(edm::ProductID id, size_t clusterIndex, int hit, int layer, const TrackingRecHit* hitPtr) { +void MkFitHitIndexMap::insert(edm::ProductID id, size_t clusterIndex, int hit, int layer, const TrackingRecHit* hitPtr) { // mapping CMSSW->mkfit auto found = resizeByClusterIndexImpl(colls_, id, clusterIndex); found->infos[clusterIndex] = HitInfo(hit, layer); @@ -47,7 +47,7 @@ void MkFitIndexLayer::insert(edm::ProductID id, size_t clusterIndex, int hit, in hits_[layer][hit].clusterIndex = clusterIndex; } -const MkFitIndexLayer::HitInfo& MkFitIndexLayer::get(edm::ProductID id, size_t clusterIndex) const { +const MkFitHitIndexMap::HitInfo& MkFitHitIndexMap::get(edm::ProductID id, size_t clusterIndex) const { auto found = std::find_if(colls_.begin(), colls_.end(), [&](const auto& item) { return item.productID == id; }); if (found == colls_.end()) { auto exp = cms::Exception("Assert"); diff --git a/RecoTracker/MkFit/src/MkFitInputWrapper.cc b/RecoTracker/MkFit/src/MkFitInputWrapper.cc index 54ee931bbc5ac..6a4806efd2dd9 100644 --- a/RecoTracker/MkFit/src/MkFitInputWrapper.cc +++ b/RecoTracker/MkFit/src/MkFitInputWrapper.cc @@ -7,11 +7,11 @@ MkFitInputWrapper::MkFitInputWrapper() = default; -MkFitInputWrapper::MkFitInputWrapper(MkFitIndexLayer&& indexLayers, +MkFitInputWrapper::MkFitInputWrapper(MkFitHitIndexMap&& hitIndexMap, std::vector&& hits, mkfit::TrackVec&& seeds, mkfit::LayerNumberConverter&& lnc) - : indexLayers_{std::move(indexLayers)}, + : hitIndexMap_{std::move(hitIndexMap)}, hits_{std::move(hits)}, seeds_{std::make_unique(std::move(seeds))}, lnc_{std::make_unique(std::move(lnc))} {} From 5bee1f587a6b415d61029382fe2480dd913976aa Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Thu, 5 Sep 2019 00:31:31 +0200 Subject: [PATCH 17/20] Rename MkFitHitIndexMap internals and interface for better self-documentation --- .../MkFit/interface/MkFitHitIndexMap.h | 73 ++++++++++++++----- .../MkFit/plugins/MkFitInputConverter.cc | 9 ++- .../MkFit/plugins/MkFitOutputConverter.cc | 5 +- RecoTracker/MkFit/src/MkFitHitIndexMap.cc | 48 ++++++------ 4 files changed, 90 insertions(+), 45 deletions(-) diff --git a/RecoTracker/MkFit/interface/MkFitHitIndexMap.h b/RecoTracker/MkFit/interface/MkFitHitIndexMap.h index f59a242f6e6dd..c7b6eea50263a 100644 --- a/RecoTracker/MkFit/interface/MkFitHitIndexMap.h +++ b/RecoTracker/MkFit/interface/MkFitHitIndexMap.h @@ -7,41 +7,80 @@ class TrackingRecHit; +/** + * This class provides mappings + * - from CMSSW(ProductID, cluster index) to mkFit(layer index, hit index) + * - from mkFit(layer index, hit index) to pointer to CMSSW hit + */ class MkFitHitIndexMap { public: - struct HitInfo { - HitInfo() : index(-1), layer(-1) {} - HitInfo(int i, int l) : index(i), layer(l) {} - int index; - int layer; - }; + // This class holds the index and layer of a hit in the hit data + // structure passed to mkFit + class MkFitHit { + public: + MkFitHit() = default; + explicit MkFitHit(int i, int l) : index_{i}, layer_{l} {} - struct Coll { - explicit Coll(edm::ProductID id) : productID(id) {} - edm::ProductID productID; - std::vector infos; // indexed by cluster index + int index() const { return index_; } + int layer() const { return layer_; } + + private: + int index_ = -1; + int layer_ = -1; }; MkFitHitIndexMap() = default; + /** + * Can be used to preallocate the internal vectors for CMSSW->mkFit mapping + */ void resizeByClusterIndex(edm::ProductID id, size_t clusterIndex); + + /** + * Can be used to preallocate the internal vectors for mkFit->CMSSW mapping + * + * \param layer Layer index (in mkFit convention) + * \param additionalSize Number of additional elements to make space for + */ void increaseLayerSize(int layer, size_t additionalSize); - 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; + /** + * Inserts a new hit in the mapping + * + * \param id ProductID of the cluster collection + * \param clusterIndex Index of the cluster in the cluster collection + * \param hit Index and layer of the hit in the mkFit hit data structure + * \param hitPtr Pointer to the TrackingRecHit + */ + void insert(edm::ProductID id, size_t clusterIndex, MkFitHit hit, const TrackingRecHit* hitPtr); - const TrackingRecHit *getHitPtr(int layer, int hit) const { return hits_.at(layer).at(hit).ptr; } + /// Get mkFit hit index and layer + const MkFitHit& mkFitHit(edm::ProductID id, size_t clusterIndex) const; - size_t getClusterIndex(int layer, int hit) const { return hits_.at(layer).at(hit).clusterIndex; } + /// Get CMSSW hit pointer + const TrackingRecHit* hitPtr(MkFitHit hit) const { return mkFitToCMSSW_.at(hit.layer()).at(hit.index()).ptr; } + + /// Get CMSSW cluster index (currently used only for debugging) + size_t clusterIndex(MkFitHit hit) const { return mkFitToCMSSW_.at(hit.layer()).at(hit.index()).clusterIndex; } private: + // Helper struct to map (edm::ProductID, cluster index) to MkFitHit + struct ClusterToMkFitHit { + explicit ClusterToMkFitHit(edm::ProductID id) : productID(id) {} + edm::ProductID productID; + std::vector mkFitHits; // indexed by cluster index + }; + + // Helper struct to map MkFitHit to (TrackingRecHit *, cluster index) struct CMSSWHit { - const TrackingRecHit *ptr = nullptr; + CMSSWHit() = default; + explicit CMSSWHit(const TrackingRecHit* p, size_t i) : ptr{p}, clusterIndex{i} {} + const TrackingRecHit* ptr = nullptr; size_t clusterIndex = 0; }; - std::vector colls_; // mapping from CMSSW(ProductID, index) -> mkFit(index, layer) - std::vector > hits_; // reverse mapping, mkFit(layer, index) -> CMSSW hit + std::vector cmsswToMkFit_; // mapping from CMSSW(ProductID, cluster index) -> mkFit(index, layer) + std::vector > mkFitToCMSSW_; // reverse mapping, mkFit(layer, index) -> CMSSW hit }; #endif diff --git a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc index 7b7245623a581..499054817c6c8 100644 --- a/RecoTracker/MkFit/plugins/MkFitInputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -183,7 +183,10 @@ void MkFitInputConverter::convertHits(const HitCollection& hits, << layer << " isStereo " << isStereo << " zplus " << isPlusSide(detid) << " ilay " << ilay; - hitIndexMap.insert(hit.firstClusterRef().id(), hit.firstClusterRef().index(), mkFitHits[ilay].size(), ilay, &hit); + hitIndexMap.insert(hit.firstClusterRef().id(), + hit.firstClusterRef().index(), + MkFitHitIndexMap::MkFitHit{static_cast(mkFitHits[ilay].size()), ilay}, + &hit); mkFitHits[ilay].emplace_back(pos, err, totalHits); ++totalHits; } @@ -226,8 +229,8 @@ mkfit::TrackVec MkFitInputConverter::convertSeeds(const edm::View(*iHit).firstClusterRef(); - const auto& info = hitIndexMap.get(clusterRef.id(), clusterRef.index()); - ret.back().addHitIdx(info.index, info.layer, 0); // per-hit chi2 is not known + const auto& mkFitHit = hitIndexMap.mkFitHit(clusterRef.id(), clusterRef.index()); + ret.back().addHitIdx(mkFitHit.index(), mkFitHit.layer(), 0); // per-hit chi2 is not known } ++index; } diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index d967ddfc3127c..2ff59a2689d61 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -303,13 +303,14 @@ TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutp // them in the TrackProducer. lastHitInvalid = true; } else { - recHits.push_back(hitIndexMap.getHitPtr(hitOnTrack.layer, hitOnTrack.index)->clone()); + recHits.push_back(hitIndexMap.hitPtr(MkFitHitIndexMap::MkFitHit{hitOnTrack.index, hitOnTrack.layer})->clone()); LogTrace("MkFitOutputConverter") << " pos " << recHits.back().globalPosition().x() << " " << recHits.back().globalPosition().y() << " " << recHits.back().globalPosition().z() << " mag2 " << recHits.back().globalPosition().mag2() << " detid " << recHits.back().geographicalId().rawId() << " cluster " - << hitIndexMap.getClusterIndex(hitOnTrack.layer, hitOnTrack.index); + << hitIndexMap.clusterIndex( + MkFitHitIndexMap::MkFitHit{hitOnTrack.index, hitOnTrack.layer}); lastHitInvalid = false; } } diff --git a/RecoTracker/MkFit/src/MkFitHitIndexMap.cc b/RecoTracker/MkFit/src/MkFitHitIndexMap.cc index f68dee6a83c85..89260eadbe565 100644 --- a/RecoTracker/MkFit/src/MkFitHitIndexMap.cc +++ b/RecoTracker/MkFit/src/MkFitHitIndexMap.cc @@ -7,59 +7,61 @@ namespace { template auto resizeByClusterIndexImpl(T& cmsswToMkFit, edm::ProductID id, size_t clusterIndex) -> typename T::iterator { - auto found = std::find_if(cmsswToMkFit.begin(), cmsswToMkFit.end(), [&](const auto& item) { return item.productID == id; }); + auto found = + std::find_if(cmsswToMkFit.begin(), cmsswToMkFit.end(), [&](const auto& item) { return item.productID == id; }); if (found == cmsswToMkFit.end()) { found = cmsswToMkFit.emplace(cmsswToMkFit.end(), id); } - if (found->infos.size() <= clusterIndex) { - found->infos.resize(clusterIndex + 1); + if (found->mkFitHits.size() <= clusterIndex) { + found->mkFitHits.resize(clusterIndex + 1); } return found; } -} +} // namespace void MkFitHitIndexMap::resizeByClusterIndex(edm::ProductID id, size_t clusterIndex) { - resizeByClusterIndexImpl(colls_, id, clusterIndex); + resizeByClusterIndexImpl(cmsswToMkFit_, id, clusterIndex); } void MkFitHitIndexMap::increaseLayerSize(int layer, size_t additionalSize) { - if (layer >= static_cast(hits_.size())) { - hits_.resize(layer + 1); + if (layer >= static_cast(mkFitToCMSSW_.size())) { + mkFitToCMSSW_.resize(layer + 1); } - hits_[layer].resize(hits_[layer].size() + additionalSize); + mkFitToCMSSW_[layer].resize(mkFitToCMSSW_[layer].size() + additionalSize); } -void MkFitHitIndexMap::insert(edm::ProductID id, size_t clusterIndex, int hit, int layer, const TrackingRecHit* hitPtr) { +void MkFitHitIndexMap::insert(edm::ProductID id, size_t clusterIndex, MkFitHit hit, const TrackingRecHit* hitPtr) { // mapping CMSSW->mkfit - auto found = resizeByClusterIndexImpl(colls_, id, clusterIndex); - found->infos[clusterIndex] = HitInfo(hit, layer); + auto found = resizeByClusterIndexImpl(cmsswToMkFit_, id, clusterIndex); + found->mkFitHits[clusterIndex] = hit; // mapping mkfit->CMSSW // when client calls increaseLayerSize() the two checks below are // redundant, but better to keep them - if (layer >= static_cast(hits_.size())) { - hits_.resize(layer + 1); + if (hit.layer() >= static_cast(mkFitToCMSSW_.size())) { + mkFitToCMSSW_.resize(hit.layer() + 1); } - if (hit >= static_cast(hits_[layer].size())) { - hits_[layer].resize(hit + 1); + auto& layer = mkFitToCMSSW_[hit.layer()]; + if (hit.index() >= static_cast(layer.size())) { + layer.resize(hit.index() + 1); } - hits_[layer][hit].ptr = hitPtr; - hits_[layer][hit].clusterIndex = clusterIndex; + layer[hit.index()] = CMSSWHit(hitPtr, clusterIndex); } -const MkFitHitIndexMap::HitInfo& MkFitHitIndexMap::get(edm::ProductID id, size_t clusterIndex) const { - auto found = std::find_if(colls_.begin(), colls_.end(), [&](const auto& item) { return item.productID == id; }); - if (found == colls_.end()) { +const MkFitHitIndexMap::MkFitHit& MkFitHitIndexMap::mkFitHit(edm::ProductID id, size_t clusterIndex) const { + auto found = + std::find_if(cmsswToMkFit_.begin(), cmsswToMkFit_.end(), [&](const auto& item) { return item.productID == id; }); + if (found == cmsswToMkFit_.end()) { auto exp = cms::Exception("Assert"); exp << "Encountered a seed with a hit having productID " << id << " which is not any of the input hit collections: "; - for (const auto& elem : colls_) { + for (const auto& elem : cmsswToMkFit_) { exp << elem.productID << " "; } throw exp; } - const HitInfo& ret = found->infos[clusterIndex]; - if (ret.index < 0) { + const MkFitHit& ret = found->mkFitHits.at(clusterIndex); + if (ret.index() < 0) { throw cms::Exception("Assert") << "No hit index for cluster " << clusterIndex << " of collection " << id; } return ret; From c26b3d9454bf5360b186e7c5168fb7c56d35ba50 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Mon, 9 Sep 2019 23:21:42 +0200 Subject: [PATCH 18/20] Simplify the DetLayer vector creation even more --- .../MkFit/plugins/MkFitOutputConverter.cc | 79 +++++-------------- 1 file changed, 19 insertions(+), 60 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 2ff59a2689d61..8406853d02f9b 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -183,70 +183,29 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: const TrackerTopology& ttopo) const { std::vector dets(lnc.nLayers(), nullptr); - auto setDet = [&](unsigned int index, DetId id) { - auto layer = tracker.idToLayer(id); - if (layer == nullptr) { - throw cms::Exception("LogicError") << "No layer for DetId " << id.rawId(); - } - LogTrace("MkFitOutputConverter") << "Setting DetLayer for index " << index << " subdet " << id.subdetId() - << " layer " << ttopo.layer(id) << " ptr " << layer; - - dets[index] = layer; - }; - - auto setBarrel = [&](int det, auto detIdFunc, int layer, int stereo) { - setDet(lnc.convertLayerNumber(det, layer, false, stereo, true), detIdFunc(layer, stereo)); + auto isPlusSide = [&ttopo](const DetId& detid) { + return ttopo.side(detid) == static_cast(TrackerDetSide::PosEndcap); }; - auto setForward = [&](int det, auto detIdFunc, int disk, int stereo) { - // minus side - setDet(lnc.convertLayerNumber(det, disk, false, stereo, false), - detIdFunc(static_cast(TrackerDetSide::NegEndcap), disk, stereo)); - // plus side - setDet(lnc.convertLayerNumber(det, disk, false, stereo, true), - detIdFunc(static_cast(TrackerDetSide::PosEndcap), disk, stereo)); - }; - - auto pxbDetId = [&ttopo](int layer, int stereo) { return ttopo.pxbDetId(layer, 0, 0); }; - auto tibDetId = [&ttopo](int layer, int stereo) { return ttopo.tibDetId(layer, 0, 0, 0, 0, stereo); }; - auto tobDetId = [&ttopo](int layer, int stereo) { return ttopo.tobDetId(layer, 0, 0, 0, stereo); }; - - auto pxfDetId = [&ttopo](int side, int disk, int stereo) { return ttopo.pxfDetId(side, disk, 0, 0, 0); }; - auto tidDetId = [&ttopo](int side, int disk, int stereo) { return ttopo.tidDetId(side, disk, 0, 0, 0, stereo); }; - auto tecDetId = [&ttopo](int side, int disk, int stereo) { return ttopo.tecDetId(side, disk, 0, 0, 0, 0, stereo); }; - - // TODO: detector structure currently hardcoded... - // BPix - for (int layer = 1; layer <= 4; ++layer) { - setBarrel(PixelSubdetector::PixelBarrel, pxbDetId, layer, 0); - } - // TIB - for (int layer = 1; layer <= 4; ++layer) { - setBarrel(StripSubdetector::TIB, tibDetId, layer, 0); - if (layer == 1 or layer == 2) { - setBarrel(StripSubdetector::TIB, tibDetId, layer, 1); + constexpr int isMono = 0; + constexpr int isStereo = 1; + for (const DetLayer* lay : tracker.allLayers()) { + const auto& comp = lay->basicComponents(); + if (UNLIKELY(comp.empty())) { + throw cms::Exception("LogicError") << "Got a tracker layer (subdet " << lay->subDetector() + << ") with empty basicComponents."; } - } - // TOB - for (int layer = 1; layer <= 6; ++layer) { - setBarrel(StripSubdetector::TOB, tobDetId, layer, 0); - if (layer == 1 or layer == 2) { - setBarrel(StripSubdetector::TOB, tobDetId, layer, 1); + // First component is enough for layer and side information + const auto& detId = comp.front()->geographicalId(); + const auto subdet = detId.subdetId(); + const auto layer = ttopo.layer(detId); + + // TODO: mono/stereo structure is still hardcoded for phase0/1 strip tracker + dets[lnc.convertLayerNumber(subdet, layer, false, isMono, isPlusSide(detId))] = lay; + if (((subdet == StripSubdetector::TIB or StripSubdetector::TOB) and (layer == 1 or layer == 2)) or + subdet == StripSubdetector::TID or subdet == StripSubdetector::TEC) { + dets[lnc.convertLayerNumber(subdet, layer, false, isStereo, isPlusSide(detId))] = lay; } } - // FPix - for (int disk = 1; disk <= 3; ++disk) { - setForward(PixelSubdetector::PixelEndcap, pxfDetId, disk, 0); - } - // TID - for (int disk = 1; disk <= 3; ++disk) { - setForward(StripSubdetector::TID, tidDetId, disk, 0); - setForward(StripSubdetector::TID, tidDetId, disk, 1); - } - // TEC - for (int disk = 1; disk <= 9; ++disk) { - setForward(StripSubdetector::TEC, tecDetId, disk, 0); - setForward(StripSubdetector::TEC, tecDetId, disk, 1); - } return dets; } From 55840b2bc377e9471103f5f8ee391d2a06ef499d Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Thu, 12 Sep 2019 12:30:22 -0700 Subject: [PATCH 19/20] Fix subdet == TOB case --- RecoTracker/MkFit/plugins/MkFitOutputConverter.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 8406853d02f9b..94e90936121f6 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -201,7 +201,7 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: // TODO: mono/stereo structure is still hardcoded for phase0/1 strip tracker dets[lnc.convertLayerNumber(subdet, layer, false, isMono, isPlusSide(detId))] = lay; - if (((subdet == StripSubdetector::TIB or StripSubdetector::TOB) and (layer == 1 or layer == 2)) or + if (((subdet == StripSubdetector::TIB or subdet == StripSubdetector::TOB) and (layer == 1 or layer == 2)) or subdet == StripSubdetector::TID or subdet == StripSubdetector::TEC) { dets[lnc.convertLayerNumber(subdet, layer, false, isStereo, isPlusSide(detId))] = lay; } From d4140cf587de6f358408ababc01fbf9036d0d9a0 Mon Sep 17 00:00:00 2001 From: Matti Kortelainen Date: Fri, 13 Sep 2019 16:55:17 +0200 Subject: [PATCH 20/20] Throw an exception for an invalid mkFit layer index --- .../MkFit/plugins/MkFitOutputConverter.cc | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc index 94e90936121f6..87ac3e4023677 100644 --- a/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -186,8 +186,17 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: auto isPlusSide = [&ttopo](const DetId& detid) { return ttopo.side(detid) == static_cast(TrackerDetSide::PosEndcap); }; - constexpr int isMono = 0; - constexpr int isStereo = 1; + auto setDet = [&lnc, &dets, &isPlusSide]( + const int subdet, const int layer, const int isStereo, const DetId& detId, const DetLayer* lay) { + const int index = lnc.convertLayerNumber(subdet, layer, false, isStereo, isPlusSide(detId)); + if (index < 0 or static_cast(index) >= dets.size()) { + throw cms::Exception("LogicError") << "Invalid mkFit layer index " << index << " for DetId " << detId + << " subdet " << subdet << " layer " << layer << " isStereo " << isStereo; + } + dets[index] = lay; + }; + constexpr int monoLayer = 0; + constexpr int stereoLayer = 1; for (const DetLayer* lay : tracker.allLayers()) { const auto& comp = lay->basicComponents(); if (UNLIKELY(comp.empty())) { @@ -200,10 +209,10 @@ std::vector MkFitOutputConverter::createDetLayers(const mkfit:: const auto layer = ttopo.layer(detId); // TODO: mono/stereo structure is still hardcoded for phase0/1 strip tracker - dets[lnc.convertLayerNumber(subdet, layer, false, isMono, isPlusSide(detId))] = lay; + setDet(subdet, layer, monoLayer, detId, lay); if (((subdet == StripSubdetector::TIB or subdet == StripSubdetector::TOB) and (layer == 1 or layer == 2)) or subdet == StripSubdetector::TID or subdet == StripSubdetector::TEC) { - dets[lnc.convertLayerNumber(subdet, layer, false, isStereo, isPlusSide(detId))] = lay; + setDet(subdet, layer, stereoLayer, detId, lay); } }