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', 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..11bc66edf6800 --- /dev/null +++ b/RecoTracker/MkFit/README.md @@ -0,0 +1,30 @@ +# 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. + +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` + * 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/MkFitHitIndexMap.h b/RecoTracker/MkFit/interface/MkFitHitIndexMap.h new file mode 100644 index 0000000000000..c7b6eea50263a --- /dev/null +++ b/RecoTracker/MkFit/interface/MkFitHitIndexMap.h @@ -0,0 +1,86 @@ +#ifndef RecoTracker_MkFit_MkFitHitIndexMap_h +#define RecoTracker_MkFit_MkFitHitIndexMap_h + +#include "DataFormats/Provenance/interface/ProductID.h" + +#include + +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: + // 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} {} + + 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); + + /** + * 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); + + /// Get mkFit hit index and layer + const MkFitHit& mkFitHit(edm::ProductID id, size_t clusterIndex) const; + + /// 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 { + CMSSWHit() = default; + explicit CMSSWHit(const TrackingRecHit* p, size_t i) : ptr{p}, clusterIndex{i} {} + const TrackingRecHit* ptr = nullptr; + size_t clusterIndex = 0; + }; + + 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/interface/MkFitInputWrapper.h b/RecoTracker/MkFit/interface/MkFitInputWrapper.h new file mode 100644 index 0000000000000..ae17cf6e32584 --- /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/MkFitHitIndexMap.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(MkFitHitIndexMap&& hitIndexMap, + std::vector&& hits, + mkfit::TrackVec&& seeds, + mkfit::LayerNumberConverter&& lnc); + ~MkFitInputWrapper(); + + MkFitInputWrapper(MkFitInputWrapper const&) = delete; + MkFitInputWrapper& operator=(MkFitInputWrapper const&) = delete; + MkFitInputWrapper(MkFitInputWrapper&&); + MkFitInputWrapper& operator=(MkFitInputWrapper&&); + + 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: + MkFitHitIndexMap hitIndexMap_; + 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..8297b7c511f77 --- /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&&); + MkFitOutputWrapper& operator=(MkFitOutputWrapper&&); + + 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..7574ba951ba0c --- /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..499054817c6c8 --- /dev/null +++ b/RecoTracker/MkFit/plugins/MkFitInputConverter.cc @@ -0,0 +1,240 @@ +#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/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" +#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, + MkFitHitIndexMap& hitIndexMap, + 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 MkFitHitIndexMap& hitIndexMap, + 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_; + const float minGoodStripCharge_; +}; + +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()}, + minGoodStripCharge_{static_cast( + iConfig.getParameter("minGoodStripCharge").getParameter("value"))} {} + +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"}); + + 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 { + 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()); + 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, 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_), hitIndexMap, ttrhBuilder, iSetup.getData(mfToken_)); + + 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 { + return (siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster()) > minGoodStripCharge_); +} + +bool MkFitInputConverter::passCCC(const SiPixelRecHit& hit, const DetId hitId) const { return true; } + +template +void MkFitInputConverter::convertHits(const HitCollection& hits, + std::vector& mkFitHits, + MkFitHitIndexMap& hitIndexMap, + int& totalHits, + 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(); + hitIndexMap.resizeByClusterIndex(lastClusterRef.id(), lastClusterRef.index()); + hitIndexMap.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)); + hitIndexMap.increaseLayerSize(ilay, detset.size()); // to minimize memory allocations + + 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(); + + LogTrace("MkFitInputConverter") << "Adding hit detid " << detid.rawId() << " subdet " << subdet << " layer " + << layer << " isStereo " << isStereo << " zplus " << isPlusSide(detid) << " ilay " + << ilay; + + 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; + } + } +} + +mkfit::TrackVec MkFitInputConverter::convertSeeds(const edm::View& seeds, + const MkFitHitIndexMap& hitIndexMap, + 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) { + if (not trackerHitRTTI::isFromDet(*iHit)) { + throw cms::Exception("Assert") << "Encountered a seed with a hit which is not trackerHitRTTI::isFromDet()"; + } + const auto& clusterRef = static_cast(*iHit).firstClusterRef(); + const auto& mkFitHit = hitIndexMap.mkFitHit(clusterRef.id(), clusterRef.index()); + ret.back().addHitIdx(mkFitHit.index(), mkFitHit.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..87ac3e4023677 --- /dev/null +++ b/RecoTracker/MkFit/plugins/MkFitOutputConverter.cc @@ -0,0 +1,513 @@ +#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 "DataFormats/TrackerCommon/interface/TrackerDetSide.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 MkFitHitIndexMap& hitIndexMap, + 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 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.hitIndexMap(), + 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 isPlusSide = [&ttopo](const DetId& detid) { + return ttopo.side(detid) == static_cast(TrackerDetSide::PosEndcap); + }; + 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())) { + throw cms::Exception("LogicError") << "Got a tracker layer (subdet " << lay->subDetector() + << ") with empty basicComponents."; + } + // 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 + 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) { + setDet(subdet, layer, stereoLayer, detId, lay); + } + } + + return dets; +} + +TrackCandidateCollection MkFitOutputConverter::convertCandidates(const MkFitOutputWrapper& mkFitOutput, + const MkFitHitIndexMap& hitIndexMap, + 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; + // 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) { + // 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). + // + // 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!"; + } + // 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(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.clusterIndex( + MkFitHitIndexMap::MkFitHit{hitOnTrack.index, hitOnTrack.layer}); + 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 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; + } + + 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") + << "Candidate " << candIndex << " from seed " << seedIndex << " hit " << i + << " has different layer in candidate (" << candHitOnTrack.layer << ") and seed (" << hitOnTrack.layer + << ")." + << " Hit indices are " << candHitOnTrack.index << " and " << hitOnTrack.index << ", respectively"; + } + if (hitOnTrack.index != candHitOnTrack.index) { + throw cms::Exception("LogicError") << "Candidate " << candIndex << " from seed " << seedIndex << " hit " << i + << " has different hit index in candidate (" << candHitOnTrack.index + << ") and seed (" << hitOnTrack.index << ") on layer " << hitOnTrack.layer; + } + } + + // state + auto state = cand.state(); // copy because have to modify + state.convertFromCCSToCartesian(); + 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, // mkFit does not produce loopers, so set nLoops=0 + 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{}); + } + } + + // 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); + + // assume for now that the propagation in mkfit always alongMomentum + PropagationDirection backFitDirection = oppositeToMomentum; + + // only direction matters in this context + TrajectorySeed fakeSeed(PTrajectoryStateOnDet(), edm::OwnVector(), backFitDirection); + + // 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."; + + 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(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..33a2ac0f33853 --- /dev/null +++ b/RecoTracker/MkFit/plugins/MkFitProducer.cc @@ -0,0 +1,137 @@ +#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_); + + 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, [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 + // 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..9c4e0ef7f4e3d --- /dev/null +++ b/RecoTracker/MkFit/python/customizeInitialStepOnly.py @@ -0,0 +1,148 @@ +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/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') + ) +) diff --git a/RecoTracker/MkFit/src/MkFitHitIndexMap.cc b/RecoTracker/MkFit/src/MkFitHitIndexMap.cc new file mode 100644 index 0000000000000..89260eadbe565 --- /dev/null +++ b/RecoTracker/MkFit/src/MkFitHitIndexMap.cc @@ -0,0 +1,68 @@ +#include "RecoTracker/MkFit/interface/MkFitHitIndexMap.h" + +#include "FWCore/Utilities/interface/Exception.h" + +#include + +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->mkFitHits.size() <= clusterIndex) { + found->mkFitHits.resize(clusterIndex + 1); + } + return found; + } +} // namespace + +void MkFitHitIndexMap::resizeByClusterIndex(edm::ProductID id, size_t clusterIndex) { + resizeByClusterIndexImpl(cmsswToMkFit_, id, clusterIndex); +} + +void MkFitHitIndexMap::increaseLayerSize(int layer, size_t additionalSize) { + if (layer >= static_cast(mkFitToCMSSW_.size())) { + mkFitToCMSSW_.resize(layer + 1); + } + mkFitToCMSSW_[layer].resize(mkFitToCMSSW_[layer].size() + additionalSize); +} + +void MkFitHitIndexMap::insert(edm::ProductID id, size_t clusterIndex, MkFitHit hit, const TrackingRecHit* hitPtr) { + // mapping CMSSW->mkfit + 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 (hit.layer() >= static_cast(mkFitToCMSSW_.size())) { + mkFitToCMSSW_.resize(hit.layer() + 1); + } + auto& layer = mkFitToCMSSW_[hit.layer()]; + if (hit.index() >= static_cast(layer.size())) { + layer.resize(hit.index() + 1); + } + layer[hit.index()] = CMSSWHit(hitPtr, clusterIndex); +} + +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 : cmsswToMkFit_) { + exp << elem.productID << " "; + } + throw exp; + } + 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; +} diff --git a/RecoTracker/MkFit/src/MkFitInputWrapper.cc b/RecoTracker/MkFit/src/MkFitInputWrapper.cc new file mode 100644 index 0000000000000..6a4806efd2dd9 --- /dev/null +++ b/RecoTracker/MkFit/src/MkFitInputWrapper.cc @@ -0,0 +1,24 @@ +#include "RecoTracker/MkFit/interface/MkFitInputWrapper.h" + +// mkFit includes +#include "Hit.h" +#include "LayerNumberConverter.h" +#include "Track.h" + +MkFitInputWrapper::MkFitInputWrapper() = default; + +MkFitInputWrapper::MkFitInputWrapper(MkFitHitIndexMap&& hitIndexMap, + std::vector&& hits, + mkfit::TrackVec&& seeds, + mkfit::LayerNumberConverter&& lnc) + : hitIndexMap_{std::move(hitIndexMap)}, + hits_{std::move(hits)}, + seeds_{std::make_unique(std::move(seeds))}, + lnc_{std::make_unique(std::move(lnc))} {} + +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 new file mode 100644 index 0000000000000..0c670e22a7660 --- /dev/null +++ b/RecoTracker/MkFit/src/MkFitOutputWrapper.cc @@ -0,0 +1,14 @@ +#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; + +MkFitOutputWrapper::MkFitOutputWrapper(MkFitOutputWrapper&&) = default; +MkFitOutputWrapper& MkFitOutputWrapper::operator=(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",