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..239fcae26c385 --- /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([&geom](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",