diff --git a/RecoTracker/IterativeTracking/python/InitialStepPreSplitting_cff.py b/RecoTracker/IterativeTracking/python/InitialStepPreSplitting_cff.py index 4397d57b1a4eb..57990c085216e 100644 --- a/RecoTracker/IterativeTracking/python/InitialStepPreSplitting_cff.py +++ b/RecoTracker/IterativeTracking/python/InitialStepPreSplitting_cff.py @@ -156,7 +156,8 @@ import RecoTracker.MkFit.mkFitProducer_cfi as _mkFitProducer_cfi import RecoTracker.MkFit.mkFitOutputConverter_cfi as _mkFitOutputConverter_cfi mkFitSiPixelHitsPreSplitting = _mkFitSiPixelHitConverter_cfi.mkFitSiPixelHitConverter.clone( - hits = 'siPixelRecHitsPreSplitting' + hits = 'siPixelRecHitsPreSplitting', + clusters = 'siPixelClustersPreSplitting' ) mkFitSiStripHits = _mkFitSiStripHitConverter_cfi.mkFitSiStripHitConverter.clone() # TODO: figure out better place for this module? mkFitEventOfHitsPreSplitting = _mkFitEventOfHitsProducer_cfi.mkFitEventOfHitsProducer.clone( diff --git a/RecoTracker/MkFit/plugins/MkFitEventOfHitsProducer.cc b/RecoTracker/MkFit/plugins/MkFitEventOfHitsProducer.cc index 8f87957922585..9b0e5f79b9b81 100644 --- a/RecoTracker/MkFit/plugins/MkFitEventOfHitsProducer.cc +++ b/RecoTracker/MkFit/plugins/MkFitEventOfHitsProducer.cc @@ -18,7 +18,6 @@ #include "Geometry/CommonTopologies/interface/StripTopology.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -#include "RecoTracker/MkFit/interface/MkFitClusterIndexToHit.h" #include "RecoTracker/MkFit/interface/MkFitEventOfHits.h" #include "RecoTracker/MkFit/interface/MkFitGeometry.h" #include "RecoTracker/MkFit/interface/MkFitHitWrapper.h" @@ -39,15 +38,13 @@ class MkFitEventOfHitsProducer : public edm::global::EDProducer<> { private: void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; - void fill(const std::vector& hits, - mkfit::EventOfHits& eventOfHits, - const MkFitGeometry& mkFitGeom) const; + void fill(const std::vector& hits, mkfit::EventOfHits& eventOfHits, const MkFitGeometry& mkFitGeom) const; const edm::EDGetTokenT beamSpotToken_; const edm::EDGetTokenT pixelHitsToken_; const edm::EDGetTokenT stripHitsToken_; - const edm::EDGetTokenT pixelClusterIndexToHitToken_; - const edm::EDGetTokenT stripClusterIndexToHitToken_; + const edm::EDGetTokenT> pixelLayerIndexToHitToken_; + const edm::EDGetTokenT> stripLayerIndexToHitToken_; const edm::ESGetToken mkFitGeomToken_; edm::ESGetToken pixelQualityToken_; edm::ESGetToken stripQualityToken_; @@ -61,8 +58,8 @@ MkFitEventOfHitsProducer::MkFitEventOfHitsProducer(edm::ParameterSet const& iCon : beamSpotToken_{consumes(iConfig.getParameter("beamSpot"))}, pixelHitsToken_{consumes(iConfig.getParameter("pixelHits"))}, stripHitsToken_{consumes(iConfig.getParameter("stripHits"))}, - pixelClusterIndexToHitToken_{consumes(iConfig.getParameter("pixelHits"))}, - stripClusterIndexToHitToken_{consumes(iConfig.getParameter("stripHits"))}, + pixelLayerIndexToHitToken_{consumes(iConfig.getParameter("pixelHits"))}, + stripLayerIndexToHitToken_{consumes(iConfig.getParameter("stripHits"))}, mkFitGeomToken_{esConsumes()}, putToken_{produces()}, usePixelQualityDB_{iConfig.getParameter("usePixelQualityDB")}, @@ -171,8 +168,8 @@ void MkFitEventOfHitsProducer::produce(edm::StreamID iID, edm::Event& iEvent, co mkfit::StdSeq::loadDeads(*eventOfHits, deadvectors); } - fill(iEvent.get(pixelClusterIndexToHitToken_).hits(), *eventOfHits, mkFitGeom); - fill(iEvent.get(stripClusterIndexToHitToken_).hits(), *eventOfHits, mkFitGeom); + fill(iEvent.get(pixelLayerIndexToHitToken_), *eventOfHits, mkFitGeom); + fill(iEvent.get(stripLayerIndexToHitToken_), *eventOfHits, mkFitGeom); mkfit::StdSeq::cmssw_LoadHits_End(*eventOfHits); @@ -183,13 +180,12 @@ void MkFitEventOfHitsProducer::produce(edm::StreamID iID, edm::Event& iEvent, co iEvent.emplace(putToken_, std::move(eventOfHits)); } -void MkFitEventOfHitsProducer::fill(const std::vector& hits, +void MkFitEventOfHitsProducer::fill(const std::vector& layerToHits, mkfit::EventOfHits& eventOfHits, const MkFitGeometry& mkFitGeom) const { - for (unsigned int i = 0, end = hits.size(); i < end; ++i) { - const auto* hit = hits[i]; - if (hit != nullptr) { - const auto ilay = mkFitGeom.mkFitLayerNumber(hit->geographicalId()); + for (unsigned int i = 0, end = layerToHits.size(); i < end; ++i) { + const auto ilay = layerToHits[i]; + if (ilay >= 0) { eventOfHits[ilay].registerHit(i); } } diff --git a/RecoTracker/MkFit/plugins/MkFitPhase2HitConverter.cc b/RecoTracker/MkFit/plugins/MkFitPhase2HitConverter.cc index 91e50345d3848..bf5213198f274 100644 --- a/RecoTracker/MkFit/plugins/MkFitPhase2HitConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitPhase2HitConverter.cc @@ -20,11 +20,16 @@ namespace { class ConvertHitTraitsPhase2 { public: + using Clusters = Phase2TrackerCluster1DCollectionNew; + using Cluster = Clusters::data_type; + static constexpr bool applyCCC() { return false; } - static std::nullptr_t clusterCharge(const Phase2TrackerRecHit1D& hit, DetId hitId) { return nullptr; } + static float chargeScale(DetId id) { return 0; } + static const Cluster& cluster(const Clusters& prod, unsigned int index) { return prod.data()[index]; } + static std::nullptr_t clusterCharge(const Cluster&, float) { return nullptr; } static bool passCCC(std::nullptr_t) { return true; } - static void setDetails(mkfit::Hit& mhit, const Phase2TrackerCluster1D& cluster, int shortId, std::nullptr_t) { - mhit.setupAsStrip(shortId, (1 << 8) - 1, cluster.size()); + static void setDetails(mkfit::Hit& mhit, const Cluster& clu, int shortId, std::nullptr_t) { + mhit.setupAsStrip(shortId, (1 << 8) - 1, clu.size()); } }; } // namespace @@ -40,31 +45,34 @@ class MkFitPhase2HitConverter : public edm::global::EDProducer<> { void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; const edm::EDGetTokenT siPhase2RecHitToken_; + const edm::EDGetTokenT siPhase2ClusterToken_; const edm::ESGetToken ttrhBuilderToken_; const edm::ESGetToken ttopoToken_; const edm::ESGetToken mkFitGeomToken_; const edm::EDPutTokenT wrapperPutToken_; const edm::EDPutTokenT clusterIndexPutToken_; + const edm::EDPutTokenT> layerIndexPutToken_; const edm::EDPutTokenT> clusterChargePutToken_; const ConvertHitTraitsPhase2 convertTraits_; }; MkFitPhase2HitConverter::MkFitPhase2HitConverter(edm::ParameterSet const& iConfig) - : siPhase2RecHitToken_{consumes( - iConfig.getParameter("siPhase2Hits"))}, - ttrhBuilderToken_{esConsumes( - iConfig.getParameter("ttrhBuilder"))}, - ttopoToken_{esConsumes()}, - mkFitGeomToken_{esConsumes()}, - wrapperPutToken_{produces()}, - clusterIndexPutToken_{produces()}, - clusterChargePutToken_{produces>()}, + : siPhase2RecHitToken_{consumes(iConfig.getParameter("hits"))}, + siPhase2ClusterToken_{consumes(iConfig.getParameter("clusters"))}, + ttrhBuilderToken_{esConsumes(iConfig.getParameter("ttrhBuilder"))}, + ttopoToken_{esConsumes()}, + mkFitGeomToken_{esConsumes()}, + wrapperPutToken_{produces()}, + clusterIndexPutToken_{produces()}, + layerIndexPutToken_{produces()}, + clusterChargePutToken_{produces()}, convertTraits_{} {} void MkFitPhase2HitConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; - desc.add("siPhase2Hits", edm::InputTag{"siPhase2RecHits"}); + desc.add("hits", edm::InputTag{"siPhase2RecHits"}); + desc.add("clusters", edm::InputTag{"siPhase2Clusters"}); desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"}); descriptions.add("mkFitPhase2HitConverterDefault", desc); @@ -77,6 +85,7 @@ void MkFitPhase2HitConverter::produce(edm::StreamID iID, edm::Event& iEvent, con MkFitHitWrapper hitWrapper; MkFitClusterIndexToHit clusterIndexToHit; + std::vector layerIndexToHit; std::vector clusterCharge; edm::ProductID stripClusterID; @@ -85,18 +94,22 @@ void MkFitPhase2HitConverter::produce(edm::StreamID iID, edm::Event& iEvent, con if (not phase2Hits.empty()) { stripClusterID = mkfit::convertHits(ConvertHitTraitsPhase2{}, phase2Hits, + iEvent.get(siPhase2ClusterToken_), hitWrapper.hits(), clusterIndexToHit.hits(), + layerIndexToHit, dummy, ttopo, ttrhBuilder, - mkFitGeom); + mkFitGeom, + phase2Hits.dataSize()); } hitWrapper.setClustersID(stripClusterID); iEvent.emplace(wrapperPutToken_, std::move(hitWrapper)); iEvent.emplace(clusterIndexPutToken_, std::move(clusterIndexToHit)); + iEvent.emplace(layerIndexPutToken_, std::move(layerIndexToHit)); iEvent.emplace(clusterChargePutToken_, std::move(clusterCharge)); } diff --git a/RecoTracker/MkFit/plugins/MkFitProducer.cc b/RecoTracker/MkFit/plugins/MkFitProducer.cc index 756a2c06bb93a..12a16f8b4392a 100644 --- a/RecoTracker/MkFit/plugins/MkFitProducer.cc +++ b/RecoTracker/MkFit/plugins/MkFitProducer.cc @@ -177,7 +177,7 @@ void MkFitProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::Ev stripContainerMask.copyMaskTo(stripMask); } } else { - if (mkFitGeom.isPhase1()) + if (mkFitGeom.isPhase1() && minGoodStripCharge_ > 0) stripClusterChargeCut(iEvent.get(stripClusterChargeToken_), stripMask); } diff --git a/RecoTracker/MkFit/plugins/MkFitSiPixelHitConverter.cc b/RecoTracker/MkFit/plugins/MkFitSiPixelHitConverter.cc index 39d14c66c2d1c..a89067b68466d 100644 --- a/RecoTracker/MkFit/plugins/MkFitSiPixelHitConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitSiPixelHitConverter.cc @@ -19,11 +19,16 @@ namespace { struct ConvertHitTraits { + using Clusters = SiPixelClusterCollectionNew; + using Cluster = Clusters::data_type; + static constexpr bool applyCCC() { return false; } - static std::nullptr_t clusterCharge(const SiPixelRecHit& hit, DetId hitId) { return nullptr; } + static float chargeScale(DetId) { return 0; } + static const Cluster& cluster(const Clusters& prod, unsigned int index) { return prod.data()[index]; } + static std::nullptr_t clusterCharge(const Cluster&, float) { return nullptr; } static bool passCCC(std::nullptr_t) { return true; } - static void setDetails(mkfit::Hit& mhit, const SiPixelCluster& cluster, int shortId, std::nullptr_t) { - mhit.setupAsPixel(shortId, cluster.sizeX(), cluster.sizeY()); + static void setDetails(mkfit::Hit& mhit, const Cluster& clu, int shortId, std::nullptr_t) { + mhit.setupAsPixel(shortId, clu.sizeX(), clu.sizeY()); } }; } // namespace @@ -39,26 +44,30 @@ class MkFitSiPixelHitConverter : public edm::global::EDProducer<> { void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override; const edm::EDGetTokenT pixelRecHitToken_; + const edm::EDGetTokenT pixelClusterToken_; const edm::ESGetToken ttrhBuilderToken_; const edm::ESGetToken ttopoToken_; const edm::ESGetToken mkFitGeomToken_; const edm::EDPutTokenT wrapperPutToken_; const edm::EDPutTokenT clusterIndexPutToken_; + const edm::EDPutTokenT> layerIndexPutToken_; }; MkFitSiPixelHitConverter::MkFitSiPixelHitConverter(edm::ParameterSet const& iConfig) - : pixelRecHitToken_{consumes(iConfig.getParameter("hits"))}, - ttrhBuilderToken_{esConsumes( - iConfig.getParameter("ttrhBuilder"))}, - ttopoToken_{esConsumes()}, - mkFitGeomToken_{esConsumes()}, - wrapperPutToken_{produces()}, - clusterIndexPutToken_{produces()} {} + : pixelRecHitToken_{consumes(iConfig.getParameter("hits"))}, + pixelClusterToken_{consumes(iConfig.getParameter("clusters"))}, + ttrhBuilderToken_{esConsumes(iConfig.getParameter("ttrhBuilder"))}, + ttopoToken_{esConsumes()}, + mkFitGeomToken_{esConsumes()}, + wrapperPutToken_{produces()}, + clusterIndexPutToken_{produces()}, + layerIndexPutToken_{produces()} {} void MkFitSiPixelHitConverter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add("hits", edm::InputTag{"siPixelRecHits"}); + desc.add("clusters", edm::InputTag{"siPixelClusters"}); desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"}); descriptions.addWithDefaultLabel(desc); @@ -71,21 +80,27 @@ void MkFitSiPixelHitConverter::produce(edm::StreamID iID, edm::Event& iEvent, co MkFitHitWrapper hitWrapper; MkFitClusterIndexToHit clusterIndexToHit; + std::vector layerIndexToHit; std::vector dummy; + auto const& hits = iEvent.get(pixelRecHitToken_); auto pixelClusterID = mkfit::convertHits(ConvertHitTraits{}, - iEvent.get(pixelRecHitToken_), + hits, + iEvent.get(pixelClusterToken_), hitWrapper.hits(), clusterIndexToHit.hits(), + layerIndexToHit, dummy, ttopo, ttrhBuilder, - mkFitGeom); + mkFitGeom, + hits.dataSize()); hitWrapper.setClustersID(pixelClusterID); iEvent.emplace(wrapperPutToken_, std::move(hitWrapper)); iEvent.emplace(clusterIndexPutToken_, std::move(clusterIndexToHit)); + iEvent.emplace(layerIndexPutToken_, std::move(layerIndexToHit)); } DEFINE_FWK_MODULE(MkFitSiPixelHitConverter); diff --git a/RecoTracker/MkFit/plugins/MkFitSiStripHitConverter.cc b/RecoTracker/MkFit/plugins/MkFitSiStripHitConverter.cc index d15dff6788f59..632e74a20a7f5 100644 --- a/RecoTracker/MkFit/plugins/MkFitSiStripHitConverter.cc +++ b/RecoTracker/MkFit/plugins/MkFitSiStripHitConverter.cc @@ -22,14 +22,16 @@ namespace { class ConvertHitTraits { public: ConvertHitTraits(float minCharge) : minGoodStripCharge_(minCharge) {} + using Clusters = edmNew::DetSetVector; + using Cluster = Clusters::data_type; static constexpr bool applyCCC() { return true; } - static float clusterCharge(const SiStripRecHit2D& hit, DetId hitId) { - return siStripClusterTools::chargePerCM(hitId, hit.firstClusterRef().stripCluster()); - } + static float chargeScale(DetId id) { return siStripClusterTools::sensorThicknessInverse(id); } + static const Cluster& cluster(const Clusters& prod, unsigned int index) { return prod.data()[index]; } + static float clusterCharge(const Cluster& clu, float scale) { return clu.charge() * scale; } bool passCCC(float charge) const { return charge > minGoodStripCharge_; } - static void setDetails(mkfit::Hit& mhit, const SiStripCluster& cluster, int shortId, float charge) { - mhit.setupAsStrip(shortId, charge, cluster.amplitudes().size()); + static void setDetails(mkfit::Hit& mhit, const Cluster& clu, int shortId, float charge) { + mhit.setupAsStrip(shortId, charge, clu.size()); } private: @@ -49,25 +51,28 @@ class MkFitSiStripHitConverter : public edm::global::EDProducer<> { const edm::EDGetTokenT stripRphiRecHitToken_; const edm::EDGetTokenT stripStereoRecHitToken_; + const edm::EDGetTokenT> stripClusterToken_; const edm::ESGetToken ttrhBuilderToken_; const edm::ESGetToken ttopoToken_; const edm::ESGetToken mkFitGeomToken_; const edm::EDPutTokenT wrapperPutToken_; const edm::EDPutTokenT clusterIndexPutToken_; + const edm::EDPutTokenT> layerIndexPutToken_; const edm::EDPutTokenT> clusterChargePutToken_; const ConvertHitTraits convertTraits_; }; MkFitSiStripHitConverter::MkFitSiStripHitConverter(edm::ParameterSet const& iConfig) - : stripRphiRecHitToken_{consumes(iConfig.getParameter("rphiHits"))}, - stripStereoRecHitToken_{consumes(iConfig.getParameter("stereoHits"))}, - ttrhBuilderToken_{esConsumes( - iConfig.getParameter("ttrhBuilder"))}, - ttopoToken_{esConsumes()}, - mkFitGeomToken_{esConsumes()}, - wrapperPutToken_{produces()}, - clusterIndexPutToken_{produces()}, - clusterChargePutToken_{produces>()}, + : stripRphiRecHitToken_{consumes(iConfig.getParameter("rphiHits"))}, + stripStereoRecHitToken_{consumes(iConfig.getParameter("stereoHits"))}, + stripClusterToken_{consumes(iConfig.getParameter("clusters"))}, + ttrhBuilderToken_{esConsumes(iConfig.getParameter("ttrhBuilder"))}, + ttopoToken_{esConsumes()}, + mkFitGeomToken_{esConsumes()}, + wrapperPutToken_{produces()}, + clusterIndexPutToken_{produces()}, + layerIndexPutToken_{produces()}, + clusterChargePutToken_{produces()}, convertTraits_{static_cast( iConfig.getParameter("minGoodStripCharge").getParameter("value"))} {} @@ -76,6 +81,7 @@ void MkFitSiStripHitConverter::fillDescriptions(edm::ConfigurationDescriptions& desc.add("rphiHits", edm::InputTag{"siStripMatchedRecHits", "rphiRecHit"}); desc.add("stereoHits", edm::InputTag{"siStripMatchedRecHits", "stereoRecHit"}); + desc.add("clusters", edm::InputTag{"siStripClusters"}); desc.add("ttrhBuilder", edm::ESInputTag{"", "WithTrackAngle"}); edm::ParameterSetDescription descCCC; @@ -92,16 +98,29 @@ void MkFitSiStripHitConverter::produce(edm::StreamID iID, edm::Event& iEvent, co MkFitHitWrapper hitWrapper; MkFitClusterIndexToHit clusterIndexToHit; + std::vector layerIndexToHit; std::vector clusterCharge; - auto convert = [&](auto& hits) { - return mkfit::convertHits( - convertTraits_, hits, hitWrapper.hits(), clusterIndexToHit.hits(), clusterCharge, ttopo, ttrhBuilder, mkFitGeom); - }; - edm::ProductID stripClusterID; const auto& stripRphiHits = iEvent.get(stripRphiRecHitToken_); const auto& stripStereoHits = iEvent.get(stripStereoRecHitToken_); + const auto maxSizeGuess(stripRphiHits.dataSize() + stripStereoHits.dataSize()); + auto const& clusters = iEvent.get(stripClusterToken_); + + auto convert = [&](auto& hits) { + return mkfit::convertHits(convertTraits_, + hits, + clusters, + hitWrapper.hits(), + clusterIndexToHit.hits(), + layerIndexToHit, + clusterCharge, + ttopo, + ttrhBuilder, + mkFitGeom, + maxSizeGuess); + }; + if (not stripRphiHits.empty()) { stripClusterID = convert(stripRphiHits); } @@ -119,6 +138,7 @@ void MkFitSiStripHitConverter::produce(edm::StreamID iID, edm::Event& iEvent, co iEvent.emplace(wrapperPutToken_, std::move(hitWrapper)); iEvent.emplace(clusterIndexPutToken_, std::move(clusterIndexToHit)); + iEvent.emplace(layerIndexPutToken_, std::move(layerIndexToHit)); iEvent.emplace(clusterChargePutToken_, std::move(clusterCharge)); } diff --git a/RecoTracker/MkFit/plugins/convertHits.h b/RecoTracker/MkFit/plugins/convertHits.h index 6560f50e208e6..b2702f94f17d8 100644 --- a/RecoTracker/MkFit/plugins/convertHits.h +++ b/RecoTracker/MkFit/plugins/convertHits.h @@ -21,52 +21,43 @@ #include "RecoTracker/MkFitCore/interface/HitStructures.h" namespace mkfit { - template + template edm::ProductID convertHits(const Traits& traits, const HitCollection& hits, + const ClusterCollection& clusters, mkfit::HitVec& mkFitHits, std::vector& clusterIndexToHit, + std::vector& layerIndexToHit, std::vector& clusterChargeVec, const TrackerTopology& ttopo, const TransientTrackingRecHitBuilder& ttrhBuilder, - const MkFitGeometry& mkFitGeom) { + const MkFitGeometry& mkFitGeom, + std::size_t maxSizeGuess = 0) { if (hits.empty()) return edm::ProductID{}; - edm::ProductID clusterID; - { - const auto& lastClusterRef = hits.data().back().firstClusterRef(); - clusterID = lastClusterRef.id(); - if (lastClusterRef.index() >= mkFitHits.size()) { - auto const size = lastClusterRef.index(); - mkFitHits.resize(size); - clusterIndexToHit.resize(size, nullptr); - if constexpr (Traits::applyCCC()) { - clusterChargeVec.resize(size, -1.f); - } + const auto& lastClusterRef = hits.data().back().firstClusterRef(); + edm::ProductID clusterID = lastClusterRef.id(); + auto const size = std::max(static_cast(lastClusterRef.index() + 1), maxSizeGuess); + if (mkFitHits.size() < size) { + mkFitHits.resize(size); + clusterIndexToHit.resize(size, nullptr); + layerIndexToHit.resize(size, -1); + if constexpr (Traits::applyCCC()) { + clusterChargeVec.resize(size, -1.f); } } for (const auto& detset : hits) { + if (detset.empty()) + continue; const DetId detid = detset.detId(); const auto ilay = mkFitGeom.mkFitLayerNumber(detid); + const auto uniqueIdInLayer = mkFitGeom.uniqueIdInLayer(ilay, detid.rawId()); + const auto chargeScale = traits.chargeScale(detid); + const auto& surf = detset.begin()->det()->surface(); for (const auto& hit : detset) { - const auto charge = traits.clusterCharge(hit, detid); - if (!traits.passCCC(charge)) - 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(); - auto clusterRef = hit.firstClusterRef(); if UNLIKELY (clusterRef.id() != clusterID) { throw cms::Exception("LogicError") @@ -74,6 +65,22 @@ namespace mkfit { << ", but encountered Ref to product " << clusterRef.id() << " on detid " << detid.rawId(); } const auto clusterIndex = clusterRef.index(); + + const auto& clu = traits.cluster(clusters, clusterIndex); + const auto charge = traits.clusterCharge(clu, chargeScale); + if (!traits.passCCC(charge)) + continue; + + const auto& gpos = surf.toGlobal(hit.localPosition()); + SVector3 pos(gpos.x(), gpos.y(), gpos.z()); + const auto& gerr = ErrorFrameTransformer::transform(hit.localPositionError(), surf); + SMatrixSym33 err{{float(gerr.cxx()), + float(gerr.cyx()), + float(gerr.cyy()), + float(gerr.czx()), + float(gerr.czy()), + float(gerr.czz())}}; + LogTrace("MkFitHitConverter") << "Adding hit detid " << detid.rawId() << " subdet " << detid.subdetId() << " layer " << ttopo.layer(detid) << " isStereo " << ttopo.isStereo(detid) << " zplus " @@ -82,18 +89,19 @@ namespace mkfit { if UNLIKELY (clusterIndex >= mkFitHits.size()) { mkFitHits.resize(clusterIndex + 1); clusterIndexToHit.resize(clusterIndex + 1, nullptr); + layerIndexToHit.resize(clusterIndex + 1, -1); if constexpr (Traits::applyCCC()) { clusterChargeVec.resize(clusterIndex + 1, -1.f); } } mkFitHits[clusterIndex] = mkfit::Hit(pos, err); clusterIndexToHit[clusterIndex] = &hit; + layerIndexToHit[clusterIndex] = ilay; if constexpr (Traits::applyCCC()) { clusterChargeVec[clusterIndex] = charge; } - const auto uniqueIdInLayer = mkFitGeom.uniqueIdInLayer(ilay, detid.rawId()); - traits.setDetails(mkFitHits[clusterIndex], *(hit.cluster()), uniqueIdInLayer, charge); + traits.setDetails(mkFitHits[clusterIndex], clu, uniqueIdInLayer, charge); } } return clusterID; diff --git a/RecoTracker/MkFit/python/customizeHLTIter0ToMkFit.py b/RecoTracker/MkFit/python/customizeHLTIter0ToMkFit.py index 8a123556109c7..ac3d7b13fe663 100644 --- a/RecoTracker/MkFit/python/customizeHLTIter0ToMkFit.py +++ b/RecoTracker/MkFit/python/customizeHLTIter0ToMkFit.py @@ -24,7 +24,40 @@ def customizeHLTIter0ToMkFit(process): return process # mkFit needs all clusters, so switch off the on-demand mode + process.hltSiStripRawToClustersFacility = cms.EDProducer( + "SiStripClusterizerFromRaw", + ProductLabel = cms.InputTag( "rawDataCollector" ), + ConditionsLabel = cms.string( "" ), + onDemand = cms.bool( True ), + DoAPVEmulatorCheck = cms.bool( False ), + LegacyUnpacker = cms.bool( False ), + HybridZeroSuppressed = cms.bool( False ), + Clusterizer = cms.PSet( + ConditionsLabel = cms.string( "" ), + MaxClusterSize = cms.uint32( 32 ), + ClusterThreshold = cms.double( 5.0 ), + SeedThreshold = cms.double( 3.0 ), + Algorithm = cms.string( "ThreeThresholdAlgorithm" ), + ChannelThreshold = cms.double( 2.0 ), + MaxAdjacentBad = cms.uint32( 0 ), + setDetId = cms.bool( True ), + MaxSequentialHoles = cms.uint32( 0 ), + RemoveApvShots = cms.bool( True ), + clusterChargeCut = cms.PSet( refToPSet_ = cms.string( "HLTSiStripClusterChargeCutNone" ) ), + MaxSequentialBad = cms.uint32( 1 ) + ), + Algorithms = cms.PSet( + Use10bitsTruncation = cms.bool( False ), + CommonModeNoiseSubtractionMode = cms.string( "Median" ), + useCMMeanMap = cms.bool( False ), + TruncateInSuppressor = cms.bool( True ), + doAPVRestore = cms.bool( False ), + SiStripFedZeroSuppressionMode = cms.uint32( 4 ), + PedestalSubtractionFedMode = cms.bool( True ) + ) + ) process.hltSiStripRawToClustersFacility.onDemand = False + process.hltSiStripRawToClustersFacility.Clusterizer.MaxClusterSize = 8 process.hltSiStripRecHits = SiStripRecHitConverter_cfi.siStripMatchedRecHits.clone( ClusterProducer = "hltSiStripRawToClustersFacility", @@ -35,15 +68,17 @@ def customizeHLTIter0ToMkFit(process): # Use fourth hit if one is available process.hltIter0PFLowPixelSeedsFromPixelTracks.includeFourthHit = cms.bool(True) - process.hltMkFitGeometryESProducer = mkFitGeometryESProducer_cfi.mkFitGeometryESProducer.clone() + process.load("RecoTracker.MkFit.mkFitGeometryESProducer_cfi") process.hltIter0PFlowCkfTrackCandidatesMkFitSiPixelHits = mkFitSiPixelHitConverter_cfi.mkFitSiPixelHitConverter.clone( hits = "hltSiPixelRecHits", + clusters = "hltSiPixelClusters", ttrhBuilder = ":hltESPTTRHBWithTrackAngle", ) process.hltIter0PFlowCkfTrackCandidatesMkFitSiStripHits = mkFitSiStripHitConverter_cfi.mkFitSiStripHitConverter.clone( rphiHits = "hltSiStripRecHits:rphiRecHit", stereoHits = "hltSiStripRecHits:stereoRecHit", + clusters = "hltSiStripRawToClustersFacility", ttrhBuilder = ":hltESPTTRHBWithTrackAngle", minGoodStripCharge = dict(refToPSet_ = 'HLTSiStripClusterChargeCutLoose'), ) @@ -66,7 +101,7 @@ def customizeHLTIter0ToMkFit(process): eventOfHits = "hltIter0PFlowCkfTrackCandidatesMkFitEventOfHits", seeds = "hltIter0PFlowCkfTrackCandidatesMkFitSeeds", config = ('', 'hltIter0PFlowTrackCandidatesMkFitConfig'), - minGoodStripCharge = dict(refToPSet_ = 'HLTSiStripClusterChargeCutLoose'), + minGoodStripCharge = dict(refToPSet_ = 'HLTSiStripClusterChargeCutNone'), ) process.hltIter0PFlowCkfTrackCandidates = mkFitOutputConverter_cfi.mkFitOutputConverter.clone( seeds = "hltIter0PFLowPixelSeedsFromPixelTracks", @@ -95,4 +130,136 @@ def customizeHLTIter0ToMkFit(process): if not path.contains(process.HLTIterativeTrackingIteration0) and path.contains(process.hltIter0PFlowCkfTrackCandidates): path.replace(process.hltIter0PFlowCkfTrackCandidates, replaceWith) + process.hltIter0PFlowTrackCandidatesMkFitConfig.config = 'RecoTracker/MkFit/data/mkfit-phase1-hltiter0.json' + + process.hltIter0PFlowTrackCutClassifier.mva.maxChi2 = cms.vdouble( 999.0, 25.0, 99.0 ) + + process.hltIter0PFlowTrackCutClassifier.mva.maxChi2n = cms.vdouble( 1.2, 1.0, 999.0 ) + + process.hltIter0PFlowTrackCutClassifier.mva.dr_par = cms.PSet( + d0err = cms.vdouble( 0.003, 0.003, 0.003 ), + dr_par1 = cms.vdouble( 3.40282346639E38, 0.6, 0.6 ), + dr_par2 = cms.vdouble( 3.40282346639E38, 0.45, 0.45 ), + dr_exp = cms.vint32( 4, 4, 4 ), + d0err_par = cms.vdouble( 0.001, 0.001, 0.001 ) + ) + process.hltIter0PFlowTrackCutClassifier.mva.dz_par = cms.PSet( + dz_par1 = cms.vdouble( 3.40282346639E38, 0.6, 0.6 ), + dz_par2 = cms.vdouble( 3.40282346639E38, 0.51, 0.51 ), + dz_exp = cms.vint32( 4, 4, 4 ) + ) + + if hasattr(process, 'hltIter0PFlowTrackCutClassifierSerialSync'): + process.hltIter0PFlowTrackCutClassifierSerialSync.mva.maxChi2 = cms.vdouble( 999.0, 25.0, 99.0 ) + process.hltIter0PFlowTrackCutClassifierSerialSync.mva.maxChi2n = cms.vdouble( 1.2, 1.0, 999.0 ) + process.hltIter0PFlowTrackCutClassifierSerialSync.mva.dr_par = cms.PSet( + d0err = cms.vdouble( 0.003, 0.003, 0.003 ), + dr_par1 = cms.vdouble( 3.40282346639E38, 0.6, 0.6 ), + dr_par2 = cms.vdouble( 3.40282346639E38, 0.45, 0.45 ), + dr_exp = cms.vint32( 4, 4, 4 ), + d0err_par = cms.vdouble( 0.001, 0.001, 0.001 ) + ) + process.hltIter0PFlowTrackCutClassifierSerialSync.mva.dz_par = cms.PSet( + dz_par1 = cms.vdouble( 3.40282346639E38, 0.6, 0.6 ), + dz_par2 = cms.vdouble( 3.40282346639E38, 0.51, 0.51 ), + dz_exp = cms.vint32( 4, 4, 4 ) + ) + + return process + +def customizeHLTSiStripClusterizerOnDemandFalse(process): + + # if any of the following objects does not exist, do not apply any customisation + for objLabel in [ + 'hltSiStripRawToClustersFacility', + ]: + if not hasattr(process, objLabel): + print(f'# WARNING: customize command failed (object with label "{objLabel}" not found) - no customisation applied !') + return process + + # mkFit needs all clusters, so switch off the on-demand mode + process.hltSiStripRawToClustersFacility.onDemand = False + return process + +def customizeHLTSiStripClusterizerOnDemandFalseMaxClusterSize8(process): + + for objLabel in [ + 'hltSiStripRawToClustersFacility', + ]: + if not hasattr(process, objLabel): + print(f'# WARNING: customize command failed (object with label "{objLabel}" not found) - no customisation applied !') + return process + + process.hltSiStripRawToClustersFacility = cms.EDProducer( + "SiStripClusterizerFromRaw", + ProductLabel = cms.InputTag( "rawDataCollector" ), + ConditionsLabel = cms.string( "" ), + onDemand = cms.bool( True ), + DoAPVEmulatorCheck = cms.bool( False ), + LegacyUnpacker = cms.bool( False ), + HybridZeroSuppressed = cms.bool( False ), + Clusterizer = cms.PSet( + ConditionsLabel = cms.string( "" ), + MaxClusterSize = cms.uint32( 32 ), + ClusterThreshold = cms.double( 5.0 ), + SeedThreshold = cms.double( 3.0 ), + Algorithm = cms.string( "ThreeThresholdAlgorithm" ), + ChannelThreshold = cms.double( 2.0 ), + MaxAdjacentBad = cms.uint32( 0 ), + setDetId = cms.bool( True ), + MaxSequentialHoles = cms.uint32( 0 ), + RemoveApvShots = cms.bool( True ), + clusterChargeCut = cms.PSet( refToPSet_ = cms.string( "HLTSiStripClusterChargeCutNone" ) ), + MaxSequentialBad = cms.uint32( 1 ) + ), + Algorithms = cms.PSet( + Use10bitsTruncation = cms.bool( False ), + CommonModeNoiseSubtractionMode = cms.string( "Median" ), + useCMMeanMap = cms.bool( False ), + TruncateInSuppressor = cms.bool( True ), + doAPVRestore = cms.bool( False ), + SiStripFedZeroSuppressionMode = cms.uint32( 4 ), + PedestalSubtractionFedMode = cms.bool( True ) + ) + ) + process.hltSiStripRawToClustersFacility.onDemand = False + process.hltSiStripRawToClustersFacility.Clusterizer.MaxClusterSize = 8 + + return process + +def modifyMinOutputModuleForTrackingValidation(process, filename="output.root"): + + for objLabel in [ + 'hltOutputMinimal', + ]: + if not hasattr(process, objLabel): + print(f'# WARNING: customize command failed (object with label "{objLabel}" not found) - no customisation applied !') + return process + + process.hltOutputMinimal.outputCommands = cms.untracked.vstring( + 'drop *', + 'keep edmTriggerResults_*_*_*', + 'keep triggerTriggerEvent_*_*_*', + 'keep GlobalAlgBlkBXVector_*_*_*', + 'keep GlobalExtBlkBXVector_*_*_*', + 'keep l1tEGammaBXVector_*_EGamma_*', + 'keep l1tEtSumBXVector_*_EtSum_*', + 'keep l1tJetBXVector_*_Jet_*', + 'keep l1tMuonBXVector_*_Muon_*', + 'keep l1tTauBXVector_*_Tau_*', + 'keep *_*_*_HLTX', + 'drop *_hltHbherecoLegacy_*_*', + 'drop *_hlt*Pixel*SoA*_*_*', + 'keep recoGenParticles_genParticles_*_*', + 'keep TrackingParticles_*_*_*', + 'keep *_*_MergedTrackTruth_*', + 'keep *_simSiPixelDigis_*_*', + 'keep *_simSiStripDigis_*_*', + 'keep PSimHits_g4SimHits_*_*', + 'keep SimTracks_g4SimHits_*_*', + 'keep SimVertexs_g4SimHits_*_*', + 'keep PileupSummaryInfos_addPileupInfo_*_*', + ) + process.hltOutputMinimal.fileName = filename + process.schedule.remove( process.DQMOutput ) return process diff --git a/RecoTracker/MkFitCore/interface/Hit.h b/RecoTracker/MkFitCore/interface/Hit.h index efead1df7fcfd..4ace95ee5af43 100644 --- a/RecoTracker/MkFitCore/interface/Hit.h +++ b/RecoTracker/MkFitCore/interface/Hit.h @@ -5,6 +5,12 @@ #include "RecoTracker/MkFitCore/interface/MatrixSTypes.h" #include +#include +#include +#include +#include +#include +#include #include #include @@ -31,17 +37,17 @@ namespace mkfit { inline float getInvRad2(float x, float y) { return 1.0f / (x * x + y * y); } - inline float getPhi(float x, float y) { return std::atan2(y, x); } + inline float getPhi(float x, float y) { return vdt::fast_atan2f(y, x); } - inline float getTheta(float r, float z) { return std::atan2(r, z); } + inline float getTheta(float r, float z) { return vdt::fast_atan2f(r, z); } - inline float getEta(float r, float z) { return -1.0f * std::log(std::tan(getTheta(r, z) / 2.0f)); } + inline float getEta(float r, float z) { return -1.0f * vdt::fast_logf(vdt::fast_tanf(getTheta(r, z) / 2.0f)); } - inline float getEta(float theta) { return -1.0f * std::log(std::tan(theta / 2.0f)); } + inline float getEta(float theta) { return -1.0f * vdt::fast_logf(vdt::fast_tanf(theta / 2.0f)); } inline float getEta(float x, float y, float z) { - const float theta = std::atan2(std::sqrt(x * x + y * y), z); - return -1.0f * std::log(std::tan(theta / 2.0f)); + const float theta = vdt::fast_atan2f(std::sqrt(x * x + y * y), z); + return -1.0f * vdt::fast_logf(vdt::fast_tanf(theta / 2.0f)); } inline float getHypot(float x, float y) { return std::sqrt(x * x + y * y); } @@ -82,22 +88,24 @@ namespace mkfit { inline float getPxPxErr2(float ipt, float phi, float vipt, float vphi) { // ipt = 1/pT, v = variance const float iipt2 = 1.0f / (ipt * ipt); //iipt = 1/(1/pT) = pT - const float cosP = std::cos(phi); - const float sinP = std::sin(phi); + float cosP; + float sinP; + vdt::fast_sincosf(phi, sinP, cosP); return iipt2 * (iipt2 * cosP * cosP * vipt + sinP * sinP * vphi); } inline float getPyPyErr2(float ipt, float phi, float vipt, float vphi) { // ipt = 1/pT, v = variance const float iipt2 = 1.0f / (ipt * ipt); //iipt = 1/(1/pT) = pT - const float cosP = std::cos(phi); - const float sinP = std::sin(phi); + float cosP; + float sinP; + vdt::fast_sincosf(phi, sinP, cosP); return iipt2 * (iipt2 * sinP * sinP * vipt + cosP * cosP * vphi); } inline float getPzPzErr2(float ipt, float theta, float vipt, float vtheta) { // ipt = 1/pT, v = variance const float iipt2 = 1.0f / (ipt * ipt); //iipt = 1/(1/pT) = pT - const float cotT = 1.0f / std::tan(theta); - const float cscT = 1.0f / std::sin(theta); + const float cotT = 1.0f / vdt::fast_tanf(theta); + const float cscT = 1.0f / vdt::fast_sinf(theta); return iipt2 * (iipt2 * cotT * cotT * vipt + cscT * cscT * cscT * cscT * vtheta); } @@ -121,20 +129,13 @@ namespace mkfit { struct MeasurementState { public: MeasurementState() {} - MeasurementState(const SVector3& p, const SVector6& e) : pos_(p), err_(e) {} - MeasurementState(const SVector3& p, const SMatrixSym33& e) : pos_(p) { - for (int i = 0; i < 6; ++i) - err_[i] = e.Array()[i]; - } + MeasurementState(const SVector3& p, const SMatrixSym33& e) : pos_(p), err_(e) {} + const SVector3& parameters() const { return pos_; } - SMatrixSym33 errors() const { - SMatrixSym33 result; - for (int i = 0; i < 6; ++i) - result.Array()[i] = err_[i]; - return result; - } + const SMatrixSym33& errors() const { return err_; } + SVector3 pos_; - SVector6 err_; + SMatrixSym33 err_; }; class Hit { @@ -146,14 +147,14 @@ namespace mkfit { const SVector3& position() const { return state_.parameters(); } const SVector3& parameters() const { return state_.parameters(); } - const SMatrixSym33 error() const { return state_.errors(); } + const SMatrixSym33& error() const { return state_.err_; } const float* posArray() const { return state_.pos_.Array(); } const float* errArray() const { return state_.err_.Array(); } // Non-const versions needed for CopyOut of Matriplex. SVector3& parameters_nc() { return state_.pos_; } - SVector6& error_nc() { return state_.err_; } + SMatrixSym33& error_nc() { return state_.err_; } float r() const { return sqrtf(state_.parameters().At(0) * state_.parameters().At(0) + diff --git a/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.cc b/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.cc index bd4223b2da009..13e0a6f4669be 100644 --- a/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.cc +++ b/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.cc @@ -934,7 +934,7 @@ namespace mkfit { MPlexQF msRad; #pragma omp simd for (int n = 0; n < NN; ++n) { - msRad.At(n, 0, 0) = std::hypot(msPar.constAt(n, 0, 0), msPar.constAt(n, 1, 0)); + msRad.At(n, 0, 0) = hipo(msPar.constAt(n, 0, 0), msPar.constAt(n, 1, 0)); } propagateHelixToRMPlex(psErr, psPar, Chg, msRad, propErr, propPar, outFailFlag, N_proc, propFlags); @@ -983,7 +983,7 @@ namespace mkfit { #pragma omp simd for (int n = 0; n < NN; ++n) { if (n < N_proc) { - msRad.At(n, 0, 0) = std::hypot(msPar.constAt(n, 0, 0), msPar.constAt(n, 1, 0)); + msRad.At(n, 0, 0) = hipo(msPar.constAt(n, 0, 0), msPar.constAt(n, 1, 0)); } else { msRad.At(n, 0, 0) = 0.0f; } @@ -1053,7 +1053,7 @@ namespace mkfit { MPlexQF rotT01; for (int n = 0; n < NN; ++n) { if (n < N_proc) { - const float r = std::hypot(msPar.constAt(n, 0, 0), msPar.constAt(n, 1, 0)); + const float r = hipo(msPar.constAt(n, 0, 0), msPar.constAt(n, 1, 0)); rotT00.At(n, 0, 0) = -(msPar.constAt(n, 1, 0) + psPar.constAt(n, 1, 0)) / (2 * r); rotT01.At(n, 0, 0) = (msPar.constAt(n, 0, 0) + psPar.constAt(n, 0, 0)) / (2 * r); } else { diff --git a/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.icc b/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.icc index b6d5b73469a87..8134af5bce816 100644 --- a/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.icc +++ b/RecoTracker/MkFitCore/src/KalmanUtilsMPlex.icc @@ -54,7 +54,7 @@ template static inline void ConvertToCCS_imp(const TfLV1& a, TfLV2& b, TfLL& c, const int nmin, const int nmax) { #pragma omp simd for (int n = nmin; n < nmax; ++n) { - const float pt = getHypot(a(n, 0, 3), a(n, 0, 4)); + const float pt = hipo(a(n, 0, 3), a(n, 0, 4)); const float p2 = pt * pt + a(n, 0, 5) * a(n, 0, 5); // b(n, 0, 0) = a(n, 0, 0); diff --git a/RecoTracker/MkFitCore/src/MiniPropagators.cc b/RecoTracker/MkFitCore/src/MiniPropagators.cc index 66396bfbaaa0e..0c03258df6ca1 100644 --- a/RecoTracker/MkFitCore/src/MiniPropagators.cc +++ b/RecoTracker/MkFitCore/src/MiniPropagators.cc @@ -39,7 +39,7 @@ namespace mkfit::mini_propagators { for (int i = 0; i < Config::Niter; ++i) { // compute tangental and ideal distance for the current iteration. // 3-rd order asin for symmetric incidence (shortest arc lenght). - float r0 = std::hypot(c.x, c.y); + float r0 = hipo(c.x, c.y); float td = (R - r0) * curv; float id = oo_curv * td * (1.0f + 0.16666666f * td * td); // This would be for line approximation: @@ -67,7 +67,7 @@ namespace mkfit::mini_propagators { } } // should have some epsilon constant / member? relative vs. abs? - c.fail_flag = std::abs(std::hypot(c.x, c.y) - R) < 0.1f ? 0 : 1; + c.fail_flag = std::abs(hipo(c.x, c.y) - R) < 0.1f ? 0 : 1; return c.fail_flag; } diff --git a/RecoTracker/MkFitCore/src/MkBase.h b/RecoTracker/MkFitCore/src/MkBase.h index 2244ea1f8cccb..569c06d43a113 100644 --- a/RecoTracker/MkFitCore/src/MkBase.h +++ b/RecoTracker/MkFitCore/src/MkBase.h @@ -46,7 +46,7 @@ namespace mkfit { MPlexQF msRad; #pragma omp simd for (int n = 0; n < NN; ++n) { - msRad.At(n, 0, 0) = std::hypot(par.constAt(n, 0, 0), par.constAt(n, 1, 0)); + msRad.At(n, 0, 0) = hipo(par.constAt(n, 0, 0), par.constAt(n, 1, 0)); } propagateHelixToRMPlex( @@ -84,9 +84,9 @@ namespace mkfit { #pragma omp simd for (int n = 0; n < NN; ++n) { const float slope = std::tan(m_Par[iC].constAt(n, 5, 0)); - // msZ.At(n, 0, 0) = ( Config::beamspotz0 + slope * ( Config::beamspotr0 - std::hypot(m_Par[iC].constAt(n, 0, 0), m_Par[iC].constAt(n, 1, 0))) + slope * slope * m_Par[iC].constAt(n, 2, 0) ) / ( 1+slope*slope); // PCA w.r.t. z0, r0 + // msZ.At(n, 0, 0) = ( Config::beamspotz0 + slope * ( Config::beamspotr0 - hipo(m_Par[iC].constAt(n, 0, 0), m_Par[iC].constAt(n, 1, 0))) + slope * slope * m_Par[iC].constAt(n, 2, 0) ) / ( 1+slope*slope); // PCA w.r.t. z0, r0 msZ.At(n, 0, 0) = (slope * (slope * m_Par[iC].constAt(n, 2, 0) - - std::hypot(m_Par[iC].constAt(n, 0, 0), m_Par[iC].constAt(n, 1, 0)))) / + hipo(m_Par[iC].constAt(n, 0, 0), m_Par[iC].constAt(n, 1, 0)))) / (1 + slope * slope); // PCA to origin } diff --git a/RecoTracker/MkFitCore/src/MkBuilder.cc b/RecoTracker/MkFitCore/src/MkBuilder.cc index 71c7b717f368e..de1c9c6501321 100644 --- a/RecoTracker/MkFitCore/src/MkBuilder.cc +++ b/RecoTracker/MkFitCore/src/MkBuilder.cc @@ -120,14 +120,14 @@ namespace { const float pt = 1.f / fir->getPar(0, 0, 3); std::cout << "propagate to lay=" << ilay << " start from x=" << fir->getPar(0, 0, 0) << " y=" << fir->getPar(0, 0, 1) << " z=" << fir->getPar(0, 0, 2) - << " r=" << getHypot(fir->getPar(0, 0, 0), fir->getPar(0, 0, 1)) + << " r=" << hipo(fir->getPar(0, 0, 0), fir->getPar(0, 0, 1)) << " px=" << pt * std::cos(fir->getPar(0, 0, 4)) << " py=" << pt * std::sin(fir->getPar(0, 0, 4)) << " pz=" << pt / std::tan(fir->getPar(0, 0, 5)) << " pT=" << pt << std::endl; } void post_prop_print(int ilay, MkBase *fir) { std::cout << "propagate to lay=" << ilay << " arrive at x=" << fir->getPar(0, 1, 0) << " y=" << fir->getPar(0, 1, 1) - << " z=" << fir->getPar(0, 1, 2) << " r=" << getHypot(fir->getPar(0, 1, 0), fir->getPar(0, 1, 1)) + << " z=" << fir->getPar(0, 1, 2) << " r=" << hipo(fir->getPar(0, 1, 0), fir->getPar(0, 1, 1)) << std::endl; } diff --git a/RecoTracker/MkFitCore/src/MkFinder.cc b/RecoTracker/MkFitCore/src/MkFinder.cc index 07f42007ff373..ccc60233194dc 100644 --- a/RecoTracker/MkFitCore/src/MkFinder.cc +++ b/RecoTracker/MkFitCore/src/MkFinder.cc @@ -20,7 +20,9 @@ #include "RecoTracker/MkFitCore/standalone/RntDumper/MkFinder_selectHitIndices.icc" #endif -#include "vdt/atan2.h" +#include +#include +#include #include #include @@ -407,15 +409,15 @@ namespace mkfit { const float z = m_Par[iI].constAt(itrack, 2, 0); const float dz = std::abs(nSigmaZ * std::sqrt(m_Err[iI].constAt(itrack, 2, 2))); - const float edgeCorr = - std::abs(0.5f * (L.layer_info().rout() - L.layer_info().rin()) / std::tan(m_Par[iI].constAt(itrack, 5, 0))); + const float edgeCorr = std::abs(0.5f * (L.layer_info().rout() - L.layer_info().rin()) / + vdt::fast_tanf(m_Par[iI].constAt(itrack, 5, 0))); // XXX-NUM-ERR above, m_Err(2,2) gets negative! m_XWsrResult[itrack] = L.is_within_z_sensitive_region(z, std::sqrt(dz * dz + edgeCorr * edgeCorr)); assignbins(itrack, z, dz, phi, dphi, min_dq, max_dq, min_dphi, max_dphi); // Relax propagation-fail detection to be in line with pre-43145. - if (m_FailFlag[itrack] && std::sqrt(r2) >= L.layer_info().rin()) { + if (m_FailFlag[itrack] && r2 >= sqr(L.layer_info().rin())) { m_FailFlag[itrack] = 0; } } @@ -446,11 +448,10 @@ namespace mkfit { const float r2Inv = 1.f / r2; const float dphidx = -y * r2Inv, dphidy = x * r2Inv; const float phi = getPhi(x, y); - const float dphi2 = - calcdphi2(itrack, dphidx, dphidy) - //range from finite layer thickness - + std::pow(layerD * std::tan(m_Par[iI].At(itrack, 5, 0)) * std::sin(m_Par[iI].At(itrack, 4, 0) - phi), 2) * - r2Inv; + const float tanT = vdt::fast_tanf(m_Par[iI].At(itrack, 5, 0)); + const float dphi2 = calcdphi2(itrack, dphidx, dphidy) + //range from finite layer thickness + + std::pow(layerD * tanT * vdt::fast_sinf(m_Par[iI].At(itrack, 4, 0) - phi), 2) * r2Inv; #ifdef HARD_CHECK assert(dphi2 >= 0); #endif @@ -462,8 +463,7 @@ namespace mkfit { y * y * m_Err[iI].constAt(itrack, 1, 1) + 2 * x * y * m_Err[iI].constAt(itrack, 0, 1)) / r2); - const float edgeCorr = std::abs(0.5f * (L.layer_info().zmax() - L.layer_info().zmin()) * - std::tan(m_Par[iI].constAt(itrack, 5, 0))); + const float edgeCorr = std::abs(0.5f * (L.layer_info().zmax() - L.layer_info().zmin()) * tanT); m_XWsrResult[itrack] = L.is_within_r_sensitive_region(r, std::sqrt(dr * dr + edgeCorr * edgeCorr)); assignbins(itrack, r, dr, phi, dphi, min_dq, max_dq, min_dphi, max_dphi); @@ -627,8 +627,8 @@ namespace mkfit { float hx = thishit.x(); float hy = thishit.y(); float hz = thishit.z(); - float hr = std::hypot(hx, hy); - float hphi = std::atan2(hy, hx); + float hr = hipo(hx, hy); + float hphi = vdt::fast_atan2f(hy, hx); float hex = ngr( std::sqrt(thishit.exx()) ); float hey = ngr( std::sqrt(thishit.eyy()) ); float hez = ngr( std::sqrt(thishit.ezz()) ); @@ -640,7 +640,7 @@ namespace mkfit { float tx = m_Par[iI].At(itrack, 0, 0); float ty = m_Par[iI].At(itrack, 1, 0); float tz = m_Par[iI].At(itrack, 2, 0); - float tr = std::hypot(tx, ty); + float tr = hipo(tx, ty); float tphi = std::atan2(ty, tx); // float tchi2 = ngr( m_Chi2(itrack, 0, 0) ); // unused float tex = ngr( std::sqrt(m_Err[iI].At(itrack, 0, 0)) ); @@ -652,7 +652,7 @@ namespace mkfit { float tephi = ngr( std::sqrt( (ty * ty * tex * tex + tx * tx * tey * tey - 2.0f * tx * ty * m_Err[iI].At(itrack, 0, 1)) / (tr * tr * tr * tr)) ); - float ht_dxy = std::hypot(hx - tx, hy - ty); + float ht_dxy = hipo(hx - tx, hy - ty); float ht_dz = hz - tz; float ht_dphi = cdist(std::abs(hphi - tphi)); @@ -1041,7 +1041,7 @@ namespace mkfit { } } else { prop_fail = mp_is.propagate_to_z(mp::PA_Exact, L.hit_qbar(hi), mp_s, true); - new_q = std::hypot(mp_s.x, mp_s.y); + new_q = hipo(mp_s.x, mp_s.y); new_ddq = std::abs(new_q - L.hit_q(hi)); } @@ -1265,7 +1265,7 @@ namespace mkfit { } dprint("ADD FAKE HIT FOR TRACK #" << itrack << " withinBounds=" << (fake_hit_idx != Hit::kHitEdgeIdx) - << " r=" << std::hypot(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0))); + << " r=" << hipo(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0))); m_msErr.setDiagonal3x3(itrack, 666); m_msPar(itrack, 0, 0) = m_Par[iP](itrack, 0, 0); @@ -1347,16 +1347,17 @@ namespace mkfit { msErr.constAt(itrack, 1, 1) * hitT2inv}; const bool detXY_OK = std::abs(proj[0] * proj[2] - proj[1] * proj[1]) < 0.1f; //check that zero-direction is close - const float cosP = cos(pPar.constAt(itrack, 4, 0)); - const float sinP = sin(pPar.constAt(itrack, 4, 0)); - const float sinT = std::abs(sin(pPar.constAt(itrack, 5, 0))); + float sinP; + float cosP; + vdt::fast_sincosf(pPar.constAt(itrack, 4, 0), sinP, cosP); + const float sinT = std::abs(vdt::fast_sinf(pPar.constAt(itrack, 5, 0))); //qSF = sqrt[(px,py)*(1-proj)*(px,py)]/p = sinT*sqrt[(cosP,sinP)*(1-proj)*(cosP,sinP)]. qSF = detXY_OK ? sinT * std::sqrt(std::abs(1.f + cosP * cosP * proj[0] + sinP * sinP * proj[2] - 2.f * cosP * sinP * proj[1])) : 1.f; } else { //project on z // p_zLocal/p = p_z/p = cosT - qSF = std::abs(cos(pPar.constAt(itrack, 5, 0))); + qSF = std::abs(vdt::fast_cosf(pPar.constAt(itrack, 5, 0))); } const float qCorr = pcm * qSF; @@ -1623,7 +1624,7 @@ namespace mkfit { } dprint("ADD FAKE HIT FOR TRACK #" << itrack << " withinBounds=" << (fake_hit_idx != Hit::kHitEdgeIdx) - << " r=" << std::hypot(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0))); + << " r=" << hipo(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0))); // QQQ as above, only create and add if score better TrackCand newcand; @@ -2151,8 +2152,8 @@ namespace mkfit { #ifdef DEBUG_BACKWARD_FIT_BH // Dump per hit chi2 for (int i = 0; i < N_proc; ++i) { - float r_h = std::hypot(m_msPar.At(i, 0, 0), m_msPar.At(i, 1, 0)); - float r_t = std::hypot(m_Par[iC].At(i, 0, 0), m_Par[iC].At(i, 1, 0)); + float r_h = hipo(m_msPar.At(i, 0, 0), m_msPar.At(i, 1, 0)); + float r_t = hipo(m_Par[iC].At(i, 0, 0), m_Par[iC].At(i, 1, 0)); // if ((std::isnan(tmp_chi2[i]) || std::isnan(r_t))) // if ( ! std::isnan(tmp_chi2[i]) && tmp_chi2[i] > 0) // && tmp_chi2[i] > 30) @@ -2182,7 +2183,7 @@ namespace mkfit { m_Par[ti].At(i, 5, 0), // pt, phi, theta std::atan2(m_msPar.At(i, 1, 0), m_msPar.At(i, 0, 0)), // phi_h std::atan2(m_Par[ti].At(i, 1, 0), m_Par[ti].At(i, 0, 0)), // phi_t - 1e4f * std::hypot(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0), + 1e4f * hipo(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0), m_msPar.At(i, 1, 0) - m_Par[ti].At(i, 1, 0)), // d_xy 1e4f * (m_msPar.At(i, 2, 0) - m_Par[ti].At(i, 2, 0)) // d_z // e2s((m_msErr.At(i,0,0) + m_msErr.At(i,1,1)) / (r_h * r_h)), // ephi_h @@ -2415,7 +2416,7 @@ namespace mkfit { bb.pT(), beg_cur_sep, 1.0f / m_Par[ti].At(i, 3, 0), bb.posEta(), bb.posPhi(), beg_cur_sep, std::atan2(m_Par[ti].At(i, 1, 0), m_Par[ti].At(i, 0, 0)), - std::hypot(m_Par[ti].At(i, 0, 0), m_Par[ti].At(i, 1, 0)), + hipo(m_Par[ti].At(i, 0, 0), m_Par[ti].At(i, 1, 0)), m_Par[ti].At(i, 2, 0), chi_prnt, std::isnan(chi), std::isfinite(chi), chi > 0, @@ -2427,7 +2428,7 @@ namespace mkfit { e2s(std::abs(m_Err[ti].At(i, 0, 0))), e2s(std::abs(m_Err[ti].At(i, 1, 1))), e2s(std::abs(m_Err[ti].At(i, 2, 2))), // sx_t sy_t sz_t -- track errors - 1e4f * std::hypot(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0), + 1e4f * hipo(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0), m_msPar.At(i, 1, 0) - m_Par[ti].At(i, 1, 0)), // d_xy 1e4f * (m_msPar.At(i, 2, 0) - m_Par[ti].At(i, 2, 0)) // d_z ); diff --git a/RecoTracker/MkFitCore/src/MkFitter.cc b/RecoTracker/MkFitCore/src/MkFitter.cc index b345751ffe705..7efe276646fc9 100644 --- a/RecoTracker/MkFitCore/src/MkFitter.cc +++ b/RecoTracker/MkFitCore/src/MkFitter.cc @@ -22,7 +22,7 @@ namespace mkfit { void MkFitter::printPt(int idx) { for (int i = 0; i < NN; ++i) { - printf("%5.2f ", std::hypot(m_Par[idx].At(i, 3, 0), m_Par[idx].At(i, 4, 0))); + printf("%5.2f ", getHypot(m_Par[idx].At(i, 3, 0), m_Par[idx].At(i, 4, 0))); } } diff --git a/RecoTracker/MkFitCore/src/PropagationMPlex.cc b/RecoTracker/MkFitCore/src/PropagationMPlex.cc index 197e73c8b4600..dcdb9312d708e 100644 --- a/RecoTracker/MkFitCore/src/PropagationMPlex.cc +++ b/RecoTracker/MkFitCore/src/PropagationMPlex.cc @@ -1062,8 +1062,8 @@ namespace mkfit { if (fabs(sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1]) - msRad[0]) > 0.0001) { std::cout << "DID NOT GET TO R, FailFlag=" << failFlag[0] - << " dR=" << msRad[0] - std::hypot(outPar[0],outPar[1]) - << " r=" << msRad[0] << " rin=" << std::hypot(inPar[0],inPar[1]) << " rout=" << std::hypot(outPar[0],outPar[1]) + << " dR=" << msRad[0] - hipo(outPar[0],outPar[1]) + << " r=" << msRad[0] << " rin=" << hipo(inPar[0],inPar[1]) << " rout=" << hipo(outPar[0],outPar[1]) << std::endl; // std::cout << " pt=" << pt << " pz=" << inPar.At(n, 2) << std::endl; } diff --git a/RecoTracker/MkFitCore/src/PropagationMPlexCommon.cc b/RecoTracker/MkFitCore/src/PropagationMPlexCommon.cc index 767589f648ad6..98964b4773c61 100644 --- a/RecoTracker/MkFitCore/src/PropagationMPlexCommon.cc +++ b/RecoTracker/MkFitCore/src/PropagationMPlexCommon.cc @@ -5,6 +5,9 @@ #include "PropagationMPlex.h" +#include +#include + //#define DEBUG #include "Debug.h" @@ -93,26 +96,31 @@ namespace mkfit { const float ipt = outPar.constAt(n, 3, 0); const float pt = 1.f / ipt; //fixme, make sure it is positive? const float ipt2 = ipt * ipt; - const float p = pt / std::sin(theta); - const float pz = p * std::cos(theta); + float sT; + float cT; + vdt::fast_sincosf(theta, sT, cT); + const float p = pt / sT; + const float pz = p * cT; const float p2 = p * p; constexpr float mpi = 0.140; // m=140 MeV, pion constexpr float mpi2 = mpi * mpi; // m=140 MeV, pion const float beta2 = p2 / (p2 + mpi2); const float beta = std::sqrt(beta2); //radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum) - const float invCos = - p / std::abs(pt * std::cos(outPar.constAt(n, 4, 0)) * plNrm.constAt(n, 0, 0) + - pt * std::sin(outPar.constAt(n, 4, 0)) * plNrm.constAt(n, 1, 0) + pz * plNrm.constAt(n, 2, 0)); + float sinP; + float cosP; + vdt::fast_sincosf(outPar.constAt(n, 4, 0), sinP, cosP); + const float invCos = p / std::abs(pt * cosP * plNrm.constAt(n, 0, 0) + pt * sinP * plNrm.constAt(n, 1, 0) + + pz * plNrm.constAt(n, 2, 0)); radL = radL * invCos; //fixme works only for barrel geom // multiple scattering //vary independently phi and theta by the rms of the planar multiple scattering angle // XXX-KMD radL < 0, see your fixme above! Repeating bailout if (radL < 1e-13f) continue; - // const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15 + // const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*vdt::fast_logf(radL))/(beta*p);// eq 32.15 // const float thetaMSC2 = thetaMSC*thetaMSC; - const float thetaMSC = 0.0136f * (1.f + 0.038f * std::log(radL)) / (beta * p); // eq 32.15 + const float thetaMSC = 0.0136f * (1.f + 0.038f * vdt::fast_logf(radL)) / (beta * p); // eq 32.15 const float thetaMSC2 = thetaMSC * thetaMSC * radL; if /*constexpr*/ (Config::usePtMultScat) { outErr.At(n, 3, 3) += thetaMSC2 * pz * pz * ipt2 * ipt2; @@ -135,12 +143,13 @@ namespace mkfit { constexpr float me = 0.0005; // m=0.5 MeV, electron const float wmax = 2.f * me * beta2 * gamma2 / (1.f + 2.f * gamma * me / mpi + me * me / (mpi * mpi)); constexpr float I = 16.0e-9 * 10.75; - const float deltahalf = std::log(28.816e-9f * std::sqrt(2.33f * 0.498f) / I) + std::log(beta * gamma) - 0.5f; + const float deltahalf = + vdt::fast_logf(28.816e-9f * std::sqrt(2.33f * 0.498f) / I) + vdt::fast_logf(beta * gamma) - 0.5f; const float dEdx = - beta < 1.f - ? (2.f * (hitsXi.constAt(n, 0, 0) * invCos * - (0.5f * std::log(2.f * me * beta2 * gamma2 * wmax / (I * I)) - beta2 - deltahalf) / beta2)) - : 0.f; //protect against infs and nans + beta < 1.f ? (2.f * (hitsXi.constAt(n, 0, 0) * invCos * + (0.5f * vdt::fast_logf(2.f * me * beta2 * gamma2 * wmax / (I * I)) - beta2 - deltahalf) / + beta2)) + : 0.f; //protect against infs and nans // dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg //std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.constAt(n,0,0) << std::endl; const float dP = propSign.constAt(n, 0, 0) * dEdx / beta; diff --git a/RecoTracker/MkFitCore/src/PropagationMPlexEndcap.cc b/RecoTracker/MkFitCore/src/PropagationMPlexEndcap.cc index a424ffd9b6405..657f4c0f58a68 100644 --- a/RecoTracker/MkFitCore/src/PropagationMPlexEndcap.cc +++ b/RecoTracker/MkFitCore/src/PropagationMPlexEndcap.cc @@ -137,7 +137,7 @@ namespace mkfit { hitsRl(n, 0, 0) = 0.f; hitsXi(n, 0, 0) = 0.f; } else { - const auto hypo = std::hypot(outPar(n, 0, 0), outPar(n, 1, 0)); + const float hypo = hipo(outPar(n, 0, 0), outPar(n, 1, 0)); const auto mat = tinfo.material_checked(std::abs(msZ(n, 0, 0)), hypo); hitsRl(n, 0, 0) = mat.radl; hitsXi(n, 0, 0) = mat.bbxi; @@ -456,7 +456,7 @@ namespace mkfit { // Are we close to apex? Same condition as in propToR, 12.5 deg, cos(78.5deg) = 0.2 float dotp = (outPar.At(n, 0, 0) * std::cos(outPar.At(n, 4, 0)) + outPar.At(n, 1, 0) * std::sin(outPar.At(n, 4, 0))) / - std::hypot(outPar.At(n, 0, 0), outPar.At(n, 1, 0)); + hipo(outPar.At(n, 0, 0), outPar.At(n, 1, 0)); if (dotp < 0.2 || dotp < 0) { dprintf("helixAtZ: dot product bad, dotp = %f\n", dotp); outFailFlag[n] = 1; diff --git a/RecoTracker/MkFitCore/src/PropagationMPlexPlane.cc b/RecoTracker/MkFitCore/src/PropagationMPlexPlane.cc index 1861957aa21db..906c730250ed7 100644 --- a/RecoTracker/MkFitCore/src/PropagationMPlexPlane.cc +++ b/RecoTracker/MkFitCore/src/PropagationMPlexPlane.cc @@ -659,7 +659,7 @@ namespace mkfit { hitsRl(n, 0, 0) = 0.f; hitsXi(n, 0, 0) = 0.f; } else { - const float hypo = std::hypot(outPar(n, 0, 0), outPar(n, 1, 0)); + const float hypo = hipo(outPar(n, 0, 0), outPar(n, 1, 0)); const auto mat = tinfo.material_checked(std::abs(outPar(n, 2, 0)), hypo); hitsRl(n, 0, 0) = mat.radl; hitsXi(n, 0, 0) = mat.bbxi; diff --git a/RecoTracker/MkFitCore/src/Track.cc b/RecoTracker/MkFitCore/src/Track.cc index 4fdd34437b01f..1980794c21949 100644 --- a/RecoTracker/MkFitCore/src/Track.cc +++ b/RecoTracker/MkFitCore/src/Track.cc @@ -197,11 +197,11 @@ namespace mkfit { // center of helix in x,y plane const float x_center = x() - k * py(); const float y_center = y() + k * px(); - return std::hypot(x_center - x_bs, y_center - y_bs) - abs_ooc_half; + return hipo(x_center - x_bs, y_center - y_bs) - abs_ooc_half; } } float TrackBase::swimPhiToR(const float x0, const float y0) const { - const float dR = getHypot(x() - x0, y() - y0); + const float dR = hipo(x() - x0, y() - y0); // XXX-ASSUMPTION-ERROR can not always reach R, should see what callers expect. // For now return PI to signal apex on the ohter side of the helix. const float v = dR / 176.f / pT() * charge(); @@ -213,7 +213,7 @@ namespace mkfit { bool TrackBase::canReachRadius(float R) const { const float k = ((charge() < 0) ? 100.0f : -100.0f) / (Const::sol * Config::Bfield); const float ooc = 2.0f * k * pT(); - return std::abs(ooc) > R - std::hypot(x(), y()); + return std::abs(ooc) > R - hipo(x(), y()); } float TrackBase::maxReachRadius() const { @@ -222,7 +222,7 @@ namespace mkfit { // center of helix in x,y plane const float x_center = x() - k * py(); const float y_center = y() + k * px(); - return std::hypot(x_center, y_center) + abs_ooc_half; + return hipo(x_center, y_center) + abs_ooc_half; } float TrackBase::zAtR(float R, float* r_reached) const { @@ -240,14 +240,14 @@ namespace mkfit { const float lambda = pz() * ipt; //printf("Track::zAtR to R=%f: k=%e, ipt=%e, c=%e, ooc=%e -- can hit = %f (if > 1 can)\n", - // R, k, ipt, c, ooc, ooc / (R - std::hypot(xc,yc))); + // R, k, ipt, c, ooc, ooc / (R - hipo(xc,yc))); float D = 0; for (int i = 0; i < Config::Niter; ++i) { // compute tangental and ideal distance for the current iteration. // 3-rd order asin for symmetric incidence (shortest arc lenght). - float r0 = std::hypot(xc, yc); + float r0 = hipo(xc, yc); float td = (R - r0) * c; float id = ooc * td * (1.0f + 0.16666666f * td * td); // This would be for line approximation: @@ -270,7 +270,7 @@ namespace mkfit { } if (r_reached) - *r_reached = std::hypot(xc, yc); + *r_reached = hipo(xc, yc); return z() + lambda * D; @@ -325,7 +325,7 @@ namespace mkfit { // pxc = pxc * cosa - pyc * sina; // pyc = pyc * cosa + pxo * sina; - return std::hypot(xc, yc); + return hipo(xc, yc); } const char* TrackBase::algoint_to_cstr(int algo) {