diff --git a/CommonTools/RecoAlgos/plugins/PrimaryVertexSorter.h b/CommonTools/RecoAlgos/plugins/PrimaryVertexSorter.h index 518b4c3434807..300ce6cb2c14b 100644 --- a/CommonTools/RecoAlgos/plugins/PrimaryVertexSorter.h +++ b/CommonTools/RecoAlgos/plugins/PrimaryVertexSorter.h @@ -301,27 +301,28 @@ void PrimaryVertexSorter::produce(edm::Event& iEvent, const } template <> -void PrimaryVertexSorter>::doConsumesForTiming( +inline void PrimaryVertexSorter>::doConsumesForTiming( const edm::ParameterSet& iConfig) { tokenTrackTimeTag_ = consumes>(iConfig.getParameter("trackTimeTag")); tokenTrackTimeResoTag_ = consumes>(iConfig.getParameter("trackTimeResoTag")); } template <> -void PrimaryVertexSorter>::doConsumesForTiming(const edm::ParameterSet& iConfig) {} +inline void PrimaryVertexSorter>::doConsumesForTiming(const edm::ParameterSet& iConfig) { +} template <> -bool PrimaryVertexSorter>::needsProductsForTiming() { +inline bool PrimaryVertexSorter>::needsProductsForTiming() { return true; } template <> -bool PrimaryVertexSorter>::needsProductsForTiming() { +inline bool PrimaryVertexSorter>::needsProductsForTiming() { return false; } template <> -std::pair +inline std::pair PrimaryVertexSorter>::runAlgo(const reco::VertexCollection& vertices, const reco::RecoChargedRefCandidate& pf, const edm::ValueMap* trackTimeTag, @@ -332,7 +333,7 @@ PrimaryVertexSorter>::runAlgo(const r } template <> -std::pair PrimaryVertexSorter>::runAlgo( +inline std::pair PrimaryVertexSorter>::runAlgo( const reco::VertexCollection& vertices, const reco::PFCandidate& pf, const edm::ValueMap* trackTimeTag, diff --git a/DataFormats/DTRecHit/interface/DTDriftTimeParameters.icc b/DataFormats/DTRecHit/interface/DTDriftTimeParameters.icc index 0ca987d0d803e..029ee93a2538b 100644 --- a/DataFormats/DTRecHit/interface/DTDriftTimeParameters.icc +++ b/DataFormats/DTRecHit/interface/DTDriftTimeParameters.icc @@ -404,7 +404,8 @@ const double THIS_CLASS::fun_sigma_t[N_alpha][N_By][N_Bz][N_Sigma_t] = { /*** Auxiliary functions ***/ -unsigned short THIS_CLASS::MB_DT_Check_boundaries(double distime, double alpha, double by, double bz, short ifl) const { +inline unsigned short THIS_CLASS::MB_DT_Check_boundaries( + double distime, double alpha, double by, double bz, short ifl) const { unsigned short status = 1; std::string name = "MB_DT_drift_time"; @@ -458,7 +459,7 @@ unsigned short THIS_CLASS::MB_DT_Check_boundaries(double distime, double alpha, return (status); } -void THIS_CLASS::MB_DT_Get_grid_values( +inline void THIS_CLASS::MB_DT_Get_grid_values( double Var, unsigned short *pi, unsigned short *pj, short Initial, unsigned short N, const double *Values) const { unsigned short i, iValue, jValue; @@ -482,15 +483,15 @@ void THIS_CLASS::MB_DT_Get_grid_values( } } -void THIS_CLASS::MB_DT_Get_grid_points(double alpha, - double by, - double bz, - unsigned short *p_alpha, - unsigned short *p_By, - unsigned short *p_Bz, - unsigned short *q_alpha, - unsigned short *q_By, - unsigned short *q_Bz) const { +inline void THIS_CLASS::MB_DT_Get_grid_points(double alpha, + double by, + double bz, + unsigned short *p_alpha, + unsigned short *p_By, + unsigned short *p_Bz, + unsigned short *q_alpha, + unsigned short *q_By, + unsigned short *q_Bz) const { MB_DT_Get_grid_values(fabs(alpha), p_alpha, q_alpha, 4, N_alpha, alpha_value); MB_DT_Get_grid_values(by, p_By, q_By, 0, N_By, By_value); MB_DT_Get_grid_values(bz, p_Bz, q_Bz, 0, N_Bz, Bz_value); @@ -503,7 +504,7 @@ void THIS_CLASS::MB_DT_Get_grid_points(double alpha, /*** Multidimensional linear interpolation ***/ -double THIS_CLASS::MB_DT_MLInterpolation(double *al, double *by, double *bz, double *f) const { +inline double THIS_CLASS::MB_DT_MLInterpolation(double *al, double *by, double *bz, double *f) const { double q1, q2, q3, p1, p2, p3; double fx11, fx21, fxy1, fx12, fx22, fxy2, fxyz; @@ -528,7 +529,7 @@ double THIS_CLASS::MB_DT_MLInterpolation(double *al, double *by, double *bz, dou return (fxyz); } -double THIS_CLASS::MB_DT_sigma_t_m(double dist, double *par) const { +inline double THIS_CLASS::MB_DT_sigma_t_m(double dist, double *par) const { double x = fabs(dist); // the parametrisations are symmetric under 'distance' if (x > 20.5) @@ -537,7 +538,7 @@ double THIS_CLASS::MB_DT_sigma_t_m(double dist, double *par) const { return (par[6] * x); } -double THIS_CLASS::MB_DT_sigma_t_p(double dist, double *par) const { +inline double THIS_CLASS::MB_DT_sigma_t_p(double dist, double *par) const { double x2, x = fabs(dist); // the parametrisations are symmetric under 'distance' if (x > 20.5) diff --git a/EventFilter/HcalRawToDigi/plugins/PackerHelp.h b/EventFilter/HcalRawToDigi/plugins/PackerHelp.h index ccf166ad0395b..b9d7dcd4803eb 100644 --- a/EventFilter/HcalRawToDigi/plugins/PackerHelp.h +++ b/EventFilter/HcalRawToDigi/plugins/PackerHelp.h @@ -617,10 +617,10 @@ class UHTRpacker { // converts HE QIE digies to HB data format -QIE11DataFrame convertHB(QIE11DataFrame qiehe, - std::vector const& tdc1, - std::vector const& tdc2, - const int tdcmax) { +inline QIE11DataFrame convertHB(QIE11DataFrame qiehe, + std::vector const& tdc1, + std::vector const& tdc2, + const int tdcmax) { QIE11DataFrame qiehb = qiehe; HcalDetId did = HcalDetId(qiehb.detid()); int adc, tdc; diff --git a/MagneticField/GeomBuilder/src/buildBox.icc b/MagneticField/GeomBuilder/src/buildBox.icc index 15027724473d1..b09c4235b9b12 100644 --- a/MagneticField/GeomBuilder/src/buildBox.icc +++ b/MagneticField/GeomBuilder/src/buildBox.icc @@ -4,7 +4,7 @@ * \author N. Amapane - INFN Torino */ -void volumeHandle::buildBox(double halfX, double halfY, double halfZ) { +inline void volumeHandle::buildBox(double halfX, double halfY, double halfZ) { LogTrace("MagGeoBuilder") << "Building box surfaces...: "; // Global vectors of the normals to X, Y, Z axes diff --git a/MagneticField/GeomBuilder/src/buildCons.icc b/MagneticField/GeomBuilder/src/buildCons.icc index 3dc381fbc398e..5699ac6de8b00 100644 --- a/MagneticField/GeomBuilder/src/buildCons.icc +++ b/MagneticField/GeomBuilder/src/buildCons.icc @@ -6,13 +6,13 @@ #include "DataFormats/GeometrySurface/interface/SimpleConeBounds.h" -void volumeHandle::buildCons(double zhalf, - double rInMinusZ, - double rOutMinusZ, - double rInPlusZ, - double rOutPlusZ, - double startPhi, - double deltaPhi) { +inline void volumeHandle::buildCons(double zhalf, + double rInMinusZ, + double rOutMinusZ, + double rInPlusZ, + double rOutPlusZ, + double startPhi, + double deltaPhi) { LogTrace("MagGeoBuilder") << "Building cons surfaces...: "; LogTrace("MagGeoBuilder") << "zhalf " << zhalf << newln << "rInMinusZ " << rInMinusZ << newln << "rOutMinusZ " << rOutMinusZ << newln << "rInPlusZ " << rInPlusZ << newln << "rOutPlusZ " << rOutPlusZ diff --git a/MagneticField/GeomBuilder/src/buildPseudoTrap.icc b/MagneticField/GeomBuilder/src/buildPseudoTrap.icc index 3c57e627b090f..2f0e0c40e5c1f 100644 --- a/MagneticField/GeomBuilder/src/buildPseudoTrap.icc +++ b/MagneticField/GeomBuilder/src/buildPseudoTrap.icc @@ -12,7 +12,7 @@ * \author N. Amapane - INFN Torino */ -void volumeHandle::buildPseudoTrap( +inline void volumeHandle::buildPseudoTrap( double x1, double x2, double y1, double y2, double halfZ, double radius, bool atMinusZ) { LogTrace("MagGeoBuilder") << "Building PseudoTrap surfaces...: "; LogTrace("MagGeoBuilder") << "halfZ " << halfZ << newln << "x1 " << x1 << newln << "x2 " << x2 diff --git a/MagneticField/GeomBuilder/src/buildTrap.icc b/MagneticField/GeomBuilder/src/buildTrap.icc index 094f8c6012ab3..1b662b98c5ec5 100644 --- a/MagneticField/GeomBuilder/src/buildTrap.icc +++ b/MagneticField/GeomBuilder/src/buildTrap.icc @@ -4,17 +4,17 @@ * \author N. Amapane - INFN Torino */ -void volumeHandle::buildTrap(double x1, - double x2, - double x3, - double x4, - double y1, - double y2, - double theta, - double phi, - double halfZ, - double alpha1, - double alpha2) { +inline void volumeHandle::buildTrap(double x1, + double x2, + double x3, + double x4, + double y1, + double y2, + double theta, + double phi, + double halfZ, + double alpha1, + double alpha2) { LogTrace("MagGeoBuilder") << "Building trapezoid surfaces...: "; LogTrace("MagGeoBuilder") << "x1 " << x1 << newln << "x2 " << x2 << newln << "x3 " << x3 << newln << "x4 " << x4 << newln << "y1 " << y1 << newln << "y2 " << y2 << newln << "theta " << theta << newln diff --git a/MagneticField/GeomBuilder/src/buildTruncTubs.icc b/MagneticField/GeomBuilder/src/buildTruncTubs.icc index 068efe9e87e71..75b190578df5f 100644 --- a/MagneticField/GeomBuilder/src/buildTruncTubs.icc +++ b/MagneticField/GeomBuilder/src/buildTruncTubs.icc @@ -9,14 +9,14 @@ * \author N. Amapane - INFN Torino */ -void volumeHandle::buildTruncTubs(double zhalf, - double rIn, - double rOut, - double startPhi, - double deltaPhi, - double cutAtStart, - double cutAtDelta, - bool cutInside) { +inline void volumeHandle::buildTruncTubs(double zhalf, + double rIn, + double rOut, + double startPhi, + double deltaPhi, + double cutAtStart, + double cutAtDelta, + bool cutInside) { LogTrace("MagGeoBuilder") << "Building TruncTubs surfaces...: "; LogTrace("MagGeoBuilder") << "zhalf " << zhalf << newln << "rIn " << rIn << newln << "rOut " << rOut << newln << "startPhi " << startPhi << newln << "deltaPhi " << deltaPhi << newln diff --git a/MagneticField/GeomBuilder/src/buildTubs.icc b/MagneticField/GeomBuilder/src/buildTubs.icc index 9e9872d4e1286..3b68c996b0546 100644 --- a/MagneticField/GeomBuilder/src/buildTubs.icc +++ b/MagneticField/GeomBuilder/src/buildTubs.icc @@ -1,4 +1,4 @@ -void volumeHandle::buildTubs(double zhalf, double rIn, double rOut, double startPhi, double deltaPhi) { +inline void volumeHandle::buildTubs(double zhalf, double rIn, double rOut, double startPhi, double deltaPhi) { LogTrace("MagGeoBuilder") << "Building tubs surfaces...: "; LogTrace("MagGeoBuilder") << "zhalf " << zhalf << newln << "rIn " << rIn << newln << "rOut " << rOut << newln << "startPhi " << startPhi << newln << "deltaPhi " << deltaPhi; diff --git a/PhysicsTools/PatUtils/interface/PATJetCorrExtractor.h b/PhysicsTools/PatUtils/interface/PATJetCorrExtractor.h index e1448ab80b670..abe5d5d5f1ca1 100644 --- a/PhysicsTools/PatUtils/interface/PATJetCorrExtractor.h +++ b/PhysicsTools/PatUtils/interface/PATJetCorrExtractor.h @@ -29,25 +29,6 @@ #include #include -namespace { - std::string format_vstring(const std::vector& v) { - std::string retVal; - - retVal.append("{ "); - - unsigned numEntries = v.size(); - for (unsigned iEntry = 0; iEntry < numEntries; ++iEntry) { - retVal.append(v[iEntry]); - if (iEntry < (numEntries - 1)) - retVal.append(", "); - } - - retVal.append(" }"); - - return retVal; - } -} // namespace - class PATJetCorrExtractor { public: reco::Candidate::LorentzVector operator()( @@ -83,11 +64,20 @@ class PATJetCorrExtractor { } catch (cms::Exception const&) { throw cms::Exception("InvalidRequest") << "The JEC level " << jetCorrLabel << " does not exist !!\n" - << "Available levels = " << format_vstring(jet.availableJECLevels()) << ".\n"; + << "Available levels = { " << format_vstring(jet.availableJECLevels()) << " }.\n"; } return corrJetP4; } + +private: + static std::string format_vstring(const std::vector& v) { + std::string retVal; + auto ss = std::accumulate(v.begin(), v.end(), 0, [](int a, std::string const& s) { return a + s.length() + 2; }); + retVal.reserve(ss); + for_each(v.begin(), v.end(), [&](std::string const& s) { retVal += (retVal.empty() ? "" : ", ") + s; }); + return retVal; + } }; #endif diff --git a/RecoBTag/FeatureTools/interface/SeedingTracksConverter.h b/RecoBTag/FeatureTools/interface/SeedingTracksConverter.h index 9fed81b17c46c..d02407f6153c2 100644 --- a/RecoBTag/FeatureTools/interface/SeedingTracksConverter.h +++ b/RecoBTag/FeatureTools/interface/SeedingTracksConverter.h @@ -29,7 +29,7 @@ namespace btagbtvdeep { bool computeProbabilities, std::vector& seedingT_features_vector); - float logWithOffset(float v, float logOffset = 0) { + inline float logWithOffset(float v, float logOffset = 0) { if (v == 0.) return 0.; return logOffset + log(std::fabs(v)) * std::copysign(1.f, v); diff --git a/RecoBTag/ImpactParameter/plugins/IPProducer.h b/RecoBTag/ImpactParameter/plugins/IPProducer.h index 752ca49cdd564..56bead210899c 100644 --- a/RecoBTag/ImpactParameter/plugins/IPProducer.h +++ b/RecoBTag/ImpactParameter/plugins/IPProducer.h @@ -441,7 +441,7 @@ void IPProducer::checkEventSetup(const edm::EventSetup& // Specialized templates used to fill 'descriptions' // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ template <> -void IPProducer::fillDescriptions( +inline void IPProducer::fillDescriptions( edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add("maximumTransverseImpactParameter", 0.2); @@ -462,7 +462,8 @@ void IPProducer -void IPProducer, reco::JetTagInfo, IPProducerHelpers::FromJetAndCands>::fillDescriptions( +inline void +IPProducer, reco::JetTagInfo, IPProducerHelpers::FromJetAndCands>::fillDescriptions( edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add("maximumTransverseImpactParameter", 0.2); diff --git a/RecoBTag/SecondaryVertex/interface/CombinedSVSoftLeptonComputer.h b/RecoBTag/SecondaryVertex/interface/CombinedSVSoftLeptonComputer.h index 42bcf34520458..20ac814dc35d6 100644 --- a/RecoBTag/SecondaryVertex/interface/CombinedSVSoftLeptonComputer.h +++ b/RecoBTag/SecondaryVertex/interface/CombinedSVSoftLeptonComputer.h @@ -27,7 +27,9 @@ class CombinedSVSoftLeptonComputer : public CombinedSVComputer { bool SoftLeptonFlip; }; -double CombinedSVSoftLeptonComputer::flipSoftLeptonValue(double value) const { return SoftLeptonFlip ? -value : value; } +inline double CombinedSVSoftLeptonComputer::flipSoftLeptonValue(double value) const { + return SoftLeptonFlip ? -value : value; +} template reco::TaggingVariableList CombinedSVSoftLeptonComputer::operator()(const IPTI &ipInfo, diff --git a/RecoEgamma/EgammaTools/interface/EGExtraInfoModifierFromValueMaps.h b/RecoEgamma/EgammaTools/interface/EGExtraInfoModifierFromValueMaps.h index 3465fe253b5bb..2cd2e3aa17d4e 100644 --- a/RecoEgamma/EgammaTools/interface/EGExtraInfoModifierFromValueMaps.h +++ b/RecoEgamma/EgammaTools/interface/EGExtraInfoModifierFromValueMaps.h @@ -281,7 +281,7 @@ void EGXtraModFromVMObjFiller::addValueToObject( template <> template <> -void EGXtraModFromVMObjFiller::addValuesToObject( +inline void EGXtraModFromVMObjFiller::addValuesToObject( pat::Electron& obj, const edm::Ptr& ptr, const std::unordered_map>>& vmaps_token, @@ -299,7 +299,7 @@ void EGXtraModFromVMObjFiller::addValuesToObject( template <> template <> -void EGXtraModFromVMObjFiller::addValuesToObject( +inline void EGXtraModFromVMObjFiller::addValuesToObject( pat::Photon& obj, const edm::Ptr& ptr, const std::unordered_map>>& vmaps_token, diff --git a/RecoEgamma/EgammaTools/interface/MVAValueMapProducer.h b/RecoEgamma/EgammaTools/interface/MVAValueMapProducer.h index 6f15997541252..9d61fdd95eb10 100644 --- a/RecoEgamma/EgammaTools/interface/MVAValueMapProducer.h +++ b/RecoEgamma/EgammaTools/interface/MVAValueMapProducer.h @@ -32,6 +32,35 @@ class MVAValueMapProducer : public edm::global::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); private: + static auto getMVAEstimators(const edm::VParameterSet& vConfig) { + std::vector> mvaEstimators; + + // Loop over the list of MVA configurations passed here from python and + // construct all requested MVA estimators. + for (auto& imva : vConfig) { + // The factory below constructs the MVA of the appropriate type based + // on the "mvaName" which is the name of the derived MVA class (plugin) + if (!imva.empty()) { + mvaEstimators.emplace_back( + AnyMVAEstimatorRun2Factory::get()->create(imva.getParameter("mvaName"), imva)); + + } else + throw cms::Exception(" MVA configuration not found: ") + << " failed to find proper configuration for one of the MVAs in the main python script " << std::endl; + } + + return mvaEstimators; + } + + static std::vector getValueMapNames(const edm::VParameterSet& vConfig, std::string&& suffix) { + std::vector names; + for (auto& imva : vConfig) { + names.push_back(imva.getParameter("mvaName") + imva.getParameter("mvaTag") + suffix); + } + + return names; + } + void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; const edm::EDGetTokenT> srcToken_; @@ -65,35 +94,6 @@ namespace { iEvent.put(std::move(valMap), label); } - auto getMVAEstimators(const edm::VParameterSet& vConfig) { - std::vector> mvaEstimators; - - // Loop over the list of MVA configurations passed here from python and - // construct all requested MVA estimators. - for (auto& imva : vConfig) { - // The factory below constructs the MVA of the appropriate type based - // on the "mvaName" which is the name of the derived MVA class (plugin) - if (!imva.empty()) { - mvaEstimators.emplace_back( - AnyMVAEstimatorRun2Factory::get()->create(imva.getParameter("mvaName"), imva)); - - } else - throw cms::Exception(" MVA configuration not found: ") - << " failed to find proper configuration for one of the MVAs in the main python script " << std::endl; - } - - return mvaEstimators; - } - - std::vector getValueMapNames(const edm::VParameterSet& vConfig, std::string&& suffix) { - std::vector names; - for (auto& imva : vConfig) { - names.push_back(imva.getParameter("mvaName") + imva.getParameter("mvaTag") + suffix); - } - - return names; - } - template auto getKeysForValueMapsToken(edm::InputTag const& keysForValueMapsTag, edm::ConsumesCollector&& cc) { const bool tagGiven = !keysForValueMapsTag.label().empty(); diff --git a/RecoMuon/MuonIdentification/plugins/MuonSelectionTypeValueMapProducer.h b/RecoMuon/MuonIdentification/plugins/MuonSelectionTypeValueMapProducer.h index 89d4ad05dd4c6..2f0e4208d2bcc 100644 --- a/RecoMuon/MuonIdentification/plugins/MuonSelectionTypeValueMapProducer.h +++ b/RecoMuon/MuonIdentification/plugins/MuonSelectionTypeValueMapProducer.h @@ -36,7 +36,7 @@ class MuonSelectionTypeValueMapProducer : public edm::stream::EDProducer<> { muon::SelectionType selectionType_; }; -void MuonSelectionTypeValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { +inline void MuonSelectionTypeValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { // input muon collection edm::Handle muonsH; iEvent.getByToken(muonToken_, muonsH); diff --git a/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h index f651657612e10..461d010199c3f 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h +++ b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h @@ -10,22 +10,6 @@ #include #include "RealisticCluster.h" -namespace { - - float getDecayLength(unsigned int layer, unsigned int fhOffset, unsigned int bhOffset) { - constexpr float eeDecayLengthInLayer = 2.f; - constexpr float fhDecayLengthInLayer = 1.5f; - constexpr float bhDecayLengthInLayer = 1.f; - - if (layer <= fhOffset) - return eeDecayLengthInLayer; - else if (layer > fhOffset && layer <= bhOffset) - return fhDecayLengthInLayer; - else - return bhDecayLengthInLayer; - } -} // namespace - class RealisticHitToClusterAssociator { using Hit3DPosition = std::array; @@ -283,6 +267,19 @@ class RealisticHitToClusterAssociator { const std::vector& realisticClusters() const { return realisticSimClusters_; } private: + static float getDecayLength(unsigned int layer, unsigned int fhOffset, unsigned int bhOffset) { + constexpr float eeDecayLengthInLayer = 2.f; + constexpr float fhDecayLengthInLayer = 1.5f; + constexpr float bhDecayLengthInLayer = 1.f; + + if (layer <= fhOffset) + return eeDecayLengthInLayer; + else if (layer > fhOffset && layer <= bhOffset) + return fhDecayLengthInLayer; + else + return bhDecayLengthInLayer; + } + // the vector of the Realistic SimClusters std::vector realisticSimClusters_; // the vector of the Realistic SimClusters diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cc b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cc index 96381673388ca..5978ef8851c73 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cc +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cc @@ -1 +1,50 @@ -#include "CAHitNtupletGeneratorKernelsAlloc.h" +#include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h" + +#include "CAHitNtupletGeneratorKernels.h" + +template <> +#ifdef __CUDACC__ +void CAHitNtupletGeneratorKernelsGPU::allocateOnGPU(int32_t nHits, cudaStream_t stream) { +#else +void CAHitNtupletGeneratorKernelsCPU::allocateOnGPU(int32_t nHits, cudaStream_t stream) { +#endif + ////////////////////////////////////////////////////////// + // ALLOCATIONS FOR THE INTERMEDIATE RESULTS (STAYS ON WORKER) + ////////////////////////////////////////////////////////// + + device_theCellNeighbors_ = Traits::template make_unique(stream); + device_theCellTracks_ = Traits::template make_unique(stream); + +#ifdef GPU_DEBUG + std::cout << "Allocation for tuple building. N hits " << nHits << std::endl; +#endif + + nHits++; // storage requires one more counter; + assert(nHits > 0); + device_hitToTuple_ = Traits::template make_unique(stream); + device_hitToTupleStorage_ = Traits::template make_unique(nHits, stream); + hitToTupleView_.assoc = device_hitToTuple_.get(); + hitToTupleView_.offStorage = device_hitToTupleStorage_.get(); + hitToTupleView_.offSize = nHits; + + device_tupleMultiplicity_ = Traits::template make_unique(stream); + + device_storage_ = Traits::template make_unique(3, stream); + + device_hitTuple_apc_ = (cms::cuda::AtomicPairCounter*)device_storage_.get(); + device_hitToTuple_apc_ = (cms::cuda::AtomicPairCounter*)device_storage_.get() + 1; + device_nCells_ = (uint32_t*)(device_storage_.get() + 2); + + // FIXME: consider collapsing these 3 in one adhoc kernel + if constexpr (std::is_same::value) { + cudaCheck(cudaMemsetAsync(device_nCells_, 0, sizeof(uint32_t), stream)); + } else { + *device_nCells_ = 0; + } + cms::cuda::launchZero(device_tupleMultiplicity_.get(), stream); + cms::cuda::launchZero(hitToTupleView_, stream); // we may wish to keep it in the edm +#ifdef GPU_DEBUG + cudaDeviceSynchronize(); + cudaCheck(cudaGetLastError()); +#endif +} diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cu b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cu index 96381673388ca..68ee08d657e81 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cu +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.cu @@ -1 +1 @@ -#include "CAHitNtupletGeneratorKernelsAlloc.h" +#include "CAHitNtupletGeneratorKernelsAlloc.cc" diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.h b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.h deleted file mode 100644 index 5978ef8851c73..0000000000000 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsAlloc.h +++ /dev/null @@ -1,50 +0,0 @@ -#include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h" - -#include "CAHitNtupletGeneratorKernels.h" - -template <> -#ifdef __CUDACC__ -void CAHitNtupletGeneratorKernelsGPU::allocateOnGPU(int32_t nHits, cudaStream_t stream) { -#else -void CAHitNtupletGeneratorKernelsCPU::allocateOnGPU(int32_t nHits, cudaStream_t stream) { -#endif - ////////////////////////////////////////////////////////// - // ALLOCATIONS FOR THE INTERMEDIATE RESULTS (STAYS ON WORKER) - ////////////////////////////////////////////////////////// - - device_theCellNeighbors_ = Traits::template make_unique(stream); - device_theCellTracks_ = Traits::template make_unique(stream); - -#ifdef GPU_DEBUG - std::cout << "Allocation for tuple building. N hits " << nHits << std::endl; -#endif - - nHits++; // storage requires one more counter; - assert(nHits > 0); - device_hitToTuple_ = Traits::template make_unique(stream); - device_hitToTupleStorage_ = Traits::template make_unique(nHits, stream); - hitToTupleView_.assoc = device_hitToTuple_.get(); - hitToTupleView_.offStorage = device_hitToTupleStorage_.get(); - hitToTupleView_.offSize = nHits; - - device_tupleMultiplicity_ = Traits::template make_unique(stream); - - device_storage_ = Traits::template make_unique(3, stream); - - device_hitTuple_apc_ = (cms::cuda::AtomicPairCounter*)device_storage_.get(); - device_hitToTuple_apc_ = (cms::cuda::AtomicPairCounter*)device_storage_.get() + 1; - device_nCells_ = (uint32_t*)(device_storage_.get() + 2); - - // FIXME: consider collapsing these 3 in one adhoc kernel - if constexpr (std::is_same::value) { - cudaCheck(cudaMemsetAsync(device_nCells_, 0, sizeof(uint32_t), stream)); - } else { - *device_nCells_ = 0; - } - cms::cuda::launchZero(device_tupleMultiplicity_.get(), stream); - cms::cuda::launchZero(hitToTupleView_, stream); // we may wish to keep it in the edm -#ifdef GPU_DEBUG - cudaDeviceSynchronize(); - cudaCheck(cudaGetLastError()); -#endif -} diff --git a/RecoPixelVertexing/PixelVertexFinding/interface/FindPeakFastPV.h b/RecoPixelVertexing/PixelVertexFinding/interface/FindPeakFastPV.h index 75b3d496a638b..7a5cbe488d0f7 100644 --- a/RecoPixelVertexing/PixelVertexFinding/interface/FindPeakFastPV.h +++ b/RecoPixelVertexing/PixelVertexFinding/interface/FindPeakFastPV.h @@ -10,12 +10,12 @@ #include -float FindPeakFastPV(const std::vector &zProjections, - const std::vector &zWeights, - const float oldVertex, - const float m_zClusterWidth, - const float m_zClusterSearchArea, - const float m_weightCut) { +inline float FindPeakFastPV(const std::vector &zProjections, + const std::vector &zWeights, + const float oldVertex, + const float m_zClusterWidth, + const float m_zClusterSearchArea, + const float m_weightCut) { float centerWMax = oldVertex; if (m_zClusterWidth > 0 && m_zClusterSearchArea > 0) { std::vector::const_iterator itCenter = zProjections.begin(); diff --git a/RecoTracker/ConversionSeedGenerators/plugins/Conv4HitsReco.cc b/RecoTracker/ConversionSeedGenerators/plugins/Conv4HitsReco.cc deleted file mode 100644 index 267d2d13b7b48..0000000000000 --- a/RecoTracker/ConversionSeedGenerators/plugins/Conv4HitsReco.cc +++ /dev/null @@ -1,276 +0,0 @@ -#define Conv4HitsReco_cxx -#include "Conv4HitsReco.h" -#include "FWCore/MessageLogger/interface/MessageLogger.h" -#include -// 1 -> 4 -// 2 -> 3 -// 3 -> 2 -// 4 -> 1 - -double Conv4HitsReco::GetDm() { - return (Tq * O0 - T0 * Oq) / (Tq * Om - Tm * Oq); //dm -} - -double Conv4HitsReco::GetDq() { - return (Tm * O0 - T0 * Om) / (Tm * Oq - Tq * Om); //dq -} - -void Conv4HitsReco::SetLinSystCoeff(double m, double q) { - // dq*Tq + dm*Tm = T0 // Tangent condition - // dq*Oq + dm*Om = O0 // Orthogonality condition - - double sqrtEta2mm = sqrt(_eta2 + m * m); - double sqrtPi2qq = sqrt(_pi2 + q * q); - - double signT = 1.; - - Tq = -2. * pPN + 2. * m * pn + signT * 2. * q * sqrtEta2mm / sqrtPi2qq; - Tm = 2. * nPN + 2. * q * pn + signT * 2. * m * sqrtPi2qq / sqrtEta2mm; - T0 = PN2 - _eta2 - _pi2 - 2. * q * m * pn + 2. * q * pPN - 2. * m * nPN - signT * 2. * sqrtEta2mm * sqrtPi2qq; - - TVector3 vQminusM = q * unitVp - m * unitVn + vPminusN; - double QM = vQminusM.Mag(); - double pQM = unitVp * vQminusM; - double nQM = unitVn * vQminusM; - double NVQM = vNminusV * vQminusM; - - double signO = 1.; - - Oq = sqrtEta2mm * pQM / QM + m * pn + pNV; - Om = m * QM / sqrtEta2mm - signO * sqrtEta2mm * nQM / QM + nQM - nNV - m; - O0 = -signO * sqrtEta2mm * QM - m * nQM - NVQM; -} - -std::string Conv4HitsReco::qFromM_print(double m) { - std::stringstream ss; - TVector3 vPminusM = vPminusN - m * unitVn; - double m2 = m * m; - double nPM = unitVn * vPminusM; - double pPM = unitVp * vPminusM; - double NVPM = vNminusV * vPminusM; - - double alpha = (m * pn + pNV) * (m * pn + pNV) - _eta2 - m2; - double beta = m2 * pn * nPM + m * pn * NVPM + m * nPM * pNV - pPM * (_eta2 + m2) + pNV * NVPM; - double gamma = m2 * nPM * nPM + NVPM * NVPM + 2. * m * nPM * NVPM - vPminusM.Mag2() * (_eta2 + m2); - - double disc = sqrt(beta * beta - alpha * gamma); - - double q01 = (-beta + disc) / alpha; - double q02 = (-beta - disc) / alpha; - - ss << " m: " << m << " q01: " << std::setw(20) << q01 << " q02: " << std::setw(20) << q02 << "/n"; - return ss.str(); -} - -double Conv4HitsReco::qFromM(double m) { - LogDebug("Conv4HitsReco") << qFromM_print(m); - return 0.; -} - -TVector3 Conv4HitsReco::GetPlusCenter(double &radius) { - radius = plusRadius; - return plusCenter; -} - -TVector3 Conv4HitsReco::GetMinusCenter(double &radius) { - radius = minusRadius; - return minusCenter; -} - -// -//Point of intersection between two lines (each identified by a vector and a point) -TVector3 Conv4HitsReco::GetIntersection(TVector3 &V1, TVector3 &p1, TVector3 &V2, TVector3 &p2) { - TVector3 v1 = V1.Unit(); - TVector3 v2 = V2.Unit(); - - double v1v2 = v1 * v2; - return v2 * ((p1 - p2) * (v1v2 * v1 - v2) / (v1v2 * v1v2 - 1.)) + p2; -} - -double Conv4HitsReco::GetPtFromParamAndHitPair(double &m, double &eta) { - return 0.01 * 0.3 * 3.8 * sqrt(m * m + eta * eta); -} - -double Conv4HitsReco::GetPtMinusFromParam(double &m) { return GetPtFromParamAndHitPair(m, _eta); } - -double Conv4HitsReco::GetPtPlusFromParam(double &q) { return GetPtFromParamAndHitPair(q, _pi); } - -TVector3 Conv4HitsReco::GetConvVertexFromParams(double &m, double &q) { - TVector3 unitVQminusM = (plusCenter - minusCenter).Unit(); - TVector3 vtxViaPlus = vP + q * unitVp - plusRadius * unitVQminusM; - TVector3 vtxViaMinus = vN + m * unitVn + minusRadius * unitVQminusM; - - // return 0.5*(vN+m*unitVn+m*unitVQminusM+vP+q*unitVp-q*unitVQminusM); - LogDebug("Conv4HitsReco") << ">>>>>>>> Conversion vertex computed via Plus pair\n" - << vtxViaPlus.x() << "," << vtxViaPlus.y() << "," << vtxViaPlus.z() - << ">>>>>>>> Conversion vertex computed via Minus pair\n" - << vtxViaMinus.x() << "," << vtxViaMinus.y() << "," << vtxViaMinus.z(); - - return 0.5 * (vtxViaPlus + vtxViaMinus); -} - -int Conv4HitsReco::ConversionCandidate(TVector3 &vtx, double &ptplus, double &ptminus) { - double m; - double q; - int nits = 0; - - int isNotValidBefore = ComputeMaxLimits() + ComputeMinLimits(); - int isNotValidAfter = 0; - if (!isNotValidBefore) { - GuessStartingValues(m, q); - nits = mqFindByIteration(m, q); - - if (q > qMaxLimit || q < qMinLimit) { - isNotValidAfter = 1; - LogDebug("Conv4HitsReco") << ">>>>>>>> quad result not valid for q: qMin= " << qMinLimit << " q= " << q - << " qMax= " << qMaxLimit << "\n"; - } - if (m > mMaxLimit || m < mMinLimit) { - isNotValidAfter = 1; - LogDebug("Conv4HitsReco") << ">>>>>>>> quad result not valid for m: mMin= " << mMinLimit << " m= " << m - << " mMax= " << mMaxLimit << "\n"; - } - - ptminus = GetPtMinusFromParam(m); - ptplus = GetPtPlusFromParam(q); - minusCenter = vN + m * unitVn; - minusRadius = (hit4 - minusCenter).Mag(); - plusCenter = vP + q * unitVp; - plusRadius = (hit1 - plusCenter).Mag(); - convVtx = GetConvVertexFromParams(m, q); - vtx = convVtx; - } - - if (isNotValidBefore) - return 0; - if (IsNotValidForPtLimit(m, q, ptLegMinCut, ptLegMaxCut)) - return -1000; - if (IsNotValidForVtxPosition(maxVtxDistance)) - return -2000; - if (isNotValidAfter) - return -1 * nits; - return nits; -} - -int Conv4HitsReco::mqFindByIteration(double &m, double &q) { - int maxIte = maxNumberOfIterations; - double err = iterationStopRelThreshold; - double edm = 1.; - double edq = 1.; - int i = 0; - while (((edq > err) || (edm > err)) && (i < maxIte)) { - SetLinSystCoeff(m, q); - double dm = GetDm(); - double dq = GetDq(); - /* - while( m+dm > mMaxLimit || m+dm < mMinLimit || q+dq > qMaxLimit || q+dq < qMinLimit ){ - - LogDebug("Conv4HitsReco")<<">>>>>>>> Going outside limits, reducing increments \n"; - dm=dm/2.; - dq=dq/2.; - } - */ - m += dm; - q += dq; - edm = fabs(dm / m); - edq = fabs(dq / q); - LogDebug("Conv4HitsReco") << ">>>>>>>> Iteration " << i << " m: " << m << " q: " << q << " dm: " << dm - << " dq: " << dq << " edm: " << edm << " edq: " << edq << "\n"; - i++; - } - - return i; -} - -int Conv4HitsReco::ComputeMaxLimits() { - // Q max limit - TVector3 vVminusHit2Orthogonal = (vV - hit2).Orthogonal(); - TVector3 medianPointHit2V = 0.5 * (vV + hit2); - vQMaxLimit = GetIntersection(vVminusHit2Orthogonal, medianPointHit2V, unitVp, vP); - qMaxLimit = (vQMaxLimit - vP).Mag(); - - // M max limit - TVector3 vVminusHit3Orthogonal = (vV - hit3).Orthogonal(); - TVector3 medianPointHit3V = 0.5 * (vV + hit3); - vMMaxLimit = GetIntersection(vVminusHit3Orthogonal, medianPointHit3V, unitVn, vN); - mMaxLimit = (vMMaxLimit - vN).Mag(); - - LogDebug("Conv4HitsReco") << " >>>>>> Limits: qMax= " << qMaxLimit << " ==>pt " - << GetPtFromParamAndHitPair(qMaxLimit, _pi) << " mMax= " << mMaxLimit << " ==>pt " - << GetPtFromParamAndHitPair(mMaxLimit, _eta) << "\n"; - - return IsNotValidForPtLimit(mMaxLimit, qMaxLimit, ptLegMinCut, 100000.); //Max limit not applied here -} - -int Conv4HitsReco::IsNotValidForPtLimit(double m, double q, double ptmin, double ptmax) { - if (GetPtFromParamAndHitPair(q, _pi) < ptmin) - return 1; - if (GetPtFromParamAndHitPair(m, _eta) < ptmin) - return 1; - if (GetPtFromParamAndHitPair(q, _pi) > ptmax) - return 1; - if (GetPtFromParamAndHitPair(m, _eta) > ptmax) - return 1; - return 0; -} - -int Conv4HitsReco::IsNotValidForVtxPosition(double &maxDist) { - TVector3 hitAve = 0.25 * (hit1 + hit2 + hit3 + hit4); - if ((convVtx - hitAve).Mag() > maxDist) - return 1; - return 0; -} - -int Conv4HitsReco::ComputeMinLimits() { - //Evaluate if quad is valid and compute min limits - if (((vV - vQMaxLimit).Cross(vMMaxLimit - vQMaxLimit)).Z() > 0.) { - // - //Quad is invalid - LogDebug("Conv4HitsReco") << " >>>>>> Quad is invalid\n"; - return 1; - } else { - // - // Compute q and m Min limits - TVector3 vQMinLimit = GetIntersection(v1minus2, vMMaxLimit, unitVp, vP); - qMinLimit = (vQMinLimit - vP) * unitVp; - TVector3 vMMinLimit = GetIntersection(v3minus4, vQMaxLimit, unitVn, vN); - mMinLimit = (vMMinLimit - vN) * unitVn; - if (mMinLimit > mMaxLimit || qMinLimit > qMaxLimit) { - LogDebug("Conv4HitsReco") << " >>>>>> Quad is invalid. qMin= " << qMinLimit << " mMin= " << mMinLimit << "\n"; - return 2; - } - if (IsNotValidForPtLimit(mMinLimit, qMinLimit, -1000., ptLegMaxCut)) { //Min limit not applied here - return 2; - } - - LogDebug("Conv4HitsReco") << " >>>>>> Quad is valid. qMin= " << qMinLimit << " mMin= " << mMinLimit << "\n"; - return 0; - } -} - -int Conv4HitsReco::GuessStartingValues(double &m, double &q) { - /* - m = 0.5*(mMinLimit+mMaxLimit); - q = 0.5*(qMinLimit+qMaxLimit); - */ - - m = mMaxLimit; - q = qMaxLimit; - - LogDebug("Conv4HitsReco") << " >>>>>> Starting values: q= " << q << " m= " << m << "\n"; - - return 0; -} - -void Conv4HitsReco::Dump() { - LogDebug("Conv4HitsReco") << " ======================================= " - << "\n" - << " Photon Vertex: " << vV.x() << "," << vV.y() << "," << vV.z() << " Hit1: " << hit1.x() - << "," << hit1.y() << "," << hit1.z() << " Hit2: " << hit2.x() << "," << hit2.y() << "," - << hit2.z() << " Hit3: " << hit3.x() << "," << hit3.y() << "," << hit3.z() - << " Hit4: " << hit4.x() << "," << hit4.y() << "," << hit4.z() << " N: " << vN.x() << "," - << vN.y() << "," << vN.z() << " P: " << vP.x() << "," << vP.y() << "," << vP.z() - << " P-N: " << vPminusN.x() << "," << vP.y() << "," << vP.z() << " n: " << unitVn.x() << "," - << unitVn.y() << "," << unitVn.z() << " p: " << unitVp.x() << "," << unitVp.y() << "," - << unitVp.z() << " eta: " << _eta << " pi: " << _pi << "\n"; -} diff --git a/RecoTracker/ConversionSeedGenerators/plugins/Conv4HitsReco.h b/RecoTracker/ConversionSeedGenerators/plugins/Conv4HitsReco.h deleted file mode 100644 index d7e258b54ec8f..0000000000000 --- a/RecoTracker/ConversionSeedGenerators/plugins/Conv4HitsReco.h +++ /dev/null @@ -1,149 +0,0 @@ -#ifndef Conv4HitsReco_h -#define Conv4HitsReco_h - -#include -#include -#include - -#include - -class Conv4HitsReco { -public: - //Needed for construction - TVector3 vV; - TVector3 hit4; - TVector3 hit3; - TVector3 hit2; - TVector3 hit1; - - //Important quantities - TVector3 v3minus4; - TVector3 v1minus2; - TVector3 vN; - TVector3 vP; - TVector3 vPminusN; - TVector3 vPminusV; - TVector3 vNminusV; - TVector3 unitVn; - TVector3 unitVp; - double pn; - double pPN; - double nPN; - double PN; - double PN2; - double pNV; - double nNV; - double _eta; - double _pi; - double _eta2; - double _pi2; - TVector3 vMMaxLimit; - TVector3 vQMaxLimit; - double qMaxLimit; - double mMaxLimit; - double qMinLimit; - double mMinLimit; - TVector3 plusCenter; - TVector3 minusCenter; - TVector3 convVtx; - double plusRadius; - double minusRadius; - - double iterationStopRelThreshold; - int maxNumberOfIterations; - double maxVtxDistance; - double ptLegMinCut; - double ptLegMaxCut; - double ptPhotMaxCut; - - std::string qFromM_print(double m); - double qFromM(double); //For debugging purposes - - //Elements of the linearized system matrix - double Tq, Tm, T0; - double Oq, Om, O0; - int ComputeMaxLimits(); - int ComputeMinLimits(); - int GuessStartingValues(double &, double &); - int mqFindByIteration(double &, double &); - TVector3 GetIntersection(TVector3 &, TVector3 &, TVector3 &, TVector3 &); - TVector3 GetPlusCenter(double &); - TVector3 GetMinusCenter(double &); - void SetPtLegMinCut(double val) { ptLegMinCut = val; }; - void SetPtLegMaxCut(double val) { ptLegMaxCut = val; }; - void SetPtPhotMaxCut(double val) { ptPhotMaxCut = val; }; - void SetLinSystCoeff(double, double); - void Set(double); - void SetIterationStopRelThreshold(double val) { iterationStopRelThreshold = val; }; - void SetMaxNumberOfIterations(int val) { maxNumberOfIterations = val; }; - void SetMaxVtxDistance(int val) { maxVtxDistance = val; }; - double GetDq(); - double GetDm(); - double GetPtFromParamAndHitPair(double &, double &); - double GetPtPlusFromParam(double &); - double GetPtMinusFromParam(double &); - TVector3 GetConvVertexFromParams(double &, double &); - int IsNotValidForPtLimit(double, double, double, double); - int IsNotValidForVtxPosition(double &); - - int ConversionCandidate(TVector3 &, double &, double &); - void Dump(); - - Conv4HitsReco(TVector3 &, TVector3 &, TVector3 &, TVector3 &, TVector3 &); - ~Conv4HitsReco(); -}; - -#endif - -#ifdef Conv4HitsReco_cxx -Conv4HitsReco::Conv4HitsReco(TVector3 &vPhotVertex, TVector3 &h1, TVector3 &h2, TVector3 &h3, TVector3 &h4) { - vV = vPhotVertex; - hit4 = h4; - hit3 = h3; - hit2 = h2; - hit1 = h1; - - //Let's stay in the transverse plane... - vV.SetZ(0.); - hit4.SetZ(0.); - hit3.SetZ(0.); - hit2.SetZ(0.); - hit1.SetZ(0.); - - //Filling important quantities - vN.SetXYZ(0.5 * hit4.X() + 0.5 * hit3.X(), 0.5 * hit4.Y() + 0.5 * hit3.Y(), 0.5 * hit4.Z() + 0.5 * hit3.Z()); - vP.SetXYZ(0.5 * hit2.X() + 0.5 * hit1.X(), 0.5 * hit2.Y() + 0.5 * hit1.Y(), 0.5 * hit2.Z() + 0.5 * hit1.Z()); - vPminusN = vP - vN; - vPminusV = vP - vV; - vNminusV = vN - vV; - v3minus4 = hit3 - hit4; - v1minus2 = hit1 - hit2; - unitVn = v3minus4.Orthogonal().Unit(); - unitVp = v1minus2.Orthogonal().Unit(); - pPN = unitVp * vPminusN; - nPN = unitVn * vPminusN; - pNV = unitVp * vNminusV; - nNV = unitVn * vNminusV; - pn = unitVp * unitVn; - PN = vPminusN.Mag(); - PN2 = vPminusN.Mag2(); - - _eta = 0.5 * (hit3 - hit4).Mag(); - _pi = 0.5 * (hit1 - hit2).Mag(); - _eta2 = _eta * _eta; - _pi2 = _pi * _pi; - - //Default values for parameters - SetIterationStopRelThreshold(0.0005); - SetMaxNumberOfIterations(10); - SetPtLegMinCut(0.2); - SetPtLegMaxCut(20.); - SetPtPhotMaxCut(30.); - SetMaxVtxDistance(20.); -} - -Conv4HitsReco::~Conv4HitsReco() { - // std::cout << " Bye..." << std::endl; -} - -#endif // #ifdef Conv4HitsReco_h diff --git a/RecoTracker/DebugTools/plugins/CkfDebugger.cc b/RecoTracker/DebugTools/plugins/CkfDebugger.cc index f047b8b70ab68..58fcae6ac243f 100644 --- a/RecoTracker/DebugTools/plugins/CkfDebugger.cc +++ b/RecoTracker/DebugTools/plugins/CkfDebugger.cc @@ -15,7 +15,6 @@ #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" #include "DataFormats/SiStripDetId/interface/StripSubdetector.h" -#include "GluedDetFromDetUnit.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" #include "DataFormats/SiStripDetId/interface/StripSubdetector.h" @@ -34,6 +33,11 @@ using namespace std; +inline DetId gluedId(const DetId& du) { + unsigned int mask = ~3; // mask the last two bits + return DetId(du.rawId() & mask); +} + CkfDebugger::CkfDebugger(edm::EventSetup const& es, edm::ConsumesCollector&& iC) : trackerHitAssociatorConfig_(std::move(iC)), totSeeds(0) { file = new TFile("out.root", "recreate"); diff --git a/RecoTracker/DebugTools/plugins/GluedDetFromDetUnit.h b/RecoTracker/DebugTools/plugins/GluedDetFromDetUnit.h deleted file mode 100644 index b0da7ce37a0f5..0000000000000 --- a/RecoTracker/DebugTools/plugins/GluedDetFromDetUnit.h +++ /dev/null @@ -1,6 +0,0 @@ -#include "DataFormats/DetId/interface/DetId.h" - -DetId gluedId(const DetId& du) { - unsigned int mask = ~3; // mask the last two bits - return DetId(du.rawId() & mask); -} diff --git a/RecoVertex/AdaptiveVertexFinder/interface/TrackVertexArbitration.h b/RecoVertex/AdaptiveVertexFinder/interface/TrackVertexArbitration.h index db6b491148052..85f079bb72daf 100644 --- a/RecoVertex/AdaptiveVertexFinder/interface/TrackVertexArbitration.h +++ b/RecoVertex/AdaptiveVertexFinder/interface/TrackVertexArbitration.h @@ -119,10 +119,10 @@ bool TrackVertexArbitration::trackFilterArbitrator(const reco::TransientTra return true; } -float trackWeight(const reco::Vertex &sv, const reco::TransientTrack &track) { +inline float trackWeight(const reco::Vertex &sv, const reco::TransientTrack &track) { return sv.trackWeight(track.trackBaseRef()); } -float trackWeight(const reco::VertexCompositePtrCandidate &sv, const reco::TransientTrack &tt) { +inline float trackWeight(const reco::VertexCompositePtrCandidate &sv, const reco::TransientTrack &tt) { const reco::CandidatePtrTransientTrack *cptt = dynamic_cast(tt.basicTransientTrack()); if (cptt == nullptr) diff --git a/RecoVertex/AdaptiveVertexFinder/plugins/TemplatedVertexArbitrator.h b/RecoVertex/AdaptiveVertexFinder/plugins/TemplatedVertexArbitrator.h deleted file mode 100644 index 00b9d08e9c918..0000000000000 --- a/RecoVertex/AdaptiveVertexFinder/plugins/TemplatedVertexArbitrator.h +++ /dev/null @@ -1,153 +0,0 @@ -#ifndef TemplatedVertexArbitrator_h -#define TemplatedVertexArbitrator_h -#include -#include - -#include "TrackingTools/TransientTrack/interface/TransientTrack.h" -#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" -#include "TrackingTools/Records/interface/TransientTrackRecord.h" -#include "TrackingTools/IPTools/interface/IPTools.h" - -#include "FWCore/Framework/interface/stream/EDProducer.h" -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/Utilities/interface/InputTag.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/Utilities/interface/ESGetToken.h" - -#include "DataFormats/Common/interface/Handle.h" -#include "DataFormats/TrackReco/interface/Track.h" -#include "DataFormats/TrackReco/interface/TrackFwd.h" -#include "DataFormats/VertexReco/interface/Vertex.h" -#include "DataFormats/VertexReco/interface/VertexFwd.h" -#include "RecoVertex/VertexTools/interface/VertexDistance3D.h" - -#include "DataFormats/TrackReco/interface/Track.h" -#include "DataFormats/TrackReco/interface/TrackFwd.h" -#include "DataFormats/VertexReco/interface/Vertex.h" -#include "DataFormats/VertexReco/interface/VertexFwd.h" -#include "DataFormats/BeamSpot/interface/BeamSpot.h" - -#include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h" -#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h" -#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackCompatibilityEstimator.h" -#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h" -#include "DataFormats/Math/interface/deltaR.h" - -#include "RecoVertex/ConfigurableVertexReco/interface/ConfigurableVertexReconstructor.h" -#include "RecoVertex/AdaptiveVertexFinder/interface/TrackVertexArbitration.h" -#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h" - -#include "RecoVertex/AdaptiveVertexFinder/interface/TTHelpers.h" - -//#define VTXDEBUG - -const unsigned int nTracks(const reco::Vertex &sv) { return sv.nTracks(); } -const unsigned int nTracks(const reco::VertexCompositePtrCandidate &sv) { return sv.numberOfSourceCandidatePtrs(); } - -template -class TemplatedVertexArbitrator : public edm::stream::EDProducer<> { -public: - typedef std::vector Product; - TemplatedVertexArbitrator(const edm::ParameterSet ¶ms); - - static void fillDescriptions(edm::ConfigurationDescriptions &cdesc) { - edm::ParameterSetDescription pdesc; - pdesc.add("beamSpot", edm::InputTag("offlineBeamSpot")); - pdesc.add("primaryVertices", edm::InputTag("offlinePrimaryVertices")); - if (std::is_same::value) { - pdesc.add("tracks", edm::InputTag("generalTracks")); - pdesc.add("secondaryVertices", edm::InputTag("vertexMerger")); - } else if (std::is_same::value) { - pdesc.add("tracks", edm::InputTag("particleFlow")); - pdesc.add("secondaryVertices", edm::InputTag("candidateVertexMerger")); - } else { - pdesc.add("tracks", edm::InputTag("generalTracks")); - pdesc.add("secondaryVertices", edm::InputTag("vertexMerger")); - } - pdesc.add("dLenFraction", 0.333); - pdesc.add("dRCut", 0.4); - pdesc.add("distCut", 0.04); - pdesc.add("sigCut", 5.0); - pdesc.add("fitterSigmacut", 3.0); - pdesc.add("fitterTini", 256); - pdesc.add("fitterRatio", 0.25); - pdesc.add("trackMinLayers", 4); - pdesc.add("trackMinPt", 0.4); - pdesc.add("trackMinPixels", 1); - pdesc.add("maxTimeSignificance", 3.5); - if (std::is_same::value) { - cdesc.add("trackVertexArbitratorDefault", pdesc); - } else if (std::is_same::value) { - cdesc.add("candidateVertexArbitratorDefault", pdesc); - } else { - cdesc.addDefault(pdesc); - } - } - - void produce(edm::Event &event, const edm::EventSetup &es) override; - -private: - bool trackFilter(const reco::TrackRef &track) const; - - edm::EDGetTokenT token_primaryVertex; - edm::EDGetTokenT token_secondaryVertex; - edm::EDGetTokenT token_tracks; - edm::EDGetTokenT token_beamSpot; - edm::ESGetToken token_trackBuilder; - std::unique_ptr > theArbitrator; -}; - -template -TemplatedVertexArbitrator::TemplatedVertexArbitrator(const edm::ParameterSet ¶ms) { - token_primaryVertex = consumes(params.getParameter("primaryVertices")); - token_secondaryVertex = consumes(params.getParameter("secondaryVertices")); - token_beamSpot = consumes(params.getParameter("beamSpot")); - token_tracks = consumes(params.getParameter("tracks")); - token_trackBuilder = - esConsumes(edm::ESInputTag("", "TransientTrackBuilder")); - produces(); - theArbitrator.reset(new TrackVertexArbitration(params)); -} - -template -void TemplatedVertexArbitrator::produce(edm::Event &event, const edm::EventSetup &es) { - using namespace reco; - - edm::Handle secondaryVertices; - event.getByToken(token_secondaryVertex, secondaryVertices); - Product theSecVertexColl = *(secondaryVertices.product()); - - edm::Handle primaryVertices; - event.getByToken(token_primaryVertex, primaryVertices); - - auto recoVertices = std::make_unique(); - if (!primaryVertices->empty()) { - const reco::Vertex &pv = (*primaryVertices)[0]; - - edm::Handle tracks; - event.getByToken(token_tracks, tracks); - - edm::ESHandle trackBuilder = es.getHandle(token_trackBuilder); - - edm::Handle beamSpot; - event.getByToken(token_beamSpot, beamSpot); - - std::vector selectedTracks; - for (typename InputContainer::const_iterator track = tracks->begin(); track != tracks->end(); ++track) { - selectedTracks.push_back(tthelpers::buildTT(tracks, trackBuilder, track - tracks->begin())); - } - - // const edm::RefVector< TrackCollection > tracksForArbitration= selectedTracks; - Product theRecoVertices = theArbitrator->trackVertexArbitrator(beamSpot, pv, selectedTracks, theSecVertexColl); - - for (unsigned int ivtx = 0; ivtx < theRecoVertices.size(); ivtx++) { - if (!(nTracks(theRecoVertices[ivtx]) > 1)) - continue; - recoVertices->push_back(theRecoVertices[ivtx]); - } - } - event.put(std::move(recoVertices)); -} - -#endif diff --git a/RecoVertex/AdaptiveVertexFinder/plugins/VertexArbitrators.cc b/RecoVertex/AdaptiveVertexFinder/plugins/VertexArbitrators.cc index b154d24246a84..68d982afced75 100644 --- a/RecoVertex/AdaptiveVertexFinder/plugins/VertexArbitrators.cc +++ b/RecoVertex/AdaptiveVertexFinder/plugins/VertexArbitrators.cc @@ -1,19 +1,155 @@ +#include +#include + +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/IPTools/interface/IPTools.h" + #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/Utilities/interface/InputTag.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/ESGetToken.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "RecoVertex/VertexTools/interface/VertexDistance3D.h" -#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h" #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/VertexReco/interface/Vertex.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" -#include "RecoVertex/AdaptiveVertexFinder/plugins/TemplatedVertexArbitrator.h" -#include "DataFormats/Common/interface/View.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" + +#include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h" +#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h" +#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackCompatibilityEstimator.h" +#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h" +#include "DataFormats/Math/interface/deltaR.h" + +#include "RecoVertex/ConfigurableVertexReco/interface/ConfigurableVertexReconstructor.h" +#include "RecoVertex/AdaptiveVertexFinder/interface/TrackVertexArbitration.h" +#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h" + +#include "RecoVertex/AdaptiveVertexFinder/interface/TTHelpers.h" //#define VTXDEBUG +inline const unsigned int nTracks(const reco::Vertex &sv) { return sv.nTracks(); } +inline const unsigned int nTracks(const reco::VertexCompositePtrCandidate &sv) { + return sv.numberOfSourceCandidatePtrs(); +} + +template +class TemplatedVertexArbitrator : public edm::stream::EDProducer<> { +public: + typedef std::vector Product; + TemplatedVertexArbitrator(const edm::ParameterSet ¶ms); + + static void fillDescriptions(edm::ConfigurationDescriptions &cdesc) { + edm::ParameterSetDescription pdesc; + pdesc.add("beamSpot", edm::InputTag("offlineBeamSpot")); + pdesc.add("primaryVertices", edm::InputTag("offlinePrimaryVertices")); + if (std::is_same::value) { + pdesc.add("tracks", edm::InputTag("generalTracks")); + pdesc.add("secondaryVertices", edm::InputTag("vertexMerger")); + } else if (std::is_same::value) { + pdesc.add("tracks", edm::InputTag("particleFlow")); + pdesc.add("secondaryVertices", edm::InputTag("candidateVertexMerger")); + } else { + pdesc.add("tracks", edm::InputTag("generalTracks")); + pdesc.add("secondaryVertices", edm::InputTag("vertexMerger")); + } + pdesc.add("dLenFraction", 0.333); + pdesc.add("dRCut", 0.4); + pdesc.add("distCut", 0.04); + pdesc.add("sigCut", 5.0); + pdesc.add("fitterSigmacut", 3.0); + pdesc.add("fitterTini", 256); + pdesc.add("fitterRatio", 0.25); + pdesc.add("trackMinLayers", 4); + pdesc.add("trackMinPt", 0.4); + pdesc.add("trackMinPixels", 1); + pdesc.add("maxTimeSignificance", 3.5); + if (std::is_same::value) { + cdesc.add("trackVertexArbitratorDefault", pdesc); + } else if (std::is_same::value) { + cdesc.add("candidateVertexArbitratorDefault", pdesc); + } else { + cdesc.addDefault(pdesc); + } + } + + void produce(edm::Event &event, const edm::EventSetup &es) override; + +private: + bool trackFilter(const reco::TrackRef &track) const; + + edm::EDGetTokenT token_primaryVertex; + edm::EDGetTokenT token_secondaryVertex; + edm::EDGetTokenT token_tracks; + edm::EDGetTokenT token_beamSpot; + edm::ESGetToken token_trackBuilder; + std::unique_ptr > theArbitrator; +}; + +template +TemplatedVertexArbitrator::TemplatedVertexArbitrator(const edm::ParameterSet ¶ms) { + token_primaryVertex = consumes(params.getParameter("primaryVertices")); + token_secondaryVertex = consumes(params.getParameter("secondaryVertices")); + token_beamSpot = consumes(params.getParameter("beamSpot")); + token_tracks = consumes(params.getParameter("tracks")); + token_trackBuilder = + esConsumes(edm::ESInputTag("", "TransientTrackBuilder")); + produces(); + theArbitrator.reset(new TrackVertexArbitration(params)); +} + +template +void TemplatedVertexArbitrator::produce(edm::Event &event, const edm::EventSetup &es) { + using namespace reco; + + edm::Handle secondaryVertices; + event.getByToken(token_secondaryVertex, secondaryVertices); + Product theSecVertexColl = *(secondaryVertices.product()); + + edm::Handle primaryVertices; + event.getByToken(token_primaryVertex, primaryVertices); + + auto recoVertices = std::make_unique(); + if (!primaryVertices->empty()) { + const reco::Vertex &pv = (*primaryVertices)[0]; + + edm::Handle tracks; + event.getByToken(token_tracks, tracks); + + edm::ESHandle trackBuilder = es.getHandle(token_trackBuilder); + + edm::Handle beamSpot; + event.getByToken(token_beamSpot, beamSpot); + + std::vector selectedTracks; + for (typename InputContainer::const_iterator track = tracks->begin(); track != tracks->end(); ++track) { + selectedTracks.push_back(tthelpers::buildTT(tracks, trackBuilder, track - tracks->begin())); + } + + // const edm::RefVector< TrackCollection > tracksForArbitration= selectedTracks; + Product theRecoVertices = theArbitrator->trackVertexArbitrator(beamSpot, pv, selectedTracks, theSecVertexColl); + + for (unsigned int ivtx = 0; ivtx < theRecoVertices.size(); ivtx++) { + if (!(nTracks(theRecoVertices[ivtx]) > 1)) + continue; + recoVertices->push_back(theRecoVertices[ivtx]); + } + } + event.put(std::move(recoVertices)); +} + typedef TemplatedVertexArbitrator TrackVertexArbitrator; typedef TemplatedVertexArbitrator, reco::VertexCompositePtrCandidate> CandidateVertexArbitrator; diff --git a/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianEXT.icc b/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianEXT.icc index 29b5a9cdd5f27..84b7e171dd10b 100644 --- a/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianEXT.icc +++ b/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianEXT.icc @@ -1,13 +1,9 @@ - -void -AnalyticalCurvilinearJacobian::computeFullJacobian -(const GlobalTrajectoryParameters& globalParameters, - const GlobalPoint& x, - const GlobalVector& p, - const GlobalVector& h, - const double& s) -{ - +// NOLINTNEXTLINE(misc-definitions-in-headers) +void AnalyticalCurvilinearJacobian::computeFullJacobian(const GlobalTrajectoryParameters& globalParameters, + const GlobalPoint& x, + const GlobalVector& p, + const GlobalVector& h, + const double& s) { //GlobalVector p1 = fts.momentum().unit(); GlobalVector p1 = globalParameters.momentum().unit(); GlobalVector p2 = p.unit(); @@ -20,54 +16,56 @@ AnalyticalCurvilinearJacobian::computeFullJacobian //double qbp = fts.signedInverseMomentum(); double qbp = globalParameters.signedInverseMomentum(); double absS = s; - + // calculate transport matrix // Origin: TRPRFN - double cosl0 = p1.perp(); double cosl1 = 1./p2.perp(); + double cosl0 = p1.perp(); + double cosl1 = 1. / p2.perp(); - // define average magnetic field and gradient + // define average magnetic field and gradient // at initial point - inlike TRPRFN GlobalVector hnf = h.unit(); double qp = -h.mag(); // double q = -h.mag()*qbp; - double q = qp*qbp; + double q = qp * qbp; - double theta = q*absS; double sint,cost; vdt::fast_sincos(theta,sint,cost); + double theta = q * absS; + double sint, cost; + vdt::fast_sincos(theta, sint, cost); Vec3D t1; Vec3D t2; - + Vec3D hn; Vec3D dx; - - for ( int i=0; i!=4;++i) { - t1[i]=p1.basicVector().v[i]; - t2[i]=p2.basicVector().v[i]; - hn[i]=hnf.basicVector().v[i]; - dx[i]=dxf.basicVector().v[i]; + + for (int i = 0; i != 4; ++i) { + t1[i] = p1.basicVector().v[i]; + t2[i] = p2.basicVector().v[i]; + hn[i] = hnf.basicVector().v[i]; + dx[i] = dxf.basicVector().v[i]; } - double gamma = dot(hn,t2); - Vec3D an = cross3(hn,t2); + double gamma = dot(hn, t2); + Vec3D an = cross3(hn, t2); - typedef unsigned long long v4sl __attribute__ ((vector_size (32))); - v4sl mask{1,0,5,4}; - Vec4D t = __builtin_shuffle(t1,t2,mask); - Vec4D tt = t*t; + typedef unsigned long long v4sl __attribute__((vector_size(32))); + v4sl mask{1, 0, 5, 4}; + Vec4D t = __builtin_shuffle(t1, t2, mask); + Vec4D tt = t * t; double au1 = std::sqrt(tt[0] + tt[1]); double au2 = std::sqrt(tt[2] + tt[3]); - Vec4D au{-au1,au1,-au2,au2}; - Vec4D uu = t/au; + Vec4D au{-au1, au1, -au2, au2}; + Vec4D uu = t / au; - Vec2D u1{uu[0],uu[1]}; - Vec2D u2{uu[2],uu[3]}; - - Vec3D u13{uu[0],uu[1],0,0}; - Vec3D v1 = cross3(t1,u13); - Vec3D u23{uu[2],uu[3],0,0}; - Vec3D v2 = cross3(t2,u23); + Vec2D u1{uu[0], uu[1]}; + Vec2D u2{uu[2], uu[3]}; + Vec3D u13{uu[0], uu[1], 0, 0}; + Vec3D v1 = cross3(t1, u13); + Vec3D u23{uu[2], uu[3], 0, 0}; + Vec3D v2 = cross3(t2, u23); // now prepare the transport matrix // pp only needed in high-p case (WA) @@ -76,122 +74,84 @@ AnalyticalCurvilinearJacobian::computeFullJacobian // moved up (where -h.mag() is needed() // double qp = q*pp; - double anv = -dot(hn,u23); - double anu = dot(hn,v2); + double anv = -dot(hn, u23); + double anu = dot(hn, v2); + + double omcost = 1. - cost; + double tmsint = theta - sint; - double omcost = 1. - cost; double tmsint = theta - sint; - - - Vec3D hu = cross3(hn,u13); - Vec3D hv = cross3(hn,v1); + Vec3D hu = cross3(hn, u13); + Vec3D hv = cross3(hn, v1); - // 1/p - doesn't change since |p1| = |p2| - theJacobian(0,0) = 1.; for (auto i=1;i<5; ++i) theJacobian(0,i)=0.; + theJacobian(0, 0) = 1.; + for (auto i = 1; i < 5; ++i) + theJacobian(0, i) = 0.; // lambda - - theJacobian(1,0) = -qp*anv*dot(t2,dx); - - theJacobian(1,1) = - cost*dot(v1,v2) + - sint*dot(hv,v2) + - omcost*dot(hn,v1) * dot(hn,v2) + - anv*(-sint*dot(v1,t2) + - omcost*dot(v1,an) - - tmsint*gamma*dot(hn,v1) - ) - ; - - theJacobian(1,2) = - cost*dot2(u1,v2) + - sint*dot(hu,v2) + - omcost*dot2(hn,u1) * dot(hn,v2) + - anv*(-sint*dot2(u1,t2) + - omcost*dot2(u1,an) - - tmsint*gamma*dot2(hn,u1) - ) - ; - theJacobian(1,2) *= cosl0; - - theJacobian(1,3) = -q*anv*dot2(u1,t2); - - theJacobian(1,4) = -q*anv*dot(v1,t2); + + theJacobian(1, 0) = -qp * anv * dot(t2, dx); + + theJacobian(1, 1) = cost * dot(v1, v2) + sint * dot(hv, v2) + omcost * dot(hn, v1) * dot(hn, v2) + + anv * (-sint * dot(v1, t2) + omcost * dot(v1, an) - tmsint * gamma * dot(hn, v1)); + + theJacobian(1, 2) = cost * dot2(u1, v2) + sint * dot(hu, v2) + omcost * dot2(hn, u1) * dot(hn, v2) + + anv * (-sint * dot2(u1, t2) + omcost * dot2(u1, an) - tmsint * gamma * dot2(hn, u1)); + theJacobian(1, 2) *= cosl0; + + theJacobian(1, 3) = -q * anv * dot2(u1, t2); + + theJacobian(1, 4) = -q * anv * dot(v1, t2); // phi - theJacobian(2,0) = -qp*anu*cosl1*dot(t2,dx); - - theJacobian(2,1) = - cost*dot(xy(v1),u2) + - sint*dot(xy(hv),u2) + - omcost*dot(hn,v1)*dot2(hn,u2) + - anu*(-sint*dot(v1,t2) + - omcost*dot(v1,an) - - tmsint*gamma*dot(hn,v1) - ) - ; - theJacobian(2,1) *= cosl1; - - theJacobian(2,2) = - cost*dot(u1,u2) + - sint*dot2(hu,u2) + - omcost*dot2(hn,u1)*dot2(hn,u2) + - anu*(-sint*dot2(u1,t2) + - omcost*dot2(u1,an) - - tmsint*gamma*dot2(hn,u1) - ) - ; - theJacobian(2,2) *= cosl1*cosl0; - - theJacobian(2,3) = -q*anu*cosl1*dot2(u1,t2); - - theJacobian(2,4) = -q*anu*cosl1*dot(v1,t2); + theJacobian(2, 0) = -qp * anu * cosl1 * dot(t2, dx); + + theJacobian(2, 1) = cost * dot(xy(v1), u2) + sint * dot(xy(hv), u2) + omcost * dot(hn, v1) * dot2(hn, u2) + + anu * (-sint * dot(v1, t2) + omcost * dot(v1, an) - tmsint * gamma * dot(hn, v1)); + theJacobian(2, 1) *= cosl1; + + theJacobian(2, 2) = cost * dot(u1, u2) + sint * dot2(hu, u2) + omcost * dot2(hn, u1) * dot2(hn, u2) + + anu * (-sint * dot2(u1, t2) + omcost * dot2(u1, an) - tmsint * gamma * dot2(hn, u1)); + theJacobian(2, 2) *= cosl1 * cosl0; + + theJacobian(2, 3) = -q * anu * cosl1 * dot2(u1, t2); + + theJacobian(2, 4) = -q * anu * cosl1 * dot(v1, t2); // yt + double overQ = 1. / q; - double overQ = 1./q; + theJacobian(3, 1) = (sint * dot(xy(v1), u2) + omcost * dot2(hv, u2) + tmsint * dot(hn, v1) * dot2(hn, u2)) * overQ; - theJacobian(3,1) = (sint*dot(xy(v1),u2) + - omcost*dot2(hv,u2) + - tmsint*dot(hn,v1)*dot2(hn,u2) - )*overQ; + theJacobian(3, 2) = + (sint * dot(u1, u2) + omcost * dot2(hu, u2) + tmsint * dot2(hn, u1) * dot2(hn, u2)) * (cosl0 * overQ); - theJacobian(3,2) = (sint*dot(u1,u2) + - omcost*dot2(hu,u2) + - tmsint*dot2(hn,u1)*dot2(hn,u2) - )*(cosl0*overQ); + theJacobian(3, 3) = dot(u1, u2); - theJacobian(3,3) = dot(u1,u2); - - theJacobian(3,4) = dot(xy(v1),u2); + theJacobian(3, 4) = dot(xy(v1), u2); // zt - theJacobian(4,1) = (sint*dot(v1,v2) + - omcost*dot(hv,v2) + - tmsint*dot(hn,v1)*dot(hn,v2) - )*overQ; + theJacobian(4, 1) = (sint * dot(v1, v2) + omcost * dot(hv, v2) + tmsint * dot(hn, v1) * dot(hn, v2)) * overQ; - theJacobian(4,2) = (sint*dot(u1,xy(v2)) + - omcost*dot(hu,v2) + - tmsint*dot2(hn,u1)*dot(hn,v2) - )*(cosl0*overQ); + theJacobian(4, 2) = + (sint * dot(u1, xy(v2)) + omcost * dot(hu, v2) + tmsint * dot2(hn, u1) * dot(hn, v2)) * (cosl0 * overQ); - theJacobian(4,3) = dot2(u1,v2); + theJacobian(4, 3) = dot2(u1, v2); - theJacobian(4,4) = dot(v1,v2); + theJacobian(4, 4) = dot(v1, v2); //double cutCriterion = abs(s/fts.momentum().mag()); // double cutCriterion = fabs(s/globalParameters.momentum().mag()); - double cutCriterion = std::abs(s*qbp); - const double limit = 5.; // valid for propagations with effectively float precision + double cutCriterion = std::abs(s * qbp); + const double limit = 5.; // valid for propagations with effectively float precision if (cutCriterion > limit) { - double pp = 1./qbp; - theJacobian(3,0) = pp*dot2(u2,dx); -// theJacobian(4,0) = -pp*dot(v2,dx); // was wrong sign??? -/* + double pp = 1. / qbp; + theJacobian(3, 0) = pp * dot2(u2, dx); + // theJacobian(4,0) = -pp*dot(v2,dx); // was wrong sign??? + /* It seems so. The effect was noticed analysing the distribution of reduced chi2 of general tracks vs eta. A +3% difference was noticed whem using the - sign instead of the plus, @@ -199,34 +159,31 @@ in the region .75 < |eta| < 1.5, in particular for <1GeV tracks. Indeed, the reduced chi2 is only one of the symptoms: many other quantities were affected (momentum for example). In addition additional tracks were reconstructed. */ - theJacobian(4,0) = pp*dot(v2,dx); + theJacobian(4, 0) = pp * dot(v2, dx); - } - else { - Vec3D hp1 = cross3(hn,t1); - double temp1 = dot2(hp1,u2); - Vec3D ghnmp = gamma*hn - t1; - double temp2 = dot(xy(ghnmp),u2); + } else { + Vec3D hp1 = cross3(hn, t1); + double temp1 = dot2(hp1, u2); + Vec3D ghnmp = gamma * hn - t1; + double temp2 = dot(xy(ghnmp), u2); - double qps = qp*s; - double h2 = qps*qbp; - double h3 = (-1./8.) * h2; + double qps = qp * s; + double h2 = qps * qbp; + double h3 = (-1. / 8.) * h2; double secondOrder41 = 0.5 * temp1; - double thirdOrder41 = (1./3.) *temp2; + double thirdOrder41 = (1. / 3.) * temp2; double fourthOrder41 = h3 * temp1; - theJacobian(3,0) = (s*qps)*(secondOrder41 + h2*(thirdOrder41 + fourthOrder41)); + theJacobian(3, 0) = (s * qps) * (secondOrder41 + h2 * (thirdOrder41 + fourthOrder41)); - double temp3 = dot(hp1,v2); - double temp4 = dot(ghnmp,v2); + double temp3 = dot(hp1, v2); + double temp4 = dot(ghnmp, v2); double secondOrder51 = 0.5 * temp3; - double thirdOrder51 = (1./3.) *temp4; - double fourthOrder51 = h3 * temp3; - theJacobian(4,0) = (s*qps)*(secondOrder51 + h2*(thirdOrder51 + fourthOrder51)); + double thirdOrder51 = (1. / 3.) * temp4; + double fourthOrder51 = h3 * temp3; + theJacobian(4, 0) = (s * qps) * (secondOrder51 + h2 * (thirdOrder51 + fourthOrder51)); } - - // end of TRPRFN } diff --git a/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianSSE.icc b/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianSSE.icc index 8a5562658ab9e..1dd4f27f55868 100644 --- a/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianSSE.icc +++ b/TrackingTools/AnalyticalJacobians/src/AnalyticalCurvilinearJacobianSSE.icc @@ -1,4 +1,4 @@ - +// NOLINTNEXTLINE(misc-definitions-in-headers) void AnalyticalCurvilinearJacobian::computeFullJacobian(const GlobalTrajectoryParameters& globalParameters, const GlobalPoint& x, const GlobalVector& p, diff --git a/TrackingTools/TrackFitters/plugins/KFFittingSmoother.h b/TrackingTools/TrackFitters/plugins/KFFittingSmoother.h deleted file mode 100644 index 8a3e89bb670a8..0000000000000 --- a/TrackingTools/TrackFitters/plugins/KFFittingSmoother.h +++ /dev/null @@ -1,423 +0,0 @@ -/** \class KFFittingSmoother - * A TrajectorySmoother that rpeats the forward fit before smoothing. - * This is necessary e.g. when the seed introduced a bias (by using - * a beam contraint etc.). Ported from ORCA - * - * \author todorov, cerati - */ - -#include "TrackingTools/PatternTools/interface/TrajectorySmoother.h" -#include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h" -#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" - -#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" -#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" - -#include "FWCore/MessageLogger/interface/MessageLogger.h" -#include "FWCore/Utilities/interface/isFinite.h" -#include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h" -#include "CommonTools/Utils/interface/DynArray.h" - -namespace { - - struct KFFittingSmootherParam { - explicit KFFittingSmootherParam(const edm::ParameterSet& conf) - : theEstimateCut(conf.getParameter("EstimateCut")), - theMaxFractionOutliers(conf.getParameter("MaxFractionOutliers")), - theMaxNumberOfOutliers(conf.getParameter("MaxNumberOfOutliers")), - theNoOutliersBeginEnd(conf.getParameter("NoOutliersBeginEnd")), - theMinDof(conf.getParameter("MinDof")), - theMinNumberOfHits(conf.getParameter("MinNumberOfHits")), - rejectTracksFlag(conf.getParameter("RejectTracks")), - breakTrajWith2ConsecutiveMissing(conf.getParameter("BreakTrajWith2ConsecutiveMissing")), - noInvalidHitsBeginEnd(conf.getParameter("NoInvalidHitsBeginEnd")) {} - - double theEstimateCut; - - float theMaxFractionOutliers; - int theMaxNumberOfOutliers; - bool theNoOutliersBeginEnd; - int theMinDof; - - int theMinNumberOfHits; - bool rejectTracksFlag; - bool breakTrajWith2ConsecutiveMissing; - bool noInvalidHitsBeginEnd; - }; - - class KFFittingSmoother final : public TrajectoryFitter, private KFFittingSmootherParam { - public: - ~KFFittingSmoother() override {} - - KFFittingSmoother(const TrajectoryFitter& aFitter, - const TrajectorySmoother& aSmoother, - const edm::ParameterSet& conf) - : KFFittingSmootherParam(conf), theFitter(aFitter.clone()), theSmoother(aSmoother.clone()) {} - - private: - static void fillDescriptions(edm::ParameterSetDescription& desc) { - desc.add("EstimateCut", -1); - desc.add("MaxFractionOutliers", 0.3); - desc.add("MaxNumberOfOutliers", 3); - desc.add("MinDof", 2); - desc.add("NoOutliersBeginEnd", false); - desc.add("MinNumberOfHits", 5); - desc.add("RejectTracks", true); - desc.add("BreakTrajWith2ConsecutiveMissing", true); - desc.add("NoInvalidHitsBeginEnd", true); - desc.add("LogPixelProbabilityCut", 0); - } - - Trajectory fitOne(const Trajectory& t, fitType type) const override; - Trajectory fitOne(const TrajectorySeed& aSeed, - const RecHitContainer& hits, - const TrajectoryStateOnSurface& firstPredTsos, - fitType type) const override; - Trajectory fitOne(const TrajectorySeed& aSeed, const RecHitContainer& hits, fitType type) const override; - - const TrajectoryFitter* fitter() const { return theFitter.get(); } - const TrajectorySmoother* smoother() const { return theSmoother.get(); } - - std::unique_ptr clone() const override { - return std::unique_ptr(new KFFittingSmoother(*theFitter, *theSmoother, *this)); - } - - void setHitCloner(TkCloner const* hc) override { - theFitter->setHitCloner(hc); - theSmoother->setHitCloner(hc); - } - - KFFittingSmoother(const TrajectoryFitter& aFitter, - const TrajectorySmoother& aSmoother, - KFFittingSmootherParam const& other) - : KFFittingSmootherParam(other), theFitter(aFitter.clone()), theSmoother(aSmoother.clone()) {} - - Trajectory smoothingStep(Trajectory&& fitted) const { - if (theEstimateCut > 0) { - // remove "outlier" at the end of Traj - while ( - !fitted.empty() && fitted.foundHits() >= theMinNumberOfHits && - (!fitted.lastMeasurement().recHitR().isValid() || (fitted.lastMeasurement().recHitR().det() != nullptr && - fitted.lastMeasurement().estimate() > theEstimateCut))) - fitted.pop(); - if (fitted.foundHits() < theMinNumberOfHits) - return Trajectory(); - } - return theSmoother->trajectory(fitted); - } - - private: - const std::unique_ptr theFitter; - const std::unique_ptr theSmoother; - - /// Method to check that the trajectory has no NaN in the states and chi2 - static bool checkForNans(const Trajectory& theTraj); - - friend class KFFittingSmootherESProducer; - }; - - // #define VI_DEBUG - -#ifdef VI_DEBUG -#define DPRINT(x) std::cout << x << ": " -#define PRINT std::cout -#else -#define DPRINT(x) LogTrace(x) -#define PRINT LogTrace("") -#endif - - Trajectory KFFittingSmoother::fitOne(const Trajectory& t, fitType type) const { - if (!t.isValid()) - return Trajectory(); - return smoothingStep(theFitter->fitOne(t, type)); - } - - bool KFFittingSmoother::checkForNans(const Trajectory& theTraj) { - if (edm::isNotFinite(theTraj.chiSquared())) - return false; - auto const& vtm = theTraj.measurements(); - for (auto const& tm : vtm) { - if (edm::isNotFinite(tm.estimate())) - return false; - auto const& v = tm.updatedState().localParameters().vector(); - for (int i = 0; i < 5; ++i) - if (edm::isNotFinite(v[i])) - return false; - if (!tm.updatedState().curvilinearError().posDef()) - return false; //not a NaN by itself, but will lead to one - auto const& m = tm.updatedState().curvilinearError().matrix(); - for (int i = 0; i < 5; ++i) - for (int j = 0; j < (i + 1); ++j) - if (edm::isNotFinite(m(i, j))) - return false; - } - return true; - } - - namespace { - inline void print(const std::string& header, const TrajectoryStateOnSurface& tsos) { - DPRINT("TrackFitters") << header << tsos.globalPosition().perp() << ' ' << tsos.globalPosition().z() << ' ' - << 1. / tsos.signedInverseMomentum() << ' ' << 1. / tsos.transverseCurvature() << ' ' - << tsos.globalMomentum().eta() << std::endl; - } - } // namespace - - Trajectory KFFittingSmoother::fitOne(const TrajectorySeed& aSeed, - const RecHitContainer& hits, - const TrajectoryStateOnSurface& firstPredTsos, - fitType type) const { - LogDebug("TrackFitters") << "In KFFittingSmoother::fit"; - - print("firstPred ", firstPredTsos); - - if (hits.empty()) - return Trajectory(); - - RecHitContainer myHits = hits; - Trajectory tmp_first; - - //call the fitter - Trajectory smoothed = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos)); - - do { -#ifdef EDM_ML_DEBUG - //if no outliers the fit is done only once - for (unsigned int j = 0; j < myHits.size(); j++) { - if (myHits[j]->det()) - LogTrace("TrackFitters") << "hit #:" << j + 1 << " rawId=" << myHits[j]->det()->geographicalId().rawId() - << " validity=" << myHits[j]->isValid(); - else - LogTrace("TrackFitters") << "hit #:" << j + 1 << " Hit with no Det information"; - } -#endif - -#if defined(VI_DEBUG) || defined(EDM_ML_DEBUG) - if (smoothed.isValid()) { - print("first state ", smoothed.firstMeasurement().updatedState()); - print("last state ", smoothed.lastMeasurement().updatedState()); - } -#endif - - bool hasNaN = false; - if (!smoothed.isValid() || (hasNaN = !checkForNans(smoothed)) || (smoothed.foundHits() < theMinNumberOfHits)) { - if (hasNaN) - edm::LogWarning("TrackNaN") << "Track has NaN or the cov is not pos-definite"; - if (smoothed.foundHits() < theMinNumberOfHits) - LogTrace("TrackFitters") << "smoothed.foundHits() trajectory rejected with nhits/chi2 " << smoothed.foundHits() - << '/' << smoothed.chiSquared() << "\n"; - if (rejectTracksFlag) { - return Trajectory(); - } else { - std::swap(smoothed, tmp_first); // if first attempt, tmp_first would be invalid anyway - DPRINT("TrackFitters") << "smoothed invalid => returning orignal trajectory with nhits/chi2 " - << smoothed.foundHits() << '/' << smoothed.chiSquared() << "\n"; - } - break; - } -#ifdef EDM_ML_DEBUG - else { - LogTrace("TrackFitters") << "dump hits after smoothing"; - Trajectory::DataContainer meas = smoothed.measurements(); - for (Trajectory::DataContainer::iterator it = meas.begin(); it != meas.end(); ++it) { - LogTrace("TrackFitters") << "hit #" << meas.end() - it - 1 << " validity=" << it->recHit()->isValid() - << " det=" << it->recHit()->geographicalId().rawId(); - } - } -#endif - - if (myHits.size() != smoothed.measurements().size()) - DPRINT("TrackFitters") << "lost hits. before/after: " << myHits.size() << '/' << smoothed.measurements().size() - << "\n"; - - if (theEstimateCut <= 0) - break; - - // Check if there are outliers - - auto msize = smoothed.measurements().size(); - declareDynArray(unsigned int, msize, bad); - unsigned int nbad = 0; - unsigned int ind = 0; - unsigned int lastValid = smoothed.measurements().size(); - for (auto const& tm : smoothed.measurements()) { - if (tm.estimate() > theEstimateCut && - tm.recHitR().det() != nullptr // do not consider outliers constraints and other special "hits" - ) - bad[nbad++] = ind; - if (ind < lastValid && tm.recHitR().det() != nullptr && tm.recHitR().isValid()) - lastValid = ind; - ++ind; - } - - if (0 == nbad) - break; - - DPRINT("TrackFitters") << "size/found/outliers list " << smoothed.measurements().size() << '/' - << smoothed.foundHits() << ' ' << nbad << ": "; - for (auto i = 0U; i < nbad; ++i) - PRINT << bad[i] << ','; - PRINT << std::endl; - - if ( - // (smoothed.foundHits() == theMinNumberOfHits) || - int(nbad) > theMaxNumberOfOutliers || float(nbad) > theMaxFractionOutliers * float(smoothed.foundHits())) { - DPRINT("TrackFitters") << "smoothed low quality => trajectory with nhits/chi2 " << smoothed.foundHits() << '/' - << smoothed.chiSquared() << "\n"; - PRINT << "try to remove " << lastValid << std::endl; - nbad = 0; // try to short the traj... (below lastValid will be added) - - // do not perform outliers rejection if track is already low quality - /* - if ( rejectTracksFlag && (smoothed.chiSquared() > theEstimateCut*smoothed.ndof()) ) { - DPRINT("TrackFitters") << "smoothed low quality => trajectory rejected with nhits/chi2 " << smoothed.foundHits() << '/' << smoothed.chiSquared() << "\n"; - return Trajectory(); - } else { - DPRINT("TrackFitters") << "smoothed low quality => return original trajectory with nhits/chi2 " << smoothed.foundHits() << '/' << smoothed.chiSquared() << "\n"; - } - break; - */ - } - - // always add last valid hit as outlier candidate - bad[nbad++] = lastValid; - - // if ( (smoothed.ndof()::max(); - - auto loc = nbad; - for (auto i = 0U; i < nbad; ++i) { - auto j = NHits - bad[i] - 1; - assert(myHits[j]->geographicalId() == smoothed.measurements()[bad[i]].recHitR().geographicalId()); - auto removedHit = myHits[j]; - myHits[j] = std::make_shared(*removedHit->det(), TrackingRecHit::missing); - smoothedCand[i] = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos)); - myHits[j] = removedHit; - if (smoothedCand[i].isValid() && smoothedCand[i].chiSquared() < minChi2) { - minChi2 = smoothedCand[i].chiSquared(); - loc = i; - } - } - - if (loc == nbad) { - DPRINT("TrackFitters") << "New trajectories all invalid" - << "\n"; - return Trajectory(); - } - - DPRINT("TrackFitters") << "outlier removed " << bad[loc] << '/' << minChi2 << " was " << smoothed.chiSquared() - << "\n"; - - if (minChi2 > smoothed.chiSquared()) { - DPRINT("TrackFitters") << "removing outlier makes chi2 worse " << minChi2 << '/' << smoothed.chiSquared() - << "\nOri: "; - for (auto const& tm : smoothed.measurements()) - PRINT << tm.recHitR().geographicalId() << '/' << tm.estimate() << ' '; - PRINT << "\nNew: "; - for (auto const& tm : smoothedCand[loc].measurements()) - PRINT << tm.recHitR().geographicalId() << '/' << tm.estimate() << ' '; - PRINT << "\n"; - - // return Trajectory(); - // break; - } - - std::swap(smoothed, tmp_first); - myHits[NHits - bad[loc] - 1] = - std::make_shared(*myHits[NHits - bad[loc] - 1]->det(), TrackingRecHit::missing); - std::swap(smoothed, smoothedCand[loc]); - // firstTry=false; - - DPRINT("TrackFitters") << "new trajectory with nhits/chi2 " << smoothed.foundHits() << '/' - << smoothed.chiSquared() << "\n"; - - // Look if there are two consecutive invalid hits FIXME: take into account split matched hits!!! - if (breakTrajWith2ConsecutiveMissing) { - unsigned int firstinvalid = myHits.size(); - for (unsigned int j = 0; j < myHits.size() - 1; ++j) { - if (((myHits[j]->type() == TrackingRecHit::missing) && (myHits[j]->geographicalId().rawId() != 0)) && - ((myHits[j + 1]->type() == TrackingRecHit::missing) && (myHits[j + 1]->geographicalId().rawId() != 0)) && - ((myHits[j]->geographicalId().rawId() & (~3)) != - (myHits[j + 1]->geographicalId().rawId() & (~3))) // same gluedDet - ) { - firstinvalid = j; - DPRINT("TrackFitters") << "Found two consecutive missing hits. First invalid: " << firstinvalid << "\n"; - break; - } - } - - //reject all the hits after the last valid before two consecutive invalid (missing) hits - //hits are sorted in the same order as in the track candidate FIXME?????? - if (firstinvalid != myHits.size()) { - myHits.erase(myHits.begin() + firstinvalid, myHits.end()); - smoothed = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos)); - DPRINT("TrackFitters") << "Trajectory shortened " << smoothed.foundHits() << '/' << smoothed.chiSquared() - << "\n"; - } - } - - } // do - while (true); - - if (smoothed.isValid()) { - if (noInvalidHitsBeginEnd && !smoothed.empty() //should we send a warning ? - ) { - // discard latest dummy measurements - if (!smoothed.empty() && !smoothed.lastMeasurement().recHitR().isValid()) - LogTrace("TrackFitters") << "Last measurement is invalid"; - - while (!smoothed.empty() && !smoothed.lastMeasurement().recHitR().isValid()) - smoothed.pop(); - - //remove the invalid hits at the begin of the trajectory - if (!smoothed.empty() && !smoothed.firstMeasurement().recHitR().isValid()) { - LogTrace("TrackFitters") << "First measurement is in`valid"; - Trajectory tmpTraj(smoothed.seed(), smoothed.direction()); - Trajectory::DataContainer& meas = smoothed.measurements(); - auto it = meas.begin(); - for (; it != meas.end(); ++it) - if (it->recHitR().isValid()) - break; - //push the first valid measurement and set the same global chi2 - tmpTraj.push(std::move(*it), smoothed.chiSquared()); - - for (auto itt = it + 1; itt != meas.end(); ++itt) - tmpTraj.push(std::move(*itt), 0); //add all the other measurements - - std::swap(smoothed, tmpTraj); - - } // if ( !smoothed.firstMeasurement().recHit()->isValid() ) - - } // if ( noInvalidHitsBeginEnd ) - - LogTrace("TrackFitters") << "end: returning smoothed trajectory with chi2=" << smoothed.chiSquared(); - - //LogTrace("TrackFitters") << "dump hits before return"; - //Trajectory::DataContainer meas = smoothed.measurements(); - //for (Trajectory::DataContainer::iterator it=meas.begin();it!=meas.end();++it) { - //LogTrace("TrackFitters") << "hit #" << meas.end()-it-1 << " validity=" << it->recHit()->isValid() - //<< " det=" << it->recHit()->geographicalId().rawId(); - //} - } - - return smoothed; - } - - Trajectory KFFittingSmoother::fitOne(const TrajectorySeed& aSeed, const RecHitContainer& hits, fitType type) const { - throw cms::Exception("TrackFitters", - "KFFittingSmoother::fit(TrajectorySeed, ) not implemented"); - - return Trajectory(); - } - -} // namespace diff --git a/TrackingTools/TrackFitters/plugins/KFFittingSmootherESProducer.cc b/TrackingTools/TrackFitters/plugins/KFFittingSmootherESProducer.cc index 8c54a80e9eaa1..696b5c277fd16 100644 --- a/TrackingTools/TrackFitters/plugins/KFFittingSmootherESProducer.cc +++ b/TrackingTools/TrackFitters/plugins/KFFittingSmootherESProducer.cc @@ -1,7 +1,428 @@ +/** \class KFFittingSmoother + * A TrajectorySmoother that rpeats the forward fit before smoothing. + * This is necessary e.g. when the seed introduced a bias (by using + * a beam contraint etc.). Ported from ORCA + * + * \author todorov, cerati + */ +#include "TrackingTools/PatternTools/interface/TrajectorySmoother.h" +#include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" -// to be included only here... -#include "KFFittingSmoother.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/isFinite.h" +#include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h" +#include "CommonTools/Utils/interface/DynArray.h" + +namespace { + + struct KFFittingSmootherParam { + explicit KFFittingSmootherParam(const edm::ParameterSet& conf) + : theEstimateCut(conf.getParameter("EstimateCut")), + theMaxFractionOutliers(conf.getParameter("MaxFractionOutliers")), + theMaxNumberOfOutliers(conf.getParameter("MaxNumberOfOutliers")), + theNoOutliersBeginEnd(conf.getParameter("NoOutliersBeginEnd")), + theMinDof(conf.getParameter("MinDof")), + theMinNumberOfHits(conf.getParameter("MinNumberOfHits")), + rejectTracksFlag(conf.getParameter("RejectTracks")), + breakTrajWith2ConsecutiveMissing(conf.getParameter("BreakTrajWith2ConsecutiveMissing")), + noInvalidHitsBeginEnd(conf.getParameter("NoInvalidHitsBeginEnd")) {} + + double theEstimateCut; + + float theMaxFractionOutliers; + int theMaxNumberOfOutliers; + bool theNoOutliersBeginEnd; + int theMinDof; + + int theMinNumberOfHits; + bool rejectTracksFlag; + bool breakTrajWith2ConsecutiveMissing; + bool noInvalidHitsBeginEnd; + }; + + class KFFittingSmoother final : public TrajectoryFitter, private KFFittingSmootherParam { + public: + ~KFFittingSmoother() override {} + + KFFittingSmoother(const TrajectoryFitter& aFitter, + const TrajectorySmoother& aSmoother, + const edm::ParameterSet& conf) + : KFFittingSmootherParam(conf), theFitter(aFitter.clone()), theSmoother(aSmoother.clone()) {} + + private: + static void fillDescriptions(edm::ParameterSetDescription& desc) { + desc.add("EstimateCut", -1); + desc.add("MaxFractionOutliers", 0.3); + desc.add("MaxNumberOfOutliers", 3); + desc.add("MinDof", 2); + desc.add("NoOutliersBeginEnd", false); + desc.add("MinNumberOfHits", 5); + desc.add("RejectTracks", true); + desc.add("BreakTrajWith2ConsecutiveMissing", true); + desc.add("NoInvalidHitsBeginEnd", true); + desc.add("LogPixelProbabilityCut", 0); + } + + Trajectory fitOne(const Trajectory& t, fitType type) const override; + Trajectory fitOne(const TrajectorySeed& aSeed, + const RecHitContainer& hits, + const TrajectoryStateOnSurface& firstPredTsos, + fitType type) const override; + Trajectory fitOne(const TrajectorySeed& aSeed, const RecHitContainer& hits, fitType type) const override; + + const TrajectoryFitter* fitter() const { return theFitter.get(); } + const TrajectorySmoother* smoother() const { return theSmoother.get(); } + + std::unique_ptr clone() const override { + return std::unique_ptr(new KFFittingSmoother(*theFitter, *theSmoother, *this)); + } + + void setHitCloner(TkCloner const* hc) override { + theFitter->setHitCloner(hc); + theSmoother->setHitCloner(hc); + } + + KFFittingSmoother(const TrajectoryFitter& aFitter, + const TrajectorySmoother& aSmoother, + KFFittingSmootherParam const& other) + : KFFittingSmootherParam(other), theFitter(aFitter.clone()), theSmoother(aSmoother.clone()) {} + + Trajectory smoothingStep(Trajectory&& fitted) const { + if (theEstimateCut > 0) { + // remove "outlier" at the end of Traj + while ( + !fitted.empty() && fitted.foundHits() >= theMinNumberOfHits && + (!fitted.lastMeasurement().recHitR().isValid() || (fitted.lastMeasurement().recHitR().det() != nullptr && + fitted.lastMeasurement().estimate() > theEstimateCut))) + fitted.pop(); + if (fitted.foundHits() < theMinNumberOfHits) + return Trajectory(); + } + return theSmoother->trajectory(fitted); + } + + private: + const std::unique_ptr theFitter; + const std::unique_ptr theSmoother; + + /// Method to check that the trajectory has no NaN in the states and chi2 + static bool checkForNans(const Trajectory& theTraj); + + friend class KFFittingSmootherESProducer; + }; + + // #define VI_DEBUG + +#ifdef VI_DEBUG +#define DPRINT(x) std::cout << x << ": " +#define PRINT std::cout +#else +#define DPRINT(x) LogTrace(x) +#define PRINT LogTrace("") +#endif + + inline Trajectory KFFittingSmoother::fitOne(const Trajectory& t, fitType type) const { + if (!t.isValid()) + return Trajectory(); + return smoothingStep(theFitter->fitOne(t, type)); + } + + inline bool KFFittingSmoother::checkForNans(const Trajectory& theTraj) { + if (edm::isNotFinite(theTraj.chiSquared())) + return false; + auto const& vtm = theTraj.measurements(); + for (auto const& tm : vtm) { + if (edm::isNotFinite(tm.estimate())) + return false; + auto const& v = tm.updatedState().localParameters().vector(); + for (int i = 0; i < 5; ++i) + if (edm::isNotFinite(v[i])) + return false; + if (!tm.updatedState().curvilinearError().posDef()) + return false; //not a NaN by itself, but will lead to one + auto const& m = tm.updatedState().curvilinearError().matrix(); + for (int i = 0; i < 5; ++i) + for (int j = 0; j < (i + 1); ++j) + if (edm::isNotFinite(m(i, j))) + return false; + } + return true; + } + + namespace { + inline void print(const std::string& header, const TrajectoryStateOnSurface& tsos) { + DPRINT("TrackFitters") << header << tsos.globalPosition().perp() << ' ' << tsos.globalPosition().z() << ' ' + << 1. / tsos.signedInverseMomentum() << ' ' << 1. / tsos.transverseCurvature() << ' ' + << tsos.globalMomentum().eta() << std::endl; + } + } // namespace + + inline Trajectory KFFittingSmoother::fitOne(const TrajectorySeed& aSeed, + const RecHitContainer& hits, + const TrajectoryStateOnSurface& firstPredTsos, + fitType type) const { + LogDebug("TrackFitters") << "In KFFittingSmoother::fit"; + + print("firstPred ", firstPredTsos); + + if (hits.empty()) + return Trajectory(); + + RecHitContainer myHits = hits; + Trajectory tmp_first; + + //call the fitter + Trajectory smoothed = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos)); + + do { +#ifdef EDM_ML_DEBUG + //if no outliers the fit is done only once + for (unsigned int j = 0; j < myHits.size(); j++) { + if (myHits[j]->det()) + LogTrace("TrackFitters") << "hit #:" << j + 1 << " rawId=" << myHits[j]->det()->geographicalId().rawId() + << " validity=" << myHits[j]->isValid(); + else + LogTrace("TrackFitters") << "hit #:" << j + 1 << " Hit with no Det information"; + } +#endif + +#if defined(VI_DEBUG) || defined(EDM_ML_DEBUG) + if (smoothed.isValid()) { + print("first state ", smoothed.firstMeasurement().updatedState()); + print("last state ", smoothed.lastMeasurement().updatedState()); + } +#endif + + bool hasNaN = false; + if (!smoothed.isValid() || (hasNaN = !checkForNans(smoothed)) || (smoothed.foundHits() < theMinNumberOfHits)) { + if (hasNaN) + edm::LogWarning("TrackNaN") << "Track has NaN or the cov is not pos-definite"; + if (smoothed.foundHits() < theMinNumberOfHits) + LogTrace("TrackFitters") << "smoothed.foundHits() trajectory rejected with nhits/chi2 " << smoothed.foundHits() + << '/' << smoothed.chiSquared() << "\n"; + if (rejectTracksFlag) { + return Trajectory(); + } else { + std::swap(smoothed, tmp_first); // if first attempt, tmp_first would be invalid anyway + DPRINT("TrackFitters") << "smoothed invalid => returning orignal trajectory with nhits/chi2 " + << smoothed.foundHits() << '/' << smoothed.chiSquared() << "\n"; + } + break; + } +#ifdef EDM_ML_DEBUG + else { + LogTrace("TrackFitters") << "dump hits after smoothing"; + Trajectory::DataContainer meas = smoothed.measurements(); + for (Trajectory::DataContainer::iterator it = meas.begin(); it != meas.end(); ++it) { + LogTrace("TrackFitters") << "hit #" << meas.end() - it - 1 << " validity=" << it->recHit()->isValid() + << " det=" << it->recHit()->geographicalId().rawId(); + } + } +#endif + + if (myHits.size() != smoothed.measurements().size()) + DPRINT("TrackFitters") << "lost hits. before/after: " << myHits.size() << '/' << smoothed.measurements().size() + << "\n"; + + if (theEstimateCut <= 0) + break; + + // Check if there are outliers + + auto msize = smoothed.measurements().size(); + declareDynArray(unsigned int, msize, bad); + unsigned int nbad = 0; + unsigned int ind = 0; + unsigned int lastValid = smoothed.measurements().size(); + for (auto const& tm : smoothed.measurements()) { + if (tm.estimate() > theEstimateCut && + tm.recHitR().det() != nullptr // do not consider outliers constraints and other special "hits" + ) + bad[nbad++] = ind; + if (ind < lastValid && tm.recHitR().det() != nullptr && tm.recHitR().isValid()) + lastValid = ind; + ++ind; + } + + if (0 == nbad) + break; + + DPRINT("TrackFitters") << "size/found/outliers list " << smoothed.measurements().size() << '/' + << smoothed.foundHits() << ' ' << nbad << ": "; + for (auto i = 0U; i < nbad; ++i) + PRINT << bad[i] << ','; + PRINT << std::endl; + + if ( + // (smoothed.foundHits() == theMinNumberOfHits) || + int(nbad) > theMaxNumberOfOutliers || float(nbad) > theMaxFractionOutliers * float(smoothed.foundHits())) { + DPRINT("TrackFitters") << "smoothed low quality => trajectory with nhits/chi2 " << smoothed.foundHits() << '/' + << smoothed.chiSquared() << "\n"; + PRINT << "try to remove " << lastValid << std::endl; + nbad = 0; // try to short the traj... (below lastValid will be added) + + // do not perform outliers rejection if track is already low quality + /* + if ( rejectTracksFlag && (smoothed.chiSquared() > theEstimateCut*smoothed.ndof()) ) { + DPRINT("TrackFitters") << "smoothed low quality => trajectory rejected with nhits/chi2 " << smoothed.foundHits() << '/' << smoothed.chiSquared() << "\n"; + return Trajectory(); + } else { + DPRINT("TrackFitters") << "smoothed low quality => return original trajectory with nhits/chi2 " << smoothed.foundHits() << '/' << smoothed.chiSquared() << "\n"; + } + break; + */ + } + + // always add last valid hit as outlier candidate + bad[nbad++] = lastValid; + + // if ( (smoothed.ndof()::max(); + + auto loc = nbad; + for (auto i = 0U; i < nbad; ++i) { + auto j = NHits - bad[i] - 1; + assert(myHits[j]->geographicalId() == smoothed.measurements()[bad[i]].recHitR().geographicalId()); + auto removedHit = myHits[j]; + myHits[j] = std::make_shared(*removedHit->det(), TrackingRecHit::missing); + smoothedCand[i] = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos)); + myHits[j] = removedHit; + if (smoothedCand[i].isValid() && smoothedCand[i].chiSquared() < minChi2) { + minChi2 = smoothedCand[i].chiSquared(); + loc = i; + } + } + + if (loc == nbad) { + DPRINT("TrackFitters") << "New trajectories all invalid" + << "\n"; + return Trajectory(); + } + + DPRINT("TrackFitters") << "outlier removed " << bad[loc] << '/' << minChi2 << " was " << smoothed.chiSquared() + << "\n"; + + if (minChi2 > smoothed.chiSquared()) { + DPRINT("TrackFitters") << "removing outlier makes chi2 worse " << minChi2 << '/' << smoothed.chiSquared() + << "\nOri: "; + for (auto const& tm : smoothed.measurements()) + PRINT << tm.recHitR().geographicalId() << '/' << tm.estimate() << ' '; + PRINT << "\nNew: "; + for (auto const& tm : smoothedCand[loc].measurements()) + PRINT << tm.recHitR().geographicalId() << '/' << tm.estimate() << ' '; + PRINT << "\n"; + + // return Trajectory(); + // break; + } + + std::swap(smoothed, tmp_first); + myHits[NHits - bad[loc] - 1] = + std::make_shared(*myHits[NHits - bad[loc] - 1]->det(), TrackingRecHit::missing); + std::swap(smoothed, smoothedCand[loc]); + // firstTry=false; + + DPRINT("TrackFitters") << "new trajectory with nhits/chi2 " << smoothed.foundHits() << '/' + << smoothed.chiSquared() << "\n"; + + // Look if there are two consecutive invalid hits FIXME: take into account split matched hits!!! + if (breakTrajWith2ConsecutiveMissing) { + unsigned int firstinvalid = myHits.size(); + for (unsigned int j = 0; j < myHits.size() - 1; ++j) { + if (((myHits[j]->type() == TrackingRecHit::missing) && (myHits[j]->geographicalId().rawId() != 0)) && + ((myHits[j + 1]->type() == TrackingRecHit::missing) && (myHits[j + 1]->geographicalId().rawId() != 0)) && + ((myHits[j]->geographicalId().rawId() & (~3)) != + (myHits[j + 1]->geographicalId().rawId() & (~3))) // same gluedDet + ) { + firstinvalid = j; + DPRINT("TrackFitters") << "Found two consecutive missing hits. First invalid: " << firstinvalid << "\n"; + break; + } + } + + //reject all the hits after the last valid before two consecutive invalid (missing) hits + //hits are sorted in the same order as in the track candidate FIXME?????? + if (firstinvalid != myHits.size()) { + myHits.erase(myHits.begin() + firstinvalid, myHits.end()); + smoothed = smoothingStep(theFitter->fitOne(aSeed, myHits, firstPredTsos)); + DPRINT("TrackFitters") << "Trajectory shortened " << smoothed.foundHits() << '/' << smoothed.chiSquared() + << "\n"; + } + } + + } // do + while (true); + + if (smoothed.isValid()) { + if (noInvalidHitsBeginEnd && !smoothed.empty() //should we send a warning ? + ) { + // discard latest dummy measurements + if (!smoothed.empty() && !smoothed.lastMeasurement().recHitR().isValid()) + LogTrace("TrackFitters") << "Last measurement is invalid"; + + while (!smoothed.empty() && !smoothed.lastMeasurement().recHitR().isValid()) + smoothed.pop(); + + //remove the invalid hits at the begin of the trajectory + if (!smoothed.empty() && !smoothed.firstMeasurement().recHitR().isValid()) { + LogTrace("TrackFitters") << "First measurement is in`valid"; + Trajectory tmpTraj(smoothed.seed(), smoothed.direction()); + Trajectory::DataContainer& meas = smoothed.measurements(); + auto it = meas.begin(); + for (; it != meas.end(); ++it) + if (it->recHitR().isValid()) + break; + //push the first valid measurement and set the same global chi2 + tmpTraj.push(std::move(*it), smoothed.chiSquared()); + + for (auto itt = it + 1; itt != meas.end(); ++itt) + tmpTraj.push(std::move(*itt), 0); //add all the other measurements + + std::swap(smoothed, tmpTraj); + + } // if ( !smoothed.firstMeasurement().recHit()->isValid() ) + + } // if ( noInvalidHitsBeginEnd ) + + LogTrace("TrackFitters") << "end: returning smoothed trajectory with chi2=" << smoothed.chiSquared(); + + //LogTrace("TrackFitters") << "dump hits before return"; + //Trajectory::DataContainer meas = smoothed.measurements(); + //for (Trajectory::DataContainer::iterator it=meas.begin();it!=meas.end();++it) { + //LogTrace("TrackFitters") << "hit #" << meas.end()-it-1 << " validity=" << it->recHit()->isValid() + //<< " det=" << it->recHit()->geographicalId().rawId(); + //} + } + + return smoothed; + } + + inline Trajectory KFFittingSmoother::fitOne(const TrajectorySeed& aSeed, + const RecHitContainer& hits, + fitType type) const { + throw cms::Exception("TrackFitters", + "KFFittingSmoother::fit(TrajectorySeed, ) not implemented"); + + return Trajectory(); + } + +} // namespace #include "FWCore/Framework/interface/ESHandle.h"